From 276602af3a5360387523e13cf44e0bb3d893e0b9 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 15 Sep 2009 15:01:39 +0000 Subject: [PATCH] Moved the water viscosity calculation to WaterProps.cpp. Changed calcDensity to be a virtual function started at VPStandardStateTP(). It's a protected function as well. At the VPStandardStateTP level, calcDensity() produces an error when called, because the EOS isn't specified. Derived classes must define this function, supplying the EOS. WaterPropsIAPW: I added documentation to this routine, defining what functions changes the state of the object. IdealMolalSoln: I had forgotten to update this object to employ the newer setState_TP() strategy that I had used in VPStandardStateTP(). --- Cantera/src/thermo/DebyeHuckel.cpp | 56 ++++---- Cantera/src/thermo/DebyeHuckel.h | 14 +- Cantera/src/thermo/HMWSoln.h | 2 +- Cantera/src/thermo/IdealMolalSoln.cpp | 49 ++----- Cantera/src/thermo/IdealMolalSoln.h | 11 ++ Cantera/src/thermo/IdealSolnGasVPSS.h | 4 +- Cantera/src/thermo/VPStandardStateTP.cpp | 5 + Cantera/src/thermo/VPStandardStateTP.h | 26 ++++ Cantera/src/thermo/WaterProps.cpp | 128 ++++++++++++++++- Cantera/src/thermo/WaterProps.h | 29 ++++ Cantera/src/thermo/WaterPropsIAPWS.cpp | 166 ++++++++++++++++++++--- Cantera/src/thermo/WaterPropsIAPWS.h | 61 +++++++-- 12 files changed, 446 insertions(+), 105 deletions(-) diff --git a/Cantera/src/thermo/DebyeHuckel.cpp b/Cantera/src/thermo/DebyeHuckel.cpp index f5b5aeb3d..9c4d78919 100644 --- a/Cantera/src/thermo/DebyeHuckel.cpp +++ b/Cantera/src/thermo/DebyeHuckel.cpp @@ -330,8 +330,37 @@ namespace Cantera { */ _updateStandardStateThermo(); + /* + * Calculate all of the other standard volumes + * -> note these are constant for now + */ + calcDensity(); + } + + /* + * Calculate the density of the mixture using the partial + * molar volumes and mole fractions as input + * + * The formula for this is + * + * \f[ + * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}} + * \f] + * + * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are + * the molecular weights, and \f$V_k\f$ are the pure species + * molar volumes. + * + * Note, the basis behind this formula is that in an ideal + * solution the partial molar volumes are equal to the pure + * species molar volumes. We have additionally specified + * in this class that the pure species molar volumes are + * independent of temperature and pressure. + * + */ + void DebyeHuckel::calcDensity() { if (m_waterSS) { - + /* * Store the internal density of the water SS. * Note, we would have to do this for all other @@ -339,42 +368,19 @@ namespace Cantera { */ m_densWaterSS = m_waterSS->density(); } - /* - * Calculate all of the other standard volumes - * -> note these are constant for now - */ - /* - * Get the partial molar volumes of all of the - * species. -> note this is a lookup for - * water, here since it was done above. - */ double *vbar = &m_pp[0]; getPartialMolarVolumes(vbar); - - /* - * Get mole fractions of all species. - */ double *x = &m_tmpV[0]; getMoleFractions(x); - - /* - * Calculate the solution molar volume and the - * solution density. - */ doublereal vtotal = 0.0; for (int i = 0; i < m_kk; i++) { vtotal += vbar[i] * x[i]; } doublereal dd = meanMolecularWeight() / vtotal; - - /* - * Now, update the State class with the results. This - * stores the density. - */ State::setDensity(dd); - } + /* * The isothermal compressibility. Units: 1/Pa. * The isothermal compressibility is defined as diff --git a/Cantera/src/thermo/DebyeHuckel.h b/Cantera/src/thermo/DebyeHuckel.h index 6c9f2f2ee..1467ec860 100644 --- a/Cantera/src/thermo/DebyeHuckel.h +++ b/Cantera/src/thermo/DebyeHuckel.h @@ -766,10 +766,10 @@ namespace Cantera { */ virtual void setPressure(doublereal p); - /** - * Calculate the density of the mixture using the partial - * molar volumes and mole fractions as input - * + protected: + //! Calculate the density of the mixture using the partial + //! molar volumes and mole fractions as input + /*! * The formula for this is * * \f[ @@ -785,12 +785,10 @@ namespace Cantera { * species molar volumes. We have additionally specified * in this class that the pure species molar volumes are * independent of temperature and pressure. - * - * NOTE: This is a non-virtual function, which is not a - * member of the ThermoPhase base class. */ - void calcDensity(); + virtual void calcDensity(); + public: //! Set the internally storred density (gm/m^3) of the phase. /*! * Overwritten setDensity() function is necessary because the diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h index ff9b5c0a5..e476d3069 100644 --- a/Cantera/src/thermo/HMWSoln.h +++ b/Cantera/src/thermo/HMWSoln.h @@ -1456,7 +1456,7 @@ namespace Cantera { */ virtual void setPressure(doublereal p); - private: + protected: /** * Calculate the density of the mixture using the partial * molar volumes and mole fractions as input diff --git a/Cantera/src/thermo/IdealMolalSoln.cpp b/Cantera/src/thermo/IdealMolalSoln.cpp index 6823f3384..280e0060a 100644 --- a/Cantera/src/thermo/IdealMolalSoln.cpp +++ b/Cantera/src/thermo/IdealMolalSoln.cpp @@ -288,55 +288,20 @@ namespace Cantera { * The mass density is not a function of pressure. */ void IdealMolalSoln::setPressure(doublereal p) { + setState_TP(temperature(), p); + } -#ifdef DEBUG_MODE - //printf("setPressure: %g\n", p); -#endif - /* - * Store the current pressure - */ - m_Pcurrent = p; - - /* - * update the standard state thermo - * -> This involves calling the water function and setting the pressure - */ - updateStandardStateThermo(); - - /* - * Calculate all of the other standard volumes - * -> note these are constant for now - */ - /* - * Get the partial molar volumes of all of the - * species. -> note this is a lookup for - * water, here since it was done above. - */ + void IdealMolalSoln::calcDensity() { double *vbar = &m_pp[0]; getPartialMolarVolumes(vbar); - - /* - * Get mole fractions of all species. - */ double *x = &m_tmpV[0]; getMoleFractions(x); - - /* - * Calculate the solution molar volume and the - * solution density. - */ doublereal vtotal = 0.0; for (int i = 0; i < m_kk; i++) { vtotal += vbar[i] * x[i]; } doublereal dd = meanMolecularWeight() / vtotal; - - /* - * Now, update the State class with the results. This - * stores the density. - */ State::setDensity(dd); - } /* @@ -409,6 +374,14 @@ namespace Cantera { } } + void IdealMolalSoln::setState_TP(doublereal temp, doublereal pres) { + State::setTemperature(temp); + m_Pcurrent = pres; + updateStandardStateThermo(); + //m_densWaterSS = m_waterSS->density(); + calcDensity(); + } + // // ------- Activities and Activity Concentrations // diff --git a/Cantera/src/thermo/IdealMolalSoln.h b/Cantera/src/thermo/IdealMolalSoln.h index 752178ab7..4a5118762 100644 --- a/Cantera/src/thermo/IdealMolalSoln.h +++ b/Cantera/src/thermo/IdealMolalSoln.h @@ -278,6 +278,7 @@ namespace Cantera { */ virtual void setPressure(doublereal p); + protected: /** * Calculate the density of the mixture using the partial * molar volumes and mole fractions as input @@ -303,6 +304,7 @@ namespace Cantera { */ void calcDensity(); + public: /** * Overwritten setDensity() function is necessary because the * density is not an indendent variable. @@ -336,6 +338,15 @@ namespace Cantera { */ void setMolarDensity(const doublereal rho); + //! Set the temperature (K) and pressure (Pa) + /*! + * Set the temperature and pressure. + * + * @param t Temperature (K) + * @param p Pressure (Pa) + */ + virtual void setState_TP(doublereal t, doublereal p); + //! The isothermal compressibility. Units: 1/Pa. /*! * The isothermal compressibility is defined as diff --git a/Cantera/src/thermo/IdealSolnGasVPSS.h b/Cantera/src/thermo/IdealSolnGasVPSS.h index 2003e055f..be8623de0 100644 --- a/Cantera/src/thermo/IdealSolnGasVPSS.h +++ b/Cantera/src/thermo/IdealSolnGasVPSS.h @@ -130,7 +130,7 @@ namespace Cantera { */ virtual doublereal isothermalCompressibility() const; - private: + protected: /** * Calculate the density of the mixture using the partial * molar volumes and mole fractions as input @@ -154,7 +154,7 @@ namespace Cantera { * NOTE: This is a non-virtual function, which is not a * member of the ThermoPhase base class. */ - void calcDensity(); + virtual void calcDensity(); public: diff --git a/Cantera/src/thermo/VPStandardStateTP.cpp b/Cantera/src/thermo/VPStandardStateTP.cpp index 5bd0a5a7e..6a79674aa 100644 --- a/Cantera/src/thermo/VPStandardStateTP.cpp +++ b/Cantera/src/thermo/VPStandardStateTP.cpp @@ -347,6 +347,10 @@ namespace Cantera { updateStandardStateThermo(); } + void VPStandardStateTP::calcDensity() { + err("VPStandardStateTP::calcDensity() called, but EOS for phase is not known"); + } + void VPStandardStateTP::setState_TP(doublereal t, doublereal pres) { /* @@ -371,6 +375,7 @@ namespace Cantera { */ //setTemperature(t); //setPressure(pres); + calcDensity(); } diff --git a/Cantera/src/thermo/VPStandardStateTP.h b/Cantera/src/thermo/VPStandardStateTP.h index 72f7cfebd..f9097565d 100644 --- a/Cantera/src/thermo/VPStandardStateTP.h +++ b/Cantera/src/thermo/VPStandardStateTP.h @@ -277,7 +277,33 @@ namespace Cantera { */ virtual void setPressure(const doublereal p); +protected: + /** + * Calculate the density of the mixture using the partial + * molar volumes and mole fractions as input + * + * The formula for this is + * + * \f[ + * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}} + * \f] + * + * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are + * the molecular weights, and \f$V_k\f$ are the pure species + * molar volumes. + * + * Note, the basis behind this formula is that in an ideal + * solution the partial molar volumes are equal to the pure + * species molar volumes. We have additionally specified + * in this class that the pure species molar volumes are + * independent of temperature and pressure. + * + * NOTE: This is a non-virtual function, which is not a + * member of the ThermoPhase base class. + */ + virtual void calcDensity(); + public: //! Set the temperature and pressure at the same time /*! * Note this function triggers a reevalulation of the standard diff --git a/Cantera/src/thermo/WaterProps.cpp b/Cantera/src/thermo/WaterProps.cpp index 9e748600d..848dfaf43 100644 --- a/Cantera/src/thermo/WaterProps.cpp +++ b/Cantera/src/thermo/WaterProps.cpp @@ -434,11 +434,26 @@ namespace Cantera { return pres; } - + // Returns the density of water + /* + * This function sets the internal temperature and pressure + * of the underlying object at the same time. + * + * @param T Temperature (kelvin) + * @param P pressure (pascal) + */ double WaterProps::density_IAPWS(double temp, double press) { + double dens = m_waterIAPWS->density(temp, press, WATER_LIQUID); + return dens; + } - double dens; - dens = m_waterIAPWS->density(temp, press, WATER_LIQUID); + // Returns the density of water + /* + * This function uses the internal state of the + * underlying water object + */ + double WaterProps::density_IAPWS() const { + double dens = m_waterIAPWS->density(); return dens; } @@ -463,4 +478,111 @@ namespace Cantera { } + + + + // Parameters for the viscosityWater() function + + const double H[4] = {1., + 0.978197, + 0.579829, + -0.202354}; + + const double Hij[6][7] = + { + { 0.5132047, 0.2151778, -0.2818107, 0.1778064, -0.04176610, 0., 0.}, + { 0.3205656, 0.7317883, -1.070786 , 0.4605040, 0., -0.01578386, 0.}, + { 0., 1.241044 , -1.263184 , 0.2340379, 0., 0., 0.}, + { 0., 1.476783 , 0., -0.4924179, 0.1600435, 0., -0.003629481}, + {-0.7782567, 0.0 , 0., 0. , 0., 0., 0.}, + { 0.1885447, 0.0 , 0., 0. , 0., 0., 0.}, + }; + const double TStar = 647.27; // Kelvin + const double rhoStar = 317.763; // kg / m3 + const double presStar = 22.115E6; // Pa + const double muStar = 55.071E-6; //Pa s + + + // Returns the viscosity of water at the current conditions + // (kg/m/s) + /* + * This function calculates the value of the viscosity of pure + * water at the current T and P. + * + * The formulas used are from the paper + * + * J. V. Sengers, J. T. R. Watson, "Improved International + * Formulations for the Viscosity and Thermal Conductivity of + * Water Substance", J. Phys. Chem. Ref. Data, 15, 1291 (1986). + * + * The formulation is accurate for all temperatures and pressures, + * for steam and for water, even near the critical point. + * Pressures above 500 MPa and temperature above 900 C are suspect. + */ + double WaterProps::viscosityWater() const { + + double temp = m_waterIAPWS->temperature(); + double dens = m_waterIAPWS->density(); + + //WaterPropsIAPWS *waterP = new WaterPropsIAPWS(); + //m_waterIAPWS->setState_TR(temp, dens); + //double pressure = m_waterIAPWS->pressure(); + //printf("pressure = %g\n", pressure); + //dens = 18.02 * pressure / (GasConstant * temp); + //printf ("mod dens = %g\n", dens); + + double rhobar = dens/rhoStar; + double tbar = temp / TStar; + // double pbar = pressure / presStar; + + double tbar2 = tbar * tbar; + double tbar3 = tbar2 * tbar; + + double mu0bar = std::sqrt(tbar) / (H[0] + H[1]/tbar + H[2]/tbar2 + H[3]/tbar3); + + //printf("mu0bar = %g\n", mu0bar); + //printf("mu0 = %g\n", mu0bar * muStar); + + double tfac1 = 1.0 / tbar - 1.0; + double tfac2 = tfac1 * tfac1; + double tfac3 = tfac2 * tfac1; + double tfac4 = tfac3 * tfac1; + double tfac5 = tfac4 * tfac1; + + double rfac1 = rhobar - 1.0; + double rfac2 = rfac1 * rfac1; + double rfac3 = rfac2 * rfac1; + double rfac4 = rfac3 * rfac1; + double rfac5 = rfac4 * rfac1; + double rfac6 = rfac5 * rfac1; + + double sum = (Hij[0][0] + Hij[1][0]*tfac1 + Hij[4][0]*tfac4 + Hij[5][0]*tfac5 + + Hij[0][1]*rfac1 + Hij[1][1]*tfac1*rfac1 + Hij[2][1]*tfac2*rfac1 + Hij[3][1]*tfac3*rfac1 + + Hij[0][2]*rfac2 + Hij[1][2]*tfac1*rfac2 + Hij[2][2]*tfac2*rfac2 + + Hij[0][3]*rfac3 + Hij[1][3]*tfac1*rfac3 + Hij[2][3]*tfac2*rfac3 + Hij[3][3]*tfac3*rfac3 + + Hij[0][4]*rfac4 + Hij[3][4]*tfac3*rfac4 + + Hij[1][5]*tfac1*rfac5 + Hij[3][6]*tfac3*rfac6 + ); + double mu1bar = std::exp(rhobar * sum); + + // Apply the near-critical point corrections if necessary + double mu2bar = 1.0; + if ((tbar >= 0.9970) && tbar <= 1.0082) { + if ((rhobar >= 0.755) && (rhobar <= 1.290)) { + double drhodp = 1.0 / m_waterIAPWS->dpdrho(); + drhodp *= presStar / rhoStar; + double xsi = rhobar * drhodp; + if (xsi >= 21.93) { + mu2bar = 0.922 * std::pow(xsi, 0.0263); + } + } + } + + double mubar = mu0bar * mu1bar * mu2bar; + + return mubar * muStar; + } + + + } diff --git a/Cantera/src/thermo/WaterProps.h b/Cantera/src/thermo/WaterProps.h index a52d65ad1..560baac3d 100644 --- a/Cantera/src/thermo/WaterProps.h +++ b/Cantera/src/thermo/WaterProps.h @@ -265,11 +265,22 @@ namespace Cantera { //! Returns the density of water /*! + * This function sets the internal temperature and pressure + * of the underlying object at the same time. + * * @param T Temperature (kelvin) * @param P pressure (pascal) */ double density_IAPWS(double T, double P); + //! Returns the density of water + /*! + * This function uses the internal state of the + * underlying water object + */ + double density_IAPWS() const; + + //! returns the coefficient of thermal expansion /*! * @param T Temperature (kelvin) @@ -284,6 +295,24 @@ namespace Cantera { */ double isothermalCompressibility_IAPWS(double T, double P); + //! Returns the viscosity of water at the current conditions + //! (kg/m/s) + /*! + * This function calculates the value of the viscosity of pure + * water at the current T and P. + * + * The formulas used are from the paper + * J. V. Sengers, J. T. R. Watson, "Improved International + * Formulations for the Viscosity and Thermal Conductivity of + * Water Substance", J. Phys. Chem. Ref. Data, 15, 1291 (1986). + * + * The formulation is accurate for all temperatures and pressures, + * for steam and for water, even near the critical point. + * Pressures above 500 MPa and temperature above 900 C are suspect. + */ + double viscosityWater() const; + + protected: //! Pointer to the WaterPropsIAPWS object diff --git a/Cantera/src/thermo/WaterPropsIAPWS.cpp b/Cantera/src/thermo/WaterPropsIAPWS.cpp index f86793130..6ba8de17c 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWS.cpp @@ -49,6 +49,7 @@ static const doublereal Rgas = 8.314371E3; // Joules kmol-1 K-1 #endif //@} +// Base constructor WaterPropsIAPWS:: WaterPropsIAPWS() : m_phi(0), tau(-1.0), @@ -58,6 +59,10 @@ WaterPropsIAPWS:: WaterPropsIAPWS() : m_phi = new WaterPropsIAPWSphi(); } +// Copy constructor +/* + * @param b Object to be copied + */ WaterPropsIAPWS::WaterPropsIAPWS(const WaterPropsIAPWS &b) : m_phi(0), tau(b.tau), @@ -68,6 +73,10 @@ WaterPropsIAPWS::WaterPropsIAPWS(const WaterPropsIAPWS &b) : m_phi->tdpolycalc(tau, delta); } +// assignment constructor +/* + * @param right Object to be copied + */ WaterPropsIAPWS & WaterPropsIAPWS::operator=(const WaterPropsIAPWS &b) { if (this == &b) return *this; tau = b.tau; @@ -77,12 +86,20 @@ WaterPropsIAPWS & WaterPropsIAPWS::operator=(const WaterPropsIAPWS &b) { return *this; } +// destructor WaterPropsIAPWS::~WaterPropsIAPWS() { delete (m_phi); m_phi = 0; } - +/* + * Calculate the dimensionless temp and rho and store internally. + * + * @param temperature input temperature (kelvin) + * @param rho density in kg m-3 + * + * this is a private function + */ void WaterPropsIAPWS::calcDim(doublereal temperature, doublereal rho) { tau = T_c / temperature; delta = rho / Rho_c; @@ -100,6 +117,8 @@ void WaterPropsIAPWS::calcDim(doublereal temperature, doublereal rho) { } } +// Calculate the Helmholtz free energy in mks units of J kmol-1 K-1, +// using the last temperature and density doublereal WaterPropsIAPWS::helmholtzFE() const { doublereal retn = m_phi->phi(tau, delta); doublereal temperature = T_c/tau; @@ -194,7 +213,32 @@ doublereal WaterPropsIAPWS::density(doublereal temperature, doublereal pressure, return density_retn; } - +// 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 WaterPropsIAPWS::density_const(doublereal pressure, int phase, doublereal rhoguess) const { doublereal temperature = T_c / tau; @@ -256,12 +300,24 @@ doublereal WaterPropsIAPWS::density_const(doublereal pressure, return density_retn; } - - +// Returns the density (kg m-3) +/* + * The density is an independent variable in the underlying equation of state + * + * @return Returns the density (kg m-3) + */ doublereal WaterPropsIAPWS::density() const { return (delta * Rho_c); } +// Returns the temperature (Kelvin) +/* + * @return Returns the internally storred temperature + */ +doublereal WaterPropsIAPWS::temperature() const { + return (T_c / tau); +} + /* * psat_est provides a rough estimate of the saturation * pressure given the temperature. This is used as an initial @@ -319,6 +375,13 @@ doublereal WaterPropsIAPWS::isothermalCompressibility() const { return (1.0 / (dens * dpdrho_val)); } +// Returns the value of dp / drho at constant T at the current +// state of the object +/* + * units - Joules / kg + * + * @return returns dpdrho + */ doublereal WaterPropsIAPWS::dpdrho() const { doublereal retn = m_phi->dimdpdrho(tau, delta); doublereal temperature = T_c/tau; @@ -326,11 +389,26 @@ doublereal WaterPropsIAPWS::dpdrho() const { return val; } +// Returns the isochoric pressure derivative wrt temperature +/* + * 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()) + */ doublereal WaterPropsIAPWS:: coeffPresExp() const { doublereal retn = m_phi->dimdpdT(tau, delta); return (retn); } +// Returns the coefficient of thermal expansion. +/* + * alpha = d (ln V) / dT at constant P. + * + * @return Returns the coefficient of thermal expansion + */ doublereal WaterPropsIAPWS:: coeffThermExp() const { doublereal kappa = isothermalCompressibility(); doublereal beta = coeffPresExp(); @@ -338,15 +416,27 @@ doublereal WaterPropsIAPWS:: coeffThermExp() const { return (kappa * dens * Rgas * beta / M_water); } +// Calculate the Gibbs free energy in mks units of J kmol-1 K-1. +// using the last temperature and density doublereal WaterPropsIAPWS::Gibbs() const { doublereal gRT = m_phi->gibbs_RT(); doublereal temperature = T_c/tau; return (gRT * Rgas * temperature); } + +// Utility routine in the calculation of the saturation pressure /* + * Private routine + * * Calculate the Gibbs free energy in mks units of * J kmol-1 K-1. + * + * @param temperature temperature (kelvin) + * @param pressure pressure (Pascal) + * @param densLiq Output density of liquid + * @param densGas output Density of gas + * @param delGRT output delGRT */ void WaterPropsIAPWS:: corr(doublereal temperature, doublereal pressure, doublereal &densLiq, @@ -373,6 +463,16 @@ corr(doublereal temperature, doublereal pressure, doublereal &densLiq, delGRT = gibbsLiqRT - gibbsGasRT; } +// Utility routine in the calculation of the saturation pressure +/* + * Private routine + * + * @param temperature temperature (kelvin) + * @param pressure pressure (Pascal) + * @param densLiq Output density of liquid + * @param densGas output Density of gas + * @param pcorr output corrected pressure + */ void WaterPropsIAPWS:: corr1(doublereal temperature, doublereal pressure, doublereal &densLiq, doublereal &densGas, doublereal &pcorr) { @@ -401,13 +501,27 @@ corr1(doublereal temperature, doublereal pressure, doublereal &densLiq, pcorr = rhs * Rgas * temperature / M_water; } -/** - * Calculate the saturation pressure given the temperature. - * p : Pascals : Newtons/m**2 - */ -static int method = 1; +// This function returns the saturation pressure given the +// temperature as an input parameter, and sets the internal state to the saturated +// conditions. +/* + * Note this function will return the saturation pressure, given the temperature. + * It will then set the state of the system to the saturation condition. The input + * parameter waterState is used to either specify the liquid state or the + * gas state at the desired temperatue and saturated pressure. + * + * If the input temperature, T, is above T_c, this routine will set the internal + * state to T and the pressure to P_c. Then, return P_c. + * + * @param temperature input temperature (kelvin) + * @param waterState integer specifying the water state + * + * @return Returns the saturation pressure + * units = Pascal + */ doublereal WaterPropsIAPWS::psat(doublereal temperature, int waterState) { + static int method = 1; doublereal densLiq = -1.0, densGas = -1.0, delGRT = 0.0; doublereal dp, pcorr; if (temperature >= T_c) { @@ -450,6 +564,16 @@ doublereal WaterPropsIAPWS::psat(doublereal temperature, int waterState) { return p; } +// Returns the Phase State flag for the current state of the object +/* + * @param checkState If true, this function does a complete check to see where + * in paramters space we are + * + * 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_SUPERCRIT above the critical temperature + */ int WaterPropsIAPWS::phaseState(bool checkState) const { if (checkState) { if (tau <= 1.0) { @@ -494,8 +618,11 @@ int WaterPropsIAPWS::phaseState(bool checkState) const { return iState; } - -// Find the water spinodal density +// Return the value of the density at the water spinodal point (on the liquid side) +// for the current temperature. +/* + * @return returns the density with units of kg m-3 + */ doublereal WaterPropsIAPWS::densSpinodalWater() const { doublereal temperature = T_c/tau; doublereal delta_save = delta; @@ -587,8 +714,11 @@ doublereal WaterPropsIAPWS::densSpinodalWater() const { return dens_new; } - -// Find the steam spinodal density +// Return the value of the density at the water spinodal point (on the gas side) +// for the current temperature. +/* + * @return returns the density with units of kg m-3 + */ doublereal WaterPropsIAPWS::densSpinodalSteam() const { doublereal temperature = T_c/tau; doublereal delta_save = delta; @@ -682,9 +812,7 @@ doublereal WaterPropsIAPWS::densSpinodalSteam() const { return dens_new; } - - -/** +/* * Sets the internal state of the object to the * specified temperature and density. */ @@ -693,7 +821,6 @@ void WaterPropsIAPWS::setState_TR(doublereal temperature, doublereal rho) { m_phi->tdpolycalc(tau, delta); } - /* * Calculate the enthalpy in mks units of * J kmol-1 K-1. @@ -704,7 +831,6 @@ doublereal WaterPropsIAPWS::enthalpy() const { return (hRT * Rgas * temperature); } - /* * Calculate the internal Energy in mks units of * J kmol-1 K-1. @@ -733,11 +859,15 @@ doublereal WaterPropsIAPWS::cv() const { return (cvR * Rgas); } +// Calculate the constant pressure heat capacity in mks units of J kmol-1 K-1 +// at the last temperature and density doublereal WaterPropsIAPWS::cp() const { doublereal cpR = m_phi->cp_R(); return (cpR * Rgas); } +// Calculate the molar volume (kmol m-3) +// at the last temperature and density 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 51b741e13..5c85ebfd5 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.h +++ b/Cantera/src/thermo/WaterPropsIAPWS.h @@ -134,6 +134,25 @@ * - WATER_LIQUID * - WATER_SUPERCRIT * + * There are only three functions which actually change the value of the internal + * state of this object after it's been instantiated + + * - setState_TR(temperature, rho) + * - density(temperature, pressure, phase, rhoguess) + * - psat(temperature, waterState); + * + * The setState_TR() is the main function that sets the temperature and rho value. + * The density() function serves as a setState_TP() function, in that it sets + * internal state to a temperature and pressure. However, note that this is potentially + * multivalued. Therefore, we need to supply in addition a phase guess and a rho guess + * to the input temperature and pressure. + * The psat() function sets the internal state to the saturated liquid or saturated gas + * state, dependeing on the waterState parameter. + * + * Because the underlying object WaterPropsIAPWSphi is privately held, you can be + * sure that the underlying state of this object doesn't change except due to the + * three function calls listed above. + * * @ingroup thermoprops * */ @@ -243,6 +262,7 @@ public: * 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. @@ -264,9 +284,17 @@ public: //! Returns the density (kg m-3) /*! * The density is an independent variable in the underlying equation of state + * + * @return Returns the density (kg m-3) */ doublereal density() const; + //! Returns the temperature (Kelvin) + /*! + * @return Returns the internally storred temperature + */ + doublereal temperature() const; + //! Returns the coefficient of thermal expansion. /*! * alpha = d (ln V) / dT at constant P. @@ -276,7 +304,7 @@ public: */ doublereal coeffThermExp() const; - //! Returns the isochoric pressure-temperature coefficient + //! Returns the isochoric pressure derivative wrt temperature /*! * * beta = M / (rho * Rgas) (d (pressure) / dT) at constant rho @@ -322,14 +350,22 @@ public: */ doublereal psat_est(doublereal temperature) const; - //! This function returns the saturation pressure given the - //! temperature as an input parameter. + //! This function returns the saturation pressure given the + //! temperature as an input parameter, and sets the internal state to the saturated + //! conditions. /*! + * Note this function will return the saturation pressure, given the temperature. + * It will then set the state of the system to the saturation condition. The input + * parameter waterState is used to either specify the liquid state or the + * gas state at the desired temperatue and saturated pressure. + * + * If the input temperature, T, is above T_c, this routine will set the internal + * state to T and the pressure to P_c. Then, return P_c. + * * @param temperature input temperature (kelvin) * @param waterState integer specifying the water state * - * @return - * Returns the saturation pressure + * @return Returns the saturation pressure * units = Pascal */ doublereal psat(doublereal temperature, int waterState = WATER_LIQUID); @@ -364,23 +400,24 @@ public: /*! * This is hard coded to the value 647.096 Kelvin */ - doublereal Tcrit() { return 647.096;} + doublereal Tcrit() const { return 647.096;} //! Returns the critical pressure of water (22.064E6 Pa) /*! * This is hard coded to the value of 22.064E6 pascals */ - doublereal Pcrit() { return 22.064E6;} + doublereal Pcrit() const { return 22.064E6;} //! Return the critical density of water (kg m-3) /*! * This is equal to 322 kg m-3. */ - doublereal Rhocrit() { return 322.;} + doublereal Rhocrit() const { return 322.;} private: - /** - * Calculate the dimensionless temp and rho and store internally. + //! Calculate the dimensionless temp and rho and store internally. + /*! + * Private routine * * @param temperature input temperature (kelvin) * @param rho density in kg m-3 @@ -389,6 +426,8 @@ private: //! Utility routine in the calculation of the saturation pressure /*! + * Private routine + * * @param temperature temperature (kelvin) * @param pressure pressure (Pascal) * @param densLiq Output density of liquid @@ -400,6 +439,8 @@ private: //! Utility routine in the calculation of the saturation pressure /*! + * Private routine + * * @param temperature temperature (kelvin) * @param pressure pressure (Pascal) * @param densLiq Output density of liquid