From 8a1d49e40274941fd64af017dcc3c5ad9351f7ca Mon Sep 17 00:00:00 2001
From: Philipp Nuske <philipp.nuske@mailbox.org>
Date: Wed, 9 Apr 2014 16:18:35 +0000
Subject: [PATCH] adding a new specialization to the mpnc models: - it is now
possible to specify the number of solved energy equations (property
NumEnergyEquations). Either 0 isothermal 1 standard-nonisothermal
2 *new* one energy equation for the void-space and one for the solid 3 one
energy equation for each phase (wetting, non-wetting, solid)
Therefore, the boolean property EnableKineticEnergy had to be renamed to the integer property NumEnergyEquations
- solve bug in test_boxmpnckinetic
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@12741 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
CHANGELOG | 3 +
.../mpnc/energy/mpncfluxvariablesenergy.hh | 9 +-
.../energy/mpncfluxvariablesenergykinetic.hh | 177 +++++++-
.../implicit/mpnc/energy/mpncindicesenergy.hh | 9 +-
.../mpnc/energy/mpncindicesenergykinetic.hh | 55 ++-
.../mpnc/energy/mpnclocalresidualenergy.hh | 23 +-
.../energy/mpnclocalresidualenergykinetic.hh | 406 +++++++++++++++++-
.../mpnc/energy/mpncvolumevariablesenergy.hh | 11 +-
.../mpncvolumevariablesenergykinetic.hh | 217 +++++++++-
.../mpnc/energy/mpncvtkwriterenergy.hh | 6 +-
.../mpnc/energy/mpncvtkwriterenergykinetic.hh | 240 ++++++++++-
.../mpnc/mass/mpnclocalresidualmass.hh | 43 +-
.../mpnc/mass/mpnclocalresidualmasskinetic.hh | 129 ++----
.../mass/mpncvolumevariablesmasskinetic.hh | 49 ++-
.../mpnc/mass/mpncvtkwritermasskinetic.hh | 23 +
dumux/implicit/mpnc/mpncfluxvariables.hh | 4 +-
dumux/implicit/mpnc/mpncindices.hh | 6 +-
dumux/implicit/mpnc/mpnclocalresidual.hh | 17 +-
dumux/implicit/mpnc/mpncmodel.hh | 16 +-
dumux/implicit/mpnc/mpncmodelkinetic.hh | 226 +++++-----
dumux/implicit/mpnc/mpncproperties.hh | 4 +-
dumux/implicit/mpnc/mpncpropertieskinetic.hh | 4 +-
dumux/implicit/mpnc/mpncpropertydefaults.hh | 6 +-
dumux/implicit/mpnc/mpncvolumevariables.hh | 11 +-
dumux/implicit/mpnc/mpncvolumevariablesia.hh | 4 +-
.../mpnc/mpncvolumevariablesiakinetic.hh | 148 ++++++-
.../thermalconductivitysimplefluidlumping.hh | 81 ++++
.../2p/vangenuchtenoftemperature.hh | 1 +
.../mpnc/evaporationatmosphereproblem.hh | 20 +-
.../evaporationatmospherespatialparams.hh | 4 +-
test/implicit/mpnc/test_boxmpnckinetic.cc | 2 +-
31 files changed, 1641 insertions(+), 313 deletions(-)
create mode 100644 dumux/material/fluidmatrixinteractions/2p/thermalconductivitysimplefluidlumping.hh
diff --git a/CHANGELOG b/CHANGELOG
index 9771d4d32b..65088aaaf8 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -6,6 +6,9 @@ Differences Between DuMuX 2.5 and DuMuX 2.6
* IMPROVEMENTS and ENHANCEMENTS:
* IMMEDIATE INTERFACE CHANGES not allowing/requiring a deprecation period:
+- in the fully implicit mpnc model a further specialization allows now to describe two-phase flow with
+two energy equations. The boolean property EnableKineticEnergy has therefore been changed to the
+integer property NumEnergyEquations
* Deprecated CLASSES/FILES, to be removed after 2.6:
diff --git a/dumux/implicit/mpnc/energy/mpncfluxvariablesenergy.hh b/dumux/implicit/mpnc/energy/mpncfluxvariablesenergy.hh
index 1a8c8dcf35..0371ce3aed 100644
--- a/dumux/implicit/mpnc/energy/mpncfluxvariablesenergy.hh
+++ b/dumux/implicit/mpnc/energy/mpncfluxvariablesenergy.hh
@@ -31,19 +31,18 @@
namespace Dumux
{
-
/*!
* \ingroup MPNCModel
* \ingroup ImplicitFluxVariables
* \brief Variables for the enthalpy fluxes in the MpNc model
*/
-template <class TypeTag, bool enableEnergy/*=false*/, bool kineticEnergyTransfer/*=false*/>
+template <class TypeTag, bool enableEnergy/*=false*/, int numEnergyEquations/*=0*/>
class MPNCFluxVariablesEnergy
{
- static_assert(!(kineticEnergyTransfer && !enableEnergy),
+ static_assert(!(numEnergyEquations && !enableEnergy),
"No kinetic energy transfer may only be enabled "
"if energy is enabled in general.");
- static_assert(!kineticEnergyTransfer,
+ static_assert(!numEnergyEquations,
"No kinetic energy transfer module included, "
"but kinetic energy transfer enabled.");
@@ -83,7 +82,7 @@ public:
};
template <class TypeTag>
-class MPNCFluxVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*kineticEnergyTransfer=*/false>
+class MPNCFluxVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*numEnergyEquations=*/1>
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
diff --git a/dumux/implicit/mpnc/energy/mpncfluxvariablesenergykinetic.hh b/dumux/implicit/mpnc/energy/mpncfluxvariablesenergykinetic.hh
index 06a4bd7e64..a8c699add9 100644
--- a/dumux/implicit/mpnc/energy/mpncfluxvariablesenergykinetic.hh
+++ b/dumux/implicit/mpnc/energy/mpncfluxvariablesenergykinetic.hh
@@ -33,8 +33,11 @@
namespace Dumux
{
+/*!
+ * \brief Specialization for the case of *3* energy balance equations.
+ */
template <class TypeTag>
-class MPNCFluxVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*kineticEnergyTransfer=*/true>
+class MPNCFluxVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*numEnergyEquations=*/3>
{
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
@@ -110,8 +113,180 @@ public:
return temperatureGradient_[energyEqIdx];
}
+private:
+ DimVector temperatureGradient_[numEnergyEqs];
+};
+
+
+/*!
+ * \brief Specialization for the case of *2* energy balance equations.
+ */
+template <class TypeTag>
+class MPNCFluxVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*numEnergyEquations=*/2>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+ typedef typename GridView::ctype CoordScalar;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ enum {dim = GridView::dimension};
+ enum {numEnergyEqs = Indices::numPrimaryEnergyVars};
+ enum {wPhaseIdx = FluidSystem::wPhaseIdx};
+ enum {nPhaseIdx = FluidSystem::nPhaseIdx};
+
+
+ typedef Dune::FieldVector<CoordScalar, dim> DimVector;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename FVElementGeometry::SubControlVolume SCV;
+ typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+ typedef typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel) ThermalConductivityModel;
+
+
+public:
+ /*!
+ * \brief The constructor
+ */
+ MPNCFluxVariablesEnergy()
+ {}
+ /*!
+ * \brief update
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry in the fully implicit scheme
+ * \param face The SCV (sub-control-volume) face
+ * \param fluxVars The flux variables
+ * \param elemVolVars The volume variables of the current element
+ */
+ void update(const Problem & problem,
+ const Element & element,
+ const FVElementGeometry & fvGeometry,
+ const SCVFace & face,
+ const FluxVariables & fluxVars,
+ const ElementVolumeVariables & elemVolVars)
+ {
+ // calculate temperature gradient using finite element
+ // gradients
+ DimVector tmp ;
+
+ for(int energyEqIdx=0; energyEqIdx<numEnergyEqs; energyEqIdx++)
+ temperatureGradient_[energyEqIdx] = 0.;
+
+ for (unsigned int idx = 0;
+ idx < face.numFap;
+ idx++){
+ // FE gradient at vertex idx
+ const DimVector & feGrad = face.grad[idx];
+
+ for (int energyEqIdx =0; energyEqIdx < numEnergyEqs; ++energyEqIdx){
+ // index for the element volume variables
+ int volVarsIdx = face.fapIndices[idx];
+
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].temperature(energyEqIdx);
+ temperatureGradient_[energyEqIdx] += tmp;
+ }
+ }
+
+ lambdaEff_ = 0;
+ calculateEffThermalConductivity_(problem,
+ element,
+ fvGeometry,
+ face,
+ elemVolVars);
+
+
+ }
+
+ /*!
+ * \brief The lumped / average conductivity of solid plus phases \f$[W/mK]\f$.
+ */
+ Scalar lambdaEff() const
+ { return lambdaEff_; }
+
+ /*!
+ * \brief The total heat flux \f$[J/s]\f$ due to heat conduction
+ * of the rock matrix over the sub-control volume's face.
+ *
+ * \param energyEqIdx The index of the energy equation
+ */
+ DimVector temperatureGradient(const unsigned int energyEqIdx) const
+ {
+ return temperatureGradient_[energyEqIdx];
+ }
+
+protected:
+ /*!
+ * \brief Calculate the effective thermal conductivity of
+ * the porous medium plus residing phases \f$[W/mK]\f$.
+ * This basically means to access the model for averaging
+ * the individual conductivities, set by the property ThermalConductivityModel.
+ * Except the adapted arguments, this is the same function
+ * as used in the implicit TwoPTwoCNIFluxVariables.
+ */
+ void calculateEffThermalConductivity_(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry & fvGeometry,
+ const SCVFace & face,
+ const ElementVolumeVariables &elemVolVars)
+ {
+ const unsigned i = face.i;
+ const unsigned j = face.j;
+ Scalar lambdaI, lambdaJ;
+
+ if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
+ {
+ lambdaI =
+ ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].fluidState().saturation(wPhaseIdx),
+ elemVolVars[i].thermalConductivity(wPhaseIdx),
+ elemVolVars[i].thermalConductivity(nPhaseIdx),
+ problem.spatialParams().thermalConductivitySolid(element, fvGeometry, i),
+ problem.spatialParams().porosity(element, fvGeometry, i));
+ lambdaJ =
+ ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].fluidState().saturation(wPhaseIdx),
+ elemVolVars[j].thermalConductivity(wPhaseIdx),
+ elemVolVars[j].thermalConductivity(nPhaseIdx),
+ problem.spatialParams().thermalConductivitySolid(element, fvGeometry, j),
+ problem.spatialParams().porosity(element, fvGeometry, j));
+ }
+ else
+ {
+ const Element & elementI = *fvGeometry.neighbors[i];
+ FVElementGeometry fvGeometryI;
+ fvGeometryI.subContVol[0].global = elementI.geometry().center();
+
+ lambdaI =
+ ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].fluidState().saturation(wPhaseIdx),
+ elemVolVars[i].thermalConductivity(wPhaseIdx),
+ elemVolVars[i].thermalConductivity(nPhaseIdx),
+ problem.spatialParams().thermalConductivitySolid(elementI, fvGeometryI, 0),
+ problem.spatialParams().porosity(elementI, fvGeometryI, 0));
+
+ const Element & elementJ = *fvGeometry.neighbors[j];
+ FVElementGeometry fvGeometryJ;
+ fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
+
+ lambdaJ =
+ ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].fluidState().saturation(wPhaseIdx),
+ elemVolVars[j].thermalConductivity(wPhaseIdx),
+ elemVolVars[j].thermalConductivity(nPhaseIdx),
+ problem.spatialParams().thermalConductivitySolid(elementJ, fvGeometryJ, 0),
+ problem.spatialParams().porosity(elementJ, fvGeometryJ, 0));
+ }
+
+ // -> arithmetic mean, open to discussion
+ lambdaEff_ = 0.5 * (lambdaI+lambdaJ);
+ }
+
private:
+ Scalar lambdaEff_ ;
DimVector temperatureGradient_[numEnergyEqs];
};
diff --git a/dumux/implicit/mpnc/energy/mpncindicesenergy.hh b/dumux/implicit/mpnc/energy/mpncindicesenergy.hh
index ec3cfe751e..fbe135b8b0 100644
--- a/dumux/implicit/mpnc/energy/mpncindicesenergy.hh
+++ b/dumux/implicit/mpnc/energy/mpncindicesenergy.hh
@@ -33,13 +33,14 @@ namespace Dumux
*
* This is a dummy class for the isothermal case.
*/
-template <int PVOffset, bool enableEnergy/*=false*/, bool kineticEnergyTransfer/*=false*/>
+template <int PVOffset, bool enableEnergy/*=false*/, int numEnergyEquations/*=0*/>
struct MPNCEnergyIndices
{
- static_assert(!(kineticEnergyTransfer && !enableEnergy),
+ static_assert(((numEnergyEquations<1) and not enableEnergy),
"No kinetic energy transfer may only be enabled "
"if energy is enabled in general.");
- static_assert(!kineticEnergyTransfer,
+
+ static_assert( (numEnergyEquations < 1) ,
"No kinetic energy transfer module included, "
"but kinetic energy transfer enabled.");
public:
@@ -56,7 +57,7 @@ public:
* \brief The indices required for the energy equation.
*/
template <int PVOffset>
-struct MPNCEnergyIndices<PVOffset, /*isNonIsothermal=*/true, /*kineticEnergyTransfer*/false>
+struct MPNCEnergyIndices<PVOffset, /*isNonIsothermal=*/true, /*numEnergyEquations=*/ 1 >
{
public:
/*!
diff --git a/dumux/implicit/mpnc/energy/mpncindicesenergykinetic.hh b/dumux/implicit/mpnc/energy/mpncindicesenergykinetic.hh
index 7a4cef1a1a..f69892fdd2 100644
--- a/dumux/implicit/mpnc/energy/mpncindicesenergykinetic.hh
+++ b/dumux/implicit/mpnc/energy/mpncindicesenergykinetic.hh
@@ -29,10 +29,11 @@ namespace Dumux
{
/*!
- * \brief The indices required for the energy equation.
+ * \brief The indices required for the energy equation. Specialization for the case of
+ * *3* energy balance equations.
*/
template <int PVOffset>
-struct MPNCEnergyIndices<PVOffset, /*enableEnergy=*/true, /*kineticEnergyTransfer=*/true>
+struct MPNCEnergyIndices<PVOffset, /*enableEnergy=*/true, /*numEnergyEquations=*/3>
{
public:
/*!
@@ -60,12 +61,58 @@ public:
static const unsigned int energyEqIdx = energyEq0Idx;
};
+/*!
+ * \brief The indices required for the energy equation. Specialization for the case of
+ * *2* energy balance equations.
+ */
+template <int PVOffset>
+struct MPNCEnergyIndices<PVOffset, /*enableEnergy=*/true, /*numEnergyEquations=*/2>
+{
+public:
+ /*!
+ * \brief This module defines one new primary variable.
+ */
+ static const unsigned int numPrimaryVars = 2;
+
+ /*!
+ * \brief Index for the temperature of the wetting phase in a vector of primary
+ * variables.
+ */
+ static const unsigned int temperature0Idx = PVOffset + 0;
+
+ /*!
+ * \brief Compatibility with non kinetic models
+ */
+ static const unsigned int temperatureIdx = temperature0Idx;
+ /*!
+ * \brief Equation index of the energy equation.
+ */
+ static const unsigned int energyEq0Idx = PVOffset + 0;
+ /*!
+ * \brief Equation index of the energy equation.
+ */
+ static const unsigned int energyEqSolidIdx = energyEq0Idx + numPrimaryVars - 1 ;
+
+ /*!
+ * \brief Index for storing e.g. temperature fluidState
+ */
+ static const unsigned int temperatureFluidIdx = 0 ;
+ static const unsigned int temperatureSolidIdx = 1 ;
+
+
+ /*!
+ * \brief Compatibility with non kinetic models
+ */
+ static const unsigned int energyEqIdx = energyEq0Idx;
+
+
+};
/*!
* \brief The indices required for the energy equation.
*/
template <int PVOffset, bool isNonIsothermal>
-struct MPNCEnergyIndices<PVOffset, isNonIsothermal, /*kineticEnergyTransfer=*/false >
+struct MPNCEnergyIndices<PVOffset, isNonIsothermal, /*numEnergyEquations=*/1 >
{
public:
/*!
@@ -90,7 +137,7 @@ public:
* This is a dummy class for the isothermal case.
*/
template <int PVOffset>
-struct MPNCEnergyIndices<PVOffset, /*isNonIsothermal=*/false, /*kineticEnergyTransfer*/false>
+struct MPNCEnergyIndices<PVOffset, /*isNonIsothermal=*/false, /*numEnergyEquations*/0>
{
public:
/*!
diff --git a/dumux/implicit/mpnc/energy/mpnclocalresidualenergy.hh b/dumux/implicit/mpnc/energy/mpnclocalresidualenergy.hh
index e1bec1099b..2ad7657a34 100644
--- a/dumux/implicit/mpnc/energy/mpnclocalresidualenergy.hh
+++ b/dumux/implicit/mpnc/energy/mpnclocalresidualenergy.hh
@@ -34,13 +34,13 @@ namespace Dumux {
*
* This class just does nothing.
*/
-template <class TypeTag, bool enableEnergy/*=false*/, bool kineticEnergyTransfer /*=false*/>
+template <class TypeTag, bool enableEnergy/*=false*/, int numEnergyEquations /*=0*/>
class MPNCLocalResidualEnergy
{
- static_assert(!(kineticEnergyTransfer && !enableEnergy),
+ static_assert(!(numEnergyEquations && !enableEnergy),
"No kinetic energy transfer may only be enabled "
"if energy is enabled in general.");
- static_assert(!kineticEnergyTransfer,
+ static_assert(!numEnergyEquations,
"No kinetic energy transfer module included, "
"but kinetic energy transfer enabled.");
@@ -121,7 +121,7 @@ public:
template <class TypeTag>
-class MPNCLocalResidualEnergy<TypeTag, /*enableEnergy=*/true, /*kineticenergyTransfer=*/false>
+class MPNCLocalResidualEnergy<TypeTag, /*enableEnergy=*/true, /*numEnergyEquations=*/ 1 >
{
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
@@ -186,6 +186,11 @@ public:
* fs.internalEnergy(phaseIdx)
* fs.saturation(phaseIdx)
* volVars.porosity();
+
+#ifndef NDEBUG
+if (!std::isfinite(storage[energyEqIdx]))
+ DUNE_THROW(NumericalProblem, "Calculated non-finite energy storage");
+#endif
}
/*!
* \brief Evaluates the total flux of all conservation quantities
@@ -244,6 +249,10 @@ public:
// the enthalpy transport
const VolumeVariables &up = elemVolVars[upIdx];
flux[energyEqIdx] += up.fluidState().enthalpy(phaseIdx) * massFlux;
+#ifndef NDEBUG
+if (!std::isfinite(flux[energyEqIdx]) )
+ DUNE_THROW(NumericalProblem, "Calculated non-finite energy flux");
+#endif
}
/*!
* \brief The heat conduction in the phase
@@ -260,7 +269,11 @@ public:
Scalar lumpedConductivity = fluxVars.fluxVarsEnergy().lambdaEff() ;
Scalar temperatureGradientNormal = fluxVars.fluxVarsEnergy().temperatureGradientNormal() ;
Scalar lumpedHeatConduction = - lumpedConductivity * temperatureGradientNormal ;
- flux[energyEqIdx] += lumpedHeatConduction ;
+ flux[energyEqIdx] += lumpedHeatConduction;
+#ifndef NDEBUG
+if (!std::isfinite(flux[energyEqIdx]) )
+ DUNE_THROW(NumericalProblem, "Calculated non-finite energy flux");
+#endif
}
/*!
diff --git a/dumux/implicit/mpnc/energy/mpnclocalresidualenergykinetic.hh b/dumux/implicit/mpnc/energy/mpnclocalresidualenergykinetic.hh
index c235fb6f99..f2f3c35b7b 100644
--- a/dumux/implicit/mpnc/energy/mpnclocalresidualenergykinetic.hh
+++ b/dumux/implicit/mpnc/energy/mpnclocalresidualenergykinetic.hh
@@ -26,12 +26,17 @@
#define DUMUX_MPNC_LOCAL_RESIDUAL_ENERGY_KINETIC_HH
#include <dumux/implicit/mpnc/mpnclocalresidual.hh>
+#include <dumux/common/spline.hh>
+
namespace Dumux
{
+/*!
+ * \brief Specialization for the case of *3* energy balance equations.
+ */
template <class TypeTag>
-class MPNCLocalResidualEnergy<TypeTag, /*enableEnergy=*/true, /*kineticEnergyTransfer=*/true>
+class MPNCLocalResidualEnergy<TypeTag, /*enableEnergy=*/true, /*numEnergyEquations=*/3>
{
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
@@ -123,10 +128,9 @@ public:
else
DUNE_THROW(Dune::NotImplemented,
"wrong index");
- #ifndef NDEBUG
+
if (!std::isfinite(storage[energyEq0Idx+phaseIdx]))
DUNE_THROW(NumericalProblem, "Calculated non-finite storage");
- #endif
}
/*!
@@ -161,10 +165,8 @@ public:
fluxVars,
elemVolVars,
energyEqIdx);
- #ifndef NDEBUG
if (!std::isfinite(flux[energyEq0Idx + energyEqIdx]))
- DUNE_THROW(NumericalProblem, "Calculated non-finite flux");
- #endif
+ DUNE_THROW(NumericalProblem, "Calculated non-finite flux in phase " << energyEqIdx);
}
}
/*!
@@ -185,26 +187,37 @@ public:
// calculate the mass flux in the phase i.e. make mass flux out
// of mole flux and add up the fluxes of a phase
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx){
massFlux += molarComponentValuesMassTransport[compIdx]
* FluidSystem::molarMass(compIdx) ;
+ }
unsigned int upIdx = fluxVars.face().i;
+ unsigned int dnIdx = fluxVars.face().j;
if (massFlux < 0){
upIdx = fluxVars.face().j;
+ dnIdx = fluxVars.face().i;
}
// use the phase enthalpy of the upstream vertex to calculate
// the enthalpy transport
const VolumeVariables & up = elemVolVars[upIdx];
+ const VolumeVariables & dn = elemVolVars[dnIdx];
/* todo
* CAUTION: this is not exactly correct: does diffusion carry the upstream phase enthalpy?
* To be more precise this should be the components enthalpy.
* In the same vein: Counter current diffusion is not accounted for here.
*/
- const Scalar enthalpy = up.fluidState().enthalpy(phaseIdx) ;
- flux[energyEq0Idx + phaseIdx] += enthalpy * massFlux ;
+ const Scalar transportedThingUp = up.fluidState().enthalpy(phaseIdx) ;
+ const Scalar transportedThingDn = dn.fluidState().enthalpy(phaseIdx) ;
+
+
+ const Scalar massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
+ flux[energyEq0Idx + phaseIdx] += massFlux *
+ (massUpwindWeight_ * transportedThingUp
+ +
+ (1.-massUpwindWeight_) * transportedThingDn ) ;
}
/*!
* \brief The heat conduction in the phase
@@ -236,7 +249,7 @@ public:
const Scalar klambda = kVolVar.thermalConductivity(phaseIdx);
// todo which average?
// Using a harmonic average is justified by its properties: if one phase does not conduct energy, there is no transfer
- const Scalar barLambda = Dumux::harmonicMean(ilambda, klambda);
+ const Scalar barLambda = Dumux::harmonicMean(ilambda, klambda) ;
const Scalar gradientNormal = fluxVars.fluxVarsEnergy().temperatureGradient(phaseIdx)
* fluxVars.face().normal ;
@@ -336,10 +349,8 @@ public:
- #ifndef NDEBUG
if (!std::isfinite(source[energyEq0Idx + phaseIdx]))
DUNE_THROW(NumericalProblem, "Calculated non-finite source, " << "Tw="<< Tw << " Tn="<< Tn<< " Ts="<< Ts);
- #endif
}// end phases
#define MASS_ENERGY_TRANSPORT 1
@@ -389,6 +400,377 @@ public:
}// end source
};
+
+/*!
+ * \brief Specialization for the case of *2* energy balance equations.
+ */
+template <class TypeTag>
+class MPNCLocalResidualEnergy<TypeTag, /*enableEnergy=*/true, /*numEnergyEquations=*/2>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+ enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+ enum { energyEq0Idx = Indices::energyEq0Idx };
+ enum { numEnergyEqs = Indices::numPrimaryEnergyVars};
+ enum { wPhaseIdx = FluidSystem::wPhaseIdx};
+ enum { nPhaseIdx = FluidSystem::nPhaseIdx};
+ enum { nCompIdx = FluidSystem::nCompIdx};
+ enum { wCompIdx = FluidSystem::wCompIdx};
+ enum { sPhaseIdx = FluidSystem::sPhaseIdx};
+ enum { energyEqSolidIdx = Indices::energyEqSolidIdx};
+ enum { temperatureFluidIdx = Indices::temperatureFluidIdx};
+ enum { temperatureSolidIdx = Indices::temperatureSolidIdx};
+ enum { dim = GridView::dimension}; // Grid and world dimension
+
+
+ typedef typename Dune::FieldVector<Scalar, numComponents> ComponentVector;
+ typedef typename Dune::FieldMatrix<Scalar, numPhases, numComponents> PhaseComponentMatrix;
+ typedef typename FluidSystem::ParameterCache ParameterCache;
+
+ typedef typename GET_PROP_TYPE(TypeTag, BaseLocalResidual) ParentType;
+ typedef Dumux::Spline<Scalar> Spline;
+
+
+public:
+ /*! \brief Evaluate the amount all conservation quantities
+ * (e.g. phase mass) within a sub-control volume.
+ *
+ * The result should be averaged over the volume (e.g. phase mass
+ * inside a sub-control volume divided by the volume)
+ *
+ * \param storage The mass of the component within the sub-control volume
+ * \param volVars the Volume Variables
+ */
+ static void computeStorage(PrimaryVariables & storage,
+ const VolumeVariables & volVars)
+ {
+ for(int energyEqIdx=0; energyEqIdx< numEnergyEqs; energyEqIdx++)
+ storage[energyEq0Idx+energyEqIdx] = 0;
+
+ // energy of the fluids
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ addPhaseStorage(storage, volVars, phaseIdx);
+ }
+ // heat stored in the rock matrix
+ storage[energyEqSolidIdx] +=
+ volVars.temperature(temperatureSolidIdx)
+ * volVars.densitySolid()
+ * (1.0 - volVars.porosity())
+ * volVars.heatCapacity();
+
+ if (!std::isfinite(storage[energyEqSolidIdx]))
+ DUNE_THROW(NumericalProblem, "Calculated non-finite storage");
+ }
+
+ /*! \brief Calculate the storage for all mass balance equations
+ * within a single fluid phase
+ *
+ * \param storage The mass of the component within the sub-control volume
+ * \param volVars the Volume Variables
+ * \param phaseIdx The local index of the phases
+ */
+ static void addPhaseStorage(PrimaryVariables & storage,
+ const VolumeVariables & volVars,
+ const unsigned int phaseIdx)
+ {
+
+ const FluidState & fs = volVars.fluidState();
+
+ if (phaseIdx not_eq sPhaseIdx) {
+ storage[energyEq0Idx ] +=
+ fs.density(phaseIdx) *
+ fs.internalEnergy(phaseIdx) *
+ fs.saturation(phaseIdx) *
+ volVars.porosity();
+ }
+ else
+ DUNE_THROW(Dune::NotImplemented,
+ "wrong index");
+
+ if (!std::isfinite(storage[energyEq0Idx]))
+ DUNE_THROW(NumericalProblem, "Calculated non-finite storage");
+ }
+
+ /*!\brief Evaluates the total flux of all conservation quantities
+ * over a face of a sub-control volume.
+ *
+ * \param flux The flux over the SCV (sub-control-volume) face for each component
+ * \param fluxVars The flux Variables
+ * \param elemVolVars The volume variables of the current element
+ * \param molarPhaseComponentValuesMassTransport
+ */
+ static void computeFlux(PrimaryVariables & flux,
+ const FluxVariables & fluxVars,
+ const ElementVolumeVariables & elemVolVars,
+ const ComponentVector molarPhaseComponentValuesMassTransport[numPhases])
+ {
+ // reset all energy fluxes
+ for(int energyEqIdx=0; energyEqIdx<numEnergyEqs; ++energyEqIdx)
+ flux[energyEq0Idx + energyEqIdx] = 0.0;
+
+ // only the fluid phases transport enthalpy
+ for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx)
+ computePhaseEnthalpyFlux(flux,
+ fluxVars,
+ elemVolVars,
+ phaseIdx,
+ molarPhaseComponentValuesMassTransport[phaseIdx]);
+
+ // all phases take part in conduction
+ for(int energyEqIdx=0; energyEqIdx<numEnergyEqs; ++energyEqIdx){
+ computeHeatConduction(flux,
+ fluxVars,
+ elemVolVars,
+ energyEqIdx);
+
+ if (!std::isfinite(flux[energyEq0Idx + energyEqIdx]))
+ DUNE_THROW(NumericalProblem, "Calculated non-finite flux in phase " << energyEqIdx);
+ }
+ }
+
+ /*! \brief the advective Flux of the enthalpy
+ * \param flux The flux over the SCV (sub-control-volume) face for each component
+ * \param fluxVars The flux Variables
+ * \param elemVolVars The volume variables of the current element
+ * \param phaseIdx The local index of the phases
+ * \param molarComponentValuesMassTransport
+ */
+ static void computePhaseEnthalpyFlux(PrimaryVariables & flux,
+ const FluxVariables & fluxVars,
+ const ElementVolumeVariables & elemVolVars,
+ const unsigned int phaseIdx,
+ const ComponentVector & molarComponentValuesMassTransport)
+ {
+ Scalar massFlux = 0; // mass flux is not the perfect term: sum_kappa (rho_alpha v_alpha x_alpha^kappa) = v_alpha rho_alpha
+
+ // calculate the mass flux in the phase i.e. make mass flux out
+ // of mole flux and add up the fluxes of a phase
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx){
+ massFlux += molarComponentValuesMassTransport[compIdx]
+ * FluidSystem::molarMass(compIdx) ;
+ }
+
+ unsigned int upIdx = fluxVars.face().i;
+ unsigned int dnIdx = fluxVars.face().j;
+ if (massFlux < 0){
+ upIdx = fluxVars.face().j;
+ dnIdx = fluxVars.face().i;
+ }
+
+ // use the phase enthalpy of the upstream vertex to calculate
+ // the enthalpy transport
+ const VolumeVariables & up = elemVolVars[upIdx];
+ const VolumeVariables & dn = elemVolVars[dnIdx];
+
+ /* todo
+ * CAUTION: this is not exactly correct: does diffusion carry the upstream phase enthalpy? To be more precise this should be the components enthalpy. In the same vein: Counter current diffusion is not accounted for here.
+ */
+ const Scalar transportedThingUp = up.fluidState().enthalpy(phaseIdx) ;
+ const Scalar transportedThingDn = dn.fluidState().enthalpy(phaseIdx) ;
+
+ const Scalar massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
+ flux[energyEq0Idx] += massFlux *
+ (massUpwindWeight_ * transportedThingUp
+ +
+ (1.-massUpwindWeight_) * transportedThingDn ) ;
+ }
+ /*! \brief The heat conduction in the phase
+ *
+ * \param flux The flux over the SCV (sub-control-volume) face for each component
+ * \param fluxVars The flux Variables
+ * \param elemVolVars The volume variables of the current element
+ * \param phaseIdx The local index of the phases
+ */
+ static void computeHeatConduction(PrimaryVariables & flux,
+ const FluxVariables & fluxVars,
+ const ElementVolumeVariables & elemVolVars,
+ const unsigned int energyEqIdx)
+ {
+ const unsigned int iIdx = fluxVars.face().i;
+ const unsigned int kIdx = fluxVars.face().j; // k can be better distinguished from i
+
+ const VolumeVariables & iVolVar = elemVolVars[iIdx];
+ const VolumeVariables & kVolVar = elemVolVars[kIdx];
+
+ const Scalar iPorosity = iVolVar.porosity();
+ const Scalar kPorosity = kVolVar.porosity();
+ const Scalar barPorosity = Dumux::harmonicMean(iPorosity, kPorosity);
+
+ const Scalar iSolidLambda = iVolVar.thermalConductivity(sPhaseIdx);
+ const Scalar kSolidLambda = kVolVar.thermalConductivity(sPhaseIdx);
+
+ // Using a harmonic average is justified by its properties: if one phase does not conduct energy, there is no transfer
+ const Scalar barSolidLambda = (iSolidLambda+kSolidLambda) / 2.0;
+
+ const Scalar lumpedConductivity = fluxVars.fluxVarsEnergy().lambdaEff() ;
+
+ const Scalar gradientNormal = fluxVars.fluxVarsEnergy().temperatureGradient(energyEqIdx)
+ * fluxVars.face().normal ;
+
+ // heat conduction of the rock matrix and the fluid phases
+ if (energyEqIdx == temperatureSolidIdx){
+ flux[energyEqSolidIdx] -= barSolidLambda * (1.-barPorosity) * gradientNormal ;
+ }
+ else if (energyEqIdx == temperatureFluidIdx){
+ flux[energyEq0Idx ] -= lumpedConductivity /*already includes porosity*/ * gradientNormal ;
+ }
+ else
+ DUNE_THROW(Dune::NotImplemented,
+ "wrong index");
+ }
+
+ /*! \brief Calculate the source term of the equation
+ *
+ * \param source The source/sink in the sub-control volume for each component
+ * \param volVars The volume variables
+ * \param componentIntoPhaseMassTransfer
+ */
+ static void computeSource(PrimaryVariables & source,
+ const VolumeVariables & volVars,
+ const ComponentVector componentIntoPhaseMassTransfer[numPhases])
+ {
+ const FluidState & fs = volVars.fluidState() ;
+
+ const Scalar characteristicLength = volVars.characteristicLength() ;
+
+ // Shi & Wang, Transport in porous media (2011)
+ const Scalar as = 6.0 * (1.0-volVars.porosity()) / characteristicLength ;
+
+ const Scalar TFluid = volVars.temperature(temperatureFluidIdx);
+ const Scalar TSolid = volVars.temperature(temperatureSolidIdx);
+
+ const Scalar satW = fs.saturation(wPhaseIdx) ;
+ const Scalar satN = fs.saturation(nPhaseIdx) ;
+
+ const Scalar eps = 1e-6 ;
+ Scalar solidToFluidEnergyExchange ;
+
+ Scalar fluidConductivity ;
+ if (satW < 1.0 - eps)
+ fluidConductivity = volVars.thermalConductivity(nPhaseIdx) ;
+ else if (satW >= 1.0 - eps)
+ fluidConductivity = volVars.thermalConductivity(wPhaseIdx) ;
+ else
+ DUNE_THROW(Dune::NotImplemented,
+ "wrong range");
+
+ const Scalar lambdaSolid = volVars.thermalConductivity(sPhaseIdx);
+ const Scalar factorEnergyTransfer = volVars.factorEnergyTransfer() ;
+
+ solidToFluidEnergyExchange = factorEnergyTransfer * (TSolid - TFluid) / characteristicLength * as * fluidConductivity ;
+
+ const Scalar epsRegul = 1e-3 ;
+
+ if (satW < (0 + eps) )
+ solidToFluidEnergyExchange *= volVars.nusseltNumber(nPhaseIdx) ;
+
+ else if ( (satW >= 0 + eps) and (satW < 1.0-eps) ){
+ solidToFluidEnergyExchange *= (volVars.nusseltNumber(nPhaseIdx) * satN );
+ Scalar qBoil ;
+
+ if (satW<=epsRegul){// regularize
+ Spline sp(0.0, epsRegul, // x1, x2
+ QBoilFunc(volVars, 0.0), QBoilFunc(volVars, epsRegul), // y1, y2
+ 0.0, dQBoil_dSw(volVars, epsRegul) ); // m1, m2
+
+ qBoil = sp.eval(satW) ;
+ }
+
+ else if (satW>= (1.0-epsRegul) ){// regularize
+ Spline sp(1.0-epsRegul, 1.0, // x1, x2
+ QBoilFunc(volVars, 1.0-epsRegul), 0.0, // y1, y2
+ dQBoil_dSw(volVars, 1.0-epsRegul), 0.0 ); // m1, m2
+
+ qBoil = sp.eval(satW) ;
+ }
+ else
+ qBoil = QBoilFunc(volVars, satW) ;
+
+ solidToFluidEnergyExchange += qBoil;
+ }
+ else if (satW >= 1.0-eps)
+ solidToFluidEnergyExchange *= volVars.nusseltNumber(wPhaseIdx) ;
+ else
+ DUNE_THROW(Dune::NotImplemented,
+ "wrong range");
+
+ for(int energyEqIdx =0; energyEqIdx<numEnergyEqs; ++energyEqIdx){
+ switch (energyEqIdx){
+ case 0 :
+ source[energyEq0Idx + energyEqIdx] = solidToFluidEnergyExchange;
+ break;
+ case 1 :
+ source[energyEq0Idx + energyEqIdx] = - solidToFluidEnergyExchange;
+ break;
+ default:
+ DUNE_THROW(Dune::NotImplemented,
+ "wrong index");
+ } // end switch
+
+ if (!std::isfinite(source[energyEq0Idx + energyEqIdx]))
+ DUNE_THROW(NumericalProblem, "Calculated non-finite source, " << "TFluid="<< TFluid << " TSolid="<< TSolid );
+
+ }// end energyEqIdx
+ Valgrind::CheckDefined(source);
+ }// end source
+
+
+ /*! \brief Calculate the energy transfer during boiling, i.e. latent heat
+ *
+ * \param volVars The volume variables
+ * \param satW The wetting phase saturation. Not taken from volVars, because we regularize.
+ */
+ static Scalar QBoilFunc(const VolumeVariables & volVars,
+ const Scalar satW)
+ {
+ // using saturation as input (instead of from volVars)
+ // in order to make regularization (evaluation at different points) easyer
+ const FluidState & fs = volVars.fluidState() ;
+ const Scalar g( 9.81 ) ;
+ const Scalar gamma(0.0589) ;
+ const Scalar TFluid = volVars.temperature(temperatureFluidIdx);
+ const Scalar TSolid = volVars.temperature(temperatureSolidIdx);
+ const Scalar characteristicLength = volVars.characteristicLength() ;
+
+ const Scalar as = 6.0 * (1.0-volVars.porosity()) / characteristicLength ;
+ const Scalar mul = fs.viscosity(wPhaseIdx) ;
+ const Scalar deltahv = fs.enthalpy(nPhaseIdx) - fs.enthalpy(wPhaseIdx);
+ const Scalar deltaRho = fs.density(wPhaseIdx) - fs.density(nPhaseIdx) ;
+ const Scalar firstBracket = std::pow(g * deltaRho / gamma, 0.5);
+ const Scalar cp = FluidSystem::heatCapacity(fs, wPhaseIdx) ;
+ const Scalar Tsat = FluidSystem::vaporTemperature(fs, nPhaseIdx ) ;
+ const Scalar deltaT = TSolid - Tsat ;
+ const Scalar secondBracket = std::pow( (cp *deltaT / (0.006 * deltahv) ) , 3.0 ) ;
+ const Scalar Prl = volVars.prandtlNumber(wPhaseIdx) ;
+ const Scalar thirdBracket = std::pow( 1/Prl , (1.7/0.33) );
+ const Scalar QBoil = satW * as * mul * deltahv * firstBracket * secondBracket * thirdBracket ;
+ return QBoil;
+ }
+
+ /*! \brief Calculate the derivative of the energy transfer function during boiling. Needed for regularization.
+ *
+ * \param volVars The volume variables
+ * \param satW The wetting phase saturation. Not taken from volVars, because we regularize.
+ */
+ static Scalar dQBoil_dSw(const VolumeVariables & volVars,
+ const Scalar satW)
+ {
+ // on the fly derive w.r.t. Sw.
+ // Only linearly depending on it (directly)
+ return (QBoilFunc(volVars, satW) / satW ) ;
+ }
+};
+
+
} // end namespace Dumux
#endif // DUMUX_MPNC_LOCAL_RESIDUAL_ENERGY_KINETIC_HH
diff --git a/dumux/implicit/mpnc/energy/mpncvolumevariablesenergy.hh b/dumux/implicit/mpnc/energy/mpncvolumevariablesenergy.hh
index 450eafcebd..11a123a1f6 100644
--- a/dumux/implicit/mpnc/energy/mpncvolumevariablesenergy.hh
+++ b/dumux/implicit/mpnc/energy/mpncvolumevariablesenergy.hh
@@ -37,13 +37,13 @@ namespace Dumux
* only isothermal in the sense that the temperature at a location and
* a time is specified outside of the model!
*/
-template <class TypeTag, bool enableEnergy/*=false*/, bool kineticEnergyTransfer /*=don't care*/>
+template <class TypeTag, bool enableEnergy/*=false*/, int numEnergyEquations /*=don't care*/>
class MPNCVolumeVariablesEnergy
{
- static_assert(!(kineticEnergyTransfer && !enableEnergy),
+ static_assert(not (numEnergyEquations and not enableEnergy),
"No kinetic energy transfer may only be enabled "
"if energy is enabled in general.");
- static_assert(!kineticEnergyTransfer,
+ static_assert(numEnergyEquations < 1,
"No kinetic energy transfer module included, "
"but kinetic energy transfer enabled.");
@@ -123,7 +123,7 @@ public:
* finite volume in the two-phase, N-component model.
*/
template <class TypeTag>
-class MPNCVolumeVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*kineticEnergyTransfer=*/false>
+class MPNCVolumeVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*numEnergyEquations=*/ 1 >
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
@@ -227,9 +227,6 @@ public:
Scalar densitySolid() const
{ return densitySolid_; }
- DUNE_DEPRECATED_MSG("use densitySolid() instead")
- Scalar soilDensity() const
- { return densitySolid(); }
/*!
* \brief If running under valgrind this produces an error message
diff --git a/dumux/implicit/mpnc/energy/mpncvolumevariablesenergykinetic.hh b/dumux/implicit/mpnc/energy/mpncvolumevariablesenergykinetic.hh
index 69540b4369..2e6f83d1e1 100644
--- a/dumux/implicit/mpnc/energy/mpncvolumevariablesenergykinetic.hh
+++ b/dumux/implicit/mpnc/energy/mpncvolumevariablesenergykinetic.hh
@@ -32,9 +32,10 @@ namespace Dumux
/*!
* \brief Contains the volume variables of the kinetic energy transfer
* module of the M-phase, N-component model.
+ * Specialization for the case of *3* energy balance equations.
*/
template <class TypeTag>
-class MPNCVolumeVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*kineticEnergyTransfer=*/true>
+class MPNCVolumeVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*numEnergyEquations=*/3>
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
@@ -110,6 +111,7 @@ public:
temperature_[phaseIdx]= T;
fluidState.setTemperature(phaseIdx, T);
}
+
temperature_[sPhaseIdx] = priVars[temperature0Idx + sPhaseIdx];
Valgrind::CheckDefined(temperature_);
@@ -182,10 +184,6 @@ public:
Scalar densitySolid() const
{ return densitySolid_; }
- DUNE_DEPRECATED_MSG("use densitySolid() instead")
- Scalar soilDensity() const
- { return densitySolid(); }
-
/*!
* \brief Returns the conductivity of the given solid phase [kg / m^3] in
* the sub-control volume.
@@ -193,10 +191,6 @@ public:
Scalar thermalConductivitySolid() const
{ return thermalConductivitySolid_; }
- DUNE_DEPRECATED_MSG("use thermalConductivitySolid() instead")
- Scalar soilThermalConductivity() const
- { return thermalConductivitySolid(); }
-
/*!
* \brief Returns the conductivity of the given fluid [kg / m^3] in
* the sub-control volume.
@@ -238,6 +232,211 @@ protected:
Scalar thermalConductivityFluid_[numPhases];
};
+
+
+/*!
+ * \brief Contains the volume variables of the kinetic energy transfer
+ * module of the M-phase, N-component model.
+ * Specialization for the case of *2* energy balance equations.
+ */
+template <class TypeTag>
+class MPNCVolumeVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*numEnergyEquations=*/2>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+ enum { temperature0Idx = Indices::temperature0Idx };
+ enum { nPhaseIdx = FluidSystem::nPhaseIdx };
+ enum { wPhaseIdx = FluidSystem::wPhaseIdx };
+ enum { sPhaseIdx = FluidSystem::sPhaseIdx };
+ enum { numEnergyEqs = Indices::numPrimaryEnergyVars};
+
+ /*!
+ * \brief The fluid state which is used by the volume variables to
+ * store the thermodynamic state.
+ *
+ * If chemical equilibrium is not considered, we use the most
+ * generic fluid state.
+ */
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename FluidSystem::ParameterCache ParameterCache;
+
+public:
+ /*!
+ * \brief Update the temperature of the sub-control volume.
+ *
+ * \param priVars The primary variables
+ * \param element The finite Element
+ * \param fvGeometry The finite-volume geometry in the fully implicit scheme
+ * \param scvIdx The index of the sub-control volume
+ * \param problem The problem
+ * \param temperatureIdx The temperature Index
+ */
+ Scalar getTemperature(const PrimaryVariables & priVars,
+ const Element & element,
+ const FVElementGeometry & fvGeometry,
+ const unsigned int scvIdx,
+ const Problem & problem,
+ const unsigned int temperatureIdx) const
+ {
+ // retrieve temperature from solution vector
+ return priVars[temperature0Idx + temperatureIdx]; // could also be solid phase
+ }
+
+ /*!
+ * \brief Update the temperature of the sub-control volume.
+ *
+ * \param fluidState Container for all the secondary variables concerning the fluids
+ * \param paramCache Container for cache parameters
+ * \param priVars The primary Variables
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry in the fully implicit scheme
+ * \param scvIdx The index of the sub-control volume
+ * \param problem The problem
+ */
+ void updateTemperatures(FluidState & fluidState,
+ ParameterCache & paramCache,
+ const PrimaryVariables & priVars,
+ const Element & element,
+ const FVElementGeometry & fvGeometry,
+ const unsigned int scvIdx,
+ const Problem & problem)
+ {
+ assert(2 == numEnergyEqs);
+
+ for(int energyEqIdx=0; energyEqIdx < numPhases; ++energyEqIdx){
+ // retrieve temperature from solution vector
+ const Scalar T = priVars[temperature0Idx + energyEqIdx];
+ temperature_[energyEqIdx]= T;
+ }
+
+ fluidState.setTemperature(priVars[temperature0Idx]);
+
+ Valgrind::CheckDefined(temperature_);
+ }
+
+ /*!
+ * \brief Update the enthalpy and the internal energy for a given
+ * control volume.
+ *
+ * \param fluidState Container for all the secondary variables concerning the fluids
+ * \param paramCache Container for cache parameters
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry in the fully implicit scheme
+ * \param scvIdx The index of the sub-control volume
+ * \param problem The problem
+ */
+ void update(FluidState & fluidState,
+ ParameterCache & paramCache,
+ const Element & element,
+ const FVElementGeometry & fvGeometry,
+ const unsigned int scvIdx,
+ const Problem & problem)
+ {
+ heatCapacity_ =
+ problem.spatialParams().heatCapacity(element, fvGeometry, scvIdx);
+ Valgrind::CheckDefined(heatCapacity_);
+
+ for(int phaseIdx =0; phaseIdx<numPhases; ++phaseIdx){
+ thermalConductivityFluid_[phaseIdx] =
+ FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
+ }
+ Valgrind::CheckDefined(thermalConductivityFluid_);
+
+
+ densitySolid_ =
+ problem.spatialParams().densitySolid(element, fvGeometry, scvIdx);
+ Valgrind::CheckDefined(densitySolid_);
+
+ thermalConductivitySolid_ =
+ problem.spatialParams().thermalConductivitySolid(element, fvGeometry, scvIdx);
+ Valgrind::CheckDefined(thermalConductivitySolid_);
+
+ // set the enthalpies
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ Scalar h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
+ Valgrind::CheckDefined(h);
+ fluidState.setEnthalpy(phaseIdx, h);
+ }
+ }
+
+ /*!
+ * \brief Returns the total heat capacity [J/(K m^3)] of the rock matrix in
+ * the sub-control volume.
+ */
+ Scalar heatCapacity() const
+ { return heatCapacity_; }
+
+ /*!
+ * \brief Returns the temperature in fluid / solid phase(s)
+ * the sub-control volume.
+ * \param phaseIdx The local index of the phases
+ */
+ Scalar temperature(const unsigned int phaseIdx) const
+ { return temperature_[phaseIdx]; }
+
+ /*!
+ * \brief Returns the total density of the given solid phase [kg / m^3] in
+ * the sub-control volume.
+ */
+ Scalar densitySolid() const
+ { return densitySolid_; }
+
+ /*!
+ * \brief Returns the conductivity of the given solid phase [kg / m^3] in
+ * the sub-control volume.
+ */
+ Scalar thermalConductivitySolid() const
+ { return thermalConductivitySolid_; }
+
+ /*!
+ * \brief Returns the conductivity of the given fluid [kg / m^3] in
+ * the sub-control volume.
+ *
+ * \param phaseIdx The local index of the phases
+ */
+ Scalar thermalConductivity(const unsigned int phaseIdx) const
+ {
+ if(phaseIdx == wPhaseIdx or phaseIdx == nPhaseIdx )
+ return thermalConductivityFluid_[phaseIdx];
+ else if (phaseIdx == sPhaseIdx )
+ return thermalConductivitySolid_;
+ else
+ DUNE_THROW(Dune::NotImplemented,
+ "wrong index");
+ }
+
+ void checkDefinedTemp() const
+ { Valgrind::CheckDefined(temperature_); }
+
+ /*!
+ * \brief If running under valgrind this produces an error message
+ * if some of the object's attributes is undefined.
+ */
+ void checkDefined() const
+ {
+ Valgrind::CheckDefined(temperature_);
+ Valgrind::CheckDefined(thermalConductivityFluid_);
+ Valgrind::CheckDefined(thermalConductivitySolid_);
+ Valgrind::CheckDefined(densitySolid_);
+ Valgrind::CheckDefined(heatCapacity_);
+ }
+
+protected:
+ Scalar temperature_[numEnergyEqs];
+ Scalar heatCapacity_;
+ Scalar densitySolid_;
+ Scalar thermalConductivitySolid_;
+ Scalar thermalConductivityFluid_[numPhases];
+};
+
} // end namespace
#endif
diff --git a/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh b/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh
index 208e97ef44..c4c64d7399 100644
--- a/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh
+++ b/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh
@@ -39,10 +39,10 @@ namespace Dumux
*/
template<class TypeTag,
bool enableEnergy /* = false */,
- bool enableKineticEnergy /* = false */>
+ int numEnergyEquations/*=0*/>
class MPNCVtkWriterEnergy : public MPNCVtkWriterModule<TypeTag>
{
- static_assert(enableKineticEnergy == false,
+ static_assert(numEnergyEquations < 1,
"If you enable kinetic energy transfer between fluids, you"
"also have to enable the energy in general!");
@@ -130,7 +130,7 @@ private:
* local thermal equilibrium. (i.e. no kinetic energy transfer)
*/
template<class TypeTag>
-class MPNCVtkWriterEnergy<TypeTag, /* enableEnergy = */ true, /* enableKineticEnergy = */ false >
+class MPNCVtkWriterEnergy<TypeTag, /* enableEnergy = */ true, /* numEnergyEquations = */ 1 >
: public MPNCVtkWriterModule<TypeTag>
{
typedef MPNCVtkWriterModule<TypeTag> ParentType;
diff --git a/dumux/implicit/mpnc/energy/mpncvtkwriterenergykinetic.hh b/dumux/implicit/mpnc/energy/mpncvtkwriterenergykinetic.hh
index c134c84f5d..ff7581b2b1 100644
--- a/dumux/implicit/mpnc/energy/mpncvtkwriterenergykinetic.hh
+++ b/dumux/implicit/mpnc/energy/mpncvtkwriterenergykinetic.hh
@@ -36,11 +36,11 @@ namespace Dumux
* \brief VTK writer module for the energy related quantities of the
* MpNc model.
*
- * This is the specialization for the case with an energy equation and
+ * This is the specialization for the case with *3* energy balance equations and
* no local thermal equilibrium. (i.e. _including_ kinetic energy transfer)
*/
template<class TypeTag>
-class MPNCVtkWriterEnergy<TypeTag, /*enableEnergy = */ true, /* enableKineticEnergy = */ true >
+class MPNCVtkWriterEnergy<TypeTag, /*enableEnergy = */ true, /* numEnergyEquations = */ 3 >
: public MPNCVtkWriterModule<TypeTag>
{
typedef MPNCVtkWriterModule<TypeTag> ParentType;
@@ -267,6 +267,242 @@ private:
ScalarVector ans_;
};
+/*!
+ * \ingroup MPNCModel
+ *
+ * \brief VTK writer module for the energy related quantities of the
+ * MpNc model.
+ *
+ * This is the specialization for the case with *2* energy balance equations and
+ * no local thermal equilibrium. (i.e. _including_ kinetic energy transfer)
+ */
+template<class TypeTag>
+class MPNCVtkWriterEnergy<TypeTag, /*enableEnergy = */ true, /* numEnergyEquations = */ 2 >
+ : public MPNCVtkWriterModule<TypeTag>
+{
+ typedef MPNCVtkWriterModule<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ enum { dim = GridView::dimension };
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+ enum { numEnergyEqs = Indices::numPrimaryEnergyVars};
+ enum { velocityAveragingInModel = GET_PROP_VALUE(TypeTag, VelocityAveragingInModel) };
+
+ enum { temperatureOutput = GET_PROP_VALUE(TypeTag, VtkAddTemperatures) };
+ enum { enthalpyOutput = GET_PROP_VALUE(TypeTag, VtkAddEnthalpies) };
+ enum { internalEnergyOutput = GET_PROP_VALUE(TypeTag, VtkAddInternalEnergies) };
+ enum { reynoldsOutput = GET_PROP_VALUE(TypeTag, VtkAddReynolds) };
+ enum { prandtlOutput = GET_PROP_VALUE(TypeTag, VtkAddPrandtl) };
+ enum { nusseltOutput = GET_PROP_VALUE(TypeTag, VtkAddNusselt) };
+ enum { interfacialAreaOutput= GET_PROP_VALUE(TypeTag, VtkAddInterfacialArea) };
+ enum { velocityOutput = GET_PROP_VALUE(TypeTag, VtkAddVelocities) };
+
+ typedef typename ParentType::ScalarVector ScalarVector;
+ typedef typename ParentType::PhaseVector PhaseVector;
+ typedef typename ParentType::ComponentVector ComponentVector;
+ typedef Dune::array<ScalarVector, numEnergyEqs> EnergyEqVector;
+
+ typedef Dune::FieldVector<Scalar, dim> DimVector;
+ typedef Dune::BlockVector<DimVector> DimField;
+ typedef Dune::array<DimField, numPhases> PhaseDimField;
+
+
+public:
+ MPNCVtkWriterEnergy(const Problem &problem)
+ : ParentType(problem)
+ {
+ }
+
+ /*!
+ * \brief Allocate memory for the scalar fields we would like to
+ * write to the VTK file.
+ */
+ template <class MultiWriter>
+ void allocBuffers(MultiWriter &writer)
+ {
+ resizeTemperaturesBuffer_(temperature_);
+ this->resizePhaseBuffer_(enthalpy_);
+ this->resizePhaseBuffer_(internalEnergy_);
+ this->resizePhaseBuffer_(reynoldsNumber_);
+ this->resizePhaseBuffer_(prandtlNumber_);
+ this->resizePhaseBuffer_(nusseltNumber_);
+
+
+ if (velocityAveragingInModel and not velocityOutput/*only one of the two output options, otherwise paraview segfaults due to two times the same field name*/) {
+ Scalar numVertices = this->problem_.gridView().size(dim);
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ velocity_[phaseIdx].resize(numVertices);
+ velocity_[phaseIdx] = 0;
+ }
+ }
+ }
+
+ /*!
+ * \brief Modify the internal buffers according to the volume
+ * variables seen on an element
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry in the fully implicit scheme
+ * \param elemVolVars The volume variables of the current element
+ * \param elemBcTypes
+ */
+ void processElement(const Element & element,
+ const FVElementGeometry & fvGeometry,
+ const ElementVolumeVariables & elemVolVars,
+ const ElementBoundaryTypes & elemBcTypes)
+ {
+ int numLocalVertices = element.geometry().corners();
+ for (int localVertexIdx = 0; localVertexIdx < numLocalVertices; ++localVertexIdx) {
+ const unsigned int globalIdx = this->problem_.vertexMapper().map(element, localVertexIdx, dim);
+ const VolumeVariables &volVars = elemVolVars[localVertexIdx];
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ enthalpy_[phaseIdx][globalIdx] = volVars.fluidState().enthalpy(phaseIdx);
+ internalEnergy_[phaseIdx][globalIdx] = volVars.fluidState().internalEnergy(phaseIdx);
+ reynoldsNumber_[phaseIdx][globalIdx] = volVars.reynoldsNumber(phaseIdx);
+ prandtlNumber_[phaseIdx][globalIdx] = volVars.prandtlNumber(phaseIdx);
+ nusseltNumber_[phaseIdx][globalIdx] = volVars.nusseltNumber(phaseIdx);
+ }
+
+ // because numPhases only counts liquid phases
+ for (int phaseIdx = 0; phaseIdx < numEnergyEqs; ++ phaseIdx) {
+ temperature_[phaseIdx][globalIdx] = volVars.temperature(phaseIdx);
+ Valgrind::CheckDefined(temperature_[phaseIdx][globalIdx]);
+ }
+
+ if (velocityAveragingInModel and not velocityOutput/*only one of the two output options, otherwise paraview segfaults due to two times the same field name*/){
+ int numVertices = this->problem_.gridView().size(dim); // numVertices for vertexCentereed, numVolumes for volume centered
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int I = 0; I < numVertices; ++I)
+ velocity_[phaseIdx][I] = this->problem_.model().volumeDarcyVelocity(phaseIdx, I);
+ }
+ }
+ }
+
+ /*!
+ * \brief Add all buffers to the VTK output writer.
+ */
+ template <class MultiWriter>
+ void commitBuffers(MultiWriter & writer)
+ {
+
+ if (temperatureOutput){
+ this->commitTemperaturesBuffer_(writer, "T_%s", temperature_);
+ }
+
+ if (enthalpyOutput)
+ this->commitPhaseBuffer_(writer, "h_%s", enthalpy_);
+ if (internalEnergyOutput)
+ this->commitPhaseBuffer_(writer, "u_%s", internalEnergy_);
+ if (reynoldsOutput)
+ this->commitPhaseBuffer_(writer, "reynoldsNumber_%s", reynoldsNumber_);
+ if (prandtlOutput)
+ this->commitPhaseBuffer_(writer, "prandtlNumber_%s", prandtlNumber_);
+ if (nusseltOutput)
+ this->commitPhaseBuffer_(writer, "nusseltNumber_%s", nusseltNumber_);
+ if (velocityAveragingInModel and not velocityOutput/*only one of the two output options, otherwise paraview segfaults due to two timies the same field name*/){
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ // commit the phase velocity
+ std::ostringstream oss;
+ oss << "velocity_" << FluidSystem::phaseName(phaseIdx);
+ writer.attachVertexData(velocity_[phaseIdx],
+ oss.str(),
+ dim);
+ }
+ }
+ }
+
+private:
+ /*!
+ * \brief Allocate the space for a buffer storing temperatures.
+ * This is one more entry than (fluid) phases.
+ */
+ void resizeTemperaturesBuffer_(EnergyEqVector & buffer,
+ bool vertexCentered = true)
+ {
+ Scalar n; // numVertices for vertexCentereed, numVolumes for volume centered
+ if (vertexCentered)
+ n = this->problem_.gridView().size(dim);
+ else
+ n = this->problem_.gridView().size(0);
+
+ for (int energyEqIdx = 0; energyEqIdx < numEnergyEqs; ++energyEqIdx) {
+ buffer[energyEqIdx].resize(n);
+ std::fill(buffer[energyEqIdx].begin(), buffer[energyEqIdx].end(), 0.0);
+ }
+ }
+
+/*!
+ * \brief Add a buffer for the three tmeperatures (fluids+solid) to the VTK result file.
+ */
+ template <class MultiWriter>
+ void commitTemperaturesBuffer_(MultiWriter & writer,
+ const char *pattern,
+ EnergyEqVector & buffer,
+ bool vertexCentered = true)
+ {
+ static const char *name[] = {
+ "fluid",
+ "solid"
+ };
+
+ for (int energyEqIdx = 0; energyEqIdx < numEnergyEqs; ++energyEqIdx) {
+ std::ostringstream oss;
+ oss << "T_" << name[energyEqIdx];
+
+ if (vertexCentered)
+ writer.attachVertexData(buffer[energyEqIdx], oss.str(), 1);
+ else
+ writer.attachCellData(buffer[energyEqIdx], oss.str(), 1);
+ }
+ }
+
+// static const char *phaseName(int phaseIdx)
+// {
+// static const char *name[] = {
+// "w",
+// "n"
+// };
+//
+// assert(0 <= phaseIdx && phaseIdx < numPhases);
+// return name[phaseIdx];
+// }
+
+ EnergyEqVector temperature_ ;
+ PhaseVector enthalpy_ ;
+ PhaseVector internalEnergy_ ;
+ PhaseVector reynoldsNumber_ ;
+ PhaseVector prandtlNumber_ ;
+ PhaseVector nusseltNumber_ ;
+
+ PhaseDimField velocity_;
+};
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
}
#endif
diff --git a/dumux/implicit/mpnc/mass/mpnclocalresidualmass.hh b/dumux/implicit/mpnc/mass/mpnclocalresidualmass.hh
index edd70b4a77..c2bdf90d59 100644
--- a/dumux/implicit/mpnc/mass/mpnclocalresidualmass.hh
+++ b/dumux/implicit/mpnc/mass/mpnclocalresidualmass.hh
@@ -51,11 +51,11 @@ protected:
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
- enum { enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy) };
+ enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations) };
typedef typename Dune::FieldVector<Scalar, numComponents> ComponentVector;
typedef MPNCDiffusion<TypeTag, enableDiffusion> Diffusion;
- typedef MPNCLocalResidualEnergy<TypeTag, enableEnergy, enableKineticEnergy> EnergyResid;
+ typedef MPNCLocalResidualEnergy<TypeTag, enableEnergy, numEnergyEquations> EnergyResid;
public:
/*!
@@ -78,9 +78,14 @@ public:
storage[compIdx] +=
volVars.fluidState().saturation(phaseIdx)*
volVars.fluidState().molarity(phaseIdx, compIdx);
+#ifndef NDEBUG
+if (!std::isfinite(storage[compIdx]))
+ DUNE_THROW(NumericalProblem, "Calculated non-finite storage");
+#endif
}
storage *= volVars.porosity();
+
}
/*!
@@ -98,10 +103,7 @@ public:
{
const Scalar volumeFlux = fluxVars.volumeFlux(phaseIdx) ;
- #ifndef NDEBUG
- if (!std::isfinite(volumeFlux))
- DUNE_THROW(NumericalProblem, "Calculated non-finite normal flux");
- #endif
+
// retrieve the upwind weight for the mass conservation equations. Use the value
// specified via the property system as default, and overwrite
@@ -118,6 +120,8 @@ public:
const VolumeVariables &up = fluxVars.volVars(upIdx);
const VolumeVariables &dn = fluxVars.volVars(dnIdx);
+if (!std::isfinite(volumeFlux))
+ DUNE_THROW(NumericalProblem, "Calculated non-finite normal flux in phase" << phaseIdx);
////////
// advective fluxes of all components in the phase
////////
@@ -179,9 +183,9 @@ public:
(( massUpwindWeight)*up.fluidState().molarity(phaseIdx, compIdx)
+
( 1. - massUpwindWeight)*dn.fluidState().molarity(phaseIdx, compIdx) );
-
-
- }
+ if (!std::isfinite(flux[compIdx]))
+ DUNE_THROW(NumericalProblem, "Calculated non-finite normal flux in phase " << phaseIdx << " comp " << compIdx << "T: "<< up.fluidState().temperature(phaseIdx) << "S "<<up.fluidState().saturation(phaseIdx) ) ;
+ }
}
}
@@ -250,10 +254,10 @@ class MPNCLocalResidualMass
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
enum { conti0EqIdx = Indices::conti0EqIdx };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
- enum { enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy) };
+ enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations) };
typedef typename Dune::FieldVector<Scalar, numComponents> ComponentVector;
- typedef MPNCLocalResidualEnergy<TypeTag, enableEnergy, enableKineticEnergy> EnergyResid;
+ typedef MPNCLocalResidualEnergy<TypeTag, enableEnergy, numEnergyEquations> EnergyResid;
public:
/*!
@@ -323,15 +327,21 @@ public:
MassCommon::computeDiffusivePhaseFlux(phaseComponentValuesDiffusion, fluxVars, phaseIdx);
Valgrind::CheckDefined(phaseComponentValuesDiffusion);
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx){
flux[conti0EqIdx + compIdx] +=
phaseComponentValuesDiffusion[compIdx];
+ }
+
+
// Right now I think that adding the two contributions individually into the flux is best for debugging and understanding.
// The Energy module needs both contributions.
phaseComponentValuesMassTransport[phaseIdx] = phaseComponentValuesDiffusion + phaseComponentValuesAdvection ;
+
Valgrind::CheckDefined(flux);
}
+// std::cout<< "KPN Mass: flux: " << flux << endl;
+
// \todo
//
@@ -361,10 +371,11 @@ public:
static void computeSource(PrimaryVariables &source,
const VolumeVariables &volVars)
{
- static_assert(not enableKineticEnergy, // enableKinetic is disabled, in this specialization
- "In the case of kinetic energy transfer the advective energy transport between the phases has to be considered. "
- "It is hard (technically) to say how much mass got transfered in the case of chemical equilibrium. "
- "Therefore, kineticEnergy and no kinetic mass does not fit (yet).");
+// static_assert(not enableKineticEnergy, // enableKinetic is disabled, in this specialization
+// "In the case of kinetic energy transfer the advective energy transport between the phases has to be considered. "
+// "It is hard (technically) to say how much mass got transfered in the case of chemical equilibrium. "
+// "Therefore, kineticEnergy and no kinetic mass does not fit (yet).");
+
// mass transfer is not considered in this mass module
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
diff --git a/dumux/implicit/mpnc/mass/mpnclocalresidualmasskinetic.hh b/dumux/implicit/mpnc/mass/mpnclocalresidualmasskinetic.hh
index 4a7534c0fc..61b4dbd1e1 100644
--- a/dumux/implicit/mpnc/mass/mpnclocalresidualmasskinetic.hh
+++ b/dumux/implicit/mpnc/mass/mpnclocalresidualmasskinetic.hh
@@ -26,15 +26,6 @@
#include <dumux/implicit/mpnc/mass/mpnclocalresidualmass.hh>
-// The idea here is to calculate the mass transfer into both phases (i.e. capture the resistivities of both sides of the interface),
-// under the supplementary constrain, that mass fluxes add up to zero. From this the composition of both phases is to be calculated.
-// This way the equality of fluxes into both phases does not have to be assumed, but can be calculated.
-// NOT YET WORKING!
-#if FUNKYMASSTRANSFER
-#include <dumux/material/constraintsolvers/compositionfrommasstransfer.hh>
-#endif
-
-
namespace Dumux
{
/*!
@@ -56,11 +47,6 @@ class MPNCLocalResidualMass<TypeTag, /*enableKinetic=*/true>
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
-#if FUNKYMASSTRANSFER
- // here, we need a constraint solver, which gives me the composition of the phases as determined by the mass transfer constraints
- typedef CompositionFromMassTransfer<Scalar, FluidSystem> NonEquilConstraintSolver;
-#endif
-
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
enum { conti0EqIdx = Indices::conti0EqIdx };
@@ -69,10 +55,10 @@ class MPNCLocalResidualMass<TypeTag, /*enableKinetic=*/true>
enum { wPhaseIdx = FluidSystem::wPhaseIdx} ;
enum { nPhaseIdx = FluidSystem::nPhaseIdx} ;
enum { sPhaseIdx = FluidSystem::sPhaseIdx} ;
- enum { enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy) };
+ enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
- typedef MPNCLocalResidualEnergy<TypeTag, enableEnergy, enableKineticEnergy> EnergyResid;
+ typedef MPNCLocalResidualEnergy<TypeTag, enableEnergy, numEnergyEquations> EnergyResid;
typedef typename Dune::FieldVector<Scalar, numComponents> ComponentVector;
public:
@@ -123,10 +109,8 @@ public:
storage[conti0EqIdx + phaseIdx*numComponents + compIdx]
+= phaseComponentValues[compIdx];
- #ifndef NDEBUG
if (!std::isfinite(storage[conti0EqIdx + phaseIdx*numComponents + compIdx]))
DUNE_THROW(NumericalProblem, "Calculated non-finite storage");
- #endif
}
Valgrind::CheckDefined(storage);
@@ -160,11 +144,8 @@ public:
phaseComponentValuesDiffusion[compIdx];
Valgrind::CheckDefined(flux);
- #ifndef NDEBUG
if (!std::isfinite(flux[conti0EqIdx + phaseIdx*numComponents + compIdx]))
DUNE_THROW(NumericalProblem, "Calculated non-finite flux");
- #endif
-
}
// Right now I think that adding the two contributions individually into the flux is best for debugging and understanding.
@@ -200,93 +181,58 @@ public:
// between phases is realized via source terms there is a
// balance equation for each component in each phase
ComponentVector componentIntoPhaseMassTransfer[numPhases];
-
+#define FUNKYMASSTRANSFER 0
#if FUNKYMASSTRANSFER
- // calculate the mass transfer coefficient
- // Here mass transfer coefficient is everything, of the mass
- // flux from one phase to another except the difference in mole fraction
- Scalar massTransferCoefficient[numPhases][numComponents];
- for(int smallLoopPhaseIdx=0; smallLoopPhaseIdx<numPhases; ++smallLoopPhaseIdx){
- for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
- static_assert(numComponents == 2,
- "this works only for binary mixtures, where D_AB = D_BA in the case of isothermal, isobaric conditions for diffusion");
- const Scalar rhoN = volVars.fluidState().molarDensity(smallLoopPhaseIdx);
- const Scalar D = FluidSystem::binaryDiffusionCoefficient(volVars.fluidState(),
- smallLoopPhaseIdx,
- wCompIdx,
- nCompIdx);
- Scalar Sh(0.) ;
- Sh = volVars.sherwoodNumber(smallLoopPhaseIdx) ;
- massTransferCoefficient[smallLoopPhaseIdx][compIdx] = volVars.interfacialArea(wPhaseIdx, nPhaseIdx)
- * volVars.fluidState().molarDensity(smallLoopPhaseIdx)
- * FluidSystem::binaryDiffusionCoefficient(volVars.fluidState(),
- smallLoopPhaseIdx,
- wCompIdx,
- nCompIdx)// see static assert above
- * volVars.sherwoodNumber(smallLoopPhaseIdx)
- / volVars.characteristicLength();
-
- Valgrind::CheckDefined(massTransferCoefficient[smallLoopPhaseIdx][compIdx]);
- }
- }
-
-
- Scalar xTransferNonEquil[numPhases][numComponents];
- {
- FluidState nonEquilFluidState;
- nonEquilFluidState.assign(volVars.fluidState()) ;
- NonEquilConstraintSolver::solve(nonEquilFluidState, // for storing the non-equilibrium mole fractions determined by mass transfer
- massTransferCoefficient,
- volVars,
- /*setViscosity=*/false,
- /*setEnthalpy=*/false) ;
-
- for(int smallLoopPhaseIdx=0; smallLoopPhaseIdx<numPhases; ++smallLoopPhaseIdx){
- for (int compIdx=0; compIdx< numComponents; ++ compIdx){
- xTransferNonEquil[smallLoopPhaseIdx][compIdx] = nonEquilFluidState.moleFraction(smallLoopPhaseIdx, compIdx);
- }
- }
- }
+ const Scalar mu_nPhaseNCompEquil = volVars.chemicalPotentialEquil(nPhaseIdx, nCompIdx) ; // very 2p2c
+ const Scalar mu_wPhaseWCompEquil = volVars.chemicalPotentialEquil(wPhaseIdx, wCompIdx); // very 2p2c
+ const Scalar mu_wPhaseNComp = volVars.chemicalPotential(wPhaseIdx, nCompIdx) ; // very 2p2c
+ const Scalar mu_nPhaseWComp = volVars.chemicalPotential(nPhaseIdx, wCompIdx); // very 2p2c
+ Valgrind::CheckDefined(mu_nPhaseNCompEquil);
+ Valgrind::CheckDefined(mu_wPhaseWCompEquil);
+ Valgrind::CheckDefined(mu_wPhaseNComp);
+ Valgrind::CheckDefined(mu_nPhaseWComp);
-
-
- for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
- for (int compIdx=0; compIdx<numComponents; ++compIdx){
- const Scalar deltaX = xTransferNonEquil[phaseIdx][compIdx] - volVars.xEquil(phaseIdx, compIdx);
- const Scalar massTransferCoeff = massTransferCoefficient[phaseIdx][compIdx];
- componentIntoPhaseMassTransfer[phaseIdx][compIdx] = (-deltaX*massTransferCoeff );
- }
- }
+ const Scalar characteristicLength = volVars.characteristicLength() ;
+ const Scalar temperature = volVars.fluidState().temperature(wPhaseIdx);
+ const Scalar pn = volVars.fluidState().pressure(nPhaseIdx);
+ const Scalar henry = FluidSystem::henry(temperature) ;
+ const Scalar gradNinWApprox = ( mu_wPhaseNComp - mu_nPhaseNCompEquil) / characteristicLength; // very 2p2c // 1. / henry *
+ const Scalar gradWinNApprox = ( mu_nPhaseWComp - mu_wPhaseWCompEquil) / characteristicLength; // very 2p2c // 1. / pn *
#else // FUNKYMASSTRANSFER
- // diffusion coefficients in wetting phase
- const Scalar diffCoeffNinW = volVars.diffCoeff(wPhaseIdx, wCompIdx, nCompIdx) ;
- Valgrind::CheckDefined(diffCoeffNinW);
- // diffusion coefficients in non-wetting phase
- const Scalar diffCoeffWinN = volVars.diffCoeff(nPhaseIdx, wCompIdx, nCompIdx) ;
- Valgrind::CheckDefined(diffCoeffWinN);
- Scalar phaseDensity[numPhases];
Scalar x[numPhases][numComponents]; // mass fractions in wetting phase
for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
- phaseDensity[phaseIdx] = volVars.fluidState().molarDensity(phaseIdx);
for (int compIdx=0; compIdx< numComponents; ++ compIdx){
x[phaseIdx][compIdx] = volVars.fluidState().moleFraction(phaseIdx, compIdx);
}
}
Valgrind::CheckDefined(x);
-// // "equilibrium" values: calculated in volume variables
+// "equilibrium" values: calculated in volume variables
const Scalar x_wPhaseNCompEquil = volVars.xEquil(wPhaseIdx, nCompIdx) ; // very 2p2c
const Scalar x_nPhaseWCompEquil = volVars.xEquil(nPhaseIdx, wCompIdx); // very 2p2c
Valgrind::CheckDefined(x_wPhaseNCompEquil);
Valgrind::CheckDefined(x_nPhaseWCompEquil);
const Scalar characteristicLength = volVars.characteristicLength() ;
- const Scalar factorMassTransfer = volVars.factorMassTransfer() ;
- const Scalar awn = volVars.interfacialArea(wPhaseIdx, nPhaseIdx);
const Scalar gradNinWApprox = (x[wPhaseIdx][nCompIdx] - x_wPhaseNCompEquil) / characteristicLength; // very 2p2c
const Scalar gradWinNApprox = (x[nPhaseIdx][wCompIdx] - x_nPhaseWCompEquil) / characteristicLength; // very 2p2c
+#endif
+ Scalar phaseDensity[numPhases];
+ for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
+ phaseDensity[phaseIdx] = volVars.fluidState().molarDensity(phaseIdx);
+ }
+
+ // diffusion coefficients in wetting phase
+ const Scalar diffCoeffNinW = volVars.diffCoeff(wPhaseIdx, wCompIdx, nCompIdx) ;
+ Valgrind::CheckDefined(diffCoeffNinW);
+ // diffusion coefficients in non-wetting phase
+ const Scalar diffCoeffWinN = volVars.diffCoeff(nPhaseIdx, wCompIdx, nCompIdx) ;
+ Valgrind::CheckDefined(diffCoeffWinN);
+
+ const Scalar factorMassTransfer = volVars.factorMassTransfer() ;
+ const Scalar awn = volVars.interfacialArea(wPhaseIdx, nPhaseIdx);
const Scalar sherwoodWPhase = volVars.sherwoodNumber(wPhaseIdx);
const Scalar sherwoodNPhase = volVars.sherwoodNumber(nPhaseIdx);
@@ -305,8 +251,6 @@ public:
componentIntoPhaseMassTransfer[nPhaseIdx][wCompIdx] = wCompIntoNPhase;
componentIntoPhaseMassTransfer[wPhaseIdx][wCompIdx] = wCompIntoWPhase;
-#endif
-
#if MASS_TRANSFER_OFF
#warning MASS TRANSFER OFF
for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
@@ -330,10 +274,9 @@ public:
for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
const unsigned int eqIdx = conti0EqIdx + compIdx + phaseIdx*numComponents;
source[eqIdx] += componentIntoPhaseMassTransfer[phaseIdx][compIdx] ;
- #ifndef NDEBUG
- if (!std::isfinite(source[eqIdx]))
- DUNE_THROW(NumericalProblem, "Calculated non-finite source");
- #endif
+
+ if (!std::isfinite(source[eqIdx]))
+ DUNE_THROW(NumericalProblem, "Calculated non-finite source");
}
}
Valgrind::CheckDefined(source);
diff --git a/dumux/implicit/mpnc/mass/mpncvolumevariablesmasskinetic.hh b/dumux/implicit/mpnc/mass/mpncvolumevariablesmasskinetic.hh
index c06ade3f8b..036e020749 100644
--- a/dumux/implicit/mpnc/mass/mpncvolumevariablesmasskinetic.hh
+++ b/dumux/implicit/mpnc/mass/mpncvolumevariablesmasskinetic.hh
@@ -130,6 +130,24 @@ public:
}
}
+ // Setting the equilibrium chemical potential (in a kinetic model not necessarily the same as the actual chemical potential)
+ for(int smallLoopPhaseIdx=0; smallLoopPhaseIdx<numPhases; ++smallLoopPhaseIdx){
+ for (int compIdx=0; compIdx< numComponents; ++ compIdx){
+ chemicalPotentialEquil_[smallLoopPhaseIdx][compIdx] = FluidSystem::chemicalPotential(equilFluidState,
+ smallLoopPhaseIdx,
+ compIdx);
+ }
+ }
+
+ // Setting the actual chemical potential (in a kinetic model not necessarily the same as the equilibrium chemical potential)
+ for(int smallLoopPhaseIdx=0; smallLoopPhaseIdx<numPhases; ++smallLoopPhaseIdx){
+ for (int compIdx=0; compIdx< numComponents; ++ compIdx){
+ chemicalPotential_[smallLoopPhaseIdx][compIdx] = FluidSystem::chemicalPotential(actualFluidState,
+ smallLoopPhaseIdx,
+ compIdx);
+ }
+ }
+
// compute densities of all phases
for(int smallLoopPhaseIdx=0; smallLoopPhaseIdx<numPhases; ++smallLoopPhaseIdx){
const Scalar rho = FluidSystem::density(actualFluidState, paramCache, smallLoopPhaseIdx);
@@ -141,7 +159,8 @@ public:
}
/*!
- * \brief The mole fraction we would have in the case of chemical equilibrium
+ * \brief The mole fraction we would have in the case of chemical equilibrium /
+ * on the interface.
*
* \param phaseIdx The index of the fluid phase
* \param compIdx The local index of the component
@@ -151,6 +170,29 @@ public:
return xEquil_[phaseIdx][compIdx] ;
}
+ /*!
+ * \brief The chemical potential we would have in the case of chemical equilibrium / on the interface
+ *
+ * \param phaseIdx The index of the fluid phase
+ * \param compIdx The local index of the component
+ */
+ const Scalar chemicalPotential(const unsigned int phaseIdx, const unsigned int compIdx) const
+ {
+ return chemicalPotential_[phaseIdx][compIdx] ;
+ }
+
+ /*!
+ * \brief The actual chemical potential we currently have. In the case of non-equilibrium this is not
+ * necessarily equilibrium.
+ *
+ * \param phaseIdx The index of the fluid phase
+ * \param compIdx The local index of the component
+ */
+ const Scalar chemicalPotentialEquil(const unsigned int phaseIdx, const unsigned int compIdx) const
+ {
+ return chemicalPotentialEquil_[phaseIdx][compIdx] ;
+ }
+
/*!
* \brief If running in valgrind this makes sure that all
* quantities in the volume variables are defined.
@@ -158,13 +200,16 @@ public:
void checkDefined() const
{
#if HAVE_VALGRIND && !defined NDEBUG
-// Valgrind::CheckDefined(interfacialArea_); // only include if we do not have the interfacial area from the energy volume variables
Valgrind::CheckDefined(xEquil_);
+ Valgrind::CheckDefined(chemicalPotentialEquil_);
+ Valgrind::CheckDefined(chemicalPotential_);
#endif
}
private:
Scalar xEquil_[numPhases][numComponents];
+ Scalar chemicalPotentialEquil_[numPhases][numComponents];
+ Scalar chemicalPotential_[numPhases][numComponents];
};
} // end namespace
diff --git a/dumux/implicit/mpnc/mass/mpncvtkwritermasskinetic.hh b/dumux/implicit/mpnc/mass/mpncvtkwritermasskinetic.hh
index 5f8f6a95b3..c1fcdf89f0 100644
--- a/dumux/implicit/mpnc/mass/mpncvtkwritermasskinetic.hh
+++ b/dumux/implicit/mpnc/mass/mpncvtkwritermasskinetic.hh
@@ -84,8 +84,13 @@ public:
template <class MultiWriter>
void allocBuffers(MultiWriter &writer)
{
+// this->resizePhaseComponentBuffer_(drivingThingyMu_);
+// this->resizePhaseComponentBuffer_(drivingThingyMole_);
this->resizePhaseComponentBuffer_(moleFracEquil_);
+
+// this->resizePhaseComponentBuffer_(percentageEquilMoleFrac_);
if (deltaPOutput_) this->resizeScalarBuffer_(deltaP_, isBox);
+
}
/*!
@@ -112,6 +117,15 @@ public:
moleFracEquil_[phaseIdx][compIdx][globalIdx] = volVars.xEquil(phaseIdx, compIdx);
}
}
+
+// for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+// for (unsigned int compIdx = 0; compIdx < numComponents; ++ compIdx) {
+// drivingThingyMole_[phaseIdx][compIdx][globalIdx] = volVars.xEquil(phaseIdx, compIdx) - volVars.fluidState().moleFraction(phaseIdx, compIdx);
+// drivingThingyMu_[phaseIdx][compIdx][globalIdx] = volVars.chemicalPotentialEquil(phaseIdx, compIdx) - volVars.chemicalPotential(phaseIdx, compIdx);
+//// percentageEquilMoleFrac_[phaseIdx][compIdx][globalIdx] = (volVars.fluidState().moleFraction(phaseIdx, compIdx)) / volVars.xEquil(phaseIdx, compIdx);
+// }
+// }
+
if (deltaPOutput_) {
deltaP_[globalIdx] = volVars.fluidState().pressure(nPhaseIdx) - 100000.;
}
@@ -124,6 +138,10 @@ public:
template <class MultiWriter>
void commitBuffers(MultiWriter &writer)
{
+// this->commitPhaseComponentBuffer_(writer, "dirivingThingyMole_%s^%s", drivingThingyMole_);
+// this->commitPhaseComponentBuffer_(writer, "dirivingThingyMu_%s^%s", drivingThingyMu_);
+// this->commitPhaseComponentBuffer_(writer, "percentageEquilMoleFrac_%s^%s", percentageEquilMoleFrac_);
+
if(xEquilOutput){
this->commitPhaseComponentBuffer_(writer, "xEquil_%s^%s", moleFracEquil_);
}
@@ -133,6 +151,11 @@ public:
private:
PhaseComponentMatrix moleFracEquil_;
+// PhaseComponentMatrix drivingThingyMole_;
+// PhaseComponentMatrix drivingThingyMu_;
+// PhaseComponentMatrix percentageEquilMoleFrac_;
+
+
ScalarVector deltaP_;
};
diff --git a/dumux/implicit/mpnc/mpncfluxvariables.hh b/dumux/implicit/mpnc/mpncfluxvariables.hh
index 8a81fe410a..f5dcb7005a 100644
--- a/dumux/implicit/mpnc/mpncfluxvariables.hh
+++ b/dumux/implicit/mpnc/mpncfluxvariables.hh
@@ -64,12 +64,12 @@ class MPNCFluxVariables
enum {enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion)};
enum {enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy)};
enum {enableKinetic = GET_PROP_VALUE(TypeTag, EnableKinetic)};
- enum {enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy)};
+ enum {numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations)};
typedef Dune::FieldVector<Scalar, dim> DimVector;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef MPNCFluxVariablesDiffusion<TypeTag, enableDiffusion> FluxVariablesDiffusion;
- typedef MPNCFluxVariablesEnergy<TypeTag, enableEnergy, enableKineticEnergy> FluxVariablesEnergy;
+ typedef MPNCFluxVariablesEnergy<TypeTag, enableEnergy, numEnergyEquations> FluxVariablesEnergy;
public:
/*
diff --git a/dumux/implicit/mpnc/mpncindices.hh b/dumux/implicit/mpnc/mpncindices.hh
index 96e6b96a8d..2320a0481b 100644
--- a/dumux/implicit/mpnc/mpncindices.hh
+++ b/dumux/implicit/mpnc/mpncindices.hh
@@ -57,17 +57,17 @@ struct MPNCIndices :
public MPNCEnergyIndices<BasePVOffset +
MPNCMassIndices<0, TypeTag, GET_PROP_VALUE(TypeTag, EnableKinetic) >::numPrimaryVars,
GET_PROP_VALUE(TypeTag, EnableEnergy),
- GET_PROP_VALUE(TypeTag, EnableKineticEnergy)>
+ GET_PROP_VALUE(TypeTag, NumEnergyEquations)>
{
private:
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
enum { enableKinetic = GET_PROP_VALUE(TypeTag, EnableKinetic) }; //mass transfer
- enum { enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy) }; // energy transfer
+ enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations) }; // energy transfer
enum { numPhases = FluidSystem::numPhases };
typedef MPNCMassIndices<BasePVOffset, TypeTag, enableKinetic> MassIndices;
- typedef MPNCEnergyIndices<BasePVOffset + MassIndices::numPrimaryVars, enableEnergy, enableKineticEnergy> EnergyIndices;
+ typedef MPNCEnergyIndices<BasePVOffset + MassIndices::numPrimaryVars, enableEnergy, numEnergyEquations> EnergyIndices;
public:
/*!
diff --git a/dumux/implicit/mpnc/mpnclocalresidual.hh b/dumux/implicit/mpnc/mpnclocalresidual.hh
index c88a7748ae..fd7ccd0f99 100644
--- a/dumux/implicit/mpnc/mpnclocalresidual.hh
+++ b/dumux/implicit/mpnc/mpnclocalresidual.hh
@@ -53,12 +53,13 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
enum {numPhases = GET_PROP_VALUE(TypeTag, NumPhases)};
+ enum {numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
enum {numEq = GET_PROP_VALUE(TypeTag, NumEq)};
enum {enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy)};
- enum {enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy)};
+ enum {numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations)};
enum {enableKinetic = GET_PROP_VALUE(TypeTag, EnableKinetic)};
enum {phase0NcpIdx = Indices::phase0NcpIdx};
-
+
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
@@ -69,7 +70,7 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef MPNCLocalResidualEnergy<TypeTag, enableEnergy, enableKineticEnergy> EnergyResid;
+ typedef MPNCLocalResidualEnergy<TypeTag, enableEnergy, numEnergyEquations> EnergyResid;
typedef MPNCLocalResidualMass<TypeTag, enableKinetic> MassResid;
public:
@@ -238,16 +239,16 @@ public:
*/
void eval(const Element &element,
const FVElementGeometry &fvGeometry,
- const ElementVolumeVariables &prevVolVars,
- const ElementVolumeVariables &curVolVars,
+ const ElementVolumeVariables &prevElemVolVars,
+ const ElementVolumeVariables &curElemVolVars,
const ElementBoundaryTypes &bcType)
{
ParentType::eval(element,
fvGeometry,
- prevVolVars,
- curVolVars,
+ prevElemVolVars,
+ curElemVolVars,
bcType);
-
+
if (GET_PROP_VALUE(TypeTag, ImplicitIsBox)
|| !bcType.hasDirichlet())
{
diff --git a/dumux/implicit/mpnc/mpncmodel.hh b/dumux/implicit/mpnc/mpncmodel.hh
index 00344c4b5c..244a8bb693 100644
--- a/dumux/implicit/mpnc/mpncmodel.hh
+++ b/dumux/implicit/mpnc/mpncmodel.hh
@@ -118,11 +118,17 @@ class MPNCModel : public GET_PROP_TYPE(TypeTag, BaseModel)
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef Dumux::MPNCVtkWriter<TypeTag> MPNCVtkWriter;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
enum {enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy)};
enum {enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion)};
enum {enableKinetic = GET_PROP_VALUE(TypeTag, EnableKinetic)};
- enum {enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy)};
+ enum {numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations)};
enum {enableSmoothUpwinding = GET_PROP_VALUE(TypeTag, ImplicitEnableSmoothUpwinding)};
enum {enablePartialReassemble = GET_PROP_VALUE(TypeTag, ImplicitEnablePartialReassemble)};
enum {enableJacobianRecycling = GET_PROP_VALUE(TypeTag, ImplicitEnableJacobianRecycling)};
@@ -130,6 +136,12 @@ class MPNCModel : public GET_PROP_TYPE(TypeTag, BaseModel)
enum {numPhases = GET_PROP_VALUE(TypeTag, NumPhases)};
enum {numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
enum {numEq = GET_PROP_VALUE(TypeTag, NumEq)};
+ enum { numEnergyEqs = Indices::numPrimaryEnergyVars};
+ enum { dimWorld = GridView::dimensionworld};
+ enum { dim = GridView::dimension};
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
public:
void init(Problem &problem)
@@ -144,7 +156,7 @@ public:
<< " components: " << numComponents << "\n"
<< " equations: " << numEq << "\n"
<< " kinetic mass transfer: " << enableKinetic<< "\n"
- << " kinetic energy transfer: " << enableKineticEnergy<< "\n"
+ << " number of energy equations: " << numEnergyEquations<< "\n"
<< " diffusion: " << enableDiffusion << "\n"
<< " energy equation: " << enableEnergy << "\n"
<< " smooth upwinding: " << enableSmoothUpwinding << "\n"
diff --git a/dumux/implicit/mpnc/mpncmodelkinetic.hh b/dumux/implicit/mpnc/mpncmodelkinetic.hh
index 40a2d860d4..2f8fb91993 100644
--- a/dumux/implicit/mpnc/mpncmodelkinetic.hh
+++ b/dumux/implicit/mpnc/mpncmodelkinetic.hh
@@ -63,7 +63,7 @@ class MPNCModelKinetic : public MPNCModel<TypeTag>
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy)};
enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion)};
enum { enableKinetic = GET_PROP_VALUE(TypeTag, EnableKinetic)};
- enum { enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy)};
+ enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations)};
enum { enableSmoothUpwinding = GET_PROP_VALUE(TypeTag, ImplicitEnableSmoothUpwinding)};
enum { enablePartialReassemble = GET_PROP_VALUE(TypeTag, ImplicitEnablePartialReassemble)};
enum { enableJacobianRecycling = GET_PROP_VALUE(TypeTag, ImplicitEnableJacobianRecycling)};
@@ -193,110 +193,126 @@ public:
}// end all phases
}// end calcVelocity
- /*!
- * \brief Check whether the current solution makes sense.
- */
- void checkPlausibility() const
- {
- // Looping over all elements of the domain
- ElementIterator eEndIt = this->problem_().gridView().template end<0>();
- for (ElementIterator eIt = this->problem_().gridView().template begin<0>() ; eIt not_eq eEndIt; ++eIt)
- {
- ElementVolumeVariables elemVolVars;
- FVElementGeometry fvGeometry;
-
- // updating the volume variables
- fvGeometry.update(this->problem_().gridView(), *eIt);
- elemVolVars.update(this->problem_(), *eIt, fvGeometry, false);
-
- std::stringstream message ;
- // number of scv
- const unsigned int numScv = fvGeometry.numScv; // box: numSCV, cc:1
-
- for (unsigned int scvIdx = 0; scvIdx < numScv; ++scvIdx) {
-
- const FluidState & fluidState = elemVolVars[scvIdx].fluidState();
-
- // energy check
- for(unsigned int energyEqIdx=0; energyEqIdx<numEnergyEqs; energyEqIdx++){
- const Scalar eps = 1e-6 ;
- const Scalar temperatureTest = elemVolVars[scvIdx].temperature(energyEqIdx);
- if (not std::isfinite(temperatureTest) or temperatureTest < 0.-eps ){
- message <<"\nUnphysical Value in Energy: \n";
- message << "\tT" <<"_"<<FluidSystem::phaseName(energyEqIdx)<<"="<< temperatureTest <<"\n";
- }
- }
-
- // mass Check
- for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++){
- const Scalar eps = 1e-6 ;
- for (int compIdx=0; compIdx< numComponents; ++ compIdx){
- const Scalar xTest = fluidState.moleFraction(phaseIdx, compIdx);
- if (not std::isfinite(xTest) or xTest < 0.-eps or xTest > 1.+eps ){
- message <<"\nUnphysical Value in Mass: \n";
-
- message << "\tx" <<"_"<<FluidSystem::phaseName(phaseIdx)
- <<"^"<<FluidSystem::componentName(compIdx)<<"="
- << fluidState.moleFraction(phaseIdx, compIdx) <<"\n";
- }
- }
- }
-
- // interfacial area check (interfacial area between fluid as well as solid phases)
- for(int phaseIdxI=0; phaseIdxI<numPhases+1; phaseIdxI++){
- const Scalar eps = 1e-6 ;
- for (int phaseIdxII=0; phaseIdxII< numPhases+1; ++ phaseIdxII){
- if (phaseIdxI == phaseIdxII)
- continue;
- assert(numEnergyEqs == 3) ; // otherwise this ia call does not make sense
- const Scalar ia = elemVolVars[scvIdx].interfacialArea(phaseIdxI, phaseIdxII);
- if (not std::isfinite(ia) or ia < 0.-eps ) {
- message <<"\nUnphysical Value in interfacial area: \n";
- message << "\tia" <<FluidSystem::phaseName(phaseIdxI)
- <<FluidSystem::phaseName(phaseIdxII)<<"="
- << ia << "\n" ;
- message << "\t S[0]=" << fluidState.saturation(0);
- message << "\t S[1]=" << fluidState.saturation(1);
- message << "\t p[0]=" << fluidState.pressure(0);
- message << "\t p[1]=" << fluidState.pressure(1);
- }
- }
- }
-
- // General Check
- for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++){
- const Scalar eps = 1e-6 ;
- const Scalar saturationTest = fluidState.saturation(phaseIdx);
- if (not std::isfinite(saturationTest) or saturationTest< 0.-eps or saturationTest > 1.+eps ){
- message <<"\nUnphysical Value in Saturation: \n";
- message << "\tS" <<"_"<<FluidSystem::phaseName(phaseIdx)<<"=" << std::scientific
- << fluidState.saturation(phaseIdx) << std::fixed << "\n";
- }
- }
-
- // Some check wrote into the error-message, add some additional information and throw
- if (not message.str().empty()){
- // Getting the spatial coordinate
- const GlobalPosition & globalPosCurrent = fvGeometry.subContVol[scvIdx].global;
- std::stringstream positionString ;
-
- // Add physical location
- positionString << "Here:";
- for(int i=0; i<dim; i++)
- positionString << " x"<< (i+1) << "=" << globalPosCurrent[i] << " " ;
- message << "Unphysical value found! \n" ;
- message << positionString.str() ;
- message << "\n";
-
- message << " Here come the primary Variables:" << "\n" ;
- for(unsigned int priVarIdx =0 ; priVarIdx<numEq; ++priVarIdx){
- message << "priVar[" << priVarIdx << "]=" << elemVolVars[scvIdx].priVar(priVarIdx) << "\n";
- }
- DUNE_THROW(NumericalProblem, message.str());
- }
- } // end scv-loop
- } // end element loop
- }
+// /*!
+// * \brief Check whether the current solution makes sense.
+// */
+// void checkPlausibility() const
+// {
+// // Looping over all elements of the domain
+// ElementIterator eEndIt = this->problem_().gridView().template end<0>();
+// for (ElementIterator eIt = this->problem_().gridView().template begin<0>() ; eIt not_eq eEndIt; ++eIt)
+// {
+// ElementVolumeVariables elemVolVars;
+// FVElementGeometry fvGeometry;
+//
+// // updating the volume variables
+// fvGeometry.update(this->problem_().gridView(), *eIt);
+// elemVolVars.update(this->problem_(), *eIt, fvGeometry, false);
+//
+// std::stringstream message ;
+// // number of scv
+// const unsigned int numScv = fvGeometry.numScv; // box: numSCV, cc:1
+//
+// for (unsigned int scvIdx = 0; scvIdx < numScv; ++scvIdx) {
+//
+// const FluidState & fluidState = elemVolVars[scvIdx].fluidState();
+//
+// // energy check
+// for(unsigned int energyEqIdx=0; energyEqIdx<numEnergyEqs; energyEqIdx++){
+// const Scalar eps = 1e-6 ;
+//// const Scalar temperatureTest = elemVolVars[scvIdx].fluidState().temperature();
+// const Scalar temperatureTest = elemVolVars[scvIdx].temperature(energyEqIdx);
+//// const Scalar temperatureTest = 42;
+//
+// if (not std::isfinite(temperatureTest) or temperatureTest < 0. ){
+// message <<"\nUnphysical Value in Energy: \n";
+// message << "\tT" <<"_"<<FluidSystem::phaseName(energyEqIdx)<<"="<< temperatureTest <<"\n";
+// }
+// }
+//
+// // mass Check
+// for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++){
+// const Scalar eps = 1e-6 ;
+// for (int compIdx=0; compIdx< numComponents; ++ compIdx){
+// const Scalar xTest = fluidState.moleFraction(phaseIdx, compIdx);
+// if (not std::isfinite(xTest) or xTest < 0.-eps or xTest > 1.+eps ){
+// message <<"\nUnphysical Value in Mass: \n";
+//
+// message << "\tx" <<"_"<<FluidSystem::phaseName(phaseIdx)
+// <<"^"<<FluidSystem::componentName(compIdx)<<"="
+// << fluidState.moleFraction(phaseIdx, compIdx) <<"\n";
+// }
+// }
+// }
+//
+// // interfacial area check (interfacial area between fluid as well as solid phases)
+// for(int phaseIdxI=0; phaseIdxI<numPhases+1; phaseIdxI++){
+// const Scalar eps = 1e-6 ;
+// for (int phaseIdxII=0; phaseIdxII< numPhases+1; ++ phaseIdxII){
+// if (phaseIdxI == phaseIdxII)
+// continue;
+// assert(numEnergyEqs == 3) ; // otherwise this ia call does not make sense
+// const Scalar ia = elemVolVars[scvIdx].interfacialArea(phaseIdxI, phaseIdxII);
+// if (not std::isfinite(ia) or ia < 0.-eps ) {
+// message <<"\nUnphysical Value in interfacial area: \n";
+// message << "\tia" <<FluidSystem::phaseName(phaseIdxI)
+// <<FluidSystem::phaseName(phaseIdxII)<<"="
+// << ia << "\n" ;
+// message << "\t S[0]=" << fluidState.saturation(0);
+// message << "\t S[1]=" << fluidState.saturation(1);
+// message << "\t p[0]=" << fluidState.pressure(0);
+// message << "\t p[1]=" << fluidState.pressure(1);
+// }
+// }
+// }
+//
+// // General Check
+// for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++){
+// const Scalar eps = 1e-6 ;
+// const Scalar saturationTest = fluidState.saturation(phaseIdx);
+// if (not std::isfinite(saturationTest) or saturationTest< 0.-eps or saturationTest > 1.+eps ){
+// message <<"\nUnphysical Value in Saturation: \n";
+// message << "\tS" <<"_"<<FluidSystem::phaseName(phaseIdx)<<"=" << std::scientific
+// << fluidState.saturation(phaseIdx) << std::fixed << "\n";
+// }
+// }
+//
+// // velocity Check
+// const unsigned int globalVertexIdx = this->problem_().vertexMapper().map(*eIt, scvIdx, dim);
+// for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++){
+// const Scalar eps = 1e-6 ;
+// const Scalar velocityTest = volumeDarcyMagVelocity(phaseIdx, globalVertexIdx);;
+// if (not std::isfinite(velocityTest) ){
+// message <<"\nUnphysical Value in Velocity: \n";
+// message << "\tv" <<"_"<<FluidSystem::phaseName(phaseIdx)<<"=" << std::scientific
+// << velocityTest << std::fixed << "\n";
+// }
+// }
+//
+//
+// // Some check wrote into the error-message, add some additional information and throw
+// if (not message.str().empty()){
+// // Getting the spatial coordinate
+// const GlobalPosition & globalPosCurrent = fvGeometry.subContVol[scvIdx].global;
+// std::stringstream positionString ;
+//
+// // Add physical location
+// positionString << "Here:";
+// for(int i=0; i<dim; i++)
+// positionString << " x"<< (i+1) << "=" << globalPosCurrent[i] << " " ;
+// message << "Unphysical value found! \n" ;
+// message << positionString.str() ;
+// message << "\n";
+//
+// message << " Here come the primary Variables:" << "\n" ;
+// for(unsigned int priVarIdx =0 ; priVarIdx<numEq; ++priVarIdx){
+// message << "priVar[" << priVarIdx << "]=" << elemVolVars[scvIdx].priVar(priVarIdx) << "\n";
+// }
+// DUNE_THROW(NumericalProblem, message.str());
+// }
+// } // end scv-loop
+// } // end element loop
+// }
/*!
* \brief Access to the averaged (magnitude of) velocity for each vertex.
diff --git a/dumux/implicit/mpnc/mpncproperties.hh b/dumux/implicit/mpnc/mpncproperties.hh
index 5ca631ddc4..db25803211 100644
--- a/dumux/implicit/mpnc/mpncproperties.hh
+++ b/dumux/implicit/mpnc/mpncproperties.hh
@@ -110,8 +110,8 @@ NEW_PROP_TAG(EnableDiffusion);
//! Enable kinetic resolution of mass transfer processes?
NEW_PROP_TAG(EnableKinetic);
-//! Enable kinetic resolution of energy transfer processes?
-NEW_PROP_TAG(EnableKineticEnergy);
+//! Property for the definition of the number of energy equations (0,1,2,3)
+NEW_PROP_TAG(NumEnergyEquations);
//! Enable Maxwell Diffusion? (If false: use Fickian Diffusion) Maxwell incorporated the mutual
//! influences of multiple diffusing components. However, Fick seems to be more robust.
diff --git a/dumux/implicit/mpnc/mpncpropertieskinetic.hh b/dumux/implicit/mpnc/mpncpropertieskinetic.hh
index 8beb452418..84013280c6 100644
--- a/dumux/implicit/mpnc/mpncpropertieskinetic.hh
+++ b/dumux/implicit/mpnc/mpncpropertieskinetic.hh
@@ -55,8 +55,8 @@ NEW_PROP_TAG(MPNCModelKinetic);
//! Enable kinetic resolution of mass transfer processes?
NEW_PROP_TAG(EnableKinetic);
-//! Enable kinetic resolution of energy transfer processes?
-NEW_PROP_TAG(EnableKineticEnergy);
+//! Property for the definition of the number of energy equations (0,1,2,3)
+NEW_PROP_TAG(NumEnergyEquations);
//! average the velocity in the model
NEW_PROP_TAG(VelocityAveragingInModel);
diff --git a/dumux/implicit/mpnc/mpncpropertydefaults.hh b/dumux/implicit/mpnc/mpncpropertydefaults.hh
index 5f6f539237..2434da4859 100644
--- a/dumux/implicit/mpnc/mpncpropertydefaults.hh
+++ b/dumux/implicit/mpnc/mpncpropertydefaults.hh
@@ -144,7 +144,7 @@ SET_BOOL_PROP(MPNC, EnableDiffusion, false);
SET_BOOL_PROP(MPNC, EnableKinetic, false);
//! disable kinetic energy transfer by default
-SET_BOOL_PROP(MPNC, EnableKineticEnergy, false);
+SET_INT_PROP(MPNC, NumEnergyEquations, 0);
//! disable Maxwell Diffusion by default: use Fick
SET_BOOL_PROP(MPNC, UseMaxwellDiffusion, false);
@@ -197,9 +197,9 @@ SET_PROP(MPNC, MPNCVtkEnergyModule)
{
private:
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
- enum { enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy) };
+ enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations) };
public:
- typedef MPNCVtkWriterEnergy<TypeTag, enableEnergy, enableKineticEnergy> type;
+ typedef MPNCVtkWriterEnergy<TypeTag, enableEnergy, numEnergyEquations> type;
};
//! Somerton is used as default model to compute the effective thermal heat conductivity
diff --git a/dumux/implicit/mpnc/mpncvolumevariables.hh b/dumux/implicit/mpnc/mpncvolumevariables.hh
index 73e59f0228..61f138d95d 100644
--- a/dumux/implicit/mpnc/mpncvolumevariables.hh
+++ b/dumux/implicit/mpnc/mpncvolumevariables.hh
@@ -45,10 +45,10 @@ namespace Dumux
template <class TypeTag>
class MPNCVolumeVariables
: public ImplicitVolumeVariables<TypeTag>
- , public MPNCVolumeVariablesIA<TypeTag, GET_PROP_VALUE(TypeTag, EnableKinetic), GET_PROP_VALUE(TypeTag, EnableKineticEnergy)>
+ , public MPNCVolumeVariablesIA<TypeTag, GET_PROP_VALUE(TypeTag, EnableKinetic), GET_PROP_VALUE(TypeTag, NumEnergyEquations)>
, public MPNCVolumeVariablesMass<TypeTag, GET_PROP_VALUE(TypeTag, EnableKinetic)>
, public MPNCVolumeVariablesDiffusion<TypeTag, GET_PROP_VALUE(TypeTag, EnableDiffusion) || GET_PROP_VALUE(TypeTag, EnableKinetic)>
- , public MPNCVolumeVariablesEnergy<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergy), GET_PROP_VALUE(TypeTag, EnableKineticEnergy)>
+ , public MPNCVolumeVariablesEnergy<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergy), GET_PROP_VALUE(TypeTag, NumEnergyEquations)>
{
typedef ImplicitVolumeVariables<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
@@ -76,15 +76,15 @@ class MPNCVolumeVariables
enum {numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
enum {enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy)};
enum {enableKinetic = GET_PROP_VALUE(TypeTag, EnableKinetic)};
- enum {enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy)};
+ enum {numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations)};
enum {enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) || enableKinetic};
enum {s0Idx = Indices::s0Idx};
enum {p0Idx = Indices::p0Idx};
typedef typename GridView::template Codim<0>::Entity Element;
typedef MPNCVolumeVariablesMass<TypeTag, enableKinetic> MassVolumeVariables;
- typedef MPNCVolumeVariablesEnergy<TypeTag, enableEnergy, enableKineticEnergy> EnergyVolumeVariables;
- typedef MPNCVolumeVariablesIA<TypeTag, enableKinetic, enableKineticEnergy> IAVolumeVariables;
+ typedef MPNCVolumeVariablesEnergy<TypeTag, enableEnergy, numEnergyEquations> EnergyVolumeVariables;
+ typedef MPNCVolumeVariablesIA<TypeTag, enableKinetic, numEnergyEquations> IAVolumeVariables;
typedef MPNCVolumeVariablesDiffusion<TypeTag, enableDiffusion> DiffusionVolumeVariables;
public:
@@ -214,6 +214,7 @@ public:
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
// viscosities
Scalar mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
+
fluidState_.setViscosity(phaseIdx, mu);
}
diff --git a/dumux/implicit/mpnc/mpncvolumevariablesia.hh b/dumux/implicit/mpnc/mpncvolumevariablesia.hh
index c49ef0d7db..7c45da6c56 100644
--- a/dumux/implicit/mpnc/mpncvolumevariablesia.hh
+++ b/dumux/implicit/mpnc/mpncvolumevariablesia.hh
@@ -39,10 +39,10 @@ namespace Dumux
* This is the specialization for the cases which do _not_ require
* specific interfacial area.
*/
-template <class TypeTag, bool enableKinetic /* = false */, bool enableKineticEnergy /* = false */>
+template <class TypeTag, bool enableKinetic /* = false */, bool numEnergyEquations /* = false */>
class MPNCVolumeVariablesIA
{
- static_assert(not enableKinetic and not enableKineticEnergy,
+ static_assert(not enableKinetic and not numEnergyEquations,
"The kinetic energy modules need specific interfacial area "
"but no suitable specialization of the IA volume variables module "
"has been included.");
diff --git a/dumux/implicit/mpnc/mpncvolumevariablesiakinetic.hh b/dumux/implicit/mpnc/mpncvolumevariablesiakinetic.hh
index d802f52160..279fccf8d3 100644
--- a/dumux/implicit/mpnc/mpncvolumevariablesiakinetic.hh
+++ b/dumux/implicit/mpnc/mpncvolumevariablesiakinetic.hh
@@ -41,7 +41,7 @@ namespace Dumux
// specialization for the case of kinetic mass AND energy transfer
////////////////////////////////////////////////////////////////////////////////////////////////////
template <class TypeTag, bool enableKinetic >
-class MPNCVolumeVariablesIA<TypeTag, enableKinetic, /*bool enableKineticEnergy=*/true>
+class MPNCVolumeVariablesIA<TypeTag, enableKinetic, /*numEnergyEqs=*/3>
{
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
@@ -69,6 +69,8 @@ class MPNCVolumeVariablesIA<TypeTag, enableKinetic, /*bool enableKineticEnergy=*
enum { numEnergyEqs = Indices::numPrimaryEnergyVars};
enum { nusseltFormulation = GET_PROP_VALUE(TypeTag, NusseltFormulation)} ;
enum { sherwoodFormulation = GET_PROP_VALUE(TypeTag, SherwoodFormulation)} ;
+ enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion)} ;
+
typedef DimensionlessNumbers<Scalar> DimLessNum ;
@@ -223,9 +225,11 @@ if (AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams,
const Scalar heatCapacity = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
const Scalar thermalConductivity = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
+ // If Diffusion is not enabled, Sherwood is divided by zero
+ assert( enableDiffusion ) ;
+
// diffusion coefficient of non-wetting component in wetting phase
const Scalar diffCoeff = volVars.diffCoeff(phaseIdx, wCompIdx, nCompIdx) ;
-
const Scalar porosity = problem.spatialParams().porosity(element,
fvGeometry,
scvIdx);
@@ -332,11 +336,149 @@ private:
};
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// specialization for the case of NO kinetic mass but kinteic energy transfer of a fluid mixture and solid
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////
+template <class TypeTag>
+class MPNCVolumeVariablesIA<TypeTag, /*enableKinetic=*/false, /*numEnergyEqs=*/2>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename FluidSystem::ParameterCache ParameterCache;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+ enum { dim = GridView::dimension};
+ enum { dimWorld = GridView::dimensionworld};
+ enum { numEnergyEqs = Indices::numPrimaryEnergyVars};
+ enum { nusseltFormulation = GET_PROP_VALUE(TypeTag, NusseltFormulation)} ;
+
+ typedef DimensionlessNumbers<Scalar> DimLessNum ;
+public:
+ /*!
+ * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
+ */
+ void update(const VolumeVariables & volVars,
+ const FluidState & fluidState,
+ const ParameterCache ¶mCache,
+ const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element & element,
+ const FVElementGeometry & fvGeometry,
+ const unsigned int scvIdx)
+ {
+ factorMassTransfer_ = problem.spatialParams().factorMassTransfer(element,
+ fvGeometry,
+ scvIdx);
+
+ factorEnergyTransfer_ = problem.spatialParams().factorEnergyTransfer(element,
+ fvGeometry,
+ scvIdx);
+
+ characteristicLength_ = problem.spatialParams().characteristicLength(element,
+ fvGeometry,
+ scvIdx);
+
+ // setting the dimensionless numbers.
+ // obtaining the respective quantities.
+ const unsigned int globalVertexIdx = problem.vertexMapper().map(element, scvIdx, dim);
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ const Scalar darcyMagVelocity = problem.model().volumeDarcyMagVelocity(phaseIdx, globalVertexIdx);
+ const Scalar dynamicViscosity = fluidState.viscosity(phaseIdx);
+ const Scalar density = fluidState.density(phaseIdx);
+ const Scalar kinematicViscosity = dynamicViscosity / density;
+ const Scalar heatCapacity = FluidSystem::heatCapacity(fluidState,
+ paramCache,
+ phaseIdx);
+ const Scalar thermalConductivity = FluidSystem::thermalConductivity(fluidState,
+ paramCache,
+ phaseIdx);
+
+ const Scalar porosity = problem.spatialParams().porosity(element,
+ fvGeometry,
+ scvIdx);
+
+ reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity,
+ characteristicLength_,
+ kinematicViscosity);
+
+ prandtlNumber_[phaseIdx] = DimLessNum::prandtlNumber(dynamicViscosity,
+ heatCapacity,
+ thermalConductivity);
+
+
+ nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
+ prandtlNumber_[phaseIdx],
+ porosity,
+ nusseltFormulation);
+ }
+ }
+
+
+ //! access function Reynolds Number
+ const Scalar reynoldsNumber(const unsigned int phaseIdx) const
+ { return reynoldsNumber_[phaseIdx]; }
+
+ //! access function Prandtl Number
+ const Scalar prandtlNumber(const unsigned int phaseIdx) const
+ { return prandtlNumber_[phaseIdx]; }
+
+ //! access function Nusselt Number
+ const Scalar nusseltNumber(const unsigned int phaseIdx) const
+ { return nusseltNumber_[phaseIdx]; }
+
+ //! access function characteristic length
+ const Scalar characteristicLength() const
+ { return characteristicLength_; }
+
+ //! access function pre factor energy transfer
+ const Scalar factorEnergyTransfer() const
+ { return factorEnergyTransfer_; }
+
+ //! access function pre factor mass transfer
+ const Scalar factorMassTransfer() const
+ { return factorMassTransfer_; }
+
+ /*!
+ * \brief If running in valgrind this makes sure that all
+ * quantities in the volume variables are defined.
+ */
+ void checkDefined() const
+ {
+#if !defined NDEBUG && HAVE_VALGRIND
+ Valgrind::CheckDefined(reynoldsNumber_);
+ Valgrind::CheckDefined(prandtlNumber_);
+ Valgrind::CheckDefined(nusseltNumber_);
+ Valgrind::CheckDefined(characteristicLength_);
+ Valgrind::CheckDefined(factorEnergyTransfer_);
+ Valgrind::CheckDefined(factorMassTransfer_);
+#endif
+ }
+
+private:
+ //! dimensionless numbers
+ Scalar reynoldsNumber_[numPhases];
+ Scalar prandtlNumber_[numPhases];
+ Scalar nusseltNumber_[numPhases];
+ Scalar characteristicLength_;
+ Scalar factorEnergyTransfer_;
+ Scalar factorMassTransfer_;
+ Scalar solidSurface_ ;
+};
+
+
////////////////////////////////////////////////////////////////////////////////////////////////////
// specialization for the case of (only) kinetic mass transfer
////////////////////////////////////////////////////////////////////////////////////////////////////
template <class TypeTag>
-class MPNCVolumeVariablesIA<TypeTag, /*enableKinetic=*/true, /*bool enableKineticEnergy=*/false>
+class MPNCVolumeVariablesIA<TypeTag, /*enableKinetic=*/true, /*numEnergyEqs=*/0>
{
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
diff --git a/dumux/material/fluidmatrixinteractions/2p/thermalconductivitysimplefluidlumping.hh b/dumux/material/fluidmatrixinteractions/2p/thermalconductivitysimplefluidlumping.hh
new file mode 100644
index 0000000000..32da03ac8a
--- /dev/null
+++ b/dumux/material/fluidmatrixinteractions/2p/thermalconductivitysimplefluidlumping.hh
@@ -0,0 +1,81 @@
+// -*- 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 THERMALCONDUCTIVITY_SIMPLE_FLUID_LUMPING_HH
+#define THERMALCONDUCTIVITY_SIMPLE_FLUID_LUMPING_HH
+
+#include <algorithm>
+
+namespace Dumux
+{
+/*!
+ * \ingroup fluidmatrixinteractionslaws
+ *
+ */
+template<class TypeTag>
+class ThermalConductivitySimpleFluidLumping
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations)};
+
+
+public:
+ /*!
+ * \brief Returns the effective thermal conductivity \f$[W/(m K)]\f$.
+ *
+ * \param sw The saturation of the wetting phase
+ * \param lambdaW the thermal conductivity of the wetting phase
+ * \param lambdaN the thermal conductivity of the non-wetting phase
+ * \param lambdaSolid the thermal conductivity of the solid phase
+ * \param porosity The porosity
+ *
+ * \return Effective thermal conductivity of the fluid phases
+ */
+ static Scalar effectiveThermalConductivity(const Scalar sw,
+ const Scalar lambdaW,
+ const Scalar lambdaN,
+ const Scalar lambdaSolid,
+ const Scalar porosity)
+ {
+ assert(numEnergyEquations != 3) ;
+
+ // Franz Lindner / Shi & Wang 2011
+ const Scalar satW = std::max<Scalar>(0.0, sw);
+
+ const Scalar kfeff = porosity *((1.-satW)*lambdaN + satW*lambdaW) ; // arithmetic
+
+ Scalar keff ;
+
+ if ( numEnergyEquations == 2){ // solid dealed with individually (extra balance equation)
+ keff = kfeff ;
+ }
+ else{
+ const Scalar kseff = (1.0-porosity) * lambdaSolid ;
+ keff = kfeff + kseff;
+ }
+
+ return keff ;
+ }
+};
+}
+#endif
diff --git a/dumux/material/fluidmatrixinteractions/2p/vangenuchtenoftemperature.hh b/dumux/material/fluidmatrixinteractions/2p/vangenuchtenoftemperature.hh
index 4f00adba50..aedd98cfff 100644
--- a/dumux/material/fluidmatrixinteractions/2p/vangenuchtenoftemperature.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/vangenuchtenoftemperature.hh
@@ -43,6 +43,7 @@ template <class ScalarT, class ParamsT = RegularizedVanGenuchtenParams<ScalarT>
class RegularizedVanGenuchtenOfTemperature : public Dumux::RegularizedVanGenuchten<ScalarT, ParamsT>
{
typedef Dumux::RegularizedVanGenuchten<ScalarT, ParamsT> RegularizedVanGenuchten;
+ // Data is in /home/pnuske/paper/pcOfT/
public:
typedef ParamsT Params;
diff --git a/test/implicit/mpnc/evaporationatmosphereproblem.hh b/test/implicit/mpnc/evaporationatmosphereproblem.hh
index 30c2de4e8a..78e6691ed8 100644
--- a/test/implicit/mpnc/evaporationatmosphereproblem.hh
+++ b/test/implicit/mpnc/evaporationatmosphereproblem.hh
@@ -62,9 +62,9 @@
#include <dumux/io/plotoverline2d.hh>
#include <dumux/material/fluidstates/nonequilibriumfluidstate.hh>
-#include <dumux/material/fluidstates/nonequilibriumenergyfluidstate.hh>
-#include <dumux/material/fluidstates/nonequilibriummassfluidstate.hh>
-#include <dumux/material/fluidstates/compositionalfluidstate.hh>
+//#include <dumux/material/fluidstates/nonequilibriumenergyfluidstate.hh>
+//#include <dumux/material/fluidstates/nonequilibriummassfluidstate.hh>
+//#include <dumux/material/fluidstates/compositionalfluidstate.hh>
#include "evaporationatmospherespatialparams.hh"
@@ -161,7 +161,7 @@ SET_BOOL_PROP(EvaporationAtmosphereProblem, NewtonWriteConvergence, false);
// Specify whether there is any energy equation
SET_BOOL_PROP(EvaporationAtmosphereProblem, EnableEnergy, true );
// Specify whether the kinetic energy module is used
-SET_BOOL_PROP(EvaporationAtmosphereProblem, EnableKineticEnergy, true);
+SET_INT_PROP(EvaporationAtmosphereProblem, NumEnergyEquations, 3);
// Specify whether the kinetic mass module is use
SET_BOOL_PROP(EvaporationAtmosphereProblem, EnableKinetic, true);
//#################
@@ -236,7 +236,7 @@ class EvaporationAtmosphereProblem
enum { nCompIdx = FluidSystem::N2Idx};
enum { enableKinetic = GET_PROP_VALUE(TypeTag, EnableKinetic)};
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy)};
- enum { enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy)};
+ enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations)};
enum { velocityOutput = GET_PROP_VALUE(TypeTag, VtkAddVelocities)};
enum { velocityAveragingInModel = GET_PROP_VALUE(TypeTag, VelocityAveragingInModel)};
@@ -539,7 +539,7 @@ public:
for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++)
fluidState.setPressure(phaseIdx, pnInitial_);
- if(enableKineticEnergy){
+ if(numEnergyEquations){
fluidState.setTemperature(nPhaseIdx, TInject_ );
fluidState.setTemperature(wPhaseIdx, TInitial_ ); // this value is a good one, TInject does not work
}
@@ -559,7 +559,7 @@ public:
nPhaseIdx);
fluidState.setDensity(nPhaseIdx, density);
- if(enableKineticEnergy){
+ if(numEnergyEquations){
for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++){
const Scalar h = FluidSystem::enthalpy(fluidState,
dummyCache,
@@ -583,7 +583,7 @@ public:
priVars[conti00EqIdx + nPhaseIdx * numComponents + nCompIdx]
= -molarFlux * fluidState.moleFraction(nPhaseIdx, nCompIdx);
// energy equations are specified mass specifically
- if(enableKineticEnergy){
+ if(numEnergyEquations){
priVars[energyEq0Idx + nPhaseIdx] = - massFluxInjectedPhase
* fluidState.enthalpy(nPhaseIdx) ;
}
@@ -673,7 +673,7 @@ private:
for (int phaseIdx = 0; phaseIdx < numPhases ; ++phaseIdx) {
equilibriumFluidState.setSaturation(phaseIdx, S[phaseIdx]);
- if(enableKineticEnergy){
+ if(numEnergyEquations){
equilibriumFluidState.setTemperature(phaseIdx, TInitial_ );
}
else
@@ -727,7 +727,7 @@ private:
else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << pressureFormulation << " is invalid.");
// temperature
- if(enableEnergy or enableKineticEnergy)
+ if(enableEnergy or numEnergyEquations)
for (int energyEqIdx=0; energyEqIdx< numEnergyEqs; ++energyEqIdx)
priVars[energyEq0Idx + energyEqIdx] = T;
diff --git a/test/implicit/mpnc/evaporationatmospherespatialparams.hh b/test/implicit/mpnc/evaporationatmospherespatialparams.hh
index 5eec2849a8..2be1019741 100644
--- a/test/implicit/mpnc/evaporationatmospherespatialparams.hh
+++ b/test/implicit/mpnc/evaporationatmospherespatialparams.hh
@@ -151,7 +151,7 @@ class EvaporationAtmosphereSpatialParams : public ImplicitSpatialParams<TypeTag>
enum {wPhaseIdx = FluidSystem::wPhaseIdx};
enum {nPhaseIdx = FluidSystem::nPhaseIdx};
enum {sPhaseIdx = FluidSystem::sPhaseIdx};
- enum { enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy)};
+ enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations)};
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)};
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy)};
@@ -267,7 +267,7 @@ public:
fluidState.setSaturation(phaseIdx, S[phaseIdx]);
if(enableEnergy){
- if (enableKineticEnergy)
+ if (numEnergyEquations>1)
fluidState.setTemperature(phaseIdx,TInitial);
else
fluidState.setTemperature(TInitial);
diff --git a/test/implicit/mpnc/test_boxmpnckinetic.cc b/test/implicit/mpnc/test_boxmpnckinetic.cc
index 9bf299945f..538e1e313e 100644
--- a/test/implicit/mpnc/test_boxmpnckinetic.cc
+++ b/test/implicit/mpnc/test_boxmpnckinetic.cc
@@ -336,7 +336,7 @@ int start(int argc, char **argv)
int main(int argc, char** argv)
{
-#if !HAVE_UG || !HAVE_ALUGRID
+#if !HAVE_UG && !HAVE_ALUGRID
std::cout<<"Evaporation Atmosphere not built, needs either UG or ALU for the log mesh." << std::endl;
return 77;
#else
--
GitLab