Cantera  3.2.0a2
Loading...
Searching...
No Matches
MixtureFugacityTP.cpp
Go to the documentation of this file.
1/**
2 * @file MixtureFugacityTP.cpp
3 * Methods file for a derived class of ThermoPhase that handles
4 * non-ideal mixtures based on the fugacity models (see @ref thermoprops and
5 * class @link Cantera::MixtureFugacityTP MixtureFugacityTP@endlink).
6 */
7
8// This file is part of Cantera. See License.txt in the top-level directory or
9// at https://cantera.org/license.txt for license and copyright information.
10
14#include "cantera/base/global.h"
15
16using namespace std;
17
18namespace Cantera
19{
20
22{
24}
25
27{
28 forcedState_ = solnBranch;
29}
30
32{
33 return forcedState_;
34}
35
37{
38 return iState_;
39}
40
41// ---- Molar Thermodynamic Properties ---------------------------
43{
44 double h_ideal = RT() * mean_X(m_h0_RT);
45 double h_nonideal = hresid();
46 return h_ideal + h_nonideal;
47}
48
49
51{
52 double s_ideal = GasConstant * (mean_X(m_s0_R) - sum_xlogx()
53 - std::log(pressure()/refPressure()));
54 double s_nonideal = sresid();
55 return s_ideal + s_nonideal;
56}
57
58// ----- Thermodynamic Values for the Species Standard States States ----
59
61{
62 copy(m_g0_RT.begin(), m_g0_RT.end(), g);
63 double tmp = log(pressure() / refPressure());
64 for (size_t k = 0; k < m_kk; k++) {
65 g[k] = RT() * (g[k] + tmp);
66 }
67}
68
69void MixtureFugacityTP::getEnthalpy_RT(double* hrt) const
70{
72}
73
74void MixtureFugacityTP::getEntropy_R(double* sr) const
75{
76 copy(m_s0_R.begin(), m_s0_R.end(), sr);
77 double tmp = log(pressure() / refPressure());
78 for (size_t k = 0; k < m_kk; k++) {
79 sr[k] -= tmp;
80 }
81}
82
83void MixtureFugacityTP::getGibbs_RT(double* grt) const
84{
85 copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
86 double tmp = log(pressure() / refPressure());
87 for (size_t k = 0; k < m_kk; k++) {
88 grt[k] += tmp;
89 }
90}
91
93{
94 scale(m_g0_RT.begin(), m_g0_RT.end(), g, RT());
95 double tmp = log(pressure() / refPressure()) * RT();
96 for (size_t k = 0; k < m_kk; k++) {
97 g[k] += tmp;
98 }
99}
100
102{
103 copy(m_h0_RT.begin(), m_h0_RT.end(), urt);
104 for (size_t i = 0; i < m_kk; i++) {
105 urt[i] -= 1.0;
106 }
107}
108
109void MixtureFugacityTP::getCp_R(double* cpr) const
110{
111 copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
112}
113
115{
116 for (size_t i = 0; i < m_kk; i++) {
117 vol[i] = RT() / pressure();
118 }
119}
120
121// ----- Thermodynamic Values for the Species Reference States ----
122
124{
125 copy(m_h0_RT.begin(), m_h0_RT.end(), hrt);
126}
127
129{
130 copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
131}
132
134{
135 const vector<double>& gibbsrt = gibbs_RT_ref();
136 scale(gibbsrt.begin(), gibbsrt.end(), g, RT());
137}
138
139const vector<double>& MixtureFugacityTP::gibbs_RT_ref() const
140{
141 return m_g0_RT;
142}
143
145{
146 copy(m_s0_R.begin(), m_s0_R.end(), er);
147}
148
149void MixtureFugacityTP::getCp_R_ref(double* cpr) const
150{
151 copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
152}
153
155{
156 for (size_t i = 0; i < m_kk; i++) {
157 vol[i]= RT() / refPressure();
158 }
159}
160
161bool MixtureFugacityTP::addSpecies(shared_ptr<Species> spec)
162{
163 bool added = ThermoPhase::addSpecies(spec);
164 if (added) {
165 if (m_kk == 1) {
166 moleFractions_.push_back(1.0);
167 } else {
168 moleFractions_.push_back(0.0);
169 }
170 m_h0_RT.push_back(0.0);
171 m_cp0_R.push_back(0.0);
172 m_g0_RT.push_back(0.0);
173 m_s0_R.push_back(0.0);
174 }
175 return added;
176}
177
179{
182 // depends on mole fraction and temperature
183 updateMixingExpressions();
184 iState_ = phaseState(true);
185}
186
188{
189 // A pretty tricky algorithm is needed here, due to problems involving
190 // standard states of real fluids. For those cases you need to combine the T
191 // and P specification for the standard state, or else you may venture into
192 // the forbidden zone, especially when nearing the triple point. Therefore,
193 // we need to do the standard state thermo calc with the (t, pres) combo.
194
195 double t = temperature();
196 double rhoNow = density();
197 if (forcedState_ == FLUID_UNDEFINED) {
198 double rho = densityCalc(t, p, iState_, rhoNow);
199 if (rho > 0.0) {
200 setDensity(rho);
201 iState_ = phaseState(true);
202 } else {
203 if (rho < -1.5) {
204 rho = densityCalc(t, p, FLUID_UNDEFINED , rhoNow);
205 if (rho > 0.0) {
206 setDensity(rho);
207 iState_ = phaseState(true);
208 } else {
209 throw CanteraError("MixtureFugacityTP::setPressure",
210 "neg rho");
211 }
212 } else {
213 throw CanteraError("MixtureFugacityTP::setPressure",
214 "neg rho");
215 }
216 }
217 } else if (forcedState_ == FLUID_GAS) {
218 // Normal density calculation
219 if (iState_ < FLUID_LIQUID_0) {
220 double rho = densityCalc(t, p, iState_, rhoNow);
221 if (rho > 0.0) {
222 setDensity(rho);
223 iState_ = phaseState(true);
224 if (iState_ >= FLUID_LIQUID_0) {
225 throw CanteraError("MixtureFugacityTP::setPressure",
226 "wrong state");
227 }
228 } else {
229 throw CanteraError("MixtureFugacityTP::setPressure",
230 "neg rho");
231 }
232 }
233 } else if (forcedState_ > FLUID_LIQUID_0) {
234 if (iState_ >= FLUID_LIQUID_0) {
235 double rho = densityCalc(t, p, iState_, rhoNow);
236 if (rho > 0.0) {
237 setDensity(rho);
238 iState_ = phaseState(true);
239 if (iState_ == FLUID_GAS) {
240 throw CanteraError("MixtureFugacityTP::setPressure",
241 "wrong state");
242 }
243 } else {
244 throw CanteraError("MixtureFugacityTP::setPressure",
245 "neg rho");
246 }
247 }
248 }
249}
250
252{
255 updateMixingExpressions();
256}
257
259{
261 double p_RT = pressure() / RT();
262 for (size_t k = 0; k < m_kk; k++) {
263 c[k] *= moleFraction(k)*p_RT;
264 }
265}
266
268{
269 return pressure() * meanMolecularWeight() / (density() * RT());
270}
271
273{
274 throw NotImplementedError("MixtureFugacityTP::sresid");
275}
276
278{
279 throw NotImplementedError("MixtureFugacityTP::hresid");
280}
281
282double MixtureFugacityTP::psatEst(double TKelvin) const
283{
284 double pcrit = critPressure();
285 double tt = critTemperature() / TKelvin;
286 if (tt < 1.0) {
287 return pcrit;
288 }
289 double lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918;
290 return pcrit*exp(lpr);
291}
292
293double MixtureFugacityTP::liquidVolEst(double TKelvin, double& pres) const
294{
295 throw NotImplementedError("MixtureFugacityTP::liquidVolEst");
296}
297
298double MixtureFugacityTP::densityCalc(double TKelvin, double presPa,
299 int phase, double rhoguess)
300{
301 double tcrit = critTemperature();
302 double mmw = meanMolecularWeight();
303 if (rhoguess == -1.0) {
304 if (phase != -1) {
305 if (TKelvin > tcrit) {
306 rhoguess = presPa * mmw / (GasConstant * TKelvin);
307 } else {
308 if (phase == FLUID_GAS || phase == FLUID_SUPERCRIT) {
309 rhoguess = presPa * mmw / (GasConstant * TKelvin);
310 } else if (phase >= FLUID_LIQUID_0) {
311 double lqvol = liquidVolEst(TKelvin, presPa);
312 rhoguess = mmw / lqvol;
313 }
314 }
315 } else {
316 // Assume the Gas phase initial guess, if nothing is specified to
317 // the routine
318 rhoguess = presPa * mmw / (GasConstant * TKelvin);
319 }
320 }
321
322 double molarVolBase = mmw / rhoguess;
323 double molarVolLast = molarVolBase;
324 double vc = mmw / critDensity();
325
326 // molar volume of the spinodal at the current temperature and mole
327 // fractions. this will be updated as we go.
328 double molarVolSpinodal = vc;
329 bool conv = false;
330
331 // We start on one side of the vc and stick with that side
332 bool gasSide = molarVolBase > vc;
333 if (gasSide) {
334 molarVolLast = (GasConstant * TKelvin)/presPa;
335 } else {
336 molarVolLast = liquidVolEst(TKelvin, presPa);
337 }
338
339 // OK, now we do a small solve to calculate the molar volume given the T,P
340 // value. The algorithm is taken from dfind()
341 for (int n = 0; n < 200; n++) {
342 // Calculate the predicted reduced pressure, pred0, based on the current
343 // tau and dd. Calculate the derivative of the predicted pressure wrt
344 // the molar volume. This routine also returns the pressure, presBase
345 double presBase;
346 double dpdVBase = dpdVCalc(TKelvin, molarVolBase, presBase);
347
348 // If dpdV is positive, then we are in the middle of the 2 phase region
349 // and beyond the spinodal stability curve. We need to adjust the
350 // initial guess outwards and start a new iteration.
351 if (dpdVBase >= 0.0) {
352 if (TKelvin > tcrit) {
353 throw CanteraError("MixtureFugacityTP::densityCalc",
354 "T > tcrit unexpectedly");
355 }
356
357 // TODO Spawn a calculation for the value of the spinodal point that
358 // is very accurate. Answer the question as to whether a
359 // solution is possible on the current side of the vapor dome.
360 if (gasSide) {
361 if (molarVolBase >= vc) {
362 molarVolSpinodal = molarVolBase;
363 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
364 } else {
365 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
366 }
367 } else {
368 if (molarVolBase <= vc) {
369 molarVolSpinodal = molarVolBase;
370 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
371 } else {
372 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
373 }
374 }
375 continue;
376 }
377
378 // Check for convergence
379 if (fabs(presBase-presPa) < 1.0E-30 + 1.0E-8 * presPa) {
380 conv = true;
381 break;
382 }
383
384 // Dampen and crop the update
385 double dpdV = dpdVBase;
386 if (n < 10) {
387 dpdV = dpdVBase * 1.5;
388 }
389
390 // Formulate the update to the molar volume by Newton's method. Then,
391 // crop it to a max value of 0.1 times the current volume
392 double delMV = - (presBase - presPa) / dpdV;
393 if ((!gasSide || delMV < 0.0) && fabs(delMV) > 0.2 * molarVolBase) {
394 delMV = delMV / fabs(delMV) * 0.2 * molarVolBase;
395 }
396 // Only go 1/10 the way towards the spinodal at any one time.
397 if (TKelvin < tcrit) {
398 if (gasSide) {
399 if (delMV < 0.0 && -delMV > 0.5 * (molarVolBase - molarVolSpinodal)) {
400 delMV = - 0.5 * (molarVolBase - molarVolSpinodal);
401 }
402 } else {
403 if (delMV > 0.0 && delMV > 0.5 * (molarVolSpinodal - molarVolBase)) {
404 delMV = 0.5 * (molarVolSpinodal - molarVolBase);
405 }
406 }
407 }
408 // updated the molar volume value
409 molarVolLast = molarVolBase;
410 molarVolBase += delMV;
411
412 if (fabs(delMV/molarVolBase) < 1.0E-14) {
413 conv = true;
414 break;
415 }
416
417 // Check for negative molar volumes
418 if (molarVolBase <= 0.0) {
419 molarVolBase = std::min(1.0E-30, fabs(delMV*1.0E-4));
420 }
421 }
422
423 // Check for convergence, and return 0.0 if it wasn't achieved.
424 double densBase = 0.0;
425 if (! conv) {
426 molarVolBase = 0.0;
427 throw CanteraError("MixtureFugacityTP::densityCalc",
428 "Process did not converge");
429 } else {
430 densBase = mmw / molarVolBase;
431 }
432 return densBase;
433}
434
435void MixtureFugacityTP::updateMixingExpressions()
436{
437}
438
439int MixtureFugacityTP::corr0(double TKelvin, double pres, double& densLiqGuess,
440 double& densGasGuess, double& liqGRT, double& gasGRT)
441{
442 int retn = 0;
443 double densLiq = densityCalc(TKelvin, pres, FLUID_LIQUID_0, densLiqGuess);
444 if (densLiq <= 0.0) {
445 retn = -1;
446 } else {
447 densLiqGuess = densLiq;
448 setState_TD(TKelvin, densLiq);
449 liqGRT = gibbs_mole() / RT();
450 }
451
452 double densGas = densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess);
453 if (densGas <= 0.0) {
454 if (retn == -1) {
455 throw CanteraError("MixtureFugacityTP::corr0",
456 "Error occurred trying to find gas density at (T,P) = {} {}",
457 TKelvin, pres);
458 }
459 retn = -2;
460 } else {
461 densGasGuess = densGas;
462 setState_TD(TKelvin, densGas);
463 gasGRT = gibbs_mole() / RT();
464 }
465 return retn;
466}
467
468int MixtureFugacityTP::phaseState(bool checkState) const
469{
470 int state = iState_;
471 if (checkState) {
472 double t = temperature();
473 double tcrit = critTemperature();
474 double rhocrit = critDensity();
475 if (t >= tcrit) {
476 return FLUID_SUPERCRIT;
477 }
478 double tmid = tcrit - 100.;
479 if (tmid < 0.0) {
480 tmid = tcrit / 2.0;
481 }
482 double pp = psatEst(tmid);
483 double mmw = meanMolecularWeight();
484 double molVolLiqTmid = liquidVolEst(tmid, pp);
485 double molVolGasTmid = GasConstant * tmid / pp;
486 double densLiqTmid = mmw / molVolLiqTmid;
487 double densGasTmid = mmw / molVolGasTmid;
488 double densMidTmid = 0.5 * (densLiqTmid + densGasTmid);
489 double rhoMid = rhocrit + (t - tcrit) * (rhocrit - densMidTmid) / (tcrit - tmid);
490
491 double rho = density();
492 int iStateGuess = FLUID_LIQUID_0;
493 if (rho < rhoMid) {
494 iStateGuess = FLUID_GAS;
495 }
496 double molarVol = mmw / rho;
497 double presCalc;
498
499 double dpdv = dpdVCalc(t, molarVol, presCalc);
500 if (dpdv < 0.0) {
501 state = iStateGuess;
502 } else {
503 state = FLUID_UNSTABLE;
504 }
505 }
506 return state;
507}
508
510{
511 throw NotImplementedError("MixtureFugacityTP::densSpinodalLiquid");
512}
513
515{
516 throw NotImplementedError("MixtureFugacityTP::densSpinodalGas");
517}
518
519double MixtureFugacityTP::satPressure(double TKelvin)
520{
521 double molarVolGas;
522 double molarVolLiquid;
523 return calculatePsat(TKelvin, molarVolGas, molarVolLiquid);
524}
525
526double MixtureFugacityTP::calculatePsat(double TKelvin, double& molarVolGas,
527 double& molarVolLiquid)
528{
529 // The algorithm for this routine has undergone quite a bit of work. It
530 // probably needs more work. However, it seems now to be fairly robust. The
531 // key requirement is to find an initial pressure where both the liquid and
532 // the gas exist. This is not as easy as it sounds, and it gets exceedingly
533 // hard as the critical temperature is approached from below. Once we have
534 // this initial state, then we seek to equilibrate the Gibbs free energies
535 // of the gas and liquid and use the formula
536 //
537 // dp = VdG
538 //
539 // to create an update condition for deltaP using
540 //
541 // - (Gliq - Ggas) = (Vliq - Vgas) (deltaP)
542 //
543 // @todo Suggestions for the future would be to switch it to an algorithm
544 // that uses the gas molar volume and the liquid molar volumes as the
545 // fundamental unknowns.
546
547 // we need this because this is a non-const routine that is public
548 setTemperature(TKelvin);
549 double densSave = density();
550 double tempSave = temperature();
551 double pres;
552 double mw = meanMolecularWeight();
553 if (TKelvin < critTemperature()) {
554 pres = psatEst(TKelvin);
555 // trial value = Psat from correlation
556 double volLiquid = liquidVolEst(TKelvin, pres);
557 double RhoLiquidGood = mw / volLiquid;
558 double RhoGasGood = pres * mw / (GasConstant * TKelvin);
559 double delGRT = 1.0E6;
560 double liqGRT, gasGRT;
561
562 // First part of the calculation involves finding a pressure at which
563 // the gas and the liquid state coexists.
564 double presLiquid = 0.;
565 double presGas;
566 double presBase = pres;
567 bool foundLiquid = false;
568 bool foundGas = false;
569
570 double densLiquid = densityCalc(TKelvin, presBase, FLUID_LIQUID_0, RhoLiquidGood);
571 if (densLiquid > 0.0) {
572 foundLiquid = true;
573 presLiquid = pres;
574 RhoLiquidGood = densLiquid;
575 }
576 if (!foundLiquid) {
577 for (int i = 0; i < 50; i++) {
578 pres = 1.1 * pres;
579 densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
580 if (densLiquid > 0.0) {
581 foundLiquid = true;
582 presLiquid = pres;
583 RhoLiquidGood = densLiquid;
584 break;
585 }
586 }
587 }
588
589 pres = presBase;
590 double densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
591 if (densGas <= 0.0) {
592 foundGas = false;
593 } else {
594 foundGas = true;
595 presGas = pres;
596 RhoGasGood = densGas;
597 }
598 if (!foundGas) {
599 for (int i = 0; i < 50; i++) {
600 pres = 0.9 * pres;
601 densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
602 if (densGas > 0.0) {
603 foundGas = true;
604 presGas = pres;
605 RhoGasGood = densGas;
606 break;
607 }
608 }
609 }
610
611 if (foundGas && foundLiquid && presGas != presLiquid) {
612 pres = 0.5 * (presLiquid + presGas);
613 bool goodLiq;
614 bool goodGas;
615 for (int i = 0; i < 50; i++) {
616 densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
617 if (densLiquid <= 0.0) {
618 goodLiq = false;
619 } else {
620 goodLiq = true;
621 RhoLiquidGood = densLiquid;
622 presLiquid = pres;
623 }
624 densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
625 if (densGas <= 0.0) {
626 goodGas = false;
627 } else {
628 goodGas = true;
629 RhoGasGood = densGas;
630 presGas = pres;
631 }
632 if (goodGas && goodLiq) {
633 break;
634 }
635 if (!goodLiq && !goodGas) {
636 pres = 0.5 * (pres + presLiquid);
637 }
638 if (goodLiq || goodGas) {
639 pres = 0.5 * (presLiquid + presGas);
640 }
641 }
642 }
643 if (!foundGas || !foundLiquid) {
644 warn_user("MixtureFugacityTP::calculatePsat",
645 "could not find a starting pressure; exiting.");
646 return 0.0;
647 }
648 if (presGas != presLiquid) {
649 warn_user("MixtureFugacityTP::calculatePsat",
650 "could not find a starting pressure; exiting");
651 return 0.0;
652 }
653
654 pres = presGas;
655 double presLast = pres;
656 double RhoGas = RhoGasGood;
657 double RhoLiquid = RhoLiquidGood;
658
659 // Now that we have found a good pressure we can proceed with the algorithm.
660 for (int i = 0; i < 20; i++) {
661 int stab = corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
662 if (stab == 0) {
663 presLast = pres;
664 delGRT = liqGRT - gasGRT;
665 double delV = mw * (1.0/RhoLiquid - 1.0/RhoGas);
666 double dp = - delGRT * GasConstant * TKelvin / delV;
667
668 if (fabs(dp) > 0.1 * pres) {
669 if (dp > 0.0) {
670 dp = 0.1 * pres;
671 } else {
672 dp = -0.1 * pres;
673 }
674 }
675 pres += dp;
676 } else if (stab == -1) {
677 delGRT = 1.0E6;
678 if (presLast > pres) {
679 pres = 0.5 * (presLast + pres);
680 } else {
681 // we are stuck here - try this
682 pres = 1.1 * pres;
683 }
684 } else if (stab == -2) {
685 if (presLast < pres) {
686 pres = 0.5 * (presLast + pres);
687 } else {
688 // we are stuck here - try this
689 pres = 0.9 * pres;
690 }
691 }
692 molarVolGas = mw / RhoGas;
693 molarVolLiquid = mw / RhoLiquid;
694
695 if (fabs(delGRT) < 1.0E-8) {
696 // converged
697 break;
698 }
699 }
700
701 molarVolGas = mw / RhoGas;
702 molarVolLiquid = mw / RhoLiquid;
703 // Put the fluid in the desired end condition
704 setState_TD(tempSave, densSave);
705 return pres;
706 } else {
707 pres = critPressure();
708 setState_TP(TKelvin, pres);
709 molarVolGas = mw / density();
710 molarVolLiquid = molarVolGas;
711 setState_TD(tempSave, densSave);
712 }
713 return pres;
714}
715
716double MixtureFugacityTP::dpdVCalc(double TKelvin, double molarVol, double& presCalc) const
717{
718 throw NotImplementedError("MixtureFugacityTP::dpdVCalc");
719}
720
722{
723 double Tnow = temperature();
724
725 // If the temperature has changed since the last time these
726 // properties were computed, recompute them.
727 if (m_tlast != Tnow) {
728 m_spthermo.update(Tnow, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]);
729 m_tlast = Tnow;
730
731 // update the species Gibbs functions
732 for (size_t k = 0; k < m_kk; k++) {
733 m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
734 }
735 double pref = refPressure();
736 if (pref <= 0.0) {
737 throw CanteraError("MixtureFugacityTP::_updateReferenceStateThermo",
738 "negative reference pressure");
739 }
740 }
741}
742
744{
745 double pc, tc, vc;
746 calcCriticalConditions(pc, tc, vc);
747 return tc;
748}
749
751{
752 double pc, tc, vc;
753 calcCriticalConditions(pc, tc, vc);
754 return pc;
755}
756
758{
759 double pc, tc, vc;
760 calcCriticalConditions(pc, tc, vc);
761 return vc;
762}
763
765{
766 double pc, tc, vc;
767 calcCriticalConditions(pc, tc, vc);
768 return pc*vc/tc/GasConstant;
769}
770
772{
773 double pc, tc, vc;
774 calcCriticalConditions(pc, tc, vc);
775 double mmw = meanMolecularWeight();
776 return mmw / vc;
777}
778
779void MixtureFugacityTP::calcCriticalConditions(double& pc, double& tc, double& vc) const
780{
781 throw NotImplementedError("MixtureFugacityTP::calcCriticalConditions");
782}
783
784int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b,
785 double aAlpha, double Vroot[3], double an,
786 double bn, double cn, double dn, double tc, double vc) const
787{
788 fill_n(Vroot, 3, 0.0);
789 if (T <= 0.0) {
790 throw CanteraError("MixtureFugacityTP::solveCubic",
791 "negative temperature T = {}", T);
792 }
793
794 // Derive the center of the cubic, x_N
795 double xN = - bn /(3 * an);
796
797 // Derive the value of delta**2. This is a key quantity that determines the number of turning points
798 double delta2 = (bn * bn - 3 * an * cn) / (9 * an * an);
799 double delta = 0.0;
800
801 // Calculate a couple of ratios
802 // Cubic equation in z : z^3 - (1-B) z^2 + (A -2B -3B^2)z - (AB- B^2- B^3) = 0
803 double ratio1 = 3.0 * an * cn / (bn * bn);
804 double ratio2 = pres * b / (GasConstant * T); // B
805 if (fabs(ratio1) < 1.0E-7) {
806 double ratio3 = aAlpha / (GasConstant * T) * pres / (GasConstant * T); // A
807 if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) {
808 // A and B terms in cubic equation for z are almost zero, then z is near to 1
809 double zz = 1.0;
810 for (int i = 0; i < 10; i++) {
811 double znew = zz / (zz - ratio2) - ratio3 / (zz + ratio1);
812 double deltaz = znew - zz;
813 zz = znew;
814 if (fabs(deltaz) < 1.0E-14) {
815 break;
816 }
817 }
818 double v = zz * GasConstant * T / pres;
819 Vroot[0] = v;
820 return 1;
821 }
822 }
823
824 int nSolnValues = -1; // Represents number of solutions to the cubic equation
825 double h2 = 4. * an * an * delta2 * delta2 * delta2; // h^2
826 if (delta2 > 0.0) {
827 delta = sqrt(delta2);
828 }
829
830 double h = 2.0 * an * delta * delta2;
831 double yN = 2.0 * bn * bn * bn / (27.0 * an * an) - bn * cn / (3.0 * an) + dn; // y_N term
832 double disc = yN * yN - h2; // discriminant
833
834 //check if y = h
835 if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
836 if (disc > 1e-10) {
837 throw CanteraError("MixtureFugacityTP::solveCubic",
838 "value of yN and h are too high, unrealistic roots may be obtained");
839 }
840 disc = 0.0;
841 }
842
843 if (disc < -1e-14) {
844 // disc<0 then we have three distinct roots.
845 nSolnValues = 3;
846 } else if (fabs(disc) < 1e-14) {
847 // disc=0 then we have two distinct roots (third one is repeated root)
848 nSolnValues = 2;
849 // We are here as p goes to zero.
850 } else if (disc > 1e-14) {
851 // disc> 0 then we have one real root.
852 nSolnValues = 1;
853 }
854
855 double tmp;
856 // One real root -> have to determine whether gas or liquid is the root
857 if (disc > 0.0) {
858 double tmpD = sqrt(disc);
859 double tmp1 = (- yN + tmpD) / (2.0 * an);
860 double sgn1 = 1.0;
861 if (tmp1 < 0.0) {
862 sgn1 = -1.0;
863 tmp1 = -tmp1;
864 }
865 double tmp2 = (- yN - tmpD) / (2.0 * an);
866 double sgn2 = 1.0;
867 if (tmp2 < 0.0) {
868 sgn2 = -1.0;
869 tmp2 = -tmp2;
870 }
871 double p1 = pow(tmp1, 1./3.);
872 double p2 = pow(tmp2, 1./3.);
873 double alpha = xN + sgn1 * p1 + sgn2 * p2;
874 Vroot[0] = alpha;
875 Vroot[1] = 0.0;
876 Vroot[2] = 0.0;
877 } else if (disc < 0.0) {
878 // Three real roots alpha, beta, gamma are obtained.
879 double val = acos(-yN / h);
880 double theta = val / 3.0;
881 double twoThirdPi = 2. * Pi / 3.;
882 double alpha = xN + 2. * delta * cos(theta);
883 double beta = xN + 2. * delta * cos(theta + twoThirdPi);
884 double gamma = xN + 2. * delta * cos(theta + 2.0 * twoThirdPi);
885 Vroot[0] = beta;
886 Vroot[1] = gamma;
887 Vroot[2] = alpha;
888
889 for (int i = 0; i < 3; i++) {
890 tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
891 if (fabs(tmp) > 1.0E-4) {
892 for (int j = 0; j < 3; j++) {
893 if (j != i && fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) {
894 warn_user("MixtureFugacityTP::solveCubic",
895 "roots have merged for T = {}, p = {}: {}, {}",
896 T, pres, Vroot[i], Vroot[j]);
897 }
898 }
899 }
900 }
901 } else if (disc == 0.0) {
902 //Three equal roots are obtained, that is, alpha = beta = gamma
903 if (yN < 1e-18 && h < 1e-18) {
904 // yN = 0.0 and h = 0 (that is, disc = 0)
905 Vroot[0] = xN;
906 Vroot[1] = xN;
907 Vroot[2] = xN;
908 } else {
909 // h and yN need to figure out whether delta^3 is positive or negative
910 if (yN > 0.0) {
911 tmp = pow(yN/(2*an), 1./3.);
912 // In this case, tmp and delta must be equal.
913 if (fabs(tmp - delta) > 1.0E-9) {
914 throw CanteraError("MixtureFugacityTP::solveCubic",
915 "Inconsistency in solver: solver is ill-conditioned.");
916 }
917 Vroot[1] = xN + delta;
918 Vroot[0] = xN - 2.0*delta; // liquid phase root
919 } else {
920 tmp = pow(yN/(2*an), 1./3.);
921 // In this case, tmp and delta must be equal.
922 if (fabs(tmp - delta) > 1.0E-9) {
923 throw CanteraError("MixtureFugacityTP::solveCubic",
924 "Inconsistency in solver: solver is ill-conditioned.");
925 }
926 delta = -delta;
927 Vroot[0] = xN + delta;
928 Vroot[1] = xN - 2.0*delta; // gas phase root
929 }
930 }
931 }
932
933 // Find an accurate root, since there might be a heavy amount of roundoff error due to bad conditioning in this solver.
934 double res, dresdV = 0.0;
935 for (int i = 0; i < nSolnValues; i++) {
936 for (int n = 0; n < 20; n++) {
937 res = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
938 if (fabs(res) < 1.0E-14) { // accurate root is obtained
939 break;
940 }
941 dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn; // derivative of the residual
942 double del = - res / dresdV;
943 Vroot[i] += del;
944 if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) {
945 break;
946 }
947 double res2 = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
948 if (fabs(res2) < fabs(res)) {
949 continue;
950 } else {
951 Vroot[i] -= del; // Go back to previous value of Vroot.
952 Vroot[i] += 0.1 * del; // under-relax by 0.1
953 }
954 }
955 if ((fabs(res) > 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i]))) {
956 throw CanteraError("MixtureFugacityTP::solveCubic",
957 "root failed to converge for T = {}, p = {} with "
958 "V = {}", T, pres, Vroot[i]);
959 }
960 }
961
962 if (nSolnValues == 1) {
963 // Determine the phase of the single root.
964 // nSolnValues = 1 represents the gas phase by default.
965 if (T > tc) {
966 if (Vroot[0] < vc) {
967 // Supercritical phase
968 nSolnValues = -1;
969 }
970 } else {
971 if (Vroot[0] < xN) {
972 //Liquid phase
973 nSolnValues = -1;
974 }
975 }
976 } else {
977 // Determine if we have two distinct roots or three equal roots
978 // nSolnValues = 2 represents 2 equal roots by default.
979 if (nSolnValues == 2 && delta > 1e-14) {
980 //If delta > 0, we have two distinct roots (and one repeated root)
981 nSolnValues = -2;
982 }
983 }
984 return nSolnValues;
985}
986
987}
Header file for a derived class of ThermoPhase that handles non-ideal mixtures based on the fugacity ...
#define FLUID_UNSTABLE
Various states of the Fugacity object.
Base class for exceptions thrown by Cantera classes.
int iState_
Current state of the fluid.
int reportSolnBranchActual() const
Report the solution branch which the solution is actually on.
double enthalpy_mole() const override
Molar enthalpy. Units: J/kmol.
int standardStateConvention() const override
This method returns the convention used in specification of the standard state, of which there are cu...
vector< double > m_g0_RT
Temporary storage for dimensionless reference state Gibbs energies.
double critPressure() const override
Critical pressure (Pa).
double critDensity() const override
Critical density (kg/m3).
void getEntropy_R(double *sr) const override
Get the array of nondimensional Enthalpy functions for the standard state species at the current T an...
vector< double > m_h0_RT
Temporary storage for dimensionless reference state enthalpies.
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...
void getStandardChemPotentials(double *mu) const override
Get the array of chemical potentials at unit activity.
double critTemperature() const override
Critical temperature (K).
void getCp_R(double *cpr) const override
Get the nondimensional Heat Capacities at constant pressure for the standard state of the species at ...
virtual double densSpinodalLiquid() const
Return the value of the density at the liquid spinodal point (on the liquid side) for the current tem...
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
double satPressure(double TKelvin) override
Calculate the saturation pressure at the current mixture content for the given temperature.
void getActivityConcentrations(double *c) const override
This method returns an array of generalized concentrations.
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.
double critCompressibility() const override
Critical compressibility (unitless).
void setPressure(double p) override
Set the internally stored pressure (Pa) at constant temperature and composition.
const vector< double > & gibbs_RT_ref() const
Returns the vector of nondimensional Gibbs free energies of the reference state at the current temper...
void getStandardVolumes_ref(double *vol) const override
Get the molar volumes of the species reference states at the current T and P_ref of the solution.
virtual double densSpinodalGas() const
Return the value of the density at the gas spinodal point (on the gas side) for the current temperatu...
void getPureGibbs(double *gpure) const override
Get the pure Gibbs free energies of each species.
double calculatePsat(double TKelvin, double &molarVolGas, double &molarVolLiquid)
Calculate the saturation pressure at the current mixture content for the given temperature.
vector< double > moleFractions_
Storage for the current values of the mole fractions of the species.
void getEnthalpy_RT(double *hrt) const override
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
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.
vector< double > m_s0_R
Temporary storage for dimensionless reference state entropies.
virtual double dpdVCalc(double TKelvin, double molarVol, double &presCalc) const
Calculate the pressure and the pressure derivative given the temperature and the molar volume.
void getGibbs_RT(double *grt) const override
Get the nondimensional Gibbs functions for the species at their standard states of solution at the cu...
double entropy_mole() const override
Molar entropy. Units: J/kmol/K.
virtual double densityCalc(double TKelvin, double pressure, int phaseRequested, double rhoguess)
Calculates the density given the temperature and the pressure and a guess at the density.
double critVolume() const override
Critical volume (m3/kmol).
virtual double psatEst(double TKelvin) const
Estimate for the saturation pressure.
void getCp_R_ref(double *cprt) const override
Returns the vector of nondimensional constant pressure heat capacities of the reference state at the ...
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...
virtual double sresid() const
Calculate the deviation terms for the total entropy of the mixture from the ideal gas mixture.
void getIntEnergy_RT(double *urt) const override
Returns the vector of nondimensional internal Energies of the standard state at the current temperatu...
int forcedState_
Force the system to be on a particular side of the spinodal curve.
virtual double liquidVolEst(double TKelvin, double &pres) const
Estimate for the molar volume of the liquid.
int forcedSolutionBranch() const
Report the solution branch which the solution is restricted to.
void compositionChanged() override
Apply changes to the state which are needed after the composition changes.
vector< double > m_cp0_R
Temporary storage for dimensionless reference state heat capacities.
int corr0(double TKelvin, double pres, double &densLiq, double &densGas, double &liqGRT, double &gasGRT)
Utility routine in the calculation of the saturation pressure.
void setForcedSolutionBranch(int solnBranch)
Set the solution branch to force the ThermoPhase to exist on one branch or another.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
virtual double hresid() const
Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture.
void getGibbs_RT_ref(double *grt) const override
Returns the vector of nondimensional Gibbs Free Energies of the reference state at the current temper...
double z() const
Calculate the value of z.
int phaseState(bool checkState=false) const
Returns the Phase State flag for the current state of the object.
void getEnthalpy_RT_ref(double *hrt) const override
Returns the vector of nondimensional enthalpies of the reference state at the current temperature of ...
virtual void update(double T, double *cp_R, double *h_RT, double *s_R) const
Compute the reference-state properties for all species.
An error indicating that an unimplemented function has been called.
size_t m_kk
Number of species in the phase.
Definition Phase.h:855
void setState_TD(double t, double rho)
Set the internally stored temperature (K) and density (kg/m^3)
Definition Phase.cpp:377
double temperature() const
Temperature (K).
Definition Phase.h:563
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:656
virtual void setDensity(const double density_)
Set the internally stored density (kg/m^3) of the phase.
Definition Phase.cpp:586
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition Phase.cpp:434
double sum_xlogx() const
Evaluate .
Definition Phase.cpp:626
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:588
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
Definition Phase.cpp:923
virtual void setTemperature(double temp)
Set the internally stored temperature of the phase (K).
Definition Phase.h:624
double mean_X(const double *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition Phase.cpp:616
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:581
virtual void setState_TP(double t, double p)
Set the temperature (K) and pressure (Pa)
double RT() const
Return the Gas Constant multiplied by the current temperature.
double m_tlast
last value of the temperature processed by reference state
virtual void getActivityCoefficients(double *ac) const
Get the array of non-dimensional molar-based activity coefficients at the current solution temperatur...
MultiSpeciesThermo m_spthermo
Pointer to the calculation manager for species reference-state thermodynamic properties.
virtual double refPressure() const
Returns the reference pressure in Pa.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
virtual double gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
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 Pi
Pi.
Definition ct_defs.h:68
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Definition global.h:263
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
const int cSS_CONVENTION_TEMPERATURE
Standard state uses the molar convention.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...