# Parallel Nonoverlapping DDM Combined with the Characteristic Method for Incompressible Miscible Displacements in Porous Media

Authors: Keying Ma; Tongjun Sun
Journal: Advances in Numerical Analysis
 (2013)

---

## Abstract

Two types of approximation schemes are established for incompressible miscible
displacements in porous media. First, standard mixed finite element method is used
to approximate the velocity and pressure. And then parallel non-overlapping domain
decomposition methods combined with the characteristics method are presented for
the concentration. These methods use the characteristic method to handle the material
derivative term of the concentration equation in the subdomains and explicit flux
calculations on the interdomain boundaries by integral mean method or extrapolation
method to predict the inner-boundary conditions. Thus, the velocity and pressure can
be approximated simultaneously, and the parallelism can be achieved for the concentration
equation. The explicit nature of the flux prediction induces a time step limitation
that is necessary to preserve stability. These schemes hold the advantages of nonoverlapping
domain decomposition methods and the characteristic method. Optimal
error estimates inL2-norm are derived for these two schemes, respectively.


---

## Body

## 1. Introduction

The two-phase fluid displacements in porous media is one of the most important basic problems in the oil reservoir numerical simulation. It is governed by a nonlinear coupled system of partial differential equations with initial and boundary values. In this paper, we will consider the following incompressible miscible case: the pressure is governed by an elliptic equation and the concentration is governed by a convection-diffusion equation [1–5]. (1a)-∇·(k(x)μ(c)(∇p-γ(c)∇d(x)))≡∇·u=q,x∈Ω,t∈J,(1b)ϕ∂c∂t+u·∇c-∇·(D(u)∇c)=(c~-c)q~,x∈Ω,t∈J,(1c)u·ν=(D(u)∇c)·ν=0,x∈∂Ω,t∈J,(1d)c(x,0)=c0(x),x∈Ω,whereΩ is a bounded domain in R2, J=(0,T], and q~=max{q,0} is nonzero at injection wells only. The variables in (1a)–(1d) are the pressure p(x,t) in the fluid mixture, the Darcy velocity u=(u1,u2)′, and the relative concentration c(x,t) of the injected fluid. The ν is the unit outward normal vector on boundary ∂Ω.The coefficients and data in (1a)–(1d) are k(x): the permeability of the porous media; μ(c), the viscosity of the fluid mixture; q(x,t): representing flow rates at wells; γ(c) and d(x), the gravity coefficient and vertical coordinate; ϕ(x): the porosity of the rock; c~(x,t), the injected concentration at injection wells (q>0) and the resident concentration at production wells (q<0). Here, D(u) is a tensor 2×2 matrix and generally has the form
(2)D(u)=ϕ(x){dmI+|u|(dlE(u)+dtE⊥(u))},
where 2×2 matrix E=(eij) satisfies eij=uiuj/|u|2, E⊥=I-E, dm is the molecular diffusivity, and dl, dt are longitudinal and transverse dispersivities, respectively. Furthermore, a compatibility condition ∫Ω‍q(x,t)dx=0 must be imposed to determine the pressure.The pressure equation is elliptic and easily handled by standard mixed finite element method, which has been proven to be an effective numerical method for solving fluid problems. It has an advantage to approximate the unknown variable and its diffusive flux simultaneously. There are many research articles on this method [6–9]. The concentration equation is parabolic and normally convection dominated. It is well known that standard Galerkin scheme applied to the convection dominated problems does not work well and produces excessive numerical diffusion or nonphysical oscillation. A variety of numerical techniques have been introduced to obtain better approximations for (1a)–(1d), such as characteristic finite difference method [10], characteristic finite element method [11], the modified of characteristic finite element method (MMOC-Galerkin) [12], and the Eulerian-Lagrangian localized adjoint method (ELLAM) [13].It is well known that parallel algorithms, based upon overlapping or non-overlapping domain decompositions, are effective ways to solve the large scale of PDE systems for most practical problems in engineering (e.g., see [14–17]). We have presented parallel Galerkin domain decomposition procedures for parabolic equation [18–20]. These procedures use implicit Galerkin methods in the subdomains and explicit flux calculations on the interdomain boundaries by integral mean method or extrapolation-integral mean method. Some constraints for time step are still needed for these procedures to preserve stability, but less severe than that for fully explicit methods. With respect to the accuracy order of h, L2-norm error estimates are optimal for higher-order finite element spaces and almost optimal for linear finite element space in twodimensional domain. Compared with Dawson-Dupont's schemes [21], these L2-norm error estimates avoid the loss of H-1/2 factor.We also have considered using the procedures in [18] for wave equation [22] and convection-diffusion equation [23]. This paper is one of our sequent research papers. The main purpose is to use parallel Galerkin domain decomposition procedures in [18] combined with the characteristic method for the concentration equation of incompressible miscible displacements in porous media. This paper is organized as follows. In Section 2, we first present mixed finite element method for the velocity and pressure and then formulate parallel non-overlapping domain decomposition methods combined with the characteristics method for the concentration. We establish the combined approximation schemes. Auxiliary lemmas are listed in Section 3, which show some properties of finite element spaces and projections. In Section 4, we derive the optimal-order L2-norm error estimates. Finally in Section 5, we extend the consideration for another approximate scheme by using extrapolation method. It is worthwhile to specially emphasize that the research of this paper is creative. No former researchers discussed it. These schemes not only hold the advantages of non-overlapping domain decomposition methods, but also hold the advantages of the characteristic method.


## 2.  Formulation of the Methods

In this paper, we adopt notations and norms of usual Sobolev spaces [3]:
(3)Hm(Ω)={f:∂|α|f∂xα∈L2(Ω)for|α|≤m},∥f∥m2=∑|α|≤m‍∥∂|α|f∂xα∥2,m≥0,W∞m(Ω)={f:∂|α|f∂xα∈L∞(Ω)for|α|≤m},∥f∥W∞m=max|α|≤m∥∂|α|f∂xα∥∞,m≥0.
In particular, H0(Ω) and W∞0(Ω)=L∞(Ω). The inner product on L2(Ω) is denoted by (·,·).We also use the following spaces that incorporate time dependence. Let[a,b]⊂J and let X be any of the spaces just defined. For f(x,t) suitably smooth on Ω×[a,b], we let
(4)Hm(a,b;X)={f:∫ab‍∥∂αf(·,t)∂tα∥Xdt<∞,forα≤m},∥f∥Hm(a,b;X)=(∑α=0m‍∫ab‍∥∂αf(·,t)∂tα∥Xdt)1/2,m≥0.
Similarly, W∞m(a,b;X) and the norm ∥f∥W∞m(a,b;X) are defined. If [a,b]=J-, we simplify our notation and write L∞(W∞1) for L∞(0,T;W∞1(Ω)), and so forth.Iff=(f1,f2) is a vector function, we note that f∈X if f1∈X and f2∈X. We also use the vector-function spaces and norms:
(5)Hm(div;Ω)={f:f1,f2,∇·f∈Hm(Ω)},∥f∥Hm(div;Ω)2=∥f1∥m2+∥f2∥m2+∥∇·f∥m2,m≥0,H(div;Ω)=H0(div;Ω).We need some assumptions. The regularity assumptions on the solution of (1a)–(1d) are noted by [4]:
(6)(P):{c∈L∞(H2)∩H1(H1)∩L∞(W∞1)∩H2(L2),p∈L∞(H1),u∈L∞(H1(div))∩L∞(W∞1)∩W∞1(L∞)∩H2(L2).We also require the following assumptions on the coefficients in (1a)–(1d) [4]. Let a0, a1, ϕ0, ϕ1, d0, d1, K1, and K2 be positive constants such that
(7)(Q):{0<a0≤k(x)μ(c)≤a1,0<ϕ0≤ϕ(x)≤ϕ1,|∂(k/μ)∂c(x,c)|+|∂γ∂c(x,c)|+|∇ϕ(x)|+|q~(x,t)|+|∂q~∂t(x,t)|≤K1,|∂D∂u(u)|+|∂D-1∂u(u)|≤K2,d0∥ξ∥2≤∑i,j=12D(u)ξiξj≤d1∥ξ∥2,∀ξ∈ℛ2.
Further assumptions will be made in individual theorems as necessary.For convenience, we assume that (1a)–(1d) is Ω periodic (see [4]); that is, all functions will be assumed to be spatially Ω-periodic throughout the rest of this paper. This is physically reasonable, because no-flow conditions (1c) are generally treated by reflection, and in general, interior flow patterns are much more important than boundary effects in reservoir simulation. Thus, the boundary conditions (1c) can be dropped.Throughout the analysis,C and Ki(i=3,…,6) will denote generic positive constants, independent of hc, hp, Δtc, and Δtp, but possibly depending on constants in (Q), norms in (P). Similarly, ε will denote a generic small positive constant.


### 2.1. A Mixed Finite Element Method for the Pressure and Velocity

Let(8)H={χ∈H(div;Ω)∣χ·ν=0on∂Ω},(9)W=L2(Ω){w≡constantonΩ}.As in [4], the pressure equation (1a) is equal to the following saddle-point problem of finding a map (u,p):J→H×W such that
(10)(a)A(c;u,v)+B(v,p)=(γ(c)∇d,v),∀v∈H(div;Ω),(b)B(u,w)=-(q,w),∀w∈L2(Ω),
where the bilinear forms
(11)A(c;u,v)=(μ(c)k(x)u,v),B(v,p)=-(∇·v,p),B(u,w)=-(∇·u,w).Forhp>0, we discrete (10) in space on a quasi-uniform mesh 𝒯hp of Ω with diameter of element ≤hp. Let Vhp×Whp⊂H×W be Raviart-Thomas spaces [6, 7] of index l=0 for this mesh.The mixed method for pressure and velocity, given a concentration approximationC at a time t∈J, consists of U∈Vhp and P∈Whp such that
(12)A(C;U,v)+B(v,P)=(γ(c)∇d,v),∀v∈Vhp,B(U,w)=-(q,w),∀w∈Whp.
Existence and uniqueness of U and P are proved in [1], based on ideas of [24].


### 2.2. A Characteristic DDM for the Concentration

In this section, we assumeu in (1a)–(1d) is given. Define
(13)M={v∈L2(Ω)∣vispiecewiseconstant},MisdenseinL2(Ω).Let(14)ψ(x)=[|u(x)|2+(ϕ(x))2]1/2=[(u1(x))2+(u2(x))2+(ϕ(x))2]1/2,
and let the characteristic direction associated with the operator ϕ∂c/∂t+u·∇c be denoted by τ, where
(15)ψ(x)∂c∂τ=ϕ(x)∂c∂t+u·∇c.The weak form of (1a)–(1d) is finding a map c:J→L2(Ω) such that
(16)(a)(ψ∂c∂τ,φ)+((D(u)∇c),∇φ)=((c~-c)q~,φ),∀φ∈L2(Ω),(b)c(x,0)=c0(x),∀x∈Ω.PartJ into 0=t0<t1<⋯<tN=T, with Δtcn=tn-tn-1. The analysis is valid for variable time steps, but we drop the superscript from Δtc for convenience. For functions f on Ω×J, we write fn(x) for f(x,tn). Approximate (∂cn/∂τ)(x)=(∂c/∂τ)(x,tn) by a backward difference quotient in the τ-direction,
(17)∂cn∂τ(x)≃cn(x)-cn-1(x-(u(x)/ϕ(x))Δtc)Δtc1+|u(x)|2/(ϕ(x))2.If we letx-=x-(u(x)/ϕ(x))Δtc and f-(x)=f(x-), then we get
(18)ψ∂cn∂τ≃ϕcn-c-n-1Δtc.
Since the problem is Ω-periodic [4], c-n-1 is always defined and the tangent to the characteristic (i.e., the τ segment) cannot cross a boundary to an undefined location. The difference quotient relates the concentration at a given x at time tn to the concentration that would flow to x from time tn-1 if the problem were purely hyperbolic.The time difference (18) will be combined with a characteristic DDM in the space variables. We recall the domain decomposition procedures in [18] here. For simplicity and without losing generality, we only discuss the case of two subdomains. But the algorithms and theories can be extended to the case of many subdomains. Divide Ω into two subdomains  Ωj (j=1,2) by an interdomain boundary Γ, which is a surface of dimension d-1, see Figure 1. We denote by Γj=∂Ωj⋂‍∂Ω the part of the boundary of the subdomains which coincides with ∂Ω. Denote the unit vector normal to Γ as νΓ, which points from Ω1 toward Ω2.Figure 1
The domainΩ with the interdomain boundary Γ.Let𝒯j,hc be quasi-uniform partitions of Ωj (j=1,2) and 𝒯hc=𝒯1,hc⋃‍𝒯2,hc. Here, h denotes the maximal element diameter of 𝒯hc. We construct the finite element space ℳhc on 𝒯hc which satisfies the following condition (I)(1)
Forj=1,2, let ℳj,hc be a finite element subspace of H1(Ωj), and let ℳhc⊂L2(Ω) such that if v∈ℳhc, then v∣Ωj∈ℳj,hc.(2)
Forj=1,2, Pr(Ωj)⊂ℳj,hc, where Pr(Ωj) is a polynomial space of degree at most r.(3)
Forj=1,2, h∈(0,1], the integer k≥1, and u∈Hk(Ωj), there exists a positive constant C independent of h such that
(19)infv∈ℳj,hc∥u-v∥Hs(Ωj)≤Chσ∥u∥Hk(Ωj),0≤s≤1,where σ=min(r+1-s,k-s).From the definition above, we note that functionsv in ℳhc have a well-defined jump [v] on Γ:
(20)[v](x)=v(x+)-v(x-),∀xonΓ,
where v(x±):=limλ→0±v(x+λνΓ).To construct parallel algorithm, for a small constant0<H<min{diam(Ω1),diam(Ω2)}, we introduce an integral mean value of a given function V∈L2(Ω) on the interdomain boundary Γ as
(21)V¯H(x)=12H∫-HHV(x+λνΓ)dλ,∀xonΓ.
Furthermore, we define the extrapolation of  V¯H(x) on Γ as
(22)V^H(x)=4V¯H/2(x)-V¯H(x)3,∀xonΓ.Generally, near the intersection of boundary∂Ω and inner-boundary Γ, the value of V outside Ω is needed for V¯H and V^H. If x∉Ω, let x~∈Ω denote the symmetric point of x with respect to ∂Ω. For a given function u∈L2(Ω), we define
(23)Eu(x)={u(x),ifx∈Ω,u(x~),ifx∉Ω.
By (23), we know V¯H and V^H have the values on a strip domain G={y∣y=x+λνΓ,λ∈[-H,H],∀xonΓ}, see Figure 2.Figure 2
The strip domainG with width 2H.Now, we present the non-overlapping characteristic DDM for the concentration equation:(24)(ϕCn-C^n-1Δtc,v)+(D∇Cn,∇v)+(qnCn,v)+(C¯μ,Hn-1,[v])Γ+(v¯μ,H,[Cn-1])Γ+d1KH-1([Cn-1],[v])Γ=(qnc~n,v),∀v∈ℳhc,
where
(25)C^n-1(x)=Cn-1(x-Un(x)ϕ(x)Δtc),Cμn-1=(D∇ECn-1)·νΓ.
The coefficient K will be given in Lemma 3.


### 2.3. The Combined Approximate Schemes

We now present our sequential time-stepping procedure that combines (24) and (12). As in [4], let us part J into pressure time steps 0=t0<t1<⋯<tM=T, with Δtpm=tm-tm-1. Each pressure step is also a concentration step; that is, for each m there exists n such that tm=tn; in general, Δtp>Δtc. We may vary Δtp, but except for Δtp1 we drop the superscript. For functions f on Ω×J, we write fm(x) for f(x,tm); thus, subscripts refer to pressure steps and superscripts to concentration steps.If concentration steptn relates to pressure steps by tm-1<tn≤tm, we require a velocity approximation for (25) based on Um-1 and earlier values. If m≥2, take the linear extrapolation of Um-1 and Um-2 defined by
(26)FUn=(1+tn-tm-1tm-1-tm-2)Um-1-tn-tm-1tm-1-tm-2Um-2.
If m=1, set
(27)FUn=U0.
We retain the superscript on Δtp1 because FUn is first order correct in time during the first pressure step and second order during later steps.The combined time-stepping procedure is finding a mapC:{t0,t1,…,tN}→ℳhc and a map (U,P):{t0,t1,…,tM}→Vhp×Whp such that(28a)C0=C~0,(28b)A(Cm;Um,z)+B(z,Pm)=(γ(Cm)∇d,z),∀z∈Vhp,(28c)B(Um,w)=-(qm,w),∀w∈Whp,m≥0,(28d)(ϕCn-Cˇn-1Δtc,v)+(D∇Cn,v)+(qnCn,v)+(C¯μ,Hn-1,[v])Γ+(v¯μ,H,[Cn-1])Γ+d1KH-1([Cn-1],[v])Γ=(qnc~n,v),∀v∈ℳhc,n≥1,
where (29)Cˇn-1(x)=Cn-1(xˇ)=Cn-1(x-FUn(x)ϕ(x)Δtc),C~0 is given by (45).The steps of the above calculation are as follows:Step  1.C0 known →  solve (U0,P0) by  (28b) and  (28c);Step  2. On each domain, use  (28d) to parallelly compute {Cj}j=1n1 until tn1=t1;Step  3. Then by  (28b) and  (28c) for (U1,P1);Step  4. Calculate the approximations in turn analogically to get the pressure, velocity, and concentration at other time-step, respectively.The convergence analysis will make use of an analogue ofx^ defined for the exact velocity un. If f is a function on Ω, set
(30)fˇ(x)=f(xˇ)=f(x-Fun(x)ϕ(x)Δtc).
The time step tn will be clear from the context.Remark 1.
In the scheme (28a)–(28d), the numerical flux on Γ is computed explicitly from Cn-1, so that Cn can be computed on Ω1 and Ω2 fully parallel once Cn-1 has been got.


## 2.1. A Mixed Finite Element Method for the Pressure and Velocity

Let(8)H={χ∈H(div;Ω)∣χ·ν=0on∂Ω},(9)W=L2(Ω){w≡constantonΩ}.As in [4], the pressure equation (1a) is equal to the following saddle-point problem of finding a map (u,p):J→H×W such that
(10)(a)A(c;u,v)+B(v,p)=(γ(c)∇d,v),∀v∈H(div;Ω),(b)B(u,w)=-(q,w),∀w∈L2(Ω),
where the bilinear forms
(11)A(c;u,v)=(μ(c)k(x)u,v),B(v,p)=-(∇·v,p),B(u,w)=-(∇·u,w).Forhp>0, we discrete (10) in space on a quasi-uniform mesh 𝒯hp of Ω with diameter of element ≤hp. Let Vhp×Whp⊂H×W be Raviart-Thomas spaces [6, 7] of index l=0 for this mesh.The mixed method for pressure and velocity, given a concentration approximationC at a time t∈J, consists of U∈Vhp and P∈Whp such that
(12)A(C;U,v)+B(v,P)=(γ(c)∇d,v),∀v∈Vhp,B(U,w)=-(q,w),∀w∈Whp.
Existence and uniqueness of U and P are proved in [1], based on ideas of [24].


## 2.2. A Characteristic DDM for the Concentration

In this section, we assumeu in (1a)–(1d) is given. Define
(13)M={v∈L2(Ω)∣vispiecewiseconstant},MisdenseinL2(Ω).Let(14)ψ(x)=[|u(x)|2+(ϕ(x))2]1/2=[(u1(x))2+(u2(x))2+(ϕ(x))2]1/2,
and let the characteristic direction associated with the operator ϕ∂c/∂t+u·∇c be denoted by τ, where
(15)ψ(x)∂c∂τ=ϕ(x)∂c∂t+u·∇c.The weak form of (1a)–(1d) is finding a map c:J→L2(Ω) such that
(16)(a)(ψ∂c∂τ,φ)+((D(u)∇c),∇φ)=((c~-c)q~,φ),∀φ∈L2(Ω),(b)c(x,0)=c0(x),∀x∈Ω.PartJ into 0=t0<t1<⋯<tN=T, with Δtcn=tn-tn-1. The analysis is valid for variable time steps, but we drop the superscript from Δtc for convenience. For functions f on Ω×J, we write fn(x) for f(x,tn). Approximate (∂cn/∂τ)(x)=(∂c/∂τ)(x,tn) by a backward difference quotient in the τ-direction,
(17)∂cn∂τ(x)≃cn(x)-cn-1(x-(u(x)/ϕ(x))Δtc)Δtc1+|u(x)|2/(ϕ(x))2.If we letx-=x-(u(x)/ϕ(x))Δtc and f-(x)=f(x-), then we get
(18)ψ∂cn∂τ≃ϕcn-c-n-1Δtc.
Since the problem is Ω-periodic [4], c-n-1 is always defined and the tangent to the characteristic (i.e., the τ segment) cannot cross a boundary to an undefined location. The difference quotient relates the concentration at a given x at time tn to the concentration that would flow to x from time tn-1 if the problem were purely hyperbolic.The time difference (18) will be combined with a characteristic DDM in the space variables. We recall the domain decomposition procedures in [18] here. For simplicity and without losing generality, we only discuss the case of two subdomains. But the algorithms and theories can be extended to the case of many subdomains. Divide Ω into two subdomains  Ωj (j=1,2) by an interdomain boundary Γ, which is a surface of dimension d-1, see Figure 1. We denote by Γj=∂Ωj⋂‍∂Ω the part of the boundary of the subdomains which coincides with ∂Ω. Denote the unit vector normal to Γ as νΓ, which points from Ω1 toward Ω2.Figure 1
The domainΩ with the interdomain boundary Γ.Let𝒯j,hc be quasi-uniform partitions of Ωj (j=1,2) and 𝒯hc=𝒯1,hc⋃‍𝒯2,hc. Here, h denotes the maximal element diameter of 𝒯hc. We construct the finite element space ℳhc on 𝒯hc which satisfies the following condition (I)(1)
Forj=1,2, let ℳj,hc be a finite element subspace of H1(Ωj), and let ℳhc⊂L2(Ω) such that if v∈ℳhc, then v∣Ωj∈ℳj,hc.(2)
Forj=1,2, Pr(Ωj)⊂ℳj,hc, where Pr(Ωj) is a polynomial space of degree at most r.(3)
Forj=1,2, h∈(0,1], the integer k≥1, and u∈Hk(Ωj), there exists a positive constant C independent of h such that
(19)infv∈ℳj,hc∥u-v∥Hs(Ωj)≤Chσ∥u∥Hk(Ωj),0≤s≤1,where σ=min(r+1-s,k-s).From the definition above, we note that functionsv in ℳhc have a well-defined jump [v] on Γ:
(20)[v](x)=v(x+)-v(x-),∀xonΓ,
where v(x±):=limλ→0±v(x+λνΓ).To construct parallel algorithm, for a small constant0<H<min{diam(Ω1),diam(Ω2)}, we introduce an integral mean value of a given function V∈L2(Ω) on the interdomain boundary Γ as
(21)V¯H(x)=12H∫-HHV(x+λνΓ)dλ,∀xonΓ.
Furthermore, we define the extrapolation of  V¯H(x) on Γ as
(22)V^H(x)=4V¯H/2(x)-V¯H(x)3,∀xonΓ.Generally, near the intersection of boundary∂Ω and inner-boundary Γ, the value of V outside Ω is needed for V¯H and V^H. If x∉Ω, let x~∈Ω denote the symmetric point of x with respect to ∂Ω. For a given function u∈L2(Ω), we define
(23)Eu(x)={u(x),ifx∈Ω,u(x~),ifx∉Ω.
By (23), we know V¯H and V^H have the values on a strip domain G={y∣y=x+λνΓ,λ∈[-H,H],∀xonΓ}, see Figure 2.Figure 2
The strip domainG with width 2H.Now, we present the non-overlapping characteristic DDM for the concentration equation:(24)(ϕCn-C^n-1Δtc,v)+(D∇Cn,∇v)+(qnCn,v)+(C¯μ,Hn-1,[v])Γ+(v¯μ,H,[Cn-1])Γ+d1KH-1([Cn-1],[v])Γ=(qnc~n,v),∀v∈ℳhc,
where
(25)C^n-1(x)=Cn-1(x-Un(x)ϕ(x)Δtc),Cμn-1=(D∇ECn-1)·νΓ.
The coefficient K will be given in Lemma 3.


## 2.3. The Combined Approximate Schemes

We now present our sequential time-stepping procedure that combines (24) and (12). As in [4], let us part J into pressure time steps 0=t0<t1<⋯<tM=T, with Δtpm=tm-tm-1. Each pressure step is also a concentration step; that is, for each m there exists n such that tm=tn; in general, Δtp>Δtc. We may vary Δtp, but except for Δtp1 we drop the superscript. For functions f on Ω×J, we write fm(x) for f(x,tm); thus, subscripts refer to pressure steps and superscripts to concentration steps.If concentration steptn relates to pressure steps by tm-1<tn≤tm, we require a velocity approximation for (25) based on Um-1 and earlier values. If m≥2, take the linear extrapolation of Um-1 and Um-2 defined by
(26)FUn=(1+tn-tm-1tm-1-tm-2)Um-1-tn-tm-1tm-1-tm-2Um-2.
If m=1, set
(27)FUn=U0.
We retain the superscript on Δtp1 because FUn is first order correct in time during the first pressure step and second order during later steps.The combined time-stepping procedure is finding a mapC:{t0,t1,…,tN}→ℳhc and a map (U,P):{t0,t1,…,tM}→Vhp×Whp such that(28a)C0=C~0,(28b)A(Cm;Um,z)+B(z,Pm)=(γ(Cm)∇d,z),∀z∈Vhp,(28c)B(Um,w)=-(qm,w),∀w∈Whp,m≥0,(28d)(ϕCn-Cˇn-1Δtc,v)+(D∇Cn,v)+(qnCn,v)+(C¯μ,Hn-1,[v])Γ+(v¯μ,H,[Cn-1])Γ+d1KH-1([Cn-1],[v])Γ=(qnc~n,v),∀v∈ℳhc,n≥1,
where (29)Cˇn-1(x)=Cn-1(xˇ)=Cn-1(x-FUn(x)ϕ(x)Δtc),C~0 is given by (45).The steps of the above calculation are as follows:Step  1.C0 known →  solve (U0,P0) by  (28b) and  (28c);Step  2. On each domain, use  (28d) to parallelly compute {Cj}j=1n1 until tn1=t1;Step  3. Then by  (28b) and  (28c) for (U1,P1);Step  4. Calculate the approximations in turn analogically to get the pressure, velocity, and concentration at other time-step, respectively.The convergence analysis will make use of an analogue ofx^ defined for the exact velocity un. If f is a function on Ω, set
(30)fˇ(x)=f(xˇ)=f(x-Fun(x)ϕ(x)Δtc).
The time step tn will be clear from the context.Remark 1.
In the scheme (28a)–(28d), the numerical flux on Γ is computed explicitly from Cn-1, so that Cn can be computed on Ω1 and Ω2 fully parallel once Cn-1 has been got.


## 3. Auxiliary Lemmas

We adopt some auxiliary lemmas about the finite element spaces, which will be used in the next section.


### 3.1.   For the Pressure and Velocity

The Raviart-Thomas spacesVhp×Whp in Section 2.1 possess the following approximation and inverse properties [7]:
(31)(Ap):{infv∈Vhp∥f-v∥≤K3hp∥f∥1,infv∈Vhp∥f-v∥H(div)≤K3hp∥f∥H1(div),infw∈Whp∥g-w∥≤K3hp∥g∥1,(Ip):{∥v∥L∞≤K3hp-1∥v∥,∀v∈Vhp,∥v∥W∞1(Tp)≤K3hp-1∥v∥L∞(Tp),∀v∈Vhp,
where K3 is a positive constant independent of hp and Tp is an element of the mesh 𝒯hp.Define the map(U~,P~):J→Vhp×Whp by (see [7])
(32)A(c;U~,v)+B(v,P~)=(γ(c)∇d,v),∀v∈Vhp,B(U~,w)=-(q,w),∀w∈Whp,
where c is the exact solution of (1a)–(1d).By arguments in [1, 24], the map (U~,P~) exists and (Ap) implies that
(33)∥u-U~∥L∞(H(div))+∥p-P~∥L∞(L2)≤K4{infv∈Vhp∥u-v∥H(div)+infw∈Whp∥p-w∥}≤K4(∥u∥L∞(H1(div)),∥p∥L∞(H1))hp.
The positive constant K4 depends on constant in (Q) but independent of hp, u, p, and c.The estimate (33) and (Ip) imply that
(34)∥U~∥L∞(L∞)≤K5.
As in [1], comparison of (32) and (12) implies that
(35)∥U-U~∥H(div)+∥P-P~∥≤K6(1+∥U~∥L∞)∥c-C∥.
The estimates (33) and (35) will handle the coupling of concentration and velocity errors in the convergence analysis.


### 3.2. For the Concentration

In this section, we adopt some following lemmas [18].Lemma 2.
For smooth enough functionW, there holds estimates:
(36)∥W¯H-W∥L2(Γ)≤2H∥∇W∥L2(Ω),∥W¯H-W∥L∞(Γ)≤CH2∥W∥W2,∞(Ω),W(x)-W¯H(x)=-16H2WνΓ2(x)-1120H4WνΓ4(x)+o(H6),∀xonΓ,
where WνΓ2 and WνΓ4 are the second and fourth normal derivative of  W on Γ, respectively.Lemma 3.
LetG={y∣y=x+tνΓ,t∈[-H,H],∀xonΓ}. If  ψ∈H1(Ω) and H>0 is small, one has
(37)∥Eψ∥L2(G)≤K∥ψ∥L2(Ω),∥∇(Eψ)·νΓ∥L2(G)≤K∥∇ψ∥L2(Ω),
where
(38)K={1,ifG⊂Ω,2,ifG⊄Ω.Define an elliptic projectionC~~∈ℳhc for c:  ∀v∈ℳj,hc, j=1,2,
(39)(D∇C~~(t),∇v)+(q(t)C~~(t),v)=((D∇c(t),∇v)+(q(t)c(t),v))=(q(t)c̿(t)-ϕ∂c∂t(t)-u(t)·∇c(t),v)+(-1)j+1(g,v)Γ,
where g=(D∇c)·νΓ. Let η=c-C~~. From [25–28], we see the following.Lemma 4.
Forη defined by (39), there holds the following.
L
2-norm error estimate
(40)max0≤t≤T∥η∥L2(Ω)+∥ηt∥L2(0,T;L2(Ω))≤Chcr+1{∥u∥L∞(0,T;Hr+1(Ω))+∥ut∥L2(0,T;Hr+1(Ω))},L∞-norm error estimate
(41)∥η∥L∞(Ω)+∥ηt∥L∞(Ω)≤Chc2|lnhc|{∥u∥W2,∞(Ω)+∥ut∥W2,∞(Ω)},ifr=1,d=2,∥η∥L∞(Ω)+∥ηt∥L∞(Ω)≤Chcr+1{∥u∥Wr+1,∞(Ω)+∥ut∥Wr+1,∞(Ω)},ifr>1.For functionsψ with restrictions in H1(Ω1)∪H1(Ω2), we define a norm
(42)|||ψ|||2=(D∇ψ,∇ψ)+d1KH-1([ψ],[ψ])Γ.Lemma 5.
There exists a positive constantC0=1-(2/2) such that for small H>0,
(43)C0|||ψ|||2≤(D∇ψ,∇ψ)+d1KH-1([ψ],[ψ])Γ+2(ψ¯μ,H,[ψ])Γ,∀ψ∈ℳhc.As we have shown, the combined approximation scheme (28d) includes two terms on the interdomain boundary Γ by integral mean method to present explicit flux calculation. These terms are distinct ones different from Dawson-Dupont's schemes such that the standard elliptic projection (39) is insufficient for optimal error estimates. To get optimal error estimates, we need a new elliptic projection including terms on Γ. This new elliptic projection C~∈ℳhc of the solution c is defined as
(44)(D∇(c-C~),∇v)+((c-C~)¯μ,H,[v])Γ+(v¯μ,H,[c-C~])Γ+d1KH-1([c-C~],[v])Γ=0,∀v∈ℳhc.
It follows from Lemma 5 that the projection problem (44) has a unique solution for small H. Choose the initial C~0 to be the projections of c0(x) by (44) and take the initial conditions
(45)C0=C~0.
Let
(46)ξn=cn-C~n,θn=Cn-C~n.
The following lemma gives the bounds of ξn.Lemma 6.
There holds a priori estimates:(47)∥ξ∥L2(Ω)≤C{hcr+1+H1/2∥η∥L∞(Ω)},(48)∥ξt∥L2(Ω)≤C{hcr+1+H1/2∥ηt∥L∞(Ω)}.


## 3.1.   For the Pressure and Velocity

The Raviart-Thomas spacesVhp×Whp in Section 2.1 possess the following approximation and inverse properties [7]:
(31)(Ap):{infv∈Vhp∥f-v∥≤K3hp∥f∥1,infv∈Vhp∥f-v∥H(div)≤K3hp∥f∥H1(div),infw∈Whp∥g-w∥≤K3hp∥g∥1,(Ip):{∥v∥L∞≤K3hp-1∥v∥,∀v∈Vhp,∥v∥W∞1(Tp)≤K3hp-1∥v∥L∞(Tp),∀v∈Vhp,
where K3 is a positive constant independent of hp and Tp is an element of the mesh 𝒯hp.Define the map(U~,P~):J→Vhp×Whp by (see [7])
(32)A(c;U~,v)+B(v,P~)=(γ(c)∇d,v),∀v∈Vhp,B(U~,w)=-(q,w),∀w∈Whp,
where c is the exact solution of (1a)–(1d).By arguments in [1, 24], the map (U~,P~) exists and (Ap) implies that
(33)∥u-U~∥L∞(H(div))+∥p-P~∥L∞(L2)≤K4{infv∈Vhp∥u-v∥H(div)+infw∈Whp∥p-w∥}≤K4(∥u∥L∞(H1(div)),∥p∥L∞(H1))hp.
The positive constant K4 depends on constant in (Q) but independent of hp, u, p, and c.The estimate (33) and (Ip) imply that
(34)∥U~∥L∞(L∞)≤K5.
As in [1], comparison of (32) and (12) implies that
(35)∥U-U~∥H(div)+∥P-P~∥≤K6(1+∥U~∥L∞)∥c-C∥.
The estimates (33) and (35) will handle the coupling of concentration and velocity errors in the convergence analysis.


## 3.2. For the Concentration

In this section, we adopt some following lemmas [18].Lemma 2.
For smooth enough functionW, there holds estimates:
(36)∥W¯H-W∥L2(Γ)≤2H∥∇W∥L2(Ω),∥W¯H-W∥L∞(Γ)≤CH2∥W∥W2,∞(Ω),W(x)-W¯H(x)=-16H2WνΓ2(x)-1120H4WνΓ4(x)+o(H6),∀xonΓ,
where WνΓ2 and WνΓ4 are the second and fourth normal derivative of  W on Γ, respectively.Lemma 3.
LetG={y∣y=x+tνΓ,t∈[-H,H],∀xonΓ}. If  ψ∈H1(Ω) and H>0 is small, one has
(37)∥Eψ∥L2(G)≤K∥ψ∥L2(Ω),∥∇(Eψ)·νΓ∥L2(G)≤K∥∇ψ∥L2(Ω),
where
(38)K={1,ifG⊂Ω,2,ifG⊄Ω.Define an elliptic projectionC~~∈ℳhc for c:  ∀v∈ℳj,hc, j=1,2,
(39)(D∇C~~(t),∇v)+(q(t)C~~(t),v)=((D∇c(t),∇v)+(q(t)c(t),v))=(q(t)c̿(t)-ϕ∂c∂t(t)-u(t)·∇c(t),v)+(-1)j+1(g,v)Γ,
where g=(D∇c)·νΓ. Let η=c-C~~. From [25–28], we see the following.Lemma 4.
Forη defined by (39), there holds the following.
L
2-norm error estimate
(40)max0≤t≤T∥η∥L2(Ω)+∥ηt∥L2(0,T;L2(Ω))≤Chcr+1{∥u∥L∞(0,T;Hr+1(Ω))+∥ut∥L2(0,T;Hr+1(Ω))},L∞-norm error estimate
(41)∥η∥L∞(Ω)+∥ηt∥L∞(Ω)≤Chc2|lnhc|{∥u∥W2,∞(Ω)+∥ut∥W2,∞(Ω)},ifr=1,d=2,∥η∥L∞(Ω)+∥ηt∥L∞(Ω)≤Chcr+1{∥u∥Wr+1,∞(Ω)+∥ut∥Wr+1,∞(Ω)},ifr>1.For functionsψ with restrictions in H1(Ω1)∪H1(Ω2), we define a norm
(42)|||ψ|||2=(D∇ψ,∇ψ)+d1KH-1([ψ],[ψ])Γ.Lemma 5.
There exists a positive constantC0=1-(2/2) such that for small H>0,
(43)C0|||ψ|||2≤(D∇ψ,∇ψ)+d1KH-1([ψ],[ψ])Γ+2(ψ¯μ,H,[ψ])Γ,∀ψ∈ℳhc.As we have shown, the combined approximation scheme (28d) includes two terms on the interdomain boundary Γ by integral mean method to present explicit flux calculation. These terms are distinct ones different from Dawson-Dupont's schemes such that the standard elliptic projection (39) is insufficient for optimal error estimates. To get optimal error estimates, we need a new elliptic projection including terms on Γ. This new elliptic projection C~∈ℳhc of the solution c is defined as
(44)(D∇(c-C~),∇v)+((c-C~)¯μ,H,[v])Γ+(v¯μ,H,[c-C~])Γ+d1KH-1([c-C~],[v])Γ=0,∀v∈ℳhc.
It follows from Lemma 5 that the projection problem (44) has a unique solution for small H. Choose the initial C~0 to be the projections of c0(x) by (44) and take the initial conditions
(45)C0=C~0.
Let
(46)ξn=cn-C~n,θn=Cn-C~n.
The following lemma gives the bounds of ξn.Lemma 6.
There holds a priori estimates:(47)∥ξ∥L2(Ω)≤C{hcr+1+H1/2∥η∥L∞(Ω)},(48)∥ξt∥L2(Ω)≤C{hcr+1+H1/2∥ηt∥L∞(Ω)}.


## 4. A Priori Error Estimates inL2-Norm

Now, we turn to derive an optimal priori error estimate inL2-norm for the concentration of approximation (28d). Optimal error estimates for velocity in H(div;Ω) and pressure in L2(Ω) follow at once from (33) and (35).It follows from trace theorem in [29] that
(49)∥ψ∥L2(Γ)2≤C2∥ψ∥∥ψ∥1,∀ψ∈H1(Ω)
which will be used in the following proof. Then, we can derive an L2(Ω)-norm error estimate for θn.Lemma 7.
Forθn defined by (46), there holds the following error estimate:
(50)max0≤n≤N∥θn∥2≤C{H[max1≤n≤N(∥ηn∥L∞(Ω)2+∥ηn∥2)+∥ηt∥L2(0,T;L∞(Ω))2+∥ηt∥L2(0,T;L2(Ω))2]+(Δtc)2+H5+hp2l+2+hc2r+2+(Δtp1)3+(Δtp)4{[max1≤n≤N(∥ηn∥L∞(Ω)2+∥ηn∥2)]}},
provided
(51)Δt≤C1H2,C1=d0(1-δ)2C0316d12K2C22,∀0<δ≪1.Proof.
Combining (39) and (44), we have ∀v∈ℳ,
(52)(ϕC~n-C~n-1Δtc,v)+(D∇C~n,∇v)+(C~¯μ,Hn-1,[v])Γ+(v¯μ,H,[C~n-1])Γ+d1KH-1([C~n-1],[v])Γ=(C~¯μ,Hn-1-C~¯μ,Hn,[v])Γ+(v¯μ,H,[C~n-1-C~n])Γ+d1KH-1([C~n-1-C~n],[v])Γ+(qc~n,v)-(qncn,v)-(u·∇c+ϕ∂cn∂t,v)+(ϕC~n-C~n-1Δtc,v)+(c¯μ,Hn-cμn,[v])Γ.
Subtracting (52) from (28d) and taking v=θn, we adopted the skill of [4] to obtain
(53)(ϕθn-θn-1Δtc,θn)+(D∇θn,∇θn)+d1KH-1([θn],[θn])Γ+2(θ¯μ,Hn,[θn])Γ=(ϕ∂cn∂t+Fun·∇cn-ϕcn-cˇn-1Δtc,θn)+((un-Fun)·∇cn,θn)+(ϕξn-ξn-1Δtc,θn)+{(θ¯μ,Hn-θ¯μ,Hn-1,[θn])Γ+(θ¯μ,Hn,[θn-θn-1])Γ+d1KH-1([θn-θn-1],[θn])Γ+(D∇(ξn-ξn-1),∇θn)+(c¯μ,Hn-c¯μ,Hn-1,[θn])Γ+(cμn-c¯μ,Hn,[θn])Γ{(θ¯μ,Hn-θ¯μ,Hn-1,[θn])Γ+(θ¯μ,Hn,[θn-θn-1])Γ}+(ϕc^n-1-cˇn-1Δtc,θn)-(ϕξ^n-1-ξˇn-1Δtc,θn)+(ϕθ^n-1-θˇn-1Δtc,θn)-(ϕξˇn-1-ξn-1Δtc,θn)+(ϕθˇn-1-θn-1Δtc,θn)+(qn(ξn-θn),θn).
Since
(54)(θ¯μ,Hn,[θn-θn-1])Γ=(θ¯μ,Hn,[θn])Γ-(θ¯μ,Hn-1,[θn-1])Γ-(θ¯μ,Hn-θ¯μ,Hn-1,[θn-1])Γ,
by summing (53) over n we have
(55)12(ϕθn,θn)+Δtc[{(θ¯μ,Hn,[θn])Γ}(D∇θn,∇θn)+d1KH-1([θn],[θn])Γ+(θ¯μ,Hn,[θn])Γ]+Δtc∑k=1n-1‍[{2(θ¯μ,Hk,[θk])Γ}(D∇θk,∇θk)+d1KH-1([θk],[θk])Γ+2(θ¯μ,Hk,[θk])Γ]=12(ϕθ0,θ0)-Δtc(θ¯μ,H0,[θ0])Γ+Δtc∑k=1n‍{(θ¯μ,Hk-θ¯μ,Hk-1,[θk-θk-1])Γ+d1KH([θk-θk-1],[θk])Γ+(D∇(ξk-ξk-1),∇θk)+(c¯μ,Hk-c¯μ,Hk-1,[θk])Γ+(cμk-c¯μ,Hk,[θk])Γ{(θ¯μ,Hk-θ¯μ,Hk-1,[θk-θk-1])Γ}}-∑k=1n‍(qk(ξk-θk),θk)Δtc+∑k=1n([ϕ∂cn∂t+Fun·∇cn]-ϕcn-cˇn-1Δtc,θk)Δtc‍+∑k=1n‍((uk-Fuk)·∇ck,θk)Δtc+∑k=1n‍(ϕξk-ξk-1Δtc,θk)Δtc+∑k=1n‍(ϕc^n-1-cˇn-1Δtc,θn)Δtc-∑k=1n‍(ϕξ^k-1-ξˇk-1Δtc,θk)Δtc+∑k=1n‍(ϕθ^k-1-θˇk-1Δtc,θk)Δtc-∑k=1n‍(ϕξˇk-1-ξk-1Δtc,θk)Δtc+∑k=1n‍(ϕθˇk-1-θk-1Δtc,θk)Δtc.
Furthermore, noting that
(56)2(θ¯μ,Hn,[θn])Γ≤ε∥θn∥1,D2+d1K2εH-1∥[θn]∥L2(Γ)2,
where ∥v∥1,D2=(D∇v,∇v), and taking 0<ε=2/2<1, we have
(57)12∥θn∥2+Δtc2(1-C0)((D∇θn,∇θn)+d1KH-1∥[θn]∥L2(Γ)2)+Δtc∑k=1n‍[Δtc2∥∂tθk∥2+C0({∥[θk]∥L2(Γ)2}(D∇θk,∇θk)+d1KH-1∥[θk]∥L2(Γ)2){Δtc2}]≤12∥θ0∥2+Δtc|(θ¯μ,H0,[θ0])Γ|+Δtc∑k=1n‍[(θ¯μ,Hk-θ¯μ,Hk-1,[θk-θk-1])Γ+d1KH([θk-θk-1],[θk])Γ+(D∇(ξk-ξk-1),∇θk)+(u¯μ,Hk-u¯μ,Hk-1,[θk])Γ+(uμk-u¯μ,Hk,[θk])Γ]+[∑k=1n‍|(qkθk,θk)|Δtc+∑k=1n‍|(qkξk,θk)|Δtc]+∑k=1n‍([ϕ∂ck∂t+Fuk·∇ck]-ϕck-cˇk-1Δtc,θk)Δtc+∑k=1n‍((uk-Fuk)·∇ck,θk)Δtc+∑k=1n‍(ϕξk-ξk-1Δtc,θk)Δtc+∑k=1n‍(ϕc^k-1-cˇk-1Δtc,θk)Δtc+∑k=1n‍|(ϕξ^k-1-ξˇk-1Δtc,θk)|Δtc+∑k=1n‍(ϕθ^k-1-θˇk-1Δtc,θk)Δtc+∑k=1n‍(ϕξˇk-1-ξk-1Δtc,θk)Δtc+∑k=1n‍(ϕθˇk-1-θk-1Δtc,θk)Δtc.
We estimate the terms on the right-hand side of (57) one by one. We denote the terms as T1,T2,T3,…,T10 from the third term on the right-hand side of (57). We turn to analyze them one by one.
For the termT1, similarly as that of [18], by taking
(58)Δtc≤d0(1-δ)2C03H216d12K2C22,∀0<δ≪1,
we can obtain
(59)|T1|≤(Δtc)22∑k=1n‍∥∂tθk∥2+C0Δtc∑k=1n‍[∥θk∥1,D2+d0∥θk∥2]+(1-δ)C0Δtc∑k=1n‍∥θk∥1,D2+(1-δ)C0Δtcd1KH∑k=1n‍∥[θk]∥L2(Γ)2+Δtc∑k=1n‍[Δtc2(∥θk∥1,D2+∥∂tξk∥1,D2)+C∥ck-ck-1∥1,D2+δC0d1KH-1∥[θk]∥L2(Γ)2+CH5∥ck∥W2,∞(Ω)2].
The analyses of the terms T2,…,T10 are similar to that of [4]
(60)|T2|≤CΔtc∑k=1n‍(∥θk∥2+∥ξk∥2),|T3|≤CΔtc∑k=1n‍(∥∂2c∂τ2∥L2(tk-1,tk;L2)Δtc+∥θk∥2),|T4|≤CΔtc(Δtp)3∥∂2u∂t2∥L2(t0,tn;L2)+CΔtc∑k=1n‍∥θk∥2,|T5|≤CΔtc∑k=1n‍(hc2r+2∥c∥H1(tk-1,tk;Hr+1)2+∥θk∥2),|T6|+|T7|+|T8|≤CΔtc(hp2l+2+hc2r+2+(Δtc)2(Δtp1)3+(Δtp)4)+ϵΔtc∑k=1n‍∥θk∥12,|T9|+|T10|≤CΔtc(hc2r+2+∑k=1n‍∥θk-1∥2)+ϵΔtc∑k=1n‍∥θk∥12.
Combining the above analyses, we have
(61)12∥θn∥2+Δtc2(1-C0)(∥θn∥1,D2+d1KH-1∥[θn]∥L2(Γ)2)≤12∥θ0∥2+Δtc|(θ¯μ,H0,[θ0])Γ|+CΔtc∑k=1n‍[Δtc(∥θk∥1,D2+∥∂tξk∥1,D2)+H5∥uk∥H2,∞(Ω)2]+C(hp2l+2+hc2r+2+(Δtc)2+(Δtp1)3+(Δtp)4).
Taking θ0=0, applying the Gronwall Lemma and Lemmas 5 and 6, we can derive (50).Applying the triangle inequality, Lemmas 4 and 7, we have the following error theore.Theorem 8.
Suppose that the assumptions(P), (Q), (I), (Ap), (Ip) and r≥1,l≥0 hold. Assume that the discretization parameters obey the relations
(62)Δtc=o(hp),hcr+1=O(hp),(Δtp1)3/2=O(hp),(Δtp)2=O(hp).
Then for hc, hp, and Δtc sufficiently small, the errors of the approximation (28a)–(28d) for (1a)–(1d) satisfy
(63)max0≤n≤N∥cn-Cn∥≤C{Δtc+(Δtp)2+(Δtp1)3/2+hcr+1+hpl+1+H5/2{(Δtp1)3/2}},
provided
(64)Δtc≤C1H2,
where C1 is given by (51). Here C is a positive constant dependent on (P), (Q), (I), (Ap), (Ip), ∥∂2c/∂τ2∥, ∥∂u/∂t∥, and ∥∂2u/∂t2∥, but independent of hc, hp, Δtc, and Δtp.By combining Theorem8 with (33) and (35), we obtain at once the following result.Corollary 9.
Under the assumption of Theorem8, the errors for velocity u, and pressure p are obtained by
(65)max0≤m≤M(∥um-Um∥H(div;Ω)+∥pm-Pm∥)≤C{Δtc+(Δtp)2+(Δtp1)3/2+hcr+1+hpl+1+H5/2}.Remark 10.
As for the accuracy order ofhc and hp, we know that (63) and (65) are optimal for the concentration c, the velocity u and the pressure p.


## 5.   Extensions

To get higher-order accuracy with respect toH, we apply the definition (22) to an extrapolation value C^H(x) on the interdomain boundary Γ. Similar to the analysis of [18], we can present another approximate scheme for (1a)–(1d):
(66)(a)C0=C~0,(b)A(Cm;Um,z)+B(z,Pm)=(γ(Cm)∇d,z),∀z∈Vhp,(c)B(Um,w)=-(qm,w),∀w∈Whp,m≥0,(d)(ϕCn-Cˇn-1Δtc,v)+(D∇Cn,v)+(qnCn,v)+(C^μ,Hn-1,[v])Γ+(v^μ,H,[Cn-1])Γ+d1KH-1([Cn-1],[v])Γ=(qnc~n,v),∀v∈ℳhc,n≥1.Since the differences between two schemes (28a)–(28d) and (66) are the second and third term to calculate the flux on the inner-domain boundary Γ, the convergence analysis of the scheme (66) is similar to that of the scheme (28a)–(28d). For the sake of brevity, we describe the processes of analysis for (66) simply.The proof of Theorem15 is based on the following four lemmas. The former three lemmas are adopted from [18]. We omit the proof of Lemma 14, which is similar to that of Lemma 7 in Section 4.Lemma 11.
For sufficiently smooth functionW, there holds estimates:
(67)∥W^H-W∥L2(Γ)≤22+132H∥∇W∥L2(Ω),∥W^H-W∥L∞(Γ)≤CH4∥W∥W4,∞(Ω),W(x)-W^H(x)=1480H4WνΓ4(x)+o(H6),∀xonΓ,
where WνΓ4 is the fourth normal derivative of  W on Γ.For functionsψ with restrictions in H1(Ω1)∪H1(Ω2), we use the definition of the norm
(68)|||ψ|||2=(D∇ψ,∇ψ)+d1KH-1([ψ],[ψ])Γ.Lemma 12.
There exists a positive constantC~0=1-(22/3) such that for small H>0,
(69)C~0|||ψ|||2≤(D∇ψ,∇ψ)+d1KH-1([ψ],[ψ])Γ+2(ψ^μ,H,[ψ])Γ,∀v∈ℳhc.Similarly as the elliptic projection (44), in order to get optimal error estimates, we introduce an elliptic projection C~∈ℳhc of the solution u as follows:
(70)(D∇(c-C~),∇v)+(c^μ,H-C~^μ,H,[v])Γ+(v^μ,H,[c-C~])Γ+d1KH-1([c-C~],[v])Γ=0,∀v∈ℳhc.
It follows from Lemma 12 that the projection problem (70) has a unique solution for small H.Let(71)ξn=cn-C~n,θn=Cn-C~n.
The following lemma gives the bounds of ξn.Lemma 13.
There holds a priori estimates:(72)∥ξ∥L2(Ω)≤C{hcr+1+H1/2∥η∥L∞(Ω)},∥ξt∥L2(Ω)≤C{hcr+1+H1/2∥ηt∥L∞(Ω)}.Now, we turn to derive anL2(Ω)-norm error estimate for θn.Lemma 14.
Forθn defined by (71), there exists the following error estimate
(73)max0≤n≤N∥θn∥2≤C{H[max1≤n≤N(∥ηn∥L∞(Ω)2+∥ηn∥2)+∥ηt∥L2(0,T;L∞(Ω))2+∥ηt∥L2(0,T;L2(Ω))2{[max1≤n≤N(∥ηn∥L∞(Ω)2+∥ηn∥2)]}]+(Δtc)2+H9+hp2l+2+hc2r+2+(Δtp1)3+(Δtp)4{[max1≤n≤N(∥ηn∥L∞(Ω)2+∥ηn∥2)]}},
provided
(74)Δt≤C~1H2,C~1=a0(1-δ)2C~0316a12K2C22,∀0<δ≪1.Theorem 15.
Suppose that the assumptions(P), (Q), (I), (Ap), (Ip) and r≥1, l≥0 hold. Assume that the discretization parameters obey the relations
(75)Δtc=o(hp),hcr+1=O(hp),(Δtp1)3/2=O(hp),(Δtp)2=O(hp).
Then for hc, hp and Δtc sufficiently small, the errors of the approximation (66) for (1a)–(1d) satisfy
(76)max0≤n≤N∥cn-Cn∥≤C{Δtc+(Δtp)2+(Δtp1)3/2+hcr+1+hpl+1+H9/2},
provided
(77)Δtc≤C~1H2,
where C~1 is given by (74). Here C is a positive constant dependent on (P), (Q), (I), (Ap), (Ip), ∥∂2c/∂τ2∥, ∥∂u/∂t∥, ∥∂2u/∂t2∥, but independent of hc, hp, Δtc and Δtp.By combining Theorem15 with (33) and (35), we obtain at once the following result.Corollary 16.
Under the assumption of Theorem15, the errors for velocity u and pressure p are obtained by
(78)max0≤m≤M(∥um-Um∥H(div;Ω)+∥pm-Pm∥)≤C{Δtc+(Δtp)2+(Δtp1)3/2+hcr+1+hpl+1+H9/2}.Remark 17.
As for the accuracy order ofhc and hp, we know that (76) and (78) are optimal for the concentration c, the velocity u and the pressure p.Remark 18.
From Theorem15, we can know that the scheme (66) has the same accuracy order as that of the scheme (28a)–(28d) with respect to Δtc, hc and hp, except that has an accuracy of higher order for H. This shows that the scheme (66) can use larger width H of middle strip domain than that of the scheme (28a)–(28d) so that the time step constraint is more weaker.


## 6. Numerical Experiments

In this section, we present some numerical experiments for the procedures described above. All computer programs below are written by Fortran 90 code and run on a Lenovo PC with Intel(R) Core i5 3.2 GHz CPU and 4 GB memory. The resulting linear systems of algebraic equations are solved by banded Gaussian elimination. Single precision is used for all calculations.The main purpose of this paper is to analyze the integral mean non-overlapping DDM combined with the characteristic method for the concentration equation. There were many experimental results for the integral mean non-overlapping DDM in [18–20, 22, 23]. For simplicity, we consider the following convection dominated diffusion equation
(79)ϕ∂c∂t+u·∇c-∇·(D∇c)+qc=f,(x,t)∈Ω×(0,T],∂c∂ν=0,(x,t)∈∂Ω×(0,T],c0=0,x∈Ω,t=0,
where, Ω=[0,1]×[0,1], u1=2+x12, u2=1+x22, D=0.05I, ϕ=0.2, q=2. We choose f(x,y,t) suitably so that the exact solution of (79) is  c(x,y,t)=100tx3(1-x)2cos(2πy), see Figure 3.Figure 3
The exact solutionc(x,y,t)=100tx3(1-x)2cos(2πy) at time t=0.5.We consider two scenarios:  (1) the characteristic implicit Galerkin scheme (Charac-teristic-FEM scheme) on uniform mesh; that is, no domain decomposition;   (2) the characteristic integral mean DDM scheme (Characteristic-DDM scheme) on global uniform mesh with two equal sub-domains Ω1=(0,1/2)×(0,1), Ω2=(1/2,1)×(0,1), with the interdomain boundary Γ={1/2}×(0,1). In these runs, we approximate (79) by using linear finite element with 4-node quadrilateral mesh on 20×20, 40×40, 80×80 and 160×160 grids, respectively. For each dynamic domain decomposition case, we still take H5/2=h2 to balance error accuracy with respect to h, H and mesh ration Δt=H5/2 in order to satisfy the conditional stability due to the explicit calculation of flux on the interface.The following Figure4 presents the approximate solution of the characteristic integral mean DDM scheme at time t=0.5 with space step size h=1/40.Figure 4
Approximate solution of the characteristic-DDM scheme att=0.5, h=1/40.Table1 shows L2-norm error estimate of eh=c-C at time t=0.5, where the error order is choosen as α=log(ehl/ehl-1)/log(hl/hl-1).Table 1
L
2-norm error estimate of eh=c-C at time t=0.5.
Grids
h
H
/
h
Characteristic-FEM scheme
Characteristic-DDM scheme
∥
e
h
∥
L
2
α
∥
e
h
∥
L
2
α
50
×
50
.
5000
e
-
01
.
1826
e
+
01
.
0921
e
-
01
.
0878
e
-
01
40
×
40
.
2500
e
-
01
.
2091
e
+
01
.
2448
e
-
02
1.9100
.
2282
e
-
02
1.9439
80
×
80
.
1250
e
-
01
.
2402
e
+
01
.
6245
e
-
03
1.9708
.
5781
e
-
03
1.9809
160
×
160
.
6250
e
-
02
.
2760
e
+
01
.
1583
e
-
03
1.9800
.
1455
e
-
03
1.9903Figure5 shows the convergence order of the characteristic integral mean DDM scheme.Figure 5
Convergence order of the characteristic-DDM scheme.Remark 19.
From Table1, we see that the error estimate of the characteristic integral mean DDM scheme is better than that of the characteristic implicit Galerkin scheme. From Table 1 and Figure 5, we see that the convergence order is near 2 which is consisting with the analytical results.



---

*Source: 303952-2013-02-19.xml*