diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index e1b02f19381f1af1e6fc07883ad6379d755f0d71..490c0ca489d69d8e3168153e2fd6e81c8ffada30 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -87,7 +87,8 @@ template<class TypeTag, class MyTypeTag>
struct ReplaceCompEqIdx { using type = UndefinedProperty; }; //!< The component balance index that should be replaced by the total mass/mole balance
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
/////////////////////////////////////////////
// Properties used by finite volume schemes:
/////////////////////////////////////////////
@@ -140,10 +141,14 @@ struct SolutionDependentAdvection { using type = UndefinedProperty; };
template<class TypeTag, class MyTypeTag>
struct MolecularDiffusionType { using type = UndefinedProperty; }; //!< The type for the calculation of the molecular diffusion fluxes
template<class TypeTag, class MyTypeTag>
+struct DispersionFluxType { using type = UndefinedProperty; }; //!< The type for the calculation of the dispersive fluxes
+template<class TypeTag, class MyTypeTag>
struct SolutionDependentMolecularDiffusion { using type = UndefinedProperty; }; //!< specifies if the parameters for the diffusive fluxes depend on the solution
template<class TypeTag, class MyTypeTag>
struct HeatConductionType { using type = UndefinedProperty; }; //!< The type for the calculation of the heat conduction fluxes
template<class TypeTag, class MyTypeTag>
+struct DispersionTensorType { using type = UndefinedProperty; }; //!< The type for the calculation of the dispersion tensor
+template<class TypeTag, class MyTypeTag>
struct SolutionDependentHeatConduction { using type = UndefinedProperty; }; //!< specifies if the parameters for the heat conduction fluxes depend on the solution
template<class TypeTag, class MyTypeTag>
diff --git a/dumux/flux/box/dispersionflux.hh b/dumux/flux/box/dispersionflux.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d9bd77b7fcc4e708c4872b88935669c6bbb5be8e
--- /dev/null
+++ b/dumux/flux/box/dispersionflux.hh
@@ -0,0 +1,141 @@
+// -*- 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 3 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
+ * \ingroup BoxFlux
+ * \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
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/discretization/extrusion.hh>
+#include <dumux/flux/traits.hh>
+
+namespace Dumux {
+
+// forward declaration
+template<class TypeTag, class DiscretizationMethod, ReferenceSystemFormulation referenceSystem>
+class OnePDispersionFluxImplementation;
+
+/*!
+ * \ingroup BoxFlux
+ * \brief Specialization of a dispersion flux for the box method
+ */
+template <class TypeTag, ReferenceSystemFormulation referenceSystem>
+class OnePDispersionFluxImplementation<TypeTag, DiscretizationMethods::Box, referenceSystem>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using Extrusion = Extrusion_t<GridGeometry>;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
+ using FluxVarCache = GetPropType<TypeTag, Properties::FluxVariablesCache>;
+ using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
+ using FluxTraits = typename Dumux::FluxTraits<FluxVariables>;
+ using BalanceEqOpts = GetPropType<TypeTag, Properties::BalanceEqOpts>;
+ using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+
+ enum { dim = GridView::dimension} ;
+ enum { dimWorld = GridView::dimensionworld} ;
+ enum
+ {
+ numPhases = ModelTraits::numFluidPhases(),
+ numComponents = ModelTraits::numFluidComponents()
+ };
+
+ using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+ using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>;
+
+ static constexpr bool stationaryVelocityField = FluxTraits::hasStationaryVelocityField();
+
+public:
+
+ //return the reference system
+ static constexpr ReferenceSystemFormulation referenceSystemFormulation()
+ { return referenceSystem; }
+
+ /*!
+ * \brief Returns the dispersive fluxes of all components within
+ * a fluid phase across the given sub-control volume face.
+ * The computed fluxes are given in mole/s or kg/s, depending
+ * on the template parameter ReferenceSystemFormulation.
+ */
+ static ComponentFluxVector compositionalDispersionFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int phaseIdx,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ ComponentFluxVector componentFlux(0.0);
+
+ // evaluate gradX at integration point and interpolate density
+ const DimWorldMatrix dispersionTensor = VolumeVariables::DispersionTensorType::dispersionTensor(problem, scvf, fvGeometry,
+ elemVolVars, elemFluxVarsCache);
+
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+ const auto& shapeValues = fluxVarsCache.shapeValues();
+
+ // density interpolation
+ Scalar rhoMassOrMole(0.0);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto rho = massOrMolarDensity(elemVolVars[scv], referenceSystem, phaseIdx);
+ rhoMassOrMole += rho * shapeValues[scv.indexInElement()][0];
+ }
+
+ for (int compIdx = 0; compIdx < numComponents; compIdx++)
+ {
+ // the mole/mass fraction gradient
+ Dune::FieldVector<Scalar, dimWorld> gradX(0.0);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto x = massOrMoleFraction(elemVolVars[scv], referenceSystem, phaseIdx, compIdx);
+ gradX.axpy(x, fluxVarsCache.gradN(scv.indexInElement()));
+ }
+
+ // compute the dispersion flux
+ componentFlux[compIdx] = -1.0 * rhoMassOrMole * vtmv(scvf.unitOuterNormal(), dispersionTensor, gradX) * scvf.area();
+ if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem())
+ componentFlux[phaseIdx] -= componentFlux[compIdx];
+ }
+ return componentFlux;
+ }
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/flux/dispersionflux.hh b/dumux/flux/dispersionflux.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f1d09969eea082eb422ef43999d33c793d49d76c
--- /dev/null
+++ b/dumux/flux/dispersionflux.hh
@@ -0,0 +1,30 @@
+// -*- 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 3 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
+ * \ingroup Flux
+ * \brief Dispersion flux for different discretization schemes
+ */
+#ifndef DUMUX_FLUX_DISPERSION_FLUX_HH
+#define DUMUX_FLUX_DISPERSION_FLUX_HH
+
+#include <dumux/flux/dispersionflux_fwd.hh>
+#include <dumux/flux/box/dispersionflux.hh>
+
+#endif
diff --git a/dumux/flux/dispersionflux_fwd.hh b/dumux/flux/dispersionflux_fwd.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8112a15a0bc1d238ebb186c0777b40258785c823
--- /dev/null
+++ b/dumux/flux/dispersionflux_fwd.hh
@@ -0,0 +1,47 @@
+// -*- 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 3 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
+ * \ingroup Flux
+ * \brief Dispersion flux for different discretization schemes
+ */
+#ifndef DUMUX_FLUX_DISPERSION_FWD_HH
+#define DUMUX_FLUX_DISPERSION_FWD_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/method.hh>
+
+namespace Dumux {
+
+// declaration of primary template
+template<class TypeTag, class DiscretizationMethod, ReferenceSystemFormulation referenceSystem>
+class OnePDispersionFluxImplementation;
+
+/*!
+ * \ingroup Flux
+ * \brief Evaluates the dispersive flux
+ */
+template<class TypeTag, ReferenceSystemFormulation referenceSystem = ReferenceSystemFormulation::massAveraged>
+using OnePDiffusiveDispersionFlux = OnePDispersionFluxImplementation<TypeTag,
+ typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod,
+ referenceSystem>;
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/material/fluidmatrixinteractions/CMakeLists.txt b/dumux/material/fluidmatrixinteractions/CMakeLists.txt
index ddc38f14ff65ea3a4fca6531c6f9f4bcedd9382c..5c00e429803c3670bc0af7ff878f28d98f46b966 100644
--- a/dumux/material/fluidmatrixinteractions/CMakeLists.txt
+++ b/dumux/material/fluidmatrixinteractions/CMakeLists.txt
@@ -3,6 +3,7 @@ add_subdirectory(1pia)
add_subdirectory(2p)
add_subdirectory(3p)
add_subdirectory(frictionlaws)
+add_subdirectory(dispersiontensors)
add_subdirectory(mp)
add_subdirectory(porenetwork)
diff --git a/dumux/material/fluidmatrixinteractions/dispersiontensors/CMakeLists.txt b/dumux/material/fluidmatrixinteractions/dispersiontensors/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e04467da7e102539ca2847f67f75893f76d2d864
--- /dev/null
+++ b/dumux/material/fluidmatrixinteractions/dispersiontensors/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_MATERIAL_FLUIDMATRIXINTERACTIONS_DISPERSIONTENSORS_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_MATERIAL_FLUIDMATRIXINTERACTIONS_DISPERSIONTENSORS_HEADERS}
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/material/fluidmatrixinteractions/dispersiontensors)
diff --git a/dumux/material/fluidmatrixinteractions/dispersiontensors/fulltensor.hh b/dumux/material/fluidmatrixinteractions/dispersiontensors/fulltensor.hh
new file mode 100644
index 0000000000000000000000000000000000000000..851133c6434f2c81417ffd8a702d9ce2ffe86dce
--- /dev/null
+++ b/dumux/material/fluidmatrixinteractions/dispersiontensors/fulltensor.hh
@@ -0,0 +1,57 @@
+// -*- 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 3 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
+ * \ingroup Fluidmatrixinteractions
+ * \copydoc Dumux::FullDispersionTensor
+ */
+#ifndef DUMUX_MATERIAL_FLUIDMATRIX_DISPERSIONTENSORS_FULLTENSOR_HH
+#define DUMUX_MATERIAL_FLUIDMATRIX_DISPERSIONTENSORS_FULLTENSOR_HH
+
+namespace Dumux {
+
+template<class TypeTag>
+class FullDispersionTensor
+{
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ static const int dimWorld = GridView::dimensionworld;
+ using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+
+public:
+
+ template <class ElementFluxVariablesCache>
+ static DimWorldMatrix dispersionTensor([[maybe_unused]] const Problem& problem,
+ [[maybe_unused]] const SubControlVolumeFace& scvf,
+ [[maybe_unused]] const FVElementGeometry& fvGeometry,
+ [[maybe_unused]] const ElementVolumeVariables& elemVolVars,
+ [[maybe_unused]] const ElementFluxVariablesCache& elemFluxVarsCache)
+ { return problem.spatialParams().dispersionTensor(scvf.center()); }
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh b/dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh
new file mode 100644
index 0000000000000000000000000000000000000000..18a2b00f6fb3d81534b3bd3e0d99db0c582ff16c
--- /dev/null
+++ b/dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh
@@ -0,0 +1,150 @@
+// -*- 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 3 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
+ * \ingroup Fluidmatrixinteractions
+ * \copydoc Dumux::ScheideggersDispersionTensor
+ */
+#ifndef DUMUX_MATERIAL_FLUIDMATRIX_DISPERSIONTENSORS_SCHEIDEGGER_HH
+#define DUMUX_MATERIAL_FLUIDMATRIX_DISPERSIONTENSORS_SCHEIDEGGER_HH
+
+#include <algorithm>
+#include <cmath>
+#include <dune/common/math.hh>
+
+namespace Dumux {
+
+namespace Detail {
+template <class Problem, class SubControlVolumeFace>
+using HasVelocityInSpatialParams = decltype(std::declval<Problem>().spatialParams().velocity(std::declval<SubControlVolumeFace>()));
+
+template<class Problem, class SubControlVolumeFace>
+static constexpr bool hasVelocityInSpatialParams()
+{ return Dune::Std::is_detected<HasVelocityInSpatialParams, Problem, SubControlVolumeFace>::value; }
+}
+
+template<class TypeTag>
+class ScheideggersDispersionTensor
+{
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ static const int dimWorld = GridView::dimensionworld;
+ using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+
+ using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
+ using FluxTraits = typename Dumux::FluxTraits<FluxVariables>;
+ static constexpr bool stationaryVelocityField = FluxTraits::hasStationaryVelocityField();
+
+public:
+ template <class ElementFluxVariablesCache>
+ static DimWorldMatrix dispersionTensor(const Problem& problem,
+ const SubControlVolumeFace& scvf,
+ [[maybe_unused]] const FVElementGeometry& fvGeometry,
+ [[maybe_unused]] const ElementVolumeVariables& elemVolVars,
+ [[maybe_unused]] const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ DimWorldMatrix dispersionTensor(0.0);
+
+ // Calculate Darcy's velocity
+ Dune::FieldVector<Scalar, dimWorld> velocity(0.0);
+ if constexpr (stationaryVelocityField)
+ {
+ if constexpr (!Detail::hasVelocityInSpatialParams<Problem,SubControlVolumeFace>() )
+ DUNE_THROW(Dune::NotImplemented, "\n Please provide the stationary velocity field in the spatialparams via a velocity function.");
+ else
+ velocity = problem.spatialParams().velocity(scvf);
+ }
+ else
+ {
+ if constexpr (FVElementGeometry::GridGeometry::discMethod == DiscretizationMethods::box)
+ {
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+ const auto& shapeValues = fluxVarsCache.shapeValues();
+
+ // get inside and outside permeability tensors and calculate the harmonic mean
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+ const auto K = problem.spatialParams().harmonicMean(insideVolVars.permeability(),
+ outsideVolVars.permeability(),
+ scvf.unitOuterNormal());
+
+ // evaluate gradP - rho*g at integration point
+ Dune::FieldVector<Scalar, dimWorld> gradP(0.0);
+ Scalar rho(0.0);
+ static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+
+ if (enableGravity)
+ rho += volVars.density(0)*shapeValues[scv.indexInElement()][0];
+ // the global shape function gradient
+ gradP.axpy(volVars.pressure(0), fluxVarsCache.gradN(scv.indexInElement()));
+ }
+
+ if (enableGravity)
+ gradP.axpy(-rho, problem.spatialParams().gravity(scvf.center()));
+
+ velocity = gradP;
+ velocity *= K;
+
+ velocity /= -0.5 * (insideVolVars.viscosity() + outsideVolVars.viscosity());
+ }
+ else
+ DUNE_THROW(Dune::NotImplemented, "\n Scheidegger Dispersion for compositional models is only implemented using the Box method.");
+ }
+
+ //calculate dispersion tensor
+
+ // collect the dispersion alphas at this location
+ std::array<Scalar,2> dispersivity = problem.spatialParams().dispersionAlphas(scvf.center());
+
+ //matrix multiplication of the velocity at the interface: vv^T
+ for (int i=0; i < dimWorld; i++)
+ for (int j = 0; j < dimWorld; j++)
+ dispersionTensor[i][j] = velocity[i]*velocity[j];
+
+ //normalize velocity product --> vv^T/||v||, [m/s]
+ Scalar vNorm = velocity.two_norm();
+
+ dispersionTensor /= vNorm;
+ if (vNorm < 1e-20)
+ dispersionTensor = 0;
+
+ //multiply with dispersivity difference: vv^T/||v||*(alphaL - alphaT), [m^2/s] --> alphaL = longitudinal disp., alphaT = transverse disp.
+ dispersionTensor *= (dispersivity[0] - dispersivity[1]);
+
+ //add ||v||*alphaT to the main diagonal:vv^T/||v||*(alphaL - alphaT) + ||v||*alphaT, [m^2/s]
+ for (int i = 0; i < dimWorld; i++)
+ dispersionTensor[i][i] += vNorm*dispersivity[1];
+
+ return dispersionTensor;
+ }
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/1pnc/model.hh b/dumux/porousmediumflow/1pnc/model.hh
index d35809893fe5199a57367fdeec9d21e9538c45c5..18024763d0a89a365f32536f610fd4f9727b6838 100644
--- a/dumux/porousmediumflow/1pnc/model.hh
+++ b/dumux/porousmediumflow/1pnc/model.hh
@@ -339,20 +339,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 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 DTL, class DT, class EDM, class ETCM>
struct NCNITraits : public BaseTraits
{
+ using DispersionTensorLaw = DTL;
using DiffusionType = DT;
using EffectiveDiffusivityModel = EDM;
using EffectiveThermalConductivityModel = ETCM;
};
- using EquilibriumVolVars = OnePNCVolumeVariables<NCNITraits<BaseTraits, DT, EDM, ETCM>>;
+ using EquilibriumVolVars = OnePNCVolumeVariables<NCNITraits<BaseTraits, DTL, DT, EDM, ETCM>>;
public:
- using type = NonEquilibriumVolumeVariables<NCNITraits<BaseTraits, DT, EDM, ETCM>, EquilibriumVolVars>;
+ using type = NonEquilibriumVolumeVariables<NCNITraits<BaseTraits, DTL, DT, EDM, ETCM>, EquilibriumVolVars>;
};
} // end namespace Properties
diff --git a/dumux/porousmediumflow/1pncmin/model.hh b/dumux/porousmediumflow/1pncmin/model.hh
index d0bf189ceddc30077687871a6fc0811d5cab7825..dbeb9e94104637ff33dcbd3b9c3a4b3fb01cfda0 100644
--- a/dumux/porousmediumflow/1pncmin/model.hh
+++ b/dumux/porousmediumflow/1pncmin/model.hh
@@ -115,18 +115,20 @@ private:
static_assert(FST::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid state");
using BaseTraits = OnePVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
+ using DTL = GetPropType<TypeTag, Properties::DispersionTensorLaw>;
using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
- template<class BaseTraits, class DT, class EDM>
+ template<class BaseTraits, class DTL, class DT, class EDM>
struct NCTraits : public BaseTraits
{
+ using DispersionTensorLaw = DTL;
using DiffusionType = DT;
using EffectiveDiffusivityModel = EDM;
};
- using NonMinVolVars = OnePNCVolumeVariables<NCTraits<BaseTraits, DT, EDM>>;
+ using NonMinVolVars = OnePNCVolumeVariables<NCTraits<BaseTraits, DTL, DT, EDM>>;
public:
- using type = MineralizationVolumeVariables<NCTraits<BaseTraits, DT, EDM>, NonMinVolVars>;
+ using type = MineralizationVolumeVariables<NCTraits<BaseTraits, DTL, DT, EDM>, NonMinVolVars>;
};
// Use the mineralization local residual
@@ -201,19 +203,21 @@ private:
static_assert(FST::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid state");
using BaseTraits = OnePVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
+ using DTL = GetPropType<TypeTag, Properties::DispersionTensorLaw>;
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 DTL, class DT, class EDM, class ETCM>
struct NCNITraits : public BaseTraits
{
+ using DispersionTensorLaw = DTL;
using DiffusionType = DT;
using EffectiveDiffusivityModel = EDM;
using EffectiveThermalConductivityModel = ETCM;
};
- using NonMinVolVars = OnePNCVolumeVariables<NCNITraits<BaseTraits, DT, EDM, ETCM>>;
+ using NonMinVolVars = OnePNCVolumeVariables<NCNITraits<BaseTraits, DTL, DT, EDM, ETCM>>;
public:
- using type = MineralizationVolumeVariables<NCNITraits<BaseTraits, DT, EDM, ETCM>, NonMinVolVars>;
+ using type = MineralizationVolumeVariables<NCNITraits<BaseTraits, DTL, DT, EDM, ETCM>, NonMinVolVars>;
};
//! Use the average for effective conductivities
template<class TypeTag>
diff --git a/dumux/porousmediumflow/fluxvariables.hh b/dumux/porousmediumflow/fluxvariables.hh
index 3a3478f3138e39948cb1886702ab77332ffc8c2c..fea40c60b3406ed9a613773b032960032663466e 100644
--- a/dumux/porousmediumflow/fluxvariables.hh
+++ b/dumux/porousmediumflow/fluxvariables.hh
@@ -63,11 +63,13 @@ class PorousMediumFluxVariables
public:
using UpwindScheme = UpScheme;
using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
+ using DispersionFluxType = GetPropType<TypeTag, Properties::DispersionFluxType>;
using MolecularDiffusionType = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
using HeatConductionType = GetPropType<TypeTag, Properties::HeatConductionType>;
static constexpr bool enableAdvection = ModelTraits::enableAdvection();
static constexpr bool enableMolecularDiffusion = ModelTraits::enableMolecularDiffusion();
+ static constexpr bool enableCompositionalDispersion = ModelTraits::enableCompositionalDispersion();
static constexpr bool enableEnergyBalance = ModelTraits::enableEnergyBalance();
static constexpr bool enableThermalNonEquilibrium = getPropValue<TypeTag, Properties::EnableThermalNonEquilibrium>();
@@ -123,6 +125,23 @@ public:
return Dune::FieldVector<Scalar, numComponents>(0.0);
}
+ /*!
+ * \brief Returns the compositional dispersion flux computed by the respective law.
+ */
+ Dune::FieldVector<Scalar, numComponents> compositionalDispersionFlux([[maybe_unused]] const int phaseIdx) const
+ {
+ if constexpr (enableCompositionalDispersion)
+ return DispersionFluxType::compositionalDispersionFlux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->elemVolVars(),
+ this->scvFace(),
+ phaseIdx,
+ this->elemFluxVarsCache());
+ else
+ return Dune::FieldVector<Scalar, numComponents>(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/properties.hh b/dumux/porousmediumflow/properties.hh
index 3b5a58ba99c915f536a49bdb8df44362dd79771a..c69a546c6d5e649a8352bfca1abb5c991c78a651 100644
--- a/dumux/porousmediumflow/properties.hh
+++ b/dumux/porousmediumflow/properties.hh
@@ -38,10 +38,12 @@
#include <dumux/flux/darcyslaw.hh>
#include <dumux/flux/fickslaw.hh>
#include <dumux/flux/fourierslaw.hh>
+#include <dumux/flux/dispersionflux.hh>
#include <dumux/material/solidstates/inertsolidstate.hh>
#include <dumux/material/solidsystems/1csolid.hh>
#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh>
namespace Dumux {
namespace Properties {
@@ -73,6 +75,19 @@ struct AdvectionType<TypeTag, TTag::PorousMediumFlow> { using type = DarcysLaw<T
template<class TypeTag>
struct MolecularDiffusionType<TypeTag, TTag::PorousMediumFlow> { using type = FicksLaw<TypeTag>; };
+//! Per default, we do not include compositional dispersion effects
+template<class TypeTag>
+struct EnableCompositionalDispersion<TypeTag, TTag::PorousMediumFlow> { static constexpr bool value = false; };
+
+
+//! By default, we use a diffusive flux for the dispersive fluxes
+template<class TypeTag>
+struct DispersionFluxType<TypeTag, TTag::PorousMediumFlow> { using type = OnePDiffusiveDispersionFlux<TypeTag>; };
+
+//! By default, we use Scheideggers's law for the dispersive tensor calculation
+template<class TypeTag>
+struct DispersionTensorType<TypeTag, TTag::PorousMediumFlow> { using type = ScheideggersDispersionTensor<TypeTag>; };
+
//! By default, we use fourier's law as the default for heat conduction fluxes
template<class TypeTag>
struct HeatConductionType<TypeTag, TTag::PorousMediumFlow> { using type = FouriersLaw<TypeTag>; };