© C. Tripodo et al., published by EDP Sciences, 2019

1 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 fourth-generation nuclear reactors known as Generation-IV International Forum [1]. The demonstration of the load-following 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 object-oriented 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 component-oriented 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 one-dimensional 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.

2 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].

Fig. 1 MSFR plant conceptual scheme.

2.1 Reactor fuel circuit and core

The main conceptual feature that distinguish the MSFR is the nuclear fuel that is dissolved in a liquid fluoride- or 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³, 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.

Fig. 2 MSFR fuel circuit conceptual scheme.

Fig. 3 MSFR fuel circuit layout.

2.2 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-pressure-related 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].

3 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 real-time 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).

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 β_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 neutron-importance 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.(1) (2) (3) (4) (5)

The RHS source terms and 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 ψ 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 and decay heat concentration are computed from equations (7) and (8).(6) (7) (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.(9) (10)

The total reactivity − equation (11) − is the sum of the externally inserted reactivity δρ_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 neutron-importance-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).(11) (12) (13) (14) (15)

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.(16) (17) (18)

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.(19)

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.

Fig. 4 Conceptual scheme adopted for the MSFR neutron-kinetics.

Fig. 5 Object-oriented Modelica model of the fuel circuit.

Fig. 6 Fuel circuit geometry.

3.2 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:(20) (21) (22)

The geometrical parameters that can be specified in the component are the length L, the cross-section area A and the heat transfer perimeter ω, which for a PCHE are expressed as(23) (24)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. One-dimensional finite-volume discretization with counter-current 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.(25) (26) (27)

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 semi-circular ducts [22] were implemented. In equations, it reads:(28) (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 (30) (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.

Fig. 7 Object-oriented Modelica model of the IHX.

3.3 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 η_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 Δp_v and the bypass mass flow rate G_bypass:(32)where K_v is a hydraulic conductance parameter set to 10⁻² 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 zero-order 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⁴ and Re_D ≈1.5 × 10⁵, 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.

Fig. 8 Object-oriented Modelica model of the intermediate loop.

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

Fig. 9 Object-oriented Modelica model of the ECS.

3.4.1 Compressor

The compressor was modelled by considering an energy balance. Since Helium can be considered as perfect gas, the following relations hold:(33) (34)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, η_c is the compressor efficiency, T_iso,c is the isentropic outlet temperature of the compressor and γ 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).

3.4.2 Turbine

The turbine was modelled by considering an energy balance similar to that used in compressor component. In particular,(35) (36)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, η_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.

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

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

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

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

Fig. 10 Object-oriented Dymola model of the coupled plant simulator.

4 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 β_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. (37)

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 τ_out. It was assumed that, in the out-of-core 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:(38)where(39) (40)where with υ(r) is indicated the velocity profile shape, which for uniform velocity is υ(r) = 1 and for parabolic velocity is υ(r) = 2(1 − r ²/R ²). 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 β_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. (41) (42) (43)

Figure 11 shows the comparison between the trend against the fuel salt mass flow rate of the steady-state effective delayed fractions β_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 β_eff,stat and those predicted by a lumped-parameter model [25] according to the expression:(44)

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 β_eff,g for high flow velocity, when the fuel circuit recirculation time is far less than the DNPs time constants − 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 β_eff,tot as predicted by the different models. An important fact that can be deduced from the figures is that, in nominal conditions (Γ ≈ 30 t/s, τ_recirc ≈ 4 s), β_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 β_eff,tot is equivalent to a reactivity insertion, this fact has an impact on reactor dynamic behavior during transients (see Sect. 5.1).

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.

Fig. 12 Total effective delayed neutron fraction versus fuel mass flow rate and fuel circuit recirculation time as predicted from the different models.

5 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 full-power 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).

5.1 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 ΔT. 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 β_eff due to the reduction of the fuels salt velocity. The final value of δρ_tot of about −4 pcm (Fig. 13b) therefore correspond to −Δβ_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 β _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 δρ = Δβ_eff. For a zero-power reactor, i.e. neglecting feedback effects of temperature and density on the reactivity, this δρ > 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 β_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).

Fig. 13 Controlled variables evolution for 20% decrease of the fuel salt mass flow rate: (a) reactor thermal power; (b) net reactivity variation; (c) fuel inlet temperature; (d) fuel outlet temperature; (e) fuel average temperature; (f) heat transfer rates in the heat exchangers; (g) intermediate salt minimum temperature; (h) intermediate salt maximum temperature; (i) turbine mechanical power; (m) helium turbine inlet temperature.

5.2 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).

Fig. 14 Controlled variables evolution for 20% decrease of the intermediate salt mass flow rate: (a) reactor thermal power; (b) net reactivity variation; (c) fuel inlet temperature; (d) fuel outlet temperature; (e) fuel average temperature; (f) heat transfer rates in the heat exchangers; (g) intermediate salt minimum temperature; (h) intermediate salt maximum temperature; (i) turbine mechanical power; (m) helium turbine inlet temperature.

5.3 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 post-transient steady-state power level, suggesting that the helium flow rate might be a suitable control variable for reactor power regulation.

Fig. 15 Controlled variables evolution for 20% increase of the turbine helium mass flow rate: (a) reactor thermal power; (b) net reactivity variation; (c) fuel inlet temperature; (d) fuel outlet temperature; (e) fuel average temperature; (f) heat transfer rates in the heat exchangers; (g) intermediate salt minimum temperature; (h) intermediate salt maximum temperature; (i) turbine mechanical power; (m) helium turbine inlet temperature.

5.4 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 ΔT, 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.

Fig. 16 Controlled variables evolution for 0.1 $ external reactivity step insertion: (a) reactor thermal power; (b) net reactivity variation; (c) fuel inlet temperature; (d) fuel outlet temperature; (e) fuel average temperature; (f) heat transfer rates in the heat exchangers; (g) intermediate salt minimum temperature; (h) intermediate salt maximum temperature; (i) turbine mechanical power; (m) helium turbine inlet temperature.

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

• Settling times to post-transient equilibrium values are of the order of 100–150 s for all the considered transients. 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.

6 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 object-oriented 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 reactor-follows-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.

Author contribution statement

C. Tripodo developed the core and fuel circuit models, and implemented the power plant simulator. A. Di Ronco developed the intermediate loop and power conversion loop models. C. Tripodo and S. Lorenzi carried out the numerical simulations, and analyzed and interpreted the results. C. Tripodo wrote the manuscript with input and feedback from S. Lorenzi. A. Cammi and S. Lorenzi conceived the presented idea, supervised the work and offered assistance throughout the project.

Disclaimer

The content of this paper does not reflect the official opinion of the European Union. Responsibility for the information and/or views expressed therein lies entirely with the authors.

Nomenclature

Acronyms

ALFRED: Advanced Lead Fast Reactor European Demonstrator

BoP: Balance of Plant

DNP: Delayed Neutron Precursor

FP: Fission Product

GIF-IV: Generation-IV International Forum

IHX: Intermediate Heat Exchanger

MSFR: Molten Salt Fast Reactor

PCHE: Printed Circuit Heat Exchanger

SAMOFAR: Safety Assessment of the Molten Salt Fast Reactor

SHX: Secondary Heat Exchanger

Latin symbols

A : Cross sectional area (m²)

C_f : Friction coefficient (−)

c_g : Normalized precursor density for gth group (−)

C_g : Precursor density for gth group (m⁻³)

cmd: Bypass valve command signal (−)

d : Density (kg m⁻³)

d_ch : Channel diameter of heat exchanger (m)

D_hyd : Hydraulic equivalent diameter (m)

f : Fanning friction factor (−)

f_Darcy : Darcy friction factor (−)

f_k : Decay heat fraction for kth group (−)

F_k : Decay heat density for kth group (W m⁻³)

g : Gravitational acceleration (m s⁻²)

G_bypass : Bypass mass flow rate (kg s⁻¹)

h : Specific enthalpy (J kg⁻¹)

h_cold : Cold side convective heat transfer coefficient (W m⁻² K⁻¹)

h_hot : Hot side convective heat transfer coefficient (W m⁻² K⁻¹)

H : Core height (m)

H_e : Core extrapolated height (m)

J₀ : Zero order Bessel function (−)

k : Thermal conductivity (W m⁻¹ K⁻¹)

K_v : Hydraulic conductance parameter (kg s⁻¹ Pa⁻¹)

L_e : Core extrapolated length (m)

n₀ : Neutron density at core center (m⁻³)

n_fiss : Normalized fission power (−)

Nu: Nusselt number (−)

p : Pressure (Pa)

P₀ : Nominal reactor power (W)

Pr: Prandtl number (−)

q ^″ : Heat flux (W m⁻²)

q ^‴ : Power density (W m⁻³)

Q : Power (W)

r : Radial coordinate (m)

R : Core radius (m)

R_e : Core extrapolated radius (m)

Re: Reynolds number (−)

s : HX wall thiness (m)

S : External neutron source (s⁻¹)

t : Time (s)

T : Temperature (K)

u : Axial velocity (m s⁻¹)

w : Mass flow rate (kg s⁻¹)

w_f : Average fission energy (J fiss⁻¹)

x : Coordinate along fuel flow (m)

z : Vertical coordinate (m)

Greek symbols

α : Reactivity feedback coefficient (K⁻¹)

β : Delayed neutron fraction (−)

γ : Specific heat ratio (−)

Γ: Fuel mass flow rate (kg s⁻¹)

η : Isentropic efficiency (−)

: Angular velocity (rad s⁻¹)

λ_DH,k : Decay heat decay constant for kth group (s⁻¹)

λ_g : DNP decay constant for gth group (s⁻¹)

Λ : Effective neutron lifetime (s)

ν : Emitted neutron per fission (fiss⁻¹)

φ : Neutron flux 2-D profile shape (−)

δρ : Reactivity variation (−)

τ_core : Core transit time (s)

τ_out : Out-of-core transit time (s)

τ_recirc : Fuel circuit recirculation time (s)

υ : Radial velocity profile shape (−)

ψ : Neutron flux axial profile shape (−)

ω : Tube perimeter (m)

Subscripts

0: Starting steady-state value

c : Compressor

cold: Cold side

dens: Density

DH: Decay Heat

eff: Effective

exch: Exchange

ext: External

fiss: Fission

g : DNP group

hot: Hot side

in: Inlet

int: Intermediate

k : Decay heat group

lump: Lumped

m : Metal

out: Outlet

t : Turbine

tot: Total

vol: Volume

w : Wall

Acknowledgments

This project has received funding from the EURATOM research and training programme 2014-2018 under grant agreement No 661891.

References

All Tables

All Figures

	Fig. 1 MSFR plant conceptual scheme.
In the text

	Fig. 2 MSFR fuel circuit conceptual scheme.
In the text

	Fig. 3 MSFR fuel circuit layout.
In the text

	Fig. 4 Conceptual scheme adopted for the MSFR neutron-kinetics.
In the text

	Fig. 5 Object-oriented Modelica model of the fuel circuit.
In the text

	Fig. 6 Fuel circuit geometry.
In the text

	Fig. 7 Object-oriented Modelica model of the IHX.
In the text

	Fig. 8 Object-oriented Modelica model of the intermediate loop.
In the text

	Fig. 9 Object-oriented Modelica model of the ECS.
In the text

	Fig. 10 Object-oriented Dymola model of the coupled plant simulator.
In the text

	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.
In the text

	Fig. 12 Total effective delayed neutron fraction versus fuel mass flow rate and fuel circuit recirculation time as predicted from the different models.
In the text

Fig. 13 Controlled variables evolution for 20% decrease of the fuel salt mass flow rate: (a) reactor thermal power; (b) net reactivity variation; (c) fuel inlet temperature; (d) fuel outlet temperature; (e) fuel average temperature; (f) heat transfer rates in the heat exchangers; (g) intermediate salt minimum temperature; (h) intermediate salt maximum temperature; (i) turbine mechanical power; (m) helium turbine inlet temperature.

In the text

Fig. 14 Controlled variables evolution for 20% decrease of the intermediate salt mass flow rate: (a) reactor thermal power; (b) net reactivity variation; (c) fuel inlet temperature; (d) fuel outlet temperature; (e) fuel average temperature; (f) heat transfer rates in the heat exchangers; (g) intermediate salt minimum temperature; (h) intermediate salt maximum temperature; (i) turbine mechanical power; (m) helium turbine inlet temperature.

In the text

Fig. 15 Controlled variables evolution for 20% increase of the turbine helium mass flow rate: (a) reactor thermal power; (b) net reactivity variation; (c) fuel inlet temperature; (d) fuel outlet temperature; (e) fuel average temperature; (f) heat transfer rates in the heat exchangers; (g) intermediate salt minimum temperature; (h) intermediate salt maximum temperature; (i) turbine mechanical power; (m) helium turbine inlet temperature.

In the text

Fig. 16 Controlled variables evolution for 0.1 $ external reactivity step insertion: (a) reactor thermal power; (b) net reactivity variation; (c) fuel inlet temperature; (d) fuel outlet temperature; (e) fuel average temperature; (f) heat transfer rates in the heat exchangers; (g) intermediate salt minimum temperature; (h) intermediate salt maximum temperature; (i) turbine mechanical power; (m) helium turbine inlet temperature.

In the text