How to validate Cantera OD Ideal Gas Reactor Manually

[Partth]

Hi All,

I have set up a OD reactor using Cantera's IdealGas reactor function and it works well. I wanted to write an explicit OD solver where I manually integrated the rate equations myself to validate/compare against the Cantera results for species mass fraction and temperature. Currently, I believe I put in all rate information and setup my rate equations, but my ODE solver seems to diverge.

I simply chose the base Zel'dovich mechanism for NOx production using wikipedia (https://en.wikipedia.org/wiki/Zeldovich_mechanism) for simplicity but unsure if it right. I also assumed a high enough initial temperature and pressure to allow species evolution over a short time. I typically do not work with such mixtures, but I just needed a simple mechanism to try and resolve my manual ODE solver; being said if my initial concentrations, temperature or pressure are unrealistic, please feel free to correct me/change it.

I have attached my mechanism .yaml file, cantera reactor and my manual 0D reactor codes here or in subsequent messages. 

[Partth]

For reference, I calculated the change in enthalpy for each reaction by taking the difference of 1) the sum of enthalpies of formation of products and 2) sum of enthalpies of formation of reactants at the specified T. These values were evaluated from the NIST Table: https://janaf.nist.gov/. If someone can see what I am doing incorrect in my manual approach (man_zel.py) I would appreciate the feedback as I am not getting the results shown in the Cantera approach (ct_zel.py). I am doing this exercise to see if I am understanding the steps correctly. I tried to look through the Cantera documentation and understand they use an adaptive step integration method, but I was not sure if there was any preconditioner/other steps taken to ensure solver stability. Otherwise, I could be over-complicating it in my approach so any feedback would be helpful and appreciated. Thank you for your time!

-Partth

[Partth]

I would like to first see if I can understand what Cantera is doing for a simple mechanism; I have some mechanisms which involve reactions with a generic species (+M); how would one similarly write out the explicitly rate equation/calculation using the collision efficiencies? I just want to better understand how Cantera would evaluate this if someone can shed light on this as well. Also, please excuse any typos in the prior messages if present as I have written this in haste.

[Ray Speth]

Hi,

Explicit ODE integration methods generally do not work for solving combustion chemistry. You need to use an implicit numerical method such as the BDF method implemented by the SUNDIALS Cvode solver. You’re also very unlikely to have any success without a method that implements an adaptive time step size. I would not recommend trying to code up such a solver yourself. Instead, you could use an existing ODE solver. There are several options available in SciPy. We even have an existing example of doing this: custom.py.

Regarding third body reactions, the calculation of the effective third body concentration is done as weighted sum of all the species concentrations multiplied by their corresponding efficiency. You can get this value for each reaction using the gas.third_body_concentrations property in Python. The way this concentration is used depends on whether you’re talking about a simple third body reaction or a falloff reaction. Both cases are explained in the Cantera Science Documentation.

Regards,
Ray

​

[Partth]

Hi Ray,

Thanks for this. Actually, I figured as much when looking into the the SUNDIALS solver but I figured it was a decent exercise. In any case, with the example above, I found issues using my explicit integration approach as you said; for a 1 step mechanism, I had rewritten the code and it somewhat worked, but even this is sensitive to inputs. In practice, I would use the Cantera reactor setup/option natively, but I was asked to verify this manually for my research group to see if it would be possible without any SciPy package, hence trying it manually. Thanks again!

-Partth