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