Cantera  3.1.0a4
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 m_tmpV.push_back(0.0);
175 }
176 return added;
177}
178
180{
183 // depends on mole fraction and temperature
184 updateMixingExpressions();
185 iState_ = phaseState(true);
186}
187
189{
190 // A pretty tricky algorithm is needed here, due to problems involving
191 // standard states of real fluids. For those cases you need to combine the T
192 // and P specification for the standard state, or else you may venture into
193 // the forbidden zone, especially when nearing the triple point. Therefore,
194 // we need to do the standard state thermo calc with the (t, pres) combo.
195
196 double t = temperature();
197 double rhoNow = density();
198 if (forcedState_ == FLUID_UNDEFINED) {
199 double rho = densityCalc(t, p, iState_, rhoNow);
200 if (rho > 0.0) {
201 setDensity(rho);
202 iState_ = phaseState(true);
203 } else {
204 if (rho < -1.5) {
205 rho = densityCalc(t, p, FLUID_UNDEFINED , rhoNow);
206 if (rho > 0.0) {
207 setDensity(rho);
208 iState_ = phaseState(true);
209 } else {
210 throw CanteraError("MixtureFugacityTP::setPressure",
211 "neg rho");
212 }
213 } else {
214 throw CanteraError("MixtureFugacityTP::setPressure",
215 "neg rho");
216 }
217 }
218 } else if (forcedState_ == FLUID_GAS) {
219 // Normal density calculation
220 if (iState_ < FLUID_LIQUID_0) {
221 double rho = densityCalc(t, p, iState_, rhoNow);
222 if (rho > 0.0) {
223 setDensity(rho);
224 iState_ = phaseState(true);
225 if (iState_ >= FLUID_LIQUID_0) {
226 throw CanteraError("MixtureFugacityTP::setPressure",
227 "wrong state");
228 }
229 } else {
230 throw CanteraError("MixtureFugacityTP::setPressure",
231 "neg rho");
232 }
233 }
234 } else if (forcedState_ > FLUID_LIQUID_0) {
235 if (iState_ >= FLUID_LIQUID_0) {
236 double rho = densityCalc(t, p, iState_, rhoNow);
237 if (rho > 0.0) {
238 setDensity(rho);
239 iState_ = phaseState(true);
240 if (iState_ == FLUID_GAS) {
241 throw CanteraError("MixtureFugacityTP::setPressure",
242 "wrong state");
243 }
244 } else {
245 throw CanteraError("MixtureFugacityTP::setPressure",
246 "neg rho");
247 }
248 }
249 }
250}
251
253{
256 updateMixingExpressions();
257}
258
260{
262 double p_RT = pressure() / RT();
263 for (size_t k = 0; k < m_kk; k++) {
264 c[k] *= moleFraction(k)*p_RT;
265 }
266}
267
269{
270 return pressure() * meanMolecularWeight() / (density() * RT());
271}
272
274{
275 throw NotImplementedError("MixtureFugacityTP::sresid");
276}
277
279{
280 throw NotImplementedError("MixtureFugacityTP::hresid");
281}
282
283double MixtureFugacityTP::psatEst(double TKelvin) const
284{
285 double pcrit = critPressure();
286 double tt = critTemperature() / TKelvin;
287 if (tt < 1.0) {
288 return pcrit;
289 }
290 double lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918;
291 return pcrit*exp(lpr);
292}
293
294double MixtureFugacityTP::liquidVolEst(double TKelvin, double& pres) const
295{
296 throw NotImplementedError("MixtureFugacityTP::liquidVolEst");
297}
298
299double MixtureFugacityTP::densityCalc(double TKelvin, double presPa,
300 int phase, double rhoguess)
301{
302 double tcrit = critTemperature();
303 double mmw = meanMolecularWeight();
304 if (rhoguess == -1.0) {
305 if (phase != -1) {
306 if (TKelvin > tcrit) {
307 rhoguess = presPa * mmw / (GasConstant * TKelvin);
308 } else {
309 if (phase == FLUID_GAS || phase == FLUID_SUPERCRIT) {
310 rhoguess = presPa * mmw / (GasConstant * TKelvin);
311 } else if (phase >= FLUID_LIQUID_0) {
312 double lqvol = liquidVolEst(TKelvin, presPa);
313 rhoguess = mmw / lqvol;
314 }
315 }
316 } else {
317 // Assume the Gas phase initial guess, if nothing is specified to
318 // the routine
319 rhoguess = presPa * mmw / (GasConstant * TKelvin);
320 }
321 }
322
323 double molarVolBase = mmw / rhoguess;
324 double molarVolLast = molarVolBase;
325 double vc = mmw / critDensity();
326
327 // molar volume of the spinodal at the current temperature and mole
328 // fractions. this will be updated as we go.
329 double molarVolSpinodal = vc;
330 bool conv = false;
331
332 // We start on one side of the vc and stick with that side
333 bool gasSide = molarVolBase > vc;
334 if (gasSide) {
335 molarVolLast = (GasConstant * TKelvin)/presPa;
336 } else {
337 molarVolLast = liquidVolEst(TKelvin, presPa);
338 }
339
340 // OK, now we do a small solve to calculate the molar volume given the T,P
341 // value. The algorithm is taken from dfind()
342 for (int n = 0; n < 200; n++) {
343 // Calculate the predicted reduced pressure, pred0, based on the current
344 // tau and dd. Calculate the derivative of the predicted pressure wrt
345 // the molar volume. This routine also returns the pressure, presBase
346 double presBase;
347 double dpdVBase = dpdVCalc(TKelvin, molarVolBase, presBase);
348
349 // If dpdV is positive, then we are in the middle of the 2 phase region
350 // and beyond the spinodal stability curve. We need to adjust the
351 // initial guess outwards and start a new iteration.
352 if (dpdVBase >= 0.0) {
353 if (TKelvin > tcrit) {
354 throw CanteraError("MixtureFugacityTP::densityCalc",
355 "T > tcrit unexpectedly");
356 }
357
358 // TODO Spawn a calculation for the value of the spinodal point that
359 // is very accurate. Answer the question as to whether a
360 // solution is possible on the current side of the vapor dome.
361 if (gasSide) {
362 if (molarVolBase >= vc) {
363 molarVolSpinodal = molarVolBase;
364 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
365 } else {
366 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
367 }
368 } else {
369 if (molarVolBase <= vc) {
370 molarVolSpinodal = molarVolBase;
371 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
372 } else {
373 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
374 }
375 }
376 continue;
377 }
378
379 // Check for convergence
380 if (fabs(presBase-presPa) < 1.0E-30 + 1.0E-8 * presPa) {
381 conv = true;
382 break;
383 }
384
385 // Dampen and crop the update
386 double dpdV = dpdVBase;
387 if (n < 10) {
388 dpdV = dpdVBase * 1.5;
389 }
390
391 // Formulate the update to the molar volume by Newton's method. Then,
392 // crop it to a max value of 0.1 times the current volume
393 double delMV = - (presBase - presPa) / dpdV;
394 if ((!gasSide || delMV < 0.0) && fabs(delMV) > 0.2 * molarVolBase) {
395 delMV = delMV / fabs(delMV) * 0.2 * molarVolBase;
396 }
397 // Only go 1/10 the way towards the spinodal at any one time.
398 if (TKelvin < tcrit) {
399 if (gasSide) {
400 if (delMV < 0.0 && -delMV > 0.5 * (molarVolBase - molarVolSpinodal)) {
401 delMV = - 0.5 * (molarVolBase - molarVolSpinodal);
402 }
403 } else {
404 if (delMV > 0.0 && delMV > 0.5 * (molarVolSpinodal - molarVolBase)) {
405 delMV = 0.5 * (molarVolSpinodal - molarVolBase);
406 }
407 }
408 }
409 // updated the molar volume value
410 molarVolLast = molarVolBase;
411 molarVolBase += delMV;
412
413 if (fabs(delMV/molarVolBase) < 1.0E-14) {
414 conv = true;
415 break;
416 }
417
418 // Check for negative molar volumes
419 if (molarVolBase <= 0.0) {
420 molarVolBase = std::min(1.0E-30, fabs(delMV*1.0E-4));
421 }
422 }
423
424 // Check for convergence, and return 0.0 if it wasn't achieved.
425 double densBase = 0.0;
426 if (! conv) {
427 molarVolBase = 0.0;
428 throw CanteraError("MixtureFugacityTP::densityCalc",
429 "Process did not converge");
430 } else {
431 densBase = mmw / molarVolBase;
432 }
433 return densBase;
434}
435
436void MixtureFugacityTP::updateMixingExpressions()
437{
438}
439
440int MixtureFugacityTP::corr0(double TKelvin, double pres, double& densLiqGuess,
441 double& densGasGuess, double& liqGRT, double& gasGRT)
442{
443 int retn = 0;
444 double densLiq = densityCalc(TKelvin, pres, FLUID_LIQUID_0, densLiqGuess);
445 if (densLiq <= 0.0) {
446 retn = -1;
447 } else {
448 densLiqGuess = densLiq;
449 setState_TD(TKelvin, densLiq);
450 liqGRT = gibbs_mole() / RT();
451 }
452
453 double densGas = densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess);
454 if (densGas <= 0.0) {
455 if (retn == -1) {
456 throw CanteraError("MixtureFugacityTP::corr0",
457 "Error occurred trying to find gas density at (T,P) = {} {}",
458 TKelvin, pres);
459 }
460 retn = -2;
461 } else {
462 densGasGuess = densGas;
463 setState_TD(TKelvin, densGas);
464 gasGRT = gibbs_mole() / RT();
465 }
466 return retn;
467}
468
469int MixtureFugacityTP::phaseState(bool checkState) const
470{
471 int state = iState_;
472 if (checkState) {
473 double t = temperature();
474 double tcrit = critTemperature();
475 double rhocrit = critDensity();
476 if (t >= tcrit) {
477 return FLUID_SUPERCRIT;
478 }
479 double tmid = tcrit - 100.;
480 if (tmid < 0.0) {
481 tmid = tcrit / 2.0;
482 }
483 double pp = psatEst(tmid);
484 double mmw = meanMolecularWeight();
485 double molVolLiqTmid = liquidVolEst(tmid, pp);
486 double molVolGasTmid = GasConstant * tmid / pp;
487 double densLiqTmid = mmw / molVolLiqTmid;
488 double densGasTmid = mmw / molVolGasTmid;
489 double densMidTmid = 0.5 * (densLiqTmid + densGasTmid);
490 double rhoMid = rhocrit + (t - tcrit) * (rhocrit - densMidTmid) / (tcrit - tmid);
491
492 double rho = density();
493 int iStateGuess = FLUID_LIQUID_0;
494 if (rho < rhoMid) {
495 iStateGuess = FLUID_GAS;
496 }
497 double molarVol = mmw / rho;
498 double presCalc;
499
500 double dpdv = dpdVCalc(t, molarVol, presCalc);
501 if (dpdv < 0.0) {
502 state = iStateGuess;
503 } else {
504 state = FLUID_UNSTABLE;
505 }
506 }
507 return state;
508}
509
511{
512 throw NotImplementedError("MixtureFugacityTP::densSpinodalLiquid");
513}
514
516{
517 throw NotImplementedError("MixtureFugacityTP::densSpinodalGas");
518}
519
520double MixtureFugacityTP::satPressure(double TKelvin)
521{
522 double molarVolGas;
523 double molarVolLiquid;
524 return calculatePsat(TKelvin, molarVolGas, molarVolLiquid);
525}
526
527double MixtureFugacityTP::calculatePsat(double TKelvin, double& molarVolGas,
528 double& molarVolLiquid)
529{
530 // The algorithm for this routine has undergone quite a bit of work. It
531 // probably needs more work. However, it seems now to be fairly robust. The
532 // key requirement is to find an initial pressure where both the liquid and
533 // the gas exist. This is not as easy as it sounds, and it gets exceedingly
534 // hard as the critical temperature is approached from below. Once we have
535 // this initial state, then we seek to equilibrate the Gibbs free energies
536 // of the gas and liquid and use the formula
537 //
538 // dp = VdG
539 //
540 // to create an update condition for deltaP using
541 //
542 // - (Gliq - Ggas) = (Vliq - Vgas) (deltaP)
543 //
544 // @todo Suggestions for the future would be to switch it to an algorithm
545 // that uses the gas molar volume and the liquid molar volumes as the
546 // fundamental unknowns.
547
548 // we need this because this is a non-const routine that is public
549 setTemperature(TKelvin);
550 double densSave = density();
551 double tempSave = temperature();
552 double pres;
553 double mw = meanMolecularWeight();
554 if (TKelvin < critTemperature()) {
555 pres = psatEst(TKelvin);
556 // trial value = Psat from correlation
557 double volLiquid = liquidVolEst(TKelvin, pres);
558 double RhoLiquidGood = mw / volLiquid;
559 double RhoGasGood = pres * mw / (GasConstant * TKelvin);
560 double delGRT = 1.0E6;
561 double liqGRT, gasGRT;
562
563 // First part of the calculation involves finding a pressure at which
564 // the gas and the liquid state coexists.
565 double presLiquid = 0.;
566 double presGas;
567 double presBase = pres;
568 bool foundLiquid = false;
569 bool foundGas = false;
570
571 double densLiquid = densityCalc(TKelvin, presBase, FLUID_LIQUID_0, RhoLiquidGood);
572 if (densLiquid > 0.0) {
573 foundLiquid = true;
574 presLiquid = pres;
575 RhoLiquidGood = densLiquid;
576 }
577 if (!foundLiquid) {
578 for (int i = 0; i < 50; i++) {
579 pres = 1.1 * pres;
580 densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
581 if (densLiquid > 0.0) {
582 foundLiquid = true;
583 presLiquid = pres;
584 RhoLiquidGood = densLiquid;
585 break;
586 }
587 }
588 }
589
590 pres = presBase;
591 double densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
592 if (densGas <= 0.0) {
593 foundGas = false;
594 } else {
595 foundGas = true;
596 presGas = pres;
597 RhoGasGood = densGas;
598 }
599 if (!foundGas) {
600 for (int i = 0; i < 50; i++) {
601 pres = 0.9 * pres;
602 densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
603 if (densGas > 0.0) {
604 foundGas = true;
605 presGas = pres;
606 RhoGasGood = densGas;
607 break;
608 }
609 }
610 }
611
612 if (foundGas && foundLiquid && presGas != presLiquid) {
613 pres = 0.5 * (presLiquid + presGas);
614 bool goodLiq;
615 bool goodGas;
616 for (int i = 0; i < 50; i++) {
617 densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
618 if (densLiquid <= 0.0) {
619 goodLiq = false;
620 } else {
621 goodLiq = true;
622 RhoLiquidGood = densLiquid;
623 presLiquid = pres;
624 }
625 densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
626 if (densGas <= 0.0) {
627 goodGas = false;
628 } else {
629 goodGas = true;
630 RhoGasGood = densGas;
631 presGas = pres;
632 }
633 if (goodGas && goodLiq) {
634 break;
635 }
636 if (!goodLiq && !goodGas) {
637 pres = 0.5 * (pres + presLiquid);
638 }
639 if (goodLiq || goodGas) {
640 pres = 0.5 * (presLiquid + presGas);
641 }
642 }
643 }
644 if (!foundGas || !foundLiquid) {
645 warn_user("MixtureFugacityTP::calculatePsat",
646 "could not find a starting pressure; exiting.");
647 return 0.0;
648 }
649 if (presGas != presLiquid) {
650 warn_user("MixtureFugacityTP::calculatePsat",
651 "could not find a starting pressure; exiting");
652 return 0.0;
653 }
654
655 pres = presGas;
656 double presLast = pres;
657 double RhoGas = RhoGasGood;
658 double RhoLiquid = RhoLiquidGood;
659
660 // Now that we have found a good pressure we can proceed with the algorithm.
661 for (int i = 0; i < 20; i++) {
662 int stab = corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
663 if (stab == 0) {
664 presLast = pres;
665 delGRT = liqGRT - gasGRT;
666 double delV = mw * (1.0/RhoLiquid - 1.0/RhoGas);
667 double dp = - delGRT * GasConstant * TKelvin / delV;
668
669 if (fabs(dp) > 0.1 * pres) {
670 if (dp > 0.0) {
671 dp = 0.1 * pres;
672 } else {
673 dp = -0.1 * pres;
674 }
675 }
676 pres += dp;
677 } else if (stab == -1) {
678 delGRT = 1.0E6;
679 if (presLast > pres) {
680 pres = 0.5 * (presLast + pres);
681 } else {
682 // we are stuck here - try this
683 pres = 1.1 * pres;
684 }
685 } else if (stab == -2) {
686 if (presLast < pres) {
687 pres = 0.5 * (presLast + pres);
688 } else {
689 // we are stuck here - try this
690 pres = 0.9 * pres;
691 }
692 }
693 molarVolGas = mw / RhoGas;
694 molarVolLiquid = mw / RhoLiquid;
695
696 if (fabs(delGRT) < 1.0E-8) {
697 // converged
698 break;
699 }
700 }
701
702 molarVolGas = mw / RhoGas;
703 molarVolLiquid = mw / RhoLiquid;
704 // Put the fluid in the desired end condition
705 setState_TD(tempSave, densSave);
706 return pres;
707 } else {
708 pres = critPressure();
709 setState_TP(TKelvin, pres);
710 molarVolGas = mw / density();
711 molarVolLiquid = molarVolGas;
712 setState_TD(tempSave, densSave);
713 }
714 return pres;
715}
716
717double MixtureFugacityTP::dpdVCalc(double TKelvin, double molarVol, double& presCalc) const
718{
719 throw NotImplementedError("MixtureFugacityTP::dpdVCalc");
720}
721
723{
724 double Tnow = temperature();
725
726 // If the temperature has changed since the last time these
727 // properties were computed, recompute them.
728 if (m_tlast != Tnow) {
729 m_spthermo.update(Tnow, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]);
730 m_tlast = Tnow;
731
732 // update the species Gibbs functions
733 for (size_t k = 0; k < m_kk; k++) {
734 m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
735 }
736 double pref = refPressure();
737 if (pref <= 0.0) {
738 throw CanteraError("MixtureFugacityTP::_updateReferenceStateThermo",
739 "negative reference pressure");
740 }
741 }
742}
743
745{
746 double pc, tc, vc;
747 calcCriticalConditions(pc, tc, vc);
748 return tc;
749}
750
752{
753 double pc, tc, vc;
754 calcCriticalConditions(pc, tc, vc);
755 return pc;
756}
757
759{
760 double pc, tc, vc;
761 calcCriticalConditions(pc, tc, vc);
762 return vc;
763}
764
766{
767 double pc, tc, vc;
768 calcCriticalConditions(pc, tc, vc);
769 return pc*vc/tc/GasConstant;
770}
771
773{
774 double pc, tc, vc;
775 calcCriticalConditions(pc, tc, vc);
776 double mmw = meanMolecularWeight();
777 return mmw / vc;
778}
779
780void MixtureFugacityTP::calcCriticalConditions(double& pc, double& tc, double& vc) const
781{
782 throw NotImplementedError("MixtureFugacityTP::calcCriticalConditions");
783}
784
785int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b,
786 double aAlpha, double Vroot[3], double an,
787 double bn, double cn, double dn, double tc, double vc) const
788{
789 fill_n(Vroot, 3, 0.0);
790 if (T <= 0.0) {
791 throw CanteraError("MixtureFugacityTP::solveCubic",
792 "negative temperature T = {}", T);
793 }
794
795 // Derive the center of the cubic, x_N
796 double xN = - bn /(3 * an);
797
798 // Derive the value of delta**2. This is a key quantity that determines the number of turning points
799 double delta2 = (bn * bn - 3 * an * cn) / (9 * an * an);
800 double delta = 0.0;
801
802 // Calculate a couple of ratios
803 // Cubic equation in z : z^3 - (1-B) z^2 + (A -2B -3B^2)z - (AB- B^2- B^3) = 0
804 double ratio1 = 3.0 * an * cn / (bn * bn);
805 double ratio2 = pres * b / (GasConstant * T); // B
806 if (fabs(ratio1) < 1.0E-7) {
807 double ratio3 = aAlpha / (GasConstant * T) * pres / (GasConstant * T); // A
808 if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) {
809 // A and B terms in cubic equation for z are almost zero, then z is near to 1
810 double zz = 1.0;
811 for (int i = 0; i < 10; i++) {
812 double znew = zz / (zz - ratio2) - ratio3 / (zz + ratio1);
813 double deltaz = znew - zz;
814 zz = znew;
815 if (fabs(deltaz) < 1.0E-14) {
816 break;
817 }
818 }
819 double v = zz * GasConstant * T / pres;
820 Vroot[0] = v;
821 return 1;
822 }
823 }
824
825 int nSolnValues = -1; // Represents number of solutions to the cubic equation
826 double h2 = 4. * an * an * delta2 * delta2 * delta2; // h^2
827 if (delta2 > 0.0) {
828 delta = sqrt(delta2);
829 }
830
831 double h = 2.0 * an * delta * delta2;
832 double yN = 2.0 * bn * bn * bn / (27.0 * an * an) - bn * cn / (3.0 * an) + dn; // y_N term
833 double disc = yN * yN - h2; // discriminant
834
835 //check if y = h
836 if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
837 if (disc > 1e-10) {
838 throw CanteraError("MixtureFugacityTP::solveCubic",
839 "value of yN and h are too high, unrealistic roots may be obtained");
840 }
841 disc = 0.0;
842 }
843
844 if (disc < -1e-14) {
845 // disc<0 then we have three distinct roots.
846 nSolnValues = 3;
847 } else if (fabs(disc) < 1e-14) {
848 // disc=0 then we have two distinct roots (third one is repeated root)
849 nSolnValues = 2;
850 // We are here as p goes to zero.
851 } else if (disc > 1e-14) {
852 // disc> 0 then we have one real root.
853 nSolnValues = 1;
854 }
855
856 double tmp;
857 // One real root -> have to determine whether gas or liquid is the root
858 if (disc > 0.0) {
859 double tmpD = sqrt(disc);
860 double tmp1 = (- yN + tmpD) / (2.0 * an);
861 double sgn1 = 1.0;
862 if (tmp1 < 0.0) {
863 sgn1 = -1.0;
864 tmp1 = -tmp1;
865 }
866 double tmp2 = (- yN - tmpD) / (2.0 * an);
867 double sgn2 = 1.0;
868 if (tmp2 < 0.0) {
869 sgn2 = -1.0;
870 tmp2 = -tmp2;
871 }
872 double p1 = pow(tmp1, 1./3.);
873 double p2 = pow(tmp2, 1./3.);
874 double alpha = xN + sgn1 * p1 + sgn2 * p2;
875 Vroot[0] = alpha;
876 Vroot[1] = 0.0;
877 Vroot[2] = 0.0;
878 } else if (disc < 0.0) {
879 // Three real roots alpha, beta, gamma are obtained.
880 double val = acos(-yN / h);
881 double theta = val / 3.0;
882 double twoThirdPi = 2. * Pi / 3.;
883 double alpha = xN + 2. * delta * cos(theta);
884 double beta = xN + 2. * delta * cos(theta + twoThirdPi);
885 double gamma = xN + 2. * delta * cos(theta + 2.0 * twoThirdPi);
886 Vroot[0] = beta;
887 Vroot[1] = gamma;
888 Vroot[2] = alpha;
889
890 for (int i = 0; i < 3; i++) {
891 tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
892 if (fabs(tmp) > 1.0E-4) {
893 for (int j = 0; j < 3; j++) {
894 if (j != i && fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) {
895 writelog("MixtureFugacityTP::solveCubic(T ={}, p ={}):"
896 " WARNING roots have merged: {}, {}\n",
897 T, pres, Vroot[i], Vroot[j]);
898 }
899 }
900 }
901 }
902 } else if (disc == 0.0) {
903 //Three equal roots are obtained, that is, alpha = beta = gamma
904 if (yN < 1e-18 && h < 1e-18) {
905 // yN = 0.0 and h = 0 (that is, disc = 0)
906 Vroot[0] = xN;
907 Vroot[1] = xN;
908 Vroot[2] = xN;
909 } else {
910 // h and yN need to figure out whether delta^3 is positive or negative
911 if (yN > 0.0) {
912 tmp = pow(yN/(2*an), 1./3.);
913 // In this case, tmp and delta must be equal.
914 if (fabs(tmp - delta) > 1.0E-9) {
915 throw CanteraError("MixtureFugacityTP::solveCubic",
916 "Inconsistency in solver: solver is ill-conditioned.");
917 }
918 Vroot[1] = xN + delta;
919 Vroot[0] = xN - 2.0*delta; // liquid phase root
920 } else {
921 tmp = pow(yN/(2*an), 1./3.);
922 // In this case, tmp and delta must be equal.
923 if (fabs(tmp - delta) > 1.0E-9) {
924 throw CanteraError("MixtureFugacityTP::solveCubic",
925 "Inconsistency in solver: solver is ill-conditioned.");
926 }
927 delta = -delta;
928 Vroot[0] = xN + delta;
929 Vroot[1] = xN - 2.0*delta; // gas phase root
930 }
931 }
932 }
933
934 // Find an accurate root, since there might be a heavy amount of roundoff error due to bad conditioning in this solver.
935 double res, dresdV = 0.0;
936 for (int i = 0; i < nSolnValues; i++) {
937 for (int n = 0; n < 20; n++) {
938 res = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
939 if (fabs(res) < 1.0E-14) { // accurate root is obtained
940 break;
941 }
942 dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn; // derivative of the residual
943 double del = - res / dresdV;
944 Vroot[i] += del;
945 if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) {
946 break;
947 }
948 double res2 = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
949 if (fabs(res2) < fabs(res)) {
950 continue;
951 } else {
952 Vroot[i] -= del; // Go back to previous value of Vroot.
953 Vroot[i] += 0.1 * del; // under-relax by 0.1
954 }
955 }
956 if ((fabs(res) > 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i]))) {
957 writelog("MixtureFugacityTP::solveCubic(T = {}, p = {}): "
958 "WARNING root didn't converge V = {}", T, pres, Vroot[i]);
959 writelogendl();
960 }
961 }
962
963 if (nSolnValues == 1) {
964 // Determine the phase of the single root.
965 // nSolnValues = 1 represents the gas phase by default.
966 if (T > tc) {
967 if (Vroot[0] < vc) {
968 // Supercritical phase
969 nSolnValues = -1;
970 }
971 } else {
972 if (Vroot[0] < xN) {
973 //Liquid phase
974 nSolnValues = -1;
975 }
976 }
977 } else {
978 // Determine if we have two distinct roots or three equal roots
979 // nSolnValues = 2 represents 2 equal roots by default.
980 if (nSolnValues == 2 && delta > 1e-14) {
981 //If delta > 0, we have two distinct roots (and one repeated root)
982 nSolnValues = -2;
983 }
984 }
985 return nSolnValues;
986}
987
988}
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.
vector< double > m_tmpV
Temporary storage - length = m_kk.
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:854
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:562
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:655
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:587
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
Definition Phase.cpp:922
virtual void setTemperature(double temp)
Set the internally stored temperature of the phase (K).
Definition Phase.h:623
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:580
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 writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:175
void writelogendl()
Write an end of line character to the screen and flush output.
Definition global.cpp:41
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:267
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...