Abstract
REV-scale LBM (representative elementary volume scale lattice Boltzmann method) provides a possible way for the unified microscopic simulation of transscale seepage combined with pore-scale LBM. However, shown by numerical results, when simulated fluid flows in tight porous media with discrete microfractures or dispersed dissolution pore, the nonphysical oscillation phenomenon will occur in the region of relatively greater velocity gradient, due to the insufficient computational node caused by computational capability constraints. And for that reason, the increase in simulation scale under certain computational capability will be restricted. In order to solve this problem, this paper extends the application of the original LBM local refinement algorithm into REV-LBM and corrects the original REV-LBM velocity processing method, such that REV-scale LBM has a more reasonable local refinement algorithm. Through conventional computation example, this paper proves that the local refinement algorithm can maintain well computational accuracy while at least double the time efficiency and could save more than 30% of CPU time in the practical application. This algorithm has a great potential of application in simulating flows in larger size porous media with complex micro-nano structures and can provide a new idea to transscale simulation in porous media.
1. Introduction
The computational fluid dynamics (CFD) applied in fluid flow and transportation simulation in porous media includes two types of calculation methods—the conventional N-S equation and the particle solver, including three scale simulation methods: macro continuous media simulation, mesoscopic kinetic simulation, and microscopic molecular dynamics simulation [1]. Among them, the mesoscopic kinetic simulation method based on the mesoscopic kinetic model not only has free constraints in continuity assumptions as the microscopic molecular dynamics simulation does but also maintains the simplicity of the continuous media simulation method which ignores the description of individual molecular motion. Therefore, it is suitable for simulating flow patterns in porous media with micro-nano pore throat structure (e.g., tight sandstone [2, 3], shale [4, 5], fibrous porous media [6], and porous electrode [7]). This simulation method mainly has lattice Boltzmann method (LBM) [8], direct simulation Monte-Carlo (DSMC) [9], dissipative particle dynamics (DPD) [10], etc. Whereas due to the simple and advantageous boundary conditions processing of the LBM method, it is more suitable to simulate flows in complex porous media.
The LBM simulation method for flows in porous media includes pore-scale LBM and REV-scale LBM. The pore-scale combined with the numerical model of porous media [11–13] can be used to study the microscopic flow of complex fluid such as multiphase flow, multicomponent fluid flow [14, 15], and non-Newtonian flow and is the popular trend of micro-nano pore-network seepage simulation. However, when the pore-scale LBM simulates the transscale porous media model, which could have the matrix of large-scale multicontinuum representation and discrete microfractures or dispersed dissolution pore, it will be insufficient due to the excessive number of nodes required for calculation. Unlike pore-scale LBM, REV-scale LBM avoids pore-throat detail characterizations and simulates the flow in the porous media which must have a REV’s size firstly. This method not only includes changing the classical well rebound form to fit the partial rebound model of Darcy equation [16] but also includes the resistance model constructed by adding external force term to transform into Stokes or Brinkman function. Among them, the generalized Lattice Boltzmann equation (GLBE) resistance model proposed by Guo and Zhao [17] is used most extensively. By only adjusting the porosity of the calculated grid, GLBE can accommodate all types of porous medium node, including adjusting the porosity and permeability to nearly achieve a blank or structure node which both referred in pore-scale LBM [18]. Thus, GLBE-based REV-scale LBM has good compatibility with pore-scale LBM. By simulating the fluid transportation in micrometer-scale shale containing microfractures, Zhang et al. [19, 20] confirmed that the combined REV-scale LBM and pore-scale LBM can simulate flows in homogeneous porous media, complex heterogeneous porous media, and single-phase flow in porous media with discrete microfractures. Therefore, GLBE-based REV-scale LBM provides an effective simulation method for the unified simulation of transscale seepage in porous media with micro-nano pore-throat structure. However, even so, shown by numerical simulation results, when simulating the seepage problem related to tight porous media with low permeability, the large velocity gradient generated by the wall boundary (near the fixed structure) or the fracture boundary (junction of the matrix-fracture system) will cause nonphysical oscillation in the boundary region when the calculation nodes are insufficient, thus leading to false perceptions of flow field and restricting the simulation size shown as chapter 4.1 in the reference [21]. Usually, the nonphysical oscillations caused by a lack of nodes can be solved by refining blocks or locally refining nodes to increase the nodes in rapidly varied boundary areas, in which the local refinement approach has more practical benefits. The local refinement algorithm applied to LBM was first proposed by Filippova and Hanel (FH algorithm) [22]. After several improvements [23, 24], it was mainly applied to simulate single-phase and multiphase flow [25, 26] or to explore particle transportation law [27], etc. However, the application to REV-scale LBM in simulating fluid seepage in porous media remains blank, and the REV-scale LBM has been used in the tight porous medium contained microfractures or dispersed dissolution pore with no verification of accuracy at the boundary layer having a large velocity gradient [20].
In order to make up the gap of the existing local refinement algorithm application and solve the nonphysical oscillation problem in the boundary region caused by insufficient nodes, this paper extends the application of the original local refinement algorithm to REV-LBM and adjusts the equilibrium velocity model of Guo’s REV-LBM to make the local refinement algorithm more suitable for REV-scale LBM. In this way, the local refinement algorithm can ensure the calculation accuracy as much as possible under certain calculation capacity constraints, saving calculation nodes, reducing the time consumption of at least one time to reach the simulated equilibrium condition, and enlarging the simulation size as well as the efficiency of the porous media flow simulation with complex microfractures or dispersed dissolution pore, when combined with pore-scale LBM. Therefore, it can provide a direction for future flow simulation of large size tight porous media and develop a new idea for simulations of transscale flow in porous media.
2. Porous Media Flow Simulation by Using REV-LBM with Multigrid Refinement Approach
2.1. REV-Scale Lattice Boltzmann Method
In the LBM method, fluid is assumed to be elements with finite discrete velocities distributed along a regular grid; through the collision and migration evolution, simulation of fluid motion is realized. Similar to the traditional LBM model, the REV-scale LBM model of porous media can establish the seepage resistance model by modifying the traditional rebound format to fit the seepage law or adding the external force term.
Taking the D2Q9 model [28] LBM as an example, its discrete velocity model is shown in Figure 1.
The corresponding parameters are shown in Table 1.
For BGK model or LBM with single relaxation time [29], its specific evolution equation is as follows:
Collision:
Streaming:
Wherein the equilibrium distribution function and the force source term are as follows:
Guo realizes the REV-LBM algorithm by adding an appropriate resistance term of porous media. This algorithm integrates porosity into equilibrium state distribution function and force source term to meet the porous media seepage model, i.e.,
To ensure the computational accuracy, the macroscopic fluid density and macroscopic velocity brought into the equation of equilibrium state are and , respectively. The fluid force in porous media includes the resistance term and the external force term, which is (Where, the first term on the right is the Darcy flow term, the second term is the nonlinear seepage term, and the third term is the external force term. Considering that there are various models of nonlinear seepage terms and external force terms, in this paper, the most universal algorithm is focused on; thus, the most common model is only discussed for Darcy flow to demonstrate the local refinement algorithm of REV-LBM). Since the force model is related to velocity, it is an implicit equation of velocity. Guo’s velocity model treats it as a simple quadratic equation.
When the seepage is classified as Darcy flow without external force term, the equilibrium macroscopic velocity can be simplified as .
To facilitate calculation, all physical quantities are uniformly nondimensionalized, and let , , , and thus, the equation of REV-LBM is transformed into as follows:
Collision:
Streaming:
Equilibrium distribution function:
Force source term:
Force model of porous media:
2.2. Multiblock Local Grid Refinement Approach
2.2.1. Parameter Settings
REV-LBM local refinement mesh model is the same as the conventional LBM local refinement mesh model. For double mesh, the local refinement mesh model has two block systems with different block distances. The schematic model is shown in Figure 2. The lattice spacing of the coarse block is , and the boundary nodes that intersect the fine blocks are represented by blue squares. The lattice spacing of the fine block is , and boundary nodes that intersect the coarse blocks are represented by red circles.
Similar to the conventional LBM local refinement algorithm data processing method, when it denotes the ratio between fine and coarse blocks as to ensure the consistency in the lattice speed of sound between two different block systems, let then:
By using the kinematic viscosity of the macroscopic fluid and the porous media’s local permeability that is independent of block size, it can be obtained that:
The relation between fine block parameters and coarse block parameters can be obtained as follows:
After adjusting the basic physical parameters, it can be found that if the REV-LBM local refinement algorithm still uses the macro velocity processing in the equilibrium equation without external force in Guo’s model, then after a single step duration in the coarse block, due to , the equilibrium velocities of the same block node on coarse and fine block systems are not matched, and the velocity on the coarse block node is slightly higher than that on the fine block node, making the simulation result of local refinement algorithm to produce unreasonable convex phenomenon in coarse block velocity. The convex velocity marked by blue can be observed in the simulation results in Figure 3. To fix this problem, the second term on the right side of the quadratic equation in the original Equation (7) is converted to the summation term, namely:
Similarly, considering that the flow in porous media is Darcy flow without external force, the equilibrium macroscopic velocity can be expressed as follows:
After nondimensionalized, the macroscopic velocity model in the equilibrium equation becomes . After a single step duration in the coarse block, would satisfy the velocity matching of nodes on coarse and fine block systems. Based on the macroscopic equilibrium velocity expression, the coarse block velocity convex problem caused by the local refinement algorithm simulation due to the original velocity processing equation can be improved to a certain degree. Figure 3 illustrates that the macroscopic equilibrium velocity expression is more scientific. Though certain velocity convex in coarse block zone’s nodes can still be observed, the convex is resulted from different cumulative errors on the fine and coarse blocks with different accuracy, and compared with the original equilibrium state velocity processing model, the modified equilibrium velocity processing model in this paper reduces the original velocity calculation error from 1.6% to about 0.8%.
2.2.2. Multigrid Algorithm
To ensure the regular collision, migration, and evolution of the coarse and fine block systems, the nodes on the coarse and fine block need to exchange data. Same as the conventional LBM local refinement algorithm, take distribution functions into consideration , where . Since , . The data exchange method on common points of coarse and fine blocks can thus be obtained.
Coarse to fine:
Fine to coarse:
To avoid the singularity problem occurred in of FH algorithm, the idea of data exchange before the collision is adopted, which based on the DC algorithm [24]. The specific algorithm flow is shown in Figure 4.
Among them, the nodes that need spatial interpolation on the fine block boundary adopt the serendipity element as shown in Figure 5 for interpolation calculation, that is . Since the interpolation element only depends on known nodes, which means different interpolation nodes are independent of each other. Moreover, the interpolation requirements on the four boundaries can be satisfied by rotation, and the calculation formulas before and after rotation are similar, so the parallel ability of the whole algorithm can be improved.
The original REV-LBM algorithm features high parallel computing capability, which mainly focuses on collision and migrations, whereas the local refinement algorithm has multiple parallel parts besides collision and migration parts on different block systems, which also includes the delivery of node values on coarse and fine blocks as well as the calculation of interpolated nodes on fine blocks. Therefore, the local refinement algorithm is capable of turning a single major task assignment into multiple small parallel assignments. Due to the peak performance of a single logical CPU, the local refinement algorithm can take advantage of computational capability for each line by disassembling major assignments under certain computational capability restrictions.
3. Numerical Case
3.1. Poiseuille Flow in Homogeneous Porous Media
The accuracy and efficiency of the local encryption algorithm were verified by simulating the Poiseuille flow in homogeneous porous media. The model size of the porous media is , the local permeability is , the local porosity is 10%, the fluid viscosity is m2/s, and the fluid density is kg/m3. The fluid is driven by fixed pressure difference at both ends, and the centerline flow rate is always maintained at . Nonequilibrium extrapolation scheme is adopted for the inlet and outlet boundary, and the halfway bounce-back scheme is adopted for the upper and lower walls.
The local refinement method of the model block is shown in Figure 6, and the parameters of each simulation method are shown in Table 2.
The simulation results are shown in Figure 7. It can be found that when only coarse blocks are used for simulation, the lack of calculation nodes on the solid wall boundary will cause about 8% abrupt error, while the local refinement algorithm can obviously eliminate the nonphysical oscillation on the boundary caused by the insufficient calculation nodes in a single coarse block system. Besides, although the number of nodes used in the local refinement method of Refinement-1 and Refinement-2 is different, the calculated results of the two methods are similar. Also, either Refinement-1 or Refinement-2 can get results close to fine block accuracy, which shows that as long as the refinement region covers large velocity gradient nodes, the error control of such abrupt velocity nodes can be realized, and the appropriate reduction of the local refinement region will not cause obvious changes in the overall results or cause nonphysical oscillations again.
In order to fully illustrate the efficiency of the local refinement algorithm, this paper doubled the simulation size with the local refinement block model as shown in Figure 6. At this point, the local refinement algorithm will further highlight the advantages of saving nodes and will be more prominent when under constraints of certain computational capability limits; the local refinement algorithm spends a shorter time to obtain a relatively good simulation result. The simulation parameters of each method are as shown in Table 3.
OpenMP was adopted for parallel acceleration in the above three simulation methods, and the maximum number of logical CPU cores for parallel programs was set as 8 cores. After repeated running several times, the running time results are shown in Figure 8. It can be found that when the same number of time steps is simulated, the single duration of each local refinement method is smaller than that of the fine block, and the average duration of the local refinement method is 431.101 s, less than that of the fine block 473.182 s, which exemplifies the time advantage of local refinement method. Moreover, the time required by the local refinement method is about 4.5 times that of a single coarse block. This is because, under a single step duration of the coarse block, the local refinement algorithm needs to carry out two node data exchanges between coarse and fine boundary, two collision migrations in the refinement region, and one collision migration in the coarse block region. Combined with the simulation parameters, the computational task load is found out to be about 4 to 5 times that of the coarse block, meaning that the local refinement algorithm has made full use of the computational capacity of a single logical processor. However, the computational task load of the overall refined block is 4 times that of the coarse block, but the time required is about 5.0 times that of the coarse block only, which shows that the overall refined block approach has exceeded the optimal peak of a single logical processor under the computational force constraints. The local refinement algorithm’s ability to make full use of each line’s computational capability by dissembling a major task under certain computational constraints is illustrated. In addition, the local refinement method has a similar relaxation time with coarse block, thus enabling the local refinement method to reach equilibrium faster and obtain about 2 times the calculation efficiency. At the same time shows the local encryption algorithm compared with the overall elaboration grid in the same stability condition by the simulation speed of at least one times, which also shows that the local refinement algorithm improves the simulation speed at least twice as fast as the overall refined block method under the same stability condition.
(a)
(b)
3.2. Flow in Fractured Porous Media
This method can also eliminate the nonphysical oscillation located at fracture boundary (REV-scale LBM is only used for the matrix, but the fracture flow is calculated by pore-scale LBM. The junction of the matrix-fracture system has a large velocity gradient which will lead to the nonphysical oscillation.). The size of porous media model with real fractures, shown as Figure 9, is . In addition, this model is only contained 28800 nodes which will be cost about twice by means of FEM (finite element method) under the same structure accuracy. And the local permeability of matrix is μm2; the local porosity of matrix is 10%. With the same simulation parameters and settings as the previous case, the simulation result is as follow in Figure 10.
(a)
(b)
From the velocity distribution colormap of the simulation results (Figure 10(a)), it can be clearly observed that many discontinuous color speckles appear at the fracture edge, and discontinuous color bands appear at the solid wall boundary, which is due to the numerical oscillation at the fracture edge and solid wall boundary caused by the original single coarse grid algorithm. Compared with Figure 10(a), it can be found that the discontinuous color speckles and color strip disappear in Figure 10(b).
As the partial enlarged view of Figure 10(a) is shown in Figure 11, it can be discovered that the bright speckles (green circles) appear at a certain distance from the fracture edge, while the dark speckles (red dash circles) are almost close to the fractures edge. This phenomenon is consistent with the previous case (Section 3.1), which shows that the single coarse grid algorithm will cause extremely serious nonphysical oscillation when the velocity gradient is large and will cause great interference to the simulation results of the porous media model. The smooth transition of the velocity distribution near the fractures edge in Figure 10(b) also proves the superiority of the numerical accuracy of the local refinement algorithm. What is more, it could save more than 30% of CPU time in this practical application, which is similar as the previous conventional example mentioned in Section 3.1.
4. Conclusion
In this paper, the existing LBM local refinement algorithm is extended to combine with the REV-scale LBM algorithm, and because of that, REV-scale LBM is more applicable to solve the complex porous media seepage problem independently or with the cooperation of pore-scale LBM. By adjusting the processing of macroscopic velocity in the equilibrium state equation in Guo’s original model, the problem of velocity mismatch on coarse and fine block caused by the application of the original model to the local refinement algorithm is eliminated. The simulation results show that the local refinement algorithm based on REV-scale LBM can achieve high computational accuracy with limited computational capability and less computational time (nearly 9 times error reduction at the boundary layer with only 70% CPU time), which can provide an idea for the flow simulation of transscale porous media with matrix meeting the REV’s size and complex microfractures or dissolution pore structures.
Symbols
| : | Discrete velocity set, m/s |
| : | The lattice speed of sound, m/s |
| : | Conversion factor for length, m |
| : | Conversion factor for time, s |
| : | Conversion factor for density, kg/m3 |
| : | Total derivative, kg/(s·m3) |
| : | Particle distribution function with velocity at position and time , kg/m3 |
| : | after collision, kg/m3 |
| : | Body force, N/m3 |
| : | External force, N/m3 |
| : | Porous channel height, m |
| : | Local permeability, m2 |
| : | The ratio of lattice spacing between two blocks, dimensionless |
| : | Velocity, m/s |
| : | Initial velocity, m/s |
| : | Centerline velocity, m/s |
| : | Lattice spacing, m |
| : | Time step, s |
| : | Local porosity, dimensionless |
| : | Turbulence coefficient of high-speed non-Darcy flow, m−1 |
| : | Dynamic viscosity, Pa·s |
| : | Kinematic viscosity, m2/s |
| : | Relaxation time, s |
| : | Density, kg/m3 |
| : | Weight with velocity , dimensionless. |
| : | Equilibrium state |
| : | Nonequilibrium state |
| : | Dimensionless parameter. |
| : | Coarse grid |
| : | Fine grid |
| : | Parameter with velocity |
| : | Tensor index. |
Data Availability
The numerical data and algorithm, except the information presented in the manuscript, used to support the funding of this study is restricted by the safety law of Petrol China, while the source code of the local refinement algorithm written in the manuscript is not open source before commercialization.
Conflicts of Interest
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Acknowledgments
This study is financially supported by the National Science and Technology Major Project (2017ZX05013-001), the Research on Comprehensive Evaluation Method of Unconventional Reservoir Based on Artificial Intelligence (2019D-500809), the Study on the Seepage Law of Typical Low-Grade Oil Reservoirs and New Methods for Enhancing Oil Recovery (2021DJ1102), and the Research on Tight Oil Physical Simulation and Production Mechanism (2021DJ2204).