- Email: [email protected]
THEORY OF VIBRATION | Impulse Response Function
THEORY OF VIBRATION/Impulse Response Function 1335
Tl
f 0
ÿKÿ1 AA KAB IBB
27
l
The generalized displacement vector xl in eqn [26] couples modal coordinates with boundary displacements. The original EOM for the lth substructure are assumed as: l Kl xl f l Ml x
28
The order of this equation is reduced by substituting eqn [27] and premultiplying the resulting equation by TTl . The reduced-order EOM are now expressed as: l Kl xl f l Ml x
xA xB xM l
number of substrictures in structure magnitude of jth component of structural response nodal vector of internal generalized displacements nodal vector of generalized displacements at boundary vector of modal coordinates eigenvector eigenvalue
See also: Basic principles; Theory of vibration, Energy methods; Theory of vibration, Equations of motion; Theory of vibration, Fundamentals; Theory of vibration, Variational methods.
29
where:
Further Reading Ml TTl Ml Tl Kl TTl Kl Tl f l Tl f l
30
The reduced-order EOM for the whole structure are given by: K x f M x
2
.. .
6 6f 6 lÿ1 6 0 6 F6 6 fl 6 0 6 6 .. 4 . fn
.. . 0
Ilÿ1;l 0 0 .. . 0
Bathe K±J (1982) Finite Element Procedures in Engineering Analysis. New Jersey, USA: Prentice-Hall. Craig RR (1981) Structural Dynamics. New York, USA: Wiley. Meirovitch L (1980) Computational Methods in Structural Dynamics. Maryland, USA: Sijthoff and Noordhoff. Weaver W Jr, Johnston PR (1987) Structural Dynamics by Finite Elements. New Jersey, USA: Prentice-Hall.
31
where M ; K and f are found via comformable summations as in eqn [24]. Eqn [31] for the reduced-order structure can be solved for x . Thereafter, xA for each substructure is computed from eqn [25]. The component mode synthesis is also a Rayleigh± Ritz method. The Ritz basis vectors are found as in eqn [7], where F now contains eigenvectors due to the inclusion of vibrational mode shapes. F is given by: 3 7 7 7 7 7 7 7 7 7 7 5
Nomenclature fj I
n uj
magnitude of jth force component identity matrix
32
Impulse Response Function R K Kapania, Virginia Polytechnic Institute & State University, Blacksburg, VA, USA Copyright # 2001 Academic Press doi:10.1006/rwvb.2001.0118
Response of structures to short-duration loads (i.e., the time duration of the external force is much shorter than the time period), is of great interest to a designer. Especially important is the limiting case of these loads, the loads whose time duration approaches zero but their total impulse, area under the force± time curve, stays finite. This happens because the response of a structure to such a load plays an important role both in describing the system and in acting as a building block towards finding the response of the system to any arbitrary excitation. Here we will describe the methods to determine the impulse response function, defined as the response of a single-degree-of-freedom spring-mass-damper system with zero initial displacement and velocity to a unit impulse.
1336 THEORY OF VIBRATION/Impulse Response Function
The governing equation of a single-degree-of-freedom system of mass m, stiffness k, damping c, and subjected to an external force f t is given as: m x cx_ kx f t
1
where x denotes the displacement from the static equilibrium, and a dot denotes the derivative with respect to t. Mathematically, the impulse applied at a given time t t is described using Dirac's delta function, and is represented as d t ÿ t. Dirac's delta function belongs to a class of functions called generalized functions. To understand this function, consider a uniform force f t acting at t t over duration Dt such that the area under the curve, termed impulse I and given as I f Dt, remains unity. Now, as the value of the load duration Dt is reduced, the value of the force f will increase so that the product f Dt remains unity (Figure 1). As Dt approaches zero, the force f approaches infinity, but the area under the force±time curve remains unity. This function, with value approaching infinity and acting over a vanishingly small time (the independent variable) is called Dirac's delta function. This function, denoted by d t, has the following important properties: Z1 d t ÿ t dt 1
2
ÿ1
Z1 g td t ÿ t dt g t
3
ÿ1
Here g t is an arbitrary function of time. Consider the single degree of freedom system of eqn [1], and subjected to an impulse of magnitude I at t 0. The governing equation then becomes: m x cx_ kx Id t
The response of the system to a unit impulse, termed the impulse response function, is represented as h t. It is the solution of the above equation with _ vimp 1=m and initial conditions x 0 x 0 0, i.e. both initial velocity and initial displacement are zero. In terms of h t, the above equation becomes: 1 h 2zon h_ o2n h d t m h 0 0; and h_ 0 0
6
A number of methods can be used to obtain h t. Here, we present two classical methods: first, the Laplace transformation method; and second the Fourier series and integral method. Both methods are extensively used in determining the response of a system to an arbitrary excitation.
Laplace Transformation Method The Laplace transform of a function g t, written as G s is defined as:
4
Dividing both sides of the above equation by the mass m, we obtain: x 2zon x_ o2n x vimp d t
Figure 1 Dirac's delta function. The area under the load±time duration remains unity.
5
p where on k=m is the natural frequency of the system in radiansp per second, z c=cc is the damping ratio with cc 2 km; vimp I=m. Note that cc is called critical damping. Physically vimp is the change in the velocity of a mass m when subjected to an impulse Id t.
Z1 G s Lg t
eÿst g t dt
7
0
The function g t is defined for t < 0 and s is a complex variable. It is assumed that the function g t is such that the above integral exists. The Laplace transform is a linear operator, i.e.: Lag t br t aLg t bLr t aG s bR s
8
THEORY OF VIBRATION/Impulse Response Function 1337
Here a and b are constants, and g and r both are functions of t. Other properties of importance to us here are: Lg_ t sG s ÿ g 0 Lg t s2 G s ÿ sg 0 ÿ g_ 0
9
Knowing the Laplace transform of a function, the original function is obtained by taking the inverse Laplace transform: 1 g t L G s 2pi ÿ1
gi1 Z
G s est ds
10
gÿi1
p Here i ÿ1. Often tables of Laplace transform pairs, widely available in many books on engineering or operational mathematics, are used to obtain the Laplace transform of a function of interest. Tables are also used to obtain the inverse transform. Determining the inverse Laplace transform of a function, not given in the tables of Laplace transform pairs, may be a rather difficult task. Laplace transform of some functions that are important in the theory of vibrations are given in Table 1. The impulse response function h t can now be obtained by taking the Laplace transform of both sides of eqn [6]: s2 H s ÿ sh 0 ÿ h_ 0 2zon sH s ÿ h 0 o2n H s
1 11 m
_ Substituting h 0 0, and h 0 0, and combining all the coefficients of H s, we obtain: ÿ2 1 s 2zon s o2n H s m
12
Note that the function H s is called the transfer function for the given system and it represents the ratio of the output of the system to the input to the system in the Laplace domain. It is given as: 1 H s ÿ 2 m s 2zon s o2n
13
The impulse response function h t is obtained by taking the inverse Laplace transform of H s, i.e. h t Lÿ1 H s. We can obtain the inverse transform of H s by using partial fractions. To that end, we need to factor the denominator, i.e. represent the
denominator in the form (s ÿ s1 )(s ÿ s2 ), where s1 and s2 , the two roots of the quadratic equation, s2 2zon s w2n 0, are called the poles of the transfer function. The poles, s1 , and s2 are given as: s1;2 ÿzon on
qÿ z2 ÿ 1
14
Thus, the nature of the poles depends upon the value of the nondimensional damping parameter z. The two roots would be complex and distinct if z < 1 (underdamped case), equal and real negative if z 1 (critically damped), and distinct and real negative if z < 1 (overdamped). These three cases lead to three completely distinct type of systems as is shown in the following. Underdamped case (z < 1): The two, complex, poles are: s1;2 ÿzon ion
qÿ
1 ÿ z2 ÿzon iod
15
p where od on 1 ÿ z2 is the so-called damped natural frequency. The transfer function in terms of the partial fractions is given as: 1 1 A B ÿ m s ÿ s1 s ÿ s2 m s2 2zon s o2n
16
Table 1 Laplace transform of some functions that are important in the theory of vibrations Ld tÿt eÿts Lu t L t n
1 s
n! sn1
L eat
L sin ot
o s2 o2
L cosh ot
n 0; 1; 2; . . .
1 s ÿ a 3
L sinh ot
s s2 ÿ o2
o s 2 ÿ o2
L eat cos ot
sÿa s ÿ a2 o2
2 t Z 1 L4 g tdt5 G s s
L eat sin ot
Leat g t G sÿa
L eat cosh ot
0
L cos ot
s s2 o2
o s ÿ a2 o2
L eat sinh ot
sÿa s ÿ a2 ÿo2 o
s ÿ a2 ÿo2
1338 THEORY OF VIBRATION/Impulse Response Function
Multiplying both sides by the denominator in the above equation, we obtain: A s ÿ s2 B s ÿ s1 1
17
Substituting s s1 and s s2 , respectively, we obtain: A ÿB
1 1 i ÿ s1 ÿ s2 2iod 2od
18
p where od on z2 ÿ 1. Note that the impulse response function for the critical and overdamped case are nonoscillatory in nature. This can be seen in Figure 2B. In both Figures 2A and 2B, the plot of nondimensional displacement mon h t=I is given as a function of nondimensional time on t. Observe that as the value of the damping is increased from the critical value, the value of the nondimensional peak amplitude decreases but the time it takes for the response to die down increases.
The impulse response function, for the underdamped case, is thus given as: h t Lÿ1 H s; t > 0 ÿi 1 1 ÿ1 L ÿ 2mod s ÿ s1 s ÿ s2 Recalling Lÿ1 1= s ÿ a eat , and s1; 2 ÿzon iod , we obtain: ÿi es1 t ÿ es2 t 2mod ÿi ÿzon t iod t e e ÿ eÿiod t 2mod 1 ÿzon t e sin od t; t > 0 mod
h t
19
Here we have made use of the Euler formula eiy cos y i sin y. Note that the impulse function for the underdamped case, because of the two poles being complex conjugate for this case, is oscillatory in nature. This is shown in Figure 2A. Critically damped case (z 1): Both the roots are equal s1; 2 ÿon . The impulse response function is given as: " h t Lÿ1
1 m s on 2
#
teÿon t ; m
t>0
20
Overdamped case (z < 1): Both the roots are distinct and are given as: s1; 2 ÿzon p on z2 ÿ 1 ÿzon od . The impulse response function becomes: 1 h t L m s ÿ s1 s ÿ s2 1 ÿzon t e sin hod t ; t > 0 mod ÿ1
21
Figure 2 Nondimensional response of a single-degree-of-freedom spring-mass-damper system of mass m, and natural frequency on, to an impulse I applied at t 0 for different values of nondimensional damping z. (A) underdamped systems z < 1; (B) critically z 1 and overdamped systems z > 1. Note that the impulse response function h t is simply x t=I.
THEORY OF VIBRATION/Impulse Response Function 1339
Overdamped case
Response to Initial Conditions In the preceding section, we obtained the response of a dynamic system to a unit impulse by treating it as a special forcing function applied at t 0. Since the effect of applying a unit impulse is to impart a sudden velocity vimp 1=m at t 0, the impulse response function can also be obtained by studying the response of the system under an initial velocity _ x 0 1=m but keeping both the force f t and the initial displacement x0 equal to zero. In this section, we present the response of a second-order system under the influence of nonzero initial conditions: _ v0 . The impulse response x 0 x0 , and x 0 function can then be obtained by substituting x0 0, and v0 1=m. The response under these initial conditions, in the absence of f t, can be obtained by taking the Laplace transform of both sides of eqn [1]: s2 X s ÿ sx 0 ÿ x_ 0 2zon sX s ÿ x 0 o2n X s 0
sx0 v0 2zon x0 s2 2zon s o2n
26 For the special case of x0 0 and v0 1=m, all the three equations reduce to the respective expressions for the impulse response function given in eqns [19]± [21], respectively. Most of the systems of practical interest are underdamped systems where the nondimensional damping factor z is of the order of 0.02±0.05. The response due to the initial conditions (as was also seen for the case of the impulse response function in Figures 2A and 2B) are oscillatory in nature only for the underdamped case.
Fourier Series and Transform 22
The Laplace transform, X s, of the desired response becomes
X s
v0 zon x0 x t eÿzon t x0 cos od t sinh o t d od
23
The response x t is obtained using inverse Laplace transform by keeping in mind that three distinct cases (underdamped, critically damped, and overdamped) arise, depending on the value of the nondimensional damping parameter z. The response for the three cases, obtained using inverse Laplace transform (see Table 1), is given as Underdamped case v0 zon x0 sin od t x t eÿzon t x0 cos od t od 24
For many practical cases, we are interested in finding the response of the system to external excitation, such as buildings subjected to wind and earthquake loads, an aircraft wing subjected to gusts, an automobile on uneven pavements and so on. The impulse response function derived above can be used to obtain the response of any system to any external excitation, say f t. The response can be written as a sum of a complementary solution and a particular solution: x t xc t xp t
The complementary part xc is the solution of the governing equation with the right-hand side equal to zero, the so-called homogeneous equation, as was done in the previous section on determining the response of the solution to initial conditions. The complementary part with arbitrary constants, due to the presence of eÿzon t term, is transitory in nature and approaches zero as t increases. While determining the response due to loads that are of longer duration, the complementary solution is often ignored and emphasis is placed only on the particular solution xp t. This part of the solution, which does not have any arbitrary constants, can be obtained using the convolution (Duhamel's) integral:
Critically damped case x t eÿon t x0 v0 on x0 t
27
Zt 25
xp t f h t
f th t ÿ t dt 0
28
1340 THEORY OF VIBRATION/Impulse Response Function
Performing the above integral, at times, may be quite a difficult task, especially when the forcing function f t is not a simple function. Numerical and transform methods, Laplace and Fourier, are used in such cases. Using transform methods we can transform the above integral into an algebraic product which is easy to compute. Then, the solution in the time domain is obtained by performing the inverse transform. In previous years, the Laplace transform was the transform of choice for many engineers. However, its use is limited to only those functions whose Laplace transform is easy to obtain and also the inverse transform of the product [H(s)F(s)] is easily available. Moreover, determining the Laplace transform of a function and determining the response as function in time by performing inverse Laplace transformation does not lend itself easily to the tremendous power of a modern digital computer. It should, however, be mentioned that the availability of symbolic manipulators like Mathematica and MAPLE has considerably improved our ability to obtain inverse Laplace transform. (Caution: an effective way to compute inverse Laplace transforms of long or complex expressions in Mathematica is to first produce the partial fraction expansion in the Laplace domain and then to take the inverse of each resulting term separately. Also, it has been our experience that if inverse Laplace transforms of long expressions are taken without first performing the partial fraction expansion, one may get incorrect results.) It is because of its ability in using the power of modern computers that the Fourier transform has gained a widespread popularity for determining the system response as well as in determining the system model using experimental methods. Moreover, since the Fourier transform works in the frequency domain, it provides an additional advantage: the system parameters for many systems of practical interest can be described easier in the frequency domain than in the time domain. To understand the use of Fourier transform in vibration theory it is important to understand the use of Fourier sine and cosine series in expressing periodic, not necessarily harmonic, functions, and to understand an extension of these series for aperiodic functions, called Fourier integrals. The Fourier series for a function f t of period 2T can be written as: 1 1 X a0 X ak cos ok t bk sin ok t; 2 k1 k1 2p k ok 2T
f t
29
In general, one only needs few terms in the series to achieve a good accuracy with respect to the given function. Here ok , called harmonics, are the discrete frequencies, and the Fourier coeficients ak and bk in eqn [29] are given by: 1 ak T 1 bk T
ZT f t cos ok t dt
k 0; 1; 2 . . .
ÿT ZT
30 f t sin ok t dt
k 1; 2; 3 . . .
ÿT
Note that the coefficients ak vanish for odd functions f t ÿf ÿ t and the coefficients bk vanish for even functions f t f ÿ t. The Fourier series converges to f t for all t if f t is continuous for 0 t 2T, and converges to f t f tÿ =2 if there is a discontinuity at t. It can be shown that the coefficients obtained by eqn [30] yield the best approximation for a given number of terms used in the Fourier expansion. That is, for a given number of terms used in the Fourier series, the square of the error between the function and its Fourier representation is a minimum if the coefficients given by eqn [30] are used in eqn [29]. This implies that the accuracy of the Fourier series can only be improved by adding additional terms in the Fourier series. In vibration theory, especially in experimental structural dynamics, a complex representation is often used. It can be shown that for a periodic function f t, the Fourier representation in terms of complex exponential harmonic function, eiok t can be written as: f t
X 1 k1 f ok eiok t Do 2p kÿ1
31
where, Do, the difference between two consecutive discrete frequencies ok and ok 1 is given as: Do ok1 ÿ ok
2p 2T
32
In eqn [31], f ok is given as: f ok
ZT
f t eÿiok t dt
33
ÿT
Note that in the complex representation, the integer k varies from ÿ1 to 1. The negative frequencies are thus also included. But this representation is convenient in extending the notion of Fourier series to Fourier transform as described below.
THEORY OF VIBRATION/Impulse Response Function 1341
The Fourier series described above can be used to represent an aperiodic function, since such a function can be considered as a periodic function with the period 2T approaching 1. For this case, Do!0. Thus, instead of a discrete representation in the frequency domain, one gets a continuous variation over o. The Fourier series, in eqn [31], becomes an integral, called the Fourier integral, as given below: 1 f t 2p
Z1
f oeiot do
f o
The transfer function in the frequency domain H o is the Fourier transform of the impulse response function. For the underdamped case the result obtained by, using eqns [34] and [19], is given as: Z1 H o
In the above equation, f o is given as: f teÿiot dt
ÿ1
Functions f t and f o are considered to form a transform pair, the Fourier transform pair. Note that, in some references on Fourier transform, the factor 2p, seen here in the denominator of eqn [34] may be placed in the denominator of eqn [35] instead p and still in some other sources it may be placed as 2p in the denominator of both eqns [34] and [35]. Also note that, unlike the case here, some references define eqn [34] as the Fourier transform and eqn [35] as the inverse Fourier transform whereas in other references the names of the two equations are reversed. As is the case for the Laplace transform, the Fourier transform, denoted as F: t, is a linear operation: Fag t br t aFg t bFr t aG o bR o
36
Also, as is the case for Laplace transform, the Fourier transform of a convolution integral in the time domain becomes an algebraic product of the respective transforms in the frequency domain. That is: 2 t 3 Z F g h t F4 g th t ÿ t dt5 0
Fg tFh t G oH o
h teÿiot dt
ÿ1
35
37
where G o is the Fourier transform of g t and H o is the Fourier transform of h t as defined by eqn [35]. To obtain gh t from G oH o, one needs to apply the inverse Fourier transform given by eqn [34]. That is:
38
ÿ1
34
ÿ1
Z1
g h t Fÿ1 G oH o Z1 1 G oH oeiot do 2p
1 mod
Z1
39 eÿzon t sin od teÿiot dt
0
Here we have used the fact that h t 0 for t 0. The integral on the right-hand side of this equation can be easily obtained as: 1 ÿ m o2n 2izoon ÿ o2 "ÿ # 1 1 ÿ b2 ÿ i 2zb k ÿ1 ÿ b2 2 2zb2
H o
40
where b o=on . It is interesting to note that the transfer function in the frequency domain is the same as the steady-state (i.e. the transient part of the response containing eÿzon t term has died down) response of the system to a harmonic exponential excitation of unit amplitude and excitation frequency o (i.e. f t eiot in eqn [1]) and is called the frequency response function (FRF). The frequency of the steady-state response will be the same as that of the excitation. The complex nature of H o implies that there exists a phase difference (f) between the harmonic excitation and the resulting harmonic steadystate response. This steady-state response, xss o, is given as: xss o H oeiot jH ojei otÿf " #1=2 1 1 ei otÿf k ÿ1 ÿ b2 2 2zb2
41
Here f represents the phase difference between the harmonic excitation and the resulting steady-state response. This phase angle is given as: ÿ f tanÿ1 2zb 1 ÿ b2
1342 THEORY OF VIBRATION/Impulse Response Function
Also, since 1=k represents the response of the system under a unit static force, the term inside the square brackets in eqn [40] can be considered as the magnification of the static response if the unit load is instead applied in a harmonic manner with frequency o. The term, inside the square brackets, is called the magnification factor. In some references, the transfer function (or the FRF) given above is denoted as H io. The impulse response function can be obtained by taking the inverse Fourier transform of H o as shown in the following: h t Fÿ1 H o Z1 42 1 1 ÿ eiot do 2p m o2n 2izoon ÿ o2
z o and the semicircular arc (denoted as S) of radius r with r ! 1; as: Ic
1 2pm
1 2pm
Z1 ÿ1
Z S
ÿ
eiot do o2n 2i zoon ÿ o2
eizt ÿ dz 2 on 2i zzon ÿ z2
The second integral on the right-hand side of the above integral approaches zero as the integrand contains z2 term in the denominator. Comparing eqns [44] and [45], we obtain the impulse response function to be given as:
ÿ1
The integral on the right of the preceding equation can be performed by using Cauchy's relation, namely: I C
f z dz 2pi f z0 z ÿ z0
Here C is any simple closed path around z0 which is traversed in the anticlockwise direction and f z is analytic within C. To apply this relation to the integral in eqn [42], consider the following integral, obtained by replacing o with z, a complex variable, in eqn [42]: 1 Ic 2pm
I C
ÿ
eizt dz o2n 2i zzon ÿ z2
43
Consider the path C to be the upper half plane, i.e., a semicircle of radius r as r ! 1 (Figure 3). Realizing that the denominator in the integrand in eqn [43] when equated to zero has two complex roots: z1; 2 izon od , the integral becomes: Ic ÿ
1 2pm
I C
eizt dz 2od z ÿ i zon ÿ od
I 1 eizt dz 2od z ÿ i zon od 2pm C h i 2pi ÿei izon od t ei izon ÿod t 4pmod 1 ÿzon t e sin od t mod
h t
1 ÿzon t e sin od t mod
46
The transfer function in the frequency domain (or the FRF) H o and the impulse response function h t thus form a Fourier transform pair. In determining the impulse response function, the inverse Fourier transform was performed using an analytical approach that used the so-called Cauchy relation (or contour integral) to perform the improper integral. Such an approach is limited only to simple functions for which such integrals can be easily obtained either using Cauchy's relation or from the tables of Fourier and inverse Fourier transform (similar to the case for the Laplace transform). However, the Fourier transform, using its discretized version, can be performed numerically for all sorts of functions using the fast Fourier transform (FFT).
Fast Fourier Transform Let fn represent the value of a function at a discrete point tn nDt where Dt 2T=N, N being the total number of discrete points in the time period 2T or the time-history duration. Similarly fm represents the Fourier coefficient corresponding to the mth discrete frequency om given as om 2pm= 2T. Values of fn and fm are given as:
44
The integral in eqn [43] can also be written, by splitting the contour C in two parts: the real axis
45
fn
ÿ1 X 1N fm exp i 2pmn=N N m0
n 0; 1; 2; 3 . . . ; N ÿ 1 Nÿ1 X fn exp i 2pmn=N fm n0
m 0; 1; 2; 3 . . . ; N ÿ 1
47
THEORY OF VIBRATION/Impulse Response Function 1343
The preceding set of equations constitutes the discrete Fourier transform (DFT) pair and can easily be computed. However, if the computations are done using these equations in a direct manner, a transform will take N 2 complex multiplications. A fast method to evaluate the DFT was proposed by Cooley and Tukey. This algorithm called the fast Fourier transform (FFT), needs only O Nlog2 N complex multiplications if N 2g , where g is an integer. This results in a substantial reduction in computational time. For example, for N 1024 210 , the direct evaluation to calculate the Fourier transform will take 10242 1 048 576 complex multiplications. On the other hand, one requires only O 10 240 multiplications if FFT is employed, around 1% of the effort required otherwise. To understand the logic behind FFT, let us consider any one of the two equations in eqn [47], say the second one. It can be rewritten as: fm
Nÿ1 X
fn exp ÿi 2pmn=N
n0
Nÿ1 X
fn W mn
n0
48
m 0; 1; 2; 3 . . . ; N ÿ 1 Here W expÿi 2p=N. The FFT is based on the fact that W 0 W N W 2N . . . 1, and W N m W m . As a result, for a given value of g N 2g , the DFT can be performed in g stages. This is easily understood for N 4 i:e: g 2. For this case, eqn [48], becomes: 8 9 2 W0 > f > > = 6 0 < 0 > f1 6W 2 > 4 W 0 > f > ; :> W0 f3
W0 W1 W2 W3
W0 W2 W4 W6
3 W0 W3 7 7 W6 5 W9
8 9 f0 > > > = < > f1 f > > > ; : 2> f3
N 4 case, we calculate the DFT, once again using N complex multiplications, as: f0 g0 W 0 g1 f2 W 2 g1 g0 f1 g2 W 1 g3 f3 W 3 g3 g2 Note that the values of f are scrambled. The same approach can be easily implemented for other values of g. For example, for N 8 g 3 we will require three stages and each stage will have eight complex multiplications. The computer programs capable of performing FFT are readily available. For example, the commercially available packages Mathematica and MATLAB both contain commands for obtaining the FFT.
Nomenclature C f I v X(s) d(t) z
closed path force impulse velocity Laplace transform Dirac's delta function phase angle damping factor
See also: Critical damping; Transform methods.
49
Obviously a straightforward matrix-vector multiplication will take 42 complex multiplications. If, however, we make use of the fact that W 0 W 4 1 and W 6 W 2 , etc., we can reduce the number of multiplications to eight by performing the above multiplication in two g stages. In the first stage, using N complex multiplications, we calculate an intermediatory vector g whose components are given as: g0 f 0 W 0 f 2 g1 f 1 W 0 f 3 g2 W 2 f2 f0 g3 W 2 f3 f1 In the second (gth or last) stage, also the last stage for
Further Reading Clough RW, Penzien, J (1993) Dynamics of Structures, 2nd edn. New York: McGraw Hill. Cooley JW, Tukey JW (1965) An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation 19(90): 297±301. Craig RR (1981) Structural Dynamics, An Introduction to Computer Methods. New York: Wiley. Gerald CF, Wheatley PO (1989) Applied Numerical Analysis, 4th edn. Reading, MA: Addison Wesley. Inman DJ (1996) Engineering Vibration. Englewood Cliffs: Prentice Hall. Kelly SG (1993) Fundamentals of Mechanical Vibrations. New York: McGraw Hill. Kreyszig E (1993) Advanced Engineering Mathematics, 7th edn. New York: Wiley. Meirovitch L (1997) Principles and Techniques of Vibrations. Upper Saddle River: Prentice Hall. Rao, SS (1995) Mechanical Vibrations, 3rd edn. Reading, MA: Addison Wesley.
Tl
f 0
ÿKÿ1 AA KAB IBB
27
l
The generalized displacement vector xl in eqn [26] couples modal coordinates with boundary displacements. The original EOM for the lth substructure are assumed as: l Kl xl f l Ml x
28
The order of this equation is reduced by substituting eqn [27] and premultiplying the resulting equation by TTl . The reduced-order EOM are now expressed as: l Kl xl f l Ml x
xA xB xM l
number of substrictures in structure magnitude of jth component of structural response nodal vector of internal generalized displacements nodal vector of generalized displacements at boundary vector of modal coordinates eigenvector eigenvalue
See also: Basic principles; Theory of vibration, Energy methods; Theory of vibration, Equations of motion; Theory of vibration, Fundamentals; Theory of vibration, Variational methods.
29
where:
Further Reading Ml TTl Ml Tl Kl TTl Kl Tl f l Tl f l
30
The reduced-order EOM for the whole structure are given by: K x f M x
2
.. .
6 6f 6 lÿ1 6 0 6 F6 6 fl 6 0 6 6 .. 4 . fn
.. . 0
Ilÿ1;l 0 0 .. . 0
Bathe K±J (1982) Finite Element Procedures in Engineering Analysis. New Jersey, USA: Prentice-Hall. Craig RR (1981) Structural Dynamics. New York, USA: Wiley. Meirovitch L (1980) Computational Methods in Structural Dynamics. Maryland, USA: Sijthoff and Noordhoff. Weaver W Jr, Johnston PR (1987) Structural Dynamics by Finite Elements. New Jersey, USA: Prentice-Hall.
31
where M ; K and f are found via comformable summations as in eqn [24]. Eqn [31] for the reduced-order structure can be solved for x . Thereafter, xA for each substructure is computed from eqn [25]. The component mode synthesis is also a Rayleigh± Ritz method. The Ritz basis vectors are found as in eqn [7], where F now contains eigenvectors due to the inclusion of vibrational mode shapes. F is given by: 3 7 7 7 7 7 7 7 7 7 7 5
Nomenclature fj I
n uj
magnitude of jth force component identity matrix
32
Impulse Response Function R K Kapania, Virginia Polytechnic Institute & State University, Blacksburg, VA, USA Copyright # 2001 Academic Press doi:10.1006/rwvb.2001.0118
Response of structures to short-duration loads (i.e., the time duration of the external force is much shorter than the time period), is of great interest to a designer. Especially important is the limiting case of these loads, the loads whose time duration approaches zero but their total impulse, area under the force± time curve, stays finite. This happens because the response of a structure to such a load plays an important role both in describing the system and in acting as a building block towards finding the response of the system to any arbitrary excitation. Here we will describe the methods to determine the impulse response function, defined as the response of a single-degree-of-freedom spring-mass-damper system with zero initial displacement and velocity to a unit impulse.
1336 THEORY OF VIBRATION/Impulse Response Function
The governing equation of a single-degree-of-freedom system of mass m, stiffness k, damping c, and subjected to an external force f t is given as: m x cx_ kx f t
1
where x denotes the displacement from the static equilibrium, and a dot denotes the derivative with respect to t. Mathematically, the impulse applied at a given time t t is described using Dirac's delta function, and is represented as d t ÿ t. Dirac's delta function belongs to a class of functions called generalized functions. To understand this function, consider a uniform force f t acting at t t over duration Dt such that the area under the curve, termed impulse I and given as I f Dt, remains unity. Now, as the value of the load duration Dt is reduced, the value of the force f will increase so that the product f Dt remains unity (Figure 1). As Dt approaches zero, the force f approaches infinity, but the area under the force±time curve remains unity. This function, with value approaching infinity and acting over a vanishingly small time (the independent variable) is called Dirac's delta function. This function, denoted by d t, has the following important properties: Z1 d t ÿ t dt 1
2
ÿ1
Z1 g td t ÿ t dt g t
3
ÿ1
Here g t is an arbitrary function of time. Consider the single degree of freedom system of eqn [1], and subjected to an impulse of magnitude I at t 0. The governing equation then becomes: m x cx_ kx Id t
The response of the system to a unit impulse, termed the impulse response function, is represented as h t. It is the solution of the above equation with _ vimp 1=m and initial conditions x 0 x 0 0, i.e. both initial velocity and initial displacement are zero. In terms of h t, the above equation becomes: 1 h 2zon h_ o2n h d t m h 0 0; and h_ 0 0
6
A number of methods can be used to obtain h t. Here, we present two classical methods: first, the Laplace transformation method; and second the Fourier series and integral method. Both methods are extensively used in determining the response of a system to an arbitrary excitation.
Laplace Transformation Method The Laplace transform of a function g t, written as G s is defined as:
4
Dividing both sides of the above equation by the mass m, we obtain: x 2zon x_ o2n x vimp d t
Figure 1 Dirac's delta function. The area under the load±time duration remains unity.
5
p where on k=m is the natural frequency of the system in radiansp per second, z c=cc is the damping ratio with cc 2 km; vimp I=m. Note that cc is called critical damping. Physically vimp is the change in the velocity of a mass m when subjected to an impulse Id t.
Z1 G s Lg t
eÿst g t dt
7
0
The function g t is defined for t < 0 and s is a complex variable. It is assumed that the function g t is such that the above integral exists. The Laplace transform is a linear operator, i.e.: Lag t br t aLg t bLr t aG s bR s
8
THEORY OF VIBRATION/Impulse Response Function 1337
Here a and b are constants, and g and r both are functions of t. Other properties of importance to us here are: Lg_ t sG s ÿ g 0 Lg t s2 G s ÿ sg 0 ÿ g_ 0
9
Knowing the Laplace transform of a function, the original function is obtained by taking the inverse Laplace transform: 1 g t L G s 2pi ÿ1
gi1 Z
G s est ds
10
gÿi1
p Here i ÿ1. Often tables of Laplace transform pairs, widely available in many books on engineering or operational mathematics, are used to obtain the Laplace transform of a function of interest. Tables are also used to obtain the inverse transform. Determining the inverse Laplace transform of a function, not given in the tables of Laplace transform pairs, may be a rather difficult task. Laplace transform of some functions that are important in the theory of vibrations are given in Table 1. The impulse response function h t can now be obtained by taking the Laplace transform of both sides of eqn [6]: s2 H s ÿ sh 0 ÿ h_ 0 2zon sH s ÿ h 0 o2n H s
1 11 m
_ Substituting h 0 0, and h 0 0, and combining all the coefficients of H s, we obtain: ÿ2 1 s 2zon s o2n H s m
12
Note that the function H s is called the transfer function for the given system and it represents the ratio of the output of the system to the input to the system in the Laplace domain. It is given as: 1 H s ÿ 2 m s 2zon s o2n
13
The impulse response function h t is obtained by taking the inverse Laplace transform of H s, i.e. h t Lÿ1 H s. We can obtain the inverse transform of H s by using partial fractions. To that end, we need to factor the denominator, i.e. represent the
denominator in the form (s ÿ s1 )(s ÿ s2 ), where s1 and s2 , the two roots of the quadratic equation, s2 2zon s w2n 0, are called the poles of the transfer function. The poles, s1 , and s2 are given as: s1;2 ÿzon on
qÿ z2 ÿ 1
14
Thus, the nature of the poles depends upon the value of the nondimensional damping parameter z. The two roots would be complex and distinct if z < 1 (underdamped case), equal and real negative if z 1 (critically damped), and distinct and real negative if z < 1 (overdamped). These three cases lead to three completely distinct type of systems as is shown in the following. Underdamped case (z < 1): The two, complex, poles are: s1;2 ÿzon ion
qÿ
1 ÿ z2 ÿzon iod
15
p where od on 1 ÿ z2 is the so-called damped natural frequency. The transfer function in terms of the partial fractions is given as: 1 1 A B ÿ m s ÿ s1 s ÿ s2 m s2 2zon s o2n
16
Table 1 Laplace transform of some functions that are important in the theory of vibrations Ld tÿt eÿts Lu t L t n
1 s
n! sn1
L eat
L sin ot
o s2 o2
L cosh ot
n 0; 1; 2; . . .
1 s ÿ a 3
L sinh ot
s s2 ÿ o2
o s 2 ÿ o2
L eat cos ot
sÿa s ÿ a2 o2
2 t Z 1 L4 g tdt5 G s s
L eat sin ot
Leat g t G sÿa
L eat cosh ot
0
L cos ot
s s2 o2
o s ÿ a2 o2
L eat sinh ot
sÿa s ÿ a2 ÿo2 o
s ÿ a2 ÿo2
1338 THEORY OF VIBRATION/Impulse Response Function
Multiplying both sides by the denominator in the above equation, we obtain: A s ÿ s2 B s ÿ s1 1
17
Substituting s s1 and s s2 , respectively, we obtain: A ÿB
1 1 i ÿ s1 ÿ s2 2iod 2od
18
p where od on z2 ÿ 1. Note that the impulse response function for the critical and overdamped case are nonoscillatory in nature. This can be seen in Figure 2B. In both Figures 2A and 2B, the plot of nondimensional displacement mon h t=I is given as a function of nondimensional time on t. Observe that as the value of the damping is increased from the critical value, the value of the nondimensional peak amplitude decreases but the time it takes for the response to die down increases.
The impulse response function, for the underdamped case, is thus given as: h t Lÿ1 H s; t > 0 ÿi 1 1 ÿ1 L ÿ 2mod s ÿ s1 s ÿ s2 Recalling Lÿ1 1= s ÿ a eat , and s1; 2 ÿzon iod , we obtain: ÿi es1 t ÿ es2 t 2mod ÿi ÿzon t iod t e e ÿ eÿiod t 2mod 1 ÿzon t e sin od t; t > 0 mod
h t
19
Here we have made use of the Euler formula eiy cos y i sin y. Note that the impulse function for the underdamped case, because of the two poles being complex conjugate for this case, is oscillatory in nature. This is shown in Figure 2A. Critically damped case (z 1): Both the roots are equal s1; 2 ÿon . The impulse response function is given as: " h t Lÿ1
1 m s on 2
#
teÿon t ; m
t>0
20
Overdamped case (z < 1): Both the roots are distinct and are given as: s1; 2 ÿzon p on z2 ÿ 1 ÿzon od . The impulse response function becomes: 1 h t L m s ÿ s1 s ÿ s2 1 ÿzon t e sin hod t ; t > 0 mod ÿ1
21
Figure 2 Nondimensional response of a single-degree-of-freedom spring-mass-damper system of mass m, and natural frequency on, to an impulse I applied at t 0 for different values of nondimensional damping z. (A) underdamped systems z < 1; (B) critically z 1 and overdamped systems z > 1. Note that the impulse response function h t is simply x t=I.
THEORY OF VIBRATION/Impulse Response Function 1339
Overdamped case
Response to Initial Conditions In the preceding section, we obtained the response of a dynamic system to a unit impulse by treating it as a special forcing function applied at t 0. Since the effect of applying a unit impulse is to impart a sudden velocity vimp 1=m at t 0, the impulse response function can also be obtained by studying the response of the system under an initial velocity _ x 0 1=m but keeping both the force f t and the initial displacement x0 equal to zero. In this section, we present the response of a second-order system under the influence of nonzero initial conditions: _ v0 . The impulse response x 0 x0 , and x 0 function can then be obtained by substituting x0 0, and v0 1=m. The response under these initial conditions, in the absence of f t, can be obtained by taking the Laplace transform of both sides of eqn [1]: s2 X s ÿ sx 0 ÿ x_ 0 2zon sX s ÿ x 0 o2n X s 0
sx0 v0 2zon x0 s2 2zon s o2n
26 For the special case of x0 0 and v0 1=m, all the three equations reduce to the respective expressions for the impulse response function given in eqns [19]± [21], respectively. Most of the systems of practical interest are underdamped systems where the nondimensional damping factor z is of the order of 0.02±0.05. The response due to the initial conditions (as was also seen for the case of the impulse response function in Figures 2A and 2B) are oscillatory in nature only for the underdamped case.
Fourier Series and Transform 22
The Laplace transform, X s, of the desired response becomes
X s
v0 zon x0 x t eÿzon t x0 cos od t sinh o t d od
23
The response x t is obtained using inverse Laplace transform by keeping in mind that three distinct cases (underdamped, critically damped, and overdamped) arise, depending on the value of the nondimensional damping parameter z. The response for the three cases, obtained using inverse Laplace transform (see Table 1), is given as Underdamped case v0 zon x0 sin od t x t eÿzon t x0 cos od t od 24
For many practical cases, we are interested in finding the response of the system to external excitation, such as buildings subjected to wind and earthquake loads, an aircraft wing subjected to gusts, an automobile on uneven pavements and so on. The impulse response function derived above can be used to obtain the response of any system to any external excitation, say f t. The response can be written as a sum of a complementary solution and a particular solution: x t xc t xp t
The complementary part xc is the solution of the governing equation with the right-hand side equal to zero, the so-called homogeneous equation, as was done in the previous section on determining the response of the solution to initial conditions. The complementary part with arbitrary constants, due to the presence of eÿzon t term, is transitory in nature and approaches zero as t increases. While determining the response due to loads that are of longer duration, the complementary solution is often ignored and emphasis is placed only on the particular solution xp t. This part of the solution, which does not have any arbitrary constants, can be obtained using the convolution (Duhamel's) integral:
Critically damped case x t eÿon t x0 v0 on x0 t
27
Zt 25
xp t f h t
f th t ÿ t dt 0
28
1340 THEORY OF VIBRATION/Impulse Response Function
Performing the above integral, at times, may be quite a difficult task, especially when the forcing function f t is not a simple function. Numerical and transform methods, Laplace and Fourier, are used in such cases. Using transform methods we can transform the above integral into an algebraic product which is easy to compute. Then, the solution in the time domain is obtained by performing the inverse transform. In previous years, the Laplace transform was the transform of choice for many engineers. However, its use is limited to only those functions whose Laplace transform is easy to obtain and also the inverse transform of the product [H(s)F(s)] is easily available. Moreover, determining the Laplace transform of a function and determining the response as function in time by performing inverse Laplace transformation does not lend itself easily to the tremendous power of a modern digital computer. It should, however, be mentioned that the availability of symbolic manipulators like Mathematica and MAPLE has considerably improved our ability to obtain inverse Laplace transform. (Caution: an effective way to compute inverse Laplace transforms of long or complex expressions in Mathematica is to first produce the partial fraction expansion in the Laplace domain and then to take the inverse of each resulting term separately. Also, it has been our experience that if inverse Laplace transforms of long expressions are taken without first performing the partial fraction expansion, one may get incorrect results.) It is because of its ability in using the power of modern computers that the Fourier transform has gained a widespread popularity for determining the system response as well as in determining the system model using experimental methods. Moreover, since the Fourier transform works in the frequency domain, it provides an additional advantage: the system parameters for many systems of practical interest can be described easier in the frequency domain than in the time domain. To understand the use of Fourier transform in vibration theory it is important to understand the use of Fourier sine and cosine series in expressing periodic, not necessarily harmonic, functions, and to understand an extension of these series for aperiodic functions, called Fourier integrals. The Fourier series for a function f t of period 2T can be written as: 1 1 X a0 X ak cos ok t bk sin ok t; 2 k1 k1 2p k ok 2T
f t
29
In general, one only needs few terms in the series to achieve a good accuracy with respect to the given function. Here ok , called harmonics, are the discrete frequencies, and the Fourier coeficients ak and bk in eqn [29] are given by: 1 ak T 1 bk T
ZT f t cos ok t dt
k 0; 1; 2 . . .
ÿT ZT
30 f t sin ok t dt
k 1; 2; 3 . . .
ÿT
Note that the coefficients ak vanish for odd functions f t ÿf ÿ t and the coefficients bk vanish for even functions f t f ÿ t. The Fourier series converges to f t for all t if f t is continuous for 0 t 2T, and converges to f t f tÿ =2 if there is a discontinuity at t. It can be shown that the coefficients obtained by eqn [30] yield the best approximation for a given number of terms used in the Fourier expansion. That is, for a given number of terms used in the Fourier series, the square of the error between the function and its Fourier representation is a minimum if the coefficients given by eqn [30] are used in eqn [29]. This implies that the accuracy of the Fourier series can only be improved by adding additional terms in the Fourier series. In vibration theory, especially in experimental structural dynamics, a complex representation is often used. It can be shown that for a periodic function f t, the Fourier representation in terms of complex exponential harmonic function, eiok t can be written as: f t
X 1 k1 f ok eiok t Do 2p kÿ1
31
where, Do, the difference between two consecutive discrete frequencies ok and ok 1 is given as: Do ok1 ÿ ok
2p 2T
32
In eqn [31], f ok is given as: f ok
ZT
f t eÿiok t dt
33
ÿT
Note that in the complex representation, the integer k varies from ÿ1 to 1. The negative frequencies are thus also included. But this representation is convenient in extending the notion of Fourier series to Fourier transform as described below.
THEORY OF VIBRATION/Impulse Response Function 1341
The Fourier series described above can be used to represent an aperiodic function, since such a function can be considered as a periodic function with the period 2T approaching 1. For this case, Do!0. Thus, instead of a discrete representation in the frequency domain, one gets a continuous variation over o. The Fourier series, in eqn [31], becomes an integral, called the Fourier integral, as given below: 1 f t 2p
Z1
f oeiot do
f o
The transfer function in the frequency domain H o is the Fourier transform of the impulse response function. For the underdamped case the result obtained by, using eqns [34] and [19], is given as: Z1 H o
In the above equation, f o is given as: f teÿiot dt
ÿ1
Functions f t and f o are considered to form a transform pair, the Fourier transform pair. Note that, in some references on Fourier transform, the factor 2p, seen here in the denominator of eqn [34] may be placed in the denominator of eqn [35] instead p and still in some other sources it may be placed as 2p in the denominator of both eqns [34] and [35]. Also note that, unlike the case here, some references define eqn [34] as the Fourier transform and eqn [35] as the inverse Fourier transform whereas in other references the names of the two equations are reversed. As is the case for the Laplace transform, the Fourier transform, denoted as F: t, is a linear operation: Fag t br t aFg t bFr t aG o bR o
36
Also, as is the case for Laplace transform, the Fourier transform of a convolution integral in the time domain becomes an algebraic product of the respective transforms in the frequency domain. That is: 2 t 3 Z F g h t F4 g th t ÿ t dt5 0
Fg tFh t G oH o
h teÿiot dt
ÿ1
35
37
where G o is the Fourier transform of g t and H o is the Fourier transform of h t as defined by eqn [35]. To obtain gh t from G oH o, one needs to apply the inverse Fourier transform given by eqn [34]. That is:
38
ÿ1
34
ÿ1
Z1
g h t Fÿ1 G oH o Z1 1 G oH oeiot do 2p
1 mod
Z1
39 eÿzon t sin od teÿiot dt
0
Here we have used the fact that h t 0 for t 0. The integral on the right-hand side of this equation can be easily obtained as: 1 ÿ m o2n 2izoon ÿ o2 "ÿ # 1 1 ÿ b2 ÿ i 2zb k ÿ1 ÿ b2 2 2zb2
H o
40
where b o=on . It is interesting to note that the transfer function in the frequency domain is the same as the steady-state (i.e. the transient part of the response containing eÿzon t term has died down) response of the system to a harmonic exponential excitation of unit amplitude and excitation frequency o (i.e. f t eiot in eqn [1]) and is called the frequency response function (FRF). The frequency of the steady-state response will be the same as that of the excitation. The complex nature of H o implies that there exists a phase difference (f) between the harmonic excitation and the resulting harmonic steadystate response. This steady-state response, xss o, is given as: xss o H oeiot jH ojei otÿf " #1=2 1 1 ei otÿf k ÿ1 ÿ b2 2 2zb2
41
Here f represents the phase difference between the harmonic excitation and the resulting steady-state response. This phase angle is given as: ÿ f tanÿ1 2zb 1 ÿ b2
1342 THEORY OF VIBRATION/Impulse Response Function
Also, since 1=k represents the response of the system under a unit static force, the term inside the square brackets in eqn [40] can be considered as the magnification of the static response if the unit load is instead applied in a harmonic manner with frequency o. The term, inside the square brackets, is called the magnification factor. In some references, the transfer function (or the FRF) given above is denoted as H io. The impulse response function can be obtained by taking the inverse Fourier transform of H o as shown in the following: h t Fÿ1 H o Z1 42 1 1 ÿ eiot do 2p m o2n 2izoon ÿ o2
z o and the semicircular arc (denoted as S) of radius r with r ! 1; as: Ic
1 2pm
1 2pm
Z1 ÿ1
Z S
ÿ
eiot do o2n 2i zoon ÿ o2
eizt ÿ dz 2 on 2i zzon ÿ z2
The second integral on the right-hand side of the above integral approaches zero as the integrand contains z2 term in the denominator. Comparing eqns [44] and [45], we obtain the impulse response function to be given as:
ÿ1
The integral on the right of the preceding equation can be performed by using Cauchy's relation, namely: I C
f z dz 2pi f z0 z ÿ z0
Here C is any simple closed path around z0 which is traversed in the anticlockwise direction and f z is analytic within C. To apply this relation to the integral in eqn [42], consider the following integral, obtained by replacing o with z, a complex variable, in eqn [42]: 1 Ic 2pm
I C
ÿ
eizt dz o2n 2i zzon ÿ z2
43
Consider the path C to be the upper half plane, i.e., a semicircle of radius r as r ! 1 (Figure 3). Realizing that the denominator in the integrand in eqn [43] when equated to zero has two complex roots: z1; 2 izon od , the integral becomes: Ic ÿ
1 2pm
I C
eizt dz 2od z ÿ i zon ÿ od
I 1 eizt dz 2od z ÿ i zon od 2pm C h i 2pi ÿei izon od t ei izon ÿod t 4pmod 1 ÿzon t e sin od t mod
h t
1 ÿzon t e sin od t mod
46
The transfer function in the frequency domain (or the FRF) H o and the impulse response function h t thus form a Fourier transform pair. In determining the impulse response function, the inverse Fourier transform was performed using an analytical approach that used the so-called Cauchy relation (or contour integral) to perform the improper integral. Such an approach is limited only to simple functions for which such integrals can be easily obtained either using Cauchy's relation or from the tables of Fourier and inverse Fourier transform (similar to the case for the Laplace transform). However, the Fourier transform, using its discretized version, can be performed numerically for all sorts of functions using the fast Fourier transform (FFT).
Fast Fourier Transform Let fn represent the value of a function at a discrete point tn nDt where Dt 2T=N, N being the total number of discrete points in the time period 2T or the time-history duration. Similarly fm represents the Fourier coefficient corresponding to the mth discrete frequency om given as om 2pm= 2T. Values of fn and fm are given as:
44
The integral in eqn [43] can also be written, by splitting the contour C in two parts: the real axis
45
fn
ÿ1 X 1N fm exp i 2pmn=N N m0
n 0; 1; 2; 3 . . . ; N ÿ 1 Nÿ1 X fn exp i 2pmn=N fm n0
m 0; 1; 2; 3 . . . ; N ÿ 1
47
THEORY OF VIBRATION/Impulse Response Function 1343
The preceding set of equations constitutes the discrete Fourier transform (DFT) pair and can easily be computed. However, if the computations are done using these equations in a direct manner, a transform will take N 2 complex multiplications. A fast method to evaluate the DFT was proposed by Cooley and Tukey. This algorithm called the fast Fourier transform (FFT), needs only O Nlog2 N complex multiplications if N 2g , where g is an integer. This results in a substantial reduction in computational time. For example, for N 1024 210 , the direct evaluation to calculate the Fourier transform will take 10242 1 048 576 complex multiplications. On the other hand, one requires only O 10 240 multiplications if FFT is employed, around 1% of the effort required otherwise. To understand the logic behind FFT, let us consider any one of the two equations in eqn [47], say the second one. It can be rewritten as: fm
Nÿ1 X
fn exp ÿi 2pmn=N
n0
Nÿ1 X
fn W mn
n0
48
m 0; 1; 2; 3 . . . ; N ÿ 1 Here W expÿi 2p=N. The FFT is based on the fact that W 0 W N W 2N . . . 1, and W N m W m . As a result, for a given value of g N 2g , the DFT can be performed in g stages. This is easily understood for N 4 i:e: g 2. For this case, eqn [48], becomes: 8 9 2 W0 > f > > = 6 0 < 0 > f1 6W 2 > 4 W 0 > f > ; :> W0 f3
W0 W1 W2 W3
W0 W2 W4 W6
3 W0 W3 7 7 W6 5 W9
8 9 f0 > > > = < > f1 f > > > ; : 2> f3
N 4 case, we calculate the DFT, once again using N complex multiplications, as: f0 g0 W 0 g1 f2 W 2 g1 g0 f1 g2 W 1 g3 f3 W 3 g3 g2 Note that the values of f are scrambled. The same approach can be easily implemented for other values of g. For example, for N 8 g 3 we will require three stages and each stage will have eight complex multiplications. The computer programs capable of performing FFT are readily available. For example, the commercially available packages Mathematica and MATLAB both contain commands for obtaining the FFT.
Nomenclature C f I v X(s) d(t) z
closed path force impulse velocity Laplace transform Dirac's delta function phase angle damping factor
See also: Critical damping; Transform methods.
49
Obviously a straightforward matrix-vector multiplication will take 42 complex multiplications. If, however, we make use of the fact that W 0 W 4 1 and W 6 W 2 , etc., we can reduce the number of multiplications to eight by performing the above multiplication in two g stages. In the first stage, using N complex multiplications, we calculate an intermediatory vector g whose components are given as: g0 f 0 W 0 f 2 g1 f 1 W 0 f 3 g2 W 2 f2 f0 g3 W 2 f3 f1 In the second (gth or last) stage, also the last stage for
Further Reading Clough RW, Penzien, J (1993) Dynamics of Structures, 2nd edn. New York: McGraw Hill. Cooley JW, Tukey JW (1965) An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation 19(90): 297±301. Craig RR (1981) Structural Dynamics, An Introduction to Computer Methods. New York: Wiley. Gerald CF, Wheatley PO (1989) Applied Numerical Analysis, 4th edn. Reading, MA: Addison Wesley. Inman DJ (1996) Engineering Vibration. Englewood Cliffs: Prentice Hall. Kelly SG (1993) Fundamentals of Mechanical Vibrations. New York: McGraw Hill. Kreyszig E (1993) Advanced Engineering Mathematics, 7th edn. New York: Wiley. Meirovitch L (1997) Principles and Techniques of Vibrations. Upper Saddle River: Prentice Hall. Rao, SS (1995) Mechanical Vibrations, 3rd edn. Reading, MA: Addison Wesley.