Theoretical Investigation on the Fully Developed Turbulent Heat Transfer Characteristics of Liquid Sodium

Liquid sodium has been used as the working medium for the high-temperature heat pipe and the generation IV sodium cooled fast reactor due to its extremely high conductivity. The heat transfer characteristics of sodium in a circular pipe is one of the most essential focuses in engineering applications. In this paper, a model to predict the heat transfer coeffcient of fully developed sodium flowing in tube was developed based on universal velocity, turbulent eddy diffusivity, and the linear law inside the thermal boundary layer. The Kays correlation for turbulent Prandtl number was used to predict the turbulent Prandtl in the bulk flow with y + larger than 60. This model was validated by experiment data of Hg, NaK, and sodium, showing superior accuracy than existing models. Besides, the dependence of the accuracy on the model parameters was also analyzed, demonstrating the universal applicability of the current model.


INTRODUCTION
Liquid metals, especially liquid sodium, have been widely used in high-temperature heat transfer process in industry and power applications, such as sodium cooled fast reactors in France, China, and Russia, and the high temperature heat pipe used in space facilities. To ensure the safety of the facilities, the heat transfer coefficient should be accurately estimated in the thermal design process of industrial applications. However, the heat transfer characteristics of liquid metal are quite different from traditional fluids. The molecular Prandtl numbers of traditional fluids vary from 1 to 10, while the Prandtl numbers of liquid metals are always <0.1, even 0.01. Furthermore, as for liquid sodium, the Prandtl number is in the range of 0.01-0.001 due to extremely large thermal conductivities. The difference in Prandtl number changes the heat transfer characteristics in the near wall and bulk regions. The thermal boundary layers of liquid metal are much thicker than the momentum boundary layer, resulting in the failure of correlations for heat transfer of traditional fluids in the applications of liquid metals. Besides, the discrepancy in the thermal and momentum boundary layers leads to the failure of Reynolds analogy. Thus, the traditional computational fluids dynamics treatment on the energy equation based on Reynolds analogy and the unit turbulent Prandtl number is not validated in liquid metals.
In the past few decades, researchers have devoted themselves to theoretical and experimental investigations. Lyon (1951) proposed the first empirical correlation for turbulent heat transfer of liquid sodium. Then Lubarsky and Kaufman (1955) corrected the Lyon correlation based on experiment data and published a new correlation for channels with specified geometry. Schleisiek (1970) carried out a series of experiments on liquid sodium heat transfer in a 9-mm inner diameter pipe and derived a formula for Nusselt number based on 24 sets of data. Shi et al. (1981) also gave correlations for liquid sodium, but these correlations have limited applicable range of Peclet number <3,000. Except for the correlations aiming at liquid sodium heat transfer, several correlations were also proposed based on data from other liquid metals, such as mercury (Hg) (Chang and Akins, 1972), sodium-potassium (NaK) eutectic (Sleicher et al., 1973), and lead-bismuth eutectic (LBE) (Ibragimov et al., 1961). In addition, Dwyer (1966) proposed a heat transfer correlation for the universal applications of liquid metal. However, these correlations predict the Nusselt number with large differences. Furthermore, the results show a noticeable divergence even for the same working fluid, especially for sodium and NaK eutectic, which means the investigation on the heat transfer characteristics of sodium has still not matured.
In this work, a theoretical model was proposed to predict the turbulent heat transfer coefficient of liquid sodium. This model employed the correlations for universal velocity and eddy diffusivity to predict the distributions of velocity and turbulence. The linear law was used to predict the nondimensional temperature in the near wall linear sublayer inside the thermal boundary layer. Beyond this sublayer, the turbulent Prandtl number was predicted using the Kays correlation (Kays, 1994), which is only related to the local turbulent parameters. The governing equations for non-dimensional temperature and heat flux can be solved after being discretized on the radial direction. The model was validated by experiment data, showing superior accuracy compared with other available models. The uncertainty from the model parameters was also investigated based on the parametric analysis.

Governing Equations for Heat Transfer
The governing equation for energy conservation in a fully developed pipe can be simplified into the following format: Due to the wall heat flux, the temperature in the pipe cannot be described by a fully developed field with zero gradient in the flow direction since the temperature rises along the pipe. However, the temperature increases linearly with the flow direction in a pipe when the flow and thermal parameters are fully developed. Thus, the temperature T(x, r) in the pipe can be divided into two parts plus a fixed wall temperature at the inlet, that is, where T w,x=0 is the wall temperature at the channel inlet; θ (r) is the temperature along the radial direction, which is independent from the axial location, that is, fully developed just as the velocity field. ∂T ∂x x is the temperature increase in the pipe, which is linear with the flow distance or heated length, and can be calculated based on the energy balance, i.e., After substituting Equations (2) and (3) into Equation 1, we can obtain: Let R = r r w and q + = q q w . Then Equation 4 can be transformed into: According to Fourier's law on heat transfer, the radial heat flux at any radial location can be written as Given Pr = c p µ k l and Pr t = ε m ε q , Equation 6 can be changed into: Define the following non-dimensional parameters: where u τ = √ τ w /ρ and θ τ = q w /ρc p u τ . Then, Equations (5) and (8) can be written as: The boundary conditions for Equations (10) and (11) are listed as follows: Until now, we obtained the governing equations and boundary conditions for the dimensionless temperature and heat flux along the radial direction. These equations can be solved numerically by the finite difference method with limited discrete grid along the radius of the pipe. The total length of the radius can be divided into n − 1 identical intervals with n points, including the points at the pipe center (i = 0) and on the wall (i = n). Thus, the interval is: and the radial location referred from the center of the tube is: The discrete format for Equations (10) and (11) on the n grid in the radial direction of the pipe can be written as: and where the mean universal velocity u + m is calculated by: Theoretically, Equations (18) and (19) can be solved iteratively under the following conditions: Auxiliary Models and Strategy to Solve the Governing Equations As mentioned earlier, the numerical solutions for governing equations (Equations 10, 11) for energy balance in the radial direction could be obtained by solving Equations (18) and (19) with the boundary conditions of Equations (21) and (22). However, the radial profiles for the velocity, turbulent eddy diffusivity, and turbulent Prandtl number are still unknown and should be given to solve the above governing equations.

Universal Velocity and Eddy Diffusivity
As for the fully developed turbulent flow in the circular tube, various non-dimensional universal velocity profiles have been proposed in the last century, such as Prandtl's linear law for the sublayer (Johnson, 2016) and Taylor's logarithmic law for the logarithmic regime beyond the buffer layer (Taylor, 1916). Most of these correlations were only validated for the boundary layer, except the universal velocity correlation proposed by Reichardt (1957), that is, where the non-dimensional radial R locates in the range 0 ≤ R ≤ 1. The non-dimensional distance to the wall is calculated by Reichardt (1957) also proposed an empirical correlation for the eddy diffusivity for the momentum transfer caused by the stochastic movement of fluid from the wall to the center of the duct, which can be expressed in terms of the fraction of eddy diffusivity and kinematic viscosity, that is, where κ = 0.4.

Turbulent Prandtl Number
As for the liquid metal with the molecular Prandtl number much <1, the velocity boundary layer is much thinner than the thermal boundary layer, which means the Reynolds analogy with the turbulent Prandtl number near to unit is not validated under these conditions. Researchers proposed several correlations for the turbulent Prandtl number (Pr t ) with a different molecular Prandtl number for the prediction of heat transfer characteristics. These correlations are only related to the Reynolds number, the molecular Prandtl number, the Peclet number, or the maximum eddy diffusivity in the tube, which means the calculated turbulent Prandtl number is kept as a constant along the radial direction. However, the essential laminar law for the non-dimensional temperature (θ + = Pr · y + ) in the laminar regime inside the thermal boundary layer cannot be satisfied if Equation 19 is directly solved with a constant turbulent Prandtl number. This discrepancy is caused by the improper assumption of a uniform turbulent Prandtl number in the radial direction. As a matter of fact, the turbulent Prandtl number is highly associated to y + , especially in the near wall region with y + less than about 100 (Kawamura et al., 1999;Duponcheel et al., 2014). Thus, the turbulent Prandtl number relevant to the local y + should be specified for an exact prediction of the temperature profile in the radial direction. However, as can be noted from the limited DNS and LES simulations (Kawamura et al., 1998(Kawamura et al., , 1999Tiselj and Cizelj, 2012;Duponcheel et al., 2014), the turbulent Prandtl number in the near wall region is highly related to the local turbulent Peclet number non-linearly (Kawamura et al., 1999;Tiselj and Cizelj, 2012;Duponcheel et al., 2014), and there is no universal correlation for the localized turbulent Prandtl number in the near wall region so far (Grötzbach, 2013; Roelofs et al., 2015). As for the bulk flow region, Bricteux et al. (2012) assessed the DNS and LES data and pointed out that the turbulent Prandtl number should be about 2 for the bulk regime with y + larger than 200. Duponcheel et al. (2014) assessed the available models with results from large eddy simulation (LES) and found that the Kays correlation (Kays, 1994) is capable of predicting the turbulent Prandtl number in the bulk region, that is, where Pe t is the turbulent Peclet number defined as In Duponcheel et al. (2014), it is claimed that the Kays correlation (Kays, 1994) is validated for y+ < 100; nevertheless, it can be noted from the comparison between the Kays correlation and the LES data that the Kays correlation agrees well with the LES data for y+ < 65, which fortunately gives the opportunity to model the heat transfer of the low Prandtl number liquid metal by using the approach introduced in the next section.

Treatment on the Near Wall Heat Transfer Without Turbulent Prandtl Number in the Thermal Boundary Layer
As noted in the investigations on the universal velocity profile vs. the non-dimensional wall distance y + , linear, and logarithmic laws are universal for all kinds of flow conditions in the turbulent flow regime, independent from the Reynolds number of the flow. The non-dimensional temperature, θ + , which is defined in Equation 9, also satisfies the linear law and the logarithmic law in the boundary layer. Thus, we can use the conception of the wall function to calculate θ + in the thermal boundary layer. The near wall profile of θ + could be solved by the wall function based on y + , instead of the governing equation in the thermal boundary layer if the wall function for temperature is known priorly. However, due to the failure of the Reynolds analogy in the temperature distribution and the velocity distribution, the universal wall function for temperature in liquid metal with a low molecular Prandtl number is unknown. Besides, the thickness of the thermal boundary layer is much larger than that of the momentum boundary layer. The node in the velocity logarithmic regime may still be in the linear regime for energy transfer. Thus, the approach to set the first node in both the thermal and momentum logarithmic regimes, which is commonly used when dealing with water flow and heat transfer, is not validated. Fortunately, from the DNS and LES studies on the heat transfer of liquid metal with a low Prandtl number, it can be observed that the thickness of the linear regime in the thermal boundary layer increases with decreasing the molecular Prandtl number (Kader, 1981;Kawamura lab, 2008;Duponcheel et al., 2014). Besides, the thermal linear regime is beyond y + of 60 when the molecular Prandtl number is <0.01, as noted by Bricteux et al. (2012). As for the liquid sodium, the molecular Prandtl number is <0.01 for most of the conditions. That means, the thickness of the thermal linear regime should be larger than about 60. The non-dimensional temperature θ + could be predicted by the linear law in the thermal linear regime. Thus, the radial non-dimensional temperature distribution θ + could be calculated by: if y + is <y + c , and by Equation 19 if y + is larger than y + c , where y + c is the critical nondimensional wall distance set to be 60 in this work.

Analysis on the Non-dimensional Parameters Along the Radial Direction
After solving the governing equations for non-dimensional heat flux and temperature (Equations 18, 19, 30), the non-dimensional heat transfer coefficient can be obtained by: The  Figure 1, including Reynolds numbers of 5,000, 50,000, and 100,000 and molecular Prandtl numbers of 0.001 and 0.01 for sodium under different thermal conditions. The results from the model proposed by Taler (2016) are presented as references. For low Reynolds number conditions, the effects of the molecular Prandtl number are quite small, and our model predicts higher non-dimensional temperature than Taler's model, which means a lower heat transfer coefficient. Under the conditions with a moderate to high Reynolds number, our model will present a higher Nusselt number than Taler's model, and the difference between these two models increases with the molecular Prandtl number. It should be mentioned that the linear approach on the near wall temperature distribution affects the temperature distribution not only in the linear regime but also in the bulk regime, since our model will degenerate into Taler's model if the critical y + for the upper limit of the linear regime is set to be zero.

Model Validation With Experiment Data
In this section, the Nusselt number predicted by Equations (18) and (19) based on the turbulent Prandtl number given by Kays (1994) is compared with the experiment data to validate the accuracy of the models. First, the heat transfer data from Hg and NaK (Skupinski et al., 1965) were used as the benchmark database. Besides, the results predicted by correlations for liquid metal heat transfer are also given as references. As can be noted from the comparisons in Figure 2, the Nusselt number predicted by Taler correlation (Taler, 2016) and the current model with Pr t obtained from Kays correlation (Kays, 1994) shows superior agreement with the experiment data. The heat transfer coefficient data for liquid sodium (Fuchs, 1974) are also used to validate the model. As can be noted from Figure 3, the accuracy of the current model is much better than any other models when predicting the liquid metal with an extremely low molecular Prandtl number.
Parameter Sensitivity Analysis on the Critical y + for Linear Sublayer As mentioned above, the upper limit y + of the linear regime in the thermal boundary layer is set to be 60 in the current model. However, the validation region of the Kays correlation (Kays, 1994) is beyond y + of about 65. Besides, the upper limit of the linear regime may vary with the molecular Prandtl number. Thus, it is necessary to analyze the effects of critical y + (denoted as y + crit ) on the model accuracy. In this section, y + crit varies from 50 to 70, and the results are analyzed to measure the potential uncertainties introduced by the selection of y + crit . The results are shown in Figure 4, which demonstrates negligible effects of y + crit on the predicted Nusselt number, especially for flow with a high Reynolds number.

CONCLUSIONS
In this paper, a new model for turbulent heat transfer of liquid sodium in a circular pipe was proposed. This model was developed based on the universal profiles for velocity and eddy diffusivity, generally applicable for turbulent flow in pipes. The Kays correlation was employed to predict the turbulent Prandtl number in the bulk flow region with non-dimensional wall distance y + larger than 60. In the near wall region with y + <60, the linear law was used instead of the turbulent Prandtl number to predict the temperature distribution. The profiles of non-dimensional temperature and heat flux can be obtained by solving the governing equations for temperature and heat flux. The selection of y + equaling to 60 as the upper limit of the linear regime has little effects on the heat transfer characteristics, which is proved by a sensitivity study related to the influences of critical y + on the predicted Nusselt number. The new model developed in this paper has superior accuracy when compared with other available correlations by the validation with the experiment data.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary files.

AUTHOR CONTRIBUTIONS
RZ contributed significantly to the analysis and manuscript preparation. ZhenyW performed the data analysis and wrote the manuscript. ZhenhW checked the English writing thoroughly. TR and JS approved the final version.