INTRODUCTION
Pressurized Water Reactors (PWRs) are the most common type of power producing nuclear reactor that are widely used all over the world to generate electric power. In the PWR, water under pressure as coolant is used in closed system of reactor vessel. A lot of assembles of the fuel rods is housed in a larger vessel, the shell. The heat created by the nuclear fission reaction that takes place in the fuel rods is then extracted by bulk water on the shell side. Generally the limitation of the heat transfer is found on the fuel rod side of the process and especially in the region close to the fuel rod wall (Espinoza, 2006).
In PWR reactors facilitating faster heat transfer in these extremely exothermic processes is usually achieved by designing the reactor as a vessel with narrow fuel rod assemblies. For proper design and operation, an accurate knowledge of the heat transfer properties is required because of high sensitivity of reactor behaviour to some operating parameters, such as coolant mixing temperature (Rohde, 2007). The fact that industrial power reactor operation has been driven as close as possible to runaway conditions for maximum capacity we need a complete understanding of the heat transfer in these reactors. In particular, for reasons of more safety importance in nuclear reactors, specification of the fluid distribution helps us to prevent and predict the possible break loss of coolant accident (Cartland et al., 2007).
Accurate modeling of these PWRs is complicated, especially in high fuel rod to shell diameter ratios and for the large number of fuel rods. The simplest method, which is probably most commonly used, is modelling heat transfer by using some empirical correlations. Due to the empirical nature, this approach might become highly inaccurate. This is because those correlations cannot account for the complex nature of the fluid dynamics and the geometry for specific situation. Consequently, while it is preferable for initial estimate in design, this empirical approach might not be sufficient for acquiring an accurate knowledge of the heat transfer. The most expensive and timeconsuming method is obviously experimental study. However, there are major difficulties inherent in this approach including obtaining temperature profiles (especially for the inside fuel rods).
Fortunately, with new methods such as Computational Fluid Dynamics, (CFD) it is possible to get a detailed view of the fluid flow and heat transfer phenomena in the vessel side of the reactors which recently is becoming of increasing interest. This might be because the threedimensional effects in the flows distribution cannot be predicted well by one dimensional system codes (Scheuerera, 2005).
Although CFD is a relatively young tool in the field of nuclear engineering; it is applied most commonly in mixing technology. The recent developments of the incorporation of chemical and nuclear reaction in the major CFD packages opened up a large area of application of CFD in reactor engineering (Bieder and Graffard, 2007).
In general, CFD involves numerically solving conservation equations for mass,
momentum and energy in flow geometry of interest, together with additional sets
of equations reflecting the problem at hand. This study at first presents the
development of a mathematical model for the shell side of the reactor (where
the coolant mixing exist), which is solved numerically by ANSYS CFX5.7.1, specific
to typical VVER1000 (V320) reactors. Next, the approach to turbulence modelling
and the numerical method implemented in CFX5.7.1 will be described. Then, specification
of the numerical simulation in CFX will be given and the results will then be
presented and discussed.
For this model, we use an Intel core dual processor of 4.8 GHz and 2 GB RAM. Tests with CFX 5.7.1 have shown that this is sufficient for the treatment of single phase flows with standard kε turbulence model using a grid of 352x10^{3} unstructured cells. For grid generation the tool ICEM CFD is used that allows the treatment of 1.7x10^{6} cells. Nevertheless, the response time for a single grid manipulation action for such a case often reaches 15 h. CFX5.7.1 allows the treatment of computational domains whose boundary nodes do not need to match exactly. For testing and handling such an approach is advisable and in the present case even necessary because the complete grid cannot be handled as one part due to limitations of the grid generation processing. Some other computational limits are given by the preprocessing of CFX5, where the computational domains and the flow physics are defined (Bottcher, 2006).
MATERIALS AND METHODS
The starting point for modelling the fluid flow in the vessel side of the reactor is the set of the fundamental equations that can be found in many well known text books (Fox and McDonald, 1978):
Continuity equation:
Momentum equation:
Energy equation:
Where:
and
The second term in E is the kinetic energy per unit mass of a material particle.
Inspection of these equations reveals the background for a few of the common
reactor model assumptions. First, these systems are low velocity flows and the
fluid mass density, ρ can be treated as (incompressible flow) uniform in
the flow field so that the terms containing the time or spatial derivative of
ρ can be neglected. The mass density is not dependent on the pressure changes
due to the flow and the viscous dissipation and pressure terms in the energy
equation can be neglected. Second, the density and thermodynamic coefficients
are not generally constants and may be functions of temperature. The governing
equations for low Mach number flow derived based on the dimensional analysis
can then be expressed as:
Modelling turbulent flows range from Direct Numerical Simulation, (DNS) to
the Reynolds Averaged NavierStokes, (RANS) approach. When RANS approach is
applied to the standard equations, the result is:
Where:
and
The nonlinear terms involved the turbulent fluctuations are called the Reynolds
stress and the Reynolds flux These terms have to be modelled to enable solution
of the NavierStokes equations.
There are many turbulence models for this purpose including zero equation model, kε model, RNG kε model and differential Reynolds Stress model (Muralidhar and Sundararajan, 1995). Only standard kε model, which is of eddy viscosity model type, was used for the CFD simulation. As a result of turbulence, in summary, the mathematical model to be solved for the shell side of the reactor consists of the following equations:
Continuity equation:
Momentum conservation:
Where:
and
Energy conservation:
Turbulence kinetic energy k conservation:
Where:
Turbulence kinetic energy dissipation ε conservation:
Where:
and
Table 1: 
Physical properties of water 

In our model so far there are nine unknown including fluid mass density ρ,
velocity v (including three coordinates), total pressure P, total energy E,
temperature T, turbulence kinetic energy k and turbulence dissipation rate ε.
However, there are only seven Eq. 16, 17,
20, 21, 23 and therefore
to completely specify the model, two more equations are required.
The first Eq. 16, involves the fluid density ρ. For water
coolant, the density was specified as a constant and hence independent of time,
pressure and temperature:
The second Eq. 17, including three equations usually called
constitutive equation relates the enthalpy change to the temperature and pressure.
For constant ρ and c_{p}, it follows:
For water coolant the thermal and physical properties (at pressure 155 bar)
are in Table 1.
Analytical solutions to the NavierStokes equations are impossible to obtain for any systems but the simplest flows under ideal conditions. For real flows, a numerical approach must be adopted whereby a discretization method involves replacing the NavierStokes equations by their algebraic approximations, which can then be solved using a numerical method. The CFD approach uses NavierStokes equations and energy balances over control volumes, small volumes within the geometry at a defined location representing the reactor internals. The size and number of control volumes (mesh density) is user determined and will influence the accuracy of the solutions to a degree. After boundary conditions have been introduced in the system the flow and energy balances are solved numerically. An iteration process decreases the error in the solution until a satisfactory result has been reached. By using CFD in the simulation of coolant of the nuclear reactors a detailed description of the flow behavior within the barrel can be established, which can then be used in more accurate modelling.
The CFX5.7.1 software is based on a Finite Volume, (FV) approach, where the
solution domain i.e., the fluid domain is subdivided into a finite number of
small Control Volumes, (CVs) by meshing. All of the solution variables and fluid
properties are stored at the computational nodes which are assigned at the centre
of the CVs or arranged so that CV faces lie midway between nodes. To complete
the approximation, it is now necessary to estimate the information for each
node in terms of known variables. Several schemes for interpolation practices
have been used including Upwind Interpolation, (UDS) which CFX5.7.1 implements
a modified version of it, where an additional term named Numerical Advection
Correction, (NAC) is included in the interpolation. This makes the approximations
secondorder accurate but at the same time less robust (ANSYS CFX, 2005).
Geometry and technical data of VVER1000: The VVER1000 is a four loop
pressurized water reactor with hexagonal fuel assembly design and horizontal
steam generators. The ANSYS ICEM used to generate the geometrical details; most
of these are modelled accurately, like: inlet nozzles, outlet nozzles, downcomer,
perforated elliptical sieve plate, upper core support plate, spacer ring. The
general characteristic of the reactor is given in Table 2.
Table 2: 
The general data of VVER1000 


Fig. 1: 
Threedimensional image of the reactor
pressure vessel from model in ANSYS ICEMCFD code 
In these nuclear reactors the coolant enters the vessel by the inlets, flows
downwards through the downcomer and enters the lower plenum by passing a perforated
elliptical bottom plate. Then the flow crossing the core bottom plate and enters
the core. The flow is heated up by the core exits from outlets. In this study
we assume that the PWR consists of vessel and 64 fuel rod assemblies. The basic
geometry of considered reactor is given in Fig. 1.
In this study, we only model the shell side of the reactor. From plant data the temperature profile along the fuel rod has been obtained and used in this model.
Calculations
The CFX5.7.1 code: As noted, CFX5.7.1 is a CFDcode using an elementbased finitevolume numerical method with secondorder discretization schemes in space and time. It works with unstructured hybrid grids consisting of tetrahedral, hexahedral, prism and pyramid elements. The other CFX5 options are: (1) Solution of the NavierStokesEquations for steady and transient flows for compressible and incompressible fluids, (2) Modelling of heat transfers and (3) Use of different coordinate systems.
Input deck
General assumption: The following assumptions for the modelling of the coolant flow in pressurized water reactor is made: (1) incompressible fluid, (2) use of the Standard kε turbulence model and (3) pressure boundary condition at the outlet.
Geometrical simplifications, local details: The geometric details of the construction internals have a strong influence on the flow field and on the mixing. Therefore, an exact representation of the inlet region, the downcomer below the inlet region, the eight spacer elements in the downcomer and the lower plenum structures are necessary (Kliem, 2006).
Grid model: In order to receive an optimal net griding for the later flow simulation one must consider the following items: Checking grid number in special regions to minimize numerical diffusion, refinement of the griding in fields with strong changes of the dependent variables, adaptation of the griding to estimated flow lines, generation of nets as orthogonal as possible. In this study, the mesh contained 1.65x10^{6} tetrahedral elements and 1.67x10^{6} nodes.
Boundary conditions: At VVER1000 (V320) nuclear reactor type the inlet boundary conditions (mass flow rate and temperature) were set at the inlet nozzles. No specific velocity profile was given. The wall was modelled using adiabatic conditions (Hohne and Kliem, 2006).
There were three boundary conditions imposed on the model including inlet, outlet and at the fuel rod wall.
Initial conditions for t>= 0:
At inlets:
At outlets:
At fuel rods:
By specifying the fuel rod`s heat transfer coefficient obtained from experiment,
the fuel rod`s temperature can be calculated in CFX5.7.1 using
Where:
q_{FR} 
= 
Heat transfer obtained from solving the heat balance 
T 
= 
Fluid temperature near the tube rods 
T_{FR} 
= 
Fuel rod`s temperature 
RESULTS AND DISCUSSION
Prior to the transient calculations, the steadystate was simulated. The parallel transient calculation with 10 iteration per time step took 18 h of computation time using two processor (dual CPU computer nodes, containing 2 GB RAM). The convergence criteria were set to 1x10^{4} for RMS residuals (mass, momentum and temperature). In the calculated cases the time step of 1s was taken into account.
Steady state simulation
Velocity distribution: A snapshot of the velocity distribution from
elliptic RPV bottom is shown in Fig. 2. It is clearly to be
seen that the flow from four loops covers the space between elliptic bottom
vessel and barrel. Figure 2 shows the streamlines in the region
of the barrel holes, where the flow is passing through in order to reduce the
effect of sector formation on the reactor core.
Table 3: 
Comparison of the temperature at the hotlegs in CFX and experiment 


Fig. 2: 
Distribution of velocity streamlines
in the elliptic RPV bottom at the steady state 
Outlet temperature distribution: A comparison of the predicted temperature
at the Hotlegs (outlet nozzles) with the plant data is shown in Table
3.
The temperature differences in the outlet nozzles are only in the range
up to 2.15K which is not very significant. On the other hand, the temperature
rise at the outlet nozzles is overpredicted by CFX5. This discrepancy could
be due to several reasons including simplification of the geometry, the grid
used for the numerical simulations and/or inaccuracy in the computational model
due to fluid leakage flows not being taken into account. In particular, the
geometric details of the construction internals have a strong influence on the
flow field and on the mixing. In here, the flow field was computed on a threedimensional
structured grid; however, the lower plenum structures and also spacers were
not included in the model and therefore, their effects on the fluid flow were
not observed.
Flow field at the downcomer: We calculate the velocity distribution
at the downcomer by CFX and this result can be seen at the Fig.
3. It can be seen that the flow fields in the downcomer is not very homogenous
and also no recirculation vortices are found. However, a minimum velocity exists
on azimuthal positions approximately below the inlet nozzles.
In Fig.
4, the velocity at azimuthal position at the end of the downcomer is shown.
As can be seen minimum velocities exist at the positions below the inlet nozzles
at the end of the downcomer. There are positions between the inlet nozzles at
the end of the downcomer that the fluid flow has maximum velocity.

Fig. 3: 
Flow field in the downcomer at nominal conditions (steady
state) 

Fig. 4: 
Velocity distribution at the end of downcomer of VVER1000 
Transient simulation: Examining the streamlines presented in Fig.
58 shows streamlines of water flowing in the downcomer
and lower plenum of the PWR. There are four plots in Fig. 4
describing the flow state at 1, 10, 40 and 80 sec. At 1 s, the flow is evenly
distributed around the downcomer. However, while the pump is operational and
the flow rate is at 100%, the streamlines flow around the circumference of the
reactor to recombine opposite the inlet position at a similar height before
moving down through the diffuser and into the downcomer. Note that streamlines
that move directly into the downcomer after entering through the inlet loop
also move around the circumference of the reactor and that there is virtually
no flow down the downcomer in the region below inlet. Figure 58
shows the streamlines in the region of the lower plenum, where the flow is passing
through and around the perforated drum in order to reduce the effect of sector
formation on the reactor core.

Fig. 5: 
Snapshots of the velocity streamlines
in the downcomer by ANSYSCFX code, 1 sec after startup 

Fig. 6: 
Snapshots of the velocity streamlines in the downcomer by
ANSYSCFX code, 10 sec after startup 

Fig. 7: 
Snapshots of the velocity streamlines
in the downcomer by ANSYSCFX code, 40 sec after startup 

Fig. 8: 
Snapshots of the velocity streamlines
in the downcomer by ANSYSCFX code, 80 sec after startup 
CONCLUSION
In this study, a detailed CFD model for a whole reactor pressure vessel of a PWRreactor of VVER1000 type for the simulation of a coolant mixing is presented. The huge computer memory requirements of such a detailed model forced us to find a compromise between the degree of spatial resolution of some design details of the reactor components and our computational limitation. Therefore some elements in the reactor such as fuel rod assemblies are modeled in a simplified way. Nevertheless the final complete modular RPVmodel consists of approximately 2 million cells. The model has been validated by the outlet temperature obtained from experiment. The mathematical background of fluid dynamics and also the CFX simulation result were discussed.
Future work is directed to the development of a model for VVER1000 type reactors to improve the geometry and to study in more details the transient mixing behavior.
ACKNOWLEDGMENTS
The authors like to thank Guilan`s Fanavary Park and also research deputy of the faculty of science of the University of Guilan for the support with computer facilities.
NOMENCLATURE
C_{p} 
= 
Specific heat (J kg^{1} K^{1}) 
C_{ε1} 
= 
Model constant 
C_{ε2} 
= 
Model constant 
C_{μ}=0.09 
= 
kε turbulent model constant 
D/Dt=d/dt + v.d/dx 
= 
Total derivative 
E 
= 
Total energy per unit mass (J kg^{1}) 
e 
= 
Internal energy per unit mass of a material particle (m^{2} sec^{2}) 
g 
= 
Gravitational acceleration (m sec^{2}) 
I 
= 
Unit tensor 
k 
= 
Turbulent kinetic energy (m^{2} sec^{2}) 
P 
= 
Shear production due to turbulence (incompressible flows) (kg m^{1}
sec^{3}) 
P 
= 
Static pressure (kg m^{1} sec^{2}) 
p 
= 
Modified pressure (kg m^{1} sec^{2}) 
q 
= 
Heat added to each material particle is at a rate per unit of Mass 
q_{ε} 
= 
Sources for ε 
q_{k} 
= 
Sources for k 
S_{E} 
= 
Energy sources or sinks (kg m^{1} sec^{3}) 
S_{M} 
= 
Momentum sources or sinks (kg m^{2} sec^{2}) 
T 
= 
Temperature (K) 
t 
= 
Time (sec) 
v 
= 
Velocity (m sec^{1}) 
z 
= 
Axial (m) 
Greek symbols