Sie sind hier: ICP » R. Hilfer » Publikationen

2 Formulation of the Model

2.1 Balance Laws

[page 509, §1]   
[509.1.1] Consider a onedimensional, homogeneous, isotropic and incompressible porous column filled with two immiscible Newtonian fluids. [509.1.2] In a one dimensional model transversal variations and column thickness are neglected. [509.1.3] Let {S_{{\mathbb{W}}}}={S_{{\mathbb{W}}}}(x,t) denote the saturation of wetting fluid (called water), and {S_{{\mathbb{O}}}}={S_{{\mathbb{O}}}}(x,t) the saturation of nonwetting fluid (called oil). [509.1.4] Here time is t\geq 0, and x\in[0,L] is the position in a column of length L. [509.1.5] Each of the two fluid phases is considered to consist of a continuous, mobile, percolating subphase, and of a discontinuous, isolated, trapped or nonpercolating subphase as discussed in [23, 28, 25, 24, 26]. [509.1.6] Following the notation of [26] the percolating phase of water is [page 510, §0]    indexed by i=1 and its nonpercolating phase is indexed by i=2. [510.0.1] The water saturation is then {S_{{\mathbb{W}}}}=S_{1}+S_{2}. [510.0.2] The percolating oil phase is indexed as i=3 and its nonpercolating phase by i=4. [510.0.3] The oil saturation is {S_{{\mathbb{O}}}}=S_{3}+S_{4}. [510.0.4] The volume fraction \phi _{i} of phase i is defined as \phi _{i}=\phi S_{i} where the volume fraction of the pore space, also called porosity, is \phi. [510.0.5] The volume fraction of the solid matrix is denoted as \phi _{5}=1-\phi. [510.0.6] Volume conservation requires

\displaystyle\phi _{1}+\phi _{2}+\phi _{3}+\phi _{4}+\phi _{5} \displaystyle=1 (1a)
\displaystyle S_{1}+S_{2}+S_{3}+S_{4} \displaystyle=1 (1b)

to hold.

[510.1.1] The mass balance of fluid phase i can be expressed in differential form as

\frac{\partial(\phi _{i}\varrho _{i})}{\partial t}+\frac{\partial(\phi _{i}\varrho _{i}{v}_{i})}{\partial x}=M_{i} (2)

where \varrho _{i}(x,t),\phi _{i}(x,t),{v}_{i}(x,t) are mass density, volume fraction and velocity of phase i as functions of position x\in\mathbb{S}=[0,L]\subset\mathbb{R} and time t\in\mathbb{R}_{+}. [510.1.2] M_{i} is the mass transfer rate from all other phases into phase i.

[510.2.1] The momentum balance is written as (i=1,2,3,4)

\phi _{i}\varrho _{i}\frac{{\rm D}^{i}}{{\rm D}t}{v}_{i}-\phi _{i}\frac{\partial\Sigma _{i}}{\partial x}-\phi _{i}F_{i}=m_{i}-{v}_{i}M_{i} (3)

where \Sigma _{i} is the stress tensor in the ith phase, F_{i} is the body force per unit volume acting on the ith phase, m_{i} is the momentum transfer into phase i from all the other phases, and {\rm D}^{i}/{\rm D}t=\partial/\partial t+{v}_{i}\partial/\partial x denotes the material derivative for phase i. [510.2.2] The stress tensor for the nonpercolating phases is defined as the momentum flux across surfaces in the threedimensional continuum (see [29] for more discussion).

2.2 Constitutive Assumptions

[510.3.1] For a macroscopically homogeneous porous medium

\phi(x)=\phi={\rm const} (4)

is assumed. [510.3.2] Incompressible fluids are assumed so that their densities

\displaystyle\varrho _{1}(x,t) \displaystyle=\varrho _{{\mathbb{W}}} (5a)
\displaystyle\varrho _{2}(x,t) \displaystyle=\varrho _{{\mathbb{W}}} (5b)
\displaystyle\varrho _{3}(x,t) \displaystyle=\varrho _{{\mathbb{O}}} (5c)
\displaystyle\varrho _{4}(x,t) \displaystyle=\varrho _{{\mathbb{O}}} (5d)

are independent of x,t.

[510.4.1] The percolating and the nonpercolating phases are able to exchange mass through breakup and coalescence of droplets, ganglia and clusters. [510.4.2] The mass transfer rates must depend on rates of saturation change. [510.4.3] They are here assumed to be

\displaystyle M_{1} \displaystyle=-M_{2}={\eta _{{2}}}\phi\varrho _{{\mathbb{W}}}\left(\frac{S_{2}-{{S_{2}}^{*}}}{{{S_{{\mathbb{W}}}}^{*}}-{S_{{\mathbb{W}}}}}\right)\frac{\partial{S_{{\mathbb{W}}}}}{\partial t} (6a)
\displaystyle M_{3} \displaystyle=-M_{4}={\eta _{{4}}}\phi\varrho _{{\mathbb{O}}}\left(\frac{S_{4}-{{S_{4}}^{*}}}{{{S_{{\mathbb{O}}}}^{*}}-{S_{{\mathbb{O}}}}}\right)\frac{\partial{S_{{\mathbb{O}}}}}{\partial t} (6b)

[page 511, §0]    where {\eta _{{2}}},{\eta _{{4}}} are constants. [511.0.1] The parameters {{S_{{\mathbb{W}}}}^{*}},{{S_{{\mathbb{O}}}}^{*}}, {{S_{2}}^{*}},{{S_{4}}^{*}} are defined by

\displaystyle{{S_{{\mathbb{W}}}}^{*}} \displaystyle=(1-S_{{\mathbb{O}\,\rm r}})\Theta\left(\partial _{t}{S_{{\mathbb{W}}}}\right)+S_{{\mathbb{W}\,\rm i}}\left[1-\Theta\left(\partial _{t}{S_{{\mathbb{W}}}}\right)\right] (7a)
\displaystyle{{S_{{\mathbb{O}}}}^{*}} \displaystyle=1-{{S_{{\mathbb{W}}}}^{*}} (7b)
\displaystyle{{S_{2}}^{*}} \displaystyle=S_{{\mathbb{W}\,\rm i}}\left[1-\Theta\left(\partial _{t}{S_{{\mathbb{W}}}}\right)\right] (7c)
\displaystyle{{S_{4}}^{*}} \displaystyle=S_{{\mathbb{O}\,\rm r}}\Theta\left(\partial _{t}{S_{{\mathbb{W}}}}\right) (7d)

where S_{{\mathbb{W}\,\rm i}},S_{{\mathbb{O}\,\rm r}} are limiting saturations for S_{2},S_{4}. [511.0.2] In eq. (7) the shorthand \partial _{t}=\partial/\partial t is used and

\Theta(x)=\begin{cases}1&\text{for~}x>0,\\ 0&\text{for~}x\leq 0,\end{cases} (8)

denotes the Heaviside unit step function. [511.0.3] Equation (7) follows from the form used in [26] for small rates of saturation change. [511.0.4] The mass exchange depends on the sign of \partial _{t}{S_{{\mathbb{W}}}}. [511.0.5] The sign determines the type of process. [511.0.6] It can switch locally between drainage and imbibition. [511.0.7] This results in hysteresis. [511.0.8] The structure of the mass exchange term was chosen with hindsight such that theoretical results obtained in the residual decoupling approximation agree with experimental measurements of capillary pressure. [511.0.9] The mass exchange terms have recently been further generalized to reproduce not only capillary pressure, but also experimental capillary desaturation curves [27].

[511.1.1] Turning to the momentum balance, note first that the inertial term will not be neglected in this paper. [511.1.2] The stress tensor for the four phases are specified as

\displaystyle\Sigma _{1} \displaystyle=-P_{1} (9a)
\displaystyle\Sigma _{2} \displaystyle=-P_{3}+\gamma P^{*}_{2}S_{2}^{{\gamma-1}} (9b)
\displaystyle\Sigma _{3} \displaystyle=-P_{3} (9c)
\displaystyle\Sigma _{4} \displaystyle=-P_{1}+\delta P^{*}_{4}S_{4}^{{\delta-1}} (9d)

where P_{1} and P_{3} are the fluid pressures in the percolating phases. [511.1.3] The constants P^{*}_{2},P^{*}_{4} and exponents \gamma,\delta are constitutive parameters. [511.1.4] The assumptions for the nonpercolating phases reflect their modified pressure. [511.1.5] This phenomenon seems to have been observed in experiment [2, Fig. 2, p. 233]. a (This is a footnote:) aFig. 2 in [2] shows that the pressure measured by the pore pressure transducers PPT3 and PPT4 rebounds after the end of the infiltration, i.e. when the DNAPL has passed and the transducers measure the pressure of water and residual PCE [511.1.6] In applications, the parameters P^{*}_{2},P^{*}_{4} and \gamma,\delta are determined by measuring capillary pressure curves (see below). [511.1.7] The body forces are assumed to be given by gravity and capillarity. [511.1.8] They are specified as

\displaystyle F_{1} \displaystyle=\varrho _{1}g\sin\vartheta (10a)
\displaystyle F_{2} \displaystyle=\varrho _{2}g\sin\vartheta+\Pi _{{\rm a}}\frac{\partial S_{1}^{{-\alpha}}}{\partial x} (10b)
\displaystyle F_{3} \displaystyle=\varrho _{3}g\sin\vartheta (10c)
\displaystyle F_{4} \displaystyle=\varrho _{4}g\sin\vartheta+\Pi _{{\rm b}}\frac{\partial S_{3}^{{-\beta}}}{\partial x} (10d)

with constitutive constants \Pi _{{\rm a}},\Pi _{{\rm b}} and exponents \alpha,\beta>0. [511.1.9] The angle 0\leq\vartheta\leq\pi/2 is the angle between the direction of the column and the direction of gravity with \vartheta=\pi/2 [page 512, §0]    corresponding to alignment. [512.0.1] In applications the parameters \Pi _{{\rm a}},\Pi _{{\rm b}} and \alpha,\beta are determined by measuring capillary pressure curves. [512.0.2] The capillary body forces in eqs. (10) reflect the wetting properties of the medium. [512.0.3] They oppose gravity and reduce buoyancy driven flows of the disconnected phases. [512.0.4] This is illustrated in the figures below.

[512.1.1] Finally, the momentum transfer terms are assumed to be given by linear viscous drag characterized by constitutive resistance coefficients R_{{ij}} through the equations

\displaystyle m_{1} \displaystyle=R_{{13}}({v}_{3}-{v}_{1})+R_{{14}}({v}_{4}-{v}_{1})-R_{{15}}{v}_{1} (11a)
\displaystyle m_{2} \displaystyle=R_{{23}}({v}_{3}-{v}_{2})+R_{{24}}({v}_{4}-{v}_{2})-R_{{25}}{v}_{2} (11b)
\displaystyle m_{3} \displaystyle=R_{{31}}({v}_{1}-{v}_{3})+R_{{32}}({v}_{2}-{v}_{3})-R_{{35}}{v}_{3} (11c)
\displaystyle m_{4} \displaystyle=R_{{41}}({v}_{1}-{v}_{4})+R_{{42}}({v}_{2}-{v}_{4})-R_{{45}}{v}_{4} (11d)

where R_{{12}}=0 and R_{{34}}=0 was used because there is no common interface and hence no direct viscous interaction between these phase pairs. [512.1.2] Remember that the index i=5 represents the rock matrix. [512.1.3] For more details on these constitutive assumptions the reader is referred to the original papers [24, 25, 26].

[512.2.1] The balance laws (1b), (2) and (3) combined with the constitutive assumptions given above provide 9 equations for the 10 unknowns S_{i},{v}_{i},P_{1},P_{3} with i=1,2,3,4. [512.2.2] To close the system of equations the conditions {v}_{2}=0 or {v}_{4}=0 could be used. [512.2.3] These conditions apply when the nonpercolating phases are immobile as it is often observed in experiment. [512.2.4] It turns out, however, that there exists a less restrictive and, in our opinion, more natural selfconsistent closure.

[512.3.1] The selfconsistent closure condition used in this paper follows naturally from many limiting cases. [512.3.2] One such limit is the residual decoupling approximation close to hydrostatic equilibrium described in detail in [26, Section 5., p. 216ff]. [512.3.3] A second, more general limiting case is the limit of vanishing velocities, i.e. {v}_{i}\to 0 for i=1,2,3,4. [512.3.4] Here we formulate the selfconsistent closure condition in its most general form as

\displaystyle 0= \displaystyle\left(\frac{R_{{13}}}{\phi _{1}}+\frac{R_{{14}}}{\phi _{1}}+\frac{R_{{15}}}{\phi _{1}}+\frac{R_{{31}}}{\phi _{3}}-\frac{R_{{41}}}{\phi _{4}}+\frac{M_{1}}{\phi _{1}}\right){v}_{1}+\varrho _{1}\frac{{\rm D}^{1}}{{\rm D}t}{v}_{1}
\displaystyle+\left(-\frac{R_{{23}}}{\phi _{2}}-\frac{R_{{24}}}{\phi _{2}}-\frac{R_{{25}}}{\phi _{2}}+\frac{R_{{32}}}{\phi _{3}}-\frac{R_{{42}}}{\phi _{4}}+\frac{M_{1}}{\phi _{2}}\right){v}_{2}-\varrho _{2}\frac{{\rm D}^{2}}{{\rm D}t}{v}_{2}
\displaystyle+\left(-\frac{R_{{13}}}{\phi _{1}}+\frac{R_{{23}}}{\phi _{2}}-\frac{R_{{31}}}{\phi _{3}}-\frac{R_{{32}}}{\phi _{3}}-\frac{R_{{35}}}{\phi _{3}}-\frac{M_{3}}{\phi _{3}}\right){v}_{3}-\varrho _{3}\frac{{\rm D}^{3}}{{\rm D}t}{v}_{3}
\displaystyle+\left(-\frac{R_{{14}}}{\phi _{1}}+\frac{R_{{24}}}{\phi _{2}}+\frac{R_{{41}}}{\phi _{4}}+\frac{R_{{42}}}{\phi _{4}}+\frac{R_{{45}}}{\phi _{4}}-\frac{M_{3}}{\phi _{4}}\right){v}_{4}+\varrho _{4}\frac{{\rm D}^{4}}{{\rm D}t}{v}_{4}. (12)

[page 513, §1]    [513.2.1] This condition follows selfconsistently from the constitutive theory. [513.2.2] It expresses the experimental observation that the pressure difference P_{3}-P_{1} depends more strongly on saturations than on velocities, and that it remains nonzero even for vanishing velocities. [513.2.3] Adding eqs. (3) for i=2 and i=3 and subtracting eqs. (3) with i=1 and i=4 from the result gives

\frac{\partial P_{3}}{\partial x}=\frac{\partial P_{1}}{\partial x}+\frac{\partial}{2\partial x}\left[\Pi _{{\rm a}}S_{1}^{{-\alpha}}-\Pi _{{\rm b}}S_{3}^{{-\beta}}+\gamma P^{*}_{2}S_{2}^{{\gamma-1}}-\delta P^{*}_{4}S_{4}^{{\delta-1}}\right]. (13)

[513.2.4] In this form the selfconsistent closure has been used in the numerical calculation below.

[513.3.1] The mathematical model defined above was introduced and discussed in [25, 24, 26]. [513.3.2] It was recently extended to include surface tension [27]. [513.3.3] Notation and model formulation in this paper follow [26]. [513.3.4] Approximations and analytical solutions for some special cases were given in [25, 24, 26]. [513.3.5] Here the system of equations will be solved by numerical methods. [513.3.6] To this end initial and boundary conditions are needed and will be discussed next.