Uncertainties Analysis on the Prediction of Flow and Heat Transfer of Liquid Sodium with CFD Technology

)e 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. )e 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.)e 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. )e 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. )is 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.


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. e investigations on flow and heat transfer characteristics of liquid sodium began from 1940s [1][2][3][4][5][6].
e conclusion reached was that the flow characteristics of liquid sodium are similar to that of water [7][8][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. e 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. ere 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]. e 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. e 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. e 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]. is 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. ey 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. e 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. e 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]. e 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.

Mathematical and Physical Models
e 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][29][30] C μ � 0.09, 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. e 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, − 〈u ′ v ′ 〉 (d〈u〉/dy) � ε. According to vortex viscosity model: In the linear layer and logarithmic layer,〈u ′ v ′ 〉/k ≈ 0.3 [31], so C μ ≈ 0.09. (2) Determination of C ε2 : e equation of k and ε in uniform turbulence can be simplified as follows: (4) e attenuation of uniform turbulence conforms to exponential law: k ∼ t − n [32], and ε ∼ − nt − n− 1 can be obtained from the equation of k. e above formulas can be simplified to e value of n was measured in the grille turbulence experiment [33]; n ∈ (1.2, 1.3). However, it was adjusted later [34]: (3) Determination of C ε1 and σ ε : In the logarithmic region, the average velocity u obeys 〈u〉 � u τ ln y/κ + c. erefore, the expression of Reynolds stress is e 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. e 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. erefore, the equation of turbulent kinetic energy dissipation ε transport is d dy 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 − 〈u ′ v ′ 〉(d〈u〉/dy). 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]. e constant κ can be calculated by logarithmic law: 〈u〉 � u t ln y/κ + c. 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. is set of coefficients was first proposed by [28] and is known as the standard k-ε Science and Technology of Nuclear Installations 3 constants.
ey are determined by three measured parameters (P/ε, C ε2 , and κ) and two assumed parameters (〈u ′ v ′ 〉/k and σ k ), which are dependent on theories and experiment. e 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.
e measured quantities, such as boundary conditions and physical properties, were subject to Gaussian distribution. e 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]. e standard deviation s reflects the degree of fluctuation. e 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. e combination of hydraulic diameter and turbulence intensity was input into the FLUENT solution, and the standard wall function was used to solve the nearwall region.

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. e 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.    Figure 6: Backward-facing step geometry from Dunn [23].

Analysis of the Circular
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:

Science and Technology of Nuclear Installations
For C f , the classical experimental correlation formula of Blasius [42] was used for comparison: Re ∈ 3 × 10 4 , 10 6 .
e 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 . e small input uncertainty would be amplified in the output.  Science and Technology of Nuclear Installations 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.

Analysis of the Adiabatic Backward-Facing
Step. Boundary layer separation and reattachment in backwardfacing 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]. e Reynolds number with outlet height H as the characteristic length is 132000 (Re H � 132000). e uncertainty of the normalized flow field is shown in Figure 7. e 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][8][9]. However, in the vortex region (x/h < 5.33), the mean value of velocity was a little larger than the experimental value. e experimental value was almost coincident with the lower limit of the confidence interval. e 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. e results of uncertainty quantification are shown in Figure 10. e C f represented the friction coefficient on the side with the step. e 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:  Science and Technology of Nuclear Installations

Analysis of the Heating Backward-Facing
Step. e geometry of backward-facing step referred to Da Vià [21] is shown in Figure 11. e Reynolds number Re h took the height h of the step as the characteristic length, and the value was 4805. e height of the channel was three times that of the step; that is, E/h � 3. e L h is the heating surface with constant heat flux. e 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. e 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. e result of quantifying uncertainty of friction coefficient C f and Nu along the flow direction is shown in Figures 13 and 14, respectively. e calculated results of C f were not precise enough, especially before the vortex region. e  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. e 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. ere were similarities between C f and Nu in the rule of uncertainty propagation. e 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.
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. e 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 κ.

Conclusion
In this work, the processes of turbulent flow and heat transfer of liquid sodium in circular tube and backwardfacing step were simulated with the standard k-ε model. e results of uncertainty quantification and sensitivity analysis are as follows. e 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. e 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. e 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.
e 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 c:
Constant C f : Friction coefficient C fout : Friction coefficient of fully developed region C p : Specific heat C ε1 : k-ε model constant C ε2 : k-ε model constant C μ : k-ε model constant h: Step height k: Turbulence kinetic energy n: Decay rate Nu: Nusselt number P: Rate of production of turbulent kinetic energy Pe: Berkeley number Pr: Prandtl number Pr t : Turbulent Prandtl number Q: Wall heat flux Re:

Reynolds number s:
Deviation of the sample t: Time T: Mainstream temperature U: Inlet velocity u: Velocity along the flow channel u ′ : Fluctuation velocity along the flow channel u k : Velocity in the direction of the coordinate axis u τ : Local wall-shear velocity V: Mainstream velocity V in : Inlet velocity v ′ : Fluctuation velocity perpendicular to the flow channel x: Coordinate along the flow channel X * : Normalized reattachment point Xr: Reattachment point y: Coordinate perpendicular to the flow channel

Data Availability
All datasets supporting the results are included in the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding this article.