From 5fc2c72c39155c57ad4f9f47c5101b120e1d5a7f Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 2 Sep 2008 17:23:40 +0000 Subject: [PATCH] Made the isothermal compressibility function a real derivative function instead of a numerical derivative. --- Cantera/src/thermo/PDSS_Water.cpp | 4 +- Cantera/src/thermo/WaterProps.cpp | 7 +- Cantera/src/thermo/WaterPropsIAPWS.cpp | 78 ++++++----------------- Cantera/src/thermo/WaterPropsIAPWS.h | 10 +-- Cantera/src/thermo/WaterPropsIAPWSphi.cpp | 16 ++++- Cantera/src/thermo/WaterPropsIAPWSphi.h | 10 +++ 6 files changed, 50 insertions(+), 75 deletions(-) diff --git a/Cantera/src/thermo/PDSS_Water.cpp b/Cantera/src/thermo/PDSS_Water.cpp index 477e20e53..85dd7972c 100644 --- a/Cantera/src/thermo/PDSS_Water.cpp +++ b/Cantera/src/thermo/PDSS_Water.cpp @@ -15,7 +15,6 @@ #include "ctml.h" #include "PDSS_Water.h" #include "WaterPropsIAPWS.h" -//#include "importCTML.h" #include "ThermoFactory.h" #include @@ -451,8 +450,7 @@ namespace Cantera { } doublereal PDSS_Water::isothermalCompressibility() const { - doublereal pres = pressure(); - doublereal val = m_sub->isothermalCompressibility(m_temp, pres); + doublereal val = m_sub->isothermalCompressibility(); return val; } diff --git a/Cantera/src/thermo/WaterProps.cpp b/Cantera/src/thermo/WaterProps.cpp index 4c4e15f7e..2eac58007 100644 --- a/Cantera/src/thermo/WaterProps.cpp +++ b/Cantera/src/thermo/WaterProps.cpp @@ -417,7 +417,7 @@ namespace Cantera { double depsRelWaterdP = relEpsilon(T, P, 3); dAdP -= A_Debye * (1.5 * depsRelWaterdP / epsRelWater); - double kappa = isothermalCompressibility_IAPWS(T, P); + double kappa = isothermalCompressibility_IAPWS(T,P); //double ddwdP = density_T_new(T, P, 3); dAdP += A_Debye * (0.5 * kappa); @@ -449,9 +449,8 @@ namespace Cantera { } double WaterProps::isothermalCompressibility_IAPWS(double temp, double press) { - - double kappa; - kappa = m_waterIAPWS->isothermalCompressibility(temp, press); + m_waterIAPWS->density(temp, press, WATER_LIQUID); + double kappa = m_waterIAPWS->isothermalCompressibility(); return kappa; } diff --git a/Cantera/src/thermo/WaterPropsIAPWS.cpp b/Cantera/src/thermo/WaterPropsIAPWS.cpp index cbaa191d6..6f3e16346 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWS.cpp @@ -117,7 +117,14 @@ double WaterPropsIAPWS::pressure(double temperature, double rho) { double retn = m_phi->pressureM_rhoRT(tau, delta); return (retn * rho * Rgas * temperature/M_water); } -double WaterPropsIAPWS::pressure() const{ + +/* + * Calculate the pressure (Pascals), using the + * current internally storred temperature and density + * 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; @@ -286,65 +293,26 @@ double WaterPropsIAPWS::coeffThermExp(double temperature, double pressure) { * of temperature and pressure. * kappa = - d (ln V) / dP at constant T. * - * Currently this function is calculated using an inaccurate - * one-sided differencing scheme. */ -double WaterPropsIAPWS:: -isothermalCompressibility(double temperature, double pressure) { - /* - * Difference amount is large, because we are solving for - * density underneath - */ - double deltaP; - double psat_at=0.0; - double rhoguess = -1; - int phase = -1; - if (temperature > T_c) { - rhoguess = pressure * M_water / (Rgas * temperature); - phase = WATER_SUPERCRIT; - } else { - psat_at = psat(temperature); - if (pressure >= psat_at) { - phase = WATER_LIQUID; - } else { - phase = WATER_GAS; - } - } - double dens_base = density(temperature, pressure, phase, rhoguess); - if (dens_base == -1.0) { - printf("problems\n"); - exit(-1); - } - - if (iState == WATER_GAS) { - deltaP = -0.0001 * pressure; - } else if (iState == WATER_LIQUID) { - deltaP = +0.0001 * pressure; - } else { - deltaP = +0.0001 * pressure; - } - - - double pres_del = pressure + deltaP; - double dens_del = density(temperature, pres_del, phase, dens_base); - double Vavg = 0.5 * (1./dens_del + 1./dens_base); - double retn = -1.0 / Vavg * (1./dens_del - 1.0/dens_base)/deltaP; - return retn; +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)); } - /* * Calculate the Gibbs free energy in mks units of * J kmol-1 K-1. */ -double WaterPropsIAPWS:: -Gibbs(double temperature, double rho) { +double WaterPropsIAPWS::Gibbs(double temperature, double rho) { setState(temperature, rho); double gRT = m_phi->gibbs_RT(); return (gRT * Rgas * temperature); } -double WaterPropsIAPWS:: -Gibbs() const { + +double WaterPropsIAPWS::Gibbs() const { double gRT = m_phi->gibbs_RT(); double temperature = T_c/tau; return (gRT * Rgas * temperature); @@ -437,14 +405,12 @@ 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(double temperature, double rho) { calcDim(temperature, rho); m_phi->tdpolycalc(tau, delta); } - /* * Calculate the enthalpy in mks units of * J kmol-1 K-1. @@ -495,25 +461,21 @@ entropy(double temperature, double rho) { * Calculate the enthalpy in mks units of * J kmol-1 K-1. */ -double WaterPropsIAPWS:: -entropy() const { +double WaterPropsIAPWS::entropy() const { double sR = m_phi->entropy_R(); return (sR * Rgas); } - /* * Calculate heat capacity at constant volume * J kmol-1 K-1. */ -double WaterPropsIAPWS:: -cv(double temperature, double rho) { +double WaterPropsIAPWS::cv(double temperature, double rho) { setState(temperature, rho); double cvR = m_phi->cv_R(); return (cvR * Rgas); } - /* * Calculate heat capacity at constant pressure * J kmol-1 K-1. diff --git a/Cantera/src/thermo/WaterPropsIAPWS.h b/Cantera/src/thermo/WaterPropsIAPWS.h index 92496b9d5..bfed57800 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.h +++ b/Cantera/src/thermo/WaterPropsIAPWS.h @@ -169,7 +169,6 @@ public: //! using the last temperature and density double helmholtzFE() const; - //! Calculate the Gibbs free energy in mks units of J kmol-1 K-1. /*! * @param temperature temperature (kelvin) @@ -312,20 +311,17 @@ public: */ double coeffThermExp(double temperature, double pressure); - - //! Returns the coefficient of isothermal compressibility as a function - //! of temperature and pressure. + //! Returns the coefficient of isothermal compressibility for the + //! state of the object /*! * kappa = - d (ln V) / dP at constant T. * * units - 1/Pascal * - * @param temperature Input temperature (Kelvin) - * @param pressure Input pressure (Pa) * @return * returns the isothermal compressibility */ - double isothermalCompressibility(double temperature, double pressure); + double isothermalCompressibility() const; //! This function returns an estimated value for the saturation pressure. /*! diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp index 1cfb7b7fb..0e5ee73b3 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp @@ -646,7 +646,7 @@ double WaterPropsIAPWSphi::pressureM_rhoRT(double tau, double delta) { } /* - * Calculate d_phiR_d(delta), the second derivative of phiR + * Calculate d2_phiR_dd(delta), the second derivative of phiR * wrt delta * * tau = dimensionless temperature @@ -752,7 +752,7 @@ double WaterPropsIAPWSphi::phiR_dd() const { } /* - * Calculate d_phi0_d(delta), the first derivative of phi0 + * Calculate d2_phi0_dd(delta), the second derivative of phi0 * wrt delta * * tau = dimensionless temperature @@ -764,7 +764,7 @@ double WaterPropsIAPWSphi::phi0_dd() const { } /* - * Calculate the dPhidDelta function, which is basically the derivative + * Calculate the d2_PhidDelta2 function, which is the second derivative * of helmholtz free energy wrt delta * Eqn. (6.4) */ @@ -776,6 +776,16 @@ double WaterPropsIAPWSphi::phi_dd(double tau, double delta) { return retn; } +double WaterPropsIAPWSphi::dimdpdrho(double tau, double delta) { + tdpolycalc(tau, delta); + double res1 = phiR_d(); + double res2 = phiR_dd(); + double retn = 1.0 + delta * (2.0*res1 + delta*res2); + return retn; +} + + + /* * Calculate d_phi0/d(tau) */ diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.h b/Cantera/src/thermo/WaterPropsIAPWSphi.h index a1110c768..a2a11f8ce 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.h +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.h @@ -98,6 +98,16 @@ public: */ double pressureM_rhoRT(double tau, double delta); + //! Dimensionless derivative of p wrt rho at constant T + /*! + * dp/drho * 1/RT = (2. delta phi_d() + delta**2 phi_dd()) + * (1.0 + 2. delta phiR_d() + delta**2 phiR_dd()) + * + * @param tau Dimensionless temperature = T_c/T + * @param delta Dimensionless density = delta = rho / Rho_c + */ + double dimdpdrho(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.