© A. Nalbandyan et al., Published by EDP Sciences, 2022

1 Introduction

Passive systems operating without active driving components, such as pumps, are of interest for engineering applications where system autonomy and reliability have to be ensured. This is particularly relevant in the nuclear sector, where power plant safety functions (reactivity control, core and containment cooling, prevention of radioactive release) should be ensured in every different situations, ranging from operational conditions to accidental scenarios [1]. For example, the cooling function, and related passive decay heat removal systems, play an important role in increasing the reliability of a nuclear reactor, as highlighted by the Fukushima accident.

Some advanced Generation III and III+ water reactors such as the AP1000 and the ESBWR as well as the Generation IV reactors foresee passive systems for core decay heat removal after reactor shutdown [2–4]. These kind of passive systems may rely on natural circulation features in order to ensure the cooling capability. Being the absence of a pump the main advantage - in terms of passive safety but also in terms of simplicity and cost reduction -, as main drawback, the natural circulation systems are prone to unstable behavior due to the strong coupling between the buoyancy forces (that depends on density differences) and friction forces [5]. Whereas in the maeopaority of the cases the heat source and the heat sink are localized, for some specific applications, distributed heat source might be present in the loop, e.g. due to the internal heat generation by the working fluid. This situation is particularly relevant to the Generation IV Molten Salt Reactors (MSRs), where the liquid fuel acts also as coolant causing the fission heat generation to be internally delivered in the fluid [6]. The unstable behaviour is characterized by large oscillations of the main thermal-hydraulic parameters such as the mass flow rate and temperature. This is not a desirable situation for any engineering application in particular if the system can reach conditions that may jeopardize the safe operation of the structural components. Thus, it is important to be able to correctly describe the flow regime of natural circulation systems. In this context, a useful equipment is represented by a Natural Circulation Loop (NCL), i.e., a rectangular or a toroidal loop with a heat sink and a heat source where the fluid flow is driven by natural circulation as a function of friction and buoyancy forces [5,7,8]. The simplicity of these systems allows focusing on the physical and phenomenological basis of the buoyancy-driven circulation and represents an ideal validation benchmark for testing modelling capabilities. Among the NCLs, the DYNASTY facility [9] located at Politecnico di Milano (Fig. 1) is aimed at investigating the stability and the dynamics of natural circulation in presence of distributed heating, in particular to investigate the decay heat removal mechanisms for the Molten Salt Fast Reactor [6]. Previous studies on the stability of singlephase NCLs with both localized and distributed heat sources applying analytical and numerical methods [7,10], indicate, that natural circulation with internal heat generation is potentially more prone to the instability with respect to the case with localized heat sources. Thus it is of interest and practical importance to further investigate the stability of NCLs with distributed heating. In order to gain in-depth information on the flow behavior inside the facility, in this paper a CFD model of DYNASTY - based on the OpenFOAM CFD framework - is presented with a comparison of several Reynolds-Averaged Navier Stokes (RANS) turbulent models. The CFD approach allows investigating 3D phenomena as well as the radial profile of the temperature. As a major outcome, the CFD simulation results are validated against experimental results carried out using water as working fluid. Despite the thermophysical differences between water and molten salts, it has been decided to start the experimental campaign with water to limit the requirements on temperatures. In addition, these experimental data represent the first data obtained with DYNASTY in distributed heating mode and the knowledge acquired with water will be employed for future experiments with molten salts. As a further verification, the modeling results are also compared to previously developed and tested methods, e.g. the stability maps and the 1D Object Oriented Modeling approach [7]. As such, this work presents a unique study of an NCL with distributed heat source comparing not only analytical and 1D modelling methods but also a full-scale detailed CFD model against an experimental case.

The paper is organized as follows: in Section 2, the description of DYNASTY facility is provided, followed by a detailed description of the CFD model in Section 3, such as the model geometry, meshing approach, the turbulence modeling, modeling of pressure and heat losses as well as the initial and boundary conditions together with the numerical schemes used by the solver. Section 4 refers to the analytical and Object Oriented modeling approaches and finally Section 5 presents and discusses the simulation results compared to the experiment, followed by a conclusion given in Section 6.

Fig. 1 DYNASTY facility.

2 Description of DYNASTY facility

DYNASTY is a large NCL facility designed to operate with molten salt as the thermal carrier but with the flexibility to run with multiple thermal carriers as water or glycol. All the components are made of stainless steel (AISI 316) to withstand the operational temperatures. The main geometrical dimensions of the facility and operational parameters for water are summarized in Tables 1 and 2, respectively.

A schematic view of DYNASTY is presented in Figure 2. The facility is divided into five sections, namely the cooler, the downcomer (pipe 1), the riser (pipe 3), the horizontal leg (pipe 2), and the pump leg. The two vertical legs are labelled riser and downcomer, assuming (just for naming purposes) a clockwise flow of fluid inside the loop. The bottom part of DYNASTY presents two parallel sections (horizontal leg and pump leg) to be used for natural circulation (NC, as in this work) and forced circulation experiments respectively. Each of the two sections can be isolated from the rest of the loop through valves placed at both ends of the section. The DYNASTY cooler is the top horizontal section, which cannot be heated and it is a finned pipe coupled to a fan which generates air, in cross-flow with the tube, at ambient temperature.

Heat is provided by electrical strips installed on the exterior of the pipes to mimic the internal heat generation. On the outside of the pipes and the heating system, a Rockwool thermal insulator is applied to reduce thermal dispersion to the environment. The heating elements are grouped in four sections GV1, GV2, GO1 and GO2 that are responsible for heating the riser and top part adjacent to it, the downcomer and the top part adjacent to it, the horizontal leg and the pump leg (Fig. 3). Each of these sections can be powered independently from one another, and the provided power can be distributed between the sections. The cooler section cannot be heated and thus does not have any heating stripes. Thanks to the independent powering of each leg of the system, it is possible to operate the facility with a wide range of power distributions ranging from localized configuration (e.g., Horizontal Heating Horizontal Cooling or Vertical Heating Horizontal Cooling) to distributed heating condition powering all three legs.

As for the instrumentation, DYNASTY is equipped with temperature sensors for the fluid and for the pipes, and with a Coriolis-effect mass flow meter on the bottom horizontal leg. The fluid temperature is measured with four resistance temperature detectors (TC1, TC2, TC3 and TC4) that are placed in positions relevant to the heat exchange section, i.e., cooler inlet, cooler outlet, down-comer outlet and riser inlet (Fig. 3). These locations are relevant since they allow computing the fluid temperature difference, which is the driving force of NC for the different heating configurations. In case of distributed heating, the driving fluid temperature difference would be TC2-TC1 as the thermal gap across the cooler.

Fig. 2 Schematic view of DYNASTY.

Fig. 3 Heating and instrumentation in DYNASTY facility.

3 The CFD model

As pointed out in [11], the stability of natural circulation is influenced by the wall thermal inertia and the heat exchange between the fluid and the solid wall. The modelling of both the fluid and the solid domain (represented by the pipe) is then of paramount importance to correctly reproduce the nature of the equilibrium state of the system. To this aim, a Conjugate Heat Transfer (CHT) approach is selected for the CFD model. In particular, the chtMultiRegionFoam heat transfer solver [12], available in the OpenFOAM v1806 version, is used in this work.

Note, that the CFD models steel pipes but for the sake of simplicity the heating stripes and the insulator are only mimicked via appropriate boundary conditions. For modeling the all-external heat flux in the heated regions of the loop, the so called external Wall HeatFlux Temperature [12] boundary condition is used on the external surfaces of heated pipe sections. This boundary condition applied on the pipe external wall surface can be imposed in terms of power, heat flux or heat transfer coefficient.

The power mode is used for the heated sections, whereas the cooler outer surface is set to coefficient mode with the overall heat transfer coefficient and the ambient temperature as input parameters. The thermal heat losses to the environment are treated by reducing the net power provided to the heated sections of the loop.

The chtMultiRegionFoam solver considers polynomial dependence on the temperature for the thermo-physical properties. The main thermo-physical properties for the water and the piping steel used in this work are presented in Tables 3 and 4. As for the pressure losses in the loop, in addition to the frictional pressure losses calculated by the CFD model, the localized pressure losses are present in the loop due to the water filling tank and the mass flow rate meter. In the OpenFOAM model the localized pressure losses are modeled as a porous zone and the Darcy-Forchheimer model is applied [21] (1)

where D and F are the friction coefficients; D is for the viscous losses (also known as the Darcy component) and F is for inertial losses (also known as the Forchheimer component).

For e.g. pressure loss in x direction: (2)

For the water filling tank the pressure loss is calculated using empirical correlations available in hydraulic handbooks as a function of the flow characteristics, such as the Re number as well as the roughness of the pipes. For the mass flow rate meter, the mass flow rate dependent pressure drop correlation is provided by the vendor.

The D and F coefficients are subsequently calculated as follows: (3) (4)

where D and F are the viscous loss and the inertial loss coefficients. Based on the pressure loss calculation for the mass flow rate meter and for the water filling tank, the D and F coefficients in the flow direction are calculated as 1.7 × 10⁵ (m⁻²) and 3.6 × 10³(m⁻¹) respectively.

3.1 Model geometry and Mesh

The geometrical model of DYNASTY (depicted in Fig. 4) features several simplifications compared to the actual geometry presented in Figure 2. The pump leg is not modeled because in the scope of this work only natural circulation regime is considered. The mass flow rate meter is included by means of porous medium instead of modeling the actual complex structure of the device as the only impact on the flow is the pressure loss caused by the mass flow rate meter and this can be modeled using the vendor provided information. The draining tank is not modeled as well. The finned structure of the cooler is also not included in the CFD model for the sake of simplicity. The impact of the fins is taken into account in the calculations of the heat transfer coefficient between the pipe and the ambient. The meshing is facilitated using ANSYS 19.2 Workbench [14]. The cooler section which includes the water filling tank, is sliced, so that the largest part of it can be meshed using hexahedral structured grid, meanwhile the tank itself is meshed using an unstructured tetrahedral mesh (Fig. 5).

Meshes of different coarseness are generated to test the model for mesh independence. The mesh average orthogonal quality and the maximum skewness are reported in Table 5 and are complying with standard meshing guidelines for ANSYS Workbench suggesting an average orthogonal quality above 0.2 and a maximal skewness less than 0.95. These two characteristics are very important to control the mesh quality, as they show how close the mesh elements are to the optimal size and shape.

Another important meshing criterion is the dimension-less wall distance otherwise known as the wall y+ value.

It is defined as (5)

where Δy is the first cell height of the mesh boundary layer, ρ is the density, µ is the dynamic viscosity and U_T is the frictional velocity which can be calculated as: (6)

where U is the free stream velocity and C_f is the skin friction factor; an empirical coefficient determined for internal flows as [14]: (7)

The dimensionless wall distance helps describing the near wall flow, which is subject to numerical and modelling challenges due to the viscosity induced effects. Usually, the near wall flow can either be modelled by resolving the viscous sub-layer, or by adopting wall functions to approximate the flow behavior across it. If the viscous sub-layer is to be resolved, a very fine mesh is usually required, with y⁺ ≈1. This is frequently not viable for a large industrial model and thus wall functions are generally adopted for near wall flow modelling. For the wall functions to yield valid results, the wall y⁺ value lower limit should be chosen correctly, because too low values may lead to incorrect modelling results. In this work we kept the wall y⁺ value between 30 and 200. With this y⁺ value, the OpenFOAM k, ω and ε wall functions are employed.

Fig. 4 CFD model of DYNASTY.

Fig. 5 DYNASTY mesh: closeup at inflation layers applied on water body.

3.2 Turbulence models

The shear stress in Navier-Stokes equations is usually represented as a sum of viscous and turbulent components. The turbulent component, also known as Reynolds stress, is usually modelled rather than solved directly for, in order to avoid large computational burden. A common way is to adopt the Reynolds Averaged Numerical Simulation (RANS) approach, where the Reynolds stress variable is represented as a sum of the mean steady value and the unsteady fluctuations. These unsteady fluctuations are represented by the Reynolds stress tensor and require additional closure models to relate the stress tensor to the mean values of the flow variables. These closure models are known as turbulence models and they commonly rely on the introduction of a parameter that relates the Reynolds stress to the mean flow, called the eddy viscosity. Similar to the molecular viscosity, the eddy viscosity is responsible for internal momentum transfer via eddies that are formed in turbulent flow. The eddy-viscosity turbulence models are classified as zero up to five equation models depending on the number of transport equations they solve. In this work we use the RANS approach to model shear stress and we test three RANS turbulence models, namely: k − ω SST, k − ω SSTLM and realizable k − ε models.

The k − ω Shear Stress Transport (k − ω SST) turbulence model is a two-equation model that is solved for the turbulent kinetic energy k and the specific turbulence dissipation rate ω. The implementation adopted in Open-FOAM is based on the solution of the transport equations for k and ω with adoption of a blending function which depends on the distance to the wall [15]. When moving away from the walls, this blending function has a value of 0, which corresponds to applying the k − ε model, while next to the wall the value of the blending function is 1 and subsequently the k − ω turbulence model is used. A second blending function is used as well, which prevents buildup of turbulence in stagnation zones. The k − ω SST method thus combines the advantages of normal k − ω model in the near wall region with the accurate performance of k − ε model in the free stream region. This model is used in many industrial flow simulations especially when flow separation or adverse pressure gradients are present. The second turbulence model tested in this work is the Langtry-Menter four equation k − ω SSTLM model [13].

The model introduces a transition Re number as: (8)

where the Tu factor depends on the free stream velocity and the turbulent kinetic energy and is defined as: (9)

The model performs well at the low-Re numbers and in predicting flow transition [16]. Based on the experimental data and the main geometrical parameters of the loop, the transition Re number is calculated to be 580, which points to the fact that in the case of natural circulation with distributed heating the transition Re number can be much lower than e.g. for an infinitely long straight pipe, for which the transition Re is around 2300 (for water).

Finally, the third model tested in this work is the realizable k − ε turbulence model. In contrast to the standard k − ε model this model adopts a new, more exact transport equation for calculating the turbulence energy dissipation rate and assumes the turbulent viscosity being a function of the mean flow parameters. The realizable k − ε model is more accurate and reliable for many application than the standard k − ε model. It is especially successful in describing complex flows with rotation, vortexes and stagnation zones.

3.3 Case initialization

A transient experimental test case is modelled in this work, with a CFD time step of 0.1s wherein the total power provided to the loop is 450 W, distributed 2/3 to the hot leg and 1/3 to the cold leg. The heat losses are taken into account by decreasing the net power provided to the heated legs: the losses have been calibrated through an energy balance performed on the facility at the end of a heating transient, in order to assess the net power provided to the fluid. The cooler is modeled by providing the convective heat transfer coefficient calculated based on semi-empirical correlations [17]. The case is initialized from a thermal equilibrium between the pipes and the water, the starting temperature of both being 333 K. A transient with a total duration of 5000 s is modeled, being the time required to reach the steady-state equilibrium from experimental evidences. The goal of the simulations is to determine the stability of the equilibrium point (at given power and cooler heat transfer conditions), and to validate the capability of the CFD analysis to reproduce the experimental results.

The solver accounts for the conjugate heat transfer between solid and fluid regions as well as the buoyancy and the turbulence effects adopting a segregated solution strategy; first the equations for the fluid region are solved followed by the solution for a solid region. For the fluid region, the PIMPLE algorithm [18] with three correctors is adopted for pressure, velocity and energy equations, meaning that the pressure is corrected three times within the PIMPLE loop. The outer correctors are two, meaning that within a time step the PIMPLE loop is performed two times for the whole set of the equations before moving to the next time step. As for the numerical schemes, first order implicit Euler scheme is used for the time derivative terms, second order unbounded Gauss linear and first order bounded Gauss upwind schemes are used for the divergence terms.

The boundary conditions are as follows:

• Contact regions between solid and fluid - A special temperature coupling boundary condition available in the OpenFOAM, TemperatureCoupledBaf-fleMixed [12], is used for the temperature on all contact zones between the fluid and the solid. This boundary condition represents the continuity of the temperature across the fluid-wall interface.

• On the tank outlet an inletOutlet boundary condition is used. This boundary condition usually behaves as a zero gradient Neumann boundary condition, except when there is backflow into the domain; then the inletOutlet boundary condition changes to a fixed value to prevent a non-physical flow re-entry situation. The boundary condition for the velocity in the liquid zone is set to uniform zero fixed value on all boundaries and to pressureIn-letOutletVelocity on the tank outlet boundary; this boundary condition applies a zero gradient condition on the outflow and for the inflow a velocity derived from an internal cell values is applied.

• Cooler external surface - externalWallHeatFluxTemperature boundary condition in coefficient mode is used: the overall heat transfer coefficient and the ambient temperature are specified.

• External surfaces of the heated sections of the loop-externalWallHeatFluxTemperature boundary condition in power mode is used.

Both the fluid and the solid zones temperature are initially set to 333. 15 K, and the fluid velocity is set to 0 m/s.

4 Analytical and object oriented modeling approaches

An additional verification of the CFD model is performed by comparing the results to the simulations performed with 1-D object oriented DYNASTY model and semi analytical calculations based on stability maps. In this section the two methods are described briefly.

4.1 Stability maps

Stability maps are a simple but powerful tool aimed to provide information about the stability of natural circulation over a large range of conditions, useful in the design phase of circulation loops. They are graphs drawn in the space spanned by two parameters (e.g., Reynolds and Prandtl number) where a neutral stability curve separates the region of asymptotically stable equilibrium points of the system from the unstable ones (Fig. 6). They rely on the modal linear analysis approach and on a semi-analytical treatment based on some simplifying assumptions (incompressible and mono-directional fluid, constant flow regime in any point of the loop, Boussinesq approximation). As opposed to CFD simulations, the stability maps require pressure drop and heat transfer correlations for the analysis. Starting from the conservation equations for mass, momentum and energy for the fluid and energy equation for the solid wall, a steady state can be found (in terms of mass flow rate and temperatures) which corresponds to the equilibrium point investigated by the stability map. The stability analysis involves a linearization of the time-dependent version of the aforementioned equation with the solution of every state variable represented as a steady-state solution plus a time dependent perturbation. The perturbation of a generic state variable can be written in the following form:

(10)

ω being the complex pulsation of the perturbation. If equation (10) is substituted into the governing equation, the system can be solved for ω. If the real part of ω is positive, the amplitude of a perturbation grows exponentially in time, hence the internal dynamics of the NCL amplifies perturbations and the system is unstable. On the contrary, if the real part of ω is negative, the exponential is decreasing in time, the perturbation is dampened, and the system is stable. ω = 0 represent the neutrally stable curve. Detailed description of stability maps derivation process can be found in e.g. [7,10].

Fig. 6 DYNASTY stability map example: distributed heating configuration, water as working fluid.

4.2 1D object oriented modeling

One drawback of the stability map approach is the limitation to an asymptotic analysis. To overcome this problem, and to be able to track the time evolution of the main variable of interest, a 1D modelling has been developed in the past, based on the Modelica object oriented language using the Dymola environment [19]. This approach consists in the solution of the one-dimensional, time dependent, non-linear governing equation (mass, momentum and energy balance for the fluid, and energy balance for the solid wall). Similar to the stability maps, also the one dimensional approach relies on pressure drop and heat transfer correlations. Modelica is an object-oriented, acausal, equation-based language used to simulate physical systems. It relies on a wide range of validated libraries containing the model of different multi-engineering components (thermal hydraulics, electrical, mechanical, …). The 1D object-oriented DYNASTY model is presented in Figure 7.

The thermal hydraulics components of the model are taken from the ThermoPower library [20] along with other components specifically developed for distributed heating and collected in the ThermoPowerIHG library. The development and validation of the model is described in [7]. The main component in the model is represented by a pipe, in which the mass, momentum and energy balance equations for the fluid are implemented, encapsulated in a metal tube and subject to a volumetric heat source. Additional components are added in order to represent the cooler and the fan, and the localized pressure drop caused by the mass flow meter, the elbow and the junctions. Pressure drops and heat transfer are modeled through dedicated semi-empirical correlations, e.g. the Darcy friction factor for the transition zone is modeled by employing the correlation derived in [10], whereas the convective heat transfer coefficient is determined using the correlation for Nu number according to Churchill-Bernstein [17].

Fig. 7 1-D model of the facility.

5 Results and discussion

The main results of the numerical simulation of the DYNASTY loop are provided and discussed in this section. The flow evolution as a function of time is investigated with a main focus on analyzing how the different turbulence models predict the flow stabilization i.e., when the mass flow rate and the temperature difference across the cooler stabilize around some mean value. This mean value is then compared to the experimental results in order to evaluate the behavior of various turbulence models and to establish which one predicts the steady-state equilibrium closest to the experimental results.

Mesh independence is checked by comparing the average mass flow rate in the loop for different meshes and is reported in Table 6. The fine mesh with 6×10⁶ elements is chosen for the simulations.

The CFD simulations are run on a HP Proliant SL230 Gen. 8 cluster with Intel Ivy-Bridge Xeon e-5-2880v2 2.GHz processors. The wall time is reported in Table 7 per each turbulent model, for the same mesh.

In Figure 8 the mass flow rate of the water is shown as a function of time, with the last 1000 s when the flow is stabilized, depicted in Figure 9 for the purpose of more clear comparison. The CFD results are compared to the experimental data and Modelica results. Both CFD and Modelica simulations reflect the initial oscillation of the mass flow rate due to the typical initial transient of the natural circulation where hot and cold fluid plug start circulating in the circuit. However, the detailed comparison of the initial transient oscillation regime of the flow is out of the scope of this paper. It is worth mentioning though, that the initial oscillations are affected by the turbulence parameters of the chosen turbulence models, as well as by the overall startup and initial conditions of the system, including but not limited to the power, ambient temperature and fluid initial temperature. In this light, a more appropriate figure of merit for the validation of the CFD model are the steady-state values. The mass flow rate stabilizes starting from t = 2000 s, with the realizable k − ε turbulence model further exhibiting some small oscillations and finally stabilizing after 4000 s with the average stable mass flow rate being 0.034 kg/s. This is ≈16 % higher than the experimental average of 0.029 kg/s. The k − ω models are in better agreement with the experimental results, overestimating the mass flow rate less than 10%. The relative error calculation is based on the results from 2000 s to 5000 s interval, where the mass flow rate does not exhibit large oscillations (see Tab. 8 for a result overview).

In Figure 10 the temperature difference between the cooler inlet and outlet sections is depicted, with a close-up view on the last 1000 s provided in Figure 11. The k − ω SST model predicts for the stabilized flow a temperature difference between the cooler inlet and the outlet of 4.5 K, which is on 1.8 K higher than the temperature predicted by the experiment. Similarly, the k − ω SSTLM turbulence model has a difference of 2.3 K with the experiment, while the realizable k − ε has a slightly lower difference of 1.6 K (see Tab. 9 for a result overview). These results are in agreement with the mass flow rate prediction by each turbulence model for the stabilized flow. The discrepancy is attributed to the technique adopted for the modelling of the distributed heat as a constant power applied on the heated pipes which is most likely falling short to capture the heat losses during the transient. A more precise heat loss model implementation can help improve of the current results, however for the first assessment of the CFD model an overall agreement in the temporal behavior of the heat transfer across the cooler is achieved.

Modelica simulations on the other hand, deliver the closest results to the experiment based on the detailed heat loss modeling, which takes into account the convection and irradiation that occur on the outer shell of the pipe.

In Figures 12 and 13, the velocity and temperature profiles (at 1000 s) obtained with the k − ω SST CFD model are shown. The CFD approach is able to highlight the strong stratification occurring in the cooler that may impact the heat exchange in this relevant part of the circuit.

In Figure 14, the adiabatic mixing temperature (weighted average of the equilibrium) in the loop is shown for the stabilized flow regime. Some deviation between the CFD model and Modelica is present, the CFD predicting around 2 K higher temperature in the cold and the hot legs. However, both models show that stability in terms of energy balance is established. The deviation can be attributed to the initial conditions and the heat loss modeling, as well as the impact of the turbulence models in the CFD part. The Re is equal to 1145, whose calculation is based on the average stabilized MFR values, whereas the Pr number is a property of the fluid, equal to 3.45. This corresponds to the stable regime on the stability map presented in Figure 15.

Fig. 8 Mass flow rate.

Fig. 9 Mass flow rate for the last 1000 s of the simulation.

Fig. 10 Temperature difference across the cooler.

Fig. 11 Temperature difference across the cooler for the last 1000 s.

Fig. 12 Velocity (sx) and temperature (dx) profiles at the cooler inlet at t = 1000 s for k-wSST turbulence model.

Fig. 13 Velocity (sx) and temperature (dx) profiles at the cooler outlet at t = 1000 s.

Fig. 14 Adiabatic mixing temperature for stabilized flow.

Fig. 15 Stability map.

6 Conclusion

In this paper the experimental results with distributed heat source obtained from DYNASTY natural circulation loop are used to set up and validate CFD simulations of the facility in OpenFOAM. The simulation results are also compared to a 1-D modeling approach (using Modelica object-oriented language) and analytical stability maps. The paper provides first complete assessment of a novel natural circulation facility, DYNASTY, aimed at analyzing the stability of fluids with internal heat generation. The simulation results in terms of the temporal behaviour of the mean thermal-hydraulic parameters, e.g. the mass flow rate and the temperature difference across the cooler section are compared to the experimental outcomes. The objective of the study is to assess the capability of the CFD model to reflect correctly the equilibrium steady - state of a natural circulation loop in presence of distributed heating. For the CFD part, three turbulence models are considered and it is shown that for the modeling of a natural circulation loop with distributed heating the k − ω SST turbulence model provides more accurate results in terms of mass flow rate prediction. The comparison of the CFD results to the experiment suggests, that the realizable k − ε turbulence model is not a suitable turbulence model for this case, most likely due to the fact that it is intended for high Re number flows, it is sensitive for initial conditions and in general less stable than other turbulence models; it overestimates the mass flow rate for the stabilized flow significantly. The adia-batic mixing temperature profiles provide further details to the thermal balance established in the loop. The temperature difference across the cooler is overestimated by the CFD model, due to the heat losses not being captured in the modelling.

Conflict of interests

The authors declare that they have no competing interests to report.

Funding

This research did not receive any specific funding.

Data availability statement

Data associated with this article cannot be disclosed.

Author contribution statement

All the authors were involved in the preparation of the manuscript. The first author carried out modelling and simulation part of the work, results interpretation, comparison to the experimental data and wrote the parts of the paper concerning the CFD model. The co-authors contributed substantially to the writing of the sections regarding the facility as well as the additional numerical modelling tools. All the authors have read and approved the final manuscript.

References

All Tables

All Figures

	Fig. 1 DYNASTY facility.
In the text

	Fig. 2 Schematic view of DYNASTY.
In the text

	Fig. 3 Heating and instrumentation in DYNASTY facility.
In the text

	Fig. 4 CFD model of DYNASTY.
In the text

	Fig. 5 DYNASTY mesh: closeup at inflation layers applied on water body.
In the text

	Fig. 6 DYNASTY stability map example: distributed heating configuration, water as working fluid.
In the text

	Fig. 7 1-D model of the facility.
In the text

	Fig. 8 Mass flow rate.
In the text

	Fig. 9 Mass flow rate for the last 1000 s of the simulation.
In the text

	Fig. 10 Temperature difference across the cooler.
In the text

	Fig. 11 Temperature difference across the cooler for the last 1000 s.
In the text

	Fig. 12 Velocity (sx) and temperature (dx) profiles at the cooler inlet at t = 1000 s for k-wSST turbulence model.
In the text

	Fig. 13 Velocity (sx) and temperature (dx) profiles at the cooler outlet at t = 1000 s.
In the text

	Fig. 14 Adiabatic mixing temperature for stabilized flow.
In the text

	Fig. 15 Stability map.
In the text