From f8e88f6d05bcef97a1e9b3606a7f64608dff8fd4 Mon Sep 17 00:00:00 2001
From: Klaus Mosthaf <klmos@env.dtu.dk>
Date: Mon, 21 Jan 2013 16:54:30 +0000
Subject: [PATCH] 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 | 53 +++++++----
dumux/implicit/2p2cni/2p2cniproperties.hh | 6 ++
.../implicit/2p2cni/2p2cnipropertydefaults.hh | 5 +
.../fluidmatrixinteractions/2p/somerton.hh | 91 +++++++++++++++++++
test/implicit/2p2cni/waterairspatialparams.hh | 62 +++----------
5 files changed, 152 insertions(+), 65 deletions(-)
create mode 100644 dumux/material/fluidmatrixinteractions/2p/somerton.hh
diff --git a/dumux/implicit/2p2cni/2p2cnifluxvariables.hh b/dumux/implicit/2p2cni/2p2cnifluxvariables.hh
index 499d7250a6..be4cda8133 100644
--- a/dumux/implicit/2p2cni/2p2cnifluxvariables.hh
+++ b/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_;
diff --git a/dumux/implicit/2p2cni/2p2cniproperties.hh b/dumux/implicit/2p2cni/2p2cniproperties.hh
index 92fb6fab82..059bc5315d 100644
--- a/dumux/implicit/2p2cni/2p2cniproperties.hh
+++ b/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
}
}
diff --git a/dumux/implicit/2p2cni/2p2cnipropertydefaults.hh b/dumux/implicit/2p2cni/2p2cnipropertydefaults.hh
index a9c7e3063f..dca073ed54 100644
--- a/dumux/implicit/2p2cni/2p2cnipropertydefaults.hh
+++ b/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>);
+
}
}
diff --git a/dumux/material/fluidmatrixinteractions/2p/somerton.hh b/dumux/material/fluidmatrixinteractions/2p/somerton.hh
new file mode 100644
index 0000000000..17e441c429
--- /dev/null
+++ b/dumux/material/fluidmatrixinteractions/2p/somerton.hh
@@ -0,0 +1,91 @@
+// -*- 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
diff --git a/test/implicit/2p2cni/waterairspatialparams.hh b/test/implicit/2p2cni/waterairspatialparams.hh
index 3961adc21c..df3b8e4dc0 100644
--- a/test/implicit/2p2cni/waterairspatialparams.hh
+++ b/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_;
};
--
GitLab