www.mooseframework.org
PorousFlowVanGenuchten.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
9 #include "libmesh/utility.h"
10 
11 namespace PorousFlowVanGenuchten
12 {
13 Real
14 effectiveSaturation(Real p, Real alpha, Real m)
15 {
16  Real n, seff;
17 
18  if (p >= 0.0)
19  return 1.0;
20  else
21  {
22  n = 1.0 / (1.0 - m);
23  seff = 1.0 + std::pow(-alpha * p, n);
24  return std::pow(seff, -m);
25  }
26 }
27 
28 Real
29 dEffectiveSaturation(Real p, Real alpha, Real m)
30 {
31  if (p >= 0.0)
32  return 0.0;
33  else
34  {
35  Real n = 1.0 / (1.0 - m);
36  Real inner = 1.0 + std::pow(-alpha * p, n);
37  Real dinner_dp = -n * alpha * std::pow(-alpha * p, n - 1.0);
38  Real dseff_dp = -m * std::pow(inner, -m - 1) * dinner_dp;
39  return dseff_dp;
40  }
41 }
42 
43 Real
44 d2EffectiveSaturation(Real p, Real alpha, Real m)
45 {
46  if (p >= 0.0)
47  return 0.0;
48  else
49  {
50  Real n = 1.0 / (1.0 - m);
51  Real inner = 1.0 + std::pow(-alpha * p, n);
52  Real dinner_dp = -n * alpha * std::pow(-alpha * p, n - 1.0);
53  Real d2inner_dp2 = n * (n - 1.0) * alpha * alpha * std::pow(-alpha * p, n - 2.0);
54  Real d2seff_dp2 = m * (m + 1.0) * std::pow(inner, -m - 2.0) * std::pow(dinner_dp, 2.0) -
55  m * std::pow(inner, -m - 1.0) * d2inner_dp2;
56  return d2seff_dp2;
57  }
58 }
59 
60 Real
61 capillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
62 {
63  if (seff >= 1.0)
64  return 0.0;
65  else if (seff <= 0.0)
66  return pc_max;
67  else
68  {
69  Real a = std::pow(seff, -1.0 / m) - 1.0;
70  return std::min(std::pow(a, 1.0 - m) / alpha, pc_max);
71  }
72 }
73 
74 Real
75 dCapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
76 {
77  if (seff <= 0.0 || seff >= 1.0)
78  return 0.0;
79  else
80  {
81  Real a = std::pow(seff, -1.0 / m) - 1.0;
82  // Return 0 if pc > pc_max
83  if (std::pow(a, 1.0 - m) / alpha > pc_max)
84  return 0.0;
85  else
86  return (m - 1.0) * std::pow(a, -m) * std::pow(seff, -1.0 - 1.0 / m) / m / alpha;
87  }
88 }
89 
90 Real
91 d2CapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
92 {
93  if (seff <= 0.0 || seff >= 1.0)
94  return 0.0;
95  else
96  {
97  Real a = std::pow(seff, -1.0 / m) - 1.0;
98  // Return 0 if pc > pc_max
99  if (std::pow(a, 1.0 - m) / alpha > pc_max)
100  return 0.0;
101  else
102  {
103  Real d2pc = -std::pow(a, -1.0 - m) * std::pow(seff, -2.0 - 2.0 / m) +
104  ((1.0 + m) / m) * std::pow(a, -m) * std::pow(seff, -1.0 / m - 2.0);
105  d2pc *= (1.0 - m) / m / alpha;
106  return d2pc;
107  }
108  }
109 }
110 
111 Real
112 relativePermeability(Real seff, Real m)
113 {
114  if (seff <= 0.0)
115  return 0.0;
116  else if (seff >= 1.0)
117  return 1.0;
118 
119  const Real a = 1.0 - std::pow(seff, 1.0 / m);
120  const Real b = 1.0 - std::pow(a, m);
121 
122  return std::sqrt(seff) * Utility::pow<2>(b);
123 }
124 
125 Real
126 dRelativePermeability(Real seff, Real m)
127 {
128  // Guard against division by zero
129  if (seff <= 0.0 || seff >= 1.0)
130  return 0.0;
131 
132  const Real a = 1.0 - std::pow(seff, 1.0 / m);
133  const Real da = -1.0 / m * std::pow(seff, 1.0 / m - 1.0);
134  const Real b = 1.0 - std::pow(a, m);
135  const Real db = -m * std::pow(a, m - 1.0) * da;
136 
137  return 0.5 * std::pow(seff, -0.5) * Utility::pow<2>(b) + 2.0 * std::sqrt(seff) * b * db;
138 }
139 
140 Real
141 d2RelativePermeability(Real seff, Real m)
142 {
143  // Guard against division by zero
144  if (seff <= 0.0 || seff >= 1.0)
145  return 0.0;
146 
147  const Real a = 1.0 - std::pow(seff, 1.0 / m);
148  const Real da = -1.0 / m * std::pow(seff, 1.0 / m - 1.0);
149  const Real d2a = -(1.0 / m) * (1.0 / m - 1.0) * std::pow(seff, 1.0 / m - 2.0);
150  const Real b = 1.0 - std::pow(a, m);
151  const Real db = -m * std::pow(a, m - 1.0) * da;
152  const Real d2b = -m * (m - 1.0) * std::pow(a, m - 2.0) * da * da - m * std::pow(a, m - 1.0) * d2a;
153 
154  return -0.25 * std::pow(seff, -1.5) * Utility::pow<2>(b) + 2.0 * std::pow(seff, -0.5) * b * db +
155  2.0 * std::sqrt(seff) * db * db + 2.0 * std::sqrt(seff) * b * d2b;
156 }
157 }
Real relativePermeability(Real seff, Real m)
Relative permeability as a function of effective saturation.
van Genuchten effective saturation, capillary pressure and relative permeability functions.
Real d2CapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
Second derivative of capillary pressure wrt effective saturation.
Real dRelativePermeability(Real seff, Real m)
Derivative of relative permeability with respect to effective saturation.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
Real effectiveSaturation(Real p, Real alpha, Real m)
Effective saturation as a function of porepressure.
Real dCapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
Derivative of capillary pressure wrt effective saturation.
Real d2EffectiveSaturation(Real p, Real alpha, Real m)
Second derivative of effective saturation wrt porepressure.
Real d2RelativePermeability(Real seff, Real m)
Second derivative of relative permeability with respect to effective saturation.
Real dEffectiveSaturation(Real p, Real alpha, Real m)
Derivative of effective saturation wrt porepressure.
Real capillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
Capillary pressure as a function of effective saturation.