minires.fluids

 1"""The two-phase fluid, held at `minires.ResSim.fluid`.
 2
 3`Fluid` holds the phase viscosities and the parameters of the Corey relative
 4permeability curves, and computes what the model needs of them: the mobilities
 5(`RelPerm`, Listing 6 of the reference paper) and the fractional flow
 6(`fractional_flow`), each with its derivative (`dRelPerm`, `dfractional_flow`).
 7Curves of another shape (e.g. tabulated) are a subclass overriding `RelPerm` and
 8`dRelPerm`.
 9"""
10
11from dataclasses import dataclass
12
13import numpy as np
14
15from minires._repr import AlignedRepr
16
17
18@dataclass
19class Fluid(AlignedRepr):
20    """A two-phase (water/oil) fluid: viscosities and Corey relative permeabilities.
21
22    The relative permeabilities are Corey (power-law) curves of the normalized
23    saturation $ S $ (`rescale_sat`),
24
25    $$ k_{rw} = k_{rw}^0 \\, S^{n_w} \\,, \\qquad k_{ro} = k_{ro}^0 \\, (1 - S)^{n_o} \\,,
26       \\qquad S = \\frac{s - S_{wc}}{1 - S_{wc} - S_{or}} \\,, $$
27
28    with $ S $ clipped to $[0, 1]$, so that a phase below its residual saturation
29    is immobile (rather than mobile with the wrong sign, as an odd power would
30    make it). The defaults give the quadratic curves of the reference paper
31    (Listing 6), $ S^2 $ and $ (1 - S)^2 $, and unit viscosities.
32
33    >>> fluid = Fluid(vo=2, swc=.2, sor=.2)
34    >>> Mw, Mo = fluid.RelPerm(np.array([.1, .2, .5, .8, .9]))
35    >>> Mw.round(4), Mo.round(4)
36    (array([0.  , 0.  , 0.25, 1.  , 1.  ]), array([0.5  , 0.5  , 0.125, 0.   , 0.   ]))
37    """
38
39    __repr__ = AlignedRepr.__repr__
40
41    vw: float = 1.0
42    """Viscosity for water."""
43    vo: float = 1.0
44    """Viscosity for oil."""
45    swc: float = 0.0
46    """Irreducible saturation, water: $ k_{rw} = 0 $ below it."""
47    sor: float = 0.0
48    """Irreducible saturation, oil: $ k_{ro} = 0 $ above $ s = 1 - S_{or} $."""
49    nw: float = 2.0
50    """Corey exponent of the water relative permeability."""
51    no: float = 2.0
52    """Corey exponent of the oil relative permeability."""
53    krw0: float = 1.0
54    """End-point (maximal) water relative permeability, at $ s = 1 - S_{or} $."""
55    kro0: float = 1.0
56    """End-point (maximal) oil relative permeability, at $ s = S_{wc} $."""
57
58    def rescale_sat(self, s: np.ndarray) -> np.ndarray:
59        """The normalized saturation $ S $ (unclipped). Ref paper, p. 32."""
60        return (s - self.swc) / (1 - self.swc - self.sor)
61
62    # RelPerm() -- listing 6
63    def RelPerm(self, s: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
64        """Rel. permeabilities of water and oil, as mobilities (perm/viscosity)."""
65        S = np.clip(self.rescale_sat(s), 0, 1)
66        Mw = self.krw0 * S**self.nw / self.vw
67        Mo = self.kro0 * (1 - S) ** self.no / self.vo
68        return Mw, Mo
69
70    def dRelPerm(self, s: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
71        """Derivatives of `RelPerm` wrt `s` (zero where the curves are clipped)."""
72        S = self.rescale_sat(s)
73        inside = (0 <= S) & (S <= 1)  # one-sided at the ends, as the reference code
74        S = np.clip(S, 0, 1)
75        w = 1 - self.swc - self.sor
76        dMw = np.where(inside, self.krw0 * self.nw * S ** (self.nw - 1) / w, 0) / self.vw
77        dMo = -np.where(inside, self.kro0 * self.no * (1 - S) ** (self.no - 1) / w, 0) / self.vo
78        return dMw, dMo
79
80    def fractional_flow(self, s: np.ndarray) -> np.ndarray:
81        """The water fractional flow, $ f_w = λ_w / λ_t $."""
82        Mw, Mo = self.RelPerm(s)
83        return Mw / (Mw + Mo)
84
85    def dfractional_flow(self, s: np.ndarray) -> np.ndarray:
86        """Derivative of `fractional_flow` wrt `s`.
87
88        Separate from `fractional_flow` (rather than a second output of it) because
89        the explicit transport scheme wants $ f_w $ alone, in its innermost loop,
90        where the derivative would cost as much again as everything else.
91        """
92        Mw, Mo = self.RelPerm(s)
93        dMw, dMo = self.dRelPerm(s)
94        return (dMw * Mo - Mw * dMo) / (Mw + Mo) ** 2