Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "MooseTypes.h"
13 : #include "libmesh/utility.h"
14 :
15 : /**
16 : * van Genuchten effective saturation, capillary pressure and relative
17 : * permeability functions.
18 : * Note: effective saturation is provided as a function of porepressure, not
19 : * capillary pressure.
20 : * Note: capillary pressure and relative permeability are functions of effective
21 : * saturation. The derivatives are therefore given wrt effective saturation. These
22 : * derivatives must be multiplied by the derivative of effective saturation wrt
23 : * the true saturation in objects using these relations.
24 : *
25 : * Based on van Genuchten, M. Th., A closed for equation for
26 : * predicting the hydraulic conductivity of unsaturated soils,
27 : * Soil Sci. Soc., 44, 892-898 (1980).
28 : */
29 :
30 : namespace PorousFlowVanGenuchten
31 : {
32 : /**
33 : * Effective saturation as a function of porepressure.
34 : * Note: seff = 1 for p >= 0
35 : * @param p porepressure
36 : * @param alpha van Genuchten parameter
37 : * @param m van Genuchten exponent
38 : * @return effective saturation
39 : */
40 : Real effectiveSaturation(Real p, Real alpha, Real m);
41 :
42 : /**
43 : * Derivative of effective saturation wrt porepressure
44 : * @param p porepressure
45 : * @param alpha van Genuchten parameter
46 : * @param m van Genuchten exponent
47 : * @return derivative of effective saturation wrt porepressure
48 : */
49 : Real dEffectiveSaturation(Real p, Real alpha, Real m);
50 :
51 : /**
52 : * Second derivative of effective saturation wrt porepressure
53 : * @param p porepressure
54 : * @param alpha van Genuchten parameter
55 : * @param m van Genuchten exponent
56 : * @return second derivative of effective saturation wrt porepressure
57 : */
58 : Real d2EffectiveSaturation(Real p, Real alpha, Real m);
59 :
60 : /**
61 : * Capillary pressure as a function of effective saturation
62 : *
63 : * @param seff effective saturation
64 : * @param alpha van Genuchten parameter
65 : * @param m van Genuchten exponent
66 : * @param pc_max maximum capillary pressure (Pa)
67 : * @return capillary pressure (Pa)
68 : */
69 : Real capillaryPressure(Real seff, Real alpha, Real m, Real pc_max);
70 :
71 : /**
72 : * Derivative of capillary pressure wrt effective saturation
73 : *
74 : * @param seff effective saturation
75 : * @param alpha van Genuchten parameter
76 : * @param m van Genuchten exponent
77 : * @param pc_max maximum capillary pressure (Pa)
78 : * @return derivative of capillary pressure wrt effective saturation
79 : */
80 : Real dCapillaryPressure(Real seff, Real alpha, Real m, Real pc_max);
81 :
82 : /**
83 : * Second derivative of capillary pressure wrt effective saturation
84 : *
85 : * @param seff effective saturation
86 : * @param alpha van Genuchten parameter
87 : * @param m van Genuchten exponent
88 : * @param pc_max maximum capillary pressure (Pa)
89 : * @return second derivative of capillary pressure wrt effective saturation
90 : */
91 : Real d2CapillaryPressure(Real seff, Real alpha, Real m, Real pc_max);
92 :
93 : /**
94 : * Relative permeability as a function of effective saturation
95 : * @param seff effective saturation
96 : * @param m van Genuchten exponent
97 : * @return relative permeability
98 : */
99 : template <typename T>
100 : T
101 751015 : relativePermeability(const T & seff, Real m)
102 : {
103 751015 : if (MetaPhysicL::raw_value(seff) <= 0.0)
104 0 : return 0.0;
105 751014 : else if (MetaPhysicL::raw_value(seff) >= 1.0)
106 0 : return 1.0;
107 :
108 : using std::pow, std::sqrt;
109 :
110 750450 : const T a = 1.0 - pow(seff, 1.0 / m);
111 750450 : const T b = 1.0 - pow(a, m);
112 :
113 750450 : return sqrt(seff) * Utility::pow<2>(b);
114 : }
115 :
116 : /**
117 : * Derivative of relative permeability with respect to effective saturation
118 : * @param seff effective saturation
119 : * @param m van Genuchten exponent
120 : * @return derivative of relative permeability wrt effective saturation
121 : */
122 : Real dRelativePermeability(Real seff, Real m);
123 :
124 : /**
125 : * Second derivative of relative permeability with respect to effective saturation
126 : * @param seff effective saturation
127 : * @param m van Genuchten exponent
128 : * @return second derivative of relative permeability wrt effective saturation
129 : */
130 : Real d2RelativePermeability(Real seff, Real m);
131 :
132 : /**
133 : * Relative permeability for a non-wetting phase as a function of effective saturation
134 : * @param seff effective saturation
135 : * @param m van Genuchten exponent
136 : * @return relative permeability
137 : */
138 : template <typename T>
139 : T
140 9322 : relativePermeabilityNW(const T & seff, Real m)
141 : {
142 9322 : if (MetaPhysicL::raw_value(seff) <= 0.0)
143 0 : return 0.0;
144 9321 : else if (MetaPhysicL::raw_value(seff) >= 1.0)
145 0 : return 1.0;
146 :
147 : using std::pow, std::sqrt;
148 :
149 9189 : const T a = pow(1.0 - seff, 1.0 / m);
150 9189 : const T b = pow(1.0 - a, 2.0 * m);
151 :
152 9189 : return sqrt(seff) * b;
153 : }
154 :
155 : /**
156 : * Derivative of relative permeability for a non-wetting phase with respect to effective saturation
157 : * @param seff effective saturation
158 : * @param m van Genuchten exponent
159 : * @return derivative of relative permeability wrt effective saturation
160 : */
161 : Real dRelativePermeabilityNW(Real seff, Real m);
162 :
163 : /**
164 : * Second derivative of relative permeability for a non-wetting phase with respect to effective
165 : * saturation
166 : * @param seff effective saturation
167 : * @param m van Genuchten exponent
168 : * @return second derivative of relative permeability wrt effective saturation
169 : */
170 : Real d2RelativePermeabilityNW(Real seff, Real m);
171 :
172 : /**
173 : * Parameters associated with the extension of the hysteretic capillary pressure function to low
174 : * saturation values
175 : * @ ExtensionStrategy the type of extension used
176 : * @ S liquid saturation at the point of extension
177 : * @ Pc capillary pressure at the point of extension
178 : * @ dPc d(Pc)/dS at the point of extension
179 : */
180 : struct LowCapillaryPressureExtension
181 : {
182 : enum ExtensionStrategy
183 : {
184 : NONE,
185 : QUADRATIC,
186 : EXPONENTIAL
187 : };
188 : ExtensionStrategy strategy;
189 : Real S;
190 : Real Pc;
191 : Real dPc;
192 :
193 : LowCapillaryPressureExtension()
194 954360 : : strategy(LowCapillaryPressureExtension::NONE),
195 954360 : S(0.0),
196 954360 : Pc(std::numeric_limits<Real>::max()),
197 14544 : dPc(std::numeric_limits<Real>::lowest()){};
198 :
199 : LowCapillaryPressureExtension(const ExtensionStrategy & strategy, Real S, Real Pc, Real dPc)
200 454368 : : strategy(strategy), S(S), Pc(Pc), dPc(dPc){};
201 : };
202 :
203 : /**
204 : * Parameters associated with the extension of the hysteretic wetting capillary pressure function to
205 : * high saturation values
206 : * @ ExtensionStrategy the type of extension used
207 : * @ S liquid saturation at the point of extension
208 : * @ Pc capillary pressure at the point of extension
209 : * @ dPc d(Pc)/dS at the point of extension
210 : */
211 : struct HighCapillaryPressureExtension
212 : {
213 : enum ExtensionStrategy
214 : {
215 : NONE,
216 : POWER
217 : };
218 : ExtensionStrategy strategy;
219 : Real S;
220 : Real Pc;
221 : Real dPc;
222 :
223 : HighCapillaryPressureExtension()
224 5962965 : : strategy(HighCapillaryPressureExtension::NONE),
225 5962965 : S(1.0),
226 5962965 : Pc(0.0),
227 4353711 : dPc(std::numeric_limits<Real>::lowest()){};
228 :
229 : HighCapillaryPressureExtension(const ExtensionStrategy & strategy, Real S, Real Pc, Real dPc)
230 454368 : : strategy(strategy), S(S), Pc(Pc), dPc(dPc){};
231 : };
232 :
233 : /**
234 : * Hysteretic capillary pressure function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of
235 : * Doughty2008).
236 : * NOTE: this function is undefined for sl < 0 and sl > 1, so you MUST ensure 0 <= sl <= 1
237 : * NOTE: this returns a non-negative quantity.
238 : * @param sl liquid saturation. 0 <= sl <= 1
239 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
240 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
241 : * 1
242 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
243 : * @param n van Genuchten n parameter. n > 1
244 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
245 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
246 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
247 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
248 : */
249 : Real capillaryPressureHys(
250 : Real sl,
251 : Real slmin,
252 : Real sgrdel,
253 : Real alpha,
254 : Real n,
255 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
256 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
257 :
258 : /**
259 : * Derivative of capillaryPressureHys with respect to sl.
260 : * NOTE: this function is undefined for sl < 0 and sl > 1, so you MUST ensure 0 <= sl <= 1
261 : * NOTE: this returns a negative quantity.
262 : * @param sl liquid saturation. 0 <= sl <= 1
263 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
264 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
265 : * 1
266 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
267 : * @param n van Genuchten n parameter. n > 1
268 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
269 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
270 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
271 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
272 : */
273 : Real dcapillaryPressureHys(
274 : Real sl,
275 : Real slmin,
276 : Real sgrdel,
277 : Real alpha,
278 : Real n,
279 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
280 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
281 :
282 : /**
283 : * Second derivative of capillaryPressureHys with respect to sl.
284 : * NOTE: this function is undefined for sl < 0 and sl > 1, so you MUST ensure 0 <= sl <= 1
285 : * @param sl liquid saturation. 0 <= sl <= 1
286 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
287 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
288 : * 1
289 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
290 : * @param n van Genuchten n parameter. n > 1
291 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
292 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
293 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
294 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
295 : */
296 : Real d2capillaryPressureHys(
297 : Real sl,
298 : Real slmin,
299 : Real sgrdel,
300 : Real alpha,
301 : Real n,
302 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
303 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
304 :
305 : /**
306 : * Hysteretic saturation function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of
307 : * Doughty2008), which is the inverse of capillaryPressureHys
308 : * @param pc capillary pressure. 0 <= pc
309 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
310 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
311 : * 1
312 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
313 : * @param n van Genuchten n parameter. n > 1
314 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
315 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
316 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
317 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
318 : */
319 : Real
320 : saturationHys(Real pc,
321 : Real slmin,
322 : Real sgrdel,
323 : Real alpha,
324 : Real n,
325 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
326 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
327 :
328 : /**
329 : * Derivative of Hysteretic saturation function with respect to pc
330 : * @param pc capillary pressure. 0 <= pc
331 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
332 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
333 : * 1
334 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
335 : * @param n van Genuchten n parameter. n > 1
336 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
337 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
338 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
339 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
340 : */
341 : Real
342 : dsaturationHys(Real pc,
343 : Real slmin,
344 : Real sgrdel,
345 : Real alpha,
346 : Real n,
347 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
348 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
349 :
350 : /**
351 : * Second derivative of Hysteretic saturation function with respect to pc
352 : * @param pc capillary pressure. 0 <= pc
353 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
354 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
355 : * 1
356 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
357 : * @param n van Genuchten n parameter. n > 1
358 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
359 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
360 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
361 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
362 : */
363 : Real
364 : d2saturationHys(Real pc,
365 : Real slmin,
366 : Real sgrdel,
367 : Real alpha,
368 : Real n,
369 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
370 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
371 :
372 : /**
373 : * Hysteretic relative permeability for liquid
374 : * @param sl liquid saturation
375 : * @param slr residual liquid saturation. For sl < slr, this function will always return 0
376 : * @param sgrdel value of gas saturation where van Genuchten wetting capillary-pressure expression
377 : * -> 0. This depends on the turning-point saturation when drying became wetting, using the Land
378 : * equation
379 : * @param sgrmax maximum value possible for sgrdel. This will be equal to sgrdel if the
380 : * turning-point saturation is small
381 : * @param sldel value of the turning-point saturation when drying became wetting
382 : * @param m van-Genuchten m parameter
383 : * @param upper_liquid_param cubic modification of the wetting relative permeability will occur
384 : * between upper_liquid_param * (1 - sgrdel) and 1 - 0.5 * sgrdel. 0 < upper_liquid_param <= 1.
385 : * Usually upper_liquid_param is close to 1 (eg 0.9)
386 : * @param y0 value of the unmodified wetting relative permeability at sl = upper_liquid_param * (1 -
387 : * sgrdel)
388 : * @param y0p value of the derivtaive of the unmodified wetting relative permeability at sl =
389 : * upper_liquid_param * (1 - sgrdel)
390 : * @param y1 value of the unmodified wetting relative permeability at sl = 1 - 0.5 * sgrdel
391 : * @param y1p value of the derivtaive of the unmodified wetting relative permeability at sl = 1 -
392 : * 0.5 * sgrdel
393 : */
394 : Real relativePermeabilityHys(Real sl,
395 : Real slr,
396 : Real sgrdel,
397 : Real sgrmax,
398 : Real sldel,
399 : Real m,
400 : Real upper_liquid_param,
401 : Real y0,
402 : Real y0p,
403 : Real y1,
404 : Real y1p);
405 :
406 : /**
407 : * Derivative of Hysteretic relative permeability for liquid, with respect to liquid saturation
408 : * @param sl liquid saturation
409 : * @param slr residual liquid saturation. For sl < slr, this function will always return 0
410 : * @param sgrdel value of gas saturation where van Genuchten wetting capillary-pressure expression
411 : * -> 0. This depends on the turning-point saturation when drying became wetting, using the Land
412 : * equation
413 : * @param sgrmax maximum value possible for sgrdel. This will be equal to sgrdel if the
414 : * turning-point saturation is small
415 : * @param sldel value of the turning-point saturation when drying became wetting
416 : * @param m van-Genuchten m parameter
417 : * @param upper_liquid_param cubic modification of the wetting relative permeability will occur
418 : * between upper_liquid_param * (1 - sgrdel) and 1 - 0.5 * sgrdel. 0 < upper_liquid_param <= 1.
419 : * Usually upper_liquid_param is close to 1 (eg 0.9)
420 : * @param y0 value of the unmodified wetting relative permeability at sl = upper_liquid_param * (1 -
421 : * sgrdel)
422 : * @param y0p value of the derivtaive of the unmodified wetting relative permeability at sl =
423 : * upper_liquid_param * (1 - sgrdel)
424 : * @param y1 value of the unmodified wetting relative permeability at sl = 1 - 0.5 * sgrdel
425 : * @param y1p value of the derivtaive of the unmodified wetting relative permeability at sl = 1 -
426 : * 0.5 * sgrdel
427 : */
428 : Real drelativePermeabilityHys(Real sl,
429 : Real slr,
430 : Real sgrdel,
431 : Real sgrmax,
432 : Real sldel,
433 : Real m,
434 : Real upper_liquid_param,
435 : Real y0,
436 : Real y0p,
437 : Real y1,
438 : Real y1p);
439 :
440 : /**
441 : * Hysteretic relative permeability for gas
442 : * @param sl liquid saturation
443 : * @param slr residual liquid saturation. For sl < slr, this function will always return 0
444 : * @param sgrdel value of gas saturation where van Genuchten wetting capillary-pressure expression
445 : * -> 0. This depends on the turning-point saturation when drying became wetting, using the Land
446 : * equation
447 : * @param sgrmax maximum value possible for sgrdel. This will be equal to sgrdel if the
448 : * turning-point saturation is small
449 : * @param sldel value of the turning-point saturation when drying became wetting
450 : * @param m van-Genuchten m parameter
451 : * @param gamma index satisfying gamma > 0. Usually gamma = 1/3.
452 : * @param k_rg_max Maximum value of unextended gas relative permeability. If no low-saturation
453 : * extension is used then gas relative permeability = k_rg_max for sl <= slr. Satisfies 0 <
454 : * k_rg_max <= 1. Frequently k_rg_max = 1 is used
455 : * @param y0p Value of the derivative of the low-saturation extension at sl = slr. If an extension
456 : * is used then the gas relative permeability in the region sl <= slr is a cubic whose value is
457 : * unity at sl = 0, and derivative is zero at sl = 0
458 : */
459 : Real relativePermeabilityNWHys(Real sl,
460 : Real slr,
461 : Real sgrdel,
462 : Real sgrmax,
463 : Real sldel,
464 : Real m,
465 : Real gamma,
466 : Real k_rg_max,
467 : Real y0p);
468 :
469 : /**
470 : * Derivative of hysteretic relative permeability for gas with respect to the liquid saturation
471 : * @param sl liquid saturation
472 : * @param slr residual liquid saturation. For sl < slr, this function will always return 0
473 : * @param sgrdel value of gas saturation where van Genuchten wetting capillary-pressure expression
474 : * -> 0. This depends on the turning-point saturation when drying became wetting, using the Land
475 : * equation
476 : * @param sgrmax maximum value possible for sgrdel. This will be equal to sgrdel if the
477 : * turning-point saturation is small
478 : * @param sldel value of the turning-point saturation when drying became wetting
479 : * @param m van-Genuchten m parameter
480 : * @param gamma index satisfying gamma > 0. Usually gamma = 1/3.
481 : * @param k_rg_max Maximum value of unextended gas relative permeability. If no low-saturation
482 : * extension is used then gas relative permeability = k_rg_max for sl <= slr. Satisfies 0 <
483 : * k_rg_max <= 1. Frequently k_rg_max = 1 is used
484 : * @param y0p Value of the derivative of the low-saturation extension at sl = slr. If an extension
485 : * is used then the gas relative permeability in the region sl <= slr is a cubic whose value is
486 : * unity at sl = 0, and derivative is zero at sl = 0
487 : */
488 : Real drelativePermeabilityNWHys(Real sl,
489 : Real slr,
490 : Real sgrdel,
491 : Real sgrmax,
492 : Real sldel,
493 : Real m,
494 : Real gamma,
495 : Real k_rg_max,
496 : Real y0p);
497 : }
|
