Commit 8a1d49e4 authored by Philipp Nuske's avatar Philipp Nuske
Browse files

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
parent 99d49adc
......@@ -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:
......
......@@ -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;
......
......@@ -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];
};
......
......@@ -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:
/*!
......
......@@ -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:
/*!
......
......@@ -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
}
/*!
......
......@@ -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
......
......@@ -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