Geofluids

Volume 2017 (2017), Article ID 1740693, 12 pages

https://doi.org/10.1155/2017/1740693

## Upscaled Lattice Boltzmann Method for Simulations of Flows in Heterogeneous Porous Media

^{1}Center for Numerical Porous Media, King Abdullah University of Science and Technology, Thuwal, Saudi Arabia^{2}Center for Integrative Petroleum Research, College of Petroleum Engineering and Geosciences, King Fahd University of Petroleum & Minerals, Dhahran, Saudi Arabia^{3}School of Mathematical Sciences, University of Nottingham, Nottingham, UK

Correspondence should be addressed to Jun Li

Received 30 June 2016; Accepted 5 January 2017; Published 16 February 2017

Academic Editor: Andri Stefansson

Copyright © 2017 Jun Li and Donald Brown. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

An upscaled Lattice Boltzmann Method (LBM) for flow simulations in heterogeneous porous media at the Darcy scale is proposed in this paper. In the Darcy-scale simulations, the Shan-Chen force model is used to simplify the algorithm. The proposed upscaled LBM uses coarser grids to represent the average effects of the fine-grid simulations. In the upscaled LBM, each coarse grid represents a subdomain of the fine-grid discretization and the effective permeability with the reduced-order models is proposed as we coarsen the grid. The effective permeability is computed using solutions of local problems (e.g., by performing local LBM simulations on the fine grids using the original permeability distribution) and used on the coarse grids in the upscaled simulations. The upscaled LBM that can reduce the computational cost of existing LBM and transfer the information between different scales is implemented. The results of coarse-grid, reduced-order, simulations agree very well with averaged results obtained using a fine grid.

#### 1. Introduction

Detailed flow simulations in porous media are often modeled using the Darcy or Brinkman approximations. In these models, effective parameters, such as absolute and relative permeabilities, depend on the pore-scale geometry. To compute these effective parameters, pore-scale simulations accounting for relevant geometric features in a Representative Elementary Volume (REV) are commonly used as in [1]. The Lattice Boltzmann Method (LBM) [2–4] is well developed for pore-scale flow simulations and extended to model two-phase systems or two immiscible fluids [5–7]. After computing the effective parameters, we are able to perform Darcy-scale simulations using traditional finite volume or element methods used in commercial reservoir simulators. However, these computations are limited to small REVs (compared to the computational grid) and rely on two distinct idealized scale concepts.

Flows at the Darcy scale can also be modeled by LBM with a modified algorithm. The model described in [8] allows particles to partially bounce back at the cells (points) with small permeability. In [9, 10], an external body force, which increases with decreasing permeability, is employed to represent the resistance effect of the porous media to the fluid, where LBM is considered as a unified framework for simulations at all scales. However, these simulations require significant computational resources to converge since the permeability distribution usually has drastic changes in space, which requires a very fine grid for high spatial resolution. To overcome this difficulty, we propose an upscaled LBM scheme.

Following the work in [10] where the generalized Navier-Stokes equation [11] is solved, we simplify the equilibrium distribution function to solve the governing equation at the Darcy scale. In addition, we replace the original Guo et al. force model [12] used in [10] by the simpler Shan-Chen force model [5, 6] to improve the efficiency. Then, an upscaled LBM scheme is proposed to improve computational efficiency by using a coarse grid (each coarse point represents a subdomain) with effective permeability. For each subdomain, the effective permeability is computed by a local scheme, which is based on the conservation principle for the average fluxes (see [13] for general overview of multiscale methods). To avoid the iterative process of finding the unknown effective permeability that satisfies the equation for the average flux, we derive an analytical formula. This analytical formula allows finding the average flux in terms of the effective permeability, which is then inversely determined from the computed average flux using the original permeability distribution in the subdomain concerned.

The computed effective permeabilities are verified in several benchmark problems, where analytical solutions are known. We implement upscaled LBM simulations using the computed effective permeabilities at coarse grids. Agreement between the coarse and fine-grid LBM simulations demonstrates the validity of the upscaled LBM scheme. The average effects of the fine-grid simulations are preserved in the coarse-grid simulations in solving the flow equation at any intermediate scale. Our numerical results show that one can achieve a substantial gain in CPU time by using coarse-grid models. In this paper, the upscaled LBM approach is applied to single-phase flows; however, this approach can be used for modeling multiphase flow phenomena.

#### 2. LBM Algorithms for Simulating Flows in Porous Media

In this section, we discuss LBM algorithms that will be used in our microscale simulations. We first present the LBM algorithm based on the force model proposed by Guo et al. in [12], where an additional term is used in the particle evolution equation to represent the force contribution. Then, we present the general Shan-Chen force model for multicomponent and multiphase systems. In our upscaling algorithm, we focus on the single-phase and single-component model for Brinkman flows. The Shan-Chen force model allows for a more efficient upscaling procedure and overall cleaner presentation. We refer to [14] for more general discussions on LBM algorithms.

First, we will introduce some basic notations associated with LBM. The grid (Lattice) points are uniformly distributed inside the computational domain and , where is the Lattice spacing and is the time step. For two-dimensional problems, we use the D2Q9 Lattice model [4], where is the total number of Lattice velocities, , , , , for 1 to 4, , , and for 5 to 8. For three-dimensional problems, the D3Q19 Lattice model [4] with different and is used, but the algorithms are unchanged.

##### 2.1. LBM Algorithm with the Guo Force Model

Following the algorithm presented in [10], the evolution algorithm of the distribution function iswhere the normalized relaxation time is appropriately selected to match the desired effective kinematic viscosity , where is the sound speed. We use the simplified truncated form of the equilibrium distribution function aswhere the density is given by and the equilibrium velocity, , is defined aswhere is the force per unit mass. Similarly, is simplified to

In the force model proposed by Guo et al. [12], the flow velocity is equal to . If is constant, and are computed by and then and of (1) are determined explicitly. For solving the Darcy-scale equation, we consider the following expressions for as a linear function of (see [10]):where is the porosity, is the physical kinematic viscosity of the fluid, is a scalar for the permeability, and is the external body force per unit mass. The force introduced above incorporates the porous media heterogeneities through the permeability function and depends on the microstructure. If has a high value in the region, then one can assume that this region is highly permeable, while if has a very low value, then this region is almost impermeable. One can also use a forcing that is nonlinear in as an extension to cases with nonlinear Forchheimer effects that are discussed in [10, 15].

In expressions (3) and (5) we have . Using this fact and solving for in (5), we obtain the explicit formula as in [10]: is computed by (5). Then, and of (1) are determined by (2) and (4), respectively.

In the incompressible limit with , the analysis [10] based on the Chapman-Enskog expansion shows that the computed pressure and flow velocity converge to the solutions of the following equation:where is the initial mass density used in LBM simulations. Here, needs not be the real density of the incompressible fluid; then, the computed is used as the pressure of the physical problem. The steady state results of LBM simulations are used as the solutions of the Brinkman equation. The parameters of , , , and can be set independently such that the steady state LBM results converge to the solutions of the continuum Darcy and Stokes equations, respectively.

##### 2.2. Simplified LBM Algorithm with the Shan-Chen Force Model

We now present the general Shan-Chen model and its application to our upscaling scheme. In the original Shan-Chen model [5, 6], which is proposed to simulate multiphase and multicomponent flows, the number of molecules of component having the velocity at and time is denoted by , where and is the total number of components. The general updating algorithm of iswhere and the equilibrium distribution function is defined aswhere ,, and is computed aswhere is related to the total volume force acting on th component. Generally speaking [16], contains three parts: the fluid-fluid interaction , fluid-solid interaction , and external force . For example, for the contribution by the external body force per unit mass. In (10), is defined as follows to conserve momentum:The flow velocity of the whole fluid is equal to the mean velocity before and after implementing the force term and is computed as follows:Recently [15], phase separation process in a fiber geometry and flow of two immiscible fluids in a cross channel are modeled using the Shan-Chen model, which shows the convenience of the LBM in dealing with complex geometries and manipulating the contact angle.

As upscaling in the multiphase problems is a very difficult and often nonlinear procedure, we focus our algorithm first to single-component and single-phase models. For flows of single-component and single-phase, the evolution of without notation isIn order to recover the Brinkman equation, the equilibrium distribution function of (9) is simplified to (2). but of (10) is modified as follows:where we use to replace the original notation, which is equal to the momentum increase per unit volume after due to the force effect through the relaxation process. Correspondingly, the flow velocity of (12) is modified toWhen solving the Brinkman equation, is a function of defined by (5). A comparison of (3) and (15) shows that the explicit formula of (6) to compute is also valid here. Then, is computed aswhich is obtained by solving (14) and (15).

As we can see, the computations of by (5) and by (4) using the computed in the original algorithm [10] are avoided in the current simplified algorithm and, therefore, the efficiency is improved. In the incompressible limit with , the computed pressure and flow velocity also converge to the solutions of the above Brinkman-like equation (7). The same idea can be implemented to (8)–(12) to solve multiphase flows at the Darcy scale. In this way, it may be possible to develop multiphase upscaling techniques based on the upscaling scheme presented below. This is a topic for future work.

#### 3. Upscaling Scheme

For many practical cases, the number of fine discretization points in the whole computational domain due to heterogeneities is very large, making the memory usage and computational time unaffordable. We use an upscaling simulation scheme to reduce the number of points in the fine grid by using a coarse grid with an effective permeability . The upscaled quantities are a tensor quantity even though the input permeability, , is assumed to be a scalar. With this approach we are able to capture fine-grid information on the coarse grid by solving many parallel local problems.

In our proposed algorithm, the computational domain is divided into many subdomains and each subdomain is represented by a coarse point (see Figure 1). This substantially reduces the degrees of freedom in the coarse-grid simulation. To compute the effective for each subdomain, we impose different external forces to drive flows in different directions in the local LBM simulations, which use a fine grid located inside the corresponding subdomain and the distribution of on the fine grid. Then, the similar local LBM simulations usually need to be run with a constant tensor as shown in (18). We seek such that the average velocities from local fine-grid simulations with the heterogeneous and homogeneous , respectively, are equal (see (21)-(22)). The onerous seeking process by adjusting the unknown to match the fluxes computed using is avoided in our simulations since can be computed explicitly by (17).