45 double h_nonideal =
hresid();
46 return h_ideal + h_nonideal;
54 double s_nonideal =
sresid();
55 return s_ideal + s_nonideal;
64 for (
size_t k = 0; k <
m_kk; k++) {
65 g[k] =
RT() * (g[k] + tmp);
78 for (
size_t k = 0; k <
m_kk; k++) {
87 for (
size_t k = 0; k <
m_kk; k++) {
96 for (
size_t k = 0; k <
m_kk; k++) {
104 for (
size_t i = 0; i <
m_kk; i++) {
116 for (
size_t i = 0; i <
m_kk; i++) {
136 scale(gibbsrt.begin(), gibbsrt.end(), g,
RT());
156 for (
size_t i = 0; i <
m_kk; i++) {
183 updateMixingExpressions();
219 if (
iState_ < FLUID_LIQUID_0) {
224 if (
iState_ >= FLUID_LIQUID_0) {
234 if (
iState_ >= FLUID_LIQUID_0) {
255 updateMixingExpressions();
262 for (
size_t k = 0; k <
m_kk; k++) {
289 double lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918;
290 return pcrit*exp(lpr);
299 int phase,
double rhoguess)
303 if (rhoguess == -1.0) {
305 if (TKelvin > tcrit) {
308 if (phase == FLUID_GAS || phase == FLUID_SUPERCRIT) {
310 }
else if (phase >= FLUID_LIQUID_0) {
312 rhoguess = mmw / lqvol;
322 double molarVolBase = mmw / rhoguess;
323 double molarVolLast = molarVolBase;
328 double molarVolSpinodal = vc;
332 bool gasSide = molarVolBase > vc;
341 for (
int n = 0; n < 200; n++) {
346 double dpdVBase =
dpdVCalc(TKelvin, molarVolBase, presBase);
351 if (dpdVBase >= 0.0) {
352 if (TKelvin > tcrit) {
354 "T > tcrit unexpectedly");
361 if (molarVolBase >= vc) {
362 molarVolSpinodal = molarVolBase;
363 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
365 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
368 if (molarVolBase <= vc) {
369 molarVolSpinodal = molarVolBase;
370 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
372 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
379 if (fabs(presBase-presPa) < 1.0E-30 + 1.0E-8 * presPa) {
385 double dpdV = dpdVBase;
387 dpdV = dpdVBase * 1.5;
392 double delMV = - (presBase - presPa) / dpdV;
393 if ((!gasSide || delMV < 0.0) && fabs(delMV) > 0.2 * molarVolBase) {
394 delMV = delMV / fabs(delMV) * 0.2 * molarVolBase;
397 if (TKelvin < tcrit) {
399 if (delMV < 0.0 && -delMV > 0.5 * (molarVolBase - molarVolSpinodal)) {
400 delMV = - 0.5 * (molarVolBase - molarVolSpinodal);
403 if (delMV > 0.0 && delMV > 0.5 * (molarVolSpinodal - molarVolBase)) {
404 delMV = 0.5 * (molarVolSpinodal - molarVolBase);
409 molarVolLast = molarVolBase;
410 molarVolBase += delMV;
412 if (fabs(delMV/molarVolBase) < 1.0E-14) {
418 if (molarVolBase <= 0.0) {
419 molarVolBase = std::min(1.0E-30, fabs(delMV*1.0E-4));
424 double densBase = 0.0;
428 "Process did not converge");
430 densBase = mmw / molarVolBase;
435void MixtureFugacityTP::updateMixingExpressions()
440 double& densGasGuess,
double& liqGRT,
double& gasGRT)
443 double densLiq =
densityCalc(TKelvin, pres, FLUID_LIQUID_0, densLiqGuess);
444 if (densLiq <= 0.0) {
447 densLiqGuess = densLiq;
452 double densGas =
densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess);
453 if (densGas <= 0.0) {
456 "Error occurred trying to find gas density at (T,P) = {} {}",
461 densGasGuess = densGas;
476 return FLUID_SUPERCRIT;
478 double tmid = tcrit - 100.;
486 double densLiqTmid = mmw / molVolLiqTmid;
487 double densGasTmid = mmw / molVolGasTmid;
488 double densMidTmid = 0.5 * (densLiqTmid + densGasTmid);
489 double rhoMid = rhocrit + (t - tcrit) * (rhocrit - densMidTmid) / (tcrit - tmid);
492 int iStateGuess = FLUID_LIQUID_0;
494 iStateGuess = FLUID_GAS;
496 double molarVol = mmw / rho;
499 double dpdv =
dpdVCalc(t, molarVol, presCalc);
522 double molarVolLiquid;
527 double& molarVolLiquid)
557 double RhoLiquidGood = mw / volLiquid;
558 double RhoGasGood = pres * mw / (
GasConstant * TKelvin);
559 double delGRT = 1.0E6;
560 double liqGRT, gasGRT;
564 double presLiquid = 0.;
566 double presBase = pres;
567 bool foundLiquid =
false;
568 bool foundGas =
false;
570 double densLiquid =
densityCalc(TKelvin, presBase, FLUID_LIQUID_0, RhoLiquidGood);
571 if (densLiquid > 0.0) {
574 RhoLiquidGood = densLiquid;
577 for (
int i = 0; i < 50; i++) {
579 densLiquid =
densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
580 if (densLiquid > 0.0) {
583 RhoLiquidGood = densLiquid;
590 double densGas =
densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
591 if (densGas <= 0.0) {
596 RhoGasGood = densGas;
599 for (
int i = 0; i < 50; i++) {
601 densGas =
densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
605 RhoGasGood = densGas;
611 if (foundGas && foundLiquid && presGas != presLiquid) {
612 pres = 0.5 * (presLiquid + presGas);
615 for (
int i = 0; i < 50; i++) {
616 densLiquid =
densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
617 if (densLiquid <= 0.0) {
621 RhoLiquidGood = densLiquid;
624 densGas =
densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
625 if (densGas <= 0.0) {
629 RhoGasGood = densGas;
632 if (goodGas && goodLiq) {
635 if (!goodLiq && !goodGas) {
636 pres = 0.5 * (pres + presLiquid);
638 if (goodLiq || goodGas) {
639 pres = 0.5 * (presLiquid + presGas);
643 if (!foundGas || !foundLiquid) {
644 warn_user(
"MixtureFugacityTP::calculatePsat",
645 "could not find a starting pressure; exiting.");
648 if (presGas != presLiquid) {
649 warn_user(
"MixtureFugacityTP::calculatePsat",
650 "could not find a starting pressure; exiting");
655 double presLast = pres;
656 double RhoGas = RhoGasGood;
657 double RhoLiquid = RhoLiquidGood;
660 for (
int i = 0; i < 20; i++) {
661 int stab =
corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
664 delGRT = liqGRT - gasGRT;
665 double delV = mw * (1.0/RhoLiquid - 1.0/RhoGas);
666 double dp = - delGRT *
GasConstant * TKelvin / delV;
668 if (fabs(dp) > 0.1 * pres) {
676 }
else if (stab == -1) {
678 if (presLast > pres) {
679 pres = 0.5 * (presLast + pres);
684 }
else if (stab == -2) {
685 if (presLast < pres) {
686 pres = 0.5 * (presLast + pres);
692 molarVolGas = mw / RhoGas;
693 molarVolLiquid = mw / RhoLiquid;
695 if (fabs(delGRT) < 1.0E-8) {
701 molarVolGas = mw / RhoGas;
702 molarVolLiquid = mw / RhoLiquid;
710 molarVolLiquid = molarVolGas;
732 for (
size_t k = 0; k <
m_kk; k++) {
737 throw CanteraError(
"MixtureFugacityTP::_updateReferenceStateThermo",
738 "negative reference pressure");
746 calcCriticalConditions(pc, tc, vc);
753 calcCriticalConditions(pc, tc, vc);
760 calcCriticalConditions(pc, tc, vc);
767 calcCriticalConditions(pc, tc, vc);
774 calcCriticalConditions(pc, tc, vc);
779void MixtureFugacityTP::calcCriticalConditions(
double& pc,
double& tc,
double& vc)
const
785 double aAlpha,
double Vroot[3],
double an,
786 double bn,
double cn,
double dn,
double tc,
double vc)
const
788 fill_n(Vroot, 3, 0.0);
791 "negative temperature T = {}", T);
795 double xN = - bn /(3 * an);
798 double delta2 = (bn * bn - 3 * an * cn) / (9 * an * an);
803 double ratio1 = 3.0 * an * cn / (bn * bn);
805 if (fabs(ratio1) < 1.0E-7) {
807 if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) {
810 for (
int i = 0; i < 10; i++) {
811 double znew = zz / (zz - ratio2) - ratio3 / (zz + ratio1);
812 double deltaz = znew - zz;
814 if (fabs(deltaz) < 1.0E-14) {
824 int nSolnValues = -1;
825 double h2 = 4. * an * an * delta2 * delta2 * delta2;
827 delta = sqrt(delta2);
830 double h = 2.0 * an * delta * delta2;
831 double yN = 2.0 * bn * bn * bn / (27.0 * an * an) - bn * cn / (3.0 * an) + dn;
832 double disc = yN * yN - h2;
835 if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
838 "value of yN and h are too high, unrealistic roots may be obtained");
846 }
else if (fabs(disc) < 1e-14) {
850 }
else if (disc > 1e-14) {
858 double tmpD = sqrt(disc);
859 double tmp1 = (- yN + tmpD) / (2.0 * an);
865 double tmp2 = (- yN - tmpD) / (2.0 * an);
871 double p1 = pow(tmp1, 1./3.);
872 double p2 = pow(tmp2, 1./3.);
873 double alpha = xN + sgn1 * p1 + sgn2 * p2;
877 }
else if (disc < 0.0) {
879 double val = acos(-yN / h);
880 double theta = val / 3.0;
881 double twoThirdPi = 2. *
Pi / 3.;
882 double alpha = xN + 2. * delta * cos(theta);
883 double beta = xN + 2. * delta * cos(theta + twoThirdPi);
884 double gamma = xN + 2. * delta * cos(theta + 2.0 * twoThirdPi);
889 for (
int i = 0; i < 3; i++) {
890 tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
891 if (fabs(tmp) > 1.0E-4) {
892 for (
int j = 0; j < 3; j++) {
893 if (j != i && fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) {
894 warn_user(
"MixtureFugacityTP::solveCubic",
895 "roots have merged for T = {}, p = {}: {}, {}",
896 T, pres, Vroot[i], Vroot[j]);
901 }
else if (disc == 0.0) {
903 if (yN < 1e-18 && h < 1e-18) {
911 tmp = pow(yN/(2*an), 1./3.);
913 if (fabs(tmp - delta) > 1.0E-9) {
915 "Inconsistency in solver: solver is ill-conditioned.");
917 Vroot[1] = xN + delta;
918 Vroot[0] = xN - 2.0*delta;
920 tmp = pow(yN/(2*an), 1./3.);
922 if (fabs(tmp - delta) > 1.0E-9) {
924 "Inconsistency in solver: solver is ill-conditioned.");
927 Vroot[0] = xN + delta;
928 Vroot[1] = xN - 2.0*delta;
934 double res, dresdV = 0.0;
935 for (
int i = 0; i < nSolnValues; i++) {
936 for (
int n = 0; n < 20; n++) {
937 res = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
938 if (fabs(res) < 1.0E-14) {
941 dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn;
942 double del = - res / dresdV;
944 if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) {
947 double res2 = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
948 if (fabs(res2) < fabs(res)) {
952 Vroot[i] += 0.1 * del;
955 if ((fabs(res) > 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i]))) {
957 "root failed to converge for T = {}, p = {} with "
958 "V = {}", T, pres, Vroot[i]);
962 if (nSolnValues == 1) {
979 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.
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.
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...
MultiSpeciesThermo m_spthermo
Pointer to the calculation manager for species reference-state thermodynamic properties.
virtual double refPressure() const
Returns the reference pressure in Pa.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
virtual double gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
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...