1 #include "gtest/gtest.h"
2 #include "cantera/thermo.h"
3 #include "cantera/kinetics.h"
4 #include "cantera/thermo/IdealGasPhase.h"
5 #include "cantera/thermo/SurfPhase.h"
6 #include "cantera/kinetics/GasKinetics.h"
7 #include "cantera/base/Solution.h"
8
9 namespace Cantera
10 {
11
12 #ifndef CT_NO_PYTHON
TEST(FracCoeff,ConvertFracCoeff)13 TEST(FracCoeff, ConvertFracCoeff)
14 {
15 IdealGasPhase thermo1("../data/frac.cti", "gas");
16 std::vector<ThermoPhase*> phases1 { &thermo1 };
17 GasKinetics kinetics1;
18 importKinetics(thermo1.xml(), phases1, &kinetics1);
19
20 IdealGasPhase thermo2("../data/frac.xml", "gas");
21 std::vector<ThermoPhase*> phases2 { &thermo2 };
22 GasKinetics kinetics2;
23 importKinetics(thermo2.xml(), phases2, &kinetics2);
24
25 ASSERT_EQ(thermo2.nSpecies(), thermo1.nSpecies());
26 ASSERT_EQ(kinetics2.nReactions(), kinetics1.nReactions());
27
28 for (size_t i = 0; i < kinetics1.nReactions(); i++) {
29 for (size_t k = 0; k < thermo1.nSpecies(); k++) {
30 EXPECT_EQ(kinetics1.reactantStoichCoeff(k,i),
31 kinetics2.reactantStoichCoeff(k,i));
32 EXPECT_EQ(kinetics1.productStoichCoeff(k,i),
33 kinetics2.productStoichCoeff(k,i));
34 }
35 }
36 }
37 #endif
38
39 class FracCoeffTest : public testing::Test
40 {
41 public:
FracCoeffTest()42 FracCoeffTest() :
43 therm("frac.yaml", "gas")
44 {
45 kin = newKinetics({&therm}, "frac.yaml", "gas");
46 therm.setState_TPX(2000, 4*OneAtm,
47 "H2O:0.5, OH:.05, H:0.1, O2:0.15, H2:0.2");
48 kH2O = therm.speciesIndex("H2O");
49 kH = therm.speciesIndex("H");
50 kOH = therm.speciesIndex("OH");
51 kO2 = therm.speciesIndex("O2");
52 kH2 = therm.speciesIndex("H2");
53 }
54 IdealGasPhase therm;
55 unique_ptr<Kinetics> kin;
56 size_t kH2O, kH, kOH, kO2, kH2;
57 };
58
TEST_F(FracCoeffTest,StoichCoeffs)59 TEST_F(FracCoeffTest, StoichCoeffs)
60 {
61 EXPECT_DOUBLE_EQ(1.0, kin->reactantStoichCoeff(kH2O, 0));
62 EXPECT_DOUBLE_EQ(1.4, kin->productStoichCoeff(kH, 0));
63 EXPECT_DOUBLE_EQ(0.6, kin->productStoichCoeff(kOH, 0));
64 EXPECT_DOUBLE_EQ(0.2, kin->productStoichCoeff(kO2, 0));
65
66 EXPECT_DOUBLE_EQ(0.7, kin->reactantStoichCoeff(kH2, 1));
67 EXPECT_DOUBLE_EQ(0.6, kin->reactantStoichCoeff(kOH, 1));
68 EXPECT_DOUBLE_EQ(0.2, kin->reactantStoichCoeff(kO2, 1));
69 EXPECT_DOUBLE_EQ(1.0, kin->productStoichCoeff(kH2O, 1));
70 }
71
TEST_F(FracCoeffTest,RateConstants)72 TEST_F(FracCoeffTest, RateConstants)
73 {
74 vector_fp kf(kin->nReactions(), 0.0);
75 vector_fp kr(kin->nReactions(), 0.0);
76 kin->getFwdRateConstants(&kf[0]);
77 kin->getRevRateConstants(&kr[0]);
78
79 // sum of reaction orders is 1.0; kf has units of 1/s
80 EXPECT_DOUBLE_EQ(1e13, kf[0]);
81
82 // sum of reaction orders is 3.8.
83 // kf = 1e13 (mol/cm^3)^-2.8 s^-1 = 1e13*1000^-2.8 (kmol/m^3)^-2.8 s^-1
84 EXPECT_NEAR(1e13*pow(1e3, -2.8), kf[1], 1e-2);
85
86 // Reactions are irreversible
87 EXPECT_DOUBLE_EQ(0.0, kr[0]);
88 EXPECT_DOUBLE_EQ(0.0, kr[1]);
89 }
90
TEST_F(FracCoeffTest,RatesOfProgress)91 TEST_F(FracCoeffTest, RatesOfProgress)
92 {
93 vector_fp kf(kin->nReactions(), 0.0);
94 vector_fp conc(therm.nSpecies(), 0.0);
95 vector_fp ropf(kin->nReactions(), 0.0);
96 therm.getConcentrations(&conc[0]);
97 kin->getFwdRateConstants(&kf[0]);
98 kin->getFwdRatesOfProgress(&ropf[0]);
99
100 EXPECT_DOUBLE_EQ(conc[kH2O]*kf[0], ropf[0]);
101 EXPECT_DOUBLE_EQ(pow(conc[kH2], 0.8)*conc[kO2]*pow(conc[kOH],2)*kf[1],
102 ropf[1]);
103 }
104
TEST_F(FracCoeffTest,CreationDestructionRates)105 TEST_F(FracCoeffTest, CreationDestructionRates)
106 {
107 vector_fp ropf(kin->nReactions(), 0.0);
108 vector_fp cdot(therm.nSpecies(), 0.0);
109 vector_fp ddot(therm.nSpecies(), 0.0);
110 kin->getFwdRatesOfProgress(&ropf[0]);
111 kin->getCreationRates(&cdot[0]);
112 kin->getDestructionRates(&ddot[0]);
113
114 EXPECT_DOUBLE_EQ(ropf[0], ddot[kH2O]);
115 EXPECT_DOUBLE_EQ(1.4*ropf[0], cdot[kH]);
116 EXPECT_DOUBLE_EQ(0.6*ropf[0], cdot[kOH]);
117 EXPECT_DOUBLE_EQ(0.2*ropf[0], cdot[kO2]);
118
119 EXPECT_DOUBLE_EQ(0.7*ropf[1]+ropf[2], ddot[kH2]);
120 EXPECT_DOUBLE_EQ(0.6*ropf[1], ddot[kOH]);
121 EXPECT_DOUBLE_EQ(0.2*ropf[1]+0.5*ropf[2], ddot[kO2]);
122 EXPECT_DOUBLE_EQ(ropf[1]+ropf[2], cdot[kH2O]);
123
124 EXPECT_DOUBLE_EQ(0.0, cdot[therm.speciesIndex("O")]);
125 EXPECT_DOUBLE_EQ(0.0, ddot[therm.speciesIndex("O")]);
126 }
127
TEST_F(FracCoeffTest,EquilibriumConstants)128 TEST_F(FracCoeffTest, EquilibriumConstants)
129 {
130 vector_fp Kc(kin->nReactions(), 0.0);
131 vector_fp mu0(therm.nSpecies(), 0.0);
132
133 kin->getEquilibriumConstants(&Kc[0]);
134 therm.getGibbs_ref(&mu0[0]); // at pRef
135
136 double deltaG0_0 = 1.4 * mu0[kH] + 0.6 * mu0[kOH] + 0.2 * mu0[kO2] - mu0[kH2O];
137 double deltaG0_1 = mu0[kH2O] - 0.7 * mu0[kH2] - 0.6 * mu0[kOH] - 0.2 * mu0[kO2];
138
139 double pRef = therm.refPressure();
140 double RT = therm.RT();
141
142 // Net stoichiometric coefficients are 1.2 and -0.5
143 EXPECT_NEAR(exp(-deltaG0_0/RT) * pow(pRef/RT, 1.2), Kc[0], 1e-13 * Kc[0]);
144 EXPECT_NEAR(exp(-deltaG0_1/RT) * pow(pRef/RT, -0.5), Kc[1], 1e-13 * Kc[1]);
145 }
146
147 class NegativePreexponentialFactor : public testing::Test
148 {
149 public:
NegativePreexponentialFactor()150 NegativePreexponentialFactor() {}
setup(const std::string & infile)151 void setup(const std::string& infile) {
152 soln = newSolution(infile);
153 soln->thermo()->setState_TPX(2000, OneAtm,
154 "H2O:1.0, H:0.2, O2:0.3, NH:0.05, NO:0.05, N2O:0.05");
155 nSpec = soln->thermo()->nSpecies();
156 nRxn = soln->kinetics()->nReactions();
157 }
158
testNetProductionRates()159 void testNetProductionRates() {
160 const double wdot_ref[] = {0.44705, -0.0021443, 0, -279.36, 0.0021432, 278.92, 0.4449, -279.36, 279.36, 0, 0, 0};
161 ASSERT_EQ(12, (int) nSpec);
162 ASSERT_EQ(12, (int) nRxn);
163 vector_fp wdot(nSpec);
164 soln->kinetics()->getNetProductionRates(&wdot[0]);
165 for (size_t i = 0; i < nSpec; i++) {
166 EXPECT_NEAR(wdot_ref[i], wdot[i], std::abs(wdot_ref[i])*2e-5 + 1e-9);
167 }
168
169 const double ropf_ref[] = {479.305, -128.202, 0, -0, 0, 0, 0, 0, 0, 0.4449, 0, 0};
170 const double ropr_ref[] = {97.94, -26.1964, 0, -0, 1.10334e-06, 0, 0, 0, 6.58592e-06, 0, 0, 0.00214319};
171
172 vector_fp ropf(nRxn);
173 vector_fp ropr(nRxn);
174 soln->kinetics()->getFwdRatesOfProgress(&ropf[0]);
175 soln->kinetics()->getRevRatesOfProgress(&ropr[0]);
176 for (size_t i = 0; i < nRxn; i++) {
177 EXPECT_NEAR(ropf_ref[i], ropf[i], std::abs(ropf_ref[i])*2e-5 + 1e-9);
178 EXPECT_NEAR(ropr_ref[i], ropr[i], std::abs(ropr_ref[i])*2e-5 + 1e-9);
179 }
180 }
181
182 shared_ptr<Solution> soln;
183 size_t nRxn, nSpec;
184 };
185
186 #ifndef CT_NO_PYTHON
TEST_F(NegativePreexponentialFactor,fromCti)187 TEST_F(NegativePreexponentialFactor, fromCti)
188 {
189 setup("../data/noxNeg.cti");
190 testNetProductionRates();
191 }
192 #endif
193
TEST_F(NegativePreexponentialFactor,fromYaml)194 TEST_F(NegativePreexponentialFactor, fromYaml)
195 {
196 setup("noxNeg.yaml");
197 testNetProductionRates();
198 }
199
TEST(InterfaceReaction,CoverageDependency)200 TEST(InterfaceReaction, CoverageDependency) {
201 IdealGasPhase gas("ptcombust.yaml", "gas");
202 SurfPhase surf("ptcombust.yaml", "Pt_surf");
203 shared_ptr<Kinetics> kin(newKinetics({&surf, &gas}, "ptcombust.yaml", "Pt_surf"));
204 ASSERT_EQ(kin->nReactions(), (size_t) 24);
205
206 double T = 500;
207 surf.setState_TP(T, 101325);
208 surf.setCoveragesByName("PT(S):0.7, H(S):0.3");
209 vector_fp kf(kin->nReactions());
210 kin->getFwdRateConstants(&kf[0]);
211 EXPECT_NEAR(kf[0], 4.4579e7 * pow(T, 0.5), 1e-14*kf[0]);
212 // Energies in XML file are converted from J/mol to J/kmol
213 EXPECT_NEAR(kf[1], 3.7e20 * exp(-(67.4e6-6e6*0.3)/(GasConstant*T)), 1e-14*kf[1]);
214 }
215
216 }
217