Development of a multiphysics model for the study of fuel compressibility effects in the Molten Salt Fast Reactor

Compressible fluid dynamics is of great practical interest in many industrial applications, ranging from chemistry to aeronautical industry, and to nuclear field as well. At the same time, modelling and simulation of compressible flows is a very complex task, requiring the development of specific approaches, in order to describe the effect of pressure on the fluid velocity field. Compressibility effects become even more important in the study of two-phase flows, due to the presence of a gaseous phase. In addition, compressibility is also expected to have a significant impact on other physics, such as chemical or nuclear reactions occurring in the mixture. In this perspective, multiphysics represents a useful approach to address this complex problem, providing a way to catch all the different physics that come into play as well as the coupling between them. In this work, a multiphysics model is developed for the analysis of the generation IV Molten Salt Fast Reactor (MSFR), with a specific focus on the compressibility effects of the fluid that acts as fuel in the reactor. The fuel mixture compressibility is expected to have an important effect on the system dynamics, especially in very rapid super-prompt-critical transients. In addition, the presence of a helium bubbling system used for online fission product removal could modify the fuel mixture compressibility, further affecting the system transient behaviour. Therefore, the MSFR represents an application of concrete interest, inherent to the analysis of compressibility effects and to the development of suitable modelling approaches. An OpenFOAM solver is developed to handle the fuel compressibility, the presence of gas bubbles in the reactor as well as the coupling between the system neutronics and fluid dynamics. The outcomes of this analysis point out that the fuel compressibility plays a crucial role in the evolution of fast transients, introducing delays in the expansion feedbacks that strongly affect the system dynamics. Moreover, it is found that the gas bubbles significantly alter the fuel compressibility, yielding even larger differences compared to the incompressible approximation usually adopted in the current MSFR solvers. 2018 The Authors. Published by Elsevier Ltd. This is an openaccess article under the CCBY license (http:// creativecommons.org/licenses/by/4.0/).


Introduction
Compressibility plays a crucial role in the propagation of waves within a fluid, since pressure and density perturbations travel at a finite velocity through the medium. The relative propagation velocity of the waves with respect to the fluid represents the local speed of sound. Although all the real fluids are compressible, this property can be often neglected, introducing the incompressibility assumption and assuming an infinite speed of sound (Thompson, 1972). However, there are many industrial applications, ranging from high-pressure chemistry to supersonic aerodynamics, in which the incompressible approximation is not suitable, since the fluid density is strongly affected by pressure. Compared to incompressible fluid dynamics, the study of compressible flows requires a specific treatment both from a theoretical and numerical point of view, calling for the coupled solution of the continuity, momentum and energy equations (Moukalled et al., 2016). In addition, the density variation may have important feedbacks on other physics involved in the problem, further complicating the analysis of many complex systems of industrial interest (e.g., chemical and nuclear reactors).
Compressibility effects become of particular interest in the analysis of two-phase liquid-gas flows. In this situation, the presence of gas bubbles may have a relevant impact on the average mixture compressibility. In addition, local effects may also arise in case of strongly heterogeneous spatial distributions of the gaseous phase, leading to phenomena that cannot be caught with single-phase or homogeneous-mixture approaches.
In this sense, the multiphysics approach constitutes a valuable tool to address the problem, providing an efficient way to describe all the different physical phenomena occurring in an industrial process (Cammi et al., 2011(Cammi et al., , 2012Aufiero et al., 2014aAufiero et al., , 2014bFiorina et al., 2014). In the present work, a multiphysics modelling approach is presented for the analysis of the impact of fuel compressibility during super-prompt-critical transients in the generation IV Molten Salt Fast Reactor (GIF, 2016).
The Molten Salt Fast Reactor (MSFR), developed in the framework of the H2020 SAMOFAR Project (http://samofar.eu/), is a circulating fuel nuclear reactor in which a mixture of molten thorium and uranium fluorides acts as fuel and coolant simultaneously (Serp et al., 2014;Dolan, 2017;Gerardin et al., 2017). From a computational point of view, the simulation of nuclear reactor dynamics is a complex multiphysics task, needing accurate solution for both neutronics and thermal hydraulics and considering the coupling between them. This is even more important in circulating fuel nuclear reactors, in which the velocity field of the fuel salt mixture has a significant impact on the distribution of the delayed neutron precursors, affecting the reactor kinetics (a brief descrip-tion of the MSFR is provided in Section 2.1). Given this tight coupling between neutronics and fluid dynamics, the effect of the fuel mixture compressibility may play a relevant role on the dynamics of the system and on the transient behaviour of the reactor.
Due to the negative temperature feedback coefficient of the MSFR (Gerardin et al., 2017), temperature increases during power excursions lead to a reduction of the system reactivity. This feedback is partly due to the Doppler effect, related to the neutron captures by the fertile nuclides, and partly to the thermal expansion of the fuel, which increases neutron leakages. While the Doppler effect acts promptly to reduce the system reactivity, the expansion feedback is delayed, since a density perturbation takes a finite time to propagate through the reactor.
In the MSFR, the speed of sound in the fuel mixture is about 1200 m/s, hence sufficiently high to consider the pressure wave propagation connected to the fluid compressibility as ''instantaneous" in most transient scenarios (Aufiero et al., 2017). On the other hand, this may not be the case for very rapid superprompt-critical transients, which could be a reason of concern during the reactor start-up, due to unwanted fuel injections. In fact, the characteristic times of these transients are in the order of a few milliseconds, comparable to the propagation time of pressure waves in the reactor. This could lead to a delay of the expansion mechanism, resulting in an overall weaker feedback. For this reason, the adoption of incompressible approximation (Aufiero et al., 2014a(Aufiero et al., , 2014bFiorina et al., 2014) may result in significant underestimations of the energy released in super-prompt-critical transients. The analysis of these strongly coupled transients provides meaningful information also from a safety point of view, as highlighted by Qiu et al. (2016) and Zhang et al. (2018).
In addition, a helium bubbling system is envisaged in the MSFR to efficiently remove the gaseous fission products in the salt, but also as a possible option for the reactivity control of the reactor, exploiting the void reactivity feedback of the air bubbles in the fuel  (Gerardin et al., 2017). The helium bubbles are injected at the core inlet and are separated from the fuel mixture at the core outlet. The presence of bubbles in the fuel mixture, as well as their spatial distribution, may have a significant impact on the fuel compressibility and on the pressure wave propagation, possibly leading to additional effects that cannot be described with the incompressible approximation. In particular, the bubbles may increase (both locally and globally) the fuel compressibility, further delaying the thermal expansion feedback. In the light of this, the MSFR constitutes an application of industrial interest, inherent at the same time to the modelling of two-phase, compressible fluid dynamics and its coupling with other physics. The purpose of the present work is to model the compressibility phenomena in the pure liquid salt, and in the bubbly mixture as well, highlighting the coupling between these phenomena and the system neutronics. This problem is addressed by means of a multiphysics OpenFOAM (http://www.openfoam.org/docs/) solver for the MSFR developed in a previous work (Cervi et al., 2017), implementing a two-phase Euler-Euler model for the system thermal hydraulics and the multi-group formulation of the diffusion equation for neutronics. Thanks to the Euler-Euler approach, the motion of the helium bubbles through the fuel mixture can be modelled. In addition, both the liquid and the gaseous phases are treated as compressible fluids. Coupling strategies between the different physics are also developed, in order to account for the presence of the bubbly flow in the reactor from both a neutronics and a thermal hydraulics point of view. A previous analysis carried out with this solver (Cervi et al., 2017) pointed out that an accurate description of the bubble distribution in the reactor is fundamental for the correct evaluation of the reactor multiplication factor and of the void reactivity coefficient. However, the impact of the helium bubbling system and of compressibility effects on the MSFR transient behaviour is still to be investigated. This work aims at filling this gap, highlighting the effect of the fuel compressibility on super-prompt-critical transients, in the pure liquid salt as well as in presence of the helium bubbles. Compressibility effects in super-prompt-critical transients are pointed out by Aufiero et al. (2017), assuming a uniformly distributed gas phase. Compared to the available literature, the novelty of the present work is that the local effect of the bubble spatial distribution on the fuel mixture compressibility is highlighted. To this aim, pressure wave propagation through the reactor is studied (i) considering the bubble distribution calculated by the multiphysics solver, and (ii) assuming a uniform distribution with the same average void fraction. In addition, the impact of the pure liquid salt compressibility, without helium bubbles, is also investigated.
The paper is organised as follows. Section 2 provides a brief description of the Molten Salt Fast Reactor. In Section 3, the main features of the OpenFOAM solver are presented. In Section 4.1, the solver is tested by simulating a super-prompt-critical reactivity insertion in a simplified model of the MSFR, considering only the liquid fuel without bubbles in it. In Section 4.2, the presence of the helium bubbling system is also considered, in order to show the effect of the bubbles on the mixture compressibility, and on the system dynamics as well. The comparison and discussion of results is given in Section 4.3, while the main conclusions of the present work and its future developments are drawn in Section 5.

Description of the Molten Salt Fast Reactor
In this Section, the MSFR is briefly described in order to highlight some peculiarities of interest for the present work. A schematic representation of the MSFR layout is shown in Fig. 1, while its main design parameters are listed in Table 1 (Gerardin et al., 2017). The salt mixture, which serves simultaneously as fuel and coolant, circulates through sixteen external circuits, where the power produced by fissions is removed by heat exchangers. Compared to traditional solid-fuelled reactors, the circulating fuel is a distinguishing feature of the MSFR (and other Molten Salt Reactors as well -e.g., see Serp et al. (2014), Dolan (2017) -which arises completely new design and technological challenges. Notably, the delayed neutron precursors do not decay in the place they are produced, as in conventional nuclear systems, but they are dragged by the fuel mixture through the reactor and the external circuits. For this reason, delayed neutrons can be emitted in peripheral regions of the reactor, where the neutron importance is lower, or even in the external circuit, where they do not contribute to fissions. As a consequence, the coupling between the system neutronics and thermal fluid dynamics is even stronger than in traditional reactors, since the fuel mixture velocity field directly influences the precursor distribution.

Modelling approach
In this Section, a multiphysics solver for the analysis of the MSFR is presented, aimed at addressing the issue of the fuel compressibility and its dependence on the bubble spatial distribution. More details on the solver can be found in (Cervi et al., 2017).
At each time step, the system thermal hydraulics and neutronics are solved in two different cycles, as indicated in Fig. 2. The thermal hydraulics sub-solver is based on the standard OpenFOAM solver ''twoPhaseEulerFoam" for the compressible fluid and the bubble modelling. The neutronics sub-solver implements the multi-group diffusion equation for the neutron flux and the transport equations of the delayed neutron precursors.
At the beginning of the time step, the thermal hydraulics cycle solves for the phase fractions, for the velocity of both phases, for the pressure and for the temperature. Picard iterations are performed until convergence is reached for the solution of the thermal hydraulic part of the problem. Then, the neutronics cycle begins, solving for the flux, for the delayed neutron precursors and for the decay heat. Once the flux (the fission power in turn) and the decay heat are known, the volumetric power source field is updated and the energy equation is solved again. Once the new temperature and density fields of the fuel are calculated, the cross sections are updated, also correcting for the void fraction (see Section 3.3, Eq. (18), and the cycle is repeated with Picard iterations until convergence is reached. In addition, a certain number of external iterations between the thermal hydraulics and the neutronics sub-solvers is performed. The external iterations are particularly important in fast transients, in which the large thermal expansions due to steep power excursions have a strong impact on fuel velocity field.
Summarizing, the following information is exchanged between the two sub-solvers. The void fraction, fuel temperature and density are passed from thermal-hydraulics to neutronics, in order to evaluate the cross sections. On the other hand, the power density distribution is passed from the neutronics to the thermalhydraulics solver in order to update temperature. Then, once convergence is reached in the external iterations between neutronics and thermal-hydraulics, the governing equations are integrated in time, moving to the following step.

Neutronics model
The multi-group diffusion equation is adopted for the estimation of the neutron flux. Despite some limitations of this modelling choice, it is well employed in the transient analysis of nuclear reactors (Fiorina et al., 2015) due to the easiness in the implementation and the relatively low computational time. For each energy group, the neutron diffusion equation reads: In the previous equation, S n;i , S d , S s;i are explicit terms which represent the neutrons released in fissions induced by neutrons from other energy groups, the delayed neutrons and the neutrons scattered from the other groups, respectively: Due to these explicit terms, an iterative procedure among the several groups is required to achieve convergence for the neutronics description. The value of k eff can be suitably changed to simulate a given reactivity insertion. Albedo boundary conditions are adopted at the top and bottom walls of the reactor (axial reflectors) and at the radial wall (blanket salt), in order to limit the solution of the neutron diffusion equation to the fuel salt circuit only (Aufiero et al., 2014b).
The precursor balance equations include the diffusion and the transport term to allow for the fuel motion (neglecting the precursor mass transfer from the liquid to the gas phase): The turbulent Schmidt number Sc t is set to 0.85, even if no data are specifically available for the diffusion of species in the MSFR salt (Aufiero et al., 2014b).
In order to properly consider the decay heat during accidental transients, the solver is provided with equations that consider the behaviour of the isotopes responsible for the decay heat, subdivided in ''decay heat groups" in a manner similar to the precursor  groups. Actually, the equations implement the balance for the precursor concentration multiplied by the average energy released by that decay group: In addition, a power iteration routine, based on the k-eigenvalue method, is implemented in the neutronics module of the solver for the calculation of the multiplication factor. For a more detailed description, the reader is referred to (Cervi et al., 2017). At the beginning of the simulation, the user can choose between the time-dependent solver or the steady-state eigenvalue solver.

Thermal hydraulics model
The need of a solver for a two-phase compressible fluid is due to the role that the salt compressibility plays in fast transients, and to the presence of the online bubbling system foreseen for fission product removal and reactivity control. To this aim, the ''twoPhaseEulerFoam" solver available in the OpenFOAM library is used, which implements a two-fluid (or Euler-Euler) approach (Rusche, 2002). Each phase is treated as a continuum interpenetrating each other, and is described by means of averaged conservation equations. Due to the averaging process, phase fractions are introduced into the governing equations.
Compared to other approaches for bubbly flows (e.g., Lagrangian-Lagrangian and Eulerian-Lagrangian), the Euler-Euler approach is characterized by lower computational requirements, and therefore is suitable to the simulation of high-Reynolds and large-scale systems, which is the case for the MSFR. In fact, considering an average fuel density of 4125 kg/m 3 , an average fuel velocity of about 1.2 m/s, a diameter of 2.26 m and a dynamic viscosity of 10 À2 Pa s (Gerardin et al., 2017), the Reynolds number is in the order of 10 6 , implying a fully turbulent flow regime. For these reasons, the Euler-Euler approach is the preferred method for many practical applications, and is adopted also in this work. On the other hand, a VOF (volume of fluid) approach would not be suitable for this analysis. In fact, in dispersed bubbly flows the gas-liquid interfaces are characterized by small scales compared to the system dimensions, and surface tracking at such high level of detail would result in an excessive computational burden.
Due to the loss of information related to the averaging process, several closure relations appear in the macroscopic balance equations. Therefore, the predictive capabilities of the adopted model crucially rely on the accuracy of these constitutive equations, representing a possible modelling limit of the Euler-Euler approach.
The mass and momentum conservation for the two phases reads: A mass source term S j is considered in the continuity equation of the gas phase to model the bubble injection/extraction. The term M j appears in the averaged momentum equations of each phase due to non-linearity, which requires closure equations. This term takes into account the momentum transfer between the two phases, due to the forces acting at the liquid-gas interface, namely the drag, virtual mass forces, the lift and turbulent dispersion.
Several models are available in literature to describe the interphase terms and to close the momentum equation. For details, the reader is referred to Enwald et al. (1996) and Gidaspow (1994).
In this work, the Schiller-Naumann correlation (Schiller and Naumann, 1933;Rusche, 2002;Xiao et al., 2017) is selected for the evaluation of the drag coefficient 1 : The bubble Reynolds number Re b is defined as: where d b is the average bubble diameter, estimated by means of the following isothermal power law: in which the reference diameter and pressure are required as input parameters. Typically, the bubble diameter ranges from 1 to 5 mm (Lance et al., 1996) In this work, it is assumed that d 0 ¼ 3 mm at p 0 ¼ 1 atm. An isothermal approximation is acceptable for the purpose of this analysis, since the pressure variations encountered in this work are much higher than temperature variations (see Sections 4.1, 4.2 and Appendix B).
A constant coefficient correlation is adopted for virtual mass forces, with C VM ¼ 0:5 (Paladino and Maliska, 2002;Rusche, 2002). 2 Lift is not considered in the present work, assuming that the bubble size is sufficiently small to neglect the effect of vorticity on the momentum exchange between phases. Moreover, turbulent dispersion is also neglected. The energy equations for ''twoPhaseEulerFoam" reads: @q j a j h j @t þ r Á q j a j u j h j þ @q j a j k j @t þ r Á q j a j u j k j where L is an inter-phase heat transfer coefficient resulting from the averaging process, and DT is the temperature difference between the two phases. In the present work, the Ranz-Marshall correlation (Ranz and Marshall, 1952;Mimouni et al., 2010). is adopted to evaluate L: in which Pr is the Prandtl number for the liquid phase. The provided references show good agreement between the selected correlations and experimental data. In addition, the authors carried out a preliminary sensitivity analysis to other correlations implemented in the ''twoPhaseEulerFoam" solver (also including lift and turbulent dispersion), pointing out a negligible dependence for the present case studies.
The eddy viscosity models are also modified with the addition of new terms to incorporate the effect of the dispersed phase on turbulence. In particular, the Lahey k-e model (Lahey, 2005) is selected, in which the turbulent transport equations are solved for the continuous phase while the turbulent viscosity is corrected by the gas fraction to account for the bubble effect. 1 The force exerted by the liquid on a gas bubble is given by The virtual mass force is given by F VM ¼ agq l Dugas Dt À Dul Dt .

Pressure correction equation and treatment of the compressibility
The continuity and the momentum equations are solved by means of a pressure correction algorithm -for details, see (Patankar, 1980). This method involves the solution of the momentum equation using an initial guessed pressure field. The obtained velocity field is a solution of the momentum equation but it is not granted that it satisfies the continuity equation. Thus, mass conservation is enforced by solving a correction equation for pressure, obtained by a combination of continuity and momentum equations (Moukalled et al., 2016). Then, using the corrected pressure, an updated velocity field is calculated from the momentum equation, and the procedure is iterated until the solution converges.
The continuity equation, Eq. (5), can be rewritten as follows: Combining together Eq. (12) and the momentum equation, Eq. (6), the pressure equation can be written in the following compact form: Expressing the density q j as a function of T and p, its differential is given by: Replacing Eq. (14) in Eq. (13): The first summation in Eq. (15) is simply a source term and can be obtained by solving the energy equation, given the equation of state. However, the second summation depends on pressure, introducing an implicit term into the equation. Expanding the material derivative of pressure, Eq. (15) can be rewritten as: where u j ¼ @q j =@p is the isothermal compressibility of the i-th phase. For incompressible fluids, u j ¼ 0 and the pressure term disappears from Eq. (16). On the other hand, for compressible fluids, u j -0 and the pressure term must be taken into account.
In ''twoPhaseEulerFoam", a distinction is made upon the flow regime of each phase. In particular, in subsonic motion (0:3 < Ma < 0:75), the convective term u j Á rp is small (Thompson, 1972) and therefore is neglected in Eq. (16), only considering the partial time derivative of pressure. On the other side, in transonic (Ma % 1) or supersonic (Ma > 1) motion, the velocity magnitude is comparable to the speed of sound (or even larger, for supersonic flows) and the convective term is included into the pressure equation.
In the case study considered in the present work, the velocity of molten salt and of helium bubbles is far below the sound barrier. However, due to strong density gradients at liquid-gas interfaces, the product u j Á rp is large, and the convective term becomes important even if the velocity is lower than the speed of sound. For this reason, the complete form of Eq. (16) is employed in this work, in order to properly describe the pressure field in case of strongly heterogeneous bubble distributions.
Finally, concerning the boundary conditions, the structure around the liquid domain is considered as ideally rigid, since a detailed thermo-mechanical analysis of the vessel and of the external circuit is out of the scope of the present work. Following this assumption, pressure waves reaching the boundary are totally reflected by the rigid walls (Thompson, 1972).

Coupling between neutronics and thermal hydraulics
Neutronics affects thermal hydraulics through the power source term in the energy equation. The source term is formed by the contributions of fission power, Q f and, decay heat, Q d : On the other hand, thermal hydraulics affects neutronics through temperature and density effects, and the void feedback as well. The expression of the macroscopic cross section for a generic reaction j and for the energy group i can be expressed as follows: The expression is constituted by a constant term, a logarithmic term (accounting for the temperature effect), the ratio between the actual density and a reference density (accounting for the density effect), and a correction term for the gas fraction (accounting for the presence of bubbles). The constant term and the coefficient of the logarithmic term are evaluated by means of the Monte Carlo reactor physics and burnup code SERPENT 2 (Leppänen et al., 2015), using the JEFF-3.1 cross section library (Koning et al., 2006).

Results
The multiphysics solver described in Section 3 is used to simulate a transient in the MSFR, highlighting the impact of the fuel compressibility on the system dynamics. In particular, a superprompt-critical transient at the reactor start-up is chosen as case study in the present work. The power ramp procedure from low initial power is a possible initiator of a reactivity accident, due to an unwanted fuel injection.
A simplified axial-symmetric 2D geometry (Fiorina et al., 2014) is adopted for the present study (Fig. 3a). At t ¼ 0, the system is in zero-power condition, the fuel temperature is 923 K (equal to the design inlet temperature) and the initial reactivity, evaluated with the power iteration routine, is equal to 500 pcm. The power iteration routine is employed on the initial configuration, i.e., at t ¼ 0, before starting the transient simulation. Therefore, the criticality and the time-dependent calculations are carried out independently each from the other. Also note that a 500 pcm reactivity is superprompt-critical in the MSFR, whose delayed neutron fraction is around 310 pcm (Aufiero et al., 2014a).
In all the following cases, the neutron flux is estimated using six energy groups. Six-group cross sections and albedo coefficients are generated using a Serpent model of the same simplified geometry, considering also the reflectors and the blanket (see Fig. 3b).

Effect of the pure salt compressibility on the MSFR dynamics
In this Section, the presence of bubbles in the fuel mixture is not taken into account, in order to highlight only the effects due to the pure salt compressibility. The transient is simulated in two different cases: Compressible fuel (Case I-a): the fuel density depends on temperature and pressure according to the following linear law: where q 0 = 4125 kg/m 3 , b th = 0.882 kg/m 3 K and T ref = 973 K (Ignatiev et al., 2012). The isothermal compressibility w is evaluated as: where K fuel is the bulk modulus of the pure salt, for which a value of 6.3 GPa can be adopted (Aufiero et al., 2017), while at 923 K q fuel = 4169 kg/m 3 . In addition, it is assumed that p ref ¼ 2 bar in Eq. (19).
Incompressible fuel (Case I-b): the fuel density is evaluated by means of Eq. (19), assuming w ¼ 0. In this way, thermal expansion, as well as its effects on neutronics and thermal hydraulics, is taken into account. On the other hand, the dependence of density on pressure is neglected, meaning that pressure waves cannot propagate through the fuel.
A discussion on the result sensitivity to the values of b th and w is presented in Appendix A.
Time evolution of core power and pressure field in Cases I-a and I-b is presented in Figs. 4 and 5, respectively. For the sake of completeness, the temperature fields of the different cases are presented in Appendix B. The power peak value predicted in the compressible case (Case I-a) is 28% higher, compared to the incompressible case (Case I-b). The time evolution of the pressure field is also completely different in the two cases. While in the compressible case (Case I-a) a pressure wave propagates from the centre to the reactor walls, in the incompressible case (Case I-b) no wave is observed. Since in Case I-a the pressure wave takes a finite time to propagate through the reactor, the fuel expansion feedback is delayed, and only the Doppler effect acts promptly to counterbalance the reactivity insertion. For this reason, the power release is higher in Case I-a, compared to Case I-b.

Effect of the fuel mixture compressibility on the MSFR dynamics with bubbles
The same transient of Section 4.1 is now simulated also considering the presence of the bubbles and the helium bubbling system. The injection point as well as the bubble spatial distribution (evaluated with the multiphysics solver) are represented in Fig. 6. The bubbles are injected at the bottom of the reactor and are removed in the heat exchanger region, by considering suitable source terms in Eq. (5). The core-averaged value of the gas fraction is 0.67%. The following three cases are considered: Compressible fuel (Case II-a): the fuel density depends on temperature and pressure according to Eq. (19), while the helium density is evaluated using the equation of state for ideal gases, depending in turn on pressure and temperature. The spatial bubble distribution is calculated by the multiphysics solver and is shown in Fig. 6. Compressible fuel with uniform bubble distribution (Case II-b): as in Case II-a, both the fuel and the helium are treated as compressible fluids. However, a uniform bubble distribution is assumed instead of the one calculated with the solver and presented in Fig. 6. The uniform gas fraction is 0.67%, equal to the core-averaged value of Case II-a. Incompressible fuel (Case II-c): again, the fuel density is evaluated by means of Eq. (19), assuming w ¼ 0, while a density of 0.1 kg/m 3 is assumed for the helium, corresponding to a temperature of 923 K and a pressure of 2 bar. The bubble distribution calculated by the solver is assumed in this case (see Fig. 6).
Due to the void reactivity feedback introduced by the bubbles, the system multiplication factor (evaluated with the power iteration routine implemented into the solver) decreases, compared to Section 4.1, in which the bubbles are not considered. Moreover, the reactivity of Case II-b is different from Cases II-a and II-c, due to the different distribution of the void fraction -for more details on the effect of the bubble spatial distribution on the system reactivity, the reader is referred to (Cervi et al., 2017). For this reason, the value of k eff in Eq. (1) is adjusted so that the initial reactivity is equal to 500 pcm in all the cases. This is to ensure that the differences observed in the various cases are only due to the fuel compressibility, and not to a different initial reactivity insertion. Time evolution of core power and pressure field in Cases II-a, II-b and II-c is presented in Figs. 7 and 8, respectively.
The power peak value predicted in the compressible cases (Cases II-a and II-b) is 64% and 33% higher, respectively, compared to the incompressible one (Case II-c). This indicates that, compared to the pure salt, the helium bubbles increase (both locally and globally) the fuel mixture compressibility, further delaying the  expansion feedback and leading to a larger energy release. This effect is stronger in Case II-a, compared to Case II-b, suggesting that the bubble spatial distribution also plays a fundamental role in determining the system dynamics, locally affecting the fuel mixture compressibility.
The effect of the bubble distribution is also evident from Fig. 8. When a uniform bubble distribution is assumed (Case II-b) the pressure wave has a spherical shape, similarly to the compressible case with pure salt only (Case I-a). However, when the spatial bubble distribution calculated by the solver is considered (Case II-a), the wave-front shape is completely modified by the presence of the bubbly flow. It is also interesting to note that, in Case II-a, wave reflections occur at the reactor wall. This suggests that the wavefront path is also dependent on the reactor geometry, which may have an important role on the system behaviour in superprompt-critical scenarios. As a further note, the elasticity of the reactor walls (here considered as completely rigid) could also have an important effect on wave reflections (this issue is out of the scope of the present work and will be investigated in the future).

Comparison and discussion of results
In all the cases, the power increases of 10 decades, going from zero to about 40-70 GW. Then, the reactor is shut down by the reactivity feedback related to the Doppler effect and to the salt expansion introduced during the power excursion. However, the peak value of the power evaluated in the compressible cases strongly differs, compared to the incompressible ones, both in Sections 4.1 and 4.2. In fact, in the incompressible cases, the thermally-expanded fuel rigidly translates to the expansion tank. On the other hand, in the compressible cases, the pressure wave  due to the thermal expansion travels at the speed of sound, requiring a finite time to propagate through the reactor. Thus, while in the incompressible cases both the Doppler and the expansion feedbacks act promptly to counterbalance the reactivity insertion, in the compressible cases the expansion feedback is delayed, leading to a higher reactivity and, as a consequence, to a larger power excursion. The results of Sections 4.1 and 4.2 are summarized in Table 2.  II-c). Again, pressure waves are eliminated by assumption in this last case. Note that Case II-b is very similar to Case I-a, since in both cases the speed of sound is uniform over the whole medium. However, when the calculated bubble distribution is considered (Case II-a), the wave propagation changes completely.
Even in the pure salt, without considering the helium bubbling system (Case I-a vs. Case I-b), compressibility has a sensible effect on the energy release. Therefore, in the case of super-promptcritical reactivity insertion, the incompressible approximation is not acceptable not only for the gas phase, but also for the liquid salt. Nevertheless, the bubbles play an important role, since they increase the fuel mixture compressibility both locally and globally, further delaying the expansion feedback and leading to a stronger power excursion.
In addition, results are strongly dependent on the bubble spatial distribution. When a uniform bubble distribution is considered, the power increase in the compressible case (Case II-b) with respect to the incompressible case (Case II-c) is 33%. This increase is only slightly larger compared to the pure liquid salt, in which the power peak in the compressible case (Case I-a) is 28% higher with respect to the incompressible case (Case I-b). In fact, both in Case I-a and IIb, the speed of sound through the medium is uniform and, as a consequence, the shape of the wave-front is similar. In addition, the speed of sound in helium is about 1000 m/s, not much smaller than speed of sound in the salt (about 1200 m/s). For this reason, the presence of gas bubbles only yields an additional power increase of 5%, compared to the pure liquid salt. Therefore, even if the bubbles have a visible effect, the most important contribution is due to the liquid compressibility. Again, this confirms that the liquid compressibility should be properly modelled, for an accurate simulation of fast transients.
On the other hand, in Case II-a (compressible fuel mixture with calculated bubble distribution), most of the bubbles are concentrated in the centre of the reactor, as calculated by the multiphysics solver. As a consequence, the mixture compressibility becomes non-uniform, and increases in the central region of the core, where the neutron importance is higher. For this reason, the effect of compressibility on the system reactivity is stronger, leading to a significantly higher power increase (64%). In this case, the presence of bubbles inside the reactor yields an additional power increase of 36%, compared to the pure liquid. Therefore, accurate models are required to predict the bubble spatial distribution. In fact, simply assuming a uniform distribution would result in a strong underestimation of the energy release.

Conclusions and future developments
The Molten Salt Fast Reactor represents an inherent case study for the analysis of two-phase, compressible fluid dynamics and its coupling with other physics, being at the same time an application of scientific and industrial interest. As pointed out in the present work, the helium bubble motion, the propagation of pressure waves in a compressible medium, such as the fuel mixture of the MSFR, and neutronics are strongly coupled each with the other. This is particularly true for fast transients, whose characteristic times are similar to the pressure wave propagation times, at the speed of sound.
In this paper, particular focus is dedicated to the investigation of the impact of the fuel mixture compressibility during accidental super-prompt-critical scenarios, adopting a multiphysics model developed on purpose and implemented in OpenFOAM. A power excursion resulting from a 500 pcm reactivity insertion is evaluated (i) considering the fuel mixture compressibility, and (ii) approximating the fuel mixture as incompressible. The analysis is carried out for the pure salt, and considering the presence of helium bubbles in the reactor as well. Results of this case study point out that approaches neglecting the fuel compressibility may significantly underestimate the power excursions in superprompt-critical reactivity insertions. Therefore, the incompressible approximation turns out to be unsuitable to simulate fast transients in the MSFR (and other liquid-fuelled Molten Salt Reactors as well), not only for the gas phase, but also for the liquid one. In particular, even when the helium bubbling system is not taken into account, the pure salt compressibility has a significant effect on the transient evolution. The system dynamics is also influenced by the bubble spatial distribution, which affects the fuel mixture compressibility as well as the propagation of pressure waves in the medium.
Future efforts will be devoted to the extension of the current analysis to the real 3D geometry of the MSFR, focusing in particular on the geometric effect on the pressure wave reflection. As a further development, the influence of the wall elasticity on the wave reflection could be considered, as well as the wave propagation through the solid reflector, which may have an effect on neutron leakages. As a further note, the fission gas production during the fuel irradiation is also expected to increase the compressibility effects pointed out in this work. In fact, compared to helium, these gases typically have a larger molecular mass and, as a consequence, a lower speed of sound and a higher compressibility (Thompson, 1972). Therefore, considering the fission gas production could be an interesting extension of the present analysis. Moreover, a detailed thermo-mechanical analysis is required to assess whether the pressure waves could challenge the structural integrity of the walls.
As far as neutronics is concerned, the thermal expansion feedback is related to neutron leakages, which in turn depend on the flux behaviour at the reactor walls. In this sense, the neutronics description could benefit from the adoption of a transport model, in order to overcome the limitations of diffusion theory in the flux estimation at material interfaces. of pressure waves, as temperature increases, while the compressibility coefficient makes the waves propagate with a finite velocity through the reactor. Therefore, both these properties must be considered to study the compressibility effects presented in the manuscript. Ignatiev et al. (2012) found the following experimental correlation between fuel density and temperature, without considering the effect of pressure: where q 0 = 4125 kg/m 3 , b th = 0.882 kg/m 3 K and T ref = 973 K. The measurement error is estimated as 0.9% (Ignatiev et al., 2012), even if no specific information is given on the uncertainty of b th . On the other hand, to the best of the authors' knowledge, no information is available on the uncertainty of the compressibility coefficient. Therefore, a variation of ±10% is considered for both the parameters, in order to assess their weight independently from their different experimental uncertainties. Case II-a (compressible fuel mixture with calculated bubble distribution) is selected for this analysis. In fact, this case is the most interesting and complete, involving both compressibility and two-phase flow, as well as their coupling with neutronics.
The power transients obtained by variation of one parameter at a time are shown in Fig. A.1, while the corresponding power peak values are listed in Table A.1.
Results point out that the peak power increases with the compressibility coefficient. In fact, if compressibility increases, the speed of sound decreases, leading to a larger delay of the thermal expansion feedback and, as a consequence, to a higher energy release. On the other hand, concerning the thermal expansion coefficient, two opposite effects come into play: 1. A higher b th causes a stronger expansion feedback on neutronics, leading to a lower energy release; 2. A higher b th triggers stronger pressure waves and, as a consequence, larger fuel compressions, leading to a higher energy release.
Results point out that the first effect prevails on the second, since the power peak decreases with b th . More in general, it can be observed that transients are slightly more sensitive to b th then w, even if sensitivity to both parameters is quite small (±10% perturbations of the coefficients lead to variations of the peak power values that are always smaller than 1.3%).

Appendix B. Temperature fields
In this Appendix, the temperature fields of all the considered case studies are presented (see Figs. B.1-B.5).