Development and Assessment of an Isotropic Four-Equation Model for Heat Transfer of Low Prandtl Number Fluids

In the simple gradient diffusion hypothesis, the turbulent Prandtl number ( P r t ) with a constant of 0.85 is difficult to accurately predict for liquid metals having low Prandtl numbers ( P r ), while a four-equation model can improve this solution by introducing the turbulence time-scale into the calculation of turbulent thermal diffusivity. However, the four-equation model’s transport form and numerical stability are so complex that suitable commercial code is lacking. Therefore, an isotropic four-equation model with simple Dirichlet wall boundary conditions is built in the present work. Based on the open-source computational fluid dynamics program OpenFOAM, the fully developed velocity, temperature, Reynolds stress, and heat flux of low P r fluids ( P r = 0.01–0.05) in the parallel plane are obtained by numerical simulation. The results show that the time-average statistics predicted using the present four-equation model are in good agreement with the direct numerical simulation data. Then, the isotropic four-equation model is used to analyze the flow and heat of liquid metal ( P r = 0.01) in a quadrilateral infinite rod bundle. The numerical results are compared with the various and available experimental relationships. The Nusselt numbers calculated using the isotropic four-equation model are betweenness the available correlations, while the turbulent Prandtl number model using a constant of 0.85 over predicts heat transfer. More detailed local heat transfer phenomena and distribution of low Pr fluids are obtained using the present isotropic four-equation model.


INTRODUCTION
Liquid metals are widely considered coolants with good thermal-hydraulic characteristics for many energy systems, such as fast reactors and subcritical reactors Wang et al., 2021). However, liquid metals' low molecular Prandtl number Pr leads to special heat transfer compared to traditional fluids. Typically, the Pr of lead-bismuth is 0.03-0.01 and that of sodium is 0.01-0.006 (Vià et al., 2020). It is significant to study liquid metals' turbulent heat transfer characteristics, which will affect the economy and security of the energy system. It is difficult, dangerous, and costly to experiment with liquid metals (Schroer et al., 2012;Ejenstam and Szakálos, 2015), while the computational fluid dynamics (CFD) method is more commonly applied to investigate the thermal-hydraulic properties of these liquids. However, in CFD, the computing cost required by the direct numerical simulation (DNS) method and the large eddy simulations (LES) method is too high to rely on this technique for a quick and economic calculation of heat transfer in complex geometries (Kawamura et al., 1999), while the Reynolds Averaged Navier-Stokes (RANS) method can be promoted. The velocity boundary layer and the temperature layer of traditional fluid (Pr ≈ 1) are generally considered to be similar in the RANS framework. In this way, one can obtain a constant turbulent Prandtl number Pr t to simplify the calculation of the energy equation after using the simple gradient diffusion hypothesis (SGDH) (Groetzbach, 2013). However, it is invalid for low Pr fluids (Reynolds, 1975).
In the SGDH framework, Cheng (Cheng and Tak, 2006) derived a Pr t correlation through employing global Reynolds numbers, while Kays (Kays, 1994) introduced the local turbulence effect into the Pr t relation. These nonlinear Pr t relations can effectively improve this problem in some simple geometries. However, these relations still need to be further verified in complex geometries (Duponcheel et al., 2014). Unlike the SGDH method, the differential or algebraic heat flux model (D/AHFM) establishes differential or algebraic transport equations to consider the dissimilarity between velocity and the temperature field to improve the heat transfer accuracy of liquid metals. Assessment and calibration results of D/AHFM models for low Pr fluids completed in some simple geometries show that the heat transport of second-moment closure is very sensitive to the model coefficients and functions (Lai and So, 1990;Shikazono and Kasagi, 1996;Choi and Kim, 2007;Shams et al., 2019).
Another popular model in the SGDH work, called the fourequation k-ε-k θ -ε θ turbulent heat transfer model, is introduced into turbulent and thermal time-scales for the simulation of the explicit first-order turbulent heat diffusivity. The literature (Nagano and Kim, 1988;Abe et al., 1995;Nagano and Shimada. 1996) has contributed extensively to near-wall model closure and thermal turbulence effect for a four-equation model. In recent years, Manservisi (Manservisi and Menghini, 2014a) improved the four-equation k-ε-k θ -ε θ model proposed by Abe and Nagano. The model predicted the heat transfer process of the plane, circular tube, triangular rod bundles, and quadrilateral rod bundles for Pr 0.025 fluids (Manservisi and Menghini, 2014b;Manservisi and Menghini, 2015). However, the numerical stability of the four-equation model is affected by its near-wall boundary conditions. To improve the problem, a four-equation k-ω-k θ -ω θ model with the specific dissipation rates ω = ε/(C μ k) and ω θ = ε θ /(C μ k θ ) was developed based on Manservisi's k-ε-k θ -ε θ model (Cerroni et al., 2015). Then by the following work of the literatures (Vià et al., 2016;Chierici et al., 2019;Vià and Manservivi, 2019), logarithmic specific dissipation rates of ω and ω θ , that is, Ω = ln(ω) and Ω θ = ln (ω θ ), have been used to simplify the near-wall boundary conditions and have been introduced into the logarithmic four-equation k-Ω-k θ -Ω θ model. Other state variables as proposed by Youssef (Youssef, 2006) are the velocity τ u = k/ε and temperature time-scale τ θ = k θ / ε θ to improve its numerical stability. The proposed k θ -τ θ model does not suffer from numerical stiffness problems since natural boundary conditions for the variables k θ and τ θ are used (k θ = τ θ = 0, at walls). In addition, the static variablesε and ε θ , called the isotropic dissipation rates, which are linked to the "true" dissipation rates ε and ε θ , can also provide Dirichlet wall boundary conditions and simpler transport modes than the k θ -τ θ modes (Nagano and Shimada, 1996;Nagano et al., 1997).
However, the application codes of an isotropic four-equation model in complex geometries are still lacking for liquid metals, and its reliability needs to be further verified and evaluated. So, in the present work, the four-equation model in an isotropic dissipation rate k −ε − k θ − ε θ formulation for liquid metals to improve its numerical robustness was presented by Taylor series expansion and near-wall turbulence analysis based on Abe's model and Manservisi's model. The turbulent heat transfer process of Pr = 0.01~0.05 fluids in parallel planes is numerically studied based on the finite volume method and the open-source CFD program OpenFOAM. The validity of the present four-equation model was verified by DNS data. To evaluate the present model's applicability in complex geometries, the turbulent heat transfer of Pr 0.01 fluids in a bare quadrilateral infinite rod bundle with different pitch-todiameter ratios is predicted on the present isotropic fourequation model. The numerical results obtained from two SGDH options, the Pr t 0.85 model and the present isotropic k −ε − k θ − ε θ model, are compared and analyzed with the available experimental relations and CFD results. The local distributions of dimensionless temperature, temperature fluctuation, turbulent heat diffusivity, and turbulent Prandtl numbers are analyzed.

Isotropic Four-Equation Turbulence Model
The incompressible RANS equations with no gravity and constant physical properties for the calculation of velocity, pressure, and temperature fields are considered as follows: The turbulent viscosity ] t can be solved using two-equation k − ε or k − ω models, while the turbulent thermal diffusivity α t is considered as follows: Due to low Pr fluids' physical properties, there are great differences between the momentum and heat. A constant Pr t ≈ 0.85-0.9 cannot give acceptable results for liquid metals. It is possible to consider α t as a function of the turbulence variables, such as two-equation k θ − ε θ or k θ − ω θ analogy to dynamic twoequation k − ε or k − ω. To smooth the isotropic variables at the wall, when one makesε and ε θ zero at the wall, an isotropic fourequation k −ε − k θ − ε θ model can be written using the Taylor series expansion as follows: whereε and ε θ call the isotropic dissipation rate. ε and ε θ are the "true" dissipation rate, respectively. f ε , f w , f d2 , and f wθ are wallproximity effect functions. Based on the analysis of the literature and the DNS database (Abe et al., 1995;Kawamura et al., 1998;Manservisi and Menghini, 2014a), the isotropic four-equation model coefficients for liquid metals are used and corrected as shown in Table 1. Here, the k −ε − k θ − ε θ model functions and constants recommended in this study, such as f ε , f w , f d2 , f wθ , and C d1 , are different from the k − ε − k θ − ε θ model values of Manservisi (Manservisi and Menghini, 2014b). It is noted that P kθ is the production term of temperature fluctuations k θ , which represents the energy transferred to the turbulent flux by the average temperature change rate of the RANS flow. That is, the temperature fluctuations k θ is the result of the combined action of the average temperature change rate and Reynolds heat flux. Through the transport of Eq. 11, after the temperature fluctuation occurs, it will experience convection, molecular diffusion, and turbulent diffusion until it is dissipated.
Both the turbulent Reynolds number R t k 2 /]ε and R ε are introduced into the isotropic model. Let δ be the wall distance. R ε can be recommended by Abe et al. (Abe et al., 1995), as follows: where η (] 3 /ε) 1/4 is the Kolmogorov length scale which can calculate the separated flow well. When isotropic variables are used, R t and R ε in Eqs 14-17 should become the following:

Turbulent Viscosity and Turbulent Thermal Diffusivity
Some dynamic and thermal time-scales should be applied to calculate the turbulent viscosity ] t and turbulent thermal diffusivity α t . Using dynamic and thermal time-scales, that is, τ u k/ε and τ θ k θ /ε θ , turbulent viscosity and turbulent thermal diffusivity can be defined (Abe et al., 1995;Manservisi and Menghini, 2015) as follows: and and TABLE 1 | Values for the k −ε − k θ − ε θ model constants in Eqs 9-12.
The above model has been proved to be effective in producing correct turbulence behavior of liquid metal, which draws into the mixing time-scale τ m τ u R/(C m + R) and the ratio of the velocity time-scale to the thermal time-scale R τ θ /τ u (Abe et al., 1995). The values of Pr t∞ 0.9 and C m 0.3 are recommended for liquid metals . In order to establish a complete isotropic four-equation turbulent heat transfer model k −ε − k θ − ε θ , we used the isotropic variablesε and ε θ in Eqs 20-26 instead of the "true" variables ε and ε θ .

Wall Boundary Conditions for Turbulence Models
Appropriate boundary conditions should be applied to solve the four-equation turbulent models without the available wall functions. By using series expansions as shown in Table 2, the near-wall behaviors of dynamical and thermal turbulence variables k, ε, k θ , and ε θ can be obtained (Deng et al., 2001) as follows: According to the definitions of Eqs 7, 8, when δ tends to zero, k,ε, k θ , and ε θ tend to zero. Thus, the simple Dirichlet wall boundary conditions with the isotropic k −ε − k θ − ε θ model can be imposed.

Modified Navier-Stokes Equations for Periodic Boundary Conditions
For a fully developed turbulent field, the pressure gradient term in Eq. 2 is divided into the mean of the pressure gradient dP f /dz which is constant along the flow direction z and the fluctuating pressure gradient zP/zx i . The mean term of the pressure gradient is given as the known power source term.
For a fully developed temperature field, temperature T is written as follows: where T in is the inlet temperature, and ΔT b is a constant wall temperature gradient, while θ satisfies periodic boundary conditions along the flow direction. According to the energy conservation for thermal fully developed flow with a wall flux q w , ΔT b can be set as 4q w /(ρC p u b D h ), where u b is average velocity, and D h is an equivalent diameter. After introducing periodic temperature θ, Eq. 3 can be written as follows: 3 NUMERICAL RESULTS AND ANALYSIS

Numerical Solver of the Isotropic Four-Equation Model
For the modified RANS system and the isotropic four-equation k −ε − k θ − ε θ model in Section 2, the appropriate discrete format could be selected according to the needs. Here, the convection term adopts the Gauss upwind format, the diffusion term uses the Gauss linear format, and the time term uses the steady-state format. The SIMPLE algorithm solves the continuity and momentum simultaneously, and coupled multigrid iteration is used for the matrix solution. The convergence conditions of residual error are as follows: where Q stands for u i , θ, k,ε, k θ and ε θ . The index i denotes the steps of calculation. All calculations were realized on OpenFOAM (Weller et al., 1998;Moukalled et al., 2016). The detailed analysis of the four-equation model solver by Gu and Su (Gu and su, 2021) based on OpenFOAM can be further referred to. The height of the grid point closest to the wall is approximately 10 −3 mm to meet the requirements of the isotropic four-equation model (y + < 1).

Numerical Verification of the Plane
The full development process of different Pr fluids in a constant heat flow heated plane was studied using DNS (Kawamura et al., 1998;Kawamura et al., 1999;Tiselj and Cizelj, 2012). This study compares the DNS data with Pr = 0.01 0.05. The plane geometry and mesh parameters are  consistent with those in the literature . Here, the three-dimensional parallel plane with a flow length of 6.4 h, a spanwise width of 3.2 h, and a plane height of 2 h were selected as the calculation domain, where h is the half-height of the plane, 30.25 mm. The wall heat flux is set as q w = 360 kW/m 2 . Finally, 80, 50, and 100 nodes are divided into three directions, respectively: Flow direction, span direction, and height. Due to the symmetry of the plane, half of the plane is taken as the calculation domain, and the structured grid is used. The grid height of the first boundary layer is set to 3 × 10 −6 m to meet the requirements of the low Reynolds number turbulence model, y + < 1. The physical properties of Pr = 0.01-0.05 are reported in Table 3 where u τ is the friction velocity on the wall side, h is the halfheight of the plane, and the numbers A~G represent different Reynolds number calculation conditions. For the convenience of analysis, u τ and the friction temperature T τ are taken as the dimensionless reference quantity, where T τ is set as q w ρ −1 Cp −1 u −1 τ , and q w is the wall heat flux.

Velocity Field Verification
The dimensionless velocity u + u/u τ distribution of cases A and F along the dimensionless wall-normal distance y + yu τ /v is shown in Figure 1. At a low Reynolds number and a high Reynolds number, the velocity distribution predicted by the present isotropic four-equation model is in good agreement with the DNS results, and the velocity distribution meets u + ＝ y + in the linear region (u + ＝y + ) while it meets u + ＝2.5ln y + +5 in the logarithmic rate region (y + ＞30).
The distribution of dimensionless Reynolds stress τ + R uυu −2 τ and total shear stress τ + total τ + R + ]υ ,y u −2 τ along y + under cases A and F is shown in Figure 2. The stress calculation results are in good agreement with DNS results. In the near-wall region (1＜ y + ＜7), this study predicted τ + R to be lower than the DNS value. However, it has little effect on the calculation results of the velocity field in the linear region, mainly because the molecular diffusion is dominant in the near-wall region. With the increase in wall distance, the turbulent diffusion effect increases gradually. When y + ＞30, the total shear stress coincides with the Reynolds stress curve.

Temperature Field Verification
The dimensionless temperature θ + θ/T τ distribution of the Pr = 0.01 fluid along y + is shown in Figure 3A. Cases A, B, and C agree with DNS and meet θ + ＝Pry + in the linear region (1＜ y + ＜40).
The distribution of dimensionless Reynolds heat flux q + The heat flux calculation results of cases A-C are in good agreement with DNS results. Due to the thick thermal boundary layer and strong molecular heat conduction of low Pr fluid, although the predicted Reynolds heat flux deviates from the DNS results at 1＜y + ＜30, it does not affect the linear distribution of the temperature field in this region. With the increase in wall distance, although the effect of turbulent heat diffusion is gradually increasing, there is still a difference between the total heat flow and Reynolds heat flux, mainly caused by the strong molecular heat conduction of low Pr fluid.
2) Pr = 0.025 The distribution of θ + /Pr of the Pr 0.025 fluid along y + is shown in Figure 4A. Cases A, B, and C are in good agreement with DNS. The peak value of θ + /Pr increases with the increase in the Reynolds number. The distribution of q + R along y + is shown in Figure 4B. Cases A and B agree with DNS. With the increase in Re, the q + R peak increases and moves to the turbulent core.

3) Pr=0.05
As shown in Figure 5A, the dimensionless temperature θ + of case A of Pr = 0.05 fluid is in good agreement with the DNS results. At the same Reynolds number, with the decrease in Pr, the molecular thermal conductivity increases and the average θ + decreases. As shown in Figure 5B, the dimensionless Reynolds

Numerical Analysis of the Square Bundle
To further analyze the applicability of the current isotropic fourequation model in complex geometries, the flow and heat of liquid metal with a Prandtl value of 0.01 in the quadrilateral infinite rod bundle were studied. Figure 6 shows the schematic diagram of the computational domain and local hexahedral mesh of the quadrilateral infinite rod bundle. Lines ab, bo, oc, and arc ca are the local distributions to be analyzed next. The flow parameters are reported in Table 4. The surfaces of rods in contact with the fluids are heated by uniform wall fluxes q w . Symmetry condition is applied on other faces to simulate the infinite rod bundle region with a quadrilateral arrangement. The heat fully developed flow length L z is set to 10D h (Ge et al., 2017). The rod diameter D and X = P/D ratio are the same as in Zhukov's experiment (Zhukov et al., 2002).

Heat Transfer Evaluation and Analysis
In      reviewed in the study by (Mikityuk 2009). Zhukov (Zhukov et al., 2002) obtained 36 Nu(Pe) data from the heat transfer experiment of 22%Na-78%K in a quadrilateral rod bundle. The correlations obtained by Subbotin (Subbotin et al., 1965) and BREST (Adamov 2001) were derived from triangular bundles, while Mikityuk, from a wide database, conducted experimental results of liquid metals (Mikityuk 2009). It is worth noting that there are no more available experimental data. The correlations should be used with care and be paid more attention to for analysis of the reliability. The Nu can be calculated as follows: where n i is a vector of a unit perpendicular to a sliced face. First, the simulations for Pr = 0.01 are characterized by X = 1.25, 1.34, and 1.46 with a corresponding hydraulic diameter D h of 11.87, 15.43, and 20.57 mm. Nine simulations with bulk Reynolds numbers of 250-400,000 were operated. The corresponding Peclet number is approximately 250-4,000. Figure 7 shows the Nusselt number results for X 1.25 -1.46 cases. In   between these empirical relations is too poor to relate because these relationships have very different slopes.
In Figures 7A-C, the prediction Nu of the isotropic k −ε − k θ − ε θ model almost lies between the experimental relationship. The prediction of the Pr t 0.85 model is higher than Subbotin and Mikityuk correlations. At low Peclet numbers, the Nu results having a special slope predicted using the isotropic four-equation model is closer to the Subbotin and Mikityuk correlations, while the calculation results at high Pe are conservative, similar to the Zhukov and BREST correlations. The overall numerical Nu seems to be more conservative than these experimental correlations. Poor experimental relevance makes this problem's numerical verification difficult to go on. Figure 7D, the numerical slope in the Nu(Pe) line obtained using the k −ε − k θ − ε θ model is very similar to that predicted using Manservisi's k − ε − k θ − ε θ model, which has successfully predicted heat transfer calculation for Pr 0.025 fluids and can provide an additional CFD relationship reference. Therefore, we still use the present isotropic fourequation model for thermal analysis and local temperature calculations.

As shown in
This present work focuses on the comparison of the effects of two turbulent thermal diffusion models (the Pr t 0.85 model and the present isotropic k −ε − k θ − ε θ model) on the temperature field, without considering the influence of different RANS turbulence k − ε or k − w models. The detailed   Figure 8 shows the dimensionless temperature θ + λ(T − T b )/(q w · D h ) on a sliced surface of square bundle X = 1.25. The overall dimensionless temperature in the channel decreases with the increase in Pe. The temperature at the rod surface is not constant. The wall temperature near the center of the channels is lower than that near the center of gaps. Figure 6 has annotated the line ab, line bo, line oc, and arc ca, where points b and c represent the center of the gap and the channel, respectively. The dimensionless temperature field with Pe = 2000 is selected for analysis. In Figure 9A, with the decrease in P/D, the gap width decreases, and the overall dimensionless temperature distribution θ + at the gap line ab increases. In Figures 9B,D, the coolant dimensionless temperature distribution from gap center point b to channel center point o and the wall dimensionless temperature distribution over arc ca are smoother in pace with strengthening of P/D. Due to the  Figure 9C.

2) Fluctuation and isotropic dissipation
The distribution and shape of k θ and ε θ for Pe = 100, 2000, and 4,000 are shown in Figure 10. k θ has a lower local value near the center point o. The maximum value appears near the wall. The maximum value of k θ decreases when Pe increases, while that of ε θ increases. The numerical shape and distribution of k θ and its isotropic dissipation ε θ obtained from the present k −ε − k θ − ε θ model is similar to that predicted by Manservisi's k − ε − k θ − ε θ model . Further details on over line ab, line bo, and line oc can be found in Figure 11. The dimensionless root square means temperature θ rms λ 2k θ /(qD h ) is set. From wall point a to wall point b, θ rms first increases to a small peak and then begins to decrease to point b. Then along the bo line, the maximum peak value of θ rms is reached near the midpoint of the bo section. Three peaks of different sizes were formed on lines ab, bo, and oc, respectively. When it goes straight toward the wall point, it will eventually approach 0. When the value is X, the overall value of θ rms decreases.
3) Turbulence and molecular thermal diffusivity Figure 12 shows the proportional relationship of turbulent to molecular thermal diffusivity, that is, dimensionless thermal diffusivity a + α t /α, over line ab, line bo, and line oc with the change of P/D for square bundle X = 1.25 at Pe = 2000. Reynolds heat flux q t and molecular heat flux q m can be defined as follows: Thus, the specific value of turbulent to molecular thermal diffusivity represents the ratio of Reynolds to molecular heat flux. Reynolds heat flux is mainly caused by thermal diffusion caused by turbulent flow, while molecular heat conduction generates molecular heat flux. In Figure 12, in the area closest to the wall point a or c, a + tends to 0, which means turbulent heat diffusion tends to be ignored here. With the increase in wall distance from wall point a to gap center point b, the turbulent heat diffusion gradually increases, so the value of a + increases. However, a + is still less than one over line ab, which means that the molecular heat conduction of liquid metal always affects the heat transfer at the gap area. At the gap center point b, a + reaches a local maximum over  line ab. It increases along the bo segment and reaches its maximum near the middle of the bo segment. In most regions far from the wall over lines bo and oc, due to the development of the turbulent flow, turbulent diffusion is greater than the molecular heat conduction, so a + > 1. Interestingly, the maximum value of a + does not appear at the channel center point o but near point o. At the gap line ab, the smaller the gap length, the greater the influence of molecular heat conduction, so the value of a + increases over the line ab when the valve of P/D increases. But it decreases over the line bo and line oc when the valve of P/D increases because at the same Pe, the larger the P/D, the larger the flow area, the smaller the overall velocity, and the smaller the turbulent heat diffusion over those areas. In Figure 13, the dimensionless thermal diffusivity on a slice for Pe = 250, 1,000, 2000, and 4,000 is shown. At the same P/D, when the Pe increases, turbulent heat diffusion will be enhanced, and the grid point with a + > 1 begins to appear.

4) Turbulent Prandtl number
In Figure 14, the Pr t with the change of P/D is shown over the lines ab, bo, and oc for the square bundle with X = 1.25 at Pe =2000. The Pr t is not constant and linear. It is larger than the general value of 0.85 for liquid metal. The Pr t distribution on the line bo is to be smoother than that of lines ab and co. The value of Pr t decreases over line ab when the value of P/D increases, while it increases over lines bo and oc. This phenomenon echoes the changing trend of dimensionless thermal diffusivity. In Figure 15, the turbulent Prandtl number on a slice for Pe = 1,000, 2000, and 4,000 for the square bundle is shown. When Pe increases, Pr t near the wall and the center decreases gradually.
The average turbulent Prandtl number Pr tm can be defined as follows:

Pr tm
A Pr t dA A dA (34) Table 6 summarizes the calculation results of average turbulent Prandtl number Pr tm under different P/D and Pe using the present k −ε − k θ − ε θ model. One can see that at the same P/D, the average turbulent Prandtl number decreases with the increase in Pe, while at the same Pe, it increases when the value of P/D increases.

5) Performance analysis of the present four-equation model
When the same turbulence k −ε model is kept to solve the turbulent viscosity, the present four-equation model needs to solve two more equations k θ − ε θ than the Pr t = 0.85 model. To analyze the calculation speed of the k −ε − k θ − ε θ model, the single-step iterative calculation time and calculation memory of these two models are compared and listed in Table 7. The current test platform is HP ProDesk 680 G 4 MT with an Intel(R) Core(TM) i7-9700 CPU at 3.00 GHz and 16 GB system memory at 2,667 MHz. At present, the model solver is written based on OpenFOAM V6 version, compiled and integrated in Ubuntu 18.04 system environment. As shown in Table 7, the calculation time and memory required by the  k −ε − k θ − ε θ model are greater than those of the k −ε − Pr t 0.85 model.

CONCLUSION
The present work studied an isotropic four-equation model for low Pr number that uses simple Dirichlet wall boundaries. The isotropic four-equation model can provide more references for calculating the thermal-hydraulic phenomena of liquid metals. But its wider applicability needs further verification.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
XS: concept, research, writing, editing, code, and data processing. XL: modification, concept, research, and code. XW: research, modification, and code. YL: concept, editing, and research. QC: editing and research. QS: editing and research. XS: editing and research. GL: funding, project management, concept, and research.

FUNDING
This work was supported by the Research on key technology and safety verification of primary circuit, Grant No. 2020YFB1902104, and the Experimental study on thermalhydraulics of fuel rod bundle, Grant No. Y828020XZ0.