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++) {
184 updateMixingExpressions();
220 if (
iState_ < FLUID_LIQUID_0) {
225 if (
iState_ >= FLUID_LIQUID_0) {
235 if (
iState_ >= FLUID_LIQUID_0) {
256 updateMixingExpressions();
263 for (
size_t k = 0; k <
m_kk; k++) {
290 double lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918;
291 return pcrit*exp(lpr);
300 int phase,
double rhoguess)
304 if (rhoguess == -1.0) {
306 if (TKelvin > tcrit) {
309 if (phase == FLUID_GAS || phase == FLUID_SUPERCRIT) {
311 }
else if (phase >= FLUID_LIQUID_0) {
313 rhoguess = mmw / lqvol;
323 double molarVolBase = mmw / rhoguess;
324 double molarVolLast = molarVolBase;
329 double molarVolSpinodal = vc;
333 bool gasSide = molarVolBase > vc;
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;
436void MixtureFugacityTP::updateMixingExpressions()
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;
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;
477 return FLUID_SUPERCRIT;
479 double tmid = tcrit - 100.;
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);
493 int iStateGuess = FLUID_LIQUID_0;
495 iStateGuess = FLUID_GAS;
497 double molarVol = mmw / rho;
500 double dpdv =
dpdVCalc(t, molarVol, presCalc);
523 double molarVolLiquid;
528 double& molarVolLiquid)
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;
711 molarVolLiquid = molarVolGas;
733 for (
size_t k = 0; k <
m_kk; k++) {
738 throw CanteraError(
"MixtureFugacityTP::_updateReferenceStateThermo",
739 "negative reference pressure");
747 calcCriticalConditions(pc, tc, vc);
754 calcCriticalConditions(pc, tc, vc);
761 calcCriticalConditions(pc, tc, vc);
768 calcCriticalConditions(pc, tc, vc);
775 calcCriticalConditions(pc, tc, vc);
780void MixtureFugacityTP::calcCriticalConditions(
double& pc,
double& tc,
double& vc)
const
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.
int iState_
Current state of the fluid.
int reportSolnBranchActual() const
Report the solution branch which the solution is actually on.
double enthalpy_mole() const override
Molar enthalpy. Units: J/kmol.
int standardStateConvention() const override
This method returns the convention used in specification of the standard state, of which there are cu...
vector< double > m_g0_RT
Temporary storage for dimensionless reference state Gibbs energies.
double critPressure() const override
Critical pressure (Pa).
double critDensity() const override
Critical density (kg/m3).
void getEntropy_R(double *sr) const override
Get the array of nondimensional Enthalpy functions for the standard state species at the current T an...
vector< double > m_h0_RT
Temporary storage for dimensionless reference state enthalpies.
void getGibbs_ref(double *g) const override
Returns the vector of the Gibbs function of the reference state at the current temperature of the sol...
void getStandardChemPotentials(double *mu) const override
Get the array of chemical potentials at unit activity.
double critTemperature() const override
Critical temperature (K).
void getCp_R(double *cpr) const override
Get the nondimensional Heat Capacities at constant pressure for the standard state of the species at ...
virtual double densSpinodalLiquid() const
Return the value of the density at the liquid spinodal point (on the liquid side) for the current tem...
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
double satPressure(double TKelvin) override
Calculate the saturation pressure at the current mixture content for the given temperature.
vector< double > m_tmpV
Temporary storage - length = m_kk.
void getActivityConcentrations(double *c) const override
This method returns an array of generalized concentrations.
int solveCubic(double T, double pres, double a, double b, double aAlpha, double Vroot[3], double an, double bn, double cn, double dn, double tc, double vc) const
Solve the cubic equation of state.
double critCompressibility() const override
Critical compressibility (unitless).
void setPressure(double p) override
Set the internally stored pressure (Pa) at constant temperature and composition.
const vector< double > & gibbs_RT_ref() const
Returns the vector of nondimensional Gibbs free energies of the reference state at the current temper...
void getStandardVolumes_ref(double *vol) const override
Get the molar volumes of the species reference states at the current T and P_ref of the solution.
virtual double densSpinodalGas() const
Return the value of the density at the gas spinodal point (on the gas side) for the current temperatu...
void getPureGibbs(double *gpure) const override
Get the pure Gibbs free energies of each species.
double calculatePsat(double TKelvin, double &molarVolGas, double &molarVolLiquid)
Calculate the saturation pressure at the current mixture content for the given temperature.
vector< double > moleFractions_
Storage for the current values of the mole fractions of the species.
void getEnthalpy_RT(double *hrt) const override
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
void getEntropy_R_ref(double *er) const override
Returns the vector of nondimensional entropies of the reference state at the current temperature of t...
void setTemperature(const double temp) override
Set the temperature of the phase.
vector< double > m_s0_R
Temporary storage for dimensionless reference state entropies.
virtual double dpdVCalc(double TKelvin, double molarVol, double &presCalc) const
Calculate the pressure and the pressure derivative given the temperature and the molar volume.
void getGibbs_RT(double *grt) const override
Get the nondimensional Gibbs functions for the species at their standard states of solution at the cu...
double entropy_mole() const override
Molar entropy. Units: J/kmol/K.
virtual double densityCalc(double TKelvin, double pressure, int phaseRequested, double rhoguess)
Calculates the density given the temperature and the pressure and a guess at the density.
double critVolume() const override
Critical volume (m3/kmol).
virtual double psatEst(double TKelvin) const
Estimate for the saturation pressure.
void getCp_R_ref(double *cprt) const override
Returns the vector of nondimensional constant pressure heat capacities of the reference state at the ...
void getStandardVolumes(double *vol) const override
Get the molar volumes of each species in their standard states at the current T and P of the solution...
virtual double sresid() const
Calculate the deviation terms for the total entropy of the mixture from the ideal gas mixture.
void getIntEnergy_RT(double *urt) const override
Returns the vector of nondimensional internal Energies of the standard state at the current temperatu...
int forcedState_
Force the system to be on a particular side of the spinodal curve.
virtual double liquidVolEst(double TKelvin, double &pres) const
Estimate for the molar volume of the liquid.
int forcedSolutionBranch() const
Report the solution branch which the solution is restricted to.
void compositionChanged() override
Apply changes to the state which are needed after the composition changes.
vector< double > m_cp0_R
Temporary storage for dimensionless reference state heat capacities.
int corr0(double TKelvin, double pres, double &densLiq, double &densGas, double &liqGRT, double &gasGRT)
Utility routine in the calculation of the saturation pressure.
void setForcedSolutionBranch(int solnBranch)
Set the solution branch to force the ThermoPhase to exist on one branch or another.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
virtual double hresid() const
Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture.
void getGibbs_RT_ref(double *grt) const override
Returns the vector of nondimensional Gibbs Free Energies of the reference state at the current temper...
double z() const
Calculate the value of z.
int phaseState(bool checkState=false) const
Returns the Phase State flag for the current state of the object.
void getEnthalpy_RT_ref(double *hrt) const override
Returns the vector of nondimensional enthalpies of the reference state at the current temperature of ...
virtual void update(double T, double *cp_R, double *h_RT, double *s_R) const
Compute the reference-state properties for all species.
An error indicating that an unimplemented function has been called.
size_t m_kk
Number of species in the phase.
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 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...