Is Cantera ODE solver really trustworthy?

[Marc Fischer]

Hi folks,
I'm facing a problem that makes me call the reliability of Cantera ODE solver into question.

I'm writing functions in order to know whether an explosion (exponential increase of the temperature) takes place during the interval [0;t_end].
Currently, I'm checking whether Cantera is able to predict both the species concentration and the temperature profiles.

The function "assess_explosion" simulate the temporal evolution of the system with Cantera ODE solver.
I defined

delta_t = t_end/n_t; # Time step
r = ct.Reactor(gas, volume = V);
sim = ct.ReactorNet([r])

and then in a while loop:

time += delta_t # Next time point.
sim.advance(time)


I also created the function "homemade_assess_explosion" that solves the ODE system numerically through a one-step explicit method.
The function "analytical_assess_explosion" computes the species mole fractions using analytical formula under the assumption that the temperature remains constant.

The reaction I'm considering all along is:
reaction('C3H8 + 5 O2 => 3 CO2 + 4 H2O', [8.6E+05, 0, 125520],order="C3H8:0.1 O2:1.65")


I first tested these functions by considering a highly diluted mixture, so that the temperature remains pretty much constant:

X = 'C3H8:0.00004058,O2:0.0002029,N2:0.9999975652' #'C3H8:0.000407,O2:0.002029,N2:0,AR:0.99'
T0 =  1178.92;

As you can see in Figure "X_C3H8_isotherm.png", the three methods predict the same mole fraction profile.
At a constant temperature, I was able to get such a good agreement for other values of the reaction orders and initial temperatures as well.
I am thus confident that my home-made numerical solving of the concentration ODE is correct.

I then considered a situation where the temperature would no longer stay constant:

T0 =  878.92; # Initial temperature
X = 'C3H8:0.04058,O2:0.2029,N2:0.9975652' #'C3H8:0.000407,O2:0.002029,N2:0,AR:0.99'

The results of Cantera and of my home-made numerical solver are very similar (see "T_878.92.png" and "X_C3H8_878.92.png").

Nevertheless, if I consider a lower temperature (and the same composition),
T0 =  778.92; # Initial temperature
the two models considerably differ: while my home-made ODE solver predicts an exponential increase of the temperature in the time interval, Cantera ODE solver only predicts a slow linear increase thereof (see "T_778.92_n_t1E+04.png" and "X_C3H8_778.92_n_t1E+04.png").
If, however, I strongly increase the number of time points (and thus strongly decrease the time step)

n_t = 1E+05;
T0 =  778.92;

my solver and Cantera solver agree again (see "T_778.92_n_t1E+05.png" and "X_C3H8_778.92_n_t1E+05.png").

This shows that my home-made explicit solver was right.

I face the same problem all the time: to get accurate results with Cantera, I must use very small time steps, even though I only consider one reaction!

Do you know where this problem stems from?
Is it due to the "unusual" reaction orders of the reactants C3H8 and O2 that results in a very stiff ODE system?

How can I trust Cantera for computations involving large reaction mechanisms if it isn't even able to simulate ONE reaction without having to resort to very small time steps?


P.S: I could send you the source code I wrote but it is complicated and I'm not sure it'd help you better understand the situation...

[Bryan W. Weber]

Hi Marc,

As you've seen, Cantera uses SUNDIALS (specifically, CVODES) to solve the ODEs for the 0-D (specifically, instances of ReactorNet) chemical systems. I would say that CVODES is quite reliable in general, although like all software I'm sure it has bugs, and there are probably problems for which it is not suitable. I doubt this is one of those cases though, considering the prevalence of this kind of problem in the field.

One thing you've noticed is that the results depend on the time step you choose. When you use ReactorNet.advance, CVODES internally uses an interpolation step to get you to the time that you specified. If, instead, you use ReactorNet.step, CVODES will take a single, internal time step, so you get the maximum resolution out of the solver. This then typically involves a while loop or something:

time = 0

n_steps = 0

while reac.T < XXXX and time < 0.1 and n_steps < 10000:

    time = netw.step()

    n_steps += 1

Realistically, you may not need all those conditions on the while loop, I just put them there to show a few different options that you have. Really, you can be as creative as you want with that condition (derivatives, etc.).

Without seeing your homemade ODE solver, it's difficult to do a direct comparison. What was the step size used in your own ODE solver when the discrepancy was present? Was it the same as the step size for advance?

Best,

Bryan

[Marc Lüttingen]

Thanks for your answer Bryan!

I'm conducting further tests thanks to your useful pieces of information and will post the results here.

Cheers.

--
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/8de1da26-190f-469c-b044-bfab55fbc2f0n%40googlegroups.com.