CSTR convergence error!

[PrayShah]

I am trying to run CSTR with Cantera v3.0.0 for two reaction mechanism, namely Mechanism-1 and Mechanism-2. Mechanism-1 is subset of Mechanism-2. The CSTR works fine with Mechanims-1 however Mechanism-2 gives convergence error. I have setup CSTR job as shown in one of Cantera's example on the website Both mechanisms are working prefect with Chemkin! I have attached the mechanism files and jupyter notebooks I am using to run Cantera. Could you please help me to fix the error I am facing?

CSTR_Reactors.zip

[speth]

Using the provided notebooks, I saw integrator errors with both mechanisms. There are three issues I noticed:

• Relaxing the default integrator tolerances, especially atol can be lead to erratic integrator behavior, since the solution does depend on the behavior of species with small concentrations.
• Probably the most common issue leading to integrator errors is reactions that have non-physical rate constants (by many orders of magnitude), which introduce artificial stiffness. You can find such rates (at a particular temperature) by running something like the following:

  gas.TP = 750, 101325
  kf = gas.forward_rate_constants
  kr = gas.reverse_rate_constants
  print('Forward:')
  for i in range(gas.n_reactions):
  if kf[i] > 1e13:
      print(f"{i:>5d}: {gas.reaction(i).equation:30s}  {kf[i]:.2e}")
  print('Reverse:')
  for i in range(gas.n_reactions):
  if kr[i] > 1e13:
      print(f"{i:>5d}: {gas.reaction(i).equation:30s}  {kr[i]:.2e}")

• In the second notebook, the integration loop is using the wrong time when calling ReactorNet.advance:

  for i in range(N):
      SimTime += dt
      temp_dependence.append(stirred_reactor.thermo.state)
      reactor_network.advance(max_simulation_time)

[PrayShah]

Thank you for the response. I corrected the time used in integration loop for second mechanism and got the same error as first mechanism. So next I checked reaction rates that have high rate constants by the approach you suggested and here is what I got:

Forward:
  108: CH2O2H <=> CH2O + OH            2.56e+14
 1602: CCYCCOOC-I2 <=> CHOIC3H6O       2.58e+13
Reverse:
  862: C3H3 + H <=> C3H2 + H2          1.05e+13
  873: C3H2(S) + M <=> C3H2 + M        3.76e+16
 1851: C4H63,1-2OH <=> C4H612 + OH     4.22e+15
 2245: C4H5-I + O2 <=> C2H3 + CH2CO + O  1.55e+14
 2685: IC3H5CHO + IC4H7 <=> H15DE25DM-SO  8.47e+13
 2699: C2H3COCH3 + C4H71-3 <=> C8H131-5,3-4,TAO  1.54e+15
 2722: C6H10D24 <=> 2 C3H5-S           2.07e+13
 2726: C2H3COCH3 + C4H71-3 <=> C8H131-5,3,TAO  5.87e+18
 2727: C4H71-3 + SC3H5CHO <=> C8H131-5,3,SAO  1.28e+15
 2747: C4H71-3 + SC3H5CHO <=> C8H132-6,SAO  1.55e+15
 2995: CH3OCH2OCH3 <=> CH3OCH + CH3OH  1.06e+18
 3117: CH3OCH2OCH2OCH3 <=> CH3OCH + CH3OCH2OH  1.19e+19
 3172: CH3OCH2OCHO + HO2 <=> CH3OCH2OCO + H2O2  3.83e+13
 3175: CH3OCH2OCHO + O2 <=> CH3OCH2OCO + HO2  4.62e+15
 3180: CH3OCHOCHO <=> CH3OCHO + HCO    1.48e+13

Reaction 2995, 3117 and 3180 are of importance and other reaction will anyway not happen at that temperature. Now if I pass the same mechanism in Chemkin with their PSR model it works fine so I am assuming that Chemkin is somehow managing the stiffness. Similarly is there a way I can modify some parameter of Cantera so that the integrator can manage such stiff reaction model?

[speth]

I noticed another issue with your reactor network construction. The exhaust reservoir is being created before setting the state of the gas object to your initial condition, which means that the initial condition for stirred_reactor is higher by 0.03 * 101325 = 3039.75 Pa. Further, the PressureController is being set up without a value for its pressure_coeff, which means that it defaults to 1.0 kg/s/Pa, creating an initial mass flow rate out of the system of 3039 kg/s. Given the initial mass in your system is ~ 1e-4 kg, this creates a very short time constant for the system that gives the integrator a lot of trouble. Furthermore, the valve pressure coefficient will drive the pressure in the stirred_reactor to the pressure of the exhaust, rather than maintaining the initial pressure you specified.

I would recommend the following changes to your code:

1. Create the ehxaust object after setting the initial state of the gas in your loop over the reactor temperature
2. Set a much smaller valve coefficient, for example by adding K=1e-3 to the keyword arguments when creating the pressure_regulator object.
3. Use a tighter value for atol, and leave the default rtol, for example rtol = 1.e-9; atol = 1.e-20.

These changes were sufficient for me to get reasonably efficient solutions using Mechanism 2.

Regarding the reaction rates, if those reactions are important, then it's presumably important to get their rates right. For reactions 2995 and 3117, those rates are physically impossible. The problem of non-physical rate constants in combustion mechanisms has been discussed in the literature. I would refer you to two papers in particular:

• https://doi.org/10.1016/j.combustflame.2017.08.005 suggests an upper bound for bimolecular reactions of $k_{coll} = 2.9 \times 10^{15} \rm cm^3/mol/s$ (or $2.9 \times 10^{12} \rm m^3/kmol/s$ in Cantera's units) with stricter limits that can be calculated if some transport properties are available.
• For unimolecular decomposition, https://doi.org/10.1016/j.combustflame.2020.06.010 suggests limits ranging from $6\times 10^{12}\ 1/\rm s$ at 300 K to $6\times 10^{13}\ 1/\rm s$ at 3000 K. In this case reaction 2726's reverse rate is the biggest outlier.

[PrayShah]

I implemented the changes to my reactor code and it is working fine. Thanks you for giving me insights. Further I look into reaction 3117 and 2995 rates. I compared kinetic model generated and Master Equation rates and there is a huge error with my model rates. Thank you for pointing me in the direction of non physical rates. I will be able to fix this now and will keep in mind to do these checks from now onwards!