From 68f6c5ed9c6bffaee0ba40cc68ce8cdd6fc99cdc Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sat, 31 Oct 2020 00:39:29 +0100 Subject: [PATCH 01/15] [remediationscenario] Adapt to new material law --- .../remediationscenariosexercise.input | 26 ++++++ .../remediationscenariosspatialparams.hh | 92 ++++--------------- 2 files changed, 44 insertions(+), 74 deletions(-) diff --git a/lecture/mm/remediationscenarios/remediationscenariosexercise.input b/lecture/mm/remediationscenarios/remediationscenariosexercise.input index cce99be..b9071d9 100644 --- a/lecture/mm/remediationscenarios/remediationscenariosexercise.input +++ b/lecture/mm/remediationscenarios/remediationscenariosexercise.input @@ -45,6 +45,32 @@ AirFlux = -0.02 ContaminantFlux = -0.0000001 HeatFlux = -5555.0 +[SpatialParams.Fine] +Permeability = 6.28e-12 +Porosity = 0.42 +Swr = 0.12 +Snr = 0.07 +Sgr = 0.01 +ParkerVanGenuchtenAlpha = 0.0005 +ParkerVanGenuchtenN = 4.0 +ParkerVanGenuchtenRegardSnrForKrn = true +ParkerVanGenuchtenBetaNw = 1.83 +ParkerVanGenuchtenBetaGn = 2.2 +ParkerVanGenuchtenBetaGw = 1.0 + +[SpatialParams.Coarse] +Permeability = 9.14e-10 +Porosity = 0.42 +Swr = 0.12 +Snr = 0.07 +Sgr = 0.01 +ParkerVanGenuchtenAlpha = 0.0015 +ParkerVanGenuchtenN = 4.0 +ParkerVanGenuchtenRegardSnrForKrn = true +ParkerVanGenuchtenBetaNw = 1.83 +ParkerVanGenuchtenBetaGn = 2.2 +ParkerVanGenuchtenBetaGw = 1.0 + [Newton] MaxRelativeShift = 1e-4 TargetSteps = 4 diff --git a/lecture/mm/remediationscenarios/remediationscenariosspatialparams.hh b/lecture/mm/remediationscenarios/remediationscenariosspatialparams.hh index 793559e..1eaabaa 100644 --- a/lecture/mm/remediationscenarios/remediationscenariosspatialparams.hh +++ b/lecture/mm/remediationscenarios/remediationscenariosspatialparams.hh @@ -23,16 +23,11 @@ #ifndef DUMUX_KUEVETTE3P3CNI_SPATIAL_PARAMS_HH #define DUMUX_KUEVETTE3P3CNI_SPATIAL_PARAMS_HH - -//TODO Do we need this #include - #include #include #include -#include -#include -#include +#include namespace Dumux { @@ -52,11 +47,9 @@ class KuevetteSpatialParams using GlobalPosition = typename SubControlVolume::GlobalPosition; using ParentType = FVSpatialParams>; - using EffectiveLaw = RegularizedParkerVanGen3P; + using ThreePhasePcKrSw = FluidMatrix::ParkerVanGenuchten3PDefault; public: - using MaterialLaw = EffToAbsLaw; - using MaterialLawParams = typename MaterialLaw::Params; using PermeabilityType = Scalar; /*! @@ -66,46 +59,16 @@ public: */ KuevetteSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) + , pcKrSwCurveFine_("SpatialParams.Fine") + , pcKrSwCurveCoarse_("SpatialParams.Coarse") { // intrinsic permeabilities - fineK_ = 6.28e-12; - coarseK_ = 9.14e-10; + fineK_ = getParam("SpatialParams.Fine.Permeability"); + coarseK_ = getParam("SpatialParams.Coarse.Permeability"); // porosities - finePorosity_ = 0.42; - coarsePorosity_ = 0.42; - - // residual saturations - fineMaterialParams_.setSwr(0.12); - fineMaterialParams_.setSnr(0.07); - fineMaterialParams_.setSgr(0.01); - coarseMaterialParams_.setSwr(0.12); - coarseMaterialParams_.setSnr(0.07); - coarseMaterialParams_.setSgr(0.01); - - // parameters for the scaling of capillary pressure (GW = 1); - fineMaterialParams_.setBetaNw(1.83); - fineMaterialParams_.setBetaGn(2.2); - fineMaterialParams_.setBetaGw(1.0); - coarseMaterialParams_.setBetaNw(1.83); - coarseMaterialParams_.setBetaGn(2.2); - coarseMaterialParams_.setBetaGw(1.0); - - - // parameters for the 3phase van Genuchten law - fineMaterialParams_.setVgAlpha(0.0005); - coarseMaterialParams_.setVgAlpha(0.0015); - fineMaterialParams_.setVgn(4.0); - coarseMaterialParams_.setVgn(4.0); - - coarseMaterialParams_.setKrRegardsSnr(true); - fineMaterialParams_.setKrRegardsSnr(true); - - // parameters for adsorption - coarseMaterialParams_.setKdNAPL(0.); - coarseMaterialParams_.setRhoBulk(1500.); - fineMaterialParams_.setKdNAPL(0.); - fineMaterialParams_.setRhoBulk(1500.); + finePorosity_ = getParam("SpatialParams.Fine.Porosity"); + coarsePorosity_ = getParam("SpatialParams.Coarse.Porosity"); } /*! @@ -156,49 +119,30 @@ public: } } -/* - Scalar porosityAtPos(const GlobalPosition& globalPos) const - { - if (isFineMaterial_(globalPos)) - return finePorosity_; - else - return coarsePorosity_; - } -*/ - - /*! - * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.). + * \brief Returns the fluid-matrix interaction law (kr-sw, pc-sw, etc.). * * \param element The current element * \param scv The sub-control volume inside the element. * \param elemSol The solution at the dofs connected to the element. - * \return the material parameters object */ template - const MaterialLawParams& materialLawParams(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol) const + auto fluidMatrixInteraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const { const auto& globalPos = scv.dofPosition(); if (isFineMaterial_(globalPos)) - return fineMaterialParams_; + return makeFluidMatrixInteraction(pcKrSwCurveFine_); else - return coarseMaterialParams_; + return makeFluidMatrixInteraction(pcKrSwCurveCoarse_); } private: bool isFineMaterial_(const GlobalPosition &globalPos) const { - if (0.13 < globalPos[0] + eps_ && 1.20 > globalPos[0] - eps_ && - 0.32 < globalPos[1] + eps_ && globalPos[1] < 0.57 - eps_) - { - return true; - } - else - { - return false; - } + return (0.13 < globalPos[0] + eps_ && 1.20 > globalPos[0] - eps_ && + 0.32 < globalPos[1] + eps_ && globalPos[1] < 0.57 - eps_); } Scalar fineK_; @@ -207,8 +151,8 @@ private: Scalar finePorosity_; Scalar coarsePorosity_; - MaterialLawParams fineMaterialParams_; - MaterialLawParams coarseMaterialParams_; + const ThreePhasePcKrSw pcKrSwCurveFine_; + const ThreePhasePcKrSw pcKrSwCurveCoarse_; static constexpr Scalar eps_ = 1.5e-7; }; -- GitLab From 4575b4d949728cc3314bfe0151c7817ed73a80f9 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sat, 31 Oct 2020 00:43:03 +0100 Subject: [PATCH 02/15] [sagd][zerohyst] Adapt to new material law --- .../3p/parkervangenuchtenzerohysteresis.hh | 434 ++++++++++++------ .../heavyoil/sagdcyclichyst/spatialparams.hh | 70 +-- 2 files changed, 332 insertions(+), 172 deletions(-) diff --git a/lecture/mm/heavyoil/3p/parkervangenuchtenzerohysteresis.hh b/lecture/mm/heavyoil/3p/parkervangenuchtenzerohysteresis.hh index b7f54ac..aeb2c8c 100644 --- a/lecture/mm/heavyoil/3p/parkervangenuchtenzerohysteresis.hh +++ b/lecture/mm/heavyoil/3p/parkervangenuchtenzerohysteresis.hh @@ -27,9 +27,9 @@ #include -#include "parkervangenuchtenzerohysteresisparams.hh" +#include -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 ParkerVanGenZeroHyst3P +template +class ParkerVanGenZeroHyst3P : Adapter, 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 ¶ms, 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 ¶ms, Scalar swe) - { - return pc(params, swe); - } + Params(Scalar vgAlpha, + Scalar vgn, + Dune::FieldVector 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 ¶ms, Scalar sw) - { - return 0; - } + void setTrappedSatN(Scalar v) + { + TrappedSatN_ = v; + } - DUNE_DEPRECATED_MSG("use pcgw() (uncapitalized 'cgw') instead") - static Scalar pCGW(const Params ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 residualSaturation) + { + setSwr(residualSaturation[0]); + setSnr(residualSaturation[1]); + setSgr(residualSaturation[2]); + } - DUNE_DEPRECATED_MSG("use sw() (uncapitalized 's') instead") - static Scalar Sw(const Params ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms) - { - return params.rhoBulk() * params.KdNAPL(); - } - +private: + Params params_; }; -} // end namespace Dumux +} // end namespace Dumux::FluidMatrix #endif diff --git a/lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh b/lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh index ef7a842..acae9de 100644 --- a/lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh +++ b/lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh @@ -30,8 +30,6 @@ #include #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; + + using ThreePhasePcKrSw = FluidMatrix::ParkerVanGenZeroHyst3P; public: - //using MaterialLaw = EffToAbsLaw; - using MaterialLaw = ParkerVanGenZeroHyst3P; - 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 + baseParamsFine.setSgr(0.01); + baseParamsCoarse.setSwr(0.1); + baseParamsCoarse.setSnr(0.09); + baseParamsCoarse.setSgr(0.01); // parameters for the 3phase van Genuchten law - fineMaterialParams_.setVgn(4.0); - coarseMaterialParams_.setVgn(4.0); - fineMaterialParams_.setVgAlpha(0.0005); - coarseMaterialParams_.setVgAlpha(0.0015); + baseParamsFine.setVgn(4.0); + baseParamsCoarse.setVgAlpha(0.0005); + baseParamsCoarse.setVgn(4.0); + baseParamsCoarse.setVgAlpha(0.0015); + + baseParamsFine.setKrRegardsSnr(false); + baseParamsCoarse.setKrRegardsSnr(false); - coarseMaterialParams_.setKrRegardsSnr(false); - fineMaterialParams_.setKrRegardsSnr(false); + pcKrSwCurveFine_ = std::make_unique(baseParamsFine); + pcKrSwCurveCoarse_ = std::make_unique(baseParamsCoarse); } /*! @@ -142,33 +144,32 @@ public: } /*! - * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.). + * \brief Returns the fluid-matrix interaction law (kr-sw, pc-sw, etc.). * * \param element The current element * \param scv The sub-control volume inside the element. * \param elemSol The solution at the dofs connected to the element. - * \return the material parameters object */ template - const MaterialLawParams& materialLawParams(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol) const + auto fluidMatrixInteraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const { - return materialLawParamsAtPos(scv.dofPosition()); + return fluidMatrixInteractionAtPos(scv.dofPosition()); } /*! - * \brief Returns the parameter object for the capillary-pressure/ - * saturation material law + * \brief Returns the fluid-matrix interaction law (kr-sw, pc-sw, etc.). * * \param globalPos The global position */ - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const + auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const { + assert(pcKrSwCurveFine_ && pcKrSwCurveCoarse_); if (isFineMaterial_(globalPos)) - return fineMaterialParams_; + return makeFluidMatrixInteraction(*pcKrSwCurveFine_); else - return coarseMaterialParams_; + return makeFluidMatrixInteraction(*pcKrSwCurveCoarse_); } struct MaxSaturations @@ -235,8 +236,11 @@ public: trapSats_[dofIdxGlobal].trappedSatN = (maxSats_[dofIdxGlobal].MaxSatN)/(1.0 + landParameter*maxSats_[dofIdxGlobal].MaxSatN); - coarseMaterialParams_.setTrappedSatN(trapSats_[dofIdxGlobal].trappedSatN); - fineMaterialParams_.setTrappedSatN(trapSats_[dofIdxGlobal].trappedSatN); + // TODO this makes no sense, you are overwriting the value with each dof. + // you probably want a vector of pcSwcurves for each dof? + // or update the trapped sn just before you return it from fluidMatrixInteractionAtPos? + pcKrSwCurveFine_->updateTappedSn(trapSats_[dofIdxGlobal].trappedSatN); + pcKrSwCurveCoarse_->updateTappedSn(trapSats_[dofIdxGlobal].trappedSatN); } } } @@ -258,8 +262,8 @@ private: Scalar finePorosity_; Scalar coarsePorosity_; - MaterialLawParams fineMaterialParams_; - MaterialLawParams coarseMaterialParams_; + std::unique_ptr pcKrSwCurveFine_; + std::unique_ptr pcKrSwCurveCoarse_; Scalar eps_; }; -- GitLab From f674b23ece8331069731a080de9fac3f16a25dd6 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sat, 31 Oct 2020 00:43:58 +0100 Subject: [PATCH 03/15] [efm][1p2c_2p_2p2c] Adapt to new material laws --- lecture/efm/1p2c_2p_2p2c/exercise3.input | 42 ++++++++-------- .../efm/1p2c_2p_2p2c/lens1p2cspatialparams.hh | 9 ++-- .../efm/1p2c_2p_2p2c/lens2pspatialparams.hh | 48 +++++++------------ lecture/efm/1p2cvs2p/exercise1.input | 25 +++++----- lecture/efm/2p/lens2pexercise2.input | 28 ++++++----- 5 files changed, 72 insertions(+), 80 deletions(-) diff --git a/lecture/efm/1p2c_2p_2p2c/exercise3.input b/lecture/efm/1p2c_2p_2p2c/exercise3.input index 902a0be..3cf360d 100755 --- a/lecture/efm/1p2c_2p_2p2c/exercise3.input +++ b/lecture/efm/1p2c_2p_2p2c/exercise3.input @@ -1,30 +1,32 @@ [TimeLoop] -MaxTimeStepSize = 5.0e6 # maximal time step size [s] -TEnd = 5.0e7 # end time of the simulation [s] -DtInitial = 0.0625 # initial time step for the simulation [s] +MaxTimeStepSize = 5.0e6 # maximal time step size [s] +TEnd = 5.0e7 # end time of the simulation [s] +DtInitial = 0.0625 # initial time step for the simulation [s] EpisodeLength = 5.0e6 [Grid] -UpperRight = 4 2 # upper right corner coordinates (x,y) [m] -Cells = 40 20 # grid resolution in (x,y) direction [-] +UpperRight = 4 2 # upper right corner coordinates (x,y) [m] +Cells = 40 20 # grid resolution in (x,y) direction [-] -[SpatialParams] -FinePermeability = 3.1e-11 # intrinsic permeability of the fine porous medium [m^2] -CoarsePermeability = 4.6e-9 # intrinsic permeability of the coarse porous medium [m^2] -FinePorosity = 0.2 # porosity of the fine porous medium [-] -CoarsePorosity = 0.6 # porosity of the coarse porous medium [-] -FineBrooksCoreyLambda = 3.5 # pore size distribution parameter for the Brooks-Corey capillary pressure-saturation relationship in the fine soil [-] -FineBrooksCoreyEntryPressure = 400 # entry pressure for the Brooks-Corey capillary pressure-saturation relationship in the fine soil [Pa] -CoarseBrooksCoreyLambda = 2.0 # pore size distribution parameter for the Brooks-Corey capillary pressure-saturation relationship in the coarse soil [-] -CoarseBrooksCoreyEntryPressure = 200 # entry pressure for the Brooks-Corey capillary pressure-saturation relationship in the coarse soil [Pa] -FineResidualSaturationWetting = 0.05 # residual saturation of the wetting phase in the fine soil [-] -FineResidualSaturationNonwetting = 0.3 # residual saturation of the nonwetting phase in the fine soil [-] -CoarseResidualSaturationWetting = 0.05 # residual saturation of the wetting phase in the coarse soil [-] -CoarseResidualSaturationNonwetting = 0.1 # residual saturation of the nonwetting phase in the coarse soil [-] +[SpatialParams.Fine] +Permeability = 3.1e-11 # intrinsic permeability of the fine porous medium [m^2] +Porosity = 0.2 # porosity of the fine porous medium [-] +BrooksCoreyLambda = 3.5 # pore size distribution parameter for the Brooks-Corey capillary pressure-saturation relationship in the fine soil [-] +BrooksCoreyPcEntry = 400 # entry pressure for the Brooks-Corey capillary pressure-saturation relationship in the fine soil [Pa] +Swr = 0.05 # residual saturation of the wetting phase in the fine soil [-] +Snr = 0.3 # residual saturation of the nonwetting phase in the fine soil [-] + +[SpatialParams.Coarse] +Permeability = 4.6e-9 # intrinsic permeability of the coarse porous medium [m^2] +Porosity = 0.6 # porosity of the coarse porous medium [-] +BrooksCoreyLambda = 2.0 # pore size distribution parameter for the Brooks-Corey capillary pressure-saturation relationship in the coarse soil [-] +BrooksCoreyPcEntry = 200 # entry pressure for the Brooks-Corey capillary pressure-saturation relationship in the coarse soil [Pa] +Swr = 0.05 # residual saturation of the wetting phase in the coarse soil [-] +Snr = 0.1 # residual saturation of the nonwetting phase in the coarse soil [-] [Boundary] -LowerPressure = 2.0e5 # Dirichlet pressure value for the boundary condition at the lower boundary [Pa] -UpperPressure = 4.0e5 # Dirichlet pressure value for the boundary condition at the upper boundary [Pa] +LowerPressure = 2.0e5 # Dirichlet pressure value for the boundary condition at the lower boundary [Pa] +UpperPressure = 4.0e5 # Dirichlet pressure value for the boundary condition at the upper boundary [Pa] [Problem] EnableGravity = false diff --git a/lecture/efm/1p2c_2p_2p2c/lens1p2cspatialparams.hh b/lecture/efm/1p2c_2p_2p2c/lens1p2cspatialparams.hh index 1888161..1d9b051 100644 --- a/lecture/efm/1p2c_2p_2p2c/lens1p2cspatialparams.hh +++ b/lecture/efm/1p2c_2p_2p2c/lens1p2cspatialparams.hh @@ -80,17 +80,16 @@ public: Lens1p2cSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) { - lensLowerLeft_[0] = 0.0; lensLowerLeft_[1] = 0.0; lensUpperRight_[0]= 4.1; lensUpperRight_[1]= 1.0; - lensPorosity_ = getParam("SpatialParams.FinePorosity"); - outerPorosity_ = getParam("SpatialParams.CoarsePorosity"); + lensPorosity_ = getParam("SpatialParams.Fine.Porosity"); + outerPorosity_ = getParam("SpatialParams.Coarse.Porosity"); - lensK_ = getParam("SpatialParams.FinePermeability"); - outerK_ = getParam("SpatialParams.CoarsePermeability"); + lensK_ = getParam("SpatialParams.Fine.Permeability"); + outerK_ = getParam("SpatialParams.Coarse.Permeability"); } /*! diff --git a/lecture/efm/1p2c_2p_2p2c/lens2pspatialparams.hh b/lecture/efm/1p2c_2p_2p2c/lens2pspatialparams.hh index ced425e..0921a38 100644 --- a/lecture/efm/1p2c_2p_2p2c/lens2pspatialparams.hh +++ b/lecture/efm/1p2c_2p_2p2c/lens2pspatialparams.hh @@ -26,8 +26,7 @@ #ifndef DUMUX_LENS2P_SPATIALPARAMS_HH #define DUMUX_LENS2P_SPATIALPARAMS_HH -#include -#include +#include #include namespace Dumux { @@ -49,11 +48,10 @@ class Lens2pSpatialParams using ParentType = FVSpatialParams; static constexpr int dimWorld = GridView::dimensionworld; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - using EffectiveLaw = RegularizedBrooksCorey; + + using PcKrSw = FluidMatrix::BrooksCoreyDefault; public: - using MaterialLaw = EffToAbsLaw; - using MaterialLawParams = typename MaterialLaw::Params; using PermeabilityType = Scalar; /*! @@ -61,27 +59,17 @@ public: */ Lens2pSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) + , pcKrSwFine_("SpatialParams.Fine") + , pcKrSwCoarse_("SpatialParams.Coarse") { lensLowerLeft_ = {0.0, 0.0}; lensUpperRight_= {4.0, 1.0}; - lensPorosity_ = getParam ("SpatialParams.FinePorosity"); - outerPorosity_ = getParam("SpatialParams.CoarsePorosity"); - - lensK_ = getParam("SpatialParams.FinePermeability"); - outerK_ = getParam("SpatialParams.CoarsePermeability"); - - // residual saturations - lensMaterialParams_.setSwr( getParam("SpatialParams.FineResidualSaturationWetting") ); - lensMaterialParams_.setSnr( getParam("SpatialParams.FineResidualSaturationNonwetting") ); - outerMaterialParams_.setSwr( getParam("SpatialParams.CoarseResidualSaturationWetting") ); - outerMaterialParams_.setSnr( getParam("SpatialParams.CoarseResidualSaturationNonwetting") ); + lensPorosity_ = getParam ("SpatialParams.Fine.Porosity"); + outerPorosity_ = getParam("SpatialParams.Coarse.Porosity"); - // parameters for the Brooks-Corey law - lensMaterialParams_.setPe( getParam("SpatialParams.FineBrooksCoreyEntryPressure") ); - lensMaterialParams_.setLambda( getParam("SpatialParams.FineBrooksCoreyLambda") ); - outerMaterialParams_.setPe( getParam("SpatialParams.CoarseBrooksCoreyEntryPressure") ); - outerMaterialParams_.setLambda( getParam("SpatialParams.CoarseBrooksCoreyLambda") ); + lensK_ = getParam("SpatialParams.Fine.Permeability"); + outerK_ = getParam("SpatialParams.Coarse.Permeability"); } /*! @@ -121,20 +109,20 @@ public: } /*! - * \brief return the brooks-corey context depending on the position + * \brief Returns the fluid-matrix interaction law * * \param element The current element * \param scv The sub-control volume inside the element. * \param elemSol The solution at the dofs connected to the element. */ template - const MaterialLawParams& materialLawParams(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol) const + auto fluidMatrixInteraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const { if (isInLens_(scv.dofPosition())) - return lensMaterialParams_; - return outerMaterialParams_; + return makeFluidMatrixInteraction(pcKrSwFine_); + return makeFluidMatrixInteraction(pcKrSwCoarse_); } /*! @@ -178,9 +166,9 @@ private: Scalar outerK_; Scalar lensPorosity_; Scalar outerPorosity_; - MaterialLawParams lensMaterialParams_; - MaterialLawParams outerMaterialParams_; - Scalar eps_ = 1e-6; + const PcKrSw pcKrSwFine_; + const PcKrSw pcKrSwCoarse_; + const Scalar eps_{1e-6}; }; } // end namespace Dumux diff --git a/lecture/efm/1p2cvs2p/exercise1.input b/lecture/efm/1p2cvs2p/exercise1.input index b67974d..b36f041 100644 --- a/lecture/efm/1p2cvs2p/exercise1.input +++ b/lecture/efm/1p2cvs2p/exercise1.input @@ -12,21 +12,22 @@ UpperRight = 5 4 # upper right corner coordinates (x, Cells = 40 64 # grid resolution in (x,y) direction [-] [SpatialParams] -FinePermeability = 3.1e-11 # intrinsic permeability of the fine porous medium [m^2] -CoarsePermeability = 3.1e-10 # intrinsic permeability of the coarse porous medium [m^2] -FinePorosity = 0.1 # porosity of the fine porous medium [-] -CoarsePorosity = 0.2 # porosity of the coarse porous medium [-] + +Fine.Permeability = 3.1e-11 # intrinsic permeability of the fine porous medium [m^2] +Coarse.Permeability = 3.1e-10 # intrinsic permeability of the coarse porous medium [m^2] +Fine.Porosity = 0.1 # porosity of the fine porous medium [-] +Coarse.Porosity = 0.2 # porosity of the coarse porous medium [-] ######## Parameters only relevant for two-phase simulations: ######## -FineBrooksCoreyLambda = 3.5 # pore size distribution parameter for the Brooks-Corey capillary pressure - saturation relationship in the fine soil [-] -FineBrooksCoreyEntryPressure = 400 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship in the fine soil [Pa] -CoarseBrooksCoreyLambda = 2.0 # pore size distribution parameter for the Brooks-Corey capillary pressure - saturation relationship in the coarse soil [-] -CoarseBrooksCoreyEntryPressure = 200 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship in the coarse soil [Pa] -FineResidualSaturationWetting = 0.05 # residual saturation of the wetting phase in the fine soil [-] -FineResidualSaturationNonwetting = 0.3 # residual saturation of the nonwetting phase in the fine soil [-] -CoarseResidualSaturationWetting = 0.05 # residual saturation of the wetting phase in the coarse soil [-] -CoarseResidualSaturationNonwetting = 0.1 # residual saturation of the nonwetting phase in the coarse soil [-] +Fine.BrooksCoreyLambda = 3.5 # pore size distribution parameter for the Brooks-Corey capillary pressure-saturation relationship in the fine soil [-] +Fine.BrooksCoreyPcEntry = 400 # entry pressure for the Brooks-Corey capillary pressure-saturation relationship in the fine soil [Pa] +Fine.Swr = 0.05 # residual saturation of the wetting phase in the fine soil [-] +Fine.Snr = 0.3 # residual saturation of the nonwetting phase in the fine soil [-] +Coarse.BrooksCoreyLambda = 2.0 # pore size distribution parameter for the Brooks-Corey capillary pressure-saturation relationship in the coarse soil [-] +Coarse.BrooksCoreyPcEntry = 200 # entry pressure for the Brooks-Corey capillary pressure-saturation relationship in the coarse soil [Pa] +Coarse.Swr = 0.05 # residual saturation of the wetting phase in the coarse soil [-] +Coarse.Snr = 0.1 # residual saturation of the nonwetting phase in the coarse soil [-] ######## [Boundary] diff --git a/lecture/efm/2p/lens2pexercise2.input b/lecture/efm/2p/lens2pexercise2.input index 2c66ce1..3c75032 100644 --- a/lecture/efm/2p/lens2pexercise2.input +++ b/lecture/efm/2p/lens2pexercise2.input @@ -9,19 +9,21 @@ LowerLeft = 0 0 # lower left corner coordinates (x,y UpperRight = 3 2 # upper right corner coordinates (x,y) [m] Cells = 75 50 # grid resolution in (x,y) direction [-] -[SpatialParams] -FinePermeability = 9e-12 # intrinsic permeability of the fine porous medium [m^2] -CoarsePermeability = 4.6e-10 # intrinsic permeability of the coarse porous medium [m^2] -FinePorosity = 0.38 # porosity of the fine porous medium [-] -CoarsePorosity = 0.40 # porosity of the coarse porous medium [-] -FineBrooksCoreyLambda = 3.5 # pore size distribution parameter for the Brooks-Corey capillary pressure - saturation relationship in the fine soil [-] -FineBrooksCoreyEntryPressure = 500 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship in the fine soil [Pa] -CoarseBrooksCoreyLambda = 2.0 # pore size distribution parameter for the Brooks-Corey capillary pressure - saturation relationship in the coarse soil [-] -CoarseBrooksCoreyEntryPressure = 200 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship in the coarse soil [Pa] -FineResidualSaturationWetting = 0.18 # residual saturation of the wetting phase in the fine soil [-] -FineResidualSaturationNonwetting = 0.0 # residual saturation of the nonwetting phase in the fine soil [-] -CoarseResidualSaturationWetting = 0.05 # residual saturation of the wetting phase in the coarse soil [-] -CoarseResidualSaturationNonwetting = 0.0 # residual saturation of the nonwetting phase in the coarse soil [-] +[SpatialParams.Fine] +Permeability = 9e-12 # intrinsic permeability of the fine porous medium [m^2] +Porosity = 0.38 # porosity of the fine porous medium [-] +BrooksCoreyLambda = 3.5 # pore size distribution parameter for the Brooks-Corey capillary pressure - saturation relationship in the fine soil [-] +BrooksCoreyPcEntry = 500 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship in the fine soil [Pa] +Swr = 0.18 # residual saturation of the wetting phase in the fine soil [-] +Snr = 0.0 # residual saturation of the nonwetting phase in the fine soil [-] + +[SpatialParams.Coarse] +Permeability = 4.6e-10 # intrinsic permeability of the coarse porous medium [m^2] +Porosity = 0.40 # porosity of the coarse porous medium [-] +BrooksCoreyLambda = 2.0 # pore size distribution parameter for the Brooks-Corey capillary pressure - saturation relationship in the coarse soil [-] +BrooksCoreyPcEntry = 200 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship in the coarse soil [Pa] +Swr = 0.05 # residual saturation of the wetting phase in the coarse soil [-] +Snr = 0.0 # residual saturation of the nonwetting phase in the coarse soil [-] [Boundary] LowerPressure = 1.19612e5 # Dirichlet pressure value for the boundary condition at the lower boundary [Pa] -- GitLab From 61c7d025fab7bf2ea4e35d8667c12d2193d7f659 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Sat, 31 Oct 2020 19:37:11 +0100 Subject: [PATCH 04/15] [co2] adapt to new material laws --- .../mm/co2plume/co2plumeshapeexercise.input | 6 ++-- .../co2plumeshapespatialparameters.hh | 34 +++++-------------- 2 files changed, 10 insertions(+), 30 deletions(-) diff --git a/lecture/mm/co2plume/co2plumeshapeexercise.input b/lecture/mm/co2plume/co2plumeshapeexercise.input index 8eb6029..26e4965 100644 --- a/lecture/mm/co2plume/co2plumeshapeexercise.input +++ b/lecture/mm/co2plume/co2plumeshapeexercise.input @@ -21,12 +21,10 @@ Salinity = 0.1 Permeability = 1e-13 # [m^2] intrinsic permeability DipAngle = 0.0 # [deg] dip angle for the domain Porosity = 0.3 # porosity - -[MaterialLaw] Swr = 0.2 # [-] residual wetting phase sat. Snr = 0.2 # [-] residual nonwetting phase sat. -Pe = 5e3 # [Pa] capillary entry pressure -Lambda = 2 # [-] Brooks Corey parameter +BrooksCoreyPcEntry = 5e3 # [Pa] capillary entry pressure +BrooksCoreyLambda = 2 # [-] Brooks Corey parameter [Problem] Name = co2plumeshape # the name of the output files diff --git a/lecture/mm/co2plume/co2plumeshapespatialparameters.hh b/lecture/mm/co2plume/co2plumeshapespatialparameters.hh index 049e253..509bafb 100644 --- a/lecture/mm/co2plume/co2plumeshapespatialparameters.hh +++ b/lecture/mm/co2plume/co2plumeshapespatialparameters.hh @@ -28,8 +28,7 @@ #include #include -#include -#include +#include #include @@ -54,12 +53,9 @@ class PlumeShapeSpatialParams using ParentType = FVSpatialParams>; using GlobalPosition = typename SubControlVolume::GlobalPosition; - - using EffectiveLaw = RegularizedBrooksCorey; + using PcKrSwCurve = FluidMatrix::BrooksCoreyDefault; public: - using MaterialLaw = EffToAbsLaw; - using MaterialLawParams = typename MaterialLaw::Params; using PermeabilityType = Scalar; /*! @@ -67,20 +63,13 @@ public: */ PlumeShapeSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) + , pcKrSwCurve_("SpatialParams") { // intrinsic permeability permeability_ = getParam("SpatialParams.Permeability"); // porosity porosity_ = getParam("SpatialParams.Porosity"); - - // residual saturations - materialParams_.setSwr(getParam("MaterialLaw.Swr")); - materialParams_.setSnr(getParam("MaterialLaw.Snr")); - - // parameters for the Brooks-Corey law - materialParams_.setPe(getParam("MaterialLaw.Pe")); - materialParams_.setLambda(getParam("MaterialLaw.Lambda")); } /*! @@ -96,9 +85,7 @@ public: PermeabilityType permeability(const Element& element, const SubControlVolume& scv, const ElementSolution& elemSol) const - { - return permeability_; - } + { return permeability_; } /*! * \brief Returns the porosity \f$[-]\f$ @@ -112,9 +99,7 @@ public: Scalar porosity(const Element& element, const SubControlVolume& scv, const ElementSolution& elemSol) const - { - return porosity_; - } + { return porosity_; } /*! @@ -123,10 +108,8 @@ public: * \return the material parameters object * \param globalPos The position of the center of the element */ - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const - { - return materialParams_; - } + const auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const + { return pcKrSwCurve_; } /*! * \brief Function for defining which phase is to be considered as the wetting phase. @@ -142,8 +125,7 @@ private: Scalar permeability_; Scalar porosity_; - MaterialLawParams materialParams_; - + const PcKrSwCurve pcKrSwCurve_; }; } // end namespace Dumux -- GitLab From cdbfcec5a25fc1d3ada131697043282927e47b2b Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Mon, 2 Nov 2020 09:24:14 +0100 Subject: [PATCH 05/15] [mm][columnxylene] adapt to new material laws --- .../columnxylene/columnxyleneexercise.input | 29 +++++++- .../columnxylene/columnxylenespatialparams.hh | 66 +++++-------------- 2 files changed, 41 insertions(+), 54 deletions(-) diff --git a/lecture/mm/columnxylene/columnxyleneexercise.input b/lecture/mm/columnxylene/columnxyleneexercise.input index 026202a..a20898a 100644 --- a/lecture/mm/columnxylene/columnxyleneexercise.input +++ b/lecture/mm/columnxylene/columnxyleneexercise.input @@ -9,9 +9,32 @@ File = ./grids/column.dgf [SpatialParams] Permeability = 1.4e-11 PermeabilityExtDomain = 1.e-8 -BetaNw = 1.83 -BetaGn = 2.2 -BetaGw = 1.0 +ThreePNAPLAdsorptionRhoBulk = 1500.0 +ThreePNAPLAdsorptionKdNAPL = 0.0 + +[SpatialParams.FineMaterial] +Swr = 0.12 +Snr = 0.10 +Sgr = 0.01 +ParkerVanGenuchtenAlpha = 0.0005 +ParkerVanGenuchtenN = 4.0 +ParkerVanGenuchtenRegardSnrForKrn = true +ParkerVanGenuchtenBetaNw = 1.83 +ParkerVanGenuchtenBetaGn = 2.2 +ParkerVanGenuchtenBetaGw = 1.0 +VanGenuchtenConstantRegularization = true + +[SpatialParams.CoarseMaterial] +Swr = 0.12 +Snr = 0.10 +Sgr = 0.01 +ParkerVanGenuchtenAlpha = 0.5 +ParkerVanGenuchtenN = 4.0 +ParkerVanGenuchtenRegardSnrForKrn = false +ParkerVanGenuchtenBetaNw = 1.83 +ParkerVanGenuchtenBetaGn = 2.2 +ParkerVanGenuchtenBetaGw = 1.0 +VanGenuchtenConstantRegularization = true [Problem] Name = columnxylol diff --git a/lecture/mm/columnxylene/columnxylenespatialparams.hh b/lecture/mm/columnxylene/columnxylenespatialparams.hh index 1d26bc7..b4b6eb1 100644 --- a/lecture/mm/columnxylene/columnxylenespatialparams.hh +++ b/lecture/mm/columnxylene/columnxylenespatialparams.hh @@ -27,9 +27,8 @@ #include #include #include -#include -#include -#include +#include +#include namespace Dumux { @@ -50,11 +49,10 @@ class ColumnSpatialParams ColumnSpatialParams>; using GlobalPosition = typename SubControlVolume::GlobalPosition; - using EffectiveLaw = RegularizedParkerVanGen3P; + using ThreePhasePcKrSw = FluidMatrix::ParkerVanGenuchten3PDefault; + using AdsorptionModel = FluidMatrix::ThreePNAPLAdsorption; public: - using MaterialLaw = EffToAbsLaw; - using MaterialLawParams = typename MaterialLaw::Params; using PermeabilityType = Scalar; /*! @@ -62,6 +60,9 @@ public: */ ColumnSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) + , fineThreePhasePcKrSw_("SpatialParams.FineMaterial") + , coarseThreePhasePcKrSw_("SpatialParams.CoarseMaterial") + , adsorption_("SpatialParams") { // intrinsic permeabilities fineK_ = getParam("SpatialParams.Permeability"); @@ -70,41 +71,6 @@ public: // porosities finePorosity_ = 0.46; coarsePorosity_ = 0.46; - - // residual saturations - fineMaterialParams_.setSwr(0.12); - fineMaterialParams_.setSnr(0.10); - fineMaterialParams_.setSgr(0.01); - coarseMaterialParams_.setSwr(0.12); - coarseMaterialParams_.setSnr(0.10); - coarseMaterialParams_.setSgr(0.01); - - // parameters for the 3phase van Genuchten law - fineMaterialParams_.setVgAlpha(0.0005); - coarseMaterialParams_.setVgAlpha(0.5); - fineMaterialParams_.setVgn(4.0); - coarseMaterialParams_.setVgn(4.0); - - coarseMaterialParams_.setKrRegardsSnr(false); - fineMaterialParams_.setKrRegardsSnr(true); - - // parameters for the scaling of capillary pressure (GW = 1); - fineMaterialParams_.setBetaNw(getParam("SpatialParams.BetaNw")); - fineMaterialParams_.setBetaGn(getParam("SpatialParams.BetaGn")); - fineMaterialParams_.setBetaGw(getParam("SpatialParams.BetaGw")); - coarseMaterialParams_.setBetaNw(getParam("SpatialParams.BetaNw")); - coarseMaterialParams_.setBetaGn(getParam("SpatialParams.BetaGn")); - coarseMaterialParams_.setBetaGw(getParam("SpatialParams.BetaGw")); - - // parameters for adsorption - coarseMaterialParams_.setKdNAPL(0.); - coarseMaterialParams_.setRhoBulk(1500.); - fineMaterialParams_.setKdNAPL(0.); - fineMaterialParams_.setRhoBulk(1500.); - - // use capping as regularization - coarseMaterialParams_.useConstRegularization(true); - fineMaterialParams_.useConstRegularization(true); } /*! @@ -168,19 +134,16 @@ public: * \param elemSol The solution at the dofs connected to the element */ template - const MaterialLawParams& materialLawParams(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol) const + auto fluidMatrixInteraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const { const auto& globalPos = scv.dofPosition(); if (isFineMaterial_(globalPos)) - return fineMaterialParams_; - - else - return coarseMaterialParams_; + return makeFluidMatrixInteraction(fineThreePhasePcKrSw_, adsorption_); + return makeFluidMatrixInteraction(coarseThreePhasePcKrSw_, adsorption_); } - private: bool isFineMaterial_(const GlobalPosition &pos) const { @@ -200,8 +163,9 @@ private: Scalar fineHeatCap_; Scalar coarseHeatCap_; - MaterialLawParams fineMaterialParams_; - MaterialLawParams coarseMaterialParams_; + const ThreePhasePcKrSw fineThreePhasePcKrSw_; + const ThreePhasePcKrSw coarseThreePhasePcKrSw_; + const AdsorptionModel adsorption_; static constexpr Scalar eps_ = 1.5e-7; }; -- GitLab From f1b9f520f522ee798a9dfd97cb86382387df49b9 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Mon, 2 Nov 2020 09:42:29 +0100 Subject: [PATCH 06/15] [mm][fracture] adapt to new material laws --- lecture/mm/fractures/fracture_exercise.input | 45 +++++++++------- lecture/mm/fractures/fracturespatialparams.hh | 26 +++------ lecture/mm/fractures/matrixspatialparams.hh | 54 +++++++------------ 3 files changed, 53 insertions(+), 72 deletions(-) diff --git a/lecture/mm/fractures/fracture_exercise.input b/lecture/mm/fractures/fracture_exercise.input index d74facd..9fe6a9c 100644 --- a/lecture/mm/fractures/fracture_exercise.input +++ b/lecture/mm/fractures/fracture_exercise.input @@ -18,28 +18,33 @@ InjectionLength = 10 # [m] The x-extent of the injection well [Matrix] Problem.Name = matrix -SpatialParams.Permeability = 1e-14 -SpatialParams.Porosity = 0.1 -SpatialParams.VGAlpha = 1e-3 -SpatialParams.VGN = 2 -SpatialParams.Snr = 0.0 -SpatialParams.Swr = 0.0 -SpatialParams.OverburdenPermeability = 5e-10 -SpatialParams.OverburdenPorosity = 0.3 -SpatialParams.OverburdenVGAlpha = 5e-3 -SpatialParams.OverburdenVGN = 2 -SpatialParams.OverburdenSnr = 0.0 -SpatialParams.OverburdenSwr = 0.0 - [Fracture] Problem.Name = fractures -SpatialParams.Aperture = 1e-1 -SpatialParams.Permeability = 1e-7 -SpatialParams.Porosity = 0.85 -SpatialParams.VGAlpha = 1e-2 -SpatialParams.VGN = 2 -SpatialParams.Snr = 0.0 -SpatialParams.Swr = 0.0 + +[SpatialParams] +Aperture = 1e-1 +Permeability = 1e-14 +Porosity = 0.1 +VanGenuchtenAlpha = 1e-3 +VanGenuchtenN = 2 +Snr = 0.0 +Swr = 0.0 + +[SpatialParams.Overburden] +Permeability = 5e-10 +Porosity = 0.3 +VanGenuchtenAlpha = 5e-3 +VanGenuchtenN = 2 +Snr = 0.0 +Swr = 0.0 + +[SpatialParams.Fracture] +Permeability = 1e-7 +Porosity = 0.85 +VanGenuchtenAlpha = 1e-2 +VanGenuchtenN = 2 +Snr = 0.0 +Swr = 0.0 [Newton] MaxRelativeShift = 1e-6 diff --git a/lecture/mm/fractures/fracturespatialparams.hh b/lecture/mm/fractures/fracturespatialparams.hh index 8c4ed6f..8f6277b 100644 --- a/lecture/mm/fractures/fracturespatialparams.hh +++ b/lecture/mm/fractures/fracturespatialparams.hh @@ -25,8 +25,7 @@ #define DUMUX_COURSE_FRACTURESEXERCISE_FRACTURE_SPATIALPARAMS_HH #include -#include -#include +#include namespace Dumux { @@ -47,7 +46,7 @@ class FractureSpatialParams using GlobalPosition = typename Element::Geometry::GlobalCoordinate; // use a regularized van-genuchten material law - using EffectiveLaw = RegularizedVanGenuchten; + using PcKrSw = FluidMatrix::VanGenuchtenDefault; // we identify those fractures as barriers, that have a domain marker // of 2. This is what is set in the grid file (see grids/complex.geo) @@ -57,23 +56,14 @@ public: //! export the type used for permeabilities using PermeabilityType = Scalar; - //! export the material law and parameters used - using MaterialLaw = EffToAbsLaw< EffectiveLaw >; - using MaterialLawParams = typename MaterialLaw::Params; - //! the constructor FractureSpatialParams(std::shared_ptr fvGridGeometry, const std::string& paramGroup) : ParentType(fvGridGeometry) + , pcKrSwFracture_("SpatialParams.Fracture") { - porosity_ = getParamFromGroup(paramGroup, "SpatialParams.Porosity"); - permeability_ = getParamFromGroup(paramGroup, "SpatialParams.Permeability"); - - // set the material law parameters - materialLawParams_.setSnr(getParamFromGroup(paramGroup, "SpatialParams.Snr")); - materialLawParams_.setSnr(getParamFromGroup(paramGroup, "SpatialParams.Swr")); - materialLawParams_.setVgAlpha(getParamFromGroup(paramGroup, "SpatialParams.VGAlpha")); - materialLawParams_.setVgn(getParamFromGroup(paramGroup, "SpatialParams.VGN")); + porosity_ = getParamFromGroup(paramGroup, "SpatialParams.Fracture.Porosity"); + permeability_ = getParamFromGroup(paramGroup, "SpatialParams.Fracture.Permeability"); } //! Function for defining the (intrinsic) permeability \f$[m^2]\f$. @@ -91,8 +81,8 @@ public: { return porosity_; } //! Return the material law parameters - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const - { return materialLawParams_; } + auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const + { return pcKrSwFracture_; } //! Water is the wetting phase template< class FluidSystem > @@ -106,7 +96,7 @@ public: private: Scalar porosity_; PermeabilityType permeability_; - MaterialLawParams materialLawParams_; + const PcKrSw pcKrSwFracture_; }; } // end namespace Dumux diff --git a/lecture/mm/fractures/matrixspatialparams.hh b/lecture/mm/fractures/matrixspatialparams.hh index 737950e..29c37e2 100644 --- a/lecture/mm/fractures/matrixspatialparams.hh +++ b/lecture/mm/fractures/matrixspatialparams.hh @@ -25,8 +25,7 @@ #define DUMUX_COURSE_FRACTURESEXERCISE_MATRIX_SPATIALPARAMS_HH #include -#include -#include +#include namespace Dumux { @@ -46,61 +45,48 @@ class MatrixSpatialParams using GlobalPosition = typename Element::Geometry::GlobalCoordinate; // use a regularized van-genuchten material law - using EffectiveLaw = RegularizedVanGenuchten; + using PcKrSw = FluidMatrix::VanGenuchtenDefault; public: //! export the type used for permeabilities using PermeabilityType = Scalar; - //! export the material law and parameters used - using MaterialLaw = EffToAbsLaw< EffectiveLaw >; - using MaterialLawParams = typename MaterialLaw::Params; - //! the constructor MatrixSpatialParams(std::shared_ptr fvGridGeometry, const std::string& paramGroup) : ParentType(fvGridGeometry) + , pcKrSw_("SpatialParams") + , pcKrSwOverburden_("SpatialParams.Overburden") { porosity_ = getParamFromGroup(paramGroup, "SpatialParams.Porosity"); - porosityOverburden_ = getParamFromGroup(paramGroup, "SpatialParams.OverburdenPorosity"); permeability_ = getParamFromGroup(paramGroup, "SpatialParams.Permeability"); - permeabilityOverburden_ = getParamFromGroup(paramGroup, "SpatialParams.OverburdenPermeability"); - - // set the material law parameters - materialLawParams_.setSnr(getParamFromGroup(paramGroup, "SpatialParams.Snr")); - materialLawParams_.setSwr(getParamFromGroup(paramGroup, "SpatialParams.Swr")); - materialLawParams_.setVgAlpha(getParamFromGroup(paramGroup, "SpatialParams.VGAlpha")); - materialLawParams_.setVgn(getParamFromGroup(paramGroup, "SpatialParams.VGN")); - - materialLawParamsOverburden_.setSnr(getParamFromGroup(paramGroup, "SpatialParams.OverburdenSnr")); - materialLawParamsOverburden_.setSwr(getParamFromGroup(paramGroup, "SpatialParams.OverburdenSwr")); - materialLawParamsOverburden_.setVgAlpha(getParamFromGroup(paramGroup, "SpatialParams.OverburdenVGAlpha")); - materialLawParamsOverburden_.setVgn(getParamFromGroup(paramGroup, "SpatialParams.OverburdenVGN")); + porosityOverburden_ = getParamFromGroup(paramGroup, "SpatialParams.Overburden.Porosity"); + permeabilityOverburden_ = getParamFromGroup(paramGroup, "SpatialParams.Overburden.Permeability"); } //! Function for defining the (intrinsic) permeability \f$[m^2]\f$. PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const { - if (globalPos[1] > 35.0) - return permeabilityOverburden_; - else if ((globalPos[1] > 32)||(globalPos[1] < 30)) - return permeability_; - else return permeability_/10000.; + if (globalPos[1] > 35.0) + return permeabilityOverburden_; + else if ((globalPos[1] > 32)||(globalPos[1] < 30)) + return permeability_; + else return permeability_/10000.; } //! Return the porosity Scalar porosityAtPos(const GlobalPosition& globalPos) const { - if (globalPos[1] > 35.0) - return porosityOverburden_; - else if ((globalPos[1] > 32)||(globalPos[1] < 30)) - return porosity_; - else return 0.05; + if (globalPos[1] > 35.0) + return porosityOverburden_; + else if ((globalPos[1] > 32)||(globalPos[1] < 30)) + return porosity_; + else return 0.05; } //! Return the material law parameters - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const - { return globalPos[1] > 35.0 ? materialLawParamsOverburden_ : materialLawParams_; } + auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const + { return globalPos[1] > 35.0 ? pcKrSwOverburden_ : pcKrSw_; } //! Water is the wetting phase template< class FluidSystem > @@ -116,8 +102,8 @@ private: Scalar porosityOverburden_; PermeabilityType permeability_; PermeabilityType permeabilityOverburden_; - MaterialLawParams materialLawParams_; - MaterialLawParams materialLawParamsOverburden_; + const PcKrSw pcKrSw_; + const PcKrSw pcKrSwOverburden_; }; } // end namespace Dumux -- GitLab From 1888dc9eee61debf4c381571930767d4e5540d48 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Mon, 2 Nov 2020 09:59:19 +0100 Subject: [PATCH 07/15] [mm][heavyoil][sagd] adapt to new material laws --- lecture/mm/heavyoil/sagd/sagd.input | 16 ++++++ lecture/mm/heavyoil/sagd/spatialparams.hh | 63 +++++++---------------- 2 files changed, 34 insertions(+), 45 deletions(-) diff --git a/lecture/mm/heavyoil/sagd/sagd.input b/lecture/mm/heavyoil/sagd/sagd.input index 948a0b9..76255c8 100644 --- a/lecture/mm/heavyoil/sagd/sagd.input +++ b/lecture/mm/heavyoil/sagd/sagd.input @@ -21,3 +21,19 @@ RelTolerance = 1e-02 SolidDensity = 2650 SolidThermalConductivity = 2.8 SolidHeatCapacity = 850 + +[SpatialParams.Fine] +Swr = 0.1 +Snr = 0.09 +Sgr = 0.01 +ParkerVanGenuchtenN = 4.0 +ParkerVanGenuchtenAlpha = 0.0005 +ParkerVanGenuchtenRegardSnrForKrn = false + +[SpatialParams.Coarse] +setSwr = 0.1 +setSnr = 0.09 +setSgr = 0.01 +ParkerVanGenuchtenN =4.0 +ParkerVanGenuchtenAlpha = 0.0015 +ParkerVanGenuchtenRegardSnrForKrn = false diff --git a/lecture/mm/heavyoil/sagd/spatialparams.hh b/lecture/mm/heavyoil/sagd/spatialparams.hh index 620e037..a118507 100644 --- a/lecture/mm/heavyoil/sagd/spatialparams.hh +++ b/lecture/mm/heavyoil/sagd/spatialparams.hh @@ -25,13 +25,9 @@ #define DUMUX_SAGD_SPATIAL_PARAMS_HH #include -#include - #include - -#include -#include -#include +#include +#include namespace Dumux { @@ -55,12 +51,10 @@ class SagdSpatialParams using Element = typename GridView::template Codim<0>::Entity; using ParentType = FVSpatialParams>; - - using EffectiveLaw = RegularizedParkerVanGen3P; + using ThreePhasePcKrSw = FluidMatrix::ParkerVanGenuchten3PDefault; public: - using MaterialLaw = EffToAbsLaw; - using MaterialLawParams = typename MaterialLaw::Params; + using PermeabilityType = Scalar; /*! @@ -69,7 +63,10 @@ public: * \param gridView The grid view */ SagdSpatialParams(std::shared_ptr fvGridGeometry) - : ParentType(fvGridGeometry), eps_(1e-6) + : ParentType(fvGridGeometry) + , threePhasePcKrSwFine_("SpatialParams.Fine") + , threePhasePcKrSwCoarse_("SpatialParams.Coarse") + , eps_(1e-6) { layerBottom_ = 35.0; @@ -80,23 +77,6 @@ public: // porosities finePorosity_ = 0.10; coarsePorosity_ = 0.1; - - // 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); - - // parameters for the 3phase van Genuchten law - fineMaterialParams_.setVgn(4.0); - coarseMaterialParams_.setVgn(4.0); - fineMaterialParams_.setVgAlpha(0.0005); - coarseMaterialParams_.setVgAlpha(0.0015); - - coarseMaterialParams_.setKrRegardsSnr(false); - fineMaterialParams_.setKrRegardsSnr(false); } /*! @@ -148,12 +128,10 @@ public: * \return the material parameters object */ template - const MaterialLawParams& materialLawParams(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol) const - { - return materialLawParamsAtPos(scv.dofPosition()); - } + auto fluidMatrixInteraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { return fluidMatrixInteractionAtPos(scv.dofPosition()); } /*! * \brief Returns the parameter object for the capillary-pressure/ @@ -161,30 +139,25 @@ public: * * \param globalPos The global position */ - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const + auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const { if (isFineMaterial_(globalPos)) - return fineMaterialParams_; - else - return coarseMaterialParams_; + return threePhasePcKrSwFine_; + return threePhasePcKrSwCoarse_; } private: bool isFineMaterial_(const GlobalPosition &pos) const - { - return pos[dimWorld-1] > layerBottom_ - eps_; - }; + { return pos[dimWorld-1] > layerBottom_ - eps_; }; Scalar layerBottom_; - Scalar fineK_; Scalar coarseK_; - Scalar finePorosity_; Scalar coarsePorosity_; - MaterialLawParams fineMaterialParams_; - MaterialLawParams coarseMaterialParams_; + const ThreePhasePcKrSw threePhasePcKrSwFine_; + const ThreePhasePcKrSw threePhasePcKrSwCoarse_; Scalar eps_; }; -- GitLab From 62a573da57027126475d26e0d49b28adeba13f6d Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Mon, 2 Nov 2020 10:31:04 +0100 Subject: [PATCH 08/15] [mm][heavyoil][sagdcyclic] adapt to new material laws --- .../mm/heavyoil/sagdcyclic/sagd_cyclic.input | 16 ++++++ .../mm/heavyoil/sagdcyclic/spatialparams.hh | 56 +++++-------------- 2 files changed, 31 insertions(+), 41 deletions(-) diff --git a/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.input b/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.input index 30a59f2..a854c1e 100644 --- a/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.input +++ b/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.input @@ -21,3 +21,19 @@ MaxSteps = 8 SolidDensity = 2650 SolidThermalConductivity = 2.8 SolidHeatCapacity = 850 + +[SpatialParams.Fine] +Swr = 0.1 +Snr = 0.09 +Sgr = 0.01 +ParkerVanGenuchtenN = 4.0 +ParkerVanGenuchtenAlpha = 0.0005 +ParkerVanGenuchtenRegardSnrForKrn = false + +[SpatialParams.Coarse] +setSwr = 0.1 +setSnr = 0.09 +setSgr = 0.01 +ParkerVanGenuchtenN =4.0 +ParkerVanGenuchtenAlpha = 0.0015 +ParkerVanGenuchtenRegardSnrForKrn = false diff --git a/lecture/mm/heavyoil/sagdcyclic/spatialparams.hh b/lecture/mm/heavyoil/sagdcyclic/spatialparams.hh index 9e0a7c6..0fc1bd1 100644 --- a/lecture/mm/heavyoil/sagdcyclic/spatialparams.hh +++ b/lecture/mm/heavyoil/sagdcyclic/spatialparams.hh @@ -25,13 +25,9 @@ #define DUMUX_SAGD_SPATIAL_PARAMS_HH #include -#include - #include - -#include -#include -#include +#include +#include namespace Dumux { @@ -54,12 +50,9 @@ class SagdSpatialParams using Element = typename GridView::template Codim<0>::Entity; using ParentType = FVSpatialParams>; - - using EffectiveLaw = RegularizedParkerVanGen3P; + using ThreePhasePcKrSw = FluidMatrix::ParkerVanGenuchten3PDefault; public: - using MaterialLaw = EffToAbsLaw; - using MaterialLawParams = typename MaterialLaw::Params; using PermeabilityType = Scalar; /*! @@ -68,7 +61,10 @@ public: * \param gridView The grid view */ SagdSpatialParams(std::shared_ptr fvGridGeometry) - : ParentType(fvGridGeometry), eps_(1e-6) + : ParentType(fvGridGeometry) + , threePhasePcKrSwFine_("SpatialParams.Fine") + , threePhasePcKrSwCoarse_("SpatialParams.Coarse") + , eps_(1e-6) { layerBottom_ = 35.0; @@ -79,23 +75,6 @@ public: // porosities finePorosity_ = 0.10; coarsePorosity_ = 0.1; - - // 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); - - // parameters for the 3phase van Genuchten law - fineMaterialParams_.setVgn(4.0); - coarseMaterialParams_.setVgn(4.0); - fineMaterialParams_.setVgAlpha(0.0005); - coarseMaterialParams_.setVgAlpha(0.0015); - - coarseMaterialParams_.setKrRegardsSnr(false); - fineMaterialParams_.setKrRegardsSnr(false); } /*! @@ -147,12 +126,10 @@ public: * \return the material parameters object */ template - const MaterialLawParams& materialLawParams(const Element& element, + auto fluidMatrixInteraction(const Element& element, const SubControlVolume& scv, const ElementSolution& elemSol) const - { - return materialLawParamsAtPos(scv.dofPosition()); - } + { return fluidMatrixInteractionAtPos(scv.dofPosition()); } /*! * \brief Returns the parameter object for the capillary-pressure/ @@ -160,19 +137,16 @@ public: * * \param globalPos The global position */ - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const + auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const { if (isFineMaterial_(globalPos)) - return fineMaterialParams_; - else - return coarseMaterialParams_; + return threePhasePcKrSwFine_; + return threePhasePcKrSwCoarse_; } private: bool isFineMaterial_(const GlobalPosition &pos) const - { - return pos[dimWorld-1] > layerBottom_ - eps_; - }; + { return pos[dimWorld-1] > layerBottom_ - eps_; }; Scalar layerBottom_; @@ -182,8 +156,8 @@ private: Scalar finePorosity_; Scalar coarsePorosity_; - MaterialLawParams fineMaterialParams_; - MaterialLawParams coarseMaterialParams_; + const ThreePhasePcKrSw threePhasePcKrSwFine_; + const ThreePhasePcKrSw threePhasePcKrSwCoarse_; Scalar eps_; }; -- GitLab From 393ec3a5018074eb859abfc0b225a258a1f80ab8 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Mon, 2 Nov 2020 16:02:10 +0100 Subject: [PATCH 09/15] [mm][naplinfiltration] adapt to new material laws --- .../3p/naplinfiltration3p.input | 12 ++- lecture/mm/naplinfiltration/3p/problem.hh | 18 ++-- .../3p3c/naplinfiltration3p3c.input | 12 ++- lecture/mm/naplinfiltration/3p3c/problem.hh | 32 +++---- lecture/mm/naplinfiltration/spatialparams.hh | 94 +++++++++---------- 5 files changed, 87 insertions(+), 81 deletions(-) diff --git a/lecture/mm/naplinfiltration/3p/naplinfiltration3p.input b/lecture/mm/naplinfiltration/3p/naplinfiltration3p.input index 8b9da4e..770afb7 100644 --- a/lecture/mm/naplinfiltration/3p/naplinfiltration3p.input +++ b/lecture/mm/naplinfiltration/3p/naplinfiltration3p.input @@ -11,8 +11,16 @@ Cells = 250 10 # number of cells in (x, y) direction [-] [SpatialParams] permeability = 1.e-11 # [m^2] porosity = 0.40 # [-] -vanGenuchtenAlpha = 0.0005 # [-] -vanGenuchtenN = 4.0 # [-] +Swr = 0.12 +Snr = 0.07 +Sgr = 0.03 +ParkerVanGenuchtenN = 4.0 +ParkerVanGenuchtenAlpha = 0.0005 +ParkerVanGenuchtenBetaNw = 1.83 +ParkerVanGenuchtenBetaGn = 2.2 +ParkerVanGenuchtenBetaGw = 1.0 +ThreePNAPLAdsorptionKdNAPL = 0. +ThreePNAPLAdsorptionRhoBulk = 1500. [Problem] Name = naplinfiltration3p diff --git a/lecture/mm/naplinfiltration/3p/problem.hh b/lecture/mm/naplinfiltration/3p/problem.hh index e102c47..bb3fd9d 100644 --- a/lecture/mm/naplinfiltration/3p/problem.hh +++ b/lecture/mm/naplinfiltration/3p/problem.hh @@ -232,41 +232,37 @@ private: PrimaryVariables initial_(const GlobalPosition &globalPos) const { PrimaryVariables values(0.0); + Scalar swr=0.12, sgr=0.03; Scalar y = globalPos[1]; Scalar x = globalPos[0]; - const auto& materialLawParams = this->spatialParams().materialLawParamsAtPos(globalPos); - Scalar swr = materialLawParams.swr(); - Scalar sgr = materialLawParams.sgr(); - if(y > (-1e-3*x+5) - eps_) { Scalar pc = 9.81 * 1000.0 * (y - (-5e-4*x+5)); if (pc < 0.0) pc = 0.0; - Scalar sw = invertPcgw_(pc, this->spatialParams().materialLawParamsAtPos(globalPos)); + Scalar sw = invertPcgw_(pc, this->spatialParams().fluidMatrixInteractionAtPos(globalPos)); if(sw < swr) sw = swr; if(sw > 1.-sgr) sw = 1.-sgr; values[pressureIdx] = 1e5; values[swIdx] = sw; values[snIdx] = 0.; - }else + } + else { values[pressureIdx] = 1e5 + 9.81 * 1000.0 * ((-5e-4*x+5) - y); values[swIdx] = 1.-sgr; values[snIdx] = 0.; } - return values; } // small solver inverting the pc curve - template - static Scalar invertPcgw_(Scalar pcIn, const MaterialLawParams &pcParams) + template + static Scalar invertPcgw_(Scalar pcIn, const FMInteraction& fluidMatrixInteraction) { - using MaterialLaw = typename ParentType::SpatialParams::MaterialLaw; Scalar lower(0.0); Scalar upper(1.0); const unsigned int maxIter = 25; @@ -276,7 +272,7 @@ private: for (unsigned int k = 1; k <= maxIter; k++) { sw = 0.5*(upper + lower); - pcgw = MaterialLaw::pcgw(pcParams, sw); + pcgw = fluidMatrixInteraction.pcgw(sw, 0.0); const Scalar delta = std::abs(pcgw - pcIn); if (delta < bisLimit) return sw; diff --git a/lecture/mm/naplinfiltration/3p3c/naplinfiltration3p3c.input b/lecture/mm/naplinfiltration/3p3c/naplinfiltration3p3c.input index fb9f999..b5500a5 100644 --- a/lecture/mm/naplinfiltration/3p3c/naplinfiltration3p3c.input +++ b/lecture/mm/naplinfiltration/3p3c/naplinfiltration3p3c.input @@ -11,8 +11,16 @@ Cells = 250 10 # number of cells in (x, y) direction [-] [SpatialParams] permeability = 1.e-11 # [m^2] porosity = 0.40 # [-] -vanGenuchtenAlpha = 0.0005 # [-] -vanGenuchtenN = 4.0 # [-] +Swr = 0.12 +Snr = 0.07 +Sgr = 0.03 +ParkerVanGenuchtenN = 4.0 +ParkerVanGenuchtenAlpha = 0.0005 +ParkerVanGenuchtenBetaNw = 1.83 +ParkerVanGenuchtenBetaGn = 2.2 +ParkerVanGenuchtenBetaGw = 1.0 +ThreePNAPLAdsorptionKdNAPL = 0. +ThreePNAPLAdsorptionRhoBulk = 1500. [Output] PlotFluidMatrixInteractions = true diff --git a/lecture/mm/naplinfiltration/3p3c/problem.hh b/lecture/mm/naplinfiltration/3p3c/problem.hh index f843be7..a5c1291 100644 --- a/lecture/mm/naplinfiltration/3p3c/problem.hh +++ b/lecture/mm/naplinfiltration/3p3c/problem.hh @@ -225,38 +225,38 @@ private: { PrimaryVariables values(0.0); values.setState(wgPhaseOnly); + Scalar swr=0.12, sgr=0.03; Scalar y = globalPos[1]; Scalar x = globalPos[0]; - - const auto& materialLawParams = this->spatialParams().materialLawParamsAtPos(globalPos); - Scalar swr = materialLawParams.swr(); - Scalar sgr = materialLawParams.sgr(); - - if(y >(-1.E-3*x+5) - eps_) + if (y > (-1e-3*x + 5) - eps_) { - Scalar pc = 9.81 * 1000.0 * (y - (-5E-4*x+5)); + Scalar pc = 9.81 * 1000.0 * (y - (-5e-4*x + 5)); if (pc < 0.0) pc = 0.0; - Scalar sw = invertPcgw_(pc, this->spatialParams().materialLawParamsAtPos(globalPos)); - if (sw < swr) sw = swr; - if (sw > 1.-sgr) sw = 1.-sgr; + Scalar sw = invertPcgw_(pc, this->spatialParams().fluidMatrixInteractionAtPos(globalPos)); + if (sw < swr) + sw = swr; + if (sw > 1.-sgr) + sw = 1.-sgr; values[pressureIdx] = 1e5 ; values[switch1Idx] = sw; values[switch2Idx] = 0.0; - }else { - values[pressureIdx] = 1e5 + 9.81 * 1000.0 * ((-5E-4*x+5) - y); + } + else + { + values[pressureIdx] = 1e5 + 9.81 * 1000.0 * ((-5e-4*x + 5) - y); values[switch1Idx] = 1.-sgr; values[switch2Idx] = 0.0; } return values; } - template - static Scalar invertPcgw_(Scalar pcIn, const MaterialLawParams &pcParams) + // small solver inverting the pc curve + template + static Scalar invertPcgw_(Scalar pcIn, const FMInteraction& fluidMatrixInteraction) { - using MaterialLaw = typename ParentType::SpatialParams::MaterialLaw; Scalar lower(0.0); Scalar upper(1.0); const unsigned int maxIter = 25; @@ -266,7 +266,7 @@ private: for (unsigned int k = 1; k <= maxIter; k++) { sw = 0.5*(upper + lower); - pcgw = MaterialLaw::pcgw(pcParams, sw); + pcgw = fluidMatrixInteraction.pcgw(sw, 0.0/*sn*/); const Scalar delta = std::abs(pcgw - pcIn); if (delta < bisLimit) return sw; diff --git a/lecture/mm/naplinfiltration/spatialparams.hh b/lecture/mm/naplinfiltration/spatialparams.hh index ac7326f..794a701 100644 --- a/lecture/mm/naplinfiltration/spatialparams.hh +++ b/lecture/mm/naplinfiltration/spatialparams.hh @@ -26,11 +26,10 @@ #include #include -#include -#include -#include +#include +#include #include -#include +#include namespace Dumux { /*! @@ -47,19 +46,14 @@ class InfiltrationSpatialParams using Element = typename GridView::template Codim<0>::Entity; using ParentType = FVSpatialParams>; - using GlobalPosition = typename SubControlVolume::GlobalPosition; - - using EffectiveLaw = RegularizedParkerVanGen3P; + using ThreePhasePcKrSw = FluidMatrix::ParkerVanGenuchten3PDefault; + using AdsorptionModel = FluidMatrix::ThreePNAPLAdsorption; public: // export permeability type using PermeabilityType = Scalar; - //get the material law from the property system - using MaterialLaw = EffToAbsLaw; - using MaterialLawParams = typename MaterialLaw::Params; - /*! * \brief The constructor * @@ -67,32 +61,14 @@ public: */ InfiltrationSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) + , threePhasePcKrSw_("SpatialParams") + , adsorption_("SpatialParams") { // intrinsic permeabilities permeability_ = getParam("SpatialParams.permeability"); // porosities porosity_ = getParam("SpatialParams.porosity"); - vGAlpha_ = getParam("SpatialParams.vanGenuchtenAlpha"); - vGN_ = getParam("SpatialParams.vanGenuchtenN"); - // residual saturations - materialParams_.setSwr(0.12); - materialParams_.setSnr(0.07); - materialParams_.setSgr(0.03); - - // parameters for the scaling of capillary pressure (GW = 1); - materialParams_.setBetaNw(1.83); - materialParams_.setBetaGn(2.2); - materialParams_.setBetaGw(1.0); - - // parameters for the 3phase van Genuchten law - materialParams_.setVgAlpha(vGAlpha_); - materialParams_.setVgn(vGN_); - materialParams_.setKrRegardsSnr(false); - - // parameters for adsorption - materialParams_.setKdNAPL(0.); - materialParams_.setRhoBulk(1500.); plotFluidMatrixInteractions_ = getParam("Output.PlotFluidMatrixInteractions"); } @@ -105,15 +81,41 @@ public: { GnuplotInterface gnuplot(plotFluidMatrixInteractions_); gnuplot.setOpenPlotWindow(plotFluidMatrixInteractions_); - PlotMaterialLaw plotMaterialLaw(plotFluidMatrixInteractions_); + + const Scalar sg = 0.2; // assume a fixed gas saturation + auto swRange = linspace(sg, 1.0, 1000); + + // assume fixed sn = 0.2 for pcgw curve + auto pcgw = swRange; + std::transform(swRange.begin(), swRange.end(), pcgw.begin(), [&](auto x){ return threePhasePcKrSw_.pcgw(x, sg); }); gnuplot.resetAll(); - plotMaterialLaw.addpc(gnuplot, materialParams_); - gnuplot.plot("pc"); + gnuplot.setXlabel("Sw"); + gnuplot.setYlabel("pcgw"); + gnuplot.addDataSetToPlot(swRange, pcgw, "pcgw", "w l"); + gnuplot.plot("pcgw-sw"); + + // plot kr + swRange = linspace(sg, 0.8, 1000); + auto krw = swRange; + auto krn = swRange; + auto krg = swRange; + for (const auto& [i, sw] : enumerate(swRange)) + { + const Scalar sn = 1.0 - sg - sw; + krw[i] = threePhasePcKrSw_.krw(sw, sn); + krn[i] = threePhasePcKrSw_.krn(sw, sn); + krg[i] = threePhasePcKrSw_.krg(sw, sn); + } gnuplot.resetAll(); - plotMaterialLaw.addkr(gnuplot, materialParams_); - gnuplot.plot("kr"); + gnuplot.setXlabel("Sw"); + gnuplot.setYlabel("kr"); + gnuplot.addDataSetToPlot(swRange, krw, "krw", "w l"); + gnuplot.addDataSetToPlot(swRange, krn, "krn", "w l"); + gnuplot.addDataSetToPlot(swRange, krg, "krg", "w l"); + gnuplot.plot("kr-sw"); + } /*! @@ -128,9 +130,7 @@ public: PermeabilityType permeability(const Element& element, const SubControlVolume& scv, const ElementSolution& elemSol) const - { - return permeability_; - } + { return permeability_; } /*! * \brief Returns the porosity \f$[-]\f$ @@ -138,19 +138,15 @@ public: * \param globalPos The global position */ Scalar porosityAtPos(const GlobalPosition& globalPos) const - { - return porosity_; - } + { return porosity_; } /*! * \brief Returns the parameter object for the Brooks-Corey material law * * \param globalPos The global position */ - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const - { - return materialParams_; - } + auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const + { return makeFluidMatrixInteraction(threePhasePcKrSw_, adsorption_); } private: /* bool isFineMaterial_(const GlobalPosition &globalPos) const @@ -161,12 +157,10 @@ private: */ Scalar permeability_; - Scalar porosity_; - Scalar vGN_; - Scalar vGAlpha_; - MaterialLawParams materialParams_; + const ThreePhasePcKrSw threePhasePcKrSw_; + const AdsorptionModel adsorption_; bool plotFluidMatrixInteractions_; -- GitLab From 594ef5a103e82701e217aad70a9852cc10aac00f Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Tue, 3 Nov 2020 11:54:38 +0100 Subject: [PATCH 10/15] [mm][henryproblem] adapt to new material laws --- .../henry1p2c/henry1p2cspatialparameters.hh | 9 +--- .../henry2p/henry2pspatialparams.hh | 45 ++++++------------- 2 files changed, 15 insertions(+), 39 deletions(-) diff --git a/lecture/mm/henryproblem/henry1p2c/henry1p2cspatialparameters.hh b/lecture/mm/henryproblem/henry1p2c/henry1p2cspatialparameters.hh index 6f9e74a..daa1efd 100644 --- a/lecture/mm/henryproblem/henry1p2c/henry1p2cspatialparameters.hh +++ b/lecture/mm/henryproblem/henry1p2c/henry1p2cspatialparameters.hh @@ -27,7 +27,6 @@ #include #include -#include namespace Dumux { @@ -68,9 +67,7 @@ public: * \param globalPos The global position */ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const - { - return perm_; - } + { return perm_; } /*! * \brief Define the porosity \f$\mathrm{[-]}\f$. @@ -78,9 +75,7 @@ public: * \param globalPos The global position */ double porosityAtPos(const GlobalPosition& globalPos) const - { - return porosity_; - } + { return porosity_; } private: diff --git a/lecture/mm/henryproblem/henry2p/henry2pspatialparams.hh b/lecture/mm/henryproblem/henry2p/henry2pspatialparams.hh index 0530541..f55a5a5 100644 --- a/lecture/mm/henryproblem/henry2p/henry2pspatialparams.hh +++ b/lecture/mm/henryproblem/henry2p/henry2pspatialparams.hh @@ -27,8 +27,6 @@ #include #include -#include - #include namespace Dumux { @@ -51,12 +49,10 @@ public FVSpatialParams; + public: using PermeabilityType = Scalar; - //get the material law from the property system - using EffectiveLaw = LinearMaterial; - using MaterialLaw = EffToAbsLaw; - using MaterialLawParams = typename MaterialLaw::Params; Henry2pSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) @@ -72,15 +68,8 @@ public: std::cerr << "Unknown exception thrown!\n"; exit(1); } - - // residual saturations - materialParams_.setSwr(0.0); - materialParams_.setSnr(0.0); - - // parameters for the linear law - // entry pc and max pc - materialParams_.setEntryPc(0); - materialParams_.setMaxPc(0); + typename PcKrSwCurve::BasicParams params(0.0/*pcEntry*/, 0.0/*pcMax*/); + pcKrSwCurve_ = std::make_unique(params); K_ = 1.019368e-9; porosity_=0.35; @@ -92,9 +81,7 @@ public: * \param globalPos The global position */ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const - { - return K_; - } + { return K_; } /*! * \brief Define the porosity \f$\mathrm{[-]}\f$. @@ -102,22 +89,17 @@ public: * \param globalPos The global position */ double porosityAtPos(const GlobalPosition& globalPos) const - { - return porosity_; - } + { return porosity_; } /*! * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.). * - * \param element The current element - * \param fvElemGeom The current finite volume geometry of the element - * \param scvIdx The index of the sub-control volume. + * \param globalPos The global position * \return the material parameters object */ - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const - { - return materialParams_; - } + auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const + { return makeFluidMatrixInteraction(*pcKrSwCurve_); } + /*! * \brief Function for defining which phase is to be considered as the wetting phase. @@ -127,14 +109,13 @@ public: */ template int wettingPhaseAtPos(const GlobalPosition& globalPos) const - { - return FluidSystem::phase0Idx; - } + { return FluidSystem::phase0Idx; } private: Scalar K_; Scalar porosity_; - MaterialLawParams materialParams_; + std::unique_ptr pcKrSwCurve_; + }; -- GitLab From 2d76ea3e366a66e3192b56d8559c1b39291a2b55 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Tue, 3 Nov 2020 13:33:50 +0100 Subject: [PATCH 11/15] [mm][heatpipe] adapt to new material laws --- lecture/mm/heatpipe/heatpipespatialparams.hh | 31 +++---- lecture/mm/heatpipe/krpcheatpipe.hh | 75 +++++++++++++--- lecture/mm/heatpipe/krpcheatpipeparams.hh | 93 -------------------- 3 files changed, 74 insertions(+), 125 deletions(-) delete mode 100644 lecture/mm/heatpipe/krpcheatpipeparams.hh diff --git a/lecture/mm/heatpipe/heatpipespatialparams.hh b/lecture/mm/heatpipe/heatpipespatialparams.hh index 65d9daf..e05d26d 100644 --- a/lecture/mm/heatpipe/heatpipespatialparams.hh +++ b/lecture/mm/heatpipe/heatpipespatialparams.hh @@ -38,47 +38,40 @@ class HeatPipeSpatialParams using Element = typename GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + using PcKrSwCurve = FluidMatrix::KrPcHeatPipeDefault; + public: using PermeabilityType = Scalar; - using MaterialLaw = KrPcHeatpipe; - using MaterialLawParams = typename MaterialLaw::Params; HeatPipeSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) { permeability_ = getParam("Problem.Permeability"); porosity_ = 0.4; + Scalar p0 = std::pow((porosity_/permeability_), 0.5); - materialParams_.setSwr(0.15); - materialParams_.setSnr(0.0); - materialParams_.setP0(std::pow((porosity_/permeability_), 0.5)); + typename PcKrSwCurve::BasicParams params(0.15, 0.0, p0); + pcKrSwCurve_ = std::make_unique(params); } PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const - { - return permeability_; - } + { return permeability_; } Scalar porosityAtPos(const GlobalPosition& globalPos) const - { - return porosity_; - } + { return porosity_; } + + auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const + { return makeFluidMatrixInteraction(*pcKrSwCurve_); } - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const - { - return materialParams_; - } template int wettingPhaseAtPos(const GlobalPosition& globalPos) const - { - return FluidSystem::phase0Idx; - } + { return FluidSystem::phase0Idx; } private: PermeabilityType permeability_; Scalar porosity_; - MaterialLawParams materialParams_; + std::unique_ptr pcKrSwCurve_; }; } diff --git a/lecture/mm/heatpipe/krpcheatpipe.hh b/lecture/mm/heatpipe/krpcheatpipe.hh index 0183150..67b5397 100644 --- a/lecture/mm/heatpipe/krpcheatpipe.hh +++ b/lecture/mm/heatpipe/krpcheatpipe.hh @@ -19,13 +19,15 @@ #ifndef DUMUX_KR_PC_HEATPIPE_HH #define DUMUX_KR_PC_HEATPIPE_HH -#include -#include "krpcheatpipeparams.hh" - #include +#include -namespace Dumux -{ +#include +#include +#include +#include + +namespace Dumux::FluidMatrix { /*! * * \brief Implementation of the capillary pressure <-> saturation @@ -34,18 +36,57 @@ namespace Dumux * cap-press based on the function of Leverett. * */ -template > -class KrPcHeatpipe +class KrPcHeatPipe { public: - using Params = ParamsT; - using Scalar = ScalarT; + template + struct Params + { + Params(Scalar swr, Scalar snr, Scalar p0) + : swr_(swr) + , snr_(snr) + , p0_(p0) + {} + + Scalar swr() const{ return swr_; } + void setSwr(Scalar swr){ swr_ = swr; } + + Scalar snr() const { return snr_; } + void setSnr(Scalar snr) { snr_ = snr; } + + Scalar p0() const { return p0_; } + void setp0(Scalar p0) { p0_ = p0; } + + bool operator== (const Params& p) const + { + return Dune::FloatCmp::eq(snr(), p.snr(), 1e-6) + && Dune::FloatCmp::eq(snr(), p.snr(), 1e-6) + && Dune::FloatCmp::eq(p0(), p.p0(), 1e-6); + } + + private: + Scalar swr_, snr_, p0_; + }; + + /*! + * \brief Construct from a subgroup from the global parameter tree + * \note This will give you nice error messages if a mandatory parameter is missing + */ + template + static Params makeParams(const std::string& paramGroup) + { + const auto swr = getParamFromGroup(paramGroup, "Swr"); + const auto snr = getParamFromGroup(paramGroup, "Snr"); + const auto p0 = getParamFromGroup(paramGroup, "P0"); + return {swr, snr, p0}; + } /*! * \brief The capillary pressure-saturation curve according to Leverett. * */ - static Scalar pc(const Params ¶ms, Scalar sw) + template + static Scalar pc(Scalar sw, const Params& params) { if(sw<0.) sw=0.; /* effective values */ @@ -72,7 +113,8 @@ public: * the medium according to the Fatt-Klikoff * parameterization. */ - static Scalar krw(const Params ¶ms, Scalar sw) + template + static Scalar krw(Scalar sw, const Params& params) { /*** effective saturation ***/ Scalar Se = (sw-params.swr())/(1-params.snr()-params.swr()); @@ -88,7 +130,8 @@ public: * the medium according to the Fatt-Klikoff * parameterization. */ - static Scalar krn(const Params ¶ms, Scalar sw) + template + static Scalar krn(Scalar sw, const Params& params) { /*** effective saturation ***/ Scalar Se = (sw-params.swr())/(1-params.snr()-params.swr()); @@ -98,9 +141,15 @@ public: /* compute and return value */ return std::pow(1.-Se,3); }; +}; +/*! + * \ingroup Fluidmatrixinteractions + * \brief A default configuration for using the VanGenuchten material law + */ +template +using KrPcHeatPipeDefault = TwoPMaterialLaw; -}; } #endif diff --git a/lecture/mm/heatpipe/krpcheatpipeparams.hh b/lecture/mm/heatpipe/krpcheatpipeparams.hh deleted file mode 100644 index d0d9bec..0000000 --- a/lecture/mm/heatpipe/krpcheatpipeparams.hh +++ /dev/null @@ -1,93 +0,0 @@ -// -*- 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 2 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 . * - *****************************************************************************/ - -/*! - * \file - * - * \brief Specification of the material parameters - * for the Brooks Corey constitutive relations. - */ -#ifndef DUMUX_KR_PC_HEATPIPE_PARAMS_HH -#define DUMUX_KR_PC_HEATPIPE_PARAMS_HH - -namespace Dumux { - -/*! - * \brief Specification of the material parameters - * for the kr-pc constitutive relations for the heatpipe problem. - * - */ -template -class KrPcHeatpipeParams -{ -public: - typedef ScalarT Scalar; - - KrPcHeatpipeParams() = default; - - KrPcHeatpipeParams(Scalar swr, Scalar snr) - : swr_(swr), snr_(snr) - { - } - - /*! - * \brief Returns the entry pressure [Pa] - */ - Scalar swr() const - { return swr_; } - - /*! - * \brief Set the entry pressure [Pa] - */ - void setSwr(Scalar v) - { swr_ = v; } - - - /*! - * \brief Returns the lambda shape parameter - */ - Scalar snr() const - { return snr_; } - - /*! - * \brief Set the lambda shape parameter - */ - void setSnr(Scalar v) - { snr_ = v; } - - /*! - * \brief Returns the lambda shape parameter - */ - Scalar p0() const - { return p0_; } - - /*! - * \brief Set the lambda shape parameter - */ - void setP0(Scalar v) - { p0_ = v; } - -private: - Scalar swr_; - Scalar snr_; - Scalar p0_; -}; -} // namespace Dumux - -#endif -- GitLab From d587cb1ee994499b4c91484fa2fb45557b4d9e01 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Mon, 2 Nov 2020 23:43:36 +0100 Subject: [PATCH 12/15] [mm][buckleyLeverett] adapt to new material laws --- .../buckleyleverettanalyticsolution.hh | 76 +++++-------------- .../buckleyleverettexercise.cc | 6 +- .../buckleyleverettexercise.input | 9 +-- .../buckleyleverett/buckleyleverettproblem.hh | 4 +- .../buckleyleverettspatialparams.hh | 73 +++++++----------- 5 files changed, 58 insertions(+), 110 deletions(-) diff --git a/lecture/mm/buckleyleverett/buckleyleverettanalyticsolution.hh b/lecture/mm/buckleyleverett/buckleyleverettanalyticsolution.hh index 012b7de..5dec1ec 100644 --- a/lecture/mm/buckleyleverett/buckleyleverettanalyticsolution.hh +++ b/lecture/mm/buckleyleverett/buckleyleverettanalyticsolution.hh @@ -21,8 +21,7 @@ #include #include -#include -#include +#include /*! * \file @@ -36,33 +35,6 @@ namespace Dumux { * the Buckley-Leverett problem */ -template -struct CheckMaterialLaw -{ - static bool isLinear() - { - return false; - } -}; - -template -struct CheckMaterialLaw > -{ - static bool isLinear() - { - return true; - } -}; - -template -struct CheckMaterialLaw > > -{ - static bool isLinear() - { - return true; - } -}; - template class BuckleyLeverettAnalytic { using Problem = GetPropType; @@ -70,8 +42,6 @@ template class BuckleyLeverettAnalytic using GridView = typename GetPropType::GridView; using Scalar = GetPropType; using SpatialParams = GetPropType; - using MaterialLaw = typename SpatialParams::MaterialLaw; - using MaterialLawParams = typename MaterialLaw::Params; using FluidSystem = GetPropType; using FluidState = GetPropType; using Indices = GetPropType; @@ -87,7 +57,6 @@ template class BuckleyLeverettAnalytic using GlobalPosition = Dune::FieldVector; public: - // functions needed for analytical solution void initializeAnalytic() @@ -106,11 +75,11 @@ public: */ void prepareAnalytic() { - const MaterialLawParams& materialLawParams(problem_.spatialParams().materialLawParams(dummyElement_)); - - swr_ = materialLawParams.swr(); - snr_ = materialLawParams.snr(); - porosity_ = problem_.spatialParams().porosity(dummyElement_); + const auto& dummyElement = *problem_.gridView().template begin<0>(); + const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(dummyElement.geometry().center()); + swr_ = fluidMatrixInteraction.effToAbsParams().swr(); + snr_ = fluidMatrixInteraction.effToAbsParams().snr(); + porosity_ = problem_.spatialParams().porosity(dummyElement); time_ = 0; satVec_ = swr_; @@ -128,12 +97,11 @@ public: for (int i = 0; i < pointNum_; i++) { - fractionalW_[i] = MaterialLaw::krw(materialLawParams, satVec_[i])/viscosityW; - fractionalW_[i] /= (fractionalW_[i] + MaterialLaw::krn(materialLawParams, satVec_[i])/viscosityNW); + fractionalW_[i] = fluidMatrixInteraction.krw(satVec_[i])/viscosityW; + fractionalW_[i] /= (fractionalW_[i] + fluidMatrixInteraction.krn(satVec_[i])/viscosityNW); } dfwdsw_ = 0; - for (int i = 1; i < intervalNum_; i++) { dfwdsw_[i] = (fractionalW_[i + 1] - fractionalW_[i - 1]) / (satVec_[i + 1] - satVec_[i - 1]); @@ -177,12 +145,13 @@ public: { xf_[i] = vTot_ * time_ / porosity_ * dfwdsw_[i]; } + + // position of maximum xf_ + int xfmax = 0; int xhelp = pointNum_ / 3; int xhelpold = 0; int xhelpoldold = 0; - int xfmax = 0; - - // position of maximum xf_ + int xhelp2 = 0; for (int i = 0; i < pointNum_; i++) { if (xf_[i] > xf_[i + 1]) @@ -197,13 +166,9 @@ public: Scalar A1; Scalar A2; Scalar b; - int xhelp2 = 0; while (a) { - if (CheckMaterialLaw::isLinear()) - break; - A1 = 0; for (int i = 0; i <= xhelp - 1; i++) @@ -262,7 +227,7 @@ public: int index = problem_.variables().index(*eIt); // account for linear material law - if (CheckMaterialLaw::isLinear()) + if constexpr (SpatialParams::pcSwCurveIsLinear()) { if (globalPos[0] <= xf_[1]) { @@ -313,9 +278,7 @@ public: } BlockVector AnalyticSolution() const - { - return analyticSolution_; - } + { return analyticSolution_; } //Write saturation and pressure into file template @@ -332,9 +295,13 @@ public: } //! Construct an IMPES object. - BuckleyLeverettAnalytic(Problem& problem, Scalar totalVelocity = 3e-7) : - problem_(problem), analyticSolution_(0), error_(0), elementVolume_(0), size_(problem.gridView().size(0)), vTot_(totalVelocity), dummyElement_( - *(problem_.gridView().template begin<0> ())) + BuckleyLeverettAnalytic(Problem& problem, Scalar totalVelocity = 3e-7) + : problem_(problem) + , analyticSolution_(0) + , error_(0) + , elementVolume_(0) + , size_(problem.gridView().size(0)) + , vTot_(totalVelocity) { dummyGlobal_ = 0.0; dummyGlobal_[0] = 1.0; @@ -365,7 +332,6 @@ private: Dune::FieldVector dfwdsw_; Dune::FieldVector xf_; int dfwdswmax_; - const Element& dummyElement_; GlobalPosition dummyGlobal_; }; diff --git a/lecture/mm/buckleyleverett/buckleyleverettexercise.cc b/lecture/mm/buckleyleverett/buckleyleverettexercise.cc index 6673a30..7688e91 100644 --- a/lecture/mm/buckleyleverett/buckleyleverettexercise.cc +++ b/lecture/mm/buckleyleverett/buckleyleverettexercise.cc @@ -48,10 +48,10 @@ void usage(const char *progName, const std::string &errorMsg) "\t-SpatialParams.Porosity The porosity of the porous medium [-]\n" "\t-SpatialParams.BrooksCoreyLambda The pore size distribution parameter for the \n" "\t \t Brooks-Corey capillary pressure - saturation relationship [-]\n" - "\t-SpatialParams.BrooksCoreyEntryPressure The entry pressure for the \n" + "\t-SpatialParams.BrooksCoreyPcEntry The entry pressure for the \n" "\t \t Brooks-Corey capillary pressure - saturation relationship [Pa]\n" - "\t-SpatialParams.ResidualSaturationWetting The residual saturation of the wetting phase [-]\n" - "\t-SpatialParams.ResidualSaturationNonwetting The residual saturation of the nonwetting phase [-]\n" + "\t-SpatialParams.Swr The residual saturation of the wetting phase [-]\n" + "\t-SpatialParams.Snr The residual saturation of the nonwetting phase [-]\n" "\t-Fluid.DensityW The density of the wetting phase [kg/m^3]\n" "\t-Fluid.DensityNW The density of the nonwetting phase [kg/m^3]\n" "\t-Fluid.ViscosityW The dynamic viscosity of the wetting phase [kg/(ms)]\n" diff --git a/lecture/mm/buckleyleverett/buckleyleverettexercise.input b/lecture/mm/buckleyleverett/buckleyleverettexercise.input index b5c3e9f..b4defe1 100644 --- a/lecture/mm/buckleyleverett/buckleyleverettexercise.input +++ b/lecture/mm/buckleyleverett/buckleyleverettexercise.input @@ -6,15 +6,12 @@ DtInitial = 1e3 # initial time step for the simulati EnableGravity = false [SpatialParams] - Permeability = 1.01936799e-14 # intrinsic permeability of the porous medium [m^2] Porosity = 0.2 # porosity of the porous medium [-] - BrooksCoreyLambda = 4.0 # pore size distribution parameter for the Brooks-Corey capillary pressure - saturation relationship [-] -BrooksCoreyEntryPressure = 0 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship [Pa] - -ResidualSaturationWetting = 0.2 # residual saturation of the wetting phase [-] -ResidualSaturationNonwetting = 0.2 # residual saturation of the nonwetting phase [-] +BrooksCoreyPcEntry = 0 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship [Pa] +Swr = 0.2 # residual saturation of the wetting phase [-] +Snr = 0.2 # residual saturation of the nonwetting phase [-] [Fluid] DensityW = 1e3 # density of the wetting phase [kg/m^3] diff --git a/lecture/mm/buckleyleverett/buckleyleverettproblem.hh b/lecture/mm/buckleyleverett/buckleyleverettproblem.hh index 9a43ca8..b7eac91 100644 --- a/lecture/mm/buckleyleverett/buckleyleverettproblem.hh +++ b/lecture/mm/buckleyleverett/buckleyleverettproblem.hh @@ -90,8 +90,8 @@ public: densityNonwetting_ = getParam("Fluid.DensityNW"); - swr_ = getParam("SpatialParams.ResidualSaturationWetting"); - snr_ = getParam("SpatialParams.ResidualSaturationNonwetting"); + swr_ = getParam("SpatialParams.Swr"); + snr_ = getParam("SpatialParams.Snr"); paraviewOutput_ = getParam("Output.paraviewOutput", true); } diff --git a/lecture/mm/buckleyleverett/buckleyleverettspatialparams.hh b/lecture/mm/buckleyleverett/buckleyleverettspatialparams.hh index 9c3f90f..6bc65d1 100644 --- a/lecture/mm/buckleyleverett/buckleyleverettspatialparams.hh +++ b/lecture/mm/buckleyleverett/buckleyleverettspatialparams.hh @@ -23,9 +23,7 @@ #include #include -#include -#include -#include +#include namespace Dumux { @@ -44,19 +42,12 @@ struct BuckleyLeverettSpatialParamsTypeTag {}; template struct SpatialParams { using type = BuckleyLeverettSpatialParams; }; -// Set the material law -template -struct MaterialLaw -{ -private: - using Scalar = GetPropType; - using RawMaterialLaw = RegularizedBrooksCorey; -public: - using type = EffToAbsLaw; -}; - } // end namespace Properties +// forward declaration +class LinearMaterialDefault; +class LinearMaterial; + template class BuckleyLeverettSpatialParams: public SequentialFVSpatialParams { @@ -71,50 +62,44 @@ class BuckleyLeverettSpatialParams: public SequentialFVSpatialParams using Element = typename Grid::Traits::template Codim<0>::Entity; using GlobalPosition = Dune::FieldVector; using FieldMatrix = Dune::FieldMatrix; + using PcKrSwCurve = FluidMatrix::BrooksCoreyDefault; public: - using MaterialLaw = GetPropType; - using MaterialLawParams = typename MaterialLaw::Params; - Scalar intrinsicPermeability (const Element& element) const + BuckleyLeverettSpatialParams(const Problem& problem) + : ParentType(problem) + , pcKrSwCurve_("SpatialParams") { - return constPermeability_; + Scalar permFactor = 1.0; //0.001/(1000*9.81); + constPermeability_ = getParam("SpatialParams.Permeability")*permFactor; + porosity_ = getParam("SpatialParams.Porosity"); } + static constexpr bool pcSwCurveIsLinear() + { return (std::is_same_v || std::is_same_v); } + + Scalar intrinsicPermeability (const Element& element) const + { return constPermeability_; } + Scalar porosity(const Element &element) const - { - return porosity_; - } + { return porosity_; } /*! - * \brief DOC ME! - * \param element DOC ME! + * \brief Returns the parameters for the material law at a given location + * + * \param globalPos The global coordinates for the given location */ - - // return the parameter object for the Brooks-Corey material law which depends on the position - const MaterialLawParams& materialLawParams(const Element &element) const - { - return materialLawParams_; - } - - BuckleyLeverettSpatialParams(const Problem& problem) - :ParentType(problem) - { - Scalar permFactor = 1.0; //0.001/(1000*9.81); - - constPermeability_ = getParam("SpatialParams.Permeability")*permFactor; - materialLawParams_.setSwr( getParam("SpatialParams.ResidualSaturationWetting") ); - materialLawParams_.setSnr( getParam("SpatialParams.ResidualSaturationNonwetting") ); - //set Brooks-Corey parameters - materialLawParams_.setPe( getParam("SpatialParams.BrooksCoreyEntryPressure") ); - materialLawParams_.setLambda( getParam("SpatialParams.BrooksCoreyLambda") ); - porosity_ = getParam("SpatialParams.Porosity"); - } + auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const + { return pcKrSwCurve_; } private: - MaterialLawParams materialLawParams_; + const PcKrSwCurve pcKrSwCurve_; Scalar constPermeability_; Scalar porosity_; + Scalar swr_; + Scalar snr_; + Scalar pcEntry_; + Scalar lambda_; }; -- GitLab From 59925954a9b8d2c8b7635a2b20ac852a68f1aa27 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Mon, 2 Nov 2020 23:43:48 +0100 Subject: [PATCH 13/15] [mm][mcwhorter] adapt to new material laws --- .../mm/mcwhorter/mcwhorteranalyticsolution.hh | 33 ++++++++++--------- lecture/mm/mcwhorter/mcwhorterexercise.cc | 6 ++-- lecture/mm/mcwhorter/mcwhorterexercise.input | 6 ++-- lecture/mm/mcwhorter/mcwhorterproblem.hh | 4 +-- 4 files changed, 25 insertions(+), 24 deletions(-) diff --git a/lecture/mm/mcwhorter/mcwhorteranalyticsolution.hh b/lecture/mm/mcwhorter/mcwhorteranalyticsolution.hh index ecb7600..6166e7c 100644 --- a/lecture/mm/mcwhorter/mcwhorteranalyticsolution.hh +++ b/lecture/mm/mcwhorter/mcwhorteranalyticsolution.hh @@ -60,8 +60,6 @@ class McWhorterAnalytic using ElementIterator = typename GridView::template Codim<0>::Iterator; public: - using MaterialLaw = typename SpatialParams::MaterialLaw; - using MaterialLawParams = typename MaterialLaw::Params; // functions needed for analytical solution void initializeAnalytic() @@ -97,11 +95,12 @@ public: void prepareAnalytic() { - const MaterialLawParams& materialLawParams( problem_.spatialParams().materialLawParams(dummyElement_) ); - swr_ = materialLawParams.swr(); - snr_ = materialLawParams.snr(); - porosity_ = problem_.spatialParams().porosity(dummyElement_); - permeability_ = problem_.spatialParams().intrinsicPermeability(dummyElement_); + const auto& dummyElement = *problem_.gridView().template begin<0>(); + const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(dummyElement.geometry().center()); + swr_ = fluidMatrixInteraction.effToAbsParams().swr(); + snr_ = fluidMatrixInteraction.effToAbsParams().snr(); + porosity_ = problem_.spatialParams().porosity(dummyElement); + permeability_ = problem_.spatialParams().intrinsicPermeability(dummyElement); PrimaryVariables initVec; problem_.initialAtPos(initVec, dummyGlobal_); sInit_ = initVec[saturationIdx]; @@ -126,8 +125,8 @@ public: // get fractional flow function vector for (int i=0; i ())), tolAnalytic_(tolAnalytic) + McWhorterAnalytic(Problem& problem, Scalar tolAnalytic = 1e-14) + : problem_(problem) + , analyticSolution_(0) + , error_(0) + , elementVolume_(0) + , size_(problem.gridView().size(0)) + , tolAnalytic_(tolAnalytic) { dummyGlobal_ = 0.0; dummyGlobal_[0] = 1.0; @@ -318,8 +321,6 @@ private: BlockVector elementVolume_; int size_; - - const Element& dummyElement_; GlobalPosition dummyGlobal_; Scalar tolAnalytic_; diff --git a/lecture/mm/mcwhorter/mcwhorterexercise.cc b/lecture/mm/mcwhorter/mcwhorterexercise.cc index 9406595..8996631 100644 --- a/lecture/mm/mcwhorter/mcwhorterexercise.cc +++ b/lecture/mm/mcwhorter/mcwhorterexercise.cc @@ -48,10 +48,10 @@ void usage(const char *progName, const std::string &errorMsg) "\t-SpatialParams.Porosity The porosity of the porous medium [-]\n" "\t-SpatialParams.BrooksCoreyLambda The pore size distribution parameter for the \n" "\t \t Brooks-Corey capillary pressure - saturation relationship [-]\n" - "\t-SpatialParams.BrooksCoreyEntryPressure The entry pressure for the \n" + "\t-SpatialParams.BrooksCoreyPcEntry The entry pressure for the \n" "\t \t Brooks-Corey capillary pressure - saturation relationship [Pa]\n" - "\t-SpatialParams.ResidualSaturationWetting The residual saturation of the wetting phase [-]\n" - "\t-SpatialParams.ResidualSaturationNonwetting The residual saturation of the nonwetting phase [-]\n" + "\t-SpatialParams.Swr The residual saturation of the wetting phase [-]\n" + "\t-SpatialParams.Snr The residual saturation of the nonwetting phase [-]\n" "\t-Fluid.DensityW The density of the wetting phase [kg/m^3]\n" "\t-Fluid.DensityNW The density of the nonwetting phase [kg/m^3]\n" "\t-Fluid.ViscosityW The dynamic viscosity of the wetting phase [kg/(ms)]\n" diff --git a/lecture/mm/mcwhorter/mcwhorterexercise.input b/lecture/mm/mcwhorter/mcwhorterexercise.input index 80e6ba5..486b3e8 100644 --- a/lecture/mm/mcwhorter/mcwhorterexercise.input +++ b/lecture/mm/mcwhorter/mcwhorterexercise.input @@ -9,9 +9,9 @@ EnableGravity = false Permeability = 1.01936799e-13 # intrinsic permeability of the porous medium [m^2] Porosity = 0.2 # porosity of the porous medium [-] BrooksCoreyLambda = 3.0 # pore size distribution parameter for the Brooks-Corey capillary pressure - saturation relationship [-] -BrooksCoreyEntryPressure = 8000 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship [Pa] -ResidualSaturationWetting = 0.2 # residual saturation of the wetting phase [-] -ResidualSaturationNonwetting = 0.2 # residual saturation of the nonwetting phase [-] +BrooksCoreyPcEntry = 8000 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship [Pa] +Swr = 0.2 # residual saturation of the wetting phase [-] +Snr = 0.2 # residual saturation of the nonwetting phase [-] [Fluid] DensityW = 1e3 # density of the wetting phase [kg/m^3] diff --git a/lecture/mm/mcwhorter/mcwhorterproblem.hh b/lecture/mm/mcwhorter/mcwhorterproblem.hh index 34af635..ed5975c 100644 --- a/lecture/mm/mcwhorter/mcwhorterproblem.hh +++ b/lecture/mm/mcwhorter/mcwhorterproblem.hh @@ -79,8 +79,8 @@ public: PseudoOil::setDensity( getParam("Fluid.DensityW") ); PseudoH2O::setDensity( getParam("Fluid.DensityNW") ); - swr_ = getParam("SpatialParams.ResidualSaturationWetting"); - snr_ = getParam("SpatialParams.ResidualSaturationNonwetting"); + swr_ = getParam("SpatialParams.Swr"); + snr_ = getParam("SpatialParams.Snr"); paraviewOutput_ = getParam("Output.paraviewOutput", true); } -- GitLab From 766cb8d22880ac5e05b4caf8c635e46db763f398 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Mon, 2 Nov 2020 16:30:55 +0100 Subject: [PATCH 14/15] [mm][fuelcell] Adapt to new material laws --- lecture/mm/fuelcell/fuelcell.input | 3 + lecture/mm/fuelcell/fuelcellspatialparams.hh | 53 +-- lecture/mm/fuelcell/material/acosta.hh | 319 ++++++++++++++++-- lecture/mm/fuelcell/material/acostaparams.hh | 113 ------- .../mm/fuelcell/material/regularizedacosta.hh | 284 ---------------- .../material/regularizedacostaparams.hh | 76 ----- 6 files changed, 303 insertions(+), 545 deletions(-) delete mode 100644 lecture/mm/fuelcell/material/acostaparams.hh delete mode 100644 lecture/mm/fuelcell/material/regularizedacosta.hh delete mode 100644 lecture/mm/fuelcell/material/regularizedacostaparams.hh diff --git a/lecture/mm/fuelcell/fuelcell.input b/lecture/mm/fuelcell/fuelcell.input index 21258b5..669178c 100644 --- a/lecture/mm/fuelcell/fuelcell.input +++ b/lecture/mm/fuelcell/fuelcell.input @@ -56,3 +56,6 @@ SolidDensity = 1430 # [kg/m^3] SolidThermalConductivity = 15.6 # [W/(m*K)] SolidHeatCapacity = 710 # [J/(kg*K)] +[SpatialParams] +Swr = 0.05 +Snr = 0.05 diff --git a/lecture/mm/fuelcell/fuelcellspatialparams.hh b/lecture/mm/fuelcell/fuelcellspatialparams.hh index 930c3fe..a454a1f 100644 --- a/lecture/mm/fuelcell/fuelcellspatialparams.hh +++ b/lecture/mm/fuelcell/fuelcellspatialparams.hh @@ -28,12 +28,7 @@ #include #include -#include -#include -#include -#include - -#include "./material/regularizedacosta.hh" +#include "./material/acosta.hh" #include "./material/thermalconductivityconstant.hh" namespace Dumux { @@ -80,11 +75,8 @@ class FuelCellLectureSpatialParams: public FVSpatialParams; - + using PcKrSwCurve = FluidMatrix::AcostaPcSwDefault; public: - using MaterialLaw = EffToAbsLaw; - using MaterialLawParams = typename MaterialLaw::Params; using PermeabilityType = Scalar; /*! @@ -98,23 +90,14 @@ public: { // intrinsic permeabilities K_ = 5.2e-11; // Acosta absolute permeability diffusion layer - // porosities porosity_ = 0.78; // Acosta porosity diffusion layer - // thermalconductivity lambdaSolid_ = 15.6; // [W/(m*K)] Acosta thermal conductivity used in capillary pressure-saturation - // residual saturations - materialParams_.setSwr(0.05); // value for imbibition from Acosta - materialParams_.setSnr(0.05); // assumption + typename PcKrSwCurve::BasicParams params(-1168.75, 8.5, -0.2, -700, 0.0); + pcKrSwCurve_ = std::make_unique(params); - // parameters for the Acosta saturation capillary pressure relation - materialParams_.setAcA(-1168.75); // imbition -1168.75; drainage -600; - materialParams_.setAcB(8.5); // imbition 8.5; drainage 25; - materialParams_.setAcC(-0.2); // imbition -0.2; drainage -16; - materialParams_.setAcD(-700); // imbition -700; drainage -3300; - materialParams_.setAcE(0.0); // imbition 0.0; drainage 800; } ~FuelCellLectureSpatialParams() @@ -126,9 +109,7 @@ public: * \param globalPos The global position */ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const - { - return K_; - } + { return K_; } /*! * \brief Defines the porosity \f$[-]\f$ of the spatial parameters. @@ -154,19 +135,15 @@ public: */ template int wettingPhaseAtPos(const GlobalPosition& globalPos) const - { - return FluidSystem::phase0Idx; - } + { return FluidSystem::phase0Idx; } /*! * \brief Returns the parameter object for the Brooks-Corey material law which depends on the position. * * \param globalPos The global position */ - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const - { - return materialParams_; - } + auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const + { return makeFluidMatrixInteraction(*pcKrSwCurve_); } /*! * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix. @@ -181,9 +158,7 @@ public: Scalar solidHeatCapacity(const Element &element, const FVElementGeometry &fvGeometry, const int scvIdx) const - { - return 710; // specific heat capacity of diffusion layer Acosta [J / (kg K)] - } + { return 710; }// specific heat capacity of diffusion layer Acosta [J / (kg K)] /*! * \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix. @@ -198,9 +173,7 @@ public: Scalar solidDensity(const Element &element, const FVElementGeometry &fvGeometry, const int scvIdx) const - { - return 1430; // density of ELAT [kg/m^3] Wöhr - } + { return 1430; } // density of ELAT [kg/m^3] Wöhr /*! * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid. @@ -214,15 +187,13 @@ public: Scalar solidThermalConductivity(const Element &element, const FVElementGeometry &fvGeometry, const SubControlVolume& scv) const - { - return lambdaSolid_; - } + { return lambdaSolid_; } private: Scalar K_; Scalar porosity_; Scalar eps_ = 1e-6; - MaterialLawParams materialParams_; + std::unique_ptr pcKrSwCurve_; Scalar lambdaSolid_; }; diff --git a/lecture/mm/fuelcell/material/acosta.hh b/lecture/mm/fuelcell/material/acosta.hh index e5e5650..dc3a5b7 100644 --- a/lecture/mm/fuelcell/material/acosta.hh +++ b/lecture/mm/fuelcell/material/acosta.hh @@ -25,17 +25,15 @@ #ifndef ACOSTA_HH #define ACOSTA_HH -#include "acostaparams.hh" - #include #include -#include -#include -#include +#include +#include +#include +#include -namespace Dumux -{ +namespace Dumux::FluidMatrix { /*! * @@ -48,14 +46,62 @@ namespace Dumux * * \see AcostaParams */ -template > class AcostaPcSw { public: + template + struct Params + { + Params(Scalar acA, Scalar acB, Scalar acC, Scalar acD, Scalar acE) + : acA_(acA) + , acB_(acB) + , acC_(acC) + , acD_(acD) + , acE_(acE) + {} + + Scalar acA() const { return acA_; } + void setAcA(Scalar a) { acA_ = a; } + + Scalar acB() const { return acB_; } + void setAcB(Scalar b) { acB_ = b; } + + Scalar acC() const { return acC_; } + void setAcC(Scalar c) { acC_ = c; } + + Scalar acD() const { return acD_; } + void setAcD(Scalar d) { acD_ = d; } - typedef ParamsT Params; - typedef typename Params::Scalar Scalar; + Scalar acE() const { return acE_; } + void setAcE(Scalar e) { acE_ = e; } + + bool operator== (const Params& p) const + { + return Dune::FloatCmp::eq(acA(), p.acA(), 1e-6) + && Dune::FloatCmp::eq(acB(), p.acB(), 1e-6) + && Dune::FloatCmp::eq(acC(), p.acC(), 1e-6) + && Dune::FloatCmp::eq(acD(), p.acD(), 1e-6) + && Dune::FloatCmp::eq(acE(), p.acE(), 1e-6); + } + + private: + Scalar acA_, acB_, acC_, acD_, acE_; + }; + /*! + * \brief Construct from a subgroup from the global parameter tree + * \note This will give you nice error messages if a mandatory parameter is missing + */ + template + static Params makeParams(const std::string& paramGroup) + { + const auto acA = getParamFromGroup(paramGroup, "AcA"); + const auto acB = getParamFromGroup(paramGroup, "AcB"); + const auto acC = getParamFromGroup(paramGroup, "AcC"); + const auto acD = getParamFromGroup(paramGroup, "AcD"); + const auto acE = getParamFromGroup(paramGroup, "AcE"); + return {acA, acB, acC, acD, acE}; + } /*! * \brief The capillary pressure-saturation curve according to Acosta. * @@ -71,7 +117,8 @@ public: * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too. */ - static Scalar pc(const Params ¶ms, Scalar swe) + template + static Scalar pc(Scalar swe, const Params& params) { assert(0 <= swe && swe <= 1); Scalar pc = exp(params.acB() * swe + params.acC()); @@ -94,9 +141,10 @@ public: * is constructed accordingly. Afterwards the values are set there, too. * \return The effective saturation of the wetting phase */ - static Scalar sw(const Params ¶ms, Scalar pc) + template + static Scalar sw(Scalar pc, const Params& params) { - DUNE_THROW(Dune::NotImplemented, "Acosta::sw(params, pc)"); + DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::sw(params, pc)"); } /*! @@ -111,7 +159,8 @@ public: * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too. */ - static Scalar dpc_dsw(const Params ¶ms, Scalar swe) + template + static Scalar dpc_dsw(Scalar swe, const Params& params) { assert(0 <= swe && swe <= 1); @@ -134,9 +183,10 @@ public: * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too. */ - static Scalar dsw_dpc(const Params ¶ms, Scalar pc) + template + static Scalar dsw_dpc(Scalar pc, const Params& params) { - DUNE_THROW(Dune::NotImplemented, "Acosta::dsw_dpc(params, pc)"); + DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dsw_dpc(params, pc)"); } /*! @@ -148,7 +198,8 @@ public: * \param params A container object that is populated with the appropriate coefficients for the respective law. * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too. */ - static Scalar krw(const Params ¶ms, Scalar swe) + template + static Scalar krw(Scalar swe, const Params& params) { assert(0 <= swe && swe <= 1); Scalar krw; @@ -180,9 +231,10 @@ public: * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too. */ - static Scalar dkrw_dsw(const Params ¶ms, Scalar swe) + template + static Scalar dkrw_dsw(Scalar swe, const Params& params) { - DUNE_THROW(Dune::NotImplemented, "Acosta::dkrw_dsw(params, swe)"); + DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dkrw_dsw(params, swe)"); }; /*! @@ -195,7 +247,8 @@ public: * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too. */ - static Scalar krn(const Params ¶ms, Scalar swe) + template + static Scalar krn(Scalar swe, const Params& params) { assert(0 <= swe && swe <= 1); @@ -229,13 +282,14 @@ public: * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too. */ - static Scalar dkrn_dsw(const Params ¶ms, Scalar swe) + template + static Scalar dkrn_dsw(Scalar swe, const Params& params) { - DUNE_THROW(Dune::NotImplemented, "Acosta::dkrn_dsw(params, swe)"); + DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dkrn_dsw(params, swe)"); } private: - + template static Scalar evaluateIntegral_(const std::vector& intVector, const Scalar swe) { Scalar satPos = (integrationSteps_ - 1) * swe; @@ -257,23 +311,226 @@ private: static std::vector intDraSw1_; }; -template -std::vector AcostaPcSw::intImb0Sw_ = {0.0000000000e+00, 1.4448654771e-08, 2.2819077599e-08, 2.7262869143e-08, 2.9465634517e-08, 3.0504079857e-08, + +template +class AcostaRegularization +{ +public: + template + struct Params + { + Params(S thresholdSw = 1e-3) + : thresholdSw_(thresholdSw) + {} + + S thresholdSw() const + { return thresholdSw_; } + + bool operator== (const Params& p) const + { return Dune::FloatCmp::eq(thresholdSw(), p.thresholdSw(), 1e-6);} + + private: + S thresholdSw_; + }; + + template + void init(const MaterialLaw* m, const BaseParams& bp, const EffToAbsParams& etap, const Params& p) + { + thresholdSw_ = p.thresholdSw(); + + initParameters_(m, bp, thresholdSw_); + } + + /*! + * \brief A regularized Acosta capillary pressure-saturation + * curve. + * + * regularized part: + * - low saturation: extend the \f$p_c(S_w)\f$ curve with the slope at the regularization point (i.e. no kink). + * - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$ by a straight line (yes, there is a kink :-( ). + * + * For the non-regularized part: + * + * \copydetails AcostaPcSw::pc() + */ + OptionalScalar pc(Scalar swe) const + { + // make sure that the capilary pressure observes a derivative != 0 for 'illegal' saturations. This is + // required for example by newton solvers (if the derivative is calculated numerically) in order to get the + // saturation moving to the right direction if it temporarily is in an 'illegal' range. + + if (acE_ == 0.0) //imbibition + { + if (swe < thresholdSw_) + return pcsweLow_ + mLow_*( swe - thresholdSw_); + else if (swe >= (1.0 - thresholdSw_)) + return pcsweHigh_ + mHigh_*(swe - (1.0 - thresholdSw_)); + } + else //drainage + { + if (swe <= thresholdSw_) + return pcsweLow_ + mLow_*(swe - thresholdSw_); + else if (swe >= (1.0 - thresholdSw_)) + return pcsweHigh_ + mHigh_*(swe - (1.0 - thresholdSw_)); + } + + // if the effective saturation is in an 'reasonable' + // range, we use the real Acosta law... + return {}; + } + + /*! + * \brief A regularized Acosta saturation-capillary pressure curve. + * + * function does not exist yet! + */ + OptionalScalar sw(Scalar pc) const + { + DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::sw(bp, pc)"); + } + + /*! + * \brief A regularized version of the partial derivative + * of the \f$p_c(\overline S_w)\f$ w.r.t. effective saturation + * according to Acosta. + * + * regularized part: + * - low saturation: use the slope of the regularization point (i.e. no kink). + * - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$ + * by a straight line and use that slope (yes, there is a kink :-( ). + * + * For the non-regularized part: + * + * \copydetails AcostaPcSw::dpc_dsw() + */ + OptionalScalar dpc_dsw(Scalar swe) const + { + // derivative of the regualarization + if (swe < 0) + return m0_; + else if (swe >= (1.0 - thresholdSw_)) + return mHigh_; + + return {}; + } + + /*! + * \brief The partial derivative of the effective + * saturation to the capillary pressure according to Acosta. + * + * function does not exist yet! + */ + OptionalScalar dsw_dpc(Scalar pc) const + { DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dsw_dpc(pc, bp)"); } + + /*! + * \brief Regularized version of the relative permeability + * for the wetting phase of + * the medium implied by the Acosta + * parameterization. + * + * regularized part: + * - below \f$ \overline S_w =0\f$: set relative permeability to zero + * - above \f$ \overline S_w =1\f$: set relative permeability to one + * - between \f$ 0.95 \leq \overline S_w \leq 1\f$: use a spline as interpolation + * + * For not-regularized part: \copydetails AcostaPcSw::krw() + */ + OptionalScalar krw(Scalar swe) const + { + if (swe <= 0.0) + return 0.0; + else if (swe >= 1.0) + return 1.0; + + return {}; + } + + /*! + * \brief The derivative of the relative permeability for the + * wetting phase in regard to the wetting saturation of the + * medium implied by the Acosta parameterization. + */ + OptionalScalar dkrw_dsw(Scalar swe) const + { DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dkrw_dsw(swe, bp)"); }; + + /*! + * \brief Regularized version of the relative permeability + * for the nonwetting phase of + * the medium implied by the Acosta + * parameterization. + * + * regularized part: + * - below \f$ \overline S_w =0\f$: set relative permeability to zero + * - above \f$ \overline S_w =1\f$: set relative permeability to one + * - for \f$ 0 \leq \overline S_w \leq 0.05 \f$: use a spline as interpolation + * \copydetails AcostaPcSw::krn() + */ + OptionalScalar krn(Scalar swe) const + { + if (swe >= 1.0) + return 0.0; + else if (swe <= 0.0) + return 1.0; + + return {}; + } + + /*! + * \brief The derivative of the relative permeability for the + * nonwetting phase in regard to the wetting saturation of + * the medium as implied by the Acosta + * parameterization. + */ + OptionalScalar dkrn_dsw(Scalar swe) const + { DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dkrn_dsw(swe, bp)"); }; + +private: + template + void initParameters_(const MaterialLaw* m, const BaseParams& bp, const Scalar thresholdSw) + { + thresholdSw_ = thresholdSw; + m0_ = AcostaPcSw::dpc_dsw(0.0, bp); + + mLow_ = AcostaPcSw::dpc_dsw(thresholdSw_, bp); + pcsweLow_ = AcostaPcSw::pc(thresholdSw_, bp); + + mHigh_ = AcostaPcSw::dpc_dsw((1.0-thresholdSw_), bp); + pcsweHigh_ = AcostaPcSw::pc((1.0-thresholdSw_), bp); + + acE_ = bp.acE(); + } + + Scalar thresholdSw_; + Scalar m0_; + Scalar mLow_; + Scalar pcsweLow_; + Scalar mHigh_; + Scalar pcsweHigh_; + Scalar acE_; +}; + +/*! + * \ingroup Fluidmatrixinteractions + * \brief A default configuration for using the VanGenuchten material law + */ +template +using AcostaPcSwDefault = TwoPMaterialLaw, TwoPEffToAbsDefaultPolicy>; + + +std::vector AcostaPcSw::intImb0Sw_ = {0.0000000000e+00, 1.4448654771e-08, 2.2819077599e-08, 2.7262869143e-08, 2.9465634517e-08, 3.0504079857e-08, 3.0976893747e-08, 3.1187237078e-08, 3.1279418104e-08, 3.1319431548e-08, 3.1336696665e-08, 3.1344118639e-08, 3.1347301932e-08, 3.1348665344e-08, 3.1349248803e-08, 3.1349498362e-08, 3.1349605071e-08, 3.1349650691e-08, 3.1349670192e-08, 3.1349678527e-08, 3.1349682090e-08}; -template -std::vector AcostaPcSw::intImbSw1_ = {3.1349682090e-08, 1.6901027319e-08, 8.5306044914e-09, 4.0868129471e-09, 1.8840475735e-09, 8.4560223338e-10, +std::vector AcostaPcSw::intImbSw1_ = {3.1349682090e-08, 1.6901027319e-08, 8.5306044914e-09, 4.0868129471e-09, 1.8840475735e-09, 8.4560223338e-10, 3.7278834340e-10, 1.6244501261e-10, 7.0263986151e-11, 3.0250542321e-11, 1.2985425390e-11, 5.5634511134e-12, 2.3801587239e-12, 1.0167462066e-12, 4.3328717089e-13, 1.8372876392e-13, 7.7019533374e-14, 3.1399697336e-14, 1.1898564029e-14, 3.5629091876e-15, 0.0000000000e+00}; -template -std::vector AcostaPcSw::intDra0Sw_ = {0.0000000000e+00, 9.0286303155e-11, 1.0440347151e-09, 5.3755677032e-09, 2.0922261595e-08, 7.4676307329e-08, +std::vector AcostaPcSw::intDra0Sw_ = {0.0000000000e+00, 9.0286303155e-11, 1.0440347151e-09, 5.3755677032e-09, 2.0922261595e-08, 7.4676307329e-08, 2.8041525597e-07, 1.3662599439e-06, 2.4661875855e-05, 9.2186966645e-01, 9.2188585863e-01, 9.2189413795e-01, 9.2189728035e-01, 9.2189773354e-01, 9.2189777056e-01, 9.2189777338e-01, 9.2189777360e-01, 9.2189777360e-01, 9.2189777359e-01, 9.2189777356e-01, 9.2189777362e-01}; -template -std::vector AcostaPcSw::intDraSw1_ = {9.2189777362e-01, 9.2189777350e-01, 9.2189777259e-01, 9.2189776823e-01, 9.2189775271e-01, 9.2189769898e-01, +std::vector AcostaPcSw::intDraSw1_ = {9.2189777362e-01, 9.2189777350e-01, 9.2189777259e-01, 9.2189776823e-01, 9.2189775271e-01, 9.2189769898e-01, 9.2189749322e-01, 9.2189640736e-01, 9.2187311174e-01, 2.8107160428e-05, 1.1914958378e-05, 3.6356677731e-06, 4.9325835742e-07, 4.0074681584e-08, 3.0250536841e-09, 2.3593717317e-10, 1.8918448614e-11, 1.5374841055e-12, 1.2497552223e-13, 9.4684316652e-15, 0.0000000000e+00}; } diff --git a/lecture/mm/fuelcell/material/acostaparams.hh b/lecture/mm/fuelcell/material/acostaparams.hh deleted file mode 100644 index 11834ba..0000000 --- a/lecture/mm/fuelcell/material/acostaparams.hh +++ /dev/null @@ -1,113 +0,0 @@ -// -*- 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 2 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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Specification of the material parameters - * for the Acosta constitutive relations. - */ -#ifndef ACOSTA_PARAMS_HH -#define ACOSTA_PARAMS_HH - -namespace Dumux -{ -/*! - * - * \brief Specification of the material parameters - * for the Acosta constitutive relations. - */ -template -class AcostaPcSwParams -{ -public: - typedef ScalarT Scalar; - - AcostaPcSwParams() - {} - - AcostaPcSwParams(Scalar acA, Scalar acB, Scalar acC, Scalar acD, Scalar acE) - { setAcA(acA); setAcB(acB); setAcC(acC); setAcD(acD); setAcE(acE); }; - - /*! - * \ Acosta's curve. - */ - Scalar acA() const - { return acA_; } - - /*! - * \ Acosta's curve. - */ - void setAcA(Scalar a) - { acA_ = a; } - - /*! - * \Acosta's curve. - */ - Scalar acB() const - { return acB_; } - - void setAcB(Scalar b) - { acB_ = b; } - - /*! - * \ Acosta's curve. - */ - Scalar acC() const - { return acC_; } - - /*! - * \ Acosta's curve. - */ - void setAcC(Scalar c) - { acC_ = c; } - - /*! - * \ Acosta's curve. - */ - Scalar acD() const - { return acD_; } - - /*! - * \ Acosta's curve. - */ - void setAcD(Scalar d) - { acD_ = d; } - - /*! - * \ Acosta's curve. - */ - Scalar acE() const - { return acE_; } - - /*! - * \ Acosta's curve. - */ - void setAcE(Scalar e) - { acE_ = e; } - -private: - Scalar acA_; - Scalar acB_; - Scalar acC_; - Scalar acD_; - Scalar acE_; -}; -} // namespace Dumux - -#endif diff --git a/lecture/mm/fuelcell/material/regularizedacosta.hh b/lecture/mm/fuelcell/material/regularizedacosta.hh deleted file mode 100644 index b83f20a..0000000 --- a/lecture/mm/fuelcell/material/regularizedacosta.hh +++ /dev/null @@ -1,284 +0,0 @@ -// -*- 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 2 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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Implementation of the capillary pressure and - * water phase saturation according to Acosta. - */ -#ifndef REGULARIZED_ACOSTA_HH -#define REGULARIZED_ACOSTA_HH - -#include "acosta.hh" -#include "regularizedacostaparams.hh" - -//#include -//#include -//#include - -#include - -namespace Dumux -{ -/*! - * \brief Implementation of the regularized Acosta - * capillary pressure / relative permeability <-> saturation relation. - * This class bundles the "raw" curves as - * static members and doesn't concern itself converting - * absolute to effective saturations and vice versa. - * - * In order to avoid very steep gradients the marginal values are "regularized". - * This means that in stead of following the curve of the material law in these regions, some linear approximation is used. - * Doing this is not worse than following the material law. E.g. for very low wetting phase values the material - * laws predict infinite values for \f$p_c\f$ which is completely unphysical. In case of very high wetting phase - * saturations the difference between regularized and "pure" material law is not big. - * - * Regularizing has the additional benefit of being numerically friendly: Newton's method does not like infinite gradients. - * - * The implementation is accomplished as follows: - * - check whether we are in the range of regularization - * - yes: use the regularization - * - no: forward to the standard material law. - * - * For an example figure of the regularization: RegularizedVanGenuchten - * - * \see Acosta - */ - -/*! - * - * For general info: EffToAbsLaw - * - * \see AcostaParams - */ - -template > -class RegularizedAcosta -{ - typedef Dumux::AcostaPcSw Acosta; - -public: - typedef ParamsT Params; - typedef typename Params::Scalar Scalar; - - /*! - * \brief A regularized Acosta capillary pressure-saturation - * curve. - * - * regularized part: - * - low saturation: extend the \f$p_c(S_w)\f$ curve with the slope at the regularization point (i.e. no kink). - * - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$ by a straight line (yes, there is a kink :-( ). - * - * For the non-regularized part: - * - * \copydetails Acosta::pc() - */ - static Scalar pc(const Params ¶ms, Scalar swe) - { -// return 0.0; - const Scalar sThres = params.thresholdSw(); - - // make sure that the capilary pressure observes a - // derivative != 0 for 'illegal' saturations. This is - // required for example by newton solvers (if the - // derivative is calculated numerically) in order to get the - // saturation moving to the right direction if it - // temporarily is in an 'illegal' range. - - //imbibition - if (params.acE() == 0.0) - { - if (swe < sThres) { - Scalar m = Acosta::dpc_dsw(params, sThres); - Scalar pcsweLow = Acosta::pc(params, sThres); - return pcsweLow + m*( swe - sThres); - } - else if (swe >= (1.0 - sThres)) { - Scalar m = Acosta::dpc_dsw(params, (1.0 - sThres)); - Scalar pcsweHigh = Acosta::pc(params, (1.0 - sThres)); - return pcsweHigh + m*(swe - (1.0 - sThres)); - } - } - //drainage - else - { - if (swe <= sThres) { - Scalar m = Acosta::dpc_dsw(params, sThres); - Scalar pcsweLow = Acosta::pc(params, sThres); - return pcsweLow + m*(swe - sThres); - } - else if (swe >= (1.0 - sThres)) { - Scalar m = Acosta::dpc_dsw(params, (1.0 - sThres)); - Scalar pcsweHigh = Acosta::pc(params, (1.0 - sThres)); - return pcsweHigh + m*(swe - (1.0 - sThres)); - } - } - // if the effective saturation is in an 'reasonable' - // range, we use the real Acosta law... - return Acosta::pc(params, swe); - } - - /*! - * \brief A regularized Acosta saturation-capillary pressure curve. - * - * regularized part: - * - low saturation: extend the \f$p_c(S_w)\f$ curve with the slope at the regularization point (i.e. no kink). - * - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$ by a straight line (yes, there is a kink :-( ). - * - * The according quantities are obtained by exploiting theorem of intersecting lines. - * - * For the non-regularized part: - * - * \copydetails Acosta::sw() - */ - static Scalar sw(const Params ¶ms, Scalar pc) - { - DUNE_THROW(Dune::NotImplemented, "Acosta::sw(params, pc)"); - } - - /*! - * \brief A regularized version of the partial derivative - * of the \f$p_c(\overline S_w)\f$ w.r.t. effective saturation - * according to Acosta. - * - * regularized part: - * - low saturation: use the slope of the regularization point (i.e. no kink). - * - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$ by a straight line and use that slope (yes, there is a kink :-( ). - * - * For the non-regularized part: - * - * \copydetails Acosta::dpc_dsw() - */ - static Scalar dpc_dsw(const Params ¶ms, Scalar swe) - { - const Scalar sThres = params.thresholdSw(); - - // derivative of the regualarization - if (swe < 0) { - // calculate the slope of the straight line used in pc() - Scalar m = Acosta::dpc_dsw(params, 0); - return m; - } - else if (swe >= (1.0 - sThres)) { - // calculate the slope of the straight line used in pc() - Scalar m = Acosta::dpc_dsw(params, (1.0 - sThres)); - return m; - } - - return Acosta::dpc_dsw(params, swe); - } - - /*! - * \brief The partial derivative of the effective - * saturation to the capillary pressure according to Acosta. - * - * function does not exist yet! - * - * \param pc Capillary pressure \f$p_C\f$ - * \param params A container object that is populated with the appropriate coefficients for the respective law. - * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container - * is constructed accordingly. Afterwards the values are set there, too. - */ - static Scalar dsw_dpc(const Params ¶ms, Scalar pc) - { - DUNE_THROW(Dune::NotImplemented, "Acosta::dsw_dpc(params, pc)"); - } - - /*! - * \brief Regularized version of the relative permeability - * for the wetting phase of - * the medium implied by the Acosta - * parameterization. - * - * regularized part: - * - below \f$ \overline S_w =0\f$: set relative permeability to zero - * - above \f$ \overline S_w =1\f$: set relative permeability to one - * - between \f$ 0.95 \leq \overline S_w \leq 1\f$: use a spline as interpolation - * - * For not-regularized part: - \copydetails Acosta::krw() - */ - static Scalar krw(const Params ¶ms, Scalar swe) - { - if (swe <= 0.0) - return 0.0; - else if (swe >= 1.0) - return 1.0; - - return Acosta::krw(params, swe); - } - - - /*! - * \brief The derivative of the relative permeability for the - * wetting phase in regard to the wetting saturation of the - * medium implied by the Acosta parameterization. - * - * \param swe The mobile saturation of the wetting phase. - * \param params A container object that is populated with the appropriate coefficients for the respective law. - * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container - * is constructed accordingly. Afterwards the values are set there, too. - */ - static Scalar dkrw_dsw(const Params ¶ms, Scalar swe) - { - DUNE_THROW(Dune::NotImplemented, "Acosta::dkrw_dsw(params, swe)"); - }; - - /*! - * \brief Regularized version of the relative permeability - * for the nonwetting phase of - * the medium implied by the Acosta - * parameterization. - * - * regularized part: - * - below \f$ \overline S_w =0\f$: set relative permeability to zero - * - above \f$ \overline S_w =1\f$: set relative permeability to one - * - for \f$ 0 \leq \overline S_w \leq 0.05 \f$: use a spline as interpolation - * - \copydetails Acosta::krn() - * - */ - static Scalar krn(const Params ¶ms, Scalar swe) - { - if (swe >= 1.0) - return 0.0; - else if (swe <= 0.0) - return 1.0; - - return Acosta::krn(params, swe); - } - - /*! - * \brief The derivative of the relative permeability for the - * nonwetting phase in regard to the wetting saturation of - * the medium as implied by the Acosta - * parameterization. - * - * \param swe The mobile saturation of the wetting phase. - * \param params A container object that is populated with the appropriate coefficients for the respective law. - * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container - * is constructed accordingly. Afterwards the values are set there, too. - */ - static Scalar dkrn_dsw(const Params ¶ms, Scalar swe) - { - DUNE_THROW(Dune::NotImplemented, "Acosta::dkrn_dsw(params, swe)"); - }; - }; -} -#endif diff --git a/lecture/mm/fuelcell/material/regularizedacostaparams.hh b/lecture/mm/fuelcell/material/regularizedacostaparams.hh deleted file mode 100644 index 006c240..0000000 --- a/lecture/mm/fuelcell/material/regularizedacostaparams.hh +++ /dev/null @@ -1,76 +0,0 @@ -// -*- 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 2 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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Specification of the material parameters - * for the Acosta constitutive relations. - */ -#ifndef REGULARIZED_ACOSTA_PARAMS_HH -#define REGULARIZED_ACOSTA_PARAMS_HH - -#include "acostaparams.hh" - -namespace Dumux -{ -/*! - * - * \brief Parameters that are necessary for the \em regularization of - * the Acosta capillary pressure model. - */ -template -class RegularizedAcostaParams : public Dumux::AcostaPcSwParams -{ - typedef Dumux::AcostaPcSwParams AcostaParams; - -public: - typedef ScalarT Scalar; - - RegularizedAcostaParams() : AcostaParams() {} - - RegularizedAcostaParams(Scalar acA, Scalar acB, Scalar acC, Scalar acD, Scalar acE) - : AcostaParams(acA, acB, acC, acD, acE) - {} - - /*! - * \brief Threshold saturation below which the capillary pressure - * is regularized. - * - * This is just 1%. If you need a different value, overload this - * class. - */ - Scalar thresholdSw() const - { - // Some problems are very sensitive to this value - // (e.g. makeing it smaller might result in negative - // pressures), if you change it here, you will almost - // certainly break someone's code! - // - // If you want to use a different regularization threshold, - // overload this class and supply the new class as second - // template parameter for the RegularizedVanGenuchten law! - return /* PLEASE DO _NOT_ */ 1e-3; /* CHANGE THIS VALUE. READ - * COMMENT ABOVE! */ - } - -}; -} // namespace Dumux - -#endif - -- GitLab From 4c5967006f2c9ca7ec3dc0013151b6abf30c47e2 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Thu, 5 Nov 2020 13:57:26 +0100 Subject: [PATCH 15/15] [fuelcells][spatialParams] Use correct residual saturation --- lecture/mm/fuelcell/fuelcell.input | 4 ---- lecture/mm/fuelcell/fuelcellspatialparams.hh | 3 ++- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/lecture/mm/fuelcell/fuelcell.input b/lecture/mm/fuelcell/fuelcell.input index 669178c..464cb9d 100644 --- a/lecture/mm/fuelcell/fuelcell.input +++ b/lecture/mm/fuelcell/fuelcell.input @@ -55,7 +55,3 @@ FreqOutput = 1 # frequency of VTK output SolidDensity = 1430 # [kg/m^3] SolidThermalConductivity = 15.6 # [W/(m*K)] SolidHeatCapacity = 710 # [J/(kg*K)] - -[SpatialParams] -Swr = 0.05 -Snr = 0.05 diff --git a/lecture/mm/fuelcell/fuelcellspatialparams.hh b/lecture/mm/fuelcell/fuelcellspatialparams.hh index a454a1f..44e0f65 100644 --- a/lecture/mm/fuelcell/fuelcellspatialparams.hh +++ b/lecture/mm/fuelcell/fuelcellspatialparams.hh @@ -96,7 +96,8 @@ public: lambdaSolid_ = 15.6; // [W/(m*K)] Acosta thermal conductivity used in capillary pressure-saturation typename PcKrSwCurve::BasicParams params(-1168.75, 8.5, -0.2, -700, 0.0); - pcKrSwCurve_ = std::make_unique(params); + typename PcKrSwCurve::EffToAbsParams effToAbsParams(0.05, 0.05); + pcKrSwCurve_ = std::make_unique(params, effToAbsParams); } -- GitLab