1. society
2. crime
3. property crime
4. piracy

Simplified modelling of rotor cracks
*
†
J.E.T. Penny and M.I. Friswell
*
School of Engineering and Applied Science, Aston University, Birmingham B4 7ET, UK
†
Department of Mechanical Engineering, University of Wales Swansea, Swansea SA2 8PP, UK
e-mail: [email protected], [email protected]
Abstract
There are a number of approaches reported in the literature for modelling cracks in shafts. Although two and
three dimensional finite element models may be used, the most popular approaches are based on beam
models. This paper considers a number of different approaches to modelling cracks in shafts of rotating
machines. Of particular concern is the modelling of breathing cracks, which open and close due to the selfweight of the rotor, producing a parametric excitation. In the simplest model the shaft stiffness is bilinear,
depending on whether the crack is open or closed. More complicated functions that relate the shaft stiffness
to angular position have been proposed, and two possibilities will be considered. This paper will demonstrate
the influence of the crack model on the response of a Jeffcott rotor, using the harmonic balance approach.
to gravity. In this case the crack will open and close
(or breathe) due to the rotation of the rotor. This
problem was initially studied by Gasch [3,4], Mayes
and Davies [5,6], Grabowski [7] and others. Various
crack opening and closing models have been
proposed. In this paper the vibration response of a
simple Jeffcott rotor with three different crack
models is compared.
1. Introduction
The notion that changes in a rotor’s dynamic
behaviour could be used for on-line crack detection
and monitoring was first proposed in the 1970s and
much of the early work was presented in technical
reports of the General Electric Co, USA. To
implement this method of crack detection either a
model based method or a “pattern recognition”
based method is used. The pattern recognition
approach looks for changes in the vibration data
using some suitable analysis method, for example
Prabhakar et al. [1] used wavelet transforms to
detect and monitor cracks. The alternative model
based approach allows vibration data from the rotor
to be compared with simulated data from the rotor
and crack model. The crack location and severity
can then be determined by using an appropriate
mathematical technique to adjust these crack model
parameters to minimise the error between the actual
and simulated data [2].
2. Models of opening and closing
cracks
The simplest model of an opening and closing crack
is the hinge model [3,4]. In this model, the crack is
assumed to change from its closed to open state (and
visa-versa) abruptly as the shaft rotates. The
stiffness matrix in rotating coordinates (ξ, η) , due
to a crack is
0 
k ξ (θ )
kR = 
k η (θ)
 0
If the vibration due to any out-of-balance forces
acting on a rotor is greater than the static deflection
of the rotor due to gravity, then the crack will
remain either open or closed depending on the size
and location of the unbalance masses. If the crack
remains open, the rotor is then asymmetric and this
condition can lead to stability problems. Of more
importance in large machines is the situation where
the vibration due any out-of-balance forces acting
on a rotor is less than the deflection of the rotor due
(1)
Thus, in the hinge model the variable stiffnesses
k ξ (θ) and k η (θ) are defined by

k0
kξ ( θ ) = 
k
 ξ
607
−
π
π
≤θ<
2
2
3π
π
≤θ<
2
2
(2a)
608
P ROCEEDINGS OF ISMA2002 - VOLUME II
and
of the shaft material, the depth of the crack and the
lateral force.

k0
kη ( θ ) = 
k
 η
−
π
π
≤θ<
2
2
π
3π
≤θ<
2
2
(2b)
where θ is the angular position of the crack, k 0 is
the stiffness of the shaft when the crack is closed
and k ξ and k η are the stiffnesses in the ξ and η
directions when the crack is open. When the crack is
closed its stiffness is equal to that of an equivalent
undamaged shaft.
Whilst the hinge model might adequately
represent very shallow cracks Mayes and Davies
[5,6] proposed a model in which the opening and
closing of the crack was described by a cosine
function. Thus in the Mayes model, the stiffnesses
k ξ (θ) and k η (θ) are defined as follows:
(
) (
)
kξ ( θ ) = 12 k0 + kξ + 12 k0 − kξ cos θ
(3a)
and
(
) (
)
kη ( θ ) = 12 k0 + kη + 12 k0 − kη cos θ .
(3b)
Both the hinge model and Mayes model suffer from
the defect that there is no direct relationship
between the shaft stiffness and the depth of crack.
Jun et al. [8] proposed a crack model based on
fracture mechanics. The cross coupling stiffnesses,
as well as the direct stiffnesses are estimated as the
crack opens and closes. Thus, the stiffness matrix in
rotating coordinates becomes
 k ξξ (θ) k ξη (θ )
kR = 

k ηξ (θ ) k ηη (θ)
(4)
To determine the stiffnesses Jun et al. assumed
the crack to be at the mid-span of the shaft. The total
stress intensity factor is derived from the sum of the
stress intensity factors due to forces acting in the ξ
and η directions. Jun assumes that the sign of total
stress intensity at each point along the crack front
determines whether the crack is open or closed at
that point. The additional deflections in the ξ and
η directions due to the open crack are then
determined using Castigliano’s theorem. This
involves numerical integration over the area of the
proportion of the crack that is open. For a particular
angle of rotation, the Jun model is a function the
shaft length and diameter, the modulus of elasticity
3. Frequency components in a
crack model
3.1. The hinge model
The hinge model is defined by equation (2).
Expanding k ξ (θ ) and k η (θ) in a Fourier series, we
have
 k 0 + kξ 
kξ ( t ) = 

 2 
C3 C5 C7
4  k 0 − kξ  

+ 
+
−
+ …
  C1 −
π  2 
3
5
7

(5a)
 k0 + kη 
kη ( t ) = 

 2 
C3 C5 C7
4  k0 − kη  

+ 
+
−
+ …
  C1 −
π  2 
3
5
7

(5b)
where C1 = cos θ , C3 = cos ( 3θ ) , etc., and θ = Ωt .
Letting,
 k0 + kξ 
 k + kη 
 , k Mη =  0

k Mξ = 

 2 
2




(6a)
 k0 − kξ 
 k − kη 
 , k Dη =  0

k Dξ = 

 2 ,
 2 


(6b)
and
we have
(
kξ ( θ ) = k M ξ + π4 k Dξ C1 − 13 C3 + 15 C5 − 71 C7 + …
)
(7)
(
kη ( θ ) = k M η + 4 k Dη C1 − 13 C3 + 15 C5 − 17 C7 + …
π
)
This is the variation of stiffness in rotating
coordinates. Now the stiffness matrix in fixed
coordinates, K F is given by
K F = TT K R T
(8)
DYNAMIC MODELLING AND NUMERICAL MODELS
609
( K F )22 = 12 ( kM ξ + kM η ) − 14 ( kξ − kη ) C2
where the transformation matrix T is
S1 
C
T = 1

− S1 C1 
( ){( 34 kDξ + 83 kDη ) C1 + ( − 1528 kDξ + 158 kDη ) C3
92 k
8
+ ( 105
Dξ − 105 k Dη ) C5
+ ( − 188
k + 8 k
C + …}
315 Dξ 315 Dη ) 7
+ π1
(9)
Thus
K F = TT K R T
 kξ ( θ )
+ kη ( θ )
=
 k ( θ) − k ( θ) S C
1 1
η
 ξ
C12
S12
(
)
( kξ ( θ ) − kη ( θ ) ) S1C1 
kξ ( θ ) S12 + kη ( θ ) C12 

(10)
where S1 = sin θ , etc. Now, considering each of the
fixed co-ordinate stiffness matrix terms in turn,
beginning with the off diagonal one we have,
( K F )12 =
1
2
(
)
(
)
1
kM ξ − k M η S2 + k Dξ − k Dη ×
π
( 43 S1 + 54 S3 − 214 S5 + 454 S7 − 774 S9 …)
Since k Mξ − k Mη = 12 (k ξ − k η ) and
k Dξ − k Dη = − 12 (k ξ − k η ) , then
( K F )12 = ( kξ − kη ){ 14 S2 −
(
1 4
S + 4S
2π 3 1 5 3
4 S + 4 S − 4 S …
− 21
5 45 7 77 9
)}
(11)
Considering the diagonal terms:
(13)
3.2. The Mayes model
This model is defined, in rotating coordinates, by
equations (3a) and (3b). Thus
k η (θ) = k Mη + k DηC1
and
k ξ (θ) = k Mξ + k Dξ C1
The stiffnesses in fixed coordinates can be
developed using the same approach as for the hinge
model. Thus using equation (8) we have the
following stiffness expressions in fixed coordinates,
(K F )12 = 12 (k Mξ − k Mη )S 2 + 14 (k Dξ − k Dη )(S1 + S 3 )
or
( K F )12
)( 83 C1 + 158 C3 − 1058 C5 + …)
+ 12 k M η − 12 k M ηC2
+
(
1k
π Dη
)(
4 C − 28 C
3 1 15 3
92 C
+ 105
5
−…
)
and simplifying,
( K F )11 = 12 ( kM ξ + kM η ) + 14 ( kξ − kη ) C2
+ π1
{( 83 kDξ + 34 kDη ) C1 + ( 158 kDξ − 1528 kDη ) C3
(
)
8 k
188
+ ( 315
Dξ − 315 k Dη ) C7 + …}
8 k
92
+ − 105
Dξ + 105 k Dη C5
}
( S1 + S3 )
(14)
1
4
( (3kDξ + kDη ) C1 + ( kDξ − kDη ) C3 )
(15)
and
( K F )22 = 12 ( kM η + kM ξ ) − 14 ( kξ − kη ) C2
+
1
4
((3kDη + kDξ ) C1 + ( kDη − kDξ ) C3 )
(16)
Thus in fixed coordinates the Mayes model
generates a constant term plus 1X, 2X and 3X rotor
speed components in the diagonal stiffness terms
and 1X, 2X and 3X rotor speed components in the
off-diagonal stiffness terms.
3.3. The model of Jun et al.
(12)
Similarly,
1
2
( K F )11 = 12 ( kM ξ + kM η ) + 14 ( kξ − kη ) C2 +
= 12 k M ξ + 12 k M ξC2
(
){
It may be shown that
( K F )11 = kξ ( t ) C12 + kη ( t ) S12
+ π1 k Dξ
(
= 14 kξ − kη S2 −
In the model of Jun et al. [8] the stiffnesses k ξξ (θ) ,
k ηη (θ) and k ξη (θ) are determined numerically for
P ROCEEDINGS OF ISMA2002 - VOLUME II
a particular set of conditions. Thus the coefficients
of the Fourier expansion representing the stiffnesses
must be determined numerically using the least
squares approach. In rotating coordinates
k ξξ (θ) can be represented by the series
kξξ ( θ ) = b0 + b1 cos θ + b2 cos ( 2θ )
(17)
+ b3 cos ( 3θ ) + b4 cos ( 4θ ) +…
Normalised stiffness
610
0
180
θ (degrees)
270
360
are still cosine series and the off-diagonal
stiffness, (k F )12 is still a sine series.
Figures 1 and 2 show the stiffness variation with
angle of rotation for a 40% crack in rotating
coordinates, and also Figures 3 and 4 in stationary
coordinates. The markers indicate the simulated
values and the lines represent the fit of the Fourier
Series. Tables 1 to 4 show the coefficients of the
corresponding Fourier Series. The Fourier series
provides a good fit to the functions. In these figures
and tables, the stiffnesses have been normalised by
dividing them by the stiffness of the undamaged
shaft, k0.
1.0
1.0
0.95
k11
k22
0.9
0
90
180
θ (degrees)
270
360
Figure 3. Normalised stiffnesses (k F )11 and
(k F )22 - fixed coordinates
Normalised stiffness
The process of converting to fixed coordinates
using (8) changes the values of the coefficients in
these series but the diagonal stiffnesses, (k F )11 and
Normalised stiffness
(18)
+ a3 sin ( 3θ ) + a4 sin ( 4θ ) +…
Normalised stiffness
90
Figure 2: Normalised stiffness kηξ - rotating
coordinates
represented by the series
(k F )22
0.0
-0.02
The stiffness k ηη (θ) can be represented by a
similar series. The stiffness k ξη (θ) can be
kξη ( θ ) = a1 sin θ + a2 sin ( 2θ )
0.02
0.02
0.0
-0.02
0.95
kξξ
k
ηη
0.9
0
90
180
θ (degrees)
270
360
Figure 1: Normalised stiffnesses kξξ and kηη rotating coordinates
0
90
180
θ (degrees)
270
360
Figure 4: Normalised stiffness (k F )12 - fixed
coordinates
DYNAMIC MODELLING AND NUMERICAL MODELS
611
Crack depth = 10% diameter
Fourier
coefficient
b0
b1
b2
b3
b4
b5
Crack depth = 40% diameter
kξξ
kηη
kξξ
kηη
0.998186
0.002241
0.000001
−0.000583
−0.000001
0.000203
0.999880
0.000142
0.000000
−0.000023
0.000000
−0.000006
0.953773
0.051221
0.000845
−0.005069
−0.000201
−0.000711
0.979907
0.017422
−0.000034
0.004754
0.000170
−0.001401
Table 1: Fourier coefficients for the Jun model, rotating coordinates
Fourier
coefficient
a1
a2
a3
a4
a5
Crack depth = 10% dia
Crack depth = 40% dia
kηξ
kηξ
−0.000141
0.000000
0.000117
0.000000
−0.000078
−0.016949
−0.000364
0.006663
0.000232
0.000148
Table 2: Fourier coefficients for the Jun model, rotating coordinates
Crack depth = 10% diameter
Fourier
coefficient
b0
b1
b2
b3
b4
b5
Crack depth = 40% diameter
(k F )11
(k F )22
(k F )11
(k F )22
0.999034
0.001588
−0.000847
0.000242
0.000000
−0.000019
0.999033
0.000794
0.000848
−0.000848
−0.000008
0.000215
0.967241
0.045458
−0.012870
−0.000084
0.000022
0.000198
0.966438
0.023185
0.013681
−0.000231
−0.000053
−0.002311
Table 3: Fourier coefficients for the Jun model, fixed coordinates
Fourier
coefficient
a1
a2
a3
a4
a5
Crack depth = 10% dia.
Crack depth = 40% dia
(k F )12
(k F )12
0.000794
−0.000847
0.000363
0.000000
−0.000047
0.022711
−0.012859
−0.000123
0.000038
0.000497
Table 4: Fourier coefficients for the Jun model, fixed coordinates
612
P ROCEEDINGS OF ISMA2002 - VOLUME II
4. Equations of motion
The analysis may be performed in fixed or rotating
coordinates. If the bearings and foundations are axisymmetric then the stator dynamic stiffness will
appear constant in the rotating frame, and there is
some benefit in analysing the machine in rotating
coordinates. Typically foundations will be stiffer
vertically than horizontally, and in this case the
advantage in using rotating coordinates is
significantly reduced. Thus in this paper the
machine is analysed in fixed coordinates.
In
fixed
coordinates
we
have
K (θ) = K 0 − K D (θ) where K 0 is the diagonal
stiffness matrix for the undamaged beam and
K D ( θ ) is the stiffness change due to damage. Let
the deflection the system be q = q st + q u where
q st is the static deflection and qu is the deflection
due to whirl. Thus q = qu , and the equations of
motion in fixed coordinates is
Mqu + ( K 0 − K D ( θ ) ) ( q st + qu ) = Qu + W (19)
where Qu and W are the out of balance forces and
gravitational force respectively. Assuming that
K 0 K D ( θ ) , then
( K 0 − K D ( θ ) ) q st ≈ K 0q st = W
Fourier Series. Such parametrically excited systems
can exhibit a wide variety of responses, from
periodic with either super or sub harmonics, and
even a chaotic response [9-12]. Indeed the steady
state response may even be unstable [13]. However
most detection schemes for rotor cracks look for the
2X response, and orbits from experimental rigs
predominantly show super harmonic response. Thus
the response will be assumed to be periodic, with
the same frequency as the rotational speed. The
response may then be written as a Fourier series,
and the method of harmonic balance applied. The
steady state response is completely defined by the
coefficients of the Fourier Series, and the harmonic
balance approach [14,15] generates equations for
these coefficients by equating the sine and cosine
coefficients in the equations of motion, at all
frequencies. Since the Fourier Series is, in theory,
infinite, the series is truncated, and the coefficients
estimated using the lower frequency harmonic
balance equations.
The only slightly tricky part in applying the
harmonic balance method is to simplify the product
K D ( θ ) qu . Assuming that we know the Fourier
Series of K D ( θ ) ,
K D ( θ) = K D0
∞
+
∑
K Dci cos ( iθ ) + K Dsi sin ( iθ ) 
(21)
i =1
and the unbalance response is written as
and
Mqu + ( K 0 − K D ( θ ) ) qu = Qu
(20)
Thus far damping and gyroscopic effects have
been ignored. These effects may be included in
equation (20) and the assumption made here is that
the effect may be modelled as an additional velocity
term with a constant coefficient matrix. Thus the
skew-symmetric contribution to the stiffness matrix
from damping in the rotor is ignored, and this
damping is assumed to be axi-symmetric.
5. Harmonic balance analysis
The equations of motion (20) for the system are
linear with parametric excitation through the
variable stiffness K D ( θ ) , where θ = Ωt and Ω is
the rotor speed. With the crack opening and closing
based on the static beam deflection the stiffness is a
periodic function of time, and has been written as a
qu ( θ ) = q 0
∞
+
∑ q
j =1
cj cos (
jθ ) + q sj sin ( jθ )  ,
(22)
then K D ( θ ) qu involves products of sines and
cosines that are easily simplified using the standard
trigonometric formulae. Of course complex
exponentials could be used instead of sines and
cosines.
6. A Jeffcott rotor example
A Jeffcott rotor modelled using two degrees of
freedom will be used to demonstrate the effect of the
different crack models on the response of a rotor.
The rotor is 700 mm long, with a diameter of
15 mm, and has a 1 kg disc at the centre. The crack
is assumed to be located close to the disc and only
613
responses that are symmetrical about the shaft centre
are considered. The natural frequency of the rotor is
approximately 42 Hz the damping ratio is taken as
1%. Gyroscopic effects are neglected. A crack of
40% of the diameter is simulated, using the model
of Jun et al. [8]. This gives rise to a stiffness
reduction of approximately 9% when the crack is
open, and figures 1 to 4 give the stiffness variation
with rotor orientation. In order to compare the crack
modelling approaches, for the hinge and Mayes
models the reduction in direct stiffnesses due to the
open crack is taken directly from the results of the
model of Jun et al. [8]. An unbalance force arising
from a 10 g mass at 100 mm radius is used to excite
the rotor, and at t=0 the crack is closed and the
unbalance force is in the x direction.
Figure 5 shows the steady state synchronous
response magnitude for the three crack models.
Fourier Series terms up to 21 times the running
speed are included to ensure convergence of the
coefficients. All the responses are very close, and
only minor differences in magnitude occur near the
resonance frequency, which has reduced to
approximately 41 Hz. Figure 6 shows the steady
state 2X response as a function of rotor speed and is
clearly excited when the running speed is
approximately half of the resonance frequency.
Figure 7 shows the shaft orbit at 1240 rpm, where
the 2X component is well excited, and shows that
there is a significant quantitative difference between
the models. However the form of the orbits are very
similar.
x response (dB)
DYNAMIC MODELLING AND NUMERICAL MODELS
-20
-60
-100
y response (dB)
Hinge
Mayes
Jun et al.
0
-20
1000
2000
3000
1000
2000
Rotor speed (rpm)
3000
Hinge
Mayes
Jun et al.
-60
-100
0
x response (dB)
Figure 5: The steady state synchronous response (dB relative to 1 m) as a function of rotor speed
Hinge
Mayes
Jun et al.
-60
-80
-100
-120
y response (dB)
0
1000
2000
3000
1000
2000
Rotor speed (rpm)
3000
Hinge
Mayes
Jun et al.
-60
-80
-100
-120
0
Figure 6: The steady state 2X response (dB relative to 1 m) as a function of rotor speed
614
P ROCEEDINGS OF ISMA2002 - VOLUME II
y (mm)
1
Hinge
Mayes
Jun et al.
0
Machinery, Paper C178/76, IMechE, (1976), pp.
123-128.
[4] R. Gasch, A survey of the dynamic behaviour of
a simple rotating shaft with a transverse crack,
Journal of Sound and Vibration, Vol. 160, No. 2,
(1993), pp. 313-332.
[5] I.W. Mayes, W.G.R. Davies, A method of
calculating the vibrational behaviour of coupled
rotating shafts containing a transverse crack,
International Conference on Vibrations in
Rotating Machinery, Paper C254/80, IMechE,
(1980), pp. 17-27.
-1
-1
0
x (mm)
1
Figure 7. The steady state orbits at 1240 rpm for the
different crack models
7. Conclusions
The three crack models examined have relatively
little effect on the predicted steady state 1X
response (Figure 5) but they do have some influence
on the predicted whirl orbit and the steady state 2X
response (Figure 6 and 7). However, in any crack
identification scheme, these differences are not
likely to have a significant effect.
The advantage of using the model of Jun et al. to
represent a crack is that although it requires some
numerical integration, the stiffness variation is a
function of the crack depth.
This study is limited to a Jeffcott rotor and the
foundations are therefore inherently isotropic. The
next stage in the investigation is to incorporate the
crack models into a beam type element so that more
complex rotor models can be studied using finite
element analysis.
References
[1] S. Prabhakar, A.S. Sekhar and A.R. Mohanty,
Detection of cracks in a rotor-bearing system
using wavelet transforms. Mechanical Systems
and Signal Processing, Vol. 15, No. 2, (2001),
pp. 447-450.
[2] N. Bachschmid, P. Pennacchi, E. Tanzi, A.
Vania, Identification of transverse crack
position and depth in rotor systems, Meccanica,
Vol. 35, (2000), pp. 563-582.
[3] R. Gasch, Dynamic behaviour of a simple rotor
with a cross sectional crack, International
Conference on Vibrations in Rotating
[6] I.W. Mayes, W.G.R. Davies, Analysis of the
response of a multi-rotor-bearing system
containing a transverse crack in a rotor, Journal
of Vibration, Acoustics, Stress and Reliability in
Design, Vol. 106, No. 1, (1984), pp. 139-145
[7] B. Grabowski, Vibration behaviour of a turbine
rotor containing a transverse crack, ASME
Design Engineering Technical Conference, St
Louis, (1979), paper 79-DET-67.
[8] O.S. Jun, H.J. Eun, Y.Y. Earmme, C.-W. Lee,
Modelling and vibration analysis of a simple
rotor with a breathing crack, Journal of Sound
and Vibration, Vol. 155, No. 2, (1992), pp. 273290.
[9] I.W. Mayes, W.G.R. Davies, The vibration
behaviour of a rotating shaft system containing
a transverse crack, International Conference on
Vibrations in Rotating Machinery, Paper
C168/76, IMechE, (1976), pp. 53-64.
[10] N. Bachschmid, A. Vania, E. Tanzi, P.
Pennacchi, Statical and dynamical behaviour
of rotors with transversal crack - part I:
comparison between different models,
Proceedings of the XIV Congresso Nazionale
AIMETA, Como, (1999).
[11] P.C. Muller, J. Bajowski, D. Soffker, Chaotic
motions and fault-detection in a cracked rotor
chaotic rotor response, Nonlinear Dynamics,
Vol. 5, No. 2, (1994), pp. 233-254.
[12] P.N. Saavedra, L.A. Cuitino, Vibration
analysis of rotor for crack identification,
Journal of Vibration and Control, Vol. 8,
(2002), pp. 51-67.
[13] S.C. Huang, Y.M. Huang, S.M. Shieh,
Vibration and stability of a rotating shaft
containing a transverse crack, Journal of
Sound and Vibration, Vol. 162, No. 3, (1993),
pp. 387-401.
DYNAMIC MODELLING AND NUMERICAL MODELS
[14] N. Pugno, C. Surace, R. Ruotolo, Evaluation
of the non-linear dynamic response to
harmonic excitation of a beam with several
breathing cracks, Journal of Sound and
Vibration, Vol. 235, No. 5, (2000), pp. 749762.
615
[15] A. Tondl, T. Ruijgrok, F. Verhulst, R.
Nabergoj, Autoparametric resonance in
mechanical systems, Cambridge University
Press, 2000.
616
P ROCEEDINGS OF ISMA2002 - VOLUME II