3-D numerical studies of single control rod withdrawal transients in an LBE cooled critical reactor

3-D numerical studies of single control rod withdrawal transients in an LBE cooled critical reactor ☆ Xue-Nong Chen a,* , Yoshiharu Tobita a , Andrei Rineiski a , Barbara Kędzierska a , Guy Scheveneels b,1 , Matteo Zanetti b,2 , Bogdan Yamaji b a Institute for Neutron Physics and Reactor Technology (INR), Karlsruhe Institute of Technology (KIT), Hermann-von-Helmholtz-Platz 1, Eggenstein-Leopoldshafen D76344, Germany b SCK-CEN Belgian Nuclear Research Centre, Boerentang 200, 2400 Mol, Belgium ARTICLE INFO Keywords: LBE cooled reactor Control rod withdrawal 3-D thermal hydraulics-neutronics coupled simulations Dynamic reactivity Fuel pin degradation ABSTRACT The presented studies are carried out within the EU project ANSELMUS. A recent design version of MYRRHA (Multi-purpose hYbrid Research Reactor for High-tech Applications), which is a Lead-Bismuth-Eutectic (LBE) cooled reactor developed at SCK-CEN (Belgian Nuclear Research Centre), is investigated. The SIMMER-IV code is employed for 3-D simulations of single control rod withdrawal (CRWD) transients at the critical operation mode. A new CRWD model for the SIMMER-IV code is developed, so that the CRWD can be simulated for any constant withdrawal speed and from any initial position. The basic case is the complete withdrawal of a control rod (CR) filled with B4C absorber with the natural boron, where the reactivity worth is about 0.9 $, from a fully inserted position within 3 s. Cases with and without scram after 3 s are considered. The spatial kinetics effects on the power distribution are evaluated by comparing relative variations in time of local power densities and of the total one. The dynamic reactivity values during CRWD have been confirmed to be close to those obtained by static calculations. The transient with the scram at 3 s results in nothing severe, but that without the scram leads to local fuel melting. A further example, where the CR is filled with an enriched by B-10 absorber, its reactivity worth being 1.7 $, is calculated and shown as well. The withdrawal leads to a severe accident with fuel pin degradation, but without prompt supercritical power excursion. The numerical scenarios are presented and investigated by means of parametric studies. 1. Introduction The MYRRHA reactor (Baeten P, MYRRHA—Multipurpose hYbrid Research Reactor for High-tech Applications, SEARCH/MAXSIMA, 2014) designed by the Belgian Nuclear Research Centre is investigated currently in the EU project ANSELMUS (ANSELMUS Project. https:// www.anselmus.eu., 2024). ANSELMUS (Advanced Nuclear Safety Evaluation of Liquid Metal Using Systems) is a EURATOM project on safety assessment of heavy-liquid–metal-cooled reactor designs, which are currently developed in Europe. The present work is dedicated to the analyses of single control rod (CR) withdrawal transients in MYRRHA in its critical core configuration. The purpose of the project task is to estimate the 3-D effect of the single CR withdrawal and the deviation of the dynamic reactivity from the static one due to the CR withdrawal, in order to support the development of models for safety analysis, in terms of completeness of the model, its verification, and validation, required by the regulatory body (Courtois and Destin, 2013). An additional purpose of the paper is to understand what can happen, if the control rod withdrawal reactivity is appreciably larger than 1 $. To our knowledge, this kind of scenario has not been addressed before with the level of detail adopted in this work, in which transient simulations based on thermal hydraulics coupled to spatial kinetics with explicit CR movement are done. In general, the CR withdrawal is routinely considered in safety analysis, and examples for fast and other reactors can be found elsewhere. Whilst the studies demonstrate interest in the topic, they do not ☆ This article is part of a special issue entitled: ‘LFR’ published in Nuclear Engineering and Design. * Corresponding author. E-mail address: [email protected] (X.-N. Chen). 1 Current status: Retired 2 Current affiliation: Universit´ e libre de Bruxelles, Avenue FD Roosevelt 50, 1050 Brussels, Belgium Contents lists available at ScienceDirect Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes https://doi.org/10.1016/j.nucengdes.2025.114419 Received 31 March 2025; Received in revised form 10 July 2025; Accepted 20 August 2025 Nuclear Engineering and Design 444 (2025) 114419 Available online 30 August 2025 0029-5493/© 2025 Karlsruhe Institute of Technology. Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). address our concerns, reported above, but only the consequences of reactivity insertions, most often using point-kinetics approaches. There have been attempts over the years (Rineiski, 1997; Jeannin et al., 2024; Hirakawa et al., 1984) to address the spatial effects of a CR movement with a certain degree of explicitness rather than relying exclusively on a simple point-kinetics approach. This allows in some of these works to have a time-dependent directly evaluated reactivity, including nonseparable effect with respect to the change of the flux shape. However, these studies fall short from addressing in full our concerns, and tend to focus more on the techniques than on the specific issues that we want to address. This is particularly true for dynamic reactivity and what governs it. A Control Rod Withdrawal (CRWD) transient is a possible accident initiator in a nuclear reactor and should be studied for reactor safety assessment. In this paper we simulate such a transient with the SIMMERIV code (Yamano et al., 2008, 2008) for an LBE-cooled reactor design, MYRRHA, based on the numerical model that was made at KIT in the past (Kriventsev et al., 2014). We start from a 3-D model of the subcritical core configuration and update it to a critical configuration with a new CRWD model. The intention here is not to represent a specific design but to create the condition to emphasize the effect of the CRWD on the spatial changes of the neutron flux and the dynamics of the transients, that would otherwise be more similar to the ones predicted in cylindrically-symmetric cases if the rods are grouped in one ring. Proceeding on this basis, the updated critical core configuration model will be used in the analysis. The assumption is that one of the six control rods is withdrawn with a certain velocity. A new CRWD model is developed for SIMMER-IV. The aim is to study the dynamic reactivity evolution during such transient, and the dynamics and extent of the changes of power distribution in the core. After that, we explore the consequences of the CRWD on the core integrity. The SIMMER III/IV code (Sn Implicit Multifield Multicomponent Eulerian Recriticality), based on Advanced Fluid Dynamics Model (AFDM) (Bohl et al., 1990) and mainly developed by JAEA (Yamano et al., 2008, 2008,; Kondo S, Tobita Y, Morita K, Shirakawa N, SIMMERIII: an advanced computer program for LMFBR severe accident analysis, Proceedings of the International Conference on Design and Safety of Advanced Nuclear Power Plant (ANP’92), vol. IV, Tokyo, Japan, October 25–29, 1992), is particularly suited for studies of the transition phase of accidents, i.e. from the degradation of the pins to full core melting and core expansion. SIMMER IV is a three-dimensional fluiddynamics code coupled with a structure model including fuel pins, hexcans, etc., and a space–time and energy-dependent neutron transport Fig. 1. Geometry of the simplified CRWD model, where the area in orange is absorber in particle form. Fig. 2. SIMMER-IV simulation model, with the CR to be withdrawn, the neighboring FA and the peaking FA in the core center marked in yellow, blue and red, respectively. X.-N. Chen et al. Nuclear Engineering and Design 444 (2025) 114419 2 model. The overall fluid-dynamics solution algorithm is based on a timefactorization approach, in which intra-cell interfacial area source terms, heat and mass transfers, and the momentum exchange functions are determined separately from inter-cell fluid convection. In addition, an analytical equation-of-state (EOS) model is available to close and complete the fluid-dynamics conservation equations. The code has originally been developed and applied for safety analyses in the severe accident domain of sodium-cooled fast reactors (SFRs) (Yamano et al., 2008, 2008,; Kondo S, Tobita Y, Morita K, Shirakawa N, SIMMER-III: an advanced computer program for LMFBR severe accident analysis, Proceedings of the International Conference on Design and Safety of Advanced Nuclear Power Plant (ANP’92), vol. IV, Tokyo, Japan, October 25–29, 1992). However, the philosophy behind the SIMMER development was to generate a versatile and flexible deterministic tool, applicable for safety analyses of various reactor types with different neutron spectra and coolants, including critical reactors and accelerator driven systems (ADSs) for waste transmutation (Suzuki et al., 2005; Chen et al., 2017). SIMMER-IV includes a 3-D neutron transport solver based on the Sn method. The cross-section preparation model is based on the Bondarenko f-factor method. An isotope-wise 11-group library is employed for this study. With respect to time, an improved quasistatic scheme is used: the cross-sections and reactivity are recalculated at each reactivity step, the flux shape is recalculated at each shape step, which includes several reactivity steps. The fluid-dynamics meshes are treated as homogeneous material mixtures in neutronics. In the described simulations, we employ two fluid-dynamics meshes per subassembly in plane in order to represent the hexagonal arrangement with rectangular meshes. In this paper, the CRWD model development and the reactor simulation are described in §2.1 and §2.2. Then, we consider the following scenarios for the CRWD studies: (1) A single control rod withdrawal, where the control rod is filled with the B4C absorber of natural boron (for a worth of 0.9$) and it is completely withdrawn from an initially fully inserted position within 3 s. Afterwards, the reactor is considered to be scrammed. In addition, we study the consequences of the transient without scram, as well. The results are presented in §3.1. (2) A similar transient as the first one, except the control rod worth is increased to a worth of 1.7 $, by enriching the B4C absorber to 80 % B10, in order to make the scenario more severe and highlight the possible consequences. This transient leads to a power increase of about 18 times with pin degradation, but without prompt supercritical power excursion. After pin’s failure, the fuel is partly present in the coolant channels at the core centre in the form of particles. In addition, parametric studies on the blockage formation by fuel particles are performed and the results are presented and discussed in §3.2. 2. SIMMER-IV CRWD model and simulation 2.1. CRWD model development In order to study the spatial effects of the CRWD, we want to model the rod movement explicitly, rather than providing a reactivity as input. The SIMMER code is particularly suited for this study. As a multi-phase and multi-velocity field code, it can model material movement, including the movement of solids in the form of particles. This is achieved by treating particles as a “special” fluid. In the new CRWD model developed for SIMMER-IV, the absorber and structure materials, i.e. the B4C and the steel, in the withdrawn control rod are defined as particles, so that they can move. Then, the new model forces the control and steel particles to move uniformly at a certain user-defined velocity within a certain period at a given control rod location. This is an important development, since the standard code would allow modelling control and steel as particles, but they would be subject to interactions with surrounding fluids and structure, so that the coherent movement needed to represent the CRWD would be hampered. Fig. 1 shows the geometry of the simplified CRWD model. The new model is based on this geometry description. The volume fractions of absorber and steel particles are given in the conventional way as part of the SIMMER input file. Data for the CRWD model are provided as additional part of the input file, including the CR bottom position and its velocity. 2.2. Modeling of fuel blockage formation If the reactivity insertion due to CRWD is large enough, the fuel pin failure with possible clad and fuel melting can take place during the Fig. 3. Static CR withdrawal reactivity as a function of the axial position of the absorber upper boundary. X.-N. Chen et al. Nuclear Engineering and Design 444 (2025) 114419 3 Fig. 4. The power and reactivity transients in the natural B4C cases with and without scram. X.-N. Chen et al. Nuclear Engineering and Design 444 (2025) 114419 4 transient. Even after the fuel pin failure, in a lead or LBE cooled reactor, thanks to its high-boiling-point the coolant is still present as a liquid and rapidly cools the fuel, forming fuel particles. It is possible that these fuel particles, upon entering the fuel pin bundle flow path above the fuelfailure location, would form a blockage downstream, with a significant impact on the subsequent progression of events. SIMMER simulates blockage formation by particles using a ‘particle viscosity model’. This model represents the phenomenon of pressure-drop increase in the blockage, by considering that the effective viscosity of a mixture of liquid and dispersed solid particles increases exponentially with the volume fraction of solid particles. SIMMER handles two types of fuel particles. One is the standard ’fuel particle,’ which is formed when molten fuel cools and solidifies. The other is the ’fuel chunk,’ which consists of unmolten broken fuel pellets released into the flow path due to fuel pin failure. The particles are, by themselves, treated in the same manner by the code. The main difference, which is a user choice, is that they can be assigned to different velocity fields. In the SIMMER’s default model, molten fuel and fuel particles are assigned to the same velocity field, and the effective viscosity is determined solely based on the volume fraction of ‘fuel particle’, whilst the fuel chunk is assigned to another velocity field (say, “field two”). This stems from the need to analyse, in SFRs, the behaviour in which the mixture of molten fuel and fuel particles, formed in the central core region, infiltrates the gaps of the blockage composed of fuel chunks-fuel pellets fragments released due to the failure and collapse of fuel pins in the upper and lower core regions. Now, the effective viscosity is associated to the velocity field. Hence, in the default assignment of materials to the velocity fields, the effective viscosity is determined for “field one” solely based on the volume fraction of ‘fuel particle’, and in “field two” solely by the fraction of chunks. In leador LBE-cooled reactors, the issue arises of how to handle these two types of particles concerning the modelling of the effective viscosity of the mixture of molten fuel, fuel particles, and fuel chunks in a blockage. Compared to the known SFR case, in situations where failed fuel enters the pin bundle flow path in a leador LBE-cooled reactor, the analysis involves a phenomenon for which no experimental data exists, namely, the penetration of disrupted mixture of molten fuel and particles into a flow path filled with a high-boiling-point coolant. In this case, it is necessary to conduct analyses not only using the default model but also with conservative assumptions to evaluate the possibility of blockage formation. In the default particle viscosity model, only fuel particles are considered for the blockage model in “field one”. However, as a conservative assumption, an additional analysis is conducted in which fuel chunks are moved to “field one”. Under this condition, the solid particle volume fraction increases in “field one”, leading to a higher effective viscosity, which makes blockage formation more likely. 2.3. SIMMER-IV reactor simulation Fig. 2 shows the MYRRHA core map, with the CR to be withdrawn, the neighbouring fuel assembly (FA) near the CR and the peaking FA in the core centre marked in yellow, blue and red, respectively. The original 3-D SIMMER-IV model was prepared to simulate the sub-critical equilibrium core with 72 fuel assemblies at end-of-cycle (EOC) conditions, with different enrichments for each batch (Kriventsev et al., 2014). But now the core is loaded uniformly with MOX fuel with an overall Pu enrichment of 33 %, so that the reactor becomes critical. Consequently, the hottest fuel assembly (peaking FA) is the closest one to the centre. The effective delayed neutron fraction is 343 pcm and the CR to be withdrawn, marked in yellow, has an integral reactivity worth of about 0.9 $. The control rod is withdrawn downwards from a fully inserted position to a position outside the active core. On the fluid dynamics side the current SIMMER-IV model uses the original model (Kriventsev et al., 2014). The core inlet and outlet pressure and the core inlet temperature are given as boundary conditions, which are constant during all the transients studied in this paper, the mass flow rate and the coolant core outlet temperature are calculated. The mass flow rate has a radial distribution, which was adjusted by orifice coefficients (Kriventsev et al., 2014). 3. Simulation results and Discussion In the ANSELMUS project, the case of a control rod withdrawal from initially fully inserted position is considered. The control rod is filled with natural B4C and will be withdrawn completely within 3 s. This timeframe is considered representative of the typical response time of a protection system, but its actual value will need to be determined on the basis of safety analysis. It was intended to provide the maximum rod worth during such timeframe to amplify the neutron-kinetics effects under study. Afterwards in the protected case, the reactor is scrammed, which lasts also 3 s, that is just an arbitrary choice adopted in this study, the speed of insertion of the rod will depend on the protection system. In the unprotected case, the reactor will run continually without scram. Both protected and unprotected cases are considered and the results are shown in §3.1. In order to investigate the potential of the reactor to incur in and resist against a severe accident, we consider withdrawal of the CR filled with B-10-enriched B4C, which worth is almost doubled compared to the natural B4C control rod. For this case only the unprotected transients are calculated. Results are presented and discussed in §3.2. 3.1. CRWD transients with natural B4C First, we show the CRWD reactivity insertion as function of its position in Fig. 3, which is calculated statically as its position changes. The CR position is indicated with reference to the absorber upper boundary location. The fully inserted control rod is withdrawn at a constant speed of 0.31 m/s during a steady-state operation. It lasts 3 s, corresponding to an average reactivity insertion rate of about 0.3 $/s. There are two options afterwards (in 3 s after the CRWD): (1) the reactor core is scrammed and (2) not scrammed. The scram lasts another 3 s with introduction of a negative reactivity of 4 $. The power and reactivity transients in both Fig. 5. Temperature transients in both cases, where the fuel and coolant temperatures are the fuel pellet central temperature and the coolant temperature at the midplane of the hottest FA. X.-N. Chen et al. Nuclear Engineering and Design 444 (2025) 114419 5 Fig. 6. Material distribution at time points 20, 23, 30 s, where the CRWD transient starts at 20 s. No SCRAM performed. The radial mesh cell 3 is the CR location and the interior fuel melting in the core centre can be observed. X.-N. Chen et al. Nuclear Engineering and Design 444 (2025) 114419 6 Fig. 7. Transient results of FA power, where the CR withdraw starts at 20 s and no scram occurs. X.-N. Chen et al. Nuclear Engineering and Design 444 (2025) 114419 7 options are presented in Fig. 4. The total reactivity insertion in both cases is about 0.9 $, which can result in a fuel melting in the fuel pellet center at the hottest FAs in the case without scram, while no fuel melting takes place and the core stays intact in the case with scram. Fig. 5 shows the fuel center temperature and the coolant temperature in Kelvin at the midplane of the hottest FA. Fig. 6 shows the material distributions at 20 s (beginning of CRDW), 23 s (maximum fuel temperature is reached in the case with scram) and 30 s in the case without scram, where one can see the fuel melting at the radial meshes 8, 9 and 10. The plots in Fig. 6 present volume fractions of different material components for the core cross-section at the horizontal mesh K =16 (Y direction). This cross-section includes 17 plane (X direction) and 45 axial (Z-direction) meshes. If several components are present in a mesh, such as steel, fuel, coolant, the mesh includes several sub-regions of different color, their areas being proportional to the component volume fractions. Fig. 7 shows the power evolutions of the peaking FA, the neighboring FA (see Fig. 1) and the average for all FAs. Because of the positive reactivity insertion during the CR withdrawal, the FA power increases and then, due to Doppler feedback, stabilizes at a higher power level. The FA power is normalized by its initial value at 0 s. The power of the peaking FA clearly behaves almost in the same way as the average one, while the one of the neighboring to CR location FA increases by 15 % more than the average one. This is the local 3-D CRWD effect on power distribution. In order to compare the static and dynamic reactivity, we performed calculations for a new case: starting from a very low power of 10 kW, so that the temperatures of all materials remain almost unchanged up to about 3 s. The calculated total reactivity is the dynamic reactivity due to CRWD, this reactivity curve is presented in Fig. 8 as the red line with small dots. The dynamic reactivity is computed as a weighted reaction rate balance per fission neutron. We compare this dynamic reactivity with the static CRDW reactivity curve, computed with static k-eff calculations for the initial material and temperature distributions, but with different CR locations, as shown in Fig. 3, we see they are quite similar, but there is a small difference between them. To study it, we performed additional static k-eff calculations with the transient temperature and material distributions and computed then “static” reactivity values. The transient distributions were saved during transient calculations and then used for computing the cross-sections, then k-eff values. The result is shown in Fig. 8 as well with square red points denoted as “Static Calculation” in the legend. The differences between the static and the dynamic reactivity values do not exceed few cents, so the consistency of dynamic and static reactivity values is confirmed. The reactivity values computed in the transient and from static k-eff calculations are not exactly the same, because of the following reasons: •The transient reactivity values are obtained from the balance between the neutron production and losses, weighted with the scalar adjoint flux computed at the beginning of the transient. The static reactivity values are eigenvalues computed with the same crosssections as transient ones. Almost the same static reactivity values could be also obtained from the neutron balance, provided that the angular adjoint fluxes were computed at the transient times and used for weighting. Thus, one reason for the observed deviations is the use of the scalar adjoint flux computed at the beginning of the transient for dynamic reactivity calculations. •Next, the transient flux shape calculations are performed using an implicit scheme vs. time, thus by introducing a (1/v*dt) terms multiplied by flux in the neutron balance equation. This affects the flux shape if the flux amplitude is not constant. •The delayed neutron precursors are not in equilibrium with the prompt fission source if the flux amplitude is not constant. This influences the spectrum of the total, i.e. prompt plus delayed, fission source. Thus, another reason is that the transient flux is computed with a different effective fission spectrum. •The last reason is that the deviations between the transient and static values are comparable to the convergence criteria employed in the iterative procedures for the flux calculations. 3.2. CRWD transients with enriched B4C and blockage studies As already mentioned, we replace the absorber material of natural B4C with B-10 80 % enriched B4C in the control rod. The effective delayed neutron fraction is changed slightly to 345 pcm and the CR withdrawal worth becomes 1.69 $, which is almost doubled compared to the previous case. Fig. 9 shows the CR withdrawal reactivity for both cases of natural and enriched B4C. We consider the cases without scram, i.e. unprotected. As expected, the fuel pin failure will take place in the unprotected transients. Because of the large uncertainty in the numerical modeling, we consider two modeling options in this paper, namely the default (Case 1) and the conservative (Case 2). The case 2 is conservative with respect to the blockage formation. This is because both fuel particles (of smaller size) and fuel chunks (of larger size) are taken into for computing the viscosity in the particle viscosity model. Therefore, flow blockage takes place in Case 2. Due to the flow blockage, more pellets are broken in Case 2. Consequently, more fuel in the broken form is relocated from the core in Case 2. The reactivity and relative power are shown in Fig. 10 for both Case 1 and Case 2. The power transients are very similar for both cases. The power increases up to 18 times of the nominal power due to the CRWD reactivity insertion and then decreases because of the fuel dispersion after the fuel pin failure. The reactivity transients for different CR compositions are very similar before the fuel pin failure but then they Fig. 8. CRWD dynamic and static reactivities. Fig. 9. CRWD reactivities of natural and enriched B4C. X.-N. Chen et al. Nuclear Engineering and Design 444 (2025) 114419 8 Fig. 10. Reactivity and relative power in enriched B4C CRWD transients, where Case 1 denotes the default model and Case 2 the conservative one. X.-N. Chen et al. Nuclear Engineering and Design 444 (2025) 114419 9