Detailed Simulation of the Nominal Flow and Temperature Conditions in a Pre-Konvoi PWR Using Coupled CFD and Neutron Kinetics

: The aim of the numerical study was the detection of possible vortices in the upper part of the core of a Pre-Konvoi Pressurized Water Reactor (PWR) which could lead to temperature cycling. In addition, the practical application of this Computational Fluid Dynamic (CFD) simulation exists in the full 3D analysis of the coolant ﬂow behavior in the reactor pressure vessel of a nuclear PWR. It also helps to improve the design of future reactor types. Therefore, a CFD simulation of the ﬂow conditions was carried out based on a complex 3D model. The geometry of the model includes the entire Reactor Pressure Vessel (RPV) plus all relevant internals. The core is modelled using the porous body approach, the di ﬀ erent pressure losses along and transverse to the main ﬂow direction were considered. The spacer-grid levels were taken into account to the extent that in these areas no cross-ﬂow is possible. The calculation was carried out for nominal operating conditions, i.e., for full load operation. Furthermore, a prototypical End of Cycle (EOC) power distribution was assumed. For this, a power distribution was applied as obtained from a stationary full-core calculation with the 3D neutron kinetics code DYN3D. In order to be able to adequately reproduce ﬂow vortexes, the calculation was performed transiently with suitable Detached Eddy Simulations (DES) turbulence models. The calculation showed ﬂuctuating transverse ﬂow in the upper part of the core, starting at the 8th spacer grid but also revealed that no large dominant vortices exists in this region. It seems that the core acts as a rectiﬁer attenuating large-scale vortices. The analyses included several spacer grid levels in the core and showed that in some areas of the core cross-section an upward increasingly directed transversal ﬂow to the outlet nozzle occurs. In other areas of the core cross-section, on the other hand, there is nearly any cross-ﬂow. However, the following limitations of the model apply: In the model all fuel elements are treated identical and cross ﬂows due to di ﬀ erent axial pressure losses for di ﬀ erent FA types cannot be displayed. The complex structure of the FAs (eg. ﬂow vanes in spacer grids) could also inﬂuence the formation of large-scale vortices. Also, the possible inﬂuence of two-phase ﬂows was not considered.


Introduction
Since 2005, events have occurred involving the increased operational oxidation of M5 fuel rod claddings at several German pressurized water reactors (PWRs). The conspicuous corrosion was observed mainly in the area between the two uppermost spacer grids (SGs), i.e., the eighth and the ninth SGs (Figures 1 and 2). In this area, the transition from the active fuel rod area (filled with fuel) to the fuel rod plenum takes place. In some cases, the increased oxide layer was even found in the region of the fuel rod plenum, where no appreciable power is transferred from the fuel rod to the coolant. and large-scale turbulence. Due to different coolant temperatures in neighboring fuel assemblies (FAs), the turbulent flow could lead to a cyclic temperature load on the cladding material and the initial oxide layer, especially in the upper core region. According to this hypothesis, it is assumed that the oxide layer ruptures due to the load from alternating temperatures, creating pathways where oxidizing species reach the metallic cladding rod surface. As a result, the oxide layer loses its protective effect and thicker oxide layers can develop.

From earlier experimental investigations at ROCOM (Rossendorf Coolant Mixing Test Facility)
it is known that strong large-scale vortex structures are present in the upper plenum of the RPV [1]. Numerical simulations of the flow conditions in a German type PWR were part of this work, using a complex 3D model.   From earlier experimental investigations at ROCOM (Rossendorf Coolant Mixing Test Facility) it is known that strong large-scale vortex structures are present in the upper plenum of the RPV [1]. Numerical simulations of the flow conditions in a German type PWR were part of this work, using a complex 3D model.  The oxidation phenomenon is not fully understood. However, it seems that the oxidation is caused by a combination of cladding material characteristics, which will not be discussed here, as well as specific operational influences. There are various hypotheses about the operational influences.
One of the hypotheses assumes that an initial oxide layer first arises more or less over the entire length of the fuel rod. Depending on the height position, different thermal conditions may occur in the form of temperature fluctuations. These might be due to the fact that the upper core area could be affected by increased cross-flow, as it is expected that the flow starts "orientation" towards the reactor pressure vessel (RPV) outlet even below the uppermost SG. This flow can show small-scale and large-scale turbulence. Due to different coolant temperatures in neighboring fuel assemblies (FAs), the turbulent flow could lead to a cyclic temperature load on the cladding material and the initial oxide layer, especially in the upper core region. According to this hypothesis, it is assumed that the oxide layer ruptures due to the load from alternating temperatures, creating pathways where oxidizing species reach the metallic cladding rod surface. As a result, the oxide layer loses its protective effect and thicker oxide layers can develop.
From earlier experimental investigations at ROCOM (Rossendorf Coolant Mixing Test Facility) it is known that strong large-scale vortex structures are present in the upper plenum of the RPV [1]. Numerical simulations of the flow conditions in a German type PWR were part of this work, using a complex 3D model.

Geometry of the CFD Model
The geometry of the CFD model includes the entire RPV with all the relevant details. These include the RPV inlet and outlet nozzles, the hot leg and cold leg parts of the main coolant lines, the downcomer, the lower RPV plenum with internals, the lower core support structure (LCSS), the FAs with the spacer grids and the FA head (Figure 1), the upper core support structure (UCSS) and the installations in the upper RPV plenum.
In Figure 3, the primary circuit of a PWR is shown schematically. The model geometry is outlined in red. Figure 4 shows a top view of the model geometry with the four cold and four hot legs. In Figure 5, the CFD domain of the Pre-Konvoi plant is marked in yellow. Based on the detailed design, the CFD grid model ( Figure 6) was created. The complete grid comprises 19.4 million elements, with 13,137,993 tetrahedral elements, 245,024 wedges, 32,318 pyramids and 6,071,267 hexahedral elements. The grid was created using the grid generation program ICEM CFD. The wall area has been greatly refined. Figure 7 shows the core as a simplified porous medium with the spacer planes and the FA head or FA foot. A section through the inlet nozzle plane (z = 0 m) is given in Figure 8.

Geometry of the CFD Model
The geometry of the CFD model includes the entire RPV with all the relevant details. These include the RPV inlet and outlet nozzles, the hot leg and cold leg parts of the main coolant lines, the downcomer, the lower RPV plenum with internals, the lower core support structure (LCSS), the FAs with the spacer grids and the FA head (Figure 1), the upper core support structure (UCSS) and the installations in the upper RPV plenum.
In Figure 3, the primary circuit of a PWR is shown schematically. The model geometry is outlined in red. Figure 4 shows a top view of the model geometry with the four cold and four hot legs. In Figure 5, the CFD domain of the Pre-Konvoi plant is marked in yellow. Based on the detailed design, the CFD grid model ( Figure 6) was created. The complete grid comprises 19.4 million elements, with 13,137,993 tetrahedral elements, 245,024 wedges, 32,318 pyramids and 6,071,267 hexahedral elements. The grid was created using the grid generation program ICEM CFD. The wall area has been greatly refined. Figure 7 shows the core as a simplified porous medium with the spacer planes and the FA head or FA foot. A section through the inlet nozzle plane (z = 0 m) is given in Figure 8.

Geometry of the CFD Model
The geometry of the CFD model includes the entire RPV with all the relevant details. These include the RPV inlet and outlet nozzles, the hot leg and cold leg parts of the main coolant lines, the downcomer, the lower RPV plenum with internals, the lower core support structure (LCSS), the FAs with the spacer grids and the FA head (Figure 1), the upper core support structure (UCSS) and the installations in the upper RPV plenum.
In Figure 3, the primary circuit of a PWR is shown schematically. The model geometry is outlined in red. Figure 4 shows a top view of the model geometry with the four cold and four hot legs. In Figure 5, the CFD domain of the Pre-Konvoi plant is marked in yellow. Based on the detailed design, the CFD grid model ( Figure 6) was created. The complete grid comprises 19.4 million elements, with 13,137,993 tetrahedral elements, 245,024 wedges, 32,318 pyramids and 6,071,267 hexahedral elements. The grid was created using the grid generation program ICEM CFD. The wall area has been greatly refined. Figure 7 shows the core as a simplified porous medium with the spacer planes and the FA head or FA foot. A section through the inlet nozzle plane (z = 0 m) is given in Figure 8.                   The lower plenum internal (Schemel), which is typical for Pre-Konvoi systems, is shown in Figure 10. This design is intended to equalize the flow in the lower plenum and must be resolved in detail in the modeling in order to correctly capture the flow distribution at the core inlet.     The lower plenum internal (Schemel), which is typical for Pre-Konvoi systems, is shown in Figure 10. This design is intended to equalize the flow in the lower plenum and must be resolved in detail in the modeling in order to correctly capture the flow distribution at the core inlet. The lower plenum internal (Schemel), which is typical for Pre-Konvoi systems, is shown in Figure 10. This design is intended to equalize the flow in the lower plenum and must be resolved in detail in the modeling in order to correctly capture the flow distribution at the core inlet. The lower plenum internal (Schemel), which is typical for Pre-Konvoi systems, is shown in Figure 10. This design is intended to equalize the flow in the lower plenum and must be resolved in detail in the modeling in order to correctly capture the flow distribution at the core inlet.

Technological Boundary Conditions
The calculations were carried out for normal operating conditions, i.e., for full load operation. The mass flow rate through the core and the resulting pressure losses over the core and over the entire RPV, respectively, were set according to the information obtained from the vendor.
Furthermore, a typical end-of-cycle (EOC) power distribution was assumed in the core. For this purpose, the nodal power distribution was implemented into the CFD model as obtained from a stationary full-core calculation with the neutron kinetics code DYN3D for a generic Pre-Konvoi core loading pattern [2].
In summary, the following technological assumptions were made: -Typical EOC power distribution in the core (calculation with DYN3D for a generic Pre-Konvoi core loading pattern, off-line coupling with CFD); -Mass flow rate per loop 4910.5 (kg/s); -Temperature at the RPV inlet: 293.8 • C; -RPV internal pressure: 157 bar (Reactor Inlet).

Physical and Numerical Model Properties
The used CFD code ANSYS CFX [3] is a finite-volume CFD program system of state-of-the-art technology for solving the Navier-Stokes equations. The postprocessor enables the extensive qualitative and quantitative analysis of the calculation results. ANSYS CFX-18 relies on unstructured, hybrid computing grids in which tetrahedral, hexahedral, prism and pyramid elements can be combined in any order. Due to the flexible grid generation an optimum representation of complex geometries is possible.
As quality assurance measure grid studies were performed according to the best practice guidelines [4] for CFD applications in nuclear engineering.
The CFD simulations were based on the following assumptions: -Adiabatic walls with standard logarithmic wall function; -Pressure boundary condition at the hot main coolant pipe legs; -Incompressible coolant; -Transient CFD calculation with initially quiescent coolant; -Second order spatial and temporal discretization scheme; -Turbulence model: stress blended eddy simulation (SBES) In order to reflect the assumed turbulence, the calculations were performed with a hybrid turbulence model. This method combines the advantages of Reynolds Averaged Navier-Stokes equations (RANS) and Large Eddy Simulations (LES), which reduces the modelling part of the RANS and automatically integrates with the grid resolution. As a result, the RANS part near the wall is very high, whereas in the free stream, the vortices are resolved by means of LES.
The SBES turbulence model belongs to the hybrid RANS-LES turbulence model family [5]. SBES allows combining different RANS and LES models with a blending feature to automatically switch between them.
This turbulence model was successfully used to model density-driven flows in ROCOM. The ROCOM mixing test facility provided a unique experimental basis for the validation of CFD tools for the description of turbulent flows [6,7]. Benchmark tests [8] were used to test the turbulence models under different flow conditions to study the influence of reactor geometry, and to analyze the choice of constraints, discretization scheme, time step size, and grid size.

Porous Body Modeling of the Core
The core is represented as a heterogeneous anisotropic porous body. This means that the head losses are different for axial (streamwise) and transverse flow. The SG zones have been taken into account in such a way that in these areas no transverse flow is possible. Moreover, the axial head losses are different for SG zones, FA foot, FA head, and the open FA regions between the SGs. The head losses along the axial flow direction, as usually applied in thermal hydraulic system code calculations, are provided in Table 1. The axial drag coefficients for the CFD calculations are matched to these pressure drops. The following sink terms are inserted in the momentum equation: The permeabilities K S perm and the resistance coefficients K S loss for streamwise (S) flow are defined as follows: The resistance coefficients for transverse (T) flow direction are set according to the streamwise coefficient multiplier option offered by ANSYS CFX [3]: For the SG area, the LCSS and the UCSS, the multiplier was set to 10 in order to widely block the transverse flow in these areas.
The velocities U X , U Y , U Z are superficial velocities. In order to convert them into real flow velocities, they need to be divided by the respective porosity (about 0.4 in the vertical and about 0.25 in the horizontal direction). For example, the calculated streamwise superficial velocity of 2 m/s corresponds to a flow velocity of approx. 5 m/s in the core. Similarly, the superficial transversal velocity of 0.05 m/s is equivalent to a transverse flow velocity of 0.2 m/s.

Power Distribution within the Core
The power distribution in the reactor core was calculated with the help of a stationary core calculation using the neutron kinetics code DYN3D [9]. DYN3D is a nodal 3D reactor dynamics code for stationary and time-dependent core analyses of light water reactors. The DYN3D burn-up module was used to perform a burn-up calculation for a generic KONVOI equilibrium core loading [10] until the end of the cycle. This generic core contains 64 mixed oxide (MOX) fuel assemblies. The stationary core calculation for determining the power distribution was carried out under the following boundary conditions: As it is usual with nodal codes, the radial neutron kinetic resolution is one node per fuel assembly. In the thermal hydraulic and fuel rod module, a mean fuel rod per fuel assembly is modeled. The axial resolution is 16 layers and is identical for the neutron kinetics and thermal hydraulics modules of DYN3D. This includes one reflector node at the core bottom and the top, respectively. Figures 11 and 12 present a horizontal cross-section of the rod linear power rate in the reactor core and the axial distribution of the rod linear power rate of one fuel assembly. As it is usual with nodal codes, the radial neutron kinetic resolution is one node per fuel assembly. In the thermal hydraulic and fuel rod module, a mean fuel rod per fuel assembly is modeled. The axial resolution is 16 layers and is identical for the neutron kinetics and thermal hydraulics modules of DYN3D. This includes one reflector node at the core bottom and the top, respectively. Figures 11 and 12 present a horizontal cross-section of the rod linear power rate in the reactor core and the axial distribution of the rod linear power rate of one fuel assembly.    The power distribution calculated with the above boundary conditions was extracted from DYN3D and mapped to 13 equivalent height layers for use in CFX. As a result, 2509 nodal power densities (values for 193 fuel elements in 13 height layers each) were provided. The volume sources thus extracted from the DYN3D calculation were transferred to CFX and interpolated according to the CFX grid. Finally, the power distribution was re-calibrated to the total thermal output of 3950 MW. The resulting mean core heat-up in both calculations (DYN3D and CFX) was approx. 35 K. Figures 13 and 14 show the distribution of the heat sources in the CFX calculation. The power distribution calculated with the above boundary conditions was extracted from DYN3D and mapped to 13 equivalent height layers for use in CFX. As a result, 2509 nodal power densities (values for 193 fuel elements in 13 height layers each) were provided. The volume sources thus extracted from the DYN3D calculation were transferred to CFX and interpolated according to the CFX grid. Finally, the power distribution was re-calibrated to the total thermal output of 3950 MW. The resulting mean core heat-up in both calculations (DYN3D and CFX) was approx. 35 K. Figures 13  and 14 show the distribution of the heat sources in the CFX calculation.

Results
Figures 15 and 16 present a streamline view of the flow through one sector of the RPV (cold leg sector 1 CL 1). The flow in the downcomer forms sector patterns as can be seen in these figures. The flow is bounded laterally by the coolant flows from the adjacent sectors. Low velocity recirculation areas are located below the cold leg nozzles and high velocity areas between the cold leg nozzles. Figures 17 and 18 show the temperature distribution in the RPV. The coolant is heated up by nuclear fission (in the CFD model by heat source distribution from DYN3D) while flowing through the core (Figure 18) towards the upper RPV plenum. A temperature stratification emerges in the hot legs (hotter medium at the top, see Figure 17).

Results
Figures 15 and 16 present a streamline view of the flow through one sector of the RPV (cold leg sector 1 CL 1). The flow in the downcomer forms sector patterns as can be seen in these figures. The flow is bounded laterally by the coolant flows from the adjacent sectors. Low velocity recirculation areas are located below the cold leg nozzles and high velocity areas between the cold leg nozzles. Figures 17 and 18 show the temperature distribution in the RPV. The coolant is heated up by nuclear fission (in the CFD model by heat source distribution from DYN3D) while flowing through the core (Figure 18) towards the upper RPV plenum. A temperature stratification emerges in the hot legs (hotter medium at the top, see Figure 17). Figure 19 shows a 3D vector plot of the velocity field in the lower part of the downcomer and in the lower RPV plenum.
While flowing through the lower RPV plenum and the lower plenum internals (Schemel), it is deflected by 180 • before it enters the core. Here are recognizable strong vortices. Figure 20 shows the velocity vector views in the upper part of the core between SG levels 8 and 9 (upper part of figure) and in the FA foot area (bottom part of figure), respectively. Different absolute velocities per FA position can be recognized in the FA foot zone, but have almost disappeared in the upper part of the core.            Figure 19 shows a 3D vector plot of the velocity field in the lower part of the downcomer and in the lower RPV plenum.   Figure 19 shows a 3D vector plot of the velocity field in the lower part of the downcomer and in the lower RPV plenum.    In order to see what really happens to the transverse flow components, a closer view of the velocity field is necessary. Figure 21 shows the nomenclature (SIEMENS KWU) of the fuel element positions in plane view. Additionally, the positions of the cold (KS) and hot legs are indicated. This scheme was used in the following figures. In order to see what really happens to the transverse flow components, a closer view of the velocity field is necessary. Figure 21 shows the nomenclature (SIEMENS KWU) of the fuel element positions in plane view. Additionally, the positions of the cold (KS) and hot legs are indicated. This scheme was used in the following figures.    Figure 26 shows the velocity component w. Figure 27 displays the temperature distribution. This is determined by the heat source distribution of the DYN3D calculation.  Figure 26 shows the velocity component w. Figure 27 displays the temperature distribution. This is determined by the heat source distribution of the DYN3D calculation.                 In the peripheral zones, the temperature increase is lower as well as in the control rod positions. Figures 28-31 show the time-dependent transverse velocity components for different planes in three exposed positions (O10, N13 and E2). Starting with Figure 28 in the FA head, these transverse components are the largest. There are also larger fluctuations recognizable. These components become smaller with increasing distance in the direction of the center of the core (downwards) ( Figure 29, the eighth and ninth spacer grid-plane, Figure 30, the seventh and eighth spacer gridplane and Figure 31, the fifth and sixth spacer grid plane) and the fluctuations become smaller. It can also be concluded from the Figures 28-31 that transversal velocities are getting larger towards the core outlet depending on the positions of the hot legs (see Figure 21). Closer to the hot legs, the flow is already directed towards it and therefore the transversal velocity components are getting larger. In the peripheral zones, the temperature increase is lower as well as in the control rod positions. Figures 28-31 show the time-dependent transverse velocity components for different planes in three exposed positions (O10, N13 and E2). Starting with Figure 28 in the FA head, these transverse components are the largest. There are also larger fluctuations recognizable. These components become smaller with increasing distance in the direction of the center of the core (downwards) (Figure 29, the eighth and ninth spacer grid-plane, Figure 30, the seventh and eighth spacer grid-plane and Figure 31, the fifth and sixth spacer grid plane) and the fluctuations become smaller. It can also be concluded from the Figures 28-31 that transversal velocities are getting larger towards the core outlet depending on the positions of the hot legs (see Figure 21). Closer to the hot legs, the flow is already directed towards it and therefore the transversal velocity components are getting larger.

Securing the Results with the Wanninger Model Approach [11]
The results of our approach were compared with an approach from the literature [11]. The physical and numerical model properties according to [11] are shown below. As in our approach, the source and sink terms were used in the momentum equation in the CFD code.
The sink terms in the momentum equation are: The approach in the flow direction is as follows: For K s loss Wanninger uses the McAdams correlation. However, this provides too small pressure losses. This is why it was chosen for this value: as in our approach. The approach in the transverse direction is determined as follows: K T loss = 1.85 Re −0.2 /p FR (δ/(δ 1)) 2 ( Figure 32) with: p FR : Fuel rod pitch, δ: Pitch-to-diameter ratio: p FR /d FR . The Reynolds number is determined according: with d FR the rod diameter, v the gap velocity, and σ the kinematic viscosity. Figure 32 shows the loss coefficient as a function of the Reynolds number according to Peybernes. The results of the physical and numerical model properties according to [11]  For K s loss Wanninger uses the McAdams correlation. However, this provides too small pressure losses. This is why it was chosen for this value: as in our approach. The approach in the transverse direction is determined as follows: K T loss = 1.85 Re −0.2 /pFR(δ/(δ 1)) 2 ( Figure 32) with: pFR: Fuel rod pitch, δ: Pitch-to-diameter ratio: p /d . The Reynolds number is determined according: with dFR the rod diameter, v the gap velocity, and σ the kinematic viscosity.  This corresponds to a factor of 68 over our approach (multiplier approach). Our option for more restricted areas is K T loss × 2: K T loss (spacer level, OKG, UKG) = 544 [m −1 ].
The negation of the permeability term in the momentum sink in [11] may lead to an underestimation of the pressure loss at the core: Δptotal (core): 0.78 bar. This corresponds to a factor of 68 over our approach (multiplier approach). Our option for more restricted areas is K T loss × 2: The negation of the permeability term in the momentum sink in [11] may lead to an underestimation of the pressure loss at the core: ∆p total (core): 0.78 bar.
A comparison of results is shown in Figures 33-36. Figure 33 shows the absolute velocities of the two approaches in the FA foot over the core diameter. As the resistance coefficients in the transversal area are larger according to [11], these influences can already be seen in Figure 34 in the FA foot.  Figure 33 shows the absolute velocities of the two approaches in the FA foot over the core diameter. As the resistance coefficients in the transversal area are larger according to [11], these influences can already be seen in Figure 34 in the FA foot.

Conclusions
The aim of the study was the detection of possible vortices between the eighth and ninth spacer grid plane of a Pre-Konvoi core, which could lead to temperature cycling. The results of the CFD calculation do not show large dominant vortices in the upper core region. The core acts as a rectifier and dampens large-scale vortices. The analyses included several spacer grid levels in the core and showed that in some areas of the core cross-section, an upward increasingly directed transversal flow to the outlet nozzle occurs. In other areas of the core cross-section, on the other hand, there is nearly any cross-flow.
However, the following limitations of the model apply: in the model, all fuel elements are treated as identical and cross flows due to different axial pressure losses for different FA types cannot be displayed. The complex structure of the FAs (e.g., flow vanes in spacer grids) could also influence the formation of large-scale vortices. This effect could only be resolved with a very high number of

Conclusions
The aim of the study was the detection of possible vortices between the eighth and ninth spacer grid plane of a Pre-Konvoi core, which could lead to temperature cycling. The results of the CFD calculation do not show large dominant vortices in the upper core region. The core acts as a rectifier and dampens large-scale vortices. The analyses included several spacer grid levels in the core and showed that in some areas of the core cross-section, an upward increasingly directed transversal flow to the outlet nozzle occurs. In other areas of the core cross-section, on the other hand, there is nearly any cross-flow.
However, the following limitations of the model apply: in the model, all fuel elements are treated as identical and cross flows due to different axial pressure losses for different FA types cannot be displayed. The complex structure of the FAs (e.g., flow vanes in spacer grids) could also influence the formation of large-scale vortices. This effect could only be resolved with a very high number of grid elements (several billions). This is currently not possible from a computational point of view. In addition, the possible influence of two-phase flows was not considered.