Skip to content
Permalink
Browse files

Merge pull request #13132 from WilkAndy/sat_not_effsat_13131

PorousFlow uses saturation not effective saturation
  • Loading branch information...
permcody committed Mar 28, 2019
2 parents 85f7bd2 + 27cb709 commit 0e86a30c118bcfeb3d783c917d645e53ff7c3155
@@ -7,7 +7,7 @@ and the porous medium. Capillary pressure, $P_c$, is commonly defined as [citep!
P_c = P_{nw} - P_w,
\end{equation}
where $P_{nw}$ is the pressure of the non-wetting phase (typically the gas phase), and $P_w$ is the
pressure of the wetting phase (typically the liquid phase).
pressure of the wetting phase (typically the liquid phase). In the case of a single, unsaturated phase, it is common to use $P_{c} = -P_{w}$, which is the convention employed in PorousFlow. In this case, the capillary pressure is only relevant when $P_{w}\leq 0$, and this results in $S_{\mathrm{eff}}\leq 1$ using the formulae below, otherwise, for $P_{w}>0$, $S_{\mathrm{eff}}=1$.

The capillary pressure is given by the Young-Laplace equation [citep!bear1972]
\begin{equation}
@@ -22,9 +22,15 @@ semi-impirical formulations for capillary pressure have been proposed that relat
to effective saturation.

Several capillary pressure formulations are available in the Porous Flow module suitable for either
porepressure formulations (where effective saturation is calculated using capillary pressure), or
porepressure formulations (where effective saturation and then saturation is calculated using capillary pressure), or
porepressure-saturation formulations (where capillary pressure is calculated using the saturation).

In the following, the effective saturation is
\begin{equation}
S_{\mathrm{eff}} = \frac{S - S_{\mathrm{l,res}}}{1 - S_{\mathrm{l,res}}} \ .
\end{equation}
This is the wetting-phase (liquid) effective saturation, and $S$ is the wetting phase saturation. The quantity $S_{\mathrm{l,res}}$ is the wetting-phase residual saturation, and is termed `sat_lr` in MOOSE input files. The porepressure formulations compute $S_{\mathrm{eff}}$ as a function of capillary pressure, and hence $S$ which then appears in the fluid mass, the relative permeability, etc.

## Constant

[`PorousFlowCapillaryPressureConst`](/PorousFlowCapillaryPressureConst.md)
@@ -165,15 +171,17 @@ Only effective saturation as a function of capillary pressure is available in
## Logarithmic extension at low liquid saturations

Several of the capillary pressure formulations have capillary pressure and its derivative approach
infinity while the liquid saturation approaches zero. While this is desirable for calculation of
infinity while the liquid effective saturation approaches zero. While this is desirable for calculation of
effective saturation as a function of capillary pressure, it is undesirable when calculating the
capillary pressure numerically using saturation, often leading to numerical convergence issues.

By default, the numerical implementations of the capillary pressure curves implement a hard maximum
when the effective saturation decreases below 0 in order to avoid unphysical values of capillary
pressure and its derivatives. While this approach does avoid infinite values, it can lead to
numerical difficulties for saturations close to residual due to the discontinuous derivative of
capillary pressure with respect to saturation.
By default, the numerical implementations of the capillary pressure
curves implement a hard maximum when the effective saturation
decreases below 0 in order to avoid unphysical values of capillary
pressure and its derivatives. While this approach does avoid infinite
values, it can lead to numerical difficulties for saturations close to
residual (effective saturation close to zero) due to the discontinuous
derivative of capillary pressure with respect to saturation.

To overcome this, the logarithmic extension detailed by [cite!webb2000] is implemented for low
saturations in formulations where capillary pressure approaches infinity for small liquid
@@ -184,7 +192,7 @@ P_c = P_{c,max} 10^{m(S - S^*)}
is used for saturations less than a value $S^*$. The value of $s^*$ is calculated so that the
capillary pressure curve is continuous and smooth up to the maximum capillary pressure $P_{c,max}$,
see [pc_vg_logext] for an example for the van Genuchten capillary pressure, and [pc_bc_logext] for
the Brooks-Corey capillary pressure.
the Brooks-Corey capillary pressure. For these figures, $S_{\mathrm{l,res}}=0.1$, and $S^{*}\approx 0.11$ (VG) or $S^{*}\approx 0.101$ (BC).
!media media/porous_flow/pc_vg_logext.png
id=pc_vg_logext
@@ -2,6 +2,8 @@

!syntax description /Materials/PorousFlow1PhaseP

Given the porepressure, this Material computes the saturation for 1-phase fluid simulations. See [capillary pressure](capillary_pressure.md) for further details.

!syntax parameters /Materials/PorousFlow1PhaseP

!syntax inputs /Materials/PorousFlow1PhaseP
@@ -2,6 +2,8 @@

!syntax description /Materials/PorousFlow2PhasePP

Given the liquid and gas porepressures, this Material computes the liquid and gas saturations for 2-phase fluid simulations. See [capillary pressure](capillary_pressure.md) for further details.

!syntax parameters /Materials/PorousFlow2PhasePP

!syntax inputs /Materials/PorousFlow2PhasePP
@@ -100,6 +100,33 @@ class PorousFlowCapillaryPressure : public DiscreteElementUserObject
*/
virtual Real d2EffectiveSaturation(Real pc, unsigned qp = 0) const = 0;

/**
* Saturation as a function of capillary pressure
* @param pc capillary pressure (Pa)
* @param qp quadpoint to use (when saturation depends on coupled variables, not just
* pc)
* @return saturation
*/
Real saturation(Real pc, unsigned qp = 0) const;

/**
* Derivative of saturation wrt capillary pressure
* @param pc capillary pressure (Pa)
* @param qp quadpoint to use (when saturation depends on coupled variables, not just
* pc)
* @return derivative of saturation wrt capillary pressure
*/
Real dSaturation(Real pc, unsigned qp = 0) const;

/**
* Second derivative of saturation wrt capillary pressure
* @param pc capillary pressure
* @param qp quadpoint to use (when saturation depends on coupled variables, not just
* pc)
* @return second derivative of saturation wrt capillary pressure
*/
Real d2Saturation(Real pc, unsigned qp = 0) const;

protected:
/**
* Effective saturation of liquid phase given liquid saturation and residual
@@ -59,26 +59,26 @@ PorousFlow1PhaseP::computeQpProperties()
PorousFlowVariableBase::computeQpProperties();

buildQpPPSS();
const Real dseff = _pc_uo.dEffectiveSaturation(_porepressure_var[_qp]);
const Real ds = _pc_uo.dSaturation(_porepressure_var[_qp]);

if (!_nodal_material)
{
(*_gradp_qp)[_qp][0] = _gradp_qp_var[_qp];
(*_grads_qp)[_qp][0] = dseff * _gradp_qp_var[_qp];
(*_grads_qp)[_qp][0] = ds * _gradp_qp_var[_qp];
}

// _porepressure is only dependent on _porepressure, and its derivative is 1
if (_dictator.isPorousFlowVariable(_porepressure_varnum))
{
// _porepressure is a PorousFlow variable
_dporepressure_dvar[_qp][0][_p_var_num] = 1.0;
_dsaturation_dvar[_qp][0][_p_var_num] = dseff;
_dsaturation_dvar[_qp][0][_p_var_num] = ds;
if (!_nodal_material)
{
(*_dgradp_qp_dgradv)[_qp][0][_p_var_num] = 1.0;
(*_dgrads_qp_dgradv)[_qp][0][_p_var_num] = dseff;
(*_dgrads_qp_dgradv)[_qp][0][_p_var_num] = ds;
(*_dgrads_qp_dv)[_qp][0][_p_var_num] =
_pc_uo.d2EffectiveSaturation(_porepressure_var[_qp]) * _gradp_qp_var[_qp];
_pc_uo.d2Saturation(_porepressure_var[_qp]) * _gradp_qp_var[_qp];
}
}
}
@@ -87,5 +87,5 @@ void
PorousFlow1PhaseP::buildQpPPSS()
{
_porepressure[_qp][0] = _porepressure_var[_qp];
_saturation[_qp][0] = _pc_uo.effectiveSaturation(_porepressure_var[_qp]);
_saturation[_qp][0] = _pc_uo.saturation(_porepressure_var[_qp]);
}
@@ -72,13 +72,13 @@ PorousFlow2PhasePP::computeQpProperties()
PorousFlowVariableBase::computeQpProperties();

const Real pc = buildQpPPSS();
const Real dseff = _pc_uo.dEffectiveSaturation(pc); // d(seff)/d(pc)
const Real ds = _pc_uo.dSaturation(pc); // dS/d(pc)

if (!_nodal_material)
{
(*_gradp_qp)[_qp][0] = _phase0_gradp_qp[_qp];
(*_gradp_qp)[_qp][1] = _phase1_gradp_qp[_qp];
(*_grads_qp)[_qp][0] = dseff * ((*_gradp_qp)[_qp][0] - (*_gradp_qp)[_qp][1]);
(*_grads_qp)[_qp][0] = ds * ((*_gradp_qp)[_qp][0] - (*_gradp_qp)[_qp][1]);
(*_grads_qp)[_qp][1] = -(*_grads_qp)[_qp][0];
}

@@ -99,35 +99,31 @@ PorousFlow2PhasePP::computeQpProperties()

if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
{
_dsaturation_dvar[_qp][0][_p0var] = dseff;
_dsaturation_dvar[_qp][1][_p0var] = -dseff;
_dsaturation_dvar[_qp][0][_p0var] = ds;
_dsaturation_dvar[_qp][1][_p0var] = -ds;
}
if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
{
_dsaturation_dvar[_qp][0][_p1var] = -dseff;
_dsaturation_dvar[_qp][1][_p1var] = dseff;
_dsaturation_dvar[_qp][0][_p1var] = -ds;
_dsaturation_dvar[_qp][1][_p1var] = ds;
}

if (!_nodal_material)
{
const Real d2seff_qp = _pc_uo.d2EffectiveSaturation(pc); // d^2(seff_qp)/d(pc_qp)^2
const Real d2s_qp = _pc_uo.d2Saturation(pc); // d^2(S_qp)/d(pc_qp)^2
if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
{
(*_dgrads_qp_dgradv)[_qp][0][_p0var] = dseff;
(*_dgrads_qp_dv)[_qp][0][_p0var] =
d2seff_qp * (_phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp]);
(*_dgrads_qp_dgradv)[_qp][1][_p0var] = -dseff;
(*_dgrads_qp_dv)[_qp][1][_p0var] =
-d2seff_qp * (_phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp]);
(*_dgrads_qp_dgradv)[_qp][0][_p0var] = ds;
(*_dgrads_qp_dv)[_qp][0][_p0var] = d2s_qp * (_phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp]);
(*_dgrads_qp_dgradv)[_qp][1][_p0var] = -ds;
(*_dgrads_qp_dv)[_qp][1][_p0var] = -d2s_qp * (_phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp]);
}
if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
{
(*_dgrads_qp_dgradv)[_qp][0][_p1var] = -dseff;
(*_dgrads_qp_dv)[_qp][0][_p1var] =
-d2seff_qp * (_phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp]);
(*_dgrads_qp_dgradv)[_qp][1][_p1var] = dseff;
(*_dgrads_qp_dv)[_qp][1][_p1var] =
d2seff_qp * (_phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp]);
(*_dgrads_qp_dgradv)[_qp][0][_p1var] = -ds;
(*_dgrads_qp_dv)[_qp][0][_p1var] = -d2s_qp * (_phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp]);
(*_dgrads_qp_dgradv)[_qp][1][_p1var] = ds;
(*_dgrads_qp_dv)[_qp][1][_p1var] = d2s_qp * (_phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp]);
}
}
}
@@ -138,8 +134,8 @@ PorousFlow2PhasePP::buildQpPPSS()
_porepressure[_qp][0] = _phase0_porepressure[_qp];
_porepressure[_qp][1] = _phase1_porepressure[_qp];
const Real pc = _phase0_porepressure[_qp] - _phase1_porepressure[_qp]; // this is <= 0
const Real seff = _pc_uo.effectiveSaturation(pc);
_saturation[_qp][0] = seff;
_saturation[_qp][1] = 1.0 - seff;
const Real sat = _pc_uo.saturation(pc);
_saturation[_qp][0] = sat;
_saturation[_qp][1] = 1.0 - sat;
return pc;
}
@@ -93,6 +93,24 @@ PorousFlowCapillaryPressure::effectiveSaturationFromSaturation(Real saturation)
return (saturation - _sat_lr) / (1.0 - _sat_lr);
}

Real
PorousFlowCapillaryPressure::saturation(Real pc, unsigned qp) const
{
return effectiveSaturation(pc, qp) * (1.0 - _sat_lr) + _sat_lr;
}

Real
PorousFlowCapillaryPressure::dSaturation(Real pc, unsigned qp) const
{
return dEffectiveSaturation(pc, qp) * (1.0 - _sat_lr);
}

Real
PorousFlowCapillaryPressure::d2Saturation(Real pc, unsigned qp) const
{
return d2EffectiveSaturation(pc, qp) * (1.0 - _sat_lr);
}

Real
PorousFlowCapillaryPressure::capillaryPressureLogExt(Real saturation) const
{
@@ -163,11 +163,15 @@
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.11
[../]
[./relperm1]
type = PorousFlowRelativePermeabilityCorey
n = 3
phase = 1
s_res = 0.01
sum_s_res = 0.11
[../]
[./relperm0_nodal]
type = PorousFlowRelativePermeabilityCorey
@@ -135,11 +135,15 @@
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.11
[../]
[./relperm1]
type = PorousFlowRelativePermeabilityCorey
n = 3
phase = 1
s_res = 0.01
sum_s_res = 0.11
[../]
[./porosity]
type = PorousFlowPorosityConst
@@ -91,6 +91,8 @@
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.1
[../]
[./darcy_velocity_qp]
type = PorousFlowDarcyVelocityMaterial
@@ -92,6 +92,8 @@
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.1
[../]
[./darcy_velocity_qp]
type = PorousFlowDarcyVelocityMaterial
@@ -100,6 +100,8 @@
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.1
[../]
[./darcy_velocity_qp]
type = PorousFlowDarcyVelocityMaterial
@@ -129,11 +129,15 @@
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.1
[../]
[./relperm1_qp]
type = PorousFlowRelativePermeabilityCorey
n = 3
phase = 1
s_res = 0.0
sum_s_res = 0.1
[../]
[./darcy_velocity_qp]
type = PorousFlowDarcyVelocityMaterial
@@ -80,7 +80,7 @@
m = 0.5
alpha = 1
pc_max = 10
sat_lr = 0.1
sat_lr = 0.01
[../]
[]

@@ -85,7 +85,7 @@
m = 0.5
alpha = 1
pc_max = 10
sat_lr = 0.1
sat_lr = 0.01
[../]
[]

@@ -166,12 +166,16 @@
at_nodes = true
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.1
[../]
[./relperm1]
type = PorousFlowRelativePermeabilityCorey
at_nodes = true
n = 3
phase = 1
s_res = 0.0
sum_s_res = 0.1
[../]
[]

0 comments on commit 0e86a30

Please sign in to comment.
You can’t perform that action at this time.