A comparative study of implicit and explicit methods using unstructured voronoi meshes in petroleum reservoir simulation

A comparative study of implicit and explicit methods using unstructured voronoi meshes in petroleum reservoir simulation

Francisco Marcondes^I; Clovis R. Maliska^II; Mário C. Zambaldi^III

^IMember, ABCM, [e-mail], Federal University of Ceará - UFC Dep. of Metallurgical Eng. and Material Sci. 60341-790 Fortaleza, CE, Brazil

^[e-mail], Federal University of Santa Catarina - UFSC Department of Mechanical Engineering 88040-900 Florianópolis, SC, Brazil

^[e-mail], Federal University of Santa Catarina - UFSC Department of Mathematics 88040-900 Florianópolis, SC, Brazil

ABSTRACT

This work presents a comparative study among three linearization schemes widely used in petroleum reservoir simulation, namely the Implicit Pressure Explicit Saturation (IMPES), Fully Implicit (FI) and Adaptive Implicit Method (AIM). Attention is given to the criterion used for switching from IMPES to FI and vice-versa. The black-oil model considering the oil-water flow is employed. The equations are discretized using an unstructured Voronoi mesh in a finite-volume framework. The numerical results are analyzed based on the number of iterations in the Newton and solver procedures, number of implicit volumes, and mass balance error of components. The effect of the stopping criterion on Newton's iterations is also investigated. The outcome of the study reveals that adaptive implicit methodologies can be a good choice when unstructured grids are present, allowing the use of time steps with the same order of magnitude as the ones used in FI procedures, but with reduction in the computational effort.

Keywords: adaptive implicit method, voronoi grids, finite-volume method, petroleum reservoir simulation

Introduction

IMPES is one of the most widely used procedures in petroleum reservoir simulation. In this scheme the mobilities are evaluated with saturation and pressure from the previous time level. This fact decouples the pressure equation from the saturation equations, allowing the saturation to be calculated explicitly, while pressure is kept implicit. Such a scheme requires a small computational effort per timestep, because the pressure is the only unknown to be calculated through a linear system of equations. Another advantage is that the procedure for advancing the saturation is easily vectorized. The main disadvantage of IMPES is that the Courant-Friedrichs-Lewy (CFL) number must be smaller than unity, to avoid spurious oscillations in the solution. On the other hand, a methodology which works well with CFL larger than unity is the Fully Implicit (FI) one, which solves pressure and saturation simultaneously. However, the computational effort per timestep increases considerably when compared with IMPES. According to Thomas and Thurnau (1983), except in some critical regions of the reservoir, such as in the vicinity of the wells or near saturation fronts, the IMPES method can work with reasonable large time steps. From this observation, Thomas and Thurnau (1983) proposed the use of the Adaptive Implicit Method (AIM), whose concept was later improved by Forsyth Jr. and Sammon (1986). The goal of the Adaptive Implicit Method is to advance saturation and pressure implicitly in the regions where the CFL is larger than unity or where large variations in the saturations occur. In other regions of the reservoir only the pressure is considered implicitly. One of the difficulties of the AIM implementation is the switching criterion to change the procedure to FI and vice-versa.

Thomas and Thurnau (1983) used as switching criterion a specified saturation threshold from the previous iteration level. If an IMPES volume experiences a variation in saturation larger than the specified threshold, this volume is calculated as being FI in the next Newton's iteration. Forsyth Jr. and Sammon (1986) found that this criterion does not allow a new cell to become IMPES again once it was switched to FI, because a FI cell can have large throughput and the variation in saturation can be small. In this case, if this volume is switched to IMPES large instabilities may result. Russel (1989) proposed a criterion based on the CFL number to determine which volume may switch from FI to IMPES and vice-versa. He provides several reasons for working with the criterion based on the CFL and several ways to approximately calculate it. Nevertheless, it does not present any result for two and three phase flows. Young and Russel (1993) proposed a new relation for calculating the CFL number and suggest applying it in conjunction with the calculation of the maximum variation in the saturation. In this case, solutions using the IMPES method with the CFL larger than unity were obtained without oscillations. Fung et al. (1989) proposed a criterion based on the local analysis of the approximate equations. In each volume the eigenvalues are calculated and recommended that the spectral radius be smaller than unity for treating the cell with the IMPES method. Otherwise, the volume must be considered by the FI procedure. This scheme was applied to three-phase flow in Cartesian grids only.

The objective of this work is to present a comparison among the above discussed methodologies employing unstructured Voronoi meshes, broading the range of application of adaptive implicit methods. The results are presented in terms of CPU time, number of Newton's iterations, number of time iterations, and the degree of implicitness. The stability criterion proposed by Fung et al. (1989) added to a specified variation in the saturation field was adapted to deal with arbitrary number of connections between control volumes. The black-oil model is employed and two-phase (oil-water) in two-dimensional case studies are investigated.

Nomenclature

A = reservoir area, m ² AIM = Adaptive Implicit Method b_ij = width of the face located between the gridpoints i and j, m c = compressibity of the fluid, Pa ^-1 B = volumetric formation factor d_ij = distance between the grid points i and j, m MBE = mass balance error FI = Fully Implicit method h = height of the reservoir, m IMPES = Implicit Pressure Explicit Saturation method J = Jacobian matrix K = absolute permeability, m ² k_r = relative permeability M = viscosity ratio N _ν = number of neighbors of control volume i NIN = total Newtonian iterations NIT = number of time steps P = pressure, Pa P_c = capillarity pressure, Pa PI = percentage of implicit control volumes P_i = initial pressure, Pa PVI = porous volume injected = volumetric flow rate at storage conditions per unit of reservoir volume, (m ³ /s)/m ³ q_pi = injected volumetric flow rate, m ³ /s q_pp = produced volumetric flow rate, m ³ /s r = vector of residues r_rw = well radius, m R = residue of Newton's equation R ^* = R multiplied by the inverse of the diagonal matrix block S_or = residual oil saturation S = saturation t = time, s T = transformation matrix TCPU = total CPU time T_ij = transmissibility factor V_i = Volume of the control volume, m ³ V_f,p = final volume of phase p, m ³ V_i,p = initial volume of phase p, m ³ x',y' = local coordinate system X = represent P or S_w

Greek Symbols

= maximum pressure variation, Pa = maximum saturation variation ΔS_w = threshold change in the water saturation ΔX = increment or decrement in X Δt = timestep, s α = eigenvalues = porosity λ = mobility λ₁ = bound for forward switching λ₂ = bound for backward switching µ = viscosity, Pa.s ε = user tolerance εⁿ = error in the time level n

εⁿ⁺¹

= error in the time level n+1 ρ = spectral radius

Subscripts

w denotes water phase or well o denotes oil phase p denotes phase p

Mathematical Model

Assuming there are only two immiscible phases in the reservoir, oil (o) and water (w), and neglecting the gravitational and the capillarity pressure effects, the volumetric conservation equation for the phase p can be written as

where is the porosity and Bp is the volumetric formation factor, P is the pressure in the reservoir, _p is the volumetric flow rate at stock tank conditions per unit volume of the reservoir. The mobility, λ_p, is defined as

where K is the absolute permeability, k_rp is the relative permeability and µ_p is the viscosity of phase p.

In Equation (2), the volumetric formation factor is included in the mobility phase because the volumetric balance is evaluated at stock tank conditions. Writing Eq. (1) for the oil and water phases there are three unknowns (S_o, S_wand P) and only two equations. The closure equation comes from the global mass conservation, given by

The velocity field can be found from the Darcy's equation. Details can be found in Aziz and Setari (1979).

Integration of the Governing Equations

Figure 1 shows a control volume of Voronoi type, where i indicates the gridpoint and the j's are its neighbors. For each gridpoint j, it is possible to align a local Cartesian system x'-y' in so that the x' axis (the line that joins gridpoints i and j) be perpendicular to the face of the control volume, allowing the evaluation of the fluxes using only two gridpoints, a property of local orthogonal grid systems (Palagi and Aziz, 1994).

Integrating Eq. (1) in space and in time, one obtains:

The superscript θ in Eq. (4) identifies the type of the procedure: θ = 1, FI; θ = 0, IMPES; and if θ assumes zero or one in different regions of the domain, the procedure is AIM. Nν is the number of neighbors of control volume i, and T_ijis the transmissibility, given by

where b and d are shown in Fig. 1 and h is the thickness of the control volume.

FI, IMPES and AIM Formulations

FI Formulation

In this formulation, the unknowns pressure P and the water saturation S_ware obtained simultaneously. For the mobilities, , k_rp,ij is computed using UDS (upstream differencing

scheme) and (µ_pB_p)_ijusing an arithmetic averaging of the corresponding values at the gridpoints ij. The system of equations is solved iteratively using Newton's method. Making θ = 1 in Eq. (4), for all control volumes, the residual form of the equation reads

The Taylor series expansion of the residue gives

where ν is the iteration level and X represents the unknowns (P and Sw). When convergence is achieved the residue in the iteration ν+1 should be zero. Therefore,

The unknowns (P and S_w) are calculated after each Newtonian iteration, as

and the procedure is stopped when the following tolerances are satisfied:

where and are input parameters.

IMPES Formulation

In this procedure, the terms which contain the fluid, mobility and the production or injection flow rates are based on the saturation and pressure from the previous time level. Making θ = 0 in Eq. (4), a similar equation to Eq. (6) is obtained. All terms depending on the saturation are explicitly calculated, except the term which involves the saturation at the level (n+1)-th, which is, in fact, the required transient term. Therefore, when Newton's linearization is applied in each Jacobian block line, the only full block is the diagonal one, as shown in Fig. 2a. To obtain the pressure equation, each Jacobian block line is multiplied by its inverse diagonal block resulting in the block line of Fig. 2b. Also, the inverse of each diagonal block is multiplied by the original vector , resulting in the vector. As shown in Fig. 2c, the pressure equation is decoupled from the saturation equation. When this procedure is applied to each control volume, the pressure equation is obtained using the first line of each block of Fig. 2c, yielding the system

After solving Eq. (11) for Δ^ν+1 , P is obtained by

and S_w by

where J_ij are the coefficients of the second line of control volume i, according to Fig. 2c.

AIM Formulation

The IMPES scheme requires considerably less computational effort per timestep than the FI. However, it is often limited to small timesteps due to instabilities inherent to explicit schemes. This occurs in regions where the velocities or gradients of saturation are high. On the other hand, in the FI scheme the size of the timestep is bounded only by the accuracy of the solution. The goal of the AIM method is to put together the advantages of each methodology, using the FI method in the regions where the IMPES method shows instabilities and IMPES in the rest of the reservoir. Recall that pressure is evaluated implicitly all over the reservoir.

Evaluating the production/injection terms implicitly and considering that the mobilities of the phases can be considered implicit (FI) or explicit (IMPES), a similar equation to Eq. (6) is obtained. Using the same procedure used for the pressure equation in the IMPES method, a similar equation to Eq. (11) is obtained, where the only difference is that the entries of the Jacobian matrix can be blocks (1x2), (2x1), (2x2) or (1x1).

Switching from an IMPES to a FI Volume and Vice-versa

The success of an AIM scheme depends strongly on the efficiency of the switching criterion. A sensitive and practical criterion is not easily found. Here a local stability criterion of the Jacobian matrix akin to the one proposed by Fung et al. (1989) is employed. Denoting the residue of the conservation equations by R and the variables (P and S_w) by X, Newton's iteration can be written, in the index notation, as

where N is the total number of equations. In Eq. (14), is the solution in the (ν+1)-th iteration, and is the residual at the ν-th iteration of the i-th equation. converges when . The term is written as

where m refers to the last iteration at the time level n. If the error at the time levels n and n+1 are denoted by εⁿ and εⁿ⁺¹ , respectively, Eq. (15) can be written as

where the signal of ε is neglected, since only its norm is of interest. Considering

Eq. (16) can be written as

where T is the transformation matrix. The condition of stability is

which guarantees that the error introduced in a specified time level is not amplified in the next one. This condition is satisfied if T is normal and

where αj's are the eigenvalues of T, and ρ(T) is the spectral radius of T.

These are the stability conditions of the whole system. Computing the spectral radius is obviously impracticable. One can however identify the control volumes that yield instability by computing easily the spectral radius of its respective small submatrix. Eq. (18) can be written for each control volume and the maximum eigenvalue of T calculated. If the Eq. (20) holds, the cell can be IMPES, otherwise it is FI (see Fung et al. 1989). Two important simplifications are convenient and produce good results. The accumulation term at time level n can be excluded from Ri. For the evaluation of the derivatives, the flow is considered incompressible, which simplifies the calculation of the eigenvalues. Tests without this assumption produced similar results.

In the original proposal of Fung et al. (1989), three eigenvalues are required for all control volumes. To avoid that, they suggested that only the eigenvalues at the boundary between implicit and explicit regions be calculated. This implies that when there is a switching in the tested control volume, all the explicit neighbors must be checked, if the switching is from explicit to implicit (forward switching). Otherwise, if the switching is from implicit to explicit (backward switching), all the implicit neighbors must be checked. This procedure should be repeated until no volumes should be switched. Recalling that the compressibility effects are small in oil-water flows, the incompressibility condition was considered for calculating the eigenvalues. This implies that only one eigenvalue need to be calculated for each control volume, allowing, due to the minor computational effort required, the calculation of the eigenvalues for all control volumes. It was not possible to reproduce the implicit front reported by Fung et al. (1989) using these stability criteria. However, considering the compressibility effects and calculating the two eigenvalues for all control volumes, the results were identical. A possible reason for that is because one is solving two-phase flows, while in the work of Fung et al. (1989) three-phase flows were treated.

Searching for a more robust switching criterion, the following results where obtained with the method just described in conjunction with a switching criterion based on the variation of the saturation. The errors in the mass balance of oil and water are also reported on the tables of results. The mass balance error of each component, MBEp, is calculated by

where V_i,p and V_f,p are the initial and final volumes of phase p in a specified time level, q_wi is the volumetric flow rate of injected water and ql is the volumetric flow rate of the producing liquid (oil+water). The numerator of Eq. (21) represents the error in the volume of the phase p, in each time level, and the denominator the volume of injected water, both integrated from the starting.

Results

All cases analyzed in this section were obtained with = 68,95 kPa and = 510^-3. Similar CPU time were observed between the AIM and FI methodology for tests with = 10 and = 68,95 kPa.

For the AIM and FI methodologies, Newton's iteration is performed until the variation in saturation and pressure satisfies the tolerance for every control volume, regardless whether the saturation was treated explicitly or implicitly. As expected, the number of Newton's iterations increased considerably per timestep in the IMPES method. For the sake of comparison among the methodologies, all of them used the Jacobian matrix structure. Therefore, no criterion based on pre-fixed values was used for the IMPES method. For all cases the GMRES solver (Saad and Schultz, 1986) with ILU(1) preconditioning was employed.

In all tables of results, the eigenvalues λ₁ is the bound for the forward switching, λ₂ is the bound for the backward switching, and ΔS_wis the threshold change in the water saturation used for the implicit switching. They are listed as (λ₁, λ₂, ΔS_w), and the missing value, denoted by (*) means that the respective criterion was not used. Using the threshold change in the saturation, whenever a cell is switched, all its neighbors are checked out. For the backward switching, some tests were done to find out how frequent it must be checked. The switching after 20 iterations in time provides good results. Switching after 10, 20, 30 and 50 iterations is also considered. For all cases the results are reported for a dimensionless time level (PVI) equal to 1.82, that is, the volume of liquid injected is 1.82 times the porous volume of the reservoir. Also NIN is the total Newtonian iterations, NIT and PI are, respectively, the number of timesteps and the percentage of implicit control volumes required to reach 1.82 PVI. The CPU time were normalized by the minimum CPU time among the methods. For evaluating the pressure and saturation, an integration process with a variable timestep is used as in Aziz and Settari (1979). Except for the IMPES method, the timestep used was up to 50 days. The first test case, CASE 1, considers a flow whose relative permeability curves are quadratic functions of the saturation. The hybrid grid is shown in Fig. 3, the physical and geometric data are listed in Tab. 1, and the relative permeability and viscosity ratio curves are given by

The results are presented in Tabs. 2 and 3 for M = 10 and 50, respectively. From Tab. 2, it can be seen that the criterion based on the stability analysis requires only a small part of the volumes to be implicit. There is no effect on the tested eigenvalue bounds for the backward and forward switching. A value of ΔS_w = 0.01 in conjunction with the stability criterion provide the best results. Using only the switching criterion based on the variation of the saturation, a large part of the volumes turn to FI at the end of the simulation, especially for small ΔS_w values. There is a situation where the AIM method required more CPU time than the FI approach. In this case, there is a large difference in the number of Newtonian iterations, meaning that the criterion used by the AIM method does not detect the regions of appearance of instabilities. Also, the errors of the IMPES mass balance are larger than the AIM and FI. This fact can be explained considering that only one Newtonian iteration is performed, and the same convergence criterion in the solver for the three schemes was employed.

For M = 10, Figure 4 presents the implicit fronts for four PVI values. In this case λ₁ = 1.0, λ₂ = 0.8, and ΔS_w= 0.01. It is verified that for higher PVI values only a small number of volumes were kept implicit. There is a considerable reduction in the number of implicit cells from 0.2 to 0.6 PVI. This occurred because the water front at 0.6 PVI had already reached the producing wells. After the water breakthrough, the number of implicit cells decreases proportionally to the increasing of the saturation in the wells. The water breakthrough occurs at 0.3 PVI approximately, according to Fig. 6.

Figures 5 and 6 show the oil recovery and the water-cut curves for the AIM, FI, and IMPES methodologies for the parallel wells. The water cut is defined as the ratio of produced water volume and the total volume of produced liquid (oil+water). A good agreement among the methods is observed. Figure 7 presents the water cut for the IMPES using a timestep of 50 days. Although the oscillations were not observed in the oil recovery curves with the IMPES method, the oscillations on the water cut were detected with timesteps ranged from 2 to 50 days. This fact suggests that the oil recovery cannot be the only parameter for comparison, because it represents an integral value. In this case, the oscillations are smoothed out by the integration procedure and cannot be detected.

In Table 3, for M = 50, the stability criterion did not capture the appearance of instability, making the AIM performance inferior compared to FI. Also, this flow is more prone to instabilities than the one in Tab. 2, because of the larger viscosity ratio.

The percentages of implicit volumes were larger than those in Tab. 2. It is verified that the criterion was not able to detect the appearance of stability to keep robustness in Newton's method and, therefore, an excessive number of Newton iterations are required. Even dealing with the instability criteria in conjunction with the variation in the saturation, it is not detected instabilities for ΔS_w = 0.02 during some period of the simulation, resulting an excessive number of Newtonian iterations when compared to the same case with ΔS_w = 0.01.

A second test problem, again using the five-spot configuration, with the mesh presented in Fig. 8 is solved. In spite of being a hexagonal mesh, it resembles the Cartesian mesh used by Fung et al. (1989) in developing the stability criterion. The reservoir data are listed in Tab. 4, the relative permeabilities in Tab. 5 and viscosity curves given by

Table 6 shows that a small amount of implicit control volumes are required, even when only the threshold switching is used. Except for the IMPES method, the number of timesteps for all the simulations is coincident. The AIM scheme with the switching criterion based on the threshold change for ΔS_w= 0.02 provide the best results.

The last test problem, Case 3, involves a reservoir with eight wells using again a hybrid mesh as shown in Fig. 9.

The physical and geometrical data are listed in Tab. 7, the viscosity curves are given by Eq. (23), and the relative permeabilities by

The results are presented in Tab. 8. By decreasing λ₁ and λ₂, the local switching criterion capture the region where the implicit treatment is needed. Again, the local switching linked to the threshold criterion ΔS_w provides good results.

Conclusions

This work presented a comparison among three widely used methodologies employed in petroleum reservoir simulation. The tests were realized for two phase flows (oil+water) for areal reservoirs. The equations were discretized in a control volume framework in conjunction with unstructured Voronoi meshes. For switching a control volume from IMPES to FI and vice-versa, a stability criterion based on the eigenvalues, added to a threshold variation in water saturation was used. From the analysis, it was concluded that the AIM method can considerably save CPU time and memory if an efficient switching criterion is used. It was also observed that the stability criterion based on the eigenvalues alone was not able to detect regions where the saturation should be treated implicitly. A possible reason for this is the application of the criterion, in this work, to two-phase flow (oil+water) problems, in which the physical instabilities are much smaller than those of three phase flows (oil, gas, and water). It became clear, for most of the cases, that a stability criterion added to a water saturation threshold gives the best results in terms of CPU time, percentage of implicit volumes, and Newtonian iterations. Due to time step limitations, the IMPES method shows the worst performance when compared to FI and AIM methods. In spite of the better performance of AIM when compared to IMPES and FI, one believes that more effort need to be devoted in finding an efficient criterion for switching a control volume from IMPES to FI.

Acknowledgements

This work was partially supported by CNPq-Brazilian Council of Research through grants to the authors.

Paper accepted July, 2009.

Technical Editor: Celso K. Morooka

References

• Aziz, K., Settari, A., 1979, "Petroleum Reservoir Simulation", Applied Science Publishers, London.
• Forsyth Jr., P.A., Sammon, P.H., 1986, "Practical considerations for adaptive implicit methods in reservoir simulation", Journal of Computational Physics, Vol. 62, pp. 265-281.
• Fung, L.S.K., Collins, D.A., Nghiem, L.X., 1989, "An adaptive-implicit switching criterion based on numerical stability analysis", SPERE, February, pp. 45-51.
• Palagi, C.L., Aziz, K., 1994, "Use of Voronoi grid in reservoir simulation", SPE Advanced Technology Series, Vol. 2, No. 2, pp. 69-77.
• Peaceman, D.W., 1977, "Fundamentals of Numerical Reservoir Simulation", Elsevier Scientific Publishing Company, Amsterdam.
• Russell, T.F., 1989, "Stability analysis and switching criteria for adaptive implicit methods based on the CFL condition", paper SPE 18416, SPE Symposium on Reservoir Simulation, Houston, Texas.
• Saad, Y., Schultz, M.H., 1986, "GMRES: A general minimal residual algorithm for solving nonsymmetric linear systems", SIAM Journal Sci. Stat. Comput., Vol. 7, pp. 856-869.
• Thomas, G.W., Thurnau, D.H., 1983, "Reservoir simulating using an adaptive implicit method", SPEJ, Vol. 23, No. 5, pp. 759-768.
• Young, L.C., Russel, T.F., 1993, "Implementation of an adaptive implicit method", paper SPE 25245, SPE Symposium on Reservoir Simulation, New Orleans, Louisiana.

Publication Dates

• Publication in this collection
  19 Feb 2010
• Date of issue
  Dec 2009

History

• Accepted
  July 2009
• Received
  July 2009