According to the survey, the turbulence Prandtl number model is mainly divided into two types: the model related to the wall distance and the model unrelated to the wall distance, which are called the local turbulence Prandtl number model and the global turbulence Prandtl number model in this paper, respectively. The following representative models are selected for subsequent evaluation and analysis.

The global turbulent Prandtl number model mainly includes the following:

Aoki (1963): Pr t − 1 ⁡ = 0.014 R e 0.45 ⁡ Pr 0.2 1 − ⁡ exp − 1 0.014 R e 0.45 ⁡ Pr 0.2 . (1)

Reynolds (1974): Pr t ⁡ = 1 + 100 P e − 1 / 2 1 1 + 120 R e − 1 / 2 − 0.15 . (2)

Jischa and Rieke (1979): Pr t ⁡ = 0.9 + 182.4 P r R e 0.888 . (3)

Cheng and Tak (2006): Pr t ⁡ = 4.12 f o r   P e ⩽ 1000 0.01 P e 0.018 P e 0.8 − 7.0 − A 1.25 f o r   1000 < P e ⩽ 6000 = 4.5   f o r   P e ⩽ 1000 5.4 − 9 × 10 − 4 P e   f o r   1000 < P e ⩽ 2000 3.6   f o r   P e > 2000 . (4)

The local turbulent Prandtl number model mainly includes the following:

Kays (1994) found that in the logarithmic region, Pr _t seems to be mainly a function of Pe _t . When Pe _t is large, Pr _t generally approaches a constant of 0.85. When Pe _t decreases, Pr _t increases with uncertainty. For liquid metals with low Pr, the Pr _t value is relatively large, and the following calculation formula is proposed. The results are in good agreement with DNS and the experimental data at the center of the pipe. Pr t ⁡ = 0.85 + 0.7 P e t , (5) P e t = ε M υ Pr ⁡ = ε M κ . (6)

Weigand et al. (1997) obtained an extended relationship based on the Kays turbulent Pr _t model, and the calculated Nu can be in good agreement with the experimental data on a fully developed and thermally developed pipeline flow. This formula considers the spatial distribution of turbulent Pr _t , that is, the influence of the wall distance, and can express the turbulent Pr _t situation in the entire velocity boundary layer. It is more reasonable to deal with the near wall. 1 Pr t = 1 2 ⁡ Pr t ∞ + C P e t 1 Pr t ∞ − C P e t 2 1 − ⁡ exp − 1 C P e t Pr t ∞ , (7) Pr t ∞ ⁡ = 0.85 + 100 Pr ⁡ R e 0.888 , C = 0.3 . (8)

Li et al. (2019) put forward the three-zone theory, deduced the lead–bismuth Pr _t model, and established the Pr _t zoning theoretical model suitable for lead–bismuth. The first zone (Y+ < 50) adopts Weigand’s Pr _t model, and the second zone (50 ≤ Y+ ≤ 500) establishes a new Pr _t model by combining the renormalization group analysis method and turbulent heat transfer equation method and further using the physical property Pr correction. The third zone (Y+ > 500) adopts a Pr _t value of 1, as shown in Table 1.

This study will evaluate the existing turbulent Prandtl number models mentioned above. A set of experimental data with a relatively wide range of operating conditions and representative geometric conditions was selected as a reference, which was an experiment conducted by Crye et al. (2002) in the United States to measure the heat transfer coefficient of liquid Hg in a narrow rectangular channel. Tables 2 and 3 are the specific channel dimensions and operating conditions.

A three-dimensional, narrow rectangular channel model is established. The geometric size and boundary conditions fully refer to the above experimental conditions. In previous studies (Zhang et al., 2022b), after comparison and analysis, it was found that the turbulence model that selected the standard k-epsilon model, combined with the enhanced wall treatment model, was more reasonable. After grid independence analysis, the grid division method is as follows: 300 nodes are evenly divided along the flow direction of the channel. A total of 120 segments are divided in the long side direction of the cross section, and densification is performed near the wall. The thickness of the first layer of the grid is 0.0005 mm, with a growth rate of 1.2. The narrow edge of the cross section is divided into a total of 40 segments, with a first-layer grid thickness of 0.0005 mm and a growth rate of 1.1. There are approximately 920 thousand grids in total.

The numerical calculation is based on Fluent 18.0, which is a typical commercial CFD software application. For the solution method, SIMPLE is adopted as the pressure–velocity coupling scheme. The spatial discretization of pressure adopts the second-order format and that of momentum adopts the bounded central differencing format. Using the steady-state calculation method, the residual convergence standard is set to 10⁻⁶. The physical properties of the working fluid Hg in the flow channel are calculated by the following polynomial equations embedded in the software application, which are transformed from nonlinear equations (Jaeger et al., 2015).

Density, kg/m³, is expressed as follows: ρ = − 9.61751 × 10 − 7 T 3 + 9.89515 × 10 − 4 T 2 − 2.84542 T + 14320.20376 . (9)

Specific heat capacity, J/kg·K, is expressed as follows: C p = − 6.11564 × 10 − 9 T 3 + 8.03019 × 10 − 5 T 2 − 0.08138 T + 156.68834 . (10)

Thermal conductivity, W/m·K, is expressed as follows: λ = − 6.378 × 10 − 8 T 2 + 0.01363 T + 4.4584 . (11)

Dynamic viscosity, kg/(m·s), is expressed as follows: μ = 7.98384 × 10 − 15 T 4 − 2.54132 × 10 − 11 T 3 + 3.03574 × 10 − 8 T 2 − 1.65123 × 10 − 5 T + 0.00431 . (12)

“T" represents the fluid temperature in K. The calculation method of the Nusselt number refers to the calculation method in the experimental data processing process. The temperature at a certain point on the wall is selected as the wall temperature, and the average temperature of the cross section is obtained based on energy conservation of the fluid temperature. In this study, the position in the middle of the heating section along the axial direction was selected as the calculation position.

First, the existing global turbulence Prandtl number model was evaluated based on experimental data, as shown in Table 4.

Two sets of typical operating conditions with different Peclet numbers were selected, and the Nusselt numbers calculated using different models were compared. Among them, Nu-cal represents the calculated value of CFD; Nu-exp represents the Nusselt number calculated by the Lubarski and Kaufman relationship, which is highly consistent with the experimental data (Crye et al., 2002); and Nu-J represents the empirical value of the Nusselt number calculated by the relationship of Jaeger et al. (2015). The results show that all four models have certain errors in predicting Nu. The smaller the Pe, the greater the error. Among them, the calculation results of Cheng et al.'s model are closest to the empirical values, but the error is still greater than 20%.

Then, the applicability of three typical local turbulence Prandtl number models was evaluated. The three models were implanted into FLUENT software using the UDF method, and the Nusselt numbers were calculated and compared with the empirical values. The results are shown in Table 5:

The results show that the above three local turbulent Prandtl number models all have significant errors in calculating the heat transfer process of liquid metal flow in narrow rectangular channels, which should be because the fitting data source of the above semi-empirical model formula is mostly circular tubes or rod bundles.

V. Karman divides the flow field into three parts: a viscous bottom layer (Y+<5), a transition layer (5<Y+<30), and a turbulent wake layer (Y+>30). In the intermediate transition layer, it is assumed that the molecular and turbulent physical quantities are the same. Li et al. (2019) drew inspiration from the idea of V. Karman’s analogy and proposed a partition theory for low Pr fluids. According to the heat transfer mechanism, the entire turbulent region was divided into three parts: (1) in the velocity boundary layer, due to the high thermal conductivity of liquid metal molecules (about 20 times that of water) and the relatively small flow rate in the velocity boundary layer, the molecular thermal conductivity dominates heat transfer compared to the eddy current diffusion rate. (2) The transition zone between the velocity boundary layer and the heat flux boundary layer, where the flow velocity has reached the mainstream velocity, enhances eddy heat transfer, and the molecular viscosity is much smaller than the eddy viscosity. Due to the good thermal conductivity and high thermal diffusion coefficient of liquid metals, it can be assumed that the molecular and turbulent physical quantities in terms of heat transfer are the same. (3) The fully mainstream region outside the heat flow boundary layer, with fully developed velocity and temperature.

Drawing on the above analytical approach, a typical numerical example is used to analyze the temperature field and flow field in the narrow gap of a narrow rectangular channel in this paper. Figure 1 shows the velocity and temperature distribution along the narrow slot direction during the turbulent heat transfer process of a narrow rectangular liquid metal. It can be seen that the velocity boundary layer is very thin, and most areas have reached the mainstream velocity. Due to the good thermal conductivity of the liquid metal, the temperature boundary layer is very thick, and the gap size of the narrow rectangular channel is very small, resulting in the temperature boundary layers of the two heating surfaces overlapping. Therefore, the entire flow region only includes the near-wall region and the transition region.

It can be observed that the transition region accounts for the main proportion. Therefore, this article focuses on the turbulent Prandtl number model in the turbulent region. For the near-wall region, based on existing research experience (Li et al., 2019), the model improved by Weigand et al. (1997) is selected. For the transition region that accounts for the majority, the basic form of the model is determined by referring to Jischa’s method for solving turbulent transport equations (Jischa and Rieke, 1979). Based on the experimental data, the unknown coefficients of the semi-empirical model are iteratively solved and finally compared with other experiments to verify the accuracy of the model. The following will introduce the process of establishing the semi-empirical model.

The definition formula of the turbulent Prandtl number is shown in Equation 13. Pr t ⁡ = u ′ v ′ ¯ ∂ T ∂ y v ′ T ′ ¯ ∂ u ∂ y . (13)

In the equation, u represents the velocity in the x direction (perpendicular to the heating wall), and v represents the velocity in the y direction (along the direction of flow). In order to obtain the Pr t , u ′ v ′ ¯ and v ′ T ′ ¯ need to be derived by solving the transport equation. The following will introduce the derivation steps of Jischa and Rieke (1979).

The work is based on the work of Prandtl, who modeled the transport equation for the turbulent kinetic energy as follows: k = 1 2 u ′ 2 ¯ + v ′ 2 ¯ . (14)

For Prandtl’s original modeling assumptions:

Turbulent shear stress − u ′ v ′ ¯ = ε m ∂ u ∂ y = C τ k L ∂ u ∂ y . (15)

Production is expressed as follows: − u ′ v ′ ¯ ∂ u ∂ y = C τ k L ∂ u ∂ y 2 . (16) Dissipation is expressed as follows: ε = C D k 3 / 2 L . (17) Diffusion is expressed as follows: − ∂ ∂ y … . . = C ∂ ∂ y k L ∂ k ∂ y . (18)

Here, L is the turbulence length scale. In the following process of derivation, only Eq 15 will be used. The others are given to lead over to the similar modeling of the transport equation.

To obtain v ′ T ′ ¯ , one multiplies the equation for the instantaneous value of temperature (T + T′) with the velocity fluctuation v j ′ and adds it to the x j component of the Navier–Stokes equations multiplied by the temperature fluctuation T ′ . j and k represent the x and y directions, respectively: v k ∂ T ′ v j ′ ¯ ∂ x k = − v j ′ v k ′ ¯ ∂ T ∂ x k − v k ′ T ′ ¯ ∂ v j ∂ x k + 1 ρ C p v j ′ v k ′ ¯ ∂ p ∂ x k ⏟ P r o d u c t i o n + 1 ρ C p ρ ε ′ v j ′ ¯ + λ v j ′ ∂ 2 T ′ ∂ x k 2 ¯ + C p η T ′ ∂ 2 v j ′ ∂ x k 2 ¯ ⏟ D i s s i p a t i v e   t e r m s + 1 ρ C p − C p T ′ ∂ p ′ ∂ x j ¯ + v j ′ ∂ p ′ ∂ t ¯ + v k v j ′ ∂ p ′ ∂ x k ¯ + v j ′ v k ′ ∂ p ′ ∂ x k ¯ ⏟ R e d i s t r i b u t i v e   t e r m s − ∂ ∂ x k T ′ v j ′ v k ′ ¯ ⏟ D i f f u s i o n . (19)

Eq 19 is expanded and simplified to obtain the following Eq 20. u ∂ ∂ x T ′ v ′ ¯ + v ∂ ∂ y T ′ v ′ ¯ ⏟ C o n v e c t i o n + v ′ 2 ∂ T ¯ ∂ y − u ′ v ′ ¯ ρ C p d p d x ⏟ P r o d u c t i o n − 1 C p v ′ ε ¯ + ν Pr ⁡ + 1 Pr ∂ v ′ ∂ x ∂ T ′ ∂ x ¯ ⏟ D i s s i p a t i o n − p ′ ρ ∂ T ′ ∂ y + … ⏟ R e d i s t r i b u t i o n + ∂ ∂ y T ′ v ′ 2 ¯ − … ⏟ D i f f u s i o n = 0 . (20)

For the research object between the velocity boundary layer and the temperature boundary layer, the velocity of the flow has been fully developed; therefore, the convection term is zero in Eq 20.

In analogy to the modeling of Eq 16, neglecting the pressure gradient production, the production term is transformed as follows: v ′ 2 ∂ T ∂ y = a k ∂ T ∂ y 2 . (21)

The diffusion term can be neglected in the boundary layer. Usually, all dissipation terms in the transport equations can be neglected to expect the turbulent dissipation ε . This is because when the turbulence Reynolds number is high, regardless of isotropic or anisotropic turbulence, the dissipation correlation is zero, which means the redistribution term limits the growth of the Reynolds heat flux. So, the dissipation term and redistribution terms both can be modeled using the same assumption in analogy to Eq 17, as follows: Dissipation and  / or redistribution = K D k L T ′ v ′ ¯ . (22)

Therefore, the following can be obtained: a k d T d y + K D k L T ′ v ′ ¯ = 0 , (23) − T ′ v ′ ¯ = a K D k L d T d y , (24) Pr t ⁡ = C K D a . (25)

The above derivation ignores the second term of the dissipation, which depends on the molecular Prandtl number. However, when Pr is particularly small, the impact of this term cannot be ignored.

Need to make a new assumption Dissipation and  / or redistribution = K D + K p Pr ⁡ + 1 Pr k L T ′ v ′ ¯ . (26)

Then, a k d T d y + K D + K p Pr ⁡ + 1 Pr k L T ′ v ′ ¯ = 0 , (27) Pr t ⁡ = C a K D + K p Pr ⁡ + 1 Pr = A + B Pr ⁡ + 1 Pr = E + B Pr . (28)

Here, “E" is a constant number. When the Pr is large enough, it can be considered that the fluid has the same characteristics as conventional fluids such water so that “E" can be equal to 1.0. The coefficient “B" seems to be a function of the Reynolds number under the effect of the diffusion term, which was not considered in the previous derivation. After taking it into account, the above equation can be converted into the following equation: Pr t ⁡ = 1.0 + F P e t . (28a)

Finally, based on CFD and using the experimental data as a benchmark for an iterative solution, the value of the constant number "F" is selected to obtain the turbulent Prandtl number model for the transition region, which is expressed as follows: Pr t ⁡ = 1.0 + 8.0 P e t . (29)

The applicable range of the model for Pe: 1,200–3,000. Figure 2 shows the comparison of the Nu between that calculated by the CFD implanted the new Pr t model with the experimental data by Crye et al. (2002). Figure 3 shows the error between the Nu calculated by the CFD implanted the new Pr t model with that calculated by the Lubarski and Kaufman relationship (Lubarsky and Kaufman, 1956). The error is within the range of ±20% and is acceptable.

In order to explore the applicability of the newly proposed turbulent Prandtl number model, the obtained model was validated with the experimental data of liquid NaK in parallel narrow rectangular channels measured by Duchatelle and Vautrey (1964) from France. Tables 6 and 7 show the characteristics of the test section and the operating conditions.

The working fluid NaK consists of 78% potassium and 22% sodium. The physical properties can be calculated by the following equations.

Density, kg/m³, is expressed as follows: ρ = 953.8505 − 0.27 T . (30)

Specific heat capacity, J/kg·K, is expressed as follows: C p = 1238.86485 − 1.5126 T + 0.002 T 2 . (31)

Thermal conductivity, W/m·K, is expressed as follows: λ = 141.16476 − 0.86103 T + 0.00156 T 2 . (32)

Dynamic viscosity, kg/(m·s), is expressed as follows: μ = 0.00267 − 9.38492 × 10 − 6 T + 9.4 × 10 − 9 T 2 . (33)

When Pe is greater than 200, the experimental data are highly consistent with an empirical relationship: N u = 5.58 + 0.0003417 P e 1.29 . (34)

Figure 4 shows the comparison of Nu calculated by the proposed local Pr t model with experimental data, empirical values, and the previously proposed global model (Zhang et al., 2022b). The results show that the simulation results are in good agreement with the experimental data within approximately 20% error, and the calculation results of the local Pr t model proposed in this article are closer to the experimental values than those of the previously proposed global Pr t model. This indicates that the local Pr t model has better applicability to different runner sizes and different liquid metal types.