Made the isothermal compressibility function a real derivative function

instead of a numerical derivative.

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 <math.h>
 
@@ -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; 
   }
 

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;
   }
 

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.

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. 
   /*!

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)
  */

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.