Recent neutronics developments for reactor safety studies with SIMMER code at KIT

The SIMMER family of codes is applied for safety studies of sodium fast reactors and reactors of other types. Both neutronics and fluid-dynamics parts of SIMMER are under development. In the paper new neutronics capabilities are presented. In particular developments for neutron transport solvers and a new technique for taking into account thermal expansion effects are described. These new capabilities facilitate 3D simulations and improve accuracy of modelling for the initiation transient phase during a hypothetical severe accident.


Introduction
The SIMMER family of codes is developed by JAEA (Japan) in cooperation with CEA (France), KIT (Germany) and other partners [1]. SIMMER-III is a 2D code version and SIMMER-IV is a 3D one. SIMMER includes neutronics, fluid-dynamics and structure parts. They are employed together to simulate reactor behavior. In particular sodium-cooled fast reactors (SFRs), lead-cooled, LBE-cooled, He-cooled fast reactors, experimental water-cooled reactors, accelerator driven systems and molten salt reactors were studied with SIMMER at KIT. SFR analyses at KIT are mainly done for 2D (RZ) models representing the reactor vessel. 3D (XYZ) simulations are time-consuming, but possible for parts of the vessel that include the core and its surroundings.
Transient simulations of hypothetical severe accidents in SFRs may include several phases, starting from the initiation phase, during which the fuel subassembly can-walls remain intact, that prevents the radial material movement in the core after coolant boiling, clad and fuel failure. The initiation phase may end after a power excursion, the amplitude of which depends on the reactor design and accident initiator. Then a transition to massive core melting may occur. The melting enables radial movements of molten/broken materials that may lead to subsequent re-criticalities and appreciable mechanical energy releases.
The neutronics part of SIMMER gets material and temperature distributions from the other parts and computes the power shape and amplitude employed in the other parts. The set of neutronics subroutines includes a subset that generates macroscopic cross-sections from SIMMER nuclear data libraries, a subset that performs reactivity and power amplitude computations and a subset that calculates time-dependent neutron flux and power shapes.

SIMMER cross-section libraries and validation of neutronics models at nominal conditions
Several cross-section libraries in the extended CCCC format [2] are employed in SIMMER studies at KIT. The libraries include data in ISOTXS format, in particular group-average cross-sections and group-to group scattering cross-sections, their angular Legendre expansion coefficients, and data in BRKOXS format, in particular Bondarenko f-factors. An 11-group library is usually applied for SFR production calculations, this library being produced directly from evaluated nuclear data files by employing a weighting function related to a SFR neutron spectrum. For treating water-cooled systems, a 40-group library is used. It includes cross-sections for 10 neutron groups below 1 eV and data on neutron up-scattering. Also a 72-group library is applied for fast reactor calculations for investigating effects related to the limited number of groups in the 11-group library and different nuclear data. The 40-group and 72-group libraries are obtained from a "master" 560 group library [3]. It was shown that for simple geometries, the 560-group library provides results very close to those obtained with pointwise cross-sections while employing a Monte-Carlo code [4]. Before transient simulations, the power shapes and major reactivity effects, including coolant density/void and Doppler effects, computed by SIMMER at steady-state are compared to those computed with reference neutronics codes, such as ERANOS [5]. In most cases the deviations between the compared values are acceptable in view of existing uncertainties, in particular due to nuclear data and modeling assumptions.
Note that XYZ geometry models are applied currently in SIMMER for 3D simulations, in plane each HEX being represented by two rectangles with the same area as HEX together. For the fluiddynamics part, the node interface areas are provided by the user in addition to other geometry data, while in the neutronics part, these areas between neighboring nodes of the same axial location are approximate. For an experimental benchmark on the EBR-II reactor with metal fuel, a very good agreement between XYZ and HEX-Z models was demonstrated [6] after applying VARIANT [7], a code that includes both XYZ and HEX-Z options. VARIANT is a 3D ERANOS neutron transport solver. Effects of using of XYZ geometry for neutron transport simulations in fast reactors with hexagonal subassemblies will be studied also in the future.

Improved quasi-static scheme in SIMMER
In SIMMER, the improved quasi-static scheme [8], IQS, is employed. The neutron flux is represented by a product of (1) a neutron flux amplitude function depending on time, and (2) a neutron flux shape function, depending on the time, spatial, energy (group) and angular variables. The shape function is assumed to be normalized, while establishing the IQS equations. The amplitude is computed by solving the point-kinetics equation. The equation for the shape is obtained from the neutron transport equation by substituting the neutron flux by the product of the amplitude and shape functions. Therefore the shape calculations can fully or partially compensate inaccuracies in amplitude computation, in particular when these inaccuracies are due to approximate weighting functions used for reactivity computations. The adjoint flux at steady-state is usually used as the weighting function.
The time-dependent equation for the shape function is solved by employing a fully implicit first order scheme for discretization in time. This approach transforms the time-dependent equation into a sequence of steady-state-like problems with external-source-like terms related to the values of the shape function at the beginning of time step and to delayed neutron precursors. In case of using the same time step for shape and amplitude computations, IQS is not much more efficient than the "direct" fully implicit method for solving the time-dependent transport equation for the neutron flux. When time-consuming shape calculations are performed less frequently as compared to fast-running amplitude calculations, the IQS method is faster than the direct one.
The point-kinetics parameters, such as reactivity and effective delayed neutron fractions, are recalculated during the transient by employing the time-dependent shape function. These calculations are done for so called reactivity steps, which are smaller than shape steps, but larger than the amplitude ones. The user may provide time-dependent "external reactivity" values, which are added to the computed transient reactivity for adjusting its value. One of the reasons is the existence of some phenomena affecting the reactivity, but not yet considered in SIMMER, such as e.g. thermal vessel expansion. If "external" reactivity values are used, corresponding corrections are applied also to the neutron production term in the equation for the shape function. For molten-salt reactors, the scheme was extended, the effects related to the movement of precursors being treated as an integrally negative neutron source term [4].

PARTISN-based neutron transport solver
The current SIMMER version includes a neutron solver based on DANTSYS, a static neutronics code [9]. The set of DANTSYS subroutines was introduced into SIMMER and extended -to solve timedependent equations for the shape function -at KIT in the past. Also the treatment of delayed neutrons was included. Unlike DANTSYS, PARTISN [10] includes an option for transient calculations. It is based on the semi-implicit scheme. The latter scheme is efficient for small time steps, but less suited for IQS if large shape steps are used. Therefore the time-discretization scheme in PARTISN was extended by implementing a theta-scheme, that is semi-implicit for theta=0.5 and fully implicit for theta=1, the theta parameter being provided by the user. A delayed neutron treatment was first introduced [11] via coupling of PARTISN with the KIN3D code [12]. KIN3D is a kinetics extension for VARIANT. The extended code was validated by simulations of fast transients in zero-power subcritical experiments performed earlier and also by comparisons with other codes [11]. Then the extended PARTISN code was coupled to SIMMER framework as an external executable, and it was demonstrated that SIMMER/PARTISN could be used to provide similar results as SIMMER/DANTSYS. [13]. Unlike DANTSYS, PARTISN is a parallelized code. A high efficiency of parallelization was demonstrated for 3D cases for a large number of neutronics nodes, the computation time per shape step being reduced by about a factor of 50 while employing 64 processors. For smaller number of nodes, in particular for 2D problems, smaller gains were observed employing a smaller number of processors, up to 8, while with more processors the computation time was not reduced. Further studies on parallelization are planned.

VARIANT-based neutron transport solver
An effort was also done to extend VARIANT capabilities and make possible its application as a neutron transport solver in SIMMER. Unlike DANTSYS and PARTISN, which are both based on the Sn method, VARIANT is based on the Variational Nodal Method and cannot be employed for modeling neutron transport in vacuum sub-regions. Though a complete vacuum does not occur in the systems currently considered with SIMMER, regions with a very low material density may appear after coolant elimination due to its boiling or a fission gas release from failed pins. For such regions, the VARIANT calculations may be less fast and accurate. A spatial kinetics extension for VARIANT was developed in the past [12]. More recently, options for application of large spatial nodes were studied at KIT. A technique for treating heterogeneous nodes was developed [14,15]. It allows application of nodes with different cross-sections for axial sub-nodes; e.g. a node with two sub-nodes, with and without a control rod absorber, can be treated. Using of larger nodes reduces the node number and accelerates calculations. More effort is needed to improve the VARIANT-performance, in particular to implement parallelization options.

Treatment of core thermal expansion effects in SIMMER
A conventional scheme for performing transient analyses for SFRs in Europe and Japan includes two steps: (1) SAS4A [16] or a similar code, based on the channel approach, is used for taking into account fuel irradiation effects and is employed for initiation phase simulation and (2) SIMMER is used for later transient phases. A special coupling code is used to transfer the values obtained at the end of the first step to SIMMER. The coupling introduces appreciable uncertainties because of quite different models used in different codes. The scheme applicability to the most recent SFR designs, in particular to those containing an axial heterogeneous fertile layer, is limited. Therefore an effort is under way to perform once-through analyses with SIMMER starting from nominal conditions, while SAS4A or special fuel performance codes should be used for irradiation simulations.
In view of the two-step scheme mainly applied in the past, SIMMER did not include before recently a model for taking into account reactivity effects due to core thermal expansion. These effects are important mainly for the initiation phase. For calculations with SIMMER starting from nominal conditions, these effects should be taken into account.
In the recently implemented in SIMMER-III scheme [17,18], the reactivity effects due to core thermal expansion can be taken into account. The contribution to the effect from each node is computed as the sum of two ones: (1) due to homogeneous expansion, i.e. same for all materials, including fuel, structure, coolant, etc., and (2) due to a relative variation of the coolant density compared to solid material densities. The second contribution is calculated on the basis of the assumption that the fuel mass, steel mass and coolant density do not vary due to expansion in any node. The SIMMER framework is already established for computing the second contributions.
For taking into account the first contribution, an "external reactivity"-like approach is used. The reactivity and neutron production term in the equation for the shape are modified, but the mesh is kept constant in fluid-dynamics and neutron transport calculations. The "external reactivity" due to expansion is computed in two steps. A technique based on the equivalence method [19] is applied in these computations. According to this method, instead of an isotropic system expansion determined by a factor f, one can consider an equivalent -with respect to the reactivity and other neutronics valuesnon-expanded system with densities divided by f.
First, the equivalence principle is applied node-wise, i.e. with different factors for each node (that is an approximation), so that one obtains an "equivalent" system with initial dimensions in plane, but with "equivalent" densities and axial dimensions. Second, additional transformations are done to bring also the axial dimensions to the original mesh, while making further modifications with material and related densities. Several options have been developed for the second step. One of the options is to establish equivalence coefficients between axial dimension and density variations globally or ringwise by performing k-eff calculations at the beginning of the transient with modified axial dimensions.
Alternatively, an "equivalent" system with initial axial dimensions can be obtained first and then additional transformations are done to take into account the radial mesh variations.
Thus the task is finally reduced to modifications of material/nuclear densities employed in the calculations and computing the reactivity effects induced by these modifications. In this way rather complex variations, e.g. related to bowing of fuel assemblies could be treated. After a structure mechanics code is coupled to SIMMER, the full potential of the incorporated scheme can be investigated. Currently effects due to expansions in plane leading to a cylindrical or conical core shape can be computed. Axially expanded dimensions are computed independently for different radial locations.

Conclusions
The SIMMER code is used and developed in Japan and Europe for simulating hypothetical core disruptive accidents in sodium cooled fast reactors and other systems. The code neutronics part has been extended recently to accelerate 3D calculations via employing a new neutron transport solver based on the PARTISN code. Another neutron transport code, VARIANT, is extended to treat axially heterogeneous nodes that facilitates its application as neutron solver in SIMMER as a smaller number of nodes can be treated. A new technique has been implemented in SIMMER for taking into account core thermal expansion effects. The extensions improve the SIMMER performance and applicability range.