Equilibrium state-model hangs when computing temperature

[jbvos]

We coupled the Mutation++ equilibrium solver to a CFD solver. We try to make a calculation for 5 species air, start from free stream conditions, but the solver hangs in the 2nd step when sending the state to Mutation++. In addition there is the issue of error messages in EquilStateModel.cpp (around line 112) , it fills the disk. I also added a break to quit the iteration loop if it is not converging.

The same happens when we start from a perfect gas solution.

[jbscoggi]

Hi @jbvos, can you provide the exact input to the setState() function that is causing it to hang? That way, I can try to reproduce the problem.

[jbvos]

This is the print statement in the fortan program:

  if (n == 7485) then
  print'(a,i6,2e21.12)','rho=',n,rho(n),rhoe
  print'(a,i6,2e21.12)','p,t=',n,p(n),t(n)
  endif

           call mpp_set_state(rho(n),rhoe,impp)

This is for the first block

rho= 7485 0.510484261795E-03 0.977054025998E+04
p,t= 7485 0.448499999125E+02 0.379949999259E+03

This is for the 2nd block and here it hangs

rho= 7485 0.510671756357E-03 0.973033904733E+04
p,t= 7485 0.448499999125E+02 0.379949999259E+03

[jbscoggi]

Can you please explain how the pressure and temperature are computed/used here? They are unchanged between between the two calls.

[jbvos]

I just print them for info, they are from the previous time step. The way I compute them:

           call mpp_get_temperatures(tnew)
           pnew  = mpp_pressure()

I played around with relaxation:
!
delt = told - tnew
if (abs(delt) > 500.) then
delt = sign( min(abs(delt),0.5*told) , delt)
endif
!
t(n) = told - delt
t(n) = min(max(t(n),50.0),8000.0)

This helps a little, but not much

[jbvos]

I have to leave in about 5 minutes, but I will check again this evening

[jbscoggi]

OK, no problem. I will not have a chance to get into this tonight, but I just wanted to get a few details to help when I have time to take a look.

[jbscoggi]

Hi @jbvos, I have tried to reproduce your problem with no success. Can you please specify exactly which state model name you are using and the value of impp in the code snippet above?

[jbscoggi]

To be more specific, here is the test case I have tried

    MixtureOptions opts;
    opts.setSpeciesDescriptor("N O NO N2 O2");
    opts.setStateModel("Equil");
    Mixture mix(opts);

    double rho  = 0.510671756357E-03;
    double rhoe = 0.973033904733E+04;

    mix.setState(&rho, &rhoe);

    std::cout << mix.T() << '\n';
    std::cout << mix.P() << '\n';

This gives T = 5301.53 K and P = 1280.28 Pa.

[jbvos]

I am not surprised, I tried to do the same but then it works. I was wondering whether there is an accumulation of errors or a problem with an overwrite.

I use the Statemodel air5, impp equals zero.

[jbscoggi]

air5 should be the mixture name, I was wondering what you chose as the state model name. It should be either Equil or EquilTP. If its EquilTP, then it could explain the problem...

[jbvos]

Equil

[jbscoggi]

Hi @jbvos, I still have not been able to reproduce the problem. Do you think you could provide the state that you set just before the one before it hangs on? In other words, it looks like you are printing only for cell 7485, but can you also provide 7484? Until I can reproduce the problem, it is not really possible to debug. Thanks.

[jbvos]

I know, I tried to isolate the problem too, but did not manage. I am traveling this week, but I will try to see what I can do this evening.
About the GSI model, I tried to run the different cases we used for the Ablacat project, but I get the error message

M++ error: invalid input.
input: key = gamma_energy
Provider <gamma_energy> for type SurfaceProperties is not registered. Possible providers are:
gamma

(I used the gsi files we used from Ablacat)

Thanks for all you help

[jbscoggi]

@grgbellasvki, can you take a look at the last comment by @jbvos? Thanks.

[grgbellasvki]

Considering the gsi problems there are two things:

1. I think you need to update your Mutation++. GSI was till very recently only in a personal branch, and was moved to master one month ago, when we added the energy balance. For this reason it is complaining that it doesn't find "gamma_energy", but only "gamma".
2. Since the version you were using during ablacat some things have changed for clarity. For example the key gamma has been changed to phenomenological.

Feel free to send me your GSI file and I will update it to the new format.
Can you please remind me if you are using the C++ functions or the fortran wrapper?

[jbvos]

Ok, thanks, please let me know exactly which version of Mutation++ I need to have/download. I also attach the different gsi files we used in Ablacat. I use the fortran wrapper. Many thanks for your help
gsi.tar.gz

[jbvos]

I am for a few hours in my office, and I printed the input for a few other points in the grid:

rho= 7480 0.456073574855E-03 0.617581118526E+04
rho= 7481 0.456083566995E-03 0.617609636132E+04
rho= 7482 0.456145930175E-03 0.615630267017E+04
rho= 7483 0.456199742381E-03 0.613609481969E+04
rho= 7484 0.456141749894E-03 0.615371133006E+04
rho= 7485 0.510484261795E-03 0.977054025998E+04
rho= 7486 0.409605406221E-03 -0.484543062766E+02
rho= 7487 0.456138613203E-03 0.615585471262E+04
rho= 7488 0.456146489417E-03 0.615669547124E+04
rho= 7489 0.510548345834E-03 0.977537146738E+04
rho= 7490 0.510525621144E-03 0.977457690101E+04
rho= 7480 0.456060713159E-03 0.617563493130E+04
rho= 7481 0.456062127776E-03 0.617859464097E+04
rho= 7482 0.456094073892E-03 0.617449703052E+04
rho= 7483 0.456143676864E-03 0.615250467545E+04
rho= 7484 0.456241475056E-03 0.611579964293E+04
rho= 7485 0.510671756357E-03 0.973033904733E+04

The interesting thing is that the state when it hangs is not far from the state at other points for which it does not hang.

Is there an option to print debug info from mutation++ so that you get a better idea what is going on?

I have to travel again, back on Thursday so I am probably not very reactive

[grgbellasvki]

The latest version of Mutation++ should work fine. Note, that the fortran wrapper has a few functions slightly changed (mainly name).
Steady state pyrolysis is not available for the time being. Since nobody is currently using it, I was planning to implement it only on request. Do you need it?

Can you please send me the single gsi file you want to use? I cannot modify all of them.
When you update M++, you can find many examples in /tests/data/gsi.
The xml files seb_carbon_ablation_park, smb_carbon_ablation_park and smb_air_catalysis_gamma_ions are the ones you need to refer to.

[jbvos]

I attach one of the gsi files. For steady state pyrolysis, I think we never used it, but I need to check

gsi.zip

Thanks a lot for your help

[jbscoggi]

Hi @jbvos, I was wondering if there is any updates to this issue? I'm not sure if you have tried using the latest version of the library and still are having problems. If everything is working, I will close this issue.

[jbvos]

I just downloaded the latest version of Mutation++ and I will try it out. I will come back to you
Best regards, Jan

[jbvos]

Dear James

I downloaded and compiled the new Mutation++ version, but get problems with linking:

undefined reference to mpp_set_wall_state_' undefined reference to mpp_get_wall_state_'

Have these routines been changed, or are they absent in the Fortran version?
Thanks for your help,
Best regards Jan

[jbscoggi]

Hi @jbvos, I think these functions were renamed a while back when the GSI branch was more fully supported. If you replace wall with surface in both function calls in your code, then I believe it should work.

[jbvos]

Thanks, this worked. I will now go to the testing

[jbvos]

It took some time. I had to copy the gsi data directory of an older distribution to the new one, and then I ran probably into another compatibility issue:

M++ error: error parsing file.
file: /soft/mutationpp-master/data/gsi/gsi_oxi.xml
line: 3
Error in surface_properties for the gsi imput file. Surface properties should provided. If no properties are needed they should set to 'none'.

The gsi_oxi.xml states:

<surface_properties name = "none" >
</surface_properties>

<diffusion_model model="first_order">
    <mass distance= "1"> </mass>
</diffusion_model>

<production_terms>
    <surface_chemistry>
        <!-- 1 -->
        <reaction type= "ablation" formula="C-s + O => CO">
            <gamma_T pre_exp="0.63" T="1160.0" />
        </reaction>

    </surface_chemistry>

<surface_radiation emissivity="0.86" T_environment="0."/>

</production_terms>

[jbscoggi]

@grgbellasvki can you help here please?

[grgbellasvki]

The final format for the GSI file for carbon oxidation is shown below:

<surface_properties type = "ablation" >
    <surface label="b" species="C"
             enthalpy_surface="0." />
</surface_properties>

<solid_properties 
    virgin_to_surf_density_ratio = "1."
    enthalpy_virgin = "0.">
</solid_properties>

<surface_features 
    solid_conduction = "none"
    surface_in_thermal_equil = "true"
    gas_radiation = "false">
</surface_features>

<surface_chemistry>
    <!-- 1 -->
    <reaction type= "ablation" formula="C-b + O => CO">
        <gamma_T pre_exp="0.63" T="1160.0" />
    </reaction>
</surface_chemistry>

<surface_radiation emissivity="0.86" T_env="0."/>

Additional examples can be found in "mutation++/tests/data/gsi/". For the simulations, please pay attention to the option solid_conduction to whether "none" or "steady_state" is more physical for your testcase.

Let me know if you have any other issues.

[jbvos]

I took the file you provided:

<surface_properties type = "ablation" >
    <surface label="b" species="C"
             enthalpy_surface="0." />
</surface_properties>

<solid_properties
    virgin_to_surf_density_ratio = "1."
    enthalpy_virgin = "0.">
</solid_properties>

<surface_features
    solid_conduction = "none"
    surface_in_thermal_equil = "true"
    gas_radiation = "false">
</surface_features>

<surface_chemistry>
    <!-- 1 -->
    <reaction type= "ablation" formula="C-b + O => CO">
        <gamma_T pre_exp="0.63" T="1160.0" />
    </reaction>
</surface_chemistry>

<surface_radiation emissivity="0.86" T_env="0."/>

but when trying to run it I get the message:

M++ error: logic error.
file: /soft/mutationpp-master/src/gsi/Surface.h
line: 148
setConductiveHeatFluxModel can be called only when solving the surface energy balance!

So I do something wrong. Is there a conversion program which converts the GSI files we used during the Ablacat project to the new format?

Many thanks

[grgbellasvki]

Unfortunately no, but the changes are very simple. Most of the xml elements are completely unchanged, the only thing different is the overall structure of the file, which is more consistent with the current structure of the code.

It seems the comments here remove the parent xml element so I attach the whole files.
Is both the mass and energy balance solved on the surface? If yes, the following file should be used:
seb_carbon_oxidation_park.xml.zip

If no (only mass is solved):
smb_carbon_oxidation_park.xml.zip

In the latter case you cannot call function:
mpp_set_cond_heat_flux, because it is a logic error. This is what is causing Mutation++ to exit. If you could post a code snippet of how you invoke the gsi routines it would help.

[jbvos]

This was helpful many thanks. The code started to run the case using both the smb_carbon_oxidation_park and seb_carbon_oxidation_park files
It will probably take some time before I have a result, I will let you know probably early next week.
Again many thanks

[jbvos]

Air_P250_T280_Tw2747.zip
Air_P250_T280.zip

Sorry, it took some time but the calculations are rather long. I attach 2 results from GSI cases of the Ablacat project, one with imposed wall temperature, in the other one the temperature at the wall comes from the GSI calculation. From the CFD solver, the only modification compared to the old M++ version is the change from: call mpp_set_wall_state(rhoi,twl,1) to call mpp_set_surface_state(rhoi,twl,1)

The zipped files contain both the stagnation and wall values of different variables

The good thing is that both calculation ran. The calculation with imposed wall temperature (and using smb_carbon_oxidation_park.xml) is not too bad, although there are surprises on the wall around x=0.01, not clear what is happening there. The calculation in which the wall temperature comes from the GSI routines pose more problems.

I also attach the routine from which I compute the ablation

bcvis.zip

Any help is appreciated. Note that in the Ablacat project we also computed with Asterm, I do hope that this is still possible

Many thanks