Skip to content
Permalink
Browse files

Merge pull request #13141 from WilkAndy/nw_vg_relperm_13137

 Nonwetting vanGenuchten relperm in PorousFlow
  • Loading branch information...
permcody committed Mar 27, 2019
2 parents 60b7eaf + 06da7ee commit 027cbad9eb85b6707d8433bd5a697a70479e2aae
Binary file not shown.
@@ -3,14 +3,18 @@
The relative permeability of a phase is a function of its effective
saturation:
\begin{equation}
S_{\mathrm{eff}}(S) = \frac{S - S_{\mathrm{res}}^{\beta}}{1 -
S^{\beta}_{\mathrm{eff}}(S^{\beta}) = \frac{S^{\beta} - S_{\mathrm{res}}^{\beta}}{1 -
\sum_{\beta'}S_{\mathrm{res}}^{\beta'}}
\end{equation}
In this equation $S_{\mathrm{res}}^{\beta}$ is the residual
saturation for phase $\beta$. If $S_{\mathrm{eff}} < 0$ then the
relative permeability is zero, while if $S_{\mathrm{eff}}>1$ then the
relative permeability is unity. Otherwise, the relative permeability
is given by the expressions below.
In this equation:

- $S^{\beta}_{\mathrm{eff}}$ is the effective saturation for phase $\beta$;
- $S^{\beta}$ is the saturation for phase $\beta$;
- $S_{\mathrm{res}}^{\beta}$ is the residual saturation for phase $\beta$.

In the MOOSE input file, $S_{\mathrm{res}}$ is termed `s_res`, and the entire sum $\sum_{\beta'}S_{\mathrm{res}}^{\beta'}$ is termed `sum_s_res`. MOOSE deliberately does not check that `sum_s_res` equals $\sum$`s_res`, so that you may choose `sum_s_res` independently of your `s_res` quantities.

If $S_{\mathrm{eff}} < 0$ then the relative permeability is zero, while if $S_{\mathrm{eff}}>1$ then the relative permeability is unity. Otherwise, the relative permeability is given by the expressions below.

## Constant

@@ -42,7 +46,7 @@ where $n$ is a user-defined quantity. Originally, [cite!corey1954] used $n = 4$

[`PorousFlowRelativePermeabilityVG`](/PorousFlowRelativePermeabilityVG.md)

The relative permeability of the phase is given by [cite!vangenuchten1980]
The relative permeability of the wetting phase is given by [cite!vangenuchten1980]
\begin{equation}
k_{\mathrm{r}} = \sqrt{S_{\mathrm{eff}}} \left(1 - (1 -
S_{\mathrm{eff}}^{1/m})^{m} \right)^{2}.
@@ -65,6 +69,14 @@ Here the cubic is chosen so that its value and derivative match the
van Genuchten expression at $S=S_{c}$, and so that it is unity at
$S_{\mathrm{eff}}=1$.

For the non-wetting phase, the van-Genuchten expression is
\begin{equation}
k_{\mathrm{r}} = \sqrt{S_{\mathrm{eff}}} \left(1 - (1 -
S_{\mathrm{eff}})^{1/m} \right)^{2m} \ .
\end{equation}
As always in this page, here $S_{\mathrm{eff}}$ is the effective saturation of the appropriate phase, which is the non-wetting phase in this case.


!media media/porous_flow/relperm_vg.png style=width:100%;margin-left:10px; caption=van Genuchten relative permeability id=relperm_vg

## Brooks-Corey
@@ -23,7 +23,7 @@ def corey(s, sr, sum_s_r, n):

return relperm

# Analytical form of van Genuchten's relative permeability
# Analytical form of van Genuchten's wetting relative permeability
def vg(s, sr, sls, m):

# Effective saturation
@@ -35,6 +35,18 @@ def vg(s, sr, sls, m):
# Relative permeability is then
return relperm

# Analytical form of van Genuchten's nonwetting relative permeability
def vg_nw(s, sr, sls, m):

# Effective saturation
seff = np.clip((s - sr) / (sls - sr), 0, 1)

# The gas relative permeability is
relperm = np.sqrt(seff) * np.power(1 - np.power(1 - seff, 1.0 / m), 2.0 * m)

# Relative permeability is then
return relperm

# Saturation of phase 0 varies linearly from 0 to 1
s0 = np.linspace(0, 1, 200)

@@ -45,7 +57,7 @@ def vg(s, sr, sls, m):
# Case 1: residual saturation set to 0 for both phases, n = 1
#
# Read MOOSE simulation data
data = np.genfromtxt('../../tests/relperm/corey1_out_vpp_0001.csv', delimiter = ',', names = True, dtype = None)
data = np.genfromtxt('../../test/tests/relperm/gold/corey1_out_vpp_0001.csv', delimiter = ',', names = True, dtype = None)

plt.figure(1)
plt.plot(s0, corey(s0, 0, 0, 1), label = 'kr0')
@@ -61,7 +73,7 @@ def vg(s, sr, sls, m):
# Case 2: residual saturation set to 0 for both phases, and n = 2
#
# Read MOOSE simulation data
data = np.genfromtxt('../../tests/relperm/corey2_out_vpp_0001.csv', delimiter = ',', names = True, dtype = None)
data = np.genfromtxt('../../test/tests/relperm/gold/corey2_out_vpp_0001.csv', delimiter = ',', names = True, dtype = None)

plt.figure(2)
plt.plot(s0, corey(s0, 0, 0, 2), label = 'kr0')
@@ -77,7 +89,7 @@ def vg(s, sr, sls, m):
# Case 3: residual saturation set to 0.2 for phase 0, 0.3 for phase 1 and n = 2
#
# Read MOOSE simulation data
data = np.genfromtxt('../../tests/relperm/corey3_out_vpp_0001.csv', delimiter = ',', names = True, dtype = None)
data = np.genfromtxt('../../test/tests/relperm//gold/corey3_out_vpp_0001.csv', delimiter = ',', names = True, dtype = None)

plt.figure(3)
plt.plot(s0, corey(s0, 0.2, 0.5, 2), label = 'kr0')
@@ -98,33 +110,34 @@ def vg(s, sr, sls, m):
# Case 1: residual saturation set to 0 for both phases, m = 0.5, sls = 1
#
# Read MOOSE simulation data
data = np.genfromtxt('../../tests/relperm/vangenuchten1_out_vpp_0001.csv', delimiter = ',', names = True, dtype = None)
data = np.genfromtxt('../../test/tests/relperm/gold/vangenuchten1_out_vpp_0001.csv', delimiter = ',', names = True, dtype = None)

plt.figure(4)
plt.plot(s0, vg(s0, 0, 1, 0.5), label = 'kr0')
plt.plot(data['s0aux'], data['kr0aux'], 'ob', label = 'kr0 (MOOSE)')
plt.plot(s0, corey(1 - s0, 0, 0, 2), label = 'kr1')
plt.plot(data['s0aux'], data['kr1aux'], 'og', label = 'kr1 (MOOSE)')
plt.xlabel('Phase 0 saturation (-)')
plt.plot(s0, vg(s0, 0, 1, 0.5), label = 'kr_w')
plt.plot(data['s0aux'], data['kr0aux'], 'ob', label = 'kr_w (MOOSE)')
plt.plot(s0, vg_nw(1 - s0, 0, 1, 0.5), label = 'kr_nw')
plt.plot(data['s0aux'], data['kr1aux'], 'og', label = 'kr_nw (MOOSE)')
plt.xlabel('Wetting phase saturation (-)')
plt.ylabel('Relative permeability (-)')
plt.legend(loc = 'best')
plt.title('van Genuchten relative permeability: $S_{0r} = 0, S_{1r} = 0, m = 0.5$')
plt.title('van Genuchten relative permeability: $S_{w,res} = 0, S_{nw,res} = 0, m = 0.5$')
plt.ylim([-0.01, 1.01])
plt.savefig("vg1_fig.pdf")

# Case 2: residual saturation set to 0.25 for phase 0, 0 for phase 1, m = 0.4, sls = 1
# Case 2: residual saturation set to 0.1 for phase 0, 0.2 for phase 1, m = 0.4
#
# Read MOOSE simulation data
data = np.genfromtxt('../../tests/relperm/vangenuchten2_out_vpp_0001.csv', delimiter = ',', names = True, dtype = None)
data = np.genfromtxt('../../test/tests/relperm/gold/vangenuchten2_out_vpp_0001.csv', delimiter = ',', names = True, dtype = None)

plt.figure(5)
plt.plot(s0, vg(s0, 0.25, 1, 0.4), label = 'kr0')
plt.plot(data['s0aux'], data['kr0aux'], 'ob', label = 'kr0 (MOOSE)')
plt.plot(s0, corey(1 - s0, 0, 0.25, 2), label = 'kr1')
plt.plot(data['s0aux'], data['kr1aux'], 'og', label = 'kr1 (MOOSE)')
plt.plot(s0, vg(s0, 0.1, 0.8, 0.4), label = 'kr_w')
plt.plot(data['s0aux'], data['kr0aux'], 'ob', label = 'kr_w (MOOSE)')
plt.plot(s0, vg_nw(1 - s0, 0.2, 0.9, 0.4), label = 'kr_nw')
plt.plot(data['s0aux'], data['kr1aux'], 'og', label = 'kr_nw (MOOSE)')
plt.xlabel('Phase 0 saturation (-)')
plt.ylabel('Relative permeability (-)')
plt.legend(loc = 'best')
plt.title('van Genuchten relative permeability: $S_{0r} = 0.25, S_{1r} = 0, m = 0.4$')
plt.title('van Genuchten relative permeability: $S_{w,res} = 0.1, S_{nw,res} = 0.2, m = 0.4$')
plt.ylim([-0.01, 1.01])
plt.savefig("vg2_fig.pdf")
plt.savefig("../../doc/content/media/porous_flow/relperm_vg.png")
@@ -42,6 +42,9 @@ class PorousFlowRelativePermeabilityVG : public PorousFlowRelativePermeabilityBa
/// van Genuchten exponent m for the specified phase
const Real _m;

/// Whether to use the wetting or non-wetting van Genuchten expression
const bool _wetting;

/// Start of cubic smoothing
const Real _cut;

@@ -113,6 +113,31 @@ Real dRelativePermeability(Real seff, Real m);
* @return second derivative of relative permeability wrt effective saturation
*/
Real d2RelativePermeability(Real seff, Real m);

/**
* Relative permeability for a non-wetting phase as a function of effective saturation
* @param seff effective saturation
* @param m van Genuchten exponent
* @return relative permeability
*/
Real relativePermeabilityNW(Real seff, Real m);

/**
* Derivative of relative permeability for a non-wetting phase with respect to effective saturation
* @param seff effective saturation
* @param m van Genuchten exponent
* @return derivative of relative permeability wrt effective saturation
*/
Real dRelativePermeabilityNW(Real seff, Real m);

/**
* Second derivative of relative permeability for a non-wetting phase with respect to effective
* saturation
* @param seff effective saturation
* @param m van Genuchten exponent
* @return second derivative of relative permeability wrt effective saturation
*/
Real d2RelativePermeabilityNW(Real seff, Real m);
}

#endif // POROUSFLOWVANGENUCHTEN_H
@@ -28,23 +28,31 @@ validParams<PorousFlowRelativePermeabilityVG>()
params.addParam<bool>("zero_derivative",
false,
"Employ a cubic for seff>seff_turnover that has zero derivative at seff=1");
params.addClassDescription(
"This Material calculates relative permeability of a phase using the van Genuchten model");
params.addParam<bool>("wetting",
true,
"If true, use the van Genuchten form appropriate for a wetting (liquid) "
"phase. If false, use the non-wetting (gas) expression.");
params.addClassDescription("This Material calculates relative permeability of a phase "
"using the van Genuchten model");
return params;
}

PorousFlowRelativePermeabilityVG::PorousFlowRelativePermeabilityVG(
const InputParameters & parameters)
: PorousFlowRelativePermeabilityBase(parameters),
_m(getParam<Real>("m")),
_wetting(getParam<bool>("wetting")),
_cut(getParam<Real>("seff_turnover")),
_cub0(PorousFlowVanGenuchten::relativePermeability(_cut, _m)),
_cub1(PorousFlowVanGenuchten::dRelativePermeability(_cut, _m)),
_cub0(_wetting ? PorousFlowVanGenuchten::relativePermeability(_cut, _m)
: PorousFlowVanGenuchten::relativePermeabilityNW(_cut, _m)),
_cub1(_wetting ? PorousFlowVanGenuchten::dRelativePermeability(_cut, _m)
: PorousFlowVanGenuchten::dRelativePermeabilityNW(_cut, _m)),
_cub2(_cut < 1.0
? (getParam<bool>("zero_derivative")
? 3.0 * (1.0 - _cub0 - _cub1 * (1.0 - _cut)) / Utility::pow<2>(1.0 - _cut) +
_cub1 / (1.0 - _cut)
: PorousFlowVanGenuchten::d2RelativePermeability(_cut, _m))
: (_wetting ? PorousFlowVanGenuchten::d2RelativePermeability(_cut, _m)
: PorousFlowVanGenuchten::d2RelativePermeabilityNW(_cut, _m)))
: 0.0),
_cub3(_cut < 1.0
? (getParam<bool>("zero_derivative")
@@ -60,7 +68,12 @@ Real
PorousFlowRelativePermeabilityVG::relativePermeability(Real seff) const
{
if (seff < _cut)
return PorousFlowVanGenuchten::relativePermeability(seff, _m);
{
if (_wetting)
return PorousFlowVanGenuchten::relativePermeability(seff, _m);
else
return PorousFlowVanGenuchten::relativePermeabilityNW(seff, _m);
}

return _cub0 + _cub1 * (seff - _cut) + _cub2 * Utility::pow<2>(seff - _cut) +
_cub3 * Utility::pow<3>(seff - _cut);
@@ -70,7 +83,12 @@ Real
PorousFlowRelativePermeabilityVG::dRelativePermeability(Real seff) const
{
if (seff < _cut)
return PorousFlowVanGenuchten::dRelativePermeability(seff, _m);
{
if (_wetting)
return PorousFlowVanGenuchten::dRelativePermeability(seff, _m);
else
return PorousFlowVanGenuchten::dRelativePermeabilityNW(seff, _m);
}

return _cub1 + 2.0 * _cub2 * (seff - _cut) + 3.0 * _cub3 * Utility::pow<2>(seff - _cut);
}
@@ -156,4 +156,51 @@ d2RelativePermeability(Real seff, Real m)
return -0.25 * std::pow(seff, -1.5) * Utility::pow<2>(b) + 2.0 * std::pow(seff, -0.5) * b * db +
2.0 * std::sqrt(seff) * db * db + 2.0 * std::sqrt(seff) * b * d2b;
}

Real
relativePermeabilityNW(Real seff, Real m)
{
if (seff <= 0.0)
return 0.0;
else if (seff >= 1.0)
return 1.0;

const Real a = std::pow(1.0 - seff, 1.0 / m);
const Real b = std::pow(1.0 - a, 2.0 * m);

return std::sqrt(seff) * b;
}

Real
dRelativePermeabilityNW(Real seff, Real m)
{
// Guard against division by zero
if (seff <= 0.0 || seff >= 1.0)
return 0.0;

const Real a = std::pow(1.0 - seff, 1.0 / m);
const Real da = -1.0 / m * a / (1.0 - seff);
const Real b = std::pow(1.0 - a, 2.0 * m);
const Real db = -2.0 * m * b / (1.0 - a) * da;

return 0.5 * std::pow(seff, -0.5) * b + std::sqrt(seff) * db;
}

Real
d2RelativePermeabilityNW(Real seff, Real m)
{
// Guard against division by zero
if (seff <= 0.0 || seff >= 1.0)
return 0.0;

const Real a = std::pow(1.0 - seff, 1.0 / m);
const Real da = -1.0 / m * a / (1.0 - seff);
const Real d2a = 1.0 / m * (1.0 / m - 1) * std::pow(1.0 - seff, 1.0 / m - 2.0);
const Real b = std::pow(1.0 - a, 2.0 * m);
const Real db = -2.0 * m * b / (1.0 - a) * da;
const Real d2b =
-2.0 * m * (db / (1.0 - a) * da + b * Utility::pow<2>(da / (1.0 - a)) + b / (1.0 - a) * d2a);

return -0.25 * std::pow(seff, -1.5) * b + std::pow(seff, -0.5) * db + std::sqrt(seff) * d2b;
}
}
@@ -1,22 +1,21 @@
x,y,z,id,s0aux,s1aux,kr0aux,kr1aux
0,0,0,0,0.995,0.005,0.81767127272099,4.3064800892356e-11
0.052631578947368,0,0,0.052631578947368,0.945,0.055,0.44044507951165,5.4900521340933e-07
0.10526315789474,0,0,0.10526315789474,0.895,0.105,0.2903772133122,9.9607515954477e-06
0.15789473684211,0,0,0.15789473684211,0.845,0.155,0.19901105764177,5.7663441077274e-05
0.21052631578947,0,0,0.21052631578947,0.785,0.215,0.12830647743282,0.0002539594790149
0.26315789473684,0,0,0.26315789473684,0.735,0.265,0.088872460313874,0.00065862579390779
0.31578947368421,0,0,0.31578947368421,0.685,0.315,0.061001763100158,0.0014555937764247
0.36842105263158,0,0,0.36842105263158,0.635,0.365,0.041248157468635,0.0028772971835911
0.42105263157895,0,0,0.42105263157895,0.575,0.425,0.02508173623612,0.0058620826279179
0.47368421052632,0,0,0.47368421052632,0.525,0.475,0.016068838953146,0.0099302872835256
0.52631578947368,0,0,0.52631578947368,0.475,0.525,0.0099302872835256,0.016068838953146
0.57894736842105,0,0,0.57894736842105,0.425,0.575,0.0058620826279179,0.025081736236119
0.63157894736842,0,0,0.63157894736842,0.365,0.635,0.0028772971835911,0.041248157468635
0.68421052631579,0,0,0.68421052631579,0.315,0.685,0.0014555937764247,0.061001763100158
0.73684210526316,0,0,0.73684210526316,0.265,0.735,0.00065862579390779,0.088872460313874
0.78947368421053,0,0,0.78947368421053,0.215,0.785,0.0002539594790149,0.12830647743282
0.84210526315789,0,0,0.84210526315789,0.155,0.845,5.7663441077274e-05,0.19901105764177
0.89473684210526,0,0,0.89473684210526,0.105,0.895,9.9607515954477e-06,0.2903772133122
0.94736842105263,0,0,0.94736842105263,0.055,0.945,5.4900521340933e-07,0.44044507951165
1,0,0,1,0.005,0.995,4.3064800892356e-11,0.81767127272099

0,0,0,0,0.995,0.005,0.81767127272099,0.00079468337511365
0.052631578947368,0,0,0.052631578947368,0.945,0.055,0.44044507951165,0.025110851800172
0.10526315789474,0,0,0.10526315789474,0.895,0.105,0.2903772133122,0.06448949485448
0.15789473684211,0,0,0.15789473684211,0.845,0.155,0.19901105764177,0.11259819392439
0.21052631578947,0,0,0.21052631578947,0.785,0.215,0.12830647743282,0.17795538116674
0.26315789473684,0,0,0.26315789473684,0.735,0.265,0.088872460313874,0.23668776512578
0.31578947368421,0,0,0.31578947368421,0.685,0.315,0.061001763100158,0.29789909636912
0.36842105263158,0,0,0.36842105263158,0.635,0.365,0.041248157468635,0.36054389329128
0.42105263157895,0,0,0.42105263157895,0.575,0.425,0.02508173623612,0.43637851184429
0.47368421052632,0,0,0.47368421052632,0.525,0.475,0.016068838953146,0.49923931543403
0.52631578947368,0,0,0.52631578947368,0.475,0.525,0.0099302872835256,0.56108529784636
0.57894736842105,0,0,0.57894736842105,0.425,0.575,0.0058620826279179,0.62131825072922
0.63157894736842,0,0,0.63157894736842,0.365,0.635,0.0028772971835911,0.69070140912018
0.68421052631579,0,0,0.68421052631579,0.315,0.685,0.0014555937764247,0.74551858725947
0.73684210526316,0,0,0.73684210526316,0.265,0.735,0.00065862579390779,0.79710990847103
0.78947368421053,0,0,0.78947368421053,0.215,0.785,0.0002539594790149,0.84504001337092
0.84210526315789,0,0,0.84210526315789,0.155,0.845,5.7663441077274e-05,0.89714653900473
0.89473684210526,0,0,0.89473684210526,0.105,0.895,9.9607515954477e-06,0.93560608146977
0.94736842105263,0,0,0.94736842105263,0.055,0.945,5.4900521340933e-07,0.96916170874622
1,0,0,1,0.005,0.995,4.3064800892356e-11,0.99746260954264
@@ -3,20 +3,19 @@ x,y,z,id,s0aux,s1aux,kr0aux,kr1aux
0.052631578947368,0,0,0.052631578947368,0.945,0.055,1,0
0.10526315789474,0,0,0.10526315789474,0.895,0.105,1,0
0.15789473684211,0,0,0.15789473684211,0.845,0.155,1,0
0.21052631578947,0,0,0.21052631578947,0.785,0.215,1,0
0.26315789473684,0,0,0.26315789473684,0.735,0.265,1,0
0.31578947368421,0,0,0.31578947368421,0.685,0.315,0.56661038192749,4.5449082513583e-08
0.36842105263158,0,0,0.36842105263158,0.635,0.365,0.23996954587445,2.6372366439369e-05
0.42105263157895,0,0,0.42105263157895,0.575,0.425,0.099352541677511,0.00050637103510007
0.47368421052632,0,0,0.47368421052632,0.525,0.475,0.046506467919796,0.0023723093152781
0.52631578947368,0,0,0.52631578947368,0.475,0.525,0.020172644076885,0.0076876821741651
0.57894736842105,0,0,0.57894736842105,0.425,0.575,0.0076876821741651,0.020172644076885
0.63157894736842,0,0,0.63157894736842,0.365,0.635,0.0018074052303307,0.054380952711592
0.68421052631579,0,0,0.68421052631579,0.315,0.685,0.0003464682754216,0.11506193157389
0.73684210526316,0,0,0.73684210526316,0.265,0.735,2.6372366439369e-05,0.23996954587445
0.78947368421053,0,0,0.78947368421053,0.215,0.785,4.5449082513583e-08,0.56661038192749
0.84210526315789,0,0,0.84210526315789,0.155,0.845,0,1
0.89473684210526,0,0,0.89473684210526,0.105,0.895,0,1
0.21052631578947,0,0,0.21052631578947,0.785,0.215,0.47604805967375,0.013992735535288
0.26315789473684,0,0,0.26315789473684,0.735,0.265,0.20004065161298,0.089524470056521
0.31578947368421,0,0,0.31578947368421,0.685,0.315,0.10226261608964,0.17960810080408
0.36842105263158,0,0,0.36842105263158,0.635,0.365,0.05408762339136,0.27407835829833
0.42105263157895,0,0,0.42105263157895,0.575,0.425,0.02486349852954,0.387113860877
0.47368421052632,0,0,0.47368421052632,0.525,0.475,0.012511029810272,0.47804683204256
0.52631578947368,0,0,0.52631578947368,0.475,0.525,0.0059345605693258,0.56423831346722
0.57894736842105,0,0,0.57894736842105,0.425,0.575,0.0025855072999061,0.64456565920285
0.63157894736842,0,0,0.63157894736842,0.365,0.635,0.00080986520515468,0.73217663813379
0.68421052631579,0,0,0.68421052631579,0.315,0.685,0.00025084507084325,0.79736762215455
0.73684210526316,0,0,0.73684210526316,0.265,0.735,5.7679593442309e-05,0.85530492679017
0.78947368421053,0,0,0.78947368421053,0.215,0.785,7.8738051594707e-06,0.90615376154244
0.84210526315789,0,0,0.84210526315789,0.155,0.845,1.3903791841047e-07,0.95857273053566
0.89473684210526,0,0,0.89473684210526,0.105,0.895,1.5427817391226e-12,0.99641446294669
0.94736842105263,0,0,0.94736842105263,0.055,0.945,0,1
1,0,0,1,0.005,0.995,0,1

@@ -124,6 +124,7 @@
type = PorousFlowRelativePermeabilityVG
phase = 1
m = 0.5
wetting = false
[../]
[]

Oops, something went wrong.

0 comments on commit 027cbad

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