COMPUTATIONAL FLUID DYNAMICS MODELING OF TWO-PHASE FLOW IN A BOILING - - PDF document

Mathematics and Computation, Supercomputing, Reactor Physics and Nuclear and Biological Applications Palais des Papes, Avignon, France, September 12-15, 2005, on CD-ROM, American Nuclear Society, LaGrange Park, IL (2005)

COMPUTATIONAL FLUID DYNAMICS MODELING OF TWO-PHASE FLOW IN A BOILING WATER REACTOR FUEL ASSEMBLY

Adrian Tentner Nuclear Engineering Division Argonne National Laboratory 9700 South Cass Avenue, Argonne, Illinois 60439, USA tentner@anl.gov Simon Lo CD-adapco 200 Shepherds Bush Road, London W6 7NY, UK simon.lo@uk.cd-adapco.com Andrey Ioilev, Maskhud Samigulin, Vasily Ustinenko Institute of Theoretical and Mathematical Physics Russian Federal Nuclear Center (VNIIEF) Mir ave. 37, 607190 Sarov, Nizhnii Novgorod Region, Russian Federation ioilev@vniief.ru, sms@vniief.ru, ustinenko@vniief.ru ABSTRACT Two-phase flow phenomena inside a BWR fuel bundle include coolant phase changes and multiple flow regimes which directly influence the coolant interaction with fuel assembly and, ultimately, the reactor performance. The CFD-BWR code is being developed as a specialized module built on the foundation of the commercial CFD code STAR-CD which provides general two-phase flow modeling capabilities. New models describing two-phase flow and heat transfer phenomena specific for BWRs are developed and implemented in the CFD-BWR module. A set

• f experiments focused on two-phase flow and phase-change phenomena has been identified for

the validation of the CFD-BWR code and results of several experiment analyses are presented. The close agreement between the computed results, the measured data and the correlation results provides confidence in the accuracy of the models. KEYWORDS: Boiling Water Reactor, Two-Phase Flow, Computational Fluid Dynamics, Flow Regimes.

• 1. INTRODUCTION

It is highly desirable to understand the detailed two-phase flow phenomena inside a Boiling Water Reactor (BWR) fuel bundle. These phenomena include coolant phase changes and multiple flow regimes which directly influence the coolant interaction with fuel assembly and, ultimately, the reactor performance. Traditionally, the best analysis tools for the analysis of two-phase flow phenomena inside the BWR fuel assembly have been the sub-channel codes.

However, the resolution of these codes is too coarse for analyzing the detailed intra-assembly flow patterns, such as flow around a spacer element and it has been recognized that their basic modeling approach and computational methods no longer represent the state-of-art in the field of numerical simulation [1,2]. Recent progress in Computational Fluid Dynamics (CFD), coupled with the rapidly increasing computational power of massively parallel computers, shows promising potential for the fine-mesh, detailed simulation of fuel assembly two-phase flow

• phenomena. However, the phenomenological models available in the 3-dimensional CFD

programs are not as advanced as those currently being used in the sub-channel codes used in the nuclear industry. In particular, there are no models currently available which are able to reliably predict the nature of the flow regimes, and use the appropriate sub-models for those flow regimes. The CFD-BWR code is being developed as a specialized module built on the foundation of the commercial CFD code STAR-CD which provides general two-phase flow modeling capabilities. New models describing specific for BWR two-phase flow and heat transfer phenomena are developed and implemented in the CFD-BWR module, which interacts closely with the standard STAR-CD code to allow the study of BWR fuel assembly performance.

• 2. PHENOMENOLOGICAL MODELS

The key phenomenological models included in the CFD-BWR code are focused on the prediction of local two-phase flow regimes in a BWR fuel bundle and the definition of the appropriate mass, momentum and energy inter-phase and phase-to-boundary transfer terms. The strategy for the prediction of the local flow regime includes the use of local flow regime maps and flow regime specific phenomenological models in conjunction with an interface transport and topology transport approach. During the first phase of the project specific mass, momentum, and energy inter-phase exchange and phase-to-boundary exchange terms have been developed for the bubbly boiling flow regime, which is typical for the initial two-phase region in BWR channels and is also of interest in the analysis of PWR channels. These models are described

• below. A procedure for the identification of the local flow regime using a flow regime map based
• n the local void fraction and void fraction gradients has also been developed and its initial

implementation is described. 2.1. Transport Equations The STAR-CD Eulerian two-phase solver tracks the mass, momentum, and energy of the liquid and vapor phases in each cell. Full details of the Eulerian two-phase flow models in STAR-CD can be found in [3] and [4]. The main equations solved are the conservation of mass, momentum and energy for each phase. Mass Conservation The conservation of mass equation for phase is: k

( ) ( ) (

∑

=

− = ∇ + ∂ ∂

N i ik ki k k k k k

m m u t

1

. & & ρ α ρ α

)

(1)

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 2/17

where

k

α is the volume fraction of phase , k

k

ρ is the phase density, is the phase velocity, and are mass transfer rates to and from the phase, and is the total number of phases. The sum of the volume fractions is clearly equal to unity.

k

u

ki

m &

ik

m & N 1 =

∑

k k

α (2) Momentum Conservation The conservation of momentum equation for phase is: k

( ) ( )

( ) ( )

M g p u u u t

k k k t k k k k k k k k k k

+ + ∇ − = + ∇ − ∇ + ∂ ∂ ρ α α τ τ α ρ α ρ α . . (3) where

k

τ and are the laminar and turbulence shear stresses respectively,

t k

τ p is pressure, g is gravitational acceleration and M is the sum of the inter-phase forces. Energy Conservation The conservation of energy equation for phase is: k

( ) ( ) ( )

Q T e u e t

k k k k k k k k k k

= ∇ ∇ − ∇ + ∂ ∂ λ α ρ α ρ α . . (4) where is the phase enthalpy,

k

e

k

λ is the thermal conductivity, is the phase temperature and is the inter-phase heat transfer.

k

T Q 2.2. Turbulence Equations To calculate the continuous and dispersed phase turbulence stresses used in equation 3 above values for and k ε are required. These can be computed using the extended

• k ε equations

containing extra source terms that arise from the inter-phase forces present in the momentum

• equations. The additional terms account for the effect of bubbles on the turbulence field. The

relevant equations are:

( )

( )

2

. .

k c c k t c c c c c c c c

S G k k u k t + − + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∇ + ∇ = ∇ + ∂ ∂ ε ρ α σ µ µ α ρ α ρ α (5)

( )

( )

2 2 1

. .

ε ε

ε ρ α ε σ µ µ α ε ρ α ε ρ α S C G C u t

c c t c c c c c c c c

+ − + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∇ + ∇ = ∇ + ∂ ∂ (6) where

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 3/17

( ) (

k C A u u A S

t D d c d d c t c D k

1 2 .

2

− + ∇ − − = α σ α α ν

α

)

(7)

( )ε

ε

1 2

2

− =

t D C

A S (8)

( )

c T c c c

u u u G ∇ ∇ + ∇ = : µ (9) In the above equations, is a response coefficient defined as the ratio of the dispersed phase velocity fluctuations to those of the continuous phase:

t

C

c d t

u u C ′ ′ = (10) The dispersed-phase turbulent stress is correlated to the continuous-phase turbulent stresses via the response coefficient such that

t d

τ

t c

τ

t

C

t c t c d t d

C τ ρ ρ τ

2

= (11) Further details of the response coefficient and the turbulence model can be found in [3] and [4].

t

C 2.3. Inter-phase Transfer for the Bubbly Flow Regime The inter-phase mass, momentum, and energy exchanges depend on the local geometry and thermo-hydrodynamic conditions. During the first phase of this work the inter-phase transfer terms for the bubbly flow regime have been developed and implemented in the CFD-BWR

• module. For computational cells where the bubbly flow regime is present, the vapor is assumed

to exist in the form of spherical bubbles with a variable diameter. While the bubble diameter can vary from cell to cell, all bubbles in one cell are assumed to have the same diameter. Vapor bubbles with a prescribed diameter are generated near the heated surfaces and are entrained in the coolant stream, their trajectories being determined by the inter-phase forces. The inter-phase forces included in the computations are buoyancy, drag, turbulence drag, lift, and virtual mass

• forces. As they exchange energy and mass with the surrounding liquid the bubbles can condense

and decrease in size and number or, under certain conditions, grow due to additional liquid

• evaporation. A diagram illustrating the heat and mass exchanges between a vapor bubble and the

surrounding liquid is presented in Fig. 1. The inter-phase heat and mass transfer models were

• btained by considering the heat transfers from the gas and the liquid to the gas/liquid interface.

The net heat transfer to the interface is used to compute the mass transfer rate between the two phases.

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 4/17

Figure 1. Heat and mass transfer between a vapor bubble and liquid Heat transfer rate from the liquid to the interface is:

( )

sat l d l l

T T A h q − = & (12) Heat transfer rate from the gas to the interface is:

( )

sat g d g g

T T A h q − = & (13) where is heat transfer coefficient, is interfacial area, h

d

A T is temperature, is the saturation temperature, subscript denotes liquid and

sat

T l g denotes gas phase. Assuming that all the heat transferred to the interface is used in mass transfer (i.e. evaporation or condensation), the mass transfer rate can be written as:

fg g l

h q q m & & & + = (14) where is the latent heat.

fg

h To determine the interfacial area used in Eqs. 12 and 13 we need to specify the bubble

• diameter. Since bubbles are generated due to boiling and reduced due to condensation, the

bubble diameter is expected to be function of the liquid temperature. Kurul and Podowski [6] defined the local bubble diameter using linear interpolation between measured bubble diameters at two specified values of liquid sub-cooling:

d

A

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 5/17

( ) ( )

1 1 1

T T T T d T T d d

sub sub

∆ − ∆ ∆ − ∆ − ∆ − ∆ = ;

1

T T T

sub

∆ ≥ ∆ ≥ ∆ (15) d d = ; T Tsub ∆ ≥ ∆ (16)

1

d d = ;

1

T Tsub ∆ ≤ ∆ (17) The liquid sub-cooling is defined as:

l sat sub

T T T − = ∆ (18) The following values obtained from experiments were used:

4

10 5 . 1

−

= x d at 5 . 13

0 =

∆T

3 1

10 2

−

= x d at 5

1

− = ∆T The inter-phase forces considered in the model are: drag, turbulent drag, virtual mass and lift forces, and momentum transfer associated with mass transfer, hence

(

∑

=

− + + + + =

N i k ik i ki L M T D

u m u m F F F F M

1

& &

)

(19) Further details of the inter-phase forces can be found in [3] and [4]. 2.4. Wall Heat Partitioning Model A model describing the heat transfer between the heated wall and the coolant has also been

• developed. The heat flux from the wall is divided into three parts according to a wall heat

partitioning model which includes convective heat for the liquid, evaporative heat for generation

• f steam and quench heat for heating of liquid in the nucleation sites. If the wall heat flux is

specified, rather than the wall temperature, this model allows the calculation of the wall temperature that corresponds to the specified heat flux. The details of the wall heat partitioning model are given below. As we know, evaporation starts from nucleation sites at the heated surface, we can divide the heated surface into two parts: (1)

• fraction of wall area not covered by nucleation site and therefore subjected to convective

heat transfer.

c

A (2)

• fraction of wall area covered by nucleation sites and therefore subjected to evaporative

heat transfer.

e

A n d F A

w A e

′ ′ = 4

2

π (20)

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 6/17

e c

A A − =1 (21) where is the model constant in which a value of 2 has produced satisfactory results in our validation exercises, is the bubble departure diameter and

A

F

w

d n′ ′ is the nucleation site density. The nucleation site density can be obtained from:

( )

p

T m n

sup

∆ = ′ ′ (22) where is the wall superheat, and according to Kurul and Podowski [6]

sat wall

T T T − = ∆

sup

185 = m and . 805 . 1 = p The bubble departure diameter is obtained from: ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∆ ∆ − =

0 exp

T T d d

sub w

(23) where and are model constants. mm d 6 .

0 =

45

0=

∆T Over the area not covered by nucleation sites we have convective heating of the liquid. The convective heat flux can be modeled as:

( )

l wall c c c

T T A h q − = ′ ′ & (24) Over the area covered by nucleation sites, we have evaporation heat transfer. The evaporation heat flux at the wall is proportional to the nucleation sites density, the bubble departure diameter and bubble departure frequency, , as follow: f n f h d q

fg g w e

′ ′ = ′ ′ ρ π 6

3

& (25) The bubble departure frequency can be obtained from:

( )

l w g l

d g f ρ ρ ρ − = 3 4 (26) As the bubble detaches from the wall, the space it occupied is filled by cooler water. Part of the wall heat flux is used in heating this replacement water. We call this heating the quenching heat transfer, Del Valle and Kenning [7] modeled this heat transfer by transient heat conduction in a semi-infinite slab:

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 7/17

( )

l wall e q q

T T A h q − = ′ ′ & (27) π λ ρ

l l l w q

Cp t f h 2 = (28) where is the waiting time between the bubble departure and the activation of the next bubble,

w

t f tw 8 . = (29) The wall heat flux is therefore made up of three components as follow:

q e c w

q q q q ′ ′ + ′ ′ + ′ ′ = ′ ′ & & & & (30) 2.5. Flow Regime Map In a BWR fuel assembly the coolant exhibits a wide range of flow regimes as it heats up and

• vaporizes. In order to model correctly the inter-phase transfer terms we must determine first the

local flow regime. An initial procedure for the identification of the local flow regimes characteristic for a BWR fuel assembly has been developed; it is presented in this section. The flow regime identification is based on a flow regime map which is used for all inner cells. Additional assumptions must be used in the cells adjacent to the walls in order to separate the liquid film when the film thickness is smaller than boundary cell size in the direction orthogonal to the wall, and these assumptions are not used for inner cells. So there are different flow regime maps for inner cells and for cells adjacent to rigid walls. The flow regime map for inner points of the numerical mesh is shown in Fig. 2. The flow regime is defined by two quantities: void fraction

g

α and void fraction increment

g

α δ γ ∇ ⋅ = at a distance of characteristic cell size δ. Currently, the cell size is defined as minimum of the cell dimension in three directions:

{ }

z y x

∆ ∆ ∆ = , , min δ . Three flow regimes are distinguished in the map:

• 1. Regimes with multi-connected interface inside a cell:
• Bubbly flow – gaseous phase is distributed in liquid in the form of bubbles with diameter

smaller than the cell size;

• Mist flow – liquid phase is distributed in gas in the form of droplets with diameter

smaller than the mesh cell size;

• Froth flow is an intermediate regime between bubbly and mist flows.
• 2. Regimes with simply-connected interface in a cell that can be approximated by a plane – this

is sharp interface.

• 3. Intermediate regimes between the first and the second cases.

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 8/17

This map is used to determine inter-phase interactions. In sharp interface regime it is assumed that interfacial area is flat and is oriented orthogonally to the vector

g

∇α . Interfacial area can be found from this condition. Interfacial drag in case of sharp interface can be represented as a total

• f components along and orthogonally to interfacial area, which differ from each other.

Correlations available from published literature will be used for inter-phase interactions in bubbly and mist regimes. In transitional regimes, required quantities will be obtained by interpolating between the basic regimes.

Sharp interface (3) Transition (21) Transition (22) Transition (23) Bubbly (11) Transition (churn) (12) Mist (13)

αg,bc αg,cm γ γ2 γ1 αg Wall w w -1

g

α

f

δ

g,w 1 −

α

g,w

α

Figure 2. Flow regime map. Figure 3. Schematic of flow restructuring in a boundary cell. In cases when the liquid film thickness is smaller than near-wall cell size (normal to the wall) the above described flow regime map is unsuitable. A special flow regime is used in such cases - regime 31. This regime is defined in a boundary cell if:

• the void fraction in the cell is greater than zero,
• the void fraction in the adjacent cell (in the direction normal to the wall) is higher than that in

the boundary cell, and

• the void fraction in the adjacent cell (in the direction normal to the wall) is higher than a

predetermined value (currently, a value of 0.7 is used). If all these conditions are matched, it is assumed that in this boundary cell a liquid film is present

• n the wall and coexists with the flow regime corresponding to the void fraction in the adjacent

cell (in the direction normal to the wall). The flow configuration in the boundary cell is restructured according to these assumptions as shown in Fig. 3. Solid horizontal lines in this figure show void fraction values in the boundary cell and the adjacent cell. Dashed line designates flow restructuring in the boundary cell. A liquid film is located near the wall (void fraction is zero), while void fraction in the rest of the cell is equal to that in the adjacent cell. According to this restructuring, the film thickness δf , which is required for closure relationships, can be calculated.. To verify the flow regime map algorithm in a stand-alone mode, it was applied to a specified void fraction field in a heated vertical channel with upward coolant flow. Fig. 4a presents specified void fraction values on the mesh and Fig. 4b shows the flow regimes determined by the

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 9/17

flow regime map program. The following flow regime map parameters have been chosen for the current version: 3 .

,

=

bc g

α , 8 .

,

=

cm g

α , 2 .

1 =

γ , 4 .

2 =

γ . In the future these parameters will be adjusted based on experimental data. Flow regimes are marked in accordance with the flow regime map notation given in Fig. 2. It is seen that the program determines flow regimes rather accurately.

z/r 1 2 3 4 5 6 31 0.9 0.9 0.9 0.9 0.9 0.9 30 0.9 0.9 0.9 0.9 0.9 0.9 29 0.9 0.9 0.9 0.9 0.9 0.9

mist

28 0.9 0.9 0.9 0.9 0.9 0.8 27 0.9 0.9 0.9 0.9 0.8 0.7 26 0.9 0.9 0.9 0.9 0.7 0.6 25 0.9 0.9 0.9 0.8 0.7 0.6

annular-mist

24 0.9 0.8 0.8 0.7 0.5 0.4 23 0.8 0.8 0.7 0.6 0.5 0.4 22 0.8 0.7 0.5 0.5 0.4 0.3 21 0.7 0.6 0.5 0.5 0.4 0.3

transition

20 0.6 0.5 0.4 0.4 0.4 0.3 19 0.0 0.1 0.1 0.2 0.2 0.3 18 0.0 0.0 0.1 0.2 0.2 0.3 17 0.0 0.1 0.1 0.2 0.2 0.3 16 0.5 0.5 0.4 0.4 0.3 0.3 15 0.4 0.5 0.5 0.4 0.3 0.3 14 0.4 0.4 0.5 0.5 0.3 0.3

slug

13 0.3 0.4 0.5 0.5 0.3 0.3 12 0.3 0.3 0.4 0.5 0.5 0.3 11 0.2 0.2 0.3 0.3 0.5 0.3 10 0.1 0.2 0.2 0.2 0.4 0.3 9 0.1 0.1 0.1 0.2 0.2 0.3

bubbly

8 0.1 0.1 0.1 0.1 0.2 0.3 7 0.2 0.2 0.1 0.1 0.1 0.2 6 0.0 0.0 0.0 0.0 0.1 0.2 5 0.0 0.0 0.0 0.0 0.0 0.1

bubbly sub- cooled

4 0.0 0.0 0.0 0.0 0.0 0.0 3 0.0 0.0 0.0 0.0 0.0 0.0 2 0.0 0.0 0.0 0.0 0.0 0.0 1 0.0 0.0 0.0 0.0 0.0 0.0

water

z/r 1 2 3 4 5 6 31 13 13 13 13 13 13 30 13 13 13 13 13 13 29 13 13 13 13 13 13 28 13 13 13 13 13 31 27 13 13 13 13 12 31 26 13 13 13 13 12 31 25 13 13 13 12 12 31 24 13 12 12 12 12 12 23 12 12 12 12 12 12 22 12 12 12 12 12 11 21 12 12 12 12 12 11 20 22 22 22 12 12 11 19 11 21 11 11 11 11 18 11 11 11 11 11 11 17 11 21 11 11 11 11 16 22 12 12 12 11 11 15 12 12 12 12 11 11 14 12 12 12 12 11 11 13 11 12 12 12 11 11 12 11 11 12 12 12 11 11 11 11 11 11 12 11 10 11 11 11 11 12 11 9 11 11 11 11 11 11 8 11 11 11 11 11 11 7 11 11 11 11 11 11 6 00 00 00 00 11 11 5 00 00 00 00 00 11 4 00 00 00 00 00 00 3 00 00 00 00 00 00 2 00 00 00 00 00 00 1 00 00 00 00 00 00

Figure 4. a. Void fraction distribution (left) and b. Flow regime distribution (right) in the model problem: 00 - subcooled liquid, 11-bubbles, 13 – mist, 30 – sharp interface, 31 – film, 12, 21, 22, 23 – transition.

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 10/17

• 3. MODEL VALIDATION AND VERIFICATION

A set of experiments designed for the study of two-phase flow and phase-change phenomena has been identified for the validation of the CFD-BWR code. The initial analysis has been focused

• n the validation of the bubbly flow regime models described above. This section presents the

results of three such experiment analyses. 3.1. Validation Case 1 – Bartolomei and Chanturiya experiments The model described above was used for the analysis of the Bartolomei and Chanturiya [8] experiments involving upward water flow with boiling in a 2 m long heated pipe of 15.4 mm

• diameter. Pressure in the pipe was maintained at p=45 bar. The inlet water has a sub-cooling of

∆Ts=60°K and a mass flux of G=900 kg/m2s. A heat flux of q=0.57 MW/m2 was applied uniformly along the pipe wall. Axial distributions of the wall temperature, bulk liquid temperature and bulk void fraction were measured and therefore available for comparison with the computed results. Figs. 5 and 6 show that the calculated wall temperature, bulk temperature and void fraction profiles are in good agreement with the measured data. The calculated bulk temperature and void fraction were obtained at each axial location from the more detailed radial distributions calculated by the code. Although these calculations were performed using a specified heat flux boundary condition, the wall-heat partitioning model described in Section 2.4 calculates the heated wall temperature in addition to the vapor source. The calculated wall temperatures are compared in Fig. 5 with the measured wall temperatures and with the wall temperatures given by Dittus-Boelter correlation [9] for the single-phase entrance region and Thom correlation [9] for the two-phase region.

Thom correlation Dittus-Boelter

Figure 5. Bulk coolant temperature and wall temperature axial distributions. Figure 6. Bulk void fraction axial distribution.

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 11/17

3.2. Validation Case 2 – Bartolemei and Gorburov experiments The same bubbly flow boiling model was used for the analysis of the Bartolemei and Gorburov [10] experiments which studied upward flow with condensation of steam in a 2m-long adiabatic (none heated) pipe with diameter D=32 mm. A steam-water mixture was created in a mixing chamber attached to the bottom of the pipe. Due to the subcooling of the water, steam is condensed as the mixture moves up the pipe. A schematic representation of the experimental section is presented in Fig. 7.

Steam Mixing chamber Subcooled water Void fraction measurement section

0.25 0.5 0.75 1 1.25 1.5 1.75 2 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4

α l, m

experim ent 5-a S tar-C D experim ent 5-a experim ent 7-b S tar-C D experim ent 7-b

Figure 7. Scheme of experimental section for study of condensation. Figure 8. Void fraction as a function of length along the pipe in test-cases 5-a and 7-b. In the experiments the average void fraction along the pipe was measured, and these data are used for comparison with calculated results. Calculations were performed for two experiments:

• a. Experiment 5-a: pressure p=20 bar, mass flux of water G=250 kg/m2⋅s; water subcooling and

void fraction at the mixing chamber outlet ∆Ts=15°K and α=0.25.

• b. Experiment 7-b: pressure p=30 bar, mass flux of water is G=420 kg/m2⋅s; water subcooling

and void fraction at the mixing chamber outlet ∆Ts=35°K and α=0.4. Inter-phase heat transfer was calculated with constant Nusselt number Nu=1, bubble diameter was assumed constant, 0.0013 m (according to the diameter of holes in the mixing chamber used in the experiment). The adiabatic wall boundary condition is used: q=0. Results of comparison between numerical and experimental results are shown in Fig. 8. Good agreement between numerical and experimental results has been achieved. Discrepancy is

• bserved only at low void fractions. A possible way to improve the results at low void fractions

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 12/17

is to introduce void-fraction-dependent bubble diameter during condensation below a specified threshold void fraction. 3.3. Validation Case 3 – Avdeev, Pekhterev and Bartolemei experiments The bubbly flow boiling model was also used for the analysis of experiments reported by Avdeev and Pekhterev [11] and Bartolemei, et al. [12]. In these experiments illustrated schematically in Fig. 8 the average void fraction was measured over the pipe length in upward water flow. In the heated lower section of the pipe sub-cooled boiling occurs and steam is generated. The section above is adiabatic and vapor condensation occurs due to the mixing of the vapor generated near the wall with the still subcooled liquid core. The pipe diameter is D=12.03 mm, the heated section length is L0=1 m and the total pipe length is L=1.4 m. Scheme of experimental section is presented in Fig. 9. In our calculations we used the same models as in the test-cases described above except for the bubble size, for which we used Kurul-Podowski correlation [6] with bubble size dependent on liquid subcooling.

L0 L

Wall heat Flux, q Subcooled water Water + steam

z

Figure 9. Scheme of experimental section for study of boiling and condensation. Water at temperature Tl enters the pipe at pressure Pl=6.89 MPa and mass flux Gl. Heat flux at the channel wall in the heated section of the pipe is q=const. The degree of subcooling of water at the inlet is given by ∆Tsub. Saturation temperature for the conditions is Tsat=558°K. Values of Tl, Gl, q and ∆Tsub for the experiments discussed are given in Table I. The wall heat boundary condition in the non-heated upper section of the pipe is q=0. Table I. Experimental parameters Experiment # q, MW/m2 Gl, kg/(m2⋅s) Tl, °K ∆Tsub, °K 2 1.2 1500 495 63 3 0.8 1500 519 39 5 0.8 1000 503 55 Typical distributions of the water temperature and void fraction calculated for experiment # 2 aree presented in Figures. 10a and 10b, respectively. Figures 11a and 11b show the radial void fraction and water temperature distribution, respectively, calculated for experiment # 2 at three axial locations. The elevation of 0.6 m corresponds to the beginning of boiling, elevation of

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 13/17

0.95 m is located near the end of the heated section in the boiling region, and elevation of 1.3 m is located near the end of the experimental section in the condensation region.

• Figs. 10 and 11 illustrate development of radial heating of water, changing of void fraction and

transport of steam due to lift force and turbulent drag. By the end of the heated section the water near the heated wall reaches saturation temperature, while the water at the center of the pipe remains approximately 30°K sub-cooled. As illustrated in Fig. 11a, the vapor fraction decreases in the adiabatic section of the pipe due to condensation caused by turbulent mixing of the two- phase mixture from the near-wall region with the sub-cooled liquid in the central region. Correspondingly, the radial distribution of water temperature in the condensation section flattens mainly due to turbulent transport, as illustrated in Fig. 11b. Figure 10. a. Distribution of water temperature and b. Void fraction in the flow (simulation of experiment # 2).

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.003 0.0035 0.004 0.0045 0.005 0.0055 0.006

Radial distance (m) V oid fraction Axial distance = 0.6 m Axial distance = 0.95 m Axial distance = 1.3 m

510 515 520 525 530 535 540 545 550 555 560 565 0.001 0.002 0.003 0.004 0.005 0.006

Radial distance (m) T liquid, K Axial distance = 0.6 m Axial distance = 0.95 m Axial distance = 1.3 m Tsaturation

Figure 11. Radial distribution of : a. void fraction and b. water temperature at three axial locations (simulation of experiment # 2).

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 14/17

• Fig. 12 presents heat flux to vapor and

water due to different phenomena described by the wall heat partitioning model described in Section 2.4. Heat transfer in the lower section of the heated pipe occurs mainly due to convective heating of water then, as water is heated, vaporization starts to

• dominate. Starting from elevation of

~0.6 m, the area occupied by bubbles becomes equal to entire pipe surface, and convective heat flux to water becomes zero. Water is heated only through quenching heat transfer, i.e. in the interval between bubble departure and formation of a new bubble. In simulations of experiments #3 and #5, at lower heat fluxes into water, convective heat flux does not decrease to zero.

0.0E+00 2.0E+05 4.0E+05 6.0E+05 8.0E+05 1.0E+06 1.2E+06 1.4E+06 0.2 0.4 0.6 0.8 1

Axial distance (m) Heat flux, W /m 2 Q (conv) Q (quen) Q (evap)

Figure 12. Wall heat flux components: Qtotal=Qconvective+Qquenching+Qevaporation (simulation of experiment # 2). Comparisons between the calculated bulk void fraction distributions and the corresponding experimental measurements are presented in Fig. 13. Qualitatively, the experimental and numerical profiles of void fraction are in good agreement. Quantitatively, the discrepancy in different sections is up to ~30%. Good qualitative agreement was obtained for simulation of experiment #2. Simulations of experiments #3 and 5 overpredict and underpredict void fraction,

• respectively. This suggests that vaporization and condensation models should be improved to

handle a wide range of flow parameters in the same manner. Overall rather good agreement between numerical and experimental data was obtained in calculations for boiling and condensation processes. 3.4. Future Work At present, numerical analyses of experiments with a wider range of flow rate, density and wall heat flux than addressed in this paper are under way. The following is studied in these experiments: vapor-water flows with different void fractions, boiling and condensation in pipes with different water sub-cooling and wall heat release profiles, pressure drop along channels. The model development work also continues. The following major modeling efforts required for the simulation of coolant flow in BWR fuel bundle are planned:

• Implementation of the full flow regime map in the code allowing the identification and

treatment of all flow regimes typical for the coolant flow in BWR fuel bundles.

• Improved models for individual phenomena specific for various flow regimes, with special

attention given to the vapor/liquid mass, momentum, and energy transfers and to the fluid/wall interactions.

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 15/17

Experiment #2

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.2 0.4 0.6 0.8 1 1.2 1.4

Axial distance (m) void fraction

Expt STAR-CD

Experiment #3

0.05 0.1 0.15 0.2 0.25 0.3 0.2 0.4 0.6 0.8 1 1.2 1.4

Axial distance (m) void fraction Expt STAR-CD

Experiment #5

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.2 0.4 0.6 0.8 1 1.2 1.4

Axial distance (m) void fraction Expt STAR-CD

Figure 13. Average void fraction as a function of distance along the pipe.

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 16/17

• 4. CONCLUSIONS

Two-phase flow phenomena inside a BWR fuel bundle include coolant phase changes and multiple flow regimes which directly influence the coolant interaction with fuel assembly and, ultimately, the reactor performance. The CFD-BWR code is being developed as a specialized module built on the foundation of the commercial CFD code STAR-CD which provides general two-phase flow modeling capabilities. New models describing two-phase bubbly flow and heat transfer phenomena specific for BWRs have been developed and implemented in the CFD-BWR module. A set of experiments focused on two-phase flow and phase-change phenomena (boiling and condensation in upward flow in heated and adiabatic vertical straight pipes) has been identified for the validation of the CFD-BWR code and results of three experiment analyses have been

• presented. The close agreement between the computed results, the measured data and the

correlation results provides confidence in the accuracy of the initial models. ACKNOWLEDGMENTS The work was performed under partial financial support of ISTC, Project # 2601p, sponsored by the U.S. Department of Energy. REFERENCES

• 1. Proceedings of the OECD/CSNI Workshop on Transient Thermal-Hydraulic and Neutronic

Codes, Annapolis, Maryland, USA, November 5-8, 1996. NUREG/CP-0159, NEA/CSNI/R (97) 4.

• 2. Proceedings of the OECD/CSNI Workshop on Transient Thermal-Hydraulic and Neutronic

Codes, Barcelona, Spain, April 10-13, 2000.

• 3. STAR-CD Version 3.20 Methodology Manual, Chapter 13, CD-adapco, UK (2004).
• 4. S. Lo, “Modelling multiphase flow with an Eulerian approach”, von Karman Institute Lecture

Series - Industrial Two-Phase Flow CFD, May 23-27, von Karman Institute, Belgium (2005).

• 5. M. Lance, J. Bateille, “Turbulence in liquid and phase of a uniform bubbly air-water flow”,
• J. Fluid Mechanics, 222, pp 95-118 (1991).
• 6. N. Kurul, M.Z. Podowski, “Multidimensional effects in sub-cooled boiling”, Proc. 9'th Heat

Transfer Conference, Jerusalem (1990).

• 7. M. V. H. Del Valle, D.B.R. Kenning, “Subcooled flow boiling at high heat flux”, Int. J. Heat

Mass Transfer, 28, pp.1907-1920 (1985).

• 8. G. G. Bartolemei, V. M. Chanturiya, “Experimental study of true void fraction when boiling

subcooled water in vertical tubes”, Thermal Engineering, 14, pp.123-128 (1967).

• 9. R. T. Lahey, Jr., F. J. Moody, The Thermal-Hydraulics of a Boiling Water Nuclear Reactor,

2nd Edition, American Nuclear Society, La Grange Park, Illinois USA (1993).

• 10. G. G. Bartolemei, V. I. Gorburov, “Experimental study of vapour phase condensation in liquid

subcooled below saturation temperature”, Heat Production, No.12, p.58 (1969).

• 11. A. A. Avdeev, V. P. Pekhterev, “Vapour condensation in non-equilibrium bubbly flows”, High

Temperature, 1986, 24, No.6, p.1125.

• 12. G. G. Bartolemei, G. N. Batashova, V. G. Brantov, et.al. In: Teplomassoobmen-IV. Vol. 5.
• Minsk. ITMO AN BSSR Press, p.38 (1980).

American Nuclear Society Topical Meeting in Mathematics & Computations, Avignon, France, 2005 17/17