Fully Discrete Finite Element Analysis of Multiphase Flow in Groundwater Hydrology.
This paper deals with the development and analysis of a fully discrete finite element method for a nonlinear differential system for describing an air-water system in groundwater hydrology. The nonlinear system is written in a fractional flow formulation, i.e., in terms of a saturation and a global pressure. The saturation equation is approximated by a finite element method, while the pressure equation is treated by a mixed finite element method. The analysis is carried out first for the case where the capillary diffusion coefficient is assumed to be uniformly positive, and is then extended to a degenerate case where the diffusion coefficient can be zero. It is shown that error estimates of optimal order in the $L^2$-norm and almost optimal order in the $L^\infty$-norm can be obtained in the nondegenerate case. In the degenerate case we consider a regularization of the saturation equation by perturbing the diffusion coefficient. The norm of error estimates depends on the severity of the degeneracy in diffusivity, with almost optimal order convergence for nonsevere degeneracy. Implementation of the fractional flow formulation with various nonhomogeneous boundary conditions is also discussed. Results of numerical experiments using the present approach for modeling groundwater flow in porous media are reported.
Introduction
. In this paper we develop and analyze a fully-discrete finite element
procedure for solving the flow equations for an air-water system in groundwater
hydrology,
@t
kk rff
porous medium, OE and k are the porosity and absolute
permeability of the porous system, ae ff , s ff , p ff , u ff , and - ff are the density, saturation,
pressure, volumetric velocity, and viscosity of the ff-phase, f ff is the source/sink term,
k rff is the relative permeability of the ff-phase, and g is the gravitational, downward-
pointing, constant vector.
Flow simulation in groundwater reservoirs has been extensively studied in past
years (see, e.g., [26], [28] and the bibliographies therein). However, in most previous
works the air-phase equation is eliminated by the assumption that the air-phase
remains essentially at atmospheric pressure. This assumption, as mentioned in [13],
is reasonable in most cases because the mobility of air is much larger than that of
water, due to the viscosity difference between the two fluids. When the air-phase
pressure is assumed constant, the air-phase mass balance equation can be eliminated
and thus only the water-phase equation remains. Namely, the Richards equation is
used to model the movement of water in groundwater reservoirs. However, it provides
* Department of Mathematics and the Institute for Scientific Computation, Texas A&M Uni-
versity, College Station, Partly supported by the Department of Energy under contract
DE-ACOS-840R21400. email: zchen@isc.tamu.edu, ewing@ewing.tamu.edu.
no information on the motion of air. If contaminant transport is the main concern
and the contaminant can be transported in the air-phase, the air-phase needs to be
included to determine the advective component of air-phase contaminant transport
[7]. Furthermore, the dynamic interaction between the air and water phases is also
important in vapor extraction systems. Hence in these cases the coupled system of
nonlinear equations for the air-water system must be solved. It is the purpose of this
paper that is to develop and analyze a finite element procedure for approximating
the solution of the coupled system of nonlinear equations for the air-water system in
groundwater hydrology.
In petroleum reservoir simulation the governing equations that describe fluid flow
are usually written in a fractional flow formulation, i.e., in terms of a saturation and
a global pressure [1], [8]. The main reason for this fractional flow approach is that
efficient numerical methods can be devised to take advantage of many physical properties
inherent in the flow equations. However, this pressure-saturation formulation
has not yet achieved application in groundwater hydrology. In petroleum reservoirs
total flux type boundary conditions are conveniently imposed and often used, but in
groundwater reservoirs boundary conditions are very complicated. The most commonly
encountered boundary conditions for a groundwater reservoir are of first-type
(Dirichlet), second-type (Neumann), third-type (mixed), and "well" type [8]. The
problem of incorporating these nonhomogeneous boundary conditions into the fractional
flow formulation has been a challenge [12]. In particular, in using the fractional
flow approach a difficulty arises when the Dirichlet boundary condition is imposed for
one phase (e.g. air) and the Neumann type is used for another phase (e.g. water).
This paper follows the fractional flow formulation. Based on this approach, we
develop a fully-discrete finite element procedure for the saturation and pressure equa-
tions. The saturation equation is approximated by a Galerkin finite element method,
while the pressure equation is treated by a mixed finite element method. It is well
known that the physical transport dominates the diffusive effects in incompressible
flow in petroleum reservoirs. In the air-water system studied here, the transport
again dominates the entire process. Hence it is important to obtain good approximate
velocities. This motivates the use of the parabolic mixed method, as in [17], in the
computation of the pressure and the velocity. Also, due to its convection-dominated
feature, more efficient approximate procedures should be used to solve the saturation
equation. However, since this is the first time to carry out an analysis for the
present problem, it is of some importance to establish that the standard finite element
method for this model converges at an asymptotically optimal rate for smooth
problems. Characteristic Petrov-Galerkin methods based on operator splitting [20],
transport diffusion methods [32], and other characteristic based methods will be considered
in forthcoming papers.
The main part of this paper deals with an asymptotical analysis for the fully
discrete finite element method for the first-type and second-type boundary conditions
where p ffD and d ff are given functions,
being disjoint, and
- is the outer unit normal to @ We point out that petroleum reservoir simulation
is different from groundwater reservoir simulation. The flow of two incompressible
fluids (e.g. water and oil) is usually considered in the former case, while the latter
system consists of the air and water phases. Consequently, the finite element analyses
for these two cases differ. As shown here, compressibility and combination of the
boundary conditions (1.3) and (1.4) complicate error analyses. Indeed, if optimality is
to be preserved for the finite element method, the standard error argument just fails
unless we work with higher order time-differentiated forms of error equations, which
require properly scaling initial conditions. Also, we mention that a slightly compressible
miscible displacement problem was treated in [14], [18], [23], [33]; however, only
the single phase was handled, gravitational terms were omitted, and total flux type
boundary conditions were assumed. Furthermore, the so-called "quadratic" terms in
velocity were neglected. The dropping of these quadratic terms may not be valid near
wells, and so the miscible displacement model was oversimplified both physically and
mathematically. The analysis of this paper includes these terms. Finally, only the
Raviart-Thomas mixed finite element spaces [34] have been considered in these earlier
papers. We are here able to discuss all existing mixed spaces.
The error analysis is given first for the case where the capillary diffusion coefficient
is assumed to be uniformly positive. In this case, we show error estimates of optimal
order in the L 2 -norm and almost optimal order in the L 1 -norm. Then we treat a
degenerate case where the diffusion coefficient vanishes for two values of saturation.
In the degenerate case we consider a regularization of the saturation equation by
perturbing the diffusion coefficient to obtain a nondegenerate problem with smooth
solutions. It is shown that the regularized solutions converge to the original solution
as the perturbation parameter goes to zero with specific convergence rates given. The
norm of error estimates depends on the severity of the degeneracy in diffusivity, with
almost optimal order convergence for the degeneracy under consideration.
The rest of this paper is concerned with implementation of the fractional flow
formulation with various nonhomogeneous boundary conditions. We show that all
the commonly encountered boundary conditions can be incorporated in the fractional
flow formulation. Normally the "global" boundary conditions are highly nonlinear
functions of the physical boundary conditions for the original two flow phases. This
means that we have to iterate on these global boundary conditions as part of the solution
process. We here develop a general solution approach to handle these boundary
conditions. Results of numerical experiments using the present approach for modeling
groundwater flow are reported here.
The paper is organized as follows. In x2, we define a fractional flow formulation for
equations (1.1)-(1.4). Then, in x3 we introduce weak forms of the pressure-saturation
equations, and in x4 a fully-discrete finite element procedure for solving these equa-
tions. An asymptotical analysis is given in x5 and x6 for the nondegenerate case
and the degenerate case, respectively. Finally, in x7 we discuss implementation of
various nonhomogeneous boundary conditions and present the results of numerical
experiments.
2. A pressure-saturation formulation. In addition to (1.1)-(1.4), we impose
the customary property that the fluid fills the volume:
and define the capillary pressure function p c by
Introduce the phase mobilities
and the total mobility
To devise our numerical method, it is important to choose a reasonable set of dependent
variables. Since equal to the water residual saturation [3],
pw cannot generally be expected to lie in any Sobolev space. Air being a continuous
phase implies that p a is well behaved. Hence, as mentioned in the introduction, we
define the global pressure [1] with
sc
-w
d-
d-
Z pc (s)\Gamma -w
c (-)
The integral in the right-hand side of (2.3) is well defined [1], [8].
As usual, assume that ae ff depends on p [8]. Then we define the total velocity
where
g:
Now it can be easily seen that
(2.5a)
where q
Consequently,
Equations (1.1) and (1.2) can be manipulated using (2.1)-(2.6) to have the pressure
equation
@t
a
ff=wae ff
OEs ff
@ae ff
@t
and the saturation equation
OE @sw
@t
\Gammas w
@t
ae w
OEs w
@t
Terms of the form u ff \Delta rae ff , neglected in compressible miscible
displacement problems [14], [18], [23], [33]. The dropping of these terms may not
be valid near wells. Also, if they are neglected, the model may not be qualitatively
equivalent to the usual formulation of two phase flow. Hence we keep them in this
paper. However, the water phase is usually assumed to be incompressible. With the
incompressibility of the water phase and the following notation:
c(s;
ae a
dae a
dp
ds
ae a
dae a
dp
~
ae w
ae a
dae a
dp
ae)
f(s; p) =ae a
dae a
dp k-w q a (rp c \Gamma ~
f a
ae a
ae w
@t
equations (2.7) and (2.8) can be now written as
c(s; p) @p
@t
OE @s
@t
@t
The boundary conditions for the pressure-saturation equations become
(D(s)rs
where s D and pD are the transforms of pwD and paD by (2.2) and (2.3), and ~
The model is completed by specifying the initial conditions
The later analysis for the nondegenerate case in x5 is given under a number of
assumptions. First, the solution is assumed smooth; i.e., the external source terms
are smoothly distributed, the coefficients are smooth, the boundary and initial data
satisfy the compatibility condition, and the domain has at least the regularity required
for a standard elliptic problem to have H
2(\Omega\Gamma671/64748/ y and more if error estimates
of order bigger than one are required. Second, the coefficients a(s), OE, and c(s; p) are
assumed bounded below positively:
Finally, the capillary diffusion coefficient D(s) is assumed to satisfy
While the phase mobilities can be zero, the total mobility is always positive [31].
The assumptions (2.18) and (2.19) are physically reasonable. Also, the present analysis
obviously applies to the incompressible case where c(s; In this case, the analysis
is simpler since we have an elliptic pressure equation instead of the parabolic equation
(2.9). Thus we assume condition (2.20) for the compressible case under consideration.
Next, although the reasonableness of the assumption (2.21) is discussed in [16], the
diffusion coefficient D(s) can be zero in reality. It is for this reason that section six
is devoted to consideration of the case where the solution is not required smooth and
the assumption (2.21) is removed. As a final remark, we mention that for the case
where point sources and sinks occur in a porous medium, an argument was given in
[22] for the incompressible miscible displacement problem and can be extended to the
present case.
3. Weak forms. To handle the difficulty associated with the inhomogeneous
Neumann boundary condition (2.13) in the analysis of the mixed finite element method,
let d be such that d \Delta
d and introduce the change of variable equations
(2.9)-(2.11). Then the homogeneous Neumann boundary condition holds for ~ u. Thus,
without loss of generality, we assume that ~
To be compatible, we also require
that this homogeneous condition holds when
In the two-dimensional case, let
while it is accordingly defined in the three-dimensional case as follows:
Also, set
The weak form of (2.9)-(2.11) on which the finite element procedure is based is given
below. Let is the time interval of interest. The mixed formulation
for the pressure is defined by seeking a pair of maps
2(\Omega\Gamma such that
(3.1a)
(c(s; p) @p
@t
the inner products (\Delta; \Delta) are to be interpreted to be in L
or (L 2
(\Omega\Gamma4 d , as appropriate, and h\Delta; \Deltai
denotes the duality between H 1=2
H \Gamma1=2 (\Gamma 1 ). The weak form for the saturation s : J
OE
@s
@t
where the boundary condition (2.15) is used. Finally, to treat the nonzero initial
conditions imposed on s and p in (2.16) and (2.17), we introduce the following transformations
7in (3.1) and (3.2):
where
we have zero initial conditions for s, p, and u. Hence, without loss
of generality again, we assume that
The reason for introducing these transformations to have zero initial conditions is to
validate equation (5.15) later.
4. Fully-discrete finite element procedures.
Let\Omega be a polygonal domain.
For partitions into ele-
ments, say, simplexes, rectangular parallelepipeds, and/or prisms. In both partitions,
we also need that adjacent elements completely share their common edge or face. Let
be a standard C 0 -finite element space associated with T h such
that
where hK is the norm in the Sobolev space W k;q (K)
(we omit K when K
=\Omega and kvk
(\Omega\Gamma be the Raviart-Thomas-Nedelec [34], [29], the Brezzi-Douglas-Fortin-Marini
[5], the Brezzi-Douglas-Marini [6] (if 2), the Brezzi-Douglas-Dur'an-Fortin [4] (if
or the Chen-Douglas [11] mixed finite element space associated with the
partition T hp of index such that the approximation properties below are satisfied:
kr
where h p;K for the first two spaces,
for the second two spaces, and both cases are included in the last space.
Finally, let ft n g nT
n=0 be a quasi-uniform partition of J with t
set \Deltat
We are now in a position to introduce our finite element procedure.
The fully-discrete finite element method is given as follows. The approximation
procedure for the pressure is defined by the mixed method for a pair of maps fu n
(ff(s
(4.5a)
(c(s
(4.5b)
and the finite element method for the saturation is given for s n
h )rs n
@t
The initial conditions satisfy
After startup, for equations (4.5) and (4.6) are computed as
follows. First, using s
h , and (4.5), evaluate fu n
g. Since it is linear, (4.5)
has a unique solution for each n [10], [27]. Next, using s
h g, and (4.6),
calculate s n
h . Again, (4.6) has a unique solution for \Deltat n sufficiently small for each n
[39].
We end this section with a remark. While the backward Euler scheme is used
in (4.5b) and (4.6), the Crank-Nicolson scheme and more accurate time stepping
procedures (see, e.g., [21]) can be used. The present analysis applies to these schemes.
5. An error analysis for the fully-discrete scheme. In this section we give a
convergence analysis for the finite element procedure (4.5) and (4.6) under assumption
(2.21). As usual, it is convenient to use an elliptic projection of the solution of
into the finite element space M h . Let ~ defined by
Then it follows from standard results of the finite element method [15], [30], [37] that
(5.3a)
1. The same result applies to the time-
differentiated forms of (5.1) [40]:
@t
@t
@s
@t
As for the analysis of the mixed finite element method, we use the the following
two projections instead of the elliptic projections introduced in [14] and [18]. So the
present analysis is different from and in fact simpler than those in [14] and [18]. Each
of our mixed finite element spaces [4]-[6], [11], [29], [34] has the property that there
are projection operators \Pi
such that
kr
and (see, e.g., [9], [19])
(r
Note that, by (3.3) and (4.7),
Finally, we prove some bounds of the projections ~ s and ~
p. Let be the
interpolant of s in M h . Then we see, by (4.1), (5.3b), the approximation property of
s, and an inverse inequality in M h , that
ks
ks
ks \Gamma sk 0;1
ks
where fl is given as in (5.3b). This implies that k~sk 1;1 is bounded for sufficiently
smooth solutions since k - 1. The same argument applies to k@~s=@tk 1;1 . Next, note
that, by the approximation property of the projection P h [27],
These bounds on ~
are used below.
We are now ready to prove some results. Below " is a generic positive constant
as small as we please.
5.1. Analysis of the mixed method. We first analyze the mixed method
(4.5). We set The following error
equation is obtained by subtracting (4.5) from (3.1) at applying (5.8) and
(ff(s
@t
@t
Below C i indicates a generic constant with the given dependencies.
Lemma 5.1. Let (u; p) and solve (3.1) and (4.5), respectively. Then
ks
@t
@t
@t
Proof. Set add the resulting equations at
use (3.3), (4.7), and (5.12) to see that
where
@t
@t
Then (5.15) can be easily seen.
Lemma 5.2. Let (u; p) and and (4.5), respectively. Then
ae
@t
ks
\Theta k@'
@t

\Psi \Deltat n
oe
@t
@t
@t
Proof. Difference equations (5.13) and (5.14) with respect to n, set
in the resulting equations, divide by \Deltat n , and add to obtain
where
\Deltat n
\Deltat n
\Deltat n
\Deltat n
@t
@t
@t
@t
\Deltat n
\Deltat n
\Deltat n
\Deltat n
\Deltat n
Observe that the left-hand side of (5.16) is larger than the quantity
(5.17)2\Deltat n
where
We estimate the new term T n
2 in detail. Other terms can be bounded by a simpler
argument. To estimate T n
2 , we write
\Gamma\Phi
\Deltat n
\Deltat n
\Deltat n
Note that
[A(s
@s
(bs
@A
@s
(bs
@A
where
and similar inequalities hold for b
s
Consequently, with -
see that
@s
(bs
@p@s
(p
\Deltat
(bp
\Deltat
so that
ks n\Gamma2
where
and an analogous inequality holds for b s
Also, we see that
[A(s
\Phi @A
@s
(p
which implies that
ks
Next, it can be easily seen that
ks
Finally, since
we find that
Hence T n
2 can be bounded in terms of T n
Other terms are bounded as
follows:
ks n\Gamma2
ks n\Gamma2
ks
ks
@t
@t
ks
ks n\Gamma2
can be bounded as in (5.19), e.g.,
kbs
ks
ks
ks
ks
apply these inequalities and (5.17)-(5.20), multiply (5.16) by \Deltat n , sum n, and
properly arrange terms to complete the proof of the lemma.
The error equations (5.13) and (5.14) are usually exploited to derive error estimates
in the parabolic mixed finite element method [18], [27]. To handle the difficulty
arising from the combination of the Dirichlet boundary condition (1.3) and the non-linearity
of the differential system (2.9)-(2.11), we must use their time-differentiated
forms, as mentioned before. Also, the three terms T n
care of the quadratic
terms in the velocities, which require more regularity on u than those without
these quadratic terms, as seen from Lemma 5.2.
5.2. Analysis of the saturation equation. We now turn to analyzing the
finite element method (4.6).
Lemma 5.3. Let s and s h solve (3.2) and (4.6), respectively. Then
ae
ks
@t
ks
0;1 \Deltat n
oe
@t
@t
Proof. Subtract (4.6) from (3.2) at use (5.1) at set the test
function to see that
where
@t
@t
The left-hand side of (5.21) is bigger than the quantity
\Gamma2\Deltat n (D(s n\Gamma2
defined by
and is bounded by
Next, it can be easily seen that
@t
To avoid an apparent loss of a factor h in B n
use summation by
parts on these items. We work on B n
3 in detail, and other quantities can be estimated
similarly. Applying summation by parts in n and the fact that - we see that
\Psi
so that, using the same argument as for (5.18),
3 \Deltat n
ks
ks
oe
where kbs n
can be estimated as in (5.20). The term
7 \Deltat n has the
same bound as in (5.25). Also, we find that
4 \Deltat n
oe
and
5 \Deltat n
ks
ks
oe
multiply (5.21) by \Deltat n , sum n, and use (5.22)-(5.27) to complete the proof of
the lemma.
5.3. estimates. We now prove the main result in this section. Define
K2Thp
@t
K2Thp
@t
@s
@t
Theorem 5.4. Let (u; p; s) and
respectively. Then, if the parameters \Deltat, h p , and h satisfy
we have
ks
@t
@t
Proof. Take a 1)-multiple of the inequality in Lemma 5.3, add the resulting
inequality and the inequality in Lemma 5.2, and use (5.3)-(5.7), (5.15), and the
extension of the solution for t - 0 to obtain
ae
\Theta (k@-
oe
where In deriving (5.29), we required that the " appearing in Lemma
5.3 be sufficiently small that (C 1 increases C 2 , but not C 1 . Observe
that, by (5.12),
The same result holds for - fl and fi fl . Combine (5.29), (5.30), and an inverse inequality
to see that
ae
\Theta (k@-
oe
We now make the induction hypothesis that
. Note that, by (5.12), (5.32) holds trivially for
(5.32), (5.31) becomes
ae
\Theta (k@-
oe
Using (5.28), we choose the discretization parameters so small that
Then it follows from (5.33) that
ae
oe
which, together with Gronwall's inequality, implies that
where
for \Deltat not too large. Consequently, the induction argument is completed and the
theorem follows.
We remark that, if h and h p are of the same order as they tend to zero, then
since k    - k   + 1. Since k - 1,
3:
Also, if k    - 2, we see that
3:
Thus, for (5.28) to be satisfied, we assume that k    - 2. This excludes the mixed
finite element spaces of lowest order, i.e., k   1. The lowest order case has to be
treated using different techniques. If the nonlinear coefficients ff(s) and c(s; p) in (4.5)
are projected into the finite element space W h , the technique developed in [10] can be
used to handle the lowest order case. We shall not pursue this here.
5.4. estimates. The main objective of this paper is to establish the
estimates given in Theorem 5.4. For completeness, we end this section with
a statement of L 1 -estimates for the errors in the two-dimensional
case.
Theorem 5.5. Assume that (p; s) and (p h ; s h ) satisfy (3.1), (3.2) and (4.5),
(4.6), respectively, and the parameters h p and h satisfy (5.28). Then
(\Omega\Gamma/
ks
Proof. First, it follows from the approximation property of the projection P h [27]
that
Also, from [27, Lemma 1.2] and (5.13), we see that
so that, by Theorem 5.4,
This, together with (5.37), implies (5.35). Finally, apply the embedding inequality
(5.3b), and (5.34) to obtain (5.36).
6. Finite elements for a degenerate problem. In this section we consider
a degenerate case where the diffusion coefficient D(s) can be zero. Since the pressure
equation is the same as before, we here focus on the saturation equation. For simplicity
we neglect gravity. Then the saturation equation (2.11) can be written as
@s
@t
@t
2\Omega \Theta J:
For technical reasons we only consider the Neumann boundary condition (2.15):
@\Omega \Theta J;
and the initial condition is given by
We impose the following conditions on the degeneracy
of D(s):
where the fi i are positive constants and ff j and - satisfy the conditions:
2:
Difficulties arise when trying to derive error estimates for the approximate solution
of (6.1) and (6.2) with D(s) satisfying the condition (6.3). To get around this
problem, we consider the perturbed diffusion coefficient D - (s) defined by [13], [24],
[35], [38]
g. Since the coefficient D - (s) is bounded away from
zero, the previous error analysis applies to the perturbed problem:
OE
@t
@t
2\Omega \Theta J;
(6.4a)
(D
@\Omega \Theta J;
(6.4c)
We now state a result on the convergence of s - to s as - tends to zero. Its proof
is given in [24] for the case where dw j 0 and the right-hand side of (6.1) is zero, and
can be easily extended to the present case.
Theorem 6.1. Assume that D(s) satisfies (6.3) and there is a constant C   ? 0
such that
where
Z sD(-)d-:
Then there is C independent of -, s, and - such that
As shown in [24], the requirement (6.5) is reasonable. We now consider a fully-
discrete finite element method for (6.4). Let M h be the standard C 0 piecewise linear
polynomial space associated with T h ; due to the roughness of the solution to (6.1) and
(6.2), no improvements in the asymptotic convergence rates result from taking higher
order finite element spaces. Also, we extend the domain of D - and q w as follows:
ae D - (1) if -
and
Now the finite element solution s n
to (6.4) is given by
(6.7a)
where P h is the L 2 -projection onto M h . The following theorem states the convergence
of s h to s. For (6.8) below to be satisfied, we see from (6.6) that the perturbation
parameter - need to satisfy the relation
Theorem 6.2. Let s and s h solve (6.1), (6.2) and (6.7), respectively, and let the
hypotheses of Theorem 6.1 be satisfied. Then there is C independent of -, s, and -
such that
The proof can be carried out as in [25], [35], and [38]; we omit the details.
7. Simulation with various boundary conditions. Let
@\Omega be a set of four
disjoint regions As
mentioned in the introduction, the most commonly encountered boundary conditions
for the two-pressure equations are of first-type, second-type, third-type, and well type.
Then we consider for
(7.3a)
(7.4a)
are given functions, d j is an arbitrary scaling constant,
and - is the outer unit normal to @ Note that \Gamma 1 is of the first type, \Gamma 2 is of the
third type (it reduces to the second type as - ff j 0), \Gamma 3 is of the well type, and on \Gamma 4
we have the Dirichlet condition for the air phase and the Neumann condition for the
water phase. Let \Gamma Then the
global boundary conditions for the pressure-saturation equations (2.9)-(2.11) become
(7.7a)
where pD and s D are the transforms of pwD and paD by (2.2) and (2.3), and
Z pc (s)q a
c (-)
Z pc (s)q a
c (-)
Z pc (s)q w
c (-)
We now incorporate the boundary conditions (7.5)-(7.10) in the finite element
scheme given in (4.5) and (4.6). The constraint V h ae V says that the normal components
of the members of V h are continuous across the interior boundaries in T hp .
Following [2], [9], we relax this constraint on V h by introducing Lagrange multipliers
over interior boundaries. Since the mixed space V h is finite dimensional and defined
locally on each element K in T hp , let V h . Then we define
~
for each K 2 T hp
ae
e
for each j;
oe
and W h and M h are given as before. The mixed finite element solution of the pressure
equation is fu n
(c(s
(ff(s
and the finite element method for the saturation is given for s n
satisfying
h )rs n
@t
. The computation of these equations can be carried out as in
(4.5) and (4.6). Note that the last equation in the unconstrained mixed formulation
above enforces the continuity requirement on u h , so in fact . It is well known
[2], [9] that the linear system arising from this unconstrained mixed formulation leads
to a symmetric, positive definite system for the Lagrange multipliers, which can be
easily solved. Also, the introduction of the Lagrange multipliers makes it easier to
incorporate the boundary conditions (7.5)-(7.10).
We now present a numerical example. The relative permeability functions are
taken as follows:
where s rw and s ra are the irreducible saturations of the water and air phases, respec-
tively. The capillary pressure function is of the form
where fl and \Theta are functions of the irreducible saturations. The water and air viscosities
and densities are set to be 1cP and 0:8cP , and 100kg=m 3 and 1:3kg=m 3 ,
respectively. The permeability rate is two-dimensional domain of 4m
width by 1m depth is simulated. Finally, the boundary of the domain is divided into
the following segments:
A uniform partition
of\Omega into rectangles with \Deltay is taken, and the time
step \Deltat is required to satisfy (5.28). The Raviart-Thomas space of lowest-order over
rectangles is chosen. Tables 1 and 2 describe the errors and convergence orders for the
pressure and saturation at time respectively. Experiments at other times
and on finer meshes are also carried out; similar results are observed and not reported
here.

Table

1. Convergence of p h at

Table

2. Convergence of s h at
From

Table

1, we see that the scheme is first-order accurate both in L 2 and L 1
norms for the pressure, i.e., optimal order. Table 2 shows that the scheme is almost
optimal order for the saturation. Thus the numerical experiments in the two tables
are in agreement with our earlier analytic results.



--R

On the solvability of boundary value problems for degenerate two-phase porous flow equations
Mixed and nonconforming finite element methods: implementation
Dynamics of Fluids in Porous Media



one dimensional simulation and air phase velocities
Mathematical Models and Finite Elements for Reservoir Simula- tion
Analysis of mixed methods using conforming and nonconforming finite element methods


Multiphase flow simulation with various boundary conditions

Mixed finite element methods for compressible miscible displacement in porous media
The Finite Element Method for Elliptic Problems

The approximation of the pressure by a mixed method in the simulation of miscible displacement


Characteristic Petrov-Galerkin subdomain methods for two phase immiscible flow
Efficient time-stepping methods for miscible displacement problems in porous media
Galerkin methods for miscible displacement problems with point sources and sinks-unit mobility ratio case
Timestepping along characteristics for a mixed finite element approximation for compressible flow of contamination from nuclear waste in porous media
A priori estimates and regularization for a class of porous medium equations

Fundamentals of Soil Physics
estimates for some mixed finite element methods for parabolic type problems

Mixed finite elements in

Fundamentals of Numerical Reservoir Simulation
On the transport-diffusion algorithm and its application to the Navier-Stokes equations
An implicit diffusive numerical procedure for a slightly compressible miscible displacement problem in porous media
A mixed finite element method for second order elliptic problems
Numerical Methods for flow through porous media I
Maximum norm stability and error estimates in parabolic finite element equations
Optimal L 1 estimates for the finite element method on irregular meshes
A near optimal order approximation to a class of two sided nonlinear degenerate parabolic partial differential equations
Galerkin Finite Element Methods for Parabolic Problems
A priori L 2 error estimates for Galerkin approximation to parabolic partial differential equations
--TR

--CTR
M. Afif , B. Amaziane, On convergence of finite volume schemes for one-dimensional two-phase flow in porous media, Journal of Computational and Applied Mathematics, v.145 n.1, p.31-48, 1 August 2002
E. Abreu , J. Douglas, Jr. , F. Furtado , D. Marchesin , F. Pereira, Three-phase immiscible displacement in heterogeneous petroleum reservoirs, Mathematics and Computers in Simulation, v.73 n.1, p.2-20, 6 November 2006
Z. Chen , R. E. Ewing, Degenerate Two-Phase Incompressible Flow IV: Local Refinement and Domain Decomposition, Journal of Scientific Computing, v.18 n.3, p.329-360, June
error estimate;air-water system;finite element;porous media;numerical experiments;mixed method;compressible flow;time discretization