Abstract

An efficient combination method of Laplace transform and mixed multiscale finite-element method for coupling partial differential equations of flow in a dual-permeability system is present. First, the time terms of parabolic equation with unknown pressure term are removed by the Laplace transform. Then the transformed equations are solved by mixed FEMs which can provide the numerical approximation formulas for pressure and velocity at the same time. With some assumptions, the multiscale basis functions are constructed by utilizing the effects of fine-scale heterogeneities through basis functions formulation computed from local flow problems. Without time step in discrete process, the present method is efficient when solving spatial discrete problems. At last, the associated pressure transform is inverted by the method of numerical inversion of the Laplace transform.

1. Introduction

In recent years, there have been some important developments in numerical methods of flow in fractured porous media. And an important numerical method is mixed finite-element method. Due to the heterogeneous properties and the multiscale nature of media, a complete analysis of these problems is extremely difficult, and efficient numerical solvers usually require an extremely large amount of computer storage and CPU time. Developing multiscale methods allow us to overcome above difficulty while retaining a satisfactory accuracy. The popular multiscale methods include multiscale finite-element method, numerical upscaling method, and multiscale finite volume method [1, 2]. In general, multiscale finite-element methods are regarded as numerical methods and strategies in which basis functions are computed by solving local homogeneous PDEs subject to special boundary conditions [1, 3–11].

Multiscale FEMs are efficiency and convenience for elliptic equation of steady flow model. Despite their advantage, they still have some practical limitations in solving parabolic equations for nonsteady flow model [12]. First, their major drawback is that it is necessary to take small time steps. A severe limitation on the time step may require an excessive amount of computer time. Second, all the interior velocity or pressure must be computed at each time step. We can see that numerical methods for parabolic equation with general FEMs require the solution of some simultaneous algebraic equations at each time step. On the other hand, these approaches tend to increase cost when the solutions must be carried out over long time periods [9].

The present method, the combined use of the Laplace transform and the multiscale finite-element method, is used to solve some problems of flow in fractured porous media. It is efficient to overcome the above difficulty.

2. Mathematical Model and Laplace Transform

Incompressible flow in fractured porous media can be described by some coupled partial differential equations. We consider the following mathematical model consisting of four coupling equations for pressure 𝑝 1 , 𝑝 2 and Darcy velocity 𝑢 1 , 𝑢 2 in a spatial domain [13, 14]: Φ 1 𝐶 1 𝜕 𝑝 1 ( 𝑥 , 𝑡 ) 𝑢 𝜕 𝑡 − d i v 1 𝑝 ( 𝑥 , 𝑡 ) − 𝜆 𝛼 2 ( 𝑥 , 𝑡 ) − 𝑝 1 Φ ( 𝑥 , 𝑡 ) = 0 , 𝑥 ∈ Ω , 𝑡 > 0 , ( 2 . 1 ) 2 𝐶 2 𝜕 𝑝 2 ( 𝑥 , 𝑡 ) 𝑢 𝜕 𝑡 − d i v 2 𝑝 ( 𝑥 , 𝑡 ) + 𝜆 𝛼 2 ( 𝑥 , 𝑡 ) − 𝑝 1 𝑢 ( 𝑥 , 𝑡 ) = 0 , 𝑥 ∈ Ω , 𝑡 > 0 , ( 2 . 2 ) 1 ( 𝑥 , 𝑡 ) = − 𝜆 𝐾 1 ( 𝑥 ) ∇ 𝑝 1 ( 𝑢 𝑥 , 𝑡 ) , ( 2 . 3 ) 2 ( 𝑥 , 𝑡 ) = − 𝜆 𝐾 2 ( 𝑥 ) ∇ 𝑝 2 𝑝 ( 𝑥 , 𝑡 ) , ( 2 . 4 ) 1 ( 𝑥 , 0 ) = 𝑝 0 1 , 𝑝 ( 2 . 5 ) 2 ( 𝑥 , 0 ) = 𝑝 0 2 , ( 2 . 6 ) where Ω is a domain in 𝑅 𝑑 ( 𝑑 = 2 ) ,    𝐶 𝑖 ( 𝑖 = 1 , 2 ) denote the compressibility, 𝜙 𝑖 ( 𝑖 = 1 , 2 ) are the porosity, and 𝑘 𝑖 ( 𝑥 ) ( 𝑖 = 1 , 2 ) denote the heterogeneous field permeability tensor which have scale separation and periodicity. We assume that only 𝑢 𝑖 , 𝑝 𝑖 ( 𝑖 = 1 , 2 ) are unknown expressions. And here we assume that the equations are equipped with Neumann boundary conditions → 𝑢 1 ⋅ → 𝑛 = 0 o n 𝜕 Ω , ( 2 . 7 ) → 𝑢 2 ⋅ → 𝑛 = 0 o n 𝜕 Ω , ( 2 . 8 ) where → 𝑛 is the outward normal of 𝜕 Ω .

In order to remove time dependence from the governing equation and boundary conditions, the scheme of the Laplace transform will be utilized. The Laplace transform of the function 𝑝 and its inversion formulas are defined as [ ] = 𝑃 ( 𝑥 , 𝑠 ) = 𝐿 𝑝 ( 𝑥 , 𝑡 ) ∞ 0 𝑒 − 𝑠 𝑡 𝑝 ( 𝑥 , 𝑡 ) 𝑑 𝑡 , 𝑝 ( 𝑥 , 𝑡 ) = 𝐿 − 1 = 1 𝑝 ( 𝑥 , 𝑠 ) 2 𝜋 𝑖 𝑐 + 𝑖 ∞ 𝑐 − 𝑖 ∞ 𝑒 𝑠 𝑡 𝑝 ( 𝑥 , 𝑠 ) 𝑑 𝑠 , ( 2 . 9 ) where 𝑠 = 𝑐 + 𝑖 𝜔 , 𝑐 , 𝜔 ∈ 𝑅 .

The Laplace transform of (2.1)-(2.2) is: Φ 1 𝐶 1 𝑠 𝑝 ( 𝑥 , 𝑠 ) − 𝑝 1 ( 𝑥 , 0 ) − d i v 𝑢 1 ( 𝑥 , 𝑠 ) − 𝜆 𝛼 𝑝 2 − 𝑝 1 Φ = 0 , ( 2 . 1 0 ) 2 𝐶 2 𝑠 𝑝 2 ( 𝑥 , 𝑠 ) − 𝑝 2 ( 𝑥 , 0 ) − d i v 𝑢 2 ( 𝑥 , 𝑠 ) + 𝜆 𝛼 𝑝 2 − 𝑝 1 = 0 . ( 2 . 1 1 )

The expression (2.10)-(2.11) can be written as Φ 1 𝐶 1 𝑠 + 𝜆 𝛼 𝑝 1 − 𝜆 𝛼 𝑝 2 − Φ 1 𝐶 1 𝑝 0 1 = d i v 𝑢 1 , Φ 2 𝐶 2 𝑠 + 𝜆 𝛼 𝑝 2 − 𝜆 𝛼 𝑝 1 − Φ 2 𝐶 2 𝑝 0 2 = d i v 𝑢 2 . ( 2 . 1 2 )

Now the unknown function 𝑢 𝑖 ( 𝑥 , 𝑠 ) and 𝑝 𝑖 ( 𝑥 , 𝑠 ) , 𝑖 = 1 , 2 , will be solved.

By the Laplace transform of (2.3)-(2.4), we can get the following: 𝑢 1 ( 𝑥 , 𝑡 ) = − 𝜆 𝐾 1 ( 𝑥 ) ∇ 𝑝 1 ( 𝑥 , 𝑡 ) , ( 2 . 1 3 ) 𝑢 2 ( 𝑥 , 𝑡 ) = − 𝜆 𝐾 2 ( 𝑥 ) ∇ 𝑝 2 ( 𝑥 , 𝑡 ) . ( 2 . 1 4 )

By the Laplace transform of Neumann boundary condition (2.7)-(2.8), we can know → 𝑢 𝑖 ⋅ → 𝑛 = 0 , o n 𝜕 Ω , 𝑖 = 1 , 2 . ( 2 . 1 5 )

Next we consider mixed multiscale finite-element method to (2.10)–(2.15).

3. Mixed Finite-Element Method

The mixed finite-element method is based on a mixed formulation [6–8, 15]. For our model problem (2.10)–(2.15), the mixed formulation is as follows.

Find ( 𝑢 𝑖 , 𝑝 𝑖 ) ( 𝑖 = 1 , 2 ) ∈ 𝐻 0 1 , d i v ( Ω ) × 𝐿 2 ( Ω ) , such that Ω 𝑘 𝑣 ⋅ 𝑖 𝜆 − 1 𝑢 𝑖 𝑑 𝑥 + Ω 𝑣 ⋅ ∇ 𝑝 𝑖 𝑑 𝑥 = 0 , ( 3 . 1 ) Ω d i v 𝑢 𝑖 Φ 𝑤 𝑑 𝑥 − 𝑖 𝐶 𝑖 𝑠 + 𝜆 𝛼 𝑝 𝑖 − 𝜆 𝛼 𝑝 3 − 𝑖 − Φ 𝑖 𝐶 𝑖 𝑝 0 𝑖 𝑤 𝑑 𝑥 = 0 , ( 3 . 2 ) which hold for all 𝑣 ∈ 𝐻 0 1 , d i v ( Ω ) and 𝑤 ∈ 𝐿 2 ( Ω ) , 𝐻 0 1 , d i v 𝐿 ( Ω ) = 𝑣 ∈ 2 ( Ω ) 𝑑 ∶ d i v ( 𝑣 ) ∈ 𝐿 2 ( Ω ) , 𝑣 ⋅ 𝑛 = 0 o n 𝜕 Ω , ( 3 . 3 ) where 𝐿 2 ( Ω ) is the space of square integral function in Ω . The above formulas (3.1)-(3.2) can provides the numerical results of 𝑝 𝑖 ( 𝑥 , 𝑠 ) and 𝑢 𝑖 ( 𝑥 , 𝑠 ) ( 𝑖 = 1 , 2 ) at the same time.

In mixed finite-element method, 𝐿 2 ( Ω ) and 𝐻 0 1 , d i v ( Ω ) are approximated by finite-dimensional subspaces 𝑊   and 𝑉 , respectively [1, 6]. For instance, 𝐿 2 ( Ω ) is replaced by 𝑊 : 𝑊 = 𝑤 ∈ 𝐿 2 ( Ω ) ∶ 𝑤 | Ω 𝑘 = c o n s t ∀ Ω 𝑘 . ∈ Ω ( 3 . 4 )

And 𝐻 0 1 , d i v ( Ω ) is replaced by 𝑉 : 𝑉 = 𝑣 ∈ 𝐻 0 1 , d i v ( Ω ) ∶ d i v ( 𝑣 ) ∈ 𝐿 2 ( Ω ) , ∀ Ω 𝑖 ∈ Ω , 𝑣 ⋅ 𝑛 𝑖 𝑗 | 𝛾 𝑖 𝑗 = c o n s t , ∀ 𝛾 𝑖 𝑗 ∈ Ω , a n d 𝑣 ⋅ 𝑛 𝑖 𝑗 i s c o n t i n u o u s a r o s s 𝛾 𝑖 𝑗 . ( 3 . 5 ) Here 𝑛 𝑖 𝑗 is the unit normal to 𝛾 𝑖 𝑗 pointing from Ω 𝑖 to Ω 𝑗 .

Denoting the approximation of ( 𝑝 𝑖 , 𝑢 𝑖 ) ( 𝑖 = 1 , 2 ) individually by ( 𝑝 𝑖 ℎ , 𝑢 𝑖 ℎ ) ( 𝑖 = 1 , 2 ) , the discrete formulation reads as follows.

Find ( 𝑝 𝑖 ℎ , 𝑢 𝑖 ℎ ) ( 𝑖 = 1 , 2 ) ∈ 𝑊 × 𝑉 , such that (3.1)-(3.2) hold for all 𝑤 ∈ 𝑊 and 𝑣 ∈ 𝑉 .

By applying Green formulation for (3.1), we can conclude the following: Ω 𝑘 𝑣 ⋅ 1 𝜆 − 1 𝑢 ℎ 1 𝑑 𝑥 − Ω 𝑝 ℎ 1 ∇ ⋅ 𝑣 𝑑 𝑥 = 0 , ( 3 . 6 ) Ω 𝑘 𝑣 ⋅ 2 𝜆 − 1 𝑢 ℎ 2 𝑑 𝑥 − Ω 𝑝 ℎ 2 ∇ ⋅ 𝑣 𝑑 𝑥 = 0 . ( 3 . 7 ) Now letting 𝑉 = s p a n { 𝜓 𝑗 } and 𝑤 = s p a n { 𝜑 𝑗 } , we let 𝑣 be equal to velocity basis function 𝜓 𝑘 and 𝑤 equal pressure basis function 𝜑 𝑘 ( 𝑘 = 1 , 2 , … , 𝑁 .). And the approximations of 𝑢 1 ℎ , 𝑢 2 ℎ and 𝑝 1 ℎ , 𝑝 2 ℎ can be written as follows: 𝑢 ℎ 1 = 𝑁 𝑗 = 1 𝑈 1 𝑗 ( 𝑠 ) ⋅ 𝜓 𝑗 ( 𝑥 ) , 𝑢 ℎ 2 = 𝑁 𝑗 = 1 𝑈 𝑗 2 ( 𝑠 ) ⋅ 𝜓 𝑗 ( 𝑥 ) , 𝑝 ℎ 1 = 𝑁 𝑗 = 1 𝑃 𝑗 1 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) , 𝑝 ℎ 2 = 𝑁 𝑗 = 1 𝑃 𝑗 2 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) . ( 3 . 8 ) Then formulas (3.2)–(3.7) can show that Ω 𝜓 𝑘 ⋅ 𝑘 1 𝜆 − 1 𝑁 𝑗 = 1 𝑈 1 𝑗 ( 𝑠 ) ⋅ 𝜓 𝑗 ( 𝑥 ) 𝑑 𝑥 − Ω 𝑁 𝑗 = 1 𝑃 1 𝑗 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) ∇ ⋅ 𝜓 𝑘 𝑑 𝑥 = 0 , Ω 𝜓 𝑘 ⋅ 𝑘 2 𝜆 − 1 𝑁 𝑗 = 1 𝑈 2 𝑗 ( 𝑠 ) ⋅ 𝜓 𝑗 ( 𝑥 ) 𝑑 𝑥 − Ω 𝑁 𝑗 = 1 𝑃 2 𝑗 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) ∇ ⋅ 𝜓 𝑘 𝑑 𝑥 = 0 , Ω d i v 𝑁 𝑗 = 1 𝑈 1 𝑗 ( 𝑠 ) ⋅ 𝜓 𝑗 𝜑 ( 𝑥 ) 𝑘 = 𝑑 𝑥 Ω Φ 1 𝐶 1 𝑠 + 𝜆 𝛼 𝑁 𝑗 = 1 𝑃 𝑗 1 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) − 𝜆 𝛼 𝑁 𝑗 = 1 𝑃 𝑗 2 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) − Φ 1 𝐶 1 𝑝 0 1 𝜑 𝑘 𝑑 𝑥 , Ω d i v 𝑁 𝑗 = 1 𝑈 2 𝑗 ( 𝑠 ) ⋅ 𝜓 𝑗 𝜑 ( 𝑥 ) 𝑘 = 𝑑 𝑥 Ω Φ 2 𝐶 2 𝑠 + 𝜆 𝛼 𝑁 𝑗 = 1 𝑃 𝑗 2 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) − 𝜆 𝛼 𝑁 𝑗 = 1 𝑃 𝑗 1 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) − Φ 2 𝐶 2 𝑝 0 2 𝜑 𝑘 𝑑 𝑥 , 𝑘 = 1 , 2 , … , 𝑁 . ( 3 . 9 )

The above four formulas can be rearranged as follows: Ω 𝜓 𝑘 ⋅ 𝑘 1 𝜆 − 1 𝑁 𝑗 = 1 𝑈 𝑗 1 ( 𝑠 ) ⋅ 𝜓 𝑗 ( 𝑥 ) 𝑑 𝑥 − Ω 𝑁 𝑗 = 1 𝑃 𝑗 1 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) ∇ ⋅ 𝜓 𝑘 𝑑 𝑥 = 0 , Ω d i v 𝑁 𝑗 = 1 𝑈 𝑗 1 ( 𝑠 ) ⋅ 𝜓 𝑗 𝜑 ( 𝑥 ) 𝑘 𝑑 𝑥 − Ω Φ 1 𝐶 1 𝑠 + 𝜆 𝛼 𝑁 𝑗 = 1 𝑃 𝑗 1 ( 𝑠 ) ⋅ 𝜑 𝑗 𝜑 ( 𝑥 ) 𝑘 + 𝑑 𝑥 Ω 𝜆 𝛼 𝑁 𝑗 = 1 𝑃 𝑗 2 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) 𝑑 𝑥 = − Ω Φ 1 𝐶 1 𝑝 0 1 𝜑 𝑘 𝑑 𝑥 , Ω 𝜆 𝛼 𝑁 𝑗 = 1 𝑃 𝑗 1 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) 𝑑 𝑥 − Ω Φ 2 𝐶 2 𝑠 + 𝜆 𝛼 𝑁 𝑗 = 1 𝑃 𝑗 2 ( 𝑠 ) ⋅ 𝜑 𝑗 𝜑 ( 𝑥 ) 𝑘 + 𝑑 𝑥 Ω d i v 𝑁 𝑗 = 1 𝑈 𝑗 2 ( 𝑠 ) ⋅ 𝜓 𝑗 𝜑 ( 𝑥 ) 𝑘 𝑑 𝑥 = − Ω Φ 2 𝐶 2 𝑝 0 2 𝜑 𝑘 𝑑 𝑥 − Ω 𝑁 𝑗 = 1 𝑃 𝑗 2 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) ∇ ⋅ 𝜓 𝑘 + 𝑑 𝑥 Ω 𝜓 𝑘 ⋅ 𝑘 2 𝜆 − 1 𝑁 𝑗 = 1 𝑈 𝑗 2 ( 𝑠 ) ⋅ 𝜓 𝑗 ( 𝑥 ) 𝑑 𝑥 = 0 , 𝑘 = 1 , 2 , … , 𝑁 .

Then, we obtain the approximations: 𝑢 ℎ 1 = 𝑁 𝑗 = 1 𝑈 𝑗 1 ( 𝑠 ) ⋅ 𝜓 𝑗 ( 𝑥 ) , 𝑢 ℎ 2 = 𝑁 𝑗 = 1 𝑈 𝑗 2 ( 𝑠 ) ⋅ 𝜓 𝑗 ( 𝑥 ) , 𝑝 ℎ 1 = 𝑁 𝑗 = 1 𝑃 𝑗 1 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) , 𝑝 ℎ 2 = 𝑁 𝑗 = 1 𝑃 𝑗 2 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) , ( 3 . 1 1 ) where the coefficients 𝑈 𝑗 1 , 𝑃 𝑗 1 , 𝑃 𝑗 2 and 𝑈 𝑗 2 are obtained by solving the following linear system: ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐴 1 𝐵 1 𝐵 0 0 2 𝐷 1 𝐷 2 0 0 𝐷 3 𝐷 4 𝐵 3 0 0 𝐵 4 𝐴 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑈 𝑗 1 𝑃 𝑗 1 𝑃 𝑗 2 𝑈 𝑗 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 𝐹 1 𝐹 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 3 . 1 2 ) where 𝐴 1 = Ω 𝜓 𝑘 ⋅ 𝜆 𝐾 1 − 1 𝜓 𝑗 𝑑 𝑥 , 𝐴 2 = Ω 𝜓 𝑘 ⋅ 𝜆 𝐾 2 − 1 𝜓 𝑗 , 𝐵 𝑑 𝑥 1 = 𝐵 4 = − Ω 𝜑 𝑗 𝜓 d i v 𝑘 𝑑 𝑥 , 𝐵 2 = 𝐵 3 = Ω 𝜑 𝑘 𝜓 d i v 𝑗 , 𝐷 𝑑 𝑥 1 = − Ω Φ 1 𝐶 1 𝜑 s + 𝜆 𝛼 𝑗 𝜑 𝑘 𝑑 𝑥 , 𝐷 2 = 𝐷 3 = Ω 𝜆 𝛼 𝜑 𝑗 , 𝐷 𝑑 𝑥 4 = − Ω Φ 2 𝐶 2 𝜑 s + 𝜆 𝛼 𝑗 𝜑 𝑘 , 𝐹 𝑑 𝑥 1 = − Ω Φ 1 𝐶 1 𝑝 0 1 𝜑 𝑘 𝑑 𝑥 , 𝐹 2 = − Ω Φ 2 𝐶 2 𝑝 0 2 𝜑 𝑘 d x ( 𝑘 , 𝑗 = 1 , 2 , … , 𝑁 ) . ( 3 . 1 3 )

Then the linear system (3.12) can be written as ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐴 1 𝐵 1 𝐵 0 0 2 𝐷 1 𝐷 2 0 0 𝐷 2 𝐷 4 𝐵 2 0 0 𝐵 1 𝐴 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑈 𝑗 1 𝑃 𝑗 1 𝑃 𝑗 2 𝑈 𝑗 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 𝐹 1 𝐹 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 3 . 1 4 )

4. Multiscale Basis Functions

The main idea of the MMsFEM (Mixed Multiscale Finite-Element Method) is to construct special local basis functions that are adaptive to the local properties of the elliptic differential operator [6]. In this paper, two sets of grids are considered: a fine grid and a coarsened grid in which each coarse block Ω 𝐻 consists of connected cells from the underlying fine grid. Local basis functions 𝜑 𝐻 𝑖 for pressure and 𝜓 𝐻 𝑖 for velocity in each coarse block Ω 𝐻 generally are constructed by solving ∇ ⋅ 𝜓 𝑖 𝐻 = 𝑇 𝐻 ( 𝑥 ) , f o r 𝑥 ∈ Ω 𝐻 , 𝜓 0 , e l s e , 𝐻 𝑖 = − 𝜆 𝐾 ∇ 𝜑 𝐻 𝑖 , 𝜓 𝐻 𝑖 ⋅ 𝑛 = 0 , o n 𝜕 Ω 𝐻 , ( 4 . 1 ) where subscript 𝑖 denotes the order number of the basis function.

The weighting function 𝑇 𝐻 ( 𝑥 ) plays an important role to distribute d i v ( 𝑢 𝑖 ) ( 𝑖 = 1 , 2 ) onto the fine grid appropriately. It can be chosen on the form 𝑇 𝐻 ∫ ( 𝑥 ) = 𝑙 ( 𝑥 ) / ( Ω 𝐻 𝑙 ( 𝑥 ) 𝑑 𝑥 ) with various choices of 𝑙 ( 𝑥 ) [6].

Now denote 𝑉 𝑚 𝑠 = s p a n { 𝜓 𝐻 𝑖 } the approximation space of 𝑉 and 𝑊 = s p a n { 𝜑 𝐻 𝑗 } , we let 𝑣 𝐻 equal velocity basis function 𝜓 𝐻 𝑖 and 𝑤 𝐻 equal pressure basis function 𝜑 𝑖 𝐻 ( 𝑖 , 𝑗 = 1 , 2 , … 𝑁 .). And the approximations of 𝑢 𝐻 1 , 𝑢 𝐻 2 and 𝑝 𝐻 1 , 𝑝 𝐻 2 can be written as follows: 𝑢 𝐻 1 = 𝑁 ′ j = 1 𝑈 1 𝐻 , 𝑗 ( 𝑠 ) ⋅ 𝜓 𝑗 ( 𝑥 ) , 𝑢 𝐻 2 = 𝑁 ′ 𝑗 = 1 𝑈 2 𝐻 , 𝑗 ( 𝑠 ) ⋅ 𝜓 𝑗 ( 𝑥 ) , 𝑝 𝐻 1 = 𝑁 ′ 𝑗 = 1 𝑃 1 𝐻 , 𝑗 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) , 𝑝 𝐻 2 = 𝑁 ′ 𝑗 = 1 𝑃 2 𝐻 , 𝑗 ( 𝑠 ) ⋅ 𝜑 𝑗 ( 𝑥 ) . ( 4 . 2 )

Then the MMsFEM seeks ( 𝑝 𝐻 𝑖 , 𝑢 𝐻 𝑖 ) ( 𝑖 = 1 , 2 ) ∈ 𝑊 × 𝑉 𝑚 𝑠 ,such that (3.1)-(3.2) hold for all 𝑣 𝐻 ∈ 𝑉 𝑚 𝑠 and 𝑤 𝐻 ∈ 𝑊 . We can get the following formulas on a similar conclusion from (3.6)-(3.7) to expression (3.14): ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐴 1 𝐵 1 𝐵 0 0 2 𝐷 1 𝐷 2 0 0 𝐷 4 𝐷 4 𝐵 2 0 0 𝐵 1 𝐴 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑈 1 𝐻 , 𝑗 𝑃 1 𝐻 , 𝑗 𝑃 2 𝐻 , 𝑗 𝑈 2 𝐻 , 𝑗 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 𝐹 1 𝐹 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ 𝐴 , ( 4 . 3 ) 1 = Ω 𝜓 𝑘 𝐻 ⋅ 𝜆 𝐾 1 − 1 𝜓 𝐻 𝑗 𝑑 𝑥 , 𝐴 2 = Ω 𝜓 𝐻 𝑘 ⋅ 𝜆 𝐾 2 − 1 𝜓 𝐻 𝑗 , 𝐵 𝑑 𝑥 1 = 𝐵 4 = − Ω 𝜑 𝐻 𝑗 𝜓 d i v 𝐻 𝑘 𝑑 𝑥 , 𝐵 2 = 𝐵 3 = Ω 𝜑 𝐻 𝑘 𝜓 d i v 𝐻 𝑗 , 𝐷 𝑑 𝑥 1 = − Ω Φ 1 𝐶 1 𝜑 𝑠 + 𝜆 𝛼 𝐻 𝑗 𝜑 𝐻 𝑘 𝑑 𝑥 , 𝐷 2 = 𝐷 3 = Ω 𝜆 𝛼 𝜑 𝐻 𝑗 , 𝐷 𝑑 𝑥 4 = − Ω Φ 2 𝐶 2 𝜑 𝑠 + 𝜆 𝛼 𝐻 𝑗 𝜑 𝐻 𝑘 𝐹 𝑑 𝑥 1 = − Ω Φ 1 𝐶 1 𝑝 0 1 𝜑 𝐻 𝑘 , 𝐹 𝑑 𝑥 2 = − Ω Φ 2 𝐶 2 𝑝 0 2 𝜑 𝐻 𝑘 , 𝑑 𝑥 ( 𝑘 , 𝑗 = 1 , 2 , … , 𝑁 ) . ( 4 . 4 )

5. Numerical Inversion of the Pressure and Velocity Term

In this section we compute the numerical value of 𝑈 𝑖 𝐻 , 𝑗 ( 𝑡 ) ( 𝑖 = 1 , 2 ; 𝑗 = 1 , 2 , … , 𝑁 ) and 𝑃 𝑖 𝐻 , 𝑗 ( 𝑡 ) ( 𝑖 = 1 , 2 ; 𝑗 = 1 , 2 … , 𝑁 ) by 𝑈 𝑖 𝐻 , 𝑗 ( 𝑠 ) ( 𝑖 = 1 , 2 ; 𝑗 = 1 , 2 … 𝑁 ) and 𝑃 𝑖 𝐻 , 𝑗 ( 𝑠 ) ( 𝑖 = 1 , 2 ; 𝑗 = 1 , 2 … 𝑁 ) using the numerical inversion formula of Laplace transform. In accordance with the method of Durbin [16], the Fourier series expansion of 𝑈 𝑖 𝐻 , 𝑗 ( 𝑡 ) can be derived as follows [9, 17]: 𝑈 𝑖 𝐻 , 𝑗 1 ( 𝑡 ) = 𝑇 𝑒 𝑣 𝑡 − 1 2 𝑈 R e 𝑖 𝐻 , 𝑗 + ( 𝑣 ) ∞ ∑ 𝑘 = 0 𝑈 R e 𝑖 𝐻 , 𝑗 𝑆 𝑘 c o s 𝑘 𝜋 𝑡 𝑇 𝑈 − I m 𝑖 𝐻 , 𝑗 𝑆 𝑘 s i n 𝑘 𝜋 𝑡 𝑇 − 𝑈 0 𝑆 ( 𝑣 , 𝑡 , 𝑇 ) , 𝑘 = 𝑣 + 𝑖 𝑘 𝜋 𝑇 , 𝑘 = 0 , 1 , … , 𝑁 . ( 5 . 1 )

And here the discrete error is given by 𝑈 0 ( 𝑣 , 𝑡 , 𝑇 ) = ∞ 𝑘 = 1 𝑒 − 2 𝑣 𝑘 𝑇 𝑈 𝑖 𝐻 , 𝑗 ( 2 𝑘 𝑇 + 𝑡 ) . ( 5 . 2 ) So the numerical approximate value of 𝑈 𝑖 𝐻 , 𝑗 ( 𝑡 ) is 𝑈 𝐻 , 𝑗 𝑖 , 𝑁 1 ( 𝑡 ) = 𝑇 𝑒 𝑣 𝑡 − 1 2 𝑈 R e 𝑖 𝐻 , 𝑗 + 1 ( 𝑣 ) 𝑇 𝑒 𝑁 𝑣 𝑡 𝑘 = 0 𝑈 R e 𝑖 𝐻 , 𝑗 𝑆 𝑘 c o s 𝑘 𝜋 𝑡 𝑇 𝑈 − I m 𝑖 𝐻 , 𝑗 𝑆 𝑘 s i n 𝑘 𝜋 𝑡 𝑇 , ( 5 . 3 ) where 2 𝑇 is the period of Fourier series approximate inversion function on the interval [ 0 , 2 𝑇 ] and 𝑁 and 𝑣 * 𝑇 are free parameters. All 𝑈 𝑖 𝐻 , 𝑗 ( 𝑠 ) are computed by formulation (4.3).

Obviously, there occurs a truncation error given by 𝑈 𝐻 , 𝑗 𝑖 , 𝐴 1 ( 𝑡 ) = 𝑇 𝑒 𝑣 𝑡 ∞ 𝑘 = 𝑁 + 1 𝑈 R e 𝑖 𝐻 , 𝑗 𝑆 𝑘 c o s 𝑘 𝜋 𝑡 𝑇 𝑈 − I m 𝑖 𝐻 , 𝑗 𝑆 𝑘 s i n 𝑘 𝜋 𝑡 𝑇 . ( 5 . 4 )

Because numerical inversion formula requires the value of transformed function at each 𝑆 𝑘 ( 𝑘 = 0 , 1 , 2 , 3 , … , 𝑁 ) , it is needed to solve the equations for any 𝑆 = 𝑆 𝑘 .

The approximate value of 𝑃 𝑖 𝐻 , 𝑗 ( 𝑡 ) ( 𝑖 = 1 , 2 ; 𝑗 = 1 , 2 , … , 𝑁 ) can be computed in the same way. 𝑃 𝐻 , 𝑗 𝑖 , 𝑁 1 ( 𝑡 ) = 𝑇 𝑒 𝑣 𝑡 − 1 2 𝑃 R e 𝑖 𝐻 , 𝑗 + 1 ( 𝑣 ) 𝑇 𝑒 𝑁 𝑣 𝑡 𝑘 = 0 𝑃 R e 𝑖 𝐻 , 𝑗 𝑆 𝑘 c o s 𝑘 𝜋 𝑡 𝑇 𝑃 − I m 𝑖 𝐻 , 𝑗 𝑆 𝑘 s i n 𝑘 𝜋 𝑡 𝑇 . ( 5 . 5 )

The solution at a specific node in the given domain can be obtained from (4.3) and (5.3)–(5.5). This strategy is different from that of classical FEM, which must compute all nodal values for each time step until the specific time is reached. The present method can compute the specific nodal value at a specific time. In addition, it is obvious that the present method takes less computer time if lengthy solutions are required.

6. Conclusion

This paper focuses on efficient numerical method which combines the numerical inversion of Laplace transform with mixed FEM for coupling partial differential equations of flow in a dual-permeability system. The present combined method can save computing cost by multiscale method in spatial domain, and overcome some disadvantages coming from numerical processing of different time step size by the numerical inversion of Laplace transform. However, the present method can be extended to other two- or three-dimensional linear time-dependent problems.

Acknowledgments

This research was supported by Natural Science-Technology Major Special Task of China (2011ZX05008-004-44), Natural Science Foundation of China (11161002), and Natural Science Foundation of Jiangxi province (20114BAB201016).