Commit 393ec3a5 authored by Ned Coltman's avatar Ned Coltman
Browse files

[mm][naplinfiltration] adapt to new material laws

parent 62a573da
......@@ -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
......
......@@ -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<class MaterialLawParams>
static Scalar invertPcgw_(Scalar pcIn, const MaterialLawParams &pcParams)
template<class FMInteraction>
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;
......
......@@ -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
......
......@@ -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<class MaterialLawParams>
static Scalar invertPcgw_(Scalar pcIn, const MaterialLawParams &pcParams)
// small solver inverting the pc curve
template<class FMInteraction>
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;
......
......@@ -26,11 +26,10 @@
#include <dumux/porousmediumflow/properties.hh>
#include <dumux/material/spatialparams/fv.hh>
#include <dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3p.hh>
#include <dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3pparams.hh>
#include <dumux/material/fluidmatrixinteractions/3p/efftoabslaw.hh>
#include <dumux/material/fluidmatrixinteractions/3p/parkervangenuchten.hh>
#include <dumux/material/fluidmatrixinteractions/3p/napladsorption.hh>
#include <dumux/io/gnuplotinterface.hh>
#include <dumux/io/plotmateriallaw3p.hh>
#include <dumux/common/enumerate.hh>
namespace Dumux {
/*!
......@@ -47,19 +46,14 @@ class InfiltrationSpatialParams
using Element = typename GridView::template Codim<0>::Entity;
using ParentType = FVSpatialParams<FVGridGeometry, Scalar,
InfiltrationSpatialParams<FVGridGeometry, Scalar>>;
using GlobalPosition = typename SubControlVolume::GlobalPosition;
using EffectiveLaw = RegularizedParkerVanGen3P<Scalar>;
using ThreePhasePcKrSw = FluidMatrix::ParkerVanGenuchten3PDefault<Scalar>;
using AdsorptionModel = FluidMatrix::ThreePNAPLAdsorption<Scalar>;
public:
// export permeability type
using PermeabilityType = Scalar;
//get the material law from the property system
using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
using MaterialLawParams = typename MaterialLaw::Params;
/*!
* \brief The constructor
*
......@@ -67,32 +61,14 @@ public:
*/
InfiltrationSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
, threePhasePcKrSw_("SpatialParams")
, adsorption_("SpatialParams")
{
// intrinsic permeabilities
permeability_ = getParam<Scalar>("SpatialParams.permeability");
// porosities
porosity_ = getParam<Scalar>("SpatialParams.porosity");
vGAlpha_ = getParam<Scalar>("SpatialParams.vanGenuchtenAlpha");
vGN_ = getParam<Scalar>("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<bool>("Output.PlotFluidMatrixInteractions");
}
......@@ -105,15 +81,41 @@ public:
{
GnuplotInterface<Scalar> gnuplot(plotFluidMatrixInteractions_);
gnuplot.setOpenPlotWindow(plotFluidMatrixInteractions_);
PlotMaterialLaw<Scalar, MaterialLaw> 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_;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment