search

jcaiuwyo / cantera

Automatically exported from code.google.com/p/cantera
0 stars 0 forks source link

pdep_arrhenius rate expressions with 0 A coefficients give 'nan' rate constants (but don't raise an exception) #98

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago 
The master chemical mechanism developed by Peter Glarborg and coworkers at the 
Technical University of Denmark contains several rate expressions for CO + OH = 
HOCO  for use at different pressures.

The tricky thing is that at low pressures, Glarborg et al use only one 
Arrhenius expression to represent the rate,
at 10 bar they use the sum of two, and at 20+ bar they use the sum of three 
Arrhenius expressions.

For use in Chemkin (and RMG) I translated these into a series of duplicate PLOG 
expressions:

// (0.001-2000 bar, 300<T<2000K)
CO + OH = HOCO              5.5E44 -11.000    7948       0.0 0.0 0.0
  PLOG  /         0.001       1.0E25  -6.000    2981  /
  PLOG  /            1      6.0E26  -5.600    2881    /
  PLOG  /          3      2.2E27  -5.600    3239      /
  PLOG  /          10     1.5E25  -5.000    1987      /
  PLOG  /           20      4.2E26  -5.700    1927    /
  PLOG  /          50      4.9E25  -5.200    1987     /
  PLOG  /           80      5.2E25  -5.200    1987    /
  PLOG  /          100      1.1E28  -6.000    2384    /
  PLOG  /          650      3.2E41 -10.000    6955    /
  PLOG  /          2000      5.5E44 -11.000    7948   /
DUPLICATE                                             

// (0.001 bar, 300<T<2000K)
CO + OH = HOCO                  2.7E67 -17.000   22851      0.0 0.0 0.0
  PLOG  /           0.001       0.0      0.0       0  /
  PLOG  /            1        0.0      0.0       0    /
  PLOG  /           3         0.0      0.0       0    /
  PLOG  /            10      1.3E37  -8.400    7948   /
  PLOG  /            20      7.5E28  -6.000    3775   /
  PLOG  /            50     4.0E38  -9.000    6955    /
  PLOG  /          80       7.7E35  -8.000    6557    /
  PLOG  /          100       1.8E36  -8.000    7153   /
  PLOG  /          650       2.9E66 -17.100   19870   /
  PLOG  /          2000       2.7E67 -17.000   22851  /
DUPLICATE                                               

// (0.001 bar, 300<T<2000K)
CO + OH = HOCO                 1.0E74 -18.000   37157      0.0 0.0 0.0
  PLOG  /            0.001        0.0      0.0       0    /
  PLOG  /             1         0.0      0.0       0      /
  PLOG  /             3         0.0      0.0       0      /
  PLOG  /             10        0.0      0.0       0      /
  PLOG  /             20       4.0E39  -9.000    9935     /
  PLOG  /              50      5.0E43 -10.000   13015     /
  PLOG  /             80       9.0E47 -11.200   15499     /
  PLOG  /            100       2.0E54 -13.000   19671     /
  PLOG  /            650        2.0E63 -15.200   27421    /
  PLOG  /              2000      1.0E74 -18.000   37157   /
DUPLICATE                                    

However, when I try this in Cantera :

# Reaction 378
pdep_arrhenius('CO(15) + OH(26) <=> HOCO(30)',
               [(0.001, 'atm'), 1.000000e+25, -6.0, 2.981],
               [(1.0, 'atm'), 6.000000e+26, -5.6, 2.881],
               [(3.0, 'atm'), 2.200000e+27, -5.6, 3.239],
               [(10.0, 'atm'), 1.500000e+25, -5.0, 1.987],
               [(20.0, 'atm'), 4.200000e+26, -5.7, 1.927],
               [(50.0, 'atm'), 4.900000e+25, -5.2, 1.987],
               [(80.0, 'atm'), 5.200000e+25, -5.2, 1.987],
               [(100.0, 'atm'), 1.100000e+28, -6.0, 2.384],
               [(650.0, 'atm'), 3.200000e+41, -10.0, 6.955],
               [(2000.0, 'atm'), 5.500000e+44, -11.0, 7.948],
               options='duplicate')

# Reaction 379
# THIS REACTION GIVES A NAN FWD RATE CONSTANT
pdep_arrhenius('CO(15) + OH(26) <=> HOCO(30)',
               [(0.001, 'atm'), 0.000000e+00, 0.0, 0.0],
               [(1.0, 'atm'), 0.000000e+00, 0.0, 0.0],
               [(3.0, 'atm'), 0.000000e+00, 0.0, 0.0],
               [(10.0, 'atm'), 1.300000e+37, -8.4, 7.948],
               [(20.0, 'atm'), 7.500000e+28, -6.0, 3.775],
               [(50.0, 'atm'), 4.000000e+38, -9.0, 6.955],
               [(80.0, 'atm'), 7.700000e+35, -8.0, 6.557],
               [(100.0, 'atm'), 1.800000e+36, -8.0, 7.153],
               [(650.0, 'atm'), 2.900000e+66, -17.1, 19.87],
               [(2000.0, 'atm'), 2.700000e+67, -17.0, 22.851],
               options='duplicate')

# Reaction 380
# THIS REACTION GIVES A NAN FWD RATE CONSTANT
pdep_arrhenius('CO(15) + OH(26) <=> HOCO(30)',
               [(0.001, 'atm'), 0.000000e+00, 0.0, 0.0],
               [(1.0, 'atm'), 0.000000e+00, 0.0, 0.0],
               [(3.0, 'atm'), 0.000000e+00, 0.0, 0.0],
               [(10.0, 'atm'), 0.000000e+00, 0.0, 0.0],
               [(20.0, 'atm'), 4.000000e+39, -9.0, 9.935],
               [(50.0, 'atm'), 5.000000e+43, -10.0, 13.015],
               [(80.0, 'atm'), 9.000000e+47, -11.2, 15.499],
               [(100.0, 'atm'), 2.000000e+54, -13.0, 19.671],
               [(650.0, 'atm'), 2.000000e+63, -15.2, 27.421],
               [(2000.0, 'atm'), 1.000000e+74, -18.0, 37.157],
               options='duplicate')

This loads fine, but if you evaluate the rate coefficients at 1 atm, it does 
not complain with "Procedure: getRateCoefficient
Error:   negative or zero A coefficient for reaction __" (like normal arrhenius 
expressions) but instead just gives "NAN" as the rate constant, which 
propagates through everything else and you end up with NANs everywhere.

Cantera 2.0.0b3 on MacOS X accessed via Python 2.7.3

Original issue reported on code.google.com by r.h.w...@gmail.com on 18 Jul 2012 at 7:15

GoogleCodeExporter commented 9 years ago 
You can give multiple rate expressions at the same pressure, which will be 
added before doing the logarithmic interpolation in pressure, e.g.:

pdep_arrhenius('CO(15) + OH(26) <=> HOCO(30)',
               [(0.001, 'atm'), 1.000000e+25, -6.0, 2.981],
               [(1.0, 'atm'), 6.000000e+26, -5.6, 2.881],
               [(3.0, 'atm'), 2.200000e+27, -5.6, 3.239],
               [(10.0, 'atm'), 1.500000e+25, -5.0, 1.987],
               [(10.0, 'atm'), 1.300000e+37, -8.4, 7.948],
               [(20.0, 'atm'), 4.200000e+26, -5.7, 1.927],
               [(20.0, 'atm'), 7.500000e+28, -6.0, 3.775],
               [(20.0, 'atm'), 4.000000e+39, -9.0, 9.935],
               [(50.0, 'atm'), 4.900000e+25, -5.2, 1.987],
               [(50.0, 'atm'), 4.000000e+38, -9.0, 6.955],
               [(50.0, 'atm'), 5.000000e+43, -10.0, 13.015],
               [(80.0, 'atm'), 5.200000e+25, -5.2, 1.987],
               [(80.0, 'atm'), 7.700000e+35, -8.0, 6.557],
               [(80.0, 'atm'), 9.000000e+47, -11.2, 15.499],
               [(100.0, 'atm'), 1.100000e+28, -6.0, 2.384],
               [(100.0, 'atm'), 1.800000e+36, -8.0, 7.153],
               [(100.0, 'atm'), 2.000000e+54, -13.0, 19.671],
               [(650.0, 'atm'), 3.200000e+41, -10.0, 6.955],
               [(650.0, 'atm'), 2.900000e+66, -17.1, 19.87],
               [(650.0, 'atm'), 2.000000e+63, -15.2, 27.421],
               [(2000.0, 'atm'), 5.500000e+44, -11.0, 7.948],
               [(2000.0, 'atm'), 2.700000e+67, -17.0, 22.851],
               [(2000.0, 'atm'), 1.000000e+74, -18.0, 37.157])

This has the advantage of also working when some of the Arrhenius expressions 
at a particular pressure have negative A factors. I haven't worked out a 
general way to detect whether a particular combination of rate expressions will 
lead to a negative rate coefficient for some value of T (and subsequently break 
the logarithmic interpolation), but checking that A != 0.0 wouldn't hurt.

Original comment by yarmond on 19 Jul 2012 at 1:16

GoogleCodeExporter commented 9 years ago 
Fixed in r1713 in the 2.0 maintenance branch; will be merged back into trunk 
later.

Original comment by yarmond on 24 Jul 2012 at 11:04