Paper received: 09.02.2009 Paper accepted: 00.00.200x Spatial-Mode-Shape Identification using a Continuous Wavelet Transform Martin Cesnik - Janko Slavic - Miha Boltezar* University of Ljubljana, Faculty of Mechanical Engineering, Ljubljana, Slovenia This paper presents an experimental modal analysis of a damped multi-degree-of-freedom mechanical system using a continuous wavelet transform. An approximation of the wavelet transform of the impulse response function is deduced, which serves as a basis for the extraction of the natural frequencies, the damping ratios and the mode shapes. Due to an approximation with a finite Taylor series, a computational error in the identified oscillatory amplitude occurs and is observed for the simulated system response. The presented approach of modal identification is applied to real mechanical systems, such as a steel beam and the horizontal tail of an ultralight aircraft. Using the proposed measurement methodology, it is possible to reconstruct the spatial mode shapes of any dynamic linear system with an arbitrary geometry. © 2009 Journal of Mechanical Engineering. All rights reserved. Keywords: modal parameters, continuous wavelet transform, spatial vibrations 0 INTRODUCTION One of the key steps in the development of a numerical model is its validation, where the modal analysis of a numerical model and the experimental modal analysis (EMA) of the real system are compared. With the results of the comparison it is possible to determine the adequateness of the numerical model. The most common way to undertake a modal analysis is by measuring the frequency response-functions (FRF) [1] and [2]. This is done by signal transformation from the time domain to the frequency domain using a Fourier transformation. With the use of a fast Fourier transform (FFT) algorithm, the calculation of the FRF is a simple and useful method used to perform an EMA [3] to [5]. The described method is suitable for analyzing stationary signals, since it averages the signal in time domain. For an analysis of non-stationary signals it is convenient to use transformations in the time or time-frequency domain, for instance, least-squares complex exponential (LSCE) [1], continuous wavelet transform (CWT) [6], etc. The impulse response of a mechanical system is a typical non-stationary process where the oscillation amplitude decreases exponentially with time. Generally, the modal identification, based on non-stationary signals, is usually done with time-domain methods [1], e.g. LSCE for impulse response function or Ibrahim time domain (ITD) for free decay response. The idea of applying a CWT to a dynamic system response for damping identification was first introduced by Staszewski [7]. While Staszewski used the Morlet wavelet function, Slavic et al. [8, 9, 10] developed their idea using the Gabor wavelet function and focusing it on the edge effect and the relatively short signals. The properties of the Gabor wavelet have been studied in detail by Simonovski and Boltezar [11], who also used it for fault detection in DC electro-motors [12]. Le and Argoul [13] extended the work of precendent researchers and identified all the modal parameters by applying a CWT to the theoretical impulse response of a spatial 4-degree-of-freedom (DOF) model with the Morlet, Cauchy and harmonic wavelet functions. They also explored the problem of the edge-effect and time-frequency localization. Furthermore, Lardies and Gouttebroze [14] applied a CWT to analyze real signals measured on a tower excited by the wind and determined its natural frequencies, damping ratios and one-dimensional mode shapes. They compared the results, obtained with a CWT to those obtained by an autocorrelation method. Their work was extended by Huang and Su [15], who applied a CWT to discrete equations of motion in order to identify the modal parameters of a steel frame subjected to an earthquake excitation. As a result, they obtained the frame mode shapes by a determination of the planar displacements in a horizontal plane. To continue their work, Cesnik et al. [16] tried to extract the spatial mode shapes of a mechanical system with an arbitrary geometry. This work represents a continuation and enhancement of that research. In this paper, a procedure for the identification of the modal parameters is introduced, based on a CWT of the impulse response function. Due to the low-order approximation, a computation error of the oscillatory amplitude occurs when dealing with systems with high damping, which is presented in a numerical example. The method is also tested on a real system of a steel beam The method is also tested on a real system of a steel beam where the estimated mode shapes are compared with the theoretical ones. With the identification of the spatial mode shapes of the horizontal tail of an ultralight aircraft, the feasibility of the developed method is proven. 1 CONTINUOUS WAVELET TRANSFORM The introductory part of CWT theory will be summarized according to Mallat [6] and Tchamitchian et al. [17]. The continuous wavelet transformation of the time function f(t), that satisfies f(t)e L2 (ffl) is defined as Wf(u,s) = J f(t)-vis(t)dt, (1) where v(t) represents the wavelet function, the superscript * denotes the complex conjugate, u is the translation parameter related to time and s is the scale parameter that serves as an inverse of the frequency. Fig. 1 shows the influence of the parameters u and s on the wavelet function. The function f(t) e L 2 (ffl) may be treated as a mother wavelet function when satisfying oscillatory, energy preservation and admissibility conditions [6]. A modified wavelet function Wxjt), translated in the time domain and dilated in the frequency domain also has to satisfy the energy-preservation condition, and so we obtain 1 (t — i i s Vu,s (t) = -r-V Vs v„,s H = ^ - e-^" -v(sv), (2) (3) where vu s (®) denotes the integral Fourier transformation of the wavelet function Vxjt) in the frequency domain. The impulse response of a multi-degree-of-freedom (MDOF) system is expected to be a summation of N mutually independent function pairs [7] f (t )=£ a (t y«(t), (4) where A i(t) denotes the amplitude and tpi (t) the phase (9). Therefore, the linearity property of the CWT [6] can be used N \ N W^af \(u,s) = Za Wf )(",s). (5) 1.1 Gabor Wavelet Function The Gabor mother wavelet function is defined as [6] VGabor (t __e a . (6) i ^V/4 6 6 . i™ ) By introducing the translation and dilation parameters into (6), we obtain a family of Gabor wavelet functions, defined in the time and frequency domains as [11] (t —u V i 2 2\V4 (na s ) is (a),a,n=(4 na2 s2- e e - e (7) 2 - e— mu. (8) 1.0 0.5 0.0 -0.5 -1.0 » ill 1 A _A L " \ 1 1 / v^V r i[S] Fig. 1. (—) mother VGabor (t, u = 0, s = 1) and —) modified VGabor,u,s (t,u = 4,s = 0.5) 1=1 1.2 Approximation of the CWT of the Asymptotic Function The impulse response of a viscously damped single-degree-of-freedom (SDOF) system can be written as [8] f (t)= Ae^0' • cos (®0d t + %), (9) where Z denotes the damping ratio, cp° is the phase shift, co° is the undamped angular frequency and o0d = o0^1-Z2 is the damped angular frequency. In order to separate the amplitude and the phase information of the response function it can be written in a more general form as an analytic function [17] fa (t) = Af (t) ep(t\ (10) Similarly, a wavelet function may also be given in the form of an analytic function as Wa (t) = Aw (t)• ePr(t). With regards to the relation Wf (u, s) = 1/2 • Wfa (u, s) [6] we can deduce that +, - 1 .....1........ ........ T 1 1 1 ................... ---► -CO, (Oj co2 CO, CO Fig. 2. Ridge detection using the phase method with the recursive equation (14) Due to the linear relation between frequency and scale for the Gabor wavelet '(u,s) = q/s we can conclude that the scale variable on the ridge s(u)is constant s(u)=s 0. By introducing the Gabor wavelet function into equation (12), regarding the frequency match and neglecting the error function, it follows that [6] \Wf (u,s0)| = 1Aezu (a2s2). (15) Finally, we can deduce equations for the direct identification of the damping ratio Z and the initial amplitude A d( ln |Wf (u, s0 )|) — «const., (16) z=- du a>0 A = - 2 e z®0d u Wf ( u, ^0 )| / ^V4 (4^cT2 (/®0d) ) i const. (17) 1.4 A (t) Variation Influence on the Identified Amplitude An approximation for the theoretical CWT of the impulse response is deduced by approximating an impulse response (9) as a finite Taylor series. An approximation of A(t) = const. is exact only for asymptotic signals, like the undamped response, where the phase variations are considerably faster than the amplitude variations Af '(t) Af (t ) << V'(t )|. (18) When considering signals by which the asymptotic condition is not entirely satisfied, e.g., a highly damped impulse response, an error occurs during the amplitude A identification. The presented error represents a deviation of the reconstructed amplitude Areconstr from the real amplitude in the system response Areal. It can be observed from a numerical example of a highly damped impulse response and depends on the damping ratio Z and the parameter a (Fig.3 and Fig.4); however, it is independent of the real amplitude value Areal. Consequently, the ratio between the real amplitude Areal and the reconstructed amplitude Areconstr which serves for the mode-shape reconstruction, is constant 1 real 2 real Air 2 reconstr • = const. (19) This finding is important when a reconstruction of the mode shape is made only by comparing the amplitudes of the displacements at different measurement positions (Section 2.1). Fig. 3. Reconstructed initial amplitude at variable parameter a and Areal(t=0)=1. 0 0.02 0.04 0.06 0.08 0.1 Fig. 4. Reconstructed initial amplitude at A1real(t = 0) = 1 and A2 rea} (t = 0) = 2. 1.5 Spatial Mode Shapes - Geometry In order to reconstruct the spatial mode shapes the spatial vibrations have to be measured with accelerometers at pre-defined measurement positions. When performing a modal analysis on simple geometry such as a steel beam (Section 2.1), the response accelerometers can be oriented in the same direction and no issues arise from the geometry point of view. When analyzing a system with three-dimensional geometry, the orientation of the response accelerometers in the same direction as the reference one is not always possible (Fig. 5). To overcome this problem, the local coordinate system of each accelerometer is referenced to the global coordinate system of the measured sample. One possibility of describing this linkage is with the Euler angles (Eq. 20 to 22). Rz ,y - Fig. 5. Problem of spatial vibration measurements To reconstruct a single spatial mode shape, the identified displacements need to be transformed from the local coordinate systems (as measured) into the global coordinate system with the use of the rotational matrices (20) 1 0 0 Rx,a 0 cos (a) - sin (a) 0 sin (a) cos (a) " cos (0) 0 sin (0) Ry ,P - 0 1 0 - sin (0) 0 cos (0) (21) cos (y) - sin (Y 0 sin (Y cos (y) 0 0 0 1 where a, P and y denote the Euler rotation angles about the x, y or z axis, respectively. 1.6 Identification of the Spatial Mode Shapes The procedure to perform the experiment and reconstruct the single spatial mode shape is presented in the following steps: a) Experiment 1. Definition of the excitation position. 2. Definition of the global coordinate system, n response positions with n local coordinate systems on the measured sample's geometry. 3. Selection of the reference position and the reference direction (this should be the position/direction where high oscillations of the identified modes are expected). 4. Excitation of the system with the impulse hammer, acceleration measurement on the reference position/direction and on the response position(s) in the chosen direction(s). b) Reconstruction 5. Identification of the ridge (14) and of the initial amplitudes Aref j and Aresp i,j (17) for the reference and response signals, respectively; j denotes the consecutive number of the impulse excitation and i is the number of the response position (i=1,...,n). 6. Computation of the normalized response amplitude Anorm i with respect to the reference amplitude Aref j with the ratio A norm i = A resp i,j/ Aref j. 7. Transformation of the normalized displacements Anorm i from the local to the global coordinate system (Eq. 20 to 22). 8. Graphical presentation of the spatial-mode shape. 2 EXPERIMENT For a relatively simple steel beam, the presented approach for modal identification is experimentally compared to the beam's theoretical mode shapes. Furthermore, to demonstrate an experimental identification of the spatial modal analysis, the horizontal tail of an ultralight aircraft is presented. 2.1 Experimental Validation In order to perform a specific numerical manipulation, custom programs and algorithms were developed by the authors. Custom-made programs offer a higher flexibility and provide the author with an insight into the computation process. To check the correctness of the developed programs a simple experiment with well-defined theoretical mode shapes was carried out. With algorithms that are verified on simple models, it is possible to perform a reliable modal analysis on more complex spatial structures. A simple validation experiment for method verification was performed on a free-free supported steel beam with dimensions of 5x40x1000 mm, excited by an impulse hammer at a pre-defined "excitation" position. The experiment was made with several impulse excitations; during each excitation the response accelerometer was at a different position, as described in the following. For the sake of experimental simplicity, only two accelerometers were used for the measurements. One accelerometer was denoted as the "reference" sensor and was always placed in the "reference" position. The second accelerometer, denoted as the "response" accelerometer was placed at a different "response" position during each impulse excitation. When carrying out the described experimental procedure a precaution should be taken; due to changing of position of the response accelerometer a mass matrix of the whole beam system is changing. As a result some errors of the reconstructed mode shapes and belonging natural frequencies could occur, especially when using accelerometers with considerable mass and when identifying high-order mode shapes with local nature. With application of a lightweight accelerometer with a mass of 1g and at identification of global mode shapes only negligible errors are expected at the validation experiment. Eleven measuring positions were set along the beam: position 1 was defined as the "reference" position, and the remaining positions were defined as the "response" positions, as shown in Fig. 6. 1 2 3 4 5 6 7 8 9 10 11 10 * * * i * * X X X X X J. J. „ 100 100 Fig. 6. Measuring-position placement on the analyzed beam After the whole set of measurements was carried out, we could obtain the impulse response from any pre-defined position along the beam. Fig. 7. The validation-experiment configuration set-up With the use of a CWT on the measured responses it is possible to identify the modal parameters according to the theoretical backgrounddescribed in Section 1.3. The equation (14) is used to define the time development of an observed natural frequency. Once the value of the natural frequency is known, equations (16) and (17) are used to define the damping ratio Zof the system and the initial amplitude A at the location of the measurement. Based on the ratio between the displacements at the reference and selected positions (Section 1.6), the first nine normalized mode shapes were reconstructed and compared to the theoretical mode shape by computing the Modal Assurance Criterion MAC (®X, ®A) [1], shown in Figure 8 and stated in Table 1. Table 2 shows the associated natural frequencies f and the damping ratios Z. Table 1. MAC(Qx, QA matrix Mode 1 2 3 4 5 6 7 8 9 0.9994 0.0001 0.0815 0.0000 0.0765 0.0001 0.0742 0.0001 0.0553 0.0001 0.9996 0.0000 0.0750 0.0000 0.0751 0.0000 0.0708 0.0000 0.0695 0.0001 0.9995 0.0001 0.0751 0.0000 0.0751 0.0000 0.0588 MAC 0.0002 0.0700 0.0000 0.9992 0.0000 0.0762 0.0002 0.0758 0.0000 (®x, ®A) 0.0672 0.0000 0.0709 0.0000 0.9995 0.0001 0.0781 0.0000 0.0668 0.0001 0.0701 0.0000 0.0728 0.0001 0.9992 0.0002 0.0839 0.0000 0.0683 0.0000 0.0706 0.0000 0.0757 0.0001 0.9997 0.0001 0.0885 0.0001 0.0692 0.0001 0.0716 0.0002 0.0807 0.0000 0.9992 0.0000 0.0554 0.0000 0.0582 0.0000 0.0639 0.0000 0.0851 0.0001 0.9997 Table 2. Beam's natural frequencies and damping ratios i 1. 2. 3. 4. 5. 6. 7. 8. 9. .fi ÍHzl 26.08 71.60 140.63 232.73 347.88 486.25 647.93 832.62 1040.7 ¿i .104 3.567 1.895 20.649 1.311 1.035 0.9178 1.502 0.8018 0.5266 Fig. 8. Modal assurance criterion matrix 2.2 Horizontal Tail Experiment An experiment similar to the one described in Section 2.1. was carried out on the horizontal tail of an ultralight aircraft. Since the tail geometry is described using spatial curves (Fig. 9), a spatial vibration measurement had to be performed. In order to reduce the experimental complexity and the number of nonlinearities in the system, the tail was supported free-free. Twenty-eight measuring positions were defined on the tail's surface, as shown in Fig.10. Altogether, 28 local coordinate systems were defined in such a way that the z axis was normal to the surface and the x axis followed obvious lines in the tail geometry. As a reference, measurement position 11 was chosen (Figs. 10 and 11), oriented in the local direction z. Fig. 9. Horizontal tail in isometric view 15 Fig. 10. Measuring positions' placement and the global coordinate system of the horizontal tail In each measurement two signals were sampled: - the acceleration at the reference position in the reference direction, - the acceleration at the response position in one of the defined local axes. Altogether, 84 separate measurements were made (28 measuring positions in 3 directions). Fig. 11. A horizontal tail with accelerometers placed at the reference position (left) and at the i-th response position (right) In order to transform the identified displacements from the local to the global coordinate system, special software was developed by the authors, to which an arbitrary geometry can be imported. The geometrical characteristics are imported into the software in the form of a finite-element-method (FEM) mesh that consists of volume elements. The user defines the positions and the orientations of the local coordinate systems and the corresponding identified displacements of the single mode shape. The software generates a script which can be run in ANSYS and gives a reconstructed mode shape in the form of a static problem solution. The identified natural frequencies and the damping ratios are shown in Table 2, and the experimentally reconstructed normalized mode shapes are shown in Fig.12. Table 2. Horizontal tail's natural frequencies and damping ratios_ i f [Hz] £ .io4 1. 22.8 63.17 2. 43.4 83.89 3. 75.5 141.8 a) First mode shape b) Second mode shape c) Third mode shape Fig. 12. Reconstructed normalized horizontal tail's mode shapes 3 DISCUSSION From the validation experiment, which served as an ideal, real system, it is obvious that the identified mode shapes agree with theoretical mode shapes to a high degree, although some changes in mass matrix occured due to the displacement of a response accelerometer. With the simple beam experiment the applicability of the continuous wavelet transform to determine the modal parameters is confirmed. A measurement approach with the introduction of a reference accelerometer was shown to be suitable. According to the experiment realization and the results, performed and obtained in Section 2.2, the presented approach provides a simple and feasible method for an experimental determination of the modal parameters. Although a mechanical system, such as a horizontal tail, is not an ideal system without nonlinearities, the extraction of its modal parameters with a CWT was shown to be a robust and effective method. 4 CONCLUSIONS In this paper a modal parameter identification of a dynamic system using a continuous wavelet transformation of system impulse response was presented. The computation error due to the finite Taylor series was indicated and its influence on the mode-shape reconstruction was shown. Steel beam mode shapes were reconstructed and compared to the theoretical mode shapes. The spatial mode shapes of a free-free supported horizontal tail of an ultralight aircraft were reconstructed. A comparison among the theoretical and experimental mode shapes of the beam confirmed the suitability of the proposed method for mode-shape reconstruction. The main contribution of this paper is an enhancement of known methods for EMA using a CWT for one-dimensional cases with a new method for the reconstruction of spatial mode shapes for an arbitrary geometry. In this case it is a horizontal tail. The experimental results can be used for the validation of the numerical model. 5 ACKNOWLEDGEMENT The authors would like to thank Pipistrel d.o.o. for providing the horizontal tail of the ultralight aircraft. 6 REFERENCES [1] Maia, N.M.M., Silva, J.M.M., Theoretical and experimental modal analysis, Research Studies Press, Somerset, 1997. [2] Ewins, D.J., Modal Testing: Theory and Practice, Research Studies Press, Lechtworth, 1986. [3] Cermelj, P., Boltezar, M., An indirect approach to investigating the dynamics of a structure containing ball bearings, Journal of Sound and Vibration, 2004, vol. 276, no. 1/2, p. 401-417. [4] Cermelj, P., Boltezar, M., Modelling localised nonlinearities using the harmonic nonlinear super model, Journal of Sound and Vibration, 2006, vol. 298, no. 4/5, p. 10991112. [5] Boltezar, M., Cermelj, P., Dynamics of complex structures - valid modeling of industrial products, In: Korelc, J., Zupan, D. (Eds.): Kuhljevi dnevi 2007, Snovik, p. 1724, 2007. [6] Mallat, S., A wavelet tour of signal processing, Academic Press, 1999, second edition. [7] Staszewski, W.J., Identification of damping in mdof systems using time-scale decomposition, Journal of Sound and Vibration, 1997, vol. 203, p. 283-305. [8] Slavic, J., Simonovski, I., Boltezar, M., Damping identification using continuous wavelet transform: application to real data, Journal of Sound and Vibration, 2003, vol. 262, p. 291-307. [9] Slavic, J., Boltezar, M., Enhanced identification of damping using continuous wavelet transform, Journal of Mechanical Engineering, 2002, vol. 48, no. 11, p. 621631. [10] Boltezar, M., Slavic, J., Enhancements to the continuous wavelet transform for damping identifications on short signals, Mechanical systems and signal processing, 2004, vol. 18, no. 5, p. 1065-1076. [11] Simonovski, I., Boltezar, M., The norms and variances of the gabor, morlet and general harmonic wavelet functions, Journal of Sound and Vibration, 2003, vol. 264, p. 545557. [12] Boltezar, M., Simonovski, I., Furlan, M., Fault detection in DC electro motors using the continuous wavelet transform, Meccanica, 2003, vol. 38, no. 2, p. 251-264. [13] Le, T.P., Argoul, P., Continuous wavelet transform for modal identification using free decay response, Journal of Sound and Vibration, 2004, vol. 277, p. 73-100. [14] Lardies, J., Gouttebroze, S., Identification of modal parameters using free decay response, International Journal of Mechanical Sciences, 2002, vol. 44, p. 2263-2283. [15] Huang, C.S., Su, W.C., Identification of modal parameters of a time invariant linear system by continuous wavelet transform, Mechanical Systems and Signal Processing, 2006, vol. 21, p. 1642-1664. [16] Cesnik, M., Slavic, J., Boltezar, M., Modal parameter identification using continuous wavelet transform, In: Boltezar, M., Slavic, J. (Eds.): Kuhljevi dnevi 2008, Cerklje na Gorenjskem, p. 41-48, 2008. [17] Tchamitcian, P., Torresani, B., Asymptotic wavelet and gabor analysis: Extraction of instantaneous frequencies, IEEE transactions on information technology, 1992, vol. 38, p. 644-664.