Evaluation of the Turbulence Model Influence on the Numerical Simulation of Cavitating Flow with Emphasis on Temperature Effect

Yilin Deng 1,2, Jian Feng 2, Fulai Wan 2, Xi Shen 2 and Bin Xu 2,* 1 Institute for Energy Research, Jiangsu University, 301 Xuefu Road, Zhenjiang 212013, China; yldeng@ujs.edu.cn 2 Research Center of Fluid Machinery Engineering and Technology, Jiangsu University, 301 Xuefu Road, Zhenjiang 212013, China; 2211711042@stmail.ujs.edu.cn (J.F.); 2111711013@stmail.ujs.edu.cn (F.W.); shen.xi.01@outlook.com (X.S.) * Correspondence: norkistar@ujs.edu.cn; Tel.: +86-1358-400-3613


Introduction
Cavitation is an important issue that affects the operation and working life of fluid machinery. Cavitation refers to the process of the formation, development, and collapse of vapor cavities in the liquid or at the liquid-solid interface when the local pressure in the liquid decreases [1]. The pressure reduction and unsteady characteristics of the cavitation area usually lead to mechanical vibration and noise, resulting in negative effects such as equipment surface fatigue damage, fracture, and mechanical performance degradation [2][3][4].
When cavitation occurs, the fluid medium will absorb heat from the surrounding fluid due to the latent heat of vaporization during the transition from the liquid phase to the vapor phase, resulting in a decrease in the temperature of the cavitation zone-that is, the thermal effect on cavitation [5]. At normal temperature, the cavitation of water can ignore the influence of thermodynamic effects-that is, when the cavitation occurs, it is assumed that the temperature of the cavitation area does not change [6,7]. However, the physical properties of thermosensitive fluid such as liquid hydrogen, liquid nitrogen, and liquid oxygen are sensitive to temperature changes. During the cavitating flow, the thermodynamic effect is significant, which causes the temperature of the cavitation region to change, which changes the cavitating flow field dramatically [8][9][10].
Fluids such as liquid hydrogen and liquid oxygen are often used as propellants for liquid carrier rocket engines. Because the weight and size of the engine turbo pump are subject to strict design requirements, it is necessary to increase the engine thrust by increasing the power density of the turbo pump. At the same time, in order to ensure effective carrying capacity, reducing the volume of the propellant container will cause the pressure at the entrance of the inducer to decrease. Under the conditions of high-speed rotation and lower inlet pressure, cavitation can easily occur around the

Governing Equations
The thermosensitive cavitating flow is treated as the homogeneous multiphase flow. The continuity, momentum, and energy equations in the Cartesian coordinates are shown as follows: where ρ m = ρ v α v + ρ l α l is the density of multiphase flow, subscripts v and l mean vapor phase and liquid phase, respectively; subscripts i and j represent the coordinates in the Cartesian coordinate system; α is the volume fraction; u is the flow velocity; p is the pressure in the flow field; µ m is the dynamic viscosity; µ t is the turbulent dynamic viscosity; k e f f is the effective thermal conductivity; T is temperature; C p is specific heat; and L is the latent heat.

Turbulence Model
For the turbulence model, how to predict the turbulent eddy viscosity in the cavitation area is a key problem, affecting the accuracy of numerical simulation results. In this study, DCM and FBM are employed to modify the turbulent eddy viscosity for three turbulence models (k − ε, RNG k − ε, and SST k − ω). The implementation in detail is shown below.

k − ω Turbulence Model
This model is a classic two-equation turbulence model with good applicability and high utilization frequency. It has been embedded in many CFD solvers. The specific form of the equation is as follows: For the turbulent kinetic energy equation, For the turbulent dissipation rate equation, For the turbulent viscosity, In the formula, C ε1 , C ε2 , σ k , σ ε , and C µ are all empirical coefficients, taking the values of 1.44, 1.92, 1.0, 1.3, and 0.09, respectively. P kb and P εb are the influence source terms of buoyancy, and P k is the turbulent kinetic energy generation term. When there is a strong vortex motion and a large fluid pressure gradient in the flow field, the prediction accuracy of k-ε turbulence model is limited. Compared with the k-ε turbulence model, the RNG k-ε model considers the rotating flow. The time-averaged strain rate is also considered in the turbulent energy dissipation rate equation.
For the turbulent kinetic energy equation, For the turbulent dissipation rate equation, In the formula, the value of each constant is: The standard k-ε model and its improved RNG k-ε model have higher simulation accuracy for the turbulent flow that has been fully developed on the far wall surface, but the accuracy of flow simulation with flow separation phenomena and large backpressure gradients needs to be improved. Menter [28] proposes a two-equation SST k-ω model, which is obtained by modifying the (baseline) BSL k-ω model. The simulation accuracy of separated flow and strong curvature flow is greatly improved, but there are problems such as high requirements on grid resolution, difficulty in convergence, and sensitivity to boundary conditions. The equation expression is as follows: For the turbulent kinetic energy equation, For the turbulence frequency equation, For the turbulent viscosity equation, where P kb is the turbulent kinetic energy term produced by buoyance; P k is the turbulent kinetic energy term produced by viscous force; F 1 and F 2 are the blending functions; and S is the shear tensor. The constant coefficients are taken separately: a 1 = 0.31, σ k = 2, β' = 0.09, σ ω = 2, α = 5/9, β = 0.075, and σ ω2 = 1/0.856. The expressions of F 1 and F 2 are as follows: where V is kinematic viscosity.

Modification of Turbulence Model
Considering the influence of turbulence scale and the compressibility of gas-liquid mixed phase, the filter-based model and density correction method are employed to modify the turbulent viscosity of three turbulence models, which was used for the following numerical study: where C 3 is the empirical coefficients, which is assigned to be 1.0.

Cavitation Model for Thermosensitive Fluids
With the evolution of the cavitation process, the mass transfer rate between liquid phase and vapor phase is controlled by the combination of source term m + and sink term m − , respectively. To capture the cavitation characteristics, liquid volume fraction and vapor volume fraction needs to be obtained. A representative method is to employ a transport equation to determine the liquid volume fraction and vapor volume fraction, respectively. For the Merkle cavitation model [29], the transfer rate between the liquid phase and vapor phase is supposed to be proportional to the local pressure difference. The governing equation for the Merkle cavitation model is shown below: where C vap and C cond are the characteristic constant coefficients for the evaporation and condensation, respectively; U ∞ is the reference speed at infinity; t ∞ is the reference time; and p v is the saturated vapor pressure. The default values for C vap and C cond are listed as follows: The original Merkle cavitation model [29] is established based on the isothermal assumption. The saturated vapor pressure in the model is a fixed value, while the saturated vapor pressure of thermosensitive fluid is actually a function of temperature. In our previous study, the B-factor theory and Antoine equation were taken into consideration. Considering the influence of turbulent kinetic energy on the saturated vapor pressure, the modified source term is shown as follows: .
For the implementation of the corrected turbulence model and cavitation model in the numerical simulation, the CEL language assembly function is employed to embed them in CCL files into the CFX software.

Numerical Setup and Validation
The numerical results were validated by a recognized cavitation experiment conducted by Cervone et al. [30]. Water at different temperatures around the NACA0015 hydrofoil were investigated experimentally in the Cavitating Pump Rotordynamic Test Facility (CPRTF) laboratory. The influence of thermodynamic effect on the cavitation performance was analyzed. The experimental data was published in the American Society of Mechanical Engineers (ASME) journal, which attracted great attention from the research community.

NACA0015 Hydrofoil Geometry Model
The structure of NACA0015 hydrofoil used in the experiment is shown in Figure 1. The hydrofoil chord has a length of 115 mm and a width of 80 mm. The attack angle of the experimental hydrofoil ranges from 4 • to 8 • . There are three holes for pressure measurement at the inlet and outlet of the experimental section, ten holes for pressure measurement along the suction surface of the hydrofoil, and two holes for pressure measurement along the pressure surface of the hydrofoil. The model is built using the UG software. The 3D model of NACA0015 hydrofoil is shown in Figure 2.

Mesh Implementation
The ICEM CFD software is employed to mesh the NACA0015 hydrofoil into a hexahedral structured grid. The mesh grid around the NACA0015 hydrofoil is shown in Figure 3. The mesh quality of the boundary layer around the hydrofoil has a great influence on the accuracy of numerical simulation results. The y + value is usually employed to judge the mesh quality [31][32]. There is no specific range for the distribution of the y + value for different turbulence models. The y + value is usually maintained under 60 to keep the accuracy of the simulation results in the research community. In order to improve the quality of the mesh grid, a refined C-shaped structure grid is used around the leading edge and trailing edge of the hydrofoil. The distribution of the y + value on the upper and lower surfaces of the hydrofoil is shown in Figure 4.

Mesh Implementation
The ICEM CFD software is employed to mesh the NACA0015 hydrofoil into a hexahedral structured grid. The mesh grid around the NACA0015 hydrofoil is shown in Figure 3. The mesh quality of the boundary layer around the hydrofoil has a great influence on the accuracy of numerical simulation results. The y + value is usually employed to judge the mesh quality [31,32]. There is no specific range for the distribution of the y + value for different turbulence models. The y + value is usually maintained under 60 to keep the accuracy of the simulation results in the research community. In order to improve the quality of the mesh grid, a refined C-shaped structure grid is used around the leading edge and trailing edge of the hydrofoil. The distribution of the y + value on the upper and lower surfaces of the hydrofoil is shown in Figure 4.

Mesh Implementation
The ICEM CFD software is employed to mesh the NACA0015 hydrofoil into a hexahedral structured grid. The mesh grid around the NACA0015 hydrofoil is shown in Figure 3. The mesh quality of the boundary layer around the hydrofoil has a great influence on the accuracy of numerical simulation results. The y + value is usually employed to judge the mesh quality [31][32]. There is no specific range for the distribution of the y + value for different turbulence models. The y + value is usually maintained under 60 to keep the accuracy of the simulation results in the research community. In order to improve the quality of the mesh grid, a refined C-shaped structure grid is used around the leading edge and trailing edge of the hydrofoil. The distribution of the y + value on the upper and lower surfaces of the hydrofoil is shown in Figure 4.

Boundary Conditions
The boundary conditions for the numerical setup is shown in Figure 2. The inlet is set to be the velocity inlet boundary condition, the outlet is set to be the pressure outlet boundary condition, and the wall is set to be the non-slip boundary condition. The specific values at different temperatures are listed in Table 1.

Mesh Independence Study
Three sets of meshes with different densification levels, namely mesh 1, mesh 2, and mesh 3, are employed to perform the mesh independence study. The mesh information is shown in Table 2. To facilitate comparison and verification, the dimensionless pressure coefficient is defined as follows: (29) where is pressure at outlet; is speed at inlet; and is the density of the working fluid. The number of grids in the three sets (mesh1, mesh2 and mesh3) are 1.1 million, 2.85 million, and 4.56 million, respectively. As shown in Figure 5, the simulated pressure coefficients along the

Boundary Conditions
The boundary conditions for the numerical setup is shown in Figure 2. The inlet is set to be the velocity inlet boundary condition, the outlet is set to be the pressure outlet boundary condition, and the wall is set to be the non-slip boundary condition. The specific values at different temperatures are listed in Table 1.

Mesh Independence Study
Three sets of meshes with different densification levels, namely mesh 1, mesh 2, and mesh 3, are employed to perform the mesh independence study. The mesh information is shown in Table 2.
To facilitate comparison and verification, the dimensionless pressure coefficient p c is defined as follows: where p out is pressure at outlet; U in is speed at inlet; and ρ l is the density of the working fluid. The number of grids in the three sets (mesh 1, mesh 2 and mesh 3) are 1.1 million, 2.85 million, and 4.56 million, respectively. As shown in Figure 5, the simulated pressure coefficients along the upper surface of NACA0015 hydrofoil match well with the experimental results. Taking the calculation cost into consideration, mesh 2 is selected in the numerical simulation.

Influence of Different Turbulence Models on NACA0015 at 25°С
Using three turbulence models and their corresponding modified turbulence models, the numerical simulation about the cavitating flow around the NACA0015 hydrofoil is performed. The distribution of pressure coefficient along the suction side of the NACA0015 hydrofoil is obtained as shown in Figure 6. The results show that there is a certain difference between the simulated pressure coefficient and experimentally measured pressure coefficient in the low-pressure area of the first 30% of the hydrofoil length. At the last 70% of the hydrofoil length, the calculation results of the RNG k-ε and SST k-ω turbulence models are closer to the experimental values, and the simulation results of the k-ε model differ greatly from the experimental results. Based on the modified turbulence model, the simulation results of the revised RNG k-ε model and the revised SST k-ω showed significant improvement, which is closer to the experimental value. The simulation results of the revised k-ε model are significantly smaller than that of the uncorrected turbulence model. The correction effect is obvious.

Influence of Different Turbulence Models on NACA0015 at 25 • C
Using three turbulence models and their corresponding modified turbulence models, the numerical simulation about the cavitating flow around the NACA0015 hydrofoil is performed. The distribution of pressure coefficient along the suction side of the NACA0015 hydrofoil is obtained as shown in Figure 6. The results show that there is a certain difference between the simulated pressure coefficient and experimentally measured pressure coefficient in the low-pressure area of the first 30% of the hydrofoil length. At the last 70% of the hydrofoil length, the calculation results of the RNG k-ε and SST k-ω turbulence models are closer to the experimental values, and the simulation results of the k-ε model differ greatly from the experimental results. Based on the modified turbulence model, the simulation results of the revised RNG k-ε model and the revised SST k-ω showed significant improvement, which is closer to the experimental value. The simulation results of the revised k-ε model are significantly smaller than that of the uncorrected turbulence model. The correction effect is obvious.

Influence of Different Turbulence Models on NACA0015 at 25°С
Using three turbulence models and their corresponding modified turbulence models, the numerical simulation about the cavitating flow around the NACA0015 hydrofoil is performed. The distribution of pressure coefficient along the suction side of the NACA0015 hydrofoil is obtained as shown in Figure 6. The results show that there is a certain difference between the simulated pressure coefficient and experimentally measured pressure coefficient in the low-pressure area of the first 30% of the hydrofoil length. At the last 70% of the hydrofoil length, the calculation results of the RNG k-ε and SST k-ω turbulence models are closer to the experimental values, and the simulation results of the k-ε model differ greatly from the experimental results. Based on the modified turbulence model, the simulation results of the revised RNG k-ε model and the revised SST k-ω showed significant improvement, which is closer to the experimental value. The simulation results of the revised k-ε model are significantly smaller than that of the uncorrected turbulence model. The correction effect is obvious.  The root-mean-square (RMS) error and maximum deviation of the pressure coefficient for the numerical and the experimental results were analyzed to demonstrate the improvement of the calculation results of the three turbulence models and their modified turbulence models. The RMS error and the maximum deviation are shown in Figure 7. The results showed that the RMS error and the maximum deviation of the RNG k-ε model and the SST k-ω model are less than that of the k-ε model; meanwhile, the RMS error and the maximum deviation of the modified turbulence model showed good improvement.  At room temperature, the distribution of vapor volume fractions on the surface of NACA0015 hydrofoil is shown in Figure 8. The results showed that the k-ε model has a large-scale vortex at the tail of the cavity, and the development process of the cavitation core area is longer. Meanwhile, the development process of RNG k-ε model and SST k-ω model in the cavitation core area is shorter. The modified k-ε model eliminates the vortex at the tail of the cavity, but the cavitation core area is significantly shorter. The cavitation area is significantly expanded based on the modified RNG k-ε model and the modified SST k-ω model and the vapor volume fraction in the cavitation core area is larger.  At room temperature, the distribution of vapor volume fractions on the surface of NACA0015 hydrofoil is shown in Figure 8. The results showed that the k-ε model has a large-scale vortex at the tail of the cavity, and the development process of the cavitation core area is longer. Meanwhile, the development process of RNG k-ε model and SST k-ω model in the cavitation core area is shorter. The modified k-ε model eliminates the vortex at the tail of the cavity, but the cavitation core area is significantly shorter. The cavitation area is significantly expanded based on the modified RNG k-ε model and the modified SST k-ω model and the vapor volume fraction in the cavitation core area is larger.  At room temperature, the distribution of vapor volume fractions on the surface of NACA0015 hydrofoil is shown in Figure 8. The results showed that the k-ε model has a large-scale vortex at the tail of the cavity, and the development process of the cavitation core area is longer. Meanwhile, the development process of RNG k-ε model and SST k-ω model in the cavitation core area is shorter. The modified k-ε model eliminates the vortex at the tail of the cavity, but the cavitation core area is significantly shorter. The cavitation area is significantly expanded based on the modified RNG k-ε model and the modified SST k-ω model and the vapor volume fraction in the cavitation core area is larger.

Influence of Different Turbulence Models on NACA0015 at 50 • C
Through the discussion in Section 4.1, the modified RNG k-ε model and the modified SST k-ω model show good applicability. This section focuses on the modified RNG k-ε model and the SST k-ω model for the numerical simulation of cavitating flow around the NACA0015 hydrofoil at 50 • C and compares it with the experimental results.
At 50 • C, the pressure coefficient distribution of the NACA0015 hydrofoil suction surface simulated by the modified RNG k-ε model and the modified SST k-ω model is shown in Figure 9. The revised RNG k-ε model calculates the low-pressure region in a small range, while the revised SST k-ω model calculates the low-pressure region in a large range. The RMS error and the maximum deviation are shown in Figure 10. The results show that at 50 • C, the RMS error and maximum deviation of the modified SST k-ω model are both less than that of the modified RNG k-ε model.

Influence of Different Turbulence Models on NACA0015 at 50°С
Through the discussion in Section 4.1, the modified RNG k-ε model and the modified SST k-ω model show good applicability. This section focuses on the modified RNG k-ε model and the SST kω model for the numerical simulation of cavitating flow around the NACA0015 hydrofoil at 50°C and compares it with the experimental results.
At 50°С, the pressure coefficient distribution of the NACA0015 hydrofoil suction surface simulated by the modified RNG k-ε model and the modified SST k-ω model is shown in Figure 9. The revised RNG k-ε model calculates the low-pressure region in a small range, while the revised SST kω model calculates the low-pressure region in a large range. The RMS error and the maximum deviation are shown in Figure 10. The results show that at 50°С, the RMS error and maximum deviation of the modified SST k-ω model are both less than that of the modified RNG k-ε model.   Figure 11 shows the distribution of the pressure coefficient along the suction surface of the NACA0015 hydrofoil simulated by the modified RNG k-ε model and the modified SST k-ω model at 70°С. The results show that the simulation using revised RNG k-ε model matches well with the experimental results. The modified SST k-ω model has a longer cavitation development process. Figure 12 shows that the RMS error and maximum deviation of the modified RNG k-ε model are less than that of the modified SST k-ω model. The corrected RNG k-ε model is more accurate at 70°С.

Influence of Different Turbulence Models on NACA0015 at 50°С
Through the discussion in Section 4.1, the modified RNG k-ε model and the modified SST k-ω model show good applicability. This section focuses on the modified RNG k-ε model and the SST kω model for the numerical simulation of cavitating flow around the NACA0015 hydrofoil at 50°C and compares it with the experimental results.
At 50°С, the pressure coefficient distribution of the NACA0015 hydrofoil suction surface simulated by the modified RNG k-ε model and the modified SST k-ω model is shown in Figure 9. The revised RNG k-ε model calculates the low-pressure region in a small range, while the revised SST kω model calculates the low-pressure region in a large range. The RMS error and the maximum deviation are shown in Figure 10. The results show that at 50°С, the RMS error and maximum deviation of the modified SST k-ω model are both less than that of the modified RNG k-ε model.   Figure 11 shows the distribution of the pressure coefficient along the suction surface of the NACA0015 hydrofoil simulated by the modified RNG k-ε model and the modified SST k-ω model at 70°С. The results show that the simulation using revised RNG k-ε model matches well with the experimental results. The modified SST k-ω model has a longer cavitation development process. Figure 12 shows that the RMS error and maximum deviation of the modified RNG k-ε model are less than that of the modified SST k-ω model. The corrected RNG k-ε model is more accurate at 70°С.  Figure 11 shows the distribution of the pressure coefficient along the suction surface of the NACA0015 hydrofoil simulated by the modified RNG k-ε model and the modified SST k-ω model at 70 • C. The results show that the simulation using revised RNG k-ε model matches well with the experimental results. The modified SST k-ω model has a longer cavitation development process. Figure 12 shows that the RMS error and maximum deviation of the modified RNG k-ε model are less than that of the modified SST k-ω model. The corrected RNG k-ε model is more accurate at 70 • C.   Figure 13 shows the vapor volume fraction distribution on the surface of NACA0015 hydrofoil calculated by the modified RNG k-ε model at different water temperatures. The results show that as the temperature increases, the vapor volume fraction decreases, the cavity area decreases, and the vapor-liquid interface becomes blurred. The water vapor content decreases with the increase of temperature.   Figure 13 shows the vapor volume fraction distribution on the surface of NACA0015 hydrofoil calculated by the modified RNG k-ε model at different water temperatures. The results show that as the temperature increases, the vapor volume fraction decreases, the cavity area decreases, and the vapor-liquid interface becomes blurred. The water vapor content decreases with the increase of temperature.   Figure 13 shows the vapor volume fraction distribution on the surface of NACA0015 hydrofoil calculated by the modified RNG k-ε model at different water temperatures. The results show that as the temperature increases, the vapor volume fraction decreases, the cavity area decreases, and the vapor-liquid interface becomes blurred. The water vapor content decreases with the increase of temperature.   Figure 13 shows the vapor volume fraction distribution on the surface of NACA0015 hydrofoil calculated by the modified RNG k-ε model at different water temperatures. The results show that as the temperature increases, the vapor volume fraction decreases, the cavity area decreases, and the vapor-liquid interface becomes blurred. The water vapor content decreases with the increase of temperature.

Conclusions
This paper carries out numerical investigation of different turbulence model effect on the cavitating flow around NACA0015 hydrofoil with water at different temperatures. The k-ε, RNG k-ε, and SST k-ω turbulence model and their revised turbulence model are studied systematically. The following conclusions can be drawn from this research paper: (1) At 25 • C, the correction effect is significant for the modified k-ε model, and the vortex is eliminated in the closed area of the cavity tail. The simulation results obtained from the modified RNG k-ε model and the SST k-ω model showed reasonably good agreement with the experimental results. (2) At 50 • C, the modified RNG k-ε model and the modified SST k-ω model have a small difference between numerical results and experimental results for the RMS error and the maximum deviation. (3) At 70 • C, the modified RNG k-ε model is smaller than the result of the modified SST k-ω model in terms of the RMS error and the maximum deviation. The turbulent kinetic energy of the modified SST k-ω model near the wall is significantly larger than that obtained by the modified RNG k-ε model, and the cavitation is more serious, which is quite different from the experimental results.