From 7bb80f7b9fb953da15c1f3abd1ab19d21e2194e5 Mon Sep 17 00:00:00 2001
From: Benjamin Faigle <benjamin.faigle@posteo.de>
Date: Wed, 12 Feb 2014 10:14:29 +0000
Subject: [PATCH] Add functions for the heat capacity to the components air and
 co2, that will be required by the decoupled NI-model. Enabled acess to these
 methods in the co2-brine-fluidsystem (only pure-phase heat capacities are
 used). reviewed by Klaus

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@12462 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
 dumux/material/components/air.hh              | 38 +++++++++++++++++++
 dumux/material/components/co2.hh              | 25 ++++++++++++
 dumux/material/components/h2o.hh              |  2 +-
 .../fluidsystems/brineco2fluidsystem.hh       | 17 +++++++--
 4 files changed, 78 insertions(+), 4 deletions(-)

diff --git a/dumux/material/components/air.hh b/dumux/material/components/air.hh
index 23600c5f03..90da3c6575 100644
--- a/dumux/material/components/air.hh
+++ b/dumux/material/components/air.hh
@@ -203,6 +203,44 @@ public:
             / molarMass(); // conversion from [J/(mol K)] to [J/(kg K)]
     }
 
+    /*!
+     * \brief Specific isobaric heat capacity \f$[J/(kg K)]\f$ of pure
+     *        air.
+     *
+     *  This methods uses the formula for "zero-pressure" heat capacity that
+     *  is only dependent on temperature, because the pressure dependence is rather small.
+     *  This one should be accurate for a pressure of 1 atm.
+     *  Values taken from NASA Contractor Report 4755, Real-Gas Flow Properties for NASA
+     *  Langley Research Center Aerothermodynamic Facilities Complex Wind Tunnels
+     *  using data from
+     *  Hilsenrath et al 1955, "Tables of Thermal Properties of Gases"
+     * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+     * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+     */
+    static const Scalar gasHeatCapacity(Scalar temperature,
+                                        Scalar pressure)
+    {
+        // scale temperature with referenence temp of 100K
+        Scalar phi = temperature/100;
+
+        Scalar c_p = 0.661738E+01
+                -0.105885E+01 * phi
+                +0.201650E+00 * pow(phi,2.)
+                -0.196930E-01 * pow(phi,3.)
+                +0.106460E-02 * pow(phi,4.)
+                -0.303284E-04 * pow(phi,5.)
+                +0.355861E-06 * pow(phi,6.);
+        c_p +=   -0.549169E+01 * pow(phi,-1.)
+                +0.585171E+01* pow(phi,-2.)
+                -0.372865E+01* pow(phi,-3.)
+                +0.133981E+01* pow(phi,-4.)
+                -0.233758E+00* pow(phi,-5.)
+                +0.125718E-01* pow(phi,-6.);
+        c_p *= IdealGas::R / (molarMass() * 1000); // in J/mol/K * mol / kg / 1000 = kJ/kg/K
+
+
+        return  c_p;
+    }
 };
 
 } // end namespace
diff --git a/dumux/material/components/co2.hh b/dumux/material/components/co2.hh
index b13a72a955..6813e41919 100644
--- a/dumux/material/components/co2.hh
+++ b/dumux/material/components/co2.hh
@@ -249,11 +249,36 @@ public:
         DUNE_THROW(NumericalProblem, "CO2::liquidPressure()");
     }
 
+    /*!
+     * \brief Specific isobaric heat capacity of the component [J/kg] as a liquid.
+     * USE WITH CAUTION! Exploits enthalpy function with artificial increment
+     * of the temperature!
+     * Equation with which the specific heat capacity is calculated : \f$ c_p = \frac{dh}{dT}\f$
+     * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+     * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+     */
+    static Scalar liquidHeatCapacity(Scalar temperature, Scalar pressure)
+    {
+        //temperature difference :
+        Scalar dT = 1.; // 1K temperature increment
+        Scalar temperature2 = temperature+dT;
+
+        // enthalpy difference
+        Scalar hold = liquidEnthalpy(temperature, pressure);
+        Scalar hnew = liquidEnthalpy(temperature2, pressure);
+        Scalar dh = hold-hnew;
+
+        //specific heat capacity
+        return dh/dT ;
+    }
+
 
     /*!
      * \brief The dynamic viscosity [N/m^3*s] of CO2.
      * Equations given in: - Vesovic et al., 1990
      *                     - Fenhour et al., 1998
+     * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+     * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
      */
     static Scalar gasViscosity(Scalar temperature, Scalar pressure)
     {
diff --git a/dumux/material/components/h2o.hh b/dumux/material/components/h2o.hh
index 83fe5cdfd3..29ea091236 100644
--- a/dumux/material/components/h2o.hh
+++ b/dumux/material/components/h2o.hh
@@ -449,7 +449,7 @@ public:
     }
 
     /*!
-     * \brief Specific isochoric heat capacity of liquid water \f$\mathrm{[J/kg]}\f$.
+     * \brief Specific isochoric heat capacity of liquid water \f$\mathrm{[J/m^3]}\f$.
      *
      * \param temperature temperature of component in \f$\mathrm{[K]}\f$
      * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
diff --git a/dumux/material/fluidsystems/brineco2fluidsystem.hh b/dumux/material/fluidsystems/brineco2fluidsystem.hh
index 63f528dbe3..3fe4b8f163 100644
--- a/dumux/material/fluidsystems/brineco2fluidsystem.hh
+++ b/dumux/material/fluidsystems/brineco2fluidsystem.hh
@@ -560,14 +560,25 @@ public:
 
     /*!
      * \copydoc BaseFluidSystem::heatCapacity
+     * 
+     * We employ the heat capacity of the pure phases.
+     * Todo: Include compositional effects.
+     *
+     * \param fluidState An arbitrary fluid state
+     * \param phaseIdx The index of the fluid phase to consider
+     * \tparam FluidState the fluid state class
      */
+    using Base::heatCapacity;
     template <class FluidState>
     static Scalar heatCapacity(const FluidState &fluidState,
-                               const ParameterCache &paramCache,
                                int phaseIdx)
     {
-        // TODO!
-        DUNE_THROW(Dune::NotImplemented, "Heat capacities");
+        if(phaseIdx == wPhaseIdx)
+            return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx),
+                                           fluidState.pressure(phaseIdx));
+        else
+            return CO2::liquidHeatCapacity(fluidState.temperature(phaseIdx),
+                                       fluidState.pressure(phaseIdx));
     }
 
 private:
-- 
GitLab