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

1. Introduction

Detailed flow simulations in porous media are often modeled using the Darcy or Brinkappaman 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 up-scaled 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 c = [increment of x]/[DELTA]t, where [increment of x] is the Lattice spacing and [DELTA]t is the time step. For two-dimensional problems, we use the D2Q9 Lattice model [4], where Q = 9 is the total number of Lattice velocities, [[??].sub.0] = (0,0), [[omega].sub.0] = 4/9, [[??].sub.[alpha]] = ([cos[theta].sub.[alpha]],[sin[theta].sub.[alpha]]c, [[theta].sub.[alpha]] = ([alpha] - 1)[pi]/2, [[omega].sub.[alpha]] = 1/9 for [alpha] = 1 to 4, [[??].sub.[alpha]] = ([cos.[theta].sub.[alpha]], [sin.[theta].sub.[alpha]])[square root of 2c], [[theta].sub.[alpha]] = ([alpha] - 5)[pi]/2 + [pi]/4, and [[omega].sub.[alpha]] = 1/36 for [alpha] = 5 to 8. For three-dimensional problems, the D3Q19 Lattice model [4] with different [[??].sub.[alpha]] and [[omega].sub.[alpha]] 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 [f.sub.[alpha]]([??], t) is

[mathematical expression not reproducible] (1)

where the normalized relaxation time [tau] is appropriately selected to match the desired effective kinematic viscosity [v.sub.eff] = [c.sup.2.sub.s] ([tau] - 0.5)[DELTA]t, where [c.sub.s] = c/[square root of 3] is the sound speed. We use the simplified truncated form of the equilibrium distribution function as

[mathematical expression not reproducible] (2)

where the density is given by [rho]([??], t) = [summation.sup.Q-1.sub.[alpha]=0 [f.sub.[alpha]] ([??], t) and the equilibrium velocity, [[??].sup.(eq)]([??], t), is defined as

[mathematical expression not reproducible] (3)

where [mathematical expression not reproducible] is the force per unit mass. Similarly, [F.sub.[alpha]]([??], t) is simplified to

[mathematical expression not reproducible] (4)

In the force model proposed by Guo et al. [12], the flow velocity [??] is equal to [[??].sup.(eq)]. If [[??].sub.m] is constant, p and [[??].sup.(eq)] are computed by [f.sub.[alpha]] and then [f.sup.(eq).sub.[alpha]] and [F.sub.[alpha]] of (1) are determined explicitly. For solving the Darcy-scale equation, we consider the following expressions for [[??].sub.m] as a linear function of u (see [10]):

[mathematical expression not reproducible] (5)

where [epsilon] is the porosity, v is the physical kinematic viscosity of the fluid, [kappa]([??]) is a scalar for the permeability, and [mathematical expression not reproducible] is the external body force per unit mass. The force introduced above incorporates the porous media heterogeneities through the permeability function [kappa]([??]) and depends on the microstructure. If [kappa]([??]) has a high value in the region, then one can assume that this region is highly permeable, while if [kappa]([??]) 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 [??] = [[??].sup.(eq)]. Using this fact and solving for [??] in (5), we obtain the explicit formula as in [10]:

[mathematical expression not reproducible]. (6)

[[??].sub.m] is computed by (5). Then, [f.sup.(eq).sub.[alpha]] and [F.sub.[alpha]] of (1) are determined by (2) and (4), respectively.

In the incompressible limit with [absolute value of ([??])] [much less than] [c.sub.s], the analysis [10] based on the Chapman-Enskog expansion shows that the computed pressure p = [c.sup.2.sub.s][rho] and flow velocity [??] converge to the solutions of the following equation:

[mathematical expression not reproducible] (7)

where [[rho].sub.0] is the initial mass density used in LBM simulations. Here, [[rho].sub.0] needs not be the real density [[rho].sub.real] of the incompressible fluid; then, the computed p[[rho].sub.real]/[[rho].sub.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 [v.sub.eff], [epsilon], v, and [kappa]([??]) 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 oth component having the velocity [[??].sub.[alpha]] at [??] and time t is denoted by [f.sup.o.sub.[alpha]] ([??], t), where o = 1, ..., S and S is the total number of components. The general updating algorithm of [f.sup.o.sub.[alpha]] ([??], t) is

[mathematical expression not reproducible] (8)

where [sigma] = 1, ..., S and the equilibrium distribution function is defined as

[mathematical expression not reproducible] (9)

where [sigma] = 1, ..., S, [[rho].sup.[sigma]] = [summation.sup.Q-1.sub.[alpha]=0][f.sup.[sigma].sub.[alpha]], and [[??].sup.[sigma](eq)] is computed as

[mathematical expression not reproducible] (10)

where [mathematical expression not reproducible] is related to the total volume force acting on oth component. Generally speaking [16], [[??].sup.o] contains three parts: the fluid-fluid interaction [[??].sup.1,[sigma]], fluid-solid interaction [[??].sup.2,o], and external force [[??].sup.3,[sigma]]. For example, [[??].sup.3,[sigma]] = [DELTA]t[[rho].sup.[sigma]][??] for the contribution by the external body force [??] per unit mass. In (10), [??]' is defined as follows to conserve momentum:

[mathematical expression not reproducible] (11)

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:

[mathematical expression not reproducible] (12)

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 [f.sub.[alpha]]([??], t) without notation [sigma] is

[mathematical expression not reproducible] (13)

In order to recover the Brinkman equation, the equilibrium distribution function [f.sup.(eq).sub.[alpha]] of (9) is simplified to (2). [rho] = [summation.sup.Q-1.sub.[alpha]=0] [f.sub.[alpha]] but [[??].sup.(eq)] of (10) is modified as follows:

[mathematical expression not reproducible] (14)

where we use [DELTA]t[rho][[??].sub.m] to replace the original notation, which is equal to the momentum increase per unit volume after [DELTA]t due to the force effect through the relaxation process. Correspondingly, the flow velocity [??] of (12) is modified to

[mathematical expression not reproducible] (15)

When solving the Brinkman equation, [[??].sub.m] 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, [[??].sup.(eq)] is computed as

[mathematical expression not reproducible] (16)

which is obtained by solving (14) and (15).

As we can see, the computations of [[??].sub.m] by (5) and [F.sub.[alpha]] by (4) using the computed [[??].sub.m] in the original algorithm [10] are avoided in the current simplified algorithm and, therefore, the efficiency is improved. In the incompressible limit with [absolute value of ([??])] << [c.sub.s], the computed pressure p = [c.sup.2.sub.s][rho] 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 [[kappa].sup.*]([??]). The upscaled quantities are a tensor quantity even though the input permeability, [kappa]([??]), 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 [[kappa].sup.*] for each subdomain, we impose different external forces [[??].sub.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 [kappa]([??]) on the fine grid. Then, the similar local LBM simulations usually need to be run with a constant tensor [[kappa].sup.*] as shown in (18). We seek [[kappa].sup.*] such that the average velocities from local fine-grid simulations with the heterogeneous [kappa](x) and homogeneous [[kappa].sup.*], respectively, are equal (see (21)-(22)). The onerous seeking process by adjusting the unknown [[kappa].sup.*] to match the fluxes computed using [kappa]([??]) is avoided in our simulations since [[kappa].sup.*] can be computed explicitly by (17).

We discuss two-dimensional problems as example. In the local LBM simulations using [kappa]([??]), we drive flow in x direction by [[??].sup.(1).sub.const] = ([G.sub.const], 0) and compute the average velocity [mathematical expression not reproducible] where [bar.*] is defined as a volume average over a subdomain. We also compute [mathematical expression not reproducible] by using [[??].sup.(2).sub.const] = (0, [G.sub.const]) in another local simulation. Then, [[kappa].sup.*] is computed as follows:

[mathematical expression not reproducible] (17)

Now, we validate that the computed [[kappa].sup.*] satisfies the conservation principle of average fluxes. Assuming that we run local LBM simulations using the constant [[kappa].sup.*] computed by (17), (5) is modified to be

[mathematical expression not reproducible] (18)

where [[kappa].sup.*-1] is the inverse matrix of [[kappa].sup.*]. The evolution of [f.sub.[alpha]]([??], t) is described by (2), (13), (14), (15), and (18). As [[kappa].sup.*] and [??] are constant and the periodic boundary conditions are used in local simulations, the relation

[mathematical expression not reproducible] (19)

holds at steady state. For arbitrary [increment of x], [DELTA]t, [tau], [epsilon], v, [[kappa].sup.*] , and [??] = [[??].sub.const], the steady state solution of [f.sub.[alpha]] is independent of [??] and equal to

[mathematical expression not reproducible] (20)

which implies that the uniform density is [rho] = [[summation].sup.Q-1.sub.[alpha]=0][f.sub.[alpha]] = [[rho].sub.o).] We validate the solution of (20) by the following verification: substituting (20) into (15) and considering (18), we get the uniform velocity [mathematical expression not reproducible]. In addition, we get [[??].sub.m] [equivalent to] 0 by (18) using [??]. Then, substituting [f.sub.[alpha]] and [[??].sub.m] into (14), we get [[??].sup.(eq)] = [??], which implies that (13) is satisfied at steady state since we have [f.sub.[alpha]] = [f.sup.(eq).sub.[alpha]] according to (2) and [mathematical expression not reproducible]. According to the uniform solution of [mathematical expression not reproducible] and (17), the average velocity [mathematical expression not reproducible], of the local simulation using constant [[kappa].sup.*] and [[??].sup.(1).sub.const] = ([G.sub.const], 0) satisfies

[mathematical expression not reproducible] (21)

which implies that the average flux is conserved when using the same external force [[??].sup.(1).sub.const] but different permeability distributions, namely, using the heterogeneous [kappa]([??]) and homogeneous [[kappa].sup.*], respectively. When driving flow by [[??].sup.(2).sub.const] = (0, [G.sub.const]), the average flux is also conserved:

[mathematical expression not reproducible] (22)

After getting the value of [[kappa].sup.*] ([??]) 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. [[kappa].sup.*] ([??]), [DELTA][x.sub.coarse], and [DELTA][t.sub.coarse] in the coarse-grid simulation are different from [kappa]([??]), [DELTA][x.sub.fine], and [DELTA][t.sub.fine], respectively. The boundary conditions and the parameters [[rho].sub.0], [v.sub.eff], [epsilon], v, and [mathematical expression not reproducible] 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 that occur when using fixed pressures, for example, at the two ends along x direction because fixed quantities are numerically imposed at the initial and last points along x direction and their spatial positions are different when using different [increment of x].

4. Numerical Results

4.1. 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 et al. force model. In the two simulations using different force models, the number of grid points is 100 x 100 and [increment of x] = 0.01 m, [DELTA]t = 0.0001s, and [tau] = 0.53 making [v.sub.eff] = 0.01 [m.sup.2] [s.sup.-1], v = 2 x [10.sup.-6] [m.sup.2] [s.sup.-1], [[rho].sub.0] = 1000 [kappa][m.sup.-3], [epsilon] = 0.8, and [[??].sub.const] = (2,0) m [s.sup.-2]. The periodic boundary conditions are used and the permeability assigned to each point with index (i, j) is

[mathematical expression not reproducible] (23)

The average pressure over the whole computational domain is subtracted from the computed pressure p = [c.sup.2.sub.s][rho] in all figures of the pressure distributions. The transient results at 5000th[DELTA]t and the steady state results at 200000th[DELTA]t are given in Figure 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.

4.2. Verifications of the Computed [[kappa].sup.*]. We use the above setting of parameters but choose different values for [tau], [[??].sub.const], and [kappa] for different problems in Section 4.2. We run simulations in the whole computational domain with prescribed distribution of [kappa]([??]) and verify the computed effective permeability [[kappa].sup.*] against the analytical solutions.

4.2.1. Layered Distribution of [kappa]. Here, we uniformly divide the whole domain into 10 layers parallel to the y-axis. The odd number layers have [[kappa].sub.1] = [10.sup.-12] [m.sup.2] and [[kappa].sub.2] in the even number layers is constant with values shown in Table 1 for different cases. Flow is driven in x direction by a uniform [[??].sub.const] = (2,0) m[s.sup.-2]. [tau] = 0.53 and so [v.sub.eff] = 0.01 [m.sup.2] [s.sup.-1]. The results in Table 1 show that the computed 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 [v.sub.eff] actually yield the solutions of the Darcy equation at steady state.

When driving flow in y direction by setting [[??].sub.const] = (0,2) m[s.sup.-2], the velocity distribution along x direction is nonuniform. We set [tau] = 0.5 such that [v.sub.eff] = 0 [m.sup.2] [s.sup.-1] to recover the Darcy equation. As we can see in Table 2, the computed [[kappa].sup.*.sub.yy] by LBM simulations agrees exactly with the analytical solution.

4.2.2. Checkerboard Distribution of kappa. As on a checkerboard, we divide the whole computational domain uniformly into 10 x 10 squares with each square containing 10 x 10 points. The black squares of the checkerboard have [[kappa].sup.1] = [10.sup.-12] [m.sup.2] and [kappa.sup.2] in the white squares takes different values for different cases as shown in Table 3. Flow is driven by a uniform [[??].sub.const] = (2, 0) m[s.sup.-2] and we set [v.sup.eff] = 0 [m.sup.2] [s.sup.-1] to get the solution of the Darcy equation. The representative distributions of p, [micro], and v are given in Figure 3. The results in Table 3 show that the computed [[kappa].sup.*.sub.xx] by LBM simulations agrees well with the analytical solution when [kappa.sub.2]/[[kappa].sub.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 x 100 to 1000 x 1000 to show improving accuracy. [DELTA]x and [DELTA]t are changed to [10.sup.-3] m and [10.sup.-5] s, respectively. The results given in Table 3 show that the computed [[kappa].sup.*.sub.xx] becomes very close to the analytical solution when the permeability ratio is up to 100 but still significantly deviates from the correct value if the permeability ratio is very high, where more points are required to achieve good spatial resolution.

4.3. Verifications of the Upscaled Simulation Scheme

4.3.1. Simulations of Darcy Flows. We choose a two-dimensional 1 m x 1 m domain with periodic boundary conditions and [mathematical expression not reproducible], and [mathematical expression not reproducible] We set [v.sub.eff] = 0 [m.sup.2] [s.sup.-1] by using [tau] = 0.5 in the simulations of both fine and coarse grids and also in the calculation of [[kappa].sup.*] ([??]). In order to have obvious variations in the results of the coarse-grid simulation, a nonuniform external force [mathematical expression not reproducible] is used and the distribution of permeability [kappa]([??])in Figure 4 is set according to (24) such that the distribution of [[kappa].sup.*] ([??]) is nonuniform.

[mathematical expression not reproducible] (24)

where [[kappa].sub.const] = [10.sup.-13] [m.sup.2.] [DELTA][x.sub.fine] = 0.0025 m and [DELTA][t.sub.fine] = 0.000025 s in the fine-grid simulation. The number of fine points is 400 x 400 inside the whole computational domain that is divided uniformly into 40 x 40 subdomains. The averaged results over each set of 10 x 10 fine points located inside the same subdomain are computed and used to verify the results of the coarse-grid simulation.

We have [mathematical expression not reproducible] inside the subdomains, where kappa = [[kappa].sub.const]. For subdomains with kappa = [mathematical expression not reproducible], we use [[??].sub.const] = (2,0) [ms.sup.-2] to drive flow and get [[kappa].sup.*.sub.xx] = 8.485[[kappa].sub.const] and [[kappa].sup.*.sub.yx] = 0. The symmetric property of [kappa]([??]) inside the subdomain implies that [mathematical expression not reproducible]. We define a scalar [[kappa].sup.*] as the average value over all diagonal components of [[kappa].sup.*] and for all subdomains we have [[kappa].sup.*] = [[kappa].sup.*] I, where I is the identity tensor. Now, (18), which is a general formula in the coarse-grid simulations, can be replaced by (5), where we change kappa to [[kappa].sup.*]. In the case of [[kappa].sup.*] = [[kappa].sup.*] 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 [[kappa].sup.*] and kappa, respectively. In the coarse-grid simulation, [DELTA][x.sub.coarse] = 0.025 m and [DELTA][t.sub.coarse] = 0.00025 s. The number of coarse points is 40 x 40 and the value of [[kappa].sup.*] assigned to each coarse point with index (7, 7) is

[mathematical expression not reproducible] (25)

The distributions of the fine and coarse grids inside a representative area are given in Figure 1. Figures 5-6 show that the agreement is very good between the two simulations using the fine and coarse grids, respectively.

4.3.2. 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 [[kappa].sub.const] to be [[kappa].sub.const] = [10.sup.-7] [m.sup.2] (cf. [[kappa].sub.const] = [10.sup.-13] [m.sup.2] in Section 4.3.1) and set [v.sub.eff] = [10.sup.-5] [m.sup.2] [s.sup.-1] such that the contribution by the viscosity term [v.sub.eff] [DELTA][??] is large enough to be noticed as shown in Figure 7, which shows the comparison between the averaged results of two fine-grid simulations using [v.sub.eff] = [10.sup.-5] and 0 [m.sup.2] [s.sup.-1], respectively. In this regime, flow in some regions is close to Stokes flow while, in other regions, flow is close to Darcy flow. Since [v.sub.eff] is nonzero here, we let [v.sub.eff] = [10.sup.-5] [m.sup.2] [s.sup.-1] when computing [[kappa].sup.*] and get [[kappa].sup.*.sub.xx] (kappa, [v.sub.eff]) = [[kappa].sup.*.sub.yy](kappa, [v.sub.eff]) = 7..367[[kappa].sub.const] and [[kappa].sup.*.sub.yx](kappa, [v.sub.eff]) = [[kappa].sup.*.sub.xy](kappa, [v.sub.eff]) = 0 for subdomains with [mathematical expression not reproducible]. For subdomains with [mathematical expression not reproducible]. Thus, we have [[kappa].sup.*] (kappa, [v.sub.eff]) equal to 7.367[[kappa].sub.const]I or [[kappa].sub.const]I. As discussed in Section 4.3.1, we can use a scalar distribution of [[kappa].sup.*] (kappa, [v.sub.eff]), which is equal to 7.367[[kappa].sub.const] or [[kappa].sub.const], in the coarse-grid simulation. We use [DELTA][x.sub.fine] = 0.0025 m, [DELTA][t.sub.fine] = 0.000025 s, and [[tau].sub.fine] = 0.50012 in the fine-grid simulation and use [DELTA][x.sub.coarse] = 0.025 m, [DELTA][t.sub.coarse] = 0.00025 s, and [[tau].sub.course] = 0.500012 in the coarse-grid simulation. Figures 8-9 show that the agreement of the coarse-grid simulation using [[kappa].sup.*] (kappa, [v.sub.eff]) with the fine-grid simulation is very good. In addition, we set [v.sub.eff] = 0 [m.sup.2] [s.sup.-1] when computing [[kappa].sup.*,err] and get [mathematical expression not reproducible] for subdomains with [mathematical expression not reproducible]. For subdomains with [mathematical expression not reproducible], where we still use the superscript "err" since the local simulation procedure with [v.sub.eff] = 0 [m.sup.2] [s.sup.-1] is wrong although the obtained [[kappa].sup.*,err](kappa) is the same as the above [[kappa].sup.*] (kappa, [v.sub.eff]). Now, we have [[kappa].sup.*,err](kappa) equal to 8.485[[kappa].sub.const]I or [[kappa].sub.const]I. We use a scalar distribution of [[kappa].sup.*,err](kappa), which is equal to 8.485[[kappa].sub.const] or [[kappa].sub.const], in another coarse-grid simulation. Note that the difference between [[kappa].sup.*,err](kappa) and [[kappa].sup.*] (kappa, [v.sub.eff]) is distinct in subdomains with nonuniform [kappa]([??]). The results of the coarse-grid simulation using [[kappa].sup.*,err]([kappa]) are also given in Figure 8, which shows that the deviation of the coarse-grid simulation using [[kappa].sup.*,err]([kappa]) from the fine-grid simulation is noticeable. Thus, the previous computation of [[kappa].sup.*] (kappa, [v.sub.eff]) as the effective permeability using nonzero [v.sub.eff] is accurate for upscaling the Brinkman equation.

5. 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 an upscaled LBM algorithm to model flows at the Darcy scale using coarse grids 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 using the original permeability distribution. 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 challenging heterogeneous cases with high-contrast permeability. It is known (e.g., [13, 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. In real reservoir or groundwater applications, the boundary condition with fixed pressure is more realistic than the periodic boundary condition used here and then additional treatment is also needed to accurately compute the effective permeability for the subdomains (i.e., coarse grids) close to the boundaries with fixed pressures.

http://dx.doi.org/10.1155/2017/1740693

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to acknowledge Victor Calo for his helpful discussion on the physical implications of the models. Also, they would like to thank Oleg Iliev for his helpful insights on the upscaling of porous media and validation of the computational results.

References

[1] F. Khan, F. Enzmann, M. Kersten, A. Wiegmann, and K. Steiner, "3D simulation of the permeability tensor in a soil aggregate on basis of nanotomographic imaging and LBE solver," Journal of Soils and Sediments, vol. 12, no. 1, pp. 86-96, 2012.

[2] G. R. McNamara and G. Zanetti, "Use of the boltzmann equation to simulate lattice-gas automata," Physical Review Letters, vol. 61, no. 20, pp. 2332-2335, 1988.

[3] H. Chen, S. Chen, and W. H. Matthaeus, "Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method," Physical Review A, vol. 45, no. 8, pp. R5339-R5342, 1992.

[4] Y. H. Qian, D. D'Humieres, and P. Lallemand, "Lattice BGKAPPA models for Navier-Stokes equation," Europhysics Letters, vol. 17, no. 6, pp. 479-484, 1992.

[5] X. Shan and H. Chen, "Lattice Boltzmann model for simulating flows with multiple phases and components," Physical Review E, vol. 47, no. 3, pp. 1815-1819, 1993.

[6] X. Shan and H. Chen, "Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation," Physical Review E, vol. 49, no. 4, pp. 2941-2948, 1994.

[7] X. Shan and G. Doolen, "Diffusion in a multicomponent lattice Boltzmann equation model," Physical Review E, vol. 54, no. 4, pp. 3614-3620, 1996.

[8] J. Zhu and J. Ma, "An improved gray lattice Boltzmann model for simulating fluid flow in multi-scale porous media," Advances in Water Resources, vol. 56, pp. 61-76, 2013.

[9] Q. Kang, D. Zhang, and S. Chen, "Unified lattice Boltzmann method for flow in multiscale porous media," Physical Review E--Statistical, Nonlinear, and Soft Matter Physics, vol. 66, no. 5, Article ID 056307, 2002.

[10] Z. Guo and T. S. Zhao, "Lattice Boltzmann model for incompressible flows through porous media," Physical Review E, vol. 66, no. 3, Article ID 036304, 2002.

[11] P. Nithiarasu, K. N. Seetharamu, and T. Sundararajan, "Natural convective heat transfer in a fluid saturated variable porosity medium," International Journal of Heat and Mass Transfer, vol. 40, no. 16, pp. 3955-3967, 1997.

[12] Z. Guo, C. Zheng, and B. Shi, "Discrete lattice effects on the forcing term in the Lattice Boltzmann method," Physical Review E, vol. 65, no. 4, Article ID 046308, 2002.

[13] Y. Efendiev and T. Y. Hou, Multiscale Finite Element Methods. Theory and Applications, vol. 4 of Surveys and Tutorials in the Applied Mathematical Sciences, Springer, New York, NY, USA, 2009.

[14] M. C. Sukop and D. T. Thorne Jr., Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers, Springer, 2006.

[15] J. Li, D. Brown, V. Calo, Y. Efendiev, and O. Iliev, "Multiscale Lattice Boltzmann method for flow simulations in highly heterogenous porous media," in Proceedings of the SPE Reservoir Characterization and Simulation Conference and Exhibition, SPE-165985-MS, Abu Dhabi, UAE, September 2013.

[16] Q. Kang, D. Zhang, and S. Chen, "Displacement of a two-dimensional immiscible droplet in a channel," Physics of Fluids, vol. 14, no. 9, pp. 3203-3214, 2002.

[17] G. Qin, L. Bi, P. Popov, Y. Efendiev, and S. M. Espedal, "An efficient upscaling procedure based on Stokes-Brinkman model and discrete fracture network method for naturally fractured carbonate karst reservoirs," in Proceedings of the International Oil and Gas Conference and Exhibition in China, SPE-132236MS, June 2010.

Jun Li (1, 2) and Donald Brown (1,3)

(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, UKAPPA

Correspondence should be addressed to Jun Li; lijun04@gmail.com

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

Caption: FIGURE 1: Schematic models of the fine and coarse grids.

Caption: FIGURE 2: Comparisons of p, u, and v between two LBM simulations using different force models, transient results at 5000th [DELTA]t (a) and the steady state results (b), [v.sub.eff] = 0.01 [m.sup.2] [s.sup.-1] and [[??].sub.const] = (2,0) m [s.sup.-2].

Caption: FIGURE 3: Distributions of p, w, and v, [v.sub.eff] = 0 [m.sup.2] [s.sup.-1], [[??].sub.const] = (2,0) m [s.sup.-2], [[kappa].sub.1] = [10.sup.-12] [m.sup.2], and [kappa.sub.2]/[[kappa].sub.1] = 2.

Caption: FIGURE 4: Distribution of the permeability [mathematical expression not reproducible].

Caption: FIGURE 5: Comparisons of p, u, and v between the fine-grid results (a), fine-grid averaged results (b), and coarse-grid results using [[kappa].sup.*] (c), [mathematical expression not reproducible].

Caption: FIGURE 6: Detailed comparisons of p, u, and v between the fine-grid averaged results and coarse-grid results using [mathematical expression not reproducible]

Caption: FIGURE 7: Comparisons of p, u, and v between two fine-grid simulations using [mathematical expression not reproducible].

Caption: FIGURE 8: Comparisons of p, u, and v between the fine-grid averaged results (a) and coarse-grid results using [mathematical expression not reproducible].

Caption: FIGURE 9: Detailed comparisons of p, u, and v between the fine-grid averaged results and coarse-grid results using [mathematical expression not reproducible]
```
TABLE 1: Verification of computed [[kappa].sup.*.sub.xx],
[[kappa].sub.1] = [10.sup.-12] [m.sup.2] and [v.sub.eff] = 0.01
[m.sup.2] [s.sup.-1].

[[kappa].sub.2]/            [(1/2)             [[kappa].sup.*.sub.xx]
[[kappa].sub.1]     (1/[[kappa].sub.1] + 1/            by LBM
[[kappa].sub.2])].sup.-1]

2                   1.33333 x [10.sup.-12]     1.33333 x [10.sup.-12]
10                  1.81818 x [10.sup.-12]     1.81818 x [10.sup.-12]
50                  1.96078 x [10.sup.-12]     1.96078 x [10.sup.-12]
100                 1.98019 x [10.sup.-12]     1.98019 x [10.sup.-12]
1000                1.99800 x [10.sup.-12]     1.99800 x [10.sup.-12]
10000               1.99980 x [10.sup.-12]     1.99979 x [10.sup.-12]
100000              1.99998 x [10.sup.-12]     1.99998 x [10.sup.-12]

TABLE 2: Verification of computed [[kappa].sup.*.sub.yy],
[[kappa].sub.1] = [10.sup.-12] [m.sup.2] and [v.sub.eff] = 0
[m.sup.2] [s.sup.?1].

[[kappa].sub.2]/            (1/2)            [[kappa].sup.*.sub.yy]
[[kappa].sub.1]       ([[kappa].sub.1]               by LBM
[[kappa].sub.2])

2                  1.500000 x [10.sup.-12]   1.499999 x [10.sup.-12]
10                 5.500000 x [10.sup.-12]   5.499999 x [10.sup.-12]
50                 25.50000 x [10.sup.-12]   25.49999 x [10.sup.-12]
100                50.50000 x [10.sup.-12]   50.49999 x [10.sup.-12]
1000               500.5000 x [10.sup.-12]   500.4999 x [10.sup.-12]
10000              5000.500 x [10.sup.-12]   5000.499 x [10.sup.-12]
100000             50000.50 x [10.sup.-12]   50000.49 x [10.sup.-12]

TABLE 3: Verification of computed [k.sup.*.sub.xx],
[[kappa].sub.1] = [10.sup.-12] [m.sup.2] and
[v.sub.eff] = 0 [m.sup.2 [s.sup.-1].

[[kappa].sub.2]/      [square root of       [[kappa].sup.*.sub.xx]
[[kappa].sub.1]       ([[kappa].sub.1]              by LBM
[[kappa].sub.2])]

2                  1.41421 x [10.sup.-12]   1.41418 x [10.sup.-12]
10                 3.16227 x [10.sup.-12]   3.14081 x [10.sup.-12]
50                 7.07106 x [10.sup.-12]   6.45938 x [10.sup.-12]
100                10.0000 x [10.sup.-12]   8.25393 x [10.sup.-12]
1000               31.6227 x [10.sup.-12]   12.2496 x [10.sup.-12]
10000              100.000 x [10.sup.-12]   13.0133 x [10.sup.-12]

[[kappa].sub.2]/   [[kappa].sup.*.sub.xx]
[[kappa].sub.1]            by LBM

(1000 x 1000 points)
2
10
50                 7.01357 x [10.sup.-12]
100                9.70489 x [10.sup.-12]
1000               19.8897 x [10.sup.-12]
10000              23.2777 x [10.sup.-12]
```