Cantera  3.1.0b1
Loading...
Searching...
No Matches
RedlichKwongMFTP.cpp
Go to the documentation of this file.
1//! @file RedlichKwongMFTP.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 RedlichKwongMFTP::omega_a = 4.27480233540E-01;
21const double RedlichKwongMFTP::omega_b = 8.66403499650E-02;
22const double RedlichKwongMFTP::omega_vc = 3.33333333333333E-01;
23
24RedlichKwongMFTP::RedlichKwongMFTP(const string& infile, const string& id_)
25{
26 initThermoFile(infile, id_);
27}
28
29void RedlichKwongMFTP::setSpeciesCoeffs(const string& species,
30 double a0, double a1, double b)
31{
32 size_t k = speciesIndex(species);
33 if (k == npos) {
34 throw CanteraError("RedlichKwongMFTP::setSpeciesCoeffs",
35 "Unknown species '{}'.", species);
36 }
37
38 if (a1 != 0.0) {
39 m_formTempParam = 1; // expression is temperature-dependent
40 }
41
42 size_t counter = k + m_kk * k;
43 a_coeff_vec(0, counter) = a0;
44 a_coeff_vec(1, counter) = a1;
45
46 // standard mixing rule for cross-species interaction term
47 for (size_t j = 0; j < m_kk; j++) {
48 if (k == j) {
49 continue;
50 }
51
52 // a_coeff_vec(0) is initialized to NaN to mark uninitialized species
53 if (isnan(a_coeff_vec(0, j + m_kk * j))) {
54 // The diagonal element of the jth species has not yet been defined.
55 continue;
56 } else if (isnan(a_coeff_vec(0, j + m_kk * k))) {
57 // Only use the mixing rules if the off-diagonal element has not already been defined by a
58 // user-specified crossFluidParameters entry:
59 double a0kj = sqrt(a_coeff_vec(0, j + m_kk * j) * a0);
60 double a1kj = sqrt(a_coeff_vec(1, j + m_kk * j) * a1);
61 a_coeff_vec(0, j + m_kk * k) = a0kj;
62 a_coeff_vec(1, j + m_kk * k) = a1kj;
63 a_coeff_vec(0, k + m_kk * j) = a0kj;
64 a_coeff_vec(1, k + m_kk * j) = a1kj;
65 }
66 }
67 a_coeff_vec.getRow(0, a_vec_Curr_.data());
68 b_vec_Curr_[k] = b;
69}
70
71void RedlichKwongMFTP::setBinaryCoeffs(const string& species_i, const string& species_j,
72 double a0, double a1)
73{
74 size_t ki = speciesIndex(species_i);
75 if (ki == npos) {
76 throw CanteraError("RedlichKwongMFTP::setBinaryCoeffs",
77 "Unknown species '{}'.", species_i);
78 }
79 size_t kj = speciesIndex(species_j);
80 if (kj == npos) {
81 throw CanteraError("RedlichKwongMFTP::setBinaryCoeffs",
82 "Unknown species '{}'.", species_j);
83 }
84
85 if (a1 != 0.0) {
86 m_formTempParam = 1; // expression is temperature-dependent
87 }
88
89 m_binaryParameters[species_i][species_j] = {a0, a1};
90 m_binaryParameters[species_j][species_i] = {a0, a1};
91 size_t counter1 = ki + m_kk * kj;
92 size_t counter2 = kj + m_kk * ki;
93 a_coeff_vec(0, counter1) = a_coeff_vec(0, counter2) = a0;
94 a_coeff_vec(1, counter1) = a_coeff_vec(1, counter2) = a1;
95 a_vec_Curr_[counter1] = a_vec_Curr_[counter2] = a0;
96}
97
98// ------------Molar Thermodynamic Properties -------------------------
99
101{
103 double TKelvin = temperature();
104 double sqt = sqrt(TKelvin);
105 double mv = molarVolume();
106 double vpb = mv + m_b_current;
108 double cpref = GasConstant * mean_X(m_cp0_R);
109 double dadt = da_dt();
110 double fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
111 double dHdT_V = (cpref + mv * dpdT_ - GasConstant - 1.0 / (2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv) * fac
112 +1.0/(m_b_current * sqt) * log(vpb/mv) * (-0.5 * dadt));
113 return dHdT_V - (mv + TKelvin * dpdT_ / dpdV_) * dpdT_;
114}
115
117{
119 double TKelvin = temperature();
120 double sqt = sqrt(TKelvin);
121 double mv = molarVolume();
122 double vpb = mv + m_b_current;
123 double cvref = GasConstant * (mean_X(m_cp0_R) - 1.0);
124 double dadt = da_dt();
125 double fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
126 return (cvref - 1.0/(2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv)*fac
127 +1.0/(m_b_current * sqt) * log(vpb/mv)*(-0.5*dadt));
128}
129
131{
133
134 // Get a copy of the private variables stored in the State object
135 double T = temperature();
136 double molarV = meanMolecularWeight() / density();
137 double pp = GasConstant * T/(molarV - m_b_current) - m_a_current/(sqrt(T) * molarV * (molarV + m_b_current));
138 return pp;
139}
140
142{
144 return 1.0 / m_tmpV[k];
145}
146
148{
149 double mv = molarVolume();
150 double sqt = sqrt(temperature());
151 double vpb = mv + m_b_current;
152 double vmb = mv - m_b_current;
153
154 for (size_t k = 0; k < m_kk; k++) {
155 m_pp[k] = 0.0;
156 for (size_t i = 0; i < m_kk; i++) {
157 size_t counter = k + m_kk*i;
158 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
159 }
160 }
161 double pres = pressure();
162
163 for (size_t k = 0; k < m_kk; k++) {
164 ac[k] = (- RT() * log(pres * mv / RT())
165 + RT() * log(mv / vmb)
166 + RT() * b_vec_Curr_[k] / vmb
167 - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
168 + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
169 - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
170 );
171 }
172 for (size_t k = 0; k < m_kk; k++) {
173 ac[k] = exp(ac[k]/RT());
174 }
175}
176
177// ---- Partial Molar Properties of the Solution -----------------
178
180{
181 getGibbs_ref(mu);
182 for (size_t k = 0; k < m_kk; k++) {
183 double xx = std::max(SmallNumber, moleFraction(k));
184 mu[k] += RT()*(log(xx));
185 }
186
187 double mv = molarVolume();
188 double sqt = sqrt(temperature());
189 double vpb = mv + m_b_current;
190 double vmb = mv - m_b_current;
191
192 for (size_t k = 0; k < m_kk; k++) {
193 m_pp[k] = 0.0;
194 for (size_t i = 0; i < m_kk; i++) {
195 size_t counter = k + m_kk*i;
196 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
197 }
198 }
199 double pres = pressure();
200 double refP = refPressure();
201
202 for (size_t k = 0; k < m_kk; k++) {
203 mu[k] += (RT() * log(pres/refP) - RT() * log(pres * mv / RT())
204 + RT() * log(mv / vmb)
205 + RT() * b_vec_Curr_[k] / vmb
206 - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
207 + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
208 - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
209 );
210 }
211}
212
214{
215 // First we get the reference state contributions
216 getEnthalpy_RT_ref(hbar);
217 scale(hbar, hbar+m_kk, hbar, RT());
218
219 // We calculate dpdni_
220 double TKelvin = temperature();
221 double mv = molarVolume();
222 double sqt = sqrt(TKelvin);
223 double vpb = mv + m_b_current;
224 double vmb = mv - m_b_current;
225 for (size_t k = 0; k < m_kk; k++) {
226 m_pp[k] = 0.0;
227 for (size_t i = 0; i < m_kk; i++) {
228 size_t counter = k + m_kk*i;
229 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
230 }
231 }
232 for (size_t k = 0; k < m_kk; k++) {
233 dpdni_[k] = RT()/vmb + RT() * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb)
234 + m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb);
235 }
236 double dadt = da_dt();
237 double fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
238
239 for (size_t k = 0; k < m_kk; k++) {
240 m_tmpV[k] = 0.0;
241 for (size_t i = 0; i < m_kk; i++) {
242 size_t counter = k + m_kk*i;
243 m_tmpV[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter];
244 }
245 }
246
248 double fac2 = mv + TKelvin * dpdT_ / dpdV_;
249 for (size_t k = 0; k < m_kk; k++) {
250 double hE_v = (mv * dpdni_[k] - RT() - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac
251 + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_tmpV[k]
252 + b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac);
253 hbar[k] = hbar[k] + hE_v;
254 hbar[k] -= fac2 * dpdni_[k];
255 }
256}
257
259{
260 getEntropy_R_ref(sbar);
261 scale(sbar, sbar+m_kk, sbar, GasConstant);
262 double TKelvin = temperature();
263 double sqt = sqrt(TKelvin);
264 double mv = molarVolume();
265 double refP = refPressure();
266
267 for (size_t k = 0; k < m_kk; k++) {
268 double xx = std::max(SmallNumber, moleFraction(k));
269 sbar[k] += GasConstant * (- log(xx));
270 }
271 for (size_t k = 0; k < m_kk; k++) {
272 m_pp[k] = 0.0;
273 for (size_t i = 0; i < m_kk; i++) {
274 size_t counter = k + m_kk*i;
275 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
276 }
277 }
278 for (size_t k = 0; k < m_kk; k++) {
279 m_tmpV[k] = 0.0;
280 for (size_t i = 0; i < m_kk; i++) {
281 size_t counter = k + m_kk*i;
282 m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
283 }
284 }
285
286 double dadt = da_dt();
287 double fac = dadt - m_a_current / (2.0 * TKelvin);
288 double vmb = mv - m_b_current;
289 double vpb = mv + m_b_current;
290 for (size_t k = 0; k < m_kk; k++) {
291 sbar[k] -=(GasConstant * log(GasConstant * TKelvin / (refP * mv))
293 + GasConstant * log(mv/vmb)
294 + GasConstant * b_vec_Curr_[k]/vmb
295 + m_pp[k]/(m_b_current * TKelvin * sqt) * log(vpb/mv)
296 - 2.0 * m_tmpV[k]/(m_b_current * sqt) * log(vpb/mv)
297 + b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) * fac
298 - 1.0 / (m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac
299 );
300 }
301
303 getPartialMolarVolumes(m_partialMolarVolumes.data());
304 for (size_t k = 0; k < m_kk; k++) {
305 sbar[k] -= -m_partialMolarVolumes[k] * dpdT_;
306 }
307}
308
310{
311 // u_k = h_k - P * v_k
312 getPartialMolarVolumes(m_partialMolarVolumes.data());
314 double p = pressure();
315 for (size_t k = 0; k < nSpecies(); k++) {
316 ubar[k] -= p * m_partialMolarVolumes[k];
317 }
318}
319
321{
322 for (size_t k = 0; k < m_kk; k++) {
323 m_pp[k] = 0.0;
324 for (size_t i = 0; i < m_kk; i++) {
325 size_t counter = k + m_kk*i;
326 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
327 }
328 }
329 for (size_t k = 0; k < m_kk; k++) {
330 m_tmpV[k] = 0.0;
331 for (size_t i = 0; i < m_kk; i++) {
332 size_t counter = k + m_kk*i;
333 m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
334 }
335 }
336
337 double sqt = sqrt(temperature());
338 double mv = molarVolume();
339 double vmb = mv - m_b_current;
340 double vpb = mv + m_b_current;
341 for (size_t k = 0; k < m_kk; k++) {
342 double num = (RT() + RT() * m_b_current/ vmb + RT() * b_vec_Curr_[k] / vmb
343 + RT() * m_b_current * b_vec_Curr_[k] /(vmb * vmb)
344 - 2.0 * m_pp[k] / (sqt * vpb)
345 + m_a_current * b_vec_Curr_[k] / (sqt * vpb * vpb)
346 );
347 double denom = (pressure() + RT() * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb)
348 );
349 vbar[k] = num / denom;
350 }
351}
352
353bool RedlichKwongMFTP::addSpecies(shared_ptr<Species> spec)
354{
355 bool added = MixtureFugacityTP::addSpecies(spec);
356 if (added) {
357 a_vec_Curr_.resize(m_kk * m_kk, 0.0);
358
359 // Initialize a_vec and b_vec to NaN, to screen for species with
360 // pureFluidParameters which are undefined in the input file:
361 b_vec_Curr_.push_back(NAN);
362 a_coeff_vec.resize(2, m_kk * m_kk, NAN);
363
364 m_pp.push_back(0.0);
365 m_coeffSource.push_back(CoeffSource::EoS);
366 m_partialMolarVolumes.push_back(0.0);
367 dpdni_.push_back(0.0);
368 }
369 return added;
370}
371
373{
374 // Contents of 'critical-properties.yaml', loaded later if needed
375 AnyMap critPropsDb;
376 std::unordered_map<string, AnyMap*> dbSpecies;
377
378 for (auto& [name, species] : m_species) {
379 auto& data = species->input;
380 size_t k = speciesIndex(name);
381 if (!isnan(a_coeff_vec(0, k + m_kk * k))) {
382 continue;
383 }
384 bool foundCoeffs = false;
385
386 if (data.hasKey("equation-of-state") &&
387 data["equation-of-state"].hasMapWhere("model", "Redlich-Kwong"))
388 {
389 // Read a and b coefficients from species 'input' information (that is, as
390 // specified in a YAML input file) specific to the Redlich-Kwong model
391 auto eos = data["equation-of-state"].getMapWhere(
392 "model", "Redlich-Kwong");
393
394 if (eos.hasKey("a") && eos.hasKey("b")) {
395 double a0 = 0, a1 = 0;
396 if (eos["a"].isScalar()) {
397 a0 = eos.convert("a", "Pa*m^6/kmol^2*K^0.5");
398 } else {
399 auto avec = eos["a"].asVector<AnyValue>(2);
400 a0 = eos.units().convert(avec[0], "Pa*m^6/kmol^2*K^0.5");
401 a1 = eos.units().convert(avec[1], "Pa*m^6/kmol^2/K^0.5");
402 }
403 double b = eos.convert("b", "m^3/kmol");
404 foundCoeffs = true;
405 setSpeciesCoeffs(name, a0, a1, b);
406 m_coeffSource[k] = CoeffSource::EoS;
407 }
408
409 if (eos.hasKey("binary-a")) {
410 AnyMap& binary_a = eos["binary-a"].as<AnyMap>();
411 const UnitSystem& units = binary_a.units();
412 for (auto& [name2, coeff] : binary_a) {
413 double a0 = 0, a1 = 0;
414 if (coeff.isScalar()) {
415 a0 = units.convert(coeff, "Pa*m^6/kmol^2*K^0.5");
416 } else {
417 auto avec = coeff.asVector<AnyValue>(2);
418 a0 = units.convert(avec[0], "Pa*m^6/kmol^2*K^0.5");
419 a1 = units.convert(avec[1], "Pa*m^6/kmol^2/K^0.5");
420 }
421 setBinaryCoeffs(name, name2, a0, a1);
422 }
423 }
424
425 if (foundCoeffs) {
426 continue;
427 }
428 }
429
430 // Coefficients have not been populated from model-specific input
431 double Tc = NAN, Pc = NAN;
432 if (data.hasKey("critical-parameters")) {
433 // Use critical state information stored in the species entry to
434 // calculate a and b
435 auto& critProps = data["critical-parameters"].as<AnyMap>();
436 Tc = critProps.convert("critical-temperature", "K");
437 Pc = critProps.convert("critical-pressure", "Pa");
438 m_coeffSource[k] = CoeffSource::CritProps;
439 } else {
440 // Search 'crit-properties.yaml' to find Tc and Pc. Load data if needed
441 if (critPropsDb.empty()) {
442 critPropsDb = AnyMap::fromYamlFile("critical-properties.yaml");
443 dbSpecies = critPropsDb["species"].asMap("name");
444 }
445
446 // All names in critical-properties.yaml are upper case
447 auto ucName = boost::algorithm::to_upper_copy(name);
448 if (dbSpecies.count(ucName)) {
449 auto& spec = *dbSpecies.at(ucName);
450 auto& critProps = spec["critical-parameters"].as<AnyMap>();
451 Tc = critProps.convert("critical-temperature", "K");
452 Pc = critProps.convert("critical-pressure", "Pa");
453 m_coeffSource[k] = CoeffSource::Database;
454 }
455 }
456
457 // Check if critical properties were found in either location
458 if (!isnan(Tc)) {
459 // Assuming no temperature dependence (that is, a1 = 0)
460 double a = omega_a * pow(GasConstant, 2) * pow(Tc, 2.5) / Pc;
461 double b = omega_b * GasConstant * Tc / Pc;
462 setSpeciesCoeffs(name, a, 0.0, b);
463 } else {
464 throw InputFileError("RedlichKwongMFTP::initThermo", data,
465 "No critical property or Redlich-Kwong parameters found "
466 "for species {}.", name);
467 }
468 }
469}
470
472 AnyMap& speciesNode) const
473{
475 size_t k = speciesIndex(name);
477 if (m_coeffSource[k] == CoeffSource::EoS) {
478 auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
479 "model", "Redlich-Kwong", true);
480
481 size_t counter = k + m_kk * k;
482 if (a_coeff_vec(1, counter) != 0.0) {
483 vector<AnyValue> coeffs(2);
484 coeffs[0].setQuantity(a_coeff_vec(0, counter), "Pa*m^6/kmol^2*K^0.5");
485 coeffs[1].setQuantity(a_coeff_vec(1, counter), "Pa*m^6/kmol^2/K^0.5");
486 eosNode["a"] = std::move(coeffs);
487 } else {
488 eosNode["a"].setQuantity(a_coeff_vec(0, counter),
489 "Pa*m^6/kmol^2*K^0.5");
490 }
491 eosNode["b"].setQuantity(b_vec_Curr_[k], "m^3/kmol");
492 } else if (m_coeffSource[k] == CoeffSource::CritProps) {
493 auto& critProps = speciesNode["critical-parameters"];
494 double a = a_coeff_vec(0, k + m_kk * k);
495 double b = b_vec_Curr_[k];
496 double Tc = pow(a * omega_b / (b * omega_a * GasConstant), 2.0/3.0);
497 double Pc = omega_b * GasConstant * Tc / b;
498 critProps["critical-temperature"].setQuantity(Tc, "K");
499 critProps["critical-pressure"].setQuantity(Pc, "Pa");
500 }
501
502 if (m_binaryParameters.count(name)) {
503 auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
504 "model", "Redlich-Kwong", true);
505 AnyMap bin_a;
506 for (const auto& [name2, coeffs] : m_binaryParameters.at(name)) {
507 if (coeffs.second == 0) {
508 bin_a[name2].setQuantity(coeffs.first, "Pa*m^6/kmol^2*K^0.5");
509 } else {
510 vector<AnyValue> C(2);
511 C[0].setQuantity(coeffs.first, "Pa*m^6/kmol^2*K^0.5");
512 C[1].setQuantity(coeffs.second, "Pa*m^6/kmol^2/K^0.5");
513 bin_a[name2] = std::move(C);
514 }
515 }
516 eosNode["binary-a"] = std::move(bin_a);
517 }
518}
519
521{
522 // note this agrees with tpx
523 double rho = density();
524 double mmw = meanMolecularWeight();
525 double molarV = mmw / rho;
526 double hh = m_b_current / molarV;
527 double zz = z();
528 double dadt = da_dt();
529 double T = temperature();
530 double sqT = sqrt(T);
531 double fac = dadt - m_a_current / (2.0 * T);
532 double sresid_mol_R = log(zz*(1.0 - hh)) + log(1.0 + hh) * fac / (sqT * GasConstant * m_b_current);
533 return GasConstant * sresid_mol_R;
534}
535
537{
538 // note this agrees with tpx
539 double rho = density();
540 double mmw = meanMolecularWeight();
541 double molarV = mmw / rho;
542 double hh = m_b_current / molarV;
543 double zz = z();
544 double dadt = da_dt();
545 double T = temperature();
546 double sqT = sqrt(T);
547 double fac = T * dadt - 3.0 *m_a_current / (2.0);
548 return GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqT * m_b_current);
549}
550
551double RedlichKwongMFTP::liquidVolEst(double TKelvin, double& presGuess) const
552{
553 double v = m_b_current * 1.1;
554 double atmp;
555 double btmp;
556 calculateAB(TKelvin, atmp, btmp);
557 double pres = std::max(psatEst(TKelvin), presGuess);
558 double Vroot[3];
559 bool foundLiq = false;
560 int m = 0;
561 while (m < 100 && !foundLiq) {
562 int nsol = solveCubic(TKelvin, pres, atmp, btmp, Vroot);
563 if (nsol == 1 || nsol == 2) {
564 double pc = critPressure();
565 if (pres > pc) {
566 foundLiq = true;
567 }
568 pres *= 1.04;
569 } else {
570 foundLiq = true;
571 }
572 }
573
574 if (foundLiq) {
575 v = Vroot[0];
576 presGuess = pres;
577 } else {
578 v = -1.0;
579 }
580 return v;
581}
582
583double RedlichKwongMFTP::densityCalc(double TKelvin, double presPa, int phaseRequested,
584 double rhoguess)
585{
586 // It's necessary to set the temperature so that m_a_current is set correctly.
587 setTemperature(TKelvin);
588 double tcrit = critTemperature();
589 double mmw = meanMolecularWeight();
590 if (rhoguess == -1.0) {
591 if (phaseRequested >= FLUID_LIQUID_0) {
592 double lqvol = liquidVolEst(TKelvin, presPa);
593 rhoguess = mmw / lqvol;
594 }
595 } else {
596 // Assume the Gas phase initial guess, if nothing is specified to the routine
597 rhoguess = presPa * mmw / (GasConstant * TKelvin);
598 }
599
600
601 double volguess = mmw / rhoguess;
602 NSolns_ = solveCubic(TKelvin, presPa, m_a_current, m_b_current, Vroot_);
603
604 double molarVolLast = Vroot_[0];
605 if (NSolns_ >= 2) {
606 if (phaseRequested >= FLUID_LIQUID_0) {
607 molarVolLast = Vroot_[0];
608 } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
609 molarVolLast = Vroot_[2];
610 } else {
611 if (volguess > Vroot_[1]) {
612 molarVolLast = Vroot_[2];
613 } else {
614 molarVolLast = Vroot_[0];
615 }
616 }
617 } else if (NSolns_ == 1) {
618 if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) {
619 molarVolLast = Vroot_[0];
620 } else {
621 return -2.0;
622 }
623 } else if (NSolns_ == -1) {
624 if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) {
625 molarVolLast = Vroot_[0];
626 } else if (TKelvin > tcrit) {
627 molarVolLast = Vroot_[0];
628 } else {
629 return -2.0;
630 }
631 } else {
632 molarVolLast = Vroot_[0];
633 return -1.0;
634 }
635 return mmw / molarVolLast;
636}
637
639{
640 double Vroot[3];
641 double T = temperature();
642 int nsol = solveCubic(T, pressure(), m_a_current, m_b_current, Vroot);
643 if (nsol != 3) {
644 return critDensity();
645 }
646
647 auto resid = [this, T](double v) {
648 double pp;
649 return dpdVCalc(T, v, pp);
650 };
651
652 boost::uintmax_t maxiter = 100;
653 auto [lower, upper] = bmt::toms748_solve(
654 resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
655
656 double mmw = meanMolecularWeight();
657 return mmw / (0.5 * (lower + upper));
658}
659
661{
662 double Vroot[3];
663 double T = temperature();
664 int nsol = solveCubic(T, pressure(), m_a_current, m_b_current, Vroot);
665 if (nsol != 3) {
666 return critDensity();
667 }
668
669 auto resid = [this, T](double v) {
670 double pp;
671 return dpdVCalc(T, v, pp);
672 };
673
674 boost::uintmax_t maxiter = 100;
675 auto [lower, upper] = bmt::toms748_solve(
676 resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
677
678 double mmw = meanMolecularWeight();
679 return mmw / (0.5 * (lower + upper));
680}
681
682double RedlichKwongMFTP::dpdVCalc(double TKelvin, double molarVol, double& presCalc) const
683{
684 double sqt = sqrt(TKelvin);
685 presCalc = GasConstant * TKelvin / (molarVol - m_b_current)
686 - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
687
688 double vpb = molarVol + m_b_current;
689 double vmb = molarVol - m_b_current;
690 double dpdv = (- GasConstant * TKelvin / (vmb * vmb)
691 + m_a_current * (2 * molarVol + m_b_current) / (sqt * molarVol * molarVol * vpb * vpb));
692 return dpdv;
693}
694
696{
698 return -1 / (molarVolume() * dpdV_);
699}
700
702{
704 return -dpdT_ / (molarVolume() * dpdV_);
705}
706
708{
710 return molarVolume() * sqrt(-cp_mole() / cv_mole() * dpdV_ / meanMolecularWeight());
711}
712
714{
715 double TKelvin = temperature();
716 double mv = molarVolume();
717 double pres;
718
719 dpdV_ = dpdVCalc(TKelvin, mv, pres);
720 double sqt = sqrt(TKelvin);
721 double vpb = mv + m_b_current;
722 double vmb = mv - m_b_current;
723 double dadt = da_dt();
724 double fac = dadt - m_a_current/(2.0 * TKelvin);
725 dpdT_ = (GasConstant / vmb - fac / (sqt * mv * vpb));
726}
727
729{
730 double temp = temperature();
731 if (m_formTempParam == 1) {
732 for (size_t i = 0; i < m_kk; i++) {
733 for (size_t j = 0; j < m_kk; j++) {
734 size_t counter = i * m_kk + j;
735 a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
736 }
737 }
738 }
739
740 m_b_current = 0.0;
741 m_a_current = 0.0;
742 for (size_t i = 0; i < m_kk; i++) {
743 m_b_current += moleFractions_[i] * b_vec_Curr_[i];
744 for (size_t j = 0; j < m_kk; j++) {
745 m_a_current += a_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j];
746 }
747 }
748 if (isnan(m_b_current)) {
749 // One or more species do not have specified coefficients.
750 fmt::memory_buffer b;
751 for (size_t k = 0; k < m_kk; k++) {
752 if (isnan(b_vec_Curr_[k])) {
753 if (b.size() > 0) {
754 fmt_append(b, ", {}", speciesName(k));
755 } else {
756 fmt_append(b, "{}", speciesName(k));
757 }
758 }
759 }
760 throw CanteraError("RedlichKwongMFTP::updateMixingExpressions",
761 "Missing Redlich-Kwong coefficients for species: {}", to_string(b));
762 }
763}
764
765void RedlichKwongMFTP::calculateAB(double temp, double& aCalc, double& bCalc) const
766{
767 bCalc = 0.0;
768 aCalc = 0.0;
769 if (m_formTempParam == 1) {
770 for (size_t i = 0; i < m_kk; i++) {
771 bCalc += moleFractions_[i] * b_vec_Curr_[i];
772 for (size_t j = 0; j < m_kk; j++) {
773 size_t counter = i * m_kk + j;
774 double a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
775 aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
776 }
777 }
778 } else {
779 for (size_t i = 0; i < m_kk; i++) {
780 bCalc += moleFractions_[i] * b_vec_Curr_[i];
781 for (size_t j = 0; j < m_kk; j++) {
782 size_t counter = i * m_kk + j;
783 double a_vec_Curr = a_coeff_vec(0,counter);
784 aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
785 }
786 }
787 }
788}
789
790double RedlichKwongMFTP::da_dt() const
791{
792 double dadT = 0.0;
793 if (m_formTempParam == 1) {
794 for (size_t i = 0; i < m_kk; i++) {
795 for (size_t j = 0; j < m_kk; j++) {
796 size_t counter = i * m_kk + j;
797 dadT+= a_coeff_vec(1,counter) * moleFractions_[i] * moleFractions_[j];
798 }
799 }
800 }
801 return dadT;
802}
803
804void RedlichKwongMFTP::calcCriticalConditions(double& pc, double& tc, double& vc) const
805{
806 double a0 = 0.0;
807 double aT = 0.0;
808 for (size_t i = 0; i < m_kk; i++) {
809 for (size_t j = 0; j <m_kk; j++) {
810 size_t counter = i + m_kk * j;
811 a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
812 aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
813 }
814 }
815 double a = m_a_current;
816 double b = m_b_current;
817 if (m_formTempParam != 0) {
818 a = a0;
819 }
820 if (b <= 0.0) {
821 tc = 1000000.;
822 pc = 1.0E13;
823 vc = omega_vc * GasConstant * tc / pc;
824 return;
825 }
826 if (a <= 0.0) {
827 tc = 0.0;
828 pc = 0.0;
829 vc = 2.0 * b;
830 return;
831 }
832 double tmp = a * omega_b / (b * omega_a * GasConstant);
833 double pp = 2./3.;
834 double sqrttc, f, dfdt, deltatc;
835
836 if (m_formTempParam == 0) {
837 tc = pow(tmp, pp);
838 } else {
839 tc = pow(tmp, pp);
840 for (int j = 0; j < 10; j++) {
841 sqrttc = sqrt(tc);
842 f = omega_a * b * GasConstant * tc * sqrttc / omega_b - aT * tc - a0;
843 dfdt = 1.5 * omega_a * b * GasConstant * sqrttc / omega_b - aT;
844 deltatc = - f / dfdt;
845 tc += deltatc;
846 }
847 if (deltatc > 0.1) {
848 throw CanteraError("RedlichKwongMFTP::calcCriticalConditions",
849 "didn't converge");
850 }
851 }
852
853 pc = omega_b * GasConstant * tc / b;
854 vc = omega_vc * GasConstant * tc / pc;
855}
856
857int RedlichKwongMFTP::solveCubic(double T, double pres, double a, double b, double Vroot[3]) const
858{
859
860 // Derive the coefficients of the cubic polynomial to solve.
861 double an = 1.0;
862 double bn = - GasConstant * T / pres;
863 double sqt = sqrt(T);
864 double cn = - (GasConstant * T * b / pres - a/(pres * sqt) + b * b);
865 double dn = - (a * b / (pres * sqt));
866
867 double tmp = a * omega_b / (b * omega_a * GasConstant);
868 double pp = 2./3.;
869 double tc = pow(tmp, pp);
870 double pc = omega_b * GasConstant * tc / b;
871 double vc = omega_vc * GasConstant * tc / pc;
872
873 int nSolnValues = MixtureFugacityTP::solveCubic(T, pres, a, b, a, Vroot, an, bn, cn, dn, tc, vc);
874
875 return nSolnValues;
876}
877
878}
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:431
const UnitSystem & units() const
Return the default units that should be used to convert stored values.
Definition AnyMap.h:640
bool empty() const
Return boolean indicating whether AnyMap is empty.
Definition AnyMap.cpp:1468
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:1595
static AnyMap fromYamlFile(const string &name, const string &parent_name="")
Create an AnyMap from a YAML file.
Definition AnyMap.cpp:1841
A wrapper for a variable whose type is determined at runtime.
Definition AnyMap.h:87
void getRow(size_t n, double *const rw)
Get the nth row and return it in a vector.
Definition Array.cpp:79
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:749
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 getEntropy_R_ref(double *er) const override
Returns the vector of nondimensional entropies of the reference state at the current temperature of t...
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 ...
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
Definition Phase.cpp:153
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:231
size_t m_kk
Number of species in the phase.
Definition Phase.h:854
double temperature() const
Temperature (K).
Definition Phase.h:562
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:655
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
map< string, shared_ptr< Species > > m_species
Map of Species objects.
Definition Phase.h:867
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
virtual double density() const
Density (kg/m^3).
Definition Phase.h:587
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:873
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 dpdT_
The derivative of the pressure wrt the temperature.
void calculateAB(double temp, double &aCalc, double &bCalc) const
Calculate the a and the b parameters given the temperature.
double thermalExpansionCoeff() const override
Return the volumetric thermal expansion coefficient. Units: 1/K.
void setBinaryCoeffs(const string &species_i, const string &species_j, double a0, double a1)
Set values for the interaction parameter between two species.
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.
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).
int m_formTempParam
Form of the temperature parameterization.
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...
RedlichKwongMFTP(const string &infile="", const string &id="")
Construct a RedlichKwongMFTP object from an input file.
static const double omega_b
Omega constant for b.
vector< double > m_pp
Temporary storage - length = m_kk.
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...
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.
void pressureDerivatives() const
Calculate dpdV and dpdT 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.
static const double omega_a
Omega constant for a -> value of a in terms of critical properties.
double liquidVolEst(double TKelvin, double &pres) const override
Estimate for the molar volume of the liquid.
double dpdVCalc(double TKelvin, double molarVol, double &presCalc) const override
Calculate the pressure and the pressure derivative given the temperature and the molar volume.
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 and b coefficients.
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 a and b parameters.
double cp_mole() const override
Molar heat capacity at constant pressure. Units: J/kmol/K.
double m_a_current
Value of a in the equation of state.
double standardConcentration(size_t k=0) const override
Returns the standard concentration , which is used to normalize the generalized concentration.
double dpdV_
The derivative of the pressure wrt the volume.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
double m_b_current
Value of b in the equation of state.
void getActivityCoefficients(double *ac) const override
Get the array of non-dimensional activity coefficients at the current solution temperature,...
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.
vector< double > dpdni_
Vector of derivatives of pressure wrt mole number.
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
map< string, map< string, pair< double, double > > > m_binaryParameters
Explicitly-specified binary interaction parameters.
void setSpeciesCoeffs(const string &species, double a0, double a1, double b)
Set the pure fluid interaction parameters for a species.
int solveCubic(double T, double pres, double a, double b, double Vroot[3]) const
Prepare variables and call the function to solve the cubic equation of state.
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 fmt_append(fmt::memory_buffer &b, const std::string &tmpl, Args... args)
Versions 6.2.0 and 6.2.1 of fmtlib do not include this define before they include windows....
Definition fmt.h:29
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
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
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...