Characterization of non-linear bearings using the Hilbert–Huang transform

Abstract

Changes in the performance of bearings can significantly vary the distribution of internal forces and moments in a structure as a result of environmental or operational loads. The response of a bearing has been traditionally idealized using a linear model but a non-linear representation produces a more accurate picture at the expense of modelling complexity and computational time. In this article, a lead rubber bearing is idealized using the hysteretic Bouc–Wen model. The Hilbert–Huang transform is then employed to characterize the features of the non-linear system from the instantaneous frequencies of the bearing response to a time-varying force. Instantaneous frequencies are also shown to be a useful tool in detecting sudden damage to the bearings simulated by a reduction in the effective stiffness of the force-deformation loop.

Introduction

Generally, it is reasonable to assume that the response of a structure to moderate dynamic loads will remain in the elastic region; however, large deformations or specific components with non-linear physical properties are represented best by non-linear models. ¹ An example of a non-linear structural component is an elastomeric bearing. Elastomeric bearings are the most commonly used type of seismic isolator ² and their main purpose is to allow lateral movement due to shrinkage, temperature variation and earthquakes. In addition to this lateral flexibility, bearings can have a degree of compressibility in the vertical direction, typically modelled with translational springs. In some cases, rotational springs are also needed to accurately resemble field measurements from static and dynamic loading tests. For instance, there is experimental evidence that changes in joints and elastomeric bearings are one of the main causes of changes in bridge frequencies.^3–5

Of relevance to this investigation is the non-linear nature of elastomeric bearings due their inherent damping properties. ⁶ The nonlinearity can be described by the restoring force and displacement behaviour, that is, a hysteresis loop describing the recoverable and permanent deformations. A number of mathematical hysteresis models that predict the bearing response for different levels of stress are presented in the literature.^7–11 Hwang et al. ⁷ experimentally validate a mathematical model that replicates the hysteretic behaviour of hard rubber bearings for a variety of external loading conditions. Matsagar and Jangid ¹² investigate the influence of different isolator characteristics, such as the shape of the isolator force-deformation loop and yield displacement of the isolator, on the seismic response of the structure. They find that the seismic response is significantly influenced by the shape of the loop and that lower yield displacements increase the acceleration response of the structure. In this article, elastomeric bearings are assumed to be lead rubber bearings modelled as vertical springs. Lead rubber bearings present several convenient features with fundamental drawbacks related to the negligible increase in damping and high deformability for low static loads. The lead plug increases the horizontal stiffness significantly compared to other bearings (e.g. laminated rubber bearings). With the placing of lead plugs (Figure 1(a) where $f(t)$ is the external vertical load applied to the bearing), aside from energy dissipation for seismic response and stiffness for static loads, a single compact device is obtained which is able to satisfy most of the requirements associated with a good bearing. For that reason, rubber bearings are used extensively in practical applications in bridge structures. ¹³

Figure 1.

Typical lead elastomeric bearing used in bridges: (a) components and (b) force–displacement hysteresis loops.

Figure 1(b) represents a bilinear hysteretic model for simulating the response of a lead rubber bearing. In this figure, $u(t)$ is the time-varying vertical displacement of the bearing and $f_s(t)$ is the restoring force in the bearing. The stiffness of the bearing $(\mathrm{slope}{K}_v)$ is related to post-yield stiffness by α. The lead responds with elastic perfectly plastic loops; therefore, after yielding (f_s(t) > f_sy in Figure 1(b)), the stiffness is equal to the stiffness of the rubber bearing alone (αK_v) and the loops are almost bilinear. ¹⁴ $K_{\mathrm{eff}}$ is the effective stiffness value given by the ratio of the maximum restoring force ( $f_{\mathrm{smax}}$ ) to the maximum displacement ( $u_{\max }$ ). ¹⁵

Deterioration in elastomeric bearings is caused by external restraint forces which are generated in the movement of concrete in response to temperature, creep and shrinkage. The restraint may consist of friction at the bearings, bonding to already hardened concrete or by attachment to other components of the structure. Cracks resulting from the actions of external restraint forces develop in a similar manner to those caused by external load. Concern about the reliability of elastomeric bearings is increasing and studies on experimental and in-service bearings reveal the occurrence of permanent changes in the engineering properties of the bearings.^4,5 Cyclic loading is found to cause progressive damage, analogous to fatigue in metals. ¹⁶ Most damage cases are in the form of progressive debonding between the rubber and steel shims. A degree of stiffness loss can be expected due to repeated cyclic loading causing breading of molecular bonds between polymers. ¹⁷ A method to characterize and monitor the structural parameters characterizing the bearing response is clearly needed.

Many methods have been proposed for parameter identification of single degree-of-freedom (SDOF) and multiple degree-of-freedom (MDOF) structural systems. Recently, wavelet ridges of continuous wavelet transform (CWT) have been employed to identify the instantaneous frequency (IF) of a time-varying structure, where the frequency change is gradual representing a linear system with time-varying frequency. ¹⁸ Nagarajaiah and Basu ¹⁹ investigate modal identification of linear time-variant systems and they conclude that the wavelet analysis is more appropriate than the Hilbert–Huang transform (HHT) in identifying changes in the frequency of linear systems due to the number of scales and coefficients that can be applied. Feldman ²⁰ applies the Hilbert Transform (HT) to identify the non-linear instantaneous modal parameters of SDOF systems in forced and free vibration. However, the application of HT is limited as it only defines a ‘mono-component’ signal (i.e. obtain only one frequency value using a single component). ²¹ For ‘multi-component’ signals, the concept of IF becomes meaningless, given that the HT is inherently limited to ‘mono-component’ signals and breakdown of the signal components is needed. The signal needs to be defined in terms of its single ‘mono-components’ each having a different IF ⁸ as addressed by Huang et al. ²² They decompose the signal using Empirical Mode Decomposition (EMD) to obtain the Intrinsic Mode Functions (IMFs) and then use the HT to obtain the IF of the signal. This methodology is called the HHT and it is one of the most popular methods for analysis of non-linear and non-stationary data.^22–24 In the following sections, the HHT is applied to the characterization of the response of non-linear systems with hysteresis that define lead rubber bearings. First, basics on the HHT are reviewed. Second, the mathematical model of the bearing (a SDOF) is introduced. Then, the response of the model to a time-varying force and to a sudden change in the force–displacement hysteresis loop is investigated via the HHT.

The HHT

There are two main steps in the HHT analysis process:^22,25 EMD and the HT. The HT characterizes a time-domain signal $x(t)$ by extracting its envelope $a(t)$ and its phase $\theta (t)$ which can be expressed as $x\left(t\right)=a(t)\cos \mathrm{\theta }(t)$ . Huang et al. ²² apply the EMD process first to overcome the limitation of the HT to only ‘mono-component’ signals.

First step: EMD

The signal is decomposed into a number of IMFs using the ‘sifting’ process. The IMFs can be defined as a representation of oscillatory modes of variable frequency and amplitude that are time dependent, as opposed to simple harmonic functions of constant amplitude and frequency. The IMFs must satisfy two conditions:

• The number of extrema and the number of zero crossings must differ by no more than 1.

• The mean value of the envelope defined by the maxima and the minima at any given time must be 0.

The sifting process to obtain the IMFs is as follows:

1. Identify the local maxima and minima.

2. Find the mean m₁(t) of the upper and lower envelopes.

3. Obtain the difference $h_1(t)$ between the signal x(t) and the mean m₁(t), that is, $h_1(t)=x\left(t\right)-m_1(t)$ .

4. Check if $h_1$ satisfies the two conditions ((a) and (b) above) to be an IMF. If the two conditions are met, then h₁ will be the first IMF. However, this rarely occurs in the first iteration, and therefore, $h_1$ must be treated as the signal and the process (1)–(3) repeated again (i.e. $h_{11}(t)=h_1(t)-m_{11}(t)$ ), until the sifted signal satisfies the two conditions for an IMF.

After the sifting process, (1)–(3) has been done a number of times, that is, k, the two conditions will be met and the first IMF, $c_1(t)$ , will be obtained (equation (1))

$$\begin{array}{cc}\multicolumn{1}{c}{}& h_{1k}(t)=h_{1(k-1)}(t)-m_{1k}(t)\hfill \\ \multicolumn{1}{c}{}& c_1(t)=h_{1k}(t)\hfill \end{array}$$

(1)

However, too many sifting cycles could reduce all the components to a constant-amplitude signal with frequency modulation only, that is, only frequency variation is retained. To guarantee that the IMF components retain enough physical sense of amplitude and frequency modulations, the number of times, k, the sifting process repeats has to be limited. There are two main mathematical stoppage criteria proposed in the literature that are used to ensure that the sifting process stops when conditions (a) and (b) above are met:^22,23

• A first criterion based on limiting the normalized Square Difference (SD) between two successive sifting operations to a small value, that is

  $$\mathrm{SD}=\sum _{t=0}^T\left[\frac{{|h_{1(k-1)}(t)-h_{1k}(t)|}^2}{h_{1(k-1)}^2(t)}\right]\le 0.2-0.3$$

  (2)

  at time $t$ and for a total duration $T$ of the signal.

• A second criterion based on assuming the sifting process has ended after a predefined number of times. For instance, if after $S$ sifting times, there is no difference between the number of zero crossings and the extrema, or differ at most by 1, then the process is stopped.

The first criterion is difficult to implement in practice due to the need to select an accurate value of $\mathrm{SD}$ to stop the sifting process. A small value of SD does not necessarily satisfy the number of zero crossings and extrema. For that reason, the second criterion is adopted in this article, although there also exists difficulty in predefining the value of $S$ which is generally taken between 4 and 8. ²³ The value of S adopted here is 4 to save computational time given that choosing a higher value has not shown to significantly enhance the results of the analysed signals.

The process of sifting starts by separating the finest local mode $c_1(t)$ (equation (1)) from the data x(t). This gives the residual $r_1(t)=x(t)-c_1(t)$ . The residual is then used as the input signal to obtain the second IMF $(c_2(t))$ and so on. The decomposition into higher order IMFs is repeated $m$ times until the value of $r_m(t)$ becomes a ‘monotonic’ component from which no more IMFs are extracted.

Second step: the HT

Once the EMD process has been finalized, the second step of the HHT analysis is the HT of the IMFs.^24,25 A complex analytical signal is formed for each IMF consisting of a real part defined as the $i\mathrm{th}$ IMF ( $c_i(t)$ ) and an imaginary part defined by the HT ( $y_i(t)$ in equation (3))

$$y_i(t)=\frac{1}{\pi }\underset{y\to \infty }{lim}\underset{-y}{\overset{y}{\int }}\frac{c_i(\tau )}{t-\tau }d\tau$$

(3)

The complex signal can be expressed in exponential format as in equation (4)

$$z_i(t)=c_i(t)+{\mathrm{jy}}_i(t)=a_i(t)e^{j{\theta }_i(t)}$$

(4)

where $j=\sqrt{-1}$ , $a_i(t)$ is the amplitude at time $t$ and ${\theta }_i(t)$ is the instantaneous phase. In another words, the analytical signal $z_i(t)$ represents a rotation in the complex plane with the radius of rotation $a_i(t)$ and the instantaneous phase angle ${\theta }_i(t)$ , where

$${\theta }_i(t)=\arctan \left\{\frac{y_i(t)}{c_i(t)}\right\}\mathrm{and}{a}_i(t)=\sqrt{{|c_i(t)|}^2+{|y_i(t)|}^2}$$

(5)

The IF ${\mathrm{IF}}_i(t)$ associated with the ith IMF at time t can be obtained from

$${\mathrm{IF}}_i(t)=\frac{d{\theta }_i(t)}{\mathrm{dt}}$$

(6)

Using the HHT, it is possible to characterize the frequency of a non-linear system such as a non-linear support. The potential of the HHT in detecting changes in bearing stiffness is examined in the following sections.

Description of mathematical model of bearing

An elastomeric bearing is idealized as a SDOF system consisting of a mass (M) supported on a spring ( $k_v(t)$ ) and damper (C_v) system (Figure 2(a)). The response of a SDOF system with a non-linear hysteretic stress–strain relationship has been analysed by Dobson et al., ⁹ who provide a mathematical model based on Bouc–Wen capable of simulating the asymmetric hysteretic response found in bearings. The asymmetric hysteresis is composed of four different curves that describe the force–displacement response of the bearings. Here, the non-linear stiffness in the spring, K_v(t), is given by the slope of the curve in Figure 2(b), which is a smoother and more realistic representation of a bearing response than the simplified bilinear mathematical model consisting of straight lines shown in Figure 1(b). ²⁶ In Figure 2(b), z is the hysteretic displacement which defines the permanent displacement in the hysteresis, that is, at unloading the displacement does not return to the original value.

Figure 2.

Bearing model: (a) SDOF and (b) force–displacement relationship.

The response of the SDOF system in Figure 2 is governed by the dynamic equilibrium between the internal inertial ( $M\stackrel{··}{u}\left(t\right)$ ), damping ( $C_v\stackrel{·}{u}\left(t\right)$ ) and stiffness ( $k_v(t)u\left(t\right)$ ) forces and the external applied forces ( $f(t)$ ) given in equation (7)

$$M\stackrel{··}{u}\left(t\right)+C_v\stackrel{·}{u}\left(t\right)+k_v(t)u\left(t\right)=f(t)$$

(7)

where $\stackrel{··}{u}(t)$ , $\stackrel{·}{u}(t)$ and $u(t)$ are the acceleration, velocity and displacement, respectively, of the SDOF, which can be obtained by integration of equation (7) using the Wilson-θ method. The response of the non-linear SDOF system in equation (7) can be re-formulated ²⁷ as

$$\stackrel{··}{u}\left(t\right)+2\xi \omega (t)\stackrel{·}{u}\left(t\right)+f_s(t)=f\left(t\right)/M$$

(8)

where $\xi$ is damping, $\omega (t)$ is the natural frequency given by $\omega (t)=\sqrt{k_v(t)/M}$ , and $f_s(t)$ is the restoring force for the non-linear system given by

$$f_s(t)=\frac{K_v(t)u(t)}{M}=\alpha {\omega }^2(t)u(t)+(1-\alpha ){\omega }^2(t)z(t)$$

(9)

In equation (9), $\alpha$ is the ratio of post-yield to pre-yield stiffness and $z(t)$ is the hysteretic displacement (i.e. the difference between the displacement at loading and unloading (permanent displacement)). The post-yield to pre-yield stiffness ratio ( $\alpha$ ) regulates the shape of the hysteresis. If α = 1, the force–displacement relationship reverts to linear case, while for α = 0.2 ( $0\le \alpha \le 1$ ), a hysteresis force–displacement relationship, as illustrated in Figure 2(b), is formed.

The Bouc–Wen hysteretic displacement ( $z(t)$ , which is given by the horizontal axis in Figure 2(b)) is used to guarantee a smooth transition between pre-yield and post-yield stiffness that represents the real force-deformation behaviour of the bearing. $z(t)$ is governed by equations (10) and (11) where the influences of the characteristic strength, damping, post- to pre-yield stiffness ratio and other parameters are taken into account via pre-established coefficients ²⁷

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{z}\left(t\right)=A\stackrel{·}{u}\left(t\right)-\beta \left|\stackrel{·}{u}\left(t\right)\right|{\left|z(t-\Delta t)\right|}^{n-1}z(t-\Delta t)}\end{array}-\gamma \stackrel{·}{u}\left(t\right){\left|z(t-\Delta t)\right|}^n$$

(10)

$$z(t)=\stackrel{·}{z}\left(t\right)+z(t-\Delta t)$$

(11)

where $\stackrel{·}{z}(t)$ is the incremental hysteretic displacement; $\stackrel{·}{u}(t)$ is the velocity of the SDOF system; $A$ , $n$ , $\beta$ and $\gamma$ are the coefficients that define the shape of the hysteresis ( $A$ , $\beta$ and $\gamma$ influence the loop size and $n$ influences its smoothness); $\Delta t$ is the incremental time step and $z(t-\Delta t)$ is the hysteresis displacement at the previous time step. In general, by using equation (10), the coefficients are dependent on each other to produce a hysteretic loop that is representative of the hysteresis bearing response. In order to have more flexibility and better control over the coefficients that define the hysteresis shape, equation (12) is used instead of equation (11) to obtain the force–displacement loops presented in Figure 2(b) as described by Dobson et al. ⁹

$$\begin{array}{cc}\multicolumn{1}{c}{}& \stackrel{·}{z}(t)=A\stackrel{·}{u}\left(t\right)+B(\stackrel{·}{u}\left(t\right)+\left|\stackrel{·}{u}\left(t\right)\right|){(z(t-\Delta t)+|\stackrel{·}{z}(t)|)}^n\hfill \\ \multicolumn{1}{c}{}& +C(\stackrel{·}{u}\left(t\right)+\left|\stackrel{·}{u}\left(t\right)\right|){(z(t-\Delta t)-|\stackrel{·}{z}(t)|)}^n\hfill \\ \multicolumn{1}{c}{}& +D(\stackrel{·}{u}(t)+|\stackrel{·}{u}(t)|){(z(t-\Delta t)+|\stackrel{·}{z}(t)|)}^n\hfill \\ \multicolumn{1}{c}{}& +E(\stackrel{·}{u}(t)-|\stackrel{·}{u}(t)|){(z(t-\Delta t)-|\stackrel{·}{z}(t)|)}^n\hfill \end{array}$$

(12)

Equation (12) decouples the hysteretic loop in four loading regimes each of which being regulated by a different constant shape coefficient: $B$ , $C$ , $D$ or $E$ (shown in Figure 2(b)). The constants A and n define the shape and level of nonlinearity, respectively, in the restoring force-deformation behaviour, and the constants B, C, D and E shape the hysteresis loop. The constant shape parameter to employ will depend on the location of the force–displacement curve and the signs of the velocity ( $\stackrel{·}{u}(t)$ ) and the hysteretic displacement ( $z(t)$ ). Three of the non-linear terms in equation (12) will perish for each of the four possible combinations of signs for $\stackrel{·}{u}(t)$ and $z(t)$ . For example, C( $\stackrel{·}{u}(t)>0,z(t)<0$ ) is considered if the velocity is positive and the hysteretic displacement is negative, while $B$ , $D$ and $E$ become zero. D( $\stackrel{·}{u}(t)<0,z(t)>0$ ) is considered if the velocity is negative and the hysteretic displacement is positive (the opposite of the condition for C) and then, B, C and E are zero. Similar for the regions B( $\stackrel{·}{u}(t)>0,z(t)>0$ ) and E( $\stackrel{·}{u}(t)<0,z(t)<0$ ) illustrated in Figure 2(b). When $\stackrel{·}{u}(t)=0$ and $z(t)=0$ , the system is linear and there are no constants that need to be defined for the force–displacement curve.

The main advantages of the model in equation (12) over the Bouc–Wen model in equation (11) is the added flexibility to the hysteretic cycle as each curve is defined separately. Setting the non-linear parameters B, C, D and E equal to 0 results in a linear relationship with a slope equal to that of the initial tangent slope. Depending on the selected coefficients, choosing a positive value shall bow the curve inwards (or outwards) from the loop, with the negative value producing the opposite effect.

Characterization of non-linear SDOF bearing model using the HHT

The response of the SDOF in Figure 2 ( $u(t)$ ) is obtained for a sinusoidal load $f(t)$ given in Newton by equation (13) for testing the ability of the HHT-based method in characterizing the structural system

$$f\left(t\right)={10}^6\sin (10t)$$

(13)

In this theoretical test, the SDOF bearing model has a mass of M = 3135.8 kg and a viscous damping of 6% (ξ = 0.06).^5,28,29 A vertical effective stiffness of 1.14375 × 10⁹ N/m is assumed which leads to an ‘effective’ natural frequency of 30.7 Hz. The following values are assigned to the constants of the hysteresis model of equation (12) to define the shape of Figure 2(b): $A=1$ , $\alpha =0.2$ , $n=1$ , $B=-500$ , $C=-500$ , $D=500$ and $E=500$ . These values are selected to create a force–displacement curve that resembles previous experimental investigations on the response of hysteretic bearings.^30–32

Figures 3(a) and (b) show the acceleration ( $\stackrel{··}{u}(t)$ ) and displacement ( $u(t)$ ), respectively, of the SDOF system under the sinusoidal load specified in equation (13). The four stages (B, C, D and E) of the hysteresis in Figure 2(b) are shown by rounded dashed circles in Figure 3 for the response between 1.1 and 1.73 s (which covers one full hysteresis loop 0.63-s long, that is, the time from start of loading at B to unloading at D and E to reloading at C and back to B in the curve). When the curve is at C, the stiffness is relatively high and the acceleration response (from 1.1 to 1.225 s) has lower amplitude, which is similar to the stage when the curve is at D (from 1.425 to 1.54 s). When the curve is at B (from 1.225 to 1.425 s), the lower stiffness produces higher amplitude and a relatively linear force–displacement relationship (Figure 2(b)), similar to stage E from 1.54 to 1.73 s.

Figure 3.

Response of the non-linear SDOF system to a sinusoidal load: (a) acceleration and (b) displacement.

The HHT is applied to the acceleration signal from Figure 3(a). Figure 4(a) and (b) show the first four IMFs resulting from the EMD process and the associated IFs via the HT, respectively. Although the total number of IMFs extracted from the acceleration is seven, only the first four IMFs are shown in the figure. The last three IMFs (IMF5–7) are residues of small magnitude that do not provide meaningful information on the system behaviour. Lower IMFs contain higher frequency components of the signal, and IMF1 appears to have higher coefficient values than other IMFs and to dominate the response.

Figure 4.

HHT: (a) IMFs 1–4 using EMD; (b) associated instantaneous frequencies (IF₁(t)–IF₄(t)) in Hertz using the HT.

The frequency of the applied sinusoidal load (equation (13)) is 1.59 Hz and it is captured by IMF4 in Figure 4. The natural frequency of the SDOF system, NF(t), is time varying and determined in Hertz by

$$\mathrm{NF}(t)=\frac{1}{2\pi }\sqrt{\frac{k_v(t)}{M}}$$

(14)

where the stiffness $k_v(t)$ is obtained from the force–displacement curve in Figure 2(b) at each point in time $t$ . Figure 5 shows the close correlation between the IF₁(t) extracted from the IMF1 and the natural frequency of the system, NF(t). Discrepancies between IF₁(t) and NF(t) are noticed for the lowest stiffness in the loop (stages B and E), that is, IF₁(t) predicts about 12 Hz for a minimum NF of 14 Hz. The latter can be partially attributed to mode mixing in the EMD process.

Figure 5.

Comparison of instantaneous frequency of IMF1 using HT and actual instantaneous natural frequency (NF(t)) of the system.

Figure 6 shows the variation of stiffness with time of the SDOF system using $K_v(t)$ calculated using equation (14) with the predicted IF₁(t) in Figure 5, and the ‘actual $K_v$ ’ obtained from the force–displacement curve at each point in time. Figure 6 shows how the variations of stiffness due to frequency changes are captured with a reasonable degree of accuracy for stages B, C, D and E of a full hysteresis cycle. The high-frequency content ( $\cong 41-16\mathrm{Hz}$ ) in Figure 5 from 1.1 to 1.225 s corresponds to the points in the left portion of the C stage (associated with high slopes and stiffness in Figures 2(b) and 6, respectively). Stiffness is significantly lower when entering stage B from 1.225 to 1.425 s where the frequency range is $\cong 16-14\mathrm{Hz}$ . Stiffness is again higher in stage D from 1.425 to 1.54 s where the frequency range increases to $\cong 41-16\mathrm{Hz}$ . Finally, stage E takes place between 1.54 and 1.73 s ( $\cong 16-14\mathrm{Hz}$ ) where the stiffness reduces to similar levels found before in stage B.

Figure 6.

Comparison of stiffness prediction ( $K_v$ ) based on IF by HHT and true stiffness (‘Actual $K_v$ ’) of the system.

Although the IF₁ and $K_v$ using HHT do not exactly match those of the ‘actual values’ of the system, they closely approximate the exact frequency and stiffness of the system. Figures 7 and 8 show the analysis of the acceleration in Figure 3 using Fast Fourier Transform and CWT, respectively, to compare with the HHT. In Figure 7, the Power Spectral Density (PSD) shows the range of possible frequencies of the system with the highest peaks taking place at 1.6 and 4.8 Hz. These two large peaks are followed by others of smaller magnitude up to 42 Hz representing the continuous change in frequency of the system with time. It is evident that the PSD is very limited in characterizing a non-linear system due to the lack of time information. There is no way of knowing if these frequencies correspond to stages B, C, D or E on the force–displacement hysteresis curve.

Figure 7.

PSD of acceleration in Figure 3(a).

Figure 8.

Wavelet analysis of acceleration in Figure 3(a): (a) contour plot of wavelet coefficients versus frequency and time, (b) wavelet coefficients versus frequency for times of 1.1, 1.225, 1.425, 1.54 and 1.73 s, and (c) wavelet coefficients versus time for frequencies of 41 and 14 Hz.

The CWT uses inner products to measure the similarity between a signal and a function known as the mother wavelet. In the CWT, the signal is compared to shifted and compressed or stretched (scaled) versions of the mother wavelet. By comparing the signal to the wavelet at various scales and points in time, a frequency–time map of wavelet coefficients is obtained. ³³ Figure 8(a) shows the wavelet coefficients versus pseudo-frequency and time that result from applying the CWT to the acceleration response of Figure 3(a) using the Mexican Hat mother wavelet. The pseudo-frequency is directly proportional to the centre frequency of the wavelet in Hertz (for the Mexican Hat wavelet, this centre frequency is 0.25 Hz) and inversely proportional to the scale of the wavelet and to the sampling period. It is difficult to interpret the behaviour of the system from Figure 8 alone, so sections at times 1.1, 1.225, 1.425, 1.54 and 1.73 s (A-A, B-B, C-C, D-D and E-E, respectively) are plotted versus frequency in Figure 8(b) and sections at frequencies 14 and 41 Hz (F-F and G-G, respectively) are plotted versus time in Figure 8(c). In Figure 8(b), section A-A corresponds to the initial stage of C, while section B-B corresponds to the end of stage C and beginning of stage B. Sections through C-C and D-D correspond to the end of stage B and beginning of stage D. It can be seen that the sections C-C and D-D have the same magnitude than sections A-A and B-B, but they are in the negative region. Section E-E corresponds to the end of stage E and it is almost equal in magnitude to section A-A suggesting that the stiffness is returned to its original value and a full hysteresis cycle is complete. Figure 8(c) shows the variation in content of the pseudo-frequencies 14 and 41 Hz (lower and upper values respectively of the range of variation of the frequency of the system) with time. However, the range of variation of the frequency of the system, the frequency at which the system is responding for each point in time, and its duration are not as clearly defined by wavelets as in Figure 5 simply using the IF₁(t) by HHT.

It has been shown that HHT can be used to characterize frequency and stiffness changes due to nonlinearity properties of a SDOF system, although one problem that may arise with the EMD process is over-decomposition of the signal resulting into excessive IMFs and mode mixing, that is, IMFs sharing similar frequencies. It is also important to note that this is a theoretical investigation of a SDOF system with a specific hysteretic loop and that experimental results may contain other complexities involving external factors such as temperature and/or inconsistencies of the applied load and the specimen and/or internal factors such as friction between steel shims and rubber (for steel-reinforced rubber bearings).

Identification of changes in the bearing behaviour using the HHT

Although the HHT was initially developed to characterize non-linear systems, its application to structural damage detection has been limited mostly to sudden losses in stiffness using simple periodic functions, experimental beams with crack or linear time-variant systems.^19,34–37 The reason for this possibly lays on the lack of data to characterize how damage affects the overall performance of non-linear structural systems. In this investigation, damage is theoretically simulated through the introduction of a reduction in the effective stiffness of a hysteretic loop as shown in Figure 9(a). ³⁸ The acceleration response of the SDOF system is illustrated in Figure 9(b) where a 10% reduction in effective stiffness is introduced after 2.5 s. This change in the effective stiffness of the hysteretic loop is visualized as an increase in the maximum amplitude of the acceleration signal of 0.282 m/s² with respect to the original acceleration (i.e. maximum acceleration is 28.8% higher after the stiffness reduction). This large increase in acceleration is attributed to the low stiffness values experienced by the bearing at loading and unloading. In comparison, at higher stiffness, the change in acceleration is much lower (6.2%).

Figure 9.

(a) Force–displacement relationship for the SDOF system before and after 10% loss in effective stiffness and (b) acceleration response before and after 10% loss in effective stiffness.

Four IMFs of interest can be obtained from Figure 9(b) using EMD, and the IFs of these IMFs are shown in Figure 10(a). Edge effects (before 0.5 s and after 4.5 s) appear in the IFs as a result of the decomposition process. ²² IF₁(t)–IF₄(t) show a sudden peak at the instant of change in effective stiffness. This peak can be clearly visualized in IF₃(t) and IF₄(t) of Figure 10(a), due to the relative high-frequency content at this location with respect to the rest of the IMF. The dominant IMF is IMF1 and Figure 10(b) zooms in on the IF₁(t) shown in Figure 10(a). The average range of IF₁(t) is 42.35–10.64 Hz prior to the simulated damage and 38.35–9.46 Hz after damage has occurred. Once edge effects are ignored, the reduction in frequency as a result of the stiffness change is 9.6% (−4 Hz) at the upper range (cases C and D in Figure 2(b)) and 11.1% (−1.18 Hz) at the lower range (cases B and E in Figure 2(b)). The average drop in stiffness is (9.6% + 11.1%)/2 = 10.35%, that is, approximately the 10% reduction in effective stiffness.

Figure 10.

Instantaneous frequencies in Hertz before and after 10% effective stiffness loss of (a) IMFs 1–4 and (b) IMF1.

This section has demonstrated the capability of the HHT to extract important time–frequency information about the non-linear behaviour of a SDOF system. IF₁(t) has been able to capture a small change in effective stiffness of 10%. It is acknowledged that experimental setups and more complex systems will require IFs from more than one IMF to gather sufficient information to characterize the system.

Influence of noise on the performance of the HHT

Noise can be caused by a number of factors which include the natural intermittent instabilities in the system, the concurrent phenomena in the environment where investigations are conducted, and/or the sensors and recording systems. As a result, data are always an amalgamation of signal and noise, that here they are simulated using the additive model

$$\left\{{\stackrel{··}{u}}_n\left(t\right)\right\}=\left\{\stackrel{··}{u}\left(t\right)\right\}+n_l\{N\}\sigma$$

(15)

where $\{{\stackrel{··}{u}}_n(t)\}$ is the acceleration corrupted by noise, $\{\stackrel{··}{u}(t)\}$ is the acceleration of the system at time $t$ , $n_l$ is the noise level, $\{N\}$ is the standard normal distribution vector with zero mean and unit standard deviation, and $\sigma$ is the standard deviation of the noise-free signal. For testing purposes, the original acceleration data of Figure 3(a) are corrupted with 1% noise level (i.e. $n_l=0.01$ in equation (15)) which results into Figure 11.

Figure 11.

Acceleration response of the non-linear SDOF system to a sinusoidal load with 1% noise.

Following EMD of the corrupted acceleration signal, the corresponding IMFs 1–8 are represented in Figure 12(a) and (b). Compared to the noise-free signal (Figure 4), there is a higher number of IMFs that need to be considered as a result of the presence of noise. It can be seen that even low noise levels can have a significant impact on the decomposition process. Figure 12(c) and (d) shows the IFs associated with each IMF. IMFs 1–3 reveal high variability of frequency content between 0–500, 0–300 and 0–150 Hz for IMF1, IMF2 and IMF3, respectively, which can be attributed to noise. IMFs 4–8 contain a smoother variation of frequencies that appears to be in the range of IFs obtained in Figure 4(b) for noise-free acceleration.

Figure 12.

HHT for signal with 1% noise: (a) IMFs 1–4 using EMD. (b) IMFs 5–8 using EMD. (c) Associated instantaneous frequencies (IF₁(t)–IF₄(t)) in Hertz using the HT. (d) Associated instantaneous frequencies (IF₅(t)–IF₈(t)) in Hertz using the HT.

However, when plotting the IFs of IMF4 and IMF5 (Figure 12), together with the time-varying natural frequency of the system in Figure 13, the influence of noise is clear. The extraction of a clean IF that matches the frequency of the system is rendered unsatisfactory.

Figure 13.

Comparison of the instantaneous frequency of IMF using HT on a 1% noisy signal and actual instantaneous natural frequency (NF(t)) of the system: (a) instantaneous frequency of IMF4 and NF(t) and (b) instantaneous frequency of IMF5 and NF(t).

For cases where noise has distinct time or frequency scales from those of the true signal, Fourier transform filters can be applied to separate the noise from the signal. Once the signal is filtered from noise, the HT has been shown to be capable of extracting phases and IFs for damage identification of linear scenarios.^39,40 Similarly, sections ‘Characterization of non-linear SDOF bearing model using the HHT’ and ‘Identification of changes in the bearing behaviour using the HHT’ have demonstrated that the HHT is a valid tool to characterize instantaneous changes using the response of a non-linear SDOF bearing system where noise is non-existing or has been safely removed. However, traditional filtering methods may not be as reliable in separating mixed noise and harmonics in non-linear systems, and section ‘Influence of noise on the performance of the HHT’ has shown that the direct application of EMD and the HT to a corrupted signal can produce misleading results. In this regard, Ensemble Empirical Mode Decomposition (EEMD) ⁴¹ is an alternative to EMD that attenuates mode mixing and noise to some extent.

Conclusion

A hysteretic Bouc–Wen model has been implemented in a SDOF system to simulate the response of a lead rubber bearing. An IF analysis of the response of the SDOF via the HHT has provided information that allows characterizing the different stages of stiffness in this non-linear system (i.e. stages $B$ , $C$ , $D$ and $E$ of the hysteretic loop). The HHT has overcome the inability of the PSD to distinguish the multiple frequencies associated with one mode of vibration in a non-linear system. It has also outperformed wavelet analysis in characterizing the system. A sudden change in effective stiffness has been successfully identified via the IF. Finally, it has been shown that the accuracy of the method can be severely limited by noise. Further research is needed to address the impact of noise and of integrating the bearing as part of a more sophisticated MDOF system on the performance of the HHT.