Simplified 0-D semi-analytical model for fuel draining in molten salt reactors

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 zero-dimensional semi-analytical 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 (temperature-related internal surface damages) and neutronic (sub-critical states dominate the transient) view points and show which are the main characteristics of the fuel salt draining transient.


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 * e-mail: antonio.cammi@polimi.it 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.

Molten salt 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.

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 temperature and 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.

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 which is a balance on the volumetric flow rate at the CS-A and CS-B. 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 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 sourceQ f (t), linked to neutronics (prompt fission power and decay heat). The energy equation can be written as follows (for more detail, see [8]): where T (t) is the average temperature of the control volume and T B (t) is the fuel temperature at CS-B. 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 0-D 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 inlet-to-core temperature difference) and 0 (the mean salt temperature coincides with the outlet temperature at the draining end) respectively, it yields: 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.

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)/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: .
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 modification of 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 diffusion area is equal to the diffusion coefficient divided by the absorption cross section, it is possible to express the variation of ρ as: 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.

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: where Σ x is the macroscopic cross section for 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 where α * T (T ) is defined as the temperature reactivity feedback coefficient related to the variation of k ∞ : Two effects are present in the α * 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 α νΣ f , along with an increase of absorptions, represented by the positive value of α Σa . As general effect, the infinite multiplication factor decreases following a temperature increase, and therefore the temperature coefficient α * T (T ) is negative. Figure 4a depicts the dependence of this coefficient on temperature, using the cross section values listed in Table 3.

Variation of buckling
The geometrical buckling of the cylindrical multiplying domain is defined as: where j 0,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 respect to h, the following result is obtained: where α h (h, T ) is defined as the reactivity coefficient due to the geometry modification: 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.

Variation of diffusion area
The last term in equation (5)  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: where α * * T (h, T ) is the temperature reactivity feedback coefficient due to the diffusion area variation: Temperature increase implies on one hand the increase of absorptions of neutrons in the medium, represented by the positive value of α Σa , but on the other the reduction of neutron diffusion coefficient, stated by the slight negative value of α Dn . 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 α * * 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 α * * 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: 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: 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:

Heat source from fissions and decay heat
The power source term in equation (2),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 powerq fis (t) at a certain time instant is expressed by the following formulȧ 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 0-D model,q 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: 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 time-dependent:

Molten salt level sub-model
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 dh/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 temperature draining 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] where θ is a dimensionless constant that describes the friction losses along the core cavity and the shaft. The mathematical expression for θ is 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 cavity-to-shaft sudden contraction, which is assumed to be 0.5 [16][17][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: 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: 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].

Initial conditions for the coupled simulations
Once the function h(t) is known, the coupled equations describing the draining can be solved, provided initial Table 4. Initial conditions for the coupled calculations.
Average fuel temperature T0 700 • C Normalized neutron density η(0) 1 Normalized precursor concentration ξ(0) 1 Total system thermal powerQ f (0) 3.5 × 10 9 W Neutron population n0 1.33 × 10 13 neutron/m 3 Neutron precursors concentration conditions are defined. As for the reference MSFR operating parameters [1], the average salt temperature is set to 700 • C and the initial powerQ 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.

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 temperature difference (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, as the 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 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 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 temperature ramp 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 temperature starts decreasing to reach 1025 K at the transient end. It should be noticed that the time derivative of the temperature depends 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.

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 temperature variations 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  Table 4.  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.

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 quasisteady 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 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.