Journal of Chemistry
Volume 2015 (2015), Article ID 815051, 13 pages
http://dx.doi.org/10.1155/2015/815051

Research Article

Finite Element Method Using a Characteristic-Based Split for Numerical Simulation of a Carbonate Fracture-Cave Reservoir

Liehui Zhang, Yuhui Zhou, Lei Zhao, and Deliang Zhang

State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, Sichuan 610500, China

Received 25 July 2014; Accepted 5 September 2014

Academic Editor: Jianchao Cai

Copyright © 2015 Liehui Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 by

where 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 is

where the quantities identified by ~ symbols indicate nodal values, and

where 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 expression

The matrices are

Step 3. The value of is determined from

where

3.2. Derivation Process of the Free-Flow Model

From (5), we derive

From 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 as

The first term on the right side of (19) can be expressed as

where

The second term on the right side of (19) can be expressed as

The third term on the right side of (19) can be expressed as

where

Letting

then

Similarly, (16) can be expressed as

Letting

we can derive the result of , and (5) becomes

According to CBS polymorphisms,

and combining (31) and (30) gives

The weak form of (32) is written using the characteristic-Galerkin approximation as

The first term on the left side of (33) can be expressed as

The second term on the left side of (33) can be expressed as

The third term on the left side of (33) can be expressed as

The fourth term on the left side of (33) can be expressed as

where

Thus, (32) can be simplified to

We use a semi-implicit scheme to solve the incompressible flow problem, so , .

At the completion of this stage, the value of is fully determined, and using the characteristic-Galerkin approximation, the weak form of (31) can be expressed as

The left side of (40) can be expressed as

The first term on the right side of (40) can be expressed as

The second term on the right side of (40) can be expressed as

The third term on the right side of (40) can be expressed as

The fourth term on the right side of (40) can be expressed as

where

The fifth term on the right side of (40) can be expressed as

where

Then, we derive

The values of and can be obtained from Steps 1 and 2, and can be calculated in Step 3 from (49). Then, can be substituted into Step 1 as an initial value for the next iterative computation.

4. Case Study and Discussion

An example calculation was performed for the initial condition of a cave with a reservoir pressure of 10 MPa, a fluid density of 800 Kg/m³, a fluid viscosity of 1 mPa s, and an axial inflow velocity of 0.01 m/s. The entrance was in the middle of the left boundary of the cave, and the exit was opposite to the entrance.

The boundary conditions are nonslip, so that . The grid at the cave boundary is refined. The number of nodes is 700, and the number of units is 1282. The influence of mass force is ignored, and the physical parameters of the flow in the reservoir are constant with time. Figure 2 shows a two-dimensional single-phase free-flow model and Delaunay triangle subdivision. Figure 3 shows a single-grid figure.

Figure 2: Two-dimensional single-phase free-flow model and Delaunay triangle subdivision.

Figure 3: A single-grid figure.

Based on the research results above, we used MATLAB to develop software for numerical simulation of two-dimensional single-phase free flow in a cave. The single-phase free flow of a large-scale fracture-cave reservoir can thus be simulated, with the simulation performed on any computer with MATLAB installed. The functions of the software include data input and storage, grid partitioning, simulation of single-phase flow, data output, and graphical plotting. Figures 4, 5, and 6 show outputs of the numerical simulation software.

Figure 4: Velocity field of the free-flow region at different times.

Figure 5: Velocity vector streamline diagram of the free-flow region at different times.

Figure 6: Pressure profile patterns of the free-flow region at different times.

Figures 4–6 represent the temporal evolution of the velocity field, velocity vector streamlines, and profile patterns of the free-flow region, respectively. The figures show that two vortices form at the two corners, beside the entrance to the cave, which is in agreement with the actual characteristics of free flow. Simulations indicate that the CBS algorithm has the advantage of flexibility of grid subdivision, and high efficiency and precision of calculations for the single-phase free-flow numerical simulation of fracture-cave carbonate reservoirs.

This work provides a theoretical basis for future studies on the FEM approach for multiphase flow and Darcy flow-free flow in a discrete media network model of a fracture-cave reservoir; this is an important development in the modeling of fracture-cave flow regimes and should be of considerable use in future research in this field.

5. Conclusions

Based on a single-phase free-flow numerical simulation of a fracture-cave reservoir, the CBS algorithm was adopted to obtain numerical solutions accurately and quickly, in which each step only requires solution of a system of linear equations. The velocity field, velocity vector streamline diagram, and profile pattern of the free-flow region at different times were determined, and they demonstrated good conformity with the actual characteristics of free flow. The results confirm that the CBS algorithm is applicable to the numerical simulation of a discrete media model of a fracture-cave reservoir.

Similarly to the traditional method, we can use statistical information and data on fractures and caves obtained from an oil field to construct a discrete fracture-cave geology geometry model. Then, we can use the approach in this study to model realistic single-phase cave free flow in fractured cave media. This paper provides an effective approach to developing a fracture-cave reservoir, including development plan design, dynamic prediction, reservoir evaluation, optimization of engineering technical problems, and decision-making.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Project is supported by the National Science Fund for Distinguished Young Scholars of China (Grant no. 51125019) and the Natural Science Foundation of China (Grant no. 51374181).

References

C. Pérez-Rosales, A. Mercado-Díaz, J. Cruz-Hernández, and R. Islas-Juárez, “Electric Resistivity of Vuggy Fractured Media,” in Proceedings of the SPE International Petroleum Conference and Exhibithion in Mexico, Paper SPE 74340, pp. 11–23, Villahermosa, Mexico, February 2002.

Z. Chen, Finite Element Methods and their Applications, vol. 5, Springer, Berlin, Germany, 2005.

V. Dalen, “Simplified finite-element models for reservoir flow problems,” Society of Petroleum Engineers Journal, vol. 19, no. 5, pp. 333–343, 1979.

J. G. Kim and M. D. Deo, “Finite element, discrete-fracture model for multiphase flow in porous media,” AIChE Journal, vol. 46, no. 6, pp. 1120–1130, 2000.

J. Cai, X. Hu, D. C. Standnes, and L. You, “An analytical model for spontaneous imbibition in fractal porous media including gravity,” Colloids and Surfaces A: Physicochemical and Engineering Aspects, vol. 414, pp. 228–233, 2012.

J. Cai, E. Perfect, C.-L. Cheng, and X. Hu, “Generalized modeling of spontaneous imbibition based on hagen-poiseuille flow in tortuous capillaries with variably shaped apertures,” Langmuir, vol. 30, no. 18, pp. 5142–5151, 2014.

Q. Liu, T. Stone, G. Han, R. Marsden, and G. Shaw, “Coupled stress and fluid flow using a finite element method in a commercial reservoir simulator,” in SPE Asia Pacific Oil and Gas Conference and Exhibition, Perth, Australia, October 2004, Paper SPE 88616.

M. Sheng, G. Li, S. N. Shah, and X. Jin, “Extended finite element modeling of multi-scale flow in fractured shale gas reservoirs,” in Proceedings of the SPE Annual Technical Conference and Exhibition (ATCE '12), Paper SPE 159919, San Antonio, Tex, USA, October 2012.

Y. Wei, T. V. Carlos, A. Ridvan, A. D. Saleh, A. T. Abdullah, and E. Haluk, “Interpretation of frequency-dependent dual-laterolog measurements acquired in Middle-East carbonate reservoirs using a second-order finite-element method,” in Proceedings of the 48th Annual Logging Symposium, Society of Petrophysicists and Well Log Analysts (SPWLA), Austin, Tex, USA, June 2007.

G. Qin, L. Bi, P. Popov, Y. Efendiev, and M. Espedal, “An efficient upscaling procedure based on Stokes-Brinkman model and discrete fracture network method for naturally fractured carbonate karst reservoirs,” in Proceedings of the International Oil and Gas Conference and Exhibition in China, Beijing, China, June 2010.

A. F. Gulbransen, V. L. Hauge, and K.-A. Lie, “A multiscale mixed finite-element method for vuggy and naturally-fractured reservoirs,” in Proceedings of the SPE Reservoir Simulation Symposium, Paper SPE 119104, pp. 898–908, The Woodlands, Tex, USA, February 2009.

M. S. Jouini and S. Vega, “Elastic properties estimation of carbonate samples using X-ray computed tomography images,” in Proceedings of the EAGE Workshop on Applications & Challenges of Rock Physics for Quantitative Geophysical Interpretation, Dubai, UAE, January 2012.

M. D. Jackson, A. H. Muggeridge, M. J. Blunt et al., “Reservoir modeling for flow simulation using surfaces adaptive unstructured meshes and control-volume-finite-element methods,” in Proceedings of the SPE Reservoir Simulation Symposium, Paper SPE 163633, The Woodlands, Tex, USA, February 2013.

Z. Chen, “An ABAQUS implementation of the XFEM for hydraulic fracture problems,” in Paper Presented at the ISRM International Conference for Effective and Sustainable Hydraulic Fracturing, 2013.

Y. Cinar, S. Rahman, and N. G. Doonechaly, “A new finite-element numerical model for analyzing transient pressure response of naturally-fractured reservoirs,” in Proceedings of the SPE Annual Technical Conference and Exhibition, New Orleans, La, USA, September-October 2013.

Z. J. Luo and F. Zeng, “Finite element numerical simulation of land subsidence and groundwater exploitation based on visco-elastic-plastic Biot's consolidation theory,” Journal of Hydrodynamics, vol. 23, no. 5, pp. 615–624, 2011.

W. Lai and A. A. Khan, “Time stepping in discontinuous Galerkin method,” Journal of Hydrodynamics Series B, vol. 25, no. 3, pp. 321–329, 2013.

Y. L. Zhou, H. W. Tang, and X. H. Liu, “A split-characteristic finite element model for 1-D unsteady flows,” Journal of Hydrodynamics, vol. 19, no. 1, pp. 54–61, 2007.

T. Liu, H. Xu, and X. Qiu, “A combination method of mixed multiscale finite-element and laplace transform for flow in a dual-permeability system,” ISRN Applied Mathematics, vol. 2012, Article ID 202893, 10 pages, 2012.

D. Gupta, L. Kumar, and B. Singh, “Finite element solution of unsteady mixed convection flow of micropolar fluid over a porous shrinking sheet,” The Scientific World Journal, vol. 2014, Article ID 362351, 11 pages, 2014.

J. Yao, Z. Huang, Y. Li, C. Wang, and X. Lv, “Discrete fracture-vug network model for modeling fluid flow in fractured vuggy porous media,” in Proceedings of the International Oil and Gas Conference and Exhibition in China, Beijing, China, June 2010, Paper SPE 130287.

Z. Chen, G. Huan, and Y. Ma, Computational Methods for Multiphase Flows in Porous Media, vol. 2, SIAM, Philadelphia, Pa, USA, 2006.

Z.-Y. Liu and H.-Z. Chen, “Discontinuous Galerkin immersed finite volume element method for anisotropic flow models in porous medium,” Abstract and Applied Analysis, vol. 2014, Article ID 520404, 10 pages, 2014.

Z. Li, T. Lin, and X. Wu, “New Cartesian grid methods for interface problems using the finite element formulation,” Numerische Mathematik, vol. 96, no. 1, pp. 61–98, 2003.

Z. Li, “The immersed interface method using a finite element formulation,” Applied Numerical Mathematics, vol. 27, no. 3, pp. 253–267, 1998.

Y. Gong, B. Li, and Z. Li, “Immersed-interface finite-element methods for elliptic interface problems with nonhomogeneous jump conditions,” SIAM Journal on Numerical Analysis, vol. 46, no. 1, pp. 472–495, 2007.

Y.-L. Zhao, L.-H. Zhang, J.-X. Luo, and B.-N. Zhang, “Performance of fractured horizontal well with stimulated reservoir volume in unconventional gas reservoir,” Journal of Hydrology, vol. 512, pp. 447–456, 2014.

J. Cai and B. Yu, “A discussion of the effect of tortuosity on the capillary imbibition in porous media,” Transport in Porous Media, vol. 89, no. 2, pp. 251–263, 2011.

J. Donea and A. Huerta, Finite Element Methods for Flow Problems, John Wiley & Sons, New York, NY, USA, 2003.

J. N. Reddy and D. K. Gartling, The Finite Element Method in Heat Transfer and Fluid Dynamics, CRC Press, New York, NY, USA, 2010.