Abstract

The prediction of flow and heat transfer characteristics of liquid sodium with CFD technology is of significant importance for the design and safety analysis of sodium-cooled fast reactor. The accuracies and uncertainties of the CFD models should be evaluated to improve the confidence of the numerical results. In this work, the uncertainties from the turbulent model, boundary conditions, and physical properties for the flow and heat transfer of liquid sodium were evaluated against the experimental data. The results of uncertainty quantization show that the maximum uncertainties of the Nusselt number and friction coefficient occurred in the transition zone from the inlet to the fully developed region in the circular tube, while they occurred near the reattachment point in the backward-facing step. Furthermore, in backward-facing step flow, the maximum uncertainty of temperature migrated from the heating wall to the geometric center of the channel, while the maximum uncertainty of velocity occurred near the vortex zone. The results of sensitivity analysis illustrate that the Nusselt number was negatively correlated with the thermal conductivity and turbulent Prandtl number, while the friction coefficient was positively correlated with the density and Von Karman constant. This work can be a reference to evaluate the accuracy of the standard k-ε model in predicting the flow and heat transfer characteristics of liquid sodium.

1. Introduction

Liquid sodium is widely used as a coolant in nuclear reactor engineering, aerospace, chemical engineering, and industrial waste heat utilization, due to the high thermal conductivity and low kinematic viscosity. The investigations on flow and heat transfer characteristics of liquid sodium began from 1940s [1–6]. The conclusion reached was that the flow characteristics of liquid sodium are similar to that of water [7–9]. However, Jenkins [10] indicated that the heat transfer characteristics of liquid sodium and water, such as the turbulent Prandtl number, were different. Wu et al. [11] made an overview of the flow and heat transfer characteristics of liquid sodium in different channel geometries including circular tube, annular tube, and rod bundle. However, conducting experimental research is difficult and expensive in some specific flow regions, such as narrow channels, helical wire spacer [12, 13], and sudden enlargement pipe. To predict the flow and heat transfer characteristics of liquid sodium more efficiently, numerical simulation method has been adopted. Meanwhile with the development of turbulence theory, a variety of models have been proposed for various flow and heat transfer conditions. The standard k-ε model has been improved in previous studies [14, 15] and widely used in commercial CFD codes, which is efficient in prediction of the flow and heat transfer characteristics of liquid sodium when it flows in tubes or complex channels with high Reynolds number in nuclear engineering normally. However, the prediction based on standard k-ε model may not be accurate, since this model is developed from conventional fluids such as water and air and widely used in commercial CFD codes. However, this model is developed based on conventional fluids such as water and air.

In the calculation process of turbulent heat transfer by standard k-ε, the method to solve the energy equation is similar to that dealing with the momentum equation. There are usually two methods [16]. One is using the turbulence Prandtl number to close energy equation. Kays [17] and Weigand et al. [18] proposed turbulent Prandtl number model and corresponding extended model. Taler [19] adjusted the turbulent Prandtl number of different models, which made the calculated values in agreement with the experimental data better. Ge et al. [20] assessed the applicability of different turbulent Prandtl number models in the bundle flow and recommended models of Kays [17] and Aoki [4]. The other method is to use the four-parameter equation model to describe the turbulent energy propagation process. Da Vià and Manservisi [21] simulated the turbulent flow of liquid sodium in backward-facing step with the four-parameter turbulence model, and the calculated results were in good agreement with the DNS data. However, there is still no uniform standard for the value of turbulent Prandtl number or the coefficients in the four-parameter equation. The uncertainty of these coefficients not only leads to the lack of universality but also reduces the accuracy of models. Meanwhile, many uncertainties will also generate from input boundary conditions, model coefficients, grid division, and so on. The accuracy of the standard k-ε model still needs an evaluation.

Uncertainty and sensitivity analysis has become a common practice in CFD to evaluate the accuracy of algorithm [22]. This method can quantify the accuracy of the turbulence model and find the key parameters which have a greater impact on the results. Dunn et.al [23] analyzed the uncertainty of the standard k-ε model in predicting the flow characteristic of water in backward-facing step. They found that the largest uncertainty appears near the recirculating region and the reattachment point. Wang et al. [24] analyzed the uncertainty of laminar and turbulent aeroheating predictions for Mars entry. It showed that, in the leeside region, the key parameters, which produce significant uncertainties in laminar and turbulent cases, are evidently different. Qureshi and Chan [25] discussed the influence of vorticity-to-strain ratio on turbulence parameters and adjusted vorticity-to-strain ratio reduces uncertainties, which improved the performance of standard k-ε model. Rakhimov et al. [26] proposed an uncertainty quantification method for CFD validated for turbulent mixing experiments from GEMIX. The results showed that uncertainty in the turbulence model input parameters has a significant contribution to the overall uncertainty in the simulation results. So far, there has not been enough evaluation on the uncertainty of the standard k-ε model, especially in terms of predicting the heat transfer characteristics of liquid sodium.

In this work, in order to provide a reference for the optimization of k-ε model in predicting turbulent flow and heat transfer characteristics of liquid sodium, uncertainty quantification and sensitivity analysis were adopted. The CFD software FLUENT was used to simulate the turbulent flow of sodium in the circular tube and the backward-facing step, and the uncertainty was imported by the Latin hypercube sampling (LHS) method [27]. The aim of this work is to obtain the rules of uncertainties in the standard k-ε model. Moreover, several key factors which influence the flow and heat transfer characteristics of liquid sodium were identified.

2. Mathematical and Physical Models

The k-turbulence model can be used to describe turbulence at high Reynolds number, which is composed of the transport equations of turbulent kinetic energy k and turbulent dissipation rate ε:

Source, diffusion, and dissipation of turbulent kinetic energy and turbulent kinetic energy dissipation are contained in the k-ε equations. To close the equation, the model usually needs to specify five coefficients, and the best combination is [28–30]

Besides, turbulent Prandtl number Pr_t will be introduced, which reflects the turbulent heat transport process. For liquid metals, the value of Pr_t is different from that for conventional fluids.

The principle for determining the coefficients is that the prediction by the k-ε model should be consistent with the results of DNS or experiment [31].(1)Determination of C_μ: In simple steady and parallel shear turbulence, the generation and dissipation of turbulent kinetic energy k are in equilibrium; that is, . According to vortex viscosity model: , we have In the linear layer and logarithmic layer,[31], so .(2)Determination of C_ε2: The equation of k and ε in uniform turbulence can be simplified as follows:

 The attenuation of uniform turbulence conforms to exponential law: [32], and can be obtained from the equation of k. The above formulas can be simplified to The value of n was measured in the grille turbulence experiment [33]; . However, it was adjusted later [34]:(3)Determination of C_ε1 and σ_ε: In the logarithmic region, the average velocity u obeys . Therefore, the expression of Reynolds stress is The expression of dissipation rate ε is On the other hand, the source and dissipation of turbulent kinetic energy are equal: Multiplying (8) and (10), In the logarithmic layer, Reynolds stress and turbulent kinetic energy k are approximately constant. The gradient of turbulent kinetic energy dissipation ε along the flow direction is approximately equal to zero, and the source term is in equilibrium with the dissipative term. Therefore, the equation of turbulent kinetic energy dissipation ε transport is An equation is still needed to determine C_ε1 and σ_ε. Some experimental results of uniform shear flow can be adopted to obtain their relation [31]. In the uniform shear turbulent field, the spatial derivative of the turbulence statistic is zero, so its standard k-ε equation can be simplified as follows: P is the turbulent energy . Experiments and DNS have confirmed that k/ε in uniform shear turbulence tends to be constant gradually [31]; that is, . To sum up, the following formula will be obtained: In the early studies, the measured value of P/ε in the uniform shear turbulence experiment was about 1.4, which was later revised to 2.1 [34]. The constant κ can be calculated by logarithmic law: . Usually, κ = 0.41, which was obtained by DNS in turbulent channel flow from Kim et al. [35]. Combined with (13), (15), and κ = 0.41, the relation of C_ε1 and σ_ε is(4)Determination of σ_k: Because the turbulent kinetic energy diffusion ε occurs in the heterogeneous turbulent field, it is difficult to determine the turbulent kinetic energy diffusion coefficient σ_k. Assuming that the transfer of turbulent kinetic energy k and momentum proceed in the same mechanism, In conclusion, the coefficients of k-ε are obtained as shown in Table 1. This set of coefficients was first proposed by [28] and is known as the standard k-ε constants. They are determined by three measured parameters (P/ε, C_ε2, and κ) and two assumed parameters ( and σ_k), which are dependent on theories and experiment. The model coefficients cannot be changed arbitrarily because of their relationship in the mathematical derivation.

In this work, the input uncertainties came from three sources: (1) boundary conditions, such as the input mainstream temp T, the mainstream velocity V, and heat flux Q; (2) physical property parameters, such as density ρ, dynamic viscosity μ_f, thermal conductivity λ, and thermal capacity C_p; (3) model coefficients, such as the five coefficients of standard k-ε and turbulent Prandtl number Pr_t.

Table 2 summarizes the input samples, which were independent from each other. Probability density functions of the sampled parameters were constructed by statistical information. The measured quantities, such as boundary conditions and physical properties, were subject to Gaussian distribution. The model coefficients also obeyed Gaussian distribution, which fluctuated around the nominal value. However, the turbulent Prandtl number was uniformly distributed within the range of 0.85–4.2 due to the absence of a unified nominal value.

In this work, interval estimation was carried out for samples conforming to Gaussian distribution [36]. The standard deviation s reflects the degree of fluctuation. The confidence interval of the expectation μ was μ ± 3s, and the uncertainty was defined as s/μ.

In terms of samples with non-Gaussian distribution, ±95% of the data was considered to be within the confidence interval.

Latin hypercube sampling (LHS) method was adopted to quantify the uncertainties, which can reduce the number of samples to achieve a reasonable random distribution compared with traditional Monte Carlo method [27]. Based on Wilks formula [37, 38], the sample size was 153, which satisfied the requirements of interval.

In sensitivity analysis, the Pearson coefficients [39] and Spearman coefficients [40] were used to evaluate the correlation. If the value of the correlation coefficient between two variables is greater than +0.3, they have a strong positive relationship, and vice versa [39, 40].

Uncertainty analysis software Dakota was used for sampling, and CFD software FLUENT was adopted for simulation. The combination of hydraulic diameter and turbulence intensity was input into the FLUENT solution, and the standard wall function was used to solve the near-wall region.

3. Results and Discussion

In this work, the rules of uncertainty propagation and key factors affecting flow and heat transfer characteristics in the standard k-ε model were analyzed and discussed. Liquid sodium flowed in the circular tube and backward-facing step. The output parameters are Nusselt number Nu, friction coefficient C_f, velocity, and temperature field. In terms of the backward-facing step, the reattachment point Xr was also considered.

3.1. Analysis of the Circular Tube

Liquid sodium flowed in a circular tube with uniform diameter, which was heated with constant heat flux. The inlet temperature was 673K. The Reynolds number was 93500. Given the enhanced heat transfer characteristics of the inlet segment, the modified parameters were taken into account [41]:

For Nu, the experimental correlation formula of [4] was used for comparison:

For C_f, the classical experimental correlation formula of Blasius [42] was used for comparison:

The results of uncertainty quantification are shown in Figures 1 and 2. In the fully developed region, the standard k-ε model had a divergence effect on Nu and C_f. The small input uncertainty would be amplified in the output. The maximum uncertainty of Nu and C_f occurred in the transition region from the inlet to the fully developed region (x/D ≈ 35–40).

Figure 1 

Axial uncertainty quantification of Nu with 6 s uncertainty.

Figure 2 

Axial uncertainty quantification of C_f with 6 s uncertainty.

The results of sensitivity analysis about Nu and C_f are shown in Figures 3–5 which show the scatter diagram. There was a strong positive relationship between Nu and the thermal conductivity λ. While the turbulent Prandtl number Pr_t was significantly strongly negatively correlated with Nu. For C_f, the density ρ and the Von Karman constant κ were significantly positively correlated with it.

Figure 3 

Correlation coefficients of Nu.

Figure 4 

Correlation coefficients of C_f.

Figure 5 

Scatter diagrams of Nu and C_f.

3.2. Analysis of the Adiabatic Backward-Facing Step

Boundary layer separation and reattachment in backward-facing step make the predicted results of standard k-ε model distorted. In the adiabatic condition, the configuration of step referred to Dunn [23] (Figure 6), and the boundary conditions were consistent with the experiment results of Kim [35]. The Reynolds number with outlet height H as the characteristic length is 132000 (Re_H = 132000).

Figure 6 

Backward-facing step geometry from Dunn [23].

The uncertainty of the normalized flow field is shown in Figure 7. The confidence interval of velocity completely included the experimental point when x/h > 5.33, which proved that the flow characteristics of liquid sodium are similar to that of water [7–9]. However, in the vortex region (x/h < 5.33), the mean value of velocity was a little larger than the experimental value. The experimental value was almost coincident with the lower limit of the confidence interval.

Figure 7 

Normalized flow field with 6 s uncertainty.

The standard k-ε model underestimates Xr by about 20% [43], and for this reason, many studies have adjusted the model coefficients to better predict the reattachment point [23]. Table 3 shows the results of K–S test for the Xr, which obeyed Gaussian distribution. Figure 8 shows the sensitivity analysis results of Xr. Figure 9 shows the scatter between the inputs and outputs which were significantly correlated. It illustrates that to improve Xr to approach the experiment data, the coefficient with positive relationship σ_ε needs to be improved, and the coefficient with negative relationship C_μ, C_ε1, C_ε2, and P/ε needs to be reduced.

Figure 8 

Correlation coefficients of Xr.

Figure 9 

Scatter diagrams of Xr.

The results of uncertainty quantification are shown in Figure 10. The C_f represented the friction coefficient on the side with the step. The maximum uncertainty occurred rear the region that was close to the reattachment point Xr (X = 0). At four times of Xr, the uncertainty gradually decreased to zero:

Figure 10 

Friction coefficient along the heating wall with 6 s uncertainty.

3.3. Analysis of the Heating Backward-Facing Step

The geometry of backward-facing step referred to Da Vià [21] is shown in Figure 11. The Reynolds number Re_h took the height h of the step as the characteristic length, and the value was 4805. The height of the channel was three times that of the step; that is, E/h = 3. The L_h is the heating surface with constant heat flux.

Figure 11 

Backward-facing step geometry from Da Vià [21].

The result of uncertainty quantification of normalized velocity and temperature field is shown in Figure 12. Before y/h≈6, the calculated mean values of temperature were more consistent with DNS data. However, after y/h≈6, the calculated mean values of the velocity field matched the DNS data better. On the whole, the uncertainty of temperature and velocity increased along the flow direction. The maximum uncertainty of the velocity always occurred about the vortex region, while the maximum uncertainty of the temperature field would migrate from near the wall to the geometric center of the flow channel.

Figure 12 

Normalized flow and temperature field with 6 s uncertainty.

The result of quantifying uncertainty of friction coefficient C_f and Nu along the flow direction is shown in Figures 13 and 14, respectively. The calculated results of C_f were not precise enough, especially before the vortex region. The uncertainty tends to be the maximum at the minimum value of C_f, which jumped from −10% to 20%. After y/h≈12, the mean value of C_f was about 30% higher than DNS data, and the confidence interval did not contain any DNS data points. The calculated results of Nu were a little bit more accurate, and there would be some DNS data points in the confidence interval after y/h≈6. There were similarities between C_f and Nu in the rule of uncertainty propagation. The uncertainties of C_f and Nu increased near the vortex region. As the flow fields became more stable, their uncertainty decreases. In the fully developed region, their uncertainties were less than those from the inlet.

Figure 13 

Friction coefficient along the bottom wall with 6 s uncertainty.

Figure 14 

Axial uncertainty quantification of Nu with 6 s uncertainty.

On the whole, according to the relationship from Table 1 and the conclusion of sensitivity analysis from Figures 12–14, appropriate adjustment on model coefficients should be considered to make the calculated results of C_f and Nu conform to DNS data better. The effective way to increase Nu was to reduce the turbulent Prandtl number reasonably, while the calculated value of C_f can be reduced by reasonably selecting a smaller κ.

4. Conclusion

In this work, the processes of turbulent flow and heat transfer of liquid sodium in circular tube and backward-facing step were simulated with the standard k-ε model. The results of uncertainty quantification and sensitivity analysis are as follows.

The uncertainties of Nu and C_f in the fully developed region are larger than those of the input. In the circular tube, their uncertainty first increases and then decreases along the flow direction. The maximum uncertainty occurs in the transition zone from the inlet to the fully developed region. In the backward-facing step, the maximum uncertainties of Nu and C_f occur near the reattachment point. The maximum uncertainty of the temperature field migrates from near the heating wall to the geometric center of the channel, while the maximum uncertainty of the velocity field occurs near the vortex region.

The results of sensitivity analysis show that, for Nusselt number, the thermal conductivity is significantly positively correlated with it, while the turbulent Prandtl number is significantly strongly negatively correlated. For friction coefficient, the density and the Von Karman constant are significantly positively correlated with it. For reattachment point Xr, the coefficient σ_ε is positively correlated with it, while the coefficient C_μ, C_ε1, C_ε2, and P/ε are negatively correlated. If the standard k-ε model is used in engineering calculation to predict flow and heat transfer of liquid sodium, attention should be paid to in some regions where the boundary layer separates or reattaches. Meanwhile, the relationship between the model coefficients should be taken into account. Besides, some physical properties have impacted with the CFD results, such as density and thermal conductivity.

Nomenclature

Data Availability

All datasets supporting the results are included in the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding this article.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (no. 11705035) and Natural Science Foundation of Heilongjiang (no. LH2019A009), which are gratefully acknowledged.