*2.1. Turbulence Model and Cavitation Model*

During the development of cavitation flow, the fluid was treated as a vapor liquid mixture. The gas-liquid two-phase mixture model (Mixture model) was used for two-phase flow calculations in this paper. The mixture model assumes that the fluid is homogeneous, and the two-phase fluid components are assumed to share the same velocity and pressure. The continuity equation

of vapor/liquid, momentum equation, and vapor volume fraction mass transport equation for the two-phase flow are generally expressed as follows:

$$\frac{\partial \rho\_m}{\partial t} + \frac{\partial (\rho\_m u\_j)}{\partial x\_j} = 0 \tag{1}$$

$$\frac{\partial \rho\_{\text{m}} u\_{i}}{\partial t} + \frac{\partial (\rho\_{\text{m}} u\_{i} u\_{j})}{\partial \mathbf{x}\_{j}} = \rho f\_{i} - \frac{\partial p}{\partial \mathbf{x}\_{i}} + \frac{\partial}{\partial \mathbf{x}\_{j}} \Big[ (\mu\_{\text{m}} + \mu\_{t}) \Big\{ \frac{\partial u\_{i}}{\partial \mathbf{x}\_{j}} + \frac{\partial u\_{j}}{\partial \mathbf{x}\_{i}} - \frac{2}{3} \frac{\partial u\_{i}}{\partial \mathbf{x}\_{j}} \delta\_{ij} \Big\} \tag{2}$$

$$\frac{\partial(\alpha\_{\mathbb{P}}\rho\_{V})}{\partial t} + \frac{\partial(\alpha\_{\mathbb{P}}\rho\_{V}u\_{j})}{\partial \mathbf{x}\_{j}} = R\_{E} - R\_{\mathbb{C}} \tag{3}$$

where *p* is the mixture pressure, Pa; ρ*<sup>V</sup>* is vapor volume density, kg/m3; *ui* is the velocity in *i* direction, m/s; *uj* is the velocity in *j* direction, m/s; *fi* is the body force in the i direction; μ is the laminar viscosity; *t* is the time, s; μ*<sup>t</sup>* is the turbulent viscosity; α*<sup>v</sup>* is the volume fraction of vapor; and *RE* and *RC* respectively represent the source terms for evaporation and condensation, kg/(m·s). The mixture density ρ*<sup>m</sup>* is defined as follows:

$$
\rho\_{\rm ll} = \rho \nu \alpha\_{\rm v} + \rho\_{l} (1 - \alpha\_{\rm v}) \tag{4}
$$

where ρ*<sup>l</sup>* is the liquid volume density.

The RNG *k-*ε (Renormalization Group *k-*ε) model can simulate the region of cavitation clearly and has better adaptability to the simulation of cavitation flow [36–38]. Therefore, under the assumption of mixture homogeneous flow model, the RNG *k-*ε model was selected to calculate the cavitation characteristics of the waterjet propulsion system.

The cavitation model is a mathematical model that describes the mutual transformation between the liquid volume and the vapor volume. The transport equation model is the most commonly used cavitation model, which mainly includes three types: the Zwart–Gerber–Belamri model [39], the Kunz model [40], and the Schnerr–Sauer model [41]. In this paper, the Zwart model in ANSYS-CFX software was selected to simulate and analyze the cavitation process.

$$R\_{\epsilon} = F\_{\text{trap}} \frac{3\alpha\_{\text{Func}} (1 - \alpha\_{\upsilon}) \rho\_V}{R\_B} \sqrt{\frac{2}{3} \frac{(P\_V - P)}{\rho\_l}} \text{ P} < \text{P}\_{\text{V}} \tag{5}$$

$$R\_c = F\_{cond} \frac{3\alpha\_v \rho\_V}{R\_B} \sqrt{\frac{2}{3} \frac{(P - P\_V)}{\rho\_l}} \text{ P} > \text{P}\_V \tag{6}$$

where *PV* is the vapor pressure; *P* is the flow field pressure; *Fvap* and *Fcond* are empirical coefficients for the vaporization and condensation processes, respectively; α*ruc* is the non-condensable gas fraction in the liquid; and *RB* is the typical bubble size in the water. According to numerous literature discussions [42–44], <sup>α</sup>*ruc* takes the value of 5 <sup>×</sup> 10−4; *RB* takes the value of 1 <sup>×</sup> 10−<sup>6</sup> m; and *Fvap* and *Fcond* take the values 50 and 0.01, respectively.

## *2.2. Geometric Model and Mesh Generation*

As can be seen from Figure 1 that the propulsion pump is the core component of waterjet propulsion pump system, which consists of an impeller with six blades and a guide vane with seven vanes. Therefore, in order to ensure that the pump system has better anti-cavitation ability, it is necessary to require the pump itself to have a good anti-cavitation performance. As shown in Figure 1, the waterjet propulsion pump system is the main research object. The calculation domain of the waterjet propulsion pump system includes import extension, inlet passage, impeller, guide vane and nozzle. The inlet flow of the inlet passage is a non-uniform flow, so the inlet section needs to be extended by a distance to guarantee the accuracy and convergence of the calculation, as shown in the import extension in Figure 1a. The basic geometrical parameters of the waterjet propulsion pump

system are shown in Table 1. The rotational speed of the waterjet propulsion pump is 700 r/min. *D0* represents the inlet diameter of the impeller. Section P1 and Section P2 are respectively the inlet and outlet sections for calculating the head of the pump section.

**Figure 1.** Three-dimensional schematic diagram of the waterjet propulsion pump system. 1. Nozzle; 2. Guide vane; 3. Impeller; 4. Inlet passage; 5. Import extension; 6. Guide vane blades; 7. Impeller blades; 8. Guide vane hub; 9. Impeller hub; 10. Shaft; P1. Impeller inlet section; P2. Guide vane outlet section.



The entire computational domain was divided into three parts to generate mesh, namely, the nozzle part, the propulsion pump part and the inlet passage part. Structured mesh with hexahedral cells can improve the computational efficiency of CFD (Computational Fluid Dynamics), so the meshes of the computational domain in this paper are divided into structured mesh by ICEM CFD software (ANSYS Inc., Pittsburgh, PA, USA). ICEM CFD is the integrated computer engineering and manufacturing code for computational fluid dynamics, which is a professional preprocessing software. Since both the impeller and guide vane are periodic meshes, the generated single channel meshes are rotated and duplicated to generate the computational domain of impeller and guide vane. Since different turbulence models have different requirements for grid Y plus values, the RNG *k-*ε model requires *y*+ values between 30 and 100. The grid size of the boundary layer was controlled to ensure that Y plus meets the requirements of turbulence model in this paper. The grid of the computational domain is shown in Figure 2.

**Figure 2.** Grid of the computational domain by using ICEM CFD software.

The number of grids has a great influence on the calculation accuracy and the solution speed. In theory, the denser the calculation domain grid is, the higher the calculation accuracy is. However, in the actual calculation process, too many grids will greatly increase the calculation period and waste computing resources. In order to find the appropriate mesh size, a mesh sensitivity analysis was carried out. For this calculation model, the impeller is the core of the whole calculation domain. The mesh sensitivity of the propulsion pump section was verified by changing the size of the impeller mesh. As shown in Figure 3, when the grid number of the impeller reaches 1.6 <sup>×</sup> 106, the head of the pump section changes less. Equation (6) was used to calculate the head of the pump section. The calculated cross sections are sections P1 and P2 in Figure 1. Finally, the number of meshes of the impeller is 1.72 million, and the total number of calculation domain grids of the waterjet propulsion pump system is 5.23 million.

$$H = \frac{P\_{2-2} - P\_{1-1}}{\rho \mathcal{g}} \tag{7}$$

where P1-1 and P2-2 are the pressure of the sections P1 and P2, respectively, Pa; ρ is water density, kg/m3; and g is gravitational acceleration, m/s2.

**Figure 3.** Mesh sensitivity analysis.

## *2.3. Setting of the Boundary Condition*

The unsteady cavitating turbulent flow was simulated using the high-performance computational fluid dynamics software (ANSYS CFX 14.5, ANSYS Inc., Pittsburgh, PA, USA) based on the Reynolds-averaged Navier–Stokes (RANS) equation with the RNG *k-* ε turbulence model. The total pressure condition was applied as outlet boundary condition, and mass flow rate was adopted as the outlet boundary condition. The process of cavitation was controlled by changing the total pressure of the inlet. Assuming that there are no bubbles in the inlet fluid, the volume fractions of the liquid phase and the vapor phase at the inlet were set to 1 and 0, respectively. No slip condition was applied at solid boundaries. The effect of temperature was not considered in the calculation [45]. The mixture model was chosen as the multiphase flow model, the Zwart model was chosen as the cavitation model, and the Rayleigh–Plesset equation was used to control bubble motion. Before the calculation of the cavitation model was embedded, the pump segment under the condition of no cavitation was calculated, and the calculation result was used as the initial value of the cavitation simulation to improve the convergence speed and calculation accuracy. When applying the cavitation model, two physical parameters were given: the vaporization pressure of the liquid at normal temperature (25 ◦C) *PV* = 3574 Pa, and the surface tension of the cavitation bubble, 0.074 N/m.
