Modelling of two-phase flow inside geothermal wells

Modelling of two-phase flow inside geothermal wells P. K. Chadha and M. R. Malin CHAM Limited, Bakery House, Wimbledon,

London,

UK

A. Palacio-Perez Znstituto de Zngeneria, UNAM Cd. Universitaria,

Coyoacan, Mexico

An extended homogeneous two-phase model is presented for the prediction of the two-phase flow of liquid and its vapor in a vertical well. The one-dimensional mathematical model includes some novel features for the calculation of the pressure drop along the well. The capability of the model is demonstrated by its application to three different geothermal wells for which experimental data are available and each of which presents completely different flow characteristics. In all cases the calculations compare very well with the experimental results. Keywords:

two-phase flow, geothermal wells, finite-volume method

Introduction The general problem considered is that of predicting the flow characteristics of a two-phase mixture of liquid and its vapor flowing up a vertical circular well. That is, given the conditions at one end of the well such as pressure, mass flow rate, and enthalpy, which may be determined from field measurements, the flow characteristics along the well are determined from a mathematical model that accounts for the local flow regime. The modelling of such two-phase-flow phenomena may be divided, in general, into two main categories: the so-called single-fluid models and the multifluid models (see, for example, Boure and Delhaye’). The models in the first category employ balance equations that are formally identical to the single-phase equations, the problem being to determine the physical properties of the pseudo single-phase compressible fluid. As a consequence of this approach, some parts of the solution need to be imposed, such as nonequilibria between the phases. A further subdivision of these models leads to those based on imposed velocity protiles (see, for example, Delhaye2), which are: (a) the one-dimensional, one-velocity model or homogeneous model; (b) the two-dimensional, one-velocity model or Bankoff model, which allows a velocity distribution;

Address reprint requests to Dr. Malin at CHAM Limited, Bakery House, Wimbledon, London, SW19 SAU, UK. Received 2 April 1991; revised 11 September tember 1992

236

Appl. Math. Modelling,

1992; accepted 25 Sep-

1993, Vol. 17, May

(c) the one-dimensional, two-velocity model or Wallis model, which allows a relative velocity between the phases; and (d) the two-dimensional, two-velocity model or Zuber-and-Findlay model, which is essentially a combination of the Bankoff and Wallis models. The model presented in this paper, which was developed mainly for application to geothermal wells, follows the concept of the so-called extended homogeneous model.3 The nature of this model is such that it not only takes into account the effect of slip between the phases (which is the essence of the Wallis model), but it also provides a more accurate estimate of the total pressure gradients by using a different definition of density. This density and the expressions used to calculate the slip velocity (and therefore the friction gradient) are selected according to the prevailing phase distribution of the mixture. A graphical representation of the possible flow regimes considered in this model, together with the map used for their determination, is shown in Figure 1. Although the complete model is capable of taking into account the effect of transient heat losses between the well and the surrounding rock formation (see, for example, Palacio4), the present version neglects such effects. Therefore, it considers only those cases where the well fluid and surroundings have reached thermal equilibrium. This assumption is valid, for practical purposes, when the wells have been operating for long periods under fairly constant nominal conditions. The present contribution involves the implementation of an extended homogeneous two-phase model into the general-purpose fluid dynamics program PHOENICS (see, for example, Spalding’) for the predic-

0 1993 Butterworth-Heinemann

Two-phase

flow in geothermal

wells: P. K. Chadha et al.

.. . . . .. .. . . . . . .. .. ... . . .. .. .. . . .. . . .. .. . . .. Bubble

Slug

Mist

Transition

/

‘OL 1

Continuous

I

Bubble flow

Iiqul

phas/

Alternating

phases

Coniinuous gas phase Misl

Gas velocity

flow

number,

Ngv

Figure 1. Flow regimes and map.

tion of one-dimensional vertical two-phase flows. The model includes some novel features for the calculation of the friction gradients in the slug and mist flow regimes. Specifically, one of two existing correlations is selected for the friction gradient in the slug flow regime according to the value of a new dimensionless parameter. Further, a modified version of an existing correlation is used for the mist flow regime. The capability of the model is demonstrated by its application to three practical well configurations, namely the Cerro-Prieto M-90 well in Baja California,

Mexico; the East-Mesa 6-l well at Imperial Valley, CA; and the HGPA-1 well in Hawaii. The results of these calculations are compared with the available experimental data. The remainder of the paper consists of three main sections. The mathematical model is described in the first of these sections. The second section applies the model to the three test cases and compares the results with experimental data. The third and final section of the paper summarizes the present investigation and makes some recommendations for future work.

Appl.

Math.

Modelling,

1993,

Vol. 17, May

237

Two-phase

flow in geothermal

Mathematical Conservation

wells: P. K. Chadha

and the saturated liquid enthalpy he is computed from the following correlation:

formulation equations

The steady two-phase flows considered in this paper are simulated by assuming one-dimensional flow in the well and by treating the two-phase flow as a homogeneous mixture of liquid and vapor. The model assumes thermal equilibrium between the two phases and allows for the presence of slip and several different flow regimes. The conservation equations for the mixture are given below in one-dimensional form: Continuity:

%(m=o $ (pu2) = -E - pg- TA” Energy:

;(pUh) = -q/l”

(3)

wherein p is the mixture density, U the mixture velocity, p the fluid pressure, g the gravitational acceleration, 7 the wall shear stress, q the wall heat flux, A” the wetted area per unit volume, and h the total specific enthalpy of the mixture inclusive of kinetic and gravitational contributions. In the present study, only adiabatic systems are considered and so the mixture enthalpy equation need not be solved. Further, the kinetic and gravitational contributions to the total enthalpy may be neglected in the present applications, and so the local enthalpy h is simply equal to the inlet enthalpy hi,. In order to close the foregoing equations expressions must be supplied to calculate the mixture density p and the wall shear stress T. The actual expressions used for their determination depends on the local flow regime prevailing in the well. Flow regime detection

The mathematical model must allow for the presence of several different flow regimes, as the geometrical distribution of the liquid and vapor (gas) phases strongly influences the pressure drop along the well. In addition to pure-liquid and pure-gas flow, the generally accepted flow regimes of vertical two-phase flow are: bubble flow, slug flow, transition flow, and mist flow. They are illustrated in Figure 1 and described below. The criteria for selecting between the various flow regimes are those employed by Palacio,4 as follows: Pure-liquid flow is identified when the mixture static enthalpy h is less than that of saturated liquid at the prevailing fluid pressure. Thus, pure-liquid flow is presumed to exist when: h < h,

Appl. Math. Modelling,

he = 422.1~‘.*~‘~

(5)

where p is the pressure in bar and h, is in kJ/kg. When condition (4) is not satisfied, gas is taken to be present and one of the remaining flow regimes is selected by testing whether certain parameters fall within prescribed limits, as indicated below. Flow regime

Criteria

Bubble Slug Transition Mist

NV > q&t N, < q&q1 and N,, < RN,

(7)

RN, < N,, < RN,

(8)

(6)

N,, > RN,

(9) where qg and qr are, respectively, the volumetric flow rates of gas and mixture defined by

Momentum:

238

et al.

(4)

1993, Vol. 17, May

qx = mX/p,

(10)

41 = 9n + 9e

(11)

where the subscripts t, g, and 4 denote mixture, gas, and liquid respectively; m is the mass flow rate of the mixture; qr is the volumetric flow rate of liquid qe = m(1 - XYP,

(12)

and X is the vapor mass fraction or quality, defined by x=

h-he h, - h,

(13)

The saturated vapor enthalpy h, is computed from the following correlation h, = 2692.3~O.O’~~’

(14)

where p is the pressure in bar and h, is in kJ/kg. The various dimensionless parameters appearing in (6)-(9) above, are defined by N, = (1.071 - 7.1386U2/gD, 0.13)

1

(15)

0.25

N

=%fi '"

A [ ga

(16)

RN, = 50 + 36Nev

(17)

RN m = 75 + 84N0.” el/

(18)

wherein

1

0.25

N

t?V

=%!k A

[

gu

(19)

In these relations, A is the well flow area, D is the well hydraulic diameter, u is the liquid surface tension, U (= q,lA) is the mixture velocity, and the symbol (a, b) denotes that the maximum of a and b is taken. In the foregoing, NV is an overall dimensionless velocity number defining the bubble-slug boundary, RN, is a dimensionless number defining the slug-transition boundary, RN,,, is a dimensionless number defining the transition-mist boundary, and N,, and N,, are, respectively, dimensionless gas velocity and liquid ve-

Two-phase

locity numbers. These parameters are those reported by 0rkiszewski6 based on the boundary between bubble and slug flow proposed by Griffith and Wallis’ and by Duns and Ross for the remaining three flow regimes. Mixture-density

calculation

The calculation of the mixture density p varies according to the local flow regime prevailing in the well. The formulas used in each regime for the calculation of the density are defined below. Liquid regime. In this regime the mixture density is simply set equal to the liquid density pe, which is calculated from the following expression’: p+J= loooslp

(20)

(21)

Here, the pressure p is in bar, and the temperature difference AT is given by AT= T,.- T,

(22)

where T,. and T, are the liquid reference and ambient temperatures respectively, in degrees K. Bubble regime. In the bubble regime the density of the mixture is calculated from P=P&,77K+PP(l -77x)

N,, = (Nev + Ngv)/qg wherein the volumetric computed from

rig =

(28) fraction

of gas phase Q is

qn qr + AU,

and U, is the bubble rise velocity et al.“:

given by Nicklin

In these relations, U, is the slip velocity taken as 0.2439 m/s and U.YKis the superficial gas velocity given by

(30)

The onset of the second stage of slug flow is reached when N, 2 50 and X 2 0.03; otherwise the flow is presumed to be in the first stage of slug flow where the liquid content is high. In the first stage of slug flow the mixture density is calculated from equation (23) with pu and pn computed as in bubble flow, but with nx determined from equation (29) above. In the second stage of slug flow the mixture density is calculated from

m + pcvbA

qr + W

(23)

(24)

+ TP,

(31)

where v, is the bubble rise velocity and I the liquid distribution coefficient, both of which are defined below. The densities pe and px are computed as in bubble flow. The bubble rise velocity ub appearing in equation (31) is calculated according to the prevailing bubble Reynolds number N

=

P&bD

(32)

b

U,, = q&4

(25)

Although the use of a constant slip velocity is simple, it is a reasonably precise approximation, as discussed for example by 0rkiszewski.6 In equation (23) the gas density is calculated from the following correlation px = 0.5809p”.9588

(26)

where p is the pressure in bar. The liquid density is calculated from equation (20) with T, taken to be the saturation temperature T,, which is calculated from the following formula: T, = 255.37 + 117.5p”.225

(27)

where p is the pressure in bar and T, is in degrees K. Slug regime. In the slug flow regime, two different sets of density formulas are used to reflect the fact that larger pressure drops are observed at the onset of the regime, when the liquid content is high, than later on in the regime when the liquid content is relatively low, say when X > 0.1. This practice of dividing the slug

et al.

flow regime into two distinct stages is a new feature of the present model. Following Palacio,4v’o the boundary between the two different slug regime stages is determined through the use of the following dimensionless parameter

P=

where 71~is the gas volume or void fraction, given by the expression proposed by Griffith and Wallis’ 7& = OS[{l + U/U,} - [{I + U/U,]2 - 4U,,/U,]“*]

wells: P. K. Chadha

u, = 1.2u + 0.35(gD)O.5

where S is the specific gravity and /3 = 1 + 2.148. 10-4AT + 3.204.10-6AT2 - 4.83.10p5p

flow in geothermal

I-Q

where Ub is computed from equation (30) and pe is the viscosity of the liquid, which is determined from the following correlation: j_+ = 2.185.10-3/(0.28921

+ 0.07281T + 1.67.10p5T2)

(33)

where T is the temperature in degrees C. Once the bubble Reynolds number is known, v, is calculated from the following set of equations: When Nb 5 3000 v, = [0.546 + 8.74.10e6Re] (gD)o.5

(34)

When Nb 2 8000 v, = [0.35 + 8.74.10m6Re] (gD)o.5

(35)

When 3000 < Nb < 8000 v, = 0.5u,i + [I& + 1.1 1656.104pe/(pe~D)]o~5 (36) where vbi = [0.251 + 8.74.10p6Re] (gD)o.5

Appt. Math. Modelling,

1993, Vol. 17, May

(37)

239

Two-phase

flow in geothermal

In the foregoing set of equations, number defined by Re

wells:

P. K. Chadha et al.

Re is the Reynolds

(38)

The liquid distribution coefficient r is determined from the following equations: When U < 3.048 r = -0.2132546U

- 0.1

(39)

When U 2 3.048 m/s r = (r,, r,)

(40)

where rl = [1.7421568.10~210g,o(103~,)]/Do~799 - 0.709 - 0.162 log,,(U) - 0.888 log,,(D) - 0.541784 (41) and

r2 =

vbA 9r + ‘&A

(1

-

P/P,)

(42)

The foregoing formulas for v, and r originate from the work of Orkiszewski.6 Mist regime. The mixture density in the mist regime is computed from equation (23) as for the bubble regime, except that the gas-phase volume fraction is calculated from (43)

77s= 4&A

rather than from equation (24). Transition regime. In the case of the transition flow, the mixture density is calculated by linear interpolation of values calculated for the slug and mist regimes, as follows: P = (PA

RN,,, - Ngv RN,,, - RNs

+ (P),

- RN, gv RN,,, - RN,

(48)

where k is the roughness height. In the pure liquid regime, the liquid velocity U, is equal to the mixture velocity U, whereas in bubble flow it is calculated from (49)

Ue = qA[A(l - ~~11

where qa is given by equation (24) and the quality X is calculated from equation ( 13). Slug regime. Following the practices used earlier for the density calculation, two different sets of shear stress formulas are used in the slug flow regime. These formulas differ according to the local values of the quality X and the dimensionless parameter N,, both of which are used to determine whether the two-phase mixture is in the first or second stage of slug flow (see the discussion of slug regime in the section on mixturedensity calculation). For the first stage of slug flow the shear stress is calculated from the correlation of Chierici et a1.,13 which is given by T = (1 - qK)pefU2/8

(50)

where vr is given by equation (29) and f is determined from equation (46) wherein the Reynolds number is defined by equation (47) but is based on the mixture velocity U. For the second stage of slug flow the correlation of 0rkiszewski6 is used to calculate the wall shear stress, as follows:

1

The calculation of the wall shear stress T, as with the mixture density, varies according to the local flow regime prevailing in the well. The formulas used in each regime for the calculation of the shear stress are defined below. Liquid and bubble regimes. In the pure-liquid and bubble regimes, the shear stress is calculated as follows: 7 = pef U’,lS

(45) where f is the friction factor calculated from Swamee and Jain’s” formula:

ln,~~L_J

Appl. Math. Modelling,

E = klD

(44)

Shear stress calculation

240

and E is the relative roughness

N

where (P)~ and (p), are, respectively, the values of mixture density evaluated using the formulas appropriate to the slug flow and mist flow regimes.

f =

number Re is defined by (47)

peuD

=

and here the Reynolds

(46)

1993, Vol. 17, May

(51)

where f is evaluated from equation (46) using the mixture velocity, and vb and r are calculated according to equations (34)-(42). Mist regime. In the mist flow regime the shear stress is calculated according to the correlation of Tachimori,14 as employed by Palacio4,” in modified form: T=X*T,

+ (1

-x)Te

(52)

where Tn = fq;p,U;/8

(53)

Te = f(l - 71,)~peu?3

(54)

The modification involves replacing Tachimori’s shifting factor y with the mixture quality X, and this practice avoids the need to adjust y from one well to another. This modified form of Tachimori’s correlation is another new feature of the present model. In equations (52)-(54), Q. is calculated from equation (43), the gas and liquid velocities from equations (25) and (49), respectively, and the calculation of the friction factor

Two-phase

proceeds according to the local value of a relative roughness factor E*, as follows: For E. > 0.05: f = 0.25[log,, {0.27~}] -* + 0.268~!.‘~

(55)

For E.: 5 0.05, f is computed from equation (46) with the Reynolds number Re defined by Re

=

pk’“k?D

(56)

and the relative roughness E is taken as the maximum of E, as defined by equation (48), and E., which is calculated according to the local Weber number We, as follows: For We 5 0.005: (57)

For We > 0.005: E. = 174.8a We0.302/(pg UZD)

(58)

where (T is the surface tension, and the Weber number is defined by We = p,U~klu

(59)

and k is the roughness height. Transition regime. For the transition regime, the shear stress is calculated by linear interpolation of values calculated for the slug and mist regimes, as follows: RN,

- N,,

RN,,, - RN,

+ (T),

N,, - RN, RN,,, - RN,

(60)

where (T)~ and (T), are, respectively, the values of the wall shear stress evaluated using the formulas appropriate to the slug flow and mist flow regimes. Solution of the equations

The complete homogeneous two-phase model is introduced by utilizing the facilities provided by PHOENICS for the attachment of user-generated coding sequences. The one-dimensional differential equations for continuity and momentum are solved numerically using the parabolic finite-volume procedure embodied in the PHOENICS computer program. The parabolic option obviates the need to specify a downstream boundary condition, as the solution is obtained as a marching integration in the x-direction, starting from prescribed values of all dependent variables at the bottom of the well. At each forward step the continuity and momentum equations are solved by an iterative procedure based on the well-known Patankar-Spalding’s method. In the context of the present work, the essence of this method is that at each forward step the downstream pressure p. is determined by first calculating the mixture velocity U from the momentum equation with an estimated downstream pressure and then obtaining an appropriate correction so as to satisfy the continuity equation. This process is then repeated by iteration until convergence of the equation set.

wells: P. K. Chadha et al.

For cases in which the inertia and friction terms are negligible compared with the gravitational and pressure terms a momentum-based pressure adjustment was introduced so as to accelerate convergence at each forward step. This adjustment proved essential only for solving the second test case, namely the East-Mesa 6-l well at Imperial Valley, CA. It should be mentioned that no under-relaxation practices are necessary for any of the cases considered in the present study. In general, the number of iterations per forward step is set equal to six, although typically it is found that only three or four iterations are required for convergence of the equation set at each step. Model application

E* = 34ul(p, up,

r = (7)s

flow in geothermal

Cerro-Prieto

and results

well

The first test case is the Cerro-Prieto M-90 well in Baja California, Mexico. ” The well-bore diameter is 0.18288 m, the depth is 1298.7 m, and the roughness height is 91.44 pm. For this case the fluid enters the base of the well as saturated liquid with the following inflow conditions: Mass flow rate Pressure Temperature Density Quality Enthalpy Velocity

m = 45 kg/s pin = 88.48 bar T,,, = 577.53 K (304.38 “C) pin = 770 kg/m3 xi, = 0.0 hi, = 1.338 MJlkg Ui, = 2.225 m/s

The ambient temperature T, is taken as 288.7 K, the specific gravity S as 1.02, and the surface tension (T as 0.0729 N/m. To examine the effects of grid size on the solutions, calculations are made with four different forward-step sizes, as follows: AX = 21.65, 10.82,7.22, and 5.41 m. It is found that grid independence is achieved with a forward-step size AX = 10.82 m. Figure 2 compares the measured and predicted pres-

100

60

II -

Pressure Bar) 60

-

I /I I

I 000

I 400

Depth

Figure 2. measured

Appl.

I

II

1200 (ml

Cerro Prieto: comparison between calculated and pressure profiles along the well.

Math.

Modelling,

1993,

Vol.

17, May

241

Two-phase flow in geothermal

wells: P. K. Chadha et al.

sure profiles for the well, and the experimental data are taken from Ortiz.i6 Figures 3-5 show, respectively, computed profiles of the quality, densities, and mixture velocity. Each figure includes demarcation lines indicating the onset of different flow regimes along the wall. The density plot shows results for the gas and liquid phases, as well as for the mixture. From Figure 2 it can be seen that the predicted variation of pressure with depth is in good accord with the measurements. Further, according to the flow regime criteria utilized in the present model, the mixture evolves from saturated liquid conditions through the

0.15

0.10

Ouality 0.05

0

Figure 3. well.

I 0

Bubble

--

I I

400

800

1200

Depth(m)

Cerro Prieto: calculated evolution of quality along the

IOOOLiquid 800

Density (Kg/m’) 600 400

Bubble

-

I

-

800

12cm

Depth

Figure 4. Cerro Prieto: calculated variation and mixture density along the well.

( m)

of the gas, liquid,

I! Bubble

--

IC

U,(mh)

( Depth(m)

Figure 5. Cerro Prieto: calculated locity along the well.

242

Appl.

Math. Modelling,

variation

of the mixture ve-

1993, Vol. 17, May

bubble flow regime and finally into both types of slug flow regime. The continued reduction in fluid pressure along the well due to friction results in the formation of vapor and a progressive increase in its quality, as may be seen from Figure 3. The calculated well-head quality is 13.85%, which is within 1% of the measured quality reported by Ortiz.i6 The reduction in mixture density along the well is evident from Figure 4, and this expansion is accompanied by a corresponding increase in mixture velocity, as shown in Figure 5. The well-head mixture temperature, density, and velocity are predicted to be 527 K, 131.6 kg/m3, and 13.1 m/s, respectively. It should be mentioned that abrupt change in the mixture density variation shown in Figure 4 at points of regime transition (mainly from slug 1 to slug 2) is due to the different nature of the density equations used, which represent a parameter used to account for liquid holdup. East-Mesa well The case is the East-Mesa 6-l well at Imperial Valley, CA.17 The well-bore diameter is 0.2215 m, the roughness height is 91.44 pm, and the depth is 2133.6 m. For this case subcooled liquid enters at the base of the well with the following inflow conditions: Mass flow rate Pressure Temperature Density Quality Enthalpy Velocity

m = pin = T;,, = pin = Xi, =

12.925 kg/s 92.9655 bar 471.48 K (198.33 “C) 893 kg/m3 0.0 h, = 855.26 k.I/kg U,, = 0.3756 m/s

The inlet enthalpy corresponds to that given in “Steam Tables” for the specified inlet pressure and temperature, and the degree of subcooling is 109.66 “C. Again, the ambient temperature T, is taken as 288.7 K, the specific gravity S as 1.02, and the surface tension r as 0.0729 N/m. A grid-independent solution is achieved with a forward step size of Ax fi 8.5 m. The results are presented in Figures 6-9 in terms of longitudinal profiles of pressure, quality, mixture and phase densities, and mixture velocity. The experimental pressure data are taken from Gould,17 and the mixture density plot includes profiles of the liquid and vapor density. All of the figures include demarcation lines identifying the different flow regimes that the model predicts to prevail within the well. The predicted pressure profile shows very good agreement with the data given in Figure 6. It can be seen that the model predicts that the flashing point is located some 1248 m below the well head, above which the mixture evolves through the bubble, slug, transition and mist flow regimes. This result agrees well with the expectations from field measurements reported by Gould. I7 The formation of vapor at the flashing point is clearly evident in Figure 7, which shows the evolution of

Two-phase

flow in geothermal

wells: P. K. Chadha et al.

dominated flow. The well operates under conditions of mist flow throughout the entire length of the well. For this particular case, the well has two different diameters, as follows: 0.22555 m to a depth of 762 m from the surface and 0.17617 m from that point to the bottom of the well located at 1966 m. The roughness height of the well bore is 91.44 pm. For this case wet steam enters the bottom of the well with the following inflow conditions:

I

I 1500

2 000

Depth(m)

Figure 6.

East Mesa: comparison between sured pressure profiles along the well.

calculated

and mea-

m = 12.47 kg/s pin = 19.507 bar T;, = 484.65 K (211.5 “C) = 15.3 kg/m3 2 ~064 h;; = 2:1194 MJ/kg Ui, = 33.44 m/s

Mass flow rate Pressure Temperature Density Quality Enthalpy Velocity

As in previous test cases, the ambient temperature T, is taken as 288.7 K, the specific gravity S as 1.02,

and the surface tension u as 0.0729 N/m. For this particular well, grid independence is achieved with a forward step size of Ax 5 8.2 m. As in previous test cases, the results are presented in terms of profile plots of the variation in flow properties with depth. Figures IO-13 present, respectively, profiles of the static pressure, quality, densities, and mix-

0.10 -

--Mist

Quollty Transalmn 0.05

-

Oorn Depth (m)

Figure 7.

East Mesa:

calculated

evolution

of quality

40

along the

well.

30 U&m/s) 1000 20 a00

Densily

Ikglm’) 600

Depth

Cm)

200

Figure 9.

East Mesa: calculated ity along the well.

0 0

500

IO00

1500 Depth

Figure 8. mixture

East Mesa: calculated density along the well.

variation

of the mixture

veloc-

2 000 (m)

of the gas, liquid, and

quality along the well. The model predicts a well-head quality of 14.67% and an efflux temperature of 400 K. As may be seen from Figures 8 and 9, the onset of vapor formation is marked by a rapid expansion of the fluid during which the mixture density reduces markedly and the velocity increases. At the well head, the model predicts a mixture density of 9.44 kg/m3 and a mixture velocity of 35.53 m/s. HGPA-1

variation

Diomster

01

0

=0.2255

I 500

I I

well

The third and final test case is the HGPA-1 well in Hawaii,14 which provides a test of the model for vapor-

I

I 1000

1500 Depth

Figure 10. sured

HGPA-1: comparison between pressure profiles along the well.

Appl.

Math.

Modelling,

1993,

21 00

(m)

calculated

and mea-

Vol. 17, May

243

Two-phase

flow in geothermal

wells: P. K. Chadha et al.

reveals the sudden decrease in mixture velocity at the diameter-expansion location.

Concluding 0.66

Oiam~ter=O.2255 500

0

1000

2000

1500

Depth(m)

Figure 11. well.

HGPA-1:

calculated

evolution

of quality along the

201 I5 -

Density

( kg/m3 ) IO -

0

I 500

0

I 1000

Oiomclsr:0.17617 I 1500

2 000

Depth(m)

Figure 12. HGPA-1: calculated variation of the gas and mixture density along the well.

40

O~ometrr=0.2255

20 0

Owmrtrr

I

500

:0.17617

I IO00

I 500

2000

Depth (m)

Figure 13. HGPA-1: calculated variation of the mixture velocity along the well.

ture velocity. The pressure measurements are taken from Tachimori,14 and the density plot includes profiles of the mixture and vapor density only. It can be seen on inspection of Figure 10 that the calculated pressure profile is in excellent agreement with the experimental data. Figure 11 reveals that the steam quality increases throughout the well, exhibiting an abrupt reduction in its rate of increase at the location of the diameter change. The predicted quality and temperature at the well head are 70.4% and 422 K, respectively, which compare very well with the measured well-head quality of 70.82% reported by Tachimori. l4 The fluid expansion throughout the well results in a computed well-head mixture density of 3.656 kg/m3 and mixture velocity of 85.37 m/s, as may be noted from Figures 12 and 23, respectively. The latter figure

244

Appl.

Math. Modelling,

1993, Vol. 17, May

remarks

A one-dimensional mathematical model for the prediction of a two-phase mixture flowing up a vertical well has been implemented into the PHOENICS program. The capability of the model has been demonstrated by its successful application to three geothermal wells, each of which exhibits completely different flow characteristics. These cases were selected so as to validate the model over the whole range of flow regimes considered by the model. The good agreement found between the experimental and calculated pressure profiles attests to the capability of the model, which includes some new features for the calculation of the friction gradient in the slug flow and mist flow regimes. For the case where the fluid enters the bottom of the well as subcooled liquid (the East-Mesa 6-l well), it has been shown that the model gives a very good estimation of the flashing point location within the well. To increase the range of applicability of the present model, it is recommended that heat losses between the well and the surrounding formation be included, as under certain conditions these can have a significant effect on the properties of the mixture. This might be done to a first approximation by calculating the wall heat transfer rate on the basis of the ReynoldsColburn equation (see, for example, 0zisik18). A second option would be to follow Ramey’s approach” for the calculation of transient heat losses. This approach neglects rapid changes in the fluid temperature and also other effects such as viscous heating of the fluid and longitudinal heat conduction within the surrounding rock formation. Another option is to consider a constant overall heat transfer coefficient which, when multiplied by the product of the fluid-rock temperature difference and the corresponding cross-sectional area, represents the source or sink of heat to be included in the energy balance equation. This term may take into account the transient, two-dimensional temperature distribution within the rock formation or it may simply assume a linear depth-dependent distribution under steady-state conditions.

References Boure, J. A. and Delhaye, J. M. General equations and twophase flow modelling. Handbook of Multiphase Sysfems, ed. G. Hestroni, Hemisphere Publishing Co., 1982, pp. 1.36-1.87 Delhaye, J. M. Basic equations for two-phase flow modelling. Two-Phase Industries,

Flow and Heat

Transfer

in the Power and Process

ed. A. E. Bergles et al., Hemisphere Publishing Co., 1981, pp. 40-97 De Gance, A. and Atherton, R. Vertical and Inclined-Flow Correlations. Chemical Engineering Aspects of Two-Phase Flow, 1970, Part 6, pp. 87-94 Palacio A. Dinamica de Pozos Geotermicos. Ph.D. Thesis, National University of Mexico, 1987 Spalding, D. B. A general-purpose computer program for

Two-phase

6 7 8

9 IO

multi-dimensional one- and two-phase flow. Maths. Comput. Simul. 1981, 13, 267-276 Orkiszewski, J. Predicting two-phase pressure drops in vertical pipes. .I. Per. Tech. 1967, 19, 829-838 Griffith, P. and Wallis, G. Two-phase slug flow. _I. Hear

14

Transfer

15

Geothermal

II 12

1961, 307-320

Duns, H. Jr. and Ros, N. C. J. Vertical Flow of Gas and Liquid Mixtures in Wells. Proceedings of the 6th World Petroleum Congress, Frankfurt, Section II, Paper 22.PDG, 1967, pp. 451-465 Frick. T. Petroleum Production Handbook. SPE, Dallas, 1962 Palacio, A. Effect of heat transfer on the performance of geothermal wells. Paper accepted for publication in Inr. J. Research

and Its Applications.

13

16 17

1990. 19, 31 l-328

Nicklin, D., Wilkes. J., and Davidson, J. Two-phase flow in vertical tubes. Trans. Inst. Chum. Eng. 1962. 40, 61-68 Swamee, P. and Jain, A. Explicit equations for pipe-flow problems. J. Hydr. Div. Proc. ASCE, 1976, 657-664

18

19

flow

in geothermal

wells:

P. K. Chadha et al.

Chierici, L., Giannone, G., and Slocchi, G. A wellbore model for flow in geothermal reservoirs. SPE 10315, 1985 Tachimori, M. A numerical simulation model for vertical flow in geothermal wells. Proceedings of the 8th Workshop on Geothermal Reservoir Engineering 1982, pp. 155-160 Patankar, S. V. and Spalding, D. B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int. J. Hear Mass Transfer, 1972, 15, 17871806 Ortiz, R. J. Two-phase flow in geothermal wells: Development and use of a computer program. Report SGP-TR-66, Stanford Geothermal Program, Stanford University, 1983 Gould, T. J. Vertical two-phase steam-water flow in geothermal wells. J. Pet. Tech. 1974, 26, 833-842 Ozisik, M. N. Heat Transfer: A Basic Approach, Mechanical Engineering Series, McGraw-Hill Book Co., New York, 1989 Ramey, H. J. Wellbore heat transmission. J. Pet. Tech. 1962. pp. 427-435

Appt.

Math. Modelling,

1993, Vol. 17, May

245