DYNSUB: A high fidelity coupled code system for the evaluation of local safety parameters – Part I: Development, implementation and verification

Annals of Nuclear Energy 48 (2012) 108–122

Contents lists available at SciVerse ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

DYNSUB: A high fidelity coupled code system for the evaluation of local safety parameters – Part I: Development, implementation and verification Armando Miguel Gomez-Torres a,⇑, Victor Hugo Sanchez-Espinoza a, Kostadin Ivanov a, Rafael Macian-Juan b a b

Karlsruhe Institute of Technology, Institute for Neutron Physics and Reactor Technology, Hermann-von-Helmhotz-Platz 1, D-76344 Eggenstein-Leopoldshafen, Germany Technical University Munich, Department of Nuclear Engineering, Boltzmannstr. 15, D-85748 Garching, Germany

a r t i c l e

i n f o

Article history: Received 13 December 2011 Received in revised form 23 March 2012 Accepted 15 May 2012 Available online 26 June 2012 Keywords: Best estimate High fidelity Coupled code Simplified transport Sub-channel code Local safety parameters

a b s t r a c t DYNSUB is a novel two-way pin-based coupling of the simplified transport (SP3) version of DYN3D with the sub-channel code SUBCHANFLOW. The new coupled code system allows for a more realistic description of the core behaviour under steady state and transient conditions. The details of the developed internal coupling approach of both codes together with its implementation are presented and discussed. Comparisons of the results predicted by DYNSUB with the ones of coarser coupled solutions have shown very good agreement in the global parameters (keff and power distribution at steady state and position and magnitude of the power peak in the transient cases) validating the correctness of the coupling strategy. At local level however, important (and expected) deviations in the local safety parameters (maximal clad, fuel and moderator temperatures) have arisen. Differences up to 150 K in the centreline fuel rod temperature were found. It demonstrates the novel capabilities of the developed coupled system DYNSUB. The more detailed coupling solution had also an important impact on the convergence process, mainly in the neutronics internal convergence due to a smoother gradient on the thermal-hydraulics feedback parameters between neighbour sub-channels, compared to the gradient between assembly level channels. DYNSUB has successfully been applied to analyze the behaviour of one eight of a PWR core during a REA transient by a pin-by-pin simulation consisting of a huge amount of nodes. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction One of the main subjects related with the use of nuclear energy is Nuclear Safety. Nowadays, safety analyses are commonly carried out using a Best-Estimate (BE) approach by trying to simulate the physical phenomena taking place in the core, the coolant loops and the balance of plant as accurately as possible. In order to achieve the most realistic description of the neutron flux distribution and its coupling to the thermal-hydraulic phenomena within the reactor core, advanced multidimensional reactor dynamics codes have been developed and validated in the last decades. These state-ofthe-art codes are able to predict, for instance, non-symmetrical core power perturbations and they can calculate safety margins more accurately than the former developments (based on point kinetics) by means of three-dimensional core models with a spatial resolution at fuel assembly level for a wide range of operational transients and postulated accidents. Such methodologies do not allow the direct prediction of safety parameters base on local conditions (pin level) such as DNB and maximal cladding temperature

⇑ Corresponding author. Tel.: +49 721 608 2 3093; fax: +49 721 608 2 3718. E-mail addresses: [email protected], [email protected] (A.M. Gomez-Torres). 0306-4549/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2012.05.011

which is necessary in the frame of safety assessment of reactor cores. Some modern nodal reactor dynamics codes include own one dimensional thermal-hydraulics modules or they are coupled with best-estimate one dimensional system codes for the description of a wide range of transients taking into account important feedback mechanism between the involved domains (neutronics and thermal-hydraulics). The thermal-hydraulics model can be a one-dimensional two-phase flow model assigning parallel channels to each assembly or group of assemblies like FLOCAL (Manera et al., 2005) in DYN3D (Grundmann et al., 2005). In such coupled codes, usually fuel assembly averaged neutronics and thermalhydraulics parameters are exchanged during the solution (Ivanov et al., 1999a). In the last years, new developments are focused on more detailed spatial resolution of the neutronics and thermal-hydraulics domains using different approaches striving to have high fidelity coupled solutions based on local parameters (sub-channel/pin level). For example in Périn et al. (2010) by the coupling of a sub-channel code (COBRA-TF) with a nodal core simulator (QUABOX–CUBBOX) the cross flow between the FAs of a PWR can be taken into account. Moreover, a new generation of high fidelity codes is under development that will allow a more realistic prediction of local

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

safety parameters. Several strategies have been investigated to extend the coupled methodologies for a more detailed and physical prediction of the pin power (Ivanov and Avramova, 2007). Reactor dynamics codes like DYN3D and PARCS were extended by implementing a higher order solution of the neutron transport equations than the diffusion approximation (Beckert and Grundmann, 2008a,b; Grundmann, 2009; Downar et al., 2006), e.g. the transport approximation SP3 (Brantley and Larsen, 2000). On the other hand, discrete ordinates codes such as TORT-TD (Azmy, 1996; Seubert, 2011) are being improved for coupling with thermal-hydraulics codes (CFD, sub-channel and system codes). The coupling of advanced neutronics solutions based on SP3 or SN-approaches (allowing the prediction of the pin power) with 1D or 3D-coarse mesh thermal-hydraulics codes such as ATHLET and RELAP5 does not make possible a full coupled solution at cell/pin level since local thermal-hydraulic information is lost due to the averaging process during the calculation of the feedbacks. Consequently, a more precise and realistic description of both the thermal-hydraulics and the neutronics domains at the same spatial scale is needed in order to consider the increasing heterogeneities of modern core loadings in LWRs and future reactor designs. A new coupled code for the evaluation of local safety parameters, DYNSUB, has been developed by coupling the neutronics model of DYN3D-SP3 and SUBCHANFLOW (Imke et al., 2010) at fuel rod level. The details of the internal coupling approach of the two codes together with the implementation as well as selected results of the verification and validation work are presented and discussed. The comparison of the results predicted by DYNSUB with the ones of coarser coupled solutions have shown important deviations in predicting local safety parameters demonstrating the novel capabilities of the developed coupled system DYNSUB.

2. Description of utilized codes 2.1. DYN3D-SP3 DYN3D-SP3 is based on the well validated diffusion code DYN3D (Grundmann and Rohde, 1997; Grundmann and Kliem 2003; Grundmann et al., 2005; Grundmann, 2006; Rohde et al., 2009). The SP3 method is based on the developments of Brantley and Larsen (2000). In the PN-Method, the angular dependence of the neutron flux is expanded in spherical harmonic functions up to order N. The exact transport solution is recovered as N tends to infinity. In three dimensional geometries, the number of PN equations grows like (N + 1)2; in one dimensional planar geometry, the number of PN equations is only (N + 1). The PN equations in one dimensional planar geometry are relatively simple and can be

109

reformulated in second order form as (N + 1)/2 equations similar to the diffusion ones but coupled through the angular moments. In multidimensional geometries the case is too complicated. Although the problem can be formulated as a set of second order equations, in this case the number of equations increases and the coupling involves not only the angular moments but also mixed spatial derivatives of these moments. This increase in complexity led to the proposal of the Simplified PN approximation (Gelbard, 1960). In the SPN method, the second order derivatives in one dimensional planar geometry of the PN equations are replaced or ‘‘generalized’’ by means of the three dimensional Laplacian operator leading to a multidimensional generalization of the planar geometry PN equations that avoids the complexities of the full spherical harmonics approximation. In the DYN3D-SP3 code, a one dimensional thermal-hydraulic model (FLOCAL) is implemented. The neutronic solver is able to calculate a detailed multi-group pin power distribution based in the SP3 approximation. On the other hand, FLOCAL (Manera et al., 2005) is typically used as a 1-D assembly level TH module. 2.2. SUBCHANFLOW SUBCHANFLOW solves three mixture balance equations for mass, momentum and energy in axial direction as well as one additional lateral momentum equation at the sub-channel and fuel assembly level (Imke et al., 2010). It is based on the COBRA-family (Wheeler et al., 1976; Basile et al., 1999) of sub-channel codes but with improved capabilities. SUBCHANFLOW uses rigorously SI units internally in all modules. Coolant properties and state functions implemented for water correspond to the IAPWS-97 formulation. In addition, property functions for liquid metals (sodium and lead) and helium are available. SUBCHANFLOW calculates heat conduction in the fuel rod (pellet, gap and cladding) by means of a fully implicit finite difference method. For the heat flux transferred from the external clad wall to the fluid, the heat transfer coefficient, hcoef, is determined by using empirical correlations depending on the heat transfer mode and flow regimes. Additional constitutive relations are available for the void fraction, pressure drop, wall friction and turbulent mixing. 3. The coupling approach used in DYNSUB 3.1. Introduction The DYNSUB code is developed to overcome the limitations of the coupled codes relying on assembly based thermal-hydraulics models. A schematic view of such coupling approach is illustrated in Fig. 1. The neutronics solution is a detailed one, while the Thermal-Hydraulics (TH) solution calculates nodal averaged values at

Fig. 1. One-and-a-half way coupling.

110

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

Fig. 2. Two-way coupling in DYNSUB.

fuel assembly level. The calculation (update) of cross sections, at every iteration or time step, is done at pin level but with nodal averaged feedback parameters. By replacing the assembly averaged thermal-hydraulics model by a sub-channel code like SUBCHANFLOW, a two-way coupling at pin level is realized that aims to calculate and store the pin power in every axial node of each rod taking into account local feedback parameters in a direct manner avoiding the averaging process of the thermal hydraulic parameters. In Fig. 2, the two-way coupling scheme implemented in DYNSUB is illustrated for a fuel assembly. The pin power distribution is calculated with DYN3D-SP3 for every node in the fuel assembly geometry. The calculated pin power is transferred to the thermalhydraulics solver, where it is allocated to each axial node of the fuel rods as considered in the SUBCHANFLOW representation. In this internal coupling approach, each neutronics node is coupled directly to a thermal-hydraulics sub-channel or set of subchannels (usually one sub-channel considers four fuel rods located in the corner as it is illustrated in Section 3.4.1) allowing high fidelity detailed and direct pin-by-pin calculations of fuel assemblies (FA), FA’s clusters or cores. The development of DYNSUB and specifically the replacement of FLOCAL with SUBCHANFLOW required several changes and modifications of involved codes as well as the writing of new subroutines. In this implementation, SUBCHANFLOW acts as slave. After the compilation, a SUBCHANFLOW library is created and linked to the DYN3D-SP3 source, for a later compilation as a coupled code. 3.2. SUBCHANFLOW modifications and extensions For the internal coupling of SUBCHANFLOW with DYN3D-SP3 extensions and modifications were performed and implemented in such a way that it can be used as stand-alone code for pure thermal-hydraulics analysis of FA and cores but also for couple simulations. Hereafter the main changes of the code are described:  Source code extensions to have an additional option for direct allocation of the power distribution in every node of a 3D problem representation allowing a different axial profile for every single rod.  Implementation of a time dependent pin power map taking into account local changes. Such development replaces the former global definition of the power in which the change in power was taken into account as a global change (the original power multiplied by the fractional change in power).  Development of a generic and flexible Pre-processor for SUBCHANFLOW that is needed for the automatic input deck generation of very big problems (large number of channels

and rods in the model). This capability was also necessary for DYNSUB to represent accurately a core configuration or even a fuel assembly for which different tables with hundreds of rows (one for each sub-channel considered) have to be created. This Pre-processor is able to generate all the tables needed by SUBCHANFLOW with all the possible details e.g. the definition of an irregular cluster of assemblies each one with different inner configurations that are characterized by different type of rods with different thermal properties and control rod positions, with or without wetted boundary, etc. The Pre-processor works for the standalone version of SUBCHANFLOW and it is also directly implemented in DYNSUB as a subroutine. 3.3. DYN3D-SP3 modifications and extensions The changes of DYN3D-SP3 were mainly devoted to the replacement of the thermal-hydraulics module FLOCAL by SUBCHANFLOW. The general control of the iterative process follows almost the same logic described in (Grundmann et al., 2005). However, several modifications were necessary for the development of DYNSUB that are given hereafter:  New subroutine for the allocation of detailed 3D feedback parameters (at pin scale).  Modifications in the subroutines used for the calculation of pin-by-pin cross sections and for the interpolation with the new pin based feedback parameters.  Creation of subroutines for the handling of the detailed output of pinwise simulations keeping as much as possible the output structure of DYN3D-SP3. However, due to the more detailed output coming from SUBCHANFLOW, several modifications were implemented and are foreseen for future developments.  Modifications in the main subroutines performing the steady state and transient calculation in order to include the new thermal-hydraulics model. Initialization of new control variables and arrays, creation of SUBCHANFLOW input files by means of calling the Pre-processor subroutine, storage of old thermal-hydraulic values in case they are needed to iterate inside the time step. 3.4. The coupled code DYNSUB The development of multiphysics and multiscale coupled solutions requires the definition of a consistent spatial mapping between the involved domains. In addition, the temporal coupling and time advancement synchronization of the involved codes need to be solved in an appropriate way. Other issues such as the way of coupling (explicit, semi-implicit, full-implicit), or the check of the

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

convergence of the coupled solution are very important, too (Ivanov and Avramova 2007). In the next subsections, some details about the implemented coupling approach will be given. 3.4.1. Spatial coupling The spatial (radial/axial) mapping between the neutronics and thermal-hydraulics domains at a pin/sub-channel scale as well as the proper information exchange between these domains is one of the most challenging parts of the new coupled system. The mesh usually used in the sub-channel codes has normally a different nodalization compared with the neutronics one. Basically, they represent one sub-channel with four quarters of fuel rods in the corners. Such radial mesh is preferably used because most of the correlations used in thermal-hydraulics codes have been validated in experimental facilities with representative channels like this. SUBCHANFLOW is not limited to this channel definition, and subchannels coinciding with the neutronics nodes can be modelled; however this version of DYNSUB is fixed to traditional sub-channels and does not consider this option. In this coupling, the power coming from the neutronics code can be automatically allocated in the radial and axial rod power distribution used in SUBCHANFLOW. On the other hand, the thermal-hydraulics feedback parameters (TH-FB) are calculated in every sub-channel making necessary the use of a weighting method for having just one feedback value for each neutronics pin cell. As an example, in Fig. 3, a bundle of four rods with nine subchannels (SC) has been defined; furthermore the neutronic mesh (solid line) and the thermal-hydraulic mesh (dotted line) are shown. The radial mapping between the neutronics and thermalhydraulics domains relates four sub-channels to every neutronics node. Hence the TH-FB for each of the neutronics nodes must be calculated based on the corresponding four sub-channels before they are transferred to the neutronics part. The averaged thermal-hydraulic feedback parameter FP for every neutronic node ‘‘i’’ is given by means of the following equation: 4 X FPk  WðSC k Þ

FPi ¼

k¼1 4 X WðSC k Þ

ð1Þ

k¼1

where FPk is a given thermal-hydraulic feedback parameter calculated in the sub-channel SCk (moderator density or temperature, void fraction, etc.), W(SCk) is the weighting factor that depends on the thermal-hydraulic feedback parameter to be calculated. Thus, for the moderator density the normal weighting parameter is the volume of the sub-channel. For the temperatures different weight-

Fig. 3. Bundle of 4 fuel rods with 9 sub-channels.

111

ing parameters could be used, for the present development, the density of the sub-channels was used. The density yields more information than comparing just the volume especially in cases with changes in density in which the fixed volume does not change but the mass does. A better weighting factor could be the energy density or specific enthalpy mainly when the boiling appears in the calculation and the density becomes smaller and less representative but the energy increases. In DYNSUB the density was used, however it is foreseen to explore the use of the specific enthalpy for newer versions. FPi is the mean value of the feedback parameter transferred to the neutronics mesh. For the axial mapping, a similar weighting method can be applied. Axial weighting functions based on the geometry have to be calculated and used. However, all the results presented in this paper consider a one-to-one axial coupling (the same axial mesh is considered in both codes) making unnecessary the axial weighting. Once the feedback parameters are ready to be transferred, the communication process has to be defined. The following options were developed:  Data exchange via external files for the first version of DYNSUB based on Windows platform. Many subroutines for writing and reading of data were developed (Gomez et al., 2011). Although this way of data exchange and the weighting process were working fine, it is not computational efficient in case of large complex problems to be solved.  Data exchange of feedback parameters via memory arrays developed for a Linux platform aiming to solve real cores with huge memory requirements. 3.4.2. Temporal coupling The temporal coupling and time step selection plays also a very important role in the coupling. As a first approach in the coupling, a straightforward strategy was followed, i.e. to use a one-to-one time step selection. An explicit scheme (or also called marching scheme) was implemented. The one-to-one time coupling in DYNSUB is shown in Fig. 4. DYN3D-SP3 solves the time dependent problem with fixed cross sections (initial values or coming from the previous step). The power at the end of the time step is transferred to SUBCHANFLOW for its transient calculation. SUBCHANFLOW uses the two power maps (at the beginning and at the end of the time step), for the time calculation inside the step. In this way local changes of power inside the neutronics interval can be taken into account and even an internal division of the TH interval can be used. At the end of the thermal-hydraulics step, new feedback parameters are used for actualization of the cross sections for the next neutronics step. Later on, a nested loop iteration scheme was implemented as an approximation to an implicit scheme. In Part II of this paper, a full description of this temporal coupling together with some comparisons between the two methods is presented.

Fig. 4. Explicit coupling of DYN3D-SP3 with SUBCHANFLOW.

112

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

4. Investigations with DYNSUB

4.2. Case 1: Fast running mini-core

4.1. Introduction

4.2.1. Definition of the problem for Case 1 Fig. 6 shows the mini-core arrangement. It consists of a bundle of 2  2 fuel assemblies of UO2 at 4.2%. Reflective (North and East faces) and vacuum (South and West faces) boundary conditions are imposed. The initial operational conditions used in this case are presented in Table 1. For the upper and lower boundary conditions, in Seubert et al. (2008) and Christienne et al. (2010), water reflectors together with vacuum boundary conditions at the end were used on both sides. In this case, the lower and upper boundary conditions were adjusted by means of axial albedos in order to get a similar keff as reported in Christienne et al. (2010). At the initial steady state condition, the B2 assembly has a control rod fully inserted. Furthermore, such assembly has the largest burn-up (37.5 GWd/t). After a steady state calculation, the transient starts at time zero with a fast control rod ejection in assembly B2. The control rod must be fully ejected after 0.1 s. The transient continues until 1 s.

For testing the performance of DYNSUB two cases were defined.  Case 1: A fast transient test case (2  2 mini-core) based in the ‘‘OECD/NEA and U.S. NRC PWR MOX/UO2 core transient Benchmark’’ (Kozlowski and Downar, 2003) and previously defined and studied in (Seubert et al. 2008) and (Christienne et al., 2010) was chosen. A steady state and a control rod ejection transient were calculated and compared with the DYN3D-SP3 standalone results using FLOCAL with assembly based channels.  Case 2: In order to evaluate the practical feasibility of DYNSUB, a control rod ejection in an eighth core geometry (for the same Benchmark) is presented. Comparisons with the DYN3D-SP3 standalone using FLOCAL with assembly based channels are also performed. The cross sections library with 8 group pinwise homogenized cross-sections and kinetics parameters used in the two cases is provided in the NEMTAB-like format used in the OECD MSLB Benchmark (Ivanov et al., 1999b). The original library of the Benchmark problem was replaced with the library generated and validated in (Beckert and Grundmann, 2008a) for the development of DYN3D-SP3. In general, the 8 group structured library contains branches in fuel temperature (3), moderator density (3) and soluble boron density conditions (3) as well as 7 burn-up steps, to cover the expected range of core operating conditions. The moderator temperature effect is treated implicitly in the moderator density term assuming constant pressure of 15.5 MPa as described in (Kozlowski and Downar, 2003). The pinwise homogenized kinetics parameters are given only as a function of the burn up, i.e., they do not change with the thermal-hydraulic feedback. The configuration of the UO2 and MOX fuel assemblies is shown in Fig. 5. For control of reactivity, pins with Integral Fuel Burnable Absorber (IFBA) and Wet Annular Burnable Absorbers (WABA) are used in the UO2 and MOX fuel assembles respectively. The IFBA is a coating of zirconium diboride (ZrB2) on the fuel pellets and provides reactivity control over a relatively short burn-up period. All UO2 assemblies contain 104 IFBA pins located in the highest worth regions (vicinity of the guide tubes) and in the corners. The WABA is an annular pellet of Al2O3–B4C with wet (water filled) central region and Zircaloy cladding. In contrast to IFBA, WABA provides relatively long term reactivity control. The MOX fuel assemblies are designed with 24 WABA pins inserted in the guide tube locations (Kozlowski and Downar, 2003).

4.2.2. Steady state results for Case 1 In order to evaluate the rod worth of the assembly B2, two steady state calculations were performed with DYN3D-SP3 standalone and DYNSUB, all rods out (ARO) and all rods in (ARI). The rod worth is calculated by means of Eq. (2) (Christienne et al., 2010) and the insertion of reactivity in Dollars ($) by means of Eq. (3). ARO

rod worth ¼

q¼

ARI

keff  keff ARO ARI

keff keff

rod worth beff

ð2Þ

ð3Þ

The beta effective value ðbeff Þ for this core configuration is 695.7 pcm and it was obtained using Eq. (4). NASS X

bm eff

beff ¼ m¼1 NASS

ð4Þ

where bm eff is the beta effective for material m coming from the twogroup diffusion cross sections library of the OECD NEA PWR MOX Benchmark (Kozlowski and Downar, 2003), and NASS is the number of assemblies in the configuration. The results are presented in Table 2. Differences in keff up to 73.9 pcm for the ARO configuration and up to 88.2 for the ARI configuration were found, taken as a reference the DYN3D-SP3 calculation.

Fig. 5. Configurations for the UO2 (left) and MOX (right) fuel assemblies.

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

Fig. 6. Geometrical configuration of the 2  2 mini-core.

Table 1 Operational conditions for the 2  2 mini-core. Operational condition

Value

Number of assemblies Power (MWth) Inlet Temperature (°C) Inlet Pressure (MPa) Pressure drop over the core (kPa) Active Flow (kg/s) Fuel lattice, fuel rods per assembly Active length (cm) Assembly pitch (cm) Pin pitch (cm) Axial boundary conditions (Albedo) Number of axial nodes Boron concentration (ppm)

4 22.9 287.0 15.5 125.0 328.4849 17  17,264 365.76 21.41 1.26 a = 0.75 17 (21.515 cm) 1015.0

Table 2 keff, rod worth and reactivity (q) comparison between DYN3D-SP3 and DYNSUB for the 2  2 mini-core. DYN3D-SP3

DYNSUB

Deviation (pcm)

ARO keff ARI keff

0.980843

0.980104

73.9

0.971992

0.971110

88.2

Rod worth (pcm) q ($)

928.4 1.3344

945.0 1.3582

16.57 1.78%

Although the keff for DYN3D-SP3 in both cases is higher than the one calculated with DYNSUB, the worth of the rod for DYNSUB (945.0 pcm) is greater than the rod worth calculated with DYN3D-SP3 (928.4 pcm). The differences in keff can be explained as a result from the local differences in power due to the thermal-hydraulic refinement used in DYNSUB together with the hot full power (HFP) condition. The HFP condition results in thermal-hydraulic differences comparing DYNSUB with DYN3D-SP3 standalone as it will be illustrated and explained later in this section (Figs. 10–12). The beta effective value (beff) is in both cases smaller than the rod worth. Therefore a fast reactivity insertion leading the minicore to a super prompt critical state is expected. Fig. 7 shows for the ARI configuration, the convergence of keff (top-left), and how the deviations of the convergence parameters behave in the iterative process (the convergence criteria are defined by the user via input). Although the convergence criteria in keff (top-right) and in moderator temperature (bottom-right) were met relatively easy by the two codes, for DYN3D-SP3 standalone it was necessary much more iterations (41 in total) in order to fulfil the fuel Doppler temperature requirement (bottom-left). On the other hand, the thermalhydraulic behaviour of SUBCHANFLOW allowed DYNSUB to finish in just 10 iterations. It is important to mention that in order to validate the coupling, very strict criteria were used specially in the

113

fuel Doppler temperature (5  104), relaxing the convergence criteria would lead to a faster fulfilment of all the criteria by DYN3D-SP3 as it can be inferred in the trend of the fuel temperature deviation for DYN3D-SP3 (bottom-left plot). The axial power distribution of Fig. 8 shows a very good agreement between the two codes for the ARI configuration. The maximum is reached in the same axial node and only small differences between the two profiles can be distinguished. The radial pin power distribution for DYNSUB in the hottest layer (layer 7) for the ARI configuration is shown in Fig. 9. The effect of the control rod in the B2 assembly can be clearly observed. In Fig. 10, the deviation from the DYN3D-SP3 calculation in percent is shown. The differences range from approximately 1.0% in the coldest regions to +1.5% in the hottest regions of each assembly, with the understanding that a negative value implies an over estimation of the power by DYNSUB and therefore, a positive value results in an under estimation. Such differences arise due to the more detailed description of the thermal-hydraulic parameters coming from SUBCHANFLOW and are explained hereafter. Whereas DYN3D-SP3 predicts one value for each axial node (Fig. 11), DYNSUB is able to predict a 2-D distribution inside each axial node (Fig. 12). The fuel centreline temperature for the hottest assembly (A1 on the top-left) is in the case of DYN3D-SP3 1145.60 K. DYNSUB however, predicts the hottest temperature ranging from approx. 915 K in the coldest region of the assembly (bottom-right corner) until 1271 K in the hottest region (top-left corner). The bigger fuel centreline temperature predicted by DYNSUB in this assembly (top-left corner) results as an underestimation of the power in DYNSUB (due to the Doppler effect as well as the decrease in moderator density), likewise, the smaller temperature prediction (bottom-right corner) leads to a DYNSUB overestimation of power (see Fig. 10). Similar analysis can be done in the other fuel assemblies. It is important to stress that the small deviations in the normalized radial power distribution (from approximately 1.0% to +1.5%) shown in Fig. 10 resulted in 11% difference in the hottest centreline temperature predicted by DYNSUB (1271.0 K) in comparison with DYN3D-SP3 (1145.6 K). Further discussion related to these differences will be presented in the transient analysis results of Section 4.2.3. 4.2.3. Transient results for Case 1 As previously mentioned, the transient under study corresponds to a very fast control rod ejection. As a consequence, the rapid increase in reactivity brings the system to a prompt critical state and a significant power rise is expected in the initial few milliseconds. The strong negative fuel temperature (Doppler) feedback reactivity is the only inherent reactivity feedback responsible to avoid a bigger power excursion and to mitigate the power peak. The fuel temperature is a prompt and inherent reactivity feedback related to the Doppler broadening of resonances and is directly associated with the fuel temperature changes. The transient is started from critical condition by means of dividing the multiplication cross sections by keff. Since the difference in keff is just 88.2 pcm between the two calculations analyzed, the difference induced by such adjustment in the initial cross sections (0.0882%) can be neglected. However, if larger differences would have arisen, then it would have been necessary to start the transient with exactly the same set of cross sections in both cases in order to perform a real comparison. In order to evaluate the impact of the detailed thermal-hydraulic model several assumptions were made in order to be able to perform comparisons between DYN3D-SP3 and DYNSUB. The thermal-hydraulics models in FLOCAL and SUBCHANFLOW have different approaches and correlations implemented for the heat conduction in the fuel rod as well as for the heat transfer from

114

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

Fig. 7. keff and deviations of the convergence behaviour in the steady state calculation ARI of Case 1.

Fig. 8. Comparison of the normalized axial power profile for the ARI configuration.

the cladding to the fluid. However, for the sake of comparison, fixed and idealized properties were established for the heat conduction in the fuel for the two codes. Table 3 shows a resume of them. The heat transfer coefficient between the cladding and moderator is to be calculated using the own internal correlation of the both codes. Because this is a very fast transient, any variation may not have a great impact at the initial part of the transient. However, after the peak such variations may affect the asymptotic behaviour of the power. The effective Doppler temperature TDOPPLER is determined from the fuel rod centreline temperature Tf,C and the rod surface temperature Tf,S in agreement with the PWR MOX/UO2 core

transient Benchmark (Kozlowski and Downar, 2003) and is given by Eq. (5).

T DOPPLER ¼ ð1  aÞT f ;C þ aT f ;S

ð5Þ

The global behaviour of the power for the two calculations is presented in Fig. 13. As expected, a power peak was reached after the total ejection of the control rod. In Table 4, details about the results are given. The deviations were calculated taken DYN3D-SP3 calculation as reference with:

Dev iation ¼



ValDYNSUB  ValDYN3DSP3 ValDYN3DSP3

  100%

ð6Þ

115

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

Fig. 9. Normalized pin power distribution calculated with DYNSUB on the hottest layer (layer 7) for the ARI configuration.

Fig. 12. Thermal-hydraulic parameters in the hottest layer for the ARI configuration with DYNSUB: (a) fuel centreline temperature; (b) moderator temperature.

Table 3 Heat conduction properties for the fuel element.

Fig. 10. Percentual differences between DYN3D-SP3 and DYNSUB for the hottest layer (layer 7).

Fig. 11. Thermal-hydraulic parameters on the hottest layer for the ARI configuration with DYN3D-SP3.

Heat conduction properties

Value

Fuel thermal conductivity (kW/(m K)) Fuel density (kg/m3) Fuel specific heat (kJ/(kg K)) Cladding thermal conductivity (kW/(m K)) Cladding density (kg/m3) Cladding specific heat (kJ/(kg K)) Number of radial nodes in the fuel Heat transfer coefficient for the gas gap (kW/(m2 K)) Fraction of heat released in the fuel

3.0  103 10410.0 3.0  101 1.95  102 6504.0 3.0  101 10 10.0 1.0

The more detailed thermal-hydraulics effect in DYNSUB calculation increased the computational time in 26% comparing with DYN3D-SP3 calculation with FLOCAL on a fuel assembly based. On the other hand, although the rod worth calculated with DYNSUB was 1.78% bigger than the one calculated with DYN3DSP3, a power peak 0.46% smaller was predicted with DYNSUB. The reason can be explained by analyzing the Fig. 14 where the time evolution of the core-averaged Doppler temperature and moderator density are presented as a function of time.

116

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

Fig. 13. Global behaviour of total power during the transient.

Table 4 Comparison of DYN3D-SP3 and DYNSUB.

Power peak (MW) Peak time (milliseconds) Asymptotic power at the end of the transient (MW) Relative CPU time

DYN3DSP3

DYNSUB

Deviation (%)

879.651 105 59.5304

875.543 102 82.2091

0.46 2.85 38

1 (1530 min)

1.26 (1939 min)

26

It can be seen in the upper part of Fig. 14 that the core-averaged Doppler temperature predicted with DYNSUB increases faster than the one calculated with DYN3D-SP3. Thus, a faster insertion of negative reactivity is given in the case of DYNSUB. After reaching the power peak (102 ms), the decrease of moderator density, as a consequence of the fuel temperature increase, appears as an additional negative contribution to the reactivity. The moderator density predicted with DYNSUB decreases also faster than the one predicted with DYN3D-SP3. These two effects are responsible for the narrower shape of the DYNSUB peak in comparison with the DYN3D-SP3 peak. The faster increase of the Doppler temperature in the DYNSUB calculation comparing with DYN3D-SP3 is attributed to the impact of the HFP condition in the refined thermal-hydraulics approach used in SUBCHANFLOW. The differences between the maximal temperatures reached in some of the fuel rods predicted by DYNSUB in comparison with the fuel assembly averaged temperature predicted by FLOCAL increase significantly at the beginning of the transient. Due to this faster increase, a stronger Doppler effect is reached before in the hottest rods predicted by DYNSUB and therefore, the averaged Doppler temperature peak is reached faster and then remains under the DYN3D-SP3 curve. The discrepancies arising after the power peak and in the asymptotic part are due to the different heat transfer and thermal-hydraulics models and correlations used in SUBCHANFLOW and FLOCAL but also due to the impact of the more detailed thermal-hydraulic solution given by SUBCHANFLOW as it is discussed hereafter. In Fig. 15, the maximal fuel temperature (centreline) as a function of time for the two calculations is shown. Whereas DYN3D-SP3 calculates the assembly-averaged maximal fuel tem-

perature, DYNSUB is able to predict the rod with the maximum (rod 241) and the minimum (rod 17) fuel temperature within the hottest fuel assembly. Differences over 100 K between the assembly-averaged maximum fuel temperature (predicted with DYN3D-SP3) and the maximum fuel temperature of the hottest rod (predicted with DYNSUB) can be observed. However in the asymptotic part the difference in the maximum fuel temperature decreases. Fig. 16 shows the fuel centreline, the clad surface and the moderator axial temperature profiles at the DYNSUB peak time (Fig. 16a–c) and at the end of the transient (Fig. 16d–f) respectively. An axial shift can be seen by comparing the maximal fuel temperature profiles of Fig. 16 cases (a) and (d), where in the case of DYNSUB, the hottest axial layer has moved from axial level 6 at the peak time (102 ms) to axial level 5 at the end of the transient. At the peak time, differences up to 150 K in the centreline fuel temperature and up to 12 K in the cladding temperature are observed in all cases with an overestimation of DYNSUB. These differences become smaller as the transient progresses. At the end of the transient time there is still a DYNSUB overestimation of the fuel centreline temperature (about 60 K) but it is not the case of the maximal cladding temperature and the moderator temperature where the impact of the correlations used in the thermal-hydraulics models plays an important role. 4.3. Case 2: One-eighth sector of symmetry of PWR core 4.3.1. Definition of the problem for the Case 2 For a more realistic application, the core configuration of the OECD/NEA and U.S. NRC PWR MOX/UO2 core transient Benchmark (Kozlowski and Downar, 2003), assuming one-eighth symmetry, was modelled (Fig. 17). The benchmark was designed to provide the framework to assess the ability of modern reactor kinetics codes to predict the transient response of a core partially loaded with MOX fuel. A rod ejection may occur as a consequence of the rupture of the drive mechanism casing located on the reactor pressure vessel. This event is particular interesting for MOX fuelled cores since the delayed neutron fraction in MOX fuel is significantly smaller than in UO2 cores. The core has uniform fuel composition in axial direction. The original benchmark core considers axial reflectors (top and bottom) with the same width as the fuel assembly pitch and in

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

117

Fig. 14. Global behaviour (averaged over the whole mini-core) of the Doppler temperature and moderator density during the transient.

Fig. 15. Maximal fuel temperature during the transient.

the radial direction, the core is surrounded by a single row of reflector assemblies. In DYNSUB, since the use of reflector has not been implemented, axial boundary conditions by means of Albedos were chosen. The initial operational conditions used in this case are presented in Table 5.

4.3.2. Steady state results for the Case 2 Similarly to the mini-core, the rod worth of the assembly B5 was calculated with DYN3D-SP3 and DYNSUB. The average beta effective value ðbeff Þ for this core configuration (eighth of core) is 559.3 pcm and it was obtained by means of Eq. (4). The average

118

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

Fig. 16. Axial temperature distribution for the fuel centreline, the clad surface and the moderator at peak time (102 ms) and at the end of the transient (1000 ms).

Table 5 Operational conditions for the eighth of core.

Fig. 17. Geometrical configuration of the PWR core with one-eighth symmetry.

beta effective is in this case smaller than in the mini-core due to the introduction of fuel assemblies containing MOX. The results are presented in Table 6. Differences in keff up to 2.7 pcm for the ARO configuration and up to 3.2 pcm for the ARI

Operational condition

Value

Number of assemblies Power (MWth) Inlet Temperature (°C) Inlet Pressure (MPa) Pressure drop over the core (kPa) Active Flow (kg/s) Fuel lattice, fuel rods per assembly Active length (cm) Assembly pitch (cm) Pin pitch (cm) Axial boundary conditions (Albedo) Radial boundary conditions (Albedo) Number of axial nodes Boron concentration (ppm)

31 4.456  104 (104% rated power) 287.0 15.5 125.0 2545.69 17  17,264 365.76 21.41 1.26 aAX = 0.5 aR = 0.5 10 (36.576 cm) 1605.0

configuration were found, taken as a reference the DYN3D-SP3 calculation. These small differences, in comparison with the minicore, are due to the hot zero power (HZP) condition of the actual configuration. Therefore, it is clear that the thermal-hydraulic conditions do not play an important role at the steady state. The worth of the rod for DYNSUB (793.2 pcm) and the one of DYN3D-SP3 (792.7 pcm) are greater than the beta effective value.

119

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122 Table 6 keff, rod worth and reactivity (q) comparison between DYN3D-SP3 and DYNSUB for the eighth of core. DYN3D-SP3

DYNSUB

Deviation (pcm)

ARO

0.998378

0.998351

2.7

keff

ARI

0.990539

0.990507

3.2

Rod worth (pcm) q ($)

792.7 1.4172

793.2 1.4182

0.55 0.07%

keff

Table 7 Comparisons of DYN3D-SP3 and DYNSUB for the Case 2. DYN3D-SP3 Power peak (MW) Peak time (milliseconds) Asymptotic power at the end of the transient (MW) NK–TH iterations NK internal iterations Relative CPU time

DYNSUB 5

Deviation 5

1.28523  10 180 2.55044  103

1.27312  10 180 3.74688  103

0.94% 0.0 46.9%

524 108209 1 (9556 min)

524 103242 0.73 (7006 min)

0 4967 27%

is at HZP state and thus, the thermal-hydraulics conditions do not have a significant influence at this steady state. In Fig. 18, the percentage difference between DYN3D-SP3 and DYNSUB, for the radial power distribution in the axial layer 5, is shown. Just small differences ranging from 0.08% in the innermost region of the core (centre of the core) until 0.04% in the outer region of the core appeared.

Fig. 18. Pin power distribution, difference in percentage between DYN3D-SP3 and DYNSUB for the hottest layer (layer 5).

Thereby a fast reactivity insertion bringing the system to a super prompt critical state is expected. In this case both calculations got convergence in just 5 iterations. The iterative process was strongly related with the neutronics behaviour (deviations in keff) and not due to the thermal-hydraulics. This is explained again by the fact that the core

4.3.3. Transient results for the Case 2 In the same way that it was done for the mini-core, several assumptions were made in order to be able to perform comparisons between DYN3D-SP3 and DYNSUB. The heat conduction properties in the fuel pin were also fixed in agreement with Table 3. The transient is also started from critical condition by means of dividing the multiplication cross sections by keff. The global behaviour of the power as well as for the Doppler temperature and moderator density and temperature are presented as a function of time in Fig. 19. In Table 7, details about the results are given. The deviations were calculated taken DYN3D-SP3 calculation as reference. Similarly to the Case 1, a power peak slightly smaller was predicted with DYNSUB (0.94%). The averaged Doppler temperature in the first 200 ms increases faster with DYNSUB than with

Fig. 19. Global behaviour of total power and thermal-hydraulic parameters during the transient.

120

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

DYN3D-SP3 (top-right plot of Fig. 19). For this reason, the moderator temperature predicted with DYNSUB starts also increasing before the one predicted with DYN3D-SP3 (bottom-left plot of Fig. 19). After the first 200 ms, the Doppler temperature predicted by DYN3D-SP3 becomes greater than the one predicted with DYNSUB, thus this negative reactivity brings the power below the one predicted by DYNSUB. Important is to notice the impact of SUBCHANFLOW in the iteration process. Although the same number of Neutron-Kinetics/Thermal-Hydraulics (NK–TH) iterations was needed in the

two cases (524 iterations). The computational time was 27% larger in the case of DYN3D-SP3. The reason is that DYNSUB needed less internal iterations for the neutronic calculation (103242) than DYN3D-SP3 (108209) especially in the asymptotic region in which the power has already stabilized in its asymptotic value (Fig. 20). Similarly to the Case 1, in which a HFP state is analyzed, it can be seen that the more detailed thermal-hydraulics model in DYNSUB (by means of SUBCHANFLOW) has an impact in the CPU time. Further comparisons can be found in Gomez et al. (2011).

Fig. 20. NK internal iterations in every NK–TH iteration.

Fig. 21. Location of the hottest point in the reactor at assembly base and pin base with DYN3D-SP3 and DYNSUB respectively.

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

Like in the case of the mini-core (Case 1), the discrepancies arising after the power peak and in the asymptotic part are due to the different heat transfer and thermal-hydraulics models and correlations used in SUBCHANFLOW and FLOCAL and to the more detailed thermal-hydraulic description by means of SUBCHANFLOW. These differences resulted in a deviation of 46.9% in the power at the asymptotic part of the transient (from 250 until 1000 ms). For the prediction of the hottest point in the reactor configuration, discrepancies between the two codes were found. In the topleft graph of Fig. 21, the maximal fuel temperature (centreline) as a function of time for the two calculations is presented. Whereas DYN3D-SP3 calculates the position of the hottest fuel assembly (e.g. assembly B6 at 200 ms and C5 at the end), DYNSUB is able to predict the rod with the maximum (rod 6989) and the minimum (rod 7225) fuel centreline temperature within the hottest fuel assembly (assembly C7). DYNSUB predicted a maximum fuel temperature up to 62 K larger than DYN3D-SP3. The axial level with the hottest temperature is the same after the power peak with the two codes (top-right Fig. 21); however the radial location of the hottest fuel assembly did not match. DYNSUB predicted always the same fuel assembly (C7), whereas DYN3D-SP3 moved from assembly D6 until C5 (bottom-left Fig. 21). Finally, with DYNSUB is also possible to know exactly the hottest rod in the hottest assembly (bottom-right plot of Fig. 21). At the peak time, differences up to 15 K in the centreline fuel temperature and up to 4.3 K in the cladding temperature are observed in all cases with an overestimation of DYNSUB. The differences in the fuel centreline temperature, in opposite to the Case 1, become larger as the transient progresses (62 K at the end of the transient). However, the cladding temperature predicted by DYNSUB at the end of the transient is 5 K smaller. The differences in both cases are due to the impact of the more detailed thermal-hydraulics model. In DYNSUB a faster increase of global Doppler and moderator temperature can be observed in the first milliseconds of transient. However at the end of the transient, the global temperatures predicted by DYN3D-SP3 are larger than the ones of DYNSUB and although the fuel centreline temperature for the hottest rod remains bigger than the hottest assembly is important to notice that the difference between the hottest assembly temperature predicted by DYN3D-SP3 and the coldest rod in this assembly predicted by DYNSUB increased from approximately 100 K at peak time to 350 K at end time. Thus, at the end, although the temperature of the hottest rod (predicted by DYNSUB) remains over the hottest assembly, the average core centreline temperature is smaller.

5. Conclusions DYNSUB is a new best estimate and high fidelity coupled system, based in well validated stand alone codes. The more detailed thermal-hydraulics capability of DYNSUB avoids the loss of information due to averaging at assembly level, making possible to predict local safety parameters more accurately and directly i.e. without the use of hot channel factors. Calculations using DYNSUB have been performed and compared with the well validated code DYN3D-SP3 with its thermal-hydraulics model FLOCAL. In order to investigate the impact of the thermal-hydraulics refinement, FLOCAL has been used as a fuel assembly averaged model. The global and steady state results obtained with DYNSUB are in good agreement with the ones of FLOCAL. Comparing the local parameters predicted during the transient, it can be observed that DYNSUB predicts higher values of the temperature distributions than DYN3D-SP3. Such differences can be very significant while approaching to the operational limits. The

121

degree of details obtained with DYNSUB may allow a more accurate estimation of local safety parameters. It was also shown that the use of SUBCHANFLOW in DYNSUB has an impact in the iterative process. In the first case, the convergence was achieved in less NK–TH iterations in comparison with DYN3D-SP3, in the second case the impact was in the internal NK iterations. This impact is due to a smoother gradient on the thermal-hydraulics feedback parameters between neighbour subchannels, compared to the gradient between assembly level channels. With the LINUX version of the code it is possible to solve practical problems like the eighth of a PWR core involving a large amount of nodes. Further verification and validation studies must be done in order to assess the impact of the thermal-hydraulics models and correlations and to quantify the code’s uncertainty and sensitivity using e.g. the SUSA package.

Acknowledgement The results described in this paper are based on the valuable cooperation with the staff of the Institute of Safety Research of the Helmholtz-Zentrum Dresden-Rossendorf (HZRD), who gave us the source of DYN3D-SP3 and the cross sections library. Especial thanks to Dr. Sören Kliem, Dr. Ulrich Rohde, Dr. S. Mittag and Dr. Ulrich Grundmann.

References Azmy, Y.Y., 1996. The three-dimensional, discrete ordinates neutral particle transport code TORT: an overview. In: OECD/NEA Meeting on 3D Deterministic Radiation Transport Computer Programs, Feature, Applications and Perspectives. Paris, France, December 2–3. Basile, D., Chierici, R., Beghi, M., Salina, E., Brega, E., 1999. COBRA-EN, an Updated Version of the COBRA-3C/MIT Code for Thermal-Hydraulic Transient Analysis of Light Water Reactor Fuel Assemblies and Cores, Report 1010/1, ENEL-CRTN Compartimento di Milano. Beckert, C., Grundmann, U., 2008a. Entwicklung einer Transportnährung fürdas reaktordynamische Rechenprogramm DYN3D, Forschungzentrum DresdenRossendorf, Institut für Sicherheitsforschung, Berichts-Nr. FZD-497. Beckert, C., Grundmann, U., 2008b. Development and verification of a nodal approach for solving the multigroup SP3 equations. Annals of Nuclear Energy 35, 75–86. Brantley, P.S., Larsen, E.W., 2000. The Simplified P3 Approximation. Nuclear Science and Engineering 134, 1–21. Christienne, M., Avramova, M., Périn, Y., Seubert, A., 2010. Coupled TORT-TD/CTF capability for high-fidelity LWR core calculations. In: Proceedings of PHYSOR 2010, Pittsburg, USA, May 9–14. Downar, T.J., Xu, Y., Seker, V., 2006. PARCS v2.7U.S.NRC Core Neutronics Simulator, User Manual. University of Purdue. Gelbard, E.M., 1960. Application of Spherical Harmonics Methods to Reactor Problems, WAPD-BT-20, Bettis Atomic Power Laboratory. Gomez, A., Sanchez, V., Imke, U. and Macian, R., 2011. Pin level neutronic-thermalhydraulic two-way-coupling using DYN3D-SP3 and SUBCHANFLOW. In: Proceedings of International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, Rio de Janeiro, Brazil. Grundmann, U., Rohde, U., 1997. Verification of the code DYN3D/R with the help of international Benchmarks. Internal report Research Centre Rossendorf 195. Grundmann, U., Kliem, S., 2003. Analyses of the OECD main steam line break benchmark with the DYN3D and ATHLET codes. Nuclear Technology 142, 146– 153. Grundmann, U., Rohde, U., Mittag, S., Kliem S., 2005. DYN3D Version 3.2 Code for Calculation of Transients in LWR with Hexagonal or Quadratic Fuel Elements; Description of Models and Methods. Research Centre Dresden-Rossendorf. Grundmann, U., 2006. Calculations of a steady state of the OECD/NRC PWR MOX/ UO2 transient benchmark with DYN3D. In: Proceedings of KTG Annual Meeting on Nuclear Technology, Aachen, May 16–18. Grundmann, U., 2009. DYN3D – MG – V2.0: Code for Calculation of Steady States and Transients of Reactors by using the Multigroup Neutron Diffusion approximation for hexagonal or quadratic fuel Assemblies or the Multigroup SP3 approximation for quadratic fuel assemblies, Institute of Safety Research, Research Centre Dresden-Rossendorf. Imke, U., Sanchez, V., and Gomez, R., 2010. Subchanflow: A new empirical knowledge based subchannel code, In: Proceedings of KTG Annual Meeting on Nuclear Technology, Berlin, Germany, May 4–6.

122

A.M. Gomez-Torres et al. / Annals of Nuclear Energy 48 (2012) 108–122

Ivanov, K., Macian, R., Irani, A., Baratta, A., 1999a. Features and performance of a coupled three-dimensional thermal hydraulic/kinetics PWR analysis code TRAC-PF1/NEM. Annals of Nuclear Energy 26, 1407–1417. Ivanov, K.N., Beam, T.M., Baratta, A.J., 1999. PWR Main Steam Line Break (MSLB) Benchmark, Volume I: Final Specifications, NEA/NSC/DOC(99) 8. Ivanov, K., Avramova, M., 2007. Challenges in coupled thermal-hydraulics and neutronics simulations for LWR safety analysis. Annals of Nuclear Energy 34, 501–513. Kozlowski, T., Downar, T., 2003. OECD/NEA AND U.S. NRC PWR MOX/UO2 CORE TRANSIENT BENCHMARK – Final Specification. OECD Nuclear Energy Agency/ Nuclear Science Committee. Manera, A., Rohde, U., Prasser, H.-M., van der Hagen, T.H.J.J., 2005. Modeling of flashing-induced instabilities in the start-up phase of natural-circulation BWRs using the code FLOCAL. Nuclear Engineering and Design 235, 1517–1535. Périn, Y., Velkov, K., Pautz, A., 2010. COBRA-TF/QUABOX-CUBBOX: Code system for coupled core and subchannel analysis. In: Proceedings of PHYSOR 2010, Pittsburg, USA, May 9–14.

Rohde, U., Mittag, S., Grundmann, U., Petkov, P., Hadek, J., 2009. Application of a step-wise verification and validation procedure to the 3D neutron kinetics code DYN3D within the European NURESIM project; 17th International Conference on Nuclear Engineering ICONE-17, Brussels, Belgium, July 12–16. Seubert, A., Velkov, K., Langenbuch, S., 2008. The time-dependent 3D discrete ordinates code TORT-TD with thermal-hydraulic feedback by ATHLET models. In: Proceedings of International Conference on the Physics of Reactors, CasinoKursaal Conference Center, Interlaken, Switzerland, September 14–19. Seubert, A., 2011. The time-dependent 3-D transport code TORT-TD: Recent advances and applications, invited. Transactions of the American Nuclear Society, Vol. 104, Hollywood, Florida, June 26–30. Wheeler, C.L., et al., 1976. COBRA-IV-I: An Interim Version of COBRA for Thermal Hydraulic Analysis of Rod Bundle Nuclear Fuel Elements and Cores, BNWL1962. Pacific Northwest Laboratory.