The Influence of Fracture on the Permeability of Carbonate Reservoir Formation Using Lattice Boltzmann Method

Abstract

Due to the complex conditions of carbonate reservoir, in this paper, a unified lattice Boltzmann method was used to study the rule of the flow in carbonate reservoir as the foundation. Two group models were designed to simulate the influences of fractures and vugs. In experiments, at first, the model for carbonate reservoir considering different amounts and lengths of fractures was considered. Then, the model was improved by taking the influences of fractures and vugs into consideration. The result from the first group shows that the whole permeability enhanced a lot when multiple fractures form a big one which connects to two boundaries. At the same time, the main result of the other group shows that the main capability of vugs is flow accumulation. Through a series of experiments, the flow rule in carbonate reservoir with vugs and fractures is proved based on LBM theory, which has a huge impact on the study in LBM and researches on carbonate reservoir.

1. Introduction

Carbonate reservoirs are rich in oil and gas resources and widely distributed in the world, occupying an important position in the distribution of oil and gas. Usually, carbonate reservoirs have the characteristics of large reserves and high production, thus it is easy for them to form large oil and gas fields. However, in carbonate reservoirs, the dissolved vugs, pores, and fractures occupy a great proportion, resulting in serious heterogeneity. This has been a huge obstacle in the research and exploitation of carbonate reservoirs. Various problems not only affect carbonate reservoirs, but also the large-scale exploitation of coalbed methane (CBM) and emerging shale gas. It can be said that the research of fractures is extremely important for energy development.

The exploration of fracture theory was first proposed in 1868 by Russian scientists. It is a very important theory called cubic law, which was popularly used to describe the flow law of most fracture plane flows. In 1969, Snow et al. [1] carried out the research of the fracture system. Their team established a fracture theory using the tensor theory based on statistics which made it very simple and easy to operate, although the calculated results were inaccurate. In later studies, Kamath et al. [2] chose the boundary element method to deal with the situation of one-way flow in fractures. This method was better to match with continuous medium models; therefore, it could be used in a more complex fracture system. Then, Teimoori et al. [3] further improved the model. He made division to fracture porous medium according to the relationship between fractures and matrix, so the calculation process is greatly simplified. The above is about the original contents of the study in fracture theory. Recently, there have been many researches on fractures. Liu et al. [4] numerically studied the geometry and hydraulic characteristics of a three-dimensional intersecting fracture model. In their research, the modified successive random additions algorithm was used for the generation of rough fracture surfaces, and the effect of fracture roughness on the permeability was analyzed. IJeJe et al. [5] used TFReact, a fully implicit coupling simulator, to assess the effect of permeability anisotropy on thermal production and rock property evolution, taking into account each individual factor, including thermal, hydraulic, and mechanical fields. In the development of shale gas reservoirs, Chen et al. [6] utilized multiple interacting continua method to construct fractured reservoir models in order to determine hydraulic fracture geometry of horizontal Wells, which significantly saves computational resources and improves prediction accuracy in their optimization methods. Their works reflect the diversity of research in this field, which has strong guiding significance for our research.

Nowadays, several approaches have been developed for the issues of multiphase flow, including the CT experiment [7,8], the pore network model [9], the lattice Boltzmann method (LBM) [10], the level set method [11], the volume-of-fluid (VOF) method [12], as well as the phase field method [13,14,15]. The research based on these methods for fracture modeling has been greatly developed. Yang et al. [16] used X-ray tomography to analyze the sensitivity of pressure to volume changes in carbonate cavities, which is of great significance for the exploitation in deep carbonate reservoirs. Cai et al. [17] derived a mathematical model based on the classical Lucas–Washburn equation for spontaneous imbibition in a Y-shaped branching network model, which is significant for research on the fluid flow in some pore space of reservoirs. Chen et al. [18] have also used this technique to conduct a series of studies on the pore-scale mechanisms of immiscible and miscible gas injection in fractured carbonate rocks. In their study, fracture–matrix permeability difference determines the ultimate matrix recovery during immiscible gas injection. Wang et al. [19] proposed a fully coupled hydro-thermal model and flow transport in hydraulic fracture is described using discrete fracture model (DFM). They investigated the effect of thermal process on recovery of shale oil considering formation factors and fractal parameters. In the work of Li et al. [20], the LBM method was used for flow simulation and the results showed that the dominant channel pair of natural micro-fractures can significantly increase the permeability of coal seams. LBM is favored by many scholars because it can be programmed on parallel computers and can easily deal with complex boundary conditions. In the 1960s, Broadwell et al. [21] firstly proposed the discrete velocity model, which was the early lattice automata model. In their model, time and space were still continuous physical quantities. In 1970s, Hardy et al. [22] firstly put forward the discrete model called the HPP model, in which the fluid was set as a series of particles, while time and space were scattered to the two dimensional square lattice in the model. With the gradual development of LBM, McNamara and Zanetti [23] introduced a real parameter in 1988 to represent the local particle distribution, and the Boltzmann equation was used to replace the lattice gas automata evolution equation which made it easily apply to numerical calculation. This was the innovation of the lattice Boltzmann method. Compared with the lattice gas automata model, the new model eliminated most system noise, while keeping the advantages of lattice gas automata. The model considered the collision process which could completely ignore the influences between each particle. However, too many particles made the calculation of collision operators very difficult. In 1989, Higuera and Jimenez [24] simplified this model with the equilibrium distribution function collision operator linearization, which overcame the potential problems from the collision operator calculation. Chen et al. [25] put forward the method with single relaxation time in 1991, using the same coefficient of relaxation time to control the speed of different particles closing to grid, which simplified the collision operator. Soon afterwards, Grunau et al. [26] put forward the simplified model, which was based on the single relaxation time model but with simplified collision function, named the lattice BGK (LBGK) model. The model overcame many defects of the lattice gas automata method. It not only greatly simplified the computation of the model, but could also be transferred to the form of Navier–Stokes equations. Therefore, the LBGK model is the most widely used one at present.

Because the lattice Boltzmann method is easy to deal with complex shapes and boundaries in the simulation, the application of LBM is mainly concentrated on pore scale simulation. However, the lattice Boltzmann method can also be extended to the characterization of unit volume (REV) size of numerical simulation.

The basic idea of REV numerical simulation with LBM is to introduce a local resistance to change the local fluid velocity. Spaid and Phelan [27] put forward a single-flow model based on the Brinkman equations. Then, Freed [28] extended the resistance to the tensor theory to improve the model. In the new model, the computed nodes were no longer a pore or a determined point in the skeleton of rock. Instead, it represented the node with a certain volume of porous media. However, due to the linear dependence of velocity and dimensionless resistance, the effects of Mach number and the inertia effect of conventional Mach number are different. Then, Guo and Zhao [29] put forward a model of isothermal incompressible flow in porous media by using the porosity of the introduction of equilibrium distribution function and adding the external force to the evolution equation. The model considered linear and nonlinear partial resistance matrix (matrix drag components) as well as the inertia force and viscous force. In 2002, Kang et al. [30] extended the Freed model to a non-uniform grid, which was a unified model called the unified lattice Boltzmann model (ULBM). The model can simulate different models on diverse scales, and can also be used to simulate the models with the coexistence of multiple scales. The ULBM model can be used on micro scales as small as pores or as macro as the scales in meter or even kilometer.

However, in literatures the grid is rectangular, which makes model application restricted. For example, the fractures in the fracture system flow simulation must have smooth surface, and must be parallel to the axis, which is very difficult.

So as to extend the ULBM model application, we must improve the mesh dividing method. In fact, in order to increase calculation efficiency, many scholars have accomplished a lot of works in the aspect of treating the curve boundary or the local refinement of traditional LBM. One of the popular methods is to combine traditional LBM with quadtree grid. The quadtree grid is a kind of typical unstructured grid, which is composed of many different sizes of square grids. Quadtree grid division offers a certain division standard in a square area, so the classification is a gradual process of mesh recursive subdivision. In recent years, the quadtree grid is gradually popular because this method has the following advantages. Firstly, it is about an automatic dividing grid, which means the grid information storage structure is simple and orderly, and then local high encryption and dynamic adaptability. Crouse et al. first applied a quadtree grid to the lattice Boltzmann method. However, because of the linear interpolation processing mode adopted in the grid with different levels, the improved model cannot keep second-order accuracy of the lattice Boltzmann method. Chen et al. [31] used linear interpolation in the quadtree grid on the basis of a correct method (back-and-forth error compensation and correction, called BEFECC), making the precision corresponding to second-order accuracy.

In this paper, we choose ULBM with quadtree grid to do the research. The model we use is the D2Q9 model for carbonate reservoirs. In experiments, we established two groups of models to simulate the flow in the system with pores and fractures. In the first group, we just consider the effects of fractures. The possible variable parameters include the density of fractures and the permeability of rock matrix. In the second group, we consider the effects of holes on the model. Compared with the first group, the increased possible variable parameters are the size of holes and the quantity of holes, so the research of the second group mainly focuses on the ability of holes.

The quadtree paper only studied the influences of fractures on matrixes with diverse characters as well as the effects of fracture permeability on matrixes when the fracture positions are fixed. However, it neglected the influences of fracture number and density.

2. Numerical Method and Physical Model

2.1. LBM

In this section, we will briefly introduce the unified lattice Boltzmann model. More details can be found in Kang’s articles [30,32]. The basic idea of a unified model is to introduce a resistance to change the local fluid velocity, which not only makes nodes in different areas of the model with diverse scales able to be handled unified, but also makes the inexistence of internal boundary problem in the process of simulation.

In the Bhatnagar–Gross–Krook (BGK) lattice Boltzmann model [33] with a single relaxation time, the evolution equation is:

$$f_i\left(x+e_i{\delta }_t,t+{\delta }_t\right)=f_i\left(x,t\right)-\frac{f_i\left(x,t\right)-f_i^{eq}\left(\rho ,u,T\right)}{\tau }$$

(1)

In the equation, $f_i$ is the particle distribution function in the direction of $i$ (in the D2Q9 model, $i=0,1,2\cdots 8$, as shown in Figure 1), ${\delta }_t$ is the time step, $\tau$ represents the relaxed time, and the relationship between viscosity and relaxed time is: $v=\left(\tau -0.5\right)RT$, $f_i^{eq}$ is the equilibrium distribution function, defined as follows:

$$f_i^{eq}\left(\rho ,u,T\right)={\omega }_i\rho \left[1+\frac{e_i·u}{RT}+\frac{{\left(e_i·u\right)}^2}{2{\left(RT\right)}^2}-\frac{u^2}{2RT}\right],$$

(2)

where $R$ is the gas constant; and $u$, $\rho$, $T$ represent the fluid velocity, density, and temperature, respectively. For the D2Q9 model, the weight coefficient ${\omega }_0=4/9$ when $i=0$; ${\omega }_i=1/9$ when $i=1,2,3,4$; and ${\omega }_i$ = 1/36, when $i=5,6,7,8$. The lattice velocities $e_i$ are $e_0=(0,0)$, $e_{1,3}=(\pm 1,0)$, $e_{2,4}=(0,\pm 1)$, $e_{5,7}=(\pm 1,\pm 1)$, and $e_{6,8}=(\mp 1,\pm 1)$. The fluid density and velocity can be calculated by following the equations: $\rho ={\displaystyle \sum _i{f}_i}$ and $\rho u={\displaystyle \sum _i{e}_i{f}_i}$, respectively.

In the unified lattice Boltzmann model, two velocities, $u^{\prime }$ and $\overline{u}$, are introduced and defined as:

$$u^{\prime }=\left[\underset{\_}{\underset{\_}{I}}+\frac{1}{2}\underset{\_}{\underset{\_}{R}}\right]\left[\underset{\_}{\underset{\_}{I}}-\left(\tau -\frac{1}{2}\right)\underset{\_}{\underset{\_}{R}}\right]•u$$

(3)

$$\overline{u}=\left(1-\frac{1}{2\tau }\right)u+\frac{1}{2\tau }{u}^{\prime }$$

(4)

where $\underset{\_}{\underset{\_}{R}}$ is the resistance tensor, which is related to the permeability tensor $\underset{\_}{\underset{\_}{k}}$ and defined as $\underset{\_}{\underset{\_}{R}}=v{\underset{\_}{\underset{\_}{k}}}^{-1}$. When the resistance tensor is diagonal, Equation can be written as

$${u_{\alpha }}^{\prime }=\frac{1-\left(\tau -\frac{1}{2}\right)R_a}{1+\frac{1}{2}{R}_a}{u}_{\alpha }$$

(5)

At the same time, the equilibrium distribution function is modified as follows:

$$f_i^{eq}\left(\rho ,u,T\right)={\omega }_i\rho \left[1+\frac{e_i·u^{\prime }}{RT}+\frac{{\left(e_i·\overline{u}\right)}^2}{2{\left(RT\right)}^2}-\frac{{\overline{u}}^2}{2RT}\right].$$

(6)

According to the Chapman–Enskog expansion, Equation will recover the following macroscopic transport equations [30]:

$$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho u\right)=0,$$

(7)

$$\frac{\partial \left(\rho u\right)}{\partial t}+\nabla ·\left(\rho uu\right)=-\nabla p-\rho \underset{\_}{\underset{\_}{R}}·u+\nabla ·\left[\rho v\left(\nabla u+u\nabla \right)\right].$$

(8)

Usually, the complex reservoir contains two systems with quite different resistances, e.g., fracture and matrix in fractured reservoir, or hydraulic fracture and matrix in carbonate reservoir. In the fracture system, its permeability can be very large, and the resistance is very small. As a result, the linear velocity term is negligible compared to other terms and Equation (8) will recover the Naiver–Stokes equation. The last term in Equation (8) is called the Brinkman correction, which is usually much smaller than the linear velocity term for flow in porous media. In matrix system, the Brinkman correction term can be neglected, and Equation (8) will be reduced to Darcy’s equation.

2.2. Physical Model

(1)

Physical model I

A series of fracture models, shown in Figure 2, are generated in this work. Four parameters are used to characterize the fracture, including the coordinates of midpoint (x, y), length l, dip angle $\theta$, and fracture aperture b, which are shown in Figure 3. In the fracture model, the fracture density $D_s$ is adopted to describe the number of fractures in two-dimensional plane [34]. The fracture density is defined as:

$$D_s=\frac{{\displaystyle \sum l_i{b}_i}}{A},$$

(9)

where A is the plane area. The following steps are used to generate a model with fracture density $D_0$.

Step 1: define the size of fracture model, including the width X and the length Y for a two-dimensional model. In this work, X and Y are both set as 1000.

Step 2: the coordinates of midpoint (x, y) are generated randomly for each fracture, where $0\le x\le X$, $0\le y\le Y$. In this study, the effect of dip angle is not considered. Therefore, the dip angle $\theta$ is fixed to a certain angle firstly, such as ${30}^{°}$. Then, a small random perturbation $\Delta \theta$ is added to the dip angle, which is subject to $\Delta \theta \le {10}^{°}$ in this case. The length of each fracture is generated by the function expand(x) in Matlab. It needs a mean length $L_m$, which is set as $L_m=64$. The aperture b is set as 1 for all fractures. Then, only the fracture part, which is located inside of the area, should be kept.

Step 3: calculate the fracture density $D_s$. If fracture density is less than the estimated density $D_0$, repeat step 2.

(2)

Physical model II

Two series of fracture models with different sizes and numbers of vugs are generated based on the physical model I. Only the fracture models with a density of 0.02 are chosen from the last section. The vugs are divided into two groups by its sizes, i.e., 100 × 100 for the bigger one and 50 × 50 for the smaller one. For the first group, one, two, three, and four big vugs are added to each fracture model, respectively, while four, eight, twelve, and sixteen small vugs are added in the second group, as shown in Figure 4 and Figure 5. For each fracture density in last section, five models are created. In the figure, only the first model with holes is given while in total there are 40 models.

The locations of the holes are random. For instance, in three-hole models, the third model is produced randomly based on the two-hole models. The volume of a larger hole equals to the summation of four smaller ones, thus the porosities for both models are corresponding.

(3)

Physical model III

Two more series of fracture models vugs are generated, and the fracture density, the size, and number of vugs are different. The size of vugs are kept the same as before: 100 × 100 for the big one and 50 × 50 for the small one. For the fracture models with a density of 0.01, one, two, and three big vugs are added, and three, five, and seven small vugs are added to fracture models with a density of 0.02 randomly, as shown in Figure 6 and Figure 7. The option of this fracture density and vug number is under the consideration of keeping the porosities of fractures and vugs in a similar range. For each fracture density, five models are created. A total of 30 models are built while only the first model with vugs is shown in the figure.

3. Results and Discussion

3.1. The Influence of Fracture Density

To investigate the effect of the fracture density on the flow ability of the formation, six groups of models with different densities were created based on the physical model I. The fracture densities range from 0.01 to 0.06. There are 30 models in total, and only one model for each group is shown in Figure 8. During the simulation, the matrix is set with different permeability and the resistance in the fractures is set to zero. The permeability of matrix $k_m$ is set as 10⁻³, 10⁻², 10⁻¹, 10⁰, and 10¹ (lattice units). The permeability of the whole system $k_p$ will be calculated to study the contribution of fractures to the permeability of matrix model.

The width and length of each model have 1000 pixels. The left side is the entrance and the right side is the exit. In the lattice Boltzmann simulation, all the parameters are dimensionless. The pressure (density) boundary condition [35] is adopted on both sides, while the standard bounceback boundary condition is applied to the other two boundaries. The pressure drop (∆p) between the entrance and the exit is 0.05 (lattice units).

Figure 9 shows the effect of fracture on matrix with different permeability $k_m$, where the y-axis represents the ratio of permeability of fracture model $k_p$ to the permeability of the matrix model without fractures. We calculate the average permeability for each group with the same fracture density and matrix permeability. As shown in Figure 9, the permeability ratio decreases as the permeability of matrix increases. The fracture can change the permeability of the whole system significantly, and the permeability can be 10, even 40 times as the permeability of the matrix system. The fluid flow in the fracture can be dominant in the lower permeable formation. However, with the continuous rise in matrix permeability, the influence of fracture can be gradually ignored. As the matrix permeability increases, the matrix will make more contributions to the fluid flow in the system. At the same time, the matrix occupies more volume in the formation compared to the fractures. The permeability ratio reaches 1 when the matrix permeability is greater than 1, which indicates that the fractures have less contributions to the flowability of matrix due to the higher permeability in matrix system.

On the other side, the permeability ratio increases linearly as the fracture density increases with the same matrix permeability, as shown in Figure 10. This can be attributed to the fact that more fractures will bring more flow channels in the matrix system. It also can be seen that fractures have a more obvious influence on the lower permeability system.

3.2. The Influence of Fractures and Vugs

In some unconventional reservoirs, the properties of various locations are quite different, e.g., the fracture and the matrix in fractured formation, the organic matrix and inorganic matrix in shale, or the fracture/vug system and matrix system in carbonate reservoirs. In this section, a new fracture model with vugs is developed (physical model II). The effect of the size and total volume of vugs will be investigated according to a series of simulations.

The influence of vug size is studied first. The permeability of the matrix system is 0.01, 0.1, and 1, respectively, while other parameters remain the same as the last section. In total, there are 120 examples. By simulation, the permeability of the whole system is calculated with the addition of fractures and vugs. The results of five examples are averaged, which have the same matrix permeability as well as the same vug number and size. The result is shown in Figure 11, in which the x-axis represents the porosity of fractures and vugs while the y-axis is the permeability of the model with fractures and vugs. By comparison, the permeabilities of fracture models with big vugs and small vugs are almost the same, and no obvious differences can be observed. It can be inferred that the size of vugs has no direct relations with the improvement of permeabilities. It thus means the influences of fractures and vugs on fluid flow are different. The interconnected fracture network can improve the flowing ability of the reservoir greatly, while the vugs and disconnected fractures can only increase the fluid reservoir space with no contribution to the flowing. To further verify this conclusion, a set of new models are built for flow simulation.

3.3. The Influence of Fracture Density and Vugs Volume

In this section, the physical model III was used to simulate. The permeability of the matrix system in the simulation is 0.001, while other parameters remain the same as last section. In total, there are 30 examples. The permeability of the whole system with fractures and vugs is calculated and the result is shown in Figure 12. The x-axis represents the porosity of fractures and vugs while the y-axis is the permeability of the model with fractures and vugs. In the same porosity, the permeability of small vugs is generally higher than that of large vugs, because the connectivity between larger vugs and fractures is generally smaller than that between narrow vugs and fractures. The larger vugs corresponds to a stronger heterogeneity, which can be seen as the aggregation of small vugs in one place, while the model with small vugs has a stronger homogeneity, more dispersed vugs, and better connectivity with fractures.

4. Conclusions

In this paper, different fracture systems are simulated by ULBM based on the quadtree grid, and the D2Q9 model is used for carbonate reservoirs. This method has great feasibility in the study of carbonate reservoir and shows a good ability to deal with the complex fracture–vug system. The application of ULBM in carbonate reservoir is also very instructive. By analyzing the results of the simulation, several conclusions can be reached. With the fracture plane density increasing, the fracture percentage rise gradually. Moreover, the whole permeability will enhance a lot when fractures form big fractures connected to two boundaries. However, the number of vugs affects little to the whole permeability. The increase in the number of vugs just reduces the proportion of rock matrix, which indirectly lifts the proportion of flow in fracture, so the whole permeability rises slightly. We also find that the influence law basically remains unchanged in vugs and fractures at different scales. For various sizes of vugs, we can find that, with the increase in vug size, the porosity keeps unchanged and the whole permeability is unchanged as well, except the situation where vugs connect fractures and thus form a connected boundaries fracture. Finally, when the vugs are connected with the fractures, the permeability of the large vugs model is smaller than that of the small vugs model due to its stronger heterogeneity under the condition that the porosity remains unchanged. In future work, we will use ULBM to simulate more complex fracture–vug systems in carbonate reservoirs. For example, corrosion of fractures and vugs, the effect of fluid stress on wall surfaces, and nonisothermal modeling of surfactants on LBM models can be considered.