diff --git a/dumux/common/deprecated.hh b/dumux/common/deprecated.hh
index b3eed110cec0adf0be176ac91d3871c01f3a6d23..164a27832df9ecf14b2c976f75108666871913d8 100644
--- a/dumux/common/deprecated.hh
+++ b/dumux/common/deprecated.hh
@@ -75,6 +75,12 @@ template<class ModelTraits>
static constexpr bool hasEnableCompositionalDispersion()
{ return Dune::Std::is_detected<HasEnableCompositionalDispersionDetector, ModelTraits>::value; }
+template <typename ModelTraits>
+using HasEnableThermalDispersionDetector = decltype(ModelTraits::enableThermalDispersion());
+
+template<class ModelTraits>
+static constexpr bool hasEnableThermalDispersion()
+{ return Dune::Std::is_detected<HasEnableThermalDispersionDetector, ModelTraits>::value; }
} // end namespace Deprecated
#endif
diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index 490c0ca489d69d8e3168153e2fd6e81c8ffada30..95e1fdc9a62e946cd86a36c59d87e6f013b8da14 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -89,6 +89,9 @@ template<class TypeTag, class MyTypeTag>
struct BalanceEqOpts { using type = UndefinedProperty; }; //!< A class that collects options for the evaluation of the balance equations
template<class TypeTag, class MyTypeTag>
struct EnableCompositionalDispersion { using type = UndefinedProperty; }; //!< Property whether to include compositional dispersion
+template<class TypeTag, class MyTypeTag>
+struct EnableThermalDispersion { using type = UndefinedProperty; }; //!< Property whether to include thermal dispersion
+
/////////////////////////////////////////////
// Properties used by finite volume schemes:
/////////////////////////////////////////////
diff --git a/dumux/flux/box/dispersionflux.hh b/dumux/flux/box/dispersionflux.hh
index d9bd77b7fcc4e708c4872b88935669c6bbb5be8e..4607b3f0afb3b29eef98507faf98992e90873b08 100644
--- a/dumux/flux/box/dispersionflux.hh
+++ b/dumux/flux/box/dispersionflux.hh
@@ -22,8 +22,8 @@
* \brief This file contains the data which is required to calculate
* dispersive fluxes.
*/
-#ifndef DUMUX_DISCRETIZATION_BOX_DISPERSION_HH
-#define DUMUX_DISCRETIZATION_BOX_DISPERSION_LAW_HH
+#ifndef DUMUX_DISCRETIZATION_BOX_DISPERSION_FLUX_HH
+#define DUMUX_DISCRETIZATION_BOX_DISPERSION_FLUX_HH
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
@@ -76,6 +76,7 @@ class OnePDispersionFluxImplementation<TypeTag, DiscretizationMethods::Box, refe
using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>;
+ using HeatFluxScalar = Scalar;
static constexpr bool stationaryVelocityField = FluxTraits::hasStationaryVelocityField();
@@ -101,7 +102,7 @@ public:
{
ComponentFluxVector componentFlux(0.0);
- // evaluate gradX at integration point and interpolate density
+ // collect the dispersion tensor, the fluxVarsCache and the shape values
const DimWorldMatrix dispersionTensor = VolumeVariables::DispersionTensorType::dispersionTensor(problem, scvf, fvGeometry,
elemVolVars, elemFluxVarsCache);
@@ -134,6 +135,31 @@ public:
return componentFlux;
}
+ /*!
+ * \brief Returns the thermal dispersive flux
+ * across the given sub-control volume face.
+ */
+ static HeatFluxScalar thermalDispersionFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int phaseIdx,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ // collect the dispersion tensor
+ const DimWorldMatrix dispersionTensor = VolumeVariables::DispersionTensorType::dispersionTensor(problem, scvf, fvGeometry,
+ elemVolVars, elemFluxVarsCache);
+ // compute the temperature gradient with the shape functions
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+ Dune::FieldVector<Scalar, GridView::dimensionworld> gradTemp(0.0);
+ for (auto&& scv : scvs(fvGeometry))
+ gradTemp.axpy(elemVolVars[scv].temperature(), fluxVarsCache.gradN(scv.indexInElement()));
+
+ // compute the heat conduction flux
+ return -1.0*vtmv(scvf.unitOuterNormal(), dispersionTensor, gradTemp)*Extrusion::area(scvf);
+ }
+
};
} // end namespace Dumux
diff --git a/dumux/porousmediumflow/1pnc/model.hh b/dumux/porousmediumflow/1pnc/model.hh
index 986b99ad67e7b0b920839be644e39755a141d03a..239424c87a8f58a4bdf29f5b7ed9f5e7aad3109c 100644
--- a/dumux/porousmediumflow/1pnc/model.hh
+++ b/dumux/porousmediumflow/1pnc/model.hh
@@ -85,7 +85,7 @@ namespace Dumux {
*
* \tparam nComp the number of components to be considered.
*/
-template<int nComp, bool useM, int enableCompDisp, int repCompEqIdx = nComp>
+template<int nComp, bool useM, int enableCompDisp, int enableThermDisp, int repCompEqIdx = nComp>
struct OnePNCModelTraits
{
using Indices = OnePNCIndices;
@@ -99,6 +99,7 @@ struct OnePNCModelTraits
static constexpr bool enableAdvection() { return true; }
static constexpr bool enableMolecularDiffusion() { return true; }
static constexpr bool enableCompositionalDispersion() { return enableCompDisp; }
+ static constexpr bool enableThermalDispersion() { return enableThermDisp; }
static constexpr bool enableEnergyBalance() { return false; }
};
@@ -134,6 +135,7 @@ public:
using type = OnePNCModelTraits<FluidSystem::numComponents,
getPropValue<TypeTag, Properties::UseMoles>(),
getPropValue<TypeTag, Properties::EnableCompositionalDispersion>(),
+ getPropValue<TypeTag, Properties::EnableThermalDispersion>(),
getPropValue<TypeTag, Properties::ReplaceCompEqIdx>()>;
};
template<class TypeTag>
@@ -251,16 +253,17 @@ private:
using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
using ETCM = GetPropType< TypeTag, Properties:: ThermalConductivityModel>;
- template<class BaseTraits, class DT, class EDM, class ETCM>
+ template<class BaseTraits, class DTT, class DT, class EDM, class ETCM>
struct NCNITraits : public BaseTraits
{
+ using DispersionTensorType = DTT;
using DiffusionType = DT;
using EffectiveDiffusivityModel = EDM;
using EffectiveThermalConductivityModel = ETCM;
};
public:
- using type = OnePNCVolumeVariables<NCNITraits<BaseTraits, DT, EDM, ETCM>>;
+ using type = OnePNCVolumeVariables<NCNITraits<BaseTraits, DTT, DT, EDM, ETCM>>;
};
} // end namespace Properties
@@ -321,6 +324,7 @@ private:
using EquilibriumTraits = OnePNCModelTraits<FluidSystem::numComponents,
getPropValue<TypeTag, Properties::UseMoles>(),
getPropValue<TypeTag, Properties::EnableCompositionalDispersion>(),
+ getPropValue<TypeTag, Properties::EnableThermalDispersion>(),
getPropValue<TypeTag, Properties::ReplaceCompEqIdx>()>;
public:
using type = OnePNCUnconstrainedModelTraits<EquilibriumTraits>;
@@ -350,22 +354,22 @@ private:
using PT = typename GetPropType<TypeTag, Properties::SpatialParams>::PermeabilityType;
using BaseTraits = OnePVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
- using DTL = GetPropType<TypeTag, Properties::DispersionTensorLaw>;
+ using DTT = GetPropType<TypeTag, Properties::DispersionTensorType>;
using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
using ETCM = GetPropType< TypeTag, Properties:: ThermalConductivityModel>;
- template<class BaseTraits, class DTL, class DT, class EDM, class ETCM>
+ template<class BaseTraits, class DTT, class DT, class EDM, class ETCM>
struct NCNITraits : public BaseTraits
{
- using DispersionTensorLaw = DTL;
+ using DispersionTensorType = DTT;
using DiffusionType = DT;
using EffectiveDiffusivityModel = EDM;
using EffectiveThermalConductivityModel = ETCM;
};
- using EquilibriumVolVars = OnePNCVolumeVariables<NCNITraits<BaseTraits, DTL, DT, EDM, ETCM>>;
+ using EquilibriumVolVars = OnePNCVolumeVariables<NCNITraits<BaseTraits, DTT, DT, EDM, ETCM>>;
public:
- using type = NonEquilibriumVolumeVariables<NCNITraits<BaseTraits, DTL, DT, EDM, ETCM>, EquilibriumVolVars>;
+ using type = NonEquilibriumVolumeVariables<NCNITraits<BaseTraits, DTT, DT, EDM, ETCM>, EquilibriumVolVars>;
};
} // end namespace Properties
diff --git a/dumux/porousmediumflow/compositional/localresidual.hh b/dumux/porousmediumflow/compositional/localresidual.hh
index d06028eeb6fe6565f36cf92e9a431a9ce6de6e0c..fefff8ed512c8ad5cc3cb2610481b323adc582d1 100644
--- a/dumux/porousmediumflow/compositional/localresidual.hh
+++ b/dumux/porousmediumflow/compositional/localresidual.hh
@@ -233,8 +233,9 @@ public:
enableCompositionalDispersionMissing_<ModelTraits>();
}
- //! Add diffusive energy fluxes. For isothermal model the contribution is zero.
+ //! Add diffusive and dispersive energy fluxes. For isothermal model the contribution is zero.
EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
+ EnergyLocalResidual::heatDispersionFlux(flux, fluxVars);
return flux;
}
diff --git a/dumux/porousmediumflow/fluxvariables.hh b/dumux/porousmediumflow/fluxvariables.hh
index fea40c60b3406ed9a613773b032960032663466e..d5e40cd33de018b581aa883ef8bdeadf7ca50289 100644
--- a/dumux/porousmediumflow/fluxvariables.hh
+++ b/dumux/porousmediumflow/fluxvariables.hh
@@ -70,6 +70,7 @@ public:
static constexpr bool enableAdvection = ModelTraits::enableAdvection();
static constexpr bool enableMolecularDiffusion = ModelTraits::enableMolecularDiffusion();
static constexpr bool enableCompositionalDispersion = ModelTraits::enableCompositionalDispersion();
+ static constexpr bool enableThermalDispersion = ModelTraits::enableThermalDispersion();
static constexpr bool enableEnergyBalance = ModelTraits::enableEnergyBalance();
static constexpr bool enableThermalNonEquilibrium = getPropValue<TypeTag, Properties::EnableThermalNonEquilibrium>();
@@ -142,6 +143,23 @@ public:
return Dune::FieldVector<Scalar, numComponents>(0.0);
}
+ /*!
+ * \brief Returns the thermal dispersion flux computed by the respective law.
+ */
+ Dune::FieldVector<Scalar, 1> thermalDispersionFlux([[maybe_unused]] const int phaseIdx = 0) const
+ {
+ if constexpr (enableThermalDispersion)
+ return DispersionFluxType::thermalDispersionFlux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->elemVolVars(),
+ this->scvFace(),
+ phaseIdx,
+ this->elemFluxVarsCache());
+ else
+ return Dune::FieldVector<Scalar, 1>(0.0);
+ }
+
/*!
* \brief Returns the conductive flux computed by the respective law.
* \note This overload is used in models considering local thermal equilibrium
diff --git a/dumux/porousmediumflow/nonequilibrium/model.hh b/dumux/porousmediumflow/nonequilibrium/model.hh
index 30301c92ce73e46b0f2bdeb2446d91d211b51670..9698e89e16c5027cb8dc05d384417d747c1b6e27 100644
--- a/dumux/porousmediumflow/nonequilibrium/model.hh
+++ b/dumux/porousmediumflow/nonequilibrium/model.hh
@@ -62,6 +62,8 @@ struct NonEquilibriumModelTraits : public ET
static constexpr int numEnergyEqSolid() { return therm ? numES : 0; }
static constexpr int numEnergyEq() { return numEnergyEqFluid()+numEnergyEqSolid(); }
+ static constexpr bool enableCompositionalDispersion() { return false; }
+ static constexpr bool enableThermalDispersion() { return false; }
static constexpr bool enableEnergyBalance() { return ET::enableEnergyBalance() || therm; }
static constexpr bool enableThermalNonEquilibrium() { return therm; }
static constexpr bool enableChemicalNonEquilibrium() { return chem; }
diff --git a/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh b/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
index 63180863a97cba6bc1eadaf01c686aaa0883ccd9..173b9ef28cc7f6412ad3bd8cb9ed3374fff2836c 100644
--- a/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
+++ b/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
@@ -100,6 +100,15 @@ public:
}
}
+ /*!
+ * \brief The dispersive energy fluxes
+ *
+ * \param flux The flux
+ * \param fluxVars The flux variables.
+ */
+ static void heatDispersionFlux(NumEqVector& flux,
+ FluxVariables& fluxVars)
+ {}
//! The advective phase energy fluxes
static void heatConvectionFlux(NumEqVector& flux,
@@ -299,6 +308,17 @@ public:
flux[energyEqSolidIdx+sPhaseIdx] += fluxVars.heatConductionFlux(numPhases + sPhaseIdx);
}
}
+
+ /*!
+ * \brief The dispersive energy fluxes
+ *
+ * \param flux The flux
+ * \param fluxVars The flux variables.
+ */
+ static void heatDispersionFlux(NumEqVector& flux,
+ FluxVariables& fluxVars)
+ {}
+
/*!
* \brief Calculates the source term of the equation.
*
diff --git a/dumux/porousmediumflow/nonisothermal/localresidual.hh b/dumux/porousmediumflow/nonisothermal/localresidual.hh
index e0490095824b3992e222d22d9950fed5d98e2201..b52e5d8684f84637867e6f66408c7b6d451fbd25 100644
--- a/dumux/porousmediumflow/nonisothermal/localresidual.hh
+++ b/dumux/porousmediumflow/nonisothermal/localresidual.hh
@@ -28,6 +28,7 @@
#include <dumux/common/properties.hh>
#include <dumux/common/numeqvector.hh>
+#include <dumux/common/deprecated.hh>
namespace Dumux {
@@ -100,6 +101,16 @@ public:
static void heatConductionFlux(NumEqVector& flux,
FluxVariables& fluxVars)
{}
+
+ /*!
+ * \brief The dispersive energy fluxes
+ *
+ * \param flux The flux
+ * \param fluxVars The flux variables.
+ */
+ static void heatDispersionFlux(NumEqVector& flux,
+ FluxVariables& fluxVars)
+ {}
};
/*!
@@ -118,8 +129,10 @@ class EnergyLocalResidualImplementation<TypeTag, true>
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ using Indices = typename ModelTraits::Indices;
+ static constexpr int numPhases = ModelTraits::numFluidPhases();
enum { energyEqIdx = Indices::energyEqIdx };
public:
@@ -189,6 +202,30 @@ public:
flux[energyEqIdx] += fluxVars.heatConductionFlux();
}
+ /*!
+ * \brief The dispersive energy fluxes
+ *
+ * \param flux The flux
+ * \param fluxVars The flux variables.
+ */
+ static void heatDispersionFlux(NumEqVector& flux,
+ FluxVariables& fluxVars)
+ {
+ if constexpr (Deprecated::hasEnableThermalDispersion<ModelTraits>())
+ {
+ if constexpr (ModelTraits::enableCompositionalDispersion())
+ {
+ if constexpr (FVElementGeometry::GridGeometry::discMethod == DiscretizationMethods::box && numPhases == 1)
+ flux[energyEqIdx] += fluxVars.thermalDispersionFlux();
+ else
+ DUNE_THROW(Dune::NotImplemented, "Dispersion Fluxes are only implemented for single phase flows using the Box method.");
+ }
+ }
+ else
+ enableThermalDispersionMissing_<ModelTraits>();
+ }
+
+
/*!
* \brief heat transfer between the phases for nonequilibrium models
*
@@ -204,6 +241,14 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv)
{}
+
+private:
+
+ template <class T = ModelTraits>
+ [[deprecated("All non-isothermal models must specifiy if thermal dispersion is enabled."
+ "Please add enableThermalDispersion to the ModelTraits in your model header.")]]
+ static void enableThermalDispersionMissing_() {}
+
};
} // end namespace Dumux
diff --git a/dumux/porousmediumflow/properties.hh b/dumux/porousmediumflow/properties.hh
index c69a546c6d5e649a8352bfca1abb5c991c78a651..a357494bbeb77fc33142fb16fb13f4fe5d158e0d 100644
--- a/dumux/porousmediumflow/properties.hh
+++ b/dumux/porousmediumflow/properties.hh
@@ -79,6 +79,9 @@ struct MolecularDiffusionType<TypeTag, TTag::PorousMediumFlow> { using type = Fi
template<class TypeTag>
struct EnableCompositionalDispersion<TypeTag, TTag::PorousMediumFlow> { static constexpr bool value = false; };
+//! Per default, we do not include thermal dispersion effects
+template<class TypeTag>
+struct EnableThermalDispersion<TypeTag, TTag::PorousMediumFlow> { static constexpr bool value = false; };
//! By default, we use a diffusive flux for the dispersive fluxes
template<class TypeTag>