EVALUATION OF RELAP5 REACTOR CORE MODELING CAPABILITY
By
VINCENT J.P. ROUX
A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
UNIVERSITY OF FLORIDA
2001
Copyright 2001
by
Vincent J.P. Roux
To Lucie
ACKNOWLEDGMENTS
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
work.
TABLE OF CONTENTS
page
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
Crossflow 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
3 PROTOCOL OF EVALUATION OF CROSSFLOW MODEL USING HIGH
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
v
Turbulence Modeling in FLUENT............. ........................ .............. 42
Oneequation SpalartAllmaras model ............. ..... .....................................46
Tw oequation kepsilon m odel ................................... ............ ....... .............. 46
Defining a CFD Method for Evaluating RELAP5 CrossFlow 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
4 CROSSFLOW CHARACTERIZATION AND EVALUATION OF RELAP5
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
Crossflow Between Two Adjacent Channels: FLUENT and RELAP5 Results..........64
Simulation in a LargeDiam eter Vertical Pipe ......... ..................... .................. 65
If no power is added to the fluid .......................... .... ................... 65
If power is added to the fluid......... ......................................... .............. 70
Simulation in a SmallDiameter 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
APPENDICES
A SUMMARY OF SCDAP/RELAP5 HYDRODYNAMIC COMPONENTS............... 80
B RELAP5/SCDAP HEAT STRUCTURES FOR THE CORE AND STEAM
GENERATOR ................................................................... .......... ........ 83
C RELAP5 INPUTS FOR THE TWOCHANNELS SYSTEM............................... 84
LIST OF REFERENCES... ...................................................... ...... .............. 86
BIOGRAPHICAL SKETCH .......................................................... ........... ..... 88
LIST OF TABLES
Table Page
21 Fuel rod pow er distribution ........................ ............. .................. .............. 27
41 Comparison of axial junction mass flow for a core with and without crossflow.....56
Ai Hydrodynamic components used to build the reactor model.................. ..............80
A2 Hydrodynamic components in the steam generator................................................. 81
A3 Hydrodynamic components of the primary loop................................. ........... 81
A4 Hydrodynamic components of the pressurizer .................................................... 82
Bl Heat structures in the core and the steam generator. ............. ...........................83
LIST OF FIGURES
Figure Page
21 Difference equation nodalization schematic. .................................. .... ............... 9
22 RELAP5 nodalization diagram of the Surry reactor core............ ................ 17
23 RELAP5 nodalization diagram of a Utube steam generator .............................20
24 RELAP5 nodalization diagram of the pressurizer on the hot leg ..........................23
25 RELAP5 nodalization diagram of the primary cooling system tubing.....................25
26 Bundle matrix for the first channel: 1 is a fuel rod 2 is a control rod................ 28
27 Evolution of the temperature distribution during the steady state run........ ........ 30
28 Axial coolant temperature distribution the steady state............... ................31
29 Radial fuel temperature profile at the sixth vertical node (height=2.012 m)............32
31 Mixinglength distribution in wall boundary layers.................... ...............45
32 Basic geometry studied with RELAP5 and FLUENT...........................................49
33 FLUENT calculation domain of the system...........................................................49
34 Nodalization diagram of the equivalent RELAP5 model., ..................................... 54
41 Crossflow magnitude versus height in the core for a total power of 2443MW.......57
42 Crossflow mass flux versus core cell height..................................................... 58
43 Ratio of crossflow to axial flow for different power levels and for different
height in the core.................................... ........................................ 59
44 Power transient simulated in the transient test............................................. 61
45 Crossflow magnitude versus height in the core for a total power of 3500MW.......61
46 Crossflow to inlet flow ratio for a 1000MWth case. Crossflow is from Channel 1
to 2 ..................................... .................................... ........... 62
47 Crossflow to inlet flow ratio for a 2443MWth case. Crossflow is from Channel 1
to 2 .................................... ................... ................ ........... 63
48 Twosub channel pipe studied with FLUENT and RELAP5 ...............................64
49 Inlet velocity profiles used in crossflow simulation ...............................................65
410 Axial evolution of the axial velocity profile. z is the height in meters................... 66
411 Radial (Vvelocity) and axial (Uvelocity) velocity distribution at the entrance of
the tube for Vi=6.0 m/s and Vo=2.0 m/s...........................................................67
412 Radial (Vvelocity) and axial (Uvelocity) velocity distribution at the entrance of
the tube for Vi=4.0 m/s and Vo=2.0 m/s...........................................................67
413 Comparison of crossflow prediction between FLUENT and RELAP for two
sim ulations............. ................. .................................................. 68
414 Evolution of the axial velocity for the two RELAP5 axial channels along the pipe,
for the two inlet boundary conditions (a) and (b). ....................... ...............69
415 Temperature distribution over the total length of a 20.0 m pipe heated in the
central channel at 66.1 M W /m3. .............................................. .................. 70
416 Radial velocity for a constant inlet velocity and a heated central channel. ..............71
417 Crossflow velocity at the artificial interface between the two channels. ...............72
418 Axial velocity in Channel 1 (heated) and Channel 2 (nonheated) ..........................72
419 Radial velocity in a small size tube in an adiabatic flow. Vi=4.0 m/s and
V o=2.0 m /s............. .. ................. .............. .. ......... .... ........ 74
420 Axial velocity in a small size tube in an adiabatic flow. Vi=4.0 m/s and
V o=2.0 m /s. ............ .. ................. .............. .. ......... .... ............ 74
421 Crossflow velocity versus distance to diameter ratio for large and small diameter
tubes. ............. .............. .. .. ........ .. ................. ........... 75
422 Crossflow velocity for unheated small diameter channels (RELAP5 calculations).75
423 Crossflow velocity between heated small diameter channels..............................77
424 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
EVALUATION OF RELAP5 REACTOR CORE MODELING CAPABILITY
By
Vincent J.P. Roux
August 2001
Chairman: Professor Samim Anghaie
Major Department: Nuclear and Radiological Engineering
RELAP5 is a onedimensional reactorsystem simulation code with additional
crossflow calculation capability to include two and threedimensional 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 threeloop 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 crossflow
model. Several inlet flow rates and core power distributions are considered and modeled.
Results of the analysis showed the significant contribution of crossflow in overall
temperature and flow distributions in the core. Results of the study also showed that the
RELAP5 predictions of crossflow, at least for singlephase cases, are not consistent with
the theory.
To evaluate the accuracy of RELAP5 crossflow 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 highresolution FLUENT model for a pair of parallel reactor core
channels. Two models were developed for FLUENT calculation of crossflow. The first
model is a simple tube with axisymmetric nonuniform 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 crossflow
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 kepsilon 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 threedimensional flow phenomena using the crossflow 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 crossflow model for predicting complicated two and
threedimensional flow phenomena.
CHAPTER 1
INTRODUCTION
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 twophase 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 incore 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
performance.
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 crossflow phenomenon. Crossflow is
a wellknown phenomenon, and extensive literature is available to show experimental
observations and attempts to model this specific type of flow. By definition, the term
"crossflow" 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 crossflow cooling. However, the designation crossflow is more
commonly adapted to the flow of a jet mixing in another flow flowing in another
direction. But these are crossflows generated by geometry conditions, as analyzed by
Wang (1). Crossflow 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 crossflow jet
mixing in an axial stream. In a nuclear reactor core, crossflow 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 twophase 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), crossflow across the core is important. Crossflow in a
nuclear reactor core occurs between the different channels (or even subchannels). 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 crossflow 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 crossflow 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 crossflow phenomenon, and a comparison of the results was
performed. Finally, the calculations showed that there was a difference between RELAP5
and FLUENT crossflow calculations. This thesis presents an actual simulation of a
reactor, enhancing occurrence of crossflow inside the core, and discusses the modeling
of crossflow in RELAP5, based on highresolution CFD simulations.
CHAPTER 2
SYSTEM MODELING WITH RELAP5
In the modeling of crossflow 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 crossflow 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
crossflow 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 highresolution 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
calculated.
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 nonnuclear 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 nonequilibrium, sixequation and twofluid model. RELAP5 allows
choosing an arbitrary number of volumes, junctions and surfaces. RELAP5 is a one
dimensional code; however, it can handle multipledimensional nodalization by using
crossflow 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 noncondensable 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 welldefined 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 onedimension 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 interrelated 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
7
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
Numerics
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 21
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 illposed 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 steadystate or selfinitialization 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 semiimplicit scheme is to replace the system of partial
differential equations with a system of finitedifference equations partially implicit in
time. For the semiimplicit 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, iuf, 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
Vg
SCL
K L
Vf
j1 J j+1
Momentum control
volume or cell
Figure 21: 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 twoequation 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
velocity.
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 orderfive 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 )
Au'
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 + f4, 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 onedimensional 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 twodimensional 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 semiimplicit 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 semiimplicit 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 semiimplicit 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 21) as in the semiimplicit 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 semiimplicit 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
semiimplicit 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).
Crossflow modeling
As RELAP5 is a onedimension 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 crossflow between reactor core channels, simplified crossflow
momentum equations are used to couple two adjacent channels linked with a crossflow
junction. In the momentum equation, the transverse momentum advection terms are
neglected, meaning that there is no transport of xdirection momentum due to flow in the
transverse direction. Crossflow 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 crossflow 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 3Loop 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.
Core
Figure 22 shows the nodalization diagram of the core. This is a fivechannel 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 3channel model ([151152153154], [161162163164] and [171172173
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 crossflow 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 branchtype 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 crossflow junctions:
Components 140149 connect the center core Channel 1 and center core
Channel 2.
Components 120129 connect the center core Channel 2 and the middle
core Channel 3.
17
Components 8089 connect the middle core Channel 3 and the middle core
Channel 4.
Components 130139 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
104
118
108
Figure 22: RELAP5 nodalization diagram of the Surry reactor core.
Figure 22: 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 crossflow, the four sets of crossflow 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 crossflow 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 upperplenum (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 151161 and 152162 and 7.0 for the junctions 161171 and 162172). The
inverted hat of the upper head has been represented using branch components (153154,
163164 and 173174). 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 subchannel, 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 23 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 Utubes 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 nonradioactive 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 Utubes 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.
S282
280 safety valve
278 286 288 290
270
Separatoi i
266 272 steam
L outlet
Steam Generator
Feedwater
274 276
d Dwncomrr er
Figure 23: RELAP5 nodalization diagram of a Utube 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 builtin 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 servovalve (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 s1 whereas the
safety valve opens at a rate 20.0 s1. 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).
Pressurizer
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 24 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
440
pressurizer
441
\ 443
Figure 24: 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 branchtype 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 s1 (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 25 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 builtin 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 2cells 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 Utube
geometry (as shown on Figure 25).
valve 4 44 4
4
440
441
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 reactorcooling 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 21 shows the power fraction for each fuel rod
component and the number of rods associated with this component.
Table 21: 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 26.
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 26: 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 "steadystate" 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 steadystate 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 (104105 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.9E5, 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.3099E08, 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 27 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
60 2nd cell
595 11102
! 3rd cell :
590 11103
E 4th cell :
585 11104
5th cell:
580 11105
66th cell:
575 11106
S 7u 7th cell
S570 11107
565 X8th cell
11108
560  9th cell:
11109
555 WM O10th cell
11110
0 50 100 150
"time advancement"
Figure 27: 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 28.
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 29. These shapes may be obtained theoretically.
These types of calculation are the major interest of RELAP5. Depending on whether the
core model includes crossflow 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.
605
600 
595
590
585
580
X 575 first channel
570 second channel
565 Athird channel
X fourth channel
560 fith channel
555 
0 0.2 0.4 0.6 0.8 1
normalized axial node elevation
Figure 28: Axial coolant temperature distribution the steady state.
The purpose of this chapter was to describe the first tool used in the investigation
of the crossflow 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.
1600
.1400 >
1200 
1000 
S800
600
400
r200
0
0 0.005 radius (m)01o 0.015
Figure 29: 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 crossflow 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.
CHAPTER 3
PROTOCOL OF EVALUATION OF CROSSFLOW MODEL USING
HIGH RESOLUTION COMPUTATIONAL FLUID DYNAMICS
About FLUENT
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.
Crossflow phenomenon in a nuclear reactor core was studied in different ways.
The effect of crossflow 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 crossflow
junction. The code solves the momentum equation between the two adjacent cells to
calculate crossflow 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 crossflow 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 NavierStokes 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
crossflow 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 NavierStokes 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 crossflow
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
twodimension 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 twodimension space, which calculates on
a twodimension 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 twodimension 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))= kak 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 twodimension 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
and
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 socalled
"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 ks 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 userfriendly 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 controlvolume 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 colocated 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 massflow 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
P
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
2
terms. As the average of the perturbation term is zero and as one takes a timeaveraged
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 meanvelocity 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 zeroequation model is also called Prandtl
mixinglength 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 31. 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 mixinglength distribution is well
described by Nikuradse's formula (15):
lm/R = 0.140.08(1y/R)20.06(1y/R)4
tm
Xx Qs5 yI6 1
Figure 31: Mixinglength 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 crossflow 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: SpalartAllmaras model, ke 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
ke models, which are respectively oneequation and twoequation models and use a
Boussinesq eddy viscosity hypothesis for the divergence of the Reynolds constraints
tensor. SpalartAllmaras model could also be used without any problem
Oneequation SpalartAllmaras 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 oneequation 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 SpalartAllmaras 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 SpalartAllmaras model is done in the FLUENT
user manual (7) and in the original paper (19).
Twoequation kepsilon model.
This ke model is used as a better alternative to the oneequation 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 ke 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 ks model are:
Dk a 9, Sk
P =a[It Ok + G, + Gb PYM
Dt ax, k 2J+x,
Ou
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 ke model without
changing the default constants. A detailed description of other twoequation 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
avg
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 CrossFlow Model
RELAP5 uses a coarse integral method to solve the conservation equations. The
crossflow 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 crossflow 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 crossflow calculation, as velocity is a
49
parameter that cause crossflow. Figure 32 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 crossflow occurrence. The velocity
profile is radius dependent, with the central part of the flow at a higher velocity.
inlet velocity profile
Figure 32: Basic geometry studied with RELAP5 and FLUENT.
FLUENT Model
The FLUENT model is supposed to be as close as possible to the real geometry
and operating conditions. Figure 33 shows the calculation domain representing the pipe
of the Figure 32.
p4 wall P3
Vo
ro
inlet outlet
iinl
P L axs 92
Figure 33: 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 crossflow 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 33, 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 2D 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 onedimensional 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
developed.
Wall boundaries. The centerline boundary was represented by an axis type
boundary condition. The outer wall ("wall" on Figure 33) 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 104, and a roughness constant approximately equal to 0.5 are acceptable
numbers.
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 33, 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 crossflow phenomenon
inside the pipe. For this, the pipe was divided axially into two pipes. The Figure 34
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 33, 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 (ro22ri). 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 crossflow 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
K.
54
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 crossflow 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
11002
111 0 I 1
11001 10001
11201
i99 819
98 10018
Figure 34: Nodalization diagram of the equivalent RELAP5 model.,
The results of the equivalent simulations with FLUENT and RELAP5 are
presented in Chapter 4.
CHAPTER 4
CROSSFLOW CHARACTERIZATION AND EVALUATION OF RELAP5 MODEL
This chapter will present the results from different simulations with FLUENT and
RELAP5. For each case, crossflow has been calculated. First, RELAP5 has been run for
a steady state case and a power surge transient to characterize crossflow phenomenon in
a real case simulation. And then, the attention is focused on the crossflow phenomenon
itself by using the two models described in Chapter 3.
Surry Reactor Calculation for a Steady State and a Simulated Power Surge Transient
Crossflow is predicted in both the steady state and transient cases. The results
will show that the impact of taking into account crossflow 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 crossflow is
considered or not. Without crossflow 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 104 for a core without crossflow. Table 41
shows the axial mass flow in the five hot core channels for the case of a purely axial core
or a core with crossflow junctions between the axial channels.
Table 41: Comparison of axial junction mass flow for a core with
and without crossflow.
Axial mass flow (kg/s) for core
with crossflow
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 crossflow
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 closedchannels core is close to the one at the
different axial junctions along the core allowing crossflow, but still a difference is
noticeable. For a given channel, the closedchannel value is roughly an average value of
the axial mass flow in the open channel case. The difference comes from the crossflow
consideration. However, as can be seen on Figure 41, in the steady state conditions,
crossflow is much smaller than the main axial mass flow, which is expected. This figure
shows the magnitude of crossflow between the five different channels that compose the
core. The sign of the calculated crossflow represents its direction. The way the channels
are numbered in the RELAP5 inputs, a positive value for crossflow 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.
8
6
 from channel 1 to 2
5 from channel 2 to 3
From channel 3 to 4
4 X from channel 4 to 5
2
4 2 0 2 4 6
mass cross flow (kg/s)
Figure 41: Crossflow magnitude versus height in the core
for a total power of 2443MW.
However, to have a better idea of the relative importance of crossflow 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 crossflow versus cell height in the core. This one is calculated by dividing
the mass flow by the crossflow junction area. Actually, the direction of the crossflow is
more interesting. According to the sign of the mass flow, and looking at the setup of the
junctions, there is a crossflow going outwards for all central channels. The outer junction
58
is a little bit different. This kind of prediction is surprising as different parameters can be
considered:
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 crossflow 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 42: Crossflow mass flux versus core cell height.
16 18
Now different power levels can be considered. RELAP5 will predict different
crossflow 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 crossflow value was calculated, and normalized to the inlet
mass flow in the reactor. Figure 43 shows that in a general manner, crossflow
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
A1st cell ch. 3 to 4
X 2nd cell ch. 3 to 4
3rd cell ch. 3 to4 4
o G4th cell ch. 3 to 4 C
u 5fth cell ch. 3 to 4
d 8th cell ch. 3 to 4
o
LO C
+ o
o o
o 2443 2000 Power (MW) 1500 1000 o
Figure 43: Ratio of crossflow to axial flow for different power levels and for different
height in the core.
W I
6 2443 2000 Power (MW) 1500 1000 (?
Figure 43: Ratio of crossflow 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
crossflow 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, crossflow 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 44. 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 44: 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 crossflow magnitude. As mass flow in the core is constant, fluid temperature
increases when increasing power in the reactor.
Crossflow magnitude is larger for a greater power level. As can be seen in Figure
45, and also as in the previous part in Figure 43, the mass flux through the crossflow
junctions is larger when power is larger (whatever the sign of the crossflow is).
10
9
8
*channel 1 to 2
04
0 *channel 2 to 3
3
 channel 3 to 4
2 channel 4 to 5
0 3 6 9 1 15 18
mass cross flow (k gs)
Figure 45: Crossflow magnitude versus height in the core for a total power of 3500MW.
To illustrate the transient phase of the crossflow, focus will be drawn on the
central channel. The crossflow ratio is defined as the ratio of crossflow 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 crossflow is normalized to the inlet flow
rate, but instead, parameters like temperature will affect the crossflow variation. Figure
46 and Figure 47 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 crossflow 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
conditions.
8th cell
C"'
O2H+  c  w first cell
0 9 5 5 ?rfl ?C 97f4th cell97.
0
transient time (s)
Figure 46: Crossflow to inlet flow ratio for a 1000MWth case. Crossflow is from the
Channel 1 to 2.
The figure 47 shows that for the seventh and the eighth cell, crossflow takes
longer to stabilize and has also more oscillations in the transient phase. However, for the
bottomcore crossflow junctions, the steady state is reached quite smoothly, but it still
requires some time before the level of crossflow stabilizes. As the temperature in the
core is increasing with certain inertia, it is understandable that crossflow, 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 crossflow junctions would be
disturbed, since it picks up perturbations from the lower part of the core.
C
8 7th
Figure 47 Crossflow to inlet flow ratio for a 2443MWth CS. COSSflOW is om
o 7th cell
uJ
o 245 250 255 260 265 270 275 280
transient time (s)
Figure 47: Crossflow to inlet flow ratio for a 2443MWth case. Crossflow is from
Channel 1 to 2.
Crossflow phenomenon is now characterized with RELAP5. An expected result
is that crossflow increases as power increases. However, the crossflow 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.
Crossflow Between Two Adjacent Channels: FLUENT and RELAP5 Results
This part presents the results of crossflow calculation in the simple system
presented Chapter 3. This model has been developed to calculate crossflow between two
axial streams of fluid. This should allow doing an evaluation of the crossflow model
used in RELAP5. FLUENT and RELAP5 were used to calculate the crossflow 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 socalled
crossflow junction. The geometry of the system is reminded in Figure 48.
inlet
Figure 48: Twosub 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 LargeDiameter 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 33.
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 49, profile (a).
velocity (m/s) velocity (m/si
6.0 
4.0
2.0 2.0 I
0 0.4 radius m) 0 0.2 0.4 radius (m)
Figure 49: Inlet velocity profiles used in crossflow 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 nonuniform 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 410, this pipe is long enough such that the flow becomes fully developed. The
evolution of axial velocity shape is displayed on Figure 410. The mixing length can be
determined as the point where the flow gets his final shape that follows a oneseventh
power law. The radial velocity is the crossflow 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 411 shows a more detailed view of the mixing process. Both
plots represent a length of 5 meters. Figure 411 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.
7.00
6.00
5.00
4.00
3.00
2.00
1.00
0.00
0.00 0.10 0.20 0.30 0.40
radius (m)
Figure 410: Axial evolution of the axial velocity profile. z is the height in meters.
z=8.3 z=6.25 z=4.2
z=1.25
z=0.8
z=O.O
I7 t '5Z .
are 41 1: Kiaalai (vvelocity) ana axial tuveloclty) velocity alsmounon
at the entrance of the tube for Vi=6.0 m/s and Vo=2.0 m/s.
Figure 412: Radial (Vvelocity) and axial (Uvelocity) 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 49, and results are shown Figure 412. With such a lower velocity, crossflow
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
412) than for an inner velocity of 6 m/s (Figure 411).
 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
5
4.5
4 [1
3.5
2
13
1 2.5 11
1.5 ]
0.5
0 El I I I ,
0.00 0.01 0.02 0.03 0.04 0.05 0.06
velocity (m/s)
Figure 413: Comparison of crossflow 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 413 shows the predictions of crossflow with two different
methods: RELAP5 and CFD with FLUENT.
There are differences between the predictions of crossflow. Results are consistent
as in average, crossflow for the case where Vi is equal to 6.0 m/s is about 2.5 times the
magnitude of crossflow 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 414.
 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
5
4.5
4
3.5
2.5
2
1.5
0.5
0
0  i i'' 
1 2 3 4 5 6 7
velocity (m/s)
Figure 414: 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. Crossflow 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 nonconstant 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 415: Temperature distribution over the total length of a 20.0 m pipe heated in the
central channel at 66.1 MW/m3.
Figure 415 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 crossflow occurring due to the heating of the fluid.
Figure 416 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 416: Radial velocity for a constant inlet velocity and a heated central channel.
On Figure 416, 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 417 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 416. This graph shows that both
FLUENT and RELAP5 predict the same behavior of the crossflow velocity at the
72
entrance of the pipe. However, there is a net difference of behavior of the predictions
after that. The following Figure 418 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.
4.5
FLUENT results
a RELAP5 results
0.025 0.02 0.015 0.01 0.005 0 0.005
fluid velocity (m/s)
Figure 417: Crossflow velocity at the artificial interface between the two channels.
5 *axial velocity in the Channel 1 (heated)
4.5 NXaxial velocity in the Chanel 2 (not heated)
4
3.5
3
Z2.5
2
E 1.5
1
0.5
0
1.6 1.8 2 2.2 2.4 2.6 2.8 3
velocity (m/s)
Figure 418: Axial velocity in Channel 1 (heated) and Channel 2 (nonheated)
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 415. 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 SmallDiameter 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 48, 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 419 and 420 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 421, crossflow 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 422 shows the
crossflow velocity prediction for the three types of unheated calculations.
Figure 419: 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.5
0 ..
0.005 0 0.005 0.01 0.015 0.02 0.025
radial velocity (m/s)
Figure 421: Crossflow velocity versus distance to diameter ratio for large and small
diameter tubes.
0.2
0.18
0.16
0.14
0.12
0.1 
0.08
'" 0.06 * case Vi=2.5 m/s.
0.04 case Vi=4.0 m/s.
0.02 Acase Vi=6.0 m/s.
0.E+00 1.E03 2.E03 3.E03 4.E03 5.E03 6.E03 7.E03 8.E03
velocity (m/s)
Figure 422: Crossflow velocity for unheated small diameter channels (RELAP5
calculations).
Crossflow 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 crossflow 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 crossflow. Also crossflow is more
important in the case where power density difference in the tube is more important. As
seen on Figure 423, crossflow 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 crossflow velocity is
higher in the case of steam flow, and the mixing is faster (Figure 424).
Pl=10 kW (600 MW/m3) and P2=0
 P1=2.7 kW (160 MW/m3) and P2=0.94 kW (20 MW/m3)
0.2
0.16
E 0.12
8 0.08
0.04
0.E+00 1.E04 2.E04 3.E04 4.E04 5.E04 6.E04 7.E04
velocity (m/s)
Figure 423: Crossflow velocity between heated small diameter channels.
 outer channel axial velocity  inner channel axial velocity
oue velocity
5
4.5
4
3.5
S3
S2.5
8 2
1.5
1
0.5
0
15 20 25 30 35 40 45
steam velocity (m/s)
Figure 424: Axial velocity distribution for superheated steam.
Conclusions and Recommendations
The study showed that crossflow is an important concern in modeling nuclear
reactor cores. A series of calculations under different flow conditions have demonstrated
that the magnitude of crossflow 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 crossflow.
Results of the study have shown the lack of sufficient accuracy of the RELAP5
crossflow model even for cases involving singlephase flow and simple geometry.
Although high resolution CFD allows for high accuracy prediction of crossflow, 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 crossflow junctions. The closechannel treatment
of the flow in the reactor core, which completely ignores the crossflow, 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 crossflow 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 crossflow 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, crossflow 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 singlephase flow where the contribution of
crossflow under uniform flow and thermal conditions is insignificant. In twophase flow
systems, the relative contribution of crossflow is much more significant. Further studies
are needed to evaluate the RELAP5 modeling capability for analysis of twophase flow
and heat transfer in open channels of light water reactors. The use of high resolution CFD
model proved valuable in calculating crossflow across parallel channels. Future studies
should address both crossflow modeling and empirical relations that are used to account
for loss of momentum in transverse direction in onedimensional codes like RELAP5.
In summary, RELAP5 is a widely used thermalhydraulic code, but this study has
shown that crossflow results need to be interpreted with caution because the model used
is not accurate. Crossflow modeling in RELAP5 has to be improved to consider actual
threedimensional effects. Models could be developed using CFD calculations.
APPENDIX A
SUMMARY OF SCDAP/RELAP5 HYDRODYNAMIC COMPONENTS
Table A1: Hydrodynamic components used to build the reactor model.
Component Number Component type Feature
Pipe
Branch
Pipe
Branch
Branch
Pipe
Pipe
Pipe
Pipe
Pipe
Pipe
Single junction
Single junction
Single junction
Single junction
Branch
Branch
Branch
Branch
Branch
Branch
Branch
Branch
Branch
Branch
Branch
Branch
Single volume
Single volume
Single volume
Branch
Upper annulus
Inlet annulus
Downcommer
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
100
102
104
106
108
111
112
113
114
115
118
14X
08X
12X
13X
151
152
153
154
161
162
163
164
171
172
173
174
181
182
183
190
Table A2: 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 A3: 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
82
Table A4: Hydrodynamic components of the pressurizer.
Component Number Component type Feature
Branch
Pipe
Single junction
Pipe
Valve
Valve
Single volume
Pressurizer dome
Pressurizer tank
Junction hot leg/surge line
Surge line
Pressurizer porv
Pressurizer safety valve
Container
APPENDIX B
RELAP5/SCDAP HEAT STRUCTURES FOR THE CORE AND
STEAM GENERATOR
Table Bl: Heat structures in the core and the steam generator.
Heat struct. Number
Feature
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
APPENDIX C
RELAP5 INPUTS FOR THE TWOCHANNELS SYSTEM
This appendix presents the RELAP5 inputs used to model the single pipe with taking into
account the crossflow phenomenon.
=two channel test
100 new stdyst
102 si si
201 400.0 1.Oe6 0.05 3 50 200 200
506 time 0 gt null 0 1000.0 1 *forpower
0880000 botin tmdpvol
0880101 7.854e5 1.Oe+02 0.0 0 90.0 1.Oe+2 0 0 0
08802003
0880201 0.0 1.56821e+7 555.0
0980000 botout tmdpvol
0980101 2.356e4 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.854e5
08902000
0890201 0.0 2.0 2.0 0.0
0890202 1.Oe+06 2.0 2.0 0.0
tmdpjun
0990101 098010002 110010001 2.356e4
09902000
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.854e5 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
10013000
1001301 2.0 2.0 0.0 19
*outlet channel
1100000 channel annulus
1100001 20
1100101 2.356e4 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
11013000
1101301 2.0 2.0 0.0 19
*" junction
1120000 crfl mtpljun
1120001 200
1120011 100010004 110010003 3.1416e4 0.0 0.0 00003
+1.01.01.0
+ 10000 10000 0 20
1121011 0.00.020
*sngljun in
1 '* '* '* '*' sngljun
1200101 100200002 130000000 .854e5 0.0 0.0 0 1.0 1.0
1.0
1200201 0 2.0 0
*sngljun out
12 . '. ' sngljun
1210101 110200002 131000000 2.356e4 0.0 0.00 1.0 1.0
1.0
1210201 0 2.00
*outletin tmdpvol
1300000 outin tmdpvol
1300101 7.854e5 1.Oe+2 0 090.0 1.Oe+2 0 0 0010
13002003
1300201 0.0 1.568e+7 555.0
*outletout tmdpvol
1310000 outout tmdpvol
1310101 2.356e4 1.Oe+2 0 0 90.0 1.0e+2 0 0 0010
13102003
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
1112010001
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 ssteel
* 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
LIST OF REFERENCES
1. W. WANG, "CFD Mixing analysis of vortex generator jets injected into confined
crossflow in rectangular duct," Tuijin Jishu/Journal of Propulsion Technology
19, 2 (1998).
2. P. AJERSCH, J.M. ZHOU, S. KETLER, M. SALCUDEAN, I. S.
GARTSHORE, "Multiple jets in a crossflow: Detailed measurements and
numerical simulations," American Society of Mechanical Engineers (Paper)
(1995).
3. V. KASSERA, "Three dimensional CFDanalysis of tube vibrations induced by
crossflow," 4 International Symposium on FluidStructure Interactions,
Aeroelasticity, FlowInduced Vibrations and Voice ASME, Aerospace Division
AD, 532 (1997).
4. SCADP/RELAP5 Development Team, "SCDAP/RELAP5/MOD3.2 Code Manual
Volume I: SCDAP RELAP5 Interface Theory," NUREG/CR6150, 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/CR5535,
INEL95/0174, INEL, Idaho Falls (June 1995).
6. SCADP/RELAP5 Development Team, "SCDAP/RELAP5/MOD3.2 Code Manual
Volume III: User's Guide and Input ManualAppendix A," NUREG/CR6150,
INEL96/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 stateofthe
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
(1994).
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 TM107, Fluent Inc., Lebanon
(1993).
18. T. H. SHIH, W. W. LIOU, A. SHABBIR and J. ZHU, "A new kepsilon eddy
viscosity model for high Reynolds number turbulent flows Model development
and validation," Computer Fluids 24, 3 (1995).
19. P. SPALART, S. ALLMARAS, "A oneequation turbulence model for
aerodynamic flows," Technical Report AIAA920439, 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).
BIOGRAPHICAL SKETCH
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
assistant.
