Development of a CFD – LES model for the dynamic analysis of the DYNASTY natural circulation loop

(cid:1) A LES model for the analysis of a natural circulation loop in presence of distributed heating is developed. (cid:1) The model takes into account the ﬂuid and the solid regions, heat generation and conduction in the pipes. (cid:1) The approach is able to reproduce stable and unstable transients of DYNASTY. (cid:1) New information gathered as stratiﬁcation and counter-current ﬂows occurring during ﬂow reversal. (cid:1) LES suitable to overcome RANS limitations in the stability and dynamic analysis of natural circulation systems. Natural circulation is exploited in nuclear systems to passively remove power in case of accident scenar- ios. In this regard, the DYNASTY experimental facility at Politecnico di Milano has been setup to increase the knowledge on single-phase, buoyancy-driven systems in the presence of distributed heating. In this paper, the development of a computational ﬂuid dynamics (CFD) model of DYNASTY is presented, focus-ing on the capability of CFD to assess the dynamic behavior of the facility. The large eddy simulation (LES) model takes into account both the ﬂuid and the solid regions, with heat generation and 3D heat conduction resolved in the pipe walls. The study, conducted using OpenFOAM, shows (i) the capability of repro- ducing stable and unstable transients of DYNASTY, (ii) new observations on the features of ﬂow reversals during unstable transients, speciﬁc circulation systems.


Introduction
Natural circulation is the result of the presence of density gradients in a fluid system, induced by temperature differences, which generate convective motion as a result of the action of buoyancy forces. In nuclear engineering, the use of natural circulation to passively extract heat from, and increase the safety of, nuclear reactors has been documented and employed in several reactor designs of Gen IIIþ, such as the AP1000 (Sutharshan et al., 2011) and the ESBWR (Rassame et al., 2017). Due to the ever-increasing concerns for improving nuclear reactor safety, passive heat extraction systems are also major players in the development of Generation-IV reactors. This is particularly true in the case of the molten salt fast reactor (MSFR), due to the peculiar characteristic of its fuel, which is fluid and also acts as the reactor coolant. This feature introduces the opportunity during accident situations to take advantage of the fluid passively flowing by natural circulation inside the system and releasing its decay heat to external heat sinks, even in absence of any external action. However, such a design inherently relies on natural circulation being able to sustain a stable, high enough mass flow-rate inside the core to extract the required thermal power without overheating the system's components or operating in an unstable oscillating regime, which could also become dangerous for the integrity of the reactor. This latter aspect is of fundamental importance, given the well-known tendency of natural circulation loops to be subjected to flow instability and mass flow rate oscillations (IAEA, 2005). In particular, for these kinds of systems, the most frequent dynamic instabilities are density wave oscillations (Belblidia and Bratianu, 1979) that arise from the dependence of the buoyancy driving forces on the pressure losses and the density and temperature distribution in the fluid (Fig. 1). Because of this https://doi.org/10.1016/j.ces.2021.116520 0009-2509/Ó 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). feedback mechanism, buoyancy-driven flows typically have an initial transient characterized by oscillations of the thermal hydraulic parameters of the system, before converging to a stable circulation state. However, for certain conditions occurring in the system, these oscillations can diverge with time (Misale, 2014), leading to limit situations such as flow reversal that can have drastic consequences on the system's reliability. Differently from the majority of the work available in the literature (Vijayan et al., 2007;Misale, 2010) an additional challenge when dealing with liquid-fuelled reactors like the molten salt reactor is caused by the heat source being distributed and the impact of this on the natural circulation stability (Krepel et al., 2014;Krepel et al., 2008).
Natural circulation was first studied, both from an experimental and modelling standpoint, using natural circulation loops (NCL) (Zvirin, 1982;Greif, 1988). These facilities are composed of pipeli-nes looped in simple geometries, through which a fluid can flow because of the buoyancy force generated by the temperaturedriven differences in the fluid's density. Temperature distributions are achieved via local or distributed external heating/cooling applied to the pipes to introduce/extract heat from the fluid. These facilities are important as they provide the possibility to observe the behaviour of natural circulation flows in a controlled environment, as well as the dependence between the system's parameters and its stability, and in providing validation data for the numerous mathematical models that have been developed to predict the behaviour of NCLs. Once insight from simplified NCLs has been obtained, this can then be applied to more complex systems, obtaining a priori information on the desired system behaviour before it is physically built. In light of this, in order to support the development of the molten salt reactor, and in particular the adoption of passive decay heat removal systems in its design, as well as assisting the development of modelling approaches able to assess the dynamical behaviour of natural circulation systems, and the effectiveness of passive cooling strategies in the presence of distributed heat sources, the DYNASTY facility (Cauzzi, 2019) was built in the DYNAMO laboratories at Politecnico di Milano (Fig. 2). DYNASTY is a natural circulation loop that can operate under distributed heating conditions and with molten salt as the working fluid.
Different modelling approaches are available to describe flow instability in natural circulation systems, according to the desired degree of fidelity required by the task at hand. Welander (1967)  A. Battistini, A. Cammi, S. Lorenzi et al. Chemical Engineering Science 237 (2021) 116520 was the first to develop an approach based on a linear stability analysis in order to draw a stability map which is function of two parameters. Starting from this pioneering work, linear stability maps (Parks, 1992) have been used to predict the asymptotic stable or unstable behaviour of natural circulation systems. However, stability maps cannot provide any information on the time-dependent behaviour of the system. To overcome this issue, 1D system codes that describe the system in terms of mass, momentum and energy balance equations have been employed. With these, it is possible to study the time-dependent evolution of any unstable transient, assess the oscillation amplitude of the mass flow rate and the temperature and also predict flow reversals. On the other hand, a mono-dimensional approach still prevents the model from accounting for local three-dimensional phenomena that can be present in buoyancy-driven transients. Consequently, the presence of different flow regimes, flow reversal and laminar-to-turbulent transition mechanisms, radial temperature distribution and other 3D phenomena can be analyzed only with more advanced approaches such as computational fluid dynamics (CFD). CFD models can be a very useful tool in the analysis of the timedependent behaviour of natural circulation systems. However, proper modelling of turbulent phenomena in these systems is crucial due to the wide range of flow conditions characterizing them and the physical complexities often involved. Most of the CFD campaigns performed on natural circulation systems rely on the use of turbulence models used in conjunction with Reynolds-averaged formulation of the Navier-Stokes equations Luzzi et al., 2017). The choice of Reynolds-averaged Navier-Stokes (RANS) models for these studies -as for most common industrial flows -is mandated by the interest in averaged quantities in terms both of the velocity field and temperature distribution and justified by the limited computational burden of RANS relative to other approaches (Pini, 2017;Cauzzi, 2019). On the other hand, this choice has direct consequences for the possibility of the model correctly reproducing turbulent phenomena during unstable transients, where both turbulent and laminar conditions, and transition between the two, are found and flow reversal may occur. Previous studies indicate that RANS models, due to the ensemble average assumption they are built on, can damp the oscillations in the system, introducing a strong bias in the evaluation of both the asymptotic stability and the dynamics of the system (Cauzzi, 2019). In addition, other assumptions often made in RANS modelling, such as the isotropy of the turbulent viscosity and the proportionality between turbulent momentum and energy transfer, as well as the near-wall scaling laws commonly employed, may fail in natural circulation and in the presence of local mixing and stratification phenomena, limiting the accuracy of the CFD model (Krepper et al., 2002;Hanjalic, 2002;Choi and Kim, 2012). In view of this, more resolved approaches to modelling turbulent flows, such as large-eddy simulation (LES), offer the potential for more detailed predictions of the three-dimensional behaviour of the flow, even if at the cost of more computational resources (Rodi et al., 2013;Blocken, 2018). In LES, only the large scales of the turbulent motion are resolved, while the smallest scales are left to be modelled with an appropriate sub-grid scale model (SGS) (Rodi et al., 2013). In natural circulation flows, this approach should be able to account for the different regimes present in the flow, contributing fully in turbulent regions and becoming negligible in the case of re-laminarizations that are typically predicted during unstable transients (Luzzi et al., 2017;Cauzzi, 2019). In this work, a 3D CFD model of the DYNASTY facility is developed and LES is applied to study the dynamics of natural circulation in the system. The model includes conjugate heat transfer with the solid walls of the flow loop and the metal temperature dynamics and, to the best of the authors' knowledge, this is the first example of an LES model with conjugate heat transfer applied to the analysis of a buoyancydriven system having the length scale of the DYNASTY facility. The main aim of the work performed is to increase current knowledge of natural circulation in DYNASTY and improve understanding of the interactions between local three-dimensional phenomena and the global dynamic stability of natural circulation loops with distributed heating. Prior to this work, DYNASTY's dynamic stability was studied through various models of increasing fidelity, from stability maps, to 1D system models (DYMOLA (Dassault Systèmes, 2019)) and CFD simulations based on RANS modelling (Cauzzi, 2019). These simulations demonstrated the possibility for the system to produce both stable and unstable transients, depending on the system's setup parameters (amount of injected heat, heating distribution, cooling temperatures). Here, the work uses two previously analysed configurations which led to a stable and an unstable transient. These configurations are simulated with the LES model in order to highlight the limitations in previous RANS modelling results (i.e., the loss of local and three-dimensional resolution due to spatial averaging) and increase confidence in the LES model's ability to predict the stability of the selected DYNASTY configurations in advance of its operation with molten salt. Because of the innovative nature of the work, the pre-processing phase of the simulations focused on defining high quality computational grids (block-structured mesh), following the criteria present in A. Battistini, A. Cammi, S. Lorenzi et al. Chemical Engineering Science 237 (2021) 116520 Battistini et al. (2020) and recapped in Section 3.3, and historically required by LES (Rodi et al., 2013). Other improvements with respect to older simulation campaigns are related to the treatment of the heat source (analysed in Section 3.2) and the boundary conditions for the outlet of the system (reviewed in Section 3.2). The LES model adopted for this analysis employs the wall adapting linear eddy-viscosity (WALE) SGS model developed by Ducros et al. (1999) for modelling the filtered part of the turbulent motion. WALE is an algebraic (0 equation) eddy-viscosity model which has shown better performance in the literature than the standard Smagorinsky model (Ma et al., 2009). This was confirmed in a validation campaign using experimental data from a conventional NCL (L2, based at Università di Genova) by Battistini (2020), where the WALE model showed good properties and the ability to perform in alternating laminar and turbulent flow regimes without the excessive damping of the turbulent structures typical of the standard Smagorinksy model (Layton, 2016), a fundamental requirement for accurately predicting NCLs' behaviour. The paper is organized as follows. In Section 2, the DYNASTY facility is described with the main geometrical and thermophysical data. In Section 3, the mathematical model is presented, in addition to the simulation setup in terms of boundary conditions, modelling of the filling tank, initial conditions, and spatial and temporal discretizations. Section 4 presents the numerical results of a stable and unstable configuration, with a focus on the analysis on the flow reversal and on the temperature distributions. Some conclusions and future perspectives are outlined in Section 5.

The experimental facility
DYNASTY is a rectangular natural circulation loop composed of AISI-316 steel pipes (with thermo-physical properties in Table 1, nominal internal diameter of 38 mm and thickness of 2 mm) and dimensions reported in the schematic of Fig. 3. The loop can operate in either natural or forced circulation conditions by switching the lower horizontal leg (GO1) to an optional section (GO2) in which the flowmeter is substituted by a pump. For the remainder of this discussion, only the natural circulation configuration will be analysed. By design, every section except the top horizontal pipe is coiled with independently operated electric resistors (GV1, GV2, GO1, GO2 in Fig. 3), which allow heating of the loop in different configurations including the distributed one that is the main focus of this paper. The independently heated sections are chosen in a way so that vertical or horizontal localized heating configurations can be employed and allow the operation of DYNASTY as a conventional NCL. However, it is possible to regulate the power output of each section to achieve uniform distributed heating conditions to emulate the conditions arising in molten salt reactors 2 . The maximum power output of the heaters is 5.3 kW and the facility is designed to operate with molten salts. For the experimental campaigns, the mixture is mainly composed of sodium and potassium nitrites and nitrates (KNO 3 ; NaNO 2 ; NaNO 3 ), the thermophysical properties of which are reported in Table 2, with these values employed in the simulations. Nonetheless, DYNASTY can also be operated with water and water-glycol mixtures, the latter of which are effective in reproducing molten salt thermophysical properties whilst reducing solidification and high temperature problems encountered when using molten salts.
Insulation is applied on the external surface of the loop's pipes, reducing thermal dispersions to a minimum, and an upper tank is installed for the loop filling and to serve as an expansion tank, whereas a lower tank is present to collect the salt at the end of the experiments. Experimental data collection is performed via thermocouples at 4 locations (shown in Fig. 3), close to the elbows of the loop, and with Coriolis mass flow meter Proline Promass F80 DN 25 in the middle of the lower horizontal leg. The cooling section, a 2.1 m portion of the upper horizontal pipe, is the only section without heaters or insulation. Instead, the outer shell is covered with copper fins to increase heat exchange and external forced convection cooling is provided by a fan placed underneath the pipe.

DYNASTY numerical model
In this work, both the fluid and the solid regions of DYNASTY are modelled (Fig. 4). The 3D computational model is taken to be as similar to the experimental facility as possible. This meant including in the model the same nominal dimensions (wall thickness, internal diameter, outer dimensions etc.) and thermo-mechanical properties reported in Section 2. The only deviations from the actual system are related to features that were believed to have limited impact on the outcome of the simulations. The discrepancies are limited to the bends' shape, which for simplicity purposes have all been modelled as circular elbows with an 0.1 m radius (instead of T junctions as in the experimental facility), the omission of the lower forced circulation branch of the facility (which is closed Table 1 AISI-316 Steel thermophysical properties (Incropera et al., 2002).
. DYNASTY schematic (lengths in mm), the thermocouples can be seen at the inlet (T1) and outlet (T2) of the cooler section and before (T3) and after (T4) the lower horizontal pipe  in the simulations and therefore does not have any effect on the results), and the positioning of the expansion tank, which is offset from its real location by 14 cm. The latter feature is a solution adopted to simplify the geometric construction of the CAD model by minimising the number of T-junctions in the structure. Indeed, as a mesh block structure was used, automatic meshing algorithms could not easily discretise penetrating geometries, therefore 4 circular bends were modelled, and only a T-junction has been modelled manually to include the tank (see Fig. 5).
The modelling of the solid walls is particularly important because it allows account to be taken of the effect of the thermal inertia of the walls on the dynamic stability of the system, as well as producing information on the operational temperatures reached in the metal. In previous simulations focused on conventional NCLs (i.e., with localized heating and cooling), the inclusion of solid regions, rather than using simplified fixed temperature or wall heat flux boundary conditions on the inside wall, was found to have a dampening effect on instability phenomena (Pini, 2017).
In the fluid domain, the filtered, weakly-compressible Navier-Stokes equations (Eqs. (1)), with gravity forcing 3 , coupled to the total energy transport equation (Eq. (3)) (Vreman et al., 1994), are solved at each time-step. The density-dependent formulation is used to correctly model the relationship between the fluid motion and the buoyancy effects caused by the density differences induced by the temperature gradients in the molten salt. A second order implicit Euler method has been used for temporal discretisation, and a second order linear-upwind scheme for both momentum and energy convection. The remaining contributions to the momentum and energy equations have been discretised with second order linear schemes. The filtering operation is required by the LES model and filters out the fluctuations at the smallest scale, which would require a very fine discretisation to be correctly resolved (Rodi et al., 2013). Instead, the impact of these on the resolved flow is modelled with a SGS model and the contribution is added to the shear stress tensor: Here, the overbar denotes a filter operation and tilde the explicit Favre-filtering operation (Erlebacher et al., 1992), defined as Þ; u, and p, which are the temperature-dependent density, velocity and pressure, respectively. Additionally, e tot is the total energy and s ij the viscous stress tensor, defined as: Finally, s SGS ij is the SGS contribution, discussed in more details in the next section. For the solid regions only a heat diffusion equation is solved. The equations are solved with a segregated approach with the OpenFOAM solver chtMultiRegionFoam (The OpenFOAM Foundation Ltd., 2018) for each region of the system, and adjacent solid-liquid regions interact via coupled boundary conditions. The solution is based on a combination of the PISO (Pressure Implicit with Splitting of Operator) algorithm (Issa, 1986), a predictor-corrector approach able to simulate transient behaviours with large time-steps, and the SIMPLE (Semi-Implicit Method for Pressure-  Linked Equations) algorithm (Spalding, 1972), devised to retrieve steady state solutions of the Navier-Stokes equations. The result is the PIMPLE algorithm (The OpenFOAM Foundation Ltd., 2018), a transient solver with improved characteristics in terms of the numerical stability of the solution, due to the iteration of the PISO algorithm -non-iterative in origin -until convergence of the solution is achieved at each time-step.

Sub-grid scale turbulence model
Because unstable transients in natural circulation loops often involve alternating laminar and turbulent regimes and complex local thermal hydraulic conditions (Luzzi et al., 2017), RANS models have been found to struggle in predicting the global behaviour of the system and the velocity and temperature distributions. Whilst RANS models allow relatively fast calculations of even complex systems by solving the Navier-Stokes equations for the mean flow only, they entirely rely on how accurately turbulent phenomena are modelled by the selected turbulence model Þ . The LES approach, instead, by resolving the largest scales of turbulent fluid motion, which are the most difficult to model due to their anisotropy, and the majority of the turbulent structures (Rodi et al., 2013), has the potential to better capture the relationship between the thermal-hydraulics of the system and its stability. The effect of the smallest scale fluctuations, i.e. the unresolved part of the turbulence spectrum, are incorporated in the resolved flow computations using an appropriate SGS model (Xu, 2003). Due to the relative simplicity of its formulation, the LES eddy-viscosity based WALE SGS model (Ducros et al., 1999) has been selected. Similarly to RANS eddy-viscosity models, the sub-grid contribution is added to the shear stress by adding a sub-grid term: where m t is the SGS turbulent viscosity. The WALE model (Ducros et al., 1999) has shown promising results for several applications (Weickert et al., 2010;Kamali-Moghadam et al., 2016), and similarly to the more well-known Smagorinsky model (Smagorinsky, 1963), calculates the turbulent kinematic viscosity from the filter width D (defined as the cube root of the volume of the cells) and the resolved strain rate tensor b S ij (Eq. (6)) and its deviatoric part b S d ij (Eq. (7) The final formulation for the turbulent viscosity is: This formulation results in an improved ability to represent the alternating turbulence regimes typical of natural circulation flows, and the contribution given by the SGS model when the flow is laminar is correctly dampened and tends to zero, differently from the standard Smagorinsky model, which has been found to overdampen the turbulent oscillations (Layton, 2016). On account of the properties of the WALE model, wall-damping functions have not been adopted and the flow is wall-resolved keeping the y þ in the first cell from the walls below unity.

Boundary and initial conditions
Due to the presence of a solid pipe wall, the boundary conditions for the fluid are limited to the interface boundary conditions between the molten salt and the solid wall and the upper outlet section. The fluid/solid interface has been modelled following the conjugated heat transfer framework, by applying continuity conditions for the temperature (Eq. (9)) and the heat flux (Eq. (10)) (proportional to the temperature gradients) between adjacent domains: For the fluid velocity, the no slip condition is applied at the solid wall, while a zero gradient condition is enforced on the pressure. The modelling of the upper tank (Fig. 4) -which serves as fluid outlet and pressure boundary -has been found to be one of the most critical aspects for modelling the DYNASTY system behaviour, as will be shown in the results section. Due to the density dependence of the compressible NS equations, the fluid expands in the facility when it is heated. Therefore, especially in the first part of a transient, the model has to allow for the exit of some of the fluid, but at the same time back-flows may be also possible should reductions in the temperature occur during the unstable transients. A set of different simulation setups (reported in Table 3) has been tested, to reveal the effect the modelling choice for this outlet has on the stability of the system. The first of the adopted conditions ("Open" case) allows the spillage of fluid, while the back-flow velocity is calculated from the flux at the outlet surface in the previous time-step. As most CFD literature suggests, pressure is arbitrarily fixed at the boundary surface, to have a reference value from which the Navier-Stokes equations can be solved. With regards to the temperature, its value at the boundary has been treated as adiabatic for out-flows and the following Robin boundary condition has been applied for back-flows: where the value of h ext is taken as 5.23 W m À2 K À1 , the result of natural convection calculations based on an external air temperature T ext ¼ 20 C. This condition allows some heat dissipation to the surrounding domain, a more moderate condition with respect to the commonly adopted inletOutlet which fixes the temperature for back-flows. A second, less-physical and more ideal condition ("Closed" configuration) has also been used. In this, no out-flows or back-flows are allowed across the outlet, and the solver compensates for any change in density by correcting the amount of mass in the system. The idea of using a closed configuration is taken in accordance with previous analyses and other examples available in the literature Cauzzi, 2019;Pini, 2017). In particular, a closed configuration is usually imposed in analytical and numerical 1D models, often employed in the analysis of the stability of natural circulation loops, and numerous works on flow instability have employed the Boussinesq approximation (Krishnani and Basu, 2016), where density changes with temperature are only accounted for in the gravitational term of the momentum equation.
To maintain the possibility for adding heat losses around the heated sections as external boundary conditions to the solid domain, the power source (given by resistance wires coiled around the pipes in the experimental facility) is applied using volumetric sources in the heat diffusion equation of the solid domain. The value of such sources (Eqs. (12) and (13)) is calculated in order to uniformly provide 5.3 kW (or 1 kW for the stable transient) of thermal power in GO1, GV1 and GV2, specified in Fig. 3: In the present work, however, it is assumed that due to the insulation applied to the solid wall the losses are negligible, and an adiabatic conditions is imposed in the heated sections between the heated solid walls and the surroundings. On the other hand, the boundary conditions for the cooled section of the wall are given as a fixed temperature, which is T cool ¼ 180 C for the stable configuration and T cool ¼ 240 C for the unstable configuration. This choice has been extracted from previous RANS simulations of the DYNASTY system for the SAMOFAR project (Safety assessment of the Molten Salt Fast Reactor), one of the Research and Innovation projects in the Horizon 2020 Euratom research programme, focused on the demonstration of the key safety features of the MSFR Cauzzi, 2019). In future works, this condition could be improved to a more realistic one, which should take into account the forced convection of air from the fan positioned below the cooler pipe. The initial conditions for the entire domain match those imposed in previous RANS simulations of DYNASTY (Pini, 2017;Cauzzi, 2019), and the start-up conditions for previous experimental analyses of similar natural circulation loops (i.e. L2, Università di Genova) (Misale, 2014;Luzzi et al., 2017). The fluid is considered initially at rest with a uniform temperature equal to the temperature of the external wall of the cooler (180 C or 240 C). Uniform pressure is also imposed, with this initial condition not showing any effect on the progression of the transients, and the pressure distribution across the facility quickly adapts to the outlet boundary value and the gravitational force (lower on top, higher at the bottom).

Computational grid and time discretisation
One of the critical aspects of LES when compared with RANS is the stricter requirement on spatial and temporal resolution. Since the majority of dynamical and three-dimensional turbulent structures are resolved in LES, the computational grid needs to be sufficiently refined for these structures to be properly captured. However, given DYNASTY's physical dimensions (12 m length, 38 mm diameter) and the great number of finite volume cells % 3 Á 10 6 cells needed to discretise the whole geometry, much care should be taken in building a computational grid with the optimum amount of elements. When using RANS, there is the necessity of reaching so called "mesh independence", i.e., the level of spatial refinement above which the results of the simulation for the same setup do not change on further refinement. In LES instead, any further mesh refinement leads to the resolution of additional turbulent structures, until direct numerical simulation (DNS) resolution is reached. Therefore, the resolution necessary to resolve the desired turbulent structures needs to be determined. In Battistini et al. (2020), a preliminary analysis of the dependence of LES accuracy on geometrical meshing in cylindrical geometries in low-Reynolds number flows was carried out. Different mesh construction approaches were tested, varying the linear step in the three cylindrical coordinates (q; h; z). The mentioned work determined guidelines on the choice of refinement required whilst constraining errors in the prediction of the pressure losses, turbulence kinetic energy and resolved shear stresses to reasonable values. These criteria have been used in the present work in the construction of the mesh. Adapting the refinement criteria found in the pipe flow analysis  and applying it in nondimensional linear steps to the present geometry, average flow parameters and the reference shear velocity typical of DYNASTY's transients have to be determined: To do so, a preliminary RANS simulation (the results of which are omitted here for brevity but can be found in Battistini (2020)) has been carried out on an preliminary mesh of DYNASTY, adopting a block-structured version of the calculation grid and refinement level following the choices found in Cauzzi (2019). To retrieve a reference value of the friction velocity, the results in El Khoury et al. (2013) can be used to relate the u bulk;avg in the system to the friction velocity u s;avg : u s ¼ 0:016359 Á u 3 bulk À 0:03529 Á u 2 bulk þ 0:071989 Á u bulk þ 4:9 Á 10 À5 ð15Þ The u bulk was estimated from the mass flow rate retrieved from the preliminary RANS simulation and equal to _ m avg % 0:2 kg s À1 : Average values of the kinematic viscosity and the density in the system were determined from the mean value of the temperature from the RANS simulation and a summary of all the average quantities used for the conversion of the geometrical mesh parameters is reported in Table 4. Using the shear velocity estimated from RANS (Battistini, 2020), and the dimensionless refinement criteria established in Battistini et al. (2020), such as the Dz þ lower than 15 and the Dr þ lower than 1 at the wall required for wall-resolved LES, the criteria used to build the DYNASTY mesh for the unstable configuration in cylindrical coordinates have been retrieved and they are reported in Table 5. Visual representations of the computational grid are shown in Figs. 6 and 7.
The computational grid has % 3:5 Á 10 6 hexahedral finite volumes and is a good compromise between accuracy and computational requirements. Despite this, around 3 months of simulation time on an average number of 70 cores (% 8 Á 10 4 CPU-hours) were required to simulate a typical transient in DYNASTY. Indeed, preliminary experimental tests have shown that DYNASTY's transient time-spans are in the order of 10 3 À 10 4 s (Pini, 2017;Cauzzi, 2019), and because of the geometric mesh requirements of the LES approach, care has also been taken for the choice of the time-step, which impacts computational resources as well as the numerical stability of the solution (Versteeg et al., 1995). To eliminate the dependence of numerical stability on the time-step, an implicit scheme has been used to discretize the Navier-Stokes equations. However, in LES models the Courant-Friedrichs-Lewy (CFL) condition (Eq. (17)), which expresses the speed at which information can propagate across the calculation domain (Courant et al., 1928), needs to be properly enforced in order to ensure the correct propagation of the turbulent structures and avoid any loss of information (Lau et al., 2012). Therefore, the time step has been controlled to achieve a CFL lower than 1 in the entire loop for the duration of the transient: This meant keeping the time-step in a range between 5 Á 10 À3 s and 10 À2 s, depending on the fluid's average velocity across the entire facility, as per (Eq. (17)).

Results and Discussion
In this section, the outcomes of two LES simulations of the DYNASTY facility are reported. The power and cooler temperature configurations for the two transients 1 kW; ð 180 C À 5:3 kW; 240 CÞ were extracted from Cammi et al. (2019). The two configurations have consistently been found stable and unstable, respectively, using different modelling approaches (stability maps, DYMOLA and CFD). The stable transient 1 kW; 180 C ð Þ shows the convergence of global parameters (such as mass flow-rates, pressure drops, temperature distributions) to steady-state values, whereas the unstable transient 5:3 kW; 240 C ð Þ is characterised by an oscillating mass flow rate that diverges with time until a periodic flow reversal with characteristic frequencies of oscillation is established.
Using predictions from previous analyses, this first configuration was tested. As reported in Fig. 8, the LES results show an initial oscillating transient that eventually converges to a steady state mass flow-rate similar to previous RANS predictions . The main differences from previous results are first the flow direction predicted in the system, which may be linked to the different propagation of perturbations that eventually lead the essentially symmetrical system (except for the outlet tank) to unidirectional flow. Moreover, the temporal behaviour of the flow rate is different, as the initial oscillations start earlier with LES for both treatments of the outlet boundary, with oscillations starting first in the Open case. Secondly, when the major part of the oscillatory transient has dampened, the fluctuations of the parameters do not reduce as much as with the RANS-modelled transient, which is to be expected as the solution of the Navier-Stokes equations is calculated from filtered instead of averaged flow variables. Nonetheless, predictions for the stability of the present configuration are in a general agreement with previous results, with the statistical steady state value of the mass flow-rate slightly increased with respect to its RANS counterpart % þ10% ð Þ . This could be Table 4 Average parameters from preliminary RANS simulation (Battistini, 2020    A. Battistini, A. Cammi, S. Lorenzi et al. Chemical Engineering Science 237 (2021) 116520 related to the evidence found in Battistini (2020), where LES predicted higher pressure losses with respect to RANS simulations. Increased pressure losses lead to increased feedback of counterforces in the loop and therefore lower flow rates and a system less prone to instability, given the same buoyancy force.
Following Cammi et al. (2019), an increase in the power injected in the system and a higher cooler external temperature should favour the unstable behaviour of the facility. However, 1D model predictions showed a completely different transient (Cauzzi, 2019) and stable operating conditions. Therefore, a further analysis with LES at these operating conditions has been deemed necessary to shed some light on the behaviour of the facility in such conditions. The same set of configurations used for the stable transient (reported in Table 3) were employed and a significant influence of the modelling of the outlet section, which is discussed in detail in the next section below, was found.

Effects of the outlet tank boundary conditions
Differently from the previous stable transient, the treatment of the outlet boundary condition was found to have a major impact on the results and the stability of the system in the P ¼ 5:3 kW; T cool ¼ 240 C configuration. Fig. 9 shows how the behaviour of the system is different depending on the boundary condition at the outlet. Specifically, the free spillage and backflow condition, coupled to the Robin condition on the fluid temperature (Open LES configuration in Table 3 and Fig. 9), which more closely simulates the entire physics of DYNASTY, including the expansion tank and its mutual interaction with the loop, tends to stabilise the mass flow rate to the value observed in the preliminary RANS simulations ( _ m ss % 0:2 kg s À1 ), made using the same settings employed for the computational grid analysis (Section 3.3) and also included in Fig. 9. In contrast the closed outlet, coupled to an adiabatic treatment of the temperature (Closed LES configuration in Table 3 and Fig. 9), which models a more ideal behaviour of the facility without the real dynamics of the fluid expansion and the related inflow/outflow, shows an oscillatory behaviour that increases with time, until periodic flow reversals are reached and maintained throughout the transient. The different behaviour observed with the two boundary treatments can be related to the role of the expansion tank acting as an accumulator, which partially absorbs the alternated hot and cold plugs that are present in unstable conditions (described in more details in the following sections), dampening flow oscillations in space and time. In Fig. 9, the Open LES transient shows some oscillations at the beginning of the transient, but the presence of the tank is sufficient to damp these oscillations and drive the flow rate to an almost stable steadystate value. Instead, without the tank (Closed LES), these initial oscillations are sufficient to trigger the instability and eventually reaching flow reversal. Clearly, if the tank is included, the correct behaviour of the facility can only be predicted if its physical behaviour, and effect on the system stability, are properly modelled. Overall, therefore, the LES model is able to detect and handle an unstable system operation, but the accurate modelling of the outlet section and the upper tank, as well as its impact on the correct prediction of the stability boundary, will need to be carefully addressed in future studies at the end of the experimental campaign, when experimental data will be available. It is also important to note how the RANS simulation predicts a stable system even without the presence of the expansion tank (Closed RANS in Fig. 9). Therefore, the additional turbulent viscosity introduced by the RANS model seems to have a stabilizing effect on the system behaviour. This is shown also near the end of the transient, where the flow rate reaches a stable value and does not show even the small oscillations found in the Open LES system.

Flow reversal analysis
It is interesting to look in more detail at an integrated-value quantity such as the mass flow rate _ m during the unstable transient. This shows continuous oscillations with alternating positive and negative peaks and zero-flow conditions between the peaks. This behaviour has already been observed in previous works Luzzi et al., 2017), and it was characterised by the alternate presence of different flow regimes, from fully turbulent to laminar. However, examining the LES predictions during   and the result of the LES simulations (uniform heating at 1kW and cooler temperature at 180°C). A. Battistini, A. Cammi, S. Lorenzi et al. Chemical Engineering Science 237 (2021) 116520 these reversals, it is observed how a purely laminar flow is never established during the unstable transient. Instead, countercurrent, superimposed flows in the horizontal sections of the facility and the coexistence of oppositely directed flows in the vertical sections are observed with LES, as shown in Fig. 10. This effect cannot be captured in 1D models, with the resolution of the local flow conditions being beyond the resolution capabilities of such approaches. At the same time, in RANS CFD details will be lost due to the averaging procedure, and the accurate modelling of complex phenomena such as counter-current buoyancy-driven effects in the transition and turbulent regions is particularly challenging. Instead, these same phenomena are mostly resolved in LES, which therefore proves to be extremely valuable for their detailed analysis and understanding. Furthermore, differently from conventional natural circulation systems, heating is distributed throughout the entire facility in DYNASTY (except for the cooling section). Therefore, in the vertical section shown in Fig. 10, the fluid particles closer to the vertical wall during flow reversal tend to be hotter, and they are pushed upwards by buoyancy forces, in opposition to the downward flow in the colder central region of the pipe. Consequently, during the inversion phase, the motion of the fluid is mainly driven by the local temperature distribution, because of the lower inertia of the bulk motion. The countercurrent circulation of a central descending region of colder fluid and a cylindrical boundary of ascending hotter fluid generates turbulence, proven by the chaotic directions of the fluid particles at the boundary between ascending and descending flow in Fig. 10. This analysis shows how, during flow reversal, the stratification and the coexistence of counter-current flows occurs, more notably in the horizontal sections of the facility but with some repercussions in the vertical legs. These stratified flows, if integrated, result in low mass flow rates, although the local flow is far from zero. An analysis based on stability maps, built on average dimensionless parameters such as the Reynolds (Re) number or the Grashof (Gr) number, or on 1D system codes such as RELAP5 3D (Idaho National Laboratory, 2015) or DYMOLA (Dassault Systèmes, 2019), both using as input the low integral mass flow rate, would not detect this phenomenon. Therefore, global stability may be empirically predicted, but an important physical aspect of the facil-ity's behaviour, as well as its impact on the stability map, is entirely neglected.

Effect of instability on the temperature distribution
Dynamic instabilities are dangerous because, due to the mass flow oscillations, hot and cold plugs of fluid circulate around the loop, as already mentioned and described in Welander (1967). These may result in large oscillations of the operational quantities, which may exceed their design limits, with consequent over-stress and eventually damage to the facility. This effect is clearly noticeable in the temperature distribution within the facility during the transient. The adiabatic mixing temperature within the loop in the stable transient, plotted in Fig. 11, shows a stationary distribution once a steady-state condition is reached (hence the limited discrepancy between the two times). On the other hand, in the case of an unstable transient (Fig. 12), the adiabatic mixing temperature is fluctuating both spatially and temporally, proving the presence of locally hot and cold plugs of fluid travelling along the facility. These plugs are local regions of the fluid in which the density is consistently different from the remaining flow. Therefore, the plugs cause buoyancy counter-forces against the bulk circulation of the fluid. Eventually, the joint effect of pressure losses and these counter-forces may cause flow reversal. In addition, the hot plugs are not sufficiently cooled down by exchanging heat in the cooling section, because of their limited residence time, and, especially in uniformly heated fluids, they tend to increase their temperature while travelling in the remaining sections of the loops, more than in the corresponding stable case. Most importantly, the effect of the dynamic instability can also be observed when sampling the temperature of the solid wall. As shown in Fig. 13, the oscillating behaviour of the mass flow-rate compromises the heat extraction, inducing higher temperatures and therefore increasing the probability of damage to the structural materials.

Buoyancy effect on the temperature distribution
It is also interesting to focus on the behaviour of the fluid's temperature distribution on the pipe cross-section in the cooler por- Fig. 10. Flow reversal behaviour. Fig. 11. Adiabatic mixing temperature along the full length of the facility in the Open LES -P ¼ 5:3 kW; T cool ¼ 240 C configuration. (Start point: at the right of the cooler and then proceeding clockwise, through the heated sections first).
A. Battistini, A. Cammi, S. Lorenzi et al. Chemical Engineering Science 237 (2021) 116520 tion of the loop, which shows an increasing asymmetry along the cooler's length (see Figs. 14 and 15). The asymmetry presents itself as a stratification of layers of fluid with different temperatures, and this was also noticed in previous works on DYNASTY  for stable transients at low power (1-2 kW). Here, the stable results obtained in the Open configuration also allow examination of the same effect at the maximum power input (5.3 kW). As shown in Fig. 15, a layer of cold fluid builds-up in the lower part of the pipe. The cold particles have lower inertia, and at the same time are heavier than the adjacent hot particles. The concurrence of these two characteristics generates the observed distribution, which increases in asymmetry along the cooling section. This can also be observed by looking at the contours of velocity entering and exiting the cooler section (Fig. 16), and streamlines of the fluid      A. Battistini, A. Cammi, S. Lorenzi et al. Chemical Engineering Science 237 (2021) 116520 particles along it (Fig. 17), where those in the layer near the wall are cooled and recirculate down along the walls and accumulate at the bottom of the pipe. Due to the longer trajectory travelled by the colder particles, their residence time in contact with the cold wall is prolonged, leading to lower velocities. This phenomenon needs to be addressed in future works, as it may add uncertainty to the evaluation of the cooling capabilities of the system.

Conclusions
The present paper focused on the development of a singlephase CFD model, including the modelling of the solid walls, for the study of natural circulation dynamics with a distributed heat source. The analysis focused on the DYNASTY facility, a natural circulation loop built at the Energy Laboratories of Politecnico di Milano. In this model, with the aim of increasing the ability to predict the physical behaviour of the system and overcoming some of the limitations of RANS models applied in previous works, LES is employed. Due to the high computational cost of LES, this study focused on a small number of configurations previously associated with stable and unstable conditions . This was undertaken to confirm previous results and at the same time highlight additional physical insight into the facility's behaviour an LES is able to provide. The results have confirmed previous predictions for a stable transient in a low power configuration (1 kW, cooler temperature = 180 C). For a purportedly unstable configuration (5.3 kW, cooler temperature = 240 C), an unstable transient was obtained using a closed boundary condition that does not allow any fluid to leave the facility, an ideal condition usually employed in 1D calculations. However, a sensitivity analysis has underlined a strong dependence on this outlet boundary condition. When a more physical open outlet boundary condition was introduced (i.e., the flow is allowed to exit and re-enter the loop following changes in the density), together with the modelling of the salt filling tank, a stable system was observed. Additional studies are therefore necessary on how to properly model this aspect of the facility and its impact, and the limitations introduced by the closed approximation, on the correct prediction of the stability boundary, and these will be carried out once the experimental studies have been completed. Important local phenomena were observed. During unstable transients, flow reversals were shown to occur with stratified counter-current flows in the horizontal cooled and heated pipes, which, when integrated, result in very low mass flow-rates. In the vertical sections, instead, a radial stratification occurs with a distinct annular upflow region around the heated wall, in opposition to the core downflow region in the centre of the pipe. This observation improves previous understanding, mostly related to re-laminarisation occurring during the flow reversal, driven by the very low measured flow-rates obtained when counter-current flows are integrated. Analysis of the temperature distributions for both stable and unstable transients confirmed the presence of hot and cold fluid plugs travelling along the facility during unstable transients, contributing to destabilisation of the flow. The alternation of clockwise and counterclockwise circulation typical of unstable transients causes the deterioration of the heat removal capabilities of the system and a distinct rise (10-15 C) in the temperature of the solid walls, detected by including conjugate heat transfer in the CFD model. Results also showed an increasingly asymmetrical temperature distribution in the horizontal cooling section, leading to thermal stratification in the pipe cross-section. This aspect of the flow and its impact on the system's heat removal capabilities will need to be further investigated, and the availability of a reliable numerical tool which includes the modelling of the solid walls, will def-initely help in supporting and corroborating experimental evidence. In addition to the modelling of the filling tank and the related boundary condition, which should be analysed in depth including more appropriate designs such as a free surface on top of the system with a movable mesh interface (ALE), further development of the model should also focus on increasing the fidelity of the boundary condition used to model the cooling fan. Given the computational costs of simulations, mesh requirements were derived in a simplified geometry, in flow conditions as close as possible to the DYNASTY facility. Therefore, additional sensitivity on the mesh resolution will also be necessary, including comparisons with more refined simulations and a specific focus on the thermal field resolution, considering the high Prandtl number of the molten salt employed. In this regard, the requirement of an ad hoc model for the sub-grid heat fluxes, in conditions where some commonly adopted assumptions such as the gradient diffusion hypothesis or the Reynolds analogy are expected to fail, will be worth investigating. Finally, the validation of the model, and its quantitative predictions of the stability boundary, will be a priority once experimental data from DYNASTY become available. Once validated, the model, by providing detailed insight on the natural circulation behaviour inside the loop, will support the development of passive heat removal systems for the molten salt reactor.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.