CFD-Calculation of Fluid Flow in VVER-1000 Reactors

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:

(1)

Momentum equation:

(2)

Energy equation:

(3)

Where:

(4)

(5)

and

(6)

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:

(7)

(8)

(9)

Modelling turbulent flows range from Direct Numerical Simulation, (DNS) to the Reynolds Averaged Navier-Stokes, (RANS) approach. When RANS approach is applied to the standard equations, the result is:

(10)

(11)

(12)

Where:

(13)

(14)

and

(15)

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 Navier-Stokes 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:

(16)

Momentum conservation:

(17)

Where:

(18)

and

(19)

Energy conservation:

(20)

Turbulence kinetic energy k conservation:

(21)

Where:

(22)

Turbulence kinetic energy dissipation ε conservation:

(23)

Where:

(24)

and

Table 1:	Physical properties of water

(25)

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:

(26)

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:

(27)

For water coolant the thermal and physical properties (at pressure 155 bar) are in Table 1.

Analytical solutions to the Navier-Stokes 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 Navier-Stokes equations by their algebraic approximations, which can then be solved using a numerical method. The CFD approach uses Navier-Stokes 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 CFX-5.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 CFX-5.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 second-order accurate but at the same time less robust (ANSYS CFX, 2005).

Geometry and technical data of VVER-1000: The VVER-1000 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 VVER-1000

Fig. 1:	Three-dimensional image of the reactor pressure vessel from model in ANSYS ICEM-CFD 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 CFX-5.7.1 code: As noted, CFX-5.7.1 is a CFD-code using an element-based finite-volume numerical method with second-order discretization schemes in space and time. It works with unstructured hybrid grids consisting of tetrahedral, hexahedral, prism and pyramid elements. The other CFX-5 options are: (1) Solution of the Navier-Stokes-Equations 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⁶ tetrahedral elements and 1.67x10⁶ nodes.

Boundary conditions: At VVER-1000 (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:

(28)

(29)

At outlets:

(30)

At fuel rods:

(31)

By specifying the fuel rod`s heat transfer coefficient obtained from experiment, the fuel rod`s temperature can be calculated in CFX-5.7.1 using

(32)

RESULTS AND DISCUSSION

Prior to the transient calculations, the steady-state 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 hot-legs 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 Hot-legs (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 over-predicted by CFX-5. 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 three-dimensional 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 VVER-1000

Transient simulation: Examining the streamlines presented in Fig. 5-8 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 5-8 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 ANSYS-CFX code, 1 sec after start-up

Fig. 6:	Snapshots of the velocity streamlines in the downcomer by ANSYS-CFX code, 10 sec after start-up

Fig. 7:	Snapshots of the velocity streamlines in the downcomer by ANSYS-CFX code, 40 sec after start-up

Fig. 8:	Snapshots of the velocity streamlines in the downcomer by ANSYS-CFX code, 80 sec after start-up

CONCLUSION

In this study, a detailed CFD model for a whole reactor pressure vessel of a PWR-reactor of VVER-1000 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 RPV-model 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 VVER-1000 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

ρ	=	Fluid mass density (kg m-3)
ε	=	Rate of dissipation (m2 sec-3)
μ	=	Dynamic molecular viscosity (kg m-1 sec-1)
	=	Turbulent viscosity (kg m-1 sec-1)
μeff = μ +μt	=	Effective viscosity (kg m-1 sec-1)
Prt	=	Turbulent Prandtl number(Dimensionless)
Γ	=	Diffusivity (kg m-1 sec-1)
	=	Turbulent diffusivity (kg m-1 sec-1)
Γeff = Γ +Γt	=	Effective diffusivity (kg m-1 sec-1)
	=	Effective diffusivity for k (kg m-1 sec-1)
	=	Effective diffusivity for ε (kg m-1 sec-1)
σk = 1	=	Model constant
σε = 1.3	=	Model constant
λ	=	Thermal conductivity (W m-1 K-1)
λeff	=	Effective thermal conductivity (W m-1 K-1)
	=	Effective thermal conductivity (kg m sec-3 K-1)
	=	Bulk viscosity (kg m-1 sec-1)