Direct numerical simulations for liquid metal applications

Direct numerical simulations for liquid metal applications

6.1.1

I. Tiselj*, E. Stalio†, D. Angeli‡, J. Oder* *Reactor Engineering Division, “Jozˇef Stefan” Institute, Ljubljana, Slovenia, †Department of Engineering “Enzo Ferrari”, University of Modena and Reggio Emilia, Modena, Italy, ‡ Department of Sciences and Methods for Engineering, University of Modena and Reggio Emilia, Reggio Emilia, Italy

6.1.1.1

Introduction

6.1.1.1.1 The Navier-Stokes equations The phenomenon of turbulence eludes any definition and is usually introduced through a list of typical flow features. The most relevant property of turbulence for the present text is that it emerges as solution to the incompressible Navier-Stokes equations: !

r  U¼ 0 ! ! ∂ U
! 1 2! + U  r U ¼ rP + r U ∂t Reτ

(6.1.1.1)

where Reτ indicates the friction Reynolds number. When these equations are solved at sufficiently high Reynolds number and using suitable numerical techniques, numerical results of unsteady nature are produced, which show excellent agreement with measurements in incompressible Newtonian fluids (Pope, 2000). This approach is called direct numerical simulation or DNS. When we talk about “agreement,” one should not expect the same temporal development of the dependent variables in a selected point of computational domain and in the equivalent point of the experimental device. As solutions of Eqs. (6.1.1.1) are chaotic, only the same statistical behavior of the numerical solution and measured signal can be observed. From the mathematical and computational point of view, any minor change in boundary or initial conditions, numerical scheme, grid refinement, a different implementation of the same algorithms, etc, will result in a different instantaneous solution of the turbulent field after a sufficiently long time. Nevertheless, statistical properties of the solutions remain unchanged as long as we do not significantly compromise the accuracy of our simulations. The good thing is that statistical behavior of all these numerical solutions shows close agreement with measurements.

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors. https://doi.org/10.1016/B978-0-08-101980-1.00016-8 Copyright © 2019 Elsevier Ltd. All rights reserved.

220

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

6.1.1.1.1.1 Historical context As referenced in Moin and Mahesh (1999) and Kasagi and Iida (1999), the first turbulent flow and heat transfer simulations date back to the end of the 1980s. Particular attention was focused on the DNS of the near-wall flows, which is most relevant also for the present chapter focused on heat transfer. Another type of turbulent flow simulations is homogeneous isotropic turbulence, see Pope (2000) for a review, which is historically also very important for investigating turbulent cascade mechanisms of energy, but is not considered in the present text. Two main types of DNS heat transfer studies, both important for liquid metal applications, must be distinguished: the first is the passive scalar hypothesis (the scalar field T does not affect the velocity field) applicable for forced convection studies, and the second type of studies represents natural and mixed convection phenomena. The focus of this chapter is on simulations of near-wall turbulence as they reveal the basic mechanisms of the convective heat transfer between the fluid and the solid wall. In the past decade, this topic became important also in the field of thermal fatigue problems (Brillant et al., 2006; Aulery et al., 2012). It is important to stress that agreements between the DNS results and measurements in turbulent heat transfer studies cannot be as good as for the pure fluid dynamics problems, because the material properties of many fluids show nonnegligible variations with temperature. Most of the state-of-the-art DNS studies assume constant material properties and only some of them are considering the buoyancy effects. Thus, most of the results in this chapter are shown for the passive scalar approximation approaches. First, DNSs of heat transfer in the channel flow geometry were made at low Reynolds numbers and at Prandtl numbers of around one by Kim and Moin (1989) and Kasagi et al. (1992). Later, Kawamura et al. (1998) and Na and Hanratty (2000) performed DNS of the turbulent channel flow at Prandtl numbers up to 10. All these studies considered temperature as a passive scalar. Low Prandtl number DNSs of channel flows, which are relevant for liquid metal heat transfer and can be found in the database of Kawamura (2017), were performed at Pr ¼ 0.025 in 2003 and 2004. Tiselj and Cizelj (2012) created a DNS database at Pr ¼ 0.01, which is closer to the lead and sodium properties. Their simulations took into account the detailed heat transfer in the solid heated walls of the channel and were performed for variable combinations of the liquid and solid material properties.

6.1.1.1.2 The scalar transport equation When the flow is not isothermal, conservation of thermal energy implies the following dimensionless transport equation for the temperature field that is solved together with Eqs. (6.1.1.1) ∂T ! 1 r2 T + S + U  rT ¼ ∂t Reτ Pr

(6.1.1.2)

Direct numerical simulations for liquid metal applications

221

The dimensionless temperature T in Eq. (6.1.1.2) may be interpreted also as the concentration of a specific species and the Prandtl number is in that case replaced with the Schmidt number. The mass transfer studies are frequently met in the publications in the field of chemical engineering. In the literature, two archetypal cases of turbulent transport of thermal energy stand out, namely infinite channel flow and boundary layer flow (Li et al., 2009). This chapter focuses primarily on channel or duct flows, given their higher relevance to the specific field of application. Fully developed flows in ducts of constant or periodic cross-sections are often studied with the DNS approach. In such cases, it is convenient to assume the homogeneity of the flow also in the streamwise direction, although it is clear that some variables, such as pressure and temperature, are actually not homogeneous. However, it is possible to solve for a set of new variables, which are periodical by definition and can be used to recover the original, nonperiodic variables (see discussion in Section 6.1.1.3). Note that the enforcement of periodic boundary conditions (BCs) to represent homogeneity is a physical model, which is more accurate when the computational domain that corresponds to the periodic length imposed to the numerical solution is large enough to contain the largest relevant turbulent eddies of the physical case (see discussion of the characteristic scales in Section 6.1.1.1.3). For a fully developed duct flow the dimensionless energy equation written in terms of the fluid excess temperature θ ¼ T  Tm(x) becomes: ! ∂θ 1 ω U r2 θ  ¼ r  ðU θÞ + ∂t Reτ Pr Ω Um

(6.1.1.3)

where U is the instantaneous x-component of the velocity vector, ω/Ω is the dimensionless ratio between the normalized duct lateral surface and volume, and Um and Tm are the dimensionless bulk velocity and temperature, respectively. Equations are normalized with the channel half width h, friction velocity uτ, kinematic viscosity ν, and friction temperature Tτ ¼ qw/(uτρfcp) (qw is the heat flux from the wall into the fluid). It is worthy to note that the last term in Eq. (6.1.1.3) is a source term arising from the definition of θ, which closes the energy balance so as to enforce fully developed conditions while maintaining the periodicity of θ. When buoyancy effects are taken into account, the Boussinesq approximation is usually enforced (Gray and Giorgini, 1976). A typical form of the momentum equation for fully developed mixed convection flows is: !

! ∂ U
! 1 2 ! Gr ω! ! + U  r U ¼ rp0m + r U + 2θ g + ı ∂t Reτ Ω Reτ !

(6.1.1.4)

The forcing term along ı is an imposed pressure gradient in the streamwise direction, ! while g is the gravity unit vector. Note that Eq. (6.1.1.4) is written in terms of a modified pressure pm and the previously introduced fluid excess temperature θ, whose governing equation remains Eq. (6.1.1.3).

222

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

6.1.1.1.3 Length and time scales in turbulent flows For a correct setup of a DNS of a turbulent flow, it is crucial to derive estimates of the smallest and the largest length and time scales. The smallest length scale of turbulence (i.e., the dimension of the smallest turbulent structures), where mechanical energy is dissipated by viscosity, was introduced by Kolmogorov for homogeneous turbulence at high Reynolds numbers (see the textbook of Pope (2000) for background): η¼

 3 1=4 ν E

(6.1.1.5)

where ν represents the kinematic viscosity and E the dissipation rate of turbulent kinetic energy per unit mass k. Eq. (6.1.1.5) is useful for a posteriori evaluations of the turbulent dissipation. For an estimate, the expression: η  LRe3=4

(6.1.1.6)

is used, where L represents the macroscopic scale used to evaluate the Reynolds number. It is generally accepted that for performing a DNS, grid spacings need to be of the same order of magnitude as the Kolmogorov length scale (Coleman and Sandberg, 2010). While far from solid walls a model for the Kolmogorov scale (Eq. 6.1.1.6) might be appropriate, the calculation of the local Kolmogorov length scale in wallbounded problems requires the value of the turbulent kinetic energy dissipation ε. However, the calculation of such a quantity can be cumbersome, and, in addition, it can only be accomplished a posteriori, as it needs well-gathered statistics. For wall-bounded flows, the Kolmogorov scale can be modeled by η 1=4  ðκy + Þ δv

(6.1.1.7)

where the von Ka´rma´n constant κ is around 0.41, y+ being the distance from the closest wall in wall units, and δv is the viscous length ν/uτ, see the book by Pope (2000, p. 290). Eq. (6.1.1.7) can be considered valid for fully developed wall-bounded flows and also as a good indication for developing boundary layers. The smallest time scale that follows from Kolmogorov theory is: τη ¼

η2 ν

(6.1.1.8)

in wall-bounded flows it can be rewritten as τη 

δ2v 1=2 ðκy + Þ ν

(6.1.1.9)

It is clear that time steps in DNS must be shorter than the Kolmogorov time scale in order to accurately capture the flow dynamics.

Direct numerical simulations for liquid metal applications

223

The smallest scale for the scalar field was given by Batchelor (Pope, 2000) for high Pr flows: ηθ 

η Pr1=2

(6.1.1.10)

while for low Pr flows the Obukhov-Corrsin scale should be considered (Sagaut, 2006): ηθ 

η Pr3=4

(6.1.1.11)

Thus, for high Prandtl (Schmidt) number heat (mass) transfer, the smallest length scales of the thermal (mass concentration) field are smaller than the dimensions of the smallest vortices (Bergant and Tiselj, 2007) and this must be taken into account in true DNS studies. At Pr ¼ 1, both scales are of the same size, although some studies (Galantucci and Quadrio, 2010) indicate that accurate simulation of a passive scalar field at Pr around unity actually requires a finer resolution than the velocity field due to the higher intermittency of the scalar field. Nevertheless, tiny differences might be observed only in less relevant higher-order statistics. Because liquid metal flows are characterized by low Prandtl numbers, the smallest length scales of the scalar field are larger than the smallest length scales of the velocity field, which consequentially determine the required resolution of the DNS studies. This can be clearly shown in Fig. 6.1.1.1, which shows the very fine structures of streamwise velocity and larger structures in temperature fluctuations in the centerplane of a turbulent channel flow with Pr ¼ 0.01 (Tiselj, 2014). Further discussion of the resolution is reported in Section 6.1.1.2. The largest length scales in a given turbulent flow are in principle determined by the geometry of the domain. However, for external flows and for flows with inlet and outlet sections, these length scales cannot be easily determined. Simulations in the simplest wall-bounded flow configuration (i.e., flow between infinite parallel plates), also known as channel flow or Poiseuille flow, have shown that the streamwise and the spanwise lengthscales might be unbounded (Tiselj, 2014).

Fig. 6.1.1.1 Instantaneous fluctuation fields of streamwise velocity (left) and temperature (right) at the center of the channel. Reτ¼180, Pr ¼ 0.01. Horizontal coordinate: streamwise direction spanning from 0 to 13,751 wall units. Vertical coordinate: spanwise direction from 0 to 4523 wall units (Tiselj, 2014).

224

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

The largest time scales are also important, because they determine the time interval, which is needed to obtain statistically acceptable averaged values of turbulent quantities. Typical averaging times in domains with inlet and outlet are several dozens flow-through times, where the flow-through time is defined as the time needed by a fluid particle to travel along the whole domain. For fully developed flows, a time scale based on the mean velocity and the reference length should be considered. More about the largest length and time scales relevant for DNS is given in Section 6.1.1.2, discussing periodic BCs, and in Section 6.1.1.3, discussing autocorrelation functions, which are used to determine the required minimum length of periodic domains in DNS.

6.1.1.2

Direct numerical simulation techniques

The first DNS studies of homogeneous and near-wall turbulence were performed with spectral schemes. Authors of the first turbulent channel DNS studies, which were performed 30 years ago, are still maintaining their databases with open access data (Moser, 2017; Jimenez, 2017). Spectral schemes are known as the most accurate techniques for DNS studies and are computationally very efficient; however, they are limited to rather simple geometries (cuboid, cylinder) and BCs. A typical DNS approach for turbulent channel flow and heat transfer is based on pseudospectral schemes, using Fourier series in the x (streamwise) and z (spanwise) directions, and Chebyshev polynomials in the wall-normal y-direction (Tiselj et al., 2001b). Second-order accurate time differencing (Crank-Nicholson scheme for diffusive terms and Adams-Bashfort scheme for other terms) is used with maximum CFL numbers of approximately 0.1. The aliasing error can be removed by computing the nonlinear terms on a number of modes 1.5 times larger in each direction (Kim and Moin, 1989). DNS heat transfer results obtained with this type of codes can be found in various papers (Kasagi and Iida, 1999; Kasagi et al., 1992; Tiselj et al., 2001b, 2004; Bergant and Tiselj, 2007). Some of the DNS studies include also conjugate heat transfer, which accounts for the penetration of the turbulent temperature fluctuations into the heated wall (Tiselj et al., 2001a, 2013; Flageul et al., 2015). The channel flow database obtained at friction Reynolds numbers Reτ ¼ 180, 395, and 590 with liquid metal Prandtl number Pr ¼ 0.01 was prepared within the EU project THINS (Tiselj and Cizelj, 2012). Typical resolution (i.e., the distance between the collocation points of the spectral scheme) in DNS studies of turbulent channel flows at friction Reynolds numbers between 180 and 2000 (Hoyas and Jimenez, 2006; Moser et al., 1999) is Δx+ ¼ 1018, Δz+ ¼ 5, and Δy+ ¼ 8 (channel center), in the streamwise, spanwise, and wall-normal directions, respectively. This resolution guarantees the accuracy of

Direct numerical simulations for liquid metal applications

225

turbulent statistics for most of the relevant quantities. Nevertheless, a detailed resolution study performed by Vreman and Kuerten (2014) has shown that some higher-order statistics of turbulent flow (e.g., turbulent dissipation) require a finer resolution, namely Δx+ ¼ 6 and Δz+ ¼ 4, in the streamwise and spanwise directions, while they recommend Δy+ ¼ 3 in the center for the channel and Δy+ ¼ 1 in the near-wall region at y+ ¼ 12 in wall-normal direction. The resolution requirements proposed by Vreman and Kuerten are actually close to the theoretical prediction from Eq. (6.1.1.5), for which Δx+, Δy+, Δz+  η. For the choice of the time step size in DNS, the physical constraints and numerical accuracy and stability constraints must be met. In order to resolve the smallest time scales, the time step size should not exceed a fraction of the Kolmogorov time scale. However, numerical stability and accuracy criteria are often more stringent. For instance, in DNS of Tiselj et al. (2001b), with second-order accurate time scheme, a time step size is typically chosen such that the Courant number is around 0.1 to ensure the sufficient accuracy. The scheme is stable (but with poorer accuracy) for Courant numbers up to 0.5. Various finite volume and finite difference schemes were later developed for DNS studies. The reason why these methods are sometimes preferred to spectral schemes is that different types of BCs can be more easily implemented than in a spectral method framework. The finite difference method is thoroughly explained and documented in a couple of books (Fletcher, 1991; Anderson, 1995). In the framework of the finite difference method, compact schemes are currently used as high-order spatial schemes that increase accuracy and resolution without widening the computational stencil. The properties and techniques of the finite difference compact schemes are explained in a paper by Lele (1992). Starting in the 1990s, compact finite difference schemes have been used for DNSs. More recently, these schemes were implemented in Incompact3d, which is one of the most important Open Source codes for performing DNSs (Laizet et al., 2010). In a finite volume framework, balance equations for mass, momentum, and energy are expressed in their integral form over each of the nonoverlapping control volumes Ωi, which form the computational domain [ni¼1 Ωi ¼ Ω. Applicability of finite volume schemes on unstructured grids in complex geometries is one of the reasons why many commercial codes for industrial applications are based on a finite volume approach. The finite volume methods are usually presented in the literature restricted to the second-order accuracy form. Implementing high-order finite volume methods requires devising specific schemes for evaluating space averaged derivatives, interpolations, and fluxes, while deconvolution schemes are needed to extract pointwise values from spatial averages (Piller and Stalio, 2004). In summary, development of high-order finite volume algorithms is more demanding than development of such schemes for finite difference methods especially in complex domains (Piller and Stalio, 2008).

226

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

Finite difference method vs.

Finite volume method

Structured grids Noninherently conservative Close to the known equations Simple geometries Derivation and interpolation Higher-order spatial accuracy High-order methods

Unstructured grids Inherently conservative Physically meaningful Suitable for complex geometries Also integration and deconvolution Lower-order spatial accuracy Difficult for high-order methods

Probably the most frequently used DNS database obtained with finite volume schemes is being maintained since 1998 by Kawamura (2017). His team performed simulations of turbulent channel flow at various friction Reynolds numbers up to Reτ ¼ 1000 and Prandtl numbers between Pr ¼ 0.025 and 10. The recent DNS study of Vreman and Kuerten (2014) shows a detailed comparison of DNS results obtained with spectral and finite difference schemes. They found that their finite difference code (second-order spatial accuracy in periodic, fourth order in wall-normal direction, second order in time) required 3/4 smaller grid spacing than their spectral code to achieve the same accuracy. A number of other approaches have been also applied to DNS. We do not address the lattice Boltzmann method in this chapter, which has been shown to be capable of DNS and is known to be extremely efficient on massive parallel computers with GPUs due to its simplicity (Mayer and Hazi, 2006). However, the method, which is based on the Boltzmann transport equations, has not been widely used and thoroughly verified in single-phase heat transfer studies. It is instead worth to mention the spectral element method, developed in 1984 by Patera (1984). The method employs finite elements to discretize the domain and highorder Chebyshev or Legendre polynomials as basis functions within each element. It has gained a significant attention as a DNS tool thanks to the Open Source Code Nek5000, developed by the Argonne National Laboratory (ANL, 2017). It is known to be of similar accuracy as spectral schemes, with the advantage of being applicable to rather complex geometries. Practically speaking, when setting up a DNS from scratch, it is advisable to refer to previous works with similar geometries, for comparable flow regimes, conducted with similar numerical schemes, for the choice of spatial and temporal resolution. A simulation with a, for example, 2 times coarser grid and longer time steps is often run first in order to reach statistically steady conditions, saving computational resources. The computed fields can then be used as initial conditions for the resolved DNS, by interpolation on the final mesh.

Direct numerical simulations for liquid metal applications

6.1.1.3

227

Boundary and initial conditions

6.1.1.3.1 Periodic boundary conditions In hydrodynamically and thermally fully (time) developed conditions and also for periodically fully developed conditions (Stalio and Piller, 2007), where homogeneity of the flow in the streamwise direction can be assumed (see Section 6.1.1.1), periodic BCs can be implemented on all variables along the streamwise direction, including pressure and temperature. One of the most frequently used geometries in DNS, where periodic BCs are imposed, is the channel flow. Simulations of a turbulent flow in a channel of infinite length (x) and width (z) are performed by approximating the infinite extent of the domain in the streamwise and spanwise directions with a computational domain of finite length (also called unit cell) and by applying periodic BCs in these directions. Such conditions theoretically entail the dependent variable fields to be periodic along the selected direction. In practice, their implementation depends on the specific method itself: for instance, periodicity is inherent to the basis functions of spectral methods, while in finite volume/finite difference methods, it is achieved by imposing a circular topology of the computational stencils at the two ends of the domain. Special care should be taken when choosing the extent of the unit cell, which should be large enough to encompass the largest scales of the flow. If these conditions are satisfied, the flow in the computational domain can be considered as an accurate representation of the flow in the infinite channel. Section 6.1.1.4 reports some qualitative criteria for the choice of the channel length and width.

6.1.1.3.2 Inflow open boundary conditions Periodic BCs can be applied in the streamwise direction to the computed variables when fully developed conditions are met. In case the interest is a spatially developing flow, the DNS approach requires resorting to an inflow turbulence generation method for both the velocity and the temperature fields. Inflow turbulence generation methods for the velocity field have been recently classified and reviewed by Wu (2017). A short auxiliary periodic domain, where turbulence is computed in parallel to the main calculation, is providing results that are used as inflow condition for the target simulation. This type of inflow BCs is known as a recycling method. Synthetic methods represent an alternative attempt to fully model the inflow.

6.1.1.3.3 Outflow open boundary conditions Devising appropriate BCs at open boundaries when the mean flow is exiting the domain can be a difficult task (Sani and Gresho, 1994). The most common problems of outflow open BCs is that they produce spurious reflection of the velocity and of transported scalars, they tend to produce wiggles at the domain boundary and might affect the stability of the solution algorithm.

228

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

As outflow BCs have to represent the continuation of the physical domain beyond the edge of the computational domain, transport equations extrapolated from the domain interior are often enforced in the form of ∂Φ ∂Φ + Uc ¼0 ∂t ∂x

(6.1.1.12)

In many cases a uniform value for Uc ensuring mass conservation is used. An expression for the convective velocity Uc can be imposed in Eq. (6.1.1.12) also based on the specific flow configuration. For example, Craske and van Reeuwijk (2013), which deal with turbulent jets and plumes, suggest the use of an exponential function for Uc(y). In the review on open outflow conditions by Hattori et al. (2013), the method by Stevens (1990) is described as a promising approach for buoyant, turbulent plumes. The method proposed by Stevens (1990) is based on a one-dimensional advectiondiffusion equation at the outflow, where a diffusion term is added to Eq. (6.1.1.12) and a phase velocity is introduced in addition to Uc.

6.1.1.3.4 Boundary conditions for the thermal field Turbulent channel geometry is also the most frequently used geometry for studies of the near-wall heat transfer. The most accurate approach is conjugate heat transfer, where heat conduction inside the realistic heated walls is taken into account. A relevant study for liquid metal applications is the DNS of Tiselj and Cizelj (2012), who performed a conjugate heat transfer DNS at a low Prandtl number, Pr ¼ 0.01. Most of the DNS heat transfer studies (Kasagi et al., 1992; Kawamura et al., 1998; Na and Hanratty, 2000) were performed with simplified thermal BCs without solid walls and the dimensionless temperature at the fluid-solid contact plane that was fixed to zero: θðy ¼ hÞ ¼ 0 ðnonfluctuating thermal BCÞ

(6.1.1.13)

This type of thermal BC is a very good approximation of reality when the thermal activqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ity ratio K ¼ ðλρcp Þf =ðλρcp Þw goes to zero, while the heated wall thickness d+, and the parameter G ¼ αf/αw remains finite (α ¼ λ/ρcp). In that case, temperature fluctuations generated in the fluid do not penetrate into the solid wall. This type of ideal thermal BC was denoted as “nonfluctuating BC” by Tiselj and Cizelj (2012). An example of such system is air/metal combination of fluid and solid. For water/steel combination the thermal BC is still close to nonfluctuating BC if the heated wall is thick enough; however, the approximation fails when the metal wall thickness is small. The thermal BC, which permits maximum penetration of the turbulent temperature fluctuations from the fluid into the solid, is defined as: 0

hθðy ¼ hÞix, z, t ¼ 0

∂θ ∂y

! ¼ 0 ðfluctuating thermal BCÞ y¼h

(6.1.1.14)

Direct numerical simulations for liquid metal applications

229

Mean dimensionless temperature for the fluctuating temperature BC at the wall is still fixed to zero, while the zero derivative is prescribed for the fluctuating part of the temperature θ0 . This type of thermal BC is achieved when the thermal activity ratio K goes to infinity; it can be reproduced in experiments with a water flume that is heated with a very thin metal foil (Mosyak et al., 2001). In 2001 when the first analysis of the BCs (6.1.1.13), (6.1.1.14) in turbulent channel flow was performed by Tiselj et al. (2001b), these BCs were denoted as “isothermal”— Eq. (6.1.1.13)—and “isoflux BC for dimensionless temperature”—Eq. (6.1.1.14). Especially the “isoflux” name caused considerable confusion, because it was frequently mixed with standard isoflux BC for the physical temperature. These two BCs should not be mistaken: configuration with a constant volumetric heat source inside the solid walls corresponds to the isoflux BC for the physical temperature. The detailed BC in this case can be nonfluctuating (Eq. 6.1.1.13) or fluctuating (Eq. 6.1.1.14). Thus new names are recommended for BCs (6.1.1.13), (6.1.1.14) by Tiselj and Cizelj (2012). An important feature of different thermal BCs (6.1.1.13), (6.1.1.14) is that in the passive scalar approximation they have almost negligible influence on the heat transfer coefficient. In other words, mean temperature profiles near the wall are practically independent of the detailed thermal BC. Of course, a significant difference can be observed in temperature fluctuations and other turbulent statistics (see results in Section 6.1.1.4). For realistic fluid-solid systems the thermal BC is always in between the two limiting BCs of Eqs. (6.1.1.13), (6.1.1.14). For the liquid sodium and steel combination with K ¼ 1, the thermal BC is roughly in the middle between the fluctuating and nonfluctuating temperature BC. Hence, in order to study the details of the temperature fluctuations and their penetration into the solid wall, a conjugate heat transfer approach is recommended for DNS of heat transfer in liquid metals.

6.1.1.3.5 Initial conditions After all the equations, computational domain and BCs have been specified, initial conditions are needed to start the simulations. Given that long-term statistics are independent of initial conditions, and provided CPU time is not a limiting factor, it is often sufficient to start with any kind of initial velocity that contains at least a few perturbations of sufficiently long wavelengths. Because these initial stages of the simulations can be run on rather coarse grids, the CPU time consumption should not be too high. In most of the cases, if the geometry remains similar, it is sufficient to use instantaneous fields obtained from simulations at other Reynolds or Prandtl numbers after appropriate scaling. In case the computational cost is a limiting factor, then the initial transient can be shortened with some more sophisticated types of initial conditions (Coleman and Sandberg, 2010).

6.1.1.3.6 Statistical treatment of numerical solutions After a time interval, which typically spans over several flow-through times and 103 to 104 viscous time units, the flow in the computational domain should develop into a turbulent flow with constant mean properties. This condition is known as

230

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

“statistical steady state,” meaning that long-term statistics are independent of time. The assessment of such a condition is all but trivial. In principle, every flow parameter—velocity, pressure, or temperature at a given point, average value in a given line, on a selected surface or in a selected volume—should exhibit constant mean values with fluctuations around the means. In the turbulent channel flow, typical parameters used to identify statistical steady state are usually global variables, such as mean velocity, mean temperature, turbulence kinetic energy, friction velocities and friction temperatures computed at the walls, etc. Typically, statistical independence must be observed over several tenths of flow-through times and several thousands of viscous time units in order to recognize the statistical steady state in a running DNS. These runs in most cases require several 10–100 thousands of time steps. Once the statistical steady state is recognized, the spatial and temporal averaging of the results can start. Typical averaging procedures include computations of mean values in a point, over the line, surface, or in the selected volume. Flow statistics can be computed upon runtime or a posteriori by using snapshots of the field variables. Both approaches have obvious advantages and drawbacks. Because the number of statistical quantities to be computed can reach an order of a hundred parameters (when various budget terms must be evaluated), the choice of the postprocessing approach is influenced by the available computational and storage resources. Instantaneous fields can offer certain information; however, statistical averaging is necessary to compare results with experiments or other simulations. It is shown in the channel flow section, which is a geometry with two homogeneous directions, that the number of the time steps needed for the acceptable statistical uncertainty is rather large. The corresponding number of time steps for the flows with a single homogeneous direction is larger, and must be even larger in the geometries without homogeneous directions, where the only available type of averaging is averaging in time or ensemble averaging with independent DNS runs performed in the same geometry. The paper of Oliver et al. (2014) offers a methodology and a tool for analyses of uncertainties in DNS simulations. Their tool, which has been recently used in the work by Flageul and Tiselj (2017) to analyze statistical uncertainty in channel flow heat transfer, is recommended for quantification of the statistical uncertainties in the DNS studies.

6.1.1.4

Results: Channel flow

Overview of the turbulent heat transfer simulations performed in channel geometry at low Prandtl number is collected in Tables 6.1.1.1 and 6.1.1.2. The fourth column of Table 6.1.1.1 specifies thermal BCs used in the simulations: simulations denoted by F and NF were performed only with the two limiting thermal BCs: fluctuating and nonfluctuating temperature BC. Simulations denoted by C were conjugate heat transfer simulations. Simulations of Tiselj and Cizelj (2012) and Tiselj (2014) were performed at very low Prandtl numbers as benchmarks for the codes that will be used for heat transfer analyses in liquid sodium fast reactors. The latest publication by Oder et al. (2015) summarizes results of channel flow conjugate heat transfer simulations using the spectral element code Nek5000.

Direct numerical simulations for liquid metal applications

231

Table 6.1.1.1 Review of the DNS data at Pr ¼ 0.01 Num. scheme Spectral Spectral Spec. elements

Reτ 180, 395, 590 180 180

Thermal BC

K

G

d+

Reference

C

0.01–100

0.001–1000

0.1–1000

NF vs. F C

/ 1

/ 1

/ 180

Tiselj and Cizelj (2012) Tiselj (2014) Oder et al. (2015)

Notes: F, fluctuating temperature BC; NF, nonfluctuating temperature BC; C, conjugate heat transfer simulation.

Fig. 6.1.1.2 shows the mean velocity and the mean temperature profiles, results of the turbulent heat transfer DNSs performed at different Reynolds numbers and at low Prandtl number Pr ¼ 0.01 (Tiselj and Cizelj, 2012). These profiles are calculated with averaging of the velocity field over the planes parallel to the wall and over time. One can easily distinguish between the velocity and temperature profiles at different Reynolds numbers (180, 395, and 590) due to their different lengths. Despite the relatively low Reynolds numbers the velocity profiles exhibit a region at y+ > 50, which can be denoted as a quasilogarithmic layer. Temperature profiles in liquid metal are rather different and are more similar to parabolic profiles characteristic for laminar flows. One would probably need results at friction Reynolds number of 104 to see at least a small region similar to logarithmic profile in the mean temperature profiles. There are two temperature profiles in Fig. 6.1.1.2 for each Reynolds number: each for a separate limiting BC described by Eqs. (6.1.1.13), (6.1.1.14). One can see that heat transfer with fluctuating BC is slightly more efficient; however, this difference is negligible in the practical applications. This is in fact an important result concerning the conjugate heat transfer: the type of thermal BC, that is, fluctuating or nonfluctuating temperature BC, does not have a significant influence on the mean temperature profile and on the heat transfer coefficient. The study (Tiselj et al., 2004) has shown that the heat transfer coefficient for the fluctuating temperature BC is up to 1% higher (typically around 0.5% higher) than for the nonfluctuating temperature BC. However, this difference is difficult to measure even in the DNS simulations, because it is comparable with statistical uncertainty of the DNS results. Nevertheless, the profiles with fluctuating temperature BC are consistently around 0.5% higher than the nonfluctuating temperature BC at the same Reynolds and Prandtl numbers: 0.6% (Pr ¼ 1), 0.3% (Pr ¼ 0.01, Reτ ¼ 180), and 0.8% (Pr ¼ 0.01, Reτ ¼ 590). It remains to be confirmed that this difference remains small also in buoyant flows. The second most important parameter in turbulent heat transfer is probably temperature fluctuation profiles shown in Fig. 6.1.1.3 and obtained with similar averaging procedure as the mean temperature profiles. The left drawing in Fig. 6.1.1.3 shows profiles of RMS temperature fluctuations at Pr ¼ 0.01 obtained at Reτ ¼ 180, 395, and 590. Two profiles are shown for each Reynolds numbers for both limiting types of BCs (6.1.1.13), (6.1.1.14). A clear difference is seen between both BCs: temperature fluctuations approach to zero at y+ ¼ 0 for nonfluctuating temperature BC while

232

DNS

“Normal” “Large” “Very large” Reτ ¼ 395 Reτ ¼ 590

Domain

Mesh

Lx × Ly × Lz

Nx × Ny × Nz

Δx+

Δy+

Δz+

4π  2  4π/3 12π  2  4π 24π  2  8π 2π  2  8π 2π  2  8π

128  129  128 384  129  384 768  129  768 256  257  256 384  257  384

17.7 17.7 17.7 9.69 9.65

0.054–4.42 0.054–4.42 0.054–4.42 0.03–4.85 0.044–7.24

5.89 5.89 5.89 4.85 4.52

Aver.

Achieved

Δt+

time t+

Reτ

0.027 0.027 0.027 0.0158 0.0195

10,800 24,300 21,600 7900 11,800

180.09 179.99 180.04 394.65 588.95

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

Table 6.1.1.2 Details of various DNS, computational domain, grid, time step, averaging times of DNS, and friction Reynolds numbers achieved in simulations at Pr ¼ 0.01

233

25

3

20

2.5

F NF

2

15

qmean

mean

Direct numerical simulations for liquid metal applications

10

1.5 1

5

0.5

0 0.1

1

100

10

1000

0 0.1

1

10

y+

100

1000

y+

Fig. 6.1.1.2 Mean velocity and temperature profiles at Reτ ¼ 180 (shortest profile), 395, 590 (longest), Pr ¼ 0.01. F, NF, fluctuating and nonfluctuating thermal boundary condition, respectively.

0.06 0.25 F NF

0.05

0.2 0.04 qRMS

qRMS

0.15 0.1

0.03

F-normal domain F-large domain F-very large domain N-normal domain N-large domain

0.02

0.05

.036

0.01

0.033 50

0 0.1

100

0 1

10 y+

100

1000

0.1

1

10

100

y+

Fig. 6.1.1.3 Temperature RMS fluctuation profiles at Pr ¼ 0.01 for fluctuating (F) and nonfluctuating (NF) temperature BC. Left: At Reτ ¼ 180, 395, 590. Right: At Reτ ¼ 180 obtained on computational domains of different sizes. (Data and figures from Tiselj, I., Cizelj, L., 2012. DNS of turbulent channel flow with conjugate heat transfer at Prandtl number 0.01. Nucl. Eng. Des. 253, 153–160; Tiselj, I., 2014. Tracking of large-scale structures in turbulent channel with DNS of low Prandtl number passive scalar. Phys. Fluids 26 (12), 125111.)

the fluctuations approach to a constant value for fluctuating BC. Temperature fluctuations at the wall for fluctuating temperature BC are comparable with the maximum of the temperature fluctuations for nonfluctuating temperature BC that is typically measured at the distance around 20/Pr1/2 wall units away from the wall. Precise value of the wall fluctuations can be predicted only by DNS. Fig. 6.1.1.3 (left) obtained at low Reynolds numbers shows that values of temperature fluctuations at y+ ¼ 0 for fluctuating BC grow with increasing Reynolds number. This dependence is believed to be a consequence of nondeveloped thermal boundary layer at such low Reynolds numbers. Results obtained at moderate Prandtl numbers of around 1 show much weaker dependence of the wall temperature fluctuations on the Reynolds number variations. It is expected that this value becomes almost Reynolds number independent at sufficiently high Re, but that remains to be confirmed by higher Reynolds number DNS or very accurate LES (Large Eddy Simulation) studies.

234

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

Second specific behavior of the passive scalar turbulence at Pr ¼ 0.01 is shown in Fig. 6.1.1.3 (right): Solid line DNS results in this figure are obtained at Reτ ¼ 180 with the same resolution and in the same computational domain size as in the DNS described in Moser et al. (1999) (solid line). Other two profiles are obtained with the same spatial resolution but in larger computational domains: the “large” domain and the “very large” domain results were obtained in 3  3 and 6  6 larger computational domains in streamwise/spanwise direction, respectively. These large domains thus allowed tracking of larger turbulent structures. Fig. 6.1.1.3 (right) shows that temperature fluctuations at the wall in the “very large” domain are 20% larger than in the “normal” computational domain. This is another result that is not observed at moderate Prandtl numbers. The explanation lies in the very large turbulent structures, which are too weak to have any kind of influence on velocity field statistics and also on temperature field statistics at Pr around unity (Tiselj, 2014). However, in low Pr number flows the relative importance of these structures grows, because high thermal diffusion eliminates smaller structures much more effectively than the larger structures with longer wave lengths. It is a rather well known fact from mathematical physics, that diffusion typically dumps higher frequencies (short wave lengths) of any kind of wave motions more efficiently than lower frequencies (long wave lengths). It is important to emphasize that velocity statistics were not influenced by the enlarged computational domain and the mean temperature profiles at Pr ¼ 0.01 also remained unaffected. The same effect of computational domain size is visible (Tiselj, 2014) also at Reτ ¼ 395, but remains to be confirmed at higher Reynolds numbers. Two of the more important statistical quantities are streamwise and spanwise turbulent heat fluxes, as they are usually needed in development of the LES/RANS (Reynolds Averaged Navier-Stokes) turbulent models. Fig. 6.1.1.4 shows that turbulent heat fluxes in liquid metal flow at low Re represent only a small fraction of the total heat flux qw ¼ 1. Profiles in Fig. 6.1.1.4 show differences between the fluctuating

0.03

0.05 F-normal F-large F-very large N-normal N-large N-very large

0.03

F-normal F-large F-very large N-normal N-large N-very large

0.025 0.02 x,z,t

x,z,t

0.04

0.02

0.015 0.01

0.01

0.005 0 0

20

40

60

80

100 y+

120

140

160

180

0 0

20

40

60

80

100 y

120

140

160

180

+

Fig. 6.1.1.4 Streamwise (left) and wall-normal (right) turbulent heat flux profiles at Pr ¼ 0.01 computed in “normal,” “large,” and “very large” computational domains. Results are shown for ideal fluctuating (F) and ideal nonfluctuating (N) temperature boundary condition. (Figures from Tiselj, I., 2014. Tracking of large-scale structures in turbulent channel with DNS of low Prandtl number passive scalar. Phys. Fluids 26 (12), 125111.)

Direct numerical simulations for liquid metal applications 1

1

0.8

0.8

Fluctuating BC Nonfluctuating BC

0.6

Fluctuating BC Nonfluctuating BC

0.6

0.4

Rzz

Rxx

235

0.4

0.2

0.2

0

0

–0.2 0

1000

2000

3000

4000

5000

6000

x+

–0.2 0

500

1000

1500

2000

z+

Fig. 6.1.1.5 Streamwise (left) and spanwise (right) autocorrelation functions of fluctuating BC (solid lines) and nonfluctuating BC (dashed lines) temperature fields at y+ ¼ 91 in “normal,” “large,” and “very large” domain (Tiselj, 2014).

and nonfluctuating thermal BCs and influence of the domain size on the results, which is visible in streamwise flux. On the other side, the wall-normal turbulent heat flux Fig. 6.1.1.4 (right) is barely influenced by the domain size. Influence of the long-lived large turbulent structures is seen also in the streamwise autocorrelation functions of the temperature fields in Fig. 6.1.1.5, which are calculated as Z

Lx+

 θ ðξÞθ ðξ  xÞdξ 0

0

0

Rxx ðx, yÞ ¼ Z

Lx+ 0

0

0

z, t



(6.1.1.15)

θ ðξÞθ ðξÞdξ z, t

where integrals of θ0 are taken at fixed wall distance y while statistical averaging denoted by h i brackets is performed in spanwise direction z and time t. Similar equation is used for spanwise autocorrelations. Autocorrelation function is a convenient quantity than is used to determine the sufficient length of the periodic DNS domain: computational domain is sufficiently long (wide) if the autocorrelation function of all quantities drops to zero at half length (width) of the domain. As shown in Fig. 6.1.1.5 in the “normal” computational domain the streamwise autocorrelation functions drop to about zero at the half way of the channel length, while the spanwise autocorrelations point to nonnegligible anticorrelation at the half way of the channel width. This anticorrelation is more pronounced for the thermal field with fluctuating temperature BC. In the “large” and the “very large” domains the streamwise autocorrelation functions drop to zero, while spanwise autocorrelations are closer to zero, but still show a slight anticorrelation. Similar profiles for autocorrelations of components of the velocity field or pressure, which can be obtained in various DNS databases, show that “normal” size domain is sufficient for most of the DNS simulations. As a final remark on the size of the computational domain it should be noted that many of the small-scale turbulent phenomena can be observed even in the very small

236

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

periodic computational domain. Jimenez and Moin (1991) have shown that even domain with dimensions Lx+ ¼ 300 and Ly+ ¼ 100 is sufficient to correctly reproduce the mean velocity profiles, profiles of the velocity fluctuations, and small scale velocity spectra. Authors of the present chapter are not aware of similar test for passive scalars at low Prandtl numbers. Many other turbulent quantities, spectra of the temperature fluctuations, and budget terms of the turbulent quantities are relevant for development of the LES or RANS models or empirical heat transfer correlations. The reader is referred to the references in Table 6.1.1.1 for details about these statistics. For the conjugate heat transfer results, the reader is referred to Tiselj and Cizelj (2012): where channel flow results are available at Reτ ¼ 180, 395, and 590, and at Pr ¼ 0.01 for variable fluid-solid properties. All turbulent statistics are somewhere between the values provided by the limiting thermal BCs discussed earlier.

6.1.1.5

Results: Nonplanar geometries

When geometries become more complex, also the flow and temperature fields gain specific, relevant features that are not observed over planar walls. Errico and Stalio (2014, 2015) performed DNS of the fully developed low Prandtl number forced convective heat transfer in a channel of high aspect ratio, having one wavy wall and one flat wall. The flow is forced and no buoyancy effects are accounted for. The simulated friction Reynolds number is Reτ ¼ 282, corresponding to a Reynolds number of the bulk velocity and the hydraulic diameter, Re  18, 900. Three fluids of different thermal conductivities, corresponding to Pr ¼ 0.71, 0.20, 0.025 are considered in the research. The fluid flow displays separation, reattachment, and a shear layer downstream the wave peak, these are conditions relevant for turbulent heat transfer and passive scalar transport applications. Low Prandtl number fluids are characterized by comparatively low Nusselt numbers, but results reveal that heat transfer is enhanced by the wall undulation. A minimum Nusselt number is found in the flow separation region, where the recirculating bubble acts as a barrier to the advection of heat, a peak heat transfer rate is instead located in the flow reattachment region for all fluids investigated. A detailed investigation of the components of heat flux in vertical direction reveals that the main contribution to heat transfer is always to be ascribed to the mean advective term, at least within the range of parameters investigated. Profiles of turbulent heat flux in Fig. 6.1.1.6 show that turbulent transport of heat almost vanishes at Pr ¼ 0.025; as a consequence mean temperature profiles at those low Peclet number values are typically laminar. Instantaneous contours of the streamwise velocity component and three temperature fields corresponding to the different Prandtl numbers simulated are shown in Fig. 6.1.1.7, together with the fluctuation fields. It appears that the mean temperature profiles of laminar characteristics at Pr ¼ 0.025 derive from an unsteady temperature field, a small range of spatial scales and weak turbulent heat flux in vertical direction. The additional discussion in Errico and Stalio (2014) reveals that the use of a uniform

Direct numerical simulations for liquid metal applications

x = 0.3

x = 0.7

237

x = 1.3

x = 1.7

2 1.6 y 1.2 0.8 0.4 0 0

2

4

0

2

4

0

2

4

0

2

4

0.2 y

0.1 0 0

0.2

0.4

0.6

0.8

1 x

1.2

1.4

1.6

1.8

2

0

Fig. 6.1.1.6 Profiles of turbulent heat fluxes, u0 θ : solid line (Pr ¼ 0.71); dashed line (Pr ¼ 0.20); dash-dot line (Pr ¼ 0.025).

turbulent Prandtl number for the modeling of turbulent heat transfer for Pr ≪ 1 is not expected to be accurate, especially in flow separation conditions. Flow over a backward-facing step (BFS) is another geometry relevant for several applications. The paper by Oder et al. (2015) describes preliminary DNS results of thermal fluctuations in walls and in turbulent flow of liquid metal flowing over a BFS with finite dimensions and solid walls. The BFS geometry is shown in Fig. 6.1.1.8 and is as close as possible to the Kasola experimental device that is planned at Karlsruhe Institute of Technology ( Jaeger et al., 2017). The Pr ¼ 0.01 fluid is flowing from the narrower part into the wider part. The realistic heated wall is assumed behind the step. The temperature field is a passive scalar, which means that the natural convection is not simulated, although it can be relevant. For the inflow BC over the BFS, a fully developed turbulent velocity field with constant temperature is used. To obtain the inflow, a recycling BC is used. The values for velocity and temperature from a plane parallel and downstream from the inflow are imposed as the inflow BC. The streamwise component of velocity is scaled at this operation to ensure a constant mass flow rate. Simulations are performed with the Open Source Code Nek5000 (ANL, 2017), based on spectral element method. Strong 3D nature of the resulting temperature field makes the recirculating vortex behind the step rather complicate, although that cannot be clearly seen from Fig. 6.1.1.8, where the instantaneous field is shown in the plane of symmetry. It turns out that one of the key problems in DNS of statistically 3D turbulence is the very long computational times that are needed to obtain reasonably low statistical uncertainty. Discussion of the BFS simulations can be concluded with a reference to a recent paper by Niemann and Froehlich (2016) where DNS of vertical buoyant flow over the step was analyzed, although the results were not supported by the measurements.

238

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

Fig. 6.1.1.7 Snapshots of the instantaneous streamwise velocity component and temperature fields, and their fluctuations for various Pr values.

Fig. 6.1.1.8 Instantaneous temperature field in a flow over a backward-facing step with the heated wall under the step.

Very few works are found in the literature for the fuel bundle flow simulations concerning low Prandtl number fluids, such as liquid metals (Govindha Rasu et al., 2014). A first attempt at carrying out DNSs of fully developed, buoyancy-induced convection flows around a vertical rod bundle is done by Angeli and Stalio (2018). Finite volume computations are performed by an original, second-order accurate discretization technique based on the representation of arbitrarily shaped cylindrical boundaries on a

Direct numerical simulations for liquid metal applications

239

nonuniform Cartesian grid (Angeli et al., 2015). A unit cell of a triangular lattice of rods with a pitch-to-diameter ratio P/D ¼ 1.4 is considered as the reference geometry. A Prandtl number Pr ¼ 0.031 is chosen as the working fluid. A friction Reynolds number Reτ ¼ 550 and a Rayleigh number value of Ra¼ 5  105 were imposed, resulting in a bulk Reynolds number of about Reb ¼ 8600. Fig. 6.1.1.9 shows the overall structure of the flow in a subchannel. Statistics were computed on a unit flow cell corresponding to one-sixth of the base subchannel. Contours of the average streamwise velocity component ux and wall-bulk temperature

Fig. 6.1.1.9 Rod bundle DNS: (A) Primary flow cell on which statistics are computed (highlighted in red) and its contour, together with the definition of local abscissa γ used for the profiles. Contours of time-averaged (B) ux and (C) θw  θ, on the unit flow cell and (D) streamlines of the mean secondary flow components.

240

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

0

Fig. 6.1.1.10 Rod bundle DNS: Profiles of k, kθ , u0x θ , and εθ across the unit flow cell border.

difference θw  θ on the unit flow cell are reported in Fig. 6.1.1.9, alongside with streamlines of the cross-flow components. The buoyancy-induced flow in these geometries is seen to be characterized by nonzero time-averaged secondary flows. The circulation is driven toward the rods in the narrowest gap, while it is directed toward the center of the subchannel in the largest gap. Profiles of relevant turbulent quantities along the curvilinear abscissa γ are reported in Fig. 6.1.1.10. Overall, it is to be noted that once again the low Pr value determines very small fluctuations of θ, and an analogously small dissipation rate of temperature fluctuations. Near-wall peaks of turbulent kinetic energy and Reynolds stresses are instead found.

6.1.1.6

Conclusions

DNS is one of the tools that are being used for studies of turbulent heat transfer at low Prandtl numbers. As shown in the present chapter, DNS (at all Prandtl numbers) is used in rather simple geometries and at low Reynolds numbers. Nevertheless, this is exactly the type of flows where DNS results provide the key added value. The amount of information in even the simplest DNS database is extremely rich and the extraction of statistics and interpretation of results is often a difficult part of the job. Aside from the computational cost, the use of DNS in real-world applications is further discouraged by the huge amount of over-abundant data with respect to quantities of interest. Furthermore, uncertainties on boundary and initial conditions, which cannot be avoided in any numerical approach, do not justify the use of DNS in these cases. Of course, DNS in simple geometries can be a very useful tool for development of subgrid LES models and RANS transport equations in controlled numerical experiments. In this respect, DNS of liquid metal flows does not make any difference with respect to the DNSs of other fluids. The simplest type of passive scalar DNS in near-wall geometry reveals some peculiarities of the low Prandtl number ( 0.01) heat transfer: near-wall mean temperature

Direct numerical simulations for liquid metal applications

241

profiles at low Reynolds numbers accessible by DNS (Fig. 6.1.1.2) are similar to parabolic laminar flow profiles due to the high diffusive heat fluxes, which are typically much larger than turbulent heat fluxes, see Fig. 6.1.1.4. Further differences are in the temperature fluctuations at the liquid-solid wall: due to the comparable liquid and solid properties, turbulent temperature fluctuations penetrate into the solid walls much more efficiently than in the liquid water-solid systems. Another property of low Prandtl number liquid heat transfer is the nonnegligible role of the large-scale structures, which are practically irrelevant in the velocity field statistics. The mean temperature profiles remain unaffected by the large-scale turbulent structures, while their contribution to the temperature RMS fluctuations (Fig. 6.1.1.3) is not negligible and can account for up to 20% of the total fluctuations. The presence of the large-scale structures becomes visible in the measurable temperature statistics only at Prandtl numbers around 0.01 or lower. The observed large-scale thermal structures can be understood as echoes of the weak large-scale structures hidden in the velocity field. They should be seen also in the statistics of the LES studies. The focus of the state-of-the-art studies in the field of low Prandtl number fluids is nowadays on buoyant flows, and on slightly more complex geometries. A few examples are mentioned in this text: wavy channel flow, flow over a BFS, and buoyant flow in the subchannel. The last two cases are under development and represent the contribution of research groups of the authors to the EU H2020 project SESAME. However, two main problems remain: extraction of information from DNS and statistical uncertainties. It is not always clear what is the useful information in the huge amount of DNS data and it is not always easy to extract it. Nevertheless, the role of the DNS benchmarks, which are performed in such geometries, remains crucial for development and improvement of turbulence models in LES and RANS approaches.

References Anderson, J.D., 1995. Computational Fluid Dynamics: The Basics With Applications. McGraw-Hill, New York, NY. Angeli, D., Stalio, E., 2018. Direct numerical simulation of rod bundle at low Prandtl number, Computers and Fluids, submitted. Angeli, D., Stalio, E., Corticelli, M.A., Barozzi, G.S., 2015. A fast algorithm for direct numerical simulation of natural convection flows in arbitrarily-shaped periodic domains. J. Phys. Conf. Ser. 655(1). https://doi.org/10.1088/1742-6596/655/1/012054. ANL, 2017. Database. Available from: https://nek5000.mcs.anl.gov/ (Accessed January 2017). Aulery, F., Toutant, A., Brillant, G., Monod, R., Bataille, F., 2012. Numerical simulations of sodium mixing in a T-junction. Appl. Therm. Eng. 37, 38–43. Bergant, R., Tiselj, I., 2007. Near-wall passive scalar transport at high Prandtl numbers. Phys. Fluids 19 (6), 065105. Brillant, G., Husson, S., Bataille, F., 2006. Subgrid-scale diffusivity: wall behaviour and dynamic methods. J. Appl. Mech. Trans. ASME 73 (3), 360–367. Coleman, G.N., Sandberg, R.D., 2010. A Primer on Direct Numerical Simulation of Turbulence—Methods, Procedures and Guidelines. School of Engineering Sciences, University of Southampton.

242

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

Craske, J., van Reeuwijk, M., 2013. Robust and accurate open boundary conditions for incompressible turbulent jets and plumes. Comput. Fluids 18 (86), 284–297. Errico, O., Stalio, E., 2014. Direct numerical simulation of turbulent forced convection in a wavy channel at low and order one Prandtl number. Int. J. Therm. Sci. 86, 374–386. Errico, O., Stalio, E., 2015. Direct numerical simulation of low-Prandtl number turbulent convection above a wavy wall. Nucl. Eng. Des. 290, 87–98. Flageul, C., Tiselj, I., 2017. Impact of unresolved smaller scales on the scalar dissipation rate in direct numerical simulations of wall bounded flows. Int. J. Heat Fluid Flow 68, 173–179. Flageul, C., Benhamadouche, S., Lamballais, E., Laurence, D., 2015. DNS of turbulent channel flow with conjugate heat transfer: effect of thermal boundary conditions on the second moments and budgets. Int. J. Heat Fluid Flow 55, 34–44. Fletcher, C.A.J., 1991. second ed. Computational Techniques for Fluid Dynamics, vol. 1 and 2. Springer, Berlin. Galantucci, L., Quadrio, M., 2010. Very fine near-wall structures in turbulent scalar mixing. Int. J. Heat Fluid Flow 31, 499–506. Govindha Rasu, N., Velusamy, K., Sundararajan, T., Chellapandi, P., 2014. Simultaneous development of flow and temperature fields in wire-wrapped fuel pin bundles of sodium cooled fast reactor. Nucl. Eng. Des. 267, 44–60. Gray, D., Giorgini, A., 1976. The validity of the Boussinesq approximation for liquids and gases. Int. J. Heat Mass Transf. 19, 545–551. Hattori, T., Norris, S.E., Kirkpatrick, M.P., Armfield, S.W., 2013. Comparison of non-reflective boundary conditions for a free-rising turbulent axisymmetric plume. Int. J. Numer. Methods Fluids 72, 1307–1320. Hoyas, S., Jimenez, J., 2006. Scaling of the velocity fluctuations in turbulent channels up to Retau¼2003. Phys. Fluids 18, 011702. Jaeger, W., Hahn, T.S., Hering, W., Otic, I., Shams, A., Oder, J., Tiselj, I., 2017. Design and preevaluation of a backward facing step experiment with liquid metal coolant. In: The 17th International Topical Meeting on Nuclear Reactor Thermal Hydraulics, Xi An, China. Jimenez, J., 2017. DNS Database. UPM. Available from: http://torroja.dmt.upm.es/ftp/ channels/data/ (Accessed January 2017). Jimenez, J., Moin, P., 1991. The minimal flow unit in near-wall turbulence. J. Fluid. Mech. 225, 213–240. Kasagi, N., Iida, O., 1999. Progress in direct numerical simulation of turbulent heat transfer. In: Proceedings of the 5th ASME/JSME Joint Thermal Engineering Conference. American Society of Mechanical Engineers, San Diego, USA. Kasagi, N., Tomita, Y., Kuroda, A., 1992. Direct numerical simulation of passive scalar field in a turbulent channel flow. J. Heat Transf. Trans. ASME 114, 598–606. Kawamura, H., 2017. DNS database. Tokyo University of Science. Available from: http:// murasun.me.noda.tus.ac.jp/turbulence/poi/poi.html (Accessed January 2017). Kawamura, H., Ohsaka, K., Abe, H., Yamamoto, K., 1998. DNS of turbulent heat transfer in channel flow with low to medium-high Prandtl number fluid. Int. J. Heat Fluid Flow 19, 482–491. Kim, J., Moin, P., 1989. Transport of passive scalars in a turbulent channel flow. In: Turbulent Shear Flows VI, Springer-Verlag, Berlin, p. 85. Laizet, S., Lamballais, E., Vassilicos, J.C., 2010. A numerical strategy to combine high-order schemes, complex geometry and parallel computing for high resolution DNS of fractal generated turbulence. Comput. Fluids 39 (3), 471–484.

Direct numerical simulations for liquid metal applications

243

Lele, S.K., 1992. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 103, 16–42. Li, Q., Schlatter, P., Brandt, L., Henningson, D.S., 2009. DNS of a spatially developing turbulent boundary layer with passive scalar transport. Int. J. Heat Fluid Flow 30, 916–929. Mayer, G., Hazi, G., 2006. Direct numerical and large eddy simulation of longitudinal flow along triangular array of rods using the lattice Boltzmann method. Math. Comput. Simul. 72, 173–178. Moin, P., Mahesh, K., 1999. Direct numerical simulation: a tool in turbulence research. Ann. Rev. Fluid Mech. 30, 539–578. Moser, D.R., 2017. Database. DNS database. University of Texas. Available from: http:// turbulence.ices.utexas.edu/MKM1999.htm (Accessed January 2017). Moser, D.R., Kim, J., Mansour, N.N., 1999. Direct numerical simulation of turbulent channel flow up to Reτ¼590. Phys. Fluids 11 (4), 943–945. Mosyak, A., Pogrebnyak, E., Hetsroni, G., 2001. Effect of constant heat flux boundary condition on wall temperature fluctuations. J. Heat Transf. Trans. ASME 123, 213–218. Na, Y., Hanratty, T.J., 2000. Limiting behavior of turbulent scalar transport close to a wall. Int. J. Heat Mass Transf. 43, 1749–1758. Niemann, M., Froehlich, J., 2016. Buoyancy-affected backward facing step flow with heat transfer at low Prandtl number. Int. J. Heat Mass Transf. 101, 1237–1250. Oder, J., Urankar, J., Tiselj, I., 2015. Spectral element direct numerical simulation of heat transfer in turbulent channel sodium flow. In: 24th International Conference Nuclear Energy for New Europe, Portoroz, Slovenia. Oliver, T.A., Malaya, N., Ulerich, R., Moser, R.D., 2014. Estimating uncertainties in statistics computed from direct numerical simulation. Phys. Fluids 26 (3), 035101. Patera, A.T., 1984. A spectral element method for fluid dynamics: laminar flow in a channel expansion. J. Comp. Phys. 54, 468. Piller, M., Stalio, E., 2004. Finite-volume compact schemes on staggered grids. J. Comput. Phys. 197 (1), 299–340. Piller, M., Stalio, E., 2008. Compact finite volume schemes on boundary-fitted grids. J. Comput. Phys. 227 (9), 4736–4762. Pope, S.B., 2000. Turbulent Flows. Cambridge University Press, Cambridge. Sagaut, P., 2006. Large Eddy Simulation for Incompressible Flows. Springer, Berlin. Sani, R.L., Gresho, P.M., 1994. Resume and remarks on the open boundary condition minisymposium. Int. J. Numer. Methods Fluids 18 (10), 983–1008. Stalio, E., Piller, M., 2007. Direct numerical simulation of heat transfer in converging-diverging wavy channels. J. Heat Transf. 129 (7), 769–777. Stevens, D.P., 1990. On open boundary conditions for three dimensional primitive equation ocean circulation models. Geophys. Astrophys. Fluid Dyn. 51, 103–133. Tiselj, I., 2014. Tracking of large-scale structures in turbulent channel with DNS of low Prandtl number passive scalar. Phys. Fluids 26 (12), 125111. Tiselj, I., Cizelj, L., 2012. DNS of turbulent channel flow with conjugate heat transfer at Prandtl number 0.01. Nucl. Eng. Des. 253, 153–160. Tiselj, I., Bergant, R., Mavko, B., Bajsic, I., Hetsroni, G., 2001. DNS of turbulent heat transfer in channel flow with heat conduction in the solid wall. J. Heat Transf. Trans. ASME 123, 849–857.

244

Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors

Tiselj, I., Pogrebnyak, E., Li, C., Mosyak, A., Hetsroni, G., 2001. Effect of wall boundary condition on scalar transfer in a fully developed turbulent flume flow. Phys. Fluids 13 (4), 1028–1039. Tiselj, I., Horvat, A., Mavko, B., Pogrebnyak, E., Mosyak, A., Hetsroni, G., 2004. Wall properties and heat transfer in near wall turbulent flow. Numer. Heat Transf. A 46, 717–729. Tiselj, I., Oder, J., Cizelj, L., 2013. Double-sided cooling of heated slab: conjugate heat transfer DNS. Int. J. Heat Mass Transf. 66, 781–790. Vreman, A.W., Kuerten, J.G.M., 2014. Comparison of direct numerical simulation databases of turbulent channel flow at Reτ¼180. Phys. Fluids 26 (1), 015102. Wu, X., 2017. Inflow turbulence generation methods. Ann. Rev. Fluid Mech. 49, 23–49.