Abstract

Fracture-cave carbonate reservoirs occur widely in source rocks and are prospects for exploitation worldwide. However, the presence of massive caves and multiscale fracture systems results in extremely complex fluid flow patterns. Therefore, in this paper, a discrete network model for fracture-cave reservoirs was established to study fluid flow characteristics and pressure distributions in complex flow regimes. In this study, the cave system was treated as a free-flow region, and the fluid flow in fracture systems followed the Navier-Stokes and Darcy equations, respectively. After discrete modeling, the Galerkin finite element method was used for numerical calculation of the single-phase free flow; the method maintains a high-precision result with low grid orientations during the simulation. In addition, because only one linear equation requires solving at each step, the solution is obtained quickly. Moreover, based on the proposed discrete media network model of fracture-cave reservoirs and the finite element numerical calculation method, a corresponding simulator was also developed. The finite element numerical simulation method based on the characteristic-based split (CBS) algorithm has proven to be applicable to complex flow problems in fracture-cave reservoirs.

1. Introduction

Fracture-cave reservoirs play a vital role in the petroleum industry, as these reservoirs contain large volumes of oil and are widely distributed in areas such as Saudi Arabia, Iran, Iraq, Mexico, Oman, and Syria [1]. Therefore, the development of reliable analytical solutions that consider the geometrical properties, storage capacities, and flow properties of porous media in fracture-cave reservoirs is of utmost importance. Many of the world’s giant fields contain naturally fractured and vuggy carbonate reservoirs, and the pore systems in these reservoirs are complex as a result of the characteristics of carbonate rocks, which are particularly sensitive to postdepositional diagenesis that can involve dissolution, dolomitization, and fracturing.

A great deal of research and model development has been devoted to flow behavior in fractured vuggy reservoirs; the focus of such research, which has been conducted worldwide, can be categorized into efforts to develop three types of models: continuum models, equivalent continuum models, and discrete medium models. Discrete medium models describe the heterogeneity of local-scale reservoirs better than the other models do; therefore, these models describe actual reservoir characteristics more accurately. Three numerical methods are commonly used to solve flow models: the finite difference method (FDM), the finite element method (FEM), and the boundary element method (BEM). The FEM has been widely used in the field of solid mechanics and is now also widely used in the petroleum research, thus providing an important approach to numerical simulation of problems involving discrete fractured networks [2, 3].

Kim and Deo [4] proposed a finite element discrete fracture model to describe multiphase flow in porous media. Their study showed that oil recovery is determined by a complex interplay of absolute matrix permeability, permeability contrasts, and flow rates; fracture capillary pressure also plays a significant role in determining water penetration into the matrix [5, 6].

Liu et al. [7] discussed the discretization of coupled finite element stress equations used in commercial reservoir simulations and compared the stress predictions of these equations with the results of a commonly used commercial stress program. Sheng et al. [8] used an extended FEM to investigate an integrated and improved shale-gas transport model and characterized the main flow mechanisms and the discrete fracture network. Wei et al. [9] presented a second-order finite element algorithm that can be used to simulate the frequency-dependent laterolog response of axially symmetric, invaded, and anisotropic formations. Qin et al. [10] developed a flow calculation method based on the Stokes-Brinkman model and the discrete fracture network method, which can be used for accurate and efficient upscaling of processes in naturally fractured carbonate karst reservoirs.

Gulbransen et al. [11] developed a multiscale mixed finite element (MsMFE) method that can be used to examine vuggy and naturally fractured reservoirs; the MsMFE is a first step toward developing a uniform multiscale multiphysics framework. In addition, a comparison of solutions of the MsMFE method and the fine-scale Stokes-Brinkman model for two different test cases, including both short- and long-range fractures, was also presented; the results demonstrate how fine-scale flow in fracture networks can be represented within a coarse-scale Darcy flow model using multiscale elements, with the solution computed using the Stokes-Brinkman equations. The results demonstrate that the MsMFE method can simulate flow in vuggy and naturally fractured reservoirs using a highly detailed geocellular approach.

Jouini and Vega [12] estimated the elastic properties of core plug samples from a Middle East carbonate reservoir and computed tomography scanner images by X-ray. Using the FEM, they simulated the elastic properties of the reservoir rocks. Using Hertzian’s and Hookean’s contact models, the discrete element method can be used to investigate and predict the elastic properties of reservoir rocks.

Jackson et al. [13] developed approaches to simulate fluid flow using the pillar-grid concept. Results showed significant improvements in the representation of multiscale geological heterogeneities and predictions of flow through those heterogeneities. In their research, which was based on more than 20 years of development of innovative numerical methods in ocean modeling, Jackson et al. refined and modified their approach to address the unique challenges associated with reservoir simulation. Chen [14] presented a new finite element in ABAQUS that incorporates the extended FEM (XFEM) in the solution of hydraulic fracture problems.

Cinar et al. [15] developed a numerical model to analyze transient pressure data from naturally fractured reservoirs; their model calculated well pressures in a single well in a naturally fractured reservoir. The characteristic properties of naturally fractured reservoirs, including fracture orientation, fracture length, fracture density, and fracture distribution, were also simulated in this model. The simulation results showed that the model can predict the transient pressure response in complex naturally fractured reservoirs; in addition, it can be used to match pressure data and diagnose characteristic properties of naturally fractured reservoirs. In addition to the studies reported above, a great deal of other research related to the FEM has been carried out [16–22].

In conclusion, the FEM has been widely studied and applied in petroleum engineering, especially for the analysis of flow in sandstone reservoirs and naturally fractured reservoirs, as well as in shale-gas reservoirs; however, the technique has rarely been applied to fracture-cave reservoirs. Because of the lack of a systematic theory and appropriate modeling techniques, research on theories and techniques of numerical simulation using discrete medium models to model giant fracture-cave reservoirs is still in an exploratory stage, with existing studies concentrating mostly on fractured media. The characterization of vuggy media with extremely complex internal geometries using theoretical models remains difficult. The characteristic-based split (CBS) algorithm has been used to solve problems in solid mechanics and fluid mechanics, but until now the algorithm has not been applied to a discrete media model of fracture-cave reservoirs. In this paper, we applied the FEM to the study of fracture-cave reservoirs based on a CBS discrete medium model.

2. Discrete Fracture-Cave Network Model Concept and Mathematical Formulation

As observed in numerous studies on the field of carbonate formation, three types of porosity are typically present in natural fracture-cave reservoirs; these are matrix porosity, fracture porosity, and cave porosity [23–26]. Fractures and caves are irregularly distributed and vary in size from microscopic to macroscopic. Figure 1 presents the conceptualized discrete fracture-cave network model. The isolated caves are connected via a discrete fracture network. Flow within the matrix and in macroscopic fracture systems follows the Darcy law and in cave systems follows a free-flow region.

Figure 1 

Discrete fracture-cave network model and simplified cave model.

In this paper, we assume that fluid flow is isothermal, single phase, and slightly compressible and that the fluid viscosity is constant. For a laminar viscous fluid flowing through an open domain (e.g., a cave system), the flow is governed by the continuity and momentum equations, which can be expressed in vector form as [27, 28]where is the fluid density (kg/m³), is the fluid velocity (m/s), is a gradient operator, is the body force per unit mass (m/s²), and is the fluid stress tensor (Pa). is given bywhere is the pressure in the fracture system (Pa), is a unit tensor, is the fluid viscosity (Pa s), and , the strain rate, is

For simplicity, in this study we treat the fluid as a Newtonian fluid. The density of a slightly compressible fluid can be assumed to be nearly constant (i.e., incompressible). Therefore, the governing fluid motion equations can be reduced to the Navier-Stokes equations, which are

The boundary conditions for (2) and (3) can be classified into two types: Dirichlet (velocity) and Neumann (traction) [29, 30], which can be expressed, respectively, as where is the specific velocity at the Dirichlet boundary , is an outward unit normal vector at the boundary, and is the specific traction.

3. Derivation of the Finite Element Method Based on CBS

3.1. CBS Procedure

The theory of the CBS algorithm is based on the split- and characteristic-Galerkin procedures, which are extended in the FEM. The CBS algorithm involves three indispensable steps: Step 1 is to compute the intermediate momentum variable; Step 2 is to calculate the density or pressure field; and Step 3 involves combining the introduction of the density or pressure fields to obtain the correct momentum variables.

Step 1. The solution for in matrix form iswhere the quantities identified by ~ symbols indicate nodal values, andwhere is a matrix, is a deformation function, , and . After integration by parts, the expressions for and are

Step 2. The value of can be easily determined using the expressionThe matrices are

Step 3. The value of is determined fromwhere

3.2. Derivation Process of the Free-Flow Model

From (5), we deriveFrom Step 1, after the split-temporal discretization, (15) becomes (omitting the third-order derivative above)From Step 2, after spatial discretization, we find that

Then,The left side of (19) can be expressed asThe first term on the right side of (19) can be expressed aswhereThe second term on the right side of (19) can be expressed asThe third term on the right side of (19) can be expressed aswhereLettingthenSimilarly, (16) can be expressed asLettingwe can derive the result of , and (5) becomes

According to CBS polymorphisms, and combining (31) and (30) givesThe weak form of (32) is written using the characteristic-Galerkin approximation asThe first term on the left side of (33) can be expressed asThe second term on the left side of (33) can be expressed asThe third term on the left side of (33) can be expressed asThe fourth term on the left side of (33) can be expressed aswhereThus, (32) can be simplified to