Two-phase flow dynamics by a five-equation drift-flux model

0094-4548/81/010057-12502.00/0 ©PergBmc~Press Ltd. Printed in theUnitedStates

IN H F A T A N D M A S S TRANSFER

Vol. 8, pp. 57-68, 1981

TWO-PHASE FLOWDYNAMICSBY A FIVE-EQUATION DRIFT-FLUXMODEL

C. Kim and R. P. Roy Nuclear Engineering Program University of I l l i n o i s Urbana, IL 61820 (OmTrmmicated by J.P. Hartnett and W.J. Minkowycz)

Introduction Analyses of transient two-phase (specifically, vapor-liquid) flow problems have been carried out by various researchers employing models ranging in complexity and detail from homogeneous, thermal equilibrium (three conservation equations) to two-fluid unequal velocity, unequal temperature (six conservation equations) [1-3].

The five-equation

d r i f t - f l u x model is one of intermediate complexity and detail.

The intent

of this work is to extend the application of a numerical solution technique termed the Newton Block Gauss Seidel (NBGS) scheme, introduced in [4] and then applied to a four-equation d r i f t - f l u x model, to a five-equation d r i f t - f l u x modelt. While the addition of the f i f t h equation (an additional phasic thermal energy equation) adds to the generality of the model by removing a thermal constraint, the computations become more elaborate. We present numerical solution results for two transient problems adopted from [4] and [5-modified] respectively.

t

The version used here is somewhat different from [3] in that in place of the mixture mass and thermal energy conservation equations, the corresponding liquid phase conservation equations are included. This choice appears to yield somewhat better computational results for some problems.

57

58

C. K i m a n d R.P. Roy

Vol.

8, No. 1

Ana~ysi~

The.. .model e.quations The time-dependent, one-dimensional conservation equations are: Liquid mass:

@ a a [~(I-~) PgP~Vrl a-t- [(l-~)p~] + a-i [(I-~)Pj~Vm] - Ti L " %1 J [

Vapor mass: Mixture , momentum :

Vapor thermal energy:

3-t [(~Pg] + ~

~

'

'

= -r ( I )

I

:

£ (2)

_ l Pm 8z ~p+ gz - K VmlVml

(3)

Pm

L

J

BVm aV l @__ [ ~(l -~) Pgp~Vr2] at + Vm@-~ + - Pm Pm az L

a

~ [~pgeg] + Tz [~pgVmeg] + ~z

+ p~

[~Vm] + P

"~(l-~)pgp~Vregq Pm "~(l-~)P~Vrl

J @~

qwg + qig + £ hsg

(4)

Liquid thermal a-t [(l-~)p~e£] + -~ [(l-~)P~Vme~] - a .~(l-~)PgP~Vre~. energy: L °m ]

+P

[(l'~)Vm] - P ~

Pm

]

a(.l-~)..

qw~ + qi~ - £ hs~ The equations Liquid thermal: Ofare:State used Vapor thermal : Liquid caloric: Vapor caloric :

p£ : p~ pg = pg e£ = e£ eg = eg

(p,T~) (p, Tg) (p,T£) (p,Tg)

(5)

(6)

The net evaporation rate per unit mixture volume, £, is: £ = (-qig-qil) (7) (hsg - hs~) After subtraction of (Vm.mixture mass conservation equation)

Vol. 8, No. 1

T~D-PHASE ~

DYNAMICS

59

The following expressions have been adopted for qtg and qtl [4,6]: qtg =

~cXc=(l-a) (hs9 - hs~) Tsl/2 (Ts-Tg)

(8)

qtl =

~eXe~(1-~) (hsg-hs~) Tsl/2 (Ts-T~)

(9)

where Cc = A pg(Rg) 1/2 and eC = A p~(Rg) 1/2.

In

Also:

Xe = O.l, T~ ~ Ts

Xc = 0.1, Tg ~ T s

0.0, T~ < Ts

0.0, Tg > Ts

A = [~ NO a] 2/3

,~<_0.5

[~ No(1-~)]2/3,~>0.5

these interfacial (between the phases)heat transfer expressions, the

area is estimated on the basis of vapor bubble/liquid droplet configurations of the distributed phase. This is merely for calculational simplicity and can be easily replaced by a flow-regime-dependent (bubbly, slug, annular, etc.) interfacial area model. Appropriate flow regime maps [3,7,8] have been included in our computational model for proper prescription of relative velocity, Vr, between the phases. Equations ( I ) - (5) have also been used here for representing pure vapor ( ~ = 1.0) and pure l i q u i d (~= 0.0) regions by a r t i f i c i a l introduction of small amounts of l i q u i d and vapor respectively.

I t is adequate (for

computations) for example to render ~= 0.9999 in case of pure vapor and~ = 0.0001 in case of pure l i q u i d . An alternative would be to use single phase conservation equations when required. Solution technique A semi-implicit in time, spatially staggered (Eulerian grid) f i n i t e difference scheme [4,9] in which Pk' p' Tk' and ~ are defined at the computational mesh cell centers and Vm is defined at the cell boundaries is used. Pk' ~ ' and Tk are prescribed at the cell boundaries by donor c e i l i n g , as is Vm at the cell centers. However, central differencing is generally used by us for prescribing Pm at the cell boundaries.~ ~This approach appears to yield better computational results for some problems (eg. the o s c i l l a t i n g manometer problem described l a t e r ) .

60

C. Kim and R.P. Roy

Vol. 8, No. 1

The NBGStechnique is then applied to the difference equations (corresponding to equations (1) - (5)) and the constitutive equations• This technique has been described in detail in [4] and only the main points will be noted here. The f i r s t step is to linearize all equations around the latest iterate values of the unknowns. This is followed by the elimination of all unknown state variables (at the new time level) except the primary variables ~, P, Vm, T~, and Tg by using the linearized equations. Finally, Vm at the new time level is eliminated by using the mixture momentum equation• We now have four linear algebraic equations for four unknown primary variables ~, p, T~, and Tg, corresponding to each mesh cell. The complete equation system, with block tridiagonal coefficient matrix is shown in Fig. 1. All

Xl

Al2

Bi

X2

Ai,i

Ai , i - I

Ai,i+l

Xi -l Xi

Bi (lO)

Xi+l

AN,N-I

AN,N

XN_I XN

BN

Fig. 1 The linearized equation system (N = no. of mesh cells) Here, Xi represents the ith-cell primary state variable column vector [~i Pi Tg,i T~,i]T. Ai, i represents a 4x4 matrix whose elements are as follows (the superscripts k and n indicate the present time step previous iteration value and the previous timestep value respectively): all : -

k (P~)i

Vol. 8, No. 1

%I~D-I:~.ASEFLOWDYNAMICS

I[ ]" ]

(l-c~)p~ Pm

n

1

61

+

(l+KnlVm~+i/2 I)

i

kn

k

i-I/2 (l+Knlvn I) Jmi-I/2

t~Tt~ cr 5 ]/2 si

k a13 =

(T' k)ll 2

:

al4

p@ -k

~---nxkn,~

\~T~] i

(Is)) 1/2

k~ be ec~itm-c~i y l

a21 : (pg)~ a22 : c~i ~@P] i + A-z-z~(~Pm/i+I/2 (l+KnIVn l) \ Pmli-ll2 mi+I/2

¢ k,~.~C.,_on. ~]p"~7 ,,t c I

\B---p-Ii (Ts~)ll2

At ~-nkk CcXpi(l-~ni) a23 : o~ \@Tg/.i - (Ts~)I/2

/.~?

At ~-nxk n(l_c k) ee i a24 : - (T k)I/2 si

k e k+ k a31 : (Pg)i ( g)i Pi a32 = c~i

~@P/i

" ~@Pl i 3

L\ Pm

I~+ii~ Pmi+I12]

l

(l+KnlVmn

i+I12

[)

LX Pm

Ii-I12 Pi t~mm)i_ll2J(l+Kn I vni-I/2 I)

62

C. Kim and R.P. Roy

=n k k _n) CcXcmi (I - At

(Ts~1)l,/2

(hsg -

Vol. 8, No. 1

~;)P ii

+At [Ce~eC~,(l-c~,,+ ~ccXcC~i(1-(~,,]\ T J i (Ts~,i/Z (esg+

i

kk n k.Be .W. L~cXc~i(l-~i) k hs~ At

k/@P\k

~,,. °~ [(e,',i~)~+ (,,',t~l,J÷ (, ~,,,, (,~....

i

c c i t - i ) (esg k _

(Tsk)I/2 i k =n kn k (e k+P~--)i At CeXe~i (l -~i ) sg Pg a34 =-

(Ts.~1/2 i

a41 = - (p£e~)~-p~ (~,~k

.At 2 t~l-~)PEe,~ n

a42 = (l-a~) (p¢)~xB-Tli * (~)i (L\ Pm

I i*I/2

,

(]+KnlVmi+i/~l ÷ - ~m -)i-I/2

- At

~n. k n(l _~i~ e*e~i

+ p~ (l-~. Pmi+l /2]

k1-.n ]

(BTs~k

si 43

=

I, CACai (1-(~)

es~k

'

Pi(~--)i-I/2 (I+~IVmi_i/2 n I}

+

t

1 1/2 ' (Ts~)

i

1

Vol. 8, No. 1

a44 (l-c~) =

+ At

qWO-PHASE FLOW DYNAMICS

(Bp~k 1

(e~) k (Bp~_~k

i ~@T~)i + (PJL)~ ~T~)I

Ce~ek~i

(l-c~k)

es +

63

_-n~k n,.

I;e eC=iU-c¢ik,J

+ At (Ts~)I/2.

k

k

(hsg - hs~)

(Ts~)I/2

Elements of the 4x4 matrices Ai,i+ l and Ai,i_ 1 w i l l not be given here due to lack of space. Similarly, elements of the 4xl column vector Bi w i l l not be given here but are easily derived. I t is of interest to note that the Ai, i matrix becomes singular when = is equal to either 0.0 or 1.0. We have mentioned earlier however that in both cases, a small amount of the other phase is a r t i f i c a l l y injected for computational purposes. The unknown primary state variables Xi (i=I . . . . ,N) can now be solved either a cell-by-ceil iteration scheme which sweeps across the flow system beginning at an appropriate end, or, by inversion of the entire NxN coefficient matrix.

Both methods were used by us to solve the two transient

flow problems discussed next. The Problems 1.

The Oscillating Manometer [4]

O.Im

I0

;= EQUILIBRIUM

LEVEL

/ / i

/8/ r .f--s~

I

I VAPOR PHASE LIQUID PHASE

Fig. 2

I n i t i a l liquid and vapor phase distribution in the manometer

64

C. K i m a n d R.P. Roy

Vol. 8, No. 1

Fig. 2 schematically depicts a U-tube manometer, its subdivision into ten computational cells, and the prescribed i n i t i a l distribution of liquid and vapor (water and steam respectively, although air can be easily substituted for steam). Given that the i n i t i a l velocities in all cells are zero, we are to calculate the variation of vapor and liquid velocities and cell vapor volume fractions as functions of time.

The i n i t i a l

calculated based on hydrostatic considerations. and wall heat fluxes qwk (k = ~, g) are zero.

pressure distribution is

The evaporation rate, £, Both zero and nonzero

wall f r i c t i o n (K = O, K > O) cases are studied.

A

oO

v

E

>-

0.4

ANALYTICAL SOLN. / ~ ( Z ~ t = l ms)/&z~aAA I

I

I

I

a

02

0

1

FIVE-EQN.

I

MODEL

__1 b.I 0

0.0

I

D --~ 0.0 N

-0.2

-0.4

i

l

i

I

I

i

0.2

0.4

0.6

0.8

,.0

,.2

t

L4

(s)

Fig. 3 Liquid axial (z) velocity for the oscillating manometer problem Fig. 3 shows the variation of liquid velocity with time for the case of zero wall f r i c t i o n and Vr set equal to two-tenths of the suggested relative velocities for appropriate flow regimes [3].

Both numerically

obtained and analytically derived (liquid inventory equation of motion solution) results are shown. Similar numerical results were obtained for Vr equal to one-tenth of the suggested values as well as zero.

However,

attempts to increase Vr beyond the two-tenths magnitude resulted in numerical i n s t a b i l i t y .

Somenumerical damping of the oscillation amplitude

should be noted in Fig. 3.

Fig. 4 shows the vapor volume fraction in some

cells at 1.5s into the transient.

The numerical results shown compare

favorably with the four-equation d r i f t - f l u x model solution [4].

Vol. 8, No. 1

T~D-PHASE FI.OW DYNAMICS

A •

2

d Z

•

3

0

t=l.5

S

& F I V E - E Q N . MODEL CALCS. • I

I

ANALYTICAL I

I

SOLN. I

0.0 0.2 0.4 0.6 0.8 1.0 VAPOR VOLUME FRACTION

Fig. 4

2.

Vapor volume fraction in various cells at 1.5 s (Vr = 0.2 , recommended steadystate values)

Instantaneous Heat Addition to Liquid Flow on a Vertical Channel [5] A vertical channel of 12.72 mm inside diameters through which,

initially,

subcooled water is flowing upward with a velocity of 1.089 m/s

and entropy of 1184.7 J/g is considered.

The channel length is subdivided

into three sections - (sections 1 to 3) of length 2.69 m, 3.281 m, and 1.345 m respectively. respectively.

The inlet and outlet pressures are 7.102 MPa and 6.984 MPa To begin the transient, a uniform heat flux of 9000 W/m is

applied instantaneously to the f l u i d in Section 2.

We are to calculate the

evolution in time of the new steadystate channel condition. Fig. 5 shows the calculated time variation of the f l u i d total mass flow rates at the channel inlet and outlet.

Results obtained in [5] by a

characteristic f i n i t e difference solution of a similar problems compare well qualitatively.

Fig. 6 shows the vapor fraction axial profile at the

new steadystate condition. §The diameter adopted by us is larger than the value mentioned in [5]. There appears to be a discrepancy between the channel diameter, i n i t i a l liquid velocity, and total i n i t i a l liquid mass flow rate in [5].

65

66

C. Kim and R.P. Roy

I

I

0.4

I

Vol. 8, No. 1

I

I

OUTLET

~ 0.3 <= )0.2 O

INLET

.J I.L

CO u) O.I < ~E 0.0 0.(

At = 0.5 I

I

1.0

2.0

I

I

ms I

3.0 4.0 TIME (s)

5.0

6.0

Fig. 5 Variation of total mass flow rates at inlet and outlet with time z

0 1

I

I

I

I

I

I 2

I 3

I 4

5

6

I

(.9 .18 rr LI. w .12 1 0> . 0 6 nO

a_ O0 0 >

7 7.316

(m)

Fig. 6 Vapor fraction distribution at new steadystate

Discussion A five-equation d r i f t - f l u x model in conjunction with the Newton Block Gauss Seidel method have been used here to analyze two transient two-phase flow problems. The work reported here is an extension of the computational technique reported in [4] to the five-equation formulation.

Stable

Vol. 8, No. 1

TAO-PHASE FI_OW DYNAMICS

67

numerical results are obtained provided, in some cases, the magnitude of the relative velocity between the phases is below a certain l i m i t .

It may be

possible however to suppress the numerical i n s t a b i l i t i e s , when they arise, by a r t i f i c i a l numerical means (for example, by increasing damping).

Use of

steadystate correlations for relative velocity in lieu of a second phasic axial momentum equation has been known to contribute to solution i n s t a b i l i t y in certain cases. Nomenclature ek

k-phase internal energy per unit mass (k=g,~)

hsk

k-phase saturation enthalpy

K

wall f r i c t i o n coefficient

NO p

no. of vapor bubbles per unit mixture volume static pressure

qik

interracial heat flow into k-phase

qwk Rg

wall heat flow into k-phase gas constant

t

time

At

timestep in the numerical scheme

Tk

k-phase temperature

Vk

k-phase axial velocity

Vm

mixture axial velocity, = [~ pgVg + (1-~)p~V~]/pm

Vr

relative velocity between the phases, = (Vg-V~)

z

axial coordinate

F

evaporation rate per unit mixture volume

Pk

k-phase density

Pm

mixture density, = [~ pg + (1-~)p£]

vapor fraction

References 1.

M. I s h i i , Thermo-Fluid Dynamic Theory of Two-Phase Flow, Eyrolles, Paris (1975).

2.

W. T. Hancox, R. L. Ferch, W. S. Liu, and R. E. Nieman, Int. Journal Multiplase Flow, 6, 25-40 (1980).

3.

TRAC-PIA, An Advanced Best-Estimate Computer Program for PWR LOCA Analysis, NUREG/CR-0665 (LA-7777-MS) (1979).

4.

D. R. Liles and W. H. Reed, Journal of Computational Physics, 26, 390-407 (1978).

68

C. K i m a n d

R.P. Roy

Vol. 8, No. 1

5.

W. T. Hancox and S. Banerjee, Nuclear Science and Engineering, 64, 106-123 (1977).

6.

W. G. Rivard and M. D. Torrey, Los Alamos S c i e n t i f i c Lab. Report LA-6104-MS (1975).

7.

C. W. Solbrig, J. H, McFadden, R. W. Lyezkowski, and E. D. Hughes, Fifteenth National Heat Transfer Conf., San Francisco (1975).

.

9.

Y. Taitel and A. E. Dukler, A.I.Ch.E. Journal, 22 (1), 47-55 (1976). F. H. Harlow and A. A. Amsden, Journal of Computational Physics, 1__77, 19-52 (1975).