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

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 in 𝐿²-norm are derived for these two schemes, respectively.

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]. 𝑘(𝑥)−∇⋅𝜇(𝑐)(𝑐)(𝑥)𝜙∇𝑝−𝛾∇𝑑≡∇⋅𝐮=𝑞,𝑥∈Ω,𝑡∈𝐽,𝜕𝑐(𝐷(𝐮))=(̃)(𝐷(𝐮))𝑐()𝜕𝑡+𝐮⋅∇𝑐−∇⋅∇𝑐𝑐−𝑐𝑞,𝑥∈Ω,𝑡∈𝐽,𝐮⋅𝜈=∇𝑐⋅𝜈=0,𝑥∈𝜕Ω,𝑡∈𝐽,𝑥,0=𝑐₀(𝑥),𝑥∈Ω, (1a)(1b)(1c)(1d)

where Ω is a bounded domain in 𝐑^𝟐, 𝐽=(0,𝑇], and 𝑞=max{𝑞,0} is nonzero at injection wells only. The variables in (1a)–(1d) are the pressure 𝑝(𝑥,𝑡) in the fluid mixture, the Darcy velocity 𝐮=(𝑢₁,𝑢₂)^′, and the relative concentration 𝑐(𝑥,𝑡) of the injected fluid. The 𝜈 is the unit outward normal vector on boundary 𝜕Ω.

The coefficients and data in (1a)–(1d) are 𝑘(𝑥): the permeability of the porous media; 𝜇(𝑐), the viscosity of the fluid mixture; 𝑞(𝑥,𝑡): representing flow rates at wells; 𝛾(𝑐) and 𝑑(𝑥), the gravity coefficient and vertical coordinate; 𝜙(𝑥): the porosity of the rock; c̃(𝑥,𝑡), the injected concentration at injection wells (𝑞>0) and the resident concentration at production wells (𝑞<0). Here, 𝐷(𝐮) is a tensor 2×2 matrix and generally has the form 𝐷(𝐮)=𝜙(𝑥)𝑑_𝑚𝐼+|𝐮|𝑑_𝑙𝐸(𝐮)+𝑑_𝑡𝐸^⟂(𝐮), (2) where 2×2 matrix 𝐸=(𝑒_𝑖𝑗) satisfies 𝑒_𝑖𝑗=𝑢_𝑖𝑢_𝑗/|𝐮|², 𝐸^⟂=𝐼−𝐸, 𝑑_𝑚 is the molecular diffusivity, and 𝑑_𝑙, 𝑑_𝑡 are longitudinal and transverse dispersivities, respectively. Furthermore, a compatibility condition ∫_Ω𝑞(𝑥,𝑡)𝑑𝑥=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 ℎ, 𝐿²-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 𝐿²-norm error estimates avoid the loss of 𝐻^−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 𝐿²-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]: 𝐻^𝑚(Ω)=𝜕𝑓∶^|𝛼|𝑓𝜕𝑥^𝛼∈𝐿²(Ω)|𝛼|,𝑓for≤𝑚_2𝑚=_{|𝛼|≤𝑚}𝜕^|𝛼|𝑓𝜕𝑥^𝛼2𝑊,𝑚≥0,_𝑚∞(Ω)=𝜕𝑓∶^|𝛼|𝑓𝜕𝑥^𝛼∈𝐿^∞(Ω)|𝛼|,𝑓for≤𝑚_𝑊𝑚∞=max_{|𝛼|≤𝑚}𝜕^|𝛼|𝑓𝜕𝑥^𝛼_∞,𝑚≥0. (3) In particular, 𝐻⁰(Ω) and 𝑊⁰_∞(Ω)=𝐿^∞(Ω). The inner product on 𝐿²(Ω) is denoted by (⋅,⋅).

We also use the following spaces that incorporate time dependence. Let [𝑎,𝑏]⊂𝐽 and let 𝑋 be any of the spaces just defined. For 𝑓(𝑥,𝑡) suitably smooth on Ω×[𝑎,𝑏], we let 𝐻^𝑚(𝑎,𝑏;𝑋)=𝑓∶^𝑏_𝑎𝜕^𝛼𝑓(⋅,𝑡)𝜕𝑡^𝛼_𝑋,𝑓𝑑𝑡<∞,for𝛼≤𝑚_{𝐻𝑚(𝑎,𝑏;𝑋)}=^𝑚_𝛼=0^𝑏_𝑎𝜕^𝛼𝑓(⋅,𝑡)𝜕𝑡^𝛼_𝑋𝑑𝑡^1/2,𝑚≥0. (4) Similarly, 𝑊^𝑚_∞(𝑎,𝑏;𝑋) and the norm ‖𝑓‖_{𝑊𝑚∞(𝑎,𝑏;𝑋)} are defined. If [𝑎,𝑏]=𝐽, we simplify our notation and write 𝐿^∞(𝑊¹_∞) for 𝐿^∞(0,𝑇;𝑊¹_∞(Ω)), and so forth.

If 𝑓=(𝑓₁,𝑓₂) is a vector function, we note that 𝑓∈𝑋 if 𝑓₁∈𝑋 and 𝑓₂∈𝑋. We also use the vector-function spaces and norms: 𝐻^𝑚(div;Ω)=𝑓∶𝑓₁,𝑓₂,∇⋅𝑓∈𝐻^𝑚(Ω,𝑓_{2𝐻𝑚(div;Ω)})=𝑓₁²_𝑚+𝑓₂²_𝑚+∇⋅𝑓²_𝑚𝐻(),𝑚≥0,div;Ω=𝐻⁰(div;Ω). (5)

We need some assumptions. The regularity assumptions on the solution of (1a)–(1d) are noted by [4]: (𝑃)∶𝑐∈𝐿^∞𝐻²∩𝐻¹𝐻¹∩𝐿^∞𝑊¹_∞∩𝐻²𝐿²,𝑝∈𝐿^∞𝐻¹,𝐮∈𝐿^∞𝐻¹(div)∩𝐿^∞𝑊¹_∞∩𝑊¹_∞𝐿^∞∩𝐻²𝐿². (6)

We also require the following assumptions on the coefficients in (1a)–(1d) [4]. Let 𝑎₀, 𝑎₁, 𝜙₀, 𝜙₁, 𝑑₀, 𝑑₁, 𝐾₁, and 𝐾₂ be positive constants such that (𝑄)∶0<𝑎₀≤𝑘(𝑥)𝜇(𝑐)≤𝑎₁,0<𝜙₀≤𝜙(𝑥)≤𝜙₁,𝜕𝑘/𝜇(+𝜕𝑐𝑥,𝑐)𝜕𝛾(+(𝑥+𝑞(+𝜕𝑞𝜕𝑐𝑥,𝑐)∇𝜙)𝑥,𝑡)(𝜕𝑡𝑥,𝑡)≤𝐾₁,𝜕𝐷(𝐮+𝜕𝐮)𝜕𝐷⁻¹(𝐮𝜕𝐮)≤𝐾₂,𝑑₀𝜉²≤∑²_𝑖,𝑗=1𝐷(𝐮)𝜉_𝑖𝜉_𝑗≤𝑑₁𝜉²,∀𝜉∈ℛ². (7) 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, 𝐶 and 𝐾_𝑖(𝑖=3,…,6) will denote generic positive constants, independent of ℎ_𝑐, ℎ_𝑝, Δ𝑡_𝑐, and Δ𝑡_𝑝, but possibly depending on constants in (𝑄), norms in (𝑃). Similarly, 𝜀 will denote a generic small positive constant.

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

Let (),𝐿𝐻=𝜒∈𝐻div;Ω∣𝜒⋅𝜈=0on𝜕Ω𝑊=²(Ω){}.𝑤≡constantonΩ (8)(9)

As in [4], the pressure equation (1a) is equal to the following saddle-point problem of finding a map (𝐮,𝑝)∶𝐽→𝐻×𝑊 such that (a)𝐴()=𝛾(𝑐)(),(b)𝐵()𝑐;𝐮,𝑣+𝐵𝑣,𝑝∇𝑑,𝑣,∀𝑣∈𝐻div;Ω𝐮,𝑤=−𝑞,𝑤,∀𝑤∈𝐿²(Ω), (10) where the bilinear forms 𝐴()=𝜇(𝑐)𝑐;𝐮,𝑣𝑘(𝑥),𝐵,𝐵()().𝐮,𝑣𝑣,𝑝=−∇⋅𝑣,𝑝𝐮,𝑤=−∇⋅𝐮,𝑤 (11)

For ℎ_𝑝>0, we discrete (10) in space on a quasi-uniform mesh 𝒯_ℎ𝑝 of Ω with diameter of element ≤ℎ_𝑝. Let 𝑉_ℎ𝑝×𝑊_ℎ𝑝⊂𝐻×𝑊 be Raviart-Thomas spaces [6, 7] of index 𝑙=0 for this mesh.

The mixed method for pressure and velocity, given a concentration approximation 𝐶 at a time 𝑡∈𝐽, consists of 𝑈∈𝑉_ℎ𝑝 and 𝑃∈𝑊_ℎ𝑝 such that 𝐴(𝐶;𝑈,𝑣)+𝐵(𝑣,𝑃)=𝛾(𝑐)∇𝑑,𝑣,∀𝑣∈𝑉_ℎ𝑝,𝐵(𝑈,𝑤)=−𝑞,𝑤,∀𝑤∈𝑊_ℎ𝑝. (12) Existence and uniqueness of 𝑈 and 𝑃 are proved in [1], based on ideas of [24].

2.2. A Characteristic DDM for the Concentration

In this section, we assume 𝑢 in (1a)–(1d) is given. Define 𝑀=𝑣∈𝐿²(Ω)∣𝑣ispiecewiseconstant,𝑀isdensein𝐿²(Ω). (13)

Let 𝜓(𝑥)=|𝐮(𝑥)|²+𝜙(𝑥)^21/2=𝑢₁(𝑥)²+𝑢₂(𝑥)²+𝜙(𝑥)^21/2, (14) and let the characteristic direction associated with the operator 𝜙𝜕𝑐/𝜕𝑡+𝐮⋅∇𝑐 be denoted by 𝜏, where 𝜓(𝑥)𝜕𝑐(𝑥)𝜕𝜏=𝜙𝜕𝑐𝜕𝑡+𝐮⋅∇𝑐. (15)

The weak form of (1a)–(1d) is finding a map 𝑐∶𝐽→𝐿²(Ω) such that (a)𝜓𝜕𝑐+(𝐷(𝐮))=(̃)𝜕𝜏,𝜑∇𝑐,∇𝜑𝑐−𝑐𝑞,𝜑,∀𝜑∈𝐿²(Ω),(b)𝑐(𝑥,0)=𝑐₀(𝑥),∀𝑥∈Ω. (16)

Part 𝐽 into 0=𝑡⁰<𝑡¹<⋯<𝑡^𝑁=𝑇, with Δ𝑡^𝑛_𝑐=𝑡^𝑛−𝑡^𝑛−1. The analysis is valid for variable time steps, but we drop the superscript from Δ𝑡_𝑐 for convenience. For functions 𝑓 on Ω×𝐽, we write 𝑓^𝑛(𝑥) for 𝑓(𝑥,𝑡^𝑛). Approximate (𝜕𝑐^𝑛/𝜕𝜏)(𝑥)=(𝜕𝑐/𝜕𝜏)(𝑥,𝑡^𝑛) by a backward difference quotient in the 𝜏-direction, 𝜕𝑐^𝑛(𝑥)≃𝑐𝜕𝜏^𝑛(𝑥)−𝑐^𝑛−1𝑥−𝐮(𝑥)/𝜙(𝑥)Δ𝑡_𝑐Δ𝑡_𝑐1+|𝐮(𝑥)|²/𝜙(𝑥)². (17)

If we let 𝑥=𝑥−(𝐮(𝑥)/𝜙(𝑥))Δ𝑡_𝑐 and 𝑓(𝑥)=𝑓(𝑥), then we get 𝜓𝜕𝑐^𝑛𝑐𝜕𝜏≃𝜙^𝑛−𝑐^𝑛−1Δ𝑡_𝑐. (18) Since the problem is Ω-periodic [4], 𝑐^𝑛−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 𝑥 at time 𝑡^𝑛 to the concentration that would flow to 𝑥 from time 𝑡^𝑛−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 Ω_𝑗 (𝑗=1,2) by an interdomain boundary Γ, which is a surface of dimension 𝑑−1, see Figure 1. We denote by Γ_𝑗=𝜕Ω_𝑗⋂𝜕Ω the part of the boundary of the subdomains which coincides with 𝜕Ω. Denote the unit vector normal to Γ as 𝝂_Γ, which points from Ω₁ toward Ω₂.

Figure 1 

The domain Ω with the interdomain boundary Γ.

Let 𝒯_𝑗,ℎ𝑐 be quasi-uniform partitions of Ω_𝑗 (𝑗=1,2) and 𝒯_ℎ𝑐=𝒯_1,ℎ𝑐⋃𝒯_2,ℎ𝑐. Here, ℎ denotes the maximal element diameter of 𝒯_ℎ𝑐. We construct the finite element space ℳ_ℎ𝑐 on 𝒯_ℎ𝑐 which satisfies the following condition (I) (1)For 𝑗=1,2, let ℳ_𝑗,ℎ𝑐 be a finite element subspace of 𝐻¹(Ω_𝑗), and let ℳ_ℎc⊂𝐿²(Ω) such that if 𝑣∈ℳ_ℎc, then 𝑣∣_Ω𝑗∈ℳ_𝑗,ℎ𝑐. (2)For 𝑗=1,2, 𝑃_𝑟(Ω_𝑗)⊂ℳ_𝑗,ℎ𝑐, where 𝑃_𝑟(Ω_𝑗) is a polynomial space of degree at most 𝑟. (3)For 𝑗=1,2, ℎ∈(0,1], the integer 𝑘≥1, and 𝑢∈𝐻^𝑘(Ω_𝑗), there exists a positive constant 𝐶 independent of ℎ such that inf_{𝑣∈ℳ𝑐𝑗,ℎ}‖𝑢−𝑣‖_𝐻^𝑠_(Ω𝑗)≤𝐶ℎ^𝜎‖𝑢‖_𝐻^𝑘_(Ω𝑗),0≤𝑠≤1, (19)where 𝜎=min(𝑟+1−𝑠,𝑘−𝑠).

From the definition above, we note that functions 𝑣 in ℳ_ℎ𝑐 have a well-defined jump [𝑣] on Γ: [𝑣](𝐱)=𝑣𝐱⁺−𝑣(𝐱⁻),∀𝐱onΓ, (20) where 𝑣(𝐱^±)∶=lim_𝜆→0^±𝑣(𝐱+𝜆𝝂_Γ).

To construct parallel algorithm, for a small constant 0<𝐻<min{diam(Ω₁),diam(Ω₂)}, we introduce an integral mean value of a given function 𝑉∈𝐿²(Ω) on the interdomain boundary Γ as 𝑉_𝐻(𝐱)=12𝐻^𝐻_−𝐻𝑉𝐱+𝜆𝝂_Γ𝑑𝜆,∀𝐱onΓ. (21) Furthermore, we define the extrapolation of 𝑉_𝐻(𝐱) on Γ as 𝑉_𝐻(𝐱)=4𝑉_𝐻/2(𝐱)−𝑉_𝐻(𝐱)3,∀𝐱onΓ. (22)

Generally, near the intersection of boundary 𝜕Ω and inner-boundary Γ, the value of 𝑉 outside Ω is needed for 𝑉_𝐻 and 𝑉_𝐻. If 𝐱∉Ω, let 𝐱∈Ω denote the symmetric point of 𝐱 with respect to 𝜕Ω. For a given function 𝑢∈𝐿²(Ω), we define (𝐱)=𝑢(𝐱)𝑢(𝐱)𝐸𝑢,if𝐱∈Ω,,if𝐱∉Ω. (23) By (23), we know 𝑉_𝐻 and 𝑉_𝐻 have the values on a strip domain 𝐺={𝐲∣𝐲=𝐱+𝜆𝝂_Γ,𝜆∈[−𝐻,𝐻],∀𝐱onΓ}, see Figure 2.

Figure 2 

The strip domain 𝐺 with width 2𝐻.

Now, we present the non-overlapping characteristic DDM for the concentration equation:𝜙𝐶^𝑛−𝐶^𝑛−1Δ𝑡_𝑐,𝑣+𝐃∇𝐶^𝑛,∇𝑣+𝑞^𝑛𝐶^𝑛+,𝑣𝐶^𝑛−1_𝜇,𝐻,[𝑣]_Γ+𝑣_𝜇,𝐻,𝐶^𝑛−1_Γ+𝑑₁𝐾𝐻⁻¹𝐶^𝑛−1,[𝑣]_Γ=𝑞^𝑛𝑐^𝑛̃,𝑣,∀𝑣∈ℳ_ℎ𝑐, (24) where 𝐶^𝑛−1(𝑥)=𝐶^𝑛−1𝑈𝑥−^𝑛(𝑥)𝜙(𝑥)Δ𝑡_𝑐,𝐶_𝜇^𝑛−1=𝐃∇𝐸𝐶^𝑛−1⋅𝝂_Γ. (25) The coefficient 𝐾 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 𝐽 into pressure time steps 0=𝑡₀<𝑡₁<⋯<𝑡_𝑀=𝑇, with Δ𝑡^𝑚_𝑝=𝑡_𝑚−𝑡_𝑚−1. Each pressure step is also a concentration step; that is, for each 𝑚 there exists 𝑛 such that 𝑡_𝑚=𝑡^𝑛; in general, Δ𝑡_𝑝>Δ𝑡_𝑐. We may vary Δ𝑡_𝑝, but except for Δ𝑡¹_𝑝 we drop the superscript. For functions 𝑓 on Ω×𝐽, we write 𝑓_𝑚(𝑥) for 𝑓(𝑥,𝑡_𝑚); thus, subscripts refer to pressure steps and superscripts to concentration steps.

If concentration step 𝑡^𝑛 relates to pressure steps by 𝑡_𝑚−1<𝑡^𝑛≤𝑡_𝑚, we require a velocity approximation for (25) based on 𝑈_𝑚−1 and earlier values. If 𝑚≥2, take the linear extrapolation of 𝑈_𝑚−1 and 𝑈_𝑚−2 defined by 𝐹𝑈^𝑛=𝑡1+^𝑛−𝑡_𝑚−1𝑡_𝑚−1−𝑡_𝑚−2𝑈_𝑚−1−𝑡^𝑛−𝑡_𝑚−1𝑡_𝑚−1−𝑡_𝑚−2𝑈_𝑚−2. (26) If 𝑚=1, set 𝐹𝑈^𝑛=𝑈₀. (27) We retain the superscript on Δ𝑡¹_𝑝 because 𝐹𝑈^𝑛 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 map 𝐶∶{𝑡⁰,𝑡¹,…,𝑡^𝑁}→ℳ_ℎ𝑐 and a map (𝑈,𝑃)∶{𝑡₀,𝑡₁,…,𝑡_𝑀}→𝑉_ℎ𝑝×𝑊_ℎ𝑝 such that𝐶⁰=𝐶⁰,𝐴𝐶_𝑚;𝑈_𝑚,𝑧+𝐵𝑧,𝑃_𝑚=𝛾𝐶_𝑚∇𝑑,𝑧,∀𝑧∈𝑉_ℎ𝑝,𝐵𝑈_𝑚,𝑤=−𝑞_𝑚,𝑤,∀𝑤∈𝑊_ℎ𝑝𝜙𝐶,𝑚≥0,_𝑛−𝐶̌^𝑛−1Δ𝑡_𝑐,𝑣+𝐃∇𝐶^𝑛,𝑣+𝑞^𝑛𝐶^𝑛,𝑣+𝐶^𝑛−1_𝜇,𝐻,[𝑣]_Γ+𝑣_𝜇,𝐻,𝐶^𝑛−1_Γ+𝑑₁𝐾𝐻⁻¹𝐶^𝑛−1,[𝑣]_Γ=𝑞^𝑛𝑐^𝑛̃,𝑣,∀𝑣∈ℳ_ℎ𝑐,𝑛≥1, (28a)(28b)(28c)(28d) where 𝐶̌^𝑛−1(𝑥)=𝐶^𝑛−1(𝑥̌)=𝐶^𝑛−1𝑥−𝐹𝑈^𝑛(𝑥)𝜙(𝑥)Δ𝑡_𝑐, (29)𝐶⁰ is given by (45).

The steps of the above calculation are as follows: Step  1. 𝐶⁰ known → solve (𝑈₀,𝑃₀) by (28b) and (28c); Step  2. On each domain, use (28d) to parallelly compute {𝐶^𝑗}^𝑛1_𝑗=1 until 𝑡^𝑛1=𝑡₁; Step  3. Then by (28b) and (28c) for (𝑈₁,𝑃₁); 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 of 𝑥 defined for the exact velocity 𝐮^𝑛. If 𝑓 is a function on Ω, set 𝑓̌(𝑥)=𝑓(𝑥̌)=𝑓𝑥−𝐹𝐮^𝑛(𝑥)𝜙(𝑥)Δ𝑡_𝑐. (30) The time step 𝑡^𝑛 will be clear from the context.

Remark 1. In the scheme (28a)–(28d), the numerical flux on Γ is computed explicitly from 𝐶^𝑛−1, so that 𝐶^𝑛 can be computed on Ω₁ and Ω₂ fully parallel once 𝐶^𝑛−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 spaces 𝑉_ℎ𝑝×𝑊_ℎ𝑝 in Section 2.1 possess the following approximation and inverse properties [7]: 𝐴_𝑝∶inf_{𝑣∈𝑉ℎ𝑝}𝑓−𝑣≤𝐾₃ℎ_𝑝𝑓₁,inf_{𝑣∈𝑉ℎ𝑝}𝑓−𝑣_𝐻(div)≤𝐾₃ℎ_𝑝𝑓_𝐻1(div),inf_{𝑤∈𝑊ℎ𝑝}𝑔−𝑤≤𝐾₃ℎ_𝑝𝑔₁,𝐼_𝑝∶‖𝑣‖_𝐿^∞≤𝐾₃ℎ_𝑝⁻¹‖𝑣‖,∀𝑣∈𝑉_ℎ𝑝,‖𝑣‖_𝑊¹_{∞(𝑇𝑝)}≤𝐾₃ℎ_𝑝⁻¹‖𝑣‖_𝐿^∞_(𝑇𝑝),∀𝑣∈𝑉_ℎ𝑝, (31) where 𝐾₃ is a positive constant independent of ℎ_𝑝 and 𝑇_𝑝 is an element of the mesh 𝒯_ℎ𝑝.

Define the map (𝑈,𝑃)∶𝐽→𝑉_ℎ𝑝×𝑊_ℎ𝑝 by (see [7]) 𝐴𝑐;𝑈,𝑣+𝐵𝑣,𝑃=𝛾(𝑐)∇𝑑,𝑣,∀𝑣∈𝑉_ℎ𝑝,𝐵𝑈,𝑤=−𝑞,𝑤,∀𝑤∈𝑊_ℎ𝑝, (32) where 𝑐 is the exact solution of (1a)–(1d).

By arguments in [1, 24], the map (𝑈,𝑃) exists and (𝐴_𝑝) implies that 𝑢−𝑈_{𝐿∞(𝐻(div))}+𝑝−𝑃_{𝐿∞(𝐿2)}≤𝐾₄inf_{𝑣∈𝑉ℎ𝑝}‖𝑢−𝑣‖_𝐻(div)+inf_{𝑤∈𝑊ℎ𝑝}𝑝−𝑤≤𝐾₄‖𝑢‖_𝐿^∞_(𝐻¹_(div)),𝑝_{𝐿∞(𝐻1)}ℎ_𝑝. (33) The positive constant 𝐾₄ depends on constant in (𝑄) but independent of ℎ_𝑝, 𝑢, 𝑝, and 𝑐.

The estimate (33) and (𝐼_𝑝) imply that 𝑈_{𝐿∞(𝐿∞)}≤𝐾₅. (34) As in [1], comparison of (32) and (12) implies that 𝑈−𝑈_𝐻(div)+𝑃−𝑃≤𝐾₆1+𝑈_𝐿∞‖𝑐−𝐶‖. (35) 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 function 𝑊, there holds estimates: 𝑊_𝐻−𝑊_𝐿2(Γ)≤√2𝐻‖∇𝑊‖_𝐿²_(Ω),𝑊_𝐻−𝑊_𝐿∞(Γ)≤𝐶𝐻²‖𝑊‖_𝑊^2,∞_(Ω),𝑊(𝐱)−𝑊_𝐻(𝐱)1=−6𝐻²𝑊_𝜈²_Γ(𝐱)−1𝐻120⁴𝑊_𝜈⁴_Γ(𝐱)+𝑜𝐻⁶,∀𝐱𝑜𝑛Γ, (36) where 𝑊_𝜈²_Γ and 𝑊_𝜈⁴_Γ are the second and fourth normal derivative of 𝑊 on Γ, respectively.

Lemma 3. Let 𝐺={𝐲∣𝐲=𝐱+𝑡𝝂_Γ,𝑡∈[−𝐻,𝐻],∀𝐱𝑜𝑛Γ}. If 𝜓∈𝐻¹(Ω) and 𝐻>0 is small, one has 𝐸𝜓_𝐿2(𝐺)≤√𝐾𝜓_𝐿2(Ω),∇𝐸𝜓⋅𝜈_{Γ𝐿2(𝐺)}≤√𝐾∇𝜓_𝐿2(Ω), (37) where 𝐾=1,𝑖𝑓𝐺⊂Ω,2,𝑖𝑓𝐺⊄Ω. (38)

Define an elliptic projection 𝐶∈ℳ_ℎ𝑐 for 𝑐: ∀𝑣∈ℳ_𝑗,ℎ𝑐, 𝑗=1,2, 𝐶(𝑡)+𝑞(𝑡)𝐶(𝑡)=((𝑡))+𝑞(𝑡)𝑐(𝑡)=𝑞(𝑡)𝐃∇,∇𝑣,𝑣𝐃∇𝑐,∇𝑣,𝑣𝑐(𝑡)−𝜙𝜕𝑐(𝑡)(𝑡)(𝑡)+()𝜕𝑡−𝐮⋅∇𝑐,𝑣−1^𝑗+1𝑔,𝑣_Γ, (39) where 𝑔=(𝐃∇𝑐)⋅𝝂_Γ. Let 𝜂=𝑐−𝐶. From [25–28], we see the following.

Lemma 4. For 𝜂 defined by (39), there holds the following.
𝐿²-norm error estimate max_{0≤𝑡≤𝑇}𝜂_𝐿2(Ω)+𝜂_{𝑡𝐿2(0,𝑇;𝐿2(Ω))}≤𝐶ℎ_𝑐^𝑟+1‖𝑢‖_𝐿^∞_(0,𝑇;𝐻^𝑟+1_(Ω))+𝑢_{𝑡𝐿2(0,𝑇;𝐻𝑟+1(Ω))}, (40)𝐿^∞-norm error estimate 𝜂_𝐿∞(Ω)+𝜂_{𝑡𝐿∞(Ω)}≤𝐶ℎ²_clnℎ_c‖𝑢‖_𝑊^2,∞_(Ω)+𝑢_{𝑡𝑊2,∞(Ω)},𝜂𝑖𝑓𝑟=1,𝑑=2,_𝐿∞(Ω)+𝜂_{𝑡𝐿∞(Ω)}≤𝐶ℎ_𝑐^𝑟+1‖𝑢‖_𝑊^𝑟+1,∞_(Ω)+𝑢_{𝑡𝑊𝑟+1,∞(Ω)},𝑖𝑓𝑟>1. (41)

For functions 𝜓 with restrictions in 𝐻¹(Ω₁)∪𝐻¹(Ω₂), we define a norm 𝜓²=𝐷∇𝜓,∇𝜓+𝑑₁𝐾𝐻⁻¹𝜓,𝜓_Γ. (42)

Lemma 5. There exists a positive constant 𝐶₀=1−(√2/2) such that for small 𝐻>0, 𝐶₀𝜓²≤𝐃∇𝜓,∇𝜓+𝑑₁𝐾𝐻⁻¹𝜓,𝜓_Γ+2𝜓_𝜇,𝐻,𝜓_Γ,∀𝜓∈ℳ_ℎ𝑐. (43)

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 𝐶∈ℳ_ℎ𝑐 of the solution 𝑐 is defined as 𝐃∇𝑐−𝐶,∇𝑣+𝑐−𝐶_𝜇,𝐻,[𝑣]_Γ+𝑣_𝜇,𝐻,𝑐−𝐶_Γ+𝑑₁𝐾𝐻⁻¹𝑐−𝐶,[𝑣]_Γ=0,∀𝑣∈ℳ_ℎ𝑐. (44) It follows from Lemma 5 that the projection problem (44) has a unique solution for small 𝐻. Choose the initial 𝐶⁰ to be the projections of 𝑐₀(𝑥) by (44) and take the initial conditions 𝐶⁰=𝐶⁰. (45) Let 𝜉^𝑛=𝑐^𝑛−𝐶^𝑛,𝜃^𝑛=𝐶^𝑛−𝐶^𝑛. (46) The following lemma gives the bounds of 𝜉^𝑛.

Lemma 6. There holds a priori estimates: 𝜉_𝐿2(Ω)≤𝐶ℎ_𝑐^𝑟+1+𝐻^1/2𝜂_𝐿∞(Ω),𝜉_t𝐿2(Ω)≤Cℎ^r+1_c+𝐻^1/2𝜂_{𝑡𝐿∞(Ω)}. (47)(48)

4. A Priori Error Estimates in 𝐿²-Norm

Now, we turn to derive an optimal priori error estimate in 𝐿²-norm for the concentration of approximation (28d). Optimal error estimates for velocity in 𝐻(div;Ω) and pressure in 𝐿²(Ω) follow at once from (33) and (35).

It follows from trace theorem in [29] that 𝜓²_𝐿2(Γ)≤𝐶₂𝜓𝜓₁,∀𝜓∈𝐻¹(Ω) (49) which will be used in the following proof. Then, we can derive an 𝐿²(Ω)-norm error estimate for 𝜃^𝑛.

Lemma 7. For 𝜃^𝑛 defined by (46), there holds the following error estimate: max_{0≤𝑛≤𝑁}𝜃^𝑛2≤𝐶𝐻max_{1≤𝑛≤𝑁}𝜂^𝑛2_𝐿∞(Ω)+𝜂^𝑛2+𝜂_𝑡²_{𝐿2(0,𝑇;𝐿∞(Ω))}+𝜂_𝑡²_{𝐿20,𝑇;𝐿2(Ω)}+Δ𝑡_𝑐²+𝐻⁵+ℎ_𝑝^2𝑙+2+ℎ_𝑐^2𝑟+2+Δ𝑡¹_𝑝³+Δ𝑡_𝑝⁴, (50) provided Δ𝑡≤𝐶₁𝐻²,𝐶₁=𝑑₀(1−𝛿)²𝐶³₀16𝑑²₁𝐾²𝐶²₂,∀0<𝛿≪1. (51)

Proof. Combining (39) and (44), we have ∀𝑣∈ℳ, 𝜙𝐶^𝑛−𝐶^𝑛−1Δ𝑡_𝑐,𝑣+𝐃∇𝐶^𝑛,∇𝑣+𝐶^𝑛−1_𝜇,𝐻,[𝑣]_Γ+𝑣_𝜇,𝐻,𝐶^𝑛−1_Γ+𝑑₁𝐾𝐻⁻¹𝐶^𝑛−1,[𝑣]_Γ=𝐶^𝑛−1_𝜇,𝐻−𝐶^𝑛_𝜇,𝐻,[𝑣]_Γ+𝑣_𝜇,𝐻,𝐶^𝑛−1−𝐶^𝑛_Γ+𝑑₁𝐾𝐻⁻¹𝐶^𝑛−1−𝐶^𝑛,[𝑣]_Γ+𝑞𝑐^𝑛̃,𝑣−𝑞^𝑛𝑐^𝑛,𝑣−𝐮⋅∇𝑐+𝜙𝜕𝑐^𝑛+𝜙𝐶𝜕𝑡,𝑣_𝑛−𝐶^𝑛−1Δ𝑡_𝑐,𝑣+𝑐^𝑛_𝜇,𝐻−𝑐^𝑛_𝜇,[𝑣]_Γ. (52) Subtracting (52) from (28d) and taking 𝑣=𝜃^𝑛, we adopted the skill of [4] to obtain 𝜙𝜃^𝑛−𝜃^𝑛−1Δ𝑡_𝑐,𝜃^𝑛+𝐃∇𝜃^𝑛,∇𝜃^𝑛+𝑑₁𝐾𝐻⁻¹𝜃^𝑛,𝜃^𝑛_Γ+2𝜃^𝑛_𝜇,𝐻,𝜃^𝑛_Γ=𝜙𝜕𝑐^𝑛𝜕𝑡+𝐹𝐮^𝑛⋅∇𝑐^𝑛𝑐−𝜙^𝑛−𝑐̌^𝑛−1Δ𝑡_𝑐,𝜃^𝑛+𝐮^𝑛−𝐹𝐮^𝑛⋅∇𝑐^𝑛,𝜃^𝑛+𝜙𝜉^𝑛−𝜉^𝑛−1Δ𝑡_𝑐,𝜃^𝑛+𝜃^𝑛_𝜇,𝐻−𝜃^𝑛−1_𝜇,𝐻,𝜃^𝑛_Γ+𝜃^𝑛_𝜇,𝐻,𝜃^𝑛−𝜃^𝑛−1_Γ+𝑑₁𝐾𝐻⁻¹𝜃^𝑛−𝜃^𝑛−1,𝜃^𝑛_Γ+𝐃∇𝜉^𝑛−𝜉^𝑛−1,∇𝜃^𝑛+𝑐^𝑛_𝜇,𝐻−𝑐^𝑛−1_𝜇,𝐻,𝜃^𝑛_Γ+𝑐^𝑛_𝜇−𝑐^𝑛_𝜇,𝐻,𝜃^𝑛_Γ+𝜙𝑐̂^𝑛−1−𝑐̌^𝑛−1Δ𝑡_𝑐,𝜃^𝑛−𝜙𝜉̂^𝑛−1−𝜉̌^𝑛−1Δ𝑡_𝑐,𝜃^𝑛+𝜙𝜃^𝑛−1−𝜃̌^𝑛−1Δ𝑡_𝑐,𝜃^𝑛−𝜙𝜉̌^𝑛−1−𝜉^𝑛−1Δ𝑡_𝑐,𝜃^𝑛+𝜙𝜃̌^𝑛−1−𝜃^𝑛−1Δ𝑡_𝑐,𝜃^𝑛+𝑞^𝑛𝜉^𝑛−𝜃^𝑛,𝜃^𝑛. (53) Since 𝜃^𝑛_𝜇,𝐻,𝜃^𝑛−𝜃^𝑛−1_Γ=𝜃^𝑛_𝜇,𝐻,𝜃^𝑛_Γ−𝜃^𝑛−1_𝜇,𝐻,𝜃^𝑛−1_Γ−𝜃^𝑛_𝜇,𝐻−𝜃^𝑛−1_𝜇,𝐻,𝜃^𝑛−1_Γ, (54) by summing (53) over 𝑛 we have 12𝜙𝜃^𝑛,𝜃^𝑛+Δ𝑡_𝑐𝐃∇𝜃^𝑛,∇𝜃^𝑛+𝑑₁𝐾𝐻⁻¹𝜃^𝑛,𝜃^𝑛_Γ+𝜃^𝑛_𝜇,𝐻,𝜃^𝑛_Γ+Δ𝑡_𝑐^𝑛−1_𝑘=1𝐃∇𝜃^𝑘,∇𝜃^𝑘+𝑑₁𝐾𝐻⁻¹𝜃^𝑘,𝜃^𝑘_Γ+2𝜃^𝑘_𝜇,𝐻,𝜃^𝑘_Γ=12𝜙𝜃⁰,𝜃⁰−Δ𝑡_𝑐𝜃⁰_𝜇,𝐻,𝜃⁰_Γ+Δ𝑡_𝑐^𝑛_𝑘=1𝜃^𝑘_𝜇,𝐻−𝜃^𝑘−1_𝜇,𝐻,𝜃^𝑘−𝜃^𝑘−1_Γ+𝑑₁𝐾𝐻𝜃^𝑘−𝜃^𝑘−1,𝜃^𝑘_Γ+𝐃∇𝜉^𝑘−𝜉^𝑘−1,∇𝜃^𝑘+𝑐^𝑘_𝜇,𝐻−𝑐^𝑘−1_𝜇,𝐻,𝜃^𝑘_Γ+𝑐^𝑘_𝜇−𝑐^𝑘_𝜇,𝐻,𝜃^𝑘_Γ−^𝑛_𝑘=1𝑞^𝑘𝜉^𝑘−𝜃^𝑘,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝜕𝑐^𝑛𝜕𝑡+𝐹𝐮^𝑛⋅∇𝑐^𝑛𝑐−𝜙^𝑛−𝑐̌^𝑛−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝐮^𝑘−𝐹𝐮^𝑘⋅∇𝑐^𝑘,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝜉^𝑘−𝜉^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝑐̂^𝑛−1−𝑐̌^𝑛−1Δ𝑡_𝑐,𝜃^𝑛Δ𝑡_𝑐−^𝑛_𝑘=1𝜙𝜉̂^𝑘−1−𝜉̌^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝜃^𝑘−1−𝜃̌^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐−^𝑛_𝑘=1𝜙𝜉̌^𝑘−1−𝜉^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝜃̌^𝑘−1−𝜃^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐. (55) Furthermore, noting that 2𝜃^𝑛_𝜇,𝐻,𝜃^𝑛_Γ≤𝜀𝜃^𝑛2_1,𝐷+𝑑₁𝐾𝐻2𝜀⁻¹𝜃^𝑛2_𝐿2(Γ), (56) where ‖𝑣‖²_1,𝐷=(𝐃∇𝑣,∇𝑣), and taking 0<𝜀=√2/2<1, we have 12𝜃^𝑛2+Δ𝑡_𝑐21−𝐶₀𝐃∇𝜃^𝑛,∇𝜃^𝑛+𝑑₁𝐾𝐻⁻¹𝜃^𝑛2_𝐿2(Γ)+Δ𝑡_𝑐^𝑛_𝑘=1Δ𝑡_𝑐2𝜕_𝑡𝜃^𝑘2+𝐶₀𝐃∇𝜃^𝑘,∇𝜃^𝑘+𝑑₁𝐾𝐻⁻¹𝜃^𝑘2_𝐿2(Γ)≤12𝜃⁰²+Δ𝑡_𝑐𝜃⁰_𝜇,𝐻,𝜃⁰_Γ+Δ𝑡_𝑐^𝑛_𝑘=1𝜃^𝑘_𝜇,𝐻−𝜃^𝑘−1_𝜇,𝐻,𝜃^𝑘−𝜃^𝑘−1_Γ+𝑑₁𝐾𝐻𝜃^𝑘−𝜃^𝑘−1,𝜃^𝑘_Γ+𝐃∇𝜉^𝑘−𝜉^𝑘−1,∇𝜃^𝑘+𝐮^𝑘_𝜇,𝐻−𝐮^𝑘−1_𝜇,𝐻,𝜃^𝑘_Γ+𝐮^𝑘_𝜇−𝐮^𝑘_𝜇,𝐻,𝜃^𝑘_Γ+^𝑛_𝑘=1𝑞^𝑘𝜃^𝑘,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝑞^𝑘𝜉^𝑘,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝜕𝑐^𝑘𝜕𝑡+𝐹𝐮^𝑘⋅∇𝑐^𝑘𝑐−𝜙^𝑘−𝑐̌^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝐮^𝑘−𝐹𝐮^𝑘⋅∇𝑐^𝑘,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝜉^𝑘−𝜉^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝑐̂^𝑘−1−𝑐̌^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝜉̂^𝑘−1−𝜉̌^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝜃^𝑘−1−𝜃̌^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝜉̌^𝑘−1−𝜉^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐+^𝑛_𝑘=1𝜙𝜃̌^𝑘−1−𝜃^𝑘−1Δ𝑡_𝑐,𝜃^𝑘Δ𝑡_𝑐. (57) We estimate the terms on the right-hand side of (57) one by one. We denote the terms as 𝑇₁,𝑇₂,𝑇₃,…,𝑇₁₀ from the third term on the right-hand side of (57). We turn to analyze them one by one.
For the term 𝑇₁, similarly as that of [18], by taking Δ𝑡_𝑐≤𝑑₀(1−𝛿)²𝐶³₀𝐻²16𝑑²₁𝐾²𝐶²₂,∀0<𝛿≪1, (58) we can obtain 𝑇₁≤Δ𝑡_𝑐²2^𝑛_𝑘=1𝜕_𝑡𝜃^𝑘2+𝐶₀Δ𝑡_𝑐^𝑛_𝑘=1𝜃^𝑘2_1,𝐷+𝑑₀𝜃^𝑘2+(1−𝛿)𝐶₀Δ𝑡_𝑐^𝑛_𝑘=1𝜃^𝑘2_1,𝐷+(1−𝛿)𝐶₀Δ𝑡_𝑐𝑑₁𝐾𝐻^𝑛_𝑘=1𝜃^𝑘2_𝐿2(Γ)+Δ𝑡_𝑐^𝑛_𝑘=1Δ𝑡_𝑐2𝜃^𝑘2_1,𝐷+𝜕_𝑡𝜉^𝑘2_1,𝐷+𝐶𝑐^𝑘−𝑐^𝑘−12_1,𝐷+𝛿𝐶₀𝑑₁𝐾𝐻⁻¹𝜃^𝑘2_𝐿2(Γ)+𝐶𝐻⁵𝑐^𝑘2_{𝑊2,∞(Ω)}. (59) The analyses of the terms 𝑇₂,…,𝑇₁₀ are similar to that of [4] 𝑇₂≤𝐶Δ𝑡_𝑐^𝑛_𝑘=1𝜃^𝑘2+𝜉^𝑘2,𝑇₃≤𝐶Δ𝑡_𝑐^𝑛_𝑘=1𝜕²𝑐𝜕𝜏²_{𝐿2(𝑡𝑘−1,𝑡𝑘;𝐿2)}Δ𝑡_𝑐+𝜃^𝑘2,𝑇₄≤𝐶Δ𝑡_𝑐Δ𝑡_𝑝³𝜕²𝐮𝜕𝑡²_{𝐿2(𝑡0,𝑡𝑛;𝐿2)}+𝐶Δ𝑡_𝑐^𝑛_𝑘=1𝜃^𝑘2,𝑇₅≤𝐶Δ𝑡_𝑐^𝑛_𝑘=1ℎ_𝑐^2𝑟+2‖𝑐‖²_{𝐻1𝑡𝑘−1,𝑡𝑘;𝐻𝑟+1}+𝜃^𝑘2,𝑇₆+𝑇₇+𝑇₈≤𝐶Δ𝑡_𝑐ℎ_𝑝^2𝑙+2+ℎ_𝑐^2𝑟+2+Δ𝑡_𝑐²Δ𝑡¹_𝑝³+Δ𝑡_𝑝⁴+𝜖Δ𝑡_𝑐^𝑛_𝑘=1𝜃^𝑘2₁,𝑇₉+𝑇₁₀≤𝐶Δ𝑡_𝑐ℎ_𝑐^2𝑟+2+^𝑛_𝑘=1𝜃^𝑘−12+𝜖Δ𝑡_𝑐^𝑛_𝑘=1𝜃^𝑘2₁. (60) Combining the above analyses, we have 12𝜃^𝑛2+Δ𝑡_𝑐21−𝐶₀𝜃^𝑛2_1,𝐷+𝑑₁𝐾𝐻⁻¹𝜃^𝑛2_𝐿2(Γ)≤12𝜃⁰²+Δ𝑡_𝑐𝜃⁰_𝜇,𝐻,𝜃⁰_Γ+𝐶Δ𝑡_𝑐^𝑛_𝑘=1Δ𝑡_𝑐𝜃^𝑘2_1,𝐷+𝜕_𝑡𝜉^𝑘2_1,𝐷+𝐻⁵𝐮^𝑘2_{𝐻2,∞(Ω)}+𝐶ℎ_𝑝^2𝑙+2+ℎ_𝑐^2𝑟+2+Δ𝑡_𝑐²+Δ𝑡¹_𝑝³+Δ𝑡_𝑝⁴. (61) Taking 𝜃⁰=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 (𝑃), (𝑄), (𝐼), (𝐴_𝑝), (𝐼_𝑝) and 𝑟≥1,𝑙≥0 hold. Assume that the discretization parameters obey the relations Δ𝑡_𝑐=𝑜ℎ_𝑝,ℎ_𝑐^𝑟+1=𝑂ℎ_𝑝,Δ𝑡¹_𝑝^3/2=𝑂ℎ_𝑝,Δ𝑡_𝑝²=𝑂ℎ_𝑝. (62) Then for ℎ_𝑐, ℎ_𝑝, and Δ𝑡_𝑐 sufficiently small, the errors of the approximation (28a)–(28d) for (1a)–(1d) satisfy max_{0≤𝑛≤𝑁}𝑐^𝑛−𝐶^𝑛≤𝐶Δ𝑡_𝑐+Δ𝑡_𝑝²+Δ𝑡¹_𝑝^3/2+ℎ_𝑐^𝑟+1+ℎ_𝑝^𝑙+1+𝐻^5/2, (63) provided Δ𝑡_𝑐≤𝐶₁𝐻², (64) where 𝐶₁ is given by (51). Here 𝐶 is a positive constant dependent on (𝑃), (𝑄), (𝐼), (𝐴_𝑝), (𝐼_𝑝), ‖𝜕²𝑐/𝜕𝜏²‖, ‖𝜕𝐮/𝜕𝑡‖, and ‖𝜕²𝐮/𝜕𝑡²‖, but independent of ℎ_𝑐, ℎ_𝑝, Δ𝑡_𝑐, and Δ𝑡_𝑝.

By combining Theorem 8 with (33) and (35), we obtain at once the following result.

Corollary 9. Under the assumption of Theorem 8, the errors for velocity 𝐮, and pressure 𝑝 are obtained by max_{0≤𝑚≤𝑀}𝐮_𝑚−𝑈_{𝑚𝐻(div;Ω)}+𝑝_𝑚−𝑃_𝑚≤𝐶Δ𝑡_𝑐+Δ𝑡_𝑝²+Δ𝑡¹_𝑝^3/2+ℎ_𝑐^𝑟+1+ℎ_𝑝^𝑙+1+𝐻^5/2. (65)

Remark 10. As for the accuracy order of ℎ_𝑐 and ℎ_𝑝, we know that (63) and (65) are optimal for the concentration 𝑐, the velocity 𝐮 and the pressure 𝑝.

5. Extensions

To get higher-order accuracy with respect to 𝐻, we apply the definition (22) to an extrapolation value 𝐶_𝐻(𝐱) on the interdomain boundary Γ. Similar to the analysis of [18], we can present another approximate scheme for (1a)–(1d): (a)𝐶⁰=𝐶⁰,(b)𝐴𝐶_𝑚;𝑈_𝑚,𝑧+𝐵𝑧,𝑃_𝑚=𝛾𝐶_𝑚∇𝑑,𝑧,∀𝑧∈𝑉_ℎ𝑝,(c)𝐵𝑈_𝑚,𝑤=−𝑞_𝑚,𝑤,∀𝑤∈𝑊_ℎ𝑝(d𝜙𝐶,𝑚≥0,)^𝑛−𝐶̌^𝑛−1Δ𝑡_𝑐,𝑣+𝐃∇𝐶^𝑛,𝑣+𝑞^𝑛𝐶^𝑛+𝐶,𝑣_{𝑛−1𝜇,𝐻},[𝑣]_Γ+𝑣_𝜇,𝐻,𝐶^𝑛−1_Γ+𝑑₁𝐾𝐻⁻¹𝐶^𝑛−1,[𝑣]_Γ=𝑞^𝑛𝑐^𝑛̃,𝑣,∀𝑣∈ℳ_ℎ𝑐,𝑛≥1. (66)

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 Theorem 15 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 function 𝑊, there holds estimates: 𝑊_𝐻−𝑊_𝐿2(Γ)≤2√2+13√2𝐻‖∇𝑊‖_𝐿²_(Ω),𝑊_𝐻−𝑊_𝐿∞(Γ)≤𝐶𝐻⁴‖𝑊‖_𝑊^4,∞_(Ω),𝑊(𝐱)−𝑊_𝐻(𝐱)=1𝐻480⁴𝑊_𝜈⁴_Γ(𝐱)+𝑜𝐻⁶,∀𝐱𝑜𝑛Γ, (67) where 𝑊_𝜈⁴_Γ is the fourth normal derivative of 𝑊 on Γ.

For functions 𝜓 with restrictions in 𝐻¹(Ω₁)∪𝐻¹(Ω₂), we use the definition of the norm 𝜓²=𝐷∇𝜓,∇𝜓+𝑑₁𝐾𝐻⁻¹𝜓,𝜓_Γ. (68)

Lemma 12. There exists a positive constant 𝐶₀=1−(2√2/3) such that for small 𝐻>0, 𝐶₀𝜓²≤𝐷∇𝜓,∇𝜓+𝑑₁𝐾𝐻⁻¹𝜓,𝜓_Γ+2𝜓_𝜇,𝐻,𝜓_Γ,∀𝑣∈ℳ_ℎ𝑐. (69)

Similarly as the elliptic projection (44), in order to get optimal error estimates, we introduce an elliptic projection 𝐶∈ℳ_ℎ𝑐 of the solution 𝑢 as follows: 𝐷∇𝑐−𝐶,∇𝑣+𝑐̂_𝜇,𝐻−𝐶_𝜇,𝐻,[𝑣]_Γ+𝑣_𝜇,𝐻,𝑐−𝐶_Γ+𝑑₁𝐾𝐻⁻¹𝑐−𝐶,[𝑣]_Γ=0,∀𝑣∈ℳ_ℎ𝑐. (70) It follows from Lemma 12 that the projection problem (70) has a unique solution for small 𝐻.

Let 𝜉^𝑛=𝑐^𝑛−𝐶^𝑛,𝜃^𝑛=𝐶^𝑛−𝐶^𝑛. (71) The following lemma gives the bounds of 𝜉^𝑛.

Lemma 13. There holds a priori estimates: 𝜉_𝐿2(Ω)≤𝐶ℎ_𝑐^𝑟+1+𝐻^1/2𝜂_𝐿∞(Ω),𝜉_{𝑡𝐿2(Ω)}≤𝐶ℎ_𝑐^𝑟+1+𝐻^1/2𝜂_{𝑡𝐿∞(Ω)}. (72)

Now, we turn to derive an 𝐿²(Ω)-norm error estimate for 𝜃^𝑛.

Lemma 14. For 𝜃^𝑛 defined by (71), there exists the following error estimate max_{0≤𝑛≤𝑁}𝜃ⁿ²≤𝐶𝐻max_{1≤𝑛≤𝑁}𝜂^𝑛2_𝐿∞(Ω)+𝜂^𝑛2+𝜂_𝑡²_{𝐿2(0,𝑇;𝐿∞(Ω))}+𝜂_𝑡²_{𝐿20,𝑇;𝐿2(Ω)}+Δ𝑡_𝑐²+𝐻⁹+ℎ_𝑝^2𝑙+2+ℎ_𝑐^2𝑟+2+Δ𝑡¹_𝑝³+Δ𝑡_𝑝⁴, (73) provided Δ𝑡≤𝐶₁𝐻²,𝐶₁=𝑎₀(1−𝛿)²𝐶³₀16𝑎²₁𝐾²𝐶²₂,∀0<𝛿≪1. (74)

Theorem 15. Suppose that the assumptions (𝑃), (𝑄), (𝐼), (𝐴_𝑝), (𝐼_𝑝) and 𝑟≥1, 𝑙≥0 hold. Assume that the discretization parameters obey the relations Δ𝑡_𝑐=𝑜ℎ_𝑝,ℎ_𝑐^𝑟+1=𝑂ℎ_𝑝,Δ𝑡¹_𝑝^3/2=𝑂ℎ_𝑝,Δ𝑡_𝑝²=𝑂ℎ_𝑝. (75) Then for ℎ_𝑐, ℎ_𝑝 and Δ𝑡_𝑐 sufficiently small, the errors of the approximation (66) for (1a)–(1d) satisfy max_{0≤𝑛≤𝑁}𝑐^𝑛−𝐶^𝑛≤𝐶Δ𝑡_𝑐+Δ𝑡_𝑝²+Δ𝑡¹_𝑝^3/2+ℎ_𝑐^𝑟+1+ℎ^l_𝑝⁺¹+𝐻^9/2, (76) provided Δ𝑡_𝑐≤𝐶₁𝐻², (77) where 𝐶₁ is given by (74). Here 𝐶 is a positive constant dependent on (𝑃), (𝑄), (𝐼), (𝐴_𝑝), (𝐼_𝑝), ‖𝜕²𝑐/𝜕𝜏²‖, ‖𝜕𝐮/𝜕𝑡‖, ‖𝜕²𝐮/𝜕𝑡²‖, but independent of ℎ_𝑐, ℎ_𝑝, Δ𝑡_𝑐 and Δ𝑡_𝑝.

By combining Theorem 15 with (33) and (35), we obtain at once the following result.

Corollary 16. Under the assumption of Theorem 15, the errors for velocity 𝐮 and pressure 𝑝 are obtained by max_{0≤𝑚≤𝑀}𝐮_𝑚−𝑈_{𝑚𝐻(div;Ω)}+𝑝_𝑚−𝑃_𝑚≤𝐶Δ𝑡_𝑐+Δ𝑡_𝑝²+Δ𝑡¹_𝑝^3/2+ℎ_𝑐^𝑟+1+ℎ_𝑝^𝑙+1+𝐻^9/2. (78)

Remark 17. As for the accuracy order of ℎ_c and ℎ_𝑝, we know that (76) and (78) are optimal for the concentration 𝑐, the velocity 𝐮 and the pressure 𝑝.

Remark 18. From Theorem 15, we can know that the scheme (66) has the same accuracy order as that of the scheme (28a)–(28d) with respect to Δ𝑡_𝑐, ℎ_𝑐 and ℎ_𝑝, except that has an accuracy of higher order for 𝐻. This shows that the scheme (66) can use larger width 𝐻 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 𝜙𝜕𝑐()()(],𝜕𝑡+𝐮⋅∇𝑐−∇⋅𝐃∇𝑐+𝑞𝑐=𝑓,𝑥,𝑡∈Ω×0,T𝜕𝑐()(],𝑐𝜕𝝂=0,𝑥,𝑡∈𝜕Ω×0,T₀=0,𝑥∈Ω,𝑡=0, (79) where, Ω=[0,1]×[0,1], 𝑢₁=2+𝑥²₁, 𝑢₂=1+𝑥²₂, 𝐃=0.05𝐈, 𝜙=0.2, 𝑞=2. We choose 𝑓(𝑥,𝑦,𝑡) suitably so that the exact solution of (79) is 𝑐(𝑥,𝑦,𝑡)=100𝑡𝑥³(1−𝑥)²cos(2𝜋𝑦), see Figure 3.

Figure 3 

The exact solution 𝑐(𝑥,𝑦,𝑡)=100𝑡𝑥3(1−𝑥)2cos(2𝜋𝑦) at time 𝑡=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 Ω₁=(0,1/2)×(0,1), Ω₂=(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 𝐻^5/2=ℎ² to balance error accuracy with respect to ℎ, 𝐻 and mesh ration Δ𝑡=𝐻^5/2 in order to satisfy the conditional stability due to the explicit calculation of flux on the interface.

The following Figure 4 presents the approximate solution of the characteristic integral mean DDM scheme at time 𝑡=0.5 with space step size ℎ=1/40.

Figure 4 

Approximate solution of the characteristic-DDM scheme at 𝑡=0.5, ℎ=1/40.

Table 1 shows 𝐿²-norm error estimate of 𝑒_ℎ=𝑐−𝐶 at time 𝑡=0.5, where the error order is choosen as 𝛼=log(𝑒_ℎ𝑙/𝑒_ℎ𝑙−1)/log(ℎ_𝑙/ℎ_𝑙−1).

Table 1 

𝐿²-norm error estimate of 𝑒_ℎ=𝑐−𝐶 at time 𝑡=0.5.

Figure 5 shows the convergence order of the characteristic integral mean DDM scheme.

Figure 5 

Convergence order of the characteristic-DDM scheme.

Remark 19. From Table 1, 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.

Acknowledgments

The authors thank the referees for their valuable suggestions and constructive comments which led to great improvement of this paper. This paper was supported by the NSF of China (no. 11271231), Tian Yuan Special Foundation (no. 11126206), the Scientific Research Award Fund for Excellent Middle-Aged and Young Scientists of Shandong Province (no. BS2009NJ003) and the NSF of Shandong Province (no. ZR2010AQ019), China.