© T. Le Meute et al., Published by EDP Sciences, 2025

1. Introduction

Stability of nuclear reactors is due to physical phenomena that induce a decrease in neutronic power when the power increases and an increase in neutronic power when it decreases. In homogeneous MSRs, these neutronic feedbacks are the Doppler effect and the density effect. When power increases, the temperature of the fuel rises, causing a cross-section broadening, which in turn induces a reactivity decrease for most of the nuclei (the increase of the absorption cross-section is stronger than the fission cross-section increase); this phenomenon is known as the Doppler effect. In a similar way, when the temperature rises, the density of the salt decreases, leading to a decrease in reactivity and power; this is referred to as the density effect.

The density reactivity feedback is an important aspect of chain reaction stability in MSRs due to the increase in leakage when temperature rises, particularly in Chloride MSRs [1, 2], where the Doppler effect is minimal (in [2], the Doppler feedback is estimated at ≈ − 0.25 pcm/K in Fig. 4.68 and the density feedback at ≈ − 19 pcm/K in Fig. 4.70). Understanding the behavior of these reactors during power transients is crucial for safety studies and modeling. During fast (short or comparable time to the pressure wave propagation in the system) power insertion transient, there is a pressure transient before the salt expands when considering compressibility. This induces a delay on the density feedback.

This work aims to describe the physical effects that cause a delay in the density reactivity feedback when considering the compressibility of the salt during variations in neutronic power. As previously emphasized in other works [3–5], neglecting compressibility leads to an underestimation of neutronic power computation during transients. The focus of this article is on a methodology to quantify the effects of compressibility on density reactivity feedback using a transfer function, facilitating computation in non-compressible calculation tools. To establish this transfer function, compressible calculations are conducted to model pressure wave effects. The tools used to implement the methodology are presented in Section 2. The methodology itself is described indetail in Section 2.3. The development and verification of this methodology rely on computations performed for two different cases, as detailed in Section 3. The results of the computation of these cases are detailed in Section 4.

The primary illustration for this paper of Molten Salt Reactors (MSR) is the Molten Salt Fast Reactor (MSFR) concept, a fast spectrum version of Molten Salt Reactors (MSRs) [6]. The MSFR is a 3 GWth breeder reactor in which heavy nuclei are dissolved in a fuel salt (Lithium Fluoride salt in the reference version, with a chloride version also under development [1]). This salt acts as both fuel and coolant, circulating from the core’s center, where the chain reaction is sustained, to the 16 sectors to transfer energy to the intermediate salt in heat exchangers (Fig. 1).

Fig. 1. MSFR concept [7].

2. Methodology

2.1. Thermalhydraulics computation

To perform compressible flow computations, the BlastFOAM solver [8], based on OpenFOAM [9], is utilized to solve the thermohydraulic problem. OpenFOAM has compressible equations of state for gas and for liquids with constant temperature but lacks equations of state that allows for considering pressure and temperature variations in liquids, that’s why BlastFOAM was chosen for this purpose as it is dedicated to detonation physics simulation. Additionally, BlastFOAM offers more tools for computing compressible flow, such as flux schemes. The solver has been adapted to introduce a time-linear energy source (constant power) in the fluid with an externally defined spatial shape. Notably, the neutronic code is not coupled with the thermo-hydraulic code in this first study. BlastFOAM introduces numerous elements to OpenFOAM, including various compressible equations of state (EOS). For this study, the linear Tillotson EOS [10, 11] is selected due to its linearity which is easier to parametrize when we don’t have robust experimental studies:

$$\begin{array}{c}\hfill p_t=p_0+\omega {\rho }_t(e_t-e_0)+A\mu +B{\mu }^2+C{\mu }^3\end{array}$$(1)

where p corresponds to the pressure, e to the internal energy, ρ to the density, and $\mu =\frac{{\rho }_t}{{\rho }_0}-1$, where the subscript 0 corresponds to the reference and t denotes the value at time t. A, B, C and ω are material coefficients. The coefficients used in this work are given in Section 3. The Riemann solver used in this work is the HLLC [12] flux scheme and the integration scheme the strong-stability-preserving Runge-Kutta method of order 4 (RK4SSP) [13], these are blastFOAM’s built-in methods.

2.2. Neutronic methods

In order to compute density reactivity feedback, we use the OpenMC code [14] associated with the Correlated Sampling (CS) [15] to estimate the reactivity variation due to a local fuel density variation.

These tools allow for the computation of spatial profiles of neutronic power shape (in the MSFR case) and the density reactivity feedback shape to be computed at each time step using the blastFOAM solver. The methodology looks like a ‘one-way coupling’, first computing for the case the shapes of the power and density reactivity feedback with the OpenMC code, then computing the thermohydraulic evolution with blastFOAM, and calculating the density reactivity feedback at each time step of the blastFOAM computation.

This methodology assumes that the neutronic feedback and power spatial shapes do not evolve with time. This hypothesis is valid for smooth cases, but for more extreme cases such as those studied here, its validity requires further research. However, adjustments can be made to this methodology in future works if it is shown that such evolution is an essential aspect of the transient behavior in order to capture the impact of the spatial deformation of the neutronic feedback and power spatial shapes.

2.3. Proposed transfer function methodology

The transfer function used in this work is inspired by the linearity of the wave equation, which permits to describe the pressure wave propagation. The density feedback effect is modeled as:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\rho }^{\ast }}{\mathrm{d}t}(t)={\int }_{-\infty }^{t}P(X)F_{\mathrm{insert}}(t-X)\mathrm{d}X\end{array}$$(2)

with ρ ^* [kg m⁻³] representing the density that affects the neutron population computed during an incompressible transient, P(X) [W] the power deposit at time X, and F_insert [kg m⁻³ s⁻¹ J⁻¹] denoting the transfer function computed for an elementary energy insertion.

To derive this equation, we assumed that a power deposed at time X during time dX induces a density reactivity feedback F_insert(t − X) at time t. Because the density reactivity feedback evolves with time, we must sum on all the power deposits during the computation which induces the integral in equation (2).

In order to compute a density reactivity feedback, we first need to compute (thanks to the ‘one-way coupling’ of OpenMC and blastFOAM) the response in density reactivity feedback (using the spatial feedback profile obtained with the OpenMC code) of an energy insertion E_insert during a time t_insert (short compare to the time of propagation of pressure wave), this response will be called t → FB_insert(t).

To validate this approach, we compute a less complex case: a constant power insertion. When an elementary computation of density reactivity feedback with time FB_insert is given for an energy insertion E_insert during the time t_insert, we will compute the density reactivity feedback FB^*(t) for an energy insertion E during the time t = N × t_insert (constant power insertion). The calculation of the integral of equation (2) becomes a sum and so the feedback is written:

$$\begin{array}{c}\hfill F{B}^{\ast }(t)=\frac{E}{E_{\mathrm{insert}}}\sum _{i=0}^{N}F{B}_{\mathrm{insert}}(t-i\times t_{\mathrm{insert}}).\end{array}$$(3)

This implies that the transfer function is a sum over time of energy insertions of the elementary FB_insert during the elementary time t_insert.

In future, the goal of this methodology is to model the compressibility impact on density reactivity feedback during reactor transient. The method represented above is applied at each time step to compute the density reactivity feedback during a transient. An algorithmic scheme proposing to implement this method is displayed in Figure 2.

Fig. 2. Algorithm proposed for the modelling of compressible impact on reactivity in incompressible codes.

3. Presentation of the cases

To describe the consequences of compressible flow on density reactivity feedback, two systems are solved: a simplified model consisting in a 1D tube reactor, and the MSFR. The end of this section presents the transfer function approach; this article aims to validate it.

For the following part, the parameters of the linear Tillotson model used in the computation to model the fuel salt of the reactor are listed in Table 1. These parameters were chosen to yield a dilatability of $\frac{\partial \rho }{\partial T}=-0.882$ kg m⁻³ K⁻¹ [16] and a speed of sound c  =  1200 m/s. The speed of sound is an estimation due to the lack of experimental data.

3.1. Tube reactor

The first case represents a long parallelipedic reactor of 10 m in length. The spatial power deposition is represented by a rectangular function cutting the reactor into three equal pieces. A side of the channel’s boundary conditions is set to outlet, the other one as a perfect wall and all the others are symmetry. The density reactivity feedback was calculated using the CS approach. These characteristics are illustrated in Figure 3. The spatial feedback is fitted using the following form: $FB(x)=Asin{\left(\frac{x\pi }{L}\right)}^B+C$, where x represents the position and L the length of the tube. One of the tube’s ends is closed by a wall, while the other is left open.

Fig. 3. Presentation of the cube reactor. Power shape, density and doppler feedback.

3.2. MSFR

The reference version of the MSFR is utilized in this article. The spatial shapes of power and neutronic feedback are depicted in Figure 4. In the MSFR, the walls are considered perfectly rigid, causing pressure waves to be reflected. The boundary condition on the expansion tank, i.e. the pipe on the upside of the critical zone of the core, is set as an outlet.

Fig. 4. Presentation of the MSFR. Power shape, density and doppler feedback.

The mesh is quite coarse for the MSFR case for two main reasons. The first is that we chose to use the same mesh for neutronic shape computation, to converge on a reasonable time; the mesh can thus not be more refined. The second is that we made a lot of tests, and we wanted the blastFOAM computation also in a reasonable time. Our aim in this paper is not to compute a high-fidelity transfer function for these cases, but to test and validate the proposed methodology.

4. Results

In this section, we discuss the results of energy insertion in both the tube reactor and the MSFR reactor. This work does not focus on the origin of energy insertion, but this is important to underline that reactivity insertion simulation is an aim of this work. A systematic study of reactivity insertion transient initiators can be found in [17]. These computations include the presence of compressible density reactivity feedback initially, and the discussion revolves around the validity of the transfer function described in Section 2.3.

4.1. Tube reactor

4.1.1. Transfer function generation

For the tube reactor with a length of 10 m and a speed of sound c = 1200 m/s, the time it takes for waves to travel from the center of the tube to the extremities is t_travel = 4.1 ms. To ensure that the energy insertion occurs much faster than t_travel, a time of t_insert = 10 μs is chosen, which is the same order of magnitude than the prompt neutron lifetime in reactors such as the MSFR. The evolution of density reactivity feedback is displayed in Figure 5. The energy insertion corresponds to a temperature increase of 1 K. In this figure, the dotted lines represent the incompressible evolution of each neutronic feedback, the dashed lines depict the compressible transient, and the solid line represents the total neutronic feedback.

Fig. 5. The density reactivity feedback from this reactivity insertion will be used for different temperature variation amplitudes and time scales.

These results properly illustrate the influence of pressure wave on density feedback, here the peak results from the reflection on the tube reactor wall of the pressure wave induced by heating (the peak between 5 ms and 12 ms corresponds to the travel of the pressure wave on the maximum value of the density feedback shape). At 15 ms, which corresponds approximately to 3 × t_travel (since the wave bounces, it takes 3 times for the wave to travel from the center to the limit of the reactor), the neutronic feedback computed, considering the compressibility of the salt, reaches the incompressible value. The density reactivity feedback induced by this energy insertion will be used as a transfer function to compute density feedback for energy insertion with higher time values such as 20 μs, 3 ms, and 20 ms.

4.1.2. Application of the transfer function

The results of these transients are displayed in Figure 6. The same transients are performed, but for an energy insertion corresponding to a temperature increase of 2 K. The results are displayed in Figure 7.

Fig. 6. The density reactivity feedback for the energy insertion equivalent to 1 K increase, time of the insertion of 20 μs, 3 ms and 20 ms in the tube reactor.

Fig. 7. The density reactivity feedback for the energy insertion equivalent to 2 K increase, time of the insertion of 20 μs, 3 ms and 20 ms in the tube reactor.

Depending on the energy insertion time scale, the feedback variation is smoothed for an insertion time of 20 ms. This is expected since it is like the time required for the wave to reach the boundary of the reactor. The impact of the bounce of the pressure wave on the wall is visible on each case.

In all of these results, the transfer function constructed in Section 2.3 yields exactly the same results as the direct computation of the density reactivity feedback, with a nearly instantaneous computation.

4.2. MSFR

4.2.1. Transfer function generation

The MSFR geometry is much more complicated than that of the tube reactor presented previously. Consequently, the propagation of pressure waves is much more complex, and the evolution of FB with time will be correspondingly more difficult to interpret. The time chosen for t_insert is the same as for the tube reactor, t_insert = 10 μs. The evolution of different neutronic feedbacks is displayed in Figure 8.

Fig. 8. The density reactivity feedback in the MSFR for the elementary energy insertion.

Interpreting the evolution of time in the density reactivity feedback computed by the neutronic codes is much more challenging than interpreting the results of the tube reactor. This difficulty arises because the walls of the core are modeled as perfect rigid walls, causing pressure waves to be completely reflected, the real reflection will depend on the material. As a result, these waves bounce and travel in the core for a long time. Furthermore, compared to the tube reactor, the tube leading to the expansion tank is very narrow, it takes time for the pressure waves to reach the free level and dissipate.

Even if the transfer of mass from the core center to the expansion tank is slow, the waves propagate, bounce and dilute in the recirculation loop of the reactor. Then the density drop due to the released energy is distributed throughout the entire circuit. For this reason, most of the feedback appears quickly due to the feedback shape in the system and the lower neutron importance near the core wall and in the recirculation loop. Note that, at 0.9 ms, the wave bouncing from the wall to the core center increases the density and induces a notable impact on the reactivity.

4.2.2. Application of the transfer function

In order to validate the methodology of the transfer function on a more complex configuration (here the MSFR), the same transients of energy deposition as those for the tube reactor are displayed in Figures 9 and 10.

Fig. 9. The density reactivity feedback for the energy insertion equivalent to 1 K increase, time of the insertion of 20 μs, 1 ms and 15 ms in the MSFR.

Fig. 10. The density reactivity feedback for the energy insertion equivalent to 2 K increase, time of the insertion of 20 μs, 1 ms and 15 ms in the MSFR.

The conclusions of these transients are the same as those for the tube reactor: the transfer function yields the same results for the evolution of time of the density reactivity feedback as the computed density reactivity feedback computed with compressible computation, even with much more complicated shape of the system.

5. Conclusions and perspectives

This work focuses on computing compressible flow transients resulting from energy insertion in MSRs, particularly examining the impact on the density neutronic feedback due to flow compressibility. These computations are conducted using the blastFOAM solver developed with OpenFOAM-9 coupled with the neutronic parameters computed by the Correlated Sampling approach.

A transfer function was developed in this work to compute the density reactivity feedback for different energy insertion times and values from a single energy insertion. This methodology was validated with two systems: a tube reactor and the reference version of the MSFR. In both cases, the results obtained using the transfer function were consistent with those obtained through direct compressible computation.

Based on this work, future developments will focus on integrating these results into incompressible CFD codes and system codes. These results have led to the development of a new modeling approach to incorporate the impact of the delay in density feedback due to compressibility into any non-compressible neutronic-thermohydraulic code without increasing computational time. With this methodology, a single compressible computation can be used to provide input to the transfer function for computing reactivity in incompressible codes, accounting for compressible transients. While this approach does not allow for the computation of parameters such as damages due to pressure waves on reactor structures, it enables the accurate computation of neutronic power and, consequently, temperature increase in the fuel salt. Another approach for the computation of this kind of transient is the coupling of a neutronic solver with a mechanical computation code such as Europlexus [3]. The coupling approach needs much more computation time but is more accurate. A validation of the method presented in this paper with that kind of coupling would be important. Especially on the interaction between pressure waves and neutrons population variation are the characteristic times of these phenomena can be very close.

Funding

This research did not receive any specific funding.

Conflicts of interest

The authors declare that they have no competing interests to report.

Data availability statement

This article has no associated data.

Author contribution statement

T. Le Meute: Conceptualization, Investigation, Methodology, Software, Writing – original draft. A. Laureau: Conceptualization, Validation, Supervision, Methodology, Software,Writing – review and editing. R. Moretti: Validation, Writing – review and editing. E. Merle: Validation, Supervision, Project administration, Funding acquisition, Writing – review and editing. T. Kooyman: Validation, Project administration, Writing – review and editing

References

All Tables

All Figures

	Fig. 1. MSFR concept [7].
In the text

	Fig. 2. Algorithm proposed for the modelling of compressible impact on reactivity in incompressible codes.
In the text

	Fig. 3. Presentation of the cube reactor. Power shape, density and doppler feedback.
In the text

	Fig. 4. Presentation of the MSFR. Power shape, density and doppler feedback.
In the text

	Fig. 5. The density reactivity feedback from this reactivity insertion will be used for different temperature variation amplitudes and time scales.
In the text

	Fig. 6. The density reactivity feedback for the energy insertion equivalent to 1 K increase, time of the insertion of 20 μs, 3 ms and 20 ms in the tube reactor.
In the text

	Fig. 7. The density reactivity feedback for the energy insertion equivalent to 2 K increase, time of the insertion of 20 μs, 3 ms and 20 ms in the tube reactor.
In the text

	Fig. 8. The density reactivity feedback in the MSFR for the elementary energy insertion.
In the text

	Fig. 9. The density reactivity feedback for the energy insertion equivalent to 1 K increase, time of the insertion of 20 μs, 1 ms and 15 ms in the MSFR.
In the text

	Fig. 10. The density reactivity feedback for the energy insertion equivalent to 2 K increase, time of the insertion of 20 μs, 1 ms and 15 ms in the MSFR.
In the text