An FFT-based spectral analysis method for linear discrete dynamic systems with non-proportional damping

This paper proposes a fast Fourier transforms (FFT)-based spectral analysis method for the dynamic analysis of linear discrete dynamic systems which have non-proportional viscous damping and are subjected to non-zero initial conditions. To evaluate the proposed FFT-based spectral analysis method, the forced vibration of a three degree-of-freedom (DOF) system is considered as an illustrative problem. The accuracy of the proposed FFT-based spectral analysis method is evaluated by comparing the forced vibration responses obtained by the present FFT-based spectral analysis method with those obtained by using the well-known Runge-Kutta method and modal analysis method.


Introduction
A large number of analytical and numerical solution methods have been well developed for the linear discrete dynamic systems.By virtue of impressive progress in computer technologies during the last three decades, there have been developed diverse computer-based numerical methods to obtain satisfactory approximate solutions, especially for the coupled large DOF systems.They may include the direct integration methods, the modal analysis methods, the discrete-time system methods, and the spectral analysis methods in which the FFT technique is utilized.The first three are the time-domain methods [1,2], while the FFT-based spectral element method (SAM) is a frequency-domain method [3][4][5][6].
In the FFT-based SAM, the dependent variables of a set of ordinary differential equations are all transformed into the frequency-domain by using the discrete Fourier transforms (DFT) to transform the ordinary differential equations into a set of algebraic equations with frequency as the parameter.The algebraic equations are then solved for the Fourier (or spectral) components of dependent variables at each discrete frequency.As the final step, the time-domain responses are reconstructed from the Fourier components by using the inverse discrete Fourier transforms (IDFT).In practice, the FFT is used to carry out the DFT or IDFT.As the FFT is a remarkably efficient computer algorithm, it can offer an enormous reduction in computer time and also can increase solution accuracy [4,7].
The FFT-based SAM has been known to be very useful especially in the following situations [4,6,7]: (1) when the modern data acquisition systems are used, as in most experimental measurements, to store digitized data through the analogue-to-digital converters, (2) when the excitation forces are so complicated that one has to use numerical integration to obtain the dynamic responses by using the excitation values at a discrete set of instants, (3) when it is significantly easier to derive the constitutive equation of a material in the frequency-domain rather than in the time-domain, and (4) when the frequency-dependent spectral element (or dynamic stiffness matrix) model is used as a structure model.So it may be worthy to mention that the development of an FFT-based SAM is not necessary to compete with the exiting time-domain methods such as Runge-Kutta method, for instance.The FFT-based SAM will be especially useful in the aforementioned situations in which the frequency-domain analysis is more desirable.
The FFT-based SAM has been well applied to the prediction of the steady-state responses of dynamic systems [3][4][5][6]8].However the application of the FFT-based SAM to the transient responses has been limited to the dynamic systems with all zero initial conditions.As an effort to deal with dynamic systems with nonzero initial conditions, Veletsos and Ventura [9,10] introduced a DFT-based approach to calculate the transient responses of a linear 1-DOF system.Their procedure involves the superposition of a corrective, free vibration solution which effectively transforms the steady-state response to the desired transient response.Later Mansur et al. [11,12] used the pseudo-force concept to take into account non-zero initial conditions in the DFT-based frequency-domain analysis of continuous media discretized by FEM.Recently Lee et al. [13] developed an FFT-based spectral analysis method for the linear discrete dynamic systems with nonzero initial conditions.As the proportional viscous damping was considered in reference [13], the modal decomposition approach could be used to develop a modal-coordinates -based SAM.However the modal-coordinates-based SAM proposed in reference [13] can not be directly applied to the dynamic systems with non-proportional viscous damping.
Thus, the purpose of this paper is to develop an FFT-based SAM for the linear discrete dynamic systems with non-proportional viscous damping.The present FFT-based SAM is unique because it does not use the superposition of corrective, free vibration solution to match the initial conditions as in references [9,10], or the pseudo-force concept to take into account the non-zero initial conditions as in references [11,12], or the modal decomposition method as in reference [13].To evaluate the FFT-based SAM developed in this paper, the forced vibration of a viscously damped three-DOF vibration system is considered as an illustrative example.

Theory of discrete Fourier transforms
Because the theory of discrete Fourier transforms (DFT) is one of the key mathematical tools used to develop the present FFT-based SAM, a brief review on the DFT theory will be given in the following.A periodic function of time x(t), with period T , can be always expressed as a Fourier series of the form where i = √ −1, ω n = n(2π /T) = nω 1 are the discrete frequencies, and X n are constant Fourier components given by Equations ( 1) and ( 2) are the continuous Fourier transforms pair for a periodic function.Although x(t) is a continuous function of time, it is often the case that only sampled values of the function are available, in the form of a discrete time series {x(t r )}.If N is the number of samples, all equally spaced with a time interval equal to ∆ = T /N , the discrete time series are given by x r = x(t r ), where t r = r∆ and r = 0, 1, 2, . .., N − 1.The integral in Eq. ( 2) may be replaced approximately by the summation which is the discrete Fourier transforms (DFT) of the discrete time series {x r }.Any typical value x r of the series {x r } can be given by the inverse formula which is the inverse discrete Fourier transforms (IDFT).Thus, Eqs (3) and ( 4) represent the DFT-IFFT pair.Even though Eq. ( 3) is an approximation of Eq. ( 2), it is important to note that it allows all discrete time series {x r } to be regained exactly [6,7].The Fourier components X n in Eq. ( 4) are arranged as X N −n = X * n (n = 0, 1, 2, . .., N/2), where X * n represents the complex conjugate of X n .Note that X N/2 corresponds to the highest frequency ω N/2 = (N/2) ω 1 , which is called the Nyquist frequency.
The fast Fourier transforms (FFT) is an ingenious highly efficient computer algorithm developed to perform the numerical operations required for a DFT or IDFT, reducing the computing time drastically by the order of N/log 2 N. It should be pointed out that while the FFT-based spectral analysis uses a computer, it is not a numerical method in the usual sense, because the analytical descriptions of Eqs ( 3) and ( 4) are still retained.Further details of DFT and FFT can be found in reference [7].

Development of spectral analysis method
The forced vibration of a viscously damped m-DOFs dynamic system can be represented by the matrix equation of motion and the initial conditions where u(t) is the nodal DOFs vector and {f (t)} is the nodal forces vector.The matrices . ., m) are the mass matrix, the non-proportional viscous damping matrix, and the stiffness matrix, respectively.The total dynamic response of the system can be obtained by the sum of the steady-state response, {u p (t)}, and the transient response, {u h (t)}, as follows: The approach to compute the steady-state response {u p (t)} is basically the same as that introduced in the previous work [13] and it will be briefly summarized in the following, while the approach to compute the transient response {u h (t)} will be newly developed in this paper.

Steady-state responses
Assume that the nodal force vector {f (t)} and the steady-state response vector {u p (t)} can be represented in the spectral forms as Applying Eq. ( 8) into Eq.( 5) yields where [D (ω)] is the dynamic stiffness matrix defined by Once the Fourier components {P n } are computed from Eq. ( 9) for a given nodal force vector {f (t)}, the steady-state response vector in the time-domain can be reconstructed by using the IFFT algorithm as From Eq. (8b), the time derivative of {u p (t)} can be obtained as where

Transient responses
The transient response vector {u h (t)} must satisfy the homogeneous matrix equation of motion which can be reduced from Eq. ( 5) by enforcing {f(t)} = 0 as follows: Because [C] is the non-proportional damping matrix, Eq. ( 14) cannot be decoupled by using the modal decomposition analysis.Thus assume the general solution of Eq. ( 14) in the form or Substituting Eq. ( 15) into Eq.( 14) gives For the existence of non-trivial solution {A}, it follows that Equation ( 19) yields a 2m-degree algebraic equation with λ as an unknown.In general the roots (eigenvalues) of the algebraic equation are of complex form.As discussed in reference [14], the complex roots will appear in the complex conjugate pairs for the underdamped system, because all the coefficients of the algebraic equation are real.Thus, without loss of generality, the eigenvalues of Eq. (17) or Eq. ( 18) can be written in a complex conjugate form as where ω j represents the natural frequency and ξ j represents the rate of exponential decay of the jth damped vibration mode.
By substituting an eigenvalue λ = λ j into Eq.( 17) or Eq. ( 18), the corresponding eigenvector {A} can be computed.The ratio between the components a k of the jth eigenvector are given by [14,15].
where C jk is the co-factor of the jth row of the determinant in Eq. ( 19) for a particular λ j , and z j is an arbitrary complex number.From Eq. ( 16), the kth component of the jth transient vibration mode corresponding to the jth eigenvalue λ j is given by Since there are 2m eigenvalues λ j (j = 1, 2, . .., 2m), the kth component of the total transient response vector {u h (t)} can be compounded from Eq. ( 23) as where bj are constants.Substitute Eq. ( 23) into Eq.( 24) and use the fact that the 2m eigenvalues are in the m-pairs of complex conjugate to obtain where B j and B * j (j = 1, 2, . . ., m) are constants to be determined by initial conditions.The time derivative of Eq. ( 25) is given by The functions u hk (t) and uhk (t) can be represented into the spectral forms as H kn e iωntr (27) By the DFT theory, the Fourier components H kn and H kn can be expressed as Substituting Eqs (25) and (26) into Eq.( 28) gives Substitute Eq. (29) into Eq.( 27) and collect u hk (t) and uhk (t) (k = 1, 2, . . ., m) to form vectors as Hn e iωntr where Now consider the initial conditions.The total dynamic response given by Eq. ( 7) should satisfy the initial conditions given by Eq. ( 6).Thus, substituting Eq. ( 7) into Eq.( 6) gives Applying Eqs (8b), ( 12) and (31) into Eq.(34) gives By using Eq.(29), Eq. ( 33) can be rewritten as Substituting Eq. (36) into Eq.(35) gives The constants vectors {B} or {B * } can be solved from Eq. (38) as Once the constants vector {B} is computed from Eq. ( 41) by using given initial conditions, the Fourier components {H n } are computed first from Eq. (36).Then one can readily compute the transient responses {u h (t)} by using the IFFT algorithm as follows:

Numerical examples and discussion
To evaluate the present (FFT-based) SAM, two viscously damped three-DOFs dynamic systems are considered as example problems.The first one is the case with the proportional viscous damping, for which the modal analysis method can be readily applied to get the analytical solutions, and the second one is the case with the non-proportional viscous damping.
(a) Case 1: proportional viscous damping [13] where m = 10 kg, c = 40 N • s/m, k = 1 kN/m, and 2)] N , where s(t) represents the unit step function.The initial conditions are given by Figure 1 compares the dynamic responses (u 1 , u 2 and u 3 ) obtained by three different solution methods for the dynamic system with proportional viscous damping: the present SAM, the modal analysis method and the (4th order) Runge-Kutta method.It is quite straightforward to apply the modal analysis method to obtain exact analytical solutions for the linear dynamic systems with proportional viscous damping.Thus the dynamic responses exactly obtained by the modal analysis method are used as the exact reference solutions to evaluate the other solutions obtained by the present SAM.The DFT period T = 4.8 seconds and the number of samples N = 2 11 are used for the present SAM to obtain the dynamic responses within 0.1% time averaged error with respect to the exact reference solutions, whereas the time increment ∆t = 0.00234 seconds is used for Runge-Kutta method.Figure 1 shows that the present SAM provides accurate solutions which are very close to the exact reference solutions and also to the numerical solutions obtained by Runge-Kutta method.
Figure 2 compares the dynamic responses obtained by the present SAM and the Runge-Kutta method for the dynamic system with non-proportional viscous damping.The same values of T , N , and ∆t as used for Fig. 1 are consistently used for Fig. 2. One may find from Fig. 2 that the present SAM certainly provides the dynamic responses which are very close to those obtained by Runge-Kutta method.
Figure 3 shows the convergence of the dynamic responses of the dynamic system with non-proportional viscous damping obtained by the present SAM as the sampling number N is increased.As expected, Fig. 3 shows that more accurate solutions can be achieved by increasing N for a fixed time window T = 4.8 seconds.When the sampling number N = 2 12 is used, the time average of the absolute difference between the dynamic responses obtained by SAM and Runge-Kutta method is found to be within 0.1% of the time average of the absolute dynamic response obtained by Runge-Kutta method.
Similarly Fig. 4 shows the convergence of the dynamic responses obtained by the present SAM as the time window size T is reduced for fixed sampling number N = 2 11 .As also expected from DFT theory, Fig. 4 shows that more accurate solutions can be achieved by decreasing N for a fixed sampling number.
Finally, it is worthwhile to confirm from Fig. 1 through Fig. 4 that the present SAM certainly captures all non-zero initial conditions accurately in the dynamic responses, which might be one of motivations of the present paper.

Conclusions
In this paper, an (FFT-based) SAM is developed to obtain the dynamic responses of a linear discrete dynamic system with non-proportional damping, subjected to non-zero initial conditions.In due course, the SAM is evaluated by comparing the dynamic responses of an example three DOFs vibration system with non-proportional damping obtained by using the present SAM with the exact analytical solutions obtained by the modal analysis method and also with the numerical solutions obtained by the Runge-Kutta method.It is shown that sufficiently accurate solutions can be achieved by the present SAM, when compared with the solutions by the modal analysis method and Runge-Kutta method, by properly choosing the sampling number for a given time window in the FFT analysis process.But, it may be worthy to mention that the development of the present SAM is not necessary to compete with the exiting well-known time-domain solution methods such as Runge-Kutta method, for instance.Instead, it has been developed to be useful in some special situations in which such frequency-domain method can be more easily applied: the special situations may include when the excitation forces are measured as the digitized data and when the structural properties such as the stiffness and damping coefficients are measured or provided as the frequency-dependent properties.

Fig. 1 .
Fig.1.Comparison of the dynamic responses of the dynamic system with proportional viscous damping obtained by the present SAM, modal analysis method (exact) and Runge-Kutta method.

Fig. 2 .
Fig. 2. Comparison of the dynamic responses of the dynamic system with non-proportional viscous damping obtained by the present SAM and Runge-Kutta method.

Fig. 3 .
Fig. 3. Convergence of the dynamic responses of the dynamic system with non-proportional viscous damping obtained by the present SAM as the sampling number N is increased for fixed time window T = 4.8 seconds.

Fig. 4 .
Fig. 4. Convergence of the dynamic responses of the dynamic system with non-proportional viscous damping obtained by the present SAM as the time window size T is reduced for fixed sampling number N = 2 11 .