mutationpp / Mutationpp

The MUlticomponent Thermodynamic And Transport library for IONized gases in C++
GNU Lesser General Public License v3.0
108 stars 58 forks source link

Get T, Tv temperatures #108

Closed CatarinaGarbacz closed 4 years ago

CatarinaGarbacz commented 4 years ago

Hello,

I am using Mutation++ to initialize a thermal bath domain in a CFD solver, with the following conditions:

To debug the problem, I was reproducing this problem just in Mutation++ and I have the following results:

case 1: when I input the energies and do mix.getTemperatures, I get T = 57K , Tv = 82K

case 2: when I input the temperatures T = Tv = 57K, and I compute the energies, getting the same energy_density_mixture = -6166.675192 and vibrational_energy_density_mixture =1.461811921e-14, and then I use these same values to do mix.getTemperatures, I get T = 57K , Tv = 57K

And this happens for all values temperatures in equilibrium up to T = Tv = 95K. Above this value, things go smoothly in both cases and in the CFD solver I get a flow initialized and staying in thermal equilibrium. Things also go well when a mixture of CO2 is used, for any values of temperature.

Does anyone have an idea of what is the issue here? Is it that up to 95K we physically cannot have thermal equilibrium in a mixture of standard air, and in case 2 I am kinda forcing it? Or am I missing something else?

Any help is appreciated

jbscoggi commented 4 years ago

Hi @CatarinaGarbacz, can you post a minimum working example showing this behavior? I will use that to try and debug what's going on. Thanks

CatarinaGarbacz commented 4 years ago

Sure JB thanks. Just comment in or out the lines corresponding to case 1 or 2 to obtain the two different behaviors.

include "mutation++.h"

using namespace Mutation; using namespace Mutation::Thermodynamics; using namespace std;

// Main entry point int main() {

MixtureOptions opts("air5");
opts.setMechanism("air5");                 
opts.setStateModel("ChemNonEqTTv");       
Mixture mix(opts);                       

int ns = mix.nSpecies(); int i;

double cs0[ns], rhos0[ns], T[2], E[2], rhoE[2], Tnew[2];

double rho = 0.02373253141;

cs0 [0] = 0.0;
cs0 [1] = 0.0;
cs0 [2] = 0.0;
cs0 [3] = 0.77;
cs0 [4] = 0.23;

for (i=0; i<ns; i++) rhos0[i] = cs0[i]*rho;

// ---------- Case 2 ---------- //

T[0] = 57.0; T[1] = 57.0;

std::cout << "T = " << T[0] << std::endl << std::endl; std::cout << "Tve = " << T[1] << std::endl << std::endl;

mix.setState(rhos0, T, 1);

mix.mixtureEnergies(E);

rhoE[0] = rhoE[0]; rhoE[1] = rhoE[1];

// ---------- Case 1 ---------- //

// rhoE[0] = -6166.675192; // rhoE[1] = 1.461811921e-14;

// ---------- Common to case 1 and 2 ---------- //

std::cout << "rhoE = " << rhoE[0] << std::endl << std::endl; std::cout << "rhoEve = " << rhoE[1] << std::endl << std::endl;

mix.setState(rhos0, rhoE, 0);

mix.getTemperatures(Tnew);

std::cout << "Tnew = " << Tnew[0] << std::endl << std::endl; std::cout << "Tvnew = " << Tnew[1] << std::endl << std::endl;

return 0; }

jbscoggi commented 4 years ago

Hey @CatarinaGarbacz, I looked into this and it seems that it's not really a bug, it's just the tolerance that is set for the Newton iteration to solve for the temperature. What's happening is when you explicitly set the temperatures first (case 2) those temperatures are used as the initial guess for the Newton iteration when you call the second setState with energies. Since those temperatures are exact, it just returns with no iterations. On the other hand, if you set the energies first (case 1) then it will use whatever the the StateModel was initialized (300K I think) and perform a Newton iteration to get the temperatures. The iterations stop once the 2-norm of the difference in given energies and those computed from the temperatures satisfies a relative and absolute error tolerance criterion, which are both quite strict already at 1e-12. It just so happens that in the region below ~90 K, this tolerance is satisfied with pretty much any Tv because the vibronic energy is essentially zero below this temperature. The only energy playing a role in the iterations is the total which is used to set T=57. That's why you get this behavior. Other than decreasing the relative tolerance further, I'm not really sure what to do. If you decrease the tolerance, you will get a "better" Tv but you will pay for it by having to perform more iterations every time you set the state with energies. If you want to play with that, you can change the rtol parameter in the solveEnergies method in ChemNonEqTTvStateModel.cpp.

The bottom line is that from a thermodynamic standpoint at least, the choice of temperature below ~90K will have essentially no impact on the mixture thermodynamics. I would suggest that you just use a minimum temperature of 100K and you will avoid this problem most likely. Is it crucial to use 57 for some reason?

CatarinaGarbacz commented 4 years ago

Thanks JB for looking into it and for such a detailed explanation. In my debugging process I had already got to some of those conclusions, but you provided me with a more complete understanding of what's happening.

To answer your question, indeed it is not mandatory to have 57K and I had already replaced this for 100K test-cases. But I thought it would still be good to understand the problem.

I guess I have one remaining question about your comment "The bottom line is that from a thermodynamic standpoint at least, the choice of temperature below ~90K will have essentially no impact on the mixture thermodynamics." Does this mean that having T=57K with Tv=57K or Tv=82K would result in the same total energy and same vibrational energy, hence with no impact in the freestream and, consequently, disturbed flow? cause in the end the flow is solved through the computation of energies. and from what I understand, you're telling that, up to Tv ~ 90K, the vibrational energy is always the same equal to zero, despite of the value of Tv.

Cause if this is the case, I would keep my T=57K , Tv=57K/82K work going and have the T=Tv=100K as an additional test-case

jbscoggi commented 4 years ago

Hey @CatarinaGarbacz, to answer your last question, yes. :) I don't think you will see any difference using T=57 or T=100K. For all practical purposes, the vibronic energy is zero in both cases. Setting T=100 will just "look" nicer because the temperatures you get back in the free-stream will probably still be 100K.

CatarinaGarbacz commented 4 years ago

ok thanks a lot for the feedback ! :)

jbscoggi commented 4 years ago

Of course! I will close this issue then.