# Dynamic Instability Analysis of a Rotating Ship Shaft under a Periodic Axial Force by Discrete Singular Convolution

Authors: Wei Li; Zhiwei Song; Xuexia Gao; Zhigang Chen
Journal: Shock and Vibration
 (2015)

---

## Abstract

Dynamic instability of a rotating ship shaft subjected to a periodic axial force is studied by using discrete singular convolution (DSC) with regularized Shannon’s delta kernel. The excitation frequency is related to the spinning speed and the number of blades on the propeller. Effects of number of blades, constant term in the periodic force, and damping on dynamic instability regions are investigated. The results have shown that the increase of number of blades and damping could improve the dynamic stability of rotating shaft with damping. The increase of constant term in the periodic force leads to dynamic instability regions shifting to lower frequencies, making the shaft more sensitive to periodic force. Those dynamic instability regions obtained by DSC method have been compared with those by Floquet’s method to verify the application of DSC method to dynamic instability analysis of rotating ship shaft.


---

## Body

## 1. Introduction

When a ship sails at sea, the propulsion shaft is subjected to a periodic axial force induced by a rotating screw propeller. The excitation frequency of the periodic axial force has relation with spinning speed and the number of blades on the propeller. When the excitation frequency and the natural frequency of the shaft are in a certain condition, that may cause serious dynamic instability phenomenon [1]. In order to avoid this phenomenon happening and protect the safety of the rotating shaft, much attention should be paid to the dynamic instability analysis of such rotating ship shafts.The dynamic instability analysis of rotating systems subjected to axial forces has been studied for many years. Chen and Ku [2, 3] studied the dynamic stability of a rotating Timoshenko shaft and a rotating Timoshenko shaft-disk system, respectively. Lee [4, 5] analyzed the dynamic stability of spinning pretwisted cantilever beams subjected to axial pulsating loads and axial base excitations, respectively. Then Lee et al. [6] studied dynamic stability of a spinning Timoshenko shaft with time-dependent spin rate by the use of the multiple scales method. And also Tan et al. [7] analyzed parametric instability of spinning pretwisted beams subjected to spin speed perturbation. Liao and Huang [8, 9] addressed parametric instability of a spinning pretwisted cantilever beam under periodic axial force and a cantilever pretwisted beam with time-dependent spinning rate, respectively. Chen and Peng [10] investigated the dynamic stability of a rotating composite shaft under axial periodic forces. Sheu and Chen [11] used a lumped mass model to study the parametric instability of a cantilever shaft-disk system subjected to axial and follower loads. Young and Gau [12, 13] studied the dynamic stability of a pretwisted cantilever beam with constant spin rates and nonconstant spin rates under axial random forces at the free end. Lin and Chen [14] presented dynamic stability problems of spinning pretwisted sandwich beams with a constrained damping layer subjected to periodic axial loads. Recently, Pei [15] analyzed the differences in the instability regions for gyroscopic system by Bolotin’s method and Floquet’s method. Pavlović et al. [16] studied the stochastic stability of an elastic balanced rotating shaft subjected to action of axial forces at the ends.In these published works, the periodic axial forces are the typeP(t)=PS+PDcos⁡(θt) (PS is the constant part, PD is the amplitude of dynamic part with time t, and θ is the excitation frequency.) But, for the rotating ship shaft [17], the periodic axial forces made by rotating propeller can be expressed in a series; the type is (t)=PS+PD1cos⁡(ZpΩt)+PD2cos⁡(2ZpΩt)+PD3cos⁡(3ZpΩt)+⋯. The excitation frequencies of periodic axial forces are equal to the blade frequency ZpΩ and its multiplies (iZpΩ,i=2,3,…). The blade frequency is equal to the shaft spinning speed Ω multiplied by the number of blades Zp. PS is the constant term; PD1, PD2, and PD3 are the amplitudes of dynamic terms. Compared to PS and PD1, the values of PD2 and PD3 are very small and could be neglected. So the axial force considered in this work is expressed as P(t)=PS+PDcos⁡(ZpΩt). Therefore, it is necessary to analyze the dynamic instability problem of a rotating ship shaft subjected to such type of periodic axial force. In this work, a relatively new method, discrete singular convolution (DSC) [18], is introduced for the solutions of dynamic responses and dynamic stability analysis of rotating shafts.The discrete singular convolution method proposed by Wei [18] has emerged as a new technique for numerical solutions of differential equations. Compared to the other conventional numerical methods, the discrete singular convolution (DSC) has the global methods’ high accuracy and the local methods’ flexibility for handling complex geometry and boundary conditions in structural analysis. And this method is also able to accurately predict very high order vibration modes of plates and beams without numerical instability. Many studies have indicated that DSC has high accuracy and reliability for free vibration and buckling analysis of beams, plates, and shells. Wei and his cooperators applied the DSC method to vibration analysis of beams [19–22] and plates [23–29]. Civalek [30–33] gave numerical solutions of free vibration problems of truncated conical shells, laminated conical and cylindrical shells, and isotropic and orthotropic rectangular plates by using DSC method. To extend the research and application scope of the DSC method, Lim et al. [34] and Hou et al. [35] developed DSC-Ritz method for vibration analysis of thick plates and shallow shells. Xiang et al. [36] explored DSC-element method for vibration analysis of rectangular Mindlin plates. And then Xiang et al. [37] also explored DSC-Ritz element method for vibration analysis of rectangular Mindlin plates with mixed edge supports. These studies indicate that the DSC method works very well for the vibration analyses of beams, plates and shells, and so forth. Recently, Wang et al. [38] applied the DSC method to solve elastic wave propagation problem. Song et al. [39] used DSC method to analyze the dynamic stability and nonstationary vibration of beams.The aim of the present paper is to propose a relatively new numerical approach using the method of DSC for the dynamic instability analysis of a rotating ship shaft under periodic axial force induced by propeller. The algorithm for solving equations of rotating shafts has been presented. The dynamic instability regions of rotating ship shaft are obtained through judging the dynamic stability of dynamic responses of transverse displacements obtained by using DSC method. Effects of number of blades, constant term in the periodic force, and damping on dynamic instability regions are discussed. And all these results by using DSC are compared with those obtained by Floquet’s method to verify the application of DSC method to dynamic stability analysis of rotating shafts.


## 2. Discrete Singular Convolution

The method of discrete singular convolution (DSC) is an effective and useful approach for the numerical solutions of differential equations, which occur commonly in many science and engineering fields. The DSC method was originally introduced by Wei [18] as a simple and highly efficient numerical technique. This paper does not tell more details about this method. For more details of the mathematical background and application of the DSC method, interested readers may refer to the works of Wei et al. [19–23, 26, 27]. The mathematical foundation of the DSC algorithm is the theory of distributions and wavelet analysis. For brevity, consider a distribution T and η(t) as an element of the space of the test function. A singular convolution can be defined by Wei [19]
(1)Ft=T*ηt=∫-∞∞Tt-xηxdx,
where T(t-x) is a singular kernel. Depending on the form of the kernel T, the singular convolution is the central issue for a wide range of science and engineering problems, for example, singular kernels of delta type [20]
(2)Tx=δnxn=0,1,2,…,
where δ is the delta distribution and the superscript n denotes the nth-order distribution. Here, kernel T(x)=δ(x) is important for interpolation and T(x)=δ(n)(x) (n=1,2,…) are essential for differentiations. However, since these kernels are singular, they cannot be directly digitized in computer. To avoid the difficulty of using singular expressions directly in computer, sequences of approximations {Tα} of the distribution T can be constructed:
(3)limα→α0Τα⟶Tx,
where α0 is a generalized limit. Obviously, in the case of T(x)=δ(x), the sequence Tα(x) is a delta sequence. Moreover, with a sufficiently smooth approximation, it is useful to consider a discrete singular convolution [19]
(4)Fαt=∑kTαt-xkfxk,
where Fα(t) is an approximation to F(t) and {xk} is an appropriate set of discrete points on which the DSC is well defined. Note that the original test function η(x) has been replaced by f(x). This new discrete expression is suitable for computer realization. Recently, the use of some new kernels and regularizer such as delta regularizer [19, 20, 23, 26] was proposed to solve applied mechanics problem. Shannon’s kernel is regularized as [23]
(5)δΔ,σx-xk=sinπ/Δx-xkπ/Δx-xkexp⁡-x-xk22σ2,σ>0,
where Δ is the grid spacing. The parameter σ determines the width of the Gaussian envelope and often varies in association with the grid spacing, σ=rΔ. It is also known that the truncation error is very small due to the use of the Gaussian regularizer; the above formulation given by (5) is practically and has an essentially compact support for numerical interpolation. In the DSC method, the function f(x) and its derivatives with respect to the x coordinate at a grid point xi are approximated by a linear sum of discrete values f(xk) in a narrow bandwidth [x-xM,x+xM]. This can be expressed as [28]
(6)dnfxdxnx=xi=fnx≈∑k=-MMδΔ,σnxi-xkfxkn=0,1,2,…,
where superscript (n) denotes the nth-order derivative with respect to x. When the regularized Shannon’s delta (RSD) kernel is used, the detailed expressions for δΔ,σ(n) can be easily obtained. Readers can refer to some published references [20, 23].


## 3. Theory and Algorithm

### 3.1. Equation of Motion

Figure1 illustrates a uniform simply supported Euler shaft rotating about its longitudinal axis with a constant spinning speed Ω, subjected to a periodic axial force P(t)=PS+PDcos⁡(ZpΩt), and (XYZ) is the fixed coordinate system while (xyz) is the coordinate system attached to the shaft with the x-axis aligned with the x-axis. Using Euler-Bernoulli beam theory, the effect of rotatory inertia is neglected. The equations [9, 12] of motion can be written as
(7)m∂2uy∂t2-2mΩ∂uz∂t-mΩ2uy+C∂uy∂t+EI∂4uy∂x4+Pt∂2uy∂x2=0,m∂2uz∂t2+2mΩ∂uy∂t-mΩ2uz+C∂uz∂t+EI∂4uz∂x4+Pt∂2uz∂x2=0,
where uy, uz are the transverse displacements in the y- and z-direction, m is the mass per unit length of shaft, and C, E, I, and t are the viscous damping coefficient, Young’s modulus, axial moment of inertia, and time, respectively. PS is the constant term, PD is the amplitude, and Zp is the number of the blades on the screw propeller.Figure 1
The rotating ship shaft and coordinate systems.For simplicity, introduce the following dimensionless quantities:(8)X=xlUY=uylUZ=uzl,α=PSPcrβ=PDPcrΩ¯=Ωω,ζ=C2mωτ=ωtω*=ω1-α,Ω*¯=Ωω*=Ω¯1-αU=UY2+UZ2,
where Pcr=π2EI/l2 is the buckling force of a nonspinning simply supported shaft, ω=π2EI/ml4 is the natural bending frequency for the first mode, ω*=ω1-α is the natural bending frequency with constant axial compressive force, l is the length of the shaft, ζ is damping ratio, U is composition of UY and UZ, and Ω¯ is dimensionless spinning speed.Then, (7) can be written as
(9)∂2UY∂τ2-2Ω¯∂UZ∂τ-Ω¯2UY+2ζ∂UY∂τ+1π4∂4UY∂X4+α+βcos⁡ZpΩ¯τπ2∂2UY∂X2=0,∂2UZ∂τ2+2Ω¯∂UY∂τ-Ω¯2UZ+2ζ∂UZ∂τ+1π4∂4UZ∂X4+α+βcos⁡ZpΩ¯τπ2∂2UZ∂X2=0.And (9) can be changed into
(10)∂2UY∂τ2=2Ω¯∂UZ∂τ+Ω¯2UY-2ζ∂UY∂τ-1π4∂4UY∂X4-α+βcos⁡⁡ZpΩ¯τπ2∂2UY∂X2,∂2UZ∂τ2=-2Ω¯∂UY∂τ+Ω¯2UZ-2ζ∂UZ∂τ-1π4∂4UZ∂X4-α+βcos⁡⁡ZpΩ¯τπ2∂2UZ∂X2.The boundary conditions are(11)U0,τ=U1,τ=0,∂2U0,τ∂X2=∂2U1,τ∂X2=0.


### 3.2. Discretization and Algorithm

The DSC discrete scheme (6) is utilized for the spatial discretization and the fourth-order Runge–Kutta (RK4) scheme is used for the time discretization [39, 40]. The computational domain of coordinate X is [0,1], and the coordinate X is equally spaced; the grid sizes are denoted by ΔX=(1-0)/N(N is the total number of partition grids on the computational domain [0,1]); the grid points are denoted by Xj=(j-1)ΔX  (j=1,2,…,N+1). So Xj-Xj+k=-kΔX. The approximate values of UY and UZ at the grid point Xj are expressed as UY,j and UZ,j Then, (10) can be expressed as
(12)∂∂τUY,jUZ,j=VY,jVZ,j,∂∂τVY,jVZ,j=-2ζ2Ω¯-2Ω¯-2ζVY,jVZ,j+Ω¯200Ω¯2UY,jUZ,j-1π4∂4∂X4+α+βcos⁡ZpΩ¯τπ2∂2∂X2UY,jUZ,jj=1,2,…,N+1.Let(13)Y1,j=Y1,1,Y1,2,…,Y1,2N+2=UY,1,UY,2,…,UY,N+1,VY,1,VY,2,…,VY,N+1,Y2,j=Y2,1,Y2,2,…,Y2,2N+2=UZ,1,UZ,2,…,UZ,N+1,VZ,1,VZ,2,…,VZ,N+1j=1,2,…,2N+2.Whenj=1,2,…,N+1(14)Fi,j=Yi,j+N+1i=1,2
and when j=N+2,N+3,…,2N+2(15)F1,jF2,j=-2ζ2Ω¯-2Ω¯-2ζY1,jY2,j+Ω¯200Ω¯2Y1,j-N-1Y2,j-N-1-1π4∂4∂X4+α+βcos⁡ZpΩ¯τπ2∂2∂X2Y1,j-N-1Y2,j-N-1.Then, we have a unified semidiscretized equation:(16)∂Yi,j∂τ=Fi,ji=1,2,j=1,2,…,2N+2.The temporal discretization expressions of (16) by using fourth-order Runge–Kutta method are given as
(17)Yi,jn+1=Yi,jn+ΔτYi,j+N+1n+Δτ26Li,j+N+1,1+Li,j+N+1,2+Li,j+N+1,3i=1,2,j=1,2,…,N+1,Yi,jn+1=Yi,jn+Δτ6Li,j,1+2Li,j,2+2Li,j,3+Li,j,4i=1,2,j=N+2,N+3,…,2N+2,
where Li,j,1=Fi,j,1n, Li,j,2=Fi,j,2n, Li,j,3=Fi,j,3n, and Li,j,4=Fi,j,4n (i=1,2,j=N+2,N+3,…,2N+2) and here superscript n denotes time level; Δτ is the time step, so τ=nΔτ. Using the DSC discrete scheme (6), the Fi,j,1n, Fi,j,2n, Fi,j,3n, and Fi,j,4n discretization expressions are
(18)L1,j,1L2,j,1=F1,j,1nF2,j,1n=-2ζ2Ω¯-2Ω¯-2ζY1,jnY2,jn+Ω¯200Ω¯2Y1,j-N-1nY2,j-N-1n-α+βcos⁡ZpΩ¯τλl4π2η2∑k=-W+WδΔ,σ2-kΔX1λl4∑k=-W+WδΔ,σ4-kΔX+α+βcos⁡ZpΩ¯τλl4π2η2∑k=-W+WδΔ,σ2-kΔX1λl4∑k=-W+WδΔ,σ4-kΔX·Y1,j+k-N-1nY2,j+k-N-1n,L1,j,2L2,j,2=F1,j,2nF2,j,2n=-2ζ2Ω¯-2Ω¯-2ζY1,jn+Δτ2L1,j,1Y2,jn+Δτ2L2,j,1+Ω¯200Ω¯2Y1,j-N-1n+Δτ2Y1,jnY2,j-N-1n+Δτ2Y2,jn-1λl4∑k=-W+WδΔ,σ4-kΔX+α+βcos⁡ZpΩ¯τλl4π2η2·∑k=-W+WδΔ,σ2-kΔXY1,j+k-N-1n+Δτ2Y1,j+knY2,j+k-N-1n+Δτ2Y2,j+knL1,j,3L2,j,3=F1,j,3nF2,j,3n=-2ζ2Ω¯-2Ω¯-2ζY1,jn+Δτ2L1,j,2Y2,jn+Δτ2L2,j,2+Ω¯200Ω¯2Y1,j-N-1n+Δτ2Y1,jn+Δτ24L1,j,1Y2,j-N-1n+Δτ2Y2,jn+Δτ24L2,j,1-1λl4∑k=-W+WδΔ,σ4-kΔX+α+βcos⁡ZpΩ¯τλl4π2η2·∑k=-W+WδΔ,σ2-kΔX·Y1,j+k-N-1n+Δτ2Y1,j+kn+Δτ24L1,j+k,1Y2,j+k-N-1n+Δτ2Y2,j+kn+Δτ24L2,j+k,1L1,j,4L2,j,4=F1,j,4nF2,j,4n=-2ζ2Ω¯-2Ω¯-2ζY1,jn+ΔτL1,j,3Y2,jn+ΔτL2,j,3+Ω¯200Ω¯2Y1,j-N-1n+ΔτY1,jn+Δτ22L1,j,2Y2,j-N-1n+ΔτY2,jn+Δτ22L2,j,2-1λl4∑k=-W+WδΔ,σ4-kΔX+α+βcos⁡ZpΩ¯τλl4π2η2·∑k=-W+WδΔ,σ2-kΔX·Y1,j+k-N-1n+ΔτY1,j+kn+Δτ22L1,j+k,2Y2,j+k-N-1n+ΔτY2,j+kn+Δτ22L2,j+k,2j=N+2,N+3,…,2N+2,
where [-W,+W] is the computational bandwidth. Kernels δΔ,σ(2) and δΔ,σ(4) can be easily obtained [20, 23]. All of these coefficients are only dependent on grid size. When the grid point distribution is given, the coefficients can be computed once and stored for use during the computation.The overall calculation procedure can be summarized as follows [39, 40]:(a)
the initial values forYi,j0 or the values of Yi,jn (i=1,2,j=1,2,…,2N+2) at previous time level n can be obtained;(b)
the values ofFi,j,1n, Fi,j,2n, Fi,j,3n, and Fi,j,4n (i=1,2,j=N+2,N+3,…,2N+2) at the time level n can be obtained;(c)
substituting the values in (18) into (17), the values of Yi,jn+1 (i=1,2,j=1,2,…,2N+2) at new time level n+1 can be calculated;(d)
the computational time is advanced (i.e.,τ=τ+Δτ, n=n+1), and the whole procedure above is repeated, until calculation precision is reached.


## 3.1. Equation of Motion

Figure1 illustrates a uniform simply supported Euler shaft rotating about its longitudinal axis with a constant spinning speed Ω, subjected to a periodic axial force P(t)=PS+PDcos⁡(ZpΩt), and (XYZ) is the fixed coordinate system while (xyz) is the coordinate system attached to the shaft with the x-axis aligned with the x-axis. Using Euler-Bernoulli beam theory, the effect of rotatory inertia is neglected. The equations [9, 12] of motion can be written as
(7)m∂2uy∂t2-2mΩ∂uz∂t-mΩ2uy+C∂uy∂t+EI∂4uy∂x4+Pt∂2uy∂x2=0,m∂2uz∂t2+2mΩ∂uy∂t-mΩ2uz+C∂uz∂t+EI∂4uz∂x4+Pt∂2uz∂x2=0,
where uy, uz are the transverse displacements in the y- and z-direction, m is the mass per unit length of shaft, and C, E, I, and t are the viscous damping coefficient, Young’s modulus, axial moment of inertia, and time, respectively. PS is the constant term, PD is the amplitude, and Zp is the number of the blades on the screw propeller.Figure 1
The rotating ship shaft and coordinate systems.For simplicity, introduce the following dimensionless quantities:(8)X=xlUY=uylUZ=uzl,α=PSPcrβ=PDPcrΩ¯=Ωω,ζ=C2mωτ=ωtω*=ω1-α,Ω*¯=Ωω*=Ω¯1-αU=UY2+UZ2,
where Pcr=π2EI/l2 is the buckling force of a nonspinning simply supported shaft, ω=π2EI/ml4 is the natural bending frequency for the first mode, ω*=ω1-α is the natural bending frequency with constant axial compressive force, l is the length of the shaft, ζ is damping ratio, U is composition of UY and UZ, and Ω¯ is dimensionless spinning speed.Then, (7) can be written as
(9)∂2UY∂τ2-2Ω¯∂UZ∂τ-Ω¯2UY+2ζ∂UY∂τ+1π4∂4UY∂X4+α+βcos⁡ZpΩ¯τπ2∂2UY∂X2=0,∂2UZ∂τ2+2Ω¯∂UY∂τ-Ω¯2UZ+2ζ∂UZ∂τ+1π4∂4UZ∂X4+α+βcos⁡ZpΩ¯τπ2∂2UZ∂X2=0.And (9) can be changed into
(10)∂2UY∂τ2=2Ω¯∂UZ∂τ+Ω¯2UY-2ζ∂UY∂τ-1π4∂4UY∂X4-α+βcos⁡⁡ZpΩ¯τπ2∂2UY∂X2,∂2UZ∂τ2=-2Ω¯∂UY∂τ+Ω¯2UZ-2ζ∂UZ∂τ-1π4∂4UZ∂X4-α+βcos⁡⁡ZpΩ¯τπ2∂2UZ∂X2.The boundary conditions are(11)U0,τ=U1,τ=0,∂2U0,τ∂X2=∂2U1,τ∂X2=0.


## 3.2. Discretization and Algorithm

The DSC discrete scheme (6) is utilized for the spatial discretization and the fourth-order Runge–Kutta (RK4) scheme is used for the time discretization [39, 40]. The computational domain of coordinate X is [0,1], and the coordinate X is equally spaced; the grid sizes are denoted by ΔX=(1-0)/N(N is the total number of partition grids on the computational domain [0,1]); the grid points are denoted by Xj=(j-1)ΔX  (j=1,2,…,N+1). So Xj-Xj+k=-kΔX. The approximate values of UY and UZ at the grid point Xj are expressed as UY,j and UZ,j Then, (10) can be expressed as
(12)∂∂τUY,jUZ,j=VY,jVZ,j,∂∂τVY,jVZ,j=-2ζ2Ω¯-2Ω¯-2ζVY,jVZ,j+Ω¯200Ω¯2UY,jUZ,j-1π4∂4∂X4+α+βcos⁡ZpΩ¯τπ2∂2∂X2UY,jUZ,jj=1,2,…,N+1.Let(13)Y1,j=Y1,1,Y1,2,…,Y1,2N+2=UY,1,UY,2,…,UY,N+1,VY,1,VY,2,…,VY,N+1,Y2,j=Y2,1,Y2,2,…,Y2,2N+2=UZ,1,UZ,2,…,UZ,N+1,VZ,1,VZ,2,…,VZ,N+1j=1,2,…,2N+2.Whenj=1,2,…,N+1(14)Fi,j=Yi,j+N+1i=1,2
and when j=N+2,N+3,…,2N+2(15)F1,jF2,j=-2ζ2Ω¯-2Ω¯-2ζY1,jY2,j+Ω¯200Ω¯2Y1,j-N-1Y2,j-N-1-1π4∂4∂X4+α+βcos⁡ZpΩ¯τπ2∂2∂X2Y1,j-N-1Y2,j-N-1.Then, we have a unified semidiscretized equation:(16)∂Yi,j∂τ=Fi,ji=1,2,j=1,2,…,2N+2.The temporal discretization expressions of (16) by using fourth-order Runge–Kutta method are given as
(17)Yi,jn+1=Yi,jn+ΔτYi,j+N+1n+Δτ26Li,j+N+1,1+Li,j+N+1,2+Li,j+N+1,3i=1,2,j=1,2,…,N+1,Yi,jn+1=Yi,jn+Δτ6Li,j,1+2Li,j,2+2Li,j,3+Li,j,4i=1,2,j=N+2,N+3,…,2N+2,
where Li,j,1=Fi,j,1n, Li,j,2=Fi,j,2n, Li,j,3=Fi,j,3n, and Li,j,4=Fi,j,4n (i=1,2,j=N+2,N+3,…,2N+2) and here superscript n denotes time level; Δτ is the time step, so τ=nΔτ. Using the DSC discrete scheme (6), the Fi,j,1n, Fi,j,2n, Fi,j,3n, and Fi,j,4n discretization expressions are
(18)L1,j,1L2,j,1=F1,j,1nF2,j,1n=-2ζ2Ω¯-2Ω¯-2ζY1,jnY2,jn+Ω¯200Ω¯2Y1,j-N-1nY2,j-N-1n-α+βcos⁡ZpΩ¯τλl4π2η2∑k=-W+WδΔ,σ2-kΔX1λl4∑k=-W+WδΔ,σ4-kΔX+α+βcos⁡ZpΩ¯τλl4π2η2∑k=-W+WδΔ,σ2-kΔX1λl4∑k=-W+WδΔ,σ4-kΔX·Y1,j+k-N-1nY2,j+k-N-1n,L1,j,2L2,j,2=F1,j,2nF2,j,2n=-2ζ2Ω¯-2Ω¯-2ζY1,jn+Δτ2L1,j,1Y2,jn+Δτ2L2,j,1+Ω¯200Ω¯2Y1,j-N-1n+Δτ2Y1,jnY2,j-N-1n+Δτ2Y2,jn-1λl4∑k=-W+WδΔ,σ4-kΔX+α+βcos⁡ZpΩ¯τλl4π2η2·∑k=-W+WδΔ,σ2-kΔXY1,j+k-N-1n+Δτ2Y1,j+knY2,j+k-N-1n+Δτ2Y2,j+knL1,j,3L2,j,3=F1,j,3nF2,j,3n=-2ζ2Ω¯-2Ω¯-2ζY1,jn+Δτ2L1,j,2Y2,jn+Δτ2L2,j,2+Ω¯200Ω¯2Y1,j-N-1n+Δτ2Y1,jn+Δτ24L1,j,1Y2,j-N-1n+Δτ2Y2,jn+Δτ24L2,j,1-1λl4∑k=-W+WδΔ,σ4-kΔX+α+βcos⁡ZpΩ¯τλl4π2η2·∑k=-W+WδΔ,σ2-kΔX·Y1,j+k-N-1n+Δτ2Y1,j+kn+Δτ24L1,j+k,1Y2,j+k-N-1n+Δτ2Y2,j+kn+Δτ24L2,j+k,1L1,j,4L2,j,4=F1,j,4nF2,j,4n=-2ζ2Ω¯-2Ω¯-2ζY1,jn+ΔτL1,j,3Y2,jn+ΔτL2,j,3+Ω¯200Ω¯2Y1,j-N-1n+ΔτY1,jn+Δτ22L1,j,2Y2,j-N-1n+ΔτY2,jn+Δτ22L2,j,2-1λl4∑k=-W+WδΔ,σ4-kΔX+α+βcos⁡ZpΩ¯τλl4π2η2·∑k=-W+WδΔ,σ2-kΔX·Y1,j+k-N-1n+ΔτY1,j+kn+Δτ22L1,j+k,2Y2,j+k-N-1n+ΔτY2,j+kn+Δτ22L2,j+k,2j=N+2,N+3,…,2N+2,
where [-W,+W] is the computational bandwidth. Kernels δΔ,σ(2) and δΔ,σ(4) can be easily obtained [20, 23]. All of these coefficients are only dependent on grid size. When the grid point distribution is given, the coefficients can be computed once and stored for use during the computation.The overall calculation procedure can be summarized as follows [39, 40]:(a)
the initial values forYi,j0 or the values of Yi,jn (i=1,2,j=1,2,…,2N+2) at previous time level n can be obtained;(b)
the values ofFi,j,1n, Fi,j,2n, Fi,j,3n, and Fi,j,4n (i=1,2,j=N+2,N+3,…,2N+2) at the time level n can be obtained;(c)
substituting the values in (18) into (17), the values of Yi,jn+1 (i=1,2,j=1,2,…,2N+2) at new time level n+1 can be calculated;(d)
the computational time is advanced (i.e.,τ=τ+Δτ, n=n+1), and the whole procedure above is repeated, until calculation precision is reached.


## 4. Numerical Results and Discussions

### 4.1. Dynamic Instability Regions for Shafts without Damping

In this section, first, dynamic instability regions with different number of bladesZp are carried out. And then effect of constant term in the periodic axial force on dynamic instability region is discussed. In order to discuss the effects of number of blades Zp on dynamic instability regions, the value of α is chosen as 0.1. Then the frequency is ω*=ω1-α=0.9ω. The initial conditions and parameters in DSC algorithm are set as follows:U
Y
(
i
Δ
X
,
0
)
=
0.001
sin
(
i
Δ
X
π
), UZ(iΔX,0)=0, ∂UY(iΔX,0)/∂τ=0, ∂UZ(iΔX,0)/∂τ=0 (i=0,1,2,…,N), N=16, W=15, r=2.5, Δτ = 3.0 × 10−6. For the choice of N, W, and r, readers may refer to Wei et al. [20, 26].Figure2 shows dynamic responses of transverse dimensionless displacements at midpoint for the rotating shaft with Zp = 5 and ζ = 0. It is seen that the dynamic responses for ZpΩ/2ω* = 1.0 β = 0.05, ZpΩ/2ω* = 0.9 β = 0.36, and ZpΩ/2ω* = 1.1 β = 0.38 in Figures 2(a), 2(c), and 2(e) are dynamically unstable, while the dynamic responses for ZpΩ/2ω* = 0.9 β = 0.34 and ZpΩ/2ω* = 1.1 β = 0.36 in Figures 2(b) and 2(d) are dynamically stable. The whole procedure above is repeated, the dynamic instability region for the first mode shown as Figure 3 is obtained. In Figure 3, the dynamic instability region obtained by DSC algorithm is coincident with the instability region obtained by Floquet’s method [15] very well. It is also found that the dynamic instability region is “V” shape, and when ZpΩ/2ω* = 1.0, that is to say, the excitation frequency is twice the natural frequency with constant axial compressive force, the dynamic instability is the most serious, which is similar to the dynamic instability region when a nonspinning shaft is subjected to periodic force [1].Dynamic responses of transverse displacements withZp=5. (a) ZpΩ/2ω*=1.0  β = 0.05, (b) ZpΩ/2ω*=0.9  β = 0.34, (c) ZpΩ/2ω*=0.9  β = 0.36, (d) ZpΩ/2ω*=1.1  β = 0.36, and (e) ZpΩ/2ω*=1.1  β = 0.38.
(a)
(b)
(c)
(d)
(e)Figure 3
Dynamic instability region of a rotating ship shaft withZp = 5. “—,” instability region boundary by Floquet’s method; “●” and “○,” stable and unstable results by DSC method.The dynamic instability region withZp = 3 shown as Figure 4 is also obtained by using DSC algorithm. It is found that the dynamic instability region is “V” shape, and when ZpΩ/2ω* = 1.0, the dynamic instability is the most serious. And also the instability region obtained by DSC algorithm agrees well with that obtained by Floquet’s method. The comparisons of dynamic instability regions with different number of blades Zp are shown in Figure 5. It is seen that the dynamic instability regions obtained by DSC algorithm and Floquet’s method are “V” shapes and do not change with the number of blades Zp increasing. That is to say, the increase of number of blades Zp has little influence on dynamic instability region of rotating ship shafts without damping.Figure 4
Dynamic instability region of a rotating ship shaft withZp=3. Key as for Figure 3.Figure 5
Dynamic instability regions with different number of bladesZp. “—,” “- -,” instability region boundary with Zp=5, Zp=3 by Floquet’s method; “●” “○,” “▼” “▽,” stable and unstable results with Zp=5, Zp=3 by DSC method.In order to discuss the effect of constant term in the periodic force on dynamic instability region, the value ofα is increased to 0.3. And the corresponding frequency ω*=ω1-α = 0.7ω = 0.8367ω. When α = 0.1, the corresponding frequency is ω* = 0.9ω = 0.9487ω. The comparisons of dynamic instability regions with different α are shown in Figure 6. It is found that the dynamic instability regions with α = 0.1 and 0.3 are “V” shape. And when the excitation frequency ZpΩ is near twice the natural frequency with constant axial compressive force, the dynamic instability is the most serious. The increase in the value of α causes the dynamic instability regions shifting to lower exciting frequencies, which makes the shaft more sensitive to periodic force.Figure 6
Dynamic instability regions with differentα. “—,” “- -,” instability region boundary with α = 0.1, 0.3 by Floquet’s method; “●” “○,” “▼” “▽,” stable and unstable results with α = 0.1, 0.3 by DSC method.


### 4.2. Dynamic Instability Regions for Shafts with Damping

In this section, dynamic instability regions with damping effects are obtained, and effects of different number of bladesZp are also discussed. The value of ζ is chosen as 0.05. And the value of α is chosen as 0.1. The initial conditions and other parameters in DSC algorithm are the same as the previous computation.Dynamic responses of transverse dimensionless displacements at midpoint for the rotating shaft withZp = 5 and ζ = 0.05 are shown in Figure 7. It is seen that the dynamic responses for ZpΩ/2ω* = 1.0 β = 0.18 and ZpΩ/2ω* = 1.1 β = 0.40 in Figures 7(b) and 7(d) are dynamically unstable, while the dynamic responses for ZpΩ/2ω* = 1.0 β = 0.16 and ZpΩ/2ω* = 1.1 β = 0.38 in Figures 7(a) and 7(c) are dynamically stable. The whole procedure above is repeated; the dynamic instability region shown in Figure 8 is obtained. In Figure 8, the dynamic instability region obtained by DSC algorithm is coincident with that obtained by Floquet’s method very well. And it is also found that damping improves dynamic stability of the region with ZpΩ/2ω* near 1.0 obviously.Dynamic responses of transverse displacements withZp=5 and ζ=0.05. (a) ZpΩ/2ω*=1.0β=0.16, (b) ZpΩ/2ω*=1.0β=0.18, (c) ZpΩ/2ω*=1.1β=0.38, and (d) ZpΩ/2ω*=1.1β=0.40.
(a)
(b)
(c)
(d)Figure 8
Dynamic instability region of a rotating ship shaft withZp = 5 and ζ = 0.05. “—,” “- -,” instability region boundary with Zp = 5 ζ = 0.05 and Zp = 5 ζ = 0 by Floquet’s method, “●,” “○,” stable and unstable data with Zp = 5 ζ = 0.05 by DSC method.The dynamic instability region withZp = 3 and ζ = 0.05 shown in Figure 9 is also obtained by using DSC algorithm. It is found that the dynamic stability of the region with ZpΩ/2ω* near 1.0 is improved obviously due to the existent of damping, while dynamic stability of other regions decreases. The comparisons of dynamic instability regions with different number of blades Zp and ζ = 0.05 are shown in Figure 10. It is also seen that the dynamic instability regions obtained by DSC algorithm agree well with those obtained by Floquet’s method, respectively. And the dynamic instability regions expand with decreasing of the number of blades Zp; that is to say, the more blades on the propeller, the better of dynamic stability, which is different from that obtained from Figure 5. Moreover, the dynamic instability of region with ZpΩ/2ω* near 1.0 is also the most serious, although the dynamic stability of this region is better than that shown in Figure 5. Therefore, in order to prevent dynamic instability occurring, the spinning speed Ω should be far away from 2ω*/Zp.Figure 9
Dynamic instability region of a rotating ship shaft withZp = 3 and ζ = 0.05. “—,” “- -,” instability region boundary with Zp = 3, ζ = 0.05 and Zp = 3, ζ = 0 by Floquet’s method; “●,” “○,” stable and unstable data with Zp = 3, ζ = 0.05 by DSC method.Figure 10
Dynamic instability regions with different numbers of bladesZp and ζ = 0.05. “—,” “- -,” instability region boundary with Zp = 5, Zp = 3 by Floquet’s method, “●” “○,” “▼” “▽,” stable and unstable data with Zp = 5, Zp = 3 by DSC method.


## 4.1. Dynamic Instability Regions for Shafts without Damping

In this section, first, dynamic instability regions with different number of bladesZp are carried out. And then effect of constant term in the periodic axial force on dynamic instability region is discussed. In order to discuss the effects of number of blades Zp on dynamic instability regions, the value of α is chosen as 0.1. Then the frequency is ω*=ω1-α=0.9ω. The initial conditions and parameters in DSC algorithm are set as follows:U
Y
(
i
Δ
X
,
0
)
=
0.001
sin
(
i
Δ
X
π
), UZ(iΔX,0)=0, ∂UY(iΔX,0)/∂τ=0, ∂UZ(iΔX,0)/∂τ=0 (i=0,1,2,…,N), N=16, W=15, r=2.5, Δτ = 3.0 × 10−6. For the choice of N, W, and r, readers may refer to Wei et al. [20, 26].Figure2 shows dynamic responses of transverse dimensionless displacements at midpoint for the rotating shaft with Zp = 5 and ζ = 0. It is seen that the dynamic responses for ZpΩ/2ω* = 1.0 β = 0.05, ZpΩ/2ω* = 0.9 β = 0.36, and ZpΩ/2ω* = 1.1 β = 0.38 in Figures 2(a), 2(c), and 2(e) are dynamically unstable, while the dynamic responses for ZpΩ/2ω* = 0.9 β = 0.34 and ZpΩ/2ω* = 1.1 β = 0.36 in Figures 2(b) and 2(d) are dynamically stable. The whole procedure above is repeated, the dynamic instability region for the first mode shown as Figure 3 is obtained. In Figure 3, the dynamic instability region obtained by DSC algorithm is coincident with the instability region obtained by Floquet’s method [15] very well. It is also found that the dynamic instability region is “V” shape, and when ZpΩ/2ω* = 1.0, that is to say, the excitation frequency is twice the natural frequency with constant axial compressive force, the dynamic instability is the most serious, which is similar to the dynamic instability region when a nonspinning shaft is subjected to periodic force [1].Dynamic responses of transverse displacements withZp=5. (a) ZpΩ/2ω*=1.0  β = 0.05, (b) ZpΩ/2ω*=0.9  β = 0.34, (c) ZpΩ/2ω*=0.9  β = 0.36, (d) ZpΩ/2ω*=1.1  β = 0.36, and (e) ZpΩ/2ω*=1.1  β = 0.38.
(a)
(b)
(c)
(d)
(e)Figure 3
Dynamic instability region of a rotating ship shaft withZp = 5. “—,” instability region boundary by Floquet’s method; “●” and “○,” stable and unstable results by DSC method.The dynamic instability region withZp = 3 shown as Figure 4 is also obtained by using DSC algorithm. It is found that the dynamic instability region is “V” shape, and when ZpΩ/2ω* = 1.0, the dynamic instability is the most serious. And also the instability region obtained by DSC algorithm agrees well with that obtained by Floquet’s method. The comparisons of dynamic instability regions with different number of blades Zp are shown in Figure 5. It is seen that the dynamic instability regions obtained by DSC algorithm and Floquet’s method are “V” shapes and do not change with the number of blades Zp increasing. That is to say, the increase of number of blades Zp has little influence on dynamic instability region of rotating ship shafts without damping.Figure 4
Dynamic instability region of a rotating ship shaft withZp=3. Key as for Figure 3.Figure 5
Dynamic instability regions with different number of bladesZp. “—,” “- -,” instability region boundary with Zp=5, Zp=3 by Floquet’s method; “●” “○,” “▼” “▽,” stable and unstable results with Zp=5, Zp=3 by DSC method.In order to discuss the effect of constant term in the periodic force on dynamic instability region, the value ofα is increased to 0.3. And the corresponding frequency ω*=ω1-α = 0.7ω = 0.8367ω. When α = 0.1, the corresponding frequency is ω* = 0.9ω = 0.9487ω. The comparisons of dynamic instability regions with different α are shown in Figure 6. It is found that the dynamic instability regions with α = 0.1 and 0.3 are “V” shape. And when the excitation frequency ZpΩ is near twice the natural frequency with constant axial compressive force, the dynamic instability is the most serious. The increase in the value of α causes the dynamic instability regions shifting to lower exciting frequencies, which makes the shaft more sensitive to periodic force.Figure 6
Dynamic instability regions with differentα. “—,” “- -,” instability region boundary with α = 0.1, 0.3 by Floquet’s method; “●” “○,” “▼” “▽,” stable and unstable results with α = 0.1, 0.3 by DSC method.


## 4.2. Dynamic Instability Regions for Shafts with Damping

In this section, dynamic instability regions with damping effects are obtained, and effects of different number of bladesZp are also discussed. The value of ζ is chosen as 0.05. And the value of α is chosen as 0.1. The initial conditions and other parameters in DSC algorithm are the same as the previous computation.Dynamic responses of transverse dimensionless displacements at midpoint for the rotating shaft withZp = 5 and ζ = 0.05 are shown in Figure 7. It is seen that the dynamic responses for ZpΩ/2ω* = 1.0 β = 0.18 and ZpΩ/2ω* = 1.1 β = 0.40 in Figures 7(b) and 7(d) are dynamically unstable, while the dynamic responses for ZpΩ/2ω* = 1.0 β = 0.16 and ZpΩ/2ω* = 1.1 β = 0.38 in Figures 7(a) and 7(c) are dynamically stable. The whole procedure above is repeated; the dynamic instability region shown in Figure 8 is obtained. In Figure 8, the dynamic instability region obtained by DSC algorithm is coincident with that obtained by Floquet’s method very well. And it is also found that damping improves dynamic stability of the region with ZpΩ/2ω* near 1.0 obviously.Dynamic responses of transverse displacements withZp=5 and ζ=0.05. (a) ZpΩ/2ω*=1.0β=0.16, (b) ZpΩ/2ω*=1.0β=0.18, (c) ZpΩ/2ω*=1.1β=0.38, and (d) ZpΩ/2ω*=1.1β=0.40.
(a)
(b)
(c)
(d)Figure 8
Dynamic instability region of a rotating ship shaft withZp = 5 and ζ = 0.05. “—,” “- -,” instability region boundary with Zp = 5 ζ = 0.05 and Zp = 5 ζ = 0 by Floquet’s method, “●,” “○,” stable and unstable data with Zp = 5 ζ = 0.05 by DSC method.The dynamic instability region withZp = 3 and ζ = 0.05 shown in Figure 9 is also obtained by using DSC algorithm. It is found that the dynamic stability of the region with ZpΩ/2ω* near 1.0 is improved obviously due to the existent of damping, while dynamic stability of other regions decreases. The comparisons of dynamic instability regions with different number of blades Zp and ζ = 0.05 are shown in Figure 10. It is also seen that the dynamic instability regions obtained by DSC algorithm agree well with those obtained by Floquet’s method, respectively. And the dynamic instability regions expand with decreasing of the number of blades Zp; that is to say, the more blades on the propeller, the better of dynamic stability, which is different from that obtained from Figure 5. Moreover, the dynamic instability of region with ZpΩ/2ω* near 1.0 is also the most serious, although the dynamic stability of this region is better than that shown in Figure 5. Therefore, in order to prevent dynamic instability occurring, the spinning speed Ω should be far away from 2ω*/Zp.Figure 9
Dynamic instability region of a rotating ship shaft withZp = 3 and ζ = 0.05. “—,” “- -,” instability region boundary with Zp = 3, ζ = 0.05 and Zp = 3, ζ = 0 by Floquet’s method; “●,” “○,” stable and unstable data with Zp = 3, ζ = 0.05 by DSC method.Figure 10
Dynamic instability regions with different numbers of bladesZp and ζ = 0.05. “—,” “- -,” instability region boundary with Zp = 5, Zp = 3 by Floquet’s method, “●” “○,” “▼” “▽,” stable and unstable data with Zp = 5, Zp = 3 by DSC method.


## 5. Conclusions

In the present study, the method of discrete singular convolution (DSC) has been applied to solve vibration governing equations of motion for rotating ship shafts under periodic axial forces with type ofP(t)=PS+PDcos⁡(ZpΩt). The excitation frequency of the periodic axial force is related to the spinning speed and the number of blades. The dynamic instability regions of rotating ship shafts are obtained through judging the dynamic stability of dynamic responses of transverse displacements. Effects of number of blades, constant term in the periodic force, and damping on dynamic instability regions also have been discussed. The following can be concluded.(1)
The dynamic instability regions for different cases obtained by using DSC method are coincident with the instability regions obtained by Floquet’s method very well, respectively, which verifies the application of DSC method to this class of problems.(2)
For the rotating ship shaft without damping, the dynamic instability regions are “V” shape and do not change with the number of blades increasing. The increase in the value of the constant term in the periodic force causes the dynamic instability regions shifting to lower exciting frequencies, which makes the shaft more sensitive to periodic force.(3)
For the rotating ship shaft with damping, the dynamic stability of regions where the excitation frequency is in the neighborhood of twice the natural frequency with constant axial compressive force is improved obviously, due to damping existing, while the dynamic stability of other regions decrease. In addition, the dynamic instability regions expand with the number of blades decreasing; that is to say, the increase of number of blades improves the dynamic stability, which is very different from the result for the rotating ship shaft without damping. Therefore, it is suggested that the screw propeller should have more blades for rotating ship shafts.(4)
When the excitation frequency is twice the natural frequency with constant axial compressive force, the dynamic instability is the most serious for the system with no damping and with damping. Therefore, the spinning speedΩ should be far away from 2ω*/Zp, to avoid serious dynamic instability happening.



---

*Source: 482607-2015-03-02.xml*