J Sci Comput https://doi.org/10.1007/s10915-018-0752-4

Local Discontinuous Galerkin Method with Implicit–Explicit Time Marching for Incompressible Miscible Displacement Problem in Porous Media Haijin Wang1 · Jingjing Zheng2 · Fan Yu2 · Hui Guo2 · Qiang Zhang3

Received: 24 December 2017 / Revised: 9 May 2018 / Accepted: 29 May 2018 © Springer Science+Business Media, LLC, part of Springer Nature 2018

Abstract In this paper, we shall present two fully-discrete local discontinuous Galerkin methods, coupled with multi-step implicit–explicit time discretization up to second order, for solving the two-dimensional incompressible miscible displacement problem. To avoid the solving of nonlinear algebraic systems, the extrapolation linearization is adopted to diffusion– dispersion tensor. Under weak temporal-spatial conditions, the optimal error estimates in L ∞ (L 2 ) norm for both concentration and velocity are derived. Numerical experiments are also given to demonstrate the theoretical results.

H. Wang: Research supported by NSFC Grants 11601241 and 11671199, Natural Science Foundation of Jiangsu Province Grant BK20160877, and NUPTSF Grant NY215067. J. Zheng, F. Yu and H. Guo: Research supported by NSFC Grant 11571367 and the Fundamental Research Funds for the Central Universities. Q. Zhang: Research supported by NSFC Grants 11671199 and 11571290.

B

Qiang Zhang [email protected] Haijin Wang [email protected] Jingjing Zheng [email protected] Fan Yu [email protected] Hui Guo [email protected]

1

School of Science, Nanjing University of Posts and Telecommunications, Nanjing 210023, Jiangsu Province, People’s Republic of China

2

College of Science, China University of Petroleum, Qingdao 266580, Shandong Province, People’s Republic of China

3

Department of Mathematics, Nanjing University, Nanjing 210093, Jiangsu Province, People’s Republic of China

123

J Sci Comput

Keywords Local discontinuous Galerkin method · Implicit–explicit time-marching scheme · Incompressible miscible displacement problem · Extrapolation linearization · Error estimates Mathematics Subject Classification 65M12 · 65M15 · 65M60

1 Introduction The incompressible miscible displacement in porous media have wide applications in various engineering areas, such as reservoir simulations and exploration of underground water and oil. The classical equations in two-dimensional space are given as follows: ∇ · u = q, −κ(x, y) u= ∇ p, μ(c) ˜ φ(x, y)ct + ∇ · (uc) = ∇ · (D(u)∇c) + cq,

(1.1a) (1.1b) (1.1c)

which is a coupled system of nonlinear partial differential equations. Here ct is the firstorder time derivative of c. In the above system, (1.1a) together with (1.1b) are called flow equation, and (1.1c) is called transport equation. The unknown functions p, u and c are the pressure, the velocity, and the concentration of interested species, respectively. The given functions φ and κ are porosity and permeability of the medium, respectively. Moreover, μ is the concentration-dependent viscosity, and q is the external volumetric flow rate. The concentration of fluid in the external flow c˜ depends on q, specifically,  c , if q > 0, c˜ = ∗ (1.2) c, if q < 0, where c∗ is a given injection concentration. The second-order matrix-valued function D(u) is the diffusion–dispersion tensor which may be given in different forms (see [4] for details). For simplicity, we assume in this paper that the system is defined in a bounded rectangle domain  subject to the initial condition c(x, y, 0) = c0 (x, y), (x, y) ∈ ,



(1.3)

and periodic boundary conditions. Moreover, a compatibility condition  pdxdy = 0 must be imposed to uniquely determine the pressure. Numerical simulations and analysis for system (1.1) were investigated extensively in the last several decades, such as the standard Galerkin–Galerkin methods [13,19], the Galerkinmixed finite element method [10,18], characteristic-mixed methods [11,12], upstream weighting and mixed finite elements [17], Eulerian–Lagrangian localized adjoint methods [23], discontinuous Galerkin (DG) methods [1,20,27], mixed finite element method with DG [22], and so on. Recently, the local discontinuous Galerkin (LDG) methods [15,16] is considered also. The LDG method was introduced by Cockburn and Shu in [8] for solving nonlinear convection diffusion equations, motivated by the work of Bassi and Rebay [2] for the compressible Navier–Stokes equations. The idea of LDG method is to rewrite high order equations into an equivalent first order system, then apply the DG method to the system. The LDG method shares the advantages of DG methods, it can easily handle meshes with hanging nodes, elements of general shapes and local spaces of different types, thus it is flexible for hp-adaptivity.

123

J Sci Comput

Besides, a key advantage of this method is the local solvability, that is, the auxiliary variables approximating the derivatives of the solution can be locally eliminated [5,8]. The methods were further developed in [28–31] for solving many nonlinear wave equations with higher order derivatives. In [15], the LDG method for solving (1.1) was just presented and analyzed in the semidiscrete framework. Actually, the efficient time discretization is also a very important aspect. As is well-known, explicit time-marching often suffer from stringent time-step restriction. Hence, we would like in this paper to consider the implicit–explicit (IMEX) time-marching. Roughly speaking, the main advantage of IMEX methods is that larger time step can be taken, compared with the explicit time-marchings. The IMEX time-marching, coupled with the LDG methods for convection–diffusion problems, has been studied in [24–26], in which the diffusion part was given in a linear form, and the unconditional stability and optimal error estimates in L ∞ (L 2 ) norm were presented. Actually, even for convection–diffusion problems with nonlinear diffusion terms, the IMEX time-marching would also show their advantages in obtaining an elliptic algebraic system, which is easy to solve by many iterative methods. If both convection term and diffusion term are treated implicitly, the resulting algebraic system will be far from elliptic, and convergence of many iterative solvers will suffer. In this paper, we will adopt the multi-step IMEX time-marching [14] and present two fully-discrete LDG methods to solve (1.1). However, a direct application of those algorithms in [14] will lead to a nonlinear algebraic system, in which the velocity and concentration are coupled, since the diffusion–dispersion tensor in transport equation depends on the velocity in flow equation. To improve the efficiency in each time-marching, we would like to decouple the numerical computation of flow equation and transport equation and avoid the solving of nonlinear algebraic system. To this purpose, we adopt the extrapolation technique to linearize the diffusion–dispersion tensor in transport equation. Similar idea can be found in [10,13], where a linearized semi-implicit Euler scheme was adopted. Based on the above discussion, we will design two schemes in this paper, where the explicit time discretization is adopted for the convection part of transport equation, meanwhile the implicit time discretization as well as the linearization techniques are adopted for the nonlinear diffusion part of transport equation. By virtue of energy analysis, the optimal error estimate in L ∞ (L 2 ) norm are obtained, for both the concentration and the velocity, under weak temporal-spatial conditions. Namely, the order of accuracy can achieve optimal in both space and time; seeing Theorem 3.1. However, the linearization techniques will cause some difficulties in the theoretical analysis, especially when high-order time-marching algorithm is used. To overcome this issue, we would like to introduce a series of suitable reference functions at each time level. The paper is organized as follows. In Sect. 2, we present some preliminaries, including hypotheses, the basic notations, finite element space, the used projections and so on. The next three sections are the main body of the paper. In Sect. 3, we present two IMEX time-marching of LDG schemes and the main theorem with respect to optimal error estimates in L ∞ (L 2 ) norm. In Sects. 4 and 5, we prove the kernel lemma for the first-order and second-order time-marching scheme, respectively. In Sect. 6, numerical results are given to demonstrate the accuracy of the two proposed methods. Finally, some concluding remarks are presented in Sect. 7.

123

J Sci Comput

2 Preliminaries In this section we give some assumptions on the incompressible miscible displacement problem, the discontinuous finite element space, basic notations and some useful projections.

2.1 Hypotheses Following [15], we assume that the exact solution of (1.1) exists uniquely and is smooth enough. Furthermore, we also make the following hypotheses (H): 1. There exist five positive constants κ∗ , κ ∗ , φ∗ , φ ∗ and q ∗ , such that κ∗ ≤ κ(x, y) ≤ κ ∗ , φ∗ ≤ φ(x, y) ≤ φ ∗ , |q(x, y)| ≤ q ∗ , for any (x, y) ∈ . 2. The functions c, ∇c and u are bounded for any t ∈ [0, T ] and (x, y) ∈ , where T > 0 is the final time. 3. There exist two positive constants μ∗ and μ∗ , such that μ∗ ≤ μ(c) ≤ μ∗ holds for any considered c. Moreover, μ(c) is Lipschitz continuous with respect to c. 4. The diffusion–dispersion tensor D(u) is smooth enough, such that each element and its derivatives up to the second order are bounded uniformly [3] for any considered u. Moreover, D(u) is symmetric and positive definite uniformly, namely, there exists a positive constant D∗ such that v  D(u)v ≥ D∗ v 2 , v ∈ R2 , holds for any considered u. Let A(c) := A(x, y, c) = μ(c)/κ(x, y). From hypotheses 1 and 3, it is easy to see that μ∗ /κ ∗ ≤ A(c) ≤ μ∗ /κ∗ ,

(2.1)

and A(c) is also Lipschitz continuous with respect to c.

2.2 Notations Some standard notations of Sobolev space will be used. For any domain D and non-negative integer s, let H s (D) = W s,2 (D) and W s,∞ (D) respectively represent the space equipped with the norm · s,D and · s,∞,D , in which the function and its derivatives up to the s-th order are square-integrable and essentially bounded, respectively. Similarly, we define the norms 1/2  u s,D = u 1 2s,D + u 2 2s,D , u s,∞,D = max u i s,∞,D , 1≤i≤2

D s,D =

2  i, j=1

di j 2s,D

1/2

(2.2) ,

D s,∞,D = max di j s,∞,D , 1≤i, j≤2

for vector-valued function u = (u 1 , u 2 ) and matrix-valued function D = {di j }2×2 , respectively. If s = 0 or D = , the related subscripts are omitted. We also use L 2 (H s ), L ∞ (H s ), L ∞ (L ∞ ) and L ∞ (W 1,∞ ) to denote the usual space-time Sobolev spaces. Here the time interval [0, T ] and the space domain  are omitted.

123

J Sci Comput

2.3 Discontinuous Finite Element Space j=1,2,...,N

Let h = {K i j }i=1,2,...,Nxy be a partition of  with rectangular element K i j = Ii × J j , where Ii = (xi−1/2 , xi+1/2 ) and J j = (y j−1/2 , y j+1/2 ). Here h ix = xi+1/2 − xi−1/2 is the y y element length, h j = y j+1/2 − y j−1/2 is the element width, and h = maxi, j {h ix , h j } is the mesh parameter. We also assume that h is quasi-uniform in this paper, namely, the quantity y maxi {h/ h ix } and max j {h/ h j } are upper bounded by a positive constant, as h goes to zero. Associated with this partition, we define the discontinuous finite element space Whk = {v : v| K ∈ Q k (K ), ∀K ∈ h }, k ≥ 0,

(2.3)

which is contained in the broken Sobolev space H 1 (h ) = {v : v ∈ L 2 (), v| K ∈ H 1 (K ), ∀K ∈ h }.

(2.4)

Here Q k (K ) denotes the space of polynomials of degree at most k in each variable on K . The functions in H 1 (h ) are permitted to be discontinuous across element interfaces. On each vertical edge, v + and v − represent the values taken from right and left, respectively. On each horizontal edge, v + and v − represent the values taken from above and below, respectively. The corresponding jump is defined by [[v]] = v + − v − . Let w and v be scalar functions in H 1 (h ). The L 2 inner products on each element and its boundaries are defined as 



(w, v) K = K

wvdxdy, w, v ∂ K =

∂K

wvds.

The sum over all elements are denoted by (w, v) =



(w, v) K , w, v =

K ∈h



w, v ∂ K .

K ∈h

The associated norm with ·, · is defined by w h =



w − 2e + w + 2e

1/2

,

(2.5)

e∈h

where h is the set of all element interfaces, and · e means the L 2 norm on edge e. The classical inverse properties from [6] will be used in this paper. For any function w ∈ Whk , there holds h w ∞ + h 1/2 w h ≤ C w , (2.6) where the inverse constant C > 0 is independent of h and w. The above notations and properties can be extended to vector-valued functions. To save space, we omit them.

2.4 Projections Several special projections will be used in this paper.

123

J Sci Comput

1. The Gauss–Radau [5] projection πh+ : for any scalar function w ∈ H 1 (), the projection is the unique element in Whk , satisfying  

Jj

(πh+ w − w, v) K i j = 0,   (πh+ w − w) xi− 1 , y v(y)dy = 0, 2

  (πh+ w − w) x, y j− 1 v(x)dx = 0, 2 Ii   (πh+ w − w) xi− 1 , y j− 1 = 0, 2

∀v ∈ Q k−1 (K i j ),

(2.7a)

∀v ∈ P k−1 (J j ),

(2.7b)

∀v ∈ P k−1 (Ii ),

(2.7c) (2.7d)

2

for any i = 1, 2, . . . , N x and j = 1, 2, . . . , N y . 2. The L 2 -projection P h : for any w ∈ [L 2 ()]2 , the projection P h w is the unique element . in W kh = [Whk ]2 such that ( P h w − w, v) = 0, ∀v ∈ W kh .

(2.8)

3. The projection h following [9]: for any vector-valued function w ∈ [H 1 ()]2 , the projection h w is the unique element in W kh , such that  

Jj

Ii

(h w − w, ∇v) K i j = 0,   (h w − w) xi+ 1 , y · nv(y)dy = 0,

∀v ∈ Q k (K i j ),

(2.9a)

∀v ∈ P k (J j ),

(2.9b)

  (h w − w) x, y j+ 1 · nv(x)dx = 0,

∀v ∈ P k (Ii ),

(2.9c)

2

2

hold for any i = 1, 2, . . . , N x and j = 1, 2, . . . , N y , where n is the unit normal on the interfaces, pointing to the right or above. For the special projections mentioned above, there hold the standard approximation properties [5,9]. Lemma 2.1 Suppose w ∈ H k+1 () and w ∈ [H k+1 ()]2 , then there hold πh+ w − w + h 1/2 πh+ w − w h + h πh+ w − w ∞ ≤ Ch k+1 w k+1 ,

(2.10)

w k+1 ,

(2.11)

w k+1 ,

(2.12)

P h w − w + h

1/2

P h w − w h + h P h w − w ∞ ≤ Ch

h w − w + h

1/2

h w − w h + h h w − w ∞ ≤ Ch

k+1

k+1

where the bounding constant C > 0 solely depends on the regularity of w or w. Moreover, the projection πh+ on cartesian meshes has the following superconvergence property (see [7], Lemma 3.6). It is very important in obtaining the optimal error estimate in the L 2 -norm. Lemma 2.2 Suppose w ∈ H k+2 (), then for any K and w ∈ W kh there holds |(w − πh+ w, ∇ · w) K − w − πh+ w, w · n K ∂ K | ≤ Ch k+1 w k+2,K w K ,

(2.13)

where n K is the outward normal of K , and the bounding constant C > 0 is independent of K and h.

123

J Sci Comput

3 IMEX-LDG Schemes and Main Results In this section we present two fully-discrete LDG schemes. To this end, we start from the semi-discrete scheme, following [15] and using the same notations.

3.1 Semi-discrete LDG Scheme As usual, we would like to introduce two auxiliary variables s = −∇c and z = D(u)s, and represent (1.1) into an equivalent first order system: φct + ∇ · (uc) + ∇ · z = cq, ˜

(3.1a)

s = −∇c,

(3.1b)

z = D(u)s,

(3.1c)

A(c)u = −∇ p,

(3.1d)

∇ · u = q.

(3.1e)

Multiply (3.1a)–(3.1e) by test functions v, ζ ∈ Whk , and w, ψ, θ ∈ W kh , respectively, and integrate by parts in each element K . Then we get (φct , v) K = (uc + z, ∇v) K − (uc + z) · n K , v ∂ K + (cq, ˜ v) K , (s, w) K = (c, ∇ · w) K − c, w · n K ∂ K , (z, ψ) K = (D(u)s, ψ) K , (A(c)u, θ ) K = ( p, ∇ · θ ) K −  p, θ · n K ∂ K , −(u, ∇ζ ) K = − u · n K , ζ ∂ K + (q, ζ ) K , where n K is the outward normal of K . Replacing the exact solutions by the numerical solutions, and then introducing reasonable numerical fluxes at every element interface, we can define the semi-discrete LDG scheme as follows: find for t ∈ (0, T ] the numerical solutions ch , ph ∈ Whk and s h , z h , uh ∈ W kh , such that, in any element K , the variational forms (φch t , v) K = LcK (uh , ch , v) + LdK (z h , v) + (c˜h q, v) K ,

(3.2a)

(s h , w) K = D K (ch , w),

(3.2b)

(z h , ψ) K = (D(uh )s h , ψ) K ,

(3.2c)

(A(ch )uh , θ ) K = D K ( ph , θ ), (q, ζ ) K =

− LdK (uh , ζ ),

(3.2d) (3.2e)

hold for any test functions v, ζ ∈ Whk , and w, ψ, θ ∈ W kh . In the above formulations,  c , q > 0, (3.3) c˜h = ∗ ch , q < 0, due to (1.2), and three functionals LcK (s, c, v) = (sc, ∇v) K −  sc · n K , v ∂ K , LdK (s, v)

= (s, ∇v) K −  s · n K , v ∂ K ,

D K (c, w) = (c, ∇ · w) K −  c , w · n K ∂ K ,

(3.4a) (3.4b) (3.4c)

123

J Sci Comput

are used for notational convenience. In this paper we adopt the alternating numerical fluxes [5,29] for the diffusion term, namely + h = u− zh = z − ph = ph+ . h = ch , u h , c h,

(3.5)

For the convection term, we use the Lax-Friedrichs-type flux [21], namely u  h ch =

1 + + − (u c + u− h ch − αne [[ch ]]), 2 h h

(3.6)

where α ≥ 0 can be chosen as any constant and ne is the normal of e ∈ h pointing to the right or above. Obviously, the above numerical fluxes are consistent and Lipschitz continuous. The initial condition is given as ch (x, y, 0) = πh+ c0 (x, y),

(3.7)

where πh+ is the Gauss–Radau projection defined in (2.7), and c0 (x, y) is the given initial solution. Till now we complete the definition of the semi-discrete LDG scheme. Remark 3.1 For the above semi-discrete LDG scheme, Guo et al. [15] derived the following optimal error estimate c − ch L ∞ (L 2 ) + s − s h L 2 (L 2 ) + u − uh L ∞ (L 2 ) ≤ Ch k+1 ,

(3.8)

if the exact solutions are smooth enough, for example, c ∈ L ∞ (H k+2 ), ct ∈ L ∞ (H k+1 ), p ∈ L ∞ (H k+2 ) and u ∈ [L ∞ (H k+1 )]2 . Here the bounding constant C > 0 is independent of the mesh size h. For convenience of notations and analysis, we would like to write the semi-discrete LDG scheme into a global form. By denoting    Lc (s, c, v) = LcK (s, c, v), Ld (s, v) = LdK (s, v), D(c, w) = D K (c, w), K ∈h

K ∈h

K ∈h

the global form is almost the same as (3.2), except that the subscript K is dropped.

3.2 Implicit–Explicit Fully-Discrete LDG Schemes M be a uniform partition of the time interval [0, T ], with time step τ = T /M. Let {t n = nτ }n=0 The time step could actually change from step to step, but in this paper we take the time step as a constant for simplicity. Two multi-step IMEX time marching scheme coupled with LDG spatial discretization will be considered. For convenience, we use the same notations for two schemes.

3.2.1 First Order Time-Marching Scheme The first order time-marching scheme, denoted by IMEX-LDG(k,1) in this paper, is given as follows. For any n ≥ 0, suppose that the numerical solutions chn , phn , snh , z nh and unh are and z n+1 available, we would like to find the numerical solutions chn+1 , sn+1 h h , such that

 chn+1 − chn n n (3.9) φ , v = Lc (unh , chn , v) + Ld (z n+1 h , v) + (c˜h q , v), τ

123

J Sci Comput

together with n+1 (sn+1 h , w) = D (ch , w),

(z n+1 h , ψ)

= (D(un+1, )sn+1 h h , ψ),

(3.10a) (3.10b)

hold for any v ∈ Whk , and w, ψ ∈ W kh . In (3.9), the forward Euler discretization is used for convection term and backward Euler discretization is used for diffusion term. Noting that the diffusion–dispersion tensor D(·) in (3.10b) is explicitly linearized in the form un+1, = unh , n ≥ 0. h

(3.11a)

For convenience of notations, we would like to supplement a definition for n = −1, namely, 0 u0, h = uh .

(3.11b)

After that, we compute the numerical solution un+1 and phn+1 , by the following variational h forms n+1 (A(chn+1 )un+1 h , θ ) = D ( ph , θ ),

(q

n+1

,ζ) = − L

d

(un+1 h , ζ ),

(3.12a) (3.12b)

which hold for any ζ ∈ Whk and θ ∈ W kh . To start this scheme, we would like to take ch0 = πh+ c0 , same as the semi-discrete case, and then compute u0h and ph0 by the variation form (3.12) with n = −1. Note that we do not need to compute s0h and z 0h for this scheme.

3.2.2 Second Order Time-Marching Scheme The second order time-marching scheme, denoted by IMEX-LDG(k,2) in this paper, is defined as follows. For any n ≥ 1, suppose the numerical solutions at t n and t n−1 have been known, we would like to find the numerical solutions chn+1 , sn+1 and z n+1 h h , such that 

cn+1 − chn 1 3 d n+1 3 n−1 , v = Lc (unh , chn , v) − Lc (un−1 φ h h , ch , v) + L (z h , v) τ 2 2 4 (3.13) 1 d n−1 3 n n 1 n−1 n−1 + L (z h , v) + (c˜h q , v) − (c˜h q , v), 4 2 2 together with (3.10), hold for any v ∈ Whk , and w, ψ ∈ W kh . In (3.13), the standard AdamsBashforth extrapolation and a kind of modified Adams-Moulton interpolation are used for the explicit discretization of convection and the implicit discretization of diffusion, respectively; see [14]. To maintain the second order accuracy in time, the linear extrapolation un+1, = 2unh − un−1 h h ,

(3.14a)

is used in (3.10), instead of (3.11a). After that, we use the variation form (3.12) to compute and phn+1 . the numerical solutions un+1 h To start this scheme, we need to provide initial solutions at two time levels. Same as that in the previous subsection, we take ch0 = πh+ c0 , and compute u0h and ph0 by the variation form (3.12) with n = −1. As a main difference, now we need to compute the numerical solutions s0h and z 0h by the variation form (3.10) with n = −1 and 0 u0, h = uh .

(3.14b)

123

J Sci Comput

To obtain the initial solution at t 1 , we would like to employ the first order time-marching scheme introduced in the previous subsection. Namely, we firstly compute the numerical solutions ch1 , s1h and z 1h by the variation forms (3.9) and (3.10) with n = 0 and 0 u1, h = uh .

(3.14c)

After that, we compute u1h and ph1 by the variation form (3.12) with n = 0. Till now we complete the definition of the second order time-marching scheme.

3.2.3 Main Result To obtain the optimal error estimates of the IMEX-LDG(k, r ) scheme, where r = 1, 2, we assume that the exact solutions of (3.1) are smooth enough, for example, c, p ∈ L ∞ (H k+2 ) and u, s, z ∈ [L ∞ (H k+1 )]2 ,

(3.15a)

as well as ( )

ct

(r +1)

∈ L ∞ (H k+1 ) for = 1, · · · , r, and ct

∈ L ∞ (L 2 ).

(3.15b)

( )

Here ct means the -th order derivative of c in the time direction. For r = 2, we require in addition that both ut and utt belong to [L ∞ (W 1,∞ )]2 . For any n ≥ 0, denote the numerical error at t n by en = (ecn , enu , ens , enz , enp ) = (cn − chn , un − unh , sn − snh , z n − z nh , p n − phn ),

(3.16)

where (cn , un , sn , z n , p n ) are the reference functions at the time level t n , which are defined as cn = c(t n ), un = u(t n ), sn = s(t n ),

z n = D(un, )s(t n ),

p n = p(t n ).

Here un, is defined corresponding to the adopted linearization technique, namely un−1 , n ≥ 1 n, u = n = 0, u0 , for r = 1, and

u

n,

=

2un−1 − un−2 , n ≥ 2 u0 , n = 0, 1

(3.17)

(3.18)

(3.19)

for r = 2. Note that z n = D(un )s(t n ), namely, the reference function is not defined as the exact solution. It is an important trick to obtain the optimal error estimate in time, especially for the second order time-marching scheme, since the reference function z n preserves the same structure as that in the used time-marching. In this paper we would like to use the symbol C to denote a generic constant independent of n, h and τ , which may have different values at different occurrences. Now we present the main result as follows. Theorem 3.1 Let k ≥ 1. Assume the temporal-spatial condition, for r = 1, 2, that  1+2δ  τ =O h r ,

123

(3.20)

J Sci Comput

where δ > 0 is any given constant. Then the IMEX-LDG(k, r ) scheme satisfies the optimal error estimate ecn 2 + enu 2 + τ

n 

2 2k+2 em + τ 2r ), ∀n ≥ 1, s ≤ C(h

(3.21)

m=1

provided the mesh parameter h is small enough. Here the bounding constant C > 0 is independent of n, h and τ . Below we present a snapshot of the proof. As the traditional treatment in the finite element analysis, we split the numerical error into two parts, namely en = ξ n − η n ,

(3.22)

for n = 0, 1, . . . , M, where ξ n = (ξcn , ξ nu , ξ ns , ξ nz , ξ pn ) = (πh+ cn − chn , h un − unh , P h sn − snh , h z n − z nh , πh+ p n − phn ), η = n

=

(ηcn , ηnu , ηns , ηnz , ηnp ) (πh+ cn − cn , h un −

(3.23a)

un , P h sn − sn , h z n − z n , πh+ p n − p n ),

(3.23b)

with projections πh+ , P h and h being defined in Sect. 2.4. From Lemma 2.1 and the linear structure of these projections, we have the following approximation properties. Lemma 3.1 The projection errors satisfy, for any n ≥ 0, the following properties ηcn + ηnu + ηns + ηnz + ηnp ≤ Ch k+1 , ηcn h + ηnu h ≤ Ch

k+ 21

ηcn ∞

k

+ ηnu ∞

,

≤ Ch ,

(3.24a) (3.24b) (3.24c)

and ηcn+1 − ηcn ≤ Ch k+1 τ,

(3.25)

where the bounding constant C > 0 is independent of n, h and τ . Thanks to the triangle inequality, Theorem 3.1 can be proved by establishing the optimal estimate of ξ n . In this process, in order to conveniently deal with the difficulties resulted from nonlinearity, we follow [32,33] and make the a priori assumption: if h is small enough, there holds ecm ≤ h 1+δ , 0 ≤ m ≤ n, (3.26) for any n ≥ 0. After setting up the error equations with respect to ξ n ∈ [Whk ]8 , we can obtain the following lemma by carrying out the energy analysis. Lemma 3.2 Assume the temporal-spatial condition (3.20) holds. Under the a prior assumption (3.26), the IMEX-LDG(k, r ) scheme (r = 1, 2) satisfies ξcn 2 + ξ nu 2 + τ

n 

2 2k+2 ξ m + τ 2r ), ∀n ≥ 1, s ≤ C(h

(3.27)

m=1

where the bounding constant C > 0 is independent of n, h and τ .

123

J Sci Comput

The proof of this lemma for two schemes will be given respectively in the next two sections. To that end, we present here some elemental conclusions about the DG spatial discretization, which will be used many times. 1. Integrating by parts on each cell, it is easy to check the following identity: for any functions v ∈ H 1 (h ) and w ∈ [H 1 (h )]2 , there holds Ld (w, v) + D(v, w) = 0.

(3.28)

2. Based on the definitions of the projection h and the operator Ld , it is easy to see that Ld (ηnu , v) = Ld (ηnz , v) = 0, ∀v ∈ Whk ,

(3.29)

holds for any n ≥ 0. 3. Moreover, it follows from Lemma 2.2 and the definition of D(·, ·) that |D(ηcn , w)| ≤ Ch k+1 cn k+2 w , |D(ηnp , w)| ≤ Ch k+1 p n k+2 w ,

(3.30)

for any w ∈ W kh and any n ≥ 0. We refer the readers to [15] for more details. Finally, we need to show the reasonability of the a priori assumption (3.26), in order to complete the proof of Theorem 3.1. In fact, it can be easily verified by a simple induction. Since ξc0 = 0, due to the setting of initial solution, it follows from k ≥ 1 and the approximation property (3.24) that ec0 = ηc0 ≤ Ch k+1 ≤ h 1+δ , if h is small enough. Assume (3.26) holds for any given n, then (3.27) holds with the bounding constant independent of n and h, which together with (3.24) and (3.20) imply ecn+1 ≤ C(h k+1 + τ r ) ≤ h 1+δ , if h is small enough. Hence, the a priori assumption (3.26) is reasonable. Remark 3.2 The condition (3.20) is acceptable, if the order of accuracy in space and in time are balanced for k ≥ 1. This condition is just used for the a priori assumption. And the requirement k ≥ 1 is just for the theoretical analysis, we could also choose k = 0 in numerical experiments.

4 Proof of Lemma 3.2 with r = 1 In this section, we consider the first order time-marching scheme.

123

J Sci Comput

4.1 Error Equation Owing to the smoothness assumption (3.15), it can be verified that the reference functions satisfy, for any n ≥ 0, the following variational forms

n+1  − cn c φ (4.1a) , v = Lc (un , cn , v) + Ld (z n+1 , v) + (c˜n q n , v) + (ς n , v), τ (sn+1 , w) = D(cn+1 , w), (z

n+1

, ψ) = (D(u

n+1,

)s

(4.1b)

n+1

, ψ),

(4.1c)

(A(c )u , θ ) = D( p , θ ), n

n

n

(4.1d)

(q , ζ ) = − L (u , ζ ), n

for any v, ζ ∈

Whk ,

d

and w, ψ, θ ∈

n

W kh .

(4.1e)

Here ς n is the local truncation error satisfying

ς n ≤ Cτ, ∀n ≥ 0,

(4.2)

where the bounding constant C > 0 depends on the regularity of c, ct and ctt . Here the boundedness of D(u) and its first order derivative (hypothesis 4) is used. Subtracting (4.1) from the first order time-marching scheme, we get the following error equations

n+1  − ecn e n n n φ c , v = F (un , cn ; unh , chn ; v) + Ld (en+1 z , v) + (e˜c q , v) + (ς , v), τ (4.3a) , w) = D(ecn+1 , w), (en+1 s (A(cn )un −

(en+1 z , ψ) A(chn )unh , θ ) Ld (enu , ζ )

which hold for any v, ζ ∈

= (D(u

n+1,

)s

(4.3b)

n+1

− D(un+1, )sn+1 h h , ψ),

(4.3c)

= D(enp , θ ),

(4.3d)

= 0,

Whk ,

(4.3e)

and w, ψ, θ ∈

W kh .

Here

F (un , cn ; unh , chn ; v) = Lc (un , cn , v) − Lc (unh , chn , v)

and

 n

e˜c =

(4.4)

0, q n > 0, ecn , q n < 0.

(4.5)

Furthermore, according to the error splitting and property (3.29), we derive, for any n ≥ 0, the following error equations  n+1 

n+1 − ξcn − ηcn η ξ ,v = φ c , v + F (un , cn ; unh , chn ; v) φ c τ τ n n n + Ld (ξ n+1 z , v) + (e˜c q , v) + (ς , v),

(A(cn )un −

, w) (ξ n+1 s n+1 (ξ z , ψ) A(chn )unh , θ ) Ld (ξ nu , ζ )

which hold for any v, ζ ∈

(4.6a)

= (ηn+1 , w) + D(ξcn+1 , w) − D(ηcn+1 , w), s n+1, n+1 = (ηn+1 )s − D(un+1, )sn+1 z , ψ) + (D(u h , ψ), h n n = D(ξ p , θ ) − D(η p , θ ),

(4.6d)

= 0,

(4.6e)

Whk ,

and w, ψ, θ ∈

(4.6b) (4.6c)

W kh .

123

J Sci Comput

4.2 Elemental Inequalities We present some elemental inequalities in this subsection. The next lemma establishes the important relationship between the gradient and jump of numerical solution with the numerical solution of the gradient. It can be easily derived by following [26] with minor changes, so the detailed proof is omitted here. Lemma 4.1 The variation form (4.6b) implies for n ≥ 0 that 1

∇ξcn+1 + h − 2 [[ξcn+1 ]] h ≤ C( ξ n+1 + h k+1 ), s

(4.7)

where the bounding constant C > 0 is independent of n and h. Owing to the error Eq. (4.6), we have the following lemmas. Lemma 4.2 Under the a priori error assumption (3.26), we have ξ nu ≤ C(h k+1 + ξcn ), ∀n ≥ 0,

(4.8)

where the bounding constant C > 0 is independent of n, h and τ . Proof Taking θ = ξ nu and ζ = ξ pn in (4.6d) and (4.6e), respectively, and using (3.28), we obtain   A(cn )ξ nu , ξ nu = T1 + T2 + T3 + T4 , (4.9) where

T1 = (A(cn )ηnu , ξ nu ), T2 = −((A(cn ) − A(chn ))un , ξ nu ), T3 = − D(ηnp , ξ nu ), T4 = −((A(cn ) − A(chn ))(unh − un ), ξ nu ).

(4.10)

Using Cauchy–Schwarz inequality and Lemma 3.1, we can get |T1 | ≤ C ηnu ξ nu ≤ Ch k+1 ξ nu ,

(4.11)

since A(·) is bounded (see (2.1)). By hypothesis 2, Cauchy–Schwarz inequality, the Lipschitz continuity of A(·) and Lemma 3.1, we have   |T2 | ≤ C A(cn ) − A(chn ) ξ nu ≤ C ecn ξ nu ≤ C h k+1 + ξcn ξ nu . (4.12) For T3 , we use (3.30) to obtain |T3 | ≤ Ch k+1 p n k+2 ξ nu .

(4.13)

By the Lipschitz continuity of A(·), the triangle inequality, the inverse inequality, Lemma 3.1 and the a priori assumption (3.26), we can derive A(cn ) − A(chn ) ∞ ≤ C ecn ∞ ≤ C( ξcn ∞ + ηcn ∞ ) ≤ C(h −1 ξcn + h k ) ≤ C(h −1 ecn + h k ) ≤ Ch δ .

(4.14)

Hence, letting h be small enough, we can get from (2.1) that 1 1 μ∗ /κ ∗ ≤ A(cn ). 4 4 An application of Cauchy–Schwarz inequality and Lemma 3.1 yields 1  |T4 | ≤ A(cn )ξ nu 2 + Ch 2k+2 . 2 A(cn ) − A(chn ) ∞ ≤

123

(4.15)

J Sci Comput

Substituting (4.11), (4.12), (4.13) and (4.15) into (4.9), we have the estimate    A(cn )ξ nu 2 ≤ C h k+1 + ξcn ξ nu + Ch 2k+2 ,  

which further yields (4.8) owing to (2.1). Together with Lemma 4.2, the a priori assumption (3.26) implies δ ecm ∞ + em u ∞ ≤ Ch , for 0 ≤ m ≤ n, ∀n ≥ 0,

(4.16)

by a similar argument as (4.14). It further implies chm ∞ + um h ∞ ≤ C, for 0 ≤ m ≤ n, ∀n ≥ 0,

(4.17)

by the boundedness of c and u (hypothesis 2). Then it follows from hypothesis 4 that D(unh ) ∞ ≤ C,

(4.18)

where D ∞ has been defined in (2.2) with s = 0 and D = . Lemma 4.3 Under the a priori error assumption (3.26), we have for any n ≥ 0 that n n+1 ξ n+1 + h k+1 ), z ≤ C( ξ u + ξ s

(4.19)

where the bounding constant C > 0 is independent of n, h and τ . Proof Taking ψ = ξ n+1 in (4.6c) we get z 2 n+1 n+1 n+1, n+1 n+1 ξ n+1 )s − D(un+1, )sn+1 z = (η z , ξ z ) + (D(u h h , ξz ) n+1 n n n+1 n+1 n+1 , ξ z ) + (D(unh )(sn+1 − sn+1 = (ηn+1 z , ξ z ) + (D(u ) − D(uh ))s h ), ξ z ).

Then we have

  2 n+1 n n n+1 n+1 ξ n+1 + ξ n+1 ) ξ n+1 z ≤ η z + D(u ) − D(uh ) ξ z + C( η s s z   ξ n+1 ≤ C h k+1 + ξ nu + ξ n+1 s z ,

where in the first step we have used Cauchy–Schwarz inequality, the boundedness of s = −∇c assumed in hypothesis 2 and (4.18), the second step follows from the Lipschitz continuity of diffusion–dispersion tensor (implied by hypothesis 4) and Lemma 3.1. Thus canceling ξ n+1   z in the above equation yields (4.19).

4.3 Final Estimate We continue to derive the final estimate. Taking the test function v = ξcn+1 , w = ξ n+1 z , and in (4.6a), (4.6b) and (4.6c), respectively, adding them together, and using (3.28), ψ = −ξ n+1 s we obtain   1 )ξ n+1 , ξ n+1 (φ(ξcn+1 − ξcn ), ξcn+1 ) + D(un+1, s s h τ (4.20) n+1 n+1 n+1 + Rp2 + Rp3 , = Rcn+1 + Rdn+1 + Rp1

123

J Sci Comput

where Rcn+1 = F (un , cn ; unh , chn ; ξcn+1 ),     n+1, n+1 n+1 n+1, n+1 n+1 Rdn+1 = D(un+1, − (D(u , )η , ξ ) − D(u ))s , ξ s s s h h

(4.21b)

n+1 Rp1

(4.21c)

n+1 Rp2 n+1 Rp3

n+1 n+1 n+1 = − D(ηcn+1 , ξ n+1 , ξ z ) − (ηn+1 ), z ) + (η s z , ξs

1 = (φ(ηcn+1 − ηcn ), ξcn+1 ) + (e˜c n q n , ξcn+1 ), τ = (ς n , ξcn+1 ).

(4.21a)

(4.21d) (4.21e)

Next we estimate the above terms one by one. Firstly, we notice that  n cn ) · n , [[ξ n+1 ]] ,  Rcn+1 = (un cn − unh chn , ∇ξcn+1 ) − (un cn − u e e c h h

(4.22)

e∈h \∂+

where ∂+ means the right and top boundary of . Proceeding the same argument as the estimate for R3 in [15], we can derive    k+1 , (4.23) + h |Rcn+1 | ≤ C h k+1 + ξ nu + ξcn ξ n+1 s where Lemma 4.1 plays an important role. For more details please refer to [15]. Using Cauchy–Schwarz inequality, (4.18), the Lipschitz continuity of diffusion–dispersion tensor (implied by hypothesis 4), and Lemma 3.1, we can get ξ n+1 + C D(un ) − D(unh ) ξ n+1 |Rdn+1 | ≤ C ηn+1 s s s ≤ Ch k+1 ξ n+1 + C un − unh ξ n+1 s s   . ≤ C h k+1 + ξ nu ξ n+1 s

(4.24)

By Cauchy–Schwarz inequality, Lemma 3.1, (3.30), hypothesis 1 and (4.2), we have n+1 n+1 n+1 + Rp2 | ≤ Ch k+1 ( ξ n+1 ) + C(h k+1 + ξcn ) ξcn+1 , |Rp1 z + ξ s n+1 | |Rp3

≤ Cτ ξcn+1 .

(4.25a) (4.25b)

Substituting (4.23)–(4.25) into (4.20) and using Lemmas 4.2 and 4.3, we obtain  1  1  n+1 2 1  φξc + φ(ξcn+1 − ξcn ) 2 − φξcn 2 + τ D(un+1, )ξ n+1 2 s h 2 2 2 ≤ Cτ (h k+1 + ξcn )(h k+1 + ξ n+1 ) + Cτ (h k+1 + τ + ξcn ) ξcn+1 s  1 )ξ n+1 2 + Cτ ( ξcn 2 + ξcn+1 2 ) + Ch 2k+2 τ + Cτ 3 , (4.26) ≤ τ D(un+1, s h 2 by using Young’s inequality. Here we have used the uniformly positive definite property (hypothesis 4) of diffusion–dispersion tensor. Summing (4.26) and then using the discrete Gronwall’s inequality, we get ξcn 2 + τ

n−1 

ξ m+1 2 ≤ C(h 2k+2 + τ 2 ), s

(4.27)

m=0

if τ is small enough, where we have used hypotheses 1. It further yields ξ nu ≤ C(h k+1 + τ ), by Lemma 4.2, and we complete the proof of this lemma for r = 1.

123

(4.28)

J Sci Comput

5 The Proof of Lemma 3.2 with r = 2 The proof line is similar as but a little complicated than that in the previous section.

5.1 Error Equation and Some Elemental Inequalities Owing to the smoothness assumption (3.15), it can be verified that the exact solution satisfies 

n+1 3 c 1 3 − cn , v = Lc (un , cn , v) − Lc (un−1 , cn−1 , v) + Ld (z n+1 , v) φ τ 2 2 4 1 3 1 + Ld (z n−1 , v) + (c˜n q n , v) − (c˜n−1 q n−1 , v) + (ς n , v), (5.1) 4 2 2 for any n ≥ 1, and satisfies (4.1b)–(4.1e) for any n ≥ 0. Here ς n is the local truncation error which satisfies Cτ 2 , n ≥ 3 and n = 1, n ς ≤ (5.2) Cτ, n = 2, where the bounding constant C > 0 depends on the regularity of ut , utt and c, up to the third order time derivative. Here the boundedness of D(u) up to the second order derivatives (hypothesis 4) is used. Note that the order-reduction for n = 2 is resulted from the first order linearization at time t 1 . Then we can get the following error equation  n+1 

n+1 ηc − ηcn 3 ξc − ξcn ,v = φ , v + F (un , cn ; unh , chn ; v) φ τ τ 2 1 n−1 − F (un−1 , cn−1 ; un−1 h , ch ; v) 2 1 d n−1 3 + Ld (ξ n+1 z , v) + L (ξ z , v) 4 4 1 n−1 n−1 3 n n q , v) + (ς n , v), (5.3) + (e˜c q , v) − (e˜c 2 2 for n ≥ 1, and the error Eq. (4.6a) for n = 0, with ς 0 ≤ Cτ.

(5.4)

Along the same line, we also have for m ≥ 0 the following error equations m m m (ξ m s , w) = (η s , w) + D (ξc , w) − D (ηc , w),

(5.5a)

m m, m m )s − D(um, (ξ m z , ψ) = (η z , ψ) + (D(u h )s h , ψ),

(5.5b)

m m ((A(cm )um − A(chm )um h ), θ ) = D (ξ p , θ ) − D (η p , θ ), d

L

(ξ m u ,ζ)

= 0,

(5.5c) (5.5d)

m, have been defined in (3.14) and (3.19), respectively. where um, h and u For the IMEX-LDG(k, 2) scheme, Lemma 4.2 still holds under the a priori assumption (3.26), namely k+1 ξ m + ξcm ), ∀m ≥ 0. (5.6) u ≤ C(h

Hence, the inequalities from (4.16) to (4.18) are also true when h is small enough, owing to the a priori assumption (3.26). However, the estimate to ξ m z has a little modification, because the linearization methods are not the same. The detailed conclusion is presented in the following lemma without proof, since the proof line is the same as Lemma 4.3.

123

J Sci Comput

Lemma 5.1 Under the a priori error assumption (3.26), we have m−1 k+1 ξ m + ξ m−2 + ξ m ), m ≥ 2, z ≤ C( ξ u u s +h

(5.7a)

m−1 k+1 ξ m + ξ m ), m = 1, z ≤ C( ξ u s +h

(5.7b)

m m k+1 ξ m ), m = 0. z ≤ C( ξ u + ξ s + h

(5.7c)

and

and

5.2 Energy Estimate In what follows, let n ≥ 1. We first take the test function v = ξcn+1 in (5.3), and then take the test function w = 43 ξ n+1 in (5.5a) and ψ = − 43 ξ n+1 in (5.5b), with m = n + 1. Adding z s the resulted equations together, we obtain  3  1 φ(ξcn+1 − ξcn ), ξcn+1 + D(un+1, )ξ n+1 , ξ n+1 s s h τ 4 n+1 n+1 n+1 n+1 n+1 + Rd2 + Rp1 + Rp2 + Rp3 , = Rcn+1 + Rd1

(5.8)

where 3 1 n−1 n+1 Rcn+1 = F (un , cn ; unh , chn ; ξcn+1 ) − F (un−1 , cn−1 ; un−1 ), h , ch ; ξc 2 2     3 3 n+1, n+1 n+1 n+1 n+1, n+1 n+1 D(un+1, (D(u − , = )η , ξ ) − D(u ))s , ξ Rd1 s s s h h 4 4 1 n+1 n+1 = Ld (ξ n−1 ), Rd2 z , ξc 4 3 3 n+1 n+1 3 n+1 Rp1 = − D(ηcn+1 , ξ n+1 , ξ z ) − (ηn+1 , ξ n+1 ), z ) + (η s s 4 4 4 z 1 3 1 n+1 = (φ(ηcn+1 − ηcn ), ξcn+1 ) + (e˜c n q n , ξcn+1 ) − (e˜c n−1 q n−1 , ξcn+1 ), Rp2 τ 2 2 n+1 = (ς n , ξcn+1 ). Rp3 Below we will estimate them separately. n+1 are similar as those in Sect. 4.2. Hence, we omit the The estimates for Rcn+1 and Rd1 detailed process and directly present the final results:    n n−1 n+1 k+1 |Rcn+1 | ≤ C h k+1 + ξ nu + ξ n−1 + ξ + ξ ξ + h u c c s    k+1 n n−1 n+1 k+1 ≤C h (5.9a) + ξc + ξc ξ s + h   n+1 n+1 | ≤ C h k+1 + ξ nu + ξ n−1 |Rd1 u ξ s   , (5.9b) ≤ C h k+1 + ξcn + ξcn−1 ξ n+1 s

123

J Sci Comput

where (5.6) is used. Also, we have n+1 n+1 |Rp1 | ≤ Ch k+1 ( ξ n+1 ) ≤ Ch k+1 ( ξcn + ξcn−1 + ξ n+1 + h k+1 ), z + ξ s s

(5.9c) n+1 |Rp2 |

≤ C(h

k+1

+ ξcn + ξcn−1 ) ξcn+1 ,

(5.9d)

Cτ 2 ξcn+1 , n ≥ 3 and n = 1, n+1 n n+1 |Rp3 | ≤ ς ξc ≤ Cτ ξcn+1 , n = 2.

(5.9e)

where Lemma 5.1, (5.6) and (5.2) are used. n+1 is much more complex. Using property (3.28), and taking the test The estimate to Rd2 n−1 function w = ξ z in (5.5a) with m = n + 1, we can get that  1 1  n+1 n−1 n+1 n+1 n−1 Rd2 (ξ s , ξ z ) − (ηn+1 = − D(ξcn+1 , ξ n−1 , ξ n−1 z )=− s z ) + D (ηc , ξ z ) . 4 4 in (5.5b) with m = n − 1, we have Taking ψ = ξ n+1 s n+1 n+1 n+1 ) = (ηn−1 ) + (D(un−1, )sn−1 − D(un−1, )sn−1 ). (ξ n−1 z , ξs z , ξs h h , ξs

As a consequence, we have n+1 Rd2 = n+1 + n+1 1 2 ,

(5.10)

where

 1  n+1 n−1 n+1 n+1 n−1 (ηs , ξ z ) − (ηn−1 , ξ ) − D (η , ξ ) , (5.11a) z s c z 4 1 n+1 = − (D(un−1, )sn−1 − D(un−1, )sn−1 ). (5.11b) n+1 2 h h , ξs 4 Together with the superconvergence property (3.30), a simple use of Cauchy–Schwarz inequality and Lemma 3.1 yields n+1 = 1

k+1 n+1 ( ξ n−1 ). |n+1 z + ξ s 1 | ≤ Ch

Hence it follows from Lemma 5.1, inequality (5.6), and ξc0 = 0 that Ch k+1 ( ξcn−2 + ξcn−3 + ξ n−1 + ξ n+1 + h k+1 ), n ≥ 3 n+1 s s |1 | ≤ + ξ n+1 + h k+1 ), n = 1, 2, Ch k+1 ( ξ n−1 s s ≤ Ch

k+1

(5.12)

( ξcn−2 + ξcn−3 + ξ n−1 + ξ n+1 + h k+1 ), s s

here and below we use ξc−1 = ξc−2 = 0 for convenience of notations. The estimate to n+1 2 starts from the following separation n+1 n+1 = n+1 n+1 2 21 + 22 + 23 ,

(5.13)

where 1 n−1, n+1 ) − D(un−1, )]sn−1 , ξ n+1 ) s 21 = − ([D(u h 4 1 n−1, n−1 n+1 )ηs , ξ s ), n+1 22 = (D(uh 4 1 n−1, n−1 n+1 )ξ s , ξ s ). n+1 23 = − (D(uh 4

(5.14a) (5.14b) (5.14c)

123

J Sci Comput

By the boundedness assumption of s = −∇c (hypothesis 2) and the Lipschitz continuity of diffusion–dispersion tensor (implied by hypothesis 4), we have k+1 + ξcn−2 + ξcn−3 ) ξ n+1 , |n+1 s 21 | ≤ C(h

(5.15)

where (5.6) is used. Applying (4.18) and Lemma 3.1 to n+1 22 , we get k+1 n+1 |n+1 ξ s . 22 | ≤ Ch

(5.16)

Applying Cauchy–Schwarz inequality and Young’s inequality to n+1 23 , we get |n+1 23 | ≤

1 1  n+1 , (5.17) )ξ n−1 , ξ n−1 ) + (D(un−1, )ξ n+1 , ξ n+1 ) = n+1 +  (D(un−1, s s s s h h 8 8

where 

n+1

 n+1 

    1 n−1, n−1 2 n+1, n+1 2 D(uh = )ξ s + D(uh )ξ s , 8  1 . = (D(un−1, ) − D(un+1, ))ξ n+1 , ξ n+1 s s h h 8

(5.18a) (5.18b)

Noticing D(un−1, ) = D(un+1, ) − [D(un+1, ) − D(un−1, )] h h )] − [D(un−1, ) − D(un−1, )], + [D(un+1, ) − D(un+1, h h we can further get  n+1 | ≤ C τ + |

3 

 en−i u ∞

2 ≤ C(τ + h δ ) ξ n+1 2 , ξ n+1 s s

i=0

under the a priori assumption (3.26), where the Lipschitz continuity of diffusion–dispersion tensor (implied by hypothesis 4) and the inequality (4.16) related to the a priori assumption are used. Therefore, if h is small enough, we can get  1 n+1 | ≤ D(un+1, )ξ n+1 2 , (5.19) | s h 8 due to the uniformly positive definite property (hypothesis 4) of diffusion–dispersion tensor. Collecting up the above estimations yields  1 n+1 |Rd2 | ≤ D(un+1, )ξ n+1 2 + n+1 s 8  h (5.20)   n+1 + C h k+1 + ξcn−2 + ξcn−3 h k+1 + ξ n−1 + ξ . s s Consequently, we obtain from (5.8), (5.9) and (5.20) that  1  n+1 2 1  n 2 5 )ξ n+1 2 φξc − φξc + τ D(un+1, s h 2 2 8 n+1 ≤ τ n+1 + τ ϒ n+1 + τ |Rp3 |,

(5.21)

for any n ≥ 1, where 3     ϒ n+1 = C h k+1 + ξcn−i h k+1 + ξcn+1 + ξ n−1 + ξ n+1 . s s i=0

123

(5.22)

J Sci Comput

Using Young’s inequality and the uniformly positive definite property (hypothesis 4) of diffusion–dispersion tensor, we can achieve 4    ϒ n+1 ≤ n+1 + C h 2k+2 + ξcn+1−i 2 .

(5.23)

i=0

Let m ≥ 3 be any integer. Summing the inequality (5.21) from n = 3 to n = m, and noticing n+1 |Rp3 | ≤ Cτ 2 ξcn+1 by (5.9e), we have m  1  m+1 2 1  D(un+1, )ξ n+1 2 φξc + τ s h 2 8 n=3

≤ R0 + C(h 2k+2 + τ 4 ) + Cτ

m+1 

ξcn 2 ,

(5.24)

n=0

where

    1  3 2 1 3, 3 2 2 2 )ξ + D(u )ξ φξc + τ D(u2, s s h h 2 4

  3   i 2 φξci 2 + τ D(ui, . ≤ Rstart := h )ξ s

R0 =

(5.25)

i=1

In the next subsection, we will prove Rstart ≤ C(h 2k+2 + τ 4 ).

(5.26)

As a consequence, noticinging hypotheses 1 and 4 again, and using the discrete Gronwall’s inequality yields ξcn 2 + τ

n−1 

ξ m+1 2 ≤ C(h 2k+2 + τ 4 ), ∀n ≥ 1, s

(5.27)

m=0

if τ is small enough. It further yields ξ nu ≤ C(h k+1 + τ 2 ),

(5.28)

by (5.6) and thus we get the desired result (3.27) with r = 2. Now we complete the proof of this lemma.

5.3 Proof of (5.26) In this subsection we prove (5.26). Firstly, since the second order scheme (3.13) starts with the first order scheme (3.9), taking n = 0 in (4.26) and noticing ξc0 = 0 leads to  1  1 2 1 2 k+1 k+1 φξc + τ D(u1, (h + ξ 1s ) + Cτ (h k+1 + τ ) ξc1 h )ξ s ≤ Cτ h 2 ≤ ε ξc1 2 + ετ ξ 1s 2 + C(h 2k+2 τ + τ 4 ), for arbitrary positive ε, where Young’s inequality is used in the last step. Due to hypotheses 1 and 4, we choose ε small enough, and get   1 2 2k+2 φξc1 2 + τ D(u1, + τ 4 ). (5.29) h )ξ s ≤ C(h

123

J Sci Comput

Secondly, letting n = 1 in (5.21) and noticing a rough estimate 2 τ |Rp3 | ≤ Cτ 2 ξc2 ≤ Cτ 4 + ε ξc2 2 ,

(5.30)

by (5.9e), where ε is any positive constant, we get  1  2 2 1  1 2 3 2 2 φξc − φξc + τ D(u2, h )ξ s 2 2 8  1 0 2 2 2 1 2 2k+2 ≤ τ D(u0, τ + Cτ 4 ). h )ξ s + C(τ + ε) ξc + Cτ ξc + C(h 4 Taking w = ξ 0s in (5.5a) with m = 0, by Cauchy–Schwarz inequality, Lemma 3.1 and (3.30) we can easily get ξ 0s ≤ Ch k+1 . Hence, taking ε small enough and substituting (5.29) into the above estimate, we have   2 2 2k+2 + τ 4 ), (5.31) φξc2 2 + τ D(u2, h )ξ s ≤ C(h if τ is small enough. Thirdly, we also have   3 2 2k+2 φξc3 2 + τ D(u3, + τ 4 ), h )ξ s ≤ C(h

(5.32)

along the same analysis line. Now we complete the proof of (5.26). 3 |. Remark 5.1 In (5.30), we only use a rough estimate that also holds for τ |Rp3

6 Numerical Experiments In this section, we provide some numerical examples to illustrate the accuracy and capability of the proposed methods. We will simulate problem (1.1) on  = (0, 1) × (0, 1) up to T = 1. In Examples 6.1 and 6.2, we present the optimal accuracy of the IMEX-LDG(k, r ) schemes. In Example 6.3, we display the evolution of concentration at different time for a practical problem. In all experiments, if not otherwise stated, we take κ(x, y) = 1 and μ(c) = 1, and let h = 1/N be the mesh size, where N = N x = N y is the number of rectangles in x- and y-direction, respectively. Example 6.1 Let φ(x, y) = 1, consider the problem with D(u) = γ I, where γ is a positive constant. In addition, we take the initial condition c0 = sin(2π(x + y)), and add a source term f = 2π[4π(x + y + t) + 1] cos(2π(x + y + t)) + 8γ π 2 sin(2π(x + y + t) in the concentration equation. The exact solution is c = sin(2π(x + y + t)), u = (4π x + 2πt, 4π y + 2πt) . The L 2 -norm errors and orders of accuracy for IMEX-LDG(k, r ) schemes are listed in k+1 Tables 1 and 2, for γ = 1 and 0.01. Here we take time step τ = 6h r . For both schemes, we can observe the optimal orders of accuracy for both c and u.

123

2

8.9543E−04

1.2193E−04

1.5991E−05

2.0389E−06

40

80

160

8.5112E−05

160

20

3.3971E−04

80

6.3254E−03

1.3512E−03

40

10

5.3213E−03

80

160

20

4.5916E−03

2.2968E−03

40

2.0481E−02

9.1787E−03

20

10

1.8330E−02

10

0

1

3.6572E−02

N

k

γ =1 c − ch

2.971

2.931

2.877

2.821

–

1.997

1.992

1.977

1.944

–

0.999

0.999

0.998

0.997

–

Order

2.9500E−06

2.3412E−05

1.8048E−04

1.3597E−03

9.9958E−03

2.0757E−04

8.2357E−04

3.2572E−03

1.2756E−02

4.8571E−02

1.1362E−02

2.2682E−02

4.5213E−02

8.9215E−02

1.7581E−01

u − uh

Table 1 Example 6.1: L 2 -norm errors and orders for IMEX-LDG(k, 1) schemes

2.988

2.946

2.913

2.878

–

1.988

1.984

1.969

1.929

–

0.997

0.995

0.981

0.979

–

Order

1.6567E−07

1.3105E−06

9.8843E−06

7.2474E−05

4.9179E−04

3.6512E−05

1.4532E−04

5.7164E−04

2.1971E−03

8.2681E−03

6.2209E−04

1.2430E−03

2.4812E−03

4.9468E−03

9.8571E−03

γ = 0.01 c − ch

2.984

2.915

2.874

2.763

–

1.993

1.976

1.942

1.912

–

0.999

0.997

0.995

0.995

–

Order

2.5465E−06

2.0244E−05

1.5867E−04

1.2123E−03

9.1476E−03

1.8326E−04

7.2959E−04

2.8795E−03

1.1206E−02

4.2869E−02

1.0871E−02

2.1686E−02

4.3197E−02

8.5569E−02

1.6813E−01

u − uh

2.991

2.971

2.934

2.916

–

1.993

1.981

1.960

1.936

–

0.996

0.994

0.986

0.974

–

Order

J Sci Comput

123

123

3

2.6168E−04

2.0105E−05

1.3297E−06

8.5473E−08

40

80

160

2.5141E−06

160

20

1.9569E−05

80

2.5869E−03

1.4955E−04

40

10

1.0998E−03

80

160

20

2.1090E−04

5.3011E−05

40

7.9712E−03

8.3748E−04

20

10

3.2987E−03

10

1

2

1.2848E−02

N

k

γ =1 c − ch

3.960

3.918

3.702

3.305

–

2.960

2.934

2.879

2.858

–

1.992

1.989

1.978

1.962

–

Order

6.7123E−08

1.0459E−06

1.6200E−05

2.3959E−04

3.4815E−03

2.6246E−06

2.0917E−05

1.6573E−04

1.3019E−03

1.0095E−02

3.2918E−04

1.3126E−03

5.2256E−03

2.0649E−02

8.1249E−02

u − uh

Table 2 Example 6.1: L 2 -norm errors and orders for IMEX-LDG(k, 2) schemes

3.962

3.953

3.887

3.861

–

2.994

2.986

2.974

2.955

–

1.995

1.993

1.982

1.976

–

Order

5.7628E−08

8.7924E−07

1.2679E−05

1.5989E−04

1.7582E−03

1.9637E−06

1.5341E−05

1.1632E−04

8.6593E−04

6.1647E−03

3.3256E−05

1.3272E−04

5.2857E−04

2.0826E−03

8.1857E−03

γ = 0.01 c − ch

3.931

3.850

3.657

3.459

–

2.966

2.923

2.896

2.832

–

1.997

1.994

1.978

1.975

–

Order

5.1514E−08

8.0623E−07

1.2400E−05

1.8823E−04

2.8518E−03

5.2683E−06

4.1546E−05

3.2237E−04

2.4739E−03

1.8712E−02

3.8311E−04

1.5285E−03

6.0268E−03

2.3577E−02

9.0185E−02

u − uh

3.968

3.943

3.924

3.921

–

2.979

2.956

2.940

2.919

–

1.996

1.979

1.968

1.935

–

Order

J Sci Comput

J Sci Comput Table 3 Example 6.2: L 2 -norm errors and orders for IMEX-LDG(k, 1) schemes k=0 Error

N c − ch

u − uh

Order

k=1 Error

Order

k=2 Error

Order

10

5.7159E−02

–

3.8471E−02

–

9.4716E−03

–

20

2.9851E−02

0.937

1.1321E−02

1.765

1.4754E−03

2.682

40

1.5179E−02

0.976

3.1512E−03

1.845

2.1219E−04

2.798

80

7.6916E−03

0.981

8.3971E−04

1.908

2.9041E−05

2.869

160

3.8868E−03

0.985

2.1511E−04

1.965

3.7389E−06

2.957

10

2.6872E−01

–

5.7164E−02

–

6.5718E−03

–

20

1.3892E−01

0.952

1.5756E−02

1.859

1.0597E−03

2.633

40

7.1241E−02

0.963

4.2057E−03

1.905

1.5215E−04

2.800

80

3.6192E−02

0.977

1.0936E−03

1.943

2.0151E−05

2.917

160

1.8198E−02

0.992

2.7757E−04

1.978

2.5776E−06

2.967

Order

k=3 Error

Order

Table 4 Example 6.2: L 2 -norm errors and orders for IMEX-LDG(k, 2) schemes k=1 Error

N c − ch

u − uh

Order

k=2 Error

10

1.8764E−02

–

9.1457E−03

–

4.7582E−03

–

20

5.2987E−03

1.824

1.2598E−03

2.860

4.0168E−04

3.566

40

1.4737E−03

1.846

1.6955E−04

2.893

2.9105E−05

3.787

80

3.9090E−04

1.915

2.2569E−05

2.909

1.9342E−06

3.911

160

1.0301E−04

1.924

2.9141E−06

2.953

1.2543E−07

3.947

10

9.0128E−02

–

1.7462E−02

–

6.8125E−03

–

20

2.3649E−02

1.930

2.3019E−03

2.923

5.8859E−04

3.533

40

6.1562E−03

1.942

2.9573E−04

2.960

4.3200E−05

3.768

80

1.5813E−03

1.961

3.7517E−05

2.979

2.8459E−06

3.924

160

4.0176E−04

1.977

4.7096E−06

2.994

1.8281E−07

3.960

 u 21 + 1 u 1 u 2 . In 2 u1u2 u2 + 1 addition, we choose the initial condition c0 = sin(2π(x + y)), and add a source term

Example 6.2 Let φ(x, y) = 1, consider the problem with D(u) =

f = 2π[4π(x + y + t))(1 − 12π) + 1] cos(2π(x + y + t)) + 4π 2 [16π 2 (x + y + t)2 + 2] sin(2π(x + y + t)) in the concentration equation. The exact solution is still c = sin(2π(x + y + t)), u = (4π x + 2πt, 4π y + 2πt) . The L 2 -norm errors and orders of accuracy of IMEX-LDG(k, r ) schemes are contained in k+1 Tables 3 and 4. Here we take time step τ = 10h r . We can also observe the optimal orders of accuracy for both schemes. Example 6.3 Let us consider a petroleum production scenario which happens in the real world with injection well and production well, treated as delta functions. The model is

123

J Sci Comput

Fig. 1 Example 6.3: Distribution of oil at different time. a T = 0.1, b T = 0.2, c T = 0.5, d T = 1.0

so-called quarter-of-a-five-spot pattern, which can be described by Eq. (1.1) subject to the non-permissible boundary condition. Let φ(x, y) = 0.1, and take D(u) = φ|u|I. The initial condition is taken as c0 = 1 which means there are only oil at the beginning time; we also take c˜ = 0 which means we inject the pure water into underground. The source term is defined as q = 0.05δ(x, y) − 0.05δ(x − 1, y − 1), to simulate injection well and production well located in the lower-left corner and upper-right corner respectively. In the numerical simulation, we only care about the concentration c and investigate the time evolution of oil. We adopt the IMEX-LDG(1,2) scheme to compute the concentration at T = 0.1, 0.2, 0.5, 1.0, with N = N x = N y = 50 and τ = 0.1h. The distribution of oil at different time are shown in Fig. 1a–d, respectively. From these pictures we can observe the water is injected from the lower-left corner and the oil is pushed toward production well with time evolution.

7 Concluding Remarks The LDG spatial discretization coupled with two multi-step IMEX schemes are applied to flow and transport equations in two dimensional incompressible miscible displacement

123

J Sci Comput

problem. In the presented methods, explicit and implicit time discretization are used for the convection and diffusion part in the transport equation, respectively. To decouple the numerical computation of the flow equation from the transport equation, and to avoid the solving of a nonlinear algebraic system, we use the extrapolation linearization for the diffusion–dispersion tensor in transport equation. We also proved in this paper the optimal convergence order in L ∞ (L 2 ) norm for both the concentration and the velocity, under weak temporal-spatial conditions. In the future work, we will design more efficient and higher order time-marching methods for this problem.

References 1. Bartels, S., Jensen, M., Müller, R.: Discontinuous Galerkin finite element convergence for incompressible miscible displacement problem of low regularity. SIAM J. Numer. Anal. 47, 3720–3743 (2009) 2. Bassi, F., Rebay, S.: A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations. J. Comput. Phys. 131, 267–279 (1997) 3. Bear, J., Bachmat, Y.: A generalized theory of hydrodynamic dispersion in porous media. Symposium of Haifa. Int. Assoc. Sci. Hydrol. 72, 7–16 (1967) 4. Bear, J., Bachmat, Y.: Introduction to Modeling of Transport Phenomena in Porous Media. Springer, New York (1990) 5. Castillo, P., Cockburn, B., Schötzau, D., Schwab, C.: Optimal a priori error estimates for the hp-version of the LDG method for convection diffusion problems. Math. Comput. 71, 455–478 (2002) 6. Ciarlet, P.: The Finite Element Method for Elliptic Problem. North Holland, Englewood Cliffs (1975) 7. Cockburn, B., Kanschat, G., Perugia, I., Schötzau, D.: Superconvergence of the local discontinuous Galerkin method for elliptic problems on cartesian grids. SIAM J. Numer. Anal. 39, 264–285 (2001) 8. Cockburn, B., Shu, C.-W.: The local discontinuous Galerkin method for time-dependent convection– diffusion systems. SIAM J. Numer. Anal. 35, 2440–2463 (1998) 9. Dong, B., Shu, C.-W.: Analysis of a local discontinuous Galerkin methods for linear time-dependent fourth-order problems. SIAM J. Numer. Anal. 47, 3240–3268 (2009) 10. Douglas Jr., J., Ewing, R.E., Wheeler, M.F.: A time-discretization procedure for a mixed finite element approximation of miscible displacement in porous media. R.A.I.R.O. Analyse numérique 17, 249–256 (1983) 11. Durán, R.G.: On the approximation of miscible displacement in porous media by a method of characteristics combined with a mixed method. SIAM J. Numer. Anal. 25, 989–1001 (1988) 12. Ewing, R.E., Russell, T.F., Wheeler, M.F.: Convergence analysis of an approximation of miscible displacement in porous media by mixed finite elements and a modified method of characteristics. Comput. Methods Appl. Mech. Eng. 47, 73–92 (1984) 13. Ewing, R.E., Wheeler, M.F.: Galerkin methods for miscible displacement problems in porous media. SIAM J. Numer. Anal. 17, 351–365 (1980) 14. Gottlieb, S., Wang, C.: Stability and convergence analysis of fully discrete Fourier collocation spectral method for 3-D viscous Burgers’ equation. J. Sci. Comput. 53, 102–128 (2012) 15. Guo, H., Yu, F., Yang, Y.: Local discontinuous Galerkin method for incompressible miscible displacement problem in porous media. J. Sci. Comput. 71, 615–633 (2017) 16. Guo, H., Zhang, Q., Yang, Y.: A combined mixed finite element method and local discontinuous Galerkin method for miscible displacement problem in porous media. Sci. China Math. 57, 2301–2320 (2014) 17. Jaffre, J., Roberts, J.E.: Upstream weighting and mixed finite elements in the simulation of miscible displacements. ESAIM: M2AN 19, 443–460 (1985) 18. Li, B., Sun, W.: Unconditional convergence and optimal error estimates of a Galerkin-mixed FEM for incompressible miscible flow in porous media. SIAM J. Numer. Anal. 51, 1959–1977 (2013) 19. Li, B., Wang, J., Sun, W.: The stability and convergence of fully Galerkin–Galerkin FEMs for porous medium flows. Commun. Comput. Phys. 15, 1141–1158 (2014) 20. Rivière, B.: Discontinuous Galerkin Finite Element Methods for Solving the Miscible Displacement Problem in Porous Media. Ph.D. Thesis, The University of Texas at Austin (2000) 21. Shu, C.-W.: Discontinuous Galerkin methods: general approach and stability. Numerical solutions of partial differential equations. In: Bertoluzza, S., Falletta, S., Russo, G., Shu, C.-W. (eds.) Advanced Courses in Mathematics CRM Barcelona, pp. 149–201. Basel, Birkhäuser (2009)

123

J Sci Comput 22. Sun, S., Rivière, B., Wheeler, M.F.: A combined mixed finite element and discontinuous Galerkin method for miscible displacement problem in porous media. In: Chan, T., et al. (eds.) Recent Progress in Computational and Applied PDEs, pp. 323–351. Kluwer Academic Publishers, Dordrecht (2002) 23. Wang, H., Liang, D., Ewing, R.E., Lyons, S.L., Qin, G.: An approximation to miscible fluid flows in porous media with point sources and sinks by an Eulerian–Lagrangian localized adjoint method and mixed finite element methods. SIAM J. Sci. Comput. 22, 561–581 (2000) 24. Wang, H., Shu, C.-W., Zhang, Q.: Stability and error estimates of local discontinuous Galerkin methods with implicit–explicit time-marching for advection–diffusion problems. SIAM J. Numer. Anal. 53, 206– 227 (2015) 25. Wang, H., Shu, C.-W., Zhang, Q.: Stability analysis and error estimates of local discontinuous Galerkin methods with implicit–explicit time-marching for nonlinear convection–diffusion problems. Appl. Math. Comput. 272, 237–258 (2016) 26. Wang, H., Wang, S., Zhang, Q., Shu, C.-W.: Local discontinuous Galerkin methods with implicit–explicit time-marching for multi-dimensional convection–diffusion problems. ESAIM: M2AN 50, 1083–1105 (2016) 27. Wheeler, M.F., Darlow, B.L.: Interiori Penalty Galerkin Methods for Miscible Displacement Problems in Porous Media. Computational Methods in Nonlinear Mechanics, pp. 458–506. North-Holland, Amsterdam (1980) 28. Xu, Y., Shu, C.-W.: Local discontinuous Galerkin methods for nonlinear Schrödinger equations. J. Comput. Phys. 205, 72–97 (2005) 29. Xu, Y., Shu, C.-W.: Local discontinuous Galerkin methods for high-order time-dependent partial differential equations. Commun. Comput. Phys. 7, 1–46 (2010) 30. Yan, J., Shu, C.-W.: A local discontinuous Galerkin method for KdV type equations. SIAM J. Numer. Anal. 40, 769–791 (2002) 31. Yan, J., Shu, C.-W.: Local discontinuous Galerkin methods for partial differential equations with higher order derivatives. J. Sci. Comput. 17, 27–47 (2002) 32. Zhang, Q., Shu, C.-W.: Error estimates to smooth solutions of Runge–Kutta discontinuous Galerkin methods for scalar conservation laws. SIAM J. Numer. Anal. 42, 641–666 (2004) 33. Zhang, Q., Shu, C.-W.: Stability analysis and a priori error estimates of the third order explicit Runge– Kutta discontinuous Galerkin method for scalar conservation laws. SIAM J. Numer. Anal. 48, 1038–1063 (2010)

123