1 #include "gtest/gtest.h"
2 #include "cantera/thermo/PengRobinson.h"
3 #include "cantera/thermo/ThermoFactory.h"
4
5
6 namespace Cantera
7 {
8
9 class PengRobinson_Test : public testing::Test
10 {
11 public:
PengRobinson_Test()12 PengRobinson_Test() {
13 test_phase.reset(newPhase("../data/thermo-models.yaml", "CO2-PR"));
14 }
15
16 //vary the composition of a co2-h2 mixture:
set_r(const double r)17 void set_r(const double r) {
18 vector_fp moleFracs(7);
19 moleFracs[0] = r;
20 moleFracs[2] = 1-r;
21 test_phase->setMoleFractions(&moleFracs[0]);
22 }
23
24 std::unique_ptr<ThermoPhase> test_phase;
25 };
26
TEST_F(PengRobinson_Test,construct_from_yaml)27 TEST_F(PengRobinson_Test, construct_from_yaml)
28 {
29 PengRobinson* peng_robinson_phase = dynamic_cast<PengRobinson*>(test_phase.get());
30 EXPECT_TRUE(peng_robinson_phase != NULL);
31 }
32
TEST_F(PengRobinson_Test,chem_potentials)33 TEST_F(PengRobinson_Test, chem_potentials)
34 {
35 test_phase->setState_TP(298.15, 101325.);
36 /* Chemical potential should increase with increasing co2 mole fraction:
37 * mu = mu_0 + RT ln(gamma_k*X_k).
38 * where gamma_k is the activity coefficient. Run regression test against values
39 * calculated using the model.
40 */
41 const double expected_result[9] = {
42 -457338129.70445037,
43 -457327078.87912911,
44 -457317214.31077951,
45 -457308354.65227401,
46 -457300353.74028891,
47 -457293092.45485628,
48 -457286472.73969948,
49 -457280413.14238912,
50 -457274845.44186872
51 };
52
53 double xmin = 0.6;
54 double xmax = 0.9;
55 int numSteps = 9;
56 double dx = (xmax-xmin)/(numSteps-1);
57 vector_fp chemPotentials(7);
58 for(int i=0; i < numSteps; ++i)
59 {
60 set_r(xmin + i*dx);
61 test_phase->getChemPotentials(&chemPotentials[0]);
62 EXPECT_NEAR(expected_result[i], chemPotentials[0], 1.e-6);
63 }
64 }
65
TEST_F(PengRobinson_Test,chemPotentials_RT)66 TEST_F(PengRobinson_Test, chemPotentials_RT)
67 {
68 double T = 410.0;
69 test_phase->setState_TP(T, 130 * OneAtm);
70
71 // Test that chemPotentials_RT*RT = chemPotentials
72 const double RT = GasConstant * T;
73 vector_fp mu(7);
74 vector_fp mu_RT(7);
75 double xmin = 0.6;
76 double xmax = 0.9;
77 int numSteps = 9;
78 double dx = (xmax-xmin)/(numSteps-1);
79
80 for(int i=0; i < numSteps; ++i)
81 {
82 const double r = xmin + i*dx;
83 set_r(r);
84 test_phase->getChemPotentials(&mu[0]);
85 test_phase->getChemPotentials_RT(&mu_RT[0]);
86 EXPECT_NEAR(mu[0], mu_RT[0]*RT, 1.e-6);
87 EXPECT_NEAR(mu[2], mu_RT[2]*RT, 1.e-6);
88 }
89 }
90
TEST_F(PengRobinson_Test,activityCoeffs)91 TEST_F(PengRobinson_Test, activityCoeffs)
92 {
93 double T = 330.0;
94 test_phase->setState_TP(T, 120 * OneAtm);
95
96 // Test that mu0 + RT log(activityCoeff * MoleFrac) == mu
97 const double RT = GasConstant * T;
98 vector_fp mu0(7);
99 vector_fp activityCoeffs(7);
100 vector_fp chemPotentials(7);
101 double xmin = 0.6;
102 double xmax = 0.9;
103 int numSteps = 9;
104 double dx = (xmax-xmin)/(numSteps-1);
105
106 for(int i=0; i < numSteps; ++i)
107 {
108 const double r = xmin + i*dx;
109 set_r(r);
110 test_phase->getChemPotentials(&chemPotentials[0]);
111 test_phase->getActivityCoefficients(&activityCoeffs[0]);
112 test_phase->getStandardChemPotentials(&mu0[0]);
113 EXPECT_NEAR(chemPotentials[0], mu0[0] + RT*std::log(activityCoeffs[0] * r), 1.e-6);
114 EXPECT_NEAR(chemPotentials[2], mu0[2] + RT*std::log(activityCoeffs[2] * (1-r)), 1.e-6);
115 }
116 }
117
TEST_F(PengRobinson_Test,standardConcentrations)118 TEST_F(PengRobinson_Test, standardConcentrations)
119 {
120 EXPECT_DOUBLE_EQ(test_phase->pressure()/(test_phase->temperature()*GasConstant),
121 test_phase->standardConcentration(0));
122 EXPECT_DOUBLE_EQ(test_phase->pressure()/(test_phase->temperature()*GasConstant),
123 test_phase->standardConcentration(1));
124 }
125
TEST_F(PengRobinson_Test,activityConcentrations)126 TEST_F(PengRobinson_Test, activityConcentrations)
127 {
128 // Check to make sure activityConcentration_i == standardConcentration_i * gamma_i * X_i
129 vector_fp standardConcs(7);
130 vector_fp activityCoeffs(7);
131 vector_fp activityConcentrations(7);
132 double xmin = 0.6;
133 double xmax = 0.9;
134 int numSteps = 9;
135 double dx = (xmax-xmin)/(numSteps-1);
136 test_phase->setState_TP(350, 100 * OneAtm);
137
138 for(int i=0; i < numSteps; ++i)
139 {
140 const double r = xmin + i*dx;
141 set_r(r);
142 test_phase->getActivityCoefficients(&activityCoeffs[0]);
143 standardConcs[0] = test_phase->standardConcentration(0);
144 standardConcs[2] = test_phase->standardConcentration(2);
145 test_phase->getActivityConcentrations(&activityConcentrations[0]);
146
147 EXPECT_NEAR(standardConcs[0] * r * activityCoeffs[0], activityConcentrations[0], 1.e-6);
148 EXPECT_NEAR(standardConcs[2] * (1-r) * activityCoeffs[2], activityConcentrations[2], 1.e-6);
149 }
150 }
151
TEST_F(PengRobinson_Test,setTP)152 TEST_F(PengRobinson_Test, setTP)
153 {
154 // Check to make sure that the phase diagram is accurately reproduced for a few select isobars
155
156 // All sub-cooled liquid:
157 const double rho1[6] = {
158 6.6603507723749249e+002,
159 1.6824762614489907e+002,
160 1.6248354709581241e+002,
161 1.5746729362032696e+002,
162 1.5302217175386241e+002,
163 1.4902908974486667e+002
164 };
165 // Phase change between temperatures 4 & 5:
166 const double rho2[6] = {
167 7.5732259810273172e+002,
168 7.2766981078381912e+002,
169 6.935475475396446e+002,
170 6.5227027102964917e+002,
171 5.9657442842753153e+002,
172 3.9966973266966875e+002
173 };
174 // Supercritical; no discontinuity in rho values:
175 const double rho3[6] = {
176 8.0601205067780199e+002,
177 7.8427655940884574e+002,
178 7.6105347579146576e+002,
179 7.3605202492828505e+002,
180 7.0887891410210011e+002,
181 6.7898591969734434e+002
182 };
183
184 for(int i=0; i<6; ++i)
185 {
186 const double temp = 294 + i*2;
187 set_r(0.999);
188 test_phase->setState_TP(temp, 5542027.5);
189 EXPECT_NEAR(test_phase->density(),rho1[i],1.e-8);
190
191 test_phase->setState_TP(temp, 7388370.);
192 EXPECT_NEAR(test_phase->density(),rho2[i],1.e-8);
193
194 test_phase->setState_TP(temp, 9236712.5);
195 EXPECT_NEAR(test_phase->density(),rho3[i],1.e-8);
196 }
197 }
198
TEST_F(PengRobinson_Test,getPressure)199 TEST_F(PengRobinson_Test, getPressure)
200 {
201 // Check to make sure that the P-R equation is accurately reproduced for a few selected values
202
203 /* This test uses CO2 as the only species.
204 * Values of a_coeff, b_coeff are calculated based on the the critical temperature
205 * and pressure values of CO2 as follows:
206 * a_coeff = 0.457235(RT_crit)^2/p_crit
207 * b_coeff = 0.077796(RT_crit)/p_crit
208 * The temperature dependent parameter in P-R EoS is calculated as
209 * \alpha = [1 + \kappa(1 - sqrt{T/T_crit}]^2
210 * kappa is a function calulated based on the accentric factor.
211 */
212
213 double a_coeff = 3.958095109E+5;
214 double b_coeff = 26.62616317/1000;
215 double acc_factor = 0.228;
216 double pres_theoretical, kappa, alpha, mv;
217 const double rho = 1.0737;
218
219 //Calculate kappa value
220 kappa = 0.37464 + 1.54226*acc_factor - 0.26992*acc_factor*acc_factor;
221
222 for (int i = 0; i < 15; i++)
223 {
224 const double temp = 296 + i * 50;
225 set_r(1.0);
226 test_phase->setState_TR(temp, rho);
227 const double Tcrit = test_phase->critTemperature();
228 mv = 1 / rho * test_phase->meanMolecularWeight();
229 //Calculate pressure using Peng-Robinson EoS
230 alpha = pow(1 + kappa*(1 - sqrt(temp / Tcrit)), 2.0);
231 pres_theoretical = GasConstant*temp / (mv - b_coeff)
232 - a_coeff*alpha / (mv*mv + 2*b_coeff*mv - b_coeff*b_coeff);
233 EXPECT_NEAR(test_phase->pressure(), pres_theoretical, 3);
234 }
235 }
236
TEST_F(PengRobinson_Test,gibbsEnergy)237 TEST_F(PengRobinson_Test, gibbsEnergy)
238 {
239 // Test that g == h - T*s
240 const double T = 360.;
241 double xmin = 0.6;
242 double xmax = 0.9;
243 int numSteps = 9;
244 double dx = (xmax - xmin) / (numSteps - 1);
245 double gibbs_theoretical;
246
247 for (int i = 0; i < numSteps; ++i)
248 {
249 const double r = xmin + i * dx;
250 test_phase->setState_TP(T, 150e5);
251 set_r(r);
252 gibbs_theoretical = test_phase->enthalpy_mole() - T * (test_phase->entropy_mole());
253 EXPECT_NEAR(test_phase->gibbs_mole(), gibbs_theoretical, 1.e-6);
254 }
255 }
256
TEST_F(PengRobinson_Test,totalEnthalpy)257 TEST_F(PengRobinson_Test, totalEnthalpy)
258 {
259 // Test that hbar = \sum (h_k*x_k)
260 double hbar, sum = 0.0;
261 vector_fp partialEnthalpies(7);
262 vector_fp moleFractions(7);
263 double xmin = 0.6;
264 double xmax = 0.9;
265 int numSteps = 9;
266 double dx = (xmax - xmin) / (numSteps - 1);
267
268 for (int i = 0; i < numSteps; ++i)
269 {
270 sum = 0.0;
271 const double r = xmin + i * dx;
272 test_phase->setState_TP(430., 120e5);
273 set_r(r);
274 hbar = test_phase->enthalpy_mole();
275 test_phase->getMoleFractions(&moleFractions[0]);
276 test_phase->getPartialMolarEnthalpies(&partialEnthalpies[0]);
277 for (int k = 0; k < 7; k++)
278 {
279 sum += moleFractions[k] * partialEnthalpies[k];
280 }
281 EXPECT_NEAR(hbar, sum, 1.e-6);
282 }
283 }
284
TEST_F(PengRobinson_Test,cpValidate)285 TEST_F(PengRobinson_Test, cpValidate)
286 {
287 // Test that cp = dH/dT at constant pressure using finite difference method
288
289 double p = 200e5;
290 double Tmin = 298;
291 int numSteps = 20;
292 double dT = 1e-4;
293 test_phase->setMoleFractionsByName("CO2: 0.7, H2O: 0.1, H2: 0.2");
294
295 for (int i = 0; i < numSteps; ++i) {
296 const double T = Tmin + 10 * i;
297 test_phase->setState_TP(T - dT, p);
298 double h1 = test_phase->enthalpy_mole(); // J/kmol
299 test_phase->setState_TP(T, p);
300 double cp = test_phase->cp_mole(); // unit is J/kmol/K
301 test_phase->setState_TP(T + dT, p);
302 double h2 = test_phase->enthalpy_mole();
303
304 double dh_dT = (h2 - h1) / (2 * dT);
305 EXPECT_NEAR(cp, dh_dT, 1e-6 * cp);
306 }
307 }
308
TEST_F(PengRobinson_Test,CoolPropValidate)309 TEST_F(PengRobinson_Test, CoolPropValidate)
310 {
311 // Validate the P-R EoS in Cantera with P-R EoS from CoolProp
312
313 const double rhoCoolProp[10] = {
314 9.067928191884574,
315 8.318322900591179,
316 7.6883521740498155,
317 7.150504298001246,
318 6.685330199667018,
319 6.278630757480957,
320 5.919763108091383,
321 5.600572727499541,
322 5.314694056926007,
323 5.057077678380463
324 };
325
326 double p = 5e5;
327
328 // Calculate density using Peng-Robinson EoS from Cantera
329 for(int i=0; i<10; i++)
330 {
331 const double temp = 300 + i*25;
332 set_r(1.0);
333 test_phase->setState_TP(temp, p);
334 EXPECT_NEAR(test_phase->density(),rhoCoolProp[i],1.e-5);
335 }
336 }
337
TEST_F(PengRobinson_Test,partialMolarPropertyIdentities)338 TEST_F(PengRobinson_Test, partialMolarPropertyIdentities)
339 {
340 // unique_ptr<ThermoPhase> phase(newPhase("co2_PR_example.yaml"));
341 vector_fp hk(test_phase->nSpecies());
342 vector_fp uk(test_phase->nSpecies());
343 vector_fp sk(test_phase->nSpecies());
344 vector_fp gk(test_phase->nSpecies());
345 vector_fp vk(test_phase->nSpecies());
346 vector_fp X(test_phase->nSpecies());
347
348 test_phase->setMoleFractionsByName("CO2: 0.7, H2O: 0.1, H2: 0.2");
349 test_phase->getMoleFractions(X.data());
350 double P = 100 * OneAtm;
351
352 for (int i = 0; i < 5; i++) {
353 double T = 300 + i*60;
354 test_phase->setState_TP(T, P);
355 test_phase->getPartialMolarEnthalpies(hk.data());
356 test_phase->getPartialMolarIntEnergies(uk.data());
357 test_phase->getPartialMolarEntropies(sk.data());
358 test_phase->getChemPotentials(gk.data());
359 test_phase->getPartialMolarVolumes(vk.data());
360
361 double h_mix = test_phase->enthalpy_mole();
362 double u_mix = test_phase->intEnergy_mole();
363 double s_mix = test_phase->entropy_mole();
364 double g_mix = test_phase->gibbs_mole();
365 double v_mix = test_phase->molarVolume();
366
367 double h = dot(X.begin(), X.end(), hk.begin());
368 EXPECT_NEAR(h, h_mix, 1e-11 * std::abs(h_mix));
369 double u = dot(X.begin(), X.end(), uk.begin());
370 EXPECT_NEAR(u, u_mix, 1e-11 * std::abs(u_mix));
371 double s = dot(X.begin(), X.end(), sk.begin());
372 EXPECT_NEAR(s, s_mix, 1e-11 * std::abs(s_mix));
373 double g = dot(X.begin(), X.end(), gk.begin());
374 EXPECT_NEAR(g, g_mix, 1e-11 * std::abs(g_mix));
375 double v = dot(X.begin(), X.end(), vk.begin());
376 EXPECT_NEAR(v, v_mix, 1e-11 * std::abs(v_mix));
377
378 for (size_t k = 0; k < test_phase->nSpecies(); k++) {
379 EXPECT_NEAR(uk[k] + P * vk[k], hk[k], 1e-11 * std::abs(h_mix));
380 EXPECT_NEAR(hk[k] - T * sk[k], gk[k], 1e-11 * std::abs(g_mix));
381 }
382 }
383 }
384
TEST(PengRobinson,lookupSpeciesProperties)385 TEST(PengRobinson, lookupSpeciesProperties)
386 {
387 AnyMap phase_def = AnyMap::fromYamlString(
388 "{name: test, species: [{gri30.yaml/species: [CO2, CH4, N2]}],"
389 " thermo: Peng-Robinson}"
390 );
391 unique_ptr<ThermoPhase> test(newPhase(phase_def));
392
393 // Check for correspondence to values in critical properties "database"
394 test->setState_TPX(330, 100 * OneAtm, "CH4: 1.0");
395 EXPECT_NEAR(test->critPressure(), 4.63e+06, 1e-2);
396 EXPECT_NEAR(test->critTemperature(), 190.7, 1e-3);
397
398 test->setState_TPX(330, 180 * OneAtm, "N2: 1.0");
399 EXPECT_NEAR(test->critPressure(), 3.39e+06, 1e-2);
400 EXPECT_NEAR(test->critTemperature(), 126.2, 1e-3);
401 }
402
TEST(PengRobinson,lookupSpeciesPropertiesMissing)403 TEST(PengRobinson, lookupSpeciesPropertiesMissing)
404 {
405 AnyMap phase_def = AnyMap::fromYamlString(
406 "{name: test, species: [{gri30.yaml/species: [CO2, CH3, N2]}],"
407 " thermo: Peng-Robinson}"
408 );
409
410 // CH3 is not in the critical properties database, so this should be
411 // detected as an error
412 EXPECT_THROW(newPhase(phase_def), CanteraError);
413 }
414
415 };
416