Fix watertransport

include/cantera/thermo/WaterProps.h

@@ -236,40 +236,6 @@ class WaterProps
      */
     doublereal isothermalCompressibility_IAPWS(doublereal T, doublereal 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.
-     */
-    doublereal viscosityWater() const;
-
-    //! Returns the thermal conductivity of water at the current conditions
-    //! (W/m/K)
-    /*!
-     * This function calculates the value of the thermal conductivity of
-     * 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.
-     */
-    doublereal thermalConductivityWater() const;
-
 protected:
     //! Pointer to the WaterPropsIAPWS object
     WaterPropsIAPWS* m_waterIAPWS;

include/cantera/thermo/WaterSSTP.h

@@ -65,9 +65,14 @@ class WaterProps;
  *
  * ## Instantiation of the Class
  *
- * The constructor for this phase is NOT located in the default ThermoFactory
- * for %Cantera. However, a new WaterSSTP object may be created by the following
- * code snippets, combined with an XML file given in the XML example section.
+ * A new WaterSSTP object may be created by the following code snippets,
+ * combined with an XML file given in the XML example section.
+ *
+ * @code
+ *      ThermoPhase* w = newPhase("waterSSTPphase.xml");
+ * @endcode
+ *
+ * or
  *
  * @code
  *      WaterSSTP *w = new WaterSSTP("waterSSTPphase.xml","");

include/cantera/transport/WaterTransport.h

@@ -81,11 +81,9 @@ class WaterTransport : public Transport
      */
     virtual doublereal thermalConductivity();
 
-private:
-    //! Routine to do some common initializations at the start of using
-    //! this routine.
-    void initTP();
+    virtual void init(thermo_t* thermo, int mode=0, int log_level=0);
 
+private:
     //! Pointer to the WaterPropsIAPWS object, which does the actual calculations
     //! for the real equation of state
     /*!

interfaces/cython/cantera/liquidvapor.py

@@ -1,11 +1,19 @@
 # This file is part of Cantera. See License.txt in the top-level directory or
 # at http://www.cantera.org/license.txt for license and copyright information.
 
-from . import PureFluid
+from . import PureFluid, _cantera
+
+
 
 def Water():
-    """Create a `PureFluid` object using the equation of state for water."""
-    return PureFluid('liquidvapor.xml','water')
+    """
+    Create a `PureFluid` object using the equation of state for water and the
+    `WaterTransport` class for viscosity and thermal conductivity."""
+    class WaterWithTransport(PureFluid, _cantera.Transport):
+        __slots__ = ()
+
+
+    return WaterWithTransport('liquidvapor.xml', 'water', transport_model='Water')
 
 def Nitrogen():
     """Create a `PureFluid` object using the equation of state for nitrogen."""

interfaces/cython/cantera/test/test_transport.py

@@ -268,6 +268,59 @@ def test_molar_fluxes(self):
         # self.assertNear(sum(fluxes1) / sum(abs(fluxes1)), 0.0)
 
 
+class TestWaterTransport(utilities.CanteraTest):
+    """
+    Comparison values are taken from the NIST Chemistry WebBook. Agreement is
+    limited by the use of a different equation of state here (Reynolds) than
+    in the Webbook (IAPWS95), as well as a different set of coefficients for
+    the transport property model. Differences are largest in the region near
+    the critical point.
+    """
+    @classmethod
+    def setUpClass(cls):
+        cls.water = ct.Water()
+
+    def check_viscosity(self, T, P, mu, rtol):
+        self.water.TP = T, P
+        self.assertNear(self.water.viscosity, mu, rtol)
+
+    def check_thermal_conductivity(self, T, P, k, rtol):
+        self.water.TP = T, P
+        self.assertNear(self.water.thermal_conductivity, k, rtol)
+
+    def test_viscosity_liquid(self):
+        self.check_viscosity(400, 1e6, 2.1880e-4, 1e-3)
+        self.check_viscosity(400, 8e6, 2.2061e-4, 1e-3)
+        self.check_viscosity(620, 1.6e7, 6.7489e-5, 2e-3)
+        self.check_viscosity(620, 2.8e7, 7.5684e-5, 2e-3)
+
+    def test_thermal_conductivity_liquid(self):
+        self.check_thermal_conductivity(400, 1e6, 0.68410, 1e-3)
+        self.check_thermal_conductivity(400, 8e6, 0.68836, 1e-3)
+        self.check_thermal_conductivity(620, 1.6e7, 0.45458, 2e-3)
+        self.check_thermal_conductivity(620, 2.8e7, 0.49705, 2e-3)
+
+    def test_viscosity_vapor(self):
+        self.check_viscosity(600, 1e6, 2.1329e-5, 1e-3)
+        self.check_viscosity(620, 5e6, 2.1983e-5, 1e-3)
+        self.check_viscosity(620, 1.5e7, 2.2858e-5, 2e-3)
+
+    def test_thermal_conductivity_vapor(self):
+        self.check_thermal_conductivity(600, 1e6, 0.047636, 1e-3)
+        self.check_thermal_conductivity(620, 5e6, 0.055781, 1e-3)
+        self.check_thermal_conductivity(620, 1.5e7, 0.10524, 2e-3)
+
+    def test_viscosity_supercritical(self):
+        self.check_viscosity(660, 2.2e7, 2.7129e-5, 2e-3)
+        self.check_viscosity(660, 2.54e7, 3.8212e-5, 1e-2)
+        self.check_viscosity(660, 2.8e7, 5.3159e-5, 1e-2)
+
+    def test_thermal_conductivity_supercritical(self):
+        self.check_thermal_conductivity(660, 2.2e7, 0.14872, 1e-2)
+        self.check_thermal_conductivity(660, 2.54e7, 0.35484, 2e-2)
+        self.check_thermal_conductivity(660, 2.8e7, 0.38479, 1e-2)
+
+
 class TestTransportData(utilities.CanteraTest):
     @classmethod
     def setUpClass(cls):

src/thermo/ThermoFactory.cpp

@@ -39,6 +39,7 @@
 #include "cantera/thermo/IdealMolalSoln.h"
 #include "cantera/thermo/MolarityIonicVPSSTP.h"
 #include "cantera/thermo/IdealSolnGasVPSS.h"
+#include "cantera/thermo/WaterSSTP.h"
 #include "cantera/base/stringUtils.h"
 
 using namespace std;
@@ -77,6 +78,7 @@ ThermoFactory::ThermoFactory()
     reg("RedlichKwong", []() { return new RedlichKwongMFTP(); });
     m_synonyms["RedlichKwongMFTP"] = "RedlichKwong";
     reg("MaskellSolidSolnPhase", []() { return new MaskellSolidSolnPhase(); });
+    reg("PureLiquidWater", []() { return new WaterSSTP(); });
 }
 
 ThermoPhase* ThermoFactory::newThermoPhase(const std::string& model)

src/thermo/WaterProps.cpp

@@ -263,153 +263,4 @@ doublereal WaterProps::isothermalCompressibility_IAPWS(doublereal temp, doublere
     return m_waterIAPWS->isothermalCompressibility();
 }
 
-static const doublereal H[4] = {1., 0.978197, 0.579829, -0.202354};
-
-static const doublereal 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.},
-};
-
-static const doublereal rhoStar = 317.763; // kg / m3
-static const doublereal presStar = 22.115E6; // Pa
-
-doublereal WaterProps::viscosityWater() const
-{
-    static const doublereal TStar = 647.27; // Kelvin
-    static const doublereal muStar = 55.071E-6; //Pa s
-    doublereal temp = m_waterIAPWS->temperature();
-    doublereal dens = m_waterIAPWS->density();
-
-    doublereal rhobar = dens/rhoStar;
-    doublereal tbar = temp / TStar;
-    doublereal tbar2 = tbar * tbar;
-    doublereal tbar3 = tbar2 * tbar;
-    doublereal mu0bar = std::sqrt(tbar) / (H[0] + H[1]/tbar + H[2]/tbar2 + H[3]/tbar3);
-
-    doublereal tfac1 = 1.0 / tbar - 1.0;
-    doublereal tfac2 = tfac1 * tfac1;
-    doublereal tfac3 = tfac2 * tfac1;
-    doublereal tfac4 = tfac3 * tfac1;
-    doublereal tfac5 = tfac4 * tfac1;
-
-    doublereal rfac1 = rhobar - 1.0;
-    doublereal rfac2 = rfac1 * rfac1;
-    doublereal rfac3 = rfac2 * rfac1;
-    doublereal rfac4 = rfac3 * rfac1;
-    doublereal rfac5 = rfac4 * rfac1;
-    doublereal rfac6 = rfac5 * rfac1;
-
-    doublereal 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
-                     );
-    doublereal mu1bar = std::exp(rhobar * sum);
-
-    // Apply the near-critical point corrections if necessary
-    doublereal mu2bar = 1.0;
-    if (tbar >= 0.9970 && tbar <= 1.0082 && rhobar >= 0.755 && rhobar <= 1.290) {
-        doublereal drhodp = 1.0 / m_waterIAPWS->dpdrho();
-        drhodp *= presStar / rhoStar;
-        doublereal xsi = rhobar * drhodp;
-        if (xsi >= 21.93) {
-            mu2bar = 0.922 * std::pow(xsi, 0.0263);
-        }
-    }
-
-    doublereal mubar = mu0bar * mu1bar * mu2bar;
-    return mubar * muStar;
-}
-
-doublereal WaterProps::thermalConductivityWater() const
-{
-    static const doublereal Tstar = 647.27;
-    static const doublereal rhostar = 317.763;
-    static const doublereal lambdastar = 0.4945;
-    static const doublereal presstar = 22.115E6;
-    static const doublereal L[4] = {
-        1.0000,
-        6.978267,
-        2.599096,
-        -0.998254
-    };
-    static const doublereal Lji[6][5] = {
-        { 1.3293046,    1.7018363,   5.2246158,   8.7127675, -1.8525999},
-        {-0.40452437,  -2.2156845, -10.124111,   -9.5000611,  0.93404690},
-        { 0.24409490,   1.6511057,   4.9874687,   4.3786606,  0.0},
-        { 0.018660751, -0.76736002, -0.27297694, -0.91783782, 0.0},
-        {-0.12961068,   0.37283344, -0.43083393,  0.0,        0.0},
-        { 0.044809953, -0.11203160,  0.13333849,  0.0,        0.0},
-    };
-
-    doublereal temp = m_waterIAPWS->temperature();
-    doublereal dens = m_waterIAPWS->density();
-
-    doublereal rhobar = dens/rhostar;
-    doublereal tbar = temp / Tstar;
-    doublereal tbar2 = tbar * tbar;
-    doublereal tbar3 = tbar2 * tbar;
-    doublereal lambda0bar = sqrt(tbar) / (L[0] + L[1]/tbar + L[2]/tbar2 + L[3]/tbar3);
-
-    doublereal tfac1 = 1.0 / tbar - 1.0;
-    doublereal tfac2 = tfac1 * tfac1;
-    doublereal tfac3 = tfac2 * tfac1;
-    doublereal tfac4 = tfac3 * tfac1;
-
-    doublereal rfac1 = rhobar - 1.0;
-    doublereal rfac2 = rfac1 * rfac1;
-    doublereal rfac3 = rfac2 * rfac1;
-    doublereal rfac4 = rfac3 * rfac1;
-    doublereal rfac5 = rfac4 * rfac1;
-
-    doublereal sum = (Lji[0][0]       + Lji[0][1]*tfac1        + Lji[0][2]*tfac2       + Lji[0][3]*tfac3       + Lji[0][4]*tfac4       +
-                      Lji[1][0]*rfac1 + Lji[1][1]*tfac1*rfac1  + Lji[1][2]*tfac2*rfac1 + Lji[1][3]*tfac3*rfac1 + Lji[1][4]*tfac4*rfac1 +
-                      Lji[2][0]*rfac2 + Lji[2][1]*tfac1*rfac2  + Lji[2][2]*tfac2*rfac2 + Lji[2][3]*tfac3*rfac2 +
-                      Lji[3][0]*rfac3 + Lji[3][1]*tfac1*rfac3  + Lji[3][2]*tfac2*rfac3 + Lji[3][3]*tfac3*rfac3 +
-                      Lji[4][0]*rfac4 + Lji[4][1]*tfac1*rfac4  + Lji[4][2]*tfac2*rfac4 +
-                      Lji[5][0]*rfac5 + Lji[5][1]*tfac1*rfac5  + Lji[5][2]*tfac2*rfac5
-                     );
-    doublereal lambda1bar = exp(rhobar * sum);
-    doublereal mu0bar = std::sqrt(tbar) / (H[0] + H[1]/tbar + H[2]/tbar2 + H[3]/tbar3);
-    doublereal tfac5 = tfac4 * tfac1;
-    doublereal rfac6 = rfac5 * rfac1;
-
-    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
-          );
-    doublereal mu1bar = std::exp(rhobar * sum);
-    doublereal t2r2 = tbar * tbar / (rhobar * rhobar);
-    doublereal drhodp = 1.0 / m_waterIAPWS->dpdrho();
-    drhodp *= presStar / rhoStar;
-    doublereal xsi = rhobar * drhodp;
-    doublereal xsipow = std::pow(xsi, 0.4678);
-    doublereal rho1 = rhobar - 1.;
-    doublereal rho2 = rho1 * rho1;
-    doublereal rho4 = rho2 * rho2;
-    doublereal temp2 = (tbar - 1.0) * (tbar - 1.0);
-
-    //     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 beta = m_waterIAPWS->coeffPresExp();
-    doublereal dpdT_const_rho = beta * GasConstant * dens / 18.015268;
-    dpdT_const_rho *= Tstar / presstar;
-    doublereal lambda2bar = 0.0013848 / (mu0bar * mu1bar) * t2r2 * dpdT_const_rho * dpdT_const_rho *
-                             xsipow * sqrt(rhobar) * exp(-18.66*temp2 - rho4);
-    return (lambda0bar * lambda1bar + lambda2bar) * lambdastar;
-}
-
 }

src/transport/TransportFactory.cpp

@@ -8,6 +8,7 @@
 #include "cantera/transport/MixTransport.h"
 #include "cantera/transport/UnityLewisTransport.h"
 #include "cantera/transport/IonGasTransport.h"
+#include "cantera/transport/WaterTransport.h"
 #include "cantera/transport/SolidTransport.h"
 #include "cantera/transport/DustyGasTransport.h"
 #include "cantera/transport/SimpleTransport.h"
@@ -52,6 +53,7 @@ TransportFactory::TransportFactory()
     reg("Mix", []() { return new MixTransport(); });
     reg("Multi", []() { return new MultiTransport(); });
     reg("Ion", []() { return new IonGasTransport(); });
+    reg("Water", []() { return new WaterTransport(); });
     m_synonyms["CK_Mix"] = "Mix";
     m_synonyms["CK_Multi"] = "Multi";
     reg("HighP", []() { return new HighPressureGasTransport(); });