Error estimates for the finite volume discretization for the porous medium equation

Journal of Computational and Applied Mathematics 234 (2010) 2135–2142

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Error estimates for the finite volume discretization for the porous medium equation I.S. Pop a,∗ , M. Sepúlveda b , F.A. Radu c,d , O.P. Vera Villagrán e a b c

CASA, TU Eindhoven, P.O. Box 513, 5600 MB Eindhoven, The Netherlands CI2 MA, Universidad de Concepción, Casilla 160-C, Concepción, Chile UFZ-Helmholtz Center for Environmental Research, Permoserstr. 15, D-04318 Leipzig, Germany

d

University of Jena, Wöllnitzerstr. 7, D-07749, Jena, Germany

e

Universidad del Bío-Bío, Collao 1202, Casilla 5-C, Concepción, Chile

article

info

Article history: Received 12 September 2008 Received in revised form 13 January 2009 Keywords: Degenerate parabolic equation Porous medium Finite volumes Error estimates

abstract We analyze the convergence of a numerical scheme for a class of degenerate parabolic problems modelling reactions in porous media, and involving a nonlinear, possibly vanishing diffusion. The scheme involves the Kirchhoff transformation of the regularized nonlinearity, as well as an Euler implicit time stepping and triangle based finite volumes. We prove the convergence of the approach by giving error estimates in terms of the discretization and regularization parameter. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Degenerate parabolic equations appear in the mathematical modeling of numerous real life processes. A well known example in this sense is the porous medium equation, describing the flow of an ideal gas in a homogeneous porous medium. More complex situations are encountered in petroleum reservoir and groundwater aquifer simulations. Compared to regular parabolic problems, and in particular to the heat equation, the diffusive term in the degenerate case depends on the unknown solution, and may vanish or blow up. Thus the parabolic character of the equation may change into an elliptic, or even hyperbolic one. The interfaces separating the domains of regularity – also called free boundaries – have finite speed of propagation. Generally, these are not known in advance and have to be determined together with the solution. Typically, the solutions of such problems are lacking regularity. Eventual singularities do not smooth out in time; these may even develop in time, giving the problem a strongly nonlinear character. Consequently, the numerical approximation of such solutions require adequate algorithms that are able to deal with both the free boundary, as well as the singularities of the solution. This paper is motivated by a combined mixed finite element (MFE) - finite volume (FV) scheme of a two phase flow model for the heap leaching of copper ores [1]. The convergence of such schemes has been investigated in [2] and [3], where the MFE method is employed for the flow component, whereas the saturation is discretized by using a FV scheme. The convergence results there, are obtained under several simplifying assumptions that rule out the degeneracy of the model. For the convergence of numerical schemes for degenerate parabolic problems, we refer to [4] for the finite element discretization, [5–8] for MFE schemes, [9] for a DG approach, and [10] for a multipoint flux approximation method. FV

∗

Corresponding address: TU Eindhoven, Mathematics and Computer Science, Den Dolech 2, 5612AZ, Eindhoven, The Netherlands. E-mail address: [email protected] (I.S. Pop).

0377-0427/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2009.08.071

2136

I.S. Pop et al. / Journal of Computational and Applied Mathematics 234 (2010) 2135–2142

schemes for porous media models are analyzed rigorously in [11,12], whereas a MFEM-FV approach is considered in [13]. There, convergence is obtained by compactness, and no estimates for the approximation error are given. Here, we consider the FV discretization of the degenerate parabolic equation

∂t u − ∆β(u) = r (u),

in QT ≡ (0, T ) × Ω .

(1)

Initially we have u(0) = u0 in Ω , whereas u = 0 on ∂ Ω . In the above 0 < T < ∞ is fixed, Ω is a bounded polygonal domain in R2 with a Lipschitz continuous boundary. The function β : R → R is non-decreasing and differentiable. Specifically, we assume the following: (A1) β is Lipschitz and differentiable, β(0) = 0, 0 ≤ β 0 (u) ≤ Lβ . (A2) u0 ∈ L2 (Ω ) and β(u0 ) ∈ H01 (Ω ). (A3) r : R → R is continuous in u; furthermore,

|r (u) − r (v)|2 ≤ C (u − v)(β(u) − β(v)) for any u, v ∈ R, where C > 0 does not depend on x, t , u and v . Moreover, r (0) = 0. By degeneracy we mean a vanishing diffusion, namely β 0 (u) = 0 for some u. An important example that can be written in the above form is the porous medium equation, where β(u) = um for some m > 1 whenever u ≥ 0, while r = 0. Another example is a simple model for melting and solidification, where β is increasing, piecewise linear, and vanishes on a certain interval, say [0, 1]. More complex is the Richards equation, which models unsaturated flow in porous media and involves nonzero convection. For ease of presentation, we restrict this paper to homogeneous Dirichlet boundary conditions; considering more general ones is straightforward. We use standard notations for function spaces, norms and scalar products: L2 (Ω ), H01 (Ω ), or its dual H −1 (Ω ). With X being one of the spaces above, L2 (0, T ; X ) extends the square integrability to time dependent functions. We let (·, ·) stand for the inner product on L2 (Ω ), or the duality pairing between H01 (Ω ) and H −1 (Ω ), k · k for the norm in L2 (Ω ), whereas k · kk denotes the norm in H k (Ω ). Moreover, we write u or u(t ) instead of u(t , x) and use C to denote a positive constant independent of the discretization parameters or the function itself. Saying this, we look for the weak solution of Problem P, solving: Problem WP. Find u ∈ H 1 (0, T ; H −1 (Ω )) such that β(u) ∈ L2 (0, T ; H01 (Ω )), and

(∂t u, ϕ) + (∇β(u), ∇ϕ) = (r (u), ϕ), for all ϕ ∈

H01

(2)

(Ω ) and for all t ∈ (0, T ], whereas u(0) = u0 in H

−1

(Ω ).

Existence and uniqueness for Problem WP is proven, e. g. in [14] and [15]. Notice that β(u) is more regular than u, this being a property that we exploit below. We also employ a regularization step in constructing the numerical scheme: with ε > 0 a given parameter, the nonlinear function β is approximated by βε satisfying β 0 ≥ ε . For simplicity, we consider the global perturbation

βε (u) ≡ β(u) + ε u.

(3)

Other perturbations can be considered and the analysis is similar. In general, βε should be invertible satisfying

ε ≤ βε0 (u) ≤ C ,

0 ≤ βε0 (u) − β 0 (u) ≤ ε,


0

and C ≤ βε−1 (βε (u)) ≤ 1/ε.

(4)

2. The time discretization Due to the lack of regularity, we employ the Euler implicit scheme to discretize the regularized Problem WP in time. This idea is used for constructing effective numerical algorithms; see e.g. [16], where compactness arguments are considered for showing the convergence in a general setting. Further, since β(u) is more regular than u, we first approximate θ ≈ β(u), and then u = βε−1 (θ ). With n ∈ N and τ = T /n denoting the (fixed) time step, we let tk = kτ . For k ∈ {1, . . . , n} we define the time discrete approximation θ k of β(u(tk )) as the solution of Problem WTk . Given θ k−1 ∈ H01 (Ω ), find θ k ∈ H01 (Ω ) such that for all ϕ ∈ H01 (Ω )

(βε−1 (θ k ) − βε−1 (θ k−1 ), ϕ) + τ (∇θ k , ∇ϕ) = τ (r (βε−1 (θ k )), ϕ).

(5)

The scheme is completed by the initial data. A straightforward choice is θ = βε (u ). Whenever u ∈ (Ω ), this gives θ 0 0 2 in the same space. However, in (A2) we have only required u ∈ L (Ω ). Following the discussion in [17], Chapter 3, one can consider θ 0 = β(u0 ) + ερµ ∗ u0 , where ρµ is a mollifier having a compact support in B(0, µ). With µ = O(ε), θ 0 is bounded uniformly in H 1 , whereas ku0 − βε−1 (θ 0 )k vanishes as ε & 0. It is worth mentioning that the convolution can be replaced by the solution of the heat equation at a (small) time, where the initial data is u0 . 0

0

0

H01

I.S. Pop et al. / Journal of Computational and Applied Mathematics 234 (2010) 2135–2142

2137

Remark 2.1. Inverting βε may be tedious. Further, since function calls increase the computing time significantly, for implementing the scheme we construct a look-up table of values for βε . Requiring little computer memory, as well as linear interpolation for the values not included in the table, this leads to a reduction of the computing time. Moreover, the monotonicity of βε allows a fast searching in this table, and thus the values of βε or its inverse can be obtained efficiently. Existence and uniqueness for the time discrete problems WPk is provided by standard results for nonlinear elliptic equations. Furthermore, assuming that the initial data is essentially bounded, the sequence of solutions θ k remain essentially bounded uniformly in k. The following estimates are obtained in [17], Chapter 3 (see also [18–20]). Lemma 2.1. Assume (A1)–(A3). With θ k solving the Problem WTk , for all 0 ≤ p ≤ n we have

kθ p k2 +

n X 

(βε−1 (θ k ) − βε−1 (θ k−1 ), θ k − θ k−1 ) + kθ k − θ k−1 k2 + τ k∇θ k k2 ≤ C .

(6)

k=1

Remark 2.2. The estimate (6) immediately implies n X

kβε−1 (θ k ) − βε−1 (θ k−1 )k2−1 ≤ C τ .

(7)

k=1

To prove the above, we use (5) and the Poincaré inequality to obtain

|(βε−1 (θ k ) − βε−1 (θ k−1 ), ϕ)| ≤ τ k∇θ k k + C kr (βε−1 (θ k ))k




for all ϕ ∈ H01 (Ω ) s.t. k∇ϕk = 1. Now 2.2 follows by (A3) and the above estimates. We use the following notations. Given a function f : QT → R integrable in time, define 1 f¯ k := τ

Z

kτ

(k−1)τ

f (s, ·)ds,

if k ≥ 1.

Further, f¯ 0 := f (0, ·). The errors are obtained in terms of eku and ekθ defined as k

ekθ := β(u) − θ k ,

eku := u¯ k − βε−1 (θ k ),

(8)

where k ≥ 0. Given a sequence {f ∈ (Ω ), k = 1, n}, the piecewise constant extension in time f∆ is defined as f∆ (t ) = θ k for t ∈ (tk−1 , tk ]. Further, G : H −1 (Ω ) → H01 (Ω ) stands for the Green operator defined by H01

k

(∇ Gψ, ∇ϕ) = (ψ, ϕ), for all ϕ ∈

H01

(Ω ), where ψ ∈ H

k∇ Gψk = kψk−1 ,

(9) −1

(Ω ). Therefore

kψk−1 ≤ C kψk,

(10)

where the last inequality applies only if ψ ∈ L (Ω ). We have ([17], Chapter 3): 2

Theorem 2.1. Assume (A1)–(A3). If u and θ k (k = 1, n) solve the problems given above, then eku 2−1

sup k k k=1,n

T

Z + 0

(βε (u(t )) − θ∆ (t ), u(t ) − βε−1 (θ∆ (t )))dt + kβ(u) − θ∆ k2L2 (Q ) ≤ C {τ + ε} . T

These estimates hold in a more  general setting, where convection terms can also be included. Using the results in [21], the estimates become optimal, C τ 2 + ε 2 . This holds in a simplified case (e.g. in the absence of convection and if β is maximal monotone having the range R). 3. The finite volume discretization Here we refer to the framework in [22] (see also [23,28,29] and [3]) and let T := {Ti , i ∈ I ⊂ N} be a regular and acute decomposition of Ω into triangles. We assume that the diameter of any triangle T ∈ T does not exceed h. Further, E and P stand for the set of triangle edges, respectively the set of nodes. Since Ω is assumed polygonal, such a decomposition is possible without introducing additional errors occurring when discretizing curved boundaries. In this case we also have E = Eint ∪ Eext , where Eext = E ∩ ∂ Ω and Eint = E \ Eext . In what follows we use the following notation:

|T |- the area of T ∈ T , |`|- the length of ` ∈ E , Ni - the triangles adjacent to Ti ∈ T , Ei - the edges of Ti ,

2138

I.S. Pop et al. / Journal of Computational and Applied Mathematics 234 (2010) 2135–2142

xi - the center of the circumcircle of Ti , `ij - the edge between Ti and Tj (where Tj ∈ Ni ), d(xi , `ij )- the distance from xi to `ij , and dij = d(xi , `ij ) + d(xj , `ij ) if `ij ∈ Eint , σij = |`ij |/dij - the ‘‘transmissibility’’ through `ij nij - the outward unit normal to `ij pointing into Tj (Tj ∈ Ni ). The assumptions on T ensure that xi ∈ int (Ti ) for all i. Furthermore, if Tj ∈ Ni , then the line through xi and xj is orthogonal to `ij . Given T we define the finite dimensional subspace of L2 (Ω ) Wh := {v ∈ L2 (Ω )/ v is constant on any T ∈ T },

(11)

which is spanned by the triangle indicator functions {χT / T ∈ T }. Furthermore we define Ph : L2 (Ω ) → Wh ,

(Ph w − w, wh ) = 0

(12)

for any wh ∈ Wh . With s ∈ {0, 1}, a constant C > 0 exists such that

kPh w − wk ≤ Chs kwks

(13)

for all w ∈ H (Ω ). Moreover, for any w ∈ L (Ω ) and Ti ∈ T we have s

2

w ¯ i := (Ph w)|Ti =

1

Z

|Ti |

w(x)dx.

(14)

Ti

As in the spatially continuous case, for any u, v ∈ L2 (Ω ) we define the discrete inner products

(u, v)h :=

X

|Ti |¯ui v¯ i ,

(u, v)1,h :=

X

σij (¯ui − u¯ j )(¯vi − v¯ j ),

(15)

`ij ∈E

Ti ∈T

where the value of u¯ and v¯ are extended by 0 outside Ω in view of the homogeneous Dirichlet boundary conditions. The associated discrete norms are denoted by k · kh , respectively k · k1,h . It is easy to see that for all u, v ∈ L2 (Ω ),

(u, v)h = (u, Ph v) = (Ph u, v).

(16)

Furthermore, in [22] the following discrete Poincaré inequality is given:

kuk ≤ C kuk1,h ,

(17)

for all u ∈ Wh , where C > 0 does not depend on h or u. Using (16) and (13), one obtains

|(u, v) − (u, v)h | = |(u − Ph u, v)| = |(u − Ph u, v − Ph v)| ≤ Chs+p kuks kvkp , for all u ∈ H (Ω ) and v ∈ H (Ω ), s, p ∈ {0, 1}. Further we define the discrete H s

kuk−1,h =

p

sup

wh ∈Wh ,kwh k1,h =1

−1

|(u, wh )h |.

(18)

norm (19)

Following [22] and [23] we write the finite volume scheme for the time discrete Problem WTk . With θh,i = θh |Ti and given

θhk−1 ∈ Wh , we seek θ k ∈ Wh such that for all Ti ∈ T it holds X |Ti |(βε−1 (θhk,i ) − βε−1 (θhk,−i 1 )) + τ σij (θhk,i − θhk,j ) = τ |Ti |r (βε−1 (θhk,i )).

(20)

`ij ∈Ei

To give a weak form of the above scheme, for any wh ∈ Wh we multiply (20) by wi = wh |Ti and sum up the resulting for all Ti ∈ T . Recalling the definitions in (15), after changing the summation order in the second term on the left we obtain the following Problem WDk . Given θhk−1 ∈ Wh , find θhk ∈ Wh such that for all wh ∈ Wh

(βε−1 (θhk ) − βε−1 (θhk−1 ), wh )h + τ (θhk , wh )1,h = τ (r (βε−1 (θhk )), wh )h . To complete the scheme, we define the initial data θ = 0 h

βε (Ph (βε−1 (θ 0 )))

(21)

∈ Wh .

The stability properties below are similar to the ones in the time discrete case. Lemma 3.1. Assume (A1)–(A3), and let θhk solving (21). For any 1 ≤ p ≤ n we have

kβε−1 (θhp )k2h +

n X

kβε−1 (θhk ) − βε−1 (θhk−1 )k2h + τ

k=1

n X k=1

kθhk k21,h ≤ C ,

n n X X (βε−1 (θhk ) − βε−1 (θhk−1 ), θhk − θhk−1 )h + kθhk − θhk−1 k2h ≤ C . k=1

k=1

(22)

I.S. Pop et al. / Journal of Computational and Applied Mathematics 234 (2010) 2135–2142

2139

Proof. We start by noticing that for any wh ∈ Wh , by (A3) and (4) we get 1

C kwh k21,h ≤ (βε−1 (wh ), wh )1,h ,

and

|(r (wh ), wh )h | ≤ C (wh , β(wh ))h2 kwh kh ≤ C kwh k2h .

(23)

βε−1 (θhk )

Next we take wh = into (21), sum the resulting up for k = 1, . . . , p and use the elementary identity 2a(a − b) = a2 − b2 + (a − b)2 and obtain 1

p h

kβε−1 (θ )k2h − kβε−1 (θh0 )k2h +

2

p X

! kβε−1 (θhk ) − βε−1 (θhk−1 )k2h

k=1

+τ

p p X X (θhk , βε−1 (θhk ))1,h = τ (r (βε−1 (θhk )), βε−1 (θhk ))h . k=1

k=1

By (23), the first part of the estimates is a direct consequence of the Gronwall lemma. Using the assumptions on β and r, testing (21) with wh = θhk − θhk−1 completes the proof.



Remark 3.1. As in the spatially continuous case, the estimates above immediately imply n X

kβε−1 (θhk ) − βε−1 (θhk−1 )k2−1,h ≤ C τ .

(24)

k=1

The fully discrete problems have unique solutions, as follows from: Theorem 3.1. Assume (A1)–(A3), and let θhk−1 ∈ Wh be given. Then the fully discrete problem WD k has a unique solution θhk , at least for moderately small time steps τ . Proof. We start with the uniqueness, which is a direct consequence of the monotonicity of β . To see this we consider two piecewise constant functions θh , θ¯h ∈ Wh satisfying (21) for any wh ∈ Wh . Subtracting the resulting two equalities, we obtain

(βε−1 (θh ) − βε−1 (θ¯h ), wh )h + τ kwh k21,h = τ (r (βε−1 (θh )) − r (βε−1 (θ¯h )), wh )h .

(25)

With wh = θh − θ¯h , using (A3), the Cauchy inequality and the inequality of means, we obtain

τ (r (βε−1 (θh )) − r (βε−1 (θ¯h )), θh − θ¯h )h ≤

1 2

(βε−1 (θh ) − βε−1 (θ¯h ), θh − θ¯h )h +

τ2 kθh − θ¯h k2h . 2C¯

By the discrete Poincaré inequality, uniqueness follows whenever τ ≤ C , where C > 0 is a constant that does not depend on the parameters τ , h, or ε . For the existence, we use Lemma 1.4, p. 140 in [24], which is an abstract result for finite dimensional Hilbert spaces. In this sense we define the continuous mapping P : Wh → Wh

Pθ = Φ =

X

αi χTi ,

Ti ∈T

where χTi is the indicator function of the triangle Ti , while αi ∈ R are given by

αi = |Ti |(βε−1 (θi ) − βε−1 (θhk,−i 1 )) + τ

X `ij ∈Ei

σij (θi − θj ) − τ |Ti |r (βε−1 (θi )).

For any θ ∈ Wh , we use (21) to estimate the inner product (P θ , θ )h :

(P θ , θ )h = (βε−1 (θ ), θ )h + τ kθk21,h − (βε−1 (θhk−1 ), θ )h − τ (r (βε−1 (θ )), θ )h . The first term on the right is bounded by

(βε−1 (θ ), θ )h ≥

1 2

(βε−1 (θ ), θ )h +

C 2

kθ k2h ,

whereas for the third term we use (4) and the Cauchy inequality to obtain

|(βε−1 (θhk−1 ), θ )h | ≤ C kθhk−1 kh kθ kh ≤

C0

δ

kθhk−1 k2h + δkθ k2h .

Proceeding as for the a priori estimates (22), the last term yields

τ |(r (βε−1 (θ )), θ )h | ≤

1 4

(βε−1 (θ ), θ )h + C τ 2 kθ k2h .

(26)

2140

I.S. Pop et al. / Journal of Computational and Applied Mathematics 234 (2010) 2135–2142

Choosing δ properly, for moderately small τ the above inequalities as well as (17) imply

(P θ , θ )h ≥

1 4

τ

(βε−1 (θ ), θ )h +

2

kθ k21,h − K ,

with K = C˜ kθhk−1 k2h for some fixed C˜ . Now existence follows by the result mentioned above.



For obtaining the error estimates we proceed as in the time discrete case and use the discrete Green operator Gh : L2 (Ω ) → Wh defined by

(Gh ψ, ϕ)1,h = (ψ, ϕ)h ,

(27)

for all ϕ ∈ Wh . As in the spatially continuous case, for any ψ ∈ Wh one gets

kGh ψk1,h = kψk−1,h ,

kψk−1,h ≤ C kψkh .

(28)

Gh is well defined as the FV approximation of the Poisson equation with homogeneous Dirichlet boundary conditions and an L2 right hand side (see [22] and [23]). As shown in [25,22,28] and [23], for any ψ ∈ L2 (Ω ) one has the estimates

k(G − Gh )Ψ k1,h ≤ ChkΨ k.

(29)

To estimate the error due to the spatial discretization, we define for each time step t = kτ k,h

euk,h := βε−1 (θ k ) − βε−1 (θhk ), eθ

:= θ k − θhk ,

(30)

see also (8). The errors defined above are estimated in the following lemma: Lemma 3.2. Assume (A1)–(A3), let θ k and θhk solving (5), respectively (21). We have sup keuk,h k2−1 + C τ

k=1,n

n X

kekθ,h k2 + C τ

k=1

n X


(eku,h , ekθ,h ) ≤ C ke0u,h k2−1 + h2 /ε ,

k=1

provided τ is reasonably small. Proof. We take ϕ = Geku,h ∈ H01 in (5) and ϕ = Gh eku,h ∈ Wh in (21), subtract the resulting and use (16) to obtain

  (euk,h − euk−1,h , Geku,h ) + τ (∇θ k , ∇ Geku,h ) − (θhk , Gh eku,h )1,h = −(β −1 (θhk ) − β −1 (θhk−1 ), (G − Gh )eku,h )h + τ (r (β −1 (θhk )), (G − Gh )eku,h )h + τ (r (β −1 (θ k )) − r (β −1 (θhk )), Geku,h ).

(31)

We sum up the above for k = 1, . . . , p, denote the resulting terms by S1 , . . . , S5 , and proceed by estimating them separately. By (9), for S1 we have 2S1 = 2

p X (∇ G(eku,h − eku−1,h ),

∇ Geku,h ) = kepu,h k2−1 − ke0u,h k2−1 +

k=1

p X

keku,h − euk−1,h k2−1 .

(32)

k=1

Using (4), (9), (16) and (27), S2 becomes S2 = τ

p p X X (θ k , eku,h ) − (θhk , eku,h )h = τ (ekθ,h , eku,h ) k=1

k=1

p

≥

p X k,h τ X k,h k,h τε X (eθ , eu ) + C τ keku,h k2 . ke θ k2 + p

3 k=1

k=1

3

(33)

k=1

To estimate S3 we use the estimates (24) and (29), as well as (19) and (28) to obtain

|S3 | ≤

p X

kβ −1 (θhk ) − β −1 (θhk−1 )k−1,h k(G − Gh )eku,h k1,h

k=1

≤ δ1

p X

kβ −1 (θhk ) − β −1 (θhk−1 )k2−1,h +

k=1

≤C

h2

ε

+

p τε X

12 k=1

keku,h k2 ,

p Ch2 X

4δ1 k=1

keku,h k2 (34)

where in the above we have taken δ1 = O(h2 /(τ ε)). Alternatively, the L2 estimates for β −1 (θh ) and β −1 (θhk ) imply P τ pk=1 keuk,h k2 ≤ C . For δ1 = h/τ this yields |S3 | ≤ Ch.

I.S. Pop et al. / Journal of Computational and Applied Mathematics 234 (2010) 2135–2142

2141

For S4 we use (A3) and (17), and proceed in a similar manner to get

|S4 | ≤ τ

p X

kr (β −1 (θhk ))k−1,h k(G − Gh )eku,h k1,h

k =1

≤ Cτ h

p X

kr (β −1 (θhk ))kh keku,h k ≤ C

k=1

h2

ε

+

p τε X

12 k=1

keku,h k2 .

(35)

As above, the alternative estimate is |S4 | ≤ Ch. By (A3), the last term gives

|S5 | ≤ τ

p X

kr (β −1 (θhk )) − r (β −1 (θhk−1 ))kkGeku,h k

k =1

≤ C˜ τ δ2

p X

(ekθ,h , eku,h ) +

k=1

p τ X kek,h k2 . δ2 k=1 u −1

Taking δ2 = 1/(6C˜ ), and using (31)–(36), the discrete Gronwall lemma provides the result.

(36) 

Remark 3.2. As following from the proof, the ratio h2 /ε in the estimates can be replaced by h. Furthermore, the initial error can be made arbitrarily small. To see this, notice that e0u,h = βε−1 (θ 0 ) − Ph (βε−1 (θ 0 )). As mentioned in the beginning of Section 2, a mollifying step is involved in constructing a θ 0 that is H 1 , having an ε -uniformly bounded norm. By (13) this gives ke0u,h k ≤ Ch. The FV scheme is convergent, as follows from Theorem 2.1 and of Lemma 3.1: Theorem 3.2. Assuming (A1)–(A3), the FV approximation converges to the solution of Problem WP as τ , h and ε & 0. The following estimates hold k

sup k¯u − k=1,n

βε−1 (θhk )k2−1

T

Z + 0

(βε (u(t )) − θ∆,h (t ), u(t ) − βε−1 (θ∆,h (t )))dt

  h2 + kβ(u) − θ∆,h k2L2 (Q ) ≤ C τ + ε + , T ε where (similarly to the time discrete case) θ∆,h (t ) = θhk whenever t ∈ (tk−1 , tk ]. Remark 3.3. The above estimates are sub-optimal when compared to the ones for the heat equation. As mentioned before, in a certain framework, one can obtain optimal (first order) estimates for the time discretization. To extend such a result to the FV discretization, one needs higher order estimates in (29), as suggested, for example, in [22,26] and [27].

Acknowledgments MS and ISP acknowledge the support of Fondecyt through the project #7070129, by which this joint work has been initiated. MS has been supported by Fondecyt project # 1070694, FONDAP and BASAL projects CMM, Universidad de Chile, and CI2 MA, Universidad de Concepción. OVV acknowledges the support of Dirección de Investigación de la Universidad del Bío-Bío, DIPRODE project #061008 1/R. The work of ISP was supported by the Dutch government through the national program BSIK: knowledge and research capacity, in the ICT project BRICKS (http://www.bsik-bricks.nl), theme MSV1, as well as the International Research Training Group NUPUS. References [1] E. Cariaga, F. Concha, M. Sepúlveda, Flow through porous media with applications to heap leaching of copper ores, Chem. Eng. J. 111 (2005) 151–165. [2] E. Cariaga, F. Concha, M. Sepúlveda, Convergence of a MFE-FV method for two phase flow with applications to heap leaching of copper ores, Comput. Methods Appl. Mech. Engrg. 196 (2007) 2541–2554. [3] M. Ohlberger, Convergence of a mixed finite element - finite volume method for the two phase flow in porous media, East-West J. Numer. Math. 5 (1997) 183–210. [4] R.H. Nochetto, C. Verdi, Approximation of degenerate parabolic problems using numerical integration, SIAM J. Numer. Anal. 25 (1988) 784–814. [5] T. Arbogast, M.F. Wheeler, N.Y. Zhang, A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media, SIAM J. Numer. Anal. 33 (1996) 1669–1687. [6] F.A. Radu, I.S. Pop, P. Knabner, Order of convergence estimates for an Euler implicit, mixed finite element discretization of Richards’ equation, SIAM J. Numer. Anal. 42 (2004) 1452–1478. [7] F.A. Radu, I.S. Pop, P. Knabner, Error estimates for a mixed finite element discretization of some degenerate parabolic equations, Numer. Math. 109 (2008) 285–311.

2142

I.S. Pop et al. / Journal of Computational and Applied Mathematics 234 (2010) 2135–2142

[8] E. Schneid, P. Knabner, F.A. Radu, A priori error estimates for a mixed finite element discretization of the Richards’ equation, Numer. Math. 98 (2004) 353–370. [9] B. Riviere, M.F. Wheeler, Discontinuous Galerkin methods for flow and transport problems in porous media, Comm. Numer. Methods Engrg. 18 (2002) 63–68. [10] R.A. Klausen, F.A. Radu, G.T. Eigestad, Convergence of MPFA on triangulations and for Richards’ equation, Internat. J. Numer. Methods Fluids 58 (2008) 1327–1351. doi:10.1002/fld.1787. [11] M. Afif, B. Amaziane, Convergence of finite volume schemes for a degenerate convection-diffusion equation arising in flow in porous media, Comput. Methods Appl. Mech. Engrg. 191 (2002) 5265–5286. [12] R. Eymard, T. Gallouët, R. Herbin, A. Michel, Convergence of a finite volume scheme for nonlinear degenerate parabolic equations, Numer. Math. 92 (2002) 41–82. [13] R. Eymard, D. Hilhorst, M. Vohralík, A combined finite volume-nonconforming/mixed-hybrid finite element scheme for degenerate parabolic problems, Numer. Math. 105 (2006) 73–131. [14] H.W. Alt, S. Luckhaus, Quasilinear elliptic-parabolic differential equations, Math. Z. 183 (1983) 311–341. [15] F. Otto, L1 -contraction and uniqueness for quasilinear elliptic-parabolic equations, J. Differential Equations 131 (1996) 20–38. [16] W. Jäger, J. Kačur, Solution of doubly nonlinear and degenerate parabolic problems by relaxation schemes, M2 AN (Math. Model. Numer. Anal.) 29 (1995) 605–627. [17] I.S. Pop, Regularization Methods in the Numerical Analysis of Some Degenerate Parabolic Equations, Ph.D. Thesis, Universitatea ‘‘Babeş-Bolyai’’ ClujNapoca, Romania, (1998). [18] I.S. Pop, Error estimates for a time discretization method for the Richards’ equation, Comput. Geosci. 6 (2002) 141–160. [19] I.S. Pop, W.A. Yong, A numerical approach to degenerate parabolic equations, Numer. Math. 92 (2002) 357–381. [20] I.S. Pop, Numerical schemes for degenerate parabolic problems, in: Progress in Industrial Mathematics at ECMI 2004, in: A.Di Bucchianico, R.M.M. Mattheij, M.A. Peletier (Eds.), Mathematics in Industry, Vol. 8, Springer-Verlag, Heidelberg, 2006, pp. 513–517. [21] J. Rulla, Error analysis for implicit approximations to solutions to Cauchy problems, SIAM J. Numer. Anal. 33 (1996) 68–87. [22] R. Eymard, T. Gallouët, R. Herbin, Finite volume methods, in: Handbook of Numerical Analysis, Vol. VII, North-Holland, Amsterdam, 2000, pp. 713–1020. [23] T. Gallouët, R. Herbin, M.H. Vignal, Error estimates on the approximate finite volume solution of convection diffusion equations with general boundary conditions, SIAM J. Numer. Anal. 37 (2000) 1935–1972. [24] R. Temam, Navier–Stokes Equations: Theory and Numerical Analysis, AMS Chelsea Publishing, Providence, RI, 2001. [25] Y. Coudiére, T. Gallouët, R. Herbin, Discrete Sobolev inequalities and Lp error estimates for approximate finite volume solutions of convection diffusion equations, M2 AN (Math. Model. Numer. Anal.) 35 (2001) 767–778. [26] J. Baranger, J.F. Maître, F. Oudin, Connection between finite volume and mixed finite element methods, RAIRO - Modél. Math. Anal. Numér. 30 (1996) 445–465. [27] Z. Chen, Expanded mixed finite element methods for quasilinear second order elliptic problems, M2 AN (Math. Model. Numer. Anal.) 32 (1998) 501–520. [28] R. Herbin, An error estimate for a finite volume scheme for a diffusion convection problem on a triangular mesh, Numer. Methods Partial Differential Equations 11 (1995) 165–173. [29] D. Kröner, Numerical Schemes for Conservation Laws, Wiley-Teubner, Stuttgart, 1997.