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 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 was , 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.5, 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.

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 was ANSYS 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 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) -based model
The -based SST model accounts for the transport of the turbulent shear stress. The BSL model combines the advantages of the Wilcox and the model via a blending function. For free shear flows, the SST model is identical to the model. One of the advantages of the formulation is the near wall treatment for low Reynolds number computations, where it is more accurate and more robust. The convergence behavior of the model is often similar to that of the model. Since the zonal 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, , 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):

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 (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.

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 of and a transverse multiplier coefficient of 1000 and the inflow into the support columns with an isotropic loss coefficient of .

3.4. 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 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 270 and 16 MPa. (i)Density : .(ii)Dynamic viscosity : .(iii)Thermal conductivity : .(iv)Heat capacity : .

Table 1: Inlet boundary conditions.

4. Computational Results

4.1. Flow Field in the RPV

Figures 5 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 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).

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 Figure 11. Red color represents maximum temperature changes; blue color describes minimum changes.

Figure 11: Relative temperature rise at the core inlet calculated with the three different turbulence models.

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 in the experiment. This displacement could be partly reproduced by ANSYS CFX; a difference of 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.

Acknowledgments

The Author would like to thank the Benchmark Team, especially U. Bieder (CEA) and N. Kolev (IRNE), for their support and fruitful discussions. The work reported about in this paper was supported by the German Federal Ministry of Economics and Trade within Project no. 150 1260 on scientific-technical cooperation between Germany and Russian Federation in the field of nuclear reactor safety.