Detection and removal of harmonic components in operational modal analysis

1. Introduction

During the last decade, the modal analysis becomes a key technology in structural dynamics behavior study [1, 2]. Operational modal analysis (OMA) is a technology for the identification of the linear structural modal parameters by using only the output responses [3, 4]. It is used to derive an experimental dynamics model from vibration measurements on a structure in operational conditions. In many operational cases, it is impossible to measure the excitation force, and the structural responses are the only information that can be collected for the system identification algorithms. A new method based on the transmissibility function is proposed to process the parameters in OMA with unrestricted types load [5, 6], but it needs twice tests with different loads and may generate fake resonance peaks. Therefore, in the OMA, the system is generally assumed to be excited by a broadband random signal [7, 8]. However, in many applications, the presence of periodic excitation in the random loads is unavoidable. The periodic forces generated by unbalanced masses of rotating parts or aerodynamic or electric actuators, may cause harmonic components in the responses [9]. These harmonic components will cause uncertainty in extraction of modal parameters [10], and need to be detected and removed before modal identification.

A direct approach for processing the harmonic components in outputs is to take them as zero-damping modes. In practice, this method would not be effective when the structural modes have very low damping, or the harmonic frequencies are very close to the structural modal frequencies [11]. A method based on random decrement is proposed for distinguishing the harmonic frequencies from the natural frequencies [12, 13]. It argues that the harmonic components are un-attenuated. The harmonic components, in other words, are also assumed to be zero-damping modes in this method. Another method based on the probability density function (PDF) of the output for distinguishing a harmonic response from narrowband stochastic responses is widely used [9, 14, 15]. This method, however, is effective when the harmonic frequencies are well separated from the structural frequencies, since it would imply that temporal filter must be used. Furthermore, it will require tedious calculations and be restricted by bandwidth of the filter. A kurtosis-based method was initially proposed in time domain for distinguishing the random and period signals [16-18], and processed in a way similar to the PDF method mentioned above. In order to avoid tedious data calculation, the spectral kurtosis (SK) of a signal is defined as the kurtosis of its frequency component, and has been applied to harmonic components detection in the OMA [19-21]. Another advantage of this method is that the information of each frequency component can be indicated promptly in the frequency band of interest. However, the method is not perfect, since only one channel signal can be dealt with at a time. It brings a time-consuming problem and may lead to a failure of detection.

In order to improve the quality and efficiency of the harmonic components detection, an enhanced spectral kurtosis (ESK) based on the enhanced PSD function is proposed in this paper. Unlike the original SK method, the ESK calculation utilizes all multiple-channel signals at each resonance frequency of interest. Moreover, due to the application of a PSD enhancement technique, the band-pass filter is not required any more. Starting from the enhanced PSD obtained from the frequency spatial domain decomposition (FSDD) [22, 23], the ESK algorithm is derived in detail. Then, a new method named the virtual excitation technique (VET) is presented for the removal of harmonic components even when they are close to the natural frequencies. By using multiple input & multiple output (MIMO) model and virtual excitation assumption after detection of the harmonic frequencies, the responses which do not contain the harmonic components are obtained, and the structural model parameters are identified with these responses by the FSDD method. The result from the proposed method is also compared with that from the previous method by a numerical simulation. Finally, an application of the proposed method to detect and remove the harmonic component of a GARTEUR (Group of Aeronautical Research and Technology in EURope) aircraft model is presented to demonstrate its procedures and credibility.

2. Theory background

2.2. VET algorithm

The unknown mixed signal $x$ can be considered as the linear superposition of a random signal and a harmonic signal:

10

$x\left(n\right)=x_b\left(n\right)+x_s\left(n\right),$

where $x_b$ is a random signal and $x_s$ is a harmonic signal. In the frequency domain, the random signal ${\mathbf{X}}_b\left(\omega \right)$ can be assumed as the product of a virtual white noise ${\mathbf{X}}_p\left(\omega \right)$ (constant spectrum) and a function of coefficient $C_1\left(\omega \right)$. Similarly, the harmonic signal ${\mathbf{X}}_s\left(\omega \right)$ can be assumed as the product of a virtual harmonic signal ${\mathbf{X}}_h\left(\omega \right)$ and a constant coefficient $c_2$. The DFT of the mixed signal $x$ can be expressed with the virtual excitations as:

11

$\mathbf{X}\left(\omega \right)={\mathbf{X}}_b\left(\omega \right)+{\mathbf{X}}_s\left(\omega \right)=C_1{\left(\omega \right)\mathbf{X}}_p\left(\omega \right)+c_2{X}_h\left(\omega \right).$

Indeed, the spectrum of the random input is not always flat in practice, but can be generally considered as constant spectrum in a narrow band [23]. The MIMO model can be written at each frequency in the narrow band as:

12

$\underset{\omega \to {\omega }_k}{\lim }\mathbf{Y}\left(\omega \right)=\mathbf{H}\left(\omega \right)\mathbf{X}\left(\omega \right)=\mathbf{H}\left(\omega \right)\left({\stackrel{-}{c}}_1{{\mathbf{X}}_p\left(\omega \right)+c}_2{\mathbf{X}}_h\left(\omega \right)\right)$
$=\left[\begin{array}{cc}{\stackrel{-}{c}}_1\mathbf{H}\left(\omega \right)& c_2\mathbf{H}\left(\omega \right)\end{array}\right]\left[\begin{array}{c}{\mathbf{X}}_p\left(\omega \right)\\ {\mathbf{X}}_h\left(\omega \right)\end{array}\right],$

where $\mathbf{H}$ is the frequency response function (FRF) matrix of the system, and is a constant coefficient. Using the Least Squares algorithm, the proportional virtual responses of this system can be obtained as:

13

$\left[\begin{array}{cc}\stackrel{-}{c_1}\mathbf{H}\left(\omega \right)& c_2\mathbf{H}\left(\omega \right)\end{array}\right]=\mathbf{Y}\left(\omega \right){\left[\begin{array}{c}{\mathbf{X}}_p\left(\omega \right)\\ {\mathbf{X}}_h\left(\omega \right)\end{array}\right]}^{+},$

where the superscript + denotes pseudo inversion, and $\omega ⟶{\omega }_k$ Therefore, the response by the virtual excitations can be obtained as:

14

$\tilde{\mathbf{Y}}\left(\omega \right)=\stackrel{-}{c_1}\mathbf{H}\left(\omega \right).$

It is obvious that the $\tilde{\mathbf{Y}}\left(\omega \right)$ does not contain the harmonic modes any more, which means the harmonic components are removed successfully.

In order to improve the understandability of the method, a flow chart is given as follows:

Fig. 1Flow chart of identification based on ESK and VET

3. Numerical case

A two-storey building structure with six DOFs (regarding two floors as two lumped masses, each lumped mass has two translational DOFs and a rotational DOF), shown in Fig. 2, is used to validate the proposed algorithm. The simulation technique based on the FRF modes is employed to generate six responses with five uncorrelated random noise signals and a mixed signal with a periodic frequency at 33.4 Hz in the $X$-direction. The simulated time domain responses have a length of 204, 800 points, with a sample frequency of 256 Hz. The PSD matrix of responses will be computed by the Welch method, using Hanning windowing and resolution of 0.2 Hz. Then, the mode indicator function (MIF) curves can be computed by singular value decomposition (SVD) in 0-100 Hz, as shown in Fig. 3. The original SK method is used to detect a harmonic component and the results are shown in Fig. 4.

Fig. 2Two-storey building model

Fig. 3Mode indicator curves of two-storey building

Fig. 4SK curves of responses (two-storey building)

a)

b)

c)

d)

e)

f)

g)

h)

In Fig. 3, the first (solid) curve, second (dash) curve and third (dotted) curve consist of the 1st, 2nd and 3rd singular value at each spectral line, respectively. Six structural modes and a harmonic component (at the frequency of 33.4 Hz) are apparently indicted by the mode indicator curves. Moreover, the harmonic component is very close, actually only three spectral lines are far, to the 3rd natural frequency. Fig. 4 shows that the harmonic is not always indicated by SK of response at each DOF, since the harmonic component(s) may not, in general, be strongly contained in all measured responses. Therefore, the SK of all responses should be calculated for a good harmonic detection. Moreover, although the harmonic is indicated in several SK curves, like Fig. 4(e) and Fig. 4(f), the SK values at the harmonic frequency are not exact –1 as expected. Since the frequency band is too narrow to separate harmonic frequency from nearby structural frequency by band-pass filtering, the computational accuracy is influenced.

The $k$th ($k=$1, 2,…, 7) ESK curve is computed by Eq. (6) and Eq. (8) after picking the corresponding peak in the MIF curve, as shown in Fig. 5. It is clearly observed that the 4th (harmonic component) ESK value at 33.4 Hz is very close to –1, which means that the harmonic component has been indicated successfully. In contrast with the SK method, the ESK method is more reliable to detect the harmonic components, because all of the measured responses are completely considered by the calculation of ESK from the enhanced PSD.

Fig. 5ESK curves at each peak (two-storey building)

Furthermore, Fig. 5 shows that the harmonic frequency is also indicated by other ESK curves (such as the first and third one), since the harmonic component is not always completely decoupled from the structural modes. As a result, the enhanced PSD of a structural mode may contain the harmonic component as well. Accordingly, the harmonic component should be removed to avoid the errors during the identification of modal parameters, especially when the harmonic frequency is close to a modal frequency. The virtual harmonic input signal can be assumed as: $x_h=\sin (2\pi \times 33.4t)$ and the random input signal is assumed as virtual white noise $x_p$. According to the Eq. (12) and Eq. (13), the third enhanced power spectral density of virtual response is proportional to that of the real response except at harmonic frequencies in a narrow band, as shown in Fig. 6(a). For comparison, the constant of proportionality is normalized, as shown in Fig. 6(b).

Fig. 6Enhanced PSD curve in narrow band before and after harmonic component removal

a)

b)

It should be noted that they do not match completely (shown in Fig. 6(b)), because the assumed excitation is uncorrelated with the real excitation (see Eq. (10)), which means the assumed excitation is indeed noise. Fortunately, the noise level is acceptable to identify the modal parameters by FSDD method. The modal frequencies and damping ratios identified before and after harmonic components removal are shown in Table 1 and Table 2, respectively.

As shown in Table 1 and Table 2, the relative error of the 3rd modal damming ratio is up to 24 % before the harmonic removal, but it is corrected after the harmonic removal. The modal assurance criterion (MAC) is introduced to check the orthogonality between the theoretical and identified mode shapes, as shown in Fig. 7. The mode shapes identified after harmonic component removal are shown in Fig. 8.

It is seen in Fig. 7 that the orthogonality of the mode shapes is greatly improved after the harmonic removal. The analytical results indicate that the removal of harmonic components is an efficient way to adjust the structural modes. It also proves that the proposed method is an effective tool to remove the harmonic components in OMA.

Fig. 7Modal assurance criterion (MAC)

a) MAC before harmonic components removal

b) MAC after harmonic components removal

Fig. 8Structural modes of two-storey building after harmonic component removal

4. Experimental case

The technique discussed above is applied to the GARTEUR aircraft model in the laboratory. Two shakers are fixed in two directions respectively, as shown in Fig. 9. A random noise of 0-128 Hz is applied on the bottom of the structure by one shaker, and a harmonic load of 39.15 Hz is applied on the middle of the structure by another shaker. The experiment is performed with 102, 400 points at a sampling frequency of 128 Hz. The acceleration responses of thirty-four DOFs are obtained in time domain, as shown in Fig. 10. The PSD matrix of responses is computed by the Welch method, with Hanning window applied. The frequency resolution is 0.08 Hz. The FSDD mode indicator curves of response signals, with nine resonance peaks indicated, are showed in Fig. 10. The first peak is the rigid body mode caused by rubber cord suspension, and the fifth peak is the harmonic component which is very close to the fourth structural mode. The original SKs of the thirty-four measured responses are computed to detect the harmonic component, and part of them is shown in Fig. 11.

Fig. 9Experimental scene

Fig. 10Typical time response (top) and mode indicator curves (bottom) of GARTEUR

Fig. 11SK curves of responses (GARTEUR)

The SK curves in Fig. 11 are not always effective enough to indicate the harmonic component, since the computational accuracy is affected by the band-pass filtering. Therefore, the SKs of all measured responses should be calculated and examined to make sure that there are no missing harmonics in practice. By comparison, there are only nine ESKs that need to be computed, and the results are shown in Fig. 12. It is observed that the harmonic is clearly indicated by the 5th ESK value at 39.15 Hz. It is also noticed that the computation of the SKs takes about 12.68 seconds, while the computation of the ESKs (including SVD) needs only about 3.03 seconds on the same computing platform. It turns out that the ESK method is more reliable and efficient to detect the harmonic components than the original SK method.

Fig. 12ESK curves at each peak (GARTEUR)

The harmonic component may affect the modal parameters identified by FSDD, since it is close to the 5th mode (6th peak) and 6th mode (7th peak) and cannot be decoupled in the enhanced PSD. Therefore, the harmonic component, which has been detected, should be removed by using the VET. The comparison of the enhanced PSD before and after harmonic component removal is shown in Fig. 13.

Fig. 13Enhanced PSD of 5th and 6th mode enhanced PSD curves in narrow band before and after harmonic component removal

As shown in Table 3 and Table 4, the relative error of the 4th and 5th modal damming ratios are up to 57.14 % and 50 %, respectively, before the harmonic removal, but they are modified after the harmonic removal. The MAC is used to check the orthogonality between the reference and identified mode shapes, as shown in Fig. 14.

Fig. 14Modal assurance criterion (MAC)

a) MAC before harmonic components removal

b) MAC after harmonic components removal

The improvement of the orthogonality (note the MAC value of the 4th mode) can also be observed after the harmonic removal. The experimental results prove that the proposed method is an effective tool to remove the harmonic components in the OMA. Then, the structural mode shapes in 0-50 Hz are obtained after harmonic component removal as shown in Fig. 15.

Fig. 15Mode shapes of GARTUER after harmonic components removal

5. Conclusions

Harmonic components make a difference in the dynamical response of many structural or mechanical systems which are excited by periodic loads. It is a challenge for many OMA techniques to detect and remove the harmonic components.

The ESK method is proposed in this paper to optimize the process and result of the harmonic component detection. It takes all multiple-channel signals at each resonance frequency into account by the enhanced PSD function. Less computation than the original SK method is asked, because the number of measured DOFs is generally larger than the number of modes in the OMA. It can also improve the accuracy of computation and avoids the risk of the harmonic component(s) absence by computing the enhanced mode.

The VET method is presented in this paper for the removal of harmonic components. By assuming the Multiple Input Multiple Output (MIMO) model and virtual noise excitation, the responses without harmonic components are obtained by the Least Squares (LS) algorithm. Then the accuracy of identified modal parameters can be greatly improved for the modes highly coupled by the harmonic component.

The results of the simulation and experiment demonstrate that the presented harmonic detection and removal methods are efficient and credible.