Abstract
Propeller exciting forces are the main causes of stern vibrations. In this paper, three-dimension exciting bearing forces of one blade and the whole propeller under nonuniform ship’s wake were predicted, and the influence of cross flows on these exciting forces was studied. All simulations were carried out using a commercial solver, STAR-CCM+. To obtain the nominal wake for studying propeller exciting forces, flow field around a bare hull was simulated. Numerical results were widely validated by measured data, especially the velocity field at the propeller plane. Harmonic characteristics of the nonuniform ship’s wake were studied. Then, a propeller under uniform inflow and nonuniform ship’s wake with/without cross flows was simulated. Free-water surface and hull boundary were considered using a specially designed dummy stern. Results show that the influence of cross flows on propeller exciting forces is obvious. As for the exciting forces of one blade, the cross flows have greater influence on the axial force. As for the exciting forces of the whole propeller, the cross flows have greater influence on the transverse and vertical forces, and if the cross flows in ship’s wake are not considered, the amplitudes of the main harmonics of transverse and vertical forces increase obviously.
1. Introduction
With the increasing pressure of market competition, both the load capacity of the ship and the main engine power are increased. The propeller is operated under heavy loading conditions. Besides, under the influence of nonuniform ship’s wake, propeller exciting forces, including bearing forces and surface forces, may cause severe structural vibrations. In this paper, only the bearing forces were studied, which may cause blade and shaft vibrations.
In experiments, the nonuniform ship’s wake is needed for measuring propeller exciting forces. The entire ship model with a model-sized propeller can be used in large cavitation tunnels or depressurized tanks to perform the measurement experiments. But in small or medium-sized cavitation tunnels, the most common method is still using the wire grid alone or using the wire grid after a dummy stern to simulate the wake [1]. If the wire grid is used alone, then only the axial velocity is simulated, and the measured results need to be revised to consider the cross flows. If the wire grid is used with a dummy stern, the simulated cross flows may still be different from that in the actual ship wake. Besides, the flow field to the propeller is very likely to be affected by the wall boundaries of small tanks or tunnels. Therefore, it is important to study the influence of cross flows in the wake field on propeller exciting forces [2].
Along with the fast development of computational fluid dynamics, numerical methods become an important auxiliary tool for experimental studies. There are two commonly used methods for numerically studying propeller exciting forces under nonuniform ship’s wake. One is that the propeller is simulated separately from the hull using measured ship’s wake [3–5] or simulated ship wake as velocity inlet boundary condition. Here, the measured ship’s wake usually refers to the nominal wake because the effective wake is difficult to measure in experiments for now. In [6, 7], the three velocity components of the measured wake were directly used in the inlet boundary for propeller simulations. In this method, the accuracy of predicted propeller exciting forces is mainly dependent on the measured (simulated) wake data. Actually, small turbulent structures in ship’s wake are difficult to measure using a limited number of measurement points. But using this method, the influence of cross flows on propeller exciting forces can be investigated. The other is to simulate the rotating propeller right after the hull [8–10]. This is the way to predict the most accurate propeller exciting forces. However, the only shortcoming of this method for the present study is that the influence of cross flows cannot be separated from the predicted propeller exciting forces. Therefore, in this paper, the propeller simulations for studying propeller exciting forces were separated from the flow field simulation of the ship, and the predicted nonuniform ship’s wake was used as a velocity inlet boundary condition for propeller simulations. Literatures [11–13] provided references for the numerical simulations of the ship, and the measured results of the flow field around the KCS ship [14] were used to validate the numerical predictions of the bare hull in this paper. To consider the influence of the free-water surface and hull boundary on propeller exciting forces, a specially designed dummy stern was used for the propeller simulation.
The present study focused on predicting three-dimension exciting bearing forces of one blade and the whole propeller under nonuniform ship’s wake using numerical methods and studying the influence of cross flows on these exciting forces. All simulations were carried out using a professional commercial solver, the STAR-CCM+. First, the unsteady Reynolds-averaged Navier-Stokes (URANS) and Detached Eddy Simulation (DES) were compared for the flow field simulations of the ship to obtain the nonuniform ship’s wake for propeller simulations. Grid independent verification and validation of the numerical methods for the flow field simulations of the ship were conducted. Ship resistance, global ship waves, wave profile, and velocity distributions at the propeller plane were compared with the measured results. Harmonics of the nonuniform ship’s wake was studied. Second, the predicted ship’s wake with/without cross flows was applied for propeller simulations to predict three-dimension exciting forces of one blade and the whole propeller. Based on the experimental data, numerical methods for propeller simulations under uniform inflow and nonuniform ship’s wake were validated. Last, frequency analysis of propeller exciting forces was conducted. The influence of cross flows in ship’s wake on the mean values and the harmonics of the axial, transverse, and vertical exciting forces of one blade and the whole propeller was studied.
2. Numerical Methods for the Flow Field Simulations of the Ship
2.1. Physical Models
A model-sized 3600-TEU container ship (KCS) designed by the Korea Research Institute of Ships and Ocean Engineering (KRISO) was simulated to obtain the nonuniform ship’s wake. Large numbers of experiments of this ship were conducted by the Maritime and Ocean Engineering Research Institute (MOERI), National Maritime Research Institute (NMRI), and Schiffbau Versuchsanstalt Potsdam GmbH (SVA). Related researches about the ship can be found in the Tokyo 2005/2015 CFD Workshop and the Gothenburg 2010 Workshop. The principal particulars of the model-sized KCS and the model propeller KP505 are presented in Table 1. The geometrical model of the ship and the propeller are shown in Figure 1.

2.2. Numerical Setup
Dimensions of the computational domain in the - plane and the main boundary conditions are presented in Figure 2. The length between the side boundaries and the ship centerline is , where is the ship length. Uniform inflow with the speed of 2.196 m/s is assigned to the velocity inlet boundary. The entire ship hull is set as no-slip boundaries. For a better simulation of the turbulence after the ship, vortices shedding from the hull and ship waves, multiple blocks are designed to improve local mesh resolutions in these areas of interest (Figure 3). Automatic meshing approach is applied to discretize the computational domain (Figure 4). There are six layers of boundary-layer meshes near the hull and the average wall + value is about 60. A simple multiphase model, volume of fluid (VOF) [13], is applied to resolve the interface (waves) between water and air. When the waves reach the boundary of the virtual tank, they will be reflected back and interfere with the newly generated waves; then the simulated ship waves are different from the actual waves. Therefore, the wave damping function [15] is used to block the reflected waves.



Both URANS and DES solvers were applied for the flow field simulations around the ship. All simulations were performed with the STAR-CCM+ and the version number is 10.04.011. The code solves continuity equations in integral form by means of the finite volume technique. The second-order temporal scheme was used. The simulation is first carried out using URANS solver with the SST (shear stress transport) -Omega model [16]. The time step size is set to 0.04 second, and the maximum inner iterations within each time step is set to 10. The residuals monitoring plot was presented in Figure 5(a). It is clear that each residual quantity decreases to a small number within each time step. From the physical quantities monitoring plot presented in Figure 5(b), the simulation needs about 60 more seconds of physical time before the pressure resistance coefficient () of the ship showing the sign of convergence. It is clear that the friction resistance coefficient () of the ship is much easier to converge than . After the calculation running for 100 seconds of physical time, DES solver is selected for the simulation. The time step size is set to seconds. About 20 seconds of physical time is needed before these coefficients showing a clear sign of convergence. Detailed information about the numerical setup and mathematical principles can be found in the User’s Guide of the STAR-CCM+ [17].

(a)

(b)
3. Grid Refinement Study and Validation
3.1. Ship Resistance
Based on the predicted and of the ship, the grid refinement study is done by increasing the mesh resolution. Three sets of meshes with different mesh resolutions (the increasing rate of the base element size is 1.2) are simulated using both URANS and DES solvers. It should be noted that the distribution of the meshes within the boundary-layer in the spanwise direction keeps basically the same no matter the mesh is coarse or fine, because the wall + value should be kept within the recommended range to use the standard wall function to predict wall shear stress accurately. Thus, increasing the mesh resolution refines the meshes in the streamwise direction within the boundary-layer. The aspect ratio of each grid cell within the boundary-layer is reduced and the shape is close to a cube, which will yield more accurate predictions. Here, a scalar is defined for the grid refinement study (here and after):
Results from Table 2 show that the differences in and are less than 1% when the cell number increases from 3.12 million to 4.39 million. It suggests that the numerical results are little affected by the mesh resolution when the cell number is up to 4.39 million.
The predicted and of the ship by URANS and DES (the 4.39-million grid case) are validated by experimental data [14]. The following formula is used to quantify how well the predicted results agree with experimental data (here and after):where is the numerical results and is the experimental data. It is clear that the predicted ship resistances by URANS and DES are close to the experimental data, and all deviations are within 6% (Table 3). It should be noted that the viscous force of the ship predicted by DES is more accurate than that by URANS.
3.2. Global Ship Waves and Wave Profile
Global ship waves predicted by URANS of the cases with different mesh resolutions are presented in Figure 6. Generally, it is quite difficult to distinguish the difference among the three wave figures, while, when looking closer, you can still find slightly more details after the ship stern area in the finest mesh case (more contours).

Global ship waves predicted by URANS and DES (using the finest mesh case) are compared in Figure 7. The measured global ship waves [14] are presented in Figure 7(a). It is quite clear that the waves predicted by DES contain much more details than that predicted by URANS. For a better comparison with the measured results, Table 4 lists the wave heights corresponding to different levels defined in Figure 7(a). Qualitatively and quantitatively, the global ship waves predicted by DES are closer to the measured results than that predicted by URANS.

(a)

(b)
The predicted wave profiles along the hull by URANS and DES with different mesh resolutions are presented in Figure 8. Generally, they are very similar and agree well with the measured results [14].

3.3. Velocity Distribution at the Propeller Plane
The predicted axial velocity () contours at the propeller plane right after the ship with different mesh resolutions are present in Figure 9. With the benefit of the grid type and the meshing techniques (all + values are kept within the recommended range and the meshes near the propeller plane are greatly refined), it is clear that the predicted results are very little affected by mesh resolutions.

The definition of the circumferential degree (), tangential velocity (), and radial velocity () at the propeller plane is present in Figure 10. A comparison between the predicted three-dimension velocities at the propeller plane by URANS and DES (the finest case) is presented in Figure 11. The tangential and radial velocities (cross flows) are shown as vectors. It is clear that the tangential and radial velocities are much smaller than the axial velocity. Qualitatively, the axial, tangential, and radial velocities at the propeller plane predicted by URANS and DES have similar patterns. From the axial velocity in Figure 11, the viscous effect of the hull is obvious, and there is a clear sign of vortices near the hub.


A quantitative comparison of the axial, tangential, and radial velocities of , , and predicted by URANS and DES is presented in Figure 12. By comparing with the measured results [14], the axial velocities predicted by DES are slightly better than that predicted by URANS. The axial velocity agrees well with the measured results except that of , where large difference appears below the hub (0~60 degrees). The main reason is that the flow close to the hub belongs to the TBL (turbulent boundary-layer) problems, and the velocities within the boundary-layer are difficult to predict accurately, especially for the flow around a hull. Considering that these places of relatively large differences are close to the blade root, the contribution of these velocities to propeller exciting forces is small (the main contribution is from the region of 0.75R~0.9R of the blade). Therefore, the influence of the prediction error of the axial velocity on propeller exciting forces is negligible. From Figure 12, although the tangential and radial velocities are much smaller than the axial velocity, there are considerable values at and , and the cross flows may have influence on the harmonics of propeller exciting forces.

From the radial velocity () in Figure 12, some significant changes occur at the vicinity of 150 degrees: the radial velocity changes direction in these places. This phenomenon is clearly caused by vortices. To validate the estimated vortices from Figure 12, the commonly used -criteria [18] is adopted for the visualization of vortex structures in the field. The -criteria are mathematically described aswhere is vorticity of the fluid and is shear strain rate of the fluid. A field function of these -criteria is provided by the STAR-CCM+, and it can be directly applied for this study. The simulated vortices are shown in Figure 13. Clearly, two pairs of strong vortices (A and B) are found near the hub. One pair of vortices (A) appears above the hub at the vicinity and the other pair of vortices (B) is below and very close to the hub (about 0.2R~0.3R). The positions of the pair of vortices (A) explain the reason why the radial velocity changes directions at the vicinity of 150 degrees.

3.4. Harmonic Wake Analysis
The circumferential distribution of the axial wake coefficient of profile at the propeller plane is analyzed harmonically:where denotes the mean wake coefficient and are always equal to zero theoretically because the wake flow is symmetric.
Results of the harmonic analysis of two predicted wakes and the measured wake of the ship are presented in Figure 14, where the number 0 in the horizontal axis denotes the mean wake coefficient. It is clear that the wake of the ship is primarily composed of low harmonics, especially the first and second harmonics, which are decided by the symmetrical hull form. Amplitudes of the high harmonics are quite small, which are caused by turbulent structures in the field (including vortices shedding from the hull). Results show that the hull form plays a decisive role on the harmonics of the wake. Thus, a well-design hull form is very important for reducing propeller exciting forces. Generally, the high harmonics can be ignored if there is no special need. However, if high-frequency structural vibration is the most concerned problem, then the high order harmonics of the wake should be studied. From Figure 14, the amplitudes of low harmonics of the wakes predicted by URANS and DES are close to that based on the measured wake. However, the amplitudes of the high harmonics of the predicted wake by DES are much larger than that of the measured wake. The possible reason is that the number of points used for measuring velocity is not enough to capture those small turbulent structures in the wake.

(a)

(b)
4. Numerical Studies of Propeller Exciting Forces
4.1. Grid Refinement Study and Validation of Numerical Methods for Propeller Simulations
Based on the measured open water performance of KP505, validation of the numerical methods for propeller simulation in open water was carried out. Dimensions of the computational domain and the main boundary conditions are shown in Figure 15, where is the diameter of the propeller. The computational domain is decomposed into two main parts: internal rotating domain and external static domain. The Moving Reference Frame model [19], a steady-state approach, is applied in the rotating domain for simulating the rotating propeller. A stationary reference is used in the static domain. Several working conditions (, known as advance coefficient) are simulated by changing the inflow speed () based on the following equation: . The rate of the propeller is kept the same ( Hz). The uniform inflow speed () is applied to the velocity inlet boundary. Automatic meshing approach is applied to discretize the computational domain. Fifteen layers of boundary-layer meshes are generated near the blade (Figure 16). Low + wall treatment is applied to resolve the viscous sublayer of the blade. All wall + values on the blade surface are basically within 1. Only five layers of boundary-layer meshes are generated near the hub, because the viscous force of the hub is not interested. The URANS method with SST -omega model is used to model turbulence around the propeller.


Grid refinement study is performed based on the predicted propeller thrust coefficients (). Three sets of meshes with different mesh resolutions are simulated. Numerical results of are presented in Table 5, where are calculated based on (1). When the cell number increases from 3.69 million to 4.12 million, the relative percent difference of is within 3%. It suggests that the numerical results are little affected by the mesh resolution when the cell number is up to 4.12 million.
The predicted of the propeller (the 4.12 million grid case) with different advance coefficients () are validated by the measured results from NMRI, shown in Table 6, where the deviations between the predicted results and the measured data are calculated based on (2). The predicted are close to the measured data when the advance coefficient is less than 0.8 (the deviations are within 3%). It should be noted that there are large deviations when the advance coefficient is great. This is because the is too small when the propeller works at such a great advance coefficient (, ), which makes the relative percentage difference look large.
4.2. Numerical Methods for Propeller Simulations under Nonuniform Ship’s Wake
Numerical simulations of the propeller under nonuniform ship’s wake were carried out to predict propeller exciting forces. The computational domain is presented in Figure 17 and the main dimensions are similar to that shown in Figure 15. To accurately predict propeller exciting forces, the influence of the free-water surface and wall boundary was considered, and the flow field that the propeller comes across needs to be the same as that at the propeller plane right after the ship. Thus, a dummy stern is designed in the computational domain (Figure 17). The section of the stern at the velocity inlet boundary for propeller simulations is the exact section of the ship’s hull at the propeller plane, and each section of the dummy stern is the same. In this way, the wall boundary right above the rotating propeller is just like that right after the ship. The velocity field at the velocity inlet boundary is delivered from a circular region (the diameter is ) at the propeller plane right after the ship, as shown in Figure 18. To study the influence of cross flows on these exciting forces, two cases were simulated: one is that only the axial velocity is assigned to the velocity inlet boundary, and the other is that all velocity components, axial, tangential, and radial velocities, are assigned to the velocity inlet boundary. It should be noted that, in the propeller simulation using all three velocity components, the inlet boundary is away from the propeller plane. Sliding mesh technique [11] is used in the rotating domain for simulating the rotating propeller. A transient simulation is needed for obtaining propeller exciting forces. The time step size is set as s, which is the time of the propeller rotating one degree. The maximum iterations value within each time step is 10. The simulation is performed for one second of physical time and the propeller rotates 9.5 revolutions ( Hz). The URANS method with SST -Omega model is used to predict propeller exciting forces.


Figure 19 shows the axial velocity contour near the propeller. It is clear that the interaction between the propeller and the dummy hull is strong. But as for propeller bearing forces, the velocity field within the propeller plane disc is the key. From Figure 19, the velocity condition set at the inlet boundary kept the same for a distance; then a clear acceleration appears near the propeller, which is caused by the rotating propeller (induced velocities).

Figure 20 shows the gauge (static) pressure contour of the entire domain. It is clear the gauge pressure in the air is basically the same (the pressure is close to zero relatively to standard atmospheric pressure), while that in the water is increased linearly along the negative direction. Complex pressure fluctuations appear near the rotating propeller.

Figure 21 shows the three-dimension velocity distributions at the velocity inlet boundary and a section upstream of the propeller plane. It is clear that the velocities at the section are only slightly different from that at the velocity inlet boundary. These small differences are caused by the induced velocities of the rotating propeller and also some numerical dissipation effects. More importantly, the main characteristics of the nonuniform ship’s wake are kept.

Time history curves of the axial force of one blade () under the wake with/without cross flows are shown in Figure 22. It is clear that when the propeller rotates the third revolution, the simulation is already convergent.

Figure 23 shows the predicted fluctuating thrusts of three cases with different mesh resolutions. It is clear that the mesh resolution has less effect on the numerical results (the relative percent difference between the numerical results of the middle case and the fine case is about 0.7%). Table 7 shows a comparison of the numerical result of the mean thrust (the finest case) with the measured data in the self-propulsion experiment [20]. The deviation is −2.24%, which is much smaller than that of the numerical study in [20], where the deviation is 11.65%.

4.3. Propeller Exciting Forces under Ship’s Wake considering Cross Flows
The following simulations for studying propeller exciting forces were performed using the fine mesh case.
As for the blade structural vibration, the exciting forces of one blade were studied. The predicted nondimensionalized axial force, , transverse force, , and vertical force, , of one blade under ship’s wake with/without cross flows are shown in Figure 24 (transverse and vertical forces are also called the lateral forces). The 0 degrees in Figure 24 and the following figures points to the negative direction (Figure 10). It is clear that the mean values of the transverse and vertical forces of one blade are close to zero. That is the reason why the transverse and vertical forces of the whole propeller (five blades) are much smaller than the axial force. From Figure 24, the cross flows have more pronounced impact on the axial force than the lateral forces.

As for one blade, the exciting forces complete a periodic change when the propeller rotates one cycle because of the bilaterally symmetric ship wake. Consequently, the base frequency of the exciting forces of one blade is the shaft frequency of the propeller (9.5 Hz). The frequency analysis of the exciting bearing forces of one blade is performed and the results are shown in Figure 25. First, the three exciting forces of one blade are primarily composed of the shaft frequency (SF) and multiple SF components. Second, the amplitudes of SF of these two lateral forces are greater than that of the axial force. Third, the cross flows have little influence on the main harmonics of the lateral forces. It should be noted that the amplitude of the axial force is increased obviously if the cross flows are not considered.

The shafting structural vibration is caused by exciting bearing forces of the whole propeller. The predicted nondimensionalized axial force, , transverse force, , and vertical force, , of the whole propeller under ship’s wake with/without cross flows are shown in Figure 26. It is clear that the mean values of the transverse force and vertical force are much smaller than that of the axial force. Generally, these two lateral forces are not considered. However, if shafting bending vibration is the most concerned problem, then these two lateral forces should be studied. The influence of cross flows on the transverse force should cause attention. Ignoring the cross flows, the transverse force of the whole propeller is increased obviously.

As for the whole propeller, the exciting forces complete a periodic change when the propeller rotates 72 degrees (360/5) because the propeller has five blades. Consequently, the base frequency of the exciting forces of the whole propeller is the blade passing frequency ( Hz). The frequency analysis of the exciting forces of the whole propeller is performed and the results are shown in Figure 27. First, the three exciting forces of the whole propeller are primarily composed of the blade passing frequency (BPF) component. Second, it is clear that the cross flows have little influence on the axial force. Third, it should be noted that the cross flows have strong influence on the harmonics of the transverse force and vertical force. Ignoring the cross flows, the amplitudes of BPF of the transverse and vertical forces are increased significantly.

5. Conclusions
In this paper, numerical methods were applied for predicting propeller exciting forces under nonuniform ship’s wake. The influence of cross flows on the three-dimension exciting forces of one blade and the whole propeller was studied. Results show that the influence of cross flows in ship’s wake on the propeller exciting bearing forces is considerable. Some important conclusions were made as follows.
Numerical results, including ship resistances, wave patterns, and velocity distributions, predicted by DES and URANS are very close and agree well with the measured data. Circumferential distributions of the axial, tangential, and radial velocities of several profiles at the propeller plane were studied, and the axial velocity was validated by measured results. Numerical results show that the radial and tangential velocities are much smaller than the axial velocity. The mean wake coefficient and the amplitudes of the low harmonics of the predicted wake are very close to that of the measured wake. The predicted ship’s wake by DES contains larger amplitudes of the high harmonics than the measured data.
As for the exciting forces of one blade, comparatively, the cross flows have greater influence on the axial force than that of the other two lateral forces. When the cross flows are not considered, the amplitude of SF of the axial force increases significantly. It should be noted that, considering the cross flows or not, the amplitudes of SF of the transverse and vertical forces are even larger than that of the axial force.
As for the exciting forces of the whole propeller, first, the mean values of the transverse and vertical forces are much smaller than the axial force. Second, as exciting forces, the cross flows have strong influence on the transverse and vertical forces than the axial force. When the cross flows are not considered, the amplitudes of BPF of the transverse and vertical forces increase significantly. Last, the cross flows also have strong influence on the mean value of the transverse force of the whole propeller.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (no. 51479116 and no. 11272213).