Development of a control-oriented power plant simulator for the molten salt fast reactor

In this paper, modelling and simulation of a control-oriented plant-dynamics tool for the molten salt fast reactor (MSFR) is presented. The objective was to develop a simulation tool aimed at investigating the plant response to standard control transients, in order to support the system design finalization and the definition of control strategies. The simulator was developed employing the well tested, flexible and open-source objectoriented Modelica language. A one-dimensional modelling approach was used for thermal-hydraulics and heat transfer. Standard and validated thermal-hydraulic Modelica libraries were employed for various plant components (tubes, pumps, turbines, etc.). An effort was spent in developing a new MSR library modelling the 1D flow of a liquid nuclear fuel, including an ad-hoc neutron-kinetics model which properly takes into consideration the motion of the Delayed Neutron Precursors along the fuel circuit and the consequent reactivity insertion due to the variation of the effective delayed fractions. An analytical steady-state 2-D model of the core and the fuel circuit was developed using MATLAB in order to validate the Decay Neutron Precursors model implemented in the plant simulator. The plant simulator was then employed to investigate the plant dynamics in response to three transients (variation of fuel flow rate, intermediate flow rate and turbine gas flow rate) that are relevant to control purposes. Simulation outcomes highlight the typical reactor-follows-turbine behavior of the MSFR, and they show the small influence of fuel and intermediate flow rate on the reactor power and their strong effects on the temperatures in their respective circuits. Starting from the insights on the reactor behavior gained from the analysis of its free dynamics, the plant simulator here developed will provide a valuable tool in support to the finalization of the design phase, the definition of control strategies and the identification of controlled operational procedures for reactor startup and shutdown.


Introduction
The objective of this work was to develop a fast-running, control-oriented plant-dynamics simulation tool for the molten salt fast reactor (MSFR) and the associated Balance of Plant, and to use it to investigate and analyze the plant dynamics.
The MSFR is the circulating-fuel fast-neutron-spectrum reactor concept currently object of research under the EU SAMOFAR project (http://samofar.eu/), within the international framework for the development of fourthgeneration nuclear reactors known as Generation-IV International Forum [1]. The demonstration of the loadfollowing capabilities and the control operability of the reactor is one of the objectives of the SAMOFAR project. In this view, it is important to rely on a power plant simulator to study the system dynamics and to develop and test the control strategies. Due to the dynamic and control related purposes of the power plant simulator, an objectoriented modelling approach is selected as suitable choice for the model-based control design. Due to its features in terms of hierarchical structure, abstraction and encapsulation, this approach allows developing a model that satisfies the requirements of modularity, openness and efficiency [2]. A viable path to achieve the above-mentioned goals is constituted by the adoption of the Modelica language [3]. Modelica is an object-oriented, declarative, equation-based language developed for the componentoriented modelling of complex physical and engineering systems [2]. It allows a description of single system components (or objects) directly in terms of physical equations and principles, and to connect different components through standardized interfaces (or connectors). In addition, his acausal component-based modelling strategy provides high reusability of the models and flexibility of the plant configuration, as well as a more realistic description of the plant, since several validated libraries of power plant components exist (e.g. the ThermoPower library [4]). Modelica is open-source and it has already been successfully adopted in different fields, such as automotive, robotics, thermo-hydraulic and mechatronic systems, but also in nuclear simulation field. Plant simulators were developed for control purposes for the ALFRED (Advanced Lead Fast Reactor European Demonstrator) reactor [5] and the IRIS reactor [6]. As simulation environment, Dymola (Dynamic Modelling Laboratory) [7] was adopted, even if open-source implementation can be considered as alternative option, e.g. OpenModelica [8].
In developing the power plant simulator for the MSFR, it is essential to consider the peculiar features of this reactor, firstly the presence of a liquid circulating nuclear fuel that acts contemporarily as coolant. The strong physical coupling of thermo-fluid-dynamics and neutronics which characterizes the MSFR indeed required to take into account the motion of the Delayed Neutron Precursors (DNPs), which circulate along the fuel circuit. A onedimensional modelling approach was therefore employed for the reactor (as well as for the remaining of the plant) as the best compromise between, on one hand, the need to consider the spatial dependence of the DNPs concentration, and, on the other hand, the need to have a computationally efficient, fast running simulation tool suitable to be employed for plant dynamics investigation and subsequently in support to the design of the plant control system. An ad-hoc point-kinetics model, which is able to take into account the DNPs position in the core, was implemented using a hybrid 0D-1D approach. To verify the DNPs model employed in the MSFR power plant simulator and the corresponding predicted values of the effective delayed neutron fractions for the various delayed groups, an analytical steady-state 2-D model of the reactor core was developed by using the MATLAB ® software [9], under suitable simplifying assumptions. The plant simulator was then employed to investigate the plant free dynamics (i.e., the plant response with no control actions) in response to different transients that are relevant for the development of the control strategy. Four different transients were simulated and analyzed: (i) a reduction of the fuel mass flow rate; (ii) a reduction of the intermediate salt mass flow rate; (iii) an increase of the helium mass flow rate in the turbine unit; and (iv) an external reactivity insertion.
These transients were selected since they involved three of the possible control variables that can be chosen in the control strategy of the reactor for the full power mode, i.e., the operational mode from 50% to 110% of the power. The possibility to control the reactor in this operational mode acting only on the mass flow rates of the different circuit is relevant since the MSFR does not foresee the use of control rods for load-following operation.
The paper is organized as follows. In Section 2, the MSFR reference design is briefly presented, whereas the modelling approach employed for the description of the reactor and the Balance of Plant is described in Section 3. Section 4 illustrates an analytical 2-D benchmark model and its results compared with those of the simulator, in Section 5 the simulator is used to investigate the MSFR plant free dynamics, and in Section 6 some conclusions are drawn.

Reference plant and reactor description
The conceptual scheme of the MSFR BoP is shown in Figure 1. As it can be seen from the figure, the non-nuclear part of the plant consists of a conventional circuit with two loops in series: (i) the Intermediate Loop, through which a fluoride-based molten salt circulates, serves to extract the heat generated in the reactor À through an Intermediate Heat Exchanger (IHX) [10] À and to transport it to the power conversion system; (ii) the Power Conversion Loop, which consists of a conventional Joule-Brayton gas-turbine cycle [11].

Reactor fuel circuit and core
The main conceptual feature that distinguish the MSFR is the nuclear fuel that is dissolved in a liquid fluorideor chloride-based salt which acts contemporarily as fuel and coolant. The reference MSFR design [12] is a 3000 MW th reactor with a total fuel salt volume of 18 m 3 , operated at a mean fuel temperature of 700°C. The reactor schemes are shown in Figures 2 and 3. The fuel circuit, defined as the circuit containing the fuel salt during power generation, includes the core cavity and the recirculation loops (also called 'sectors') including the inlet and outlet pipes, a gas injection system, salt-bubble separators, pumps and fuel heat exchangers. Sixteen cooling sectors are arranged circumferentially around the vessel. Due to the liquid nature of the nuclear fuel, which does not require the presence of any solid fuel-element, and the fast neutron spectrum which does not require any moderating materials, the MSFR core is constituted by a simple, empty cavity, surrounded by an axial reflector and a radial blanket. The fuel salt, with an inlet temperature of about 650°C, enters radially from the bottom into the active zone, where it temporarily reaches criticality and its heated to the outlet temperature of about 750°C. The fuel then exits from the top of the core and it is recirculated through the 16 fuel sectors.

MSFR potentialities
Thanks to its peculiar features, the MSFR presents numerous advantages that make it attractive for the long-term perspective of the nuclear energy. It can operate with very flexible fuel-cycle strategies, reaching high breeding ratios with the thorium cycle, and it is capable to operate as a waste-burner for transuranic waste produced in traditional once-through nuclear reactors, thereby allowing a significant reduction in radiotoxicity [13]. The liquid nature of the fuel allows adjusting on-line the fissile content, with the consequence that no excess reactivity is required in the core at any time to compensate for temperature and power defects, or to compensate fission-products-related reactivity losses. This means that neither burnable poisons nor long-term-adjustment control rods are needed in the core. The continuous removal of fission products allows a better chemical control and allows removing any FPs-related negative reactivity. In particular, the removal of the main nuclear poison Xenon eliminates the reactor dead-time following shutdowns or power reductions, paving the way to much more flexible reactor operation and load-following applications. Great advantages are also present looking at the intrinsic safety aspects of the MSFR. Since the fuel is in a fluid state, the core meltdown scenario is eliminated by-design and no limits exist for the attainable fuel burnup with respect to rods cladding damage and fission gas release. The low vapor tension of the molten salt allows operating the reactor at atmospheric pressure, reducing mechanical stresses on structural components and excluding all high-pressurerelated accidental scenarios. Besides, in case of accidents an emergency fuel-draining system allows to automatically and passively drain the whole fuel content of the reactor, to assure its sub-criticality, and to passively cool it long-term [14]. Finally, the dual fuel/coolant role of the salt, together with its neutronics characteristics, implies that the MSFR has very large, negative, prompt temperature and void reactivity feedback coefficients, making the reactor extremely stable [15].

MSFR plant simulator modelling
In the perspective of identifying effective plant control strategies for an innovative reactor concept like the MSFR, an essential preliminary step was to acquire sufficiently accurate knowledge and understanding of both the reactor system dynamics and the whole Balance of Plant dynamics. To this aim, a control-oriented plant-dynamics simulator was developed and then used to study the MSFR dynamics. A proper dynamic simulation tool for control-oriented purposes, especially in a preliminary design phase, should satisfy some basic requirements [4,5]. In particular it should be modular and extensible, in order to be easily modified and updated to follow the design evolutions; readable, to allow an easy understanding of the equations implemented; computationally efficient, to allow fast-running and realtime simulations; be easily integrable with the control system model.
With the above requirements to be fulfilled, the modelling choice fell on the Modelica language [3]. Modelica is an object-oriented, acausal, equation-based language which offers great advantages in terms of modularity, extensibility, readability and integrability with control-dedicated software (e.g. MATLAB control  toolbox). The simulator was implemented within the Dymola simulation environment [7], which is equipped with state-of-the-art implicit numerical integration algorithms (e.g. DASSL) to handle non-linear differential-algebraic equations sets and with effective homotopy-based model-initialization algorithms [16], and which provides powerful model-linearization tools potentially useful in the future control system design phase. The tested and validated ThermoPower thermal-hydraulic Modelica library [4] has been used for the simulator modelling, and it has been significantly modified and extended into an MSR library to account for all the balance equations pertaining the various nuclear variables (see Sect. 3.1).

Fuel circuit and core
The usual approach employed for dynamics and control in conventional solid-fueled reactors is the so-called Point-Kinetics (i.e., zero-dimensional kinetics) [5,6]. In a circulating-fuel reactor like the MSFR the DNPs move along the fuel circuit, and a proper neutronics modelling needed to take into account the position of emission of the delayed neutrons in the core. Besides, a fraction of the delayed neutrons are emitted in the out-of-core portion of the primary circuit, thereby reducing the effective delayed neutron fraction b eff [17], with a clear impact on the reactor dynamics. An ad-hoc neutronics model, which is able to take into account the DNPs position in the core, was therefore developed using a hybrid 0D-1D approach. Similar approaches have been proposed in previous works on circulating-fuel reactors' dynamics [18]. The conceptual scheme adopted for the fuel circuit modelling is shown in Figure 4.
The circuit thermal-hydraulics determines the spatial distribution of the DNPs concentration along the fuel circuit. The DNPs spatial profile is then used to compute an effective core-averaged value of the DNPs concentration in the core, suitable to be used in the reactor kinetics equation [19]. To correctly account for the drift of the DNPs, i.e., the fact that they are created in a different location with respect to the emission of the corresponding delayed neutron, in the averaging procedure the delayed neutron source intensity can be weighted with a neutronimportance function that can be both the direct flux or more properly the adjoint neutron flux [20]. Similarly, the average temperature used for the feedback reactivity evaluation is obtained as weighted-average of the temperature profile in the core multiplied by the importance function. The decay heat distribution was modelled using the same 1-D modelling approach. The total reactor power is the sum of the fission power in the core and the decay power throughout the whole fuel circuit.
The Modelica model of the fuel circuit is shown in Figure 5. The thermal-hydraulics of the reactor core was modelled in the MSR_Core component (Fig. 5). It is described by the mass (Eq. (1)), X-momentum (Eq. (2)), energy (Eq. (3)) conservation equation and the balance for the DNPs concentration for the 8 DNPs groups (Eq. (4)) and Decay Heat (DH) concentration for the 3 decay-heat groups (Eq. (5)). In the last three equations, the generation term due to the fission process is included. Longitudinal heat and species diffusion were neglected.
The RHS source terms q 000 fiss and q 000 DH in equation (3) are the fission power density and the decay-heat generation density, respectively. The friction coefficient C f appearing in the momentum equation (2) is evaluated using the Colebrook hydraulic correlation [21].
The term c is the neutron importance function and it was assumed to be fixed and equal to the fundamental eigenfunction of the single-energy diffusion theory for bare uniform reactor À i.e., a sinusoidal profile with a proper extrapolation length (Eq. (6)). The values of the fission heat concentration q 000 fiss and decay heat concentration q 000 DH are computed from equations (7) and (8).
The time evolution of the normalized core fission power n fiss (t) = Q fiss (t)/Q fiss,0 is determined in the Neutron_ Kinetics component by the reactor-kinetics equation (9), in which the effective, neutron-importance-weighted averages of the DNPs concentrations À equation (10) À are used, noting that, in the single-energy diffusion theory approximation, the neutron-importance function is taken as the neutron flux profile.
The total reactivity À equation (11) À is the sum of the externally inserted reactivity dr ext and the feedback reactivity of fuel salt temperature and density. The latter two are determined by equations (12) and (13), where the effective temperature T eff is determined as a neutronimportance-weighted core average À equation (14) À and T eff,0 is the reference temperature with respect to which the reactivity defects are calculated. The effective delayed neutron fractions, which take into account the spatial distributions of the DNPs and the importance of the emitted neutrons, are evaluated according to equation (15).
dr T ðtÞ ¼ a T ½T eff ðtÞ À T eff;0 ð 12Þ dr dens ðtÞ ¼ a dens ½T eff ðtÞ À T eff;0 ð 13Þ The 16 external loops forming the fuel circuit were modelled as a single equivalent loop formed by a hot leg section, representing the piping from the core outlet to the IHX inlet, the IHX and a cold leg section representing the piping from the IHX outlet to the core inlet (Fig. 6). The HotLeg and Cold leg tube components implement the single-state, one-dimensional, finite-volume-discretized conservation equations for mass (Eq. (1)) and momentum (Eq. (2)), whereas energy (Eq. (16)), DNPs (Eq. (17)), and DH (Eq. (18)) equations are modified to consider only the decay term.
Ideal, mass-flow-rate-controlled pumps (PumpFuel component) establish the salt flow through the circuit.
The reactor total power is the sum of the fission power in the core and the decay heat generated along the whole fuel circuit À equation (19). Reactor geometrical, operational, physical and neutronic data used in the following are shown in Tables 1 and 2 (with reference to Fig. 6). All the parameters of the simulator are easily modifiable at runtime to allow for model modification and update throughout the various design phases.
Since the fuel circuit forms a closed loop, it was important to provide an expansion tank to avoid strong pressure variations caused by temperature transients. The SinkPressure component allows handling any mass insurge or outsurge transient, with no associated dynamic effect. When mass flows from the sink to the loop, the outsurge fluid was assumed to be at the same temperature of the cold leg.

Intermediate heat-exchanger (IHX)
Due to its non-conventional design, an effort was spent to set up a specific component representing the MSFR intermediate heat exchangers [10]. The heat exchangers were modelled as counterflow heat exchangers, with particular reference to the Printed Circuit Heat Exchanger À a proposed technology for the MSFR, for more detail see [11] À but any other counterflow arrangement based on parallel flow pipes subject to heat transfer through their lateral surface can be modelled as well with little modification. The Intermediate_HX model (Fig. 5) is based on components from the ThermoPower library, especially the Flow1DFV component, which describes the fluid flow in a rigid tube. It is based on a 1D finite volume discretization of the mass (Eq. (20)), momentum (Eq. (21)) and energy transport (Eq. (22)) equations: The geometrical parameters that can be specified in the component are the length L, the cross-section area A and the heat transfer perimeter v, which for a PCHE are expressed as where d ch is the channel diameter. Figure 7 shows the Modelica model of the IHX whereas the geometric and operational parameters are shown in Table 3. Onedimensional finite-volume discretization with countercurrent flows was employed for the heat transfer in the heat exchanger. A single, equivalent heat exchanger component, representative of the 16 parallel ones (one for each of the parallel fuel circuit loops), was used. Longitudinal heat transfer along the flow direction was neglected, while the heat capacity of the metal walls of the heat exchanger was accounted for. Equations (25) and (26) are the heat exchange equations on the hot (fuel salt) and cold (intermediate salt) sides, respectively. Equation (27) is the energy balance equation for the IHX metal wall.
Due to the small channel size, the resulting flow is laminar in most of the cases for the fuel salt side. This simplifies considerably the heat transfer modelling (even if it restricts the heat transfer coefficients to quite low values). The average Fanning friction factor (Eq. (28)) and Nusselt number (Eq. (29)) for fully developed laminar flow in semicircular ducts [22] were implemented. In equations, it reads: Nu ¼ 4:089: ð29Þ On the cold (intermediate salt) side the flow regime is in the transition zone (Re ≈ 5000 Ä 7000), and the Gnielinski [21] correlation (Eq. (30)) is used. f Darcy (Re) is the Darcy friction factor, for which the Petukhov [21] correlation for smooth tubes (Eq. (31)) is used f Darcy ¼ ð0:79 lnðReÞ À 1:64Þ À2 : ð31Þ The thermohydraulic correlations to be used in the IHX component are selectable at runtime, to allow for design variations in geometrical and/or operational IHX parameters.

Intermediate loop and secondary heat-exchanger (SHX)
The four intermediate loops were modelled as a single equivalent loop formed by a hot leg section, representing the piping from the IHX outlet to the SHX inlet, a bypass line and a cold leg section representing the piping from the SHX outlet to the IHX inlet (Fig. 8). The intermediate loop model was assembled by using standard components from the ThermoPower library. The adopted scheme is represented in Figure 8. The two basic components are the hotLeg and coldLeg components, which are modelled by Flow1DFV objects. The transport delay associated with the hot/cold leg, with the geometrical parameters indicated in Table 4, is of the order of some seconds. In addition, the dynamic effect associated with thermal inertia is not negligible, hence the total volume of the intermediate loop has a significant influence on dynamics.
The loop includes two active components, a pump and a bypass valve (Fig. 8). The pump class models a simple centrifugal pump with no energy or momentum dynamics and the power consumption was simply estimated through a constant pump efficiency h p . The pump has an external input port which can be used to control the rotational speed and thus the mass flow rate. The valve component was modelled by the ValveLin class, which simply provides a linear constitutive equation to relate the pressure drop Dp v and the bypass mass flow rate G bypass : where K v is a hydraulic conductance parameter set to 10 À2 and cmd is the command signal, provided by an external input port. The valve can be used to control the mass flow rate flowing in the secondary heat exchanger, providing another way to control heat extracted from the intermediate loop. As explained in Section 3.1 for the fuel circuit, an expansion tank was provided to avoid the strong pressure variations related to temperature transients of an incompressible liquid in closed loop and to establish a reference pressure level in the cold leg (1 bar). The expansion tank was modelled using the expansionTank component (Fig. 8).
The other components appearing in Figure 8 are simple temperature and mass flow rate sensors, which model zeroorder sensors providing ideal measurements. Geometric and operational parameters of the intermediate loop are shown in Table 4.
In the SHX, heat is transferred from the intermediate salt to the helium in the Energy Conversion System (ECS). The modelling approach employed for the SHX was identical to that used for the IHX (see Sect. 3.2). Geometric and operational parameters are shown in Table 4. The flow regime in the SHX is fully turbulent on both the salt and gas sides (Re D ≈ 4 Â 10 4 and Re D ≈ 1.5 Â 10 5 , respectively). The Gnielinski [21] correlation is used to evaluate the convective heat transfer coefficients. Also in this case, the correlations to be used in the SHX component are selectable at runtime, to allow for design variations in geometrical and/or operational SHX parameters.

Energy conversion system (ECS)
The energy conversion system model was assembled by using standard components from the ThermoPower library. In particular, a Helium Joule-Brayton cycle with regeneration and three stages of reheating and intercooling was considered. This configuration turned out to ensure a gas temperature at secondary heat exchanger inlet that can avoid salt solidification problem [11]. The adopted scheme is represented in Figure 9.
There are five main components in the model, namely, the compressor, the turbine, the intercooler, the reheater, and the recuperator (Fig. 9). The cycle was modelled as open, i.e., disregarding the final heat sink section. This is a common choice to simplify the modelling of the cycle [5] and it has no impact on the dynamics results since the final sink acts as an infinite heat sink.

Compressor
The compressor was modelled by considering an energy balance. Since Helium can be considered as perfect gas, the following relations hold: where T in,c and p in,c are the gas temperature and pressure at the inlet of the compressor, T out,c and p out,c are the gas temperature and pressure at the outlet of the compressor, h c is the compressor efficiency, T iso,c is the isentropic outlet temperature of the compressor and g is the specific heat ratio of the gas. The efficiency and the pressure ratio can be set by the user in order to adapt the component to the cycle parameters. The component can be connected to a shaft in order to calculate the compressor work (and hence the cycle efficiency).

Turbine
The turbine was modelled by considering an energy balance similar to that used in compressor component. In particular, where T in,t and p in,t are the gas temperature and pressure at the inlet of the turbine, T out,t and p out,t are the gas temperature and pressure at the outlet of the turbine, h t is the turbine efficiency, T iso,t is the isentropic outlet temperature of the turbine. Also, in this case, the efficiency and the pressure ratio are user-selectable parameters and the turbine work can be calculated.

Intercooler
The intercooler is a heat exchanger with the gas and an infinite sink at prescribed temperature (T_intercooling). It is adopted to improve efficiency by decreasing the average specific volume of the gas in the compression stages.

Reheater
The reheater is a heat exchanger with the gas at the outlet turbine and a hotter source. It is adopted to improve efficiency by increasing the average specific volume of the gas in the expansion stages. In the present model, a fraction of the gas at the secondary heat exchanger outlet is extracted to reheat the colder gas at the turbine outlet. An alternative option is to employ the hot intermediate molten salt in the reheaters.

Recuperator
The recuperator is a heat exchanger aimed at performing regeneration, i.e., preheating the gas at the inlet of the secondary heat exchanger with the high temperature gas at the turbine outlet. This improves the efficiency and can also avoid the problem of the salt solidification in the secondary heat exchanger.

Full plant simulator
The MSFR plant simulator model was built assembling the various sub-models illustrated in Sections 3.1 through 3.4. The full, coupled, Modelica model is shown in Figure 10, while Table 5 shows the various model input and output variables. Table 6 summarizes various meaningful plant variables in steady-state nominal operating conditions, as obtained with the present plant simulator.

Analytical benchmark model
An analytical model was developed to verify the DNPs model implemented in the MSFR power plant simulator and the corresponding values of the effective delayed neutron fractions b eff,g for the various delayed groups. The analytical model, implemented in MATLAB, is able to calculate the two-dimensional, steady-state DNPs spatial distribution in the core and the corresponding effective delayed fractions. In order to obtain an analytical solution, several simplifying assumptions were made. A 2-D r-z axial-symmetric cylindrical geometry was assumed for the core, and the neutronic flux shape was assumed to be the single-group diffusion fundamental eigenfunction for bare cylindrical reactors, i.e. a sinusoidal axial dependence and a 0th-order Bessel radial dependence À equation (37). The presence of a reflector was allowed for through the introduction of a proper extrapolation length.
The fuel salt velocity profile was assumed fixed and axially directed, the motion occurs in the form of parallel streamlines and DNPs turbulence mixing and molecular diffusion were neglected. Two different velocity profiles were considered, a uniform one and a parabolic one. The out-of-core portion of the circuit was modelled with a 0-D geometry, i.e., with a simple, mass-flow-rate-dependent, out-of-core time t out . It was assumed that, in the out-ofcore portion of the fuel circuit, complete fluid mixing occurs. The DNPs re-entry boundary condition at core inlet was therefore assumed to be uniform, and equal to the average outlet concentration, reduced by the fraction of DNP which decays in out-of-core portion of the circuit. All physical properties were considered constant and evaluated at the core average temperature. Geometric, operational and physical parameters used in the following are shown in Table 7.
Under the above modelling assumptions, the DNPs balance equation is: where where with y(r) is indicated the velocity profile shape, which for uniform velocity is y(r) = 1 and for parabolic velocity is y(r) = 2(1 À r 2 /R 2 ). Under the hypothesis of parallel streamline flow, equation (38) can be solved analytically, resulting in equation (41) for the DNPs concentrations. The core inlet boundary condition C g,in under the complete mixing assumption is expressed by equation (42). Having obtained the DNPs spatial distributions, the effective delayed neutron fractions b eff,g are calculated taking into account the spatial neutron-importance of the emitted neutrons, according to equation (43) [17], approximated with the direct neutron flux.
See equations (41)-(43) below. Figure 11 shows the comparison between the trend against the fuel salt mass flow rate of the steady-state effective delayed fractions b eff,g as predicted by the Dymola 1-D plant simulator and those predicted by the analytical MATLAB 2-D model for the two velocity profiles. Figure 11 also shows the values of the delayed fractions for static fuel b eff,stat and those predicted by a lumped-parameter model [25] according to the expression: C g ðr; zÞ ¼ C g;in ðrÞe ! 8 > < > : A core e l g t out À f r; z ð Þl g C g r; z ð Þ2prdrdz  Fig. 11. Trend of the effective delayed neutron fractions versus fuel salt mass flow rate as predicted from the different models. Delayed groups 1-8. Results shown in the following refers to parameters values indicated in Table 7. Note that, for reasons of comparison with previous works [26], the values of the static-fuel delayed fractions used in this section for purposes of verification of the DNPs model are slightly different from those employed for plant free dynamics simulation (Tab. 2). As can be seen from Figure 11, the MSFR plant simulator predicts with sufficient accuracy the trend of the effective delayed fractions as a function of the fuel salt mass flow rate, when compared to the analytical 2-D model, in particular for the uniform velocity profile (for the parabolic velocity profile, the MATLAB model underestimates the effective delayed fractions in most of the mass flow rate range: the reason is the higher fuel salt velocity near the core axis, which leads to lower DNP concentrations in a region of high neutron importance which are not compensated by the higher concentrations in the annular region due to the lower neutron importance of this zone). A small offset of about 10% is present for the saturation values À i.e., the limiting values of the b eff,g for high flow velocity, when the fuel circuit recirculation time is far less than the DNPs time constants t recirc ≪ l À1 g À due to the different modelling dimensionalities [18]. For small mass flow rates, the effective delayed fractions tend to their corresponding static-fuel values as can be seen for delayed groups 6 to 8 (the "fastest" groups). For the slowest groups (1 to 4), the static-fuel values are only reached at very low mass flow rates, out of the simulated range of Figure 11. For these groups, in the simulated mass flow rate range the effective delayed fractions are essentially constant and equal to the corresponding saturation values. Figure 12 shows the dependency on the fuel mass flow rate of the total delayed fraction b eff,tot as predicted by the different models. An important fact that can be deduced from the figures is that, in nominal conditions (G ≈ 30 t/s, t recirc ≈ 4 s), b eff,tot is almost insensitive to small mass flow rate variations, due to the fact that groups 1 through 5 are essentially at saturation. Since a variation of b eff,tot is equivalent to a reactivity insertion, this fact has an impact on reactor dynamic behavior during transients (see Sect. 5.1).

Plant free dynamics simulation
The study of the plant free dynamics is a fundamental step in order to understand the behavior of a reactor in response to different transient initiators. It also provides the necessary insights on the plant to support the definition of suitable control strategies and operating procedures. In this section, the simulation results for four different transients are presented and analyzed: (i) fuel salt mass flow rate reduction, (ii) intermediate salt mass flow rate reduction, (iii) turbine helium mass flow rate increase, and (iv) external reactivity insertion. All the transients were simulated starting from nominal fullpower steady-state operating conditions. The plant simulator developed in the present work allows simulating transients of 200 s with computational times of less than 15 s (2.80 GHz with 16 GB memory laptop).

Reduction of the fuel salt mass flow rate
The plant response to a transient consisting of a 20% exponential reduction of the fuel salt mass flow rate was investigated. An exponential time constant of 5 s was chosen to take into account the pumps' inertia. Simulation results are shown in Figure 13. As soon as the flow reduction begins, the slowing down of the fuel salt in the core immediately leads to an increase in its average temperature (Fig. 13e), and therefore to a prompt negative reactivity insertion (Fig. 13b), which leads to a core power reduction (Fig. 13a). Meanwhile, the fuel salt in the IHX experiences a temperature decrease, due to the reduced mass flow rate. When this cooled salt enters the core (Fig. 13c) after a delay corresponding to the cold leg transit time, it provides positive reactivity, as we can see from the small relative peak at about 8 s. The reduced fuel flow rate leads to a decrease in the heat transfer rate in the IHX  ( Fig. 13f), and the intermediate salt outlet temperature start to decrease (Fig. 13h). The reduced intermediate salt temperature then causes a corresponding decrease in the temperature of the helium at the turbine admission (Fig. 13m), which leads to a reduction of the turbine mechanical power output (Fig. 13i). When the transient initiator ends, the tradeoff between the increased fuel heating in the core and the increased fuel cooling in the IHX leads to a new equilibrium with a larger DT. At the end of the transient, the new average core temperature will be such that it exactly compensates the variation in the effective delayed neutron precursor fraction b eff due to the reduction of the fuels salt velocity. The final value of dr tot of about À4 pcm (Fig. 13b) therefore correspond to ÀDb eff . This transient feature is peculiar to circulating fuel reactors and of the MSFR in particular [25]. From any starting condition, a primary circuit flow velocity reduction causes an increase in the effective delayed neutron fraction b eff , because less delayed neutrons precursors decayed outside in the core, and the delayed neutron are created in core positions of higher neutron importance (Sect. 4, [17]). This corresponds to a positive reactivity insertion dr = Db eff . For a zero-power reactor, i.e. neglecting feedback effects of temperature and density on the reactivity, this dr > 0 causes reactor power to start increasing according to the inhour equation. If temperature feedback effects are included in the analysis, core temperature variations cause reactivity insertions, and, depending on the relative magnitudes of the negative "temperature-reactivity" and the positive "precursor-reactivity", the reactor power might either decrease or increase. The exact features of the power evolution during the transient will depend on the starting power level (which determines the magnitudes of the temperature variations), the starting mass flow rate (which determines the b eff dependency on the mass flow rate À see Sect. 4), and on the operating parameters of the heat exchangers and of the energy conversion system. Starting from nominal full-power operating condition, it is clear from Figure 13 that a variation of the fuel mass flow rate has an almost negligible impact on the reactor power level, with power variation below 1% at the end of the transient (Fig. 13a). The small differences (a few MWs, see Figs. 13a and 13f) between the reactor thermal power and the heat transfer rates in the two heat exchangers are due to the pumping powers in the fuel and intermediate circuits. It is instead evident the strong impact it has on the minimum and maximum temperatures of the fuel salt, with variations of about 12°C in opposite directions (Figs. 13c and 13d).

Reduction of the intermediate salt mass flow rate
A 20% intermediate mass flow rate reduction is considered. An exponential flow reduction with a time constant of 5 s was chosen to take into account the pumps' inertia. Simulation results are shown in Figure 14. As soon as the flow reduction begins, the reduced flow rate causes an increase of the intermediate salt maximum temperature, at the IHX outlet (Fig. 14h), and a decrease of the minimum temperature, at the outlet of the SHX (Fig. 14g), due to the increased residence time in the two heat exchangers. The reduced heat transfer rates ( Fig. 14f), due to the reduced velocity, cause core inlet temperature to increase (less power is removed, Fig. 14c), and gas temperature to decrease (less power is ceded, Fig. 14m). When the hotter fuel salt starts filling the core (Fig. 14e), the negative reactivity it provides (Fig. 14b) starts reducing the core power (Fig. 14a). Core inlet temperature reaches a maximum of about +5°C, and then decreases when the cooler intermediate salt reaches the IHX after its transit time. This maximum corresponds to the max negative reactivity value, and in turn to the lowest reactor power, at about t = 18 s. Due to this power minimum, the core outlet temperature reaches in turn a minimum at about t = 25 s. After all temperatures have reached corresponding local maxima/minima at t = 18-25 s, the system slowly stabilizes to its final new equilibrium state when the power extracted from the SHX matches the reduced reactor power (plus the pumping powers). The mechanical power of the turbine unit correspondingly decreases due to the reduction of the helium admission temperature. Intermediate flow rate has a small impact on reactor power and fuel salt temperatures, while it strongly influences the temperatures of the intermediate salt, especially its minimum value (Figs. 14g and 14h).

Increase of the helium turbine mass flow rate
The third simulated transient consists of a 20% increase of the helium mass flow rate in the turbine unit. An exponential flow variation with a time constant of 5 s was considered to take into account the compressors' inertia. Simulation results are shown in Figure 15. As the transient begins, the increased gas flow rate causes an increase in the heat disposal from the SHX, due to the increased heat transfer coefficient at higher gas velocities and to the reduction of the average helium temperature in the SHX (Fig. 15f). For constant reheating gas flow rate, an increase in the flow rate from the compressors' line causes a reduction of the helium temperature at turbine admission (Fig. 15m). This increased heat transfer rate causes a rapid cooling of the intermediate salt exiting the SHX (Fig. 15g), and when this cooled salt reaches the IHX after its transit delay time, it starts to extract more heat from the fuel salt, whose temperature decreases (Fig. 15c). When, in turn, this cooled fuel salt begins to fill the core, it injects positive reactivity (Fig. 15b), determining a corresponding increase of the reactor power (Fig. 15a). This, in turn, causes the fuel outlet temperature to increase (Fig. 15d), and when this hot salt reaches the IHX, the heat transfer to the intermediate salt increases (Fig. 15f), limiting the reduction of the intermediate salt maximum temperature (Fig. 15h) with respect to its minimum temperature (Fig. 15g). At about 140 seconds the system reaches a new equilibrium steady-state, characterized by a higher reactor power level that matches (considering the pumping powers) the augmented heat exchange in the SHX at increased helium flow rate. This transient shows large variations for all the plant's relevant temperatures, and, more importantly, a large variation in the final posttransient steady-state power level, suggesting that the helium flow rate might be a suitable control variable for reactor power regulation.

External reactivity insertion
The fourth and last simulated transient consists of a 0.1 $ external reactivity step insertion. Simulation results are shown in Figure 16. As can be seen from Figure 16a, the reactor thermal power undergoes a prompt jump of about 20%. The average fuel temperature immediately increases of about 6°C in fractions of a second (Fig. 16e) due to the higher fission power, leading to a fast injection of negative reactivity thanks to the large negative prompt feedback coefficient (−1.46 pcm/K from Doppler feedback, −2.91 pcm/K from density feedback) and the power rapidly decreases from the peak prompt-jump value. When the heated fuel salt (Fig. 16d) reaches the IHX, the heat transfer rate in the IHX increases (Fig. 16f) due to the higher DT, and the intermediate salt maximum temperature starts increasing as a consequence. Similarly, when the hotter intermediate salt enters the SHX at about t = 9 s, the power extracted from the SHX increases and the helium temperature at turbine admission increases in turn. For a constant reheating gas mass flow rate (see Sect. 3.4), the average helium temperature in the SHX also increases. When the hot fuel salt produced in the power peak re-enters the core, the average fuel temperature has a slight increase that causes a second, sudden power decrease.
All the plant relevant temperatures are "dragged up" by the augmented fuel temperature. The system simply stabilizes at a higher temperature level, determined by the effective core temperature at which the feedback effects exactly compensate the external reactivity. In particular, for a 0.1 $ external reactivity insertion and a total feedback coefficient of À4.37 pcm/K, this corresponds to about 11.5°C temperature increase in the core, other temperatures being determined by heat transfer in the HXs. Since all the temperatures shift upward following the fuel temperature, the final post-transient value of the latter can be reached at a reactor power level that is very similar to the starting value: a 0.1 $ insertion causes a variation of less than 2% in the reactor thermal power (Fig. 16a) and of about 3% in the turbine mechanical power (Fig. 16i) (the discrepancy is due to the varying thermodynamic efficiency of the cycle at varying temperatures). The practical consequence is that external reactivity is not a suitable input variable for power regulation. On the other hand, it has a strong impact on all plant temperatures.

Outcomes of the MSFR dynamic simulations
From the results of the simulated transients, valuable insights in the reactor behavior, useful for the definition of the normal operating procedures and the selection of reactor control strategies, can be obtained: -The plant dynamics analysis for a reactivity insertion transient clearly shows that an insertion of external reactivity has a small effect on the reactor power. In this perspective, an external reactivity system is not required for the load-following capabilities in terms of power variation. -The MSFR shows a typical reactor-follows-turbine behavior in which is the ultimate heat loop extraction (i.e., the energy conversion system) that drives the core power. This confirms the need to include the modelling of this loop in a power plant simulator aimed at studying the control strategies. In this view, the most suitable candidate for controlling the reactor power output seems to be the helium mass flow rate. -The fuel mass flow rate has small effect on the power. This comes from the design choice of using the printed circuit heat exchangers for the intermediate HX.
In particular, the laminar condition imposed by this type of HX does not allow strong variation in the heat transfer properties when changing the fuel velocity. The fuel mass flow rate variation has a remarkable impact on the inlet and outlet core fuel temperature. The dynamics of the plant is mainly governed by the heat capacity in the HXs rather than physical recirculation time. The characteristics time of the MSFR is strongly influenced by the choice of the metal HX material and its configuration, and different design options for the heat exchangers will lead to different dynamics feature of the reactor.

Conclusions
In this paper, a plant-dynamics simulator oriented to the control design of the MSFR power plant was developed by employing the well tested, flexible and open-source objectoriented Modelica language. Components from validated thermal-hydraulic libraries were used to model the intermediate circuit, the energy conversion system and the heat exchangers, and a new library was created to model the 1-D flow of a liquid nuclear fuel, with the associated finite-volume-discretized balance equations for the DNPs concentration and the decay heat density. In particular, an effort was spent in implementing an ad-hoc hybrid 0D-1D neutron-kinetics which properly takes into account the position of the DNPs along the fuel circuit and the consequent reactivity insertion due to the variation of the effective delayed fractions. In addition, an analytical steady-state 2-D model of the core and the fuel circuit was developed using MATLAB in order to verify the DNPs model. The simulator was then employed to investigate the MSFR power plant free dynamics in response to four typical design-basis transient initiators. Computational times for all the considered transients are of the order of a few seconds, proving the simulator to be a very computationally-efficient tool.
Simulation outcomes highlight the typical reactorfollows-turbine behavior, in which the ultimate heat extraction loop (i.e., the energy conversion system) determines the evolution of the reactor power, confirming the MSFR load-following capabilities without requiring an external reactivity insertion system. The mass flow rates in the fuel and intermediate circuits are shown to have small effects on reactor power, whereas they strongly influence the salt temperatures in the respective circuits. The plant dynamics characteristic times are mainly due to the characteristics of the heat exchangers. Starting from the insights on the reactor behavior gained from the analysis of its free dynamics here presented, and thanks to the characteristics of the Modelica language in terms of flexibility and integrability with control design software (e.g. MATLAB Control System Toolbox), the plant simulator here developed will provide a valuable tool in support to the finalization of the design phase and to the definition of model-based plant control strategies. This project has received funding from the EURATOM research and training programme 2014-2018 under grant agreement No 661891.