Numerical simulation of the turbulent convective buoyant flow of sodium over a backward- facing step

A forced convective and a buoyancy-aided turbulent liquid sodium flow over a backward-facing step with a constant heat flux applied on the indented wall is simulated. Linear eddy viscosity models are used for the Reynolds stresses. Turbulent heat fluxes are modelled with a single gradient diffusion hypotheses with two different approaches to evaluate the turbulent Prandtl number. Moreover, the influence of turbulence on heat transfer to sodium is also assessed through simulations with zero turbulent thermal diffusivity. The results are compared with DNS data from literature. The velocity and turbulent kinetic energy profiles predicted by all models are in good agreement with the DNS data. The local Nusselt number trend is qualitatively well captured, however, its magnitude is underestimated by all models for the mixed convection case. For forced convection, the heat transfer is overestimated by all heat flux models. The simulation with neglected turbulent heat transfer shows the best overall agreement for the forced convection case. For the mixed convection best agreement is obtained using a correlation to locally evaluate the turbulent thermal diffusivity.


Introduction
The influence of buoyancy on the heat transfer for very low Prandtl number, Pr, fluids like liquid metals -especially liquid sodium -differs from that in the case of fluids with a Prandtl number around unity or higher. Experiments of mixed convection with liquid sodium in a vertical pipe show an increase of heat transfer for buoyancy-aided mixed convection due to the increased advection of heat [1]. The same trend has been found by Niemann and Fröhlich [2] in a numerical investigation of a turbulent flow over a backward-facing step under buoyancy-aided mixed convection conditions. The present contribution investigates the influence of buoyancy on the flow structure and the heat transfer characteristic for liquid sodium by means of steady state RANS simulations. The turbulent heat flux is modelled with a Simple Gradient Diffusion Hypothesis (SGDH). Simulations are performed with different approaches to evaluate the turbulent thermal diffusivity, α t , including the case of α t = 0. The results are compared with DNS data [2]. The choice of the turbulence models for the Reynolds stresses is based on previous studies of Cotton and Jackson [1], who used the model of Launder and Sharma [3], including an additional source term from Yap [4] in the equation of the dissipation rate of turbulent kinetic energy. The choice of the turbulent heat flux model is based on a literature research focused on low Prandtl number fluids. The Reynolds analogy, which assumes a constant turbulent Prandtl number, Pr t , close to unity, cannot be applied to low Prandtl number fluids. Therefore, Pr t is locally evaluated on one side with the correlation of Kays [5] and on the other side with the model of Manservisi and Menghini (MM) [6]. This recently developed four-equation model takes into account the local difference between the turbulent thermal and hydrodynamical time scales. However, both approaches have been developed and validated for forced convection only. They have already been applied to a turbulent forced convection sodium flow over a backward-facing step at Re H = 10000 and Er = 1.5 [7].

Numerical approach
The present calculations consider an incompressible Newtonian fluid with constant thermophysical properties and no viscous dissipation. The effect of density variation has been thus neglected except in the body force of the momentum equations, where it has been accounted for through the Boussinesq approximation. As already mentioned, two linear k − ε models have been used to calculate the turbulent viscosity, ν t , namely the one of Launder and Sharma [3] with the correction of Yap [4] (LSY) and that of Manservisi and Menghini [6] (MM). The time averaged turbulent heat fluxes, u i θ, are evaluated with a SGDH as follows: where T stands for the time averaged temperature and Pr t denotes the turbulent Prandtl number. For the simulations with the LSY model, Pr t is computed by means of the Kays [5] correlation: This accounts for the diminishing influence of the turbulent heat transfer in the near wall region. The turbulent Péclet number in the above equation is defined as Pe t = ν t /ν · Pr, with the kinematic viscosity ν. The model of Manservisi and Menghini [6] solves two additional transport equations, namely one for the temperature variance, k θ , and one for its dissipation rate, ε θ . The ratio between the characteristic turbulent thermal,τ θ , and hydrodynamic time scale, τ u , can thus be locally evaluated, R = τ θ /τ u = (k θ )/(k θ ), where k stands for the turbulent kinetic energy and for its dissipation rate. For higher Prandtl number fluids the latter assumes a constant value of approximately R = 0.5. Three different local thermal time scales are used to compute the variation of the thermal diffusivity, α t , and thus of Pr t , over the wall distance: with the damping function f 1 and the thermal time scales for the near wall region τ θ,nw , intermediate region τ θ,int and asymptotic region τ θ,asym being: τ θ,asym = 0.9 · k/ .  Table 1: Calculation of the discretization error for the streamwise velocity and temperature based on the second grid level (medium grid).
In the above equations the turbulent Reynolds number, Re t = k 2 /(ν ) and the wall distance, based on the Kolmogorov velocity scale, y * = y · (ν · ) 0.25 /ν are used, where y is the distance in wall normal direction. For the forced convection simulations a coarser grid than for the mixed convection case has been used. The coarsening of the mesh only affects the streamwise direction of the transition from the end of the heater to the relaxation zone, and the relaxation zone itself. A non-dimensional distance from the wall of y + = yuτ ν ≈ 1, with the friction velocity, u τ has been used, as required by the low-Reynolds turbulence models and in order to achieve high quality results. Grid convergence has been checked on the basis of the skin friction coefficient and the local Nusselt number. For the mixed convection cases the Grid Convergence Index (GCI) method of Roache [8] has been used to quantify the numerical discretization errors in the domain, as shown in Table 1. The refinement ratio between the fine and the medium grid is r 12 , while r 23 stands for the refinement ratio between the medium and coarse one. The mesh is equally refined in streamwise and wall normal direction. The global observed order of accuracy is expressed by p. The standard deviation σ is given for the discretization error e and the grid convergence index GCI. The results shown in the paper refer to the second grid level with approximately 110000 cells.

Computational setup
The steady-state, two-dimensional simulations have been performed using OpenFOAM. The SIMPLE algorithm for the pressure-velocity coupling has been used. Blended schemes with an order of accuracy between one and two have been used for the spatial discretization. The inlet profile of a fully developed turbulent channel flow has been generated via a separate simulation, defined by the friction velocity Reynolds number of Re τ = uτ h ν ≈ 395 with identical cell distribution in wall normal direction and of second order accuracy in space. The Reynolds number based on the inlet bulk velocity, U b , and the step height, h, is The inlet section is located 4h upstream of the step. Downstream of it the channel height, H, equals 2h, which corresponds to an expansion ratio of Er = H/(H − h) = 2. The heated surface at the indented wall is L h = 20h followed by an adiabatic one of length L r = 20h. The other walls are adiabatic and the non-slip condition has been enforced at all surfaces. The molecular Prandtl number is Pr = 0.0088. Two Richardson numbers, have been considered, namely Ri = 0 and Ri = 0.2. In the previous equation,q is the applied wall heat flux, g the acceleration due to gravity, β the volume expansion coefficient and λ the molecular thermal conductivity. For the MM model following wall boundary conditions have  Figure 1: Sketch of the numerical domain.
x/h = 0 x/h = 10 been applied: At the inlet section k θ ≈ 0 and ε θ ≈ 0, while at the outlet section zero gradient condition for all calculated variables have been imposed. Figure 2 shows the streamlines for the forced convection, Ri = 0 (top), and the buoyancy-aided mixed convection, Ri = 0.2 (bottom) case. For the former case, the corner counter-clockwise rotating eddy extends up to x/h ≈ 1.5h while the clockwise rotating one up to x/H ≈ 9 [2]. For the latter case, the influence of buoyancy on the velocity profile reduces the clockwise recirculation zone and shortens the reattachment length behind the step, as already observed [2,9]. Moreover, the clockwise rotating vortex has approximately the same extension of the counter-clockwise rotating one. It lies above the latter, thus detached from the heated wall, which is in agreement with previous investigations [2,9,10]. As can be seen in Figure 3, already at x/h = 6 the velocity profile for Ri = 0.2 is inverted with  respect to that for Ri = 0. The reason is the buoyancy force that accelerates the fluid close to the heated wall, while the fluid away from it slows down due to continuity. Successively, on the adiabatic wall, the fluid close to the wall decelerates and the velocity profile tends to that of a fully developed channel flow. For both cases the velocity fields predicted by the different models are in good agreement with the DNS data [2]. However, along the adiabatic wall differences between the simulations using a turbulent heat flux model and the one with α t = 0 are apparent. The reduced total heat transfer of the LSY model with α t = 0 results in higher fluid temperatures close to the wall. Consequently, the higher buoyancy force implies a slower deceleration of the fluid. The profiles of turbulent kinetic energy across the channel at different x/h positions are shown in Figure 4. The differences between the RANS predictions and the DNS data are more pronounced than for the velocity field. For the case with Ri = 0, turbulent kinetic energy is generated in the shear layer between the flow in streamwise direction and the clockwise rotating vortex. Behind the step both turbulence models overestimate the magnitude of k in the near wall region up to approximately y/h = 1, resulting in a stronger mixing within the shear layer. This leads to a shorter reattachment length, x r , compared to the DNS (x r,MM ≈ 7.5h, x r,LSY ≈ 8h, x r,DNS ≈ 9h). With increasing distance from the separation point, the turbulent kinetic energy decreases because of the reduced velocity gradients, as shown in Figure 3.

Flow field
For the buoyancy-aided case a different development of k is shown in Figure 4. Due to the strong shear stresses turbulent kinetic energy is generated downstream of the step up to x/h ≈ 4. From there on, the buoyancy forces strongly modify the velocity profile compared to forced convection, by accelerating the fluid close to the wall. The resulting very small velocity gradients away from it cause the flow to laminarize, as shown in Figure 4 by the very small values of turbulent kinetic energy. Thus, the production of turbulent kinetic energy is suppressed due to the negligible strain rate magnitude. Towards the end of the heated wall, the increased buoyancy influence causes a sign inversion of the shear stresses away from the wall and a recover in the production of turbulent kinetic energy, as shown in Figure 4 by the profile of k. Downstream of the heater, the turbulent kinetic energy further increases due to the still changing shear stresses. After x/h ≈ 24 it decreases due to the reducing velocity gradients shown in Figure 3.   force leads to a flow relaminarization, as already shown in Figure 4, leading to negligible c f differences between both turbulence models. However, starting from x/H = 8, where turbulence production recovers, differences in c f values become apparent. The neglected turbulent diffusion of heat in the LSY model with α t = 0 results in a higher near wall temperature of the fluid and therefore in higher buoyancy and increased skin friction. The influence of the turbulent kinetic energy is reflected in the strong increase of c f towards the end of the heater. The magnitude of the increase in c f is directly linked to the heat transfer model used. According to Figure  5 (right) the model predicting the highest local Nusselt number shows the lowest increase in skin friction because the high thermal diffusion results in lower near wall temperatures and hence lower buoyancy forces. Downstream of the heater the skin friction decreases due to the deceleration of the flow near the indented wall and tends towards the c f value for Ri = 0.

Thermal performance
The heat transfer predicted by the RANS simulations are compared with DNS data [2] in Figure 5 (right). The local Nusselt number, Nu x , is defined in terms of the step height h and the difference between the local wall temperature and the inlet temperature, ∆T : For the forced convection case with Ri = 0 the LSY model with α t = 0 predicts the local Nusselt number in very close agreement with the DNS data in the region of the corner eddy as well as in the recirculation zone. Further downstream, where the thermal boundary layer develops, the turbulent heat fluxes become important, as shown by the differences between the DNS data and the LSY model with α t = 0. The simulations using a turbulent heat flux model are not capable to compute the heat transfer correctly. One reason for this behaviour is the underestimation of the corner eddy, which leads to an increased convective transport closer to the step and thus to larger Nu x . Furthermore, the higher level of turbulent kinetic energy within the recirculation zone directly enhances the heat transfer and leads to an overestimation of the turbulent heat fluxes with respect to the DNS. Downstream of the reattachment point the computed Nu x converges towards the DNS values for both the MM and the LSY&Kays model. In the case of buoyancy-aided mixed convection the total heat transfer increases compared to the forced convection case, even where turbulent production is impaired. The same trend has been experimentally observed for a pipe flow of liquid sodium with buoyancy-aided mixed convection [11] as well as buoyancy-aided laminar mixed convection of air over a backward-facing step [9,10,12]. Compared to the case with Ri = 0, the initial decrease of Nu x extends to a higher distance from the step, according to the greater extension of the corner eddy. From there, Nu x increases as a result of the additional momentum and energy transported from the non-negligible wall normal velocity towards the heated surface. Nu x reaches its maximum value where the wall normal velocity becomes negligible and thereafter decreases because of the growing thermal boundary layer. However, the steep increase and the peak value are not captured in their full extent by the models. One reason for the underestimation of the turbulent heat fluxes results from the near wall behaviour of the models. The Kays correlation underestimates the turbulent heat diffusion in the near wall region due to its asymptotic behaviour towards infinity. Therefore, the increasing heat transfer due to the impinging flow is not captured by this model. For the MM model the turbulent heat transfer depends on both, the turbulent kinetic energy and the damping function f 1 , where f 1 itself is dependant on the Kolmogorov time scale and the Prandtl number. Due to the relaminarization of the flow the turbulent kinetic energy, and therefore the Kolmogorov time scale is reduced. As a result, the damping function tends to zero and the turbulent heat transfer becomes negligible. The other reason for the lower heat transfer arises from the single gradient diffusion hypothesis, which provides an insufficient approximation of the streamwise turbulent heat fluxes. However, using the generalized gradient diffusion hypothesis suggested by Ince and Launder [13] in combination with the LSY&Kays only model showed minor improvements on the Nu x distribution. Using a more sophisticated heat flux model to further improve the prediction of the streamwise heat fluxes is therefore expected to have a positive influence on the heat transfer predictions.

Conclusion
Two linear k − models in combination with two heat flux models have been compared with DNS data [2] for the forced and buoyancy-aided mixed convection of a turbulent liquid sodium flow over a backward-facing step. Furthermore, the influence of turbulence on heat transfer has been examined through simulations with zero turbulent thermal diffusion. For the forced convection case, the eddy diffusivity has only a small influence on the local Nusselt number just downstream of the step. As soon as the flow reattaches the turbulent heat transfer increases. Both investigated turbulence and heat flux models have failed to predict the correct local Nusselt number distribution for the considered Reynolds number. By increasing the buoyancy force, the heat transfer has shown the same behaviour as for laminar flows for buoyancy-aided mixed convection, i.e. increased Nusselt number values compared to forced convection. Due to the acceleration the flow has laminarized downstream of the separation point. The turbulence has recovered with increasing distance from the step causing a steep increase in skin friction. For Ri = 0.2, the recirculation zone has been detached from the heated wall. At the same time the heat transfer has increased due to the additional convection and the impinging effect downstream of the recirculation zone. All models investigated predicted the influence of buoyancy on the flow field and the heat transfer characteristics qualitatively correct. However, the impinging effect is underestimated by both heat flux models used. For these cases, the turbulence models show good agreement with DNS data for the streamwise velocity as well as for the turbulent kinetic energy. At the relatively low Reynolds number considered here, the LSY&Kays model has shown similar results compared to the more complex model of Manservisi and Menghini. Further investigations at higher Reynolds numbers have to be performed to assess the performance of the models when turbulence plays a more prominent role and to further investigate the buoyancy-aided turbulent convection of sodium for this geometry.