INTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL METHODS IN GEOMECHANICS. VOL. 3.355-366 (1979) FINITE DIFFERENCE TECHNIQUES FOR MODELLING GEOTHERMAL RESERVOIRS G. A. ZYVOLOSKI, M. J. O’SULLIVAN AND D. E. KROL Department of Theoretical and Applied Mechanics, University of Auckland, Auckland, New Zeoland SUMMARY Two finite difference algorithms suitable for long-time simulation of the exploitation of a two-phase geothermal reservoir are presented. One is based on the hopscotch method proposed by Saul’yev’ and analysed further by Gordon’ and Gourlay.’ The other is based on the well-known AD1 method. Both methods use a Newton-Raphson iterative technique in order to obtain accurate solutions of the non-linear difference equations involved. Rapid convergence of the iterative schemes occurs both for single-phase and two-phase reservoir problems. One- and two-dimensional model problems are presented. INTRODUCTION The recent interest in the use of geothermal energy has created a need for mathematical models of geothermal reservoirs. The most successful modelling of one- and two-phase geothermal flows previously reported has been carried out b Brownell et a1.: Garg et ul.,’ Pritchett et a1.,6 and by Mercer and Faust7 and Faust and MercerJ9Several other authors have considered either one-phase or two-phase flows but not flows involving the transition from compressed water to steam and water (see Toronyi and Farouq Ali” for example). Pritchett et aL6 use the density and internal energy of the two-phase fluid as dependent variables. This approach enables accurate representation of the mass and energy accumulation terms in the basic mass and energy conservation equations but these variables are not convenient for dealing with the flux terms. These authors use the AD1 method for simplifying the solution of the resulting non-linear difference equations. However, their use of density and internal energy means that the efficient and rapidly convergent procedures used in the present work cannot be implemented. Mercer and Faust7 and Faust and use the same dependent variables, namely pressure and mixture enthalpy, as in the present work. However, they do not use the AD1 splitting of the problem and therefore are forced to deal directly with large sparse matrices. The most troublesome problem in the numerical modelling of geothermal reservoirs is the transition from one-phase to two-phase conditions. In the course of this transition properties of water change dramatically. To obtain a numerical scheme which remains stable during a phase transition it is necessary to take account of the change of fluid properties by including a weighted average of flux and accumulation terms at the beginning and end of the time step (that is, to use implicit time differencing). Also because of the non-linear dependence of mass and energy on pressure and enthalpy it is necessary to iterate to achieve sufficiently accurate mass and energy balances. Because of this requirement for an implicit iterative procedure it is desirable to choose a method which is arithmetically efficient. For this reason finite difference methods were selected, after some experimentation with finite element methods, because they lead to sparse matrix equations of a particularly simple form. 0363-9061/79/0403-0355$01.00 Received 8 August 1978 @ 1979 by John Wiley & Sons, Ltd. Revised 15 December 1978 355
Transcript
INTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL METHODS IN GEOMECHANICS. VOL. 3.355-366 (1979)
FINITE DIFFERENCE TECHNIQUES FOR MODELLING GEOTHERMAL RESERVOIRS
G . A. ZYVOLOSKI, M. J. O’SULLIVAN AND D. E. KROL
Department of Theoretical and Applied Mechanics, University of Auckland, Auckland, New Zeoland
SUMMARY
Two finite difference algorithms suitable for long-time simulation of the exploitation of a two-phase geothermal reservoir are presented. One is based on the hopscotch method proposed by Saul’yev’ and analysed further by Gordon’ and Gourlay.’ The other is based on the well-known AD1 method. Both methods use a Newton-Raphson iterative technique in order to obtain accurate solutions of the non-linear difference equations involved. Rapid convergence of the iterative schemes occurs both for single-phase and two-phase reservoir problems. One- and two-dimensional model problems are presented.
INTRODUCTION
The recent interest in the use of geothermal energy has created a need for mathematical models of geothermal reservoirs. The most successful modelling of one- and two-phase geothermal flows previously reported has been carried out b Brownell et a1.: Garg et ul.,’ Pritchett et a1.,6 and by Mercer and Faust7 and Faust and MercerJ9 Several other authors have considered either one-phase or two-phase flows but not flows involving the transition from compressed water to steam and water (see Toronyi and Farouq Ali” for example). Pritchett et aL6 use the density and internal energy of the two-phase fluid as dependent variables. This approach enables accurate representation of the mass and energy accumulation terms in the basic mass and energy conservation equations but these variables are not convenient for dealing with the flux terms. These authors use the AD1 method for simplifying the solution of the resulting non-linear difference equations. However, their use of density and internal energy means that the efficient and rapidly convergent procedures used in the present work cannot be implemented. Mercer and Faust7 and Faust and use the same dependent variables, namely pressure and mixture enthalpy, as in the present work. However, they do not use the AD1 splitting of the problem and therefore are forced to deal directly with large sparse matrices.
The most troublesome problem in the numerical modelling of geothermal reservoirs is the transition from one-phase to two-phase conditions. In the course of this transition properties of water change dramatically. To obtain a numerical scheme which remains stable during a phase transition it is necessary to take account of the change of fluid properties by including a weighted average of flux and accumulation terms at the beginning and end of the time step (that is, to use implicit time differencing). Also because of the non-linear dependence of mass and energy on pressure and enthalpy it is necessary to iterate to achieve sufficiently accurate mass and energy balances.
Because of this requirement for an implicit iterative procedure it is desirable to choose a method which is arithmetically efficient. For this reason finite difference methods were selected, after some experimentation with finite element methods, because they lead to sparse matrix equations of a particularly simple form.
0363-9061/79/0403-0355$01.00 Received 8 August 1978 @ 1979 by John Wiley & Sons, Ltd. Revised 15 December 1978
355
356 G. A. ZYVOLOSKI, M. J. O’SULLIVAN AND D. E. KROL
The methods used in this paper have two advantages. Firstly, they involve the solution of small and sparse matrix equations. In the AD1 method tridiagonal matrix equations are solved and the hopscotch method is completely explicit. The second advantage of these methods is that they require a small number of function evaluations of the non-linear terms in the difference equations. Because they work either line by line in the case of the AD1 method or point by point in the case of the hopscotch method the number of iterations required and hence the number of function evaluations is matched to the local accuracy of the solution. The procedure used by Mercer and Faust’ and Faust and Mercer8” requires global iteration and therefore many function evaluations. This point is particularly important since typically the point by point convergence of the solution of the non-linear difference scheme is very fast except at a few points where the flow is changing from one- to two-phase.
MATHEMATICAL FORMULATION
Detailed derivations of the governing equations for two-phase geothermal reservoirs have been presented by several investigators (Mercer et al.,” Mercer and Faust’ and Brownell et al.4 for example) and therefore only a brief development will be presented here.
Conservation of mass is expressed by the equation
aAm -+V. f,+q,=O at
where the mass per unit volume A, is given by
and the mass flux f, is given by fm = pvVv + PIVI (3)
Here (b is the porosity of the matrix, S, and S, are saturations, p, and p1 are densities, V, and VI are flow rates per unit volume with the subscripts v and 1 indicating quantities for the vapour phase or liquid phase, respectively. Sources or sinks (that is, bores or reinjection wells) are represented by the term 9,.
Conservation of energy is expressed by the equation
aAe -+v . r ,+q, = 0 at (4)
where the energy per unit volume A, is given by
Ae = (1 - 4)prur + (Svpvuv + SIPIUI) ( 5 )
and the energy flux f, is given by
f, =p,h,V,+plhlVl-KVT (6) Here the subscript r refers to the rock matrix, ur, uv and u1 are specific internal energies, h, and hl are specific enthalpies, K is an effective thermal conductivity, T is the temperature and 4. is the energy contributed from sources and sinks. The energy equation used by Pritchett et a1.6 ignores the small contribution to the energy flux term arising from the power input to a typical control volume by the pV term and therefore they use internal energy instead of enthalpy in the flux term. Mercer and Faust’ and Faust and Mercersv9 make an approximation of a similar order by replacing internal energies in the energy content A, by enthalpies.
MODELLING GEOTHERMAL RESERVOIRS 357
To complete the governing equations it is assumed that Darcy’s Law applies to the movement of each phase:
k R , v,= - - ( V p - p , g ) PV
(7)
Here k is the permeability, R, and R1 are the relative permeabilities, P,, and pl are viscosities, p is the pressure and g represents the acceleration due to gravity. The relative permeabilities used here are a version of a form suggested by Carey'* and adopted by Faust and M e r ~ e r : ~
R, = (1 - Sf2) (1 -SF)* (9)
RI = SF4 (10)
where S: = (Sl - Sh- Svr)/(l - S1, - Svr) and Slr, S,, are the residual saturations at which the liquid and vapour phases respectively become immobile. In the present work values of SIr = 0.05 and S, = 0.05 were used to obtain comparison with the work of Faust and Mercer.’ A value of Sh = 0.3 has been used by other authors and is probably closer to the true field value.
Using Darcy’s Law the basic conservation equations (1) and (4) can be rewritten in the form
and aA at
V . ( D e V p ) + V . ( k V T ) - q e - 2 = 0
Here the transmissibilities D, and D, are given by
and
with
PI PY
The conduction term V . (KVT) in equation (12) is negligible in most geothermal reservoir problems. However, it is significant in large-scale slow convective motion and therefore is included for completeness.
The functional dependence of the thermodynamic and transport quantities in equations (1 1) and (12) on p and h is not given in detail here. The formulae given by Mercer and Faust’ were used and the reader can refer to that paper for details.
The mixture enthalpy h is related to the other variables by the formula
h = ( P V h V S V +PlhlSI)I(P,S, + PISI) (16)
The source and sink terms in equations (1) and (4) arise from bores, and assuming that the total mass withdrawal qm for each bore is specified, then the energy withdrawal qe is determined as follows:
4 e = q v h v + q l h l (17)
358 G . A. ZYVOLOSKI, M. J . OSULLIVAN AND D. E. KROL
The form of equation (18) shows the importance of the ratio of relative permeabilities R I / R v in controlling the discharge composition. The limited field data available (see G r a d 3 ) show that there may be significant differences between the relative permeability formulae used here in equations (4) and (10) (and by other authors) and field permeabilities.
NUMERICAL SOLUTION-AD1 METHOD
The first method considered in this paper is based on the alternating direction implicit method (see Richtmyer and Morton14). Application of AD1 splitting to the finite difference approxima- tion of equations (11) and (12) gives the following formulae for a two-dimensional problem:
1 1 n + 1 / 2 n + l n + l n+1/2 n + l F::' .r (A::' - AZi,J -- [Da,+1/2, (P i+1,j -P i,j - Da,-,/,.j - P i-1.j)I A x
for LY = m and CY = e.
Here n + 1 / 2 -1 P:, = P ( ~ A X , JAY, n ~ t ) , D ~ , . ~ + ~ , ~ - 2 (Z,::112 +D:, .~+~/~)
D:1.1+1/2 = m:,.,+l +D:,.,)
and
for example. In equations (19) only half of the overall AD1 scheme is shown. At this time step the variables in the x derivatives are treated implicitly and the y derivatives are explicit. At the next time step the implicit and explicit directions are reversed as follows:
for CY = m and CY =e. The equations (19) and (20) for all the nodes considered can each be written in the form
where {F,} and {F,} are vectors containing elements F::,' and Fz,zl respectively for equation (19) and FkT; and F,",:' for equation (20).
To solve the non-linear difference equations (21) arising at the n + 1 th time step (similarly for n + 2 th time step) the Newton-Raphson method is used to obtain corrections Ap:,, Ah:, to
h:'lVk +Ah:,. This process is started with the initial estimates =p: , , hl., = h:, and
n + l . k + l - n + l . k +A~:] , h;,+l.k+l = iteratively adjust the values of p:T1 and h:,+' using p I . , - PI., n c l . 0
MODELLING GEOTHERMAL RESERVOIRS 359
terminated when convergence to a sufficient level of accuracy is obtained ((F,,,'sl < IFe's( < lo-'). The equationssolved at each iteration can be represented in block form as follows:
Here matrices such as [aF,,,/dp] contain the derivatives of F,,,., with respect to all the variables pimi. The conduction terms V . (KVT) have been omitted from equations (19) and (20). However, their inclusion is a straightforward matter. In equation (22) the contribution to the coefficient matrix from the conduction terms is neglected and therefore the effect on the solution process of the including of the conduction terms merely corresponds to additional terms added to the energy source term 4e.
To obtain a computationally efficient scheme for solving (22) the derivatives of the trans- missibilities with respect to pressure and enthalpy are neglected. With this very important modification the matrices [@',,,/ah] and [aF,/ah] become diagonal while [aFm/ap] and [dFe/ap] are tridiagonal. This simplification means that the equations for the correction {Ah} and {Ap} can be separated, giving the equation
{Ah} = [ %] -' ( - {Fe} + [ 3] {Ap}] aP
and
As a consequence of the diagonal nature of [dF,/dh] and [aFe/ah], equation (24) represents a tridiagonal system which is easily solved using the Thomas algorithm. Once {Ap} is known {Ah} is obtained by direct evaluation using equation (24). Although there is an approximation made in the solution process by neglecting the derivatives of the transmissibilities with respect t o p and h it is important to note that this involves no approximation in the final solution of (21). The coefficient matrix and the right-hand side of (22) are updated at each iteration using the latest information to calculate transmissibilities and other quantities.
Some modification is possible to the above scheme to improve it. For example, upstream weighting of the transmissibilities is often recommended (see C r i c h l o ~ ' ~ ) and does not significantly change the equations (19) and (20) or the overall solution procedure. Also some economy is possible in solving equation (22) by updating the coefficient matrix only every few iterations. However, the saving in time for function calculations in the derivatives is often offset by a requirement for an increased number of iterations to reach convergence, and in any case the procedure usually converges in two or three iterations (except on rows or columns including points when a phase change occurs when four or five iterations are required).
NUMERICAL SOLUTION-HOPSCOTCH METHOD
The second method used is a variation of the hopscotch method. The hopscotch technique involves alternate applications of the explicit and implicit finite difference approximations of equations (1 1) and (12). The second implicit step ensures that the hopscotch method does not have a restriction on the time step size for stability, and a special scanning of the mesh (suggested
360 G . A. ZWOLOSKI, M. J. O'SULLIVAN AND D. E. KROL
by the name), alternated at each time step, means that the equations involve unknowns from only one mesh point at a time. A comprehensive description of the technique is given by Gourlay3 and all details are not repeated here. For equations (11) and (12) the process involves first a solution of the explicit equations below on points for which n + i + j is odd:
for a = m and a =e. This is followed by a solution of an implicit difference approximation of (1 1) and (12) on points for which n + i + j is even:
for a = m and a =e. Because of the 'hopscotch' design of the scanning, the grid equations (26) above involve only two unknowns p;;', hyf' . All the other advanced time level quantities are already known from equations (25). Thus the solution of (25) and (26) involves solving only two coupled non-linear equations in two unknowns at each node. This is accomplished using the Newton-Raphson procedure (the scalar version of equation (22)).
The one disadvantage of the hopscotch method is that it introduces systematic errors in the overall balance of mass and energy in the system. This arises when the solution is changing rapidly and occurs because energy and mass withdrawals from sinks tend to be half a time step ahead of the redistribution of the remaining quantities. However, as the results in the next section show this problem does not produce an accumulating error. In the problems solved so far it was found that the Newton-Raphson iteration converged in two or three iterations at all points except when a phase change was occurring, in which case four or five were required.
As for the AD1 method inclusion of the conduction terms presents no difficulty. The only important point to note is that their contribution to derivatives of F:,:' with respect to the p's and h's is neglected in the Newton-Raphson solution process. However, their contribution to the calculation of FZ,:' at each iteration is included.
MODEL PROBLEMS
The first example to which these models are applied is a hypothetical one-dimensional reservoir previously analysed by Mercer and Faust7 and Faust and Mercer.' The data for this problem are given in Table I and a sketch of the reservoir is shown in Figure 1. The results obtained by using the AD1 method for a five year simulation are presented in Table I1 and similar results for the hopscotch method are given in Table 111. As can be seen from Tables I1 and I11 the results for the hopscotch method and the AD1 method are nearly identical after five years and very close throughout the simulation. The results for the AD1 method are also very close to those of Mercer and Faust7 and Faust and Mercer.'
MODELLING GEOTHERMAL RESERVOIRS 361
The second example considered is a hypothetical two-dimensional reservoir proposed by Faust and Mercer.' Data used for the simulation is given in Table IV and a sketch of the reservoir is given in Figure 2. The grid layout and values of pressure, enthalpy and saturation after 510 days are given in Table V. Comparisons with the results of Faust and Mercer' are not given because their results appear to be in error showing a drying out of the field adjacent to the boundaries furthest from the sink rather than adjacent to the boundaries closest to the sink. In
Table I. Parameters for the one-dimensional model problem
Parameter Symbol Value
Permeability k 1.0 x 10-'*m2
Porosity 4 0.10 Discharge q m 20 kg/s
Thermal conductivity K 3-20 W/m.K Rock density Pr 2500 kg/m'
Aquifer length - 2000 m Aquifer cross-section area - 2.5 x los m2 Initial pressure P? 43.8 bar Initial enthalpy h? 1-02 MJ/kg Initial temperature T? 238°C Mesh spacing Ax 200 m Time step At 1 day (until 80 days)
2 days (after 80 days)
Rock specific heat C, 1.01 kJ1kg.K
~~
Table 11. Results for the one-dimensional model problem-AD1 method
Table VI the results of a simulation of the same problem using the hopscotch method are given and in Table VII results for upstream weighting of transmissibilities are shown. The results shown in Table VIII were obtained by using a finer mesh (12 X 8 rather than 6 X 4). The results are in good agreement.
In each case the total mass balance error is given. This was calculated as (Q(t)-MR(O)+ MR( t ) ) /Q( t ) , where Q(t) is the total mass discharge from the reservoir up to time t and MR(f) is the total mass in the reservoir at time f.
Table IV. Parameters for the two-dimensional model problem
Parameter Symbol Value
Permeability Discharge Porosity Thermal conductivity Rock density Rock specific heat Aquifer length Aquifer width Aquifer thickness Initial pressure Initial enthalpy Initial temperature Mesh spacing Time step
k 1 . 0 ~ lo-’* m2 4 m 4 kg/s 4 0.1 K 2.5 W/m.K P r 2500 kg/m3 C, 1.01 kJ/kg.K - 600 m - 400 m - 100 m d, 69.0 bar h :j 1.05 MJ/kg T11, 242.8”C
Ax, AY At
100 m or 50 m 1 day (until 80 days) 2 days (after 80 days)
MODELLING GEOTHERMAL RESERVOIRS 363
qln* - 5
‘ 5 _ i l _ , 1 - 1 r.2 1 - 3 i-4 1 - 5 1-6 r = 7 1.0 1-9 i-10
Figure 1. The grid layout for the one-dimensional model problem. The total discharge, qm = 20 kg/s, is taken from the first cell
Table V. Results after 5 10 days using the AD1 method. Total mass balance error -0.6 x
p ( W 35.16 35.16 35.12 34-04 35.12 35.16 h(MJ/kg) 1-048 1.048 1,049 1.053 1.049 1.048 j = 4 SI 1.000 0.999 0.957 0.723 0.957 0.999
The test problems considered here both correspond to the relatively slow exploitation of a reservoir. The methods have also been used for short-time pressure transient calculations and there the advantage of upstream weighting of the transmissibilities is more noticeable.
DISCUSSION
Of the two methods presented the hopscotch technique is more efficient (approximately 40 per cent faster than the AD1 method for the 2-D model problem) but because of its errors in overall mass and energy balances it should be used with caution when an accurate simulation of a short-term transient effect is required. Further work is being carried out by the authors on developing procedures for systematically adjusting the mass balance errors in the hopscotch method.
Table VII. Results after 510 days using the AD1 method and upstream weighting of trans- missibilities. Total mass balance error -0.5 x lo-’
d b a r ) 35.16 35.16 35.14 34.65 35.14 35.16 h (MJ/kg) 1.048 1.048 1.048 1.063 1.048 1.048 j = 4 SI 1.000 0.999 0.974 0.657 0.974 0.999
Table VIII. Results after 510 days using the AD1 method and upstream weighting of trans- missibilities on a 12 x 8 grid. Results shown are calculated by averaging over four blocks. Total
mass balance error -0.3 x lo-’
p ( W 35.16 35.16 35.11 34.50 35.11 35.16 h(MJ/kg) 1.048 1.048 1.050 1.054 1.050 1.048 j = 4 SI 1.000 1.000 0.938 0.715 0.938 1.000
There is no difficulty in extending the methods presented here to non-homogeneous and non-isotropic porous media. However, because of the nature of the finite difference methods used the reservoir must be represented by a grid of rectangular cells (the overall shape of the reservoir need not be rectangular).
The use of the Newton-Raphson correction to achieve accuracy together with the use of the AD1 method or the hopscotch method to achieve efficiency in handling a large number of nodes means that the methods described above are capable of handling the long-term simulation of large-scale, complex geothermal flows. Although only two small model problems are considered here they are of sufficient generality to provide a reasonable test of the methods, particularly since they involve a phase change from one-phase to two-phase conditions.
REFERENCES
1. V. K. Saul’yev, Integration of Equations of Parabolic Type by the Method of Nets, MacMillan, New York, 1964. 2. P. Gordon, ‘Nonsymmetric difference equations’, J. SOC. Indust. Appl . Math. 13,667-673 (1965). 3. A. R. Gourlay, ‘Hopscotch: a fast second-order partial differential equation solver’, J. Inst. Maths. Applics. 6,
4. D. H. Brownell, Jr., S. K. Garg and J. W. Pritchett, ‘Computer simulation of geothermal reservoirs’, paperSPE5381 presented at 45th California Regional Meeting of the SOC. of Pet. Eng. of AIME, Ventura (1975).
5. S. K. Garg, J. W. Pritchett and D. H. Brownell, Jr., ‘Transport of mass and energy in porous media’, paper presented at the Second United Nations Symposium on the Development and Use of Geothermal Resources, San Francisco (1975).
6. J. W. Pritchett, S. K. Garg, D. H. Brownell, Jr. and H. B. Levine, ‘Geohydrological environmental effects of geothermal power production-phase I’, Report No. SSS-R- 75-2733, Systems, Science and Sofrware, La Jolla, California (1975).
7. J. W. Mercer and C. R. Faust, ‘Simulation of water- and vapor-dominated hydrothermal reservoirs’, paper SPE 5520 presented at 50th Annual Fall Meeting of the SOC. of Pet. Eng. of AIME, Dallas, Texas (1975).
8. C. R. Faust and J. W. Mercer, ‘Mathematical modeling of geothermal systems’, in Proc. Second United Nations Symposium on the Development and Use of Geothermal Resources, San Francisco (1975).
9. C. R. Faust and J. W. Mercer, ‘An analysis of finite-difference and finite-element techniques for geothermal reservoir simulation’, paper SPE 5742 presented at the Fourth Symposium of Numerical Simulation of Reservoir Performance of the SOC. of Pet. Eng. of AIME, Los Angeles, Calif. (1976).
375-390 (1970).
366 G. A. ZYVOLOSKI. M. J. O’SULLIVAN AND D. E. KROL
10. R. M. Toronyi and S. M. Farouq Ali, ‘Two-phase two-dimensional simulation of a geothermal reservoir and
11. J. W. Mercer, Jr., C. R. Faust and G. F. Pinder, ‘Geothermal reservoir simulation’; Proc. Conf. on Research for the
12. A. T. Corey, ‘The interrelation between gas and oil relative permeabilities’, Producers Monrhly, 19,3841 (1954). 13. M. A. Grant, ‘Permeability reduction factors at Wairakei’, paper 77-HT-52, presented at the A I M - A I M E Hear
14. R. D. Richtrnyer and K. W. Morton, Difference Merhods for Inirial- Value Problems, 2nd edition, Interscience, New
15. H. B. Crichlow, Modern Reseruoir Engineering-a Simulation Approach, Prentice-Hall, Englewood Cliffs, New
wellbore system’, Soc. Pet. Eng. J. 17, 171-183 (1977).
Development of Geothermal Energy Resources, Pasadena, California (1974).
