Next Article in Journal / Special Issue
Molecular Simulation of Shale Gas Adsorption and Diffusion in Clay Nanopores
Previous Article in Journal
Optical Properties of Silicon-Rich Silicon Nitride (SixNyHz) from First Principles
Previous Article in Special Issue
An Incompressible, Depth-Averaged Lattice Boltzmann Method for Liquid Flow in Microfluidic Devices with Variable Aperture
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multiscale Simulations for Coupled Flow and Transport Using the Generalized Multiscale Finite Element Method

1
Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong Kong, China
2
Center for Numerical Porous Media (NumPor), King Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia
3
Department of Mathematics, Texas A&M University, College Station, TX 77843, USA
*
Author to whom correspondence should be addressed.
Submission received: 9 October 2015 / Revised: 19 November 2015 / Accepted: 1 December 2015 / Published: 11 December 2015
(This article belongs to the Special Issue Advances in Modeling Flow and Transport in Porous Media)

Abstract

:
In this paper, we develop a mass conservative multiscale method for coupled flow and transport in heterogeneous porous media. We consider a coupled system consisting of a convection-dominated transport equation and a flow equation. We construct a coarse grid solver based on the Generalized Multiscale Finite Element Method (GMsFEM) for a coupled system. In particular, multiscale basis functions are constructed based on some snapshot spaces for the pressure and the concentration equations and some local spectral decompositions in the snapshot spaces. The resulting approach uses a few multiscale basis functions in each coarse block (for both the pressure and the concentration) to solve the coupled system. We use the mixed framework, which allows mass conservation. Our main contributions are: (1) the development of a mass conservative GMsFEM for the coupled flow and transport; (2) the development of a robust multiscale method for convection-dominated transport problems by choosing appropriate test and trial spaces within Petrov-Galerkin mixed formulation. We present numerical results and consider several heterogeneous permeability fields. Our numerical results show that with only a few basis functions per coarse block, we can achieve a good approximation.

1. Introduction

Many porous media problems occur over multiple scales. Because of scale disparity, some type of model reduction or averaging are often employed. Typical flow-based upscaling approaches ([1,2,3]) compute the permeabilities on a coarse (computational) grid and solve the flow and transport equations together on a coarse grid. These approaches are easy to implement with a minimal code modification; however, these approaches can result in large errors for complex heterogeneities because limited information is used in upscaling. In this paper, we consider a multiscale method for a coupled flow and transport system, where the flow equation for the pressure field is described by a steady state elliptic equation (derived by Darcy law) and the transport equation for the concentration field is described by a convection-dominated parabolic equation.
Our approach derives its foundation from rigorous multiscale techniques developed recently [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]. In our approach, multiscale basis functions are computed and used to solve flow equations. The main contributions of this paper are (1) to use multiscale basis functions for both flow and transport equations and (2) design of a novel mixed multiscale methods for convection-dominated transport equations within the context of Generalized Multiscale Finite Element Method (GMsFEM). The GMsFEM is introduced in [16], where the objective is to construct multiscale basis functions in each coarse domain. The main distinction of this approach from some of existing multiscale methods is that the GMsFEM systematically constructs multiscale basis functions using local spectral problems and the local snapshot spaces. Both of these concepts (local spectral problem and local snapshot spaces) are important in a design of the method.
In our proposed approach, we use a mixed formulation for both the flow and the transport equations. The use of a mixed formulation for the flow equation guarantees the mass conservation [7,26,27,28,29,30]. The use of a mixed formulation is important for the transport equation, guarantees the mass conservation and helps with the stabilization of it. In particular, by adding more flux basis functions, we can achieve a better stability based on our numerical studies [31]. In the mixed formulation for a coupled system, we first define snapshot spaces for the appropriately defined fluxes in the flow and the concentration equations. For the pressure and the concentration fields, we use piecewise constant basis functions. The snapshot space represents the solution space in each coarse region and provides a mass conservative approach. In the snapshot space, we perform a local spectral decomposition to identify multiscale basis functions. The functions corresponding to dominant modes are used to construct coarse spaces for the fluxes. In particular, we use only a few degrees of freedom (maximum two) in each coarse region. For the test functions for the convection-dominated transport, we use the solutions of adjoint problem. As a result, we have a full coarse-grid approximation for the coupled flow and transport system.
We would like to discuss a relation between our proposed approaches and the flow-based upscaling [2,3]. In the latter, the upscaled permeabilities are computed and, then, the coupled flow and transport equation is solved on a coarse grid. The computations of upscaled permeabilities can be compared to the computations of multiscale basis functions, which are performed locally in each coarse region by solving local problems. The resulting coarse-grid problem in the upscaled model is similar to our resulting coarse-grid system. There are several main advantages of the proposed method: (1) it can use more multiscale basis functions in each region and achieve a desired accuracy; (2) the number of multiscale basis functions can be adaptively selected in space and time; (3) our approach can easily reconstruct the fine-scale distributions of the solution; (4) our approaches can handle any coarse-grid geometries.
We consider several numerical examples by considering two different hydraulic conductivity fields. In our numerical examples, we use a relatively coarse grid with the coarsening factor of 20. We study the numerical errors in the concentration field and depict concentration profiles at different time instants. Our numerical results show that the use of only one basis function per coarse edge (which is similar to numerical upscaling of hydraulic conductivity) does not produce an accurate solution, while if we use two basis functions, we can obtain a much more accurate approximation of the solution. These preliminary numerical results show that our fully coarse-grid model for the coupled flow and transport equations can accurately predict the solution using only a few degrees of freedom in each coarse region.
Novelty. The novelty of our work can be summarized as follows:
  • The development of Generalized Multiscale Finite Element Method for a mixed formulation of the convection-diffusion equation using different test and trial spaces.
  • The development and implementation of a mass conservative mixed GMsFEM for the coupled flow and transport equations.
The paper is organized as follows. In the next section, we present preliminaries and describe coarse and fine grids as well as a mixed formulation. Section 3 is devoted to the construction of multiscale basis functions. Numerical results are presented in Section 4. Finally, we present conclusions.

2. Preliminaries and Notations

We consider a coupled flow and transport system that arises in many hydrological applications
c t + v · c - D Δ c = f c in Ω × ( 0 , T ) , v + κ p = 0 in Ω , · v = f q in Ω .
Here Ω is the computational domain, T > 0 is a fixed time, κ is a hydraulic conductivity, c is the concentration, p is the pressure, D is a diffusivity of the medium, and f c , f q are source terms. The equations are subject to the boundary conditions c ( x , t ) = 0 on Ω × [ 0 , T ] , p ( x ) = g on Ω and the initial condition c ( x , 0 ) = 0 in Ω. In order to have a mass conservation, we re-write the system in a mixed formulation.
c t + · q = f c + c f q in Ω × ( 0 , T ) , q + D c - v c = 0 in Ω × ( 0 , T ) , v + κ p = 0 in Ω , · v = f q in Ω .
Here q is an auxiliary variable, we call it flux or flux-concentration. It is well known that the Galerkin method inherits the stability of the continuous problem, and it yields to spurious oscillation when the convection coefficient is larger than the diffusive one ( v D ). The mixed finite element methods can be used to achieve a mass conservation. In this paper, we aim to construct multiscale basis functions for the flux q and the velocity v. As for the concentration c and the pressure p, we choose piecewise constant functions on a coarse grid.
Next, we introduce the discretized form of the Problem (2). We first comment on some notations of the coarse and the fine mesh. We assume that Ω is partitioned in a usual way into a union of rectangles (in 2-D) or hexadedrals (in 3-D) denoted by T H , where H is the coarse mesh size. The fine mesh is obtained by splitting each coarse block in T H into a connected union of fine-grid blocks, and we denote it by T h with mesh size h, h < H . In addition, we denote by E H the set of all edges/faces of elements in T H , i.e., E H = { E i : 1 i N e } , where N e is the number of coarse edges/faces. Note that each E i E H is the union of some fine-grid edges/faces. That is, E i = { e j : e j is on E i } , where e j are the fine-grid edges/faces lying on E i . We also define the coarse neighborhood of the edge/face E i by
ω i = { K T H : K E i ϕ } .
See Figure 1 for an illustration of a coarse neighborhood associated with one coarse edge, where the coarse mesh and the fine mesh are represented by black lines and gray lines, respectively.
Figure 1. An illustration of a coarse neighborhood ω i (in green) corresponding to the coarse edge E i .
Figure 1. An illustration of a coarse neighborhood ω i (in green) corresponding to the coarse edge E i .
Computation 03 00670 g001
Next, we define the coarse trial spaces for q , v , c , and p. The concentration c and the pressure p are approximated in the space of piecewise constant functions with respect to the coarse grid T H . We denote this space by Q H , and
Q H = { ψ : ψ | K i P 0 ( K i ) for all K i T H } .
As for the flux-concentration q and the velocity v, the basis functions are associated with coarse edges and are supported in the corresponding coarse neighborhoods. In particular, to obtain the basis functions for a coarse edge E i , we will solve a local problem in the coarse neighborhood ω i . Let Φ i q and Φ i v be the set of multiscale basis functions for q and v with respect to E i , respectively. Then, we define the multiscale solution space for q and v as
V H q = E i Φ i q and V H v = E i Φ i v .
Furthermore, we discretize the time interval [ 0 , T ] uniformly by the points
t n = n τ , n = 0 , 1 , , N
where τ is the time step and N = T / τ . Using the above spaces, the GMsFEM approximation of (2) can be given as follows. For all n 1 , find ( c H n , q H n , p H , v H ) Q H × V H q × Q H × V H v , such that
Ω c H n - c H n - 1 τ c ˜ + Ω ( · q H n - θ ) c ˜ = Ω ( f c n - θ + c H n - θ f q ) c ˜ for c ˜ Q H , Ω 1 D q H n - θ · q ˜ = Ω ( · q ˜ + v H D · q ˜ ) c H n - θ for q ˜ W H q , Ω 1 κ v H · v ˜ = Ω ( · v ˜ ) p H for v ˜ V H v , Ω ( · v H ) p ˜ = Ω f q p ˜ for p ˜ Q H ,
subject to the boundary condition p H = g . Here, we let z n - θ = θ z n - 1 + ( 1 - θ ) z n , where z = c H or q H . For θ = 0 , 1 / 2 and 1, we get the Backward Euler, Crank-Nicolson and Forward Euler Methods respectively. Note that W H q is the testing space for q H n . We emphasize that in general, W H q V H q , which allows a better stability. The construction of V H q , V H v , and W H q will be discussed in the next section.
In addition, to get the fine-scale solution, we use the standard lowest-order Raviart-Thomas space Q h × V h q × Q h × V h v for the approximation of (2) on the fine grid T h . The fine-scale solution ( c h , q h ) will be considered as a reference solution for the comparison with the multiscale solution ( c H , q H ) .

3. Generalized Multiscale Finite Element Method

In this section, we will discuss the construction of the vector solution spaces V H q , V H v , and the testing space W H q . Note that solving q H in (6) requires v H . Therefore, we will first compute v H using the last two equations of (6).
We briefly introduce the basic idea of the GMsFEM. In the framework of GMsFEM, it follows an offline-online procedure to generate the multiscale solution space. At the offline stage, we first construct the snapshot space obtained by solutions of local problems with all possible boundary conditions, which are resolvable by the fine-grid. Then, we choose a subspace of the snapshot space by selecting dominant modes via some local spectral decompositions. The resulting reduced dimensional space, called the offline space, will be used at the online stage to get a multiscale solution. Note that the snapshot space can be taken as the solution space to obtain a solution with a good accuracy. However, it is expensive to use the snapshot space due to its high dimensionality. Thus, there is a need to perform a space reduction within the snapshot space, and choose only few dominant modes to form the offline space. We remark that for parameter dependent problems (in which the coefficients depend on some extra parameters), one needs to construct an online space using the offline space and an appropriate spectral problem.

3.1. Multiscale Solution Space V H v

We form the multiscale solution space V H v for v H . In order to build V H v , we first construct the snapshot space, denoted by V snap v . For each coarse-grid edge E i E H , a local problem is defined on each of its neighboring coarse cells. Because the snapshot space is supposed to capture the multiscale information within this local domain, the corresponding boundary conditions of the local problem should cover all possibilities. Formally, the local problem is defined as the following: For each E i E H , on its neighboring coarse cell K = K r or K s (see Figure 2), find ( p j i , ψ v , j i ) Q h × V h v , such that
ψ v , j i + κ p j i = 0 in K , · ψ v , j i = α j i in K , ψ v , j i · n = 0 on K E i , ψ v , j i · n = δ j i on E i .
Here, j varies over all fine-grid edges on E i . For the convenience, we use the same notations as the solution to denote local snapshot fields. Here α j i is a constant defined on K, n is a fixed unit-normal vector on K , and δ j i ( x ) is defined on E i with respect to the fine-grid edges. Note that E i = { e j i : e j i in E i } . We require δ j i has value 1 on e j i and 0 on the other fine-grid edges, 1 j J i , where J i is the number of the fine-grid edges on E i . That is,
δ j i = 1 , on e j i , 0 , on the other fine-grid edges.
Figure 2. Neighboring coarse cells K r and K s corresponding to the coarse-grid edge E i . The fine-grid edge e j i and e j i are in red.
Figure 2. Neighboring coarse cells K r and K s corresponding to the coarse-grid edge E i . The fine-grid edge e j i and e j i are in red.
Computation 03 00670 g002
Note that the constant α j i in (7) is chosen so that the compatibility condition K α j i = ( ± ) E i δ j i is satisfied, for K = K r or K s . The above local Problem (7) can be solved numerically on the underlying fine grid of K by the lowest-order Raviart-Thomas element method.
After solving the local problems (on K r and K s ) with respect to E i , a local snapshot V snap v , i is generated by the solutions of (7). That is
V snap v , i = span { ψ v , j i : 1 j J i } .
Then, the snapshot space V snap v is formed as
V snap v = E i V snap v , i = span { ψ v , j i : 1 j J i , 1 i N e } .
We re-numerate the snapshot functions by a single index to create the snapshot matrix
R snap v = ψ 1 v , ψ 2 v , , ψ L snap v ,
where L snap = i = 1 N e J i is the dimension of the snapshot space V snap v . Note that the matrix R snap v maps from the coarse space V snap v to the fine space. The next step is to construct the offline space V off v for v H . For a coarse-grid edge E i , we define the following eigenvalue problem in V snap v , i :
A v off Ψ k v = λ k M v off Ψ k v ,
where
A v off = [ a s t ] = E i ( ψ v , s i · n ) ( ψ v , t i · n ) ,
and
M v off = [ m s t ] = ω i 1 κ ψ v , s i · ψ v , t i .
The eigenvalue Problem (9) will produce a set of eigenvectors for generating the local offline space V off v , i by first selecting L off i eigenvectors corresponding to the smallest L off i eigenvalues, and then computing the element ϕ v , k i = j = 1 J i Ψ k j v ψ v , j i for k = 1 , 2 , , L off i , where Ψ k j v is the j-th component of Ψ k v . Then, the local offline space with respect to E i is formed as
V off v , i = span { ϕ v , k i : 1 k L off i }
and the offline space V off v is formed as
V off v = E i V off v , i = span { ϕ v , k i : 1 k L off i , 1 i N e } .
We re-numerate the basis functions by a single index to create the offline matrix
R off v = ϕ 1 v , ϕ 2 v , , ϕ L off v ,
where L off = i = 1 N e L off i is the dimension of the offline space V off v .
We take the multiscale solution space as
V H v = V off v .
Then, v H is approximated by solving the problem: Find ( p H , v H ) Q H × V H v such that
Ω 1 κ v H · v ˜ = Ω ( · v ˜ ) p H for v ˜ V H v , Ω ( · v H ) q ˜ = Ω f q q ˜ for q ˜ Q H ,
subject to the boundary condition p H = g .
Remark 1. The design of a ”good” eigenvalue problem is a key ingredient of the GMsFEM. It plays a critical role in reducing the dimension of the coarse solution space. After applying such eigenvalue problem, a few dominant modes are selected to capture the multiscale behavior of the multiscale solution up to some certain extent.
Remark 2. Generally, the offline space should contain a constant basis function. Therefore, the local snapshot space V snap v , i , as a solution space of the eigenvalue Problem (9), can be decomposed into V 1 v , i V 2 v , i , where V 1 v , i = { ( 1 , 1 , , 1 ) } and V 2 v , i = { ψ : E i ψ · n i = 0 , ψ V snap v , i } . We note that the element with { ( 1 , 1 , , 1 ) } in the snapshot space corresponds to the local solution with the unit normal flux on the edge and the elements of V 2 v , i have average zero normal flux traces. We eliminate the constant boundary condition basis function in the snapshot space and put it, ( 1 , 1 , , 1 ) V 1 v , i , as the first basis function of the offline space, and identify the other dominant modes in the space V 2 v , i .

3.2. Multiscale Solution Space V H q

We construct the solution space V H q and the testing space W H q for q H . They follow a similar procedure to the one for v H . In order to build V H q , the snapshot space V snap q is to be constructed. We define a local problem corresponding to each coarse-grid edge E i E H as follows: Find ( c j i , ψ q , j i ) Q h × V h q , such that
ψ q , j i + D c j i - v H c j i = 0 in K , · ψ q , j i = α j i in K , ψ q , j i · n = 0 on K E i , ψ q , j i · n = δ j i on E i ,
where K = K r or K s depicted in Figure 2. The local snapshot space V snap q , i , the snapshot space V snap q , and the snapshot matrix R snap q are formed as
V snap q , i = span { ψ q , j i : 1 j J i } ,
V snap q = E i V snap q , i = span { ψ q , j i : 1 j J i , 1 i N e } .
R snap q = ψ 1 q , ψ 2 q , , ψ L snap q ,
respectively.
To construct the testing space W H q , a snapshot space W snap q is also needed. We solve the adjoint problem of (12)
w q , j i - D p j i = 0 in K , · w q , j i + 1 D v H · w q , j i = α j i in K , w q , j i · n = 0 on K E i , w q , j i · n = δ j i on E i .
Then, the snapshot space W snap q is formed by the snapshot functions ω q , j i for all E i T H .
The local offline space V off q , i for u H is obtained by solving the following eigenvalue problem in V snap q , i :
A q off Ψ k q = λ k M q off Ψ k q ,
where
A q off = [ a s t ] = E i ( ψ q , s i · n ) ( ψ q , t i · n ) ,
and
M q off = [ m s t ] = ω i D ˜ ψ q , s i · ψ q , t i .
Here, we let D ˜ = 1 + | v H | D , where | v H | denotes the vector l 2 norm of v H . Then, the local offline space V off q , i , the offline space V off q , and the offline matrix R off q are formed as:
V off q , i = span { ϕ q , k i : 1 k L off i } ,
V off q = E i V off q , i = span { ϕ q , k i : 1 k L off i , 1 i N e } ,
R off q = ϕ 1 q , ϕ 2 q , , ϕ L off q ,
respectively. Here ϕ q , k i = j = 1 J i Ψ k j q ψ q , j i , k = 1 , 2 , , L off i , is the dominant mode in the space V snap q , i .
We let V H q = V off q and define the testing space as
W H q = { w W snap q : ( w - q ) · n | E = 0 , for all E E H , for some q V H q } ,
where n is the unit-normal vector of the coarse-grid edge E. Now, we have the multiscale solution space V H q and the testing space W H q . The GMsFEM solution ( c H , q H ) can be computed through (6).
Remark 3. We remark that the fine-scale velocity field is constructed and used in computing multiscale basis functions for the concentration. In future, we plan to study joint multiscale basis construction procedures, which can compute multiscale basis functions for the concentration without solving the velocity field. Because the concentration field is a nonlinear function of the velocity field, an online procedure for computing multiscale basis functions for the concentration will be employed.

4. Numerical Results

We present a representative set of numerical experiments that demonstrate the performance of the mixed GMsFEM for approximating the coupled flow and transport Equation (2). The computational domain Ω = [ 0 , 1 ] × [ 0 , 1 ] . We consider two different permeability fields κ, as depicted in Figure 3a,b. The source term f c to be used are plotted in Figure 3c,d, and f q is chosen to be 0. In all experiments below, the fine grid T h is fixed to be 200 × 200 uniform mesh, i.e., h = 1 / 200 . The coarse grid consists of 20 × 20 uniform meshes, i.e., H = 1 / 20 . As for the time discretization, the time step is chosen to be τ = 1 × 10 - 4 and the Crank-Nicolson scheme is used. We denote c h ¯ as the average value of c h in each coarse cell, i.e., c h ¯ = 1 | K | K c h , for every K T H . The relative numerical errors are measured in L 2 norm for the concentration.
Figure 3. Two permeability fields (a,b) and two source terms (c,d).
Figure 3. Two permeability fields (a,b) and two source terms (c,d).
Computation 03 00670 g003
In our numerical experiments, we approximate the velocity v H and the flux q H (at any time instant) in the offline space V H v and V H q in which each coarse edge is equipped with L off i basis functions. We choose at most L off i = 2 and compare concentration profiles with L off i = 1 . Note that there are H / h basis functions associated with each coarse edge. We only select up to the two dominant modes per edge to form the offline space.
The first experiment is for the Problem (2), where the diffusivity D = 1 , the source term f c = f 1 , the permeability field κ = κ 1 , and the boundary condition g = 7 x y . The boundary condition represents the flow from lower left corner to upper right corner. The Peclet number in the entire domain is 180 . 18 , while the coarse-grid Peclet numbers vary with the maximum one being 11 . 01 . In Figure 4, we depict the averaged fine-scale solutions and the GMsFEM solutions at the time T = 0 . 015 , 0 . 04 , and 0 . 1 . We choose the averaged fine-scale solution since piecewise constant basis functions are used in the approximation of the concentration on a coarse grid. The averaged fine-scale solution is depicted on the top of Figure 4. In the middle row, we depict the GMsFEM solution, where only one basis function per edge is used for the approximation of the velocity and one basis function is used for the approximation of the flux-concentration field. The L 2 errors in the concentration field are 21 . 22 %, 30 . 16 %, and 27 . 19 % at T = 0 . 015 , 0 . 04 , and 0 . 1 . At the bottom row of Figure 4, we depict the concentrations computed using two basis functions for the velocity and two basis functions for the flux-concentration field. As we observe, the results look much better when two basis functions are used. For example, at T = 0 . 1 , we can observe that the concentration profile is similar to the averaged fine-scale concentration profile. The L 2 errors in the concentration field are 8 . 1 %, 11 . 1 %, and 8 . 44 % at T = 0 . 015 , 0 . 04 , and 0 . 1 . Note that since we use 20 × 20 coarse grid and piecewise constant basis functions for the concentration, one expects a first order convergence with respect to the coarse mesh (i.e., about 5% error in our case).
Figure 4. First row: Fine-scale solutions c h ¯ . Second row: GMsFEM solutions c H (one multiscale basis function per edge for both velocity and flux-concentration). Third row: GMsFEM solutions c H (two multiscale basis function per edge for both velocity and flux-concentration). From left column to right column: T = 0 . 015 , 0 . 04 , 0 . 1 . We use κ = κ 1 , f c = f 1 , H = 1/20.
Figure 4. First row: Fine-scale solutions c h ¯ . Second row: GMsFEM solutions c H (one multiscale basis function per edge for both velocity and flux-concentration). Third row: GMsFEM solutions c H (two multiscale basis function per edge for both velocity and flux-concentration). From left column to right column: T = 0 . 015 , 0 . 04 , 0 . 1 . We use κ = κ 1 , f c = f 1 , H = 1/20.
Computation 03 00670 g004
In our second set of numerical examples, we consider the permeability κ = κ 2 , a source term f c = f 2 , and a boundary condition g = 8 ( x - y ) . We let the diffusivity to be D = 1 . In this case, the problem becomes convection-dominated, where the coarse-grid Peclet number has a maximum value around 8.1. The numerical results are shown in Figure 5, where the top row shows the averaged fine-scale concentration profiles, the middle row shows the GMsFEM solution with one velocity basis per edge and one flux-concentration basis function, and the bottom row shows the GMsFEM solution, which uses two velocity basis functions and two flux-concentration basis functions. In all numerical results, we use piecewise constant pressure and piecewise constant concentration basis functions, as mentioned earlier. Thus, for the concentration field, there is an irreducible error of order H (5% in our case). First, we study the numerical results using only one multiscale basis function (the middle row of Figure 5). The L 2 errors in the concentration field are 77 . 6 %, 82 . 53 %, and 77 . 96 % at T = 0 . 01 , 0 . 02 , and 0 . 1 . We observe from the figure that the results with one basis function per edge do not capture the averaged fine-scale solution behavior correctly. On the other hand, we observe a good agreement when two multiscale basis functions are used (the bottom row of Figure 5). We observe that the main features of the concentration profile are captured. The L 2 errors in the concentration field are 14 . 9 %, 18 . 24 %, and 15 . 71 % at T = 0 . 01 , 0 . 02 , and 0 . 1 .
In our numerical simulations, we limit ourselves to piecewise constant (lowest order) concentration basis functions. This allows a low computational cost and avoids a design of more complex concentration basis functions, which can satisfy inf-sup conditions (see e.g., [32]). Since we use lowest order basis functions for the concentration, the accuracy is limited to O ( H ) . This is reasonable as in large-scale simulations (due to very detailed κ), the coarse grid sizes can be taken sufficiently small. In general, the numerical results shown above can be improved using more basis functions for the flux-concentration and designing multiscale basis functions for the concentration field. One of our main objectives is to show that the results with one basis function per edge (which is similar to flow-based upscaling) can be substantially improved if one uses more multiscale basis functions, and we also provide a procedure for computing the multiscale basis functions and suitable global formulations.
Figure 5. First row: Fine-scale solutions c h ¯ . Second row: GMsFEM solutions c H (one multiscale basis function per edge for both velocity and flux-concentration). Third row: GMsFEM solutions c H (two multiscale basis function per edge for both velocity and flux-concentration). From left column to right column: T = 0 . 01 , 0 . 02 , 0 . 1 . We use κ = κ 2 , f c = f 2 , H = 1/20.
Figure 5. First row: Fine-scale solutions c h ¯ . Second row: GMsFEM solutions c H (one multiscale basis function per edge for both velocity and flux-concentration). Third row: GMsFEM solutions c H (two multiscale basis function per edge for both velocity and flux-concentration). From left column to right column: T = 0 . 01 , 0 . 02 , 0 . 1 . We use κ = κ 2 , f c = f 2 , H = 1/20.
Computation 03 00670 g005

5. Discussions on Time Dependent Velocity

In our examples above, the velocity field was time-independent, and thus, the coupling between the flow and transport is weak. In above examples, we derived a multiscale coupling framework and demonstrated a robust framework for multiscale convection-diffusion equation. When the velocity field becomes time-dependent, one needs to re-calculate multiscale basis functions as velocity varies. There are several procedures, which we will mention, after a brief discussion on time varying velocity field scenarios. First, we remark on several scenarios, when the velocity field can depend on time (we call it “time-dependent” in our further discussions, even though the time may enter implicitly). The velocity field can depend on time via time-varying concentration field as in miscible or immiscible flows (e.g., [33] and the referenes therein) or the velocity field can, in addition, have a time-dependent component due to compressibility. We plan to address these application problems in our future works.
Computations of multiscale basis functions for the time-dependent velocity field can be performed offline [16,34] or using the residuals in the online stage [12,13]. For the offline computations, one needs a richer class of snapshot vectors, which can span possible velocity fields and perform the same local spectral decompositions and global coupling framework, as we discussed above. “Rich” snapshot spaces are discussed in [16] for parameter-dependent problem (the time can be treated as a parameter). Another approach for identifying multiscale basis functions, when the snapshot space is very large, is the use of l 1 -minimization techniques as proposed in [34]. All these approaches can be used adaptively in space and time. Finally, one can use residual-based online basis functions [12,13], where one additional multiscale basis function is computed adaptively using the residual information. The success of these approaches depends on the offline space as it is shown, thus a good offline space is needed to achieve a high accuracy.
In some time-dependent velocity cases, one can get a good accuracy without modifying multiscale basis functions. For simplicity, we consider a “miscible” type flow and transport, which are described by
c t + v · c - D Δ c = f c in Ω × ( 0 , T ) , v + ζ ( c ) κ p = 0 in Ω , · v = f q in Ω .
We take ζ ( c ) = ( c 2 + 0 . 2 ( 1 - c ) 2 ) . We use the same setup as in the first example. The source term is the same, that is, f c = f 1 and the boundary condition is g = 28 x y . The numerical results are shown in Figure 6. The first row shows the averaged fine-scale concentration profiles and the second row shows the GMsFEM solution with two velocity basis and two flux-concentration basis functions. The L 2 errors in the concentration field are 5 . 81 %, 8 . 63 %, and 7 . 72 % at T = 0 . 015 , 0 . 04 , and 0 . 1 , respectively. These errors are comparable to the case with time-independent velocity field. We note that the errors are close to the irreducible error due to the coarse-mesh size.
Figure 6. First row: Fine-scale solutions c h ¯ . Second row: GMsFEM solutions c H (two multiscale basis function per edge for both velocity and flux-concentration). From left column to right column: T = 0 . 015 , 0 . 04 , 0 . 1 . We use κ = ( c 2 + 0 . 2 ( 1 - c ) 2 ) κ 1 , f c = f 1 , H = 1 / 20 . Time-dependent velocity case.
Figure 6. First row: Fine-scale solutions c h ¯ . Second row: GMsFEM solutions c H (two multiscale basis function per edge for both velocity and flux-concentration). From left column to right column: T = 0 . 015 , 0 . 04 , 0 . 1 . We use κ = ( c 2 + 0 . 2 ( 1 - c ) 2 ) κ 1 , f c = f 1 , H = 1 / 20 . Time-dependent velocity case.
Computation 03 00670 g006

6. Conclusions

In this paper, we develop a Generalized Multiscale Finite Element Method (GMsFEM) for coupled flow and transport equations. The transport equation is convection dominated and we choose an appropriate test space to achieve a stability and improve numerical results. We study a mixed formulation for both flow and transport, which guarantees mass conservation. The multiscale spaces for the flux and the velocity fields are constructed by appropriately choosing the snapshot spaces and performing local spectral decompositions. The design of the local test and trial multiscale spaces is one of our novel contributions. From the numerical experiments, we observe that one needs to have more than one multiscale basis function to capture the concentration profile on the coarse grid. In particular, using only two multiscale basis functions for the pressure equation and two multiscale basis functions for the concentration equation, we can achieve good accuracy. Our future studies will include a systematic design of multiscale basis functions for the concentration that can provide a stable discretization.

Acknowledgments

Eric Chung’s research is supported by the Hong Kong RGC General Research Fund project 400411 and CUHK Faculty of Science Research Incentive Scheme 2015-16. This publication was made possible by NPRP grant NPRP 7-1482-1-278 from the Qatar National Research Fund (a member of Qatar Foundation). The statements made herein are solely the responsibility of the authors. Yalchin Efendiev’s work is partially supported by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award Number DE-FG02-13ER26165 and the DoD Army ARO Project.

Author Contributions

Jun Ren was a lead in the paper. Wing Tat Leung, Eric T. Chung, Yalchin Efendiev contributed to the computations and method developments.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Chen, Y.; Durlofsky, L. An Ensemble Level Upscaling Approach for eFficient Estimation of fIne-Scale Production Statistics Using Coarse-Scale Simulations. In SPE Reservoir Simulation Symposium; Society of Petroleum Engineers: Houston, TX, USA, 2007. [Google Scholar]
  2. Chen, Y.; Durlofsky, L.; Gerritsen, M.; Wen, X. A coupled local-global upscaling approach for simulating flow in highly heterogeneous formations. Adv. Water Resour. 2003, 26, 1041–1060. [Google Scholar] [CrossRef]
  3. Durlofsky, L.J. Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media. Water Resour. Res. 1991, 27, 699–708. [Google Scholar] [CrossRef]
  4. Aarnes, J.E.; Krogstad, S.; Lie, K.-A. A hierarchical multiscale method for two-phase flow based upon mixed finite elements and nonuniform grids. Multiscale Model. Simul. 2006, 5, 337–363. [Google Scholar] [CrossRef]
  5. Arbogast, T.; Pencheva, G.; Wheeler, M.F.; Yotov, I. A multiscale mortar mixed finite element method. Multiscale Model. Simul. 2007, 6, 319–346. [Google Scholar] [CrossRef]
  6. Brown, D.L.; Peterseim, D. A multiscale method for porous microstructures. ArXiv E-Prints 2014. arXiv:1411.1944. [Google Scholar]
  7. Bush, L.; Ginting, V.; Presho, M. Application of a conservative, generalized multiscale finite element method to flow models. J. Comput. Appl. Math. 2014, 260, 395–409. [Google Scholar] [CrossRef]
  8. Chan, H.Y.; Chung, E.T.; Efendiev, Y. Adaptive mixed gmsfem for flows in heterogeneous media. ArXiv E-Prints 2015. arXiv:1507.01659. [Google Scholar]
  9. Chung, E.; Efendiev, Y. Reduced-contrast approximations for high-contrast multiscale flow problems. Multiscale Model. Simul. 2010, 8, 1128–1153. [Google Scholar] [CrossRef]
  10. Chung, E.; Leung, W.T. A sub-grid structure enhanced discontinuous galerkin method for multiscale diffusion and convection-diffusion problems. Comput. Phys. 2013, 14, 370–392. [Google Scholar] [CrossRef]
  11. Chung, E.T.; Efendiev, Y.; Fu, S. Generalized multiscale finite element method for elasticity equations. GEM-Int. J. Geomath. 2014, 5, 225–254. [Google Scholar] [CrossRef]
  12. Chung, E.T.; Efendiev, Y.; Leung, W.T. An online generalized multiscale discontinuous galerkin method (gmsdgm) for flows in heterogeneous media. ArXiv E-Prints 2015. arXiv:1504.04417. [Google Scholar]
  13. Chung, E.T.; Efendiev, Y.; Leung, W.T. Residual-driven online generalized multiscale finite element methods. ArXiv E-Prints 2015. arXiv:1501.04565. [Google Scholar] [CrossRef]
  14. Chung, E.T.; Efendiev, Y.; Li, G.; Vasilyeva, M. Generalized multiscale finite element methods for problems in perforated heterogeneous domains. Appl. Anal. 2015, in press. [Google Scholar] [CrossRef]
  15. Chung, E.T.; Efendiev, Y.; Li, G. An adaptive GMsFEM for high-contrast flow problems. J. Comput. Phys. 2014, 273, 54–76. [Google Scholar] [CrossRef]
  16. Efendiev, Y.; Galvis, J.; Hou, T. Generalized multiscale finite element methods. J. Comput. Phys. 2013, 251, 116–135. [Google Scholar] [CrossRef]
  17. Efendiev, Y.; Galvis, J.; Lazarov, R.; Willems, J. Robust domain decomposition preconditioners for abstract symmetric positive definite bilinear forms. ESAIM Math. Model. Numer. Anal. 2012, 46, 1175–1199. [Google Scholar] [CrossRef]
  18. Efendiev, Y.; Galvis, J.; Vassilevski, P.S. Spectral Element Agglomerate Algebraic Multigrid Methods for Elliptic Problems with High-Contrast Coefficients. In Domain Decomposition Methods in Science and Engineering XIX; Springer: Heidelberg, Germany, 2011. [Google Scholar]
  19. Efendiev, Y.; Hou, T.; Wu, X.H. Convergence of a nonconforming multiscale finite element method. SIAM J. Numer. Anal. 2000, 37, 888–910. [Google Scholar] [CrossRef]
  20. Hajibeygi, H.; Kavounis, D.; Jenny, P. A hierarchical fracture model for the iterative multiscale finite volume method. J. Comput. Phys. 2011, 230, 8729–8743. [Google Scholar] [CrossRef]
  21. Henning, P.; Peterseim, D. Oversampling for the multiscale finite element method. Multiscale Model. Simul. 2013, 11, 1149–1175. [Google Scholar] [CrossRef]
  22. Hou, T.; Wu, X.H. A multiscale finite element method for elliptic problems in composite materials and porous media. J. Comput. Phys. 1997, 134, 169–189. [Google Scholar] [CrossRef]
  23. Jenny, P.; Lee, S.H.; Tchelepi, H. Multi-scale finite volume method for elliptic problems in subsurface flow simulation. J. Comput. Phys. 2003, 187, 47–67. [Google Scholar] [CrossRef]
  24. Kaulmann, S.; Flemisch, B.; Haasdonk, B.; Lie, K.A.; Ohlberger, M. The localized reduced basis multiscale method for two-phase flows in porous media. Int. J. Numer. Methods Eng. 2015, 102, 1018–1040. [Google Scholar] [CrossRef]
  25. Målqvist, A.; Peterseim, D. Localization of elliptic multiscale problems. Math. Comput. 2014, 83, 2583–2603. [Google Scholar] [CrossRef]
  26. Presho, M.; Galvis, J. A mass conservative generalized multiscale finite element method applied to two-phase flow in heterogeneous porous media. ArXiv E-Prints 2015. arXiv:1504.02033. [Google Scholar] [CrossRef]
  27. Arbogast, T.; Wheeler, M.F.; Yotov, I. Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences. SIAM J. Numer. Anal. 1997, 34, 828–852. [Google Scholar] [CrossRef]
  28. Douglas, J., Jr.; Ewing, R.E.; Wheeler, M.F. The approximation of the pressure by a mixed method in the simulation of miscible displacement. RAIRO-Anal. Numérique 1983, 17, 17–33. [Google Scholar]
  29. Durlofsky, L.J. A triangle based mixed finite element-finite volume technique for modeling two phase flow through porous media. J. Comput. Phys. 1993, 105, 252–266. [Google Scholar] [CrossRef]
  30. 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. 1984, 47, 73–92. [Google Scholar] [CrossRef]
  31. Calo, V.M.; Chung, E.T.; Efendiev, Y.; Leung, W.T. Multiscale stabilization for convection-dominated diffusion in heterogeneous media. ArXiv E-Prints 2015. arXiv:1509.06833. [Google Scholar]
  32. Chung, E.T.; Efendiev, Y.; Lee, C.S. Mixed generalized multiscale finite element methods and applications. Multiscale Model. Simul. 2015, 13, 338–366. [Google Scholar] [CrossRef]
  33. Ewing, R.; Efendiev, Y.; Ginting, V.; Wang, H. Upscaling of Transport Equations for Multiphase and Multicomponent Flows; Springer: Heidelberg, Germany, 2008. [Google Scholar]
  34. Chung, E.; Efendiev, Y.; Leung, T.W.; Li, G. Sparse generalized multiscale finite element methods and their applications. ArXiv E-Prints 2015. arXiv:1506.08509. [Google Scholar]

Share and Cite

MDPI and ACS Style

Chung, E.T.; Efendiev, Y.; Leung, W.T.; Ren, J. Multiscale Simulations for Coupled Flow and Transport Using the Generalized Multiscale Finite Element Method. Computation 2015, 3, 670-686. https://0-doi-org.brum.beds.ac.uk/10.3390/computation3040670

AMA Style

Chung ET, Efendiev Y, Leung WT, Ren J. Multiscale Simulations for Coupled Flow and Transport Using the Generalized Multiscale Finite Element Method. Computation. 2015; 3(4):670-686. https://0-doi-org.brum.beds.ac.uk/10.3390/computation3040670

Chicago/Turabian Style

Chung, Eric T., Yalchin Efendiev, Wing Tat Leung, and Jun Ren. 2015. "Multiscale Simulations for Coupled Flow and Transport Using the Generalized Multiscale Finite Element Method" Computation 3, no. 4: 670-686. https://0-doi-org.brum.beds.ac.uk/10.3390/computation3040670

Article Metrics

Back to TopTop