Coupled vibration characteristics of marine shafting with complex change of inclination under ice loads

Abstract

With the intensifying exploitation of the resources and the waterway of the Arctic region in worldwide, the variable inclination propulsion system are gradually used in the ice-sailing ship, its design and strength study are of great engineering significance and application value. A variable inclination propulsion shafting with universal joint is introduced and its mechanical model is established in this paper. By combine torsional stress and inclination angle calculation, a decoupling algorithm based on Newmark is designed, and then the analysis of ice loads response in both frequency domain and time domain are completed. The results show that due to the complex triangle function relation among the various components of the shafting, there is an interaction between the propeller and the transmission shafting. In particular, when considering the dynamic change of inclination, nonlinear bending-torsional coupling vibration may occur and further lead to the enhancement effect of torsional vibration. Among the three working conditions defined by DNV, the bending-torsional coupling effect vary with rotation speed and inclination angle, and is mainly concentrated in high frequency bands. In the ice loads condition 2, with the largest time domain excitation amplitude and the smallest high order component, it’s vibration response is the strongest, and the change law of dangerous phase angle is similar to the bending moment. With the increase of initial inclination angle, the amplitude of low frequency vibration decreases, yet the amplitude of high frequency vibration increase. This trend is more obvious with the increase of rotation speed, and more high frequency resonance points appear in the condition 1 and condition 3.

Keywords

Propulsion system marine shafting universal joint coupled vibration ice loads

Introduction

In recent years, some transmission devices with certain inclination angle are being widely used in ships. In a special boat structures, the propeller shaft is installed in the submerged bodies, and there is a large difference between the output center line of motor shaft and propeller shaft. Therefore, the two shafts need to be linked up by the coupling at a certain angular, which can deform greatly to prevent propulsion system from being overloaded. Besides, along with the need of path tracking attitude keeping, and capable of travelling at low speed, this device is being applied various sea vehicles. Typical cases include the semi-immersed propulsion system, the vector thruster and so on. In this structure, there are an interaction between the propeller and the transmission shafting.¹ The fluctuating force created by the loads on propeller will cause the vibration of shafting and further cause the deformation of shafting, which create a feedback loop that would increase vibration and deformation further. Therefore, the motion characteristics and load response characteristics of the propulsion shafting will be different from the traditional technical specifications.

There are some basic researches on this kind of propulsion system. Applying finite analysis software, Li² used the finite element method to complete the modal analysis of elastic coupling at rated torque, and obtain the first ten steps natural frequencies and the modes of vibration. Li³ performed the kinematic and dynamic simulation of transmission shafting with universal coupling in the right-hand coordinate system according to the multi-body theory. Xia⁴ analyzed the motion of the single universal coupling by using the UG software. Petrescu and Petrescu⁵ analyzed the geometric, structural and kinematic properties of the universal joints aimed at the drive system with universal. For better application of this kind of propulsion system, lots of scholars have studied its response characteristics under external loads and its structure optimization method. Using frequency method and time domain method, Zhao et al.⁶ obtained the stress distribution of an universal coupling under the action of impact loads, which provided a basis for shock resistance performance evaluation of the corresponding device. Gagandeep et al.⁷ established a 3D model of universal joint by the Solidworks, and used the structural steel material to reduce the stress of the universal joint under the impact loads. An and Wand⁸ analyzed the influence of phase angle of loads on the transmission performance of the universal coupling. Daskalov⁹ analyzed the reasons for the unequal speed of the driving and driven shaft, and further discovered that this phenomenon can cause fluctuations in the additional torque, resulting in forced vibration of the system.

However, when the driven system with certain inclination angle is applied to ship propeller, in the process of inclination adjustment, the driven shaft occur whirling motion along with rotating around its center line. Accordingly, the propeller is driven to make complex space rotation motion that influence its inflow characteristics, and even shafting vibration and hull posture. Therefore, chosen the rotor-coupling-bearing-vibration absorbed of marine power plant as object. Xu et al.¹⁰ established the nonlinear torsional vibration equations of universal joint, and came to a conclusion that the angle can be omitted in torsional calculation when it is less than 20 deg. Liu et al.¹¹ studied the influence of bending-torsional vibration of universal joint on the rotational speed in y and z directions of the driving shaft. Liu et al.¹² developed the theory and test research on the dynamic characteristics of shafting system and its main components, which revealed the nonlinear torsional vibration generated by the coupling joint. This phenomenon has become one of technical problems to be solved, that is particularly the case in ice-sailing ships. Owing to including multiple frequency excitation in the ice loads, the propulsion shafting torsional is a very complicated process and have high risk, which is one of the research focuses currently. Among them, Waal¹³ established the milling and impact model of propeller blades under ice loads, and obtained the laws of thrust and torque of propeller under ice loads. Baik¹⁴ studied the ship propulsion system with a six-cylinder diesel engine, analyzed the change of loads during the interaction between propeller and ices, and then obtained the transient response of propulsion shafting synthesized the excitation of diesel engine and ice loads. Persson¹⁵ introduced the GTORSI software of MAN B&W company, and analyzed the calculation of transient torsional vibration of propulsion shafting under the ice loads.

In the current research, because the vibration fluctuation of propulsion shafting is very small with small inclination, the vibration influence caused by shafting inclination is often ignored in engineering applications. Accordingly the vibration coupling effect caused by the change of the inclination is rarely studied at present. In fact, the inclination angle change complicated with the change of loads in the process of transferring torque. The inclination angle also have a very complex influence on the vibration of the shafting system, in which nonlinear vibration enhancement effect may occur. Scholars have proposed many mathematical models and test methods to address this issue.^16,17 Mazzei and Scott¹⁸ used Galerkin’s technique to obtain a set of coupled linear differential equations with time-dependent coefficients, and then some effects of internal viscous damping on parametric and flutter instability zones are investigated by the monodromy matrix technique. Yu¹⁹ tested the vibration sound radiation of the underwater vehicle under different shafting angles and pedestal installation angles. To amalgamate the torsional vibration and inclination angle calculation, these models consists of a set of coupled, linear partial differential equations with time-dependent coefficients, that is difficult to solve directly.²⁰

Therefore, a mechanical model of marine shafting with complex changes of inclination angle is established in this paper firstly. Furthermore, based on three working conditions refer to the regulations of classification societies, the torque characteristics of propeller under difference ice loads are analyzed, and then the influence of inclination angle and external loads are comprehensively considered and applied to dynamic loads response calculation. Finally, according the computational requirements of DNV, the torsional vibration response are calculated in the time domain and the frequency domain respectively, and the dangerous phase angles are analyzed under different ice loads. This study can play a guiding role in solving practical engineering problems in ships design and ships handling.

Kinematic analysis of shafting with universal joint

As shown in Figure 1 is the application of variable inclination propulsion system in an electrical propulsion ship. In general, both motor and propulsion system are mounted in the submerged body, however, due to the small size of the ship design, the motor is moved to install on the deck. They are connected through a cross shaft universal joint, which can make full use of the upper space of boat and greatly increase the flexibility of design, maintenance and overhaul of propulsion system. Besides, to enhance the robustness of navigation control and autonomy, the variable inclination propulsion, such as the vector thruster are also gradually used and its design and strength study are of great engineering significance and application value.

Figure 1.

Application of variable inclination propulsion system.

Due to the small size and high speed of shafting, two cross shaft of universal joint are taken into account to calculate the torsional vibration characteristics accurately. The equivalent model of the universal joint is established as shown in Figure 2, in which the $o X_{1} Y_{1} Z_{1}$ and $o X_{3} Y_{3} Z_{3}$ are the fixed coordinate of the driving shaft and the driven shaft respectively. The model consists of six nodes with their own rotating inertia prospectively and free degrees of angular displacement. As shown in Figure 2, $J_{1}$ , $J_{2}$ , $J_{3}$ , $J_{4}$ , $J_{5}$ , $J_{6}$ are the moment of inertia of the driving shaft, driving cross shaft, intermediate shaft, driven cross shaft and driven shaft respectively.

Figure 2.

The equivalent torsion model of shafting with universal joint.

The Lagrange method is a universal method for establishing dynamic equations of multibody systems, in which the second type of Lagrange equation is applicable to tree shaped multibody systems. So, the dynamic equation of propulsion system can be established according to the second type of Lagrange equation in this paper.

Kinetic energy of the shafting can be calculated as follows

E_{p} = (J_{1} \cdot θ_{1}^{2} + J_{2} {\overset{\cdot}{θ}}_{2^{2}} + J_{4} {\overset{\cdot}{θ}}_{4^{2}} + J_{5} {\overset{\cdot}{θ}}_{5^{2}} + J_{6} {\overset{\cdot}{θ}}_{6^{2}} + J_{7} {\overset{\cdot}{θ}}_{7^{2}}) / 2

(1)

Potential energy of the shafting can be calculated as follows.

\begin{array}{l} E_{v} = [c_{1} {(θ_{1} - θ_{2 x})}^{2} + c_{2 x} {(θ_{2 x} - θ_{4} \cos r)}^{2} + c_{2 y} {(θ_{2 y} - θ_{4} \sin r)}^{2} \\ + . c_{3} {(θ_{4} - θ_{5})}^{2} + c_{4 y} {(θ_{5} \cos r - θ_{6 x})}^{2} + c_{4 x} {(θ_{5} \sin r - θ_{6 y})}^{2} \\ + c_{5} {(θ_{6 x} - θ_{7})}^{2}] / 2 \end{array}

(2)

Where in: $θ$ are the angular displacement, and in its lower corner mark, the numbers and characters are the symbols corresponding to the part number of Figure 2 and the coordinate axis. $r$ is the inclination angle. $C_{3}$ is the torsional stiffness of the middle shaft; $C_{1}$ , $C_{2 x}$ , $C_{2 y}$ , $C_{4 x}$ , $C_{4 y}$ , $C_{5}$ are the contact stiffness between the cross shaft and the driving or driven shaft; $C_{2 z}$ , $C_{5 z}$ are the rotational stiffness of the cross shaft, which are negligible compared with the contact stiffness.

The kinetic energy of the cross shaft can be expressed as follows.

J_{6} {\overset{\cdot}{θ}}_{6^{2}} = (J_{6 x} {\overset{\cdot}{θ}}_{{6 x}^{2}} + J_{6 y} {\overset{\cdot}{θ}}_{{6 y}^{2}} + J_{6 z} {\overset{\cdot}{θ}}_{{6 z}^{2}})

(3)

According to the symmetry of the cross shaft as shown in equation (4) and the motion equation of universal joint, the torsional angular velocity of cross shaft can be obtained as equation (5).

J_{6 y} = J_{6 z}

(4)

{\begin{matrix} {\dot{θ}}_{6 y} = \frac{\sin 2 r \cdot {(\cos θ_{6 x})}^{2}}{1 - {(\sin r \cdot \cos θ_{6 x})}^{2}} \cdot {\dot{θ}}_{6 x} \\ {\dot{θ}}_{6 z} = \frac{\sin 2 r \cdot \sin θ_{6 x} \cdot \cos θ_{6 x}}{1 - {(\sin r \cdot \cos θ_{6 x})}^{2}} \cdot {\dot{θ}}_{6 x} \end{matrix}

(5)

Put equations (4) and (5) into the equation (3), which can be expressed as equation (6).

J_{6} {\dot{θ}}_{6^{2}} = [J_{6 x} + J_{6 y} {(\frac{\sin 2 r \cdot \cos θ_{6 x}}{1 - {(\sin r \cdot \cos θ_{6 x})}^{2}})}^{2}] \cdot {\dot{θ}}_{6 x}^{2}

(6)

The equation (7) can be obtained by integrating two sides of equation (6), and the same calculation for $θ_{2 y}$ , $θ_{2 z}$ can be performed.

{\begin{matrix} θ_{6 y} = \frac{2}{\sin r} \cdot \arctan (\frac{\tan θ_{6 x}}{\cos r}) - \frac{2 θ_{6 x}}{\tan r} \\ θ_{6 z} = - \frac{1}{\tan r} \cdot \ln ({(\sin r \cdot \cos θ_{6 x})}^{2} - 1) \end{matrix}

(7)

$θ_{1}$ , $θ_{2 x}$ , $θ_{4}$ , $θ_{5}$ , $θ_{6 x}$ and $θ_{7}$ are selected as generalized coordinates, the differential equations of vibration system can be established as equation (8):

{\begin{array}{l} J_{1} \dot{θ_{1}} + C_{1} (θ_{1} - θ_{2 x}) = 0 \\ [J_{2 x} + J_{2 y} {(\frac{\sin 2 rcos θ_{2 x}}{1 - {(sinr)}^{2} {(\cos θ_{2 x})}^{2}})}^{2}] \dot{θ_{2 x}} - C_{1} (θ_{1} - θ_{2 x}) + C_{2 x} (θ_{2 x} - θ_{4} cosr) \\ + C_{2 y} [\frac{2}{sinr} \arctan (\frac{\tan θ_{2 x}}{cosr}) - \frac{2 θ_{2 x}}{tanr} - θ_{4} sinr] \frac{\sin 2 r {(\cos θ_{2 x})}^{2}}{1 - {(sinr)}^{2 *} {(\cos θ_{2 x})}^{2}} = 0 \\ J_{3} \dot{θ_{4}^{2}} + C_{3} (θ_{4} - θ_{5}) - C_{2 x} cosr (θ_{2 x} - θ_{4} cosr) - C_{2 y} sinr [\frac{2}{sinr} \arctan (\frac{\tan θ_{2 x}}{cosr}) \\ - \frac{2 θ_{2 x}}{tanr} - θ_{4} sinr] = 0 \\ J_{4} \dot{θ_{5}^{2}} - C_{3} (θ_{4} - θ_{5}) + C_{4 y} cosr (θ_{5} cosr - θ_{6 x}) + C_{4 x} sinr [- \frac{2}{sinr} \arctan (\frac{\tan θ_{6 x}}{cosr}) \\ + \frac{2 θ_{6 x}}{tanr} + θ_{4} sinr] = 0 \\ [J_{5 x} + J_{5 y} {(\frac{\sin 2 ros θ_{6 x}}{1 - {(\sin r)}^{2} * {(\cos θ_{6 x})}^{2}})}^{2}] θ_{6 x} - C_{4 y} (θ_{5} cosr - θ_{6 x}) + C_{5} (θ_{6 x} - θ_{7}) \\ - C_{4 x} [- \frac{2}{\sin r} \arctan (\frac{\tan θ_{6 x}}{cosr}) + \frac{2 θ_{6 x}}{tanr} + θ_{5} sinr] \frac{\sin 2 r {(\cos θ_{6 x})}^{2}}{1 - {(\sin r)}^{2} {(\cos θ_{6 x})}^{2}} = 0 \\ J_{6} \dot{θ_{7}} + C_{5} (θ_{6 x} - θ_{7}) = 0 \end{array}

(8)

Initial parameters are set as: $J_{1} = J_{4} = 4.77 kg . m^{2}$ , $J_{2} = J_{3} = 8.50 kg . m^{2}$ . Geometrical parameters of intermediate shaft $d = 0.3 m, l = 1.2 m$ ; Both the driving shaft and the driven shaft are $d = 0.2 m, l = 0.6 m$ ; $J_{2 x} = J_{5 x} = 1.9 kg . m^{2}$ ; $J_{2 y} = J_{5 y} = J_{2 z} = J_{5 z} = 1.14 kg . m^{2}$ ; All contact stiffness are $2 \times 10^{10} Nm / rad$ . The equation (8) contain a large number of complex trigonometric functions and is difficult to obtain its analysis solution. Take $r = 1 \deg$ and numerical simulation is conducted, the angular velocity curve of the driven shaft can be obtained as shown in Figure 3.

Figure 3.

Time domain diagram of angular velocity. (a) Time domain of 0–0.01 s and(b) time domain of 0–0.05 s.

Figure 3(a) is the time domain diagram of the first 0.01 s, it can be seen that due to there are complex trigonometric functions between the motion characteristics of the driving and the driven shaft, the angular velocity has fluctuation according to time while the overall change trend remains unchanged. In Figure 3(b), from 0.01 to 0.03 s the curve began to show its various facets, but the change law is similar to Figure 3(a). The fluctuation of angular velocity is further enhanced and gradually display the change law of quasi-periodic. So it can be considered that the low-frequency fluctuation is determined by the inertia and stiffness characteristics of the driving shaft, intermediate shaft and driven shaft. Instead, the performances in the high frequency band are determined by the inertia and stiffness characteristics of the cross shaft, in which the additional elastic torsional vibration is easy to form under external loads.

Bending-torsional coupling vibration under ice-loads

Torque characteristics of propeller under ice-loads

The ice loads are the main excitation of propeller when ships operating in ice area that can be calculated refer to the regulations of classification societies such as ABS, CCS, DNV.²¹ According to the computing formula of the ice loads provided by DNV, the process of collision between propeller blades and ices can be divided into three stages as shown in Figure 4: (1) contact phase; (2) penetration phase; (3) separation phase. Based on these process, the torque excitation of a single blade under ice loads are divided into three working conditions as shown in Table 1, and the torque function can be calculated by equation (9).

{\begin{matrix} T_{t} = C_{b} . T_{max} . \sin [θ (180 / α_{n})], 0 \leq θ \leq a_{n} \\ T_{t} = 0, α_{n} \leq θ \leq 360 \end{matrix}

(9)

Where in: $θ$ is the rotation angle of blades; $C_{b}$ is constant value that correlates with $α_{n}$ as shown in Table 1. $T_{max}$ is the maximum torque of the propeller under ice loads, taking the four-blade propeller as the object and set $D \leq D_{lim it}$ , that can be calculated by equation (10)²² :

T_{max} = 10.9 - (1 - \frac{d}{D}) {(\frac{p_{0.7}}{D})}^{0.16} \overset{0.17}{(nD)} D^{3}

(10)

Where in: $d$ is the external diameter of propeller boss, and is set as 0.33 m; $D$ is diameter of propeller, and is set as 0.6 m; $D_{lim it}$ is the critical value of propeller diameter, and $D_{lim it} = 1.8 H_{ice}$ ; $p_{0.7}$ is the pitch at 0.7R, and is set as 0.5 m; $n$ is revolution speed of propeller.

Figure 4.

The interaction between blades and ices.

Table 1.

Parameters of three working conditions.

Conditions	Characteristic of ices	$α_{n}$	$C_{b}$
CASE1	single ice	90°	0.75
CASE2	single ice	135°	1.0
CASE3	two ices (phase difference 45°)	45°	0.5

According the definition of ice loads in classification society, the excitation torque is treated as a linear ramp function in 270 degrees rotation angle scope from integration to separation between the blades and ices, and then the time domain curves of excitation torque can be fitted as shown in Figure 5(a), (c), and (e). After harmonic analysis, the excitation torque of three working conditions are transformed into a series of frequency components²³ as shown in Figure 5(b), (d), and (f).

Figure 5.

Torque characteristics under different ice load conditions. (a) Torque at condition 1, (b) inclination angle at condition 1, (c) torque at condition 2, (d) inclination angle at condition 2, (e) torque at condition 3, and (f) inclination angle at condition 3.

It can be found that the instantaneous torque acting on the propeller shaft are almost the same as a series of semi-sinusoidal waves in all three conditions, and the difference mainly lie in the peak size. The maximum instantaneous torque of CASE1 is 11.6 kN.m, the effective component is 7.1 kN.m, the main order are 4, 8, 12, …; The maximum instantaneous torque of CASE2 is 19.8 kN.m, the effective component is 13.5 kN.m, the maximum instantaneous torque of CASE3 is 12 kN.m, the effective component is 6.4 kN.m, the main order of both of CASE2 and of CASE3 are 8, 16, 24, …. In the CASE2, the amplitude in time domain is the biggest, but the amplitude of each order is the minimum after the frequency-domain transformation. The fourth order of CASE1, and the eighth order CASE3 are large that needs to be focused on when calculating the torsional vibration.

Estimation of bending moment and its influence on inclination angle

When ships navigate near surface, the pulse power caused by ice impact on propeller will cause the driven shaft to be subjected to bending moment. And then the inclination angle is changed, that will cause the angular displacement and angular velocity of the driven shaft to change, and lead to the change of torsional vibration characteristics, that is, the bending-torsional coupling vibration.

According to the analysis of literature,²⁴ when the driving shaft and the driven shaft are parallel, the torque acting on the intermediate shaft is 0. In addition, the intermediate shaft usually has the characteristic of scalability, that is, it is without axial vibration and lateral vibration. By set the bending moment of driven and driving shaft are M1 and M2 respectively and expand along the X, Y, and Z directions, the equation (11) can be obtained²⁵:

{\begin{matrix} M_{2} \cos θ_{3} - M_{1} \sin θ_{1} = 0 \\ M_{1} \cos θ_{1} + M_{2} \sin θ_{3} \cos r + T_{2} \sin r = 0 \\ - T_{1} - M_{2} \sin θ_{3} \sin r + T_{2} \cos r = 0 \end{matrix}

(11)

Furthermore, according to motion characteristics of universal joint: $\tan θ_{1} = \cos r' \tan θ_{3}$ , the total bending moment on the driven shaft under ice loads can be expressed as equation (12).

M_{1} = - \cos θ_{1} \tan r \cdot T_{1}

(12)

To obtain the influence of bending moment on the torsional vibration, its influence on the inclination angle must be calculated firstly. The motion model of universal joint is established as shown in Figure 6, the solid line represents the initial position of the intermediate shaft, and the dotted line represents its position under the additional bending moment $M_{1}$ , in which the inclination angle changes from $r$ to $r'$ .

Figure 6.

The change of inclination angle under bending moment.

According to the mechanism of movement as shown in Figure 6 and the difference in speed between the driving and the driven shaft, equation (13) can be listed.

\frac{H - \frac{M \sin ({\overset{\cdot}{θ}}_{1} tL)}{EI}}{H} = \frac{\sin r'}{\sin r}

(13)

In order to observe the variation relationship of angular velocity between the driven shaft and the cross shaft, according to the motion equation of the rigid body of rotation about a fixed point and the symmetry of cross shaft, the equation (14) can be listed.

{\overset{\cdot}{θ}}_{2 y} = 2 \overset{\cdot}{α} \cos θ_{1} = \frac{\sin 2 r \cdot {(\cos θ_{1})}^{2}}{1 - {(\sin r \cdot \cos θ_{1})}^{2}} \cdot {\overset{\cdot}{θ}}_{1}

(14)

Where in, $α$ is rotation angle of coordinate system $o X_{1} Y_{1} Z_{1}$ . Put equation (14) in equation (13), the expression for $r'$ can be concluded as equation (15).

\frac{H - \frac{M \sin ({\overset{\cdot}{θ}}_{2 y} \cdot \frac{1 - {(\sin r \cdot \cos θ_{1})}^{2}}{\sin 2 r \cdot {(\cos θ_{1})}^{2}} \cdot t \cdot L)}{EI}}{H} = \frac{\sin r'}{\sin r}

(15)

Taking the inclination angle 3, 5, 8 degrees respectively, the bending moment generated at the propeller at no-loads state can be made as shown in Figure 7(a). It can be found that the change of bending moment with rotation angle are the same as those of torque excitation.²⁶ And then, according to three ice load conditions, the inclination angle changing with time are calculated by equation (15) as shown in Figure 7(b). It can be found that the total bending moment acting on the shaft is different in different ice loads conditions, which correspondingly leads to complex change of inclination angle, even the driving shaft rotates at a constant angular velocity. Compared with condition 1 and 3, the bending moment has a more influence on the inclination angle in condition 2. However, this kind of fluctuation is often ignored in engineering applications due to its small scope, that will be calculated for improve the accuracy of vibration response in this paper.

Figure 7.

Change of bending moment and inclination angle. (a) Bending moment and (b) inclination angle under ice loads.

Calculation process based on decoupling algorithm

According to the equivalent torsion model of shafting in chapter 1, when calculating the bending-torsional coupling vibration under the ice loads, the vibration characteristic of shafting and the dynamic change of inclination angle should be calculated simultaneously. But the two calculation process are coupled to each other, it is difficult to decouple directly.²⁷ Therefore, we propose a dynamic decoupling algorithm, in which the two calculation process are performed separately in the time domain, and then combine them by using Newmark method to guarantee convergence. The iterative steps of this algorithm are described as follows:

1) Assign each parameter of the shafting: set the initial inclination angle $r'$ , the time t = 0; and the angular acceleration and angular displacement of the driving shaft, driven shaft and intermediate shaft are all 0.

2) Calculate bending moment according to the ice loads.

3) set $t = t_{0} + Δ t$ , in which the current value of bending moment, angular acceleration are updated and put in equation (8), so the inclination angle also is obtained by $r = r' + Δ r$ .

4) Suppose the acceleration varies linearly within the time $[t, t + Δ t]$ , the state parameters $(θ, \overset{\cdot}{θ})$ can be established as follows:

{\begin{matrix} {\overset{\cdot}{θ}}_{t + Δ t} = {\overset{\cdot}{θ}}_{t} + [(1 - δ) {\overset{\cdot\cdot}{θ}}_{t} + δ {\overset{\cdot\cdot}{θ}}_{t + Δ t}] Δ t \\ θ_{t + Δ t} = θ_{t} + {\overset{\cdot}{θ}}_{t} Δ t + [(1 / 2 - β) {\overset{\cdot\cdot}{θ}}_{t} + β {\overset{\cdot\cdot}{θ}}_{t + Δ t}] Δ t^{2} \end{matrix}

(16)

After variable substitution according to Newmrk method, equation (16) becomes equation (17)

{\begin{matrix} {\overset{\cdot\cdot}{θ}}_{t + Δ t} = \frac{1}{β Δ t^{2}} ({θ}_{t + Δ t} - {θ}_{t}) - \frac{1}{β Δ t} {\overset{\cdot}{θ}}_{t} - (\frac{1}{2 β} - 1) Δ t {\overset{\cdot\cdot}{θ}}_{t} \\ {\overset{\cdot}{θ}}_{t + Δ t} = \frac{δ}{β Δ t} (θ_{t + Δ t} - θ_{t}) + (1 - \frac{δ}{β}) {\overset{\cdot}{θ}}_{t} + (1 - \frac{δ}{2 β}) Δ t {\overset{\cdot\cdot}{θ}}_{t} \end{matrix}

(17)

Put equation (17) in equation (8)

\begin{matrix} [M] {{\overset{\cdot}{θ}}_{t} + [(1 - δ) {\overset{\cdot\cdot}{θ}}_{t} + δ {\overset{\cdot\cdot}{θ}}_{t + Δ t}] Δ t} + [K] θ_{t + Δ t} + [C] \\ {\frac{δ}{β Δ t} (θ_{t + Δ t} - θ_{t}) + (1 - \frac{δ}{β}) {\overset{\cdot}{θ}}_{t} + (1 - \frac{δ}{2 β}) Δ t {\overset{\cdot\cdot}{θ}}_{t}} \\ = - [M] {p} {\overset{\cdot\cdot}{u}}_{t + Δ t} \end{matrix}

(18)

So, $θ_{t + Δ t}$ can be obtained by solving equation (18).

5) Repeat the steps 2–4 to obtain the new state parameters $(θ, \overset{\cdot}{θ})$ after the interval $Δ τ$ until the convergence of phase trajectory is realized, which in turn are used to calculate other vibration characteristics.

In order to verify the feasibility of the above calculation method, set $ω = 3600 \deg / s$ , $r = 1 \deg$ , set the relative error of calculation is 1e-6, 1e-9 respectively, the angular displacement of the driven shaft in the time domain calculated by our algorithm and ODE45 algorithm can be obtained as Figure 8(a) and (b). It can be seen that with reduce of the relative error, the two curves gradually coincide with each other, especially the deviation of peak value greatly reduce as shown in the A and B while the calculation time of the ODE45 algorithm greatly increase with the decrease of relative error. So our algorithm can improve the calculation efficiency on the condition of calculation accuracy regardless of the load characteristics, and can also be applied to other calculation of vibration characteristics of propulsion shafting.

Figure 8.

Computation comparison between two algorithms. (a) Relative error 1e-6 and (b) relative error 1e-9.

To observe the bending-torsional coupling effect by this decoupling algorithm, it is necessary to measure the vibration characteristics of the no-load state firstly. So set the initial inclination angle to 3 degrees, the angular displacement changing with time in two conditions considering and not considering the dynamic change of inclination angle are obtained as Figure 9(a). Further, three frequency bands are divided, that is, low frequency band (0–150 Hz), intermediate frequency band (150–500 Hz) and high frequency band (over 500 Hz), and the amplitude of the three frequency bands under different inclination angles are calculated as shown in Figure 9(b).

Figure 9.

Angular displacement considering and not considering variation of inclination. (a) Time domain at two condition and (b) amplitude under different inclination.

In Figure 9(a), there are certain error existence at value between two conditions, in which the angular displacement considering the dynamic change of inclination angle is slightly smaller than that not considering. This trend is more obvious after 0.03 s, in which the third and fourth peaks of each cycle increase as the excitation continued. In Figure 9(b), with the inclined angle increasing, the change rule under two conditions of three frequency bands are basically unchanged, but the deviation of high-frequency band is more remarkable than of low-frequency band. This result is consistent with the Figure 9(a). The influences of the inclination angle on the torsional vibration are mainly concentrated in the high frequency region. So although the influence of inclination angle can be ignored in the allowable error range of engineering, the high order excitation of external loads can not be negligible, especially in the case of the inclination angle is above 8 degrees.

Vibration characteristics under action of ice loads

Spectrum analysis under different ice loads

According to the above calculation process, the torsional vibration response under ice loads are calculated and analyzed in the time domain and the frequency domain respectively below. Since the torque of propeller is proportional to the square of rotation speed, the ice loads excitation can be calculated by the power and rotation speed of propeller, and ignores the hydrodynamic torque caused by dynamic speed drop generated by the ice loads. Choosing the conditions of CASE1 and making calculation simulation, the vibration response frequency spectrum are obtained by fast Fourier transform as shown in Figure 10(a). Furthermore, according to the rule of frequency division, the vibration response frequency spectrum in different ranges are magnified and displayed as shown in Figure 10(b), (c), and (d).

Figure 10.

Spectrogram of angular displacement response of driven shaft. (a) Full band, (b) low frequency band, (c) medium frequency band, and (d) high frequency band.

In the low frequency band of Figure 10(b), the driven shaft angular displacement amplitude had a decreasing trend with the growth of inclination angle. In the condition of inclination angle 3 degrees, there are two obvious peaks, but with the increase of inclination angle, the two peaks gradually decrease and approach, showing a trend of merging, and the corresponding frequency move to the right. This can be seen as the inertia of the driving shaft, intermediate shaft and driven shaft are effected by the increase of inclination angle.

In the middle frequency band of Figure 10(c), due to the relatively small moment of inertia of the cross shaft, the angular displacement vary little. However, in the high frequency band of Figure 10(d), the angular displacement fluctuations become more volatile with the increase of inclination angle. Meanwhile, due to the principle of conservation of energy, the corresponding spectrum will shift to the left. These phenomenon are caused by the complex trigonometric relationship between the driven and the driving shaft become stronger and stronger, which are the origin of the influence of the bending-torsional coupling on propulsion shafting. So the influence of the increased inclination angle on angular displacement are mainly concentrated in the high frequency band,²⁸ which will be emphatically analyzed in detail farther below.

Torsional vibration calculation under different ice loads

In engineering practice, although the calculated value of torque are within the normal range in frequency domain, it may exceed the instantaneous allowable value at a certain moment.²⁹ So the calculations of transient torsional vibration induced by ice loads have become mandatory since January 2011 in DNV specifications. According to three different conditions as shown in Table 1, the transient torque of driven shaft at speed $ω = 3600 \deg / s$ with inclination angles of 1 deg are calculated, in which the calculation time step is taken as that the propeller take to rotate 1 deg, and the loads are applied at the 720th time steps. After a total of 4000 sub-steps are calculated, the time-domain curves of torque can be obtained as shown in Figure 11(a), (c), and (e). Furthermore, the multiple speeds with inclination angles of 3,5,8 degrees under each condition are selected for time domain calculation, and the maximum amplitude of them are extracted for comparison as shown in Figure 11(b), (d), and (f) accordingly.

Figure 11.

Torsional amplitude at different conditions. (a) Condition 1 under rated speed, (b) condition 1 under all rotation speeds, (c) condition 2 under rated speed, (d) condition 2 under all rotation speeds, (e) condition 3 under rated speed, and (f) condition 3 under all rotation speeds.

By comparing the transient torque under three different ice loads conditions, it is found that the condition 2 has the strongest response with the condition 1 the second and the condition 3 the weakest. For example, at the same inclination angle of 3 degrees, the torque amplitude is 30 kN.m and the duration time is 8 s in Figure 11(c), that is average 23 kN.m and 5 s respectively in Figure 11(a), and is 34 kN.m but discontinue respectively in Figure 11(d). From the point of view of waveforms of the condition 2 and 3, the transient torque reach the maximum when the propeller get in touch with an ice, but in condition 1, the maximum transient torque generated is delayed until the propeller contacts the ice for a period time.

As shown in Figure 11(d), (e), and (f), it can be found a strong resonance point near the speed of 430 r/min, and two resonance point at 250 and 650 r/min with small amplitudes in all conditions, that is consistent with the variation trend of each frequency band shown in Figure 8. However, there are some resonance points with small amplitude near the speed of 1000 r/min in condition 3, that can be thought to result from the high order components of ice loads as shown in Figure 5. On the other hand, with the increasing rotation speed, the vibration amplitude changes more obviously with the inclination angle in all working conditions. It shows that the higher the excitation frequency is, the stronger the effect of bending-torsional coupling caused by the inclination angle of shafting can be, even the torque may exceed the instantaneous allowable value at a certain moment. That is a important factor in the reliability analysis of marine shafting with the universal joint in engineering practice.

Dangerous phase angle analysis under different ice loads

On the demands of the DNV, the time-domain calculation of ice loads must be performed in the worst phase angle condition. As described in the above vibration response calculation, the ice loads are all applied at the 720th load step, that is, the phase angles between the bending moment and ice loads excitation are considered to be zero. However, the bending moments are affected by the rotation angle and the inclination angle as can be seen from equation (12), when the inclination angle is unchanged, there is an arigonometric function relationship between the bending moment and the rotation angle. But if taking into account the dynamic changes of inclination angle, the function relationship become irregular which result in the initial phase angle between the bending moment and ice loads excitation has a great influence on the calculation results.

Due to the complex trigonometric relationship between the driven and driving shaft, it is difficult to directly calculate the phase angle of ice loads corresponding to the worst effects by analytic method, so the torsional vibration response with phase angle of ice loads from 0 to 359 degrees at intervals of 1 deg is calculated respectively, that is, the ice loads are applied at every one load step from 720th to 1079th load steps. The maximum torsional stress of the driven shaft can be obtained as shown in Figure 12, and conclusions can be are summarized as follows.

Figure 12.

Dangerous phase angle under different ice loads conditions. (a) Ice load condition 1, (b) ice load condition 2, and (c) ice load condition 3.

The most dangerous phase angles are different in different ice load conditions. It can be seen that the amplitude-phase curve of condition 1 and condition 3 are similar and irregular, but the change law of condition 2 is relatively regular and consistent with the change law of bending moment, which are related to the large amplitude of primary harmonic and the small component of high harmonic in condition 2. As can be seen from Figure 12, at the same inclination angle of 3 degrees, the most dangerous phase angle is 275, 97, and 105 degrees in condition 1, 2, 3 respectively, in which condition 3 is 170 degrees different from condition 1.

With the increase of inclination angle, the amplitude-phase curve of all conditions change, which are related to the change of natural frequency of shafting, and further lead to change of the bending-torsional coupling vibration under ice loads. Compared with ice loads condition 1 and condition 2, the amplitude-phase angle curve of condition 3 have a wider variation range with increasing inclination angle. As shown in Figure 12(c), there is 90 degrees different in the most dangerous phase angle between the inclination angle 3 and 8 degrees. Due to this complexity, experimental studies on relationship between the phase angle and the inclination angle are necessary. Although it is feasible and reliable to conduct ice loading experiments based on the propeller ice cutting process, there may be many shortcomings in measurement accuracy, and further improvement is needed.

Conclusion

Marine shafting with universal joint shows complex nonlinear dynamics properties due to the trigonometric function relationship between the driving and driven shaft. To obtain the accurately vibration characteristics of this shafting under ice loads, the dynamic change of inclination angle must be considered to calculate the bending-torsional coupling vibration. By using the Newmark method to amalgamate the torsional vibration and inclination angle calculation, the results of the numerical simulation are summarized as follows:

(1) The torsional fluctuating caused by ice loads on the propeller will cause driven shaft subject additional bending moment, that can induces complex change of inclination angle and the further strengthening of bending-torsional coupling vibration. Meanwhile the influencing extent of bending moment on inclination angle is different in different frequency bands, that in the high frequency range is greater than in the low frequency range, which may play a role in the occurrence of vibration enhancement effect. Due to including multiple frequency component in ice loads, although the calculated value of torque are within the normal range in frequency domain, it may exceed the instantaneous allowable value at a certain moment.

(2) Due to there are strong nonlinear characteristic in vibration response of the universal joint, when calculating bending-torsional coupling vibration under the ice loads, the shafting system torsion vibration characteristic and the dynamic change of inclination angle should be calculated simultaneously. The dynamic decoupling algorithm based on the Newmark method can be well applied for ice loads response calculation, in which can improve the calculation efficiency on the condition of calculation accuracy, and can also be applied to the design of propulsion shafting with universal joint under other loads.

(3) The torsional vibration response of the propulsion shafting shows different changes with the different ice loads conditions and operating conditions. According to the regulations of the classification society on ice-sailing ships, the torque excitation of a single blade under ice loads are divided into three working conditions, in which have similar time waveform but defense frequency spectrum. The propulsion system has the strongest vibration response in the condition 2 with the second in the condition 1 and the weakest in the condition 3. Moreover, in higher rotate speed range, there are additional resonance point with small amplitude in the condition 3. The dangerous phase angle is not equal under different working conditions. There is a wider variation rang under working condition 3, which is something to consider in the reliability analysis of the propulsion shafting. Therefore, when sailing in the ice area, the operating condition of propulsion shafting should be reasonably selected according to the different loads.

Footnotes

Handling Editor: Chenhui Liang

Authorship contributions

The first author: Conceptualization, Methodology, Software, Writing Original Draft. The second author: Formal analysis, Supervision, Writing-Review & Editing.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by Zhejiang Provincial Natural Science Foundation of China under Grant No. LY20E090002.

Ethical approval

This article does not contain any studies with human participants or animals performed by any of the authors.

Informed consent

Informed consent was obtained from all individual participants included in the study.

ORCID iDs

Liangxiong Dong

Jubao Li

Data availability

The data used to support the findings of this study are available.

References

Evernden

HIF

. The propeller shaft or hooke’s coupling and the cardan joint. Proc Inst Mech Eng Autom Div 2006; 2: 100–110.

. Modification of the speed ratio formula of universal coupling and the calculation formula of belt drive pressure. World Sci Techn Res Dev 2008: 30; 38–41.

Wang

Lou

, et al. Study on rigid-flexible coupling dynamics of double cross joint. Mech Trans 2016; 40: 33–37.

Xia

. Motion analysis and simulation of universal coupling based on UG software. Equip Manu Techn 2013; 4: 45–46.

Petrescu

. The structure, geometry, and kinematics of a universal joint. Indep J Manag Prod 2019; 10: 1713.

Zhao

Liu

Chen

, et al. Experimental study on misalignment and frequency locking trouble of rotor-coupling-bearing-isolator system. J Dyn Control 2009; 7: 171–176.

Gagandeep

Singh

. Finite element methods based structural analysis of universal joint using structural steel material. Indian J Sci Techn 2017; 10: 1–7.

Wang

. Transmission performance and fault analysis of a vehicle universal joint. Adv Mech Eng 2017; 9: 1–10.

Daskalov

. Kinematic analysis of cardan drives. Mech Mach Theory 1990; 25: 479–486.

10.

Zhou

Wang

. Simulation of nonlinear torsional vibration of universal joint based on MATLAB. J Wuhan Univ Techn 2011; 33: 135–138.

11.

Liu

Dong

Gan

. Kinematic characteristics analysis of an unmanned boat shafting based on universal joint. Ship Sci Techn 2023; 45: 80–85.

12.

Liu

Bose

Veitch

. Evaluation, design and optimization for strength and integrity of polar class propellers. Cold Reg Sci Techn 2015; 113: 10421–10425.

13.

Waal

. An investigation of shaft line torsional vibration during ice impacts on PSRVs. Stellenbosch University, Western Cape, 2017.

14.

Baik

. The study for stress calculation of slip damage between propeller boss and shaft on the large vessel. J Korean Soc Mar Environ Safety 2011; 17: 291–294.

15.

Persson

. Ice impact simulation for propulsion machinery. MTZ Ind 2015; 5: 34–41.

16.

Masoud

Ghazavi

Asghar

. Parametric resonances for torsional vibration of excited rotating machineries with nonconstant velocity joints. J Vib Control 2018; 24: 3262–3277.

17.

Huang

. Motion simulation of double universal coupling. J Donghua Univ 2011; 37: 506–508.

18.

Mazzei

Scott

. Effects of internal viscous damping on the stability of a rotating shaft driven through a universal joint. J Sound Vib 2003; 265: 863–885.

19.

. Influence of shaft angle of cardan joint on acoustic radiation characteristics of underwater vehicle. Digit Ocean Underw Warfare 2021; 4: 395–400.

20.

Zhou

, et al. Coupling matching research and torsional vibration characteristics prediction of MTU diesel engine propulsion shafting. China Ship Res 2019; 14: 120–127.

21.

DNV. Ice strengthening of propulsion machinery, Classification Notes No.51.1, DNV (January 2011).

22.

Xiao

Zhou

. Calculation of transient torsional vibration of ship shafting in ice area and software development. China Shipbuild 2019; 60: 138–149.

23.

Barro

Lee

. Transient torsional vibration analysis of ice-class propulsion shafting system driven by electric motor. Trans Korean Soc Noise Vib Eng 2014; 24: 667–674.

24.

Vijayan

Barik

Sha

. Shock transmission through universal joint of cutter suction dredger. Ocean Eng 2021; 233: 109185.

25.

Cui

Zhou

Wand

. Research on shock resistance of ship propulsion shafting with flexible coupling. Shipbuild China 2023; 64: 61–72.

26.

Batrak

Serdjuchenko

Tarasenko

. Calculation of torsional vibration response propulsion shafting system caused by ice impacts. Proc 1st Torsi Vibra Sympo 2014; 999: 1–8.

27.

Wang

. Transmission performance and fault analysis of a vehicle universal joint. Adv Mech Eng 2017; 9: 1–10.

28.

Wen

Cai

Wang

. Numerical simulation of ship propulsion shafting under cyclic load. China Navig 2021; 44: 13–19.

29.

Barro

Eom

Lee

. Transient torsional vibration response due to ice impact torque excitation on Marine diesel engine propulsion shafting. Trans Korean Soc Noise Vib Eng 2015; 25: 321–328.