- Email: [email protected]
A pressure-enthalpy finite element model for simulating hydrothermal reservoirs
Mathematics and Computers 0 North-Holland Publishing
in Simulation Company
A PRESSUREENTHALPY
XX (1978)
167-178
FINITE ELEMENT MODEL FOR SIMULATING HYDROTHERMAL
RESERVOIRS Pongsarl HUYAKORN New Mexico Institute of Mining and Technology,
Socorro, New Mexico 87801,
USA
and George F. PINDER Ainceton
University, Princeton, New Jersey 08540,
USA
A numerical model is presented for simulating single or two-phase flow and energy transport in hydrothermal reservoirs. The model is formulated via two non-linear equations for fluid pressure and enthalpy. Both equations are solved simultaneously using a new finite element technique which employs asymmetric weighting functions to overcome numerical oscillation. Non-linearity is treated by a modified Newton-Raphson scheme which takes into account derivative discontinuities in the non-linear coefficients. This scheme also treats unknown flux boundary conditions inplicitly, thus allowing larger time steps to be taken without inducing instability. The proposed model is applied to two test examples involving one-dimensional flow in both hot water and steam dominated reservoirs. Results indicate that the numerical technique presented is efficient and the model can be used to simulate both types of reservoirs.
1. Introduction
problems he encountered when attempting the standard Galerkin finite element technique in conjunction with a simple iterative procedure to such a situation. These difficulties stem from the high non-linearities in the strongly coupled partial differential equations and discontinuities in the derivatives of their non-linear coefficients. The purpose of this paper is to present a pressureenthalpy model which employs a new finite element technique to overcome the above numerical difficulties. The new technique is formulated via a general weighted residual procedure and non-linear iteration is performed using a modified Newton-Raphson scheme which allows for slope discontinuities in the non-linear coefficients.
In recent years, there has been a number of attempts to develop numerical models for simulating the behavior of hydrothermal reservoirs. Several workers [l-8] have derived different sets of governing equations and have employed finite difference and finite element. methods to obtain approximate solutions. Among those who have developed single-phase flow models for simulating hot water reservoirs are Harlow and Pracht [ 11, Mercer [2] and Sorey [3]. Their models have limited use and are not applicable to vapor-dominated reservoirs. Recognizing the limitation of the single-phase models, Mercer and Faust [4] developed a more flexible model capable of simulating both hot-water and vapordominated reservoirs. In their model, the conventional Gale&in finite element method was employed to solve the non-linear equations written in terms of fluid pressure and enthalpy. Other similar models were also developed by Garg et al. [S] and Lasseter et al. 161. It was soon realized that serious numerical difficulties could be encountered in dealing with a situation involving rapid phase conversion. Faust [7] reported several
2. Basic equations The equations employed to simulate a hydrothermal reservoir comprising both steam and water or either water or steam have been derived by Mercer and Faust [4] with an assumption that capillary effects can be 167
P. Huyakorn,
168
G.F. Pinder / Model for hydrothermal
neglected when steam and water coexist. These equations can be written in the form
reservoirs I
0.6 -
I
I
h =669.3
I
JOULE/qm
- 0.9 I A
(1)
0.6 N P
B
: -0.4
(2)
$ COMPRESSED WATER ZONE
in which p is fluid pressure, h is specific enthalpy of the fluid mixture, 7, F, A, /3 and Hare non-linear coeffl. cients defined as follows: 1.4
r=rwtr8,
(3a)
F=(bp.
WI
: - 0.2
I.6 pXlO-’
I
0.4
I 2.6
30
dyne/cm2
I
h=689.3
ForhL
I 2.2
I
I
JOULE/qm
0.4
/===~o.3
..[
Forhh: h=Kg+rh,
(34
P=K$
(3e)
01
H=@~h+U -@PA.
(3f)
In eqs. (3a) and (3b), 7, and 7, are mobilities of the water and steam phases respectively, $ is the effective porosity of the porous medium and p is the density of the fluid mixture. The phase mobilities and mixture density are defined as cm = w, 4
Trn = Pmkk,mIAn P = &VP,
+ SSP,
9
,
(44 (4b)
where k is the specific permeability of the porous medium and pm, k,m, pm and S, are density, relative permeability, dynamic viscosity and saturation of fluid phase m respectively. In eqs. (3c-3f), K is thermal conductivity of the medium, T is temperature, hk and h,* are saturated enthalpies of water and steam respectively, and pr and h, are density and enthalpy of the rock matrix. The two dependent variables, p and h, specified as unknowns in eqs. (1) and (2) uniquely define the thermodynamic state of the system. Once p and h are known, the non-linear coefficients in these equations can be evaluated using the procedure described in
I
I
I
I
8
I
1.4
LB
2.2
2.6
3.:
pxlo-‘dyne/cm2
Fig. 1. Variation (h = 889.3 J/g).
of non-linear
coefficients
with pressure
Appendix A. Typical plots of the functions r, /3 and X and the functions F and H are illustrated in figs. 1 and 2 respectively. As can be seen, these functions are highly non-linear in the two-phase zone of the thermodynamic region. Furthermore, there exist sharp discontinuities in their derivatives at the boundary of the two-phase and single-phase zones. It is this feature which makes the geothermal problem very difficult to treat by conventional finite difference or Galerkin finite element methods when steam and water coexist and there is also significant phase conversion. In such a case, both eqs. (1) and (2) become highly non-linear and display hyperbolic character when the values of coefficients approach zero. The combination of high non-linearities and slope discontinuities usually results in convergence difficulties and places severe restrictions on the size of allowable time steps. On the other hand, the hyperbolic character of the two equations often leads to oscillatory behavior of the solution obtained. The numerical
P. Huyakorn,
I
169
reservoirs
1
I
~‘I.84X 107dy”e
0.6 -
G.F. Pinder /Model for hydrothermal
- 0.6
/cm2
-0.6
N 0 x
E-PHASE
ZONE
- 0.4
P x r(
c o,2_
0 650
COMPRESSED WATER ZONE
I
I
I
850
950
1050
-
h (JOULE /gm)
1
I
I p=l.64X107
4
dyne/cm2
2-PHASE
- 1.0
ZONE
-0.9
n x 0 0.1 -
0 650
COMPRESSED WATER ZONE
I 750
1
i
Fig. 3. Asymmetric
weighting
functions
for linear
1-D elements.
-0.2
750
04-
_I_
L
”
linear convection-diffusion equation and the non-linear two-phase flow equations (see refs. [9,10]). These functions can be expressed as a combination of the linear shape functions and quadratic functions that are dependent on the flow direction along each side of an element. To apply the proposed formulation, let trial solutions for pressure and enthalpy be expressed in the form P_&) >
(5)
h(Xj, t) = NJ(Xj) h_r(t) .
(6)
PC%, t) = N&i)
I
I
950
a50
I 1050
0.6 1150
h(JOULE/gm)
Fig. 2. Variation of non-linear @ = 1.84 X 10’ dyne/cm2).
coefficients
with enthalpy
oscillations are unacceptable because they can be greatly amplified through the calculation of the nonlinear coefficients.
3. Weighted residual finite element formulation In order to overcome the aforementioned difticulties, we propose a new finite element technique. In this technique, eqs. (1) and (2) are discretized via a general weighted residual formulation and the resulting set of non-linear algebraic equations are solved by a modified fully implicit Newton-Raphson scheme. In the proposed weighted residual formulation, the spatial terms in eqs. (1) and (2) are weighted using asymmetric weighting functions and the accumulation time derivative terms are weighted using the shape functions. For linear iso-parametric elements, the asymmetric functions (see fig. 3) have been successfully applied to the
Let us also assume that the non-linear be represented by
coefficients
can
P(xi, t) = NJ(xi>
PAtI
>
(7)
%G, t) = N&i)
F,(t)
9
(8)
Vxi,
t)
= NJ(Xi)
hJ(t>
>
(9)
PCxi,
f)
= NJ(Xi>
Pdt>
>
(10)
HXi, f) = NJ(%) ffJ(f) >
(11)
where NJ denotes a set of basis functions, pJ and hJ are nodal values of pressure and enthalpy respectively, TJ, FJ, AJ, PJ and HJ are nodal values of coefficients 7, F, X, /3and H respectively, and the nodal summation convection for subscript J, i.e.,
NJcxi> HJ(~)
= J$
NJ(Xt)
HAtI
.
Let WI be a complete set of asymmetric weighting functions. Application of the proposed weighted residual procedure to eqs. (1) and (2) leads to
170
P. Huyakorn,
G.F. Pinder /Model for hydrothermal
reservoirs
treated carefully to avoid numerical instability. On substituting eq. (16) into eq. (1 S), one obtains
s
_
iV$dR=O;
(13)
R
where R is the entire flow region. Using Greens’s theorem to transform the second derivative terms, one obtains
dR
s
_
B
I-W %.dB=O ‘&Xi ’
’
(18)
4. Treatment of non-linear equations
(14)
(ApJ + /3hJ) + NI E]dR
Eqs. (14) and (18) represent a system of non-linear equations which can be efficiently treated by the Newton-Raphson technique. This technique is modified slightly to allow for discontinuities in derivatives of the non-linear coefficients when phase conversion occurs. To employ the proposed technique, it is convenient to write eqs. (14) and (18) in residual form as follows:
(15) I
where B is the flow boundary and Yliis the outward unit normal vector. The boundary integral in eq. (14) can be identified as the fluid flux at boundary node 1. This flux may be either known and prescribed or unknown. The boundary integral term in eq. (15) can also be written in the form
1
T~$~~+~(~+~‘-F’)
RI-
dR=O,
I
(19)
GI =
+$
(Ht+At -ZIP)
dR=o, 1
t
s
W# $.
nidB , I
B
(16)
where the time derivative terms are approximated using finite differences and the boundary integrals are dropped. The non-linear boundary terms will be treated separately. On applying the Newton-Raphson procedure, one obtains
where aRI
@-gt!K
ap.
k+l
-p$~+~(lzJ (17)
The first and second terms on the RHS of eq. (16) can be identified as the conductive and convective heat fluxes at boundary node I, respectively. The conductive nodal heat flux is usually known but the convective nodal heat flux is always unknown and has to be
.I
-h3=-RF,
(21)
aGI
-P:)+-(h$+‘-h3=-Gf, ahJ
(22)
where k and k + 1 refer to the previous and current iterations at time level t + At. The derivatives of RI and GI are given by
P. Huyakorn,
W -=
G.F. Pinder /Model for hydrothermal
171
reservoirs
dR
@J
+
s-’ __-- 1*+A*dR NI
aFf+**
At
aPJ
R
aw,
a7
(26)
(27)
&J=d”-?k
NI aF*+** dR s nEah,
R At
where the terms in the square brackets are evaluated using the kth iterate and &J and ti J are pressure and enthalpy increments defined as
arj
i3Xi ahJ [ axi
+
a WI aNJ =-~--(?&3h~)dl-~NI(Hk-H*)dR, R axi axi
(234
(23b)
’
‘%J
=
hj+l - h$ .
(28)
Eqs. (2.5) and (26) can be put in matrix form as
=I -- _
(29)
&J where t$-[~]*~A*]dR+~~[~]*~A*dR,
(24a) fi=-jsTzp$dR-J%(Fk-F*)dR, Rax , &-
(24b)
substituting eqs. (19-20) and (23-24) into eqs. (2122) and linearizing the non-linear terms using known nodal values at the previous iteration, one obtains
s
[AIJI =
[
awZ
AaxJt axi
ah
(Fk - F*) dR , (25)
2
a@,
ab
bIJ
CIJ
dIJ
1
(3Oc)
’
The coefficients aIJ, bIJ, cIJ and dIJ can be readily extracted from eqs. (25) and (26).
R;-RIteI,
R
of flux boundary
conditions
(31)
(32)
ah
ap,axi apJaxi
aIJ
For boundary nodes which admit fluid and heat fluxes, the boundary integral terms exist and eqs. (19) and (20) have to be modified. It is convenient to write the modified equations in the form
sAt
aNJ k J -T----PJdRax axi
R
Gob)
3(Hk-H*)dR,
R At
5. Treatment
=_
(304
awI aNJ s --@&Ph$W R axi axi
‘$-[$]+Ar)dR+$$~:i_4idX.
R At
)
where
(33a)
P. Huyakorn,
172
G.F. Pinder /Model for hydrothermal
Wb) Normally, the fluid flux ~71and the conductive heat flux are prescribed. Thus, they merely appear on the RHS of the equations obtained by applying the Newton-Raphson procedure to R; and CT. On the other hand, the remaining boundary integral in eq. (32) is a non-linear term which should be treated implicitly to avoid numerical instability. Application of the Newton-Raphson procedure to eqs. (3 1) and (32) leads to the following matrix equation [A*]
(IO @
f”
=
Ah
7. Analysis of one-dimensional
where
(354
f;=fr-ciI? JWr0”
z
f’IJ ~A?rl =arJ 1 c;J
= CIJ
d;J
-
JWr g
wrae s ahJ B
case
A brief analysis is now presented for a one-dimensional case in which the equations obtained by applying the weighted residual procedure reduce to
1
dR=O,
(35c)
’
B
&J=dIJ-
Wb)
nidB , I
B
[ c;J
from these values. The increments are taken in such a manner that the possibility of crossing the phase boundary is avoided when the values pj and h$ are near saturated values. In brief, the derivatives with respect to the pressure variable are evaluated by taking a negative or backward increment if the point (p$, hJk) lies in the two-phase zone (see fig. 1) and by taking a positive or forward increment if the point ($3, h$) lies in the single-phase zone. On the other hand, the derivatives with respect to the enthalpy are evaluated by taking a positive increment if (ps, hf) lies in the two-phase zone (see fig. 2) and by taking a negative increment if it lies in the single-phase zone.
(34)
g*’
g; = gr - &+
reservoirs
2
nidB 1
2nidB aXi
The remaining coefficients
1
dwr WJ
JW#3" F
~~(XpJ+flhJ)+h$~J~
dR=O,
B
’
We)
are as given previously.
6. Solution procedure Eq. (29) or (34) represents the required set of linearized equations which can be solved by a standard Gaussian elimination scheme. After a new set of nodal values has been obtained, it is used to update the coefficient matrix and the RHS vector. The iteration cycle is then repeated for each time step until convergence to a prescribed tolerance is achieved. In evaluating the coefficient matrix, the derivatives of the nonlinear coefficients r, Fand their counterparts are approximated by chord slopes which are determined using nodal values p$ and h$ and taking very small increments
(36)
(37)
where Re is an element subregion and Xe denotes summation over all elements. For a typical linear element (see fig. 3), the weighting functions can be expressed in the form
W,(x) =Ni(x)
- 3CY$-‘i (
W,(X) = Nz(x) - 3cu $ (
1
,
(38)
,
_; )
in which (Yis a free parameter with a sign dependent on the flow direction and 1 is the element length. The sign of (Yis determined such that CYis positive if the flow velocity is in the positive direction. Substituting eqs. (38) and (39) into (36) and (37) and assuming that equal sized elements are used, one obtains the following equations for a particular node I in the network:
P. Huyakorn,
G.F. Pinder /Model for hydrothermal
173
reservoirs
8.1. Example 1
+I 6
@I--l 4dFI @I,1 (-+++at at at 1 ’ co
f (PI-PI-I)
(40)
+ :--flpI - P1+1)
+f (b - hI-I) + :@I
- &+I) (41)
This example was chosen to verify the finite element model by comparison with an analytical solution and experimental data available in the literature. The particular problem considered is illustrated in fig. 4. It concerns linear flow and thermal convection in a cylindrical sandstone core injected with hot water. The case simulated corresponds to the experimental condition modellecl by Arihara [ 121. For this case, the temperature distribution can be accurately described by the following linear equation: aT i3T PWUWCW z + -at [(I - $J)P&r + ~PpwC,l
where ?=i[(l
+(Y)71-1 +(I -a)711
7* = 1[(1 •t a) Tr + (1 - a)
,
q+11>
(42) (43)
and similar definitions hold for >;, X*, fi and p*. If the mass matrix is lumped by adding values along a row and placing the sum on its diagonal, the time derivatives in eqs. (40) and (41) become dFI/dt and dHI/dt respectively. The spatial terms in eq. (40) can be interpreted as the fluid flux towards node I, contributed by two adjacent elements. The effective value of the non-linear coefficient 7 for each element block is automatically expressed as a weighted sum of values at the upstream and downstream nodes. The role of parameter Q should be noted. For a= 0, the effective value is the average value normally obtained using the standard Galerkin procedure or conventional finite difference approximation. For (Y= 1, the effective value is the upstream value which can also be obtained from a finite difference approximation with one-point upstream weighting. Similar interpretation can be given to the corresponding terms in eq. (41).
+z(TTo
T-)=0,
(44)
subject to T(0, t) = Tj and T(x, 0) = To, where v, is the flow velocity p ,,,, c, and pr, c, are the density and specific heat at constant pressure of water and sandstone respectively, U is the overall heat transfer coefficient based on the radius r, of the core, and T, is the ambient temperature, The physical constants used by Arihara [ 121 in comparing the analytical solution with experimental data are given as follows: pr = 2.65 g/cm3 ,
c, = 0.88 J/g “C ,
Pw = 1 g/cm3,
c, = 4.187 J/g “C ,
u, = 0.0115 cm/s ,
$I = 0.22 )
2U/r, = 0.554 X 10e3 J/cm3 - s - “C , To = T, = 25.65 “C , Tj = f$ [48(1 - e-“~0387f/60) + 50 .51 2 where Ti is in “C and t is in seconds. To obtain our corresponding finite element solution, additional data on boundary conditions of the pressure
8. Numerical examples To demonstrate the applicability of the proposed numerical technique, two examples are presented. The first example concerns a single-phase flow problem and the second example concerns a more complex problem involving conversion from single-phase to twophase flow.
Fig. 4. Example
1, single-phase
flow in cylindrical
core.
P. Huyakom,
174
G.F. Pinder /Model for hydrothermal
reservoirs
8.2. Example 2
202
0
IO
20
30
40
50
60
I.ml
Fig. 5. Comparison of calculated and experimental temperature profiles.
This example was chosen to demonstrate the reservoir behavior during rapid phase conversion and to test our numerical technique under conditions of severe non-linearities in the governing equations and boundary conditions and slope discontinuities in the coefficients. The particular problem considered is illustrated in fig. 6 and is that of non-linear flow and heat transport in a cylindrical sandstone core initially filled with hot water and later subjected to rapid pressure drop at one end. It is assumed that the core is placed in an oven and that there is direct heat transfer which is given by
Q* =-F(Tand enthalpy variables were required. These are given as follows: ~(0, t)
= 1.725 X 10’ dyne/cm2
q(60, t) = 0.0115 g/(cm2s)
,
,
W, t> = M-i) J/g, where hi is an implicit function Tic -27.9144
t 0.301117hi
given by
t 3536.37/hi
- 4.31919 X lo-‘h? 1 . The flow region was discretized into 20 equally sized elements (Ax = 3 cm). The time discretization used wast=60secforO
T,).
(45)
The values of physical constants used in this simulation are as follows: pr = 2.65 g/cm3,
c, = 0.88 J/(g “C) ,
k= IO-*cm2,
$ = 0.36 ,
K = 0.021 J/(sec cm “C) ;
T, = 199°C)
2u _ = 1.04 X 10d3 J/(cmm3 set “C)
.
r0
Relative permeabilities’were of saturation as follows: kv
= S&t
assumed to be functions
k,, = (1 - S,)2
.
The initial and boundary conditions used are given in table 1. The flow region was again discretized into 20 equally sized elements. The time discretization was At = 10 set for t < 20 set, At = 2 set for 20 < t < 30 set and At = 10 set for t > 30 sec. (The smaller time step size was used during the phase conversion period.) Results were obtained by setting the upstream parameter (Yto unity for all elements and lumping the ele-
SURROUNDING
I------
60cm
TEMP.
‘& =199-C
-------A
Fig. 6. Example 2, single and two-phase flow in cylindrical core.
P. Huyakorn,
G. F. Pinder /Model for hydrothermal
Table 1 Initial and boundary
conditions
Boundary
data
Initial
Exit pressure
Initial pressure, p. = 2.17 x 10’ (dyne/cm*)
condition
for example
Time, t(sec)
pe(dyne/cm2)
0 10 20 30 40
2.17 1.731 1.25 1.07 0.97
x x x x x
reservoirs
175
2 condition
data
Initial enthalpy, ho = 812.4 - 0.2.x = 807.6 - 0.533x = 806 - 0.82(x - 27)
10’ 10’ 10’ 10’
ment mass matrix. The lumped mass matrix was used in preference to the consistent mass matrix as it was found that the latter spread the effect of the accumulation terms for one node over the adjacent nodes and eventually led to numerical ocillation and convergence difficulty. With the lumped mass matrix and full upstream weighting, no such numerical difficulty was encountered and the proposed Newton-Raphson technique led to rapid convergence. The convergence criterion used was a maximum norm of less than 0.05 %. It was found that this criterion gave satisfactory mass and energy balance checks. Fig. 7 shows saturation profiles at various times. As can be seen, the flow remains single-phase at t = 10 sets. and rapid phase conversion occurs during the period 20 < t < 30 sec. The phase conversion starts from the upstream end of the core and propagates very rapidly in the downstream direction. Figs. 8 and 9 show typical temperature and pressure profiles respectively. As can be expected, when the flow is partly single-phase and partly two-phase, the portion of a profile which corresponds to the single-
(0 < x =G24cm) (24 4 x Q 27cm) (27
Fig. 8. Calculated
0.6
I
0
IO
temperature
1
I
20
Fig. 9. Calculated
profiles.
I
30 x.cm
pressure
t
40
50
60
profiles.
c ----At=lOrec
,*I: 0
IO
-l%=Lrec
.25 0
20
Fig. 7. Calculated
30 x.cm
saturation
40
50
60
0
IO
t 20
I
I
I
30
40
50
I 60
r,cm
profiles.
Fig. 10. Comparison of saturation profiles coarse and refined time discretization.
obtained
using
176
P. Huyakorn,
G.F. Pinder /Model for hydrothermal
phase zone is merely a straight line with a negative slope. A study was also performed to assess stability of the proposed numerical scheme. Using the same mesh but a much coarser time discretization (At = 10 set) during phase conversion, a second set of results was obtained and compared with the first set. (In obtaining the second set of results, a critical stage was encountered at which the upstream node of the exit boundary element went two-phase while the downstream node still remained single-phase. When this happened, it was necessary to apply a phase constraint at the exit node by setting its enthalpy value to the saturated water enthalpy. Unless this was done, the exit nodal values would not converge or would converge to a wrong solution). Fig. 10 shows a comparison of saturation profiles at t = 30 sets and t = 40 sec. It may be noted that although there is a considerable difference between the two corresponding profiles at t = 30 set, the difference becomes much smaller at the subsequent time level. This indicates that the proposed scheme exhibits good stability characteristics. The large discrepancy in saturation values at t = 30 set is due to the fact that the time step size of At = 10 set is too large to accommodate rapid saturation change during 20 < t < 30 sec.
Charles R. Faust and James W. Mercer of the U.S. Geological Survey for helpful discussions and suggestions.
Appendix Evaluation of non-linear
h,*= 2822.82
The work reported in this paper was supported by the National Science Foundation Grant No. NSF-AER 74-01765. The authors would like to acknowledge
- 39.9521~ + 2.54342/p*
- 0.938879~~
,(A-1)
where hi is given in J/g. Knowing h& and hi, the thermodynamic state can be readily determined. (h < hk indicates compressed water zone, h; < h < hz indicates two-phase zone and h > h,” indicates superheated steam zone). 2. The known values of p and h are used to calculate water and steam densities (g/cm3). For compressed water zone pw = 0.989875
+ 4.00894
x 10-4p
+ 2.66608/h
+ 5.46283
x 10-7ph
(ifh
,
> 200)
For superheated PW
=
x 10-5h
- 1.29958
x 10-7h2
.
(A-3)
steam zone
-2.26162
x 1O-5 + 0.0438441~
+ 3.69276
- 4.00489
(A-2)
pw = 1 (for h $ 200)
x 10-8p4
+ 5.17644
- 1.79088
x 10-5ph
X lo-r3ph3.
(A-4)
For two-phase zone, h& is used in place of h in eq. (A-2) and hl is used in place of h in eq. (A-4). 3. Knowing saturated phase enthalpies and phase densities, water saturation and mixture density can be calculated from s,=
1
s,=o
P&
S,= Acknowledgement
coefficients
For specified values of pressure in (dyne/cm*) x lo-’ or J/cm3 and enthalpy in J/g, the non-linear coefficients 7, F, h, p and Hare determined in a computer subroutine which performs the following steps. 1. The known p value is used to calculate saturated enthalpy of water and steam, h; and h,*. The saturated enthalpy of water is obtained from linear interpolation using vector arrays containing values extracted from a steam table. The saturated enthalpy of steam, is calculated from the following regression formula (see ref. [4]).
9. Conclusions A pressure-enthalpy numerical model has been presented for simulating single or two-phase flow and energy transport in hydrothermal reservoirs. In this model, the equations for pressure and enthalpy are treated using a new finite element technique which employs (1) asymmetric weighting functions to overcome numerical oscillation (2) a modified NewtonRaphson scheme to treat non-linearities. The model has been applied to two one-dimensional test examples. Results indicate that the numerical technique used is efficient and effective in handling a situation involving high non-linearities and rapid phase conversion.
reservoirs
h(p, P = P\ysw
-h)
forh
(A-5)
forh
ah,*,
(A-6)
for h;
< h < h:,
(A-7)
- ps) - (&.,pw - h@s) + PJl
- S,)
4. The values of p and h are used to calculate degree celcius. For two-phase zone
(A-8) temperature
in
P. Huyakorn, 4667.0754
T=P@)
G.F. Pinder /Model for hydrothermal
273.15
12.598833-ln(lOp)
(A-9)
-
kr K
ni For compressed
water zone
T=lc,l@,h)+[p@)-~1@,h~)l, j, 1@, h) = -28.1515 + 3536.37/h
(A-l Oa)
- 0.137458~
For superheated-steam
Q
+ 0.301117h
x lo-‘h*
- 4.31919
+ 7.39386
+ 47.9921p
x 10-sh2
+ 0.03571541~~ - 2.26861
(A-l Ob)
R S T
(A-l la)
WI $J n
region
T=J,*(p,h)+[P@)-i*@,h~], ~j,*@, h) = -374.669
- 3.3372
- 1.1725
- 0.633606~~
P
7
x 106/(p2h2)
A, P
x 10-9h3p
x 1015/h4.
(A-lib)
The square bracket terms in eqs. (A-1Oa) and (A-l la) are introduced to keep the function T continuous at the boundaries of the thermodynamic zones. The derivatives of T can be derived from eqs. (A-9) to (A-l 1). 5. The value of temperature is used to compute viscosities of water and steam and the enthalpy of the solid matrix. nw=
10-e
x 239.4
/.I~ = 10-6(0.407T h, = 106(2041.
x 10[*48.37/(~+133.15)] + 80.4)
,
,
k nv =S$,
(A-15)
k,,=(l-Sw)*.
(A-16)
Eqs. (A-15) and (A-16) can be replaced by other types of functional relations depending on the properties of the porous medium. 7. Using the values calculated in steps 1 to 6, the non-linear coefficients 7, F, h, P and H are computed from eqs. (4a) and (3a) to (30. Note: The regression formulae given by eqs. (A-l), (A-2), (A-4), (A-lob), (A-lib) and (A-13) were obtained from ref. [4]. The formulae given by (A-9), (A-12) to (A-14) were obtained from ref. [ 131.
Nomenclature = = = = = =
flow boundary, non-linear coefficient specific enthalpy, saturated enthalpy, non-linear coefficient permeability,
= = = = = =
relative permeability, thermal conductivity, outward unit normal vector, shape functions, pressure, fluid flux, = heat flux, = flow region, = saturation, = temperature, = weighting functions, = porosity, = dynamic viscosity, = density, = non-linear coefficients in flow equation, = non-linear coefficients in heat equation.
i
= coordinate
r, J k m r t W
subscript
= nodal subscripts = iteration level = phase m = rock matrix = steam = time level = water
(A-14)
where clw and ps are given in g/(cm-set) and h, is given in J/g. 6. The values of Sw is used to calculate relative permeabilities of water and steam. For simplicity, it is assumed that
B F h h* H k
177
Indices
S
(A-12) (A-13)
+ 8.790,
NI P 4
reservoirs
References [l] F.H. Harlow and W.E. Pracht, “A theoretical
[2]
[3] [4]
[5]
[6]
in flow equations,
in heat equations,
[7]
study of geothermal energy extraction”, J. Geophys. Res., 77, 35, 7033-7048 (1972). J.W. Mercer, “Finite elements approach to the modelling of hydrothermal system”, Ph. D. Thesis, U. of Illinois, Urbana-Champaign (1973). M.L. Sorey, “Numerical modelling of liquid geothermal systems”, Ph. D. Thesis, U. of Calif., Berkeley (1975). J.W. Mercer and C.R. Faust, “Simulation of water- and vapor-dominated hydrothermal reservoirs”, paper SPE 5520 presented at 50th SPE Annual Fall Meeting, Dallas, Texas, (1975). S.K. Garg, J.W. Pritchett and D.H. Brownell, Jr., “Transport of mass and energy in porous media”, paper presented at 2nd United Nations Symposium on the Development and Use of Geothermal Resources, San Francisco, Calif. (1975). T.J. Lasseter, P.A. Witherspoon and M.J. Lippman, “The numerical simulation of heat and mass transfer in multidimensional two-phase geothermal reservoirs”, paper presented at 2nd United Symposium on the Development and Use of Geothermal Reservoirs, San Francisco, CA (1975). CR. Faust, “Numerical simulation of fluid flow and energy transport in liquid- and vapor-dominated hydro-
178
P. Huyakorn,
G.F. Pinder /Model for hydrothermal
thermal reservoirs”, Ph. D. Thesis, Pennsylvania State U., (1976). [8] R.M. Toronyi, “Two-phase two-dimensional simulation of a geothermal reservoir and well-bore system”, Ph. D. Thesis, Pennsylvania State U. (1974). [9] P.S. Huyakorn, “Solution of steady-state, convective transport equation using an upwind finite element scheme”, Applied Numerical Modelling, IPC Science and Technology Press, U.K., 1, 4, 187-195 (1977). [lo] P.S. Huyakorn and C. Taylor, “Solution of the transient convection-diffusion equation using extended upwind finite element scheme”, Applied Numerical Modelling, (in Press).
reservoirs
[l l] P.S. Huyakorn and G.F. Pinder, “A new finite element technique for the solution of two-phase flow through porous media”, Princeton U., Water Resources Program, Res. Report No. 76-WR-4 (1976). [ 121 N. Arihara, “A study of non-isothermal and two-phase flow through consolidated sandstone”, Ph. D. Thesis, Stanford U., Stanford, Calif., (1974). [13] H.J. Ramey (Jr.), W.E. Bragham, H.K. Chen, P.G. Aitkinson and N. Arihara, “Thermodynamic properties of hydrothermal systems”, Proc. of NSF Conf. on the Utilization of Volcano Energy, held in Hilo, Hawaii (1974).
in Simulation Company
A PRESSUREENTHALPY
XX (1978)
167-178
FINITE ELEMENT MODEL FOR SIMULATING HYDROTHERMAL
RESERVOIRS Pongsarl HUYAKORN New Mexico Institute of Mining and Technology,
Socorro, New Mexico 87801,
USA
and George F. PINDER Ainceton
University, Princeton, New Jersey 08540,
USA
A numerical model is presented for simulating single or two-phase flow and energy transport in hydrothermal reservoirs. The model is formulated via two non-linear equations for fluid pressure and enthalpy. Both equations are solved simultaneously using a new finite element technique which employs asymmetric weighting functions to overcome numerical oscillation. Non-linearity is treated by a modified Newton-Raphson scheme which takes into account derivative discontinuities in the non-linear coefficients. This scheme also treats unknown flux boundary conditions inplicitly, thus allowing larger time steps to be taken without inducing instability. The proposed model is applied to two test examples involving one-dimensional flow in both hot water and steam dominated reservoirs. Results indicate that the numerical technique presented is efficient and the model can be used to simulate both types of reservoirs.
1. Introduction
problems he encountered when attempting the standard Galerkin finite element technique in conjunction with a simple iterative procedure to such a situation. These difficulties stem from the high non-linearities in the strongly coupled partial differential equations and discontinuities in the derivatives of their non-linear coefficients. The purpose of this paper is to present a pressureenthalpy model which employs a new finite element technique to overcome the above numerical difficulties. The new technique is formulated via a general weighted residual procedure and non-linear iteration is performed using a modified Newton-Raphson scheme which allows for slope discontinuities in the non-linear coefficients.
In recent years, there has been a number of attempts to develop numerical models for simulating the behavior of hydrothermal reservoirs. Several workers [l-8] have derived different sets of governing equations and have employed finite difference and finite element. methods to obtain approximate solutions. Among those who have developed single-phase flow models for simulating hot water reservoirs are Harlow and Pracht [ 11, Mercer [2] and Sorey [3]. Their models have limited use and are not applicable to vapor-dominated reservoirs. Recognizing the limitation of the single-phase models, Mercer and Faust [4] developed a more flexible model capable of simulating both hot-water and vapordominated reservoirs. In their model, the conventional Gale&in finite element method was employed to solve the non-linear equations written in terms of fluid pressure and enthalpy. Other similar models were also developed by Garg et al. [S] and Lasseter et al. 161. It was soon realized that serious numerical difficulties could be encountered in dealing with a situation involving rapid phase conversion. Faust [7] reported several
2. Basic equations The equations employed to simulate a hydrothermal reservoir comprising both steam and water or either water or steam have been derived by Mercer and Faust [4] with an assumption that capillary effects can be 167
P. Huyakorn,
168
G.F. Pinder / Model for hydrothermal
neglected when steam and water coexist. These equations can be written in the form
reservoirs I
0.6 -
I
I
h =669.3
I
JOULE/qm
- 0.9 I A
(1)
0.6 N P
B
: -0.4
(2)
$ COMPRESSED WATER ZONE
in which p is fluid pressure, h is specific enthalpy of the fluid mixture, 7, F, A, /3 and Hare non-linear coeffl. cients defined as follows: 1.4
r=rwtr8,
(3a)
F=(bp.
WI
: - 0.2
I.6 pXlO-’
I
0.4
I 2.6
30
dyne/cm2
I
h=689.3
ForhL
I 2.2
I
I
JOULE/qm
0.4
/===~o.3
..[
Forh
(34
P=K$
(3e)
01
H=@~h+U -@PA.
(3f)
In eqs. (3a) and (3b), 7, and 7, are mobilities of the water and steam phases respectively, $ is the effective porosity of the porous medium and p is the density of the fluid mixture. The phase mobilities and mixture density are defined as cm = w, 4
Trn = Pmkk,mIAn P = &VP,
+ SSP,
9
,
(44 (4b)
where k is the specific permeability of the porous medium and pm, k,m, pm and S, are density, relative permeability, dynamic viscosity and saturation of fluid phase m respectively. In eqs. (3c-3f), K is thermal conductivity of the medium, T is temperature, hk and h,* are saturated enthalpies of water and steam respectively, and pr and h, are density and enthalpy of the rock matrix. The two dependent variables, p and h, specified as unknowns in eqs. (1) and (2) uniquely define the thermodynamic state of the system. Once p and h are known, the non-linear coefficients in these equations can be evaluated using the procedure described in
I
I
I
I
8
I
1.4
LB
2.2
2.6
3.:
pxlo-‘dyne/cm2
Fig. 1. Variation (h = 889.3 J/g).
of non-linear
coefficients
with pressure
Appendix A. Typical plots of the functions r, /3 and X and the functions F and H are illustrated in figs. 1 and 2 respectively. As can be seen, these functions are highly non-linear in the two-phase zone of the thermodynamic region. Furthermore, there exist sharp discontinuities in their derivatives at the boundary of the two-phase and single-phase zones. It is this feature which makes the geothermal problem very difficult to treat by conventional finite difference or Galerkin finite element methods when steam and water coexist and there is also significant phase conversion. In such a case, both eqs. (1) and (2) become highly non-linear and display hyperbolic character when the values of coefficients approach zero. The combination of high non-linearities and slope discontinuities usually results in convergence difficulties and places severe restrictions on the size of allowable time steps. On the other hand, the hyperbolic character of the two equations often leads to oscillatory behavior of the solution obtained. The numerical
P. Huyakorn,
I
169
reservoirs
1
I
~‘I.84X 107dy”e
0.6 -
G.F. Pinder /Model for hydrothermal
- 0.6
/cm2
-0.6
N 0 x
E-PHASE
ZONE
- 0.4
P x r(
c o,2_
0 650
COMPRESSED WATER ZONE
I
I
I
850
950
1050
-
h (JOULE /gm)
1
I
I p=l.64X107
4
dyne/cm2
2-PHASE
- 1.0
ZONE
-0.9
n x 0 0.1 -
0 650
COMPRESSED WATER ZONE
I 750
1
i
Fig. 3. Asymmetric
weighting
functions
for linear
1-D elements.
-0.2
750
04-
_I_
L
”
linear convection-diffusion equation and the non-linear two-phase flow equations (see refs. [9,10]). These functions can be expressed as a combination of the linear shape functions and quadratic functions that are dependent on the flow direction along each side of an element. To apply the proposed formulation, let trial solutions for pressure and enthalpy be expressed in the form P_&) >
(5)
h(Xj, t) = NJ(Xj) h_r(t) .
(6)
PC%, t) = N&i)
I
I
950
a50
I 1050
0.6 1150
h(JOULE/gm)
Fig. 2. Variation of non-linear @ = 1.84 X 10’ dyne/cm2).
coefficients
with enthalpy
oscillations are unacceptable because they can be greatly amplified through the calculation of the nonlinear coefficients.
3. Weighted residual finite element formulation In order to overcome the aforementioned difticulties, we propose a new finite element technique. In this technique, eqs. (1) and (2) are discretized via a general weighted residual formulation and the resulting set of non-linear algebraic equations are solved by a modified fully implicit Newton-Raphson scheme. In the proposed weighted residual formulation, the spatial terms in eqs. (1) and (2) are weighted using asymmetric weighting functions and the accumulation time derivative terms are weighted using the shape functions. For linear iso-parametric elements, the asymmetric functions (see fig. 3) have been successfully applied to the
Let us also assume that the non-linear be represented by
coefficients
can
P(xi, t) = NJ(xi>
PAtI
>
(7)
%G, t) = N&i)
F,(t)
9
(8)
Vxi,
t)
= NJ(Xi)
hJ(t>
>
(9)
PCxi,
f)
= NJ(Xi>
Pdt>
>
(10)
HXi, f) = NJ(%) ffJ(f) >
(11)
where NJ denotes a set of basis functions, pJ and hJ are nodal values of pressure and enthalpy respectively, TJ, FJ, AJ, PJ and HJ are nodal values of coefficients 7, F, X, /3and H respectively, and the nodal summation convection for subscript J, i.e.,
NJcxi> HJ(~)
= J$
NJ(Xt)
HAtI
.
Let WI be a complete set of asymmetric weighting functions. Application of the proposed weighted residual procedure to eqs. (1) and (2) leads to
170
P. Huyakorn,
G.F. Pinder /Model for hydrothermal
reservoirs
treated carefully to avoid numerical instability. On substituting eq. (16) into eq. (1 S), one obtains
s
_
iV$dR=O;
(13)
R
where R is the entire flow region. Using Greens’s theorem to transform the second derivative terms, one obtains
dR
s
_
B
I-W %.dB=O ‘&Xi ’
’
(18)
4. Treatment of non-linear equations
(14)
(ApJ + /3hJ) + NI E]dR
Eqs. (14) and (18) represent a system of non-linear equations which can be efficiently treated by the Newton-Raphson technique. This technique is modified slightly to allow for discontinuities in derivatives of the non-linear coefficients when phase conversion occurs. To employ the proposed technique, it is convenient to write eqs. (14) and (18) in residual form as follows:
(15) I
where B is the flow boundary and Yliis the outward unit normal vector. The boundary integral in eq. (14) can be identified as the fluid flux at boundary node 1. This flux may be either known and prescribed or unknown. The boundary integral term in eq. (15) can also be written in the form
1
T~$~~+~(~+~‘-F’)
RI-
dR=O,
I
(19)
GI =
+$
(Ht+At -ZIP)
dR=o, 1
t
s
W# $.
nidB , I
B
(16)
where the time derivative terms are approximated using finite differences and the boundary integrals are dropped. The non-linear boundary terms will be treated separately. On applying the Newton-Raphson procedure, one obtains
where aRI
@-gt!K
ap.
k+l
-p$~+~(lzJ (17)
The first and second terms on the RHS of eq. (16) can be identified as the conductive and convective heat fluxes at boundary node I, respectively. The conductive nodal heat flux is usually known but the convective nodal heat flux is always unknown and has to be
.I
-h3=-RF,
(21)
aGI
-P:)+-(h$+‘-h3=-Gf, ahJ
(22)
where k and k + 1 refer to the previous and current iterations at time level t + At. The derivatives of RI and GI are given by
P. Huyakorn,
W -=
G.F. Pinder /Model for hydrothermal
171
reservoirs
dR
@J
+
s-’ __-- 1*+A*dR NI
aFf+**
At
aPJ
R
aw,
a7
(26)
(27)
&J=d”-?k
NI aF*+** dR s nEah,
R At
where the terms in the square brackets are evaluated using the kth iterate and &J and ti J are pressure and enthalpy increments defined as
arj
i3Xi ahJ [ axi
+
a WI aNJ =-~--(?&3h~)dl-~NI(Hk-H*)dR, R axi axi
(234
(23b)
’
‘%J
=
hj+l - h$ .
(28)
Eqs. (2.5) and (26) can be put in matrix form as
=I -- _
(29)
&J where t$-[~]*~A*]dR+~~[~]*~A*dR,
(24a) fi=-jsTzp$dR-J%(Fk-F*)dR, Rax , &-
(24b)
substituting eqs. (19-20) and (23-24) into eqs. (2122) and linearizing the non-linear terms using known nodal values at the previous iteration, one obtains
s
[AIJI =
[
awZ
AaxJt axi
ah
(Fk - F*) dR , (25)
2
a@,
ab
bIJ
CIJ
dIJ
1
(3Oc)
’
The coefficients aIJ, bIJ, cIJ and dIJ can be readily extracted from eqs. (25) and (26).
R;-RIteI,
R
of flux boundary
conditions
(31)
(32)
ah
ap,axi apJaxi
aIJ
For boundary nodes which admit fluid and heat fluxes, the boundary integral terms exist and eqs. (19) and (20) have to be modified. It is convenient to write the modified equations in the form
sAt
aNJ k J -T----PJdRax axi
R
Gob)
3(Hk-H*)dR,
R At
5. Treatment
=_
(304
awI aNJ s --@&Ph$W R axi axi
‘$-[$]+Ar)dR+$$~:i_4idX.
R At
)
where
(33a)
P. Huyakorn,
172
G.F. Pinder /Model for hydrothermal
Wb) Normally, the fluid flux ~71and the conductive heat flux are prescribed. Thus, they merely appear on the RHS of the equations obtained by applying the Newton-Raphson procedure to R; and CT. On the other hand, the remaining boundary integral in eq. (32) is a non-linear term which should be treated implicitly to avoid numerical instability. Application of the Newton-Raphson procedure to eqs. (3 1) and (32) leads to the following matrix equation [A*]
(IO @
f”
=
Ah
7. Analysis of one-dimensional
where
(354
f;=fr-ciI? JWr0”
z
f’IJ ~A?rl =arJ 1 c;J
= CIJ
d;J
-
JWr g
wrae s ahJ B
case
A brief analysis is now presented for a one-dimensional case in which the equations obtained by applying the weighted residual procedure reduce to
1
dR=O,
(35c)
’
B
&J=dIJ-
Wb)
nidB , I
B
[ c;J
from these values. The increments are taken in such a manner that the possibility of crossing the phase boundary is avoided when the values pj and h$ are near saturated values. In brief, the derivatives with respect to the pressure variable are evaluated by taking a negative or backward increment if the point (p$, hJk) lies in the two-phase zone (see fig. 1) and by taking a positive or forward increment if the point ($3, h$) lies in the single-phase zone. On the other hand, the derivatives with respect to the enthalpy are evaluated by taking a positive increment if (ps, hf) lies in the two-phase zone (see fig. 2) and by taking a negative increment if it lies in the single-phase zone.
(34)
g*’
g; = gr - &+
reservoirs
2
nidB 1
2nidB aXi
The remaining coefficients
1
dwr WJ
JW#3" F
~~(XpJ+flhJ)+h$~J~
dR=O,
B
’
We)
are as given previously.
6. Solution procedure Eq. (29) or (34) represents the required set of linearized equations which can be solved by a standard Gaussian elimination scheme. After a new set of nodal values has been obtained, it is used to update the coefficient matrix and the RHS vector. The iteration cycle is then repeated for each time step until convergence to a prescribed tolerance is achieved. In evaluating the coefficient matrix, the derivatives of the nonlinear coefficients r, Fand their counterparts are approximated by chord slopes which are determined using nodal values p$ and h$ and taking very small increments
(36)
(37)
where Re is an element subregion and Xe denotes summation over all elements. For a typical linear element (see fig. 3), the weighting functions can be expressed in the form
W,(x) =Ni(x)
- 3CY$-‘i (
W,(X) = Nz(x) - 3cu $ (
1
,
(38)
,
_; )
in which (Yis a free parameter with a sign dependent on the flow direction and 1 is the element length. The sign of (Yis determined such that CYis positive if the flow velocity is in the positive direction. Substituting eqs. (38) and (39) into (36) and (37) and assuming that equal sized elements are used, one obtains the following equations for a particular node I in the network:
P. Huyakorn,
G.F. Pinder /Model for hydrothermal
173
reservoirs
8.1. Example 1
+I 6
@I--l 4dFI @I,1 (-+++at at at 1 ’ co
f (PI-PI-I)
(40)
+ :--flpI - P1+1)
+f (b - hI-I) + :@I
- &+I) (41)
This example was chosen to verify the finite element model by comparison with an analytical solution and experimental data available in the literature. The particular problem considered is illustrated in fig. 4. It concerns linear flow and thermal convection in a cylindrical sandstone core injected with hot water. The case simulated corresponds to the experimental condition modellecl by Arihara [ 121. For this case, the temperature distribution can be accurately described by the following linear equation: aT i3T PWUWCW z + -at [(I - $J)P&r + ~PpwC,l
where ?=i[(l
+(Y)71-1 +(I -a)711
7* = 1[(1 •t a) Tr + (1 - a)
,
q+11>
(42) (43)
and similar definitions hold for >;, X*, fi and p*. If the mass matrix is lumped by adding values along a row and placing the sum on its diagonal, the time derivatives in eqs. (40) and (41) become dFI/dt and dHI/dt respectively. The spatial terms in eq. (40) can be interpreted as the fluid flux towards node I, contributed by two adjacent elements. The effective value of the non-linear coefficient 7 for each element block is automatically expressed as a weighted sum of values at the upstream and downstream nodes. The role of parameter Q should be noted. For a= 0, the effective value is the average value normally obtained using the standard Galerkin procedure or conventional finite difference approximation. For (Y= 1, the effective value is the upstream value which can also be obtained from a finite difference approximation with one-point upstream weighting. Similar interpretation can be given to the corresponding terms in eq. (41).
+z(TTo
T-)=0,
(44)
subject to T(0, t) = Tj and T(x, 0) = To, where v, is the flow velocity p ,,,, c, and pr, c, are the density and specific heat at constant pressure of water and sandstone respectively, U is the overall heat transfer coefficient based on the radius r, of the core, and T, is the ambient temperature, The physical constants used by Arihara [ 121 in comparing the analytical solution with experimental data are given as follows: pr = 2.65 g/cm3 ,
c, = 0.88 J/g “C ,
Pw = 1 g/cm3,
c, = 4.187 J/g “C ,
u, = 0.0115 cm/s ,
$I = 0.22 )
2U/r, = 0.554 X 10e3 J/cm3 - s - “C , To = T, = 25.65 “C , Tj = f$ [48(1 - e-“~0387f/60) + 50 .51 2 where Ti is in “C and t is in seconds. To obtain our corresponding finite element solution, additional data on boundary conditions of the pressure
8. Numerical examples To demonstrate the applicability of the proposed numerical technique, two examples are presented. The first example concerns a single-phase flow problem and the second example concerns a more complex problem involving conversion from single-phase to twophase flow.
Fig. 4. Example
1, single-phase
flow in cylindrical
core.
P. Huyakom,
174
G.F. Pinder /Model for hydrothermal
reservoirs
8.2. Example 2
202
0
IO
20
30
40
50
60
I.ml
Fig. 5. Comparison of calculated and experimental temperature profiles.
This example was chosen to demonstrate the reservoir behavior during rapid phase conversion and to test our numerical technique under conditions of severe non-linearities in the governing equations and boundary conditions and slope discontinuities in the coefficients. The particular problem considered is illustrated in fig. 6 and is that of non-linear flow and heat transport in a cylindrical sandstone core initially filled with hot water and later subjected to rapid pressure drop at one end. It is assumed that the core is placed in an oven and that there is direct heat transfer which is given by
Q* =-F(Tand enthalpy variables were required. These are given as follows: ~(0, t)
= 1.725 X 10’ dyne/cm2
q(60, t) = 0.0115 g/(cm2s)
,
,
W, t> = M-i) J/g, where hi is an implicit function Tic -27.9144
t 0.301117hi
given by
t 3536.37/hi
- 4.31919 X lo-‘h? 1 . The flow region was discretized into 20 equally sized elements (Ax = 3 cm). The time discretization used wast=60secforO
T,).
(45)
The values of physical constants used in this simulation are as follows: pr = 2.65 g/cm3,
c, = 0.88 J/(g “C) ,
k= IO-*cm2,
$ = 0.36 ,
K = 0.021 J/(sec cm “C) ;
T, = 199°C)
2u _ = 1.04 X 10d3 J/(cmm3 set “C)
.
r0
Relative permeabilities’were of saturation as follows: kv
= S&t
assumed to be functions
k,, = (1 - S,)2
.
The initial and boundary conditions used are given in table 1. The flow region was again discretized into 20 equally sized elements. The time discretization was At = 10 set for t < 20 set, At = 2 set for 20 < t < 30 set and At = 10 set for t > 30 sec. (The smaller time step size was used during the phase conversion period.) Results were obtained by setting the upstream parameter (Yto unity for all elements and lumping the ele-
SURROUNDING
I------
60cm
TEMP.
‘& =199-C
-------A
Fig. 6. Example 2, single and two-phase flow in cylindrical core.
P. Huyakorn,
G. F. Pinder /Model for hydrothermal
Table 1 Initial and boundary
conditions
Boundary
data
Initial
Exit pressure
Initial pressure, p. = 2.17 x 10’ (dyne/cm*)
condition
for example
Time, t(sec)
pe(dyne/cm2)
0 10 20 30 40
2.17 1.731 1.25 1.07 0.97
x x x x x
reservoirs
175
2 condition
data
Initial enthalpy, ho = 812.4 - 0.2.x = 807.6 - 0.533x = 806 - 0.82(x - 27)
10’ 10’ 10’ 10’
ment mass matrix. The lumped mass matrix was used in preference to the consistent mass matrix as it was found that the latter spread the effect of the accumulation terms for one node over the adjacent nodes and eventually led to numerical ocillation and convergence difficulty. With the lumped mass matrix and full upstream weighting, no such numerical difficulty was encountered and the proposed Newton-Raphson technique led to rapid convergence. The convergence criterion used was a maximum norm of less than 0.05 %. It was found that this criterion gave satisfactory mass and energy balance checks. Fig. 7 shows saturation profiles at various times. As can be seen, the flow remains single-phase at t = 10 sets. and rapid phase conversion occurs during the period 20 < t < 30 sec. The phase conversion starts from the upstream end of the core and propagates very rapidly in the downstream direction. Figs. 8 and 9 show typical temperature and pressure profiles respectively. As can be expected, when the flow is partly single-phase and partly two-phase, the portion of a profile which corresponds to the single-
(0 < x =G24cm) (24 4 x Q 27cm) (27
Fig. 8. Calculated
0.6
I
0
IO
temperature
1
I
20
Fig. 9. Calculated
profiles.
I
30 x.cm
pressure
t
40
50
60
profiles.
c ----At=lOrec
,*I: 0
IO
-l%=Lrec
.25 0
20
Fig. 7. Calculated
30 x.cm
saturation
40
50
60
0
IO
t 20
I
I
I
30
40
50
I 60
r,cm
profiles.
Fig. 10. Comparison of saturation profiles coarse and refined time discretization.
obtained
using
176
P. Huyakorn,
G.F. Pinder /Model for hydrothermal
phase zone is merely a straight line with a negative slope. A study was also performed to assess stability of the proposed numerical scheme. Using the same mesh but a much coarser time discretization (At = 10 set) during phase conversion, a second set of results was obtained and compared with the first set. (In obtaining the second set of results, a critical stage was encountered at which the upstream node of the exit boundary element went two-phase while the downstream node still remained single-phase. When this happened, it was necessary to apply a phase constraint at the exit node by setting its enthalpy value to the saturated water enthalpy. Unless this was done, the exit nodal values would not converge or would converge to a wrong solution). Fig. 10 shows a comparison of saturation profiles at t = 30 sets and t = 40 sec. It may be noted that although there is a considerable difference between the two corresponding profiles at t = 30 set, the difference becomes much smaller at the subsequent time level. This indicates that the proposed scheme exhibits good stability characteristics. The large discrepancy in saturation values at t = 30 set is due to the fact that the time step size of At = 10 set is too large to accommodate rapid saturation change during 20 < t < 30 sec.
Charles R. Faust and James W. Mercer of the U.S. Geological Survey for helpful discussions and suggestions.
Appendix Evaluation of non-linear
h,*= 2822.82
The work reported in this paper was supported by the National Science Foundation Grant No. NSF-AER 74-01765. The authors would like to acknowledge
- 39.9521~ + 2.54342/p*
- 0.938879~~
,(A-1)
where hi is given in J/g. Knowing h& and hi, the thermodynamic state can be readily determined. (h < hk indicates compressed water zone, h; < h < hz indicates two-phase zone and h > h,” indicates superheated steam zone). 2. The known values of p and h are used to calculate water and steam densities (g/cm3). For compressed water zone pw = 0.989875
+ 4.00894
x 10-4p
+ 2.66608/h
+ 5.46283
x 10-7ph
(ifh
,
> 200)
For superheated PW
=
x 10-5h
- 1.29958
x 10-7h2
.
(A-3)
steam zone
-2.26162
x 1O-5 + 0.0438441~
+ 3.69276
- 4.00489
(A-2)
pw = 1 (for h $ 200)
x 10-8p4
+ 5.17644
- 1.79088
x 10-5ph
X lo-r3ph3.
(A-4)
For two-phase zone, h& is used in place of h in eq. (A-2) and hl is used in place of h in eq. (A-4). 3. Knowing saturated phase enthalpies and phase densities, water saturation and mixture density can be calculated from s,=
1
s,=o
P&
S,= Acknowledgement
coefficients
For specified values of pressure in (dyne/cm*) x lo-’ or J/cm3 and enthalpy in J/g, the non-linear coefficients 7, F, h, p and Hare determined in a computer subroutine which performs the following steps. 1. The known p value is used to calculate saturated enthalpy of water and steam, h; and h,*. The saturated enthalpy of water is obtained from linear interpolation using vector arrays containing values extracted from a steam table. The saturated enthalpy of steam, is calculated from the following regression formula (see ref. [4]).
9. Conclusions A pressure-enthalpy numerical model has been presented for simulating single or two-phase flow and energy transport in hydrothermal reservoirs. In this model, the equations for pressure and enthalpy are treated using a new finite element technique which employs (1) asymmetric weighting functions to overcome numerical oscillation (2) a modified NewtonRaphson scheme to treat non-linearities. The model has been applied to two one-dimensional test examples. Results indicate that the numerical technique used is efficient and effective in handling a situation involving high non-linearities and rapid phase conversion.
reservoirs
h(p, P = P\ysw
-h)
forh
(A-5)
forh
ah,*,
(A-6)
for h;
< h < h:,
(A-7)
- ps) - (&.,pw - h@s) + PJl
- S,)
4. The values of p and h are used to calculate degree celcius. For two-phase zone
(A-8) temperature
in
P. Huyakorn, 4667.0754
T=P@)
G.F. Pinder /Model for hydrothermal
273.15
12.598833-ln(lOp)
(A-9)
-
kr K
ni For compressed
water zone
T=lc,l@,h)+[p@)-~1@,h~)l, j, 1@, h) = -28.1515 + 3536.37/h
(A-l Oa)
- 0.137458~
For superheated-steam
Q
+ 0.301117h
x lo-‘h*
- 4.31919
+ 7.39386
+ 47.9921p
x 10-sh2
+ 0.03571541~~ - 2.26861
(A-l Ob)
R S T
(A-l la)
WI $J n
region
T=J,*(p,h)+[P@)-i*@,h~], ~j,*@, h) = -374.669
- 3.3372
- 1.1725
- 0.633606~~
P
7
x 106/(p2h2)
A, P
x 10-9h3p
x 1015/h4.
(A-lib)
The square bracket terms in eqs. (A-1Oa) and (A-l la) are introduced to keep the function T continuous at the boundaries of the thermodynamic zones. The derivatives of T can be derived from eqs. (A-9) to (A-l 1). 5. The value of temperature is used to compute viscosities of water and steam and the enthalpy of the solid matrix. nw=
10-e
x 239.4
/.I~ = 10-6(0.407T h, = 106(2041.
x 10[*48.37/(~+133.15)] + 80.4)
,
,
k nv =S$,
(A-15)
k,,=(l-Sw)*.
(A-16)
Eqs. (A-15) and (A-16) can be replaced by other types of functional relations depending on the properties of the porous medium. 7. Using the values calculated in steps 1 to 6, the non-linear coefficients 7, F, h, P and H are computed from eqs. (4a) and (3a) to (30. Note: The regression formulae given by eqs. (A-l), (A-2), (A-4), (A-lob), (A-lib) and (A-13) were obtained from ref. [4]. The formulae given by (A-9), (A-12) to (A-14) were obtained from ref. [ 131.
Nomenclature = = = = = =
flow boundary, non-linear coefficient specific enthalpy, saturated enthalpy, non-linear coefficient permeability,
= = = = = =
relative permeability, thermal conductivity, outward unit normal vector, shape functions, pressure, fluid flux, = heat flux, = flow region, = saturation, = temperature, = weighting functions, = porosity, = dynamic viscosity, = density, = non-linear coefficients in flow equation, = non-linear coefficients in heat equation.
i
= coordinate
r, J k m r t W
subscript
= nodal subscripts = iteration level = phase m = rock matrix = steam = time level = water
(A-14)
where clw and ps are given in g/(cm-set) and h, is given in J/g. 6. The values of Sw is used to calculate relative permeabilities of water and steam. For simplicity, it is assumed that
B F h h* H k
177
Indices
S
(A-12) (A-13)
+ 8.790,
NI P 4
reservoirs
References [l] F.H. Harlow and W.E. Pracht, “A theoretical
[2]
[3] [4]
[5]
[6]
in flow equations,
in heat equations,
[7]
study of geothermal energy extraction”, J. Geophys. Res., 77, 35, 7033-7048 (1972). J.W. Mercer, “Finite elements approach to the modelling of hydrothermal system”, Ph. D. Thesis, U. of Illinois, Urbana-Champaign (1973). M.L. Sorey, “Numerical modelling of liquid geothermal systems”, Ph. D. Thesis, U. of Calif., Berkeley (1975). J.W. Mercer and C.R. Faust, “Simulation of water- and vapor-dominated hydrothermal reservoirs”, paper SPE 5520 presented at 50th SPE Annual Fall Meeting, Dallas, Texas, (1975). S.K. Garg, J.W. Pritchett and D.H. Brownell, Jr., “Transport of mass and energy in porous media”, paper presented at 2nd United Nations Symposium on the Development and Use of Geothermal Resources, San Francisco, Calif. (1975). T.J. Lasseter, P.A. Witherspoon and M.J. Lippman, “The numerical simulation of heat and mass transfer in multidimensional two-phase geothermal reservoirs”, paper presented at 2nd United Symposium on the Development and Use of Geothermal Reservoirs, San Francisco, CA (1975). CR. Faust, “Numerical simulation of fluid flow and energy transport in liquid- and vapor-dominated hydro-
178
P. Huyakorn,
G.F. Pinder /Model for hydrothermal
thermal reservoirs”, Ph. D. Thesis, Pennsylvania State U., (1976). [8] R.M. Toronyi, “Two-phase two-dimensional simulation of a geothermal reservoir and well-bore system”, Ph. D. Thesis, Pennsylvania State U. (1974). [9] P.S. Huyakorn, “Solution of steady-state, convective transport equation using an upwind finite element scheme”, Applied Numerical Modelling, IPC Science and Technology Press, U.K., 1, 4, 187-195 (1977). [lo] P.S. Huyakorn and C. Taylor, “Solution of the transient convection-diffusion equation using extended upwind finite element scheme”, Applied Numerical Modelling, (in Press).
reservoirs
[l l] P.S. Huyakorn and G.F. Pinder, “A new finite element technique for the solution of two-phase flow through porous media”, Princeton U., Water Resources Program, Res. Report No. 76-WR-4 (1976). [ 121 N. Arihara, “A study of non-isothermal and two-phase flow through consolidated sandstone”, Ph. D. Thesis, Stanford U., Stanford, Calif., (1974). [13] H.J. Ramey (Jr.), W.E. Bragham, H.K. Chen, P.G. Aitkinson and N. Arihara, “Thermodynamic properties of hydrothermal systems”, Proc. of NSF Conf. on the Utilization of Volcano Energy, held in Hilo, Hawaii (1974).