20 MixtureFugacityTP::MixtureFugacityTP() :
22 forcedState_(FLUID_UNDEFINED),
52 for (
size_t k = 0; k <
m_kk; k++) {
53 muRT[k] *= 1.0 /
RT();
64 for (
size_t k = 0; k <
m_kk; k++) {
65 g[k] =
RT() * (g[k] + tmp);
79 for (
size_t k = 0; k <
m_kk; k++) {
89 for (
size_t k = 0; k <
m_kk; k++) {
99 for (
size_t k = 0; k <
m_kk; k++) {
108 for (
size_t i = 0; i <
m_kk; i++) {
122 for (
size_t i = 0; i <
m_kk; i++) {
144 scale(gibbsrt.begin(), gibbsrt.end(), g,
RT());
168 for (
size_t i = 0; i <
m_kk; i++) {
190 if (state.
hasChild(
"temperature")) {
191 t =
getFloat(state,
"temperature",
"temperature");
195 double p =
getFloat(state,
"pressure",
"pressure");
197 }
else if (state.
hasChild(
"density")) {
198 double rho =
getFloat(state,
"density",
"density");
240 void MixtureFugacityTP::setMoleFractions_NoState(
const doublereal*
const x)
244 updateMixingExpressions();
250 "called, but EOS for phase is not known");
265 updateMixingExpressions();
275 rho =
densityCalc(t, pres, FLUID_UNDEFINED , rhoNow);
280 throw CanteraError(
"MixtureFugacityTP::setState_TP()",
"neg rho");
283 throw CanteraError(
"MixtureFugacityTP::setState_TP()",
"neg rho");
288 if (
iState_ < FLUID_LIQUID_0) {
294 if (
iState_ >= FLUID_LIQUID_0) {
295 throw CanteraError(
"MixtureFugacityTP::setState_TP()",
"wrong state");
298 throw CanteraError(
"MixtureFugacityTP::setState_TP()",
"neg rho");
302 if (
iState_ >= FLUID_LIQUID_0) {
309 throw CanteraError(
"MixtureFugacityTP::setState_TP()",
"wrong state");
312 throw CanteraError(
"MixtureFugacityTP::setState_TP()",
"neg rho");
318 void MixtureFugacityTP::setState_TR(doublereal T, doublereal rho)
325 updateMixingExpressions();
331 setMoleFractions_NoState(x);
342 throw CanteraError(
"MixtureFugacityTP::sresid()",
"Base Class: not implemented");
347 throw CanteraError(
"MixtureFugacityTP::hresid()",
"Base Class: not implemented");
357 doublereal lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918;
358 return pcrit*exp(lpr);
363 throw CanteraError(
"MixtureFugacityTP::liquidVolEst()",
"unimplemented");
367 int phase, doublereal rhoguess)
371 if (rhoguess == -1.0) {
373 if (TKelvin > tcrit) {
376 if (phase == FLUID_GAS || phase == FLUID_SUPERCRIT) {
378 }
else if (phase >= FLUID_LIQUID_0) {
380 rhoguess = mmw / lqvol;
390 double molarVolBase = mmw / rhoguess;
391 double molarVolLast = molarVolBase;
396 double molarVolSpinodal = vc;
400 bool gasSide = molarVolBase > vc;
409 for (
int n = 0; n < 200; n++) {
414 double dpdVBase =
dpdVCalc(TKelvin, molarVolBase, presBase);
419 if (dpdVBase >= 0.0) {
420 if (TKelvin > tcrit) {
422 "T > tcrit unexpectedly");
429 if (molarVolBase >= vc) {
430 molarVolSpinodal = molarVolBase;
431 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
433 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
436 if (molarVolBase <= vc) {
437 molarVolSpinodal = molarVolBase;
438 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
440 molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
447 if (fabs(presBase-presPa) < 1.0E-30 + 1.0E-8 * presPa) {
453 doublereal dpdV = dpdVBase;
455 dpdV = dpdVBase * 1.5;
460 double delMV = - (presBase - presPa) / dpdV;
461 if ((!gasSide || delMV < 0.0) && fabs(delMV) > 0.2 * molarVolBase) {
462 delMV = delMV / fabs(delMV) * 0.2 * molarVolBase;
465 if (TKelvin < tcrit) {
467 if (delMV < 0.0 && -delMV > 0.5 * (molarVolBase - molarVolSpinodal)) {
468 delMV = - 0.5 * (molarVolBase - molarVolSpinodal);
471 if (delMV > 0.0 && delMV > 0.5 * (molarVolSpinodal - molarVolBase)) {
472 delMV = 0.5 * (molarVolSpinodal - molarVolBase);
477 molarVolLast = molarVolBase;
478 molarVolBase += delMV;
480 if (fabs(delMV/molarVolBase) < 1.0E-14) {
486 if (molarVolBase <= 0.0) {
487 molarVolBase = std::min(1.0E-30, fabs(delMV*1.0E-4));
492 double densBase = 0.0;
495 throw CanteraError(
"MixtureFugacityTP::densityCalc()",
"Process did not converge");
497 densBase = mmw / molarVolBase;
502 void MixtureFugacityTP::updateMixingExpressions()
507 doublereal& densGasGuess, doublereal& liqGRT, doublereal& gasGRT)
510 doublereal densLiq =
densityCalc(TKelvin, pres, FLUID_LIQUID_0, densLiqGuess);
511 if (densLiq <= 0.0) {
514 densLiqGuess = densLiq;
515 setState_TR(TKelvin, densLiq);
519 doublereal densGas =
densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess);
520 if (densGas <= 0.0) {
523 "Error occurred trying to find gas density at (T,P) = {} {}",
528 densGasGuess = densGas;
529 setState_TR(TKelvin, densGas);
543 return FLUID_SUPERCRIT;
545 double tmid = tcrit - 100.;
553 double densLiqTmid = mmw / molVolLiqTmid;
554 double densGasTmid = mmw / molVolGasTmid;
555 double densMidTmid = 0.5 * (densLiqTmid + densGasTmid);
556 doublereal rhoMid = rhocrit + (t - tcrit) * (rhocrit - densMidTmid) / (tcrit - tmid);
559 int iStateGuess = FLUID_LIQUID_0;
561 iStateGuess = FLUID_GAS;
563 double molarVol = mmw / rho;
566 double dpdv =
dpdVCalc(t, molarVol, presCalc);
578 throw CanteraError(
"MixtureFugacityTP::densSpinodalLiquid",
"unimplemented");
583 throw CanteraError(
"MixtureFugacityTP::densSpinodalGas",
"unimplemented");
588 doublereal molarVolGas;
589 doublereal molarVolLiquid;
594 doublereal& molarVolLiquid)
624 double RhoLiquidGood = mw / volLiquid;
625 double RhoGasGood = pres * mw / (
GasConstant * TKelvin);
626 doublereal delGRT = 1.0E6;
627 doublereal liqGRT, gasGRT;
631 doublereal presLiquid = 0.;
633 doublereal presBase = pres;
634 bool foundLiquid =
false;
635 bool foundGas =
false;
637 doublereal densLiquid =
densityCalc(TKelvin, presBase, FLUID_LIQUID_0, RhoLiquidGood);
638 if (densLiquid > 0.0) {
641 RhoLiquidGood = densLiquid;
644 for (
int i = 0; i < 50; i++) {
646 densLiquid =
densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
647 if (densLiquid > 0.0) {
650 RhoLiquidGood = densLiquid;
657 doublereal densGas =
densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
658 if (densGas <= 0.0) {
663 RhoGasGood = densGas;
666 for (
int i = 0; i < 50; i++) {
668 densGas =
densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
672 RhoGasGood = densGas;
678 if (foundGas && foundLiquid && presGas != presLiquid) {
679 pres = 0.5 * (presLiquid + presGas);
682 for (
int i = 0; i < 50; i++) {
683 densLiquid =
densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
684 if (densLiquid <= 0.0) {
688 RhoLiquidGood = densLiquid;
691 densGas =
densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
692 if (densGas <= 0.0) {
696 RhoGasGood = densGas;
699 if (goodGas && goodLiq) {
702 if (!goodLiq && !goodGas) {
703 pres = 0.5 * (pres + presLiquid);
705 if (goodLiq || goodGas) {
706 pres = 0.5 * (presLiquid + presGas);
710 if (!foundGas || !foundLiquid) {
711 writelog(
"error couldn't find a starting pressure\n");
714 if (presGas != presLiquid) {
715 writelog(
"error couldn't find a starting pressure\n");
720 double presLast = pres;
721 double RhoGas = RhoGasGood;
722 double RhoLiquid = RhoLiquidGood;
725 for (
int i = 0; i < 20; i++) {
726 int stab =
corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
729 delGRT = liqGRT - gasGRT;
730 doublereal delV = mw * (1.0/RhoLiquid - 1.0/RhoGas);
731 doublereal dp = - delGRT *
GasConstant * TKelvin / delV;
733 if (fabs(dp) > 0.1 * pres) {
741 }
else if (stab == -1) {
743 if (presLast > pres) {
744 pres = 0.5 * (presLast + pres);
749 }
else if (stab == -2) {
750 if (presLast < pres) {
751 pres = 0.5 * (presLast + pres);
757 molarVolGas = mw / RhoGas;
758 molarVolLiquid = mw / RhoLiquid;
760 if (fabs(delGRT) < 1.0E-8) {
766 molarVolGas = mw / RhoGas;
767 molarVolLiquid = mw / RhoLiquid;
769 setState_TR(tempSave, densSave);
775 molarVolLiquid = molarVolGas;
776 setState_TR(tempSave, densSave);
783 throw CanteraError(
"MixtureFugacityTP::pressureCalc",
"unimplemented");
788 throw CanteraError(
"MixtureFugacityTP::dpdVCalc",
"unimplemented");
802 for (
size_t k = 0; k <
m_kk; k++) {
807 throw CanteraError(
"MixtureFugacityTP::_updateReferenceStateThermo()",
"neg ref pressure");
virtual int forcedSolutionBranch() const
Report the solution branch which the solution is restricted to.
virtual doublereal satPressure(doublereal TKelvin)
Calculate the saturation pressure at the current mixture content for the given temperature.
virtual void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input...
virtual void setState_TPX(doublereal t, doublereal p, const doublereal *x)
Set the temperature (K), pressure (Pa), and mole fractions.
virtual void getEnthalpy_RT_ref(doublereal *hrt) const
virtual bool addSpecies(shared_ptr< Species > spec)
virtual doublereal sresid() const
Calculate the deviation terms for the total entropy of the mixture from the ideal gas mixture...
virtual doublereal densityCalc(doublereal TKelvin, doublereal pressure, int phaseRequested, doublereal rhoguess)
Calculates the density given the temperature and the pressure and a guess at the density.
virtual int standardStateConvention() const
This method returns the convention used in specification of the standard state, of which there are cu...
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
MultiSpeciesThermo m_spthermo
Pointer to the calculation manager for species reference-state thermodynamic properties.
virtual void getGibbs_RT_ref(doublereal *grt) const
Returns the vector of nondimensional Gibbs Free Energies of the reference state at the current temper...
std::string getChildValue(const XML_Node &parent, const std::string &nameString)
This function reads a child node with the name, nameString, and returns its XML value as the return s...
doublereal temperature() const
Temperature (K).
virtual doublereal hresid() const
Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture...
An error indicating that an unimplemented function has been called.
virtual void getChemPotentials_RT(doublereal *mu) const
Get the array of non-dimensional species chemical potentials These are partial molar Gibbs free energ...
virtual doublereal psatEst(doublereal TKelvin) const
Estimate for the saturation pressure.
virtual void update(doublereal T, doublereal *cp_R, doublereal *h_RT, doublereal *s_R) const
Compute the reference-state properties for all species.
virtual doublereal dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal &presCalc) const
Calculate the pressure and the pressure derivative given the temperature and the molar volume...
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
doublereal z() const
Calculate the value of z.
int forcedState_
Force the system to be on a particular side of the spinodal curve.
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
const int cSS_CONVENTION_TEMPERATURE
Standard state uses the molar convention.
vector_fp m_g0_RT
Temporary storage for dimensionless reference state Gibbs energies.
Class XML_Node is a tree-based representation of the contents of an XML file.
virtual doublereal density() const
Density (kg/m^3).
doublereal calculatePsat(doublereal TKelvin, doublereal &molarVolGas, doublereal &molarVolLiquid)
Calculate the saturation pressure at the current mixture content for the given temperature.
#define FLUID_UNSTABLE
Various states of the Fugacity object.
virtual void getPureGibbs(doublereal *gpure) const
Get the pure Gibbs free energies of each species.
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
virtual int reportSolnBranchActual() const
Report the solution branch which the solution is actually on.
virtual bool addSpecies(shared_ptr< Species > spec)
virtual void invalidateCache()
Invalidate any cached values which are normally updated only when a change in state is detected...
virtual void invalidateCache()
Invalidate any cached values which are normally updated only when a change in state is detected...
vector_fp m_s0_R
Temporary storage for dimensionless reference state entropies.
virtual void getGibbs_ref(doublereal *g) const
Returns the vector of the Gibbs function of the reference state at the current temperature of the sol...
virtual doublereal critPressure() const
Critical pressure (Pa).
virtual void setTemperature(const doublereal temp)
Set the temperature of the phase.
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
Header file for a derived class of ThermoPhase that handles non-ideal mixtures based on the fugacity ...
virtual void getCp_R_ref(doublereal *cprt) const
Returns the vector of nondimensional constant pressure heat capacities of the reference state at the ...
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
virtual doublereal densSpinodalLiquid() const
Return the value of the density at the liquid spinodal point (on the liquid side) for the current tem...
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
void setMoleFractionsByName(const compositionMap &xMap)
Set the species mole fractions by name.
void setMassFractionsByName(const compositionMap &yMap)
Set the species mass fractions by name.
Base class for exceptions thrown by Cantera classes.
virtual void setMoleFractions(const doublereal *const x)
Set the mole fractions to the specified values.
virtual void setForcedSolutionBranch(int solnBranch)
Set the solution branch to force the ThermoPhase to exist on one branch or another.
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity.
virtual void setPressure(doublereal p)
Set the internally stored pressure (Pa) at constant temperature and composition.
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
virtual doublereal densSpinodalGas() const
Return the value of the density at the gas spinodal point (on the gas side) for the current temperatu...
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
virtual doublereal critTemperature() const
Critical temperature (K).
virtual doublereal pressureCalc(doublereal TKelvin, doublereal molarVol) const
Calculate the pressure given the temperature and the molar volume.
int iState_
Current state of the fluid.
virtual void getEntropy_R_ref(doublereal *er) const
Returns the vector of nondimensional entropies of the reference state at the current temperature of t...
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
doublereal m_Tlast_ref
The last temperature at which the reference state thermodynamic properties were calculated at...
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
virtual doublereal refPressure() const
Returns the reference pressure in Pa.
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
int phaseState(bool checkState=false) const
Returns the Phase State flag for the current state of the object.
Contains declarations for string manipulation functions within Cantera.
int corr0(doublereal TKelvin, doublereal pres, doublereal &densLiq, doublereal &densGas, doublereal &liqGRT, doublereal &gasGRT)
Utility routine in the calculation of the saturation pressure.
doublereal getFloat(const XML_Node &parent, const std::string &name, const std::string &type)
Get a floating-point value from a child element.
const vector_fp & gibbs_RT_ref() const
Returns the vector of nondimensional Gibbs free energies of the reference state at the current temper...
vector_fp moleFractions_
Storage for the current values of the mole fractions of the species.
virtual void getEnthalpy_RT(doublereal *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of each species in their standard states at the current T and P of the solution...
size_t m_kk
Number of species in the phase.
virtual void setState_TP(doublereal T, doublereal pres)
Set the temperature (K) and pressure (Pa)
virtual void getStandardVolumes_ref(doublereal *vol) const
Get the molar volumes of the species reference states at the current T and P_ref of the solution...
virtual void getGibbs_RT(doublereal *grt) const
Get the nondimensional Gibbs functions for the species at their standard states of solution at the cu...
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the standard state of the species at ...
Namespace for the Cantera kernel.
virtual doublereal critDensity() const
Critical density (kg/m3).
virtual void getIntEnergy_RT(doublereal *urt) const
Returns the vector of nondimensional internal Energies of the standard state at the current temperatu...
virtual void setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase.
vector_fp m_h0_RT
Temporary storage for dimensionless reference state enthalpies.
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Enthalpy functions for the standard state species at the current T an...
vector_fp m_cp0_R
Temporary storage for dimensionless reference state heat capacities.
virtual void setStateFromXML(const XML_Node &state)
Set the initial state of the phase to the conditions specified in the state XML element.
virtual doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const
Estimate for the molar volume of the liquid.