Numerical implementation of structural dynamics analysis

Numerical implementation of structural dynamics analysis

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...

855KB Sizes 1 Downloads 99 Views

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).