problem of PDSS module

[junfeng bai]

Hi, 

I'm trying to change the ideal gas law based on PDSS_IdealGas.cpp considering real gas effects. So I wrote a new module inheriting from PDSS and compiled successfully. 

However, no difference in results of my speciation simulation in ConstPressureReactor when I changed the equation-of-state model in yaml file. 

And then I changed the  equation-of-state model to constant-volume in yaml file, 

  equation-of-state:

    model: constant-volume

    molar-volume: 1.3 cm^3/mol

 still, no difference with the results using default  equation-of-state.

So I assumed that the PDSS module is not called at all.

And I found this introduction in 

http://web.mit.edu/2.62/cantera/doc/html/group__pdssthermo.html

Additionally, sometimes the VPSSMgr object will not call the PDSS object at all to calculate thermodynamic properties, instead relying on its own determination/knowledge for how to calculate thermo quantities quickly given what it knows about the PDSS objects under its control.

So when the PDSS module will be called and specified equation-of-state will be used?

The problem may be a little complicated, could anyone help me? Thank you.

Best,

Junfeng  

[Ray Speth]

Hi Junfeng,

First, I would suggest looking at the current documentation for Cantera available at https://cantera.org/documentation/index.html, rather than documentation that’s 10 years out of date. For instance, the VPSSMgr class mentioned there doesn’t even exist anymore, so the discussion of how it works is pretty much moot.

Second, I’m not sure that the PDSS classes are a very good way to try implementing real gas behavior - what ThermoPhase model are you putting them in? The one existing real gas equation of state implemented in Cantera (Redlich-Kwong) does not use this VPSS system, nor does the Peng-Robinsion implementation that’s being developed.

Last, if you’re trying to verify that your implementation is doing what you think it should, I would start with much simpler calculations than what effect it has on a reactor simulation, where the impact of the non-ideal terms could be quite small. I would focus on making sure that it’s giving you values you expect for simple properties like the pressure - density - temperature surface itself.

Regards,
Ray

​

[junfeng bai]

Hi Ray, 

Thank you very much for your advice. And I could try implementing real gas behavior based on Redlich-Kwong method EOS in cantera. 

However, I have another question when I tried to compare with Redlich-Kwong method. According to my test using 0-D constant-pressure reactor, the correction code of entropy is not actually called while the corrections of enthalpy and Cp are called. 

I have attached the code of RedlichKwongMFTP::entropy_mole() here. 

doublereal RedlichKwongMFTP::entropy_mole() const

{

    _updateReferenceStateThermo();

    doublereal sr_ideal = GasConstant * (mean_X(m_s0_R)

                                          - sum_xlogx() - std::log(pressure()/refPressure()));

    doublereal sr_nonideal = sresid();

    return sr_ideal + sr_nonideal;

}

Actually, I found 

    doublereal entropy_mass() const {

        return entropy_mole()/meanMolecularWeight();

    }

 

in ThermoPhase.cpp is also not called. 

Could this part be right without calling the correction of entropy? I don't think so. 

Could you help me? Thanks a lot. 

Sincerely, 

Junfeng

[Ray Speth]

HI Junfeng,

The mixture entropy isn't needed to evaluate the reactor governing equations, so it isn't really surprising that this method doesn't get called. For that matter, the specific heat doesn't appear in the governing equations either, but it is used by the iterative solver when setting the state at constant enthalpy and pressure, which does occur in the constant pressure reactor.

Regards,

Ray

[junfeng bai]

Hi Ray, 

Thank you for your explainations. However now I also find 

void RedlichKwongMFTP::getPartialMolarEnthalpies(doublereal* hbar) const

void RedlichKwongMFTP::getPartialMolarEntropies(doublereal* sbar) const

are not called either.  

Is this part right?

Thank you very much.

Sincerely,

Junfeng

[Ray Speth]

Junfeng,

Again, those quantities just don't appear in the governing equations. If you look at the implementation of the Redlich-Kwong model, you'll see that the total mixture enthalpy (which is used in the governing equations) isn't calculated by summing the partial molar enthalpies, but rather through other equivalent equations (although following the math to see that they are equivalent isn't necessarily straightforward).

Regards,

Ray

[junfeng bai]

Hi Ray, 

I understand that  partial molar enthalpies and entropies are not appearing in the governing equations. But shouldn't enthalpy and entropy for each species be used in determining reaction rate and equilibrium? As long as entropy for each species is 

used, RK real gas behavior should be considered, and correction of entropy should be called. Am I right?

Actually, I found the results of speciation simulation using RK method in Cantera and Chemkin are showing  reverse real gas effects. Inhibiting in chemkin while promoting in Cantera. (Using same Tc and Pc data). So I am checking the codes. Thank you again for your patient reply. 

Sincerely,

Junfeng

[Ray Speth]

Hi Junfeng,

The calculation of the equilibrium constant, to determine the reverse reaction rate, is done using the species chemical potentials. These are calculated directly, rather than using the identity . This eliminates some redundant calculations that would be performed if you computed the chemical potentials using the partial molar enthalpies and entropies.

Regards,
Ray

​

[junfeng bai]

Hi Ray,

I assumed that  the species chemical potentials with real gas corrections should be calculated by calling 

void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const

However, I found this part is not called either. Again, is this right?

Thank you very much. 

Junfeng

[Steven DeCaluwe]

Hi Junfeng,

Those calculations are done on the back-end, so to speak. The C++ routine that calculates the reaction rate will check to see if the state has changed, update the thermo and reaction rate coefficients as necessary, and then calculate whatever production rates have been requested. 

This is precisely one of the value propositions of Cantera—these steps are all automated for you, so that your high-level Python code does not have to. It makes the coding easier, more generalized, and less prone to bugs and errors. 

Best,

Steven

Sent from my iPhone

--
You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/125d2716-8330-46c3-ac73-aa46f1d925a8n%40googlegroups.com.

[junfeng bai]

Hi Steven, 

However, when I use RK EOS method to calculate, I think the program need to update thermo with corrections of real gas to calculate reaction rate. According to explaination from Ray, this is done using the species chemical potentials. So I believe the corrections code of  species chemical potentials in RedlichKwongMFTP.cpp should be called during calculations. 

However, according to my test, 

void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const

is not called at all. So I believe the real gas behavior is not correctly introduced into reaction rate calculations. 

Please correct me if I am wrong. Thank you very much!

Sincerely,

Junfeng

[Steven DeCaluwe]

Junfeng,

Please see the ‘updateROP()’ function, which begins on line 176 of GasKinetics.cpp:

https://github.com/Cantera/cantera/blob/4945645867d52d3a3792dcd7cec6ce152c518c79/src/kinetics/GasKinetics.cpp#L176

It first calls two routines: ‘update_rate_C’ and ‘update_rates_T’, which update the rate constants due to any change in relevant state variables.  `update_rates_T` also updates the equilibrium constant, Kc, via a call to `updateKc`.

To follow up on Ray’s post, the species chemical potentials are not needed at all, actually.  The equilibrium constant is a function of the standard state chemical potentials, which are called on line 120 of GasKinetics.cpp.  The real gas effects show up via the activity concentrations, which are equal to the molar concentrations, C_k (kmol/m3), multiplied by the activity coefficients, gamma_k.  For an ideal gas, gamma_k = 1 for all species, and for a real gas they are calculated according to the EoS.

I hate being the person who says “read my paper” to answer a question, but: this is all spelled out pretty carefully in the appendix (section A.3) to a paper we published in Combustion & Flame on real-gas effects in shock-tube simulations:

Kogekar, Karakaya, Liskovich, Oehschlaeger, Kee, DeCaluwe, “Impact of non-ideal behavior on ignition delay and chemical kinetics in high-pressure shock tube reactors,” Comb. Flame, 189, 2018.

https://doi.org/10.1016/j.combustflame.2017.10.014

You can find a preprint here: 

https://people.mines.edu/decaluwe/wp-content/uploads/sites/49/2021/05/KogekarEtAl_CombFlame_2018_NTC_Preprint.pdf

To directly answer your question, though: all real-gas effects in Cantera’s kinetics implementation show up in the activity concentrations, the calculation of which are equivalent to calculation of the chemical potentials, but do not require direct calculation of same.

Best,

Steven

To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/c7cc0095-2c74-4ada-94bc-71777f2b4036n%40googlegroups.com.

[Steven DeCaluwe]

Also, note that there are a few slight differences in some of the equations in that paper and in Cantera’s Redlich-Kwong implementation.  These have to do with how the standard state is defined (i.e. how certain energy terms are split between the “standard state” and the “activity”) but do not result in any differences in the resulting calculations.

SD

To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/66210D43-7980-4D1E-9BE1-B787F2EF1CCB%40gmail.com.

[junfeng bai]

Hi Steven,

Thank you very much for your explainations. I am learning your paper now and I want to know more details for rate of progress equation.

However, it seems not appeared in [50] G. Froment , K. Bischoff,Chemical reactor analysis and design, Wiley, New York, 1979 .

And PDF of 

[51] C. Eckert , M. Boudart , On the use of fugacities in gas kinetics, Chem. Eng. Sci. 18 (1963) 144–147 .

is not available on scholar.google. 

Could you send me this paper? Thank you again.

Sincerely,

Junfeng

[junfeng bai]

Hi Ray, 

I have successfully changed EOS equations, H and Cp corrections and activity concentration calculations code using new method based on RK method. Timehistory results could be outputted normally at 100 atm and results are reasonable. However, at 200 atm and 300 atm, it seems like  t = reactorNetwork.step()  is very hard to finish because the time interval between steps become shorter and shorter. 

Do you have any idea what could cause this problem? Could you tell me how the time interval is decided in Cantera and in which part of code? Or could time interval be changed manually?

Thank you very much!

Sincerely, 

Junfeng

[Ray Speth]

Hi Junfeng,

Time step sizes in reactor network integration are determined by the Sundials integrator. The time step will be changed during the course of integration to make sure that the solution satisfies the error tolerances specified on the reactor network. If you find that the integration is becoming so short that you can't complete integration, two possible causes are (a) that you have some reactions where the rates in one direction or the other are too high (and non-physical), or some sort of discontinuity/asyptote in your governing equations. The best general advice I can offer is to save all of your reactor state variables to after each timestep so you can examine them and see if there are any specific variables that are not behaving as you would expect.

Regards,

Ray

[junfeng bai]

Hi Ray，

Thank you very much for your advice. I tried to save all variables, including temperatures, pressures, molarvolume, non-ideal H and Cp corrections values, activity coefficients. However, they are all  behaving as I would expect. According to my test,  the time interval is becoming too short only when I used non-ideal Cp corrections, I think this may affect the solution of   governing equations? The  Cp corrections values I used are similar to RK method, but RK method could always finish normally. Real-gas-effect is making Cp bigger than ideal gas, so I do not understand how the change of  Cp  could cause discontinuity/asyptote in my governing equations. 

Could you give me some more advice according to my new tests? Or, since I am assuming my code is right actually, is there any simple method could lower the standard, sacrifice some accuracy to make time interval larger? Thank you very much!

Regards,

Junfeng
 

[William R. Smith]

Hi Junfeng/Ray,

If all the input is correct, what Junfeng is probably seeing is evidence of a well-known problem that occurs when integrating certain differential equations (DEs) with certain algorithms.  This occurs when the DE's are "stiff"; as evidenced by  some require very small time steps in comparison with others; since only a single time-step is possible in an integration algorithm,  this means that a very small time step must be used for all the equations, and it becomes very slow.  Thie "stiffness problem" was noticed some time ago by a chemist (Hirschfelder) and ultimately remedied remedied by using a different numerical integration scheme.  C.W. Gear in 1971 was one of the early developers of methods that alleviated the problem.

Ray, I don't know the details of the "Sundial" algorithm, but perhaps it doesn't perform well for stiff DEs?

William

[Ray Speth]

Junfeng,

Without seeing the details of your code or the full results of your simulations up to the point where they fail, it's pretty much It's impossible to give much more concrete advice. If you have made your changes to a branch of the Cantera Git repository, you could push that back to your fork on Github and someone might have the time to look at your work.

William,

The Sundials library contains several ODE solvers. The one used by Cantera is CVODES, which implements a variable-order, variable timestep BDF method that is closely related to Gear's method and is well suited for the stiff problems that arise in chemical kinetics. However, even this algorithm has its limits. In my experience, the most common failure modes for the CVODES solver in Cantera are when the problems posed have non-physical sources of stiffness, such as reaction mechanisms with rates that exceed collision limits by 10 orders of magnitude or more.

Regards,

Ray

[junfeng bai]

Hi Ray, 

Thank you very much. Could I enlarge error tolerances to make simulation faster? Is there any option about accuracy in Cantera or could I change tolerances in source code?

Sincerely,

Junfeng

[Bryan Weber]

Hi Junfeng,

Sure, this is documented here: https://cantera.org/documentation/docs-2.5/doxygen/html/d7/def/classCantera_1_1ReactorNet.html#a4d7bd81259e8942445eab714e947490b The rtol and atol methods of the ReactorNet control the tolerances.

Best,

Bryan