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
6 #include "cantera/thermo/PengRobinson.h"
7 #include "cantera/thermo/ThermoFactory.h"
8 #include "cantera/thermo/Species.h"
9 #include "cantera/base/stringUtils.h"
10 #include "cantera/base/ctml.h"
11
12 #include <boost/math/tools/roots.hpp>
13
14 #define _USE_MATH_DEFINES
15 #include <math.h>
16
17 using namespace std;
18 namespace bmt = boost::math::tools;
19
20 namespace Cantera
21 {
22
23 const double PengRobinson::omega_a = 4.5723552892138218E-01;
24 const double PengRobinson::omega_b = 7.77960739038885E-02;
25 const double PengRobinson::omega_vc = 3.07401308698703833E-01;
26
PengRobinson(const std::string & infile,const std::string & id_)27 PengRobinson::PengRobinson(const std::string& infile, const std::string& id_) :
28 m_b(0.0),
29 m_a(0.0),
30 m_aAlpha_mix(0.0),
31 m_NSolns(0),
32 m_dpdV(0.0),
33 m_dpdT(0.0)
34 {
35 fill_n(m_Vroot, 3, 0.0);
36 initThermoFile(infile, id_);
37 }
38
setSpeciesCoeffs(const std::string & species,double a,double b,double w)39 void PengRobinson::setSpeciesCoeffs(const std::string& species, double a, double b, double w)
40 {
41 size_t k = speciesIndex(species);
42 if (k == npos) {
43 throw CanteraError("PengRobinson::setSpeciesCoeffs",
44 "Unknown species '{}'.", species);
45 }
46
47 // Calculate value of kappa (independent of temperature)
48 // w is an acentric factor of species
49 if (w <= 0.491) {
50 m_kappa[k] = 0.37464 + 1.54226*w - 0.26992*w*w;
51 } else {
52 m_kappa[k] = 0.374642 + 1.487503*w - 0.164423*w*w + 0.016666*w*w*w;
53 }
54
55 //Calculate alpha (temperature dependent interaction parameter)
56 double critTemp = speciesCritTemperature(a, b); // critical temperature of individual species
57 double sqt_T_r = sqrt(temperature() / critTemp);
58 double sqt_alpha = 1 + m_kappa[k] * (1 - sqt_T_r);
59 m_alpha[k] = sqt_alpha*sqt_alpha;
60
61 m_a_coeffs(k,k) = a;
62 double aAlpha_k = a*m_alpha[k];
63 m_aAlpha_binary(k,k) = aAlpha_k;
64
65 // standard mixing rule for cross-species interaction term
66 for (size_t j = 0; j < m_kk; j++) {
67 if (k == j) {
68 continue;
69 }
70 double a0kj = sqrt(m_a_coeffs(j,j) * a);
71 double aAlpha_j = a*m_alpha[j];
72 double a_Alpha = sqrt(aAlpha_j*aAlpha_k);
73 if (m_a_coeffs(j, k) == 0) {
74 m_a_coeffs(j, k) = a0kj;
75 m_aAlpha_binary(j, k) = a_Alpha;
76 m_a_coeffs(k, j) = a0kj;
77 m_aAlpha_binary(k, j) = a_Alpha;
78 }
79 }
80 m_b_coeffs[k] = b;
81 }
82
setBinaryCoeffs(const std::string & species_i,const std::string & species_j,double a0)83 void PengRobinson::setBinaryCoeffs(const std::string& species_i,
84 const std::string& species_j, double a0)
85 {
86 size_t ki = speciesIndex(species_i);
87 if (ki == npos) {
88 throw CanteraError("PengRobinson::setBinaryCoeffs",
89 "Unknown species '{}'.", species_i);
90 }
91 size_t kj = speciesIndex(species_j);
92 if (kj == npos) {
93 throw CanteraError("PengRobinson::setBinaryCoeffs",
94 "Unknown species '{}'.", species_j);
95 }
96
97 m_a_coeffs(ki, kj) = m_a_coeffs(kj, ki) = 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
cp_mole() const105 double PengRobinson::cp_mole() const
106 {
107 _updateReferenceStateThermo();
108 double T = temperature();
109 double mv = molarVolume();
110 double vpb = mv + (1 + M_SQRT2)*m_b;
111 double vmb = mv + (1 - M_SQRT2)*m_b;
112 calculatePressureDerivatives();
113 double cpref = GasConstant * mean_X(m_cp0_R);
114 double dHdT_V = cpref + mv * m_dpdT - GasConstant
115 + 1.0 / (2.0 * M_SQRT2 * m_b) * log(vpb / vmb) * T * d2aAlpha_dT2();
116 return dHdT_V - (mv + T * m_dpdT / m_dpdV) * m_dpdT;
117 }
118
cv_mole() const119 double PengRobinson::cv_mole() const
120 {
121 _updateReferenceStateThermo();
122 double T = temperature();
123 calculatePressureDerivatives();
124 return (cp_mole() + T * m_dpdT * m_dpdT / m_dpdV);
125 }
126
pressure() const127 double PengRobinson::pressure() const
128 {
129 _updateReferenceStateThermo();
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 double pp = GasConstant * T / (mv - m_b) - m_aAlpha_mix / denom;
135 return pp;
136 }
137
standardConcentration(size_t k) const138 double PengRobinson::standardConcentration(size_t k) const
139 {
140 getStandardVolumes(m_tmpV.data());
141 return 1.0 / m_tmpV[k];
142 }
143
getActivityCoefficients(double * ac) const144 void PengRobinson::getActivityCoefficients(double* ac) const
145 {
146 double mv = molarVolume();
147 double vpb2 = mv + (1 + M_SQRT2)*m_b;
148 double vmb2 = mv + (1 - M_SQRT2)*m_b;
149 double vmb = mv - m_b;
150 double pres = pressure();
151
152 for (size_t k = 0; k < m_kk; k++) {
153 m_pp[k] = 0.0;
154 for (size_t i = 0; i < m_kk; i++) {
155 m_pp[k] += moleFractions_[i] * m_aAlpha_binary(k, i);
156 }
157 }
158 double num = 0;
159 double denom = 2 * M_SQRT2 * m_b * m_b;
160 double denom2 = m_b * (mv * mv + 2 * mv * m_b - m_b * m_b);
161 double RTkelvin = RT();
162 for (size_t k = 0; k < m_kk; k++) {
163 num = 2 * m_b * m_pp[k] - m_aAlpha_mix * m_b_coeffs[k];
164 ac[k] = (-RTkelvin * log(pres * mv/ RTkelvin) + RTkelvin * log(mv / vmb)
165 + RTkelvin * m_b_coeffs[k] / vmb
166 - (num /denom) * log(vpb2/vmb2)
167 - m_aAlpha_mix * m_b_coeffs[k] * mv/denom2
168 );
169 }
170 for (size_t k = 0; k < m_kk; k++) {
171 ac[k] = exp(ac[k]/ RTkelvin);
172 }
173 }
174
175 // ---- Partial Molar Properties of the Solution -----------------
176
getChemPotentials(double * mu) const177 void PengRobinson::getChemPotentials(double* mu) const
178 {
179 getGibbs_ref(mu);
180 double RTkelvin = RT();
181 for (size_t k = 0; k < m_kk; k++) {
182 double xx = std::max(SmallNumber, moleFraction(k));
183 mu[k] += RTkelvin * (log(xx));
184 }
185
186 double mv = molarVolume();
187 double vmb = mv - m_b;
188 double vpb2 = mv + (1 + M_SQRT2)*m_b;
189 double vmb2 = mv + (1 - M_SQRT2)*m_b;
190
191 for (size_t k = 0; k < m_kk; k++) {
192 m_pp[k] = 0.0;
193 for (size_t i = 0; i < m_kk; i++) {
194 m_pp[k] += moleFractions_[i] * m_aAlpha_binary(k, i);
195 }
196 }
197 double pres = pressure();
198 double refP = refPressure();
199 double denom = 2 * M_SQRT2 * m_b * m_b;
200 double denom2 = m_b * (mv * mv + 2 * mv * m_b - m_b * m_b);
201
202 for (size_t k = 0; k < m_kk; k++) {
203 double num = 2 * m_b * m_pp[k] - m_aAlpha_mix * m_b_coeffs[k];
204
205 mu[k] += (RTkelvin * log(pres/refP) - RTkelvin * log(pres * mv / RTkelvin)
206 + RTkelvin * log(mv / vmb)
207 + RTkelvin * m_b_coeffs[k] / vmb
208 - (num /denom) * log(vpb2/vmb2)
209 - m_aAlpha_mix * m_b_coeffs[k] * mv/denom2
210 );
211 }
212 }
213
getPartialMolarEnthalpies(double * hbar) const214 void PengRobinson::getPartialMolarEnthalpies(double* hbar) const
215 {
216 // First we get the reference state contributions
217 getEnthalpy_RT_ref(hbar);
218 scale(hbar, hbar+m_kk, hbar, RT());
219 vector_fp tmp;
220 tmp.resize(m_kk,0.0);
221
222 // We calculate m_dpdni
223 double T = temperature();
224 double mv = molarVolume();
225 double vmb = mv - m_b;
226 double vpb2 = mv + (1 + M_SQRT2)*m_b;
227 double vmb2 = mv + (1 - M_SQRT2)*m_b;
228 double daAlphadT = daAlpha_dT();
229
230 for (size_t k = 0; k < m_kk; k++) {
231 m_pp[k] = 0.0;
232 tmp[k] = 0.0;
233 for (size_t i = 0; i < m_kk; i++) {
234 double grad_aAlpha = m_dalphadT[i]/m_alpha[i] + m_dalphadT[k]/m_alpha[k];
235 m_pp[k] += moleFractions_[i] * m_aAlpha_binary(k, i);
236 tmp[k] +=moleFractions_[i] * m_aAlpha_binary(k, i) * grad_aAlpha;
237 }
238 }
239
240 double denom = mv * mv + 2 * mv * m_b - m_b * m_b;
241 double denom2 = denom * denom;
242 double RTkelvin = RT();
243 for (size_t k = 0; k < m_kk; k++) {
244 m_dpdni[k] = RTkelvin / vmb + RTkelvin * m_b_coeffs[k] / (vmb * vmb) - 2.0 * m_pp[k] / denom
245 + 2 * vmb * m_aAlpha_mix * m_b_coeffs[k] / denom2;
246 }
247
248 double fac = T * daAlphadT - m_aAlpha_mix;
249 calculatePressureDerivatives();
250 double fac2 = mv + T * m_dpdT / m_dpdV;
251 double fac3 = 2 * M_SQRT2 * m_b * m_b;
252 double fac4 = 0;
253 for (size_t k = 0; k < m_kk; k++) {
254 fac4 = T*tmp[k] -2 * m_pp[k];
255 double hE_v = mv * m_dpdni[k] - RTkelvin - m_b_coeffs[k] / fac3 * log(vpb2 / vmb2) * fac
256 + (mv * m_b_coeffs[k]) /(m_b * denom) * fac
257 + 1/(2 * M_SQRT2 * m_b) * log(vpb2 / vmb2) * fac4;
258 hbar[k] = hbar[k] + hE_v;
259 hbar[k] -= fac2 * m_dpdni[k];
260 }
261 }
262
getPartialMolarEntropies(double * sbar) const263 void PengRobinson::getPartialMolarEntropies(double* sbar) const
264 {
265 // Using the identity : (hk - T*sk) = gk
266 double T = temperature();
267 getPartialMolarEnthalpies(sbar);
268 getChemPotentials(m_tmpV.data());
269 for (size_t k = 0; k < m_kk; k++) {
270 sbar[k] = (sbar[k] - m_tmpV[k])/T;
271 }
272 }
273
getPartialMolarIntEnergies(double * ubar) const274 void PengRobinson::getPartialMolarIntEnergies(double* ubar) const
275 {
276 // u_i = h_i - p*v_i
277 double p = pressure();
278 getPartialMolarEnthalpies(ubar);
279 getPartialMolarVolumes(m_tmpV.data());
280 for (size_t k = 0; k < m_kk; k++) {
281 ubar[k] = ubar[k] - p*m_tmpV[k];
282 }
283 }
284
getPartialMolarCp(double * cpbar) const285 void PengRobinson::getPartialMolarCp(double* cpbar) const
286 {
287 throw NotImplementedError("PengRobinson::getPartialMolarCp");
288 }
289
getPartialMolarVolumes(double * vbar) const290 void PengRobinson::getPartialMolarVolumes(double* vbar) const
291 {
292 for (size_t k = 0; k < m_kk; k++) {
293 m_pp[k] = 0.0;
294 for (size_t i = 0; i < m_kk; i++) {
295 m_pp[k] += moleFractions_[i] * m_aAlpha_binary(k, i);
296 }
297 }
298
299 double mv = molarVolume();
300 double vmb = mv - m_b;
301 double vpb = mv + m_b;
302 double fac = mv * mv + 2 * mv * m_b - m_b * m_b;
303 double fac2 = fac * fac;
304 double RTkelvin = RT();
305
306 for (size_t k = 0; k < m_kk; k++) {
307 double num = (RTkelvin + RTkelvin * m_b/ vmb + RTkelvin * m_b_coeffs[k] / vmb
308 + RTkelvin * m_b * m_b_coeffs[k] /(vmb * vmb)
309 - 2 * mv * m_pp[k] / fac
310 + 2 * mv * vmb * m_aAlpha_mix * m_b_coeffs[k] / fac2
311 );
312 double denom = (pressure() + RTkelvin * m_b / (vmb * vmb)
313 + m_aAlpha_mix/fac
314 - 2 * mv* vpb * m_aAlpha_mix / fac2
315 );
316 vbar[k] = num / denom;
317 }
318 }
319
speciesCritTemperature(double a,double b) const320 double PengRobinson::speciesCritTemperature(double a, double b) const
321 {
322 if (b <= 0.0) {
323 return 1000000.;
324 } else if (a <= 0.0) {
325 return 0.0;
326 } else {
327 return a * omega_b / (b * omega_a * GasConstant);
328 }
329 }
330
addSpecies(shared_ptr<Species> spec)331 bool PengRobinson::addSpecies(shared_ptr<Species> spec)
332 {
333 bool added = MixtureFugacityTP::addSpecies(spec);
334 if (added) {
335 m_a_coeffs.resize(m_kk, m_kk, 0.0);
336 m_b_coeffs.push_back(0.0);
337 m_aAlpha_binary.resize(m_kk, m_kk, 0.0);
338 m_kappa.push_back(0.0);
339
340 m_alpha.push_back(0.0);
341 m_dalphadT.push_back(0.0);
342 m_d2alphadT2.push_back(0.0);
343
344 m_pp.push_back(0.0);
345 m_partialMolarVolumes.push_back(0.0);
346 m_dpdni.push_back(0.0);
347 }
348 return added;
349 }
350
getCoeff(const std::string & iName)351 vector<double> PengRobinson::getCoeff(const std::string& iName)
352 {
353 vector_fp spCoeff{ NAN, NAN, NAN };
354
355 // Get number of species in the database
356 // open xml file critProperties.xml
357 XML_Node* doc = get_XML_File("critProperties.xml");
358 size_t nDatabase = doc->nChildren();
359
360 // Loop through all species in the database and attempt to match supplied
361 // species to each. If present, calculate pureFluidParameters a_k and b_k
362 // based on crit properties T_c and P_c:
363 for (size_t isp = 0; isp < nDatabase; isp++) {
364 XML_Node& acNodeDoc = doc->child(isp);
365 std::string iNameLower = toLowerCopy(iName);
366 std::string dbName = toLowerCopy(acNodeDoc.attrib("name"));
367
368 // Attempt to match provided species iName to current database species
369 // dbName:
370 if (iNameLower == dbName) {
371 // Read from database and calculate a and b coefficients
372 double vParams;
373 double T_crit = 0.0, P_crit = 0.0, w_ac = 0.0;
374
375 if (acNodeDoc.hasChild("Tc")) {
376 vParams = 0.0;
377 XML_Node& xmlChildCoeff = acNodeDoc.child("Tc");
378 if (xmlChildCoeff.hasAttrib("value")) {
379 std::string critTemp = xmlChildCoeff.attrib("value");
380 vParams = strSItoDbl(critTemp);
381 }
382 if (vParams <= 0.0) { //Assuming that Pc and Tc are non zero.
383 throw CanteraError("PengRobinson::getCoeff",
384 "Critical Temperature must be positive");
385 }
386 T_crit = vParams;
387 }
388 if (acNodeDoc.hasChild("Pc")) {
389 vParams = 0.0;
390 XML_Node& xmlChildCoeff = acNodeDoc.child("Pc");
391 if (xmlChildCoeff.hasAttrib("value")) {
392 std::string critPressure = xmlChildCoeff.attrib("value");
393 vParams = strSItoDbl(critPressure);
394 }
395 if (vParams <= 0.0) { //Assuming that Pc and Tc are non zero.
396 throw CanteraError("PengRobinson::getCoeff",
397 "Critical Pressure must be positive");
398 }
399 P_crit = vParams;
400 }
401 if (acNodeDoc.hasChild("omega")) {
402 vParams = 0.0;
403 XML_Node& xmlChildCoeff = acNodeDoc.child("omega");
404 if (xmlChildCoeff.hasChild("value")) {
405 std::string acentric_factor = xmlChildCoeff.attrib("value");
406 vParams = strSItoDbl(acentric_factor);
407 }
408 w_ac = vParams;
409 }
410
411 spCoeff[0] = omega_a * (GasConstant * GasConstant) * (T_crit * T_crit) / P_crit; //coeff a
412 spCoeff[1] = omega_b * GasConstant * T_crit / P_crit; // coeff b
413 spCoeff[2] = w_ac; // acentric factor
414 break;
415 }
416 }
417 // If the species is not present in the database, throw an error
418 if(isnan(spCoeff[0]))
419 {
420 throw CanteraError("PengRobinson::getCoeff",
421 "Species '{}' is not present in the database", iName);
422 }
423 return spCoeff;
424 }
425
initThermo()426 void PengRobinson::initThermo()
427 {
428 for (auto& item : m_species) {
429 // Read a and b coefficients and acentric factor w_ac from species input
430 // information, specified in a YAML input file.
431 if (item.second->input.hasKey("equation-of-state")) {
432 auto eos = item.second->input["equation-of-state"].getMapWhere(
433 "model", "Peng-Robinson");
434 double a0 = 0;
435 if (eos["a"].isScalar()) {
436 a0 = eos.convert("a", "Pa*m^6/kmol^2");
437 }
438 double b = eos.convert("b", "m^3/kmol");
439 // unitless acentric factor:
440 double w = eos["acentric-factor"].asDouble();
441
442 setSpeciesCoeffs(item.first, a0, b, w);
443 if (eos.hasKey("binary-a")) {
444 AnyMap& binary_a = eos["binary-a"].as<AnyMap>();
445 const UnitSystem& units = binary_a.units();
446 for (auto& item2 : binary_a) {
447 double a0 = 0;
448 if (item2.second.isScalar()) {
449 a0 = units.convert(item2.second, "Pa*m^6/kmol^2");
450 }
451 setBinaryCoeffs(item.first, item2.first, a0);
452 }
453 }
454 } else {
455 // Check if a and b are already populated for this species (only the
456 // diagonal elements of a). If not, then search 'critProperties.xml'
457 // to find critical temperature and pressure to calculate a and b.
458 size_t k = speciesIndex(item.first);
459 if (m_a_coeffs(k, k) == 0.0) {
460 vector<double> coeffs = getCoeff(item.first);
461
462 // Check if species was found in the database of critical
463 // properties, and assign the results
464 if (!isnan(coeffs[0])) {
465 setSpeciesCoeffs(item.first, coeffs[0], coeffs[1], coeffs[2]);
466 }
467 }
468 }
469 }
470 }
471
sresid() const472 double PengRobinson::sresid() const
473 {
474 double molarV = molarVolume();
475 double hh = m_b / molarV;
476 double zz = z();
477 double alpha_1 = daAlpha_dT();
478 double vpb = molarV + (1.0 + M_SQRT2) * m_b;
479 double vmb = molarV + (1.0 - M_SQRT2) * m_b;
480 double fac = alpha_1 / (2.0 * M_SQRT2 * m_b);
481 double sresid_mol_R = log(zz*(1.0 - hh)) + fac * log(vpb / vmb) / GasConstant;
482 return GasConstant * sresid_mol_R;
483 }
484
hresid() const485 double PengRobinson::hresid() const
486 {
487 double molarV = molarVolume();
488 double zz = z();
489 double aAlpha_1 = daAlpha_dT();
490 double T = temperature();
491 double vpb = molarV + (1 + M_SQRT2) * m_b;
492 double vmb = molarV + (1 - M_SQRT2) * m_b;
493 double fac = 1 / (2.0 * M_SQRT2 * m_b);
494 return GasConstant * T * (zz - 1.0) + fac * log(vpb / vmb) * (T * aAlpha_1 - m_aAlpha_mix);
495 }
496
liquidVolEst(double T,double & presGuess) const497 double PengRobinson::liquidVolEst(double T, double& presGuess) const
498 {
499 double v = m_b * 1.1;
500 double atmp;
501 double btmp;
502 double aAlphatmp;
503 calculateAB(atmp, btmp, aAlphatmp);
504 double pres = std::max(psatEst(T), presGuess);
505 double Vroot[3];
506 bool foundLiq = false;
507 int m = 0;
508 while (m < 100 && !foundLiq) {
509 int nsol = solveCubic(T, pres, atmp, btmp, aAlphatmp, Vroot);
510 if (nsol == 1 || nsol == 2) {
511 double pc = critPressure();
512 if (pres > pc) {
513 foundLiq = true;
514 }
515 pres *= 1.04;
516 } else {
517 foundLiq = true;
518 }
519 }
520
521 if (foundLiq) {
522 v = Vroot[0];
523 presGuess = pres;
524 } else {
525 v = -1.0;
526 }
527 return v;
528 }
529
densityCalc(double T,double presPa,int phaseRequested,double rhoGuess)530 double PengRobinson::densityCalc(double T, double presPa, int phaseRequested, double rhoGuess)
531 {
532 // It's necessary to set the temperature so that m_aAlpha_mix is set correctly.
533 setTemperature(T);
534 double tcrit = critTemperature();
535 double mmw = meanMolecularWeight();
536 if (rhoGuess == -1.0) {
537 if (phaseRequested >= FLUID_LIQUID_0) {
538 double lqvol = liquidVolEst(T, presPa);
539 rhoGuess = mmw / lqvol;
540 }
541 } else {
542 // Assume the Gas phase initial guess, if nothing is specified to the routine
543 rhoGuess = presPa * mmw / (GasConstant * T);
544 }
545
546 double volGuess = mmw / rhoGuess;
547 m_NSolns = solveCubic(T, presPa, m_a, m_b, m_aAlpha_mix, m_Vroot);
548
549 double molarVolLast = m_Vroot[0];
550 if (m_NSolns >= 2) {
551 if (phaseRequested >= FLUID_LIQUID_0) {
552 molarVolLast = m_Vroot[0];
553 } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
554 molarVolLast = m_Vroot[2];
555 } else {
556 if (volGuess > m_Vroot[1]) {
557 molarVolLast = m_Vroot[2];
558 } else {
559 molarVolLast = m_Vroot[0];
560 }
561 }
562 } else if (m_NSolns == 1) {
563 if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) {
564 molarVolLast = m_Vroot[0];
565 } else {
566 return -2.0;
567 }
568 } else if (m_NSolns == -1) {
569 if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) {
570 molarVolLast = m_Vroot[0];
571 } else if (T > tcrit) {
572 molarVolLast = m_Vroot[0];
573 } else {
574 return -2.0;
575 }
576 } else {
577 molarVolLast = m_Vroot[0];
578 return -1.0;
579 }
580 return mmw / molarVolLast;
581 }
582
densSpinodalLiquid() const583 double PengRobinson::densSpinodalLiquid() const
584 {
585 double Vroot[3];
586 double T = temperature();
587 int nsol = solveCubic(T, pressure(), m_a, m_b, m_aAlpha_mix, Vroot);
588 if (nsol != 3) {
589 return critDensity();
590 }
591
592 auto resid = [this, T](double v) {
593 double pp;
594 return dpdVCalc(T, v, pp);
595 };
596
597 boost::uintmax_t maxiter = 100;
598 std::pair<double, double> vv = bmt::toms748_solve(
599 resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
600
601 double mmw = meanMolecularWeight();
602 return mmw / (0.5 * (vv.first + vv.second));
603 }
604
densSpinodalGas() const605 double PengRobinson::densSpinodalGas() const
606 {
607 double Vroot[3];
608 double T = temperature();
609 int nsol = solveCubic(T, pressure(), m_a, m_b, m_aAlpha_mix, Vroot);
610 if (nsol != 3) {
611 return critDensity();
612 }
613
614 auto resid = [this, T](double v) {
615 double pp;
616 return dpdVCalc(T, v, pp);
617 };
618
619 boost::uintmax_t maxiter = 100;
620 std::pair<double, double> vv = bmt::toms748_solve(
621 resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
622
623 double mmw = meanMolecularWeight();
624 return mmw / (0.5 * (vv.first + vv.second));
625 }
626
dpdVCalc(double T,double molarVol,double & presCalc) const627 double PengRobinson::dpdVCalc(double T, double molarVol, double& presCalc) const
628 {
629 double denom = molarVol * molarVol + 2 * molarVol * m_b - m_b * m_b;
630 double vpb = molarVol + m_b;
631 double vmb = molarVol - m_b;
632 double dpdv = -GasConstant * T / (vmb * vmb) + 2 * m_aAlpha_mix * vpb / (denom*denom);
633 return dpdv;
634 }
635
calculatePressureDerivatives() const636 void PengRobinson::calculatePressureDerivatives() const
637 {
638 double T = temperature();
639 double mv = molarVolume();
640 double pres;
641
642 m_dpdV = dpdVCalc(T, mv, pres);
643 double vmb = mv - m_b;
644 double denom = mv * mv + 2 * mv * m_b - m_b * m_b;
645 m_dpdT = (GasConstant / vmb - daAlpha_dT() / denom);
646 }
647
updateMixingExpressions()648 void PengRobinson::updateMixingExpressions()
649 {
650 double temp = temperature();
651
652 // Update individual alpha
653 for (size_t j = 0; j < m_kk; j++) {
654 double critTemp_j = speciesCritTemperature(m_a_coeffs(j,j), m_b_coeffs[j]);
655 double sqt_alpha = 1 + m_kappa[j] * (1 - sqrt(temp / critTemp_j));
656 m_alpha[j] = sqt_alpha*sqt_alpha;
657 }
658
659 //Update aAlpha_i, j
660 for (size_t i = 0; i < m_kk; i++) {
661 for (size_t j = 0; j < m_kk; j++) {
662 m_aAlpha_binary(i, j) = sqrt(m_alpha[i] * m_alpha[j]) * m_a_coeffs(i,j);
663 }
664 }
665 calculateAB(m_a,m_b,m_aAlpha_mix);
666 }
667
calculateAB(double & aCalc,double & bCalc,double & aAlphaCalc) const668 void PengRobinson::calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const
669 {
670 bCalc = 0.0;
671 aCalc = 0.0;
672 aAlphaCalc = 0.0;
673 for (size_t i = 0; i < m_kk; i++) {
674 bCalc += moleFractions_[i] * m_b_coeffs[i];
675 for (size_t j = 0; j < m_kk; j++) {
676 aCalc += m_a_coeffs(i, j) * moleFractions_[i] * moleFractions_[j];
677 aAlphaCalc += m_aAlpha_binary(i, j) * moleFractions_[i] * moleFractions_[j];
678 }
679 }
680 }
681
daAlpha_dT() const682 double PengRobinson::daAlpha_dT() const
683 {
684 double daAlphadT = 0.0, k, Tc, sqtTr, coeff1, coeff2;
685 for (size_t i = 0; i < m_kk; i++) {
686 // Calculate first derivative of alpha for individual species
687 Tc = speciesCritTemperature(m_a_coeffs(i,i), m_b_coeffs[i]);
688 sqtTr = sqrt(temperature() / Tc); //we need species critical temperature
689 coeff1 = 1 / (Tc*sqtTr);
690 coeff2 = sqtTr - 1;
691 k = m_kappa[i];
692 m_dalphadT[i] = coeff1 * (k*k*coeff2 - k);
693 }
694 //Calculate mixture derivative
695 for (size_t i = 0; i < m_kk; i++) {
696 for (size_t j = 0; j < m_kk; j++) {
697 daAlphadT += moleFractions_[i] * moleFractions_[j] * 0.5 * m_aAlpha_binary(i, j)
698 * (m_dalphadT[i] / m_alpha[i] + m_dalphadT[j] / m_alpha[j]);
699 }
700 }
701 return daAlphadT;
702 }
703
d2aAlpha_dT2() const704 double PengRobinson::d2aAlpha_dT2() const
705 {
706 for (size_t i = 0; i < m_kk; i++) {
707 double Tcrit_i = speciesCritTemperature(m_a_coeffs(i, i), m_b_coeffs[i]);
708 double sqt_Tr = sqrt(temperature() / Tcrit_i); //we need species critical temperature
709 double coeff1 = 1 / (Tcrit_i*sqt_Tr);
710 double coeff2 = sqt_Tr - 1;
711 // Calculate first and second derivatives of alpha for individual species
712 double k = m_kappa[i];
713 m_dalphadT[i] = coeff1 * (k*k*coeff2 - k);
714 m_d2alphadT2[i] = (k*k + k) * coeff1 / (2*sqt_Tr*sqt_Tr*Tcrit_i);
715 }
716
717 //Calculate mixture derivative
718 double d2aAlphadT2 = 0.0;
719 for (size_t i = 0; i < m_kk; i++) {
720 double alphai = m_alpha[i];
721 for (size_t j = 0; j < m_kk; j++) {
722 double alphaj = m_alpha[j];
723 double alphaij = alphai * alphaj;
724 double term1 = m_d2alphadT2[i] / alphai + m_d2alphadT2[j] / alphaj;
725 double term2 = 2 * m_dalphadT[i] * m_dalphadT[j] / alphaij;
726 double term3 = m_dalphadT[i] / alphai + m_dalphadT[j] / alphaj;
727 d2aAlphadT2 += 0.5 * moleFractions_[i] * moleFractions_[j] * m_aAlpha_binary(i, j)
728 * (term1 + term2 - 0.5 * term3 * term3);
729 }
730 }
731 return d2aAlphadT2;
732 }
733
calcCriticalConditions(double & pc,double & tc,double & vc) const734 void 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
solveCubic(double T,double pres,double a,double b,double aAlpha,double Vroot[3]) const753 int PengRobinson::solveCubic(double T, double pres, double a, double b, double aAlpha, double Vroot[3]) const
754 {
755 // Derive the coefficients of the cubic polynomial (in terms of molar volume v) to solve.
756 double bsqr = b * b;
757 double RT_p = GasConstant * T / pres;
758 double aAlpha_p = aAlpha / pres;
759 double an = 1.0;
760 double bn = (b - RT_p);
761 double cn = -(2 * RT_p * b - aAlpha_p + 3 * bsqr);
762 double dn = (bsqr * RT_p + bsqr * b - aAlpha_p * b);
763
764 double tc = a * omega_b / (b * omega_a * GasConstant);
765 double pc = omega_b * GasConstant * tc / b;
766 double vc = omega_vc * GasConstant * tc / pc;
767
768 int nSolnValues = MixtureFugacityTP::solveCubic(T, pres, a, b, aAlpha, Vroot, an, bn, cn, dn, tc, vc);
769
770 return nSolnValues;
771 }
772
773 }
774