From 0f57869240d44496cc9b638e62e5c52153c90169 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 3 Sep 2008 20:29:29 +0000 Subject: [PATCH] Changed a setState function to a setState_TR function for consistency with the rest of Cantera Added an analytical derivative for dpdT for the water object. Started adding states for the water object that refer to unstable conditions within the spinodal curve. --- Cantera/src/thermo/PDSS_Water.cpp | 14 +++--- Cantera/src/thermo/ThermoFactory.cpp | 6 +-- Cantera/src/thermo/WaterPropsIAPWS.cpp | 59 +++++++++++++---------- Cantera/src/thermo/WaterPropsIAPWS.h | 25 ++++++++-- Cantera/src/thermo/WaterPropsIAPWSphi.cpp | 8 ++- Cantera/src/thermo/WaterPropsIAPWSphi.h | 17 ++++--- Cantera/src/thermo/WaterSSTP.cpp | 10 ++-- 7 files changed, 88 insertions(+), 51 deletions(-) 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; }