Error estimates for the finite volume discretization for the
porous medium equation
Citation for published version (APA):
Sepúlveda, M., Vera Villagrán, O. P., Pop, I. S., & Radu, F. A. (2008). Error estimates for the finite volume discretization for the porous medium equation. (CASA-report; Vol. 0813). Technische Universiteit Eindhoven.
Document status and date: Published: 01/01/2008
Document Version:
Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)
Please check the document version of this publication:
• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.
• The final author version and the galley proof are versions of the publication after peer review.
• The final published version features the final layout of the paper including the volume, issue and page numbers.
Link to publication
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal.
If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:
www.tue.nl/taverne
Take down policy
If you believe that this document breaches copyright please contact us at:
openaccess@tue.nl
Error estimates for the finite volume discretization for the porous
medium equation
M. Sep´ulveda1, O. P. Vera Villagr´an2, I. S. Pop3 and F. A. Radu4,5
1CI2MA and Departamento de Ingenier´ıa Matem´atica, Universidad de Concepci´on, Casilla 160-C, Concepci´on, Chile 2Universidad del B´ıo-B´ıo, Collao 1202, Casilla 5-C, Concepci´on, Chile
3CASA, Eindhoven University of Technology, P. O. Box 513, 5600 MB Eindhoven, The Netherlands 4UFZ-Helmholtz Center for Environmental Research, Permoserstr. 15, D-04318 Leipzig, Germany
5University of Jena, W¨ollnitzerstr. 7, D-07749, Jena, Germany
Abstract
In this paper we analyze the convergence of a numerical scheme for a class of degenerate parabolic problems. Such problems are often used to model reactions in porous media, and in-volve a nonlinear, possibly vanishing diffusion. The scheme considered here inin-volves the Kirchhoff transformation coupled with the regularization of the nonlinearity, and is based on the Euler implicit time stepping and the triangle based finite volume spatial discretization. We prove the convergence of the approach by giving estimates for the error in terms of the discretization and regularization parameter.
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 (see, e.g. [13]). When 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. In this way 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. Conse-quently, 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 [5]. The convergence of such schemes has been investigated in [6] and [20], 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 [19] for the finite element discretization, [3], [24], [25], [28] for MFE schemes, [26] for a DG approach, and [16] for a multipoint flux approximation method. FV schemes for porous media models are analyzed rigorously in [1], [10], whereas a MFEM-FV approach is considered in [11]. In these cases the convergence is obtained by means of compactness arguments, and no estimates for the approximation error are given.
Here we consider the FV discretization of the degenerate parabolic equation
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) ∈ H1
0(Ω).
(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) = umfor 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. This involves nonzero convection, and possibly reaction.
For the ease of presentation we have considered homogeneous Dirichlet boundary conditions; these can be extended straightforwardly to other types. In what follows we use standard notations for the spaces of functions, norms and scalar products: L2(Ω), H01(Ω), or its dual H−1(Ω). With X being one of the spaces before, 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·kkdenotes the norm in Hk(Ω). Moreover, we often write u or u(t) instead of u(t, x) and use C to denote a generic positive constant independent of the discretization parameters or the function itself. Having this in mind we seek for the weak solution of Problem P, solving:
Problem WP.
Find u ∈ H1(0, T ; H−1(Ω)) such that β(u) ∈ L2(0, T ; H1
0(Ω)), and for all ϕ ∈ H01(Ω) we have
(∂tu, ϕ) + (∇β(u), ∇ϕ) = (r(u), ϕ), (2)
for all t ∈ (0, T ], whereas u(0) = u0in H−1(Ω).
The existence and uniqueness of a solution for Problem WP is proven, for example, in [2] and [21]. Notice that β(u) has a better regularity than u, this being a property that we will exploit in what follows. Furthermore, we employ a regularization step in constructing the numerical scheme. Specifically, with ε > 0 a given parameter, the nonlinear function β is approximated by βε satisfying β0 ≥ ε. For simplicity we consider the global perturbation of β
βε(u) ≡ β(u) + εu. (3)
Other possibilities can be considered as well, and the analysis is absolutely similar. Then βεis invertible and we have
ε ≤ βε0(u) ≤ C, 0 ≤ β0ε(u) − β0(u) ≤ ε, and C ≤ β−1ε 0
(βε(u)) ≤ 1
ε. (4)
2
The time discretization
To discretize Problem WP in time we consider the Euler implicit scheme, combined with regularization. This idea has been used for constructing effective numerical algorithms. In this sense we mention [15], where compactness arguments are considered for showing the convergence in a general setting. We also mention [29] for a linear regularization based scheme, as well as [24], where (2) is considered in a mixed form. Linearized approaches are also discussed in [22], in connection with finite elements. In all these situations first order schemes are considered, this choice being determined by the lacking regularity of the solution.
As noticed before, β(u) is more regular than u itself. Therefore we first calculate θ ≈ β(u), and then obtain 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 θkof β(u(tk)) as the solution of
Problem WTk. Given θk−1∈ H1
0(Ω), find θk∈ H01(Ω) such that for all ϕ ∈ H01(Ω)
(βε−1(θk) − βε−1(θk−1), ϕ) + τ (∇θk, ∇ϕ) = τ (r(βε−1(θk)), ϕ). (5) At this point we have to specify the initial data. A straightforward choice is θ0 = βε(u0). Whenever u0 ∈ H1
0(Ω), this gives θ0 in the same space. However, in (A2) we have only required u0 ∈ L2(Ω). Following the discussion in [22], Chapter 3, one can consider θ0 = β(u0) + ερµ∗ u0, where ρµ is a mollifier having a compact support in B(0, µ). With µ = O(ε), θ0is bounded uniformly in H1, 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.
Remark 2.1 Inverting βε may be tedious. Furthermore, since function calls increase the computing time significantly, in implementing the scheme we have constructed a look-up table of values forβε. At the expense of some computer memory and involving a simple (linear) interpolation step for obtaining values which are not included in the table, the computing time is significantly reduced, whereas the errors are controlled by an adaptive choice of the interpolation knots. Moreover, the monotonicity of βε allows a fast searching in this table and therefore the values ofβε or the inverse can be obtained efficiently.
Existence and uniqueness for the time discrete problems WPk is provided by standard results for non-linear elliptic equations. Furthermore, assuming that the initial data is essentially bounded, the sequence of solutions θkremain essentially bounded uniformly in k. The following estimates for the time discrete approximation are obtained in [22], Chapter 3 (see also [23]).
Lemma 2.1 Assume (A1) - (A3). With θksolving the problems WTk, for all0 ≤ p ≤ n we have kθpk2+ n X k=1 n (βε−1(θk) − βε−1(θk−1), θk− θk−1) + kθk− θk−1k2+ τ k∇θkk2o≤ C. (6)
Remark 2.2 The estimate (6) immediately implies n
X k=1
kβε−1(θk) − β−1ε (θk−1)k2−1 ≤ Cτ. (7)
To prove the above we use (5) and the Poincar´e inequality to obtain
|(βε−1(θk) − βε−1(θk−1), φ)| ≤ τ k∇θkk + kr(βε−1(θk))k
for allφ ∈ H01(Ω) s.t. k∇φk = 1. The inequality now follows by (A3) and the above estimates.
To give the error estimates introduce some notations. For any function f : QT → R that is integrable in time we define ¯ fk := 1 τ Z kτ (k−1)τ f (s, ·)ds, if k ≥ 1. Further, ¯f0 := f (0, ·). The errors are obtained in terms of ekuand ekθ defined as
eku := ¯uk− βε−1(θk), ekθ := β(u)k− θk, (8) where k ≥ 0. Given a sequence {fk ∈ H1
0(Ω), k = 1, n}, the piecewise constant extension in time f∆ is defined as f∆(t) = θkfor t ∈ (tk−1, tk].
Further, G : H−1(Ω) → H01(Ω) stands for the Green operator defined by
(∇Gψ, ∇ϕ) = (ψ, ϕ), (9)
for all ϕ ∈ H01(Ω), where ψ ∈ H−1(Ω). Therefore
k∇Gψk = kψk−1, kψk−1 ≤ Ckψk, (10)
where the last inequality applies only if ψ ∈ L2(Ω). With this notations we have ([22], Chapter 3): Theorem 2.1 Assuming (A1) -(A3), if u solves Problem WP and θk solves Problem WTk (k = 1, n), then sup k=1,n kek uk2−1+ RT 0 (βε(u(t)) − θ∆(t), u(t) − β −1 ε (θ∆(t)))dt + kβ(u) − θ∆k2L2(Q T) ≤ C {τ + ε} .
The estimates above are valid for a more general setting, where convection terms can also be included. Using the results in [27], these estimates become optimal, Cτ2+ ε2 . This holds in a simplified setting, e.g. in the absence of convection and if β is maximal monotone having the range R.
3
The finite volume discretization
In this section we refer to the framework in [9] (see also [12] and [20]) 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, xi- the center of the circumcircle of Ti,
`ij - the edge between Tiand Tj (where Tj ∈ Ni),
d(xi, `ij) - the distance from xi to `ij, and dij = d(xi, `ij) + d(xi, `ij) if `ij ∈ Eint, σij = |`ij|/dij - the ”transmissibility” through `ij
nij - the outward unit normal to `ij pointing into Tj (Tj ∈ Ni).
Notice that the assumptions on the triangularization 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 can define the finite dimensional subspace of L2(Ω)
Wh := {v ∈ L2(Ω)/ v is constant on any T ∈ T }, (11) which is a space spanned by the triangle indicator functions {χT / T ∈ T }. Furthermore we define the projector
Ph: L2(Ω) → Wh, (Phw − w, wh) = 0 (12) for any wh ∈ Wh. With s ∈ {0, 1}, a constant C > 0 exists such that
for all w ∈ Hs(Ω). Furthermore, for any w ∈ L2(Ω) and Ti ∈ T we have ¯ wi := (Phw)|Ti = 1 |Ti| Z Ti w(x)dx. (14)
In analogy to the spatially continuous case, for any u, v ∈ L2(Ω) we can define the following discrete inner products
(u, v)h := X Ti∈T |Ti|¯uiv¯i, (u, v)1,h := X `ij∈E σij(¯ui− ¯uj)(¯vi− ¯vj), (15)
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, Phv)h = (Phu, v)h. (16) Furthermore, in [9] the following discrete Poincar´e inequality is given:
kuk ≤ Ckuk1,h, (17)
for all u ∈ Wh, where C > 0 does not depend on h or u. Using (16) and (13), one immediately obtains |(u, v) − (u, v)h| = |(u − Phu, v)| = |(u − Phu, v − Phv)| ≤ Chs+pkukskvkp, (18) for all u ∈ Hs(Ω) and v ∈ Hp(Ω), s, p ∈ {0, 1}. Moreover, as in the continuous case we can define the discrete H−1norm
kuk−1,h = sup wh∈Wh,kwhk1,h=1
|(u, wh)h|. (19)
At this point we can follow [9] and [12] to introduce the finite volume scheme for the time discrete problem WTk. Specifically, given θk−1h ∈ Wh, we seek for θk ∈ Whsuch that for all Ti∈ T it holds
|Ti|(βε−1(θkh,i) − βε−1(θh,ik−1)) + τ X `ij∈Ei
σij(θkh,i− θkh,j) = τ |Ti|r(βε−1(θh,ik )), (20)
where θh,i= θh|Ti.
The above scheme can be brought to a weak form. To this aim, for any wh ∈ Whwe 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∈ Whsuch that for all wh ∈ Wh
(βε−1(θhk) − βε−1(θk−1h ), wh)h+ τ (θhk, wh)1,h= τ (r(βε−1(θhk)), wh)h. (21) To complete the scheme, we define the initial data θ0h = βε(Ph(βε−1(θ0))) ∈ Wh.
First we give some stability properties of the FV approximation that are similar to the time discrete case.
Lemma 3.1 Assume (A1) - (A3), and let θhksolving (21). For any1 ≤ p ≤ n we have kβ−1ε (θhp)k2h+Pn k=1kβ −1 ε (θhk) − β −1 ε (θk−1h )k2h+ τ Pn k=1kθhkk21,h ≤ C, Pn k=1(βε−1(θkh) − β −1 ε (θhk−1), θkh− θ k−1 h )h+ Pn k=1kθhk− θk−1k2h ≤ C. (22) 5
Proof. We start by noticing that for any wh∈ Wh, by (A3) and (4) we get Ckwhk21,h ≤ (β −1 ε (wh), wh)1,h, as well as |(r(wh), wh)h| ≤ C(w¯ h, β(wh))h 1/2 kwhkh ≤ Ckwhk2 h. (23)
Next we take wh = βε−1(θhk) 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 2 kβ −1 ε (θ p h)k 2 h− kβε−1(θ0h)k2h+ Pp k=1kβ −1 ε (θkh) − βε−1(θhk−1)k 2 h +τPp k=1(θhk, β −1 ε (θkh))1,h= τ Pp k=1(r(β −1 ε (θhk)), β −1 ε (θhk))h.
Using now the inequalities in (23), the first part of the estimates follow are a direct consequence of the Gronwall lemma.
The second estimates are obtained by testing (21) with wh= θhk− θk−1h , and using the assumptions
on β and r.
Remark 3.1 As in the spatially continuous case, the estimates above immediately imply n
X k=1
kβ−1ε (θhk) − βε−1(θhk−1)k2−1,h≤ Cτ. (24)
The above estimates are obtained assuming that the fully discrete problems are uniquely solvable: Theorem 3.1 Assume (A1) - (A3), and let θhk−1 ∈ Wh be given. Then the fully discrete problem WDk 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, both satisfying (21) for any wh ∈ Wh. Taking then wh = θh− ¯θh and subtracting the resulting two identities we obtain
(βε−1(θh) − βε−1(¯θh), wh)h+ τ kwhk21,h= τ (r(β −1
ε (θh)) − r(β−1ε (¯θh)), wh)h. (25) Using (A3), the Cauchy inequality as well as the inequality of means, we estimate the term on the right by τ (r(βε−1(θh)) − r(βε−1(¯θh)), wh)h ≤ 1 2(β −1 ε (θh) − βε−1(¯θh), wh)h+ τ2 2 ¯Ckwhk 2 h.
By the discrete Poincar´e 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 see Lemma 1.4, p. 140 in [30], which is an abstract result for finite dimensional Hilbert spaces. In this sense we define the continuous mapping P : Wh→ Whas
Pθ = Φ = X Ti∈T
αiχTi,
where χTiis the indicator function of the triangle Ti, while αi ∈ R are given by
αi = |Ti|(βε−1(θi) − βε−1(θk−1h,i )) + τ X `ij∈Ei
σij(θi− θj) − τ |Ti|r(βε−1(θi)). (26)
Given a θ ∈ Wh, we use (21) to estimate the inner product (Pθ, θ)h: (Pθ, θ)h = (β−1ε (θ), θ)h+ τ kθk1,h2 − (βε−1(θk−1h ), θ)h− τ (r(β
−1
The first term on the right is bounded by (βε−1(θ), θ)h≥ 1 2(β −1 ε (θ), θ)h+ C 2kθ| 2 h, whereas for the third term we use (4) and the Cauchy inequality to obtain
|(βε−1(θk−1h ), θ)h| ≤ Ckθk−1h khkθkh≤ C0 δ kθ k−1 h k 2 h+ δkθ|2h. Proceeding as for the apriori estimates (22), the last term yields
τ |(r(βε−1(θ)), θ)h| ≤ 1 4(β
−1
ε (θ), θ)h+ Cτ2kθ|2h.
Choosing δ properly, for moderately small τ the above inequalities as well as (17) imply (Pθ, θ)h ≥ 1 4(β −1 ε (θ), θ)h+ τ 2kθk 2 1,h− K, with K = ˜Ckθk−1h k2
h for some fixed ˜C. The existence of a solution is provided now by the abstract
result mentioned above.
For obtaining the error estimates we proceed as in the time discrete case and use the discrete Green operator Gh : L2(Ω) → Whdefined by
(Ghψ, ϕ)1,h= (ψ, ϕ)h, (27)
for all ϕ ∈ Wh. As in the spatially continuous case, for any ψ ∈ Whone gets
kGhψk1,h = kψk−1,h, kψk−1,h ≤ Ckψkh. (28) Ghis well defined as the FV approximation of the Poisson equation with homogeneous Dirichlet bound-ary conditions and an L2 right hand side (see [9] and [12]). Furthermore, as follows from these papers (see also [8] and [18]), for any ψ ∈ L2(Ω) one has the error estimates
k(G − Gh)Ψk1,h ≤ ChkΨk. (29)
To determine the error added by the spatial discretization to the time discretization we define for each time step t = kτ (k = 0, . . . , n) the errors
ek,hu := βε−1(θk) − βε−1(θkh), ek,hθ := θk− θkh, (30) see also (8). The errors defined above are estimated in the following lemma:
Lemma 3.2 Assume (A1) - (A3), let θkandθk
hsolving (5), respectively (21). We have sup k=1,n kek,hu k2 −1+ Cτ Pn k=1ke k,h θ k2+ Cτ Pn k=1(e k,h u , ek,hθ ) ≤ C kek,hu k2−1+ h2 ε ! , providedτ is reasonably small.
Proof. We take ϕ = Gek,hu ∈ H01 in (5) and ϕ = Ghek,hu ∈ Wh in (21), subtract the resulting and use (16) to obtain yields (ek,hu − ek−1,hu , Gek,hu ) + τ h (∇θk, ∇Gek,hu ) − (θhk, Ghek,hu )1,h i = −(β−1(θhk) − β−1(θk−1h ), (G − Gh)ek,hu )h +τ (r(β−1(θkh)), (G − Gh)ek,hu )h+ τ (r(β−1(θk)) − r(β−1(θkh)), Ge k,h u ). (31) 7
We sum up the above for k = 1, . . . , p, denote the resulting terms by S1, . . . , S5, and proceed by estimating them separately. By (27), S1gives
S1 = Ppk=1(G(e k,h u − ek−1,hu ), Gek,hu )−1 = 12nkep,hu k2−1− ke0,hu k2−1+ Pp k=1ke k,h u − ek−1,hu k2−1 o . (32)
Using (4), (9), (16) and (27), S2becomes
S2 = τPpk=1(θk, ek,hu ) − (θhk, ek,hu )h = τ Pp k=1(e k,h θ , e k,h u ) ≥ τ3 Pp k=1(e k,h θ , e k,h u ) + CτPpk=1kek,hθ k2+τ ε3 Ppk=1kek,hu k2. (33)
To estimate S3we use the estimates (24) and (29), as well as (19) and (28) to obtain |S3| ≤ Ppk=1kβ−1(θhk) − β−1(θk−1h )k−1,hk(G − Gh)ek,hu k1,h ≤ δ1Ppk=1kβ−1(θhk) − β−1(θk−1h )k2−1,h+Ch 2 4δ1 Pp k=1ke k,h u k2 ≤ Chε2 +τ ε12Pp k=1ke k,h u k2, (34)
where in the above we have taken δ1 = O(h2/(τ ε)). Alternatively, one may use the L2 estimates for β−1(θh) and β−1(θkh) to show that τ
Pp k=1ke
k,h
u k2 ≤ C, yielding for δ1 = h/τ : |S3| ≤ Ch. For S4we use (A3) and (17), and proceed in a similar manner to get
|S4| ≤ τPp k=1kr(β −1(θk h))k−1,hk(G − Gh)ek,hu k1,h ≤ Cτ hPp k=1kr(β −1(θk h))khkek,hu k ≤ Ch 2 ε + τ ε 12 Pp k=1ke k,h u k2. (35)
As above, the alternative estimate is |S4| ≤ Ch. By (A3), the last term gives |S5| ≤ τPpk=1kr(β−1(θhk)) − r(β−1(θhk−1))kkGek,hu k
≤ Cτ δ˜ 2Ppk=1(ek,hθ , ek,hu ) +δτ2Ppk=1kek,hu k2−1.
(36)
Taking now δ2 = 1/(6 ˜C), using (32) - (36) into (31), the result is a direct consequence of the discrete
Gronwall lemma.
Remark 3.2 As following from the proof, the ratio h2/ε in the above estimates can be replaced by h. Furthermore, the initial error can be made arbitrarily small. To see this we first notice that e0,hu = βε−1(θ0) − Ph(βε−1(θ0)). As mentioned in the beginning of Section 2, a mollifying step is involved in constructing aθ0 that isH1, having anε-uniformly bounded norm. By (13) this gives ke0,hu k ≤ Ch.
The convergence of the scheme is a straightforward consequence of 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
sup k=1,n k¯uk− βε−1(θkh)k2−1+R0T(βε(u(t)) − θ∆,h(t), u(t) − β−1ε (θ∆,h(t)))dt +kβ(u) − θ∆,hk2L2(Q T) ≤ C n τ + ε + hε2 o , where (similarly to the time discrete case)θ∆,h(t) = θhkwhenevert ∈ (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 e.g. in [9], [4] and [7].
Acknowledgements. 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 CI2MA, Universidad de Con-cepci´on and CMM, Universidad de Chile, through the Fondecyt project #1070694 and Fondap in Ap-plied Mathematics (Project #15000001). OVV acknowledges the support of Direcci´on de Investigaci´on 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] M. Afif, B. Amaziane Convergence of finite volume schemes for a degenerate convectiondiffu-sion equation arising in flow in porous media, Comput. Methods Appl. Mech. Eng. 191 (2002) 52655286.
[2] H.W. Alt, S. Luckhaus, Quasilinear elliptic-parabolic differential equations, Math. Z. 183 (1983) 311-341.
[3] T. Arbogast, M. F. Wheeler and N. Y. Zhang, A nonlinear mixed finite element method for a de-generate parabolic equation arising in flow in porous media, SIAM J. Numer. Anal. 33 (1996) 1669-1687.
[4] J. Baranger, J.F. Maˆıtre and F. Oudin, Connection between finite volume and mixed finite element methods, RAIRO - Mod´el. Math. Anal. Num´er. 30 (1996) 445-465.
[5] E. Cariaga, F. Concha, M. Sep´ulveda, Flow through porous media with applications to heap leach-ing of copper ores, Chemical Engineerleach-ing Journal 111 (2005) 151-165.
[6] E. Cariaga, F. Concha, M. Sep´ulveda, 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.
[7] Z. Chen, Expanded mixed finite element methods for quasilinear second order elliptic problems, M2AN 32 (1998) 501-520.
[8] Y. Coudi´ere, T. Gallou¨et, R. Herbin, Discrete Sobolev inequalities and Lp error estimates for ap-proximate finite volume solutions of convection diffusion equations, M2AN (Math. Model. Numer. Anal.) 35 (2001) 767-778.
[9] R. Eymard, T. Gallou¨et, R. Herbin, Finite volume methods, in Handbook of numerical analysis, Vol. VII North-Holland, Amsterdam, 2000, 713-1020.
[10] R. Eymard, T. Gallou¨et, R. Herbin, A. Michel, Convergence of a finite volume scheme for nonlinear degenerate parabolic equations, Numer. Math. 92 (2002) 4182.
[11] 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.
[12] T. Gallou¨et, 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.
[13] R. Helmig, Multiphase Flow and Transport Processes in the Subsurface, Springer-Verlag, Heidel-berg, 1997.
[14] R. Herbin, An error estimate for a finite volume scheme for a diffusion convection problem on a triangular mesh, Numer. Meth. Part. Diff. Equ. 11 (1995) 165-173.
[15] W. J¨ager J. Kaˇcur, Solution of doubly nonlinear and degenerate parabolic problems by relaxation schemes, M2AN (Math. Model. Numer. Anal.) 29 (1995) 605-627.
[16] R. A. Klausen, F. A. Radu, G. T. Eigestad, Convergence of MPFA on triangulations and for Richards’ equation, Internat. J. Numer. Methods Fluids (2008), DOI:10.1002/fld.1787.
[17] D. Kr¨oner, Numerical Schemes for Conservation Laws, Wiley-Teubner, Stuttgart, 1997.
[18] I.D. Mishev, Q.-Y. Chen, A mixed finite volume method for elliptic problems, Numer. Meth. Part. Diff. Equ. 23 (2007) 1122-1138.
[19] R.H. Nochetto, C. Verdi, Approximation of degenerate parabolic problems using numerical inte-gration, SIAM J. Numer. Anal. 25 (1988) 784-814.
[20] 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.
[21] F. Otto, L1-contraction and uniqueness for quasilinear elliptic-parabolic equations, J. Differ. Equations 131 (1996) 20-38.
[22] I.S. Pop, Regularization Methods in the Numerical Analysis of Some Degenerate Parabolic Equa-tions, PhD Thesis, Universitatea ”Babes¸-Bolyai” Cluj-Napoca, Romania, (1998).
[23] I. S. Pop, W. A. Yong, A numerical approach to degenerate parabolic equations, Numer. Math. 92 (2002) 357-381.
[24] 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. [25] 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 .
[26] B. Riviere, M. F. Wheeler, Discontinuous Galerkin methods for flow and transport problems in porous media, Comm. Numer. Methods Engrg. 18 (2002) 63-68.
[27] J. Rulla, Error analysis for implicit approximations to solutions to Cauchy problems, SIAM J. Numer. Anal. 33 (1996) 68-87
[28] 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.
[29] M. Slodiˇcka, Solution of nonlinear parabolic problems by linearization, Preprint M3-92, Comenius University Bratislava, Slovakia (1992).
[30] R. TEMAM, Navier-Stokes equations: theory and numerical analysis, AMS Chelsea Publishing, Providence, RI, 2001.
[31] C. Woodward, C. Dawson, Analysis of expanded mixed finite element methods for a nonlinear parabolic equation modeling flow into variably saturated porous media, SIAM J. Numer. Anal. 37 (2000) 701-724.
