research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Kinetics] Fix Chebyshev rate evaluation with only 1 point in T or P

Browse source
This commit is contained in:
Ray Speth 2019-09-13 14:24:07 -04:00
parent b44f189569
commit ecce98c1dc
2 changed files with 49 additions and 8 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

16
include/cantera/kinetics/RxnRates.h
View file

@ -364,13 +364,13 @@ public:
//! @param c base-10 logarithm of the pressure in Pa
void update_C(const doublereal* c) {
double Pr = (2 * c[0] + PrNum_) * PrDen_;
double Cnm1 = 1;
double Cn = Pr;
double Cnm1 = Pr;
double Cn = 1;
double Cnp1;
for (size_t j = 0; j < nT_; j++) {
dotProd_[j] = chebCoeffs_[nP_*j] + Pr * chebCoeffs_[nP_*j+1];
dotProd_[j] = chebCoeffs_[nP_*j];
}
for (size_t i = 2; i < nP_; i++) {
for (size_t i = 1; i < nP_; i++) {
Cnp1 = 2 * Pr * Cn - Cnm1;
for (size_t j = 0; j < nT_; j++) {
dotProd_[j] += Cnp1 * chebCoeffs_[nP_*j + i];
@ -387,11 +387,11 @@ public:
*/
doublereal updateRC(doublereal logT, doublereal recipT) const {
double Tr = (2 * recipT + TrNum_) * TrDen_;
double Cnm1 = 1;
double Cn = Tr;
double Cnm1 = Tr;
double Cn = 1;
double Cnp1;
double logk = dotProd_[0] + Tr * dotProd_[1];
for (size_t i = 2; i < nT_; i++) {
double logk = dotProd_[0];
for (size_t i = 1; i < nT_; i++) {
Cnp1 = 2 * Tr * Cn - Cnm1;
logk += Cnp1 * dotProd_[i];
Cnm1 = Cn;

41
interfaces/cython/cantera/test/test_kinetics.py
View file

@ -943,6 +943,47 @@ class TestReaction(utilities.CanteraTest):
self.assertNear(gas2.net_rates_of_progress[0],
gas1.net_rates_of_progress[4])
def test_chebyshev_single_P(self):
species = ct.Species.listFromFile('pdep-test.cti')
r = ct.ChebyshevReaction()
r.reactants = 'R5:1, H:1'
r.products = 'P5A:1, P5B:1'
r.set_parameters(Tmin=300.0, Tmax=2000.0, Pmin=1000, Pmax=10000000,
coeffs=[[ 5.28830e+00],
[ 1.97640e+00],
[ 3.17700e-01],
[-3.12850e-02]])
gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics',
species=species, reactions=[r])
# rate constant should be pressure independent
for T in [300, 500, 1500]:
gas.TP = T, 1e4
k1 = gas.forward_rate_constants[0]
gas.TP = T, 1e6
k2 = gas.forward_rate_constants[0]
self.assertNear(k1, k2)
def test_chebyshev_single_T(self):
species = ct.Species.listFromFile('pdep-test.cti')
r = ct.ChebyshevReaction()
r.reactants = 'R5:1, H:1'
r.products = 'P5A:1, P5B:1'
r.set_parameters(Tmin=300.0, Tmax=2000.0, Pmin=1000, Pmax=10000000,
coeffs=[[ 5.28830e+00, -1.13970e+00, -1.20590e-01, 1.60340e-02]])
gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics',
species=species, reactions=[r])
# rate constant should be temperature independent
for P in [1e4, 2e5, 8e6]:
gas.TP = 400, P
k1 = gas.forward_rate_constants[0]
gas.TP = 1700, P
k2 = gas.forward_rate_constants[0]
self.assertNear(k1, k2)
def test_chebyshev_rate(self):
gas1 = ct.Solution('pdep-test.cti')
gas1.TP = 800, 2*ct.one_atm