From d33490f16b932c6b3d63b5d7a54a852fde0a6687 Mon Sep 17 00:00:00 2001
From: Thomas Fetzer <thomas.fetzer@iws.uni-stuttgart.de>
Date: Thu, 16 Nov 2017 12:14:55 +0100
Subject: [PATCH] [components] Introduce heat capacity functions, use them to
 calculate enthalpy

Nitrogen already had a heat capacity function, however, it was calculated differently for the other gases.
---
 dumux/material/components/ch4.hh            | 33 ++++++++++-------
 dumux/material/components/h2.hh             | 29 +++++++++------
 dumux/material/components/n2.hh             | 40 +++++----------------
 dumux/material/components/o2.hh             | 31 ++++++++++------
 dumux/material/fluidsystems/h2on2kinetic.hh |  2 +-
 5 files changed, 70 insertions(+), 65 deletions(-)

diff --git a/dumux/material/components/ch4.hh b/dumux/material/components/ch4.hh
index fe095eec15..65c7be37ff 100644
--- a/dumux/material/components/ch4.hh
+++ b/dumux/material/components/ch4.hh
@@ -132,13 +132,26 @@ public:
     /*!
      * \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of pure methane gas.
      *
-     * \param T temperature of component in \f$\mathrm{[K]}\f$
+     * \param temperature temperature of component in \f$\mathrm{[K]}\f$
      * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
-     *
-     * See: R. Reid, et al. (1987, pp 154, 657, 671) \cite reid1987
      */
-    static const Scalar gasEnthalpy(Scalar T,
+    static const Scalar gasEnthalpy(Scalar temperature,
                                     Scalar pressure)
+    {
+        return gasHeatCapacity(temperature, pressure) * temperature;
+    }
+
+    /*!
+     * \brief Specific isobaric heat capacity \f$\mathrm{[J/(kg*K)]}\f$ of pure
+     *        methane gas.
+     *
+     * This is equivalent to the partial derivative of the specific
+     * enthalpy to the temperature.
+     *
+     * See: R. Reid, et al. (1987, pp 154, 657, 665) \cite reid1987
+     */
+    static Scalar gasHeatCapacity(Scalar T,
+                                  Scalar pressure)
     {
         // method of Joback
         const Scalar cpVapA = 19.25;
@@ -146,16 +159,12 @@ public:
         const Scalar cpVapC = 1.197e-5;
         const Scalar cpVapD = -1.132e-8;
 
-        //Scalar cp =
-        //    cpVapA + T*(cpVapB + T*(cpVapC + T*cpVapD));
-
-        // calculate: \int_0^T c_p dT
         return
-            1/molarMass()* // conversion from [J/mol] to [J/kg]
-            T*(cpVapA + T*
-               (cpVapB/2 + T*
+            1/molarMass()* // conversion from [J/(mol*K)] to [J/(kg*K)]
+            (cpVapA + T*
+              (cpVapB/2 + T*
                 (cpVapC/3 + T*
-                 (cpVapD/4))));
+                  (cpVapD/4))));
     }
 
     /*!
diff --git a/dumux/material/components/h2.hh b/dumux/material/components/h2.hh
index 65af6a93c0..32e7d43b9b 100644
--- a/dumux/material/components/h2.hh
+++ b/dumux/material/components/h2.hh
@@ -146,13 +146,26 @@ public:
     /*!
      * \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of pure hydrogen gas.
      *
-     * \param T temperature of component in \f$\mathrm{[K]}\f$
+     * \param temperature temperature of component in \f$\mathrm{[K]}\f$
      * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+     */
+    static const Scalar gasEnthalpy(Scalar temperature,
+                                    Scalar pressure)
+    {
+        return gasHeatCapacity(temperature, pressure) * temperature;
+    }
+
+    /*!
+     * \brief Specific isobaric heat capacity \f$\mathrm{[J/(kg*K)]}\f$ of pure
+     *        hydrogen gas.
+     *
+     * This is equivalent to the partial derivative of the specific
+     * enthalpy to the temperature.
      *
      * See: R. Reid, et al. (1987, pp 154, 657, 665) \cite reid1987
      */
-    static const Scalar gasEnthalpy(Scalar T,
-                                    Scalar pressure)
+    static const Scalar gasHeatCapacity(Scalar T,
+                                        Scalar pressure)
     {
         // method of Joback
         const Scalar cpVapA = 27.14;
@@ -160,14 +173,10 @@ public:
         const Scalar cpVapC = -1.381e-5;
         const Scalar cpVapD = 7.645e-9;
 
-        //Scalar cp =
-        //    cpVapA + T*(cpVapB + T*(cpVapC + T*cpVapD));
-
-        // calculate: \int_0^T c_p dT
         return
-            1/molarMass()* // conversion from [J/mol] to [J/kg]
-            T*(cpVapA + T*
-               (cpVapB/2 + T*
+            1/molarMass()* // conversion from [J/(mol*K)] to [J/(kg*K)]
+            (cpVapA + T*
+              (cpVapB/2 + T*
                 (cpVapC/3 + T*
                  (cpVapD/4))));
     }
diff --git a/dumux/material/components/n2.hh b/dumux/material/components/n2.hh
index b9da95dce7..b0d5724008 100644
--- a/dumux/material/components/n2.hh
+++ b/dumux/material/components/n2.hh
@@ -160,36 +160,13 @@ public:
     /*!
      * \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of pure nitrogen gas.
      *
-     * \param T temperature of component in \f$\mathrm{[K]}\f$
+     * \param temperature temperature of component in \f$\mathrm{[K]}\f$
      * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
-     *
-     * See: R. Reid, et al. (1987, pp 154, 657, 665) \cite reid1987
      */
-    static const Scalar gasEnthalpy(Scalar T,
+    static const Scalar gasEnthalpy(Scalar temperature,
                                     Scalar pressure)
     {
-        // method of Joback
-        const Scalar cpVapA = 31.15;
-        const Scalar cpVapB = -0.01357;
-        const Scalar cpVapC = 2.680e-5;
-        const Scalar cpVapD = -1.168e-8;
-
-        // calculate: \int_0^T c_p dT
-        return
-            1/molarMass()* // conversion from [J/(mol K)] to [J/(kg K)]
-
-            T*(cpVapA + T*
-               (cpVapB/2 + T*
-                (cpVapC/3 + T*
-                 (cpVapD/4))));
-
-//#warning NIST DATA STUPID INTERPOLATION
-//        Scalar T2 = 300.;
-//        Scalar T1 = 285.;
-//        Scalar h2 = 311200.;
-//        Scalar h1 = 295580.;
-//        Scalar h = h1+ (h2-h1) / (T2-T1) * (T-T1);
-//        return h ;
+        return gasHeatCapacity(temperature, pressure) * temperature;
     }
 
     /*!
@@ -220,6 +197,8 @@ public:
      *
      * This is equivalent to the partial derivative of the specific
      * enthalpy to the temperature.
+     *
+     * See: R. Reid, et al. (1987, pp 154, 657, 665) \cite reid1987
      */
     static const Scalar gasHeatCapacity(Scalar T,
                                         Scalar pressure)
@@ -232,11 +211,10 @@ public:
 
         return
             1/molarMass()* // conversion from [J/(mol K)] to [J/(kg K)]
-
-            cpVapA + T*
-            (cpVapB + T*
-             (cpVapC + T*
-              (cpVapD)));
+            (cpVapA + T*
+              (cpVapB/2 + T*
+                (cpVapC/3 + T*
+                  (cpVapD/4))));
     }
 
     /*!
diff --git a/dumux/material/components/o2.hh b/dumux/material/components/o2.hh
index 0ecb4a4daa..6b7ab93f2b 100644
--- a/dumux/material/components/o2.hh
+++ b/dumux/material/components/o2.hh
@@ -158,13 +158,26 @@ public:
     /*!
      * \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of pure oxygen gas.
      *
-     * \param T temperature of component in \f$\mathrm{[K]}\f$
+     * \param temperature temperature of component in \f$\mathrm{[K]}\f$
      * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+     */
+    static Scalar gasEnthalpy(Scalar temperature,
+                              Scalar pressure)
+    {
+        return gasHeatCapacity(temperature, pressure) * temperature;
+    }
+
+    /*!
+     * \brief Specific isobaric heat capacity \f$\mathrm{[J/(kg*K)]}\f$ of pure
+     *        oxygen gas.
+     *
+     * This is equivalent to the partial derivative of the specific
+     * enthalpy to the temperature.
      *
      * See: R. Reid, et al. (1987, pp 154, 657, 665) \cite reid1987
      */
-    static Scalar gasEnthalpy(Scalar T,
-                                    Scalar pressure)
+    static Scalar gasHeatCapacity(Scalar T,
+                                  Scalar pressure)
     {
         // method of Joback
         const Scalar cpVapA = 28.11;
@@ -172,16 +185,12 @@ public:
         const Scalar cpVapC = 1.746e-5;
         const Scalar cpVapD = -1.065e-8;
 
-        //Scalar cp =
-        //    cpVapA + T*(cpVapB + T*(cpVapC + T*cpVapD));
-
-        // calculate: \int_0^T c_p dT
         return
-            1/molarMass()* // conversion from [J/mol] to [J/kg]
-            T*(cpVapA + T*
-               (cpVapB/2 + T*
+            1/molarMass()* // conversion from [J/(mol*K)] to [J/(kg*K)]
+            (cpVapA + T*
+              (cpVapB/2 + T*
                 (cpVapC/3 + T*
-                 (cpVapD/4))));
+                  (cpVapD/4))));
     }
 
     /*!
diff --git a/dumux/material/fluidsystems/h2on2kinetic.hh b/dumux/material/fluidsystems/h2on2kinetic.hh
index 3559b675e6..f635e7ab17 100644
--- a/dumux/material/fluidsystems/h2on2kinetic.hh
+++ b/dumux/material/fluidsystems/h2on2kinetic.hh
@@ -90,7 +90,7 @@ public:
                 case ParentType::H2OIdx:
                     return ParentType::H2O::liquidEnthalpy(T, p);
                 case ParentType::N2Idx:
-                    return ParentType::N2::gasEnthalpy(T, p);
+                    return ParentType::N2::gasEnthalpy(T, p); // TODO: should be liquid enthalpy
                 default:
                     DUNE_THROW(Dune::NotImplemented,
                                "wrong index");
-- 
GitLab