Title: Evaluation of RELAP5 reactor core modeling capability
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00100820/00001
 Material Information
Title: Evaluation of RELAP5 reactor core modeling capability
Physical Description: Book
Language: English
Creator: Roux, Vincent J.-P., 1977-
Publisher: University of Florida
Place of Publication: Gainesville Fla
Gainesville, Fla
Publication Date: 2001
Copyright Date: 2001
Subject: Nuclear and Radiological Engineering thesis, M.S   ( lcsh )
Dissertations, Academic -- Nuclear and Radiological Engineering -- UF   ( lcsh )
Genre: government publication (state, provincial, terriorial, dependent)   ( marcgt )
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )
Summary: ABSTRACT: RELAP5 is a one-dimensional reactor-system simulation code with additional cross-flow calculation capability to include two- and three-dimensional effects in light water nuclear reactor cores. The code is used to model the core, steam generator, and the balance of the Surry reactor, which is a three-loop Westinghouse Pressurized Water Reactor (PWR) system. A detailed RELAP5 model including full nodalization of the system is developed and implemented for this study. The RELAP5 Surry core model uses one or several parallel channels to compare and assess the performance of the cross-flow model. Several inlet flow rates and core power distributions are considered and modeled. Results of the analysis showed the significant contribution of cross-flow in overall temperature and flow distributions in the core. Results of the study also showed that the RELAP5 predictions of cross-flow, at least for single-phase cases, are not consistent with the theory. To evaluate the accuracy of RELAP5 cross-flow model, an industry standard Computational Fluid Dynamics code, FLUENT, is used to perform two- and three-dimensional calculations. Initial and boundary conditions for the RELAP5 model are used to develop a high-resolution FLUENT model for a pair of parallel reactor core channels. Two models were developed for FLUENT calculation of cross-flow. The first model is a simple tube with axisymmetric non-uniform inlet flow velocities. The second model included different power generation rates in the inner tube and the outer annulus.
Summary: ABSTRACT (cont.): Results of flow rate in the transverse direction were compared with RELAP 5 cross-flow calculations under similar geometry and boundary conditions. All FLUENT calculations were performed using simple orthogonal mesh, however with a finer mesh in the regions of interest. The k-epsilon turbulence model is used to calculate turbulent viscosity. The power density is adjusted to account for the uniform heat generation inside the core. Results of this study show the limitations and shortcomings of the RELAP5 to account for three-dimensional flow phenomena using the cross-flow model. The coarse mesh nodalization used in RELAP5 model limits the accuracy and resolution of cross-flow calculation. Results of RELAP5 calculation for different inlet flow conditions show the inadequacy of the simplified cross-flow model for predicting complicated two- and three-dimensional flow phenomena.
Summary: KEYWORDS: RELAP5, cross-flow, CFD, core modeling
Thesis: Thesis (M.S.)--University of Florida, 2001.
Bibliography: Includes bibliographical references (p. 86-87).
System Details: System requirements: World Wide Web browser and PDF reader.
System Details: Mode of access: World Wide Web.
Statement of Responsibility: by Vincent J.-P. Roux.
General Note: Title from first page of PDF file.
General Note: Document formatted into pages; contains xi, 88 p.; also contains graphics.
General Note: Vita.
 Record Information
Bibliographic ID: UF00100820
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: oclc - 49053197
alephbibnum - 002763559
notis - ANP1581


This item has the following downloads:

master ( PDF )

Full Text







Copyright 2001


Vincent J.-P. Roux

To Lucie


This research would not be completed without the contributions of several people.

I would like to thank Dr. Samim Anghaie for his support and his wise advice throughout

the project. His trust and guidance made my master's degree more than an academic

experience. I would like to express my thanks to everyone who took the time to help me;

especially Dr Gary Chen, whose experience and availability I value. Dr. Chris Allison

and Dr. Alexandre Ezzidi were also ready to help me every time I needed it. I would also

like to thank Dr. Edward T. Dugan and Dr. Raj at Mittal for their quality teaching and for

taking the time to be on my supervisory committee.

I would like to thank my parents and my sisters Helene and Cecile for their

understanding and their support. I would also like to express my thanks to the French

team. Their support and encouragement helped me to complete this research. I also thank

all my friends in the department and at INSPI who made me enjoy my stay at the

University of Florida. I am also grateful to Lynne, Denielle and Helene for reviewing my




A C K N O W L E D G M EN T S ......... ................................................................................ iv

LIST OF TABLES ................ ......... .................. vii

LIST OF FIGURES ................ ......... ...................viii

1 IN TR O D U C T IO N ................................................................... ........................ 1

2 SYSTEM MODELING WITH RELAP5 ............................. .......................... 4

A bout R E L A P 5 ...................................................................................... . . 4
Capabilities of RELAP 5 .................................................................... ..............5
Physical Models and Numerical Methods ............... .........................................6
Conservation equations ................................. ............................. ..... 6
N um erics ................................................................................................... 8
Cross-flow m odeling..................... ............. ... ....................................... 14
Comparing RELAP5 Core Model with and without Crossflow ............................. 14
Hydrodynam ic Com ponents ........................................................... .............. 15
C o re ............................................................ 1 5
Steam generator .............................................................................................. 18
Pressurizer ................................................................. .. .... ........ 22
Prim ary cooling system tubing .................................................................... 24
H eat Structures ............................................................................................ ........ 26
H eat structures in the core ....................................................................... .......... 26
H eat structure in the loop ....................................................................... ........... 28
Heat structures in the steam generator ............ ........ ................28
Steady State and Initialization Process ...................... ..................... 29

RESOLUTION COMPUTATIONAL FLUID DYNAMICS ..................................... 33

A bout FLU EN T ............................................................................................... .......33
G en eral O v erv iew ............................................................................................. 3 4
Basic Physical Models ................................................................................................35
N u m erical S o lv ers ............................................................................................ 3 7
Numerical schemes and discretization........................ ..... ......... 37
Boundary conditions .............................................................................. 39
U ser defin ed fu n action s .................................................................................. 4 1


Turbulence Modeling in FLUENT............. ........................ .............. 42
One-equation Spalart-Allmaras model ............. ..... .....................................46
Tw o-equation k-epsilon m odel ................................... ............ ....... .............. 46
Defining a CFD Method for Evaluating RELAP5 Cross-Flow Model........................48
F L U E N T M o d el .............................................................................................. 4 9
GAM BIT grid generation........................................................... .......... .... 50
Boundary and operating conditions ......... ............................... ............. 51
R E L A P 5 E equivalent M odel ..................................................................... ....... 53

M O D EL ............................................................... ..... ....... ........ 55

Surry Reactor Calculation for a Steady State and a Simulated Power Surge Transient 55
Steady State ........ ............................. ... .. ......... ...... ........... 55
Transient Case: Simulation of a Power Surge Transient .................................... 60
Cross-flow Between Two Adjacent Channels: FLUENT and RELAP5 Results..........64
Simulation in a Large-Diam eter Vertical Pipe ......... ..................... .................. 65
If no power is added to the fluid .......................... .... ................... 65
If power is added to the fluid......... ......................................... .............. 70
Simulation in a Small-Diameter Pipe ......................... ..............73
If no power is added to the fluid .......................... .... ................... 73
If power is added to the fluid......... ......................................... .............. 76
Steam flow ........................ ... ................... .............. 76
Conclusions and Recommendations..... ...................... ................78



GENERATOR ................................................................... .......... ........ 83

C RELAP5 INPUTS FOR THE TWO-CHANNELS SYSTEM............................... 84

LIST OF REFERENCES... ...................................................... ...... .............. 86

BIOGRAPHICAL SKETCH .......................................................... ........... ..... 88


Table Page

2-1 Fuel rod pow er distribution ........................ ............. .................. .............. 27

4-1 Comparison of axial junction mass flow for a core with and without cross-flow.....56

A-i Hydrodynamic components used to build the reactor model.................. ..............80

A-2 Hydrodynamic components in the steam generator................................................. 81

A-3 Hydrodynamic components of the primary loop................................. ........... 81

A-4 Hydrodynamic components of the pressurizer .................................................... 82

B-l Heat structures in the core and the steam generator. ............. ...........................83


Figure Page

2-1 Difference equation nodalization schematic. .................................. .... ............... 9

2-2 RELAP5 nodalization diagram of the Surry reactor core............ ................ 17

2-3 RELAP5 nodalization diagram of a U-tube steam generator .............................20

2-4 RELAP5 nodalization diagram of the pressurizer on the hot leg ..........................23

2-5 RELAP5 nodalization diagram of the primary cooling system tubing.....................25

2-6 Bundle matrix for the first channel: 1 is a fuel rod 2 is a control rod................ 28

2-7 Evolution of the temperature distribution during the steady state run........ ........ 30

2-8 Axial coolant temperature distribution the steady state............... ................31

2-9 Radial fuel temperature profile at the sixth vertical node (height=2.012 m)............32

3-1 Mixing-length distribution in wall boundary layers.................... ...............45

3-2 Basic geometry studied with RELAP5 and FLUENT...........................................49

3-3 FLUENT calculation domain of the system...........................................................49

3-4 Nodalization diagram of the equivalent RELAP5 model., ..................................... 54

4-1 Cross-flow magnitude versus height in the core for a total power of 2443MW.......57

4-2 Cross-flow mass flux versus core cell height..................................................... 58

4-3 Ratio of cross-flow to axial flow for different power levels and for different
height in the core.................................... ........................................ 59

4-4 Power transient simulated in the transient test............................................. 61

4-5 Cross-flow magnitude versus height in the core for a total power of 3500MW.......61

4-6 Cross-flow to inlet flow ratio for a 1000MWth case. Cross-flow is from Channel 1
to 2 ..................................... .................................... ........... 62

4-7 Cross-flow to inlet flow ratio for a 2443MWth case. Cross-flow is from Channel 1
to 2 .................................... ................... ................ ........... 63

4-8 Two-sub channel pipe studied with FLUENT and RELAP5 ...............................64

4-9 Inlet velocity profiles used in cross-flow simulation ...............................................65

4-10 Axial evolution of the axial velocity profile. z is the height in meters................... 66

4-11 Radial (V-velocity) and axial (U-velocity) velocity distribution at the entrance of
the tube for Vi=6.0 m/s and Vo=2.0 m/s...........................................................67

4-12 Radial (V-velocity) and axial (U-velocity) velocity distribution at the entrance of
the tube for Vi=4.0 m/s and Vo=2.0 m/s...........................................................67

4-13 Comparison of cross-flow prediction between FLUENT and RELAP for two
sim ulations............. ................. .................................................. 68

4-14 Evolution of the axial velocity for the two RELAP5 axial channels along the pipe,
for the two inlet boundary conditions (a) and (b). ....................... ...............69

4-15 Temperature distribution over the total length of a 20.0 m pipe heated in the
central channel at 66.1 M W /m3. .............................................. .................. 70

4-16 Radial velocity for a constant inlet velocity and a heated central channel. ..............71

4-17 Cross-flow velocity at the artificial interface between the two channels. ...............72

4-18 Axial velocity in Channel 1 (heated) and Channel 2 (non-heated) ..........................72

4-19 Radial velocity in a small size tube in an adiabatic flow. Vi=4.0 m/s and
V o=2.0 m /s............. .. ................. .............. .. ......... .... ........ 74

4-20 Axial velocity in a small size tube in an adiabatic flow. Vi=4.0 m/s and
V o=2.0 m /s. ............ .. ................. .............. .. ......... .... ............ 74

4-21 Cross-flow velocity versus distance to diameter ratio for large and small diameter
tubes. ............. .............. .. .. ........ .. ................. ........... 75

4-22 Cross-flow velocity for unheated small diameter channels (RELAP5 calculations).75

4-23 Cross-flow velocity between heated small diameter channels..............................77

4-24 Axial velocity distribution for superheated steam. ...........................................77

Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science



Vincent J.-P. Roux

August 2001

Chairman: Professor Samim Anghaie
Major Department: Nuclear and Radiological Engineering

RELAP5 is a one-dimensional reactor-system simulation code with additional

cross-flow calculation capability to include two- and three-dimensional effects in light

water nuclear reactor cores. The code is used to model the core, steam generator, and the

balance of the Surry reactor, which is a three-loop Westinghouse Pressurized Water

Reactor (PWR) system. A detailed RELAP5 model including full nodalization of the

system is developed and implemented for this study. The RELAP5 Surry core model uses

one or several parallel channels to compare and assess the performance of the cross-flow

model. Several inlet flow rates and core power distributions are considered and modeled.

Results of the analysis showed the significant contribution of cross-flow in overall

temperature and flow distributions in the core. Results of the study also showed that the

RELAP5 predictions of cross-flow, at least for single-phase cases, are not consistent with

the theory.

To evaluate the accuracy of RELAP5 cross-flow model, an industry standard

Computational Fluid Dynamics code, FLUENT, is used to perform two- and three-

dimensional calculations. Initial and boundary conditions for the RELAP5 model are

used to develop a high-resolution FLUENT model for a pair of parallel reactor core

channels. Two models were developed for FLUENT calculation of cross-flow. The first

model is a simple tube with axisymmetric non-uniform inlet flow velocities. The second

model included different power generation rates in the inner tube and the outer annulus.

Results of flow rate in the transverse direction were compared with RELAP 5 cross-flow

calculations under similar geometry and boundary conditions. All FLUENT calculations

were performed using simple orthogonal mesh, however with a finer mesh in the regions

of interest. The k-epsilon turbulence model is used to calculate turbulent viscosity. The

power density is adjusted to account for the uniform heat generation inside the core.

Results of this study show the limitations and shortcomings of the RELAP5 to

account for three-dimensional flow phenomena using the cross-flow model. The coarse

mesh nodalization used in RELAP5 model limits the accuracy and resolution of cross-

flow calculation. Results of RELAP5 calculation for different inlet flow conditions show

the inadequacy of the simplified cross-flow model for predicting complicated two- and

three-dimensional flow phenomena.


Research and development in the nuclear field is always active to create a new

generation of safer and safer reactors. For simplicity of operating procedures and also to

optimize the efficiency of the reactors, people try to understand more and more in detail

all the phenomena that can occur in a reactor core. Researchers develop models to

represent and calculate the behavior of the core in any situation to be able to know

exactly the state of the reactor, should a problem occur; or even during normal operating

conditions. A nuclear reactor is a complex entity that requires knowledge in neutronics,

thermalhydraulics, chemistry, but also automated systems, informatics to develop the

control system of the reactor. This study focuses on the thermalhydraulic characteristics

and behavior of the reactor core. Many codes have been developed to predict temperature

distribution in the core and pressure field intensity. An accurate mapping of the flow

paths is useful to calculate heat removal capabilities, stresses inside the reactor, efficiency

of the reactor, and safety features. Also from a neutronics point of view, flow path can be

important, especially for two-phase flows. This will help in the design of new reactors. It

is important to make sure that no component in the core reaches its melting point. Among

all these tools, RELAP5 is one of the most popular and most widely used codes to do

calculation of in-core parameters and to simulate reactor transients in Light Water

Reactors (LWR). It allows for the modeling of the entire reactor and performs

calculations with an excellent accuracy. The entire nuclear community recognizes its


However, numerical codes like RELAP5 need closure relations to solve the

system of governing equations. Physical models adapted to a numerical treatment need to

be used. These models, based on theoretical relations or empirical observations, are

continuously improving. This study focuses on the cross-flow phenomenon. Crossflow is

a well-known phenomenon, and extensive literature is available to show experimental

observations and attempts to model this specific type of flow. By definition, the term

"cross-flow" designates the fluid moving in a direction that is not the one of the main

stream. Both streams do not necessarily have a direct interaction with each other. Cross-

flow can be found in heat exchangers, for example. If the system is designed so that the

cooling fluid is flowing around the tube bundle of the primary system in a normal

direction, this is cross-flow cooling. However, the designation cross-flow is more

commonly adapted to the flow of a jet mixing in another flow flowing in another

direction. But these are cross-flows generated by geometry conditions, as analyzed by

Wang (1). Cross-flow has been extensively studied considering jet mixers. This particular

case concerns injection of fluid through a jet inside the main flow. Mixing velocity and

turbulences generated are then studied (2). Computational Fluid Dynamics (CFD) tools

allow calculating flow patterns, fluid velocities and turbulence parameters. By comparing

experimental observations, turbulence models have been validated. A lot of work has

been done to simulate flow path and the steady state velocity field of a cross-flow jet

mixing in an axial stream. In a nuclear reactor core, cross-flow is due to grids and mixing

devices, but also to inhomogeneous properties of the fluid. This phenomenon has a low

magnitude in a typical Pressurized Water Reactor (PWR) because the pumps, which

impose a mass flux in the loop, drive the flow axially.

However, in case of natural circulation conditions, there is not the same strength

anymore on the fluid, and then the flow is more sensitive to heterogeneity in the fluid

properties. Sometimes, the problem can be much more difficult to handle. If the interest

is in two-phase flow, then the modeling must be done very carefully because the

difference of density between the two phases present is so large that buoyancy cannot be

neglected and stress interaction between the two phases has to be carefully modeled. In a

Boiling Water Reactor (BWR), cross-flow across the core is important. Cross-flow in a

nuclear reactor core occurs between the different channels (or even sub-channels). If the

peaking factors in the two adjacent channels, or assemblies are different, or if there is a

grid mixer that makes the flow diverging from its path, it is important to understand what

the effects of the cross-flow can be on the system. As illustrated by Kassera (3), cross-

flow can induce tube vibration, which would be a serious concern in a nuclear reactor or

in a steam generator, where the flow in the secondary loop is a mixture of two phases.

The purpose of this study is to analyze the efficiency of cross-flow modeling in

RELAP5. The CFD calculation code FLUENT was used for this purpose. It has been

useful to work on a simple system, and develop a model adapted for each of the codes.

These models enhanced the cross-flow phenomenon, and a comparison of the results was

performed. Finally, the calculations showed that there was a difference between RELAP5

and FLUENT cross-flow calculations. This thesis presents an actual simulation of a

reactor, enhancing occurrence of cross-flow inside the core, and discusses the modeling

of cross-flow in RELAP5, based on high-resolution CFD simulations.


In the modeling of cross-flow phenomenon, the first tool used was the

thermalhydraulics code RELAP5. This code was developed to simulate a nuclear reactor

coolant system during accidents as well as for normal operating conditions of the reactor.

RELAP5 was used to simulate cross-flow phenomenon in a nuclear reactor core. First we

present a brief description of RELAP5 capabilities. The model used for the calculation is

based on the Surry reactor. The RELAP5 nodalization diagram of this reactor is

explained in details in this chapter. Finally, the results of our simulation, emphasizing

cross-flow behavior in the reactor core is presented.

The code RELAP5 was developed to simulate the behavior of a whole cooling

system of the reactor. Even with a detailed description of the core's thermalhydraulic

behavior, it is impossible to have a very high-resolution representation of the core. The

whole reactor system was run to achieve a steady state and to achieve initial conditions

for the CFD simulation.

About RELAP5

RELAP5 was used to simulate the whole core behavior during steady state and

transient conditions. This code was developed for the U.S. Nuclear Regulatory

Commission (US NRC) for transient analysis of light water reactors (LWR). The

RELAP5 version used was provided by Dr. Chris Allison, from Innovative Software

System, Idaho Falls, Idaho (4). This version of the code is equipped with the severe core

damage package SCDAP. The whole SCDAP/RELAP5/MOD3.2 package has a plot

option that can plot evolution of any variable versus time, as requested. In the core

model, SCDAP components will be used, however, no accident scenario will be


Capabilities of RELAP5

A specific use of RELAP5 is the simulation of transients in LWR systems.

RELAP5 could, however, be applied for other kind of non-nuclear related

thermalhydraulic system. This part shows the capabilities of the code RELAP5. The

system's thermalhydraulic models are fully coupled at each time step. Hydrodynamics,

heat conduction, reactor kinetics, control systems, trip logic and special component

models, such as pumps or turbines will govern a RELAP5 case. Hydrodynamics is

modeled by non-equilibrium, six-equation and two-fluid model. RELAP5 allows

choosing an arbitrary number of volumes, junctions and surfaces. RELAP5 is a one-

dimensional code; however, it can handle multiple-dimensional nodalization by using

cross-flow junctions. Orthogonal direction then has to be specified. Heat transfer in

RELAP5 is handled by ID heat conduction theory for rectangular, cylindrical, or

spherical geometry. An extensive heat transfer package is included in RELAP5 with the

consideration of the effect of non-condensable gases. RELAP5 also takes radiative heat

transfers into account. A reactor kinetics option is available in RELAP5. This task is

performed by point kinetics theory or space dependent kinetics.

A control system in RELAP5 is performed by using different tools. The basic

arithmetic operations are used to build the control system. It is also possible to use

differentiation and integration of variables to have a well-defined control system. Other

direct control tools, such as standard or tabular functions can be used. To have an

automated control system, trip logic operates the system. The variable trips are defined

using arithmetic comparisons like "greater or equal" or "less than"... Logical trips were

built with boolean operators "and," "or" and "xor." For instance, more developed logical

trips are built in the code to control valve opening or pump motor velocity.

Physical Models and Numerical Methods

RELAP5 is a one-dimension thermalhydraulic code adapted for transient, two-

fluid model. Thus, the code solves six equations of conservation as it simulates the two

phases, liquid water and steam. Conservation equations of mass, momentum and energy

are solved by RELAP5 for each phase. These equations are inter-related with each other.

The eight primary independent variables considered in these equations are pressure (P),

phasic specific internal energies (Ug and Uf), void fraction (a ), phasic velocities (vg and

vf), noncondensable quality (Xn) and boron density ( pb). Independent variables are time

(t) and distance (x). There are also secondary dependent variables such as phasic densities

( pg and pf), phasic temperatures (Tg and Tf), saturation temperature (Ts) and

noncondensable mass fraction in noncondensable gas phase (Xni). The nomenclature used

is the same as in the RELAP5 Code Manual (5).

Conservation equations

The following are the conservation equations solved in RELAP5. The momentum

equations are not modified when noncondensable species are present, however energy


equations are slightly modified and consider heat transfer at the noncondensable gas-

liquid interface. These are the terms in {-} in the equations.

1. Vapor mass conservation: (, p)+ 1 a pa vA)=Fg
at A ax

2. Liquid mass conservation: -(a p)+ -- a (apvA)= Ff
at A Ox

3. Vapor momentum conservation:

ap9,A + -apgA v
Ot 2 Ox
AP + ag pBxA- A(agpA)FWG(vg)+ FA(, -vg )

(ag pA)FIG(v, -vf)-CaapmA[ v) +Vf -

4. Liquid momentum conservation:

avf 1 av2
a, pA At-+-a A fpfA
Qt 2 8x
-afA +afpfBA- (afpfA)FWF(vf)- A(v -vf)
a a )g

-(apA)FIF(vf -v,)-CaafpmA A v -v) +v V
at 9 a x f

5. Vapor energy conservation:

a (agpgUg)+-l (agpgUgVgA)=

P- -au a(vA)+Q, +Qg + ghg + +Fg +DISSg -QgfJ
at Aa

6. Liquid energy conservation:

(a p U, )+ (apUvA)=
at A+
aP f a(afVfA)+Q f Qf Fh* 17 h' DJ
at Aa f

7. Noncondensable component mass conservation:

a (ogpgXn)r1 ( gpgXnvgA)= 0

8. Boron concentration in the liquid field:

Pb +l (a fbv f A)= 0
at A 8x


Spatial discretization of the set of differential equations is performed by

integrating the differential equations over the cell volume for the volume quantities (that

define scalar properties like pressure, energies and void fraction), and between cell

centers for junction quantities (representing vector quantities like velocities). Figure 2-1

shows the typical spatial noding resulting from this definition (5). This approach results

in a numerical scheme having a staggered spatial mesh. Numerically there should not be

an ill-posed problem; RELAP5 uses different techniques to get rid of instability problems

for mesh sizes of practical interest. This is done by avoiding explicit treatment of time

advancement. Two possible time advancement schemes are available in RELAP5. Each

of them is more or less adapted to a specific case. For instance, it will be more adapted to

use a nearly implicit scheme for a steady-state or self-initialization case problem where

the time step is limited by the Courant limit, but also for slow phases of a transient

problem. Basically, a semi implicit time advancement scheme is adapted for transient

runs. The principle of the semi-implicit scheme is to replace the system of partial

differential equations with a system of finite-difference equations partially implicit in

time. For the semi-implicit solution strategy, seven simultaneous equations are to be

solved: two per junction and five per volume. The five density and energy variables are

expressed as differences (in time) and are ordered, AXn,i=Xn,im+ -Xn,im, Aug,i=Ug,i -

ug,im, Auf Uiml A, = m, a, and A =im+-Pim, m is the time index.
Ug,i uf, i-uf, -u i g, g1 -{z g,1

Junction unknowns are velocities vg,jm+ and vfjm+for the junction j.

Vector node
or junction
Vg, Vf Mass and energy control
Scalar node volume or cell
P, og, Ug, Uf


j-1 J j+1
Momentum control
volume or cell
Figure 2-1: Difference equation nodalization schematic.

For each junction, the only unknowns involved are the velocities and the pressure

of the volumes connected by the junction. The tilde (-) sign is used to specify that the

term is a matrix type term. All terms in these matrices are calculated using "old time"

values. Components of the vector v"+1 are the gaseous and liquid phase velocities at

junction j at time m+1.

Btm 1=(pm +AP) (pm+ AP,)+r,

For each junction, we have a two-equation system easily solved in term of the pressure

change, however it is not yet possible to get a value for the new pressure, nor the new


Vm1 B (A1 (A +1 1)+B r1
v h = B k + A/m+l )+1fi b

The five other volume equations can be written in an order-five matrix form. Fortunately,

the unknowns for each volume equation involve only the unknown quantities for this

given volume and the velocities for the junctions attached to the volume. And so, for each

volume i,

4A, =Ab+IV(91+/;1 )

with A = Auf, and A is the order five matrix of the system. The detail of the terms in
A ,

these equations can be obtained in the RELAP5 user manual (5). Thus, the five equations

allow the volume quantity for each volume to be written in terms of the velocities in the

attached junctions. Then from the previous comment, at this point only the Ag solution

is needed. The matrix Ai is then factorized into lower and upper triangular matrix

(Ai=LU), and from these, the bottom row of the inverse A'i is computed. It can be shown

that this particular row represents an equation involving only the pressure in the volume

and the velocities in the attached junction:

A 51 =A,, + g,Jg, +A521, 2,z + i.g2,jVg,+ +A, bz, +
1 b m \ j3 1 j V7i^

A51, b4,, +- -(4, 1'J + f-4, jV+) 5+A 1 5, +1 (s, v"'J + f, v"f )

This equation involves only the pressure in the volume and the velocities in the attached

junctions. Substituting the velocity expressions previously obtained from the junction

equations into this expression, the resulting expression involves the AP for the volume i

and the AP's for each volume connected to this volume by the attached junctions. So

combining the AP equation for each volume leads to a system of simultaneous

equations, one per volume. And,

CAP = r

For a strictly one-dimensional pipe, the matrix C is tridiagonal. A typical reactor system

with loops and branches would not have a regular pattern, but a proper ordering of the

equations would show sub matrices of tridiagonal and two-dimensional form. This

system is solved and the AP 's obtained from the pressure equations are directly

substituted into the two velocity equations for each junction, thus obtaining the velocities.

The AP 's and velocities are substituted into the remaining four equations (out of the

original five) for each volume. These four equations are actually solved for A Xn,i,

Au g,i, Au f,i and Aa *, which are the unknown values at a provisional "new time". The

actual "new time" variables AXn,i, Aug,i, Aufi and Ag,, are determined using the

unexpanded difference equations (5). Note that instead of solving a set of two equations

per junction plus five equations per volume, the semi-implicit technique allows solving

an order two system for each junction, an order five system for each volume, and one

system of one equation per volume. This was possible by using new or old time values

for hydrodynamic unknowns. This technique is characterized by using only old time

values for convection quantities. The semi-implicit method limits the time step to the

Courant limit. The Courant limit is a parameter equal to the ratio of the length of the cell

to the velocity of the fluid in the cell. This numerical limitation has a physical meaning.

By limiting the time step to the Courant limit, it does not allow fluid to travel a path

longer than the cell itself during the time step. The nearly implicit method is not limited

by Courant limit. The basic idea of the nearly implicit scheme is to split the equations

into fractional steps based upon physical phenomena.

This scheme consists of two major steps to avoid solving the complete set of two

equations per junction and five equations per volume. The first solves all conservation

equation in the expanded form, as for the semi-implicit method. It results in one equation

per volume involving the Ag for the volume i and the velocities of the junctions attached

to this volume. Convective terms in the momentum conservation equations (junction

equations) are treated implicitly, meaning that the scheme uses new time values. The

resultant equations involve the velocities in junction j and the AP 's from the volume k

and 1 connected by junction j (Figure 2-1) as in the semi-implicit technique, but also

involve the velocities from each junction connected to volumes k and 1. The simple

solution of the order two matrices for the velocities in terms of the pressure changes is no

longer possible and the substitution of the velocity expressions into pressure equation is

likewise not possible. Instead, substitute the pressure equation into the velocity equations.

The resultant two equations for each junction involve vgj and vfj and the velocities for

each junction connected to the volumes k and 1. Combining these two equations per

junction leads to a set of simultaneous equations involving only the velocities. The shapes

of the matrix for the semi-implicit technique and this stage of the nearly implicit step

have similarities. The junction equations are solved using the same sparse matrix solver

as used in the semi implicit technique. The velocities are substituted into the five

equations for each volume, and each order five system is solved for the volume quantities

AX, ,, Au,, Au,, ,,,A, and APi. At this point, we do not obtain "new time" values. The

tilde (-) symbol on the top of the first four quantities indicates that these are not the final

value for these quantities.

The second step consists of solving the unexpanded form of the mass and energy

equation. Using the final values for the velocities and pressure changes, and the

intermediate values for the noncondensable mass fractions, internal energies, and void

fractions, the phasic conservation of mass and conservation of energy equations are

written with the convective quantities now at new time values. Interphase quantities such

as interphase mass transfer are evaluated using the linearized equations with intermediate

results. For a volume, each phasic equation involves only one type of unknown. Vapor

noncondensable equation involves only (a pX )m+l, vapor mass equation involves only

(a p, )m+1, liquid mass equation involves only (a pf )m+1, vapor energy equation

involves only (a, pu,)m+l and liquid energy equation involves only (af puf )m+1. An

equation for a volume involves the convected terms for that volume and convected terms

from the volumes connected by the attached junctions. As five sets of simultaneous

equations result from gathering the same type of equation for every volume, the shape of

these five matrices are the same and have the same shape as the pressure matrix C in the

semi-implicit scheme. The terms Xnm+1, Um, ufm+1 and a7 are obtained by dividing

between each other's previous five terms. p7'1 is obtained from linearized equation of

state using Pm+1 and ufm+1

A detailed description of the two methods is extensively presented in the first

volume of the RELAP5 code manual (5).

Cross-flow modeling

As RELAP5 is a one-dimension code, a specific model is used to calculate cross-

flow, which by definition, is occurring in the normal direction to the main direction of the

flow. For small cross-flow between reactor core channels, simplified cross-flow

momentum equations are used to couple two adjacent channels linked with a cross-flow

junction. In the momentum equation, the transverse momentum advection terms are

neglected, meaning that there is no transport of x-direction momentum due to flow in the

transverse direction. Cross-flow areas can be entered or calculated by the code. If the

code has to calculate the junction area, it considers the geometry of a cylindrical pipe.

Detailed governing equations for modeling of cross-flow phenomenon are presented in

the user's manual (5).

Comparing RELAP5 Core Model with and without Crossflow

Calculations were performed with the model of the Surry reactor. The Surry

Nuclear Plant is located in Norfolk, Virginia. It has two Westinghouse 3-Loop PWR

reactors. Each unit is rated 823 MWe (2443 MWth). Surry is owned and operated by

Virginia Power. Each loop is composed of a steam generator, a pump and a tubing

system. One of them has a pressurizer on the hot leg. An accumulator is also present on

every loop on the cold leg, for the emergency cooling system. The following description

presents these elements: the core, the steam generator, the pressurizer and the tubing

related to the primary loop. For each element, a picture of the nodalization diagram and

the initial conditions for temperatures and flow rates in the cells and junctions are given.

Also this paper describes the heat structures that are used to model the heat transfer

characteristics inside the core and along the loops and also are used to model the power

generation in the core. A detailed description of the input requirements can be found in

the Appendix A of the SCDAP/RELAP5/MOD3.2 user's guide and input manual (6).

Hydrodynamic Components

The reactor model is mainly composed of hydrodynamic components that

basically represent the parts of the reactor where coolant passes through and heat

structures that represent solid parts of the reactor (there is no flow in the heat structures)

where also heat can be generated or conducted to other regions of the reactor or the

environment. The following is a presentation of the different parts of the reactor modeled

in this problem.


Figure 2-2 shows the nodalization diagram of the core. This is a five-channel core

model. The core is composed of a downcommer (104) coming to the lower head (106)

and the lower plenum (108). The fueled part of the core is composed of five channels

(111, 112, 113, 114 and 115) and a bypass region (118). The upper plenum of the core is

only a 3-channel model ([151-152-153-154], [161-162-163-164] and [171-172-173-

174]). The control assembly housing is represented in three parts (181, 182 and 183). The

upper head (190) is covering the core. Upper annulus (100) and inlet annulus (102) are

also represented in this model. This is a quite detailed nodalization scheme. It can be

useful to describe and analyze cross-flow phenomenon. Inlet annulus (component 102) is

the hydrodynamic volume where cold coolant enters the core. Junctions with upper

annulus (100) and upper plenum (172) are considered, but surface junctions are much

smaller than with the downcommer (104), meaning that most of the flow will go down

the downcommer. Also initial junction coolant velocities are set in consequence. The

downcommer is represented by a pipe and divided into seven cells. Initial conditions for

these cells are pressure and temperature. Pressure is calculated to take into account the

gravity effects (from 2286.0 psi for the upper cell to 2290.5 psi for the lower cell).

Temperature is essentially constant throughout the downcommer and about equal to

543.0 F. The lower head (106) is a branch-type hydrodynamic component. Pressure and

temperature are given in this branch and equal to 2294.4 psi and 543.02 F. Initial liquid

and vapor velocities in the junctions are given. The lower plenum is described the same

way as the lower head. Inlet loss coefficients are set to 8.0 for each channel (bypass and

core channels). Inlet liquid velocity is 13.4 ft/s in the core channels and 8.4 ft/s in the

bypass region. The core channels are modeled with pipe components. The only difference

between them is the volume flow area of the cells. Initial pressure is 2275.0 psi in all the

cells and 577.0 F. Fluid velocities are given at the junctions (14 ft/s). Ten cells of

identical size are used to model the length of the channel (12 feet long heated rods). Heat

structures associated to these hydrodynamic volumes will be described later. As

represented on the diagram, cells of different channels are interconnected with each other

using single cross-flow junctions:

Components 140-149 connect the center core Channel 1 and center core

Channel 2.

Components 120-129 connect the center core Channel 2 and the middle

core Channel 3.


Components 80-89 connect the middle core Channel 3 and the middle core

Channel 4.

Components 130-139 connect the middle core Channel 4 and the outer

core Channel.

174 164 154
100 183 ~ 182 s
173 [- IJ 163 Y 153
from cold I to hot legs
legs 102 : 172 162 152

171 161 151

113 115 112 114 111



Figure 2-2: RELAP5 nodalization diagram of the Surry reactor core.
Figure 2-2: RELAP5 nodalization diagram of the Surry reactor core.

Initial velocities for all junctions between the channels are set to zero. There is no

junction between the bypass region and the center channels. In the case of calculation

without cross-flow, the four sets of cross-flow junctions are not included. The flow is not

able to pass from one channel to the other. The bypass region is modeled by a five cell-

pipe. Pressure in assumed to be 2276.0 psi and temperature 543.2 F. Initial velocity at

the four inner junctions is 0.8124 ft/s. The upper plenum is composed of branch

components, which allow cross-flow in the upper part of the core and drive the flow of

coolant to the outlet of the core hydrodynamicc component 172). Junction between the

heated channels and the upper-plenum (volume 151, 161 or 171) induces energy flow

losses. The loss coefficient is equal to 8.0 for the outlet effect of the channels. Energy

losses are also considered in the transversal direction (coefficient equal to 5.0 for the

junctions 151-161 and 152-162 and 7.0 for the junctions 161-171 and 162-172). The

inverted hat of the upper head has been represented using branch components (153-154,

163-164 and 173-174). In each of these components, initial pressure is equal to 2250.0

psi and temperature to 606.0 F. Housing for the control assembly is represented by

single volume components (one for each upper head sub-channel, 181, 182 and 183).

Pressure is set to 2250.0 psi and temperature to 550.0 F. To cover all this and the core

component, the Upper Head (190) is modeled using a branch component in order to take

into account all the junctions with control assembly housing volumes. Initial pressure in

the upper head is 2250.0 psi and temperature is 543.0 F. It is considered that there is no

loss between the upper head and the upper plenum and initial liquid velocity at the

junction upper head/upper plenum is 1.0 ft/s. The flow is going down in the reactor as

some flow is coming up from the upper annulus (100). Actually, this flow is really low.

Steam generator

Figure 2-3 represents the nodalization of the secondary side of one of the steam

generators used in each of the loops. The primary side, represented only by U-tubes is

described with the primary side tubing system. Radioactive water from the primary side

enters the steam generator at the bottom and flows through small diameter reversed U-

tubes. This coolant looses its heat, which is transferred to the water of the secondary side

coming from the steam generator feedwater and flowing outside the tubes. This system

allows this non-radioactive water to get turned into steam and then, moist steam is

produced and goes through the steam separator which provides dry steam to the steam

outlet (this steam is then sent to drive the turbines) and remaining water is sent directly to

the downcommer of the steam generator to be reheated. Finally, the superheated steam is

distributed through valves present at the outlet of the steam generator dome. The steam

generator feedwater is represented by a time dependent volume (266) where the initial

conditions for water supply are 785.0 psi and 430.0 F. The downcommer of the

secondary side of the steam generator is comprised of three components, the

hydrodynamic volumes (270, 272 and 274). The inlet volume is the branch 272 which

receive water from the inlet feedwater, the upper volume 270 and the separator 278.

Water coming from the separator is the result of the separation of steam and moisture to

provide only dry steam to the turbines. Initial properties for the water in the volume 272

are a pressure of 832.13 psi, and specific internal energies of 495.85 Btu/lb for the liquid

phase and 1114.42 Btu/lb for the vapor phase. A pipe divided into four cells models the

main downcommer part. For these volumes, initial conditions are given by Pressure,

liquid specific internal energy and vapor specific internal energy. Pressure calculation

includes a gravity effect. Change in internal energy (for liquid or vapor) is not significant

and is about 495.8 Btu/lb for the liquid and 1114.3 Btu/lb for the vapor phase. The

pressure is 834.24 psi at the top of the pipe and 839.92 psi at the bottom of the pipe. The

vertical pipe used to model the cavity surrounding the inverted U-tubes where water will

be boiled (276) is divided into five parts. For each of them, initial conditions on Pressure,

specific internal energy for liquid and vapor phase are given. For this region, it is

necessary to specify the hydraulic diameter, taking into account the number of tubes. This

has been calculated and the hydraulic diameter of the component is 0.14 ft.


280 safety valve

278 286 288 290
Separatoi i

266 272 steam
L- outlet
Steam Generator

274 276

d Dwncomrr er

Figure 2-3: RELAP5 nodalization diagram of a U-tube steam generator.

Also energy loss coefficients are specified and equal to 16.9 at the first junction

(which can represent the bottom of the bundle) and 6.0 for the middle region (the second

and third junctions). No loss is considered in the last junction. The swirl vane moisture

separator has a built-in model in the code, and this model is used to model the steam

separator (278). The purpose of this component is to separate steam from liquid water

and send the latter back to the downcommer. Drying of the steam is performed in two

stages so that an additional dryer is added (branch 280 in this model). The dryer is also

connected to the downcommer to get rid of the excess liquid water. At the top of the unit,

the steam generator Dome is represented by a single volume component. Initial

conditions for the fluid in the dome are a Pressure of 830.16 psi, specific internal energies

of 511.85 Btu/lb for the liquid phase and 1114.49 Btu/lb for the vapor phase. The vapor

void fraction is very close to 1.0. The steam line connecting the steam dome to the

different valves managing the steam distribution and the safety of the system is

represented by the pipe 284. Three valves are present in the system. The purpose of each

of them is well defined.

The main valve (285) controlling the steam flow through the main steam

line is represented with a servo-valve (model built in the code. The valve

flow area is calculated by a control system). The initial velocities are set to

102.2 ft/s for the liquid phase (even if the quantity of liquid is much lower

than the vapor one) and 115.6 ft/s for the vapor phase.

The second valve (287) is called the porv "power operated relief valve". A

"motor valve" type valve models this valve. This valve opens if pressure

in the dome steam line 284 gets greater than 1050.0 psi and closes when

pressure in 284 has become less than 1000.0 psi. This valve is initially

closed and the valve change rate is specified in the inputs. Initial velocity

has to be set to zero as the valve is initially closed.

The third one is a safety valve (289), which is also a motor valve. It opens

when the pressure in 284 is higher than 1184.0 psi and closes when

pressure in 284 drops down below 1092.0 psi. Again, this valve is initially

closed, but has a much higher junction area (which allows the pressure to

decrease faster) and the valve change rate is about 40 times larger than the

porv valve, which means faster opening ability. In term of rate of change

of normalized valve area, the porv opens at a rate 0.556 s-1 whereas the

safety valve opens at a rate 20.0 s-1. Here again, the initial velocity is set to

zero as the safety valve is closed in normal conditions of operation.

The opening of the valves is commanded by trips defined in the input file. With

these trips, it is possible to build control variables (using basically logic statements to

elaborate different trips).


The Pressurizer is used to maintain a constant pressure in the primary cooling

system. It is composed of a large tank. It is about half full of liquid water and the other

half contains steam. The whole tank was designed to remain at constant pressure, and

because it is connected to the primary system, it will regulate the pressure of the primary

cooling system. Temperature in the pressurizer is a bit higher than the saturation pressure

corresponding to the set pressure point. Electric heaters will be activated to vaporize a

fraction of the liquid to increase the pressure in the pressurizer and in the loop, whereas

sprays are designed to cool the vapor and condense it to decrease the pressure if needed.

The tank of the pressurizer should never become totally full of water (it would then be

impossible to decrease the pressure), or on the contrary full of steam (it would then be

impossible to increase the pressure). Figure 2-4 shows the pressurizer on the hot leg of

one of the loops. 172 is the outlet of the core. The Pipe 400 is the hot leg of the loop. The

branch will connect the pressurizer to the circuit. These elements will be described more

precisely later. The Pipe 443 represents the surge line. Pressure in this pipe is initially

taken as 2260.0 psi.

valve 444



\ 443

Figure 2-4: RELAP5 nodalization diagram of the pressurizer on the hot leg.

The pipe 441 represents the tank. The upper part of this tank, the pressurizer dome

(440), is represented by a branch-type element. Two valves can connect the tank to a

container represented by a single volume (449). The first valve (444), which is a

motorvalve, would represent the spray nozzle intended to condense vapor in case of a

need to reduce the pressure. This event would occur if the pressure becomes higher than

2350.0 psi and would stop when pressure goes lower than 2280.0 psi. The safety valve of

the pressurizer is numbered 445 in the model. The safety valve is also considered to be a

motorvalve. This valve opens when pressure in the tank becomes greater than 2575.0 psi

and closes when the pressure drops back below 2375.0 psi. Although the safety valve

diameter is larger (about three times) than the pressurizer porv, both of them open at the

same rate: 5.0 s-1 (the definition of the rate is also based on the normalized valve area).

Initial mass flow rates through the valves are of course equal to zero.

Primary cooling system tubing

The primary loop shown on Figure 2-5 contains the primary side of the steam

generator and a pressurizer (only one is necessary for the whole plant) on the hot leg and

the pump and an accumulator (used as emergency cooling system) on the cold leg. The

pressurizer has been previously described. The pump included in this system is the same

for each loop. It is possible to use built in models of pumps usually used (Westinghouse

for instance) or describe the characteristics of head and torque, pump head and torque

multiplier, and two phase difference curve. The characteristics were entered in this input

deck. The purpose of the accumulator is to act as an Emergency Core Cooling System

(ECCS) if needed. The accumulator has a built-in model in the code, as it requires special

numerical treatment. The geometry of the tank is needed, pressure (set at 615.0 psi),

temperature (set at 120.0 F) and tank initial fill conditions. Conditions are the same for

the three systems. As stated previously, the hot leg is mostly composed of pipes. The pipe

200 going from the vessel to the pressurizer is a 2-cells 9.16 ft long pipe. Initial pressure

is 2234.9 psi. The Branch 202 is used to connect the pressurizer to the hot leg. This

branch is included on each leg, but only one of them, the 4XX has a pressurizer plugged

on it. The Single Volume 204 is the inlet tube of the primary side of the steam generator.

Initial conditions in all these tubes are given with the pressure, the specific internal

energies for liquid and vapor phases and vapor void fraction. A single pipe represents the

tube sheet of the steam generator on the primary side. The Branch 206 is used to connect

the Tube Sheet 208 to the inlet of the Steam Generator 204. The tube sheet is represented

by the Pipe 208, which is divided into eight cells, so that it can represent the U-tube

geometry (as shown on Figure 2-5).

valve 4 44 4



For all these cells, initial conditions are set with pressure (from 2232.2 psi for the

inlet to 2170.4 psi for the outlet), taking into account the energy loss coefficients (0.048

at the junctions and 0.00015 for the wall roughness). Also, specific internal energies are

given and the vapor void fraction is set to zero. A Branch 210 is used to connect the

outlet of the steam generator primary side to the suction pipe of the reactor-cooling pump

(212). This simple pipe is composed of four cells, which allow giving a realistic

representation of the actual geometry of the system. The cold leg of the loop is

represented by a series of branches and pipes. The pump is plugged right after the

previously described Pipe 212 and sends water through a series of single volume

components 216, 218 and 220 (where is located the ECCS). The inlet of the reactor is

done with a single volume component that joins the Cell 200 in the reactor. A complete

summary of hydrodynamic components can be found in Appendix A.

Heat Structures

Heat structures in the core are of two types. The version of the code used is a

modified version of RELAP5. This includes a specific version to treat severe accident.

SCDAP/RELAP5 MOD3.2 allows fuel to melt and can model the relocalization of the

molten fuel and can predict the quantity of hydrogen produced in case of oxidation of

steam during LOCA. We obtained the code from ISS, Idaho Falls, Idaho. Different

materials are used to model all the heat structures in the code. We entered them in some

specific cards 201MMMNN for Heat Structures Material Properties.

Heat structures in the core

There are a lot of heat structures in the core. Especially, all the walls have

conduction properties that cannot be neglected. All vessel walls, thermal shields (in

addition to the vessel wall), different structures and plates present in the core to maintain

the assemblies and the integrity of the vessel, the core barrel and the core baffle are

modeled with heat structures. Geometrical models are to be chosen to model the heat

conduction (rectangular, cylindrical or spherical geometry). Other information essentially

needed are the heat structure mesh interval data, the composition of the structure, the

boundary conditions on the right and left side of the structure and the initial temperature

distribution. For all these heat structures no source of heat is used because the only power

source of the reactor is the fuel. The fuel can be modeled by RELAP5 type inputs, or

SCDAP inputs. SCDAP components were chosen. The fuel rods and control rods are

modeled using SCDAP components. For each channel (total of five), one component

"Fuel" for the fuel rod and one component "Control" for the control rods and

instrumentation tube are described. Information concerning radial mesh spacing in the

rod component (the code has built in models for fuel rods), initial radial temperature

distribution, power distribution, axially and radially, and information about the power

history of the fuel are specified. For each fuel rod component, the number of fuel rods

included in the model is specified. Table 2-1 shows the power fraction for each fuel rod

component and the number of rods associated with this component.

Table 2-1: Fuel rod power distribution.
Component Type # of rods Power fraction
1 Fuel 1020 0.04
2 Control 105 0
3 Fuel 4080 0.15
4 Control 420 0
5 Fuel 7344 0.24
6 Control 756 0
7 Fuel 12360 0.37
8 Control 1260 0
9 Fuel 7344 0.20
10 Control 756 0

The power fraction for each fuel element can be calculated knowing the number

of rods modeled in each element and the power peaking factors associated to each fuel

rod (or assembly). Radiation heat transfer option is activated in this study, so we have to

input the bundle matrix so that the code can calculate the view factor for every channel.

A typical bundle (here for the first channel) is described Figure 2-6.

1 1 111 1 1 1 1 1 1 11 1
111111111111 iii

1 1 1 1 1 1 1 1 1 2 1 1 1 1
1 1 1 I I 1 1 1 1 1 1 1 1
1 1 1 1 1 11 1 1 1 1 1 1 1 1
1 1 1 I l 1 1 1 1. 1 1
1 1 11 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 1 1 1 1 1 1 1 2 1 1

1 1 1 2 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Figure 2-6: Bundle matrix for the first channel: 1 is a fuel rod 2 is a control rod.

Heat structure in the loop

Heat structure concerning hot leg piping is simply a cylindrical geometry with left

boundary condition taken as the hydrodynamic volumes associated to the hot leg piping

and "0" as right boundary condition, meaning that we consider zero gradient of

temperature in the insulation of the pipe. Steam generator heat structures will be

discussed in the following part even if it is really part of the primary loop. Cold leg

piping and pump suction piping are also cylindrical geometry heat structures.

Heat structures in the steam generator

The heat structure of the steam generator is very complex due to the large number

of support plates for the tubes, and also, heat transfer between primary and secondary

sides need to be accurately modeled to represent the heat removal capability of the steam

generator. It is a heat structure that does the thermal connection between Pipe 208 (the

tube sheet of the primary side of the loop) and Pipe 276 (boiler part of the secondary side

of the steam generator). A cylindrical geometry is used with no heat generation inside the

heat structure. For this particular case, heat transfer is treated as convective phenomenon

at the boundary. The heat transfer coefficient can be calculated using a table that the user

can enter. In this study, the convective heat transfer coefficients were calculated with the

package included in the code. The secondary side of the steam generator also contains

heat structure components for walls (shell top, upper shell, middle shell, shell base, lower

shell) and for the wrappers (lower wrapper between 274 and 276) and upper wrapper and

swirl vanes (between 276 and 278, and 272 and 270). Calculations of heat transfers

through these components are also based on a cylindrical geometry. A summary of the

heat structures used in this model can be found in Appendix B.

Steady State and Initialization Process

The initialization run is a "steady-state" type run. Nearly implicit time

advancement has been chosen for the numerical time incrementation method. The power

for this experiment is set and controlled by a power table and equal to 2443.0 MWth.

Maximum time step allowed is 0.01s. The denomination time step does not have a real

sense here, because the steady-state option uses special numerical methods to reach

steady state more rapidly. So, saying that a time step is 0.01s means that the

incrementation step is 0.01. This time step was chosen low enough as a larger value

brings oscillations that prevent obtaining a good stability of the solution. Steady state was

reached after "time" advancement of 194.72, meaning 19472 iterations.

We can verify the initialization of our problem by looking the mass error m /M .

If this is low enough (104-10-5 or less), then we can trust the results of the steady state

run and run some transients. From the output file, we can read that for a total system

mass of 3.306E+5 kg, the mass error cumulatedd) is only 13.125 kg, which correspond to

a relative error of 3.9E-5, which is acceptable. And in the output file, "err.est" which

corresponds to the estimated truncation error fraction at the last advancement is only

equal to 2.3099E-08, which is very low. It is possible to adapt the problem to obtain a

good initialization level by playing with the friction factors for instance. Because it is a

parameter that is not exactly known, we can modify them a little bit from a previous

guessing hand calculation so that the problem is well initialized, meaning that the mass

error is as low as possible (or at least below an acceptable limit). Figure 2-7 shows the

evolution with time advancement of the solution for the steady state calculation. Once we

get this, we can start to set up transients.

605 -- 1st cell:
600 11101
6--0- -2nd cell
595 11102
-!- 3rd cell :
590 11103
-E --4th cell :
585 11104
5th cell:
580 11105
6--6th cell:
575 11106
S 7u 7th cell
S570 11107
565 X--8th cell
560 ---- 9th cell:
555 WM O10th cell
0 50 100 150
"time advancement"

Figure 2-7: Evolution of the temperature distribution during the steady state run.

The purpose of this picture is to show how, starting from the initial conditions

(temperatures all equal to 577 F=575.94 K), a final temperature distribution is reached.

The steady state temperature distribution in the core channel is presented in Figure 2-8.

The temperature distribution has an expected behavior. Whereas the central channels

have temperatures close to each other's, the outer channel has the lowest temperature and

the middle channels are in a kind of intermediate range. And the temperature profile

inside a fuel rod is given in Figure 2-9. These shapes may be obtained theoretically.

These types of calculation are the major interest of RELAP5. Depending on whether the

core model includes cross-flow junctions or not, the steady state conditions obtained are

not exactly the same, because coolant is allowed to flow in the transversal direction and is

not limited to only axial flow. Coolant flow remains, however, mainly axial. This is

treated in the Chapter 3.


600 -





X 575 -----first channel
570 ---second channel

565 -A-third channel
X-- fourth channel
560 fith channel

555 -
0 0.2 0.4 0.6 0.8 1
normalized axial node elevation

Figure 2-8: Axial coolant temperature distribution the steady state.

The purpose of this chapter was to describe the first tool used in the investigation

of the cross-flow calculation. However, results obtained with RELAP5 were not always

convincing, and RELAP5 is not designed to give the level of definition that was needed

to achieve a detailed calculation.

.1400 >
1200 -
1000 -
0 0.005 radius (m)01o 0.015

Figure 2-9: Radial fuel temperature profile at the sixth vertical node (height=2.012 m).

Results of the investigation will be discussed in more detail in Chapter 3. Another

method was then used to calculate cross-flow occurrence in the core in a more detailed

way, by using a Computational Fluid Dynamics (CFD) code. This method, more time

consuming may not allow modeling systems as large as RELAP5, but it allows having a

precise distribution of fluid velocity in an adapted core model.



RELAP5 is a proven tool for estimating the general properties of a large size

system, including all components contained in the reactor loop like pumps for example.

Cross-flow phenomenon in a nuclear reactor core was studied in different ways.

The effect of cross-flow can be of a major concern in a case of natural circulation of the

coolant in the core. RELAP5 calculates the radial mass flow rates by using a cross-flow

junction. The code solves the momentum equation between the two adjacent cells to

calculate cross-flow between these two cells. This model in RELAP5 needs to be

improved. In order to validate the models used in RELAP5, a similar model was created

to calculate cross-flow phenomenon using CFD methods. FLUENT (5) is a code that

performs CFD calculations. As the code solves the fundamental conservation equations,

results are supposed to be better. It solves the Navier-Stokes equations using a finite

element method on a grid that can be generated by the grid generator software GAMBIT.

FLUENT4 allows building simple grids with an integrated grid building option.

However, another version, FLUENT5, was used as options not available in FLUENT4

were needed, like User Defined Functions (UDF) that can be coded in C. UDF principles

will be discussed later in this chapter. GAMBIT was used as an external grid builder,

which is feasible as this problem is a simple case. The procedure of grid building with

GAMBIT is explained later on. Boundary conditions have a critical importance in the

modelisation of a system and it will be discussed that the major problem in achieving

cross-flow evaluation is to setup an adequate model. The boundary conditions available

in FLUENT, and those chosen in the model will be detailed in this chapter.

General Overview

Solving the Navier-Stokes equations is a real challenge, as they do not have an

exact analytical solution. Thus, numerical tools are necessary to solve numerically the

conservation equations. FLUENT is a Computational Fluid Dynamics (CFD) code that

solves Navier Stokes equations using finite differences methods. Numerical methods are

numerous to solve partial differential equations (PDE), and the one used in FLUENT is

described in this chapter. Although discretization of PDE's is powerful, it has important

limitations concerning stability of the solution and also from the size of the grid that is

chosen. The larger the grid is (or the finer the mesh size is), the more time is needed to

solve the problem. Also, the CPU time needed is generally in the order ofN2 for an NxN

grid. So, one can not perform calculations for the exact core model as modeling the

whole core with all the channels would require too many calculation cells and a quantity

of CPU and memory that a computer can not handle. To solve this problem, an equivalent

model is used based on a modelisation of the core with different annulus all entered on

the centerline of the core. This model will be described in this part too.

The basis idea of using FLUENT is to perform an evaluation of the cross-flow

models included in RELAP5. CFD methods, as they calculate the flow patterns and

properties at every point of a detailed grid, are more accurate than models used in

RELAP. Development of a model requires assumptions to match as well as possible the

actual conditions in the reactor. Also, FLUENT is solving the conservation equations for

two dimensions or three dimensions cases. This study will detail the advantage to do a

two-dimension run.

Basic Physical Models

For all types of problems, as shown in the user's guide (7), FLUENT solves

continuity and momentum conservation equations. If one chooses to include the "energy"

option in the setup of the problem, which is required to calculate temperature

distributions, heat transfer phenomenon or include compressibility properties, then an

additional energy equation for energy conservation is solved. The problem of interest has

a cylindrical geometry. Then only a 2D calculation is necessary. There are two types of

spaces for this model. The first one is a regular two-dimension space, which calculates on

a two-dimension axisymmetric domain. The last one has special properties that will be

described later. The conservation equations are written in a regular Cartesian coordinates

system for both spaces. For the regular two-dimension domain, the mass conservation

equation is written for each phase as follows:

+ a (p)= Sm
at ax

The term Sm is the mass source term, which is generally a phase change source term. In

the case of a single element, two phase problem, condensation or evaporation of the

element is included in this source term.

The momentum conservation equation is written:

a ( 8 + a aP vI
(p,) (pu ) + + pg, +-,
at al ax x

a ru + u 2 (u j
and z = L +K + - 3c,

cij is the stress tensor and Fi represents the external body forces.

The energy equation solved when needed is:

a(pE)+ (u,(pE+P))=- ka-k -hJ,,+ (,) +Sh
at ax Ox, a I

keff is the effective conductivity, Jj, is the diffusion flux of species j' and Sh represents any

volumetric heat source. This can be a User Defined Function.

It could be useful, in this case, to use a slightly different two-dimension space

model to perform our calculations. FLUENT uses a solver in a 2D axisymmetric space.

The conservation equations are a bit modified, in which axial and radial coordinates are

used (it is still a 2D problem) and axial (u) and radial (v) velocities are calculated. The

continuity equation for this coordinates system is written as:

p+ a(pu)+ (pv)+ PV= S
at Ox Or r

Momentum conservation equations for axial and radial momentum are more complicated

than for a Cartesian mesh:

8 1 8 1 8
-(pu ) + (rpuu ) + (rpvu )=
6t r ox r or
p 1 8 u 2 1 8 u 8v
---+--- r2 (V -) +--- rp -+- + F
Cx r c x Cx 3 r cr cr Cx X


8 18 1
(pv) +- (rpuv) + -- (rpvv)=

op 1 8 1v 2 Y (uv ou v 2 (V ) 2
--+-- r 2 (V ) +-- rL -+- r + -21-+ p-- +F,
dr r dr r 3 r ax x r r 3 r r
o is the swirl velocity.

Numerical Solvers

This part will present numerical solvers available in FLUENT. Different

algorithms can be used, depending on the type of problem that needs to be solved. Also,

boundary conditions are to be chosen among the list of available models presented in this

part. All types of boundary conditions are not always compatible with each other, or even

with the type of problem, depending on the geometry or the flow conditions for example.

Numerical schemes and discretization

There are two different numerical schemes in FLUENT, the so-called

"segregated" solver and "coupled" solver. The principle of both solvers is to solve the

integral form of the equations of conservation of mass, momentum and energy over a

mesh grid. Turbulence related equations (like in the k-s model) also employ these

solvers. A detailed description of both methods can be found in the FLUENT user

manual (8). This paper focuses mainly on the segregated solver.

Both solvers use a control volume technique. For this method, the domain is

divided in control volumes delimited by a computational grid. The grid choice depends

on the geometry, the type of problem to be solved (heat transfer, fluid flow), and the

resolution needed for the solution. One may choose simple orthogonal grid for simple

geometry where more resolution in certain regions of the domain is not needed. However,

if the focus is on boundary layers or some regions around obstacles, then one may need to

use a finer mesh grid in the interesting regions of the domain. The generation of this type

of grid may be very time consuming, and can need specific tools. The software GAMBIT

was used to generate a grid for two reasons. First, it is a convenient and user-friendly grid

generation software. Second, the version FLUENT5 does not generate grids by itself.

The governing equations, which are actually partial differential equations, are

integrated over a control volume to obtain a set of algebraic equations involving all

variables. These equations are then linearized and the system obtained is solved. These

two steps (linearization and resolution) are different depending on the solver.

In the segregated solution method, starting from an initial estimation or the

current solution, the fluid properties are updated. Next, the momentum equations for

every direction are solved using current values for pressure and face mass fluxes, in order

to update the velocity field. A Poisson type equation obtained from mass conservation

and linearized momentum equation is solved to correct pressure and velocity field so that

mass conservation is satisfied. When needed, scalar quantities like energy or turbulence

are calculated (by solving appropriate equations). The last step is to check the

convergence, which is tested on the residual of the unknown parameters. The

convergence level of each residual can be chosen separately. In the segregated solution

method, terms are treated explicitly and implicitly with respect to the iteration

advancement variable. The variable treated implicitly is the equation's dependent

variable. Thus, for each variable, we have a set of linear equations, one for each cell in

the domain. The principle of the segregated solver is that it solves a system of equations

for the variables one after the other.

Different discretization schemes are available for each governing equation. As

noted earlier, the code uses a control-volume technique. The governing equations are

integrated over the volume of the cell, such that the values of pressure and velocity are

needed at the cell faces. However, FLUENT uses a co-located scheme, meaning that both

pressure and velocity are stored at the center of the cell. Therefore, an interpolation

scheme is needed for pressure and velocity. Different schemes are available in FLUENT.

This study used the standard scheme for pressure described in reference (6). The

momentum interpolation utilizes a first order upwind scheme. It is assumed that the

variable face value is equal to the cell center value of the variable in the upstream cell. As

observed, the second step, after solving the momentum equations, is to perform a

pressure and velocity field correction. The pressure velocity coupling is realized with the

SIMPLE algorithm. This algorithm provides a pressure correction and is described in the

FLUENT user's guide (8). Since the energy equation must be solved, interpolation for

energy is done using a first order upwind scheme. This is also the case for k and F

transport equations. More extensive descriptions of numerical methods employed by

FLUENT can be found in the FLUENT user's guide (8).

Boundary conditions

In computational fluid dynamics, conservation equations govern the state of the

fluid and its behavior inside the computation domain. The state of the fluid at the

boundaries has to be well known to adapt the solution to the specific case to be solved.

The choice of boundary conditions and additional models will be decisive in the

matching of the FLUENT model to RELAP5 model. Basically, we need to have

boundary conditions for the outer walls, the inlet and the outlet regions of our domain.

FLUENT permits many different boundary condition (BC) types (9). A short description

of the most commonly used and those used for this model will follow. All the following

parameters can be fixed as constant or defined by a User Defined Function (UDF), or

simple piecewise linear distribution.

The inlet BC can be modeled using the following models:

The mass-flow inlet type BC. One can choose the inlet mass flux (kg/s) or mass

flow rate (kg/m2.s) for the chosen surface. One is also asked to define the inlet fluid

temperature and pressure. One should also provide the flow direction. It is possible by

writing the two components of the flow as work in two dimensions.

The pressure inlet type BC. One may choose the initial gauge pressure of the fluid

when it enters the calculation domain, also one need to provide the inlet temperature and

the turbulence model chosen.

The velocity inlet type BC. For a specified direction, one can choose the velocity

vector components, and the temperature of the inlet fluid. The turbulence model needs

again to be entered.

The outlet boundary conditions can be chosen among the following types:

The outflow type BC. One may simply choose the fraction of fluid that leaves

through this surface by choosing the flow rate weighting. This type of boundary condition

requires a fully developed flow.

The pressure outlet type BC. One needs to input the gauge pressure at the outer

surface and the temperature of the backflow with the turbulence model. The advantage on

the outflow boundary condition is that the flow does not need to be stabilized to have this

type of boundary condition.

These are basically the two interesting models of BC for the outlets. The walls

have specific types of BC that can be personalized to every kind of problem. As one

choose to solve the energy equation, one may impose boundary conditions that involve

temperature and heat flux characteristics.

The wall type BC. It can be largely personalized by choosing the type of thermal

boundary condition: it can be set by heat flux, temperature, convection, radiation or

mixed type of boundary conditions. For every kind of temperature, one will need to enter

the wall thickness and the heat generation rate, and also whether the temperature or the

heat flux at the surface or the heat transfer coefficient or the emissivity of the wall and

the temperature of emissivity. A mixed type BC will associate convection and radiative

heat transfer boundary conditions. Also in every case the user needs to input the wall

roughness and the material name. Parameters like heat flux, temperature, heat transfer

coefficient or emissivity can be modeled with a constant parameter or by using a UDF.

The symmetry type BC. If one chooses to have this type of boundary condition,

the normal velocity is equal to zero at the symmetry plane and normal gradients of all

variables are set to zero. This has to be used for a geometry that is actually symmetric, or

it can also be used to model "slip" walls, without shear stress.

The "axis" type BC. For an axisymmetric problem the physical value of a

particular variable at a point on the axis is determined by using the cell value in the

adjacent cell. For an axisymmetric geometry, the user should employ an "axis" boundary

condition type rather than the "symmetry" one.

The interior of the domain can be chosen as solid or liquid. In all these cases, the

liquid type of interior fluid is considered. Any kind of boundary condition can be

imposed at the limits, but this choice has to be coherent.

User defined functions

User defined functions (UDF) can be used to specify very detailed distributions of

boundary variables, source terms, property definitions, wall heat flux or solution

initialization. The UDF have to first be coded in the C language. The UDF is then

compiled by FLUENT (for the most basic utilization) and then information contained in

the UDF is incorporated in the FLUENT model to enhance the capabilities of the code. In

these tests, UDF were used to model a power source in the reactor, and inlet velocity

distribution, which is not elementary distributed. This type of function could also be used

to model the temperature dependence of the water density. As temperature varies slightly,

the density is expected to change only a few percent, which will have a negligible effect.

Power density is different in the channels of the core, due to the difference in peaking

factors. UDF can also be useful to set up a specific inlet velocity profile. In RELAP5, the

inlet velocity is the velocity averaged at the junction. However, for this model, a profile

following the 1/7th power law at the wall boundary can be proposed as an alternative for

inlet velocity profile. A UDF is used for this modelisation. The description of all possible

user defined functions can be found in the FLUENT user guide (9) and (10).

Turbulence Modeling in FLUENT

In this study, the flow is turbulent, as Reynolds number is very high (about 105-

106). Reynolds number is equal to Re = PVDh where Dh is the hydraulic diameter, p is

the fluid density, v is the fluid velocity and [i is the fluid viscosity. As there is a quite

high temperature fluid, viscosity is really low, which makes the Reynolds number high.

Also, there is a large diameter system and treating the problem like a laminar flow would

not make any sense. A turbulent viscosity term has to be taken into account, as seen in

the previous part. The way this parameter is calculated depends on the turbulence model

one chooses. Unfortunately, there is no single universal model for turbulence modeling.

Actually, turbulence modeling may need to be developed for each specific case, as it

depends a lot on the fluid type and properties, and on the geometry. Turbulent flows are

characterized by fluctuating velocity fields. As these fluctuations are of very small scale

or high frequency, it is difficult and too expensive to simulate them from an engineering

point of view. Governing equations can be averaged over time or ensemble so that

resulting modified equations are much easier to solve. But then, models are needed to

determine additional unknown variables. Turbulent viscosity can be calculated with

different methods that can require more or less parameters of the flow, like turbulent

kinetic energy k, or even the dissipation rate of the turbulent kinetic energy e .

The velocity field (and any scalar quantity field) can be decomposed in two

distinct components (mean and perturbation) such as:

u, = u +u'

The turbulent kinetic energy is defined as k = (v2 + v+ v2 ), using the perturbation

terms. As the average of the perturbation term is zero and as one takes a time-averaged

form of the continuity and momentum equation, one obtains, for the mean velocity,

written ui now:

ap+- a (pu)= 0
at ax,

Du, ap 8 u, u 2 8
P =- +- p- +-_8 +- pu'u'
Dt 8x, a8x ax8 ax, 3 Ix 8xK pl ,

These are the "Averaged Navier Stokes Equations". 6, is the Kronecker symbol. The

second term in the right hand side of the equation represents the divergence of the

viscous constraints tensor. The effects of turbulence appear in the Reynolds stress term

T', = -p u', u', that requires a specific modeling to close the Reynolds averaged

momentum equation. In analogy to the viscous stress term in laminar flows, the turbulent

stress term has been proposed to be proportional to the mean-velocity gradients. This

concept has been proposed by Boussinesq in 1877, and has been used as a basic to

numerous turbulence models. The Boussinesq eddy viscosity hypothesis for the

divergence of the Reynolds constraints tensor is written as:

ou, Ou 2
-pu' u', = Pt, + pk8

The term [tt is called turbulent viscosity and has to be determined to be able to

evaluate the Reynolds stress. Contrary to the molecular viscosity [t, [tt is not a fluid

property but depends strongly on the state of turbulence. So [tt may vary greatly from one

point to another in the fluid, especially close to the boundaries, as in walls for instance.

Many models have been developed to calculate turbulent viscosity. Some of them are

very simple and based on empirical correlations, but other models are more complex and

mix theoretical and empirical considerations to setup transport equations for turbulence-

related quantities used to obtain the turbulent viscosity. So the Boussinesq eddy viscosity

hypothesis does not constitute the whole turbulence model, but it provides the framework

for constructing a real model. The main problem is actually to determine the distribution

of [tt. A comprehensive review of turbulence models and their application in hydraulics is

presented in reference (12) and (13).

An algebraic model may be sufficient for this simulation, as there is a very simple

geometry. The principle of an algebraic model is to determine the turbulent viscosity only

by one single relation. The most common zero-equation model is also called Prandtl

mixing-length model. A characteristic length of the turbulence is introduced: the Prandtl

mixing length. This is basically the length along which an eddy remains identical

mechanically. The turbulent viscosity is then calculated as:

,2 SU
1)t = 12

U is the mean velocity gradient along the transversal direction. 1m is the mixing

length. The calculation of lm is then the problem. A satisfactory estimation for lm is

illustrated in Figure 3-1. The parameter d is defined as the layer width. It is calculated as

the distance from the wall to the 1% point of the outer edge, or the point where velocity

differs by 1% from the free stream velocity. The constant K is the Von Karman constant

and is equal to =0.41. A is equal to A =0.09. These constants have been determined

originally by Patankar and Spalding (14).

For developed duct flows (channels, pipes), the mixing-length distribution is well

described by Nikuradse's formula (15):

lm/R = 0.14-0.08(1-y/R)2-0.06(1-y/R)4


Xx Qs5 yI6 1
Figure 3-1: Mixing-length distribution in wall boundary layers.

However, this is a very approximate modelisation, and as a high intensity cross-

flow phenomenon is not expected, one may want to have a better model to treat

turbulence so that the order of magnitude of the error introduced by the turbulence

treatment would not affect the results of cross-flow calculation. Furthermore, such a zero

equation model is not available in FLUENT.

Six models of turbulence can be used in FLUENT, each of them being adapted to a

certain type of situation. They are: Spalart-Allmaras model, k-e model (Standard (16),

RNG (17) and Realizable (18)), Reynolds stress model and Large Eddy simulation. Two-

equation models are the most widely used in engineering simulations. This study will use

k-e models, which are respectively one-equation and two-equation models and use a

Boussinesq eddy viscosity hypothesis for the divergence of the Reynolds constraints

tensor. Spalart-Allmaras model could also be used without any problem

One-equation Spalart-Allmaras model

Although this model has been originally developed for aerospace flows involving

wall bounded flows, it has been adapted to fit with the needs of FLUENT and high

Reynolds number calculations. However, as this model is quite new, results in some very

specific cases should be evaluated with caution. This is a simple one-equation model that

uses the Boussinesq hypothesis. The main advantage of this method is that one does not

need to calculate a length scale related to the local shear layer thickness. The unknown of

the transport equation solved in the Spalart-Allmaras model is a modified form of the

turbulent kinematic viscosity. Constants used in the model are these proposed by default

in the code. A detailed description of the Spalart-Allmaras model is done in the FLUENT

user manual (7) and in the original paper (19).

Two-equation k-epsilon model.

This k-e model is used as a better alternative to the one-equation Spalart-

Allmaras model. This is a very common turbulence model as it is based on semi-

empirical considerations. Since it was proposed by Launder (16), the k-e model has

always been used and improved. In this model, the transport equations of turbulent

kinetic energy (k) and its dissipation rate (e) are solved. Basically, the transport

equations of k and s are built on the classical balance equation:

Rate of change + Convection = Diffusion + (generation destruction)

The transport equations for the Standard k-s model are:

Dk a 9, Sk
P =a[-It Ok + G, + Gb P-YM
Dt ax, k 2J+x,

where Gk is the production of turbulent kinetic energy defined as Gk = -pu', u',

determined using Boussinesq assumption, Gb is the generation term due to buoyancy,

put OT
defined as Gb = fig,, L YM is the dilatation dissipation term for high Mach number
Pr, 9x,

flows, not relevant here. ak is the turbulent Prandtl number for k.

pC + +C,C (Gk 3+G,)-C p
Dt aO, a~x, k k

where C1,, C2, and C3, are constants, O7, is the turbulent Prandtl number for e In

FLUENT, constants have the default values C1,=1.44, C2 =1.92 C3 =0.09, ak =1.0 and

O, =1.3. These values have been determined experimentally. However, the user can

change them in FLUENT. This study used only the Standard k-e model without

changing the default constants. A detailed description of other two-equation models can

be found in the user manual (7).

In this model, the turbulent viscosity is computed as [tt=pCk2/s. This value can

then be used to calculate the Reynolds stress, from the Boussinesq assumption. The main

difficulty is to choose the boundary conditions for k and e. The user guide manual

proposed the following approximations to evaluate the initial conditions for k and e .

The turbulent kinetic energy k (m2/s2) can be estimated from the

u' 3 (Uagi)2
turbulence intensity I defined as I =- 0.16(ReD : k = (u
u 2

The turbulent dissipation rate e (m2/s3) can be estimated from a length

3/4 k3/2
scale = 0.07DH. Then e is approximately equal to e = C4 k The

constant Ci has been determined equal to 0.09.

This research used these methods to determine our inlet boundary values for k and e .

Defining a CFD Method for Evaluating RELAP5 Cross-Flow Model

RELAP5 uses a coarse integral method to solve the conservation equations. The

cross-flow model is very simplified, as presented in the Chapter 2. A protocol has been

developed to determine an evaluation of this model. A very simple system has been

studied with cross-flow calculations. This part presents the two different models used to

perform the simulation with FLUENT and RELAP5. The system consists in a simple

vertical pipe where the fluid (liquid water) is moving upwards. Conditions of operation

are these of a real reactor. Inlet fluid temperature is 555 K, and the pressure is 15.68

MPa. The inlet conditions are very specific. This study used a velocity profile distribution

to analyze the impact of velocity distribution for cross-flow calculation, as velocity is a


parameter that cause cross-flow. Figure 3-2 represents the system that has been calculated

with FLUENT and RELAP5. Also, it is possible to have internal heat generation in the

models to study the effect of temperature on the cross-flow occurrence. The velocity

profile is radius dependent, with the central part of the flow at a higher velocity.

inlet velocity profile

Figure 3-2: Basic geometry studied with RELAP5 and FLUENT.


The FLUENT model is supposed to be as close as possible to the real geometry

and operating conditions. Figure 3-3 shows the calculation domain representing the pipe

of the Figure 3-2.

p4 wall P3
inlet outlet

P L axs 92

Figure 3-3: FLUENT calculation domain of the system.

Calculations were performed on a two dimensional axisymmetric space, as

boundary conditions were applied with a cylindrical geometry. The distance ri does not

represent an actual boundary, but it is used as a limit to impose internal heat generation,

and also it will serve as a basis to define velocity profiles. Also, this is the limit that has

been used to build the RELAP5 equivalent model, where the channel was divided in two

pipes linked together with the cross-flow junction of interest. The calculations were

performed on this simple surface, so a grid was generated for this domain. The software

GAMBIT was used.

GAMBIT grid generation

Grid generation is a very important step in CFD problems. Generally, grid

generation represents a large fraction of the total CPU time necessary to the resolution of

a problem. If the geometry of the domain is complex, it can take a long time to build an

adequate grid. This study used a very simple geometry. Since the calculation domain is a

rectangular area, an orthogonal grid was used. The GAMBIT software generated the grid,

as this task cannot be handled by FLUENT5 (contrary to FLUENT4 that allowed

building simple grids). Following is a description of the few steps necessary for grid

construction and the establishment of some reference points.

The first step is to set the limit points that will be used to create the edges.

GAMBIT allows creating the points with their coordinates.

The next step is to create the edges in the domain. In the present case, the edges

are the external limits of the domain. So edges are created for both walls and open

boundaries. On Figure 3-3, the four edges are created with the four points P1, P2, P3, and

P4. They correspond to the inlet zones, one outlet area, and the two side boundaries,

corresponding to the surrounding wall and the centerline of the physical domain that is

modeled by a real boundary. Once edges are created, it is possible to create faces and to

build volumes. This option was not needed in the 2-D problem.

The next step can be very time consuming in term of CPU time. However, for this

simple domain, the grid generation is very fast. For an orthogonal grid, the first step is to

generate the mesh edges. For every edge, a one-dimensional mesh is generated, each cell

having the same length. It could be possible to create a finer grid at the boundary layer if

the purpose of the study was to have detailed description of the events occurring close to

the wall surfaces. Finally, the last step is to create the final mesh of the domain. In

GAMBIT, the "mesh faces" option generates the final mesh, based on the mesh edges.

Then, after assigning specific boundary types to the boundaries, it is possible to export

the grid to FLUENT, where materials properties and boundary conditions are specified.

Boundary and operating conditions

The boundary conditions are very important, and will be the determinant in the

validity of the model. As a commercial code is used that is reliable and uses powerful

numerical methods, the only way one can adapt the solver to a problem is by setting good

and adapted boundary conditions. The inlet and outlet boundary conditions will

especially have a capital effect on the behavior of the fluid in the calculation domain.

Operating conditions. The operating conditions need also to be specified before

beginning the calculations. The reference for the pressure is at the outlet of the

rectangular calculation model, on the centerline. Also, an option can be activated to take

into account the effect of gravity by setting the acceleration of gravity (9.8 m/s2) and also

the direction to which it applies (vertical top to bottom). Operating pressure was chosen

equal to 15.68 MPa.

Outlet boundary. The outlet was set as an outlet pressure boundary with a static

gauge pressure equal to zero (which means that the absolute outlet pressure is the

operating pressure). An outflow could not be used, as it requires the flow to be fully


Wall boundaries. The centerline boundary was represented by an axis type

boundary condition. The outer wall ("wall" on Figure 3-3) was of a wall type. An

adiabatic wall was chosen for the outer boundary (meaning heat flux set to zero). Choice

of wall roughness will determine the losses due to friction on the wall and of course have

a strong effect on the velocity profile (at least on the transversal direction) for a domain

with a mall radius, where the wall boundary layer can be visible. A roughness height (ks)

equal to 10-4, and a roughness constant approximately equal to 0.5 are acceptable


Inlet boundary. The inlet boundary was set to have the velocity fixed by an inlet

velocity type BC. As velocity inlet boundary conditions was set, other parameters like

turbulence terms boundary values need to be specified.

Fluid properties. As one considers an incompressible flow, the material properties

are almost constant values. However, water properties tables (20) were used to set the

water properties respect to temperature as a piecewise linear function. This has an effect

only when heat is added to the fluid:

Two different domains were studied. The first one represents a large pipe, where

the wall effects were avoided because of the large diameter. The second one is a much

smaller one. Based on the geometry shown Figure 3-3, two models were built. The first

one had an outer radius of 0.4 m, and the pipe was divided into two sub channels such

that ri=0.2m and ro=0.4m. The length of the pipe L was equal to 5.0m. The second model

was a much smaller one. The dimensions of the pipes were ri=0.5cm, ro=1.0cm and the

length of the pipe was taken 20 cm. Results of calculations on these two models are

presented in the Chapter 4.

RELAP5 Equivalent Model

An equivalent model was built for the analysis of the cross-flow phenomenon

inside the pipe. For this, the pipe was divided axially into two pipes. The Figure 3-4

shows the nodalization diagram of the equivalent system. Pipes, which area is the same as

the area of the two sub channels that compose the pipe Figure 3-3, represent the two axial

channels. This means that the cross section area of the pipe 100 is equal to 7r ri2. The area

of the pipe 110 is then equal to 7r (ro2-2ri). Volumes 88, 98, 130 and 131 are time

dependent volumes and represent reservoirs to provide the water flow in the system.

Pressure in 130 and 131 set the pressure in the system. Pressure has then to be adjusted in

volumes 88 and 98, where the temperature of the entering fluid is set. Time dependent

junctions 89 and 99 regulate the inlet velocity for the system. The junctions 120 and 121

are simple single junctions. The set of Junctions 112XX is a multiple cross-flow junction

that allows the fluid going from a channel to the other. This is the junction of interest in

this study. The pressure in the two top Volumes 130 and 131 is set to 15.68 MPa and is

not changed. Only the pressure in the inlet time dependent volumes is adapted to each

problem. The temperature in the inlet Time Dependent Volumes 88 and 98 is set to 555



The inlet velocity regulated by the time dependent junctions 89 and 99 is the

parameter that sets the flow in the system, and is specified for each test in the Chapter 4.

The cross-flow junction area in the RELAP5 model has been set as the outer area of the

cylinder of radius ri and height equal to the height of a nodalization cell. Heat structures

are used to provide heat to the system and are attached to the volumes 110 and 100. An

additional heat structure is added to realize heat conduction between channels 1 and 2. A

copy of the RELAP5 input deck is provided in Appendix C.

11019 10019

11018 10018

11202 10002
111 0 I 1

11001 10001

i99 819

98 10018

Figure 3-4: Nodalization diagram of the equivalent RELAP5 model.,

The results of the equivalent simulations with FLUENT and RELAP5 are

presented in Chapter 4.


This chapter will present the results from different simulations with FLUENT and

RELAP5. For each case, cross-flow has been calculated. First, RELAP5 has been run for

a steady state case and a power surge transient to characterize cross-flow phenomenon in

a real case simulation. And then, the attention is focused on the cross-flow phenomenon

itself by using the two models described in Chapter 3.

Surry Reactor Calculation for a Steady State and a Simulated Power Surge Transient

Cross-flow is predicted in both the steady state and transient cases. The results

will show that the impact of taking into account cross-flow in the core is real.

Steady State

In the case of a steady state calculation, no attention is focused on the evolution of

the different variables during the calculation phase. Those values have indeed no physical

meaning, since the code uses numerical tricks to reach steady state more rapidly, and then

the iteration step size does not correspond to a real time value. By looking at the mass

flow in the axial channels, different results can be observed whether cross-flow is

considered or not. Without cross-flow junction, the mass flow is constant along the core

length (due to mass conservation). Though, velocity is increasing as temperature

increases and density decreases. The relative mass error introduced by numerical

calculations along the channel is less than 10-4 for a core without cross-flow. Table 4-1

shows the axial mass flow in the five hot core channels for the case of a purely axial core

or a core with cross-flow junctions between the axial channels.

Table 4-1: Comparison of axial junction mass flow for a core with
and without cross-flow.
Axial mass flow (kg/s) for core
with cross-flow
junction# \ channel 1 2 3 4 5
inlet 387.24 1549.0 2788.2 4647.0 2788.2
1 387.09 1548.4 2788.2 4648.0 2787.8
2 386.76 1547.1 2787.6 4648.9 2789.2
3 386.51 1546.1 2787.2 4649.9 2789.9
4 386.3 1545.2 2786.8 4650.3 2790.9
5 386.02 1544.2 2786.4 4650.9 2792.1
6 385.79 1543.2 2785.7 4651.1 2793.8
7 385.74 1543.1 2786.2 4651.8 2792.9
8 385.22 1540.9 2784.8 4653.3 2795.5
9 385.51 1542.3 2785.4 4650.4 2796.2
outlet 383.8 1535.4 2784.4 4662.1 2794.2

Axial mass flow (kg/s) for core
without cross-flow
actual data 382.26 1529.1 2762.8 4615.1 2768.2
normalized data 385.49 1542.002786.11 4654.04 2791.56

For a more accurate comparison, data have been normalized to a same total core

mass flow. The mass flow rate for a closed-channels core is close to the one at the

different axial junctions along the core allowing cross-flow, but still a difference is

noticeable. For a given channel, the closed-channel value is roughly an average value of

the axial mass flow in the open channel case. The difference comes from the cross-flow

consideration. However, as can be seen on Figure 4-1, in the steady state conditions,

cross-flow is much smaller than the main axial mass flow, which is expected. This figure

shows the magnitude of cross-flow between the five different channels that compose the

core. The sign of the calculated cross-flow represents its direction. The way the channels

are numbered in the RELAP5 inputs, a positive value for cross-flow means that the flow

is going outwards. The total cross mass flow is only about 0.5% of the total axial mass

flow. This is due to the fact that the pumps drive the flow axially.


-- from channel 1 to 2
5 from channel 2 to 3
From channel 3 to 4
4- -X- -from channel 4 to 5


-4 -2 0 2 4 6
mass cross flow (kg/s)
Figure 4-1: Cross-flow magnitude versus height in the core
for a total power of 2443MW.

However, to have a better idea of the relative importance of cross-flow from a

junction to the other, it is better to look at the mass flux at the junctions. Figure 4.2 shows

mass flux for cross-flow versus cell height in the core. This one is calculated by dividing

the mass flow by the cross-flow junction area. Actually, the direction of the cross-flow is

more interesting. According to the sign of the mass flow, and looking at the setup of the

junctions, there is a cross-flow going outwards for all central channels. The outer junction


is a little bit different. This kind of prediction is surprising as different parameters can be


Temperature: the temperature is higher in the central part of the core, which

should, by density difference, push the flow inwards.

Axial velocity: the axial velocity is changing along the core, but always the

axial velocity in the central parts of the core is higher than at the periphery,

which means that the pressure difference induced by velocity differences

should also drive the flow inwards.

If the conditions are not constant anymore, and if a transient phase is taken into

account where conditions are changing in the core, then cross-flow will be different.

- -2443MW channel 1 to 2
- -2443MW channel 3 to 4
- -3500 MW channel 1 to 2
- - 3500MW channel 3 to 4

-4 -2 0 2 4 6 8 10 12 14
mass flux (kg/m2.s)

Figure 4-2: Cross-flow mass flux versus core cell height.

16 18

Now different power levels can be considered. RELAP5 will predict different

cross-flow magnitudes for different power levels. As the mass flow is set constant in the

system, only the temperature of the fluid has an effect on the phenomenon. For different

power levels, the steady state cross-flow value was calculated, and normalized to the inlet

mass flow in the reactor. Figure 4-3 shows that in a general manner, cross-flow

magnitude is increasing when power is higher.

Sr ---1st cell ch. 1 to 2
cross flow ratio
-- 2nd cell ch. 1 to 2 cross flow ratio
channel 1 to 2 3rd cell ch. 1 to 2
---4th cell ch. 1 to 2 channel 3 to 4 ,
-o *-5th cell ch. 1 to2 o
-A- 8th cell ch. 1 to 2 o
-A-1st cell ch. 3 to 4
-X 2nd cell ch. 3 to 4
-3rd cell ch. 3 to4 4
o -G--4th cell ch. 3 to 4 C
u --5fth cell ch. 3 to 4
d ---8th cell ch. 3 to 4


+ o
o o

o 2443 2000 Power (MW) 1500 1000 o

Figure 4-3: Ratio of cross-flow to axial flow for different power levels and for different
height in the core.
W I-

6 2443 2000 Power (MW) 1500 1000 (?

Figure 4-3: Ratio of cross-flow to axial flow for different power levels and for different
height in the core.

This figure represents the ratio of the fluid flow in the radial direction to the

bottom flow in the axial direction. These four cases represent steady state runs. The

cross-flow is noticeable very low. However, a good point to remember is that the mass

flow in the loop is maintained constant, and then for each case, there is approximately the

same mass flow in the core. Therefore, cross-flow cannot be mainly caused by mass flow.

The parameters that can make the flow go radially would be the density difference, since

the temperature profile in the core is not flat, lowering the density in the central part

where flow temperature is higher. Indeed, looking at the level of the fifth cell

(Hydrodynamic Volume 11105, etc...) reveals a difference of temperature between the

center channel and the outer channel of about 1%. This makes a very small difference.

Density variation is also small. But added to this phenomenon, axial velocity is higher at

the center of the core than at the periphery.

Transient Case: Simulation of a Power Surge Transient

The transient case being considered consists in a power surge transient. RELAP5

allows using two methods to input power in the core. One is the reactor kinetics

calculation. But as the interest resides in the thermalhydraulic behavior of the core, the

other method using a power table and setting the total power in the core was preferred.

The use of the power fraction factors in the heat structures does the repartition of the

power in the different channels. The power transient is represented in Figure 4-4. The

reactor model is built to keep a constant mass flow rate in the reactor equal to 4100.0 kg/s

in each leg. Then, as the inlet mass flow in the core remains constant, the impact of the

power level on a steady state can be observed, and also during the transient phase, as

temperature increases and physical properties are consequently modified (much more

than during the temperature change across the core).

Total Power [MW]
3500 -

2443 -

240 250 260 time (s]
Figure 4-4: Power transient simulated in the transient test.

The calculations were done using the restart option of RELAP5, starting from the

steady state determined in the previous part. The general conclusion that can be made

after analysis of the results of this investigation is that the power level in the reactor has

an effect on cross-flow magnitude. As mass flow in the core is constant, fluid temperature

increases when increasing power in the reactor.

Cross-flow magnitude is larger for a greater power level. As can be seen in Figure

4-5, and also as in the previous part in Figure 4-3, the mass flux through the cross-flow

junctions is larger when power is larger (whatever the sign of the cross-flow is).


-*-channel 1 to 2
0 -*--channel 2 to 3
--- channel 3 to 4
2 ---channel 4 to 5

0 3 6 9 1 15 18
mass cross flow (k gs)

Figure 4-5: Cross-flow magnitude versus height in the core for a total power of 3500MW.

To illustrate the transient phase of the cross-flow, focus will be drawn on the

central channel. The cross-flow ratio is defined as the ratio of cross-flow to inlet flow in

the core is calculated and its variation due to a power surge can be analyzed. The flow

rate range will not be so much of a problem as cross-flow is normalized to the inlet flow

rate, but instead, parameters like temperature will affect the cross-flow variation. Figure

4-6 and Figure 4-7 show the variation of this ratio for a power surge of 1000 to 1450

MWth. The first comment is that it takes some time to reach a new steady state once the

power surge is simulated. Since it is an instantaneous large step of power, it is expected

that the flow will oscillate a little bit before reaching its new level. Results show in both

cases that the cross-flow magnitude tends to increase after the power surge. The higher in

the core, the more the oscillations last before reaching a steady state for the new


8th cell


O2H+ --- c -- w ----first cell

0 9 5 5 ?rfl ?C 97f4th cell97.


transient time (s)

Figure 4-6: Cross-flow to inlet flow ratio for a 1000MWth case. Cross-flow is from the
Channel 1 to 2.

The figure 4-7 shows that for the seventh and the eighth cell, cross-flow takes

longer to stabilize and has also more oscillations in the transient phase. However, for the

bottom-core cross-flow junctions, the steady state is reached quite smoothly, but it still

requires some time before the level of cross-flow stabilizes. As the temperature in the

core is increasing with certain inertia, it is understandable that cross-flow, which also

depends on the fluid properties, will reach its steady state conditions with a delay.

Furthermore, it is logical that the flow through the top core cross-flow junctions would be

disturbed, since it picks up perturbations from the lower part of the core.


-8 7th

Figure 4-7 Cross-flow to inlet flow ratio for a 2443MWth CS. COSS-flOW is om

o 7th cell
o 245 250 255 260 265 270 275 280
transient time (s)
Figure 4-7: Cross-flow to inlet flow ratio for a 2443MWth case. Cross-flow is from
Channel 1 to 2.

Cross-flow phenomenon is now characterized with RELAP5. An expected result

is that cross-flow increases as power increases. However, the cross-flow direction is not

going in the expected way. As seen previously, the temperature distribution in the core

and the fluid velocity distribution should both push the flow inwards, due to density

difference and pressure difference, even if these two phenomena have only a little effect

due to the small differences in the distributions.

Cross-flow Between Two Adjacent Channels: FLUENT and RELAP5 Results

This part presents the results of cross-flow calculation in the simple system

presented Chapter 3. This model has been developed to calculate cross-flow between two

axial streams of fluid. This should allow doing an evaluation of the cross-flow model

used in RELAP5. FLUENT and RELAP5 were used to calculate the cross-flow through

an imaginary boundary inside a vertical pipe. The main interest is that, as we consider a

system where the inner channel is an open channel, special model is needed in RELAP5,

which is doing calculations using closed channels. This is the function of the so-called

cross-flow junction. The geometry of the system is reminded in Figure 4-8.

Figure 4-8: Two-sub channel pipe studied with FLUENT and RELAP5.

The first set of calculations is applied to a quite large system. The system is 5.0 m

long, ri=0.2m and ro=0.4m.

Simulation in a Large-Diameter Vertical Pipe

The flow is expected to behave in a different manner whether its temperature

remains constant or changes due to a power density. Indeed, the temperature of the fluid

affects its energy, and the fluid properties may change consequently as its temperature

increases. To analyze this effect, tests have been done with and without power. The

results from FLUENT calculations presented in this paper show the solution obtained by

the code on the calculation domain. So, the domain presented in all the following

FLUENT figures is actually an elementary angular domain of the pipe like represented on

Figure 3-3.

If no power is added to the fluid

The first test performed with this geometry was with a different inlet velocity in

each of the imaginary channel. The velocity profile is described Figure 4-9, profile (a).

velocity (m/s) velocity (m/si

6.0 -

2.0 2.0- I

0 0.4 radius m) 0 0.2 0.4 radius (m)
Figure 4-9: Inlet velocity profiles used in cross-flow simulation.

Two tests with different pipe length were performed. The first one was over a

50m long pipe. The grid on this area is a non-uniform mesh of 500 x 60 mesh points. The

mesh is much finer in the regions of high gradients or of regions of interest. These are the

central region at r=0.2m, and the entrance region. For these regions, GAMBIT allows

having a finer mesh by setting a ratio between the dimensions of adjacent cells.. As seen

on Figure 4-10, this pipe is long enough such that the flow becomes fully developed. The

evolution of axial velocity shape is displayed on Figure 4-10. The mixing length can be

determined as the point where the flow gets his final shape that follows a one-seventh-

power law. The radial velocity is the cross-flow velocity, and data at r=0.2 have been

collected on a second run with the same boundary conditions, but on a shorter domain

such that only the beginning of the mixing process appears. This second run was done

only on a 0.4mx5.0m domain. The mesh grid for this simulation was a 500 x 40

orthogonal mesh. Figure 4-11 shows a more detailed view of the mixing process. Both

plots represent a length of 5 meters. Figure 4-11 shows the calculation domain. This is

actually an elementary angular domain of cylinder. The centerline is located at the bottom

of the calculation domain, as explained previously.









0.00 0.10 0.20 0.30 0.40
radius (m)
Figure 4-10: Axial evolution of the axial velocity profile. z is the height in meters.

z=8.3 z=6.25 z=4.2



I7 t '5Z- .

are 4-1 1: Kiaalai (v-velocity) ana axial tu-veloclty) velocity alsmounon
at the entrance of the tube for Vi=6.0 m/s and Vo=2.0 m/s.

Figure 4-12: Radial (V-velocity) and axial (U-velocity) velocity distribution
at the entrance of the tube for Vi=4.0 m/s and Vo=2.0 m/s.

The sign of the radial velocity shows two behaviors. First, the fluid at lower

velocity, right at the entrance is attracted by the region of low velocity by pressure

difference. Then, the fluid is going from the central part of the flow to the outer part.

There is a transfer of momentum of the liquid from the high velocity part to the lower

velocity region, and then, the fluid is moving outwards to reach the one seventh power

law distribution as seen above.

The same simulation was realized with the inlet velocity profile (b), as shown on

Figure 4-9, and results are shown Figure 4-12. With such a lower velocity, cross-flow

intensity is strongly reduced. The attraction of fluid at the entrance (the blue region where

radial velocity is negative) is less important in case of an inner velocity of 4 m/s (Figure

4-12) than for an inner velocity of 6 m/s (Figure 4-11).

-- Vo=2.0m/s and Vi=4.0m/s RELAP5 --Vo=2.0m/s and Vi=6.0m/s RELAP5
Vo=2.0m/s and Vi=4.0m/s FLUENT - - -Vo=2.0m/s and Vi=6.0m/s FLUENT



4 [1

1 2.5 11

1.5 ]

0 El I I I ,

0.00 0.01 0.02 0.03 0.04 0.05 0.06
velocity (m/s)

Figure 4-13: Comparison of cross-flow prediction between FLUENT and RELAP
for two simulations.

Data got from these runs have been plotted and compared to the same conditions

solved with RELAP5. Figure 4-13 shows the predictions of cross-flow with two different

methods: RELAP5 and CFD with FLUENT.

There are differences between the predictions of cross-flow. Results are consistent

as in average, cross-flow for the case where Vi is equal to 6.0 m/s is about 2.5 times the

magnitude of cross-flow for an inner velocity of 4.0 m/s, but still, FLUENT prediction is

about ten times larger than RELAP5 results. Also, RELAP5 does not predict the

oscillations at the beginning of the pipe. This is explained by the fact that RELAP5 uses

integrated quantities over a much larger volume than FLUENT, and these transient

phenomena are then eluded. Other results concern the axial velocity and are presented on

Figure 4-14.

-- case Vi=6.0 m/s and Vo=2.0 m/s channel 1
X case Vi=6.0 m/s and Vo=2.0 m/s channel 2
-- case Vi=4.0 m/s and Vo=2.0 m/s channel 1
-- case Vi=4.0 m/s and Vo=2.0 m/s channel 2


0 --- i------ i-----'' --

1 2 3 4 5 6 7
velocity (m/s)

Figure 4-14: Evolution of the axial velocity for the two RELAP5 axial channels along the
pipe, for the two inlet boundary conditions (a) and (b).

FLUENT predicts a mixing much faster than RELAP5. A weakness of RELAP5

is that it considers pipes as closed channels. Drag force from a channel stream to the

other needs to be taken into account. Even with a low viscosity fluid, FLUENT predicts a

mixing length much smaller than that. Cross-flow phenomenon is not treated with enough

accuracy to calculate this problem.

If power is added to the fluid

The last situation studied with this specific geometry is the case of a channel

where the fluid is entering at a constant velocity over the cross section area, but the

central part of the channel is heated. The case was studied with a non-constant power

density over the cross section of the pipe. The simulation presented here considers that

the central channel is heated with a power density equal to 66.1 MW/m3. The velocity of

the fluid (water) entering the volume is 2.0 m/s everywhere.

Figure 4-15: Temperature distribution over the total length of a 20.0 m pipe heated in the
central channel at 66.1 MW/m3.

Figure 4-15 shows the evolution of the temperature of the fluid as it goes all the

way through a 20.0 meters long pipe, heated in its central part.

But the major interest is the cross-flow occurring due to the heating of the fluid.

Figure 4-16 shows the results of the calculation of radial velocity in a shorter domain.

The 5.0 meters pipe was discretized in a very fine mesh grid (40 x 500).

Figure 4-16: Radial velocity for a constant inlet velocity and a heated central channel.

On Figure 4-16, the first three pictures (and the scale on the left) represent the

general behavior of the radial velocity along the 20.0 meters long pipe heated in the

central channel. The fourth one focuses on the inlet region. Figure 4-17 shows the radial

velocity at a radius of 20 cm (which corresponds to the interface between the two

channels). These results have been obtained with FLUENT and RELAP5. The negative

sign means that the fluid is moving from the outer channel to the inner channel. This is

illustrated with the blue and green regions on Figure 4-16. This graph shows that both

FLUENT and RELAP5 predict the same behavior of the cross-flow velocity at the


entrance of the pipe. However, there is a net difference of behavior of the predictions

after that. The following Figure 4-18 shows that axial velocity is changing significantly

as temperature increases. The density of the fluid heated in the central channel is

becoming lower and lower and the velocity of the water increases.


FLUENT results
-a- RELAP5 results

-0.025 -0.02 -0.015 -0.01 -0.005 0 0.005
fluid velocity (m/s)

Figure 4-17: Cross-flow velocity at the artificial interface between the two channels.

5 --*-axial velocity in the Channel 1 (heated)

4.5 -NX-axial velocity in the Chanel 2 (not heated)






E 1.5


1.6 1.8 2 2.2 2.4 2.6 2.8 3
velocity (m/s)

Figure 4-18: Axial velocity in Channel 1 (heated) and Channel 2 (non-heated)

Whereas RELAP5 does not predict a big change in the temperature of the fluid in

the outer channel, FLUENT predicts a much faster heating of the second channel as seen

on Figure 4-15. Water from the outer channel is attracted to the central stream, which

explains the decrease of the axial velocity in the beginning of the pipe.

Simulation in a Small-Diameter Pipe

Similar tests were conducted with another set of domains. The principle is the

same, but the diameter of the tubes was much smaller. For the geometry described on

Figure 4-8, the tubes were 2 cm outer diameter. This means that ro=l cm and ri = 0.5 cm.

The first simulations were conducted without fluid heating. For these calculations, the

mesh was composed of a 200 x 50 orthogonal grid.

If no power is added to the fluid

In this case where the power density for the heat source was zero, three different

inlet velocity profiles were chosen. The outer channel axial inlet velocity was set equal to

2.0 m/s and three tests were done for inner channel velocity: 2.5 m/s, 4.0 m/s and 6.0 m/s.

The Figures 4-19 and 4-20 show the prediction of axial and radial velocity for this

experiment. The behavior of the flow is similar to the case of a large diameter pipe.

However, as seen on Figure 4-21, cross-flow is lower for small diameter tubes. This is

explained by the wall effect that acts to push the flow inwards. Small diameter tubes are

more sensitive to this phenomenon. But this phenomenon was observed in a region close

to the inlet of the tube, as farther, the same behavior is observed. RELAP5 tests have also

been performed on the same geometry under the same conditions. Figure 4-22 shows the

cross-flow velocity prediction for the three types of unheated calculations.

Figure 4-19: Radial velocity in a small size tube in an adiabatic tlow.
Vi=4.0 m/s and Vo=2.0 m/s.

: Axial velocity in a small size tube in an a
Vi=4.0 m/s and Vo=2.0 m/s.


0 ..
-0.005 0 0.005 0.01 0.015 0.02 0.025
radial velocity (m/s)

Figure 4-21: Cross-flow velocity versus distance to diameter ratio for large and small
diameter tubes.

0.1 -
'" 0.06 -* case Vi=2.5 m/s.
0.04 --case Vi=4.0 m/s.
0.02 -A-case Vi=6.0 m/s.

0.E+00 1.E-03 2.E-03 3.E-03 4.E-03 5.E-03 6.E-03 7.E-03 8.E-03
velocity (m/s)

Figure 4-22: Cross-flow velocity for unheated small diameter channels (RELAP5

Cross-flow is increasing as inner channel velocity is increasing. This is the

general behavior expected from previous cases. However, the same comments can be

done concerning axial velocity. The RELAP5 predictions of axial velocity evolution

show again that mixing is not occurring as fast as FLUENT predicts it. This is a problem

in RELAP5 modelisation of cross-flow that has been discussed previously.

If power is added to the fluid

The effect of heating on flow evolution in the channel is obvious for small

diameter channels too. Two different cases were run to observe the effect of heating. The

first one was run with a very high power density, 600 MW/m3, which makes a total

power of 10 kW in the first channel, and no power in the second channel. Inlet velocities

were taken equal to 2.0 m/s in both channels, and inlet temperature is still 555 K for a

pressure of 15.68 MPa. And the second case has the same initial conditions, but power

density is more realistic: the inner channel has a 160 MW/m3 (2.7 kW for the total

channel) and the outer channel has a 20 MW/m3 (0.94 kW for the whole channel). As

heat is transferred to the two channels, with a power density higher in the center than in

the outer channel, the axial velocity increases in both channels. This phenomenon is

calculated by RELAP5 and FLUENT. As inlet velocity is the same across the inlet area,

the heating process and temperature increase induces cross-flow. Also cross-flow is more

important in the case where power density difference in the tube is more important. As

seen on Figure 4-23, cross-flow is much larger when the difference of power density

between the two channels is the largest.

Steam flow

Previous cases were considering pure liquid flows. This case shows an example

with superheated steam. No heating occurs along the path, but initial velocities are set

much higher. The inner channel initial velocity is equal to 40.0 m/s and the outer channel

initial velocity is 20.0 m/s. A pressure of 10 bar was used (106 Pa), and the inlet

temperature is 520 K. Calculations done with RELAP5 show that cross-flow velocity is

higher in the case of steam flow, and the mixing is faster (Figure 4-24).

--Pl=10 kW (600 MW/m3) and P2=0

-- P1=2.7 kW (160 MW/m3) and P2=0.94 kW (20 MW/m3)



E 0.12

8 0.08


0.E+00 1.E-04 2.E-04 3.E-04 4.E-04 5.E-04 6.E-04 7.E-04
velocity (m/s)

Figure 4-23: Cross-flow velocity between heated small diameter channels.

-- outer channel axial velocity -- inner channel axial velocity
oue velocity-----------------------
8 2

15 20 25 30 35 40 45
steam velocity (m/s)

Figure 4-24: Axial velocity distribution for superheated steam.

Conclusions and Recommendations

The study showed that cross-flow is an important concern in modeling nuclear

reactor cores. A series of calculations under different flow conditions have demonstrated

that the magnitude of cross-flow varies with inlet flow and temperature conditions as well

as the velocity and temperature axial gradients. It was shown that both velocity and

temperature gradients have significant effect on the magnitude of the cross-flow.

Results of the study have shown the lack of sufficient accuracy of the RELAP5

cross-flow model even for cases involving single-phase flow and simple geometry.

Although high resolution CFD allows for high accuracy prediction of cross-flow, it

cannot be used to model large portion of a reactor core. The number of computational

mesh points to resolve a full reactor core would be in excess of billions that is far beyond

the current computational capability of the most advanced computer systems.

Results of the study have shown the limitations of RELAP5 to model multi-

dimensional flow and heat transfer problems. These limitations were illustrated by the

calculation of the radial velocity in the cross-flow junctions. The close-channel treatment

of the flow in the reactor core, which completely ignores the cross-flow, is a less

desirable solution. In RELAP5 model because the advection terms are neglected in the

momentum equation, there is no transport of axial momentum in the radial direction. As a

result, using a simple loss factor RELAP5 treats losses in the transverse direction. This is

obviously a crude approximate

Along with the cross-flow junction comes the problem of stress between fluid

layers. In a problem like the one treated in this study, shear stress between streams of

fluid is important. Especially when axial velocities are different, viscosity phenomenon

will play an important role in defining the shape of the flow. In reactor cores, special

attention is needed to model cross-flow in presence of fuel rod bundles that partially

obstruct radial flow motion. The primary concern is to find a proper loss factor that

would be used as a wall loss factor to model friction between vertical layers of flow in

different channels. Another issue is the proper modeling of the turbulence effect on cross-

flow calculation. In nuclear reactor cores, cross-flow is also influenced by all the

equipment and core components that are used to divert the flow from its initial path to

provide better mixing and homogenization of flow parameters.

The study was conducted for a purely single-phase flow where the contribution of

cross-flow under uniform flow and thermal conditions is insignificant. In two-phase flow

systems, the relative contribution of cross-flow is much more significant. Further studies

are needed to evaluate the RELAP5 modeling capability for analysis of two-phase flow

and heat transfer in open channels of light water reactors. The use of high resolution CFD

model proved valuable in calculating cross-flow across parallel channels. Future studies

should address both cross-flow modeling and empirical relations that are used to account

for loss of momentum in transverse direction in one-dimensional codes like RELAP5.

In summary, RELAP5 is a widely used thermalhydraulic code, but this study has

shown that cross-flow results need to be interpreted with caution because the model used

is not accurate. Cross-flow modeling in RELAP5 has to be improved to consider actual

three-dimensional effects. Models could be developed using CFD calculations.


Table A-1: Hydrodynamic components used to build the reactor model.
Component Number Component type Feature

Single junction
Single junction
Single junction
Single junction
Single volume
Single volume
Single volume

Upper annulus
Inlet annulus
Lower head
Lower plenum
Core channel 1
Core channel 3
Outer core channel
Core channel 2
Core channel 4
Core bypass
Connect 111MM to 114MM
Connect 112MM to 115MM
Connect 114MM to 112MM
Connect 115MM to 113MM
Upper plenum center 1
Upper plenum center 2
Upper plenum center 3
Upper plenum center 4
Upper plenum middle 1
Upper plenum middle 2
Upper plenum middle 3
Upper plenum middle 4
Upper plenum outer 1
Upper plenum outer 2
Upper plenum outer 3
Upper plenum outer 4
Control assembly housing
Control assembly housing
Control assembly housing
Upper head


Table A-2: Hydrodynamic components in the steam generator.
Component Number Component type Feature
X66 Time dep. volume Sg feedwater
X67 Time dep. Junction Sg feedwater
X70 Single volume Sg upper downcommer
X72 Branch Sg inlet downcommer
X74 Pipe Sg downcommer
X75 Single junction Sg boiler
X76 Pipe Sg boiler
X78 Separator Sg swirl vanes
X80 Branch Sg dryer
X82 Single volume Sg dome
X84 Branch Sg steam line
X85 Valve Sg steam outlet
X86 Time dep. Volume Steam outlet
X87 Valve Sg porv
X88 Time dep. Volume Safety container
X89 Valve Sg safety valve
X90 Time dep. Volume Safety container

Table A-3: Hydrodynamic components of the primary loop.
Component Number Component type Feature
XOO Pipe Hot leg vessel outlet
X02 Branch Hot leg pressurizer
X03 Single junction Hot leg stop
X04 Single volume Hot leg sg inlet
X06 Branch Sg primary coolant inlet
X08 Pipe Sg tube sheet
X10 Branch Sg primary coolant outlet
X12 Pipe Pump suction pipe
X14 Pump Primary loop pump
X16 Single volume Cold leg pump outlet
X17 Single junction Cold leg stop
X18 Single volume Cold leg ECC1
X20 Branch Cold leg ECC2
X22 Single volume Cold leg vessel inlet


Table A-4: Hydrodynamic components of the pressurizer.
Component Number Component type Feature

Single junction
Single volume

Pressurizer dome
Pressurizer tank
Junction hot leg/surge line
Surge line
Pressurizer porv
Pressurizer safety valve


Table B-l: Heat structures in the core and the steam generator.

Heat struct. Number


X001 Hot leg piping
X061 s.g. inlet plenum
X062 s.g. partition plate
X081 s.g. tubes
X082 s.g. tubesheet
X101 s.g. outlet plenum
X121 s.g. pump suction plenum
X161 Cold leg piping
X701 s.g. upper shell
X721 s.g. middle shell
X741 s.g. shell base
X742 s.g. lower shell
X761 s.g. lower wrapper
X781 s.g. upper wall/swirl vane
X821 s.g. shell top
1001 Upper Head
1002 Upper vessel wall
1003 Middle vessel wall
1004 Lower vessel wall
1041 Thermal shield
1061 Lower head structures
1071 Lower support plates
1081 Lower plenum structures
1151 Core baffle
1201 Core barrel
1501 Center channel upper plenum structures
1601 Middle channel upper plenum structures
1701 Outer channel upper plenum structures
1811 Center channel control housing
1821 Middle channel control housing
1831 Outer channel control housing
1851 Upper support plate
1901 Upper head structures


This appendix presents the RELAP5 inputs used to model the single pipe with taking into
account the cross-flow phenomenon.

=two channel test
100 new stdy-st
102 si si
201 400.0 1.Oe-6 0.05 3 50 200 200

506 time 0 gt null 0 1000.0 1 *forpower

0880000 botin tmdpvol
0880101 7.854e-5 1.Oe+02 0.0 0 90.0 1.Oe+2 0 0 0
0880201 0.0 1.56821e+7 555.0

0980000 botout tmdpvol
0980101 2.356e-4 1.Oe+02 0.0 0 90.0 1.Oe+2 0 0 0
0980200 3
0980201 0.0 1.56821e+7 555.0

0890000 junin tmdpjun
0890101 088010002 100010001 7.854e-5
0890201 0.0 2.0 2.0 0.0
0890202 1.Oe+06 2.0 2.0 0.0

0990101 098010002 110010001 2.356e-4
0990201 0.0 2.0 2.0 0.0
0990202 1.Oe+6 2.0 2.0 0.0

*inlet channel
1000000 channel annulus
1000001 20
1000101 7.854e-5 20
1000301 0.01 20
1000601 90.0 20
1000801 0 0.0120
1001001 020
1001101 019
1001201 3 1.58e+7 550.0 0 0 0 20
1001301 2.0 2.0 0.0 19

*outlet channel
1100000 channel annulus
1100001 20
1100101 2.356e-4 20
1100301 0.01 20
1100601 90.0 20
1100801 0 0.0120
1101001 020
1101101 019
1101201 3 1.58e+7 550.0 0 0 0 20
1101301 2.0 2.0 0.0 19

*" junction
1120000 crfl mtpljun
1120001 200
1120011 100010004 110010003 3.1416e-4 0.0 0.0 00003
+ 10000 10000 0 20
1121011 0.00.020

*sngljun in
1 '* '* '* '*' sngljun
1200101 100200002 130000000 .854e-5 0.0 0.0 0 1.0 1.0
1200201 0 2.0 0

*sngljun out
12 . '. ' sngljun
1210101 110200002 131000000 2.356e-4 0.0 0.00 1.0 1.0
1210201 0 2.00

*outlet-in tmdpvol
1300000 outin tmdpvol
1300101 7.854e-5 1.Oe+2 0 090.0 1.Oe+2 0 0 0010
1300201 0.0 1.568e+7 555.0

*outlet-out tmdpvol
1310000 outout tmdpvol
1310101 2.356e-4 1.Oe+2 0 0 90.0 1.0e+2 0 0 0010
1310201 0.0 1.568e+7 555.0

* heat structure inner channel
11000000 20 2 2 0 0.005
11000100 0 1
11000101 1 0.006
11000201 51
11000301 1.01
11000400 0
11000401 555.02
11000501 10001000010000 1 1 0.01 20
11000601 0 0 0 1 0.01 20
11000701 10 0.05 0 0 20
11000801 0.01 10.0 10.0 0.0 0.0 0.0 0.0 1.0 20

* heat structure outer channel
11100000 20 2 2 0 0.0087
11100100 0 1
11100101 1 0.0097
11100201 51
11100301 1.01
11100400 0
11100401 555.02

11100501 11001000010000 1 1 0.01 20
11100601 0 0 0 1 0.01 20
11100701 11 0.05 0 0 20
11100801 0.015 10.0 10.0 0.0 0.0 0.0 0.0 1.0 20

*conduction between the channels
11120000 20 2 2 0 0.005
111201011 0.006
11120201 0041
11120301 0.01
11120499 555.0 2
11120501 100010000 10000 1 1 0.01 20
11120601 110010000 10000 1 1 0.01 20
11120701 00.00020
11120801 0.01 10.0 10.0 0.0 0.0 0.0 0.01.0 20
11120901 0.03 10.0 10.0 0.0 0.0 0.0 0.0 1.0 20

* general tables
*thermal properties for water

5 ,' ,','-- ''. .. 1 -1
20100401 523.0 0.635
20100402 548.0 0.6034
20100403 573.0 0.5617
20100404 598.0 0.5086

20100451 2.24e+6
20100452 2.312e+6
20100453 2.341e+6
20100454 2.313e+6
* stainless steel
20100500 s-steel

* table of time, core power
20201000 power 506
20201001 0.0 10000.01.0e+610000.0
20201100 power 506
20201101 0.0 0.0 .Oe+6 0.0


1. W. WANG, "CFD Mixing analysis of vortex generator jets injected into confined
cross-flow in rectangular duct," Tuijin Jishu/Journal of Propulsion Technology
19, 2 (1998).

GARTSHORE, "Multiple jets in a cross-flow: Detailed measurements and
numerical simulations," American Society of Mechanical Engineers (Paper)

3. V. KASSERA, "Three dimensional CFD-analysis of tube vibrations induced by
cross-flow," 4 International Symposium on Fluid-Structure Interactions,
Aeroelasticity, Flow-Induced Vibrations and Voice ASME, Aerospace Division
AD, 53-2 (1997).

4. SCADP/RELAP5 Development Team, "SCDAP/RELAP5/MOD3.2 Code Manual
Volume I: SCDAP RELAP5 Interface Theory," NUREG/CR-6150, INEL-
96/0422, INEL, Idaho Falls (October 1997).

5. RELAP5 Code Development Team, "RELAP5/MOD3 Code Manual Volume I:
Code Structure, System Models and Solution Methods," NUREG/CR-5535,
INEL-95/0174, INEL, Idaho Falls (June 1995).

6. SCADP/RELAP5 Development Team, "SCDAP/RELAP5/MOD3.2 Code Manual
Volume III: User's Guide and Input Manual-Appendix A," NUREG/CR-6150,
INEL-96/0422, INEL, Idaho Falls (November 1997).

7. "FLUENT 5 User's Guide, Volume 2," Fluent Inc., Lebanon (July 1998).

8. "FLUENT 5 User's Guide, Volume 3," Fluent Inc., Lebanon (July 1998).

9. "FLUENT 5 User's Guide, Volume 1," Fluent Inc., Lebanon (July 1998).

10. "FLUENT 5.5 UDF User's Guide," Fluent Inc., Lebanon (September 2000).

11. "FLUENT 5 User's Guide, Volume 4," Fluent Inc., Lebanon (July 1998).

12. W. RODI, Turbulence models and their application in hydraulics, A state-of-the-
art review, IAHR/AIRH, Rotterdam (1993).

13. W.K. GEORGE, R. ARNDT, Advances in Turbulences, Hemisphere Publishing
International, New York (1989).

14. S.V. PATANKAR,. D.B. SPALDING, Heat and Mass Transfer in Boundary
Layers, Applied Mathematics and Mechanics, Vol 15, Academic Press, New York

15. H. SCHLICHTING, Boundary Layer Theory, McGraw Hill, New York (1969).

16. B.E. LAUNDER, D.B. SPALDING, Lectures in Mathematical Models of
Turbulence, Academic Press, London, England (1972).

17. D. CHOUDHURY, "Introduction to the Renormalization Group Method and
Turbulence Modeling," Technical Memorandum TM-107, Fluent Inc., Lebanon

18. T. -H. SHIH, W. W. LIOU, A. SHABBIR and J. ZHU, "A new k-epsilon eddy
viscosity model for high Reynolds number turbulent flows Model development
and validation," Computer Fluids 24, 3 (1995).

19. P. SPALART, S. ALLMARAS, "A one-equation turbulence model for
aerodynamic flows," Technical Report AIAA-92-0439, American Institute of
Aeronautics and Astronautics, New York (1992).

20. U. GRIGULL, J. STRAUB, P. SCHIEBNER, Steam Tables in SI Units, Springer-
Verlag, New York (1990).


Vincent Roux was born on May 25, 1977, in Perigueux, France. He graduated

from Girault de Borneil High School in France. He entered the National School for

Physics in Grenoble in the fall of 1997. As part of an academic exchange, he entered the

University of Florida in the fall of 1999 as a graduate student, and obtained his French

engineer diploma in Summer 2000. Since then, he has been pursuing a Master of Science

degree in nuclear engineering sciences while working at INSPI as a graduate research


University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs