www.mooseframework.org
Classes | Functions
PorousFlowVanGenuchten Namespace Reference

van Genuchten effective saturation, capillary pressure and relative permeability functions. More...

Classes

struct  HighCapillaryPressureExtension
 Parameters associated with the extension of the hysteretic wetting capillary pressure function to high saturation values @ ExtensionStrategy the type of extension used @ S liquid saturation at the point of extension @ Pc capillary pressure at the point of extension @ dPc d(Pc)/dS at the point of extension. More...
 
struct  LowCapillaryPressureExtension
 Parameters associated with the extension of the hysteretic capillary pressure function to low saturation values @ ExtensionStrategy the type of extension used @ S liquid saturation at the point of extension @ Pc capillary pressure at the point of extension @ dPc d(Pc)/dS at the point of extension. More...
 

Functions

Real effectiveSaturation (Real p, Real alpha, Real m)
 Effective saturation as a function of porepressure. More...
 
Real dEffectiveSaturation (Real p, Real alpha, Real m)
 Derivative of effective saturation wrt porepressure. More...
 
Real d2EffectiveSaturation (Real p, Real alpha, Real m)
 Second derivative of effective saturation wrt porepressure. More...
 
Real capillaryPressure (Real seff, Real alpha, Real m, Real pc_max)
 Capillary pressure as a function of effective saturation. More...
 
Real dCapillaryPressure (Real seff, Real alpha, Real m, Real pc_max)
 Derivative of capillary pressure wrt effective saturation. More...
 
Real d2CapillaryPressure (Real seff, Real alpha, Real m, Real pc_max)
 Second derivative of capillary pressure wrt effective saturation. More...
 
template<typename T >
relativePermeability (const T &seff, Real m)
 Relative permeability as a function of effective saturation. More...
 
Real dRelativePermeability (Real seff, Real m)
 Derivative of relative permeability with respect to effective saturation. More...
 
Real d2RelativePermeability (Real seff, Real m)
 Second derivative of relative permeability with respect to effective saturation. More...
 
template<typename T >
relativePermeabilityNW (const T &seff, Real m)
 Relative permeability for a non-wetting phase as a function of effective saturation. More...
 
Real dRelativePermeabilityNW (Real seff, Real m)
 Derivative of relative permeability for a non-wetting phase with respect to effective saturation. More...
 
Real d2RelativePermeabilityNW (Real seff, Real m)
 Second derivative of relative permeability for a non-wetting phase with respect to effective saturation. More...
 
Real capillaryPressureHys (Real sl, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
 Hysteretic capillary pressure function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of Doughty2008). More...
 
Real dcapillaryPressureHys (Real sl, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
 Derivative of capillaryPressureHys with respect to sl. More...
 
Real d2capillaryPressureHys (Real sl, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
 Second derivative of capillaryPressureHys with respect to sl. More...
 
Real saturationHys (Real pc, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
 Hysteretic saturation function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of Doughty2008), which is the inverse of capillaryPressureHys. More...
 
Real dsaturationHys (Real pc, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
 Derivative of Hysteretic saturation function with respect to pc. More...
 
Real d2saturationHys (Real pc, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
 Second derivative of Hysteretic saturation function with respect to pc. More...
 
Real relativePermeabilityHys (Real sl, Real slr, Real sgrdel, Real sgrmax, Real sldel, Real m, Real upper_liquid_param, Real y0, Real y0p, Real y1, Real y1p)
 Hysteretic relative permeability for liquid. More...
 
Real drelativePermeabilityHys (Real sl, Real slr, Real sgrdel, Real sgrmax, Real sldel, Real m, Real upper_liquid_param, Real y0, Real y0p, Real y1, Real y1p)
 Derivative of Hysteretic relative permeability for liquid, with respect to liquid saturation. More...
 
Real relativePermeabilityNWHys (Real sl, Real slr, Real sgrdel, Real sgrmax, Real sldel, Real m, Real gamma, Real k_rg_max, Real y0p)
 Hysteretic relative permeability for gas. More...
 
Real drelativePermeabilityNWHys (Real sl, Real slr, Real sgrdel, Real sgrmax, Real sldel, Real m, Real gamma, Real k_rg_max, Real y0p)
 Derivative of hysteretic relative permeability for gas with respect to the liquid saturation. More...
 

Detailed Description

van Genuchten effective saturation, capillary pressure and relative permeability functions.

Note: effective saturation is provided as a function of porepressure, not capillary pressure. Note: capillary pressure and relative permeability are functions of effective saturation. The derivatives are therefore given wrt effective saturation. These derivatives must be multiplied by the derivative of effective saturation wrt the true saturation in objects using these relations.

Based on van Genuchten, M. Th., A closed for equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc., 44, 892-898 (1980).

Function Documentation

◆ capillaryPressure()

Real PorousFlowVanGenuchten::capillaryPressure ( Real  seff,
Real  alpha,
Real  m,
Real  pc_max 
)

Capillary pressure as a function of effective saturation.

Parameters
seffeffective saturation
alphavan Genuchten parameter
mvan Genuchten exponent
pc_maxmaximum capillary pressure (Pa)
Returns
capillary pressure (Pa)

Definition at line 63 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowCapillaryPressureVG::capillaryPressureCurve(), and TEST().

64 {
65  if (seff >= 1.0)
66  return 0.0;
67  else if (seff <= 0.0)
68  return pc_max;
69  else
70  {
71  Real a = std::pow(seff, -1.0 / m) - 1.0;
72  return std::min(std::pow(a, 1.0 - m) / alpha, pc_max);
73  }
74 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)

◆ capillaryPressureHys()

Real PorousFlowVanGenuchten::capillaryPressureHys ( Real  sl,
Real  slmin,
Real  sgrdel,
Real  alpha,
Real  n,
const LowCapillaryPressureExtension low_ext = LowCapillaryPressureExtension(),
const HighCapillaryPressureExtension high_ext = HighCapillaryPressureExtension() 
)

Hysteretic capillary pressure function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of Doughty2008).

NOTE: this function is undefined for sl < 0 and sl > 1, so you MUST ensure 0 <= sl <= 1 NOTE: this returns a non-negative quantity.

Parameters
slliquid saturation. 0 <= sl <= 1
slminvalue of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
sgrdelvalue of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <= 1
alphavan Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
nvan Genuchten n parameter. n > 1
low_extstrategy and parameters to use for the extension in the small-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)
high_extstrategy and parameters to use for the extension in the high-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)

Definition at line 180 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowHystereticCapillaryPressure::capillaryPressureQp(), PorousFlowHystereticCapillaryPressure::computeTurningPointInfo(), PorousFlowHystereticCapillaryPressure::firstOrderWettingPc(), PorousFlowHystereticCapillaryPressure::initQpStatefulProperties(), PorousFlowHystereticCapillaryPressure::secondOrderDryingPc(), and TEST().

187 {
188  Real pc = 0.0;
189  if (sl < low_ext.S) // important for initializing low_ext that this is < and not <=
190  {
191  switch (low_ext.strategy)
192  {
193  case LowCapillaryPressureExtension::QUADRATIC:
194  pc = low_ext.Pc + low_ext.dPc * 0.5 * (sl * sl - low_ext.S * low_ext.S) / low_ext.S;
195  break;
196  case LowCapillaryPressureExtension::EXPONENTIAL:
197  pc = low_ext.Pc * std::exp(low_ext.dPc * (sl - low_ext.S) / low_ext.Pc);
198  break;
199  default:
200  pc = low_ext.Pc;
201  }
202  return pc;
203  }
204  if (sl > high_ext.S) // important for initializing high_ext that this is >, not >=
205  {
206  switch (high_ext.strategy)
207  {
208  case HighCapillaryPressureExtension::POWER:
209  {
210  if (sl >= 1.0)
211  pc = 0.0;
212  else
213  {
214  const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
215  pc = high_ext.Pc * std::pow((1.0 - sl) / (1.0 - high_ext.S), expon);
216  }
217  break;
218  }
219  default:
220  pc = 0.0;
221  }
222  return pc;
223  }
224  const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
225  if (seff >= 1.0)
226  pc = 0.0; // no sensible high extension defined
227  else if (seff <= 0.0)
228  pc = low_ext.Pc; // no sensible low extension defined
229  else
230  {
231  const Real a = std::pow(seff, n / (1.0 - n)) - 1.0;
232  pc = (1.0 / alpha) * std::pow(a, 1.0 / n);
233  }
234  return pc;
235 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)

◆ d2CapillaryPressure()

Real PorousFlowVanGenuchten::d2CapillaryPressure ( Real  seff,
Real  alpha,
Real  m,
Real  pc_max 
)

Second derivative of capillary pressure wrt effective saturation.

Parameters
seffeffective saturation
alphavan Genuchten parameter
mvan Genuchten exponent
pc_maxmaximum capillary pressure (Pa)
Returns
second derivative of capillary pressure wrt effective saturation

Definition at line 93 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowCapillaryPressureVG::d2CapillaryPressureCurve(), and TEST().

94 {
95  if (seff <= 0.0 || seff >= 1.0)
96  return 0.0;
97  else
98  {
99  Real a = std::pow(seff, -1.0 / m) - 1.0;
100  // Return 0 if pc > pc_max
101  if (std::pow(a, 1.0 - m) / alpha > pc_max)
102  return 0.0;
103  else
104  {
105  Real d2pc = -std::pow(a, -1.0 - m) * std::pow(seff, -2.0 - 2.0 / m) +
106  ((1.0 + m) / m) * std::pow(a, -m) * std::pow(seff, -1.0 / m - 2.0);
107  d2pc *= (1.0 - m) / m / alpha;
108  return d2pc;
109  }
110  }
111 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)

◆ d2capillaryPressureHys()

Real PorousFlowVanGenuchten::d2capillaryPressureHys ( Real  sl,
Real  slmin,
Real  sgrdel,
Real  alpha,
Real  n,
const LowCapillaryPressureExtension low_ext = LowCapillaryPressureExtension(),
const HighCapillaryPressureExtension high_ext = HighCapillaryPressureExtension() 
)

Second derivative of capillaryPressureHys with respect to sl.

NOTE: this function is undefined for sl < 0 and sl > 1, so you MUST ensure 0 <= sl <= 1

Parameters
slliquid saturation. 0 <= sl <= 1
slminvalue of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
sgrdelvalue of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <= 1
alphavan Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
nvan Genuchten n parameter. n > 1
low_extstrategy and parameters to use for the extension in the small-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)
high_extstrategy and parameters to use for the extension in the high-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)

Definition at line 299 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowHystereticCapillaryPressure::d2capillaryPressureQp(), PorousFlowHystereticCapillaryPressure::d2firstOrderWettingPc(), PorousFlowHystereticCapillaryPressure::d2secondOrderDryingPc(), and TEST().

306 {
307  Real d2pc = 0.0;
308  if (sl < low_ext.S) // important for initializing low_ext that this is < and not <=
309  {
310  switch (low_ext.strategy)
311  {
312  case LowCapillaryPressureExtension::QUADRATIC:
313  d2pc = low_ext.dPc / low_ext.S;
314  break;
315  case LowCapillaryPressureExtension::EXPONENTIAL:
316  d2pc = std::pow(low_ext.dPc, 2) / low_ext.Pc *
317  std::exp(low_ext.dPc * (sl - low_ext.S) / low_ext.Pc);
318  break;
319  default:
320  d2pc = 0.0;
321  }
322  return d2pc;
323  }
324  if (sl > high_ext.S) // important for initializing high_ext that this is >, not >=
325  {
326  switch (high_ext.strategy)
327  {
328  case HighCapillaryPressureExtension::POWER:
329  {
330  if (sl >= 1.0)
331  d2pc = 0.0;
332  else
333  {
334  const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
335  d2pc = high_ext.dPc * (1.0 - expon) / (1.0 - high_ext.S) *
336  std::pow((1.0 - sl) / (1.0 - high_ext.S), expon - 2.0);
337  }
338  break;
339  }
340  default:
341  d2pc = 0.0;
342  }
343  return d2pc;
344  }
345  const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
346  if (seff >= 1.0)
347  d2pc = 0.0; // no sensible high extension defined
348  else if (seff <= 0.0)
349  d2pc = 0.0; // no sensible low extension defined
350  else
351  {
352  const Real a = std::pow(seff, n / (1.0 - n)) - 1.0;
353  const Real dseff = 1.0 / (1.0 - sgrdel - slmin);
354  const Real d2pc_dseff =
355  (1.0 / alpha / (1.0 - n)) *
356  (std::pow(a, 1.0 / n - 2.0) * std::pow(seff, 2.0 * (n / (1.0 - n) - 1.0)) +
357  (n / (1.0 - n) - 1.0) * std::pow(a, 1.0 / n - 1.0) * std::pow(seff, n / (1.0 - n) - 2.0));
358  d2pc = d2pc_dseff * dseff * dseff;
359  }
360  return d2pc;
361 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)

◆ d2EffectiveSaturation()

Real PorousFlowVanGenuchten::d2EffectiveSaturation ( Real  p,
Real  alpha,
Real  m 
)

Second derivative of effective saturation wrt porepressure.

Parameters
pporepressure
alphavan Genuchten parameter
mvan Genuchten exponent
Returns
second derivative of effective saturation wrt porepressure

Definition at line 46 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowCapillaryPressureVG::d2EffectiveSaturation(), and TEST().

47 {
48  if (p >= 0.0)
49  return 0.0;
50  else
51  {
52  Real n = 1.0 / (1.0 - m);
53  Real inner = 1.0 + std::pow(-alpha * p, n);
54  Real dinner_dp = -n * alpha * std::pow(-alpha * p, n - 1.0);
55  Real d2inner_dp2 = n * (n - 1.0) * alpha * alpha * std::pow(-alpha * p, n - 2.0);
56  Real d2seff_dp2 = m * (m + 1.0) * std::pow(inner, -m - 2.0) * std::pow(dinner_dp, 2.0) -
57  m * std::pow(inner, -m - 1.0) * d2inner_dp2;
58  return d2seff_dp2;
59  }
60 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)

◆ d2RelativePermeability()

Real PorousFlowVanGenuchten::d2RelativePermeability ( Real  seff,
Real  m 
)

Second derivative of relative permeability with respect to effective saturation.

Parameters
seffeffective saturation
mvan Genuchten exponent
Returns
second derivative of relative permeability wrt effective saturation

Definition at line 129 of file PorousFlowVanGenuchten.C.

Referenced by TEST().

130 {
131  // Guard against division by zero
132  if (seff <= 0.0 || seff >= 1.0)
133  return 0.0;
134 
135  const Real a = 1.0 - std::pow(seff, 1.0 / m);
136  const Real da = -1.0 / m * std::pow(seff, 1.0 / m - 1.0);
137  const Real d2a = -(1.0 / m) * (1.0 / m - 1.0) * std::pow(seff, 1.0 / m - 2.0);
138  const Real b = 1.0 - std::pow(a, m);
139  const Real db = -m * std::pow(a, m - 1.0) * da;
140  const Real d2b = -m * (m - 1.0) * std::pow(a, m - 2.0) * da * da - m * std::pow(a, m - 1.0) * d2a;
141 
142  return -0.25 * std::pow(seff, -1.5) * Utility::pow<2>(b) + 2.0 * std::pow(seff, -0.5) * b * db +
143  2.0 * std::sqrt(seff) * db * db + 2.0 * std::sqrt(seff) * b * d2b;
144 }
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
const GeochemicalDatabaseReader db("database/moose_testdb.json", true, true, false)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MooseUnits pow(const MooseUnits &, int)

◆ d2RelativePermeabilityNW()

Real PorousFlowVanGenuchten::d2RelativePermeabilityNW ( Real  seff,
Real  m 
)

Second derivative of relative permeability for a non-wetting phase with respect to effective saturation.

Parameters
seffeffective saturation
mvan Genuchten exponent
Returns
second derivative of relative permeability wrt effective saturation

Definition at line 162 of file PorousFlowVanGenuchten.C.

Referenced by TEST().

163 {
164  // Guard against division by zero
165  if (seff <= 0.0 || seff >= 1.0)
166  return 0.0;
167 
168  const Real a = std::pow(1.0 - seff, 1.0 / m);
169  const Real da = -1.0 / m * a / (1.0 - seff);
170  const Real d2a = 1.0 / m * (1.0 / m - 1) * std::pow(1.0 - seff, 1.0 / m - 2.0);
171  const Real b = std::pow(1.0 - a, 2.0 * m);
172  const Real db = -2.0 * m * b / (1.0 - a) * da;
173  const Real d2b =
174  -2.0 * m * (db / (1.0 - a) * da + b * Utility::pow<2>(da / (1.0 - a)) + b / (1.0 - a) * d2a);
175 
176  return -0.25 * std::pow(seff, -1.5) * b + std::pow(seff, -0.5) * db + std::sqrt(seff) * d2b;
177 }
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
const GeochemicalDatabaseReader db("database/moose_testdb.json", true, true, false)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MooseUnits pow(const MooseUnits &, int)

◆ d2saturationHys()

Real PorousFlowVanGenuchten::d2saturationHys ( Real  pc,
Real  slmin,
Real  sgrdel,
Real  alpha,
Real  n,
const LowCapillaryPressureExtension low_ext = LowCapillaryPressureExtension(),
const HighCapillaryPressureExtension high_ext = HighCapillaryPressureExtension() 
)

Second derivative of Hysteretic saturation function with respect to pc.

Parameters
pccapillary pressure. 0 <= pc
slminvalue of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
sgrdelvalue of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <= 1
alphavan Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
nvan Genuchten n parameter. n > 1
low_extstrategy and parameters to use for the extension in the small-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)
high_extstrategy and parameters to use for the extension in the high-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)

Definition at line 507 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowHystereticCapillaryPressure::d2firstOrderWettingSat(), PorousFlowHystereticCapillaryPressure::d2liquidSaturationQp(), PorousFlowHystereticCapillaryPressure::d2secondOrderDryingSat(), and TEST().

514 {
515  if (pc <= 0)
516  return 0.0;
517  Real d2s = 0.0;
518  if (pc > low_ext.Pc) // important for initialization of the low_ext that this is > and not >=
519  {
520  switch (low_ext.strategy)
521  {
522  case LowCapillaryPressureExtension::QUADRATIC:
523  {
524  const Real s2 = low_ext.S * low_ext.S + 2.0 * (pc - low_ext.Pc) * low_ext.S / low_ext.dPc;
525  if (s2 <= 0.0) // this occurs when we're trying to find a saturation on the wetting curve
526  // defined by sl = sgrDel at pc = Pcd_Del, if this pc is actually impossible
527  // to achieve on this wetting curve
528  d2s = 0.0;
529  else
530  {
531  const Real ds2 = 2.0 * low_ext.S / low_ext.dPc;
532  d2s = -0.25 * ds2 * ds2 / std::pow(s2, 1.5);
533  }
534  break;
535  }
536  case LowCapillaryPressureExtension::EXPONENTIAL:
537  {
538  const Real s = low_ext.S + std::log(pc / low_ext.Pc) * low_ext.Pc / low_ext.dPc;
539  if (s <= 0.0) // this occurs when we're trying to find a saturation on the
540  // wetting curve defined by sl = sgrDel at pc = Pcd_Del, if this
541  // pc is actually impossible to achieve on this wetting curve
542  d2s = 0.0;
543  else
544  d2s = -low_ext.Pc / std::pow(pc, 2.0) / low_ext.dPc;
545  break;
546  }
547  default:
548  d2s = 0.0;
549  }
550  return d2s;
551  }
552  if (pc < high_ext.Pc) // important for initialization of the high_ext that this is < and not <=
553  {
554  switch (high_ext.strategy)
555  {
556  case HighCapillaryPressureExtension::POWER:
557  {
558  const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
559  d2s = -(1.0 - high_ext.S) * (1.0 / expon) * (1.0 / expon - 1.0) /
560  std::pow(high_ext.Pc, 2.0) * std::pow(pc / high_ext.Pc, 1.0 / expon - 2.0);
561  break;
562  }
563  default:
564  d2s = 0.0;
565  }
566  return d2s;
567  }
568  if (pc == std::numeric_limits<Real>::max())
569  d2s = 0.0;
570  else
571  {
572  const Real seffpow = 1.0 + std::pow(pc * alpha, n);
573  const Real dseffpow = n * (seffpow - 1.0) / pc;
574  const Real d2seffpow = (n - 1.0) * dseffpow / pc;
575  const Real seff = std::pow(seffpow, (1.0 - n) / n);
576  const Real dseff = (1.0 - n) / n * seff / seffpow * dseffpow;
577  const Real d2seff =
578  (1.0 - n) / n *
579  (dseff * dseffpow - seff * dseffpow * dseffpow / seffpow + seff * d2seffpow) / seffpow;
580  d2s = (1.0 - sgrdel - slmin) * d2seff;
581  }
582  return d2s;
583 }
if(subdm)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)

◆ dCapillaryPressure()

Real PorousFlowVanGenuchten::dCapillaryPressure ( Real  seff,
Real  alpha,
Real  m,
Real  pc_max 
)

Derivative of capillary pressure wrt effective saturation.

Parameters
seffeffective saturation
alphavan Genuchten parameter
mvan Genuchten exponent
pc_maxmaximum capillary pressure (Pa)
Returns
derivative of capillary pressure wrt effective saturation

Definition at line 77 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowCapillaryPressureVG::dCapillaryPressureCurve(), and TEST().

78 {
79  if (seff <= 0.0 || seff >= 1.0)
80  return 0.0;
81  else
82  {
83  Real a = std::pow(seff, -1.0 / m) - 1.0;
84  // Return 0 if pc > pc_max
85  if (std::pow(a, 1.0 - m) / alpha > pc_max)
86  return 0.0;
87  else
88  return (m - 1.0) * std::pow(a, -m) * std::pow(seff, -1.0 - 1.0 / m) / m / alpha;
89  }
90 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)

◆ dcapillaryPressureHys()

Real PorousFlowVanGenuchten::dcapillaryPressureHys ( Real  sl,
Real  slmin,
Real  sgrdel,
Real  alpha,
Real  n,
const LowCapillaryPressureExtension low_ext = LowCapillaryPressureExtension(),
const HighCapillaryPressureExtension high_ext = HighCapillaryPressureExtension() 
)

Derivative of capillaryPressureHys with respect to sl.

NOTE: this function is undefined for sl < 0 and sl > 1, so you MUST ensure 0 <= sl <= 1 NOTE: this returns a negative quantity.

Parameters
slliquid saturation. 0 <= sl <= 1
slminvalue of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
sgrdelvalue of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <= 1
alphavan Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
nvan Genuchten n parameter. n > 1
low_extstrategy and parameters to use for the extension in the small-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)
high_extstrategy and parameters to use for the extension in the high-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)

Definition at line 238 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowHystereticCapillaryPressure::computeTurningPointInfo(), PorousFlowHystereticCapillaryPressure::dcapillaryPressureQp(), PorousFlowHystereticCapillaryPressure::dfirstOrderWettingPc(), PorousFlowHystereticCapillaryPressure::dsecondOrderDryingPc(), and TEST().

245 {
246  Real dpc = 0.0;
247  if (sl < low_ext.S) // important for initializing low_ext that this is < and not <=
248  {
249  switch (low_ext.strategy)
250  {
251  case LowCapillaryPressureExtension::QUADRATIC:
252  dpc = low_ext.dPc * sl / low_ext.S;
253  break;
254  case LowCapillaryPressureExtension::EXPONENTIAL:
255  dpc = low_ext.dPc * std::exp(low_ext.dPc * (sl - low_ext.S) / low_ext.Pc);
256  break;
257  default:
258  dpc = 0.0;
259  }
260  return dpc;
261  }
262  if (sl > high_ext.S) // important for initializing high_ext that this is >, not >=
263  {
264  switch (high_ext.strategy)
265  {
266  case HighCapillaryPressureExtension::POWER:
267  {
268  if (sl >= 1.0)
269  dpc = 0.0;
270  else
271  {
272  const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
273  dpc = high_ext.dPc * std::pow((1.0 - sl) / (1.0 - high_ext.S), expon - 1.0);
274  }
275  break;
276  }
277  default:
278  dpc = 0.0;
279  }
280  return dpc;
281  }
282  const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
283  if (seff >= 1.0)
284  dpc = 0.0; // no sensible high extension defined
285  else if (seff <= 0.0)
286  dpc = 0.0; // no sensible low extension defined
287  else
288  {
289  const Real a = std::pow(seff, n / (1.0 - n)) - 1.0;
290  const Real dseff = 1.0 / (1.0 - sgrdel - slmin);
291  const Real dpc_dseff = (1.0 / alpha / (1.0 - n)) * std::pow(a, 1.0 / n - 1.0) *
292  std::pow(seff, n / (1.0 - n) - 1.0);
293  dpc = dpc_dseff * dseff;
294  }
295  return dpc;
296 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)

◆ dEffectiveSaturation()

Real PorousFlowVanGenuchten::dEffectiveSaturation ( Real  p,
Real  alpha,
Real  m 
)

Derivative of effective saturation wrt porepressure.

Parameters
pporepressure
alphavan Genuchten parameter
mvan Genuchten exponent
Returns
derivative of effective saturation wrt porepressure

Definition at line 31 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowCapillaryPressureVG::dEffectiveSaturation(), and TEST().

32 {
33  if (p >= 0.0)
34  return 0.0;
35  else
36  {
37  Real n = 1.0 / (1.0 - m);
38  Real inner = 1.0 + std::pow(-alpha * p, n);
39  Real dinner_dp = -n * alpha * std::pow(-alpha * p, n - 1.0);
40  Real dseff_dp = -m * std::pow(inner, -m - 1) * dinner_dp;
41  return dseff_dp;
42  }
43 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)

◆ dRelativePermeability()

Real PorousFlowVanGenuchten::dRelativePermeability ( Real  seff,
Real  m 
)

Derivative of relative permeability with respect to effective saturation.

Parameters
seffeffective saturation
mvan Genuchten exponent
Returns
derivative of relative permeability wrt effective saturation

Definition at line 114 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowRelativePermeabilityVGTempl< is_ad >::dRelativePermeability(), and TEST().

115 {
116  // Guard against division by zero
117  if (seff <= 0.0 || seff >= 1.0)
118  return 0.0;
119 
120  const Real a = 1.0 - std::pow(seff, 1.0 / m);
121  const Real da = -1.0 / m * std::pow(seff, 1.0 / m - 1.0);
122  const Real b = 1.0 - std::pow(a, m);
123  const Real db = -m * std::pow(a, m - 1.0) * da;
124 
125  return 0.5 * std::pow(seff, -0.5) * Utility::pow<2>(b) + 2.0 * std::sqrt(seff) * b * db;
126 }
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
const GeochemicalDatabaseReader db("database/moose_testdb.json", true, true, false)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MooseUnits pow(const MooseUnits &, int)

◆ drelativePermeabilityHys()

Real PorousFlowVanGenuchten::drelativePermeabilityHys ( Real  sl,
Real  slr,
Real  sgrdel,
Real  sgrmax,
Real  sldel,
Real  m,
Real  upper_liquid_param,
Real  y0,
Real  y0p,
Real  y1,
Real  y1p 
)

Derivative of Hysteretic relative permeability for liquid, with respect to liquid saturation.

Parameters
slliquid saturation
slrresidual liquid saturation. For sl < slr, this function will always return 0
sgrdelvalue of gas saturation where van Genuchten wetting capillary-pressure expression -> 0. This depends on the turning-point saturation when drying became wetting, using the Land equation
sgrmaxmaximum value possible for sgrdel. This will be equal to sgrdel if the turning-point saturation is small
sldelvalue of the turning-point saturation when drying became wetting
mvan-Genuchten m parameter
upper_liquid_paramcubic modification of the wetting relative permeability will occur between upper_liquid_param * (1 - sgrdel) and 1 - 0.5 * sgrdel. 0 < upper_liquid_param <= 1. Usually upper_liquid_param is close to 1 (eg 0.9)
y0value of the unmodified wetting relative permeability at sl = upper_liquid_param * (1 - sgrdel)
y0pvalue of the derivtaive of the unmodified wetting relative permeability at sl = upper_liquid_param * (1 - sgrdel)
y1value of the unmodified wetting relative permeability at sl = 1 - 0.5 * sgrdel
y1pvalue of the derivtaive of the unmodified wetting relative permeability at sl = 1 - 0.5 * sgrdel

Definition at line 642 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowHystereticRelativePermeabilityLiquid::computeRelPermQp(), PorousFlowHystereticRelativePermeabilityLiquid::computeTurningPoint0Info(), and TEST().

653 {
654  if (sl <= slr) // by the definition of slr, always return 0
655  return 0.0;
656  if (sl == 1.0) // derivative is infinite at this point
657  return std::numeric_limits<Real>::max();
658  const Real sl_bar = (sl - slr) / (1.0 - slr); // effective saturation
659  const Real sl_bar_prime = 1.0 / (1.0 - slr);
660  // a and b are useful parameters. Define b along the drying curve initially, and
661  // modify a and b appropriately if the wetting result is required
662  Real a = 0;
663  Real a_prime = 0.0;
664  Real b = 0;
665  Real b_prime = 0.0;
666  if (sgrdel == 0.0 || sl <= sldel) // along the drying curve
667  {
668  const Real c = std::pow(sl_bar, 1.0 / m);
669  const Real dc_dsbar = c / m / sl_bar;
670  a = std::pow(1.0 - c, m);
671  const Real da_dsbar = -m * a / (1.0 - c) * dc_dsbar;
672  a_prime = da_dsbar * sl_bar_prime;
673  }
674  else // along the wetting curve
675  {
676  // In most cases, use sldel and sgrdel as provided to this function. However, because "there is
677  // no hysteresis along the extension" according to p6 of Doughty2008, if the turning point is
678  // less than slr, then use the expressions for the case when the turning point was slr
679  const Real my_sldel = (sldel < slr) ? slr : sldel;
680  const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
681  if (sl >= 1.0 - 0.5 * my_sgrdel)
682  {
683  // follow the drying curve. The parameter b has already been defined. It
684  // is important for initialization of the curic that the above condition is >= and not >
685  const Real c = std::pow(sl_bar, 1.0 / m);
686  const Real dc_dsbar = c / m / sl_bar;
687  a = std::pow(1.0 - c, m);
688  const Real da_dsbar = -m * a / (1.0 - c) * dc_dsbar;
689  a_prime = da_dsbar * sl_bar_prime;
690  }
691  else if (sl > upper_liquid_param * (1.0 - my_sgrdel))
692  {
693  // follow the cubic modification of the wetting curve. Immediately exit from this function by
694  // returning the cubic result
696  sl, upper_liquid_param * (1.0 - my_sgrdel), y0, y0p, 1.0 - 0.5 * my_sgrdel, y1, y1p);
697  }
698  else
699  {
700  // standard case of wetting curve outside the cubic-modification and drying-curve regions
701  const Real sl_bar_del = (my_sldel - slr) / (1.0 - slr);
702  const Real s_gt_bar =
703  my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
704  const Real s_gt_bar_prime = my_sgrdel / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
705  const Real c = std::pow(sl_bar + s_gt_bar, 1.0 / m);
706  const Real c_prime = c / m / (sl_bar + s_gt_bar) * (sl_bar_prime + s_gt_bar_prime);
707  a = (1 - s_gt_bar / (1.0 - sl_bar_del)) * std::pow(1.0 - c, m);
708  a_prime =
709  -s_gt_bar_prime / (1.0 - sl_bar_del) * std::pow(1.0 - c, m) - m * a / (1.0 - c) * c_prime;
710  b = s_gt_bar / (1.0 - sl_bar_del) * std::pow(1.0 - std::pow(sl_bar_del, 1.0 / m), m);
711  b_prime = s_gt_bar_prime * b / s_gt_bar;
712  }
713  }
714  const Real kr = std::sqrt(sl_bar) * Utility::pow<2>(1.0 - a - b);
715  return 0.5 * kr / sl_bar * sl_bar_prime -
716  std::sqrt(sl_bar) * 2.0 * (1.0 - a - b) * (a_prime + b_prime);
717 }
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
Real dcubic(Real x, Real x0, Real y0, Real y0p, Real x1, Real y1, Real y1p)
Derivative of cubic function, f(x), with respect to x.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MooseUnits pow(const MooseUnits &, int)

◆ dRelativePermeabilityNW()

Real PorousFlowVanGenuchten::dRelativePermeabilityNW ( Real  seff,
Real  m 
)

Derivative of relative permeability for a non-wetting phase with respect to effective saturation.

Parameters
seffeffective saturation
mvan Genuchten exponent
Returns
derivative of relative permeability wrt effective saturation

Definition at line 147 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowRelativePermeabilityVGTempl< is_ad >::dRelativePermeability(), and TEST().

148 {
149  // Guard against division by zero
150  if (seff <= 0.0 || seff >= 1.0)
151  return 0.0;
152 
153  const Real a = std::pow(1.0 - seff, 1.0 / m);
154  const Real da = -1.0 / m * a / (1.0 - seff);
155  const Real b = std::pow(1.0 - a, 2.0 * m);
156  const Real db = -2.0 * m * b / (1.0 - a) * da;
157 
158  return 0.5 * std::pow(seff, -0.5) * b + std::sqrt(seff) * db;
159 }
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
const GeochemicalDatabaseReader db("database/moose_testdb.json", true, true, false)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MooseUnits pow(const MooseUnits &, int)

◆ drelativePermeabilityNWHys()

Real PorousFlowVanGenuchten::drelativePermeabilityNWHys ( Real  sl,
Real  slr,
Real  sgrdel,
Real  sgrmax,
Real  sldel,
Real  m,
Real  gamma,
Real  k_rg_max,
Real  y0p 
)

Derivative of hysteretic relative permeability for gas with respect to the liquid saturation.

Parameters
slliquid saturation
slrresidual liquid saturation. For sl < slr, this function will always return 0
sgrdelvalue of gas saturation where van Genuchten wetting capillary-pressure expression -> 0. This depends on the turning-point saturation when drying became wetting, using the Land equation
sgrmaxmaximum value possible for sgrdel. This will be equal to sgrdel if the turning-point saturation is small
sldelvalue of the turning-point saturation when drying became wetting
mvan-Genuchten m parameter
gammaindex satisfying gamma > 0. Usually gamma = 1/3.
k_rg_maxMaximum value of unextended gas relative permeability. If no low-saturation extension is used then gas relative permeability = k_rg_max for sl <= slr. Satisfies 0 < k_rg_max <= 1. Frequently k_rg_max = 1 is used
y0pValue of the derivative of the low-saturation extension at sl = slr. If an extension is used then the gas relative permeability in the region sl <= slr is a cubic whose value is unity at sl = 0, and derivative is zero at sl = 0

Definition at line 764 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowHystereticRelativePermeabilityGas::computeRelPermQp(), and TEST().

773 {
774  if (sl < slr)
775  {
776  // in the extended region, so immediately return with the relevant value
777  if (k_rg_max == 1.0)
778  return 0.0;
779  return PorousFlowCubic::dcubic(sl, 0.0, 1.0, 0.0, slr, k_rg_max, y0p);
780  }
781  if (sl > 1.0 - sgrdel) // saturation is above 1.0 - residual gas saturation
782  return 0.0;
783  const Real sl_bar = (sl - slr) / (1.0 - slr);
784  const Real sl_bar_prime = 1.0 / (1.0 - slr);
785  Real s_gt_bar = 0.0; // initialize this parameter as if on the drying curve
786  Real s_gt_bar_prime = 0.0; // again, assume on drying curve
787  if (sgrdel != 0.0 && sl > sldel)
788  {
789  // On the wetting curve
790  // In most cases, use sldel and sgrdel as provided to this function. However, because "there is
791  // no hysteresis along the extension" according to p6 of Doughty2008, if the turning point is
792  // less than slr, then use the expressions for the case when the turning point was slr
793  const Real my_sldel = (sldel < slr) ? slr : sldel;
794  const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
795  s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
796  s_gt_bar_prime = my_sgrdel / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
797  }
798  Real kr_prime = 0.0;
799  if (sl_bar + s_gt_bar < 1.0) // check for the condition where sl is too big, in which case kr =
800  // 0, irrespective of hysteresis
801  {
802  const Real a = std::pow(1.0 - (sl_bar + s_gt_bar), gamma);
803  const Real a_prime = -gamma * a / (1.0 - (sl_bar + s_gt_bar)) * (sl_bar_prime + s_gt_bar_prime);
804  const Real c = std::pow(sl_bar + s_gt_bar, 1.0 / m);
805  const Real c_prime =
806  (c == 0 ? 0.0 : c / m / (sl_bar + s_gt_bar) * (sl_bar_prime + s_gt_bar_prime));
807  const Real b = std::pow(1.0 - c, 2.0 * m);
808  const Real b_prime = -2.0 * m * b / (1.0 - c) * c_prime;
809  kr_prime = k_rg_max * (a * b_prime + a_prime * b);
810  }
811  return kr_prime;
812 }
Real dcubic(Real x, Real x0, Real y0, Real y0p, Real x1, Real y1, Real y1p)
Derivative of cubic function, f(x), with respect to x.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MooseUnits pow(const MooseUnits &, int)

◆ dsaturationHys()

Real PorousFlowVanGenuchten::dsaturationHys ( Real  pc,
Real  slmin,
Real  sgrdel,
Real  alpha,
Real  n,
const LowCapillaryPressureExtension low_ext = LowCapillaryPressureExtension(),
const HighCapillaryPressureExtension high_ext = HighCapillaryPressureExtension() 
)

Derivative of Hysteretic saturation function with respect to pc.

Parameters
pccapillary pressure. 0 <= pc
slminvalue of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
sgrdelvalue of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <= 1
alphavan Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
nvan Genuchten n parameter. n > 1
low_extstrategy and parameters to use for the extension in the small-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)
high_extstrategy and parameters to use for the extension in the high-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)

Definition at line 433 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowHystereticCapillaryPressure::dfirstOrderWettingSat(), PorousFlowHystereticCapillaryPressure::dliquidSaturationQp(), PorousFlowHystereticCapillaryPressure::dsecondOrderDryingSat(), and TEST().

440 {
441  if (pc <= 0)
442  return 0.0;
443  Real ds = 0.0;
444  if (pc > low_ext.Pc) // important for initialization of the low_ext that this is > and not >=
445  {
446  switch (low_ext.strategy)
447  {
448  case LowCapillaryPressureExtension::QUADRATIC:
449  {
450  const Real s2 = low_ext.S * low_ext.S + 2.0 * (pc - low_ext.Pc) * low_ext.S / low_ext.dPc;
451  if (s2 <= 0.0) // this occurs when we're trying to find a saturation on the wetting curve
452  // defined by sl = sgrDel at pc = Pcd_Del, if this pc is actually impossible
453  // to achieve on this wetting curve
454  ds = 0.0;
455  else
456  {
457  const Real ds2 = 2.0 * low_ext.S / low_ext.dPc;
458  ds = 0.5 * ds2 / std::sqrt(s2);
459  }
460  break;
461  }
462  case LowCapillaryPressureExtension::EXPONENTIAL:
463  {
464  const Real s = low_ext.S + std::log(pc / low_ext.Pc) * low_ext.Pc / low_ext.dPc;
465  if (s <= 0.0) // this occurs when we're trying to find a saturation on the
466  // wetting curve defined by sl = sgrDel at pc = Pcd_Del, if this
467  // pc is actually impossible to achieve on this wetting curve
468  ds = 0.0;
469  else
470  ds = low_ext.Pc / pc / low_ext.dPc;
471  break;
472  }
473  default:
474  ds = 0.0;
475  }
476  return ds;
477  }
478  if (pc < high_ext.Pc) // important for initialization of the high_ext that this is < and not <=
479  {
480  switch (high_ext.strategy)
481  {
482  case HighCapillaryPressureExtension::POWER:
483  {
484  const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
485  ds = -(1.0 - high_ext.S) / pc / expon * std::pow(pc / high_ext.Pc, 1.0 / expon);
486  break;
487  }
488  default:
489  ds = 0.0;
490  }
491  return ds;
492  }
493  if (pc == std::numeric_limits<Real>::max())
494  ds = 0.0;
495  else
496  {
497  const Real seffpow = 1.0 + std::pow(pc * alpha, n);
498  const Real dseffpow = n * (seffpow - 1.0) / pc;
499  const Real seff = std::pow(seffpow, (1.0 - n) / n);
500  const Real dseff = (1.0 - n) / n * seff / seffpow * dseffpow;
501  ds = (1.0 - sgrdel - slmin) * dseff;
502  }
503  return ds;
504 }
if(subdm)
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)

◆ effectiveSaturation()

Real PorousFlowVanGenuchten::effectiveSaturation ( Real  p,
Real  alpha,
Real  m 
)

Effective saturation as a function of porepressure.

Note: seff = 1 for p >= 0

Parameters
pporepressure
alphavan Genuchten parameter
mvan Genuchten exponent
Returns
effective saturation

Definition at line 16 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowCapillaryPressureVG::effectiveSaturation(), and TEST().

17 {
18  Real n, seff;
19 
20  if (p >= 0.0)
21  return 1.0;
22  else
23  {
24  n = 1.0 / (1.0 - m);
25  seff = 1.0 + std::pow(-alpha * p, n);
26  return std::pow(seff, -m);
27  }
28 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)

◆ relativePermeability()

template<typename T >
T PorousFlowVanGenuchten::relativePermeability ( const T &  seff,
Real  m 
)

Relative permeability as a function of effective saturation.

Parameters
seffeffective saturation
mvan Genuchten exponent
Returns
relative permeability

Definition at line 101 of file PorousFlowVanGenuchten.h.

Referenced by PorousFlowRelativePermeabilityVGTempl< is_ad >::relativePermeability(), and TEST().

102 {
103  if (MetaPhysicL::raw_value(seff) <= 0.0)
104  return 0.0;
105  else if (MetaPhysicL::raw_value(seff) >= 1.0)
106  return 1.0;
107 
108  const T a = 1.0 - std::pow(seff, 1.0 / m);
109  const T b = 1.0 - std::pow(a, m);
110 
111  return std::sqrt(seff) * Utility::pow<2>(b);
112 }
auto raw_value(const Eigen::Map< T > &in)
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
MooseUnits pow(const MooseUnits &, int)

◆ relativePermeabilityHys()

Real PorousFlowVanGenuchten::relativePermeabilityHys ( Real  sl,
Real  slr,
Real  sgrdel,
Real  sgrmax,
Real  sldel,
Real  m,
Real  upper_liquid_param,
Real  y0,
Real  y0p,
Real  y1,
Real  y1p 
)

Hysteretic relative permeability for liquid.

Parameters
slliquid saturation
slrresidual liquid saturation. For sl < slr, this function will always return 0
sgrdelvalue of gas saturation where van Genuchten wetting capillary-pressure expression -> 0. This depends on the turning-point saturation when drying became wetting, using the Land equation
sgrmaxmaximum value possible for sgrdel. This will be equal to sgrdel if the turning-point saturation is small
sldelvalue of the turning-point saturation when drying became wetting
mvan-Genuchten m parameter
upper_liquid_paramcubic modification of the wetting relative permeability will occur between upper_liquid_param * (1 - sgrdel) and 1 - 0.5 * sgrdel. 0 < upper_liquid_param <= 1. Usually upper_liquid_param is close to 1 (eg 0.9)
y0value of the unmodified wetting relative permeability at sl = upper_liquid_param * (1 - sgrdel)
y0pvalue of the derivtaive of the unmodified wetting relative permeability at sl = upper_liquid_param * (1 - sgrdel)
y1value of the unmodified wetting relative permeability at sl = 1 - 0.5 * sgrdel
y1pvalue of the derivtaive of the unmodified wetting relative permeability at sl = 1 - 0.5 * sgrdel

Definition at line 586 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowHystereticRelativePermeabilityLiquid::computeRelPermQp(), PorousFlowHystereticRelativePermeabilityLiquid::computeTurningPoint0Info(), and TEST().

597 {
598  if (sl <= slr) // by the definition of slr, always return 0
599  return 0.0;
600  const Real sl_bar = (sl - slr) / (1.0 - slr); // effective saturation
601  // a and b are useful parameters. Define b along the drying curve initially, and
602  // modify a and b appropriately if the wetting result is required
603  Real a = 0;
604  Real b = 0;
605  if (sgrdel == 0.0 || sl <= sldel) // along the drying curve
606  a = std::pow(1.0 - std::pow(sl_bar, 1.0 / m), m);
607  else // along the wetting curve
608  {
609  // In most cases, use sldel and sgrdel as provided to this function. However, because "there is
610  // no hysteresis along the extension" according to p6 of Doughty2008, if the turning point is
611  // less than slr, then use the expressions for the case when the turning point was slr
612  const Real my_sldel = (sldel < slr) ? slr : sldel;
613  const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
614  if (sl >= 1.0 - 0.5 * my_sgrdel)
615  {
616  // follow the drying curve. The parameter b has already been defined. It
617  // is important for initialization of the curic that the above condition is >= and not >
618  a = std::pow(1.0 - std::pow(sl_bar, 1.0 / m), m);
619  }
620  else if (sl > upper_liquid_param * (1.0 - my_sgrdel))
621  {
622  // follow the cubic modification of the wetting curve. Immediately exit from this function by
623  // returning the cubic result
624  return PorousFlowCubic::cubic(
625  sl, upper_liquid_param * (1.0 - my_sgrdel), y0, y0p, 1.0 - 0.5 * my_sgrdel, y1, y1p);
626  }
627  else
628  {
629  // standard case of wetting curve outside the cubic-modification and drying-curve regions
630  const Real sl_bar_del = (my_sldel - slr) / (1.0 - slr);
631  const Real s_gt_bar =
632  my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
633  a = (1 - s_gt_bar / (1.0 - sl_bar_del)) *
634  std::pow(1.0 - std::pow(sl_bar + s_gt_bar, 1.0 / m), m);
635  b = s_gt_bar / (1.0 - sl_bar_del) * std::pow(1.0 - std::pow(sl_bar_del, 1.0 / m), m);
636  }
637  }
638  return std::sqrt(sl_bar) * Utility::pow<2>(1.0 - a - b);
639 }
Real cubic(Real x, Real x0, Real y0, Real y0p, Real x1, Real y1, Real y1p)
Cubic function f(x) that satisfies f(x0) = y0 f&#39;(x0) = y0p f(x1) = y1 f&#39;(x1) = y1p.
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MooseUnits pow(const MooseUnits &, int)

◆ relativePermeabilityNW()

template<typename T >
T PorousFlowVanGenuchten::relativePermeabilityNW ( const T &  seff,
Real  m 
)

Relative permeability for a non-wetting phase as a function of effective saturation.

Parameters
seffeffective saturation
mvan Genuchten exponent
Returns
relative permeability

Definition at line 138 of file PorousFlowVanGenuchten.h.

Referenced by PorousFlowRelativePermeabilityVGTempl< is_ad >::relativePermeability(), and TEST().

139 {
140  if (MetaPhysicL::raw_value(seff) <= 0.0)
141  return 0.0;
142  else if (MetaPhysicL::raw_value(seff) >= 1.0)
143  return 1.0;
144 
145  const T a = std::pow(1.0 - seff, 1.0 / m);
146  const T b = std::pow(1.0 - a, 2.0 * m);
147 
148  return std::sqrt(seff) * b;
149 }
auto raw_value(const Eigen::Map< T > &in)
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
MooseUnits pow(const MooseUnits &, int)

◆ relativePermeabilityNWHys()

Real PorousFlowVanGenuchten::relativePermeabilityNWHys ( Real  sl,
Real  slr,
Real  sgrdel,
Real  sgrmax,
Real  sldel,
Real  m,
Real  gamma,
Real  k_rg_max,
Real  y0p 
)

Hysteretic relative permeability for gas.

Parameters
slliquid saturation
slrresidual liquid saturation. For sl < slr, this function will always return 0
sgrdelvalue of gas saturation where van Genuchten wetting capillary-pressure expression -> 0. This depends on the turning-point saturation when drying became wetting, using the Land equation
sgrmaxmaximum value possible for sgrdel. This will be equal to sgrdel if the turning-point saturation is small
sldelvalue of the turning-point saturation when drying became wetting
mvan-Genuchten m parameter
gammaindex satisfying gamma > 0. Usually gamma = 1/3.
k_rg_maxMaximum value of unextended gas relative permeability. If no low-saturation extension is used then gas relative permeability = k_rg_max for sl <= slr. Satisfies 0 < k_rg_max <= 1. Frequently k_rg_max = 1 is used
y0pValue of the derivative of the low-saturation extension at sl = slr. If an extension is used then the gas relative permeability in the region sl <= slr is a cubic whose value is unity at sl = 0, and derivative is zero at sl = 0

Definition at line 720 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowHystereticRelativePermeabilityGas::computeRelPermQp(), and TEST().

729 {
730  if (sl < slr)
731  {
732  // in the extended region, so immediately return with the relevant value
733  if (k_rg_max == 1.0)
734  return 1.0;
735  return PorousFlowCubic::cubic(sl, 0.0, 1.0, 0.0, slr, k_rg_max, y0p);
736  }
737  if (sl > 1.0 - sgrdel) // saturation is above 1.0 - residual gas saturation
738  return 0.0;
739  const Real sl_bar = (sl - slr) / (1.0 - slr);
740  Real s_gt_bar = 0.0; // initialize this parameter as if on the drying curve
741  if (sgrdel != 0.0 && sl > sldel)
742  {
743  // On the wetting curve
744  // In most cases, use sldel and sgrdel as provided to this function. However, because "there is
745  // no hysteresis along the extension" according to p6 of Doughty2008, if the turning point is
746  // less than slr, then use the expressions for the case when the turning point was slr
747  const Real my_sldel = (sldel < slr) ? slr : sldel;
748  const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
749  s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
750  }
751  Real kr = 0.0;
752  if (sl_bar + s_gt_bar < 1.0) // check for the condition where sl is too big, in which case kr =
753  // 0, irrespective of hysteresis
754  {
755  const Real a = std::pow(1.0 - (sl_bar + s_gt_bar), gamma);
756  const Real c = std::pow(sl_bar + s_gt_bar, 1.0 / m);
757  const Real b = std::pow(1.0 - c, 2.0 * m);
758  kr = k_rg_max * a * b;
759  }
760  return kr;
761 }
Real cubic(Real x, Real x0, Real y0, Real y0p, Real x1, Real y1, Real y1p)
Cubic function f(x) that satisfies f(x0) = y0 f&#39;(x0) = y0p f(x1) = y1 f&#39;(x1) = y1p.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MooseUnits pow(const MooseUnits &, int)

◆ saturationHys()

Real PorousFlowVanGenuchten::saturationHys ( Real  pc,
Real  slmin,
Real  sgrdel,
Real  alpha,
Real  n,
const LowCapillaryPressureExtension low_ext = LowCapillaryPressureExtension(),
const HighCapillaryPressureExtension high_ext = HighCapillaryPressureExtension() 
)

Hysteretic saturation function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of Doughty2008), which is the inverse of capillaryPressureHys.

Parameters
pccapillary pressure. 0 <= pc
slminvalue of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
sgrdelvalue of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <= 1
alphavan Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
nvan Genuchten n parameter. n > 1
low_extstrategy and parameters to use for the extension in the small-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)
high_extstrategy and parameters to use for the extension in the high-saturation region (defaults to no extension: this default is not recommended for simulations of real phenomena)

Definition at line 364 of file PorousFlowVanGenuchten.C.

Referenced by PorousFlowHystereticCapillaryPressure::computeTurningPointInfo(), PorousFlowHystereticCapillaryPressure::d2secondOrderDryingSat(), PorousFlowHystereticCapillaryPressure::dsecondOrderDryingSat(), PorousFlowHystereticCapillaryPressure::firstOrderWettingSat(), PorousFlowHystereticCapillaryPressure::liquidSaturationQp(), PorousFlowHystereticCapillaryPressure::secondOrderDryingSat(), and TEST().

371 {
372  if (pc <= 0)
373  return 1.0;
374  Real s = 1.0;
375  if (pc > low_ext.Pc) // important for initialization of the low_ext that this is > and not >=
376  {
377  switch (low_ext.strategy)
378  {
379  case LowCapillaryPressureExtension::QUADRATIC:
380  {
381  const Real s2 = low_ext.S * low_ext.S + 2.0 * (pc - low_ext.Pc) * low_ext.S / low_ext.dPc;
382  if (s2 <= 0.0) // this occurs when we're trying to find a saturation on the wetting curve
383  // defined by sl = sgrDel at pc = Pcd_Del, if this pc is actually impossible
384  // to achieve on this wetting curve
385  s = 0.0;
386  else
387  s = std::sqrt(s2);
388  break;
389  }
390  case LowCapillaryPressureExtension::EXPONENTIAL:
391  {
392  const Real ss = low_ext.S + std::log(pc / low_ext.Pc) * low_ext.Pc / low_ext.dPc;
393  if (ss <= 0.0) // this occurs when we're trying to find a saturation on the
394  // wetting curve defined by sl = sgrDel at pc = Pcd_Del, if this
395  // pc is actually impossible to achieve on this wetting curve
396  s = 0.0;
397  else
398  s = ss;
399  break;
400  }
401  default:
402  s = low_ext.S;
403  }
404  return s;
405  }
406  if (pc < high_ext.Pc) // important for initialization of the high_ext that this is < and not <=
407  {
408  switch (high_ext.strategy)
409  {
410  case HighCapillaryPressureExtension::POWER:
411  {
412  const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
413  s = 1.0 - std::pow(pc / high_ext.Pc, 1.0 / expon) * (1.0 - high_ext.S);
414  break;
415  }
416  default:
417  s = high_ext.S;
418  }
419  return s;
420  }
421  if (pc == std::numeric_limits<Real>::max())
422  s = 0.0;
423  else
424  {
425  const Real seffpow = 1.0 + std::pow(pc * alpha, n);
426  const Real seff = std::pow(seffpow, (1.0 - n) / n);
427  s = (1.0 - sgrdel - slmin) * seff + slmin;
428  }
429  return s;
430 }
if(subdm)
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Definition: NS.h:126
MooseUnits pow(const MooseUnits &, int)