45 double h_nonideal =
hresid();
46 return h_ideal + h_nonideal;
54 double s_nonideal =
sresid();
55 return s_ideal + s_nonideal;
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();
76 for (
size_t k = 0; k <
m_kk; k++) {
77 g[k] =
RT() * (g[k] + tmp);
90 for (
size_t k = 0; k <
m_kk; k++) {
99 for (
size_t k = 0; k <
m_kk; k++) {
108 for (
size_t k = 0; k <
m_kk; k++) {
116 for (
size_t i = 0; i <
m_kk; i++) {
128 for (
size_t i = 0; i <
m_kk; i++) {
148 scale(gibbsrt.begin(), gibbsrt.end(), g,
RT());
168 for (
size_t i = 0; i <
m_kk; i++) {
196 updateMixingExpressions();
232 if (
iState_ < FLUID_LIQUID_0) {
237 if (
iState_ >= FLUID_LIQUID_0) {
247 if (
iState_ >= FLUID_LIQUID_0) {
268 updateMixingExpressions();
275 for (
size_t k = 0; k <
m_kk; k++) {
302 double lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918;
303 return pcrit*exp(lpr);
312 int phase,
double rhoguess)
316 if (rhoguess == -1.0) {
318 if (TKelvin > tcrit) {
321 if (phase == FLUID_GAS || phase == FLUID_SUPERCRIT) {
323 }
else if (phase >= FLUID_LIQUID_0) {
325 rhoguess = mmw / lqvol;
335 double molarVolBase = mmw / rhoguess;
336 double molarVolLast = molarVolBase;
341 double molarVolSpinodal = vc;
345 bool gasSide = molarVolBase > vc;
354 for (
int n = 0; n < 200; n++) {
359 double dpdVBase =
dpdVCalc(TKelvin, molarVolBase, presBase);
364 if (dpdVBase >= 0.0) {
365 if (TKelvin > tcrit) {
367 "T > tcrit unexpectedly");
374 if (molarVolBase >= vc) {
375 molarVolSpinodal = molarVolBase;
376 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
378 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
381 if (molarVolBase <= vc) {
382 molarVolSpinodal = molarVolBase;
383 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
385 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
392 if (fabs(presBase-presPa) < 1.0E-30 + 1.0E-8 * presPa) {
398 double dpdV = dpdVBase;
400 dpdV = dpdVBase * 1.5;
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;
410 if (TKelvin < tcrit) {
412 if (delMV < 0.0 && -delMV > 0.5 * (molarVolBase - molarVolSpinodal)) {
413 delMV = - 0.5 * (molarVolBase - molarVolSpinodal);
416 if (delMV > 0.0 && delMV > 0.5 * (molarVolSpinodal - molarVolBase)) {
417 delMV = 0.5 * (molarVolSpinodal - molarVolBase);
422 molarVolLast = molarVolBase;
423 molarVolBase += delMV;
425 if (fabs(delMV/molarVolBase) < 1.0E-14) {
431 if (molarVolBase <= 0.0) {
432 molarVolBase = std::min(1.0E-30, fabs(delMV*1.0E-4));
437 double densBase = 0.0;
441 "Process did not converge");
443 densBase = mmw / molarVolBase;
448void MixtureFugacityTP::updateMixingExpressions()
453 double& densGasGuess,
double& liqGRT,
double& gasGRT)
456 double densLiq =
densityCalc(TKelvin, pres, FLUID_LIQUID_0, densLiqGuess);
457 if (densLiq <= 0.0) {
460 densLiqGuess = densLiq;
465 double densGas =
densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess);
466 if (densGas <= 0.0) {
469 "Error occurred trying to find gas density at (T,P) = {} {}",
474 densGasGuess = densGas;
489 return FLUID_SUPERCRIT;
491 double tmid = tcrit - 100.;
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);
505 int iStateGuess = FLUID_LIQUID_0;
507 iStateGuess = FLUID_GAS;
509 double molarVol = mmw / rho;
512 double dpdv =
dpdVCalc(t, molarVol, presCalc);
535 double molarVolLiquid;
540 double& molarVolLiquid)
570 double RhoLiquidGood = mw / volLiquid;
571 double RhoGasGood = pres * mw / (
GasConstant * TKelvin);
572 double delGRT = 1.0E6;
573 double liqGRT, gasGRT;
577 double presLiquid = 0.;
579 double presBase = pres;
580 bool foundLiquid =
false;
581 bool foundGas =
false;
583 double densLiquid =
densityCalc(TKelvin, presBase, FLUID_LIQUID_0, RhoLiquidGood);
584 if (densLiquid > 0.0) {
587 RhoLiquidGood = densLiquid;
590 for (
int i = 0; i < 50; i++) {
592 densLiquid =
densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
593 if (densLiquid > 0.0) {
596 RhoLiquidGood = densLiquid;
603 double densGas =
densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
604 if (densGas <= 0.0) {
609 RhoGasGood = densGas;
612 for (
int i = 0; i < 50; i++) {
614 densGas =
densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
618 RhoGasGood = densGas;
624 if (foundGas && foundLiquid && presGas != presLiquid) {
625 pres = 0.5 * (presLiquid + presGas);
628 for (
int i = 0; i < 50; i++) {
629 densLiquid =
densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
630 if (densLiquid <= 0.0) {
634 RhoLiquidGood = densLiquid;
637 densGas =
densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
638 if (densGas <= 0.0) {
642 RhoGasGood = densGas;
645 if (goodGas && goodLiq) {
648 if (!goodLiq && !goodGas) {
649 pres = 0.5 * (pres + presLiquid);
651 if (goodLiq || goodGas) {
652 pres = 0.5 * (presLiquid + presGas);
656 if (!foundGas || !foundLiquid) {
657 warn_user(
"MixtureFugacityTP::calculatePsat",
658 "could not find a starting pressure; exiting.");
661 if (presGas != presLiquid) {
662 warn_user(
"MixtureFugacityTP::calculatePsat",
663 "could not find a starting pressure; exiting");
668 double presLast = pres;
669 double RhoGas = RhoGasGood;
670 double RhoLiquid = RhoLiquidGood;
673 for (
int i = 0; i < 20; i++) {
674 int stab =
corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
677 delGRT = liqGRT - gasGRT;
678 double delV = mw * (1.0/RhoLiquid - 1.0/RhoGas);
679 double dp = - delGRT *
GasConstant * TKelvin / delV;
681 if (fabs(dp) > 0.1 * pres) {
689 }
else if (stab == -1) {
691 if (presLast > pres) {
692 pres = 0.5 * (presLast + pres);
697 }
else if (stab == -2) {
698 if (presLast < pres) {
699 pres = 0.5 * (presLast + pres);
705 molarVolGas = mw / RhoGas;
706 molarVolLiquid = mw / RhoLiquid;
708 if (fabs(delGRT) < 1.0E-8) {
714 molarVolGas = mw / RhoGas;
715 molarVolLiquid = mw / RhoLiquid;
723 molarVolLiquid = molarVolGas;
745 for (
size_t k = 0; k <
m_kk; k++) {
750 throw CanteraError(
"MixtureFugacityTP::_updateReferenceStateThermo",
751 "negative reference pressure");
759 calcCriticalConditions(pc, tc, vc);
766 calcCriticalConditions(pc, tc, vc);
773 calcCriticalConditions(pc, tc, vc);
780 calcCriticalConditions(pc, tc, vc);
787 calcCriticalConditions(pc, tc, vc);
792void MixtureFugacityTP::calcCriticalConditions(
double& pc,
double& tc,
double& vc)
const
798 double aAlpha,
double Vroot[3],
double an,
799 double bn,
double cn,
double dn,
double tc,
double vc)
const
801 fill_n(Vroot, 3, 0.0);
804 "negative temperature T = {}", T);
808 double xN = - bn /(3 * an);
811 double delta2 = (bn * bn - 3 * an * cn) / (9 * an * an);
816 double ratio1 = 3.0 * an * cn / (bn * bn);
818 if (fabs(ratio1) < 1.0E-7) {
820 if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) {
823 for (
int i = 0; i < 10; i++) {
824 double znew = zz / (zz - ratio2) - ratio3 / (zz + ratio1);
825 double deltaz = znew - zz;
827 if (fabs(deltaz) < 1.0E-14) {
837 int nSolnValues = -1;
838 double h2 = 4. * an * an * delta2 * delta2 * delta2;
840 delta = sqrt(delta2);
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;
845 double disc = yN * yN - h2;
848 if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
851 "value of yN and h are too high, unrealistic roots may be obtained");
859 }
else if (fabs(disc) < 1e-14) {
863 }
else if (disc > 1e-14) {
871 double tmpD = sqrt(disc);
872 double tmp1 = (- yN + tmpD) / (2.0 * an);
878 double tmp2 = (- yN - tmpD) / (2.0 * an);
884 double p1 = pow(tmp1, 1./3.);
885 double p2 = pow(tmp2, 1./3.);
886 double alpha = xN + sgn1 * p1 + sgn2 * p2;
890 }
else if (disc < 0.0) {
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);
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]);
914 }
else if (disc == 0.0) {
916 if (yN < 1e-18 && h < 1e-18) {
924 tmp = pow(yN/(2*an), 1./3.);
926 if (fabs(tmp - delta) > 1.0E-9) {
928 "Inconsistency in solver: solver is ill-conditioned.");
930 Vroot[1] = xN + delta;
931 Vroot[0] = xN - 2.0*delta;
933 tmp = pow(yN/(2*an), 1./3.);
935 if (fabs(tmp - delta) > 1.0E-9) {
937 "Inconsistency in solver: solver is ill-conditioned.");
940 Vroot[0] = xN + delta;
941 Vroot[1] = xN - 2.0*delta;
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) {
954 dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn;
955 double del = - res / dresdV;
957 if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) {
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)) {
965 Vroot[i] += 0.1 * del;
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]);
975 if (nSolnValues == 1) {
992 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.
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.
void setState_TD(double t, double rho)
Set the internally stored temperature (K) and density (kg/m^3)
double temperature() const
Temperature (K).
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
virtual void setDensity(const double density_)
Set the internally stored density (kg/m^3) of the phase.
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
double sum_xlogx() const
Evaluate .
double moleFraction(size_t k) const
Return the mole fraction of a single species.
virtual double density() const
Density (kg/m^3).
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
virtual void setTemperature(double temp)
Set the internally stored temperature of the phase (K).
double mean_X(const double *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
virtual double pressure() const
Return the thermodynamic pressure (Pa).
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.
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.
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...