diff --git a/dumux/implicit/2pni/2pnifluxvariables.hh b/dumux/implicit/2pni/2pnifluxvariables.hh
index 5683986a41bdb19477c27e3bbe99c95d47ae24c2..3c2ed188475262af9564ebfae1ef179749a7ea76 100644
--- a/dumux/implicit/2pni/2pnifluxvariables.hh
+++ b/dumux/implicit/2pni/2pnifluxvariables.hh
@@ -53,13 +53,14 @@ class TwoPNIFluxVariables : public ImplicitDarcyFluxVariables<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;
enum { dimWorld = GridView::dimensionworld };
typedef typename GridView::ctype CoordScalar;
- typedef Dune::FieldVector<CoordScalar, dimWorld> Vector;
+ typedef Dune::FieldVector<CoordScalar, dimWorld> DimVector;
public:
@@ -82,46 +83,88 @@ public:
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
: ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+ {
+ faceIdx_ = faceIdx;
+
+ calculateValues_(problem, element, elemVolVars);
+ }
+
+ /*!
+ * \brief The total heat flux \f$\mathrm{[J/s]}\f$ due to heat conduction
+ * of the rock matrix over the sub-control volume's face in
+ * direction of the face normal.
+ */
+ 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,
+ const Element &element,
+ const ElementVolumeVariables &elemVolVars)
{
// calculate temperature gradient using finite element
// gradients
- Vector temperatureGrad(0);
- for (int idx = 0; idx < fvGeometry.numFap; idx++)
+ temperatureGrad_ = 0;
+ DimVector tmp(0.0);
+ for (int idx = 0; idx < this->fvGeometry_.numFap; idx++)
{
- Vector feGrad = this->face().grad[idx];
+ tmp = this->face().grad[idx];
// index for the element volume variables
int volVarsIdx = this->face().fapIndices[idx];
-
- feGrad *= elemVolVars[volVarsIdx].temperature();
- temperatureGrad += feGrad;
+
+ tmp *= elemVolVars[volVarsIdx].temperature();
+ temperatureGrad_ += tmp;
}
- // The spatial parameters calculates the actual heat flux vector
- Vector heatFlux;
- problem.spatialParams().matrixHeatFlux(heatFlux,
- *this,
- elemVolVars,
- temperatureGrad,
- element,
- fvGeometry,
- faceIdx);
+ lambdaEff_ = 0;
+ calculateEffThermalConductivity_(problem, element, elemVolVars);
+
// project the heat flux vector on the face's normal vector
- normalMatrixHeatFlux_ = heatFlux * 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);
}
- /*!
- * \brief The total heat flux \f$\mathrm{[J/s]}\f$ due to heat conduction
- * of the rock matrix over the sub-control volume's face in
- * direction of the face normal.
- */
- Scalar normalMatrixHeatFlux() const
- { return normalMatrixHeatFlux_; }
private:
+ Scalar lambdaEff_;
Scalar normalMatrixHeatFlux_;
+ DimVector temperatureGrad_;
+ int faceIdx_;
};
-} // end namepace
+} // end namespace
#endif
diff --git a/dumux/implicit/2pni/2pniproperties.hh b/dumux/implicit/2pni/2pniproperties.hh
index 7a77f2d0a95787cb03514280c5e65253fc11a39f..a885e82df3c37c1621efb2f4b46a2021cc7cdeca 100644
--- a/dumux/implicit/2pni/2pniproperties.hh
+++ b/dumux/implicit/2pni/2pniproperties.hh
@@ -43,6 +43,12 @@ NEW_TYPE_TAG(TwoPNI, INHERITS_FROM(TwoP));
NEW_TYPE_TAG(BoxTwoPNI, INHERITS_FROM(BoxModel, TwoPNI));
NEW_TYPE_TAG(CCTwoPNI, INHERITS_FROM(CCModel, TwoPNI));
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(ThermalConductivityModel); //!< The model for the effective thermal conductivity
+
}
}
diff --git a/dumux/implicit/2pni/2pnipropertydefaults.hh b/dumux/implicit/2pni/2pnipropertydefaults.hh
index d19cbe21b6def819f8b2b0081a1342fc63755b46..b0566265c0543c44955e380c6a7ab7d3f6572428 100644
--- a/dumux/implicit/2pni/2pnipropertydefaults.hh
+++ b/dumux/implicit/2pni/2pnipropertydefaults.hh
@@ -36,6 +36,8 @@
#include "2pnifluxvariables.hh"
#include "2pniindices.hh"
+#include <dumux/material/fluidmatrixinteractions/2p/somerton.hh>
+
namespace Dumux
{
@@ -65,6 +67,9 @@ SET_TYPE_PROP(TwoPNI, FluxVariables, TwoPNIFluxVariables<TypeTag>);
//! The indices required by the non-isothermal two-phase model
SET_TYPE_PROP(TwoPNI, Indices, TwoPNIIndices<TypeTag, 0>);
+//! Somerton is used as default model to compute the effective thermal heat conductivity
+SET_TYPE_PROP(TwoPNI, ThermalConductivityModel, Somerton<TypeTag>);
+
}
}
diff --git a/dumux/implicit/2pni/2pnivolumevariables.hh b/dumux/implicit/2pni/2pnivolumevariables.hh
index d3804ab6bec8cb72f26f36b267b771d4c123eec9..8a5a4e89f77346a53eef78006ca73c3eda2fda68 100644
--- a/dumux/implicit/2pni/2pnivolumevariables.hh
+++ b/dumux/implicit/2pni/2pnivolumevariables.hh
@@ -79,6 +79,14 @@ public:
Scalar heatCapacity() const
{ return heatCapacity_; };
+ /*!
+ * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of the fluid phase in
+ * the sub-control volume.
+ */
+ Scalar thermalConductivity(const int phaseIdx) const
+ { return FluidSystem::thermalConductivity(this->fluidState_, phaseIdx); };
+
+
protected:
// this method gets called by the parent class. since this method
// is protected, we are friends with our parent..
diff --git a/test/implicit/2p2c/injectionspatialparams.hh b/test/implicit/2p2c/injectionspatialparams.hh
index a9726fdbcc44c07a6dc9451f0a6aae315cfb96bc..193e483fc6858ba12b7db5f16ab92974279881e3 100644
--- a/test/implicit/2p2c/injectionspatialparams.hh
+++ b/test/implicit/2p2c/injectionspatialparams.hh
@@ -116,6 +116,9 @@ public:
finePorosity_ = 0.3;
coarsePorosity_ = 0.3;
+ // heat conductivity of granite
+ lambdaSolid_ = 2.8;
+
// residual saturations
fineMaterialParams_.setSwr(0.2);
fineMaterialParams_.setSnr(0.0);
@@ -196,7 +199,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
{
@@ -206,51 +209,17 @@ public:
* (1 - porosity(element, fvGeometry, scvIdx));
}
+
/*!
- * \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 tempGrad 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 &tempGrad,
- 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 sW = 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(sW) * (ldry - lsat);
-
- // the matrix heat flux is the negative temperature gradient
- // times the heat conductivity.
- heatFlux = tempGrad;
- heatFlux *= -heatCond;
+ return lambdaSolid_;
}
private:
@@ -264,6 +233,8 @@ private:
Scalar finePorosity_;
Scalar coarsePorosity_;
+ Scalar lambdaSolid_;
+
MaterialLawParams fineMaterialParams_;
MaterialLawParams coarseMaterialParams_;
};