Non-linear analysis for a natural circulation boiling channel

ELSEVIER

Nuclear Engineeringand Design 152 (1994) 349-360

Nuclear E need.ng andDesign

Non-linear analysis for a natural circulation boiling channel Y . N . Lin, Chin Pan * Department of Nuclear Engineering, National Tsing Hua University, 101, Sect. 2, Kuang Fu Road, Hsinchu 30043, Taiwan

Abstract

The non-linear dynamics of a two-phase natural circulation boiling channel were investigated in this study, by employing the Galerkin nodal approximation method based on the homogeneous flow model. A stability map, which confirms the presence of an instability region under the condition of low exit quality, was obtained. Theoretical results indicate that the oscillation frequency in this low quality region is significantly smaller than that for the density wave oscillation. The model was also utilized to study the transient behavior of the channel following a step change of power. Theoretical results further revealed that an oscillating channel can be stabilized by a change of power from the unstable region to the stable region.

1. Introduction

Two-phase natural circulation loops have received extensive attention because of their simplicity and high heat transfer capability. They are also very attractive for application in nuclear reactors, because of their passive nature. For example, the next generation boiling water reactor employs natural circulation for heat transfer either under normal operational states or abnormal situations. Stable operation of a two-phase natural circulation loop requires a delicate balance between the gravity head available and loop pressure drops. Experimental work reported in previous literature has revealed several different kinds of thermal hydraulic instability in two-phase natural circulation loops. Chexal and Bergles (1973) reported

* Corresponding author. Elsevier Science S.A. SSDI 0029-5493 (94) 00805-9

several different kinds of flow instabilities in addition to the well-known density wave oscillations for a low pressure thermosiphon reboiler. Fukuda and Kobori (1978) identified the 'type-I' instability for a parallel channel, natural circulation loop under the conditions of low exit quality, in addition to the density wave oscillations under the conditions of high exit quality. Recently, Aritomi and his coworkers (1992, 1993) and Chiang et al. (1993) reported three possible varieties of instability during the start-up of a natural circulation boiling water reactor, i.e. geysering, natural circulation oscillations and density wave oscillations. Delmastro et al. (1991) also experimentally verified the presence of instability at low powers in a natural circulation loop with a single boiling channel. More recently, Wang et al. (1994a) reported thermal hydraulic oscillations in low pressure, two-phase natural circulation loops with low powers and high inlet subcooling. Large-

350

Y.N. Lin, C. Pan / Nuclear Engineering and Design 152 (1994) 349-360

amplitude thermal hydraulic oscillations with the presence of reversed flows have been reported. The threshold of instability for a boiling channel may be obtained by linear stability analysis. Standard techniques of linear stability analysis in frequency domain are available in earlier literature (Lahey, 1989). Recently, Wang et al. (1994b) and Lee and Lee (1991) employed this technique in establishing the stability map for a two-phase, natural circulation loop. Non-linear stability analysis for a boiling channel has recently attracted a considerable amount of interest, in that it may provide detailed information regarding oscillation, e.g., frequency and amplitude (Achard, 1985; Rizwan-Uddin, 1986). Chaotic oscillations in two-phase flow systems have also been reported in earlier literature, through non-linear analysis. Rizwan-Uddin and Doming (1988) reported a strange attractor in a boiling channel which is subject to a periodically forced flow. Clausse and Lahey (1980) and Lahey (1991) also found a chaotic attractor in their non-linear analysis of autonomous density wave instabilities. In work by Lahey, a Galerkin nodal approximation was applied to transform the partial differential conservation equations into a set of non-linear ordinary differential equations. Such non-linear analysis has not yet been employed in available literature on a two-phase natural circulation loop. The methodology developed by Clausse and Lahey (1990) is adopted in this study to investigate the non-linear dynamics of a two-phase, natural circulation loop. The characteristics of limit cycle oscillations, as well as unstable oscillations in the unstable operational regions, are studied. Moreover, the transient behaviors following a step change of heating power are also reported.

2. Analysis The non-linear analysis given previously by Clausse and Lahey (1990) is applied for a twophase, natural circulation loop. For a natural circulation loop, however, the mass flow rate is not an independent variable, as in a forced circulation loop; instead, it depends on the power in

the heated section, and on the geometric and operational conditions. This circumstance complicates the analysis and solution procedure, as is discussed later. It is also well known that the loop mass flow rate for a two-phase natural circulation loop increases with an increase in heating power at low powers, where the gravitational pressure drop is dominant in the loop. Furthermore, the mass flow rate decreases with an increase in heating power at high powers, where the two-phase, frictional pressure drop is dominant. Thus, a maximum mass flow rate exists at a certain intermediate heating power. At this power, both the gravitational and two-phase, frictional pressure drops are significant. This maximum mass flow rate may be selected as the reference mass flow rate. Recently, Jeng and Pan (1994) derived the following representative mass flow rate: W s = prAHUs

(1)

where Us is the characteristic velocity in the loop and is given as ( g D u ~ '/2 = \

(2)

For simplicity, it is assumed that the two-phase, frictional factor is twice that of single-phase flow. Thus, we have f2, = 2f~6

(3)

The flow in a natural circulation loop is generally not fully developed, and the single-phase, frictional factor may be evaluated by the following equation that has been developed from a closed-square, natural circulation loop (Su, 1991): 0.213 fl~

-

-

Re0.24 !

(4)

Substituting Eqs. (3) and (4) into Eq. (2) yields the following equation for the characteristic velocity, which is used to normalize the equations: US =

1.62g°'S69D~i 7°5 1j0.137

(5)

The present non-linear analysis focuses on a simplified natural circulation boiling channel as shown in Fig. 1, which consists of a heated section

Y.N. Lin, C. Pan / Nuclear Engineering and Design 152 (1994) 349-360

f-

m

_

351

P r,,Vlt

k

Pr,2

LR/NR

LR

............. j

Pr, I

/

I 0 Q

1

--)

<._

-) -)

<--

pc,he

<-

q" --) -..) .-) -->

Ltt

~-

hi+~Ah~

~-

hi+~Ah~b

~LI¢, =~,(t) ~Ll(t:2(t )

hi

w

Ui

/

Ui

LOOP

NODALIT_ATKW

Fig. 1. Schematic diagram of the two-phase natural circulation loop and nodalization.

and a riser. The present channel analysis approximates a natural circulation loop in which the frictional and form losses in the downcomer are negligibly smaller than those in the heated section and adiabatic riser. The nodalization for the analysis follows that previously suggested by Clausse and Lahey (1990), as is also shown in Fig. 1. The following assumptions are made in the analysis: (1) constant properties are used at a given system pressure; (2) subcooled boiling is not considered; (3) the two-phase flow in the loop is treated by the homogeneous flow model; (4) uniform heat flux distribution is assumed to be in the axial direction; (5) viscous dissipation is neglected; (6) the fluid enters the heated section in a subcooled state, with constant inlet subcooling.

Assuming a linear enthalpy profile between two nodes in the single-phase region, i.e. the region in front of the boiling boundary where the fluid reaches the saturation temperature, and integrating the one-dimensional, single-phase energy equation yields (Clausse, 1990) dLn+ =2ui + - 2NsNVCh ( L + -- L +_l) dt + Nsu b

dL + dt +

n = l . . . . . Ns

(6)

where LJ- is the entrance point, an invariant, and L~vs is the boiling boundary, which is also denoted as 2 +. Eq. (6) provides the dynamics of the boiling boundary, in addition to other nodes in the single-phase region, with constant enthalpy g i v e n as Nsmn h ~+ -

- Ns

hi+

(7)

352

Y.N. Lin, C. Pan / Nuclear Engineering and Design 152 (1994) 349-360

For the two-phase part in the heated region, integrating the continuity equation from the inlet to the exit of the heated section produces (Clausse, 1990)

approximately unchanged following the area change at the exit of the heater, the following equation holds as dz tends to zero:

dM& dt +

Ur'in -~-

-

u; ~ -

p:

1 R--~ue

(14)

(8)

u:

where M& is the total mass in the heated region and is given as (Rizwan-Uddin, 1988) ln( 1/p + ) M& = 2 + + (1 - 2 +) ( l i p +) - 1

(9)

Furthermore, the fluid velocity is not a function of the axial location in the adiabatic riser, so that

1 u~+ =~AA u+

(15)

In Eq. (8), u + is the two-phase mixture velocity at the exit of the heated region, and it can be shown that (Rizwan-Uddin, 1988)

The dynamics of the inlet velocity can be obtained by balancing the dynamic pressure drop in the channel and the available gravity pressure head, which is a constant and is given as

ue+ = ui+ + Npeh(1 -- 2 +)

APext = p f g ( L r + LH)

(10)

In the riser region, the dynamics of two-phase mass in the rth node can be obtained by integrating the continuity equation. Thus, we have dMr+ dt+ - u~+ (p~+_ 1 - p r+)

P.__L_eue Ur'in -- RA Pr,in

1 ~Pr, in dz Pr,in St

(17)

I=AP~-+AP~+AP~-+AP~+AP+~

Each component of the dynamic pressure drop can be obtained by integrating the momentum equations in the heated and riser regions (Clausse, 1990), i.e.

(12)

1 d +" + N P c h ( 1 - 2 + ) " - M ~ ) } AP+ = E----udt + ui Mob+

Here, Eq. (11) is different from that given by Clausse and Lahey (1990), who implicitly assumed that u,+ = u~+ everywhere in the riser; however, Eq. (11) is based on the assumption that Ur+ = U~+/RA as is discussed in what follows. The fluid velocity in the riser is not a function of the axial location in the riser, because the riser is adiabatic, as suggested by Eq. (7a) of Clausse and Lahey (1990). Considering the mass balance in an infinitesimal control volume right next to the exit of the heated section leads to

1

Thus, we have

(11)

where 1 ln(pr+_l/Pr +) Mr+ = R A R e N a ( 1/p r+) - ( l i p +_ 1)

(16)

(13)

where u,,i. is the fluid velocity on the upper control surface of the control volume and Pr.in is the average fluid density of the control volume. The fluid density on the upper control surface is assumed to be that of the average in the control volume. The term dz is the height of the control volume. Assuming that the fluid density remains

(18) AP~-

+ + NexpNsub f 2+ Mch Eu Fr Eu Fr [

1 2Ns

x 2 (2n - 1)(L 2 - L + _ , )

(19)

n=l

1

AP+ = Euu {P+(U+)2 - (u+)Z}

(20)

AP~ A2o( + 2 ln(1/P +) =Al*(ui+Z2+lEu + E u \(ui ) (1 - 2 + ) ( l / p + ) _ 1 2u: N

+

al - 2 + ) ( 1 - M:.) 2 (-1/~+) - - i ( (1/p +) - 1 j

1 + 2Eu {ki(u+)2 - khrP+ (u+)z}

(21)

Y.N. Lin, C. Pan / Nuclear Engineering and Design 152 (1994) 349-360

1

353

0.70

1

AP+ = ~uu [~A { d ~ (U+M{~) + (P~vR--P+)(U+) 2 0.60

Dri + 2 +

1 + A2,~--~n(ue) MR +ke.~p+R(U:) 2 1

1

0.50

+-]

+~rr~AMRJ

(22)

In Eq. (19), we have

0.40 +~:~'0.30

NR

= Y M, +

(23)

r=l

0.20 0.10

3. Solution method

As noted earlier, the mass flow rate (also inlet velocity) is not an independent variable in a natural circulation loop. For the initial steady state conditions, the inlet velocity at a given power and for given geometric and operational conditions is obtained by solving Eq. (17) with time-dependent terms being set to zero, thereby resulting in a transcendental equation for u~. The equation is solved numerically by using the subroutine SNSQE of Kahauer et al. (1989). The subroutine employs the PoweU hybrid scheme--which emerges from the Newton m e t h o d - - t o solve the equation. The non-linear dynamics of the system at a given initial steady state are obtained by solving the set of non-linear ordinary differential equations with a perturbation of the inlet velocity from the steady state. The set of non-linear ordinary differential equations is solved by the subroutine SDRIV2 of Kahauer et al. (1989). This subroutine employs the Gear multi-value method to solve the set of equations. Following the suggestion of Clausse (1990), both the single-phase region and the riser are divided into three nodes.

4. Results and discussion

The non-dimensional steady state inlet velocity is shown in Fig. 2 as a function of the phase change number. In the present study, the non-dimensional phase change number (Np~h) is based

N~.b=1.326

0.00 IIlIIIIIIIIIIIIIIIIIIIIIIIF'IIIlIlIII:IIIIIIIIII 0.00 1.00 2.00 3.00

['~rLe t ~ k ~ f t 5 ['i,r~,e 2 . kc= I 0 [~,rLe 3 ~ kc: 5

4.00

5.00

Npcn Fig. 2. Inlet velocity as a function of heat flux and inlet loss coefficient in a two-phase natural circulation loop.

on the characteristic velocity given by Eq. (5) rather than by the steady state inlet velocity of Clausse and Lahey (1990). Fig. 2 reveals the typical inlet velocity dependence on the heating power. In a lower power region, the inlet velocity increases quite rapidly with an increase in power, and reaches a maximum at a certain power. Subsequently, the inlet velocity decreases with an increase in power. Notably, the non-dimensional inlet velocity is of the order of unity, indicating that the reference velocity given by Eq. (5) is representative. A stability map can be obtained by evaluating the system dynamics following a perturbation of given steady state conditions. The corresponding state is considered to be unstable if the magnitude of oscillation grows continuously following the perturbation. Furthermore, the corresponding state is a state with limit cycle oscillation if the oscillation magnitude grows to a certain value and then remains unchanged. If the magnitude of oscillation dampens out, however, the corresponding state is considered to be stable. Fig. 3 shows the stability maps corresponding to the three different inlet loss coefficients, and other conditions listed in Table 1. This figure reveals the presence of an unstable region at low powers, in addition to the well-known density

Y.N. Lin, C. Pan / Nuclear Engineering and Design 152 0994) 349-360

354

8.00

~I -

t y p e - I / n s t a b i l l t y region stable region density wave oscillation region

III-

~.=I~

6.00

~

1

4.00

III

2.00

0.00

IIII111

I l l l l l l l l l l l l l l l

0.00

1.00

I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 1 1

2.00

3.00

v,,h

¢.00

5.00

6.00

Fig. 3. Stability maps o f a two-phase natural circulation loop.

Table 1 Conditions for the reference case

and K i for these five points are listed in Table 2. A comparison of oscillating behaviors betweenpoints T1 and DI clearly indicates that the system oscillates with a longer period at low powers in the type-I instability region. Notably, TI and D~ have the same difference in Npch from the stability boundary. At an intermediate power level in the stable region, i.e. point So, the system returns to its initial steady state after the initial transient following the perturbation. The small-amplitude limit cycle oscillations exist at the marginal stability boundary, as is evident from a single close loop in the phase diagram of ui+ and 2 + shown in Fig. 6. The non-dimensional frequency along the marginal stability boundary is shown in Fig. 7 as a function of the phase change number. It should be noted that, along the stability boundary, the subeooling number varies as well as the phase change number. This figure clearly indicates that the frequency of oscillation for the type-I instability boundary is quite insensitive to the variations of the phase change (and subcooling) number. In addition, the inlet loss coefficient has an insignificant effect on the frequency of the type-I limit cycle oscillation. This is because the channel pressure drop is dominated by the gravitational drop. Npc h

P LH

68.947 bar 2.44m

khr

2{1 --(I/jI~A)}/RA

LR DH

2.44m 0.0132m

fl~s f2~,

0.213/Re°241 2fi~,

DR

1.4DH

ke~

0

wave oscillations at high powers. The area of the stable region is larger for a system with a larger inlet loss coefficient, as expected. The region of instability at lower powers is referred to as 'type-I' instability in the literature (Fukuda, 1978). In such a region, the pressure drop in the loop is dominated by gravitation, which is in contrast to the density wave oscillations at high powers, in which the two-phase, frictional pressure drop is the dominant contributor. The stability map also can be obtained by linear stability analysis. The same shape of stability map also has been reported elsewhere (Lee, 1991; Wang, 1994). The system oscillation behavior can be quantified by non-linear analyses. Figs. 4 and 5 show the time evolution of the inlet velocity and boiling boundary, respectively, for three of the five points marked in Fig. 3. The values of Nsub,

Y.N. Lin, C. Pan/ NuclearEngineeringandDesign152(1994)349-360 0.5200

0.8400

Tt

7'1

0.5000

0.8200

0.4,800

0.8000

0,4.800

0.7800

0,4400

0.7600

0.7400

0.4-200

N~=3,575

NI~,=3.575 0,4000 0

355

........

2b.b6 . . . . . . . t +

~b.b~ . . . . . . .

0.7200 0.00

~.00

40.00

60.00

f+

0.5100

0.7100

So

0,5000

0.4900

0.6800 0.4800

4.

"< 0.6700

0A.700 0.6600

0,4.600

0.6500

N~e~=4,257 0"45000

0

t÷

'

"

....

2'o.b6

. . . . . . .

,:o.b6

. . . . . . .

~o.00

t÷

0,7111111

0.7000

0.8000

°~'~°o.od

Df

1)t

0.5000 0.4000

+ 0.5000 ,<

O.3000I

e

0A.I~

0.2000I N~,,=4.940 °'~00°o.~o" .......

I07o6 .......

N~,~=4. 940 ~o.b~

'"

.....

do. o

0.3000 t+

Fig, 4. Time evolution of inlet velocity at three different phase

change numbers,

Fig. 5. Time evolution of boiling boundary at three different phase change numbers.

Y.N. Lin, C. Pan / Nuclear Engineering and Design 152 (1994) 349-360

356

The non-dimensional frequency is in the narrow range 0.1-0.3 along the type-I instability boundary. Since the subcooling number along the type-I instability boundary increases following the increase of the phase change number, the insensitivity of the oscillation frequency to the change in the phase change number suggests that the subcooling number has an opposite and competing effect on the oscillation frequency with respect to the phase change number, In contrast, significantly different characteristics of the frequency appear for the density wave oscillation region.

Table 2 Values of Nr~h for points in Fig. 3 Ns=b = 5.966

k i = 15

Point

Npch

T1 T2 Sl So D2 DI

3.575 3.642 3.845 4.257 4.873 4.940

0.39

0.21

0.38

0.20

.,°36

:

\

\

\

0.19

\

0.18 p. 0.17

0.34

0.16

0.33 : ......... , ......... , ......... , ......... , ......... , ......... 0.74 0.75 0.76 0~77 0.78 0.79 0.80

0.15

,T,,,,,,, I..... ,,,,i,, ....... i,,, ...... i,,,, ..... i,,,,,,,,, 0.50 0.52 0.54 0.56 0.58 0.60 0.62

X*

X÷

0.49

0.52

0.50

0.48

0.48 0.47 0.46

..b

+

~0.46 0.44

0.45 0.42

(a)

0.44 , , , , , , , , , i ........ 0.74 0.75

, i , , , , , , , , , i , ,

0.76

~.177

.......

r .........

0.78

0.40

i . . . . . . ,,. 0.79 0.80

(b)

......... , ......... , ................... , ......... , ......... 0.50 0.52 0.54 0.~6 0.58 0.60 0.62

X*

Fig. 6. Limit cycle oscillations (a) at point T 2 and (b) on the stability boundary.

357

Y.N. Lin, C. Pan / Nuclear Engineering and Design 152 (1994) 349-360

1

1.20

• :::: 1.00

.-" : : :

~/'pe-I ~.'tuflal~,t~,t'gregio'n de~eitli 1 were - osoitl~tiorL region 7 k:~=

,5"

/

• ,=,o

k~=15

0.80

0.60

/

/

/

0.85

0.80

I

///

0.75

0.70

0.65

So ---> 8f St ehee~l/ s t ¢ ~

0.40 0.60

......... 0

020

-

.

.

.

i ....... 50

,,i ......... 100

J,,, ...... i,.r,,,,,,i, ....... ~5.~ 200 250

'l''' 300

.

0.54

0.00

. . . . . . . . .

0.00

,

. . . . . . . . .

2.00

0.52

~ . . . . . . . . .

4.00

6.00

NpcH

Fig. 7. Non-dimensionalfrequencyalong the marginal stability boundary.

0.50

0.48

÷ .o46 ! Dramatic changes in frequency arise along the stability boundary. The inlet loss coefficient has a significant effect on the oscillation frequency, which is in contrast to the type-I instability. In addition to the stability analysis described above, the model also may be used to study the transient behavior of a natural circulation loop. Fig. 8 shows that, in the stable region, the system evolves to a new steady state following a step change in the heating power. This figure reveals that the system approaches the new steady state oscillatorily with a strong damping effect. The figure also indicates that the system moves from one point attractor to another point attractor. Fig. 9 illustrates the transient behavior for the system experiencing a step change in the heating power from a stable state to an unstable state. The full line represents the transient following the step change at t + = 100, and the broken line displays the time evolution of the system variables following a small perturbation of the steady state initially at t ÷ = 0. This figure reveals that the system evolves into the oscillating behavior of the unstable state following a change of power from a stable to an unstable state. The phase diagram reveals that the system changes from one point attractor to a repeller.

0.44 0.42 -= >

° 4°.6~"''b'.~'" 'E#C"6.'~'~""b'.'f6"~'~'£"'~3~"''~.'~'" 'E~T" '~.'~'" 'b'.'82 Fig. 8. Transient behavior followinga step change of heating power from a stable state to another stable state.

Fig. 10 shows the transient behavior of the system following a step change of the heating power from an unstable state to a stable state. This figure clearly indicates that the system evolves into the new stable state oscillatorily, with a strong damping effect. The system changes from a repeller to an attractor, as shown in the phase diagram. Therefore, the possibility that an oscillatory s y s t e m can be stabilized by a change in heating power is quite likely.

5. Conclusions

A non-linear analysis for a two-phase, natural circulation channel was performed in this study using a Galerkin nodal approximation method.

358

Y.N. Lin, C. Pan / Nuclear Engineering and Design 152 (1994) 349-360

1O0

0.90

o.go

0.80

0.80 0.70 0.70 + 0.60 ,,<

4'~

0.60

0.50

0.50

- -

0.40

0.40

~q# - - >

T,

0.30

40

0

t.a0

Tf

-

0.30

12o

..............

160

i ..............

30

i, .............

Dt $o S, steady stale - - >

-

i .......

60

90

, ......

i ..............

120

150

t*

0.70

0.60

0.60

0,55

0.50

0.50 0.40 0.45

÷ff

+~:~j'0.30

0.40

0.20 0.35

0.10 --:, 0.00

.........

0,30

~. . . . . . . . .

0.40

, .........

0.5o

~. . . . . . . . .

0.50

, .........

0.70

~. . . . . . . . .

0.80

D1

T1 ~. . . . . . . . .

o.go

0.30

1.00

.........

0.40

i .........

0.45

i .........

0.50

i .........

0.55

t .........

0.6~+

i .........

0.65

i .........

0.70

- - >

z .........

0.75

So i .........

0.80

0.85

X÷

Fig. 9. Transient behavior following a step change of heating power from a stable state to an unstable state.

Fig. 10. Transientbehavior following a step change of heating power from an unstable state to a stable state.

The following conclusions could be drawn from the present study. (1) The present non-linear study confirms the presence of 'type-I' instability under the conditions of low quality at the exit of the heated section. (2) The oscillation behavior of 'type-I' instability is fundamentally different from the conventional density wave oscillations. The frequency is lower than that in the density wave oscillations. (3) The model is capable of predicting the loop transient behavior following a change in heating power. The loop approaches the new steady state oscillatorily, with a strong damping effect in the stable region. Meanwhile, the system evolves into the oscillating behavior of the unstable state fol-

lowing a change of power from a stable state to an unstable state. However, the system evolves into the new stable state oscillatorily, with a damping effect, following a step change of the heating power from an unstable state to a stable state. This occurrence indicates the possibility that an oscillatory system can be stabilized by a change in power.

Acknowledgements This work was supported by grants from the National Science Council of the Republic of China under Contracts NSC81-0401-E007-18 and NSC82-0413-E007-081.

359

Y.N. Lin, C. Pan / Nuclear Engineering and Design 152 (1994) 349-360

Appendix A: Nomenclature Ax.s

cross-sectional area of the heated section

AR

cross-sectional area of the riser section

Cpf

liquid constant pressure specific heat (J kg -1 K -l) diameter of the heated section (m) friction factor or frequency gravitational acceleration (m s -2) enthalpy (J kg -I) latent heat of evaporation (J kg -~) loss coefficient length (m) mass (kg) number of nodes in the single-phase region number of nodes in the riser region heated perimeter (m) pressure drop (Pa) heat flux (W m -2) heating power (W) riser-to-heated section area ratio Reynolds number riser-to-heated section length ratio time (s) velocity (m s -l) specific volume of saturated liquid (m 3 kg-1) difference in specific volume of saturated liquid and vapor (m 3 kg -1) mass flow rate (kg s -1) axial coordinate (m)

(m 2) (m 2) DH f g h hfg K L M Ns NR PH AP q" Q RA

Re RL

t u vf

/)fg W z

Greek letters

P Pf

2

liquid thermal expansion coefficient (K -l) density (kg m -3) density of saturated liquid (kg m -3) boiling boundary (m)

H i I g n r R l~b 2~b 0 hr

heater inlet of heated section inertial pressure drop gravitational pressure drop nth node in the single-phase region rth node in the riser region riser single phase two phase steady state interface of heated section and riser

Dimensionless groups

Q0

hfg

Fr =

Al¢ -----flcLH single-phase friction number 2DH A2~, =f2~Lrt 2DH

Cpf

Vfg

Eu = APext pfu 2

Re-

thermal expansion number

Euler number

usDH

vf

Reynolds number (Re is based on ui during dynamic calculation)

f+ =fLu/u s h + = (h - h f ) l { Q o / ( p f a x _ s U s )

L +n = L . L n M

M +_

pfLHAx-s

t ÷ = t/(LH/us)

z + = Z/LH

ext f

two-phase friction number

flhfg o f

Nexp----

a e ex

of

Froude number

gL.

U + = U/U S

acceleration pressure drop channel exit of heated section exit of riser section external frictional pressure drop

phase change number

Nsub h f - hi Vfg subcooling number

Subscripts

ch

vfg

Nr~h -- prAx_sUs hfgvf

Ap + =

AP

APoxI

2 + = L~vs = L N s / L H P+ = P/Pf

}

360

Y.N. Lin, C. Pan / Nuclear Engineering and Design 152 (1994) 349-360

References J.L. Achard, D.A. Drew and R.T. Lahey, The analysis of nonlinear density-wave oscillations in boiling channels, J. Fluid Mech. 155 (1985) 213-232. M. Aritomi, J.H. Chiang, T. Nakahashi, M. Wataru and M. Moil, Fundamental study on thermo-hydraulics during start-up in natural circulation boiling water reactors (I), thermo-hydraulic instabilities, J. Nucl. Sci. Technol. 29 (1992) 631-641. M. Aritomi, J.H. Chiang and M. Moil, Geyseilng in parallel boiling channels, Nucl. Eng. Des. 141 (1993) 111 - 121. V.K. Chexal and A.E. Bergles, Two-phase instabilities in a low pressure natural circulation loop, AIChE. Symp. Ser. 69 (1973) 37-45. J.H. Chiang, M. Aritomi and M. Moil, Fundamental study on thermo-hydraulics during start-up in natural circulation boiling water reactors (II), natural circulation oscillation induced by hydrostatic head fluctuation, J. Nucl. Sci. Technol. 29 (1993) 203-211. A. Clausse and R.T. Lahey, Jr., An investigation of periodic and strange attractors in boiling flows using chaos theory, Proc. 9th Int. Heat Transfer Conf., Jerusalem, 1990, Vol. 2, pp. 3-8. D.F. Delmastro, A. Clausse and J. Converti, The influence of gravity on the stability of boiling flows, Nucl. Eng. Des. 127 (1991) 129-139. K. Fukuda and T. Kobori, Two-phase flow instabilities in parallel channels, Proc. Sixth Int. Heat Transfer Conf., Vol. 1, 1978, pp. 369-374. H.R. Jeng and C. Pan, Multiple stcady-state solutions for a two-phase natural circulation loop with subcooled boiling, Proc. 4th Int. Topical Meet. on Nuclear Thermal

Hydraulics, Operations and Safety, Taipei, 5-8 April 1994, paper 9-A. D. Kahaner, C. Moler and S. Nash, Numerical Methods and Software, Prentice-Hall, Englewood Cliffs, NJ, 1989. R.T. Lahey, Jr., An application of fractal and chaos theory in the field of two-phase flow and heat transfer, Wiirme- und StoJ~bertragung 26 (1991) 351-363. R.T. Lahey, Jr. and M.Z. Podowski, On the analysis of various instabilities in two-phase flows, Ch. 3, in G.F. Hewitt, J.M. Delhaye and N. Zuber (eds.), Multiphase Science and Technology, vol. 4, Hemisphere, 1989. S.Y. Lee and D.W. Lee, Linear analysis of flow instability in an open two-phase natural circulation loop, Nucl. Eng. Des. 28 (1991) 317-330. Rizwan-Uddin and J.J. Doming, Some nonlinear dynamics of a heated channel, Nucl. Eng. Des. 93 (1986) 1-14. Rizwan-Uddin and J.J. Doming, A chaotic attractor in a periodically forced two-phase flow system, Nucl. Sei. Eng. 100 (1988) 393-404. Y.J. Su, C.J. Chang and C.H. Lee, Experimental study of natural circulation in a rectangular closed loop, Proc. 2nd World Conf. on Experimental Heat Transfer, Fluid Mechanics and Thermodynamics, Dubrovnik, 1991. S.B. Wang, C.Y. Wu, C. Pan and W,K. Lin, Thermal hydraulic oscillations in a low pressure two-phase natural circulation loop at low powers and high inlet subeooling, Proc. 4th Int. Topical Meet. on Nuclear Thermal Hydraulics, Operations and Safety, Taipei, 5-8 April 1994a, paper 53-C. F.S. Wang, L.W. Hu and C. Pan, Thermal and stability analysis of a two-phase natural circulation loop, Nucl. Sci. Eng. ll7 (1994b) 33-46.