Influence evaluation of wake shedding on natural frequencies of marine propulsion shafts

1. Introduction

Propulsion shafting is an important unit of marine components. It consists of a shaft, propeller and some support bearings. Usually, vibrations are inevitable for propulsion shafts because the propeller has inescapable eccentricity due to machining accuracy or wear and tear. Meanwhile, the propeller also suffers from large fluid pulse pressure due to the non-uniform flow. Vibrations of propulsion shafts have great effect on crew comfort, ship equipment strength and consequently on the ship safety. Therefore, the dynamic analyses for propulsion shafts are very important for designers. Many studies [1, 2] have investigated about the vibration characteristics and vibration control of propulsion shafts, including the prediction of natural frequencies and mode shapes, the calculation of forced vibrations and so on. At the design stage of propulsion shafts, natural frequencies and mode shapes must be analyzed to know the possible resonance speed. And it is also very important to evaluate them to avoid the resonance phenomena that can be caused by periodic forces of propulsion engine and propeller. Generally, in the dynamic model of the propulsion shaft system, the shaft could be assumed to be a beam supported with some linear springs in bearing points and the propeller could be considered as a lumped mass attached to the end of the beam. So, the mechanical attributes of the entire system could be described with relative clarity. But there are still some problems to predict natural frequencies. One of the reasons is that the propeller immerses in water and fluid-structure interaction occurs. The fluid force has great influence on propulsion shaft vibrations. Some studies [3, 4] show that part of these forces are proportional to vibration acceleration of wetted surface. Therefore, these forces can be replaced by added mass. Added mass is usually significant because on the one hand, it reaches nearly up to the propeller self-mass in some cases, and on the other hand, the added mass will play a dominant role in the entire system, due to the propeller location at a free end of the shaft like a cantilever beam with large mass at its end. So, the added mass cannot be omitted in the dynamic analysis of the propulsion shaft system. The natural frequencies of the shaft are significantly lower than those in air due to added mass. So even the self-mass properties of the propeller can be found without much difficulty, it is the added mass that is the great challenge to calculate the natural frequencies. Some studies dealt with added mass about propellers are simple and imprecise. For example, Murawski and Charchalis [5] studied torsional vibrations of propulsion shafts and Murawski [6] also investigated axial vibrations. Zou et al. [1, 7, 8] analyzed nonlinear forced vibrations of a marine propulsion shaft. In these studies, the added mass was estimated using empirical equations. Exact calculation of propellers added mass is important to accurately predict natural frequencies of propulsion shafts. There are some studies about the calculation of added mass caused by fluid. Lamb [9] described the calculation method about added mass in detail. But it is main to non-lifting bodies and the geometry shape is simple (e.g. cylinders and spheres). For lifting bodies (e.g. propeller), there is vortex shedding (wake vortex) from a sharp edge (e.g. the trailing edge of a propeller when it rotates). This wake is usually called the wake shedding or trailing vortex wake in ship engineering. It is controversial about the effect of the wake vortex on added mass. So its effect on natural frequencies is also unclear.

On the one hand, some studies think the effect of wake on added mass could be ignored. For example, Lee et al. [10] studied hydroelastic responses of marine propellers. Hutchison et al. [11] studied ducted propellers added mass. Ghassemi and Yari [12] evaluated the added mass of marine propellers. Lin and Tsai [13] analyzed free vibrations of a composite propeller blade. In these analysis, the effect of wake shedding was ignored. On the other hand, some studies think the wake has important effect on added mass. For example, Michelin [14] thought that the separation and vortex shedding had significant influence on added mass when the body had a sharp edge. These vortices could produce the thrust and drag. He et al. [15] calculated the propeller added mass. Esch et al. [16] evaluated the added mass of propellers due to torsional and axial vibrations. Their analysis showed the wake shedding had an influence on the magnitude. Young [17, 18] calculated the added mass about the flexible composite propellers. In these analyses, the effect of wake shedding was considered. In addition, Gaschler and Abdel-Maksoud [19] studied the hydrodynamic mass coefficients for a cavitation and non-cavitation marine propeller. The non-rotating (in this case, the effect of wake was ignored) and rotating (in this case, the effect of wake was included) conditions were considered. MacPherson et al. [20] also estimated the axial and torsional added mass with rotating and non-rotating conditions.

To sum up, it is unclear about the effect of the wake shedding on added mass for marine propellers. Some studies think the wake could be ignored and some studies don’t think so. But no papers make a comparison between the case with consideration of wake shedding and the case without wake shedding. And no papers evaluate the effect of wake shedding on natural frequencies of propulsion shafts. Therefore, the comparison and evaluation of the influence of wake shedding on propellers added mass and propulsion shafts natural frequencies are the subject of this paper.

2. Mathematical formulations

In our study, the natural frequencies of propulsion shafts are calculated by FEM, and the added mass of propellers is determined by BEM. The detail description of the formulation and implementation is provided below.

2.1. Calculation of natural frequencies

A typical marine propulsion shafts is composed of a shaft, propeller, back stern bearing, front stern bearing, middle bearing and thrust bearing, which is illustrated in Fig. 1. In general, the propeller vibrates with the end of a shaft. It is assumed that the displacement of the end shaft is $\left\{q\right\}=[{\delta }_x,{\delta }_y,{\delta }_z,{\theta }_x,{\theta }_y,{\theta }_z{]}^T$. Many investigations in linear rotor dynamics deal with the problems of natural, unbalance and transient vibrations. Powerful approximation methods, e.g., the finite element method (FEM), are available for solving these problems. So, in our study, the FEM is applied to solve the natural frequencies.

Fig. 1A diagram of a typical marine propulsion shafts

In the FEM model, the shaft is simulated with Timoshenko beam element which has six degree of freedom. The bearings are modelled by linear spring element. The propeller is modelled by lumped mass element. The FEM model is showed in Fig. 2.

Fig. 2The FEM model of propeller-shaft-bearings system

The system equation of motion can be expressed as [21]:

1

$\left(\mathbf{M}+{\mathbf{M}}_a\right)\ddot{\mathbf{q}}+\mathbf{C}\dot{\mathbf{q}}+\mathbf{K}\mathbf{q}=\mathbf{f}\left(t\right),$

where $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$ are the global mass, damping and stiffness matrices respectively. The detail formulation of $\mathbf{M}$, $\mathbf{C}$ and could see the reference [21]. ${\mathbf{M}}_a$ is the added mass matrix of propellers. $\mathbf{f}\left(t\right)$ is external loads.

The natural frequencies can be obtained by solving the following equations:

2

$\left|-\left(\mathbf{M}+{\mathbf{M}}_a\right){\omega }^2+\mathbf{K}\right|=0,$

where $\omega$ is the natural frequency.

So, the added mass of propellers is an important parameter for the calculation of natural frequencies. It would decrease natural frequencies compared with those in air. The next section will describe the method for determining the added mass.

2.2. Calculation of added mass

Young [17, 18] have calculated the added mass about the flexible composite propellers with the boundary element method (BEM). But in their analysis, the propeller is elastic body and has infinite degrees of freedom. So, the added mass is distributional and should be expressed combined with FEM. In this study, the boundary element method described by Young is modified and fit for solving the added mass of rigid body which has six degrees of freedom.

2.2.1. Basic theory

The following assumptions are employed. (1) the flow is inviscid, incompressible, and irrotational; (2) the propeller immerses in deep water and the effect of free surface is ignored; (3) cavitation is not considered; (4) the geometry model of wake shedding is linear and prescribed helical surface.

The pressure on the propeller is caused by both parts, first it originates from the steady rotation of the propeller blades through the non-uniform wake field, and second it arises from the propeller vibration [18]. As long as the vibration amplitude is small, independent evaluation of these two pressure is possible. Because the added mass is only relevant to the second part, only the pressure caused by the propeller vibration is studied in this paper. This pressure can be written as (for more detailed derivation, see reference [17, 18]):

3

$P_{\phi }=\rho \left[-\frac{\partial \phi }{\partial t}-{\mathbf{V}}_{in}\cdot \nabla \phi \right],$

where $\phi$ is the perturbation velocity potential due to the propeller vibration. ${\mathbf{V}}_{in}$ is the inflow velocity of the propeller and $\rho$ is the density of water.

For the perturbation velocity potential $\phi$, combining the kinematic boundary condition and applying Green’s identity, the following expression can be obtained [22]:

4

$2\pi \phi \left(P\right)={\iint }_{S_B}\phi \left(Q\right)\frac{\partial }{\partial n_Q}\left(\frac{1}{R(P,Q)}\right)dS+{\iint }_{S_W}\mathrm{\Delta }\phi \left(Q_1\right)\frac{\partial }{\partial n_{Q_1}}\left(\frac{1}{R(P,Q_1)}\right)dS$
$+{\iint }_{S_B}\left({\mathbf{V}}_b\cdot {\mathbf{n}}_Q\right)\left(\frac{1}{R(P,Q)}\right)dS,$

where $P$, $Q$ are points on the propeller surface $S_B$. $Q_1$ is any point on the wake surface $S_W$ and $\mathrm{\Delta }\phi \left(Q_1\right)$ is potential jump across the wake surface $S_W$. ${\mathbf{n}}_Q$ is the unit normal vector of the propeller surface $S_B$. ${\mathbf{V}}_b$ is the vibration velocity of the propeller caused by the shaft vibration.

It assumes that the vortex sheets will leave the trailing edge at constant pitch angle, so the wake geometry which could be express as [23]:

5

$\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_w\left(r\right)=\left(1-C\right)\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_p+C\mathrm{t}\mathrm{a}\mathrm{n}\beta ,$

where $C$ is a constant between 0 and 1. Usually, $C$ is set equal to 0.2. ${\beta }_w\left(r\right)$ is the pitch angle of trailing vortices at radius $r$. ${\beta }_p$ is blade pitch angle which depends on the blade geometry. $\beta$ is undisturbed-flow pitch angle, which could be written as $\beta =\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(V_0/\left(\omega r\right)\right)$ determined by the advance coefficient $J$ ($J=\left(\pi V_0\right)/\left(\omega r\right)$).

The boundary integral Eq. (4) is hard to solve analytically for complex geometry. So, a numerical solution is obtained by dividing the surface into a number of small elements. It is assumed that the numbers of the panels on the propeller and wake surface are respectively $N$ and $N_W$ shown as in Fig. 3 (only one blade wake surface is shown).

Fig. 3Panel arrangement of propeller and wake

The boundary integral Eq. (4) can be expressed in matrix form as:

6

$\left\{\phi \right\}=\left[H\right]\left\{{\mathbf{V}}_b\cdot \mathbf{n}\right\},$

where $\left[H\right]$ is influence coefficient matries and solved by Morino’s analytical formulation [24].

Take the derivative to coordinate $x$ (or $y$ and $z$) on both sides of Eq. (4). Similarly, the following matrix form can be obtained:

7

$\left\{\frac{\partial \phi }{\partial x_i}\right\}=\left[H^{x_i}\right]\left\{{\mathbf{V}}_b\cdot \mathbf{n}\right\},x_i=x,y,z,$

where$\left[H^{x_i}\right]$ is the velocity influence coefficient matrixes and also calculated by Morino’s analytical formulation [24].

2.2.2. Added mass formulation

As shown in Fig. 2, the vibration velocity of the shaft is $\left\{\dot{q}\right\}=[{\dot{\delta }}_x,{\dot{\delta }}_y,{\dot{\delta }}_z,{\dot{\theta }}_x,{\dot{\theta }}_y,{\dot{\theta }}_z{]}^T$. For propeller, the coordinate of the $i$th panel center is assumed to be $(x_i,y_i,z_i)$. So, its vibration velocity could be expressed as:

8

${\overrightarrow{V}}_{bi}=\left({\dot{\delta }}_x+z_i{\dot{\theta }}_y-y_i{\dot{\theta }}_z,{\dot{\delta }}_y+x_i{\dot{\theta }}_z-z_i{\dot{\theta }}_x,{\dot{\delta }}_z+y_i{\dot{\theta }}_x-x_i{\dot{\theta }}_y\right).$

The normal vector ${\overrightarrow{n}}_i$ of the $i$th panel center is also assumed to be $(n_{xi},n_{yi},n_{zi})$. So ${\overrightarrow{V}}_{bi}\cdot {\overrightarrow{n}}_i$ could be written as:

9

${\overrightarrow{V}}_{bi}\cdot {\overrightarrow{n}}_i=\left[T_i\right]\left\{\dot{q}\right\},$

where $\left[T_i\right]=\left(n_{xi},n_{yi},n_{zi},n_{zi}{y}_i-n_{yi}{z}_i,n_{xi}{z}_i-n_{zi}{x}_i,n_{yi}{x}_i-n_{xi}{y}_i\right)$.

According to Eq. (3), then the pressure of the $i$th panel is obtained:

10

$P_{\phi }^i=-\rho \sum _{j=1}^N{H}_{ij}\left[T_j\right]\left\{\ddot{q}\right\}-\rho \sum _{j=1}^N\left(V_{in}^x{H}_{ij}^x+V_{in}^y{H}_{ij}^y+V_{in}^z{H}_{ij}^z\right)\left[T_j\right]\left\{\dot{q}\right\},$

where $V_{in}^x$, $V_{in}^y$, $V_{in}^z$ are three coordinate components of vector $V_{in}$.

Then the resultant force and moment by all panels could be written as:

11

$FM=-\rho \sum _{i=1}^N\sum _{j=1}^N{S}_i{H}_{ij}{\left[T_i\right]}^T\left[T_j\right]\left\{\ddot{q}\right\}-\rho \sum _{i=1}^N\sum _{j=1}^N{S}_{i}H{V}_{ij}{\left[T_i\right]}^T\left[T_j\right]\left\{\dot{q}\right\},$

where $S_i$ is the area of the $i$th panel and $H{V}_{ij}=\left(V_{in}^x{H}_{ij}^x+V_{in}^y{H}_{ij}^y+V_{in}^z{H}_{ij}^z\right)$.

In Eq. (11), the first term is proportional to vibration acceleration of the propeller, so it could be called added mass. Similarly, the second term is proportional to vibration velocity and could be called added damping.

So, from Eq. (11), it is known that the added mass could be written as:

12

$\left[M_a\right]=-\rho \sum _{i=1}^N\sum _{j=1}^N{S}_i{H}_{ij}{\left[T_i\right]}^T\left[T_j\right].$

In general terms, a propeller vibrates in the six rigid body modes defined in Fig. 2. So, the added mass is six order matrix and has the form [3]:

13

$\left[M_a\right]=\left\{\begin{array}{l}{m}_{ij},\text{if}i=\mathrm{1,4}\text{and}j=\mathrm{1,4}\text{or}i=\mathrm{2,3},\mathrm{5,6}\text{and}j=\mathrm{2,3},\mathrm{5,6},\\ 0,\text{if}i,j=\text{ot}\text{h}\text{ers},\end{array}\right.$

where $m_{11}$ is called axial added mass and $m_{22}$ is called lateral added mass. $m_{44}$ is called polar added mass moment of inertia and $m_{55}$ is called diametric added mass moment of inertia. For the added mass matrix, first, several of the terms, for example $m_{22}$ and $m_{33}$, have identical values. Second, the matrix is symmetrical. Third, there also are some non-diagonal coupling terms.

3. Program codes validation studies

Some numerical convergence and validation studies of the current BEM program codes (with MATLAB) have been conducted in this section. The added mass of a sphere immersed in infinite water is calculated, becasue for it exact analytical results are available. Then the exact results and the numerical results are compared. The radius of the sphere $r$ is 0.1 m. According to the analytical solutions in mass matrix, only $m_{11}$, $m_{22}$ and $m_{33}$ have values and the others are zero due to the geometry shape. The three coefficients are equal with each other. For a non-lifting body, there is no wake shedding from the body. Its panel arrangement is shown in Fig. 4. The mesh numbers on the sphere surface are 1200.

Fig. 4Panel arrangement of the sphere

The calculation values for added mass of the sphere using the BEM codes are as follows:

In Eq. (14), only $m_{11}$, $m_{22}$ and $m_{33}$ have values and the others are zero, which agrees with the analytical solution. For a sphere, the analytical values of added mass are equal to [25]:

Compared the numerical solutions with the analytical solutions, it is known that the errors between them are very small. In order to evaluate the mesh convergence, the added mass matrix is investigated for three different mesh elements (1200, 2400 and 4800 elements), using the BEM codes and the obtained results and errors are checked. The results are shown in Table 1.

It is evident that the convergence characteristics of the mass terms are very good. With the increase of mesh numbers, the numerical values converge to the analytical values.

To sum up, from this example, it can be found that there is a very good accuracy between the numerical solutions and the experiment data or the analytical solutions, suggesting the validity of our BEM program codes written in MATLAB.

4. Numerical example

4.1. The effect of wake on added mass

In this section, some numerical examples are considered to examine the influence of the wake vortex emanating from trailing edge on added mass. A series of four five-blade skewed marine propellers originate in practical engineering are evaluated. The parent propeller has a symmetric blade outline with zero skew, and the other three propellers have maximum projected skew angles of 36°, 72° and 108° at the propeller tip. The geometric parameters for all four propellers are the same except skew, pitch, and camber. The propellers each have a diameter of 5 m and the design advance coefficient $J=$0.889. The primary geometries are given in [3, 26], and are shown in Fig. 5.

Fig. 5BEM representation of propeller 4381, 4382, 4383 and 4384 and wake shedding

Two models are analyzed, in which one ignores the wake shedding and the other considers it. In the second model, the trajectories of wake shedding are defined by Eq. (5). The shape of wake shedding is shown in Fig. 5 (take one blade wake for example).

The mesh numbers are subdivided as many as possible until the numerical results have little variation and are convergent. These results have been non-dimensionalized by the divisor given in Table 2 (where $\rho$ is the density of water and $D$ is the diameter of the propeller).

Table 2Nondimensionalization of added mass

Divisor	Corresponding $i$ and $j$ number in $m_{ij}$	Type of coefficient
$\rho D^3/100$	$i,j=$1, 2, 3	Mass terms
$\rho D^4/100$	$i=$1, 2, 3; $j=$4, 5, 6 or $i=$4, 5, 6; $j=$1, 2, 3	Coupling mass terms
$\rho D^5/1000$	$i,j=$4, 5, 6	Mass moment of inertia terms

Fig. 6 is the results comparison of the main diagonal added mass between these two cases, in which one ignores wake shedding and the other considers it. The superscript $w$ denotes the added mass with consideration of wake shedding. From this figure, it is seen that, first, the wake greatly increases lateral added mass $m_{22}$ and diametric added mass moment of inertial $m_{55}$ for all skew angles. Second, the wake has minor influence on polar added mass moment of inertia $m_{44}$ when the skew angle is not large (such as propeller 4381 and 4382). As the skew angle increases, so does the effect of wake on $m_{44}$. Third, for axial added mass $m_{11}$, the wake has little effect on it when the skew angle is 72 degrees (such as propeller 4383). The wake increases $m_{11}$ when the skew angle is larger than 72 degrees (such as propeller 4384), and reduces it if the skew angle is smaller than 72 degrees (such as propeller 4381 and 4382).

Fig. 6Comparison of the main diagonal values with and without wake shedding

4.2. The effect of wake on shafts’ natural frequencies

In the previous section, the effect of wake shedding to added mass is discussed. As the introduction says, the ultimate aim of exact calculation the added mass is to predict the natural frequencies accurately. So, in this section, the effect of wake shedding on natural frequencies of propulsion shafts is also studied. Four dynamic models are established according a marine propulsion shafting in practical engineering and the structure is shown in Fig. 1. In these four models, the only difference is the propeller (the propeller is 4381 to 4384 respectively). The length of the shaft is 20 m. The external diameter of the cross section is 500 mm and the internal diameter is 320 mm (to simplify, the shaft is considered as uniform circular in cross section). The Young’s modulus is 210 GPa and the Poisson’s ratio is 0.33. The density is 7800 Kgm^-3. The axial coordinate of transverse bearing $k_1$, $k_2$ and $k_3$ is 2 m, 10.4 m and 20 m respectively. The stiffness are all 3e8 Nm^-1. The axial coordinate of thrust bearing $k_t$ is 20 m and the stiffness is 5e8 Nm^-1. Most parameters of the shaft come from the reference [27] and a small number of them, for example the thrust bearing stiffness, come from our experience. The natural frequencies of two models are calculated. In one model, the added mass is without wake shedding and in the other model, the added mass is with wake shedding. Then these results are compared with each other. $f_a$ and $f_a^w$ are defined as the axial (longitudinal) natural frequencies with or without wake respectively. Similarly, $f_l$ and $f_l^w$ are defined as the lateral (transverse) natural frequencies. $f_t$ and $f_t^w$ are defined as the torsion natural frequencies. The first three natural frequencies are studied and these results are shown in Fig. 7 to Fig. 9. Meanwhile, the relative errors are defined and shown in Table 3 to Table 5.

Fig. 7Comparison of the first natural frequencies with and without wake shedding

a) The axial natural frequencies

b) The lateral natural frequencies

c) The torsion natural frequencies

From these tables and figures, some conclusions can also be obtained. First, for lateral, axial and torision vibration, because the first natural frequencies occur near the propeller and the second, the third occur away from the propeller, the wake shedding only has effect on the first natural frequencies. Second, the first lateral natural frequency ${\mathrm{f}}_{\mathrm{l}}^{\mathrm{w}}$ is smaller than ${\mathrm{f}}_{\mathrm{l}}$ for any skew angle. The smaller the skew angle is, the larger the relative error will be. Even if the skew angle is 108 degrees, the relative error reaches up to 4.63 %. So, the wake shedding has important effect on lateral natural frequencies and should be considered. Second, the wake has great impact on the first torsion natural frequency when the skew angle is large. For example, when the skew angle is 108 degrees, the relative error amounts to 5.10 %. As the skew angle is small, the wake could be ignored. Third, for axial vibration, the influence of wake is little and may be ignored for all skew angles.

Fig. 8Comparison of the second natural frequencies with and without wake shedding

a) The axial natural frequencies

b) The lateral natural frequencies

c) The torsion natural frequencies

Fig. 9Comparison of the third natural frequencies with and without wake shedding

a) The axial natural frequencies

b) The lateral natural frequencies

c) The torsion natural frequencies

Table 3The relative errors of the first natural frequencies

	4381	4382	4383	4384
$\left|\left(f_a^w-f_a\right)/f_a^w\right|\times 100\%$	1.18 %	0.95 %	0.09 %	1.07 %
$\left|\left(f_l^w-f_l\right)/f_l^w\right|\times 100\%$	6.52 %	6.51 %	5.90 %	4.63 %
$\left|\left(f_t^w-f_t\right)/f_t^w\right|\times 100\%$	0.28 %	0.57 %	2.92 %	5.10 %

Table 4The relative errors of the second natural frequencies

	4381	4382	4383	4384
$\left|\left(f_a^w-f_a\right)/f_a^w\right|\times 100\%$	0.55 %	0.45 %	0.05 %	0.55 %
$\left|\left(f_l^w-f_l\right)/f_l^w\right|\times 100\%$	0.15 %	0.17 %	0.13 %	0.11 %
$\left|\left(f_t^w-f_t\right)/f_t^w\right|\times 100\%$	0.008 %	0.017 %	0.089 %	0.157 %

Table 5The relative errors of the third natural frequencies

	4381	4382	4383	4384
$\left|\left(f_a^w-f_a\right)/f_a^w\right|\times 100\%$	0.14 %	0.12 %	0.01 %	0.15 %
$\left|\left(f_l^w-f_l\right)/f_l^w\right|\times 100\%$	1.14 %	1.10 %	0.95 %	0.74 %
$\left|\left(f_t^w-f_t\right)/f_t^w\right|\times 100\%$	0.002 %	0.003 %	0.016 %	0.028 %

5. Conclusions

In this study, the influence of wake shedding from the propeller trail edge on natural frequencies is investigated. The most important results can be expressed as:

1) For lateral, axial and torsion vibration, the wake shedding has great effect on the first natural frequencies compared with the others. This because the first natural frequencies occur near the propeller.

2) For the first lateral natural frequency, the max relative error reaches up to 6.52 % if the wake shedding is ignored. So, for the lateral vibration, the wake should be considered.

3) For the first torsion natural frequency, the max relative error amounts to 5.10 % when the skew angle is 108 degrees. So, for the torsion vibration, the wake shedding should also be considered when the skew angle is large.

4) For the axial vibration, the influence of wake is little and may be ignored.