# A New Coupled CFD/Neutron Kinetics System for High Fidelity Simulations of LWR Core Phenomena: Proof of Concept

Authors: Jorge Pérez Mañes; Victor Hugo Sánchez Espinoza; Sergio Chiva; Robert Stieglitz
Journal: Science and Technology of Nuclear Installations
 (2014)

---

## Abstract

The Institute for Neutron Physics and Reactor Technology (INR) at the Karlsruhe Institute of Technology (KIT) is investigating the application of the meso- and microscale analysis for the prediction of local safety parameters for light water reactors (LWR). By applying codes like CFD (computational fluid dynamics) and SP3 (simplified transport) reactor dynamics it is possible to describe the underlying phenomena in a more accurate manner than by the nodal/coarse 1D thermal hydraulic coupled codes. By coupling the transport (SP3) based neutron kinetics (NK) code DYN3D with NEPTUNE-CFD, within a parallel MPI-environment, the NHESDYN platform is created. The newly developed system will allow high fidelity simulations of LWR fuel assemblies and cores. In NHESDYN, a heat conduction solver, SYRTHES, is coupled to NEPTUNE-CFD. The driver module of NHESDYN controls the sequence of execution of the solvers as well as the communication between the solvers based on MPI. In this paper, the main features of NHESDYN are discussed and the proof of the concept is done by solving a single pin problem. The prediction capability of NHESDYN is demonstrated by a code-to-code comparison with the DYNSUB code. Finally, the future developments and validation efforts are highlighted.


---

## Body

## 1. Introduction

The safety evaluation methodologies of nuclear power plants (NPP) are in continuous development and validation and must reflect the state of the art in science and technology. Recent advances in computer and engineering sciences have led to a gradual improvement of the prediction capability of numerical tools regarding key safety-relevant phenomena of nuclear reactors. Currently, conservative and best estimate (BE) safety analysis methodologies are being applied with increasing importance of the BE computer codes within the licensing process of NPP.Due to the increasing complexity of fuel assemblies (FA) in design and material compositions, a general trend to move towards multiscale and multiphysics simulations aiming to obtain high fidelity simulations is pursued [1–6].This new trend permits a better assessment of the safety margins and contributes to reduction of the conservatism in the applied methodologies for core analysis. The availability of high performance computers allows an increase of the spatial resolution of the thermal hydraulics (TH) and neutron physical phenomena at very detailed scales paving the way for the possibility to directly predict local safety parameters [4, 5, 7].In this paper, the novel system NHESDYN (NEPTUNE-SYRTHES-DYN3-SP3) is presented. It consists of the coupling of the neutron kinetic code DYN3D-SP3 [8], the CFD solver NEPTUNE-CFD [9], and the heat conduction solver SYRTHES [10]. The main goal of this paper is to present the proof of the concept for the implemented coupling approach meaning that the MPI implementation of the communication between the involved solvers, the update of the cross-sections based on the thermal hydraulic local parameters, is working appropriately and that the new code is predicting physical sound results.The far goal of this development is to have a coupled system able to describe both the thermal hydraulic and neutronic phenomena in LWR cores with fewer approximations than the legacy coupled codes such as RELAP5/PARCS [11] or ATHLET/DYN3D [11]. Since the detailed description of the thermal hydraulics phenomena at macro- and mesolevel by the CFD codes can be very time consuming, though all CFD codes are parallelized, NHESDYN is developed within a MPI-environment. Similar efforts are being followed in different places combining different code versions and types, for example, ANSYS/PARCS [11], MCNP/STAR-CD [12], or DYN3D/ANSYS CFX [13].In this paper, the main features of NHESDYN to solve one single pin problem are presented since the application of NHESDYN to solve pin clusters, fuel assemblies, or even a reactor core is very challenging and too time consuming at the moment. Hence the authors decided to restrict themselves to a simple problem where the potentials of the coupling methodologies can be clearly demonstrated. After the description of the coupling approach and implementation, the validation of NHESDYN is described and discussed in Section3. In Section 4, a summary and the conclusions are given. Finally, the areas for future work are highlighted in Section 5.


## 2.  Coupling Methodology

By coupling NEPTUNE-CFD with the heat conduction solver SYRTHES and the NK code DNY3D-SP3, a new code named NHESDYN was created. In this coupling approach, the codes are externally coupled using the Message Passing Interface (MPI) standards for the communication between them. In Figure1, a simplified communication scheme to understand the main information exchange during a time step is shown. This exchange has several stages.(i)
The SYRTHES code receives the information of the volumetric heat flux from DYN3D-SP3.(ii)
The fuel temperature gradients plus pin gap and cladding heat conduction are solved.(iii)
The heat flux is transferred from the cladding to the fluid domain, solving the conjugated heat transfer, between TH codes.(iv)
The NEPTUNE-CFD code calculates the temperatures, pressure, and velocities of the fluid.(v)
DYN3D-SP3 updates the cross-sections by taking the information of the Doppler temperature (calculated by SYRTHES) and moderator temperature and density (calculated by NEPTUNE-CFD).Figure 1
Simplified communication scheme for the coupled solution NEPTUNE_CFD/SYRTHES/DYN3D_SP3.Finally, the power calculated by the NK code is passed over to the SYRTHES code to restart again the process. This methodology is repeated at each time step. The codes involved have to wait for the information generated by its predecessor to start the calculation and each code is normally executed one time. This kind of time hierarchy in the coupled solutions is called explicit. The solution can be considered a combination of external and internal coupling. The heat conduction solver SYRTHES is coupled in NEPTUNE-CFD by means of a semi-implicit approach; therefore the thermal hydraulic and heat conduction solver are considered as one. The combination of NEPTUNE-CFD/SYRTHES performs an external coupling with the NK tool.Because the effort to validate codes used for safety analysis is very large, a preservation of the established codes is of strong importance. Therefore, coupling between existing codes representing different physics fields such as TH and NK is very attractive compared to the option of developing a complete new code. This approach is followed by the route of a combined coupling and it works quite well for the TH/NK systems where well validated codes exist. The TH tool is always one step ahead of the NK calculation.The time discretization has to be small enough to avoid numerical stability problems. Exploring studies have shown that the explicit method of updating the data is free of oscillations only if the time step is small enough. For large time steps some variables develop unphysical oscillations of numerical nature that are not related to the physics of the problem under consideration. The explicit method is conditionally stable. Normally, the time step is selected by the TH code attending to the flow characteristics like the Courant number. In case of the simulation of a transient scenario with externally coupled codes, both codes share the same time step. The time scale of the coupled solutions is always conditioned by the most restrictive time step. TH tools use a smaller time step compared to the NK code. Therefore, NK code has always a sufficient time step in the coupled solution.


### 2.1. General Coupling Approach

A general overview of the data flow between the codes and how it was implemented is described here. First of all, the first step for the communication among the codes is the creation of the MPI groups. A general scheme concerning the group distribution of the different codes is shown in Figure2. This scheme defines the direct dependencies between the two TH codes involved and the NK code with the main program. The amount of information to deal with is high and to work with three different outputs is tedious. To unify the information provided by the codes an extra program is created; this is the main program (NHESDYN). It is programmed in FORTRAN and it was especially developed to control the time synchronization for the coupled solution and to monitor the progress of the coupled simulations.Figure 2
General distribution of the three MPI groups and the codes included within them.NHESDYN generates DYN3D-SP3 as a subprocess (spawn operation) and creates a point to point communication with NEPTUNE-CFD. At the same time the NK code creates a point to point communication with SYRTHES. In the MPI, the point to point communication is a client and a server.In Figure2, the TH tools are considered like a single code. These codes were already coupled by the code developers (EDF and CEA) and the communication between them is done directly with no MPI interaction, for example, to solve the conjugated heat transfer. Additional information is shown in Figure 2 concerning the kind of information transferred (densities, temperatures, etc.) which is classified by tags.In Figure3, a simplified flow chart to understand the function of each module in the coupled solution is given. The NK/TH tools follow a specific order in the information flow while the main code (NHESDYN) controls the times and performs the convergence loops during some stages of the simulation. Considering that the NK code has solved the previous time step, the information (volumetric heat flux) is sent to the SYRTHES code. Then, the TH step is calculated and the conjugated heat transfer through the solid-liquid interface is solved.Figure 3
Simplified flow chart for the coupled NHESDYN system.The information generated takes two different ways. The information from SYRTHES (Doppler temperature) can be transferred directly to DYN3D-SP3 because both are included in the same MPI group. The second way is followed by the information from NEPTUNE-CFD (moderator temperatures and densities) which is sent to the main code. Here, if it applies some operations are performed (e.g., convergence loops for a steady state case). Afterwards, the information is resent to the NK code. By this last operation the NK has enough information to refresh the precalculated cross-sections and solve the transport equation. When the calculation is finished the heat flux is sent again to the TH code which is waiting for the next time step calculation.


### 2.2. Execution of the Coupled Code NHESDYN

In Figure4, the different steps performed during the coupled solution are illustrated. The coupled simulation is started after a converged thermal hydraulic stand-alone NEPTUNE-CFD/SYRTHES run is performed. Since no steady option for NEPTUNE-CFD exists, a null transient simulation is performed for constant boundary conditions. The thermal hydraulic time step is used later on as a common time step for the other solvers during the iterative coupling approach. During the first iteration, NEPTUNE-CFD and SYRTHES provide the initial conditions for the neutron kinetic solver (DYN3D). With the TH information DYN3D calculates the power distribution. This loop runs till the predefined convergence criteria are achieved.Figure 4
NHESDYN main steps of execution.DYN3D-SP3 first performs a steady state run and then the transient calculation is initiated. In this case, both codes run at the same time and the thermal hydraulic boundary conditions can be modified and these modifications can be taken into account by the transient neutronic solver.


### 2.3. Description of the Time Advancement

The coupling scheme depends on the type of calculation to be done, namely, steady state or transient. In case of the coupled steady state calculation, the neutronic solver and the thermal hydraulic solvers are run sequentially till a converged solution is found. At each time step, the information is exchanged between them; for example, the power is passed from DYN3D-SP3 to the thermal hydraulic solver and they provide the feedback parameters like coolant density or temperature, fuel temperature, and so forth to the neutronic solver for the cross-section update before the neutron transport problem is solved. In Figure5, the different time steps (T0,T1,…,Tn) and flow of information between the involved solvers are shown. In Step 1 the initial TH variables from the previous converged simulation are transferred to the NK-solver without passing through the convergence loop (initialization). Then, DYN3D_SP3 updates the cross-section (Step 2) and calculates the power which is passed over first to the SYRTHES solver (Step 3). In Step 4, the TH calculation is done; both TH codes advance one time step. For the next step, the convergence loop is followed. In Step 5, the information calculated by the TH codes is sent first to the main program. In case of a coupled steady state simulation, the convergence loop is executed after all the information has been collected by the main code. If the convergence criteria are achieved, then the TH parameters are sent to the NK code (Step 6). Step 7 follows which is analog to Step 2.Figure 5
Time advancement description for the codes.Due to the fact that there are three solvers coupled to each other solving the fluid dynamics (NEPTUNE-CFD), heat conduction and transfer (SYRTHES), and neutronic (DYN3D-SP3), the coupled system may need a lot of iterations to converge since the TH parameters are very sensitive to any power changes. Since SYRTHES uses the same time step as the other code solvers but it is not mandatory to do so, this time step can be modified to speed up the convergence. Hence, high values for the time step of SYRTHES (e.g., 104 seconds) were selected in order to neglect the thermal inertia from the energy equation during the heat conduction solution. By this way, the temperature changes in the fuel will rapidly accommodate to the new power distribution.As soon as the NK code changes to a transient simulation mode, the time step of SYRTHES will depend on the NEPTUNE_CFD solver, which defines the time step for all involved codes. During the coupled simulation, the NK code waits until the thermal hydraulic solver finalizes the simulation. The following thermal hydraulic variables are monitored and their relative changes will be used as convergence criteria: Doppler temperature, moderator temperature, and moderator density. In addition, the heat flux at the solid-liquid interface is also checked for as convergence criteria. For this purpose, the relative variation of these parameters between the current and the last time step is calculated. The convergence criterion is reached when the relative variation of these variables is below 10−4. Furthermore, to monitor the evolution of the steady state coupled simulation, the main program writes the convergence parameters at each axial level of the TH discretization, that is, local convergence criteria. When all the axial levels are converged then the information is sent to the NK code to perform another time step calculation.In case of a coupled transient simulation, there is no convergence loop before Step 6 and the TH parameters are sent directly to the NK solver by the main program after they have been collected. Both the TH and the NK solvers are executed in transient mode. To assure that the effects of the thermal inertia are accurately described, SYRTHES uses the same time step as the other solvers. In Figure5, the different steps of the time advancement of the involved codes are indicated. After the cross-sections have been updated in Step 7, the power is predicted by DYN3D-SP3 and it is transferred to the SYRTHES code (Step 8). With the refreshed power distribution the TH code advances one time step (Step 9). Then, the TH feedback parameters are sent to the main program (Step 10) which passes the information to the NK solver (Step 11) without any convergence loop. Based on this new information, this solver can again update the cross-sections before solving the neutron transport problem for a new time step (Step 12).For the information exchange between the involved solvers in the coupled system NHESDYN, the definition of an appropriate spatial mapping is very important. This will be described in the next subsection.


### 2.4. Description of the Spatial Mapping

Several steps must be performed using the MPI-features for the exchange of information of the different solvers considering the peculiarities of each spatial nodalization. For this purpose, the information has to be prepared before and after each data exchange between solvers. One has to make sure that the axial and radial discretization of any problem used by the different solvers are compatible with each other and to assure that the information exchange is self-consistent. For the neutronic SP3 solver, a fuel rod is radially treated as one homogenized cell while in the axial direction it can be subdivided into any number of cells. The SYRTHES solver, on the other hand, can discretize the fuel pin (fuel, gap, and cladding) in different radial rings and axial nodes. Finally, the CFD solver allows a very detailed discretization of the fluid domain including a large number of cells. The spatial refinement of the fluid domain is problem dependent. Those different discretization approaches of the involved codes represent a real challenge to be solved by the coupled system. To avoid mesh dependencies in the communication, a “virtual interface” is created before the information exchange is done between solvers. Doing so, the important values to be passed from one domain to another are approximated to continuous functions. These functions are described by 4th degree polynomial approximations. Due to the simplicity of the case a 4th degree polynomial is enough to describe the axial profiles. In case of setting a more complex problem, for example, different burn up or enrichment, other approaches need to be investigated.Then, the real information exchanged between codes is done via these approximations. By this innovative approach, the amount of information to be treated and sent by the MPI code is considerably reduced. In Figure6, the different solver discretizations are shown. The representation of the mesh for a single pin is illustrated by (a), where the fuel gap and clad are modeled by SYRTHES and the fluid domain is modeled by NEPTUNE-CFD. For each axial level of this mesh, the averaged values of the TH parameters are taken (b); these are moderator temperature and density and Doppler temperature. These TH values are approximated to polynomial functions (c). The DYN3D-SP3 code reads the TH values at the required location by interpolation of the polynomial functions (d). The NK computes the information generating a power map distribution (e). The power values are approximated by a polynomial function (f) and finally these values are the input for SYRTHES.Figure 6
Spatial mapping for a pin: (a) NEPTUNE-CFD/SYRTHES model, (b) averaged TH values at each axial level, (c) polynomial approximation of TH parameters, (d) DYN3D-SP3 model, (e) computed power, and (f) polynomial approximation of the axial power.


## 2.1. General Coupling Approach

A general overview of the data flow between the codes and how it was implemented is described here. First of all, the first step for the communication among the codes is the creation of the MPI groups. A general scheme concerning the group distribution of the different codes is shown in Figure2. This scheme defines the direct dependencies between the two TH codes involved and the NK code with the main program. The amount of information to deal with is high and to work with three different outputs is tedious. To unify the information provided by the codes an extra program is created; this is the main program (NHESDYN). It is programmed in FORTRAN and it was especially developed to control the time synchronization for the coupled solution and to monitor the progress of the coupled simulations.Figure 2
General distribution of the three MPI groups and the codes included within them.NHESDYN generates DYN3D-SP3 as a subprocess (spawn operation) and creates a point to point communication with NEPTUNE-CFD. At the same time the NK code creates a point to point communication with SYRTHES. In the MPI, the point to point communication is a client and a server.In Figure2, the TH tools are considered like a single code. These codes were already coupled by the code developers (EDF and CEA) and the communication between them is done directly with no MPI interaction, for example, to solve the conjugated heat transfer. Additional information is shown in Figure 2 concerning the kind of information transferred (densities, temperatures, etc.) which is classified by tags.In Figure3, a simplified flow chart to understand the function of each module in the coupled solution is given. The NK/TH tools follow a specific order in the information flow while the main code (NHESDYN) controls the times and performs the convergence loops during some stages of the simulation. Considering that the NK code has solved the previous time step, the information (volumetric heat flux) is sent to the SYRTHES code. Then, the TH step is calculated and the conjugated heat transfer through the solid-liquid interface is solved.Figure 3
Simplified flow chart for the coupled NHESDYN system.The information generated takes two different ways. The information from SYRTHES (Doppler temperature) can be transferred directly to DYN3D-SP3 because both are included in the same MPI group. The second way is followed by the information from NEPTUNE-CFD (moderator temperatures and densities) which is sent to the main code. Here, if it applies some operations are performed (e.g., convergence loops for a steady state case). Afterwards, the information is resent to the NK code. By this last operation the NK has enough information to refresh the precalculated cross-sections and solve the transport equation. When the calculation is finished the heat flux is sent again to the TH code which is waiting for the next time step calculation.


## 2.2. Execution of the Coupled Code NHESDYN

In Figure4, the different steps performed during the coupled solution are illustrated. The coupled simulation is started after a converged thermal hydraulic stand-alone NEPTUNE-CFD/SYRTHES run is performed. Since no steady option for NEPTUNE-CFD exists, a null transient simulation is performed for constant boundary conditions. The thermal hydraulic time step is used later on as a common time step for the other solvers during the iterative coupling approach. During the first iteration, NEPTUNE-CFD and SYRTHES provide the initial conditions for the neutron kinetic solver (DYN3D). With the TH information DYN3D calculates the power distribution. This loop runs till the predefined convergence criteria are achieved.Figure 4
NHESDYN main steps of execution.DYN3D-SP3 first performs a steady state run and then the transient calculation is initiated. In this case, both codes run at the same time and the thermal hydraulic boundary conditions can be modified and these modifications can be taken into account by the transient neutronic solver.


## 2.3. Description of the Time Advancement

The coupling scheme depends on the type of calculation to be done, namely, steady state or transient. In case of the coupled steady state calculation, the neutronic solver and the thermal hydraulic solvers are run sequentially till a converged solution is found. At each time step, the information is exchanged between them; for example, the power is passed from DYN3D-SP3 to the thermal hydraulic solver and they provide the feedback parameters like coolant density or temperature, fuel temperature, and so forth to the neutronic solver for the cross-section update before the neutron transport problem is solved. In Figure5, the different time steps (T0,T1,…,Tn) and flow of information between the involved solvers are shown. In Step 1 the initial TH variables from the previous converged simulation are transferred to the NK-solver without passing through the convergence loop (initialization). Then, DYN3D_SP3 updates the cross-section (Step 2) and calculates the power which is passed over first to the SYRTHES solver (Step 3). In Step 4, the TH calculation is done; both TH codes advance one time step. For the next step, the convergence loop is followed. In Step 5, the information calculated by the TH codes is sent first to the main program. In case of a coupled steady state simulation, the convergence loop is executed after all the information has been collected by the main code. If the convergence criteria are achieved, then the TH parameters are sent to the NK code (Step 6). Step 7 follows which is analog to Step 2.Figure 5
Time advancement description for the codes.Due to the fact that there are three solvers coupled to each other solving the fluid dynamics (NEPTUNE-CFD), heat conduction and transfer (SYRTHES), and neutronic (DYN3D-SP3), the coupled system may need a lot of iterations to converge since the TH parameters are very sensitive to any power changes. Since SYRTHES uses the same time step as the other code solvers but it is not mandatory to do so, this time step can be modified to speed up the convergence. Hence, high values for the time step of SYRTHES (e.g., 104 seconds) were selected in order to neglect the thermal inertia from the energy equation during the heat conduction solution. By this way, the temperature changes in the fuel will rapidly accommodate to the new power distribution.As soon as the NK code changes to a transient simulation mode, the time step of SYRTHES will depend on the NEPTUNE_CFD solver, which defines the time step for all involved codes. During the coupled simulation, the NK code waits until the thermal hydraulic solver finalizes the simulation. The following thermal hydraulic variables are monitored and their relative changes will be used as convergence criteria: Doppler temperature, moderator temperature, and moderator density. In addition, the heat flux at the solid-liquid interface is also checked for as convergence criteria. For this purpose, the relative variation of these parameters between the current and the last time step is calculated. The convergence criterion is reached when the relative variation of these variables is below 10−4. Furthermore, to monitor the evolution of the steady state coupled simulation, the main program writes the convergence parameters at each axial level of the TH discretization, that is, local convergence criteria. When all the axial levels are converged then the information is sent to the NK code to perform another time step calculation.In case of a coupled transient simulation, there is no convergence loop before Step 6 and the TH parameters are sent directly to the NK solver by the main program after they have been collected. Both the TH and the NK solvers are executed in transient mode. To assure that the effects of the thermal inertia are accurately described, SYRTHES uses the same time step as the other solvers. In Figure5, the different steps of the time advancement of the involved codes are indicated. After the cross-sections have been updated in Step 7, the power is predicted by DYN3D-SP3 and it is transferred to the SYRTHES code (Step 8). With the refreshed power distribution the TH code advances one time step (Step 9). Then, the TH feedback parameters are sent to the main program (Step 10) which passes the information to the NK solver (Step 11) without any convergence loop. Based on this new information, this solver can again update the cross-sections before solving the neutron transport problem for a new time step (Step 12).For the information exchange between the involved solvers in the coupled system NHESDYN, the definition of an appropriate spatial mapping is very important. This will be described in the next subsection.


## 2.4. Description of the Spatial Mapping

Several steps must be performed using the MPI-features for the exchange of information of the different solvers considering the peculiarities of each spatial nodalization. For this purpose, the information has to be prepared before and after each data exchange between solvers. One has to make sure that the axial and radial discretization of any problem used by the different solvers are compatible with each other and to assure that the information exchange is self-consistent. For the neutronic SP3 solver, a fuel rod is radially treated as one homogenized cell while in the axial direction it can be subdivided into any number of cells. The SYRTHES solver, on the other hand, can discretize the fuel pin (fuel, gap, and cladding) in different radial rings and axial nodes. Finally, the CFD solver allows a very detailed discretization of the fluid domain including a large number of cells. The spatial refinement of the fluid domain is problem dependent. Those different discretization approaches of the involved codes represent a real challenge to be solved by the coupled system. To avoid mesh dependencies in the communication, a “virtual interface” is created before the information exchange is done between solvers. Doing so, the important values to be passed from one domain to another are approximated to continuous functions. These functions are described by 4th degree polynomial approximations. Due to the simplicity of the case a 4th degree polynomial is enough to describe the axial profiles. In case of setting a more complex problem, for example, different burn up or enrichment, other approaches need to be investigated.Then, the real information exchanged between codes is done via these approximations. By this innovative approach, the amount of information to be treated and sent by the MPI code is considerably reduced. In Figure6, the different solver discretizations are shown. The representation of the mesh for a single pin is illustrated by (a), where the fuel gap and clad are modeled by SYRTHES and the fluid domain is modeled by NEPTUNE-CFD. For each axial level of this mesh, the averaged values of the TH parameters are taken (b); these are moderator temperature and density and Doppler temperature. These TH values are approximated to polynomial functions (c). The DYN3D-SP3 code reads the TH values at the required location by interpolation of the polynomial functions (d). The NK computes the information generating a power map distribution (e). The power values are approximated by a polynomial function (f) and finally these values are the input for SYRTHES.Figure 6
Spatial mapping for a pin: (a) NEPTUNE-CFD/SYRTHES model, (b) averaged TH values at each axial level, (c) polynomial approximation of TH parameters, (d) DYN3D-SP3 model, (e) computed power, and (f) polynomial approximation of the axial power.


## 3. Validation of the NHESDYN Coupled System

Since the involved codes are already in the process of validation, the coupled NHESDYN system will be verified and validated by a code-to-code comparison. For this purpose, the DYNSUB system, developed at KIT by coupling DYN3D-SP3 with the subchannel code SUBCHANFLOW [7, 14], was selected. A set of test problems were defined to check the prediction capability of the new NHESDYN system compared to DYNSUB where the thermal hydraulic conditions, for example, coolant mass flow, temperature, and system pressure, were changed in time. The geometrical and material data as well as the thermal hydraulic conditions for the test cases were derived from the OECD/NEA and US NRC PWR MOX/UO2 core transient benchmark [15]. The cross-section data were also taken from the mentioned benchmark.


### 3.1. Definition of the Single Fuel Rod Problem

The test problem consists of one pin surrounded by the coolant with the thermal conditions given in Table1. The fuel is UO2 with Zircaloy as cladding material.Table 1
Nominal conditions for the simulations according to the OECD/NEA and US NRC PWR MOX/UO2 core transient benchmark adapted for a single pin geometry.
Power (kW)
Inlet temp. (K)
Pressure (MPa)
Mass flow (kg/s)
21.68
560
15.5
0.28The dimensions of the fuel rod are exhibited in Table2.Table 2
Rod geometry considered for the validation based on the the OECD/NEA and US NRC PWR MOX/UO2 core transient benchmark.
Fuel radius (mm)
Clad inner radius (mm)
Rod radius (mm)
3.951
4.010
4.583


### 3.2. Modeling Issues for the Code-to-Code Comparison

For a consistent code-to-code comparison between NHESDYN and DYNSUB we made sure that both coupled codes use the same thermophysical material data for the fuel, gap, and cladding, for example, thermal conductivity, density, and the heat capacity as defined in [15], the same water/steam tables, for example, the IAPWS [16], the same gap width (0.05 mm) filled with helium, and the same heat transfer coefficient.In DYNSUB, a constant heat transfer coefficient (HTC) of 104 W/m2K is used. In NHESDYN, the temperature gradient over the gap is calculated by solving locally the conduction within a solid using the thermophysical properties of helium. This is a simplification of the complex heat transfer phenomena that can take place in the gap, normally a combination of conduction, convection, and radiation. The difference between these two solution approaches leads to significant changes in the Doppler temperature prediction. To reduce the difference, the same HTC is used in both codes. After a NHESDYN steady state calculation, an averaged HTC of 11710 W/m2K is calculated based on the heat transfer area, the heat flux, and gradient of the fuel rod gap temperature. This value is used in DYNSUB as the HTC over the fuel rod gap. In Figure 7, the details of the computational domains are illustrated.(a) Plant view of a 17× 17 PWR FA. (b) Plant view of a single pin with the geometric characteristics. (c) 3D view of the numerical domain simulated with the TH tools of NHESDYN.
(a)
(b)
(c)NEPTUNE-CFD describes the fluid domain where the mass flow rate is imposed at the inlet (the bottom of the subchannel) and the pressure is set at the outlet (the top of the subchannel) as boundary conditions.The fluid domain is 12.6 mm in theX and in the Y direction. The walls of the perimeter are set with a symmetry boundary condition. The computational domain is 3.7 m long in the Z direction; the fuel and clad geometrical aspects simulated by SYRTHES are summarized in Table 2.


### 3.3. Discretization Issues in SYRTHES and NEPTUNE-CFD

In the coupled code NHESDYN, the solid domain is solved by SYRTHES while the fluid domain is solved by NEPTUNE-CFD. In Figure8, the meshing for SYRTHES (fuel pin) and for NEPTUNE-CFD (fluid domain) is given. Two meshes were considered for both domains: a coarse mesh (shown in (a) and (b)) and a refined mesh (see (c) and (d)). The coarse mesh has 44 axial levels while the refined mesh consists of 60 axial levels. The pin model in DYNSUB consists of 24 axial levels. In radial direction, reflective boundary conditions are considered for this single fuel rod problem.Fuel pin with a fluid domain with different discretizations: (a) and (b) correspond to a coarse mesh and (c) and (d) to a refined mesh.
(a)
(b)
(c)
(d)


### 3.4. Testing the Mesh Dependence of NHESDYN for Steady State Simulations

Since the modelling of the heat transfer over the fuel gap is very important, a mesh sensitivity study was performed with the goal to find out the most appropriate meshing of both fuel and fluid domains for the coupled simulations. In SYRTHES, a coarser and a refined mesh of the fuel and clad are considered; the number of cells in the gap is 2; see Table3. In R1 and R2 the boundary nodes of the clad are set with the helium properties; in R2 only the nodes within the gap are defined with the helium properties.Table 3
Specifications of the test matrix for the sensitivity study of the influence of a mesh refinement on the prediction of fuel and cladding temperatures.
Case
Mesh type
Gap boundary nodes properties
Run_1 (R1)
Refined
He
Run_2 (R2)
Coarse
He
Run_3 (R3)
Coarse
UO2-ZircaloyThen, three NHESDYN steady sate simulations and one DYNSUB were performed for the three cases illustrated in Table3 to determine the influence of the mesh type on the predicted results, for example, fuel centerline temperature and coolant temperature, which are important feedback parameters between the neutronic and thermal hydraulic solvers. A comparison of the predicted axial Doppler temperature for the cases R1, R2, and R3 with the one predicted by DYNSUB is shown in Figure 9.Figure 9
Comparison of the axial Doppler temperature distribution calculated by DYNSUB and by NHESDYN using three different resolutions (R1, R2, and R3).The Doppler temperature is calculated taking into account the fuel center line and the fuel outer surface. It can be seen that the Doppler temperature for the cases R1 and R2 is very close to the one of DYNSUB (reference solution). For case 3 the Doppler temperature is largely underpredicted compared to the reference value. In addition, the radial temperature profiles along the fuel rod radius predicted by NHESDYN for the three cases (R1, R2, and R3) and by DYNSUB are compared to each other in Figure10. There it can be observed that the radial temperature profiles for cases R1 and R2 are in good agreement with the one calculated by DYNSUB, especially for case R1 (fine grid resolution). Finally, in Figure 11 axial profiles of the averaged cladding temperature for the cases R1, R2, and R3 and for DYNSUB are exhibited. The profiles slightly deviate from each.Figure 10
Radial temperature distribution over the radius of the fuel pin at 1.8 m elevation. Comparison of DYNSUB and different NHESDYN resolutions (R1, R2, and R3).Figure 11
Comparison of the axial evolution of the cladding temperatures for DYNSUB and different NHESDYN resolutions (R1, R2, and R3).The axial coolant temperature distribution calculated for the three cases (R1, R2, and R3) and by DYNSUB is shown in Figure12. There it can be seen that the deviations are negligibly small. They increase slightly with the axial height. It seems that a better match is achieved using a refined mesh. But the differences between the calculations are smaller than one degree. In Figure 13, the evolution of the multiplication factor (keff) during the iteration steps during the coupled NHESDYN steady state simulation is given.Figure 12
Comparison of the axial moderator temperature evolution predicted by DYNSUB and by NHESDYN using three resolutions (R1, R2, and R3).Figure 13
Evolution ofkeff during the iterations of the steady state coupled NHESDYN and DYNBSUB simulations for the three cases (R1, R2, and R3).In a steady state case the coupled system reaches convergence after four NK iterations. The TH needs more time to develop a temperature and density profile after a change in the power; to reach a steady scenario for the TH hundreds of iterations are needed by executing the convergence loop. Case R1 provides the best result regardingkeff compared to the reference value. In addition the coarse mesh lowers CPU times, mainly due to the reduced number of cells for the SYRTHES code in charge of solving the thermal conduction. With the coarse mesh a 60 seconds transient takes 5 minutes of CPU time to be solved; this time is multiplied by a factor of 7 for the refined mesh. The SYRTHES code in this case is the main CPU time consumer due to its serial execution; nevertheless its parallelization is onwards what could lead to a remarkable decrease of CPU time.


### 3.5. Definition of Transient Cases for the Testing of NHESDYN

To test and validate the new NHESDYN coupled system by a code-to-code comparison with the DYNSUB code, the following three transient cases are defined.(i)
Case 1: variation of the coolant temperature at fuel rod inlet.(ii)
Case 2: variation of the coolant mass flow rate.(iii)
Case 3: variation of the system pressure.These transient cases are used to check the robustness and numerical stability as well as the prediction accuracy of NHESDYN.The initial conditions for all transient cases are given in Table1 and they were taken from the previous converged steady state based on the OECD/NEA and US NRC PWR MOX/UO2 core transient benchmark [15].The simulations performed have been tested with the two mesh refinement cases R1 and R2.


#### 3.5.1. Variation of the Coolant Temperature Case

In this postulated transient case, the inlet temperature is decreased 10 degrees from the initial 560 K within few seconds and it remains at this condition for six seconds. Later on, the temperature returns to nominal conditions. Due to the temperature reduction, the moderator density increases leading to a better moderation of the fast neutrons and hence the power increases.In Figure14, the evolution of the inlet temperature as well as the averaged moderator temperature during the transient as predicted by NHESDYN (R1 and R2) and DYNSUB is exhibited. It can be observed that there is no big difference in coolant temperature predicted by NHESDYN compared to the DYNSUB prediction. For the NHESDYN, case R2, the maximal deviation from the reference value (DYNSUB) amounts to 3 degrees. The time evolution of the predicted temperature is faster in the case of DYNSUB than in the case of NHESDYN.Figure 14
Water inlet temperature evolution imposed and averaged water temperature of the domain for different NHESDYN resolutions (R1 and R2) and DYNSUB.The evolutions of the Doppler temperature as predicted by NHESDYN (R1 and R2) and by DYNSUB are compared to each other, Figure15. In general, NHESDYN follows quite well the trend of the reference solution.Figure 15
Comparison of the pin averaged Doppler temperature computed with DYNSUB and with NHESDYN using two different resolutions (R1 and R2).In Figure16, both the fuel rod power and the total reactivity provided by DYNSUB and NHESDYN, cases R1 and R2, are compared with each other. There are very small deviations only at the power peak. Compared with the DYNSUB results, R2 provides 1.55% less power and R1 1.03% more power at the power peak (second 43). Regarding the reactivity, there is a small reactivity peak that appears just before the power returns to nominal conditions for NHESDYN using the nodalization of case R2. This effect is related to the small delay in the TH moderator temperature calculated by NHESDYN.Figure 16
Comparison of the total power and reactivity evolution predicted by DYNSUB and NHESDYN using different resolutions (R1 and R2).


#### 3.5.2. Variation of the Coolant Mass Flow Case

For this postulated transient, it is assumed that the moderator mass flow rate is progressively decreased from the nominal value down to 50% of the nominal conditions, that is, from 0.28 kg/s to 0.14 kg/s within 14 seconds. It remains at this value for four seconds and then it increases with the same change rate till nominal conditions are achieved.By decreasing the moderator flow and keeping all other parameters at the nominal values, the coolant will heat up leading to a reduction of the neutron moderation and hence of the multiplication factor. The moderator temperature predicted at the axial elevation of 2 m by NHESDYN and DYNSUB is illustrated in Figure17. The prediction using the refined mesh (R1) is very close to the reference value (deviation of 0.2 degrees) but in general both NHESDYN predictions follow the same trend as DYNSUB.Figure 17
Comparison of the moderator temperature evolution for two heights in theZ direction (2 and 3.63 m) predicted by DYNSUB and NHESDYN with two different resolutions (R1and R2).In Figure18 a comparison of the Doppler temperature evolution at the axial level of 2 m calculated by NHESDYN (R1 and R2) and DYNSUB is shown. All curves are following similar trends and they are very close to each other. NHESDYN tends to overestimate the temperatures by about two degrees.Figure 18
Comparison of the local Doppler temperature evolution during the transient at the axial height 2 (2 m) predicted by DYNSUB and NHESDYN using two different resolutions (R1 and R2).The predicted total power by NHESDYN and DYNSUB is presented in Figure19. As expected the power decreases due to the reduction of the neutron moderation. All the three power evolutions are very close to each other.Figure 19
Comparison of the power evolution calculated by DYNSUB and NHESDYN using two resolutions (R1 and R2).


#### 3.5.3. The Variation of the System Pressurization Case

A depressurization scenario is defined to study the robustness of NHESDYN when two phase flow conditions appear in the fluid domain. Such conditions are expected to occur when a pipe break in the primary circuit happens. It was assumed that the outlet pressure is decreased from nominal value of 15.5 MPa to 7.1 MPa within 15 seconds (Figure20).Figure 20
Outlet pressure imposed by NHESDYN during the transient. And void fraction evolution predicted by NHESDYN R1 and R2.As it can be seen in Figure20, as soon as the saturation temperature is reached, void fraction (VF) appears in the domain.To make the problem more challenging, a pressure oscillation was defined from second 26 to second 41 as shown in Figure20. In this case, the saturation temperature will also change leading to an oscillation of the void fraction generation.When the void appears, large density gradients in the moderator are expected. The decrease of the density of the moderation leads to a reduction of the neutron moderation and a negative reactivity is produced. On the other hand if the moderation is improved an increase in the power is expected. The refined mesh provided by case R1 has an earlier onset of boiling and generates more void than the coarse mesh provided by case R2. The local density evolution at the axial height of 3.63 m is illustrated in Figure21. Differences in the density evolution between codes are small. Compared to the reference solution, it seems that better agreement is reached with the coarse resolution NEHSDYN (R2) than with the refined one (R1).Figure 21
Comparison of the local moderator density evolution for axial heights in theZ direction (3.63 m) predicted by DYNSUB and NHESDYN using two resolutions (R1 and R2).Before the boiling starts (around 24.5 seconds) the moderator density decreases due to the pressure decrease; this fact provides negative reactivity and the power decreases too; see Figure22.Figure 22
Comparison of the power evolution for the depressurization transient scenario calculated with DYNSUB and NHESDYN using two resolutions (R1 and R2).At that time, a sudden power decrease occurs due to the negative reactivity provided by the moderator density. During the pressure oscillation the void generated oscillates affecting the moderator density and power due to the strong coupling of neutronic and thermal hydraulic processes. The power peaks fit in time with the guided pressure peaks which increase the saturation temperature and the neutron moderation. In Figure23, averaged Doppler temperature as predicted by NHESYN and DYNSUB is given. There is a constant 2-3 degree over prediction by NHESDYN compared to DYNSUB during the oscillation period, where NHESDYN with refined mesh (R1) predicts a Doppler temperature closer to the reference one than with the coarse mesh (R2).Figure 23
Comparison of the total averaged Doppler temperature evolution predicted by DYNSUB and NHESDYN using two resolutions (R1 and R2).


## 3.1. Definition of the Single Fuel Rod Problem

The test problem consists of one pin surrounded by the coolant with the thermal conditions given in Table1. The fuel is UO2 with Zircaloy as cladding material.Table 1
Nominal conditions for the simulations according to the OECD/NEA and US NRC PWR MOX/UO2 core transient benchmark adapted for a single pin geometry.
Power (kW)
Inlet temp. (K)
Pressure (MPa)
Mass flow (kg/s)
21.68
560
15.5
0.28The dimensions of the fuel rod are exhibited in Table2.Table 2
Rod geometry considered for the validation based on the the OECD/NEA and US NRC PWR MOX/UO2 core transient benchmark.
Fuel radius (mm)
Clad inner radius (mm)
Rod radius (mm)
3.951
4.010
4.583


## 3.2. Modeling Issues for the Code-to-Code Comparison

For a consistent code-to-code comparison between NHESDYN and DYNSUB we made sure that both coupled codes use the same thermophysical material data for the fuel, gap, and cladding, for example, thermal conductivity, density, and the heat capacity as defined in [15], the same water/steam tables, for example, the IAPWS [16], the same gap width (0.05 mm) filled with helium, and the same heat transfer coefficient.In DYNSUB, a constant heat transfer coefficient (HTC) of 104 W/m2K is used. In NHESDYN, the temperature gradient over the gap is calculated by solving locally the conduction within a solid using the thermophysical properties of helium. This is a simplification of the complex heat transfer phenomena that can take place in the gap, normally a combination of conduction, convection, and radiation. The difference between these two solution approaches leads to significant changes in the Doppler temperature prediction. To reduce the difference, the same HTC is used in both codes. After a NHESDYN steady state calculation, an averaged HTC of 11710 W/m2K is calculated based on the heat transfer area, the heat flux, and gradient of the fuel rod gap temperature. This value is used in DYNSUB as the HTC over the fuel rod gap. In Figure 7, the details of the computational domains are illustrated.(a) Plant view of a 17× 17 PWR FA. (b) Plant view of a single pin with the geometric characteristics. (c) 3D view of the numerical domain simulated with the TH tools of NHESDYN.
(a)
(b)
(c)NEPTUNE-CFD describes the fluid domain where the mass flow rate is imposed at the inlet (the bottom of the subchannel) and the pressure is set at the outlet (the top of the subchannel) as boundary conditions.The fluid domain is 12.6 mm in theX and in the Y direction. The walls of the perimeter are set with a symmetry boundary condition. The computational domain is 3.7 m long in the Z direction; the fuel and clad geometrical aspects simulated by SYRTHES are summarized in Table 2.


## 3.3. Discretization Issues in SYRTHES and NEPTUNE-CFD

In the coupled code NHESDYN, the solid domain is solved by SYRTHES while the fluid domain is solved by NEPTUNE-CFD. In Figure8, the meshing for SYRTHES (fuel pin) and for NEPTUNE-CFD (fluid domain) is given. Two meshes were considered for both domains: a coarse mesh (shown in (a) and (b)) and a refined mesh (see (c) and (d)). The coarse mesh has 44 axial levels while the refined mesh consists of 60 axial levels. The pin model in DYNSUB consists of 24 axial levels. In radial direction, reflective boundary conditions are considered for this single fuel rod problem.Fuel pin with a fluid domain with different discretizations: (a) and (b) correspond to a coarse mesh and (c) and (d) to a refined mesh.
(a)
(b)
(c)
(d)


## 3.4. Testing the Mesh Dependence of NHESDYN for Steady State Simulations

Since the modelling of the heat transfer over the fuel gap is very important, a mesh sensitivity study was performed with the goal to find out the most appropriate meshing of both fuel and fluid domains for the coupled simulations. In SYRTHES, a coarser and a refined mesh of the fuel and clad are considered; the number of cells in the gap is 2; see Table3. In R1 and R2 the boundary nodes of the clad are set with the helium properties; in R2 only the nodes within the gap are defined with the helium properties.Table 3
Specifications of the test matrix for the sensitivity study of the influence of a mesh refinement on the prediction of fuel and cladding temperatures.
Case
Mesh type
Gap boundary nodes properties
Run_1 (R1)
Refined
He
Run_2 (R2)
Coarse
He
Run_3 (R3)
Coarse
UO2-ZircaloyThen, three NHESDYN steady sate simulations and one DYNSUB were performed for the three cases illustrated in Table3 to determine the influence of the mesh type on the predicted results, for example, fuel centerline temperature and coolant temperature, which are important feedback parameters between the neutronic and thermal hydraulic solvers. A comparison of the predicted axial Doppler temperature for the cases R1, R2, and R3 with the one predicted by DYNSUB is shown in Figure 9.Figure 9
Comparison of the axial Doppler temperature distribution calculated by DYNSUB and by NHESDYN using three different resolutions (R1, R2, and R3).The Doppler temperature is calculated taking into account the fuel center line and the fuel outer surface. It can be seen that the Doppler temperature for the cases R1 and R2 is very close to the one of DYNSUB (reference solution). For case 3 the Doppler temperature is largely underpredicted compared to the reference value. In addition, the radial temperature profiles along the fuel rod radius predicted by NHESDYN for the three cases (R1, R2, and R3) and by DYNSUB are compared to each other in Figure10. There it can be observed that the radial temperature profiles for cases R1 and R2 are in good agreement with the one calculated by DYNSUB, especially for case R1 (fine grid resolution). Finally, in Figure 11 axial profiles of the averaged cladding temperature for the cases R1, R2, and R3 and for DYNSUB are exhibited. The profiles slightly deviate from each.Figure 10
Radial temperature distribution over the radius of the fuel pin at 1.8 m elevation. Comparison of DYNSUB and different NHESDYN resolutions (R1, R2, and R3).Figure 11
Comparison of the axial evolution of the cladding temperatures for DYNSUB and different NHESDYN resolutions (R1, R2, and R3).The axial coolant temperature distribution calculated for the three cases (R1, R2, and R3) and by DYNSUB is shown in Figure12. There it can be seen that the deviations are negligibly small. They increase slightly with the axial height. It seems that a better match is achieved using a refined mesh. But the differences between the calculations are smaller than one degree. In Figure 13, the evolution of the multiplication factor (keff) during the iteration steps during the coupled NHESDYN steady state simulation is given.Figure 12
Comparison of the axial moderator temperature evolution predicted by DYNSUB and by NHESDYN using three resolutions (R1, R2, and R3).Figure 13
Evolution ofkeff during the iterations of the steady state coupled NHESDYN and DYNBSUB simulations for the three cases (R1, R2, and R3).In a steady state case the coupled system reaches convergence after four NK iterations. The TH needs more time to develop a temperature and density profile after a change in the power; to reach a steady scenario for the TH hundreds of iterations are needed by executing the convergence loop. Case R1 provides the best result regardingkeff compared to the reference value. In addition the coarse mesh lowers CPU times, mainly due to the reduced number of cells for the SYRTHES code in charge of solving the thermal conduction. With the coarse mesh a 60 seconds transient takes 5 minutes of CPU time to be solved; this time is multiplied by a factor of 7 for the refined mesh. The SYRTHES code in this case is the main CPU time consumer due to its serial execution; nevertheless its parallelization is onwards what could lead to a remarkable decrease of CPU time.


## 3.5. Definition of Transient Cases for the Testing of NHESDYN

To test and validate the new NHESDYN coupled system by a code-to-code comparison with the DYNSUB code, the following three transient cases are defined.(i)
Case 1: variation of the coolant temperature at fuel rod inlet.(ii)
Case 2: variation of the coolant mass flow rate.(iii)
Case 3: variation of the system pressure.These transient cases are used to check the robustness and numerical stability as well as the prediction accuracy of NHESDYN.The initial conditions for all transient cases are given in Table1 and they were taken from the previous converged steady state based on the OECD/NEA and US NRC PWR MOX/UO2 core transient benchmark [15].The simulations performed have been tested with the two mesh refinement cases R1 and R2.


### 3.5.1. Variation of the Coolant Temperature Case

In this postulated transient case, the inlet temperature is decreased 10 degrees from the initial 560 K within few seconds and it remains at this condition for six seconds. Later on, the temperature returns to nominal conditions. Due to the temperature reduction, the moderator density increases leading to a better moderation of the fast neutrons and hence the power increases.In Figure14, the evolution of the inlet temperature as well as the averaged moderator temperature during the transient as predicted by NHESDYN (R1 and R2) and DYNSUB is exhibited. It can be observed that there is no big difference in coolant temperature predicted by NHESDYN compared to the DYNSUB prediction. For the NHESDYN, case R2, the maximal deviation from the reference value (DYNSUB) amounts to 3 degrees. The time evolution of the predicted temperature is faster in the case of DYNSUB than in the case of NHESDYN.Figure 14
Water inlet temperature evolution imposed and averaged water temperature of the domain for different NHESDYN resolutions (R1 and R2) and DYNSUB.The evolutions of the Doppler temperature as predicted by NHESDYN (R1 and R2) and by DYNSUB are compared to each other, Figure15. In general, NHESDYN follows quite well the trend of the reference solution.Figure 15
Comparison of the pin averaged Doppler temperature computed with DYNSUB and with NHESDYN using two different resolutions (R1 and R2).In Figure16, both the fuel rod power and the total reactivity provided by DYNSUB and NHESDYN, cases R1 and R2, are compared with each other. There are very small deviations only at the power peak. Compared with the DYNSUB results, R2 provides 1.55% less power and R1 1.03% more power at the power peak (second 43). Regarding the reactivity, there is a small reactivity peak that appears just before the power returns to nominal conditions for NHESDYN using the nodalization of case R2. This effect is related to the small delay in the TH moderator temperature calculated by NHESDYN.Figure 16
Comparison of the total power and reactivity evolution predicted by DYNSUB and NHESDYN using different resolutions (R1 and R2).


### 3.5.2. Variation of the Coolant Mass Flow Case

For this postulated transient, it is assumed that the moderator mass flow rate is progressively decreased from the nominal value down to 50% of the nominal conditions, that is, from 0.28 kg/s to 0.14 kg/s within 14 seconds. It remains at this value for four seconds and then it increases with the same change rate till nominal conditions are achieved.By decreasing the moderator flow and keeping all other parameters at the nominal values, the coolant will heat up leading to a reduction of the neutron moderation and hence of the multiplication factor. The moderator temperature predicted at the axial elevation of 2 m by NHESDYN and DYNSUB is illustrated in Figure17. The prediction using the refined mesh (R1) is very close to the reference value (deviation of 0.2 degrees) but in general both NHESDYN predictions follow the same trend as DYNSUB.Figure 17
Comparison of the moderator temperature evolution for two heights in theZ direction (2 and 3.63 m) predicted by DYNSUB and NHESDYN with two different resolutions (R1and R2).In Figure18 a comparison of the Doppler temperature evolution at the axial level of 2 m calculated by NHESDYN (R1 and R2) and DYNSUB is shown. All curves are following similar trends and they are very close to each other. NHESDYN tends to overestimate the temperatures by about two degrees.Figure 18
Comparison of the local Doppler temperature evolution during the transient at the axial height 2 (2 m) predicted by DYNSUB and NHESDYN using two different resolutions (R1 and R2).The predicted total power by NHESDYN and DYNSUB is presented in Figure19. As expected the power decreases due to the reduction of the neutron moderation. All the three power evolutions are very close to each other.Figure 19
Comparison of the power evolution calculated by DYNSUB and NHESDYN using two resolutions (R1 and R2).


### 3.5.3. The Variation of the System Pressurization Case

A depressurization scenario is defined to study the robustness of NHESDYN when two phase flow conditions appear in the fluid domain. Such conditions are expected to occur when a pipe break in the primary circuit happens. It was assumed that the outlet pressure is decreased from nominal value of 15.5 MPa to 7.1 MPa within 15 seconds (Figure20).Figure 20
Outlet pressure imposed by NHESDYN during the transient. And void fraction evolution predicted by NHESDYN R1 and R2.As it can be seen in Figure20, as soon as the saturation temperature is reached, void fraction (VF) appears in the domain.To make the problem more challenging, a pressure oscillation was defined from second 26 to second 41 as shown in Figure20. In this case, the saturation temperature will also change leading to an oscillation of the void fraction generation.When the void appears, large density gradients in the moderator are expected. The decrease of the density of the moderation leads to a reduction of the neutron moderation and a negative reactivity is produced. On the other hand if the moderation is improved an increase in the power is expected. The refined mesh provided by case R1 has an earlier onset of boiling and generates more void than the coarse mesh provided by case R2. The local density evolution at the axial height of 3.63 m is illustrated in Figure21. Differences in the density evolution between codes are small. Compared to the reference solution, it seems that better agreement is reached with the coarse resolution NEHSDYN (R2) than with the refined one (R1).Figure 21
Comparison of the local moderator density evolution for axial heights in theZ direction (3.63 m) predicted by DYNSUB and NHESDYN using two resolutions (R1 and R2).Before the boiling starts (around 24.5 seconds) the moderator density decreases due to the pressure decrease; this fact provides negative reactivity and the power decreases too; see Figure22.Figure 22
Comparison of the power evolution for the depressurization transient scenario calculated with DYNSUB and NHESDYN using two resolutions (R1 and R2).At that time, a sudden power decrease occurs due to the negative reactivity provided by the moderator density. During the pressure oscillation the void generated oscillates affecting the moderator density and power due to the strong coupling of neutronic and thermal hydraulic processes. The power peaks fit in time with the guided pressure peaks which increase the saturation temperature and the neutron moderation. In Figure23, averaged Doppler temperature as predicted by NHESYN and DYNSUB is given. There is a constant 2-3 degree over prediction by NHESDYN compared to DYNSUB during the oscillation period, where NHESDYN with refined mesh (R1) predicts a Doppler temperature closer to the reference one than with the coarse mesh (R2).Figure 23
Comparison of the total averaged Doppler temperature evolution predicted by DYNSUB and NHESDYN using two resolutions (R1 and R2).


## 3.5.1. Variation of the Coolant Temperature Case

In this postulated transient case, the inlet temperature is decreased 10 degrees from the initial 560 K within few seconds and it remains at this condition for six seconds. Later on, the temperature returns to nominal conditions. Due to the temperature reduction, the moderator density increases leading to a better moderation of the fast neutrons and hence the power increases.In Figure14, the evolution of the inlet temperature as well as the averaged moderator temperature during the transient as predicted by NHESDYN (R1 and R2) and DYNSUB is exhibited. It can be observed that there is no big difference in coolant temperature predicted by NHESDYN compared to the DYNSUB prediction. For the NHESDYN, case R2, the maximal deviation from the reference value (DYNSUB) amounts to 3 degrees. The time evolution of the predicted temperature is faster in the case of DYNSUB than in the case of NHESDYN.Figure 14
Water inlet temperature evolution imposed and averaged water temperature of the domain for different NHESDYN resolutions (R1 and R2) and DYNSUB.The evolutions of the Doppler temperature as predicted by NHESDYN (R1 and R2) and by DYNSUB are compared to each other, Figure15. In general, NHESDYN follows quite well the trend of the reference solution.Figure 15
Comparison of the pin averaged Doppler temperature computed with DYNSUB and with NHESDYN using two different resolutions (R1 and R2).In Figure16, both the fuel rod power and the total reactivity provided by DYNSUB and NHESDYN, cases R1 and R2, are compared with each other. There are very small deviations only at the power peak. Compared with the DYNSUB results, R2 provides 1.55% less power and R1 1.03% more power at the power peak (second 43). Regarding the reactivity, there is a small reactivity peak that appears just before the power returns to nominal conditions for NHESDYN using the nodalization of case R2. This effect is related to the small delay in the TH moderator temperature calculated by NHESDYN.Figure 16
Comparison of the total power and reactivity evolution predicted by DYNSUB and NHESDYN using different resolutions (R1 and R2).


## 3.5.2. Variation of the Coolant Mass Flow Case

For this postulated transient, it is assumed that the moderator mass flow rate is progressively decreased from the nominal value down to 50% of the nominal conditions, that is, from 0.28 kg/s to 0.14 kg/s within 14 seconds. It remains at this value for four seconds and then it increases with the same change rate till nominal conditions are achieved.By decreasing the moderator flow and keeping all other parameters at the nominal values, the coolant will heat up leading to a reduction of the neutron moderation and hence of the multiplication factor. The moderator temperature predicted at the axial elevation of 2 m by NHESDYN and DYNSUB is illustrated in Figure17. The prediction using the refined mesh (R1) is very close to the reference value (deviation of 0.2 degrees) but in general both NHESDYN predictions follow the same trend as DYNSUB.Figure 17
Comparison of the moderator temperature evolution for two heights in theZ direction (2 and 3.63 m) predicted by DYNSUB and NHESDYN with two different resolutions (R1and R2).In Figure18 a comparison of the Doppler temperature evolution at the axial level of 2 m calculated by NHESDYN (R1 and R2) and DYNSUB is shown. All curves are following similar trends and they are very close to each other. NHESDYN tends to overestimate the temperatures by about two degrees.Figure 18
Comparison of the local Doppler temperature evolution during the transient at the axial height 2 (2 m) predicted by DYNSUB and NHESDYN using two different resolutions (R1 and R2).The predicted total power by NHESDYN and DYNSUB is presented in Figure19. As expected the power decreases due to the reduction of the neutron moderation. All the three power evolutions are very close to each other.Figure 19
Comparison of the power evolution calculated by DYNSUB and NHESDYN using two resolutions (R1 and R2).


## 3.5.3. The Variation of the System Pressurization Case

A depressurization scenario is defined to study the robustness of NHESDYN when two phase flow conditions appear in the fluid domain. Such conditions are expected to occur when a pipe break in the primary circuit happens. It was assumed that the outlet pressure is decreased from nominal value of 15.5 MPa to 7.1 MPa within 15 seconds (Figure20).Figure 20
Outlet pressure imposed by NHESDYN during the transient. And void fraction evolution predicted by NHESDYN R1 and R2.As it can be seen in Figure20, as soon as the saturation temperature is reached, void fraction (VF) appears in the domain.To make the problem more challenging, a pressure oscillation was defined from second 26 to second 41 as shown in Figure20. In this case, the saturation temperature will also change leading to an oscillation of the void fraction generation.When the void appears, large density gradients in the moderator are expected. The decrease of the density of the moderation leads to a reduction of the neutron moderation and a negative reactivity is produced. On the other hand if the moderation is improved an increase in the power is expected. The refined mesh provided by case R1 has an earlier onset of boiling and generates more void than the coarse mesh provided by case R2. The local density evolution at the axial height of 3.63 m is illustrated in Figure21. Differences in the density evolution between codes are small. Compared to the reference solution, it seems that better agreement is reached with the coarse resolution NEHSDYN (R2) than with the refined one (R1).Figure 21
Comparison of the local moderator density evolution for axial heights in theZ direction (3.63 m) predicted by DYNSUB and NHESDYN using two resolutions (R1 and R2).Before the boiling starts (around 24.5 seconds) the moderator density decreases due to the pressure decrease; this fact provides negative reactivity and the power decreases too; see Figure22.Figure 22
Comparison of the power evolution for the depressurization transient scenario calculated with DYNSUB and NHESDYN using two resolutions (R1 and R2).At that time, a sudden power decrease occurs due to the negative reactivity provided by the moderator density. During the pressure oscillation the void generated oscillates affecting the moderator density and power due to the strong coupling of neutronic and thermal hydraulic processes. The power peaks fit in time with the guided pressure peaks which increase the saturation temperature and the neutron moderation. In Figure23, averaged Doppler temperature as predicted by NHESYN and DYNSUB is given. There is a constant 2-3 degree over prediction by NHESDYN compared to DYNSUB during the oscillation period, where NHESDYN with refined mesh (R1) predicts a Doppler temperature closer to the reference one than with the coarse mesh (R2).Figure 23
Comparison of the total averaged Doppler temperature evolution predicted by DYNSUB and NHESDYN using two resolutions (R1 and R2).


## 4. Summary and Conclusions

In this paper, the main features of coupled system NHESDYN are described in detail. The peculiarities of the communication between the involved solvers in the frame of the MPI implementation are also presented. The testing and validation of the prediction capability of NHESDYN are carried out by a code-to-code comparison with DYNSUB, which has the same neutronic solver but a subchannel code instead of NEPTUNE-CFD/SYRTHES. For this purpose, a single pin problem was defined. Three postulated transient scenarios were defined to check the prediction accuracy of this new coupled system, namely, variation of the coolant inlet temperature, variation of the mass flow rate, and variation of the system pressure. These transient scenarios were then predicted by NHESDYN and DYNSUB using the same material properties for the fuel rod as defined in OECD/NEA and US NRC PWR MOX/UO2 core transient benchmark [15].Based on the code-to-code comparison, it can be stated that NHESDYN is able to predict the behaviour of the fuel rod under all postulated conditions satisfactorily. By this it is demonstrated that the coupling approach and the information exchange among the solvers are consistent and are working properly. The real advantages of the new system NHESDYN compared to DYNSUB can only be shown, when a pin cluster or a fuel assembly is simulated with a very detailed resolution as usual for CFD codes. Hence the work presented here is the first step in the development of NHESDYN as a high fidelity simulation tool. The promising results obtained till now and presented here are very encouraging for the further development and testing of NHESDYN in solving large problems as fuel rod clusters of fuel assemblies of a PWR core. By this way, NHESDYN makes feasible multiphysics simulations describing the TH phenomena in a multiscale sense, for example, at a micro- and mesoscale, and hence local safety parameters such as clad temperature and pin power can be estimated using a refined space discretization if demanded.NHESDYN showed a stable and robust behaviour dealing with rapid boundary conditions changes as defined in the three scenarios which are challenging especially for CFD codes.The steady state analyses provided satisfactory predictions of the main thermal hydraulic parameters with a good agreement with the reference code in both axial and radial spatial profiles. The convergence loop for the steady state works properly and ensures a converged solution as initial condition for transient simulations. The designed main program drives properly the transition steady state to transient.During the transients, a small delay in the evolution of the moderator temperature was calculated by NHESDYN compared to DYNSUB. The reason for it may be the different spatial discretization used in the thermal hydraulic solvers of NHESDYN (NEPTUNE-CFD) and DYNSUB (SUBCHANFLOW).The code-to-code comparison has demonstrated the prediction capability of NHESDYN for both steady state and transient test problems.Finally, the modularity of the developed coupling scheme makes it possible, for example, to replace the serial SP3-transport solver by a parallel Monte Carlo one, so that the coupled system can be run fully in a parallel mode in the frame of high fidelity simulations for reactor design and safety.


## 5. Outlook

Despite the fact that the proof of concept for NHESDYN was very satisfactory, there are several developments that are necessary to fully take advantage of this MPI-based coupled approach. The following investigations should be performed in the near future.(i)
Extension of the mapping scheme to deal with pin clusters or fuel assemblies is done. The radial mapping for these cases must be extended for a consistent consideration of the feedbacks of larger computational domains. Then, additional effects and feedbacks can be analysed such as cross flow and pins with different power.(ii)
Extension of the coupling scheme for the analysis of boron dilution transients in PWR cores or fuel pin arrangements is done.(iii)
Simulation of BWR relevant transient cases for fast (e.g., rod drop accidents) and slow transients such as pressure, coolant temperature, and coolant mass flow rate perturbations is done.The MPI-based implementation of NHESDYN facilitates high fidelity simulations of fuel rod clusters and large problems using a very detailed spatial resolution.Since the SYRTHES solver is not parallel, it limits the overall performance of NHESDYN. Hence NEPTUNE-CFD/SYRTHES should be replaced by commercial CFD codes numerically stable, robust, and more user friendly such as ANSYS CFX [17] or FLUENT [18]; therefore this will allow the analysis of larger problems making use of massive parallel computers.



---

*Source: 294648-2014-05-21.xml*