Cantera 2.6.0
PengRobinson.cpp
Go to the documentation of this file.
1//! @file PengRobinson.cpp
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
10#include "cantera/base/ctml.h"
11
12#include <boost/math/tools/roots.hpp>
13
14namespace bmt = boost::math::tools;
15
16namespace Cantera
17{
18
19const double PengRobinson::omega_a = 4.5723552892138218E-01;
20const double PengRobinson::omega_b = 7.77960739038885E-02;
21const double PengRobinson::omega_vc = 3.07401308698703833E-01;
22
23PengRobinson::PengRobinson(const std::string& infile, const std::string& id_) :
24 m_b(0.0),
25 m_a(0.0),
26 m_aAlpha_mix(0.0),
27 m_NSolns(0),
28 m_dpdV(0.0),
29 m_dpdT(0.0)
30{
31 std::fill_n(m_Vroot, 3, 0.0);
32 initThermoFile(infile, id_);
33}
34
35void PengRobinson::setSpeciesCoeffs(const std::string& species, double a, double b,
36 double w)
37{
38 size_t k = speciesIndex(species);
39 if (k == npos) {
40 throw CanteraError("PengRobinson::setSpeciesCoeffs",
41 "Unknown species '{}'.", species);
42 }
43
44 // Calculate value of kappa (independent of temperature)
45 // w is an acentric factor of species
46 if (w <= 0.491) {
47 m_kappa[k] = 0.37464 + 1.54226*w - 0.26992*w*w;
48 } else {
49 m_kappa[k] = 0.374642 + 1.487503*w - 0.164423*w*w + 0.016666*w*w*w;
50 }
51 m_acentric[k] = w; // store the original acentric factor to enable serialization
52
53 // Calculate alpha (temperature dependent interaction parameter)
54 double critTemp = speciesCritTemperature(a, b);
55 double sqt_T_r = sqrt(temperature() / critTemp);
56 double sqt_alpha = 1 + m_kappa[k] * (1 - sqt_T_r);
57 m_alpha[k] = sqt_alpha*sqt_alpha;
58
59 m_a_coeffs(k,k) = a;
60 double aAlpha_k = a*m_alpha[k];
61 m_aAlpha_binary(k,k) = aAlpha_k;
62
63 // standard mixing rule for cross-species interaction term
64 for (size_t j = 0; j < m_kk; j++) {
65 if (k == j) {
66 continue;
67 }
68 double a0kj = sqrt(m_a_coeffs(j,j) * a);
69 double aAlpha_j = a*m_alpha[j];
70 double a_Alpha = sqrt(aAlpha_j*aAlpha_k);
71 if (m_a_coeffs(j, k) == 0) {
72 m_a_coeffs(j, k) = a0kj;
73 m_aAlpha_binary(j, k) = a_Alpha;
74 m_a_coeffs(k, j) = a0kj;
75 m_aAlpha_binary(k, j) = a_Alpha;
76 }
77 }
78 m_b_coeffs[k] = b;
79}
80
81void PengRobinson::setBinaryCoeffs(const std::string& species_i,
82 const std::string& species_j, double a0)
83{
84 size_t ki = speciesIndex(species_i);
85 if (ki == npos) {
86 throw CanteraError("PengRobinson::setBinaryCoeffs",
87 "Unknown species '{}'.", species_i);
88 }
89 size_t kj = speciesIndex(species_j);
90 if (kj == npos) {
91 throw CanteraError("PengRobinson::setBinaryCoeffs",
92 "Unknown species '{}'.", species_j);
93 }
94
95 m_a_coeffs(ki, kj) = m_a_coeffs(kj, ki) = a0;
96 m_binaryParameters[species_i][species_j] = a0;
97 m_binaryParameters[species_j][species_i] = a0;
98 // Calculate alpha_ij
99 double alpha_ij = m_alpha[ki] * m_alpha[kj];
100 m_aAlpha_binary(ki, kj) = m_aAlpha_binary(kj, ki) = a0*alpha_ij;
101}
102
103// ------------Molar Thermodynamic Properties -------------------------
104
106{
108 double T = temperature();
109 double mv = molarVolume();
110 double vpb = mv + (1 + Sqrt2) * m_b;
111 double vmb = mv + (1 - Sqrt2) * m_b;
113 double cpref = GasConstant * mean_X(m_cp0_R);
114 double dHdT_V = cpref + mv * m_dpdT - GasConstant
115 + 1.0 / (2.0 * Sqrt2 * m_b) * log(vpb / vmb) * T * d2aAlpha_dT2();
116 return dHdT_V - (mv + T * m_dpdT / m_dpdV) * m_dpdT;
117}
118
120{
122 double T = temperature();
124 return (cp_mole() + T * m_dpdT * m_dpdT / m_dpdV);
125}
126
128{
130 // Get a copy of the private variables stored in the State object
131 double T = temperature();
132 double mv = molarVolume();
133 double denom = mv * mv + 2 * mv * m_b - m_b * m_b;
134 return GasConstant * T / (mv - m_b) - m_aAlpha_mix / denom;
135}
136
138{
140 return 1.0 / m_tmpV[k];
141}
142
144{
145 double mv = molarVolume();
146 double vpb2 = mv + (1 + Sqrt2) * m_b;
147 double vmb2 = mv + (1 - Sqrt2) * m_b;
148 double vmb = mv - m_b;
149 double pres = pressure();
150
151 for (size_t k = 0; k < m_kk; k++) {
152 m_pp[k] = 0.0;
153 for (size_t i = 0; i < m_kk; i++) {
154 m_pp[k] += moleFractions_[i] * m_aAlpha_binary(k, i);
155 }
156 }
157 double num = 0;
158 double denom = 2 * Sqrt2 * m_b * m_b;
159 double denom2 = m_b * (mv * mv + 2 * mv * m_b - m_b * m_b);
160 double RT_ = RT();
161 for (size_t k = 0; k < m_kk; k++) {
162 num = 2 * m_b * m_pp[k] - m_aAlpha_mix * m_b_coeffs[k];
163 ac[k] = (-RT_ * log(pres * mv/ RT_) + RT_ * log(mv / vmb)
164 + RT_ * m_b_coeffs[k] / vmb
165 - (num /denom) * log(vpb2/vmb2)
166 - m_aAlpha_mix * m_b_coeffs[k] * mv/denom2
167 );
168 }
169 for (size_t k = 0; k < m_kk; k++) {
170 ac[k] = exp(ac[k]/ RT_);
171 }
172}
173
174// ---- Partial Molar Properties of the Solution -----------------
175
177{
178 getGibbs_ref(mu);
179 double RT_ = RT();
180 for (size_t k = 0; k < m_kk; k++) {
181 double xx = std::max(SmallNumber, moleFraction(k));
182 mu[k] += RT_ * (log(xx));
183 }
184
185 double mv = molarVolume();
186 double vmb = mv - m_b;
187 double vpb2 = mv + (1 + Sqrt2) * m_b;
188 double vmb2 = mv + (1 - Sqrt2) * m_b;
189
190 for (size_t k = 0; k < m_kk; k++) {
191 m_pp[k] = 0.0;
192 for (size_t i = 0; i < m_kk; i++) {
193 m_pp[k] += moleFractions_[i] * m_aAlpha_binary(k, i);
194 }
195 }
196 double pres = pressure();
197 double refP = refPressure();
198 double denom = 2 * Sqrt2 * m_b * m_b;
199 double denom2 = m_b * (mv * mv + 2 * mv * m_b - m_b * m_b);
200
201 for (size_t k = 0; k < m_kk; k++) {
202 double num = 2 * m_b * m_pp[k] - m_aAlpha_mix * m_b_coeffs[k];
203
204 mu[k] += RT_ * log(pres/refP) - RT_ * log(pres * mv / RT_)
205 + RT_ * log(mv / vmb) + RT_ * m_b_coeffs[k] / vmb
206 - (num /denom) * log(vpb2/vmb2)
207 - m_aAlpha_mix * m_b_coeffs[k] * mv/denom2;
208 }
209}
210
212{
213 // First we get the reference state contributions
214 getEnthalpy_RT_ref(hbar);
215 scale(hbar, hbar+m_kk, hbar, RT());
216 vector_fp tmp;
217 tmp.resize(m_kk,0.0);
218
219 // We calculate m_dpdni
220 double T = temperature();
221 double mv = molarVolume();
222 double vmb = mv - m_b;
223 double vpb2 = mv + (1 + Sqrt2) * m_b;
224 double vmb2 = mv + (1 - Sqrt2) * m_b;
225 double daAlphadT = daAlpha_dT();
226
227 for (size_t k = 0; k < m_kk; k++) {
228 m_pp[k] = 0.0;
229 tmp[k] = 0.0;
230 for (size_t i = 0; i < m_kk; i++) {
231 double grad_aAlpha = m_dalphadT[i]/m_alpha[i] + m_dalphadT[k]/m_alpha[k];
232 m_pp[k] += moleFractions_[i] * m_aAlpha_binary(k, i);
233 tmp[k] +=moleFractions_[i] * m_aAlpha_binary(k, i) * grad_aAlpha;
234 }
235 }
236
237 double denom = mv * mv + 2 * mv * m_b - m_b * m_b;
238 double denom2 = denom * denom;
239 double RT_ = RT();
240 for (size_t k = 0; k < m_kk; k++) {
241 m_dpdni[k] = RT_ / vmb + RT_ * m_b_coeffs[k] / (vmb * vmb)
242 - 2.0 * m_pp[k] / denom + 2 * vmb * m_aAlpha_mix * m_b_coeffs[k] / denom2;
243 }
244
245 double fac = T * daAlphadT - m_aAlpha_mix;
247 double fac2 = mv + T * m_dpdT / m_dpdV;
248 double fac3 = 2 * Sqrt2 * m_b * m_b;
249 double fac4 = 0;
250 for (size_t k = 0; k < m_kk; k++) {
251 fac4 = T*tmp[k] -2 * m_pp[k];
252 double hE_v = mv * m_dpdni[k] - RT_
253 - m_b_coeffs[k] / fac3 * log(vpb2 / vmb2) * fac
254 + (mv * m_b_coeffs[k]) / (m_b * denom) * fac
255 + 1 / (2 * Sqrt2 * m_b) * log(vpb2 / vmb2) * fac4;
256 hbar[k] = hbar[k] + hE_v;
257 hbar[k] -= fac2 * m_dpdni[k];
258 }
259}
260
262{
263 // Using the identity : (hk - T*sk) = gk
264 double T = temperature();
267 for (size_t k = 0; k < m_kk; k++) {
268 sbar[k] = (sbar[k] - m_tmpV[k])/T;
269 }
270}
271
273{
274 // u_i = h_i - p*v_i
275 double p = pressure();
278 for (size_t k = 0; k < m_kk; k++) {
279 ubar[k] = ubar[k] - p*m_tmpV[k];
280 }
281}
282
283void PengRobinson::getPartialMolarCp(double* cpbar) const
284{
285 throw NotImplementedError("PengRobinson::getPartialMolarCp");
286}
287
289{
290 for (size_t k = 0; k < m_kk; k++) {
291 m_pp[k] = 0.0;
292 for (size_t i = 0; i < m_kk; i++) {
293 m_pp[k] += moleFractions_[i] * m_aAlpha_binary(k, i);
294 }
295 }
296
297 double mv = molarVolume();
298 double vmb = mv - m_b;
299 double vpb = mv + m_b;
300 double fac = mv * mv + 2 * mv * m_b - m_b * m_b;
301 double fac2 = fac * fac;
302 double RT_ = RT();
303
304 for (size_t k = 0; k < m_kk; k++) {
305 double num = RT_ + RT_ * m_b/ vmb + RT_ * m_b_coeffs[k] / vmb
306 + RT_ * m_b * m_b_coeffs[k] /(vmb * vmb) - 2 * mv * m_pp[k] / fac
307 + 2 * mv * vmb * m_aAlpha_mix * m_b_coeffs[k] / fac2;
308 double denom = pressure() + RT_ * m_b / (vmb * vmb) + m_aAlpha_mix/fac
309 - 2 * mv* vpb * m_aAlpha_mix / fac2;
310 vbar[k] = num / denom;
311 }
312}
313
314double PengRobinson::speciesCritTemperature(double a, double b) const
315{
316 if (b <= 0.0) {
317 return 1000000.;
318 } else if (a <= 0.0) {
319 return 0.0;
320 } else {
321 return a * omega_b / (b * omega_a * GasConstant);
322 }
323}
324
325bool PengRobinson::addSpecies(shared_ptr<Species> spec)
326{
327 bool added = MixtureFugacityTP::addSpecies(spec);
328 if (added) {
329 m_a_coeffs.resize(m_kk, m_kk, 0.0);
330 m_b_coeffs.push_back(0.0);
331 m_aAlpha_binary.resize(m_kk, m_kk, 0.0);
332 m_kappa.push_back(0.0);
333 m_acentric.push_back(0.0);
334 m_alpha.push_back(0.0);
335 m_dalphadT.push_back(0.0);
336 m_d2alphadT2.push_back(0.0);
337 m_pp.push_back(0.0);
338 m_partialMolarVolumes.push_back(0.0);
339 m_dpdni.push_back(0.0);
340 m_coeffSource.push_back(CoeffSource::EoS);
341 }
342 return added;
343}
344
346{
347 // Contents of 'critical-properties.yaml', loaded later if needed
348 AnyMap critPropsDb;
349 std::unordered_map<std::string, AnyMap*> dbSpecies;
350
351 for (auto& item : m_species) {
352 auto& data = item.second->input;
353 size_t k = speciesIndex(item.first);
354 if (m_a_coeffs(k, k) != 0.0) {
355 continue;
356 }
357 bool foundCoeffs = false;
358 if (data.hasKey("equation-of-state") &&
359 data["equation-of-state"].hasMapWhere("model", "Peng-Robinson"))
360 {
361 // Read a and b coefficients and acentric factor w_ac from species input
362 // information, specified in a YAML input file.
363 auto eos = data["equation-of-state"].getMapWhere(
364 "model", "Peng-Robinson");
365 if (eos.hasKey("a") && eos.hasKey("b") && eos.hasKey("acentric-factor")) {
366 double a0 = eos.convert("a", "Pa*m^6/kmol^2");
367 double b = eos.convert("b", "m^3/kmol");
368 // unitless acentric factor:
369 double w = eos["acentric-factor"].asDouble();
370 setSpeciesCoeffs(item.first, a0, b, w);
371 foundCoeffs = true;
372 }
373
374 if (eos.hasKey("binary-a")) {
375 AnyMap& binary_a = eos["binary-a"].as<AnyMap>();
376 const UnitSystem& units = binary_a.units();
377 for (auto& item2 : binary_a) {
378 double a0 = units.convert(item2.second, "Pa*m^6/kmol^2");
379 setBinaryCoeffs(item.first, item2.first, a0);
380 }
381 }
382 if (foundCoeffs) {
383 m_coeffSource[k] = CoeffSource::EoS;
384 continue;
385 }
386 }
387
388 // Coefficients have not been populated from model-specific input
389 double Tc = NAN, Pc = NAN, omega_ac = NAN;
390 if (data.hasKey("critical-parameters")) {
391 // Use critical state information stored in the species entry to
392 // calculate a, b, and the acentric factor.
393 auto& critProps = data["critical-parameters"].as<AnyMap>();
394 Tc = critProps.convert("critical-temperature", "K");
395 Pc = critProps.convert("critical-pressure", "Pa");
396 omega_ac = critProps["acentric-factor"].asDouble();
397 m_coeffSource[k] = CoeffSource::CritProps;
398 } else {
399 // Search 'crit-properties.yaml' to find Tc and Pc. Load data if needed.
400 if (critPropsDb.empty()) {
401 critPropsDb = AnyMap::fromYamlFile("critical-properties.yaml");
402 dbSpecies = critPropsDb["species"].asMap("name");
403 }
404
405 // All names in critical-properties.yaml are upper case
406 auto ucName = boost::algorithm::to_upper_copy(item.first);
407 if (dbSpecies.count(ucName)) {
408 auto& spec = *dbSpecies.at(ucName);
409 auto& critProps = spec["critical-parameters"].as<AnyMap>();
410 Tc = critProps.convert("critical-temperature", "K");
411 Pc = critProps.convert("critical-pressure", "Pa");
412 omega_ac = critProps["acentric-factor"].asDouble();
413 m_coeffSource[k] = CoeffSource::Database;
414 }
415 }
416
417 // Check if critical properties were found in either location
418 if (!isnan(Tc)) {
419 double a = omega_a * std::pow(GasConstant * Tc, 2) / Pc;
420 double b = omega_b * GasConstant * Tc / Pc;
421 setSpeciesCoeffs(item.first, a, b, omega_ac);
422 } else {
423 throw InputFileError("PengRobinson::initThermo", data,
424 "No Peng-Robinson model parameters or critical properties found for "
425 "species '{}'", item.first);
426 }
427 }
428}
429
430void PengRobinson::getSpeciesParameters(const std::string& name,
431 AnyMap& speciesNode) const
432{
434 size_t k = speciesIndex(name);
436
437 // Pure species parameters
438 if (m_coeffSource[k] == CoeffSource::EoS) {
439 auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
440 "model", "Peng-Robinson", true);
441 eosNode["a"].setQuantity(m_a_coeffs(k, k), "Pa*m^6/kmol^2");
442 eosNode["b"].setQuantity(m_b_coeffs[k], "m^3/kmol");
443 eosNode["acentric-factor"] = m_acentric[k];
444 } else if (m_coeffSource[k] == CoeffSource::CritProps) {
445 auto& critProps = speciesNode["critical-parameters"];
446 double Tc = speciesCritTemperature(m_a_coeffs(k, k), m_b_coeffs[k]);
447 double Pc = omega_b * GasConstant * Tc / m_b_coeffs[k];
448 critProps["critical-temperature"].setQuantity(Tc, "K");
449 critProps["critical-pressure"].setQuantity(Pc, "Pa");
450 critProps["acentric-factor"] = m_acentric[k];
451 }
452 // Nothing to do in the case where the parameters are from the database
453
454 if (m_binaryParameters.count(name)) {
455 // Include binary parameters regardless of where the pure species parameters
456 // were found
457 auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
458 "model", "Peng-Robinson", true);
459 AnyMap bin_a;
460 for (const auto& item : m_binaryParameters.at(name)) {
461 bin_a[item.first].setQuantity(item.second, "Pa*m^6/kmol^2");
462 }
463 eosNode["binary-a"] = std::move(bin_a);
464 }
465}
466
468{
469 double molarV = molarVolume();
470 double hh = m_b / molarV;
471 double zz = z();
472 double alpha_1 = daAlpha_dT();
473 double vpb = molarV + (1.0 + Sqrt2) * m_b;
474 double vmb = molarV + (1.0 - Sqrt2) * m_b;
475 double fac = alpha_1 / (2.0 * Sqrt2 * m_b);
476 double sresid_mol_R = log(zz*(1.0 - hh)) + fac * log(vpb / vmb) / GasConstant;
477 return GasConstant * sresid_mol_R;
478}
479
481{
482 double molarV = molarVolume();
483 double zz = z();
484 double aAlpha_1 = daAlpha_dT();
485 double T = temperature();
486 double vpb = molarV + (1 + Sqrt2) * m_b;
487 double vmb = molarV + (1 - Sqrt2) * m_b;
488 double fac = 1 / (2.0 * Sqrt2 * m_b);
489 return GasConstant * T * (zz - 1.0)
490 + fac * log(vpb / vmb) * (T * aAlpha_1 - m_aAlpha_mix);
491}
492
493double PengRobinson::liquidVolEst(double T, double& presGuess) const
494{
495 double v = m_b * 1.1;
496 double atmp, btmp, aAlphatmp;
497 calculateAB(atmp, btmp, aAlphatmp);
498 double pres = std::max(psatEst(T), presGuess);
499 double Vroot[3];
500 bool foundLiq = false;
501 int m = 0;
502 while (m < 100 && !foundLiq) {
503 int nsol = solveCubic(T, pres, atmp, btmp, aAlphatmp, Vroot);
504 if (nsol == 1 || nsol == 2) {
505 double pc = critPressure();
506 if (pres > pc) {
507 foundLiq = true;
508 }
509 pres *= 1.04;
510 } else {
511 foundLiq = true;
512 }
513 }
514
515 if (foundLiq) {
516 v = Vroot[0];
517 presGuess = pres;
518 } else {
519 v = -1.0;
520 }
521 return v;
522}
523
524double PengRobinson::densityCalc(double T, double presPa, int phaseRequested,
525 double rhoGuess)
526{
527 // It's necessary to set the temperature so that m_aAlpha_mix is set correctly.
529 double tcrit = critTemperature();
530 double mmw = meanMolecularWeight();
531 if (rhoGuess == -1.0) {
532 if (phaseRequested >= FLUID_LIQUID_0) {
533 double lqvol = liquidVolEst(T, presPa);
534 rhoGuess = mmw / lqvol;
535 }
536 } else {
537 // Assume the Gas phase initial guess, if nothing is specified to the routine
538 rhoGuess = presPa * mmw / (GasConstant * T);
539 }
540
541 double volGuess = mmw / rhoGuess;
542 m_NSolns = solveCubic(T, presPa, m_a, m_b, m_aAlpha_mix, m_Vroot);
543
544 double molarVolLast = m_Vroot[0];
545 if (m_NSolns >= 2) {
546 if (phaseRequested >= FLUID_LIQUID_0) {
547 molarVolLast = m_Vroot[0];
548 } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
549 molarVolLast = m_Vroot[2];
550 } else {
551 if (volGuess > m_Vroot[1]) {
552 molarVolLast = m_Vroot[2];
553 } else {
554 molarVolLast = m_Vroot[0];
555 }
556 }
557 } else if (m_NSolns == 1) {
558 if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT
559 || phaseRequested == FLUID_UNDEFINED)
560 {
561 molarVolLast = m_Vroot[0];
562 } else {
563 return -2.0;
564 }
565 } else if (m_NSolns == -1) {
566 if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED
567 || phaseRequested == FLUID_SUPERCRIT)
568 {
569 molarVolLast = m_Vroot[0];
570 } else if (T > tcrit) {
571 molarVolLast = m_Vroot[0];
572 } else {
573 return -2.0;
574 }
575 } else {
576 molarVolLast = m_Vroot[0];
577 return -1.0;
578 }
579 return mmw / molarVolLast;
580}
581
583{
584 double Vroot[3];
585 double T = temperature();
586 int nsol = solveCubic(T, pressure(), m_a, m_b, m_aAlpha_mix, Vroot);
587 if (nsol != 3) {
588 return critDensity();
589 }
590
591 auto resid = [this, T](double v) {
592 double pp;
593 return dpdVCalc(T, v, pp);
594 };
595
596 boost::uintmax_t maxiter = 100;
597 std::pair<double, double> vv = bmt::toms748_solve(
598 resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
599
600 double mmw = meanMolecularWeight();
601 return mmw / (0.5 * (vv.first + vv.second));
602}
603
605{
606 double Vroot[3];
607 double T = temperature();
608 int nsol = solveCubic(T, pressure(), m_a, m_b, m_aAlpha_mix, Vroot);
609 if (nsol != 3) {
610 return critDensity();
611 }
612
613 auto resid = [this, T](double v) {
614 double pp;
615 return dpdVCalc(T, v, pp);
616 };
617
618 boost::uintmax_t maxiter = 100;
619 std::pair<double, double> vv = bmt::toms748_solve(
620 resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
621
622 double mmw = meanMolecularWeight();
623 return mmw / (0.5 * (vv.first + vv.second));
624}
625
626double PengRobinson::dpdVCalc(double T, double molarVol, double& presCalc) const
627{
628 double denom = molarVol * molarVol + 2 * molarVol * m_b - m_b * m_b;
629 double vpb = molarVol + m_b;
630 double vmb = molarVol - m_b;
631 return -GasConstant * T / (vmb * vmb) + 2 * m_aAlpha_mix * vpb / (denom*denom);
632}
633
635{
636 double T = temperature();
637 double mv = molarVolume();
638 double pres;
639
640 m_dpdV = dpdVCalc(T, mv, pres);
641 double vmb = mv - m_b;
642 double denom = mv * mv + 2 * mv * m_b - m_b * m_b;
643 m_dpdT = (GasConstant / vmb - daAlpha_dT() / denom);
644}
645
647{
648 double temp = temperature();
649
650 // Update individual alpha
651 for (size_t j = 0; j < m_kk; j++) {
652 double critTemp_j = speciesCritTemperature(m_a_coeffs(j,j), m_b_coeffs[j]);
653 double sqt_alpha = 1 + m_kappa[j] * (1 - sqrt(temp / critTemp_j));
654 m_alpha[j] = sqt_alpha*sqt_alpha;
655 }
656
657 // Update aAlpha_i, j
658 for (size_t i = 0; i < m_kk; i++) {
659 for (size_t j = 0; j < m_kk; j++) {
660 m_aAlpha_binary(i, j) = sqrt(m_alpha[i] * m_alpha[j]) * m_a_coeffs(i,j);
661 }
662 }
664}
665
666void PengRobinson::calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const
667{
668 bCalc = 0.0;
669 aCalc = 0.0;
670 aAlphaCalc = 0.0;
671 for (size_t i = 0; i < m_kk; i++) {
672 bCalc += moleFractions_[i] * m_b_coeffs[i];
673 for (size_t j = 0; j < m_kk; j++) {
674 aCalc += m_a_coeffs(i, j) * moleFractions_[i] * moleFractions_[j];
675 aAlphaCalc += m_aAlpha_binary(i, j) * moleFractions_[i] * moleFractions_[j];
676 }
677 }
678}
679
681{
682 double daAlphadT = 0.0, k, Tc, sqtTr, coeff1, coeff2;
683 for (size_t i = 0; i < m_kk; i++) {
684 // Calculate first derivative of alpha for individual species
685 Tc = speciesCritTemperature(m_a_coeffs(i,i), m_b_coeffs[i]);
686 sqtTr = sqrt(temperature() / Tc);
687 coeff1 = 1 / (Tc*sqtTr);
688 coeff2 = sqtTr - 1;
689 k = m_kappa[i];
690 m_dalphadT[i] = coeff1 * (k*k*coeff2 - k);
691 }
692 // Calculate mixture derivative
693 for (size_t i = 0; i < m_kk; i++) {
694 for (size_t j = 0; j < m_kk; j++) {
695 daAlphadT += moleFractions_[i] * moleFractions_[j] * 0.5
696 * m_aAlpha_binary(i, j)
697 * (m_dalphadT[i] / m_alpha[i] + m_dalphadT[j] / m_alpha[j]);
698 }
699 }
700 return daAlphadT;
701}
702
704{
705 for (size_t i = 0; i < m_kk; i++) {
706 double Tcrit_i = speciesCritTemperature(m_a_coeffs(i, i), m_b_coeffs[i]);
707 double sqt_Tr = sqrt(temperature() / Tcrit_i);
708 double coeff1 = 1 / (Tcrit_i*sqt_Tr);
709 double coeff2 = sqt_Tr - 1;
710 // Calculate first and second derivatives of alpha for individual species
711 double k = m_kappa[i];
712 m_dalphadT[i] = coeff1 * (k*k*coeff2 - k);
713 m_d2alphadT2[i] = (k*k + k) * coeff1 / (2*sqt_Tr*sqt_Tr*Tcrit_i);
714 }
715
716 // Calculate mixture derivative
717 double d2aAlphadT2 = 0.0;
718 for (size_t i = 0; i < m_kk; i++) {
719 double alphai = m_alpha[i];
720 for (size_t j = 0; j < m_kk; j++) {
721 double alphaj = m_alpha[j];
722 double alphaij = alphai * alphaj;
723 double term1 = m_d2alphadT2[i] / alphai + m_d2alphadT2[j] / alphaj;
724 double term2 = 2 * m_dalphadT[i] * m_dalphadT[j] / alphaij;
725 double term3 = m_dalphadT[i] / alphai + m_dalphadT[j] / alphaj;
726 d2aAlphadT2 += 0.5 * moleFractions_[i] * moleFractions_[j]
727 * m_aAlpha_binary(i, j)
728 * (term1 + term2 - 0.5 * term3 * term3);
729 }
730 }
731 return d2aAlphadT2;
732}
733
734void PengRobinson::calcCriticalConditions(double& pc, double& tc, double& vc) const
735{
736 if (m_b <= 0.0) {
737 tc = 1000000.;
738 pc = 1.0E13;
739 vc = omega_vc * GasConstant * tc / pc;
740 return;
741 }
742 if (m_a <= 0.0) {
743 tc = 0.0;
744 pc = 0.0;
745 vc = 2.0 * m_b;
746 return;
747 }
748 tc = m_a * omega_b / (m_b * omega_a * GasConstant);
749 pc = omega_b * GasConstant * tc / m_b;
750 vc = omega_vc * GasConstant * tc / pc;
751}
752
753int PengRobinson::solveCubic(double T, double pres, double a, double b, double aAlpha,
754 double Vroot[3]) const
755{
756 // Derive the coefficients of the cubic polynomial (in terms of molar volume v)
757 double bsqr = b * b;
758 double RT_p = GasConstant * T / pres;
759 double aAlpha_p = aAlpha / pres;
760 double an = 1.0;
761 double bn = (b - RT_p);
762 double cn = -(2 * RT_p * b - aAlpha_p + 3 * bsqr);
763 double dn = (bsqr * RT_p + bsqr * b - aAlpha_p * b);
764
765 double tc = a * omega_b / (b * omega_a * GasConstant);
766 double pc = omega_b * GasConstant * tc / b;
767 double vc = omega_vc * GasConstant * tc / pc;
768
769 return MixtureFugacityTP::solveCubic(T, pres, a, b, aAlpha, Vroot,
770 an, bn, cn, dn, tc, vc);
771}
772
773}
Declaration for class Cantera::Species.
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
A map of string keys to values whose type can vary at runtime.
Definition: AnyMap.h:399
double convert(const std::string &key, const std::string &units) const
Convert the item stored by the given key to the units specified in units.
Definition: AnyMap.cpp:1508
const UnitSystem & units() const
Return the default units that should be used to convert stored values.
Definition: AnyMap.h:598
bool empty() const
Return boolean indicating whether AnyMap is empty.
Definition: AnyMap.cpp:1401
static AnyMap fromYamlFile(const std::string &name, const std::string &parent_name="")
Create an AnyMap from a YAML file.
Definition: AnyMap.cpp:1743
void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition: Array.cpp:52
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:61
Error thrown for problems processing information contained in an AnyMap or AnyValue.
Definition: AnyMap.h:702
vector_fp m_cp0_R
Temporary storage for dimensionless reference state heat capacities.
virtual bool addSpecies(shared_ptr< Species > spec)
doublereal z() const
Calculate the value of z.
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of each species in their standard states at the current T and P of the solution...
int solveCubic(double T, double pres, double a, double b, double aAlpha, double Vroot[3], double an, double bn, double cn, double dn, double tc, double vc) const
Solve the cubic equation of state.
virtual doublereal psatEst(doublereal TKelvin) const
Estimate for the saturation pressure.
virtual double critDensity() const
Critical density (kg/m3).
virtual void setTemperature(const doublereal temp)
Set the temperature of the phase.
vector_fp moleFractions_
Storage for the current values of the mole fractions of the species.
virtual double critTemperature() const
Critical temperature (K).
virtual void getGibbs_ref(doublereal *g) const
Returns the vector of the Gibbs function of the reference state at the current temperature of the sol...
vector_fp m_tmpV
Temporary storage - length = m_kk.
virtual double critPressure() const
Critical pressure (Pa).
virtual void getEnthalpy_RT_ref(doublereal *hrt) const
An error indicating that an unimplemented function has been called.
Definition: ctexceptions.h:187
virtual double cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
virtual bool addSpecies(shared_ptr< Species > spec)
void setSpeciesCoeffs(const std::string &species, double a, double b, double w)
Set the pure fluid interaction parameters for a species.
virtual double standardConcentration(size_t k=0) const
Returns the standard concentration , which is used to normalize the generalized concentration.
virtual double densityCalc(double T, double pressure, int phase, double rhoguess)
Calculates the density given the temperature and the pressure and a guess at the density.
virtual double pressure() const
Return the thermodynamic pressure (Pa).
static const double omega_b
Omega constant: b0 (= omega_b) used in Peng-Robinson equation of state.
Definition: PengRobinson.h:321
virtual void getPartialMolarEnthalpies(double *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
vector_fp m_pp
Temporary storage - length = m_kk.
Definition: PengRobinson.h:280
virtual double densSpinodalLiquid() const
Return the value of the density at the liquid spinodal point (on the liquid side) for the current tem...
virtual void getActivityCoefficients(double *ac) const
Get the array of non-dimensional activity coefficients at the current solution temperature,...
virtual double dpdVCalc(double T, double molarVol, double &presCalc) const
Calculate the pressure and the pressure derivative given the temperature and the molar volume.
virtual double cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
std::vector< CoeffSource > m_coeffSource
For each species, specifies the source of the a, b, and omega coefficients.
Definition: PengRobinson.h:308
double m_dpdT
The derivative of the pressure with respect to the temperature.
Definition: PengRobinson.h:297
PengRobinson(const std::string &infile="", const std::string &id="")
Construct and initialize a PengRobinson object directly from an input file.
virtual void updateMixingExpressions()
Update the , , and parameters.
virtual void getChemPotentials(double *mu) const
Get the species chemical potentials. Units: J/kmol.
virtual double liquidVolEst(double T, double &pres) const
Estimate for the molar volume of the liquid.
virtual double densSpinodalGas() const
Return the value of the density at the gas spinodal point (on the gas side) for the current temperatu...
int solveCubic(double T, double pres, double a, double b, double aAlpha, double Vroot[3]) const
Prepare variables and call the function to solve the cubic equation of state.
std::map< std::string, std::map< std::string, double > > m_binaryParameters
Explicitly-specified binary interaction parameters, to enable serialization.
Definition: PengRobinson.h:273
virtual void initThermo()
Initialize the ThermoPhase object after all species have been set up.
double daAlpha_dT() const
Calculate temperature derivative .
void calculatePressureDerivatives() const
Calculate and at the current conditions.
virtual void getPartialMolarIntEnergies(double *ubar) const
Return an array of partial molar internal energies for the species in the mixture.
double m_aAlpha_mix
Value of in the equation of state.
Definition: PengRobinson.h:256
double m_a
Value of in the equation of state.
Definition: PengRobinson.h:250
static const double omega_a
Omega constant: a0 (= omega_a) used in Peng-Robinson equation of state.
Definition: PengRobinson.h:315
void setBinaryCoeffs(const std::string &species_i, const std::string &species_j, double a)
Set values for the interaction parameter between two species.
static const double omega_vc
Omega constant for the critical molar volume.
Definition: PengRobinson.h:324
virtual void getPartialMolarVolumes(double *vbar) const
Return an array of partial molar volumes for the species in the mixture.
void calculateAB(double &aCalc, double &bCalc, double &aAlpha) const
Calculate the , , and parameters given the temperature.
virtual double sresid() const
Calculate the deviation terms for the total entropy of the mixture from the ideal gas mixture.
double m_b
Value of in the equation of state.
Definition: PengRobinson.h:244
virtual void getSpeciesParameters(const std::string &name, AnyMap &speciesNode) const
Get phase-specific parameters of a Species object such that an identical one could be reconstructed a...
vector_fp m_dpdni
Vector of derivatives of pressure with respect to mole number.
Definition: PengRobinson.h:304
virtual double hresid() const
Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture.
double d2aAlpha_dT2() const
Calculate second derivative .
virtual void getPartialMolarEntropies(double *sbar) const
Returns an array of partial molar entropies of the species in the solution.
virtual void getPartialMolarCp(double *cpbar) const
Calculate species-specific molar specific heats.
double m_dpdV
The derivative of the pressure with respect to the volume.
Definition: PengRobinson.h:290
virtual double speciesCritTemperature(double a, double b) const
Calculate species-specific critical temperature.
vector_fp m_acentric
acentric factor for each species, length m_kk
Definition: PengRobinson.h:262
doublereal mean_X(const doublereal *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:717
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
Definition: Phase.cpp:211
std::string name() const
Return the name of the phase.
Definition: Phase.cpp:70
size_t m_kk
Number of species in the phase.
Definition: Phase.h:943
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:751
double moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition: Phase.cpp:548
doublereal temperature() const
Temperature (K).
Definition: Phase.h:654
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:187
double molarVolume() const
Molar volume (m^3/kmol).
Definition: Phase.cpp:682
shared_ptr< Species > species(const std::string &name) const
Return the Species object for the named species.
Definition: Phase.cpp:950
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:782
void initThermoFile(const std::string &inputFile, const std::string &id)
virtual doublereal refPressure() const
Returns the reference pressure in Pa.
Definition: ThermoPhase.h:150
virtual void getSpeciesParameters(const std::string &name, AnyMap &speciesNode) const
Get phase-specific parameters of a Species object such that an identical one could be reconstructed a...
Definition: ThermoPhase.h:1726
Unit conversion utility.
Definition: Units.h:161
double convert(double value, const std::string &src, const std::string &dest) const
Convert value from the units of src to the units of dest.
Definition: Units.cpp:558
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data.
const double Sqrt2
Sqrt(2)
Definition: ct_defs.h:55
Namespace for the Cantera kernel.
Definition: AnyMap.h:29
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:192
const double SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:153
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:100
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:184
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:113
Contains declarations for string manipulation functions within Cantera.