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 750450 : const T a = 1.0 - std::pow(seff, 1.0 / m);
109 750450 : const T b = 1.0 - std::pow(a, m);
110 :
111 750450 : return std::sqrt(seff) * Utility::pow<2>(b);
112 : }
113 :
114 : /**
115 : * Derivative of relative permeability with respect to effective saturation
116 : * @param seff effective saturation
117 : * @param m van Genuchten exponent
118 : * @return derivative of relative permeability wrt effective saturation
119 : */
120 : Real dRelativePermeability(Real seff, Real m);
121 :
122 : /**
123 : * Second derivative of relative permeability with respect to effective saturation
124 : * @param seff effective saturation
125 : * @param m van Genuchten exponent
126 : * @return second derivative of relative permeability wrt effective saturation
127 : */
128 : Real d2RelativePermeability(Real seff, Real m);
129 :
130 : /**
131 : * Relative permeability for a non-wetting phase as a function of effective saturation
132 : * @param seff effective saturation
133 : * @param m van Genuchten exponent
134 : * @return relative permeability
135 : */
136 : template <typename T>
137 : T
138 9322 : relativePermeabilityNW(const T & seff, Real m)
139 : {
140 9322 : if (MetaPhysicL::raw_value(seff) <= 0.0)
141 0 : return 0.0;
142 9321 : else if (MetaPhysicL::raw_value(seff) >= 1.0)
143 0 : return 1.0;
144 :
145 9189 : const T a = std::pow(1.0 - seff, 1.0 / m);
146 9189 : const T b = std::pow(1.0 - a, 2.0 * m);
147 :
148 9189 : return std::sqrt(seff) * b;
149 : }
150 :
151 : /**
152 : * Derivative of relative permeability for a non-wetting phase with respect to effective saturation
153 : * @param seff effective saturation
154 : * @param m van Genuchten exponent
155 : * @return derivative of relative permeability wrt effective saturation
156 : */
157 : Real dRelativePermeabilityNW(Real seff, Real m);
158 :
159 : /**
160 : * Second derivative of relative permeability for a non-wetting phase with respect to effective
161 : * saturation
162 : * @param seff effective saturation
163 : * @param m van Genuchten exponent
164 : * @return second derivative of relative permeability wrt effective saturation
165 : */
166 : Real d2RelativePermeabilityNW(Real seff, Real m);
167 :
168 : /**
169 : * Parameters associated with the extension of the hysteretic capillary pressure function to low
170 : * saturation values
171 : * @ ExtensionStrategy the type of extension used
172 : * @ S liquid saturation at the point of extension
173 : * @ Pc capillary pressure at the point of extension
174 : * @ dPc d(Pc)/dS at the point of extension
175 : */
176 : struct LowCapillaryPressureExtension
177 : {
178 : enum ExtensionStrategy
179 : {
180 : NONE,
181 : QUADRATIC,
182 : EXPONENTIAL
183 : };
184 : ExtensionStrategy strategy;
185 : Real S;
186 : Real Pc;
187 : Real dPc;
188 :
189 : LowCapillaryPressureExtension()
190 954576 : : strategy(LowCapillaryPressureExtension::NONE),
191 954576 : S(0.0),
192 954576 : Pc(std::numeric_limits<Real>::max()),
193 14544 : dPc(std::numeric_limits<Real>::lowest()){};
194 :
195 : LowCapillaryPressureExtension(const ExtensionStrategy & strategy, Real S, Real Pc, Real dPc)
196 454476 : : strategy(strategy), S(S), Pc(Pc), dPc(dPc){};
197 : };
198 :
199 : /**
200 : * Parameters associated with the extension of the hysteretic wetting capillary pressure function to
201 : * high saturation values
202 : * @ ExtensionStrategy the type of extension used
203 : * @ S liquid saturation at the point of extension
204 : * @ Pc capillary pressure at the point of extension
205 : * @ dPc d(Pc)/dS at the point of extension
206 : */
207 : struct HighCapillaryPressureExtension
208 : {
209 : enum ExtensionStrategy
210 : {
211 : NONE,
212 : POWER
213 : };
214 : ExtensionStrategy strategy;
215 : Real S;
216 : Real Pc;
217 : Real dPc;
218 :
219 : HighCapillaryPressureExtension()
220 5964711 : : strategy(HighCapillaryPressureExtension::NONE),
221 5964711 : S(1.0),
222 5964711 : Pc(0.0),
223 4355133 : dPc(std::numeric_limits<Real>::lowest()){};
224 :
225 : HighCapillaryPressureExtension(const ExtensionStrategy & strategy, Real S, Real Pc, Real dPc)
226 454476 : : strategy(strategy), S(S), Pc(Pc), dPc(dPc){};
227 : };
228 :
229 : /**
230 : * Hysteretic capillary pressure function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of
231 : * Doughty2008).
232 : * NOTE: this function is undefined for sl < 0 and sl > 1, so you MUST ensure 0 <= sl <= 1
233 : * NOTE: this returns a non-negative quantity.
234 : * @param sl liquid saturation. 0 <= sl <= 1
235 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
236 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
237 : * 1
238 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
239 : * @param n van Genuchten n parameter. n > 1
240 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
241 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
242 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
243 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
244 : */
245 : Real capillaryPressureHys(
246 : Real sl,
247 : Real slmin,
248 : Real sgrdel,
249 : Real alpha,
250 : Real n,
251 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
252 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
253 :
254 : /**
255 : * Derivative of capillaryPressureHys with respect to sl.
256 : * NOTE: this function is undefined for sl < 0 and sl > 1, so you MUST ensure 0 <= sl <= 1
257 : * NOTE: this returns a negative quantity.
258 : * @param sl liquid saturation. 0 <= sl <= 1
259 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
260 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
261 : * 1
262 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
263 : * @param n van Genuchten n parameter. n > 1
264 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
265 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
266 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
267 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
268 : */
269 : Real dcapillaryPressureHys(
270 : Real sl,
271 : Real slmin,
272 : Real sgrdel,
273 : Real alpha,
274 : Real n,
275 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
276 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
277 :
278 : /**
279 : * Second derivative of capillaryPressureHys with respect to sl.
280 : * NOTE: this function is undefined for sl < 0 and sl > 1, so you MUST ensure 0 <= sl <= 1
281 : * @param sl liquid saturation. 0 <= sl <= 1
282 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
283 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
284 : * 1
285 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
286 : * @param n van Genuchten n parameter. n > 1
287 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
288 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
289 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
290 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
291 : */
292 : Real d2capillaryPressureHys(
293 : Real sl,
294 : Real slmin,
295 : Real sgrdel,
296 : Real alpha,
297 : Real n,
298 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
299 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
300 :
301 : /**
302 : * Hysteretic saturation function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of
303 : * Doughty2008), which is the inverse of capillaryPressureHys
304 : * @param pc capillary pressure. 0 <= pc
305 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
306 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
307 : * 1
308 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
309 : * @param n van Genuchten n parameter. n > 1
310 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
311 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
312 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
313 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
314 : */
315 : Real
316 : saturationHys(Real pc,
317 : Real slmin,
318 : Real sgrdel,
319 : Real alpha,
320 : Real n,
321 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
322 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
323 :
324 : /**
325 : * Derivative of Hysteretic saturation function with respect to pc
326 : * @param pc capillary pressure. 0 <= pc
327 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
328 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
329 : * 1
330 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
331 : * @param n van Genuchten n parameter. n > 1
332 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
333 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
334 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
335 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
336 : */
337 : Real
338 : dsaturationHys(Real pc,
339 : Real slmin,
340 : Real sgrdel,
341 : Real alpha,
342 : Real n,
343 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
344 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
345 :
346 : /**
347 : * Second derivative of Hysteretic saturation function with respect to pc
348 : * @param pc capillary pressure. 0 <= pc
349 : * @param slmin value of liquid sat where the van Genuchten expression -> infinity. 0 <= slmin < 1
350 : * @param sgrdel value of gas saturation where van Genuchten expression -> 0. slmin < 1 - Sgrdel <=
351 : * 1
352 : * @param alpha van Genuchten alpha parameter, with dimensions 1/Pa. alpha > 0
353 : * @param n van Genuchten n parameter. n > 1
354 : * @param low_ext strategy and parameters to use for the extension in the small-saturation region
355 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
356 : * @param high_ext strategy and parameters to use for the extension in the high-saturation region
357 : * (defaults to no extension: this default is not recommended for simulations of real phenomena)
358 : */
359 : Real
360 : d2saturationHys(Real pc,
361 : Real slmin,
362 : Real sgrdel,
363 : Real alpha,
364 : Real n,
365 : const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
366 : const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());
367 :
368 : /**
369 : * Hysteretic relative permeability for liquid
370 : * @param sl liquid saturation
371 : * @param slr residual liquid saturation. For sl < slr, this function will always return 0
372 : * @param sgrdel value of gas saturation where van Genuchten wetting capillary-pressure expression
373 : * -> 0. This depends on the turning-point saturation when drying became wetting, using the Land
374 : * equation
375 : * @param sgrmax maximum value possible for sgrdel. This will be equal to sgrdel if the
376 : * turning-point saturation is small
377 : * @param sldel value of the turning-point saturation when drying became wetting
378 : * @param m van-Genuchten m parameter
379 : * @param upper_liquid_param cubic modification of the wetting relative permeability will occur
380 : * between upper_liquid_param * (1 - sgrdel) and 1 - 0.5 * sgrdel. 0 < upper_liquid_param <= 1.
381 : * Usually upper_liquid_param is close to 1 (eg 0.9)
382 : * @param y0 value of the unmodified wetting relative permeability at sl = upper_liquid_param * (1 -
383 : * sgrdel)
384 : * @param y0p value of the derivtaive of the unmodified wetting relative permeability at sl =
385 : * upper_liquid_param * (1 - sgrdel)
386 : * @param y1 value of the unmodified wetting relative permeability at sl = 1 - 0.5 * sgrdel
387 : * @param y1p value of the derivtaive of the unmodified wetting relative permeability at sl = 1 -
388 : * 0.5 * sgrdel
389 : */
390 : Real relativePermeabilityHys(Real sl,
391 : Real slr,
392 : Real sgrdel,
393 : Real sgrmax,
394 : Real sldel,
395 : Real m,
396 : Real upper_liquid_param,
397 : Real y0,
398 : Real y0p,
399 : Real y1,
400 : Real y1p);
401 :
402 : /**
403 : * Derivative of Hysteretic relative permeability for liquid, with respect to liquid saturation
404 : * @param sl liquid saturation
405 : * @param slr residual liquid saturation. For sl < slr, this function will always return 0
406 : * @param sgrdel value of gas saturation where van Genuchten wetting capillary-pressure expression
407 : * -> 0. This depends on the turning-point saturation when drying became wetting, using the Land
408 : * equation
409 : * @param sgrmax maximum value possible for sgrdel. This will be equal to sgrdel if the
410 : * turning-point saturation is small
411 : * @param sldel value of the turning-point saturation when drying became wetting
412 : * @param m van-Genuchten m parameter
413 : * @param upper_liquid_param cubic modification of the wetting relative permeability will occur
414 : * between upper_liquid_param * (1 - sgrdel) and 1 - 0.5 * sgrdel. 0 < upper_liquid_param <= 1.
415 : * Usually upper_liquid_param is close to 1 (eg 0.9)
416 : * @param y0 value of the unmodified wetting relative permeability at sl = upper_liquid_param * (1 -
417 : * sgrdel)
418 : * @param y0p value of the derivtaive of the unmodified wetting relative permeability at sl =
419 : * upper_liquid_param * (1 - sgrdel)
420 : * @param y1 value of the unmodified wetting relative permeability at sl = 1 - 0.5 * sgrdel
421 : * @param y1p value of the derivtaive of the unmodified wetting relative permeability at sl = 1 -
422 : * 0.5 * sgrdel
423 : */
424 : Real drelativePermeabilityHys(Real sl,
425 : Real slr,
426 : Real sgrdel,
427 : Real sgrmax,
428 : Real sldel,
429 : Real m,
430 : Real upper_liquid_param,
431 : Real y0,
432 : Real y0p,
433 : Real y1,
434 : Real y1p);
435 :
436 : /**
437 : * Hysteretic relative permeability for gas
438 : * @param sl liquid saturation
439 : * @param slr residual liquid saturation. For sl < slr, this function will always return 0
440 : * @param sgrdel value of gas saturation where van Genuchten wetting capillary-pressure expression
441 : * -> 0. This depends on the turning-point saturation when drying became wetting, using the Land
442 : * equation
443 : * @param sgrmax maximum value possible for sgrdel. This will be equal to sgrdel if the
444 : * turning-point saturation is small
445 : * @param sldel value of the turning-point saturation when drying became wetting
446 : * @param m van-Genuchten m parameter
447 : * @param gamma index satisfying gamma > 0. Usually gamma = 1/3.
448 : * @param k_rg_max Maximum value of unextended gas relative permeability. If no low-saturation
449 : * extension is used then gas relative permeability = k_rg_max for sl <= slr. Satisfies 0 <
450 : * k_rg_max <= 1. Frequently k_rg_max = 1 is used
451 : * @param y0p Value of the derivative of the low-saturation extension at sl = slr. If an extension
452 : * is used then the gas relative permeability in the region sl <= slr is a cubic whose value is
453 : * unity at sl = 0, and derivative is zero at sl = 0
454 : */
455 : Real relativePermeabilityNWHys(Real sl,
456 : Real slr,
457 : Real sgrdel,
458 : Real sgrmax,
459 : Real sldel,
460 : Real m,
461 : Real gamma,
462 : Real k_rg_max,
463 : Real y0p);
464 :
465 : /**
466 : * Derivative of hysteretic relative permeability for gas with respect to the liquid saturation
467 : * @param sl liquid saturation
468 : * @param slr residual liquid saturation. For sl < slr, this function will always return 0
469 : * @param sgrdel value of gas saturation where van Genuchten wetting capillary-pressure expression
470 : * -> 0. This depends on the turning-point saturation when drying became wetting, using the Land
471 : * equation
472 : * @param sgrmax maximum value possible for sgrdel. This will be equal to sgrdel if the
473 : * turning-point saturation is small
474 : * @param sldel value of the turning-point saturation when drying became wetting
475 : * @param m van-Genuchten m parameter
476 : * @param gamma index satisfying gamma > 0. Usually gamma = 1/3.
477 : * @param k_rg_max Maximum value of unextended gas relative permeability. If no low-saturation
478 : * extension is used then gas relative permeability = k_rg_max for sl <= slr. Satisfies 0 <
479 : * k_rg_max <= 1. Frequently k_rg_max = 1 is used
480 : * @param y0p Value of the derivative of the low-saturation extension at sl = slr. If an extension
481 : * is used then the gas relative permeability in the region sl <= slr is a cubic whose value is
482 : * unity at sl = 0, and derivative is zero at sl = 0
483 : */
484 : Real drelativePermeabilityNWHys(Real sl,
485 : Real slr,
486 : Real sgrdel,
487 : Real sgrmax,
488 : Real sldel,
489 : Real m,
490 : Real gamma,
491 : Real k_rg_max,
492 : Real y0p);
493 : }
|
