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

Authors: Jun Li; Donald Brown
Journal: Geofluids
 (2017)

---

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


---

## Body

## 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 andc=Δx/Δt, where Δx is the Lattice spacing and Δ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, e→0=(0,0), ω0=4/9, e→α=(cos⁡θα,sin⁡θα)c, θα=(α-1)π/2, ωα=1/9 for α= 1 to 4, e→α=(cos⁡θα,sin⁡θα)2c, θα=(α-5)π/2+π/4, and ωα=1/36 for α= 5 to 8. For three-dimensional problems, the D3Q19 Lattice model [4] with different e→α 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 fα(x→,t) is(1)fαx→+Δte→α,t+Δt=fαx→,t+fαeqx→,t-fαx→,tτ+ΔtFαx→,t,where the normalized relaxation time τ is appropriately selected to match the desired effective kinematic viscosity νeff=cs2(τ-0.5)Δt, where cs=c/3 is the sound speed. We use the simplified truncated form of the equilibrium distribution function as(2)fαeq=ωαρ1+e→α·u→eqcs2,where the density is given by ρ(x→,t)=∑α=0Q-1fα(x→,t) and the equilibrium velocity, u→(eq)(x→,t), is defined as(3)u→eq=∑α=0Q-1e→αfα+1/2Δtρf→mρ,where f→m(x→,t) is the force per unit mass. Similarly, Fα(x→,t) is simplified to(4)Fα=ωαρ1-12τe→α·f→mcs2.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 (1) are determined explicitly. For solving the Darcy-scale equation, we consider the following expressions for f→m as a linear function of u→ (see [10]):(5)f→m=-ϵνκu→+ϵG→,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 that are discussed in [10, 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]:(6)u→=∑α=0Q-1e→αfα+Δt/2ϵρG→ρ1+ϵΔtν/2κ.f→m is computed by (5). Then, fα(eq) and Fα of (1) are determined by (2) and (4), respectively.In the incompressible limit with|u→|≪cs, the analysis [10] based on the Chapman-Enskog expansion shows that the computed pressure p=cs2ρ and flow velocity u→ converge to the solutions of the following equation:(7)∇·u→=0∂u→∂t=-1ρ0∇p+νeffΔu→-ϵνκu→+ϵG→,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.


### 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 σ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(8)fασx→+Δte→α,t+Δt=fασx→,t+fασeqx→,t-fασx→,tτσ,where σ=1,…,S and the equilibrium distribution function is defined as(9)fασeq=ρσωα1+e→α·u→σeqcs2+e→α·u→σeq22cs4-u→σeq·u→σeq2cs2,where σ=1,…,S,ρσ=∑α=0Q-1fασ, and u→σ(eq) is computed as(10)u→σeq=ρσu→′+τσF→σρσ,σ=1,…,S,where F→σ(x→,t) is related to the total volume force acting on σ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 (10), u→′ is defined as follows to conserve momentum:(11)u→′=∑σ=1S1/τσ∑α=0Q-1e→αfασ∑σ=1S1/τσ∑α=0Q-1fασ.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:(12)u→=∑σ=1S∑α=0Q-1e→αfασ+1/2∑σ=1SF→σ∑σ=1Sρσ.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 offα(x→,t) without notation σ is(13)fαx→+Δte→α,t+Δt=fαx→,t+fαeqx→,t-fαx→,tτ.In order to recover the Brinkman equation, the equilibrium distribution function fα(eq) of (9) is simplified to (2). ρ=∑α=0Q-1fα but u→(eq) of (10) is modified as follows:(14)u→eq=∑α=0Q-1e→αfα+τΔtρf→mρ,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 (12) is modified to(15)u→=∑α=0Q-1e→αfα+1/2Δtρf→mρ.When solving the Brinkman equation, f→m is a function of u→ defined by (5). A comparison of (3) and (15) shows that the explicit formula of (6) to compute u→ is also valid here. Then, u→(eq) is computed as(16)u→eq=2τu→+1-2τ∑α=0Q-1e→αfαρ,which is obtained by solving (14) and (15).As we can see, the computations off→m by (5) and Fα by (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→≪cs, the computed pressure p=cs2ρ and flow velocity u→ 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.


## 2.1. LBM Algorithm with the Guo Force Model

Following the algorithm presented in [10], the evolution algorithm of the distribution function fα(x→,t) is(1)fαx→+Δte→α,t+Δt=fαx→,t+fαeqx→,t-fαx→,tτ+ΔtFαx→,t,where the normalized relaxation time τ is appropriately selected to match the desired effective kinematic viscosity νeff=cs2(τ-0.5)Δt, where cs=c/3 is the sound speed. We use the simplified truncated form of the equilibrium distribution function as(2)fαeq=ωαρ1+e→α·u→eqcs2,where the density is given by ρ(x→,t)=∑α=0Q-1fα(x→,t) and the equilibrium velocity, u→(eq)(x→,t), is defined as(3)u→eq=∑α=0Q-1e→αfα+1/2Δtρf→mρ,where f→m(x→,t) is the force per unit mass. Similarly, Fα(x→,t) is simplified to(4)Fα=ωαρ1-12τe→α·f→mcs2.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 (1) are determined explicitly. For solving the Darcy-scale equation, we consider the following expressions for f→m as a linear function of u→ (see [10]):(5)f→m=-ϵνκu→+ϵG→,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 that are discussed in [10, 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]:(6)u→=∑α=0Q-1e→αfα+Δt/2ϵρG→ρ1+ϵΔtν/2κ.f→m is computed by (5). Then, fα(eq) and Fα of (1) are determined by (2) and (4), respectively.In the incompressible limit with|u→|≪cs, the analysis [10] based on the Chapman-Enskog expansion shows that the computed pressure p=cs2ρ and flow velocity u→ converge to the solutions of the following equation:(7)∇·u→=0∂u→∂t=-1ρ0∇p+νeffΔu→-ϵνκu→+ϵG→,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.


## 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 σ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(8)fασx→+Δte→α,t+Δt=fασx→,t+fασeqx→,t-fασx→,tτσ,where σ=1,…,S and the equilibrium distribution function is defined as(9)fασeq=ρσωα1+e→α·u→σeqcs2+e→α·u→σeq22cs4-u→σeq·u→σeq2cs2,where σ=1,…,S,ρσ=∑α=0Q-1fασ, and u→σ(eq) is computed as(10)u→σeq=ρσu→′+τσF→σρσ,σ=1,…,S,where F→σ(x→,t) is related to the total volume force acting on σ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 (10), u→′ is defined as follows to conserve momentum:(11)u→′=∑σ=1S1/τσ∑α=0Q-1e→αfασ∑σ=1S1/τσ∑α=0Q-1fασ.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:(12)u→=∑σ=1S∑α=0Q-1e→αfασ+1/2∑σ=1SF→σ∑σ=1Sρσ.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 offα(x→,t) without notation σ is(13)fαx→+Δte→α,t+Δt=fαx→,t+fαeqx→,t-fαx→,tτ.In order to recover the Brinkman equation, the equilibrium distribution function fα(eq) of (9) is simplified to (2). ρ=∑α=0Q-1fα but u→(eq) of (10) is modified as follows:(14)u→eq=∑α=0Q-1e→αfα+τΔtρf→mρ,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 (12) is modified to(15)u→=∑α=0Q-1e→αfα+1/2Δtρf→mρ.When solving the Brinkman equation, f→m is a function of u→ defined by (5). A comparison of (3) and (15) shows that the explicit formula of (6) to compute u→ is also valid here. Then, u→(eq) is computed as(16)u→eq=2τu→+1-2τ∑α=0Q-1e→αfαρ,which is obtained by solving (14) and (15).As we can see, the computations off→m by (5) and Fα by (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→≪cs, the computed pressure p=cs2ρ and flow velocity u→ 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κ∗(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 Figure1). This substantially reduces the degrees of freedom in the coarse-grid 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 (18). We seek κ∗ such that the average velocities from local fine-grid simulations with the heterogeneous κ(x→) and homogeneous κ∗, respectively, are equal (see (21)-(22)). 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 (17).Figure 1
Schematic models of the fine and coarse grids.We discuss two-dimensional problems as example. In the local LBM simulations usingκ(x→), we drive flow in x direction by G→const(1)=(Gconst,0) and compute the average velocity u→¯κ(x→)(1), where ·¯ is defined as a volume average over a subdomain. We also compute u→¯κ(x→)(2) by using G→const(2)=(0,Gconst) in another local simulation. Then, κ∗ is computed as follows:(17)κ∗=κxx∗,κxy∗κyx∗,κyy∗=νGconstu→¯κx→1·1,0,u→¯κx→2·1,0u→¯κx→1·0,1,u→¯κx→2·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 (17), (5) is modified to be(18)f→m=-ϵνκ∗-1·u→+ϵG→,where κ∗-1 is the inverse matrix of κ∗. The evolution of fα(x→,t) is described by (2), (13), (14), (15), and (18). As κ∗ and G→ are constant and the periodic boundary conditions are used in local simulations, the relation (19)fαx→+Δte→α,t+Δt=fαx→,tholds at steady state. For arbitrary Δx, Δt, τ, ϵ, ν, κ∗, and G→=G→const, the steady state solution of fα is independent of x→ and equal to(20)fα=ωαρ01+e→αcs2·κ∗·G→constν,which implies that the uniform density is ρ=∑α=0Q-1‍fα≡ρ0. We validate the solution of (20) by the following verification: substituting (20) into (15) and considering (18), we get the uniform velocity u→≡κ∗·G→const/ν. In addition, we get f→m≡0 by (18) using u→. Then, substituting fα and f→m into (14), we get u→(eq)=u→, which implies that (13) is satisfied at steady state since we have fα=fα(eq) according to (2) and fα(x→+Δte→α,t+Δt)=fα(x→,t). According to the uniform solution of u→≡κ∗·G→const/ν and (17), the average velocity u→¯κ∗1 of the local simulation using constant κ∗ and G→const(1)=(Gconst,0) satisfies(21)u→¯κ∗1=κ∗·Gconst,0ν=u→¯κx→1,which implies that the average flux is conserved when using the same external force G→const(1) but different permeability distributions, namely, using the heterogeneous κ(x→) and homogeneous κ∗, respectively. When driving flow by G→const(2)=(0,Gconst), the average flux is also conserved:(22)u→¯κ∗2=κ∗·0,Gconstν=u→¯κx→2.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. κ∗(x→), Δxcoarse, and Δtcoarse in the coarse-grid simulation are different from κ(x→), Δxfine, and Δtfine, 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 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 Δ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 is100×100 and Δx=0.01 m, Δt = 0.0001 s, and τ=0.53 making νeff=0.01 m2 s−1, ν=2×10-6 m2 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(23)κ=10-6,31≤i,j≤70,10-5,otherwise.The average pressure over the whole computational domain is subtracted from the computed pressure p=cs2ρ in all figures of the pressure distributions. The transient results at 5000thΔt and the steady state results at 200000thΔ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.Figure 2
Comparisons ofp, u, and v between two LBM simulations using different force models, transient results at 5000thΔt (a) and the steady state results (b), νeff=0.01 m2 s-1 and G→const=(2,0) m s−2.
(a) (b)


### 4.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.


#### 4.2.1. Layered Distribution ofκ

Here, we uniformly divide the whole domain into 10 layers parallel to they-axis. The odd number layers have κ1=10-12 m2 and κ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 G→const=(2,0) m s−2. τ=0.53 and so νeff=0.01 m2 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.Table 1
Verification of computedκxx∗, κ1=10-12 m2 and νeff=0.01 m2 s−1.
κ 2 / κ 1 1 / 2 1 / κ 1 + 1 / κ 2 - 1 κ x x ∗ by LBM 2 1.33333 × 1 0 - 12 1.33333 × 1 0 - 12 10 1.81818 × 1 0 - 12 1.81818 × 1 0 - 12 50 1.96078 × 1 0 - 12 1.96078 × 1 0 - 12 100 1.98019 × 1 0 - 12 1.98019 × 1 0 - 12 1000 1.99800 × 1 0 - 12 1.99800 × 1 0 - 12 10000 1.99980 × 1 0 - 12 1.99979 × 1 0 - 12 100000 1.99998 × 1 0 - 12 1.99998 × 1 0 - 12When driving flow iny direction by setting G→const=(0,2) m s−2, the velocity distribution along x direction is nonuniform. We set τ=0.5 such that νeff=0 m2 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.Table 2
Verification of computedκyy∗, κ1=10-12 m2 and νeff=0 m2 s−1.
κ 2 / κ 1 1 / 2 ( κ 1 + κ 2 ) κ y y ∗ by LBM 2 1.500000 × 1 0 - 12 1.499999 × 1 0 - 12 10 5.500000 × 1 0 - 12 5.499999 × 1 0 - 12 50 25.50000 × 1 0 - 12 25.49999 × 1 0 - 12 100 50.50000 × 1 0 - 12 50.49999 × 1 0 - 12 1000 500.5000 × 1 0 - 12 500.4999 × 1 0 - 12 10000 5000.500 × 1 0 - 12 5000.499 × 1 0 - 12 100000 50000.50 × 1 0 - 12 50000.49 × 1 0 - 12


#### 4.2.2. Checkerboard Distribution ofκ

As on a checkerboard, we divide the whole computational domain uniformly into10×10 squares with each square containing 10×10 points. The black squares of the checkerboard have κ1=10-12 m2 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 m2 s−1 to get the solution of the Darcy equation. The representative distributions of p, u, and v are given in Figure 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 deviates 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 m2 and νeff=0 m2 s−1.
κ 2 / κ 1 κ 1 κ 2 κ x x ∗ by LBM κ x x ∗ by LBM (1000×1000 points) 2 1.41421 × 1 0 - 12 1.41418 × 1 0 - 12 10 3.16227 × 1 0 - 12 3.14081 × 1 0 - 12 50 7.07106 × 1 0 - 12 6.45938 × 1 0 - 12 7.01357 × 1 0 - 12 100 10.0000 × 1 0 - 12 8.25393 × 1 0 - 12 9.70489 × 1 0 - 12 1000 31.6227 × 1 0 - 12 12.2496 × 1 0 - 12 19.8897 × 1 0 - 12 10000 100.000 × 1 0 - 12 13.0133 × 1 0 - 12 23.2777 × 1 0 - 12Figure 3
Distributions ofp, u, and v, νeff=0 m2 s−1, G→const=(2,0) m s−2, κ1=10-12 m2, and κ2/κ1=2.


### 4.3. Verifications of the Upscaled Simulation Scheme

#### 4.3.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, and ν=2×10-6 m2 s−1. We set νeff=0 m2 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 Figure 4 is set according to (24) such that the distribution of κ∗(x→) is nonuniform.(24)κ=κconst,0.45≤x,y≤0.55,κconst,0.2≤x,y≤0.3,κconst,0.7≤x,y≤0.8,κconst,0.2≤x≤0.3,0.7≤y≤0.8,κconst,0.7≤x≤0.8,0.2≤y≤0.3,101+sin⁡80xπcos⁡80yπκconst,elsewhere,where κconst=10-13 m2. Δxfine=0.0025 m and Δtfine=0.000025 s in the fine-grid simulation. The number of fine points is 400 × 400 inside the whole computational domain that 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 coarse-grid simulation.Figure 4
Distribution of the permeabilityκ(x→), κconst=10-13.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, (18), which is a general formula in the coarse-grid simulations, can be replaced by (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, Δxcoarse=0.025 m and Δtcoarse=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(25)κ∗=κconst,19≤I,J≤22,κconst,9≤I,J≤12,κconst,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 Figure1. Figures 5-6 show that the agreement is very good between the two simulations using the fine and coarse grids, respectively.Figure 5
Comparisons ofp, u, and v between the fine-grid results (a), fine-grid averaged results (b), and coarse-grid results using κ∗ (c), νeff=0 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-13 m2.
(a) (b) (c)Figure 6
Detailed comparisons ofp, u, and v between the fine-grid averaged results and coarse-grid results using κ∗, νeff=0 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-13 m2.


#### 4.3.2. Simulations of Brinkman Flows

The physical problem studied here is similar to that described in Section4.3.1. The differences are that we increase the value of κconst to be κconst=10-7 m2 (cf. κconst=10-13 m2 in Section 4.3.1) and set νeff=10-5 m2 s−1 such that the contribution by the viscosity term νeffΔu→ 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 νeff=10-5 and 0 m2 s−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 νeff is nonzero here, we let νeff=10-5 m2 s−1 when computing κ∗ and get κxx∗(κ,νeff)=κyy∗(κ,νeff)=7.367κconst and κyx∗(κ,νeff)=κxy∗(κ,νeff)=0 for subdomains with κ=10(1+sin⁡80xπcos⁡(80yπ))κconst. For subdomains with κ=κconst, κxx∗(κ,νeff)=κyy∗(κ,νeff)=κconst and κyx∗(κ,νeff)=κxy∗(κ,νeff)=0. Thus, we have κ∗(κ,νeff) equal to 7.367κconstI or κconstI. As discussed in Section 4.3.1, we can use a scalar distribution of κ∗(κ,νeff), which is equal to 7.367κconst or κconst, in the coarse-grid simulation. We use Δxfine=0.0025 m, Δtfine = 0.000025 s, and τfine=0.50012 in the fine-grid simulation and use Δxcoarse=0.025 m, Δtcoarse=0.00025 s, and τcoarse=0.500012 in the coarse-grid simulation. Figures 8-9 show that the agreement of the coarse-grid simulation using κ∗(κ,νeff) with the fine-grid simulation is very good. In addition, we set νeff=0 m2 s−1 when computing κ∗,err and get κxx∗,err(κ)=κyy∗,err(κ)=8.485κconst and κyx∗,err(κ)=κxy∗,err(κ)=0 for subdomains with κ=10(1+sin⁡80xπcos⁡(80yπ))κconst. For subdomains with κ=κconst, κxx∗,err(κ)=κyy∗,err(κ)=κconst and κyx∗,err(κ)=κxy∗,err(κ)=0, where we still use the superscript “err” since the local simulation procedure with νeff=0 m2 s−1 is wrong although the obtained κ∗,err(κ) is the same as the above κ∗(κ,νeff). Now, we have κ∗,err(κ) equal to 8.485κconstI or κconstI. We use a scalar distribution of κ∗,err(κ), which is equal to 8.485κconst or κconst, in another coarse-grid simulation. Note that the difference between κ∗,err(κ) and κ∗(κ,νeff) is distinct in subdomains with nonuniform κ(x→). The results of the coarse-grid simulation using κ∗,err(κ) are also given in Figure 8, which shows that the deviation of the coarse-grid simulation using κ∗,err(κ) from the fine-grid simulation is noticeable. Thus, the previous computation of κ∗(κ,νeff) as the effective permeability using nonzero νeff is accurate for upscaling the Brinkman equation.Figure 7
Comparisons ofp, u, and v between two fine-grid simulations using νeff=10-5 and 0 m2 s−1, respectively, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-7 m2.Figure 8
Comparisons ofp, u, and v between the fine-grid averaged results (a) and coarse-grid results using κ∗(κ,νeff) (b) and κ∗,err(κ) (c), νeff=10-5 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-7 m2.
(a) (b) (c)Figure 9
Detailed comparisons ofp, u, and v between the fine-grid averaged results and coarse-grid results using κ∗(κ,νeff), νeff=10-5 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-7 m2.


## 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 is100×100 and Δx=0.01 m, Δt = 0.0001 s, and τ=0.53 making νeff=0.01 m2 s−1, ν=2×10-6 m2 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(23)κ=10-6,31≤i,j≤70,10-5,otherwise.The average pressure over the whole computational domain is subtracted from the computed pressure p=cs2ρ in all figures of the pressure distributions. The transient results at 5000thΔt and the steady state results at 200000thΔ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.Figure 2
Comparisons ofp, u, and v between two LBM simulations using different force models, transient results at 5000thΔt (a) and the steady state results (b), νeff=0.01 m2 s-1 and G→const=(2,0) m s−2.
(a) (b)


## 4.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.


### 4.2.1. Layered Distribution ofκ

Here, we uniformly divide the whole domain into 10 layers parallel to they-axis. The odd number layers have κ1=10-12 m2 and κ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 G→const=(2,0) m s−2. τ=0.53 and so νeff=0.01 m2 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.Table 1
Verification of computedκxx∗, κ1=10-12 m2 and νeff=0.01 m2 s−1.
κ 2 / κ 1 1 / 2 1 / κ 1 + 1 / κ 2 - 1 κ x x ∗ by LBM 2 1.33333 × 1 0 - 12 1.33333 × 1 0 - 12 10 1.81818 × 1 0 - 12 1.81818 × 1 0 - 12 50 1.96078 × 1 0 - 12 1.96078 × 1 0 - 12 100 1.98019 × 1 0 - 12 1.98019 × 1 0 - 12 1000 1.99800 × 1 0 - 12 1.99800 × 1 0 - 12 10000 1.99980 × 1 0 - 12 1.99979 × 1 0 - 12 100000 1.99998 × 1 0 - 12 1.99998 × 1 0 - 12When driving flow iny direction by setting G→const=(0,2) m s−2, the velocity distribution along x direction is nonuniform. We set τ=0.5 such that νeff=0 m2 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.Table 2
Verification of computedκyy∗, κ1=10-12 m2 and νeff=0 m2 s−1.
κ 2 / κ 1 1 / 2 ( κ 1 + κ 2 ) κ y y ∗ by LBM 2 1.500000 × 1 0 - 12 1.499999 × 1 0 - 12 10 5.500000 × 1 0 - 12 5.499999 × 1 0 - 12 50 25.50000 × 1 0 - 12 25.49999 × 1 0 - 12 100 50.50000 × 1 0 - 12 50.49999 × 1 0 - 12 1000 500.5000 × 1 0 - 12 500.4999 × 1 0 - 12 10000 5000.500 × 1 0 - 12 5000.499 × 1 0 - 12 100000 50000.50 × 1 0 - 12 50000.49 × 1 0 - 12


### 4.2.2. Checkerboard Distribution ofκ

As on a checkerboard, we divide the whole computational domain uniformly into10×10 squares with each square containing 10×10 points. The black squares of the checkerboard have κ1=10-12 m2 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 m2 s−1 to get the solution of the Darcy equation. The representative distributions of p, u, and v are given in Figure 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 deviates 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 m2 and νeff=0 m2 s−1.
κ 2 / κ 1 κ 1 κ 2 κ x x ∗ by LBM κ x x ∗ by LBM (1000×1000 points) 2 1.41421 × 1 0 - 12 1.41418 × 1 0 - 12 10 3.16227 × 1 0 - 12 3.14081 × 1 0 - 12 50 7.07106 × 1 0 - 12 6.45938 × 1 0 - 12 7.01357 × 1 0 - 12 100 10.0000 × 1 0 - 12 8.25393 × 1 0 - 12 9.70489 × 1 0 - 12 1000 31.6227 × 1 0 - 12 12.2496 × 1 0 - 12 19.8897 × 1 0 - 12 10000 100.000 × 1 0 - 12 13.0133 × 1 0 - 12 23.2777 × 1 0 - 12Figure 3
Distributions ofp, u, and v, νeff=0 m2 s−1, G→const=(2,0) m s−2, κ1=10-12 m2, and κ2/κ1=2.


## 4.2.1. Layered Distribution ofκ

Here, we uniformly divide the whole domain into 10 layers parallel to they-axis. The odd number layers have κ1=10-12 m2 and κ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 G→const=(2,0) m s−2. τ=0.53 and so νeff=0.01 m2 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.Table 1
Verification of computedκxx∗, κ1=10-12 m2 and νeff=0.01 m2 s−1.
κ 2 / κ 1 1 / 2 1 / κ 1 + 1 / κ 2 - 1 κ x x ∗ by LBM 2 1.33333 × 1 0 - 12 1.33333 × 1 0 - 12 10 1.81818 × 1 0 - 12 1.81818 × 1 0 - 12 50 1.96078 × 1 0 - 12 1.96078 × 1 0 - 12 100 1.98019 × 1 0 - 12 1.98019 × 1 0 - 12 1000 1.99800 × 1 0 - 12 1.99800 × 1 0 - 12 10000 1.99980 × 1 0 - 12 1.99979 × 1 0 - 12 100000 1.99998 × 1 0 - 12 1.99998 × 1 0 - 12When driving flow iny direction by setting G→const=(0,2) m s−2, the velocity distribution along x direction is nonuniform. We set τ=0.5 such that νeff=0 m2 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.Table 2
Verification of computedκyy∗, κ1=10-12 m2 and νeff=0 m2 s−1.
κ 2 / κ 1 1 / 2 ( κ 1 + κ 2 ) κ y y ∗ by LBM 2 1.500000 × 1 0 - 12 1.499999 × 1 0 - 12 10 5.500000 × 1 0 - 12 5.499999 × 1 0 - 12 50 25.50000 × 1 0 - 12 25.49999 × 1 0 - 12 100 50.50000 × 1 0 - 12 50.49999 × 1 0 - 12 1000 500.5000 × 1 0 - 12 500.4999 × 1 0 - 12 10000 5000.500 × 1 0 - 12 5000.499 × 1 0 - 12 100000 50000.50 × 1 0 - 12 50000.49 × 1 0 - 12


## 4.2.2. Checkerboard Distribution ofκ

As on a checkerboard, we divide the whole computational domain uniformly into10×10 squares with each square containing 10×10 points. The black squares of the checkerboard have κ1=10-12 m2 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 m2 s−1 to get the solution of the Darcy equation. The representative distributions of p, u, and v are given in Figure 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 deviates 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 m2 and νeff=0 m2 s−1.
κ 2 / κ 1 κ 1 κ 2 κ x x ∗ by LBM κ x x ∗ by LBM (1000×1000 points) 2 1.41421 × 1 0 - 12 1.41418 × 1 0 - 12 10 3.16227 × 1 0 - 12 3.14081 × 1 0 - 12 50 7.07106 × 1 0 - 12 6.45938 × 1 0 - 12 7.01357 × 1 0 - 12 100 10.0000 × 1 0 - 12 8.25393 × 1 0 - 12 9.70489 × 1 0 - 12 1000 31.6227 × 1 0 - 12 12.2496 × 1 0 - 12 19.8897 × 1 0 - 12 10000 100.000 × 1 0 - 12 13.0133 × 1 0 - 12 23.2777 × 1 0 - 12Figure 3
Distributions ofp, u, and v, νeff=0 m2 s−1, G→const=(2,0) m s−2, κ1=10-12 m2, and κ2/κ1=2.


## 4.3. Verifications of the Upscaled Simulation Scheme

### 4.3.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, and ν=2×10-6 m2 s−1. We set νeff=0 m2 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 Figure 4 is set according to (24) such that the distribution of κ∗(x→) is nonuniform.(24)κ=κconst,0.45≤x,y≤0.55,κconst,0.2≤x,y≤0.3,κconst,0.7≤x,y≤0.8,κconst,0.2≤x≤0.3,0.7≤y≤0.8,κconst,0.7≤x≤0.8,0.2≤y≤0.3,101+sin⁡80xπcos⁡80yπκconst,elsewhere,where κconst=10-13 m2. Δxfine=0.0025 m and Δtfine=0.000025 s in the fine-grid simulation. The number of fine points is 400 × 400 inside the whole computational domain that 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 coarse-grid simulation.Figure 4
Distribution of the permeabilityκ(x→), κconst=10-13.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, (18), which is a general formula in the coarse-grid simulations, can be replaced by (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, Δxcoarse=0.025 m and Δtcoarse=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(25)κ∗=κconst,19≤I,J≤22,κconst,9≤I,J≤12,κconst,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 Figure1. Figures 5-6 show that the agreement is very good between the two simulations using the fine and coarse grids, respectively.Figure 5
Comparisons ofp, u, and v between the fine-grid results (a), fine-grid averaged results (b), and coarse-grid results using κ∗ (c), νeff=0 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-13 m2.
(a) (b) (c)Figure 6
Detailed comparisons ofp, u, and v between the fine-grid averaged results and coarse-grid results using κ∗, νeff=0 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-13 m2.


### 4.3.2. Simulations of Brinkman Flows

The physical problem studied here is similar to that described in Section4.3.1. The differences are that we increase the value of κconst to be κconst=10-7 m2 (cf. κconst=10-13 m2 in Section 4.3.1) and set νeff=10-5 m2 s−1 such that the contribution by the viscosity term νeffΔu→ 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 νeff=10-5 and 0 m2 s−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 νeff is nonzero here, we let νeff=10-5 m2 s−1 when computing κ∗ and get κxx∗(κ,νeff)=κyy∗(κ,νeff)=7.367κconst and κyx∗(κ,νeff)=κxy∗(κ,νeff)=0 for subdomains with κ=10(1+sin⁡80xπcos⁡(80yπ))κconst. For subdomains with κ=κconst, κxx∗(κ,νeff)=κyy∗(κ,νeff)=κconst and κyx∗(κ,νeff)=κxy∗(κ,νeff)=0. Thus, we have κ∗(κ,νeff) equal to 7.367κconstI or κconstI. As discussed in Section 4.3.1, we can use a scalar distribution of κ∗(κ,νeff), which is equal to 7.367κconst or κconst, in the coarse-grid simulation. We use Δxfine=0.0025 m, Δtfine = 0.000025 s, and τfine=0.50012 in the fine-grid simulation and use Δxcoarse=0.025 m, Δtcoarse=0.00025 s, and τcoarse=0.500012 in the coarse-grid simulation. Figures 8-9 show that the agreement of the coarse-grid simulation using κ∗(κ,νeff) with the fine-grid simulation is very good. In addition, we set νeff=0 m2 s−1 when computing κ∗,err and get κxx∗,err(κ)=κyy∗,err(κ)=8.485κconst and κyx∗,err(κ)=κxy∗,err(κ)=0 for subdomains with κ=10(1+sin⁡80xπcos⁡(80yπ))κconst. For subdomains with κ=κconst, κxx∗,err(κ)=κyy∗,err(κ)=κconst and κyx∗,err(κ)=κxy∗,err(κ)=0, where we still use the superscript “err” since the local simulation procedure with νeff=0 m2 s−1 is wrong although the obtained κ∗,err(κ) is the same as the above κ∗(κ,νeff). Now, we have κ∗,err(κ) equal to 8.485κconstI or κconstI. We use a scalar distribution of κ∗,err(κ), which is equal to 8.485κconst or κconst, in another coarse-grid simulation. Note that the difference between κ∗,err(κ) and κ∗(κ,νeff) is distinct in subdomains with nonuniform κ(x→). The results of the coarse-grid simulation using κ∗,err(κ) are also given in Figure 8, which shows that the deviation of the coarse-grid simulation using κ∗,err(κ) from the fine-grid simulation is noticeable. Thus, the previous computation of κ∗(κ,νeff) as the effective permeability using nonzero νeff is accurate for upscaling the Brinkman equation.Figure 7
Comparisons ofp, u, and v between two fine-grid simulations using νeff=10-5 and 0 m2 s−1, respectively, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-7 m2.Figure 8
Comparisons ofp, u, and v between the fine-grid averaged results (a) and coarse-grid results using κ∗(κ,νeff) (b) and κ∗,err(κ) (c), νeff=10-5 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-7 m2.
(a) (b) (c)Figure 9
Detailed comparisons ofp, u, and v between the fine-grid averaged results and coarse-grid results using κ∗(κ,νeff), νeff=10-5 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-7 m2.


## 4.3.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, and ν=2×10-6 m2 s−1. We set νeff=0 m2 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 Figure 4 is set according to (24) such that the distribution of κ∗(x→) is nonuniform.(24)κ=κconst,0.45≤x,y≤0.55,κconst,0.2≤x,y≤0.3,κconst,0.7≤x,y≤0.8,κconst,0.2≤x≤0.3,0.7≤y≤0.8,κconst,0.7≤x≤0.8,0.2≤y≤0.3,101+sin⁡80xπcos⁡80yπκconst,elsewhere,where κconst=10-13 m2. Δxfine=0.0025 m and Δtfine=0.000025 s in the fine-grid simulation. The number of fine points is 400 × 400 inside the whole computational domain that 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 coarse-grid simulation.Figure 4
Distribution of the permeabilityκ(x→), κconst=10-13.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, (18), which is a general formula in the coarse-grid simulations, can be replaced by (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, Δxcoarse=0.025 m and Δtcoarse=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(25)κ∗=κconst,19≤I,J≤22,κconst,9≤I,J≤12,κconst,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 Figure1. Figures 5-6 show that the agreement is very good between the two simulations using the fine and coarse grids, respectively.Figure 5
Comparisons ofp, u, and v between the fine-grid results (a), fine-grid averaged results (b), and coarse-grid results using κ∗ (c), νeff=0 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-13 m2.
(a) (b) (c)Figure 6
Detailed comparisons ofp, u, and v between the fine-grid averaged results and coarse-grid results using κ∗, νeff=0 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-13 m2.


## 4.3.2. Simulations of Brinkman Flows

The physical problem studied here is similar to that described in Section4.3.1. The differences are that we increase the value of κconst to be κconst=10-7 m2 (cf. κconst=10-13 m2 in Section 4.3.1) and set νeff=10-5 m2 s−1 such that the contribution by the viscosity term νeffΔu→ 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 νeff=10-5 and 0 m2 s−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 νeff is nonzero here, we let νeff=10-5 m2 s−1 when computing κ∗ and get κxx∗(κ,νeff)=κyy∗(κ,νeff)=7.367κconst and κyx∗(κ,νeff)=κxy∗(κ,νeff)=0 for subdomains with κ=10(1+sin⁡80xπcos⁡(80yπ))κconst. For subdomains with κ=κconst, κxx∗(κ,νeff)=κyy∗(κ,νeff)=κconst and κyx∗(κ,νeff)=κxy∗(κ,νeff)=0. Thus, we have κ∗(κ,νeff) equal to 7.367κconstI or κconstI. As discussed in Section 4.3.1, we can use a scalar distribution of κ∗(κ,νeff), which is equal to 7.367κconst or κconst, in the coarse-grid simulation. We use Δxfine=0.0025 m, Δtfine = 0.000025 s, and τfine=0.50012 in the fine-grid simulation and use Δxcoarse=0.025 m, Δtcoarse=0.00025 s, and τcoarse=0.500012 in the coarse-grid simulation. Figures 8-9 show that the agreement of the coarse-grid simulation using κ∗(κ,νeff) with the fine-grid simulation is very good. In addition, we set νeff=0 m2 s−1 when computing κ∗,err and get κxx∗,err(κ)=κyy∗,err(κ)=8.485κconst and κyx∗,err(κ)=κxy∗,err(κ)=0 for subdomains with κ=10(1+sin⁡80xπcos⁡(80yπ))κconst. For subdomains with κ=κconst, κxx∗,err(κ)=κyy∗,err(κ)=κconst and κyx∗,err(κ)=κxy∗,err(κ)=0, where we still use the superscript “err” since the local simulation procedure with νeff=0 m2 s−1 is wrong although the obtained κ∗,err(κ) is the same as the above κ∗(κ,νeff). Now, we have κ∗,err(κ) equal to 8.485κconstI or κconstI. We use a scalar distribution of κ∗,err(κ), which is equal to 8.485κconst or κconst, in another coarse-grid simulation. Note that the difference between κ∗,err(κ) and κ∗(κ,νeff) is distinct in subdomains with nonuniform κ(x→). The results of the coarse-grid simulation using κ∗,err(κ) are also given in Figure 8, which shows that the deviation of the coarse-grid simulation using κ∗,err(κ) from the fine-grid simulation is noticeable. Thus, the previous computation of κ∗(κ,νeff) as the effective permeability using nonzero νeff is accurate for upscaling the Brinkman equation.Figure 7
Comparisons ofp, u, and v between two fine-grid simulations using νeff=10-5 and 0 m2 s−1, respectively, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-7 m2.Figure 8
Comparisons ofp, u, and v between the fine-grid averaged results (a) and coarse-grid results using κ∗(κ,νeff) (b) and κ∗,err(κ) (c), νeff=10-5 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-7 m2.
(a) (b) (c)Figure 9
Detailed comparisons ofp, u, and v between the fine-grid averaged results and coarse-grid results using κ∗(κ,νeff), νeff=10-5 m2 s−1, G→=(sin⁡πx,sin⁡πy) m s−2, and κconst=10-7 m2.


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



---

*Source: 1740693-2017-02-16.xml*