- Email: [email protected]
A comparison of the finite-difference and the finite-element methods for simulating unstable displacements
Journal of Petroleum Science and Engineering, 5 ( 1991 ) 205-218
205
Elsevier Science Publishers B.V., Amsterdam
A comparison of the finite-difference and the finite-element methods for simulating unstable displacements Santanu Khataniar and Ekwere J. Peters Petroleum Engineering Dep., The University of Texas at Austin, Austin, TX 78712, U.S.A. (Received May 1, 1990; accepted after revision June 10, 1990)
ABSTRACT Khataniar, S. and Peters, E.J., 1991. A comparison of the finite-difference and the finite-elementmethods for simulating unstable displacements. J. Pet. Sci. Eng., 5: 205-218. The finite-difference and the finite-element methods are the two most commonly used numerical methods in reservoir simulation. A comparative study of the two methods in compositional and noncompositionai simulations is presented. The comparisons are based on grid orientation effect and computational efficiency. In the grid orientation study, the finite-elementmodel with Gaussian quadrature showed less grid orientation than either the finite-difference model or the finite-elementmodel with Lobatto quadrature. In heterogeneous media, the finite-elementmodel with Gaussian quadrature suffered from severe numerical oscillations inspite of upstream weighting. The finite-element model with Lobatto quadrature was stable in heterogeneous media; however, the results were essentially identical to those obtained by the finite-difference method. In simulating two-phase noncompositionaldisplacements, the finite-elementmodels took twice as much computing time as the finite-difference model due to the lengthy process of formulating and assembling the transmissibilitymatrices. In compositional systems, however, the finite-element method took only 50% more time than the finite-difference method because the bulk of the computational effort was concentrated in the phase behavior calculations rather than in the formulation of the system matrices.
Introduction
Traditionally, the five-point finite-difference method is used in reservoir simulation. Although the method is adequate for modelling fluid displacements in general, it suffers from certain drawbacks such as numerical dispersion and grid orientation, especially in the simulation of unstable processes. As a result, a number of papers (Langsgrnd, 1976; Spivak et al., 1976; Huyakorn and Pinder, 1978; Dalen, 1979; Douglas et al., 1979; Kukreti et al., 1988; Khataniar and Peters, 1990) have appeared in which the finite-element method and its modifications have been recommended as an alternative modelling tool for such processes. The finite-element method also is not entirely free from numerical difficulties when applied 0920-4105/91/$03.50
to multiphase flow in porous media. For example, some researchers (Langsgrud, 1976; Spivak et al., 1976; Pinder and Gray, 1977) have shown that the finite-element method may converge to incorrect solutions or suffer from numerical oscillations under certain conditions for waterflooding problems. Several numerical procedures such as upstream weighting and capacitance lumping for the finite-element method have been proposed to remedy such problems. These procedures, while stabilizing the solutions, also increase numerical dispersion in the finite-element method (Cohen, 1985). The finite-element method is also known to be computationally more expensive than the finite-difference method (Douglas et al., 1969; McMichael and Thomas, 1973). However, there has not been
© 1991 - - Elsevier Science Publishers B.V.
S. KHATANIARAND E.J. PETERS
206
Nomenclature A
cross sectional area normal to flow (m 2) fugacity of the ith hydrocarbon component in gas phase (Pa) j~,, fugacity of the ith hydrocarbon component in oil phase (Pa) n total number of nodes ~\. number of hydrocarbon components P pressure (Pa) qw molar flow rate of the aqueous component (kgmole
fg,
S -'l )
q,
molar flow rate of the ith hydrocarbon component (kgmole s -t ) qf overall hydrocarbon molar flow rate (kgmole s- L) Sg gas phase saturation (dimensionless) So oil phase saturation (dimensionless) S~ water phase saturation (dimensionless) Tg molar transmissibility of gas phase (kgmole-s kg- ~) I~, molar transmissibility of oil phase (kgmole-s kg- ~) xi mole fraction of the ith hydrocarbon component in oil phase Yi mole fraction of the ith hydrocarbon component in gas phase _-, overall mole fraction of the ith component V divergence operator 0£2 boundary of the solution domain /) outward normal 0 fractional porosity, dimensionless p~ molar density of gas phase (kgmole m - 3) po molar density ofoil phase (kgmole m -3) p~ molar density of water phase (kgmole m - 3) £2 solution domain
a detailed comparison of the finite-difference and the finite-element methods that would enable a model developer to select one method over the other. The objective of our study was to perform such a detailed comparison of the two methods for the specific purpose of simulating viscous fingering. We used the five-point finite-difference method and the Galerkin finite-element method to simulate unstable displacements in noncompositional and compositional systems in homogeneous and heterogeneous media. In order to ensure a fair comparison, the two models were developed and coded by the same individual using identical input data and solution technique. The only difference was in the method of discretization of the differential equations of the displacement to obtain the numerical model. In one case, a five-point fi-
nite-difference technique was used and in the other, a Galerkin-based finite-element technique was used. Upstream weighting was used in both methods. The finite-element method was further subdivided into two classes depending on the type of quadrature used to evaluate the integrals in the transmissibility terms: (a) the finite-element Gauss method, which used the Gauss quadrature, and (b) the finite-element Lobatto method, which used the Lobatto quadrature. The comparisons are presented in terms of grid orientation effect and computational efficiency. Qualitative comparisons of the solutions obtained by the three methods are also presented. Mathematical model
A generalized compositional model developed by Young and Stephenson (1983) was used in this study. The model for noncompositional oil-water systems was derived as a special case of the above model by suppressing the compositional part. In the following sections, the basic equations of the model and their finite-element formulations are described briefly. A more detailed description of the model is available in Young and Stephenson (1983). The compositional model consists of a set of molar material balance equations coupled with a set of phase equilibrium relations. The material balance equations are as follows: For the aqueous component: 0W
O~--V.(TwVP)=O
(1)
For the N o - 1 hydrocarbon components:
0
O(Fz,) Ot V.(Tox,+rgyi)VP=O
(2)
for i = 1,2,...,Nc- l Overall hydrocarbon balance:
OF O--~-V.(To + Tg)VP=O
(3)
207
FINITE-DIFFERENCE AND FINITE-ELEMENT METHODS FOR SIMULATING UNSTABLE DISPLACEMENTS
The terms W and F in the above equations represent the phase quantities of the aqueous component and overall hydrocarbon, respectively, and are defined as:
now be written in residual form using Eqs. 9 and 10:
R=¢F
foi =fgi, i = 1,2,...,Nc
(4)
A set of boundary conditions are used to introduce the source/sink terms into the material balances. The boundary condition at a no flow boundary is given by: VY=0
n
- E PjV.(TwVOj)
Application of Galerkin's criterion yields:
faRO~dl2=O for i= 1,2 ..... n i.e.:
(5)
The boundary conditions at a well (source/ sink) are given by: For the aqueous component:
-ATw~.Ve=qw
dW
j=l k dtJj j=l
W=pwSw, F=poSo +pgSg The phase equilibria relations are given by:
n
(6)
For the Arc- 1 hydrocarbon components:
--A ( ToXi + Tgyi)n" 17p=qi, i = 1,2,...,N¢- 1 (7)
The second term on the left-hand side of Eq. 11 is integrated by parts using the divergence theorem. Equation 11 then yields:
~ ~ [~d- ~W~) J fa OjOid£2+j= ,~PJfaa TwVOj. VOidf2 =
Overall hydrocarbon:
- A ( To + Tg)~'Ve=qf
(8)
Because all the material balance equations in the model have the same general form, only the aqueous component balance (Eq. 1 ) will be used to illustrate our formulation. The remaining material balances are treated similarly.
Galerkin formulation The primary variables in Eq. 1, P and W, are approximated by means of shape functions as follows: n
P= ~ Pj(t)Oj(x,y)
(9)
j=l
W= ~ Wj(t)Oj(x,y)
(10)
= fog2OiTwnoveda
(12)
the boundary integral on the fight-hand side of Eq. 12 is replaced by either qw or 0 to represent a source/sink or no-flow boundary according to the boundary conditions given by Eqs. 5 and 6. Equation 12 is a system of ordinary differential equations ( O D E ) w i t h P and Was the unknowns. The integrals on the left-hand side of this equation are computed numerically. The second integral on the left-hand side of Eq. 12 involves phase transmissibilities and will be referred to as the transmissibility integral. Evaluation of this integral is discussed in the following section.
Upstream weighting scheme
j=l
where O(x,y) is the bilinear shape function which is independent of time. Equation 1 can
In the finite-element method, the transmissibility integral on the left-hand side of Eq. 12
208
S. KHATANIAR AND E.J. PETERS
is evaluated numerically by using quadrature rules. In doing so, it is customary to apply an upstream weighting technique. In this study, the upstream weighting scheme proposed by Dalen (1979) is used because of its conceptual simplicity. According to this scheme, the transmissibility integral is approximated as shown below:
faTwVOj'VOidK2-~Twi~ VOj*VOid~
if I
rwjf VOj"VOidg2 if e, < ej (13) There are several advantages to using an upstream weighting scheme. Since transmissibility is a time dependent parameter, time consuming numerical integrations must be performed at each time step if the transmissibility integral shown on the left-hand side of Eq. 13 were to be evaluated directly. This would increase the computational burden of the method considerably. However, in the upstream weighted approximation shown on the right-hand side of Eq. 13, the time-dependent transmissibility term is no longer a part of the integrand. Since the shape functions are not time dependent, this integral needs to be evaluated only once for an entire simulation. Also, the transmissibility is evaluated at a node rather than at a quadrature point. This eliminates the need for interpolating nodal transmissibilities which would have been necessary had the upstream weighting not been applied. The upstream weighted finite-element method can be further subdivided into two categories depending on the quadrature rule used to evaluate the integral on the right-hand side of Eq. 13: (a) the finite-element Gauss method, which uses the Gauss quadrature, and (b) the finite-element Lobatto method, which uses the Lobatto quadrature. The essential difference between the two quadratures lies in the integral involving the diagonally opposite nodes of an element. In the Gauss quadrature,
this integral is generally nonzero which signifies fluid flow along the diagonals of an element. However, in the Lobatto quadrature, this integral is always zero which signifies no fluid flow along the diagonals of an element. Hence, the resultant computational molecule is similar to that of the finite-difference method (Young, 1981 ). After performing the integrations over all the elements and assembling the element matrices into global matrices, the following system of nonlinear ODE's is obtained:
B~+TP=R
(14)
where B and T are global coefficient matrices and R is a known vector incorporating the boundary, source and sink terms. In forming the coefficient matrix B, which arises from the first integral on the left-hand side of Eq. 12, the off diagonal terms are lumped onto the main diagonal so that matrix B is a diagonal matrix. This is known as the "lumped formulation" (Langsgrud, ! 976 ) which enhances the numerical stability of the solution. The time derivative of W in Eq. 14 is discretized using backward (implicit) difference in time. The remaining material balance equations (Eqs. 2 and 3 ) also give rise to ODE's similar to Eq. 14. After the time derivatives are discretized, a set of nonlinear algebraic equations are obtained from these ODE's. These equations, together with the phase equilibria relations are written in the residual form and solved by using the Newton-Raphson method. The details of this procedure may be found in Young and Stephenson (1983). The finite-difference model was a standard five-point central difference model with one point upstream weighting.
Results and discussion The finite-difference, the finite-element Gauss and the finite-element Lobatto methods were compared in both noncompositional and
209
FINITE-DIFFERENCE AND FINITE-ELEMENT METHODS FOR SIMULATING UNSTABLE DISPLACEMENTS
compositional systems in one and two dimensions. In l-D, the saturation or concentration profiles were qualitatively compared in a homogeneous porous medium. In 2-D, the three methods were compared on the basis of grid orientation effect and computational efficiency. Some qualitative comparison also was made among the saturation or concentration maps produced by the methods in a heterogeneous porous medium. The noncompositional study was a waterflood at a constant oil-water viscosity ratio. In the compositional study, two processes were examined: (a) an immiscible gas drive, and (b) a multi-contact miscible condensing (rich) gas drive.
A one-dimensional waterflood in a homogeneous porous m e d i u m at an oil-water viscosity ratio of 10 was simulated on a 100 grid block ( 100 element) mesh using the finite-difference, the finite-element Gauss and the finite-element Lobatto methods. Figure 1 shows the water saturation profiles at 0.35 pore volu m e injection obtained by the three methods and the analytical (Buckley-Leverett) solution to the waterflood problem. All the three 1.0
FD 8 • u
FE*Gauss BL FE-Lobatto
• u
FD FE-Gauss FE-Lobalto
0.6.
0.4' ,= o :E
I
0.2-
L 0.0 0.0
015
Dimensionless
1.0
Distance
numerical methods produced identical solutions that agreed well with the analytical solution. A one-dimensional immiscible lean gas drive process in a homogeneous porous m e d i u m was simulated using the compositional versions of the three models. A 100 grid block (100 element) mesh was used. The initial oil-injection gas viscosity ratio was 10. Figure 2 shows the concentration profiles of methane, the main component of the injection gas, at 0.40 pore volume injection. As in the noncompositional case, the three solutions are identical.
2-D non-compositional (waterflood) study
0.6"
Grid orientation
0.4
0.2
0.0 1 0.0
:E
---e---
Fig. 2. Methane concentration profile at 0.40 PV for 1-D compositionalsystem.
One-dimensional displacements
0.8"
o.8~ 1.0
015 Dimensionless
......
1.0
Distance
Fig. 1. Water saturation profiles at 0.35 PV for 1-D waterflood.
Grid orientation effect in the finite-difference, the finite-element Gauss and the finiteelement Lobatto methods was studied for an oil-water viscosity ratio of 100. The adverse viscosity ratio was chosen because grid orientation is known to be more severe for adverse viscosity ratios than for favorable viscosity ratios. Two sets of grids were used in the study: (a) a coarse 7 × 7 diagonal and 10 × 10 parallel grid, and (b) a refined 17 × 17 diagonal and 2 4 X 2 4 parallel grid. The purpose of the re-
210
S. K H A T A N I A R AND E.J. PETERS
fined grid set was to study the relative performance of the three methods under grid refinement since grid refinement is known to reduce the grid orientation effect. Note that the grid sizes are given in terms of grid blocks or elements so that a 7 X 7 grid for the finite element method consists of 8 X 8 nodes. Figures 3-5 show the oil recovery profiles for the three methods in the case of the 7 X 7 diagonal and 10 X 10 parallel grids. In the ab-
0.4
0.3
P o
>
o
o n-
0.2
O > Q. 0.1 t
~"
10x10 parallel(Lobatto)
0.4 i
0.£
L
i
1 "o
PV Water Injected
0.3
P o
Fig. 5. Oil recovery profiles tbr finite- element Lobatto method, coarse grid.
>
o
o o n-
0.2
O
0.4
13. 0.1
410xl
0 parallel(Gauss) 10
0.3
P
o.0"
I
0
1 PV
Water
¢D 0 0 0.2
Injected
Fig. 3. Oil recovery profiles for finite-element Gauss method, coarse grid,
O
0.1 r
0,4
0.0
0,3i
24x24 parallel(Gauss)
~
I I
PV W a t e r
2
,
2
Injected
F i g . 6. Oil recovery profiles for finite-element Gauss method, refined grid.
> O o Q (E
o
0.2
0 0. 0.1 t
410x10
o.o ""
i
o
I
parallel (FD)
PV Water Injected
Fig. 4, Oil recovery profiles for finite-difference method, coarse grid.
sence o f grid orientation effect, the recovery profiles for the diagonal and parallel grids should be identical. The larger the separation between the two recovery profiles, the more severe is the grid orientation effect. Thus, the finite-element method with Gauss quadrature (Fig. 3 ) shows the least grid orientation among the three methods. The finite-difference method (Fig. 4) shows more grid orientation than the finite-element Gauss method. The fi-
FINITE-DIFFERENCEAND FINITE-ELEMENTMETHODSFOR SIMULATINGUNSTABLEDISPLACEMENTS
2 11
0.4
I REAO,.PUT I "0
Io
0.3
2q)
I GENERATE TRANSMISSIBILITY MATRICES I
oo
t FORM RHS VECTOR I
> 0
n-
®1
0.2
~
O
YES
> Q. 0.1
,'
24x24 parallel (FD) I COMPUTE TIME DEPENDENT JACOBIAN ELEMENTS
I
I
0.0
i
!
I SOLVE AND UPDATE VARIABLES I
=
/
1 PV Water Injected
I COMPUTE PHASE DENSITY AND SATURATION
®1
Fig. 7. Oil recovery profiles for finite-difference method, refined grid.
I OPOATET'MEt
0.4
"O
0.3
2a)
Fig. 9. Flow chart showing the sequence of major computational tasks.
>
o
o o n-
0.2 0.10
O
1> a. 0.1
= •
0.0
,
0
0.08
24x24 parallel(Lobatto)
I
1
PV Water Injected
.~_ 0.06 O Q '~
0.O4
D.. (~
0.O2
E
Fig. 8. Oil recovery profiles for finite-element Lobatto method, refined grid.
nite-element Lobatto method (Fig. 5) shows slightly more grid orientation than the finitedifference method. The superior behavior of the finite-element Gauss method is because it accounts for fluid flow along the diagonal of an element. This diagonal flow is neglected in the finite-difference and the finite-element Lobatto methods. With the 17 × 17 diagonal and 24 × 24 parallel grid system, the grid orientation effect was
0.00 matrix
iteration
total
Fig. 10. Computation times for 20 × 20 waterflood.
significantly reduced in all three methods due to the grid refinement. The recovery profiles are shown in Figs. 6-8. As in the previous case, the finite-element Gauss method showed the least, almost negligible, level of grid orienta-
212
S. KHATANIARAND E.J. PETERS o.~
Computational efficiency
0.4
E
o.~
G
0.2
¢~
0.1
E
0.0
matrix
iteration
total
Fig. 11. C o m p u t a t i o n times for 40 × 40 waterflood.
tion. The finite-difference and the finite-element Lobatto methods showed about the same amount of grid orientation in this case. Thus, the finite-element Gauss m e t h o d is the least susceptible of the three methods to the grid orientation effect. The performance of the finite-difference and the finite-element Lobatto methods is also satisfactory provided a fine mesh is used.
The performance of the three methods was compared on a 20 × 20 and a 40 X 40 grids in a heterogeneous medium. Linear displacements were simulated at an oil-water viscosity ratio of 100. To ensure a fair comparison, the number of nodes in the finite-element methods was kept equal to the number of grid blocks in the finite-difference method to ensure that the same number of equations were solved in each method. Thus, the finite-element counterpart of a 20 × 20 block finite-difference grid consistent of a 19 X 19 element (20 X 20 node ) grid. For a clearer understanding of the computational procedure, a flow chart showing the major computational tasks is presented in Fig. 9. The task categories are the same for both the compositional and noncompositional models. The computational tasks between points 1 and 4 in Fig. 9 represent one time step while those between points 2 and 3 represent one iteration. It may be noted that the transmissibility matrices are formed only once in a time step. Figure 10 shows the computation time requirements of the three methods for the 20 X 20 grid. The finite-element methods took more than twice as long as the finite-difference method in the overall computation time. It can
Fig. 12. Water saturation m a p at 0.10 PV injection for finite-difference m e t h o d with 40 x 40 grid.
FINITE-DIFFERENCEAND FINITE-ELEMENTMETHODSFOR SIMULATINGUNSTABLEDISPLACEMENTS
213
Fig. 13. Water saturation map at 0.10 PV injection for finite-element Lobatto method with 40 × 40 grid.
Fig. 14. Water saturation map at 0. l 0 PV injection for finite-element Gauss method with 40 × 40 grid.
be seen that about one half of the total computational effort in the finite-element methods is concentrated in the formation of the transmissibility matrices. By contrast, this matrix formation is about fourteen times faster in the finite-difference method. Since the same number of equations are being solved in all three methods, the iteration times are nearly equal. The finite-element Gauss m e t h o d takes slightly more iteration time than the other two because its transmissibility matrices have
larger bandwidths which arise from the diagonal flow terms. Similar results were obtained for the 40 X 40 grid as shown in Fig. 11. However, the computation times for the finite-element Gauss m e t h o d are not available because it broke down soon after 0.10 pore volume injection due to numerical instability. The finite-elem e n t Lobatto m e t h o d took slightly less than twice as long as the finite-difference m e t h o d in this case. As in the previous case, the forma-
214
~ v ¢D > O O n"
S. K H A T A N I A R A N D E.J. P E T E R S 80
80
60
80
0 0
4O
c 0= O ¢3
40
n" C ¢,}
20
0
----=---
20
2 0 x 2 0 parallel (FD) ~
i
0
0
1
0
2
= 0
PV Injected
2 0 x 2 0 parallel
(Lobatto)
i 1
2
PV Injected
Fig. 15. Decane recovery profiles for the finite-difference method.
Fig. 17. Decane recovery profiles for the finite-element Lobatto method. I.;"
80
I.(
so Q
> 0 0
0.(
E 40
0 n"
(U
0 Q
0.(
E
0.,~
O. 0
O..C
0
~
20
) ~
i
20=(20 parallel ( G a u s s )
]
i
1
0.0
matrix
iteration
total
PV Injected
Fig. 16. Decane recovery profiles for the finite-element Gauss method.
tion of the transmissibility matrices was fourteen times faster in the finite-difference method than in the finite-element method. In order to make a qualitative comparison of the simulation results obtained from the three models, the water saturation contours at 0.10 pore volume injection for the 40 X 40 grid are plotted in Figs. 12-14. From Figs. 12 and 13, it can be seen that the finite-difference and the finite-element Lobatto methods produce
Fig. 18. Computation times for rich gas drive.
nearly identical saturation maps. The use of Gauss quadrature in the finite-element method, however, results in severe oscillations (Fig. 14) in the saturation values and the eventual breakdown of the solution. Such oscillations were also observed in the case of the 20 × 20 grid, but not in the homogeneous cases (grid orientation study). Thus, the finite-element Gauss method always produced oscillatory solutions in heterogeneous systems. The only difference between this method and the
FINITE-DIFFERENCE AND FINITE-ELEMENT METHODS FOR SIMULATING UNSTABLE DISPLACEMENTS 0.8"
Compositional study []
A la, G)
0.6
FD FE-Lobatto
Grid orientation
0
E O
0.4
O
E
¢., el
215
0.2
o
0.0
matrix
iteration
total
Fig. 19. Computation times for lean gas drive.
finite-element Lobatto method is the quadrature rule used to evaluate the transmissibility integrals. In the Lobatto quadrature, the diagonal flow is omitted whereas in the Gauss quadrature, the diagonal flow is incorporated. Therefore, it can be inferred that the numerical instability of the finite element-Gauss method is related in some way to the incorporation of the diagonal flow in the model. Further study regarding this aspect is necessary.
For the grid orientation study in compositional systems, an immiscible gas injection process at an oil-injection gas viscosity ratio of 100 was simulated. The injection gas was primarily methane and the resident fluid was primarily decane. A 14 X 14 diagonal and 20 X 20 parallel grid systems were used in this study. Results of the grid orientation study are presented in Figs. 15-17. Figure 15 shows the decane recovery profiles for the finite-difference method using the diagonal and parallel grids. The large separation of the two recovery profiles indicates severe grid orientation. Figure 16 shows the recovery profiles for the finiteelement Gauss method. Although there still is substantial grid orientation, the effect is significantly less than that in the finite-difference method. Recovery profiles from the finite-element Lobatto method are shown in Fig. 17. This method shows the most grid orientation of the three methods.
Computational efficiency The performance of the three models was compared on a 20 X 20 blocks (20 X 20 nodes)
Fig. 20. Methane concentration map at 0.25 PV injection for finite-difference method, rich gas drive.
2 16
S. KHATANIAR AND E.J. PETERS
Fig. 21. Methane concentration map at 0.25 PV injection for finite-element Lobatto method, rich gas drive.
Fig. 22. Methane concentration map at 0.25 PV injection for finite-element Gauss method, rich gas drive. grid in a heterogeneous medium. Two linear displacement processes were simulated: (a) a multi-contact miscible rich gas drive, and (b) an immiscible lean gas drive. Oil-injection gas viscosity ratio was kept at 100 for both processes. Computational time requirements for the finite-difference and the finite-element Lobatto methods for the rich gas drive are shown in Fig. 18. Data for the finite-element Gauss m e t h o d are not available because it failed soon after
0.25 pore volume injection due to numerical instability. In this case, the finite-element Lobatto m e t h o d took about 20% more overall computation time than the finite-difference method. The improved performance o f the finite-element m e t h o d relative to the finite-difference m e t h o d is because the iteration time is now large compared to the transmissibility matrix formation time. The large iteration time in the compositional models results from the lengthy phase behavior calculations. The ma-
FINITE-DIFFERENCE A N D F I N I T E - E L E M E N T M E T H O D S FOR S I M U L A T I N G UNSTABLE DISPLACEMENTS
217
( \
o ,,,I
Fig. 23. Methane concentration map at 0.25 PV for immiscible gas drive using the finite-difference method.
Fig. 24. Methane concentration map at 0.25 PV for immiscible gas drive using the finite-element Lobatto method.
trix formation is still about nine times slower than that of the finite-difference method. Figure 19 shows the computation times for the immiscible displacement. Once again, the finite-element m e t h o d broke down rapidly when Gauss quadrature was used. The finiteelement Lobatto method took nearly 45% more overall computation time than the finite-difference method. This is because the iteration time was not as large as it was in the rich gas
drive so that the computational burden shifted from the iterations to formulating the matrices, a task in which the finite-element method is much slower than the finite-difference method. Figures 20-22 show the methane mole fraction contours for the rich gas displacement at 0.25 pore volume injection for the three methods. The finite-element Lobatto method yields somewhat longer fingers than the finite-differ-
2 18
ence method. Similar results were observed in the immiscible gas drive as well (Figs. 23 and 24). This small difference could be due to the way the heterogeneous permeability field is modelled in the two methods. In the finite-difference model, an entire grid block has the same permeability whereas in the finite-element model, the permeability varies from node to node in a piecewise linear manner. The concentration map shown in Fig. 22 for the finiteelement Gauss method suffers from numerical oscillations. Conclusions (1) The finite-difference method and the upstream weighted finite-element methods produce identical results for linear displacements in a homogeneous porous medium. (2) The use of Gauss quadrature in the finite-element method produces less grid orientation than Lobatto quadrature and the finitedifference method. ( 3 ) The finite-element method is computationally less efficient than the finite-difference method. The degree of this inefficiency depends on the relative computational effort spent in formulating the space discretization matrices (transmissibility matrices). (4) The finite-element method with Gauss quadrature is not suitable for simulating multiphase flow in heterogeneous porous media because of its inherent numerical instability. ( 5 ) The finite-difference method provides a more robust and efficient means of simulating multi-phase flow in heterogeneous systems than the finite-element method. References Cohen, M.F., 1985. Finite Element Methods for Enhanced Oil Recovery Simulation. Paper SPE13512,
S. KHATANIAR AND E.J. PETERS
1985SPE Reservoir Simulation Symposium, Dallas, TX, Feb. 10-13. Dalen, V., 1979. Simplified finite-element models for reservoir flow problems. Soc. Pet. Eng. J., (Oct.) 19 ( 5 ): 333-343; Trans. AIME, 267. Douglas, J., Jr., Dupont, T. and Rachtbrd, H.H., 1969. The application of variational methods to waterflooding problems, J. Can. Pet. Technol. (July-Sept.) 8 ( 3 ): 79-85. Douglas, J., Jr., Darlow, B.L., Wheeler, M. and Kendall, R.P., 1979. Self-Adaptive Galerkin Methods for OneDimensional, Two-Phase Immiscible flow. Paper SPE7679, 1979 SPE Fifth Symposium on Reservoir Simulation, Denver, CO, Feb. 1-2. Huyakorn, P.S. and Pinder, G.F., 1978. A new finite element technique for the solution of two-phase flow through porous media. Adv. Water Resour., 1 ( 5 ): 285298. Khataniar, S. and Peters, E.J, t990. A finite element method for simulating unstable two-phase flow. J. Pet. Sci. Eng., (4): 169-181. Kukreti, A.R., Zaman, M.M., Civan, F. and Rajapaksa, Y., 1988. Finite-element analyses of waterflooding in irregularly shaped reservoirs. J. Pet. Sci. Eng., 1 (4): 277-294. Langsgrud, O., 1976. Simulation of two-phase flow by finite element method. SPE Fourth Syrup. Numerical Simulation of Reservoir Performance, Los ANgeles, CA; Pap. SPE 5725. McMichael, C.L. and Thomas, G.W., 1973. Reservoir simulation by Galerkin's method. Soc. Pet. Eng. J. (June), 125-138. Pinder, G.F. and Gray, W.G., 1977. Finite Element Simulation in Surface and Subsurface Hydrology. Academic Press, New York. N.Y., pp. 161-169. Spivak, A., Price, H.S. and Settari, A., 1976. Solutions of the equations for multi-dimensional two-phase immiscible flow by variational methods. SPE Fourth Symp. Numerical Simulation of Reservoir Performance, Los Angeles, CA; Pap. SPE 5723. Young, LC., 1981. A finite-element method for reservoir simulation. Soc. Pet. Eng. J. (Feb.), 115-128 Young, L.C. and Stephenson, R.E., 1983. A generalized compositional approach for reservoir simulation. Soc. pet. Eng. J., (Oct.): 727-742.
205
Elsevier Science Publishers B.V., Amsterdam
A comparison of the finite-difference and the finite-element methods for simulating unstable displacements Santanu Khataniar and Ekwere J. Peters Petroleum Engineering Dep., The University of Texas at Austin, Austin, TX 78712, U.S.A. (Received May 1, 1990; accepted after revision June 10, 1990)
ABSTRACT Khataniar, S. and Peters, E.J., 1991. A comparison of the finite-difference and the finite-elementmethods for simulating unstable displacements. J. Pet. Sci. Eng., 5: 205-218. The finite-difference and the finite-element methods are the two most commonly used numerical methods in reservoir simulation. A comparative study of the two methods in compositional and noncompositionai simulations is presented. The comparisons are based on grid orientation effect and computational efficiency. In the grid orientation study, the finite-elementmodel with Gaussian quadrature showed less grid orientation than either the finite-difference model or the finite-elementmodel with Lobatto quadrature. In heterogeneous media, the finite-elementmodel with Gaussian quadrature suffered from severe numerical oscillations inspite of upstream weighting. The finite-element model with Lobatto quadrature was stable in heterogeneous media; however, the results were essentially identical to those obtained by the finite-difference method. In simulating two-phase noncompositionaldisplacements, the finite-elementmodels took twice as much computing time as the finite-difference model due to the lengthy process of formulating and assembling the transmissibilitymatrices. In compositional systems, however, the finite-element method took only 50% more time than the finite-difference method because the bulk of the computational effort was concentrated in the phase behavior calculations rather than in the formulation of the system matrices.
Introduction
Traditionally, the five-point finite-difference method is used in reservoir simulation. Although the method is adequate for modelling fluid displacements in general, it suffers from certain drawbacks such as numerical dispersion and grid orientation, especially in the simulation of unstable processes. As a result, a number of papers (Langsgrnd, 1976; Spivak et al., 1976; Huyakorn and Pinder, 1978; Dalen, 1979; Douglas et al., 1979; Kukreti et al., 1988; Khataniar and Peters, 1990) have appeared in which the finite-element method and its modifications have been recommended as an alternative modelling tool for such processes. The finite-element method also is not entirely free from numerical difficulties when applied 0920-4105/91/$03.50
to multiphase flow in porous media. For example, some researchers (Langsgrud, 1976; Spivak et al., 1976; Pinder and Gray, 1977) have shown that the finite-element method may converge to incorrect solutions or suffer from numerical oscillations under certain conditions for waterflooding problems. Several numerical procedures such as upstream weighting and capacitance lumping for the finite-element method have been proposed to remedy such problems. These procedures, while stabilizing the solutions, also increase numerical dispersion in the finite-element method (Cohen, 1985). The finite-element method is also known to be computationally more expensive than the finite-difference method (Douglas et al., 1969; McMichael and Thomas, 1973). However, there has not been
© 1991 - - Elsevier Science Publishers B.V.
S. KHATANIARAND E.J. PETERS
206
Nomenclature A
cross sectional area normal to flow (m 2) fugacity of the ith hydrocarbon component in gas phase (Pa) j~,, fugacity of the ith hydrocarbon component in oil phase (Pa) n total number of nodes ~\. number of hydrocarbon components P pressure (Pa) qw molar flow rate of the aqueous component (kgmole
fg,
S -'l )
q,
molar flow rate of the ith hydrocarbon component (kgmole s -t ) qf overall hydrocarbon molar flow rate (kgmole s- L) Sg gas phase saturation (dimensionless) So oil phase saturation (dimensionless) S~ water phase saturation (dimensionless) Tg molar transmissibility of gas phase (kgmole-s kg- ~) I~, molar transmissibility of oil phase (kgmole-s kg- ~) xi mole fraction of the ith hydrocarbon component in oil phase Yi mole fraction of the ith hydrocarbon component in gas phase _-, overall mole fraction of the ith component V divergence operator 0£2 boundary of the solution domain /) outward normal 0 fractional porosity, dimensionless p~ molar density of gas phase (kgmole m - 3) po molar density ofoil phase (kgmole m -3) p~ molar density of water phase (kgmole m - 3) £2 solution domain
a detailed comparison of the finite-difference and the finite-element methods that would enable a model developer to select one method over the other. The objective of our study was to perform such a detailed comparison of the two methods for the specific purpose of simulating viscous fingering. We used the five-point finite-difference method and the Galerkin finite-element method to simulate unstable displacements in noncompositional and compositional systems in homogeneous and heterogeneous media. In order to ensure a fair comparison, the two models were developed and coded by the same individual using identical input data and solution technique. The only difference was in the method of discretization of the differential equations of the displacement to obtain the numerical model. In one case, a five-point fi-
nite-difference technique was used and in the other, a Galerkin-based finite-element technique was used. Upstream weighting was used in both methods. The finite-element method was further subdivided into two classes depending on the type of quadrature used to evaluate the integrals in the transmissibility terms: (a) the finite-element Gauss method, which used the Gauss quadrature, and (b) the finite-element Lobatto method, which used the Lobatto quadrature. The comparisons are presented in terms of grid orientation effect and computational efficiency. Qualitative comparisons of the solutions obtained by the three methods are also presented. Mathematical model
A generalized compositional model developed by Young and Stephenson (1983) was used in this study. The model for noncompositional oil-water systems was derived as a special case of the above model by suppressing the compositional part. In the following sections, the basic equations of the model and their finite-element formulations are described briefly. A more detailed description of the model is available in Young and Stephenson (1983). The compositional model consists of a set of molar material balance equations coupled with a set of phase equilibrium relations. The material balance equations are as follows: For the aqueous component: 0W
O~--V.(TwVP)=O
(1)
For the N o - 1 hydrocarbon components:
0
O(Fz,) Ot V.(Tox,+rgyi)VP=O
(2)
for i = 1,2,...,Nc- l Overall hydrocarbon balance:
OF O--~-V.(To + Tg)VP=O
(3)
207
FINITE-DIFFERENCE AND FINITE-ELEMENT METHODS FOR SIMULATING UNSTABLE DISPLACEMENTS
The terms W and F in the above equations represent the phase quantities of the aqueous component and overall hydrocarbon, respectively, and are defined as:
now be written in residual form using Eqs. 9 and 10:
R=¢F
foi =fgi, i = 1,2,...,Nc
(4)
A set of boundary conditions are used to introduce the source/sink terms into the material balances. The boundary condition at a no flow boundary is given by: VY=0
n
- E PjV.(TwVOj)
Application of Galerkin's criterion yields:
faRO~dl2=O for i= 1,2 ..... n i.e.:
(5)
The boundary conditions at a well (source/ sink) are given by: For the aqueous component:
-ATw~.Ve=qw
dW
j=l k dtJj j=l
W=pwSw, F=poSo +pgSg The phase equilibria relations are given by:
n
(6)
For the Arc- 1 hydrocarbon components:
--A ( ToXi + Tgyi)n" 17p=qi, i = 1,2,...,N¢- 1 (7)
The second term on the left-hand side of Eq. 11 is integrated by parts using the divergence theorem. Equation 11 then yields:
~ ~ [~d- ~W~) J fa OjOid£2+j= ,~PJfaa TwVOj. VOidf2 =
Overall hydrocarbon:
- A ( To + Tg)~'Ve=qf
(8)
Because all the material balance equations in the model have the same general form, only the aqueous component balance (Eq. 1 ) will be used to illustrate our formulation. The remaining material balances are treated similarly.
Galerkin formulation The primary variables in Eq. 1, P and W, are approximated by means of shape functions as follows: n
P= ~ Pj(t)Oj(x,y)
(9)
j=l
W= ~ Wj(t)Oj(x,y)
(10)
= fog2OiTwnoveda
(12)
the boundary integral on the fight-hand side of Eq. 12 is replaced by either qw or 0 to represent a source/sink or no-flow boundary according to the boundary conditions given by Eqs. 5 and 6. Equation 12 is a system of ordinary differential equations ( O D E ) w i t h P and Was the unknowns. The integrals on the left-hand side of this equation are computed numerically. The second integral on the left-hand side of Eq. 12 involves phase transmissibilities and will be referred to as the transmissibility integral. Evaluation of this integral is discussed in the following section.
Upstream weighting scheme
j=l
where O(x,y) is the bilinear shape function which is independent of time. Equation 1 can
In the finite-element method, the transmissibility integral on the left-hand side of Eq. 12
208
S. KHATANIAR AND E.J. PETERS
is evaluated numerically by using quadrature rules. In doing so, it is customary to apply an upstream weighting technique. In this study, the upstream weighting scheme proposed by Dalen (1979) is used because of its conceptual simplicity. According to this scheme, the transmissibility integral is approximated as shown below:
faTwVOj'VOidK2-~Twi~ VOj*VOid~
if I
rwjf VOj"VOidg2 if e, < ej (13) There are several advantages to using an upstream weighting scheme. Since transmissibility is a time dependent parameter, time consuming numerical integrations must be performed at each time step if the transmissibility integral shown on the left-hand side of Eq. 13 were to be evaluated directly. This would increase the computational burden of the method considerably. However, in the upstream weighted approximation shown on the right-hand side of Eq. 13, the time-dependent transmissibility term is no longer a part of the integrand. Since the shape functions are not time dependent, this integral needs to be evaluated only once for an entire simulation. Also, the transmissibility is evaluated at a node rather than at a quadrature point. This eliminates the need for interpolating nodal transmissibilities which would have been necessary had the upstream weighting not been applied. The upstream weighted finite-element method can be further subdivided into two categories depending on the quadrature rule used to evaluate the integral on the right-hand side of Eq. 13: (a) the finite-element Gauss method, which uses the Gauss quadrature, and (b) the finite-element Lobatto method, which uses the Lobatto quadrature. The essential difference between the two quadratures lies in the integral involving the diagonally opposite nodes of an element. In the Gauss quadrature,
this integral is generally nonzero which signifies fluid flow along the diagonals of an element. However, in the Lobatto quadrature, this integral is always zero which signifies no fluid flow along the diagonals of an element. Hence, the resultant computational molecule is similar to that of the finite-difference method (Young, 1981 ). After performing the integrations over all the elements and assembling the element matrices into global matrices, the following system of nonlinear ODE's is obtained:
B~+TP=R
(14)
where B and T are global coefficient matrices and R is a known vector incorporating the boundary, source and sink terms. In forming the coefficient matrix B, which arises from the first integral on the left-hand side of Eq. 12, the off diagonal terms are lumped onto the main diagonal so that matrix B is a diagonal matrix. This is known as the "lumped formulation" (Langsgrud, ! 976 ) which enhances the numerical stability of the solution. The time derivative of W in Eq. 14 is discretized using backward (implicit) difference in time. The remaining material balance equations (Eqs. 2 and 3 ) also give rise to ODE's similar to Eq. 14. After the time derivatives are discretized, a set of nonlinear algebraic equations are obtained from these ODE's. These equations, together with the phase equilibria relations are written in the residual form and solved by using the Newton-Raphson method. The details of this procedure may be found in Young and Stephenson (1983). The finite-difference model was a standard five-point central difference model with one point upstream weighting.
Results and discussion The finite-difference, the finite-element Gauss and the finite-element Lobatto methods were compared in both noncompositional and
209
FINITE-DIFFERENCE AND FINITE-ELEMENT METHODS FOR SIMULATING UNSTABLE DISPLACEMENTS
compositional systems in one and two dimensions. In l-D, the saturation or concentration profiles were qualitatively compared in a homogeneous porous medium. In 2-D, the three methods were compared on the basis of grid orientation effect and computational efficiency. Some qualitative comparison also was made among the saturation or concentration maps produced by the methods in a heterogeneous porous medium. The noncompositional study was a waterflood at a constant oil-water viscosity ratio. In the compositional study, two processes were examined: (a) an immiscible gas drive, and (b) a multi-contact miscible condensing (rich) gas drive.
A one-dimensional waterflood in a homogeneous porous m e d i u m at an oil-water viscosity ratio of 10 was simulated on a 100 grid block ( 100 element) mesh using the finite-difference, the finite-element Gauss and the finite-element Lobatto methods. Figure 1 shows the water saturation profiles at 0.35 pore volu m e injection obtained by the three methods and the analytical (Buckley-Leverett) solution to the waterflood problem. All the three 1.0
FD 8 • u
FE*Gauss BL FE-Lobatto
• u
FD FE-Gauss FE-Lobalto
0.6.
0.4' ,= o :E
I
0.2-
L 0.0 0.0
015
Dimensionless
1.0
Distance
numerical methods produced identical solutions that agreed well with the analytical solution. A one-dimensional immiscible lean gas drive process in a homogeneous porous m e d i u m was simulated using the compositional versions of the three models. A 100 grid block (100 element) mesh was used. The initial oil-injection gas viscosity ratio was 10. Figure 2 shows the concentration profiles of methane, the main component of the injection gas, at 0.40 pore volume injection. As in the noncompositional case, the three solutions are identical.
2-D non-compositional (waterflood) study
0.6"
Grid orientation
0.4
0.2
0.0 1 0.0
:E
---e---
Fig. 2. Methane concentration profile at 0.40 PV for 1-D compositionalsystem.
One-dimensional displacements
0.8"
o.8~ 1.0
015 Dimensionless
......
1.0
Distance
Fig. 1. Water saturation profiles at 0.35 PV for 1-D waterflood.
Grid orientation effect in the finite-difference, the finite-element Gauss and the finiteelement Lobatto methods was studied for an oil-water viscosity ratio of 100. The adverse viscosity ratio was chosen because grid orientation is known to be more severe for adverse viscosity ratios than for favorable viscosity ratios. Two sets of grids were used in the study: (a) a coarse 7 × 7 diagonal and 10 × 10 parallel grid, and (b) a refined 17 × 17 diagonal and 2 4 X 2 4 parallel grid. The purpose of the re-
210
S. K H A T A N I A R AND E.J. PETERS
fined grid set was to study the relative performance of the three methods under grid refinement since grid refinement is known to reduce the grid orientation effect. Note that the grid sizes are given in terms of grid blocks or elements so that a 7 X 7 grid for the finite element method consists of 8 X 8 nodes. Figures 3-5 show the oil recovery profiles for the three methods in the case of the 7 X 7 diagonal and 10 X 10 parallel grids. In the ab-
0.4
0.3
P o
>
o
o n-
0.2
O > Q. 0.1 t
~"
10x10 parallel(Lobatto)
0.4 i
0.£
L
i
1 "o
PV Water Injected
0.3
P o
Fig. 5. Oil recovery profiles tbr finite- element Lobatto method, coarse grid.
>
o
o o n-
0.2
O
0.4
13. 0.1
410xl
0 parallel(Gauss) 10
0.3
P
o.0"
I
0
1 PV
Water
¢D 0 0 0.2
Injected
Fig. 3. Oil recovery profiles for finite-element Gauss method, coarse grid,
O
0.1 r
0,4
0.0
0,3i
24x24 parallel(Gauss)
~
I I
PV W a t e r
2
,
2
Injected
F i g . 6. Oil recovery profiles for finite-element Gauss method, refined grid.
> O o Q (E
o
0.2
0 0. 0.1 t
410x10
o.o ""
i
o
I
parallel (FD)
PV Water Injected
Fig. 4, Oil recovery profiles for finite-difference method, coarse grid.
sence o f grid orientation effect, the recovery profiles for the diagonal and parallel grids should be identical. The larger the separation between the two recovery profiles, the more severe is the grid orientation effect. Thus, the finite-element method with Gauss quadrature (Fig. 3 ) shows the least grid orientation among the three methods. The finite-difference method (Fig. 4) shows more grid orientation than the finite-element Gauss method. The fi-
FINITE-DIFFERENCEAND FINITE-ELEMENTMETHODSFOR SIMULATINGUNSTABLEDISPLACEMENTS
2 11
0.4
I REAO,.PUT I "0
Io
0.3
2q)
I GENERATE TRANSMISSIBILITY MATRICES I
oo
t FORM RHS VECTOR I
> 0
n-
®1
0.2
~
O
YES
> Q. 0.1
,'
24x24 parallel (FD) I COMPUTE TIME DEPENDENT JACOBIAN ELEMENTS
I
I
0.0
i
!
I SOLVE AND UPDATE VARIABLES I
=
/
1 PV Water Injected
I COMPUTE PHASE DENSITY AND SATURATION
®1
Fig. 7. Oil recovery profiles for finite-difference method, refined grid.
I OPOATET'MEt
0.4
"O
0.3
2a)
Fig. 9. Flow chart showing the sequence of major computational tasks.
>
o
o o n-
0.2 0.10
O
1> a. 0.1
= •
0.0
,
0
0.08
24x24 parallel(Lobatto)
I
1
PV Water Injected
.~_ 0.06 O Q '~
0.O4
D.. (~
0.O2
E
Fig. 8. Oil recovery profiles for finite-element Lobatto method, refined grid.
nite-element Lobatto method (Fig. 5) shows slightly more grid orientation than the finitedifference method. The superior behavior of the finite-element Gauss method is because it accounts for fluid flow along the diagonal of an element. This diagonal flow is neglected in the finite-difference and the finite-element Lobatto methods. With the 17 × 17 diagonal and 24 × 24 parallel grid system, the grid orientation effect was
0.00 matrix
iteration
total
Fig. 10. Computation times for 20 × 20 waterflood.
significantly reduced in all three methods due to the grid refinement. The recovery profiles are shown in Figs. 6-8. As in the previous case, the finite-element Gauss method showed the least, almost negligible, level of grid orienta-
212
S. KHATANIARAND E.J. PETERS o.~
Computational efficiency
0.4
E
o.~
G
0.2
¢~
0.1
E
0.0
matrix
iteration
total
Fig. 11. C o m p u t a t i o n times for 40 × 40 waterflood.
tion. The finite-difference and the finite-element Lobatto methods showed about the same amount of grid orientation in this case. Thus, the finite-element Gauss m e t h o d is the least susceptible of the three methods to the grid orientation effect. The performance of the finite-difference and the finite-element Lobatto methods is also satisfactory provided a fine mesh is used.
The performance of the three methods was compared on a 20 × 20 and a 40 X 40 grids in a heterogeneous medium. Linear displacements were simulated at an oil-water viscosity ratio of 100. To ensure a fair comparison, the number of nodes in the finite-element methods was kept equal to the number of grid blocks in the finite-difference method to ensure that the same number of equations were solved in each method. Thus, the finite-element counterpart of a 20 × 20 block finite-difference grid consistent of a 19 X 19 element (20 X 20 node ) grid. For a clearer understanding of the computational procedure, a flow chart showing the major computational tasks is presented in Fig. 9. The task categories are the same for both the compositional and noncompositional models. The computational tasks between points 1 and 4 in Fig. 9 represent one time step while those between points 2 and 3 represent one iteration. It may be noted that the transmissibility matrices are formed only once in a time step. Figure 10 shows the computation time requirements of the three methods for the 20 X 20 grid. The finite-element methods took more than twice as long as the finite-difference method in the overall computation time. It can
Fig. 12. Water saturation m a p at 0.10 PV injection for finite-difference m e t h o d with 40 x 40 grid.
FINITE-DIFFERENCEAND FINITE-ELEMENTMETHODSFOR SIMULATINGUNSTABLEDISPLACEMENTS
213
Fig. 13. Water saturation map at 0.10 PV injection for finite-element Lobatto method with 40 × 40 grid.
Fig. 14. Water saturation map at 0. l 0 PV injection for finite-element Gauss method with 40 × 40 grid.
be seen that about one half of the total computational effort in the finite-element methods is concentrated in the formation of the transmissibility matrices. By contrast, this matrix formation is about fourteen times faster in the finite-difference method. Since the same number of equations are being solved in all three methods, the iteration times are nearly equal. The finite-element Gauss m e t h o d takes slightly more iteration time than the other two because its transmissibility matrices have
larger bandwidths which arise from the diagonal flow terms. Similar results were obtained for the 40 X 40 grid as shown in Fig. 11. However, the computation times for the finite-element Gauss m e t h o d are not available because it broke down soon after 0.10 pore volume injection due to numerical instability. The finite-elem e n t Lobatto m e t h o d took slightly less than twice as long as the finite-difference m e t h o d in this case. As in the previous case, the forma-
214
~ v ¢D > O O n"
S. K H A T A N I A R A N D E.J. P E T E R S 80
80
60
80
0 0
4O
c 0= O ¢3
40
n" C ¢,}
20
0
----=---
20
2 0 x 2 0 parallel (FD) ~
i
0
0
1
0
2
= 0
PV Injected
2 0 x 2 0 parallel
(Lobatto)
i 1
2
PV Injected
Fig. 15. Decane recovery profiles for the finite-difference method.
Fig. 17. Decane recovery profiles for the finite-element Lobatto method. I.;"
80
I.(
so Q
> 0 0
0.(
E 40
0 n"
(U
0 Q
0.(
E
0.,~
O. 0
O..C
0
~
20
) ~
i
20=(20 parallel ( G a u s s )
]
i
1
0.0
matrix
iteration
total
PV Injected
Fig. 16. Decane recovery profiles for the finite-element Gauss method.
tion of the transmissibility matrices was fourteen times faster in the finite-difference method than in the finite-element method. In order to make a qualitative comparison of the simulation results obtained from the three models, the water saturation contours at 0.10 pore volume injection for the 40 X 40 grid are plotted in Figs. 12-14. From Figs. 12 and 13, it can be seen that the finite-difference and the finite-element Lobatto methods produce
Fig. 18. Computation times for rich gas drive.
nearly identical saturation maps. The use of Gauss quadrature in the finite-element method, however, results in severe oscillations (Fig. 14) in the saturation values and the eventual breakdown of the solution. Such oscillations were also observed in the case of the 20 × 20 grid, but not in the homogeneous cases (grid orientation study). Thus, the finite-element Gauss method always produced oscillatory solutions in heterogeneous systems. The only difference between this method and the
FINITE-DIFFERENCE AND FINITE-ELEMENT METHODS FOR SIMULATING UNSTABLE DISPLACEMENTS 0.8"
Compositional study []
A la, G)
0.6
FD FE-Lobatto
Grid orientation
0
E O
0.4
O
E
¢., el
215
0.2
o
0.0
matrix
iteration
total
Fig. 19. Computation times for lean gas drive.
finite-element Lobatto method is the quadrature rule used to evaluate the transmissibility integrals. In the Lobatto quadrature, the diagonal flow is omitted whereas in the Gauss quadrature, the diagonal flow is incorporated. Therefore, it can be inferred that the numerical instability of the finite element-Gauss method is related in some way to the incorporation of the diagonal flow in the model. Further study regarding this aspect is necessary.
For the grid orientation study in compositional systems, an immiscible gas injection process at an oil-injection gas viscosity ratio of 100 was simulated. The injection gas was primarily methane and the resident fluid was primarily decane. A 14 X 14 diagonal and 20 X 20 parallel grid systems were used in this study. Results of the grid orientation study are presented in Figs. 15-17. Figure 15 shows the decane recovery profiles for the finite-difference method using the diagonal and parallel grids. The large separation of the two recovery profiles indicates severe grid orientation. Figure 16 shows the recovery profiles for the finiteelement Gauss method. Although there still is substantial grid orientation, the effect is significantly less than that in the finite-difference method. Recovery profiles from the finite-element Lobatto method are shown in Fig. 17. This method shows the most grid orientation of the three methods.
Computational efficiency The performance of the three models was compared on a 20 X 20 blocks (20 X 20 nodes)
Fig. 20. Methane concentration map at 0.25 PV injection for finite-difference method, rich gas drive.
2 16
S. KHATANIAR AND E.J. PETERS
Fig. 21. Methane concentration map at 0.25 PV injection for finite-element Lobatto method, rich gas drive.
Fig. 22. Methane concentration map at 0.25 PV injection for finite-element Gauss method, rich gas drive. grid in a heterogeneous medium. Two linear displacement processes were simulated: (a) a multi-contact miscible rich gas drive, and (b) an immiscible lean gas drive. Oil-injection gas viscosity ratio was kept at 100 for both processes. Computational time requirements for the finite-difference and the finite-element Lobatto methods for the rich gas drive are shown in Fig. 18. Data for the finite-element Gauss m e t h o d are not available because it failed soon after
0.25 pore volume injection due to numerical instability. In this case, the finite-element Lobatto m e t h o d took about 20% more overall computation time than the finite-difference method. The improved performance o f the finite-element m e t h o d relative to the finite-difference m e t h o d is because the iteration time is now large compared to the transmissibility matrix formation time. The large iteration time in the compositional models results from the lengthy phase behavior calculations. The ma-
FINITE-DIFFERENCE A N D F I N I T E - E L E M E N T M E T H O D S FOR S I M U L A T I N G UNSTABLE DISPLACEMENTS
217
( \
o ,,,I
Fig. 23. Methane concentration map at 0.25 PV for immiscible gas drive using the finite-difference method.
Fig. 24. Methane concentration map at 0.25 PV for immiscible gas drive using the finite-element Lobatto method.
trix formation is still about nine times slower than that of the finite-difference method. Figure 19 shows the computation times for the immiscible displacement. Once again, the finite-element m e t h o d broke down rapidly when Gauss quadrature was used. The finiteelement Lobatto method took nearly 45% more overall computation time than the finite-difference method. This is because the iteration time was not as large as it was in the rich gas
drive so that the computational burden shifted from the iterations to formulating the matrices, a task in which the finite-element method is much slower than the finite-difference method. Figures 20-22 show the methane mole fraction contours for the rich gas displacement at 0.25 pore volume injection for the three methods. The finite-element Lobatto method yields somewhat longer fingers than the finite-differ-
2 18
ence method. Similar results were observed in the immiscible gas drive as well (Figs. 23 and 24). This small difference could be due to the way the heterogeneous permeability field is modelled in the two methods. In the finite-difference model, an entire grid block has the same permeability whereas in the finite-element model, the permeability varies from node to node in a piecewise linear manner. The concentration map shown in Fig. 22 for the finiteelement Gauss method suffers from numerical oscillations. Conclusions (1) The finite-difference method and the upstream weighted finite-element methods produce identical results for linear displacements in a homogeneous porous medium. (2) The use of Gauss quadrature in the finite-element method produces less grid orientation than Lobatto quadrature and the finitedifference method. ( 3 ) The finite-element method is computationally less efficient than the finite-difference method. The degree of this inefficiency depends on the relative computational effort spent in formulating the space discretization matrices (transmissibility matrices). (4) The finite-element method with Gauss quadrature is not suitable for simulating multiphase flow in heterogeneous porous media because of its inherent numerical instability. ( 5 ) The finite-difference method provides a more robust and efficient means of simulating multi-phase flow in heterogeneous systems than the finite-element method. References Cohen, M.F., 1985. Finite Element Methods for Enhanced Oil Recovery Simulation. Paper SPE13512,
S. KHATANIAR AND E.J. PETERS
1985SPE Reservoir Simulation Symposium, Dallas, TX, Feb. 10-13. Dalen, V., 1979. Simplified finite-element models for reservoir flow problems. Soc. Pet. Eng. J., (Oct.) 19 ( 5 ): 333-343; Trans. AIME, 267. Douglas, J., Jr., Dupont, T. and Rachtbrd, H.H., 1969. The application of variational methods to waterflooding problems, J. Can. Pet. Technol. (July-Sept.) 8 ( 3 ): 79-85. Douglas, J., Jr., Darlow, B.L., Wheeler, M. and Kendall, R.P., 1979. Self-Adaptive Galerkin Methods for OneDimensional, Two-Phase Immiscible flow. Paper SPE7679, 1979 SPE Fifth Symposium on Reservoir Simulation, Denver, CO, Feb. 1-2. Huyakorn, P.S. and Pinder, G.F., 1978. A new finite element technique for the solution of two-phase flow through porous media. Adv. Water Resour., 1 ( 5 ): 285298. Khataniar, S. and Peters, E.J, t990. A finite element method for simulating unstable two-phase flow. J. Pet. Sci. Eng., (4): 169-181. Kukreti, A.R., Zaman, M.M., Civan, F. and Rajapaksa, Y., 1988. Finite-element analyses of waterflooding in irregularly shaped reservoirs. J. Pet. Sci. Eng., 1 (4): 277-294. Langsgrud, O., 1976. Simulation of two-phase flow by finite element method. SPE Fourth Syrup. Numerical Simulation of Reservoir Performance, Los ANgeles, CA; Pap. SPE 5725. McMichael, C.L. and Thomas, G.W., 1973. Reservoir simulation by Galerkin's method. Soc. Pet. Eng. J. (June), 125-138. Pinder, G.F. and Gray, W.G., 1977. Finite Element Simulation in Surface and Subsurface Hydrology. Academic Press, New York. N.Y., pp. 161-169. Spivak, A., Price, H.S. and Settari, A., 1976. Solutions of the equations for multi-dimensional two-phase immiscible flow by variational methods. SPE Fourth Symp. Numerical Simulation of Reservoir Performance, Los Angeles, CA; Pap. SPE 5723. Young, LC., 1981. A finite-element method for reservoir simulation. Soc. Pet. Eng. J. (Feb.), 115-128 Young, L.C. and Stephenson, R.E., 1983. A generalized compositional approach for reservoir simulation. Soc. pet. Eng. J., (Oct.): 727-742.