diff --git a/Cantera/src/thermo/PDSS_Water.cpp b/Cantera/src/thermo/PDSS_Water.cpp index 85dd7972c..dca0edb76 100644 --- a/Cantera/src/thermo/PDSS_Water.cpp +++ b/Cantera/src/thermo/PDSS_Water.cpp @@ -346,7 +346,7 @@ namespace Cantera { doublereal T = m_temp; doublereal dens0 = m_sub->density(T, m_p0); doublereal h = m_sub->enthalpy(T, dens0); - m_sub->setState(m_temp, m_dens); + m_sub->setState_TR(m_temp, m_dens); return ((h + EW_Offset - SW_Offset*T)/(T * GasConstant)); } @@ -356,7 +356,7 @@ namespace Cantera { doublereal T = m_temp; doublereal dens0 = m_sub->density(T, m_p0); doublereal h = m_sub->enthalpy(T, dens0); - m_sub->setState(m_temp, m_dens); + m_sub->setState_TR(m_temp, m_dens); return ((h + EW_Offset)/(T * GasConstant)); } @@ -365,7 +365,7 @@ namespace Cantera { doublereal T = m_temp; doublereal dens0 = m_sub->density(T, m_p0); doublereal s = m_sub->entropy(T, dens0); - m_sub->setState(m_temp, m_dens); + m_sub->setState_TR(m_temp, m_dens); return ((s + SW_Offset)/GasConstant); } @@ -374,7 +374,7 @@ namespace Cantera { doublereal T = m_temp; doublereal dens0 = m_sub->density(T, m_p0); doublereal cp = m_sub->cp(T, dens0); - m_sub->setState(m_temp, m_dens); + m_sub->setState_TR(m_temp, m_dens); return (cp/GasConstant); } @@ -383,7 +383,7 @@ namespace Cantera { doublereal T = m_temp; doublereal dens0 = m_sub->density(T, m_p0); doublereal mv = m_sub->molarVolume(T, dens0); - m_sub->setState(m_temp, m_dens); + m_sub->setState_TR(m_temp, m_dens); return (mv); } @@ -465,7 +465,7 @@ namespace Cantera { void PDSS_Water::setDensity(doublereal dens) { m_dens = dens; - m_sub->setState(m_temp, m_dens); + m_sub->setState_TR(m_temp, m_dens); } doublereal PDSS_Water::density() const { @@ -475,7 +475,7 @@ namespace Cantera { void PDSS_Water::setTemperature(doublereal temp) { m_temp = temp; doublereal dd = m_dens; - m_sub->setState(temp, dd); + m_sub->setState_TR(temp, dd); } void PDSS_Water::setState_TP(doublereal temp, doublereal pres) { diff --git a/Cantera/src/thermo/ThermoFactory.cpp b/Cantera/src/thermo/ThermoFactory.cpp index 9e8775d55..797457511 100644 --- a/Cantera/src/thermo/ThermoFactory.cpp +++ b/Cantera/src/thermo/ThermoFactory.cpp @@ -440,7 +440,7 @@ namespace Cantera { skip = true; else throw CanteraError("importPhase", - "duplicate species: "+name); + "duplicate species: \"" + name + "\""); } if (!skip) { declared[name] = true; @@ -453,8 +453,8 @@ namespace Cantera { ++k; } else { - throw CanteraError("importPhase","no data for species " - +name); + throw CanteraError("importPhase","no data for species, \"" + + name + "\""); } } } diff --git a/Cantera/src/thermo/WaterPropsIAPWS.cpp b/Cantera/src/thermo/WaterPropsIAPWS.cpp index 6f3e16346..03a16c036 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWS.cpp @@ -14,6 +14,8 @@ */ #include "WaterPropsIAPWS.h" +#include "ctexceptions.h" +#include "stringUtils.h" #include #include #include @@ -95,7 +97,7 @@ void WaterPropsIAPWS::calcDim(double temperature, double rho) { * J kmol-1 K-1. */ double WaterPropsIAPWS::helmholtzFE(double temperature, double rho) { - setState(temperature, rho); + setState_TR(temperature, rho); double retn = m_phi->phi(tau, delta); double RT = Rgas * temperature; return (retn * RT); @@ -155,14 +157,20 @@ density(double temperature, double pressure, int phase, double rhoguess) { if (temperature > T_c) { rhoguess = pressure * M_water / (Rgas * temperature); } else { - if (phase != WATER_LIQUID) { + if (phase == WATER_GAS || phase == WATER_SUPERCRIT) { rhoguess = pressure * M_water / (Rgas * temperature); - } else { + } 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 { @@ -176,7 +184,7 @@ density(double temperature, double pressure, int phase, double rhoguess) { } double p_red = pressure * M_water / (Rgas * temperature * Rho_c); deltaGuess = rhoguess / Rho_c; - setState(temperature, rhoguess); + setState_TR(temperature, rhoguess); double delta_retn = m_phi->dfind(p_red, tau, deltaGuess); double density_retn; if (delta_retn >0.0) { @@ -190,7 +198,7 @@ density(double temperature, double pressure, int phase, double rhoguess) { * Set the internal state -> this may be * a duplication. However, let's just be sure. */ - setState(temperature, density_retn); + setState_TR(temperature, density_retn); } else { @@ -302,12 +310,17 @@ double WaterPropsIAPWS::isothermalCompressibility() const { return (1.0 / (dens * dpdrho)); } +double WaterPropsIAPWS:: coeffPresExp() const { + double retn = m_phi->dimdpdT(tau, delta); + return (retn); +} + /* * Calculate the Gibbs free energy in mks units of * J kmol-1 K-1. */ double WaterPropsIAPWS::Gibbs(double temperature, double rho) { - setState(temperature, rho); + setState_TR(temperature, rho); double gRT = m_phi->gibbs_RT(); return (gRT * Rgas * temperature); } @@ -331,7 +344,7 @@ corr(double temperature, double pressure, double &densLiq, printf("error liq\n"); exit(-1); } - setState(temperature, densLiq); + setState_TR(temperature, densLiq); double gibbsLiqRT = m_phi->gibbs_RT(); densGas = density(temperature, pressure, WATER_GAS, densGas); @@ -339,7 +352,7 @@ corr(double temperature, double pressure, double &densLiq, printf("error gas\n"); exit(-1); } - setState(temperature, densGas); + setState_TR(temperature, densGas); double gibbsGasRT = m_phi->gibbs_RT(); delGRT = gibbsLiqRT - gibbsGasRT; @@ -350,11 +363,11 @@ corr1(double temperature, double pressure, double &densLiq, double &densGas, double &pcorr) { densLiq = density(temperature, pressure, WATER_LIQUID, densLiq); - setState(temperature, densLiq); + setState_TR(temperature, densLiq); double prL = m_phi->phiR(); densGas = density(temperature, pressure, WATER_GAS, densGas); - setState(temperature, densGas); + setState_TR(temperature, densGas); double prG = m_phi->phiR(); double rhs = (prL - prG) + log(densLiq/densGas); @@ -405,7 +418,7 @@ int WaterPropsIAPWS::phaseState() const { * Sets the internal state of the object to the * specified temperature and density. */ -void WaterPropsIAPWS::setState(double temperature, double rho) { +void WaterPropsIAPWS::setState_TR(double temperature, double rho) { calcDim(temperature, rho); m_phi->tdpolycalc(tau, delta); } @@ -417,7 +430,7 @@ void WaterPropsIAPWS::setState(double temperature, double rho) { */ double WaterPropsIAPWS:: enthalpy(double temperature, double rho) { - setState(temperature, rho); + setState_TR(temperature, rho); double hRT = m_phi->enthalpy_RT(); return (hRT * Rgas * temperature); } @@ -435,7 +448,7 @@ enthalpy() const { */ double WaterPropsIAPWS:: intEnergy(double temperature, double rho) { - setState(temperature, rho); + setState_TR(temperature, rho); double uRT = m_phi->intEnergy_RT(); return (uRT * Rgas * temperature); } @@ -452,7 +465,7 @@ intEnergy() const{ */ double WaterPropsIAPWS:: entropy(double temperature, double rho) { - setState(temperature, rho); + setState_TR(temperature, rho); double sR = m_phi->entropy_R(); return (sR * Rgas); } @@ -471,7 +484,7 @@ double WaterPropsIAPWS::entropy() const { * J kmol-1 K-1. */ double WaterPropsIAPWS::cv(double temperature, double rho) { - setState(temperature, rho); + setState_TR(temperature, rho); double cvR = m_phi->cv_R(); return (cvR * Rgas); } @@ -480,27 +493,23 @@ double WaterPropsIAPWS::cv(double temperature, double rho) { * Calculate heat capacity at constant pressure * J kmol-1 K-1. */ -double WaterPropsIAPWS:: -cp(double temperature, double rho) { - setState(temperature, rho); +double WaterPropsIAPWS::cp(double temperature, double rho) { + setState_TR(temperature, rho); double cpR = m_phi->cp_R(); return (cpR * Rgas); } -double WaterPropsIAPWS:: -cp() const { +double WaterPropsIAPWS::cp() const { double cpR = m_phi->cp_R(); return (cpR * Rgas); } -double WaterPropsIAPWS:: -molarVolume(double temperature, double rho) { - setState(temperature, rho); +double WaterPropsIAPWS::molarVolume(double temperature, double rho) { + setState_TR(temperature, rho); return (M_water / rho); } -double WaterPropsIAPWS:: -molarVolume() const { +double WaterPropsIAPWS::molarVolume() const { double rho = delta * Rho_c; return (M_water / rho); } diff --git a/Cantera/src/thermo/WaterPropsIAPWS.h b/Cantera/src/thermo/WaterPropsIAPWS.h index bfed57800..62bd6a85b 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.h +++ b/Cantera/src/thermo/WaterPropsIAPWS.h @@ -22,12 +22,20 @@ * @name Names for the phase regions * * These constants are defined and used in the interface - * to describe desired phases. + * to describe the location of where we are in (T,rho) space. + * + * WATER_UNSTABLELIQUID indicates that we are in the unstable region, inside the + * spinodal curve where dpdrho < 0.0 amonst other properties. The difference + * between WATER_UNSTABLELIQUID and WATER_UNSTABLEGAS is that + * for WATER_UNSTABLELIQUID d2pdrho2 > 0 and dpdrho < 0.0 + * for WATER_UNSTABLEGAS d2pdrho2 < 0 and dpdrho < 0.0 */ //@{ #define WATER_GAS 0 #define WATER_LIQUID 1 #define WATER_SUPERCRIT 2 +#define WATER_UNSTABLELIQUID 3 +#define WATER_UNSTABLEGAS 4 //@} //! Class for calculating the equation of state of water. @@ -155,7 +163,7 @@ public: * @param temperature temperature (kelvin) * @param rho density (kg m-3) */ - void setState(double temperature, double rho); + void setState_TR(double temperature, double rho); //! Calculate the Helmholtz free energy in mks units of J kmol-1 K-1. @@ -180,7 +188,6 @@ public: //! using the last temperature and density double Gibbs() const; - //! Calculate the enthalpy in mks units of J kmol-1 /*! * @param temperature temperature (kelvin) @@ -311,6 +318,18 @@ public: */ double coeffThermExp(double temperature, double pressure); + //! Returns the isochoric pressure-temperature coefficient + /*! + * + * beta = M / (rho * Rgas) (d (pressure) / dT) at constant rho + * + * Note for ideal gases this is equal to one. + * + * beta = delta (phi0_d() + phiR_d()) + * - tau delta (phi0_dt() + phiR_dt()) + */ + double coeffPresExp() const; + //! Returns the coefficient of isothermal compressibility for the //! state of the object /*! diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp index 0e5ee73b3..c05dc8ced 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp @@ -784,7 +784,13 @@ double WaterPropsIAPWSphi::dimdpdrho(double tau, double delta) { return retn; } - +double WaterPropsIAPWSphi::dimdpdT(double tau, double delta) { + tdpolycalc(tau, delta); + double res1 = phiR_d(); + double res2 = phiR_dt(); + double retn = (1.0 + delta * res1) - tau * delta * (res2); + return retn; +} /* * Calculate d_phi0/d(tau) diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.h b/Cantera/src/thermo/WaterPropsIAPWSphi.h index a2a11f8ce..e75164964 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.h +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.h @@ -73,13 +73,6 @@ public: */ double phi_tt(double tau, double delta); - //! Second derivative of phi wrt tau, then delta - /*! - * @param tau Dimensionless temperature = T_c/T - * @param delta Dimensionless density = delta = rho / Rho_c - */ - double phi_dt(double tau, double delta); - //! Internal check # 1 void check1(); @@ -108,6 +101,16 @@ public: */ double dimdpdrho(double tau, double delta); + //! Dimensionless derivative of p wrt T at constant rho + /*! + * dp/dT * M/(Rho R) = (1.0 + delta phiR_d() + * - tau delta (phiR_dt()) + * + * @param tau Dimensionless temperature = T_c/T + * @param delta Dimensionless density = delta = rho / Rho_c + */ + double dimdpdT(double tau, double delta); + /** * This program computes the reduced density, given the reduced pressure * and the reduced temperature, tau. It takes an initial guess, deltaGuess. diff --git a/Cantera/src/thermo/WaterSSTP.cpp b/Cantera/src/thermo/WaterSSTP.cpp index dc052a78b..41ca96a5e 100644 --- a/Cantera/src/thermo/WaterSSTP.cpp +++ b/Cantera/src/thermo/WaterSSTP.cpp @@ -398,7 +398,7 @@ namespace Cantera { if (dd <= 0.0) { throw CanteraError("setPressure", "error"); } - m_sub->setState(T, dd); + m_sub->setState_TR(T, dd); doublereal g = m_sub->Gibbs(T, dd); *grt = (g + EW_Offset - SW_Offset*T)/ (GasConstant * T); dd = m_sub->density(T, p, waterState, dens); @@ -427,7 +427,7 @@ namespace Cantera { if (dd <= 0.0) { throw CanteraError("setPressure", "error"); } - m_sub->setState(T, dd); + m_sub->setState_TR(T, dd); doublereal s = m_sub->entropy(T, dd); *sr = (s + SW_Offset)/ (GasConstant); @@ -445,7 +445,7 @@ namespace Cantera { waterState = WATER_LIQUID; } doublereal dd = m_sub->density(T, OneAtm, waterState, dens); - m_sub->setState(T, dd); + m_sub->setState_TR(T, dd); if (dd <= 0.0) { throw CanteraError("setPressure", "error"); } @@ -539,7 +539,7 @@ namespace Cantera { void WaterSSTP::setTemperature(double temp) { State::setTemperature(temp); doublereal dd = density(); - m_sub->setState(temp, dd); + m_sub->setState_TR(temp, dd); } @@ -549,7 +549,7 @@ namespace Cantera { doublereal tsave = temperature(); doublereal dsave = density(); doublereal pp = m_sub->psat(t); - m_sub->setState(tsave, dsave); + m_sub->setState_TR(tsave, dsave); return pp; }