A mixed nite element discretization on non-matching multiblock grids ...

East-West J. Numer. Math., Vol.0, No. 0, pp.1{21 (1997)

c VSP 1997

A mixed nite element discretization on non-matching multiblock grids for a degenerate parabolic equation arising in porous media ow IVAN YOTOV

Abstract | Mixed nite element methods on multiblock domains are considered for nonlinear degenerate parabolic equations arising in modeling multiphase ow in porous media. The subdomain grids need not match on the interfaces, where mortar nite element spaces are introduced to properly impose ux-matching conditions. The low regularity of the solution is treated through time integration, and the degeneracy of the diusion is handled analitically via the Kirchho transformation. With an appropriate choice of the mortar spaces, the error for both a semidiscrete (continuous time) scheme and a fully discrete (backward Euler) scheme is bounded entirely by approximation error terms of optimal order. Keywords. Mixed nite element, degenerate parabolic equation, nonlinear, mortar nite element, multiblock, non-matching grids, error estimates, porous media AMS subject classi cations. 65M60, 65M12, 65M15, 35K65, 76S05 Multiblock nite elementtechniques on non-matching grids have become increasingly popular in recent years. They combine the exibility of modeling irregularly shaped domains with the convenience of constructing the grids locally. In porous media applications they also allow accurate modeling of large scale geological features such as faults, layers, and fractures. We de ne a multiblock domain to be a simply connected domain 2 Rd, d = 2 or 3, that is a union of non-overlapping subdomains or blocks i, i = 1; :::; k. For the purpose of the analysis we assume that each block i is convex. In the numerical modeling of multiphase ow in porous media the coupled nonlinear system of conservation equations is often written as an equation of elliptic or parabolic type for some reference pressure and several saturation equations of advection-diusion type [7, 13]. The diusion in the saturation equations is degenerate (is zero at extreme saturation values). This causes a very low regularity of the solution and imposes diculties for the numerical method. In this paper we study mixed nite element discretizations for the saturation equation on multiblock domains. The local mass conservation property of the mixed methods is particularly important in modeling porous media ow. We allow the subdomain grids to be non-matching along the interfaces. To properly impose ux-matching conditions we Texas Institute for Computational and Applied Mathematics, The University of Texas at Austin, Austin, TX 78712; [email protected].

2

Ivan Yotov

introduce mortar nite element spaces along the subdomain boundaries. Mortar spaces have been a common tool in standard and spectral nite element methods (see e.g. [9]). Recently mortar [18, 29, 4] and non-mortar [6] techniques have been successfully applied in the context of mixed nite element methods on multiblock grids. Due to the degeneracy in the diusion the solutions can have compact support, thus behaving very dierently form the solutions of non-degenerate parabolic problems. In addition, certain error terms cannot be directly bounded in the analysis. Many authors introduce a regularized problem and then approximate it [27, 17, 24, 15, 14]. The numerical error is then a sum of the discretization error for the regularized problem and the dierence between the solutions to the regularized and the original problems. Another common technique is to handle the degeneracy analytically via the Kirchho transformation [27, 15, 14, 5, 28], which is the approach we take. We combine the use of mortar nite elements along subdomain interfaces with techniques from [5], where mixed discretizations for similar degenerate parabolic equations on a single block have been studied. We consider a continuous time scheme and a fully discrete backward Euler scheme and bound the discretization error in the two versions entirely by approximation error terms of optimal order. Critical in the analysis is the choice of mortar spaces along the non-matching interfaces. The mortar elements need to provide one order higher approximation than the traces of the subdomain velocity spaces (see Remark 2.2). At the same time, the mortar spaces should be controled by the subdomain spaces (see hypothesis (H1) and Remark 3.1). Our analysis improve the single block results from [5] in several ways. We avoid the assumption from [5] on smoothness of @ by requiring a minimal smoothness for the ux (see Lemma 2.2.). Also, we weaken the assumption ([5], A2) on the advective and the source terms and only assume physically reasonable behavior (see (2.7), (2.8)). Finally, by choosing a dierent ux test function in (3.4), we avoid a non-standard approximation error term involving the divergence of the dierence of two discrete projections of the ux, which appears in the analysis in [5]. The rest of the paper is organized as follows. In the next section we illustrate how a degenerate parabolic saturation equation arises in the fractional ow formulation of two phase incompressible ow. In Section 2 we formulate the mortar mixed method for the saturation equation. The analysis of the continuous in time and the fully discrete schemes is presented in Sections 3 and 4, respectively.

1. FRACTIONAL FLOW FORMULATION FOR TWO PHASE FLOW

Two phase immiscible ow in porous media is modeled by the system of conservation equations [7, 13] @ ('sii) + r u = q in (0; T ]; (1.1) i i i @t ui = ? ki(si)K (rpi ? i grD) in (0; T ]; (1.2) i

i = w (wetting), n (non-wetting), coupled with sw + sn = 1;

(1.3)

A multiblock mixed method for a degenerate parabolic equation

3

pc (sw ) = pn ? pw ; (1.4) where si is the phase saturation, i is the phase density, ' is the porosity, qi is the source term, ui is the Darcy velocity, pi is the phase pressure, K is the absolute permeability tensor, ki(si) is the phase relative permeability, i is the phase viscosity, g is the gravitational constant, D is the depth, and T is the nal time. We start by reformulating the problem in a standard way by writing it in a fractional

ow form (pressure and saturation equation). Let i = ki ; i = w; n; i denote the phase mobilities, and let = w + n be the total mobility. Let u = u w + un be the total velocity. For simplicity of the presentation, we assume incompressible ow and medium (constant i and ') and neglect gravity eects. Multiplying equations (1.1) by 1=i and adding them together, we get r u = q; (1.5) where q = qw =w + qn=n . Let s = sw , and de ne the global pressure [13] to be Z p (s) n ! ?1 ( ) d: p p = pw + c 0 Thus u = ?K rp: (1.6) Equation (1.5), coupled with (1.6), is referred to as the pressure equation. Since > 0 and K is a symmetric positive de nite tensor, this is an elliptic equation. For compressible ow the pressure equation is parabolic. To derive the saturation equation, we rst observe that w u = u ? w n K rp (s): w c Substituting this expression into the water conservation equation (1.1), we get the saturation equation ' @s (1.7) @t + r ( (s)u + (s)K rpc (s)) = q~w ; where (s) = w =, (s) = w n=, and q~w = qw =w . Note that pc (s) is a strictly monotone decreasing function; therefore with c (1.8) (s) = ?(s) @p @s c

4

Ivan Yotov

equation (1.7) takes the advection-diusion form

' @s @t + r ( (s)u ? (s)K rs) = q~w :

(1.9)

The diusion term vanishes at s = 0; 1, the minimum and maximum saturation values. This is due to the behavior of the relative permeability and the capillary pressure functions (see, e.g., [8]). This double degeneracy is the main source of diculties in the numerical approximation. The solutions of degenerate parabolic equations have very low regularity. It has been shown that (see [25, 16, 19, 2, 1, 3])

s 2 L1(0; T ; L1( )); @s 2 L2(0; T ; H ?1( )): @t Since the saturation s satis es 0 s(x; t) 1, (x; t) 2 [0; T ], we also have s 2 L1 (0; T ; L1( )):

(1.10) (1.11) (1.12)

We handle the degeneracy in the diusion analytically using the Kirchho transformation [27, 15, 14, 5]. Let Zs D(s) = ( ) d: 0

Then

rD(s) = (s)rs;

and (1.9) becomes

' @s @t + r ( (s)u ? K rD(s)) = q~w :

(1.13)

2. A MORTAR MIXED METHOD FOR THE SATURATION EQUATION

In this section we present a multiblock variational formulation of (1.13) that respects the low regularity of the solution, and a mixed nite element discretization using mortar elements on the interfaces. For the rest of the paper we omit the porosity ', which is assumed constant and does not eect the error analysis. The saturation s(x; t) satis es

@s + r ( (s)u ? K rD(s)) = q~ (s) w @t ( (s)u ? K rD(s)) = 0 s(x; 0) = s0(x)

in (0; T ]; on @ [0; T ]; in ;

(2.1) (2.2) (2.3)

where is the outward unit normal vector to @ . For simplicity we assume that no ow boundary conditions are imposed, although more general boundary conditions can also be


5

treated. We need several assumptions on the coecients of the above equation. We rst assume that 8 > 0 s 1 ; < 1jsj1 ; s ; (2.4) (s) > 2; : 3j1 ? sj2 ; 12 s 1;2 where i, 1 i 3, are positive constants, and i and i, i = 1; 2, satisfy 0 < 1 < 1=2 < 2 < 1; 0 < i 2: Note that (2.4) controls the rate of degeneracy of the diusion. We also assume that there exists a positive constant C such that kD(s1 ) ? D(s2)k20 C (D(s1) ? D(s2 ); s1 ? s2); for s1; s2 2 L2( ): (2.5) Here and for the reminder of the paper (; )S and h; i@S denote the L2-inner product (or duality pairing) on S Rd and @S 2 Rd?1, respectively, and k k0;S denotes the L2-norm on S . We omit S if S = . A sucient condition for (2.5) is 0 @D (2.6) @s (x; t; s) C for (x; t) 2 [0; T ]; 0 s 1: Bounds (2.4) and (2.6) follow from the physical behavior of the relative permeabilities and the capillary pressure [3, 15]. Typical relations (see, e.g., [8, 20]) are w (s) s2; pc (s) s?1 as s ! 0; o(s) (1 ? s)2; pc (s) (1 ? s)2 as s ! 1; 0 < 1 < 1; 0 < 2 < 1: Therefore, with (1.8), (s) s1?1 as s ! 0 and (s) (1 ? s)1+2 as s ! 1; which implies both (2.4) and (2.6). Finally, we assume that, for 0 s1; s2 1, j (s1) ? (s2)j2 C (D(s1) ? D(s2 ))(s1 ? s2); (2.7) 2 jq~w (s1) ? q~w (s2)j C (D(s1) ? D(s2))(s1 ? s2): (2.8) Bounds (2.7) and (2.8) are justi ed by the following lemma, proven in [15]. Lemma 2.1. Suppose satis es (2.4). If f

2 C 1[0; 1] and f 0(0) = f 0(1) = 0 with f 0

Lipschitz at 0 and 1, then there exists a positive constant C such that jf (a) ? f (b)j2 C (D(a) ? D(b))(a ? b) for 0 a; b 1:

Note that the fractional ow function (s) satis es the conditions of Lemma 2.1.. The well term q~w (s) satis es the conditions of Lemma 2.1. at the injection wells. At the production wells, q~w(s) kw (s), so q~w0 (0) = 0. Therefore (2.8) holds, if s s < 1 at the production well, which covers all cases of physical interest.

6

Ivan Yotov

Remark 2.1. The fractional ow function (s) and the integrated diusion function

D(s) are both S -shaped with zero derivatives at the end points. Bound (2.7) relates the

rates of degeneracy of the derivatives of the two functions and indicates, in a sense, that the diusion dominates the advection.

In the standard mixed variational formulation, equation (2.1) is multiplied by a test function w 2 L2( ) and integrated in space. In our case however, because of (1.11), the integral ( @s @t ; w) is not well de ned. To avoid this problem we integrate (2.1) in time from 0 to t [23, 5] to obtain the equivalent equation Zt Zt (2.9) d = q~w (s) d + s0(x); in [0; T ]; s(x; t) + r 0

0

where

= (s)u ? K rD(s): (2.10) Before presenting the variational formulation, we give the following regularity result. Lemma 2.2. Assume that there exists some 0 < " < 1=2 such that, for 1 i k ,

Zt 0

Then, for every t 2 [0; T ],

Zt 0

(s)u d 2 L2(0; T ; H "( i)):

(2.11)

d 2 (H " ( i))d \ H (div; i):

(2.12)

Proof. The argument is similar to one from [5]. The assumption (2.11) is reasonable, since the fractional ow function (s) has zero derivatives at the degeneracy values s = 0; 1. It follows from (2.9) and (1.12) that Zt d 2 H (div; i): 0

Using (2.10) we have in i

Zt

?r K r 0 D(s) d = r

Zt 0

d ? r

Zt 0

(s)u d:

By (2.11) and elliptic regularity, R0t D(s) d 2 H 1+" ( i), which, along with (2.10) and (2.11) implies that Zt d 2 (H "( i ))d: 2 0 Let ?i;j = @ i \ @ j , ? = [1i 0, we have Zt Zt 2 jT1j C 0 ks^ ? sk0 d + 0 kD(s) ? D(sh )k20 d; Z tZ Zt jT3j C 0 0 (s ? sh ; D(s) ? D(sh )) d d + 0 kD(s) ? D(sh )k20 d; Zt Zt jT4j 0 (s ? sh; D(s) ? D(sh )) d + C 0 kk20 d; Zt Zt jT5j C h?1 k ? Ph k20;? d + kk20 d ; 0

0

11


using (2.8) for the bound of T3, (2.7) for the bound of T4, and Lemma 3.1. for the bound of T5. To estimate T2 we integrate by parts in time: Z Zt Z Z ! @ a( ? h) d; @ T2 = ? d d d ? 0 0 0 0 Z t Zt Zt d ; d ? + a( ? h) d; 0

0

0

therefore,

8Z Z

2 Z t

@ Z Z

2 < t

jT2j C : 0 0 ( ? h ) d d + 0

@ 0 d ? 0 d

d 0 0

2

Z t Z t

2)

Z t d

+

( ? h) d

d ? +

0 0 0 0 0

The term T6 is the most dicult to bound. Integration by parts in time gives + Z Z k Z t * @ Z X T6 = ? d i; (Ph ? h ) d d d ? 0 0 i=1 0 @ 0 ? Z Z Z k t t t X d i d; (Ph ? h ) d : d ? + i

0

0

i=1

0

?i

Therefore,

8Z Z Z

2 < t

@ jT6j C : 0

@ 0 d ? 0 d

d 0;?

2

Z t Z t

2 ) Z t

Z

d d + (Ph ? h) d

d d ? + 0 0 0 0 0;? 0;?

Z t

2 +

(Ph ? h ) d

: 0 0;?

To bound the last two terms, we consider, for 1 i n and any xed t 2 (0; T ], the auxiliary problem

' ? ' = 0; Zt r' i = Qh;i(Ph ? h ) d; 0

in i ;

(3.11)

on @ i:

(3.12)

Note that the Neumann boundary data is in H 1=2?"( i) for any " > 0 and by elliptic regularity

Z t

; 1 m 2: (3.13) k'km?"; C 0 Qh;i(Ph ? h ) d

m?"?3=2;@

i

i

12

Ivan Yotov

We now integrate in time from 0 to t and take v = r' in (3.4) to obtain

Z t

2

Qh;i(Ph ? h ) d

0 0;@

Z t Z t d d = ? a( ? h) d; r' + (D(s) ? D(sh )) d; '

Z t 0 0 Z t + a( (s) ? (sh))u d; r' + Q ( P

?

) d; r ' h;i h i 0 0

@

= T7 + T8 + T9 + T10: We bound the terms on the right-hand side of (3.14) as follows. For any > 0,

Z t

2

jT7j C 0 ( ? h) d

+ (kr'k2"; + kr r'k20; )

2

20;

Z t

Z t

C 0 ( ? h) d + 0 Qh;i(Ph ? h) d

; 0;@

0;

i

i

i

i

i

i

(3.14)

i

i

i

i

using (3.3) for the rst inequality and (3.13) for the last bound. Note that in the above inequalities is an arbitrary small generic constant which may be dierent each time it appears. Similarly,

Z t

2 Zt jT8j C 0 kD(s) ? D(sh )k20; d +

0 Qh;i(Ph ? h ) d

; 0;@

2

Z t Zt

jT9j C 0 (s ? sh; D(s) ? D(sh )) d + 0 Qh;i(Ph ? h ) d

; 0;@

Z t

2 Zt jT10j C 0 kPh ? k20;@ d +

0 Qh;i(Ph ? h ) d

: 0;@

i

i

i

i

i

i

Combining together (3.14), the bounds on T7{T10, and (H1), we obtain

Z t

2

(Ph ? h ) d

0 ( Z t 0;? Z t 2 Zt C

0 d ? 0 d

0 + k(t)k20 + 0 kD(s) ? D(sh )k20 d Zt Zt + (s ? sh ; D(s) ? D(sh )) d + kPh ? k20;? d : 0

0

(3.15)

Combining (3.10), the bounds on T1{T6, (3.15), and using (2.5) and Gronwall's inequality, we arrive at the following result. Theorem 3.1. Assume that (2.4){(2.8) and (H1) hold. For the semi-discrete mixed -

nite element approximation (2.16){(2.18) of problem (2.1){(2.3), there exists a positive constant C such that, for every t 2 [0; T ],

2

Z t

(s ? sh; D(s) ? D(sh )) d + ( ? h) d

0 0 0

Zt

13


(Z t

Z t Z t

2

C 0 ks^ ? sk d + 0 d ? 0 d

0

Z t Z t

2 Zt d

d ? +h?1 k ? Ph k20;? d +

0 0 0 0;? Z t

@ Z Z

2 +

@ d

d d ? 0 0 0 0 Z

2 9 Z t

@ Z = + 0

@ 0 d ? 0 d

d ; : 2 0

0;?

Theorem 3.1. bounds the size of kD(s) ? D(sh )k0 by (2.5). It also allows us to derive a bound on ks ? shk?1 (see also [5]), where k k?1 is the H ?1( )-norm de ned by k k?1 = sup1 (k'; k') : 1 '2H ( ) Theorem 3.2. Assume that (2.4){(2.8) and (H1) hold. For the semi-discrete mixed nite element approximation (2.16){(2.18) of problem (2.1){(2.3), there exists a positive constant C such that, for every t 2 [0; T ],

ks(; t) ? (sh(; t)k2?1

Z t Z t

2 Zt 2 2 2

C h ks^ ? sk0 + 0 ks^ ? sk0 d + 0 d ? 0 d

0

Z t Z t

2

d

d ? h k d + 0 0 0 0;? Z t

@ Z Z

2 + 0

@ d ? 0 d

d 0 0 Z t

@ Z Z

2 9 = d

d ; : +

@ d ? 0 0 0 0;?

Zt +h?1 k ? P

2 0;?

Proof. For any ' 2 H01( ), we have (s ? sh ; ') = (s ? sh; ' ? '^) + (s ? sh; '^) = (s ? s^; ' ? '^) + (s ? sh; '^): By (3.5) we have that Zt k Zt X r 0 d ? 0 (s ? sh; '^) = ? i=1

Z t + 0 (~qw (s) ? q~w (sh)) d; '^ : h d ; '^

i

For the rst term on the right we write Zt Zt k Zt k Zt X X ? r d ? h d ; '^ = ? r d ? i=1

0

0

i

i=1

0

0

d ;' h

i

14

Ivan Yotov

X Zt Zt k Z t k Zt X d ? h d i; ' ? d ? = h d; r' 0 0 0 0 ?

i=1 i=1 Zt Zt Zt Zt k k X X d i ; ' d ? ? d ? = h d; r' 0 0 0 0 ?

i=1 i=1 k Z t X ( ? h) i d; ' ? Ph' ; ? 0 i

i

i

i

?i

i=1

using (3.6) for the last equality. Therefore j(s ? sh; '()j

Z t

Zt Zt Zt C hks ? s^k0 +

0 d ? 0 h d

+

0 d ? 0 d

0;? 0

Z t

Z t 1=2)

+ ( ? h) d + (s ? sh; D(s) ? D(sh )) d k'k1; 0 0 0 using that

k' ? '^k0 Chk'k1

for the rst term on the right, Lemma 3.1. and k' ? Ph'k0;? Ch1=2k'k1;

for the fourth term, and (2.8) for the last term. An application of Theorem 3.1. completes the proof. 2 i

i

4. DISCRETE TIME ERROR ANALYSIS

In this section we present the analysis for the fully discrete scheme (2.19){(2.21). The following error equations are obtained by subtracting (2.19){(2.21) from (2.13){(2.15) for t = tn . (a( n ? hn); v) = (D(sn ) ? D(snh ); r v)

v 2 Vh;i; (4.1) ? h n ? hn; v i i? + (a( (sn) ? (snh))un ; v) ; 1 1 0 0Z n t X j tj A ; wA (sn ? snh; w) + @r @ d ? h 0 j =1

1 0Z n t X @ = q~w (s) d ? q~w (sjh)tj ; wA + (s0 ? s0;h; w) ; w 2 Wh;i ; (4.2) i

i

i

i

n

i

i

n

k *Z t X i=1

0

0

n

i d ?

j =1

n X j =1

j j h i t ;

i

i

+

= 0; 2 Mh: (4.3) ? P (sn ) ? Dd (snh ), and = Ph n ? hn, d ? nj=1 hj tj , w = Dd

We now take v = n R0t to obtain (sn ? snh; Dd (sn ) ? Dd (snh )) + (a( n ? hn); n) n

i

15


1 0Z n t X (sn ) ? Dd (snh )A + (a( (sn) ? (snh))un; n ) = @ q~w (s) d ? q~w (sjh )tj ; Dd 0 j =1 ! + Zt k D k * Zt X X n E n n n n d ? d i; Ph ? h :(4.4) ? ? P h ; i ? + n

n

i=1

0

i=1

i

n

0

?i

We replace n by j in the above equation, multiply by tj and sum on j . The rst term on the left in (4.4) becomes n n X X (4.5) (sj ) ? Dd (sj ))tj = (sj ? sj ; D(sj ) ? D(sj ))tj ? T1; (sj ? sj ; Dd h

h

j =1

where

h

h

j =1

n X T1 = ? (sbj ? sj ; D(sj ) ? D(sjh ))tj j =1

To manipulate the terms involving n we rewrite it as ! X Zt Zt n n Zt X n j j ( j ? hj )tj : d ? t + = d ? d + ?1 n

j

n

0

0

t

j =1

j =1

j

(4.6)

Note that n is represented as a sum of an approximation error term, a time discretization error term, and an error term we are trying to bound. Using (4.6), the second term on the left in (4.4) becomes n X j =1

(a( j ? hj ); j )tj j X (a( j ? hj ); ( l ? hl )tl)tj ? T2 ? T3 j =1 l=1

2

n n

2 X

1=2 X 1 = 2 a ( j ? hj )tj

+ 21

a1=2( j ? hj )tj

0 ? T2 ? T3;

=

n X

j =1

where

j =1

0

n X

Zt

Z j ? j ); t h 0

(4.7)

!

d tj ; d ? a( 0 j =1 0 !1 j Zt n X X d ? ltl A tj : T3 = ? @a( j ? hj ); ?1

T2 = ?

j

j

l

j =1

l=1

t

l

For the second equality in (4.7) we used the well known identity, for any sequence fj g, 0 j 1 0 n 12 n n X X X @ j A + 2j = 2 @j X lA : j =1

j =1

j =1

l=1

16

Ivan Yotov

Combining (4.4), (4.5), and (4.7), we arrive at n X (sj ? sjh ; D(sj ) ? D(sjh ))tj j =1

2 n n 7

2 X X

1=2 X

1 + 2 a ( j ? hj )tj

+ 21

a1=2( j ? hj )tj

0 = Tm;

j =1

where

T4 = T5 =

0 n Zt X @

j =1 n X j =1

0

0

j

q~w (s) d ?

j X l=1

m=1

j =1

(4.8)

1 d q~w (slh )tl; Dd (sj ) ? D(sjh )A tj ;

(a( (sj ) ? (sjh))uj ; j )tj ;

k D n X X

E

j ? Ph j ; j i ? tj ; j =1 i=1 ! + Zt k * Zt n XX j j d i; Ph ? h tj : d ? T7 =

T6 = ?

i

j

j

j =1 i=1

0

0

?i

We next bound the terms Tm, m = 1; :::; 7. For any > 0 we have n n X X jT1j C ksbj ? sj k20tj + kD(sj ) ? D(sjh ))k20tj ; j =1 j =1 9 8 !

2 X > j j Zt n > = j
:

w w h h

> ; j =1 l=1 t ?1 l =1 0 n X + kD(sj ) ? D(sjh )k20tj ; l

l

j =1

n X j j j j j jT5j (s ? sh; D(s ) ? D(sh ))t + C k j k20tj ; j =1 j =1 n n o X jT6j C h?1k j ? Ph j k20;? + k j k20 tj ; j =1 n X

where we used (2.8) for the bound of T4, (2.7) for the bound of T5, and Lemma 3.1. for the bound of T6. To bound the rest of the terms we need the following discrete integration by parts identity. For sequences fAj gnj=0 and f j gnj=0, n nX ?1 X (4.9) Aj ( j+1 ? j ) = ? (Aj ? Aj?1) j + An n ? A0 0: j =1 j =0 We will apply (4.9) for Aj = Pjl=1 l, where fj gnj=1 is another sequence. In this case A0 = 0 and (4.9) leads to j nX ?1 X n n X X (4.10) j j = ? l( j+1 ? j ) + j n: j =1

j =1 l=1

j =1

17


R R To estimate T2, we apply (4.10) with j = a( j ? hj )tj and j = 0t d ? 0t d . 0j Z t +1 !1 Z t +1 nX ?1 X 1 T2 = @ a( l ? hl )tl; tj+1 d A tj+1 d ? t t j =1 l=1 1 0n Zt Zt X j ? @ a( j ? h)tj ; 0 d ? 0 d A j =1

j

2 nX ?1

X C1

a( l ? hl )tl

tj j =1 l=1 0 Z t +1 !

2 Z +1 nX ?1

1 t

+

tj+1 d

tj+1 d ? t t j =1 0

2 Z

n

2 Zt

t

X j j j

d ? d

; + a( ? h)t + C

0 0

0

j=1 0 where we assumed that there exists a constant C > 0 such that tj+1 C1tj ; 1 j N ? 1: (4.11) Similarly,

2

j nX ?1

X jT3j C1

a( l ? hl )tl

tj j =1 l=1 0 !

2 Z +1 nX ?1

1 t

tj+1 j +1 +

tj+1 d ? j +1 t

0 t j =1

2

n

2

Z n X

t

X j j j

j tj

; d ? + a( ? h)t + C

j=1

0 0 j =1 0 and

j

2 nX ?1

X jT7j

(Ph l ? hl )tl

tj j =1 l=1 0;? Z t +1 !

2 Z +1 nX ?1

1 t

+C

tj+1 d

tj+1 d ? t t j =1 0;?

n

2

!

2 Zt

Z t

X j j j

d

: d ? + (Ph ? h)t + C

0 0

j=1

0;? 0;? To further estimate the rst and the third terms, we consider, for 1 i k and any xed 1 n N , the auxiliary problem ' ? ' = 0; in i ; (4.12) n X on @ i; (4.13) r' i = Qh;i(Ph j ? hj )tj ; j

j

j

j

j

n

n

j

j

j

j

n

n

j

j

n

j

j

j

j

n

j =1

n

j

18

Ivan Yotov

Note that the Neumann boundary data is in H 1=2?"( i) for any " > 0 and by elliptic regularity

n

X j j j k'km?"; C

Qh;i(Ph ? h )t

; 1 m 2: (4.14) j =1 m?"?3=2;@

i

i

We now replace n by j in (4.1), multiply by tj , sum on 1 j n, and take v = r' to obtain

n

2

X

j

Qh;i(Ph j ? h )tj

j =1 0;@

1 1 0n 0n X X j j j j (s) ? Dd (sh ))t ; 'A = ? @ a( ? h)t ; r'A + @ (Dd j =1 j =1

1 *n 0n + X X j j j j j j j A @ Qh;i(Ph ? )t ; r' i a( (s ) ? (sh))u t ; r' + + i

i

i

j =1

i

j =1

@

= T8 + T9 + T10 + T11: We now have, for any small generic constant > 0,

n

X j j j

2 jT8j C

( ? h )t

+ kr'k20;

j =1

02;

n X C

( j ? hj )tj

+ (kr'k2; + kr r'k20;

j =1

20;

n

n

X

X j j j

C

( ? h )t

+

Qh;i(Ph j ? hj )tj

; j =1 j =1 0;@

0;

i

i

i

i

i

i

i

using (3.3) for the second inequality and (4.14) for the third inequality. Similarly,

n n X X

jT9j C kD(s) ? D(sh )k20; tj +

Qh;i(Ph j ? hj )tj

; j =1 j =1 0;@

n n X X

j jT10j C (s ? sh ; D(s) ? D(sh )) tj +

Qh;i(Ph j ? h )tj

; j =1 j =1 0;@

n n X X jT11j C kPh j ? j k20;@ tj +

Qh;i(Ph j ? hj )tj

: j =1 j =1 0;@

i

i

i

i

i

i

Combining (4.15), the bounds on T8{T11, and (H1), we conclude that 8

n

2

2 n

X > n < X

j ? j )tj C ( j ? j )tj

+ X kD(sj ) ? D(sj )k2 tj ( P

h

h

h h 0

> :

j=1 j =1 0;? 0 j =1

i

(4.15)

19


9 = + (sj ? sjh ; D(sj ) ? D(sjh ))tj + kPh j ? j k20;?tj ; : j =1 j =1 n X

n X

(4.16)

Putting together (4.8), the bounds on T1{T7, (4.16), manipulating j in T5 and T6 as in (4.6), applying (2.5) and the discrete Gronwall's inequality, we obtain the following estimate. Theorem 4.1. Assume that (2.4){(2.8) and (H1) hold, and that there exists a constant

C1 > 0 such that

tj+1 C1tj ; 1 j N ? 1: For the fully discrete mixed nite element approximation (2.19){(2.21) of (2.1){(2.3), there exists a positive constant C such that, for any 1 n N ,

2 n

n n X X

X

(sj ? sjh; D(sj ) ? D(sjh ))tj +

( j ? hj )tj

+ k( j ? hj )tj k20

0 j=1

j=1 j =1 8

!

2 j Zt n > < bj j 2

X X l )tl

q ~ ( s ) d ? q ~ ( s C >:ks ? s k0 +

w w

j =1 l=1 t ?1 0

!

2 ! 2 Z Z Z

1 t t t 1

+

tj ?1 d ? ?1 d

+

tj ?1 d ? j tj

t t t 0 9 0

! 2 Z Z =

t t +

1tj ?1 d ? ?1 d

+ h?1 kPh j ? j k20;?; tj : l

l

j

j

j

j

j

j

j

j

t

t

j

j

0;?

The following theorem can be shown as in the proof of Theorem 3.2.. Theorem 4.2. Under the assumptions of Theorem 4.1., there exists a positive constant

C such that, for any 1 n N ,

ksn ? snhk2?1 8

j Z !

2 n > < bj j 2

X t X l )tl

+ h?1 kP j ? j k2 q ~ ( s ) d ? q ~ ( s C >:ks ? s k0 +

h w w 0;?

j =1 l=1 t ?1 0

!

2

!

2 Zt Zt Zt 1 1

j +

tj ?1 d ? ?1 d

+

tj ?1 d ? j t

t t 0 0 9 t

! 2 Z Z =

t t +

1tj t ?1 d ? t ?1 d

; tj + Ch2kscn ? snk20: 0;? l

l

j

j

j

j

j

j

j

j

j

j

Remark 4.1. The estimates from Theorems 3.1.{4.2. bound the discretization error by optimal order approximation terms in time or space. The term h?1=2 k ? Ph k0;? provides approximation of order h1=2 higher then the other terms, assuming enough regularity of , since the functions in h are piece-wise polynomials of one degree higher than these in Vh .

20

Ivan Yotov

Acknowledgments. This work was partially supported by the United States Department

of Energy. The author thanks Professors Todd Arbogast and Mary F. Wheeler for useful comments and discussions.

REFERENCES 1. H. W. Alt and E. DiBenedetto, Nonsteady ow of water and oil through inhomogeneous porous media, Ann. Scuola Norm. Sup. Pisa Cl. Sci., 12 (1985), pp. 335{392. 2. H. W. Alt and S. Luckhaus, Quasilinear elliptic-parabolic dierential equations, Math. Z., 183 (1983), pp. 311{341. 3. T. Arbogast, The existence of weak solutions to single porosity and simple dual-porosity models of two-phase incompressible ow, Nonlinear Analysis, Theory, Methods and Applications, 19 (1992), pp. 1009{1031. 4. T. Arbogast, L. C. Cowsar, M. F. Wheeler, and I. Yotov, Mixed nite element methods on non-matching multiblock grids, TICAM Report 96-50, Texas Inst. Comp. Appl. Math., University of Texas at Austin, Nov. 1996. 5. T. Arbogast, M. F. Wheeler, and N. Zhang, A nonlinear mixed nite element method for a degenerate parabolic equation arising in ow in porous media, SIAM J. Numer. Anal., 33 (1996), pp. 1669{1687. 6. T. Arbogast and I. Yotov, A non-mortar mixed nite element method for elliptic problems on non-matching multiblock grids, in Symposium on Computational Mechanics Advances, Austin, Texas, 1997. To appear in Comput. Meth. Appl. Mech. Eng. 7. K. Aziz and A. Settari, Petroleum reservoir simulation, Applied Science Publishers, London and New York, 1979. 8. J. Bear, Dynamics of uids in porous media, Dover, New York, 1972. 9. C. Bernardi, Y. Maday, and A. T. Patera, A new nonconforming approach to domain decomposition: the mortar element method, in Nonlinear partial dierential equations and their applications, H. Brezis and J. L. Lions, eds., Longman Scienti c & Technical, UK, 1994. 10. F. Brezzi, J. Douglas, Jr., R. Duran, and M. Fortin, Mixed nite elements for second order elliptic problems in three variables, Numer. Math., 51 (1987), pp. 237{250. 11. F. Brezzi, J. Douglas, Jr., M. Fortin, and L. D. Marini, Ecient rectangular mixed nite elements in two and three space variables, RAIRO Model. Math. Anal. Numer., 21 (1987), pp. 581{ 604. 12. F. Brezzi, J. Douglas, Jr., and L. D. Marini, Two families of mixed elements for second order elliptic problems, Numer. Math., 88 (1985), pp. 217{235. 13. G. Chavent and J. Jaffre, Mathematical models and nite elements for reservoir simulation, NorthHolland, Amsterdam, 1986. 14. K. Fadimba and R. Sharpley, Galerkin nite element method for a class of porous medium equations. Preprint. , A priori estimates and regularization for a class of porous medium equations, Nonlinear World, 15. 2 (1995), pp. 13{41. 16. A. Friedman, The Stefan problem in several space variables, Trans. Amer. Math. Soc., 133 (1968), pp. 51{87. 17. J. Jerome and M. Rose, Error estimates for the multidimensional two-phase stefan problem, Math. Comp., 39 (1982), pp. 377{414. 18. Y. A. Kuznetsov and M. F. Wheeler, Optimal order substructuring preconditioners for mixed nite element methods on non-matching grids, East-West J. Numer. Math., 3 (1995), pp. 127{143.


21

19. O. Ladyzhenskaya, V. A. Solonnikov, and N. N. Ural'ceva, Linear and quasilinear equations of parabolic type, vol. 23 of Transl. Math. Monographs, Amer. Math. Soc., Providence, R.I., 1968. 20. L. W. Lake, Fundamentals of enhanced oil recovery, Prentice Hall, New Jersey, 1989. 21. T. P. Mathew, Domain decomposition and iterative re nement methods for mixed nite element discretizations of elliptic problems, PhD thesis, Courant Institute of Mathematical Sciences, New York University, 1989. Tech. Rep. 463. 22. J. C. Nedelec, Mixed nite elements in R3, Numer. Math., 35 (1980), pp. 315{341. 23. R. N. Nochetto, Error estimates for two-phase Stefan problems in several space variables, I: Linear boundary conditions, Calcolo, 22 (1985), pp. 457{499. 24. R. N. Nochetto, Error estimates for the multidimensional singular parabolic problems, Japan J. Appl. Math., 4 (1987), pp. 111{138. 25. O. A. Oleinik, A. S. Kalashnikov, and C. Yui-Lin, The Cauchy problem and boundary value problems for equations of type of nonstationary ltration, Izv. Akad. Nauk SSSR Ser Mat, 22 (1958), pp. 667{704. In Russian. 26. R. A. Raviart and J. M. Thomas, A mixed nite element method for 2nd order elliptic problems, in Mathematical Aspects of the Finite Element Method, Lecture Notes in Mathematics, vol. 606, SpringerVerlag, New York, 1977, pp. 292{315. 27. M. Rose, Numerical methods for ows through porous media. I, Math. Comp., 40 (1983), pp. 435{467. 28. C. S. Woodward and C. N. Dawson, Analysis of expanded mixed nite element methods for a nonlinear parabolic equation modeling ow into variably saturated porous media, TIACM Report 96-51, Texas Inst. Comp. Appl. Math.,University of Texas at Austin, Nov. 1996. 29. I. Yotov, Mixed nite element methods for ow in porous media, PhD thesis, Rice University, Houston, Texas, 1996. TR96-09, Dept. Comp. Appl. Math., Rice University and TICAM report 96-23, University of Texas at Austin.