Hysteresis in PorousFlow
Hysteretic capillary pressure and relative permeability functions are implemented in PorousFlow. Hysteresis means that capillary pressure and relative permeability depend on the saturation history as well as the current state.
The implementation of hysteresis in PorousFlow closely follows the TOUGH implementation as described in Doughty (2007) and Doughty (2008). Some other pertinent material may be found in Jaynes (1994), Niemi and Bodvarsson (1988) and Finsterle et al. (1998). In this implementation, capillary pressure at any point in the finite-element mesh is described by one of four "curves" (or functions), and relative permeability at any point is described by one of two "curves". The capillary-pressure curves are illustrated in Figure 1, which shall now be qualitatively described.
Figure 1: Illustration of a hysteretic capillary pressure function.
Assume that the point in question starts at full liquid saturation, . Its history then proceeds as follows.
The liquid saturation is gradually reduced. This is called "drying" or "draining". The capillary pressure follows the order zero curve or zeroth-order drying curve, which is also termed the "primary drying" curve.
At some time, the liquid saturation begins to increase. At this time, the liquid saturation in the figure ("TP" stands for "turning point", where ). Increasing liquid saturation is called "wetting" or "imbibing". The capillary pressure follows the order one or first-order wetting curve. This curve is an mixture of the primary drying curve and a modified version of the "primary wetting" curve. The mathematical description of the mixture is described below.
At some later time, the liquid saturation begins to decrease again. At this time, the liquid saturation in the figure. The capillary pressure follows the order two or second-order drying curve.
At some later time, the liquid saturation beings to increase again. The capillary pressure follows the order three or third-order curve.
At later times, the liquid saturation can decrease and increase, and provided the saturation obeys , the capillary pressure follows the third-order curve.
The saturation history just described means that the order increases. However, the order may also decrease, as illustrated in Figure 2.
The second-order wetting curve becomes the zeroth-order drying curve if .
The third-order curve becomes the second-order drying curve if .
The third-order curve becomes the first-order wetting curve if .
Figure 2: Hysteretic capillary pressure curves. This history is: (a) drying to ; (b) wetting to ; (c) drying to , where the second-order drying curve is followed to and then the zeroth-order drying curve is followed; (d) wetting to ; (e) drying to , where the second-order drying curve is followed to and then the zeroth-order drying curve is followed; (f) wetting to ; (g) drying to where the second-order drying curve is followed to and then the zeroth-order drying curve is followed.
The order is computed by the PorousFlowHysteresisOrder Material, as discussed in detail below.
Capillary curves
All capillary curves are based on the van Genuchten curve (1) where In these formulae:
is the capillary pressure, with SI units Pa. It obeys .
is the liquid saturation, which is dimensionless.
is a van Genuchten parameter, with SI units Pa. Two values may be entered by the user — one describing the primary drying curve and the other for the primary wetting curve — and both obey .
is a van Genuchten parameter, which is dimensionless. Two values may be entered by the user — one describing the primary drying curve and the other for the primary wetting curve – and both obey .
is the minimum saturation for which the van Genuchten expression is valid. It is a dimensionless user-input and obeys . Because , behaves like as . It is sometimes numerically advantageous and/or physically necessary to define extensions of for , which is discussed below. Two extensions are shown in Figure 3
is the residual gas saturation, and is the maximum saturation for which the van Genuchten expression is valid. It is dimensionless but it is not a user input (see more details below). It obeys and . As , . It is sometimes numerically advantageous and/or physically necessary to define extensions of for as described below, as shown in Figure 3.
The primary drying curve uses The primary wetting curve uses where is a user input (as noted below must be strictly positive for the cubic-spline modification of the water wetting relative permeability to be well-defined). The other curves use other values of as described below.
Figure 3: Various extensions of the hysteretic capillary pressure curve available in PorousFlow. The primary wetting curve is shown, which has the property that as (recall that the primary wetting curve uses the user-inputted for ).
Lower extensions
Two types of lower (small ) extensions are available (Doughty, 2008). They are based on the user inputting the value which is the value at which extension commences (that is, all will use the extension, not the original van Genuchten expression Eq. (1)). Figure 3 shows examples for . Both extensions are designed so that the resulting curve is continuous and its derivative is continuous.
Quadratic: is a quadratic in that satisfies at .
Exponential: is an exponential in .
Lower extensions could be necessary because the model could be heated to evaporate the water, or chemical reactions could occur to reduce below . They could also be necessary because the region might be explored as the MOOSE solver converges towards a solution.
Upper extensions
One type of upper extension (large ) is available (Doughty, 2008). It is based on the user inputting a value When then the upper extension is used: where the coefficient of proportionality and are chosen to ensure the result is C1-continuous. Figure 3 shows an example for . This type of extension has the property that . This type of extension is not necessary for the primary drying curve.
Upper extensions could be necessary to help the MOOSE solver converge towards a solution.
Zeroth-order drying curve
This is defined by Eq. (1) (along with an optional lower extension) with .
First-order wetting curve
The idea here is to use a wetting curve (with nonzero and the wetting and parameters) but with a value of liquid saturation that ensures continuity.
The first-order wetting curve uses a non-zero value for the quantity in Eq. (1). This is (Doughty, 2007) (2) This is the so-called "Land" expression. It depends on , which obeys and is the "liquid residual saturation" defined to be the saturation at which the liquid relative permeability tends to zero. The following may be noticed:
if then (assuming no upper extension). So is the saturation point at which the curves on Figure 1 tend to zero.
if is close to 0, then is close to 1. This means that the first-order wetting curve will be quite close to the zeroth-order drying curve.
if is close to then is close to . This means that the first-order wetting curve will be quite close to the primary wetting curve.
These points are illustrated in Figure 4.
Figure 4: Various first-order curves resulting from different turning points. Notice the different values of which depend on the turning-point saturation.
However, if Eq. (1) is used with , the result will be discontinuous at the turning point. This may be seen in Figure 5: the result would jump from the red drying curve to the blue wetting curve.
Figure 5: Construction of the first-order curve.
To ensure continuity, define (see Figure 5):
, which is the value of at , using the primary drying curve.
, which is the value of at the capillary pressure defined using Eq. (1) using the wetting and parameters and
Armed with these quantities, the first-order wetting curve is defined by Eq. (1) using the following:
the and parameters for the primary wetting curve
defined by the Land expression Eq. (2)
instead of , to ensure continuity of the result. Here
(3) This is designed so that
when then and
when then .
The result is the green curve in Figure 5 which smoothly transitions from the red drying curve and the blue wetting curve.
Second-order drying
This is drawn from Eqn(3) of Niemi and Bodvarsson (1988) and Eqn(5) of Finsterle et al. (1998) but may be slightly different because their explanations are a little opaque. The idea is similar to the first-order case: use the drying curve, but with a different saturation, , to ensure continuity of the result.
Define
to be the saturation on the primary drying curve corresponding to at
to smoothly interpolate between when , and when , ie
(4)
The second-order drying curve uses:
the and parameters of the primary drying curve
instead of , as defined by Eq. (4)
Third-order curves
The approach here is to take an appropriate average of the first-order wetting curve and the second-order drying curve, which is probably the same as Niemi and Bodvarsson (1988), Finsterle et al. (1998) and Doughty (2007), although details are sparse in those references.
Define
is the first-order wetting capillary pressure using defined by Eq. (3)
is the second-order drying capillary pressure using defined by Eq. (4)
Then, the third-order capillary pressure is defined by
This produces continuous capillary-pressure functions as shown in Figure 1.
Relative permeabilities
Hysteresis is defined for one-phase and two-phase systems only. The phase is assumed to be the liquid phase in one-phase models. The liquid and gas relative permeability functions are both hysteretic. Only the drying and first-order wetting curves are defined. This means that if the system dries (following the drying curve) and then wets (following the first-order wetting curve) and then subsequently dries, the system will move along the first-order wetting curve until the turning point is reached, when it will move along the drying curve.
The starting point for the relative permeability functions for liquid () and gas () are of the van Genuchten form Doughty (2007) (5) In these expressions Here
is the residual liquid saturation, which marks the point at which the liquid relative permeability becomes zero. It is a dimensionless user-input, which satisfies .
is the maximum gas relative permeability, which occurs on the drying curve for . It is a dimensionless user-input, which satisfies . Most frequently is chosen.
is a dimensionless index, inputted by the user, which satisfies . Most frequently is chosen.
is a dimensionless index, inputted by the user, which satisfies .
Extending and modifying the functions
An immediate problem with the functions given in Eq. (5) is that they are not defined for (meaning that ). The wetting curve for water also has an infinite derivative. These behaviors are obviously undesirable in a numerical implementation, and are probably unphysical, so the functions need to be extended and modified.
For reference, the zeroth-order (drying) are given by Eq. (5) but because in this case, , and the functions reduce to (6)
Assuming the relative permeability values are constant outside the well-defined domain results in curves of the form shown in Figure 6.
Figure 6: The basic, unextended relative permeability curves.
Doughty (2008) recommends extending and modifying the curves in the following ways.
There is no hysteresis in the region : the drying curves equal the wetting curves in this region, and if then is used in the wetting curves.
The gas relative permeability is extended to the region so that . Two types of extension are possible in MOOSE. Both are cubic functions satisfying , and :
A "linear-like" extension, where the cubic's derivative at equals . That is, the derivative is equal to the average slope in the extended region. This means the final result has a discontinuous derivative at , but the result often "looks better" to the eye (see figures below).
A "cubic" extension, where the cubic's derivative at equals the primary drying-curve's derivative at that point. This means the drying curve is C1 continuous (the wetting curve is not).
The water wetting curve is modified around the point , which is the point of infinite derivative. Firstly, it is assumed that . If this is not the case then the following modification will not work, but may be unnecessary anyway. Two points are defined
, where is a (dimensionless) user-input satisfying , and optimally should be close to 1, for example .
. Note that Doughty (2008) defines , but this can be greater than 1 and can result in a poor modification.
Then:
for the liquid wetting curve, given in Eq. (5), is used
for the liquid drying curve, given in Eq. (6), is used
otherwise, a cubic spline interpolates between these two, with parameters chosen so that the final result is continuous and has a continuous derivative.
The result is illustrated in Figure 7, using a linear-like extension for the gas relative permeability curve, and .
Figure 7: Extended and modified relative permeability curves, using a linear-like extension for the gas relative permeability curve, and being 0.9.
Examples
Figure 8: Example hysteretic relative permeability functions. The system initialises at full saturation, dries to saturation 0.5 (greater than ) and then wets back to saturation 1.
Figure 9: Example hysteretic relative permeability functions. The system initialises at full saturation, dries to saturation 0.1 (less than ) and then wets back to saturation 1. In the extended region (), the drying curve equals the wetting curve.
Exploring hysteresis using the python script
A python script that produces plots of hysteretic capillary pressure and relative permeability has been included in the MOOSE repository. This allows users to quickly explore the consequences of choosing different extension strategies to the capillary and relative-permeability curves, as well as the impact of hysteresis order on the shape of the curves. It was used to produce all the figures displayed on this page. It is
# Draws Pc curves and relperm curves
import os, sys
import math
import matplotlib.pyplot as plt
plt.rcParams.update({'figure.max_open_warning': 0})
def pc_fcn(sl, slmin, sgrdel, alpha, n, sm, pcm, pcmprime, lower_extension_type, sm_upper, pcm_upper, pcmprime_upper, upper_extension_type):
# Eqn(1) of Doughty2007
# with extensions defined in Doughty2008, page5 and Fig1
# lower extension strategy:
# User inputs Pcmax, which is a value of Pc on the main drainage curve.
# From this, sm is worked out. For sl < sm, an extension is used for the main drainage curve
# Clearly, sm > slmin because to get sm=slmin, absPcmax = infinity would have to be used
# Similarly for the main wetting curve: Pcmax gives value of sm
# Upper extension strategy:
# User inputs r < 1.0 (usually r ~ 1.0), and code calculates sm_upper = r * (1 - sgrdel).
# For sl > sm_upper, the upper extension is used. This is a power law that is zero at sl = 1, and matches the main curve's value (pcm_upper) and derivative (pcmprime_upper) at sm_upper.
if (sl < sm): # important for initialising pcm and pcmprime that this is <, not <=
if lower_extension_type == "quadratic":
return pcm + pcmprime * (sl * sl - sm * sm) / 2.0 / sm
elif lower_extension_type == "exponential":
return pcm * math.exp(pcmprime * (sl - sm) / pcm)
else:
return pcm
if (sl > sm_upper): # important for initialising pcm_upper and pcmprime_upper that this is >, not >=
if upper_extension_type == "power":
expon = -pcmprime_upper / pcm_upper * (1.0 - sm_upper)
if sl >= 1.0:
return 0.0
return pcm_upper * math.pow((1.0 - sl) / (1.0 - sm_upper), expon)
seff = (sl - slmin) / (1 - sgrdel - slmin)
if (seff >= 1.0):
return 0.0 # no valid upper-extension_type defined
elif (seff <= 0.0):
return 100.0 # no valid lower_extension_type defined
return -(1.0 / alpha) * pow(pow(seff, n / (1.0 - n)) - 1.0, 1.0 / n)
def pc_fcn_prime(sl, slmin, sgrdel, alpha, n, sm, pcm, pcmprime, lower_extension_type, sm_upper, pcm_upper, pcmprime_upper, upper_extension_type):
# derivative of pc_fcn wrt sl
if (sl < sm): # important for initialising pcm and pcmprime that this is <, not <=
if lower_extension_type == "quadratic":
return pcmprime * sl / sm
elif lower_extension_type == "exponential":
return pcmprime * math.exp(pcmprime * (sl - sm) / pcm)
if (sl > sm_upper): # important for initialising pcm_upper and pcmprime_upper that this is >, not >=
if upper_extension_type == "power":
return pcmprime_upper * math.pow((1.0 - sl) / (1.0 - sm_upper), -pcmprime_upper / pcm_upper * (1.0 - sm_upper) - 1.0)
seff = (sl - slmin) / (1.0 - sgrdel - slmin)
if (seff >= 1.0):
return 0.0 # no valid upper-extension_type defined
elif (seff <= 0.0):
return 0.0 # no valid lower_extension_type defined
dseff = 1.0 / (1.0 - sgrdel - slmin)
dpc_dseff = -(1.0 / alpha) * pow(pow(seff, n / (1.0 - n)) - 1.0, 1.0 / n - 1.0) * pow(seff, n / (1.0 - n) - 1.0) / (1.0 - n)
return dpc_dseff * dseff
def sat_fcn(pc, slmin, sgrdel, alpha, n, sm, pcm, pcmprime, lower_extension_type, sm_upper, pcm_upper, pcmprime_upper, upper_extension_type):
# inverse of pc_fcn
if (pc >= 0):
return 1.0
if (pc < pcm): # Important for initialisation that the condition is < not <=. Remember that pc < 0
if lower_extension_type == "quadratic":
s2 = sm * sm + (pc - pcm) * 2.0 * sm / pcmprime
if s2 <= 0.0: # this occurs when we're tyring to find the saturation on the wetting curve defined by self.gr_Del at pc = self.Pcd_Del, if this pc is actually impossible to achieve on this wetting curve
return 0.0
return math.sqrt(s2)
elif lower_extension_type == "exponential":
s = sm + math.log(pc / pcm) * pcm / pcmprime
if s < 0.0: # this occurs when we're tyring to find the saturation on the wetting curve defined by self.gr_Del at pc = self.Pcd_Del, if this pc is actually impossible to achieve on this wetting curve
return 0.0
return s
else:
return sm
if (pc > pcm_upper): # Important for initialisation that the condition is > not >=. Remember that pc < 0
if upper_extension_type == "power":
expon = -pcmprime_upper / pcm_upper * (1.0 - sm_upper)
return 1.0 - math.power(pc / pcm_upper, 1.0 / expon) * (1.0 - sm_upper)
seffpow = 1.0 + pow(-alpha * pc, n)
seff = pow(seffpow, (1.0 - n) / n)
return (1.0 - sgrdel - slmin) * seff + slmin
def cubic(x, x0, y0, y0p, x1, y1, y1p):
''' Cubic f(x) that satisfies
f(x0) = y0
f'(x0) = y0p
f(x1) = y1
f'(x1) = y1p
'''
d = x1 - x0
d2 = math.pow(d, 2.0)
mean = 0.5 * (x1 + x0)
sq3 = 0.5 * math.sqrt(3.0) * d
term1 = y0p * (x - x0) * math.pow(x - x1, 2) / d2 # term1(x0) = term1(x1) = term1'(x1) = 0, term1'(x0) = y0p
term2 = y1p * (x - x1) * math.pow(x - x0, 2) / d2 # term2(x0) = term2(x1) = term2'(x0) = 0, term2'(x1) = y1p
term3 = (x - mean - sq3) * (x - mean) * (x - mean + sq3)
# term3' = (x - mean) * (x - mean + sq3) + (x - mean - sq3) * (x - mean + sq3) + (x - mean - sq3) * (x - mean)
# = 3 (x - mean)^2 - sq3^2
# note term3' = 0 when x = mean +/- sq3/sqrt(3) = 0.5 * (x1 + x0) +/- 0.5 * (x1 - x0) = {x1, x0}
term3_x0 = (x0 - mean - sq3) * (x0 - mean) * (x0 - mean + sq3)
term3_x1 = (x1 - mean - sq3) * (x1 - mean) * (x1 - mean + sq3)
return (y0 * (term3 - term3_x1) + y1 * (term3_x0 - term3)) / (term3_x0 - term3_x1) + term1 + term2
def krel_liquid(sl, slr, sgrdel, sgrmax, sldel, m, upper_liquid_param, y0, y0p, y1, y1p):
if (sl < slr):
return 0.0
sl_bar = (sl - slr) / (1.0 - slr)
if (sgrdel == 0.0 or sl <= sldel): # along the drying curve
s_gt_bar = 0.0
a = math.pow(1.0 - math.pow(sl_bar, 1.0 / m), m)
b = 0.0
else:
if (sldel < slr): # turning point occured below slr, but "there is no hysteresis along the extension" according to p6 of Doughty2008. The above checked that sl >= slr. I assume the wetting relperm curve is the same as if the turning point = slr:
my_sldel = slr
my_sgrdel = sgrmax
else:
my_sldel = sldel
my_sgrdel = sgrdel
# This is Doughty2008: if sl >= (2.0 - upper_liquid_param) * (1.0 - my_sgrdel): # imporant for initialisation of cubic that this is >= and not >
if sl >= 1.0 - 0.5 * my_sgrdel: # imporant for initialisation of cubic that this is >= and not >
# follow the drying curve
s_gt_bar = 0.0
a = math.pow(1.0 - math.pow(sl_bar, 1.0 / m), m)
b = 0.0
elif sl > upper_liquid_param * (1.0 - my_sgrdel): # important for initialisation of cubic that this is > and not >=
# this is Doughty2008: return cubic(sl, upper_liquid_param * (1.0 - my_sgrdel), y0, y0p, (2.0 - upper_liquid_param) * (1.0 - my_sgrdel), y1, y1p)
return cubic(sl, upper_liquid_param * (1.0 - my_sgrdel), y0, y0p, 1.0 - 0.5 * my_sgrdel, y1, y1p)
else: # standard case
sl_bar_del = (my_sldel - slr) / (1.0 - slr)
s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel)
a = (1 - s_gt_bar / (1.0 - sl_bar_del)) * math.pow(1.0 - math.pow(sl_bar + s_gt_bar, 1.0 / m), m)
b = s_gt_bar / (1.0 - sl_bar_del) * math.pow(1.0 - math.pow(sl_bar_del, 1.0 / m), m)
return math.sqrt(sl_bar) * math.pow(1.0 - a - b, 2)
def krel_liquid_prime(sl, slr, sgrdel, sgrmax, sldel, m, upper_liquid_param):
if (sl < slr):
return 0.0
sl_bar = (sl - slr) / (1.0 - slr)
sl_bar_prime = 1.0 / (1.0 - slr)
if sgrdel == 0.0: # along the drying curve
s_gt_bar = 0.0
c = math.pow(sl_bar, 1.0 / m)
dc_dsbar = c / m / sl_bar
a = math.pow(1.0 - c, m)
da_dsbar = - m * a / (1.0 - c) * dc_dsbar
aprime = da_dsbar * sl_bar_prime
b = 0.0
bprime = 0.0
else:
if (sldel < slr): # turning point occured below slr, but "there is no hysteresis along the extension" according to p6 of Doughty2008. The above checked that sl >= slr. I assume the wetting relperm curve is the same as if the turning point = slr:
my_sldel = slr
my_sgrdel = sgrmax
else:
my_sldel = sldel
my_sgrdel = sgrdel
# this is Doughty2008: if sl >= (2.0 - upper_liquid_param) * (1.0 - my_sgrdel): # imporant for initialisation of cubic that this is >= and not >
if sl >= 1.0 - 0.5 * my_sgrdel: # imporant for initialisation of cubic that this is >= and not >
# follow the drying curve
s_gt_bar = 0.0
c = math.pow(sl_bar, 1.0 / m)
dc_dsbar = c / m / sl_bar
a = math.pow(1.0 - c, m)
da_dsbar = -m * a / (1.0 - c) * dc_dsbar
aprime = da_dsbar * sl_bar_prime
b = 0.0
bprime = 0.0
elif sl > upper_liquid_param * (1.0 - my_sgrdel): # important that this is > and not >=
return 100.0 # NOTE: this is incorrect, but it is never used
else: # standard case
sl_bar_del = (my_sldel - slr) / (1.0 - slr)
s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel)
s_gt_bar_prime = my_sgrdel / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel)
c = math.pow(sl_bar + s_gt_bar, 1.0 / m)
cprime = c / m / (sl_bar + s_gt_bar) * (sl_bar_prime + s_gt_bar_prime)
a = (1 - s_gt_bar / (1.0 - sl_bar_del)) * math.pow(1.0 - c, m)
aprime = -s_gt_bar_prime / (1.0 - sl_bar_del) * math.pow(1.0 - c, m) - m * a / (1.0 - c) * cprime
b = s_gt_bar / (1.0 - sl_bar_del) * math.pow(1.0 - math.pow(sl_bar_del, 1.0 / m), m)
bprime = s_gt_bar_prime * b / s_gt_bar
kr = math.sqrt(sl_bar) * math.pow(1.0 - a - b, 2)
return 0.5 * kr / sl_bar * sl_bar_prime - math.sqrt(sl_bar) * 2.0 * (1.0 - a - b) * (aprime + bprime)
def krel_gas(sl, slr, sgrdel, sgrmax, sldel, m, gamma, k_rg_max, y0p):
if (sl < slr):
if (k_rg_max == 1.0):
return 1.0
return cubic(sl, 0.0, 1.0, 0.0, slr, k_rg_max, y0p)
sl_bar = (sl - slr) / (1.0 - slr)
if (sgrdel == 0.0 or sl < sldel): # on the drying curve
s_gt_bar = 0.0
else:
if (sldel < slr): # turning point occured below slr, but "there is no hysteresis along the extension" according to p6 of Doughty2008. The above checked that sl >= slr. I assume the wetting relperm curve is the same as if the turning point = slr:
my_sldel = slr
my_sgrdel = sgrmax
else:
my_sldel = sldel
my_sgrdel = sgrdel
s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel)
if (sl_bar + s_gt_bar > 1.0):
return 0.0 # always zero, irrespective of any hysteresis
a = math.pow(1.0 - (sl_bar + s_gt_bar), gamma)
c = math.pow(sl_bar + s_gt_bar, 1.0 / m)
b = math.pow(1.0 - c, 2.0 * m)
return k_rg_max * a * b
def krel_gas_prime(sl, slr, sgrdel, sgrmax, sldel, m, gamma, k_rg_max):
if (sl < slr):
if (k_rg_max == 1.0):
return 0.0
return 0.0 # NOTE: this is incorrect but it is never used
if (sl > 1.0 - sgrdel):
return 0.0 # always zero, irrespective of any hysteresis
sl_bar = (sl - slr) / (1.0 - slr)
sl_bar_prime = 1.0 / (1.0 - slr)
if sgrdel == 0.0: # on the drying curve
s_gt_bar = 0.0
s_gt_bar_prime = 0.0
else:
if (sldel < slr): # turning point occured below slr, but "there is no hysteresis along the extension" according to p6 of Doughty2008. The above checked that sl >= slr. I assume the wetting relperm curve is the same as if the turning point = slr:
my_sldel = slr
my_sgrdel = sgrmax
else:
my_sldel = sldel
my_sgrdel = sgrdel
s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel)
s_gt_bar_prime = my_sgrdel / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel)
a = math.pow(1.0 - (sl_bar + s_gt_bar), gamma)
aprime = -gamma * a / (1.0 - (sl_bar + s_gt_bar)) * (sl_bar_prime + s_gt_bar_prime)
c = math.pow(sl_bar + s_gt_bar, 1.0 / m)
cprime = 1.0 / m * math.pow(sl_bar + s_gt_bar, 1.0 / m - 1.0) * (sl_bar_prime + s_gt_bar_prime)
b = math.pow(1.0 - c, 2.0 * m)
bprime = -2.0 * m * b / (1.0 - c) * cprime
return k_rg_max * (a * bprime + aprime * b)
class Hys:
def __init__(self, S_lr, m, gamma, k_rg_max, upper_liquid_param, S_lmin, S_grmax, alpha_d, alpha_w, n_d, n_w, absPcmax, lower_extension_type, upper_ratio, upper_extension_type, relperm_gas_extension):
self.S_lr = S_lr # liquid residual saturation, where liquid relperm -> 0 and gas relperm = k_rg_max
self.m = m # exponent in liquid relperm
self.gamma = gamma # exponent in gas relperm
self.k_rg_max = k_rg_max # maximum value of gas relperm (value attained S_liquid = S_lr)
self.upper_liquid_param = upper_liquid_param # cubic spline extension for liquid relperm. must be positive and <= 1. should be close to 1
self.S_lmin = S_lmin # liquid saturation where van-Genuchten -> infinity
self.S_grmax = S_grmax # gas saturation where wetting cap -> 0
self.alpha_d = alpha_d # drying alpha
self.alpha_w = alpha_w # wetting alpha
self.n_d = n_d # drying n (= 1/(1 - m) for m being MOOSE param)
self.n_w = n_w # wetting n (= 1/(1 - m) for m being MOOSE param)
self.reinit_turning_points()
self.absPcmax = absPcmax # defines point extension of capillary pressure at low saturation
self.lower_extension_type = lower_extension_type # quadratic or exponential type of extension
self.S_m_drying = sat_fcn(-self.absPcmax, self.S_lmin, 0.0, self.alpha_d, self.n_d, 0.0, -self.absPcmax, 0.0, self.lower_extension_type, 1.0, 0.0, 0.0, "none") # saturation at point of lower extension on primary drying curve
self.dPc_dSm_drying = pc_fcn_prime(self.S_m_drying, self.S_lmin, 0.0, self.alpha_d, self.n_d, 0.0, -self.absPcmax, 0.0, self.lower_extension_type, 1.0, 0.0, 0.0, "none") # deriviative at point of lower extension on primary drying curve
self.S_m_wetting = sat_fcn(-self.absPcmax, self.S_lmin, self.S_grmax, self.alpha_w, self.n_w, 0.0, -self.absPcmax, 0.0, self.lower_extension_type, 1.0, 0.0, 0.0, "none") # saturation at point of lower extension on primary wetting curve, with no upper extension at this point
self.dPc_dSm_wetting = pc_fcn_prime(self.S_m_wetting, self.S_lmin, self.S_grmax, self.alpha_w, self.n_w, 0.0, -self.absPcmax, 0.0, self.lower_extension_type, 1.0, 0.0, 0.0, "none") # deriviative at point of lower extension on primary wetting curve, with no upper extension at this point
self.upper_ratio = upper_ratio # must be positive and < 1.0
self.upper_extension_type = upper_extension_type # if "power" then power-law extension of wetting curve at upper end, otherwise no extension
self.S_m_upper = self.upper_ratio * (1.0 - self.S_grmax) # saturation at which upper extension starts on the primary wetting curve
self.Pcm_upper = pc_fcn(self.S_m_upper, self.S_lmin, self.S_grmax, self.alpha_w, self.n_w, self.S_m_wetting, -self.absPcmax, self.dPc_dSm_wetting, self.lower_extension_type, self.S_m_upper, 0.0, 0.0, self.upper_extension_type) # pc on primary wetting curve at point of upper extension
self.Pcmprime_upper = pc_fcn_prime(self.S_m_upper, self.S_lmin, self.S_grmax, self.alpha_w, self.n_w, self.S_m_wetting, -self.absPcmax, self.dPc_dSm_wetting, self.lower_extension_type, self.S_m_upper, self.Pcm_upper, 0.0, self.upper_extension_type) # dpc/ds on primary wetting curve at point of upper extension
if (relperm_gas_extension == "linear_like"):
self.krel_gas_prime = (self.k_rg_max - 1.0) / self.S_lr # this means the slope for the gas relperm at S_lr equals the slope of a straightline between (0, 1) and (S_lr, k_rg_max)
else:
self.krel_gas_prime = krel_gas_prime(self.S_lr, self.S_lr, 0.0, self.S_grmax, 1.0, self.m, self.gamma, self.k_rg_max) # This ensures the drying gas relperm is C1, but the wetting as relperm may have discontinuous deriv at self.S_lr if the turning point was less than self.S_lr
def primary_drying(self, saturations):
''' Returns the absolute value of the primary drying capillary pressure values for the saturations. There is no upper extension on this curve '''
return [abs(pc_fcn(s, self.S_lmin, 0.0, self.alpha_d, self.n_d, self.S_m_drying, -self.absPcmax, self.dPc_dSm_drying, self.lower_extension_type, 1.0, 0.0, 0.0, "none")) for s in saturations]
def primary_drying_sats(self, pcs):
''' Returns the absolute value of the saturations resulting from the pc values given, using the primary drying curve'''
return [sat_fcn(pc, self.S_lmin, 0.0, self.alpha_d, self.n_d, self.S_m_drying, -self.absPcmax, self.dPc_dSm_drying, self.lower_extension_type, 1.0, 0.0, 0.0, "none") for pc in pcs]
def primary_wetting(self, saturations):
''' Returns the absolute value of the primary wetting capillary pressure values for the saturations '''
return [abs(pc_fcn(s, self.S_lmin, self.S_grmax, self.alpha_w, self.n_w, self.S_m_wetting, -self.absPcmax, self.dPc_dSm_wetting, self.lower_extension_type, self.S_m_upper, self.Pcm_upper, self.Pcmprime_upper, self.upper_extension_type)) for s in saturations]
def relperm_liquid_drying(self, saturations):
''' Returns the liquid relperm during primary drying for the saturations. '''
return [krel_liquid(s, self.S_lr, 0.0, self.S_grmax, 1.0, self.m, self.upper_liquid_param, 0.0, 0.0, 0.0, 0.0) for s in saturations]
def relperm_gas_drying(self, saturations):
''' Returns the gas relperm during primary drying for the saturations. '''
return [krel_gas(s, self.S_lr, 0.0, self.S_grmax, 1.0, self.m, self.gamma, self.k_rg_max, self.krel_gas_prime) for s in saturations]
def s_grdel_land(self, slDel):
''' Computes S_gr^Delta as per Eqn(2) of Doughty2007. I *think* this is only for wetting (imbibition) curves, because it is zero for drying (drainage) curves '''
a = 1.0 / self.S_grmax - 1.0 / (1.0 - self.S_lr)
return (1.0 - slDel) / (1.0 + a * (1.0 - slDel))
def reinit_turning_points(self):
''' Reinitialises the history of turning points '''
self.S_lDel = [] # list of turning points
self.Pcd_Del = [] # Pc(self.S_lDel) calculated using the primary drying curve
self.gr_Del = [] # S_gr^Del(self.S_lDel) = Land S_gr^Del calculated at each turning point
self.lc_Del = [] # 1 - self.gr_Del
self.S_m_wetting_Del = [] # saturation at point of extension on the wetting curve defined by self.gr_Del
self.dPc_dS_m_wetting_Del = [] # derivative at point of extension on the wetting curve defined by self.gr_Del
self.sl_wetting_Del = [] # saturation on the wetting curve defined by self.gr_Del, at pc = self.Pcd_Del
self.Pc_Del = [] # Pc calulcated by the scanning-curve procedure at self.S_lDel
self.Sd_Del = [] # saturation on the drying curve using Pc_Del
self.S_m_upper_Del = [] # saturation on the wetting curve defined by self.gr_Del at which upper extension starts
self.Pcm_upper_Del = [] # Pc on the wetting curve defined by self.gr_Del, at saturation = self.S_m_upper_Del
self.Pcmprime_upper_Del = [] # dPc/dS on the wetting curve defined by self.gr_Del, at saturation = self.S_m_upper_Del
self.kl_begin_Del = [] # liquid relperm value at self.upper_liquid_param * self.lc_Del, used in cubic-spline transition from wetting to drying curve
self.klp_begin_Del = [] # derivative of liquid relperm value at self.upper_liquid_param * self.lc_Del, used in cubic-spline transition from wetting to drying curve
self.kl_end_Del = [] # liquid relperm value at (2 - self.upper_liquid_param) * self.lc_Del, used in cubic-spline transition from wetting to drying curve
self.klp_end_Del = [] # derivative of liquid relperm value at (2 - self.upper_liquid_param) * self.lc_Del, used in cubic-spline transition from wetting to drying curve
self.order = 0 # order of curve
self.prev_saturation = 1.0 # previous saturation encountered
self.latest_pc = 0.0 # latest value of Pc calculated
self.latest_kl = 1.0 # latest value of liquid relperm calculated
self.latest_kg = 0.0 # latest value of gas relperm calculated
def pc(self, saturations, reinit_tp):
''' Computes pc given the saturation history contained in the "saturations" list. if reinit_tp == true then the turning points and order are reinitialised to their 'primary drying' values '''
if (reinit_tp):
self.reinit_turning_points()
if (len(saturations) == 0):
return []
self.set_pc_val(saturations[0])
pc_vals = [self.latest_pc]
self.prev_saturation = saturations[0]
self.define_order(saturations[0])
for s in saturations[1:]:
self.define_order(s)
self.set_pc_val(s)
pc_vals.append(self.latest_pc)
return pc_vals
def set_pc_val(self, s):
'''
Sets self.latest_pc, which is the capillary pressure for the saturation, based on self.order
'''
if self.order == 0:
self.latest_pc = pc_fcn(s, self.S_lmin, 0.0, self.alpha_d, self.n_d, self.S_m_drying, -self.absPcmax, self.dPc_dSm_drying, self.lower_extension_type, 1.0, 0.0, 0.0, "none") # no upper extension on primary drying curve
elif self.order == 1:
# define the interpolation: s_to_use smoothly transitions from Sl_wetting_Del[0] (the value of Sl on the wetting curve defined by self.gr_Del[0]) when s = S_lDel[0] (the value of Sl on the drying curve when we transitioned to wetting), to lc_Del[0] (the value of Sl at Pc=0 on the wetting curve defined by self.gr_Del[0]) when s = lc_Del[0]
if self.upper_extension_type == "power":
s_to_use = self.sl_wetting_Del[0] + (1.0 - self.sl_wetting_Del[0]) * (s - self.S_lDel[0]) / (1.0 - self.S_lDel[0]) # in this case lc_Del[0] is effectively 1.0 since the wetting capillary curve is extended to Sl=1
else:
s_to_use = self.sl_wetting_Del[0] + (self.lc_Del[0] - self.sl_wetting_Del[0]) * (s - self.S_lDel[0]) / (self.lc_Del[0] - self.S_lDel[0])
self.latest_pc = pc_fcn(s_to_use, self.S_lmin, self.gr_Del[0], self.alpha_w, self.n_w, self.S_m_wetting_Del[0], -self.absPcmax, self.dPc_dS_m_wetting_Del[0], self.lower_extension_type, self.S_m_upper_Del[0], self.Pcm_upper_Del[0], self.Pcmprime_upper_Del[0], self.upper_extension_type)
elif self.order == 2:
# NOTE: here i think i don't use Eqn(3) from Niemi and Bodvarsson (which is Eqn(5) from Finsterle, Sonnenborg and Faybishenko). I couldn't get their idea to work, so i made up something i think is sensible
# The idea is to define capillary pressure using the primary drying curve, but use a different saturation than just "s".
# This is similar to the order=1 case, where the capillary pressure is defined using the wetting curve (with sgrdel) but the saturation is not just "s", but s_to_use.
# At the latest turning point:
sdprime = self.Sd_Del[1] # saturation on the primary drying curve corresponding to Pc calculated at the latest turning point
# At the previous turning point:
sd = self.S_lDel[0] # saturation at the previous turning point
# The idea is:
# if s = self.S_lDel[1] then s_to_use = sdprime
# if s = self.S_lDel[0] then s_to_use = sd
# An alternative idea would be to use a quadratic that obeys the above constraints and has slope = 1 at s = self.S_lDel[0] because that would give a continuous slope of the curve at that point. I don't know if that is advantagious or not
s_to_use = sd + (s - self.S_lDel[0]) * (sdprime - sd) / (self.S_lDel[1] - self.S_lDel[0])
self.latest_pc = pc_fcn(s_to_use, self.S_lmin, 0.0, self.alpha_d, self.n_d, self.S_m_drying, -self.absPcmax, self.dPc_dSm_drying, self.lower_extension_type, 1.0, 0.0, 0.0, "none") # no upper extension on drying curve
# The following commented-out stuff is Niemi et al, i think
# slc = 1.0 - sgrdel
# pcplus = pc_fcn(s, self.S_lmin, 0.0, self.alpha_d, self.n_d) # on the original drying curve
# slwplus = sat_fcn(pcplus, self.S_lmin, sgrdel, self.alpha_w, self.n_w) # on the wetting curve with S_gr = S_grdel but not the scanning version
# lhs = (s - self.S_lDel[-1]) * pow(slc - slwplus, 2.0) / (slc - s) # MOOSE: be careful with s=slc
# s_to_use = 0.5 * (slc + self.S_lDel[-1]) - 0.5 * math.sqrt(pow(slc - self.S_lDel[-1], 2.0) - 4 * lhs)
# return pc_fcn(s_to_use, self.S_lmin, sgrdel, self.alpha_w, self.n_w)
else:
# Niemi and Bodvarsson aruge that a straight-line fit between the last two reversal points is appropriate, which is the approach taken by Finsterle, Sonnenborg and Faybishenko too, as well as Doughty2007 (except that in Dought2007 Fig5 doesn't look log-linear to me)
# The following is the straight-line (in semi-log coords) but the results aren't that great
# self.latest_pc = -math.exp(math.log(-self.Pc_Del[-1]) + (s - self.S_lDel[-1]) * (math.log(-self.Pc_Del[-2]) - math.log(-self.Pc_Del[-1])) / (self.S_lDel[-2] - self.S_lDel[-1]))
# Instead, I interpolate between the previous first-order wetting and second-order trying curves
if self.upper_extension_type == "power":
s_to_use = self.sl_wetting_Del[0] + (1.0 - self.sl_wetting_Del[0]) * (s - self.S_lDel[0]) / (1.0 - self.S_lDel[0])
else:
s_to_use = self.sl_wetting_Del[0] + (self.lc_Del[0] - self.sl_wetting_Del[0]) * (s - self.S_lDel[0]) / (self.lc_Del[0] - self.S_lDel[0])
first_order_wetting = pc_fcn(s_to_use, self.S_lmin, self.gr_Del[0], self.alpha_w, self.n_w, self.S_m_wetting_Del[0], -self.absPcmax, self.dPc_dS_m_wetting_Del[0], self.lower_extension_type, self.S_m_upper_Del[0], self.Pcm_upper_Del[0], self.Pcmprime_upper_Del[0], self.upper_extension_type)
s_to_use = self.S_lDel[0] + (s - self.S_lDel[0]) * (self.Sd_Del[1] - self.S_lDel[0]) / (self.S_lDel[1] - self.S_lDel[0])
second_order_drying = pc_fcn(s_to_use, self.S_lmin, 0.0, self.alpha_d, self.n_d, self.S_m_drying, -self.absPcmax, self.dPc_dSm_drying, self.lower_extension_type, 1.0, 0.0, 0.0, "none") # no upper extension on drying curve
if (-first_order_wetting <= 0.0 or -second_order_drying <= 0.0):
self.latest_pc = 0.0;
else:
self.latest_pc = -math.exp(math.log(-first_order_wetting) + (s - self.S_lDel[1]) * (math.log(-second_order_drying) - math.log(-first_order_wetting)) / (self.S_lDel[2] - self.S_lDel[1]))
def relperms(self, saturations, reinit_tp):
''' Computes liquid and gas relperms given the saturation history contained in the "saturations" list. if reinit_tp == true then the turning points and order are reinitialised to their 'primary drying' values '''
if (reinit_tp):
self.reinit_turning_points()
if (len(saturations) == 0):
return [[], []]
self.set_relperm_vals(saturations[0])
relperm_vals = [[self.latest_kl, self.latest_kg]]
self.prev_saturation = saturations[0]
self.define_order(saturations[0])
for s in saturations[1:]:
self.define_order(s)
self.set_relperm_vals(s)
relperm_vals.append([self.latest_kl, self.latest_kg])
return relperm_vals
def set_relperm_vals(self, s):
'''
Sets self.latest_kl and self.latest_kg, which are the liquid and gas relative permeabilities for the saturation, based on self.order
'''
if self.order == 0:
self.latest_kl = krel_liquid(s, self.S_lr, 0.0, self.S_grmax, 1.0, self.m, self.upper_liquid_param, 0.0, 0.0, 0.0, 0.0)
self.latest_kg = krel_gas(s, self.S_lr, 0.0, self.S_grmax, 1.0, self.m, self.gamma, self.k_rg_max, self.krel_gas_prime)
else:
# min(self.gr_Del[0], self.S_grmax) is used to ensure the gr_Del used is never greater than grmax, which would occur if the turning point were less than S_lr
self.latest_kl = krel_liquid(s, self.S_lr, min(self.gr_Del[0], self.S_grmax), self.S_grmax, self.S_lDel[0], self.m, self.upper_liquid_param, self.kl_begin_Del[0], self.klp_begin_Del[0], self.kl_end_Del[0], self.klp_end_Del[0])
self.latest_kg = krel_gas(s, self.S_lr, min(self.gr_Del[0], self.S_grmax), self.S_grmax, self.S_lDel[0], self.m, self.gamma, self.k_rg_max, self.krel_gas_prime)
def order_and_turning_points(self, saturations):
'''
Returns: (order, list of turning points) for the list of liquid saturations in "saturations". This is used for testing, not for generating capillary pressures. Note: self.order and self.S_lDel are not reinitialised by this method, which allows successive calls to this method to be tested
'''
for s in saturations:
self.define_order(s)
return (self.order, self.S_lDel)
def define_order(self, new_saturation):
'''
Based on new_saturation and self.prev_saturation, defines self.order.
Also defines self.prev_saturation = new_saturation
Also, if the order has changed, adds to self.S_lDel and other things
Also, checks to see if the order can be reduced and self.S_lDel can be reduced
'''
dry2wet = (self.order % 2 == 0) and (new_saturation > self.prev_saturation) # have been reducing saturation ("drying" or "draining") but now are increasing saturation ("wetting" or "imbibing")
wet2dry = (self.order % 2 == 1) and (new_saturation < self.prev_saturation) # have been increasing saturation ("wetting" or "imbibing") but now are decreasing saturation ("drying" or "draining")
if dry2wet or wet2dry:
self.order += 1
self.S_lDel.append(self.prev_saturation)
self.Pcd_Del.append(pc_fcn(self.S_lDel[-1], self.S_lmin, 0.0, self.alpha_d, self.n_d, self.S_m_drying, -self.absPcmax, self.dPc_dSm_drying, self.lower_extension_type, 1.0, 0.0, 0.0, "none")) # no upper extension on drying curve
self.gr_Del.append(self.s_grdel_land(self.S_lDel[-1]))
self.lc_Del.append(1.0 - self.gr_Del[-1])
self.S_m_wetting_Del.append(sat_fcn(-self.absPcmax, self.S_lmin, self.gr_Del[-1], self.alpha_w, self.n_w, 0.0, -self.absPcmax, 0.0, self.lower_extension_type, 1.0, 0.0, 0.0, "none")) # no upper extension at this point of computation
self.dPc_dS_m_wetting_Del.append(pc_fcn_prime(self.S_m_wetting_Del[0], self.S_lmin, self.gr_Del[-1], self.alpha_w, self.n_w, 0.0, -self.absPcmax, 0.0, self.lower_extension_type, 1.0, 0.0, 0.0, "none")) # no upper extension at this point of computation
self.sl_wetting_Del.append(sat_fcn(self.Pcd_Del[-1], self.S_lmin, self.gr_Del[-1], self.alpha_w, self.n_w, self.S_m_wetting_Del[-1], -self.absPcmax, self.dPc_dS_m_wetting_Del[-1], self.lower_extension_type, 1.0, 0.0, 0.0, "none")) # no upper extension at this point of computation
self.Pc_Del.append(self.latest_pc)
self.Sd_Del.append(sat_fcn(self.Pc_Del[-1], self.S_lmin, 0.0, self.alpha_d, self.n_d, self.S_m_drying, -self.absPcmax, self.dPc_dSm_drying, self.lower_extension_type, 1.0, 0.0, 0.0, "none")) # no upper extension on drying curve
self.S_m_upper_Del.append(self.upper_ratio * (1.0 - self.gr_Del[-1]))
self.Pcm_upper_Del.append(pc_fcn(self.S_m_upper_Del[-1], self.S_lmin, self.gr_Del[-1], self.alpha_w, self.n_w, self.S_m_wetting_Del[-1], -self.absPcmax, self.dPc_dS_m_wetting_Del[-1], self.lower_extension_type, self.S_m_upper_Del[-1], 0.0, 0.0, self.upper_extension_type))
self.Pcmprime_upper_Del.append(pc_fcn_prime(self.S_m_upper_Del[-1], self.S_lmin, self.gr_Del[-1], self.alpha_w, self.n_w, self.S_m_wetting_Del[-1], -self.absPcmax, self.dPc_dS_m_wetting_Del[-1], self.lower_extension_type, self.S_m_upper_Del[-1], self.Pcm_upper_Del[-1], 0.0, self.upper_extension_type))
# Relperm stuff
# only use [0] because there is only one wetting curve
if (self.S_lDel[0] < self.S_lr): # deal with the case that "there is no hysteresis along the extension"
self.kl_begin_Del.append(krel_liquid(self.upper_liquid_param * (1.0 - self.S_grmax), self.S_lr, self.S_grmax, self.S_grmax, self.S_lr, self.m, self.upper_liquid_param, 0.0, 0.0, 0.0, 0.0))
self.klp_begin_Del.append(krel_liquid_prime(self.upper_liquid_param * (1.0 - self.S_grmax), self.S_lr, self.S_grmax, self.S_grmax, self.S_lr, self.m, self.upper_liquid_param))
# Doughty2008: self.kl_end_Del.append(krel_liquid((2.0 - self.upper_liquid_param) * (1.0 - self.S_grmax), self.S_lr, self.S_grmax, self.S_grmax, self.S_lr, self.m, self.upper_liquid_param, 0.0, 0.0, 0.0, 0.0))
# Doughty2008: self.klp_end_Del.append(krel_liquid_prime((2.0 - self.upper_liquid_param) * (1.0 - self.S_grmax), self.S_lr, self.S_grmax, self.S_grmax, self.S_lr, self.m, self.upper_liquid_param))
self.kl_end_Del.append(krel_liquid(1.0 - 0.5 * self.S_grmax, self.S_lr, self.S_grmax, self.S_grmax, self.S_lr, self.m, self.upper_liquid_param, 0.0, 0.0, 0.0, 0.0))
self.klp_end_Del.append(krel_liquid_prime(1.0 - 0.5 * self.S_grmax, self.S_lr, self.S_grmax, self.S_grmax, self.S_lr, self.m, self.upper_liquid_param))
else:
self.kl_begin_Del.append(krel_liquid(self.upper_liquid_param * self.lc_Del[0], self.S_lr, self.gr_Del[0], self.S_grmax, self.S_lDel[0], self.m, self.upper_liquid_param, 0.0, 0.0, 0.0, 0.0))
self.klp_begin_Del.append(krel_liquid_prime(self.upper_liquid_param * self.lc_Del[0], self.S_lr, self.gr_Del[0], self.S_grmax, self.S_lDel[0], self.m, self.upper_liquid_param))
# Doughty2008: self.kl_end_Del.append(krel_liquid((2.0 - self.upper_liquid_param) * self.lc_Del[0], self.S_lr, self.gr_Del[0], self.S_grmax, self.S_lDel[0], self.m, self.upper_liquid_param, 0.0, 0.0, 0.0, 0.0))
# Doughty2008: self.klp_end_Del.append(krel_liquid_prime((2.0 - self.upper_liquid_param) * self.lc_Del[0], self.S_lr, self.gr_Del[0], self.S_grmax, self.S_lDel[0], self.m, self.upper_liquid_param))
self.kl_end_Del.append(krel_liquid(1.0 - 0.5 * self.gr_Del[0], self.S_lr, self.gr_Del[0], self.S_grmax, self.S_lDel[0], self.m, self.upper_liquid_param, 0.0, 0.0, 0.0, 0.0))
self.klp_end_Del.append(krel_liquid_prime(1.0 - 0.5 * self.gr_Del[0], self.S_lr, self.gr_Del[0], self.S_grmax, self.S_lDel[0], self.m, self.upper_liquid_param))
can_reduce = self.order > 1
while (can_reduce):
can_reduce = False
reducing_and_low_sat = (self.order % 2 == 0) and (new_saturation <= self.S_lDel[-2]) # are reducing saturation and new_saturation is <= a previously-encountered turning point
increasing_and_high_sat = (self.order % 2 == 1) and (new_saturation >= self.S_lDel[-2]) # are increasing saturation and new_saturation is >= a previously-encountered turning point
if reducing_and_low_sat or increasing_and_high_sat:
self.order -= 2
for i in range(2):
self.S_lDel.pop()
self.Pcd_Del.pop()
self.gr_Del.pop()
self.lc_Del.pop()
self.sl_wetting_Del.pop()
self.Pc_Del.pop()
self.Sd_Del.pop()
self.S_m_wetting_Del.pop()
self.dPc_dS_m_wetting_Del.pop()
self.S_m_upper_Del.pop()
self.Pcm_upper_Del.pop()
self.Pcmprime_upper_Del.pop()
self.kl_begin_Del.pop()
self.klp_begin_Del.pop()
self.kl_end_Del.pop()
self.klp_end_Del.pop()
can_reduce = self.order > 1
self.prev_saturation = new_saturation
Slmin = 0.1
Slr = 0.2 # krel_liquid(Slr) = 0, krel_gas(Slr) = krgmax. Slr >= Slmin
m = 0.9 # exponent in liquid relperm < 1
gamma = 0.3 # exponent in gas relperm (0.3 -> 0.5 is reasonable)
krgmax = 0.8 # value of gas relperm when S_liquid = S_lr. <=1
relperm_gas_extension = "linear_like" # the cubic is more like a line between (0, 1) and (Slr, krgmax)
Sgrmax = 0.3 # must be > 0 to ensure cubic spline works OK on liquid relperm curve
alphad = 5.0
alphaw = 10.0
nd = 1.7
nw = 1.9 # 1.1 is good for testing
absPcmax = 1.5 # value at which extension starts. Because alphad!=alphaw and nd!=nw this causes the drying curve to exceed the wetting curve! MOOSE note: for the quadratic extension, TOUGH demands that users specify |Pc(S=0)| which might be more user-friendly
lower_extension_type = "exponential"
upper_ratio = 0.9 # upper extension will start at S_l = 0.9 * (1 - Sgrmax) for the primary wetting curve
upper_extension_type = "power"
upper_liquid_param = 0.9 # cubic spline is between upper_liquid_param * (1 - Sgr_del) and 1 - 0.5 * Sgr_del. NOTE: Doughty2008 sets this upper number to (2.0 - upper_liquid_param) * (1 - Sgr_del), but that can be >1 and i don't think it gives such a good result
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, absPcmax, lower_extension_type, upper_ratio, upper_extension_type, relperm_gas_extension)
#### Testing of calculating order and the turning points
(order, S_ldel) = hys.order_and_turning_points([1.0])
assert order == 0 and len(S_ldel) == 0
(order, S_ldel) = hys.order_and_turning_points([0.9, 0.8, 0.7, 0.6, 0.5])
assert order == 0 and len(S_ldel) == 0
(order, S_ldel) = hys.order_and_turning_points([0.6, 0.7])
assert order == 1 and len(S_ldel) == 1 and S_ldel[0] == 0.5
(order, S_ldel) = hys.order_and_turning_points([0.7, 0.8, 0.9])
assert order == 1 and len(S_ldel) == 1 and S_ldel[0] == 0.5
(order, S_ldel) = hys.order_and_turning_points([0.8, 0.7, 0.6])
assert order == 2 and len(S_ldel) == 2 and S_ldel[0] == 0.5 and S_ldel[1] == 0.9
(order, S_ldel) = hys.order_and_turning_points([0.5])
assert order == 0 and len(S_ldel) == 0
(order, S_ldel) = hys.order_and_turning_points([0.4])
assert order == 0 and len(S_ldel) == 0
(order, S_ldel) = hys.order_and_turning_points([0.9, 0.6, 0.8, 0.7])
assert order == 4 and len(S_ldel) == 4 and S_ldel[0] == 0.4 and S_ldel[1] == 0.9 and S_ldel[2] == 0.6 and S_ldel[3] == 0.8
(order, S_ldel) = hys.order_and_turning_points([0.85])
assert order == 3 and len(S_ldel) == 3 and S_ldel[0] == 0.4 and S_ldel[1] == 0.9 and S_ldel[2] == 0.6
(order, S_ldel) = hys.order_and_turning_points([0.2])
assert order == 0 and len(S_ldel) == 0
eps = 1E-8
assert cubic(1, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0) == 2.0
deriv = (cubic(1.0 + eps, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0) - cubic(1.0 - eps, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0)) / 2.0 / eps
assert deriv < 3 + 10 * eps and deriv > 3 - 10 * eps
assert cubic(4, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0) == 5.0
deriv = (cubic(4.0 + eps, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0) - cubic(4.0 - eps, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0)) / 2.0 / eps
assert deriv < 6 + 10 * eps and deriv > 6 - 10 * eps
#### Plot the primary drying and wetting curves
plt.figure(0)
svals = [0.001 * i for i in range(1000)]
plt.semilogy(svals, hys.primary_drying(svals), 'r-', label = "Drying")
plt.semilogy(svals, hys.primary_wetting(svals), 'b-', label = "Wetting")
plt.semilogy([Slmin, Slmin], [0, max(hys.primary_drying(svals))], 'k:')
plt.text(Slmin, 0.01, "$S_{lmin}$", horizontalalignment='center', verticalalignment='top')
plt.semilogy([1 - Sgrmax, 1 - Sgrmax], [0, max(hys.primary_drying(svals))], 'k:')
plt.text(1 - Sgrmax, 0.01, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 10])
plt.legend()
plt.tight_layout()
plt.title("Primary drying and wetting curves")
#### Plot hysteretic curves
plt.figure(1)
svals = [0.001 * i for i in range(1000)]
plt.semilogy(svals, hys.primary_drying(svals), 'r-', label = "Primary drying")
plt.semilogy(svals, hys.primary_wetting(svals), 'b-', label = "Primary wetting")
plt.semilogy([Slmin, Slmin], [0, max(hys.primary_drying(svals))], 'k:')
plt.text(Slmin, 0.01, "$S_{lmin}$", horizontalalignment='center', verticalalignment='top')
plt.semilogy([1 - Sgrmax, 1 - Sgrmax], [0, max(hys.primary_drying(svals))], 'k:')
plt.text(1 - Sgrmax, 0.01, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.5 is the turning point in the following
svals = [1 - 0.001 * i for i in range(500)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "hysteretic 0")
plt.semilogy(svals[-1], pcvals[-1], 'rs', label = "turning point 0")
svals = [0.5 + 0.001 * i for i in range(500)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-', label = "hysteretic 1")
plt.semilogy([1.0 - hys.s_grdel_land(0.5)], [0.01], 'gs', label = "$1 - S_{gr}^{\Delta}$")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 10])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic capillary pressure: first order")
#### Plot hysteretic curves
plt.figure(2)
svals = [0.001 * i for i in range(1000)]
plt.semilogy(svals, hys.primary_drying(svals), 'r-', label = "Primary drying")
plt.semilogy(svals, hys.primary_wetting(svals), 'b-', label = "Primary wetting")
plt.semilogy([Slmin, Slmin], [0, max(hys.primary_drying(svals))], 'k:')
plt.text(Slmin, 0.01, "$S_{lmin}$", horizontalalignment='center', verticalalignment='top')
plt.semilogy([1 - Sgrmax, 1 - Sgrmax], [0, max(hys.primary_drying(svals))], 'k:')
plt.text(1 - Sgrmax, 0.01, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.2 and 0.65 are the turning points in the following
svals = [1 - 0.001 * i for i in range(800)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "hysteretic 0")
plt.semilogy(svals[-1], pcvals[-1], 'rs', label = "turning point 0")
svals = [0.2 + 0.001 * i for i in range(450)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-', label = "hysteretic 1")
plt.semilogy(svals[-1], pcvals[-1], 'r.', label = "turning point 1")
svals = [0.65 - 0.001 * i for i in range(550)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g:', label = "hysteretic 2")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 10])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic capillary pressure: second order")
#### Plot hysteretic curves
plt.figure(3)
svals = [0.001 * i for i in range(1000)]
plt.semilogy(svals, hys.primary_drying(svals), 'r-', label = "Primary drying")
plt.semilogy(svals, hys.primary_wetting(svals), 'b-', label = "Primary wetting")
plt.semilogy([Slmin, Slmin], [0, max(hys.primary_drying(svals))], 'k:')
plt.text(Slmin, 0.01, "$S_{lmin}$", horizontalalignment='center', verticalalignment='top')
plt.semilogy([1 - Sgrmax, 1 - Sgrmax], [0, max(hys.primary_drying(svals))], 'k:')
plt.text(1 - Sgrmax, 0.01, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.2 and 0.65 and 0.25 are the turning points in the following
svals = [1 - 0.001 * i for i in range(800)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "hysteretic 0")
svals = [0.2 + 0.001 * i for i in range(450)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-', label = "hysteretic 1")
svals = [0.65 - 0.001 * i for i in range(400)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g:', label = "hysteretic 2")
svals = [0.25 + 0.001 * i for i in range(500)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-.', label = "hysteretic 3")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 10])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic capillary pressure: third order")
#plt.savefig("cp_3.png")
#### Plot hysteretic curves
plt.figure(4)
svals = [0.001 * i for i in range(1000)]
plt.semilogy(svals, hys.primary_drying(svals), 'r-', label = "Primary drying")
plt.semilogy(svals, hys.primary_wetting(svals), 'b-', label = "Primary wetting")
plt.semilogy([Slmin, Slmin], [0, max(hys.primary_drying(svals))], 'k:')
plt.text(Slmin, 0.01, "$S_{lmin}$", horizontalalignment='center', verticalalignment='top')
plt.semilogy([1 - Sgrmax, 1 - Sgrmax], [0, max(hys.primary_drying(svals))], 'k:')
plt.text(1 - Sgrmax, 0.01, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.2 and 0.65 and 0.25 and 0.55 are the turning points in the following
svals = [1 - 0.001 * i for i in range(800)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "hysteretic 0")
svals = [0.2 + 0.001 * i for i in range(450)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-', label = "hysteretic 1")
svals = [0.65 - 0.001 * i for i in range(400)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g:', label = "hysteretic 2")
svals = [0.25 + 0.001 * i for i in range(300)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-.', label = "hysteretic 3")
svals = [0.55 - 0.03 * i for i in range(15)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'ko', markerfacecolor = 'none', label = "hysteretic 4, etc")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 10])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic capillary pressure: fourth order")
#### Plot the drying relperm
plt.figure(5)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
plt.plot(svals, hys.relperm_gas_drying(svals), 'r:', label = "Drying, gas")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, 0.0, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, 0.0, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Relative permeability")
plt.grid()
plt.xlim([0, 1])
#plt.ylim([0, 1])
plt.legend()
plt.tight_layout()
plt.title("Drying relative permeability functions")
#### Plot hysteretic curves
plt.figure(6)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
plt.plot(svals, hys.relperm_gas_drying(svals), 'r:', label = "Drying, gas")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, 0.0, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, 0.0, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.5 is the turning point in the following
svals = [1 - 0.001 * i for i in range(500)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[0] for k in kvals], 'g--', label = "liquid 0")
plt.plot(svals, [k[1] for k in kvals], 'b--', label = "gas 0")
plt.plot(svals[-1], kvals[-1][0], 'rs')
plt.plot(svals[-1], kvals[-1][1], 'rs')
svals = [0.5 + 0.001 * i for i in range(500)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[0] for k in kvals], 'g-', label = "liquid 1")
plt.plot(svals, [k[1] for k in kvals], 'b-', label = "gas 1")
plt.plot([1.0 - hys.s_grdel_land(0.5)], [0.0], 'gs', label = "$1 - S_{gr}^{\Delta}$")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Relative permeability")
plt.grid()
plt.xlim([0, 1])
##plt.ylim([0, 1])
plt.legend()
plt.tight_layout()
plt.title("Relative permeability functions")
#### Plot hysteretic liquid relperm
plt.figure(7)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, 0.0, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, 0.0, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.35 and 0.9 and 0.4 and 0.55 are the turning points in the following
svals = [1 - 0.001 * i for i in range(650)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[0] for k in kvals], 'g--', label = "hysteretic 0")
svals = [0.35 + 0.001 * i for i in range(550)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[0] for k in kvals], 'y-', label = "hysteretic 1")
svals = [0.9 - 0.001 * i for i in range(500)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[0] for k in kvals], 'g:', label = "hysteretic 2")
svals = [0.4 + 0.001 * i for i in range(150)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[0] for k in kvals], 'g-.', label = "hysteretic 3")
svals = [0.55 - 0.03 * i for i in range(15)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[0] for k in kvals], 'ko', markerfacecolor = 'none', label = "hysteretic 4, etc")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Liquid relative permeability")
plt.grid()
plt.xlim([0, 1])
#plt.ylim([0, 1])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic liquid relative permeability: fourth order")
#### Plot hysteretic gas relperm
plt.figure(8)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_gas_drying(svals), 'r-', label = "Drying")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, 0.0, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, 0.0, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.25 and 0.9 and 0.4 and 0.55 are the turning points in the following
svals = [1 - 0.001 * i for i in range(750)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[1] for k in kvals], 'g--', label = "hysteretic 0")
svals = [0.25 + 0.001 * i for i in range(650)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[1] for k in kvals], 'y-', label = "hysteretic 1")
svals = [0.9 - 0.001 * i for i in range(500)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[1] for k in kvals], 'g:', label = "hysteretic 2")
svals = [0.4 + 0.001 * i for i in range(150)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[1] for k in kvals], 'g-.', label = "hysteretic 3")
svals = [0.55 - 0.03 * i for i in range(15)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[1] for k in kvals], 'ko', markerfacecolor = 'none', label = "hysteretic 4, etc")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Gas relative permeability")
plt.grid()
plt.xlim([0, 1])
#plt.ylim([0, 1])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic gas relative permeability: fourth order")
#### Plot hysteretic gas relperm
plt.figure(9)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_gas_drying(svals), 'r-', label = "Drying")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, 0.0, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, 0.0, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.6 is the turning points in the following
svals = [1 - 0.001 * i for i in range(400)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[1] for k in kvals], 'g--', label = "hysteretic 0")
svals = [0.6 + 0.001 * i for i in range(400)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[1] for k in kvals], 'y-', label = "hysteretic 1")
plt.plot([1.0 - hys.s_grdel_land(0.6)], [0.0], 'gs', label = "$1 - S_{gr}^{\Delta}$")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Gas relative permeability")
plt.grid()
plt.xlim([0, 1])
#plt.ylim([0, 1])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic gas relative permeability: first order")
#### Plot hysteretic gas relperm
plt.figure(10)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_gas_drying(svals), 'r-', label = "Drying")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, 0.0, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, 0.0, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.05 is the turning points in the following
svals = [1 - 0.001 * i for i in range(950)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[1] for k in kvals], 'g--', label = "hysteretic 0")
svals = [0.05 + 0.001 * i for i in range(950)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[1] for k in kvals], 'y-', label = "hysteretic 1")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Gas relative permeability")
plt.grid()
plt.xlim([0, 1])
#plt.ylim([0, 1])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic gas relative permeability: first order")
#### Plot hysteretic water relperm
plt.figure(11)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, 0.0, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, 0.0, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.05 is the turning points in the following
svals = [1 - 0.001 * i for i in range(950)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[0] for k in kvals], 'g--', label = "hysteretic 0")
svals = [0.05 + 0.001 * i for i in range(950)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[0] for k in kvals], 'y-', label = "hysteretic 1")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Liquid relative permeability")
plt.grid()
plt.xlim([0, 1])
#plt.ylim([0, 1])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic water relative permeability: first order")
#### Plot hysteretic curves to output nice figures for MOOSE documentation
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, 10.0, lower_extension_type, 0.999, upper_extension_type, relperm_gas_extension)
plt.figure(12)
svals = [0.001 * i for i in range(1000)]
plt.semilogy(svals, hys.primary_drying(svals), 'r-', label = "Primary drying")
plt.semilogy(svals, hys.primary_wetting(svals), 'b-', label = "Primary wetting")
# 0.18 and 0.65 and 0.25 are the turning points in the following
svals = [1 - 0.001 * i for i in range(820)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "order 0")
plt.semilogy([0.18, 0.18], [0, pcvals[-1]], 'k:')
plt.text(0.18, 0.01, "TP$_{0}\,$", horizontalalignment='right', verticalalignment='bottom', fontsize = 16)
svals = [0.18 + 0.001 * i for i in range(470)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-', label = "order 1")
plt.semilogy([0.65, 0.65], [0, pcvals[-1]], 'k:')
plt.text(0.65, 0.01, "TP$_{1}$", horizontalalignment='right', verticalalignment='bottom', fontsize = 16)
svals = [0.65 - 0.001 * i for i in range(400)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g:', label = "order 2")
plt.semilogy([0.25, 0.25], [0, pcvals[-1]], 'k:')
plt.text(0.25, 0.01, "TP$_{2}\,$", horizontalalignment='left', verticalalignment='bottom', fontsize = 16)
svals = [0.25 + 0.001 * i for i in range(300)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-.', label = "order 3")
plt.xlabel("Liquid saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 10])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic capillary pressure")
#plt.savefig("hysteretic_order.png")
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, 10.0, lower_extension_type, 0.999, upper_extension_type, relperm_gas_extension)
plt.figure(13)
svals = [0.001 * i for i in range(1000)]
plt.semilogy(svals, hys.primary_drying(svals), 'r-', label = "Primary drying")
plt.semilogy(svals, hys.primary_wetting(svals), 'b-', label = "Primary wetting")
svals = [1 - 0.001 * i for i in range(400)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "(a) order 0")
svals = [0.6 + 0.001 * i for i in range(300)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-', label = "(b) order 1")
svals = [0.9 - 0.001 * i for i in range(300)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g:', label = "(c) order 2")
svals = [0.6 - 0.001 * i for i in range(300)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g--', label = "(c) order 0")
svals = [0.3 + 0.001 * i for i in range(400)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'y-', label = "(d) order 1")
svals = [0.7 - 0.001 * i for i in range(400)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'y:', label = "(e) order 2")
svals = [0.3 - 0.001 * i for i in range(100)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'y--', label = "(e) order 0")
svals = [0.2 + 0.001 * i for i in range(200)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'c-', label = "(f) order 1")
svals = [0.4 - 0.001 * i for i in range(200)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'c:', label = "(g) order 2")
svals = [0.2 - 0.001 * i for i in range(100)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'c--', label = "(g) order 0")
plt.xlabel("Liquid saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 10])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic capillary pressure")
#plt.savefig("hysteretic_order_012.png")
plt.figure(14)
svals = [0.001 * i for i in range(1000)]
this_Slmin = 0.1
this_Sgrmax = 0.15
this_n = 3.0
plt.semilogy([this_Slmin, this_Slmin], [0, 1E5], 'k:')
plt.text(this_Slmin, 0.01, "$S_{l,min}$", horizontalalignment='center', verticalalignment='top')
plt.semilogy([1 - this_Sgrmax, 1 - this_Sgrmax], [0, 1E5], 'k:')
plt.text(1 - this_Sgrmax, 0.01, "$1 - S_{gr}^{max}$", horizontalalignment='center', verticalalignment='top')
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, this_Slmin, this_Sgrmax, alphad, alphad, this_n, this_n, 1.1, "quadratic", 0.999999, "none", relperm_gas_extension)
plt.semilogy(svals, hys.primary_wetting(svals), 'r--', label = "Quadratic lower")
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, this_Slmin, this_Sgrmax, alphad, alphad, this_n, this_n, 1.1, "exponential", 0.9999999, "none", relperm_gas_extension)
plt.semilogy(svals, hys.primary_wetting(svals), 'g-.', label = "Exponential lower")
xval = hys.S_m_wetting
plt.semilogy([0, xval], [1.1, 1.1], 'k:')
plt.text(xval, 1.1, " " + "$P_{c}^{max}$", horizontalalignment='left', verticalalignment='center')
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, this_Slmin, this_Sgrmax, alphad, alphad, this_n, this_n, 1E4, "none", 0.9, "power", relperm_gas_extension)
plt.semilogy(svals, hys.primary_wetting(svals), 'c--', label = "Power upper")
plt.semilogy([0.9 * (1 - this_Sgrmax), 0.9 * (1 - this_Sgrmax)], [0, hys.primary_wetting([0.9 * (1 - this_Sgrmax)])[0]], 'k:')
plt.text(0.9 * (1 - this_Sgrmax), hys.primary_wetting([0.9 * (1 - this_Sgrmax)])[0], " " + "$0.9 (1 - S_{gr}^{max})$", horizontalalignment='center', verticalalignment='bottom', rotation=90)
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, this_Slmin, this_Sgrmax, alphad, alphad, this_n, this_n, 1000.0, "none", 0.9999999, "none", relperm_gas_extension)
plt.semilogy(svals, hys.primary_wetting(svals), 'k-', label = "No extension")
plt.xlabel("Liquid saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xticks([0, 0.3, 0.5, 0.7, 1], ["0", "0.3", "0.5", "0.7", "1"])
plt.xlim([0, 1])
plt.ylim([0.01, 10])
plt.legend()
plt.tight_layout()
plt.title("Primary wetting capillary pressure with various extensions")
#plt.savefig("hysteretic_cap_extensions.png")
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, 10.0, lower_extension_type, 0.999, "none", relperm_gas_extension)
plt.figure(15)
svals = [0.001 * i for i in range(1000)]
plt.semilogy(svals, hys.primary_drying(svals), 'r-', label = "Primary drying")
plt.semilogy(svals, hys.primary_wetting(svals), 'b-', label = "Primary wetting")
plt.text(1 - Sgrmax, 0.01, "$1 - S_{gr}^{max}$", horizontalalignment='center', verticalalignment='top')
plt.text(Slmin, 0.01, "$S_{l,min}$", horizontalalignment='center', verticalalignment='top')
svals = [1 - 0.001 * i for i in range(100)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
svals = [0.9 + 0.001 * i for i in range(100)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-', label = "order 1")
lc = [hys.lc_Del[0]]
svals = [1 - 0.001 * i for i in range(400)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
svals = [0.6 + 0.001 * i for i in range(400)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
lc.append(hys.lc_Del[0])
plt.semilogy(svals, pcvals, 'g-', label = "order 1")
svals = [1 - 0.001 * i for i in range(800)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
svals = [0.2 + 0.001 * i for i in range(800)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
lc.append(hys.lc_Del[0])
plt.semilogy(svals, pcvals, 'g-', label = "order 1")
plt.semilogy([lc[0]], [0.01], 'cs', label = "$1 - S_{gr}^{\Delta}$")
plt.semilogy([lc[1]], [0.01], 'cs')
plt.semilogy([lc[2]], [0.01], 'cs')
plt.xlabel("Liquid saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 10])
plt.legend()
plt.tight_layout()
plt.title("Various first-order wetting curves")
#plt.savefig("hysteretic_different_tps.png")
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, 10.0, lower_extension_type, 0.999, "none", relperm_gas_extension)
plt.figure(16)
svals = [0.001 * i for i in range(1000)]
plt.semilogy(svals, hys.primary_drying(svals), 'r-', label = "Primary drying")
plt.semilogy([Slmin, Slmin], [0.01, 1E4], 'k:')
plt.text(Slmin, 0.01, "$S_{l,min}$", horizontalalignment='center', verticalalignment='top')
svals = [1 - 0.001 * i for i in range(600)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
svals = [0.4 + 0.001 * i for i in range(600)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-', label = "First-order curve")
plt.semilogy([0.4, 0.4, 0], [0.01, -hys.Pcd_Del[0], -hys.Pcd_Del[0]], 'k:')
plt.semilogy([hys.sl_wetting_Del[0], hys.sl_wetting_Del[0], 0], [0.01, -hys.Pcd_Del[0], -hys.Pcd_Del[0]], 'k:')
plt.text(0.4, 0.01, "TP$_{0}$", horizontalalignment='center', verticalalignment='top')
plt.text(hys.sl_wetting_Del[0], 0.01, "$S_{l,wet}^{\Delta}$", horizontalalignment='center', verticalalignment='top')
plt.text(0, -hys.Pcd_Del[0], "$P_{c}^{\Delta}$", horizontalalignment='right', verticalalignment='center')
gr_del = hys.gr_Del[0]
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, gr_del, alphad, alphaw, nd, nw, 10.0, lower_extension_type, 0.999, "none", relperm_gas_extension)
svals = [0.001 * i for i in range(1000)]
plt.semilogy(svals, hys.primary_wetting(svals), 'b-', label = "Wetting curve using $S_{gr}^{\Delta}$")
plt.semilogy([1 - gr_del], [0.01], 'cs', label = "$1 - S_{gr}^{\Delta}$")
plt.xlabel("Liquid saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xticks([0, 0.3, 0.5, 0.7, 1], ["0", "0.3", "0.5", "0.7", "1"])
plt.yticks([0.01, 0.1, 10], ["0.01", "0.1", "10"])
plt.xlim([0, 1])
plt.ylim([0.01, 10])
plt.legend()
plt.tight_layout()
plt.title("Construction of the first-order wetting curve")
#plt.savefig("hysteretic_order1.png")
#### Plot drying and wetting relperms with no extensions
upper_liquid_param = 1.0 # below ignores upper extension on liquid wetting curve
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, absPcmax, lower_extension_type, upper_ratio, upper_extension_type, "no_relperm_gas_extension")
plt.figure(17)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
svals = [0.001 * i for i in range(200, 1001)]
plt.plot(svals, hys.relperm_gas_drying(svals), 'r:', label = "Drying, gas")
plt.plot([0, 0.2], [krgmax, krgmax], 'r:')
plt.text(0, krgmax, "$k_{r,g}^{max}$", horizontalalignment='right', verticalalignment='center')
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, -1E-2, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, -1E-2, "$1 - S_{gr}^{max}$", horizontalalignment='center', verticalalignment='top')
svals = [1 - 0.001 * i for i in range(800)]
kvals = hys.relperms(svals, True)
plt.plot(svals[-1], kvals[-1][0], 'rs', label = "turning point")
svals = [0.2 + 0.001 * i for i in range(500)]
kvals = hys.relperms(svals, False)
svals += [0.2 + 0.001 * i for i in range(500, 800)]
rpvals = [k[0] for k in kvals] + [1.0 for s in range(500, 800)]
plt.plot(svals, rpvals, 'b-', label = "first-order, liquid")
svals = [1 - 0.001 * i for i in range(800)]
kvals = hys.relperms(svals, True)
plt.plot(svals[-1], kvals[-1][1], 'rs')
svals = [0.2 + 0.001 * i for i in range(800)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[1] for k in kvals], 'b:', label = "first-order, gas")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Relative permeability")
plt.grid()
plt.xticks([0, 0.5, 1], ["0", "0.5", "1"])
plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1], ["0", "0.2", "0.4", "0.6", "", "1"])
plt.xlim([0, 1])
plt.ylim([-1E-2, 1 + 1E-2])
plt.legend()
#plt.tight_layout()
plt.title("Unmodified relative permeability functions")
#plt.savefig("hysteretic_krel_unextended.png")
upper_liquid_param = 0.9
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, absPcmax, lower_extension_type, upper_ratio, upper_extension_type, relperm_gas_extension)
plt.figure(18)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
plt.plot(svals, hys.relperm_gas_drying(svals), 'r:', label = "Drying, gas")
plt.text(0, krgmax, "$k_{r,g}^{max}$", horizontalalignment='right', verticalalignment='center')
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, -1E-2, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, -1E-2, "$1 - S_{gr}^{max}$", horizontalalignment='center', verticalalignment='top')
svals = [1 - 0.001 * i for i in range(800)]
kvals = hys.relperms(svals, True)
plt.plot(svals[-1], kvals[-1][0], 'rs', label = "turning point")
svals = [0.2 + 0.001 * i for i in range(800)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[0] for k in kvals], 'b-', label = "first-order, liquid")
svals = [1 - 0.001 * i for i in range(800)]
kvals = hys.relperms(svals, True)
plt.plot(svals[-1], kvals[-1][1], 'rs')
svals = [0.2 + 0.001 * i for i in range(800)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[1] for k in kvals], 'b:', label = "first-order, gas")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Relative permeability")
plt.grid()
plt.xticks([0, 0.5, 1], ["0", "0.5", "1"])
plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1], ["0", "0.2", "0.4", "0.6", "", "1"])
plt.xlim([0, 1])
plt.ylim([-1E-2, 1 + 1E-3])
plt.legend()
plt.tight_layout()
plt.title("Relative permeability functions with extensions")
#plt.savefig("hysteretic_krel_extended.png")
plt.figure(19)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
plt.plot(svals, hys.relperm_gas_drying(svals), 'r:', label = "Drying, gas")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, -1E-2, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, -1E-2, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.5 is the turning point in the following
svals = [1 - 0.001 * i for i in range(500)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[0] for k in kvals], 'g--', label = "liquid 0")
plt.plot(svals, [k[1] for k in kvals], 'b--', label = "gas 0")
plt.plot(svals[-1], kvals[-1][0], 'rs')
plt.plot(svals[-1], kvals[-1][1], 'rs')
svals = [0.5 + 0.001 * i for i in range(500)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[0] for k in kvals], 'g-', label = "liquid 1")
plt.plot(svals, [k[1] for k in kvals], 'b-', label = "gas 1")
plt.plot([1.0 - hys.s_grdel_land(0.5)], [0.0], 'gs', label = "$1 - S_{gr}^{\Delta}$")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Relative permeability")
plt.xticks([0, 0.5, 1], ["0", "0.5", "1"])
plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1], ["0", "0.2", "0.4", "0.6", "", "1"])
plt.grid()
plt.xlim([0, 1])
plt.ylim([-1E-2, 1 + 1E-3])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic relative permeability functions")
#plt.savefig("hysteretic_krel_example_1.png")
plt.figure(20)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
plt.plot(svals, hys.relperm_gas_drying(svals), 'r:', label = "Drying, gas")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, -1E-2, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, -1E-2, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
# 0.1 is the turning point in the following
svals = [1 - 0.001 * i for i in range(900)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[0] for k in kvals], 'g--', label = "liquid 0")
plt.plot(svals, [k[1] for k in kvals], 'b--', label = "gas 0")
plt.plot(svals[-1], kvals[-1][0], 'rs')
plt.plot(svals[-1], kvals[-1][1], 'rs')
svals = [0.1 + 0.001 * i for i in range(900)]
kvals = hys.relperms(svals, False)
plt.plot(svals, [k[0] for k in kvals], 'g-', label = "liquid 1")
plt.plot(svals, [k[1] for k in kvals], 'b-', label = "gas 1")
plt.plot([1.0 - min(hys.gr_Del[0], hys.S_grmax)], [0.0], 'gs', label = "$1 - S_{gr}^{\Delta}$")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Relative permeability")
plt.xticks([0, 0.5, 1], ["0", "0.5", "1"])
plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1], ["0", "0.2", "0.4", "0.6", "", "1"])
plt.grid()
plt.xlim([0, 1])
plt.ylim([-1E-2, 1 + 1E-3])
plt.legend()
plt.tight_layout()
plt.title("Hysteretic relative permeability functions")
#plt.savefig("hysteretic_krel_example_2.png")
# For comparison with the vary_sat series of tests
Slmin = 0.1
Slr = 0.2
Sgrmax = 0.3
alphad = 10.0
alphaw = 7.0
nd = 1.5
nw = 1.9
absPcmax = 12
lower_extension_type = "quadratic"
upper_ratio = 0.9
upper_extension_type = "power"
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, absPcmax, lower_extension_type, upper_ratio, upper_extension_type, relperm_gas_extension)
plt.figure(21)
svals = [1 - 0.001 * i for i in range(1000)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'b-', label = "primary drying")
f = open("gold/vary_sat_1_out.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
#plt.semilogy([x[3] for x in data], [x[2] for x in data], 'b.')
svals = [1 - 0.001 * i for i in range(1000)] + [0.0 + 0.001 * i for i in range(1000)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "drying then wetting")
f = open("gold/vary_sat_1b.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.semilogy([x[3] for x in data], [x[2] for x in data], 'g+', label="MOOSE")
svals = [1 - 0.001 * i for i in range(800)] + [0.2 + 0.001 * i for i in range(800)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'y-.', label = "drying then first-order wetting")
f = open("gold/vary_sat_1c.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.semilogy([x[3] for x in data], [x[2] for x in data], 'y*', label="MOOSE")
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 40])
plt.legend()
plt.title("Hysteretic capillary pressure examples")
#plt.savefig("../../../doc/content/media/porous_flow/hys_vary_sat_1.png")
plt.figure(22)
svals = [1 - 0.001 * i for i in range(800)] + [0.2 + 0.001 * i for i in range(600)] + [0.8 - 0.001 * i for i in range(800)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'k-', label = 'Expected')
f = open("gold/vary_sat_1d.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.semilogy([x[3] for x in data], [x[2] for x in data], 'b.', label = 'MOOSE')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 40])
plt.legend()
plt.title("Hysteretic capillary pressure 2nd-order example")
#plt.savefig("../../../doc/content/media/porous_flow/hys_vary_sat_1_2ndorder.png")
plt.figure(23)
svals = [1 - 0.001 * i for i in range(800)] + [0.2 + 0.001 * i for i in range(600)] + [0.8 - 0.001 * i for i in range(500)] + [0.3 + 0.001 * i for i in range(700)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'k-', label = 'Expected')
f = open("gold/vary_sat_1e.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.semilogy([x[3] for x in data], [x[2] for x in data], 'b.', label = 'MOOSE')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 40])
plt.legend()
plt.title("Hysteretic capillary pressure 3rd-order example")
#plt.savefig("../../../doc/content/media/porous_flow/hys_vary_sat_1_3rdorder.png")
plt.figure(24)
svals = [1 - 0.001 * i for i in range(1000)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'b-', label = "primary drying")
svals = [1 - 0.001 * i for i in range(905)] + [0.096 + 0.001 * i for i in range(905)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "drying then wetting")
f = open("gold/1phase_out.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.semilogy([x[4] for x in data], [max(-x[3], 0) for x in data], 'g+', label = 'MOOSE')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 40])
plt.legend()
plt.title("Hysteretic capillary pressure examples")
#plt.savefig("../../../doc/content/media/porous_flow/hys_1phase_1.png")
plt.figure(25)
svals = [1 - 0.001 * i for i in range(1000)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'b-', label = "primary drying")
svals = [1 - 0.001 * i for i in range(905)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
svals = [0.096 + 0.001 * i for i in range(705)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'g-', label = "1st-order wetting")
svals = [0.8 - 0.001 * i for i in range(605)]
pcvals = [abs(pc) for pc in hys.pc(svals, False)]
plt.semilogy(svals, pcvals, 'k-', label = "2nd-order drying")
# do not plot hys.pc since the third-order S(Pc) is not an inverse of Pc(S)
#svals = [0.196 + 0.001 * i for i in range(605)]
#pcvals = [abs(pc) for pc in hys.pc(svals, False)]
#plt.semilogy(svals, pcvals, 'y-', label = "3rd order")
#instead, read from a by-hand calculation:
f = open("3rd_order_eg.txt", "r")
data = [list(map(float, line.strip().split())) for line in f.readlines()]
plt.semilogy([x[0] for x in data], [x[1] for x in data], 'y-', label = '3rd order')
f.close()
f = open("gold/1phase_3rd_out.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.semilogy([x[4] for x in data], [max(-x[3], 0) for x in data], 'r+', label = 'MOOSE')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.001, 40])
plt.legend()
plt.title("Hysteretic capillary pressure example: 3rd order")
#plt.savefig("../../../doc/content/media/porous_flow/hys_1phase_3.png")
plt.figure(26)
svals = [1 - 0.001 * i for i in range(1000)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'b-', label = "primary drying")
svals = [1 - 0.001 * i for i in range(465)] + [0.535 + 0.001 * i for i in range(465)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "drying then wetting")
f = open("gold/2phasePP_out.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.semilogy([x[5] for x in data], [x[4] - x[3] for x in data], 'g+', label = 'MOOSE, 2 phase PP')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0.4, 1])
plt.ylim([0.01, 1])
plt.legend()
plt.title("Hysteretic capillary pressure examples (2 phase)")
#plt.savefig("../../../doc/content/media/porous_flow/hys_2phasePP_1.png")
plt.figure(27)
svals = [1 - 0.001 * i for i in range(1000)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'b-', label = "primary drying")
svals = [1 - 0.001 * i for i in range(512)] + [0.512 + 0.001 * i for i in range(440)] + [0.952 - 0.001 * i for i in range(500)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "drying, wetting, re-drying")
f = open("gold/2phasePP_2_out.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.semilogy([x[5] for x in data], [x[4] - x[3] for x in data], 'g+', label = 'MOOSE, 2 phase PP')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 40])
plt.legend()
plt.title("Hysteretic capillary pressure examples (2 phase)")
#plt.savefig("../../../doc/content/media/porous_flow/hys_2phasePP_2.png")
plt.figure(28)
svals = [1 - 0.001 * i for i in range(1000)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'b-', label = "primary drying")
svals = [1 - 0.001 * i for i in range(465)] + [0.535 + 0.001 * i for i in range(465)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "drying then wetting")
f = open("gold/2phasePS_out.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.semilogy([x[5] for x in data], [x[4] - x[3] for x in data], 'g+', label = 'MOOSE, 2 phase PS')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0, 1])
plt.ylim([0.01, 40])
plt.legend()
plt.title("Hysteretic capillary pressure examples (2 phase)")
#plt.savefig("../../../doc/content/media/porous_flow/hys_2phasePS_1.png")
plt.figure(29)
svals = [1 - 0.001 * i for i in range(1000)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'b-', label = "primary drying")
svals = [1 - 0.001 * i for i in range(512)] + [0.512 + 0.001 * i for i in range(440)] + [0.952 - 0.001 * i for i in range(500)]
pcvals = [abs(pc) for pc in hys.pc(svals, True)]
plt.semilogy(svals, pcvals, 'g--', label = "drying, wetting, re-drying")
f = open("gold/2phasePS_2_out.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.semilogy([x[5] for x in data], [x[4] - x[3] for x in data], 'g+', label = 'MOOSE, 2 phase PS')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Capillary pressure ($|P_{c}|$)")
plt.grid("x")
plt.xlim([0.4, 1])
plt.ylim([0.001, 1])
plt.legend()
plt.title("Hysteretic capillary pressure examples (2 phase)")
#plt.savefig("../../../doc/content/media/porous_flow/hys_2phasePS_2.png")
# For comparison with the relperm series of tests
Slmin = 0.1
Sgrmax = 0.3
alphad = 10.0
alphaw = 10.0
nd = 1.5
nw = 1.5
absPcmax = 1E10
lower_extension_type = "quadratic"
upper_extension_type = "power"
Slr = 0.1
Sgrmax = 0.2
m = 0.9
upper_ratio = 0.9
gamma = 0.33
krgmax = 0.8
relperm_gas_extension = "linear_like"
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, absPcmax, lower_extension_type, upper_ratio, upper_extension_type, relperm_gas_extension)
plt.figure(30)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, -1E-2, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
svals = [1 - 0.001 * i for i in range(502)] + [0.498 + 0.001 * i for i in range(502)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[0] for k in kvals], 'g--', label = "expected")
f = open("gold/1phase_relperm_out.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.plot([x[4] for x in data], [x[3] for x in data], 'g+', label = 'MOOSE, 1-phase')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Liquid relative permeability")
plt.xticks([0, 0.4, 0.5, 0.6, 0.7, 0.9, 1], ["0", "0.4", "0.5", "0.6", "0.7", "0.9", "1"])
plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1], ["0", "0.2", "0.4", "0.6", "0.8", "1"])
plt.grid()
plt.xlim([0.4, 1])
plt.ylim([-1E-2, 1 + 1E-3])
plt.legend()
plt.tight_layout()
plt.title("1-phase hysteretic relative permeability")
#plt.savefig("../../../doc/content/media/porous_flow/hys_1phase_relperm.png")
plt.figure(31)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, -1E-2, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.text(1 - Sgrmax, -1E-2, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
svals = [1 - 0.001 * i for i in range(302)] + [0.698 + 0.001 * i for i in range(126)] + [0.824 - 0.001 * i for i in range(802)] + [0.02 + 0.001 * i for i in range(980)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[0] for k in kvals], 'g--', label = "expected")
f = open("gold/1phase_relperm_2_csv.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.plot([x[4] for x in data], [x[3] for x in data], 'g+', label = 'MOOSE, 1-phase')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Liquid relative permeability")
plt.xticks([0, 0.5, 1], ["0", "0.5", "1"])
plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1], ["0", "0.2", "0.4", "0.6", "0.8", "1"])
plt.grid()
plt.xlim([0.0, 1])
plt.ylim([-1E-2, 1 + 1E-3])
plt.legend()
plt.tight_layout()
plt.title("1-phase hysteretic relative permeability: cyclic saturation")
#plt.savefig("../../../doc/content/media/porous_flow/hys_1phase_2_relperm.png")
plt.figure(32)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
plt.plot(svals, hys.relperm_gas_drying(svals), 'r:', label = "Drying, gas")
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, -1E-2, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
svals = [1 - 0.001 * i for i in range(457)] + [0.543 + 0.001 * i for i in range(457)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[0] for k in kvals], 'g--', label = "exepcted, liquid")
plt.plot(svals, [k[1] for k in kvals], 'b--', label = "expected, gas")
plt.plot([1.0 - hys.s_grdel_land(0.543)], [0.0], 'gs', label = "$1 - S_{gr}^{\Delta}$")
f = open("gold/2phasePS_relperm_out.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.plot([x[5] for x in data], [x[4] for x in data], 'g+', label = 'MOOSE, liquid')
plt.plot([x[5] for x in data], [x[3] for x in data], 'b+', label = 'MOOSE, gas')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Relative permeability")
plt.xticks([0, 0.5, 0.6, 0.7, 0.9, 1], ["0", "0.5", "0.6", "0.7", "0.9", "1"])
plt.grid()
plt.xlim([0.5, 1])
plt.ylim([-1E-2, 1 + 1E-3])
plt.legend()
plt.tight_layout()
plt.title("2-phase hysteretic relative permeability")
#plt.savefig("../../../doc/content/media/porous_flow/hys_2phasePS_relperm.png")
Slr = 0.4
krgmax = 0.4
relperm_gas_extension = "linear_like"
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, absPcmax, lower_extension_type, upper_ratio, upper_extension_type, relperm_gas_extension)
plt.figure(33)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
plt.plot(svals, hys.relperm_gas_drying(svals), 'r:', label = "Drying, gas")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, -1E-2, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, -1E-2, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
plt.text(0, krgmax, "$k_{r,g}^{max}$", horizontalalignment='right', verticalalignment='center')
svals = [1 - 0.001 * i for i in range(725)] + [0.275 + 0.001 * i for i in range(725)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[0] for k in kvals], 'g--', label = "exepcted, liquid")
plt.plot(svals, [k[1] for k in kvals], 'b--', label = "expected, gas")
f = open("gold/2phasePS_relperm_2_linear_like.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.plot([x[5] for x in data], [x[4] for x in data], 'g+', label = 'MOOSE, liquid')
plt.plot([x[5] for x in data], [x[3] for x in data], 'b+', label = 'MOOSE, gas')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Relative permeability")
plt.xticks([0, 0.5, 1], ["0", "0.5", "1"])
plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1], ["0", "0.2", "", "0.6", "0.8", "1"])
plt.grid()
plt.xlim([0, 1])
plt.ylim([-1E-2, 1 + 1E-3])
plt.legend()
plt.tight_layout()
plt.title("2-phase hysteretic relative permeability: linear-like extension")
#plt.savefig("../../../doc/content/media/porous_flow/hys_2phasePS_relperm_2_linear_like.png")
relperm_gas_extension = "cubic"
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, absPcmax, lower_extension_type, upper_ratio, upper_extension_type, relperm_gas_extension)
plt.figure(34)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
plt.plot(svals, hys.relperm_gas_drying(svals), 'r:', label = "Drying, gas")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, -1E-2, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, -1E-2, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
plt.text(0, krgmax, "$k_{r,g}^{max}$", horizontalalignment='right', verticalalignment='center')
svals = [1 - 0.001 * i for i in range(725)] + [0.275 + 0.001 * i for i in range(725)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[0] for k in kvals], 'g--', label = "exepcted, liquid")
plt.plot(svals, [k[1] for k in kvals], 'b--', label = "expected, gas")
f = open("gold/2phasePS_relperm_2_cubic.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.plot([x[5] for x in data], [x[4] for x in data], 'g+', label = 'MOOSE, liquid')
plt.plot([x[5] for x in data], [x[3] for x in data], 'b+', label = 'MOOSE, gas')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Relative permeability")
plt.xticks([0, 0.5, 1], ["0", "0.5", "1"])
plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1], ["0", "0.2", "", "0.6", "0.8", "1"])
plt.grid()
plt.xlim([0, 1])
plt.ylim([-1E-2, 1 + 1E-3])
plt.legend()
plt.tight_layout()
plt.title("2-phase hysteretic relative permeability: cubic extension")
#plt.savefig("../../../doc/content/media/porous_flow/hys_2phasePS_relperm_2_cubic.png")
Slr = 0.4
krgmax = 1.0
relperm_gas_extension = "linear_like"
hys = Hys(Slr, m, gamma, krgmax, upper_liquid_param, Slmin, Sgrmax, alphad, alphaw, nd, nw, absPcmax, lower_extension_type, upper_ratio, upper_extension_type, relperm_gas_extension)
plt.figure(35)
svals = [0.001 * i for i in range(1001)]
plt.plot(svals, hys.relperm_liquid_drying(svals), 'r-', label = "Drying, liquid")
plt.plot(svals, hys.relperm_gas_drying(svals), 'r:', label = "Drying, gas")
plt.plot([Slr, Slr], [0, 1], 'k:')
plt.text(Slr, -1E-2, "$S_{lr}$", horizontalalignment='center', verticalalignment='top')
plt.plot([1 - Sgrmax, 1 - Sgrmax], [0, 1.0], 'k:')
plt.text(1 - Sgrmax, -1E-2, "$1 - S_{grmax}$", horizontalalignment='center', verticalalignment='top')
plt.text(0, krgmax, "$k_{r,g}^{max}$", horizontalalignment='right', verticalalignment='center')
svals = [1 - 0.001 * i for i in range(725)] + [0.275 + 0.001 * i for i in range(725)]
kvals = hys.relperms(svals, True)
plt.plot(svals, [k[0] for k in kvals], 'g--', label = "exepcted, liquid")
plt.plot(svals, [k[1] for k in kvals], 'b--', label = "expected, gas")
f = open("gold/2phasePS_relperm_2_none.csv", "r")
data = [list(map(float, line.strip().split(","))) for line in f.readlines()[2:]]
f.close()
plt.plot([x[5] for x in data], [x[4] for x in data], 'g+', label = 'MOOSE, liquid')
plt.plot([x[5] for x in data], [x[3] for x in data], 'b+', label = 'MOOSE, gas')
plt.xlabel("Aqueous saturation ($S_{l}$)")
plt.ylabel("Relative permeability")
plt.xticks([0, 0.5, 1], ["0", "0.5", "1"])
plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1], ["0", "0.2", "0.4", "0.6", "0.8", ""])
plt.grid()
plt.xlim([0, 1])
plt.ylim([-1E-2, 1 + 1E-3])
plt.legend()
plt.tight_layout()
plt.title("2-phase hysteretic relative permeability: no extension")
#plt.savefig("../../../doc/content/media/porous_flow/hys_2phasePS_relperm_2_none.png")
plt.show()
sys.exit(0)
Input file syntax for hysteretic capillarity
To include hysteretic capillarity in an existing (non-hysteretic) input file, the following changes need to be made.
Any capillary-pressure UserObjects, such as PorousFlowCapillaryPressureVG will no-longer be used, so should be removed from the input file for clarity.
A PorousFlowHysteresisOrder Material needs to be included.
A Material that computes the porepressure(s) and saturation(s) needs to be included.
For 1-phase partially-saturated situations, PorousFlow1PhaseP should be removed and replaced by PorousFlow1PhaseHysP. Note the van Genuchten parameter input is for the non-hysteretic version, but for the hysteretic version, where and .
For 2-phase situations using the two porepressures as the primary variables, PorousFlow2PhasePP should be removed and replaced by PorousFlow2PhaseHysPP.
For 2-phase situations using the liquid porepressure and gas saturation as the primary variables, PorousFlow2PhasePS should be removed and replaced by PorousFlow2PhaseHysPS.
An example is:
[hys_order_material]
type = PorousFlowHysteresisOrder
[]
[pc_calculator]
type = PorousFlow1PhaseHysP
alpha_d = 10.0
alpha_w = 7.0
n_d = 1.5
n_w = 1.9
S_l_min = 0.1
S_lr = 0.2
S_gr_max = 0.3
Pc_max = 12.0
high_ratio = 0.9
low_extension_type = quadratic
high_extension_type = power
porepressure = pp
[]
[]
Preliminary example
Before introducing hysteresis into an input file, it might be useful to assess how hysteresis evolves during model evolution using a PorousFlowHystereticInfo Material. This Material does not compute porepressures or saturations, but instead allows users to visualise capillary pressure. Hence, the following input file does not perform a usual PorousFlow simulation, but instead allows a preliminary exploration of hysteresis:
# The saturation is varied with time and the capillary pressure is computed
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
number_fluid_phases = 1
number_fluid_components = 1
porous_flow_vars = ''
[]
[]
[Variables]
[dummy]
[]
[]
[Kernels]
[dummy]
type = TimeDerivative
variable = dummy
[]
[]
[AuxVariables]
[sat]
initial_condition = 1
[]
[hys_order]
family = MONOMIAL
order = CONSTANT
[]
[pc]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[sat_aux]
type = FunctionAux
variable = sat
function = '1 - t'
[]
[hys_order]
type = PorousFlowPropertyAux
variable = hys_order
property = hysteresis_order
[]
[pc]
type = PorousFlowPropertyAux
variable = pc
property = hysteretic_info
[]
[]
[Materials]
[hys_order]
type = PorousFlowHysteresisOrder
[]
[pc_calculator]
type = PorousFlowHystereticInfo
alpha_d = 10.0
alpha_w = 7.0
n_d = 1.5
n_w = 1.9
S_l_min = 0.1
S_lr = 0.2
S_gr_max = 0.3
Pc_max = 12.0
high_ratio = 0.9
low_extension_type = quadratic
high_extension_type = power
sat_var = sat
[]
[]
[Postprocessors]
[hys_order]
type = PointValue
point = '0 0 0'
variable = hys_order
[]
[sat]
type = PointValue
point = '0 0 0'
variable = sat
[]
[pc]
type = PointValue
point = '0 0 0'
variable = pc
[]
[]
[Executioner]
type = Transient
solve_type = Linear
dt = 0.1
end_time = 1
[]
[Outputs]
csv = true
[]
By changing the FunctionAux that controls the saturation, various hysteretic curves may be constructed.
[sat_aux]
type = FunctionAux
variable = sat
function = '1 - t'
[]
Figure 10 shows the results of two hysteretic simulations. Both are initialised at full saturation and are then dried. The first (green curve) dries to before wetting, so follows the primary wetting curve. The second (yellow curve) dries to before wetting, so follows a first-order wetting curve.
Figure 10: The results of two hysteretic simulations. The lines show the expected result (from the python script) while the crosses and asterisks show the MOOSE result.
Careful examination of Figure 10 will reveal a subtle feature of the implementation of hysteresis in PorousFlow. As the system dries and then re-wets, it follows the primary drying curve for the first time-step after the turning point. This is why there is no asterisk at : instead it appears on the primary drying curve at . Similarly, there is no cross around . Consequences of this are discussed in a section below.
Figure 11 results when the FunctionAux is
if(t <= 0.4, 1 - 2 * t, if(t <= 0.7, 2 * t - 0.6, if(t <= 0.95, 0.8 - 2 * (t - 0.7), 0.3 + 2 * (t - 0.95))))
Figure 11: The result of a hysteretic simulation where the system is dried, then re-wet, then dried, then re-wet, so that it follows the zeroth, first, second and third-order curves
Single-phase example
A simulation that simply removes and adds water to a system to observe the hysteretic capillary pressure is explored in this section. The water flux is controlled by the Postprocessor
[flux]
type = FunctionValuePostprocessor
function = 'if(t <= 9, -10, if(t <= 16, 10, if(t <= 22, -10, 10)))'
[]
and the DiracKernel
[DiracKernels]
[pump]
type = PorousFlowPointSourceFromPostprocessor
mass_flux = flux
point = '0.5 0 0'
variable = pp
[]
[]
The remainder of the input file is standard, with the inclusion of the hysteretic capillary pressure:
[hys_order_material]
type = PorousFlowHysteresisOrder
[]
[pc_calculator]
type = PorousFlow1PhaseHysP
alpha_d = 10.0
alpha_w = 7.0
n_d = 1.5
n_w = 1.9
S_l_min = 0.1
S_lr = 0.2
S_gr_max = 0.3
Pc_max = 12.0
high_ratio = 0.9
low_extension_type = quadratic
high_extension_type = power
porepressure = pp
[]
[]
The result is Figure 12.
Figure 12: The result of a single-phase simulation in which an external pump removes and adds water to a porous material in order to observe the hysteretic capillary pressure.
Two-phase example using the PP formulation
A simulation that simply adds gas then removes gas from a 2-phase system to observe the hysteretic capillary pressure is explored in this section. The gas flux is controlled by the Postprocessor
[flux]
type = FunctionValuePostprocessor
function = 'if(t <= 9, 10, -10)'
[]
and the DiracKernel
[DiracKernels]
[pump]
type = PorousFlowPointSourceFromPostprocessor
mass_flux = flux
point = '0.5 0 0'
variable = pp1
[]
[]
The remainder of the input file is standard, with the inclusion of the hysteretic capillary pressure:
[hys_order_material]
type = PorousFlowHysteresisOrder
[]
[pc_calculator]
type = PorousFlow2PhaseHysPP
alpha_d = 10.0
alpha_w = 7.0
n_d = 1.5
n_w = 1.9
S_l_min = 0.1
S_lr = 0.2
S_gr_max = 0.3
Pc_max = 12.0
high_ratio = 0.9
low_extension_type = quadratic
high_extension_type = power
phase0_porepressure = pp0
phase1_porepressure = pp1
[]
[]
The result is Figure 13.
Figure 13: The result of a two-phase simulation using a PP formulation in which an external pump adds to and removes gas from a porous material in order to observe the hysteretic capillary pressure.
Two-phase example using the PS formulation
A simulation that simply adds gas, then removes gas, and adds it again to a 2-phase system to observe the hysteretic capillary pressure is explored in this section. The gas flux is controlled by the Postprocessor
[flux]
type = FunctionValuePostprocessor
function = 'if(t <= 14, 10, if(t <= 25, -10, 10))'
[]
and the DiracKernel
[DiracKernels]
[pump]
type = PorousFlowPointSourceFromPostprocessor
mass_flux = flux
point = '0.5 0 0'
variable = sat1
[]
[]
The remainder of the input file is standard, with the inclusion of the hysteretic capillary pressure:
[hys_order_material]
type = PorousFlowHysteresisOrder
[]
[pc_calculator]
type = PorousFlow2PhaseHysPS
alpha_d = 10.0
alpha_w = 7.0
n_d = 1.5
n_w = 1.9
S_l_min = 0.1
S_lr = 0.2
S_gr_max = 0.3
Pc_max = 12.0
high_ratio = 0.9
low_extension_type = quadratic
high_extension_type = power
phase0_porepressure = pp0
phase1_saturation = sat1
[]
[]
The result is Figure 14.
Figure 14: The result of a two-phase simulation using a PS formulation in which an external pump adds, removes and then adds gas to a porous material in order to observe the hysteretic capillary pressure.
Hysteretic capillary-pressure implementation remarks
As mentioned above, Figure 10 reveals a subtle feature of the implementation of hysteresis in PorousFlow. There is no asterisk around . Instead it appears on the primary drying curve at around . As the saturation is reduced along the drying curve, and then increased again, PorousFlow follows the primary drying curve for the first time-step after the turning point. At subsequent steps, it follows the correct curve.
This is because the computation of hysteresis order lags one timestep behind the MOOSE simulation. This is to ensure reasonable convergence behavior, as mentioned in Doughty (2007) and Doughty (2008).
The lagging behaviour has the unfortunate side-effect that simulations involving hysteretic capillary pressure do not conserve fluid mass. The remainder of PorousFlow conserves fluid mass exactly, at all times. However, when dealing with hysteretic capillary pressures, a small amount of mass is lost or gained whenever a turning point is encountered. This should have limited impact upon models if saturations do not change excessively during a single time-step.
Input file syntax for hysteretic relative permeability
To include hysteretic relative permeability in an existing (non-hysteretic) input file, the following changes need to be made
A PorousFlowHysteresisOrder Material needs to be included.
Materials are needed to compute the relative permeability. The existing Materials, such as PorousFlowRelativePermeabilityVG need to be replaced with PorousFlowHystereticRelativePermeabilityLiquid, and, for 2-phase simulations PorousFlowHystereticRelativePermeabilityGas.
A 1-phase example is:
[hys_order_material]
type = PorousFlowHysteresisOrder
[]
[relperm_material]
type = PorousFlowHystereticRelativePermeabilityLiquid
phase = 0
S_lr = 0.1
S_gr_max = 0.2
m = 0.9
liquid_modification_range = 0.9
[]
[]
commentnote
Models need not contain both hysteretic capillary pressures and hysteretic relative permeabilities: the hysteresis may only appear in the capillary pressure, or the relative permeabilities.
Single-phase examples
A simulation that simply removes and adds water to a system to observe the hysteretic relative permeability is explored in this section. The water flux is controlled by the Postprocessor
[flux]
type = FunctionValuePostprocessor
function = 'if(t <= 5, -10, 10)'
[]
and the DiracKernel
[DiracKernels]
[pump]
type = PorousFlowPointSourceFromPostprocessor
mass_flux = flux
point = '0.5 0 0'
variable = pp
[]
[]
The remainder of the input file is standard, with the inclusion of the hysteretic capillary pressure:
[hys_order_material]
type = PorousFlowHysteresisOrder
[]
[relperm_material]
type = PorousFlowHystereticRelativePermeabilityLiquid
phase = 0
S_lr = 0.1
S_gr_max = 0.2
m = 0.9
liquid_modification_range = 0.9
[]
[]
The result is Figure 15. By altering the flux, the system may be dried, re-wet, dried and re-wet again, to generate results such as Figure 16
Figure 15: The result of a single-phase simulation in which an external pump removes and adds water to a porous material in order to observe the hysteretic relative permeability.
Figure 16: The result of a single-phase simulation in which an external pump removes water until , then adds water until , then removes water until , then adds water until full saturation is reached.
Two-phase examples
A simulation that simply adds and removes gas from a system to observe the hysteretic relative permeability is explored in this section. The water flux is controlled by the Postprocessor
[flux]
type = FunctionValuePostprocessor
function = 'if(t <= 9, 10, -10)'
[]
and the DiracKernel
[DiracKernels]
[pump]
type = PorousFlowPointSourceFromPostprocessor
mass_flux = flux
point = '0.5 0 0'
variable = sat1
[]
[]
The remainder of the input file is standard, with the inclusion of the hysteretic relative permeabilities:
[hys_order_material]
type = PorousFlowHysteresisOrder
[]
[relperm_liquid]
type = PorousFlowHystereticRelativePermeabilityLiquid
phase = 0
S_lr = 0.1
S_gr_max = 0.2
m = 0.9
liquid_modification_range = 0.9
[]
[relperm_gas]
type = PorousFlowHystereticRelativePermeabilityGas
phase = 1
S_lr = 0.1
S_gr_max = 0.2
m = 0.9
gamma = 0.33
k_rg_max = 0.8
gas_low_extension_type = linear_like
[]
[]
The result is Figure 17. By altering the flux, k_rg_max and gas_low_extension_type, the impact of various extensions may be explored, as shown in Figure 18, Figure 19 and Figure 20.
Figure 17: The result of a two-phase simulation in which an external pump adds and removes gas from a porous material in order to observe the hysteretic relative permeability.
Figure 18: The result of a two-phase simulation in which an external pump adds and removes gas from a porous material in order to observe the hysteretic relative permeability. A cubic extension is used for the gas relative permeability.
Figure 19: The result of a two-phase simulation in which an external pump adds and removes gas from a porous material in order to observe the hysteretic relative permeability. A linear-like extension is used for the gas relative permeability.
Figure 20: The result of a two-phase simulation in which an external pump adds and removes gas from a porous material in order to observe the hysteretic relative permeability. No extension is used for the gas relative permeability, since is unity.
References
- C. Doughty.
Modeling geologic storage of carbon dioxide: comparison of non-hysteretic and hysteretic characteristic curves.
Energy Conservation and Management, 48:1768–1781, 2007.
doi:10.1016/j.enconman.2007.01.022.[BibTeX]
- C. Doughty.
User's guide for hysteretic capillary pressure and relative permeability functions in iTOUGH2, Rev3, Draft, March 2008.
Technical Report, Earth Sciences Division, Lawrence Berkeley National Laboratory, 2008.[BibTeX]
- S. Finsterle, T. O. Sonnenborg, and B. Faybishenko.
Inverse modelling of a multistep outflow experiment for determining hysteretic hydraulic properties.
In K Pruess (ed.), Proceedings of the TOUGH Workshop '98, Lawrence Berkeley National Laboratory Report LBNL-41995, pp. 250–256. 1998.[BibTeX]
- D. B. Jaynes.
Comparison of soil-water hysteresis models.
Journal of Hydrology, 75:287–299, 1994.[BibTeX]
- A. Niemi and G. S. Bodvarsson.
Preliminary capillary hysteresis simulations in fractured rocks, Yucca Mountain, Nevada.
Journal of Contaminant Hydrology, 3:277–291, 1988.[BibTeX]