www.mooseframework.org
Functions
PorousFlowBroadbridgeWhite Namespace Reference

Broadbridge-White version of relative permeability, and effective saturation as a function of capillary pressure. More...

Functions

Real LambertW (Real z)
 Provides the Lambert W function, which satisfies W(z) * exp(W(z)) = z. More...
 
Real effectiveSaturation (Real pc, Real c, Real sn, Real ss, Real las)
 Effective saturation as a function of capillary pressure If pc>=0 this will yield 1, otherwise it will yield <1. More...
 
Real dEffectiveSaturation (Real pc, Real c, Real sn, Real ss, Real las)
 Derivative of effective saturation wrt capillary pressure. More...
 
Real d2EffectiveSaturation (Real pc, Real c, Real sn, Real ss, Real las)
 Second derivative of effective saturation wrt capillary pressure. More...
 
template<typename T >
relativePermeability (const T &s, Real c, Real sn, Real ss, Real kn, Real ks)
 Relative permeability as a function of saturation. More...
 
Real dRelativePermeability (Real s, Real c, Real sn, Real ss, Real kn, Real ks)
 Derivative of relative permeability with respect to saturation. More...
 
Real d2RelativePermeability (Real s, Real c, Real sn, Real ss, Real kn, Real ks)
 Second derivative of relative permeability with respect to saturation. More...
 

Detailed Description

Broadbridge-White version of relative permeability, and effective saturation as a function of capillary pressure.

Valid for Kn small. P Broadbridge, I White `‘Constant rate rainfall infiltration: A versatile nonlinear model, 1 Analytical solution’'. Water Resources Research 24 (1988) 145–154.

Function Documentation

◆ d2EffectiveSaturation()

Real PorousFlowBroadbridgeWhite::d2EffectiveSaturation ( Real  pc,
Real  c,
Real  sn,
Real  ss,
Real  las 
)

Second derivative of effective saturation wrt capillary pressure.

Parameters
pccapillary pressure
cBW's C parameter
snBW's S_n parameter
ssBW's S_s parameter
lasBW's lambda_s parameter
Returns
second derivative of effective saturation wrt capillary pressure

Definition at line 113 of file PorousFlowBroadbridgeWhite.C.

Referenced by PorousFlowCapillaryPressureBW::d2EffectiveSaturation(), and TEST_F().

114 {
115  if (pressure >= 0.0)
116  return 0.0;
117  const Real x = (c - 1.0) * std::exp(c - 1.0 - c * pressure / las);
118  const Real lamw = LambertW(x);
119  return -(ss - sn) * Utility::pow<3>(c) / Utility::pow<2>(las) * lamw * (1.0 - 2.0 * lamw) /
120  Utility::pow<5>(1.0 + lamw);
121 }
Real LambertW(Real z)
Provides the Lambert W function, which satisfies W(z) * exp(W(z)) = z.
const std::vector< double > x
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string pressure
Definition: NS.h:56

◆ d2RelativePermeability()

Real PorousFlowBroadbridgeWhite::d2RelativePermeability ( Real  s,
Real  c,
Real  sn,
Real  ss,
Real  kn,
Real  ks 
)

Second derivative of relative permeability with respect to saturation.

Parameters
ssaturation
cBW's C parameter
snBW's S_n parameter
ssBW's S_s parameter
knBW's K_n parameter
ksBW's K_s parameter
Returns
second derivative of relative permeability wrt saturation

Definition at line 139 of file PorousFlowBroadbridgeWhite.C.

Referenced by TEST_F().

140 {
141  if (s <= sn)
142  return 0.0;
143 
144  if (s >= ss)
145  return 0.0;
146 
147  const Real coef = (ks - kn) * (c - 1.0);
148  const Real th = (s - sn) / (ss - sn);
149  const Real krelpp = coef * (2.0 / (c - th) + 4.0 * th / Utility::pow<2>(c - th) +
150  2.0 * Utility::pow<2>(th) / Utility::pow<3>(c - th));
151  return krelpp / Utility::pow<2>(ss - sn);
152 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ dEffectiveSaturation()

Real PorousFlowBroadbridgeWhite::dEffectiveSaturation ( Real  pc,
Real  c,
Real  sn,
Real  ss,
Real  las 
)

Derivative of effective saturation wrt capillary pressure.

Parameters
pccapillary pressure
cBW's C parameter
snBW's S_n parameter
ssBW's S_s parameter
lasBW's lambda_s parameter
Returns
derivative of effective saturation wrt capillary pressure

Definition at line 103 of file PorousFlowBroadbridgeWhite.C.

Referenced by PorousFlowCapillaryPressureBW::dEffectiveSaturation(), and TEST_F().

104 {
105  if (pressure >= 0.0)
106  return 0.0;
107  const Real x = (c - 1.0) * std::exp(c - 1.0 - c * pressure / las);
108  const Real lamw = LambertW(x);
109  return (ss - sn) * Utility::pow<2>(c) / las * lamw / Utility::pow<3>(1.0 + lamw);
110 }
Real LambertW(Real z)
Provides the Lambert W function, which satisfies W(z) * exp(W(z)) = z.
const std::vector< double > x
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string pressure
Definition: NS.h:56

◆ dRelativePermeability()

Real PorousFlowBroadbridgeWhite::dRelativePermeability ( Real  s,
Real  c,
Real  sn,
Real  ss,
Real  kn,
Real  ks 
)

Derivative of relative permeability with respect to saturation.

Parameters
ssaturation
cBW's C parameter
snBW's S_n parameter
ssBW's S_s parameter
knBW's K_n parameter
ksBW's K_s parameter
Returns
derivative of relative permeability wrt saturation

Definition at line 124 of file PorousFlowBroadbridgeWhite.C.

Referenced by PorousFlowRelativePermeabilityBaseTempl< is_ad >::computeQpProperties(), PorousFlowRelativePermeabilityBWTempl< is_ad >::dRelativePermeability(), and TEST_F().

125 {
126  if (s <= sn)
127  return 0.0;
128 
129  if (s >= ss)
130  return 0.0;
131 
132  const Real coef = (ks - kn) * (c - 1.0);
133  const Real th = (s - sn) / (ss - sn);
134  const Real krelp = coef * (2.0 * th / (c - th) + Utility::pow<2>(th) / Utility::pow<2>(c - th));
135  return krelp / (ss - sn);
136 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ effectiveSaturation()

Real PorousFlowBroadbridgeWhite::effectiveSaturation ( Real  pc,
Real  c,
Real  sn,
Real  ss,
Real  las 
)

Effective saturation as a function of capillary pressure If pc>=0 this will yield 1, otherwise it will yield <1.

Parameters
pccapillary pressure
cBW's C parameter
snBW's S_n parameter
ssBW's S_s parameter
lasBW's lambda_s parameter
Returns
effective saturation

Definition at line 93 of file PorousFlowBroadbridgeWhite.C.

Referenced by PorousFlowRelativePermeabilityBaseTempl< is_ad >::computeQpProperties(), PorousFlowCapillaryPressureBW::effectiveSaturation(), and TEST_F().

94 {
95  if (pressure >= 0.0)
96  return 1.0;
97  const Real x = (c - 1.0) * std::exp(c - 1.0 - c * pressure / las);
98  const Real th = c / (1.0 + LambertW(x)); // use branch 0 for positive x
99  return sn + (ss - sn) * th;
100 }
Real LambertW(Real z)
Provides the Lambert W function, which satisfies W(z) * exp(W(z)) = z.
const std::vector< double > x
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string pressure
Definition: NS.h:56

◆ LambertW()

Real PorousFlowBroadbridgeWhite::LambertW ( Real  z)

Provides the Lambert W function, which satisfies W(z) * exp(W(z)) = z.

Parameters
zz value, which must be positive
Returns
W

Definition at line 15 of file PorousFlowBroadbridgeWhite.C.

Referenced by d2EffectiveSaturation(), dEffectiveSaturation(), and effectiveSaturation().

16 {
17  /* Lambert W function.
18  Was ~/C/LambertW.c written K M Briggs Keith dot Briggs at bt dot com 97 May 21.
19  Revised KMB 97 Nov 20; 98 Feb 11, Nov 24, Dec 28; 99 Jan 13; 00 Feb 23; 01 Apr 09
20 
21  Computes Lambert W function, principal branch.
22  See LambertW1.c for -1 branch.
23 
24  Returned value W(z) satisfies W(z)*exp(W(z))=z
25  test data...
26  W(1)= 0.5671432904097838730
27  W(2)= 0.8526055020137254914
28  W(20)=2.2050032780240599705
29  To solve (a+b*R)*exp(-c*R)-d=0 for R, use
30  R=-(b*W(-exp(-a*c/b)/b*d*c)+a*c)/b/c
31 
32  Test:
33  gcc -DTESTW LambertW.c -o LambertW -lm && LambertW
34  Library:
35  gcc -O3 -c LambertW.c
36 
37  Modified trivially by Andy to use MOOSE things
38  */
39  mooseAssert(z > 0, "LambertW function in RichardsSeff1BWsmall called with negative argument");
40 
41  const Real eps = 4.0e-16; //, em1=0.3678794411714423215955237701614608;
42  Real p, e, t, w;
43 
44  /* Uncomment this stuff is you ever need to call with a negative argument
45  if (z < -em1)
46  mooseError("LambertW: bad argument ", z, "\n");
47 
48  if (0.0 == z)
49  return 0.0;
50  if (z < -em1+1e-4)
51  {
52  // series near -em1 in sqrt(q)
53  Real q=z+em1,r=std::sqrt(q),q2=q*q,q3=q2*q;
54  return
55  -1.0
56  +2.331643981597124203363536062168*r
57  -1.812187885639363490240191647568*q
58  +1.936631114492359755363277457668*r*q
59  -2.353551201881614516821543561516*q2
60  +3.066858901050631912893148922704*r*q2
61  -4.175335600258177138854984177460*q3
62  +5.858023729874774148815053846119*r*q3
63  -8.401032217523977370984161688514*q3*q; // error approx 1e-16
64  }
65  */
66  /* initial approx for iteration... */
67  if (z < 1.0)
68  {
69  /* series near 0 */
70  p = std::sqrt(2.0 * (2.7182818284590452353602874713526625 * z + 1.0));
71  w = -1.0 + p * (1.0 + p * (-0.333333333333333333333 + p * 0.152777777777777777777777));
72  }
73  else
74  w = std::log(z); /* asymptotic */
75  if (z > 3.0)
76  w -= std::log(w); /* useful? */
77  for (unsigned int i = 0; i < 10; ++i)
78  {
79  /* Halley iteration */
80  e = std::exp(w);
81  t = w * e - z;
82  p = w + 1.0;
83  t /= e * p - 0.5 * (p + 1.0) * t / p;
84  w -= t;
85  if (std::abs(t) < eps * (1.0 + std::abs(w)))
86  return w; /* rel-abs error */
87  }
88  /* should never get here */
89  mooseError("LambertW: No convergence at z= ", z, "\n");
90 }
const Real eps
void mooseError(Args &&... args)
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ relativePermeability()

template<typename T >
T PorousFlowBroadbridgeWhite::relativePermeability ( const T &  s,
Real  c,
Real  sn,
Real  ss,
Real  kn,
Real  ks 
)

Relative permeability as a function of saturation.

Parameters
ssaturation
cBW's C parameter
snBW's S_n parameter
ssBW's S_s parameter
knBW's K_n parameter
ksBW's K_s parameter
Returns
relative permeability

Definition at line 80 of file PorousFlowBroadbridgeWhite.h.

Referenced by PorousFlowRelativePermeabilityBaseTempl< is_ad >::computeQpProperties(), PorousFlowRelativePermeabilityBWTempl< is_ad >::relativePermeability(), and TEST_F().

81 {
82  if (MetaPhysicL::raw_value(s) <= sn)
83  return kn;
84 
85  if (MetaPhysicL::raw_value(s) >= ss)
86  return ks;
87 
88  const T coef = (ks - kn) * (c - 1.0);
89  const T th = (s - sn) / (ss - sn);
90  const T krel = kn + coef * Utility::pow<2>(th) / (c - th);
91  return krel;
92 }
auto raw_value(const Eigen::Map< T > &in)