A TVD SCHEME FOR 3 D UNSTRUCTURED GRIDS APPLIED TO COMPOSITIONAL RESERVOIR SIMULATION

In reducing the grid orientation effect for the numerical solution of partial differential equations, interpolation functions play an important role when the advective transport of the governing equations is considered. This is due to the fact that, in general, the unknowns are evaluated in the vertices of the elements and such properties must be extrapolated to inner parts of the elements. First-order schemes, such as upwind, are the easiest methods to use for performing the extrapolation of the properties. However, such methods introduce a large amount of numerical diffusion in the solution. A few higher-order interpolation schemes, on the other hand, are capable of providing solutions free of numerical diffusion, increasing the accuracy of the method and reducing the computational efforts required. In this work, we investigate the TVD interpolation scheme for three-dimensional unstructured grids in conjunction with Element-based Finite Volume Method (EbFVM) using four types of elements: hexahedron, tetrahedron, prism and pyramid.


INTRODUCTION
Interpolation schemes play an important role when numerical methods such as finite-volume method are used.Usually, most of the unknowns are evaluated in the center of each grid block, for structured grids, or in each vertex, for cell-vertex grids.However, in order to compute advective and diffusive terms it is necessary to extrapolate these variables from the cell's center or vertex of the grid to the interfaces of the grid blocks.This is performed through the use of interpolation functions.Unfortunately, accurate interpolation schemes cannot be used arbitrarily for advective dominant problems.The use of high-order schemes for the advective dominant problems yields spurious oscillations in the solution.In order to overcome this limitation, the use of a Total Variation Diminishing (TVD) scheme can be an alternative option.It is known that TVD schemes can preserve the monotonicity of the solution, hence avoiding spurious oscillations.
Higher-order TVD schemes usually blend the stability of the upwind schemes with the use of higher-order terms through the use of flux-limiters.Van Leer (1979) developed the first monotonicity preserving flux-limiter.However, the importance of the TVD in the monotonicity preserving was only discussed later by Harten (1983).Further, Sweby (1984) developed the TVD region for one-dimensional problems.Several flux-limiters were later introduced and studied by other authors, such as Roe (1983), Chakravarthy and Osher (1983), Roe (1986) and Koren (1993).
However, the application of these functions for unstructured grids was not straightforward.The Controlvolume Finite Element Method (CVFEM) presented by Baliga and Patankar (1980) was used in conjunction with lower-order schemes; they used upwind and exponential schemes.This exponential function, also known as Flow-Oriented Interpolation (FLO), performs the interpolation based on the exponential element shape functions for the element's average velocity and linear for the normal direction.Prakash (1986) added the source terms evaluation to the original FLO function, giving rise to the FLOS scheme.Later, Maliska (2004) renamed the CVFEM as Element-based Finite Volume Method (EbFVM), since we still have a conservative approach that borrows the idea of elements and shape functions from the finite-element method.New interpolation functions were introduced for the EbFVM, mainly for the simulation of skew grids.Hassam et al. (1983) and Schneider and Raw (1986) were the first authors to propose skew upwind schemes for twodimensional problems.Also, Swaminathan andVoller (1992a, 1992b) adapted the streamline upwind Petrov-Garlekin (SUPG) function developed by Brooks and Hughes (1995) from Cartesian to unstructured grids using EbFVM, presenting the streamline upwind control volume (SUCV).Many of those schemes may produce results with second-order accuracy in space, however at the cost of low numerical stability and possibly downwind influence, which requires a proper treatment, and in general, they are reduced to first order schemes.
Various TVD schemes have been studied and implemented for unstructured grids, such as Bruner and Walters (1995), and Hou et al. (2011, 2012).However, none of those works aimed at the application of TVD in cell-vertex approaches like as EbFVM.Just recently, Fernandes et al. (2013), based on the work of Darwish and Moukalled (2003), introduced a TVD scheme for 2D unstructured grids for the EbFVM in conjunction with compositional reservoir simulation.Later, Fernandes et al. (2015a) extended the method for solving three-dimensional problems using Hexahedron element grids.In this work, we extend the work presented by Fernandes et al. (2015a) to tetrahedron, prism, and pyramid elements as well as hybrid element grids.This implementation is performed in the in-house compositional reservoir simulator UTCOMP (1990).UTCOMP is an IMPEC-based (Implicit Pressure Explicit Compositional) compositional, multiphase/ multicomponent simulator that can take into account up to four phases developed at The University of Texas at Austin.The EbFVM implemented in the UTCOMP simulator are based on Marcondes and Sepehrnoori (2010), Marcondes et al. (2013) and Santos et al. (2013) for the General Purpose Adaptive Reservoir Simulator (GPAS).In UTCOMP, the implementation for 2D and 3D unstructured grids using the EbFVM and the IMPEC approach are presented in Fernandes et al. (2012) and Araújo et al. (2016).Additionally, IMPSAT (Implicit Pressure and Saturations) approaches (Fernandes et al., 2014;Fernandes et. al., 2015b) and two volume balance based fully implicit approaches (Fernandes, 2014;Fernandes et al., 2016) were implemented for both Cartesian and unstructured grids.In this manuscript, only IMPEC will be considered.

MATHEMATICAL MODEL
The isothermal fluid flow in porous media can be modeled through a material balance equation for each chemical species and constraint equations.The local equilibrium is assumed in order to model the hydrocarbon properties, amounts and compositions (Chang, 1990).
The material balance equation of each component considering both advection and diffusion transport and assuming multiphase Darcy's flow is written as where n c represents the number of hydrocarbon components, n c +1 refers to the water component, n p is the number of fluid phases, ij Κ is the physical dispersion tensor, k is the absolute permeability tensor, N i is the number of moles of the i-th component, ξ j and λ j represent, respectively, the molar density and the mobility of the j-th phase, φ is the formation porosity, x ij represents the mole fraction of the i-th component in the j-th phase, q i is the molar flow rate of the i-th component of a grid block that contains a well, and V b stands for the bulk volume of the grid block.The hydraulic potential of a phase, Φ j , is given by Z P j j j γ − = Φ where P j e j γ are the pressure and specific gravity of the phase j, respectively, and Z denotes the depth, which is positive in the downward direction.
Local phase equilibrium is assumed between all hydrocarbon phases in each grid block.In this work, grid block refers to the control volume that is assembled around each vertex of the grid.This condition is represented by the (1) (2) phase fugacities (f), as follows: 0 ; 1,...., ; 2,....., where the superscript r describes the reference phase, and fugacity is a function described by ( ) ln ; 1,..., ; 2,..., In Eq. ( 4), ∅ ij stands for the fugacity coefficient of the i-th component in the j-th phase.The mole fraction constraint equation, for each hydrocarbon phase, and the Rachford-Rice equation are given by where z i and K i are the overall mole fraction and the equilibrium ratio of the component i, respectively, and υ is the gas mole fraction in the absence of water.
The calculations associated with phase equilibrium and PVT properties are conducted using the Peng-Robinson Equation of State (Peng and Robinson, 1978).It is assumed that there is no mass transfer between the water and hydrocarbon phases.Therefore, water is not included in the flash calculations.Further details about this calculation can be found in Perschke (1988).
The last equation is the volume constraint equation, which is given by where S j denotes the saturation of the phase j.
The UTCOMP simulator is based in an IMPEC (Implicit Pressure Explicit Composition) formulation as described by Àcs et al. (1985).The unknown primary variables are the oil pressure and the total number of moles of each component.The oil pressure is obtained from a volume balance, giving rise to the following expression: The applied formulation evaluates the oil pressure implicitly, while all the remaining properties are explicitly evaluated from the previous time step.The mole balance and flash calculations, for instance, are only conducted after the pressure at the new time is evaluated.

APPROXIMATE EQUATIONS
In the EbFVM, the physical domain is divided into elements, which are sub-divided into sub-elements according to the number of vertices.These sub-elements are called sub-control volumes, since the material balance equations are integrated to each one of these sub-elements.Figure 1 shows the four elements investigated in this work.For each element, a local coordinate system is established.This allows the simulator to perform the flow rate at each interface of the sub-control volume, etc. Regardless of how distorted the element is.It is important to stress that except for the pyramid element, the sub-control volumes of the other three elements always have three quadrilateral integration surfaces, where the material balance is evaluated; the sub-control volumes of the pyramid element associated to the base have also three integration surfaces, but they are triangular surfaces, and the sub-control volume associated to apex of the pyramid has 4 quadrilateral integration surfaces.
In order to obtain the approximate balance equation for each component and pressure, Eqs. ( 1) and ( 7) are integrated in time and for each sub-control of each element.Integrating, for instance, Eq. (1), in space and time and applying the Gauss theorem for the advective and dispersion terms, we obtain: In order to perform the integrations of the first and second terms of Eq. ( 8), it is necessary to make use of the shape functions.The shape functions, for hexahedron, prism, tetrahedron, and pyramid elements, are respectively given by Eqs. ( 9) through ( 12).
Using the shape functions, the area, volume, and gradients that are necessary to evaluate each term of Eq. ( 8), can be easily evaluated.Further details of this calculation can be found in Marcondes et al. (2013) and Santos et al. (2013).Adopting an explicit scheme and evaluating the physical properties at each interface of the control volume, the following expressions for the first and second terms of Eq. ( 8) are obtained: where ip denotes each integration point located on the surfaces of each sub-control volume and Nv is the number of vertices of each element.A similar procedure is performed for the pressure equation.Substituting Eqs. ( 13) and ( 14) in Eq. ( 8), the following equation for each subcontrol volume of each element is obtained: , , 0 ; 1,..., ; 1,..., 1 The above equation stands for the material conservation for each sub-control volume of each element.The total material balance of each grid volume is computed assembling the contributions of all sub-control volumes that share the same vertex.This procedure is described in detail in Marcondes et al. (2013).
In order to evaluate each term of Eq. ( 14), it is necessary to calculate physical properties, such as mole fractions, phase molar densities, and mobilities at the integration points of each sub-control volume.This extrapolation is conducted through interpolation functions.In this work, two interpolation schemes are implemented for the four element types shown before: a first order upwind, and a second order TVD in conjunction with two flux-limiters.
The upwind interpolation for the EbFVM is analogous to the classic upwind used for Cartesian grids.Considering, for instance, the integration point 1, for all elements shown in Figure 1, the product of physical properties is evaluated by In Eq. ( 16), the normal to the integration surface is always positive when it is orientated in the outward direction of that interface.
For the TVD function, we will investigate two fluxlimiters to switch between lower and higher-order schemes.Although the scheme proposed by Fernandes et al. (2013) for approximating the successive slope ratio has been shown to be efficient for 2D unstructured grids, results obtained by Fernandes et al. (2014) demonstrated that, when the gravitational effects are added this scheme is more prone to spurious oscillations.Therefore, in this work we used the successive slope ratio proposed by Darwish and Moukalled (2003), that performed better than the one used by Fernandes et al. (2013) and Fernandes et al. (2014).According to Darwish and Moukalled (2003), the interpolation of any physical property for the EbFVM can be evaluated through the following equation: In Eq. ( 17), the subscript f denotes the integration point where the property needs to be evaluated, U denotes the property evaluated in the upwind vertex, D denotes the property evaluated in the downwind vertex, and r f is the successive slope ratio.In this investigation, we will employ the following equation given by Darwish and Moukalled (2003) to evaluate the slope ratio: , where U ∇∆  represents the property's gradient evaluated at the upwind vertex, and UD r ∆  is the vector distance from the upwind node to the downwind node.The gradient at the vertex is evaluated by the volumetric average of the gradient of all elements that share the same vertex (Tran et al., 2006).
Finally, the two flux-limiters investigated in this work for the TVD interpolation function are the MINMOD (Roe, 1986), or MM, as we shall call it from now on, the most diffusive limiter, and the Koren (1993), a very compressive limiter.Their expressions are given, respectively, by MINMOD: .

RESULTS
Three case studies will be presented to investigate the current implementation.The first is an advectiondispersion tracer injection in a quarter of five-spot, the second case study is a CO 2 flooding, and the third case is a CO 2 injection in an irregular reservoir.The accuracy of the TVD schemes is presented for each element individually, although the method presented is implemented in a general way and hybrid grids composed of hexahedron, tetrahedron, pyramid or prism can be considered without any further consideration.
In the first case, a non-reactive tracer is injected in a quarter of five-spot.Water is the only mobile phase and the tracer does not partition in any other phase other than water.Tracer is injected along with water until a pore volume injected (PVI) of 0.02.After 0.02 PVI, just water is injected.Gravitational effects are neglected in order to validate the solutions with the analytical solution given by Abbaszadeh-Dehghani and Brigham (1984).All data for this case are presented in Table 1.A quarter-of-five spot configuration is used with two vertical wells in the corners.The grids are formed dividing the reservoir as a Cartesian grid.For elements other than hexahedron, each grid cell is divided again to form the other elements, such as prisms, pyramids, and tetrahedrons.
The tracer's effluent profile is presented in Figure 2, where two grids were used to compare each interpolation function to the exact solution.The results obtained for hexahedron, tetrahedron, prism, and pyramid elements are presented Figures 2a and 2b, 2c and 2d, 2e and 2f, and 2g and 2h, respectively.As expected, the Koren's TVD fluxlimiter is always the most accurate while the upwind is less accurate.For the 50x50 grid using hexahedron elements, the results for the Koren's flux-limiter were almost in good agreement with the exact solution.For the grids with tetrahedron elements, the results of the TVD, although accurate, but were not as accurate as for the other elements.One possible reason for this is the grid orientation effect.Also, due to the shape functions the tetrahedron is the less accurate element.As can be observed from Figure 2f, the result obtained with Koren's flux limiter matches the exact solution.Results presented in Fig. 2h were quite good for pyramid elements as well.
In the second case study gravity effects are considered.In this case, CO 2 is injected in a quarter of five-spot configuration.We illustrate this case by running a coarse and a fine grid for hexahedrons, tetrahedrons, and prism elements.The reservoir data and fluid properties for this case are presented in Table 2. Corey's model (Corey, 1986) is used for modeling the multiphase relative permeability curves.
The oil and gas volumetric production rate curves are presented in Figure 3.It can be observed that, for the hexahedron and prism element grids, the solution for the upwind scheme converges to the higher resolution schemes as the grid is refined for both Koren and MINMOD (MM) flux limiters.However, for the tetrahedron element, it can be seen that refining the grid makes the solution for the upwind scheme deviate from the TVD schemes investigated.The reason for such behavior is that all grids employed in this case study have the same number of vertices in the Z-direction.Since the tetrahedron element is the less accurate element and few vertices in the Z-direction are commonly used in reservoir simulation, the upwind solutions for the tetrahedron elements have a poor resolution in the Z-direction.As we use the TVD schemes, the resolution in the Z-direction is highly increased and the curves converge to the correct answer.This argument was confirmed by running more refined grids in the Z-direction, which resulted in an upwind solution closer to that of the TVD solution.
The gas saturation fields at 260 days for the coarser grids are presented in Figure 3, where it can observed that the gas saturation front is more dispersed when using the upwind scheme, and less dispersed for MM and Koren's flux limiters.This shows how effective the TVD schemes are in reducing the effects of the numerical dispersion.
The third Case is similar to Case 2, but now an irregular reservoir is considered.The reservoir geometry is presented in Figure 5.We used the same data set shown in Table 2, except for the reservoir dimensions.This case shows that a better accuracy can still be obtained when using TVD even if the grid elements are distorted or vary in size.We present the gas production curves in Figure 6.From Figure 6, it can be clearly observed that the upwind curves tend to that obtained by using Koren's flux limiter.It also indicates that the TVD scheme can produce accurate results even for irregular grids.
Figure 6 presents a comparison of the gas saturation field between the upwind and TVD scheme using a Koren flux limiter, in a plane cut between two of the injecting wells, for the coarser grids at 500 days.This view was chosen for showing both gravity and displacement effects of the gas saturation front.From this figure, once again it can be observed that the fields obtained with the TVD function present less numerical dispersion than those obtained with the upwind function.The TVD scheme was not only capable to capture a sharper gas front, but also shows that the gas is expected to migrate to the higher reservoir layers.Such behavior is in accordance to that observed when the grid is refined in the Z-direction, showing that the TVD schemes also enhanced the accuracy of the gravitational effects.

CONCLUSIONS
Two 3D TVD interpolation functions (MINMOD and Koren) in conjunction with unstructured grids were implemented in a compositional reservoir simulator using the Element based Finite-Volume Method considering four element types: hexahedron, tetrahedron, prism, and pyramids and two-flux limiters.Comparison of the results using both TVD schemes and the upwind scheme shows that TVD schemes produce more accurate nonoscillatory results.In addition, the EbFVM in conjunction with the higher-order TVD limiters can easily be applied to irregular grids in order to model irregularly shaped complex reservoirs.Also, the proposed method (EbFVM with higher-order TVD limiters) allows the usage of coarse grids, while maintaining high accuracy for the computational results.

V
Total fluid volume (m 3 ) ti V Total fluid Partial molar volume (m 3 /mol) i V Phase partial molar volume (m 3 /mol)

Table 1 .
Fluid and reservoir data for Case 1.