Wavelet-based modal analysis for time-variant systems

Mechanical Systems and Signal Processing 50-51 (2015) 323–337

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Wavelet-based modal analysis for time-variant systems K. Dziedziech, W.J. Staszewski n, T. Uhl Department of Robotics and Mechatronics, AGH University of Science and Technology, Al. Mickiewicza 30, 30-059 Kraków, Poland

a r t i c l e i n f o

abstract

Article history: Received 2 August 2013 Received in revised form 28 April 2014 Accepted 2 May 2014 Available online 7 June 2014

The paper presents algorithms for modal identification of time-variant systems. These algorithms utilise the wavelet-based Frequency Response Function, and lead to estimation of all three modal parameters, i.e. natural frequencies, damping and mode shapes. The method presented utilises random impact excitation and signal post-processing based on the crazy climbers Algorithm. The method is validated using simulated and experimental data from vibration time-variant systems. The results show that the method captures correctly the dynamics of the analysed systems, leading to correct modal parameter identification. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Vibration time-variant systems Wavelet analysis Modal analysis Wavelet-based frequency response function System identification

1. Introduction Modal analysis is an important and useful tool used for various engineering applications, such as design or dynamic testing. Parameters gathered from modal analysis – i.e. natural frequencies, damping ratios and mode shapes – are important not only to achieve desirable properties and performance but also to prevent undesirable characteristics (e.g. excessive noise and vibrations) and catastrophic failures. These parameters can be obtained from the Frequency Response Function through numerical simulations and experimental analysis. For a simple case of linear time-invariant systems, methods for estimation of FRFs are well known and established. However, these classical methods often fail to deal properly with systems that are non-linear and/or time-variant. Various approaches can be used for identification of time-variant systems. Altogether these approaches can be divided into two major groups. The first group considers parametric methods that are based on parametric models. A general form of FRF is then assumed to be a ratio of two polynomials with unknown coefficients. The denominator of this ratio contains information about natural frequencies and damping, whereas the numerator gives information about the geometry of the analysed system, i.e. deflection amplitudes and phases. The Auto-Regressive (AR), Auto-Regressive Moving-Average (ARMA), Auto-Regressive Moving-Average with eXogenous inputs (ARMAX) models and many other similar approaches can be used for linear time-invariant systems. Time-variant systems are often analysed using similar – but somehow modified – approaches with coefficients that depend on time. There exist a number of time-dependent parametric methods, starting from simple models such as the Time-dependent AR (TAR) and Time-dependent ARMA (TARMA) to more complex and computationally demanding models such as the Recursive Maximum Likelihood estimated TARMA (RML-TARMA), Smoothness Priors TARMA (SP-TARMA), Functional Series TARMA (FS-TARMA) and other variations of these models, are discussed in [1,2].

n

Corresponding author. Tel.: +48 12 617 3505; fax: +48 12 617 3505. E-mail addresses: [email protected] (K. Dziedziech), [email protected] (W.J. Staszewski), [email protected] (T. Uhl).

http://dx.doi.org/10.1016/j.ymssp.2014.05.003 0888-3270/& 2014 Elsevier Ltd. All rights reserved.

324

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

The second group of methods for time-variant systems considers non-parametric methods, in which no a priori information about the analysed systems is required. Classical input–output methods – that are used for time-variant systems – consider output and input spectra obtained via the Fourier transform and leading to FRFs. Various extensions of these classical methods are based on time–frequency analysis, as discussed in [1]. The most commonly known approach uses the well-known Short Time Fourier Transform (STFT), which is both simple and relatively powerful. The STFT-based concept breaks the entire signal into small time intervals and then performs the Fourier Transform for all short intervals. This operation enables localisation of short-time events in both time and frequency domains. Another widely used timefrequency approach is based on the Wigner–Ville distribution. These two methods are considered to be members of the so-called Cohen's class of distributions [3,4]. There are many other members of this class, including for example the Rihaczek [5] and ChoiWilliams [6] distributions. The continuous wavelet transform [7,8] is a time-scale method that is qualitatively different from the other members of the Cohen's class. This approach employs different types of shifted and dilated wavelet functions to decompose signals. The analysis in the combined time–frequency domain is the first step of the entire process of modal identification. Modal models have to be extracted from the time–frequency domain. Methods for damping estimation are presented in [9–11]. Methods for the nonlinear system analysis are described in [12]. Methods for estimation of instantaneous frequency and rotational velocity from vibrations are presented in [13]. Recently online identification approaches based on adaptive wavelets were developed, as demonstrate in [14,15]. A different approach based on wavelet has been proposed in [16], where the transfer function is used for system identification. A good overview of a various wavelet-based approaches can be found in [17,18]. Other than time–frequency and time-scale approaches investigated are: pseudo-modal parameters based on eigenvalues and discrete-time state transition matrix [19], time-dependent mode shapes based on eigenvectors and Taylor expansion [20], point evolution approach based on complex output signals attached to input signals [21], partial characteristics computed for selected time periods [22], underspread theory [23] or autoregressive moving average models [24]. All above approaches are adaptive transfer functions, mainly for periodic time-variances. It is important to note that various output-only methods are also used for the analysis of time-variant systems. This includes: peak picking, frequency domain decomposition, random decrement technique, Ibrahim time-domain method, least square complex exponential method, stochastic subspace identification algorithm, various modified ARMA approaches and [25], as reviewed in [17]. A few attempts have been made to introduce time-variant FRFs. This includes the well-known concepts based on evolutionary spectra [26,27], frozen spectra [1,28–30], identification algorithms based Littlewood-Paley wavelets [14], time– frequency based methodology [31] and the wavelet-based FRF [32,33]. The latter approach has been extended in [34] to provide the theoretical background, physical interpretation based on the wavelet convolution and numerical implementation. The focus of the work was on the amplitude of the wavelet-based FRF. Only one modal parameter, i.e. natural frequency, was investigated. The need for time–frequency averaging is the major drawback of the work presented. The effect of input excitation on FRF amplitude has been investigated in [35]. The paper builds upon recent developments reported in [31,34,35]. The major objective of the work presented is to propose algorithms for modal identification of time-variant systems. These algorithms utilise the wavelet-based FRF and lead to estimation of all three modal parameters, i.e. natural frequencies, damping and mode shapes. In order to achieve the objective a number of new post-processing algorithms for wavelet-based modal analysis are proposed. A new procedure for the phase extraction from the wavelet-based FRF and an improved algorithm for natural frequency estimation are proposed. The latter is based on the so-called crazy climbers Algorithm. The work described in this paper also demonstrates how the problematic, combined time–frequency averaging – used in [34] – can be avoided with the application of the newly proposed random impact excitation. The structure of the paper is as follows: The wavelet-based FRF is briefly described in Section 2 using time-dependent auto- and cross-power spectra. Numerical implementation of the method is presented in Section 3. This part of the paper also describes the new algorithms for excitation and modal analysis post-processing algorithms, where the wavelet ridge extraction algorithm based on crazy climbers, the FRF phase extraction algorithm and algorithms for modal parameter estimation are presented. The proposed methodology for time-variant analysis is illustrated in Section 4 using simulated examples and validated in Section 5 using experimental data. The latter utilises a simple test rig with a varying mass. Finally, the paper is concluded and the future work is proposed in Section 6.

2. Wavelet-based frequency response function This section briefly introduces the wavelet-based FRF. The time-dependent auto- and cross-power spectra are firstly introduced to provide the definition. Physical interpretation is then provided.

2.1. Definition For the input x(t) and output y(t) signals, the relevant wavelet transforms can be defined as   Z 1 1 t b dt xðtÞψ Xða; bÞ ¼ pffiffiffiffiffiffi a jaj  1

ð1Þ

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

1 Yða; bÞ ¼ pffiffiffiffiffiffi jaj

Z



1

yðtÞψ 1

 t b dt a

325

ð2Þ

where b is translation indicating locality in the time domain, a is a scale parameter providing locality in the frequency domain and ψðtÞ is a continuous function – in both time and frequency domains – called the mother wavelet. Following the above definitions the relevant input and output wavelet auto-power spectra can be defined as N

^ XX ða; bÞ ¼ 1 ∑ X i ða; bÞX n ða; bÞ G i Ni¼1 N

^ YY ða; bÞ ¼ 1 ∑ Y i ða; bÞY i n ða; bÞ G Ni¼1

ð3Þ

ð4Þ

where the superscript “n” indicates a complex conjugate and N is the number of averages These two functions show how the power in a signal is distributed over the combined time–frequency plane. By analogy, the relevant wavelet cross-power spectra can be defined as N

^ XY ða; bÞ ¼ G ^ YX ða; bÞ ¼ 1 ∑ X i ða; bÞY i n ða; bÞ G Ni¼1

ð5Þ

Eqs. (3)–(5) can be used to define the wavelet-based FRF. Different classical estimators – such as H1, H2 or Hv – can be used in practice. The work presented in this paper uses the H1 estimator for the wavelet-based FRF. This estimator can be defined as H 1 ða; bÞ ¼

^ YX ða; bÞ G ða; bÞ ^ XX G

ð6Þ

By analogy to the classical modal analysis, the wavelet-based coherence can be defined as γ¼

^ XY ða; bÞj2 jG ^ XX ða; bÞ ða; bÞG ^ YY G

ð7Þ

Some comments are needed here with respect to the averaged power spectra. Numerical implementation of Eqs. (6) and (7) is very difficult in practice. Time–frequency domain averaging is not easy when time-variant systems are analysed. Numerical simulations allow one to obtain a number of different biased excitation and response signals that can be averaged. However, it is very difficult in practice to capture the same time-variant behaviour of the analysed system when experimental measurements are taken. One possible solution to this problem is to correlate the relevant input and output signals, as demonstrated in [34]. However, this is not easy in practice. Another possible solution is to use the biased (i.e. not averaged) H1 estimator for N ¼1 and to apply some additional signal post-processing to improve the results. The latter approach is used in the current paper and the extra post-processing is explained in Section 3. It is clear that Eq. (7) cannot be utilised in practice when this approach is used. It is clear that the H1 estimator – given by Eq. (6) – results in a complex function. The absolute value and argument of the H1 gives the amplitude and phase of the wavelet-based FRF as H 1a ða; bÞ ¼ absðH 1 ða; bÞÞ

ð8Þ

H 1p ða; bÞ ¼ argðH 1 ða; bÞÞ

ð9Þ

respectively. The above equations are used in the paper for the wavelet-based FRF calculations. 2.2. Physical interpretation Question now arises what the ratio in Eq. (6) represent. The convolution theorem is important to answer this question. For a given mechanical system, dynamic properties can be described in the time domain using the impulse response function h(t). Signal response can be then described as Z þ1 yðtÞ ¼ xðtÞnhðtÞ ¼ xðτÞhðt τÞdτ ð10Þ 1

where “n” is the convolution operator and τ is a convolution variable. It is well known that the continuous wavelet transform can be expressed as a convolution of analysed signal and the scaled mother wavelet, i.e. Z 1 xðτÞψ a ðb  τÞdτ ð11Þ Xða; bÞ ¼ xðbÞnψ a ðbÞ ¼ 1

where ψ a ðbÞ indicates the mother wavelet for a given scale a. Combining Eqs. (10) and (11) yields yðtÞnψ a ðtÞ ¼ ðxðtÞnhðtÞÞnψ a ðtÞ

ð12Þ

326

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

When the convolution theorem is used it is easy to show – but only for a given (i.e. fixed) scale function – that xðtÞnyðtÞ ¼ FT  1 fFTfxðtÞgFTfyðtÞgg

ð13Þ

1

indicated the Fourier and inverse Fourier transforms, respectively. As a consequence, when Eq. (13) is where FT and FT applied to Eq. (12) the following can be obtained as: FT  1 fFTfyðtÞgFTfψ a ðtÞgg ¼ FT  1 fFTfFT  1 fFTfxðtÞgFTfhðtÞgggFTfψ a ðtÞgg

ð14Þ

and in a shorter form as FT  1 fFTfyðtÞgFTfψ a ðtÞgg ¼ FT  1 fFTfxðtÞgFTfhðtÞgFTfψ a ðtÞgg

ð15Þ

Applying the Fourier transform to both sides of the above equation yields FTfyðtÞgFTfψ a ðtÞg ¼ FTfxðtÞgFTfhðtÞgFTfψ a ðtÞg

ð16Þ

As a result, the Fourier transform of the impulse response function can be written as FTfhðtÞg ¼

FTfyðtÞgFTfψ a ðtÞg FTfxðtÞgFTfψ a ðtÞg

ð17Þ

given finally – again for a fixed scale value – the frequency domain representation as FTfhðtÞg ¼

FTfyðtÞg FTfxðtÞg

and the time domain representation as   FTfyðtÞg hðtÞ ¼ FT  1 FTfxðtÞg

ð18Þ

ð19Þ

Eqs. (18) and (19) show that the classical theory holds for time-variant systems. However, the relevant analysis is fixed to only one scale (or frequency) value. When this assumption is used then the wavelet-based FRF indeed represents dynamics of the analysed mechanical system. More generic physical interpretation for all scales is not a trivial task and requires the convolution for time-scale representations. The concept of the generalised translation operator and convolution for wavelets – proposed in [36] – has been used in [34] to interpret the wavelet-based FRF. However, it is important to note that this concept is useful only for physical interpretation since analytical forms of the wavelet-based FRF for even simple systems are not easy to obtain in practice, as explained in [34]. 3. Numerical implementation This section briefly describes numerical implementation of the wavelet-based FRF used in the current investigations. The problem of wavelet function selection is briefly discussed. Input excitation is defined. Finally, signal processing – used to improve identification results from the biased H1 estimate – is described. 3.1. Wavelet function It is well known that when time–frequency analysis is used proper selection of mother wavelet function is a very important, as discussed in [37]. Wavelet desirable features are: compact support in both – i.e. time and frequency – domains and complex nature of wavelets. The need for compact support in the time and frequency domains is required by timevariant nature of analysed systems. Clearly, the ability to capture this time-variant behaviour is essential when identification is performed. Complex – rather than real – wavelet functions are needed to avoid amplitude zero crossings in the time and frequency domains. This is important when additional signal processing is applied to improve the results, as described in Section 3.3. A number of different wavelet functions – that fulfil these requirements – can be used in practice. The work presented in this paper utilises the complex Morlet wavelet function defined as ψðtÞ ¼ e 

jtj2 2

ejω0 t

ð20Þ

This function has been selected following numerous previous successful applications related to the analysis of time-variant and nonlinear dynamic systems [9,12,34]. Fig. 1 illustrates the complex Morlet wavelet function in the time and frequency domain. 3.2. Input excitation When the input–output analysis of time-variant systems is considered it is important to use broadband excitation that provides energy to the system for a relatively long period of time. This problem has been discussed in [35]. Various types of broadband excitations can be used in the classical input–output dynamic analysis. This includes: chirp, white noise and impact excitations. Chirp excitation is a commonly used broadband excitation. A pre-defined range of frequencies can be

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

327

Fig. 1. Complex Morlet wavelet function presented in the time (a) and frequency (b) domains.

easily selected to provide a broadband nature. Analysis of this excitation in the time domain can reveal some information related to time and frequency localisation. The former is very straightforward and does not need any explanation. The latter can be achieved by analysing zero-crossings and calculating periods of oscillations. It is important to note that whenever small parts of signal are considered chirp excitation is not truly broadband. This is a major problem associated with the chirp excitation when used for the input–output analysis of time-variant systems, as illustrated in [35]. White noise excitation is commonly used in modal analysis. However, white noise is considered to be broadband only when it is long enough to be ergodic. The application of white noise excitation for the input–output analysis of time-variant systems requires relatively long signals and time–frequency averaging [34,35]. Impact excitation is one of the most commonly used excitations in engineering applications. This type of excitation provides perfect localisation in the time domain and a truly broadband nature in the frequency domain. The major drawback of this excitation – when used for the analysis of time-variant systems – is the fact that the system is excited only once for a short period of time. When heavy damping is involved, the input– output analysis is limited only to a very short period of time. This is not sufficient for system identification, as demonstrated in [35]. Random impact excitation has been recently proposed for the analysis of time-variant systems in [35]. This excitation is a series of impacts spaced randomly in the time domain. Such excitation is broadband at given time instances and provides continuous inflow of energy to the system. This is advantageous when time-variant systems are analysed since continuous inflow of energy enables one to observe system dynamic properties and their variations with time. Random impact excitation – illustrated in the time, frequency and combined time–frequency domain in Fig. 2 – has been used for the input–output analysis of time-variant systems in the current investigations.

3.3. Signal post-processing for system identification Additional signal post-processing has been applied to improve the results obtained from the biased wavelet-based FRF H1 estimator, calculated for the random impact excitation. This is required for reliable identification of modal parameters. The concept of wavelet ridges – well established for time–frequency and time-scale transformations – has been used to achieve this task. This concept – previously used for the application of the wavelet-based FRF [31,34] will be also utilised in the current studies. One should note that when ridges are used, selection of proper excitation is important. When repeated impacts or white noise excitation are applied, frequency values are not distorted in theory. However, this is not the case with chirp excitation. Chirp excitation has a very high amplitude attached to a single frequency. This will directly affect ridges of the wavelet-based auto-power spectrum, calculated from the output signal. Calculations of the wavelet-based FRF need the division by the wavelet-based auto-power spectrum of the input signal. Therefore possible frequency distortions are reduced. Nevertheless, the problem cannot be fully avoided and needs further investigations. However, in contrast to previous investigations a new – and more efficient – algorithm for the ridge extraction is proposed in this section. Ridges are intuitively amplitude areas that concentrate most of the energy in the combined time–frequency plots. Various methods can be used for ridge extraction. This paper utilises the so-called crazy climbers algorithm based on Markov chain Monte Carlo simulations and explained in [8]. The main idea is as follows: A large number of particles (climbers) are initially randomly seeded in the domain D at time step zero. Then, each climber evolves according to a Markov chain in the D domain with a transition mechanism depending on local D-domain values. The amplitude of the wavelet-based FRF estimator H1(a,b) is the D domain in the case investigated. Climbers are encouraged to “climb on the hills”, to reach the ridges by the Hastings–Metropolis algorithm and a temperature schedule similar to the simulated annealing algorithm.

328

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

Fig. 2. Series of randomly spaced impacts used for excitation of time-variant systems – the time (a); frequency (b) and combined time–frequency (c) domain representations.

Climbers move independently of each other. Therefore the entire idea can be explained using the movement of only one climber. This movement is iterative and can be described using the following few steps: (1) At t ¼ 0, the random position Xð0Þ of the climber on the grid is chosen as Γ ¼ fMðj; kÞ; j ¼ 1; …; J; k ¼ 1; …; Kg

ð21Þ

where Γ is the grid of points described by the wavelet-based FRF estimator H1(a,b), j is discretised scale (frequency) operator, J is the discretized scale (frequency) maximum value, k is discretised translation (time) operator, K is discretised translation (time) maximum value. At this step the temperature T 0 , have to be initialised as well. 0 (2) Given the position XðtÞ ¼ ðj; kÞ at time t, the position Xðt þ 1Þ ¼ ðj0 ; k Þ at time tþ1 is determined according to the following rule: 0 (a) n If k Z1 and k rK  2, then k ¼ k þ1 with probability 1=2 and k 1 with probability 1=2 . 0 n If k ¼ 0 then k ¼ 1. 0 n If k ¼ K 1 then k ¼ K 2. (b) Once move in k direction is done, the possible movement in j direction is considered. The climber tries to move j' ¼ j þ1 or j 1. However, unlike the movement in k direction, movement in j direction is not necessarily made; the probability of the j direction movement depends on the temperature T t . 0 0 0 0 0 n If Mðj ; k Þ 4Mðj ; kÞ, then the j direction movement is done: Xðt þ1Þ ¼ ðj ; k Þ. 0 0 0 0 0 n If Mðj ; k Þ oMðj ; kÞ, then Xðt þ 1Þ ¼ ðj ; k Þ with the probability equal to   1 0 ð22Þ pt ¼ exp ðMðj0 ; k Þ  Mðj0 ; kÞÞ Tt and Xðt þ1Þ ¼ ðj0 ; kÞ with the probability equal to 1  pt . The temperature is then updated. (3) The second step is then repeated until climber temperature is below certain threshold. It is necessary to set a temperature schedule in such a way that when the temperature goes to zero, the trajectories concentrate on the ridges. Points extracted using the crazy climbers algorithm leads to a new time–frequency representation called the occupation measure O. This representation informs how often given time–frequency points are visited by all climbers. In other words, the occupation measure shows how attractive certain paths are for climbers. More detailed information on the crazy climbers algorithm can be find in [8]. These paths also include ridges that are important for the wavelet-based FRF analysis. The entire analysis often leads to multiple ridges. In such situations – for a given time moment – only one ridge is considered. An optimisation algorithm is used to fulfil this requirement. Ridges can be expressed as onedimensional functions that can be presented in the following form: b-φðbÞ

ð23Þ

where b is time variable, φðbÞ is the instantaneous frequency of ridge. The entire procedure leads to the problem of functional optimisation and can be described mathematically using the following expression [38]: Z ε1 ðφ Þ ¼  OðφðbÞ; bÞdb ð24Þ

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

329

where Oðb; φðbÞÞ is a path on the occupation measure where ridges are build. An important assumption – that ridges are continuous and slowly varying functions – has to be made. This is due to the fact that mechanical vibrating systems do not change their physical properties very rapidly. Following this assumption, the ridge “chaining” problem is in fact an optimisation problem. The task is to minimise the value of the following functional [38]: Z Z ε2 ðφðbÞÞ ¼ λ jφ0 ðbÞj2 db þ μ jφ″ðbÞj2 db ð25Þ where λ and μ are weight coefficients which determine smoothness of ridges. Finally Eqs. (24) and (25) can be assembled to obtain Z Z Z ð26Þ εðφðbÞÞ ¼  OðφðbÞ; bÞdb þ λ jφ0 ðbÞj2 db þμ jφ″ðbÞj2 db Once the entire functional is assembled the weight coefficients λ and μ are more meaningful. Both parameters should be selected using a trial and error approach to obtain the trade-off between ridge smoothness and following the path of occupation measure maximum amplitudes by the ridge. When short signals are analysed computational costs involved in the entire procedure are relatively small. These costs can grow very rapidly for long vibration signals. Then it is easier and more time efficient to fit some continuous curves, e.g. polynomials, to given data points. These data points can be obtained using the following iterative steps. The first step is to define parts of the occupation measure where potential ridges are located. One has to select the region defined by lower and upper bounds of time variable bl and bu and lower and upper bounds of frequency variable φl and φu . The second step is to select maxima of frequency variables for the given time instance bk . In first iteration bk have to be equal to bl , i.e. φðbk Þ ¼ maxðOðφl Cφu ; bk ÞÞ

ð27Þ

This step is repeated with changed time variable bk ¼ bk þ 1, until the end of the selected part of occupation measure is reached, i.e. bk ¼ bu . In such a way the obtained ridge is the best solution for penalty function given by Eq. (24) for given region of occupation measure. Once this initial solution is obtained, one can try to fit some polynomial curve to this data, since such a curve will assure smoothness of the finally obtained ridge. It is always advisable to exclude outliers, before fitting the polynomial curve to data. Fitting of polynomial curve should start at the 1st order polynomial and proceed to higher polynomial orders until satisfactory results are obtained. Result of curve fitting should be considered as a ridge. Once a ridge is constructed, mode shape for the analysed time-variant system can be estimated from the amplitude and phase of the wavelet-based FRF corresponding to the ridge. It is clear that when a series of random impacts is used as an excitation, the best estimates are obtained for the neighbourhood corresponding to the excitation points. In practice when experimental tests are performed it is important to avoid double impacts since these impacts can lead to erroneous results. It is clear that impact time instances bi have to be selected using the excitation signal used. Only such time instances are considered for any further analysis; values in between these time instances can be approximated. This approach is possible due to the fact that the current study deals with physical systems and these systems are assumed to be slowly varying. The combination of these time instances and previously obtained ridges define points on the wavelet based FRF which are used to obtain mode shapes, i.e. φi ¼ φðbi Þ

ð28Þ

Until now, entire procedure is global and performed for the most representative point only. However, for a given structure mode shapes are considered to be global parameters defined by local values of amplitude and phase. This is why the entire procedure has to be repeated for all measured points. For a given point P k , evolution of amplitude and phase with time can be defined by as P ka ðbi Þ ¼ H 1a ðφi ; bi ÞfP k g

ð29Þ

P kp ðbi Þ ¼ H 1p ðφi ; bi ÞfP k g

ð30Þ

Once evolutions of amplitude and phase for all points are extracted, the process of amplitude normalisation of amplitude needs to be done. For any, preferably the most representative point P m , the amplitude should be set to the value of 1 for all time variables. The amplitude of all other points is then normalised by dividing the amplitude at a given point and time instance by the amplitude of P m , i.e. P ka ðbi Þ ¼

P ka ðbi Þ P ma ðbi Þ

ð31Þ

P ma ðbi Þ ¼1 P ma ðbi Þ

ð32Þ

P ma ðbi Þ ¼

The last unknown modal model parameter is damping. Since impact excitation is used, this parameter can be estimated using the exponential decay of response. This methodology requires decoupling of modes with the use of the Continuous Wavelet Transform. Once modes are decoupled, the envelope of particular mode can be analyzed. It is well-known that

330

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

system response to impact excitation will result in impact response which can be described mathematically as qffiffiffiffiffiffiffiffiffiffiffiffi  1 ζ2 ωt þθ yðtÞ ¼ Ae  ζωt sin

ð33Þ

where Ae  ζωt is the envelope of impulse response, ζ is the damping ratio, θ is the phase. Damping ratio can be found with the use of curve fitting techniques. More details on this methodology are given in [9]. 4. System identification for simulated vibration data This section illustrates the wavelet-based FRF for the simulated vibration data. A simple 2-DOF time-variant system – illustrated in Fig. 3 – was analysed. This linear, highly dumped system included one time-varying mass element. All other relevant physical parameters were constant. The values of physical parameters used were: m2 ¼10 [kg], k1 ¼3000000 [N/m], k2 ¼800000 [N/m], c1 ¼150 [N/(m/s)] and c2 ¼450 [N/(m/s)]. The m1 mass element was decreasing linearly in time following the function m1(t) ¼16–0.7t [kg]. As a result natural frequency of the analysed system was increasing in time. The 2-DOF system was simulated using the Matlab/Simulink platform. The sampling frequency was equal to 1 kHz in these simulations. Firstly, the system was simulated – as a time-invariant system – for all time instances involved, i.e. from 0 to 10 s. Altogether, 10001 numerical simulations were performed with constant physical parameters. Impact excitation was used without any noises in order to obtain best results. These parameters were the same for all time instances, except the value of mass m1 that was different for all simulations, following the m1(t) function. Once these simulations were performed, the classical FRF were calculated for these 10001 time-invariant systems. The results – assembled for all time instances – are shown in Fig. 4. Here, a cascade of amplitudes and phases of the classical FRFs reveal the expected (theoretical) change of natural frequencies. Fig. 5 shows examples of the classical FRFs calculated for two extreme time instances, i.e. for the values of time equal to 0 and 10 s. These classical FRFs were used to estimate (theoretical) natural frequencies. The results in Figs. 4 and 5 show that the frequencies change from 27 to 32 and from 113 to 131 Hz for the first and second vibration modes respectively. Secondly, the time-variant 2-DOF system was simulated using a series of random impacts as excitation. The waveletbased FRF was calculated following the procedure described in Sections 2 and 3. The results are presented in Fig. 6. Here, the amplitude and phase of the wavelet-based FRF exhibit broken patterns corresponding to a series of impacts. These results can be compared Fig. 4. Although, the results reveal increasing natural frequencies for both vibration modes – as expected – additional signal processing is required.

Fig. 3. Schematic diagram of the simulated 2-DOF time-variant system.

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

331

Fig. 4. Cascade of classical FRFs for a series of simulated time-invariant 2-DOF systems calculated for various time instances: (a) amplitude and (b) phase.

Fig. 5. Classical FRFs calculated for the two simulated time-invariant 2-DOF systems for t ¼0 s (solid lines) and t¼ 10 s (dashed lines): (a) amplitude and (b) phase.

Fig. 6. Wavelet-based FRF for the simulated time-variant 2-DOF system (initial results): (a) amplitude and (b) phase.

332

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

The crazy climbers algorithm was applied, as described in Section 3.3. As a result the occupation measure was obtained, as shown in Fig. 7a. Then the optimisation algorithm was used to perform chaining and obtain ridges corresponding two vibration modes. The final wavelet-based FRF is presented in Fig. 7b. Two vibration modes can be clearly observed in the wavelet-based FRF. Frequencies of these modes increase with time, as expected. The frequency of the first vibration mode changes from 27 Hz for 0 s to 31 Hz for 10 s whereas the frequency of the second vibration mode changes from 113 Hz for 0 s to 130 Hz for 10 s. When Fig. 7b is compared with Fig. 4a, very good agreement with the expected results can be observed. The second vibration mode was selected as an example to perform identification. Firstly, amplitude and phase values of the wavelet-based FRF were recovered for points on the ridge corresponding to this vibration mode. The amplitude values were normalised to the value of m2. Then the second mode shape was calculated and damping estimated, following the procedure described in Section 3.3. The identification results for the second vibration mode are presented in Fig. 8. Reference (expected or theoretical) results are given using solid lines, whereas estimated results are given using circles. The results in Fig. 8a show that the ratio of vibration amplitudes of both masses changed from around 0.7 to 1.1–1.3 with the change of time from 0 to 10 s. The results in Fig. 8b show that vibrating masses remain out of phase. The results in Fig. 8a and b show that the second mode shape was identified correctly. A small discrepancy for amplitude and phase was obtained between 7 and 10 s. Fig. 8c presents the identified natural frequency. The results show that the natural frequency increases from 113 to 131 Hz for the analysed period of time, as expected. Almost perfect agreement was obtained when the expected and estimated natural frequencies are compared. Damping estimation results are presented in Fig. 8d. The results show estimated value of damping changes 0.92–1.48% for the analysed period of time. These values are slightly underestimated if compared with the reference (theoretical) damping change from 0.96% to 1.66%. Nevertheless the estimated trend is correct.

Fig. 7. Wavelet-based FRF for the simulated time-variant 2-DOF system (after post-processing): (a) occupation measure and (b) final result - FRF's ridges.

Fig. 8. Modal parameters estimated from the wavelet-based FRF (shown in Figs. 6–7) for the second vibration mode: (a) mode shape - amplitude, (b) mode shape - phase, (c) natural frequency and (d) damping. Solid line shows reference (theoretical) results; circles indicate estimated results.

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

333

In summary, the identification results presented in Fig. 8 show that good agreement was achieved between the estimated and reference modal parameters. 5. System identification for experimental vibration data A simple experiment was conducted to obtain vibration data from a time-variant system. The system analysed was a time-variant frame-like structure shown in Fig. 9. The structure consisted of two beam elements and a tank. The vertical and horizontal beams were 80 cm and 60 cm long respectively. The tank – fixed to the free end of the horizontal beam – was filled with sand. The structure was excited using a modal hammer. A series of random impacts were applied to the structure, as described in Section 3.2. Once the structure was excited, vibration measurements were taken – using accelerometers – for 12 points uniformly distributed on the structure. When measurements were taken, sand was draining from the tank. As a result the initial mass of the tank was reduced from 3 to 0 kg within the period of 65 s. All measurements – i.e. excitation signals and vibration responses – were recorded using modal hammer and 24 accelerometers evenly distributed on frame, measuring in two directions. Although a number of tests were conducted to check the repeatability, data from only one test were used for further analysis. In other words, averaging was not applied, as explained in Section 2. Firstly, the system was tested when the tank was full and then when empty. These two cases represented a timeinvariant system for two different (i.e. extreme) values of tank masses - maximum and minimum. The data gathered from these two experiments can serve as a reference data from a time-invariant system. Since the classical modal analysis was used, averaging was performed to improve signal-to-noise ratio. Modal analysis and identification was conducted using the TestLabs software. Fig. 10 presents the amplitude of the classical FRF for the two cases investigated. Here, solid and dashed lines are used to present the FRF amplitude for the structure with the full and empty tank, respectively. The analysed FRFs are different due to different masses of the tank, as expected. The major - and most interesting - difference can be observed in the 250–300 Hz range. This area in encircled in Fig. 10 and zoomed in Fig. 11. When the tank was full one vibration mode (natural frequency equal to 276 Hz) was identified in the FRF. In contrast, two vibration modes (natural frequencies equal to 261 and 289 Hz) can be observed in the FRF for the empty tank in Figs. 10 and 11. In theory, the reduced mass of the tank should result in the increased natural frequency. Therefore, one would expect that the frequency of the 276 Hz vibration mode was shifted to 289 Hz. However, further analysis – shown in Fig. 12 – shows that the mode shapes for the mode 276 Hz (full tank) and 261 Hz (empty tank) are identical. Clearly, the classical FRF does not capture the dynamics of the system properly when the tank is drained. The values of damping were also estimated for these two modes as 1.28 and 3.10% for structure with the full and empty tank. The second test involved the experimental analysis of the time-variant frame-like system with the tank. The structure was excited using a series of random impacts, as explained in Section 3.3. This time sand was draining from the tank when measurements were taken. The 70 s long response sequence was captured. Since the data were sampled at 2048 Hz frequency rate, altogether 1400000 data samples were acquired. The dynamic analysis was limited to the 250–300 Hz frequency range to

Fig. 9. Frame-like structure with the tank.

334

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

10.0 5.0 3.0 2.0 1.0 Log

0.5 0.3 0.2 0.1 50e-3 30e-3 20e-3 10e-3 0

100

200

300

400

500

600

700

800

900

1024

Hz Fig. 10. FRF amplitude for the frame-like time-invariant structure with the tank - full tank (solid line); empty tank (dashed line).

10.0 6.0 4.0 3.0

Log

2.0 1.0 0.6 0.4 0.3 0.2 0.1

261

276

289

250 255 258 260 263 265 268 270 273 275 278 280 283 285 288 290 293 295 300 Hz Fig. 11. The same as Fig. 10 but the analysis is zoomed for the 250–300 Hz frequency range.

Fig. 12. Vibration mode shapes for the time variant frame-like structure with the tank: (a) full tank - 276 Hz and (b) empty tank - 261 Hz. The results were obtained from the classical FRF analysis.

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

335

further analyse the vibration modes shown in Figs. 10 and 11. The wavelet-based FRF was calculated - following the procedure described in Sections 2 and 3. Fig. 13 shows the amplitude and phase of the wavelet-based FRF. Since only one response sequence was analysed (and averaging not used), the results are quite noisy and do not reveal anything interesting. Therefore additional post-processing was applied. The crazy climbers algorithm applied to the waveletbased FRF amplitude. This analysis led to the occupation measure shown in Fig. 14a. Two vibration modes can be already observed in the analysed frequency range. The optimisation procedure was then used to chain the results and produce ridges corresponding to the analysed vibration modes. Fig. 14b shows the final result for the amplitude of the wavelet-based FRF. Two vibration modes – with varying frequencies – can be clearly observed in the results. The results show that the 276 Hz vibration modes shifts towards the lower frequency of 258 Hz. The second observed vibration mode of 292 Hz emerged after approximately 20 s. The frequency of this mode decreased to 289 Hz after further 50 s. The wavelet-based FRF amplitude in Fig. 14b can be compared with the classical FRF in Fig. 11. The waveletbased FRF – shown in Fig. 14b – correctly captures the time-variant dynamics of the structure explaining the transfer in Fig. 11 from one vibration mode for the solid tank to two vibration modes for the empty tank. Identification procedures – described in Section 3.3 – were additionally applied for the first analysed (276 Hz) vibration mode. Fig. 15a and b presents amplitude and phase values respectively for all wavelet-based FRF points captured at time instances corresponding to excitation impacts. Relatively large amplitude fluctuations can be observed in Fig. 15a in some amplitude characteristics. Similarly large phase shifts can be observed in Fig. 15b for measuring points corresponding to small amplitudes. These points were located close to the nodes of the relevant mode shapes. Estimated natural frequencies for various time instances are shown in Fig. 15c. The results show that the natural frequency changes from 276 Hz (full tank) to 258 Hz (empty tank) in 70 s. Fig. 15d

Fig. 13. Wavelet-based FRF for the experimental time-variant frame-like structure with the draining tank (initial results): (a) amplitude and (b) phase.

Fig. 14. Wavelet-based FRF for the experimental time-variant frame-like structure with the draining tank (after post-processing): (a) occupation measure and (b) final result - FRF's ridges.

336

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

Fig. 15. Modal parameters for the experimental time-variant frame-like structure with the tank. These results were estimated from the wavelet-based FRF shown in Figs. 13–14 for the first analysed vibration mode: (a) mode shape - amplitude, (b) mode shape - phase, (c) natural frequency and (d) damping.

Fig. 16. Vibration mode shapes for the time-variant frame-like structure with the tank: (a) full tank - 276 Hz and (b) empty tank - 261 Hz. The results were obtained from the wavelet-based FRF analysis.

gives the estimated value of damping for the analysed vibration mode. The value of damping changed from 0.94 (full tank) to 3.00% (empty tank) within the analysed 70 s period. The results shown in Fig. 15a and b were used to reconstruct the relevant mode shape. The results for the full (0 s) and empty tank (70 s) are presented in Fig. 16, respectively. These results resemblance the modes shapes presented in Fig. 12. In summary, in contrast to the classical FRF analysis shown previously, the wavelet-based FRF properly captures the dynamics of the time-variant system. Also, when system identification results (i.e. natural frequencies, damping and mode shapes) obtained from the wavelet-based FRF amplitude for the two extreme cases (i.e. for the full and empty tank) are compared with the relevant results shown for the classical FRF, good agreement was achieved.

6. Conclusion A new method for identification of modal model for time-variant structures was used. The method – based on the recently proposed wavelet-based FRF – enables one to identify natural frequencies, damping ratios and mode shapes evolving in time. The method has been illustrated using simulated and experimental vibration data for MDOF systems. The method does not require any averaging that is not easy for time–frequency characteristics. The results presented in the paper show that the method used captures correctly the dynamics of the analysed time systems. All modal parameters – and their evolution in time – were identified correctly. This method has tested only for a special type of random impact excitation. Another limitation is connected with the need for measuring of excitation signals. Although, signal/data averaging is not required it is important to note that the analysis can led to uncertainties in estimated parameters. This problem needs further investigations.

K. Dziedziech et al. / Mechanical Systems and Signal Processing 50-51 (2015) 323–337

337

More research work needs to be performed to confirm the results given in this paper and establish the method presented. Further studies should involve more advanced identification techniques for the wavelet-based FRF, other types of broadband excitation. The method should be also tested for more complex structures – particularly with close modes and nonlinearities – and for operational data. Acknowledgements The work presented in this paper was supported by funding from the WELCOME research project no. 2010-3/2, sponsored by the Foundation for Polish Science (Innovative Economy, National Cohesion Programme, EU). References [1] A.G. Poulimenos, S.D. Fassois, Parametric time-domain methods for non-stationary random vibration modelling and analysis - a critical survey and comparison, Mech. Syst. Signal Process. 20 (4) (2006) 763–816. [2] M.D. Spiridonakos, S.D. Fassois, Parametric identification of a time-varying structure based on vector vibration response measurements, Mech. Syst. Signal Process. 23 (6) (2009) 2029–2048. [3] E. Wigner, On the quantum correction for thermodynamic equilibrium, Phys. Rev. 40 (1932) 749–759. [4] L. Cohen, Time–frequency analysis, vol. 778, Prentice Hall PTR, New Jersey, 1995. [5] A.W. Rihaczek, Principles of High-Resolution Radar, Artech House Radar Library. Artech House, Norwood, MA, USA, 1996. [6] H.I. Choi, W.J. Williams, Improved time–frequency representation of multicomponent signals using exponential kernels, IEEE Trans. Acoust. Speech Signal Process. 37 (6) (1989) 862–871. [7] F. Hlawatsch, F. Auger, Time–frequency analysis: concepts and methods. Digital signal and image processing series, ISTE, London, UK, 2008. [8] R. Carmona, W.L. Hwang, B. Torrésani, Practical Time–frequency Analysis: Gabor and Wavelet Transforms With an Implementation in S. Wavelet Analysis and Its Applications, 9, Academic Press Inc., San Diego, CA, USA, 1998. [9] W.J. Staszewski, Identification of damping in md of systems using time-scale decomposition, J. Sound Vib. 203 (2) (1997) 283–305. [10] S.L. Chen, J.J. Liu, H.C. Lai, Wavelet analysis for identification of damping ratios and natural frequencies, J. Sound Vib. 323 (1–2) (2009) 130–147. [11] T.P. Le, P. Argoul, Continuous wavelet transform for modal identification using free decay response, J. Sound Vib. 277 (1-2) (2004) 73–100. [12] W.J. Staszewski, Identification of non-linear systems using multi-scale ridges and skeletons of the wavelet transform, J. Sound Vib. 214 (4) (1998) 639–658. [13] J. Urbanek, T. Barszcz, J. Antoni., A two-step procedure for estimation of instantaneous rotational speed with large fluctuations, Mech. Syst. Signal Process. 38 (1) (2013) 96–102. [14] B. Basu, S. Nagarajaiah, A. Chakraborty, Online identification of linear time-varying stiffness of structural systems by wavelet analysis, Struct. Health Monit. Int. J. 7 (2008) 21–36. [15] T. Uhl, M. Petko, G. Karpiel, A. Klepka, Real time estimation of modal parameters and their quality assessment, Shock Vib. 15 (3) (2008) 299–306. [16] X. Shan, J.B. Burl, Continuous wavelet based linear time-varying system identification, Signal Process. 91 (6) (2011) 1476–1488. [17] S. Nagarajaiah, B. Basu, Output only modal identification and structural damage detection using time frequency & wavelet techniques, Earthq. Eng. Eng. Vib. 8 (4) (2009) 583–605. [18] A.N. Robertson, B. Basu., Wavelet Analysis, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2009. [19] K. Liu, Extension of modal analysis to linear time-varying systems, J. Sound Vib. 226 (1999) 149–167. [20] R. Pedersen, I. Santos, I.A. Hede, Advantages and drawbacks of applying periodic time-variant modal analysis to spur gear dynamics, Mech. Syst. Signal Process. 24 (5) (2010) 1495–1508. [21] J.A. Ball, I. Gohberg, M.A. Kaashoek, A frequency response function for linear, time-varying systems, Math. Control Signals Syst. 8 (4) (1995) 334–351. [22] T. Prothmann, Estimation of partial frequency response functions of periodic time-variant systems under broadband excitation, Proc. Appl. Math. Mech. 3 (1) (2003) 336–337. [23] G. Matz, F. Hlawatsch, Time–frequency transfer function calculus (symbolic calculus) of linear time-varying systems (linear operators) based on a generalized underspread theory, J. Math. Phys. 39 (8) (1998) 4041–4070. [24] R. Zou, K.H. Chon, Robust algorithm for estimation of time-varying transfer functions, IEEE Trans. Biomed. Eng. 51 (2) (2004) 219–228. [25] K. Dziedziech, W.J. Staszewski, T. Uhl, Time–frequency analysis of time-variant systems, Diagnostyka 14 (2013) 37–42. [26] M.B. Priestley, Non-linear and non-stationary time series analysis, Academic Press Inc., San Diego, CA, USA, 1988. [27] P. Tratskas, P.D. Spanos, Linear multi-degree-of-freedom system stochastic response by using the harmonic wavelet transform, J. Appl. Mech. 70 (5) (2003) 724–731. [28] G. Longo, B. Picinbono, International Centre for Mechanical Sciences, Time and frequency representation of signals and systems. Courses and lectures International Centre for Mechanical Sciences, Springer Verlag Wien, New York, 1989. [29] D. Tjøstheim, Spectral generating operators for non-stationary processes, Adv. Appl. Probab. (1976) 831–846. [30] G. Melard Processus Purement Indetérminables á Paramétre Discret: Approaches Fréquenttielle et Temporelle. (Ph.D. thesis), Université Libre de Bruxelles, 1975. [31] K. Dziedziech, W. J. Staszewski, and T. Uhl. Input–output time–frequency analysis of time-variant systems. In Proceedings of the International conference on Noise and Vibration Engineering, 2012. [32] W.J. Staszewski, A.N. Robertson, Time–frequency and time–scale analyses for structural health monitoring, Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci. 365 (1851) (2007) 449–477. [33] W. J. Staszewski and J. Giacomin. Application of the wavelet based frfs to the analysis of nonstationary vehicle data. In Proceedings — SPIE the International Society for Optical Engineering, pp. 425–431, 1997. [34] W.J. Staszewski, D.M. Wallace, Wavelet-based frequency response function for time-variant systems — an exploratory study, Mech. Syst. Signal Process. 47 (1-2) (2014) 35–49. [35] K. Dziedziech, W.J. Staszewski, T. Uhl, Wavelet-based frequency response function – comparative study of input excitation, Shock Vib. (2014), in press. [36] R.S. Pathak., The Wavelet Transform. Atlantis studies in mathematics for engineering and science, Atlantis Press, 2009. [37] W.J. Staszewski, Wavelets for mechanical and structural damage identification, In Polish Academy of Sciences, IMP PAN, 2000. [38] R.A. Carmona, W.L. Hwang, B. Torresani, Multiridge detection and time–frequency reconstruction, IEEE Trans. Signal Process. 47 (2) (1999) 480–492.