# CFD Simulation of Thermal-Hydraulic Benchmark V1000CT-2 Using ANSYS CFX

Authors: Thomas Höhne
Journal: Science and Technology of Nuclear Installations
 (2009)

---

## Abstract

Plant measured data from VVER-1000 coolant mixing experiments were used within the OECD/NEA
and AER coupled code benchmarks for light water reactors to test and validate computational fluid
dynamic (CFD) codes. The task is to compare the various calculations with measured data, using
specified boundary conditions and core power distributions. The experiments, which are provided for
CFD validation, include single loop cooling down or heating-up by disturbing the heat transfer in the
steam generator through the steam valves at low reactor power and with all main coolant pumps in
operation. CFD calculations have been performed using a numerical grid model of 4.7 million
tetrahedral elements. The Best Practice Guidelines in using CFD in nuclear reactor safety applications
has been used. Different advanced turbulence models were utilized in the numerical simulation. The
results show a clear sector formation of the affected loop at the downcomer, lower plenum and core
inlet, which corresponds to the measured values. The maximum local values of the relative temperature
rise in the calculation are in the same range of the experiment. Due to this result, it is now possible to
improve the mixing models which are usually used in system codes.


---

## Body

## 1. Introduction

Several mixing phenomena characterize the
various operating conditions of pressurized water reactors (PWRs) and influence
the safety analyses of the plant operating states. Computational fluid dynamics
(CFD) is the best suited tool to study such phenomena in detail. Since there
are large uncertainties in the proper application of turbulence models in
various cases, the validation of CFD codes for reactor applications requires
well-defined experiments.Also recent coupled code benchmarks [1]
have identified coolant mixing in the reactor vessel as an unresolved issue in
the analysis of complex plant transients with reactivity insertion. As a result
phase 2 of the VVER-1000 coolant transient benchmark, [2] was defined aiming at
mixing models testing and single effect analysis of main steam line break (MSLB)
transients with improved vessel thermal-hydraulic models. One purpose of the
V1000CT-2 thermal-hydraulics benchmark was in general to test the capability of
CFD codes to represent vessel thermal hydraulics and to analyze in particular
the coolant mixing in the downcomer and lower plenum of the reactor vessel.The experiment includes single-loop heating up
by disturbing the heat transfer in the steam generator (SG) through the steam
valves, at low reactor power in the range of 5–14% and with all main coolant pumps (MCPs) in operation. It was conducted during the plant commissioning phase
at Kozloduy-6.


## 2. The VVER-1000 Reactor Design

The Russian VVER-1000 reactor type(1000MWel) constructed by
Gidropress/Podolsk is a four-loop pressurized water reactor (PWR) with
hexagonal core geometry and horizontal steam generators. The core contains 163
hexagonal fuel assemblies. The geometry of the reactor vessel is presented in
Figure 1. The primary circuit coolant flows to the
core through the perforated elliptical sieve plate and perforated support
columns serving as flow distributors. The support columns are inserted into
corresponding holes of the core inlet plate and welded together at the top so that
no flow passes around the columns. The primary coolant flows through the slots
into the columns, and then further upward through the support columns into the
fuel assemblies. The location of the inlet and outlet nozzles of the reactor
vessel is nonuniform in azimuthal direction. Measurements taken on the
Kozloduy-6 reactor have shown small discrepancies in these angles with respect
to the design values. These angular differences were taken into account in the
CAD model for grid generation.Figure 1
Geometry of the VVER-1000 reactor
vessel.Coolant mixing
experiments at Kozloduy unit 6 were conducted at the beginning of cycle 1
during rise to power. The purpose of the selected experiment “Experiment 1” was
to determine the mixing coefficients (rate of mass exchange) between cold and hot legs and from
cold legs to the inlet of fuel assemblies. Additionally, the azimuthal rotation
of the loop flows relative to cold leg axes has been determined. The mixing
“Experiment 1” was initiated by disturbing loop no. 1. The experiment was
conducted in three states: a stabilized initial state, a transient state, and a
stabilized final state. Additionally, pressure losses
were measured at different locations of the reactor pressure vessel (RPV)
during nominal operation conditions. These values are used for modelling
resistance coefficients in the CFD calculation.All four main
coolant pumps and four steam generators were in operation. The thermal power of
the reactor was 281 MW, that is, 9.36% of the nominal value. The pressure above
the core was 15.59 MPa, close to the nominal value of 15.7 MPa. The coolant
temperature at the reactor inlet was268.6°C, 19.2 K below the nominal cold leg
temperature.A transient was
initiated by closing the steam isolation valve of SG-1 and isolating SG-1 from
feed water. The coolant temperature in the cold and hot legs of loop no. 1 rose
by 13–13.5C°, and the
mass flow rate decreased by 3.4%. The mass flow rate through the reactor is decreased
by 1%. The reactor power changed by 0.16% calculated from primary circuit
balance. The initially symmetric core power distribution did not change
significantly.The relative temperature
distribution at the core inlet (Figure2) has been calculated for the final
state from the measured core outlet temperatures and measured relative fuel
assembly temperature rise in the initial state. The temperatures at the fuel
assembly without thermocouples are interpolated linearly from measured values.
The temperature rise was assumed constant during the transient due to the
constant normalized power distribution. They were calculated from measured cold
leg and hot leg temperatures in the initial state, weighted by the loop mass
flows. The mean value of 3.2 K fuel assembly temperature rise is used to
estimate core inlet temperatures.Figure 2
Estimated and interpolated relative core inlet temperatures rise  [%].


## 3. CFD Code and Sensitivity Analysis According
to BPG

The CFD code for simulating the mixing
studies wasANSYS CFX release 10 [3]. ANSYS CFX is an
element-based finite-volume method with second-order discretization schemes in space and time.
It uses a coupled algebraic multigrid algorithm to solve the linear systems
arising from discretization. The discretization schemes and the multigrid
solver are scalably parallelized. ANSYS CFX works with unstructured
hybrid grids consisting of tetrahedral, hexahedral, prism, and pyramid
elements.The best practice guidelines (BPGs) by Menter [4] were used
to minimize numerical errors and to compare different advanced turbulence
models. In the current study, the CFD simulations were performed according to
these BPGs. A residual convergence
criterion for RMS mass-momentum equations of 1×10−4 was used to ensure negligible small iteration errors.


### 3.1. Advanced
Turbulence Modelling

The following turbulence models [3] were used
to describe the mixing processes.(i) Shear Stress Transport (SST)k-ω-based model
Thek-ω-based SST model
accounts for the transport of the turbulent shear stress. The BSL model
combines the advantages of the Wilcox and the k-ε model via a blending function. For free shear
flows, the SST model is identical to the k-ε model. One of the
advantages of the k-ω formulation is the near wall treatment for low Reynolds
number computations, where it is more accurate and more robust. The convergence
behavior of the k-ω model is often similar to that of the k-ε
model. Since the zonal k-ω models (BSL and SST) include blending
functions in the near wall region that are a function of wall distance, an
additional equation is solved to compute the wall distance at the start of
simulations.(ii) Large Eddy Simulation (LES) Model
Large eddies of the turbulence are
computed and only the smallest eddies are modelled. The main advantage
of LES over computationally cheaper Reynolds-averaged Navier stokes (RANSs) approaches
is the increased level of detail it can deliver. While RANS methods provide “averaged”
results, LES is able to predict instantaneous flow characteristics and resolve
turbulent flow structures. Small-scale turbulence is assumed to
be nearly isotropic and has a more universal characteristic.
Usually, the computational grid serves as a low-pass filter and only the
subgrid scale turbulent phenomena are modelled.
The subgrid scale model in industrial
applications is the one proposed by Smagorinsky; it is an eddy viscosity
model that is based on the assumption that the effect of the small scales
eddies can be accounted for by adding a contribution to the momentum diffusivity.(iii) Detached Eddy Simulation (DES) Model
DES is an attempt to combine elements of RANS
and LES formulations in order to arrive at a hybrid formulation, where RANS is
used inside attached and mildly separated boundary layers. Additionally, LES is
applied in massively separated regions. The idea behind the DES model is to
switch from the SST-RANS model to an LES model in regions, where the turbulent
length,Lt, predicted by the RANS model is larger than the local
grid spacing. In this case, the length scale used in the computation of the
dissipation rate in the equation for the turbulent kinetic energy is replaced
by the local grid spacing. The numerical formulation is switched between a
second-order upwind scheme and a central difference scheme in the RANS and LES
regions, respectively.DES is at least one order of magnitude more
computer intensive than RANS models.


### 3.2. Discretization
Schemes

In the RANS approach (SST model), a steady-state
calculation was performed using the 1st-order UPWIND advection scheme.The LES calculation requires a central difference
advection scheme and a 1st-order backward
Euler transient scheme with a time step size of 0.0001 second fulfilling the Courant-Friedrich-Levy (CFL) criteria (1):(1)Δtconv=1|Ux|/Δx+|Uy|/Δy+|Uz|/Δz


#### 3.2.1. Courant-Friedrich-Levy Criteria (CFL)

The DES calculation requires a switching
procedure of a central difference advection scheme and a higher order UPWIND
scheme and a 2nd-order backward Euler
transient scheme with a time step size of 0.001 second also fulfilling the CFL
criteria.The DES calculation on 8 processors of a
100-processor RedHat LINUX cluster (dual CPU compute nodes XEON, 3.2 GHz,~1.3 Gflops, each containing 2 GBytes RAM) took 2 weeks for 5 seconds simulation
time.


### 3.3. Grid
Generation

The model
consists of the inlets nozzles, downcomer, lower plenum, and a part of the core
and is constructed by 4.67 million unstructured tetrahedral cell elements (Figure3), and the outlines were modelled according to the real plant data. Grid refinement
was done at the spacer elements (structures for fixing the core barrel against
the RPV-wall), the perforated elliptic core barrel plate, and the core support columns
(Figure 4). The purpose of the bottom plate is to equalize the flow profile by
a large pressure loss. Additional pressure loss coefficients were introduced to
address provided design pressure drops measured for nominal steady-state
conditions.Figure 3
Flow domain.Figure 4
Lower
plenum structures.Two porous
regions were modelled, the elliptical sieve plate with a stream wise resistance
loss coefficient of0.101[m−1] and a
transverse multiplier coefficient of 1000 and the inflow into the support columns
with an isotropic loss coefficient of 0.1[m−1].


### 3.4. Boundary
Conditions

The calculation
domain and the inlet boundary conditions at the RPV nozzles of the four loops
are given in Table1. An outlet condition was imposed above the core inlet
plate at approximately 1/4 of the core height. This part of the core was
modelled as an open volume. Studies in [5] have shown that the mixing is not
affected by this simplification. Wall functions were applied on all solid
structures. The walls are treated as adiabatic. The physical properties of the
fluid are those of water at 270C° and 16 MPa.
(i)
Density(ρ): 784kg·m−3.(ii)
Dynamic
viscosity(μ): 101.2010-6Ns·m−2.(iii)
Thermal
conductivity(λ): 0.61W·m−1K−1.(iv)
Heat capacity(Cp): 4910.82J·kg−1K−1.Table 1
Inlet boundary conditions.
Loopn°Volume flux[m3/s]Velocity[m/s]Temperature[C°]16.0777810.71069282.226.0694410.69599269.936.0777810.71069269.046.1805610.89181269.2


## 3.1. Advanced
Turbulence Modelling

The following turbulence models [3] were used
to describe the mixing processes.(i) Shear Stress Transport (SST)k-ω-based model
Thek-ω-based SST model
accounts for the transport of the turbulent shear stress. The BSL model
combines the advantages of the Wilcox and the k-ε model via a blending function. For free shear
flows, the SST model is identical to the k-ε model. One of the
advantages of the k-ω formulation is the near wall treatment for low Reynolds
number computations, where it is more accurate and more robust. The convergence
behavior of the k-ω model is often similar to that of the k-ε
model. Since the zonal k-ω models (BSL and SST) include blending
functions in the near wall region that are a function of wall distance, an
additional equation is solved to compute the wall distance at the start of
simulations.(ii) Large Eddy Simulation (LES) Model
Large eddies of the turbulence are
computed and only the smallest eddies are modelled. The main advantage
of LES over computationally cheaper Reynolds-averaged Navier stokes (RANSs) approaches
is the increased level of detail it can deliver. While RANS methods provide “averaged”
results, LES is able to predict instantaneous flow characteristics and resolve
turbulent flow structures. Small-scale turbulence is assumed to
be nearly isotropic and has a more universal characteristic.
Usually, the computational grid serves as a low-pass filter and only the
subgrid scale turbulent phenomena are modelled.
The subgrid scale model in industrial
applications is the one proposed by Smagorinsky; it is an eddy viscosity
model that is based on the assumption that the effect of the small scales
eddies can be accounted for by adding a contribution to the momentum diffusivity.(iii) Detached Eddy Simulation (DES) Model
DES is an attempt to combine elements of RANS
and LES formulations in order to arrive at a hybrid formulation, where RANS is
used inside attached and mildly separated boundary layers. Additionally, LES is
applied in massively separated regions. The idea behind the DES model is to
switch from the SST-RANS model to an LES model in regions, where the turbulent
length,Lt, predicted by the RANS model is larger than the local
grid spacing. In this case, the length scale used in the computation of the
dissipation rate in the equation for the turbulent kinetic energy is replaced
by the local grid spacing. The numerical formulation is switched between a
second-order upwind scheme and a central difference scheme in the RANS and LES
regions, respectively.DES is at least one order of magnitude more
computer intensive than RANS models.


## 3.2. Discretization
Schemes

In the RANS approach (SST model), a steady-state
calculation was performed using the 1st-order UPWIND advection scheme.The LES calculation requires a central difference
advection scheme and a 1st-order backward
Euler transient scheme with a time step size of 0.0001 second fulfilling the Courant-Friedrich-Levy (CFL) criteria (1):(1)Δtconv=1|Ux|/Δx+|Uy|/Δy+|Uz|/Δz


### 3.2.1. Courant-Friedrich-Levy Criteria (CFL)

The DES calculation requires a switching
procedure of a central difference advection scheme and a higher order UPWIND
scheme and a 2nd-order backward Euler
transient scheme with a time step size of 0.001 second also fulfilling the CFL
criteria.The DES calculation on 8 processors of a
100-processor RedHat LINUX cluster (dual CPU compute nodes XEON, 3.2 GHz,~1.3 Gflops, each containing 2 GBytes RAM) took 2 weeks for 5 seconds simulation
time.


## 3.2.1. Courant-Friedrich-Levy Criteria (CFL)

The DES calculation requires a switching
procedure of a central difference advection scheme and a higher order UPWIND
scheme and a 2nd-order backward Euler
transient scheme with a time step size of 0.001 second also fulfilling the CFL
criteria.The DES calculation on 8 processors of a
100-processor RedHat LINUX cluster (dual CPU compute nodes XEON, 3.2 GHz,~1.3 Gflops, each containing 2 GBytes RAM) took 2 weeks for 5 seconds simulation
time.


## 3.3. Grid
Generation

The model
consists of the inlets nozzles, downcomer, lower plenum, and a part of the core
and is constructed by 4.67 million unstructured tetrahedral cell elements (Figure3), and the outlines were modelled according to the real plant data. Grid refinement
was done at the spacer elements (structures for fixing the core barrel against
the RPV-wall), the perforated elliptic core barrel plate, and the core support columns
(Figure 4). The purpose of the bottom plate is to equalize the flow profile by
a large pressure loss. Additional pressure loss coefficients were introduced to
address provided design pressure drops measured for nominal steady-state
conditions.Figure 3
Flow domain.Figure 4
Lower
plenum structures.Two porous
regions were modelled, the elliptical sieve plate with a stream wise resistance
loss coefficient of0.101[m−1] and a
transverse multiplier coefficient of 1000 and the inflow into the support columns
with an isotropic loss coefficient of 0.1[m−1].


## 3.4. Boundary
Conditions

The calculation
domain and the inlet boundary conditions at the RPV nozzles of the four loops
are given in Table1. An outlet condition was imposed above the core inlet
plate at approximately 1/4 of the core height. This part of the core was
modelled as an open volume. Studies in [5] have shown that the mixing is not
affected by this simplification. Wall functions were applied on all solid
structures. The walls are treated as adiabatic. The physical properties of the
fluid are those of water at 270C° and 16 MPa.
(i)
Density(ρ): 784kg·m−3.(ii)
Dynamic
viscosity(μ): 101.2010-6Ns·m−2.(iii)
Thermal
conductivity(λ): 0.61W·m−1K−1.(iv)
Heat capacity(Cp): 4910.82J·kg−1K−1.Table 1
Inlet boundary conditions.
Loopn°Volume flux[m3/s]Velocity[m/s]Temperature[C°]16.0777810.71069282.226.0694410.69599269.936.0777810.71069269.046.1805610.89181269.2


## 4. Computational Results

### 4.1. Flow Field in the RPV

Figures5 and 6 are showing the flow field which
establishes in the VVER-1000 reactor during normal operation conditions. The
temperature profile at the core outlet, relevant for the determination of the
reactor power and thus for economical plant operation, is directly influenced
by the flow distribution at the core inlet. Figures 5 and 7 show a stable flow
field in the downcomer. The coolant of loop no. 1 is basically covering the
corresponding sector of the loop in the downcomer. The highest values of the
velocity appear below the inlet nozzles. The spacer elements do only slightly
disturb the flow (Figure 5). Figure 6 shows the flow through the lower plenum
structures. The coolant is entering the lower plenum; is flowing through the
perforated plate in upwards direction, besides the support columns; is entering
the support columns; is flowing through the perforations and through the core
inlet plate into the core region (see also Figure 8).Figure 5
Stream
lines in the downcomer.Figure 6
Stream
lines in the lower plenum structures.Figure 7
Azimuthal velocity
profile in the downcomer at two different horizontal positions.Figure 8
Velocity distribution at the core inlet over the fuel element positions.


### 4.2. Temperature Distribution at the Core Inlet

It is assumed
that the flow and temperature field at the final state is independent of the
initial state as well as of the transient of the experiment. The calculated
temperature field of the SST turbulence model is shown in Figures9 and 10. In
Figure 9, the inlet nozzle plane is shown; while in Figure 10 the wall
temperature of the vessel for the stabilized period is given. It is visible
that the flow turns already in the downcomer slightly in counter-clockwise
direction due to the nonuniform and asymmetric azimuthal distribution of the
cold leg nozzles (Figure 9).Figure 9
Temperature distribution in the horizontal
inlet nozzle plane, SST model.Figure 10
Temperature
distribution at the outer wall of the downcomer, SST model.The comparison
of the relative temperature rise at the core inlet calculated with three
different turbulence models, and the measured relative temperature rise is
shown in Figure11. Red color represents maximum temperature changes; blue
color describes minimum changes.Relative temperature rise[%] at the
core inlet calculated with the three different turbulence models.
(a)(b)(c)(d)In these
Figures, the arrows represent the axes of the cold legs. The flow center maximum
of the flow coming from loop 1 is rotated in counter-clockwise direction by about34° in the experiment. This displacement could be partly reproduced by ANSYS
CFX; a difference of 6° in clockwise direction is still remaining. It is
important to note that this rotation has been calculated when the real
Kozloduy-6 geometry is used (differences in angular positions of inlet nozzles
compared to the design values). The LES and DES CFX
calculations were done with the steady-state flow field, and the transient slug
behavior (temperature rise) was modelled. At the end, the results of the relative
temperature change were
interpolated in time in both models (2–5 seconds of
simulation time). The RANS calculation was done in a steady-state mode,
therefore no time interpolation was necessary. The best agreement with the
Kozloduy Experiment 1 at the core inlet is shown by the DES simulation. The
results of all models in agreement with the experiment show a clear sector
formation of the affected loop at the downcomer, lower plenum, and core inlet.
The maximum local values of the relative temperature rise show a good agreement
at the core inlet. It amounts in the experiment 97.7% and in the DES
calculation 97.3% (Figure 12).Figure 12
Relative temperature rise at the core inlet calculated with the
three different turbulence models in comparison with the experimental values.


## 4.1. Flow Field in the RPV

Figures5 and 6 are showing the flow field which
establishes in the VVER-1000 reactor during normal operation conditions. The
temperature profile at the core outlet, relevant for the determination of the
reactor power and thus for economical plant operation, is directly influenced
by the flow distribution at the core inlet. Figures 5 and 7 show a stable flow
field in the downcomer. The coolant of loop no. 1 is basically covering the
corresponding sector of the loop in the downcomer. The highest values of the
velocity appear below the inlet nozzles. The spacer elements do only slightly
disturb the flow (Figure 5). Figure 6 shows the flow through the lower plenum
structures. The coolant is entering the lower plenum; is flowing through the
perforated plate in upwards direction, besides the support columns; is entering
the support columns; is flowing through the perforations and through the core
inlet plate into the core region (see also Figure 8).Figure 5
Stream
lines in the downcomer.Figure 6
Stream
lines in the lower plenum structures.Figure 7
Azimuthal velocity
profile in the downcomer at two different horizontal positions.Figure 8
Velocity distribution at the core inlet over the fuel element positions.


## 4.2. Temperature Distribution at the Core Inlet

It is assumed
that the flow and temperature field at the final state is independent of the
initial state as well as of the transient of the experiment. The calculated
temperature field of the SST turbulence model is shown in Figures9 and 10. In
Figure 9, the inlet nozzle plane is shown; while in Figure 10 the wall
temperature of the vessel for the stabilized period is given. It is visible
that the flow turns already in the downcomer slightly in counter-clockwise
direction due to the nonuniform and asymmetric azimuthal distribution of the
cold leg nozzles (Figure 9).Figure 9
Temperature distribution in the horizontal
inlet nozzle plane, SST model.Figure 10
Temperature
distribution at the outer wall of the downcomer, SST model.The comparison
of the relative temperature rise at the core inlet calculated with three
different turbulence models, and the measured relative temperature rise is
shown in Figure11. Red color represents maximum temperature changes; blue
color describes minimum changes.Relative temperature rise[%] at the
core inlet calculated with the three different turbulence models.
(a)(b)(c)(d)In these
Figures, the arrows represent the axes of the cold legs. The flow center maximum
of the flow coming from loop 1 is rotated in counter-clockwise direction by about34° in the experiment. This displacement could be partly reproduced by ANSYS
CFX; a difference of 6° in clockwise direction is still remaining. It is
important to note that this rotation has been calculated when the real
Kozloduy-6 geometry is used (differences in angular positions of inlet nozzles
compared to the design values). The LES and DES CFX
calculations were done with the steady-state flow field, and the transient slug
behavior (temperature rise) was modelled. At the end, the results of the relative
temperature change were
interpolated in time in both models (2–5 seconds of
simulation time). The RANS calculation was done in a steady-state mode,
therefore no time interpolation was necessary. The best agreement with the
Kozloduy Experiment 1 at the core inlet is shown by the DES simulation. The
results of all models in agreement with the experiment show a clear sector
formation of the affected loop at the downcomer, lower plenum, and core inlet.
The maximum local values of the relative temperature rise show a good agreement
at the core inlet. It amounts in the experiment 97.7% and in the DES
calculation 97.3% (Figure 12).Figure 12
Relative temperature rise at the core inlet calculated with the
three different turbulence models in comparison with the experimental values.


## 5. Conclusions

CFD calculations
have been performed for the themalhydraulic benchmark V1000CT-2. The numerical
grid model was generated with the grid generator ANSYS ICEM-CFD and contains
4.7 Mio. tetrahedral elements. Different advanced turbulence models were used
in the numerical simulation. The results of all calculations show a clear
sector formation of the affected loop at the downcomer, lower plenum, and core
inlet, which correspond to the measured values. The maximum local values of the
relative temperature rise in the experiment amount 97.7% and in the calculation
97.3%. Due to this result, it is now possible to
improve the mixing models which are usually used in system codes.



---

*Source: 835162-2009-03-30.xml*