Nonlinearity identification by time-domain-only signal processing

International Journal of Non-Linear Mechanics 54 (2013) 85–98

Contents lists available at SciVerse ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

Nonlinearity identification by time-domain-only signal processing P. Frank Pai a,n, Bao-Anh Nguyen a, Mannur J. Sundaresan b a b

Department of Mechanical and Aerospace Engineering, University of Missouri, Columbia, MO 65211, United States Department of Mechanical Engineering, North Carolina A&T State University, Greensboro, NC 27411, United States

art ic l e i nf o

a b s t r a c t

Article history: Received 26 June 2011 Received in revised form 30 March 2013 Accepted 2 April 2013 Available online 10 April 2013

A conjugate-pair decomposition (CPD) method is proposed for signal decomposition, dynamics characterization, and nonlinearity identification all in the time domain only. CPD uses the empirical mode decomposition method with signal conditioning techniques to decompose a compound signal into well separated intrinsic mode functions (IMFs) and then uses a pair of sliding conjugate functions to accurately extract the time-varying frequency and amplitude of each IMF using only three neighboring data points for each time instant. Because the variations of frequencies and amplitudes of IMFs contain system characteristics, they can be used for dynamics characterization and nonlinearity identification of discrete and continuum systems. Several discrete and continuum systems' simulated and experimental responses are used to validate CPD's accuracy and capability for system and nonlinearity identification. Experimental nonlinear free vibration of a horizontally cantilevered steel beam subject to an initial tip displacement is analyzed using CPD, and direct numerical simulations using a fully nonlinear finiteelement code are also performed and compared. Both experimental and numerical results show that the first-mode vibration contains a softening cubic nonlinearity, and the shortening-induced longitudinal inertia and nonlinear modal coupling of the first few modes make the effective modal damping of the first mode nonlinear. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Time-domain signal processing Time-frequency analysis Signal conditioning Time-varying frequency and amplitude Nonlinearity identification

1. Introduction Every mechanical system has its own resonant frequencies, and its response to an excitation is primarily determined by the closeness and commensurability (if a nonlinear system with polynomial types of nonlinearity) between the excitation frequency and the system's resonant frequencies [1,2]. A signal is mainly described by its frequency and amplitude. Hence, extraction of a signal's frequency and amplitude is the major task in signal processing for system identification, structural damage detection and on-line health monitoring, experimental verification and updating of structural models, aircraft flight tests, design of effective control schemes, and many other engineering applications [3–8]. A linear second-order system's free vibration is well described by a regular harmonic function with a constant frequency, and hence harmonic functions are often chosen as basis functions for signal decomposition of dynamical system response. However, a complex system's characteristic responses may not be harmonics but are physically meaningful basis functions for analyzing the system's responses to other compound/complicated excitations.

n

Corresponding author. Tel.: +1 573 884 1474; fax: +1 573 884 5090. E-mail address: [email protected] (P. Frank Pai).

0020-7462/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijnonlinmec.2013.04.002

For example, a nonlinear single-degree-of-freedom oscillator's free vibration is a harmonic with amplitude and frequency modulations and is called a distorted harmonic in this work. It is more meaningful to treat this distorted harmonic as a basis function, instead of decomposing it into several regular harmonics. For a nonlinear continuous or multiple-degree-of-freedom system subject to a general excitation, its response consists of multiple coupled basis functions. How to separate these basis functions and how to identify different system parameters and nonlinearities from the basis functions are challenging but very desirable for many engineering applications [1–3]. Nonlinearities can be due to large elastic deformations (i.e., geometric nonlinearity), deformation-dependent material properties (i.e., material nonlinearity), backlash between parts, clearances between mounting brackets, geometric constraints on deformation, misalignment of substructures, dry friction, and many types of nonlinear hysteretic damping (e.g., aerodynamic damping, damping of shock absorbers for vehicles, and material damping of shape memory alloys and other materials) [1,2]. Notorious nonlinear vibrations include wing and panel flutters, tail buffeting, limit-cycle oscillation of wings with attached weapons, and amplitude-dependent natural frequencies of pylon-store-wing assemblies and tail-fin assemblies in aircraft engineering; lowfrequency large-amplitude vibrations of solar panels of satellites, membrane wrinkling of space structures, complex deployment

86

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

dynamics of deployable/inflatable structures, multiple possible deployed geometries of large scientific balloons, and odd opening modes of parachutes for landing space vehicles in inner- and outer-space engineering [8]; large-amplitude vibrations of hightower buildings under gusty wind, flutter of cable bridges, and whirling of power and communication cables in civil engineering; and dynamics of leaf springs and shock absorbers of cars, brake squeal of cars, dynamics of microbeams in MEMS devices, and chattering of machining processes in mechanical engineering. Damage introduces extra nonlinearities and aging aggravates nonlinearities [2–5]. For example, dry friction between two fractured surfaces is a nonlinear effect, and the breathing of a crack causes nonlinear intermittent transient response. Also, demountable/retractable civil structures are built with nonlinearities caused by loose joints, and the usage/aging aggravates such nonlinearities. For analysis of dynamical systems, signal processing techniques are mainly developed for processing time-domain signals to filter noise, separate independent basis components, extract frequencies and amplitudes of the separated basis components, reveal linearity/nonlinearity, and/or capture transient events for dynamics characterization, system identification, damage detection, and/or health monitoring. Most of signal processing methods for dynamics characterization require both time- and frequencydomain signal processing techniques. Transforming time-domain data into the frequency domain introduces errors and/or numerical noise, and some methods even require the transformed frequency-domain data be transformed back into the time domain, which introduces even more errors. For transformation from the time domain to the frequency domain, the first step is to choose a complete set of basis functions, and the second step is to perform convolution computation between the signal under analysis and the chosen basis functions to extract components similar to the basis functions [9]. Unfortunately, this approach is not suitable for a nonlinear and/or nonstationary signal because its frequency and amplitude change with time. Hence, adaptive basis functions need to be used, and they can only be derived from the signal itself. Fourier transform (FT) is often used to transform time-domain data into frequency-domain data in order to identify a signal's characteristic frequencies and amplitudes from the obtained spectrum [9], but a spectrum cannot show the time-varying characteristics of the processed signal. To overcome this problem, short-time Fourier transform (STFT), wavelet transform (WT), and Hilbert– Huang transform (HHT) have been developed for time–frequency analysis [10–13]. STFT uses windowed regular harmonics and function conformality to simultaneously extract time-localized regular harmonics, and WT uses scaled and shifted wavelets and function conformality to simultaneously extract time-localized components similar to the wavelets [12]. Unfortunately, STFT and WT cannot accurately extract time-varying frequencies and amplitudes because of the use of pre-determined basis functions and function conformality for simultaneous extraction of components. On the other hand, HHT combines the use of empirical mode decomposition (EMD) and Hilbert transform (HT) for time–frequency decomposition. Because HHT does not use predetermined basis functions and function conformality for component extraction, it provides more concise extraction of components and more accurate time-varying amplitudes and frequencies of extracted components. EMD uses the apparent time scales revealed by the signal's local maxima and minima to sequentially sift components of different time scales, starting from high-frequency to lowfrequency ones [13–15]. Unfortunately, when a high-frequency component decays, EMD may extract the next high-frequency wave into the same component. This causes intermittence in the extracted components and difficulties in calculation of timevarying frequencies and amplitudes. Moreover, HHT uses Hilbert

transform in order to compute the time-varying frequency and amplitude of a component extracted by EMD, but Hilbert transform aggravates the discontinuity-induced Gibbs' effect introduced by FT at the two data ends when it transforms the frequency-domain data back into the time domain, as explained later in Section 2.2. In this work we present a conjugate-pair decomposition (CPD) method that processes only time-domain data for accurate time– frequency decomposition of signals, improves the sifting capability of EMD, and alleviates Gibbs' effect at the two data ends. Moreover, we present a method of using perturbation solutions and results from time–frequency analysis for system and nonlinearity identification. The capability and accuracy of CPD for system and nonlinearity identification is verified by using numerical simulations and experimental data.

2. Frequency- and time-domain signal processing methods As shown later in Section 3, a harmonic with smoothly and periodically varying amplitude and frequency is a physically meaningful and independent basis function of nonlinear dynamical systems, and it is called in this work a distorted harmonic. For a distorted harmonic uðtÞ ¼ aðtÞcos θðtÞ, if its conjugate part vðtÞ≡aðtÞcosðθðtÞ−π=2Þ ¼ aðtÞsin θðtÞ is known, its time-varying amplitude aðtÞ, phase angle θðtÞ, and frequency ωðtÞ and the corresponding analytic function zðtÞ can be defined as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tan−1 v _ zðtÞ≡u þ jv ¼ aejθ ; ω≡θ; u2 þ v2 ; θ≡ ð1Þ u pffiffiffiffiffiffi _ where θ≡dθ=dt and j≡ −1. Unfortunately, two unknown variables (i.e., θðtÞ and aðtÞ) cannot be uniquely determined from just one known function (i.e., uðtÞ). Hence, vðtÞ needs to be obtained from signal processing on uðtÞ. If uðtÞð ¼ acos θÞ is measured and _ _ uðtÞð ¼ a_ cos θ−aθsin θÞ is obtained by numerical differentiation or both are experimentally and simultaneously measured by using displacement and velocity sensors, it is still impossible to use the two known functions to calculate aðtÞ and θðtÞ because there are _ four unknowns aðtÞ, θðtÞ, a_ ðtÞ, and θðtÞ in this case. For a compound signal consisting of multiple independent basis functions, signal processing is even more difficult because basis functions need to be separated before extraction of each function's time-varying amplitude and frequency. Hence, signal processing techniques often transform a time-domain signal into the frequency domain in order to separate basis functions and obtain amplitudes and frequencies of separated basis functions. Unfortunately, transformation between different domains often introduces numerical errors. Next we compare signal processing methods that process data in different domains. a≡

2.1. Frequency-domain methods If many sampled values of uðtÞ are available and the sampled signal duration T is larger than the fundamental period of uðtÞ, the harmonic components and their amplitudes of uðtÞ can be extracted using the discrete Fourier transform (DFT) as [9] N=2

uðt k Þ≡uk ¼ a0 þ 2∑i ¼ 1 ðai cosωi t k þ bi sinωi t k Þ   N=2 ¼ a0 þ Real 2∑i ¼ 1 U i ejωi tk 1 N 1 ∑ u cosωi t k ; bi ¼ ∑N u sinωi t k ; N k¼1 k N k¼1 k 1 N U i ≡ai −jbi ¼ ∑k ¼ 1 uk e−jωi t k N

ð2aÞ

ai ¼

ð2bÞ

where ωi ≡2πi=T, t k ≡ðk−1ÞΔt, k ¼ 1; …; N, N is the total number of samples (assumed to be even here), Δt is the sampling interval, 1=Δt is the sampling frequency, Tð ¼ NΔtÞ is the sampled signal

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

duration, Δf ð ¼ 1=TÞ is the frequency resolution, the maximum (or Nyquist) frequency is 1=ð2ΔtÞð ¼ N=ð2TÞÞ, U i is the Fourier spectrum, and ai and bi are constant amplitudes of the extracted regular harmonics. Eq. (2b) shows that the conformality (or orthogonality) between uðtÞ and cosωi t and sinωi t is used to calculate the timeaveraged amplitude of each regular harmonic component. This is a brute-force method because regular harmonics of uniformly spaced frequencies are chosen as basis functions and uðtÞ is expressed as a summation of all N=2 regular harmonics of constant amplitudes and frequencies. Eq. (2a) shows that uðtÞ is assumed to be a periodic function having a period T although the original function may not be periodic. Hence, the period T needs to cover at least one fundamental period of the signal in order to identify the signal's lowest-frequency component from the spectrum. Moreover, because uNþ1 ¼ u1 is implicitly assumed in DFT, if sudden change (i.e., discontinuity) exists between uN and u1 due to sampling and/or transient effects, the ω1 harmonic is not a real component. Furthermore, in order to account for this discontinuity, many high-frequency harmonics exist and they cause the leakage problem in the frequency domain and Gibbs' effect (also called the edge effect) in the time domain at the two data ends. In order to reduce the leakage and Gibbs' effect, the sampled signal duration T needs to cover at least three cycles of the lowestfrequency component to be extracted [16]. Moreover, Shannon's sampling theory requires more than two samples per period of the component to be extracted in order to capture the component's frequency. Hence, more than 6 samples covering 3 periods are needed in order to extract the correct but averaged frequency. More seriously, a small T results in a low frequency resolution Δf ð ¼ 1=TÞ, and a small N results in a small maximum frequency ð ¼ N=ð2TÞÞ that can be revealed. Hence, T and N are often numerically increased by zero-filling in using the short-time Fourier transform (STFT) to capture time-varying amplitudes and frequencies. Unfortunately, different lengths of zero-filling result in different values for the sampled signal duration T and hence different ωi ð ¼ 2πi=TÞ for the basis functions in Eq. (2a). To match with the zero-valued section of the modified function and the discontinuity due to zero-filling, additional artificial harmonics are often introduced into Eq. (2a). Hence, this frequency-domain method is not able to provide accurate time-varying frequency and amplitude.

2.2. Dual-domain methods Hilbert transform (HT) is a typical example of dual-domain signal processing methods [17], and it is capable of accurate extraction of time-varying frequency and amplitude of a signal without any moving average. The conjugate part vðt k Þ (i.e., the Hilbert transform of u) of the uðt k Þ in Eq. (2a) can be obtained by shifting the phase of each harmonic by π=2 as −π  N=2 þ 2∑i ¼ 1 ½ai cosðωi t k −π=2Þ þ bi sinðωi t k −π=2Þ vðt k Þ≡vk ¼ a0 cos 2   N=2 N=2 ð3Þ ¼ 2∑i ¼ 1 ðai sinωi t k −bi cosωi t k Þ ¼ Imag 2∑i ¼ 1 U i ejωi t k In other words, the spectral coefficients ai and bi of uðtÞ can be used to compute vðtÞ using the inverse discrete Fourier transform (IDFT) from the frequency domain back into the time domain. After v is obtained from Eq. (3) and if a0 ¼ 0, the instantaneous amplitude a and phase angle θ of uðt k Þ can be obtained using Eqs. (2a) and (3) as a¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v N=2 u2k þ v2k ; θ ¼ tan−1 k ; uk þ jvk ¼ 2∑i ¼ 1 U i ejωi t k ¼ aejθ uk

ð4Þ

87

Moreover, it follows from Eq. (4) that N=2 u_ k þ jv_ k ¼ 2∑i ¼ 1 jωi U i ejωi tk

ð5Þ

In other words, the spectral coefficients ai and bi of uðtÞ can also be used to compute the instantaneous velocity using IDFT, instead of using finite difference in the time domain. After vk ; u_ k ; and v_ k are obtained using Eqs. (3) and (5), the instantaneous frequency ω at t ¼ t k can be computed as ω≡

dθ dðtan−1 vk =uk Þ u v_ −u_ v ¼ ¼ k k 2 k k dt dt a

ð6Þ

If significant noise exists, the ω at t ¼ t k can be obtained by averaging over 2pΔt to reduce the influence of noise as ω≈

∑pi¼ −pþ1 ½θðt k þ iΔtÞ−θðt k þ ði−1ÞΔtÞ ð2pΔtÞ

ð7Þ

Unfortunately, this dual-domain method is essentially based on DFT and can extract a signal's time-varying frequency and amplitude only by processing a large length of sampled data and only if the signal has no moving average. More seriously, the two transformations DFT and IDFT between time and frequency domains worsen the discontinuity-induced edge effect in DFT at the two data ends, and hence the extracted ωðtÞ and aðtÞ are inaccurate at ends. To avoid this edge effect caused by DFT, timedomain-only signal processing is a better approach, as shown next in Sections 2.3 and 2.4. 2.3. Time-domain-only signal processing method A signal may consist of multiple independent regular or distorted harmonics and to separate and characterize them in only the time domain is challenging but possible. Here we propose a conjugate-pair decomposition (CPD) method for time-domainonly data-driven signal processing by using the empirical mode decomposition (EMD) with signal conditioning techniques and a 3-point frequency extraction method. EMD is a data-driven signal decomposition technique that sequentially extracts regular or distorted zero-mean harmonics from a signal, starting from high- to low-frequency components. EMD is essentially different from Fourier and wavelet decomposition methods because it does not use pre-determined basis functions and the conformality between the signal and the chosen basis functions to extract components [13–15]. The EMD decomposes a signal uðtÞ into n intrinsic mode functions (IMFs) ci ðtÞ and a residual r n as [14,15] uðtÞ ¼ ∑ni¼ 1 ci ðtÞ þ r n ðtÞ

ð8Þ

where c1 is the first extracted IMF and has the shortest characteristic time scale. The characteristic time scale of c1 is defined by the time lapse between the local extrema of u. Once the local extrema are identified, compute the upper envelope by connecting all the local maxima using a natural cubic spline, compute the lower envelope by connecting all the local minima using another natural cubic spline, subtract the mean of the upper and lower envelopes, m11 , from the signal, and then treat the residuary signal as a new signal. Repeat these steps for K times until the left signal has a pair of symmetric envelopes (i.e., m1K ≈0), and then define c1 as c1 ≡u−m11 ⋯−m1K

ð9Þ

This sifting process eliminates low-frequency riding waves, makes the wave profile symmetric, and separates the highestfrequency IMF from the current residuary signal. During the sifting process for an IMF a deviation Dv is computed from any two consecutive sifting results as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 N 2 Dv ≡ ∑N ð10Þ i ¼ 1 ½c1k ðt i Þ−c1k−1 ðt i Þ =∑i ¼ 1 c1k−1 ðt i Þ

88

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

where, for example, c1k ≡u−m11 ⋯−m1k , t i ¼ ði−1ÞΔt, and Tð ¼ NΔtÞ is the sampled signal duration. A systematic way to end the iteration is to limit Dv to be a small number and/or to limit the maximum number of iterations. After c1 is obtained, define the residual r 1 ð≡u−c1 Þ, treat r 1 as a new signal, and repeat the steps shown in Eq. (9) to obtain other ci ði ¼ 2; …; nÞ as

obtained as   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi D1 a1 ¼ C 21 þ D21 ; θ1 ≡ωt s þ ϕ1 ¼ tan−1 C1 ω1 ¼ θ_ 1 ¼

ð14aÞ

    d D1 1 C1 ðt s −ΔtÞD1 ðt s Þ−C1 ðt s ÞD1 ðt s −ΔtÞ tan−1 ≈ tan−1 dt Δt C1 ðt s −ΔtÞC 1 ðt s Þ þ D1 ðt s ÞD1 ðt s −ΔtÞ C1

ð11Þ

ð14bÞ

The whole sifting process can be stopped when the residual r n ð ¼ r n−1 −cn Þ becomes a monotonic function from which no more IMF can be extracted. For data with a trend, the last IMF r n should be the trend and has at most one local extremum. Unfortunately, because EMD depends on the appearance of peaks and valleys to extract IMFs, it cannot well separate a signal's IMFs that only exist over some time segments of the sampled signal duration T. This intermittence problem and a signal conditioning technique for solving this problem will be illustrated by examples in Section 4.1. Hilbert–Huang transform (HHT) is a dual-domain method [13]. It uses the Hilbert transform shown in Eq. (3) to compute the conjugate part of each IMF and then uses Eqs. (4) and (6) to compute the time-varying amplitude and frequency. Unfortunately, the accuracy of time–frequency analysis of HHT is significantly deteriorated by the edge effect caused by HT. Teager–Kaiser algorithm (TKA) is another time-domain method for time–frequency analysis [6]. It is based on treating the signal as a regular harmonic u≡acosðωt þ ϕÞ and using finite difference to _ u€ and u; € and it requires 4 local data points approximate u; (i.e., ui−2 ; ui−1 ; ui and uiþ1 ) for frequency estimation [6]. Unfortunately, this four-point TKA cannot work for signals with moving averages, and we developed a five-point TKA to avoid this problem [18]. Even with the five-point TKA, because of the use of finite difference, results from TKA are inaccurate when the processed signal contains noise and/or high-frequency intra-wave modulation [18]. Next we propose a three-point method for extracting the timevarying frequency and amplitude of each IMF in the time domain. To initiate the extraction process for a specific IMF hðtÞ starting from t ¼ t s , we use TKA to estimate its frequency ω as [18] 0 1 1 ðhs −hs−1 Þ2 −ðhs−1 −hs−2 Þðhsþ1 −hs ÞA −1 @   cos ð12Þ 1− ω¼ 2 Δt 2 h −h h

Because there are only 3 unknowns in Eq. (13a), only 3 local data points (i.e., m ¼ 3) in Eq. (13c) need to be processed. However, more points can be used to reduce the influence of measurement noise and/or signal processing errors, but the estimated frequency is more a locally averaged value. For a noisy/complicated IMF, Eq. (12) may not provide an accurate frequency estimation. But, because the initial guess of ω from Eq. (12) is updated using Eq. (14b), this method is able to adapt and capture the actual variation of frequency. To further reduce the influence of noise on the calculated frequency, one can compute the ω at t ¼ t s by averaging over ð2p þ 1ÞΔt as

cn ¼ r n−1 −mn1 ⋯−mnK ; r n−1 ≡uðtÞ−c1 ⋯−cn−1

s

s−1 sþ1

where hs ≡hðt s Þ and t s ≡ðs−1ÞΔt: To enable frequency extraction of a signal having a moving average C 0 , we assume the signal to have the following form ^ ¼ C þ e cosðωtÞ−f sinðωtÞ ¼ C þ a cosðωt þ ϕ Þ hðtÞ 0 1 0 1 1 1 ¼ C 0 þ C 1 cosðωtÞ−D1 sinðωtÞ

ð13aÞ

where e1 ,f 1 , andC 0 are unknown constants, t ð≡t−t s Þ is a shifted local time, t s is the current time instant, and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f 2 a1 ≡ e21 þ f 1 ; ϕ1 ≡tan−1 1 ; C 1 ≡a1 cosðωt s þ ϕ1 Þ; e1 D1 ≡a1 sinðωt s þ ϕ1 Þ

ð13bÞ

The three unknowns C 1 , D1 and C 0 for the current data point at t ¼ 0ði:e:; t ¼ t s Þ can be determined by minimizing the squared error ðm−1Þ=2 Serror ≡∑i ¼ −ðm−1Þ=2 ðh^ sþi −hsþi Þ2

ð13cÞ

where mð≥3Þ is the total number of processed local data points around t ¼ t s and is assumed to be an odd number, h^ sþi denotes the function form shown in Eq. (13a) with t≡t i ¼ iΔt, and hsþi denotes the actual hðtÞ at t ¼ t s þ iΔt. After C 1 and D1 are determined, the instantaneous amplitude a1 , phase angle θ, and frequency ω can be

ω1 ≈

1 ∑p ω1 ðt s þ iΔtÞ 2p þ 1 i ¼ −p

ð15Þ

For any point other than the starting point at t ¼ t s , the previous point's frequency can be used as an initial guess of its frequency without using Eq. (12). Because the two ends of an IMF may be slightly distorted during the sifting process using natural cubic splines, it is better to start the frequency extraction process from a point away from boundaries. One can choose t ¼ t s to be around the midpoint and then process the data toward the two ends respectively. Note that the C 0 in Eq. (13a) is close to zero for an IMF extracted by EMD. This conjugate-pair decomposition (CPD) method shown in Eqs. (13a)–(15) can also be used for further decomposition of an IMF distorted by nonlinearities by assuming n

^ ¼ C þ ∑ ½e cosðω tÞ−f sinðω tÞ hðtÞ 0 i i i i i¼1 n

¼ C 0 þ ∑ ½C i cosðωi tÞ−Di sinðωi tÞ i¼1

where ei , f i and C 0 are unknown constants, and     qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Di d Di tan−1 ai ≡ C 2i þ D2i ; θi ≡tan−1 ; ωi ¼ dt Ci Ci

ð16aÞ

ð16bÞ

As shown later in Section 3, one can assume n ¼ 3 and ωi ¼ iω for an IMF of a nonlinear system with quadratic and cubic nonlinearities. 2.4. Time-domain-only noise filtering Because of the use of finite difference in Eq. (12), the accuracy of frequency estimation using TKA is easily destroyed by any noise. Noise filtering can be done in the time domain based on the use of transfer functions, but it often distorts signals [19,20]. Fortunately, the special characteristics of EMD enables the development of a time-domain-only noise filtering method. Because the displacement, velocity and acceleration are continuous at each node of a cubic spline and the m1k in Eq. (9) is the _ 1k average of the upper and lower cubic spline envelopes, m1k ,m € 1k are continuous time functions [21]. In other words, the and m _ and u€ remain in the first IMF c1 ðtÞ, and all discontinuities of u, u, other ci , c_ i , and c€ i are always continuous because ci ði≥2Þ are linear combinations of cubic splines. Hence, only c1 contains noise from the original signal. To take advantage of this special characteristics one can add a known high-frequency wave to the signal under processing and then extract the added wave and the signal's noise and discontinuities as the first IMF. After that, the discontinuities

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

and noise of the original signal can be accurately extracted by subtracting the added wave from the first IMF. Note that all the noise filtering, extraction of IMFs, and calculation of time-varying frequency and amplitude of each extracted IMF are done in the time domain, and the method is data-driven without the use of pre-determined basis functions or mathematical models. Next we show how to use the proposed CPD for parametric identification of nonlinear systems.

89

Imag a3

a a1

Real Fig. 1. Amplitude and frequency modulation of a Duffing oscillator.

3. Time-domain parametric identification of nonlinear systems We consider the following Duffing oscillator subjected to a harmonic excitation at a frequency Ω close to its linear natural frequency ω and the corresponding steady-state second-order perturbation solution [1]: u€ þ 2ςωu_ þ ω2 u þ αu3 ¼ FcosðΩtÞ

ð17aÞ

uðtÞ ¼ a1 cosðΩt−ϕÞ þ a3 cosð3Ωt−3ϕÞ ¼ aðtÞcosðΩt−ϕ þ ΘðtÞÞ; a3 ≡

αa31 32Ω2

≪a1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi aðtÞ≡ a21 þ a23 þ 2a1 a3 cosð2Ωt−2ϕÞ≈a1 þ a3 cosð2Ωt−2ϕÞ ΘðtÞ≡tan−1

a3 sinð2Ωt−2ϕÞ a3 ≈ sinð2Ωt−2ϕÞ a1 þ a3 cosð2Ωt−2ϕÞ a1

_ ¼Ωþ ^ Ω≡Ω þΘ þ 

3α 8ω

2 a61 −

2Ωa3 cosð2Ωt−2ϕÞ a1

 2 3αs 4 F a1 þ ðs2 þ ς2 ω2 Þa21 − ¼ 0; s≡Ω−ω 4ω 2ω

3α 2 a 8ω 1

_ ¼ a^ 1 sinðΩt−ϕÞ þ a^ 3 sinð3Ωt−3ϕÞ uðtÞ

ð17iÞ

where a^ 1 ¼ −Ωa1 , and a^ 3 ¼ −3Ωa3 . It represents a projection onto the vertical axis in Fig. 1. Hence, all the modulation phenomena presented in Eqs. (17a)–(17h)) for displacements are also valid for velocities. Velocities are considered here because a velocity sensor system is used in later experimental verifications.

ð17cÞ 4. Numerical and experimental results ð17dÞ

2Ωa23 þ 2Ωa1 a3 cosð2Ωt−2ϕÞ ≈Ω a21 þ a23 þ 2a1 a3 cosð2Ωt−2ϕÞ

  3αa21 s tan ϕ − ς¼ 2 ω 8ω ^ ¼ωþ ω

ð17bÞ

two harmonics in Eq. (17b) by using the proposed CPD, one can assign ω1 ¼ Ω, ω3 ¼ 3Ω, and all other harmonics to be zero in Eq. (16a). Numerical results show that all these characteristics under a harmonic excitation also exist in transient vibrations [22]. If the measured signals are velocities, it follows from Eq. (17b) that

ð17eÞ

ð17f Þ

ð17gÞ ð17hÞ

where ϕ is the phase difference between the forcing function and the response, and ς; α; F, and ϕ are constants. The amplitude a1 is a nonlinear function of F, ς, ω, and α, as shown in Eq. (17f) from first-order perturbation analysis [1]. Eq. (17b) shows that uðtÞ consists of two synchronous harmonics (i.e., reaching maxima at the same time), but, because a1 ≫a3 , uðtÞ would appear as one ^ varying distorted harmonic with an amplitude a and a frequency Ω at a frequency 2Ω. This amplitude- and frequency-modulation phenomenon can be easily explained by the vector plot shown in Fig. 1, and it can be used to determine the order (cubic or other) of ^ and a are at their nonlinearity. Moreover, if α 4 0, a3 =a1 4 0 , Ω ^ and a maxima when uðtÞ is at its maxima or minima (i.e., u_ ¼ 0), Ω _ are at their minima when uðtÞ is close to its maxima or minima ^ and a are at their minima when (i.e., u ¼ 0). If α o 0, a3 =a1 o 0, Ω ^ and a are at their uðtÞ is at its maxima or minima (i.e., u_ ¼ 0), and Ω _ maxima when uðtÞ is close to its maxima or minima (i.e., u ¼ 0). This phenomenon can be used to determine the type (hardening or softening) of nonlinearity, and the magnitude of nonlinearity can be estimated using Eq. (17b) and the time-domain data a1 ðtÞ and a3 ðtÞ asα ¼ 32Ω2 a3 =a31 . Moreover, one can use Eq. (17g) to estimate the damping ratio ς, use Eq. (17h) to estimate the ^ and use undamped amplitude-dependent natural frequency ω, Eq. (17f) to obtain F and then estimate the mass as m ¼ F 0 =F, where F 0 is the known excitation amplitude [1]. To separate the

Because vibration phenomena from numerical simulations using low-frequency signals are also valid for high-frequency signals if an appropriate scaling of the time is used [1], we will perform processing of low-frequency signals to reveal different physical and numerical phenomena. 4.1. Signal conditioning and noise filtering To show the time-domain signal decomposition using EMD we consider Eq. (18) that represents a signal consisting of three exponentially decaying harmonics, where Δt ¼ 1=300s is used and noise is a vector of normally distributed zero-mean random numbers with a standard deviation of 0:002 uðtÞ ¼ e−ζ 1 ω1 t sinðω1 tÞ þ 0:2e−ζ2 ω2 t sinðω2 tÞ þ 0:05e−ζ3 ω3 t sinðω3 tÞ þnoiseω1 ¼ 2π; ω2 ¼ 10π; ω3 ¼ 40π; ζ 1 ¼ ζ 2 ¼ ζ 3 ¼ 0:006

ð18Þ

Without the noise, Fig. 2 shows the signal and its first 3 IMFs extracted using the original EMD. Because EMD uses local peaks and valleys to sift an IMF, when an IMF is too small to create local peaks on a low-frequency riding wave, EMD shifts to extract the riding wave. This is why the beginning part of c1 is the 20 Hz component of uðtÞ; then the 5 Hz component, and then the 1 Hz component. Similarly, the beginning part of c2 is the 5 Hz component and then the 1 Hz component is extracted after the 5 Hz component is too small to create peaks on the 1 Hz riding wave. This intermittent behavior becomes worse when noise is included. Hence, signal conditioning techniques are needed for improving EMD. Fig. 3 shows the signal with the added noise and its first 3 IMFs extracted by adding 0:1sinðω3 tÞ to uðtÞ before sifting c~ 1 and then taking c~ 1 −0:1sinðω3 tÞ as c1 , then adding 0:1sinðω2 tÞ to the residue before sifting c~ 2 and then taking c~ 2 −0:1sinðω2 tÞ as c2 , and then adding 0:1sinðω1 tÞ to the residue before sifting c~ 3 and then taking c~ 3 −0:1sinðω1 tÞ as c3 . The ωi can be estimated from the signal's spectrum, and numerical simulations show that they do not need to be very accurate for the added harmonics before EMD analysis. For c2 and c1 , the minor distortion of amplitude at the beginning is caused by the EMD process. Note that the signal's noise is extracted into c1 as explained in Section 2.4.

90

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

Fig. 2. EMD decomposition of a transient response containing three exponentially decaying harmonics.

Fig. 3. EMD analysis with signal conditioning of a transient response containing three exponentially decaying harmonics.

4.2. Extraction of time-varying frequency and amplitude To demonstrate the proposed CPD for extracting time-varying amplitudes and frequencies we consider the following signal: uðtÞ ¼ ð2 þ 0:2sinωa tÞsinðωt þ sinωp tÞ þ noise; ω ¼ 2π; ωa ¼ 0:2ω; ωp ¼ 0:1ω

ð19Þ

where Δt ¼ 0:1 s is used and noise is a vector of normally distributed zero-mean random numbers with a standard deviation

of 0:005. Fig. 4 compares the original time-varying amplitude and frequency (thin blue lines) with those extracted by HHT and CPD with p ¼ 2 used in Eq. (7) for HHT and Eq. (15) for CPD. It reveals that both HHT and CPD are able to extract the time-varying frequency and amplitude, and HHT and CPD have different types of errors. The accuracy of HHT can be improved by using a smaller Δt because its accuracy depends on accurate capture of the signal's local peaks and valleys. On the other hand, the accuracy of CPD can be improved by using a larger Δt to increase the range covered by

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

91

Fig. 4. HHT and CPD analyses of Eq. (19) (a,b) from HHT analysis, and (c,d) from CPD analysis. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. HHT analysis of the nonlinear steady-state velocity of Eqs. (17a) and (20). (a) velocity, (b) spectrum, (c) amplitude, and (d) frequency.

92

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

the 3 points (i.e., 2Δt) used in Eq. (13c). This indicates that, for a fixed Δt, HHT is good for processing low-frequency signals, and CPD is good for processing high-frequency signals. 4.3. Nonlinearity identification To show the use of CPD for nonlinearity identification, we consider the nonlinear oscillator shown in Eq. (17a) with ω ¼ 2π; ς ¼ 0:01; α ¼ 0:2ω2 ; m ¼ 2; F 0 ð ¼ FmÞ ¼ 20; Ω ¼ ω

ð20Þ

Direct numerical integration is performed using the Runge– Kutta method with Δt ¼ 1=30 s for 150 s to damp out any transient effect. Fig. 5b shows that the steady-state velocity consists of several harmonics equally separated by 2Ω(¼2 Hz), indicating a frequency modulated at 2Ω. Fig. 5c,d shows that the time-varying frequency and amplitude from HHT analysis modulate at 2Ω, indicating a cubic nonlinearity (see Eqs. (17c) and (17e)). Moreover, because the frequency and amplitude are at their maxima _ is at its maxima, when u_ ¼ 0 and are around their minima when juj α 4 0. Velocities are considered here because a laser vibrometer is used to measure velocities in later experiments. If only one pair of conjugate functions are used in Eq. (16a), Fig. 6a, b shows that the time-varying amplitude and frequency from CPD do modulate at 2Ω and the amplitude is at its maxima when u_ ¼ 0, indicating a hardening cubic nonlinearity. However, the time-varying frequency does not clearly indicate a hardening nonlinearity and the variations of amplitude and frequency are not as accurate as those in Fig. 5c, d from HHT analysis. Fig. 6c, d shows that the amplitudes and frequencies from CPD using two pairs of conjugate functions in Eq. (16a) are constant, ω1 ¼ Ω, ω3 ¼ 3Ω, a1 ¼ 7:4502ð ¼ a1 ΩÞ, and a3 ¼ 0:2095ð ¼ 3Ωa3 Þ. Note that a1 and a3 are equal to the average

and varying amplitudes in Fig. 5c from HHT. Because there are five unknowns in Eq. (16a), at least 5 data points in Eq. (13c) are needed for frequency extraction. The transient effect at the beginning of Fig. 6c, d is caused by the use of TKA to estimate the frequency at t s ¼ 2Δt. One can use Eq. (17b) to estimate α ¼ 32Ω2 a3 =a31 ¼ 0:2133ω2 . Furthermore, we obtain ϕ ¼ 0:0961 from Fig. 5a and Eqs. (17a) and (17b) and ς ¼ ½3αa21 =ð8ωÞ−stan ϕ=ω ¼ 0:0108 from Eq. (17g) with sð ¼ Ω−ωÞ ¼ 0. Then Eq. (17f) with s ¼ 0 is used to obtain F ¼ 10:577ð ¼ F 0 =mÞ and hence m ¼ F 0 =F ¼ 1:891 is estimated. The key point here is that all α, ς and m can be estimated with an error less than 7% from just one steady-state response to a harmonic excitation at Ω≈ω. It follows from Eq. (17h) that the undamped natural frequency ^ ¼ 1:112ω, which indiunder a vibration amplitude a1 ¼ 1:1857 is ω cates that this is not really a weakly nonlinear problem. Because the perturbation solutions shown in Eqs. (17a)–(17h) are derived based on the assumption of weak nonlinearity, the estimation accuracy increases when the nonlinearity coefficient α decreases. Moreover, the estimations can be refined if more sets of responses to different F 0 and/or Ω are available. An alternative way to obtain Fig. 6c, d is to use EMD with the proposed signal conditioning technique by adding 3Ω and Ω harmonics, respectively, and then use CPD to extract the amplitude and frequency of each of the two IMFs. However, the EMD introduces the edge effect, and the results are not as good as those in Fig. 6c, d. Hence, it is better to use two pairs of conjugate functions for CPD without using EMD. Next we consider a transient response of the nonlinear oscillator shown in Eq. (17a) with _ ω ¼ 2π; ς ¼ 0:01; α ¼ 0:2ω2 ; m ¼ 2; uð0Þ ¼ 2; uð0Þ ¼0

a1 3Ω

Ω

a3

Fig. 6. CPD analysis of the nonlinear steady-state velocity of Eqs. (17a) and (20). (a,b) Using one conjugate pair, and (c,d) using two conjugate pairs.

ð21Þ

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

Fig. 7b shows that the response has a wideband spectrum ^ (see Eq. (17h)) changes with time. because its natural frequency ω Fig. 7c, d shows that the time-varying frequency and amplitude ^ and the frequency and from HHT analysis modulate at 2ω, amplitude are at their maxima when u_ ¼ 0 and are around their _ is at its maxima, indicating a hardening cubic minima when juj nonlinearity. However, the discontinuity-induced edge effect makes it difficult to estimate a3 from Fig. 7c. Fig. 7e, f shows the results from CPD using two pairs of conjugate functions. The hardening nonlinearity is obvious because the two frequencies decrease with the amplitudes. Moreover, the accurately extracted a1 and a3 can be used to estimate α and ς as shown in the harmonically exited case.

93

consider a horizontally cantilevered 727:0  31:75  3:175 mm3 steel beam subject to an initial tip displacement, as shown in Figs. 8 and 9a shows the transverse vibration velocity of a point at 267 mm from the clamped end measured using a Polytec PSV-200 scanning laser vibrometer with a sampling frequency of 1280 Hz (i.e., Δt ¼ 1=1280 s). The first six natural frequencies ωi are identified from the power spectral density (PSD) to be 4.45 Hz, 28.2 Hz, 78.7 Hz, 154.5 Hz, 255.4 Hz, and 382 Hz, and the mass density ρ is measured to be ρ ¼ 6951 kg=m3 . Hence, Young's

727.0mm×31.75mm×3.175mm steel beam 267mm

4.4. Experimental validation PSV-200 laser vibrometer

To demonstrate the use of CPD for time-domain-only signal processing for identification of nonlinear continuous systems, we

Fig. 8. Transient vibration of a cantilever subject to an initial tip displacement.

3ωˆ

a1 ωˆ a3

Fig. 7. Analysis of the nonlinear transient velocity of Eqs. (17a) and (21). (a) Velocity, (b) spectrum, (c,d) from HHT, and (e,f) from CPD using two pairs of conjugate functions.

94

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

Fig. 9. EMD analysis (a) the experimental velocity, and (b) IMFs from EMD analysis.

a1 ω3 ω2

a3

a2

ω1

Fig. 10. Time-varying amplitudes and frequencies of IMFs from HHT analysis (a) amplitudes, and (b) frequencies.

modulus E is estimated by matching the first five experimental natural frequencies with the theoretical linear frequencies (4.49 Hz, 28.16 Hz, 78.84 Hz, 154.49 Hz, and 255.39 Hz) to be E ¼ 149 GPa. Fig. 9b shows the three IMFs sequentially extracted using EMD with the proposed signal conditioning technique by adding and subtracting 0:5sinðωi tÞ, where i ¼ 3; 2; 1. We use 0:5sinðω3 tÞ for signal conditioning in order to extract the third and other highfrequency modal vibrations and noise into c1 and then obtain the second and first modal vibrations as c2 and c3 , respectively. Fig. 10 shows the time-varying amplitudes and frequencies of IMFs from HHT analysis. Because the first IMF c1 (mainly the third modal vibration) contains other high-frequency modes and measurement noise of the original signal, its calculated frequency ω3 shows significant variation in Fig. 10b. Moreover, Fig. 10b shows that, when an IMF's amplitude is small, its calculated frequency is noisy

because a small u in Eq. (6) makes tan−1 ðv=uÞ difficult to be accurate, especially when noise exists. For HHT, the inaccurate amplitudes and frequencies at the two data ends in Fig. 10 are caused by the discontinuity-induced Gibbs' effect discussed in Section 2.1. Fig. 11a shows that the spectrum of the first-mode vibration ^ 1 and (i.e., c3 in Fig. 9b) contains two major regular harmonics at ω ^ 1 , indicating a cubic nonlinearity (see Eq. (17b)). Fig. 11c, d shows 3ω the time-varying amplitude and frequency of the first-mode vibration extracted by using CPD to process the c3 in Fig. 9b after downsampling using Δt~ ¼ 20Δt (i.e., Three points and hence 2Δt~ cover about 14% of the period of c3 ). The amplitude and frequency in Fig. 11c,d are averaged over 3Δt~ and those in Fig. 11e, f are averaged over 13Δt~ . Fig. 11f clearly shows that the frequency increases from 4.45 Hz–4.48 Hz when the amplitude decreases, and Fig. 11d shows ^ 1 modulates at a frequency of 2ω ^ 1 . These two phenomena that ω

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

95

10 ς 1

10 ς 1

Fig. 11. Spectrum, amplitude, frequency, and damping ratio of the third IMF (a) PSD, (b) from HHT, (c,d) from CPD with averaging over 3 time steps, and (e,f) from CPD with averaging over 13 time steps.

together reveal that the first-mode vibration contains a softening cubic nonlinearity, which agrees with previous numerical results from multiple shooting analysis [23]. Note that the edge effect in Fig. 11b from HHT makes it difficult to observe the increase of frequency with the decrease of amplitude, even if 13Δt~ is used for averaging. Fig. 11e shows that the modal damping ratio decreases from 0.007–0.0035 when the amplitude decreases. This nonlinear damping effect is caused by the nonlinear longitudinal inertia due to the shortening effect, as explained later in Figs. 12–14. Figs. 10a and 11c, e show that the use of CPD avoids the discontinuity-induced edge effect on the extracted amplitude, and the obtained smooth amplitude makes it possible to have accurate estimation of damping. If the amplitude of the original signal shown in Fig. 9a is taken as the first-mode amplitude without modal decomposition (see Fig. 9b) to compute damping, erroneous results are obtained because its amplitude starts at about 1.0 but the actual amplitude a1 of c3 starts at about 0.75. CPD analyses are also done on c1 and c2 . The timevarying frequencies of both c1 and c2 also indicate the existence of a softening cubic nonlinearity.

Fig. 12a, b shows the time-varying amplitude and frequency from TKA analysis of the first-mode vibration (i.e., c3 in Fig. 9b) after down-sampling using Δt~ ¼ 20Δt and being averaged over 3Δt~ (the same conditions used in Fig. 11c, d). Although both CPD and TKA are based on the concept of circle fitting on the phase plane [18], TKA uses finite difference approximation but CPD uses curve fitting in numerical implementation. Hence, Fig. 12a, b shows that the time-varying frequency and amplitude from TKA are not as accurate as those from CPD. Note that c3 contains no noise because noise is extracted into c1 , as explained in Section 4.1. For a signal with noise (e.g., c1 ), the accuracy of TKA for time– frequency analysis is easily destroyed by the noise [18]. On the other hand, because CPD is based on curve fitting, noise filtering is an implicit capability and its accuracy increases with the number of processed data points. To explain the identified softening cubic nonlinearity and nonlinear damping we also perform nonlinear transient analysis of the beam using a geometrically exact nonlinear finite-element code GESA [23]. Six fully nonlinear Euler–Bernoulli beam elements

96

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

Fig. 12. TKA analysis of the third IMF (a) amplitude, and (b) frequency.

u, shortening

z (m)

F=25N

w

x (m)

N17

load (N)

N7

displacement (m) Fig. 13. Nonlinear static deformation (a) the deformed geometry at 25 N, and (b) the load-deflection curves of Nodes 7 and 17.

are used for x¼ 0 to 267 mm (i.e., the measured point at Node 7) and 10 elements for x ¼267 mm–727 mm and the Newmark-beta method is used for direct numerical integration. The initial deformed geometry (see Fig. 13a under F¼25 N) is obtained by applying a concentrated force of 14 N at the tip, and an integration time step Δt ¼ 1=2560 s and modal damping ratios of 0.003, 0.002, and 0.001 for the first, second, and other modes are used. The static load-deflection curve shown in Fig. 13b indicates that the system has a hardening cubic stiffness. The spectrum in Fig. 14b and the time-varying modal amplitudes in Figs. 14c and 10a show that the finite element solution is close to but not exactly the same as the experimental one. Except the first mode, other modal amplitudes in Fig. 14c are larger than the experimental ones in Fig. 10a. Moreover, the modal frequencies of the 5th, 6th, and 7th modes shown in Fig. 14b are slightly lower than the experimental ones. These are probably because the hand loading shown in Fig. 8 has an axial load and is different from the load shown in Fig. 13a, the actual modal damping might not be constant, and/or the actual mass distribution might not be uniform. Because constant modal damping ratios are used in the finite element model, a constant modal damping ratio of 0.003 for the

first mode is exactly extracted from the linear transient finiteelement solution from GESA. Hence, the nonlinear modal damping shown in Figs. 15a and 11c is believed to be caused by nonlinear modal coupling among the first few modes and the nonlinear longitudinal inertia induced by the shortening effect shown in Fig. 13a. Although the time-varying frequency in Fig. 15b does not clearly show softening effect (increase of frequency with decrease of amplitude) like that in Fig. 11f, zoom-in views show that the frequency modulates at 2ω ^ 1 and the frequency is at a local maximum when jdw=dtj is at its maximum, indicating a softening cubic nonlinearity (see Eqs. (17e) and (17b)). This softening cubic nonlinearity agrees with the experiment but is different from the hardening cubic nonlinearity predicted by perturbation analysis [23]. Perturbation analysis of weakly nonlinear vibrations of a cantilevered beam shows that the effective coefficient of cubic nonlinearity, αef f ective , is given by [23] Z L 2 ϕ½ϕ′ðϕ′ϕ″Þ′′dx; αef f ective ¼ α1 − α2 ω2 ; α1 ≡ 3 0 Z L Z xZ x α2 ≡ ϕ½ϕ′ ϕ′2 dxdx′dx ð22Þ 0

L

0

where ϕðxÞ is the first linear mode shape of the transverse displacement w and ϕ′≡∂ϕ=∂x. The α2 accounts for the shortening effect u shown in Fig. 13(a) and is based on the approximation Rx u ¼ −0:5 0 ðw0 Þ2 dx for weakly nonlinear vibration. For the first mode, αef f ective 40 is predicted by perturbation analysis, but both the experimental results and finite element simulations indicate αef f ective o 0. This indicates that the α1 in Eq. (22) does not correctly account for stiffening due to large deformation (see Fig. 13a and b), and/or the α2 in Eq. (22) does not fully account for inertia-induced softening due to shortening.

5. Conclusions A conjugate-pair decomposition (CPD) method is derived for time-domain-only signal processing of dynamic responses for dynamics characterization and system identification. CPD uses the empirical mode decomposition with signal conditioning techniques to avoid the intermittence problem in extraction of intrinsic mode functions (IMFs) and then use one or more pairs of conjugate functions to extract time-varying frequency and amplitude of each IMF. Because the variations of an IMF's frequency and amplitude can reveal system characteristics and nonlinearities, they can be used for dynamics characterization and system identification of discrete and continuum systems. CPD has an implicit noise filtering capability, and this capability can be

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

97

a1 ω3

a3

ω2

a2

ω1

Fig. 14. Nonlinear finite element solution (a) response, (b) spectrum, (c) amplitudes of IMFs from HHT analysis, and (d) frequencies of IMFs from HHT analysis.

10ς 1

Fig. 15. Amplitude, damping, and frequency of the first-mode IMF with averaging over 3 time steps (a) amplitude and damping ratio, and (b) frequency.

enhanced by processing more neighboring data points for each estimation. Furthermore, CPD is better than HHT for analysis of transient events because CPD is not affected by the edge effect caused by discontinuity. However, CPD needs the Teager–Kaiser algorithm to estimate the frequency of the first point in order to start the process, and it is accurate only if 3 consecutive data points can cover an enough portion (¼ or more) of the period of the component to be extracted. Moreover, for a signal with frequency and/or amplitude varying at a frequency close to or higher than the signal's fundamental frequency, the accuracy of CPD decreases. Discrete systems' simulated responses and experimental beam vibrations are used to validate the capability and accuracy of CPD.

Acknowledgment This work is supported by the National Science Foundation under Grants CMMI-1039433 and IIP-1026883 and the NASA Headquarters under Grant NNX09AV08A. The support is gratefully acknowledged. References [1] A.H. Nayfeh, D.T. Mook, Nonlinear Oscillations, Wiley-Interscience, New York, 1979. [2] F.C. Moon, Chaotic Vibrations, An Introduction for Applied Scientists and Engineers, Wiley-Interscience, New York, 1987.

98

P. Frank Pai et al. / International Journal of Non-Linear Mechanics 54 (2013) 85–98

[3] G. Kerschen, K. Worden, A.F. Vakakis, J.C. Golinval, Past, present and future of nonlinear system identification in structural dynamics, Mechanical Systems and Signal Processing 20 (2006) 505–592. [4] S.W. Doebling, C.R. Farrar, M.B. Prime, D.W. Shevitz, Damage Identification and Health Monitoring of Structural and Mechanical Systems from Changes in their Vibration Characteristics: A Literature Review, Report No. LA-13070-MS, Los Alamos National Laboratory, 1996. [5] H. Sohn, C.R. Farrar, F.M. Hemez, D.D. Shunk, D.W. Stinemates, B.R. Nadler, A Review of Structural Health Monitoring Literature: 1996–2001, Report No. LA13976-MS, Los Alamos National Laboratory, 2003. [6] P. Maragos, J.F. Kaiser, T.F. Quatieri, Energy separation in signal modulations with application to speech analysis, IEEE Transactions on Signal Processing 41 (1993) 3024–3051. [7] C. Mei, K. Abdel-Motagaly, R. Chen, Review of nonlinear panel flutter at supersonic and hypersonic speeds, Applied Mechanics Reviews 52 (2000) 321–332. [8] C.H.M. Jenkins (Ed.), Gossamer Spacecraft: Membrane and Inflatable Structures Technology for Space Applications, AIAA, Reston, Virginia, 2001. [9] E.O. Brigham, The Fast Fourier Transform, Prentice-Hall, Englewood Cliffs, New Jersey, 1974. [10] D.F. Shi, L.S. Qu, N.N. Gindy, General interpolated fast Fourier transform: a new tool for diagnosing large rotating machinery, Journal of Vibration and Acoustics 127 (2005) 351–361. [11] P. Flandrin, Time–frequency/Time-Scale Analysis, Academic Press, San Diego, CA, 1999. [12] I. Daubechies, Ten Lectures on Wavelets, CBMS-NSF Lecture Notes 61, SIAM, Philadelphia, 1992.

[13] N.E. Huang, N.O. Attoh-Okine (Eds.), The Hilbert–Huang Transform in Engineering, CRC Press, Boca Raton, FL, 2005. [14] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proceedings of the Royal Society of London A 454 (1998) 903–995. [15] N.E. Huang, M.C. Wu, S.R. Long, S.S.P. Shen, W. Qu, P. Gloersen, K.L. Fan, A confidence limit for the empirical mode decomposition and Hilbert spectral analysis, Proceedings of the Royal Society of London A 459 (2003) 2317–2345. [16] S. Rossignol, P. Desain H. Honing, State-of-the-art in fundamental frequency tracking, in: Proceedings of Workshop on Current Research Directions in Computer Music, Barcelona, Spain, 2001, pp. 244–254. [17] S.L. Hahn, Hilbert Transforms in Signal Processing, Artech House, Boston, 1996. [18] P.F. Pai, Circular instantaneous frequency, Advances in Adaptive Data Analysis 2 (2010) 39–64. [19] A.V. Oppenheim, R.W. Schafer, Discrete-Time Signal Processing, Prentice Hall, Englewood Cliffs, NJ, 1989. [20] R.M. Mersereau, M.J.T. Smith, Digital Filtering, Wiley-Interscience, New York, 1981. [21] R.L. Burden, J.D. Faires, Numerical Analysis, 3rd ed., Prindle, Weber & Schmidt, Boston, Massachusetts, 1985. [22] P.F. Pai, A.N. Palazotto, Detection and identification of nonlinearities by amplitude and frequency modulation analysis, Journal of Mechanical Systems and Signal Processing 22 (2008) 1107–1132. [23] P.F. Pai, Highly Flexible Structures: Modeling, Computation and Experimentation, AIAA, Reston, Virginia, 2007.