Direct Numerical Simulation of liquid metal forced and mixed convection in a square rod bundle

This paper reports results of Direct Numerical Simulations of fully developed, forced and mixed convection of Liquid Lead Bismuth Eutectic (LBE) around four vertical cylindrical rods arranged in a square lattice at Pr = 0.031. The solutions provided here are placed in a context where very few data are available for low-Prandtl flows, especially in mixed convection conditions, and can hence be useful for the development and validation of advanced turbulent heat transfer models. The equations are discretized using a Finite Volume implementation of a second order projection method. The irregular cylindrical boundaries are handled with an original Immersed Boundary technique. A single friction Reynolds number value is set (Reτ = 180) and buoyancy effects are accounted for by imposing the Rayleigh number, under the Boussinesq approximation. A Ra-value Ra = 2.5 ×104 is selected in order to investigate the main differences between forced and mixed convection at low Prandtl number values. Time-averaged velocity and temperature fields are shown, in order to discuss the main features of the flow and thermal fields. First order statistics are also presented, highlighting the effect of aiding buoyancy on turbulent heat and momentum transport.


Introduction and scope
Heat transfer around cylindrical rods is part of a highly multidisciplinary series of studies on liquid metal fast reactors, aimed at increasing their safety and improving the efficient use of resources.In this context, the application of Direct Numerical Simulations (DNS) fulfils the need for advanced tools and high-fidelity solutions.The advancements in computational methods and resources made it possible to overcome the relative complexity of the geometries involved and the challenges posed by the computational demands of DNS.Despite these efforts, several aspects concerning flow and heat transfer around bundles of heated rods still need clarification.
So far, cases of forced convection on vertical rod bundles have mainly been considered and not much has been documented about buoyancy effects.El-Genk carried out several experiments regarding convection of water in a uniformly heated square lattice, collecting data for different number of rods and variable pitch-to-diameter ratio [1,2], from which he found that the rod arrangement negligibly affects the overall Nusselt number in both forced and natural convection regimes.In his review, Meyer [3] pointed out that a natural mixing phenomenon from one subchannel to another, through the gaps between the rods, enhances the heat transfer to the cooling fluid by reducing the temperature differences.More recently, Angeli et al. [4] performed DNSs for fully-developed forced and mixed convection around a bare triangular rod bundle, under a constant pressure gradient, pointing out that buoyancy effects tend to decrease turbulent fluctuations, both increasing the bulk Reynolds number and decreasing the heat transfer rate.However, the need for further understanding of the effects of buoyancy on the flow and heat transfer in this kind of domains is still significant, because of the impact that these effects can have on the operating conditions in actual reactor cores, especially in emergency conditions such as DHR (Decay Heat Removal).Furthermore, data for liquid metals or more generally for low Prandtl numbers, relevant for Gen-IV reactors, are even scarcer.
In the following, data of fully-developed flow and heat transfer around a bare vertical, square rod bundle are presented, obtained by means of Direct Numerical Simulations in both forced and mixed convection conditions.Simulations are performed applying a friction Reynolds number Re τ = 180, considering Liquid Lead-Bismuth Eutectic (LBE) at Pr = 0.031 as the coolant fluid.Buoyancy effects are introduced by imposing a Rayleigh number Ra = 2.5 ×10 4 , corresponding to a Richardson number Ri ≈ 0.24.

Problem statement
A flow around 4 cylindrical rods is analysed.The rods are arranged in a square lattice, as shown in Figure 1(a).The rod diameter is denoted as D, while the centers are spaced by pitch P , see Figure 1(b).In the present case, the pitch-to-diameter-ratio is fixed at P/D = 1.25.Each bar is treated as uniformly heated by a constant heat flux q ′′ and the flow is considered fully developed.Periodic conditions enforced on the velocity field and the modified pressure field p m well delineate the physical situation.Furthermore, the temperature field needs to be normalized so that periodic BCs can be also set on a modified temperature-like variable θ [5].

Mathematical formulation
The incompressible Navier-Stokes formulation is adopted.The Boussinesq approximation is adopted, all the fluid properties being consistently assumed as constant, with the exception of density in the buoyancy term.
Accordingly, the continuity equation and the momentum equation can be written as: where T b (z) is the bulk temperature and the subscript "o" denotes reference condition.
In Equation 3, a denotes the temperature slope in z, which may be evaluated from an energy balance with q ′′ the imposed wall heat flux over the wetted perimeter ω, ṁ the mass flow rate and Ω the cross section.
With the gravity directed along z, (g 1 = 0, g 2 = 0, g 3 = −g), the first and the third terms in equation 3 may be incorporated into the pressure field giving place to a modified pressure: Due to fully developed conditions, the time averaged streamwise gradient of the modified pressure is constant [6], and, considering the pure natural convection requirement, such a constant equals zero.The heated rods investigated are characterised by a heat source of finite height and, as they are associated with a suitable cooling section, the flow induced by buoyancy is a non-zero mean flow [7].As a consequence, the flow regime under consideration cannot be viewed as pure natural convection, because in that case the hypothesis of full development would entail a zero mean flow rate.To overcome this issue, a forcing term is added which realises a fully developed mixed convection flow, representative of the actual flow configuration.
The modified pressure term can be split into a linear part P m , depending on the streamwise coordinate only, and a fluctuation p ′ m , which is set as periodic in the present numerical procedure: The pressure gradient term becomes: where σ denotes the constant pressure gradient due to external forcing of the flow (σ = −dP m /dz).The fluid excess temperature is defined as where T b (z) is the time averaged bulk temperature.Given the fully developed conditions, θ is independent of the streamwise coordinate [8].
Defining the temperature scale as ∆T ref = q ′′ L ref /λ and the velocity scale as u ref = τ w /ρ o , the momentum equation can be written in its dimensionless form: The term ω * /Ω * , which represents the dimensionless constant pressure gradient due to external forcing of the flow, has to be modified to take into account the possible excess or defect of the buoyancy source.In statistically steady conditions, the volume average of the fluid-excess temperature ⟨θ * ⟩ is time invariant and depends on the initial conditions only, in addition, ⟨θ * ⟩ is not zero, even though ⟨u * z θ * ⟩ = 0, by definition of θ under isoflux conditions.If ⟨θ * ⟩ ̸ = 0, the term Gr/Re 2 τ ⟨θ * ⟩ acts as an additional uniform momentum source.In order to ensure that the integral of the momentum source adds up to the same value both with or without the buoyancy source, the following constraint is imposed: where (ω * /Ω * ) ref is the reference value of the uniform momentum source when buoyancy is not active.
The energy equation, in its dimensionless form, can be written as where Pr is the Prandtl number and u * b denotes the mean bulk velocity.Values of the dimensionless parameters Re τ , Ra = Gr Pr and Pr imposed in the examined cases are reported in Table 1, together with the values of the bulk Reynolds number Re b , Richardson number Ri = Gr/Re 2 b and Nusselt number resulting from the DNS.The imposed value of Re τ = 180 has been selected.Purely forced and mixed convection conditions are considered, in order to investigate the effect of buoyancy on the mean flow and turbulent quantities.To achieve Ri ≈ 0.24 (i.e.non negligible buoyancy effects), a Rayleigh number of 2.5×10 4 has been set.The Prandtl number value is Pr = 0.031, corresponding to LBE at the reference temperature T o = 220 • C, see [4].
No-slip conditions ( u i | v = 0) and a constant heat flux (− λ∂θ/∂η| v = q ′′ ) are imposed at the rod wall, whilst periodic boundary conditions are enforced at all ends of the unit cell portrayed in Figure 1(b

Numerical method
The Navier-Stokes equations are treated in their conservative form, making use of the Finite Volume method discretization on Cartesian grids.A second-order Projection Method [9] is implemented to decouple the pressure field from velocity.The time-advancement scheme is second-order accurate and it is the result of a three-level scheme, which is semi-implicit (explicit in the streamwise direction z) for the diffusive terms, and explicit (using Adams-Bashforth) for the advective terms.For spatial derivatives a second-order central difference scheme on staggered non-uniform Cartesian grids is implemented.The modeling of irregular cylindrical boundaries on Cartesian grids is accomplished thanks to an original scheme developed for convective cases, see [10] and [11].The computational domain size and grid spacing are reported in Table 2.A constant grid spacing is set along the streamwise direction z, while the discretization is non-uniform along the two cross-flow directions (x and y).Considering the fully developed conditions outlined above, the bundle of rods can be assumed to have an infinite extension along the streamwise direction z.The conformity of the computational grid spacing to DNS standards is verified by a posteriori evaluation of the minimum Kolmogorov scale η min (see again Table 2).

Results and discussion
Some relevant results obtained are presented in this section.For the sake of conciseness, asterisks denoting dimensionless variables will be dropped starting from now on.
Firstly, contours of the average streamwise velocity component u z are displayed: in Figure 2 the different behaviour between the forced (Ri = 0) and mixed case (Ri = 0.24) is evident by observing the extent of the region of maximum velocity.Indeed, for Ri = 0 the streamwise velocity component u z peaks at the center of the subchannel and decreases its value moving towards the narrow gap.For Ri = 0.24 the contour area for the maximum velocity value expands towards the narrow gap, as in Figure 2(b).The previous observation are reinforced by the profile plots in Figure 3, for the narrowest and largest gap between two adjacent rods.The mean velocity reaches its maximum value at the center of the subchannel for both cases and it assumes almost the same value (Figure 3(b)), whereas in the narrow gap (Figure 3(a)) the maximum velocity is notably higher for the mixed case.
Another parameter that reflects this behaviour is the bulk Reynolds number: the values obtained from the simulations are Re b = 1785 for the forced case and Re b = 1850 for the mixed one, see Table 1.The two Reynolds numbers do not differ much from each other, even with the effect of aiding buoyancy, which is mainly due to the friction Reynolds number of both simulations being set at 180 (a relatively low value).Contours of the time-averaged cross-flow velocity magnitude u ⊥ = u 2 x + u 2 y are shown in Figure 4 buoyancy determines a decrease of the magnitude of u ⊥ .Moreover, in Figure 4 the streamlines of the cross-flow components are displayed: 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 for Ri = 0. Otherwise, for the mixed case, a different trend is observed: the circulation is directed towards the center of the subchannel for both large and narrow gap.
To better understand the effect of aiding buoyancy, a brief analysis of the fluid excess temperature field is reported.For conciseness, in Figure 5(a) the contour of the average flow excess temperature is shown for the forced case only, since the mixed convection case for Ri = 0.24 exhibited a similar qualitative trend.Considering the low Prandtl number of the working fluid, heat transfer is actually dominated by diffusion.The profiles for the narrow and large gap are presented as well in Figure 5(b)-(c): in general, for both cases, slightly higher values of θ are found in the narrow gap; in the large gap (Figure 5(c)) temperatures for Ri = 0.24 are visibly lower with respect to Ri = 0.These rather small differences can be appreciated in the slight decrease of the Nusselt number, from Nu = 7.60 for Ri = 0 to Nu = 7.50 Ri = 0.24 (Table 1).
In Figure 6(a) and 6(b), contours of the time-averaged turbulent kinetic energy are displayed.At first glance it is easily understood that for Ri = 0.24 there is a significant decrease of fluctuation intensity (up to ≈ 80%).The result is in agreement with other studies carried out with different Richardson numbers [12], which highlight how the presence of aiding buoyancy leads to a decrease in kinetic energy over the entire domain.Taking into account the profiles in  Figure 7, it can be pointed out that the decrease is more pronounced in the large gap region, see Figure 7(b).In addition, in Figure 7 the contributions of u ′ x u ′ x , u ′ y u ′ y and u ′ z u ′ z are reported: the decrease of fluctuation intensity can be mainly ascribed to the z-component u ′ z u ′ z .The low Prandtl value considered in this study also entails that the intensity of temperature fluctuations is generally very low, as heat transfer is dominated by diffusion.The effect of aiding buoyancy further reduces temperature fluctuations, in agreement with what has been observed so far for the fluctuating velocity field.A proof comes from the maps of temperature variance in Figure 8: a significant decrease in the values of k θ can be appreciated over the entire domain, since the peak value is reduced by ≈ 93% when aiding buoyancy is active.but even more so by focusing on the fluctuations related to turbulent kinetic energy.Moreover, the effects on the mean temperature fields is almost negligible.The Nusselt number reveals a slight deterioration of the heat transfer rate for Ri = 0.24.This can be ascribed to the reduction of turbulent fluctuations of both velocity and temperature.The data presented herein are also intended to serve as a benchmark for the development and validation of turbulent heat transfer models, both current and to come.

Figure 1 .
Figure 1.3D view (a) and cross-sectional sketch (b) of the computational domain; (c) onequarter of the 512 2 mesh in the x − y plane.

Table 1 .
Flow and thermal parameters set for the two DNSs performed and relevant integral values calculated.
The buoyancy force in Equation 2 is obtained by means of the Boussinesq approximation and the assumptions of the case as previously stated:

Table 2 .
Domain size and grid spacing for the present DNSs. ).