DIRECT NUMERICAL SIMULATION OF NATURAL, MIXED AND FORCED CONVECTION IN LIQUID METALS: SELECTED RESULTS

Selected results of three Direct Numerical Simulations are presented, on relevant test cases for the thermal hydraulics of liquid-metal-cooled nuclear reactors, encompassing a wide spectrum of turbulent convection regimes. The ﬁrst test case is a Rayleigh-Bénard cell at a Grashof number Gr = 5 × 10 7 , representative of the conditions in the unstably stratiﬁed layer of coolant in a reactor pool in both standard operating conditions and emergency situations, e.g. shutdown of the cooling system. The second case is the mixed convection in a cold-hot-cold triple jet conﬁguration, representative of liquid streams exiting from the core into the pool, and relevant for the modeling of thermal striping and thermal fatigue phenomena on the vessel containment walls. The third case is the fully-developed ﬂow in a vertical bare rod bundle with triangular arrangement and a pitch-to-diameter ratio P/D = 1 . 4 , in both forced and mixed convection conditions. These regimes respectively represent normal operation or decay heat removal conditions in reactor cores. The availability of these numerical databases will allow for an in-depth analysis of the turbulent ﬂow and heat transfer in liquid metals under different convection regimes, and is also relevant for the development, calibration and validation of turbulent heat transfer models.


INTRODUCTION
Convection heat transfer in liquid metals is a fundamental topic of relevant interest in the field of nuclear engineering.Turbulent flow and heat exchange phenomena taking place in a liquid-metal-cooled reactor are generally rather complex, due to the peculiar flow configurations.Recent studies improved the availability of reference data for this particular category of fluids, see for example ref. [1,2].Since the execution of detailed experiments is hampered by practical difficulties which can be overcome only with the design of specific experimental setups and procedures [3], the resort to highly resolved Direct Numerical Simulation is a viable alternative to provide reliable data for the development, calibration and validation of advanced turbulent heat transfer models, given the accuracy of models for nuclear applications based on the Simple Gradient Diffusion Hypothesis (SGDH) and the specific behavior at low Prandtl number is still a matter of active debate [4][5][6].
The a priori analysis in ref. [7] shows that in a moderately complex wall-bounded flow involving separation and reattachment and a one-way coupling between momentum and temperature fields (constant density), the turbulent Prandtl number distribution is far from uniform at Pr = 0.025, Pr = 0.2 and Pr = 0.71 and the irregularities are comparable.As a consequence, the uniform turbulent Prandtl number approximation is not expected to introduce errors that depend strongly on the molecular Prandtl number.However, results reported in ref. [7] and the model for the turbulent diffusivity in ref. [8] indicate that the space-averaged turbulent Prandtl number increases for decreasing molecular Prandtl numbers.Thus in case of use of a SGDH model, the selection of the most suitable turbulent Prandtl number is the main difficulty in low-Prandtl number representations of turbulent heat transfer phenomena.Instead, when a more detailed study is conducted, it should be considered that large-scale fluctuations are in general smaller for order one molecular Prandtl number with respect to the low-Prandtl number fluids, see [9].
In this work, a selection of results of Direct Numerical Simulations (DNSs) performed on three different test cases at low Prandtl number are presented, representative of different convection regimes.These are: (i) Rayleigh-Bénard convection for Pr= 0.025, representative of the conditions in the unstably stratified layer of liquid metal in a reactor pool, found for example in Lead-cooled Fast Reactors (LFRs) or in the vessels of Light Water Reactors (LWRs), see ref. [10,11]; (ii) a mixed-convection cold-hot-cold triple jet configuration at a Richardson number of 0.25, for Pr= 0.031, representative of liquid streams exiting from the core into the pool, and potentially giving rise to thermal striping phenomena on the containment walls [12]; (iii) fully-developed wall-bounded forced and mixed convection (Ri= 0.22) in a bare rod bundle configuration, representative of normal operation or decay heat removal conditions in the reactor core.DNSs in the above-described configurations and for low-Prandtl-number values are scarce or missing, at least to the authors' knowledge.For example, Rayleigh-Bénard convection has been extensively studied but mainly for Pr ∼ O(1) or higher.Numerical works on natural convection at low Prandtl number can be found for instance in ref. [13][14][15].However they are limited to cells of small aspect ratio or relatively low Rayleigh numbers.Instead, triple jet simulations reported in literature have been performed considering unsteady RANS or LES turbulence models, see for example ref. [16][17][18].These studies aim at reproducing experimental results reported in the two main experimental benchmarks [19,20].The only attempt of DNS in this configuration is reported in [21].However the spanwise grid resolution (3 points) employed is not sufficient to represent the three tridimensional features of the flow, therefore simulation in [21] can hardly be considered a DNS.Finally, the bare rod bundle configuration has been investigated in both forced and mixed/natural convection regimes in fluids with Prandtl number of unitary order or higher, see ref. [22][23][24].To the authors' knowledge no DNSs have been performed for a low-Prandtl-number fluid despite the relevance of this configuration for nuclear reactor applications.
The DNSs presented here are carried out in the frame of the Euratom H2020 SESAME project and in the PRACE project P-TURB (Prandtl number effect on turbulent Rayleigh-Bénard convection).Overall, computations took almost 25 million core hours.Results obtained represent a valuable database for further development and validation of turbulent heat transfer models by providing new and highly-accurate data for the thermal-hydraulic behavior of liquid metals in various flow configurations.

Flow configuration
Natural convection is investigated in a laterally unbounded Rayleigh-Bénard cell, see figure 1.
where H is the Rayleigh-Bénard cell height, ∆T is the temperature difference between the walls, g is gravity and β, α and ν are respectively the thermal expansion coefficient, thermal diffusivity and kinematic viscosity of the working fluid.The Reynolds number based on the free-fall velocity U f f = √ g β ∆T H and the Grashof number are defined as For the case under consideration, Ra has been set to Ra = 1.25 × 10 6 , while a value of Pr = 0.025 was enforced, representing liquid mercury [25].Velocity and time are made dimensionless respectively by U f f and H/U f f , while non-dimensional temperature θ is computed as θ = (T − T m )/∆T , where T m is the average temperature between the hot and cold walls.

Numerical method
Numerical simulations for the Rayleigh-Bénard case are performed using Incompact3d [26,27].It solves the incompressible Navier-Stokes equations together with a passive scalar transport equation on Cartesian grids.A third-order Runge-Kutta time integration scheme is used for the present study.As Incompact3d originally implements only passive scalar transport, the buoyancy term according to the Boussinesq approximation has been implemented for the present study.A validation test reported in [28] indicates that buoyancy forces have been implemented without errors.
To assess the validity of the Boussinesq approximation, the classical method by Gray and Giorgini [29] is used.The validity map relative to liquid mercury in selected conditions is reported in figure 2. It is shown that the Boussinesq approximation, when a 10% relative variation of thermophysical properties around their reference values is admitted, holds for a wide range of reference lengths L ref and temperature differences ∆T ref .
The domain size was chosen such that The number of grid points along each direction is 10 −5 .The computational grid is uniform along the homogeneous directions x and z, while in y-direction computational nodes are clustered at the walls using the clustering function embedded in Incompact3d, by setting a stretching factor β = 0.35, see ref. [26].Periodic boundary conditions are imposed along x and z directions, while no-slip and constant temperature conditions are set at y-normal walls.The discretisation in space and time has been checked a posteriori by comparing grid spacing and time-step against the turbulence smallest scales.Indicating by η K and τ K respectively the Kolmogorov length and time scales computed at the walls, grid spacing on the wall-parallel directions is ∆ π = 3.039 η K , while along the wall normal direction the maximum spacing is ∆ y = 0.904 η K ; the time step is largely smaller than the Kolmogorov time scale, ∆t = 0.129 τ K .The mesh refinement along the wall-normal direction has been set in order to reproduce instantaneous velocity and temperature profiles smoothly.Beside small-scale resolution also the capability of the computational domain to reproduce large-scale motions needs to be considered.In this work the horizontal extension of the domain has been selected by following considerations reported in ref. [15], where the same domain is used in a DNS performed at the same Prandtl number but smaller Rayleigh number, Ra = 10 5 .An a posteriori assessment of the large scale resolution by means of the twopoint correlations presented in ref. [30] suggests that domain enlargements should be taken into account for future DNS.However this task goes beyond the scope of the present work and requires computational efforts unachievable by the available computational resources.
To further assess the quality of the numerical procedure employed the Nusselt number computed as Nu = (d θ /dy) y=0 is compared against three consistency relations [31,32], which are exact equalities derived from the set of governing equations.These relations are where V indicates values averaged over the entire computational domain and ε and εθ are respectively the pseudo-dissipation rate of kinetic energy and the thermal dissipation, defined as Note that equation (4) involves velocity components and temperature, not only their fluctuating components.
The comparison between relations (3) and Nusselt number computed is reported in table I and indicates that the numerical procedure employed is suitable to study the case under investigation.Errors in the consistency relations obtained with present numerical algorithm, which employs sixth-order finite difference schemes, are comparable with respect to results obtained with spectral methods, see for example [33].Validity region of the Boussinesq's approximation for liquid mercury constructed using the method by [29].Reference conditions are T 0 = 20 • C and p 0 = 1 atm, and a 10% relative variations of thermophysical properties around reference values is admitted.

Selected results
In Figure 3  Figure 4 displays the instantaneous fields of velocity magnitude and non-dimensional temperature on a xy-plane.Due to the low Prandtl number selected, the temperature field is mainly diffusive, while the velocity field exhibits a wide-range of scales, see ref. [14].As already mentioned in the studied configuration coherent structures are found to extend from wall to wall, see figure 4(a).Rayleigh-Bénard convection is characterised by a Large Scale Circulation (LSC) regardless of Ra and Pr numbers, at least for the ranges studied so far, see ref. [34].In fluids with Pr ∼ 1 or higher the LSC originates by a cluster of smallscale plumes, as stated for example in ref. [35] and [36].In the present configuration at a very low Prandtl number, no clusters of small-scale plumes are observed and the LSC is composed by a single large-scale plume spanning the entire cell height.This change of the flow pattern for very low Prandtl numbers has a direct repercussion on the heat transfer regime as reported in [13] where a change in the scaling laws is identified for Pr = 0.35.
Statistics have been computed by considering that the flow is statistically homogeneous in time and along x and z directions.After the initial transient is washed out, velocity and temperature fields have been collected for a non-dimensional time period ∆τ = 100.In figure 5 the average temperature profile, as well as the root-mean-squared profiles of temperature and velocity fluctuations are reported.Figure 5(a) displays that the mean temperature field has a diffusive behavior and that the thermal boundary layer thickness, defined as δ θ = H/(2 Nu), covers almost one tenth of layer height.Fluctuations profiles show that outside the thermal boundary layer temperature fluctuations are essentially homogeneous, while velocity fluctuations are inhomogeneous.
Turbulence modeling is usually performed by two equations, eddy-diffusivity models, which involve the transport equations of quantities related to the turbulent kinetic energy k and its dissipation rate ε.When also temperature is a relevant flow variable the gradient-diffusion hypothesis requires the computation of a turbulent thermal diffusivity α t , and this is usually derived from the turbulent viscosity ν t by imposing a constant turbulent Prandtl number Pr t = ν t /α t .Figure 6(a) displays the profiles of k and the average pseudo-dissipation rate of turbulent kinetic energy, defined as As reported in general terms by [37] and tested in the present configuration [30], the difference between ε and the dissipation rate ε = 2/ √ Gr (s ij s ij ), where s ij represents the strain rate tensor, is negligible.Figure 6(a) reports also the turbulent viscosity ν t computed as prescribed in traditional two-equations eddyviscosity models ν t = C µ k 2 /ε, where C µ is a constant customarily set to C µ = 0.09.The profiles k and ε display that inhomogeneities are mainly confined to the region close to the walls.The turbulent kinetic energy oscillates from 0.95 to 1.12 in the region 0.07 < y < 0.93, while dissipation can be considered constant and equal to 0.02 in the range 0.04 < y < 0.96, see figure 6(a).
In Reynolds-averaged models, turbulent heat fluxes are computed assuming the gradient-diffusion hypothesis which prescribes that turbulent heat fluxes are directed down the mean scalar gradient ∇ θ and are proportional to ∇ θ through the turbulent thermal diffusivity α t .In the present case turbulent heat transfer occurs only in the vertical direction, which results in Turbulent thermal diffusivity computed by its definition ( 6) is compared in figure 6(b) against the α t profile obtained from α t = ν t /Pr t where Pr t = 1.5, as suggested by ref. [38] for fluids of Prandtl number Pr ∼ O(10 −2 ).This a priori analysis shows that in present case the combination of a two-equations eddyviscosity model and the assumption of a constant turbulent Prandtl number is expected to fail in reproducing the profile of turbulent thermal diffusivity, both in magnitude and shape.The inset of figure 6(b) displays also the profile of turbulent heat flux v θ for comparison.These data are expected to be useful for developing turbulence models specific for wall-bounded, buoyancy-driven flows at low Prandtl numbers.

Flow configuration
A sketch of the triple jet configuration is provided in figure 7. The jets have equal mean velocities U h = U c = U j but the temperature of the central stream T h is higher than that of the lateral jets T c .The jet exits are spaced by 2.5 d, where d is the width of jets slot.Liquid Lead-Bismuth Eutectic alloy (LBE) is the simulated fluid.A typical working temperature for LBE in nuclear reactor applications is   for which the Prandtl number is Pr = 0.031.Values of the Reynolds and the Grashof numbers are set to Re = U j d/ν = 5000 and Gr = g β(T h − T c ) d 3 /ν 2 = 6.25 × 10 6 , respectively.Note that Re is five times lower than in the PLAJEST experiment by Kimura et al. [20] and other previous numerical studies, see for example ref. [16,17].The lower Reynolds number selected for the present study allows for a detailed DNS investigation of the flow.Unlike the experiment reported in ref. [20] and the numerical studies in ref. [16,17], the Richardson number in the present configuration is Ri = Gr/Re 2 = 0.25, thus buoyancy effects are non-negligible.
For a reference temperature T ref = 220 • C, thermophysical properties of LBE can be calculated from correlations reported in ref. [39].By selecting the jet width d = 1.5 × 10 −2 m, which matches the order of magnitude of the PLAJEST experiment by Kimura et al. [20], the chosen value of the Grashof number issues ∆T = 72.7 • C. A validity map similar to the one displayed in figure 2 and constructed for the triple jet case shows that Boussinesq approximation holds also in this case, see ref. [30].

Numerical method
Numerical simulations for the triple jet case are also performed using Incompact3d.Besides the addition of the buoyancy term, inflow and outflow boundary conditions have been developed specifically, see ref. [28].
Periodic boundary conditions are set along the cross-flow y and the homogeneous z directions.On the inflow plane (x = 0) velocity and temperature are set equal to zero at the solid wall, while in the three slots u, v, w and θ are assigned using snapshots recorded from a precursor channel simulation.The mean dimensionless fluid excess temperature of the hot jet centreline is θ h = 0.5, while for the lateral jets θ c = −0.25.Thus, the net inflow of thermal energy is zero.The fluid excess temperature is defined so that buoyancy forces vanish for θ = 0 (T = T ref ).In the x = L x plane, outflow conditions are enforced on both the velocity and thermal fields.
Snapshots required to set inlet condition in the triple jet simulation are recorded from a DNS of the fully does not introduce significant errors as the jet slot mesh is able to reproduce wall-normal gradients with an acceptable accuracy.To decorrelate the three inflows, a phase-shift is introduced between the three jets, the time shift considered is more than twice the integral time scale of the streamwise velocity component.In this way, the inflow turbulence is realistically generated rather than reproduced by synthetic methods.Further details about precursor channel simulation can be found in paper [28].
As reported in ref. [40] the choice of an outflow condition which is both stable and accurate in spatiallydeveloping flows is a delicate subject.The outflow condition originally implemented in Incompact3d (see ref. [41]) leads to stability issues in the present configuration, as reported in the preliminary studies on the triple jet configuration [42].To overcome these stability issues, the outflow condition presented in ref. [43] and reviewed in the work [44] is selected for the present study.This condition is based upon the one dimensional advection-diffusion equation: where χ represents the kinematic or thermal diffusivity, U s is a vertical advection velocity and U p is a phase velocity, which is computed using the generic quantity φ: In definitions (8), φ indicates either a velocity component or the fluid excess temperature θ, superscripts indicate time levels and n x denotes the grid node at the outflow plane.∆t and ∆x stands for the temporal and spatial discretisation along the vertical coordinate.In addition, it is prescribed to clip the sum U p + U s so to ensure that where U j is the jets centreline velocity and U f is the maximum velocity related to the buoyancy effect, equal to Ri 0.5 in dimensionless form., where these values are computed in the location where the dissipation rate peaks.Grid spacing in the x and y directions equals ∆x = ∆y = 3.8 η min , while ∆z = 3.0 η min .The time-step used corresponds to ∆t ≈ 6.7 × 10 −3 τ min and it is selected on the basis of numerical considerations.The domain extension and in particular its ability to resolve the largest structures of motions in the present configuration has been assessed in a preliminary study on the triple jet configuration [42].A critical discussion of the domain spanwise size based upon two-point spatial correlations can be found in ref. [28].

Selected Results
The instantaneous fields of vertical velocity component and temperature on a z-normal section are displayed in figure 8.The three jets are observed to coalesce close to the inlet forming a singular ascending stream.Two regions of the flow can be identified: a mixing region and a far-field region.The mixing region is confined by x < 6 and is characterized by the bending of lateral jets towards the central one and by large instantaneous temperature differences, see figure 8.This results in a high heat transfer rate due to both the fluctuations induced by jets mixing (see also the contours of the turbulent heat fluxes in figure 10(b) and 10(c)) and the high thermal diffusivity of the fluid.On the other hand, the far-field region features a singular stream where temperature differences are almost negligible, thus this region can be considered only marginally affected by buoyancy forces.In figure 9, a selection of statistics is reported, values are given in their non-dimensional form, see ref. [28] for the mathematical formulation and reference quantities.For convenience, data are presented in a region smaller than the computational domain, (x, y) ∈ [0, 20] × [0, 10].The statistics are calculated considering that the flow is homogeneous in time and along the spanwise direction z.The time period over which statistics are computed spans 110 d/U j .In addition, symmetry or antisymmetry has been exploited depending on the variable considered: u, w and θ are symmetric, while v is antisymmetric.
Figure 9(a) displays the mean streamwise velocity component.Induced by buoyancy the central jet accelerates, while the lateral ones decelerate and deviate towards the centreline.Regions of negative streamwise velocity are observed at the sides of the cold jets.This is due both to continuity and buoyancy effects where θ < 0.About x > 6, the three jets are considered to be coalesced in a single stream.The average fluid excess temperature field is displayed in figure 9(d).Close to the inlets high temperature gradients are observed, while for x > 6 temperature variations are bounded by |∆θ| < 0.1.Thus, far enough from the inlet, buoyancy effects become negligible and the flow behaves as purely mechanical.
Also temperature fluctuations, which are presented in figure 9(e), undergo an increase in intensity together with a spreading in the cross-flow direction.Their peak is found where the lateral jets merge with the central one, see figure 9(a), while for x > 10 their intensity is largely reduced.This behavior can be ascribed to the production mechanisms, which are in turn related to the average temperature gradients.As discussed above temperature gradients are large in the mixing region and almost negligible in the far-field region.
In its conventional definition (e.g. for a boundary layer) the turbulent Prandtl number is written as Pr t becomes singular in each location where the denominator changes sign and displays an irregular behavior in separated flows [7] and complex flow fields like the one presented here.Figure 10(a) displays the distribution of log 10 |Pr t |, where the variable is drawn only in the interval 0.1 ≤ |Pr t | ≤ 10.Contours are plotted in the region (x, y) ∈ [0, 10] × [0, 5] in order to focus on the heat transfer phenomena occurring where the jets mix between each other.It appears that if by one side models including a uniform Pr t are not expected to correctly represent turbulent heat transfer in such complex flows, on the other side also non-uniform Pr t turbulent heat transfer models are expected to be difficult to devise.The irregular behavior of the turbulent Prandtl number is not observed only for specific values of the molecular Prandtl number, or when velocity and temperature fields are two-way coupled, see for example figure 13 in ref. [7] where a very irregular behavior of Pr t is calculated for Pr = 0.025, 0.2, 0.71.In addition figures 10(b) and (c) report turbulent heat fluxes in the vertical and cross-flow directions in the region (x, y) ∈ [0, 10] × [0, 5].Turbulent heat transfer is shown to substantially increase in the pool and peaks are found in the mixing region, x < 6.

Flow Configuration
The case under consideration is schematized in Figure 11.An infinite triangular lattice of rods of diameter D is considered, whose centers are spaced by a pitch P , with P/D = 1.4.Each rod is considered as uniformly heated with a constant heat flux q".In the present work, a rectangular cell consisting of 4 subchannels (highlighted in Figure 11(b)) is considered.The particular choice of the P/D parameter comes from ref. [45], where a pre-test analysis of an experimental campaign on flow blockage in an LBE-cooled triangular array of rods is performed.Based on the wide test matrix foreseen for their experiment, a single value of the friction Reynolds number has been chosen for the DNS, i.e.Re τ = 550.Forced and mixed convection conditions have been considered.In the mixed convection case, the Richardson number imposed is R i = 0.22.The Prandtl number is fixed to Pr = 0.031, corresponding to LBE at a reference temperature T ref = 220 • C, see ref. [39].

Numerical Methodology
The numerical technique adopted is based on a Finite Volume implementation of a second order Projection Method, following ref.[46].Time-discretisation in the conservation equations is performed according to a three-level scheme, which is semi-implicit (explicit in the streamwise direction x) for the diffusive terms and explicit Adams-Bashforth for the advective terms.Such a practice is second order accurate in time.
Spatial derivatives are approximated with second order central differences on staggered Cartesian grids, which can be non-uniform on the transverse y − z plane.A fast direct resolution of the discrete momentum and energy equations at each time-step is made possible by means of Approximate Factorization, while the Poisson problem associated with the pressure-velocity coupling [46] is solved through a fast Poisson solver, based on Matrix Decomposition.The modeling of arbitrarily irregular cylindrical boundaries on Cartesian grids is achieved, thanks to the original scheme developed by Barozzi et al. [47] and extended to convective problems by Angeli and Stalio [48,49].The technique involves a local modification of the computational stencil where boundary segments intersect the stencil arms.The method has been implemented in a parallel, highly efficient in-house code.A validation of the code against literature DNS data and an assessment of its computational performances can be found in ref. [49].
The computational domain size considered is 8πD h × 1.2D h × 2.1D h , where is the hydraulic diameter of a subchannel.For the present case, D h 1.16D.The domain in the cross-flow direction is such that it encompasses four subchannels.The number of control volumes along the three directions is 2048 × 448 × 256.A constant spacing has been kept along all directions, corresponding, in wall units, ∆y + = ∆z + ≈ 2.5, ∆x + ≈ 6.The intersection between the immersed circular boundaries and the Cartesian grid implies a variable distance between the wall and the closest cell center, leading, in the present case, to an average value of y + wall ≈ 0.3.An a posteriori assessment of the grid spacing yields ∆x = 3.9 η K,min and ∆y, ∆z = 1.6 η K,min , where η K,min is the Kolmogorov scale computed where turbulence dissipation rate is maximum.Concerning boundary conditions, the physical situation of a bundle of heated rods under the hypothesis of fully developed flow and heat transfer can be well represented through periodic conditions enforced on the velocity field and the modified pressure field p m .Furthermore, the temperature field needs to be normalized so that periodic boundary conditions can be also set on a modified temperature-like variable θ, see ref. [50].The complete derivation of the governing equation in fully-developed conditions and the determination of source terms in the momentum and energy equations for the present case is omitted here for brevity, and can be found in ref. [51].

Selected Results
Contours of the instantaneous wall temperature and streamwise velocity fields are shown in figure 12. Regions of low and high temperature are observed on the heated cylinder walls, corresponding to high-and low-speed streaks, respectively.Figure 11(c) schematizes the overall structure of the flow in a subchannel.As acquainted above, under the hypothesis of fully-developed flow, all the variables whose time-averaged behavior in the streamwise direction is known are normalized in order to be able to impose periodic boundary conditions in the streamwise direction and, hence, treat such direction as homogeneous [50].Therefore, statistics do not depend anymore on the streamwise coordinate and the average fields become 2D.Furthermore, exploiting the inherent symmetries of the domain, statistics can be reduced to the unit flow cell corresponding to one sixth of the base subchannel.
Contours of the average streamwise velocity component u x , wall-bulk temperature difference θ w − θ in a subchannel are reported in figure 13 for both Ri = 0 and Ri = 0.22, alongside contours with superimposed streamlines of the cross-flow components u ⊥ .The flow in these geometries is seen to be characterized by a secondary circulation whose intensity is in the order of 0.1% of the bulk velocity, this value is consistent with the prediction reported in [52] for comparable P/D ratios.The circulation is driven towards the rods in the narrowest gap, while it is directed towards the center of the subchannel in the largest gap.The effect of buoyancy is particularly visible when observing the distribution of velocity in the forced and mixed convection cases: in the mixed case, the hotter fluid in the narrow gap increases the flow rate in that region.Conversely, the intensity of cross-flow components is attenuated.
Profiles of the turbulent Prandtl number Pr t and turbulent thermal diffusivity along the border of the unit flow cell are displayed in figure 14(a): it can be observed that higher values of Pr t are reached in the case Ri = 0.22, with a slightly steeper near-wall slope.As shown in figure 14(b), the higher Pr t -values in the mixed convection case are mainly due to the decrease of turbulent diffusivity α t .Hence, this suggests that turbulent diffusion of heat moving from the wall to the bulk of the flow becomes lower when buoyancy is present.Both occurrences can be linked to the resulting average Nusselt number values, i.e.Nu = 12.82 and Nu = 12.65 for Ri = 0 and Ri = 0.22, respectively, indicating that the overall heat transfer, for the chosen parameter values, is slightly impaired by buoyancy.

CONCLUSION
Three Direct Numerical Simulations (DNSs) are presented with the aim of providing a database of statistics which may be useful for the development and validation of turbulent convection models specific for liquidmetal-cooled reactors.Resort to numerical reference data as experimental results is almost mandatory due to measurements difficulties in opaque liquids.DNSs reported are performed on configurations relevant to nuclear reactors and include forced, mixed and natural convection regimes.The cases studied are (i) natural convection in a laterally-unbounded Rayleigh-Bénard cell at Pr= 0.025, (ii) mixed convection in a cold-hotcold triple jet configuration characterized by Ri= 0.25 and Pr= 0.031, and (iii) fully-developed wall-bounded forced and mixed convection (Ri= 0.22) in a triangular bare rod bundle configuration at Pr= 0.031.Rayleigh-Bénard convection at such a low Prandtl number exhibits coherent plumes which are found to extend from wall to wall and represent the Large Scale Circulation of the cell.As expected the thermal boundary layer covers a large portion of the domain, almost one tenth per side.The a priori analysis of the turbulent thermal diffusivity shows that the combination of traditional two-equation eddy-viscosity models with the assumption of a constant turbulent Prandtl number leads to significant errors in reproducing the vertical turbulent heat flux.The three jets in mixed convection regime are observed to undergo a region of intense mixing close to the inlet, where fluctuations of both velocity components and temperature reach their maximum value.Sufficiently far from the inlets, for x > 6d, jets are coalesced in a single stream characterized by very small temperature gradients, thus the ascending stream can be treated as isothermal.
The very irregular turbulent Prandtl number distribution indicates that the development of detailed turbulence models is a difficult task in this configuration.The gap between the rod bundle setup analyzed here is characterized by a secondary flow in the order of one thousandth of the average streamwise velocity.The comparison between the mixed and forced convection cases shows that buoyancy reduces the intensity of the secondary flow while increasing the the flow rate in the narrow gap region.In addition the turbulent heat transfer between the bulk and the rods is reduced when buoyancy forces are present, resulting in a slightly lower Nusselt number.
The statistics reported are expected to provide insights into the major flow features of turbulent heat transfer in the configurations considered and support the development and validation of turbulence models specific for nuclear rectors applications.Results obtained in two out of three cases (triple jet and rod bundle) are available for the scientific community at the website [53].

Figure 2 .
Figure 2. Validity region of the Boussinesq's approximation for liquid mercury constructed using the method by [29].Reference conditions are T 0 = 20 • C and p 0 = 1 atm, and a 10% relative variations of thermophysical properties around reference values is admitted.
instantaneous fields of vertical velocity v and non-dimensional temperature θ on horizontal sections at the edge of the thermal boundary layer y/H = 1/(2 Nu) = 0.094 are reported.The temperature field in figure3(a) displays four large regions of warm fluid, connected to each other through a net of high-temperature sheet-like plumes.Beside these regions, cold fluid portions indicate the presence of cold plumes in proximity to the lower (hot) wall, suggesting that in the configuration studied plumes cross the entire fluid layer.The velocity field in figure3(b) reflects the temperature distribution as it shows that warm fluid is rising while cold plumes are falling.

Figure 7 .
Figure 7. Triple jet case: three-dimensional representation of the flow configuration.
The simulation is performed on a domain of dimensions L x × L y × L z = 30 d × 30 d × 6 d, discretised by 2049 × 2048 × 512 equally-spaced grid points.The computational time step is ∆t = 0.0005 d/U j .The smallest turbulence scales computed a posteriori on the flow domain are: Kolmogorov length scale η min /d = 3.86 × 10 −3 and Kolmogorov temporal scale τ min U j /d = 7.43 × 10 −2

Figures 9 (
Figures 9(b) and (c) show fluctuations of the velocity components.Fluctuation intensities of the flow exiting the three slots rapidly evolve to much larger values and exhibit peaks where the three jets coalesce.Velocity fluctuations in y-direction undergo a larger increase with respect to fluctuations along x.

Figure 11 .
Figure 11.Bare rod bundle DNS: (a) 3D representation of a triangular bundle; (b) cross-flow layout with periodic rectangular module highlighted, comprising four subchannels; (c) unit flow cell.

Figure 12 .
Figure 12.Bare rod bundle DNS at Re τ = 550, Ri = 0.22: instantaneous contours of dimensionless wall temperature and cross-flow contours of streamwise velocity component.