diff --git a/dumux/implicit/richards/richardsindices.hh b/dumux/implicit/richards/richardsindices.hh
index 8077abcad66fc4c5bcd74c16f4f6cea6d51f6392..6685c8e8b3609b90dc146c25368b6234d1503c31 100644
--- a/dumux/implicit/richards/richardsindices.hh
+++ b/dumux/implicit/richards/richardsindices.hh
@@ -46,7 +46,8 @@ struct RichardsIndices
//! Primary variable index for the wetting phase pressure
static const int pwIdx = 0;
-
+ //! Primary variable index for the wetting phase pressure head (used for pressure head formulation)
+ static const int hIdx = 0;
//////////
// equation indices
//////////
diff --git a/dumux/implicit/richards/richardslocalresidual.hh b/dumux/implicit/richards/richardslocalresidual.hh
index 14918e5859ad4e22d2adedb61889f01032b0ad7e..641f76ca27941943d01563459cbd000aedeac504 100644
--- a/dumux/implicit/richards/richardslocalresidual.hh
+++ b/dumux/implicit/richards/richardslocalresidual.hh
@@ -48,6 +48,7 @@ class RichardsLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
};
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ static const bool useHead = GET_PROP_VALUE(TypeTag, UseHead);
public:
/*!
@@ -87,10 +88,13 @@ public:
this->curVolVars_(scvIdx);
// partial time derivative of the wetting phase mass
+ //pressure head formulation
storage[contiEqIdx] =
- volVars.density(wPhaseIdx)
- * volVars.saturation(wPhaseIdx)
- * volVars.porosity();
+ volVars.saturation(wPhaseIdx)
+ * volVars.porosity();
+ // pressure formulation
+ if (!useHead)
+ storage[contiEqIdx] *= volVars.density(wPhaseIdx);
}
@@ -134,12 +138,14 @@ public:
const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(wPhaseIdx));
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(wPhaseIdx));
- flux[contiEqIdx] +=
- fluxVars.volumeFlux(wPhaseIdx)
- *
- (( massUpwindWeight_)*up.density(wPhaseIdx)
- +
- (1 - massUpwindWeight_)*dn.density(wPhaseIdx));
+ //pressure head formulation
+ flux[contiEqIdx] =
+ fluxVars.volumeFlux(wPhaseIdx);
+ //pressure formulation
+ if(!useHead)
+ flux[contiEqIdx] *=(( massUpwindWeight_)*up.density(wPhaseIdx)
+ +
+ (1 - massUpwindWeight_)*dn.density(wPhaseIdx));
}
/*!
diff --git a/dumux/implicit/richards/richardsmodel.hh b/dumux/implicit/richards/richardsmodel.hh
index 2f4ed78b235cb6405d9be7d49530a00d941b5234..7b772502c1f8ad55c17f6b368b685d1b56674ce7 100644
--- a/dumux/implicit/richards/richardsmodel.hh
+++ b/dumux/implicit/richards/richardsmodel.hh
@@ -111,6 +111,7 @@ class RichardsModel : public GET_PROP_TYPE(TypeTag, BaseModel)
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
@@ -126,6 +127,8 @@ class RichardsModel : public GET_PROP_TYPE(TypeTag, BaseModel)
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
+ static const bool useHead = GET_PROP_VALUE(TypeTag, UseHead);
+
public:
/*!
* \brief All relevant primary and secondary of a given
@@ -155,6 +158,10 @@ public:
ScalarField *mobN = writer.allocateManagedBuffer(numDofs);
ScalarField *poro = writer.allocateManagedBuffer(numDofs);
ScalarField *Te = writer.allocateManagedBuffer(numDofs);
+ ScalarField *ph = writer.allocateManagedBuffer(numDofs);
+ ScalarField *wc = writer.allocateManagedBuffer(numDofs);
+ ScalarField *source = writer.allocateManagedBuffer(numDofs);
+
VectorField *velocity = writer.template allocateManagedBuffer<double, dimWorld>(numDofs);
ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
@@ -199,7 +206,12 @@ public:
#else
int dofIdxGlobal = this->dofMapper().map(*eIt, scvIdx, dofCodim);
#endif
-
+ PrimaryVariables sourcevalues;
+ this->problem_().solDependentSource(sourcevalues,
+ *eIt,
+ fvGeometry,
+ scvIdx,
+ elemVolVars);
(*pw)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(wPhaseIdx);
(*pn)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(nPhaseIdx);
(*pc)[dofIdxGlobal] = elemVolVars[scvIdx].capillaryPressure();
@@ -211,6 +223,9 @@ public:
(*mobN)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(nPhaseIdx);
(*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
(*Te)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
+ (*ph)[dofIdxGlobal] = elemVolVars[scvIdx].pressureHead(wPhaseIdx);
+ (*wc)[dofIdxGlobal] = elemVolVars[scvIdx].waterContent(wPhaseIdx);
+ (*source)[dofIdxGlobal] = sourcevalues;
}
// velocity output
@@ -229,6 +244,10 @@ public:
writer.attachDofData(*mobN, "mobN", isBox);
writer.attachDofData(*poro, "porosity", isBox);
writer.attachDofData(*Te, "temperature", isBox);
+ writer.attachDofData(*ph, "pressure head", isBox);
+ writer.attachDofData(*wc, "water content", isBox);
+ writer.attachDofData(*source, "source", isBox);
+
if (velocityOutput.enableOutput())
{
writer.attachDofData(*velocity, "velocity", isBox, dim);
diff --git a/dumux/implicit/richards/richardsnewtoncontroller.hh b/dumux/implicit/richards/richardsnewtoncontroller.hh
index 37aa2a2537fd3b9e79b4feefa804621800729620..bb8fd82455e07238859aba894bde4b328075a1a7 100644
--- a/dumux/implicit/richards/richardsnewtoncontroller.hh
+++ b/dumux/implicit/richards/richardsnewtoncontroller.hh
@@ -58,6 +58,7 @@ class RichardsNewtonController : public NewtonController<TypeTag>
enum { dim = GridView::dimension };
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
+ static const bool useHead = GET_PROP_VALUE(TypeTag, UseHead);
public:
/*!
@@ -84,7 +85,7 @@ public:
{
ParentType::newtonUpdate(uCurrentIter, uLastIter, deltaU);
- if (!GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, UseLineSearch))
+ if ( (!GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, UseLineSearch)) && (!useHead) )
{
// do not clamp anything after 5 iterations
if (this->numSteps_ > 4)
diff --git a/dumux/implicit/richards/richardsproblem.hh b/dumux/implicit/richards/richardsproblem.hh
index e2b98ae91fff43ffb1e83b2912a7eaf7cbdee578..6bd16a9dbeec5da30a70b92619fc8574fa247219 100644
--- a/dumux/implicit/richards/richardsproblem.hh
+++ b/dumux/implicit/richards/richardsproblem.hh
@@ -37,7 +37,7 @@ namespace Dumux
* For a description of the Richards model, see Dumux::RichardsModel
*/
template<class TypeTag>
-class RichardsBoxProblem : public ImplicitPorousMediaProblem<TypeTag>
+class RichardsProblem : public ImplicitPorousMediaProblem<TypeTag>
{
typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
@@ -61,7 +61,7 @@ public:
* \param timeManager The TimeManager which keeps track of time
* \param gridView The GridView used by the problem.
*/
- RichardsBoxProblem(TimeManager &timeManager, const GridView &gridView)
+ RichardsProblem(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
{}
diff --git a/dumux/implicit/richards/richardsproperties.hh b/dumux/implicit/richards/richardsproperties.hh
index c390b7af808997693c1a253dabe067ac43e9a580..b587344cd343cc65b4ffaf774a7d636ae8e4e36e 100644
--- a/dumux/implicit/richards/richardsproperties.hh
+++ b/dumux/implicit/richards/richardsproperties.hh
@@ -71,6 +71,7 @@ NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the weight of the upwi
NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< The value of the weight for the upwind mobility in the velocity calculation
NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
+NEW_PROP_TAG(UseHead); //!< Defines whether pressure [Pa] (false) or pressure head [cm] (true) is used
// \}
}
diff --git a/dumux/implicit/richards/richardspropertydefaults.hh b/dumux/implicit/richards/richardspropertydefaults.hh
index aff4cd2956a220ae5d77ec2665f0fb03a2585d77..8909ce8cedffbeb100d9ebddaad8623711ec12ed 100644
--- a/dumux/implicit/richards/richardspropertydefaults.hh
+++ b/dumux/implicit/richards/richardspropertydefaults.hh
@@ -57,6 +57,9 @@ SET_INT_PROP(Richards, NumEq, GET_PROP_VALUE(TypeTag, IsothermalNumEq));
//! Number of fluid phases considered
SET_INT_PROP(Richards, NumPhases, 2);
+//! Use pressure [Pa] by default
+SET_BOOL_PROP(Richards, UseHead, false);
+
//! The local residual operator
SET_TYPE_PROP(Richards,
LocalResidual,
@@ -71,9 +74,9 @@ SET_TYPE_PROP(Richards, VolumeVariables, typename GET_PROP_TYPE(TypeTag, Isother
//! The class for the quantities required for the flux calculation
SET_TYPE_PROP(Richards, FluxVariables, typename GET_PROP_TYPE(TypeTag, IsothermalFluxVariables));
-//! The class of the newton controller
-SET_TYPE_PROP(Richards, NewtonController, RichardsNewtonController<TypeTag>);
-
+//! The class of the newton controller
+SET_TYPE_PROP(Richards, NewtonController, RichardsNewtonController<TypeTag>);
+
//! The upwind weight for the mass conservation equations
SET_SCALAR_PROP(Richards, ImplicitMassUpwindWeight, 1.0);
diff --git a/dumux/implicit/richards/richardsvolumevariables.hh b/dumux/implicit/richards/richardsvolumevariables.hh
index 0668e6e379f97ad60d03ffe4a7eec61d9c8c8873..2319ee51b6516d9da370949f6547872906748286 100644
--- a/dumux/implicit/richards/richardsvolumevariables.hh
+++ b/dumux/implicit/richards/richardsvolumevariables.hh
@@ -55,14 +55,16 @@ class RichardsVolumeVariables : public ImplicitVolumeVariables<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
- enum {
- pwIdx = Indices::pwIdx,
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx
+
+ enum{hIdx = Indices::hIdx,
+ pwIdx = Indices::pwIdx,
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx
};
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
+ static const bool useHead = GET_PROP_VALUE(TypeTag, UseHead);
public:
//! The type returned by the fluidState() method
@@ -118,20 +120,35 @@ public:
Scalar t = Implementation::temperature_(priVars, problem, element,
fvGeometry, scvIdx);
fluidState.setTemperature(t);
-
- // pressures
- Scalar pnRef = problem.referencePressure(element, fvGeometry, scvIdx);
+
const MaterialLawParams &matParams =
- problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
- Scalar minPc = MaterialLaw::pc(matParams, 1.0);
- fluidState.setPressure(wPhaseIdx, priVars[pwIdx]);
- fluidState.setPressure(nPhaseIdx, std::max(pnRef, priVars[pwIdx] + minPc));
-
- // saturations
- Scalar sw = MaterialLaw::sw(matParams, fluidState.pressure(nPhaseIdx) - fluidState.pressure(wPhaseIdx));
- fluidState.setSaturation(wPhaseIdx, sw);
- fluidState.setSaturation(nPhaseIdx, 1 - sw);
-
+ problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
+
+ pnRef_ = problem.referencePressure(element, fvGeometry, scvIdx);
+ gravity_ = problem.gravity().two_norm();
+
+ // pressure head formulation
+ if (useHead){
+ Scalar pw = (0.01*priVars[hIdx])*1000.0 * gravity_ ;
+ fluidState.setPressure(wPhaseIdx, pw);
+ fluidState.setPressure(nPhaseIdx, 0.0);
+
+ // saturations
+ Scalar sw = MaterialLaw::sw(matParams,-fluidState.pressure(wPhaseIdx));
+ fluidState.setSaturation(wPhaseIdx, sw);
+ fluidState.setSaturation(nPhaseIdx, 1 - sw);
+ }
+ else{ //pressure formulation
+
+ Scalar minPc = MaterialLaw::pc(matParams, 1.0);
+ fluidState.setPressure(wPhaseIdx, priVars[pwIdx]);
+ fluidState.setPressure(nPhaseIdx, std::max(pnRef_, priVars[pwIdx] + minPc));
+
+ // saturations
+ Scalar sw = MaterialLaw::sw(matParams, fluidState.pressure(nPhaseIdx) - fluidState.pressure(wPhaseIdx));
+ fluidState.setSaturation(wPhaseIdx, sw);
+ fluidState.setSaturation(nPhaseIdx, 1 - sw);
+ }
// density and viscosity
typename FluidSystem::ParameterCache paramCache;
paramCache.updateAll(fluidState);
@@ -141,7 +158,7 @@ public:
fluidState.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState, paramCache, wPhaseIdx));
fluidState.setViscosity(nPhaseIdx, 1e-10);
- // compute and set the enthalpy
+ // compute and set the enthalpy
fluidState.setEnthalpy(wPhaseIdx, Implementation::enthalpy_(fluidState, paramCache, wPhaseIdx));
fluidState.setEnthalpy(nPhaseIdx, Implementation::enthalpy_(fluidState, paramCache, nPhaseIdx));
@@ -245,7 +262,45 @@ public:
* \f[ p_c = p_n - p_w \f]
*/
Scalar capillaryPressure() const
- { return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx); }
+ {
+ // pressure head formulation
+ if (useHead)
+ return -fluidState_.pressure(wPhaseIdx);
+ else // pressure formulation
+ return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx);
+ }
+
+ /*!
+ * \brief Returns the pressureHead \f$\mathrm{[cm]}\f$ of a given phase within
+ * the control volume.
+ *
+ * For the non-wetting phase (i.e. the gas phase), we assume
+ * infinite mobility, which implies that the non-wetting phase
+ * pressure is equal to the finite volume's reference pressure
+ * defined by the problem.
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar pressureHead(const int phaseIdx) const
+ {
+ // pressure head formulation
+ if (useHead)
+ return (100.) *(fluidState_.pressure(phaseIdx))/ fluidState_.density(phaseIdx)/ gravity_;
+ else // pressure formulation
+ return (100.) *(fluidState_.pressure(phaseIdx) - pnRef_)/ fluidState_.density(phaseIdx)/ gravity_;
+ }
+
+ /*!
+ * \brief Returns the water content
+ * fluid phase within the finite volume.
+ *
+ * The water content is defined as the fraction of
+ * the saturation devided by the porosity
+
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar waterContent (const int phaseIdx) const
+ { return fluidState_.saturation(phaseIdx)* porosity_; }
protected:
static Scalar temperature_(const PrimaryVariables &primaryVariables,
@@ -279,6 +334,8 @@ protected:
FluidState fluidState_;
Scalar relativePermeabilityWetting_;
Scalar porosity_;
+ static Scalar pnRef_;
+ static Scalar gravity_;
private:
Implementation &asImp_()
@@ -288,6 +345,11 @@ private:
{ return *static_cast<const Implementation*>(this); }
};
+template <class TypeTag>
+typename RichardsVolumeVariables<TypeTag>::Scalar RichardsVolumeVariables<TypeTag>::pnRef_;
+template <class TypeTag>
+typename RichardsVolumeVariables<TypeTag>::Scalar RichardsVolumeVariables<TypeTag>::gravity_;
+
}
#endif
diff --git a/test/implicit/richards/richardslensproblem.hh b/test/implicit/richards/richardslensproblem.hh
index e6e24b436d3b8c791a72d9caafecc00d647796ea..0b5c916b8fb0ebac7806094ba41e93365d8b0a76 100644
--- a/test/implicit/richards/richardslensproblem.hh
+++ b/test/implicit/richards/richardslensproblem.hh
@@ -63,8 +63,11 @@ public:
typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
};
-// Set the maximum number of newton iterations of a time step
-SET_INT_PROP(RichardsLensProblem, NewtonMaxSteps, 28);
+// Enable gravity
+SET_BOOL_PROP(RichardsLensProblem, ProblemEnableGravity, true);
+
+//! Use pressure [Pa] or pressure head [cm] formulation
+SET_BOOL_PROP(RichardsLensProblem, UseHead, false);
}
/*!
@@ -73,7 +76,7 @@ SET_INT_PROP(RichardsLensProblem, NewtonMaxSteps, 28);
*
* \brief A water infiltration problem with a low-permeability lens
* embedded into a high-permeability domain which uses the
- * Richards box model.
+ * Richards model.
*
* The domain is box shaped. Left and right boundaries are Dirichlet
* boundaries with fixed water pressure (fixed Saturation \f$S_w = 0\f$),
@@ -95,9 +98,9 @@ SET_INT_PROP(RichardsLensProblem, NewtonMaxSteps, 28);
* simulation time is 10,000,000 seconds (115.7 days)
*/
template <class TypeTag>
-class RichardsLensProblem : public RichardsBoxProblem<TypeTag>
+class RichardsLensProblem : public RichardsProblem<TypeTag>
{
- typedef RichardsBoxProblem<TypeTag> ParentType;
+ typedef RichardsProblem<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
@@ -106,19 +109,24 @@ class RichardsLensProblem : public RichardsBoxProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
// copy some indices for convenience
pwIdx = Indices::pwIdx,
+ hIdx = Indices::hIdx,
contiEqIdx = Indices::contiEqIdx,
// Grid and world dimension
dimWorld = GridView::dimensionworld
};
-
+
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<0>::Entity Element;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ static const bool useHead = GET_PROP_VALUE(TypeTag, UseHead);
public:
/*!
@@ -234,7 +242,11 @@ public:
const GlobalPosition &globalPos) const
{
// use initial values as boundary conditions
- initial_(values, globalPos);
+ if (onInlet_(globalPos))
+ // inflow of water
+ values[contiEqIdx] = -0.;
+ else
+ initial_(values, globalPos);
}
/*!
@@ -253,7 +265,10 @@ public:
values = 0.0;
if (onInlet_(globalPos)) {
// inflow of water
- values[contiEqIdx] = -0.04; // kg/(m*s)
+ if (useHead)
+ values[contiEqIdx] = -0.04/1000.; // kg/(m*s) / density
+ else
+ values[contiEqIdx] = -0.04; // kg/(m*s)
}
}
@@ -284,7 +299,12 @@ private:
Scalar pc =
MaterialLaw::pc(this->spatialParams().materialLawParams(globalPos),
sw);
- values[pwIdx] = pnRef_ - pc;
+ Scalar g = this->gravity().two_norm();
+
+ if (useHead)
+ values[hIdx] = (-pc)/1000./ g *(100.);
+ else
+ values[pwIdx] = pnRef_ - pc;
}
bool onLeftBoundary_(const GlobalPosition &globalPos) const