21 int MixtureFugacityTP::standardStateConvention()
const
26 void MixtureFugacityTP::setForcedSolutionBranch(
int solnBranch)
28 forcedState_ = solnBranch;
31 int MixtureFugacityTP::forcedSolutionBranch()
const
36 int MixtureFugacityTP::reportSolnBranchActual()
const
42 double MixtureFugacityTP::enthalpy_mole()
const
44 double h_ideal = RT() * mean_X(m_h0_RT);
45 double h_nonideal = hresid();
46 return h_ideal + h_nonideal;
50 double MixtureFugacityTP::entropy_mole()
const
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;
60 void MixtureFugacityTP::getStandardChemPotentials(
double* g)
const
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);
69 void MixtureFugacityTP::getEnthalpy_RT(
double* hrt)
const
71 getEnthalpy_RT_ref(hrt);
74 void MixtureFugacityTP::getEntropy_R(
double* sr)
const
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++) {
83 void MixtureFugacityTP::getGibbs_RT(
double* grt)
const
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++) {
92 void MixtureFugacityTP::getPureGibbs(
double* g)
const
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++) {
101 void MixtureFugacityTP::getIntEnergy_RT(
double* urt)
const
103 copy(m_h0_RT.begin(), m_h0_RT.end(), urt);
104 for (
size_t i = 0; i < m_kk; i++) {
109 void MixtureFugacityTP::getCp_R(
double* cpr)
const
111 copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
114 void MixtureFugacityTP::getStandardVolumes(
double* vol)
const
116 for (
size_t i = 0; i < m_kk; i++) {
117 vol[i] = RT() / pressure();
123 void MixtureFugacityTP::getEnthalpy_RT_ref(
double* hrt)
const
125 copy(m_h0_RT.begin(), m_h0_RT.end(), hrt);
128 void MixtureFugacityTP::getGibbs_RT_ref(
double* grt)
const
130 copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
133 void MixtureFugacityTP::getGibbs_ref(
double* g)
const
135 const vector<double>& gibbsrt = gibbs_RT_ref();
136 scale(gibbsrt.begin(), gibbsrt.end(), g, RT());
139 const vector<double>& MixtureFugacityTP::gibbs_RT_ref()
const
144 void MixtureFugacityTP::getEntropy_R_ref(
double* er)
const
146 copy(m_s0_R.begin(), m_s0_R.end(), er);
149 void MixtureFugacityTP::getCp_R_ref(
double* cpr)
const
151 copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
154 void MixtureFugacityTP::getStandardVolumes_ref(
double* vol)
const
156 for (
size_t i = 0; i < m_kk; i++) {
157 vol[i]= RT() / refPressure();
161 bool MixtureFugacityTP::addSpecies(shared_ptr<Species> spec)
163 bool added = ThermoPhase::addSpecies(spec);
166 moleFractions_.push_back(1.0);
168 moleFractions_.push_back(0.0);
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);
179 void MixtureFugacityTP::setTemperature(
const double temp)
181 Phase::setTemperature(temp);
182 _updateReferenceStateThermo();
184 updateMixingExpressions();
185 iState_ = phaseState(
true);
188 void MixtureFugacityTP::setPressure(
double p)
196 double t = temperature();
197 double rhoNow = density();
198 if (forcedState_ == FLUID_UNDEFINED) {
199 double rho = densityCalc(t, p, iState_, rhoNow);
202 iState_ = phaseState(
true);
205 rho = densityCalc(t, p, FLUID_UNDEFINED , rhoNow);
208 iState_ = phaseState(
true);
218 }
else if (forcedState_ == FLUID_GAS) {
220 if (iState_ < FLUID_LIQUID_0) {
221 double rho = densityCalc(t, p, iState_, rhoNow);
224 iState_ = phaseState(
true);
225 if (iState_ >= FLUID_LIQUID_0) {
234 }
else if (forcedState_ > FLUID_LIQUID_0) {
235 if (iState_ >= FLUID_LIQUID_0) {
236 double rho = densityCalc(t, p, iState_, rhoNow);
239 iState_ = phaseState(
true);
240 if (iState_ == FLUID_GAS) {
252 void MixtureFugacityTP::compositionChanged()
254 Phase::compositionChanged();
255 getMoleFractions(moleFractions_.data());
256 updateMixingExpressions();
259 void MixtureFugacityTP::getActivityConcentrations(
double* c)
const
261 getActivityCoefficients(c);
262 double p_RT = pressure() / RT();
263 for (
size_t k = 0; k < m_kk; k++) {
264 c[k] *= moleFraction(k)*p_RT;
268 double MixtureFugacityTP::z()
const
270 return pressure() * meanMolecularWeight() / (density() * RT());
273 double MixtureFugacityTP::sresid()
const
278 double MixtureFugacityTP::hresid()
const
283 double MixtureFugacityTP::psatEst(
double TKelvin)
const
285 double pcrit = critPressure();
286 double tt = critTemperature() / TKelvin;
290 double lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918;
291 return pcrit*exp(lpr);
294 double MixtureFugacityTP::liquidVolEst(
double TKelvin,
double& pres)
const
299 double MixtureFugacityTP::densityCalc(
double TKelvin,
double presPa,
300 int phase,
double rhoguess)
302 double tcrit = critTemperature();
303 double mmw = meanMolecularWeight();
304 if (rhoguess == -1.0) {
306 if (TKelvin > tcrit) {
309 if (phase == FLUID_GAS || phase == FLUID_SUPERCRIT) {
311 }
else if (phase >= FLUID_LIQUID_0) {
312 double lqvol = liquidVolEst(TKelvin, presPa);
313 rhoguess = mmw / lqvol;
323 double molarVolBase = mmw / rhoguess;
324 double molarVolLast = molarVolBase;
325 double vc = mmw / critDensity();
329 double molarVolSpinodal = vc;
333 bool gasSide = molarVolBase > vc;
337 molarVolLast = liquidVolEst(TKelvin, presPa);
342 for (
int n = 0; n < 200; n++) {
347 double dpdVBase = dpdVCalc(TKelvin, molarVolBase, presBase);
352 if (dpdVBase >= 0.0) {
353 if (TKelvin > tcrit) {
355 "T > tcrit unexpectedly");
362 if (molarVolBase >= vc) {
363 molarVolSpinodal = molarVolBase;
364 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
366 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
369 if (molarVolBase <= vc) {
370 molarVolSpinodal = molarVolBase;
371 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
373 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
380 if (fabs(presBase-presPa) < 1.0E-30 + 1.0E-8 * presPa) {
386 double dpdV = dpdVBase;
388 dpdV = dpdVBase * 1.5;
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;
398 if (TKelvin < tcrit) {
400 if (delMV < 0.0 && -delMV > 0.5 * (molarVolBase - molarVolSpinodal)) {
401 delMV = - 0.5 * (molarVolBase - molarVolSpinodal);
404 if (delMV > 0.0 && delMV > 0.5 * (molarVolSpinodal - molarVolBase)) {
405 delMV = 0.5 * (molarVolSpinodal - molarVolBase);
410 molarVolLast = molarVolBase;
411 molarVolBase += delMV;
413 if (fabs(delMV/molarVolBase) < 1.0E-14) {
419 if (molarVolBase <= 0.0) {
420 molarVolBase = std::min(1.0E-30, fabs(delMV*1.0E-4));
425 double densBase = 0.0;
429 "Process did not converge");
431 densBase = mmw / molarVolBase;
436 void MixtureFugacityTP::updateMixingExpressions()
440 int MixtureFugacityTP::corr0(
double TKelvin,
double pres,
double& densLiqGuess,
441 double& densGasGuess,
double& liqGRT,
double& gasGRT)
444 double densLiq = densityCalc(TKelvin, pres, FLUID_LIQUID_0, densLiqGuess);
445 if (densLiq <= 0.0) {
448 densLiqGuess = densLiq;
449 setState_TD(TKelvin, densLiq);
450 liqGRT = gibbs_mole() / RT();
453 double densGas = densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess);
454 if (densGas <= 0.0) {
457 "Error occurred trying to find gas density at (T,P) = {} {}",
462 densGasGuess = densGas;
463 setState_TD(TKelvin, densGas);
464 gasGRT = gibbs_mole() / RT();
469 int MixtureFugacityTP::phaseState(
bool checkState)
const
473 double t = temperature();
474 double tcrit = critTemperature();
475 double rhocrit = critDensity();
477 return FLUID_SUPERCRIT;
479 double tmid = tcrit - 100.;
483 double pp = psatEst(tmid);
484 double mmw = meanMolecularWeight();
485 double molVolLiqTmid = liquidVolEst(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);
492 double rho = density();
493 int iStateGuess = FLUID_LIQUID_0;
495 iStateGuess = FLUID_GAS;
497 double molarVol = mmw / rho;
500 double dpdv = dpdVCalc(t, molarVol, presCalc);
510 double MixtureFugacityTP::densSpinodalLiquid()
const
515 double MixtureFugacityTP::densSpinodalGas()
const
520 double MixtureFugacityTP::satPressure(
double TKelvin)
523 double molarVolLiquid;
524 return calculatePsat(TKelvin, molarVolGas, molarVolLiquid);
527 double MixtureFugacityTP::calculatePsat(
double TKelvin,
double& molarVolGas,
528 double& molarVolLiquid)
549 setTemperature(TKelvin);
550 double densSave = density();
551 double tempSave = temperature();
553 double mw = meanMolecularWeight();
554 if (TKelvin < critTemperature()) {
555 pres = psatEst(TKelvin);
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;
565 double presLiquid = 0.;
567 double presBase = pres;
568 bool foundLiquid =
false;
569 bool foundGas =
false;
571 double densLiquid = densityCalc(TKelvin, presBase, FLUID_LIQUID_0, RhoLiquidGood);
572 if (densLiquid > 0.0) {
575 RhoLiquidGood = densLiquid;
578 for (
int i = 0; i < 50; i++) {
580 densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
581 if (densLiquid > 0.0) {
584 RhoLiquidGood = densLiquid;
591 double densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
592 if (densGas <= 0.0) {
597 RhoGasGood = densGas;
600 for (
int i = 0; i < 50; i++) {
602 densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
606 RhoGasGood = densGas;
612 if (foundGas && foundLiquid && presGas != presLiquid) {
613 pres = 0.5 * (presLiquid + presGas);
616 for (
int i = 0; i < 50; i++) {
617 densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
618 if (densLiquid <= 0.0) {
622 RhoLiquidGood = densLiquid;
625 densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
626 if (densGas <= 0.0) {
630 RhoGasGood = densGas;
633 if (goodGas && goodLiq) {
636 if (!goodLiq && !goodGas) {
637 pres = 0.5 * (pres + presLiquid);
639 if (goodLiq || goodGas) {
640 pres = 0.5 * (presLiquid + presGas);
644 if (!foundGas || !foundLiquid) {
645 warn_user(
"MixtureFugacityTP::calculatePsat",
646 "could not find a starting pressure; exiting.");
649 if (presGas != presLiquid) {
650 warn_user(
"MixtureFugacityTP::calculatePsat",
651 "could not find a starting pressure; exiting");
656 double presLast = pres;
657 double RhoGas = RhoGasGood;
658 double RhoLiquid = RhoLiquidGood;
661 for (
int i = 0; i < 20; i++) {
662 int stab = corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
665 delGRT = liqGRT - gasGRT;
666 double delV = mw * (1.0/RhoLiquid - 1.0/RhoGas);
667 double dp = - delGRT *
GasConstant * TKelvin / delV;
669 if (fabs(dp) > 0.1 * pres) {
677 }
else if (stab == -1) {
679 if (presLast > pres) {
680 pres = 0.5 * (presLast + pres);
685 }
else if (stab == -2) {
686 if (presLast < pres) {
687 pres = 0.5 * (presLast + pres);
693 molarVolGas = mw / RhoGas;
694 molarVolLiquid = mw / RhoLiquid;
696 if (fabs(delGRT) < 1.0E-8) {
702 molarVolGas = mw / RhoGas;
703 molarVolLiquid = mw / RhoLiquid;
705 setState_TD(tempSave, densSave);
708 pres = critPressure();
709 setState_TP(TKelvin, pres);
710 molarVolGas = mw / density();
711 molarVolLiquid = molarVolGas;
712 setState_TD(tempSave, densSave);
717 double MixtureFugacityTP::dpdVCalc(
double TKelvin,
double molarVol,
double& presCalc)
const
722 void MixtureFugacityTP::_updateReferenceStateThermo()
const
724 double Tnow = temperature();
728 if (m_tlast != Tnow) {
729 m_spthermo.update(Tnow, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]);
733 for (
size_t k = 0; k < m_kk; k++) {
734 m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
736 double pref = refPressure();
738 throw CanteraError(
"MixtureFugacityTP::_updateReferenceStateThermo",
739 "negative reference pressure");
744 double MixtureFugacityTP::critTemperature()
const
747 calcCriticalConditions(pc, tc, vc);
751 double MixtureFugacityTP::critPressure()
const
754 calcCriticalConditions(pc, tc, vc);
758 double MixtureFugacityTP::critVolume()
const
761 calcCriticalConditions(pc, tc, vc);
765 double MixtureFugacityTP::critCompressibility()
const
768 calcCriticalConditions(pc, tc, vc);
772 double MixtureFugacityTP::critDensity()
const
775 calcCriticalConditions(pc, tc, vc);
776 double mmw = meanMolecularWeight();
780 void MixtureFugacityTP::calcCriticalConditions(
double& pc,
double& tc,
double& vc)
const
785 int 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
789 fill_n(Vroot, 3, 0.0);
792 "negative temperature T = {}", T);
796 double xN = - bn /(3 * an);
799 double delta2 = (bn * bn - 3 * an * cn) / (9 * an * an);
804 double ratio1 = 3.0 * an * cn / (bn * bn);
806 if (fabs(ratio1) < 1.0E-7) {
808 if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) {
811 for (
int i = 0; i < 10; i++) {
812 double znew = zz / (zz - ratio2) - ratio3 / (zz + ratio1);
813 double deltaz = znew - zz;
815 if (fabs(deltaz) < 1.0E-14) {
825 int nSolnValues = -1;
826 double h2 = 4. * an * an * delta2 * delta2 * delta2;
828 delta = sqrt(delta2);
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;
833 double disc = yN * yN - h2;
836 if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
839 "value of yN and h are too high, unrealistic roots may be obtained");
847 }
else if (fabs(disc) < 1e-14) {
851 }
else if (disc > 1e-14) {
859 double tmpD = sqrt(disc);
860 double tmp1 = (- yN + tmpD) / (2.0 * an);
866 double tmp2 = (- yN - tmpD) / (2.0 * an);
872 double p1 = pow(tmp1, 1./3.);
873 double p2 = pow(tmp2, 1./3.);
874 double alpha = xN + sgn1 * p1 + sgn2 * p2;
878 }
else if (disc < 0.0) {
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);
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]);
902 }
else if (disc == 0.0) {
904 if (yN < 1e-18 && h < 1e-18) {
912 tmp = pow(yN/(2*an), 1./3.);
914 if (fabs(tmp - delta) > 1.0E-9) {
916 "Inconsistency in solver: solver is ill-conditioned.");
918 Vroot[1] = xN + delta;
919 Vroot[0] = xN - 2.0*delta;
921 tmp = pow(yN/(2*an), 1./3.);
923 if (fabs(tmp - delta) > 1.0E-9) {
925 "Inconsistency in solver: solver is ill-conditioned.");
928 Vroot[0] = xN + delta;
929 Vroot[1] = xN - 2.0*delta;
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) {
942 dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn;
943 double del = - res / dresdV;
945 if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) {
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)) {
953 Vroot[i] += 0.1 * del;
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]);
963 if (nSolnValues == 1) {
980 if (nSolnValues == 2 && delta > 1e-14) {
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.
An error indicating that an unimplemented function has been called.
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.
void writelogendl()
Write an end of line character to the screen and flush output.
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
const double GasConstant
Universal Gas Constant [J/kmol/K].
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Namespace for the Cantera kernel.
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...