Work on heat conduction and thermal conductivity (FS #172)

- added somerton.hh in fluidmatrixinteractions/2p which provides a model
  for the effective thermal conductivity (static class)
- removed *matrixHeatFlux() methods in waterairspatialparams, introduced
  thermalConductivitySolid() instead
- added property ThermalConductivityModel, which is set as default to
  Somerton
- added effThermalConductivity and calculateEffThermalConductivity
  methods, thermal conductivities are averaged harmonically

Reviewed by Melanie


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@10049 2fb0f335-1f38-0410-981e-8018bf24f1b0

dumux/implicit/2p2cni/2p2cnifluxvariables.hh

@@ -54,6 +54,7 @@ class TwoPTwoCNIFluxVariables : public TwoPTwoCFluxVariables<TypeTag>
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+    typedef typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel) ThermalConductivityModel;
 
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
     typedef typename GridView::template Codim<0>::Entity Element;
@@ -91,8 +92,17 @@ public:
     Scalar normalMatrixHeatFlux() const
     { return normalMatrixHeatFlux_; }
 
+    /*!
+     * \brief The local temperature gradient at the IP of the considered scv face.
+     */
     DimVector temperatureGradient() const
     { return temperatureGrad_; }
+    
+    /*!
+     * \brief The harmonically averaged effective thermal conductivity.
+     */
+    Scalar effThermalConductivity() const
+    { return lambdaEff_; }
 
 protected:
     void calculateValues_(const Problem &problem,
@@ -114,28 +124,37 @@ protected:
             temperatureGrad_ += tmp;
         }
 
-        // The spatial parameters calculates the actual heat flux vector
-        if (this->face().i != this->face().j)
-            problem.spatialParams().matrixHeatFlux(tmp,
-                                                   *this,
-                                                   elemVolVars,
-                                                   temperatureGrad_,
-                                                   element,
-                                                   this->fvGeometry_,
-                                                   faceIdx_);
-        else // heat flux at outflow boundaries
-            problem.spatialParams().boundaryMatrixHeatFlux(tmp,
-                                                           *this,
-                                                           elemVolVars,
-                                                           this->face(),
-                                                           element,
-                                                           this->fvGeometry_);
+        lambdaEff_ = 0;
+        calculateEffThermalConductivity_(problem, element, elemVolVars);
 
         // project the heat flux vector on the face's normal vector
-        normalMatrixHeatFlux_ = tmp*this->face().normal;
+        normalMatrixHeatFlux_ = temperatureGrad_*
+                                this->face().normal;
+        normalMatrixHeatFlux_ *= -lambdaEff_;
+    }
+
+    void calculateEffThermalConductivity_(const Problem &problem,
+                                        const Element &element,
+                                        const ElementVolumeVariables &elemVolVars)
+    {
+        const Scalar lambdaI =
+                ThermalConductivityModel::effectiveThermalConductivity(element,
+                                                                       elemVolVars,
+                                                                       this->fvGeometry_,
+                                                                       problem.spatialParams(),
+                                                                       this->face().i);
+        const Scalar lambdaJ =
+                ThermalConductivityModel::effectiveThermalConductivity(element,
+                                                                       elemVolVars,
+                                                                       this->fvGeometry_,
+                                                                       problem.spatialParams(),
+                                                                       this->face().j);
+        // -> harmonic mean
+        lambdaEff_ = harmonicMean(lambdaI, lambdaJ);
     }
 
 private:
+    Scalar lambdaEff_;
     Scalar normalMatrixHeatFlux_;
     DimVector temperatureGrad_;
     int faceIdx_;

dumux/implicit/2p2cni/2p2cniproperties.hh

@@ -43,6 +43,12 @@ namespace Properties
 NEW_TYPE_TAG(TwoPTwoCNI, INHERITS_FROM(TwoPTwoC));
 NEW_TYPE_TAG(BoxTwoPTwoCNI, INHERITS_FROM(BoxModel, TwoPTwoCNI));
 NEW_TYPE_TAG(CCTwoPTwoCNI, INHERITS_FROM(CCModel, TwoPTwoCNI));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(ThermalConductivityModel);   //!< The model for the effective thermal conductivity
 }
 }
 

dumux/implicit/2p2cni/2p2cnipropertydefaults.hh

@@ -36,6 +36,8 @@
 #include "2p2cnivolumevariables.hh"
 #include "2p2cnifluxvariables.hh"
 
+#include <dumux/material/fluidmatrixinteractions/2p/somerton.hh>
+
 namespace Dumux
 {
 
@@ -70,6 +72,9 @@ SET_PROP(TwoPTwoCNI, Indices)
     typedef TwoPTwoCNIIndices<TypeTag, formulation, 0> type;
 };
 
+//! Somerton is used as default model to compute the effective thermal heat conductivity
+SET_TYPE_PROP(TwoPTwoCNI, ThermalConductivityModel, Somerton<TypeTag>);
+
 }
 
 }

dumux/material/fluidmatrixinteractions/2p/somerton.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief   Relation for the saturation-dependent effective thermal conductivity
+ */
+#ifndef SOMERTON_HH
+#define SOMERTON_HH
+
+#include <algorithm>
+
+namespace Dumux
+{
+/*!
+ * \ingroup fluidmatrixinteractionslaws
+ *
+ * \brief Relation for the saturation-dependent effective thermal conductivity
+ *
+ *  The Somerton method computes the thermal conductivity of dry and the wet soil material
+ *  and uses a root function of the wetting saturation to compute the
+ *  effective thermal conductivity for a two-phase fluidsystem.
+ */
+template<class TypeTag>
+class Somerton
+{
+public:
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    enum { wPhaseIdx = Indices::wPhaseIdx };
+
+
+    /*!
+     * \brief Returns the effective thermal conductivity \f$[W/m^2]\f$ after Somerton (1974).
+     *
+     * The material law is:
+     * \f[
+     l_eff = l_solid + (l_wet - l_solid)
+     \f]
+     *
+     * \param element The finite element
+     * \param elemVolVars The volume variables on the element
+     * \param fvGeometry The finite volume geometry
+     * \param spatialParams The spatial parameters
+     * \param scvIdx The local index of the sub-control volume where
+     *                    the effective thermal conductivity is computed
+     *
+     * \return Effective thermal conductivity \f$[W/m^2]\f$ after Somerton (1974)
+     */
+    static Scalar effectiveThermalConductivity(const Element &element,
+                                               const ElementVolumeVariables &elemVolVars,
+                                               const FVElementGeometry &fvGeometry,
+                                               const SpatialParams &spatialParams,
+                                               const int scvIdx)
+    {
+        const Scalar lambdaSolid = spatialParams.thermalConductivitySolid(element, fvGeometry, scvIdx);
+        const Scalar poro = spatialParams.porosity(element, fvGeometry, scvIdx);
+
+        const Scalar Sw = std::max<Scalar>(0.0, elemVolVars[scvIdx].saturation(wPhaseIdx));
+        const Scalar lWater = elemVolVars[scvIdx].thermalConductivity(wPhaseIdx);
+
+        const Scalar lsat = pow(lambdaSolid, (1-poro)) * pow(lWater, poro);
+        const Scalar ldry = pow(lambdaSolid, (1-poro));
+
+        return ldry + sqrt(Sw) * (ldry - lsat);
+    }
+};
+}
+#endif

test/implicit/2p2cni/waterairspatialparams.hh

@@ -79,11 +79,6 @@ class WaterAirSpatialParams : public ImplicitSpatialParams<TypeTag>
         dimWorld=GridView::dimensionworld
     };
 
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-    enum {
-        wPhaseIdx = Indices::wPhaseIdx
-    };
-
     typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
     typedef Dune::FieldVector<CoordScalar,dimWorld> DimVector;
 
@@ -114,6 +109,9 @@ public:
         // porosities
         finePorosity_ = 0.3;
         coarsePorosity_ = 0.3;
+        
+        // heat conductivity of granite
+        lambdaSolid_ = 2.8;
 
         // residual saturations
         fineMaterialParams_.setSwr(0.2);
@@ -157,7 +155,7 @@ public:
      * \param scvIdx The local index of the sub-control volume where
      *                    the porosity needs to be defined
      */
-    double porosity(const Element &element,
+    Scalar porosity(const Element &element,
                     const FVElementGeometry &fvGeometry,
                     const int scvIdx) const
     {
@@ -198,7 +196,7 @@ public:
      * \param scvIdx The local index of the sub-control volume where
      *                    the heat capacity needs to be defined
      */
-    double heatCapacity(const Element &element,
+    Scalar heatCapacity(const Element &element,
                         const FVElementGeometry &fvGeometry,
                         const int scvIdx) const
     {
@@ -209,50 +207,15 @@ public:
     }
 
     /*!
-     * \brief Calculate the heat flux \f$[W/m^2]\f$ through the
-     *        rock matrix based on the temperature gradient \f$[K / m]\f$
-     *
-     * This is only required for non-isothermal models.
+     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
      *
-     * \param heatFlux The resulting heat flux vector
-     * \param fluxVars The flux variables
-     * \param elemVolVars The volume variables
-     * \param temperatureGrad The temperature gradient
-     * \param element The current finite element
-     * \param fvGeometry The finite volume geometry of the current element
-     * \param faceIdx The local index of the sub-control volume face where
-     *                    the matrix heat flux should be calculated
+     * \param pos The global position
      */
-    void matrixHeatFlux(DimVector &heatFlux,
-                        const FluxVariables &fluxVars,
-                        const ElementVolumeVariables &elemVolVars,
-                        const DimVector &temperatureGrad,
-                        const Element &element,
-                        const FVElementGeometry &fvGeometry,
-                        const int faceIdx) const
+    Scalar thermalConductivitySolid(const Element &element,
+                                    const FVElementGeometry &fvGeometry,
+                                    const int scvIdx) const
     {
-        static const Scalar lWater = 0.6;
-        static const Scalar lGranite = 2.8;
-
-        // arithmetic mean of the liquid saturation and the porosity
-        const int i = fvGeometry.subContVolFace[faceIdx].i;
-        const int j = fvGeometry.subContVolFace[faceIdx].j;
-        Scalar Sl = std::max<Scalar>(0.0, (elemVolVars[i].saturation(wPhaseIdx) +
-                                           elemVolVars[j].saturation(wPhaseIdx)) / 2);
-        Scalar poro = (porosity(element, fvGeometry, i) +
-                       porosity(element, fvGeometry, j)) / 2;
-
-        Scalar lsat = pow(lGranite, (1-poro)) * pow(lWater, poro);
-        Scalar ldry = pow(lGranite, (1-poro));
-
-        // the heat conductivity of the matrix. in general this is a
-        // tensorial value, but we assume isotropic heat conductivity.
-        Scalar heatCond = ldry + sqrt(Sl) * (ldry - lsat);
-
-        // the matrix heat flux is the negative temperature gradient
-        // times the heat conductivity.
-        heatFlux = temperatureGrad;
-        heatFlux *= -heatCond;
+        return lambdaSolid_;
     }
 
 private:
@@ -266,6 +229,9 @@ private:
     Scalar finePorosity_;
     Scalar coarsePorosity_;
 
+    // heat conductivity of the solid material only
+    Scalar lambdaSolid_;
+
     MaterialLawParams fineMaterialParams_;
     MaterialLawParams coarseMaterialParams_;
 };