On the fluid structure interaction of a marine cycloidal propeller

Applied Ocean Research 64 (2017) 105–127

Contents lists available at ScienceDirect

Applied Ocean Research journal homepage: www.elsevier.com/locate/apor

On the fluid structure interaction of a marine cycloidal propeller J.J. Prabhu, V. Nagarajan ∗ , M.R. Sunny, O.P. Sha Indian Institute of Technology, Kharagpur, India

a r t i c l e

i n f o

Article history: Received 2 September 2016 Received in revised form 25 December 2016 Accepted 28 January 2017 Keywords: Cycloidal propeller Unsteady load Vibration Fluid structure coupling

a b s t r a c t Marine cycloidal propulsion system is efficient in maneuvering ships like tugs, ferries, etc. It is capable of vectoring thrust in all direction in a horizontal plane. When used in pair, the system enables a vessel to perform maneuvers like moving sideways, perform rotation about a point, i.e. turning diameter of its own length, etc. In this system, the propeller blades have to change their angle of attack at different angular position of the disc. Due to this reason, the inflow velocity vector to propeller blades changes continuously. The propeller blade oscillates about a vertical axis passing through its body and at the same time rotates about a point. Superposed on these motions is the dynamics of the ship on which the propulsion system is installed. This results in a formidable and challenging hydrodynamics problem. Each of the propeller blade sections could be considered as an aerofoil operating in combined heave and pitch oscillation mode. Due to the constantly varying inflow velocity, the hydrodynamic flow is unsteady. The unsteady hydrodynamic flow is simulated by incorporating the effect of shed vortices at different time instant behind the trailing edge. Due to the kinematics of the problem, the blade is subjected to higher structural deformation and vibration load. The structural deformation and vibration when coupled with the hydrodynamic loading add another level of complexity to the problem. In this paper, the variation of hydrodynamic load on the propeller blade due to steady and unsteady flow is compared. We also model the structural dynamics of the blade and study its effect on the hydrodynamic loading. Finally, we couple the structural dynamics with hydrodynamics loading and study its influence on the propeller blade for different operating regimes. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction In marine cycloidal propulsion system, the propulsion and maneuvering units are combined together as a single unit. Therefore, no separate rudder is required to maneuver the ship. Marine Cycloidal Propellers (MCP) is fitted in pairs. Each unit consists of five to six blades. The location of the two propellers depends upon the type of ship. In ferries, the propellers may be fitted one forward and one aft, which would allow the vessel to move sideways or to turn within a circle of its own length. In a pusher tug, they may be fitted side by side at the stern, whereas in a tractor tug, the propellers may be fitted side by side at the bow. MCP has larger swept area. This makes the rotational speed of MCP about 75% less than the speed of the normal conventional propeller of comparable size and power [17]. The other benefit is reduced thrust loading which delays cavitation inception. The blades of MCP are arranged vertically on the periphery of a rotating disc. Therefore, the inflow velocity remains same throughout the length of the blade. Therefore, no twist is

∗ Corresponding author. E-mail address: vishwanath [email protected] (V. Nagarajan). http://dx.doi.org/10.1016/j.apor.2017.01.019 0141-1187/© 2017 Elsevier Ltd. All rights reserved.

required on blades, unlike the conventional marine screw propellers. However, the inflow velocity to the blades changes at the different angular position of the disc. This makes the flow unsteady. The angle of attack of the blade has to be changed continuously, due to the continuous variation in the inflow velocity. Therefore, the thrust generated by each blade varies with time. Due to this reason, the load variations on each one of the cycloidal propeller blades is high. This variation increases during maneuvering when the vessel has motions in different degrees of freedom. This time varying load can lead to vibration and fatigue failure of the propeller blade and its components. The varying inflow velocity vector also gives an opportunity to optimize the hydrodynamic performance of the propeller blade. The operating variables are very high for marine cycloidal propeller. It is difficult to achieve the above objectives by carrying out experiments. A huge number of experiments will be required for this purpose. To achieve the above objectives a good fluid structure coupled model for the cycloidal propeller system is required. Such a model will be helpful in the design and development of this propulsion system. The arrangement of a typical cycloidal propeller on a ship is shown in Fig. 1. Allan and Molyneux [1], presented the hydrodynamic performance of three different combinations of hull shape, appendages

106

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

Nomenclature A CD CL CM CP c D d d g k L N N.F n

Exposed area of blade (m2 ) Coefficient of drag Coefficient of lift Coefficient of moment Coefficient of pressure Chord length of blade (m) Drag force (kN) Depth from free surface to the propeller surface (m) Change in shed vortices Acceleration due to gravity (m/s2 ) Reduced frequency Lift force (kN) Number of panels Normal chord force on blade (kN) RPS of propeller

P

Distance vector from blade stock to midpoint of panels Atmospheric pressure (101.325 kN/m2 ) Vapour pressure (1.704 kN/m2 ) Hydrodynamic torque (kN m) Source strength per unit length (m/s-1 ) Radius of the propeller disk (m) Polar coordinates of the field point Length of panel (m) Advance velocity of ship (m/s) Onset flow to the blade (m/s) Field point Angle of attack to the blade (deg) Vortex strength per unit length (m/s-1 ) Angle of line joining the vortex shed point and panel Angular position of disc (deg) Advance coefficient of ship Load coefficient of propeller Density of sea water (1025 kg/m-3 ) Angle of line joining the stock point and panel Potential due to the pitching motion of blade Potential due to the source Potential due to the vortices Potential due to shed vortices Potential due to the uniform onset flow Inflow angle to the blade (deg) Cavitation number Pitching rate of propeller blade (rad/s) Rotational speed of propeller (rad/s)

→

PA PV Q q(s) R r,  s VA VR x, y ␣ (s)     w

B

s

v



∞ ˝ ωf ωp

and propulsion system, and discussed the effectiveness of alternative design solutions for escort tugs with specified steering force requirements. They found there is no significant difference between the hydrodynamic forces of the combined hull and propulsion systems tested for the three tug concepts tested as part of the project. The specified steering force (for example, 150 t) can be achieved by all of the concepts, provided that the hull is designed to provide the appropriate level of stability. Barrett et al. [2] experimentally computed the force and power on robotic fish. They demonstrated that the power required to propel an actively swimming, streamlined, fish-like body is significantly smaller than the power needed to tow the body straight and rigid at the same speed. The unsteady motion of a body offers possibilities for more efficient propulsion as compared to conventional steady propellers. It has been shown that unsteady motion of airfoils result in higher lift coefficient and produce propulsive thrust very efficiently. Haberman

Fig. 1. Typical set up of MCP on a ship.

and Harley [7] used Taniguchi method to compute the performance characteristics of vertical axis propellers having cycloidal blade motion and semi-elliptic blades. Taniguchi method is based on the assumption that only the induced velocities due to the trailing vortex in the direction of propeller advance contributes to the thrust and torque of the propeller. The vortex remains constant over the length of blade. The value of the induced velocity is obtained from momentum considerations with modifications based on experimental performance of a six-bladed vertical axis propeller. It is found that the thrust and torque of each blade and the maximum efficiency of the propeller decreases with increase in number of blades of same dimensions. Jurgens and Heinke [10] studied the cavitation behaviour of heavily loaded Voith Schneider Propeller (VSP) blades under bollard pull condition. Experiments were carried out with different profile shapes of the VSP blades. A new blade profile (P9659) was developed using a numerical optimisation strategy based on CFD calculation. The blade profile showed improved cavitation behaviour. Jurgens and Palm [11] studied the influence of the Voith Schneider Propeller (VSP) on stern slamming conditions, the roll damping capabilities and the impact of air ventilation. Investigations were carried out using experiments and computational fluid dynamics and compared with the results of azimuth thrusters. It is found that the VSP has a reduced impact of pressure loads due to stern slamming. The VSP uses positively the slamming behaviour of a vessel because of its vertical rotating axis. Additionally the VSP can be used to reduce the roll motion of a vessel. When used in this configuration it is called as Voith Roll Stabilization System. The VSP provides effective functionality for roll damping besides propulsion and steering. Palm et al. [12] compared the ventilation characteristics of azimuth thrusters with Voith Schneider propeller. Experiments were carried out and results supported with CFD numerical simulation. They found that the cycloidal propeller is less prone to ventilation than the azimuth thrusters. Sfakiotakis et al. [13] presented the fish swimming types and specific swimming modes. The study is made on the propulsor and the type of movements (oscillatory or undulatory) employed for thrust generation. The locomotion by angelfish is similar to rowing action consisting of two phases, the power and recovery stroke. In power stroke, the fins move perpendicular to the body at a high attack angle and with a velocity greater than the overall swimming speed. During recovery stroke, the fins are feathered to reduce resistance and brought forward. In MCP also, the blades are making the power and recovery stroke as the blades moves backward and forward respectively. In angel fish locomotion no thrust is generated in recovery stroke but in MCP a nominal thrust is generated. Stefan et al. [14] carried out a study to demonstrate the capability of a cycloidal propeller to use unsteady dynamic lift for operation. The investigation was carried out at Reynolds numbers (Re = 103 ) based on the chord length. The thrust production and energy extraction

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

107

Fig. 2. Coordinate system describing ship and MCP.

are achieved by using the unsteady dynamic lift formulation. The pitching mechanism of the cycloidal propeller blades appears to be a good example for computing the thrust based on unsteady hydrodynamics. Waclawek and Molyneux [16], conducted planar motion mechanism (PMM) captive model experiments for a tug fitted with MCP system. The influence of inflow angle on the propeller performance is studied for a range of yaw angle and thruster angle. The braking and steering forces are recorded. It is observed that steering force is independent of the inflow angle up to inflow angle of 40◦ . The experiment describes the development of a hybrid modeling process that combines captive model experiments for obtaining hydrodynamic coefficients of the tug with numerical simulation of the complete system. In this paper, we develop a numerical model to determine the time varying hydrodynamic load distribution around a cycloidal propeller blade in operation. Potential flow approximation is made. The blade thrust and torque is computed and compared with experiment results. The numerical model can compute both the disc torque and the individual blade torque. We then model the structural dynamics of the blade and couple it with the hydrodynamic loading. The influence of higher structural load on propulsive thrust and torque is demonstrated.

sient motion. The marine cycloidal propeller may be driven either by an electric motor or by diesel engine mounted above the disc. The engine shaft and propeller shaft are usually connected through gear mechanism. Usually, vibration dampers are provided on the shaft to reduce the influence of variation in hydrodynamic load on the prime mover [17]. The eccentricity point is controlled by two hydraulic rams at right angles to each other so that a variation in the magnitude and direction of the propeller thrust can be obtained. The blades of a marine cycloidal propeller are usually made of stainless steel to withstand cavitation erosion and corrosion. For a given propeller thrust there is an optimum combination of eccentricity and the angular velocity that gives the minimum delivered power. However, it is usual to run the propeller at a constant rpm and control the magnitude of thrust by controlling only the eccentricity.

2. Operating methology

The angle ‘’ for port and starboard propeller unit can be computed from Eq. (2).

Each blade of MCP makes one pitching oscillation about its own axis as the disc makes one rotation. The location of MCP on the ship and its coordinate system is shown in Fig. 2. The path of each blade of MCP is some form of cycloid – an epicycloid if VA < ωR, a cycloid if VA = ωR, and a trochoid if VA > ωR. The magnitude of the thrust depends on the ship motion dynamics and the revolution of the disc. The thrust direction can be controlled by an imaginary eccentricity point on the disc, which is the distance from the center of the disc. The eccentricity can be adjusted to any value between zero and ±R. The magnitude of the resultant thrust of the propeller can be controlled by varying the eccentricity distance. While the direction of the thrust can be controlled by varying the angle of the eccentricity point with the propeller disc axis. By changing the eccentricity point, the angles of the blades can be changed with respect to the longitudinal axis of ship. The eccentricity point is fixed for a particular magnitude and direction of thrust. Each one of the propeller blade is connected to a mechanical gear link mechanism. The hydraulic gear arrangement ensures that the blade chord remains perpendicular to the line drawn from the blade to the eccentricity point. Due to faster time response of the hydraulic system, the blade orientation changes much faster than the ship dynamics. It must be noted that blade thrust and torque depends on the inflow velocity to the blade. Therefore, even if the eccentricity point is changed, there will be some time before the thrust magnitude and the direction is stabilized. The thrust magnitude and direction will stabilize only after the body dynamics is stabilized. The intermediate stage where the blade thrust and torque is not stabilized is called tran-

Fig. 3. Flow dynamics and force vector acting on a propeller blade of MCP system.

→

The flow dynamics and force vectors are shown in Fig. 3. Vector CA in Fig. 2 is computed as shown in Eq. (1).



→







CA = R cos P ˆi + R sin P ˆj



S





=

P

ωS



−ωP



(1)

t

(2)



S in Eq. (2) is for the first blade. For the other blades, Angle P there will be a phase angle difference depending on the location of the blade on the disc. The 3-DoF (surge (u), sway (v) and yaw (r)) dynamics at ship coordinate center can be written in the vector form as Eq. (3). →

→

VO = uˆi + vˆj, r k

(3)

Velocity at blade position ‘A’ due to ship dynamics is given by Eq. (4) →

→



→

→

VA = VO + r k × OA



⎫ ⎬







(4)



⎭ OA = xP + R cos P i + yP + R sin P j →

Substituting Eq. (3) in Eq. (4), we get Eq. (5). →

VA = uˆi + vˆj + r kˆ ×









xP + R cos P i + yP + R sin P j

(5)

Simplifying Eq. (5), we get Eq. (6), the velocity vector at blade position ‘A’ due to ship dynamics. →





 





 

VA = u − yP + R sin  r ˆi + v + xP + R cos  r ˆj

(6)

108

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

Fig. 4. LHS: Details of the face and back of the propeller blade at different blade orbital angle for a particular eccentricity ratio. RHS: Illustration showing that the blade performs combined heave and pitch motions superimposed on steady inflow speed.

It is noted that the blade has an angular rotation about the point C (ωP ). The velocity at A due to angular rotation is obtained as per Eq. (7).



→

  

 

VAA = ωP k x R cos i + R sin  j

(7)





VAA = (ωP R) − sin i + cos j

(8)

→

The resultant velocity VR (only due to ship and disc dynamics) and its angle of inclination ‘ ’ with the ship’s longitudinal axis can be computed from Eq. (9). →

→

→

V R = VA + VAA



→







V R = u − yP + R sin  r − ωP R sin  ˆi









⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬

⎪ ⎪ ⎪ ⎪ ⎪     ⎪ ⎪ ␯ + xP + R cos  r + ωP R cos  ⎪ ⎪ ⎪    ⎪ =  ⎭

+ + xP + R cos  r + ωP R cos  ˆj

tan

(9)

u − yP + R sin  r − ωP R sin 

The detailed operating mechanics of MCP propellers is defined in Ghose and Gokarn [5]. An eccentricity point with coordinates (e2 , e1 ) with respect to disc center (‘c’) is defined. If we draw a line from the eccentricity point to the blade stock position, the blade chord will be perpendicular to this line. This will be the situation at any position on the circumference of the disc. The angle of the line drawn from eccentricity point to the blade stock position with respect to the longitudinal axis of the ship is given by Eq. (10). tan

=

R sin  − e 1 R cos  − e2

(10)

The propeller blade angle of attack is given by Eq. (11). ˛=

− (90 + )

(11)

The lift and drag of the blade can be obtained by Eq. (12).



  L

D

=

CL 2CD

TY = L sin



(12)

The factor 2 with CD in Eq. (12) is to account for the wetted surface on both the sides of the propeller blade. The lift is normal and drag is tangential to the effective inflow velocity V R . We are













− 90 + D cos 180 −

 

− 90 − D sin 180 −





− 90



+ D cos  −

⎫ ⎪ ⎬ 

 ⎪ ⎭

(13)

− 90

It must be noted that Eq. (13) gives the thrust and torque for a single blade on the disc. The total thrust acting on the ship will be due to ‘Z’ number of blades each on the port and starboard disc. Total torque is individually computed on each disc due to ‘Z’ number of blades. In this paper, we add the contribution of thrust and torque from each blade to get the total propeller thrust and torque. However, it must be noted that the blades do mutually interact with each other and it is well known as propeller blade in cascade. The influence of mutual interaction will be investigated in our future work. It may be noted that torque on the port and starboard disc will be in opposite direction (due to their different direction of rotation). It may be noted that the magnitudes of TX , |TY | and |Q | are same for port and starboard blades, only when the ship has straight line motion. When ship has three DoF motions, their magnitudes will be different. For cycloidal propellers, when the forces due to all the blades are summed up the resultant magnitude of TY though small remains non-zero. The magnitude of resultant TY is higher for odd number of blades as compared to even number of blades. This is because, with even number of blades Y direction force of opposite blades on the disc tends to cancel out each other. Due to this reason, cycloid propellers are installed in pairs. Each one of the unit rotates in opposite direction so that the TY due to the individual units cancel out and no net sway force acts on the ship during forward motion. Else, the propulsive efficiency will significantly reduce due to effort expended in maneuvering the ship over a straight course. The thrust direction is given by Eq. (14).



ˇ = tan−1 e2 /e1



(14)

We are also interested in computing the magnitude of the pitching angle. For this purpose, we need to know the angle that the blade chord makes with the tangent drawn to the disc at any time instant. We define this angle as ‘ı’ and compute it as shown in Eq. (15). ı= −



1 w AVR2 2

TX = L cos

Q/R = L sin  −

Eq. (7) can be simplified as Eq. (8). →

interested in the thrust, torque and transverse forces that will be acting on the ship. They are estimated as shown in Eq. (13).

(15)

It can be observed from previous equations, that angles ‘ ’ and ‘’ vary disproportionately with time. Due to this reason, pitch angle (ı) keeps varying (or oscillating) with the rotation of the disc. Therefore, angle of attack (␣) (Eq. (11)) also keeps changing from positive to negative values. Hereafter, we will discuss about the straight running motion. When blade orbit angle () is between 270◦ ∼ 0◦ ∼ 90◦ ,

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

109

resolved as per ship coordinate system (see Eq. (13)). For computing the deflection and stresses on the individual blade, the lift, drag and centrifugal forces need to be resolved as per the blade coordinate system. The normal and tangential force acting on the blade is given by Eq. (19). FN = L sin(90 − ˛) + D sin ˛ + C.F sin( FT = L cos(90 − ˛) − D cos ˛ + C.F cos(

Fig. 5. Details of free vortex shed from the trailing edge of a single propeller blade of MCP system.

angle of attack (␣) is positive. However, when blade orbit angle () is between 90◦ ∼ 180◦ ∼ 270◦ , angle of attack (␣) is negative. Even though, angle of attack (␣) changes sign, the lift force remains in positive X direction for the entire range of disc angle. This is due to reversing of the propeller blade orientation. This means that the face and back of the aerofoil are interchanged for one half of the disc rotation. This is shown in Fig. 4. Although this is a necessity from the hydrodynamic loading view point, structurally it results in severe cyclic loading on the blade. In Fig. 4, it can be observed that the aerofoil works in combined heaving and pitching mode. It must be observed that the combined 3 DoF ship dynamics will cause significant variation in the inflow speed to the propeller blade at different blade orbit angle. Additionally, the continuous variation in the pitch angle (ı) means that the blade has an additional angular velocity (and angular acceleration) about the stock. Due to pitching of the blade, besides the hydrodynamic torque there will be additional inertial torque acting on the stock of individual blades. This will be investigated later in the paper. The frequency of pitching motion can be computed as shown in Eq. (16). ωf =

dı dt

(16)

From the working condition of the blade, it can be observed that ωf = / 0. As the blade is rotating about the disc center ‘C’, there will be centrifugal force acting on each blade. The direction of centrifugal force is outward normal to the disc at all points. In case of even number of blades, the net force in ‘X’ or ‘Y’ direction due to centrifugal effects will cancel out. In case of odd number of blades, a small residual value may remain. Therefore, the centrifugal force does not significantly affect the total propulsion thrust. However, while computing the stress on individual blade or the torque on the stock of individual blade, the centrifugal force needs to be considered. The centrifugal force is given by Eq. (17).









 Px ω  P x R cos i + R sin j d (C.F) = dm ω

(17)

The centrifugal force in Eq. (17) is for each discretized blade element. the centrifugal force on the entire blade is determined by integrating Eq. (17) over the blade domain. It may be noted that the centrifugal force will also cause a moment at the propeller blade stock position. It can be computed as shown in Eq. (18). Q d(C.F) = (x − S.P)









P × ω  P × R cos i + R sin j × ω

(18)

For obtaining the thrust and transverse force acting on the ship and the torque on the rotating disc, the lift and drag forces are

−  − ˛) −  − ˛)

 (19)

It can be observed that the forces acting on the blade are in mutually perpendicular planes. Besides, above, the load distribution is not uniform and the blade support is asymmetric. Therefore, the structural problem becomes that of combined bending and torsion. For a MCP propeller moving ahead, a major component of the centrifugal force is in the direction of normal force on the propeller for one-half of the disc. For the other half of the disc, the centrifugal force is opposite to the normal force. This may cause cyclic stresses on the blade. The reversal of the face and back surface of the aerofoil, described in Eqs. (11) and (15), exacerbates the variation in the cyclic stress. It must be noted that centroid of lift and drag act on the hydrodynamic center of the blade. Centrifugal force, which is inertial, acts on the mass centroid of the blade. Now we discuss about the torque acting on the stock of individual propeller blade. The stock of propeller blade can be located on the hydrodynamic center or the mass center depending on the magnitude of hydrodynamic and centrifugal force. The hydrodynamic center varies along the chord depending on the angle of attack and aerofoil characteristics of the blade. The mass centroid of the aerofoil can be assumed constant as the magnitude of structural deflections is much less than the dimensions of the blade. Therefore, if the propeller blade stock is at mass centroid, the torque component due to centrifugal force can be made zero. This will decrease the torque on the propeller blade stock thereby reducing the load on the prime mover. It must be noted that shear force and bending moment component may remain. Similarly, an axial force on the blade and stock due to self-weight of propeller blade and stock will also be present. All these forces and moment need to be considered while computing the effective stresses on the blade. These will be demonstrated later by computation. It is noted from Eq. (19) that there is a component of the force tangential to the chord of the blade. As the blade sectional moment of inertia is high in this direction, we presently ignore this force component. Similarly, we also ignore the axial stress on the blade stock due to self-weight of the propeller blade. However, while computing the effective stress of the stock these stresses must be taken into account.

3. Numerical methodology for flow analysis In the previous section propeller blade lift and drag due to rigid body motion characteristics were discussed. In this section hydrodynamic loading on the blade will be described. We assume inviscid fluid for computing the hydrodynamic forces and moment. Panel method is used for computing the flow dynamics around the blade. The panel method employs sources and vortices as proposed by Hess and Smith [8]. The potential equation can be written as Eq. (20).

= ∞ + s + v + B + 

(20)

The flow potential is determined by integrating the strength of source and vortices over the body surface as shown in Eq. (21). The

110

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

flow potential due to onset flow, pitching frequency of blade and the shed vortices are given in Eq. (22).



s =



v =

⎫

Eq. (25) can be simplified to get the local velocity components by simple geometric interpretation as shown in Eq. (26).

q (s) ln rds ⎪ ⎪ ⎬ 2  (s) ds 2

u∗sij = − (21)

⎪ ⎪ ⎭

∞ = VR (x cos ˛ + y sin ˛) 

B = ωf P (x sin + y cos )



 = d/2 r x sin  + y cos 

⎪ ⎪ ⎭

(22)

ui = VR cos ˛ + ωf Pi sin i +

N 

2

v∗vij =

rij+1 1 ln 2 rij

q(s) (s) [ lnr − ]ds + ... 2 2 panel

M 

d/2 r(xsin + ycos)

(27)

⎪ ⎭

−ui sin i + vi cos i = 0

(28)

u1 cos 1 + v1 sin 1 = −uN cos N − vN sin N

(29)

The velocity components are computed at the middle of each panel by the contribution from the onset flow, the sources and vortices on each panel and can be written as Eq. (30). qj usij + 

N 

uvij +

j=1 N 

M 

⎫ ⎪ ⎪ ⎪ ⎪ ⎬

dk /2 rk sin k

k=1

qj vsij + 

j=1

N  



⎫ ⎪ ⎬

ˇij

The flow tangency condition is applied at all the panels and Kutta condition on panels connecting the trailing edge tip. Therefore, for ‘N’ number of panels, we get ‘N + 1 equation. The flow tangency condition and the Kutta condition can be written as in Eqs. (28) and (29) respectively.

vi = VR sin ˛ + ωf Pi cos i + (Vrel )i +

... + ωf P (xsin + ycos ) +

2

u∗vij =

j=1

i=1

(26)

Similarly, the local velocity due to a unit-strength vortex on the ‘jth’ panel can be written as Eq. (27).

In our problem, the bound circulation on the foil changes at each time step as the inflow to propeller blade changes. This makes the problem unsteady and the formulation is done by the procedure explained in [3] and [9]. We assume that at each time step, a vortex is shed tangential to the trailing edge of the blade. The free vortex shed follows the trajectory of the trailing edge as shown in Fig. 5. The strength of this free vortex is equal to the change in the bound vortex but is opposite in sign. The surface of the body is discretized by a number of points, called nodes. The nodes are connected by straight lines, which become the panels as shown in Fig. 6. There are ‘N’ number of panels and ‘M’ number of shed vortices taken into consideration. We then distribute the sources and vortices on the straight-line panels and the potential is given by Eq. (23).

= VR (xcos ˛ + ysin ˛) +

⎫ ⎪ ⎬ ⎪ ⎭

ˇij

v∗sij =

⎫ ⎪ ⎪ ⎬

rij+1 1 ln 2 rij

N 

vvij +

j=1

M 

(30)

⎪ ⎪

dk /2 rk cos k ⎪ ⎪ ⎭

k=1

Eq. (30) is substituted in flow tangency condition and Kutta condition equations (Eqs. (28) and (30)). The resultant equation in the compact form is shown in Eq. (31). (23)

N 

k=1

Ai,j qj + Ai,N+1  = Bi + Ci + Di + Si

(31)

j=1

The source strength is taken to be constant on each panel, but variable from one panel to the other. The vortex strength remains same in all the panels. Therefore, if there are ‘N’ numbers of panels we have to find ‘N’ number of source strengths and single vortex strength. By using the flow tangency condition on each panel and the Kutta condition on trailing edge panels, we can get the source and vortex strength. The ‘jth’ panel is considered to be in between ‘jth’ and ‘j + 1th’ nodes, and its inclination with the x-axis be ‘ j ’ as shown in Fig. 6. The local velocity components of the panel can be converted to the global axis by Eq. (24). u = u∗ cos j − v∗ sin j



(24)

v = u∗ sin j + v∗ cos j

The local velocity components at ‘ith’ panel due to a unitstrength source distribution on the ‘jth’ panel can be written as Eq. (25).

⎫ ⎪  ⎪ ∗ ⎪ x −t 1 ⎪ ⎪ u∗sij = dt ⎪ 2 2 ⎬ (x∗ − t) + y∗2 ⎪ lj

0 lj

v∗sij =

1 2



0

y∗ (x∗ − t)2 + y∗2

dt

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭

(25)

The LHS of Eq. (31) is the contribution due to the sources and vortex, the RHS term ‘Bi ’ corresponds to the onset flow, ‘Ci ’ due to pitching frequency, ‘Di ’ due to blade vibration velocity ‘Vrel ’ and ‘Si ’ due to the influence of shed vortices. The source and vortex strength is computed by solving the set of algebraic equation obtained from Eq. (31) by putting it in a matrix form as in Eq. (32).

⎡

A1,2

→

A1,N

A1,N+1

A2,2

→

A2,N

A2,N+1

↓

↓

↓

↓

↓

↓

↓

↓

AN,2

→

AN,N

AN,N+1

AN+1,2

→

AN+1,N

AN+1,N+1

A1,1

⎢ ⎢ A2,1 ⎢ ⎢ ↓ ⎢ ⎢ ↓ ⎢ ⎢ A ⎣ N,1 AN+1,1

⎡

⎢ ⎢ ⎢ ⎢ =⎢ ⎢ ⎢ ⎢ ⎣

B1

⎤ ⎡

⎥ ⎢ B2 ⎥ ⎢ ⎥ ⎢ ⎢ ↓ ⎥ ⎥+⎢ ⎥ ⎢ ↓ ⎥ ⎢ ⎥ ⎢ BN ⎦ ⎣ BN+1

C1

⎤ ⎡

⎥ ⎢ C2 ⎥ ⎢ ⎥ ⎢ ⎢ ↓ ⎥ ⎥+⎢ ⎥ ⎢ ↓ ⎥ ⎢ ⎥ ⎢ CN ⎦ ⎣ CN+1

D1

⎤ ⎡

⎥ ⎢ D2 ⎥ ⎢ ⎥ ⎢ ⎢ ↓ ⎥ ⎥+⎢ ⎥ ⎢ ↓ ⎥ ⎢ ⎥ ⎢ DN ⎦ ⎣ DN+1

⎤⎧ ⎫ q ⎪ ⎪ 1⎪ ⎪ ⎪ ⎪ ⎥⎪ ⎪ ⎪ q ⎪ ⎥⎪ 2 ⎪ ⎪ ⎪ ⎥⎪ ⎪ ⎥⎨ ↓ ⎬ ⎥ ⎥⎪ ↓ ⎪ ⎥⎪ ⎪ ⎪ ⎪ ⎥⎪ ⎪ qN ⎪ ⎪ ⎪ ⎦⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎭ 

S1

⎤

⎥ ⎥ ⎥ ↓ ⎥ ⎥ ⎥ ↓ ⎥ ⎥ SN ⎦ S2

SN+1

(32)

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

111

Fig. 6. (Left) Discretization of propeller blade surface with nodes and panels. (Right) Notation of panel angles and its field distance.

The LHS terms ‘A1,1 ’ to ‘AN,N ’ correspond to the sources as shown in Eq. (33). Ai,j = u∗si,j sin(i − j ) + v∗si,j cos(i − j )

(33)

The terms from ‘A1,N+1 ’ to ‘AN,N+1 ’ is due the vortex distribution and expressed as Eq. (34).

Ai,N+1 =

N 

Eq. (38) is obtained by applying Kutta condition due to onset flow, pitching motion of blade, vibration velocity of blade and the shed vortices effect on the trailing edge panels.





BN+1 = −V R cos(1 − ˛) + cos(N − ˛)





CN+1 = −ωf P1 sin( 1 − 1 ) + PN sin( N − N )

DN+1 = −(V rel )1 cos(1 − /2) − (Vrel )N cos(N − /2)

v∗vi,j cos(i − j ) − uv∗i,j sin(i − j )

(34)

j=1

⎪ ⎪ ⎪ ⎪ ⎪ M M ⎪  d  d ⎪ k k ⎪ SN+1 = − sin(k − 1 ) − sin(k − N ) ⎪ ⎭ 2 rk 2 rk k=1

The terms from ‘AN+1,1 ’ to ‘AN+1,N ’ is the coefficients of source distribution obtained by applying Kutta condition and expressed as Eq. (35).

⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬

(38)

k=1

By solving Eq. (32), we get the strength of sources and vortex distributed over the panels. The tangential velocity on each panel is computed by the vector sum of the effective inflow and the effect of sources and bound and trailing vortices. The tangential velocity at the mid of each panel is computed by Eq. (39).

⎫ ⎪ ⎪ Vti = VR cos(i − ˛) + (Vrel )i cos(i − /2) + ωf Pi sin( i − i ) + sin(k − i ) + ... ⎪ ⎪ ⎪ 2 rk ⎪ ⎬ k=1 ⎡ ⎤ ⎡ ⎤ r sin( −  )ˇ ij+1 i j ij N N  ⎪ qj ⎢ ⎪ ⎥   ⎣ sin(i − j ) ln rij ⎦ ⎪ + ... + ⎪ ⎣ ⎦ ⎪ rij+1 2 2 ⎪ ⎭ − cos(i − j ) ln j=1 j=1 + cos( −  )ˇ M  d

k

rij

AN+1,j =



v∗sk,j sin(k − j ) − u∗sk,j cos(k − j )

(35)

k=1,N

The term ‘AN+1,N+1 ’ is vortex coefficient from Kutta condition and expressed as Eq. (36).

AN+1,N+1 =

N  



v∗vk,j sin(k − j ) + u∗vk,j cos(k − j )

(36)

k=1,N j=1

The RHS terms ‘B1 ’ to ‘BN ’ (corresponding to the onset flow), ‘C1 ’ to ‘CN ’ (corresponding to the pitching motion), ‘D1 ’ to ‘DN ’ (corresponding to the blade vibration velocity) and ‘S1 ’ to ‘SN ’ (corresponding to shed vortices) are expressed in Eq. (37). Bi = VR sin(i − ˛) Ci = ωf Pi cos( i − i ) Di = (Vrel )i sin(i − /2)

⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬

⎪ ⎪ ⎪ M ⎪  ⎪ ⎪ dk Si = cos(k − i ) ⎪ ⎪ ⎭ 2 rk k=1

, i = 1, 2.....N

(37)

i

j

(39)

ij

We assume that shed vortices up to a distance of 5 times the diameter of the propeller disc (5D) behind each blade have an influence on the fluid flow around the blade. Beyond this distance, the effect of shed vortex in fluid flow around blade is negligible. In Eq. (37), the value of ‘M’ is taken such that shed vortices up to a distance of 5 times the blade diameter are considered for computation. It will be found later, that MCP propeller is used for different flow regimes (cruising (high speed) and towing (low speed)). The distance of 5D becomes suitable for all the investigated flow regimes. It must be noted that as long as the ship has forward speed, the above procedure will work well. If the MCP propeller is used for thrust vectoring at a single point, the above procedure will no longer be valid. This is because the shed vortices will lie within the circular blade trajectory. Therefore, their influence is not diminished as their distance from the blade is not increasing. We need to apply more advanced theory for such case. The MCP blade axis follows the trajectory of a cycloid, but the leading and trailing edge trajectory is not purely cycloidal. The deviation in the path of leading and trailing edge from the axis path is shown in Fig. 7. The above flow hydrodynamics have been derived by considering the blade as a rigid body. However, due to external load, the blade undergoes structural deflection, and velocities and accelerations due to structural vibration. It needs to be investigated if the structural dynamics of the system has any influence on the hydrodynamic loading of the blade. For this purpose fluid structure interaction need to be considered for the propeller blade. Each one of the blade can be considered as a cantilevered plate with

112

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

and span gives the total lift force generated by the blade. This is shown in Eq. (43).

CFY =

Tip  L.E 

2 W AVR2

pi

 Root T.E

1+



dzi dxi

dxi dzi

(43)

The component of the force (discussed in Eq. (43)) parallel to the blade chord when integrated along the blade chord and span gives the induced drag force on the blade. We add to the induced drag force the viscous drag component to get the total drag force on the blade. This is shown in Eq. (44).

 ⎛ Tip L.E ⎞ Tip  L.E dz    pi dxi i CFX = CD dxi dz i ⎠(44)   dxi dzi + 2 ⎝ W AVR2 dzi 2

Root T.E

Fig. 7. Trajectories of leading edge, trailing edge and propeller blade stock during ship’s straight line motion.

fixed support at one end. It may be noted that the fixed support is not continuous for the full length of the blade root. Due to the blade deflection, the radial distance (CA as shown in Fig. 2) varies from the blade root to the tip. The angular velocity (ωp ) in Eq. (2) is usually kept constant. However, due to structural deformation, the inflow →

velocity (VAA in Eq. (7)) will vary. This way, the blade deflection is coupled with the hydrodynamic loading. Similarly, as the support for blade is asymmetric, the blade deflections at leading and trailing edge are also asymmetric. The angular rotation of the plate due to structural deformation will change the angle of attack derived in Eq. (11). As the blade chord is less as compared to span, the twist of the blade or angular rotation can be assumed as constant. We thereby get Eq. (40), i.e. the correction term to the angle of attack mentioned in Eq. (11). −1

˛ = sin



(yLE − yTE ) c



⎫

 w  2 ⎪ VR − Vti2 Back ⎬ 2

(Pi − P∞ )Face =

 w  2 VR − Vti2 Face 2

⎪ ⎭

(41)

The pressure difference between corresponding face and back panels of the blade is computed as shown in Eq. (42). p = 0.5w



(Vti )2face

− (Vti )2back

CL = CFy cos ˛ − CFx sin ˛



(42)

The pressure difference multiplied by panel area is the hydrodynamic force acting on a single panel. The component of this force normal to the blade chord when integrated along the blade chord

Root T.E



CD = CFx cos ˛ + CFy sin ˛

(45)

We are also interested in knowing the torque on the blade stock. For this purpose, we first determine the hydrodynamic force acting on each panel by multiplying the pressure difference (in Eq. (42)) and the panel area. We take the moment of this force about the blade stock. This moment when integrated along the blade chord and span gives the total hydrodynamic moment acting on the blade stock. This is shown in Eq. (46). Tip  L.E 

(40)

(Pi − P∞ )Back =

dxi

The CD in Eq. (44) is taken as flat plate friction resistance component, 0.075/(Log10 Rn − 2)2 . The range of Reynolds number for different speed is given in Table 2. The coefficient of lift and drag is computed by resolving the forces along the flow direction as shown in Eq. (45).



QHB =

It may be noted that the magnitude of ˛ will keep varying along the span of the propeller blade. This is another instance of the coupling of structural deformation with the hydrodynamic loading. We have already described earlier how the blade vibration (velocity normal to the cord) is coupled with the hydrodynamic load. As the plate is vibrating, it has surface velocities and accelerations. The blade velocity ‘Vrel ’, due to structural vibration, is added in Eqs. 29 ∼ 32 to incorporate the structural coupling. The method of determining blade lift and drag force will be described. The coefficient of pressure is computed at mid of each panel and assumed constant over the area of the panel. The expression for the coefficient of pressure is obtained using Bernoulli’s theorem as given by Eq. (41).

1+

Root T.E

pi 1+



dzi dxi

(xi − S.P) dxi dz

(46)

In Eq. (46) S . P refers to the distance of the stock point from the leading edge of the blade. As the blade is oscillating, the torque on the blade stock is expressed as shown in Eq. (47).





IBlade+stock + ıIBlade+stock ı¨ + QHB + QFriction + QC.F = QStock

(47)

The torque QStock mentioned in Eq. (47) refers to the single individual propeller blade stock torque and must be distinguished from the disc torque mentioned in Eq. (13). In Eq. (47), ıIBlade+stock is the added mass moment of inertia of the propeller blade and stock, while QFriction is the frictional losses in bearing and contact surface of the shaft. The centrifugal force on the blade acts at its mass center. The term QC.F in Eq. (47) can be made equal to zero by locating the stock at the mass center of the blade along the blade chord. This will help in reducing the blade torque. As the flow characteristics along the blade surface is known, the cavitations behaviour of the blade can be investigated. The cavitation number of the blade is computed by Eq. (48). ˝=

PA + w gd − PV 0.5w Vti 2

(48)

By knowing the variation of ˝, the location where the blade will start cavitating can be predicted. For MCP the advance coefficient and eccentricity ratios are defined differently from that of

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

113

Table 1 Operating conditions of the Tug. Voith Water Tractor 70 t BP

Speed (knots)

Resistance (kN)

Thrust per unit (kN)

Eccentricity Ratio (e = #    

Advance coefficient ( = u/ nD)

Load Coefficient ()

2

e1 R

e2 R

+

2

)

Cruising

Case 1 Case 2 Case 3

13.5 12 10

80 66 50

40 33 25

0.596 0.528 0.440

0.566 0.503 0.419

0.8276 0.7298 0.6083

Towing

Case 4 Case 5 Case 6

6 5 3

330 250 110

165 125 55

0.330 0.270 0.155

0.251 0.209 0.125

0.4140 0.3421 0.2023

40

0.025

Converged DeltaT

Converged Number of Panels 0.02

Deflection (m)

Thrust (kN)

20

0

-20

-40

-60

0.015

0.01

0.005

0

5

10

15 20 25 30 Number of Panels

35

40

-0.005

0

0.001

0.002 DeltaT (sec)

0.003

0.004

Fig. 8. (Left) Time series variation of the thrust of a single blade for different number of discretized panels. The discretized panels are on one side of the blade surface. The entire time series variation is plotted on a single point of X axis. (Right) Time series variation of the deflection at the trailing edge of blade tip for different time steps. The entire time series variation is plotted on a single point of X axis.

conventional screw propeller [17]. In this paper, we investigate two different operating conditions, i.e. cruising and towing. In order to have a common coefficient to compare the performance at the different operating conditions, we define a new load coefficient ‘’. It is computed as shown in Eq. (49).

$% =

2 +

 e 2 1

R

+

 e 2 & 2

R

2

2

∂ y ∂ h(x, z) + ∂t 2 ∂x2 2

∂ .... + ∂x∂z

% % D

%

D (2(1 − ))

2

2

∂ y ∂ y + ∂x2 ∂z 2

%

2

∂ y ∂x∂z

&&

&&

2

∂ + ∂z 2

% % D

˙ y¨ , t) = F(x, z, y, y,

2

2

∂ y ∂ y + ∂z 2 ∂x2

&&

⎫ ⎬

+ .... ⎪ ⎪

⎪ ⎪ ⎭ (50)

(49)

4. Vibration analysis The propeller blade is considered as a cantilever plate tapered along its length. It is fixed at the stock point to the disc and free at the other end hanging downwards. Each discretized panel of the plate is considered as a trapezoidal panel and modeled as a thin Kirchhoff plate. When the structure is immersed in fluid, the vibration of the structure is transferred to the fluid and give rise to fluid motion. As a result, there is a discernible increase in the kinetic energy due to the additional kinetic energy of the fluid. Because of increase in the kinetic energy, the natural frequencies of structures which are in contact with fluid, or immersed in fluid, decrease significantly as compared to the natural frequencies in vacuum. This problem is more severe in water as compared to air due to higher density of water. This problem is referred to as the fluid-structure interaction problem or the hydro elastic vibration of structures. The forced dynamics of a Kirchhoff plate in fluid can be modeled as shown in Eq. (50).

where, ‘’ is the density of plate, ‘h(x, z)’ thickness of the plate, ˙ t)’ represents the external force distribution on the ‘F(x, z, y, y, plate and includes the added mass effect (¨y dependence) on account of plate vibration in water, ‘D = Eh3 /12(1 − 2 )’, ‘E’ is Young’s modulus, ‘’ is Poisons ratio. In this problem we used the added mass values for different modes as proposed by Vu et al. [15]. It may be noted that ‘F’ became a function of ‘y’ and y˙ due to the explanation given in the previous section. In Eq. (50), only Rayleigh structural damping is considered. As we assumed the fluid to be inviscid, we do not presently consider the viscous fluid damping effect on propeller blade vibration. It may be noted that the viscous damping on account of fluid may be significant for this problem configuration. We now discuss about the solution method of Eq. (50) to know the plate vibration characteristics, i.e. time and spatial variation of ‘y’. We assume the problem to be linear and solve it by the mode superposition method. The mode shapes are computed by finite element method. Let 1 , 2 , 3 ,. . . n be the modes shapes. Then the solution will be of the form as shown in Eq. (51). y(x, z, t) =



(x, z)p(t)

(51)

Eq. (50) can be written as Eq. (52). y(x,z,t) = p1 1 + p2 2 + p3 3 + ..... + pn n

(52)

114

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

In order to get the time dependent, ‘p1 ’, ‘p2 ’, ‘p3 ’,. . . ‘pn ’, we need to substitute Eq. (52) in Eq. (50) and then multiply with the mode shapes and integrate throughout the domain to get a set of linear algebraic equations. Substituting Eq. (52) in Eq. (50) and integrating over the domain, we get Eq. (53). a b n    i=1 0

h(x, z)p¨ i i dxdz +

a b n     

D pi i,xx + pi i,zz

i=1 0

0





.... + D (2(1 − )) pi i,xz

 ,xz

0

dxdz =

a b n    i=1 0



the (a) minimum number of mode shapes required to capture the deflection characteristics of the blade (b) minimum number of discretized panels required to capture the deflection characteristics of the blade (c) maximum permissible time step to be used for

 

,xx

+ D pi i,zz + pi i,xx



⎫ ⎪ ⎪ ⎪ + .... ⎪ ⎪ ⎪ ,zz ⎬ (53)

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭

˙ t)dxdz F(x, z, y, y,

0

Multiplying Eq. (53) by the modes individually, we get Eqs. (54) and (55). a b n    i=1 0

h(x, z)p¨ i i 1 dxdz +

a b n     

D pi i,xx + pi i,zz

i=1 0

0





.... + D (2(1 − )) pi i,xz



0

1 dxdz =

,xz

a b n    i=1 0

a b n    i=1 0

hp¨ i i n dxdz +

0





.... + D (2(1 − )) pi i,xz

 ,xz

0

n dxdz =

a b n    i=1

0

,xx

+ D pi i,zz + pi i,xx



⎫ ⎪ ⎪ ⎪ + .... ⎪ ⎪ ⎪ ,zz ⎬ (54)

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭

˙ t) 1 dxdz F(x, z, y, y,

D pi i,xx + pi i,zz

0

 

0

a b n      i=1





 

,xx

+ D pi i,zz + pi i,xx



⎫ ⎪ ⎪ ⎪ + .... ⎪ ⎪ ⎪ ,zz ⎬ (55)

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭

˙ t) n dxdz F(x, z, y, y,

0

The above equations are second order ODE. The equations are solved numerically using Newmark – Beta method. The coefficients ‘␥’ and ‘␤’ used in Newmark – Beta integration scheme are 1.00 and 0.25 respectively. The equations can be put into matrix form as shown in Eq. (56).

⎡

⎤⎧ ⎫ ⎡ ⎤⎧ ⎫ ⎡ ⎤⎧ ⎫ ⎧ ⎫ b11 b12 .. b1n ⎪ p˙ 1 ⎪ c11 c12 .. c1n ⎪ p1 ⎪ ⎪ F1 ⎪ p¨ 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎨ ⎪ ⎨ ⎪ ⎥⎪ ⎥⎪ ⎥⎪ ⎬ ⎢b ⎬ ⎢c ⎬ ⎪ ⎨F ⎪ ⎬ .. a2n ⎥ p¨ 2 ⎢ 21 b22 .. b2n ⎥ p˙ 2 ⎢ 21 c22 .. c2n ⎥ p2 2 + + = (56) ⎥ ⎢ ⎥ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ ⎪ ⎣ : ⎦⎪ ⎣ : ⎦⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ : : ⎦⎪ : : : : : : : : : : ⎪ ⎭ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎩ ⎭ ⎩ ⎭ ⎩ ⎪ ⎭

a11

a12

..

a1n

an1

an2

..

ann

⎢a ⎢ 21 a22 ⎢ ⎣ : :

p¨ n

bn1

bn2

..

bnn

p˙ n

Eq. (56) looks like an ordinary second order differential equation. However, it must be noted that the force matrix Fi is coupled with pi and p˙ i . Therefore, Eq. (56) is nonlinear. Substituting the values of ‘p1 ’, ‘p2 ’, ‘p3 ’,. . ..‘pn ’ back in Eq. (52), we get the deflections at each node as a function of time. For this problem, we used first 20 mode shapes and find that after adding more number of mode shapes there is no significant change in the results. The propeller blade is discretized into 25 number of elements along the chord and 100 number of elements along the blade span. A convergence test was also carried out to demonstrate that the selected number of discretizations and mode shapes are sufficient for this problem. 5. Results and discussion The resultant inflow to the MCP blades is a vector sum of advanced speed, the disc rotational speed, the oscillation of the blade about its own vertical shaft axis and the induced velocity by the vortex shed by the blades. As the flow is not symmetric during the forward and backward movement of the blade, the angle of attack of the blade profile need to be optimized accordingly [11]. Before carrying out the analysis for different operating conditions, using the above model a convergence study was carried out initially. The aim of the convergence test was to determine

cn1

cn2

..

cnn

pn

Fn

numerical integration, which will capture the dynamics of the system accurately. The convergence for number of panels and time is shown in Fig. 8. It was observed that with 25 numbers of panels on one side along the chord of the blade, reasonably accurate results were obtained. The aspect ratio of the blade is 4.61 with tapering ratio of 0.8. The 25 divisions along the chord lead to 25 number of panels on one side of the aerofoil section. Similarly, the time convergence was checked and it was observed that a time step of 0.001 s or lower will capture the structural dynamics of the system. This time step size is for the minimum number of panels described earlier In this study, we consider the symmetric profile of NACA0012. The validation of the panel method code is done by computing the flow characteristics over an aerofoil of 1.0 m chord length placed in a grid with values of X ranging from −2 m to 3 m and values of Y ranging from −2 m to 2 m. The aerofoil is kept in the middle of the grid; in such way, that the leading edge tip is at origin (0, 0). The grid has an ‘X’ offset 2 m from leading edge tip and trailing edge tip. The grid interval is taken as 0.05 m and 0.01 m in X and Y directions respectively. Fig. 9 shows the flow around the aerofoil at 5◦ (Left) and 15◦ (Right) angle of attack with an inflow velocity of 20 m/s. It is found that the velocity magnitude increases in the

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

115

2 NACA 0012 Cheng [4] Panel

CL

1

0

-1

-2 -20

-10

0 a(deg)

10

20

Fig. 10. Comparison of coefficient of lift with experiment result.

Fig. 9. (Left) Velocity distribution around the aerofoil at 5◦ and (Right) 15◦ angle of attack. Fluid inflow speed is 20 m/s (steady) in both cases.

upper surface of the aerofoil as the angle of attack increases from 5◦ to 15◦ . The coefficient of lift is computed by panel method and compared with the experiments results [4] in Fig. 10. The analysis is done for the angle of attack ranging from −15◦ to 15◦ . It is found that the panel method results slightly deviate from the experimental results. This is because of ignoring viscous nature of the fluid and flow separation phenomenon. Since the panel method works on the principle of potential flow theory, the flow separation effect cannot be captured. MCP blades normally encounter flow angle of attack less than 15◦ . As the flow, separation occurs at high angle of attack, the inviscid flow theory holds well for MCP blades. The analysis has been done for cruising and towing condition of a tug in forward motion. In cruising condition, only the resistance of the tug is considered for computing the required thrust. In the towing condition, the resistance of the towed ship is added to the resistance of the tug for computing the required thrust. In this paper, the tug is considered as “Voith Water Tractor 70t BP” [18] and the towed ship is a KVLCC tanker [17]. Six operating conditions of the tug at various speeds were investigated. The details are given in Table 1. The propeller rotation speed is kept constant at 65 rpm for all the investigated cases. The analysis is carried out with MCP of model ‘VSP 36R6/300-2 of Voith [20]. The particulars of the MCP propeller are given in Table 3. The eccentricity value and blade orientation for a port side propeller rotating in anticlockwise direction for a straight run case is shown in Fig. 11. The movement of the blade is such that it moves aft as the disc moves from 0◦ to 180◦ positions and moves forward as the disc moves from 180◦ to 360◦ positions. The motion of mcp blade mimics the flapping motion of a dolphin tail [19]. The thrust is generated for all the disc positions, except at 90◦ and 270◦ angle positions.

Fig. 11. Orientation of propeller blade for a particular eccentricity at various position along the disc. The ship is considered to be moving from left to right side. Six number of blades are shown as black. The red blades only show the location at which delta angle is zero, no physical blades exist at that location. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The variation of the lift coefficient with steady and unsteady flow consideration is shown in Fig. 12. On the X axis besides time, the blade orbit angle is also indicated. For considering unsteady flow effects, the magnitude and direction of incremental shed vortices must be known. This information is generated only after a few rotations of the blade. Therefore, for t < 2.0 s, the lift coefficient values have not stabilized. It is observed that the influence of unsteady flow is high for 90◦ blade orbit angle. There are multiple reasons for this. One is, at this blade orbit angle, the trajectory is such that shed vortices are geometrically near the blade (see Fig. 7). The other is, at this blade orbit angle, there is a drastic change in the intensity of shed vortices due to interchange of the face and back surface of the blade as described in Fig. 4. Finally, at this blade

116

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

CRUISING 13.5 knots unsteady 12 knots unsteady 10 knots unsteady 13.5 knots steady 12 knots steady 10 knots steady

0

Blade orbit angle (deg) 360 720

TOWING 6 knots unsteady 5 knots unsteady 3 knots unsteady 6 knots steady 5 knots steady 3 knots steady 1080

0

Blade orbit angle (deg) 360 720

1080

4

Coefficient of lift (CL)

3

2

1

0

-1 0

1

2

3

0

1

Time (sec)

2

3

Time (sec)

Fig. 12. Variation of the lift coefficient for a propeller blade when it rotates along the disc.

TOWING 6 knots 5 knots 3 knots

CRUISING 13.5 knots 12 knots 10 knots 0

Blade orbit angle (deg) 360 720

1080

0

Blade orbit angle (deg) 360 720

1080

Flapping frequency (rad/sec)

12

8

4

0

-4

0

1

2 Time (sec)

3

0

1

2

3

Time (sec)

Fig. 13. Variation of pitching frequency of the propeller blade.

orbit angle the pitching frequency is higher. This will be described next. The lift coefficient has been nondimensionalized by the resultant of the inflow velocity and the tangential velocity at the blade

position due to rotation of the disc. Unsteady effect magnitude is higher for cruising condition as compared to towing condition. The variation of pitching frequency is shown in Fig. 13. It is noted that

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

117

CRUISING 13.5 knots unsteady 12 knots unsteady 10 knots unsteady 13.5 knots steady 12 knots steady 10 knots steady

0

360

0

1

Blade orbit angle (deg) 720 1080 1440

TOWING 6 knots unsteady 5 knots unsteady 3 knots unsteady 6 knots steady 5 knots steady 3 knots steady

1800

0

360

0

1

Blade orbit angle (deg) 720 1080 1440

1800

Angle of attack (deg)

20

10

0

-10

2 3 Time (sec)

4

5

2 3 Time (sec)

Fig. 14. Variation in the angle of attack of the propeller blade.

Fig. 15. Variation of inflow velocity to the propeller blade.

4

5

118

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

CRUISING 13.5 knots unsteady 12 knots unsteady 10 knots unsteady 13.5 knots steady 12 knots steady 10 knots steady

0

Blade orbit angle (deg) 360 720

TOWING 6 knots unsteady 5 knots unsteady 3 knots unsteady 6 knots steady 5 knots steady 3 knots steady 1080

0

Blade orbit angle (deg) 360 720

1080

80 60 40

Thrust (kN)

20 0 -20 -40 -60 -80 -100 -120

0

1

2

3

0

1

Time (sec)

2

3

Time (sec)

Fig. 16. Variation of the thrust of a single propeller blade.

Table 2 Range of Reynolds number for different operating conditions of the Tug. Voith Water Tractor 70 t BP

Speed (knots)

Reynolds number (Minimum ∼ Maximum) × 106

Cruising

Case 1 Case 2 Case 3

13.5 12 10

3.68 ∼ 13.31 4.21 ∼ 12.77 4.93 ∼ 12.06

Towing

Case 4 Case 5 Case 6

6 5 3

6.35 ∼ 10.63 6.71 ∼ 10.28 7.42 ∼ 9.56

0.6

Simulated results e = 0.59 e = 0.52 e = 0.44

Thrust coefficient (ks)

0.5

Experiment results e = 0.59 e = 0.52 e = 0.44

e = 0.64 0.4 e = 0.56 0.3 e = 0.48

the pitching frequency is higher for higher ship speed. There is no change in the period of the pitching motion, as the disc revolution is constant. This means that the blade has to change its direction of rotation rapidly in case of higher forward speed. From Table 1 it can be noted that higher ship speed corresponds to higher eccentricity ratio. The variation in the angle of attack of the propeller blade is shown in Fig. 14. On the X axis, both time and the blade orbit angle are shown. The angle of attack keeps varying with the blade orbit angle. The variation in angle of attack is very high (about 20◦ ) just after 90◦ blade orbit angle in cruising condition (13.5 knots). The sudden variation in angle of attack is observed in towing condi-

0.2 e = 0.40

Table 3 Particulars of the propeller.

0.1

-0 -0

0.1

0.2 0.3 0.4 0.5 Advance coefficient (l )

0.6

0.7

Fig. 17. Comparison of numerical simulations with experiment results of Voith [19].

Diameter RPM Top chord of blade Bottom chord of blade Length of blade Area of blade (one side) Number of blades Mass per blade Mass moment of inertia about stock of each blade

3.6 m 65 0.650 m 0.520 m 3m 1.755 m2 6 661 kg 19.625 kgm2

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

CRUISING 13.5 knots unsteady 12 knots unsteady 10 knots unsteady 13.5 knots steady 12 knots steady 10 knots steady

0

Blade orbit angle (deg) 360 720

119

TOWING 6 knots unsteady 5 knots unsteady 3 knots unsteady 6 knots steady 5 knots steady 3 knots steady

1080

0

Blade orbit angle (deg) 360 720

1080

300

Load pressure (kN/m2)

200 100 0 -100 -200 -300 -400

0

1

2 Time (sec)

3

0

1

2

3

Time (sec)

Fig. 18. Variation of pressure load at the leading edge tip of the propeller blade.

Fig. 19. Net pressure (hydrodynamic + centrifugal) on the propeller blade surface at every 45◦ interval on the disc. The values correspond to 13.5 knots cruising speed.

tion also. However, its magnitude is much lower. This fluctuation is due to the high pitching frequency at that position (see Fig. 13). At most of the blade orbit angle, the angle of attack is within the range (+10◦ to −10◦ ). The magnitude of angle of attack is higher in towing condition as compared to cruising condition. It is noted that the fluctuation in angle of attack is higher in unsteady flow as compared to steady flow condition. The variation in the inflow velocity to the propeller blade is shown in Fig. 15. On the X axis both time and the blade orbit angle are shown. The major variation in the inflow velocity is due to the direction of rotation being with the

inflow and against the inflow for each half of the disc revolution. Understandably, the variation in inflow velocity is lower in towing condition as compared to cruising condition. Sudden jump in inflow velocity (in unsteady flow condition) is observed in cruising condition after 90◦ blade orbit angle. It may be noted that at this location the variation in shed vortex intensity is maximum. The variation in the blade thrust is shown in Fig. 16. For both the cruising and towing speeds, the blade thrust varies with the blade orbit angle. The variation in the thrust values is higher at the 90◦ blade orbit angle. It takes a few revolutions for the unsteady hydrodynamic load to

120

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

Fig. 20. Deflection of the propeller blade at every 45◦ interval on the disc. The values correspond to 13.5 knots cruising speed.

CRUISING 13.5 knots unsteady 12 knots unsteady 10 knots unsteady 13.5 knots steady 12 knots steady 10 knots steady

0

Blade orbit angle (deg) 360 720

TOWING 6 knots unsteady 5 knots unsteady 3 knots unsteady 6 knots steady 5 knots steady 3 knots steady

1080

0

Blade orbit angle (deg) 360 720

1080

50 40

Deflection (mm)

30 20 10 0 -10 -20 0

1

2 Time (sec)

3

0

1

2

3

Time (sec)

Fig. 21. Deflection pattern of the trailing edge tip of propeller blade for different flow and load conditions.

become stabilized. This is because only after a few revolutions, the intensity and trajectory of shed vortex is known. It is observed that, for 6 knots towing speed, the magnitude of thrust is higher as compared to 13.5 knots cruising speed. The thrust magnitude changes after applying unsteady corrections. This phenomenon is observed, both in the cruising and towing condition. The negative thrust region is more in case of cruising as compared to towing. The decrease in thrust is due to the centrifugal force being out of phase with the hydrodynamic lift force. The time varying thrust for all the six propeller blades need to be determined, to compare

the average propulsive thrust for cruising and towing condition. Only after determining the average propulsive thrust and torque, the efficiency of the propulsion system can be estimated. It is found that blade is giving the least thrust at 90◦ and 270◦ position of the disc. The angle of attack at 90◦ and 270◦ position of the disc is zero, hence no lift is generated at this position. An important conclusion is higher pitching frequency, in cruising condition, does not correspond to higher thrust. The optimal combination of disc rotation rate and forward speed needs to be investigated. It must be noted that, higher pitching frequency will increase the inertial torque load

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

TOWING 6 knots 5 knots 3 knots

CRUISING 13.5 knots 12 knots 10 knots

0.005 0.0025

121

Rotation (deg)

0 -0.0025 -0.005 -0.0075 -0.01 -0.0125 -0.015 -0.0175

0

1

2

3

0

1

Time (sec)

2

3

Time (sec)

Fig. 22. Twisting of the blade tip section due to non uniform deflection between leading and trailing edge.

CRUISING 13.5 knots coupled 12 knots coupled 10 knots coupled 13.5 knots 12 knots 10 knots

0.8

TOWING 6 knots coupled 5 knots coupled 3 knots coupled 6 knots 5 knots 3 knots

Velocity (m/s)

0.4

0

-0.4

-0.8

-0.01

0

0.01 0.02 Deflection (m)

0.03

0.04 -0.02

0

0.02 Deflection (m)

0.04

0.06

Fig. 23. Phase plot of deflection and velocity at the trailing edge of propeller blade tip. The hydrodynamic loading is for unsteady flow condition with and without structural coupling.

thereby resulting in increased blade shaft torque. The thrust coefficient (ks ) is calculated and compared with experiment results obtained from Voith [19] and shown in Fig. 17. Experiment data is available from eccentricity ratio 0.40 onwards. Therefore, only cruising conditions could be compared. Experiment data is linearly interpolated to get the eccentricity values corresponding to numerical simulation. It is observed that simulation results overestimates the thrust coefficient by 0.06 at all the three advance coefficient. The trend of experiment data is well captured by numerical simulation. The numerical model needs to be improved by incorporating viscous effects of blades, rotating disc above the blades and mutual interaction of the blades,. These effects have not been considered in this work. The variation in the hydrodynamic pressure at the leading edge tip is shown in Fig. 18. It is observed that there is smooth variation in the load while applying steady inflow condition. However, for unsteady flow condition there is a sharp fluctuation in the

load. It is found that the direction of hydrodynamic pressure acting on the propeller blades changes from positive to negative as the blade changes its orbit angle. This is because the face and back side of the aerofoil gets reversed once during each revolution of the disc. The magnitude of hydrodynamic pressure on one half of the disc is higher than the other half. This is because of the influence of the centrifugal force. The centrifugal force always acts towards the outer side of the disc. the hydrodynamic force direction is always towards the direction of motion. When the blade moves from (: 270◦ → 0◦ → 90◦ ), centrifugal force and hydrodynamic force are “in phase”. When the blade moves from (: 90◦ → 180◦ → 270◦ ), centrifugal force and hydrodynamic force are “out of phase”. We need to do the vector sum of both these forces as they change their direction at each blade orbit angle. This leads to asymmetric loading on the blade for each rotation of the disc. This may cause blade vibration, which is investigated later. The fluctuation in load is higher

122

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

CRUISING 13.5 knots coupled 12 knots coupled 10 knots coupled

Spectrum (mm2/Hz)

15

TOWING 6 knots coupled 5 knots coupled 3 knots coupled

10

5

0 0.01

0.1 1 Frequency (Hz)

10 0.01

0.1 1 Frequency (Hz)

10

Fig. 24. Spectrum plot of the deflection of the blade tip at the trailing edge. The blade structural deformation is computed considering coupling of hydrodynamic load with structural deformation.

in cruising condition (13.5 knots) as compared to towing condition (6 knots). However, the overall hydrodynamic load is higher in 6 knots towing condition as compared to cruising condition. The variation of net hydrodynamic pressure on the entire blade surface at every 45◦ blade orbit angle is shown in Fig. 19. The variation is shown for 6 knots towing condition, as the propeller loading at this condition is maximum. The values shown are the difference in pressure between the two sides of the aerofoil section after the transients died out (t > 3 s). As mentioned earlier (see Fig. 4), the face and backside are interchanged for each half of the disc. The centrifugal force, also becomes in phase and out of phase with the hydrodynamic force. Therefore, the pressure values become negative. However, the resultant of the vector component of the force will mostly be in the X axis direction of the ship. The variation of blade deflection at different blade orbit angles is shown in Fig. 20. The deflection is shown for 6 knots towing with unsteady flow condition as the load at this condition is the maximum. The values shown are for t > 3 s, to exclude the influence of transients. The blade deflection shown are ‘y’ in the blade fixed coordinate system. It is observed that due to hydrodynamic and structural dynamics coupling, the deflection pattern is complex and does not mirror the trend of hydrodynamic/centrifugal load. The deflection pattern at the trailing edge tip of the propeller blade is shown in Fig. 21. The deflection is shown for 6 knots towing with unsteady flow condition as the load at this condition is the maximum. The values shown are for t > 3 s, to exclude the influence of transients. The blade deflection shown are the values ‘y’ in the blade fixed coordinate system. It is seen that the vibration is severe while applying the unsteady flow effects. It is also found that the blade deflection values in towing condition are comparable in magnitude to the cruising condition. There is a positive blade deflection value, in the blade fixed coordinate system. This is because in (: 270◦ → 0◦ → 90◦ ) sector, centrifugal force and hydrodynamic force are in the same direction. Therefore deflections are positive. In (: 90◦ → 180◦ → 270◦ ) sector, centrifugal force and hydrodynamic force oppose each other. However, centrifugal force being higher, and in the + Y axis, the blade deflection values remain positive but their magnitude is lower. The twisting of the blade tip section about its stock point is shown in Fig. 22. The twisting is due to non uniform deflection of leading and trailing edge of the blade. It is found that the twisting of blade is marginally higher in towing condition (5 and 6 knots

speed) as compared to cruising condition. The magnitude of twisting is low as compared to deflection. This is because the aerofoil is slender (span/chord ratio: 4.61). Therefore, although the blade stock support is asymmetric (i.e. does not extend for the full chord of the blade), it does not result in higher twisting of the blade. The phase plot of the blade deflection and vibration (velocity) is shown in Fig. 23. The deflection vibration is shown for 6 knots towing with unsteady flow condition as the load at this condition is the maximum. The values shown are for t > 3 s, to exclude the influence of transients. The blade deflection shown are in the ‘Y’ axis direction of the blade fixed coordinate system. Similarly, the vibration (or velocity) is also in the ‘Y’ axis direction of the blade fixed coordinate system. It is seen that the coupled effect is slightly reducing the magnitude of both deflection and vibration (or velocity). Frequency domain analysis for blade deflection and vibration (velocity) is done using Fast Fourier Transform. The deflection and velocity is considered at the trailing edge of the blade tip. Due to the layout of the blade, the trailing edge at blade tip is expected to have highest deflections and thereby velocities. Both the deflection and vibration (velocity) are measured in the blade fixed coordinate system. The spectral analysis results for deflection are shown in Fig. 24. It is observed that spectrum for deflection becomes insignificant after 10 Hz. Spectrum is strongest for 6 knots towing condition. A principal peak corresponding to 1.08 Hz is observed. This is the frequency of rotation of the disc. Higher harmonics, of lower spectral strength, are also observed till 10 Hz. Static deflection corresponding to DC frequency is not shown. The spectral analysis results for vibration (velocity) are shown in Fig. 25. The spectral strength is maximum for 6 knots towing condition. A principal peak corresponding to 1.08 Hz is observed. This is the frequency of rotation of the disc. Higher harmonics, of lower spectral strength, are observed till 15 Hz for cruising condition. For towing condition, higher harmonics, of lower spectral strength, are observed till 10 Hz. For velocity, there is no DC frequency for both cruising and towing condition. The Poincare plot for deflection against the load coefficient is shown in Fig. 26. It is found that magnitude of deflection for towing condition is higher than cruising for similar range of load coefficient. Vibration limits for main propulsion machinery is specified by Classification Societies. If the structural vibrations are within the safe working limit specified, then it is considered that no damage will occur to the structure. The limits for vibration specified

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

TOWING 6 knots coupled 5 knots coupled 3 knots coupled

CRUISING 13.5 knots coupled 12 knots coupled 10 knots coupled

100

123

Spectrum ((mm/s) 2 /Hz)

80

60

40

20

0 0.01

0.1

1 Frequency (Hz)

10

100 0.01

0.1

1 Frequency (Hz)

10

100

Fig. 25. Spectrum plot of the velocity of the blade tip at the trailing edge. The blade structural vibration velocity is computed considering coupling of hydrodynamic load with structural deformation.

Cruising 13.5 knots Cruising 12 knots Cruising 10 knots Towing 6 knots Towing 5 knots Towing 3 knots

40

Deflection (mm)

30

20

10

0

-10

0.1

0.2

0.3 0.4 0.5 0.6 0.7 Load Coefficient (m)

0.8

0.9

Fig. 26. Poincare plot of the deflection at trailing edge of propeller blade tip. The variation is shown with respect to the load coefficient.

(18.0 mm/sec, rms value) is exceeded for the operating conditions investigated in this paper [6]. This may be because we used conventional NACA0012 series aerofoil for blade sections. The vibration limits need to be checked for the actual design aerofoil section to ensure no damage occurs to the structure. The above analysis is done by assuming the stock position of the blade at one-fourth of the chord from the leading edge i.e. hydrodynamic center. It is noted from the above analysis that centrifugal force constitutes a major part of load and it acts at the mass centroid of the blade. From the aerofoil geometry it is found that, the mass centroid is at two-fifth of the chord from leading edge. It is found that the deflection magnitude and vibration amplitude become less when the stock position is shifted from hydrodynamic center to

mass centroid. Fig. 27 shows the change in deflection and the corresponding phase plot. The torque load on the stock will be lower if the stock is located at the mass centroid of the blade. The deflection and vibration is coupled with unsteady flow and its effect is shown in Fig. 28. Due to the coupling effect, the thrust is getting low and it is found the coupling effect is more in the towing (6 knots) as the deflection magnitude is high. As the prime mover rotates the disc, the torque of all the blades is added to get the total torque. The torque contributed by a single blade during its revolution is shown in Fig. 29. It is seen that the high-speed cruising condition is giving the higher torque as compared to towing. The torque of the prime mover depends on the speed and thrust required. The individual blade stock torque is computed by Eq. (47)

124

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

13.5 knots Hydodynamic centre 13.5 knots structural CG

30

0.3

0.2

Velocity (m/s)

Deflection (mm)

20

10

0.1

0

0 -0.1

-10

0

1

2

3

-0.2 -0.01

0

Time (sec)

0.01 Deflection (m)

0.02

0.03

Fig. 27. Left: Variation of the deflection of the blade tip at trailing edge with two different blade stock positions. In “Hydrodynamic center” the blade stock is located at effective hydrodynamic center and in “Structural CG” the blade stock is located at the mass centroid. Right: Comparison of the phase plot for the two cases.

CRUISING 13.5 knots coupled 12 knots coupled 10 knots coupled 13.5 knots 12 knots 10 knots

0

Blade orbit angle (deg) 360 720

TOWING 6 knots coupled 5 knots coupled 3 knots coupled 6 knots 5 knots 3 knots

1080

0

Blade orbit angle (deg) 360 720

1080

80 60

Blade thrust (kN)

40 20 0 -20 -40 -60 -80

0

1

2 Time (sec)

3

0

1

2

3

Time (sec)

Fig. 28. Variation of blade thrust with coupled and uncoupled model. Left: Cruising condition, Right: Towing condition.

and is shown in Fig. 30. As the blade changes its direction at 90◦ and 270◦ , we can see the torque is maximum and minimum at that position. It is seen that the fluctuation between minimum and maximum decreases as the eccentricity value decreases. The pressure in face and back of the blade is computed by knowing the tangential velocity at each location. There is an interchange in the face and back of the blade as it travels from 90◦ to 270◦ and 270◦ to 90◦ . The pressure profile of face and back at critical location like 0◦ , 45◦ , 180◦ and 225◦ is shown in Fig. 31. Accordingly, the cavitation number is computed by Eq. (48) and shown in Fig. 32. It is seen that the cavitation number is less in 225◦ i.e. while blade travelling from 180◦ to 360◦ (forward stroke). It is well known that less cavitation number leads to the probability of cavitation.

In this paper, the viscous effect of the fluid is ignored while computing the hydrodynamic loading on the propeller blade. It may be noted that the blades are working below a rotating disc. Due to viscosity, the rotating motions of the disc will change the inflow characteristics to the blade. There is mutual hydrodynamic interaction of the blade. This will also change the magnitude of the hydrodynamic force acting on each blade. In case of structural deformation, only the forces normal to the blade chord have been considered. It is observed that forces parallel to the chord will also act on the blade. Additionally, the weight of the blade will contribute to axial load on the blade. To incorporate these effects, shear deformation and rotary inertia of the plate need to be considered. For higher blade loading, the plate deflections/deformations will

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

TOWING 6 knots coupled 5 knots coupled 3 knots coupled

CRUISING 13.5 knots coupled 12 knots coupled 10 knots coupled

60

125

40 20

Torque (kNm)

0 -20 -40 -60 -80 -100 -120 -140 0

1

2

3

0

1

Time (sec)

-160

2

3

Time (sec)

-180 Fig. 29. Variation of the torque at the center of the disc due to a single propeller blade in different operating conditions.

CRUISING 13.5 knots coupled 12 knots coupled 10 knots coupled 13.5 knots 12 knots 10 knots

0

Blade orbit angle (deg) 360 720

TOWING 6 knots coupled 5 knots coupled 3 knots coupled 6 knots 5 knots 3 knots

1080

0

Blade orbit angle (deg) 360 720

1080

Blade stock torque (kNm)

4

2

0

-2

0

1

2 Time (sec)

3

0

1

2

3

Time (sec)

Fig. 30. Variation of blade stock torque with coupled and uncoupled model. Left: Cruising condition, Right: Towing condition.

be in nonlinear range. In such case, finite element method need to be used for computing structural dynamics. These effects will be incorporated in our future works. 6. Conclusions In this paper, the flow characteristics around the MCP propeller are computed using panel method. The effect of unsteady flow is incorporated. The analysis is carried out for six different straight running operating conditions. The following are the main conclusions of the work. • It is found that there is a significant change in the coefficient of lift while applying the unsteady formulation. The lift coefficient

magnitude continuously changes its value in the running condition. The unsteady effect is more in the higher eccentricity ratio operating condition. Most importantly, the unsteady flow effect is concentrated at a particular location on the disc. • The panel method used here ignores the viscous effect of the fluid. The viscous effects become more critical at higher angle of attack. In some operating conditions where eccentricity ratio is high, the angle of attack gets higher. The blade stresses are also higher during these conditions. In such cases, viscous correction has to be made to accurately predict the blade stresses. • The blade vibration is severe in the towing condition. Due to higher vibration and stresses, the fluid structure coupling loads must be taken into account while designing of blades.

126

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127

Fig. 31. Variation of difference in pressure at face and back with respect to the reference pressure. The values are shown at some critical blade orbit angles on the disc.

Fig. 32. Cavitation number at face and back of the blade at critical blade orbit angles.

• From the simulation results shown in this paper, it is found that the unsteady flow and fluid structure coupling have significant impact on the characteristics of the cycloidal propeller blade. Presently the investigations have been carried out only for straight running conditions. The investigations need to be carried out for maneuvering conditions also to check the blade stresses are within limits.

References [1] R.G. Allan, D. Molyneux, Escort tug design alternatives and a comparison of their hydrodynamic performance, In Proceedings, SNAME, Transactions of the Society of Naval Architects Vol. 112 (2004) 191–205. [2] D.S. Barrett, M.S. Triantafyllou, D.K.P. Yue, M.A. Grosenbaugh, M.J. Wolfgang, Drag reduction in fish-like locomotion, J. Fluid Mech. 392 (1999) 183-212. [3] J.P. Breslin, P. Anderson, Hydrodynamics of ship propellers, in: Ocean Technology Series 3, Cambridge University Press, United Kingdom, 1994.

[4] H.M. Cheng, Propeller design based on lifting surface theory part 1 uniform chord wise load distribution, in: Hydromechanics Laboratory Research and Development Report, Department of the Navy, 1964 (Sep). [5] J.P. Ghose, R.P. Gokarn, Basic Ship Propulsion, Allied Publishers Pvt. Limited, India, 2004 (ISBN 81–7764-606-0). [6] Guidance notes on ship vibration, American Bureau of Shipping, Table 4, pp. 61, April 2006. [7] W.L. Haberman, E.E. Harley, Performance of vertical axis (Cycloidal) propellers calculated by taniguchi’s method, in: Tech Report, Report 1564, SR0090101, Hydromechanics Laboratory Research and Development Report, Virginia, 1961. [8] J.L. Hess, A.M.O. Smith, Calculation of potential flow about arbitrary bodies, Prog. Aeronaut. Sci. 8 (1966) 1–138. [9] E.L. Houghton, P.W. Carpenter, Aerodynamics for Engineering Students, 5th ed., Butterworth Heinmann, 2003 (ISBN 0–7506-5111-3). [10] D. Jurgens, H.J. Heinke, Voith schneider propeller (VSP) – investigations of the cavitation behaviour, in: First International Symposium on Marine Propulsors, SMP’09, Trondheim, Norway, 2009 (June). [11] D. Jurgens, M. Palm, Voith schneider propeller – an efficient propulsion system for DP controlled vessels, in: Dynamic Positioning Conference, Heidenheim, Germany, 2009, October 13–14.

J.J. Prabhu et al. / Applied Ocean Research 64 (2017) 105–127 [12] M. Palm, D. Jurgens, D. Bendl, Numerical and experimental study on ventilation for azimuth thrusters and cycloidal propellers, in: Second International Symposium on Marine Propulsors, Smp’11, Hamburg, Germany, 2011 (June). [13] M. Sfakiotakis, D.M. Lane, J. Bruce, C. Davies, Review of fish swimming modes for aquatic locomotion, IEEE J. Oceanic Eng. 24 (2 (April)) (1999) 237–252. [14] S. Stefan, S. Jurgen, C. Kelly, M. Thomas, A cycloidal propeller using dynamic lift, in: 37th AIAA Fluid Dynamics Conference and Exhibit, AIAA 2007–4232, Miami, 2007. [15] V.H. Vu, M. Thomas, A.A. Lakis, L. Marcouiller, Effect of added mass on submerged vibrated plates, Conference Proceedings of 25th Seminar on Machinery Vibration, Canadian Machinery Vibration Association (2007).

127

[16] P. Waclawek, D. Molyneux, Predicting the performance of a tug and tanker during escort operations using computer simulations and model tests, In Proceedings, SNAME Vol. 108 (2000) 21–43. [17] http://simman2008.dk/KVLCC/KVLCC1/kvlcc1 geometry.html (accessed 25.12.2016). [18] http://resource.voith.com/vt/publications/downloads/1371 e g1810 en voith-water-tractor-vwt.pdf (accessed 25.12.2016). [19] http://voith.com/de/1523 e g 2070 e vsp 2013-04-23 screen.pdf (accessed 25.12.2016). [20] http://voith.com/en/vsp-types-and-dimensions.pdf (accessed 25.12.2016).