Abstract
An upscaled Lattice Boltzmann Method (LBM) for flow simulations in heterogeneous porous media at the Darcy scale is proposed in this paper. In the Darcy-scale simulations, the Shan-Chen force model is used to simplify the algorithm. The proposed upscaled LBM uses coarser grids to represent the average effects of the fine-grid simulations. In the upscaled LBM, each coarse grid represents a subdomain of the fine-grid discretization and the effective permeability with the reduced-order models is proposed as we coarsen the grid. The effective permeability is computed using solutions of local problems (e.g., by performing local LBM simulations on the fine grids using the original permeability distribution) and used on the coarse grids in the upscaled simulations. The upscaled LBM that can reduce the computational cost of existing LBM and transfer the information between different scales is implemented. The results of coarse-grid, reduced-order, simulations agree very well with averaged results obtained using a fine grid.
1. Introduction
Detailed flow simulations in porous media are often modeled using the Darcy or Brinkman approximations. In these models, effective parameters, such as absolute and relative permeabilities, depend on the pore-scale geometry. To compute these effective parameters, pore-scale simulations accounting for relevant geometric features in a Representative Elementary Volume (REV) are commonly used as in [1]. The Lattice Boltzmann Method (LBM) [2–4] is well developed for pore-scale flow simulations and extended to model two-phase systems or two immiscible fluids [5–7]. After computing the effective parameters, we are able to perform Darcy-scale simulations using traditional finite volume or element methods used in commercial reservoir simulators. However, these computations are limited to small REVs (compared to the computational grid) and rely on two distinct idealized scale concepts.
Flows at the Darcy scale can also be modeled by LBM with a modified algorithm. The model described in [8] allows particles to partially bounce back at the cells (points) with small permeability. In [9, 10], an external body force, which increases with decreasing permeability, is employed to represent the resistance effect of the porous media to the fluid, where LBM is considered as a unified framework for simulations at all scales. However, these simulations require significant computational resources to converge since the permeability distribution usually has drastic changes in space, which requires a very fine grid for high spatial resolution. To overcome this difficulty, we propose an upscaled LBM scheme.
Following the work in [10] where the generalized Navier-Stokes equation [11] is solved, we simplify the equilibrium distribution function to solve the governing equation at the Darcy scale. In addition, we replace the original Guo et al. force model [12] used in [10] by the simpler Shan-Chen force model [5, 6] to improve the efficiency. Then, an upscaled LBM scheme is proposed to improve computational efficiency by using a coarse grid (each coarse point represents a subdomain) with effective permeability. For each subdomain, the effective permeability is computed by a local scheme, which is based on the conservation principle for the average fluxes (see [13] for general overview of multiscale methods). To avoid the iterative process of finding the unknown effective permeability that satisfies the equation for the average flux, we derive an analytical formula. This analytical formula allows finding the average flux in terms of the effective permeability, which is then inversely determined from the computed average flux using the original permeability distribution in the subdomain concerned.
The computed effective permeabilities are verified in several benchmark problems, where analytical solutions are known. We implement upscaled LBM simulations using the computed effective permeabilities at coarse grids. Agreement between the coarse and fine-grid LBM simulations demonstrates the validity of the upscaled LBM scheme. The average effects of the fine-grid simulations are preserved in the coarse-grid simulations in solving the flow equation at any intermediate scale. Our numerical results show that one can achieve a substantial gain in CPU time by using coarse-grid models. In this paper, the upscaled LBM approach is applied to single-phase flows; however, this approach can be used for modeling multiphase flow phenomena.
2. LBM Algorithms for Simulating Flows in Porous Media
In this section, we discuss LBM algorithms that will be used in our microscale simulations. We first present the LBM algorithm based on the force model proposed by Guo et al. in [12], where an additional term is used in the particle evolution equation to represent the force contribution. Then, we present the general Shan-Chen force model for multicomponent and multiphase systems. In our upscaling algorithm, we focus on the single-phase and single-component model for Brinkman flows. The Shan-Chen force model allows for a more efficient upscaling procedure and overall cleaner presentation. We refer to [14] for more general discussions on LBM algorithms.
First, we will introduce some basic notations associated with LBM. The grid (Lattice) points are uniformly distributed inside the computational domain and , where is the Lattice spacing and is the time step. For two-dimensional problems, we use the D2Q9 Lattice model [4], where is the total number of Lattice velocities, , , , , for 1 to 4, , , and for 5 to 8. For three-dimensional problems, the D3Q19 Lattice model [4] with different and is used, but the algorithms are unchanged.
2.1. LBM Algorithm with the Guo Force Model
Following the algorithm presented in [10], the evolution algorithm of the distribution function iswhere the normalized relaxation time is appropriately selected to match the desired effective kinematic viscosity , where is the sound speed. We use the simplified truncated form of the equilibrium distribution function aswhere the density is given by and the equilibrium velocity, , is defined aswhere is the force per unit mass. Similarly, is simplified to
In the force model proposed by Guo et al. [12], the flow velocity is equal to . If is constant, and are computed by and then and of (1) are determined explicitly. For solving the Darcy-scale equation, we consider the following expressions for as a linear function of (see [10]):where is the porosity, is the physical kinematic viscosity of the fluid, is a scalar for the permeability, and is the external body force per unit mass. The force introduced above incorporates the porous media heterogeneities through the permeability function and depends on the microstructure. If has a high value in the region, then one can assume that this region is highly permeable, while if has a very low value, then this region is almost impermeable. One can also use a forcing that is nonlinear in as an extension to cases with nonlinear Forchheimer effects that are discussed in [10, 15].
In expressions (3) and (5) we have . Using this fact and solving for in (5), we obtain the explicit formula as in [10]: is computed by (5). Then, and of (1) are determined by (2) and (4), respectively.
In the incompressible limit with , the analysis [10] based on the Chapman-Enskog expansion shows that the computed pressure and flow velocity converge to the solutions of the following equation:where is the initial mass density used in LBM simulations. Here, needs not be the real density of the incompressible fluid; then, the computed is used as the pressure of the physical problem. The steady state results of LBM simulations are used as the solutions of the Brinkman equation. The parameters of , , , and can be set independently such that the steady state LBM results converge to the solutions of the continuum Darcy and Stokes equations, respectively.
2.2. Simplified LBM Algorithm with the Shan-Chen Force Model
We now present the general Shan-Chen model and its application to our upscaling scheme. In the original Shan-Chen model [5, 6], which is proposed to simulate multiphase and multicomponent flows, the number of molecules of component having the velocity at and time is denoted by , where and is the total number of components. The general updating algorithm of iswhere and the equilibrium distribution function is defined aswhere ,, and is computed aswhere is related to the total volume force acting on th component. Generally speaking [16], contains three parts: the fluid-fluid interaction , fluid-solid interaction , and external force . For example, for the contribution by the external body force per unit mass. In (10), is defined as follows to conserve momentum:The flow velocity of the whole fluid is equal to the mean velocity before and after implementing the force term and is computed as follows:Recently [15], phase separation process in a fiber geometry and flow of two immiscible fluids in a cross channel are modeled using the Shan-Chen model, which shows the convenience of the LBM in dealing with complex geometries and manipulating the contact angle.
As upscaling in the multiphase problems is a very difficult and often nonlinear procedure, we focus our algorithm first to single-component and single-phase models. For flows of single-component and single-phase, the evolution of without notation isIn order to recover the Brinkman equation, the equilibrium distribution function of (9) is simplified to (2). but of (10) is modified as follows:where we use to replace the original notation, which is equal to the momentum increase per unit volume after due to the force effect through the relaxation process. Correspondingly, the flow velocity of (12) is modified toWhen solving the Brinkman equation, is a function of defined by (5). A comparison of (3) and (15) shows that the explicit formula of (6) to compute is also valid here. Then, is computed aswhich is obtained by solving (14) and (15).
As we can see, the computations of by (5) and by (4) using the computed in the original algorithm [10] are avoided in the current simplified algorithm and, therefore, the efficiency is improved. In the incompressible limit with , the computed pressure and flow velocity also converge to the solutions of the above Brinkman-like equation (7). The same idea can be implemented to (8)–(12) to solve multiphase flows at the Darcy scale. In this way, it may be possible to develop multiphase upscaling techniques based on the upscaling scheme presented below. This is a topic for future work.
3. Upscaling Scheme
For many practical cases, the number of fine discretization points in the whole computational domain due to heterogeneities is very large, making the memory usage and computational time unaffordable. We use an upscaling simulation scheme to reduce the number of points in the fine grid by using a coarse grid with an effective permeability . The upscaled quantities are a tensor quantity even though the input permeability, , is assumed to be a scalar. With this approach we are able to capture fine-grid information on the coarse grid by solving many parallel local problems.
In our proposed algorithm, the computational domain is divided into many subdomains and each subdomain is represented by a coarse point (see Figure 1). This substantially reduces the degrees of freedom in the coarse-grid simulation. To compute the effective for each subdomain, we impose different external forces to drive flows in different directions in the local LBM simulations, which use a fine grid located inside the corresponding subdomain and the distribution of on the fine grid. Then, the similar local LBM simulations usually need to be run with a constant tensor as shown in (18). We seek such that the average velocities from local fine-grid simulations with the heterogeneous and homogeneous , respectively, are equal (see (21)-(22)). The onerous seeking process by adjusting the unknown to match the fluxes computed using is avoided in our simulations since can be computed explicitly by (17).
We discuss two-dimensional problems as example. In the local LBM simulations using , we drive flow in direction by and compute the average velocity , where is defined as a volume average over a subdomain. We also compute by using in another local simulation. Then, is computed as follows:
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 bewhere is the inverse matrix of . The evolution of is described by (2), (13), (14), (15), and (18). As and are constant and the periodic boundary conditions are used in local simulations, the relation holds at steady state. For arbitrary , , , , , , and , the steady state solution of is independent of and equal towhich implies that the uniform density is . We validate the solution of (20) by the following verification: substituting (20) into (15) and considering (18), we get the uniform velocity . In addition, we get by (18) using . Then, substituting and into (14), we get , which implies that (13) is satisfied at steady state since we have according to (2) and . According to the uniform solution of and (17), the average velocity of the local simulation using constant and satisfieswhich implies that the average flux is conserved when using the same external force but different permeability distributions, namely, using the heterogeneous and homogeneous , respectively. When driving flow by , the average flux is also conserved:
After getting the value of 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. , , and in the coarse-grid simulation are different from , , and , respectively. The boundary conditions and the parameters , , , , and 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 direction because fixed quantities are numerically imposed at the initial and last points along direction and their spatial positions are different when using different .
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 and m, = 0.0001 s, and making m2 s−1, m2 s−1, kg m−3, , and m s−2. The periodic boundary conditions are used and the permeability assigned to each point with index isThe average pressure over the whole computational domain is subtracted from the computed pressure in all figures of the pressure distributions. The transient results at and the steady state results at 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.
(a)
(b)
4.2. Verifications of the Computed
We use the above setting of parameters but choose different values for , , and for different problems in Section 4.2. We run simulations in the whole computational domain with prescribed distribution of 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 the -axis. The odd number layers have m2 and in the even number layers is constant with values shown in Table 1 for different cases. Flow is driven in direction by a uniform m s−2. and so m2 s−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 actually yield the solutions of the Darcy equation at steady state.
When driving flow in direction by setting m s−2, the velocity distribution along direction is nonuniform. We set such that m2 s−1 to recover the Darcy equation. As we can see in Table 2, the computed by LBM simulations agrees exactly with the analytical solution.
4.2.2. Checkerboard Distribution of
As on a checkerboard, we divide the whole computational domain uniformly into squares with each square containing points. The black squares of the checkerboard have m2 and in the white squares takes different values for different cases as shown in Table 3. Flow is driven by a uniform m s−2 and we set m2 s−1 to get the solution of the Darcy equation. The representative distributions of , , and are given in Figure 3. The results in Table 3 show that the computed by LBM simulations agrees well with the analytical solution when 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 to to show improving accuracy. and are changed to m and s, respectively. The results given in Table 3 show that the computed 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 × 1 m domain with periodic boundary conditions and , kg m−3, and m2 s−1. We set m2 s−1 by using in the simulations of both fine and coarse grids and also in the calculation of . In order to have obvious variations in the results of the coarse-grid simulation, a nonuniform external force m s−2 is used and the distribution of permeability in Figure 4 is set according to (24) such that the distribution of is nonuniform.where m2. m and 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 fine points located inside the same subdomain are computed and used to verify the results of the coarse-grid simulation.
We have and inside the subdomains, where . For subdomains with , we use m s−2 to drive flow and get and . The symmetric property of inside the subdomain implies that and . We define a scalar as the average value over all diagonal components of and for all subdomains we have , where 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 , 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, m and s. The number of coarse points is and the value of assigned to each coarse point with index is
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.
(a)
(b)
(c)
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 to be m2 (cf. m2 in Section 4.3.1) and set m2 s−1 such that the contribution by the viscosity term 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 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 is nonzero here, we let m2 s−1 when computing and get and for subdomains with . For subdomains with , and . Thus, we have equal to or . As discussed in Section 4.3.1, we can use a scalar distribution of , which is equal to or , in the coarse-grid simulation. We use m, = 0.000025 s, and in the fine-grid simulation and use m, s, and in the coarse-grid simulation. Figures 8-9 show that the agreement of the coarse-grid simulation using with the fine-grid simulation is very good. In addition, we set m2 s−1 when computing and get and for subdomains with . For subdomains with , and , where we still use the superscript “err” since the local simulation procedure with m2 s−1 is wrong although the obtained is the same as the above . Now, we have equal to or . We use a scalar distribution of , which is equal to or , in another coarse-grid simulation. Note that the difference between and is distinct in subdomains with nonuniform . The results of the coarse-grid simulation using are also given in Figure 8, which shows that the deviation of the coarse-grid simulation using from the fine-grid simulation is noticeable. Thus, the previous computation of as the effective permeability using nonzero is accurate for upscaling the Brinkman equation.
(a)
(b)
(c)
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.
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.