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