Issue 
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  
DOI  https://doi.org/10.1051/epjn/2019028  
Published online  15 November 2019 
https://doi.org/10.1051/epjn/2019028
Regular Article
Simplified 0D semianalytical model for fuel draining in molten salt reactors
^{1}
Politecnico di Torino, Dipartimento Energia, NEMO group,
Torino,
Italy
^{2}
Politecnico di Milano, Department of Energy, Nuclear Engineering Division,
Milano,
Italy
^{*} email: antonio.cammi@polimi.it
Received:
23
April
2019
Received in final form:
15
July
2019
Accepted:
26
August
2019
Published online: 15 November 2019
A key feature of molten salt reactors is the possibility to reconfigure the fuel geometry (actively or passively driven by gravitational forces) in case of accidents. In this regard, the design of reference molten salt reactor of Generation IV International Forum, the MSFR, foresees the Emergency core Draining System (EDS). Therefore, the research and development of MSFRs move in the direction to study and investigate the dynamics of the fuel salt when it is drained in case of accidental situations. In case of emergency, the salt could be drained out from the core, actively or passively triggered by melting of salt plugs, and stored into a draining tank underneath the core. During the draining transient, it is relevant from a safety point of view that thermal and mechanical damages to core internal surfaces and to EDS structure – caused by the temperature increase due to the decay heat – are avoided. In addition, the subcriticality of the fuel salt should be granted during all the draining transients. A simplified zerodimensional semianalytical model is developed in this paper to capture the multiphysics interactions, to separate and analyse the different physical phenomena involved and to focus on time evolutions of temperature and system reactivity. Results demonstrate that the fuel draining occurs in safe conditions, both from the thermal (temperaturerelated internal surface damages) and neutronic (subcritical states dominate the transient) view points and show which are the main characteristics of the fuel salt draining transient.
© F. Di Lecce et al. published by EDP Sciences, 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), 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 liquidfuelled reactor concept in the frame of the Generation IV International Forum (GIV) [1]. Main fast spectrum liquidfuelled 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 timedependent 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 semianalytic zerodimensional (0D) 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 thermophysical 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 temperaturerelated nonlinearity 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.
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 timevarying control volume (in red) is adopted to describe the draining phenomenon. 
LiFThF_{4} fuel salt reference thermophysical properties.
Neutronic parameters.
3 The multiphysics draining model
The core cavity draining process is described by a 0D semianalytical 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 fluiddynamics. Due to the negative thermal feedback, both the temperatureand the geometry affect the system reactivity negatively, and therefore on the heat production. Salt thermophysical properties are in principle functions of temperature, which implies a coupling between fluiddynamics 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 submodels:

the thermalhydraulics submodel;

the neutronics submodel;

the molten salt level evolution submodel.
Each submodel is described in details in the following subsections, presenting the assumptions made and the strategies adopted for its solution.
Fig. 2 Graphical scheme depicting the multiphysics interaction of the molten salt fuel draining phenomenon. The dashed arrow identifies the coupling between fluiddynamics and energy balance through the temperaturedependent salt properties, currently not modelled. 
3.1 Thermalhydraulics submodel
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 timevarying control volume, which coincides with the salt in the cavity. The moving upper boundary is the Control Surface A (CSA), while the bottom outlet boundary represents the CSB. 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 CSB [7]. Assuming an incompressible fluid and constant properties, it is possible to demonstrate that the mass equation for this system reduces to $${v}_{A}(t)=\frac{\text{d}h(t)}{\text{d}t}={\left(\frac{d}{2R}\right)}^{2}{v}_{B}(t),$$(1)
which is a balance on the volumetric flow rate at the CSA and CSB. The functions v_{A} (t) and v_{B} (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 v_{A} and v_{B} pertains to the following submodels. 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 CSB and on the energy source ${\dot{Q}}_{f}(t)$, linked to neutronics (prompt fission power and decay heat). The energy equation can be written as follows (for more detail, see [8]): $$\begin{array}{l}\left(Hh(t)+{\left(\frac{d}{2R}\right)}^{2}L\right)\frac{\text{d}T(t)}{\text{d}t}\\ \text{\hspace{1em}}=\left(T(t){T}_{B}(t)\right)\frac{\text{d}h(t)}{\text{d}t}\frac{1}{2{c}_{p}}{\left(\frac{2R}{d}\right)}^{4}{\left(\frac{\text{d}h(t)}{\text{d}t}\right)}^{3}\\ \text{\hspace{1em}\hspace{1em}}+\frac{{\dot{Q}}_{f}(t)}{\delta {c}_{p}{R}^{2}\pi},\end{array}$$(2)
where T(t) is the average temperature of the control volume and T_{B}(t) is the fuel temperature at CSB. If the temperature was spatially homogeneous in the system (i.e., setting T(t) equal to T_{B}(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 0D model, this spatial feature is provided with the difference T(t) − T_{B}(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) − T_{B}(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 inlettocore temperature difference) and 0 (the mean salt temperature coincides with the outlet temperature at the draining end) respectively, it yields: $$\begin{array}{l}T(t){T}_{B}(t)\\ =\text{\Delta}{T}_{0}\left(2\left(1\frac{t}{{t}^{*}}\right){e}^{{\lambda}_{\text{mix}}t}2\frac{t}{{t}^{*}}{e}^{{\lambda}_{\text{mix}}({t}^{*}t)}\right)\end{array}$$(3)
where ΔT_{0} 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.
Fig. 3 Difference between average and outflow molten salt temperatures as a function of time. Red: CFD simulation of the drainingtransient; blue: analytical bestfit. 
3.2 Neutronics submodel
The molten salt system neutron kinetics is approximated by the point kinetics equations [10]. Defining the normalized neutron population η(t) = n(t)∕n_{0} and the normalized precursor concentration ξ(t) = C(t)∕C_{0}, where n_{0} and C_{0} 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: $$\{\begin{array}{c}\frac{\text{d}\eta (t)}{\text{d}t}=\frac{\rho (t)\beta}{\text{\Lambda}}\eta (t)+\frac{\beta}{\text{\Lambda}}\xi (t)\\ [1em]\frac{\text{d}\xi (t)}{\text{d}t}={\lambda}_{p}\left(\eta (t)\xi (t)\right)\end{array}.$$(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 CSB 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 timescale 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 onegroup 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: $$\begin{array}{ll}\delta \rho \hfill & =\frac{\delta k}{{k}^{2}}\simeq \frac{\delta k}{k}\hfill \\ \hfill & =\frac{\delta {k}_{\infty}(T)}{{k}_{\infty}(T)}\frac{{L}^{2}(T){B}^{2}(h)}{1+{L}^{2}(T){B}^{2}(h)}\left(\frac{\delta {B}^{2}(h)}{{B}^{2}(h)}+\frac{\delta {L}^{2}(T)}{{L}^{2}(T)}\right).\hfill \end{array}$$(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 nonleakage 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: $${\text{\Sigma}}_{x}(T)={({\Sigma}_{x})}_{0}+{\alpha}_{{\text{\Sigma}}_{x}}\mathrm{ln}\left(\frac{T}{{T}_{\text{ref}}}\right)$$(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 rewritten as $$\frac{\delta {k}_{\infty}}{{k}_{\infty}}={\alpha}_{T}^{*}(T)\delta T$$(7)
where ${\alpha}_{T}^{*}(T)$ is defined as the temperature reactivity feedback coefficient related to the variation of k_{∞} : $${\alpha}_{T}^{*}(T)=\frac{1}{T}\left(\frac{{\alpha}_{\nu {\text{\Sigma}}_{f}}}{\nu {\text{\Sigma}}_{f}(T)}\frac{{\alpha}_{{\text{\Sigma}}_{a}}}{{\text{\Sigma}}_{a}(T)}\right).$$(8)
Two effects are present in the ${\alpha}_{T}^{*}(T)$ 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 ${\alpha}_{\nu {\text{\Sigma}}_{f}}$, along with an increase of absorptions, represented by the positive value of ${\alpha}_{{\text{\Sigma}}_{a}}$. As general effect, the infinite multiplication factor decreases following a temperature increase, and therefore the temperature coefficient ${\alpha}_{T}^{*}(T)$ 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: $${B}^{2}(h)={\left(\frac{{j}_{0,1}}{R}\right)}^{2}+{\left(\frac{\pi}{Hh}\right)}^{2},$$(9)
where j_{0,1} is the first zero of the zeroth 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: $$\frac{{L}^{2}(T){B}^{2}(h)}{1+{L}^{2}(T){B}^{2}(h)}\frac{\delta {B}^{2}(h)}{{B}^{2}(h)}={\alpha}_{h}(h,T)\delta h$$(10)
where α_{h}(h, T) is defined as the reactivity coefficient due to the geometry modification: $${\alpha}_{h}(h,T)=\frac{{L}^{2}(T)}{1+{L}^{2}(T){B}^{2}(h)}\frac{2{\pi}^{2}}{{(Hh)}^{3}}.$$(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 nonleakage due to a variation of the diffusion area L^{2}. 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 L^{2} (T) = D_{n}(T)∕Σ_{a}(T), introducing the temperature relations (6) and computing the derivative with respect to temperature, it is obtained: $$\frac{{L}^{2}(T){B}^{2}(h)}{1+{L}^{2}(T){B}^{2}(h)}\frac{\delta {L}^{2}(T)}{{L}^{2}(T)}={\alpha}_{T}^{**}(h,T)\delta T,$$(12)
where ${\alpha}_{T}^{**}(h,T)$ is the temperature reactivity feedback coefficient due to the diffusion area variation: $${\alpha}_{T}^{**}(h,T)=\frac{{L}^{2}(T){B}^{2}(h)}{1+{L}^{2}(T){B}^{2}(h)}\frac{1}{T}\left(\frac{{\alpha}_{D}}{{D}_{n}(T)}\frac{{\alpha}_{{\text{\Sigma}}_{a}}}{{\text{\Sigma}}_{a}(T)}\right).$$(13)
Temperature increase implies on one hand the increase of absorptions of neutrons in the medium, represented by the positive value of ${\alpha}_{{\text{\Sigma}}_{a}},$ but on the other the reduction of neutron diffusion coefficient, stated by the slight negative value of ${\alpha}_{{D}_{n}}.$ 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 ${\alpha}_{T}^{**}(h,T)$ 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 ${\alpha}_{T}^{**}(h,T)$.
On the basis of the previous reasoning, a total temperature feedback coefficient can be defined, as the sum of the temperature coefficients aforementioned: $$\begin{array}{ll}{\alpha}_{T}(h,T)\hfill & ={\alpha}_{T}^{*}(T)+{\alpha}_{T}^{**}(h,T)\hfill \\ \hfill & =\frac{1}{T}(\frac{{\alpha}_{\nu {\text{\Sigma}}_{f}}}{\nu {\text{\Sigma}}_{f}(T)}\frac{1}{1+{L}^{2}(T){B}^{2}(h)}\frac{{\alpha}_{{\text{\Sigma}}_{a}}}{{\text{\Sigma}}_{a}(T)}\hfill \\ \hfill & \text{\hspace{1em}}\frac{{L}^{2}(T){B}^{2}(h)}{1+{L}^{2}(T){B}^{2}(h)}\frac{{\alpha}_{D}}{{D}_{n}(T)}).\hfill \end{array}$$(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: $$\delta \rho ={\alpha}_{T}(h,T)\text{\hspace{0.17em}}\delta T+{\alpha}_{h}(h,T)\text{\hspace{0.17em}}\delta h.$$(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: $$\begin{array}{ll}\rho (h,T)\hfill & ={\displaystyle {\int}_{{T}_{\text{ref}}}^{T}{\alpha}_{T}}(h,{T}^{\prime})\text{\hspace{0.17em}}\delta {T}^{\prime}+{\displaystyle {\int}_{{h}_{\text{ref}}}^{h}{\alpha}_{T}}({h}^{\prime},T)\text{\hspace{0.17em}}\delta {h}^{\prime}\hfill \\ \hfill & ={\rho}_{T}(h,T)+{\rho}_{h}(h,T),\hfill \end{array}$$(16)
where the reference conditions for temperature and core cavity height are T_{ref} and h_{ref}, respectively, and it is assumed that the reference reactivity ρ(h_{ref}, T_{ref}) 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 T_{ref} is imposed equal to 900 K [9] and the corresponding reference monitor length h_{ref} is obtained imposing criticality condition to the system.
Fig. 4 Temperature feedback coefficients as function of temperature. The coefficient ${\alpha}_{T}^{**}(h,T)$ depends alsoon the monitor length h. 
3.3 Heat source from fissions and decay heat
The power source term in equation (2), ${\dot{Q}}_{f}(t)$, 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 ${\dot{q}}_{\text{fis}}(t)$ at a certain time instant is expressed by the following formula $$\begin{array}{ll}{\dot{q}}_{\text{fis}}(t)\hfill & ={E}_{f}{\text{\Sigma}}_{f}(T)\text{UpPhi}(t)={E}_{f}{\text{\Sigma}}_{f}(T)vn(t)\hfill \\ \hfill & ={E}_{f}{\text{\Sigma}}_{f}(T)v{n}_{0}\eta (t),\hfill \end{array}$$(17)
where E_{f}Σ_{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 0D model, ${\dot{q}}_{\text{fis}}(t)$ 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: $$\frac{\text{d}{q}_{d}(t)}{\text{d}t}=f{\dot{q}}_{\text{fis}}(t){\lambda}_{d}{q}_{d}(t),$$(18)
where q_{d}(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 timedependent: $${\dot{Q}}_{f}(t)=((1f){\dot{q}}_{\text{fis}}(t)+{\lambda}_{d}{q}_{d}(t))\text{\hspace{0.17em}}{R}^{2}\pi \left(Hh(t)\right).$$(19)
3.4 Molten salt level submodel
The dynamics of the fuel salt draining depends strongly on the free surface velocity v_{A} (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 wellestablished 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 quasisteady 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]): $${v}_{B}(t)=\theta \sqrt{2g\left(Hh(t)+L\right)+2\frac{{p}_{A}{p}_{B}}{\delta}},$$(20)
where θ is a dimensionless constant that describes the friction losses along the core cavity and the shaft. The mathematical expression for θ is $$\theta =\sqrt{\frac{1}{4{C}_{p}\frac{L}{d}+{K}_{c}+1}},$$(21)
where C_{p} is the Fanning friction factor along the shaft (which is a function of the Reynolds number, and indeed, of the salt viscosity) and K_{c} the local friction factor for the cavitytoshaft sudden contraction, which is assumed to be 0.5 [16–18].
We now reconsider equation (1), expressing the unsteady mass balance of the salt draining. Substituting the relation for v_{B}(t) from equation (20), one obtains: $${v}_{A}(t)={\left(\frac{d}{2R}\right)}^{2}\theta \sqrt{2g(Hh(t)+L)+2\frac{{p}_{A}{p}_{B}}{\delta}}.$$(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: $$\begin{array}{ll}h(t)\hfill & =H\text{}+\text{}L\text{}\text{}{\left(\text{}\sqrt{H\text{}+\text{}L+\frac{{p}_{A}{p}_{B}}{\delta g}}\sqrt{\frac{g}{2}}{\left(\frac{d}{2R}\right)}^{2}\theta t\right)}^{2}\hfill \\ \hfill & \text{\hspace{1em}}+\frac{{p}_{A}{p}_{B}}{\delta g}.\hfill \end{array}$$(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 nonisothermal transients [8].
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${\dot{Q}}_{f}(0)$ 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.
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 nonlinearity. 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 semianalytical 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 variableorder 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 T_{ref} = 900 K (Tab. 3). Conversely, the molten salt volume is larger than the critical one (h_{ref} = 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× 10^{13} n/cm^{3} to about 8.3 × 10^{13}, 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 Nickelbased 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 T_{B}(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.
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 geometryrelated ρ_{h}(h, T) and the temperaturerelated ρ_{T}(h, T) system reactivities, respectively. 
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 thermophysical and some neutronics parameters of the molten salt fuel. The parametric analysis is performed with a onefactorattime method, varying a parameter at a time while keeping others fixed [21]. An uncertainty of ±10 % on thermophysical 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 frictionrelated 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.
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% (dashdot lines). 
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% (dashdot lines). 
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 fullypassive safety system concept developed for liquidfueled 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 zerodimensional semianalytical 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 quasisteady approach, to close the problem. The 0D 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 volumechange 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 semianalytical 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.
Nomenclature
${({\Sigma}_{x})}_{0}$: Macroscopic cross section at T_{ref} for Σ_{x} (1/m)
${\alpha}_{T}^{*}$: Temperature reactivity feedback coefficient related to variation of k_{∞} (1/K)
${\alpha}_{T}^{**}$: Temperature reactivity feedback coefficient due to L^{2} (1/K)
α_{h}: Reactivity coefficient due to geometry change (1/m)
${\alpha}_{{\Sigma}_{x}}$: Doppler coefficient for Σ_{x} (1/m)
${\alpha}_{{D}_{n}}$: Doppler coefficient for D_{n} (m)
β: Total effective delayed neutron fraction (1)
δ: Salt mass density (kg/m^{3})
η: 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 (m^{2}/s)
ρ_{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)
$\stackrel{.}{{Q}_{f}}$: Source of thermal power (W)
${\dot{q}}_{\text{fis}}$: Total fission volumetric power (W/m^{3})
B^{2}: Geometrical buckling of the system (m^{−2})
C_{0}: Precursor concentration at equilibrium (1/m^{3})
c_{p}: Specific heat at constant pressure (J/kg/K)
D_{n}: Diffusion coefficient (m)
f: Total fraction of decay heat precursors (1/s)
k: Effective multiplication factor (1)
k_{∞}: Infinite multiplication factor (1)
n_{0}: Neutron concentration at equilibrium (1/m^{3})
q_{d}: Decay energy stored in precursors (J/m^{3})
T_{B}: Outlet salt temperature (K)
v_{A}: Salt free surface velocity (m/s)
v_{B}: Salt outflow velocity (m/s)
References
 M. Allibert, M. Aufiero, M. Brovchenko, S. Delpech, V. Ghetta, D. Heuer, A. Laureau, E. MerleLucotte, 7Molten 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]
 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 ThU support: Fuel cycle flexibility and key material properties, Ann. Nucl. Energy 64, 408 (2014) [CrossRef] [Google Scholar]
 M. Brovchenko, E. MerleLucotte, 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 preconceptual design of the MSFR, WorkPackage WP2, Deliverable D2.2, EVOL, Evaluation and Viability of Liquid fuel fast reactor system) European FP7 project, Contract number: 249696, EVOL (2014) [Google Scholar]
 J. Serp, M. Allibert, O. Benes, S. Delpech, O. Feynberg, V. Ghetta, D. Heuer, D. Holcomb, V. Ignatiev, J.L. Kloostermang, L. Luzzi, E. MerleLucotte, 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]
 S. Wang, M. Massone, A. Rineiski, E. MerleLucotte, 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]
 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]
 Y.A. Çengel, J.M. Cimbala, Fluid Mechanics: Fundamentals and applications (McGrawHill Higher Education, 2004) [Google Scholar]
 F. Di Lecce, Neutronic and thermalhydraulic simulations for Molten Salt Fast Reactor safety assessment (unpublished master thesis), Politecnico di Torino, Torino, Italia, 2018 [Google Scholar]
 M. Aufiero, Development of advanced simulation tools for circulating fuel nuclear reactors, Doctoral dissertation, 2014 [Google Scholar]
 A.E. Waltar, A.B. Reynolds, Fast Breeder Reactors (Pergamon Press, 1981) [Google Scholar]
 D.L. Hetrick, Dynamics of Nuclear Reactors (The University of Chicago Press, 1971) [Google Scholar]
 A.E. Waltar, D.R. Todd, P.V. Tsvetkov, Fast Spectrum Reactors (Springer, 2012) [CrossRef] [Google Scholar]
 E.L. Lewis, Fundamentals of Nuclear Reactor Physics (Elsevier Science Publishing Co Inc., 2008) [Google Scholar]
 M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions (United States Department of Commerce, 1964) [Google Scholar]
 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]
 T.L. Bergman, A.S. Lavine, F.P. Incropera, D.P. DeWitt, Introduction to Heat Transfer, 6th edn. (Wiley, 2011) [Google Scholar]
 F.A. Morrison, An Introduction to Fluid Mechanics (Cambridge University Press, 2013) [CrossRef] [Google Scholar]
 Crane, Flow of Fluids through Valves, Fittings, and Pipe (Crane Co. Engineering Division, 1957) [Google Scholar]
 MATLAB and Statistics Toolbox Release R2017a, The MathWorks, Inc., Natick, Massachusetts, United States [Google Scholar]
 G.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations (PrenticeHall, Englewood Cliffs, N.J., 1971) [Google Scholar]
 V. Czitrom, OneFactorataTime Versus Designed Experiments (American Statistician, 1999) [Google Scholar]
Cite this article as: Francesco Di Lecce, Antonio Cammi, Sandra Dulla, Stefano Lorenzi, and Piero Ravetto. Simplified 0D semianalytical model for fuel draining in molten salt reactors, EPJ Nuclear Sci. Technol. 5, 14 (2019)
All Tables
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
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 timevarying control volume (in red) is adopted to describe the draining phenomenon. 

In the text 
Fig. 2 Graphical scheme depicting the multiphysics interaction of the molten salt fuel draining phenomenon. The dashed arrow identifies the coupling between fluiddynamics and energy balance through the temperaturedependent salt properties, currently not modelled. 

In the text 
Fig. 3 Difference between average and outflow molten salt temperatures as a function of time. Red: CFD simulation of the drainingtransient; blue: analytical bestfit. 

In the text 
Fig. 4 Temperature feedback coefficients as function of temperature. The coefficient ${\alpha}_{T}^{**}(h,T)$ depends alsoon the monitor length h. 

In the text 
Fig. 5 Molten salt level evolution in time. Initial value H = 2.255 m (full cavity). 

In the text 
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 geometryrelated ρ_{h}(h, T) and the temperaturerelated ρ_{T}(h, T) system reactivities, respectively. 

In the text 
Fig. 7 Simulation results of the 0D semianalytical model with input parameters as in Table 4. 

In the text 
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% (dashdot lines). 

In the text 
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% (dashdot lines). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.