The Numerical Prediction and Analysis of Propeller Cavitation Benchmark Tests of YUPENG Ship Model

: The numerical simulation of propeller cavitation benchmark tests of YUPENG ship model is studied based on OpenFOAM, an open-source CFD (Computational Fluid Dynamics) platform, and the benchmark tests are introduced as well. The propeller cavitation shape and the hull pressure ﬂuctuation are measured and predicted, respectively. The uncertainty in hull pressure ﬂuctuation measurement is also analyzed, and the analysis showed the uncertainty is below 10%. The cavitation shape and the hull pressure ﬂuctuation predicted show quite good agreement with the observations and measured data, and the inﬂuences of the grid resolutions on the unsteady propeller cavitation and the hull pressure ﬂuctuation are investigated as well. The monotonic convergence achieved and the quite small grid uncertainty illustrate the reliability of the numerical simulation methods.


Introduction
Propeller cavitation is one important aspect of hydrodynamic performance for marine propellers, and recently it has been studied based on viscous numerical simulation methods, especially the viscous RANS (Reynolds Average Navier-Stokes) approach, and more and more research articles have been published. Salvatore [1] has studied the cavitation on the standard INSEAN E779A propeller and compared the cavitation predictions using several different solvers and cavitation models. The discrepancies in cavity extent are observed in the comparative analysis. Li [2] has also predicted the cavitating flows around the INSEAN E779A propeller in non-uniform wake, using a RANS approach and Zwart cavitation model based on commercial FLUENT. The numerical results show that the key features predicted, such as the attached cavity and the collapse of the main cavity, resemble well with the experiment observations. Paik [3] has conducted the simulations of cavitating flows for a marine propeller based on FLUENT and Schnerr-Sauer cavitation model. The simulation results show reliable performance on the prediction of cavitation shape and hull pressure fluctuation. Vaz [4] has compared more cavitation and pressure pulse predictions between different solvers and cavitation models, the results indicate that the pressure fluctuations predicted by half the participants are in reasonable agreement with measurements, but the remaining calculations predict pressure levels a factor of four or five too high. Fujiyama [5] has analyzed the unsteady cavitating flow of HSP-II and CP-II propeller at behind-hull condition both in model and full scale, using commercial software SC/Tetra v13. The results show that the unsteady propeller cavitation phenomena can be captured in the numerical calculation. Long [6] has researched the propeller cavitating flow behind the hull, analyzed the vorticity distribution and particle tracks as well, using commercial software CFX and Zwart cavitation model. The cavitation patterns predicted resemble well with the experimental observations, with some over-prediction of the cavitation area. Yilmaz [7] has investigated the hull-propeller-rudder interaction in the presence of tip vortex cavitation, using commercial software STAR-CCM+ and Schnerr-Sauer cavitation model with LES (Large eddy simulation) method. The investigation results demonstrate some capability of STAR-CCM+ for investigation in the presence of tip vortex cavitation.
In addition to commercial CFD software, an open-source CFD platform, OpenFOAM has been increasingly used to conduct the numerical simulation of propeller cavitating flows. Abolfazl [8] has performed the numerical prediction of the PPTC (Potsdam Propeller Test Case) propeller cavitation size in oblique wake based on ILES (Implicit Large Eddy Simulation) approach and Schnerr-Sauer cavitation model, using OpenFOAM. The simulation results show over-prediction of the propeller cavitation extent. Bensow [9] has carried out a prediction of a marine propeller behind a chemical tanker with ILES approach and Kunz cavitation model, based on OpenFOAM. The prediction of cavity area shows good agreement with the experiments. Zheng [10,11] has simulated the unsteady propeller cavitating flows in the stern of a ship using unsteady RANS approach, based on OpenFOAM. The predicted cavitation shape change shows quite good agreement with the experimental observation.
The present study aims to research the influence of the grid resolutions on the unsteady propeller cavitation and the induced hull pressure fluctuation behind the YUPENG benchmark ship model, which will illustrate the reliability of the numerical simulation methods.

Governing Equations
The unsteady RANS method is applied in the present study for the significantly lower computational resources needed compared with LES.
The time-averaged continuity equation and momentum equation are shown in Equations (1) and where U is the averaged velocity vector, p is the mean pressure. ρ is the fluid density, µ is the dynamic viscosity, and µ t is the turbulent viscosity. F S is the surface tension force which takes place only at the free surfaces. The turbulent viscosity µ t is modeled by the SST kw [12] turbulence model combined with wall functions. According to the VOF (Volume of Fluid) method, the physical properties of the cavitating flows are scaled by the volume fraction of water γ, with γ = 1 corresponding to pure water, γ = 0 corresponding to pure vapor, 0 < γ < 1 corresponding to the interface between water and vapor.
Thus, the density and dynamic viscosity of the cavitating flows are scaled as follows. The subscripts w and v represent water and vapor, respectively.
The transport equation of the volume fraction of water γ can be written as follows.
Here, the mass transfer rate, . m, between water and vapor is to be modeled by cavitation model. Combining with Equations (3) and (5), the mass transport, Equation (1) can be re-written as follows.
It suggests that the mass transport rate, . m, between water and vapor participates in the coupled solving procedure of pressure and velocity. During the simulations of cavitating flows, . m is calculated in the PISO (Pressure Implicit with Splitting of Operators) [13] loop at first, then the transport equation of the vapor volume fraction is solved and the standard PISO procedure is progressed.

Cavitation Model
Cavitation is a complex hydrodynamic phenomenon and is generally simplified as the transition of water into vapor in the region where the local pressure is lower than the water saturation pressure during the numerical simulation. The transition is caused by the small gas nuclei in the water. In order to mimic the phase change between water and vapor, the Schnerr-Sauer cavitation model is adopted.
Here, n 0 is the number of micro gas nuclei per water volume, and R is the initial radius of micro gas nuclei. The values of n 0 and R are default, e.g., n 0 = 1.6 × 10 13 , R = 2.0 × 10 −6 . Schnerr-Sauer cavitation model can be deduced based on Rayleigh-Plesset equation describing the bubble dynamics.

Calculation Object
The investigated object is the YUPENG benchmark ship, which is a 30,000 DWT bulk carrier driven by a fixed-pitch four-bladed propeller. The benchmark propeller cavitation model test including cavitation observation and pressure fluctuation measurement is conducted in the Large Cavitation Channel (LCC) of the China Ship Scientific Research Center (CSSRC). The principal parameters of ship and propeller model are shown in Tables 1 and 2, respectively. The hull, rudder, and propeller are visualized in Figure 1. The benchmark test setup inside the LCC is shown in Figure 2. The initial rotating angle, ϕ, is defined as 0 • at twelve o'clock, and the propeller has been assembled and oriented at the initial position. J. Mar. Sci. Eng. 2019, 7, x FOR PEER REVIEW 4 of 16 The benchmark test setup inside the LCC is shown in Figure 2. The initial rotating angle, φ, is defined as 0° at twelve o'clock, and the propeller has been assembled and oriented at the initial position.

Computational Domain and Boundary Condition
The whole computational domain is a box of 7Lpp × 4Lpp × 2Lpp. The boundary conditions on the front, bottom, and sides are velocity inlet, and the boundary condition on the top is symmetry, the boundary condition on the back is pressure outlet, as shown in Figure 3. The benchmark test setup inside the LCC is shown in Figure 2. The initial rotating angle, φ, is defined as 0° at twelve o'clock, and the propeller has been assembled and oriented at the initial position.

Computational Domain and Boundary Condition
The whole computational domain is a box of 7Lpp × 4Lpp × 2Lpp. The boundary conditions on the front, bottom, and sides are velocity inlet, and the boundary condition on the top is symmetry, the boundary condition on the back is pressure outlet, as shown in Figure 3.

Computational Domain and Boundary Condition
The whole computational domain is a box of 7L pp × 4L pp × 2L pp . The boundary conditions on the front, bottom, and sides are velocity inlet, and the boundary condition on the top is symmetry, the boundary condition on the back is pressure outlet, as shown in Figure 3. The whole computational domain consists of two sub-domains, the ship sub-domain and the propeller sub-domain. The ship sub-domain is a static region which includes the boundary of the inlet, outlet, hull, and rudder. The propeller sub-domain is a rotating region surrounding the The whole computational domain consists of two sub-domains, the ship sub-domain and the propeller sub-domain. The ship sub-domain is a static region which includes the boundary of the inlet, outlet, hull, and rudder. The propeller sub-domain is a rotating region surrounding the propeller. All the grids of the two sub-domains consists of unstructured hexahedral cells of high quality. For more accurate prediction of the flow detail, several boundary layer cells are inserted near wall surfaces. The overview of the hull, rudder, and propeller surface is shown in Figure 4. The arbitrary mesh interface (AMI) method is applied to interpolate between the non-conforming interfaces of the two sub-domains. The whole computational domain consists of two sub-domains, the ship sub-domain and the propeller sub-domain. The ship sub-domain is a static region which includes the boundary of the inlet, outlet, hull, and rudder. The propeller sub-domain is a rotating region surrounding the propeller. All the grids of the two sub-domains consists of unstructured hexahedral cells of high quality. For more accurate prediction of the flow detail, several boundary layer cells are inserted near wall surfaces. The overview of the hull, rudder, and propeller surface is shown in Figure 4. The arbitrary mesh interface (AMI) method is applied to interpolate between the non-conforming interfaces of the two sub-domains.
The velocity inlet condition is set to be a fixed value of velocity, U, at the inlet boundary, U = 5.481 m/s, and the turbulence intensity at the inflow is 0.1%. The pressure outlet condition is set to be a fixed value of pressure, which is constant based on the cavitation number, σn(0.8R), at the outlet boundary. The no-slip wall condition is set at the hull, rudder, and propeller boundary.

Calculation Result Analysis
The benchmark test conditions numerically predicted in the present study are summarized in Table 3. The velocity inlet condition is set to be a fixed value of velocity, U, at the inlet boundary, U = 5.481 m/s, and the turbulence intensity at the inflow is 0.1%. The pressure outlet condition is set to be a fixed value of pressure, which is constant based on the cavitation number, σ n(0.8R) , at the outlet boundary. The no-slip wall condition is set at the hull, rudder, and propeller boundary.

Calculation Result Analysis
The benchmark test conditions numerically predicted in the present study are summarized in Table 3. The formulae are as follows.
The free surface is neglected, and the effect of the stern wave on pressure is not taken into account either in the experiment or in the numerical calculation, because the stern wave height of the low-speed vessel (Fr < 0.2) is generally less than 0.5 m at full scale, which will be negligible at model scale in the experiment.
For the computational convergence, the numerical simulation of the steady full wet flow is conducted to obtain a quasi-stable flow field using the multiple reference frame (MRF) method, then the unsteady calculation is started to simulate the propeller rotation using sliding mesh method, and the cavitation calculation is started to predict the propeller cavitating flows by activating the cavitation model, finally. Ten revolutions are simulated after the sliding mesh method is started, and 1 • per time step is used both in unsteady and cavitation calculations.

Grid Independence
In order to investigate the grid independency, three grid sizes were chosen, as shown in Table 4 and Figure 5. All the ship and propeller regions consist of unstructured hexahedral cells with the same grid topology, generated by HEXPRESS, and the grid refinement ratio is There are 120 CPU cores used in the numerical simulation, and the time cost for the unsteady cavitation prediction in the three refined grids is recorded in Table 5.  There are 120 CPU cores used in the numerical simulation, and the time cost for the unsteady cavitation prediction in the three refined grids is recorded in Table 5. The variation of the thrust coefficient K T over the last full propeller revolution is shown in Figure 6, and the results show that it is obviously periodic.  The variation of the thrust coefficient KT over the last full propeller revolution is shown in Figure  6, and the results show that it is obviously periodic. The influence of grid resolution on KT is assessed in Table 6. Based on the procedure recommended by ITTC [14]; the grid uncertainty, UG; the corrected benchmark solution, SC; and the convergence ratio, RG, are also evaluated using the grid convergence index (GCI) approach proposed by Roache [15].
In Table 6, the KT predicted using the three grids agrees well with the measured value, and the RG obtained shows that monotonic convergence was achieved for 0 < RG < 1, the relative grid uncertainty, UG/SC, is also small, so all these predicted results illustrate the reliability of the numerical simulation methods adopted in the present study.  The influence of grid resolution on K T is assessed in Table 6. Based on the procedure recommended by ITTC [14]; the grid uncertainty, U G ; the corrected benchmark solution, S C ; and the convergence ratio, R G , are also evaluated using the grid convergence index (GCI) approach proposed by Roache [15].
In Table 6, the K T predicted using the three grids agrees well with the measured value, and the R G obtained shows that monotonic convergence was achieved for 0 < R G < 1, the relative grid uncertainty, U G /S C , is also small, so all these predicted results illustrate the reliability of the numerical simulation methods adopted in the present study. Table 6. The difference between prediction and measurement of K T ; grid uncertainty, U G ; and convergence ratio, R G .

Cavitation Shape
The predicted cavity is represented by vapor iso-surface, and it is true that different value of γ v used for the iso-surface can affect the cavitation area, e.g., a larger area with the lower value and a smaller area with the higher value. The two values of γ v , 0.5 and 0.1, are used to represent the cavitation pattern in Figure 7. The results show that the dependence is limited, and the cavitation retains the main shape. The predicted propeller cavitation using three refined grids behind the stern is compared with benchmark test sketches at the corresponding rotating angles in Figure 8. The predicted propeller cavitation using three refined grids behind the stern is compared with benchmark test sketches at the corresponding rotating angles in Figure 8.
In Figure 8, the predicted cavity using three refined grids, represented by vapor iso-surface of 0.1, shows the same dynamic behavior (e.g., the change of the propeller cavitation shape) as the benchmark test observations. The main feature, the shape change of the attached cavity along with the rotating angles, correlates well with the test sketches, e.g. the cavitation appears at about the same angle, ϕ ≈ −10 • , and reaches the maximum size at ϕ ≈ 30 • . Nevertheless, the tip vortex cavity cannot be predicted well on account of the insufficient grid resolution near the blade tip. In addition, the difference of the cavity predicted among the three grids is not obvious, and it implies that Grid 2 can obtain satisfactory efficiency and precision of the unsteady propeller cavitation. The predicted propeller cavitation using three refined grids behind the stern is compared with benchmark test sketches at the corresponding rotating angles in Figure 8.  In Figure 8, the predicted cavity using three refined grids, represented by vapor iso-surface of 0.1, shows the same dynamic behavior (e.g., the change of the propeller cavitation shape) as the benchmark test observations. The main feature, the shape change of the attached cavity along with the rotating angles, correlates well with the test sketches, e.g. the cavitation appears at about the same angle, φ ≈ −10°, and reaches the maximum size at φ ≈ 30°. Nevertheless, the tip vortex cavity cannot be predicted well on account of the insufficient grid resolution near the blade tip. In addition, the difference of the cavity predicted among the three grids is not obvious, and it implies that Grid 2 can obtain satisfactory efficiency and precision of the unsteady propeller cavitation.

Hull Pressure Fluctuation Induced by Cavitation
The hull pressure fluctuation induced by propeller cavitating flows is an important index that evaluates the propeller cavitation performance, so the monitoring points on the hull surface in the stern region are arranged as shown in Figure 9.

Hull Pressure Fluctuation Induced by Cavitation
The hull pressure fluctuation induced by propeller cavitating flows is an important index that evaluates the propeller cavitation performance, so the monitoring points on the hull surface in the stern region are arranged as shown in Figure 9. The variation of the pressure at P10 over the last full propeller revolution is shown in Figure  10, and the results show it is obviously periodic as well. The variation of the pressure at P10 over the last full propeller revolution is shown in Figure 10, and the results show it is obviously periodic as well.
(a) the direction of propeller rotation (b) the location of monitoring points Figure 9. The arrangement of monitoring points.
The variation of the pressure at P10 over the last full propeller revolution is shown in Figure  10, and the results show it is obviously periodic as well. Similar to experimental analysis, the pressure fluctuation calculated at model scale is firstly analyzed by the fast Fourier transformation (FFT) method, then converted to the one at full scale. The conversion formula is as follows.
Firstly, the hull pressure fluctuation predicted in non-cavitating condition is compared with the experiment, as shown in Figure 11.The 1BF signifies the first blade frequency component of the pressure fluctuation. Similar to experimental analysis, the pressure fluctuation calculated at model scale is firstly analyzed by the fast Fourier transformation (FFT) method, then converted to the one at full scale. The conversion formula is as follows.
∆P is = K P i ·ρ sea ·n 2 s ·D 2 s (i = 1, 2, 3 . . . .) (11) Firstly, the hull pressure fluctuation predicted in non-cavitating condition is compared with the experiment, as shown in Figure 11. The 1BF signifies the first blade frequency component of the pressure fluctuation. It indicates that the maximum value of the 1BF component of the hull pressure fluctuation is at P5 in the non-cavitating condition, and the values decrease rapidly near the propeller in the vertical plane, such as from P1 to P3, from P5 to P8, and from P9 to P12. The results achieved with Grid 1 are much better than those of the other grids. It indicates that the maximum value of the 1BF component of the hull pressure fluctuation is at P5 in the non-cavitating condition, and the values decrease rapidly near the propeller in the vertical plane, such as from P1 to P3, from P5 to P8, and from P9 to P12. The results achieved with Grid 1 are much better than those of the other grids.
For the uncertainty analysis in cavitating condition, the hull pressure fluctuation measurement has been repeated three times, and the maximum value of the 1BF component of the hull pressure fluctuation is at P10 each time. The uncertainty contribution of the measurement, the combined uncertainty, and the expanded uncertainty [16] are counted in Table 7. The analysis shows that the expanded uncertainty is below 10%. The amplitudes of the 1BF component of the hull pressure fluctuation converted at full scale are compared with the benchmark test data in Figure 12. The maximum value of the 1BF component of the hull pressure fluctuation is at P10 in the cavitating condition, while it is at P5 in the non-cavitating condition. This is mainly due to the rapid change of the cavitation volume at about ϕ ≈ 30 • , which causes the levels of the pressure pulse at P10 to be higher. Moreover, the second and third blade frequency (2BF, 3BF) components are also compared with the benchmark test in Figures 13 and 14. The obvious differences can be found, and the cause may be the instability of the propeller cavitation at the same rotating angle during various periods, which can be observed in benchmark test, but cannot be simulated accurately. In order to capture the cavity shedding at high frequency, more work needs to be done in the future, such as the adoption of LES approach. The predicted results obtained show quite good agreement with the benchmark test. It is found that not only the spatial laws but also the maximum value of 1BF component of the hull pressure fluctuation, which is mostly of great concern in engineering, can be predicted well. It seems that the calculated value of Grid 2 is a little closer to the test than that of Grid 1 at P10.
Moreover, the second and third blade frequency (2BF, 3BF) components are also compared with the benchmark test in Figures 13 and 14. The obvious differences can be found, and the cause may be the instability of the propeller cavitation at the same rotating angle during various periods, which can be observed in benchmark test, but cannot be simulated accurately. In order to capture the cavity shedding at high frequency, more work needs to be done in the future, such as the adoption of LES approach. Moreover, the second and third blade frequency (2BF, 3BF) components are also compared with the benchmark test in Figures 13 and 14. The obvious differences can be found, and the cause may be the instability of the propeller cavitation at the same rotating angle during various periods, which can be observed in benchmark test, but cannot be simulated accurately. In order to capture the cavity shedding at high frequency, more work needs to be done in the future, such as the adoption of LES approach.  The influence of grid resolution on the hull pressure fluctuation at P10 is also assessed in Table  8. The predicted ΔPS values using the three grids are all close to the experimental value, and the RG obtained shows that the monotonic convergence was achieved for 0 < RG < 1, the relative grid uncertainty UG/SC is also small, so all these predicted results illustrate the reliability of the numerical simulation methods again.  The influence of grid resolution on the hull pressure fluctuation at P10 is also assessed in Table 8. The predicted ∆P S values using the three grids are all close to the experimental value, and the R G obtained shows that the monotonic convergence was achieved for 0 < R G < 1, the relative grid uncertainty U G /S C is also small, so all these predicted results illustrate the reliability of the numerical simulation methods again.

Conclusions
The numerical simulation of propeller cavitation benchmark tests of the YUPENG ship model was conducted using the unsteady RANS approach based on OpenFOAM. The influence of grid resolutions on the unsteady propeller cavitation and hull pressure fluctuation was investigated in detail, and the conclusions obtained are as follows.
(1) The thrust coefficients predicted using the three grids are all close to the measured value, the monotonic convergence was achieved, and the grid uncertainty is also quite small. (2) The difference of the cavity predicted among the three grids is not obvious, and Grid 2 can obtain satisfactory efficiency and precision of the unsteady propeller cavitation. (3) The 1BF components of the hull pressure fluctuation predicted using the three grids are all still close to the benchmark test value, the monotonic convergence was also achieved, and the grid uncertainty is quite small as well.
There are some aspects that remain to be studied further, such as the more accurate prediction of tip vortex cavity and cavity shedding at high frequency, which is the point that our subsequent study will focus on.