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 25298938 : effectiveSaturation(Real p, Real alpha, Real m)
17 : {
18 : Real n, seff;
19 :
20 25298938 : if (p >= 0.0)
21 : return 1.0;
22 : else
23 : {
24 3840613 : n = 1.0 / (1.0 - m);
25 3840613 : seff = 1.0 + std::pow(-alpha * p, n);
26 3840613 : return std::pow(seff, -m);
27 : }
28 : }
29 :
30 : Real
31 24866690 : dEffectiveSaturation(Real p, Real alpha, Real m)
32 : {
33 24866690 : if (p >= 0.0)
34 : return 0.0;
35 : else
36 : {
37 3811109 : Real n = 1.0 / (1.0 - m);
38 3811109 : Real inner = 1.0 + std::pow(-alpha * p, n);
39 3811109 : Real dinner_dp = -n * alpha * std::pow(-alpha * p, n - 1.0);
40 3811109 : Real dseff_dp = -m * std::pow(inner, -m - 1) * dinner_dp;
41 3811109 : return dseff_dp;
42 : }
43 : }
44 :
45 : Real
46 12310024 : d2EffectiveSaturation(Real p, Real alpha, Real m)
47 : {
48 12310024 : if (p >= 0.0)
49 : return 0.0;
50 : else
51 : {
52 1835268 : Real n = 1.0 / (1.0 - m);
53 1835268 : Real inner = 1.0 + std::pow(-alpha * p, n);
54 1835268 : Real dinner_dp = -n * alpha * std::pow(-alpha * p, n - 1.0);
55 1835268 : Real d2inner_dp2 = n * (n - 1.0) * alpha * alpha * std::pow(-alpha * p, n - 2.0);
56 1835268 : Real d2seff_dp2 = m * (m + 1.0) * std::pow(inner, -m - 2.0) * std::pow(dinner_dp, 2.0) -
57 1835268 : m * std::pow(inner, -m - 1.0) * d2inner_dp2;
58 1835268 : return d2seff_dp2;
59 : }
60 : }
61 :
62 : Real
63 587788 : capillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
64 : {
65 587788 : if (seff >= 1.0)
66 : return 0.0;
67 409022 : else if (seff <= 0.0)
68 19441 : return pc_max;
69 : else
70 : {
71 389581 : Real a = std::pow(seff, -1.0 / m) - 1.0;
72 389737 : return std::min(std::pow(a, 1.0 - m) / alpha, pc_max);
73 : }
74 : }
75 :
76 : Real
77 599332 : dCapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
78 : {
79 599332 : if (seff <= 0.0 || seff >= 1.0)
80 : return 0.0;
81 : else
82 : {
83 399695 : Real a = std::pow(seff, -1.0 / m) - 1.0;
84 : // Return 0 if pc > pc_max
85 399695 : if (std::pow(a, 1.0 - m) / alpha > pc_max)
86 : return 0.0;
87 : else
88 399609 : return (m - 1.0) * std::pow(a, -m) * std::pow(seff, -1.0 - 1.0 / m) / m / alpha;
89 : }
90 : }
91 :
92 : Real
93 292105 : d2CapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
94 : {
95 292105 : if (seff <= 0.0 || seff >= 1.0)
96 : return 0.0;
97 : else
98 : {
99 195074 : Real a = std::pow(seff, -1.0 / m) - 1.0;
100 : // Return 0 if pc > pc_max
101 195074 : if (std::pow(a, 1.0 - m) / alpha > pc_max)
102 : return 0.0;
103 : else
104 : {
105 195073 : Real d2pc = -std::pow(a, -1.0 - m) * std::pow(seff, -2.0 - 2.0 / m) +
106 195073 : ((1.0 + m) / m) * std::pow(a, -m) * std::pow(seff, -1.0 / m - 2.0);
107 195073 : d2pc *= (1.0 - m) / m / alpha;
108 195073 : return d2pc;
109 : }
110 : }
111 : }
112 :
113 : Real
114 502013 : dRelativePermeability(Real seff, Real m)
115 : {
116 : // Guard against division by zero
117 502013 : if (seff <= 0.0 || seff >= 1.0)
118 : return 0.0;
119 :
120 501747 : const Real a = 1.0 - std::pow(seff, 1.0 / m);
121 501747 : const Real da = -1.0 / m * std::pow(seff, 1.0 / m - 1.0);
122 501747 : const Real b = 1.0 - std::pow(a, m);
123 501747 : const Real db = -m * std::pow(a, m - 1.0) * da;
124 :
125 501747 : return 0.5 * std::pow(seff, -0.5) * Utility::pow<2>(b) + 2.0 * std::sqrt(seff) * b * db;
126 : }
127 :
128 : Real
129 55 : d2RelativePermeability(Real seff, Real m)
130 : {
131 : // Guard against division by zero
132 55 : if (seff <= 0.0 || seff >= 1.0)
133 : return 0.0;
134 :
135 53 : const Real a = 1.0 - std::pow(seff, 1.0 / m);
136 53 : const Real da = -1.0 / m * std::pow(seff, 1.0 / m - 1.0);
137 53 : const Real d2a = -(1.0 / m) * (1.0 / m - 1.0) * std::pow(seff, 1.0 / m - 2.0);
138 53 : const Real b = 1.0 - std::pow(a, m);
139 53 : const Real db = -m * std::pow(a, m - 1.0) * da;
140 53 : 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 53 : return -0.25 * std::pow(seff, -1.5) * Utility::pow<2>(b) + 2.0 * std::pow(seff, -0.5) * b * db +
143 53 : 2.0 * std::sqrt(seff) * db * db + 2.0 * std::sqrt(seff) * b * d2b;
144 : }
145 :
146 : Real
147 6189 : dRelativePermeabilityNW(Real seff, Real m)
148 : {
149 : // Guard against division by zero
150 6189 : if (seff <= 0.0 || seff >= 1.0)
151 : return 0.0;
152 :
153 6127 : const Real a = std::pow(1.0 - seff, 1.0 / m);
154 6127 : const Real da = -1.0 / m * a / (1.0 - seff);
155 6127 : const Real b = std::pow(1.0 - a, 2.0 * m);
156 6127 : const Real db = -2.0 * m * b / (1.0 - a) * da;
157 :
158 6127 : 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 2434705 : 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 2434705 : if (sl < low_ext.S) // important for initializing low_ext that this is < and not <=
190 : {
191 505921 : switch (low_ext.strategy)
192 : {
193 180497 : case LowCapillaryPressureExtension::QUADRATIC:
194 180497 : pc = low_ext.Pc + low_ext.dPc * 0.5 * (sl * sl - low_ext.S * low_ext.S) / low_ext.S;
195 180497 : break;
196 174327 : case LowCapillaryPressureExtension::EXPONENTIAL:
197 174327 : pc = low_ext.Pc * std::exp(low_ext.dPc * (sl - low_ext.S) / low_ext.Pc);
198 174327 : break;
199 151097 : default:
200 151097 : pc = low_ext.Pc;
201 : }
202 505921 : return pc;
203 : }
204 1928784 : if (sl > high_ext.S) // important for initializing high_ext that this is >, not >=
205 : {
206 401239 : switch (high_ext.strategy)
207 : {
208 248012 : case HighCapillaryPressureExtension::POWER:
209 : {
210 248012 : if (sl >= 1.0)
211 : pc = 0.0;
212 : else
213 : {
214 247814 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
215 247814 : 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 401239 : return pc;
223 : }
224 1527545 : const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
225 1527545 : if (seff >= 1.0)
226 : pc = 0.0; // no sensible high extension defined
227 1416415 : else if (seff <= 0.0)
228 5 : pc = low_ext.Pc; // no sensible low extension defined
229 : else
230 : {
231 1416410 : const Real a = std::pow(seff, n / (1.0 - n)) - 1.0;
232 1416410 : pc = (1.0 / alpha) * std::pow(a, 1.0 / n);
233 : }
234 : return pc;
235 : }
236 :
237 : Real
238 874395 : 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 874395 : if (sl < low_ext.S) // important for initializing low_ext that this is < and not <=
248 : {
249 19389 : switch (low_ext.strategy)
250 : {
251 6569 : case LowCapillaryPressureExtension::QUADRATIC:
252 6569 : dpc = low_ext.dPc * sl / low_ext.S;
253 6569 : break;
254 6405 : case LowCapillaryPressureExtension::EXPONENTIAL:
255 6405 : dpc = low_ext.dPc * std::exp(low_ext.dPc * (sl - low_ext.S) / low_ext.Pc);
256 6405 : break;
257 : default:
258 : dpc = 0.0;
259 : }
260 19389 : return dpc;
261 : }
262 855006 : if (sl > high_ext.S) // important for initializing high_ext that this is >, not >=
263 : {
264 56279 : switch (high_ext.strategy)
265 : {
266 44051 : case HighCapillaryPressureExtension::POWER:
267 : {
268 44051 : if (sl >= 1.0)
269 : dpc = 0.0;
270 : else
271 : {
272 44045 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
273 44045 : 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 56279 : return dpc;
281 : }
282 798727 : const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
283 798727 : if (seff >= 1.0)
284 : dpc = 0.0; // no sensible high extension defined
285 687934 : else if (seff <= 0.0)
286 : dpc = 0.0; // no sensible low extension defined
287 : else
288 : {
289 687927 : const Real a = std::pow(seff, n / (1.0 - n)) - 1.0;
290 687927 : const Real dseff = 1.0 / (1.0 - sgrdel - slmin);
291 687927 : const Real dpc_dseff = (1.0 / alpha / (1.0 - n)) * std::pow(a, 1.0 / n - 1.0) *
292 687927 : std::pow(seff, n / (1.0 - n) - 1.0);
293 687927 : dpc = dpc_dseff * dseff;
294 : }
295 : return dpc;
296 : }
297 :
298 : Real
299 99672 : 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 99672 : if (sl < low_ext.S) // important for initializing low_ext that this is < and not <=
309 : {
310 6491 : switch (low_ext.strategy)
311 : {
312 2217 : case LowCapillaryPressureExtension::QUADRATIC:
313 2217 : d2pc = low_ext.dPc / low_ext.S;
314 2217 : break;
315 2135 : case LowCapillaryPressureExtension::EXPONENTIAL:
316 2135 : d2pc = std::pow(low_ext.dPc, 2) / low_ext.Pc *
317 2135 : std::exp(low_ext.dPc * (sl - low_ext.S) / low_ext.Pc);
318 2135 : break;
319 : default:
320 : d2pc = 0.0;
321 : }
322 6491 : return d2pc;
323 : }
324 93181 : if (sl > high_ext.S) // important for initializing high_ext that this is >, not >=
325 : {
326 17720 : switch (high_ext.strategy)
327 : {
328 13616 : case HighCapillaryPressureExtension::POWER:
329 : {
330 13616 : if (sl >= 1.0)
331 : d2pc = 0.0;
332 : else
333 : {
334 13613 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
335 13613 : d2pc = high_ext.dPc * (1.0 - expon) / (1.0 - high_ext.S) *
336 13613 : 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 17720 : return d2pc;
344 : }
345 75461 : const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
346 75461 : if (seff >= 1.0)
347 : d2pc = 0.0; // no sensible high extension defined
348 75314 : else if (seff <= 0.0)
349 : d2pc = 0.0; // no sensible low extension defined
350 : else
351 : {
352 75311 : const Real a = std::pow(seff, n / (1.0 - n)) - 1.0;
353 75311 : const Real dseff = 1.0 / (1.0 - sgrdel - slmin);
354 : const Real d2pc_dseff =
355 75311 : (1.0 / alpha / (1.0 - n)) *
356 75311 : (std::pow(a, 1.0 / n - 2.0) * std::pow(seff, 2.0 * (n / (1.0 - n) - 1.0)) +
357 75311 : (n / (1.0 - n) - 1.0) * std::pow(a, 1.0 / n - 1.0) * std::pow(seff, n / (1.0 - n) - 2.0));
358 75311 : d2pc = d2pc_dseff * dseff * dseff;
359 : }
360 : return d2pc;
361 : }
362 :
363 : Real
364 1570545 : 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 1570545 : if (pc <= 0)
373 : return 1.0;
374 : Real s = 1.0;
375 1493091 : if (pc > low_ext.Pc) // important for initialization of the low_ext that this is > and not >=
376 : {
377 333206 : switch (low_ext.strategy)
378 : {
379 167400 : case LowCapillaryPressureExtension::QUADRATIC:
380 : {
381 167400 : const Real s2 = low_ext.S * low_ext.S + 2.0 * (pc - low_ext.Pc) * low_ext.S / low_ext.dPc;
382 167400 : 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 144608 : s = std::sqrt(s2);
388 : break;
389 : }
390 156574 : case LowCapillaryPressureExtension::EXPONENTIAL:
391 : {
392 156574 : const Real ss = low_ext.S + std::log(pc / low_ext.Pc) * low_ext.Pc / low_ext.dPc;
393 156574 : 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 9232 : default:
402 9232 : s = low_ext.S;
403 : }
404 333206 : return s;
405 : }
406 1159885 : if (pc < high_ext.Pc) // important for initialization of the high_ext that this is < and not <=
407 : {
408 115586 : switch (high_ext.strategy)
409 : {
410 115581 : case HighCapillaryPressureExtension::POWER:
411 : {
412 115581 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
413 115581 : s = 1.0 - std::pow(pc / high_ext.Pc, 1.0 / expon) * (1.0 - high_ext.S);
414 115581 : break;
415 : }
416 5 : default:
417 5 : s = high_ext.S;
418 : }
419 115586 : return s;
420 : }
421 1044299 : if (pc == std::numeric_limits<Real>::max())
422 : s = 0.0;
423 : else
424 : {
425 1044298 : const Real seffpow = 1.0 + std::pow(pc * alpha, n);
426 1044298 : const Real seff = std::pow(seffpow, (1.0 - n) / n);
427 1044298 : s = (1.0 - sgrdel - slmin) * seff + slmin;
428 : }
429 : return s;
430 : }
431 :
432 : Real
433 421750 : 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 421750 : if (pc <= 0)
442 : return 0.0;
443 : Real ds = 0.0;
444 400700 : if (pc > low_ext.Pc) // important for initialization of the low_ext that this is > and not >=
445 : {
446 72525 : switch (low_ext.strategy)
447 : {
448 34942 : case LowCapillaryPressureExtension::QUADRATIC:
449 : {
450 34942 : const Real s2 = low_ext.S * low_ext.S + 2.0 * (pc - low_ext.Pc) * low_ext.S / low_ext.dPc;
451 34942 : 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 33518 : const Real ds2 = 2.0 * low_ext.S / low_ext.dPc;
458 33518 : ds = 0.5 * ds2 / std::sqrt(s2);
459 : }
460 : break;
461 : }
462 31700 : case LowCapillaryPressureExtension::EXPONENTIAL:
463 : {
464 31700 : const Real s = low_ext.S + std::log(pc / low_ext.Pc) * low_ext.Pc / low_ext.dPc;
465 31700 : 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 30470 : ds = low_ext.Pc / pc / low_ext.dPc;
471 : break;
472 : }
473 : default:
474 : ds = 0.0;
475 : }
476 72525 : return ds;
477 : }
478 328175 : if (pc < high_ext.Pc) // important for initialization of the high_ext that this is < and not <=
479 : {
480 50487 : switch (high_ext.strategy)
481 : {
482 50481 : case HighCapillaryPressureExtension::POWER:
483 : {
484 50481 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
485 50481 : ds = -(1.0 - high_ext.S) / pc / expon * std::pow(pc / high_ext.Pc, 1.0 / expon);
486 50481 : break;
487 : }
488 : default:
489 : ds = 0.0;
490 : }
491 50487 : return ds;
492 : }
493 277688 : if (pc == std::numeric_limits<Real>::max())
494 : ds = 0.0;
495 : else
496 : {
497 277688 : const Real seffpow = 1.0 + std::pow(pc * alpha, n);
498 277688 : const Real dseffpow = n * (seffpow - 1.0) / pc;
499 277688 : const Real seff = std::pow(seffpow, (1.0 - n) / n);
500 277688 : const Real dseff = (1.0 - n) / n * seff / seffpow * dseffpow;
501 277688 : ds = (1.0 - sgrdel - slmin) * dseff;
502 : }
503 : return ds;
504 : }
505 :
506 : Real
507 155608 : 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 155608 : if (pc <= 0)
516 : return 0.0;
517 : Real d2s = 0.0;
518 146354 : if (pc > low_ext.Pc) // important for initialization of the low_ext that this is > and not >=
519 : {
520 25849 : switch (low_ext.strategy)
521 : {
522 12548 : case LowCapillaryPressureExtension::QUADRATIC:
523 : {
524 12548 : const Real s2 = low_ext.S * low_ext.S + 2.0 * (pc - low_ext.Pc) * low_ext.S / low_ext.dPc;
525 12548 : 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 11878 : const Real ds2 = 2.0 * low_ext.S / low_ext.dPc;
532 11878 : d2s = -0.25 * ds2 * ds2 / std::pow(s2, 1.5);
533 : }
534 : break;
535 : }
536 12178 : case LowCapillaryPressureExtension::EXPONENTIAL:
537 : {
538 12178 : const Real s = low_ext.S + std::log(pc / low_ext.Pc) * low_ext.Pc / low_ext.dPc;
539 12178 : 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 11768 : d2s = -low_ext.Pc / std::pow(pc, 2.0) / low_ext.dPc;
545 : break;
546 : }
547 : default:
548 : d2s = 0.0;
549 : }
550 25849 : return d2s;
551 : }
552 120505 : if (pc < high_ext.Pc) // important for initialization of the high_ext that this is < and not <=
553 : {
554 18093 : switch (high_ext.strategy)
555 : {
556 18091 : case HighCapillaryPressureExtension::POWER:
557 : {
558 18091 : const Real expon = -high_ext.dPc / high_ext.Pc * (1.0 - high_ext.S);
559 18091 : d2s = -(1.0 - high_ext.S) * (1.0 / expon) * (1.0 / expon - 1.0) /
560 18091 : std::pow(high_ext.Pc, 2.0) * std::pow(pc / high_ext.Pc, 1.0 / expon - 2.0);
561 18091 : break;
562 : }
563 : default:
564 : d2s = 0.0;
565 : }
566 18093 : return d2s;
567 : }
568 102412 : if (pc == std::numeric_limits<Real>::max())
569 : d2s = 0.0;
570 : else
571 : {
572 102412 : const Real seffpow = 1.0 + std::pow(pc * alpha, n);
573 102412 : const Real dseffpow = n * (seffpow - 1.0) / pc;
574 102412 : const Real d2seffpow = (n - 1.0) * dseffpow / pc;
575 102412 : const Real seff = std::pow(seffpow, (1.0 - n) / n);
576 102412 : const Real dseff = (1.0 - n) / n * seff / seffpow * dseffpow;
577 102412 : const Real d2seff =
578 102412 : (1.0 - n) / n *
579 102412 : (dseff * dseffpow - seff * dseffpow * dseffpow / seffpow + seff * d2seffpow) / seffpow;
580 102412 : d2s = (1.0 - sgrdel - slmin) * d2seff;
581 : }
582 : return d2s;
583 : }
584 :
585 : Real
586 30124 : 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 30124 : if (sl <= slr) // by the definition of slr, always return 0
599 : return 0.0;
600 23258 : 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 23258 : if (sgrdel == 0.0 || sl <= sldel) // along the drying curve
606 9927 : 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 13331 : const Real my_sldel = (sldel < slr) ? slr : sldel;
613 13331 : const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
614 13331 : 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 2160 : a = std::pow(1.0 - std::pow(sl_bar, 1.0 / m), m);
619 : }
620 11171 : 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 1866 : return PorousFlowCubic::cubic(
625 1866 : 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 9305 : const Real sl_bar_del = (my_sldel - slr) / (1.0 - slr);
631 9305 : const Real s_gt_bar =
632 9305 : my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
633 9305 : a = (1 - s_gt_bar / (1.0 - sl_bar_del)) *
634 9305 : std::pow(1.0 - std::pow(sl_bar + s_gt_bar, 1.0 / m), m);
635 9305 : 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 21392 : return std::sqrt(sl_bar) * Utility::pow<2>(1.0 - a - b);
639 : }
640 :
641 : Real
642 30099 : 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 30099 : if (sl <= slr) // by the definition of slr, always return 0
655 : return 0.0;
656 23235 : if (sl == 1.0) // derivative is infinite at this point
657 : return std::numeric_limits<Real>::max();
658 22095 : const Real sl_bar = (sl - slr) / (1.0 - slr); // effective saturation
659 22095 : 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 22095 : if (sgrdel == 0.0 || sl <= sldel) // along the drying curve
667 : {
668 9222 : const Real c = std::pow(sl_bar, 1.0 / m);
669 9222 : const Real dc_dsbar = c / m / sl_bar;
670 9222 : a = std::pow(1.0 - c, m);
671 9222 : const Real da_dsbar = -m * a / (1.0 - c) * dc_dsbar;
672 9222 : a_prime = da_dsbar * sl_bar_prime;
673 9222 : }
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 12873 : const Real my_sldel = (sldel < slr) ? slr : sldel;
680 12873 : const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
681 12873 : 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 1712 : const Real c = std::pow(sl_bar, 1.0 / m);
686 1712 : const Real dc_dsbar = c / m / sl_bar;
687 1712 : a = std::pow(1.0 - c, m);
688 1712 : const Real da_dsbar = -m * a / (1.0 - c) * dc_dsbar;
689 1712 : a_prime = da_dsbar * sl_bar_prime;
690 : }
691 11161 : 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 1862 : return PorousFlowCubic::dcubic(
696 1862 : 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 9299 : const Real sl_bar_del = (my_sldel - slr) / (1.0 - slr);
702 9299 : const Real s_gt_bar =
703 9299 : my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
704 9299 : const Real s_gt_bar_prime = my_sgrdel / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
705 9299 : const Real c = std::pow(sl_bar + s_gt_bar, 1.0 / m);
706 9299 : const Real c_prime = c / m / (sl_bar + s_gt_bar) * (sl_bar_prime + s_gt_bar_prime);
707 9299 : a = (1 - s_gt_bar / (1.0 - sl_bar_del)) * std::pow(1.0 - c, m);
708 9299 : a_prime =
709 9299 : -s_gt_bar_prime / (1.0 - sl_bar_del) * std::pow(1.0 - c, m) - m * a / (1.0 - c) * c_prime;
710 9299 : b = s_gt_bar / (1.0 - sl_bar_del) * std::pow(1.0 - std::pow(sl_bar_del, 1.0 / m), m);
711 9299 : b_prime = s_gt_bar_prime * b / s_gt_bar;
712 : }
713 : }
714 20233 : const Real kr = std::sqrt(sl_bar) * Utility::pow<2>(1.0 - a - b);
715 20233 : return 0.5 * kr / sl_bar * sl_bar_prime -
716 20233 : std::sqrt(sl_bar) * 2.0 * (1.0 - a - b) * (a_prime + b_prime);
717 : }
718 :
719 : Real
720 16887 : 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 16887 : if (sl < slr)
731 : {
732 : // in the extended region, so immediately return with the relevant value
733 4102 : if (k_rg_max == 1.0)
734 : return 1.0;
735 2790 : return PorousFlowCubic::cubic(sl, 0.0, 1.0, 0.0, slr, k_rg_max, y0p);
736 : }
737 12785 : if (sl > 1.0 - sgrdel) // saturation is above 1.0 - residual gas saturation
738 : return 0.0;
739 12224 : 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 12224 : 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 6158 : const Real my_sldel = (sldel < slr) ? slr : sldel;
748 6158 : const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
749 6158 : s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
750 : }
751 : Real kr = 0.0;
752 12224 : 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 11839 : const Real a = std::pow(1.0 - (sl_bar + s_gt_bar), gamma);
756 11839 : const Real c = std::pow(sl_bar + s_gt_bar, 1.0 / m);
757 11839 : const Real b = std::pow(1.0 - c, 2.0 * m);
758 11839 : kr = k_rg_max * a * b;
759 : }
760 : return kr;
761 : }
762 :
763 : Real
764 16893 : 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 16893 : if (sl < slr)
775 : {
776 : // in the extended region, so immediately return with the relevant value
777 4098 : if (k_rg_max == 1.0)
778 : return 0.0;
779 2786 : return PorousFlowCubic::dcubic(sl, 0.0, 1.0, 0.0, slr, k_rg_max, y0p);
780 : }
781 12795 : if (sl > 1.0 - sgrdel) // saturation is above 1.0 - residual gas saturation
782 : return 0.0;
783 12234 : const Real sl_bar = (sl - slr) / (1.0 - slr);
784 12234 : 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 12234 : 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 6149 : const Real my_sldel = (sldel < slr) ? slr : sldel;
794 6149 : const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
795 6149 : s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
796 6149 : s_gt_bar_prime = my_sgrdel / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
797 : }
798 : Real kr_prime = 0.0;
799 12234 : 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 11849 : const Real a = std::pow(1.0 - (sl_bar + s_gt_bar), gamma);
803 11849 : const Real a_prime = -gamma * a / (1.0 - (sl_bar + s_gt_bar)) * (sl_bar_prime + s_gt_bar_prime);
804 11849 : const Real c = std::pow(sl_bar + s_gt_bar, 1.0 / m);
805 : const Real c_prime =
806 11849 : (c == 0 ? 0.0 : c / m / (sl_bar + s_gt_bar) * (sl_bar_prime + s_gt_bar_prime));
807 11849 : const Real b = std::pow(1.0 - c, 2.0 * m);
808 11849 : const Real b_prime = -2.0 * m * b / (1.0 - c) * c_prime;
809 11849 : kr_prime = k_rg_max * (a * b_prime + a_prime * b);
810 : }
811 : return kr_prime;
812 : }
813 : }
|
