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 : #include "PorousFlowVanGenuchten.h"
11 : #include "PorousFlowCubic.h"
12 :
13 : namespace PorousFlowVanGenuchten
14 : {
15 : Real
16 37083807 : effectiveSaturation(Real p, Real alpha, Real m)
17 : {
18 : Real n, seff;
19 :
20 37083807 : if (p >= 0.0)
21 : return 1.0;
22 : else
23 : {
24 4724315 : n = 1.0 / (1.0 - m);
25 4724315 : seff = 1.0 + std::pow(-alpha * p, n);
26 4724315 : return std::pow(seff, -m);
27 : }
28 : }
29 :
30 : Real
31 36359968 : dEffectiveSaturation(Real p, Real alpha, Real m)
32 : {
33 36359968 : if (p >= 0.0)
34 : return 0.0;
35 : else
36 : {
37 4677035 : Real n = 1.0 / (1.0 - m);
38 4677035 : Real inner = 1.0 + std::pow(-alpha * p, n);
39 4677035 : Real dinner_dp = -n * alpha * std::pow(-alpha * p, n - 1.0);
40 4677035 : Real dseff_dp = -m * std::pow(inner, -m - 1) * dinner_dp;
41 4677035 : return dseff_dp;
42 : }
43 : }
44 :
45 : Real
46 17995947 : d2EffectiveSaturation(Real p, Real alpha, Real m)
47 : {
48 17995947 : if (p >= 0.0)
49 : return 0.0;
50 : else
51 : {
52 2251580 : Real n = 1.0 / (1.0 - m);
53 2251580 : Real inner = 1.0 + std::pow(-alpha * p, n);
54 2251580 : Real dinner_dp = -n * alpha * std::pow(-alpha * p, n - 1.0);
55 2251580 : Real d2inner_dp2 = n * (n - 1.0) * alpha * alpha * std::pow(-alpha * p, n - 2.0);
56 2251580 : Real d2seff_dp2 = m * (m + 1.0) * std::pow(inner, -m - 2.0) * std::pow(dinner_dp, 2.0) -
57 2251580 : m * std::pow(inner, -m - 1.0) * d2inner_dp2;
58 2251580 : return d2seff_dp2;
59 : }
60 : }
61 :
62 : Real
63 768334 : capillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
64 : {
65 768334 : if (seff >= 1.0)
66 : return 0.0;
67 580522 : else if (seff <= 0.0)
68 27479 : return pc_max;
69 : else
70 : {
71 553043 : Real a = std::pow(seff, -1.0 / m) - 1.0;
72 553360 : return std::min(std::pow(a, 1.0 - m) / alpha, pc_max);
73 : }
74 : }
75 :
76 : Real
77 796925 : dCapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
78 : {
79 796925 : if (seff <= 0.0 || seff >= 1.0)
80 : return 0.0;
81 : else
82 : {
83 564457 : Real a = std::pow(seff, -1.0 / m) - 1.0;
84 : // Return 0 if pc > pc_max
85 564457 : if (std::pow(a, 1.0 - m) / alpha > pc_max)
86 : return 0.0;
87 : else
88 564282 : return (m - 1.0) * std::pow(a, -m) * std::pow(seff, -1.0 - 1.0 / m) / m / alpha;
89 : }
90 : }
91 :
92 : Real
93 404441 : d2CapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
94 : {
95 404441 : if (seff <= 0.0 || seff >= 1.0)
96 : return 0.0;
97 : else
98 : {
99 282850 : Real a = std::pow(seff, -1.0 / m) - 1.0;
100 : // Return 0 if pc > pc_max
101 282850 : if (std::pow(a, 1.0 - m) / alpha > pc_max)
102 : return 0.0;
103 : else
104 : {
105 282848 : Real d2pc = -std::pow(a, -1.0 - m) * std::pow(seff, -2.0 - 2.0 / m) +
106 282848 : ((1.0 + m) / m) * std::pow(a, -m) * std::pow(seff, -1.0 / m - 2.0);
107 282848 : d2pc *= (1.0 - m) / m / alpha;
108 282848 : return d2pc;
109 : }
110 : }
111 : }
112 :
113 : Real
114 751015 : dRelativePermeability(Real seff, Real m)
115 : {
116 : // Guard against division by zero
117 751015 : if (seff <= 0.0 || seff >= 1.0)
118 : return 0.0;
119 :
120 750449 : const Real a = 1.0 - std::pow(seff, 1.0 / m);
121 750449 : const Real da = -1.0 / m * std::pow(seff, 1.0 / m - 1.0);
122 750449 : const Real b = 1.0 - std::pow(a, m);
123 750449 : const Real db = -m * std::pow(a, m - 1.0) * da;
124 :
125 750449 : return 0.5 * std::pow(seff, -0.5) * Utility::pow<2>(b) + 2.0 * std::sqrt(seff) * b * db;
126 : }
127 :
128 : Real
129 103 : d2RelativePermeability(Real seff, Real m)
130 : {
131 : // Guard against division by zero
132 103 : if (seff <= 0.0 || seff >= 1.0)
133 : return 0.0;
134 :
135 101 : const Real a = 1.0 - std::pow(seff, 1.0 / m);
136 101 : const Real da = -1.0 / m * std::pow(seff, 1.0 / m - 1.0);
137 101 : const Real d2a = -(1.0 / m) * (1.0 / m - 1.0) * std::pow(seff, 1.0 / m - 2.0);
138 101 : const Real b = 1.0 - std::pow(a, m);
139 101 : const Real db = -m * std::pow(a, m - 1.0) * da;
140 101 : 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 101 : return -0.25 * std::pow(seff, -1.5) * Utility::pow<2>(b) + 2.0 * std::pow(seff, -0.5) * b * db +
143 101 : 2.0 * std::sqrt(seff) * db * db + 2.0 * std::sqrt(seff) * b * d2b;
144 : }
145 :
146 : Real
147 9321 : dRelativePermeabilityNW(Real seff, Real m)
148 : {
149 : // Guard against division by zero
150 9321 : if (seff <= 0.0 || seff >= 1.0)
151 : return 0.0;
152 :
153 9187 : const Real a = std::pow(1.0 - seff, 1.0 / m);
154 9187 : const Real da = -1.0 / m * a / (1.0 - seff);
155 9187 : const Real b = std::pow(1.0 - a, 2.0 * m);
156 9187 : const Real db = -2.0 * m * b / (1.0 - a) * da;
157 :
158 9187 : return 0.5 * std::pow(seff, -0.5) * b + std::sqrt(seff) * db;
159 : }
160 :
161 : Real
162 4 : d2RelativePermeabilityNW(Real seff, Real m)
163 : {
164 : // Guard against division by zero
165 4 : if (seff <= 0.0 || seff >= 1.0)
166 : return 0.0;
167 :
168 2 : const Real a = std::pow(1.0 - seff, 1.0 / m);
169 2 : const Real da = -1.0 / m * a / (1.0 - seff);
170 2 : const Real d2a = 1.0 / m * (1.0 / m - 1) * std::pow(1.0 - seff, 1.0 / m - 2.0);
171 2 : const Real b = std::pow(1.0 - a, 2.0 * m);
172 2 : const Real db = -2.0 * m * b / (1.0 - a) * da;
173 : const Real d2b =
174 2 : -2.0 * m * (db / (1.0 - a) * da + b * Utility::pow<2>(da / (1.0 - a)) + b / (1.0 - a) * d2a);
175 :
176 2 : return -0.25 * std::pow(seff, -1.5) * b + std::pow(seff, -0.5) * db + std::sqrt(seff) * d2b;
177 : }
178 :
179 : Real
180 3727977 : capillaryPressureHys(Real sl,
181 : Real slmin,
182 : Real sgrdel,
183 : Real alpha,
184 : Real n,
185 : const LowCapillaryPressureExtension & low_ext,
186 : const HighCapillaryPressureExtension & high_ext)
187 : {
188 : Real pc = 0.0;
189 3727977 : if (sl < low_ext.S) // important for initializing low_ext that this is < and not <=
190 : {
191 789083 : switch (low_ext.strategy)
192 : {
193 281179 : case LowCapillaryPressureExtension::QUADRATIC:
194 281179 : pc = low_ext.Pc + low_ext.dPc * 0.5 * (sl * sl - low_ext.S * low_ext.S) / low_ext.S;
195 281179 : break;
196 271567 : case LowCapillaryPressureExtension::EXPONENTIAL:
197 271567 : pc = low_ext.Pc * std::exp(low_ext.dPc * (sl - low_ext.S) / low_ext.Pc);
198 271567 : break;
199 236337 : default:
200 236337 : pc = low_ext.Pc;
201 : }
202 789083 : return pc;
203 : }
204 2938894 : if (sl > high_ext.S) // important for initializing high_ext that this is >, not >=
205 : {
206 607095 : switch (high_ext.strategy)
207 : {
208 375060 : case HighCapillaryPressureExtension::POWER:
209 : {
210 375060 : if (sl >= 1.0)
211 : pc = 0.0;
212 : else
213 : {
214 374766 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
215 374766 : 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 607095 : return pc;
223 : }
224 2331799 : const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
225 2331799 : if (seff >= 1.0)
226 : pc = 0.0; // no sensible high extension defined
227 2152357 : else if (seff <= 0.0)
228 5 : pc = low_ext.Pc; // no sensible low extension defined
229 : else
230 : {
231 2152352 : const Real a = std::pow(seff, n / (1.0 - n)) - 1.0;
232 2152352 : pc = (1.0 / alpha) * std::pow(a, 1.0 / n);
233 : }
234 : return pc;
235 : }
236 :
237 : Real
238 1377179 : dcapillaryPressureHys(Real sl,
239 : Real slmin,
240 : Real sgrdel,
241 : Real alpha,
242 : Real n,
243 : const LowCapillaryPressureExtension & low_ext,
244 : const HighCapillaryPressureExtension & high_ext)
245 : {
246 : Real dpc = 0.0;
247 1377179 : if (sl < low_ext.S) // important for initializing low_ext that this is < and not <=
248 : {
249 28749 : switch (low_ext.strategy)
250 : {
251 9689 : case LowCapillaryPressureExtension::QUADRATIC:
252 9689 : dpc = low_ext.dPc * sl / low_ext.S;
253 9689 : break;
254 9525 : case LowCapillaryPressureExtension::EXPONENTIAL:
255 9525 : dpc = low_ext.dPc * std::exp(low_ext.dPc * (sl - low_ext.S) / low_ext.Pc);
256 9525 : break;
257 : default:
258 : dpc = 0.0;
259 : }
260 28749 : return dpc;
261 : }
262 1348430 : if (sl > high_ext.S) // important for initializing high_ext that this is >, not >=
263 : {
264 83647 : switch (high_ext.strategy)
265 : {
266 65539 : case HighCapillaryPressureExtension::POWER:
267 : {
268 65539 : if (sl >= 1.0)
269 : dpc = 0.0;
270 : else
271 : {
272 65533 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
273 65533 : 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 83647 : return dpc;
281 : }
282 1264783 : const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
283 1264783 : if (seff >= 1.0)
284 : dpc = 0.0; // no sensible high extension defined
285 1085846 : else if (seff <= 0.0)
286 : dpc = 0.0; // no sensible low extension defined
287 : else
288 : {
289 1085839 : const Real a = std::pow(seff, n / (1.0 - n)) - 1.0;
290 1085839 : const Real dseff = 1.0 / (1.0 - sgrdel - slmin);
291 1085839 : const Real dpc_dseff = (1.0 / alpha / (1.0 - n)) * std::pow(a, 1.0 / n - 1.0) *
292 1085839 : std::pow(seff, n / (1.0 - n) - 1.0);
293 1085839 : dpc = dpc_dseff * dseff;
294 : }
295 : return dpc;
296 : }
297 :
298 : Real
299 147620 : d2capillaryPressureHys(Real sl,
300 : Real slmin,
301 : Real sgrdel,
302 : Real alpha,
303 : Real n,
304 : const LowCapillaryPressureExtension & low_ext,
305 : const HighCapillaryPressureExtension & high_ext)
306 : {
307 : Real d2pc = 0.0;
308 147620 : if (sl < low_ext.S) // important for initializing low_ext that this is < and not <=
309 : {
310 9611 : switch (low_ext.strategy)
311 : {
312 3257 : case LowCapillaryPressureExtension::QUADRATIC:
313 3257 : d2pc = low_ext.dPc / low_ext.S;
314 3257 : break;
315 3175 : case LowCapillaryPressureExtension::EXPONENTIAL:
316 3175 : d2pc = std::pow(low_ext.dPc, 2) / low_ext.Pc *
317 3175 : std::exp(low_ext.dPc * (sl - low_ext.S) / low_ext.Pc);
318 3175 : break;
319 : default:
320 : d2pc = 0.0;
321 : }
322 9611 : return d2pc;
323 : }
324 138009 : if (sl > high_ext.S) // important for initializing high_ext that this is >, not >=
325 : {
326 26324 : switch (high_ext.strategy)
327 : {
328 20260 : case HighCapillaryPressureExtension::POWER:
329 : {
330 20260 : if (sl >= 1.0)
331 : d2pc = 0.0;
332 : else
333 : {
334 20257 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
335 20257 : d2pc = high_ext.dPc * (1.0 - expon) / (1.0 - high_ext.S) *
336 20257 : 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 26324 : return d2pc;
344 : }
345 111685 : const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
346 111685 : if (seff >= 1.0)
347 : d2pc = 0.0; // no sensible high extension defined
348 111466 : else if (seff <= 0.0)
349 : d2pc = 0.0; // no sensible low extension defined
350 : else
351 : {
352 111463 : const Real a = std::pow(seff, n / (1.0 - n)) - 1.0;
353 111463 : const Real dseff = 1.0 / (1.0 - sgrdel - slmin);
354 : const Real d2pc_dseff =
355 111463 : (1.0 / alpha / (1.0 - n)) *
356 111463 : (std::pow(a, 1.0 / n - 2.0) * std::pow(seff, 2.0 * (n / (1.0 - n) - 1.0)) +
357 111463 : (n / (1.0 - n) - 1.0) * std::pow(a, 1.0 / n - 1.0) * std::pow(seff, n / (1.0 - n) - 2.0));
358 111463 : d2pc = d2pc_dseff * dseff * dseff;
359 : }
360 : return d2pc;
361 : }
362 :
363 : Real
364 2447793 : saturationHys(Real pc,
365 : Real slmin,
366 : Real sgrdel,
367 : Real alpha,
368 : Real n,
369 : const LowCapillaryPressureExtension & low_ext,
370 : const HighCapillaryPressureExtension & high_ext)
371 : {
372 2447793 : if (pc <= 0)
373 : return 1.0;
374 : Real s = 1.0;
375 2327879 : if (pc > low_ext.Pc) // important for initialization of the low_ext that this is > and not >=
376 : {
377 520530 : switch (low_ext.strategy)
378 : {
379 261644 : case LowCapillaryPressureExtension::QUADRATIC:
380 : {
381 261644 : const Real s2 = low_ext.S * low_ext.S + 2.0 * (pc - low_ext.Pc) * low_ext.S / low_ext.dPc;
382 261644 : 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 225078 : s = std::sqrt(s2);
388 : break;
389 : }
390 245154 : case LowCapillaryPressureExtension::EXPONENTIAL:
391 : {
392 245154 : const Real ss = low_ext.S + std::log(pc / low_ext.Pc) * low_ext.Pc / low_ext.dPc;
393 245154 : 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 13732 : default:
402 13732 : s = low_ext.S;
403 : }
404 520530 : return s;
405 : }
406 1807349 : if (pc < high_ext.Pc) // important for initialization of the high_ext that this is < and not <=
407 : {
408 178170 : switch (high_ext.strategy)
409 : {
410 178165 : case HighCapillaryPressureExtension::POWER:
411 : {
412 178165 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
413 178165 : s = 1.0 - std::pow(pc / high_ext.Pc, 1.0 / expon) * (1.0 - high_ext.S);
414 178165 : break;
415 : }
416 5 : default:
417 5 : s = high_ext.S;
418 : }
419 178170 : return s;
420 : }
421 1629179 : if (pc == std::numeric_limits<Real>::max())
422 : s = 0.0;
423 : else
424 : {
425 1629178 : const Real seffpow = 1.0 + std::pow(pc * alpha, n);
426 1629178 : const Real seff = std::pow(seffpow, (1.0 - n) / n);
427 1629178 : s = (1.0 - sgrdel - slmin) * seff + slmin;
428 : }
429 : return s;
430 : }
431 :
432 : Real
433 626542 : dsaturationHys(Real pc,
434 : Real slmin,
435 : Real sgrdel,
436 : Real alpha,
437 : Real n,
438 : const LowCapillaryPressureExtension & low_ext,
439 : const HighCapillaryPressureExtension & high_ext)
440 : {
441 626542 : if (pc <= 0)
442 : return 0.0;
443 : Real ds = 0.0;
444 594804 : if (pc > low_ext.Pc) // important for initialization of the low_ext that this is > and not >=
445 : {
446 107461 : switch (low_ext.strategy)
447 : {
448 51818 : case LowCapillaryPressureExtension::QUADRATIC:
449 : {
450 51818 : const Real s2 = low_ext.S * low_ext.S + 2.0 * (pc - low_ext.Pc) * low_ext.S / low_ext.dPc;
451 51818 : 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 49764 : const Real ds2 = 2.0 * low_ext.S / low_ext.dPc;
458 49764 : ds = 0.5 * ds2 / std::sqrt(s2);
459 : }
460 : break;
461 : }
462 46960 : case LowCapillaryPressureExtension::EXPONENTIAL:
463 : {
464 46960 : const Real s = low_ext.S + std::log(pc / low_ext.Pc) * low_ext.Pc / low_ext.dPc;
465 46960 : 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 45130 : ds = low_ext.Pc / pc / low_ext.dPc;
471 : break;
472 : }
473 : default:
474 : ds = 0.0;
475 : }
476 107461 : return ds;
477 : }
478 487343 : if (pc < high_ext.Pc) // important for initialization of the high_ext that this is < and not <=
479 : {
480 75043 : switch (high_ext.strategy)
481 : {
482 75037 : case HighCapillaryPressureExtension::POWER:
483 : {
484 75037 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
485 75037 : ds = -(1.0 - high_ext.S) / pc / expon * std::pow(pc / high_ext.Pc, 1.0 / expon);
486 75037 : break;
487 : }
488 : default:
489 : ds = 0.0;
490 : }
491 75043 : return ds;
492 : }
493 412300 : if (pc == std::numeric_limits<Real>::max())
494 : ds = 0.0;
495 : else
496 : {
497 412300 : const Real seffpow = 1.0 + std::pow(pc * alpha, n);
498 412300 : const Real dseffpow = n * (seffpow - 1.0) / pc;
499 412300 : const Real seff = std::pow(seffpow, (1.0 - n) / n);
500 412300 : const Real dseff = (1.0 - n) / n * seff / seffpow * dseffpow;
501 412300 : ds = (1.0 - sgrdel - slmin) * dseff;
502 : }
503 : return ds;
504 : }
505 :
506 : Real
507 230464 : d2saturationHys(Real pc,
508 : Real slmin,
509 : Real sgrdel,
510 : Real alpha,
511 : Real n,
512 : const LowCapillaryPressureExtension & low_ext,
513 : const HighCapillaryPressureExtension & high_ext)
514 : {
515 230464 : if (pc <= 0)
516 : return 0.0;
517 : Real d2s = 0.0;
518 216566 : if (pc > low_ext.Pc) // important for initialization of the low_ext that this is > and not >=
519 : {
520 38075 : switch (low_ext.strategy)
521 : {
522 18554 : case LowCapillaryPressureExtension::QUADRATIC:
523 : {
524 18554 : const Real s2 = low_ext.S * low_ext.S + 2.0 * (pc - low_ext.Pc) * low_ext.S / low_ext.dPc;
525 18554 : 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 17590 : const Real ds2 = 2.0 * low_ext.S / low_ext.dPc;
532 17590 : d2s = -0.25 * ds2 * ds2 / std::pow(s2, 1.5);
533 : }
534 : break;
535 : }
536 17918 : case LowCapillaryPressureExtension::EXPONENTIAL:
537 : {
538 17918 : const Real s = low_ext.S + std::log(pc / low_ext.Pc) * low_ext.Pc / low_ext.dPc;
539 17918 : 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 17308 : d2s = -low_ext.Pc / std::pow(pc, 2.0) / low_ext.dPc;
545 : break;
546 : }
547 : default:
548 : d2s = 0.0;
549 : }
550 38075 : return d2s;
551 : }
552 178491 : if (pc < high_ext.Pc) // important for initialization of the high_ext that this is < and not <=
553 : {
554 26989 : switch (high_ext.strategy)
555 : {
556 26987 : case HighCapillaryPressureExtension::POWER:
557 : {
558 26987 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
559 26987 : d2s = -(1.0 - high_ext.S) * (1.0 / expon) * (1.0 / expon - 1.0) /
560 26987 : std::pow(high_ext.Pc, 2.0) * std::pow(pc / high_ext.Pc, 1.0 / expon - 2.0);
561 26987 : break;
562 : }
563 : default:
564 : d2s = 0.0;
565 : }
566 26989 : return d2s;
567 : }
568 151502 : if (pc == std::numeric_limits<Real>::max())
569 : d2s = 0.0;
570 : else
571 : {
572 151502 : const Real seffpow = 1.0 + std::pow(pc * alpha, n);
573 151502 : const Real dseffpow = n * (seffpow - 1.0) / pc;
574 151502 : const Real d2seffpow = (n - 1.0) * dseffpow / pc;
575 151502 : const Real seff = std::pow(seffpow, (1.0 - n) / n);
576 151502 : const Real dseff = (1.0 - n) / n * seff / seffpow * dseffpow;
577 151502 : const Real d2seff =
578 151502 : (1.0 - n) / n *
579 151502 : (dseff * dseffpow - seff * dseffpow * dseffpow / seffpow + seff * d2seffpow) / seffpow;
580 151502 : d2s = (1.0 - sgrdel - slmin) * d2seff;
581 : }
582 : return d2s;
583 : }
584 :
585 : Real
586 42778 : relativePermeabilityHys(Real sl,
587 : Real slr,
588 : Real sgrdel,
589 : Real sgrmax,
590 : Real sldel,
591 : Real m,
592 : Real upper_liquid_param,
593 : Real y0,
594 : Real y0p,
595 : Real y1,
596 : Real y1p)
597 : {
598 42778 : if (sl <= slr) // by the definition of slr, always return 0
599 : return 0.0;
600 33292 : 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 33292 : if (sgrdel == 0.0 || sl <= sldel) // along the drying curve
606 14159 : 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 19133 : const Real my_sldel = (sldel < slr) ? slr : sldel;
613 19133 : const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
614 19133 : 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 3192 : a = std::pow(1.0 - std::pow(sl_bar, 1.0 / m), m);
619 : }
620 15941 : 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 2760 : return PorousFlowCubic::cubic(
625 2760 : 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 13181 : const Real sl_bar_del = (my_sldel - slr) / (1.0 - slr);
631 13181 : const Real s_gt_bar =
632 13181 : my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
633 13181 : a = (1 - s_gt_bar / (1.0 - sl_bar_del)) *
634 13181 : std::pow(1.0 - std::pow(sl_bar + s_gt_bar, 1.0 / m), m);
635 13181 : 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 30532 : return std::sqrt(sl_bar) * Utility::pow<2>(1.0 - a - b);
639 : }
640 :
641 : Real
642 42753 : drelativePermeabilityHys(Real sl,
643 : Real slr,
644 : Real sgrdel,
645 : Real sgrmax,
646 : Real sldel,
647 : Real m,
648 : Real upper_liquid_param,
649 : Real y0,
650 : Real y0p,
651 : Real y1,
652 : Real y1p)
653 : {
654 42753 : if (sl <= slr) // by the definition of slr, always return 0
655 : return 0.0;
656 33269 : if (sl == 1.0) // derivative is infinite at this point
657 : return std::numeric_limits<Real>::max();
658 31505 : const Real sl_bar = (sl - slr) / (1.0 - slr); // effective saturation
659 31505 : 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 31505 : if (sgrdel == 0.0 || sl <= sldel) // along the drying curve
667 : {
668 13106 : const Real c = std::pow(sl_bar, 1.0 / m);
669 13106 : const Real dc_dsbar = c / m / sl_bar;
670 13106 : a = std::pow(1.0 - c, m);
671 13106 : const Real da_dsbar = -m * a / (1.0 - c) * dc_dsbar;
672 13106 : a_prime = da_dsbar * sl_bar_prime;
673 13106 : }
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 18399 : const Real my_sldel = (sldel < slr) ? slr : sldel;
680 18399 : const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
681 18399 : 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 2468 : const Real c = std::pow(sl_bar, 1.0 / m);
686 2468 : const Real dc_dsbar = c / m / sl_bar;
687 2468 : a = std::pow(1.0 - c, m);
688 2468 : const Real da_dsbar = -m * a / (1.0 - c) * dc_dsbar;
689 2468 : a_prime = da_dsbar * sl_bar_prime;
690 : }
691 15931 : 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
695 2756 : return PorousFlowCubic::dcubic(
696 2756 : 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 13175 : const Real sl_bar_del = (my_sldel - slr) / (1.0 - slr);
702 13175 : const Real s_gt_bar =
703 13175 : my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
704 13175 : const Real s_gt_bar_prime = my_sgrdel / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
705 13175 : const Real c = std::pow(sl_bar + s_gt_bar, 1.0 / m);
706 13175 : const Real c_prime = c / m / (sl_bar + s_gt_bar) * (sl_bar_prime + s_gt_bar_prime);
707 13175 : a = (1 - s_gt_bar / (1.0 - sl_bar_del)) * std::pow(1.0 - c, m);
708 13175 : a_prime =
709 13175 : -s_gt_bar_prime / (1.0 - sl_bar_del) * std::pow(1.0 - c, m) - m * a / (1.0 - c) * c_prime;
710 13175 : b = s_gt_bar / (1.0 - sl_bar_del) * std::pow(1.0 - std::pow(sl_bar_del, 1.0 / m), m);
711 13175 : b_prime = s_gt_bar_prime * b / s_gt_bar;
712 : }
713 : }
714 28749 : const Real kr = std::sqrt(sl_bar) * Utility::pow<2>(1.0 - a - b);
715 28749 : return 0.5 * kr / sl_bar * sl_bar_prime -
716 28749 : std::sqrt(sl_bar) * 2.0 * (1.0 - a - b) * (a_prime + b_prime);
717 : }
718 :
719 : Real
720 23791 : relativePermeabilityNWHys(Real sl,
721 : Real slr,
722 : Real sgrdel,
723 : Real sgrmax,
724 : Real sldel,
725 : Real m,
726 : Real gamma,
727 : Real k_rg_max,
728 : Real y0p)
729 : {
730 23791 : if (sl < slr)
731 : {
732 : // in the extended region, so immediately return with the relevant value
733 6040 : if (k_rg_max == 1.0)
734 : return 1.0;
735 4082 : return PorousFlowCubic::cubic(sl, 0.0, 1.0, 0.0, slr, k_rg_max, y0p);
736 : }
737 17751 : if (sl > 1.0 - sgrdel) // saturation is above 1.0 - residual gas saturation
738 : return 0.0;
739 16992 : 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 16992 : 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 8570 : const Real my_sldel = (sldel < slr) ? slr : sldel;
748 8570 : const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
749 8570 : s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
750 : }
751 : Real kr = 0.0;
752 16992 : 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 16415 : const Real a = std::pow(1.0 - (sl_bar + s_gt_bar), gamma);
756 16415 : const Real c = std::pow(sl_bar + s_gt_bar, 1.0 / m);
757 16415 : const Real b = std::pow(1.0 - c, 2.0 * m);
758 16415 : kr = k_rg_max * a * b;
759 : }
760 : return kr;
761 : }
762 :
763 : Real
764 23833 : drelativePermeabilityNWHys(Real sl,
765 : Real slr,
766 : Real sgrdel,
767 : Real sgrmax,
768 : Real sldel,
769 : Real m,
770 : Real gamma,
771 : Real k_rg_max,
772 : Real y0p)
773 : {
774 23833 : if (sl < slr)
775 : {
776 : // in the extended region, so immediately return with the relevant value
777 6036 : if (k_rg_max == 1.0)
778 : return 0.0;
779 4078 : return PorousFlowCubic::dcubic(sl, 0.0, 1.0, 0.0, slr, k_rg_max, y0p);
780 : }
781 17797 : if (sl > 1.0 - sgrdel) // saturation is above 1.0 - residual gas saturation
782 : return 0.0;
783 17038 : const Real sl_bar = (sl - slr) / (1.0 - slr);
784 17038 : 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 17038 : 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 8561 : const Real my_sldel = (sldel < slr) ? slr : sldel;
794 8561 : const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
795 8561 : s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
796 8561 : s_gt_bar_prime = my_sgrdel / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
797 : }
798 : Real kr_prime = 0.0;
799 17038 : 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 16461 : const Real a = std::pow(1.0 - (sl_bar + s_gt_bar), gamma);
803 16461 : const Real a_prime = -gamma * a / (1.0 - (sl_bar + s_gt_bar)) * (sl_bar_prime + s_gt_bar_prime);
804 16461 : const Real c = std::pow(sl_bar + s_gt_bar, 1.0 / m);
805 : const Real c_prime =
806 16461 : (c == 0 ? 0.0 : c / m / (sl_bar + s_gt_bar) * (sl_bar_prime + s_gt_bar_prime));
807 16461 : const Real b = std::pow(1.0 - c, 2.0 * m);
808 16461 : const Real b_prime = -2.0 * m * b / (1.0 - c) * c_prime;
809 16461 : kr_prime = k_rg_max * (a * b_prime + a_prime * b);
810 : }
811 : return kr_prime;
812 : }
813 : }
|
