Com+wrus & Swu~rures Vol. 65, No. I, pp. lO%125. 1997 1997 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0045-7949/97 117.00 + 0.00
,C
Pergamon 0045-7949(95)00338-x
NUMERICAL
IMPLEMENTATION OF STRUCTURAL DYNAMICS ANALYSIS K. H. Low
School of Mechanical and Production Engineering, Nanyang Technological University, Nanyang Avenue, Singapore 639798, Republic of Singapore (Recehed
9 December 1994)
Abstract-This work presents a software development and implementation on the frequency and response studies of vibratory systems, The systems can either be a discrete or continuous model, which can be represented by a set of coupled second-order differential equations. Several commonly-used numerical methods are used in the response analysis for user’s choices such as central difference, fourth-order Runge-Kutta, etc. The system frequencies are obtained by using a fast Fourier transform (FFT) method, giving data generated in the response analysis. Two types of the windows are used in the FFT analysis: rectangular and Hanning windows. The software implementation enables one to effectively analyze vibration problems on a screen environment. Examples of a discrete model and a continuous system are considered to illustrate the performance of the software. The results of response and frequencies are discussed and commented. 0 1997 Elsevier Science Ltd.
1. INTRODUCTION
system. Here [ml, [c] and [k] are the mass, damping and stiffness matrices of the system, and {Q(r)} is the external force vector. The variables {g}, {4} and {q} correspond to acceleration, velocity and displacement vectors of the system. Note that matrices are denoted by brackets and column matrices by braces. The superscript t denotes the transpose of a vector or matrix.
The present work presents a software implementation on the frequency and response analyse of vibratory systems. The systems can either be discrete or continuous model, which can be represented by a set of coupled second-order differential equations. Several numerical methods are used in the response analysis: central difference, Newmark Beta, Wilson Theta, Houbolt and fourth-order Runge-Kutta. Data generated in the response analysis are processed to obtain the system frequencies by using a fast Fourier transform (FFT) method. Two types of windows are sued and compared in the FFT analysis: rectangular and Harming windows. The software implementation with screen environment is briefly described. The simulation results of a four degrees of freedom model and a rotating beam system are presented to illustrate the performance of the software.
2.1. Four degrees
of freedom
model
As an example of mutidegree-of-freedom discrete systems, a four degrees of freedom model is considered here. A most general form of the equations of motion for the model is
2. EQlJATIONS OF MOTION If a set of generalized co-ordinated {g} = 141, q2, . . . , q,,}’ is used to describe the motion of an n-degree-of-freedom system shown in Fig. 1, the general form of the equations of motion can be written in matrix notation as
+
k,
bliB1 + [cl{41+ [Wq) = Bm.
(1) +
Equation second-order
(1) represents a set of n coupled differential equations of motion for the 109
kl2
;I-,,
1 1.
k14
kn
kx
kz,
k24
kx
kx
k,,
kx.,
k4,
ka
ka
k,
(2)
110
K. H. Low
Fig. I. Multidegree-of-freedom
2.2. Rotating
beam system
A continuous system consisting of a slender flexible beam attached to a rigid rotating shaft is shown in Fig. 2. Equations of motion of radially rotating beams were derived in Refs [l and 21. Euler-Bernoulli beam theory and the extended Hamilton’s principle [3] were utilized. The equations of motion and the associated boundary conditions for the present beam system are given by [l, 21
system
where 1 is the length of the beam, a is the radius of the rigid shaft, IR is the moment of inertia of the rigid shaft about the axis of rotation, Is is the moment of inertia of the beam, and p is the beam mass per unit length. In addition, EI denotes the flexural rigidity of the beam, and M the actuator torque applied at the motor shaft. An analytical solution to eqns (3x5) can be found by virtue of the assumed-mode method [3]. One assumes that
p {ti - w& + (a + x)0} + Eh””
M’(X,t) = i $,(x)q,(t),
+@(a
+ x)w’ = 0,
0 < x < 1,
(3)
I
-
f
’ {p(P - x2 + 2als0
+ 2p(P - x2 + 2a/w”(/, t) = 0,
(6)
/= I
- {fp8Z(I’ - x2) + @a([-x)}w”
where q,(t) are generalized coordinates, c#J~(x)are a set of shape functions, and n is the desired number of modes. Upon ignoring secondary terms, IV*, W’S’, (w’)~, w’ti’, in eqns (3) and (4), the assumed-mode method yields the following non-linear coupled set of ordinary differential equations
2ax)(w1)V
2ax)w’C’d)
dx = M,
w”‘(l, t) = 0,
(4)
,C, “fm’/ +824,(Cr, -
(5)
+S,B’=O,
Fig. 2. Rotating beam system.
m,,) + q,(k,, + k,T )) r=l,2
,...,
n,
(7)
Numerical
implementation
Fig. 3. Display
(Ze+
1,)B’+ -f
&S, = M,
of structural
dynamics
of the software
(8)
c,, =
,=I
analysis
111
window.
Jo’-
{fp(P - x2 + 2al-
2ax)cj:‘4,
‘re -p(a I
=
J
p&.4, dx,
k, =
0
+
xM,h#+}dx,
’EZ&“‘c#+dx,
i
J0
S, =
‘p(a
+ .u)$~ dx,
J0
= - EZl$:“(l)~,,l),
Fig. 4. Choices
for entering
parameters
p(a + x)? dx.
Z, =
JII
(9)
K. H. Low
112
Fig. 5. Display
It is convenient
of the response
to write eqns (7) and (8) in matrix
Fig. 6. Display
for a discrete
system
form (n = 3 for example)
of the deflection
for a rotating
as
beam system.
Numerical implementation
of structural dynamics analysis in
113
which
M,, = m,,
(i,j = l-3),
M/4 = I%, = S,,
(1 la)
(j = l-3),
(1 lb)
M44 = (ZB+ ZR),
(1 lc)
K,, = (c, - m,)B* + (k, + k$),
(i,j = l-3).
(1 ld)
Analytical integration by parts should be performed on the integrals appearing in eqn (11) in order to evaluate the matrix elements. The integration procedure is tedious and will not be produced here. General expressions for integrals containing mode shapes can be found in Ref. [4]. After evaluating the matrix elements in eqn (lo), M,, and &, the variables q,, q2, q, and t? are the four unknowns to be solved numerically. The beam deflection w(x, t) can then be evaluated by means of eqn (6). The Choleski method and the Gauss-Jordan matrix inversion scheme [S] were employed in the numerical algorithm, {q},+b, = [m-‘(Q). To overcome the difficulty of coupling, the variable 6 in eqn (1 Id) is always evaluated at one previous time point. It should be noted that this arrangement would not affect the accuracy of solutions provided that the time step used in the numerical analysis is small enough. 3. NUMERICAL
3.1. Integration
methods
A direct integration of the equations of motion provides the response of the system at discrete time intervals that are usually equally spaced.
Fig. 7. Flowchart for the FFT.
Fig.
ALGORITHMS
8. Display
of the FFT
result.
114
K. H. Low
Numerical
0
implementation
0.05
of structural
0.1
O.IJ
dynamics
0.2
analysis
115
0.25
0.3
time(s)
I 0
0.0s
0.1
aa
0.2
02S
0.3
time (s) Fig. 10. System response of clamped-free rotating beam with trapezoidal torque.
Determination of the response involves the computation of three basic parameters: displacement, velocity and acceleration. The integration algorithms are based on appropriately selected expressions that relate the response parameters at a given time interval to their values at one or more previous time points [6, 71. The following methods have been incorporated in the software implementation: (1) central difference, (2) Newmark Beta, (3) Wilson Theta, (4) Houbolt and (5) fourth-order Runge-Kutta. Algorithms of these methods, which are the standard ones used in the literature [&-lo], are omitted here for brevity. 3.2. Fast Fourier transform Natural frequencies can be determined by solving an n-degree algelbraic equation, which is the expansion of a characteristic matrix determinant.
Analytical solutions to the algebraic equation are almost impossible for n > 4. Another method to find frequencies is to express the system response in the frequency domain. A fast Fourier transform (FFT) module utilizing the Cooley-Tukey radix-2 decimation-in-frequency FFT algorithm [ 11, 121 is empoloyed to solve the frequency of the waveform. This algorithm solves for both the real and the imaginary parts of the FFT and then for the square root of the sum of the squares of the real and imaginary parts, namely the amplitude spectrum of the wave. A set of test signals is defined by the expression for the discrete Fourier transform
,4(n) =
$y ,=O
xo(t)W”‘,
(12)
116
K. H. Low
where n = 0, 1, . , N - 1. The time domain data are given by x0(O), .u,,(l), . , xo(N - I), and W is equal to ezn N. Computing the FFT of .uo(t) consists of log2 N = M stages. Each stage requires a pair of computations of the form, x,,, + ,(r) = x.,(r) + x,,,(s)
(13)
and
p=O
to N-
1;
m=O
1,
to M-
(14)
for specified integers Y and s, where x,,!(r) and x,,,(s) denote the real and imaginary values, respectively. The results at the end of mth computational stage are denoted by x.,(t), where t = 0, 1, . , N - 1. Also,
-2 ’ 0
the present algorithm is the in-place algorithm. This means that the current results replace the previous results. Thus x,,, + ,(t) overwrites x,,,(t) in going from stage m to stage m + 1. The subscript m merely defines as sequence of arrays that define the contents of an associated storage location at the end of the mth stage. After processing by the FFT algorithm, the output consists of a set of real and imaginary numbers. As such, the amplitude spectrum is defined as x,,, = Jx&)
+ x.,(s).
By taking a sample of the time domain data and then using the FFT to solve for the frequencies, a process called windowing is performed in the analysis [I 31. Two types of windows has been
I 0.2
0.15
0.05
time c”s; 1.s I
Fig.
I I. System response of pinned-free
05)
rotating
beam with trapezoidal
torque.
Numerical implementation of structural dynamics analysis
117
29.43529.4329.425-
29.4 -
0
0.05
0.1
0.1s the(s)
0.2
0.25
03
0.05
0.1
0.1s time(s)
0.2
0.25
03
116.95 116.94116.93116.92116.91-
116.63. 0
3161)
1
0
O.OS
0.1
tdl;r)
0.2
0.2S
0.3
Fig. 12. System frequencies of clamped-free rotating beam with trapezoidal torque.
K. H. Low
118
24s
24
23s
23
22s
22 21*
21 0
Fig. 13. System frequencies
___
om
of pinned-free
__
0.1
rotating
_ __
O.IJ
beam with trapezoidal
__
0.1
torque.
Numerical implementation of structural dynamics analysis
119
the (s)
0
0.01
0.02
time(s)
0.03
0.05
Fig. ‘14. Divergent system response of pinned-free rotating beam with h = 5 x 10d5 s.
incorporated into the program and their effects were studied by Low [13]. The first type is the simplest, namely the rectangular. Another is the Hanning window, which offers a good compromise between frequency resolution and reduced leakage. As illustrated in the work by Low [13], the Hanning window attenuates the results that give a narrower bandwidth around the peak if compared with those using the rectangular window. Therefore, the Hanning window is generally more suitable for capturing the vibratory waveforms. 4. SOFTWA.RE IMPLEMENTATION
Two programs were developed for the vibration analysis of discrete and continuous systems. The programs were written in Turbo Pascal language. Each program contains the same procedure and it is divided into five modules: CAS 6511-E
(1) Menu routine. This module establishes a window layout with a pop-up menu bar which ailows the user to select the option of operation. The option includes input of data, selection of numerical methods, graphic plotting and FFT (see Fig. 3). (2) Numerical routine. This module contains the procedures of different numerical methods and the FFT algorithm. (3) Data routine. This module involves the initial conditions and the input parameters. It also performs vector-matrix operations (see Fig. 4). (4) Matrix routine. This module contains matrix manipulations using the Choleski method and the Gauss-Jordan matrix inversion. (5) Graphics routine. In this module data are retrieved from the files and the maximum and minimum values are captured. The data, maximum and minimum values are transferred to the graphics subroutine to convert into screen coordinates. These
120
K. H. Low
24.5
24 I-
2x3
23
22.5
II ’ 0
OR4
0.01
time (s)
Fig. 15. Divergent system frequencies of pinned-free rotating beam with h = 5 x lO-J s.
121
Numerical implementation of structural dynamics analysis
time (s) 2.3
Fig. 16. System response of clamped-free rotating beam with sine torque.
coordinates are joined by drawing a line between them to produce a curve on the screen (see Figs 5 and 6). Also, the associated frequency spectra can be produced by using the FFT algorithm (see Figs 7 and 8). 0.1057 -0.0856 ICI =
5.1. Four degrees of, freedom model
The numerical values of the mass, damping and stiffness matrices of the four degrees of freedom system considered in Ref. [14] and are shown below: $ 0
0
0
oo+o’ [ml = [ 00 02 00
01
I
0.0904 -
5. SIMIJLATION RESULTS
-0.0856 0.0904 - 0.0563
-0.0563
1 (16
0.0742 - 0.3471 0.2762 0.3348 0.2762 -0.0942 0.0948 0.0742 -0.0942
’
The initial conditions are zero except for displacements being selected as {q(O)} = { 1.O, 0.5, 0.0, 0.1)’ and the forcing vector is zero. The initial acceleration is thus {d(O)] = [ml-‘({Q(O)} - [c]{tj(O)j - [k]{q(O)}) = { - 1.O, -0.125, 0.75, -0.125)‘. Each
K. H. Low
122
a9.a3s
a9.403
a9Mas
a9.a 29.601s
a9.401
a9a
a9.4 0
om
0.1
a15
aa
0.25
0.3
time (s)
ma46
116.661
Fig. 17. System frequencies of clamped-free rotating beam with sine torque.
123
Numerical implementation of structural dynamics analysis
a15
time (s)
0
O.OS
0.1s
0.1
0.2
0.25
0.3
time 6)
Fig. 18. System response of clamped-free rotating beam with cosine torque.
mass is given an initial displacement and released. The central difference method is used to solve for the displacements at
pinned-free models for the corresponding shape function d;(x) in eqn (6) is applied to investigate the effects of boundary conditions upon the beam vibrations. The system response was obtained by using the central difference method. Also, system frequencies for the first three modes were calculated using the scheme presented in Refs [15-171. For all the following case studies, it is assumed that I = I m, a= O.O5m, EI= 5.5Nmr and p =0.858 kg m-‘, /I = 0.3 where /I = &JIB. 5.2.1. Trapezoidal torque. The prescribed trapezoidal torque (M = I Nm) represents a servo motor being started, accelerated, moved in constant velocity, decelerated and stopped. Using the central difference method, the system was analyzed by modeling the beam attached to the motor-hub in two different models: clamped-free and pinned-free. The tip deflection and the joint velocity for the
124
K. H. Low
29.4017 29.Y)M 29.4013 29.4014 29.4013 26.4012 29.4011 29.401 29AOQ3 29au0 29AOo7
116.643
116A423
116.842
*l_/ 116.8415
I 116.841
0
O.OS
0.1
0.1s
0.2
a.23
0.3
316.848
316.8473
316.847
316.6463
Fig. 19. System frequencies of clamped-free rotating beam with cosine torque.
Numerical implementation of structural dynamics analysis and pinned-free models are shown in Figs 10 and 11, respectively. It is noted that the system frequencies are timedependent which follow the respective joint velocity profile. The different starting frequencies for the two models will thus aflect the system response. The joint velocity of the pinned-free model shown in Fig. 11 resembles closely the analytical trapezoidal torque profile. It implies that the effect of elastic vibration on the rigid body motion becomes less significant. A stability analysis in Ref. [8] suggested that numerical solutions of the central difference method for the free undamped vibration is unconditionally stable if h/z < l/n .= 0.318, where h is the time step and 7 is the period of vibration. In fact, Low [16] illustrated that the choice of time step for convergent numerical solutions of vibratory systems depends on {the highest frequency, which is associated with the shortest period. For example, the third frequency of the clamped-free case (see Fig. 12) is about 3 16.9 Hz, whereas the third frequency of the pinned-free model (see Fig. 13) is up to 7200 Hz. The lowest period (7,) is thus 3.156 x 1C~~3and 1.389 x 10-4s for the clamped-free and pinned-free cases, respectively. This is why the time step used in the pinned-free model is 1 x lo-’ s and not the 5 x lo-’ s used for the case of clamped-free. As shown in Figs 14 and 15, the solutions diverge near t = 0.0466 s with h = 5 x 1O-5s. It is seen that one can obtain different results for the same rotating beam system by using a different beam model. The question of using which beam model is beyond the scope of the present work. The problem could be studied by understanding more on the motor-beam set up [18]. 5.2.2. Sine and cosine torques. The prescribed torque to a sinusoidal function for the clamped-free beam is given by clamped-free
M(t) = A sin(2nfl),
where A = 0.5 and f=
which can be represented by a set of n coupled second-order differential equations. Several numerical methods are used in the response analysis for user’s choices: central difference, Newmark Beta, Wilson Theta, Houbolt and fourth-order RungeKutta. Data generated in the response analysis are processed to obtain the system frequencies by using a fast Fourier transform (FFT) method. The vibration problem can be analyzed on-line by virtue of the screen environment. A discrete four degrees of freedom system and a continuous rotating beam model are considered to illustrate the implementation of the software. Acknowledgement-The
author wishes to express his gratitude to Mr Ong Kian Huat for his invaluable help in carrying out the numerical simulations. REFERENCES
I. A. Yigit, R. A. Scott and G. A. Ulsoy, Flexural motion
2. 3. 4. 5. 6. 7.
9.
IO. Il.
and the prescribed lorque of the clamped-free beam to a cosine function is
12.
where A = 0.5 and f’= 10 Hz,
13.
(18)
14.
M(t) = A cos(27$),
in whichfis the forcing frequency generated by the rotating motor-hub. Numerical results associated with the two torque functions are shown in Figs 1619. It is noted from Figs 17 and 19 that the system frequencies are again similar to the corresponding joint velocity. They vary near the values of the stationary clamped-free beam: 29.4, I 16.84 and 3 16.84 Hz. 6. CONCLUDING REMARKS
A software has been developed to analyze the response and frequencies for the vibratory systems,
of a radially rotating beam attached to a rigid body. J. Sound Vi& 121, 201-210 (1988). _ K. H. Low. Eiaen-analvsis of a tie-loaded beam attached to a iota&g join; ASMEJ. Vibr. Acousr. 112, 497-500 (1990). L. Meirovitch, Analvricol Methods in Vibrations. Macmillan, New Yo& (1967). R. D. Blevins. Formulas for Natural Freauencv and Mode Shape. ‘Van Nostrand Reinhold, New ‘York (1979). C. F. Gerald, Applied Numerical Analysis. AddisonWesley, Reading, MA. (1978). R. W. Clough and J. Penzien, Dynamics of Structures. McGraw-Hill, New York (1975). K. H. Low, Convergence of the numerical methods for problems of structural dynamics. J. Sound Vibr. 150,
342-349 (1991). 8. J. L. Humar,‘Dynamics
10 Hz, (17)
125
of Structures. Prentice-Hall, Enelewood Cliffs. NJ (1990). S. 5. Rao, Mechanical Vibrations. Addison-Wesley, Reading, MA (1986). A. F. D’Souza and V. K. Garg, Advanced Dynamics: Mode&g and Analysis. Prentice-Hall, Englewood Cliffs NJ (1984). E. 0. Brigham, The Fast Fourier Transform and Its Applications. Prentice-Hall, Englewood Cliffs, NJ (1988). D. Brook and R. J. Wynne, Signal Processing: Principles and Applications. Edward Arnold, London (1988). K. H. Low, Displacement and frequency analyses of vibratory systems. Compuf. Srruct. 54, 743-755 (1995). C. Minas and D. J. Inman, Model improvement by pole placement methods. In: Proc. ASME Design Technical Conf. on Mechanical Vibration and Noise, Vibration Analysis: Techniques and Application, Montreal,
Canada, pp. 179-185 (1989). 15. F. S. Tse, E. M. Ivan and R. T. Hinkle, Mechanical Vibrarions: Theory and Applications, 2nd Edn. Allyn and Bacon, Boston (1978). 16. K. H. Low, On some numerical algorithm for the solution of mechanical vibrations. Compur. Struct. 42, 461470 (1992). 17. K. H. Low, A reliable algorithm for solving frequency equations involving transcendental functions. J. Sound Vibr. 161, 369-377 (1993).
18. K. H. Low and Michael W. S. Lau, Experimental investigation on the boundary condition of slewing beams using a high-speed camera system. Mech. Mach. Theory 30, 629-643 (1995).