Trouble getting models to converge at lower temperatures

[Carly LaGrotta]

Hi Everyone,

I am attempting to run some shock tube simulations at lower temperatures (around 300 degrees) and am having trouble getting various models to converge, although it seems that for the models I am using the thermo data should be viable at those lower temperature.  I am running Cantera version 2.4 on linux. Some of the models I have tried using are as follows:

gas = ct.Solution('lamoreaux.cti')

gas = ct.Solution('konnovo.cti')

gas  = ct.Solution('gri30.cti')

gas = ct.Solution('glarborg.cti')

gas = ct.Solution('FFCM1_custom_updated.cti')

gas = ct.Solution('metan15.cti')

gas = ct.Solution('LeedsNOx20.cti')

gas = ct.Solution('sandiego.cti')

 The code I am using to run the shock tube along with the error can be seen below. I am wondering if this error is related to the thermo data in nature? If so I was wondering if there are any suggestions about how to run these low temperature simulations or if anyone knows an models that might be viable at these temperatures and contain larger carbon molecules as well as nitrogen compounds. Thanks, any help is really appreciated!! 

##############################################################################################################

import cantera as ct

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt 

def ShockTube(pressure,temperature,conditions,initialTime,finalTime,thermalBoundary,constant,rrc_list,frc_list):

    #gas = ct.Solution('lamoreaux.cti')

    #gas = ct.Solution('konnovo.cti')

    #gas  = ct.Solution('gri30.cti')

    #gas = ct.Solution('glarborg.cti')

    #gas = ct.Solution('FFCM1_custom_updated.cti')

    #gas = ct.Solution('metan15.cti')

    #gas = ct.Solution('LeedsNOx20.cti')

    gas = ct.Solution('sandiego.cti')

    

    

    

    

    mean_molecular_weight = []

    density_list = []

    gas.TPX = temperature, pressure*101325, conditions

    if thermalBoundary == 'adiabatic' and constant == 'V': 

        shockTube = ct.IdealGasReactor(gas,name = 'R1',energy= 'on')

        #print('inside here')

    if thermalBoundary == 'adiabatic' and constant == 'P':

        shockTube = ct.IdealGasConstPressureReactor(gas,name = 'R1', energy= 'on')

    

    

    sim=ct.ReactorNet([shockTube])

    if constant =='P':

        columnNames = [shockTube.component_name(item) for item in range(shockTube.n_vars)] 

        columnNames = ['time']+['pressure'] + columnNames

        timeHistory = pd.DataFrame(columns = columnNames)

    else:

        columnNames = [shockTube.component_name(item) for item in range(shockTube.n_vars)] 

        columnNames = ['time']+['pressure'] + columnNames

        timeHistory = pd.DataFrame(columns = columnNames)        

    

   

    t=initialTime

    counter = 0

    while t < finalTime:

        t = sim.step() 

        if constant =='P':

            state = np.hstack([t, shockTube.thermo.P, shockTube.mass, 

                    shockTube.T, shockTube.thermo.X])

            mean_molecular_weight.append(gas.mean_molecular_weight)

            density_list.append(shockTube.density)

        else:

            state = np.hstack([t, shockTube.thermo.P,shockTube.mass,shockTube.volume, 

                    shockTube.T, shockTube.thermo.X])

            mean_molecular_weight.append(gas.mean_molecular_weight)

            density_list.append(shockTube.density)

        timeHistory.loc[counter] = state

        counter+=1

    return timeHistory,density_list,mean_molecular_weight

#the function I am actually running 

timeHistory,density_list,mean_molecular_weight = ShockTube(0.986923,1000,{'CH3': 1.2094157562676408e-06,

 'HO2': 7.614839946870331e-07,

 'OH': 3.041863871824672e-06,

 'H2O2': 0.0001531112203220986,

 'N2O': 0.0007737003154574132},0,.003,'adiabatic','V',frc_violation,rrc_violation_list)

##############################################################################################################

The error is as follows...

CanteraError thrown by CVodesIntegrator::step:
CVodes error encountered. Error code: -3

At t = 0.00165949 and h = 1.32238e-10, the error test failed repeatedly or with |h| = hmin.
Components with largest weighted error estimates:
15: -85.9145
16: 80.3827
6: 5.82164
2: -0.0104749
8: 0.00405638
13: 0.000791296
31: -0.000781994
21: -0.000263437
7: 0.000121084
24: -7.17314e-05
************************

[Bryan W. Weber]

Hi Carly,

Can you please post some examples of these chemistry files?

Thanks,

Bryan

[Carly LaGrotta]

Hi Bryan,

 I attached some of the cti files here. Let me know if you need anything else. 

Thanks so much for the help in advance! 

Carly

[Xiaolei ZHU]

Hi Bryan,

I just want to say that I often see similar errors, and would love to know how to diagnose this.   I understand that some rates may be much higher that cause stiffness in the equations, leading to this integration failure, but I am not sure how to use the output here to find out which reaction or species is culprit.  

It would be nice if you can explain the procedure of how to figure out which entries are causing the problem. For example, what do the "components" indexes mean in this error message?  These numbers don't seem to be species index or reaction index.   

best,

Xiaolei Zhu

[Bryan W. Weber]

Hi Carly,

A couple of things. First, I couldn't reproduce this error using the most recent version of Cantera (2.4.0) or the previous version (2.3.0); could you please specify the version you're using, and if not 2.4.0, upgrade? Second, in your code there are two undefined variables frc_violation and rrc_violation_list, I just set them to None because they don't seem to do anything but it would help to include the complete code of what you're doing. You can attach it as a script, like you did for the CTI files.

Best,

Bryan

[Bryan W. Weber]

Hi Xiaolei,

There have been a bunch of posts on the group here about understanding these errors, so I'd encourage you to search and read those posts (e.g., here) Searching for "CVodes error encountered. Error code" (with the quotes) should do the trick for more topics. The component indices are (to my understanding) the indices of the solution variables, which change slightly depending on which Reactor subclass you're using. The order is given here: https://cantera.org/documentation/docs-2.4/sphinx/html/cython/zerodim.html#cantera.ReactorNet.sensitivities So for the most part, the component indices do refer to species indices, although they are shifted by 2-3 numbers.

These errors are often caused by bad thermodynamics or kinetics information, which you can debug using our handy example Notebook: https://github.com/Cantera/ncm-2019-materials/blob/master/thermo_debugging.ipynb

Best,

Bryan

[Carly LaGrotta]

Hi Bryan,

I have attached the python file and cti file. I am trying to run at 300 degrees sorry I mistakenly posted the wrong input conditions. Thank you for your time and patience. Please let me know if you need any more information.

Carly 

[Ray Speth]

Carly,

The essential advice from the Jupyter notebook that Bryan linked to is to look at the rate constants, particularly those for the reverse reactions. In this case, there are several potentially problematic reactions. Taking a rate constant of 10^20 as a somewhat arbitrary threshold:

for i in range(gas.n_reactions):
    if kr[i] > 1e20:
        print('{:4d}  {:.3e}  {}'.format(i, kr[i], gas.reaction_equation(i)))

provides the output:

  73  4.269e+25  CH3 + M <=> CH + H2 + M
  74  5.332e+26  CH3 + M <=> CH2 + H + M
 942  1.030e+30  CH + NH3 <=> 2 H + H2CN
1337  5.780e+23  H2NCHO + M <=> HCO + NH2 + M
1338  6.219e+27  H2NCHO + M <=> H + H2NCO + M

In addition, if you look up the solution components that are called out in the error message you posted (e.g. by calling shockTube.component_name(15) etc.) you will see that they correspond to CH3, CH2, and H, which suggests that the second reaction listed is the most problematic one. According to the paper referenced for this reaction (and the preceding reaction) in the mechanism file, these rate constants are fits for temperatures in the range 2253−3527 K, so it’s not surprising for them to cause problems at 300 K. If you search for this reaction in the RMG database, you will see that it is written much more frequently in the reverse direction, generally as a falloff reaction with a much more modest rate coefficient at 300 K.

Regards,
Ray