From 55cd65fc93378425edaf384219bef6eee293b494 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Fri, 12 Sep 2008 21:51:04 +0000 Subject: [PATCH] Internal Upgrades: Added routines to calculate the spinodal curves for the water object. --- Cantera/src/thermo/PDSS.cpp | 4 + Cantera/src/thermo/PDSS.h | 7 + Cantera/src/thermo/PDSS_ConstVol.cpp | 14 +- Cantera/src/thermo/PDSS_ConstVol.h | 11 + Cantera/src/thermo/PDSS_HKFT.cpp | 4 +- Cantera/src/thermo/PDSS_IdealGas.cpp | 5 + Cantera/src/thermo/PDSS_IdealGas.h | 13 +- Cantera/src/thermo/PDSS_Water.cpp | 27 +- Cantera/src/thermo/PDSS_Water.h | 16 +- Cantera/src/thermo/WaterPropsIAPWS.cpp | 480 ++++++++++++++++++---- Cantera/src/thermo/WaterPropsIAPWS.h | 110 +++-- Cantera/src/thermo/WaterPropsIAPWSphi.cpp | 4 +- Cantera/src/thermo/WaterPropsIAPWSphi.h | 2 + 13 files changed, 565 insertions(+), 132 deletions(-) diff --git a/Cantera/src/thermo/PDSS.cpp b/Cantera/src/thermo/PDSS.cpp index 9767a5511..3d4da4b31 100644 --- a/Cantera/src/thermo/PDSS.cpp +++ b/Cantera/src/thermo/PDSS.cpp @@ -431,6 +431,10 @@ namespace Cantera { err("setState_TP()"); } + void PDSS::setState_TR(doublereal temp, doublereal rho) { + err("setState_TR()"); + } + /// saturation pressure doublereal PDSS::satPressure(doublereal t){ err("satPressure()"); diff --git a/Cantera/src/thermo/PDSS.h b/Cantera/src/thermo/PDSS.h index 66aad718c..6df1ef860 100644 --- a/Cantera/src/thermo/PDSS.h +++ b/Cantera/src/thermo/PDSS.h @@ -373,6 +373,13 @@ namespace Cantera { */ virtual void setState_TP(doublereal temp, doublereal pres); + //! Set the internal temperature and density + /*! + * @param temp Temperature (Kelvin) + * @param rho Density (kg m-3) + */ + virtual void setState_TR(doublereal temp, doublereal rho); + /** * @} * @name Miscellaneous properties of the standard state diff --git a/Cantera/src/thermo/PDSS_ConstVol.cpp b/Cantera/src/thermo/PDSS_ConstVol.cpp index 4b07333dd..fc96f7227 100644 --- a/Cantera/src/thermo/PDSS_ConstVol.cpp +++ b/Cantera/src/thermo/PDSS_ConstVol.cpp @@ -326,7 +326,7 @@ namespace Cantera { m_gss_RT_ptr[m_spindex] = m_hss_RT_ptr[m_spindex] - m_sss_R_ptr[m_spindex]; } - void PDSS_ConstVol::setTemperature(double temp) { + void PDSS_ConstVol::setTemperature(doublereal temp) { m_temp = temp; m_spthermo->update_one(m_spindex, temp, m_cp0_R_ptr, m_h0_RT_ptr, m_s0_R_ptr); @@ -342,11 +342,21 @@ namespace Cantera { } - void PDSS_ConstVol::setState_TP(double temp, double pres) { + void PDSS_ConstVol::setState_TP(doublereal temp, doublereal pres) { setTemperature(temp); setPressure(pres); } + + void PDSS_ConstVol::setState_TR(doublereal temp, doublereal rho) { + doublereal rhoStored = m_mw / m_constMolarVolume; + if (fabs(rhoStored - rho) / (rhoStored + rho) > 1.0E-4) { + throw CanteraError("PDSS_ConstVol::setState_TR", + "Inconsistent supplied rho"); + } + setTemperature(temp); + } + /// saturation pressure doublereal PDSS_ConstVol::satPressure(doublereal t){ return (1.0E-200); diff --git a/Cantera/src/thermo/PDSS_ConstVol.h b/Cantera/src/thermo/PDSS_ConstVol.h index dcae51c7e..faa5996ae 100644 --- a/Cantera/src/thermo/PDSS_ConstVol.h +++ b/Cantera/src/thermo/PDSS_ConstVol.h @@ -298,6 +298,14 @@ namespace Cantera { */ virtual void setState_TP(double temp, double pres); + + //! Set the internal temperature and density + /*! + * @param temp Temperature (Kelvin) + * @param rho Density (kg m-3) + */ + virtual void setState_TR(double temp, double rho); + /** * @} * @name Miscellaneous properties of the standard state @@ -406,6 +414,9 @@ namespace Cantera { private: //! Value of the constant molar volume for the species + /*! + * m3 / kmol + */ doublereal m_constMolarVolume; }; diff --git a/Cantera/src/thermo/PDSS_HKFT.cpp b/Cantera/src/thermo/PDSS_HKFT.cpp index f49489182..496f3d841 100644 --- a/Cantera/src/thermo/PDSS_HKFT.cpp +++ b/Cantera/src/thermo/PDSS_HKFT.cpp @@ -273,8 +273,8 @@ namespace Cantera { doublereal PDSS_HKFT::molarVolume() const { - double pbar = m_pres * 1.0E-5; - double m_presR_bar = OneAtm * 1.0E-5; + // double pbar = m_pres * 1.0E-5; + //double m_presR_bar = OneAtm * 1.0E-5; double a1term = m_a1 * 1.0E-5; diff --git a/Cantera/src/thermo/PDSS_IdealGas.cpp b/Cantera/src/thermo/PDSS_IdealGas.cpp index 6205f7723..0e191a731 100644 --- a/Cantera/src/thermo/PDSS_IdealGas.cpp +++ b/Cantera/src/thermo/PDSS_IdealGas.cpp @@ -340,6 +340,11 @@ namespace Cantera { setTemperature(temp); } + void PDSS_IdealGas::setState_TR(double temp, double rho) { + m_pres = GasConstant * temp * rho / m_mw; + setTemperature(temp); + } + /// saturation pressure doublereal PDSS_IdealGas::satPressure(doublereal t){ throw CanteraError("PDSS_IdealGas::satPressure()", "unimplemented"); diff --git a/Cantera/src/thermo/PDSS_IdealGas.h b/Cantera/src/thermo/PDSS_IdealGas.h index 963a4a50a..288c8e147 100644 --- a/Cantera/src/thermo/PDSS_IdealGas.h +++ b/Cantera/src/thermo/PDSS_IdealGas.h @@ -300,17 +300,24 @@ namespace Cantera { /*! * @param temp Temperature (Kelvin) */ - virtual void setTemperature(double temp); + virtual void setTemperature(doublereal temp); //! Return the current storred temperature - double temperature() const; + doublereal temperature() const; //! Set the internal temperature and pressure /*! * @param temp Temperature (Kelvin) * @param pres pressure (Pascals) */ - virtual void setState_TP(double temp, double pres); + virtual void setState_TP(doublereal temp, doublereal pres); + + //! Set the internal temperature and density + /*! + * @param temp Temperature (Kelvin) + * @param rho Density (Pascals) + */ + virtual void setState_TR(doublereal temp, doublereal rho); /** * @} diff --git a/Cantera/src/thermo/PDSS_Water.cpp b/Cantera/src/thermo/PDSS_Water.cpp index 5d10ac591..36c198697 100644 --- a/Cantera/src/thermo/PDSS_Water.cpp +++ b/Cantera/src/thermo/PDSS_Water.cpp @@ -366,15 +366,19 @@ namespace Cantera { } - void PDSS_Water:: - setPressure(doublereal p) { + // In this routine we must be sure to only find the water branch of the + // curve and not the gas branch + void PDSS_Water::setPressure(doublereal p) { doublereal T = m_temp; doublereal dens = m_dens; - int waterState = WATER_GAS; - doublereal rc = m_sub->Rhocrit(); - if (dens > rc) { - waterState = WATER_LIQUID; + int waterState = WATER_LIQUID; + if (T > m_sub->Tcrit()) { + waterState = WATER_SUPERCRIT; } + if (p < 1.0) { + waterState = WATER_GAS; + } + #ifdef DEBUG_HKM //printf("waterPDSS: set pres = %g t = %g, waterState = %d\n", // p, T, waterState); @@ -388,6 +392,9 @@ namespace Cantera { } m_dens = dd; m_pres = p; + + m_iState = m_sub->phaseState(true); + } // Return the volumetric thermal expansion coefficient. Units: 1/K. @@ -405,7 +412,7 @@ namespace Cantera { doublereal PDSS_Water::dthermalExpansionCoeffdT() const { doublereal pres = pressure(); doublereal dens_save = m_dens; - double tt = m_temp - 0.04; + doublereal tt = m_temp - 0.04; doublereal dd = m_sub->density(tt, pres, m_iState, m_dens); if (dd < 0.0) { throw CanteraError("PDSS_Water::dthermalExpansionCoeffdT", @@ -452,6 +459,12 @@ namespace Cantera { setPressure(pres); } + void PDSS_Water::setState_TR(doublereal temp, doublereal dens) { + m_temp = temp; + m_dens = dens; + m_sub->setState_TR(m_temp, m_dens); + } + // saturation pressure doublereal PDSS_Water::satPressure(doublereal t){ doublereal pp = m_sub->psat(t); diff --git a/Cantera/src/thermo/PDSS_Water.h b/Cantera/src/thermo/PDSS_Water.h index 0411f7213..84527bc8b 100644 --- a/Cantera/src/thermo/PDSS_Water.h +++ b/Cantera/src/thermo/PDSS_Water.h @@ -311,6 +311,14 @@ namespace Cantera { */ virtual void setState_TP(doublereal temp, doublereal pres); + + //! Set the temperature and density in the object + /*! + * @param temp Temperature (Kelvin) + * @param rho Density (kg/m3) + */ + virtual void setState_TR(doublereal temp, doublereal rho); + //! Set the density of the water phase /*! * This is a non-virtual function because it specific @@ -493,9 +501,11 @@ namespace Cantera { //! state of the fluid /*! - * 0 gas - * 1 liquid - * 2 supercrit + * 0 WATER_GAS 0 + * 1 WATER_LIQUID 1 + * 2 WATER_SUPERCRIT 2 + * 3 WATER_UNSTABLELIQUID 3 + * 4 WATER_UNSTABLEGAS */ int m_iState; diff --git a/Cantera/src/thermo/WaterPropsIAPWS.cpp b/Cantera/src/thermo/WaterPropsIAPWS.cpp index 76609924c..05fabf3aa 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWS.cpp @@ -23,13 +23,13 @@ * Critical Point values of water in mks units */ //! Critical Temperature value (kelvin) -const double T_c = 647.096; +const doublereal T_c = 647.096; //! Critical Pressure (Pascals) -static const double P_c = 22.064E6; +static const doublereal P_c = 22.064E6; //! Value of the Density at the critical point (kg m-3) -const double Rho_c = 322.; +const doublereal Rho_c = 322.; //! Molecular Weight of water that is consistent with the paper (kg kmol-1) -static const double M_water = 18.015268; +static const doublereal M_water = 18.015268; /* * Note, this is the Rgas value quoted in the paper. For consistency @@ -37,8 +37,16 @@ static const double M_water = 18.015268; * * The Ratio of R/M = 0.46151805 kJ kg-1 K-1 , which is Eqn. (6.3) in the paper. */ -//static const double Rgas = 8.314472E3; // Joules kmol-1 K-1 -static const double Rgas = 8.314371E3; // Joules kmol-1 K-1 +//static const doublereal Rgas = 8.314472E3; // Joules kmol-1 K-1 +static const doublereal Rgas = 8.314371E3; // Joules kmol-1 K-1 + +#ifndef MAX +# define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) +#endif + +#ifndef MIN +# define MIN(x,y) (( (x) < (y) ) ? (x) : (y)) +#endif WaterPropsIAPWS:: WaterPropsIAPWS() : @@ -75,7 +83,7 @@ WaterPropsIAPWS::~WaterPropsIAPWS() { } -void WaterPropsIAPWS::calcDim(double temperature, double rho) { +void WaterPropsIAPWS::calcDim(doublereal temperature, doublereal rho) { tau = T_c / temperature; delta = rho / Rho_c; /* @@ -92,10 +100,10 @@ void WaterPropsIAPWS::calcDim(double temperature, double rho) { } } -double WaterPropsIAPWS::helmholtzFE() const { - double retn = m_phi->phi(tau, delta); - double temperature = T_c/tau; - double RT = Rgas * temperature; +doublereal WaterPropsIAPWS::helmholtzFE() const { + doublereal retn = m_phi->phi(tau, delta); + doublereal temperature = T_c/tau; + doublereal RT = Rgas * temperature; return (retn * RT); } @@ -105,10 +113,10 @@ double WaterPropsIAPWS::helmholtzFE() const { * Temperature: kelvin * rho: density in kg m-3 */ -double WaterPropsIAPWS::pressure() const { - double retn = m_phi->pressureM_rhoRT(tau, delta); - double rho = delta * Rho_c; - double temperature = T_c / tau; +doublereal WaterPropsIAPWS::pressure() const { + doublereal retn = m_phi->pressureM_rhoRT(tau, delta); + doublereal rho = delta * Rho_c; + doublereal temperature = T_c / tau; return (retn * rho * Rgas * temperature/M_water); } @@ -127,10 +135,10 @@ double WaterPropsIAPWS::pressure() const { * * If a problem is encountered, a negative 1 is returned. */ -double WaterPropsIAPWS:: -density(double temperature, double pressure, int phase, double rhoguess) { +doublereal WaterPropsIAPWS::density(doublereal temperature, doublereal pressure, + int phase, doublereal rhoguess) { - double deltaGuess = 0.0; + doublereal deltaGuess = 0.0; if (rhoguess == -1.0) { if (phase != -1) { if (temperature > T_c) { @@ -161,11 +169,11 @@ density(double temperature, double pressure, int phase, double rhoguess) { } } - double p_red = pressure * M_water / (Rgas * temperature * Rho_c); + doublereal p_red = pressure * M_water / (Rgas * temperature * Rho_c); deltaGuess = rhoguess / Rho_c; setState_TR(temperature, rhoguess); - double delta_retn = m_phi->dfind(p_red, tau, deltaGuess); - double density_retn; + doublereal delta_retn = m_phi->dfind(p_red, tau, deltaGuess); + doublereal density_retn; if (delta_retn >0.0) { delta = delta_retn; @@ -186,7 +194,71 @@ density(double temperature, double pressure, int phase, double rhoguess) { return density_retn; } -double WaterPropsIAPWS::density() const { + +doublereal WaterPropsIAPWS::density_const(doublereal pressure, + int phase, doublereal rhoguess) const { + doublereal temperature = T_c / tau; + doublereal deltaGuess = 0.0; + doublereal deltaSave = delta; + if (rhoguess == -1.0) { + if (phase != -1) { + if (temperature > T_c) { + rhoguess = pressure * M_water / (Rgas * temperature); + } else { + if (phase == WATER_GAS || phase == WATER_SUPERCRIT) { + rhoguess = pressure * M_water / (Rgas * temperature); + } else if (phase == WATER_LIQUID) { + /* + * Provide a guess about the liquid density that is + * relatively high -> convergnce from above seems robust. + */ + rhoguess = 1000.; + } else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) { + throw Cantera::CanteraError("WaterPropsIAPWS::density", + "Unstable Branch finder is untested"); + } else { + throw Cantera::CanteraError("WaterPropsIAPWS::density", + "unknown state: " + Cantera::int2str(phase)); + } + } + } else { + /* + * Assume the Gas phase initial guess, if nothing is + * specified to the routine + */ + rhoguess = pressure * M_water / (Rgas * temperature); + } + + } + doublereal p_red = pressure * M_water / (Rgas * temperature * Rho_c); + deltaGuess = rhoguess / Rho_c; + + delta = deltaGuess; + m_phi->tdpolycalc(tau, delta); + // setState_TR(temperature, rhoguess); + + doublereal delta_retn = m_phi->dfind(p_red, tau, deltaGuess); + doublereal density_retn; + if (delta_retn > 0.0) { + delta = delta_retn; + + /* + * Dimensionalize the density before returning + */ + density_retn = delta_retn * Rho_c; + + } else { + density_retn = -1.0; + } + + delta = deltaSave; + m_phi->tdpolycalc(tau, delta); + return density_retn; +} + + + +doublereal WaterPropsIAPWS::density() const { return (delta * Rho_c); } @@ -201,9 +273,9 @@ double WaterPropsIAPWS::density() const { * return: * psat (Pascals) */ -double WaterPropsIAPWS::psat_est(double temperature) { +doublereal WaterPropsIAPWS::psat_est(doublereal temperature) const { - static const double A[8] = { + static const doublereal A[8] = { -7.8889166E0, 2.5514255E0, -6.716169E0, @@ -213,20 +285,20 @@ double WaterPropsIAPWS::psat_est(double temperature) { -148.39348E0, 48.631602E0 }; - double ps; + doublereal ps; if (temperature < 314.) { - double pl = 6.3573118E0 - 8858.843E0 / temperature + doublereal pl = 6.3573118E0 - 8858.843E0 / temperature + 607.56335E0 * pow(temperature, -0.6); ps = 0.1 * exp(pl); } else { - double v = temperature / 647.25; - double w = fabs(1.0-v); - double b = 0.0; + doublereal v = temperature / 647.25; + doublereal w = fabs(1.0-v); + doublereal b = 0.0; for (int i = 0; i < 8; i++) { - double z = i + 1; + doublereal z = i + 1; b += A[i] * pow(w, ((z+1.0)/2.0)); } - double q = b / v; + doublereal q = b / v; ps = 22.093*exp(q); } /* @@ -240,31 +312,35 @@ double WaterPropsIAPWS::psat_est(double temperature) { * Returns the coefficient of isothermal compressibility * of temperature and pressure. * kappa = - d (ln V) / dP at constant T. - * */ -double WaterPropsIAPWS::isothermalCompressibility() const { - double retn = m_phi->dimdpdrho(tau, delta); - double temperature = T_c/tau; - double dpdrho = retn * Rgas * temperature / M_water; - double dens = delta * Rho_c; - return (1.0 / (dens * dpdrho)); +doublereal WaterPropsIAPWS::isothermalCompressibility() const { + doublereal dpdrho_val = dpdrho(); + doublereal dens = delta * Rho_c; + return (1.0 / (dens * dpdrho_val)); } -double WaterPropsIAPWS:: coeffPresExp() const { - double retn = m_phi->dimdpdT(tau, delta); +doublereal WaterPropsIAPWS::dpdrho() const { + doublereal retn = m_phi->dimdpdrho(tau, delta); + doublereal temperature = T_c/tau; + doublereal val = retn * Rgas * temperature / M_water; + return val; +} + +doublereal WaterPropsIAPWS:: coeffPresExp() const { + doublereal retn = m_phi->dimdpdT(tau, delta); return (retn); } -double WaterPropsIAPWS:: coeffThermExp() const { - double kappa = isothermalCompressibility(); - double beta = coeffPresExp(); - double dens = delta * Rho_c; +doublereal WaterPropsIAPWS:: coeffThermExp() const { + doublereal kappa = isothermalCompressibility(); + doublereal beta = coeffPresExp(); + doublereal dens = delta * Rho_c; return (kappa * dens * Rgas * beta / M_water); } -double WaterPropsIAPWS::Gibbs() const { - double gRT = m_phi->gibbs_RT(); - double temperature = T_c/tau; +doublereal WaterPropsIAPWS::Gibbs() const { + doublereal gRT = m_phi->gibbs_RT(); + doublereal temperature = T_c/tau; return (gRT * Rgas * temperature); } @@ -273,41 +349,53 @@ double WaterPropsIAPWS::Gibbs() const { * J kmol-1 K-1. */ void WaterPropsIAPWS:: -corr(double temperature, double pressure, double &densLiq, - double &densGas, double &delGRT) { +corr(doublereal temperature, doublereal pressure, doublereal &densLiq, + doublereal &densGas, doublereal &delGRT) { densLiq = density(temperature, pressure, WATER_LIQUID, densLiq); if (densLiq <= 0.0) { - printf("error liq\n"); - exit(-1); + throw Cantera::CanteraError("WaterPropsIAPWS::corr", + "Error occurred trying to find liquid density at (T,P) = " + + Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure)); } setState_TR(temperature, densLiq); - double gibbsLiqRT = m_phi->gibbs_RT(); + doublereal gibbsLiqRT = m_phi->gibbs_RT(); densGas = density(temperature, pressure, WATER_GAS, densGas); if (densGas <= 0.0) { - printf("error gas\n"); - exit(-1); + throw Cantera::CanteraError("WaterPropsIAPWS::corr", + "Error occurred trying to find gas density at (T,P) = " + + Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure)); } setState_TR(temperature, densGas); - double gibbsGasRT = m_phi->gibbs_RT(); + doublereal gibbsGasRT = m_phi->gibbs_RT(); delGRT = gibbsLiqRT - gibbsGasRT; } void WaterPropsIAPWS:: -corr1(double temperature, double pressure, double &densLiq, - double &densGas, double &pcorr) { +corr1(doublereal temperature, doublereal pressure, doublereal &densLiq, + doublereal &densGas, doublereal &pcorr) { densLiq = density(temperature, pressure, WATER_LIQUID, densLiq); + if (densLiq <= 0.0) { + throw Cantera::CanteraError("WaterPropsIAPWS::corr1", + "Error occurred trying to find liquid density at (T,P) = " + + Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure)); + } setState_TR(temperature, densLiq); - double prL = m_phi->phiR(); + doublereal prL = m_phi->phiR(); densGas = density(temperature, pressure, WATER_GAS, densGas); + if (densGas <= 0.0) { + throw Cantera::CanteraError("WaterPropsIAPWS::corr1", + "Error occurred trying to find gas density at (T,P) = " + + Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure)); + } setState_TR(temperature, densGas); - double prG = m_phi->phiR(); + doublereal prG = m_phi->phiR(); - double rhs = (prL - prG) + log(densLiq/densGas); + doublereal rhs = (prL - prG) + log(densLiq/densGas); rhs /= (1.0/densGas - 1.0/densLiq); pcorr = rhs * Rgas * temperature / M_water; @@ -318,15 +406,15 @@ corr1(double temperature, double pressure, double &densLiq, * p : Pascals : Newtons/m**2 */ static int method = 1; -double WaterPropsIAPWS::psat(double temperature) { - double densLiq = -1.0, densGas = -1.0, delGRT = 0.0; - double dp, pcorr; - double p = psat_est(temperature); +doublereal WaterPropsIAPWS::psat(doublereal temperature) { + doublereal densLiq = -1.0, densGas = -1.0, delGRT = 0.0; + doublereal dp, pcorr; + doublereal p = psat_est(temperature); bool conv = false; for (int i = 0; i < 30; i++) { if (method == 1) { corr(temperature, p, densLiq, densGas, delGRT); - double delV = M_water * (1.0/densLiq - 1.0/densGas); + doublereal delV = M_water * (1.0/densLiq - 1.0/densGas); dp = - delGRT * Rgas * temperature / delV; } else { corr1(temperature, p, densLiq, densGas, pcorr); @@ -347,15 +435,245 @@ double WaterPropsIAPWS::psat(double temperature) { return p; } -int WaterPropsIAPWS::phaseState() const { +int WaterPropsIAPWS::phaseState(bool checkState) const { + if (checkState) { + if (tau <= 1.0) { + iState = WATER_SUPERCRIT; + } else { + doublereal T = T_c / tau; + doublereal rho = delta * Rho_c; + //doublereal psatTable = psat_est(T); + doublereal rhoMidAtm = 0.5 * (1.01E5 * M_water / (8314.472 * 373.15) + 1.0E3); + doublereal rhoMid = Rho_c + (T - T_c) * (Rho_c - rhoMidAtm) / (T_c - 373.15); + int iStateGuess = WATER_LIQUID; + if (rho < rhoMid) { + iStateGuess = WATER_GAS; + } + doublereal kappa = isothermalCompressibility(); + if (kappa >= 0.0) { + iState = iStateGuess; + } else { + // When we are here we are between the spinodal curves + doublereal rhoDel = rho * 1.000001; + + //setState_TR(T, rhoDel); + doublereal deltaSave = delta; + doublereal deltaDel = rhoDel / Rho_c; + delta = deltaDel; + m_phi->tdpolycalc(tau, deltaDel); + + doublereal kappaDel = isothermalCompressibility(); + doublereal d2rhodp2 = (rhoDel * kappaDel - rho * kappa) / (rhoDel - rho); + if (d2rhodp2 > 0.0) { + iState = WATER_UNSTABLELIQUID; + } else { + iState = WATER_UNSTABLEGAS; + } + //setState_TR(T, rho); + delta = deltaSave; + + m_phi->tdpolycalc(tau, delta); + } + } + } return iState; } + +// Find the water spinodal density +doublereal WaterPropsIAPWS::densSpinodalWater() const { + doublereal temperature = T_c/tau; + doublereal delta_save = delta; + // return the critical density if we are above or even just a little below + // the critical temperature. We just don't want to worry about the critical + // point at this juncture. + if (temperature >= T_c - 0.001) { + return Rho_c; + } + doublereal p = psat_est(temperature); + doublereal rho_low = 0.0; + doublereal rho_high = 1000; + + doublereal densSatLiq = density_const(p, WATER_LIQUID); + doublereal dens_old = densSatLiq; + delta = dens_old / Rho_c; + m_phi->tdpolycalc(tau, delta); + doublereal dpdrho_old = dpdrho(); + if (dpdrho_old > 0.0) { + rho_high = MIN(dens_old, rho_high); + } else { + rho_low = MAX(rho_low, dens_old); + } + doublereal dens_new = densSatLiq* (1.0001); + delta = dens_new / Rho_c; + m_phi->tdpolycalc(tau, delta); + doublereal dpdrho_new = dpdrho(); + if (dpdrho_new > 0.0) { + rho_high = MIN(dens_new, rho_high); + } else { + rho_low = MAX(rho_low, dens_new); + } + bool conv = false; + + for (int it = 0; it < 50; it++) { + doublereal slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old); + if (slope >= 0.0) { + slope = MAX(slope, dpdrho_new *5.0/ dens_new); + } else { + slope = -dpdrho_new; + //slope = MIN(slope, dpdrho_new *5.0 / dens_new); + // shouldn't be here for liquid spinodal + } + doublereal delta_rho = - dpdrho_new / slope; + if (delta_rho > 0.0) { + delta_rho = MIN(delta_rho, dens_new * 0.1); + } else { + delta_rho = MAX(delta_rho, - dens_new * 0.1); + } + doublereal dens_est = dens_new + delta_rho; + if (dens_est < rho_low) { + dens_est = 0.5 * (rho_low + dens_new); + } + if (dens_est > rho_high) { + dens_est = 0.5 * (rho_high + dens_new); + } + + + dens_old = dens_new; + dpdrho_old = dpdrho_new; + dens_new = dens_est; + + delta = dens_new / Rho_c; + m_phi->tdpolycalc(tau, delta); + dpdrho_new = dpdrho(); + if (dpdrho_new > 0.0) { + rho_high = MIN(dens_new, rho_high); + } else if (dpdrho_new < 0.0) { + rho_low = MAX(rho_low, dens_new); + } else { + conv = true; + break; + } + + if (fabs(dpdrho_new) < 1.0E-5) { + conv = true; + break; + } + } + + if (!conv) { + throw Cantera::CanteraError(" WaterPropsIAPWS::densSpinodalWater()", + " convergence failure"); + } + // Restore the original delta + delta = delta_save; + m_phi->tdpolycalc(tau, delta); + + return dens_new; +} + + +// Find the steam spinodal density +doublereal WaterPropsIAPWS::densSpinodalSteam() const { + doublereal temperature = T_c/tau; + doublereal delta_save = delta; + // return the critical density if we are above or even just a little below + // the critical temperature. We just don't want to worry about the critical + // point at this juncture. + if (temperature >= T_c - 0.001) { + return Rho_c; + } + doublereal p = psat_est(temperature); + doublereal rho_low = 0.0; + doublereal rho_high = 1000; + + doublereal densSatGas = density_const(p, WATER_GAS); + doublereal dens_old = densSatGas; + delta = dens_old / Rho_c; + m_phi->tdpolycalc(tau, delta); + doublereal dpdrho_old = dpdrho(); + if (dpdrho_old < 0.0) { + rho_high = MIN(dens_old, rho_high); + } else { + rho_low = MAX(rho_low, dens_old); + } + doublereal dens_new = densSatGas * (0.99); + delta = dens_new / Rho_c; + m_phi->tdpolycalc(tau, delta); + doublereal dpdrho_new = dpdrho(); + if (dpdrho_new < 0.0) { + rho_high = MIN(dens_new, rho_high); + } else { + rho_low = MAX(rho_low, dens_new); + } + bool conv = false; + + for (int it = 0; it < 50; it++) { + doublereal slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old); + if (slope >= 0.0) { + slope = dpdrho_new; + //slope = MAX(slope, dpdrho_new *5.0/ dens_new); + // shouldn't be here for gas spinodal + } else { + //slope = -dpdrho_new; + slope = MIN(slope, dpdrho_new *5.0 / dens_new); + + } + doublereal delta_rho = - dpdrho_new / slope; + if (delta_rho > 0.0) { + delta_rho = MIN(delta_rho, dens_new * 0.1); + } else { + delta_rho = MAX(delta_rho, - dens_new * 0.1); + } + doublereal dens_est = dens_new + delta_rho; + if (dens_est < rho_low) { + dens_est = 0.5 * (rho_low + dens_new); + } + if (dens_est > rho_high) { + dens_est = 0.5 * (rho_high + dens_new); + } + + + dens_old = dens_new; + dpdrho_old = dpdrho_new; + dens_new = dens_est; + + delta = dens_new / Rho_c; + m_phi->tdpolycalc(tau, delta); + dpdrho_new = dpdrho(); + if (dpdrho_new < 0.0) { + rho_high = MIN(dens_new, rho_high); + } else if (dpdrho_new > 0.0) { + rho_low = MAX(rho_low, dens_new); + } else { + conv = true; + break; + } + + if (fabs(dpdrho_new) < 1.0E-5) { + conv = true; + break; + } + } + + if (!conv) { + throw Cantera::CanteraError(" WaterPropsIAPWS::densSpinodalSteam()", + " convergence failure"); + } + // Restore the original delta + delta = delta_save; + m_phi->tdpolycalc(tau, delta); + + return dens_new; +} + + + /** * Sets the internal state of the object to the * specified temperature and density. */ -void WaterPropsIAPWS::setState_TR(double temperature, double rho) { +void WaterPropsIAPWS::setState_TR(doublereal temperature, doublereal rho) { calcDim(temperature, rho); m_phi->tdpolycalc(tau, delta); } @@ -365,9 +683,9 @@ void WaterPropsIAPWS::setState_TR(double temperature, double rho) { * Calculate the enthalpy in mks units of * J kmol-1 K-1. */ -double WaterPropsIAPWS::enthalpy() const { - double temperature = T_c/tau; - double hRT = m_phi->enthalpy_RT(); +doublereal WaterPropsIAPWS::enthalpy() const { + doublereal temperature = T_c/tau; + doublereal hRT = m_phi->enthalpy_RT(); return (hRT * Rgas * temperature); } @@ -376,9 +694,9 @@ double WaterPropsIAPWS::enthalpy() const { * Calculate the internal Energy in mks units of * J kmol-1 K-1. */ -double WaterPropsIAPWS::intEnergy() const{ - double temperature = T_c / tau; - double uRT = m_phi->intEnergy_RT(); +doublereal WaterPropsIAPWS::intEnergy() const { + doublereal temperature = T_c / tau; + doublereal uRT = m_phi->intEnergy_RT(); return (uRT * Rgas * temperature); } @@ -386,8 +704,8 @@ double WaterPropsIAPWS::intEnergy() const{ * Calculate the enthalpy in mks units of356 * J kmol-1 K-1. */ -double WaterPropsIAPWS::entropy() const { - double sR = m_phi->entropy_R(); +doublereal WaterPropsIAPWS::entropy() const { + doublereal sR = m_phi->entropy_R(); return (sR * Rgas); } @@ -395,17 +713,17 @@ double WaterPropsIAPWS::entropy() const { * Calculate heat capacity at constant volume * J kmol-1 K-1. */ -double WaterPropsIAPWS::cv() const { - double cvR = m_phi->cv_R(); +doublereal WaterPropsIAPWS::cv() const { + doublereal cvR = m_phi->cv_R(); return (cvR * Rgas); } -double WaterPropsIAPWS::cp() const { - double cpR = m_phi->cp_R(); +doublereal WaterPropsIAPWS::cp() const { + doublereal cpR = m_phi->cp_R(); return (cpR * Rgas); } -double WaterPropsIAPWS::molarVolume() const { - double rho = delta * Rho_c; +doublereal WaterPropsIAPWS::molarVolume() const { + doublereal rho = delta * Rho_c; return (M_water / rho); } diff --git a/Cantera/src/thermo/WaterPropsIAPWS.h b/Cantera/src/thermo/WaterPropsIAPWS.h index 5ae190f8d..a9601ed5b 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.h +++ b/Cantera/src/thermo/WaterPropsIAPWS.h @@ -17,7 +17,7 @@ #define WATERPROPSIAPWS_H #include "WaterPropsIAPWSphi.h" - +#include "config.h" /** * @name Names for the phase regions * @@ -163,37 +163,37 @@ public: * @param temperature temperature (kelvin) * @param rho density (kg m-3) */ - void setState_TR(double temperature, double rho); + void setState_TR(doublereal temperature, doublereal rho); //! Calculate the Helmholtz free energy in mks units of J kmol-1 K-1, //! using the last temperature and density - double helmholtzFE() const; + doublereal helmholtzFE() const; //! Calculate the Gibbs free energy in mks units of J kmol-1 K-1. //! using the last temperature and density - double Gibbs() const; + doublereal Gibbs() const; //! Calculate the enthalpy in mks units of J kmol-1 //! using the last temperature and density - double enthalpy() const; + doublereal enthalpy() const; //! Calculate the internal energy in mks units of J kmol-1 - double intEnergy() const; + doublereal intEnergy() const; //! Calculate the entropy in mks units of J kmol-1 K-1 - double entropy() const; + doublereal entropy() const; //! Calculate the constant volume heat capacity in mks units of J kmol-1 K-1 //! at the last temperature and density - double cv() const; + doublereal cv() const; //! Calculate the constant pressure heat capacity in mks units of J kmol-1 K-1 //! at the last temperature and density - double cp() const; + doublereal cp() const; //! Calculate the molar volume (kmol m-3) //! at the last temperature and density - double molarVolume() const; + doublereal molarVolume() const; //! Calculates the pressure (Pascals), given the current value of the //! temperature and density. @@ -203,7 +203,7 @@ public: * @return * returns the pressure (Pascal) */ - double pressure() const; + doublereal pressure() const; //! Calculates the density given the temperature and the pressure, //! and a guess at the density. Sets the internal state. @@ -231,14 +231,41 @@ public: * Returns the density. If an error is encountered in the calculation * the value of -1.0 is returned. */ - double density(double temperature, double pressure, - int phase = -1, double rhoguess = -1.0); + doublereal density(doublereal temperature, doublereal pressure, + int phase = -1, doublereal rhoguess = -1.0); + + //! Calculates the density given the temperature and the pressure, + //! and a guess at the density, while not changing the internal state + /*! + * Note, below T_c, this is a multivalued function. + * + * The #density() function calculates the density that is consistent with + * a particular value of the temperature and pressure. It may therefore be + * multivalued or potentially there may be no answer from this function. It therefore + * takes a phase guess and a density guess as optional parameters. If no guesses are + * supplied to density(), a gas phase guess is assumed. This may or may not be what + * is wanted. Therefore, density() should usually at leat be supplied with a phase + * guess so that it may manufacture an appropriate density guess. + * #density() manufactures the initial density guess, nondimensionalizes everything, + * and then calls #WaterPropsIAPWSphi::dfind(), which does the iterative calculation + * to find the density condition that matches the desired input pressure. + * + * @param pressure : Pressure in Pascals (Newton/m**2) + * @param phase : guessed phase of water + * : -1: no guessed phase + * @param rhoguess : guessed density of the water + * : -1.0 no guessed density + * @return + * Returns the density. If an error is encountered in the calculation + * the value of -1.0 is returned. + */ + doublereal density_const(doublereal pressure, int phase = -1, doublereal rhoguess = -1.0) const; //! Returns the density (kg m-3) /*! * The density is an independent variable in the underlying equation of state */ - double density() const; + doublereal density() const; //! Returns the coefficient of thermal expansion. /*! @@ -249,7 +276,7 @@ public: * @return * Returns the coefficient of thermal expansion */ - double coeffThermExp() const; + doublereal coeffThermExp() const; //! Returns the isochoric pressure-temperature coefficient /*! @@ -261,7 +288,7 @@ public: * beta = delta (phi0_d() + phiR_d()) * - tau delta (phi0_dt() + phiR_dt()) */ - double coeffPresExp() const; + doublereal coeffPresExp() const; //! Returns the coefficient of isothermal compressibility for the //! state of the object @@ -273,7 +300,17 @@ public: * @return * returns the isothermal compressibility */ - double isothermalCompressibility() const; + doublereal isothermalCompressibility() const; + + //! Returns the value of dp / drho at constant T for the + //! state of the object + /*! + * units - Joules / kg + * + * @return + * returns dpdrho + */ + doublereal dpdrho() const; //! This function returns an estimated value for the saturation pressure. /*! @@ -285,7 +322,7 @@ public: * @return * Returns the estimated saturation pressure */ - double psat_est(double temperature); + doublereal psat_est(doublereal temperature) const; //! This function returns the saturation pressure given the //! temperature as an input parameter. @@ -295,34 +332,43 @@ public: * Returns the saturation pressure * units = Pascal */ - double psat(double temperature); + doublereal psat(doublereal temperature); + + //! Return the value of the density at the water spinodal point + //! for the current temperature. + /*! + * + */ + doublereal densSpinodalWater() const; + + doublereal densSpinodalSteam() const; //! Returns the Phase State flag for the current state of the object /*! * There are three values: * WATER_GAS below the critical temperature but below the critical density * WATER_LIQUID below the critical temperature but above the critical density - * WATER_CRIT above the critical temperature + * WATER_SUPERCRIT above the critical temperature */ - int phaseState() const ; + int phaseState(bool checkState = false) const ; //! Returns the critical temperature of water (Kelvin) /*! * This is hard coded to the value 647.096 Kelvin */ - double Tcrit() { return 647.096;} + doublereal Tcrit() { return 647.096;} //! Returns the critical pressure of water (22.064E6 Pa) /*! * This is hard coded to the value of 22.064E6 pascals */ - double Pcrit() { return 22.064E6;} + doublereal Pcrit() { return 22.064E6;} //! Return the critical density of water (kg m-3) /*! * This is equal to 322 kg m-3. */ - double Rhocrit() { return 322.;} + doublereal Rhocrit() { return 322.;} private: /** @@ -331,7 +377,7 @@ private: * @param temperature input temperature (kelvin) * @param rho density in kg m-3 */ - void calcDim(double temperature, double rho); + void calcDim(doublereal temperature, doublereal rho); //! Utility routine in the calculation of the saturation pressure /*! @@ -341,8 +387,8 @@ private: * @param densGas output Density of gas * @param delGRT output delGRT */ - void corr(double temperature, double pressure, double &densLiq, - double &densGas, double &delGRT); + void corr(doublereal temperature, doublereal pressure, doublereal &densLiq, + doublereal &densGas, doublereal &delGRT); //! Utility routine in the calculation of the saturation pressure /*! @@ -352,8 +398,8 @@ private: * @param densGas output Density of gas * @param pcorr output corrected pressure */ - void corr1(double temperature, double pressure, double &densLiq, - double &densGas, double &pcorr); + void corr1(doublereal temperature, doublereal pressure, doublereal &densLiq, + doublereal &densGas, doublereal &pcorr); private: @@ -364,15 +410,15 @@ private: /*! * tau = T_C / T */ - double tau; + doublereal tau; //! Dimensionless density /*! * delta = rho / rho_c */ - double delta; + mutable doublereal delta; //! Current state of the system - int iState; + mutable int iState; }; #endif diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp index c05dc8ced..548bbd139 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp @@ -1111,7 +1111,7 @@ double WaterPropsIAPWSphi::dfind(double p_red, double tau, double deltaGuess) { bool conv = false; double deldd = dd; double pcheck = 1.0E-30 + 1.0E-8 * p_red; - for (int n = 0; n < 100; n++) { + for (int n = 0; n < 200; n++) { /* * Calculate the internal polynomials, and then calculate the * phi deriv functions needed by this routine. @@ -1155,7 +1155,7 @@ double WaterPropsIAPWSphi::dfind(double p_red, double tau, double deltaGuess) { if (n < 10) { dpdx = dpddelta * 1.1; } - if (dpdx < 0.1) dpdx = 0.1; + if (dpdx < 0.001) dpdx = 0.001; /* * Formulate the update to reduced density using diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.h b/Cantera/src/thermo/WaterPropsIAPWSphi.h index e75164964..ff1efd4f0 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.h +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.h @@ -17,6 +17,8 @@ #ifndef WATERPROPSIAPWSPHI_H #define WATERPROPSIAPWSPHI_H +#include "config.h" + /*! * the WaterPropsIAPSWSphi class support low level calls for * the real description of water.