2 D COMPOSITIONAL RESERVOIR SIMULATION USING UNSTRUCTURED GRIDS IN HETEROGENEOUS RESERVOIRS

This paper presents a 2D implementation of an element-based finite-volume formulation (EbFVM) using triangles, quadrilaterals, or a combination of both elements in conjunction with a compositional reservoir simulator. The formulation is implemented into the UTCOMP simulator considering a full tensor formulation for the advection and diffusion terms. The UTCOMP simulator was developed at the Center of Petroleum and Geosystems Engineering at The University of Texas at Austin for the simulation of enhanced recovery processes. UTCOMP is a compositional multiphase/multi-component, IMPEC simulator, which can handle the simulation of miscible gas flooding processes. The results of several case studies involving homogeneous and heterogeneous reservoirs are presented in terms of volumetric rates of oil and gas, and saturation fields.

There are some specific differences between the several control volume finite-element methods available in the literature.For instance, the absolute permeabilities and porosities are stored in the vertices of the elements in the approach used by [8][9]21].Therefore, in order to perform the flux evaluation, it is necessary to solve a local linear system to guarantee the flux continuity.The authors call the approach multi-point flux approximation.In most CVFEM (Control Volume Finite-Element Method) approaches used in petroleum reservoir simulation, the approximate equations for multiphase fluid flow are obtained first for a single-phase flow and then the transmissibilities are multiplied by mobilities and densities to obtain the equations for the multiphase flow.Fluid flow and heat problems were first solved using the CVFEM in [3] and [20] in conjunction with triangles and quadrilaterals, respectively.The ideas of [3] and [4] were applied to water flooding problems in [4][5][6].These authors show that if the equations are obtained from a single-phase flow equation and then transmissibilities are multiplied by densities and phase mobilities, they do not correctly approximate the equations for multiphase flow.They defined their approach as Element based Finite Volume Method (EbFVM).In this work, we will also adopt such a denomination.As explained in [17], we have a methodology that still follows the conservative principles at the discrete level and only borrows the idea of elements and shape functions from the finite element method.Compositional, multiphase, multi-component fluid-flow problems in conjunction with anisotropic and heterogeneous reservoirs were solved by EbFVM using a fully implicit procedure as in [18][19].In this work, we apply the ideas of EbFVM to an in-house simulator called UTCOMP.UTCOMP was developed at the Center for Petroleum and Geosystems Engineering at The University of Texas at Austin for the simulation of enhanced recovery processes.UTCOMP is an IMPEC (Implicit Pressure Explicit Composition), multiphase/multi-component simulator, which can handle the simulation of several enhanced oil recovery processes [1,2].

PHYSICAL MODEL
In this section, we present the formulation for discretization of the partial differential equations arising from the formulation used by the UTCOMP simulator.Then, the procedure to obtain the approximate equations using the EbFVM is presented.

Component Molar Balance
The material balance in terms of number of moles of each component present in the reservoir is given by (1) where D is the depth, which is assumed to be positive in the downward direction,  j and  j are the molar density and mobility of the j-th phase, respectively.x ij is the molar fraction of 1 0 for 1, 2, , , is molar weight of phase j, q i is the molar rate of component i, and V b is the bulk volume.and are the absolute permeability and dispersion tensor, respectively.The and components in Cartesian coordinates are given by ( The components are functions of molecular diffusion and mechanical dispersion.The expressions for each one of components can be found in [1].

Phase-Equilibrium Relationship
The phase-equilibrium calculation determines the number, amounts and composition of all equilibrium phases.The equilibrium solution must satisfy three conditions.First, the molar-balance constraint must be preserved.Second, the chemical potentials for each component must be the same in all phases.Third, the Gibbs free energy at constant temperature and pressure must be a minimum.The first partial derivative of the total Gibbs free energy with respect to the independent variables gives equality of component fugacities among all phases: (4) The phase composition constraints (5) and the equations for determining the phase amounts for the two hydrocarbon phases is given by (6) Equations ( 5)-( 6) are implicitly used in the solution of the fugacity equation, Eq. ( 4), where is the ratio of moles of gas to total moles. . .

Volume Constraint
The volume constraint states that the pore volume in each volume must be filled up completely by the total fluid volume as described below: (7) where L j is a ratio of moles in phase j to the total number of moles in the mixture, is the molar volume of phase j, and v p is the pore volume.

Pressure Equation
In the UTCOMP simulator an IMPES approach is used.In this approach, the pressure equation is solved implicitly and conservation equations (Eq. 1) are solved explicitly.The pressure equation is obtained based on the fact that the pore volume should be completely filled by the total volume of fluid: (8) where the total volume of fluid is assumed to be a function of pressure and total number of moles of each component, and the pore volume is related to pressure only.Taking the derivative of the left-hand side of Eq. ( 8) with respect to pressure and number of moles and the right-hand side with respect to pressure and substituting Eq. ( 1) in the resulting equation, the following equation for pressure is obtained: (9) where denotes the derivative of total volume of fluid relate to N i .

Material Balance Equation
In the EbFVM, each element is divided into sub-elements.These sub-elements will be called sub-control volumes.The conservation equation and pressure equation, Eq. ( 1) and ( 9), respectively need to be integrated for each one of these sub-control volumes.For a two dimensional discretization there are two possible elements to use (triangular and quadrilateral elements).Figure 1 presents a triangular and a quadrilateral element, and all of the sub-control volumes associated with each element.Integrating Eq. ( 1) in time and for each one of the subcontrol volumes, and applying the Gauss theorem for the advective and diffusion terms, we obtain: (10) To evaluate the first and second terms of Eq. ( 10), it is necessary to define the shape functions.For triangles and quadrilaterals, linear and bi-linear shape functions as defined through Eqs. ( 11) and (12), respectively, will be used.Using the shape functions, any physical properties or positions can be evaluated inside an element as (13) where N v denotes the number of vertex for each element.Elements using the same shape function for coordinates and physical properties are known as isoparametric elements [16].Using the shape functions, gradients of potentials or any other variable can be easily evaluated as (14) To evaluate the gradients, it is necessary to obtain the derivatives of shape functions relative to x and y.These derivatives are given by 1 , , , ( ) 0 for 1, 2, , .
where J t is the Jacobian matrix of the transformation, and it is given by ( 16) To perform the integral of Eq. ( 10), it is necessary to define the volumes of each subcontrol volume and the area of each interface.The volumes of each sub-control for triangles and quadrilaterals, respectively, are given by ( 17) (18) where h is the thickness of the reservoir.For a quadrilateral element, det(J t ) needs to be evaluated at the center of each sub-control volume.The area of each interface, reading counterclockwise, is given by (19) Substituting Eqs.(17) and (18) for the accumulation term, and (19) for the advective and diffusive fluxes into Eq.( 10), evaluating the fluid properties through an explicit procedure, and dividing each term by the time-step (t) the following equations for the two mentioned terms are obtained: In Equations ( 20) and ( 21), the superscript "n" and "n+1" means physical properties from the previous time-step and current time-step, respectively.By inspecting Eq. ( 21), it can be inferred that it is necessary to evaluate molar densities, molar fraction, mobilities, saturations and porosities in two interfaces of each sub-control volume.In the approach used in this work, each element has a constant porosity and absolute permeability tensor, although these properties can change from one element to another.The other physical parameters are det( ) .
based on the vertices of the elements.Therefore, in order to evaluate these properties, an upwind scheme will be used.The mobilities and other fluid properties are evaluated at integration point 1 of Fig. 1a, for instance, by Inserting Eqs. ( 20) and ( 21) into Eq.( 10), the following equation for each element is obtained: (23)

Pressure Equation
Integrating Eq. ( 9) in time and for each one of the sub-control volumes shown in Fig. 1, then applying the Gauss theorem for the advective and diffusion terms, we obtain: (24) Substituting Eq. ( 19) for the advective and diffusive fluxes into Eq.( 24) and evaluating the fluid properties in each interface of the sub-control volumes through an explicit procedure, the following equation for pressure is rendered for the accumulation, advective and diffusive, and sink/source terms: (25) where A ccp , m , F m , and B m denote the accumulation, advective and diffusive, and sink/source terms, respectively.The approximate expressions for each one of these terms are given below: .
Equations ( 23) and ( 25) denote the molar conservation of i-th component and the approximate equation for the evaluation of pressure field, for each sub-element.Now, it is necessary to assemble the equation of each control volume obtaining the contribution of each sub-control volume that shares the same vertex.This process is similar to the assembling of the stiffness global matrix in the finite element method.Figure 2 presents a control-volume around vertex 5 (dark continuous lines) that will receive contributions from scv1 of element 1, scv3 of element 2, scv4 of element 6, and svc1 of element 7. It is important to mention that each control-volume equation can have different permeabilities and porosities.

RESULTS AND DISCUSSION
In this section, we will present the three cases studies in conjunction with the EbFVM formulation.The first case refers to a three component gas injection in a quarter-of-five-spot configuration considering an isotropic and homogeneous reservoir.All reservoir data used for this case are shown in Table 1.The Corey's model relative permeabilities parameters are presented in Table 2.One of the unstructured meshes used in this case is shown in Figure 3, where the blue and red spots are the injection and production wells, respectively.The oil and gas production rates curves are shown in Figure 4 and Figure 5, respectively.From these figures, a good agreement can be verified between the unstructured investigated meshes and the refined Cartesian mesh.It is also important to mention that the number of volumes used by the unstructured grids in order to obtain the same accuracy as the refined Cartesian mesh (100x100) is much smaller.The gas saturation field at 500 days obtained with the hybrid and refined Cartesian mesh are presented in Figure 6.Although the hybrid mesh is much more coarse than the Cartesian one, we can observe a good agreement in the gas saturation front.
Figures 8 and 9 present the volumetric rates of oil and gas at standard surface, respectively.From these figures, it can be seen that a good agreement between the results obtained with the two unstructured and the two Cartesian meshes.The main goal of this case study was to verify the implementation of the dispersion tensor for unstructured meshes.It is worthwhile to mention that each element has a constant dispersion tensor; therefore, no interpolation is necessary to evaluate the dispersion fluxes along each integration point.Figure 10 shows the oil saturation field in two simulated times.Also the results for the more refined Cartesian mesh are presented.Once again, it is possible to verify a good agreement between the results produced by the two different approaches.The third case refers to the injection of immiscible gas in a saturated reservoir with irregular borders.Figure 11 presents one the employed meshes as well the position of the four wells (two injectors and two producers).The reservoir data, components initial composition and fluid injection composition are shown in Table 5, and Stone's 2 model are shown in Table 6.Due to the difficulty in adjusting a Cartesian mesh to the irregular shape of the reservoir, the results for this case are not presented for Cartesian meshes.Figures 12 through 14 present the total oil and gas rates at standard surfaces, respectively, and the average reservoir pressure.From these figures, a good agreement can be observed between all the reported results.

CONCLUSIONS
The implementation of an EbFVM in conjunction with unstructured grids considering the advection and diffusion transport terms into the UTCOMP simulator were carried out in this work.The results of select case studies were compared to the original formulation of the UTCOMP simulator using Cartesian meshes and a very good agreement between the results were obtained.

Figure 1 .
Figure 1.Triangular and quadrilateral elements and their respective sub-control volumes.