- Email: [email protected]
Hilbert transform methods for nonparametric identification of nonlinear time varying vibration systems
Mechanical Systems and Signal Processing 47 (2014) 66–77
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp
Hilbert transform methods for nonparametric identification of nonlinear time varying vibration systems Michael Feldman n Faculty of Mechanical Engineering, Technion—Israel Institute of Technology, Technion City, 32000 Haifa, Israel
a r t i c l e in f o
abstract
Article history: Received 2 January 2012 Received in revised form 31 August 2012 Accepted 2 September 2012 Available online 26 September 2012
The objective of this paper is to demonstrate a Hilbert transform (HT) method for identification of mechanical time-varying vibration systems under free and forced vibration regimes. This special kind of non-stationarity arises in experimental vibration analysis and in engineering practice. The method is based on the HT of input/output signals in a time domain to extract instantaneous dynamic structure characteristics, such as natural frequencies, stiffness, damping, and their variations in time. The HT assigns a complementary imaginary part to a given real signal part, or vice versa, by shifting each component of the signal by a quarter of a period. Thus, the HT pair provides a method for determining the instantaneous amplitude and the instantaneous frequency of a signal. For general non-stationary vibration signals, the analytic signal method does a good job of simultaneous time-frequency localization of the main signal components. The paper focuses on HT signal processing techniques and identifies three groups of dynamics time-varying SDOF systems: slow varying quasi-periodic modulation of stiffness under free and forced vibration, slow varying quasi-periodic modulation of nonlinear stiffness under free vibration, and fast inter-wave parametric stiffness modulation (Mathieu equation). & 2012 Elsevier Ltd. All rights reserved.
Keywords: Non-linear time-varying system Hilbert transform Nonparametric identification Envelope Instantaneous frequency
1. Introduction Typically, a vibration system concerned with a mechanical construction having mass, stiffness and damping elements is a physical structure that, when taking an impulse or another input force, usually initiates an output natural vibration motion. If the system has a time-varying mass or stiffness parameter, the natural frequency will depend decisively on the variation. If the system also has a time-varying damping, the solution amplitude will react in situ with the damping variation. Such frequency and amplitude variations can be clearly detected by any suitable signal processing technique in the time or frequency domain. The appropriate signal processing is mostly based on the extracted features of the level and frequency content of a varying signal in the form of the envelope and the instantaneous frequency. Thus, for example, a moving least-squares technique and a Kalman filter allow time-varying frequency tracking, while the traditional Fourier transform and wavelets technique can identify online variation in the natural frequencies and mode shapes of a vibration system [1–3].
n
Tel.: þ 972 4829 2106; fax: þ 972 4829 5711. E-mail address: [email protected]
0888-3270/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ymssp.2012.09.003
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
67
If identification of the varying dynamic system requires additional a priori information about the dynamic system and its vibration model, this identification is called parametric, and its aim is to determine the parameters of a known model. A dynamic system model is commonly created on the basis of knowledge of a mechanical construction and a functional operation of a technical object [4]. In some other more sophisticated cases, a nonparametric analysis may be required to estimate the parameter variation of the unknown physical equation of the dynamic system [5,6]. Nonparametric system identification of the modal parameters, such as frequencies, damping, and mode shapes as functions of the physical inertia, stiffness, and damping properties, is based only on analysis of the system behavior. Recent works in the area of time domain representations of vibration based on the HT show great promise for applications in nonparametric vibration system identification. An outstanding work by N. Huang, known as Empirical Mode Decomposition, adaptively decomposes a signal into the simplest intrinsic oscillatory modes (components) in the first stage. Then, in the second stage, each decomposed component forms a corresponding set of instantaneous amplitudes and frequencies called the Hilbert spectrum that can be used to obtain time-frequency vibration representation (characterization). Every significant decomposed elementary component holds exact amplitude and frequency information contained in the original signal, and the HT method makes it sensitive to any modal shift. Practically all publications dedicated to HT time domain nonlinear vibration identification with envelope and instantaneous frequency use the approach suggested in [7]. This FREEVIB approach introduces a vibration solution as a new complex variable which, when substituted into an initial differential equation, directly determines the system’s natural frequency and damping parameters. The proposed methods for identifying instantaneous modal parameters (natural frequencies, damping characteristics and their dependencies on a vibration amplitude and frequency) prove to be very simple and effective. The HT-based signal processing in parallel with other time-frequency representations in the form of Wavelet transform, Gabor Transform, or Wigner–Ville distribution provides a powerful tool for vibration analysis and identification of vibration systems [8]. Of particular interest and importance is the use of the HT to interpret nonlinear system motion. The measured instantaneous characteristics, such as the instantaneous frequency or the amplitude of the nonlinear solution, take an unusual fast oscillation (intrawave modulation) form. In the case of supplementary parameter variation, both instantaneous functions of each vibration mode – the envelope and the instantaneous frequency – become complicated modulated functions: they reflect the slow time variation of the stiffness plus fast intrawave modulations due to nonlinearity [8]. Therefore, only a smeared or blurred form of the instantaneous frequency will be obtained through direct frequency monitoring. The main purpose of this paper is to demonstrate the method of nonparametric identification of a vibration system when the linear time-varying system description fails because the non-linearity is too severe and a very accurate model is needed. In this regard, not only must the presence of a parameter variation and nonlinearity be detected, but also an adequate and readily interpretable system model must be identified. It has been suggested that HT processing involving Hilbert vibration decomposition and congruent functions can identify the time-frequency distribution of a signal with great accuracy. The digital HT signal processing the same as every other approach has some inherent limitations. The method is based on the multipoint FIR digital Hilbert transformer and therefore it requires significantly longer initial data to neglect the filter end effects. The method has a good but an ultimate frequency resolution required for separation of closely spaced harmonics. It allows separation the finite small number of valued components, which does not exceed 5–6 components. It is very sensitive to measured instrumental noise. The HT decomposition method is not effective for the separation of nonoscillation types of motion, such as random, impulse, slow aperiodic signals. 2. Hilbert transform-based identification and analysis of time varying vibration systems 2.1. The Hilbert transform in time domain approach. The HT assigns a complementary imaginary part x~ to a given real signal part x, or vice versa, by shifting each component of the signal by a quarter of period x~ ¼ H½x. Thus, the HT pair provides a method for determining the instantaneous qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ amplitude AðtÞ ¼ 7 x2 ðtÞ þ x~ 2 ðtÞ and the instantaneous frequency of a signal oðtÞ ¼ d½arctan ðxðtÞ=xðtÞÞ=dt. For a general non-stationary vibration signal, the analytic signal method does a better job of simultaneous timefrequency localizing of the main signal components [8]. The analytic signal method is equally applicable to deterministic and random processes, although, generally speaking, it does not divide them into two separate groups. Hence, it enables us to investigate any oscillating time function from a more general point of view. The method is also good for solving problems concerning analysis of stationary and non-stationary vibrations, as well as narrow- and/or wideband multicomponent signals. It also allows precise analysis of the transformation and dissipation of vibration energy, and vibration effects on machine durability. 2.2. Advantages of Hilbert transform time domain identification Taking analytic signal representation into account enables one to consider a vibration process, at any moment in time, as a quasi-harmonic oscillation, amplitude- and frequency-modulated by time-varying functions A(t) and
68
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
oðtÞ : xðtÞ ¼ AðtÞcos
Rt
0 oðtÞdt. The instantaneous parameters are functions of time and can be estimated at any point of the vibration signal. The total number of points that map the vibration is much greater than the number of peak points of the signal. Therefore the HT provides the best available time resolution of the amplitude and other estimations of vibration analysis. It also opens the way for averaging and for other statistical processing procedures, making vibration analysis more precise [8]. Classically, frequency resolution means the smallest difference or granularity with which one can distinguish details in the frequency domain of a signal. For any other digital method the spectral resolution is a discrete and terminal value that cannot be less than the ratio between the sampling frequency and the number of measured points. In contrast, the HT calculates the frequency as the rate of phase change with time, and the instantaneous frequency resolution can be as small as desired depending on the accuracy of the signal phase and time measurement. Practically all suggested HT computational algorithms are simple and fast, thus allowing online (real-time) signal analysis and identification. The HT methods are suitable for analysis and identification of both linear and nonlinear vibration systems; they produce adaptive decomposition of complicated multicomponent vibration for elementary inherent components [9]. With the help of the HT, researchers get new features, such as the envelope and the instantaneous frequency of each component and the instantaneous phase relations between time-varying components. With regard to system identification, HT methods allow estimating not only the varying stiffness but also the varying damping in vibration constructions.
3. Quasi-linear slow time varying system identification Two methods of dynamic system studying based on the HT have recently been proposed for free and forced vibration analysis, when the input signal is an impulse, shock or sweeping harmonic. The methods are suitable for testing linear and non-linear systems and for instantaneous modal parameter identification, including determining the concrete type of nonlinear spring and damping characteristics for each mode of a vibratory system. The time-varying structural system encountered in civil and mechanical engineering can be described by the following _ and K(t)x are respectively equation of motion when a SDOF system is considered: mðtÞx€ þ CðtÞx_ þ KðtÞx ¼ 0 where m(t), CðtÞx, _ and the restoring force as a function of the displacement x. the varying mass, the frictional force as a function of velocity x, The corresponding second-order differential equation of the damped system for unit mass will then take the following form: x€ þ 2hðtÞx_ þ o20 ðtÞx ¼ 0, where the term hðtÞx_ ¼ CðtÞ=2mðtÞ represents the damping force per unit mass as a function of the velocity, h(t) is the damping coefficient, and the term o20 ðtÞ ¼ KðtÞ=mðtÞ represents the time varying undamped natural frequency squared. The instantaneous natural frequency of the solution according to the FREEVIB method is estimated by the following formula [7]: 2
€ þ2A_ =A2 þ A_ o _ =Ao o20 ðtÞ ¼ o2 A=A
ð1Þ
Here A is an envelope, and o are instantaneous frequency functions of the solution. But for many real engineering structures that are under force excitation, it is useful to consider a method of modal analysis of non-linear systems with input signal excitation. The modern technique for non-linear system investigation based on the HT and called FORCEVIB [8] enables us to identify instantaneous system modal parameters (natural frequencies, damping characteristics and their dependencies on vibration amplitude and frequency) during different kinds of excitation of the dynamic system. It allows direct extraction of the linear and non-linear system parameters from the measured time signal of input and output. The HT method of forced vibration x€ þ 2hðtÞx_ þ o2x ðtÞx ¼ zðtÞ determines instantaneous modal parameters even if the input signal is a fast sweeping frequency quasi-harmonic signal: _ omA=A € þ 2A_ 2 =A2 þ A_ o _ =Ao, where o0(t) is an instantaneous modal frequency, h(t) is an o20 ðtÞ ¼ o2 þ aðtÞ=m þ bðtÞA=A instantaneous modal damping coefficient, o and A are instantaneous frequency and envelope of the vibration with their _ A, € o _ Þ, and aðtÞ ¼ ReðZ=XÞ, bðtÞ ¼ Im ðZ=XÞ are real and imaginary parts of the input–output first and second derivatives ðA, signal ratio. Below we consider some simulated examples of identification of time-varying differential equations of motion performed in SIMULINK (MATLAB) with the permanent step value ODE4 (Runge–Kutta) solver.
3.1. Free vibration regime As a first example of a quasi-linear time-varying system, we choose a SDOF vibration model with the following square wave modulated stiffness: x€ þ0:05x_ þ ½1þ 0:6sgnðcos 0:15tÞx ¼ 0; x0 ¼ 1. The free vibration solution for nonzero initial condition and the sampling frequency 13 Hz are shown in Fig. 1. The variation of the solution and its envelope reflect the time modulated stiffness. The identified instantaneous natural frequency squared function (the instantaneous modal stiffness) calculated according to (1) accurately repeats the initial square wave stiffness (Fig. 1b). The identified instantaneous modal stiffness allows presenting the static force characteristic k(x) (Fig. 2), which for the time-varying case twists around the initial linear force characteristic.
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
69
1 vibration envelope
x
0.5
0
−0.5
0
20
40
60
80
100
120
2.5
140
160
180
200
Instantaneous frequency squared stiffness
Stiffness
2 1.5 1 0.5 0 0
20
40
60
80
100
120
140
160
180
200
Time, s Fig. 1. Linear system: the free vibration (a), the stiffness variation (b).
1.5
time−invariant stiffness time−varying stiffness
1
Spring force
0.5
0
−0.5
−1
−1.5 −1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
Displacement Fig. 2. Linear system: the stiffness force characteristics.
3.2. Forced vibration regime As the next example, let us consider a structure that describes a slow harmonic modulated elastic force of the trivial dynamics system under external harmonics excitation with a fixed frequency: x€ þ0:05x_ þ ð1 þ 0:5cos 0:03tÞx ¼ cos 2p0:16t (Fig. 3).
70
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
Excitation
1 0
Displacement
−1 0
100
200
0
100
200
300
400
500
300
400
500
20 0 −20 Time, s Solution (Displacement) Power Spectrum Magnitude
Power Spectrum Magnitude
Excitation
0.2 0.15 0.1 0.05 0
6 5 4 3 2 1 0
0
0.1 0.2 Frequency, Hz
0.3
0
0.1 0.2 Frequency, Hz
0.3
Fig. 3. Forced excited system: the excitation (a), the displacement (b), the excitation spectrum (c), and the displacement spectrum (d).
Displacement and Envelope 20
10
0
−10
−20 50
100
150
200
250
300
350
400
450
500
550
400
450
500
550
Instantaneous Frequencies
Frequency, Hz
0.2 0.18 0.16 0.14 0.12 50
100
150
200
250
300
350
Time, s Fig. 4. Forced excited system: the displacement and the envelope (a) and the instantaneous frequencies (b).
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
71
The generated vibration has a rather complicated form over time and also in the frequency domain (Fig. 4). But HT identification restores this modulation in accurate detail [10]. Thus, Fig. 4b shows that the identified instantaneous natural frequency of the vibration (bold line) completely coincides pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with the varying initial elastic force modulation function f ¼ 1 þ 0:5cos 0:03t =2p (dashed line). These two examples demonstrate that the HT method enables identification of a slow modulation of the linear system parameters, in both cases the free and the forced vibration regimes. Therefore the HT-based methods proposed for identifying instantaneous modal parameters prove to be very simple and effective [10]. 4. Nonlinear time varying system identification 4.1. Nonlinear interwave system representation The nonlinear restoring force k(x) in the second-order differential equation of a vibration system can be recast into multiplication of a varying non-linear natural frequency o20 ðtÞ and system solution x : x€ þ 2h0 ðtÞx_ þ o20 ðtÞx ¼ 0. To apply the multiplication property of the HT, let us assume that the new varying nonlinear natural frequency squared could be grouped into two different parts: o2 ðtÞ ¼ o20 ðtÞ þ o21 ðtÞ. The first part o20 ðtÞ is much slower, and the second component o21 ðtÞ is faster than the system solution, so the equation of motion takes a new form: x€ þ ½o20 ðtÞ þ o21 ðtÞxðtÞ ¼ 0. Now, according to the multiplication property of the HT for overlapping functions [8], we apply the HT to both sides of the last equation ~ þo ~ 21 ðtÞxðtÞ ¼ 0, where a tilde sign indicates the HT conjugate variable. Multiplying each side of the obtained HT x~€ þ o20 ðtÞxðtÞ equation by i and adding it to the corresponding sides of the initial equation yields a new form of a differential equation of ~ 21 Þx ¼ 0, where X is a complex solution in an analytic signal form: motion in the analytic signal form X€ þ o20 X þ ðo21 þio ~ This complex equation can be transformed into a more traditional and accepted form [8] X ¼ x þ ix. X€ þ idX þ o2 X ¼ 0
ð2Þ
2 2 ~ 21 xxÞ=A ~ ~ 21 x2 o21 xxÞ=A ~ where o2 ¼ o20 þ ððo21 x2 þ o Þ; d ¼ ððo Þ. Here o2 is a varying instantaneous natural frequency and d is a fast varying instantaneous fictitious friction parameter. It should be pointed out that this equation is not a real equation of nonlinear motion; however, it produces the same varying solution of the system. It is simply an artificial fictitious equation that yields the same nonlinear primary vibration solution. The obtained varying natural frequency o2 consists of a slow o20 and a fast component. The fictitious friction parameter d, on the contrary, consists only of a fast varying component. This means that the slow part of the natural frequency forms an average period of vibration, but the fast varying fictitious friction force (induced by the zero average value) does not affect the real average friction force. It should be noted that the inherent non-linear restoring force causes fast oscillation of both the instantaneous frequency and the fictitious instantaneous damping coefficient, even in a conservative non-linear system.
4.2. Congruent modal parameters When the instantaneous modal parameters of a nonlinear system vary with time much more rapidly than the system oscillation solution, it is desirable to squeeze these fast variations into smooth initial nonlinear functions. Such precise nonlinear identification paves the way for tracing the slow time-varying parameters of nonlinear systems. Precise identification of the smooth initial nonlinear static stiffness force is based on the displacement and the natural frequency functions corresponding to the potential energy maxima. Such functions were defined as a congruent envelope and a congruent instantaneous natural frequency [8]. The congruent envelope is the envelope of the envelope function; it is equal to the algebraic sum of the solution superharmonic envelopes Ac ðtÞ ¼
N X
Al ðtÞcosfl ðtÞ
ð3Þ
l¼1
where Al(t) is the envelope of the solution l order harmonic, fl(t) is the phase angle between the primary and the l order harmonic. For example, if all high order superharmonics are congruently in phase with the primary solution, the congruent envelope is equal to the sum of harmonic envelopes. If superharmonics are in different phase relations with extreme upper and lower points of the largest energy component, this governs the symmetry – or asymmetry – of the multicomponent signal. It is clear that estimation of the congruent envelope Ac(t) requires a preliminary time domain decomposition of the system solution into a sum of high order harmonics. In turn, the varying intrawave instantaneous natural frequency o0(t) can also be decomposed into a sum of high order synchronous components. So the congruent natural frequency will be
oc ðtÞ ¼
N X
o0l ðtÞcosfol ðtÞ
ð4Þ
l¼1
estimated as an algebraic sum of the superharmonics of the instantaneous natural frequency. Both the congruent envelope and the natural frequency, being smooth functions, define the exact static stiffness force of a nonlinear system for every
72
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
moment in time: o2c Ac . In practice, the accuracy of the last expression depends on the total number N of considered synchronous vibration components. For symmetric force characteristics, the initial smooth static force function can be written in the form: ( 2 oc Ac , x 40 kðxÞ ¼ , o2c Ac , x o 0 where o2c is the congruent modal frequency squared, and Ac is the congruent envelope of the displacement. Congruent functions of the asymmetrical vibration systems are split into two different parts, separated according to the positive and the negative displacements: ( ( ocp ðtÞ, if x 4 0 Acp ðtÞ, if x 4 0 Ac ðtÞ ¼ ; oc ðtÞ ¼ : ocn ðtÞ, if x r0 Acn ðtÞ, if x r0 Therefore, the nonparametric identification of the asymmetric vibration systems considers together these two solution parts in the time domain. As a result of the asymmetric signal, the decomposition will have two sets of envelopes and instantaneous frequencies, each for its own positive and negative part of the displacement. Then, the separated congruent functions construct two different static force characteristics—separately for the positive and the negative signal parts ( 2 ocp Acp , x 40 kðxÞ ¼ ð5Þ o2cn Acn , x r 0 To estimate the congruent functions, we need to produce a decomposition of the signal first and then to construct an algebraic summation of envelopes of all the decomposed components. Such a successive signal decomposition (disassembling) and the subsequent summation (re-assembling) of the component envelopes generates desirable congruent functions as slow functions of time suitable for a precise nonparametric identification.
4.3. Nonlinear system slow stiffness variation. free vibration regime For the nonlinear time-varying system example, we can construct a vibration motion by forming a combination of the following parameters: x€ þ 0:05x_ þð1 þ0:6cos 0:15tÞð1 þ13x2 Þx ¼ 0, x0 ¼ 2 with the sampling frequency 13 Hz (Fig. 5). The corresponding solution is a free vibration of the nonlinear Duffing model whose stiffness is slowly modulated by the harmonic function. Without such slow harmonic modulation, the solution would be a classic Duffing nonlinear solution. 1.8 1.7
x
1.6 1.5 1.4 1.3 1.2 20
22
24
26
28
30 vibration envelope congruent envelope
1.5 1
x
0.5 0 −0.5 −1 −1.5 0
20
40
60
80
100
120
140
Time, s Fig. 5. Free vibration of the Duffing system.
160
180
200
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
73
Both instantaneous functions of the solution – the envelope (Fig. 5) and the instantaneous frequency (Fig. 6a) – are complicated modulated functions; they reflect the slow time variation of the stiffness plus fast intrawave modulations due to cubic nonlinearity [8]. Therefore, only a smeared or blurred form of the instantaneous frequency (Fig. 6b) is available for direct frequency monitoring. A recently developed precise nonlinear identification based on the congruent frequency function (4) squeezes the smearing and blurring spots in a single clean line that determines the real nonlinear stiffness as a function of time (Fig. 7a). By applying the congruent frequency, we eliminate the fast intrawave modulations, though the desired slow modulation will still be present in the varying stiffness function (Fig. 7a).
Fig. 6. Instantaneous frequency of Duffing system: the IF and the congruent frequency (a) and the spectrogram (b).
35 congruent stiffness smoth congruent stiffness
Stiffness, (rad/s)2
30 25 20 15 10 5 0 0
20
40
60
80
100
120
140
160
180
200
Stiffness variation
2.5 initial stiffness variation identified stiffness variation
2 1.5 1 0.5 0 0
20
40
60
80
100
120
140
160
180
200
Time, s Fig. 7. Duffing system frequency and stiffness: the congruent and the average frequencies (a) and the varying stiffness (b).
74
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
100 time−invariant stiffness time−varying stiffness
80
60
40
Spring force
20
0
−20
−40
−60
−80
−100 −2
−1.5
−1
−0.5
0
0.5
1
1.5
2
Displacement Fig. 8. Duffing varying system: the stiffness force characteristics.
Let assume that the monotonic form of the varying stiffness is an inherent feature of the nonlinear system that is not related to the desired slow varying modulation. Now we can remove the monotonic form by using a standard polynomial fitting technique in a least squares sense. The remain wave-like instantaneous stiffness function is plotted on Fig. 7b together with the initial harmonic modulation function 1 þ0.6 cos 0.15t. As can be seen, both the initial and the identified time-varying stiffness functions are in good agreement. The identified time-varying stiffness multiplied by the congruent envelope function constructs the nonlinear stiffness force characteristics, which in our case takes the form of a spiral line twisted around the initial nonlinear cubic force characteristic (Fig. 8). 4.4. Fast parametric excitation regime If a slow varying stiffness parameter is the case, the frequency of the solution will simply follow and reflect this slow variation. For an artificial fast parameter variation, in contrast, the behavior of the system solution can be quite different in principle. The fast stiffness modulation with the doubled (relative to the natural frequency) and faster multiple frequencies could even excite the system parametric resonances. Time domain identification will require a priori knowledge about the existence of the fast varying parameter and about the vibration system structure. This means that identification of the fast time-varying parameters is not a purely non-parametric problem. As an example, we describe the process of identifying a vibration system with parametric excitation in the form of linear differential equation with variable stiffness coefficient (Mathieu equation): x€ þ2h0 x_ þ o20 ð1 þ esinof ast tÞx ¼ 0, where e is the parametric modulation depth. We choose a stiffness coefficient that alternates harmonically twice over the period of the excited oscillations ofast ¼ 2o0. To identify the equation with parametric excitation, we use the non-linear system representation (2) and the R quasi-harmonic form of the solution x ¼ AðtÞcos odt. After the algebraic transform, the initial parametric equation takes a new form x€ þ ½2h0 0:5eo0 ð1cos 2o0 tÞx_ þ o20 ð1 þ0:5esin 2o0 tÞx ¼ 0
ð6Þ
This form now includes the fast varying stiffness factor 1 þ0:5esin2o0 t and also the fast varying damping component 0:5eo0 cos 2o0 t. The objective of the identification could be the unknown parametric modulation depth e. The obtained form of the equation illustrates the well-known fact that parametric excitation essentially affects system damping, reducing the real damping up to zero. The critical condition 2h0 o0.5eo0 causes loss of stability and appearance of unlimited amplitude vibration. For our Matlab simulation, we chose an initial numeric model describing parametric oscillation in the form: x€ þ 0:1x_ þ ð1 þ0:3sin 2tÞx ¼ 0, x0 ¼ 1, where o0 ¼1, e ¼0.3, 2h¼0.1. The obtained parametrically unstable excited vibration of
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
75
the system, together with the vibration envelope, is shown in Fig. 9a. The executed vibration has a frequency that is two times smaller in comparison with the parametric excitation frequency. After using the FREEVIB identification, we get the instantaneous natural frequency squared shown in Fig. 9b, together with the initial modulated stiffness function 1þ0.3 sin 2t. For the case of fast parameter modulation, the estimated instantaneous natural frequency level is equal only to half of the initial elastic parametric variation. This means that the direct time-frequency monitoring of the solution does not provide the correct value of the real stiffness variation.
100
vibration envelope
50
x
0
−50
−100 0
40
60
80
100
120
140
160
180
200
100
120
140
160
180
200
instantaneous frequency squared initial stiffness variation
1.4
Stiffness
20
1.2 1 0.8 0.6 0
20
40
60
80
Time, s Fig. 9. Parametrically excited system: the vibration (a) and the instantaneous frequency and the stiffness variation (b).
Stiffness
1.5
identified variation initial variation
1
0.5
Damping coefficiemt
0
20
40
60
80
100
120
140
160
180
200
identified variation theoretical variation initial constant value
0.2 0.1 0 −0.1 −0.2 0
20
40
60
80
100
120
140
160
180
200
Time, s Fig. 10. Parametrically excited system: the identified and the initial stiffness variation (a) and the identified and the theoretical damping variation (b).
76
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
time−invariant stiffness time−varying stiffness
100
Spring force
50
0
−50
−100
−100
−80
−60
−40
−20
0
20
40
60
80
100
Displacement Fig. 11. Parametrically excited system: stiffness force characteristics.
To identify the unknown parametric modulation depth e, we should consider the theoretical fast varying stiffness factor (6) 1 þ0:5esin2o0 t. The resultant fast varying elastic component for the estimated modulation depth e ¼0.1 is plotted together with the initial modulated stiffness function (Fig. 10a). As can be seen, both the estimated and the initial functions practically coincide. The discovered modulation depth e ¼0.1 can be also used to compare the theoretical damping variation (6) and the real time-varying damping function obtained through the FREEVIB method (Fig. 10b). Both these damping functions are in very good agreement, and they lie in the negative half-plane below the constant initial damping coefficient line (2h0 ¼ 0.1). The observed fast oscillation stiffness force can be represented as a function of the displacement in the form of traditional force characteristics (Fig. 11). Due to the fast parametric excitation, the shown force characteristics take the form of a spiral repeatedly wrapped around the initial linear force characteristic. 5. Conclusions The direct way to identify time-varying instantaneous modal parameters of nonlinear vibration systems is to use HT-based vibration decomposition and the appropriate HT time domain signal processing. The modern HT signal decomposition divides the real multicomponent solution into a principal nonstationary harmonic and a number of separate high-frequency superharmonics. A solution of nonlinear time-varying vibration systems can also include a slow variation component. By considering the high superharmonics, the congruent smooth envelope and the natural frequency, we are able to identify, precisely and non-parametrically, the time-varying nonlinear systems, including the nonlinear elastic and damping static force characteristics. Computation of the smooth time-varying characteristics, the congruent envelope and the congruent instantaneous frequency functions is based on the Hilbert vibration decomposition. According to the nonparametric identification result, the precise mathematic model of the time-varying vibration system can be established. The effectiveness of the developed identification method in the presence of time-varying parameters is investigated for the harmonically modulated Duffing differential equations of motion. The identified time-varying force function practically does not differ from the initial stiffness force modulation. References [1] C.S. Huang, S.L. Hung, W.C. Su, C.L. Wu, Identification of time-variant modal parameters using time-varying autoregressive with exogenous input and low-order polynomial function, Comput. Aided Civ. Infrastruct. Eng. 24 (7) (2009) 470–491. [2] X. Xu, W.J. Staszewski, X. Xu, Z.Y. Shi, S. Fassois, Identification of time-variant systems using wavelet analysis of force and acceleration response signals, in: Proceedings of the 4th International Operational Modal Analysis Conference (IOMAC) Istanbul, Turkey, 9–11 May, 2011. [3] Biswajit Basu, Satish Nagarajaiah, Arunasis Chakraborty, Online identification of linear time-varying stiffness of structural systems by wavelet analysis, Struct. Health Monit. (SAGE) 7 (1) (2008) 21–36.
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
77
[4] Z.Y. Shi, S.S. Law, X. Xu, Identification of linear time-varying MDOF dynamic systems from forced excitation using Hilbert transform and EMD method, J Sound Vib. 321 (2009) 572–589. [5] Young S. Lee, Stylianos Tsakirtzis, Alexander F. Vakakis, Lawrence A. Bergman, D. Michael McFarland, A time-domain nonlinear system identification method based on multiscale dynamic partitions, Meccanica 46 (2011) 625–649, http://dx.doi.org/10.1007/s11012-010-9327-7. [6] P.F. Pai, A.N. Palazotto, HHT-based nonlinear signal processing method for parametric and non-parametric identification of dynamical systems, Int. J. Mech. Sci. 50 (12) (2008) 1619–1635. [7] M. Feldman, Non-linear system vibration analysis using Hilbert transform—I. Free vibration analysis method FREEVIB, Mech. Syst. Signal Process. 8 (2) (1994) 119–127. [8] M. Feldman, Hilbert Transform Application in Mechanical Vibration, John Wiley and Sons Ltd., Chichester, United Kingdom, 2011. [9] M. Feldman, Time-varying vibration decomposition and analysis based on the Hilbert transform, J. Sound Vib. 2006 (2006) 518–530. [10] M. Feldman, Time-varying and non-linear dynamical system identification using the Hilbert transform, in: Proceedings of the 20th Biennial Conference on Mechanical Vibration and Noise, ASME, Long Beach, CA, 2005.
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp
Hilbert transform methods for nonparametric identification of nonlinear time varying vibration systems Michael Feldman n Faculty of Mechanical Engineering, Technion—Israel Institute of Technology, Technion City, 32000 Haifa, Israel
a r t i c l e in f o
abstract
Article history: Received 2 January 2012 Received in revised form 31 August 2012 Accepted 2 September 2012 Available online 26 September 2012
The objective of this paper is to demonstrate a Hilbert transform (HT) method for identification of mechanical time-varying vibration systems under free and forced vibration regimes. This special kind of non-stationarity arises in experimental vibration analysis and in engineering practice. The method is based on the HT of input/output signals in a time domain to extract instantaneous dynamic structure characteristics, such as natural frequencies, stiffness, damping, and their variations in time. The HT assigns a complementary imaginary part to a given real signal part, or vice versa, by shifting each component of the signal by a quarter of a period. Thus, the HT pair provides a method for determining the instantaneous amplitude and the instantaneous frequency of a signal. For general non-stationary vibration signals, the analytic signal method does a good job of simultaneous time-frequency localization of the main signal components. The paper focuses on HT signal processing techniques and identifies three groups of dynamics time-varying SDOF systems: slow varying quasi-periodic modulation of stiffness under free and forced vibration, slow varying quasi-periodic modulation of nonlinear stiffness under free vibration, and fast inter-wave parametric stiffness modulation (Mathieu equation). & 2012 Elsevier Ltd. All rights reserved.
Keywords: Non-linear time-varying system Hilbert transform Nonparametric identification Envelope Instantaneous frequency
1. Introduction Typically, a vibration system concerned with a mechanical construction having mass, stiffness and damping elements is a physical structure that, when taking an impulse or another input force, usually initiates an output natural vibration motion. If the system has a time-varying mass or stiffness parameter, the natural frequency will depend decisively on the variation. If the system also has a time-varying damping, the solution amplitude will react in situ with the damping variation. Such frequency and amplitude variations can be clearly detected by any suitable signal processing technique in the time or frequency domain. The appropriate signal processing is mostly based on the extracted features of the level and frequency content of a varying signal in the form of the envelope and the instantaneous frequency. Thus, for example, a moving least-squares technique and a Kalman filter allow time-varying frequency tracking, while the traditional Fourier transform and wavelets technique can identify online variation in the natural frequencies and mode shapes of a vibration system [1–3].
n
Tel.: þ 972 4829 2106; fax: þ 972 4829 5711. E-mail address: [email protected]
0888-3270/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ymssp.2012.09.003
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
67
If identification of the varying dynamic system requires additional a priori information about the dynamic system and its vibration model, this identification is called parametric, and its aim is to determine the parameters of a known model. A dynamic system model is commonly created on the basis of knowledge of a mechanical construction and a functional operation of a technical object [4]. In some other more sophisticated cases, a nonparametric analysis may be required to estimate the parameter variation of the unknown physical equation of the dynamic system [5,6]. Nonparametric system identification of the modal parameters, such as frequencies, damping, and mode shapes as functions of the physical inertia, stiffness, and damping properties, is based only on analysis of the system behavior. Recent works in the area of time domain representations of vibration based on the HT show great promise for applications in nonparametric vibration system identification. An outstanding work by N. Huang, known as Empirical Mode Decomposition, adaptively decomposes a signal into the simplest intrinsic oscillatory modes (components) in the first stage. Then, in the second stage, each decomposed component forms a corresponding set of instantaneous amplitudes and frequencies called the Hilbert spectrum that can be used to obtain time-frequency vibration representation (characterization). Every significant decomposed elementary component holds exact amplitude and frequency information contained in the original signal, and the HT method makes it sensitive to any modal shift. Practically all publications dedicated to HT time domain nonlinear vibration identification with envelope and instantaneous frequency use the approach suggested in [7]. This FREEVIB approach introduces a vibration solution as a new complex variable which, when substituted into an initial differential equation, directly determines the system’s natural frequency and damping parameters. The proposed methods for identifying instantaneous modal parameters (natural frequencies, damping characteristics and their dependencies on a vibration amplitude and frequency) prove to be very simple and effective. The HT-based signal processing in parallel with other time-frequency representations in the form of Wavelet transform, Gabor Transform, or Wigner–Ville distribution provides a powerful tool for vibration analysis and identification of vibration systems [8]. Of particular interest and importance is the use of the HT to interpret nonlinear system motion. The measured instantaneous characteristics, such as the instantaneous frequency or the amplitude of the nonlinear solution, take an unusual fast oscillation (intrawave modulation) form. In the case of supplementary parameter variation, both instantaneous functions of each vibration mode – the envelope and the instantaneous frequency – become complicated modulated functions: they reflect the slow time variation of the stiffness plus fast intrawave modulations due to nonlinearity [8]. Therefore, only a smeared or blurred form of the instantaneous frequency will be obtained through direct frequency monitoring. The main purpose of this paper is to demonstrate the method of nonparametric identification of a vibration system when the linear time-varying system description fails because the non-linearity is too severe and a very accurate model is needed. In this regard, not only must the presence of a parameter variation and nonlinearity be detected, but also an adequate and readily interpretable system model must be identified. It has been suggested that HT processing involving Hilbert vibration decomposition and congruent functions can identify the time-frequency distribution of a signal with great accuracy. The digital HT signal processing the same as every other approach has some inherent limitations. The method is based on the multipoint FIR digital Hilbert transformer and therefore it requires significantly longer initial data to neglect the filter end effects. The method has a good but an ultimate frequency resolution required for separation of closely spaced harmonics. It allows separation the finite small number of valued components, which does not exceed 5–6 components. It is very sensitive to measured instrumental noise. The HT decomposition method is not effective for the separation of nonoscillation types of motion, such as random, impulse, slow aperiodic signals. 2. Hilbert transform-based identification and analysis of time varying vibration systems 2.1. The Hilbert transform in time domain approach. The HT assigns a complementary imaginary part x~ to a given real signal part x, or vice versa, by shifting each component of the signal by a quarter of period x~ ¼ H½x. Thus, the HT pair provides a method for determining the instantaneous qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ amplitude AðtÞ ¼ 7 x2 ðtÞ þ x~ 2 ðtÞ and the instantaneous frequency of a signal oðtÞ ¼ d½arctan ðxðtÞ=xðtÞÞ=dt. For a general non-stationary vibration signal, the analytic signal method does a better job of simultaneous timefrequency localizing of the main signal components [8]. The analytic signal method is equally applicable to deterministic and random processes, although, generally speaking, it does not divide them into two separate groups. Hence, it enables us to investigate any oscillating time function from a more general point of view. The method is also good for solving problems concerning analysis of stationary and non-stationary vibrations, as well as narrow- and/or wideband multicomponent signals. It also allows precise analysis of the transformation and dissipation of vibration energy, and vibration effects on machine durability. 2.2. Advantages of Hilbert transform time domain identification Taking analytic signal representation into account enables one to consider a vibration process, at any moment in time, as a quasi-harmonic oscillation, amplitude- and frequency-modulated by time-varying functions A(t) and
68
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
oðtÞ : xðtÞ ¼ AðtÞcos
Rt
0 oðtÞdt. The instantaneous parameters are functions of time and can be estimated at any point of the vibration signal. The total number of points that map the vibration is much greater than the number of peak points of the signal. Therefore the HT provides the best available time resolution of the amplitude and other estimations of vibration analysis. It also opens the way for averaging and for other statistical processing procedures, making vibration analysis more precise [8]. Classically, frequency resolution means the smallest difference or granularity with which one can distinguish details in the frequency domain of a signal. For any other digital method the spectral resolution is a discrete and terminal value that cannot be less than the ratio between the sampling frequency and the number of measured points. In contrast, the HT calculates the frequency as the rate of phase change with time, and the instantaneous frequency resolution can be as small as desired depending on the accuracy of the signal phase and time measurement. Practically all suggested HT computational algorithms are simple and fast, thus allowing online (real-time) signal analysis and identification. The HT methods are suitable for analysis and identification of both linear and nonlinear vibration systems; they produce adaptive decomposition of complicated multicomponent vibration for elementary inherent components [9]. With the help of the HT, researchers get new features, such as the envelope and the instantaneous frequency of each component and the instantaneous phase relations between time-varying components. With regard to system identification, HT methods allow estimating not only the varying stiffness but also the varying damping in vibration constructions.
3. Quasi-linear slow time varying system identification Two methods of dynamic system studying based on the HT have recently been proposed for free and forced vibration analysis, when the input signal is an impulse, shock or sweeping harmonic. The methods are suitable for testing linear and non-linear systems and for instantaneous modal parameter identification, including determining the concrete type of nonlinear spring and damping characteristics for each mode of a vibratory system. The time-varying structural system encountered in civil and mechanical engineering can be described by the following _ and K(t)x are respectively equation of motion when a SDOF system is considered: mðtÞx€ þ CðtÞx_ þ KðtÞx ¼ 0 where m(t), CðtÞx, _ and the restoring force as a function of the displacement x. the varying mass, the frictional force as a function of velocity x, The corresponding second-order differential equation of the damped system for unit mass will then take the following form: x€ þ 2hðtÞx_ þ o20 ðtÞx ¼ 0, where the term hðtÞx_ ¼ CðtÞ=2mðtÞ represents the damping force per unit mass as a function of the velocity, h(t) is the damping coefficient, and the term o20 ðtÞ ¼ KðtÞ=mðtÞ represents the time varying undamped natural frequency squared. The instantaneous natural frequency of the solution according to the FREEVIB method is estimated by the following formula [7]: 2
€ þ2A_ =A2 þ A_ o _ =Ao o20 ðtÞ ¼ o2 A=A
ð1Þ
Here A is an envelope, and o are instantaneous frequency functions of the solution. But for many real engineering structures that are under force excitation, it is useful to consider a method of modal analysis of non-linear systems with input signal excitation. The modern technique for non-linear system investigation based on the HT and called FORCEVIB [8] enables us to identify instantaneous system modal parameters (natural frequencies, damping characteristics and their dependencies on vibration amplitude and frequency) during different kinds of excitation of the dynamic system. It allows direct extraction of the linear and non-linear system parameters from the measured time signal of input and output. The HT method of forced vibration x€ þ 2hðtÞx_ þ o2x ðtÞx ¼ zðtÞ determines instantaneous modal parameters even if the input signal is a fast sweeping frequency quasi-harmonic signal: _ omA=A € þ 2A_ 2 =A2 þ A_ o _ =Ao, where o0(t) is an instantaneous modal frequency, h(t) is an o20 ðtÞ ¼ o2 þ aðtÞ=m þ bðtÞA=A instantaneous modal damping coefficient, o and A are instantaneous frequency and envelope of the vibration with their _ A, € o _ Þ, and aðtÞ ¼ ReðZ=XÞ, bðtÞ ¼ Im ðZ=XÞ are real and imaginary parts of the input–output first and second derivatives ðA, signal ratio. Below we consider some simulated examples of identification of time-varying differential equations of motion performed in SIMULINK (MATLAB) with the permanent step value ODE4 (Runge–Kutta) solver.
3.1. Free vibration regime As a first example of a quasi-linear time-varying system, we choose a SDOF vibration model with the following square wave modulated stiffness: x€ þ0:05x_ þ ½1þ 0:6sgnðcos 0:15tÞx ¼ 0; x0 ¼ 1. The free vibration solution for nonzero initial condition and the sampling frequency 13 Hz are shown in Fig. 1. The variation of the solution and its envelope reflect the time modulated stiffness. The identified instantaneous natural frequency squared function (the instantaneous modal stiffness) calculated according to (1) accurately repeats the initial square wave stiffness (Fig. 1b). The identified instantaneous modal stiffness allows presenting the static force characteristic k(x) (Fig. 2), which for the time-varying case twists around the initial linear force characteristic.
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
69
1 vibration envelope
x
0.5
0
−0.5
0
20
40
60
80
100
120
2.5
140
160
180
200
Instantaneous frequency squared stiffness
Stiffness
2 1.5 1 0.5 0 0
20
40
60
80
100
120
140
160
180
200
Time, s Fig. 1. Linear system: the free vibration (a), the stiffness variation (b).
1.5
time−invariant stiffness time−varying stiffness
1
Spring force
0.5
0
−0.5
−1
−1.5 −1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
Displacement Fig. 2. Linear system: the stiffness force characteristics.
3.2. Forced vibration regime As the next example, let us consider a structure that describes a slow harmonic modulated elastic force of the trivial dynamics system under external harmonics excitation with a fixed frequency: x€ þ0:05x_ þ ð1 þ 0:5cos 0:03tÞx ¼ cos 2p0:16t (Fig. 3).
70
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
Excitation
1 0
Displacement
−1 0
100
200
0
100
200
300
400
500
300
400
500
20 0 −20 Time, s Solution (Displacement) Power Spectrum Magnitude
Power Spectrum Magnitude
Excitation
0.2 0.15 0.1 0.05 0
6 5 4 3 2 1 0
0
0.1 0.2 Frequency, Hz
0.3
0
0.1 0.2 Frequency, Hz
0.3
Fig. 3. Forced excited system: the excitation (a), the displacement (b), the excitation spectrum (c), and the displacement spectrum (d).
Displacement and Envelope 20
10
0
−10
−20 50
100
150
200
250
300
350
400
450
500
550
400
450
500
550
Instantaneous Frequencies
Frequency, Hz
0.2 0.18 0.16 0.14 0.12 50
100
150
200
250
300
350
Time, s Fig. 4. Forced excited system: the displacement and the envelope (a) and the instantaneous frequencies (b).
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
71
The generated vibration has a rather complicated form over time and also in the frequency domain (Fig. 4). But HT identification restores this modulation in accurate detail [10]. Thus, Fig. 4b shows that the identified instantaneous natural frequency of the vibration (bold line) completely coincides pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with the varying initial elastic force modulation function f ¼ 1 þ 0:5cos 0:03t =2p (dashed line). These two examples demonstrate that the HT method enables identification of a slow modulation of the linear system parameters, in both cases the free and the forced vibration regimes. Therefore the HT-based methods proposed for identifying instantaneous modal parameters prove to be very simple and effective [10]. 4. Nonlinear time varying system identification 4.1. Nonlinear interwave system representation The nonlinear restoring force k(x) in the second-order differential equation of a vibration system can be recast into multiplication of a varying non-linear natural frequency o20 ðtÞ and system solution x : x€ þ 2h0 ðtÞx_ þ o20 ðtÞx ¼ 0. To apply the multiplication property of the HT, let us assume that the new varying nonlinear natural frequency squared could be grouped into two different parts: o2 ðtÞ ¼ o20 ðtÞ þ o21 ðtÞ. The first part o20 ðtÞ is much slower, and the second component o21 ðtÞ is faster than the system solution, so the equation of motion takes a new form: x€ þ ½o20 ðtÞ þ o21 ðtÞxðtÞ ¼ 0. Now, according to the multiplication property of the HT for overlapping functions [8], we apply the HT to both sides of the last equation ~ þo ~ 21 ðtÞxðtÞ ¼ 0, where a tilde sign indicates the HT conjugate variable. Multiplying each side of the obtained HT x~€ þ o20 ðtÞxðtÞ equation by i and adding it to the corresponding sides of the initial equation yields a new form of a differential equation of ~ 21 Þx ¼ 0, where X is a complex solution in an analytic signal form: motion in the analytic signal form X€ þ o20 X þ ðo21 þio ~ This complex equation can be transformed into a more traditional and accepted form [8] X ¼ x þ ix. X€ þ idX þ o2 X ¼ 0
ð2Þ
2 2 ~ 21 xxÞ=A ~ ~ 21 x2 o21 xxÞ=A ~ where o2 ¼ o20 þ ððo21 x2 þ o Þ; d ¼ ððo Þ. Here o2 is a varying instantaneous natural frequency and d is a fast varying instantaneous fictitious friction parameter. It should be pointed out that this equation is not a real equation of nonlinear motion; however, it produces the same varying solution of the system. It is simply an artificial fictitious equation that yields the same nonlinear primary vibration solution. The obtained varying natural frequency o2 consists of a slow o20 and a fast component. The fictitious friction parameter d, on the contrary, consists only of a fast varying component. This means that the slow part of the natural frequency forms an average period of vibration, but the fast varying fictitious friction force (induced by the zero average value) does not affect the real average friction force. It should be noted that the inherent non-linear restoring force causes fast oscillation of both the instantaneous frequency and the fictitious instantaneous damping coefficient, even in a conservative non-linear system.
4.2. Congruent modal parameters When the instantaneous modal parameters of a nonlinear system vary with time much more rapidly than the system oscillation solution, it is desirable to squeeze these fast variations into smooth initial nonlinear functions. Such precise nonlinear identification paves the way for tracing the slow time-varying parameters of nonlinear systems. Precise identification of the smooth initial nonlinear static stiffness force is based on the displacement and the natural frequency functions corresponding to the potential energy maxima. Such functions were defined as a congruent envelope and a congruent instantaneous natural frequency [8]. The congruent envelope is the envelope of the envelope function; it is equal to the algebraic sum of the solution superharmonic envelopes Ac ðtÞ ¼
N X
Al ðtÞcosfl ðtÞ
ð3Þ
l¼1
where Al(t) is the envelope of the solution l order harmonic, fl(t) is the phase angle between the primary and the l order harmonic. For example, if all high order superharmonics are congruently in phase with the primary solution, the congruent envelope is equal to the sum of harmonic envelopes. If superharmonics are in different phase relations with extreme upper and lower points of the largest energy component, this governs the symmetry – or asymmetry – of the multicomponent signal. It is clear that estimation of the congruent envelope Ac(t) requires a preliminary time domain decomposition of the system solution into a sum of high order harmonics. In turn, the varying intrawave instantaneous natural frequency o0(t) can also be decomposed into a sum of high order synchronous components. So the congruent natural frequency will be
oc ðtÞ ¼
N X
o0l ðtÞcosfol ðtÞ
ð4Þ
l¼1
estimated as an algebraic sum of the superharmonics of the instantaneous natural frequency. Both the congruent envelope and the natural frequency, being smooth functions, define the exact static stiffness force of a nonlinear system for every
72
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
moment in time: o2c Ac . In practice, the accuracy of the last expression depends on the total number N of considered synchronous vibration components. For symmetric force characteristics, the initial smooth static force function can be written in the form: ( 2 oc Ac , x 40 kðxÞ ¼ , o2c Ac , x o 0 where o2c is the congruent modal frequency squared, and Ac is the congruent envelope of the displacement. Congruent functions of the asymmetrical vibration systems are split into two different parts, separated according to the positive and the negative displacements: ( ( ocp ðtÞ, if x 4 0 Acp ðtÞ, if x 4 0 Ac ðtÞ ¼ ; oc ðtÞ ¼ : ocn ðtÞ, if x r0 Acn ðtÞ, if x r0 Therefore, the nonparametric identification of the asymmetric vibration systems considers together these two solution parts in the time domain. As a result of the asymmetric signal, the decomposition will have two sets of envelopes and instantaneous frequencies, each for its own positive and negative part of the displacement. Then, the separated congruent functions construct two different static force characteristics—separately for the positive and the negative signal parts ( 2 ocp Acp , x 40 kðxÞ ¼ ð5Þ o2cn Acn , x r 0 To estimate the congruent functions, we need to produce a decomposition of the signal first and then to construct an algebraic summation of envelopes of all the decomposed components. Such a successive signal decomposition (disassembling) and the subsequent summation (re-assembling) of the component envelopes generates desirable congruent functions as slow functions of time suitable for a precise nonparametric identification.
4.3. Nonlinear system slow stiffness variation. free vibration regime For the nonlinear time-varying system example, we can construct a vibration motion by forming a combination of the following parameters: x€ þ 0:05x_ þð1 þ0:6cos 0:15tÞð1 þ13x2 Þx ¼ 0, x0 ¼ 2 with the sampling frequency 13 Hz (Fig. 5). The corresponding solution is a free vibration of the nonlinear Duffing model whose stiffness is slowly modulated by the harmonic function. Without such slow harmonic modulation, the solution would be a classic Duffing nonlinear solution. 1.8 1.7
x
1.6 1.5 1.4 1.3 1.2 20
22
24
26
28
30 vibration envelope congruent envelope
1.5 1
x
0.5 0 −0.5 −1 −1.5 0
20
40
60
80
100
120
140
Time, s Fig. 5. Free vibration of the Duffing system.
160
180
200
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
73
Both instantaneous functions of the solution – the envelope (Fig. 5) and the instantaneous frequency (Fig. 6a) – are complicated modulated functions; they reflect the slow time variation of the stiffness plus fast intrawave modulations due to cubic nonlinearity [8]. Therefore, only a smeared or blurred form of the instantaneous frequency (Fig. 6b) is available for direct frequency monitoring. A recently developed precise nonlinear identification based on the congruent frequency function (4) squeezes the smearing and blurring spots in a single clean line that determines the real nonlinear stiffness as a function of time (Fig. 7a). By applying the congruent frequency, we eliminate the fast intrawave modulations, though the desired slow modulation will still be present in the varying stiffness function (Fig. 7a).
Fig. 6. Instantaneous frequency of Duffing system: the IF and the congruent frequency (a) and the spectrogram (b).
35 congruent stiffness smoth congruent stiffness
Stiffness, (rad/s)2
30 25 20 15 10 5 0 0
20
40
60
80
100
120
140
160
180
200
Stiffness variation
2.5 initial stiffness variation identified stiffness variation
2 1.5 1 0.5 0 0
20
40
60
80
100
120
140
160
180
200
Time, s Fig. 7. Duffing system frequency and stiffness: the congruent and the average frequencies (a) and the varying stiffness (b).
74
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
100 time−invariant stiffness time−varying stiffness
80
60
40
Spring force
20
0
−20
−40
−60
−80
−100 −2
−1.5
−1
−0.5
0
0.5
1
1.5
2
Displacement Fig. 8. Duffing varying system: the stiffness force characteristics.
Let assume that the monotonic form of the varying stiffness is an inherent feature of the nonlinear system that is not related to the desired slow varying modulation. Now we can remove the monotonic form by using a standard polynomial fitting technique in a least squares sense. The remain wave-like instantaneous stiffness function is plotted on Fig. 7b together with the initial harmonic modulation function 1 þ0.6 cos 0.15t. As can be seen, both the initial and the identified time-varying stiffness functions are in good agreement. The identified time-varying stiffness multiplied by the congruent envelope function constructs the nonlinear stiffness force characteristics, which in our case takes the form of a spiral line twisted around the initial nonlinear cubic force characteristic (Fig. 8). 4.4. Fast parametric excitation regime If a slow varying stiffness parameter is the case, the frequency of the solution will simply follow and reflect this slow variation. For an artificial fast parameter variation, in contrast, the behavior of the system solution can be quite different in principle. The fast stiffness modulation with the doubled (relative to the natural frequency) and faster multiple frequencies could even excite the system parametric resonances. Time domain identification will require a priori knowledge about the existence of the fast varying parameter and about the vibration system structure. This means that identification of the fast time-varying parameters is not a purely non-parametric problem. As an example, we describe the process of identifying a vibration system with parametric excitation in the form of linear differential equation with variable stiffness coefficient (Mathieu equation): x€ þ2h0 x_ þ o20 ð1 þ esinof ast tÞx ¼ 0, where e is the parametric modulation depth. We choose a stiffness coefficient that alternates harmonically twice over the period of the excited oscillations ofast ¼ 2o0. To identify the equation with parametric excitation, we use the non-linear system representation (2) and the R quasi-harmonic form of the solution x ¼ AðtÞcos odt. After the algebraic transform, the initial parametric equation takes a new form x€ þ ½2h0 0:5eo0 ð1cos 2o0 tÞx_ þ o20 ð1 þ0:5esin 2o0 tÞx ¼ 0
ð6Þ
This form now includes the fast varying stiffness factor 1 þ0:5esin2o0 t and also the fast varying damping component 0:5eo0 cos 2o0 t. The objective of the identification could be the unknown parametric modulation depth e. The obtained form of the equation illustrates the well-known fact that parametric excitation essentially affects system damping, reducing the real damping up to zero. The critical condition 2h0 o0.5eo0 causes loss of stability and appearance of unlimited amplitude vibration. For our Matlab simulation, we chose an initial numeric model describing parametric oscillation in the form: x€ þ 0:1x_ þ ð1 þ0:3sin 2tÞx ¼ 0, x0 ¼ 1, where o0 ¼1, e ¼0.3, 2h¼0.1. The obtained parametrically unstable excited vibration of
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
75
the system, together with the vibration envelope, is shown in Fig. 9a. The executed vibration has a frequency that is two times smaller in comparison with the parametric excitation frequency. After using the FREEVIB identification, we get the instantaneous natural frequency squared shown in Fig. 9b, together with the initial modulated stiffness function 1þ0.3 sin 2t. For the case of fast parameter modulation, the estimated instantaneous natural frequency level is equal only to half of the initial elastic parametric variation. This means that the direct time-frequency monitoring of the solution does not provide the correct value of the real stiffness variation.
100
vibration envelope
50
x
0
−50
−100 0
40
60
80
100
120
140
160
180
200
100
120
140
160
180
200
instantaneous frequency squared initial stiffness variation
1.4
Stiffness
20
1.2 1 0.8 0.6 0
20
40
60
80
Time, s Fig. 9. Parametrically excited system: the vibration (a) and the instantaneous frequency and the stiffness variation (b).
Stiffness
1.5
identified variation initial variation
1
0.5
Damping coefficiemt
0
20
40
60
80
100
120
140
160
180
200
identified variation theoretical variation initial constant value
0.2 0.1 0 −0.1 −0.2 0
20
40
60
80
100
120
140
160
180
200
Time, s Fig. 10. Parametrically excited system: the identified and the initial stiffness variation (a) and the identified and the theoretical damping variation (b).
76
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
time−invariant stiffness time−varying stiffness
100
Spring force
50
0
−50
−100
−100
−80
−60
−40
−20
0
20
40
60
80
100
Displacement Fig. 11. Parametrically excited system: stiffness force characteristics.
To identify the unknown parametric modulation depth e, we should consider the theoretical fast varying stiffness factor (6) 1 þ0:5esin2o0 t. The resultant fast varying elastic component for the estimated modulation depth e ¼0.1 is plotted together with the initial modulated stiffness function (Fig. 10a). As can be seen, both the estimated and the initial functions practically coincide. The discovered modulation depth e ¼0.1 can be also used to compare the theoretical damping variation (6) and the real time-varying damping function obtained through the FREEVIB method (Fig. 10b). Both these damping functions are in very good agreement, and they lie in the negative half-plane below the constant initial damping coefficient line (2h0 ¼ 0.1). The observed fast oscillation stiffness force can be represented as a function of the displacement in the form of traditional force characteristics (Fig. 11). Due to the fast parametric excitation, the shown force characteristics take the form of a spiral repeatedly wrapped around the initial linear force characteristic. 5. Conclusions The direct way to identify time-varying instantaneous modal parameters of nonlinear vibration systems is to use HT-based vibration decomposition and the appropriate HT time domain signal processing. The modern HT signal decomposition divides the real multicomponent solution into a principal nonstationary harmonic and a number of separate high-frequency superharmonics. A solution of nonlinear time-varying vibration systems can also include a slow variation component. By considering the high superharmonics, the congruent smooth envelope and the natural frequency, we are able to identify, precisely and non-parametrically, the time-varying nonlinear systems, including the nonlinear elastic and damping static force characteristics. Computation of the smooth time-varying characteristics, the congruent envelope and the congruent instantaneous frequency functions is based on the Hilbert vibration decomposition. According to the nonparametric identification result, the precise mathematic model of the time-varying vibration system can be established. The effectiveness of the developed identification method in the presence of time-varying parameters is investigated for the harmonically modulated Duffing differential equations of motion. The identified time-varying force function practically does not differ from the initial stiffness force modulation. References [1] C.S. Huang, S.L. Hung, W.C. Su, C.L. Wu, Identification of time-variant modal parameters using time-varying autoregressive with exogenous input and low-order polynomial function, Comput. Aided Civ. Infrastruct. Eng. 24 (7) (2009) 470–491. [2] X. Xu, W.J. Staszewski, X. Xu, Z.Y. Shi, S. Fassois, Identification of time-variant systems using wavelet analysis of force and acceleration response signals, in: Proceedings of the 4th International Operational Modal Analysis Conference (IOMAC) Istanbul, Turkey, 9–11 May, 2011. [3] Biswajit Basu, Satish Nagarajaiah, Arunasis Chakraborty, Online identification of linear time-varying stiffness of structural systems by wavelet analysis, Struct. Health Monit. (SAGE) 7 (1) (2008) 21–36.
M. Feldman / Mechanical Systems and Signal Processing 47 (2014) 66–77
77
[4] Z.Y. Shi, S.S. Law, X. Xu, Identification of linear time-varying MDOF dynamic systems from forced excitation using Hilbert transform and EMD method, J Sound Vib. 321 (2009) 572–589. [5] Young S. Lee, Stylianos Tsakirtzis, Alexander F. Vakakis, Lawrence A. Bergman, D. Michael McFarland, A time-domain nonlinear system identification method based on multiscale dynamic partitions, Meccanica 46 (2011) 625–649, http://dx.doi.org/10.1007/s11012-010-9327-7. [6] P.F. Pai, A.N. Palazotto, HHT-based nonlinear signal processing method for parametric and non-parametric identification of dynamical systems, Int. J. Mech. Sci. 50 (12) (2008) 1619–1635. [7] M. Feldman, Non-linear system vibration analysis using Hilbert transform—I. Free vibration analysis method FREEVIB, Mech. Syst. Signal Process. 8 (2) (1994) 119–127. [8] M. Feldman, Hilbert Transform Application in Mechanical Vibration, John Wiley and Sons Ltd., Chichester, United Kingdom, 2011. [9] M. Feldman, Time-varying vibration decomposition and analysis based on the Hilbert transform, J. Sound Vib. 2006 (2006) 518–530. [10] M. Feldman, Time-varying and non-linear dynamical system identification using the Hilbert transform, in: Proceedings of the 20th Biennial Conference on Mechanical Vibration and Noise, ASME, Long Beach, CA, 2005.