Obtaining reaction rate parameters via Matlab

398 views
Skip to first unread message

charlesreid1

unread,
Oct 20, 2014, 2:12:32 PM10/20/14
to canter...@googlegroups.com
Hi,

I am using a Matlab API for a third-party program, and would like to populate species and reaction rate information using Cantera's Matlab toolbox. Specifically, I need some way of obtaining forward/reverse Arrhenius constants and activation energies. Is there a way to obtain this information via Cantera's Matlab toolbox? I didn't find anything when I looked through the Matlab toolbox code, but I figured I'd ask on the mailing list.


Charles

Ray Speth

unread,
Oct 20, 2014, 2:52:01 PM10/20/14
to canter...@googlegroups.com
Hi Charles,

There are no functions for providing this information via the Matlab toolbox or any of the other interfaces. The reason for this is, in part, because the modified Arrhenius rate expression is not general, and is not the only reaction rate parameterization supported by Cantera. For example, what would you expect such methods to return for a pressure-dependent "falloff" reaction? One approach you might consider would be to calculate locally-equivalent Arrhenius parameters indirectly by evaluating derivatives of the rate constants. Assume that k(T) has the form:

k(T) = A * T**b * exp(-Ta/T)

where (A,b,Ta) are the parameters for the modified Arrhenius equation. Then let:

    F = 1/k * dk/dT = (b*T + Ta) / T
    G = 1/k * d2k/dT2 = ((b-1)*b*T**2 + 2*(b-1)*T*Ta + Ta**2)/T**4

We can solve these two equations for Ta and b:

    Ta = T**2*(F**2*T- F -G*T)
    b = T*(2*F - F**2*T + G*T)
    A = k / (T**b * exp(-Ta/T))

F and G can be computed by finite difference (here using the Cantera Python module):

T = 1000.0
dT
= 0.1
R
= 1.9872041 # cal/mol-K
gas
.TP = T, 101325
k
= gas.forward_rate_constants
gas
.TP = T+dT, None
kp
= gas.forward_rate_constants
gas
.TP = T-dT, None
km
= gas.forward_rate_constants
F
= (kp-km) / (k*2*dT)
G
= (kp-2*k+km) / (dT**2*k)
Ea = T**2 * (F**2 * T - F - G*T) * R
b
= T * (2*F - F**2*T + G*T)
A
= k * T**-b * np.exp(Ea/(R*T))

Using dT = 0.1 seems to work better than smaller values, and gives about 5 significant figures, at least for the reactions in h2o2.xml. Working out the equivalent Matlab code should be straightforward.

Regards,
Ray

charlesreid1

unread,
Oct 20, 2014, 3:11:58 PM10/20/14
to canter...@googlegroups.com
Thanks Ray, that's a very elegant solution. I'm using an interface that is capable of handling the various Chemkin reaction types, like Cantera, so some way of accessing reaction rate parameters read from Chemkin/CTML files, whatever their form, was more what I was looking for. But this could be helpful as a way of getting towards that. Thanks for the quick response!


Charles
Reply all
Reply to author
Forward
0 new messages