The Molten Salt Fast Reactor (MSFR) has garnered much interest for its inherent safety and sustainbility features. The MSFR can adopt a closed thorium fuel cycle for sus- tainable operation through the breeding of233U from232Th. The fuel composition changes signiﬁcantly over the course of its lifespan. In this study, we investigated the steady state and transient behavior of the MSFR using Moltres, a coupled neutronics /thermal-hydraulics code developed within the Mul- tiphysics Object Oriented Simulation Environment (MOOSE) framework. Three di erent fuel compositions, start-up, early- life, and equilibrium, were examined for potentially dangerous core temperature excursions during a unprotected loss of heat sink (ULOHS) accident. The six-group and total neutron ﬂux distributions showed good agreement with SERPENT and pub- lished MSFR results, while the temperature distribution and total power showed discrepancies which can be attributed to known sources of error. For the transient behavior under the ULOHS scenario, while the transition time towards the new steady state core temperature is also in good agreement with existing MSFR simulations by Fiorina et al., Moltres under- estimated the temperature rise by a factor of ten, due to the same sources of error a ecting the steady state results. While an MSFR loaded with start-up fuel composition operates at a higher temperature than with the other two fuel compositions, all three cases were shown to be inherently safe due to the strong negative temperature feedback. I INTRODUCTION Molten Salt Reactors (MSRs) are a class of nuclear re- actors that contain nuclear fuel dissolved and circulating in a molten salt coolant loop. They potentially possess the ability to run for extended periods with minimal shutdown time due to online fuel reprocessing. Their equilibrium fuel compositions di er substantially from start-up compositions not only due to burnup of initial ﬁssile material and breeding of new ﬁssile material, but also ﬁssile material feeds and removal of ﬁssion products. Also, while early-life fuel composition is largely dependent on the initial core loading, the long-term equilib- rium fuel composition is determined by the type of feed. Since the changing fuel composition impacts safety parameters (e.g. reactivity feedback coe cients), a licensing case for this class of reactors must fully characterize those impacts. While numerous computational tools exist for conven- tional nuclear reactors, MSRs present unique computational challenges that many fail to address e ectively. MSRs di er profoundly from conventional solid-fuelled reactors, partic-ularly in their neutronics and thermal-hydraulics behaviors. New MSR simulation tools must capture strong coupling be- tween neutronics and thermal-hydraulics exhibited by delayed neutron precursor (DNP) movement as well as strong Doppler and density feedback in the fuel salt. This paper investigates the impact of changing fuel composition on fuel salt tempera- tures in the Molten Salt Fast Reactor (MSFR) concept using a new simulation tool for MSRs, Moltres [1]. Moltres is an open source coupled neutronics /thermal hy- draulics simulation application for simulating MSRs. Built on the Multiphysics Object-Oriented Simulation Environment (MOOSE) ﬁnite element framework [ 2], Moltres solves the coupled time-dependent multi-group neutron di usion, tem- perature, and DNP governing equations. The temperature and DNP equations fully account for ﬂow via fuel salt advection. The MSFR model studied in this paper is a reference de- sign for a fast-spectrum MSR developed under H2020 Safety Assessment of the Molten Salt Fast Reactor (SAMOFAR) project [ 3]. The MSFR boasts several safety and sustainability advantages over conventional reactors. Firstly, it can run on a closed thorium fuel cycle, which reduces actinide production and waste radiotoxity. The fast neutron spectrum improves 233U breeding from232Th, an isotope that is approximately three times as abundant as uranium. Lastly, the MSFR op- erates at near atmospheric pressure; this reduces the risk of containment structure failure and manufacturing costs. This paper presents results from Moltres simulations of the MSFR reference model with three fuel compositions: start- up, early-life, and equilibrium. In line with the purpose of the MSFR as a thorium breeder, the chosen start-up fuel com- position under study is a eutectic mixture of233U and232Th ﬂuorides in a lithium ﬂuoride molten salt, with a Th /233U mixture feed.233U is extracted from the blanket tank and reinserted into the core [4]. We generated group constants for each fuel composition using SERPENT [ 5], a continuous-energy Monte Carlo code for numerous reactor physics applications. The group con- stants relevant for Moltres simulations are the various macro- scopic neutron cross sections, neutron di usion coe cient, average ﬁssion energy, average neutron yield, inverse neutron speed, ﬁssion spectrum, DNP decay constants, and e ective de- layed neutron fractions. Using these group constants, Moltres then solves for the ﬂux and temperature based on the neu- tron di usion equation coupled with temperature advection due to coolant ﬂow. Transient simulations will establish core fuel temperatures during transients. These distributions will give insight into MSFR transient behavior, which will help us identify potential safety risks. Fig. 1: MSFR reactor design concept [6]. TABLE I: Speciﬁcations of the MSFR design [6]. Parameter Value Thermal /Electric output [MW th/MW e] 3000 /1500 Salt volume [m3] 18 Salt fraction in core 0.5 Number of circulation loops 16 Nominal ﬂow rate [kg s1] 18500 Nominal circulation time [s] 4.0 Inlet/outlet temperature [K] 923 /1023 Blanket volume [m3] 7.3 Section II provides the speciﬁcations and a literature re- view of the MSFR concept. Section III covers the overall modeling approach with respect to the codes used in this study: SERPENT and Moltres. Section IV presents steady state neu- tron ﬂux and temperature distribution results from Moltres for the three di erent fuel compositions studied. Section V explores the MSFR’s neutronic and thermodynamic response to a Unprotected Loss of Heat Sink (ULOHS) scenario. Lastly, Section VI provides conclusions and suggestions for future work. II MOLTEN SALT FAST REACTOR Figure 1 shows a schematic view of the MSFR. The main speciﬁcations of the MSFR are given in Table I. Lithium ﬂuoride (LiF) is the major component of the fuel and blanket molten salts used for the MSFR. Fissile and fertile isotopes are introduced into the mixture by mole fractions of 77.5%LiF-22.5%AcF 4, where AcF 4represents actinide ﬂuorides such as uranium and thorium ﬂuorides. Table II shows the mole fraction of the molecular constituents in the MSFR fuel salt, with the relevant physical properties in Table III. The MSFR supports various fuel compositions; it can run on the same235U-238U fuel used in most conventional LWRs. It can also run on a mixture of fresh and used uranium fuel containing reprocessed TRU isotopes [ 7]. However, the main conﬁguration of the MSFR is a breeder reactor running onTABLE II: Fuel salt composition used for the Serpent simulation. Isotope Mole fraction [%] 7Li 29 :0 19F 62 :6 232Th 7 :40 233U 1 :00 TABLE III: Fuel and blanket salt properties. Tdenotes temperature in Kelvins. Property Value Density,[kg m3] 4983 :560:882T Thermal cond. , k[W m1K1] 0:928+8:397105T Speciﬁc heat, cp[J kg1K1]1111 +2:78T TABLE IV: Material compositions and densities of the absorber, and reﬂector regions. All elements listed are in their natural compositions. Region Element At% Density [kg m3] AbsorberB 802520C 20 ReﬂectorNi 67 :7 10000W 25 :0 Cr 7 :0 Al 0 :3 233U-232Th [4]. Previous studies have reported breeding ratios of up to 1.1 on the MSFR [8]. The primary fuel salt ﬂows upwards through the 9 m3 central core region. At the top of the core, the ﬂow separates into 16 smaller external loops, each of which passes through a heat exchanger and is pumped back into the core. The primary heat exchangers transfer heat from the fuel salt to an interme- diate salt coolant loop. There is other instrumentation along the external loops for online reprocessing and gas sparging. The core is radially surrounded by a tank of blanket molten salt, with reﬂectors at the top and bottom of the core. The blanket salt contains fertile isotopes such as232Th for fuel breeding. There is a layer of neutron absorbing material be- hind the blanket tanks to protect the heat exchangers, pumps and other instrumentation from neutron irradiation damage. Various authors have performed a number of steady state and transient multiphysics simulations for the MSFR. A paper by Fiorina et al. [ 9] compares results between models devel- oped on the multiphysics software COMSOL, and coupled in-house codes developed at Delft University of Technology. The models used 2D axisymmetric models of the MSFR and solved multi-group neutron di usion equations for the neu- tronics, and they showed good agreement for some of the transient cases studied. A more comprehensive 3D model has been developed and studied by Auﬁero et al. [ 10] using the computional ﬂuid dynamics (CFD) toolbox OpenFOAM. Although this model relied on the one group neutron di usion equation, it enabled the study of the full 3D core geometry and 3D transient scenarios such as the failure of one of the 16 pumps in the MSFR. More recently, a coupled tool, comprising of the thermal- hydraulics code TRACE and a multi-group 3D spatial neutron- ics solver Purdue Advanced Reactor Core Simulator (PARCS), was used to run steady state and transient simulations of the MSFR [ 11]. As PARCS only supports conventional fuel as- sembly shapes, the authors approximated the 3D axisymmetric cylindrical MSFR model with hexagonal fuel assemblies. This approach allowed for a 3D model simulation of the MSFR on a coarse mesh and showed good agreement with the COMSOL and TUDelft models with the beneﬁt of much lower computa- tion times in comparison to the OpenFOAM CFD model. On top of steady state simulations, the transient scenarios stud- ied by Pettersen include loss of heat sink, pump over-speed, over-cooling and loss of ﬂow. As Moltres is developed on the MOOSE framework, it can perform implicit multiphysics coupling on an adaptive meshing scheme over multiple processing units for accurate and e cient MSR simulations. This paper aims to present Moltres simulation results and safety analysis of the MSFR, for various fuel compositions through its operational lifetime. III METHODOLOGY This section describes the overall modelling approach and the codes used in this paper. First, the early-life and equi- librium compositions were obtained from SCALE /TRITON unit cell depletion calculations by Rykhlevskii et al. [ 12] The authors used unit cell approximations for the modeling fuel cycle performance of the MSFR and three other fast-spectrum MSRs. We selected the equilibrium composition at approxi- mately 43 years after start-up because the TRU vector change between each depletion time-step from this moment onwards is less than 3% [ 12]. For the early-life composition, it is the resulting fuel composition at 300 days after start-up. We performed neutronics calculations on SERPENT us- ing these fuel compositions and a truncated, reference MSFR core model [ 11] to generate six-group neutronics data for Moltres. The material compositions and densities for the ab- sorber, and reﬂector regions are given in Table IV. We set a neutron population of 200,000 with 50 inactive cycles and 500 active cycles. For simplicity, the heat exchanger and pump region were modeled after Hastelloy N [ 13]. Isotopes in the fuel compositions, whose neutron cross-section data are not found in the JEFF-3.1.2 nuclear data library, were omitted from these calculations The energy boundaries for the six neu- tron groups were adopted from previous MSFR studies and shown in Table V. As for the DNP, there are eight DNP groups implicitly deﬁned by the JEFF-3.1.2 nuclear data library. The reactor model is a 2D axisymmetric model as shown in Fig- ure 2. This geometry focuses on the active central region of the core. Basing the geometry on the reference model facil- itates comparisons with previous MSFR studies by various authors [ 9] [11] on the reference model for code-to-code veri- ﬁcation purposes. The 2D model is extended into a 3D modelTABLE V: Neutron energy group upper bounds used in Serpent. Group number Upper bound [MeV] 1 7 :485104 2 5 :5308103 3 2 :478 75102 4 0 :4979 5 2 :2313 6 12 TABLE VI: Temperature reactivity feedback coe cients. Composition T[pcm K1] Start-up7:390:03 Early-life7:250:03 Equilibrium6:240:03 in SERPENT by rotating it around the central axis to form concentric cylinders. Table VI shows the overall temperature reactivity feedback coe cients due to the Doppler e ect and temperature-driven fuel salt expansion. The coe cients were derived by applying a least squares best ﬁt on the multiplica- tion factors at temperatures from 900 K to 1200 K with 50 K intervals. Finally, with the group constant data from SERPENT, Moltres calculated the resulting coupled neutron spectrum and temperature distribution for a given velocity proﬁle in the core, from user-speciﬁed initial conditions towards steady state, or from steady state in the case of an accident scenario. We ﬁrther simpliﬁed the mesh geometry for Moltres from the SERPENT model by omitting the absorber, heat exchanger, and reﬂector regions. This is justiﬁed as they have negligible contribution towards the neutronics. Instead, a vacuum boundary condition was imposed on the outer surface of the blanket region, along with 923 K Dirichlet boundary conditions for temperature. The 2 cm thick blanket tank structural material was also omitted as its inclusion resulted in an extra ﬁne mesh which drastically increased computational time. Figure 3 shows the Moltres MSFR mesh. Moltres has a salt loop modeling capability to account for the ﬂow and decay of DNPs outside the active core region. Moltres applies DNP boundary conditions at the inlets and outlets of the core and outer loop, determined by the inﬂow and outﬂow values of each DNP group at each time-step for each boundary. The DNP equations also govern the DNP behavior in the outer loop to allow for ﬂow, production, and decay. The velocity proﬁle is a uniform velocity of 1:1275 ms1 throughout the active core and outer loop regions. This value is derived from the MSFR loop circulation time of 4 s. Moltres also has heat exchanger modeling capabilities, represented in 3 by the heat sink term. For this paper, it is implemented at the midpoint of the outer loop region, with the heat loss value ﬁxed to induce an approximately 100 K drop in tempera- ture as speciﬁed by the MSFR design inlet /outlet temperature di erence. Fig. 2: Cross-section of the 2D axisymmetric model used in SERPENT. Derived from the MSFR reference model [11]. Figure is not drawn to scale. All dimensions are in cm. III.A Moltres Code This subsection provides the theoretical background for the coupled neutronics /thermal-hydraulics physics imple- mented in Moltres. Moltres is an application code developed in the MOOSE framework [ 2]. MOOSE application codes solve non-linear problems through the discretization of partial di erential equa- tions (PDEs) on an adaptive coarse meshing scheme provided by LibMesh [ 14] and PetSc [ 15]. Individual terms of PDEs that deﬁne the physics involved in a system are represented in MOOSE (and its applications) by kernels. For example, the various terms in the neutron di usion equation such as the dif- fusion term, time evolution term, etc. all have a corresponding physics kernel deﬁned in Moltres. Boundary conditions are also handled in a similar fashion. Moltres can solve for an arbitrary number of neutronics groups as long as the relevant group constants are provided in a compatible text format. As mentioned in the previous Moltres study [ 1], the neu- tronics in Moltres is described by the time-independent multi- group neutron di usion equation as shown in Equation 1: 1 vg@ g @tr Dgr g+ r g g =GX g,g0s g0!g g0+ p gGX g0=1(1 )f g0 g0+ d gIX iiCi;(1) where vg=average speed of neutrons in group g g=ﬂux of neutrons in group g t=time Dg=di usion coe cient of neutrons in group g r g=macroscopic cross-section for removal of neutrons from group gs g0!g=macroscopic cross-section of scattering from g0tog p g=prompt ﬁssion spectrum neutrons in group g G=number of discrete neutron groups, g =average number of neutrons produced per ﬁssion f g=macroscopic ﬁssion cross-section for neutron in group g d g=delayed ﬁssion spectrum neutrons in group g I=number of delayed neutron precursor groups =delayed neutron fraction i=average decay constant of delayed neutron precursors in precursor group i Ci=concentration of delayed neutron precursors in precursor group i: The DNPs are governed by the following equation: @Ci @t=GX g0=1 if g0 g0iCi@ @zuCi: (2) Lastly, the governing equation for temperature in the molten salt is given as: cp@T @t+rcp !uTkrT=QsQhx; (3) where =density of molten salt cp=speciﬁc heat capacity of molten salt T=temperature of molten salt  !u=velocity of molten salt k=thermal conductivity of molten salt Qhx=heat sink; and the source term Qsis given as: Qs=GX g=1gf g g: (4) IV STEADY STATE NEUTRON FLUX AND TEMPERATURE DISTRIBUTION We ﬁrst ran Moltres in steady-state solver mode to gen- erate the six-group neutron ﬂux distribution at a temperature of 973 K with the start-up fuel composition. This is veriﬁed against the six-group ﬂux generated directly from SERPENT after collapsing the ﬁne-group ﬂux. Figure 4 shows the neutron ﬂux distributions from SERPENT and Moltres, with precursor drift disabled for code veriﬁcation. There is very good agree- ment between the six group ﬂux distributions from SERPENT and Moltres. Next, the MSFR Moltres simulations were run in transient mode, with adaptive time-stepping and an initial, uniform neu- tron ﬂux of 1014cm2s1and an initial temperature of 923 Fig. 3: Mesh of the 2D axisymmetric model used in Moltres. The grey and red regions represent the fuel and blanket salt respectively. K throughout the fuel and blanket regions, for the start-up, early-life, and equilibrium fuel compositions. The simulations reached steady state after 100 s, as the volume-integrated neu- tron ﬂux values varied by less than 0.001% over the previous ten seconds. The computational time for each case averages at 30 minutes when run on four processing units on a worksta- tion. Figure 5 shows the total radial neutron ﬂux at reactor half- height, with the three fuel compositions at steady state. The peak neutron ﬂux value is highest for the start-up composition, followed by the early-life and equilibrium compositions. All three values are slightly lower than the peak value of approxi- mately 8:61015cm2s1reported by Fiorina et al. [ 7]. This TABLE VII: Average fuel inlet temperature. Composition Inlet temperature [K] Start-up 964 :85 Early-life 925 :11 Equilibrium 916 :85 TABLE VIII: Average fuel outlet temperature. Composition Outlet temperature [K] Start-up 1067 :40 Early-life 1025 :39 Equilibrium 1016 :67 Fig. 4: Fine-group and six-group neutron ﬂux distributions from SERPENT and Moltres. Fig. 5: Total radial neutron ﬂux at reactor half-height, for start-up, early-life, and equilibrium fuel compositions at steady state. Fig. 6: Neutron group ﬂuxes at reactor half-height, for start-up fuel composition at steady state. Fig. 7: Temperature distributions within the fuel salt region, for start-up (top left), early-life (top right), and equilibrium (bottom) fuel compositions at steady state. The height and radius are in the Y and X directions respectively. may be due to the truncation of the top and bottom fuel regions of the core for the model used here. The vacuum boundary conditions were imposed closer to the center and resulted in a lower peak value. Figure 6 shows the six neutron group ﬂuxes for the start-up fuel composition. Heat transfer is dominated by advection as is observed by the temperature distributions in Figure 7. Given that ﬂux is the highest at the center of the core, most of the heat is produced there. The upward ﬂow pushes the temperature peak up to the outlet boundary. Tables VII and VIII show the average fuel inlet and outlet temperatures for the three cases. The average outlet temperatures for all three compositions are close to the MSFR outlet temperature speciﬁcation of 1023 K. Due to having a Reynolds number, signiﬁcant mixing and heat transfer is expected in the fuel salt as it ﬂows out of the core and through the outer loop. However, the center of the top reﬂector region is exposed to dangerously high peak temperatures of greater than 1100 K in all three cases as observed in Figure 7. In addition, higher temperatures are expected during dangerous accident scenarios. In actuality, the peak temperatures are likely to be lower than the values seen in Figure 7 because of turbulent mixing. In choosing a uniform velocity proﬁle, we have also implied that ﬂow is laminar, which is a weak assumption given the high Reynolds number of molten salt ﬂow in the MSFR. Turbu- lent temperature mixing could potentially be approximated by incorporating it into the conduction term as turbulent di usion. The uniform velocity approximation in our models also resulted in a signiﬁcantly di erent spatial temperature distri- bution shape as compared to the results by Fiorina et al. [ 9] and Pettersen [ 11]. Their models feature signiﬁcant ﬂow stag- nation in the fuel salt near the blanket tank. This causes their temperatures distributions to peak near the blanket tank as opposed to the center of the core as seen in Figure 7. However, improved models of the MSFR, as shown by Auﬁero et al. [10], have temperature distributions closer to the results in this paper through curved walls that optimize salt ﬂow. Total power, as deﬁned by Eq. 4, at steady state is ap- proximately 2 GW. This is 1 GW less than the rated 3 GW of the MSFR. This may be mainly due to the erroneous uniform velocity proﬁle imposed in the core. From ﬂuid dynamics, velocity along the central axis of the core should be higher than the velocity at the periphral areas near the blanket tank. The higher velocity at the center results in higher advective heat transfer away from the center of the core and cools the active region. This should result in an increase in neutron ﬂux from the strong negative temperature reactivity feedback, and thus higher heat generation. The truncation of the top and bottom fuel salt region is also a factor for the lower than expected heat production. Comparing the three di erent fuel compositions, the MSFR operates at the highest temperatures with the start-up fuel composition, followed by the early-life, and the equi- librium compositions. This may be due to the generation of233Pa as an intermediate for233U breeding as mentioned by Rykhlevskii et al. [ 12] in their fuel cycle analysis of the MSFR;233Pa is a precursor for the234U neutron poison. Fig. 8: Rise in average core fuel temperature for start-up, early-life, and equilibrium fuel compositions during ULOHS. Fig. 9: Power generated, for start-up, early-life, and equilibrium fuel compositions during ULOHS. V UNPROTECTED LOSS OF HEAT SINK ULOHS may occur due to various causes of failure in the heat exchanger or secondary coolant loop such as secondary loop pump failure or loss of secondary coolant. While it is less likely in an MSFR as it consists of 16 separate coolant loops each with its own secondary cooling system, the potential consequences must still be studied. The loss of secondary cooling is simulated by exponentially decreasing the heat loss rate in the Moltres heat exchanger kernel with a time constant of 5 s. ULOHS is expected to cause temperatures to increase, followed by a decrease in neutron ﬂux and heat generation due to the strong negative temperature reactivity feedback. Eventually, the MSFR should settle on a new steady state with a higher average temperature than the operating temperature. If the high temperature is sustained for a long period of time, it could lead to severe damage in the reactor, pipes, or other instruments. As observed in Figure 8, the average fuel temperature starts rising approximately 0.7 s after initiation. The MSFR reaches a new steady state approximately 10 K higher than the initial temperature after 30 s. While the transition time to reach the new steady state is in good agreement with results reported by Fiorina et al. [ 9], the temperature rise is much smaller than the 100 K increase reported by the same authors. Power generation falls by two orders of magnitude (Figure 9), resulting in a negligible temperature di erence between the inlets and outlets for all three fuel compositions. This could also be attributed to the relatively small total power of 2 GW discussed in the previous section, and the truncated MSFR model. In comparing the three fuel compositions, an MSFR loaded with start-up fuel composition is at a higher safety risk as it has higher temperatures under normal operating con- ditions and accident scenarios. VI CONCLUSION In this study, we investigated the steady state, and tran- sient behavior of the MSFR during a ULOHS using the cou- pled neutronics /thermal-hydraulics code Moltres. Moltres boasts relatively fast computational times due to a combina- tion of parallelization, adaptive coarse meshing scheme and adaptive time-stepping. We veriﬁed the six-group neutron ﬂux distribution at 973 K from Moltres against SERPENT. The total peak neutron ﬂux at steady state agrees very closely with published results by Fiorina et al. [9]. While the temperature distribution and total power have large discrepancies, we have accounted for the main sources of error: the erroneous uniform velocity proﬁle and the reactor model truncation. We also observed that the MSFR operates at a higher temperature with the start-up fuel composition than with the early-life and equilibrium compositions, due to the presence of233Pa. The high peak temperatures near the center of the top reﬂector are a potential safety risk. While this high temperatures could have been exacerbated by absence of turbulent di usion implementation in Moltres, it warrants the need for further safety analysis. During a ULOHS, the MSFR core fuel salt temperatures rise by approximately 10 K after a transition time of 30 s. While this time duration is similar to published results, the temperature rise is relatively low compared to the 100 K rise reported by Fiorina et al. [ 9] This discrepancy is also at- tributed to the same sources of error for the steady state results. Nonetheless, a temperature rise of this magnitude is relatively benign to reactor components, and makes a strong safety case for MSRs due to their strong negative temperature feedback, with the exception of the central top reﬂector region discussed previously. Further development of the models is necessary to im- prove the accuracy and reliability of Moltres as a coupled neutronics /thermal-hydraulics code. Further work includes using an arbitrarily deﬁned velocity proﬁle, such as a parabolic proﬁle. A better option would be the development of a Navier- Stokes module to generate an accurate velocity proﬁle. A turbulent di usion model would also be an essential addition to account for turbulent temperature mixing. The heat exchanger functionality in Moltres can be im- proved beyond the current implementation of a ﬁxed heat sink. The relevant heat transfer parameters between the primary and secondary loop in the heat exchanger would be necessary fora more complicated heat exchanger system. ACKNOWLEDGEMENTS The authors thank members of the Advanced Reactors and Fuel Cycles (ARFC) group at the University of Illinois at Urbana-Champaign for helpful discussions relating to this paper. Sun Myung Park is supported by the Singapore Nu- clear Research and Safety Initiative (SNRSI) Postgraduate Scholarship.