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

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 ﬂuid dynamic (CFD) codes. The task is to compare the various calculations with measured data, using speciﬁed 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. Di ﬀ erent advanced turbulence models were utilized in the numerical simulation. The results show a clear sector formation of the a ﬀ ected 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. unrestricted


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 thermalhydraulic models.One purpose of the V1000CT-2 thermalhydraulics 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.

The VVER-1000 Reactor Design
The Russian VVER-1000 reactor type (1000 MW el ) 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.
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 was 268.6A 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.5 • C, 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 (Figure 2) 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.

CFD Code and Sensitivity Analysis According to BPG
The CFD code for simulating the mixing studies was ANSYS CFX release 10 [3].ANSYS CFX is an element-based finitevolume 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.(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, L t , 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.

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

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.

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 (Figure 3), 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.Two porous regions were modelled, the elliptical sieve plate with a stream wise resistance loss coefficient of 0.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 ].

Boundary Conditions.
The calculation domain and the inlet boundary conditions at the RPV nozzles of the four loops are given in Table 1.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    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).

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 Figures 9 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).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 Figure 11.Red color represents maximum temperature changes; blue color describes minimum changes.
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 about 34 • 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).

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.
Figure 6: Stream lines in the lower plenum structures.

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 .The k-ω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.

Figure 7 :Figure 8 :
Figure 7: Azimuthal velocity profile in the downcomer at two different horizontal positions.

Figure 9 :Figure 10 :
Figure 9: Temperature distribution in the horizontal inlet nozzle plane, SST model.

Figure 11 :
Figure 11: Relative temperature rise [%] at the core inlet calculated with the three different turbulence models.
structures.The walls are treated as adiabatic.The physical properties of the fluid are those of water at 270 • C and 16 MPa.

4. 1 .Figure 12 :
Figure 12: Relative temperature rise at the core inlet calculated with the three different turbulence models in comparison with the experimental values.
• C, 19.2 K below the nominal cold leg temperature.