EPJ Nuclear Sci. Technol.
Volume 5, 2019
Progress in the Science and Technology of Nuclear Reactors using Molten Salts
Article Number 14
Number of page(s) 14
Section Design of reactor
Published online 15 November 2019

© F. Di Lecce et al. published by EDP Sciences, 2019

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

Molten Salt Fast Reactor (MSFR) is the reference liquid-fuelled reactor concept in the frame of the Generation IV International Forum (GIV) [1]. Main fast spectrum liquid-fuelled reactor concepts are under investigation nowadays: the European Molten Salt Fast Reactor, the Russian MOlten Salt Actinide Recycler and Trasmuter (MOSART [2]) and other concepts worldwide (Terrapower MCFR, Elysium MCSFR, Indian Molten Salt Breeder Reactor, Moltex Energy Stable Salt Reactor). In this paper, the first one is considered as reference, which was studied in the frame of the EVOL (Evaluation and Viability of Liquid fuel) Euratom project and is currently being analyised within the SAMOFAR (Safety Assessment of the Molten Salt Fast Reactor) European H2020 project [3].1

The main innovative feature of molten salt reactors consists in the liquid state of the nuclear fuel. The immediate benefits that liquid fuels entail are the strong negative reactivity temperature feedback [4], the versatility in terms of composition and the possibility of reconfiguration of the fuel geometry. Specifically, the latter feature implies the opportunity of a new fully passive safety systems driven by the gravitational force, called Emergency Draining System (EDS) in the frame of SAMOFAR. The EDS could be triggered actively by operator or passively by salt plugs that melt when the temperature reaches a critical value (1755 °C) [5].

The safety assessments of the MSFR requires the analysis of accidental conditions that may occur during the reactor operation. A possible initiating event for the draining of the salt is the unintentional injection of fissile material that may bring the reactor in supercriticality conditions. In this situation, the plugs may be opened, to let the salt leaving the core and to reconfigure the system into a subcritical conditions. After the draining of the salt, a proper cooling of the fuel has to be envisaged and maintained to avoid an undesired increase in the fuel temperature inside the EDS that could damage the structure of the system.

Therefore, the accidental scenario of interest for the present analysis considers the possibility to drain the fuel out of the core cavity into a tank placed underneath it, where subcritical conditions are ensured and a proper cooling system is acting. At t = 0, the system is supercritical and the plugs open. In such a case, during the first phase of the draining process, the free surface of the fuel core volume starts to drop while the power production within the salt may still continue to increase, depending on the time evolution of the reactivity of the system. The geometry of the multiplying domain changes and this has a strong impact on neutron leakage and hence on the system reactivity. As multiplying domain geometry is intended the volume and the configuration of the fuel salt.

The criticality condition is also strongly depending on the temperature, due to the high temperature reactivity feedback coefficient of the fuel salt. Eventually, in the prosecution of the draining process, a subcritical configuration for the system is achieved, to be then maintained on longer time scales in the draining tank.

The objective of the present work is to investigate the time-dependent relation between the system reactivity and the salt temperature during the draining transient, considering the variation of the multiplying domain geometry and predict the impact on the temperature evolution. This latter aspect has a fundamental implication from a safety standpoint given the purpose of assuring the internal wall surface safety integrity during the transient. To this aim, the draining transient of a molten salt fuel is modelled with a multiphysics framework, in order to capture the more relevant physical phenomena, to analyse separately the effects of the different physics on the whole scenario and their interactions. A simplified semi-analytic zero-dimensional (0-D) mathematical model is proposed. Such integral lumped approach is able to capture the general dynamics of variable transients and to separate the main features of the phenomena. Details regarding the temperature and the reactivity evolution and their coupling are finally given in order to deduce preliminary conclusions on the reactor safety.

2 Moltensalt reactor core geometrical model

The geometrical configuration of the liquid fuel circuit considered in this work is a simplified version of the real geometry, similar to what performed in Wang et al. [5]. The system is represented by a cylinder with a hole in the center of the bottom surface and followed by a long pipe, to simulate the draining shaft (see Fig. 1a). The reference MSFR conceptual design [1], analyzed in the frame of the EVOL and SAMOFAR projects, is adopted to define the dimensions of this simplified geometrical domain (Tab. 1).

It is assumed that at the beginning of the transient the molten salt fills completely the core cavity and, during the draining transient, starts emptying it, flowing through the draining shaft of length L and exiting from the bottom outlet section of diameter d. The molten salt level is monitored by the quantity h(t) (the monitor length), defined as the distance of the salt free surface from the core cavity upper boundary. h(t) is set to zero at t = 0 (full core cavity) and reaches the value H when the cavity is emptied.

The thermo-physical properties of the molten salt are computed at the nominal operating temperature, which is 700 °C [1]. This comes from the decision to use a simplified modelling approach to understand the physics of the fuel draining to avoid temperature-related non-linearity that can be problematic for the semi analytical approach used in the paper. In addition, this choice is conservative from the temperature estimation point of view since the total heat capacity of the salt is an increasing function of the temperature given the correlation stated in [1, 6]. Tables 2 and 3 show the properties used in this paper.

thumbnail Fig. 1

Simplified geometry of the fuel circuit for the zero dimensional model. (a) The molten salt evolution is represented by the quantity h(t), i.e. the distance of the salt free surface from the cavity upper surface. (b) A time-varying control volume (in red) is adopted to describe the draining phenomenon.

Table 1

Model domain dimensions [5].

Table 2

LiF-ThF4 fuel salt reference thermo-physical properties.

Table 3

Neutronic parameters.

3 The multiphysics draining model

The core cavity draining process is described by a 0-D semi-analytical model. The lumped integral approach to the draining phenomenon allows capturing and comprehending the dynamics of main variables. In addition, thanks to the flexibility and the low computational resource of the zero dimensional approach, the model is able to deduce relevant model features by means of parametric studies. The fuel draining is a multiphysics transient problem. The energy balance require the information of the fission power and the decay heat which are an outcomes of the neutronics analysis. In turn equations for neutrons and precursors require the temperature in order to model the temperature dependency of cross sections. Moreover, the multiplying domain, i.e. the molten salt in the core cavity, changes in time as the liquid is drained. Hence, the system volume is strongly dependent on the salt fluid-dynamics. Due to the negative thermal feedback, both the temperatureand the geometry affect the system reactivity negatively, and therefore on the heat production. Salt thermo-physical properties are in principle functions of temperature, which implies a coupling between fluid-dynamics and energy balance. However this coupling is neglected since salt properties are kept temperature independent. The multiphysics interactions among the different phenomena can be graphically appreciated in Figure 2.

The model aims at capturing the multiphysics characteristics, highlighting the temperature and the system reactivity evolutions in time. The presented model can be seen as the coupling of three sub-models:

  • the thermal-hydraulics sub-model;

  • the neutronics sub-model;

  • the molten salt level evolution sub-model.

Each sub-model is described in details in the following subsections, presenting the assumptions made and the strategies adopted for its solution.

thumbnail Fig. 2

Graphical scheme depicting the multiphysics interaction of the molten salt fuel draining phenomenon. The dashed arrow identifies the coupling between fluid-dynamics and energy balance through the temperature-dependent salt properties, currently not modelled.

3.1 Thermal-hydraulics sub-model

When the draining of the molten salt starts, the lowering of the free surface level leaves an empty space on the top of the core cavity, which is assumed to be filled by air or some other inert gas. The approach adopted in this work considers a time-varying control volume, which coincides with the salt in the cavity. The moving upper boundary is the Control Surface A (CS-A), while the bottom outlet boundary represents the CS-B. The other boundaries are the peripheral walls (see Fig. 1).

The mass variation within the control volume is due only to the quantity of molten salt that is exiting from the CS-B [7]. Assuming an incompressible fluid and constant properties, it is possible to demonstrate that the mass equation for this system reduces to (1)

which is a balance on the volumetric flow rate at the CS-A and CS-B. The functions vA (t) and vB (t) in the equation (1) are the free surface and the outlet velocities respectively, while h(t) measures the increasing height of the gas head. The modelling of the velocities vA and vB pertains to the following sub-models. For additional information on the mathematical treatment, the reader may refer to [8].

The variation of energy contained in the control volume depends on the energy flux going out through CS-B and on the energy source , linked to neutronics (prompt fission power and decay heat). The energy equation can be written as follows (for more detail, see [8]): (2)

where T(t) is the average temperature of the control volume and TB(t) is the fuel temperature at CS-B. If the temperature was spatially homogeneous in the system (i.e., setting T(t) equal to TB(t)), then the salt outflow would not cause any variation in the average temperature. Anyway, the temperature within the system is not uniform and, in the frame of a 0-D model, this spatial feature is provided with the difference T(t) − TB(t) given as input. In particular, the position on the bottom surface of the draining shaft affects the molten salt extracted from the cavity and also its temperature.

For this reason, the temperature difference has been computed from a CFD analysis [8] and then employed to derive an analytical expression to model the T(t) − TB(t) term in time (see Fig. 3). Fixing the initial and the final temperature difference to be 50 K (which is assumed to be the initial inlet-to-core temperature difference) and 0 (the mean salt temperature coincides with the outlet temperature at the draining end) respectively, it yields: (3)

where ΔT0 is the initial value of the temperature difference, which is set to 50 K [1], t* is the draining time and λmix is a fitting time constant related to the temperature spatial heterogeneity and to the salt mixing. The value for this time parameter is 0.024 1/s.

thumbnail Fig. 3

Difference between average and outflow molten salt temperatures as a function of time. Red: CFD simulation of the drainingtransient; blue: analytical best-fit.

3.2 Neutronics sub-model

The molten salt system neutron kinetics is approximated by the point kinetics equations [10]. Defining the normalized neutron population η(t) = n(t)∕n0 and the normalized precursor concentration ξ(t) = C(t)∕C0, where n0 and C0 are the neutron and precursor concentrations at the beginning of the transient, the normalized point kinetics equations, assuming only one family of delayed neutron precursor for simplicity, can be written as: (4)

In the precursor equation, the terms relating to the precursor motions should be taken into account [11]. On the other hand, due to the specific configuration considered here (i.e., the absence of the circulation loop and the lack of a inflow contribution for the precursor), spatially uniform concentration of precursors is assumed, which means that the outflow contribution to the neutron escape is balanced by the system volume change. In other words, the outflow molten salt does not affect the concentration of precursors in the system because the number of precursors that exits through CS-B is perfectly balanced by the volume reduction due to the draining. Anyway, this is true if the outflow concentration and the mean concentration of precursors are equal (i.e. the precursors concentration is spatially uniform in the domain), which is not the real case but it is a consequence of the adoption of the zero dimensional approach.

The reactivity ρ(t) is the key parameter concerning the system neutron kinetics. It is affected by the fuel temperature through the feedback coefficient and the temporal modificationof the system volume. With reference to the present study, the system starts in a supercritical condition, i.e. a reactivity larger than zero, to mimic an accidental scenario in which the neutron population increases with a time-scale related to the neutron lifetime and the reactivity value. During the evolution of the draining transient, the system will experience a decrease of the reactivity, leading to a subcritical condition (negative reactivity).

An explicit definition of the reactivity as a function of the multiplication constant k and considering a one-group diffusion approximation can be drawn. In this regard, it is worth mentioning that, since the MSFR is a fast spectrum reactor, the migration area and the diffusion area can be considered equal [12, 13]. Recalling that the diffusionarea is equal to the diffusion coefficient divided by the absorption cross section, it is possible to express the variation of ρ as: (5)

Observing equation (5), we can observe that the reactivity variation δρ is composed by three contributions, namely (i) the variation of the multiplication factor of the infinite medium, which is strictly related to temperature variation of cross sections, (ii) the non-leakage probability variation, which is in turn composed by one term related to geometrical change (the larger the domain, the less the probability for neutrons to escape) and (iii) to another term associated to the neutron diffusion length, depending on temperature. These three terms are now analyzed in detail.

3.2.1 Variation of infinite multiplication factor

It is already mentioned that the variation of k = νΣfΣa can be produced by the temperature dependence of cross sections. The macroscopic cross section (Σ) dependence on temperature can be expressed by a logarithmic function [12], which is assumed to be well suited for fast reactor behavior: (6)

where Σx is the macroscopic cross sectionfor a generic reaction x, supposed to be affected by the Doppler broadening effects. Therefore, the first term in equation (5) can be re-written as (7)

where is defined as the temperature reactivity feedback coefficient related to the variation of k : (8)

Two effects are present in the trend, considering the values listed in Table 3. An increase of temperature causes a reduction of the fission cross section νΣf (T), represented by the negative value of , along with an increase of absorptions, represented by the positive value of . As general effect, the infinite multiplication factor decreases following a temperature increase, and therefore the temperature coefficient is negative. Figure 4a depicts the dependence of this coefficient on temperature, using the cross section values listed in Table 3.

3.2.2 Variation of buckling

The geometrical buckling of the cylindrical multiplying domain is defined as: (9)

where j0,1 is the first zero of the zero-th order Bessel function of the first kind [14]. The extrapolated dimensions, which the buckling refers to, are approximated to the real dimensions of the cylindrical cavity. Furthermore, the geometry variation is expressed by the monitor length h(t), that appears in the buckling relation. Rearranging the second term in equation (5), including equation (9) and developing the derivative with respectto h, the following result is obtained: (10)

where αh(h, T) is defined as the reactivity coefficient due to the geometry modification: (11)

Equation (11) represents the system reactivity variation due to the variation of leakages related to the volume modification. In general, neutrons escape with a higher probability from a small volume rather than in larger domains. Therefore, during the molten salt draining, the reduction of the volume implies the increase of the probability of leakage and hence the insertion of a negative reactivity. Figure 4b shows the geometry reactivity coefficient as a function of the monitor length for a fixed value of temperature of 900 K (see Tab. 3 for the data). A variation of the fuel temperature implies irrelevant variations of αh (h, T), assuming the same values for h, due to the limited dependence on temperature of the diffusion area.

3.2.3 Variation of diffusion area

The last term in equation (5) regards the change of the probability of non-leakage due to a variation of the diffusion area L2. The larger is the diffusion capabilities of the neutrons to travel in the medium, the larger is the probability to reach the boundaries and escape. The diffusion area is affected by temperature. Being L2 (T) = Dn(T)∕Σa(T), introducing the temperature relations (6) and computing the derivative with respect to temperature, it is obtained: (12)

where is the temperature reactivity feedback coefficient due to the diffusion area variation: (13)

Temperature increase implies on one hand the increase of absorptions of neutrons in the medium, represented by the positive value of but on the other the reduction of neutron diffusion coefficient, stated by the slight negative value of This causes the decrease of the diffusion area and indeed the decrease of leakage probability. Finally this results in a positive temperature reactivity coefficient due to diffusion area variations (Fig. 4c). Moreover, also the geometry plays a role in the definition of assumed in this work. In particular, the larger the domain (low values of h), the smaller the leakage probability and indeed the lower is the coefficient .

On the basis of the previous reasoning, a total temperature feedback coefficient can be defined, as the sum of the temperature coefficients aforementioned: (14)

It is interesting to notice that the absorption cross section has two opposite impacts on the reactivity. On one hand, the Doppler effect causes the increase in the absorptions and indeed a negative effect on reactivity. On the other hand, the absorption increase causes a reduction of the leakage rate and indeed a positive effect on the neutron economics. However, observing Figure 4d, it can be noticed that the contribution of diffusion area variation to the total temperature reactivity coefficient is smaller than the k variation contribution. At last, the total differential variation of reactivity can be written as: (15)

The temperature and geometry differential variations, weighted on the respective reactivity coefficients, influence the differential reactivity variation. To obtain the integral variation of reactivity, both sides of equation (15) have to be integrated: (16)

where the reference conditions for temperature and core cavity height are Tref and href, respectively, and it is assumed that the reference reactivity ρ(href, Tref) is null, i.e. criticality condition holds for the reference temperature and monitor length. ρT (h, T) and ρh (h, T) are the integral reactivity variations due to temperature and geometry, respectively. The reference temperature Tref is imposed equal to 900 K [9] and the corresponding reference monitor length href is obtained imposing criticality condition to the system.

thumbnail Fig. 4

Temperature feedback coefficients as function of temperature. The coefficient depends alsoon the monitor length h.

3.3 Heat source from fissions and decay heat

The power source term in equation (2), , is related to neutronics and contains two contributions, i.e., the instantaneous deposition of heat due to fissions and the decay heat. When a fission event occurs, the total energy generated is around 200 MeV. However, only a fraction of this energy amount is instantaneously deposited in the medium. The remaining part, associated to the fission products, is emitted at later times through γ, β and neutron decay emissions. The total fission volumetric power at a certain time instant is expressed by the following formula (17)

where EfΣf(T) is the energy per fission event produced times the fission cross section and v the neutron speed. In the frame of the present 0-D model, is a function time only.

As already said, the large part of the power is directly deposited in the medium, while a fraction f is stored into decay heat precursors and emitted at later times. A balance for the decay heat precursors in terms of the energy stored, where the source is represented by the fraction f of fission energy produced and the sink is the decay event, can be written as: (18)

where qd(t) is the decay energy stored in the precursors per unit volume and λd is the decay event time constant. One group of decay heat precursors is considered, without loss of generality. The total volumetric thermal power source is then obtained by summing the instantaneous fission energy production and the decay heat and multiplying them by the molten salt volume, which is time-dependent: (19)

3.4 Molten salt level sub-model

The dynamics of the fuel salt draining depends strongly on the free surface velocity vA (t) (1), i.e. on the time derivative of the monitor length h. It identifies the time evolution of the free surface of the salt that is leaving the cavity and affects the energy content in the system and mostly the neutronic dynamics. Therefore, the modeling of d h∕dt constitutes the closure to our problem.

Isothermal draining phenomena of incompressible fluids are old and well-established problems. An analytical solution is available in the current literature [15]. However, the molten salt fuel draining presents strong temperature variations during the transient and in principle it is not possible to treat it as isothermal. As a first approximation in the present study, the approach proposed in [15] for isothermal draining phenomena is adopted. This assumption is suitable to analyse the reactivity and the temperaturedraining transients.

The method assumes a quasi-steady set of equations, i.e., an unsteady mass balance in addition to a steady energy balance (the Bernoulli’s law). The main result is the expression of the outlet velocity as a function of time (for more details, please refer to [8]): (20)

where θ is a dimensionless constant that describes the friction losses along the core cavity and the shaft. The mathematical expression for θ is (21)

where Cp is the Fanning friction factor along the shaft (which is a function of the Reynolds number, and indeed, of the salt viscosity) and Kc the local friction factor for the cavity-to-shaft sudden contraction, which is assumed to be 0.5 [1618].

We now reconsider equation (1), expressing the unsteady mass balance of the salt draining. Substituting the relation for vB(t) from equation (20), one obtains: (22)

An Ordinary Differential Equation (ODE) for the variable h is obtained with initial condition h(0) = 0, i.e. the cavity full of salt. The solution of such equation is: (23)

Equation (23) represents the closure to the present problem. An analytical expression for the salt free surface time evolution allows for the equations governing the temperature and the neutronics to be fully defined and numerically solvable.

The time evolution of the molten salt level in the core is plotted in Figure 5. The draining time corresponds to the time period needed to completely empty the core cavity. The analytical expression can be derived from equation (23) imposing h(t*) = H and solving for t*. Imposing atmospheric pressures both on the free surface and the bottom pressure, the draining time results being 94.82 s, which is in the range expected in [5]. This is an indicative result for MSFRs, since it is not clear which are the pressure conditions to be expected at the draining valves. However, the draining time estimated with this approach is consistent with other literature results [5] and, moreover, CFD simulations confirm the suitability of the approach adopted even for strongly non-isothermal transients [8].

thumbnail Fig. 5

Molten salt level evolution in time. Initial value H = 2.255 m (full cavity).

3.5 Initial conditions for the coupled simulations

Once the function h(t) is known, the coupled equations describing the draining can be solved, provided initial conditions are defined. As for the reference MSFR operating parameters [1], the average salt temperature is set to 700 °C and the initial power is equal to 3.5 GW, which might be the onset reactor power of an accidental scenario. The initial conditions for the draining problem are summarized in Table 4.

Table 4

Initial conditions for the coupled calculations.

4 Simulation results

The governing equations modelling the draining of the molten salt fuel present a strong coupling among them and a high non-linearity. Therefore, such equation system needs to be solved numerically. Although all the ODEs are solved with appropriate numerical schemes, some analytical formulae are provided as input to the solution process, as the analytical expressions for the salt free surface time evolution (Eq. (1)) and for the mean outflow temperaturedifference (Eq. (3)). This aspect explains the semi-analytical characteristics of the model presented.

Since the mathematical model is stiff (i.e., highly different time scales that range from the neutron life time to the draining time scale are present), a suitable solver must be identified to treat the set of ODEs. As a result, the equation system is solved with a variable-order solver based on the numerical differentiation formulas (NDFs), which are related to the Gear’s method [19, 20].

At t = 0, the cavity is full of molten salt, the system temperature is 700 °C (973 K). The accidental conditions are represented by an initial core thermal power corresponds of 3.5 GW and an initial reactivity of 285 pcm (0.85 $), i.e, corresponding to a supercritical system. Then, the salt draining is started. Figure 6 shows the system reactivity evolution during the first 7 seconds of the transient. The ρh (h, T) and ρT (h, T) terms are depicted, as described in equation (16), as well as their sum producing the total system reactivity. The reactivity contribution due to temperature is negative as a result of temperature feedback, asthe system starts from an initial value larger than the reference temperature Tref = 900 K (Tab. 3). Conversely, the molten salt volume is larger than the critical one (href = 0.1768 m), producing an initial supercritical phase. The system stays in supercritical conditions for about 1 second, with a strong impact on the system dynamics. In particular, the initial supercritical phase induces a sudden increase in the neutron density in the system (see Fig. 7a), which increases from its initial value of 1.33× 1013 n/cm3 to about 8.3 × 1013, with the time scale of the effective neutron lifetime Λ. The increased neutron production causes both the increase of neutron precursors (see Fig. 7b) and decay heat precursors (see Fig. 7d). The system temperature increases due to the production of heat coming from the instantaneous deposition from fission events (prompt volumetric heat, which is the dominant contribution) and from the decay heat. Figures 7c and 7d show the prompt and the decay volumetric heat time evolutions. The first one follows the trend of the neutron population, since it is proportional to η(t) according to equation (17). Then, the temperature impacts strongly on the reactivity due to the negative feedback coefficient. Furthermore, the lowering of the salt free surface level causes the reduction of the multiplying system volume, the increase of neutron leakages hence adding an another negative effect on neutron economics. Therefore the reactivity goes down and the system reaches the critical state within 1 s, then it goes towards more and more negative values. The neutron population suddenly drops and the precursors approach their exponential decay form (see Figs. 7b and 7d).

The temperature time evolution is reported in Figure 7f. The first part of the transient is characterized by the temperatureramp due to power production within the salt. The mean temperature increases by 200 K, reaching a peak of 1165 K (892 °C). This value is far below the temperature that may cause damages to core inner surfaces, which is approximately 1600 K for Nickel-based walls [6]. After this ramp, the molten salt outflow enthalpy overcomes the power source term and the temperaturestarts decreasing to reach 1025 K at the transient end. It should be noticed that the time derivative of the temperaturedepends on three terms, as can be observed on the RHS of equation (2). The second term, related to the outflow of kinetic energy, is negligible during the whole transient. The heat source, involved in the last term, is positive and causes the temperature to rise up during the first part of transient. However, the contribution of this term decreases in time, since the total heat rapidly decreases (see Fig. 7e), and in the meanwhile the first term depending on the temperature difference between the mean value T(t) and the outflow value TB(t) becomes dominant. After a while, the first negative term overcomes the positive term related to heat production and the temperature begins to reduce. The heat source term is dominant for the first 20 s and then it is overcome by the term related to the molten salt outflow. As regards the safety aspects of this transient, simulations show that the strong temperature feedback coefficient (around −4 pcm/K) causes the immediate reactor shutdown in case of an hazardous temperature (or power) increase.

thumbnail Fig. 6

System reactivity as a function of time during the first 7 seconds of transient. The red, the green and the blue lines represent the total ρ(h, T), the geometry-related ρh(h, T) and the temperature-related ρT(h, T) system reactivities, respectively.

thumbnail Fig. 7

Simulation results of the 0-D semi-analytical model with input parameters as in Table 4.

4.1 Uncertainty analysis of main parameters

An uncertainty analysis has been carried out to understand and quantify the effects on the solution, in terms of temperature, reactivity and heat production, associated to variations of thermo-physical and some neutronics parameters of the molten salt fuel. The parametric analysis is performed with a one-factor-at-time method, varying a parameter at a time while keeping others fixed [21]. An uncertainty of ±10 % on thermo-physical and neutronics parameters is considered and the effects on the model solution are evaluated.

Focusing on the effects related to the fuel mass density, it is observed that an increase in the density causes a decrease in the source term of energy equation (2), as can be seen in Figures 8a and 8b. In particular, higher mass densities implies a lower temperature transient, which in turn causes a smaller temperature reactivity feedback. This results in an increase, even if limited, of the neutron population and of the power production. The viscosity of the fluoride salt, within the developed model, has an impact on the calculation of the friction-related parameter θ. The temperaturevariations caused by the viscosity uncertainties is negligible, whereas it implies a slight variation on the draining period since the larger viscosity, the larger the draining time to empty the cavity. The uncertainty on the specific heat (see Figs. 8c and 8d) has the same effects on the solution as the mass density, as can be understood looking at the energy equation.

A similar analysis has been performed on some neutronics parameters. The effective mean generation time Λ determines the slope of the neutron jump up to the peak. The neutron precursor decay constant, λp, has an impact on the transient tail. Lower decay constants mean an increase in the delayed neutron production and indeed in the power generation. This hence results in a higher temperature in the transient tail (Fig. 8). The same holds for the decay constant of decay heat precursors, where its evolution is reported in Figure 8f. The delayed neutron fraction β has a strong impact on neutronics, due to its influence on the neutron dynamics (Fig. 9b). In addition, it has a direct impact on the neutron prompt jump. Larger delayed neutron fractions enhance the neutron precursor importance and indeed imply smaller prompt neutron jumps (Fig. 9c). The temperature is anyway affected by less than the 0.6 % (Fig. 9a). Lastly, the decay heat precursor fraction, λd, has a negligible impact on the solution.

The results of this preliminary uncertainty analysis are summarised in Table 5. In this table, only the maximum variation on the temperature along the transient in correspondence to a ±10 % variation of an input variable is reported.

thumbnail Fig. 8

Uncertainty analysis. Transients of temperature and other variables of interest, where the nominal problem solution (red line) is compared with variation of parameters of ±10% (dash-dot lines).

thumbnail Fig. 9

Uncertainty analysis. Transients of temperature and other variables of interest, where the nominal problem solution (red line) is compared with variation of parameters of ±10% (dash-dot lines).

Table 5

Maximum relative variation of the temperature transient for a parameter shift of ±10% from its nominal value. The center column shows temperature maximum change for an input parameter increase of 10%, while in the right column the parameter is decreased of the same percentage.

5 Concluding remarks

The emergency core draining system is a fully-passive safety system concept developed for liquid-fueled molten salt reactors. In case of accidental situations one or more salt plugs, located at the core bottom, open and the fuel salt starts being drained, thanks to gravitational force, and stored in a draining tank underneath the core.

A simplified zero-dimensional semi-analytical model is proposed to capture the multiphysics phenomena involved in the fuel draining. Point kinetics equations are used to model the neutronics and are coupled with the energy equation. The time dependent behavior of the salt level is given with an analytical expression, adopting a quasi-steady approach, to close the problem. The 0-D model manages to describe the general dynamics of variables. The impact on the transient of the different physical phenomena can be analysed separately, understanding which are the most affecting parameters and to fully comprehend the transient dynamics. The reactivity decreases in time is due to different phenomena, namely the increase of the temperature– since the salt has a strongly negative feedback coefficient –, and the volume reduction – meaning an increased neutron leakages. The fuel salt draining gives the opportunity to study a volume-change multiplying domain problem and to derive a general analytical model to describe the reactivity feedback due to geometry variations. The analysis of the simulation results show that the temperature increases for the first 20 seconds of the transient, due to the thermal power production related to fissions and decay heat. From the initial value of 700 °C, it reaches its peak of about 900 °C, staying below the critical value to produce damages to the structures. After this ramp, the temperature starts decreasing due to the outflow enthalpy stream. The system reactivity starts from an imposed value of 285 pcm to simulate an accidental scenario, then it reaches criticality in about 1 second and then it becomes more and more subcritical. The simulation proves the intrinsic stability of the molten salt fuel thanks to the highly negative temperature feedback coefficient and the contribution on reactivity of lowering of the salt level.

Further analyses are ongoing in order to develop a CFD model of the multiphysics draining problem that may relax some of the modelling hypothesis of the semi-analytical model developed in this paper. A critical comparison of the two approaches will be the object of a future work.

Author contribution statement

The study presented in this article has been developed in the frame of Francesco Di Lecce's master thesis, which has been carried out under the supervision of Stefano Lorenzi and Antonio Cammi from PoliMi and Sandra Dulla and Piero Ravetto from PoliTo. The mathematical model and computational calculations have been developed by Francesco Di Lecce. All the authors have contributed to this work with critical and expert judgement and by providing support and proofreading. Finally, Francesco Di Lecce has mainly written the article, together with the contribution of Stefano Lorenzi and of all authors.


: Macroscopic cross section at Tref for Σx (1/m)

: Temperature reactivity feedback coefficient related to variation of k (1/K)

: Temperature reactivity feedback coefficient due to L2 (1/K)

αh: Reactivity coefficient due to geometry change (1/m)

: Doppler coefficient for Σx (1/m)

: Doppler coefficient for Dn (m)

β: Total effective delayed neutron fraction (1)

δ: Salt mass density (kg/m3)

η: Normalized neutron population (1)

Λ: Effective neutron generation lifetime (s)

λd: Total decay constant of decay heat precursors (1/s)

λp: Total decay constant of neutron precursors(1/s)

λmix: Time constant of salt mixing (s)

ν: Kinematic viscosity (m2/s)

ρ: System reactivity (1)

ρh: Reactivity variations due to geometry (1)

ρT: Reactivity variations due to temperature (1)

Σx: Macroscopic cross section of a reaction x (1/m)

ξ: Normalized precursor concentration (1)

: Source of thermal power (W)

: Total fission volumetric power (W/m3)

v: Neutron speed (m/s)

B2: Geometrical buckling of the system (m−2)

C0: Precursor concentration at equilibrium (1/m3)

cp: Specific heat at constant pressure (J/kg/K)

d: Shaft diameter (m)

Dn: Diffusion coefficient (m)

f: Total fraction of decay heat precursors (1/s)

H: Cavity height (m)

h: Monitor length (m)

k: Effective multiplication factor (1)

k: Infinite multiplication factor (1)

L: Shaft length (m)

L2: Diffusion area (m2)

n0: Neutron concentration at equilibrium (1/m3)

qd: Decay energy stored in precursors (J/m3)

R: Cavity radius (m)

T: Temperature (K)

t: Time (s)

t*: Draining time (s)

TB: Outlet salt temperature (K)

vA: Salt free surface velocity (m/s)

vB: Salt outflow velocity (m/s)


  1. M. Allibert, M. Aufiero, M. Brovchenko, S. Delpech, V. Ghetta, D. Heuer, A. Laureau, E. Merle-Lucotte, 7-Molten salt fast reactors, in Woodhead Publishing Series in Energy, Handbook of Generation IV Nuclear Reactors, edited by I.L. Pioro (Woodhead Publishing, 2016) [Google Scholar]
  2. V. Ignatiev, O. Feynberg, I. Gnidoi, A. Merzlyakov, A. Surenkov, V. Uglov, A. Zagnitko, V. Subbotin, I. Sannikov, A. Toropov, V. Afonichkin, A. Bovet, V. Khokhlov, V. Shishkin, M. Kormilitsyn, A. Lizin, Molten salt actinide recycler and transforming system without and with Th-U support: Fuel cycle flexibility and key material properties, Ann. Nucl. Energy 64, 408 (2014) [CrossRef] [Google Scholar]
  3. M. Brovchenko, E. Merle-Lucotte, H. Rouch, F. Alcaro, M. Allibert, M. Aufiero, A. Cammi, S. Dulla, O. Feynberg, L. Frima, O. Geoffroy, V. Ghetta, D. Heuer, V. Ignatiev, J.L. Kloosterman, D. Lathouwers, A. Laureau, L. Luzzi, B. Merk, P. Ravetto, A. Rineiski, P. Rubiolo, L. Rui, M. Szieberth, S. Wang, B. Yamaji, Optimization of the pre-conceptual design of the MSFR, Work-Package WP2, Deliverable D2.2, EVOL, Evaluation and Viability of Liquid fuel fast reactor system) European FP7 project, Contract number: 249696, EVOL (2014) [Google Scholar]
  4. J. Serp, M. Allibert, O. Benes, S. Delpech, O. Feynberg, V. Ghetta, D. Heuer, D. Holcomb, V. Ignatiev, J.L. Kloostermang, L. Luzzi, E. Merle-Lucotte, J. Uhlír, R. Yoshioka, D. Zhimin, The molten salt reactor (MSR) in generation IV: Overview and perspectives, Progr. Nucl. Energy 77, 308 (2014) [Google Scholar]
  5. S. Wang, M. Massone, A. Rineiski, E. Merle-Lucotte, Analytical Investigation of the Draining System for a Molten Salt Fast Reactor, in The 11th International Topical Meeting on Nuclear Rector Thermal Hydraulics, Operation and Safety, Korea, 2016 [Google Scholar]
  6. V. Ignatiev, O. Feynberg, A. Merzlyakov, A. Surenkov, A. Zagnitko, V. Afonichkin, A. Bovet, V. Khokhlov, V. Subbotin, R. Fazilov, M. Gordeev, A. Panov, A. Toropov, Progress in development of MOSART concept with Th support, in International Congress on Advances in Nuclear Power Plants 2012, 2012 [Google Scholar]
  7. Y.A. Çengel, J.M. Cimbala, Fluid Mechanics: Fundamentals and applications (McGraw-Hill Higher Education, 2004) [Google Scholar]
  8. F. Di Lecce, Neutronic and thermal-hydraulic simulations for Molten Salt Fast Reactor safety assessment (unpublished master thesis), Politecnico di Torino, Torino, Italia, 2018 [Google Scholar]
  9. M. Aufiero, Development of advanced simulation tools for circulating fuel nuclear reactors, Doctoral dissertation, 2014 [Google Scholar]
  10. A.E. Waltar, A.B. Reynolds, Fast Breeder Reactors (Pergamon Press, 1981) [Google Scholar]
  11. D.L. Hetrick, Dynamics of Nuclear Reactors (The University of Chicago Press, 1971) [Google Scholar]
  12. A.E. Waltar, D.R. Todd, P.V. Tsvetkov, Fast Spectrum Reactors (Springer, 2012) [CrossRef] [Google Scholar]
  13. E.L. Lewis, Fundamentals of Nuclear Reactor Physics (Elsevier Science Publishing Co Inc., 2008) [Google Scholar]
  14. M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions (United States Department of Commerce, 1964) [Google Scholar]
  15. D.D. Joye, B.C. Barrett, The tank drainage problem revisited: do these equations actually work? Can. J. Chem. Eng. 81, 1052 (2003) [CrossRef] [Google Scholar]
  16. T.L. Bergman, A.S. Lavine, F.P. Incropera, D.P. DeWitt, Introduction to Heat Transfer, 6th edn. (Wiley, 2011) [Google Scholar]
  17. F.A. Morrison, An Introduction to Fluid Mechanics (Cambridge University Press, 2013) [CrossRef] [Google Scholar]
  18. Crane, Flow of Fluids through Valves, Fittings, and Pipe (Crane Co. Engineering Division, 1957) [Google Scholar]
  19. MATLAB and Statistics Toolbox Release R2017a, The MathWorks, Inc., Natick, Massachusetts, United States [Google Scholar]
  20. G.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations (Prentice-Hall, Englewood Cliffs, N.J., 1971) [Google Scholar]
  21. V. Czitrom, One-Factor-at-a-Time Versus Designed Experiments (American Statistician, 1999) [Google Scholar]


A Paradigm Shift in Nuclear Reactor Safety with the Molten Salt Fast Reactor, Grant Agreement number: 661891 | SAMOFAR, Euratom research and training programme (2014–2018).

Cite this article as: Francesco Di Lecce, Antonio Cammi, Sandra Dulla, Stefano Lorenzi, and Piero Ravetto. Simplified 0-D semi-analytical model for fuel draining in molten salt reactors, EPJ Nuclear Sci. Technol. 5, 14 (2019)

All Tables

Table 1

Model domain dimensions [5].

Table 2

LiF-ThF4 fuel salt reference thermo-physical properties.

Table 3

Neutronic parameters.

Table 4

Initial conditions for the coupled calculations.

Table 5

Maximum relative variation of the temperature transient for a parameter shift of ±10% from its nominal value. The center column shows temperature maximum change for an input parameter increase of 10%, while in the right column the parameter is decreased of the same percentage.

All Figures

thumbnail Fig. 1

Simplified geometry of the fuel circuit for the zero dimensional model. (a) The molten salt evolution is represented by the quantity h(t), i.e. the distance of the salt free surface from the cavity upper surface. (b) A time-varying control volume (in red) is adopted to describe the draining phenomenon.

In the text
thumbnail Fig. 2

Graphical scheme depicting the multiphysics interaction of the molten salt fuel draining phenomenon. The dashed arrow identifies the coupling between fluid-dynamics and energy balance through the temperature-dependent salt properties, currently not modelled.

In the text
thumbnail Fig. 3

Difference between average and outflow molten salt temperatures as a function of time. Red: CFD simulation of the drainingtransient; blue: analytical best-fit.

In the text
thumbnail Fig. 4

Temperature feedback coefficients as function of temperature. The coefficient depends alsoon the monitor length h.

In the text
thumbnail Fig. 5

Molten salt level evolution in time. Initial value H = 2.255 m (full cavity).

In the text
thumbnail Fig. 6

System reactivity as a function of time during the first 7 seconds of transient. The red, the green and the blue lines represent the total ρ(h, T), the geometry-related ρh(h, T) and the temperature-related ρT(h, T) system reactivities, respectively.

In the text
thumbnail Fig. 7

Simulation results of the 0-D semi-analytical model with input parameters as in Table 4.

In the text
thumbnail Fig. 8

Uncertainty analysis. Transients of temperature and other variables of interest, where the nominal problem solution (red line) is compared with variation of parameters of ±10% (dash-dot lines).

In the text
thumbnail Fig. 9

Uncertainty analysis. Transients of temperature and other variables of interest, where the nominal problem solution (red line) is compared with variation of parameters of ±10% (dash-dot lines).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.