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

A upscaled lattice Boltzmann method (LBM) for flow simulations in heterogeneous porous media, at both pore and Darcy scales, is proposed in this paper. In the micro-scale simulations, we model flows using LBM with the modified Guo et al. algorithm where we replace the force model with a simple Shan-Chen force model. The proposed upscaled LBM uses coarser grids to represent the effects of the fine-grid (pore-scale) simulations. For the upscaled LBM, effective properties and reduced-order models are proposed as we coarsen the grid. The effective properties are computed using solutions of local problems (e.g., by performing local LBM simulations) subject to some boundary conditions. A upscaled LBM that can reduce the computational complexity 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.


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 a upscaled LBM scheme that is applicable at the pore and coarser, e.g., Darcy scales.
Following the work in [10] where the generalized Navier-Stokes equation [11] is solved, we simplify the equilibrium distribution function. 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 a coarse grid and the computed effective permeability. 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 multi-phase flow phenomena.

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 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 Brinkmann 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.

LBM algorithm with the Guo force model
Following the algorithm presented in [10], the evolution algorithm of the distribution function f α ( x, t) is: where the normalized relaxation time τ is appropriately selected to match the desired effective kinematic viscosity ν eff = c 2 s (τ −0.5)∆t, where c s = c/ √ 3 is the sound speed. We use the simplified truncated form of the equilibrium distribution function as: where the density is given by ρ( x, t) = Q−1 α=0 f α ( x, t) and the equilibrium velocity, u (eq) ( x, t), is defined as: where f m ( x, t) is the force per unit mass. Similarly, F α ( x, t) is simplified to: In the force model proposed by Guo et al. [12], the flow velocity u is equal to u (eq) . If f m is constant, ρ and u (eq) are computed by f α and then f (eq) α and F α of Eq. (1) are determined explicitly. For solving the pore-scale equation, we consider the following expressions for f m as a linear function of u (see [10]): where is the porosity, ν is the physical kinematic viscosity of the fluid, κ( x) is a scalar for the permeability and G( x) is the external body force per unit mass. The force introduced above incorporates the porous media heterogeneities through the permeability function κ( x) and depends on the microstructure. If κ( x) has a high value in the region, then one can assume that this region is highly permeable, while if κ( x) has a very low value, then this region is almost impermeable. One can also use a forcing that is nonlinear in u as an extension to cases with nonlinear Forchheimer effects which are discussed in [10] and [15]. In expressions (3) and (5) we have u = u (eq) . Using this fact and solving for u in (5), we obtain the explicit formula as in [10]: .
In the incompressible limit with | u| c s , the analysis [10] based on the Chapman-Enskog expansion shows that the computed pressure p = c 2 s ρ and flow velocity u converge to the solutions of the following equation: where ρ 0 is the initial mass density used in LBM simulations. Here, ρ 0 needs not be the real density ρ real of the incompressible fluid; then, the computed pρ real /ρ 0 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 ν eff , , ν and κ( x) can be set independently such that the steady state LBM results converge to the solutions of the continuum Darcy and Stokes equations, respectively.

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 the σ th component having the velocity e α at x and time t is denoted by f σ α ( x, t), where σ = 1, · · · , S and S is the total number of components. The general updating algorithm of f σ α ( x, t) is: where σ = 1, · · · , S and the equilibrium distribution function is defied as: where σ = 1, · · · , S, ρ σ = Q−1 α=0 f σ α and u σ(eq) is computed as: where F σ ( x, t) is related to the total volume force acting on the σ th component. Generally speaking [16], F σ contains three parts: the fluid-fluid interaction F 1,σ , fluid-solid interaction F 2,σ and external force F 3,σ . For example, F 3,σ = ∆tρ σ G for the contribution by the external body force G per unit mass. In Eq. (10), u is defined as follows to conserve momentum: The flow velocity u 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 phase 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 f α ( x, t) without notation σ is: In order to recover the Brinkman equation, the equilibrium distribution function f (eq) α of Eq. (9) is simplified to the above Eq. (2). ρ = Q−1 α=0 f α but u (eq) of Eq. (10) is modified as follows: where we use ∆tρ f m to replace the original notation which is equal to the momentum increase per unit volume after ∆t due to the force effect through the relaxation process. Correspondingly, the flow velocity u of Eq. (12) is modified to: When solving the Brinkman equation, f m is a function of u defined by Eq. (5). A comparison of Eqs. (3) and (15) shows that the explicit formula of Eq. (6) to compute u is also valid here. Then, u (eq) is computed as: which is obtained by solving Eqs. (14) and (15).
As we can see, the computations of f m by Eq. (5) and F α by Eq. (4) using the computed f m in the original algorithm [10] are avoided in the current simplified algorithm and, therefore, the efficiency is improved. In the incompressible limit with | u| c s , the computed pressure p = c 2 s ρ and flow velocity u also converge to the solutions of the above Brinkman-like equation Eq. (7). The same idea can be implemented to Eqs. (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.

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 κ * ( x). The upscaled quantities are a tensor quantity even though the input permeability, κ( x), 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 Fig. 1). This substantially reduces the degrees of freedom in the coarsegrid simulation. To compute the effective κ * for each subdomain, we impose different external forces G const 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 κ( x) on the fine grid. Then, the similar local LBM simulations usually need to be run with a constant tensor κ * as shown in Eq. (18). We seek κ * such that the average velocities from local fine-grid simulations with the heterogeneous κ( x) and homogeneous κ * , respectively, are equal (see Eqs. (20)-(21)). The onerous seeking process by adjusting the unknown κ * to match the fluxes computed using κ( x) is avoided in our simulations since κ * can be computed explicitly by Eq. (17).
We discuss two-dimensional problems as example. In the local LBM simulations using κ( x), we drive flow in the x direction by G (1) const = (G const , 0) and compute the average velocity u (1) κ( x) , where · is defined as a volume average over a subdomain. We also compute u (2) κ( x) by using G (2) const = (0, G const ) in another local simulation. Then, κ * is computed as follows: κ( x) · (0, 1) Now, we validate that the computed κ * satisfies the conservation principle of average fluxes. Assuming that we run local LBM simulations using the constant κ * computed by Eq. (17), Eq. (5) is modified to be: where κ * −1 is the inverse matrix of κ * . The evolution of f α ( x, t) is described by Eqs. (2), (13), (14), (15) and (18). As κ * and G are constant and the periodic boundary conditions are used in local simulations, the relation holds at steady state. For arbitrary ∆x, ∆t, τ , , ν, κ * and G = G const , the steady state solution of f α is independent of x and equal to: which implies that the uniform density is ρ = which implies that the average flux is conserved when using the same external force G const but different permeability distributions, namely using the heterogeneous κ( x) and homogeneous κ * , respectively. When driving flow by G (2) const = (0, G const ), the average flux is also conserved: After getting the value of κ * ( x) at each coarse point on the coarse grid, we implement two LBM simulations on the coarse and fine grids, respectively, inside the whole computational domain. The κ * ( x), ∆x coarse and ∆t coarse in the coarse-grid simulation are different from κ( x), ∆x fine and ∆t fine , respectively. The boundary conditions and the parameters ρ 0 , ν eff , , ν and G( x) in the coarse-grid simulation are the same as in the fine-grid simulation. In order to clearly verify the validity of the coarse-grid simulation of the whole computational domain, we use periodic boundary conditions to eliminate potential numerical errors which occur when using fixed pressures, for example, at the two ends along the x direction because fixed quantities are numerically imposed at the initial and last points along the x direction and their spatial positions are different when using different ∆x.

Comparison between the original and the proposed LBM algorithms
First, we verify the proposed simple LBM algorithm using the Shan-Chen force model against the original algorithm using the Guo el al. force model. In the two simulations using different force models, the number of grid points is 100×100 and ∆x = 0.01 m, ∆t = 0.0001 s and τ = 0.53 making ν eff = 0.01 m 2 s −1 , ν = 2 × 10 −6 m 2 s −1 , ρ 0 = 1000 kg m −3 , = 0.8 and G const = (2, 0) m s −2 . The periodic boundary conditions are used and the permeability assigned to each point with index (i, j) is: The average pressure over the whole computational domain is subtracted from the computed pressure p = c 2 s ρ in all figures of the pressure distributions. The transient results at the 5000 th ∆t and the steady state results at the 200000 th ∆t are given in Fig. 2 which shows the excellent agreement between the two simulations using different force models. The simulation using the Guo et al. force model takes about 23 minutes of computational time but the simulation using the Shan-Chen force model uses about 21 minutes. In the following LBM simulations, we only use the simple LBM algorithm with the Shan-Chen force model, which is described in Section 2.2.

Verifications of the computed κ *
We use the above setting of parameters but choose different values for τ , G const and κ for different problems in Section 4.2. We run simulations in the whole computational domain with prescribed distribution of κ( x) and verify the computed effective permeability κ * against the analytical solutions.

Layered distribution of κ
Here, we uniformly divide the whole domain into 10 layers parallel to the y axis. The odd number layers have κ 1 = 10 −12 m 2 and κ 2 in the even number layers is constant with values shown in Table 1 for different cases. Flow is driven in the x direction by a uniform G const = (2, 0) m s −2 . τ = 0.53 and so ν eff = 0.01 m 2 s −1 . The results in Table 1 show that the computed κ * xx by LBM agrees exactly with the analytical solution although the analytical formula is derived from the Darcy equation. This is because the steady state velocity is uniform and so the LBM simulations based on the Brinkman equation with nonzero ν eff actually yield the solutions of the Darcy equation at steady state.
When driving flow in the y direction by setting G const = (0, 2) m s −2 , the velocity distribution along the x direction is nonuniform. We set τ = 0.5 such that ν eff = 0 m 2 s −1 to recover the Darcy equation. As we can see in Table  2, the computed κ * yy by LBM simulations agrees exactly with the analytical solution.

Checkerboard distribution of κ
As on a checkerboard, we divide the whole computational domain uniformly into 10 × 10 squares with each square containing 10 × 10 points. The black squares of the checkerboard have κ 1 = 10 −12 m 2 and κ 2 in the white squares takes different values for different cases as shown in Table 3. Flow is driven by a uniform G const = (2, 0) m s −2 and we set ν eff = 0 m 2 s −1 to get the solution of the Darcy equation. The representative distributions of p, u 50.50000 × 10 −12 50.49999 × 10 −12 1000 500.5000 × 10 −12 500.4999 × 10 −12 10000 5000.500 × 10 −12 5000.499 × 10 −12 100000 50000.50 × 10 −12 50000.49 × 10 −12 and v are given in Fig. 3. The results in Table 3 show that the computed κ * xx by LBM simulations agrees well with the analytical solution when κ 2 κ 1 is not very large but deviates significantly in the case of high contrast. This deviation is due to the low spatial resolution of the grid used in the LBM simulations at high contrast of permeability. We refine the grid by increasing the total point number from 100 × 100 to 1000 × 1000 to show improving accuracy. ∆x and ∆t are changed to 10 −3 m and 10 −5 s, respectively. The results given in Table 3 show that the computed κ * xx becomes very close to the analytical solution when the permeability ratio is up to 100 but still significantly deviate from the correct value if the permeability ratio is very high, where more points are required to achieve good spatial resolution. Table 3: Verification of computed κ * xx , κ 1 = 10 −12 m 2 and ν eff = 0 m 2 s −1

. Simulations of Darcy flows
We choose a two-dimensional 1 m×1 m domain with periodic boundary conditions and = 0.8, ρ 0 = 1000 kg m −3 , ν = 2 × 10 −6 m 2 s −1 . We set ν eff = 0 m 2 s −1 by using τ = 0.5 in the simulations of both fine and coarse grids and also in the calculation of κ * ( x). In order to have obvious variations in the results of the coarse-grid simulation, a nonuniform external force G = (sin πx, sin πy) m s −2 is used and the distribution of permeability κ( x) in Fig. 4 is set according to Eq. (23) such that the distribution of κ * ( x) is nonuniform.
where κ const = 10 −13 m 2 . ∆x fine = 0.0025 m and ∆t fine = 0.000025 s in the fine-grid simulation. The number of fine points is 400×400 inside the whole computational domain which is divided uniformly into 40×40 subdomains. The averaged results over each set of 10 × 10 fine points located inside the same subdomain are computed and used to verify the results of the coarsegrid simulation. We have κ * xx = κ * yy = κ const and κ * yx = κ * xy = 0 inside the subdomains, where κ ≡ κ const . For subdomains with κ = 10(1+sin(80xπ) cos(80yπ))κ const , we use G const = (2, 0) m s −2 to drive flow and get κ * xx = 8.485κ const and κ * yx = 0. The symmetric property of κ( x) inside the subdomain implies that κ * yy = κ * xx and κ * xy = κ * yx . We define a scalar κ * as the average value over all diagonal components of κ * and for all subdomains we have κ * = κ * I, where I is the identity tensor. Now, Eq. (18), which is a general formula in the coarse-grid simulations, can be replaced by Eq. (5), where we change κ to κ * . In the case of κ * ≡ κ * I, the algorithm in the coarse-grid simulation is the same as in the fine-grid simulation (see Section 2.2) but they use different scalar permeability distributions, namely κ * and κ, respectively. In the coarse-grid simulation, ∆x coarse = 0.025 m and ∆t coarse = 0.00025 s. The number of coarse points is 40 × 40 and the value of κ * assigned to each coarse point with index (I, J) is: 29 ≤ I, J ≤ 32 κ * = κ const , 9 ≤ I ≤ 12, 29 ≤ J ≤ 32 κ * = κ const , 29 ≤ I ≤ 32, 9 ≤ J ≤ 12 κ * = 8.485κ const , otherwise. The distributions of the fine and coarse grids inside a representative area are given in Fig. 1. Figs. 5-6 show that the agreement is very good between the two simulations using the fine and coarse grids, respectively.

Simulations of Brinkman flows
The physical problem studied here is similar to that described in Section 4.3.1. The differences are that we increase the value of κ const to be κ const = 10 −7 m 2 (cf. κ const = 10 −13 m 2 in Section 4.3.1) and set ν eff = 10 −5 m 2 s −1 such that the contribution by the viscosity term ν eff ∆ u is remarkable as shown in Fig. 7, which shows the comparison between the averaged results of two fine-grid simulations using ν eff = 10 −5 and 0 m 2 s −1 , Figure 6: Detailed comparisons of p, u and v between the fine-grid averaged results and coarse-grid results using κ * , ν eff = 0 m 2 s −1 , G = (sin πx, sin πy) m s −2 , κ const = 10 −13 m 2 .

Conclusions
Pore-scale flows are routinely modeled by the LBM simulations due to their ability to handle complex geometries and physics. However, LBM simulations become very expensive as one uses large REVs. In this paper, we propose a upscaled LBM algorithm to model flows at coarse scales with a reduced computational complexity. The effective properties are computed by a local upscaling scheme. In this scheme, the local fine-grid simulations are performed and their results are averaged over the local region to compute effective properties. Effective properties are used in a coarse-grid LBM algorithm to perform the simulations at larger scales. The coarse-grid LBM simulation using the computed effective permeability agrees very well with the fine-grid LBM simulation. In addition, simulation results show that the coarse-grid LBM simulation will deviate significantly from the fine-grid LBM simulation if the effective permeability is computed by neglecting the viscosity term in modeling Brinkman flows.
Although the results presented in this paper are encouraging, there is scope for further exploration of some of the underlying approaches. As our intent here was to demonstrate that coarse scale information could be effectively used to design upscaled LBM representations, we did not consider Figure 8: Comparisons of p, u and v between the fine-grid averaged results (left), coarse-grid results using κ * (κ, ν eff ) (middle) and κ * ,err (κ) (right), ν eff = 10 −5 m 2 s −1 , G = (sin πx, sin πy) m s −2 , κ const = 10 −7 m 2 . Figure 9: Detailed comparisons of p, u and v between the fine-grid averaged results and coarse-grid results using κ * (κ, ν eff ), ν eff = 10 −5 m 2 s −1 , G = (sin πx, sin πy) m s −2 , κ const = 10 −7 m 2 . challenging heterogeneous cases with high-contrast permeability. It is known (e.g., [13] and [17]) that the presence of high heterogeneities, such as channels and high contrast, will cause a decrease in the accuracy of upscaling methods for Darcy flow problems. Similarly, we expect that our upscaled LBM algorithm will require an additional treatment to handle highly heterogeneous cases. These treatments can include oversampling, local-global, or global techniques or possibly upscaled techniques. Some of these treatments can be easily incorporated into our new upscaled LBM framework.