Modeling of Unsteady Sheet Cavitation on Marine Propeller Blades

Unsteady sheet cavitation is very common on marine propulsor blades. The authors summarize a lifting-surface and a surface-panel model to solve for the unsteady cavitating flow around a propeller that is subject to nonaxisymmetric inflow. The time-dependent extent and thickness of the cavity were determined by using an iterative method. The cavity detachment was determined by applying the smooth detachment criterion in an iterative manner. A nonzeroradius developed vortex cavity model was utilized at the tip of the blade, and the trailing wake geometry was determined using a fully unsteady wake-alignment process. Comparisons of predictions by the two models and measurements from several experiments are given.

leading edge.The VLM with the leading-edge correction was incorporated into a code named PUF-3A by Kerwin and colleagues (1986).Vortex and source lattices were placed on the mean camber surface of the blade, and a robust arrangement of singularities and control-point spacings was employed to produce accurate results (Kinnas and Fine, 1989).The method was then extended to treat supercavitating propellers subjected to steady flow (Kudo and Kinnas, 1995).Recently, the method has been renamed MPUF-3A for its added ability to search for midchord cavitation (Kinnas et al., 1998).The latest version of MPUF-3A also includes the effect of hub, wake alignment in circumferentially averaged inflow with an arbitrary shaft inclination angle (Kinnas and Pyo, 1999), and of nonlinear thicknessloading coupling (Kinnas, 1992).However, the details of the flows at the blade's leading edge and tip cannot be captured accurately due to the breakdown of either the linear cavity theory or the thickness-loading coupling corrections.In addition, the current version of MPUF-3A does not include the effect of cavity sources in the thickness-loading coupling correction.
In Kinnas and Fine (1992) and Fine and Kinnas (1993), a loworder potential-based boundary element method (BEM) was introduced for the nonlinear analysis of three-dimensional flow around cavitating propellers subjected to nonaxisymmetric inflows.The method, named PROPCAV, was later extended to predict leading-edge and midchord partial cavitation on either the face or the back of the blades (Mueller and Kinnas, 1999).
PROPCAV inherently includes the effect of nonlinear thickness-loading coupling by discretizing the blade surface instead of the mean camber surface.Thus, PROPCAV requires more Central Processing Unit (CPU) time and memory but offers a better prediction of the flow details at the propeller's leading edge and tip than does MPUF-3A.In addition, the method provides a better foundation for concurrent research efforts in the modeling of developed tip-vortex cavitation and surfacepiercing propellers.
In this study, PROPCAV was further extended to treat simultaneous face and back cavitation on conventional and supercavitating propellers as well as fully unsteady wake alignment.

FIGURE 1
A propeller is subjected to a general inflow wake.The propeller's fixed coordinate system (x, y, z) and the ship's fixed coordinate system (x s , y s , z s ) are also shown (Kinnas et al., 2002).

THE BOUNDARY ELEMENT METHOD Formulation
The BEM formulation for flow around a cavitating propeller subjected to a nonaxisymmetric inflow is given in Kinnas and Fine (1992) and in Young and Kinnas (2001).
Consider a propeller that rotates at a constant angular velocity ω and is subject to a nonaxisymmetric inflow U w (x, r, θ).* The geometry and the coordinate systems are shown in Figure 1.The total inflow velocity is defined with respect to the propeller fixed coordinate system (x, y, z): where r = y 2 + z 2 and θ = tan −1 (z/y).
For inviscid and incompressible flow, the perturbation potential φ(x, y, z, t), which satisfies the Laplace equation (∇ 2 φ = 0), can be defined as follows: where q t (x, y, z, t) is the total flow velocity.The potential φ p at arbitrary point, p, on the body must satisfy Green's third identity: where the subscript p, q corresponds to the control and variable points in the integration.Three-dimensional Green's function, G( p; q), is defined as 1/R( p; q), and R( p; q) is the distance between points p and q. n q is the unit vector normal to the integration surface at the variable point, pointing into the fluid domain.φ is the potential jump across the wake surface, S W (t). S W B (t) is the combined wetted surface, which includes the wetted blade surface (S B ), the hub surface (S H ), and the tip vortex surface (S T ).S C (t) is the cavitating surface.

Boundary Conditions
• The flow tangency condition: the fluid flow is tangent to the propeller blades and cavity surfaces: • The dynamic boundary condition on the cavity surface: the pressure (P) everywhere on the cavity surface is constant and equal to the vapor pressure (P v ).It can be shown that this is equivalent to prescribing known values of φ on the cavity, which satisfies the following relation on S C (t) (Kinnas and Fine, 1992): where s and v are the local unit vectors defined at the each panel center in the chordwise and spanwise direction, respectively.σ n = P o −P v ρ 2 n 2 D 2 is the cavitation number.n and D are rotational frequency (revolutions per sec) and diameter of propeller, respectively.
• The kinematic boundary condition on cavity surface: the kinematic boundary condition renders the partial FIGURE 2 Definition of δ and points where the induced velocity is evaluated (Lee and Kinnas, 2003).differential equation for the cavity thickness, h (Kinnas and Fine, 1992): where V s ≡ ∂φ ∂s + U in • s and V n ≡ ∂φ ∂n + U in • n are the tangential and normal components of the total velocity vector, respectively.
• The blade sheet cavity closure condition: The cavity thickness at the end of partial or super cavities should be equal to zero.
• The Kutta condition: The velocity at the propeller trailing edge is finite, ∇φ < ∞.

Wake Alignment
A potential-based low-order panel method was used to compute the velocity field induced by the dipoles and sources of the system on the trailing-wake surface.The numerical instability in the roll-up region was avoided by introducing a tip vortex with a constant circular cross-section near the tip region of the wake sheet and by calculating the induced velocity at some slightly deviated (by a distance δ normal to the wake sheet) points from the control points, as shown in Figure 2.This treatment of the roll-up region is similar to that of Krasny (1987) and Ramsey (1996) and has been found to predict two-dimensional roll-up shapes that are quite similar to those of Krasny (Lee and Kinnas, 2003).
The velocity along the trajectory of the tip vortex core, V Tip , is evaluated by using the vector sum of the velocity vectors in the circumferential direction at each streamwise location along the tip vortex.The induced velocity on the trailing-wake panels can be computed by using Green's formula because the dipole and source strength on the propeller blade and hub panels and the dipole strengths of the wake panels are already known from the previous solution.Note that the dipole strengths on the wake surface along each strip are constant in steady flow, but those strengths are convected downstream with time in unsteady flow.The induced velocity on the wake surface is given by 4π u w (t) φ w (r q , θ q , t)∇ ∂G(p; q) ∂n q (t) dS. [7] Then the total velocity on the wake surface is determined by adding the total inflow velocities, U in (x, y, z, t), and the induced velocities, u w (x, y, z, t), which are computed by using Equation ( 7).
In order to find the aligned unsteady wake geometry that satisfies the force-free condition on its surface, the following numerical procedures are implemented at each steady, unsteady aligning, and fully unsteady step (Lee and Kinnas, 2003).
Steady Mode (t = 0) 1. Solve the steady boundary value problem (BVP) with the purely helical wake and without any modeling of the contraction and the roll-up at the blade tip. 2. Once the dipole strengths on the blades and the assumed tip vortex cavity surface are known from the BVP solution, calculate the induced velocity by applying Equation ( 7) at the displaced control points.
3. Compute the mean velocity at the center of the tip's vortex cavity, and interpolate the total velocities on the wake surface, from those at the control points to those at the panel edge points.4. Find the new coordinates of the wake panels by aligning with the total local velocities by using the streamline equation.
where U x , U y , U z are the x, y, z-direction total local velocities.
The new coordinate at (n + 1)th strip is determined by the following equation: where X n = (x, y, z) n , and δθ is the angular increment of trailing wake sheet.5. Repeat solving BVP and aligning the wake geometry with the updated new wake geometries until the wake geometries converge.6. Save the wake geometry and dipole strengths on the blades (φ(x, y, z, t = 0)) and wake panels ( φ(x, y, z, t = 0)) for the unsteady wake aligning process.These steady results are the initial values for the unsteady problem, described next.
Unsteady Aligning Mode (t > 0) 1.Initially, set the wake geometries of key and other blades to be the same as those in the steady mode.2. Solve the BVP (unsteady) with the aligned wake from the steady mode.In the unsteady mode, BVP is solved only for  the potential of the key blade and the tip vortex cavity, while the potential of other blades and the potential jump of other blade wakes are assumed to be known and equal to the values on the key blade when it was located where each other blade is at the current step.3. Compute the induced velocity on the control points of the key blade wake, and align the key blade wake geometry.4. Solve the BVP again with the aligned key blade wake and the same wakes of other blades as in Equation [2], and determine the dipole strengths of key blade panels.5. Save φ(t), φ(t), and the aligned key wake geometry.6. Move to the next time step (t + 1).Update the wake geometries, φ(t + 1) and φ(t + 1), of the other blades from the previously saved data.7. Repeat the unsteady run from Equations ( 2) to (6) until the wake geometries converge.

Fully Unsteady Mode
This mode does not perform wake alignment but uses the aligned wake, as predicted in the previous mode.
1. Update the wake geometries of the key and other blades corresponding to the time step t from the results of the unsteady aligning mode run.2. Update the φ(t) and φ(t) of the other blades and wakes at the corresponding time step.3. Repeat solving the BVP by updating φ(t) and φ(t) until the last revolution.

Validation
A three-dimensional elliptic hydrofoil is first considered to validate the numerical method of predicting the wake's roll-up and contraction.The cross-section of the wing has an NACA66-415 shape with an a = 0.8 mean camber line.The maximum thickness-to-chord ratio, (t/c) max , is 15%; the aspect ratio is AR = 3.0; and the angle of attack is 10 • .Figure 3 shows the converged trailing wake sheet behind an elliptic wing, where the contraction and the three-dimensional roll-up of the trailing wake can be seen very clearly.
In Figure 4, the tip vortex cavity trajectory computed by the present method is compared with that measured in the experiment by Arndt and colleagues (1991).The thick line of experimental measurements indicates the extent of variation in the trajectory for the different physical parameters, such as angle of attack, Reynolds number, and cavitation number.Note that in the experiment it was observed that the trajectory did not depend on the cavitation number, so the tip vortex trajectory under noncavitating conditions can also be used under cavitating conditions.The tip vortex trajectory produced by the present method (this trajectory is obtained from noncavitating solution) agrees well with that measured in the experiment.
The fully unsteady wake alignment scheme on the propeller is validated by comparing the predicted forces and moments with FIGURE 6 (A) Initial cavity shape for a three-dimensional hydrofoil section with detachment locations based on the wetted pressure distributions.Also shown are the fully wetted pressure distributions.(B) The converged cavity shape and corresponding pressure distributions (Kinnas et al., 2002).the experimental measurements and those predicted by using the vortex lattice method, MPUF-3A.Boswell and colleagues (1984) performed experiments using the DTMB 4661 propeller to measure the forces and moments under an inclined inflow condition.In MPUF-3A, the wake sheet is aligned by using the circumferentially averaged inflow and is adjusted to include the effect of shaft inclination (Kinnas and Pyo, 1999).Figure 5 shows the amplitude of the first harmonic of the forces acting on one blade of the DTMB 4661 propeller, in which the inclination angle, α = 20 • ; the advance ratio, J s = 1.0; and the Froude number, F n = 4.0, are considered.The forces and moments predicted by the present method compare well with those measured in the experiments, whereas the MPUF-3A predicts fewer forces than are measured.

Face or Back Cavitation with Searched Detachment
The search for face or back cavitation is necessary because it is common for propellers to be subjected to off-design

FIGURE 7
Predicted and measured thrust (K T ) and torque (K Q ) coefficients as a function of cavitating number (σ v ) and advance ratio (J s ) for DTMB4382 propeller (Kinnas et al., 2002).
conditions.Propellers are often designed to produce a certain mean thrust.However, part or all of the blade may experience smaller loadings at certain angular positions due to nonaxisymmetric inflow.As a result, alternating or simultaneous face and

Numerical Implementation
PROPCAV searches for the cavity detachments on both sides of the blade via an iterative algorithm.First, the initial detachment lines at each time step (or blade angle) are obtained based on the fully wetted pressure distributions.The detachment lines are then adjusted iteratively at every revolution until the Villat-Brillouin smooth detachment criterion is satisfied: 1.The cavity has nonnegative thickness at its leading edge, and 2. The pressure on the wetted portion of the blade upstream of the cavity should be greater than the vapor pressure.
An example of the initial cavity shape on a three-dimensional hydrofoil section with the detachment location obtained based on the wetted pressure distribution is shown in Figure 6A.Notice that the resulting cavity has negative thickness at the leading edge due to the incorrect guess concerning the location of the cavity detachment location.Also notice the considerable underprediction of the extent and volume of the cavities, especially on the face side.The converged cavity shape and the corresponding cavitating pressure obtained by using the detachment search algorithm are shown in Figure 6B.Notice that the smooth detachment criterion is satisfied because the cavity thickness is nonnegative, and the pressure everywhere on the wetted blade surface is above the vapor pressure.

FIGURE 10
Treatment of nonzero trailing-edge sections in fully wetted, partially cavitating, and supercavitating conditions.

Predicted versus Measured Forces
In order to validate the present method, the predicted forces are compared with those measured in the experiment.Boswell (1971) performed cavitation tests on a DTMB 4382 propeller in a 24-in cavitation tunnel at the Naval Ship Research and Development Center (NSRDC) to determine the thrust breakdown due to cavitation.The predicted thrust and torque coefficients as a function of advance ratio and cavitation number are shown in Figure 7.The predicted cavitating and fully wetted forces agree

FIGURE 12
The geometry of a DTMB5168 propeller and a comparison of the thrust and torque coefficient predicted by PROPCAV with those measured (Young and Kinnas, 2001).
well with those measured in experiments.It should be noted that the algorithm for cavity detachment had to be altered for lower cavitation numbers so that very thin cavities are excluded, as described by Kinnas and colleagues (2002).

Sample Case
An example of simultaneous face and back cavitation for propeller MW1 is shown in Figure 8.The propeller geometry, given in Young and Kinnas (2001), is based on a design by Michigan

FIGURE 13
A comparison of the predicted shaft thrust and torque harmonics based on experiment, PROPCAV, and MPUF-3A for a DTMB4119 propeller.Also shown are the propeller's geometry and inflow wake (Young and Kinnas, 2001).
Wheel Corporation (Grand Rapids, MI).The flow conditions were as follows: J = 1.2, σ = 0.8, F r = 25, and the inclined inflow was at 3 • .Notice that for this propeller there is midchord supercavitation on the suction side of the blade, and there is leading partial cavitation as well as midchord supercavitation on the pressure side of the blade.The unsteady cavitating pressure contours for propeller MW1 are shown in Figure 9.

Treatment of Nonzero Trailing Edge Blade Sections
Supercavitating propellers are often believed to be the most fuel-efficient propulsive devices for high-speed vessels.However, they are difficult to model because of the unknown size and pressure in the separated region behind the thick blade's trailing edge.In the BEM, the pressure in the separated region is assumed to be constant (as suggested by measurements) and to be equal to the vapor pressure.Thus, the size and extent of the separated region can be determined in the framework of a cavity problem.For a given propeller geometry, an initial guess about the separated region boundary is assumed; then the shape of the separated region and the cavities are solved simultaneously in an iterative manner until both the kinematic and dynamic boundary conditions are satisfied on all surfaces.The treatment of nonzero trailing-edge sections in fully wetted, partially cavitating, and supercavitating conditions is depicted in Figure 10.
An example of the predicted cavity shape and cavitating pressures for supercavitating propeller M.P. No. 345 (Ship
Research Institute, Tokyo, Japan) is shown in Figure 11.(A comparison of numerical predictions with experimental measurements of a wide range of flow conditions is shown in Figure 16.)It is worth noting that at J s = 1.3, there is substantial midchord detachment.Figure 11 indicates that the detachment search criterion in PROPCAV is satisfied because the cavity thickness is nonnegative, and the pressures everywhere on the wetted blade surfaces are above the vapor pressure.

VALIDATION BY EXPERIMENT
In order to thoroughly validate PROPCAV and MPUF-3A, four different sets of experiments were carried out.

Propeller DTMB5168
Figure 12 shows a comparison between measured thrust and torque coefficients determined by experiment and predictions by PROPCAV for propeller DTMB5168 in fully wetted, uniform inflow.The geometry of the propeller is also shown in Figure 12.Notice that PROPCAV yields quite accurate force predictions for a wide range of advance ratios.

Propeller DTMB4119
Figure 13 shows a comparison of unsteady thrust and torque coefficients obtained by experiment, by PROPCAV, and by MPUF-3A for a DTMB4119 propeller.The propeller is sub-jected to a nonaxisymmetric three-cycle wake (Jessup, 1990; also shown in Figure 13) in fully wetted flow.As shown in Figure 13, both numerical codes did well in predicting the unsteady blade-force harmonics.

Propeller DTMB4148
The test geometry for the third set of experiments is the DTMB4148 propeller, as shown in Figure 14.The propeller was

FIGURE 16
Comparison of the predicted and measured K T , K Q , and η for various advance ratios in an SRI propeller (Young and Kinnas, 2002).
subjected to a screen-generated nonaxisymmetric inflow inside a cavitation tunnel (Mishima et al., 1995) under the following conditions: J s = 0.9087, F r = 9.159, and σ n = 2.576.The inflow wake used in PROPCAV, which is shown in Figures 15A, B, and C corresponds to the wake described by Mishima and colleagues (1995).The effects of the tunnel walls and vortical inflowpropeller interactions (a nonaxisymmetric "effective" wake) are accounted for by using the method of Choi (2000) and Kinnas and colleagues (2000).The equivalent J s , 0.957, for unbounded flow is obtained by matching the fully wetted thrust coefficient, K T , with the measured K T , 0.0993, from the experiment.The predicted cavity shapes using the PSF2-type alignment (Greeley and Kerwin, 1982) without the tip vortex model are shown in Figure 15B.The predicted cavity shapes using the fully unsteady wake alignment with the tip vortex model are shown in Figure 15C.Although the cavity shapes predicted by both numerical models agree well with those of experimental observations, the former has convergence problems at the blade tip because of the lack of tip vortex modeling.

Propeller 345SRI
To validate the treatment of supercavitating propellers in PROPCAV, predicted force coefficients were compared with experimental measurements (Matsuda et al., 1994) of a supercavitating propeller.The test geometry is M. P. No. 345SRI, which is designed using SSPA charts under the following conditions: J = 1.10, σ v = 0.40, and K T = 0.160.As shown in Figure 16, the predicted thrust (K T ), torque (K Q ), and efficiency (η p ) compared well with measurements made in experiments.

CONCLUSIONS
A boundary element method and a vortex-lattice method for the prediction of sheet cavitation on propellers were presented.The BEM is able to treat complex types of cavitation patterns on the back and face of conventional and supercavitating blades, as well as unsteady wake alignment with a tip vortex model.The effects of viscosity can also be included via a viscous/inviscid interactive approach, as described by Kinnas and colleagues (1994) and by Brewer and Kinnas (1997).The numerical prediction by both methods compares well with experimental measurements.

FIGURE 4
FIGURE 4Comparison of the trajectory of the tip's vortex core with that of the experiment for the elliptic wing: AR = 3.0, (t/c) max = 0.15, and α = 10 • (Lee andKinnas, 2003).

FIGURE 5
FIGURE 5Comparison of the predicted first harmonic of the forces and moments acting on one blade of the DTMB4661 propeller:inclination angle, α = 20 • , J s = 10, and F n = 4.0(Kinnas et al., 2002;Lee and Kinnas, 2003).