Abstract
A fractured-vuggy carbonate reservoir is a special reservoir formed by long-term physical, chemical, and geological processes. Its reserves are large in scale and widely distributed, showing the characteristics of free flow-seepage coupling. Conventional simulation is usually simplified by equivalent permeability, which cannot reflect the actual development characteristics. Given this, the flow in caves and fractures is treated with free flow, using the Navier–Stokes equation. The seepage simulation is used for other areas, and the Darcy formula is used. Finally, the simulation results are obtained by coupling, and the influence of oil production speed, fracture-cavity size, fracture-cavity location, dynamic viscosity, permeability, and other factors on bottom pressure is analyzed to effectively guide the field development. The results show that the production pressure of fractured-vuggy reservoirs diffuses from the central fractured-vuggy area to the surrounding matrix, and the pressure increases from the fractured-vuggy area to the surrounding matrix. The flow velocity in the seepage area is relatively stable and flows gently into the middle fracture cavity from all directions. There will be eddy current in the free-flow area. Different factors have different effects on the development. The oil production speed and oil dynamic viscosity are positively correlated with it, while the formation permeability is negatively correlated with it. The size and location distribution of fracture cavity will also have a certain impact. Simulation in advance can effectively avoid some reservoir development problems.
1. Introduction
Tahe oilfield reservoirs experienced multiphase tectonic movement and karstification, covered in Ordovician, under the influence of soluble rock seam between the growth of body piercing and fault belt; the tectonic deformation belt formed the good matching relationship, mostly occurring in deep fault belt as its core dissolution expansion of favorable reservoir trap formation, which resulted in the formation of fractured-vuggy reservoirs with fault-controlled karst characteristics [1–5]. As a new type and target of deep marine carbonate oil and gas exploration and development, it is urgent to clarify the production characteristics of this type of reservoir to provide a theoretical basis for the efficient development of crude oil in this type of reservoir [6, 7]. Fractured-vuggy reservoirs in fault-solution reservoirs are developed along the main fault zone, which is banded and diverged in some parts. The main reservoir space is a large cave, the fracture zone is wide, and the vertical faults with high angle are developed. In the early stage of mining, the main mining is exhaustion, and the phenomenon of water flooding occurs due to the channeling flow of bottom water along the fault zone [8–10]. Therefore, the understanding of the production characteristics (including dynamic variation characteristics of pressure field and dynamic variation characteristics of velocity field) in fractured-vuggy reservoirs in the early stage of depletion is helpful to provide flow field analysis for controlling water flooding caused by bottom water coning and then to take water control measures. However, different from conventional carbonate reservoirs, the main reservoir space and flow space of faulted fluid reservoirs are large caves, surrounded by fracture zones and fractures. The flow law of crude oil in faulted fluid reservoirs cannot be uniformly described by the traditional Darcy seepage theory [11–13]. A typical single cavity melt broken seam hole reservoir geological model of karst cave inside the free flow of oil displacement by N-S equation was built using COMSOL physical field coupling numerical simulation software. The fracture zone and fracture around the cave unified seepage area were described by Darcy formula, which were equivalent to two junction coupled area via BJS conditions. The evolution law of pressure field and velocity difference in the exploitation of solvent reservoir is analyzed, and the influence of exploitation speed, location, size, and viscosity of crude oil and permeability on solvent reservoir is studied.
2. Establishment of a Mathematical Model of Free Flow-Seepage Coupling in Fractures and Caves
A faulted solution reservoir is a special fractured-vuggy reservoir. Influenced by multistage geological tectonic movement and karstic action, its reservoir space is mostly large karst caves and dissolution pores located in the fracture zone and the fault is the main flow channel, as shown in Figure 1. When fluid flows in fault-solution reservoirs, its flow characteristics are a complex free flow-seepage-combined flow system, including the free-flow region of the cave system, the seepage region composed of matrix rock blocks and fractures, and the coupling region between the two regions using BJS conditions [14, 15].
(a)
(b)
(c)
2.1. Flow Equation in Free-Flow Region
When the fluid flows in a large space, the fluid has little influence on the matrix, and the flow at this time is called free flow. When the fluid flows in the fracture reservoir cave system, it is free to flow [16, 17]. The law describing the relationship between stress and strain is the conservation of dynamic quantity. The fluid motion is controlled by the continuity equation and the momentum equation, which are expressed in vector form as follows:where is the density (kg/m3), u is the fluid velocity (m/s), is the gradient operator, f is the physical force per unit mass (m/s2), and is the fluid stress tensor (Pa). For the sake of simplicity, we will consider only Newtonian fluids. The density of a slightly compressible fluid can be assumed to be almost constant, i.e., incompressible. Therefore, the controlled fluid motion equation can be summed up as the N-S equation, which is described as follows:where the pf is the pressure (Pa) and u is the fluid viscosity (Pa·s). In addition, the stress tensor σ of incompressible Newtonian fluid satisfies the following equation:where is the unit tensor and e is the strain rate, which satisfies
2.2. Flow Equation in the Seepage Zone
The flow of fluid through porous media is called percolation. A porous medium is a material consisting of a solid skeleton and interconnected pores, fractures, or capillary tubes of various types [18]. In this paper, the seepage zone includes a macroscopic fracture system and porous rock matrix. The flow of fluid through large fractures should be described as free flow because there is more space in large fractures. However, the flow model in the crack can be simplified to a parallel-plate laminar flow model and described as a Darcy formula form with equivalent flows. Therefore, the flow of the fluid in the matrix rock block and fracture system is treated as seepage flow in this paper, which follows the classical Darcy’s law and develops along the direction of mass conservation. The viscosity and compressibility of the fluid are considered, and the molecular force and physical-chemical action are not considered. The formula is shown as follows:where V is the Darcy velocity, the average volume velocity (m/s), and K is the permeability tensor (m2), pd is the average pore pressure (Pa), and Φ is the porosity of porous media. Boundary conditions of the equation include two parts: pressure condition (equation (9)) and fluid flux conditions (equation (10)):
is the specific pressure on the pressure boundary and is a specific value or expression of the inward Darcy flux (m/s) on the fluid flux boundary .
2.3. Coupled Boundary Conditions of Percolation and Free Flow
Suitable interface conditions should be introduced at the interface between the seepage zone and the free-flow zone. In this paper, the depletion and production process of a single-vuggy reservoir is simulated. It is a single first-order field problem using BJS condition coupled with Darcy equation and N-S equation [19–21]. At the junction of the seepage zone and free flow, the pressure field and velocity field satisfy where is the boundary between the seepage zone and free flow, is the pressure in the free-flow region, is the pressure in the seepage zone, is the velocity in the free-flow region, and is the velocity of the seepage zone.
Initial conditions include initial formation pressure and initial formation temperature . Boundary conditions include outlet boundary conditions satisfying bottom hole constant flow production , elastic reservoir conditions meet the compressibility of fluid and rock, and energy depletion is provided by the elastic reservoir.
3. Establishment of Geological Model of Fractured-Vuggy Reservoir
Because the reservoir is modified by multiple karst, the storage space is mostly dissolution holes and large caves formed along the fracture zone, so COMSOL software is used to establish the typical karst caves, with a geological model of single karst fracture-cave reservoir controlled by fault and a single-phase flow process of crude oil during reservoir depletion production process.
According to the formation fluid parameters, using COMSOL software, firstly, a cuboid with the length of 100 m, the width of 100 m, and height of 50 m is built to simulate the matrix of fault-controlled karst fractured-vuggy reservoir enveloping a single cave body. Then, a cube cave with a side length of 10 m was constructed in the center. Finally, a cylindrical production pipe through the internal karst cave and the external space is constructed, with an inner diameter of 0.08 m and a height of 22.5 m. The fracture-cavity area and the production pipe are selected to construct a union and form a combination. The edge of the whole mold is selected, and the boundary condition is set as no-slip geological models.
Based on the 3D geological model in Figure 2(a), fluid properties are preliminarily used at a density of 900 kg/m3, the dynamic viscosity is 0.1 Pa·s, the initial mining pressure is 40 MPa, and the mining outlet speed is controlled to 0.001 m/s. For the porous media matrix, the porosity rate is set to 0.005, and the permeability is set to 5 × 10−15 m2. Secondly, in the water storage model, set the compressibility of the fluid as 4 × 10−9[1/Pa], the effective compressibility of the matrix is 4 × 10−10(1/Pa). The tetrahedral mesh was used for 3D mesh division. The whole geometry was selected from the list of several solid layers and calibrated to fluid dynamics. The maximum element size was 9.19 m, the minimum element size was 0.2 m, the maximum element growth rate was 1.25, the curvature factor was 0.8, the narrow area resolution was 0.5, and the number of iterations was 4. After the parameters are set, grid construction is carried out to form the geological model grid of single-vuggy fault-controlled karst fractured-vuggy reservoir as shown in Figure 2(b). The complete grid contains 114,861 domain units, 4084 boundary units, and 660 edge units. In the calculation process, the output time of the 3D model is set at 300 days.
(a)
(b)
(c)
To facilitate the comparative analysis of the flow pattern and pressure variation characteristics of the porous media matrix zone in the free-flow area around the bottom of the well and in the seepage area, a longitudinal section was cut along the center of the well, as shown in Figure 2(c). The two-dimensional profile mesh was divided by a triangular mesh. The entire geometry was selected from the list of geometric solid layers and calibrated to fluid dynamics. The maximum element size was 6.7 m, the minimum element size was 6 m, the maximum element growth rate was 1.2, the curvature factor was 0.4, the narrow area resolution was 1, and the number of iterations was 8. After the parameters are set, grid construction is carried out. The complete grid contains 7684 domain units and 504 boundary units. The calculation process is analyzed according to the output time of the two-dimensional model is 6 days per unit.
4. Results and Discussion
4.1. Evolution Characteristics of Pressure Field and Velocity Field
Based on the free flow-seepage-coupled fractured-vuggy reservoir, this paper analyzes the characteristics of the development process of the fractured-vuggy reservoir from two aspects of fluid pressure and velocity:(1)For pressure analysis, select a point at the bottom of the production pipe, which is located at the upper part of the fracture hole and can typically represent the pressure change of the whole fracture hole. In addition, for the whole reservoir matrix area, the change of the overall pressure field from beginning to end can be seen through the change of tone, as well as the pressure comparison between the free-flow area and the seepage area at the same time.(2)For the analysis of velocity, streamline with arrows is used to clearly illustrate the flow direction of the fluid. Meanwhile, the density of the streamline in the velocity field can also show the velocities in different regions, and the shape of the streamline can intuitively reflect the flow state of the fluid.
To make the simulation results more convincing, the flow of the fluid in the free-flow region and the seepage region is mainly affected by the pressure brought by gravity and matrix, and the viscous force has little influence on the fluid. Only the dynamic viscosity of the extracted oil is considered, and the influence of the viscous force is not studied too much.
According to experience, in the process of reservoir development, the pressure field in the pure seepage zone will be wide at the top and narrow at the bottom, but from the actual simulation, it can be seen that the upper and lower widths of the pressure field in the fractured-vuggy reservoir are not significantly different but appear column. The lateral width of the lower pressure field is likely widened due to the relationship between the fracture-cavity region. The pressure change spreads from the bottom center of the production line to the periphery. Under the assumed physical numerical conditions, the pressure of fractured vuggy decreases with the development of reservoir production. By comparing the vuggy and porous media areas, it can be seen that the pressure in the free-flow area is always lower than that in the seepage zone at any time of reservoir production (Figure 3).
Further comparative analysis of the two-dimensional velocity fields in different periods (Figure 4) shows that the color distribution of the velocity fields is relatively uniform, indicating that the overall velocity distribution is stable. According to the second figure in Figure 3, an obvious velocity peak appears at the bottom of the oil production pipe between 0 and 1 d and then decreases and becomes stable. The flow line of the porous media area is relatively straight, the upper part is straight, the lower part is curved, and it gently merges into the middle fracture-cavity area from all directions. The flow line of the free-flow fracture-hole area is spiral, and the eddy current phenomenon is present.
From the 3D pressure diagram (Figure 5), it can be seen that the pressure distribution also presents a cup-shaped or column-shaped structure, which diffused from the center of the slit to the surrounding, and presents a petal-like at a certain moment. The pressure decreases rapidly at the center of the upper surface of the matrix, while the pressure decreases slowly at the center of the lower surface. The upper part of the side of the square matrix is not a single trend change, but a wave shape. Considering the square matrix, oil products will affect each other at the vertex of the cube during the mining process, and the closer it is to the vertex, the greater the degree of influence.
From the 3D view, a 3D drawing group was built to describe the seam hole area. The uniform density streamline was added with a 0.05 interval distance. The 3D velocity field streamline data of the free flow area at day 60 were selected for simulation analysis, as shown in Figure 6. It can be seen from the figure that the central flow line of the free flow area is more concentrated, and the edge flow line is more scattered, intuitively showing the shape of a funnel. At the same time, there is also an obvious vortex flowing into the center in the middle and upper area, which is consistent with the spiral flow line in the previous two-dimensional section.
(a)
(b)
In order to analyze the coupling situation of free flow-percolation, based on the model constructed in this paper, the central cube is selected in the three-dimensional model for analysis. With low coupling complex, the pressure decreases fast, so the coupling situation of free flow-seepage can be reflected by the change of surface pressure to some extent. Datas on day 60 were selected as the research object, and the results are shown in Figure 7. The cube corner position, especially the lower pressure, drops slowly, which indicates that the coupling complexity is relatively high. The central position pressure drops rapidly, indicating a low coupling complexity. The two views below are the upper surface and the bottom surface of the square, respectively. The area of the blue low-pressure area on the upper surface is significantly larger than that on the bottom surface, and the same result can be obtained.
(a)
(b)
4.2. Analysis of Influencing Factors of Mining Characteristics
The formation of carbonate reservoir is more complex than other reservoir types because the geological hole is the main space of the oil storage and oil speed. The subjective factors such as oil production speed and wellhead pressure, and the objective factors such as oil attribute, seam cave type, and matrix permeability all have an important influence on the development of seam cave carbonate reservoir. Song et al. based on this model studied the influence of gravity, fracture opening, water injection rate, and other factors on fluid flow characteristics and oil displacement effect [19]. This paper selected several typical factors for reference to analyze the influence law of different factors on mining characteristics [20–26].
4.2.1. Mining Rate
In the development and production process of oil and gas fields, while maintaining stable harvesting, a high mining speed should be selected as far as possible to ensure the actual production efficiency, so as to meet the national and social demand for oil. Therefore, selecting an appropriate and reasonable oil recovery speed will not only improve the reservoir development efficiency but also achieve the best overall economic effect [27–29]. To analyze the influence of oil recovery rate on the development of fractured-vuggy reservoirs, this paper sets the oil recovery rate as 172.8 m/d, 86.4 m/d, and 43.2 m/d, and the pressure field and velocity field distribution of different production rates are obtained, as shown in Figure 8(a). It can be seen from the figure that, in the production process of fractured-vuggy reservoirs, the production rate will affect the change of bottom hole pressure. At the same time, the blue low-pressure area will appear earlier when the speed is higher. As shown in Figure 8(b), the changing trend of bottom hole pressure decreases with time, and the faster the oil recovery rate, the faster the change of bottom hole pressure. The slower the recovery rate, the slower the change in bottom hole pressure. In the seepage zone, the flow rate is stable. In the free-flow region, the maximum value appears at 0.2 d and 3 d. When the oil recovery rate reaches 172.8 m/d, the bottom hole pressure will be pumped to 0 MPa in about 2.5 d. Continuing the production will cause negative pressure, which will affect the production. Therefore, the proper oil recovery rate should be selected to avoid negative pressure in the actual development process.
(a)
(b)
4.2.2. Location of Karst Cave
Both caves and fractures are formed in the stratum after long-term physical and chemical action, and their positions are random. The depth of the fracture-cave area also affects the burial depth of the reservoir [30, 31]. To analyze the influence of fracture-cavity positions on the development of fracture-cavity reservoirs, this paper selects three fracture-cavity positions in the upper part, the middle part, and the lower part for simple comparative analysis (Figure 9(a)).
(a)
(b)
The pressure distribution map on the left can be used to determine that when the oil reservoir is shallow, that is, the fracture-cavity area is close to the upper part, the streamline at the lower part of the matrix is bent more, and the pressure spreading range is smaller. When the oil reservoir is deep, that is, the fracture-cavity area is close to the lower part of the matrix, the streamline bending degree of the lower part of the matrix is smaller, and the range of pressure diffusion is larger. According to the pressure variation diagram of different fracture positions (Figure 9(b)), it can be seen that the position of the fracture has a certain influence on the reservoir development process to some extent. When the fracture is located at a higher position, the pressure drops quickly and the flow line in the fracture area is stable. When the cavity is located at a lower position, the pressure drops more slowly. As for the change of flow velocity, there is no great difference in the flow velocity of the three positions, and the flow line in the fracture-hole area is messy. Compared with the middle position, there is no obvious spiral flow line in the free-flow area of the upper and lower fractures.
4.2.3. Cave Size
To analyze the influence of the size of the slit on the development, the simulation completed the analysis by fixing the central position of the slit and changing the side length of the square slit region. The side length of the slit region was assumed to be 10 m, 20 m, and 30 m, respectively (Figure 10(a)). The larger the size of the fracture means that the contact area between the free-flow zone and the seepage zone increases, that is, the area of the coupling zone between the two increases. Therefore, the larger the size of the fracture, the faster the reservoir flow and the faster the pressure drop; the smaller the size of the fracture, the slower the reservoir flow and the slower the pressure drop. In the seepage zone, the flow rate is stable. In the free flow region, the maximum value appears at 0.2d and 0.6d, respectively. Meanwhile the larger the edge length, the more obvious the vortex flow and the larger the vortex range. According to the one-dimensional map of pressure change, in the early stage of mining, the pressure drops slower when the size of the fracture cavity is larger, and faster when the size of the fracture cavity is smaller. Late in mining, from the pressure variation (Figure 10(b)), after about 2.5 d, this trend will have a change, if seam hole size is larger, the pressure drop, and if seam hole size is small, the pressure drop, but overall, the size of the slot hole size will have little impact on the bottom hole pressure, and bottom hole pressure continues to fall over time.
(a)
(b)
4.2.4. Crude Oil Viscosity
Crude oil viscosity refers to the internal friction resistance caused by oil in the process of flow. It is a very important parameter for oil and gas transportation, accumulation, and development of oil and gas fields. The crude oil viscosity of 7 mPa·s, 5 mPa·s, and 1.28 Pa·s was selected for simulation analysis, and the results are shown in Figure 11.
(a)
(b)
Figure 11(a) shows that the bottom hole pressure drops faster at the viscosity of 5 mPa·s and 7 mPa·s, resulting in an earlier blue low-pressure zone compared with the initial standard of 1.28 Pa·s. According to the velocity field cloud diagram on the right side of Figure 11(a), the flow velocity is relatively stable in the seepage zone. In the free-flow region, the maximum value occurs at 0.2 d. Figure 11(b) shows that the viscosity of crude oil affects the bottom hole pressure. The greater the viscosity, the faster the pressure decreases. The lower the viscosity, the slower the pressure reduction, which is consistent with previous conclusions. When the viscosity of crude oil is 7 mPa·s, the pressure will reach 0 at 2.5 d. When the viscosity gradually increases, the folding point due to the elastic pressure of the matrix gradually becomes unclear. In the figure, when the dynamic viscosity is 7 mPa·s, the folding line is approximately a parabola with an upward opening.
4.2.5. Permeability of Matrix Fracture Zone
Fractured-vuggy carbonate reservoirs are widely distributed throughout the world. However, the geological strip is different in different parts of the world and therefore needs to be considered in the development process. In this paper, the influence of matrix permeability on development was considered separately. Permeability is a parameter representing the ability of rock itself to conduct liquid. After several parameter changes, 10 mD, 1 mD, and 0.5 mD permeability were finally selected for analysis. It can be seen from Figure 12(a) that the size of matrix permeability will affect the change of bottom hole pressure. At the same time, when the permeability is low, the cyan part will appear earlier in the middle fracture-hole area, indicating that the pressure drops quickly; on the contrary, when the permeability is high, the pressure drops slowly. In addition, it can be seen from the one-dimensional pressure diagram (Figure 12(b)) that the higher the permeability is, the slower the pressure drops. The lower the permeability is, the faster the pressure decreases. After about 0.4 d, the slope gradually slows down and tends to be stable. Considering the elastic pressure of the matrix mentioned above, it can be further concluded that the higher the permeability is, the smaller the elastic pressure is. The lower the permeability, the greater the elastic pressure. When the permeability is higher, the streamline distribution is more uniform. When the permeability is low, in the free-flow area, the maximum value appears at 0.2 d.
(a)
(b)
5. Conclusion
Large karst cave-fracture zone is the main oil storage space of seam cave-type carbonate reservoir. In this paper, the three-dimensional cavity model is established for the cavity reservoir simulation of free flow-seepage coupling. And reasonable fluid data and conditions are selected as simulation examples for the failure development of fracture-vug type reservoir. The results show that the simulation method can directly reflect the development characteristics of field operation under real conditions, so as to effectively guide the field development, which has certain significance for the continuous and efficient development of seam reservoir and relieve the energy tension.(1)In the development process of fractured-vuggy reservoirs, the production pressure of the reservoir spreads from the central fractured-vuggy area to the surrounding matrix. The two-dimensional section presents the cup shape and the three-dimensional petal shape. The pressure increases from the fractured-vuggy area to the surrounding area. For different oil properties and geological conditions, the declining trend of bottom hole pressure is different, so it is necessary to detect well the variation of strata pressure in actual depletion exploitation. Through the pressure curve, it can be analyzed that due to the elastic force of the matrix, the bottom hole pressure will appear at the initial bending point, and the geological conditions will have an impact on it.(2)For the change of reservoir velocity in the development process, the seepage zone is relatively stable and gently flows into the central fracture-cavity area from all directions. Vortices will appear in the free-flow area, and at the early stage of the development process, areas with increased velocity will appear at random locations. On the whole, taking the horizontal plane where the fracture-hole area is as the standard, the upper streamline is dense and the flow velocity is fast, while the lower streamline is sparse and the flow velocity is slow.(3)Many factors can affect the development of seam-type reservoir, such as oil recovery speed and oil dynamic viscosity are positively correlated, the permeability is negative, and the size and position distribution of the hole will also have a certain impact. Therefore, according to different geological conditions and reservoir attributes, different mining methods should be adopted in the actual operation process to achieve the best effect.
Data Availability
The data used to support the findings of this study are included within the supplementary information files.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this manuscript.
Acknowledgments
This research was funded by the National Science and Technology Major Project of China (2016ZX05053) and the Science and Technology Department Project of Sinopec-China Petroleum (P11089). The authors acknowledge the technical support provided by the China University of Petroleum (East China) in using COMSOL software.