*2.2. Mesh Generation and Independence Analysis*

According to the geometrical characteristics of each sub-domain, different meshing strategies are employed to construct the discrete meshes of each sub-domain. The impeller, the GV and the nozzle adopt the J-type, H-type, and O-type meshing methods, separately. By using ANSYS-ICEM, a reasonable block structure for each sub-domain is created to generate the number of structured meshes by controlling the number of mesh nodes and the growth rate. The accuracy of the calculation results will be poor if the meshes do not meet the calculation requirements. However, if the number of meshes is too large, it will occupy numerous computing resources. Therefore, the mesh independence analysis needs to be performed. The number of impeller meshes is 0.31 million, 0.5 million, 0.7 million, 0.89 million, and 1.1 million, recorded as mesh 1–5, as shown in Figure 2.

**Figure 2.** Meshes of computational domain.

The head *H* and efficiency η is calculated by Equations (1) and (2) when the flow rate is *Q*BEP, 1.13 *Q*BEP and 1.33 *Q*BEP. The relative head *H*' and relative efficiency η' are calculated by dividing the head and efficiency of mesh 4

$$H = (P\_{\rm int} - P\_{\rm outt}) / \rho \text{g} \tag{1}$$

$$
\eta = \frac{\rho \text{gQH}}{P\_{\text{shaft}}} \tag{2}
$$

where ρ is the density of water and the value is 103 kg/m3, *g* is the acceleration due to gravity and the value is 9.81 m/s2, *H* is the total head of the propulsion pump in m, *P*int is the total pressure at the entrance of the propulsion pump in Pa, *Poutt* is the total pressure on the outlet of the pump in Pa, *Pshaft* is the shaft power in kW, and η is the efficiency.

The relative head and relative efficiency are drawn in Figure 3. When the flow rate is *Q*BEP, the relative head *H*' of each mesh scheme has a certain difference, but the relative efficiency η' is basically consistent. When the flow rate is 1.33 *Q*BEP, the relative head and the relative efficiency η' has a larger difference, especially mesh 1 and mesh 2. The disparity of mesh 3 tends to decrease. As the number of meshes increases, both the relative head *H*' and the relative efficiency η' gradually enhance, and this trend exceeds to become clearer as the flow rate increases. When the number of meshes surpass 1.79 million, the relative head *H*' and the relative efficiency η' tend to be constant. Therefore, mesh 4 is chosen as the final mesh, as shown in Figure 2.

**Figure 3.** Mesh independence analysis. (**a**) Histogram of *H*' for each mesh at *Q*BEP, 1.13 *Q*BEP, and 1.33 *Q*BEP. The relative head *H*' of mesh 4 tends to be constant. (**b**) Histogram of η' for each mesh at *Q*BEP, 1.13 *Q*BEP, and 1.33 *Q*BEP. The relative head η' of mesh 4 tends to be constant.

#### *2.3. Numerical Methodolgy and Boundary Conditions*

As heat transfer does not exist in the flow process of waterjet propulsion device, Navier–Stokes equations are utilized as the governing equations to describe the flow, and the finite volume method is also applied to calculate the flow in the computational domains. Considering the incompressible flow, the inlet is set as mass flow. The impeller is the rotating domain and the spinning speed is 700 r/min or 73.27 rad/s. The reference pressure is 1 atm. The outlet is set as average static pressure. In order to guarantee the value transfer, the interfaces between each sub-domain are set, in which interfaces on the inlet and outlet of the impeller are transient rotor stator, and the remaining interfaces are general connection. Scalable wall function is processed on the wall. The standard *k*-ε turbulence model and the first-order upwind scheme are adopted. The convergence accuracy is 10<sup>−</sup>5. According to the rotating speed of the waterjet propulsion pump, the time step is 0.00023381 s and the total time is 0.685728 s. The corresponding transient turbulence convergence sample is obtained when the impeller rotates per degree. The transient turbulent convergence sample set of all time steps is obtained when the set time finishes. Usually, the numerous samples will increase the costing time and the amount of data. Therefore, this paper sets 36 samples per period; that is, the data is stored when the impeller rotated per ten degrees. Finally, 288 sample results were saved in 8 periods.

#### *2.4. Pressure Pulsation Monitoring Probe Arrangement*

Pressure pulsation can be seen as the difference between the pressure amplitude at different points in time and the average pressure amplitude over the entire time period. Pressure pulsation can usually be classified from pulsation performance and frequency. According to the pulsation performance, pressure pulsation can be divided into turbulent pressure pulsation, which ignores fluid compressibility, and pulse source pressure pulsation, which ignores fluid viscosity. Based on the pulsation frequency, the pressure pulsation can be divided into irregular random pressure pulsation, axial frequency pulsation, and blade frequency pulsation—the latter two pulsations are referred to as regular pressure pulsation. For random pressure pulsation, there are various induced factors, such as cavitation, secondary flow, non-uniform inflow, etc. In terms of the regular pressure pulsation, the main induced factors are the impeller rotation, the pump shaft rotation, and rotor–stator interaction, which are related to the axial frequency and the blade frequency, respectively. Generally, the blade frequency-related pressure pulsation is captured in the impeller. If the pressure pulsation on the upstream and downstream near the impeller shows the discipline, the rotating impeller has an effect on the flow there. The observed pressure pulsation characteristics are distinct for different pumps and for the same pump. The observed pressure pulsation characteristics are diverse under various operating conditions and monitoring positions.

Fifteen monitoring probes in the waterjet propulsion pump are shown in Figure 4, and the locations of each probe are listed in Table 1.

**Figure 4.** Diagram of pressure pulsation monitoring probe arrangement. (**a**) Lateral view. (**b**) Top view.


**Table 1.** Locations of each probe.

Fifteen monitoring probes as above are pre-set in the pre-processing. Data of monitoring parameters are stored when the impeller rotates per degree. 2880 groups of data are obtained. While obvious regular periodic diversification trend starts from the second period in the unsteady calculation process, data of the last six periods, which contains 2160 groups of data, are chosen to execute FFT transformation. The data of the eighth period on each monitoring probes are plotted into the time domain chart of pressure pulsation.

Currently, the time domain analysis method and frequency domain analysis method are the main methods in studying the pressure pulsation. The time domain chart is utilized in the time domain analysis method. In the chart, the abscissa is the time-related parameters, such as the time and period, and the ordinate is the pressure. Frequency domain analysis method converts the irregular pressure pulsation into the superposition of a simple harmonic wave with different frequencies, amplitudes, and phases by performing FFT transformation. These two methods have their own advantages. The time domain chart can intuitively reflect the change of pressure pulsation on the monitoring probe with time. The frequency domain analysis illustrates the main pulsating component of the pressure pulsation and the primary factor affecting the pressure pulsation directly. The pressure pulsation discipline of each monitoring probe on condition A, B, and C will be analyzed by using the methods above.

For the impeller rotating at 700 rev/min, the shaft frequency is 11.67 Hz by using *f* <sup>z</sup> = n/60, the blade frequency is 70 Hz by using *f* <sup>b</sup> = 6*f* z, and multiples of shaft frequency are defined as *T*f.

Pressure pulsation coefficient *Cp* is introduced to analyze the pressure pulsation characteristic of each monitoring probes, and the pressure pulsation coefficient *Cp* is calculated by applying Equation (3)

$$\mathcal{C}\_{\mathcal{P}} = \frac{(p - \overline{p})}{0.5\rho V^2} \tag{3}$$

where *p* is the instantaneous pressure in kPa, *p* is the time-averaged pressure in kPa, ρ is the density of water in kg/m3, and *V* is the blade tip speed at the entrance of propulsion pump in m/s.

#### **3. Test Arrangement and Verification**

As shown in Figure 5, double cycle waterjet propulsion test bench is established to verify the reliability of the numerical simulation method. The test device consists of two cycles—the main cycle and the secondary cycle. The main cycle, which is applied to provide the navigation speed for the waterjet propulsion pump, consists of the centrifugal auxiliary pump, electromagnetic flow meter, butterfly valve, expansion joint, rectifying device, and piping system. The secondary cycle, which is used to test the hydraulic performance, includes the test zone (mixed pump), electromagnetic flow meter, butterfly valve, and piping system. The flow rate is measured by the flow meter located in the main cycle and secondary cycle. The head is obtained by calculating the pressure measured at the beginning of the inlet flow tube and the outlet pressure measuring tube. The shaft power is calculated from the test data of the torque meter. The comprehensive error of this test bench is ±1.33%.

Numerical simulation and model test are carried out when the rotational speed is 400 rev/min. Dimensionless head *H*' and efficiency η' are obtained by dividing the head and efficiency of the best

efficiency point, plotted in Figure 6. The CFD result shows a consistent trend with the test result. Therefore, the numerical method is reliable and suitable for the following research work.

**Figure 5.** Double cycle waterjet propulsion test bench.

**Figure 6.** Contrast of CFD result and test result.

#### **4. Results**

#### *4.1. Hydraulic Performance Curve*

By diving the flow rate, head, and efficiency of the best efficiency point, dimensionless processing is performed for each condition. The dimensionless data is plotted into dimensionless a 'flow rate-head' curve and dimensionless 'flow rate-efficiency' curve, shown in Figure 7, in which the abscissa is *Q*', the left ordinate, is *H*' and the right ordinate is η'. The best efficiency point is marked as condition A. Under the condition B, the propulsion pump enters the hydraulic unstable zone. Condition C is the valley point of the hydraulic unstable zone.

**Figure 7.** External characteristic curve.

#### *4.2. Unsteady Flow Characteristics*

Condition A, condition B, and condition C are chosen to be calculated unsteadily and analyzed. Figures 8–10 list the surface streamline on the blade of three time points (0T, 1/3T, and 2/3T) in the same period for condition A, B, and C. Under the condition A, surface streamlines on the suction side of each blade are smooth and the velocity at the leading edge is high. In the same period, surface streamline on the blade is stable with the lapse of time, which indicates that the pressure field and velocity field near the blade are stable under this condition. Under condition B, the stream at the leading edge flows to the trailing edge and the tip of the blade. The pressure field and velocity field near the blade are still stable. Under condition C, an evident vortex happens on the surface of the blade and the radial location of core vortex is span = 0.65. The size and shape of the vortex change with time, which indicates that the pressure field and velocity field near the blade are unstable.

**Figure 8.** Streamline on the blade of different times for condition A (best efficiency point). (**a**) 0. (**b**) 1/3T. (**c**) 2/3T.

**Figure 9.** Streamline on the blade of different times for condition B (start point of hydraulic unstable zone). (**a**) 0. (**b**) 1/3T. (**c**) 2/3T.

**Figure 10.** Streamline on the blade of different times for condition C (valley point of hydraulic unstable zone). (**a**) 0. (**b**) 1/3T. (**c**) 2/3T.

Three turbo surfaces from the hub to the shroud are sliced and recorded as TS1 (span = 0.1) near the hub, TS2 (span = 0.65) near the intermediate surface, and TS3 (span = 0.96) near the shroud, as shown in Figure 11.

**Figure 11.** Turbo surfaces of the impeller.

Figures 12–14 are the streamlines on turbo surfaces TS1 under the condition A, B, and C, in which the velocity in the impeller is the relative velocity and the velocity in the GV is the absolute velocity. Under condition A, the streamlines are steady and vary slightly with time on the turbo surfaces of the impeller and GV. The velocity is rapid on the suction side of the impeller. The inlet attack angle is basically identical to the angle of the leading edge of the airfoil, which results in the excellent inflow condition and smooth stream in the blade-to-blade passage. A small-scale spanning vortex occurs at the tailing edge of the suction. Owing to the adjustment of the GV and the location of the vortex away from the impeller, the geometric shape and magnitude of the vortex shows no evident relationship with the time. Under condition B, slight deviation exists between the inlet attack angle and the airfoil angle of the blade in the impeller and the GV. The low-velocity region of the pressure surface of the leading edge of the impeller blade is enlarged. A large-scale spanning vortex occurs in each groove of GV, which extends from the inlet to the outlet in the axial direction, and occupies about 1/3 of the groove in the spanning direction. The streamlines of other parts in the groove are severely skewed and then gather near the outlet of the GV because of the spanning vortex. Under condition C, the smooth streamline in the groove is mildly affected. A distinct spanning vortex in the groove still exists and covers half of the groove on the spanning direction. The status and range of spanning vortex are basically maintained; however, the vortex core migrates in a small scale and the vortex status modifies, meaning the flow characteristic of GV is unstable.

**Figure 12.** Streamlines on the turbo surface TS1 of condition A (best efficiency point). (**a**) 0. (**b**) 1/3T. (**c**) 2/3T.

**Figure 13.** Streamlines on the turbo surface TS1 of condition B (start point of hydraulic unstable zone). (**a**) 0. (**b**) 1/3T. (**c**) 2/3T.

**Figure 14.** Streamlines on the turbo surface TS1 of condition C (valley point of hydraulic unstable zone). (**a**). 0 (**b**) 1/3T. (**c**) 2/3T.

Figures 15–17 are the streamlines on turbo surfaces TS2 under condition A, B, and C. Under condition A, the streamlines are smooth and no vortex occurs in the blade-to-blade passage of the impeller and GV. Under condition B, a slight deviation also occurs, but the streamlines are smooth in the groove of the impeller. As the arrows indicate in Figure 16, the streamlines deviate slightly near the tailing edge on the suction side of the impeller. The spanning vortex disappears in the groove of GV, but the shedding vortex is observed on the tailing edge of the outlet in the impeller. Part of the streamlines are severely skewed and have little impact on the mainstream. Under condition C, the inlet attack angle is consistent with the airfoil angle of the blade in the impeller. A distinct spanning vortex is observed near the tailing of different blades. The variation between the inlet attack angle of GV and the airfoil angle of the blade is apparent. Three vortexes marked as SV1, SV2, and SV3 occur on the spanning surface of GV. SV1 and SV2, located at the head edge and tailing edge of suction side, are in the channel of groove and rotate in clockwise. The range of SV1 is much larger than SV2. SV3 is the shedding vortex and located at the tailing of GV twirls on the opposite direction of the spanning vortex. As time passes, the shape and range of SV1, SV2, and SV3 varies.

**Figure 15.** Streamlines on the turbo surface TS2 of condition A (best efficiency point). (**a**) 0. (**b**) 1/3T. (**c**) 2/3T.

**Figure 16.** Streamlines on the turbo surface TS2 of condition B (start point of hydraulic unstable zone). (**a**) 0. (**b**) 1/3T. (**c**) 2/3T.

**Figure 17.** Streamlines on the turbo surface TS2 of condition C (valley point of hydraulic unstable zone). (**a**) 0. (**b**) 1/3T. (**c**) 2/3T.

Figures 18–20 are the streamlines on turbo surfaces TS2 under condition A, B, and C. Under condition A, the inlet attack angle is inconsistent with the angle of the leading edge of the airfoil. In one blade-to-blade passage, part of the stream flows from the suction side to the pressure side, then out of the blade-to-blade passage along the pressure side and the streamlines converge at the tailing edge of airfoil. The flow pattern in the blade-to-blade passage of GV is similar to the impeller. Under condition B, the stream at the entrance of the impeller flows into the neighboring groove on the opposite rotating direction instead of the GV. Part of the stream at the head edge on the suction side flows to the head edge on the pressure side of the neighboring blade and then out of the groove and into the neighboring groove after being dragged by the high-speed steam at the entrance of the groove. A low-velocity zone exists in the groove of the impeller, and the streamline is disordered. Under condition C, the distinction is huge between the inlet attack angle of the impeller and the airfoil angle of the blade, and an obvious spanning vortex appears at the head edge of the pressure side and the tailing edge of the suction side. Both spanning vortexes nearly cover the whole channel of the

groove in the spanning direction. The shape of spanning vortex will change with time, but the vortex core will not migrate and the location maintains.

**Figure 18.** Streamlines on the turbo surface TS3 of condition A (best efficiency point). (**a**) 0. (**b**) 1/3T. (**c**) 2/3T.

**Figure 19.** Streamlines on the turbo surface TS3 of condition B (start point of hydraulic unstable zone). (**a**) 0. (**b**) 1/3T. (**c**) 2/3T.

**Figure 20.** Streamlines on the turbo surface TS3 of condition C (valley point of hydraulic unstable zone). (**a**) 0. (**b**) 1/3T. (**c**) 2/3T.
