Incorrect rates for a 7,7 degree chebychev formatted rate coef.

[GoogleCodeExporter]

What steps will reproduce the problem?
1. Add this reaction to a .cti file called 'mech' (attached)

chebyshev_reaction('CH3OCH2O2 => CH2OCH2O2H',
                    Tmin=300.0, Tmax=1000.0,
                    Pmin=(0.004,'atm'), Pmax=(1363.0,'atm'),
                    coeffs=[[1.74803, 0.48491, -0.19270, 0.04531, 0.00299, -0.00661, 0.00030],
                            [-3.74931, -0.63409, 0.20864, -0.02864, -0.01588, 0.00962, 0.00221],
                            [0.00982, 0.21793, -0.02120, -0.02227, 0.01705, -0.00258, -0.00548],
                            [0.06149, -0.02515, -0.03822, 0.02761, -0.01039, -0.00312, 0.00731],
                            [-0.01192, -0.01609, 0.02664, -0.01083, -0.00013, 0.00671, -0.00652],
                            [0.00236, 0.00495, -0.00183, -0.00529, 0.00625, -0.00472, 0.00213],
                            [-0.00374, 0.00477, -0.00741, 0.00797, -0.00468, 0.00078, 0.00112]])

2. In short, with g = ct.Solution(mech), and at 1 atm, print out forward rates 
(g.forward_rate_constants) for various T.
Compared to a manual calculation of the rates, you get something like:

CH3OCH2O2 (+ M) <=> CH2OCH2O2H (+ M)
 T, K   k_fwd_cant.     k_actual     
  300  4.54e+05         1.19E-02
  400  1.12e+03         4.88E+00
  500  1.80e+01         2.76E+02
  700  2.34e-01         2.78E+04

3.What is the expected output? What do you see instead?
The expected output is the data of column 'k_actual' in the above table; 
instead, Cantera gives values in the column 'k_fwd_cant' 

What version of the product are you using? On what operating system?
dev, Ubuntu 13.10

Please provide any additional information below.

for the cti file (also, add these species to the species list):

species(name='CH3OCH2O2',
        atoms='H:5 C:2 O:3',
        thermo=(NASA([300.00, 1389.00],
                     [ 2.21029612E+00,  3.68877454E-02, -2.82561555E-05,
                       1.15730533E-08, -1.97130470E-12, -1.94940940E+04,
                       1.91463601E+01]),
                NASA([1389.00, 5000.00],
                     [ 1.24249729E+01,  1.18705986E-02, -4.07906532E-06,
                       6.35310809E-10, -3.69427867E-14, -2.29679238E+04,
                      -3.53740145E+01])),
        transport=gas_transport(geom='nonlinear',
                                diam=4.624,
                                well_depth=329.4,
                                rot_relax=1.0),
        note='7/20/98THERM')

species(name='CH2OCH2O2H',
        atoms='H:5 C:2 O:3',
        thermo=(NASA([300.00, 1393.00],
                     [ 2.52895507E+00,  4.24128290E-02, -3.73406386E-05,
                       1.66639333E-08, -2.96443312E-12, -1.44293306E+04,
                       1.76899251E+01]),
                NASA([1393.00, 5000.00],
                     [ 1.51191783E+01,  9.23718883E-03, -3.19127505E-06,
                       4.99114678E-10, -2.91162488E-14, -1.84114867E+04,
                      -4.85706618E+01])),
        transport=gas_transport(geom='nonlinear',
                                diam=4.624,
                                well_depth=329.4,
                                rot_relax=1.0),
        note='7/20/98THERM')        

Original issue reported on code.google.com by enoch.da...@gmail.com on 24 Oct 2014 at 11:02

Attachments:

• mech.cti

[GoogleCodeExporter]

How are you computing your expected rates? I wrote the following script to 
evaluate the rate constants, and get the same values produced by Cantera:

from __future__ import division

import numpy as np
import scipy.special

C=[[1.74803, 0.48491, -0.19270, 0.04531, 0.00299, -0.00661, 0.00030],
   [-3.74931, -0.63409, 0.20864, -0.02864, -0.01588, 0.00962, 0.00221],
   [0.00982, 0.21793, -0.02120, -0.02227, 0.01705, -0.00258, -0.00548],
   [0.06149, -0.02515, -0.03822, 0.02761, -0.01039, -0.00312, 0.00731],
   [-0.01192, -0.01609, 0.02664, -0.01083, -0.00013, 0.00671, -0.00652],
   [0.00236, 0.00495, -0.00183, -0.00529, 0.00625, -0.00472, 0.00213],
   [-0.00374, 0.00477, -0.00741, 0.00797, -0.00468, 0.00078, 0.00112]]

phi = [scipy.special.chebyt(i) for i in range(7)]

def k(T,P):
    Pmin = 101325.0 * 0.004
    Pmax = 101325.0 * 1363.0
    Tt = (2/T - 1/300.0 - 1/1000.0) / (1/1000.0 - 1/300.0)
    Pt = (2*np.log10(P) - np.log10(Pmin) - np.log10(Pmax)) / (np.log10(Pmax) - np.log10(Pmin))

    assert -1 <= Tt <= 1
    assert -1 <= Pt <= 1

    logk = 0.0
    for t in range(7):
        for p in range(7):
            logk += C[t][p] * phi[t](Tt) * phi[p](Pt)

    return 10**logk

for T in [300, 400, 500, 700]:
    print 'T = {:.0f}, k = {:.3e}'.format(T, k(T, 101325))

Original comment by yarmond on 27 Oct 2014 at 8:54

• Changed state: Feedback

[GoogleCodeExporter]

Hello, I've confirmed that the error was in fact in the publication (their eqn 
for calculating the rates, and that for whatever reason Tmin and Tmax values 
must be switched to get agreeable rates) from which I obtained the chebyshev 
coefficients. sorry for the misunderstanding. 

Original comment by enoch.da...@gmail.com on 5 Nov 2014 at 3:40

[GoogleCodeExporter]

Original comment by yarmond on 5 Nov 2014 at 3:55

• Changed state: Invalid