diff --git a/dumux/implicit/2p2c/2p2cfluxvariables.hh b/dumux/implicit/2p2c/2p2cfluxvariables.hh
index 64ea218dba2eef65bdcadbd8bd41b56dd7ea6716..2239ae86332539cf05dbb6bac8c79b240e9d2981 100644
--- a/dumux/implicit/2p2c/2p2cfluxvariables.hh
+++ b/dumux/implicit/2p2c/2p2cfluxvariables.hh
@@ -51,6 +51,7 @@ class TwoPTwoCFluxVariables : public GET_PROP_TYPE(TypeTag, BaseFluxVariables)
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+ typedef typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel) EffectiveDiffusivityModel;
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
@@ -72,7 +73,7 @@ class TwoPTwoCFluxVariables : public GET_PROP_TYPE(TypeTag, BaseFluxVariables)
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
-public:
+ public:
/*!
* \brief The constructor
*
@@ -89,7 +90,7 @@ public:
const int faceIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
- : BaseFluxVariables(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+ : BaseFluxVariables(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
{
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
density_[phaseIdx] = Scalar(0);
@@ -100,7 +101,7 @@ public:
calculateValues_(problem, element, elemVolVars);
}
-protected:
+ protected:
void calculateValues_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
@@ -111,7 +112,7 @@ protected:
idx < this->face().numFap;
idx++) // loop over adjacent vertices
{
- // index for the element volume variables
+ // index for the element volume variables
int volVarsIdx = this->face().fapIndices[idx];
for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
@@ -140,7 +141,7 @@ protected:
// FE gradient at vertex idx
const DimVector &feGrad = this->face().grad[idx];
- // index for the element volume variables
+ // index for the element volume variables
int volVarsIdx = this->face().fapIndices[idx];
// the mole fraction gradient of the wetting phase
@@ -177,36 +178,37 @@ protected:
const VolumeVariables &volVarsI = elemVolVars[this->face().i];
const VolumeVariables &volVarsJ = elemVolVars[this->face().j];
+ // the effective diffusion coefficients at vertex i and j
+ Scalar diffCoeffI;
+ Scalar diffCoeffJ;
+
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
// make sure to only calculate diffusion coefficients
// for phases which exist in both finite volumes
- if (volVarsI.saturation(phaseIdx) <= 0 ||
- volVarsJ.saturation(phaseIdx) <= 0)
+ if (volVarsI.saturation(phaseIdx) <= 0 || volVarsJ.saturation(phaseIdx) <= 0)
{
porousDiffCoeff_[phaseIdx] = 0.0;
continue;
}
- // calculate tortuosity at the nodes i and j needed
- // for porous media diffusion coefficient
- Scalar tauI =
- 1.0/(volVarsI.porosity() * volVarsI.porosity()) *
- pow(volVarsI.porosity() * volVarsI.saturation(phaseIdx), 7.0/3);
- Scalar tauJ =
- 1.0/(volVarsJ.porosity() * volVarsJ.porosity()) *
- pow(volVarsJ.porosity() * volVarsJ.saturation(phaseIdx), 7.0/3);
- // Diffusion coefficient in the porous medium
+ diffCoeffI = EffectiveDiffusivityModel::effectiveDiffusivity(volVarsI.porosity(),
+ volVarsI.saturation(phaseIdx),
+ volVarsI.diffCoeff(phaseIdx));
+
+ diffCoeffJ = EffectiveDiffusivityModel::effectiveDiffusivity(volVarsJ.porosity(),
+ volVarsJ.saturation(phaseIdx),
+ volVarsJ.diffCoeff(phaseIdx));
// -> harmonic mean
- porousDiffCoeff_[phaseIdx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * volVarsI.diffCoeff(phaseIdx),
- volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * volVarsJ.diffCoeff(phaseIdx));
+ porousDiffCoeff_[phaseIdx] = harmonicMean(diffCoeffI, diffCoeffJ);
}
}
-public:
+ public:
/*!
- * \brief The binary diffusion coefficient for each fluid phase.
+ * \brief The effective diffusion coefficient \f$\mathrm{[m^2/s]}\f$
+ * for each fluid phase in the porous medium.
*/
Scalar porousDiffCoeff(int phaseIdx) const
{ return porousDiffCoeff_[phaseIdx]; };
@@ -229,7 +231,7 @@ public:
const DimVector &moleFractionGrad(int phaseIdx) const
{ return moleFractionGrad_[phaseIdx]; };
-protected:
+ protected:
// mole fraction gradients
DimVector moleFractionGrad_[numPhases];
diff --git a/dumux/implicit/2p2c/2p2cproperties.hh b/dumux/implicit/2p2c/2p2cproperties.hh
index 739639c2c7a55fb8ecfe70c27c2a45a5ffdc8544..1c87bd9ea8dbfac2ffe96f2bdb1938761b0d29f7 100644
--- a/dumux/implicit/2p2c/2p2cproperties.hh
+++ b/dumux/implicit/2p2c/2p2cproperties.hh
@@ -59,6 +59,7 @@ NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters)
NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extracted from the spatial parameters)
+NEW_PROP_TAG(EffectiveDiffusivityModel); //!< The employed model for the computation of the effective diffusivity
NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the upwind weight for the mass conservation equations
diff --git a/dumux/implicit/2p2c/2p2cpropertydefaults.hh b/dumux/implicit/2p2c/2p2cpropertydefaults.hh
index 6aa279c55aca97f5503cf6016fd80acc4ad913ba..15814127ab2e5a3194377b13a24d1d17034851b0 100644
--- a/dumux/implicit/2p2c/2p2cpropertydefaults.hh
+++ b/dumux/implicit/2p2c/2p2cpropertydefaults.hh
@@ -36,6 +36,7 @@
#include "2p2clocalresidual.hh"
#include "2p2cnewtoncontroller.hh"
+#include <dumux/material/fluidmatrixinteractions/2p/diffusivitymillingtonquirk.hh>
#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
#include <dumux/material/spatialparams/implicitspatialparams.hh>
@@ -132,17 +133,25 @@ SET_SCALAR_PROP(TwoPTwoC, ImplicitMassUpwindWeight, 1.0);
SET_SCALAR_PROP(TwoPTwoC, ImplicitMobilityUpwindWeight, 1.0);
//! The indices required by the isothermal 2p2c model
-SET_PROP(TwoPTwoC, Indices)
+SET_PROP(TwoPTwoC, Indices)
{ private:
enum { Formulation = GET_PROP_VALUE(TypeTag, Formulation) };
public:
typedef TwoPTwoCIndices<TypeTag, Formulation, 0> type;
};
-//! The spatial parameters to be employed.
+//! The spatial parameters to be employed.
//! Use ImplicitSpatialParams by default.
SET_TYPE_PROP(TwoPTwoC, SpatialParams, ImplicitSpatialParams<TypeTag>);
+//! The model after Millington (1961) is used for the effective diffusivity
+SET_PROP(TwoPTwoC, EffectiveDiffusivityModel)
+{ private :
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ public:
+ typedef DiffusivityMillingtonQuirk<Scalar> type;
+};
+
// disable velocity output by default
SET_BOOL_PROP(TwoPTwoC, VtkAddVelocity, false);