# A mixed multiscale finite element method for elliptic problems with oscillating coefficients

## Abstract

The recently introduced multiscale finite element method for solving elliptic equations with oscillating coefficients is designed to capture the large-scale structure of the solutions without resolving all the fine-scale structures. Motivated by the numerical simulation of flow transport in highly heterogeneous porous media, we propose a mixed multiscale finite element method with an over-sampling technique for solving second order elliptic equations with rapidly oscillating coefficients. The multiscale finite element bases are constructed by locally solving Neumann boundary value problems. We provide a detailed convergence analysis of the method under the assumption that the oscillating coefficients are locally periodic. While such a simplifying assumption is *not* required by our method, it allows us to use homogenization theory to obtain the asymptotic structure of the solutions. Numerical experiments are carried out for flow transport in a porous medium with a random log-normal relative permeability to demonstrate the efficiency and accuracy of the proposed method.

## 1. Introduction

Many problems of fundamental and practical importance in science and engineering have multiple-scale solutions. Typical examples include composite materials with fine micro-structures and highly heterogeneous porous media. The direct numerical simulation of problems involving multiscale solutions is difficult, due to the requisite of tremendous amount of computer memory and CPU time, which can easily exceed the limit of today’s computer resources. On the other hand, in practice, it is often sufficient to predict the large scale solutions to a certain accuracy. Thus, various methods of upscaling or homogenization have been developed which replace the governing equations with multiscale solutions by the homogenized equations, whose solutions can be resolved on a coarse-scale mesh.

The recently developed multiscale finite element method Reference 18, Reference 19, Reference 13 provides an effective way to capture the large-scale structures of the solutions on a coarse mesh. The central idea of the method is to incorporate the local small-scale information of the leading-order differential operator into the finite element bases. It is through these multiscale bases and the finite element formulation that the effect of small scales on the large scales is correctly captured. Apart from the computational advantages, such as saving in computer memory and good parallel efficiency Reference 18 of the multiscale finite element method, we stress that the real significance of the method lies in its ability to solve the problems in coarse meshes. This is particularly advantageous when multiple simulations or realizations are necessary due to changes of boundary conditions or source functions for given fine microstructures of composite materials or highly heterogeneous permeability of porous media. Other relevant works on constructing special finite element bases can be found in Reference 3 for layered microstructures and in Reference 6 for convection-dominated diffusion problems.

The motivation of this paper is to explore the possibility of applying the multiscale finite element method to numerical computation of flow transport in highly heterogeneous porous media. In its simplest form, neglecting the effect of gravity, compressibility, capillary pressure, and considering constant porosity and unit mobility, the governing equations for the flow transport can be described by the following partial differential equations (see Reference 21, Reference 31, and Reference 12):

where is the pressure, is the water saturation, is the relative permeability tensor, and is the Darcy velocity. The highly heterogeneous properties of the medium are built into the permeability tensor which is generated through the use of sophisticated geological and geostatistical modeling tools. The detailed structure of the permeability coefficients makes the direct simulation of the above model infeasible. For example, it is common in real simulations to use millions of grid blocks, with each block having a dimension of tens of meters, whereas the permeability measured from cores is at a scale of centimeters ,Reference 24. This gives more than degrees of freedom per spatial dimension in the computation. This makes a direct simulation to resolve all small scales prohibitive even with today’s most powerful supercomputers. On the other hand, from an engineering perspective, it is often sufficient to predict the macroscopic properties of the solutions. Thus it is highly desirable to derive effective coarse grid models to capture the correct large solution without resolving the small-scale features. Numerical upscaling is one of the commonly used approaches in practice. There are extensive studies in the literature on the numerical upscaling of two-phase flows through porous media, see, e.g., Reference 1, Reference 21, Reference 11, Reference 14, Reference 30.

Mixed finite element approximations for second order elliptic problems, which approximate the source variable and flux simultaneously, have been studied by many authors (cf., e.g., Reference 27 and the book Reference 5). The local conservation of velocity flux is an important property in the mixed finite element methods. The violation of this local conservation property will lead to leakage of velocity flux. This will deteriorate the accuracy of the numerical solution for long-time computations. This is the reason why mixed finite element methods are attractive for porous medium simulations Reference 28 (see also the numerical results in §5.1 of this paper). In this paper, we first propose and analyze a mixed multiscale finite element method with an over-sampling technique for solving elliptic equations with oscillating coefficients, and then apply the method to compute the above flow model Equation 1.1-Equation 1.2. The use of the over-sampling technique is crucial in eliminating large resonance errors from the element boundaries. This point has already been well demonstrated in the previous study Reference 18Reference 13 for the displacement multiscale finite element method. We will demonstrate that this technique applies also to the mixed multiscale finite element method. Our computational results show that the mixed multiscale method with the over-sampling technique gives more accurate results than the corresponding method without such technique.

We have performed careful numerical experiments to validate our mixed multiscale finite element method. To test how well our method works in a realistic application, we apply our method to a flow in porous media with a random log-normal relative permeability tensor. This is one of the practical benchmark test cases for oil recovery problems. Our computational results demonstrate convincingly the importance of the local conservation property in the flow simulation. When the local conservation property is not satisfied, the fractional flow curve deviates significantly from that obtained from the well-resolved mixed finite element calculation after a short time. The mixed multiscale finite element method, on the other hand, provides an accurate numerical approximation even on a coarse grid, with accuracy comparable to that obtained using the standard fine-mesh mixed finite element method. Finally, we compare the performance of the mixed multiscale method studied in this paper with the displacement multiscale finite element method studied in Reference 18, Reference 19, Reference 13 when applied to the flow transport problem in heterogeneous porous media. The two methods give similar results at an early stage. But the fractional flow curve obtained by the displacement multiscale method deviates significantly from that obtained by the well-resolved calculation after a short time, clearly suffering from the violation of the local conservation property. The mixed multiscale method gives a more accurate solution than the displacement multiscale method for long-time calculations, due to the inherent local conservation property of the mixed multiscale method.

The ultimate goal of this study is to produce an effective coarse grid model for the two-phase flow with heterogeneous porous media. To this end, we need to upscale the saturation equation. Without capillary pressure, the saturation equation is hyperbolic. The effective equation is difficult to derive and has a nonlocal memory effect Reference 29. Using an upscaling technique recently developed in Reference 14 for the saturation equation, we show how the proposed mixed multiscale finite element method leads to a complete coarse grid algorithm. Our numerical results demonstrate convincingly that the fractional flow curve obtained from the resulting coarse grid model gives a very good approximation to that obtained from the fine grid calculation. Typically, many realizations on the same microstructure (permeability field) are made due to changes in boundary conditions and source fields. In such a case, the multiscale finite element bases are only constructed once, and can be used in the subsequent computations. Thus, the coarse grid model offers substantial saving in both memory and computing time. The saving could be as large as a factor of 10,000 if one can scale up by a factor of 10 in each space dimension (three space dimensions plus time).

The outline of the paper is as follows. In §2 we introduce the mixed multiscale finite element methods and present the main convergence results. In §3 we review the homogenization results for elliptic problems with Neumann boundary conditions. These results are the basis of our convergence analysis. In §4 we prove the error estimate for the mixed multiscale finite element method with over-sampling introduced in §2.3. The analysis depends on an abstract formulation for nonconforming mixed finite element methods, the homogenization theory, and a technique of “freezing coefficients” to deal with the local periodic coefficients. In §5 we apply our method to simulate the flow transport model Equation 1.1-Equation 1.2 for a practical random log-normal permeability.

## 2. Mixed multiscale finite element methods

### 2.1. Notation and background

Let be a polyhedral domain in with boundary whose unit outer normal is denoted by Let . be a domain in with the Lipschitz boundary For each integer . and we denote by , the standard Sobolev space of real functions having all their weak derivatives of order up to in the Lebesgue space The norm and the seminorm of . will be denoted by and respectively. As usual, when , , is denoted by with the norm and the seminorm In the following, we let . denote the subspace of whose functions have zero average over .

We consider the following second order elliptic equations with locally periodic coefficients Reference 4:

In this paper the usual Einstein summation convention for repeated indices is used. Here, is assumed to be a small parameter, and is a symmetric matrix which satisfies the uniform ellipticity condition:

for some positive constant Furthermore, we assume that . where , stands for the collection of all periodic functions with respect to the unit cube .

Let be the solution of the homogenized problem of Equation 2.1-Equation 2.2:

where with

and is the periodic solution of the cell problem

Here is the Kronecker delta, i.e., for and for Note that in .Equation 2.7 plays the role of a parameter. However, since is differentiable in we can easily show that , is also differentiable in The convergence property of . to will be studied in detail in the next section.

Throughout the paper we impose the following assumptions on the data:

(H1) , on for some .

(H2) Compatibility: .

We remark that the assumption on the boundary value is not very restrictive. For example, in the case when is a convex polygon in with sides and , such a flux , can be constructed through with being the solution of the equation in with the Neumann boundary condition on The regularity results in .Reference 17 then ensure We also stress that the explicit formula of such a flux . is not required in the definition of our methods.

Denote by the subspace of whose functions have zero average over and let , be the subspace of given by

The norm of will be denoted by Let . and the inverse matrix of then ; The mixed formulation to .Equation 2.1-Equation 2.2 then reads as follows: find a pair such that on and

Here stands for the inner product of or The existence of a unique solution of problem .Equation 2.8-Equation 2.9 and its equivalence to Equation 2.1-Equation 2.2 follow directly from the abstract Babuska-Brezzi theory Reference 5, Reference 27.

Let us assume that is a regular and quasi-uniform partition of into simplices. For any let , be its diameter, its Lebesgue measure, the unit outer normal to and , the surfaces or edges of with being the measure of The mixed multiscale finite element method that we are going to introduce is closely related to the Raviart-Thomas elements .Reference 27. For any we denote by , the lowest order Raviart-Thomas approximation of :

where is the constant element space. Let be the basis of which satisfies

Since

and the error estimate

Let

The following well-known error estimate can be obtained by using the Babuska-Brezzi theory Reference 27, Reference 5

Note that since

### 2.2. Mixed multiscale finite element method

Recalling that

Equation Equation 2.13 is the weak formulation of the following boundary value problem

Now let

Moreover, for any

imply

Let

Let

We can now introduce the following discretization of Equation 2.8-Equation 2.9: find a pair