1. Introduction

Cavitation is usually inevitable for marine propellers under working condition and the cavity generation and collapse and its unsteady characteristics can cause performance losses, structure damage and other problems. Many works have been done on such subjects for several centuries. Arndt (1981), Kerwin (1986), Kumar and Saini (2010) generalized the fundamental concepts and influence factors in their review papers. Kinnas (1996) summarized that the propeller cavitation exhibits several forms, including bubble cavitation, sheet cavitation, cloud cavitation, tip vortex cavitation and so on. The sheet cavity dynamics and tip vortex cavity dynamics are the main source of vibration and noise. In order to control or reduce the negative effects of the propeller cavitation, it is crucial to master the mechanism of the cavitation evolution and the physical characteristics of cavitating flow. Because of its complexity, many researchers studied these cavities separately. For example, Foeth (2008), Chen et al. (2016, 2017, 2019) and Ji et al. (2015) studied the sheet cavitation on various hydrofoils to master its structures evolution and unsteady characteristics. Arndt et al. (1991, 2015), Zhang et al. (2015) and Pennings et al. (2016) analyzed the behavior of the tip vortex cavitation and the cavity vibration models. For the purpose of predicting and improving the performance of the propellers, all kinds of propeller cavitation should be concerned. Therefore, the cavitating propeller in uniform or/and non-uniform wake flow gets more and more attention and is usually studied by experimental or numerical simulation means.

Pereira et al. (2002, 2004a) presented an experimental investigation on a cavitating propeller in a uniform inflow. They measured both the cavity extension and thickness, and established a complete mapping of the propeller hydrodynamic characteristics, including the cavitation pattern bucket diagram and pressure coefficient distributions on the blade suction side. Pereira et al. (2004b) also experimentally studied the cavitating propeller in non-uniformed inflow. They used an array of plates to model the wake downstream of a hull, took the harmonic analyses to identify the pressure and noise signals and found that the leading edge cavity contributed mainly to the fourth order and the tip vortex collapse acted in particular on the third harmonics. Korkut and Atlar (2012) carried out an experimental study to investigate the effects of foul release coating on the performance, cavitation and noise characteristics of a tanker model propeller. They found that cavitation development on blades was slightly reduced by the coating and the coating of blades reduced noise level under non-cavitating condition, but increased under developed cavitation condition. But they addressed that above conclusions required further investigation.

Salvatore et al. (2003) developed a theoretical and numerical methodology for marine propeller cavitation modelling at model scale, which was based on a nonlinear boundary integral formulation and boundary element method. They observed that the shape of the wake had a relevant impact on the cavity pattern prediction, even in the case of partial cavities that do not extend to the wake. And Salvatore et al. (2009) listed the results from the Rome 2008 Workshop on cavitating propeller. They presented seven sets of numerical results for INSEAN E779A propeller in uniform inflow and in a wake field, and compared them with a common experimental benchmark. The results show that in open water performance all the prediction errors for thrust and torque are within 5%, a potential flow BEM provided the most reliable results for the pressure fluctuations and LES with high gird-resolution could be used to access erosion risk, the cavity extents obtained by RANS were affected by turbulence and cavitation models. Bertetta et al. (2012) used potential panel methods to analyze the CPP propeller cavitation and noise optimization at different pitches, and such method is validated by the experiments at the towing tank and at the cavitation tunnel. Morgut et al. (2015) numerically simulated cavitating flow around a marine model scale propeller in non-uniform flow. They proved that three mass transfer models such as Zwart model, Singhal model and Kunz model could guarantee similar accuracy when the standard RANS were employed. Ji et al. (2012) employed the $ k-{\textit{ω}} $ SST turbulence model and a sliding mesh to numerically simulate the cavitation evolution around a propeller in non-uniform wake and verified that acceleration caused by the cavity volume changes was the main source of the pressure fluctuations resulted from propeller cavitation. Wu et al. (2018) used the $ k-{\textit{ω}} $ SST turbulence model with the turbulent viscosity correction to simulate the cavitating flow around a propeller in a non-uniform wake. They proved that cavity development is corresponding to the periodic large pressure fluctuation around the blade, with the dominant frequencies in accordance with the first blade passing frequency. Zhu et al. (2017) compared the tip vortex cavity shapes obtained by three different turbulent models based on Boussinesq hypothesis and proved that the turbulence model and advanced grids are important to predict tip vortex cavitation. Wang et al. (2018) simulated an oscillating propeller in cavitating flow and analyzed the effects of the non-uniform incoming flows on the propeller’s loads and the distribution of sheet cavitation and wake field.

All the studies above proved that both the experimental means and the numerical simulation method performed well on cavitating propeller problems. However, it should be noted that in most existing numerical researches, the turbulence models were based on turbulent-viscosity hypothesis, which can hardly capture the tip vortex cavitation. In order to overcome such problems, Yilmaz et al. (2019) employed the LES model to simulate the tip vortex cavitation on rotating propeller at the expense of very large mesh numbers. In the present work, the Reynolds-Averaged Naiver−Stokes (RANS) approach with Reynolds stress turbulence model (RSM) is used and its capacity on simulating the tip vortex cavitation is proved, which makes the propeller cavitation simulation proceed in a relatively acceptable mesh numbers for engineering application. The unsteady characteristics of the propeller cavitation in one uniform and two non-uniform wake flows are studied, separately. The common features and differentiating behaviors affected by the wake flows are analyzed.

2. Numerical simulation methodology

The cavitating propeller in the wake flow is a very complicated problem, which involves turbulence, phase transition and multiple coordinate systems. The mathematical models and numerical methods need special considerations.

2.1 Reynolds-averaged Navier−Stokes equations

Assumes that any physical variable is consistent on each phase component and the mixture of liquid and vapor is taken as a kind of one-fluid compressible medium. The governing equations for this problem are the instantaneous Navier−Stokes equations. In Reynolds averaging, the solution variables are decomposed into the mean and fluctuating components, and the mass and momentum conservation equations are:

where u_i is the mean velocity in the i direction, and term $- {\textit{ρ}} \overline {{u'_i}{u'_j}}$ is the Reynolds stresses. The density and viscosity are defined as a linear mixture of the liquid phase and vapor phase, and can be written as ${\textit{ρ}} = {{\textit{α}} _v}{{\textit{ρ}} _v} + \left({1 - {{\textit{α}} _v}} \right){{\textit{ρ}} _l}$ and ${\textit{μ}} = {{\textit{α}} _v}{{\textit{μ}} _v} + \left({1 - {{\textit{α}} _v}} \right){{\textit{μ}} _l}$ , respectively. Define ${{\textit{α}} _v}$ as the volume fraction of vapor phase and its transport equation can be written as:

where S is the source term that stands for the mass transfer between the liquid and vapor phase and its specific form is defined by the cavitation model.

2.2 Reynolds stress turbulence model

Since the common types of RANS turbulent models are based on Boussinesq hypothesis, which neglects the anisotropy of turbulence and makes these models perform poorly in propeller flows. Here, we adopt the Reynolds stress model (RSM), which includes several effects that are not easily handled by eddy-viscosity models. In RSM, model transport equations are solved for the individual Reynolds stresses and for the dissipation to close the Reynolds-averaged Navier−Stokes equations. The exact equations for the transport of the Reynolds stresses, $- {\textit{ρ}} \overline {{u'_i}{u'_j}}$ , can be written as follows:

Here, each term stands for special component, where C_ij is the convection, D_T,ij is the turbulent diffusion, D_L,ij is the Molecular diffusion, P_ij is the stress production, ${{\textit{φ}} _{i\!j}}$ is the pressure strain, ${{\textit{ε}} _{i\!j}}$ is the dissipation and F_ij is the production by system rotation. It should be noticed that the buoyancy production is neglected for the cavitating propeller problem considered in this paper. Among them, C_ij, D_L,ij, P_ij and F_ij can be solved directly, whereas, D_T,ij, ${{\textit{φ}} _{i\!j}}$ and ${{\textit{ε}} _{i\!j}}$ need extra models.

The turbulent diffusion term D_T,ij can be simplified as a scalar turbulent diffusivity,

where the turbulent viscosity ${{\textit{μ}} _{\rm{t}}} = {C_{\textit{μ}} }{\textit{ρ}} {{{k^2}} / {\textit{ε}} }$ , ${C_{\textit{μ}} }{\rm{ = }}0.09$ and ${{\textit{σ}} _{k1}} = 0.82$ .

The pressure strain term ${{\textit{φ}} _{ij}}$ are assumed to be linear and can be decomposed into three contributions,

where ${{\textit{φ}} _{i\!j,1}}$ is the slow pressure-strain term, ${{\textit{φ}} _{i\!j,2}}$ is the rapid pressure-strain term, and ${{\textit{φ}} _{i\!j,{\rm{w}}}}$ is the wall-reflection term. According to the work of Launder et al. (1975), Launder (1989) and Lien and Leschziner (1994), the modeled expressions are as follows:

where ${C_1} = 1.8$ , ${C_2} = 0.6$ , ${C'_1} = 0.5$ , ${C'_2} = 0.3$ , P_ij, F_ij and C_ij are defined in Eq. (4), n_k is the x_k component of the unit normal to the wall, d is the normal distance to the wall, and ${C_l}{\rm{ = }}\left.{{C_{\textit{μ}} ^{3/4}}} \right/ {\textit{κ}} $ , where ${\textit{κ}} {\rm{ = }}0.4187$ is the von Kármán constant.

For high-Reynolds-number flows, the dissipation term ${{\textit{ε}} _{i\!j}}$ is assumed to be local isotropy and can be modeled as (Pope, 2000):

Above expressions contain two unknown variables, turbulence kinetic energy k and dissipation rate ${\textit{ε}} $ , so two sets of equations are introduced to close the Reynolds stress turbulence model.

where ${{\textit{σ}} _k} = 1.0$ , ${{\textit{σ}} _{\textit{ε}} } = 1.0$ , ${C_{{\textit{ε}} 1}} = 1.44$ , and ${C_{{\textit{ε}} 2}} = 1.92$ .

2.3 Cavitation model

Mass transfer between the liquid and vapor phase includes both bubble formation (evaporation) and collapse (condensation). The source term S in Eq. (3) can be written as:

where R_e and R_c are the mass transfer terms generated by evaporation and condensation, respectively. Zwart et al. (2004) assumed that all the bubbles have the same size and proposed the Zwart−Gerber−Belamri model,

where p_v is the saturated pressure, ${{\textit{α}} _{{\rm{nuc}}}}$ is the nuclei volume fraction, ${R_{\rm{B}}}$ is the nuclei radius, ${F_{{\rm{vap}}}}$ is the evaporation coefficient and ${F_{{\rm{cond}}}}$ the condensation coefficient. In general, the model constants are ${{\textit{α}} _{\rm{nuc}}}{\rm{ = }}5 \times {10^{{\rm{ - }}4}}$ , ${R_{\rm{B}}}{\rm{ = }}1 \times {10^{{\rm{ - }}6}}$ , ${F_{{\rm{vap}}}}{\rm{ = }}50$ , and ${F_{{\rm{cond}}}}{\rm{ = }}0.01$ .

2.4 Numerical discretization method

The numerical simulation was carried out using the commercial CFD software ANSYS Fluent 14.0, which was based on the finite volume method. The SIMPLE algorithm was chosen for the calculation and the VOF model was taken to capture the liquid-vapor interface. The pressure, momentum and volume fraction terms were discretized by PRESTO! scheme, second order upwind scheme and Geo-reconstructed scheme, respectively. The body-fitted multi-block mesh strategy was adopted and the dynamic mesh model were used to simulate the rotating movement.

3. Computational conditions

A four-bladed INSEAN E779A model propeller was used for the research on the unsteady characteristics of the propeller cavitation in wake flows. This model propeller has a radium R=D/2=113.5 mm, a pitch-to-diameter of 1.1, and the hub has a diameter D_hub= 40 mm, as shown in Fig. 1. The computational domain is a cylinder with a diameter of 700 mm and a length of 1000 mm. For using the dynamic mesh, the entire domain is divided into two subdomains and the inside one is a cylinder with a diameter of 350 mm, as shown in Fig. 2.

Figure  1.  INSEAN E779A model propeller shape.

DownLoad: Full-Size Img PowerPoint

Figure  2.  Scheme of the computational domain.

DownLoad: Full-Size Img PowerPoint

In all our simulations, the advanced coefficient J=U/(nD)=0.71, and the propeller angular velocity n=36 rps. The cavitation number ${{\textit{σ}} _{\rm{n}}} = {{\left({{P_\infty } - {P_v}} \right)} / {\left[ {0.5{{\textit{ρ}} _l}{{\left({nD} \right)}^2}} \right]}}{\rm{ = }}1.763$ , where ${P_\infty }$ is the pressure in the far field and ${P_v}$ is the saturation pressure. In the uniform wake flow, the upstream inflow velocity U is 5.808 m/s in axial direction. In non-uniform wake flow, the average value of the inflow velocity is also 5.808 m/s in axial direction. Zero transversal velocity components along the inlet boundary are enforced in all cases.

For the numerical calculation, the entire domain is discretized into 6.67×10⁶ structured meshes. Fig. 3 exhibits the entire mesh on the computational domain and the local mesh on the propeller, separately. The time step size is set to 0.001 s and the convergence criterion requires that globally scaled residuals decreased to 10⁻⁴ for all equations to guarantee that the iteration steps in each time step are 100. The value of y⁺ on the four blades surfaces lies between 30 and 100, as shown in Fig. 4. Therefore, the standard wall function is adopted for the turbulent calculation.

Figure  3.  Computational meshes.

DownLoad: Full-Size Img PowerPoint

Figure  4.  y⁺ distribution on both sides of the propeller in uniform wake flow.

DownLoad: Full-Size Img PowerPoint

4. Verification of RSM turbulence model

Here, the benchmark of the INSEAN E779A model propeller in uniform wake flow is investigated by using the RSM turbulence model. The studied object and chosen parameters are the same in Salvatore’s experimental work (Salvatore et al., 2009). The cavitation shapes on the propeller blade in the uniform wake flow are given in Fig. 5, where the experimental result and numerical result, together with their comparison, are presented, respectively. The cavitation shape in our numerical result, defined as the value of vapor volume fraction of 0.5, is marked in light blue color in Fig. 5b. It can be seen that our numerical result agrees well with the experimental result in Salvatore’s work (Salvatore et al., 2009). The sheet cavitation is almost superimposed and the tip vortex cavitation is also coincident in shape and direction. However, the tip vortex cavitation in simulation is not as long as in experiment because the mesh at a distance from the blade is not fine enough to be captured.

Figure  5.  Cavitation shapes on the propeller blade in the uniform wake flow.

DownLoad: Full-Size Img PowerPoint

5. Results and discussion

Besides the uniform wake flow, two kinds of non-uniform wake flow conditions are studied here. The first one, referenced from Salvatore et al.’s (2009) work, is used to approximately simulate a propeller operating behind a single-screw hull with the inflow velocity of 6.22 m/s. Fig. 6a gives its axial velocity distribution, which is defined through wakefield measurements by LDV at a transversal plane at distance d = 0.26D upstream propeller disc. The second one is referred from Felli and Di Felice’s (2005) work where a twin-screw ship model is placed in a large circulating water channel and the propeller is under a hull wake condition, with the inflow velocity of 2.4 m/s. Fig. 6b describes the axial velocity distribution of the hull wake, which is measured by LDV without the propeller in the tunnel and the value is obtained at the position where the propeller is going to be placed. The main features of both non-uniform wake flows are absorbed in our numerical simulation, and especially for the latter one focusing on the properties that can be passed from the velocity inlet to the location of the propeller, so that some uneven distribution neither deviating from the mean velocity smaller than 8% nor far from the propeller region are neglected.

Figure  6.  Flow condition and the axial velocity distribution in non-homogeneous flow cited from (a) Salvatore et al.’s (2009) experiment and (b) Felli and Di Felice’s (2005) experiment.

DownLoad: Full-Size Img PowerPoint

The average value of the inflow velocity is 5.808 m/s in all cases. As the propagation of the wake flow is affected by the propeller rotation, the variation of the velocity distribution from velocity inlet to the propeller disc is concerned. Three slices of cross section perpendicular to the propeller hub are taken in the computational domain, as shown in Fig. 7. The first slice x = −0.3 m is at the position of velocity inlet and the third slice x = 0 m is at the middle section of the propeller. The second slice x =−0.2 m lies between these two slices. We denote the flow condition in the uniform wake flow and two kinds of non-uniform wake flows as UF, WF1 and WF2, respectively. The axial velocity distributions at three cross sections are given in Fig. 8. At x = −0.3 m, the velocity distribution is provided by the boundary conditions. In both non-uniform wake flows, the velocity at x = −0.2 m and x = 0 m is clearly influenced by the rotation of the propeller. Whereas, in the uniform wake flow the inlet velocity transports steadily to the location of x = −0.2 m and only get effected in the vicinity of the propeller at x = 0 m.

Figure  7.  Sketch for the three slices of the cross section.

DownLoad: Full-Size Img PowerPoint

Figure  8.  Axial velocity distribution at three cross sections in different wake flow (unit: m/s).

DownLoad: Full-Size Img PowerPoint

The cavity evolutions in a random period under three conditions are presented in Figs. 9−11. It can be clearly seen that the tip vortex cavities are at least partially captured in all cases. In UF, as shown in Fig. 9, the sheet cavity starting from the leading edge is rather thin and mainly stays in the vicinity of the blade tip during its development. Whereas, the root cavity appears and collapses intermittently in the rotating period. By comparing Fig. 10 with Fig. 9, it is obvious that the sheet cavity in WF1 covers more blade surface than in UF. In WF1, as shown in Fig. 10, at some time the cavity performs like an entirety, i.e. the sheet cavity and the root cavity are connected on the suction surface. It should be noted that the bubble cavity may also exist in such cavity entirety, but it is hard to tell through the RANS results based on our mesh resolution. Fig. 11 exhibits the cavity performance in WF2, which is similar to the situation in UF, but with a slightly thicker sheet cavity and a relatively more continuous root cavity. These differences in the three cases can be explained by the pressure variation, which determines the cavity generation, development and breakup. The cavity evolution exhibits strongly unsteady characteristics in all cases. Fig. 12 illustrates the pressure variation at three rotating moment for each wake flow case. The interval angle between two adjacent pictures is 32.4° for each case. It can be seen that the pressure in UF and WF2 experiences large variation in the entire surface and the increment or decrement is especially obvious in the areas around the propeller blades and sometimes in a high value level, which causes the sheet cavitation thinner and the tip-vortex cavitation shorter. Wheras, the pressure in WF1 is always within a relatively low value, compared with the other cases, therefore, the sheet cavitation remains relatively thick during the entire rotating period. The cavity volume variation with respect to time is given in Fig. 13. The cavity volume does not have an obvious characteristic of periodicity because of the various complicated cavity patterns.

Figure  9.  Cavity evolution in the uniform wake flow, and the interval angle between two adjacent pictures being 32.4°.

DownLoad: Full-Size Img PowerPoint

Figure  11.  Cavity evolution in the non-uniform wake flow of WF2, the interval angle between two adjacent pictures being 32.4°.

DownLoad: Full-Size Img PowerPoint

Figure  10.  Cavity evolution in the non-uniform wake flow of WF1, the interval angle between two adjacent pictures being 32.4°.

DownLoad: Full-Size Img PowerPoint

Figure  12.  Contour of pressure at the cross section of x=0, the color bar with the unit of kPa, and the interval angle between two adjacent pictures being 32.4° for each case.

DownLoad: Full-Size Img PowerPoint

Figure  13.  Volume variation of the cavity with respect to time.

DownLoad: Full-Size Img PowerPoint

The Q-criterion proposed by Hunt et al. (1988) is used to distinguish the vortex structure,

where ${{\textit{Ω}} _{i\!j}}$ is the vorticity tensor and ${e_{i\!j}}$ is the rate of strain tensor. The difference between the double contraction of these two tensors reflects the relative strength of the rotation and deformation of a fluid element. The region where Q>0 is defined as the vortex zone, and the rotation of the fluid plays a dominant role. Here, the iso-surface ofQ=10000 s⁻² is used to describe the vortical structures, which present a cylindrical spiral shape in each case, as shown in Fig. 14. By taking the flow characteristics in UF as a reference frame, the screw pitch of the spiral vortical structures is shorter in WF1 but slightly longer in WF2. It can be visually explained by the velocity distribution, which exhibits similar forms during a rotating period. The streamlines overlapped on the velocity contour at the cross section of x=0 are given in Fig. 15 at a random given time, when the angle of the propeller blades in the three wake flow cases are alike. In order to give a convenient analyses, the iso-surface of x=0 is divided into three zones as the propeller zone, the middle zone and the outer zone. The propeller zone is the area passed by the propeller and the velocity distribution is affected by the rotating motion, which can be observed from the shape of the streamlines, therefore, the velocity in the propeller zone is significantly larger than that in the other two zones and mainly marked by the red color in Fig. 15 for all three cases. The outer zones are confined by the blue dash lines and the velocity distribution is determined by the wake flow. The middle zone lies between them and the velocity distribution is influenced by the interaction of the rotating motion and the wake flow. Fig.15a shows that for the UF case, the streamlines are basically straight along the radial direction and the velocity distribution is almost uniform in the outer zone, whereas, the velocity obviously decreases at some areas surrounding the blades in the middle zone. Fig. 15b illustrates that for the WF1 case, the streamlines are highly curved and the velocity distribution is consistent with the coming flow in the outer zone, and the velocity in the middle region is totally different from the velocity in the outer zone, i.e., it is greater than the velocity in a small square region of the outer zone and smaller than the velocity in the other region of the outer zone, which produce a larger velocity gradient than that in the same area in UF case. This causes the spiral vortical structures shorter in WF1 than that in UF. Fig. 15c exhibits that for the WF2, the streamlines are somewhat curved and the velocity distribution exits slightly nonuniform property in the outer zone, and the velocity in the middle zone obviously decreases at a small region smaller than that in the UF case. Therefore, the flow characteristics in WF2 are very close to those in UF, the screw pitch of the spiral vortical structures in WF2 is only slightly longer than that in UF. It is worth noting that in WF1 the velocity in a relatively small region is smaller than that in the other region. However, in WF2 the contrary is the case, which may generate some components of the rate of strain tensor in the opposite signs and make the spiral vortical structures tight in WF1 but slightly loose in WF2.

Figure  14.  Iso-surface for constant value of Q-criterion (Q=10000 s⁻²) in three wake flows.

DownLoad: Full-Size Img PowerPoint

Figure  15.  Contour of velocity distribution with the streamlines (black solid lines) at the cross section of x=0, and blue dash lines being used to confine the outer region for analytical convenience.

DownLoad: Full-Size Img PowerPoint

6. Conclusion

A numerical method, based on the Reynolds-Averaged Naiver−Stokes (RANS) approach and Reynolds stress turbulence model (RSM), is used to study the propeller cavitation in the wake flow in the present study. The four-bladed INSEAN E779A model propeller rotating in one uniform and two nonuniform wake flows is analyzed here. The numerical results in the uniform flow shows a good agreement with previous experimental results on both the sheet cavitation and the tip vortex cavitation, which proves that the RSM possesses the capability of capturing the tip vortex cavitation. Two kinds of nonuniform wake flows are designed based on previous experimental researches and both of them have the same mean velocity as the uniform wake flow. The unsteady characteristics of the propeller cavitation are analyzed by comparing the results in one uniform and two nonuniform wake flows together. The results show that the effect of the wake flow on the propeller cavitation greatly depends on the distribution of the incoming velocity. The iso-surfaces of Q=10000 s⁻² are used to describe the vortical structures, which exhibits the spiral shape in the uniform wake flow, and the uneven velocity distribution causes the spiral vortical structures looser or tighter in two nonuniform wake flows.

Since in the present study the tip vortex cavitation is only partially captured, the optimizaton and refinement of the simulation mesh need to be paid more attention in the future. Meanwhile, the further research is going on the noise analysis because of the available pressure information on both the sheet cavitation and tip vortex cavitation, which are the main sources of the propeller noise.