© L. Clot et al., Published by EDP Sciences, 2025

1. Introduction

Molten salt reactors (MSRs) are innovative reactors where the fuel is a liquid circulating molten salt that acts also as coolant. The diversity of MSR concepts gives this family of reactors great versatility: the choice of the fuel salt (based on LiF for the fluoride versions or on NaCl for the chloride ones), of the neutron spectrum (from thermal to fast spectrum), the size and power allow the development of large breeder reactors as well as small burner AMRs (Advanced Modular Reactor).

MSRs have been studied at CNRS since the early 2000s and are of growing interest for the burning of Plutonium and minor actinides produced for example in pressurised water reactors.

The aim of this paper is to present the new developments made in the REM (Reactor Evolution with MCNP) code from 2023 onwards to improve the management of supplies in burner MSRs and to better simulate the physical and chemical aspects inherent to this type of reactor.

Section 2 presents the principle of neutronic depletion codes and details the former algorithm of the MSR related REM code. It presents some applications cases and model limitations too. The new developments carried out in the REM code are detailed in Section 3, followed by preliminary results on the CNRS’ reference MSFR (Molten Salt Fast Reactor). Then, Section 4 discusses the implementation of new physico-chemical constraints added to the REM code and their application in different supply strategies during the evolution.

2. Depletion simulations and the REM code

Section 2.1 presents the general principle of depletion simulations, while Section 2.2 gives a presentation of the REM code and its algorithm principle. Then Section 2.3 presents some application cases used below and model limitations that led to the new developments detailed in this paper.

2.1. Depletion simulations

The depletion simulations couple the steady state resolution of the neutron transport equation, solved in this work with a Monte-Carlo (MC) code, with an evolution module that solves the Bateman equations between two neutron transport steps. In the scope of this study, whole simulation is called depletion, whereas solving the Bateman equations is called evolution. Figure 1 shows the operating principle of the SMURE code [1], which can be generalised to depletion codes such as the one described in this work. Once the main parameters of the system under consideration and the nuclear data (as cross sections σ_i and time constants λ_i) have been defined, an initial Monte-Carlo calculation is carried out to calculate the reactivity and neutron flux ϕ at reactor start-up. These data are used to solve the Bateman Equations (1), during which certain conditions can be applied to control the evolution, such as setting the reactor power and/or the reactivity constant. Solving the Bateman equations will update the materials N_i to perform a steady state Monte-Carlo calculation again and so on until the end of the evolution.

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_i}{\mathrm{d}t}=\sum _{j\ne i}⟨{\sigma }_{j\to i}\varphi ⟩N_j-\sum _{j\ne i}⟨{\sigma }_{i\to j}\varphi ⟩N_i+\sum _{j\ne i}{\lambda }_{j\to i}{N}_j-{\lambda }_i{N}_i\end{array}$$(1)

Fig. 1. Principle of neutronic evolution calculation [1].

For molten salt reactors neutronic simulation, adaptations of the algorithm are needed. Equation (1) is established for the evolution of a core with no material addition or extraction during evolution, which is suitable for most studies and systems. For the MSR concepts based on continuous material supply and extraction during reactor operation, it is necessary to include additional terms to the Bateman equations. These terms are presented in Equation (2) [2]. Continuous processing is simulated with a time constant λ_i^processing while the supply is represented by a flow rate S​_i.

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_i}{\mathrm{d}t}=-{\lambda }_i^{\mathrm{processing}}{N}_i+S_i\end{array}$$(2)

Table A.1 in the Appendix presents the SAMOSAFER extraction schema developed during the PhD thesis of Hugo Pitois [3]. Two extractions are considered with two different time constants λ_i^processing injected in Bateman’s Equations (2):

• Chemical extraction with a λ⁻¹ of 450 days;

• Bubbling extraction with a small λ⁻¹ of 30 seconds.

2.2. The REM code

2.2.1. Overview

The in-house REM code has been under development since the early 2000s at CNRS/LPSC Grenoble [2, 4]. It is written in in-house objective C language within the DALI environment [5]. The code simulates the evolution of materials under constraints, as explained in Section 2.2.3, using the MCNP Monte-Carlo code [6] to solve the neutron transport and a Runge-Kunta 4 (RK4) algorithm for an evolution driven by precision. This code has been designed to be natively adapted for continuous loading reactors such as CANDUs and MSRs and has been tested and compared with other similar codes, notably as part of the European EVOL project [7]. Other codes have been developed since then by other actors such as MOSARELA [8], SMURE [1, 9], OpenMC [10], PYMS [11] and CEREIS [12].

2.2.2. Algorithm principle

The first step of the calculation consists of the definition of geometries, cells, materials and all the control and evolution settings and ranges for supplies. Then, during the code processing, there are different scales of calculation steps from:

• “Integration steps” where the MSR related Bateman’s Equation (2) is solved with a RK4 integration algorithm using variable steps;

• “Adjustment steps” where supplies S​_i are updated inside Bateman’s equation;

• “Control cycles” that manage adjustment steps for all supplies, as described in Section 2.2.3;

• “Neutron transport steps” where steady state transport calculations are performed to update reaction rates, as described in Section 2.2.4.

2.2.3. Adjustment steps and control cycles

In the REM code, the supplies are continuous while in most depletion codes, MSRs are supplied per batches. The flow rate F_k of each supply S_k is adjusted during the integration process at each adjustment step according to a constraint parameter C_k and its target value C_k^target. The physical parameters included in the code are as follows:

• Reactivity k_eff, calculated with interpolated fission $⟨{\overline{\nu }}_i{\sigma }_{\mathrm{f},i}\varphi ⟩$ and absorption ⟨σ_a, iϕ⟩ reaction rates (see Sect. 2.2.4) as shown in Equation (3):

  $$\begin{array}{c}\hfill k_{\mathrm{eff}}=\frac{{\sum }_i⟨{\overline{\nu }}_i{\sigma }_{\mathrm{f},i}\varphi ⟩N_i}{{\sum }_i⟨{\sigma }_{\mathrm{a},i}\varphi ⟩N_i}\end{array}$$(3)

  This estimation makes the assumption that there are no leaks, since it considers the entire system: it is therefore necessary to ensure that there are no leaks in the geometry during MCNP calculations;

• Proportion of heavy nuclei f_HN (Z ≥ 89) in the salt;

• Proportion of heavy nuclei plus Magnesium f_AnMg (Z = 12) in the salt;

• Proportion of transuranium elements f_TRU (Z ≥ 93) in the salt;

• Proportion of a specific element f_Z in the salt.

Each S​_k is supplied in a cell to keep C_k at C_k^target.

Each constraint C_k variation depends on the time t and the variation of supply flow rate F_k according to Equation (4).

$$\begin{array}{c}\hfill \mathrm{\Delta }{C}_k=\frac{\partial C_k}{\partial t}\mathrm{\Delta }t+\frac{\partial C_k}{\partial F_k}\mathrm{\Delta }{F}_k\end{array}$$(4)

So the variation in supply flow rate depends on the time and the variation in constraint according to Equation (5).

$$\begin{array}{c}\hfill \mathrm{\Delta }{F}_k=\frac{\partial F_k}{\partial C_k}(\mathrm{\Delta }{C}_k-\frac{\partial C_k}{\partial t}\mathrm{\Delta }t)\end{array}$$(5)

To calculate the variation in supply flow rate ΔF_k, each control cycle is subdivided into two adjustment steps:

• The first step is done without changing the supply variation to obtain the derivative $\frac{\partial C_k}{\partial t}=\frac{\mathrm{\Delta }{C}_k}{\mathrm{\Delta }t}$;

• The second step is used to calculate the value of $\frac{\partial F_k}{\partial C_k}$. To do this, the code will approximate the derivative $\frac{\partial C_k}{\partial t}$ with the sliding average of the last ℓ values to smooth evolution and use Formula (4) to calculate the last ℓ deviations from the target. In this way, the derivative $\frac{\partial F_k}{\partial C_k}$ is calculated as shown in Formula (6).

$$\begin{array}{c}\hfill \frac{\partial F_k}{\partial C_k}=\frac{\sum _{\ell \mathrm{last}}\frac{\mathrm{\Delta }{\mathrm{F}}_{\mathrm{k}}}{C_k^{\mathrm{target}}-C_k-\frac{\partial C_k}{\partial t}\mathrm{\Delta }t}}{\ell }\end{array}$$(6)

Once the two derivatives have been calculated, the supply flow rate is adjusted using the Formula (5) and the code moves on to the next supply in another control cycle.

2.2.4. Neutron transport/MCNP steps

Figure 2 shows the structure of the calculation and the organisation between MCNP steps (red rectangles with MCNP_n inside) and evolution steps (blue rectangles in Fig. 2 and dashed arrows between t_n and t_n + 1 in Fig. 3). Figure 3 illustrates this structure described below from the beginning to time t₂. At the beginning of the simulation, some parameters are given such as:

• The total time of simulation t_end;

• The length of the first evolution step Δt₁;

• The k_eff precision coefficient p_k;

• The interpolation function f for the reaction rates (as in the estimated k_eff calculation, Eq. (3)) among a ramp, a spline or the mean of the two values;

• The χ_max² and χ_min² coefficients used below.

Fig. 2. Principle of the organisation between MCNP steps and evolution steps in the REM code.

Fig. 3. Illustration of the organisation in Figure 2 for the first two evolution steps.

Each evolution step n has a dynamic length Δt_n adjusted as follows:

1. Starting from t_n − 1, if the code hasn’t reach the final time t_end of the simulation, a first evolution is performed to t_n = t_n − 1 + Δt_n using ⟨σϕ⟩_n − 1. At this time, REM knows the proportions N_i and the estimated k_eff (Eq. (3) and called k here) at t_n.

2. Transport is solved by running MCNP_n which gives the updated ⟨σϕ⟩_n and the MCNP k_eff with its standard deviation dk.

3. The absolute difference between MCNP k_eff and estimated k is compared with the precision coefficient p_k times the MCNP k_eff standard deviation dk: if the difference exceeds the defined criterion, this indicates that the estimated k is too far from the MCNP k_eff and its standard deviation dk, and the evolution step is remade with a reduced Δt_n (back to step 1.) then the simulation continues.

4. Back to t_n − 1, a second evolution is performed to t_n with an interpolation f between ⟨σϕ⟩_n − 1 and ⟨σϕ⟩_n. This evolution gives the proportions N_i′ at t_n.

5. A χ²-like test is performed between proportions N_i and N_i′ as shown in Equation (7).

   $$\begin{array}{c}\hfill {\chi }^2=\sqrt{\sum _{\mathrm{cells}}\frac{\sum _{i=1}^{n_{\mathrm{isotopes}}}{w}_i{(N_i^{\prime }-N_i)}^2}{{\left(\sum _{i=1}^{n_{\mathrm{isotopes}}}{N}_i\right)}^2}}\end{array}$$(7)

   The weights w_i are set to 1 by default and are currently set to 0 for the isotopes present in the supplies.

   • (a)

     If χ² >  χ_max² then another MCNP${}_n\prime$ is launched to re-update reaction rates and performed a third evolution between t_n − 1 and t_n to obtain proportions N_i′′. Then another χ² test is performed to compare N_i′ and N_i′′;

     • If the test is still unsatisfactory then the evolution step is remade (back to step 1.) with a reduced length: $\mathrm{\Delta }{t}_n/\sqrt{2}$;

     • Else, the next step is engaged with the same length: Δt_n + 1 = Δt_n.

   • (b)

     If χ² <  χ_min² and if the length of the previous step has not been reduced (as in cases 3. and 5.(a)) then the length of the next step is increased: $\mathrm{\Delta }{t}_{n+1}=\mathrm{\Delta }{t}_n\times \sqrt{2}$.

   • (c)

     Else, the next step has the same length and Δt_n + 1 = Δt_n.

2.3. Application cases and model limitations

2.3.1. Reference MSFR concept

The reference MSFR is a breeder reactor using the ²³²Th/²³³U cycle with a LiF–(HN)F₄ salt and maintaining the fraction of heavy nuclei at 22.5 mol% corresponding to the lowest melting point eutectic. This concept has been optimised during several European projects (EVOL [7], SAMOFAR [13], SAMOSAFER [14]). A more complete description of this concept can be found in [2, 7]. In the REM simulation of this reactor, there are three supplies correlated with three constraints showed in Table 1.

Reactivity is impacted both by the supply of ²³³U (direct effect) and by the breeding of ²³²Th (differentiated effect), guaranteed by a constant fraction of heavy nuclei in the salt. The 170 pcm gap between the target k_eff and criticality corresponds to prompt criticality, taking into account the circulation of precursors in the salt, a specific aspect of the considered molten salt reactors [15].

2.3.2. Chlorine MSFR concept

The Cl-MSFR is a breeder reactor developed during the PhD thesis of Pitois [3] funded by the SAMOSAFER project. Its aim is to transpose the reference MSFR but using the ²³⁸U/²³⁹Pu cycle with a NaCl–(HN)Cl₃ salt. In the performed REM simulation of this concept, there are six supplies correlated with the same amount of constraints showed in Table 2.

2.3.3. RAPTOr concept

The RAPTOr (Reactor supplied (“Alimenté” in French) with Plutonium and other Transuranic elements for Orano) concept is a burner reactor developed by Mesthiviers in her PhD thesis [16], as part of a collaboration between CNRS/LPSC Grenoble and Orano on burner MSRs. Its aim is to be supplied with Plutonium and minor actinides produced by light water reactors. At this stage of the study, a single Pu+AM supply is used to control reactivity, keeping it at −170 pcm.

2.3.4. Model limitations

Studies on burner MSRs have revealed certain limitations of the exiting model. Three main and correlated limitations have been identified:

1. Limitation on the types of constraints available as some seemed to be missing to carry out some studies;

2. To compensate the degradation of the isotopic fissile quality of their fuel, the fuel salt volume of burner reactors tends to increase, exceeding the physical volume of the core constrained by the geometry. Plus, the addition of new constraints–and therefore new supplies–tends to amplify this phenomenon, except of course for the constant ionic volume constraint (see Sect. 4.1).

3. Finally, the multiplication of constraints and supplies tends to make the model unstable.

This led to new developments in the REM code, including the implementation of new constraints and a dynamic expansion tank to manage excess salt (Sect. 4) and a new method of adjustment for the supplies (see Sect. 3 below).

3. New numerical developments

Section 3.1 presents the new algorithm implemented in the REM code, followed by preliminary results on the reference MSFR concept in Section 3.2, while Section 3.3 talks about dynamic expansion tank implementation in the code.

3.1. New algorithm implementation

The new adjustment strategy adopted for the REM code is to consider all the supplies at each control cycle. To do this, the n Equations (4)–with n number of supplies–are thought of as a single vector Equation (8) with ${\overline{\mathbf{J}}}_{\mathbf{C}}(\mathbf{F})$ the Jacobian matrix of the system defined by Formula (9), F the flow rate vector and C(F) the associated constraints vector. Now all the constraints depend on all the supplies, while the old algorithm acted as if the Jacobian matrix were diagonal.

$$\begin{array}{cc}& \mathrm{\Delta }C(F)=\frac{\partial C(F)}{\partial t}\mathrm{\Delta }t+{\overline{J}}_C(F)\mathrm{\Delta }F\hfill \end{array}$$(8)

$$\begin{array}{cc}& {\overline{J}}_C(F)=\left[\begin{array}{ccc}\frac{\partial C_1(F)}{\partial F_1}& \cdots & \frac{\partial C_1(F)}{\partial F_n}\\ ⋮& \ddots & ⋮\\ \frac{\partial C_n(F)}{\partial F_1}& \cdots & \frac{\partial C_n(F)}{\partial F_n}\end{array}\right]\hfill \end{array}$$(9)

So Equations (5) become now vectorial Equation (10).

$$\begin{array}{c}\hfill \mathrm{\Delta }F={\overline{J}}_C^{-1}(F)(C^{\mathrm{target}}-C(F)-{\langle \frac{\partial C(F)}{\partial t}\rangle }_{\mathrm{even}}\mathrm{\Delta }t)\end{array}$$(10)

Each control cycle comprises 2(n + 1) adjustment steps. The odd steps are used to calculate the derivatives forming the columns of the Jacobian matrix by modifying a supply and looking at the impact of this change on all the constraints. The even steps are used to calculate the time derivatives without adjusting the supplies (ΔF = 0) to finally take their average over the cycle to calculate Formula (10). Figure 4 shows the Chlorine supply in a RAPTOr-like simulation with three supplies: NaCl, MgCl₂ and PuCl₃, the Chlorine supply being a condensed image of these three supplies. The control cycles are shown in red and the adjustment steps in blue. Each change in Chlorine flow rate corresponds to an adjustment step and once the 8 = 2(3 + 1) steps have been completed, all the supplies are updated at the same time and a new control cycle begins.

Fig. 4. Details of the different steps during an MCNP step on a Chlorine supply (RAPTOr-like concept).

3.2. Preliminary results on the reference MSFR case

Figure 5 shows the in-core Uranium supply (flow rate and cumulative) and the reactivity in the MSFR core during an evolution, with the old controls as implemented at the beginning of 2023 (in black) and with the new controls (in red).

Fig. 5. Supply and cumulative supply of 233U and reactivity in the MSFR core. (a) Supply of 233U in the MSFR during 200 years of evolution. (b) Supply of 233U in the MSFR during 2 years of evolution. (c) Cumulative supply of 233U during 200 years of evolution. (d) Cumulative supply of 233U during 2 years of evolution. (e) Reactivity in the MSFR core during 200 years of evolution. (f) Reactivity in the MSFR core during 2 years of evolution.

Despite the radically different altitudes of the two supply curves (Fig. 5a), the supplies are equivalent on average, as shown in Figure 5c, where the cumulative supplies only deviate by 2 wt% after 200 years (from 22 419 to 21 974 kg). The new method of adjusting the supplies has tended to stabilise them out compared with the old method. Figure 5b shows that the new supplies have a continuous character that the old ones did not have. This can be seen clearly in Figure 5d, where the new cumulative supply is smoother than the old one, which is crenellated. Finally, it is possible to see in Figure 5b small peaks on the new supply: these are the small changes used to calculate the derivatives of the column of the Jacobian matrix as illustrated in Figure 4. By way of illustration, Formula (11) shows the numerical value of the Jacobian matrix during a control cycle at the end of the evolution, showing that the approximation made by the old adjustment model that this matrix was diagonal is not correct in this case.

$$\begin{array}{c}\hfill {\overline{J}}_C(F)=\left[\begin{array}{ccc}56.91& -6.27& -1.24\\ 2.38\times {10}^{-5}& 2.41\times {10}^{-5}& -2.30\times {10}^{-9}\\ -1.87\times {10}^{-9}& -1.08\times {10}^{-10}& 5.92\times {10}^{-5}\end{array}\right]\end{array}$$(11)

Making the supplies more stable has positive effects on the constraints they control, as illustrated in Figures 5e and 5f, which shows that core reactivity oscillates much less around its target value and is slightly more responsive.

3.3. Dynamic expansion tank implementation

In addition to developing new controls, the ability to dynamically modify the geometry of a core cell–known as the expansion tank–has been implemented. This process takes place at each MCNP step: when the parameter is activated, the REM code calculates the ionic volume of the fuel salt and adjusts the height of the fuel salt level to match the geometric volume when the transport is solved. Figure 6 shows two visualisations of the geometry of the Cl-MSFR concept at two instants of an evolution. The fuel salt is in purple, the fertile blanket is in pink, the B₄C neutron protection is in blue, the steel reflectors are in green and the argon sky in the empty part of the tank is in yellow. As the volume of fuel salt increased during evolution (see Sect. 2.3.4), the height of the salt level in the vessel increased accordingly.

Fig. 6. MCNP visualisation of Cl-MSFR [3] geometry at the beginning of the evolution (left) and after 50 years (right).

4. New physical models, supply strategies and applications

Section 4.1 lists the new constraints implemented in the REM code, then Section 4.2 presents and compares three supply strategies. Finally, Section 4.3 shows the impact of extraction scheme on depletion simulation.

4.1. Implementation of new constraints

As mentioned in Section 2.3, new constraints have been added to the REM code:

• Proportion of alkali metals f_A in the salt;

• Proportion of alkaline earth metals f_AT in the salt;

• Proportion of alkali and alkaline earth metals f_AAT in the salt;

• Constant ionic volume V of the salt;

• Isotopic grade of Plutonium;

• Electroneutrality of the salt.

The implementation of these new constraints was motivated by new considerations in fuel salt chemistry. For example, taking into consideration the ternary salt NaCl–MgCl₂–PuCl₃, the proportions of alkali and alkaline-earth will make it possible to control the position of the composition on the ternary, in particular to remain close to the eutectic minimum melting temperature of the salt due to safety issues. Figure A.1 in the Appendix shows this ternary salt, assuming that all the actinides behave like Plutonium in this system (i.e Pu = An) and highlighting the eutectic in red.

Controlling the constant volume ensures that the supplies are equivalent in volume to the extractions and could eliminate the need for a dynamic expansion tank. The ionic volume V is calculated with the ionic radii r_i and the atomic proportions N_i weighted with a coefficient comparing the ratio between the geometric volume V_geom and the initial estimation of the ionic volume V₀ as shown in Formula (12). For hydrogen and noble gases–which have no ionic radius–the atomic radius is taken instead as a first approximation.

$$\begin{array}{c}\hfill V=\frac{4}{3}\pi \sum _i{r}_i^3{N}_i\times \frac{V_{\mathrm{geom}}}{V_0}\end{array}$$(12)

Finally, the aim of controlling the electroneutrality of the salt is to match the theoretical valencies of the elements in the salt with those of the halogens. Formula (13) shows the result of the fission of a Plutonium nucleus into two fission products FP¹ and FP² with the distribution of the three initial Chlorine atoms. If the sum of the valencies of the two fission products is not equal to 3, then the resulting lack or excess of Chlorine can influence the reducing or oxidising nature of the salt, which has a direct impact on its corrosive nature on structural materials. To this end, the theoretical valency values of all the elements (positive for cations and negative for anions) have been statically implemented in the code to measure any discrepancies, which could lead to further studies to update them.

$$\begin{array}{c}\hfill {\mathrm{PuCl}}_3⟶({\mathrm{FP}}^1){\mathrm{Cl}}_x+({\mathrm{FP}}^2){\mathrm{Cl}}_{3-x}\end{array}$$(13)

4.2. New supply strategies for the evolution simulations

The new capabilities of the REM code are illustrated in this section on the RAPTOr case. The proposed strategies are first introduced (Sect. 4.2.1) then the processing scheme used for the simulations is presented (Sect. 4.2.2), and finally the results between one, three and four supplies strategies are compared (Sect. 4.2.3).

4.2.1. Three and four supplies strategies

Since the increase of the number of constraints–and therefore supplies–has been made possible without altering the stability of the model, it is possible to think about new supply strategies for the RAPTOr case. As well as maintaining reactivity, there are two interesting constraints to consider: staying close to the eutectic and maintaining a constant volume. Assuming the location of the initial salt, staying close to the eutectic is equivalent to maintain a constant proportion of Sodium–so a constant proportion of alkali metals. The physico-chemical problems raised another question about the supplying strategy to adopt. Again with regard to a ternary salt NaCl–MgCl₂–PuCl₃, it is possible to consider either supplying in elements (Na, Mg, Pu and Cl, see Tab. 3c) or supplying in species (NaCl, MgCl₂ and PuCl₃, see Tab. 3b). With the addition of a fourth supply, a fourth constraint is considered: ensuring electroneutrality. Both strategies have their benefits and drawbacks as presented below. The [16]’s RAPTOr is called “one supply” simulation below. This supply is made up of pure Plutonium, see Table 3a.

4.2.2. Adaptation of the SAMOSAFER extraction scheme

Because of the addition of Sodium, Magnesium and Chlorine in the supplies, it is necessary to adapt the SAMOSAFER extraction scheme in the three and four supplies simulations by adding those three elements in the chemical extraction, cross-hatched in Table A.1 in the Appendix. This adapted extraction scheme is called “with reinjection” because it behaves as if actinides and all elements displayed in white in Table A.1 were actually extracted and then reinjected instantaneously.

4.2.3. Results between one, three and four supplies strategies

Figure 7 shows the location of salt compositions on the ternary diagram during a 120 years evolution for one supply (orange), three supplies (red) and four supplies (green) strategies:

• One supply simulation remains in domain C for around 15 years before exiting while continuing to move away from the eutectic;

• Three supplies simulation remains close to the eutectic during all the evolution;

• Four supplies simulation remains close to the eutectic but begins to approach domain E after 60 years of evolution.

Fig. 7. Location of salt compositions on the ternary diagram, raw (light) and smoothed (dark) values.

In a more realistic duration of evolution of 60 years, three and four supplies simulations guarantee better safety regarding to the ternary composition of the salt.

Figure 8 shows the total salt mass for the simulations being considered. On the one hand, four supplies simulation is close to one supply simulation with a difference of −2.4 wt% at the end of the simulation. On the other hand, three supplies simulation differs from the other two by 23.7 wt% compared with one supply simulation. This difference can be explained by higher NaCl and MgCl₂ supplies to compensate for the PuCl₃ supply, while maintaining a composition close to that of the eutectic (Fig. 7). These variations are reflected in the density of the salt, with differences of −2.5% and 24.1% respectively.

Figure 9 shows the sum of all valencies in the salt for the three simulations. As expected, four supplies simulation maintains valency close to 0 during the whole evolution. In one supply simulation, the valency becomes positive, which can be interpreted as a lack of Chlorine atoms–which have a negative valency of −1. This can be explained by the constant supply of Plutonium–valency +3–instead of PuCl₃. This effect is amplified in three supplies simulations–with a value 9 times higher than in one supply simulation–but this time the explanation tends to be that too much Chlorine is extracted. The lack of experiments and data means that it is not possible to make this model “more realistic” without making too many assumptions at the moment on valence values.

Fig. 8. Mass of total salt, raw (light) and smoothed (dark) values.

Fig. 9. Valency of salt elements.

To sum up, although three supplies simulation seems more “physically realistic”, since it simulates a supply in salts rather than element by element, which is of interest in the study of deployment scenarios of such reactors, four supplies simulation has fewer drawbacks. One supply and three supplies simulations would cause a lack of Chlorine in salt which would make the salt very reductive in a real physical reactor: in reality, the elements would spontaneously change valency to match the right number of Chlorine atoms. In addition, three supplies simulation overestimates the mass balances, which can be an issue for scenario studies. Four supplies simulation seems to be the most satisfactory strategy in terms of the assumptions considered.

4.3. Impact of the extraction scheme

Figure 10 shows the total salt mass and the sum of the valences of the salt elements for a simulation of another RAPTOr-like MSR with three supplies (in black) and four supplies (in red). But in these simulations, a fraction of the whole fuel salt volume is continuously extracted without reinjecting any material: the elements displayed in white in Table A.1 are considered in blue.

Figure 10 on the top shows that this extraction scheme causes an excess of Chlorine with the three supplies strategy, which would make the salt very oxidising. But Figure 10 bottom shows that the total mass of salt is larger with four supplies simulation, mainly because the amount of Sodium and Magnesium is higher than in the three supplies simulation to match the excessive amounts of Chlorine. The four supplies strategy respects electroneutrality–modulo the assumptions made about the valences of the elements–but could overestimate the mass balances with this extraction scheme. Finally, this extraction scheme raises some question about the management of MSR spent fuel. It is more than likely that in a real MSR, a fraction of the fuel salt volume would be extracted and there will be a question of what will be reinjected. If the extracted actinides are not reinjected, the burner aim of the concept is compromised, since some actinides will end up as waste. This raises questions about the constraint on the rate of actinide reprocessing for reinjection.

Fig. 10. Valency of salt elements (top) and mass of total salt (bottom).

5. Conclusions and perspectives

A new constraint-driven supply management model has been implemented in the REM code. This model improves the algorithm by stabilising material supplies, while leaving material balances essentially unchanged compared with the previous model (differences remain within a few percent over 200 years of evolution). This improved stability makes it possible to introduce additional constraints, which are relevant for the study of a wide range of MSR concepts, including burner ones. Some of these constraints reflect specific physical and chemical requirements, such as maintaining the fraction of certain element families (alkalis and alkaline earths), controlling and matching the salt volume (through supply management or by adding a dynamic expansion vessel), and preserving the electroneutrality of the salt by accounting for element valences.

These new controls, and the flexibility they provide, enable a more detailed assessment of supply strategies during MSR operation. The results indicate that supply strategies based on chemical species yield less realistic outcomes than element-based approaches, making the latter more plausible in the context of a hypothetical real reactor. Finally, supply strategies are closely interrelated with extraction strategies. A refined management of extraction processes–whether by making them dynamic or by considering species rather than elements–may prove essential in advancing MSR depletion simulations. More broadly, these developments highlight the importance of coupling chemical constraints with reactor physics, paving the way for more realistic and operationally relevant MSR models.

Glossary

AMR: Advanced Modular Reactor

CANDU: CANadian Deuterium Uranium

FP: Fission Product

MCNP: Monte-Carlo N-Particles transport code

MSFR: Molten Salt Fast Reactor

MSR: Molten Salt Reactor

RAPTOr: Réacteur Alimenté en Plutonium et autres Transuraniens pour Orano (French)

REM: Reactor Evolution with MCNP

RK4: Runge-Kunta order 4

Funding

This research did not receive any specific funding.

Conflicts of interest

No potential conflict of interest was reported by the authors.

Data availability statement

This article has no associated data.

Author contribution statement

L. Clot: Methodology, Software, Calculations, Writing – original draft. E. Merle: Supervision, Methodology, Validation, Project administration, Writing – review and editing. A. Laureau: Supervision, Methodology, Software, Writing – review and editing. D. Heuer: Software development, Methodology, Validation. L. Tillard: Supervision, Methodology, Writing – review and editing. G. Senentz: Supervision, Project administration, Writing – review.

References

Appendix

Fig. A.1. Calculated liquidus surface of the NaCl–MgCl2–PuCl3 system with AnCl3 hypothesis and the eutectic highlighted [17].

All Tables

All Figures

	Fig. 1. Principle of neutronic evolution calculation [1].
In the text

	Fig. 2. Principle of the organisation between MCNP steps and evolution steps in the REM code.
In the text

	Fig. 3. Illustration of the organisation in Figure 2 for the first two evolution steps.
In the text

	Fig. 4. Details of the different steps during an MCNP step on a Chlorine supply (RAPTOr-like concept).
In the text

Fig. 5. Supply and cumulative supply of 233U and reactivity in the MSFR core. (a) Supply of 233U in the MSFR during 200 years of evolution. (b) Supply of 233U in the MSFR during 2 years of evolution. (c) Cumulative supply of 233U during 200 years of evolution. (d) Cumulative supply of 233U during 2 years of evolution. (e) Reactivity in the MSFR core during 200 years of evolution. (f) Reactivity in the MSFR core during 2 years of evolution.

In the text

	Fig. 6. MCNP visualisation of Cl-MSFR [3] geometry at the beginning of the evolution (left) and after 50 years (right).
In the text

	Fig. 7. Location of salt compositions on the ternary diagram, raw (light) and smoothed (dark) values.
In the text

	Fig. 8. Mass of total salt, raw (light) and smoothed (dark) values.
In the text

	Fig. 9. Valency of salt elements.
In the text

	Fig. 10. Valency of salt elements (top) and mass of total salt (bottom).
In the text

	Fig. A.1. Calculated liquidus surface of the NaCl–MgCl2–PuCl3 system with AnCl3 hypothesis and the eutectic highlighted [17].
In the text