Commit 4575b4d9 authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Ned Coltman
Browse files

[sagd][zerohyst] Adapt to new material law

parent 68f6c5ed
......@@ -27,9 +27,9 @@
#include <algorithm>
#include "parkervangenuchtenzerohysteresisparams.hh"
#include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh>
namespace Dumux {
namespace Dumux::FluidMatrix {
/*!
* \brief Implementation of van Genuchten's capillary pressure <->
* saturation relation. This class bundles the "raw" curves
......@@ -38,117 +38,290 @@ namespace Dumux {
*
* \sa VanGenuchten, VanGenuchtenThreephase
*/
template <class ScalarT, class ParamsT = ParkerVanGenZeroHyst3PParams<ScalarT> >
class ParkerVanGenZeroHyst3P
template <class ScalarT>
class ParkerVanGenZeroHyst3P : Adapter<ParkerVanGenZeroHyst3P<ScalarT>, ThreePhasePcKrSw>
{
public:
using Params = ParamsT;
using Scalar = typename Params::Scalar;
using Scalar = ScalarT;
/*!
* \brief The capillary pressure-saturation curve.
*
* \brief The parameter type
* \tparam Scalar The scalar type
*/
static Scalar pc(const Params &params, Scalar sw)
struct Params
{
DUNE_THROW(Dune::NotImplemented, "Capillary pressures for three phases is not so simple! Use pcgn, pcnw, and pcgw");
}
Params()
{
betaGw_ = betaNw_ = betaGn_ = 1.;
}
DUNE_DEPRECATED_MSG("use pc() (uncapitalized 'c') instead")
static Scalar pC(const Params &params, Scalar swe)
{
return pc(params, swe);
}
Params(Scalar vgAlpha,
Scalar vgn,
Dune::FieldVector<Scalar, 4> residualSaturation,
Scalar betaNw = 1.,
Scalar betaGn = 1.,
Scalar betaGw = 1.,
bool regardSnr = false,
Scalar TrappedSatN = 1.)
{
setVgAlpha(vgAlpha);
setVgn(vgn);
setSwr(residualSaturation[0]);
setSnr(residualSaturation[1]);
setSgr(residualSaturation[2]);
setSwrx(residualSaturation[3]);
setBetaNw(betaNw);
setBetaGn(betaGn);
setBetaGw(betaGw);
setTrappedSatN(TrappedSatN);
};
Scalar TrappedSatN() const
{
return TrappedSatN_;
}
static Scalar pcgw(const Params &params, Scalar sw)
{
return 0;
}
void setTrappedSatN(Scalar v)
{
TrappedSatN_ = v;
}
DUNE_DEPRECATED_MSG("use pcgw() (uncapitalized 'cgw') instead")
static Scalar pCGW(const Params &params, Scalar sw)
{
return pcgw(params, sw);
}
/*!
* \brief Return the \f$\alpha\f$ shape parameter of van Genuchten's
* curve.
*/
Scalar vgAlpha() const
{
return vgAlpha_;
}
static Scalar pcnw(const Params &params, Scalar sw)
{
return 0;
}
/*!
* \brief Set the \f$\alpha\f$ shape parameter of van Genuchten's
* curve.
*/
void setVgAlpha(Scalar v)
{
vgAlpha_ = v;
}
DUNE_DEPRECATED_MSG("use pcnw() (uncapitalized 'cnw') instead")
static Scalar pCNW(const Params &params, Scalar sw)
{
return pcnw(params, sw);
}
/*!
* \brief Return the \f$m\f$ shape parameter of van Genuchten's
* curve.
*/
Scalar vgm() const
{
return vgm_;
}
static Scalar pcgn(const Params &params, Scalar St)
{
return 0;
}
/*!
* \brief Set the \f$m\f$ shape parameter of van Genuchten's
* curve.
*
* The \f$n\f$ shape parameter is set to \f$n = \frac{1}{1 - m}\f$
*/
void setVgm(Scalar m)
{
vgm_ = m; vgn_ = 1/(1 - vgm_);
}
DUNE_DEPRECATED_MSG("use pcgn() (uncapitalized 'cgn') instead")
static Scalar pCGN(const Params &params, Scalar St)
{
return pcgn(params, St);
}
/*!
* \brief Return the \f$n\f$ shape parameter of van Genuchten's
* curve.
*/
Scalar vgn() const
{
return vgn_;
}
static Scalar pcAlpha(const Params &params, Scalar sn)
{
return(1);
}
/*!
* \brief Set the \f$n\f$ shape parameter of van Genuchten's
* curve.
*
* The \f$n\f$ shape parameter is set to \f$m = 1 - \frac{1}{n}\f$
*/
void setVgn(Scalar n)
{
vgn_ = n; vgm_ = 1 - 1/vgn_;
}
DUNE_DEPRECATED_MSG("use pcAlpha() (uncapitalized 'c') instead")
static Scalar pCAlpha(const Params &params, Scalar sn)
{
return pcAlpha(params, sn);
}
/*!
* \brief Return the residual saturation.
*/
Scalar satResidual(int phaseIdx) const
{
switch (phaseIdx)
{
case 0:
return swr_;
break;
case 1:
return snr_;
break;
case 2:
return sgr_;
break;
}
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
/*!
* \brief The saturation-capillary pressure curve.
*
*/
static Scalar sw(const Params &params, Scalar pc)
{
DUNE_THROW(Dune::NotImplemented, "sw(pc) for three phases not implemented! Do it yourself!");
}
/*!
* \brief Set all residual saturations.
*/
void setResiduals(Dune::FieldVector<Scalar, 3> residualSaturation)
{
setSwr(residualSaturation[0]);
setSnr(residualSaturation[1]);
setSgr(residualSaturation[2]);
}
DUNE_DEPRECATED_MSG("use sw() (uncapitalized 's') instead")
static Scalar Sw(const Params &params, Scalar pc)
{
return sw(params, pc);
}
/*!
* \brief Returns the partial derivative of the capillary
* pressure to the effective saturation.
*
*/
static Scalar dpc_dsw(const Params &params, Scalar sw)
{
DUNE_THROW(Dune::NotImplemented, "dpc/dsw for three phases not implemented! Do it yourself!");
}
/*!
* \brief Return the residual wetting saturation.
*/
Scalar swr() const
{
return swr_;
}
DUNE_DEPRECATED_MSG("use dpc_dsw() (uncapitalized 'c', 's') instead")
static Scalar dpC_dSw(const Params &params, Scalar sw)
{
return dpc_dsw(params, sw);
}
/*!
* \brief Set the residual wetting saturation.
*/
void setSwr(Scalar input)
{
swr_ = input;
}
/*!
* \brief Returns the partial derivative of the effective
* saturation to the capillary pressure.
*/
static Scalar dsw_dpc(const Params &params, Scalar pc)
{
DUNE_THROW(Dune::NotImplemented, "dsw/dpc for three phases not implemented! Do it yourself!");
}
/*!
* \brief Return the residual nonwetting saturation.
*/
Scalar snr() const
{
return snr_;
}
DUNE_DEPRECATED_MSG("use dsw_dpc() (uncapitalized 's', 'c') instead")
static Scalar dSw_dpC(const Params &params, Scalar pc)
{
return dsw_dpc(params, pc);
}
/*!
* \brief Set the residual nonwetting saturation.
*/
void setSnr(Scalar input)
{
snr_ = input;
}
/*!
* \brief Return the residual gas saturation.
*/
Scalar sgr() const
{
return sgr_;
}
/*!
* \brief Set the residual gas saturation.
*/
void setSgr(Scalar input)
{
sgr_ = input;
}
Scalar swrx() const
{
return swrx_;
}
/*!
* \brief Set the residual gas saturation.
*/
void setSwrx(Scalar input)
{
swrx_ = input;
}
/*!
* \brief defines the scaling parameters of capillary pressure between the phases (=1 for Gas-Water)
*/
void setBetaNw(Scalar input)
{
betaNw_ = input;
}
void setBetaGn(Scalar input)
{
betaGn_ = input;
}
void setBetaGw(Scalar input)
{
betaGw_ = input;
}
/*!
* \brief Return the values for the beta scaling parameters of capillary pressure between the phases
*/
Scalar betaNw() const
{
return betaNw_;
}
Scalar betaGn() const
{
return betaGn_;
}
Scalar betaGw() const
{
return betaGw_;
}
bool krRegardsSnr() const
{
return krRegardsSnr_;
}
/*!
* \brief defines if residual n-phase saturation should be regarded in its relative permeability.
*/
void setKrRegardsSnr(bool input)
{
krRegardsSnr_ = input;
}
private:
Scalar vgAlpha_;
Scalar vgm_;
Scalar vgn_;
Scalar swr_;
Scalar snr_;
Scalar sgr_;
Scalar swrx_; /* (sw+sn)_r */
Scalar betaNw_;
Scalar betaGn_;
Scalar betaGw_;
Scalar TrappedSatN_;
bool krRegardsSnr_ ;
};
ParkerVanGenZeroHyst3P(const Params& params)
: params_(params)
{}
Scalar pcgw(const Scalar sw, const Scalar sn) const
{ return 0;}
Scalar pcnw(const Scalar sw, const Scalar sn) const
{ return 0; }
Scalar pcgn(const Scalar sw, const Scalar sn) const
{ return 0; }
Scalar pcAlpha(const Scalar sw, const Scalar sn) const
{ return 1; }
void updateTappedSn(const Scalar trappedSn)
{ params_.setTrappedSatN(trappedSn); }
/*!
* \brief The relative permeability for the wetting phase of
......@@ -159,15 +332,13 @@ public:
* (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models"
* MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.)
*
* \param sw saturation saturation of the water phase.
* \param sn saturation of the NAPL phase.
* \param sg saturation of the gas phase.
* \param saturation saturation of the water phase.
* \param params Array of parameters.
*/
static Scalar krw(const Params &params, Scalar saturation, Scalar sn, Scalar sg)
Scalar krw(const Scalar sw, const Scalar sn) const
{
//transformation to effective saturation
Scalar se = (saturation - params.swr()) / (1-params.TrappedSatN()-params.swr());
const Scalar se = (sw - params_.swr()) / (1 - params_.TrappedSatN() - params_.swr());
// regularization
if (se > 1.0)
......@@ -175,7 +346,7 @@ public:
if (se < 0.0)
return 0.;
Scalar r = 1. - std::pow(1 - std::pow(se, 1/params.vgm()), params.vgm());
const Scalar r = 1. - std::pow(1 - std::pow(se, 1/params_.vgm()), params_.vgm());
return std::sqrt(se)*r*r;
};
......@@ -192,15 +363,13 @@ public:
* Journal of Contaminant Hydrology 66 (2003), 261-285
*
*
* \param sw saturation of the water phase.
* \param sg saturation of the gas phase.
* \param saturation saturation of the NAPL phase.
* \param params Array of parameters.
* \param sw saturation saturation of the water phase.
* \param sn saturation of the NAPL phase.
*/
static Scalar krn(const Params &params, Scalar sw, Scalar saturation, Scalar sg)
Scalar krn(const Scalar sw, const Scalar sn) const
{
Scalar swe = std::min((sw - params.swr()) / (1 - params.TrappedSatN()- params.swr()), 1.);
Scalar ste = std::min((sw + saturation - params.swr()) / (1 - params.swr()), 1.);
Scalar swe = std::min((sw - params_.swr()) / (1 - params_.TrappedSatN()- params_.swr()), 1.);
Scalar ste = std::min((sw + sn - params_.swr()) / (1 - params_.swr()), 1.);
// regularization
if (swe <= 0.0)
......@@ -211,19 +380,19 @@ public:
return 0.;
Scalar krn_;
krn_ = std::pow(1 - std::pow(swe, 1/params.vgm()), params.vgm());
krn_ -= std::pow(1 - std::pow(ste, 1/params.vgm()), params.vgm());
krn_ = std::pow(1 - std::pow(swe, 1/params_.vgm()), params_.vgm());
krn_ -= std::pow(1 - std::pow(ste, 1/params_.vgm()), params_.vgm());
krn_ *= krn_;
if (params.krRegardsSnr())
if (params_.krRegardsSnr())
{
// regard Snr in the permeability of the n-phase, see Helmig1997
Scalar resIncluded = std::max(std::min((saturation - params.snr()/ (1-params.swr())), 1.), 0.);
const Scalar resIncluded = std::max(std::min((sn - params_.snr()/ (1-params_.swr())), 1.), 0.);
krn_ *= std::sqrt(resIncluded );
}
else
krn_ *= std::sqrt(saturation / (1 - params.swr())); // Hint: (ste - swe) = sn / (1-Srw)
krn_ *= std::sqrt(sn / (1 - params_.swr())); // Hint: (ste - swe) = sn / (1-Srw)
return krn_;
};
......@@ -237,15 +406,13 @@ public:
* (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models"
* MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.)
*
* \param sw saturation of the water phase.
* \param sw saturation saturation of the water phase.
* \param sn saturation of the NAPL phase.
* \param saturation saturation of the gas phase.
* \param params Array of parameters.
*/
static Scalar krg(const Params &params, Scalar sw, Scalar sn, Scalar saturation)
Scalar krg(const Scalar sw, const Scalar sn) const
{
// se = (sw+sn - Sgr)/(1-Sgr)
Scalar se = std::min(((1-saturation) - params.sgr()) / (1 - params.sgr()), 1.);
const Scalar st = sw + sn;
const Scalar se = std::min(((1-st) - params_.sgr()) / (1 - params_.sgr()), 1.);
// regularization
if (se > 1.0)
......@@ -255,55 +422,44 @@ public:
Scalar scalFact = 1.;
if (saturation<=0.1)
if (st<=0.1)
{
scalFact = (saturation - params.sgr())/(0.1 - params.sgr());
scalFact = (st - params_.sgr())/(0.1 - params_.sgr());
if (scalFact < 0.)
scalFact = 0.;
}
Scalar result = scalFact * std::pow(1 - se, 1.0/3.) * std::pow(1 - std::pow(se, 1/params.vgm()), 2*params.vgm());
return result;
return scalFact * std::pow(1 - se, 1.0/3.) * std::pow(1 - std::pow(se, 1/params_.vgm()), 2*params_.vgm());
};
/*!
* \brief The relative permeability for a phase.
* \param sw saturation of the water phase.
* \param sg saturation of the gas phase.
* \param sw saturation saturation of the water phase.
* \param sn saturation of the NAPL phase.
* \param params Array of parameters.
* \param phase indicator, The saturation of all phases.
*/
static Scalar kr(const Params &params, const int phase, const Scalar sw, const Scalar sn, const Scalar sg)
Scalar kr(const int phaseIdx, const Scalar sw, const Scalar sn) const
{
switch (phase)
switch (phaseIdx)
{
case 0:
return krw(params, sw, sn, sg);
return krw(sw, sn);
break;
case 1:
return krn(params, sw, sn, sg);
return krn(sw, sn);
break;
case 2:
return krg(params, sw, sn, sg);
return krg(sw, sn);
break;
}
return 0;
};
/*
* \brief the basis for calculating adsorbed NAPL in storage term
* \param bulk density of porous medium, adsorption coefficient
*/
static Scalar bulkDensTimesAdsorpCoeff (const Params &params)
{
return params.rhoBulk() * params.KdNAPL();
}
private:
Params params_;
};
} // end namespace Dumux
} // end namespace Dumux::FluidMatrix
#endif
......@@ -30,8 +30,6 @@
#include <dumux/porousmediumflow/3pwateroil/indices.hh>
#include "../3p/parkervangenuchtenzerohysteresis.hh"
#include "../3p/parkervangenuchtenzerohysteresisparams.hh"
namespace Dumux {
......@@ -57,12 +55,10 @@ class SagdSpatialParams
static constexpr bool isBox = FVGridGeometry::discMethod == DiscretizationMethod::box;
static constexpr int dofCodim = isBox ? GridView::dimension : 0;
// using EffectiveLaw = RegularizedParkerVanGen3P<Scalar>;
using ThreePhasePcKrSw = FluidMatrix::ParkerVanGenZeroHyst3P<Scalar>;
public:
//using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
using MaterialLaw = ParkerVanGenZeroHyst3P<Scalar>;
using MaterialLawParams = typename MaterialLaw::Params;
using PermeabilityType = Scalar;
/*!
......@@ -83,22 +79,28 @@ public:
finePorosity_ = 0.10;
coarsePorosity_ = 0.1;
typename ThreePhasePcKrSw::Params baseParamsFine;
typename ThreePhasePcKrSw::Params baseParamsCoarse;
// residual saturations
fineMaterialParams_.setSwr(0.1);
fineMaterialParams_.setSnr(0.09); //Residual of NAPL if there is no water
fineMaterialParams_.setSgr(0.01);
coarseMaterialParams_.setSwr(0.1);
coarseMaterialParams_.setSnr(0.09);
coarseMaterialParams_.setSgr(0.01);
baseParamsFine.setSwr(0.1);
baseParamsFine.setSnr(0.09); //Residual of NAPL if there is no water
baseParams