15th Annual Conference of Applied Mathematics, Univ. of Central Oklahoma, Electronic Journal of Differential Equations, Conference 02, 1999, pp. 29–46. ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu (login: ftp)

Degenerate two-phase incompressible flow problems III: Perturbation analysis and

numerical experiments ^{∗}

Zhangxin Chen & Natalia L. Khlopina

Abstract

This is the third paper of a three-part series where we develop and analyze a finite element approximation for a degenerate elliptic-parabolic partial differential system which describes the flow of two incompressible, immiscible fluids in porous media. The approximation uses a mixed finite element method for the pressure equation and a Galerkin finite element method for the saturation equation. It is based on a regularization of the saturation equation. In the first paper [15] we analyzed the regularized differential system and presented numerical results. In the second paper

[16] we obtained error estimates. In the present paper we describe a perturbation analysis for the saturation equation and numerical experiments for complementing this analysis.

1 Introduction

The flow of two incompressible, immiscible fluids in a porous medium Ω ⊂ ℝ^{d}, d ≤ 3 [2, 25] is given by

\begin{matrix}ϕ∂_{t}s − ∇ · (κλ_{w}(s)(∇p_{w} + γ_{w})) = q_{w}, \\ −ϕ∂_{t}s − ∇ · (κλ_{o}(s)(∇p_{o} + γ_{o})) = q_{o}, \\ p_{c}(s) = p_{o} − p_{w}, \end{matrix}  (1.1)

where w indicates a wetting phase (e.g., water), o denotes a nonwetting phase (e.g., oil), ϕ and κ are the porosity and absolute permeability of the porous system, s is the (reduced) saturation of the wetting phase, p_{α}, λ_{α}, γ_{α}, and q_{α} are, respectively, the pressure, mobility (i.e., the relative permeability over the viscosity), gravity-density vector, and external volumetric flow rate of the α-phase (α = w, o), and p_{c} is the capillary pressure function.

∗1991 Mathematics Subject Classifications: 35K60, 35K65, 76S05, 76T05.

Key words and phrases: porous medium, degenerate elliptic-parabolic system, perturbation method, finite element, regularization, two-phase flow.

\begin{matrix}c \\ ◯ \end{matrix}1999 Southwest Texas State University and University of North Texas.

Partly supported by National Science Foundation grants DMS-9626179, DMS-9972147, and INT-9901498 and a gift grant from the Mobil Oil Company, Dallas, Texas.

Published November 24, 1999.

To separate the saturation equation from the pressure equation, we define the total mobility

λ(x, s) = λ_{w} + λ_{o}.

Also, following [1, 6] we define the global pressure

p = p_{o} − ∫_{0}^{s} (\frac{λ_{w}}{λ} \frac{∂p_{c}}{∂s} ) (x, ξ)dξ,  (1.2)

and following [8] the complementary pressure

θ = D(s) = − ∫_{0}^{s} (κ\frac{λ_{w}λ_{o}}{λ} \frac{∂p_{c}}{∂s} ) (x, ξ)dξ.  (1.3)

Then, by (1.2), (1.3), and some manipulations, it follows from (1.1) that [8, 14]

\begin{matrix}−∇ · \{κ(λ(s)∇p + γ′_{1}(s))\} = q ≡ q_{w} + q_{o}, \\ ϕ∂_{t}s − ∇ · \{∇θ + κ(λ_{w}(s)∇p + γ′_{2}(s))\} = q_{w}, \end{matrix}  (1.4)

where

γ′_{1} = −λ_{w}∇_{x}p_{c} + λ ∫_{0}^{s} ∇_{x} (\frac{λ_{w}}{λ} \frac{∂p_{c}}{∂s} ) (x, ξ)dξ + λ_{w}γ_{w} + λ_{o}γ_{o},
γ′_{2} = −λ_{w}∇_{x}p_{c} + λ_{w} ∫_{0}^{s} ∇_{x} (\frac{λ_{w}}{λ} \frac{∂p_{c}}{∂s} ) (x, ξ)dξ + λ_{w}γ_{w}
+ ∫_{0}^{s} ∇_{x} (\frac{λ_{w}λ_{o}}{λ} \frac{∂p_{c}}{∂s} ) (x, ξ)dξ.

In (1.4), s is related to θ through (1.3):

s = S(θ),  (1.5)

where S(x, θ) is the inverse of D(s) for 0 ≤ θ ≤ θ^{∗}(x) with

θ^{∗}(x) = − ∫_{0}^{1} (κ\frac{λ_{w}λ_{o}}{λ} \frac{∂p_{c}}{∂s} ) (x, ξ)dξ.

The pressure equation is given by the first equation of (1.4), while the saturation equation is described by the second equation. They determine the main unknowns p, s, and θ. The model is completed by specifying boundary and initial conditions.

For simplicity, in this paper we only consider the Neumann boundary conditions

\begin{matrix}−κ(λ(s)∇p + γ′_{1}(s)) · ν = φ_{1}(x, t), & (x, t) ∈ Γ × J, \\ −\{∇θ + κ(λ_{w}(s)∇p + γ′_{2}(s))\} · ν = φ_{2}(x, t), & (x, t) ∈ Γ × J, \end{matrix}  (1.6)

where φ_{1} and φ_{2} are given functions, J = (0, T] (T > 0), Γ is the boundary of Ω, and ν is the outer unit normal to Γ. The boundary conditions in (1.6) come from those imposed for the phase quantities via the transformations (1.2) and (1.3) [12]. For other types of boundary conditions, refer to [14]. The initial condition is given by

s(x, 0) = s_{0}(x), x ∈ Ω.  (1.7)

In recent years, the interest in the numerical simulation of two-phase fluid flow in porous media has been rising rapidly (see [18] and the bibliographies therein). In conjunction, there has been intensive research into the error analysis of numerical methods used in the simulation (see the extensive references in [12]). However, in most previous works the error analysis has been carried out under the unrealistic assumption that the capillary diffusion coefficient is uniformly positive. Namely, it has been assumed that −(λ_{w}λ_{o}/λ)(∂p_{c}/∂s) is uniformly positive. The case where this diffusion coefficient can be zero has been treated in [11, 12, 19, 20, 28]. But, in these papers only simplied saturation equations have been analyzed; i.e., the second equation in (1.4) with p (or u the total velocity; see (3.1a) later) given and q_{w} = γ′_{2} ≡ 0 has been considered [11, 12, 19, 20, 28]. More recently, in [14] the fully coupled system (1.4) has been analyzed for the finite element approximation used here. It has been shown [14] that when this approximation directly solves the degenerate system, optimal error estimates behave like O(h), where h is the discretization mesh size. The main purposes of this series of three papers are to analyze regularized versions of this fully coupled system, improve the error estimates in [14], describe a perturbation analysis, and present numerical experiments.

The error analysis presented in [16] was based on a regularization of the saturation equation. The diffusion coefficient of this equation was perturbed to obtain a nondegenerate problem with smooth solutions. The regularized solutions were shown to converge to the original solution as the perturbation parameter goes to zero with specific convergence rates given. Then a finite element approximation was used to solve the regularized solutions of the differential system. This approximation, which follows [14], combined a mixed finite element method for the pressure equation and a Galerkin finite element method for the regularized saturation equation. Then, for this approximation we proved that the norm of error estimates depended on the severity of the degeneracy in diffusivity. In particular, for the degeneracy under consideration error estimates we obtained can behave like O(h^{1+ϵ}), where 0 < ϵ < 1. This result thus improved that in [14] in some cases.

In the first paper [15] we analyzed the regularized differential system, described the finite element approximation, and presented numerical results. In the second paper [16] we obtained error estimates. Both semidiscrete (continuous in time) and fully discrete approximations were analyzed. In the present paper we describe a perturbation analysis for the saturation equation and carry out numerical experiments for complementing this analysis. We remark that since the differential system for the single-phase, miscible displacement of one incompressible fluid by another in porous media resembles that for the twophase incompressible flow studied here [13], the analysis presented in this paper extends to the miscible displacement problem. Also, due to its convectiondominated feature, more efficient approximate procedures should be used to solve the saturation equation. Characteristic-based finite element methods will be considered in forthcoming papers.

The rest of the paper is organized as follows. After preliminaries in the next section we consider a regularization of the differential system (1.4) in the third section and a perturbation method for the saturation equation in the forth section. Then we describe a finite element approximation for this regularized system in the fifth section. Numerical experiments are given in the final section.

2 Preliminaries

In this section we present preliminary results used in later sections. In particular, we collect some results for the Poisson solution operator E, which is defined below, and for the regularized diffusion coefficient.

For g ∈ H^{−1}(Ω) (the dual to H^{1}(Ω)), set

g_{Ω} = (g, 1)

in the sense of duality between H^{−1}(Ω) and H^{1}(Ω). When g is Lebesgue integrable on Ω, we have

g_{Ω} = ∫_{Ω} gdx.

For g ∈ H^{−1}(Ω), consider the Neumann boundary value problem

−Δω = g − g_{Ω}  in Ω,
∇ω · ν = 0  on Γ,
ω_{Ω} = g_{Ω}.

(2.1)

The Poisson solution operator E : H^{−1}(Ω) → H^{1}(Ω) is defined by

E(g) = ω.

It follows from (2.1) that

(∇ω, ∇v) = (g, v) − g_{Ω}v_{Ω}, v ∈ H^{1}(Ω);

i.e.,

(g, v) = (∇(E(g)), ∇v) + g_{Ω}v_{Ω}.

Especially, take v = E(g) to see that

(g, E(g)) = ∥∇(E(g))∥_{L^{2}(Ω)}^{2} + (E(g)_{Ω})^{2}.  (2.2)

It can be easily checked [19] that the operator E is linear, symmetric, and positive definite. Furthermore, for v ∈ H^{1}(Ω) note that the norm

(∥∇v∥_{L^{2}(Ω)}^{2} + (v_{Ω})^{2})^{1/2}

is equivalent to the usual norm on H^{1}(Ω). Thus, by (2.2), we can define the norm on H^{−1}(Ω)

∥g∥_{H^{−1}(Ω)} = (g, E(g))^{1/2},

which is equivalent to the usual norm on H^{−1}(Ω).

As mentioned in the introduction we improve the error estimate in [14] by studying a regularized version of the saturation equation. For this, we define the capillary diffusion coefficient

d(s) = −κ\frac{λ_{w}λ_{o}}{λ} \frac{∂p_{c}}{∂s} ,

and assume that it satisfies

d(s) ≥ \{ \begin{matrix}c_{1}|s|^{μ_{1}}, & 0 ≤ s ≤ β_{1}, \\ c_{2}, & β_{1} ≤ s ≤ β_{2}, \\ c_{3}|1 − s|^{μ_{2}}, & β_{2} ≤ s ≤ 1, \end{matrix}  (2.3)

where c_{i} (i = 1, 2, 3) are positive constants, and μ_{i} and β_{i} (i = 1, 2) satisfy

0 < β_{1} < 1/2 < β_{2} < 1, 0 < μ_{1}, μ_{2} ≤ 2.

Set

μ = max{μ_{1}, μ_{2}},

and

γ = \frac{2 + μ}{1 + μ}.

Note that γ is the conjugate to 2 + μ. As in (1.3), we also set

θ = D(s) = ∫_{0}^{s} d(ξ)dξ.

Remark that (2.3) assumes the nature of degeneracy in the coefficient d(s) near zero and one. The following lemma can be found in [19].

Lemma 2.1. Let d satisfy (2.3). Then there exists a positive constant C, depending only on the parameters in (2.3), such that

C(s_{2} − s_{1})^{1+μ} ≤ D(s_{2}) − D(s_{1}),  0 ≤ s_{1} ≤ s_{2} ≤ 1.  (2.4)

3 Regularization

For the convenience of the later analysis, we rewrite (1.4) as follows:

\begin{matrix}∇ · u = q, u = −κ(λ(s)∇p + γ′_{1}(s)) & in Ω_{T} , \\ ϕ∂_{t}s − ∇ · \{∇D(s) − f_{w}(s)u + γ_{2}(s)\} = q_{w} & in Ω_{T} , \end{matrix}  (3.1a)

where Ω_{T} = Ω × J and

f_{w}(s) = λ_{w}(s)/λ(s), γ_{2}(s) = κ{γ′_{2}(s) − f_{w}(s)γ′_{1}(s)}.

The boundary and initial conditions become

\begin{matrix}u · ν = φ_{1}(x, t), & (x, t) ∈ Γ × J, \\ −(∇D(s) − f_{w}(s)u + γ_{2}(s)) · ν = φ_{2}(x, t), & (x, t) ∈ Γ × J, \\ s(x, 0) = s_{0}(x), & x ∈ Ω. \end{matrix}  (3.1b)

Existence and uniqueness of a solution to (3.1) in the weak sense has been shown in [8] with

θ = D(s) ∈ L^{2}(J; H^{1}(Ω)), s ∈ L^{∞}(Ω_{T}), p ∈ L^{∞}(J; V̂),

where

V̂ = {v ∈ H^{1}(Ω) : v_{Ω} = 0}.

Also, it was shown under physically reasonable assumptions that u is bounded:

∥u∥_{L^{∞}(Ω_{T})} ≤ C.  (3.2)

Property (3.2) and the following assumptions are implicitly used in this paper: ϕ ∈ L^{∞}(Ω) satisfies that ϕ(x) ≥ ϕ_{∗} > 0, κ(x) is a bounded, symmetric, and uniformly positive definite matrix, i.e.,

0 < κ_{∗} ≤ |ξ|^{−2} ∑_{i,j=1}^{d} κ_{ij}(x)ξ_{i}ξ_{j} ≤ κ^{∗} < ∞, x ∈ Ω, ξ ≠ 0 ∈ ℝ^{d},

and λ(s) satisfies that

0 < λ_{∗} ≤ λ(s) ≤ λ^{∗} < ∞, s ∈ [0, 1].

Without loss of generality, we further assume that ϕ ≡ 1 (otherwise, we consider the new variable ŝ = ϕs instead of s and the subsequent analysis is the same). Also, the functions d, f_{w}, γ′_{1}, and γ_{2} are assumed to be bounded functions of s. Finally, all the functions of s need to be defined only on [0, 1].

We replace d by a positive d_{β} > 0 with d_{β} → d in some sense as β → 0; a specific example of d_{β} will be given at the end of this section. For given d_{β} > 0, define

D_{β}(s) = ∫_{0}^{s} d_{β}(ξ)dξ.

The corresponding non-degenerate differential system is given by

\begin{matrix}∇ · u_{β} = q, u_{β} = −(λ(s_{β})∇p_{β} + γ′_{1}(s_{β})) & in Ω_{T} , \\ ∂_{t}s_{β} − ∇ · \{∇D_{β}(s_{β}) − f_{w}(s_{β})u_{β} + γ_{2}(s_{β})\} = q_{w} & in Ω_{T} , \end{matrix}  (3.3a)

with the boundary and initial conditions

\begin{matrix}u_{β} · ν = φ_{1}(x, t), & (x, t) ∈ Γ × J, \\ −(∇D_{β}(s_{β}) − f_{w}(s_{β})u_{β} + γ_{2}(s_{β})) · ν = φ_{2}(x, t), & (x, t) ∈ Γ × J, \\ s_{β}(x, 0) = s_{0}(x), & x ∈ Ω. \end{matrix} (3.3b)

We now determine in what manner (p_{β}, u_{β}, s_{β}) → (p, u, s) as β → 0. Toward that end, we make the assumption

\begin{matrix}∥λ(s_{1}) − λ(s_{2})∥_{L^{2}(Ω)}^{2} + ∥γ′_{1}(s_{1}) − γ′_{1}(s_{2})∥_{L^{2}(Ω)}^{2} \\ +∥f_{w}(s_{1}) − f_{w}(s_{2})∥_{L^{2}(Ω)}^{2} + ∥γ_{2}(s_{1}) − γ_{2}(s_{2})∥_{L^{2}(Ω)}^{2} \\ ≤ C(D(s_{1}) − D(s_{2}), s_{1} − s_{2}),  0 ≤ s_{1}, s_{2} ≤ 1. \end{matrix}  (3.4)

A sufficient condition for (3.4) to hold will be described later in this section. The proof of Lemma 3.1 and Theorem 3.2 can be found in [15].

Lemma 3.1. Let (p, u, s) and (p_{β}, u_{β}, s_{β}) solve (3.1) and (3.3), respectively. Then there is a constant C independent of β such that

∥p−p_{β}∥_{L^{2}(Ω)} +∥u−u_{β}∥_{L^{2}(Ω)} ≤ C{∥λ(s)−λ(s_{β})∥_{L^{2}(Ω)} +∥γ_{1}(s)−γ_{1}(s_{β})∥_{L^{2}(Ω)}}.

Theorem 3.2. Assume that d_{β} ≥ d and conditions (2.3) and (3.4) are satisfied. Let (p, u, s) and (p_{β}, u_{β}, s_{β}) solve (3.1) and (3.3), respectively. Then

\begin{matrix}∥p − p_{β}∥_{L^{2}(Ω_{T})}^{2} + ∥u − u_{β}∥_{L^{2}(Ω_{T})}^{2} + ∥s − s_{β}∥_{L^{∞}(J;H^{−1}(Ω))}^{2} \\ + ∫_{J}(D_{β}(s) − D_{β}(s_{β}), s − s_{β})dτ ≤ C(β), \end{matrix}  (3.5)

where C(β) = C∥D(s) − D_{β}(s)∥_{L^{∞}[0,1]}^{γ}.

Corollary 3.3. Under the assumptions of Theorem 3.2, we have

\begin{matrix}∥D_{β}(s) − D_{β}(s_{β})∥_{L^{2}(Ω_{T})} ≤ C(β)^{1/2}, \\ ∥s − s_{β}∥_{L^{μ+2}(Ω_{T})} ≤ C(β)^{1/(μ+2)}. \end{matrix}

The first result follows from (3.5) and the obvious inequality

(D(s_{1}) − D(s_{2}))^{2} ≤ ∥d∥_{L^{∞}[0,1]}(D(s_{1}) − D(s_{2}))(s_{1} − s_{2}), 0 ≤ s_{1}, s_{2} ≤ 1,

while the second result follows from (2.4) and (3.5).

As in [19], we now consider a specific example of the regularization d_{β} given by

d_{β}(s) = max{d(s), β^{μ}}.  (3.6)

Note that

∥D(s) − D_{β}(s)∥_{L^{∞}[0,1]} ≤ Cβ^{μ+1},

for β small enough, so

C(β) ≤ Cβ^{μ+2}.  (3.7)

We end this section with a remark on condition (3.4). Let η represent one of the quantities λ, f_{w}, γ′_{1}, and γ_{2}. It is clear that if η satisfies that

|η(s_{1}) − η(s_{2})|^{2} ≤ C(D(s_{1}) − D(s_{2}))(s_{1} − s_{2}), a. e. 0 ≤ s_{1}, s_{2} ≤ 1, (3.8)

then assumption (3.4) is true for η. A necessary and sufficient condition for (3.8) to hold is that [14]

|η_{s}|^{2} ≤ Cd(s),  a. e. s ∈ [0, 1].  (3.9)

Inequality (3.9) means that η_{s} vanishes with d. Below we give the conditions on η so that (3.8) or (3.9) holds. The proof of the next proposition can be found in [9] or [19].

Proposition 3.4. Assume that η ∈ C^{1}[0, 1], η_{s}(0) = η_{s}(1) = 0, η_{s} is Lipschitz continuous at 0 and 1, and assumption (2.3) is satisfied. Then there is a constant C > 0 such that (3.8) holds.

4 Perturbation Analysis

In this section we report a formal application of the perturbation method for the saturation equation

∂_{t}s − ∇ · {d(s)∇s} + ∇ · {f_{w}(s)u} − ∇ · γ_{2}(s) = q_{w} .  (4.1)

The perturbation method in [21] is applied to analyze numerical solutions of this problem. We assume that the numerical method we will develop produces a solution in the asymptotic form

s_{ϵ} ∼ s_{1} + ϵs_{2} + ⋯,  (4.2)

where s_{1} is the exact solution of (4.1), s_{2} is a smooth function, and ϵ is a small constant. Substituting (4.2) into (4.1), we see that

\begin{matrix}∂_{t}s_{1} + ϵ∂_{t}s_{2} − ∇ · \{d(s_{1} + ϵs_{2})∇(s_{1} + ϵs_{2})\} \\ +∇ · \{f_{w}(s_{1} + ϵs_{2})u\} − ∇ · γ_{2}(s_{1} + ϵs_{2}) + ⋯ ∼ q_{w} . \end{matrix}

Since s_{1} satisfies (4.1), we have

ϵ∂_{t}s_{2} − ∇ · {d(s_{1} + ϵs_{2})∇(s_{1} + ϵs_{2}) − d(s_{1})∇s_{1}}
+∇ · {(f_{w}(s_{1} + ϵs_{2}) − f_{w}(s_{1})) u}
−∇ · {γ_{2}(s_{1} + ϵs_{2}) − γ_{2}(s_{1})} + ⋯ ∼ 0.

That is,

\begin{matrix}∂_{t}s_{2} − ∇ · \{s_{2} \frac{d(s_{1}+ϵs_{2})−d(s_{1})}{ϵs_{2}}∇s_{1}\} − ∇ · \{d(s_{1} + ϵs_{2})∇s_{2}\} \\ +∇ · \{s_{2} \frac{f_{w}(s_{1}+ϵs_{2})−f_{w}(s_{1})}{ϵs_{2}}u\} − ∇ · \{s_{2} \frac{γ_{2}(s_{1}+ϵs_{2})−γ_{2}(s_{1})}{ϵs_{2}}\} ∼ 0. \end{matrix} (4.3)

Assuming that d, f_{w}, γ_{2} ∈ C^{1}[0, 1] as in Proposition 3.4, it follows from (4.3) as ϵ → 0 that

\begin{matrix}∂_{t}s_{2} − ∇ · \{s_{2}d_{s}(s_{1})∇s_{1}\} − ∇ · \{d(s_{1})∇s_{2}\} \\ +∇ · \{s_{2}f_{ws}(s_{1})u\} − ∇ · \{s_{2}γ_{2s}(s_{1})\} ∼ 0; \end{matrix}

i.e.,

∂_{t}s_{2} − ∇ · {d(s_{1})∇s_{2}} + ∇ · {s_{2} (f_{ws}(s_{1})u − γ_{2s}(s_{1}) − d_{s}(s_{1})∇s_{1})} ∼ 0 . (4.4)

By (3.9), note that f_{ws} and γ_{2s} vanish with d, so (4.4) reduces to the following equation near the degeneracy of the function d(s) with some assumptions for the functions γ_{2} and f_{w}:

∂_{t}s_{2} + ∇ · {s_{2} (−d_{s}(s_{1})∇s_{1})} ∼ d_{s}(s_{1})∇s_{1} · ∇s_{2} .

Hence the above formal analysis indicates that the behavior of errors close to the degeneracy is exponential and −d_{s}(s_{1})∇s_{1} determines the gross rate of the errors in time. This is the case as shown in our numerical experiments later.

5 Finite Element Approximation

For notational convenience, we consider the case of φ_{1} ≡ 0 in the analysis below; otherwise, φ_{1} can be incorporated into the differential equation, or the later mixed finite element method can be handled by introducing the space of Lagrange multipliers [12].

For d = 2 or 3, let

H(div, Ω) = {v ∈ (L^{2}(Ω))^{d} : ∇ · v ∈ L^{2}(Ω)},

and

V = {v ∈ H(div, Ω) : v · ν = 0 on Γ}, W = {w ∈ L^{2}(Ω) : w_{Ω} = 0}.

For 0 < h_{p}, h < 1, let T_{h_{p}} and T_{h} be regular partitions into elements, say, simplexes, rectangular parallelepipeds, and/or prisms. Associated with T_{h}, let M_{h} ⊂ H^{1}(Ω) be a standard C^{0}-finite element space associated with T_{h} such that

inf_{v_{h}∈M_{h}} ∥v − v_{h}∥_{L^{r}(Ω)} ≤ Ch^{2}∥v∥_{W^{2,r}(Ω)},  1 < r < ∞.  (5.1)

In this paper we only consider lowest-order C^{0}-finite elements such that (5.1) is satisfied; due to a lack of regularity on the solution, no improvements in the asymptotic convergence rate result from taking higher order finite element spaces.

Associated with the partition T_{h_{p}}, let V_{h} × W_{h} = V_{h_{p}} × W_{h_{p}} ⊂ V × W be the Raviart-Thomas-Nedelec [26, 24], the Brezzi-Douglas-Fortin-Marini [4], the Brezzi-Douglas-Marini [5] (if d = 2), the Brezzi-Douglas-Durán-Fortin [3] (if d = 3), or the Chen-Douglas [10] mixed finite element space of index such that the approximation properties below are satisfied:

inf_{v_{h}∈V_{h}} ∥v − v_{h}∥_{L^{2}(Ω)} ≤ Ch_{p}^{l}∥v∥_{H^{l}(Ω)},  0 ≤ l ≤ k^{∗} + 1,
inf_{v_{h}∈V_{h}} ∥∇ · (v − v_{h})∥_{L^{2}(Ω)} ≤ Ch_{p}^{l}∥∇ · v∥_{H^{l}(Ω)}, 0 ≤ l ≤ k^{∗∗},
inf_{w_{h}∈W_{h}} ∥w − w_{h}∥_{L^{2}(Ω)} ≤ Ch_{p}^{l}∥w∥_{H^{l}(Ω)},  0 ≤ l ≤ k^{∗∗},

(5.2)

where k^{∗∗} = k^{∗} + 1 for the first two spaces, k^{∗∗} = k^{∗} for the second two spaces, and both cases are included in the last space.

5.1 A semidiscrete approximation

Let

a(s) = (κλ(s))^{−1}, γ_{1}(s) = −γ′_{1}(s)λ^{−1}(s).

As mentioned before, the functions of s are defined on [0, 1]. In this section the possibility that s ∉ [0, 1] is allowed. Following [11, 20, 27], the function d_{β} is extended as follows:

d_{β}(s) = \{ \begin{matrix}d_{β}(1) & if s ≥ 1, \\ d_{β}(−s) & if s ≤ 0. \end{matrix}

Let η represent one of the quantities a, f_{w}, γ_{1}, and γ_{2}. We extend η by

η(s) = \{ \begin{matrix}η(1) & if s ≥ 1, \\ η(0) & if s ≤ 0. \end{matrix}

By the above extension it follows that D_{β}(s) is now strictly increasing in s on the real line because D′_{β}(s) = d_{β}(s) > 0 for any β > 0, so D_{β} has a C^{1} inverse function S_{β}. Finally, let P_{h} indicate the L^{2}-projection operator onto M_{h}.

As remarked before, the pressure equation is approximated by the mixed finite element method. For each t ∈ J̄, the mixed finite element solution (u_{h}(·, t), p_{h}(·, t)) ∈ V_{h} × W_{h} satisfies

\begin{matrix}(∇ · u_{h}, w) = (q, w), & ∀w ∈ W_{h}, \\ (a(s_{h})u_{h}, v) − (p_{h}, ∇ · v) = (γ_{1}(s_{h}), v), & ∀v ∈ V_{h}, \end{matrix}  (5.3)

where s_{h} is determined below. First, for each t ∈ J we define θ_{h}(·, t) ∈ M_{h} by

\begin{matrix}(∂_{t}S_{β}(θ_{h}), v) + (∇θ_{h} − f_{w}(S_{β}(θ_{h}))u_{h} \\ +γ_{2}(S_{β}(θ_{h})), ∇v) + (φ_{2}, v)_{Γ} = (q_{w}, v), ∀v ∈ M_{h}, \end{matrix}  (5.4)

with the initial approximation

P_{h}S_{β}(θ_{h}(0)) = P_{h}s_{0}.  (5.5)

Now, we determine s_{h} by s_{h} = S_{β}(θ_{h}), which approximates s_{β}.

Notice that approximating D_{β}(s_{β}) by θ_{h}, then s_{β} by S_{β}(θ_{h}), yields a higher rate of convergence than approximating s_{β} by an element in M_{h} directly [14, 20]. Also, note that if bases are introduced in V_{h}, W_{h}, and M_{h}, equations (5.3)–(5.5) can be written as a nonlinear system of ordinary differential equations for s_{h} (after substituting (5.3) into (5.4)). With our assumptions on the data, this nonlinear system can be shown to have a unique solution from the fundamental theorem of ordinary differential equations [11]. An error analysis for (5.3)–(5.5) is given in the second paper [16].

5.2 A fully discrete approximation

For each positive integer N, let 0 = t^{0} < t^{1} < ⋯ < t^{N} = T be a partition of J into subintervals J^{n} = (t^{n−1}, t^{n}] with length Δt^{n} = t^{n} − t^{n−1}, 1 ≤ n ≤ N. Also, set v^{n} = v(·, t^{n}). Finally, we indicate the time difference operator by

∂v^{n} = \frac{v^{n} − v^{n−1}}{Δt^{n}} , 1 ≤ n ≤ N.

Now, the fully discrete approximation is given as follows. For any 1 ≤ n ≤ N, find (u_{h}^{n}, p_{h}^{n}) ∈ V_{h} × W_{h} such that

\begin{matrix}(∇ · u_{h}^{n}, w) = (q^{n}, w) & ∀w ∈ W_{h}, \\ (a(s_{h}^{n})u_{h}^{n}, v) − (p_{h}^{n}, ∇ · v) = (γ_{1}(s_{h}^{n}), v) & ∀v ∈ V_{h}, \end{matrix}  (5.6)

where s_{h} = S_{β}(θ_{h}) and for each n, θ_{h}^{n} ∈ M_{h} satisfies

\begin{matrix}(∂S_{β}(θ_{h}^{n}), v) + (∇θ_{h}^{n} − f_{w}(S_{β}(θ_{h}^{n}))u_{h}^{n} \\ +γ_{2}(S_{β}(θ_{h}^{n})), ∇v) + (φ_{2}^{n}, v)_{Γ} = (q_{w}^{n}, v), ∀v ∈ M_{h}, \end{matrix}  (5.7)

with the initial approximation given as in (5.5).

1/Δx p − p_{h} rate for p u − u_{h} rate for u
10  0.0585  -  0.0280  -
20  0.0277  1.07  0.0073  1.94
40  0.0135  1.04  0.0019  1.97
80  0.0066  1.02  0.0005  1.98

Table 1a. The error estimates for p and u.

The remark on existence and uniqueness of (5.6) and (5.7) can be made as above [12]. Also, an error analysis for this approximation is presented in [16].

1/Δx s − s_{h} for rate for s − s_{h} for rate for
β = 0  β = 0  β = β_{0}  β = β_{0}
10  0.0858  -  0.5443  -
20  0.0452  0.91  0.2931  0.90
40  0.0272  0.74  0.1236  1.24
80  0.0162  0.75  0.0530  1.22

Table 1b. The error estimates for s in the first example.

6 Numerical Results

The numerical experiments are presented to show convergence of our approximation scheme, to demonstrate the qualitative behavior of error estimates, to compare the present regularization technique with the un-regularized version used in [14], and to complement the perturbation analysis in the fourth section. Toward that end, we consider the pressure-saturation system of the form

\begin{matrix}−∇ · (λ(s)∇p) = q & in Ω_{T} , \\ ∂_{t}s − ∇ · \{d(s)∇s + λ_{w}(s)∇p\} = q_{w} & in Ω_{T} . \end{matrix}  (6.1)

For simplicity we focus on the Dirichlet boundary conditions

p = p_{D}, s = s_{D} on Γ × J.

The initial condition is

s(x, 0) = s_{0}(x), x ∈ Ω.

Example 1. The domain Ω is taken to be the unit square and other data are chosen as follows:

λ(s) = 1, d(s) = s(1 − s), λ_{w}(s) = 1.

In the first example the exact solution of system (6.1) is of the form

p(x, t) = sin(πx) sin(πy), s(x, t) = t sin(πx) sin(πy).

The boundary and initial data coincide with the exact solution on the boundary and at the initial time. Also, the functions in the right-hand side of system (6.1) result from the exact solution. The numerical experiments reported here are mainly to show the qualitative behavior of the finite element approximation for the degenerate saturation equation. This is why we consider the case where λ(s) = 1, so the pressure equation is decoupled from it. We also carried out experiments with λ(s) depending on s, which are not reported here and have similar results to those illustrated here. Furthermore, note that the data satisfy assumptions (2.3) with μ_{1} = μ_{2} = μ = 1 and (3.4).

Uniform partitions of Ω into triangles are used, with Δx = Δy as the lengths in the x and y directions, and the lowest-order Raviart-Thomas mixed finite element on the triangles [26] are exploited. The mixed method is used for the pressure equation and the standard finite element method with the backward Euler scheme for the time differentiation term is utilized for the saturation equation, as in (5.6) and (5.7). The time step is taken to be proportional to the space step. The mixed finite element method is implemented as in [7].

Error estimates and convergence rates in the L^{∞}-norm for the approximations to the pressure and velocity u = −λ(s)∇p are presented in Table 1a. From this table we see that the mixed method is first-order accurate for the pressure and second-order accurate for the velocity. That is, a superconvergence rate occurs for the velocity. The error estimates in the L^{∞}-norm for the saturation at t = 1 are described in Table 1b. We consider the cases where the regularization parameter β is taken to be zero or β_{0} ≡ Ch^{μ_{0}} with μ_{0} given by

μ_{0} = \frac{4}{3μ + 2},

which is chosen according to the error analysis performed in the paper [16]. The convergence rates in the case of β = β_{0} are better than those in the case of β = 0. Also, the convergence rates in the former case coincide with the theoretical results obtained in [16]. Finally, the error with Δx = 1/80 for s_{h} with β = β_{0} is shown as Figure 1 in a separate file. It can be seen that maximum errors occur in the center where d(s) is zero and near the boundary of the domain where d(s) is close to zero. This complements the perturbation analysis in the forth section.

Figure 2: The error s − s_{h} for β = β_{0} in the second example.

Example 2. In this example, we simulate a relatively simple two-phase flow problem. The data are given as follows:

k_{rw} = s,  k_{ro} = 1 − s,  μ_{w} = 0.5 cp, μ_{o} = 2 cp,
ρ_{w} = 1 g/cm^{3}, ρ_{o} = 0.7 g/cm^{3}, ϕ = 0.2,  k = 0.05 darcy.

Moreover, the function p_{c} is given by

p_{c}(s) = 1 − s .

With these, the system of equations (3.1a) now reduces to

∇ · u = q ,
u = −a(s)(∇p − G(s)) ,
ϕ∂_{t}s − ∇ · {d(s)∇s − f_{w}(s)u − b(s)} = q_{w} ,

where

a(s) = 0.0125(3s + 1),  G(s) = (0.7 + 3.3s)g̃,
d(s) = 0.1s(1 − s)/(3s + 1),  f_{w}(s) = 4s/(3s + 1),
b(s) = 0.03s(1 − s)/(3s + 1)g̃ ,

where g̃ is the gravity vector.

1/Δx p − p_{h} for rate for p − p_{h} for rate for
β = 0  β = 0  β = β_{0}  β = β_{0}
10  0.075626  -  0.069228  -
20  0.036445 1.0532 0.031686 1.1275
40  0.017600 1.0501 0.014920 1.0866
80  0.008615 1.0307 0.007110 1.0693

Table 2a. The error estimates for p in the second example.

1/Δx s − s_{h} for rate for s − s_{h} for rate for
β = 0  β = 0  β = β_{0}  β = β_{0}
10  0.2905  -  0.2423  -
20  0.1153  1.3331  0.0969  1.3222
40  0.0528  1.1268  0.0430  1.1722
80  0.0253  1.0614  0.0203  1.0829

Table 2b. The error estimates for s in the second example.

The numerical results corresponding to those in Example 1 are given in Table 2 and Figure 2. Similar observations can be made here.

Example 3. In the final example, we test a more physically adequate set of data. We simulate a two-phase flow problem [2]. The function p_{c}(s) is given by

p_{c}(s) = (1 − s) {γ(s^{−1} − 1) + θ} ,

where

γ = 20, 000 dynes/cm^{2}, θ = 100 dynes/cm^{2} .

The relative permeabilities are defined by

k_{ro} = \{ \begin{matrix}0 & if s > s_{o} , \\ s_{o}^{−2}(s_{o} − s)^{2} & if 0 ≤ s ≤ s_{o} . \end{matrix}

and

k_{rw} = \{ \begin{matrix}(s − s_{rw})^{2}(1 − s_{rw})^{−2} & if s ≥ s_{rw} , \\ 0 & if 0 ≤ s < s_{rw} , \end{matrix}

where

ϕ = 0.2,  k = 0.05 darcy, μ_{w} = 0.5 cp,
μ_{o} = 2 cp,  ρ_{w} = 1 g/cm^{3}, ρ_{o} = 0.7 g/cm^{3} ,
s_{o} = 1 − s_{ro}, s_{ro} = 0.15,  s_{rw} = 0.2 .

The domain Ω and boundary and initial conditions are taken as in Example 1. We consider the normalized water saturation

s = \frac{s_{w} − s_{rw}}{1 − s_{ro} − s_{rw}} .

The functions f_{w}(s) and d(s) are illustrated in Figure 3.

Figure 3: Normalized curves of f_{w}(s) (A), f′_{w}(s) (B), d′(s) (C), and d(s) (D).

1/Δx s − s_{βh} rate for s u − u_{βh} rate for u
10  0.0025  -  0.0321  -
20  0.0015  0.7370  0.0124  1.37
40  7.7303e-04 0.9564  0.0050  1.3
80  3.9415e-04 0.9870  0.0022  1.2

Table 3. The error estimates for s and u in the third example.

The error estimates and convergence rates in the L^{∞}-norm for the approximations to the saturation and velocity at t = 0.01 are presented in Table 3. The table shows that both are first-order accurate. Also, we can see from Figure 3 that maximum errors occur when the saturation is close to zero.

Figure 4: The error s − s_{h} for β = β_{0} in the third example.

References

[1] S. N. Antontsev, On the solvability of boundary value problems for degenerate two-phase porous flow equations, Dinamika Splos̆noĭ Sredy Vyp. 10 (1972), 28–53, in Russian.

[2] J. Bear, Dynamics of Fluids in Porous Media, Dover, New York, 1972.

[3] F. Brezzi, J. Douglas, Jr., R. Durán, and M. Fortin, Mixed finite elements for second order elliptic problems in three variables, Numer. Math. 51 (1987), 237–250.

[4] F. Brezzi, J. Douglas, Jr., M. Fortin, and L. Marini, Efficient rectangular mixed finite elements in two and three space variables, RAIRO Modèl. Math. Anal. Numér 21 (1987), 581–604.

[5] F. Brezzi, J. Douglas, Jr., and L. Marini, Two families of mixed finite elements for second order elliptic problems, Numer. Math. 47 (1985), 217– 235.

[6] G. Chavent and J. Jaffré, Mathematical Models and Finite Elements for Reservoir Simulation, North-Holland, Amsterdam, 1978.

[7] Zhangxin Chen, Equivalence between and multigrid algorithms for nonconforming and mixed methods for second order elliptic problems, East-West J. Numer. Math. 4 (1996), 1–33.

[8] Zhangxin Chen, Degenerate two-phase incompressible flow I: Existence, uniqueness and regularity of a weak solution, SMU Math Report 97–10, Dallas, Texas.

[9] Zhangxin Chen, Degenerate two-phase incompressible flow IV: Regularity, stability and stabilization, to appear.

[10] Zhangxin Chen and J. Douglas, Jr., Prismatic mixed finite elements for second order elliptic problems, Calcolo 26 (1989), 135–148.

[11] Zhangxin Chen, M. Espedal, and R. E. Ewing, Continuous-time finite element analysis of multiphase flow in groundwater hydrology, Appl. Math. 40 (1995), 203–226.

[12] Zhangxin Chen and R. E. Ewing, Fully-discrete finite element analysis of multiphase flow in groundwater hydrology, SIAM J. Numer. Anal. 34 (1997), 2228–2253.

[13] Zhangxin Chen and R. E. Ewing, Mathematical analysis for reservoir models, SIAM J. Math. Anal., 30 (1999), 431-453.

[14] Zhangxin Chen and R. E. Ewing, Degenerate two-phase incompressible flow II: Optimal error estimates, Numer. Math., to appear.

[15] Zhangxin Chen and N. L. Khlopina, Degenerate two-phase incompressible flow problems I: regularization and numerical results, Comm. Appl. Anal., to appear.

[16] Zhangxin Chen and N. L. Khlopina, Degenerate two-phase incompressible flow problems II: error estimates, Comm. Appl. Anal., to appear.

[17] P. G. Ciarlet, The finite Element Method for Elliptic Problems, North– Holland, Amsterdam, 1978.

[18] R. E. Ewing (ed.), The Mathematics of Reservoir Simulation, SIAM, Philadelphia, 1983.

[19] K. Fadimba and R. Sharpley, A priori estimates and regularization for a class of porous medium equations, Nonlinear World 2 (1995), 13–41.

[20] K. Fadimba and R. Sharpley, Galerkin finite element method for a class of porous medium equations, Preprint, 1997.

[21] M. Holmes, Introduction to Perturbation Methods, Springer-Verlag, New York, Vol. 20, 1995.

[22] J. Kevorkian and J. D. Cole, Multiple Scale and Singular Perturbation Methods, Springer-Verlag, New York, Vol. 114, 1996.

[23] N. L. Khlopina, Finite element methods for degenerate two-phase incompressible flow problems, Ph.D. Thesis, Southern Methodist University, Dallas, Texas, 1999.

[24] J. Nedelec, Mixed finite elements in ℝ^{3}, Numer. Math. 35 (1980), 315–341.

[25] D. W. Peaceman, Fundamentals of Numerical Reservoir Simulation, Elsevier, New York, 1977.

[26] P. Raviart and J. Thomas, A mixed finite element method for second order elliptic problems, Lecture Notes in Mathematics, vol. 606, Springer, Berlin, 1977, pp. 292–315.

[27] M. Rose, Numerical Methods for flow through porous media I, Math. Comp. 40 (1983), 437–467.

[28] D. Smylie, A near optimal order approximation to a class of two sided nonlinear degenerate parabolic partial differential equations, Ph. D. Thesis, University of Wyoming, Laramie, 1989.

Zhangxin Chen (e-mail: zchen@dragon.math.smu.edu) Natalia L. Khlopina (e-mail: khlopina@mail.smu.edu) Department of Mathematics Box 750156, Southern Methodist University Dallas, TX 75275-0156 USA

Note. Figure 1 is available as a PostScript file in the same directory where the TeX DVI and PDF files are.