# Quasi-Linear Convection-Dominated Transport Problem Based on Characteristics-Mixed Finite Element Method

Authors: Peixin Zhao; Hongying Man
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

## Abstract

We propose the
characteristics-mixed method for approximating the solution to a
convection-dominated transport problem. The new method is a
combination of characteristic approximation to handle the convection
part in time and a mixed finite element spatial approximation to
deal with the diffusion part. Boundary conditions are incorporated
in a natural fashion. The scheme is locally conservative since fluid
is transported along the approximate characteristics on the discrete
level and the test function can be piecewise constant. Analysis shows
that the scheme has much smaller time-truncation errors than those
of standard methods.


---

## Body

## 1. Introduction

Given an open bounded domainΩ⊂R2 with a smooth boundary Γ and a time interval (0,T], we consider in this paper the following convection-diffusion equation:
(1)(a)∂c∂t+u·∇c-∇·(a(c)∇c)=f(c)inΩ×(0,T],(b)∂c∂n=0onΓ×(0,T],(c)c(x,0)=c0(x)inΩ,
where(1)
c(x,t) denotes, for example, the concentration of a possible substance;(2)
u(x) represents the velocity of the flow;(3)
∇ and ∇· denote the gradient and the divergence operators, respectively;(4)
a(c) is sufficiently smooth and there exist constants a1 and a2 such that 0<a1≤a(c)≤a2<+∞;(5)
f denotes a source term.This equation governs such phenomena as the flow of heat within a moving fluid, the transport of dissolved nutrients or contaminants within the groundwater, and the transport of a surfactant or tracer within an incompressible oil in a petroleum reservoir.For convenience, we assume (1) is Ω-periodic, this is physically reasonable since no-flow boundaries are generally treated by reflection, and because in general interior flow patterns are much more important than boundary effects. Thus, the no-flow boundary condition above can be dropped. Because of molecular diffusion, a(c) is uniformly positive. Although this implies that the equation is uniformly parabolic, in many applications the Peclet number is quite high. Thus, convection dominates diffusion and the equation is nearly hyperbolic in nature. The concentration often develops sharp fronts that are nearly shocking. It is well known that strictly parabolic discretization schemes applied to the problem do not work well when it is convection dominated. It is especially difficult to approximate well the sharp fronts and conserve the material or mass in the system.Effective discretization schemes should recognize to some extent the hyperbolic nature of the equation. Many such schemes have been developed, such as the explicit method of characteristics, upstream-weighted finite difference schemes [1], higher-order Godunov schemes [2, 3], the streamline diffusion method [4], the least-squares-mixed finite element method [5, 6], the modified method of characteristics-Galerkin finite element procedure (MMOC-Galerkin) [7–9], and the Eulerian-Lagrangian localized adjoint method (ELLAM) [10]. Each method has its advantages and disadvantages. Explicit characteristic and Godunov schemes require that a Courant-Friedrich-Lewy (CFL) time step constraint is to be imposed. Upstream weighting tends to introduce into the solution an excessive amount of numerical diffusion near the sharp fronts. Compared with upstream weighting, the streamline diffusion method and the least-squares-mixed finite element method reduce the amount of numerical diffusion. It adds a user-defined amount of bias in the direction of streamline. In ELLAM, it can be difficult to evaluate the resulting integrals.We concentrate on MMOC-Galerkin. It is an implicit scheme, so reasonably large time steps may be used, and it does not numerically diffuse the fronts to a particularly excessive degree. Unfortunately, it has certain inherent difficulties, especially with regard to local mass balance.Mixed finite element method has been found to be very useful for solving flow equations not only for its optimal approximation of the unknown function and the vector flux, but also for its mass conservation due to the use of piecewise test function. Its theories and applications have been extensively discussed, for example, see [11–14].Therefore, to purse the high performance of an algorithm, one should take a nice combination of MMOC and the mixed finite element method. In this way, characteristics-mixed finite element method was considered. However, most integrals in this method involve the transformation images of the test function, so it is difficult to compute in practice.In this paper, we propose a new characteristics-mixed finite element scheme. It is similar to MMOC-Galerkin in that we approximate the hyperbolic part of the equation along the characteristics. We use, however, a mixed finite element spatial discretization of the equation. Since piecewise constant functions can be used as test functions, the scheme conserves mass locally and globally in the discrete level. One of our goals in the paper is to yieldK∥∂2c/∂τ2∥Δt in time truncation. Since the solution changes much less rapidly in the characteristic direction, the bound in time truncation will permit the use of larger time steps, with corresponding improvements in efficiency, at no cost in accuracy. For we use the mixed finite element, we can simultaneously approximate the scalar unknown and the vector flux. In error analysis, we introduce some new techniques to overcome difficulties in dealing with the combination of MMOC and the mixed finite element. The optimal error estimates for the unknown function and the vector flux are given.An outline of the paper is as follows. In Section2, we define an approximation to the characteristics and the characteristics-mixed finite element formulation of the problem. We give the proof of the existence and uniqueness of the discrete problem in Section 3. Several important lemmas are given in Section 4, the critical argument in this section is the proof that for any η∈L2(Ω), |(η-η¯,eh)| is bounded by a constant multiple of (h+Δt)∥η∥(∥εh∥+∥eh∥+∥ρ∥). The optimal order estimates for c-ch,p-ph in L2(Ω) are presented in Section 5.


## 2. The Characteristics-Mixed Finite Element Formulation

We begin this section by introducing some notations. We denote byWk,p(S) the standard Sobolev space of k-differential functions in Lp(S). Let ∥·∥k,p,S be its norm and let ∥·∥k,S be the norm of Hk(S)=Wk,2(S) or Hk(S)2, where we omit S if S=Ω. When k=0, we let L2 denote the corresponding space defined on Ω and ∥·∥ let denote its norm.We also use the following spaces that incorporate time dependence. If[a,b]⊂[0,T],X is a Sobolev space and f(x,t) suitably smooths on Ω×[a,b], we let
(2)Lp(a,b;X)={f:∫ab‍∥f(·,t)∥Xpdt<∞},∥f∥Lp(a,b;X)=(∫ab‍∥f(·,t)∥Xpdt)1/p,
where if p=∞, the integral is replaced by the essential supremum.We list the assumptions about the coefficients   as follows:(3)u∈W1,∞(Ω),(4)|∂a(c)∂c|+|∂2a(c)∂c2|+|∂f(c)∂c|≤K,
here and throughout this paper, K denotes different constants in different places.We also assume that the solutionc of (1) satisfies the following:
(5)(a)c∈L∞(0,T;Hk+1(Ω)),(b)ct∈L2(0,T;Hk+1(Ω))⋂‍L∞(0,T;L∞),(c)ctt∈L2(0,T;L2(Ω)).
Under the above assumptions, we begin to discretize the problem. Let ψ(x)=(1+|u|2)1/2 and let the characteristic direction associated with the operator ct+u·∇c be denoted by τ=τ(x), where
(6)∂∂τ(x)=1ψ(x)∂∂t+uψ(x)·∇,
then, (1)(a) can be put in the form the following:
(7)ψ(x)∂c∂τ-∇·(a(c)∇c)=f(c),(x,t)∈Ω×(0,T].Define the spaces:V=L2(Ω), H={χ∈H(div;Ω);χ·n→=0onΓ} and let p=a(c)∇c, then the mixed variational problem corresponding to (1) is to find (c,p):[0,T]→V×H, such that
(8)(a)(ψ(x)∂c∂τ,v)-(divp,v)=(f,v)∀v∈V,(b)(α(c)p,χ)+(c,divχ)=0∀χ∈H,(c)c(x,0)=c0(x)∀x∈Ω,
where α(c)=1/a(c)≥a0>0 and it is this form that will be discretized below.Remark 1.
The restrictionχ·n→=0 on the space H can be removed if the Dirichlet boundary condition is imposed, that is, the results presented in the following section also hold in the case of the first or the second type boundary value problems.
LetTh be a quasi-regular polygonaization of Ω and let Vh×Hh⊂V×H be the associated Raviart-Thomas-Nedelec space [15] of index k≥0. In the procedure to be used, we will consider a time step Δt>0 and approximate the solution at times tn=nΔt. The characteristic derivative will be approximated basically in the following manner.
Letx-=x-u(x)Δt and note that
(9)ψ(x)∂c∂τ|tn≈ψ(x)c(x,tn)-c(x-,tn-1)(x-x-)2+(Δt)2=c(x,tn)-c(x-,tn-1)Δt,
so, our characteristics-mixed finite element method is the determination of the map (ch,ph):{t0,t1,…,tN}→Vh×Hh satisfying the relations following:
(10)(a)(chn-c-hn-1Δt,vh)-(divphn,vh)=(f(chn),vh)∀vh∈Vh,1≤n≤N,(b)(α(chn)phn,χh)+(chn,divχh)=0∀χh∈Hh,0≤n≤N,(c)ch0=c~0,
where chn=ch(tn), and c-hn-1=chn-1(x-)=chn-1(x-u(x)Δt), c~0 are defined in (18) for t=0.Remark 2.
This scheme is an implicit scheme and it is difficult to compute in practice. For simplicity, we can usef(chn-1),α(chn-1) instead of f(chn),α(chn). Then and (10) becomes an explicit scheme, but the following results still hold.


## 3. The Existence and Uniqueness of the Solution of the Discrete Problem

In this section, we give the proof of the existence and uniqueness of the solution of the discrete problem (10).Theorem 3.
There exists a unique solution(chn,phn)∈Vh×Hh to the characteristics-mixed finite element procedure (10).Proof.
Let{ϕi}i=1M⊂Vh and {ψi}i=1N⊂Hh be two sets of bases, respectively, and let
(11)chn=∑i=1M‍Cinϕi,phn=∑j=1N‍Pjnψj,
then (10) may be written in the following matrix form.
Find(12)Cn=(C1n,C2n,…,CMn)T,Pn=(P1n,P2n,…,PNn)T,
such that
(13)(a)ACn-ΔtBPn=ΔtF(Cn)+Gn-1,(b)D(Cn)Pn+BTCn=0,
where
(14)A=((ϕi,ϕj))M×M,D(Cn)=((α(chn)ψi,ψj))N×N,B=((ϕi,divψj))M×N,F(Cn)=((f(chn),ϕ1)⋮(f(chn),ϕM)),Gn-1=((c-hn-1,ϕ1)⋮(c-hn-1,ϕM)).
Since {ϕi}i=1M and {ψi}i=1N are bases, respectively, and α(chn)≥a0, A and D are two symmetric positive definite matrices. Solve Pn from (13)(b), then substitute it into (13)(a), we have
(15)ACn+ΔtB(D(Cn))-1BTCn=ΔtF(Cn)+Gn-1,
Then,
(16)ACn=Gn-1+F~(Cn)Δt,
where
(17)F~(Cn)=F(Cn)-B(D(Cn))-1BTCn.
Since a(c) and f(c) are globally Lipschitz continuous about c, F~ is globally Lipschitz continuous about Cn too. With C0 given by (10)(c) it is clear that this nonlinear system of algebraic equations may be solved for small Δt, so that (10) defines a unique discrete solution for nΔt≤T.In order to analysis the convergence of the method, it is convenient to introduce the mixed elliptic projection associated with our equations.Let(c~,p~):[0,T]→Vh×Hh be given by the following relations
(18)(a)-(div(p-p~),vh)=0,∀vh∈Vh,0≤t≤T,(b)(α(c)(p-p~),χh)+(c-c~,divχh)=0,∀χh∈Hh,0≤t≤T.Letehn=chn-c~n,ρn=c~n-cn,εhn=phn-p~n, and σn=p~n-pn. By [11] we know that for p=2 or +∞ and 1≤s≤k+1, consider
(19)∥ρ∥Lp(0,T;L2)≤K∥c∥Lp(0,T;Hs+δ0k)hs,(20)∥σ∥Lp(0,T;L2)≤K∥c∥Lp(0,T;Hs+1)hs.
It follows easily that
(21)∥ρt∥Lp(0,T;L2)≤K(∥c∥Lp(0,T;Hs+δ0k)+∥ct∥Lp(0,T;Hs+δ0k))hs.
We also have
(22)∥p~∥L∞(0,T;L∞)≤K,∥p~t∥L∞(0,T;L∞)≤K.
By (8), (10), and (18), we obtain the error equation in the following form
(23)(a)(ehn-e-hn-1Δt,vh)-(divεhn,vh)=(ψ∂cn∂τ-cn-c-n-1Δt,vh)-(ρn-ρ-n-1Δt,vh)+(f(chn)-f(cn),vh),∀vh∈Vh,(b)(α(chn)εhn,χh)+(ehn,divχh)=((α(cn)-α(chn))p~n,χh),∀χh∈Hh.


## 4. Several Lemmas

To obtain the optimal error estimates, we introduce the following three basic lemmas. These three lemmas are crucial to our main arguments.Lemma 4.
For any functionτ∈L2(Ω), there exists a function qτ∈H1(Ω)2, such that
(24)(a)divqτ=τ,(b)∥qτ∥1≤K∥τ∥,(c)∥qτ∥≤K∥τ∥-1,
where K is a constant.Lemma 5.
For anyτ∈L2(Ω), if (eh,εh)∈Vh×Hh is satisfying
(25)(αεh,qh)+(eh,divqh)=((α(c1)-α(c2))p~,qh),∀qh∈Hh,
then, there exists a constant K>0 such that
(26)|(τ,eh)|≤K(h∥τ∥+∥τ∥-1)(∥εh∥+∥c1-c2∥).Lemma 6.
Letη∈L2(Ω), η-=η(x-u(x)Δt),u∈W1,∞(Ω). If (eh,εh)∈Vh×Hh satisfying
(27)(αεh,qh)+(eh,divqh)=((α(c1)-α(c2))p~,qh),∀qh∈Hh,
Then,
(28)|(η-η-,eh)|≤K(h+Δt)∥η∥(∥εh∥+∥c1-c2∥).Proof.
Letτ=η-η-∈L2(Ω). By Lemma 5 we have
(29)|(η-η-,ξ)|≤K(h∥η-η-∥+∥η-η-∥-1)(∥εh∥+∥c1-c2∥).
To obtain the desired estimate, what we should deduce is to demonstrate the following:
(30)∥η-η-∥-1≤K∥η∥Δt,(31)∥η-η-∥≤K∥η∥.
Thus, consider
(32)∥η-η-∥-1=supg∈H1(Ω)[{∫Ω‍(η(x)-η(x-u(x)Δt))}1∥g∥1∫Ω‍(η(x)-η(x-u(x)Δt))×g(x)dx{1∥g∥1∫Ω‍}].
About the transformation z=H(x)=x-u(x)Δt, because of the smoothness and periodicity of u,H is a differentiable homeomorphism of Ω onto itself for Δt sufficiently small (see [9]). By using the Jacobian of this transformation and a change of variables, we have
(33)∥η-η-∥-1≤supg∈H1(Ω)[{∫Ω‍(η(x)-η(x-u(x)Δt))}1∥g∥1∫Ω‍η(x)g(x)×(1-det(DH-1(x)))dx{1∥g∥1∫Ω‍η(x)}]+supg∈H1(Ω)[{∫Ω‍(η(x)-η(x-u(x)Δt))}1∥g∥1∫Ω‍η(z)×(g(z)-g(H-1(z)))det×(DH-1(z))dz{1∥g∥1∫Ω‍η(z)}]=W1+W2.
It is easy to obtain
(34)|W1|≤Ksupg∈H1(Ω)[∥η∥∥g∥∥g∥1Δt]≤K∥η∥Δt.
Next, let l→ be the unit vector in the direction of z-H-1(z); then
(35)|W2|≤Ksupg∈H1(Ω)[1∥g∥1×∫Ω‍η(z)∫01‍∂g∂l→((1-θ)H-1(z)+θz)dθ×|z-H-1(z)|dz1∥g∥1]≤Ksupg∈H1(Ω)[∫Ω‍η(z)dz∥g∥1Δt∥g∥1]≤K∥η∥Δt.
We now obtain (30) by combining (33), (34), and (36). A similar change of variables demonstrates that
(36)∥η-∥2≤(1+KΔt)∥η∥2.
Then it is easy to obtain
(37)|(η-η-,ξ)|≤K(h+Δt)∥η∥(∥εh∥+∥c1-c2∥).


## 5. Error Estimates

Under the above assumptions aboutΩ, u, c, and f, we can derive the optimal order estimates of (ch-c) and (ph-p).Theorem 7.
Let(ch,ph),(c,p) denote the solution of (10) and (8), respectively. Suppose that Δt=O(h), we have for sufficiently small Δt>0 the following:
(38)max0≤n≤N∥(ch-c)(tn)∥≤K(∥c∥L∞(0,T;Hk+1+δ0k)+∥ct∥L2(0,T;Hk+1+δ0k))hk+1+K∥∂2c∂τ2∥L2(0,T;L2)Δt,(39)max0≤n≤N∥(ph-p)(tn)∥≤K(∥c∥L∞(0,T;Hk+2)+∥ct∥L2(0,T;Hk+1+δ0k))hk+1+K∥∂2c∂τ2∥L2(0,T;L2)Δt.Proof.
Lets first show the estimate of (38). Taking vh=ehn and χh=εhn in (23), we obtain
(40)(ehn-e-hn-1Δt,ehn)+(α(chn)εhn,εhn)=(ψ∂cn∂τ-cn-c-n-1Δt,ehn)-(ρn-ρ-n-1Δt,ehn)+((α(cn)-α(chn))p~n,εhn)+(f(chn)-f(cn),ehn)=∑i=14‍Ii.
For the first term on the right-hand side of the above equation, following the treatment manner in [8], we will bound ∥ψ(∂cn/∂τ)-(cn-c-n-1)/Δt∥ through the representation below involving an integral in the parameter τ along the tangent to the characteristics from (x-,tn-1) to (x,tn). Denote the coordinates of the point on the segment by (x(τ),t(τ)). The standard backward difference quotient error equation is given by
(41)∂cn∂t-cn-cn-1Δt=1Δt∫tn-1tn‍(t-tn-1)∂2c∂t2dt,
analogously, along the tangent to the characteristics following:
(42)ψ∂cn∂τ-cn-c-n-1Δt=1Δt∫(x-,tn-1)(x,tn)‍(x(τ)-x-)2+(t(τ)-tn-1)2∂2c∂τ2dτ.
Taking the L2(Ω) norm of this error term, we obtain
(43)∥ψ∂cn∂τ-cn-c-n-1Δt∥2≤∥ψ4∥0,∞Δt∫Ω‍∫tn-1tn‍|∂2c∂τ2(tn-tΔtx-+t-tn-1Δtx,t)|2dtdx,
to relate this to a standard norm of  ∂2c/∂τ2, consider the following transformation:
(44)S:(x,t)⟶(z,t)=(tn-tΔtx-+t-tn-1Δtx,t).
By the assumption of (3), S is invertible for sufficiently small Δt and its determinant is 1+O(Δt). Since our problem is periodic, for any fixed t∈[0,T],S is a homeomorphism of Ω onto itself, so the same is true for Ω×[tn-1,tn]. It follows that
(45)∥ψ∂cn∂τ-cn-c-n-1Δt∥2≤K∥ψ4∥0,∞∥∂2c∂τ2∥L2(tn-1,tn;L2)2Δt,
and the first term on the right-hand side of (40) is bounded by
(46)|I1|≤K∥∂2c∂τ2∥L2(tn-1,tn;L2)2Δt+K∥ehn∥2.
Concerning the second term on the right-hand side of (40), write ρn-ρ¯n-1 as the sum (ρn-ρn-1)+(ρn-1-ρ-n-1), then, by Lemma 6 and (23)(b), under the assumption Δt=O(h), we have
(47)|I2|≤KΔt∫tn-1tn‍∥ρt∥2ds+K∥ehn∥2+K∥ρ∥L∞(0,T;L2)2+a02∥εhn∥2.
Now we come to the term I3, by (4)  we have the following:
(48)|I3|≤K(∥ehn∥2+∥ρn∥2)+a02∥εhn∥2.
About the last term in the right-hand side of (40), by (4), we have
(49)|I4|≤K(∥ehn∥2+∥ρn∥2),
and this completes the treatment of the right-hand side of (40).
By using (33), the left-hand side of (40) is bounded by (1/2Δt)[(ehn,ehn)-(1+KΔt)(ehn-1,ehn-1)]+a0∥εhn∥2, combining (46)–(49) with (40) to give the following recursion relation:
(50)12Δt[(ehn,ehn)-(ehn-1,ehn-1)]+a0∥εhn∥2≤KΔt∥∂2c∂τ2∥L2(tn-1,tn;L2)2+KΔt∫tn-1tn‍∥ρt∥2ds+K∥ρ∥L∞(0,T;L2)2+K∥ehn∥2+K∥ehn-1∥2+a0∥εhn∥2.
If (50) is multiplied by 2Δt and summed in time, and if it is noted that (10)(c) implies that eh0=0, then, it follows that
(51)∥ehn∥2≤K(Δt2∥∂2c∂τ2∥L2(0,T;L2)2+∥ρt∥L2(0,T;L2)2+∥ρ∥L∞(0,T;L2)2{∥∂2c∂τ2∥})+KΔt∑i=1n‍∥ehi∥2.
For sufficiently small Δt>0, by Gronwall's lemma, we obtain
(52)∥ehn∥≤K(∥ρt∥L2(0,T;L2)+∥ρ∥L∞(0,T;L2))+KΔt∥∂2c∂τ2∥L2(0,T;L2).
Recall that the error ch-c=eh+ρ, (19), (21), and (53) together imply that
(53)max0≤n≤N∥(ch-c)(tn)∥≤K(∥c∥L∞(0,T;Hk+1+δ0k)+∥ct∥L2(0,T;Hk+1+δ0k))hk+1+K∥∂2c∂τ2∥L2(0,T;L2)Δt.
This completes the first part of our conclusion.
In order to show the estimate (39), by (23)(b), we get
(54)(α(chn)εhn-α(chn-1)εhn-1Δt,χh)+(ehn-ehn-1Δt,divχh)=((α(cn)-α(chn))p~n-(α(cn-1)-α(chn-1))p~n-1Δt,χh).
Choosing vh=(ehn-ehn-1)/Δt and χh=εhn, respectively, in (23)(a) and (54) and adding them to obtain
(55)(ehn-ehn-1Δt,ehn-ehn-1Δt)+(α(chn)εhn-α(chn-1)εhn-1Δt,εhn)=(ψ∂cn∂τ-cn-c-n-1Δt,ehn-ehn-1Δt)-(ρn-ρn-1Δt,ehn-ehn-1Δt)-(ρn-1-ρ-n-1Δt,ehn-ehn-1Δt)+((α(cn)-α(chn))p~n-(α(cn-1)-α(chn-1))p~n-1Δt,εhn)-(ehn-1-e-hn-1Δt,ehn-ehn-1Δt)+(f(chn)-f(cn),ehn-ehn-1Δt).
The left-hand side satisfies the inequality:
(56)(ehn-ehn-1Δt,ehn-ehn-1Δt)+(α(chn)εhn-α(chn-1)εhn-1Δt,εhn)≥12Δt((α(chn)εhn,εhn)-(α(chn-1)εhn-1,εhn-1))+12Δt((α(chn-1)-α(chn))εhn-1,εhn-1)+1Δt((α(chn)-α(chn-1))εhn-1,εhn)∥ehn-ehn-1Δt∥2.
Substitute it into (55), we have
(57)∥ehn-ehn-1Δt∥2+12Δt((α(chn)εhn,εhn)-(α(chn-1)εhn-1,εhn-1))=(ψ∂cn∂τ-cn-c-n-1Δt,ehn-ehn-1Δt)+((α(cn)-α(chn))p~n-(α(cn-1)-α(chn-1))p~n-1Δt,εhn)+12Δt((α(chn-1)-α(chn))εhn-1,2εhn-εhn-1)+(f(chn)-f(cn),ehn-ehn-1Δt)-(ρn-ρn-1Δt,ehn-ehn-1Δt)-(ρn-1-ρ-n-1Δt,ehn-ehn-1Δt)-(ehn-1-e-hn-1Δt,ehn-ehn-1Δt)=∑i=17‍Ji.
Bound J1–J5 on the right-hand side of (56) one by one, we have
(58)|J1|≤K∥∂2c∂τ2∥L2(tn-1,tn;L2)Δt+16∥ehn-ehn-1Δt∥2,|J2|≤K(∥ehn∥2+∥ρn∥2+∥ehn-1∥2+∥ρn-1∥2)+16∥ehn-ehn-1Δt∥2+KΔt∫tn-1tn‍∥ρt∥2dt+K∥εhn∥2,|J3|≤K(∥εhn-1∥0,∞2+1)(∥εhn∥2+∥εhn-1∥2)+16∥ehn-ehn-1Δt∥2.
Concerning J4 and J5, we can easily get
(59)|J4|≤K(∥ehn∥2+∥ρn∥2)+16∥ehn-ehn-1Δt∥2,|J5|≤KΔt∫tn-1tn‍∥ρt∥2dt+16∥ehn-ehn-1Δt∥2.
Substituting (58)-(59) into (57) to obtain
(60)∥ehn-ehn-1Δt∥2+12Δt((α(chn)εhn,εhn)-(α(chn-1)εhn-1,εhn-1))≤K∥∂2c∂τ2∥L2(tn-1,tn;L2)2Δt+K(∥ehn∥2+∥ρn∥2+∥ehn-1∥2+∥ρn-1∥2)+56∥ehn-ehn-1Δt∥2+KΔt∫tn-1tn‍∥ρt∥2dt+K(1+∥εhn-1∥0,∞2)(∥εhn∥2+∥εhn-1∥2)-(ρn-1-ρ-n-1Δt,ehn-ehn-1Δt)-(ehn-1-e-hn-1Δt,ehn-ehn-1Δt).
Noting that eh0=0, we have ∥εh0∥≤K∥ρ0∥ from (23)(b). To multiply (60) by 2Δt and add it in time, we obtain
(61)Δt3∑i=1n‍∥ehn-ehn-1Δt∥2+a0∥εhn∥2≤K∥∂2c∂τ2∥L2(0,T;L2)2Δt2+KΔt(1+max0≤i≤n-1∥εhi∥0,∞)∑i=1n‍∥εhi∥2+K(∥ρ∥L∞(0,T;L2)2+∥ρt∥L2(0,T;L2)2+max1≤i≤n∥ehi∥2+∥ρ0∥2)-2(∑i=1n‍(ρi-1-ρ-i-1Δt,ehi-ehi-1)+∑i=1n‍(ehi-1-e-hi-1Δt,ehi-ehi-1)).
By (23)(b) and Lemma 6, we have
(62)∑i=1n‍(ρi-1-ρ-i-1Δt,ehi-ehi-1)≤K(max1≤i≤n∥ehi∥2+∥ρ∥L∞(0,T;L2)2+∥ρt∥L2(0,T;L2)2)+a04∥εhn∥2+KΔt∑i=1n-1‍∥εhi∥2,
where the assumption Δt=O(h) is used.
Similarly, we get(63)∑i=1n‍(ehi-1-e-hi-1Δt,ehi-ehi-1)≤K(max1≤i≤n∥ehi∥2+∥ρ∥L∞(0,T;L2)2)+a04∥εhn∥2+KΔt∑i=1n-1‍∥εhi∥2+Δt3∑i=1n-1‍∥ehi-ehi-1Δt∥2.
Substituting (62) and (63) into (61), we get
(64)Δt∑i=1n‍∥ehi-ehi-1Δt∥2+a0∥εhn∥2≤KΔt(1+max0≤i≤n∥εhi∥0,∞)∑i=1n‍∥εhi∥2+KΔt2∥∂2c∂τ2∥L2(0,T;L2)2+K(∥ρt∥L2(0,T;L2)2+∥ρ∥L∞(0,T;L2)2+max1≤i≤n∥ehi∥2).
We see that an application of Gronwall's lemma would complete the second part of our argument if we did not have the termmax0≤i≤n-1∥εhi∥0,∞ on the right hand side of (64). We will use an induction argument to prove the hypothesis as follows:
(65)∥εhi∥0,∞≤M1,i=1,2,…,N.
Since ∥εh0∥≤K∥ρ0∥, by the inverse property and (19), it can be easily verified that ∥εh0∥0,∞≤M1. Using (19), (21), (52), and (64) we see that (65) clearly holds for i=1 if Δt is sufficiently small. Assume that (65) holds for i=0,1,…,n-1, using (64), we obtain
(66)∥εhn∥2≤KΔt∑i=1n‍∥εhi∥2+K∥∂2c∂τ2∥L2(0,T;L2)2Δt2+K(max1≤i≤n∥ehi∥2+∥ρt∥L2(0,T;L2)2+∥ρ∥L∞(0,T;L2)2),
Then, an application of Gronwall's lemma gives us an estimate which together with (19), (21), and (52) and the inverse property of Hh, show that the induction hypothesis (65) holds for i=n for sufficiently small Δt. This completes the induction argument.
Finally by Gronwall's lemma, it follows from (64) that
(67)∥εhn∥≤K(∥ρt∥L2(0,T;L2)+∥ρ∥L∞(0,T;L2)+max1≤i≤n∥ehi∥)+K∥∂2c∂τ2∥L2(0,T;L2)Δt.
Note that ph-p=εh+σ, by (52) and (19)–(21), we obtain
(68)∥phn-pn∥≤K(∥ct∥L2(0,T;Hk+1+δ0k)+∥c∥L∞(0,T;Hk+2))hk+1+K∥∂2c∂τ2∥L2(0,T;L2)Δt.
This also shows that
(69)max0≤n≤N∥phn-pn∥≤K(∥ct∥L2(0,T;Hk+1+δ0k)+∥c∥L∞(0,T;Hk+2))hk+1+K∥∂2c∂τ2∥L2(0,T;L2)Δt.Note that our results in this paper do not cover the cases of nonlinear convection-dominated systems, which are of importance in some applications, particularly in reservoir simulations.



---

*Source: 264025-2013-09-09.xml*