The reactor dynamics code DYN3D – models, validation and applications

Progress in Nuclear Energy 89 (2016) 170e190

Contents lists available at ScienceDirect

Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene

The reactor dynamics code DYN3D e models, validation and applications Ulrich Rohde*, Soeren Kliem, Ulrich Grundmann 1, Silvio Baier, Yuri Bilodid, Susan Duerigen, Emil Fridman, Andre Gommlich, Alexander Grahn, Lars Holt 2, Yaroslav Kozmenkov, Siegfried Mittag 3 Helmholtz-Zentrum Dresden-Rossendorf (HZDR), Dresden, Germany

a r t i c l e i n f o

a b s t r a c t

Article history: Received 5 October 2015 Received in revised form 18 January 2016 Accepted 12 February 2016 Available online xxx

The article provides an overview of the reactor dynamics code DYN3D. The code comprises various 3D neutron kinetics solvers, a thermal-hydraulics reactor core model and a thermo-mechanical fuel rod model. The implemented models and methods and the capabilities and features of the code are described. Latest developments of models and methods are delineated. An overview on the status of verification and validation is given. Code applications for selected safety analyses are described. Furthermore, multi-physics code couplings to thermal-hydraulic system codes, CFD and sub-channel codes as well as to the fuel performance code TRANSURANUS are outlined. Developments for innovative reactor concepts, in particular Molten Salt Reactor, High Temperature Gas-cooled Reactor and Sodium Fast Reactor are delineated. The management of code maintenance is briefly described. An outlook on further code development is given. © 2016 Elsevier Ltd. All rights reserved.

Keywords: Reactor dynamics Neutron kinetics Thermal-hydraulics model Fuel rod model Code coupling Verification and validation Transient analysis Innovative reactors Code maintenance

1. Introduction The reactor dynamics code DYN3D is a three-dimensional bestestimate tool for simulating steady states and transients of Light Water Reactors (LWRs) and has been developed at the HZDR (Helmholtz-Zentrum Dresden-Rossendorf),4 Germany, and its predecessor organizations for more than 20 years (Grundmann and Rohde, 1996b; Grundmann et al., 2000). Originally, the code was developed for the analysis of reactivity-initiated transients and accidents of Russian VVER-type reactors and recommended by the IAEA as a reference code for the VVER-440/V213 type reactor (IAEA,

* Corresponding author. E-mail address: [email protected] (U. Rohde). 1 Present affiliation: U. Grundmann Physikalische Berechnungen, Dresden, Germany. 2 Present affiliation: TÜV SÜD Energietechnik GmbH Baden-Württemberg, Filderstadt, Germany. 3 Present affiliation: retired. 4 Before 2011, the Helmholtz-Zentrum Dresden-Rossendorf was officially called Forschungszentrum Dresden-Rossendorf (FZD). http://dx.doi.org/10.1016/j.pnucene.2016.02.013 0149-1970/© 2016 Elsevier Ltd. All rights reserved.

1996). The code is undergoing continuous development with respect to the improvement of physical models and numerical methods. It was coupled to different thermal-hydraulic codes and fuel performance codes. DYN3D has become an advanced simulation tool for transients in Light Water Reactors and innovative reactor designs. The user community comprises 16 users in 7 countries including universities, research organizations, commercial enterprises and nuclear regulatory authorities with their technical support organizations. A general overview of the DYN3D code in Section 2 is followed by short descriptions of the implemented models and features in Section 3. Nuclear data treatment is outlined in Section 4. DYN3D's code couplings to thermal-hydraulic system codes, CFD,5 subchannel codes and the fuel performance code TRANSURANUS are presented in Section 5. Comprehensive efforts have been put into verification and validation of DYN3D against numerical benchmarks, dedicated

5

CFD e Computational Fluid Dynamics.

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

experiments on neutron kinetics and thermal hydraulics as well as real plant data. Some examples for verification, validation and application of the code are given in Section 6. In Section 7, an overview of the development of code versions for application to transient analyses for innovative reactor concepts is given. Recently, a quality assurance and code maintenance system was established which is outlined in Section 8. Section 9 provides a summary.

171

Boundary conditions for the reactor core during transient processes, like pressure drop over the core or total core mass flow rate, pressure at core upper boundary, coolant temperature and boron concentration distribution at the core inlet, have to be provided by the DYN3D input (in form of tables in case of DYN3D stand-alone use). In the case that DYN3D is coupled with a thermal-hydraulic system code, these data have to be provided by this code. A general scheme of the DYN3D code is shown in Fig. 1. For analyzing a steady state, there are several options to make the reactor critical:

2. General overview The neutron kinetics model comprises the solution of the threedimensional two-group or multi-group neutron diffusion equations or simplified transport equations applying nodal expansion methods which is specific for the geometry of fuel assemblies. Upscattering is taken into account in the multi-group version. Cartesian as well as hexagonal fuel assembly geometries are treated. Dividing the reactor axially into layers of variable heights, the “nodes” usually represent prisms determined by the shape of the fuel assemblies. Recently, a neutronic solver on trigonal geometry was implemented, which can be used for hexagonal geometry with mesh refinement option, but also for problems without hexagonal symmetry, where different material compositions can be assigned to each trigonal node (Duerigen et al., 2013). DYN3D includes an intrinsic thermal-hydraulic model for oneand two-phase flow in the reactor core and a fuel rod model. Thermal-hydraulic parameters like fuel temperature and moderator temperature need to be known for the estimation of safety criteria to prevent fuel melting or cladding failure. Together with moderator density they are also needed for the determination of the feedback to neutronics. Based on the fission power distribution, thermal-hydraulic core parameters are calculated. The thermalhydraulic model comprises the solution of the balance equations for mass, momentum and energy of a one- or two-phase flow, heat transfer regime mapping and a thermo-mechanical fuel rod model. As a result, the distribution of fuel temperature, coolant temperature and density as well as boron concentration (in PWR) over the core are obtained for accounting the neutronic feedback. The thermal-hydraulic model also yields parameters which are relevant for safety assessment like maximum fuel temperature, fuel enthalpy, maximum cladding temperature, DNB ratio and cladding oxide thickness. An iteration loop is performed between solving neutron kinetics and thermal-hydraulic models until convergence of the feedback parameters is reached. For the solution of the neutronic equations, nuclear cross sections homogenized over the reactor core nodes are required. They depend on the nuclide composition of the fuel in each individual node, which is affected by the initial fuel composition (enrichment, presence of absorbers) and fuel burnup. Furthermore, the cross sections depend on the feedback parameters listed above. Cross sections and other neutronic parameters (e.g. kinetic parameters) are usually calculated by two-dimensional lattice codes, for example HELIOS (Casal et al., 1991). Recently, the Monte Carlo €nen, 2005) is increasingly used for this purcode Serpent (Leppa pose. Neutronic parameters are calculated for different combinations of burnup and feedback parameters and stored in macroscopic cross section libraries which are linked to the code. Cross sections for each individual node and each time step during a transient are then interpolated from these libraries. While burnup is assumed to be constant during a transient, the dependence of the cross sections from the feedback parameters is taken into account and represents the link between thermal hydraulics and neutronics in the above mentioned iteration loop.


Division of multiplication cross sections by keff
Variation of boron acid concentration
Variation of reactor power Moreover, the steady state of a sub-critical reactor with external source can be calculated. For any steady state, the node-wise fuel burnup distribution can be given as input, if it is calculated by external core design codes, or calculated by simulating a power operation history over one or even several fuel cycles. In the latter case, fuel assembly shuffling, load and unload can be taken into account. If a transient calculation is to be carried out, the following external perturbations can be treated:






Movements of single control rods or control rod banks Variation of core coolant inlet temperature Variation of boron acid concentration Changes of core pressure drop or total mass flow rates Changes of pressure

A detailed description of the models and methods can be found in the report (Grundmann et al., 2005). 3. DYN3D models 3.1. Neutronic solvers The neutronic approach implemented in DYN3D is based on nodal expansion methods (NEM). The neutronic solvers comprise transversely integrated nodal methods leading to one- or twodimensional equations. Considering the one-dimensional equations, local polynomial flux expansion ansatzes are used similarly to the nodal expansion methods in Finnemann et al., 1977, Koebke et al. 1984. Regarding the two-dimensional equations of hexagonal and trigonal geometry two-dimensional polynomial expansions are used (Grundmann, 1999; Grundmann and Holstein, 1999; Duerigen, 2013). The neutron fluxes are represented by polynomials uo to seconder and additionally, by exponential functions being the solutions of the homogeneous diffusion equations of each energy group. DYN3D is available for square and hexagonal fuel assembly geometries. The Cartesian method is applicable to square geometries on both coarse and pin-level meshes, while local refinement for selected assemblies based on flux reconstruction can be optionally defined (see Section 3.2). The hexagonal DYN3D approach (Grundmann and Holstein, 1999) can be applied only on assembly level. To ensure the capability of mesh refinement also for hexagonal problems, a trigonalgeometry approach is available in DYN3D (Duerigen et al., 2013). The trigonal model is furthermore advantageous for modeling asymmetric hexagonal fuel assemblies. Diffusion theory is commonly used and most widely accepted for performing full-core reactor calculations. It does not include direction-of-motion variables and requires low computational effort. The SP3 method

172

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

Fig. 1. General scheme of the code DYN3D.

(Gelbard, 1960; Brantley and Larsen, 2000) is able to capture anisotropic transport effects without a drastic increase in computational resources. Both approaches are imbedded in the DYN3D code (Beckert et al., 2007; Beckert and Grundmann, 2008; Duerigen and Fridman, 2012; Duerigen, 2013). An overview of the current DYN3D neutronic solvers is given in Table 1. Coarse mesh means assembly-wise resolution of the neutron flux and power distribution e which is the initial intention of nodal methods e to provide node-averaged distributions with high accuracy and minimum computational efforts. “Fine mesh solution” comprises the possibility of geometrical refinement of the coarse assembly-wise mesh, independent from the inner structure of the

fuel assembly. “Local fine mesh” stands for a pin-wise resolution of the flux and power distribution. In the DYN3D solvers, the principle unknowns of the problem are the coefficients of the polynomials and the outgoing neutron partial currents Jiþ at the faces of the nodal elements. The partial currents are shown in Fig. 2. There are two methods for hexagonal geometries: HEXNEM1 (Fig. 2(b)) and HEXNEM2, the latter including the partial currents at the corner-points (Fig. 2(c)) leading to a higher accuracy than the former, namely for greater diameters of fuel assemblies. For each homogeneous node, local responsematrix equations allow calculating the moments of the nodeinterface outgoing partial currents in terms of the flux ansatz

Table 1 Overview on current DYN3D neutronic solvers. Geometry

Diffusion/SP3 transport

Steady state/transient

Coarse mesh/fine mesh/local fine mesh

Cartesian Hexagonal Trigonal

Yes/yes Yes/no Yes/yes

Yes/yes Yes/yes Yes/no

Yes/yes/yes Yes/no/no Yes/yes/no

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

173

Fig. 2. DYN3D nodal geometries with radial (aed) and axial (e) outgoing partial currents.

coefficients and the incoming partial currents. As incoming partial currents are outgoing partial currents of adjacent nodes, the equations are solved by the inner iteration procedure. Scattering, fission source, and transverse leakage terms represented by local polynomial expansion ansatzes are updated in the outer iteration loop. Eventually, node-averaged fluxes are deduced from the neutron balance equation and successive estimates of the effective multiplication factor are computed. Concerning transient processes, time integration over a neutronic time step is performed with the help of an implicit difference scheme and exponential transformation similar to Finnemann, 1979, Grundmann et al., 2005. The exponents in each node are calculated from the neutron flux behavior in the previous time step. Delayed neutron sources are given from the integration of precursor densities over the time step, assuming exponential behavior of the fission rates. A future DYN3D option will be the local application of transport solvers using unstructured mesh with boundary conditions extracted from the global diffusion solution. An advanced methodology of current coupling collision probability with expansion of flux on orthonormal polynomials is currently under development (Litskevich and Merk, 2014), which can be used for pin power calculation.

macroscopic neutronics data library sets. However, intra-nodal flux reconstruction inside fuel assemblies is an approximation and is not equivalent to the pin-wise flux calculation. A demonstration for the application of pin-wise inner-nodal flux reconstruction is shown in Figs. 3 and 4. The Fig. 3 shows the core configuration of a generic European PWR reactor core. Fig. 4 shows a zoom of the refinement region. A DYN3D calculation was performed at nominal reactor power with partially withdrawn control rods in the left bottom corner of the core. Flux reconstruction was performed in a square region of the core in the vicinity of these control rods. Within the EU funded project NURISP in close co-operation with the Karlsruhe Institute of Technology the DYN3D version integrated into the European NURESIM software platform has been extended by a new tool for automatic mesh generation, data mapping and visualization, which uses directly the calculation results of the flux reconstruction (Gomez et al., 2014). The flux reconstruction method was verified against the MIDICORE benchmark proposed by V. Krysl et al., 2010. MIDICORE is a section of 10 fuel assemblies close to the outer boundary of a VVER1000 core. The reference solution was produced by the fine-mesh

3.2. Flux reconstruction As full core pin-wise calculations require tremendous computation times, local refinement of flux calculations can be useful in regions where significant flux gradients are expected, e.g. in the vicinity of absorber elements. Alternatively to the local mesh refinement option mentioned at the end of Section 3.1, detailed power distribution inside selected fuel assemblies can be provided based on flux reconstruction (Hadek and Grundmann, 1997; Grundmann et al., 2005). Here, twodimensional diffusion equations on the radial direction are solved inside nodes using homogenised cross sections. The neutron flux is represented by trigonometric and hyperbolic functions being the solutions of the two-group diffusion equation. The unknown expansion coefficients are determined from boundary conditions at the node edges and sides which are obtained by the nodal method. Therefore, the number of expansion coefficients which can be defined depends on the number of conditions at the node boundaries, and therefore, on the order of the nodal method and the geometry of the nodes. The functional expansion of the neutron flux inside the nodes, which is obtained by the above described ansatz, is then used in superposition with the pin-wise flux distribution obtained from the lattice code calculations. This pin-wise distribution is stored as additional information together with the cross sections in the

Fig. 3. Visualization of the power distribution in a generic PWR reactor core with inner-nodal power reconstruction in a selected region of the core (normalized).

174

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

Fig. 4. Zoom to the flux reconstruction region e left: mesh refinement, right e pin-wise power distribution (normalized).

diffusion code PERMAK (Aleshin, 2005). The reflector was modeled in real geometry. ADFs and RDFs were used in the DYN3D calculations. Cross sections were produced using the HELIOS code. The 10 assemblies of the MIDICORE sector at the boundary of the core comprise 3 types of fuel elements, two of them containing 9 Gd pins each. The following mean square deviations (in %) between pin-wise power distributions obtained by DYN3D nodal flux reconstruction and the reference solution for the different pin types (Ieremenko et al., 2012) are presented in Table 2: An acceptable error of the reconstructed pin power distribution against the reference fine mesh diffusion solution was reached. The mean square error for all types of fuel pins is less than 2.5%. For pins at the boundary to the reflector, the mean square deviation reaches 3.5%.

3.3. Adjoint flux and reactivity calculation The derivation of the reactor point kinetics equations is based on the weighted integration of the time-dependent diffusion equations over space. Using the steady-state adjoint group fluxes as weighting function, the so-called ‘dynamic reactivity’ can be calculated from changes of the neutron physical constants during a transient. An approximation of the adjoint flux in steady-state is calculated by the same procedure as the calculation of the normal flux. For this purpose, the cross sections are redefined. The adjoint equation is of the same type as the homogeneous steady-state equation with the artificial eigenvalue a instead of keff.. As keff is used for the correction of multiplication, a is close to one. Weighing the different changes of the neutron physical constants originating from the different kinds of perturbations (external perturbation by control rod movement, thermal-

hydraulic feedback effects) with the steady-state adjoint flux and with the actual normal flux, the single contributions to the dynamic reactivity can be estimated by integration over the reactor core. Components of reactivity, such as Doppler effect, moderator temperature and density effects, boron reactivity, control rod reactivity are obtained. During the considered transient, the impact of the different effects on reactivity is demonstrated. The adjoint flux as weighting function is also needed to obtain consistent definitions of the kinetic parameters of the reactor like beff and the effective neutron life time. The reactivity obtained by applying perturbation theory is called the dynamic reactivity. It must be kept in mind, that perturbation theory is valid only for small changes of the spatial shape of the neutron flux in time. An alternative algorithm to calculate reactivity in DYN3D is inverse point kinetics. The integral neutron flux obtained from the 3D neutronic calculations is interpreted as the point kinetics magnitude function, and reactivity is calculated from solving the inverse point kinetics equations. It should be noted, that in general, dynamic reactivity and inverse kinetics reactivity do not necessarily agree during transient processes. The point kinetics is valid, if the shape function of the neutron flux distribution is constant during the transient. This is, in general, not exactly fulfilled during transients.

3.4. DYN3D thermal hydraulics DYN3D allows the simulation of the neutronic and thermalhydraulic core response. A thermal-hydraulic model of the reactor core and a fuel rod model are implemented in the module FLOCAL (Manera et al., 2005) in DYN3D. The fuel rod model is described in Section 3.5. FLOCAL is a 1-D thermal hydraulics code module suitable for the

Table 2 Mean square deviations between reconstructed DYN3D pin-wise power values and reference solution for MIDICORE.

All fuel pins Periphery fuel pins (last row) Gd fuel pins Fuel pins around guide tubes Fuel pins near reflector

FA type 1 (w/o Gd)

FA type 2 (with Gd)

FA type 3 (with Gd)

0.4 0.8 e 0.2 e

1.9 2.3 2.2 1.7 3.5

1.7 2.0 1.9 1.6 3.2

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

simulation of one-phase and two-phase flows in the reactor core of PWRs, BWRs and gas-cooled reactors like HTGRs. The core is modeled as a series of parallel coolant channels associated with one or more fuel assemblies. The parallel channels are coupled by a boundary condition of equal pressure drop over their total length. Four balance equations are used for the coolant: momentum balance, energy balance and mass balance for the two-phase mixture, and a separate mass balance equation for the vapor phase. Coupling between the two mass balance equations is realized by means of additional evaporation and condensation models, allowing thermodynamic non-equilibrium between the two phases. The difference in velocity between vapor and liquid phase is taken into account in the common momentum balance equation by means of a phase-slip model. The slip model replaces the momentum equation for the second phase by providing an algebraic relation for the ratio between vapor and liquid phase velocities. By introducing two-phase mixture parameters, the mass, energy and momentum balance equations are transformed into conservative form. In this representation, the two-phase flow equations correspond formally to the one-phase balance equations with additional terms and different effective velocities for the mass, momentum and energy transport. The non-homogeneity of the flow is reflected in the definitions of the two-phase flow parameters, and the non-equilibrium between the phases is expressed by means of the condensation/evaporation rate and by an additional source term in the energy balance. Numerical methods available for conservative hyperbolic partial differential equations (PDE) can be applied. The set of four PDEs is transformed into a set of ODEs (ordinary differential equations) after subdivision of the coolant channels into a user-defined number of axial nodes. Time-dependent boundary conditions can be given for inlet temperature, core mass flow rate or pressure drop and system pressure. In the case of coupling DYN3D with a thermal-hydraulic system code (see Section 5), these boundary conditions are provided by the system code. The four balance equations are completed by constitutive relations used for the distributed and local frictional pressure drops, for the evaporation/condensation rate and for the phase slip ratio, as well as the thermodynamic state equations for phase temperature and density (depending on specific enthalpy and pressure) in the gas and liquid phases, respectively. In addition to water and steam, thermodynamic properties of sodium and gaseous Helium are implemented. A full set of correlations is implemented for heat transfer regime mapping in the range between one-phase convection in PWR normal operation and post-crisis heat transfer including superheated steam. The code switches automatically between the different heat transfer regimes according to criteria defined for the cladding surface temperature. The occurrence of the heat transfer crisis (either DNB or dry-out) is determined by using different correlations for critical heat flux or dry-out steam quality. After critical heat flux or dry-out has been reached, the heat transfer regime is switched to a transition mode, where the heat flux reduces with increasing cladding surface temperature. After the Leidenfrost temperature has been reached, the model switches to disperse flow or stable film boiling. Besides the balance equations for mass, momentum and energy of the two-phase flow, a conservation equation for the boron concentration in the coolant is solved separately for every single channel. It is taken into account that the boric acid is transported attached to the liquid phase. It is assumed that the boron is not carried by the vapor, because the solubility of boric acid in the vapor is very small. The boron conservation equation is a PDE of hyperbolic type without source and sink terms. In the case of boron dilution transients, the core response depends very sensitively on

175

the spatial- and time-dependent distribution of the boron concentration in the reactor core. Therefore, this distribution should be calculated with minimum numerical error which can be caused by artificial numerical diffusion. To avoid numerical diffusion, a special algorithm of solving the boron transport equation was implement, ed based on a Particle-In-Cell (PIC) method (Rohde and Lucas, 1997). The PIC algorithm can be activated instead of a secondorder finite difference scheme. The distribution of coolant temperature and boron acid concentration at the core inlet can be obtained from given primary circuit cold leg parameters using an intrinsic coolant mixing model for the downcomer and lower plenum of the reactor. Different options for modeling coolant mixing are available. Besides the simplest cases of homogeneous mixing and the absence of any mixing, an integrated analytical coolant mixing model can be used. It is based on the analytical solution for the diffusion in a developed flow in the potential flow approximation (Prasser, 1983; Ackermann and Draeger, 1987). The model can be adapted to realistic cases by variation of the turbulent Peclet number. It was validated against CFD calculations and mixing experiments (Kliem et al., 1999). Using the results of CFD calculations can be an alternative to obtain realistic distributions of coolant temperature and boron concentration at the core inlet. An iteration between thermal hydraulics and neutron kinetics is performed in steady state calculations as well as during transients. Analyzing transients the thermal-hydraulic time step can contain several neutron kinetic steps. If one considers transients at hot zero power (HZP) and the nuclear power undergoes rapid changes, the neutron kinetic step has to be small in comparison to thermalhydraulic step which reflects the slower changes of thermalhydraulic properties. 3.5. Fuel rod model DYN3D contains a 1D fuel rod behavior model including heat conduction in the fuel and cladding, and heat transfer through the gas gap between them (Rohde, 2001). Besides the fuel temperature, which is a reactivity feedback parameter for the neutron kinetics, and its radial distribution, the model provides data for safety estimation like fuel enthalpy, maximum cladding temperature, cladding oxide layer thickness and cladding mechanical failure. More precisely, the fuel rod model is comprised of the solution of the one-dimensional heat conduction equation for cylindrical fuel rods consisting of fuel and cladding. Thermo-physical properties (heat conductivity, specific heat) of fuel and cladding are temperature dependent; a set of correlations describing these dependencies is implemented in the code. Von-Neumann boundary conditions for the heat conduction equation are set at the fuel and cladding outer boundaries. While the heat transfer coefficient from cladding to coolant is determined in the thermal-hydraulic model, the heat transfer coefficient in the gas gap between fuel and cladding is calculated by a mechanistic model which is based on the URGAP model (Lassmann and Hohlefeld, 1987). It takes into account conduction in the filling gas, conduction by mechanical contact between fuel and cladding as well as by radiation. The general idea of the gas gap behavior modeling in DYN3D is that the parameters for the stationary reference state (e.g. geometrical gap width, gas pressure and gas composition) have to be provided as input from a detailed fuel performance code like TRANSURANUS (Lassmann, 1992), while changes of the fuel rod during the transient, in particular the gas gap parameters, are estimated by the transient fuel rod model in DYN3D. In this regard, changes of gas gap geometryical width, filling gas pressure and temperature during transients are taken into account

176

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

in DYN3D. Gas composition in the gap is assumed to be constant, i.e., no transient fission gas release is considered. For the estimation of the gas gap width, thermal expansion, elastic and plastic deformation of the cladding as well as thermal expansion and elastic deformation of the fuel are taken into account. The elastic deformation of the cladding is obtained based on the thin shell approximation. Plastic deformation of the cladding is modeled according to creeping law. In the case of gap closure, contact pressure is calculated and used to estimate contact conduction. Furthermore cladding mechanical failure is indicated conservatively if the mechanical stress will exceed the yield strength. Metal-water reaction of the cladding material at high temperatures is modeled using an Arrhenius equation. An additional heat source in the cladding from this exothermic reaction is taken into account. Therefore, the oxide layer thickness is obtained as a safety relevant parameter. The FLOCAL fuel rod model was validated against power pulse experiments at the Russian IGR and the Japanese NSRR reactor for fresh fuel rod probes (Rohde, 2001). For modeling the fuel behavior in more detail, in particular at higher burnup, fuel performance codes like TRANSURANUS are needed. For this reason, the two-way coupled code DYN3DTRANSURANUS has been developed recently (Holt et al., 2014, 2015) (see Subsection 5.3). 3.6. Burnup capabilities DYN3D can perform not only steady-state and transient calculations, but also fuel cycle simulation, taking into account fuel element shuffling operations. A core burnup calculation is done by putting a burnup loop around the stationary calculation kernel of the code: - The steady-state problem with the corresponding search for criticality by division of multiplication cross sections by keff, variation of boron acid concentration or variation of reactor power is solved for the beginning of each burnup time step. - The obtained neutron fluxes and power density distributions are assumed to be constant during a (small) burnup step. Power densities are used to calculate the increment of burnup for each node and the neutron fluxes e to obtain new distributions of the reactor poisons 135Xe and 149Sm, their precursors 135I and 149Pm and the spectral history indicators 239Pu and 239Np at the end of the burnup step. - Applying time steps that are small enough and tuned to the operation history, sufficient accuracy is reached in reasonable computing time; no predictor/corrector methods have been applied. The burnup calculation proceeds step by step until the end of cycle, or a desired burnup state is reached. The DYN3D input options allow the simulation of real fuel cycles using the following power plant data: variation of reactor power, coolant mass flow, inlet temperature and control rod movements during the cycle. The reactor poisons distributions can be dynamically calculated from the values of the previous time step by solving the Bateman equations. Alternatively, “equilibrium” values corresponding to the actual power distribution can be used. 3.7. Reactor poison dynamics To model the poison-related reactor dynamics, in particular xenon oscillations, the balance equations for the 135Xe, 149Sm, 135I and 149Pm concentrations are solved. The same approach as for the burnup calculation is used: The neutron fluxes at the beginning of the time step are obtained in the

steady-state calculation and assumed to be constant during a time step; new values of poisons concentrations are calculated from the Bateman equations using steady-state fluxes. 3.8. Decay heat model A model for calculating the space dependent decay heat power derived from the operational history as well as during the transient is integrated in DYN3D (Grundmann et al., 2005). The decay heat model takes into account contributions from fission products, the decay of actinides and other nuclides as well as of fission products that have captured neutrons. The decay heat calculation is based on the German normative document DIN Norm 25463, Part 1 (DIN90). It is suitable for decay heat calculation for UOX fuel in light water reactors with a maximum initial enrichment of 4.1%. The contribution of fission products to decay heat power is calculated from the individual contributions of four fissile isotopes 235 U, 238U, 239Pu and 241Pu. Decay heat power delivered by each of these isotopes is divided into 24 groups, each having a characteristic decay constant (DIN Norm 25463, 1990). Decay heat power is calculated individually for each node in order to take into account its spatial dependence. The decay heat from fission products is taken into account using the previous reactor operation parameters (“history”) as well as the time-dependent power release during the transient. Contribution from the power released during the transient to the total decay heat is usually small in comparison with the decay heat from the power history. The contribution of groups with longer half-lives can be neglected during the transient. Therefore only the first 10 groups with the shortest half-lives are taken into account. It was verified that this approximation is sufficient for time intervals less than 1000 s. The contribution from the power history can be taken into account by providing a power operation history diagram (operation at different power levels over certain time intervals) as input data or in a simplified way by assuming operation before shutdown or transient at the constant power level of the initial steady state. The neutron capture contribution to decay heat power is mainly caused by the actinides 239U and 239Np. Balance equations for the generation and decay of these nuclides are solved. Contributions to the decay heat power from neutron capture by other nuclides and by fission products are given by time tables in (DIN Norm 25463, 1990). 4. Treatment of neutronic data 4.1. Description of the neutronic data Nodal methods as implemented in DYN3D utilize pin- or assembly-homogenized few-group cross sections (XS). This includes neutron diffusion coefficients, reactions cross sections, as well as reactor poisons and delayed neutron data. In order to reduce homogenization errors, flux Assembly Discontinuity Factors (ADF) (Smith, 1986) in the radial direction can be employed by DYN3D (Grundmann et al., 2005; Mittag et al., 2003). In order to account for strong axial heterogeneities, ADFs in the axial core direction were also implemented in the code (Fridman et al., 2013). Furthermore, DYN3D is capable of iteratively generating and applying superhomogenization (SPH) factors (Hebert and Kavonoki, 1981; Grundmann and Mittag, 2011). XS data are generated by lattice codes for each homogenized material in the core and compiled in data libraries that can be used as DYN3D input. These data libraries represent the dependence of the material's neutronic properties on burnup and reactor

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

operation conditions such as fuel temperature, moderator temperature, boron concentration in the moderator, and moderator density. The actual distributions of burnup and operational parameters are calculated by DYN3D and then used to estimate local XS for each node based on the data from an XS library. Different types of XS libraries can be used as DYN3D input: Parametrized Data Libraries e The XS library contains a neutronic data set for one reference state point (core-averaged temperatures and densities) per burnup step. The dependences of the XS on operational parameters are described by polynomials. Polynomial coefficients dependent on burnup are also stored in the library. Multi-Dimensional-Table Libraries e The data library contains pre-calculated XS for all combinations of operational parameters. Actual local properties are obtained by multi-dimensional interpolation between closest points in the table. Approximation-and-Interpolation Libraries e a combination of approximation by polynomials and interpolation in fewdimensional tables. The library contains XS for a table of reference state points (e.g. dependence on moderator temperature and density), while the dependences on other operational parameters (fuel temperature, boron concentration, etc.) are described by polynomials. 4.2. Spectral history effects The above mentioned XS libraries describe the change of fuel properties during depletion by a dependence on a burnup parameter measured in terms of produced energy, MWd/kgHM. However, in a reactor core, fuel depletion happens under conditions varying in space and time. Different operation parameters lead to different neutron spectra, depending on the coordinates within the core (node locations) and the degree of depletion. The nature of the spectral history effects is the difference in the actinides build-up and uranium depletion rates that change with the spectrum. Fuel, depleted to a certain burnup in different spectral conditions, accumulates a difference in nuclide content. Thus, cross sections (XS) for fuel with the same burnup level could differ even if the instantaneous local conditions are the same. The most pronounced history effect in LWRs is the moderator density effect: the water density in the upper part of a core is significantly lower than in the lower part during the whole reactor operation, so the lower part of fuel assembly is depleted in a softer neutron spectrum which leads to a significant difference in the accumulated Pu, minor actinides and 235U content. Spectral history effects are taken into account in DYN3D by the correction for nodal cross sections which is proportional to the deviation of the nodal (local) 239Pu concentration from the nominal value in the lattice code depletion calculation (Bilodid and Mittag, 2010). Local 239Pu concentration for each node is calculated by DYN3D solving the Bateman equations for 239Pu and its precursor 239 Np in a fuel cycle simulation. Some additional data is necessary for a correction calculation: correction coefficients, nominal 239Pu concentration and microscopic cross sections for 238U, 239Np and 239 Pu. All these data is provided to DYN3D in a “historical sub-library” e an additional input file which is generated once for each fuel type by a lattice code. The application of the correction method significantly improves calculation results (Bilodid et al., 2012, 2014). 5. DYN3D's multi-physics coupled code systems To model the complex processes in nuclear power reactors occurring during transients, the close interaction between neutron kinetics and thermal hydraulics has to be modeled adequately. This has to be done by a tight coupling of the corresponding models.

177

Coupling core models comprising 3D neutron kinetics with thermal-hydraulic system codes became industrial practice. Such coupled codes can simulate the interaction between core and nuclear plant behavior and have limited possibilities for modeling 3D effects in the core thermal hydraulics by applying the technique of cross connection objects or coarse mesh nodalization. Engineering coefficients to model cross flows are required. Advances have been achieved by applying sub-channel models or fully 3D thermal-hydraulic modeling based on the porous body approach. These models provide more detailed and more reliable information on core thermal hydraulics but still need closure laws like turbulent mixing models or friction correlations to be derived from experiments or detailed CFD simulations. Generic CFD modeling of the coolant flow in reactor cores based on first principles is at the fore-front-end of thermal hydraulics. To adequately treat the near-wall behavior of the flow, a spatial resolution much finer then fuel rod sub-channels is required. Moreover, in the frame of the so-called “numerical nuclear reactor” approach, CFD models have been coupled with fine resolution neutron transport methods (Kochunas et al., 2012). These developments are of great methodological interest, but not practicable for applications due to tremendous computer resources required. Moreover, due to a tendency to increase the discharge burnup of nuclear fuels over the last decades, it became more and more important to properly account for fuel behavior in safety analysis. Therefore, DYN3D was coupled with the fuel performance code TRANSURANUS (Lassmann, 1992). 5.1. Coupling of DYN3D to thermal hydraulic system codes Internal, external and parallel coupling techniques (Grundmann et al., 1995, 2005; Kozmenkov et al., 2007) were used to integrate DYN3D with the thermal hydraulic system code ATHLET (Austregesilo et al., 2006), originally containing 0D (point) and 1D (one-dimensional) models of neutron kinetics. Each of these coupling techniques has its own advantages and disadvantages, depending on the problem to be solved. In the case of internal coupling, only the neutron kinetics model of DYN3D is integrated with the system code, while the core thermal hydraulics (as well as the thermal hydraulics of all other reactor components) is calculated by ATHLET. This is the most consistent approach which, however, requires significant modifications of the coupled codes. Within the external coupling approach all physical phenomena in the core, including its thermal hydraulics, are completely modeled by DYN3D, while ATHLET simulates the rest of reactor systems. This technique requires exchange of the core inlet and outlet boundary conditions between thermal-hydraulic models of DYN3D and ATHLET, which is relatively easy to implement. It helps to reduce the cost of simulations in the cases with multiple core channels, but under conditions of a strong feedback between thermal hydraulics and neutronics (e.g., in boiling water reactors) this can lead to numerical instabilities. In the third type of coupling, both ATHLET and DYN3D thermalhydraulic core models are running in parallel. ATHLET calculates the behavior of the whole plant and provides the core boundary conditions (inlet and outlet pressures, inlet coolant enthalpy) to DYN3D. Using these boundary conditions and its own thermalhydraulic model, DYN3D calculates the total core power and transfers it to ATHLET. A parallel way of coupling joins the advantages of internal and external coupling (numerical stability and simplicity of implementation), but inconsistencies might occur between two different thermal-hydraulic models of the core. The resulting coupled code DYN3D/ATHLET has advanced capabilities of modeling both steady-state spatial distributions of the core power and their evolutions during different kinds of reactor

178

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

transients. Initial stationary core power distributions can be disturbed by changes in reactor loop mass flow rates and/or temperatures, by relocations of low-temperature/diluted-boron water slugs within the primary system, or by movement of control rods. It is also highly useful in the analyses of possible reactor instabilities. Big efforts were put into the verification and validation of the coupled code system DYN3D/ATHLET. An overview can be found in Kozmenkov et al. 2015. The coupled code system DYN3D/ATHLET was used as a powerful tool for the analysis of complex accident scenarios like ATWS (Kliem et al., 2009; Mittag and Kliem, 2011) or boron dilution scenarios (Kliem and Rohde, 2007). Figs. 5 and 6 show the results of the application of DYN3D/ATHLET to the ATWS transient: “Complete failure of main feed water supply”, where the power corresponds tightly with the primary circuit parameters, especially the pressure. An internal coupling of DYN3D with the RELAP5 code (NUREG, 2001) is also available (Kozmenkov et al., 2001; Kliem et al., 2006). 5.2. Coupling of DYN3D to sub-channel and CFD codes Even the basically one-dimensional codes like RELAP5 and ATHLET make it possible to simulate certain components that have multidimensional hydraulics by including cross-flow junctions connecting pairs of control volumes of two parallel vertical or horizontal pipes, for example in the core or downcomer region of a PWR. If a series of parallel pipes is interconnected by the ATHLET cross connection objects, the individual junctions of adjacent cross connection objects may exchange momentum. An alternative approach is the coupling of sub-channel codes with neutron kinetic models. Sub-channel codes are designed to describe the flow in rod bundles, where the flow along the rods is the dominating one. Flow components in directions across the bundle are considered to be of second order. The coupled code system DYNSUB (Gomez-Torres et al., 2012) was developed by KIT comprising DYN3D and the sub-channel code SUBCHANFLOW (Imke et al., 2010). In DYNSUB, an internal coupling approach is employed. DYNSUB provides coupled solutions for square fuel lattices and allows the direct calculation of pin level safety parameters. Multi-group pin level solutions of the neutron transport problem using the diffusion approximation or simplified transport (SP3) are available. New possibilities for code coupling were opened in building the European code platform NURESIM. The platform is based on the open source software SALOME. Mesh interpolation modules have

been developed which allow for the embedding of a region with a higher level of modeling resolution (hot-channel), as the consistent interpolation of the neutron kinetic and thermal-hydraulic solution fields is required by the coupling. For VVER applications, the import and export features of the Data Exchange Model (MED) of SALOME have been extended for hexagonal geometry. Recently, within the European project NURISP, the DYN3D code was coupled with FLICA4 (Kliem et al., 2011). FLICA-4 is a fully 3D two-phase flow model which is not based on the assumption of pre-dominant flow along the bundle as it is assumed in most sub-channel codes (Toumi et al., 2000). Therefore, significantly asymmetric flow situations, even with almost stagnant flow in some regions of the core, can be simulated. Application Programming Interfaces (APIs) were developed for the coupling. DYN3D/FLICA coupling was validated against steady-state results from the OECD Main Steam Line Benchmark (Gommlich et al., 2010). Furthermore, within the mentioned NURISP project, a boron dilution benchmark was defined and solved by means of DYN3D/FLICA (Jimenez et al., 2015). A promising approach is to model the reactor core in CFD not at individual sub-channel discretization level, but based on the porous body approximation. This approach was used to implement a coupling between CFD and the DYN3D neutron kinetics (Kliem et al., 2011). The fuel rod lattice of the core is represented by a porous body with anisotropic properties with regard to porosity and friction forces. The physical interface between neutron kinetics and the flow model is the heat released in coolant per unit volume. It comprises the heat transfer from the heated surface (cladding) to the coolant as well as the heat released directly in the coolant by gamma-heating. This approach avoids detailed modeling of the heat transfer from fuel to coolant which is presently not yet solved, in particular for two-phase heat transfer regimes. Presently, the coupling is accomplished for steady-state and transient conditions, but limited to one-phase flow. Test cases covering mini and full core geometries (for PWR and VVER reactors), control rod movement and overcooling transients are presented in Grahn et al. 2015. As an example, power and coolant temperature distribution in a PWR core with 193 fuel elements calculated by ANSYS CFX/DYN3D are shown in Fig. 7. 5.3. Coupling to the fuel performance code TRANSURANUS Besides thermal fluid dynamics codes, a coupling of DYN3D with the fuel performance code TRANSURANUS (Lassmann, 1992) was developed recently (Holt et al., 2014, 2015, 2016). TRANSURANUS is

Fig. 5. Primary circuit pressure during an ATWS transient.

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

179

Fig. 6. Core power behavior during an ATWS transient.

a well validated code widely used by research organizations, safety authorities, and industry. In the frame of this development, a generalized coupling interface was created which can be used also for the coupling of TRANSURANUS with other core simulators, thermal-hydraulic system codes or sub-channel codes. The general TRANSURANUS coupling interface disposes of the following features (Holt et al., 2015):
One-way and two-way coupling. For the two-way coupling, TRANSURANUS replaces the call of the simplified DYN3D fuel behavior model, and is part of the iteration process in each timestep in DYN3D. In this case, TRANSURANUS obviously influences the entire simulation. For example, the transfer of the fuel temperature calculated by TRANSURANUS to DYN3D has direct impact on the Doppler reactivity effect simulated by the neutronics. For the one-way coupling approach, TRANSURANUS has no feedback to the simulation carried out by the DYN3D code.


Flexible time-step length control and numerical stability
Freely selectable level of detail in fuel rod behavior modeling (e.g. from complete thermal-mechanical analysis to thermal analysis with given heat transfer coefficient in the gap)
Pre- and post-calculations of fuel behavior during irradiation before and after the transient (with TRANSURANUS standalone)
Automatic switch from steady-state to transient conditions
Potential for parallelization
Coupling at pin level, at assembly level as well as a mixture of both The general TRANSURANUS coupling interface together with the modified TRANSURANUS version can be used for the simulation of several scenarios in various reactor designs (Holt et al., 2015). As a first application of the general TRANSURANUS coupling interface, the reactor dynamics code DYN3D was coupled at assembly level in order to describe the fuel behavior in more detail. This means that, identical to DYN3D, one representative/average

Fig. 7. Axial temperature distribution in the middle plane of the reactor core (left) and radial power distribution in the middle of the core (right) calculated by ANSYS CFX/DYN3D.

180

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

fuel rod is calculated in the coupled code system for each 1D vertical cooling channel (a 1D vertical cooling channel here represents one fuel assembly). In general, the coupled neutronics/thermal hydraulics code, in our case DYN3D, provides time-dependent rod power and thermal hydraulics conditions to TRANSURANUS. In one-way coupling, there is no feedback from TRANSURANUS to DYN3D, but the input parameters for the fuel performance analysis like e.g. the linear heat rate, are provided to TRANSURANUS online in the course of the transient. For two-way coupling, TRANSURANUS in turn transfers parameters like fuel and cladding temperatures back to DYN3D. The fuel temperature affects directly the neutron kinetics through the Doppler temperature, while the cladding temperature influences the thermal hydraulics through cladding-to-coolant heat transfer in DYN3D (Holt et al., 2015). A first application of the DYN3D-TRANSURANUS coupling was performed analyzing a control rod (CR) ejection scenario for a generic core of German PWR KONVOI (Holt et al., 2014). This design is characterized by a thermal power of 3750 MW, four loops and 193 fuel assemblies of the type 18  18e24. Thence 193 fuel rods were calculated with TRANSURANUS in the coupled code system, each representing the representative/average fuel rod of one fuel assembly. Comparing one-way and two-way approaches, no relevant feedback on the neutron kinetics and thermal hydraulics could be found in spite of the greater degree in detail in online fuel behavior modeling for the CR ejection scenario. Fig. 8 shows the behavior of reactor power calculated by one-way approach (values taken from the DYN3D part in the coupled code system) and twoway approach. The node-wise differences in maximum centerline fuel temperature amount between 4 K and 19 K, comparing both coupling approaches (Holt et al., 2014). Furthermore, the maximum clad surface temperature are almost identical due to absence of film boiling. Moreover, a more challenging boron dilution scenario was analyzed where relevant differences in safety relevant parameters, in particular maximum cladding temperature, occurred (Holt et al., 2016). This analysis was performed for the same core loading of the German Konvoi design as in the case of the CR ejection scenario, hence again 193 fuel rods were calculated with TRANSURANUS. In this scenario, all CR are fully inserted into the core excluding the most effective CR, which is fully withdrawn assuming a single failure. The reactor is assumed to be in a sub-critical shut down state after a steam generator tube break incident. In the course of this scenario, the coolant pressure in the primary circuit will arrive at a lower level than the pressure at the secondary side of the steam

Fig. 8. Total nuclear power over time for CR ejection starting from 30% of nominal reactor power. Taken from Holt et al. 2014.

generator. Hence, un-borated feedwater can penetrate from the secondary to the primary side in the affected loop. If then the main coolant circulation pump as first one is switched on in the affected circuit, the un-borated water will be transported to the reactor core. On the way to the core it is mixed with highly borated water in the reactor pressure vessel. The lower boron concentration in the mixed coolant leads to a positive reactivity insertion into the core which over-compensates the initial sub-criticality and induces a sharp rise of reactor power. This power rise will then be compensated by negative reactivity feedback. As input for boron concentration at each fuel assembly inlet, measured values from the Rossendorf Coolant Mixing Model (ROCOM) facility were used (Kliem, 2010). The maximum total reactivity introduced by the boron dilution achieves 2.07 $, resulting in a super-prompt critical power excursion, and the maximum total nuclear power reaches 32.144 MW in the two-way approach. The differences in the global reactor parameters are almost negligible between two-way and one-way approach. However, the picture will change if we look for local parameters in single nodes of the fuel assemblies. More precisely, in a later phase of the transient, the thermal hydraulics can be affected strongly even in fresh fuel assemblies, where film boiling appeared in one node in the two-way approach in spite of no onset of film boiling in the one-way approach (Holt et al., 2016). According to this node with the highest surface cladding temperature in Fig. 9, a higher heat transfer coefficient in the gap was calculated by TRANSURANUS leading to a higher heat flux, and therefore, to exceeding the critical heat flux in the two-way approach. The heat transfer coefficient in the gas gap changes with burnup and fuel rod enthalpy. This can be correctly modeled only with the help of detailed fuel performance codes like TRANSURANUS. A detailed analysis of the differences can be found in Ref. Holt et al., 2016. For the boron dilution transient, an additional set of calculations was accomplished performing a theoretical study in order to allow analyzing nodes in which film boiling occurs in both approaches of the coupled system DYN3D-TRANSURANUS (Holt et al., 2016). For this purpose, the ratio between the CHF and the actual heat flux, at which DNB is assumed, was increased from 1.0 (best estimate value) to 2.0 (conservative value). This value is an optional input parameter in DYN3D. Fig. 9 shows differences in node-wise clad surface temperatures between the one-way and the two-way approaches in the calculation with conservatively increased critical

Fig. 9. Maximum node clad surface temperature calculated by the two-way coupling approach versus the corresponding value of the one-way coupling approach for each node (red points). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Taken from Holt et al. 2016.

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

DNB ratio. In five nodes we observed cladding temperatures much higher than 350  C indicating post-crisis heat transfer conditions (film boiling). These differences can be mainly explained by feedback from TRANSURANUS cladding temperatures to thermal hydraulics in DYN3D which is taken into account only in the two-way approach. Thus, the one-way coupling approach can produce nonconservative results with respect to safety assessment for cladding temperature, and an advanced coupling between neutron kinetics, thermal hydraulics and detailed fuel behavior modeling is potentially required for reliable safety analysis. Nevertheless, a conclusion about the conservatism of reactor safety analysis cannot be drawn from the cladding temperature alone: Higher clad surface temperatures caused by a higher heat transfer in the gap can lead to lower fuel enthalpy. In contrary, a higher fuel enthalpy can lead to more distinct PCMI with the increased potential for fuel rod failure (Holt et al., 2016). The CPU times required to perform the transient calculation (with inclusion of the time to search for starting conditions defined by flux distribution and temperature distribution) amount 188.0 min for DYN3D standalone, 752.8 min for the one-way approach of DYN3D-TRANSURANUS and 1444.0 for the two-way approach of DYN3D-TRANSURANUS. Since the numerical performance of DYN3D-TRANSURANUS was fast and stable even for these extreme transient conditions. It is therefore concluded that the coupled code system can improve the assessment of safety criteria, at a reasonable computational cost. 6. Verification, validation and applications 6.1. Code verification For code verification, i.e. for proving the numerical correctness of the solutions or the numerical error of approximate methods against exact solution, mathematical benchmarks were used. The capability of DYN3D has been proven via a series of benchmarks. The DYN3D code was verified and validated against benchmarks for PWR, BWR and VVER type reactors. Numerical problems were solved for steady-state and transient tasks within two-group and multi-group diffusion theory as well as SP3 transport approximation. With respect to geometry, square, hexagonal and trigonal nodalization were considered. DYN3D solutions were compared with reference Monte Carlo or deterministic transport solutions. If no reference solution was available, intercode comparisons were performed. In some transient benchmarks, not only pure neutron kinetics was considered, but thermal hydraulics was included into the comparisons to verify reactor dynamics modeling. In particular, mathematical benchmarks were considered within the programs of OECD/NEA, Atomic Energy Research (AER) association on VVER reactor physics and reactor safety and the European NURESIM and NURISP projects (Rohde et al., 2009). A systematic series of benchmark exercises concerning neutron kinetics and reactor dynamics for VVER type reactors was considered within the AER (AER Benchmark Book, 1999). An overview on mathematical benchmarks solved for DYN3D verification is given in Table 3. 6.2. Code validation against measurement data DYN3D was validated against operational data and measurement data for transients in VVER reactors. Experiments in test reactors were used for code validation, too, particularly for VVER type reactors e.g. power distribution measurements for steady state as well as neutron kinetics measurements at the full size VVER-1000

181

experimental facility V-1000, which have been well documented within the EC project VALCO (Mittag et al., 2005; Hadek and Mittag, 2009; Hadek, 2012), and kinetic experiments at the LR-0 zero po  wer reactor in NRI Re z (Rypar and Hadek, 1988). Kinetic experiments performed at the LR-0 zero power reactor have also been used for the validation of the pin power reconstruction method in DYN3D (Hadek, 2012). Micro fission chambers for thermal neutron flux measurements were positioned in the center of selected fuel assemblies at different axial positions. Therefore, the detector readings were compared with reconstructed flux in the center of the assembly. Three approaches to determine the local flux at the central position were considered: (1) the node average value, (2) the value obtained by intra-nodal flux reconstruction and (3) the reconstructed value together with the pin power factors obtained from the lattice code (HELIOS) zero-net-current fuel assembly calculations. The mean deviation between calculated and measured values is less than 6%. However, single values show higher deviations (up to 19%). As expected, options (2) and (3) show a slightly better agreement with the measurement than option (1). Besides the measurements at experimental reactors, real plant data from operating NPP with VVER-1000 type reactors in the Ukraine were used for code validation. Calculational results were compared with measurement data from NPP with VVER-1000 type reactors, including critical boron concentration, reactivity coefficients and effectiveness of control rods. Full core calculations for the first three cycles of unit 6 of Zaporoshye NPP were performed by DYN3D using macroscopic XS libraries generated by use of the lattice code CASMO (Studsvik, 1995). Results of comparison between measured and calculated critical boron acid concentrations (cb) were compiled in Table 4. H10 is the position of the CR group No 10. Table 5 summarizes the results of comparison between measured and calculated reactivity coefficients for start-up experiments at the beginning of the 2nd fuel cycle at various positions of the control rod group No 10 (H10) (Kuchin et al., 2007). A benchmark for verification and validation of fuel assembly and whole core burnup calculations for VVER-1000 reactor core loadings with advanced FA (TVSA fuel assemblies containing Gd) was proposed including the calculation of four consecutive fuel cycles of the unit 2 of the Khmelnizki NPP (Loetsch et al., 2009). Task 1 of the benchmark comprises solutions and data for all fuel assembly types used in the core loadings. The second task consists of 3D core burnup calculations together with calculations of hot zero power critical states. For 3D burnup calculations, besides of other codes developed particularly for VVER type reactors, the DYN3D code was used with nuclear data libraries generated by the lattice codes CASMO and HELIOS. Significant differences between measured and calculated boron concentrations at the end of cycle occurred in the first calculations. The agreement was significantly improved by taking into account burnup spectral history (Ovdiienko et al., 2011). A comparison of calculated and measured (reconstructed) assembly-wise power distributions for the four cycles has been provided in Loetsch et al. 2013, showing that the codes and data libraries ensure a sufficient accuracy of calculation of FA power distributions. The mean square deviations between measured (reconstructed) and calculated values e for DYN3D with both CASMO and HELIOS generated macroscopic XS data libraries e are less than 2.2% with one exception. At 90 FPD of 3rd cycle, the maximum deviation was 3.66% with CASMO XS data resp. 3.17%

182

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

Table 3 List of benchmarks problems used for the validation of DYN3D. Identifier

Geometry

AER-FCM3D hex; 101 tri V1000-extra 2D hex; tri V1000-2D2D hex C1 V1000-2D2D hex C2 V1000-2D2D hex C3 V440-AER3D hex DYN2 V440-AER3D hex DYN3 V440-AER3D hex DYN4 NEACRP-CRE 3D rect NEANSC3D rect CRW OECD-MOX- 2D/3D Core rect Hebert 2D Hex; tri V1000-pin 2D Tri HTGR 2D Hex; tri

Ref.a

Type of calc.

Short description

Static, 2G diffusion

Core with partially inserted CR

Static & kinetics; 2G diffusion Dynamics 2G diffusion

Kol99; Roh09 Due13b VVER-1000 core sector against extrapolated fine mesh solution Cha95; Due13a; Due13c; VVER-1000 rodded and un-rodded core at HZP and CZP; fresh and depleted fuel; deterministic Pet06; Roh09 transport ref. solution VVER-1000 stuck rod benchmark; core with burnup; CZP or HZP; deterministic transport ref. Pet06; Roh09 solution VVER-1000 single channel boiling benchmark (reduced moderator density in one FA); Pet06; Roh09 deterministic transport ref. solution nd 2 dynamic benchmark of AER: Peripheral rod ejection in VVER-440 core at HZP; simplified AER99; Gru93; Doppler feedback only Gru94 3rd dynamic benchmark of AER: CR ejection AER99; Gru95

Dynamics; 2G diffusion

4th dynamic benchmark of AER: predefined slug of low borated coolant

AER99; Roh97

Dynamics; 2G diffusion Dynamics; 2G diffusion

NEA/CRP benchmark on CR ejection (various scenarios) NEA/NSC benchmark on CR withdrawal

Gru96a; Gru97 Gru96a; Gru97

Static & Kinetics; MG diffusion & SP3 Static; 2G Diffusion & SP3

OECD/NEA an U.S.NRC MOX Transient Benchmark

Gru06; Koz06

One-group benchmark with anisotropic scattering (academic)

Static; MG diffusion & SP3 Static; MG diffusion & SP3

Sector of one FA pinwise (VVER-1000) HTGR ring core & block with CR against Monte Carlo (Serpent code)

Heb10; Due12b; Due13c Due12a; Due13b; Due13c

Static; 2G Static, 2G diffusion Static, 2G diffusion Static 2G diffusion

a AER99¼ AER Benchmark Book, 1999; Cha95¼ Chao and Shatilla, 1995; Due12a, 12b¼ Duerigen et al., 2012a, b; Due13b¼Duerigen and Kliem, 2013; Due13c¼Duerigen, 2013; Gru93¼Grundmann and Rohde, 1993; Gru95¼Grundmann and Rohde, 1995; Gru96a¼Grundmann and Rohde, 1996a; Gru94¼Grundmann, 1994; Gru97¼ Grundmann bert, 2010; Kol99¼ Kolev et al., 1999; Koz06¼Kozlowski and Downar, 2006; Pet06¼ Petkov and Mittag, 2006, Roh97¼ and Rohde, 1997; Gru06¼Grundmann, 2006; Heb10¼He Rohde and Lucas, 1997.

Table 4 Comparison between measured and calculated boron concentrations for unit 6 of NPP Zaporoshye. Fuel cycle

Time, full power days (FPD)

Reactor power, MWth

Tinlet, ºC

H10, cm

Cb, g/kg measured

Cb, g/kg DYN3D/CASMO lib.

1

BOC, 18.8 EOC, 285.8 BOC, 18.3 EOC, 245.1 BOC, 44.9 EOC, 238.2

1236 2983 2940 2997 2991 2971

281 288 288 288 287 287

302 305 249 298 287 304

6.60 0.60 5.80 0.80 5.60 1.40

6.41 0.43 5.37 0.60 5.47 1.12

2 3

with HELIOS data. To simulate real transients in NPP, the coupled code complex DYN3D-ATHLET was used. A number of real plant transients were analyzed, in particular, for VVER type reactors (Kozmenkov et al., € m€ 2005, 2015; Kliem et al., 2006; Mittag et al., 2001; Ha al€ ainen et al., 2002).

6.3. Application examples 6.3.1. Operational and safety issues for VVER The DYN3D code has been developed for safety analyses of

reactor cores. Such analyses were performed by various nuclear reactor authorities or their technical support organizations, mostly for VVER type reactors. Naturally, not too many of these analyses were published because of their proprietary character. In a conservative approach, best-estimate codes (like DYN3D) with conservative boundary conditions were applied for the safety analysis of reactivity initiated accidents in the Ukraine using the code DYN3D. A transient scenario initiated by the drop of one cluster of the CR working group 10 at full power was calculated using DYN3D for unit 3 of the NPP South Ukraine. Due to the inadvertent drop of the

Table 5 Comparison of reactivity coefficients from DYN3D calculations and measurements. Parameter

DYN3D/CASMO

Position of CR group 10 Boron acid concentration, g/kg Boron acid concentration reactivity coefficient, acb*102,1/(g/kg) Total temperature reactivity coefficient (moderator and fuel), atotal*105,1/K Position of CR group 10 Boron acid concentration, g/kg Boron acid concentration reactivity coefficient, acb*102,1/(g/kg) Total temperature reactivity coefficient (moderator and fuel), atotal*105,1/K Efficiency of 10th group (H ¼ 80?0%), %

H10 ¼ 40% 8.76 1.39 8.78 H10 ¼ 80% 8.98 1.39 6.83 0.39

Measurement 8.74 1.43 7.2 9.03 1.43 4.34 0.32

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

cluster, the reactor power regulation system is activated, and the initial full reactor power is restored. The drop of the cluster leads to an increase of power in the upper part of assemblies situated far from the dropped cluster. It was found that, depending on the initial axial power profile, DNB may occur in the hottest rod of the corresponding assembly (Kuchin et al., 2010). Daily power maneuvering regime simulations that were proposed for 2ndunit of Khmelnitsky NPP were simulated using DYN3D. Calculation of peaking factors, linear pin power, control rod efficiency and minimal axial offset were performed for the 2nd cycle at 110 and 175 FPD. Calculation results were compared with experimental data, that were obtained during maneuvering regime's testing. The results of simulation show that a powermaneuvering regime using the 10th CR group together with the central control cluster results in a greater increase of peaking factors, compared to applying group 10 only, as initially supposed; but it requires the injection of a smaller amount of boron acid into the reactor core (Ovdiienko et al., 2006). In the following sections, three further examples are given, showing the variety of application areas of DYN3D. The first example is the control rod ejection in a VVER-1000 reactor core loaded with MOX fuel. The second one is a boron dilution analysis of a PWR core and the third one is the post-test calculation of a measured xenon dynamics transient. 6.3.2. Reactivity initiated accidents for MOX cores with weapongrade Pu A series of DYN3D calculations was performed to demonstrate the feasibility of weapon-grade plutonium burning in VVER-1000 type reactors with emphasis on safety aspects of the problem (Rohde et al., 2006). Particular objectives of weapon-grade Pu burning are the net burning of fissile Pu isotopes in general, but also the conversion of the Pu vector of the fuel from weapon-grade to reactor plutonium. While weapon-grade plutonium is almost pure fissile 239Pu, reactor plutonium is a mixture of the Pu isotopes 238, 239, 240, 241 and 242, which is not suitable for weapons. The even number isotopes are decaying spontaneously releasing neutrons, which would make a bomb un-controllable. Moreover, 238Pu decays with high heat release emitting alpha radiation. A big amount of accumulated weapon-grade plutonium makes the problem actual. A comprehensive study of Pu burning in nuclear reactors consists of the following major steps:





Fuel element design Core design Burnup calculations to get the equilibrium cycle Analyses of reactivity initiated accidents (RIAs)

Prior to RIA analyses, the first three of the above four steps were successfully carried out, including calculations aimed at proving that all design limits for the reactor core like maximum power peaking factors, maximum linear rod power or maximum coolant heat-up are met for normal operational conditions. It was demonstrated that the net burning of plutonium is achieved in the case of 100% MOX core-loading, while for the 50% MOX coreloading the weapon-grade Pu is converted to reactor plutonium. However, the focus for DYN3D is put mainly on the RIA analyses. Two scenarios of reactivity-initiated accidents (RIAs) were analyzed for three different core loadings with two fuel compositions e standard uranium dioxide and/or MOX fuel utilizing weapon-grade plutonium. Considered core-loading options are: (1) 100% loading with UO2 fuel, (2) 50% UO2 þ 50% MOX loading (shown in Fig. 10), and (3) 100% loading with MOX fuel. Fuel assemblies with profiled enrichment and gadolinium burnable absorber (6 or 12 rods with gadolinium per assembly) were used for

183

all options. The content of Gd2O3 in the UO2 fuel matrix is 5 w/o%. The uranium enrichment of rods with gadolinium is 3% for all coreloading options. The configuration of a typical profiled fuel assembly is presented in Fig. 10 B. The outer fuel pins have a lower uranium enrichment, respectively plutonium content compared with the internally located pins. Accidents were simulated for the beginning and for the specified time point (T ¼ 280 days) of the equilibrium fuel cycle (BOC and EOC). Altogether, 12 accident simulations were performed to compare the transient behavior of the VVER-1000 core loadings with and without MOX fuel. The analyzed RIA scenarios are: (a) uniformly accelerated ejection of the 6th control rods group (Fig. A) from the 284 cm position during 0.1 s and (b) unauthorized withdrawal of the same group from the same position with the speed of 2 cm/s. A failure of the reactor SCRAM was postulated in both scenarios. Case (a) was simulated up to 50 s, and case (b) - up to 100 s, starting from the accident initiation. Results of the analyses are summarized in Table 6. Evaluations show that the maximum fuel temperatures reached during the accidents are well below the fuel melting point, and the departure from nucleate boiling ratio (DNBR) is much greater than 1. Maximum reactor power is reached for the 2nd core-loading option (50% MOX fuel) at T ¼ 280 days. However, even in this case, the fuel rods were not overheated. For core-loading options 2 and 3, the highest fuel and cladding temperatures and the lowest DNBR value in MOX FA's are observed at the beginning of the cycle. Their values are kept within the permissible limits. Therefore, the nuclear safety against RIA was demonstrated for the MOX loadings, even for 100% MOX (core-loading option 3). In the simulations of scenario (a), the highest peak power at BOC was predicted for the core-loading option 3 due to the highest efficiency of the ejected control group in this case. At T ¼ 280 d (close to EOC) the peak power for option 3 is less than that for option 2 because of a lower worth of ejected control rods for option 3 due to increase of beff with increasing burnup in this case. In order to determine the maximum fuel temperature, the “hot channel” model of DYN3D was used with a power multiplication factor of 1.41. The highest fuel temperature and the lowest DNBR do not necessarily correspond to the case with the highest total peak reactor power, but with the highest local power density. Therefore, the maximum fuel temperature can be lower for EOC in comparison with BOC, even if the reactor peak power is higher, because power peaking factors at EOC are lower. In the simulations of scenario (b), only a very small power peak occurs due to the low speed of the control rod group withdrawal (2 cm/s). The inserted reactivity in this case is almost promptly compensated by negative reactivity feedback. However, as it was expected, the asymptotic (final) power levels are about the same as in scenario (a) due to the same core boundary conditions and the same initial and final positions of the control rod group. Maximum fuel temperatures and minimum critical power ratios are close to those obtained for scenario (a), but they are still a bit lower due to less prominent power excursions. 6.3.3. Boron dilution transient A boron dilution scenario has been analyzed based on the startup of the first main coolant pump (MCP) of a loop where a slug of de-borated water has accumulated. During the transport of the slug from the cold leg to the core inlet, mixing with the higher borated coolant in the reactor pressure vessel takes place. Time-dependent boron concentration curves at the inlet into each fuel assembly have been derived on the basis of experimental data from the

184

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

Fig. 10. a: Core loading pattern for option 2. b: Geometry of profiled fuel assembly.

ROCOM test facility (Prasser et al., 2003; Kliem et al., 2008). The reactor core under consideration is a generic PWR core at the beginning of an equilibrium fuel cycle. The reactor is in a hot sub-critical state with all control rods inserted except one stuck rod cluster in the sector with the highest deboration. The initial slug volume in the cold leg and the steam generator outlet chamber is 36 m3. Starting the MCP drives the slug towards the core. The boron concentration in the core decreases rapidly and the reactor becomes promptly super-critical (maximum reactivity of about 2 $) which leads to a very fast power excursion with peak power of about 8000 MW (Fig. 11). However, the power excursion is limited by the very effective Doppler feedback of the fuel temperature. The width of the power peak is only about 20 ms. Therefore, the integral energy release in the power peak is limited. This power peak occurs before the boron concentration has reached its minimum in the

core. As the positive reactivity insertion continues after the power peak appears, typical secondary power peaks emerge. The interaction of the continuing deboration with negative Doppler and void feedback effects determines the height and frequency of these power peaks. The total reactivity remains below the super-prompt margin, which is the reason why these power peaks are significantly lower. At the location of the power maximum, coolant boiling happens for a short time, with a maximum void fraction of about 70%. However, critical heat flux is not reached, so that cladding temperatures remain below 260  C and no safety-relevant limitations are violated (Kliem et al., 2004). This calculation was done using the above described PIC method for the boron transport, where practically no artificial dispersion of the deboration front occurs, as can be seen in the left part of Fig. 11, showing the time behavior of the boron concentration in a selected

Table 6 Estimated maximum values of the key safety parameters. Parameter

Initial state

Option 1

Option 2 MOX

Case (a): Ejection of control rod group Maximum reactor power (rel. to nominal) 

Max fuel temperature, C Max cladding temperature,  C DNBR Case (b): Unauthorized withdrawal of control rod group Maximum reactor power (rel. to nominal) Max fuel temperature,  C Max cladding temperature,  C DNBR

BOC 280 d BOC 280 d BOC 280 d BOC 280 d

1.23 1.55 1668 1594 350.9 350.6 1.19 1.20

1.26 1.79 1662 1488 351.0 350.6 1.21 1.28

BOC 280 d BOC 280 d BOC 280 d BOC 280 d

1.09 1.10 1667 1595 351.0 350.6 1.19 1.20

1.08 1.10 1658 1422 351.0 350.5 1.21 1.29

Option 3 UO2

1581 1610 350.8 350.8 1.29 1.16

1.33 1.58 1680 1514 351.0 350.7 1.18 1.23

1579 1612 350.9 350.7 1.28 1.16

1.06 1.07 1674 1513 351.0 350.6 1.18 1.23

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

185

Fig. 11. Behavior of boron concentration in a boron dilution transient calculated with a diffusive and non-diffusive numerical scheme and corresponding behavior of the reactor power.

fuel assembly of the considered core. Applying the usual finite difference scheme, some dispersion during the transport through the core due to numerical diffusion is observed. With this method, the first power peak is slightly higher because the dilution starts earlier, but there are no secondary power peaks as they were obtained using the PIC method. Although these secondary power peaks are smaller by magnitude, they are broader and the integral of the released energy is higher. Therefore, use of diffusive numerical schemes can lead to non-conservative results with respect to reactor safety.

moved in, starting from the insertion depth of 15%, ending at 42%, with the whole insertion procedure taking about an hour. After keeping the control rods at this lowered position for 2.5 h, they were removed to their initial depth (15%) in the same slow way as they had been inserted. During this operation, the reactor power was stabilized by varying the boron acid concentration in the coolant. The xenon dynamics triggered by this control rod motion were modeled by DYN3D and compared to the measurement. Fig. 12 depicts the axial offset (AO) of the reactor power:

AO ¼ ½ðpower in upper core halfÞ  ðpower in lower halfÞ 6.3.4. Xenon dynamics DYN3D has been designed primarily for studying possible fasttransient accident scenarios in reactor cores. A careful simulation of such transients requires the input of accurate cross section data, valid for the moment at the start of the considered event. These data depend on local burnup and operation parameters, but they are also influenced by reactor poisons, such as 135Xe and 149Sm and their (slow) dynamics. In normal smooth reactor operation at nominal power, these nuclides are present in equilibrium concentrations, determined by the spatial neutron-flux distribution. At the start of a reactor cycle, however, and generally after power (load) changes, non-equilibrium spatial distributions of these strong neutron absorbers can be observed. DYN3D is capable of modeling this poison dynamics thus producing, in a separate calculation, the actual input (cross section) information to be used in the following simulation of a fast transient. DYN3D calculations were performed for cycle No 14 of unit 3 of NPP Rovno. Xenon oscillations were induced by movement of CR group No 10, accompanied by a change of boron concentration. The oscillations are stable at low power level at BOC conditions. However, Xe oscillations may become unstable at higher burnup with increase of reactor power (at power >2200 MW at EOC). An algorithm for the damping of xenon oscillations after fast changes of reactor power close to EOC was developed and tested on DYN3D simulations (Khalimonchuk, 2008). A measurement studying poison-related reactor dynamics was carried out in the VVER-1000 of the NPP Chmelnizki-2, Ukraine (Rohde et al., 2004). Control rod bank number 10 was slowly inserted at the thermal power of about 2200 MW. The bank was

 =ðwhole core powerÞ*100%: Measured AO values are derived from the VVER-1000 in-core power-measuring system. An axial xenon oscillation is observed fading away without any plant operator action. This self-stabilizing behavior is mainly due to the thermal-hydraulic feedback (cf. Sections 3.4 And 3.5) in the VVER-1000 (PWR). Additional DYN3D test calculations, which artificially neglect this feedback, demonstrate that the core would be unstable against axial xenon oscillations without temperature feedback. Altogether, DYN3D is able to describe the measured poison dynamics; the deviation seen in Fig. 12 (period of the damped oscillation) might be connected to the fact that not the whole operation history during the days prior to starting this measurement could be taken into account in the simulation. 7. Developments and applications for innovative reactor concepts Besides of standard version of DYN3D for Light Water Reactors, special code versions for Generation IV type reactors have been developed or are under development. The considered reactor types include Molten Salt Reactors (MSR), block-type gas cooled High Temperature Reactors (HTR), and Sodium cooled Fast Reactors (SFR). 7.1. MSR developments The code version DYN3D-MSR was developed for Molten Salt

186

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

1970). The molten salt mixture circulated through a reactor vessel at temperatures above 600  C. The fuel consisted of lithium, beryllium, and zirconium fluoride salts, to which 235U was added as UF4. Later on, a small amount of 239Pu fluoride was added to the mixture. Furthermore, before the shut-down, the core was operated also with 233U.
The Molten Salt Breeder Reactor (MSBR) (Rosenthal et al., 1971) was designed as an industrial energetic successor of MSRE, but it was never realized. MSRE parameters and results of experiments are well described in ORNL reports. Therefore, it was decided by the European project MOST (Review of MOlten Salt reactor Technology) participants to define a common benchmark (Delpech et al., 2003) based on these data. Two benchmark exercises focused on the protected fuel pump start-up and coast-down transients. During these two transients a constant power was maintained. The change in the delayed neutron production was compensated by the use of control rods and the inserted reactivity was the variable to be compared with experiment. Moreover, a natural circulation experiment performed at the MSRE was calculated. The agreement of experimental data and numerical simulation is very good (Krepel et al., 2007). A number of hypothetical transients were simulated by DYN3DMSR both for the MSRE and the MSBR (Krepel et al., 2006): Fig. 12. Axial offset during a damped Xenon oscillation.

Reactor (MSR) analyses (Krepel et al., 2007, 2008). In the MSR system, the fuel is a circulating liquid mixture of sodium, zirconium, and uranium fluorides (MacPherson, 1985). A crucial point of MSR dynamics is that the fuel is dissolved in the liquid salt. The salt serves as fuel and heat carrier at the same time. Therefore, the precursor nuclides of delayed fission neutrons are moving with the salt, which leads to the effect that the delayed neutrons are released in a place different from the location where the precursors have been created by fission. This affects the reactivity contribution of the delayed neutrons, because the effective delayed neutrons fraction depends on the flow velocity of the salt. A corresponding convective transport term has to be added to the delayed neutrons balance equations in DYN3D. In MSR with thermal spectrum, the molten salt fuel flows through central channels in the hexagonal graphite blocks. For this design, the parallel flow channel concept as it is implemented in the DYN3D thermal-hydraulic model can be used. The heat exchange between graphite and salt represents a part of the heat source in the liquid salt; therefore, this source had to be modified in the thermal-hydraulic module of DYN3D for the MSR purpose. About 10% of the thermal power is released in the graphite by gamma heating. Furthermore, the thermodynamic water properties were replaced with the salt properties. Molten Salt Reactors with fast neutron spectrum (MSFR) are of particular interest for actinides transmutation. For a fast spectrum, graphite and other moderating materials must not be present in the core. The presence of structural materials inside the core has to be avoided. Therefore, the present version of DYN3D-MSR, which is limited to two neutron energy groups and modeling of parallel, channel type flow in the core, is not applicable to MSFR. For the numerical simulation, two MSR reactors studied in Oak Ridge National Laboratory (ORNL) in the 1960's and 1970's were selected:
The Molten Salt Reactor Experiment (MSRE) was an 8 MWth fluid-fueled test reactor that operated from 1965 through 1969 in support of a program to develop a thermal breeder power reactor using the thorium-uranium fuel cycle (Haubenreich,


reactivity promptly jumps up to 300 pcm for MSRE and MSBR,
unprotected fuel pump coast-down for MSRE and MSBR at nominal and zero power,
over-fueled slug flow transients initiated at zero power by fuel pump start-up (assuming an inadvertently increased fuel concentration in a slug of liquid salt entering the core) for MSRE and MSBR,
over-cooled slug flow transients for MSRE and MSBR, and
local fuel channel blockage in the MSBR core. Fig. 13 shows the behavior of reactor power after an unprotected fuel pump transient for the MSRE. The fuel pump coast-down causes gain in the effective delayed neutron precursors (DNP) concentration, which induces positive reactivity. On the other hand, the pump coast down leads to the loss of heat sink in the nominal power case. It causes a fast temperature increase and the induced reactivity almost compensates the DNP gain. There is no such a fast compensation in the zero power case and the growth caused by DN gain is stopped later, when a significant power level is reached. However, the magnitude of the power peek depends also on the neutronic parameters, especially on the average neutron generation time L and on the effective delayed neutrons fraction beff. Therefore, the power grows faster resp. decreases slower in the case of 233U fuel loading than in the 235U case. 7.2. HTR developments A code version DYN3D-HTR for application to block-type HTR has been developed (Rohde et al., 2012). One of the unique features of HTRs is associated with the use of highly robust tristructural isotropic (TRISO) coated fuel. In blocktype HTRs, the TRISO particles are compacted into fuel elements (fuel compacts), which are placed in graphite blocks with cooling channels. This fuel configuration has a double-heterogeneous structure which is a challenge for computational tools that were originally developed for LWR applications. Therefore, new procedures were developed for the homogenization of the doubleheterogeneous fuel element structures. On the one hand, the socalled Reactivity equivalent Physical Transformation (RPT), the two-step homogenization method based on 2D deterministic

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

10000

Power (kW)

time transients. Furthermore it was shown that the 3D heat conduction model is necessary to describe specific effects as warm-up of reflector assemblies accurately.

U233, Nom. Power U235, Nom. Power U233, Zero Power U235, Zero Power

8000

187

7.3. SFR developments

6000 4000 2000 0 0

2

4 6 Time (min)

8

10

Fig. 13. MSRE response to the unprotected fuel-pump coast-down at zero and nominal power for two different core loadings. Taken from Krepel et al. 2006.

lattice calculations was extended to cells with different temperatures of the materials. On the other hand, the progress in development of Monte Carlo methods for spectral calculations, in particular the development of the code SERPENT, opened a new fully consistent 3D approach, where all details of the structures at fuel particle, fuel compact and fuel block level can be taken into account within one step. A 3D heat conduction and heat transfer model including radial heat conduction between fuel assemblies was developed by taking into account the multi-level structure of the system. For every node, the model considers one representative temperature for the TRISO fuel kernel, all TRISO coatings, fuel compacts, and the graphite matrix. The coupling between TRISO and the fuel compact as well as between the fuel compact and the graphite matrix is realized by effective heat transfer coefficients taken from analytical steady state heat conduction solutions. The heat conduction model is coupled to the one-dimensional coolant channel flow model of DYN3D, developed earlier for Light Water Reactors. It is also suitable for prismatic HTRs, due to the lack of cross flow inside the reactor core. The application of the multi-group trigonal nodal solvers (Section 3) allows for more accurate modeling of the HTR block with asymmetrically located control rods and coolant channels. DYN3D-HTR was applied to a number of steady state and transient calculations. The latter comprise control rod withdrawal, helium overcooling and gas flow reduction. Fig. 14 shows the temporal evolution of the total power and the averaged graphite, compact graphite, fuel and helium temperatures for a reactivityinitiated transient (complete control rod ejection between 100s and 200s e left hand panel) and for an increase of the inlet coolant temperature of 100 (right hand panel) (Baier et al., 2014). During control rod movement, the power dramatically increases which leads to an increase in fuel temperatures. Compact, graphite and coolant temperatures follow this evolution with some time delay. Due to temperature feedback effects, the power and temperature decrease again. After a few damped oscillations the system relaxes to a new steady state with increased power and increased temperatures. The growth of the coolant temperature (right hand panel) leads to a decrease of removed heat and, thus, to an increase of the graphite temperature followed by increase of the compact and fuel temperature. Again due to the temperature feedback the power decreases and eventually a new steady state with decreased power and increased temperatures arises. The results indicate that DYN3D-HTR is able to simulate short-

Currently, the DYN3D code is being extended to the analysis of SFR cores. Thermal-physical properties of sodium (such as thermal conductivity, density, heat capacity and viscosity) were included into the thermalehydraulics module database. A two-step reactor physics analysis methodology devoted to the modeling of SFR cores was developed. According to this two-step procedure, DYN3D is used for full core calculations while the continuous-energy Monte Carlo code Serpent is considered as a tool for the generation of fewgroup XS required by DYN3D. The few-group XS generation methodology was established and verified on various simplified 2D and more realistic 3D SFR core configurations (Fridman and Shwageraus, 2013 and Rachamin et al., 2013; Nikitin et al., 2015a, 2015b). The results of the verification studies demonstrated the applicability of DYN3D to steady-state analysis of SFR cores provided that a high quality XS section library is available. In contrast to LWRs, dimensional variations of SFR core components caused by thermal expansions provide an important reactivity feedback under normal, transient and accident conditions. In order to account for thermal expansion effects (currently axial fuel rod expansion and radial diagrid expansion), a thermal-mechanical model is being developed and implemented in DYN3D. 8. Management of code maintenance The DYN3D project started in the 1980ies by U. Grundmann and U. Rohde. The amount of source code as well as the number of involved developers and supported compilers significantly increased over time. Currently the code consists of 550 files with about 85,000 lines of code. To assure the quality of the DYN3D developments, a unique system of tools for code management, development, maintenance and validation was established. One part concerns the Code Conventions which have to be kept by each active developer. They define general coding rules for the DYN3D development resulting in a more consistent, readable and maintainable source code. Another part of code maintenance is the version control system based on Subversion (SVN) from Apache. It provides a central and safe place for the DYN3D sources and the access control to it. It involves the DYN3D repository including the source code, related documents, a test matrix and a build system. By interacting with the SVN server the developers check-out local working copies, commit local changes and synchronize to the repository. This version control prevents code fragmentation and inconsistency. To support efficient DYN3D development there is a build system which covers multi-platform and multi compiler issues for code compiling and testing. It is based on CMake from Kitware and is for Intel Fortran (ifort) and GNU Fortran (gfortran) under Windows and Linux. The build system is used by the developers to validate their local developments before committing them to the SVN server. The repository includes a test matrix section to assure (permanent) code validation while the development continues. It keeps a set of test cases at different levels of validation: Real test cases (e.g. from benchmarks) as well as academic test cases to check DYN3D models and input options. Each case defines its own reference results to be compared with the development state. The automated execution and validation of the test cases against new development is managed by the build system. An analysis tool has been developed for how the DYN3D code is covered by the test matrices in

188

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

Fig. 14. Temporal evolution of temperatures and power for a control rod ejection (left hand panel) and for an overcooling transient (right hand panel).

order to identify tested and untested code parts. The results of the application of this tool serve as a basis of the extension of the test matrices; so that the matrices cover in their final states all features and models of the code. For independent code checking, a continuous integration system (CI) based on Jenkins has been implemented. It is triggered at each commit time by the SVN server. By this it invokes distributed tasks for code validation at different OS and compiler levels. Then the CI collects the task results and provides logs and statistics. Each developer can access the CI system through a web interface and also review the results. Thus, the combination of those tools guarantees a high quality assurance of the DYN3D code developments. All code management and maintenance tools are under continuous development. The availability and every-day-use are the precondition for the start of a DYN3D code refactoring program which is foreseen for the near future.

9. Summary and outlook Summarizing the status of development the core simulator code DYN3D has proven to be an effective and accurate tool for steadystate and transient core calculations thanks to the application of advanced nodal methods, the multi-group approach, consideration of discontinuity factors and intra-nodal flux reconstruction in the neutronic part of the code. Furthermore, an integrated fourequations thermal-hydraulic model for the description of oneand two-phase flow in the reactor core, and a thermo-mechanical model for the fuel behavior for fresh and low-burnup fuel which have been qualified for the simulation of reactivity initiated accidents are available and validated. For the simulation of complex transient scenarios with interaction between reactor core and full plant behavior, the coupling of DYN3D with the thermal-hydraulic system codes ATHLET and RELAP5 as well as with the advanced fuel performance code TRANSURANUS can be used as an industrial standard. Multi-physics couplings with high resolution fluid dynamics codes are being implemented as well. Furthermore, the application of DYN3D for transient analyses of advanced, innovative reactor concepts is under development. The current and foreseen developments of DYN3D are related to different code sections. Recently the development and implementation of a new method for accounting spectral history effects

hast been started. Detailed nuclide content is calculated for each region of the reactor core. Thus, a more precise modeling of fuel composition during depletion, history effects, poisonings, and decay heat release will be enabled. Furthermore, it is foreseen to continue the development of a transport solver which should be embedded into the global diffusion solution and should be able to provide a more accurate pin power distribution in selected core regions. The extension of the application area of DYN3D to fast reactor calculations will be continued. The main work is related to the development and implementation of the thermal-mechanical model for the consideration of the reactivity effects caused by thermal expansion of the core components. Additionally it is foreseen to apply the coupling of DYN3D with the system code ATHLET for fast reactor simulations. The rich HZDR experience in coupling of codes through the NURESIM SALOME platform is being used to develop a SALOME coupling interface with the fuel performance code TRANSURANUS. This work is conducted in close co-operation with the ITU Karlsruhe. As many reactor physics codes DYN3D has a long history with origins in the 1970s. There are still many parts of the codes which are more than 20 years old. However, these code fragments do not fulfill the requirements of a modern modular/flexible code structure. Thus, refactoring, restructuring and modularization of the DYN3D code structure were identified as urgent issues and will be started in the near future to facilitate the implementation of new features. Acknowledgments dek, Petko Petkov, Yuri Ovdiienko Thanks to Jiri Krepel, Jan Ha and Maksim Ieremenko for their contributions to DYN3D development and validation as well as to Dirk Lucas and Armando Gomez for their efforts on multi-physics coupling of DYN3D with other codes. Thanks to the German Federal Ministry of Economics and Technology and to the European Commission for co-funding of development and validation projects. References €ger, P., 1987. Makroskopische Kühlmittelvermischung in Ackermann, G., Dra Druckwasserreaktoren. Kernenergie 30, 454e458.

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190 AER Benchmark Book, 1999. http://aerbench.kfki.hu/aerbench/AEKI-KFKI,Hungary. Aleshin, S.S., 2005. Preparation of Licensing Documents for New Version of PERMAK-a Code Designed for Calculation of Cores with New Design Fuel Assemblies (Report of RRC Kurchatov Institute, Moscow, Russia). Austregesilo, H., Bals, C., Hora, A., Lerchl, G., Romstedt, P., 2006. ATHLET, Mod 2.1 Cycle A, Models and Methods, GRS - P - 1, vol. 4. Gesellschaft für Anlagen- und Reaktorsicherheit (GRS) mbH, Germany. Baier, S., Fridman, E., Kliem, S., Rohde, U., 2014. Extension and application of the reactor dynamics code DYN3D for block-type high temperature reactors. Nucl. Eng. Des. 271, 431e436. Beckert, C., Grundmann, U., Mittag, S., 2007. Multigroup diffusion and SP3 solutions for a PWR MOX/UO2 benchmark with the code DYN3D. Trans. Am. Nucl. Soc. Eur. Nucl. Soc. 97, 701e702. Beckert, C., Grundmann, U., 2008. Development and verification of a nodal approach for solving the multigroup SP3 equations. Ann. Nucl. Energy 35, 75e86. Bilodid, I., Mittag, S., 2010. Use of the local Pu-239 concentration as an indicator of burnup spectral history in DYN3D. Ann. Nucl. Energy 37, 1208. Bilodid, Y., Ovdiienko, I., Mittag, S., Kuchin, A., Khalimonchuk, V., Ieremenko, M., 2012. Assessment of spectral history influence on PWR and WWER core. Kerntechnik 7, 4, KT11025. Bilodid, Y., Fridman, E., Kotlyar, D., Margulis, M., Shwageraus, E., 2014. Verification of the spectral history correction method with fully coupled Monte Carlo code BGCORE. In: Proc. Of PHYSOR 2014, 28.09.-03.10.2014, Kyoto, Japan. Brantley, P.S., Larsen, E.W., 2000. The simplified P3 approximation. Nucl. Sci. Eng. 134, 1e21. Casal, R.J.J., Stamm’ler, E. A. Villarino, Ferri, A.A., 1991. HELIOS: geometric capabilities of a new fuel assembly program. In: Proc. Of the Int. Topical Meeting on Advances in Mathematics. Computations and Reactor Physics, vol. 2. Pittsburgh, Pennsylvania, USA, April 28-May 2, 1991. Chao, Y.-A., Shatilla, Y.A., 1995. Conformal mapping and hexagonal nodal methods e II: Implementation in the ANC-H code. Nucl. Sci. Eng. 121, 210e225. Delpech, M., Dulla, S., Garzenne, C., Kophazi, J., Krepel, J., Lebrun, C., Lecarpentier, D., Mattioda, F., Ravetto, P., Rineiski, A., Schikorr, M., Szieberth, M., 2003. Benchmark of dynamics simulation tool for molten salt reactors. In: Proc. Int. Conf. GLOBAL 2003, 2182e2187, New Orleans, LA, November 2003. DIN Norm 25463, 1990. Part 1, “Berechnung der Nachzerfallsleistung der Kernbrennstoffe von Leichtwasserreaktoren”. Duerigen, S., Fridman, E., 2012. The simplified P3 approach on a trigonal geometry of the nodal reactor code DYN3D. Kerntechnik 4, 226e229. Duerigen, S., Bilodid, Y., Fridman, E., Mittag, S., 2012a. A trigonal nodal SP3 method with mesh refinement capabilities - development and verification. In: PHYSOR 2012 Advances in Reactor Physics, 15.-20.04.2012, Knoxville, Tennessee, USA, on CD-ROM, American Nuclear Society, LaGrange Park, IL (2012). Duerigen, S., Nikitin, E., Mittag, S., 2012b. Verification of the trigonal-geometry diffusion and SP3 models of the code DYN3D. In: Proc. 22nd Symposium of Atomic Energy Research (AER), 01.-05.10.2012, Pruhonice, Czech Republic. Duerigen, S., Rohde, U., Bilodid, Y., Mittag, S., 2013. The reactor dynamics code DYN3D and its trigonal-geometry nodal diffusion model. Kerntechnik 310e318, 78/4, 2013. Duerigen, S., Kliem, S., 2013. The DYN3D trigonal-geometry diffusion model e verification using the AER-FCM-101 benchmark. In: Proc. 23rd Symposium of   pleso, Slovakia. Atomic Energy Research (AER), 30.09.-04.10.2013, Strbsk e Duerigen, S., 2013. Neutron Transport in Hexagonal Reactor Cores Modeled by Trigonal-geometry Diffusion and Simplified P3 Nodal Methods. Dissertation. Karlsruhe Institute of Technology, p. 175, 2013. Finnemann, H., Bennewitz, F., Wagner, M.R., 1977. Interface current techniques for multidimensional reactor calculations. Atomkernenergie 30, 123e128. Finnemann, H., 1979. Nodal expansion method for the analysis of space-time-effects in LWR's. In: Proceedings of a NEA Specialists Meeting on ‘The Calculation of 3dimensional Rating Distributions in Operating Reactors’, Paris 1979, pp. 257e281. Fridman, E., Duerigen, S., Bilodid, Y., Kotlyar, D., Shwageraus, E., 2013. Axial discontinuity factors for the nodal diffusion analysis of high conversion BWR cores. Ann. Nucl. Energy 62, 129e136. Fridman, E., Shwageraus, E., 2013. Modeling of SFR cores with SerpenteDYN3D codes sequence. Ann. Nucl. Energy 53, 354. Gelbard, E.M., 1960. Application of Spherical Harmonics Method to Reactor Problems. Tech. Report WAPD-BT-20. Bettis Atomic Power Laboratory. Gommlich, A., Kliem, S., Rohde, U., Gomez, A., Sanchez, V., 2010. Coupling of the neutron-kinetic core model DYN3D with thermal hydraulic code FLICA-4 within the NURESIM platform. In: Annual Meeting of the German Nuclear Society, Berlin, May 3-6, 2010. Gomez-Torres, A.M., Sanchez-Espinoza, V.H., Ivanov, K., Macian-Juan, R., 2012. DYNSUB: a high fidelity coupled code system for the evaluation of local safety parameters e Part I: development, implementation and verification. Ann. Nucl. Energy 48, 108. Gomez, A., Sanchez Espinosa, V.H., Kliem, S., Gommlich, A., 2014. Implementation of a fast running full core pin power reconstruction method in DYN3D. Nucl. Eng. Des. 274, 86e97. Grahn, A., Kliem, S., Rohde, U., 2015. Coupling of the 3D neutron kinetic core model DYN3D with the CFD software ANSYS CFX. Ann. Nucl. Energ 84, 197e203. Grundmann, U., Rohde, U., 1993. Definition of the second kinetic benchmark of AER. In: Proc. 3rd Symposium of Atomic Energy Research (AER). KFKI Atomic Research Institute, Budapest.

189

Grundmann, U., 1994. Results of the second kinetic AER-benchmark. In: Proc. 4th Symposium of Atomic Energy Research (AER). KFKI Atomic Research Institute, Budapest. Grundmann, U., Lucas, D., Rohde, U., 1995. Coupling of the thermo-hydraulic code ATHLET with the neutron kinetic core model DYN3D. In: Proc. Int. Conf. On Mathematics and Computations, Reactor Physics, and Environmental Analyses, Portland, Oregon, USA, April 30-May 4 1995. Grundmann, U., Rohde, U., 1995. Comparative analysis of the third threedimensional dynamic AER benchmark problem with the help of the code €, Dyn3D. In: Proc. Fifth Symposium of Atomic Energy Research (AER), Dobogoko Hungary, 15-19 October, 1995, Proc, pp. 329e343. Grundmann, U., Rohde, U., 1996a. The reactor code DYN3DR - transient calculations of NEACRP benchmarks for PWR and BWR. In: Annual Meeting on Nuclear Technology, Mannheim (Germany), May 21 e 23 1996, Proc, pp. 23e26. Grundmann, U., Rohde, U., 1996b. DYN3D - A 3-dimensional core model for steady state and transient analysis in Thermal reactors. In: Proc. Int. Conf. On the Physics of Reactors “PHYSOR 96”, Mito (Japan), 16e20 Sept. 1996. Grundmann, U., Rohde, U., 1997. Verification of the Code DYN3D/R with the Help of International Benchmarks. Forschungszentrum Rossendorf, Technical Report FZR-195, October 1997. Grundmann, U., Hollstein, F., 1999. A two-dimensional intranodal flux expansion method for hexagonal geometry. Nucl. Sci. Eng. 133, 201e212. Grundmann, U., 1999. HEXNEM e a nodal method for the solution of the neutron diffusion equation in hexagonal geometry. In: Proc. Of the Int. Conf. on Mathematics and Computation, Reactor Physics and Environmental Analysis in Nuclear Applications (M&C ’99), Madrid, Spain, September 27-30, 1999. Grundmann, U., Rohde, U., Mittag, S., 2000. DYN3D - three dimensional core model for steady-state and transient analysis of Thermal reactors. In: PHYSOR 2000Advances in Reactor Physics and Mathematics and Computation into the Next Millennium, Pittsburgh (USA), May 7-12, 2000, Proc. CD-ROM, ANS Order No. 2700281. Grundmann, U., Rohde, U., Mittag, S., Kliem, S., 2005. DYN3D Version 3.2 e Code for Calculation of Transients Light Water Reactors (LWR) with Hexagonal or Quadratic Fuel Elementse Description of Models and Methods. Tech. Report FZR-434, Forschungszentrum Rossendorf ([FZR), August 2005. Grundmann, U., 2006. Calculations of a steady state of the OECD/NRC PWR MOX/ UO2 transient benchmark with DYN3D. In: Annual Meeting on Nuclear Technology 2006, 16.-18.05.2006. CD-ROM, Aachen, Germany. Grundmann, U., Mittag, S., 2011. Super-homogenisation factors in pinwise calculations by the reactor dynamics code DYN3D. Ann. Nucl. Energy 38, 2111e2119. dek, J., Grundmann, U., 1997. Neutron flux reconstruction in a hexagonal cassette Ha e theory and implementation into the code DYN3D/H1.1. Nucleon 3, 8. dek, J., Mittag, S., 2009. Validation of DYN3D pin-power calculation against Ha experimental VVER-full-core benchmark. In: Int. Conf. On Mathematics, Computational Methods & Reactor Physics (M&C 2009), 03.-07.05.2009. Saratoga Springs, New York, USA. dek, Jan, 2012. Validation of DYN3D pin-power calculation against experimental Ha results from the LR-0 reactor. In: Proc. 22nd Symposium of Atomic Energy Research (AER), Pruhonice, Czech Republic, October 1-5, p. 735. €ma €l€ Ha ainen, A., Kyrki-Rajam€ aki, R., Mittag, S., Kliem, S., Weiss, F.-P., Langenbuch, S., Danilin, S., Hadek, J., Hegyi, G., 2002. Validation of coupled neutron kinetic/ thermal-hydraulic codes Part 2: analysis of a VVER-440 transient (Loviisa-1). Ann. Nucl. Energy 29, 255e269. Haubenreich, P.N., 1970. Molten-salt Reactor Experiments. ORNL-4396. Hebert, A., Kavenoky, A., 1981. Development of the SPH homogenization method. In: Proc. Int. Topical Meeting on Advances in Mathematical Methods for the Solution of Nuclear Engineering Problems, Munich, Germany, April 27-29 1981. bert, A., 2010. Mixed-dual implementations of the simplified PN method. Ann. He Nucl. Energy 37, 498e511. n-Juan, R., 2014. Holt, L., Rohde, U., Seidl, M., Schubert, A., Van Uffelen, P., Macia Two-way coupling between the reactor dynamics code DYN3D and the fuel performance code TRANSURANUS at assembly level. In: 22nd International Conference on Nuclear Engineering (ICONE22). ASME, Prague, Czech Republic. Holt, L., Rohde, U., Seidl, M., Schubert, A., van Uffelen, P., 2015. Development of a general coupling interface for the fuel performance code TRANSURANUS tested with the reactor dynamics code DYN3D. Ann. Nucl. Energy 84, 73e85. Holt, L., Rohde, U., Kliem, S., Baier, S., Seidl, M., Van Uffelen, P., Maci an-Juan, R., 2016. Investigation of feedback on neutron kinetics and thermal hydraulics from detailed online fuel behavior modeling during a boron dilution transient in a PWR with the two-way coupled code system DYN3D-TRANSURANUS. Nucl. Eng. Des. 297, 32e43. IAEA, 1996. Selected Safety Aspects of WWER-440 Model 213 Nuclear Power Plants, Final Report of the IAEA Project on Evaluation of Safety Aspects for WWER-440 Model 213 Nuclear Power Plants. Int. Atomic Energy Agency, Vienna, ISBN 92-0101196-2. Ieremenko, M., Ovdiienko, I., Kuchin, A., Khalimonchuk, V., Wainer, L., 2012. Validation of pin power calculations by DYN3D on MIDICORE and AER KHNPP-2 benchmarks. In: Proc. Proc. 22nd Symposium of Atomic Energy Research (AER), Pruhonice, Czech Republic, October 1-5, 2012, p. 107. Imke, U., Sanchez, V., Gomez, R., 2010. Subchanflow: a new empirical knowledge based subchannel code. In: Proc KTG Annual Meeting on Nuclear Technology 2010, Berlin, Germany, May 4e6, 2010. Jimenez, G., Herrero, J.J., Gommlich, A., Kliem, S., Cuervo, D., Jimenez, J., 2015. 2015. Boron dilution transient simulation analyses in a PWR with neutronics/thermal hydraulics coupled codes in the NURISP project. Ann. Nucl. Energy 84, 86e97.

190

U. Rohde et al. / Progress in Nuclear Energy 89 (2016) 170e190

Khalimonchuk, V.A., 2008. Dinamika jadernoho reaktora s raspredeljonnymi parametrami v issledovaniakh perekhodnykh reshimov ekspluatazii VVEP i RBMK (Dynamics of nuclear reactors with distributed parameters in investigations of transient processes in VVER and RBMK). in Russian. Osnova, Kiev, ISBN 978966-699-433-5. €hne, T., Weiß, F.-P., 1999. Main steam line break analysis of a VVER-440 Kliem, S., Ho reactor using the coupled thermo-hydraulics system/3D-neutronkinetics code DYN3D/ATHLET in combination with the CFD code CFX-4. In: 9th Int. Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-9), San Francisco, California, Oct. 3-8, 1999. ISBN -89448-650-0. Kliem, S., Rohde, U., Weiß, F.-P., 2004. Core response of a PWR to a slug of underborated water. Nucl. Eng. Des. 230, 121e132. Kliem, S., Kozmenkov, Y., Hoehne, T., Rohde, U., 2006. Analyses of the V1000CT-1 benchmark with the DYN3D/ATHLET and DYN3D/RELAP coupled code systems including a coolant mixing model validated against CFD calculations. Prog. Nucl. Energy 48, 830e848. Kliem, S., Rohde, U., 2007. DYN3D/ATHLET calculations of a boron dilution transient during natural circulation conditions. In: Proc. Of the 12th Int. Top. Meeting on Nuclear Reactor Thermal-Hydraulics. CDROM paper 053. Kliem, S., Suehnel, T., Rohde, U., Hoehne, T., Prasser, H.-M., Weiss, F.-P., 2008. Experiments at the mixing test facility ROCOM for benchmarking of CFD-codes. Nucl. Eng. Des. 238, 566e576. Kliem, S., Mittag, S., Rohde, U., Weiß, F.-P., 2009. ATWS analysis for PWR using the coupled code system DYN3D/ATHLET. Ann. Nucl. Energy 36, 1230. Kliem, S., 2010. Ein Modell zur Beschreibung der Kühlmittelvermischung und seine Anwendung auf die Analyse von Borverdünnungstransienten in Druckwasserreaktoren. PhD thesis. Technical University Dresden, Germany. Kliem, S., Gommlich, A., Grahn, A., Rohde, U., Schuetze, J., Frank, T., Gomez, A., Sanchez, V., 2011. Development of multi-physics code systems based on the reactor dynamics code DYN3D. Kerntechnik 76, 160. Kochunas, B., Stimpson, S., Collins, B., Downar, T., Brewster, R., Baglietto, E., Yan, J., 2012. Coupled full core neutron transport/CFD simulations of pressurized water reactors. In: Proc. PHYSOR 2012, Knoxville, Tennessee, USA. Koebke, K., Hetzelt, L., Wagner, M.R., Winter, H.-J., 1984. Principles and application of advanced nodal reactor analysis methods. In: Proceedings of the Topical Meeting on Reactor Physics and Shielding, Chicago, Illinois, USA, September 17-19, 1984. Kolev, N.P., Lenain, R., Fedon-Magnaud, C., 1999. AER-FCM-101Benchmark SpecificationSheet. In: AER Benchmark Book. http://aerbench.kfki.hu/aerbench/. Kozlowski, T., Downar, T.J., 2006. OECD/NEA and US NRC PWR MOX/UO2 Core Transient Benchmark e Final Report, p. 20. NEA/NSC/DOC. Kozmenkov, Y., Orekhov, Y., Grundmann, U., Kliem, S., Rohde, U., Seidel, A., 2001. Development and benchmarking of the DYN3D/RELAP code system. In: Annual Meeting of the German Nuclear Society, Dresden, Germany, May 15-17 2001. Kozmenkov, Y., Grundmann, U., Kliem, S., Rohde, U., 2005. Comparative assessment of coupled RELAP5/PARCS and DYN3D/RELAP5 codes against the Kozloduy-6 pump trip test. In: Proc. Of the Int. Top. Meeting on Mathematics and Computation, Supercomputing, Reactor Physics and Nuclear Applications. CDROM paper 170. Kozmenkov, Y., Kliem, S., Grundmann, U., Rohde, U., Weiss, F.-P., 2007. 2007. Calculation of the VVER-1000 coolant transient benchmark using the coupled code systems DYN3D/RELAP5 and DYN3D/ATHLET. Nucl. Eng. Des. 237, 1938e1951. Kozmenkov, Y., Kliem, S., Rohde, U., 2015. Validation and verification of the coupled neutron kinetic/thermal hydraulic system code DYN3D/ATHLET. Ann. Nucl. Energy 84, 153e165. Krepel, J., Rohde, U., Grundmann, U., 2006. Dynamics of molten salt reactors. In: Int. Congress on Advances in Nuclear Power Plants - ICAPP 2006, 04.-08.06.2006, Reno, United States paper 6344. Krepel, J., Rohde, U., Grundmann, U., Weiss, F.-P., 2007. DYN3D-MSR spatial dynamics code for Molten Salt Reactors. Ann. Nucl. Energy 34, 449e462. Krepel, J., Grundmann, U., Rohde, U., Weiss, F.-P., 2008. Dynamics of molten salt reactors. Nucl. Technol. 164, 34e44. Krysl, V., Mikolas, P., Sprinzl, D., Svarny, J., 2010. MIDICORE VVER-1000 core periphery power distribution benchmark proposal. In: Proc. 20th Symposium of Atomic Energy Research (AER), Hanasaari, Espoo, Finland, September 20-24, 2010. Kuchin, A., Ovdiienko, I., Loetsch, T., 2007. Comparison of CASMO and NESSEL few group cross section libraries and their usage in DYN3D. In: 17th Symposium of Atomic Energy Research (AER), Yalta, Ukraine, 2007, Proc. KFKI Atomenergia, Budapest, ISBN 978-963-372-634-1, pp. 227e244. Kuchin, A., Ovdiienko, I., Khalimonchuk, V., 2010. A study of transient connected with WWER-1000 cluster rod drop with subsequent working of automatic power controller. In: 20th Symposium of Atomic Energy Research (AER), Hanasaari, Espoo, Finland, 2010, Proc. KFKI Atomenergia, Budapest, ISBN 978963-372-645-7, pp. 527e538. Lassmann, K., Hohlefeld, F., 1987. The revised URGAP-model to describe the Gap conductance between fuel and cladding. Nucl. Eng. Des. 103, 215e221. Lassmann, K., 1992. TRANSURANUS: a fuel rod analysis code ready for use. J. Nucl. Mater. 188, 295e302. Lepp€ anen, J., 2005. A New Assembly-level Monte Carlo neutron transport code for reactor physics calculations. In: Proc. M&C 2005. Avignon, France, Sept. 12-15, 2005. Litskevich, D., Merk, B., 2014. Development and verification of a transport solver for DYN3D. In: 45th Annual Meeting on Nuclear Technology, 06.-08.05.2014, Frankfurt, Germany.

Loetsch, T., Khalimonchuk, V., Kuchin, A., 2009. Proposal of a benchmark for core burnup calculations for a VVER-1000 reactor core. In: Proc. 19th Symposium of Atomic Energy Research (AER), Varna, Bulgaria, September 21-25 2009. Loetsch, T., Khalimonchuk, V., Kuchin, A., 2013. Consolidated data and status of task 2 solution of the benchmark for core burnup calculations for a VVER-1000 reactor. In: Proc. 23rd Symposium of Atomic Energy Research (AER), Strbske Pleso, Slovak Republic, October 30-November 4 paper 2.10. MacPherson, H.G., 1985. The molten salt reactor adventure. Nucl. Sci. Eng. 90, 374e380. Manera, A., Rohde, U., Prasser, H.-M., van der Hagen, T.H.J.J., 2005. Modelling of flashing-induced instabilities in the start-up phase of natural-circulation BWRs using the code FLOCAL. Nucl. Eng. Des. 235, 1517. Mittag, S., Kliem, S., Weiss, F.-P., Kyrki-Rajamaki, R., Hamalainen, H., Langenbuch, S., Danilin, S., Hadek, J., Hegyi, G., Kuchin, A., Panayotov, D., 2001. Validation of coupled neutron kinetic/thermal-hydraulic codes Part 1: analysis of a VVER1000 transient (Balakovo-4). Ann. Nucl. Energy 28, 857. Mittag, S., Petkov, P.T., Grundmann, U., 2003. Discontinuity factors for nonmultiplying material in two-dimensional hexagonal reactor geometry. Ann. Nucl. Energy 30, 1347e1364. Mittag, S., Grundmann, U., Weiss, F.-P., Petkov, P.T., Kaloinen, E., Keresztúri, A., Panka, I., Kuchin, A., Ionov, V., Powney, D., 2005. Neutron-kinetic code validation against measurements in the Moscow V-1000 zero-power test facility. Nucl. Eng. Des. 235, 485e506. Mittag, S., Kliem, S., 2011. Burning plutonium and minimizing radioactive waste in existing PWR. Ann. Nucl. Energy 38, 98. Nikitin, E., Fridman, E., Mikityuk, K., 2015a. Solution of the OECD/NEA neutronic SFR benchmark with Serpent-DYN3D and Serpent-PARCS code systems. Ann. Nucl. Energy 75, 492e497. Nikitin, E., Fridman, E., Mikityuk, K., 2015b. On the use of the SPH method in nodal diffusion analyses of SFR cores. Ann. Nucl. Energy 85, 544e551. NUREG/CR-5535/Rev 1, 2001. RELAP5/MOD3.3 code manual. In: Code Structure, System Models, and Solution Methods, December 2001, Idaho Falls, Idaho, USA, Volume I. Ovdiienko, I., Khalimonchuk, V., Kuchin, A., 2006. The Power's maneuvering regime simulation on 2nd unit of khmelnitsky NPP. In: 16th Symposium of Atomic Energy Research (AER), Bratislava, Slovakia, 25e29 September 2006. Ovdiienko, I., Kuchin, A., Khalimonchuk, V., Ieremenko, M., 2011. Effect of burnup history by moderator density on neutron-physical characteristics of WWER1000 core. In: Proc. 21st Symposium of Atomic Energy Research (AER)Dresden, Germany, September 19-23, 2011, p. 139. Petkov, P.T., Mittag, S., 2006. Evaluation of homogenisation error in two-group nodal diffusion calculation for VVER-1000 core. In: Annual Meeting on Nuclear Technology, 16.-18.05.2006, Aachen, Germany. Prasser, H.-M., 1983. Ein Beitrag zur Ermittlung der Kühlmittelmassenstromverteilung in Druckwasserreaktoren vom Typ WWER mit eingebautem elliptischem Siebboden. Kernenergie 26, 443e448, 11. €hne, T., Kliem, S., Rohde, U., Weiss, F.-P., 2003. Prasser, H.-M., Grunwald, G., Ho Coolant mixing in a pressurized water reactor: deboration transients, steamline breaks, and emergency core cooling injection. Nucl. Technol. 143, 37e56. Rachamin, R., Wemple, C., Fridman, E., 2013. Neutronic analysis of SFR core with HELIOS-2, SERPENT, and DYN3D codes. Ann. Nucl. Energy 55, 194. Rohde, U., Lucas, D., 1997. Solution of the 4th dynamic benchmark by use of the code DYN3D with a particle-in-cell method for the description of Boron €rnitz transport. In: Proc. 7th Symposium of Atomic Energy Research (AER), Ho (Germany), September 23-26, p. 387. Rohde, U., 2001. The modelling of fuel rod behaviour under RIA conditions in the code DYN3D. Ann. Nucl. Energy 28, 1343e1363. Rohde, U., Mittag, S., Ovdiienko, I., Ieremenko, M., 2004. Final report on BMU project €ge zur sicherheitstechnischen SR2441, sub-project 861500/9 e UA-2315 „Beitra Bewertung von KKW mit WWER-1000. Reaktorphysikalische Analysen“. HZDR internal report, 2004. Rohde, U., Grundmann, U., Kozmenkov, Y., Pivovarov, V., Matveev, Yu, 2006. Core design and transient analyses for weapons plutonium burning in VVER-1000 type reactors. In: Topical Meeting of German Nuclear Society (KTG-Fachtag) “Stand der Entwicklung für LWR Brennelemente und Auslegungsmethoden”, 02.-03.03.2006, Dresden, Germany, CD-ROM Proc., p. 265. Rohde, U., Mittag, S., Grundmann, U., Petkov, P., Hadek, J., 2009. Application of a stepwise verification and validation procedure to the 3D neutron kinetics code DYN3D within the European NURESIM project. In: Proc. Of the 17th Int. Conf. on Nuclear Engineering (ICONE17), Brussels, Belgium, July 12-16 2009. Rohde, U., Baier, S., Duerigen, S., Fridman, E., Kliem, S., Merk, B., 2012. Development and verification of the coupled 3D neutron kinetics/thermal-hydraulics code DYN3D-HTR for the simulation of transients in block-type HTGR. Nucl. Eng. Des. 251, 412e422. Rosenthal, M.V., et al., 1971. Advances in the development of MSBR. In: Proc. Fourth Int. Conf. On the Peaceful Uses of Atomic Energy, A/Conf.49/P/048, Geneva. Rypar, V., H adek, J., 1988. Space kinetics Experiment on the LR-0 reactor. Nukl. ÚJV     Re z, Special Issue, Re z 30e34. Smith, K.S., 1986. Assembly homogenization techniques for light water reactor analysis. Prog. Nucl. Energy 17, 303e335. Studsvik, 1995. CASMO-4-A Fuel Assembly Burn up Program, Version 1.28.05. Studsvik/SOA-95/1, 1995. Toumi, I., et al., 2000. FLICA-4: a three-dimensional two-phase flow computer code with advanced numerical methods for nuclear applications. Nucl. Eng. Des. 200, 139e155.