A New Anisotropic Four-Parameter Turbulence Model for Low Prandtl Number Fluids

Abstract

Due to their interesting thermal properties, liquid metals are widely studied for heat transfer applications where large heat fluxes occur. In the framework of the Reynolds-Averaged Navier–Stokes (RANS) approach, the Simple Gradient Diffusion Hypothesis (SGDH) and the Reynolds Analogy are almost universally invoked for the closure of the turbulent heat flux. Even though these assumptions can represent a reasonable compromise in a wide range of applications, they are not reliable when considering low Prandtl number fluids and/or buoyant flows. More advanced closure models for the turbulent heat flux are required to improve the accuracy of the RANS models dealing with low Prandtl number fluids. In this work, we propose an anisotropic four-parameter turbulence model. The closure of the Reynolds stress tensor and turbulent heat flux is gained through nonlinear models. Particular attention is given to the modeling of dynamical and thermal time scales. Numerical simulations of low Prandtl number fluids have been performed over the plane channel and backward-facing step configurations.

1. Introduction

Liquid metals with their low Prandtl number have gained increasing attention in recent years. Compared with other coolant fluids, such as air or water, liquid metals provide large heat fluxes and can withstand high temperatures. Furthermore, some liquid metals, like sodium, can flow in the liquid phase at a wide range of temperatures without need for high pressurized systems [1]. Due to these properties, liquid metals are currently considered in a broad range of industrial applications, including the production of steel and semiconductors, in thermal solar plants [2,3] and in Generation IV nuclear power plants [4,5,6], i.e., the Lead Fast Reactor (LFR) and the Sodium Fast Reactor (SFR).

In a nuclear context, thermal-hydraulics is recognized as one of the key issues in the design and construction of liquid metal-cooled reactors. Since the possibilities for detailed measurement of local flow parameters in liquid metal cooled reactor components are challenging [7], numerical simulations of flow configurations are more important for low Prandtl number fluids than in usual cases. In this respect, Computational Fluid Dynamics (CFD) is regarded as a valuable tool to analyze the thermal-hydraulics behavior of nuclear systems. The more challenging aspects of the thermal-hydraulics of these systems are the low Prandtl number of liquid metals at operating conditions, the non-negligible buoyancy effects in the flow, and the significant turbulence anisotropy [8]. For these reasons, very sophisticated models are required to accurately simulate turbulent liquid metals flow and heat transfer.

In the Reynolds-Averaged Navier–Stokes (RANS) framework, several models have been developed in past decades for the computation of the Reynolds stress tensor $\overline{u_i^{\prime }{u}_j^{\prime }}$. The first-order models are based on the isotropic eddy diffusivity ${\nu }_t$ while the second-order models use transport equations for each component of the Reynolds stress tensor. From an academic point of view, the second-order models would be the best modeling for anisotropic momentum transfer, but the usage of these techniques requires a considerably increased numerical effort [9]. In this work, for the closure of the momentum equation, we propose an Explicit Algebraic Shear stress Model (EASM) that belongs to a class of models between first and second order. This class of models is derived from the second-order transport equation with the hypothesis of local equilibrium between production and dissipation, therefore convection and diffusion terms are neglected. The remaining closure terms are formulated in terms of the turbulent kinetic energy k and its dissipation rate $\epsilon$.

For the turbulent heat flux $\overline{u_i^{\prime }{T}^{\prime }}$, only a restricted number of models have been developed and validated. Most of them are first-order models based on the Simple Gradient Diffusion Hypothesis (SGDH) that assumes the similarity between the turbulent heat flux and the molecular heat conduction introducing the turbulent thermal diffusivity ${\alpha }_t$ and the turbulent Prandtl number $P{r}_t$, which is often set equal to a constant value in the range 0.8–1. This concept can reproduce reasonable results in the forced convection regime and for fluids with $Pr\approx 1$ whereas it is inadequate for applications involving non-unity Prandtl number fluids like liquid metals and/or non negligible buoyancy effects [10]. For these applications, the most promising models require the introduction of additional transport equations. In [11,12], an implicit Algebraic Heat Flux Model (AHFM) model has been proposed and implemented in STAR-CCM+. Its closure requires one additional transport equation for the evaluation of the temperature variance $\overline{T^{\prime 2}}$. In [11], the thermal model has been coupled with a low-Reynolds linear k-$\epsilon$ model, while in [12], the coupling with a second-order Reynolds stress model has shown better results. In [13,14,15,16], an isotropic four-parameter model has been proposed. The model introduces two additional thermal transport equations for the evaluation of the squared temperature fluctuations $k_{\theta }$ and its dissipation ${\epsilon }_{\theta }$. In the original formulation of this model, the turbulent heat flux is evaluated with an SGDH approach. In this work, we propose an anisotropic version of the above-mentioned four-parameter model and suggest an Explicit Algebraic Heat Flux Model (EAHFM) for the modeling of the turbulent heat flux. The thermal model is coupled to an Explicit Algebraic Shear stress Model (EASM) for the dynamical turbulence. For the closure of $\overline{u_i^{\prime }{u}_j^{\prime }}$ and $\overline{u_i^{\prime }{T}^{\prime }}$, we solve the four-transport model equation k-$\epsilon$-$k_{\theta }$-${\epsilon }_{\theta }$.

To validate the proposed anisotropic four-equation turbulence model (A4P), we consider two benchmarks. First, we simulate the plane channel geometry, a simple configuration widely studied in the literature. For this case, a database of DNS data is available for different $R{e}_{\tau }$ and $Pr$ number [17,18,19,20]. Then, we consider a more complex configuration, such as the flow over a backward-facing step in forced convection. Additionally, for this configuration, several studies are present in the literature and DNS data are available for $Re=9610$ and $Pr=0.088$ [21,22]. Both configurations have been tested with the isotropic four-parameter turbulence model. Plane channel simulations for several $R{e}_{\tau }$ and $Pr=0.01,0.025$ have been successfully performed [13,14,16]. Numerical simulations in forced and mixed convection regimes are very promising for the backward-facing step configuration [15], however the adoption of an anisotropic formulation is required to improve the prediction of turbulent heat flux components.

In Section 2, we present the new anisotropic four-parameter turbulence model, introducing the governing equations and the explicit algebraic expressions for both the Reynolds stress tensor and turbulent heat flux. We give special attention to dynamic and thermal time-scale modeling to reproduce the near-wall behavior of Reynolds stress tensor and turbulent heat flux. Then, we introduce the four equation model for the turbulent variables to close the system with near-wall boundary conditions. In Section 3, we illustrate the computational domain and numerical settings for the numerical simulations of the plane channel and backward-facing step configurations. We compare, then, the obtained results with the reference DNS data. Finally, conclusive remarks and future perspectives are provided in Section 4.

2. Mathematical Model

2.1. Dynamic Turbulence Modeling

The governing equations for the velocity field are written as:

$$\begin{array}{cc}\hfill & \frac{\partial u_i}{\partial x_i}=0,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & \frac{D{u}_i}{Dt}=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\nu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\overline{u_i^{\prime }{u}_j^{\prime }}\right],\hfill \end{array}$$

(2)

where $u_i$ and $u_i^{\prime }$ are the mean and fluctuating velocity components respectively, p is the mean pressure, and $\nu$ and $\rho$ are the kinematic viscosity and density. The overbar operator $\overline{\left(\right)}$ implies the mean value of a quantity, while $D/Dt$ is the substantial derivative $D/Dt=\partial /\partial t+u_j\partial /\partial x_j$. The Reynolds stress tensor $\overline{u_i^{\prime }{u}_j^{\prime }}$ is the averaged product of velocity fluctuations.

The Explicit Algebraic Stress Model (EASM) derives from the full transport equation for the Reynolds stress tensor under the hypothesis of local equilibrium between production and dissipation [23]. The hypothesis of local equilibrium is not coherent with many flow configurations. However, these are the usual hypotheses needed to close the model. The explicit expression here presented is only valid for two-dimensional flows in an inertial frame. The Reynolds stress tensor $\overline{u_i^{\prime }{u}_j^{\prime }}$ can be expressed as follows [24]:

$$\overline{u_i^{\prime }{u}_j^{\prime }}=\frac{2}{3}k{\delta }_{ij}-\frac{2{\nu }_t}{f_R}{S}_{ij}-\frac{4C_{D}k{f}_{\tau }}{f_R}\left(S_{ik}{\mathsf{\Omega }}_{kj}-{\mathsf{\Omega }}_{ik}{S}_{kj}-S_{ik}{S}_{kj}+\frac{1}{3}{S}^2{\delta }_{ij}\right),$$

(3)

where $S_{ij}$ is the strain-tensor and ${\mathsf{\Omega }}_{ij}$ the vorticity tensor:

$$\begin{array}{c}\hfill S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right),\end{array}$$

(4)

$$\begin{array}{c}\hfill {\mathsf{\Omega }}_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}-\frac{\partial u_j}{\partial x_i}\right),\end{array}$$

(5)

while $S^2=S_{ij}{S}_{ij}$ and ${\mathsf{\Omega }}^2={\mathsf{\Omega }}_{ij}{\mathsf{\Omega }}_{ij}$. The eddy viscosity ${\nu }_t$ is given by ${\nu }_t=C_{\mu }{f}_{\mu }k{\tau }_{lu}$, where $C_{\mu }$ and $f_{\mu }$ are the model constant and function, ${\tau }_{lu}$ denotes the characteristic dynamical time scale, k is the turbulent kinetic energy, and $\epsilon$ is its dissipation rate:

$$\begin{array}{c}\hfill k=\frac{1}{2}\overline{u_i^{\prime }{u}_i^{\prime }},\epsilon =\nu \overline{\frac{\partial u_i^{\prime }}{\partial x_k}\frac{\partial u_i^{\prime }}{\partial x_k}}.\end{array}$$

(6)

The value assigned to the constant $C_{\mu }$ is the standard value 0.09. The function $f_{\mu }$ and the time scale ${\tau }_{lu}$ play a key role in the eddy viscosity description and turbulence modeling. The modeling of $f_{\mu }$ is performed using the non-dimensional wall-distance $R_d$ [25], defined as $R_d=\upsilon y_d/\nu =y_d/\eta$, where $\upsilon$ is the Kolmogorov velocity scale $\upsilon ={\left(\nu \epsilon \right)}^{\frac{1}{4}}$, $\eta$ is the Kolmogorov length scale $\eta ={({\nu }^3/\epsilon )}^{\frac{1}{4}}$, and $y_d$ is the wall distance at a point, i.e., the distance between that point and nearest point on the wall surfaces. We define the function $f_{\mu }$ as follows:

$$f_{\mu }=1-exp\left[-{\left(\frac{R_d}{26}\right)}^2\right].$$

(7)

The characteristic time scale ${\tau }_{lu}$ is generally expressed with the scale of energy-containing eddies, ${\tau }_u=k/\epsilon$. However, in the proximity of the wall, there are always dissipation eddies that have to be taken into account. The effect of the dissipation eddies is added to the contribution of energy-containing eddies in the expression of the time scale:

$${\tau }_{lu}={\tau }_u(1+\frac{B_{\mu }}{R_t^{\frac{3}{4}}}{f}_{\eta }).$$

(8)

The model constant $B_{\mu }$ and the model function $f_{\eta }$ represent the effectiveness of dissipative motions and the limitation of their influence. We set $B_{\mu }=35$ and $f_{\eta }=exp(-{(R_t/30)}^{\frac{3}{4}})$, where $R_t$ is the turbulent Reynolds number $R_t=k^2/\nu \epsilon$. The function $f_R$ is given by:

$$f_R=1+\frac{22}{3}{\left(C_D{\tau }_{R_0}\right)}^2{\mathsf{\Omega }}^2+\frac{2}{3}{\left(C_D{\tau }_{R_0}\right)}^2({\mathsf{\Omega }}^2-S^2)f_B,$$

(9)

where $C_D=0.8$ and the quantity ${\tau }_{R_0}$ is the characteristic time scale of turbulence defined as ${\tau }_{R_0}={\nu }_t/k$. The model function $f_B$ is introduced to guarantee non-negative turbulent intensities when $S^2\gg {\mathsf{\Omega }}^2$. The formulation proposed is the following [23]:

$$f_B=1+C_{\eta }{\left(C_D{\tau }_{R_0}\right)}^2({\mathsf{\Omega }}^2-S^2).$$

(10)

The function $f_{\tau }$ reproduces the wall-limiting behavior and anisotropy of the Reynolds normal stress components near the wall and it is defined as:

$$f_{\tau }={\tau }_{R_0}^2+{\tau }_{R_W}^2,$$

(11)

where ${\tau }_{R_W}$ is the wall reflection time scale, defined by the expression:

$${\tau }_{R_W}=\sqrt{\frac{f_R}{6C_D{f}_{S\mathsf{\Omega }}}}\left(1-\frac{3C_{v1}{f}_{v2}}{8}\right)f_{v1}^2,$$

(12)

where $f_{v2}=1-exp(-\sqrt{R_t}/100)$, $f_{v1}=exp(-R_{tm}^2/2025)$ and $C_{v1}=0.4$. In the model function $f_{v1}$ the modified Reynolds number $R_{tm}$ is defined as follows:

$$R_{tm}=\frac{130R_d{R}_t^{\frac{1}{4}}}{130R_t^{\frac{1}{4}}+R_d}.$$

(13)

The modeling of the function $f_{S\mathsf{\Omega }}$ is a crucial aspect. In [24], the following expression is suggested for $f_{S\mathsf{\Omega }}$:

$$f_{S\mathsf{\Omega }}=\frac{{\mathsf{\Omega }}^2}{2}+\frac{S^2}{3}-{\left[\left(\sqrt{\frac{S^2}{2}}-\sqrt{\frac{{\mathsf{\Omega }}^2}{2}}\right)f_w\left(1\right)\right]}^2,$$

(14)

with $f_w\left(1\right)=exp(-R_{tm}^2)$.

Once the model for $\overline{u_i^{\prime }{u}_j^{\prime }}$ and ${\nu }_t$ is chosen, it is necessary to compute the variables appearing in the model functions, in particular the turbulent kinetic energy k and the characteristic time scale ${\tau }_u$. We propose a logarithmic turbulence model K-$\mathsf{\Omega }$ which improves the stability of a standard k-$\omega$ model since the state variables are maintained always positive during the solution process [26]. The specific dissipation rate of turbulent kinetic energy $\omega$ and the logarithmic form of k and $\omega$ are defined as:

$$\omega =\frac{\epsilon }{C_{\mu }k},\mathsf{\Omega }=ln\left(\omega \right),K=ln\left(k\right).$$

(15)

The system of equations for the K-$\mathsf{\Omega }$ model is the following:

$$\begin{array}{cc}\hfill & \frac{DK}{Dt}=\frac{\partial }{\partial x_i}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{\partial K}{\partial x_i}\right]+\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{\partial K}{\partial x_i}\frac{\partial K}{\partial x_i}+\frac{P_k}{e^K}-C_{\mu }{e}^{\mathsf{\Omega }},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & \begin{array}{c}\hfill \frac{D\mathsf{\Omega }}{Dt}=\frac{\partial }{\partial x_i}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\omega }}\right)\frac{\partial \mathsf{\Omega }}{\partial x_i}\right]+\left(\nu +\frac{{\nu }_t}{{\sigma }_{\omega }}\right)\frac{\partial \mathsf{\Omega }}{\partial x_i}\frac{\partial \mathsf{\Omega }}{\partial x_i}+2\left(\nu +\frac{{\nu }_t}{{\sigma }_{\omega }}\right)\frac{\partial K}{\partial x_i}\frac{\partial \mathsf{\Omega }}{\partial x_i}+\\ \hfill +\frac{P_k}{e^K}(C_{\epsilon 1}-1)-C_{\mu }(C_{\epsilon 2}{f}_{\epsilon }-1)e^{\mathsf{\Omega }},\end{array}\hfill \end{array}$$

(17)

where $P_k=-\overline{u_i^{\prime }{u}_j^{\prime }}\partial u_i/\partial x_j$ is the production rate of turbulent kinetic energy. The turbulent diffusion terms in Equations (16) and (17) are modeled using the Simple Gradient Diffusion Hypothesis (SGDH).

The model function $f_{\epsilon }$ has been modified since [13,14,15,25], now we set [23]:

$$f_{\epsilon }=\left\{1-0.3exp\left[-{\left(\frac{R_t}{6.5}\right)}^2\right]\right\}\left\{1-exp\left[-{\left(\frac{R_d}{3.7}\right)}^2\right]\right\},$$

(18)

and model constants ${\sigma }_k={\sigma }_{\omega }=1.4,C_{\epsilon 1}=1.5,C_{\epsilon 2}=1.9$.

2.2. Thermal Turbulence Modeling

The governing equation for the thermal field can be written as:

$$\frac{DT}{Dt}=\frac{\partial }{\partial x_i}\left(\alpha \frac{\partial T}{\partial x_i}-\overline{u_i^{\prime }{T}^{\prime }}\right),$$

(19)

where T and $T^{\prime }$ are the mean and fluctuating temperature, $\alpha$ is the thermal diffusivity, and $\overline{u_i^{\prime }{T}^{\prime }}$ is the turbulent heat flux. The explicit Algebraic Heat Flux Model (EAHFM) derives from the transport equation for the turbulent heat flux in the local equilibrium state neglecting the diffusive term [27]. Under this hypothesis, the explicit algebraic expression for the turbulent heat flux can be written as:

$$\overline{u_i^{\prime }{T}^{\prime }}=-\frac{C_{t1}}{f_{RT}}{\tau }_m\overline{u_i^{\prime }{u}_j^{\prime }}\frac{\partial T}{\partial x_j}+\frac{C_{t1}}{f_{RT}}{\tau }_m^2[(C_{t2}-C_{t3})S_{ij}+(C_{t2}-C_{t3}){\mathsf{\Omega }}_{ij}]\overline{u_j^{\prime }{u}_k^{\prime }}\frac{\partial T}{\partial x_k},$$

(20)

where $C_{t1}=0.18$, $C_{t2}=0.18$, and $C_{t3}=0.02$ are model constants. In contrast with traditional models, the components of turbulent heat flux $\overline{u_i^{\prime }{T}^{\prime }}$ and the mean temperature gradient are not necessarily in alignment due to the effects of the mean-velocity gradient and the Reynolds stress tensor $\overline{u_i^{\prime }{u}_j^{\prime }}$. To predict the heat transfer in wall flows, the characteristic time scale ${\tau }_m$ plays a key role since the turbulent heat flux (20) derives from the local equilibrium hypothesis, which does not hold in the near-wall region. The characteristic time scale ${\tau }_m$ is defined as the harmonic average of the dynamical time scale ${\tau }_u=k/\epsilon$ and the thermal time scale ${\tau }_t=k_{\theta }/{\epsilon }_{\theta }$:

$${\tau }_m\propto \frac{1}{\frac{1}{{\tau }_u}+\frac{C_m}{{\tau }_t}}={\tau }_u\frac{R}{C_m+R},$$

(21)

where we have introduced the ratio between the two scales $R={\tau }_t/{\tau }_u$. The terms appearing in the expression of the thermal time scale are $k_{\theta }$, the temperature fluctuations variance, and ${\epsilon }_{\theta }$, its dissipation:

$$\begin{array}{c}\hfill k_{\theta }=\frac{1}{2}{\overline{T^{\prime }}}^2,{\epsilon }_{\theta }=\alpha \overline{\frac{\partial T^{\prime }}{\partial x_k}\frac{\partial T^{\prime }}{\partial x_k}}.\end{array}$$

(22)

The composite time scale defined by Equation (21) is the harmonic average of the velocity and temperature time scales. The shortest time scale among ${\tau }_u$ and ${\tau }_t$ is the most important for turbulent heat flux.

In the bulk region, ${\tau }_m$ is independent of the time ratio R and the turbulent diffusion is assumed to be dominated only by velocity fluctuations. In the bulk region therefore we assume ${\tau }_m\propto {\tau }_u/P{r}_{t,\infty }$, where $P{r}_{t,\infty }$ can be assumed constant and uniform or can be modeled, for example by means of Kays model $P{r}_t=0.85+0.7/Pr{\nu }_t$ [28]. We also should introduce a model function in the ${\tau }_m$ expression to account for the wall-proximity effects. For the near-wall region, the characteristic time scale is ${\tau }_m\propto \sqrt{2R}/Pr{R}_t^{\frac{3}{4}}$. The characteristic thermal time scale is then modeled as:

$${\tau }_m={\tau }_u{f}_{1t}\left(\frac{1}{P{r}_t}+\frac{2R}{R+C_{\gamma }}{f}_{2t}+1.3\frac{\sqrt{2R}}{Pr{R}_t^{\frac{3}{4}}}{f}_{3t}\right),$$

(23)

where $C_{\gamma }=0.25/P{r}^{\frac{1}{4}}$. The model function $f_{1t}$ accounts for wall-proximity effects, and we set as in previous works:

$$f_{1t}=\left[1-exp\left(-\frac{R_d}{14}\right)\right]\left[1-exp\left(-\frac{\sqrt{Pr}{R}_d}{14}\right)\right].$$

(24)

The blending functions $f_{2t}$ and $f_{3t}$ are defined as:

$$\begin{array}{c}\hfill f_{2t}=exp\left[-{\left(\frac{R_t}{500}\right)}^2\right],f_{3t}=exp\left[-{\left(\frac{R_t}{200}\right)}^2\right].\end{array}$$

(25)

In order to evaluate $k_{\theta }$ and ${\epsilon }_{\theta }$ appearing in the model functions, we propose a logarithmic $K_{\theta }$-${\mathsf{\Omega }}_{\theta }$ turbulence model, where $K_{\theta }$ and ${\mathsf{\Omega }}_{\theta }$ represent the logarithmic values of mean temperature fluctuations $k_{\theta }$ and its dissipation rate ${\omega }_{\theta }$, defined as ${\omega }_{\theta }={\epsilon }_{\theta }/C_{\mu }{k}_{\theta }$. The transport equations for the logarithmic quantities can be written as:

$$\begin{array}{cc}\hfill & \frac{D{K}_{\theta }}{Dt}=\frac{\partial }{\partial x_i}\left[\left(\alpha +\frac{{\alpha }_t}{{\sigma }_{k\theta }}\right)\frac{\partial K_{\theta }}{\partial x_i}\right]+\left(\alpha +\frac{{\alpha }_t}{{\sigma }_{k\theta }}\right)\frac{\partial K_{\theta }}{\partial x_i}\frac{\partial K_{\theta }}{\partial x_i}+\frac{P_{k_{\theta }}}{e^{K_{\theta }}}-C_{\mu }{e}^{{\mathsf{\Omega }}_{\theta }},\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & \begin{array}{c}\hfill \frac{D{\mathsf{\Omega }}_{\theta }}{Dt}=\frac{\partial }{\partial x_i}\left[\left(\alpha +\frac{{\alpha }_t}{{\sigma }_{{\omega }_{\theta }}}\right)\frac{\partial {\mathsf{\Omega }}_{\theta }}{\partial x_i}\right]+\left(\alpha +\frac{{\alpha }_t}{{\sigma }_{{\omega }_{\theta }}}\right)\frac{\partial {\mathsf{\Omega }}_{\theta }}{\partial x_i}\frac{\partial {\mathsf{\Omega }}_{\theta }}{\partial x_i}+\frac{P_{k_{\theta }}}{e^{K_{\theta }}}(C_{p1}-1)+\\ \hfill 2\left(\alpha +\frac{{\alpha }_t}{{\sigma }_{{\omega }_{\theta }}}\right)\frac{\partial K_{\theta }}{\partial x_i}\frac{\partial {\mathsf{\Omega }}_{\theta }}{\partial x_i}+C_{p2}\frac{P_k}{e^K}-(C_{d1}-1)C_{\mu }{e}^{{\mathsf{\Omega }}_{\theta }}-C_{d2}{C}_{\mu }{e}^{\mathsf{\Omega }},\end{array}\hfill \end{array}$$

(27)

where $P_{k_{\theta }}=-\overline{u_j^{\prime }{T}^{\prime }}\partial T/\partial x_j$ and $C_{d2}$ is the following model function:

$$C_{d2}=\left\{1.9\left[1-0.3exp\left(-\frac{R_t^2}{42.25}\right)\right]-1\right\}\left[1-exp\left(-\frac{R_d^2}{25}\right)\right],$$

(28)

while $C_{p1}=1.025$, $C_{p2}=0.9$, $C_{d1}=1.1$, and ${\sigma }_{k_{\theta }}={\sigma }_{{\omega }_{\theta }}=1.4$. The eddy thermal diffusivity appearing in the diffusive terms of (26) and (27) is simplified as the scalar quantity ${\alpha }_t=C_{\theta }k{\tau }_m$, with $C_{\theta }=0.1$.

2.3. Boundary Conditions

In this subsection, we describe the boundary conditions that can be imposed on the state variables of the turbulence models. When a near-wall approach with no wall functions is used, the boundary conditions can be computed by a near-wall Taylor series expansion for the turbulence variables. For the description of the boundary conditions, we refer to the case of a plane channel where x is the wall distance, y is the streamwise coordinate, and z is the spanwise one. Moreover, u, v, and w are respectively the wall-normal, streamwise, and spanwise velocity components. In

Table 1. Near-wall Taylor expansion for the components of the mean velocity $u_i$, fluctuating velocity $u_i^{\prime }$, mean temperature T, and fluctuating temperature $T^{\prime }$.

Mean Components	Fluctuating Components
$u=A_2{x}^2+A_3{x}^3$	$u^{\prime }=a_2{x}^2+a_3{x}^3$
$v=B_{1}x+B_2{x}^2+B_3{x}^3$	$v^{\prime }=b_{1}x+b_2{x}^2+b_3{x}^3$
$w=C_{1}x+C_2{x}^2+C_3{x}^3$	$w^{\prime }=c_{1}x+c_2{x}^2+c_3{x}^3$
$T=D_0+D_{1}x+D_2{x}^2$	$T^{\prime }=d_0+d_{1}x+d_2{x}^2$

, we report the expansion for the mean and fluctuating velocity and temperature. Following the definitions (6) and (15), we obtain the following dynamical turbulence variable expansions:

$$\begin{array}{cc}\hfill & k_w\approx \frac{1}{2}(b_1^2+c_1^2)x^2=\frac{1}{2}\xi x^2,{\epsilon }_w\approx \nu (b_1^2+c_1^2)=\nu \xi ,{\omega }_w\approx \frac{2\nu }{C_{\mu }{x}^2},\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & K_w\approx ln\left(\frac{1}{2}\xi x^2\right),{\mathsf{\Omega }}_w\approx ln\left(\frac{2\nu }{C_{\mu }{x}^2}\right),\hfill \end{array}$$

(30)

where the lower-script w means the near-wall behavior. Since the value of $\xi$ depends on the components of fluctuating velocity and it is not known a priori, we transform the Dirichlet conditions (29) and (30) into Neumann conditions. By taking the derivative of k in the wall-normal direction x, we obtain ${\partial k/\partial x|}_w=\xi x=2k_w/x$, and considering the same derivative for the logarithmic variable ${\partial K/\partial x|}_w=2/x$, then for both variables it is possible to impose Neumann boundary conditions. The dissipation of turbulent kinetic energy $\epsilon$ has a constant near-wall value that can be determined from $k_w$, in particular $\epsilon =2k_w/x^2$, thus an exact Dirichlet boundary condition cannot be imposed on $\epsilon$, but the value of $\xi$ is iteratively calculated from the value of k on the wall and this can lead to convergence issues [14]. This aspect does not affect $\omega$ and $\mathsf{\Omega }$ since their values on the walls depend only on the kinematic viscosity of the fluid $\nu$, on the wall distance x, and on the model constant $C_{\mu }$. For these variables we can then impose the exact Dirichlet conditions (29) and (30).

The issue of boundary conditions on fluctuating thermal variables is still an open question [13,17,29,30]. For the energy equation, we can impose a constant wall temperature or a uniform wall heat flux. In the case of a constant wall temperature boundary condition, namely $T=T_w$, the condition must be fulfilled by both temperature and fluctuating values, so that $T^{\prime }=0$ along the wall. If a constant heat flux is applied, then temperature fluctuations can be considered null or not. If we assume that the temperature fluctuations are null (MX boundary conditions) we have $d_0=0$ and $T_w^{\prime }\approx d_{1}x$, then from definitions we obtain the following expressions:

$$\begin{array}{cc}\hfill & k_{\theta }\approx \frac{1}{2}{d}_1^2{x}^2,{\epsilon }_{\theta }\approx \alpha d_1^2,{\omega }_{\theta }\approx \frac{2\alpha }{C_{\mu }{x}^2},\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & K_{\theta }\approx ln\left(\frac{1}{2}{d}_1^2{x}^2\right),{\mathsf{\Omega }}_{\theta }\approx ln\left(\frac{2\alpha }{C_{\mu }{x}^2}\right).\hfill \end{array}$$

(32)

As in the dynamical turbulence case, the expressions of $k_{\theta }$, ${\epsilon }_{\theta }$, and $K_{\theta }$ depend on $d_1$ which is not known “a priori”. We can then reformulate (31) and (32) considering the derivative of $k_{\theta }$ and $K_{\theta }$ in the wall normal direction x as:

$$\begin{array}{cc}\hfill & \frac{\partial k_{\theta }}{\partial x}{|}_w=d_1^{2}x=\frac{2k_{\theta w}}{x},\frac{\partial K_{\theta }}{\partial x}{|}_w=\frac{2}{x},\hfill \end{array}$$

(33)

and impose Neumann boundary conditions. The quantity ${\epsilon }_{\theta w}$ is affected by the same issue of ${\epsilon }_w$ since we cannot impose an exact Dirichlet condition on this variable but only apply a Dirichlet boundary condition with a value $\alpha d_1^2$ that changes iteration by iteration. For ${\omega }_{\theta }$ and ${\mathsf{\Omega }}_{\theta }$, we can impose the exact Dirichlet conditions (31) and (32).

3. Numerical Results and Validation of the A4P Model

In this section, we aim to validate the anisotropic four-parameter turbulence model that we have illustrated in Section 2. The A4P model has been validated, simulating two different benchmark configurations. First, we have considered a plane channel configuration at different friction Reynolds numbers $R{e}_{\tau }$ for different low Prandtl numbers. The results obtained with the A4P are compared with the reference DNS data [17,19,20]. Then, a more complex configuration as a backward-facing step is considered and the numerical results are compared with DNS data [22] in a forced convection regime. The simulations have been performed using the in-house finite element multigrid code FEMuS developed at the University of Bologna [31]. The code is based on a C++ main program that handles several external open-source libraries, such as MPI and PETSc libraries. FEMuS contains solvers for Reynolds-Averaged Navier–Stokes and energy equations, four-parameter turbulence model, and explicit algebraic models for Reynolds stresses and turbulent heat flux.

3.1. Plane Channel Geometry

The plane channel flow configuration has been investigated by different authors in the framework of turbulence modeling [32,33]. This two-dimensional geometry has enabled the development of DNS reference databases, where the channel flow is characterized by different properties such as $R{e}_{\tau }$, $Pr$, and $Gr$ in the case of buoyancy, but also for the temperature boundary condition at the wall. Therefore, our simulations with RANS modeling are computed with fixed parameters in order to refer to specific DNS data. In this work the following $R{e}_{\tau }$ = 180, 395, 640, 1020 are considered for a $Pr$ number equal to 0.025 [17], meanwhile for the $Pr$ number equal to 0.01 we consider the friction Reynolds number cases for $R{e}_{\tau }$ = 180, 395, 590 [19], and 1000 [20]. In Figure 1, the computational domain is shown with the specific dimensions and the axis orientation. The problem consists of two plates located at the distance $D=2L_x=0.0605$ m characterized by the presence of a constant heat flux q equal to $3.6\times {10}^5$ W/m${}^2$. The other directions have infinite dimensions. The configuration computed is a fully developed turbulent channel flow with the presence of periodic boundary conditions on the inlet and outlet. From Figure 1, we can refer to the inlet section for the ${\mathsf{\Gamma }}_i$ section, the outlet for ${\mathsf{\Gamma }}_o$, the heated wall for ${\mathsf{\Gamma }}_w$, and the center of the channel for ${\mathsf{\Gamma }}_{sym}$. A pressure drop force $F\propto R{e}_{\tau }$, drives the channel flow, where $R{e}_{\tau }$ is defined as $u_{\tau }{L}_x/\nu$. Following this definition, we introduce the friction velocity $u_{\tau }=\sqrt{{\tau }_w/\rho }$ using the wall shear stress ${\tau }_w=\mu \frac{\partial v}{\partial x}{|}_w$, where v is the flow velocity parallel to the wall and x is the distance to the wall. In

Table 2. Physical properties employed for the numerical simulations.

Property	Symbol	Value	Units
Viscosity	$\nu$	0.001844	Pa s
Density	$\rho$	10,340	kg/m${}^3$
Thermal conductivity	$\lambda$	10.72–26.88	W/(mK)
Specific heat	$C_p$	145.75	J/(kgK)

the fluid properties are reported, with the first value of the thermal conductivity $\lambda$ referring to the case of $Pr=0.025$, and the second one for the case of $Pr=0.01$. Concerning the thermal field, the temperature has been redefined with the introduction of the variable $\theta =T-T_{w0}-L_y\Delta T_b$, where $\Delta T_b$ is the normal temperature difference, $T_{w0}$ is a constant value on ${\mathsf{\Gamma }}_w$, and $L_y$ is the axial length of the computational domain. From a computational point of view, the simulations have been computed only for one-half of the channel flow configuration, due to the symmetry of the problem. A mesh refinement near the wall ${\mathsf{\Gamma }}_w$ is performed to have the first mesh point in the viscous laminar region, in particular $x^{+}<1$, where the non-dimensional distance from the wall $x^{+}$ is defined as $x{u}_{\tau }/\nu$.

In the following graphs, all the variables are plotted in their non-dimensional forms. The variables are normalized using wall units, i.e., the friction velocity $u_{\tau }$, friction temperature $T_{\tau }$, and kinematic viscosity $\nu$. The friction velocity is used to normalize the velocity $v^{+}=v/u_{\tau }$ and the components of the Reynolds stress tensor ${\overline{u_i^{\prime }{u}_j^{\prime }}}^{+}=\overline{u_i^{\prime }{u}_j^{\prime }}/u_{\tau }^2$. The friction temperature $T_{\tau }=q/u_{\tau }\rho C_p$ is used to normalize the temperature ${\theta }^{+}=\theta /T_{\tau }$ and the components of the turbulent heat flux ${\overline{u_i^{\prime }{T}^{\prime }}}^{+}=\overline{u_i^{\prime }{T}^{\prime }}/u_{\tau }{T}_{\tau }$.

In Figure 2, the non-dimensional streamwise velocity $v^{+}$ is plotted against the non-dimensional distance from the wall $x^{+}$, for different $R{e}_{\tau }$ numbers, i.e., 180, 395, 640, and 1000, corresponding to the following Reynolds numbers, i.e., $Re\approx$ 5700, 14,100, 24,400, and 41,400. The comparison with the DNS data shows a good matching both in the linear and logarithmic region for different cases of friction Reynolds number.

In Figure 3, the components of the Reynolds stress tensor are shown. We can notice from Figure 3a a good match between the non-dimensional turbulent shear stress ${\overline{u^{\prime }{v}^{\prime }}}^{+}$ and the corresponding DNS data. On the right, in Figure 3b, the non-dimensional turbulent streamwise normal stress ${\overline{v^{\prime }{v}^{\prime }}}^{+}$ is shown. We can observe an overall good agreement with DNS for every $R{e}_{\tau }$, even though the quantities are slightly underestimated. The dimensionless wall-normal normal stress ${\overline{u^{\prime }{u}^{\prime }}}^{+}$ is shown in Figure 3c,d for $R{e}_{\tau }=180,395$ and $R{e}_{\tau }=640,1020$ respectively. In all these cases, the simulations results show an overall overestimation of the peak of the wall-normal normal stress component with respect to the DNS data.

To better understand these results, we underline that the plane channel configuration is a typical example of shear flow. Due to the symmetry of the problem, the only component of the Reynolds stress tensor that affects the mean velocity field is the turbulent shear stress, i.e., the off-diagonal component $\overline{u^{\prime }{v}^{\prime }}$. For this reason, the mean velocity is correctly estimated even though the diagonal components of the Reynolds stress tensor present some discrepancies from DNS data. These stresses are only used for the estimation of the turbulent heat flux components, according to Equation (20).

Concerning the thermal fields, the simulations results are shown for both $Pr$ numbers to validate the new anisotropic four-parameter turbulence model. The non-dimensional temperature profiles are shown in Figure 4a,b for the case $Pr=0.025$ and $Pr=0.01$ respectively. For both cases, the temperature field is in good agreement with the reference DNS data, that are available for all the test cases except for the case $R{e}_{\tau }=1020$ and $Pr=0.025$.

In Figure 5, we report the non-dimensional wall-normal turbulent heat flux ${\overline{u^{\prime }T}}^{+}$ profiles for different $R{e}_{\tau }$ and $Pr$ numbers. For shortness, only two cases of different $R{e}_{\tau }$ are shown for each $Pr$ number. In Figure 5a,b, the plots show the results for $Pr=0.025$, while in Figure 5c,d for $Pr=0.01$. We also report the effective wall-normal heat flux $q_{eff,x}^{+}$ which is defined as the sum of the molecular heat flux and the turbulent heat flux in the wall-normal direction, i.e.,

$$q_{eff,x}^{+}={\left(\alpha \frac{\partial \theta }{\partial x}\right)}^{+}-{\overline{u^{\prime }{T}^{\prime }}}^{+}.$$

(34)

We can notice some discrepancies for the ${\overline{u^{\prime }{T}^{\prime }}}^{+}$ component from DNS data in the near-wall region. However, in this region, the thermal conductivity contribution is dominant and the turbulent heat flux is almost negligible, then the total heat flux $q_{eff,x}^{+}$ is almost equal to the molecular heat flux. Thus, the bad prediction of ${\overline{u^{\prime }{T}^{\prime }}}^{+}$ in this region does not affect the total heat flux and, consequently, the mean temperature field. Moreover, increasing the distance from the wall, the turbulent heat flux contribution becomes higher but the molecular heat flux is still dominant. For this reason, the slight discrepancies in ${\overline{u^{\prime }{T}^{\prime }}}^{+}$ near the center of the channel for $Pr=0.025$ do not compromise the mean temperature estimation.

In Figure 6, the non-dimensional streamwise component of the turbulent heat flux ${\overline{v^{\prime }{T}^{\prime }}}^{+}$ is shown for $Pr=0.025$ and $Pr=0.01$. We report in Figure 6a the streamwise heat flux profiles for $Pr=0.025$ and $R{e}_{\tau }=180,640$, while in Figure 6b the profiles are shown for $Pr=0.01$ and $R{e}_{\tau }=395,1000$. The streamwise component of the turbulent heat flux is underestimated for all the cases with respect to the DNS reference data. We can then conclude that the anisotropic four-parameter turbulence model is not able to properly predict this component in the plane channel configuration. However, this bad prediction does not affect the mean temperature estimation. Indeed, due to the symmetry of the plane channel configuration, the mean temperature field is only affected by the wall-normal component of molecular and turbulent heat flux. Moreover, we underline that with the isotropic version of the presented model, the wall-normal turbulent heat flux component would be identically zero in this configuration.

3.2. Backward Facing Step Geometry

In this subsection, we report the results obtained for the simulation of a turbulent flow of liquid sodium over a vertical backward-facing step. This type of flow has been extensively studied in the literature. In [21,22,34,35], DNS simulations with different Reynolds numbers have been performed, in forced and/or mixed convection regimes. In [36], a comparison between the solutions of a RANS system of equations closed with various turbulence models is proposed for the forced convection case, showing that four-equation turbulence models, coupled with nonlinear expressions for the Reynolds stress tensor, allow for improving the predictions of the turbulent heat flux. In [11], an anisotropic three-equation turbulence model has been proposed and applied to this configuration in forced and mixed convection regimes showing a promising potential for the prediction of the turbulent heat flux. In [16], RANS simulations have been performed with an isotropic four-parameter turbulence model (4P) in forced and mixed convection regimes considering a linear expression for the Reynolds stress tensor and the turbulent heat flux. Results are promising in both regimes but the adoption of the anisotropic formulation that we are proposing could improve the prediction of the turbulent heat flux components.

The computational domain reproduces the reference DNS domain [22] and a representative sketch is reported in Figure 7. The inlet section length is $L_{in}$, the step height is h, the domain width is W, and the downstream channel height is E. The expansion ratio is $E_r=E/(E-h)$. The geometrical parameters of the simulated domain are reported in

Table 3. Backward facing step: Geometrical parameters of the simulated domain.

${\mathit{L}}_{\mathbf{in}}/\mathit{h}$	${\mathit{L}}_{\mathit{h}}/\mathit{h}$	W	${\mathit{E}}_{\mathit{r}}$	${\mathbf{Re}}_{\mathit{b}}$	$\mathbf{Pr}$
2	20	0	1.5	9610	0.0088

. The system of equations is solved with the finite element code FEMuS [31]. The mesh consists of 27,820 cells with 107,411 biquadratic nodes. Mesh cells are clustered near the corner step and reattachment zone. Mesh refinement is performed near wall boundaries to obtain a non-dimensional wall distance $x^{+}<1$ on the first mesh point near wall boundaries.

In terms of boundary conditions, a fully-developed inflow condition has been set for the velocity field and turbulent variables corresponding to $R{e}_{\tau }=h{u}_{\tau }/\nu \approx 300$ and $R{e}_b=2h{U}_b/\nu =9610$. The Reynolds number is calculated with respect to the inlet channel width $2h$ and the bulk inlet velocity $U_b$. For the temperature, a uniform value is set, i.e., $T_{ref}=423.15K$. The same temperature is used as the reference value for the evaluation of liquid sodium properties obtaining $Pr=0.0088$. At the outlet section, an outflow boundary condition is imposed on the velocity field and for all the other variables homogeneous Neumann conditions are set. All the remaining boundaries have been treated as adiabatic no-slip walls, except for the wall behind the step where a uniform heat flux $\dot{q}$ is imposed.

Numerical simulations have been performed for the forced convection calculations using the anisotropic four-parameter turbulence model (A4P) and the results are compared in the next subsections with DNS data and with the results from numerical simulations using the standard isotropic four-parameter model (4P) [4,13,14,15,16,37].

3.3. Dynamical Fields

In this subsection, the results obtained with the anisotropic four-parameter model (A4P) for flow fields are compared with DNS data and with simulations results using the isotropic four-parameter model (4P).

In Figure 8a, the contours of the non-dimensional streamwise velocity component $v^{+}=v/U_b$ and the streamlines of the velocity field are shown. The typical flow features for a backward-facing step configuration are observed, i.e., the flow separation taking place behind the step, the reattachment of the flow, and the formation of two main vortexes behind the step: A bigger one rotating in the clockwise direction and a smaller one rotating in the opposite direction.

Non-dimensional profiles of velocity taken on channel cross-section planes are reported for several streamwise coordinate $y/h$ values for the A4P results, reported with solid lines, and the 4P results, shown in dashed lines. The streamwise positions included in these plots correspond to the locations where DNS data are available [22].

The streamwise $v^{+}$ and wall-normal velocity component $u^{+}=u/U_b$ are reported respectively in Figure 9a,b. The velocity field prediction with the anisotropic four-parameter model is in good agreement with DNS results, while the isotropic model shows a slight deviation of the wall-normal velocity $u^{+}$ from DNS data at the non-dimensional height $y/h=3$. The typical features of the flow over a backward-facing step, i.e., the separation and reattachment behind the step, can be observed from the plot taken at $y/h=3$ where $v^{+}$ assumes negative values close to the heated wall.

In Figure 9c, the non-dimensional shear stress component ${\overline{u^{\prime }{v}^{\prime }}}^{+}=\overline{u^{\prime }{v}^{\prime }}/U_b^2$ is reported in comparison with DNS data. The turbulent shear stress is mainly present in the layer behind the step at $x/h\approx 0$. The prediction of the shear stress is in good agreement with reference data for the anisotropic model results. For the isotropic four-parameter model results, we have computed the shear stress as $\overline{u^{\prime }{v}^{\prime }}={\nu }_t(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})$. The prediction obtained with the isotropic model (4P) is in good agreement with DNS data, even though there are some discrepancies in the plot taken at $y/h=3$.

The wall-normal normal stress ${\overline{u^{\prime }{u}^{\prime }}}^{+}=\overline{u^{\prime }{u}^{\prime }}/U_b^2$ and the streamwise normal stress ${\overline{v^{\prime }{v}^{\prime }}}^{+}=\overline{v^{\prime }{v}^{\prime }}/U_b^2$ are reported respectively in Figure 10a,b. In Figure 10c, the turbulent kinetic energy $k^{+}=k/U_b^2$ is shown. These turbulent fields present a general good agreement with DNS data even though there are some discrepancies in the region behind the step. For the isotropic four-parameter model simulations, we have computed $\overline{u^{\prime }{u}^{\prime }}=\frac{2}{3}k+2{\nu }_t\frac{\partial u}{\partial x}$ and $\overline{v^{\prime }{v}^{\prime }}=\frac{2}{3}k+2{\nu }_t\frac{\partial v}{\partial y}$, according to Boussinesq approximation.

In Figure 11, the skin friction coefficient $c_f$ along the heated wall is reported. The skin friction profile is subjected to a double change of sign, denoting the presence of two reattachment points. The skin friction coefficient assumes negative values in the recirculation zone, which is composed of a large clockwise rotating vortex. Directly behind the step, the principal recirculating vortex causes a secondary vortex rotating in the opposite direction. The position of the first reattachment point is approximately $y/h\approx 1.26$ for the anisotropic four-parameter model and $y/h\approx 1.91$ for DNS data. The second reattachment point is located approximatively at $y/h\approx 6.52$. The DNS data give this point at $y/h\approx 7.01$, while Kasagi [38] gives this point at $y/h\approx 6.51$ through measurements.

3.4. Thermal Fields

In this subsection we propose a comparison for thermal fields between the results obtained with the anisotropic four-parameter model (A4P) and DNS data. We also report the results of simulations performed using the isotropic four-parameter model (4P). In Figure 8b, contours of the non-dimensional temperature $T^{+}=\frac{T-T_{ref}}{\Delta T}$ are reported for the simulation with the anisotropic four-parameter model (A4P). The hot fluid is located in the corner between the step and heated wall. A strong temperature increase is observed moving from the insulated wall towards the heated wall. The highest wall temperature is located in the recirculation zone and reaches a maximum closely behind the step due to the reduced heat transfer due to the backward flow.

Non-dimensional profiles of the mean temperature $T^{+}=T/\Delta T$ are reported in Figure 12a for different values of streamwise coordinate $y/h$. The temperature difference $\Delta T$ is defined using the applied heat flux $\dot{q}$ setting $\Delta T=\dot{q}h/\lambda$, where $\lambda$ is the liquid sodium thermal conductivity calculated for $T=T_{ref}$. The discrepancies with DNS results are limited to the plot taken at $y/h=0$ where the temperature is slightly overestimated. With the isotropic model, the major discrepancies with DNS values are found on the plots taken at $y/h=0$ and $y/h=3$ where an over- and under-prediction of $T^{+}$ is respectively obtained.

The turbulent heat flux components along wall-normal ${\overline{u^{\prime }{T}^{\prime }}}^{+}=\overline{u^{\prime }{T}^{\prime }}/(U_b\Delta T)$ and streamwise ${\overline{v^{\prime }{T}^{\prime }}}^{+}=\overline{v^{\prime }{T}^{\prime }}/(U_b\Delta T)$ directions are reported respectively in Figure 12b,c. The anisotropic model allows improving the prediction of the streamwise component which is completely underestimated with the isotropic model. The isotropic model assumes a unique scalar thermal diffusivity ${\alpha }_t$ for both turbulent heat flux components, i.e., $\overline{u^{\prime }{T}^{\prime }}={\alpha }_t\frac{\partial T}{\partial x}$ and $\overline{v^{\prime }{T}^{\prime }}={\alpha }_t\frac{\partial T}{\partial y}$. However, the mean temperature gradient along the streamwise direction is small and for this reason the streamwise component is totally underestimated. With the proposed anisotropic model, the wall-normal component shows a better agreement with DNS data and the streamwise component results are only slightly underestimated.

Nusselt number profiles along the heated wall are shown in Figure 13. The Nusselt number is computed as $Nu=\dot{q}h/(T-T_{ref})\lambda$. When the Nusselt value is around 1 then the heat transfer is mostly diffusive due to the low Prandtl number of the liquid metal. Inside the recirculation zone, we have $Nu<1$, then the heat transfer is prevented by the recirculating flow. As one can see in Figure 13, in the recirculation zone, the Nusselt number is slightly overestimated with the isotropic model, while the anisotropic model is in good agreement with DNS data in all the regions.

4. Conclusions

In this work, we presented a new anisotropic four-parameter turbulence model (A4P) that derives from the four-parameter turbulence model (4P) widely studied in [4,13,14,15,16] within the framework of heat transfer modeling for low-Prandtl number fluids. An Explicit Algebraic Stress Model (EASM) and an Explicit Algebraic Heat Flux Model (EAHMF) was proposed for the closure of the Reynolds stresses and turbulent heat flux instead of first-order closure relations used in the isotropic version of the model. Special attention is given to the modeling of the dynamical and thermal time scales to overcome the local equilibrium hypothesis typical of algebraic models. The closure of the model and estimation of the time scales were performed with four transport equations for the logarithmic variables K-$\mathsf{\Omega }$-$K_{\theta }$-${\mathsf{\Omega }}_{\theta }$ and suitable near-wall boundary conditions were presented. For the validation of the new anisotropic four-parameter turbulence model, we considered two different configurations, i.e., forced convection in a plane channel and over a backward-facing step, considering different $Re$ and low-$Pr$ numbers. The simulation results were compared with the available DNS data and for the backward-facing step configuration, and a comparison was also proposed with the isotropic version of the proposed model. The prediction of the velocity and temperature fields was in good agreement with DNS reference data for all the considered configurations. For the forced convection over the backward-facing step configuration, we could observe a general improvement in the prediction of dynamical and thermal fields with respect to the isotropic version of the model, above all in the estimation of the turbulent heat flux components. It can be concluded that the anisotropic four-parameter turbulence model represents a promising approach towards the accurate prediction of both turbulent momentum and heat flux for low-Prandtl number fluids. Further simulations, including complex geometries and buoyancy effects will be presented in future works.