- Email: [email protected]
Fourier spectral-based modal curvature analysis and its application to damage detection in beams
Mechanical Systems and Signal Processing 84 (2017) 763–781
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp
Fourier spectral-based modal curvature analysis and its application to damage detection in beams Zhi-Bo Yang a,b,n, Maciej Radzienski b, Pawel Kudela b, Wieslaw Ostachowicz b,c a b c
School of Mechanical Engineering, Xi’an Jiaotong University, 710049 Xi’an, PR China Institute of Fluid Flow Machinery, Polish Academy of Science, 80-231 Gdansk, Poland Faculty of Automotive and Construction Machinery, Warsaw University of Technology, 02-524 Warsaw, Poland
a r t i c l e i n f o
abstract
Article history: Received 15 February 2016 Received in revised form 10 May 2016 Accepted 2 July 2016 Available online 18 July 2016
In this paper, a simple Fourier spectral-based method is proposed to calculate the modal curvature (MC) of beams instead of the traditional central difference method. Based on the present method, damages in beam-like structures are localized. The present method provides an alternative selection to estimate MC in damage detection. There are two advantages of the present method. Firstly, the spectral calculation of spatial derivatives is conducted globally, which provides the suppression for noise. In addition, signal processing in the wavenumber domain provides an alternative choice for spatial filtering for mode shapes. Secondly, the proposed method provides a precise estimation of the MC which is related to original definition. With the absence of numerical derivative, the estimated results can be more stable and robust. Statistical analysis is conducted to show the effectiveness and noise immunity of the proposed method. In order to obtain the better identification, the MC calculated by the proposed method is employed as the input of continuous wavelet transform, and then the hybrid method is generated. The validations of the present method and comparison with the traditional central difference method are numerically and experimentally demonstrated. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Fourier transform Modal curvature Damage detection Beam Wavenumber filtering
1. Introduction Structural damage detection has been a research focus with increasing interest in mechanical, civil, aerospace and some other fields during the last few decades [1–9]. The modal curvature (MC) has been investigated in a number of studies [10– 14]. The MC is the curvature of mode shapes of structures. In beam-like structures, it can be expressed as:
w″( x) =
d2w ( x) M = − dx2 EI ( x)
(1)
where w′′(x) is the MC, w(x) is the displacement mode shape of a beam, M is the bending moment and EI(x) is the bending stiffness. A reduction of the bending stiffness induced by damage causes alteration in the MC. This change is directly related to damage location. Therefore, the MC is a theoretically feasible damage index in beam-like structures. However, only few sensor types can provide satisfactory direct measurement of the MC. Instead, some local estimation methods, such as the n
Corresponding author at: School of Mechanical Engineering, Xi’an Jiaotong University, 710049 Xi’an, PR China. E-mail addresses: [email protected], [email protected], [email protected] (Z.-B. Yang).
http://dx.doi.org/10.1016/j.ymssp.2016.07.005 0888-3270/& 2016 Elsevier Ltd. All rights reserved.
764
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Nomenclature
ε λ ψ Δ a ai cr e h i f ( ∙) k t wd w ( x) ^ ( k) w wx
the convergence error the random parameter the signal to noise ratio parameter the height of rectangular impulse the depth of crack the constant about Nuttall window the constant about damage the constant about convergence error the sampling interval imaginary number the dimensionless local compliance function wavenumber time the damaged mode shape the displacement mode shape of beam the Fourier transform of w ( x ) the discrete form of w (x )
^ w k w′( x ) w″( x ) wx′ wx″ yn(t) yn∘ ( t ) B E EI ( x ) H Kc L Lc Ll M N WNuttall ‖∙‖∞
^ (x ) the discrete form of w the first-order derivative of w (x ) the modal curvature the discrete form of w′( x ) the discrete form of w″( x ) the noisy response of measuring points n the noise-free response of measuring points n the width of beam the Young's modulus the bending stiffness the thickness of beam the stiffness of crack the length of beam crack location the width of rectangular impulse. the bending moment the number of grid point the Nuttall window the infinite norm
numerical difference techniques, are usually used in the MC approximations. Among these numerical difference techniques, the central difference approximation is usually employed in the MC estimation. For simplicity, the MC calculated by the central difference approximation is denoted as the CD-MC. The formulation of the CD-MC is:
w″( x) =
w ( x − h) − 2w ( x) + w ( x + h) h2
(2)
where h is the sampling interval. The application of the MC in damage detection can be traced back to Pandey's work in 1991 [10], in which the abrupt change of the MC is seen as the indication of damage. Following this work, some investigations have been conducted, and the MC is also used as a baseline-free damage index as presented in [12,15–17]. Although the effectiveness of Eq. (2) has been verified in these works, we are willing to point some possible deficiencies and validate them in this paper: (1) Ill-posed condition. False peaks in the CD-MC can be easily affected by noise. The problem of numerical difference (containing the central difference) is known to be ill-posed in the sense that a small measurement error/noise induces a large error in the computed derivatives. (2) Approximation order. Mode shape can be represented by trigonometric function. However, local estimation of mode shape is often done by a polynomial approximation. Thus, the approximation error is relatively large, which influences the accuracy of numerical difference further. (3) Noise suppression. Due to the local estimation, the capacity of the CD-MC for noise suppression is relatively weak. To solve these problems, some investigations have been done [13,16,18–21]. However, these studies approximate the MC by using the central difference. Therefore, the aim of this research work is to provide the MC estimation without the central difference approximation. Paying attention to problems (1–3), one can find that the Fourier spectral method is the natural choice for solving them by aid of the following properties: (1) Special process in derivative calculation. The difference of functions can be easily conducted via multiplication in wavenumber domain of the Fourier transform (FT), and reconstruction via the inverse Fourier transform (IFT). (2) The approximation function in Fourier spectral method is trigonometric function. (3) Noise suppression via averaging in the mapped domain. Fourier spectral-based methods are global, thus the parts of the noise can be suppressed in the wavenumber domain via special filters. Compared with the frequency-based method [22–24], modal shape-based method (containing MC) performs better in damage localization. But the modal shape-based method [25] is also argued by its sensitivity to small damage. Modal test often shows the low-frequency domain information compared with the high-frequency domain information provided by wave propagation-based method [26,27]. The identification of them can be complementary to each other, especially for large structures. The aim of this paper is to provide an alternative and simple method to estimate the MC used in damage detection. The following parts of this paper are organized as follows: Section 2 presents the formulas and implementation of the Fourier
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
765
spectral method for the MC estimation. Section 3 mainly focuses on the numerical comparison between the proposed method and the CD-MC method. In addition, some auxiliary methods are also introduced to improve the performance of the proposed method. Section 4 validates the proposed method via benchmark test. Finally, some discussions and conclusions are presented in Section 5.
2. Fourier spectral method for MC estimation 2.1. Theoretical basis Fourier spectral method is one of the most important techniques for numerical solution of partial differential equations. For the sake of completeness theoretical basis are explained in this section. The FT is used in signal analysis to map signals from the time domain to the frequency domain. Let's consider the mode shape w(x) as the function to be transformed by the FT and the IFT in the spatial and wavenumber domain: +∞
^ ( k) = w
∫−∞
w ( x) =
1 2π
e−ikxw ( x) dx +∞
∫−∞
(3)
^ ( k ) dk eikxw
(4)
where x ∈ and k ∈ . k is called wavenumber with the specific physical meaning. Paying attention to the distinguished characteristic of the FT:
dw ( x) 1 i = dx 2π
w′( x) =
+∞
∫−∞
^ ( k ) dk eikxkw
(5)
one can get the analytical form of the MC instead of Eq. (2):
w″( x) =
dw′( x) 1 = − dx 2π
+∞
∫−∞
^ ( k ) dk eikxk2w
(6)
However, Eq. (6) cannot be employed to estimate the MC directly because the analytical form of the mode shape w(x) is unknown in application. What is worse, the integration in Eq. (6) is defined on the interval [ −∞ , + ∞]. Given that the mode shape w(x) is discrete and bounded in numerical estimation and experiment, the discrete Fourier transform (DFT) is the natural choice for the MC calculation. Without the loss of generality, it can be assumed that w(x) is defined on the interval [0, 2π], and the interval is further divided by N 1 equal subintervals between N grid points. Let's define the length of each subinterval as h, then the related DFT synthesis can be written as:
^ = h∑ e−ikxw w k x x
wx =
1 ∑ eikxw^ k 2π k
(7)
(8)
^ are the discrete forms of w ( x ) and w ^ ( x ) respectively. Corresponding to (Eqs. (7) and 8), (Eqs. (5) and 6) are where wx and w k rewritten in the DFT frame as:
wx′ =
⎛ ⎞ h i∑ keikx ⎜⎜ ∑ e−ikxwx ⎟⎟ 2π k ⎝ x ⎠
wx″ = −
⎛ ⎞ h ∑ k2eikx ⎜⎜ ∑ e−ikxwx ⎟⎟ 2π k ⎝ x ⎠
(9)
(10)
In (Eqs. (9) and 10), the x and k are all bounded and discrete, without any technical restrictions. Based on Eq. (10), the MC can be calculated easily. This process can be seen as a trigonometric interpolation, thus Eq. (10) provides a semi-analytical method to estimate the MC. According to the characteristic of trigonometric functions, the problems (1, 2) are solved. Problem (3) will be solved by filtering in the wavenumber domain. 2.2. Damage models Damage in structures can be manifested by the singularity in mode shapes. Thus, it is necessary to investigate and classify singularities, based on which generally applicable wavenumber domain filters are provided. Before the discussion about wavenumber domain filter, the procedure for establishing damage models of beam is presented in this part of the
766
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 1. The influences of linear spring crack on the 1st–3rd mode shapes: crack features in spatial domain (a–c) and in wavenumber domain (d–f).
filter selection. In this paper, the numerical simulations are conducted using the cantilever Euler–Bernoulli beam model with the rectangular cross-section. The structure used in numerical test is a beam with the length L ¼535 mm, the width B ¼20 mm and the thickness H¼12 mm. The selection of material can be arbitrary, thus it is not highlighted here. Given that transverse vibration of the Euler–Bernoulli beam concerns only the bending vibrations of beam, the rotational spring constant is assumed to be dominant in the local flexibility matrix [1,3], and the stiffness of rotational spring is given by [9,28–30]:
Kc =
BH 2E 2
72π ( a/H ) f ( a/H )
(11)
where E is the Young's modulus, B is the width of beam, H is the depth of beam, and f(a/H) is the dimensionless local compliance function calculated as:
Fig. 2. The influences of rectangular model on the mode shape: (a) the damaged 1st mode shape, (b) damage feature in spatial domain, (c) damage feature in wavenumber domain.
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
⎛ a⎞ ⎛ a ⎞2 ⎛ a ⎞3 ⎛ a ⎞4 ⎛ a ⎞5 ⎛ a ⎞6 a f ⎜ ⎟ = 0.6384 − 1.035 + 3.7201 ⎜ ⎟ − 5.1773 ⎜ ⎟ + 7.553 ⎜ ⎟ − 7.332 ⎜ ⎟ + 2.4909 ⎜ ⎟ ⎝ H⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ H⎠ H H H H H
767
(12)
Based on (Eqs. (11) and 12), the influence of a crack on mode shape can be calculated easily by finite element method or analytical model. To illustrate effects caused by the crack model, a severe crack (a¼4.8 mm) is arranged at the location Lc ¼321 mm. The differential waveforms between cracked and intact mode shapes and the corresponding wavenumber domain representations are displayed in Fig. 1. Although crack features of 1st–3rd mode shapes are totally different in the spatial domain, the wavenumber domain representations of them are similar (belonging to the low wavenumber domain). Namely, one can conclude that the crack related features are concentrated in the low wavenumber domain via DFT. Second model, the rectangular model, is simpler compared with the linear spring model. By adding a constant component cr on the specific interval of the intact mode shape w, the damaged mode shape wd can be obtained. The formulation of this model is:
⎧ ⎪ wx + cr wd = ⎨ ⎪ ⎩ wx
( x = Lc ) ( for the others)
(13)
Fig. 2 shows the damage features with the parameter cr ¼0.02. As the damage features in every order mode shapes are the same, thus only the first mode shape is presented. Comparing the spatial domain representation with the wavenumber domain representation, one can find the feature of this damage model mainly concentrates in low wavenumber domain, similar to the linear spring model. The third model (zero-frequency component model) is formulated as:
⎧ x ∈ ( L c , L c + h, ... , L ) ⎪ wx + cr wd = ⎨ ⎪ x ∈ ( 0, L/N , ... , L c − h) ⎩ wx
(14)
In this model, the effect induced by damage is considered as a leap at the damage location Lc, and the area from damage location Lc to the free end L is also affected as show in Fig. 3. Similarly, damage features concentrates in low wavenumber domain. 2.3. Wavenumber domain filter Fig. 4 presents the representations of noise in spatial domain and wavenumber domain. The waveform shown in Fig. 4 (a) consists of the pseudo-random value generated by ‘rand’ in Matlab to simulate a white noise. We can see the noise energy distributes evenly in spatial domain and wavenumber domain. It is different with the features of the singularities mentioned above. Thus, a low-pass filter can be used in the wavenumber domain to suppress noise. In this paper, a Nuttall window WNuttall is applied in the wavenumber domain directly, which is given by:
WNuttall = a 0 − a1 cos
π ( 2hk + L ) L
+ a2 cos
2π ( 2hk + L ) L
− a3 cos
3π ( 2hk + L ) L
−
L L ≤k≤ 2h 2h
(15)
where a0 ¼ 0.3635819, a1 ¼0.4891775, a2 ¼0.1365995 and a3 ¼ 0.0106411. The comparison of the Nuttall window, the Blackman window and the Hann window is displayed in Fig. 5, in which one can observe that the Nuttall window attenuates faster than the other two windows. This is an appropriate property for noise suppression in high wavenumber domain referring to Figs. (1–4). The procedure of wavenumber domain filtering is
Fig. 3. The influences of zero-frequency component model on the mode shape: (a) the damaged 1st mode shape, (b) damage feature in spatial domain, (c) damage feature in wavenumber domain.
768
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 4. The features of noise: (a) the contaminated 1st mode shape, (b) noise in the spatial domain, (c) noise in the wavenumber domain.
Fig. 5. Comparison of window functions.
formulated as:
^ i+1 = W ^i w Nuttall *wx x
(16)
The superscript depicts the step in iteration. Given that the required times of the wavenumber domain filtering is un-
^1 known, the iteration shown in Fig. 6 is required. In Fig. 6, ε depicts the convergence error. We suggest to choose ε = e w x
‖ ‖, ∞
where e is set as 1% in this paper (the choice is flexible). To present the effectiveness of the iteration filtering, the intact mode shape, which is contaminated by the same noise level (shown in Fig. 4), before and after 5 times filtering processes is presented via spatial and wavenumber representations in Fig. 7. Noise is considerably reduced by the proposed procedure in the wavenumber domain, and the reconstructed mode shape is smoother. Further investigations and validations about the iteration will be presented in the following section.
3. Method verification 3.1. Damage scenarios The performance of the proposed method (denoted as the DFT-MC) for damage identification is investigated for a wide spectrum of damage scenarios (Table 1), and the CD-MC is employed as the numerical benchmark of the DFT-MC. In a progressive manner in terms of damage identification for different scenarios and the noise immunity tests, comprehensive comparisons are conducted between the DFT-MC and the CD-MC. In Table 1, the spectrum of damage scenarios are provided following damage models mentioned in the above section. For the linear spring model, the damage parameter vector is defined as ( a, L c ). For the rectangular model, the damage parameter vector is defined as ( Δ, Ll, L c ), where Δ is the height of rectangular impulse and Ll is the width of rectangular impulse.
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
769
w
Fig. 6. The flow chart of the wavenumber domain filtering.
L
k
Fig. 7. Illustration for the iteration result: (a) the contaminated 1st mode shape before and after filtering, (b) the noise before and after filtering.
The damage parameter vector for the zero-frequency component model is defined as ( Δ, L c ), where Δ is the height of step change. Damage scenarios are organized as follow: (1) Scenarios 1–3 are used to verify the identification capacity of DFT-MC for the three models, respectively. (2) Scenarios 4 and 5 are used to verify the capability of the present method for mixed damage identification with different damage severities. Scenarios 5 and 6 are used to verify the capability of the present method for mixed damage identification with different damage locations. (3) In scenarios 1–6, the 1st–3rd mode shapes are used in validations. The 1st, 5th and 10th mode shapes are used in scenario 7 for the purpose of comparison in low, middle and high mode shapes, respectively. (4) The further comparisons between CD-MC and DFT-MC are based on scenario 6 via noise immunity tests. (5) Scenario 8 is used for further comparisons shown in Section 3.4. It should be mentioned that cyclic continuations have been used in the following examples to avoid boundary effects, either for CD-MC or DFT-MC.
770
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Table 1 Damage scenarios used in validations. Damage scenario
Mode shape order
Linear spring model
( Δ,
( a,
[mm, mm, mm]
[mm, mm]
– (0.02, 15, 321) – (0.02, 15, 321) (0.01, 15, 321) (0.01, 15, 421) (0.01, 15, 421) –
– – (0.02, 321) (0.02, 421) (0.01, 421) (0.01, 221) (0.01, 221) –
Lc )
Rectangular model
Ll , L c )
Zero-frequency component model ( Δ, L c )
[mm, mm] 1 2 3 4 5 6 7 8
1–3 1–3 1–3 1–3 1–3 1–3 1, 5, 10 1–3
(2.4, 321) – – (2.4, 221) (1.2, 221) (1.2, 321) (1.2, 321) (0.6, 321)
Fig. 8. Scenarios 1–3: (a–c) the 1st–3rd DFT-MCs of scenarios 1, (d–f) the 1st–3rd DFT-MCs of scenarios 2, (g–i) the 1st–3rd DFT-MCs of scenarios 3. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
771
Fig. 9. Scenarios 4 and 5: (a–c) the 1st–3rd DFT-MCs of scenarios 4, (d–f) the 1st–3rd DFT-MCs of scenarios 5.
3.2. Noise-free test Fig. 8 presents the identifications related to scenarios 1–3 via DFT-MC. The present method identifies all the three kinds of damage clearly. The damage localization also agrees well with the actual locations (depicted by red triangles) The damage features displayed by modal curvature are different for these three damage models. This provides an evidence for distinguishing different type of damage. This group of simulations verifies the capacity of DFT-MC for the identification of different damages. Scenarios 4 and 5 combine three models in one mode shape, and the results are displayed in Fig. 9. The results agree with the actual conditions. It is also observed that singularity induced by the linear spring model increases with the order of mode shape. Either in scenarios 4 or 5, three singularities are all identified clearly. In higher modal curvatures, the related amplitudes are also distinct. This group of simulations verifies that DFT-MC is capable of detecting mixed damage. In scenario 6, the locations of damage are changed compared with scenario 5. However, one can see that the identification accuracy is not affected by this change (see in Fig. 10). Thus, the adaptability of the present method is validated to some extent. Fig. 11 presents a preliminary comparison between DFT-MC and CD-MC via scenario 7. This scenario is designed to verify and compare the identification capabilities of DFT-MC and CD-MC via higher mode shapes and damage mixture. These two
Fig. 10. Scenario 6: (a–c) the 1st–3rd DFT-MCs.
772
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 11. Scenario 7: (a–c) the 1st–3rd DFT-MCs, (d–f) the 1th–3rd CD-MCs.
methods obtain accurate identification for three damage models. The detected singularity via CD-MC in Fig. 11(d) at 0.321 m induced by the linear spring model is weaker than the one detected by DFT-MC in Fig. 11(a) at the same location. However, the conclusion that DFT-MC is better than CD-MC can be hardly obtained just from the above performance. Thus, further investigations are displayed via noise immunity test in the following part. 3.3. Noise immunity test Before the comparison in noise immunity test, a fact should be noticed: In the following comparisons, filtering is just applied in DFT-MC in wavenumber domain. Here, the filtering is emphasized in spatial domain, the filtering in the time domain during measurement is not included. That is to say: the results given by CD-MC is not filtered, thus comparison seems to be unfair for CD-MC in this view. However, the implementation of filtering is easier to perform in the wavenumber domain than in the spatial domain. The wavenumber domain filtering is one of the advantages of the present method. Of course, the spatial or modal filtering methods are available in the literatures [31–35], which provide the possibilities for the filtering of CD-MC. This problem is not investigated further in this paper. In this paper, noise is defined as the white Gaussian noise (WGN). The simulated WGN is added on the waveform directly
Fig. 12. (a–c) the 1st–3rd CD-MCs based on noisy mode shapes, ψ ¼0.02.
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
773
Fig. 13. DFT-MCs based on noisy mode shapes, (a) the 1st DFT-MC with ψ ¼0.02, (b) the 2nd DFT-MC with ψ ¼0.02, (c) the 3rd DFT-MC with ψ ¼0.04, (d) the 1st DFT-MC with ψ ¼0.03, (e) the 2nd DFT-MC with ψ ¼ 0.04, (f) the 3rd DFT-MC with ψ ¼0.05.
referring to [32]:
yn ( t ) = yn0 ( t ) + ψλ max ⎡⎣ yn0 ( t ) ⎤⎦
(15)
where yn(t) and ( t ) are the noisy and noise-free responses of measuring points n, λ is the random parameter, with its continuous distribution following a Gaussian distribution with zero mean and unitary standard deviation, ψ is the parameter which controls the signal to noise ratio. Fig. 12 shows the CD-MC results for scenario 6 (ψ ¼0.02). The noise amplified by the central difference severely impairs the profile of the modal curvature. From these figures, useful damage features can be hardly extracted. This example illustrates the susceptibility of CD-MC to noise. The identifications given by DFT-MC are presented in Fig. 13. The noise levels are assigned as ψ ¼ 0.02 and 0.03 for the 1st order DFT-MC, ψ ¼0.02 and 0.04 for the 2nd order DFT-MC, and ψ ¼0.04 and 0.05 for the 3rd order DFT-MC. According to the flow chart shown in Fig. 6, most of the locations of the three defects are represented clearly in each noise level. The identified damage locations agree well with the actual damage locations. It is observed in Fig. 13 that the higher order of mode shape causes lower noise susceptibility of the DFT-MCs. In addition, some false peaks (the high-wavenumber disturbances) also appear in Fig. 13, which disturb the identification of damage locations. Thus, the misjudgment and location error may emerge. In the following part, we will discuss the influences and the inhibition methods for the low/high-wavenumber disturbances by aid of statistical analysis and Monte Carlo simulation.
yn∘
3.4. Statistical analysis and post-processing method Disturbances are divided into two groups, a low- and a high-wavenumber group. The low-wavenumber disturbances are related to the waveforms of higher mode shapes, whose fluctuation will be expressed in DFT-MC too. Some post-processing damage detection methods are suggested to restrain this part of disturbance. 3.4.1. Statistical analysis Appearance of a noise may cause an error in identification and a probability of misjudgment. Influence of the noise can be hardly eliminated, but statistical analysis can help to minimize undesirable conditions. There is a requirement for the damage index, which helps to indicate the peak related to damage location among the results given by curvature methods. As most points in the measured system are related to the intact parts of structures, the
774
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 14. The illustrations of typical samples: (a) correct identification, (b) identification with locating error, (c) false dismissal, (d) failure. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
damage points become the outliers compared with the intact points distributed in a certain field. Selection of quantile to determine outlier is important: If the quantile is chosen too small, risk of erroneous judgment will sharply increase, and vice versa. An appropriate choice of quantile can be determined by the balance of risk and false dismissal probability. In this paper, quantile is selected as 3s [30]. In this part, the noise-polluted scenarios 1 and 8 are employed as the basic cases to evaluate the error-ratio/accuracy of the DFT-MC method according to Monte Carlo simulation. Some typical samples are displayed in Fig. 14, in which the filled triangles on the bottom of these figures indicate the actual location of damage, and the unfilled triangles plotted on the curves indicate the identified location (if any). The red lines parallel to the horizontal axes are the upper and lower bounds determined by 3s. The following terms should be highlighted before the simulation: (1) It can be observed in Fig. 14(a, b, and d) that the peak beyond the bounds will be deemed as damage location. If there is no point beyond the bounds, the structure will be deemed intact as show in Fig. 14(c). (2) It is noticed in Fig. 14(a) that the points beyond the bounds determine a damaged interval instead of a specific location. In this paper, we employ the peak location in this interval as the damage location instead of the center of the interval. (3) The identification as shown in Fig. 14(b) is defined as damage localization error because the difference between identified and actual location is smaller than 1% of the total length of the beam. Other identification will be seen as incorrect damage localization which means that the damage identification procedure failed. This case is indicated within the paper as failure. (4) Other complex conditions can be deemed as the combinations of the typical cases shown in Fig. 14. The Monte Carlo simulation consists of some sub-simulations: scenarios 1 with different noise levels and scenarios 8 with different noise levels. Each sub-simulation contains 1000 samples. According to the Chebyshev law of large number, the capacity as 1000 in simulation is enough to fulfill the requirement that sample distribution should cover the entire solution field. By choosing the different noise levels, some typical evaluations are listed in Table 2 and shown in Fig. 15. It should be mentioned that different noise levels are employed for the 1st–3rd order DFT-MCs to obtain similar accuracy ratio, false dismissal ratio and incorrect damage localization ratio. Among the 1st–3rd DFT-MCs, we find the lower order modes are more sensitive to noise. It is observed that the upper limitations of noise for the 1st–3rd DFT-MCs are respectively 0.03, 0.04 and 0.06 to obtain the acceptable results (larger than 70% accuracy ratio). As the increase of noise blurs the characteristic of crack, the ratio of false dismissal increases sharply while noise level rises. The incorrect damage localization (failure) ratio will ascend too, but with an approximate linear Table 2 Comparisons of different orders of DFT-MC Monte Carlo simulation results, 1000 samples, outlier identified by 3s, ( a, L c ) ¼ (2.4 mm, 321 mm). Modal order
Noise level ψ
Accuracy
False dismissal
Failure
1st
0.02 0.03 0.04
98.5% 76.9% 47.7%
0.8% 17.5% 42.1%
0.7% 5.6% 10.2%
2nd
0.02 0.04 0.06
99.7% 83.9% 47.9%
0.0% 10.0% 37.2%
0.3% 6.1% 14.9%
3rd
0.04 0.06 0.08
99.6% 87.8% 61.5%
0.2% 9.5% 29.0%
0.2% 2.7% 9.5%
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
(a)
(b)
775
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Fig. 15. The DFT-MC “Accuracy” and “False dismissal” samples for the crack case (a, Lc) ¼ (2.4 mm, 321 mm): (a–c) The 1st order DFT-MCs with ψ ¼0.02, 0.04, 0.06, (d–f) The 2nd order DFT-MCs with ψ ¼ 0.02, 0.03, 0.04, (g–i) The 3rd order DFT-MCs with ψ ¼ 0.04, 0.06, 0.08.
growth rate. As the supplement and numerical benchmark, the classical CD-MC method is also tested by the Monte Carlo simulation, and the related results are displayed in Table 3. One can see that the results given by the CD-MC method make no sense when noise exists. To obtain an acceptable accuracy rate, noise levels should be extremely small (1 10 5, 2 10 4 and 4 10 4 for 1st–3rd mode shapes, respectively). Further, we will discuss detection of smaller cracks. The results related to change of crack parameter to (a, Lc)¼ (0.6 mm, 321 mm) are shown in Table 4. The “Accuracy” and “False dismissal” samples are plotted in Fig. 16. In this case, the DFT-MC fails, the accuracy rate is 0.4–13.4%. The failure can be interpreted by the noise-free samples as shown in Fig. 16(d–f): Although crack is highlighted by spikes, the amplitudes of these spikes are not big enough to exceed the barriers given by Table 3 Comparisons of different orders of CD-MC Monte Carlo simulation results, 1000 samples, outlier identified by 3s, ( a, L c ) ¼(2.4 mm, 321 mm). Modal order
Noise level ψ
Accuracy
False dismissal
Failure
1st
1 10 5 0.02
75.3% 0.0%
24.3% 100.0%
0.4% 0.0%
2nd
2 10 4 0.02
72.7% 0.0%
27.0% 100.0%
0.3% 0.0%
3rd
4 10 4 0.02
79.3% 0.0%
20.7% 100.0%
0.0% 0.0%
776
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Table 4 Comparisons of different orders of DFT-MC Monte Carlo simulation results, 1000 samples, outlier identified by 3s, ( a, L c ) ¼ (0.6 mm, 321 mm). Modal order
Noise level ψ
Accuracy
False dismissal
Failure
1st 2nd 3rd 5th 7th 9th 10th
0.02 0.02 0.02 0.02 0.02 0.02 0.02
0.4% 13.3% 1.4% 82.2% 83.8% 0.6% 93.2%
68.9% 63.0% 93.5% 16.9% 16.2% 99.4% 6.8%
30.7% 23.7% 5.1% 0.9% 0.0% 0.0% 0.0%
3s. For the higher order DFT-MCs, one can observe the obvious improvements of damage identification results (as shown in Table 4 and Fig. 17). However, not all higher modes perform well which is indicated by the 9th order DFT-MC presented in Fig. 17(c). In application, one can hardly know which higher DFT-MC performs well, therefore, the identification obtained from different order DFT-MCs should be compared to avoid failure of damage identification. 3.4.2. Post-processing methods The characteristic of crack has been revealed in Fig. 16(d–f), but it is not strong enough to be classified as outlier. Hence, some post-processing methods could be added to intensify these characteristics. The curvature modes provided by the DFTMCs are employed as the inputs for the post-processing methods, which have been proved effective in combination with curvature modes. Then the new curves will be tested by statistical analysis to identify the outliers. In this paper, we will employ the continuous wavelet transform (CWT) as the post-processing method for DFT-MCs. Damage parameters are selected as ( a, L c ) ¼(0.6 mm, 321 mm), and the 5th, 7th, 9th and 10th orders mode shapes or curvature mode shapes are investigated due to the sensibility of higher mode shapes. Biorthogonal wavelet ‘bior6.8’ is selected in the CWT and the scale upper boundary is settled as 64. Identifications given by the CWT, the CD-MC and the present hybrid method (the DFT-MC and the CWT) for the same noise level ψ ¼0.02 are displayed in Fig. 18. It is observed in Fig. 18(a–d) that only boundary effects and some undulations induced by intact parts of mode shape are presented when the CWT is employed to decompose mode shape directly. In Fig. 18(e–h) one only can see some peaks relating to noise in low scale. However, the damage location can be easily found in Fig. 18(i–l) when the present hybrid method is used. The comparison validates the effectiveness of the hybrid method. For more stable and reasonable identification, statistical analysis is required. We propose to superpose the curves of all the scales to get a lumped curve which reflects all the information in the CWT, and enhances the characteristics of damage. Then outliers will be studied via the statistical analysis method mentioned above as shown in Fig. 19. It is seen that damage positions are clearly identified.
Fig. 16. The DFT-MCs “Accuracy” and “False dismissal” samples for the crack (a, Lc) ¼ (0.6 mm, 321 mm): (a) The 1st order DFT-MC with ψ ¼ 0.02, (b) The 2nd order DFT-MC with ψ ¼0.02, (c) The 3rd order DFT-MC with ψ ¼ 0.02, (d–f) The noise-free DFT-MC curves for the case shown in (a–c).
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
777
Fig. 17. The higher order DFT-MCs for the crack case (a, Lc) ¼ (0.6 mm, 321 mm) with ψ ¼ 0.02: (a) the 5th order, (b) the 7th order, (c) the 9th order, and (d) the 10th order.
L
L
L
L
L
L
L
L
L
L
Fig. 18. The CWT results for the crack case (a, Lc) ¼(0.6 mm, 321 mm) with ψ ¼ 0.02: (a–d) the results calculated via the 5th, 7th, 9th and 10th orders directly, (e–h) the results calculated via the 5th, 7th, 9th and 10th CD-MCs, (i–l) the results calculated via the 5th, 7th, 9th and 10th DFT-MCs.
To evaluate the accuracy, false dismissal and failure ratios of the hybrid method, the Monte Carlo simulation is employed again and the results are displayed in Table 5. By aid of the CWT, the present method gets a higher accuracy ratio and higher robustness. It should be noticed that not all the higher modes perform better than the lower ones (even with the aid of the CWT), for example the 9th order (shown in Table 5). The failure of the method is mainly caused by the boundary effect.
778
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 19. The superposition of all the scale CWT results and the outliers: (a–d) the results of 5th, 7th, 9th and 10th orders.
Table 5 Evaluation of the hybrid method, ( a, L c ) ¼ (0.6 mm, 321 mm). Case
Modal order
Noise level ψ
Accuracy
False dismissal
Failure
1 2 3 4 5 6 7 8
5th 5th 7th 7th 9th 9th 10th 10th
0.02 0.04 0.06 0.07 0.06 0.07 0.15 0.18
98.6% 70.6% 92.1% 75.4% 52.2% 45.9% 98.5% 77.4%
0.0% 0.0% 0.0% 0.0% 0.0% 0.2% 0.0% 0.0%
1.4% 29.4% 7.9% 24.6% 47.8% 53.9% 1.5% 22.6%
4. Experimental investigation The experimental investigation is based on the benchmark given by Katunin [36–38]. The benchmark can be obtained from http://ipkm.polsl.pl/wavstructdamas/, and the data we choose for the comparison are mode shapes measured for the case of beam with single crack and the baseline data (the intact structure). In this case, specimens (beam, a strip from a composite plate) are made by 24-layered glass-fiber reinforced plastic (GFRP) with unidirectional impregnated fibers and the [0°/60°/ 60°] lay-up. The dimensions of the specimens are as follows: length L ¼250 mm, width B ¼25 mm, and thickness H ¼5.28 mm. The crack is simulated via notch with the depth of 1 mm locating near 150 mm from the clamped end. The clamped length is about 20 mm. Based on the scanning LDV Polytec PSV-400, 44 sampling points are measured on each specimen to get the 1st–4th order mode shapes with a 215 mm effective length. Readers are suggested to refer to [36– 38] for more details. Firstly, the results obtained by the DFT-MC for intact (baseline) and damaged beams are presented in Fig. 20. Compared with the baseline, the identified results for the damaged beam show some singularities near 150 mm but not clear enough, especially in the 3rd and 4th order modes. Thus, the hybrid approaches are required. The damage assessments using the CD-MC, the CWT and the present hybrid method are all displayed in Figs. 21 and 22 for intact case and damaged case, respectively. The upper boundary of scale is settled as 6 considering the length of data. In assessments for intact case (Fig. 21), all three methods perform well, without any false alarm. It is observed that some small peaks and the undulations related to the waveform of mode shape in the results given by the CD-MC. The undulations related to the waveform can also be seen in the identifications generated by the DFT-MC (Fig. 21(i–l)), but not deemed as damage location. In the assessments obtained by the CWT directly (Fig. 21(e–h)), only some undulations related to the waveform and boundary effects exist. The upper and the lower bounds generated by statistical analysis in the present method provide a quantified criterion to classify the peaks and undulations, without which the normal characteristics in damage indexes may be deemed as the evidence of damage. In the damaged case, the CD-MC generates some obvious characteristics near damage location (Fig. 22(a–d)). However, these characteristics may be blurred by the irregular peaks due to numerical difference. In Fig. 22(e–h), the performances of the CWT (biorthogonal wavelet ‘bior6.8’) are better than the CD-MC, some ridges related to the damage can be easily found. However, the maximum values of these ridges appear near 160–170 mm instead of 150 mm. This localization error can be interpreted by the inappropriate choice of wavelet and parameters. It should be emphasized that the dependence of selection of wavelet type and parameters may restrict the application of the CWT in damage detection. In Fig. 22(i–l), the last assessments via present hybrid method is performed. Damage locations can be easily identified from these figures and extracted by aid of the indications of outliers.
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
779
Fig. 20. The 1st–4th order DFT-MCs of intact and cracked structures (a–d) for the intact structure, (e–h) for the cracked structure.
0.1
0.01
0.05
0.02
0.05
0
0
-0.01
-0.05
50 100 150 200
50 100 150 200
-0.1
50 100 150 200
L [mm]
L [mm] 6
5
5
4
4
4
4
50
100
150
2
1
200
3
50
100
150
1
200
50
L [mm]
L [mm] 4
100
150
1
200
-4
-10 50
100 150
L [mm]
200
w"
w"
-5
150
200
L [mm]
5
0
-2
100
10
5
0
50
L [mm]
10
2
3 2
5
w"
1
3 2
2
Scale
6
5
Scale
6
5
3
50 100 150 200
L [mm]
6
Scale
Scale
L [mm]
w"
0 -0.05
-0.02
-0.02
w"
w"
w"
w"
0
0 -5
50
100 150 200
L [mm]
0 -5
50
100 150
L [mm]
200
50
100 150 200
L [mm]
Fig. 21. Identifications using the intact structure benchmark data: (a–d) the 1st–4th order CD-MCs, (e–h) the 1st–4th order CWTs, (i–l) the results of the present hybrid method.
780
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 22. Identifications using cracked structure benchmark data: (a–d) the 1st–4th order CD-MCs, (e–h) the 1st–4th order CWTs, (i–l) the results of the present hybrid method.
5. Conclusions Modal curvature is an important dynamic feature in structural health monitoring. However, the susceptibility to noise induced by numerical derivative severely impairs the applicability in damage detection. To tackle this problem, the DFT-MC is proposed and some auxiliary methods are introduced in this paper. By aid of the characteristics of the FT, the numerical derivative in the CD-MC is accomplished by multiplication and the DFT synthesis, which enhances the numerical stability of modal curvature estimation. In addition, the wavenumber domain in the DFT provides an effective method for space filtering. Based on the iteration, statistical analysis and the post-processing method, a robust hybrid algorithm is generated. It should be mentioned that the present method is just a modification for estimating MC. Restricted by the fundamental properties of MC, the capacity of present method for small damage detection may be weak in applications. The singularity related to damage presence can be easily covered by noise. In that case, it appears few methods can conduct the satisfactory noise suppression and damage identification meanwhile based on limited time-domain measurement of mode shapes. Acknowledgments Special thanks should be given to Prof. A. Katunin, who shares the benchmark data in the website. This work is supported by the National Natural Science Foundation of China (No. 51405369), the China Scholarship Council, and the China Postdoctoral Science Foundation (No. 2014M560766).
References [1] S.P. Lele, S.K. Maiti, Modelling of transverse vibration of short beams for crack detection and measurement of crack extension, J. Sound Vib. 257 (2002)
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
781
559–583. [2] S.M. Murigendrappa, S.K. Maiti, H.R. Srirangarajan, Frequency-based experimental and theoretical identification of multiple cracks in straight pipes filled with fluid, Ndt E Int. 37 (2004) 431–438. [3] M.R. Naniwadekar, S.S. Naik, S.K. Maiti, On prediction of crack in different orientations in pipe using frequency based approach, Mech. Syst. Signal Process. 22 (2008) 693–708. [4] Z. Su, L. Ye, Y. Lu, Guided Lamb waves for identification of damage in composite structures: a review, J. Sound Vib. 295 (2006) 753–780. [5] Z. Su, L. Cheng, X. Wang, L. Yu, C. Zhou, Predicting delamination of composite laminates using an imaging approach, Smart Mater. Struct. 18 (2009) 074002. [6] M. Cao, L. Ye, L. Zhou, Z. Su, R. Bai, Sensitivity of fundamental mode shape and static deflection for damage identification in cantilever beams, Mech. Syst. Signal Process. 25 (2011) 630–643. [7] X. Qian, M. Cao, Z. Su, J. Chen, A hybrid particle swarm optimization (PSO)-simplex algorithm for damage identification of delaminated beams, Math. Probl. Eng. 2012 (2012) 607418. [8] Z.-B. Yang, X.-F. Chen, Y. Xie, H. Zuo, H.-H. Miao, X.-W. Zhang, Wave motion analysis and modeling of membrane structures using the wavelet finite element method, Appl. Math. Model. 40 (2016) 2407–2420. [9] Z. Yang, X. Chen, X. Li, Y. Jiang, H. Miao, Z. He, Wave motion analysis in arch structures via wavelet finite element method, J. Sound Vib. 333 (2014) 446–469. [10] A.K. Pandey, M. Biswas, M.M. Samman, Damage detection from changes in curvature mode shapes, J. Sound Vib. 145 (1991) 321–332. [11] Z.B. Yang, X.F. Chen, S.H. Tian, Z.J. He, Multiple damages detection in beam based approximate waveform capacity dimension, Struct. Eng. Mech. 41 (2012) 663–673. [12] M. Cao, W. Xu, W. Ostachowicz, Z. Su, Damage identification for beams in noisy conditions based on Teager energy operator-wavelet transform modal curvature, J. Sound Vib. 333 (2014) 1543–1553. [13] M. Cao, M. Radzieński, W. Xu, W. Ostachowicz, Identification of multiple damage in beams based on robust curvature mode shapes, Mech. Syst. Signal Process. 46 (2014) 468–480. [14] Z.-B. Yang, M. Radzienski, P. Kudela, W. Ostachowicz, Two-dimensional modal curvature estimation via Fourier spectral method for damage detection, Compos. Struct. 148 (2016) 155–167. [15] W. Fan, P.Z. Qiao, A strain energy-based damage severity correction factor method for damage identification in plate-type structures, Mech. Syst. Signal Process. 28 (2012) 660–678. [16] M.S. Cao, P.Z. Qiao, Novel Laplacian scheme and multiresolution modal curvatures for structural damage identification, Mech. Syst. Signal Process. 23 (2009) 1223–1242. [17] J. Xiang, T. Matsumoto, Y. Wang, Z. Jiang, Detect damages in conical shells using curvature mode shape and wavelet finite element method, Int. J. Mech. Sci. 66 (2013) 83–93. [18] R.B. Bai, M.S. Cao, Z.Q. Su, W. Ostachowicz, H. Xu, Fractal dimension analysis of higher-order mode shapes for damage identification of beam structures, Math. Probl. Eng. 2012 (2012) 454568. [19] M. Radzienski, M. Krawczuk, M. Palacz, Improvement of damage detection methods based on experimental modal parameters, Mech. Syst. Signal Process. 25 (2011) 2169–2190. [20] M. Radzienski, M. Krawczuk, W. Ostachowicz, Experimental verification and comparison of mode shape-based damage detection methods, Damage Assess. Struct. VIII 413-414 (2009) 699–706. [21] A. Katunin, P. Przystałka, Automated wavelet-based damage identification in sandwich structures using modal curvatures, J. Vibroeng. 17 (2015). [22] Z. Yang, X. Chen, J. Yu, R. Liu, Z. Liu, Z. He, A damage identification approach for plate structures based on frequency measurements, Nondestruct. Test. Eval. (2013) 1–21. [23] X. Zhang, R.X. Gao, R. Yan, X. Chen, C. Sun, Z. Yang, Multivariable wavelet finite element-based vibration model for quantitative crack identification by using particle swarm optimization, J Sound Vib. 375 (2016) 200–216. [24] J. Xiang, M. Liang, Y. He, Experimental investigation of frequency-based multi-damage detection for beams using support vector regression, Eng. Fract. Mech. 131 (2014) 257–268. [25] Z. Yang, X. Chen, Y. Jiang, Z. He, Generalised local entropy analysis for crack detection in beam-like structures, Nondestruct. Test. Eval. 29 (2014) 133–153. [26] M. Hong, Z. Su, Y. Lu, H. Sohn, X. Qing, Locating fatigue damage using temporal signal features of nonlinear Lamb waves, Mech. Syst. Signal Process. 60 (2015) 182–197. [27] Z. Su, C. Zhou, M. Hong, L. Cheng, Q. Wang, X. Qing, Acousto-ultrasonics-based fatigue damage characterization: linear versus nonlinear signal features, Mech. Syst. Signal Process. 45 (2014) 225–239. [28] B.P. Nandwana, S.K. Maiti, Detection of the location and size of a crack in stepped cantilever beams based on measurements of natural frequencies, J. Sound. Vib. 203 (1997) 435–446. [29] Z.B. Yang, X.F. Chen, Y. Xie, X.W. Zhang, The hybrid multivariate analysis method for damage detection, Struct. Control Health Monit. 23 (2016) 123–143. [30] Z.B. Yang, X.F. Chen, Y. Xie, H.H. Miao, J.J. Gao, K.Z. Qi, Hybrid two-step method of damage detection for plate-like structures, Struct. Control Health Monit. 23 (2016) 267–285. [31] A. Deraemaeker, A. Preumont, Vibration based damage detection using large array sensors and spatial filters, Mech. Syst. Signal Process. 20 (2006) 1615–1630. [32] G. Tondreau, A. Deraemaeker, Numerical and experimental analysis of uncertainty on modal parameters estimated with the stochastic subspace method, J. Sound Vib. 333 (2014) 4376–4401. [33] G. Tondreau, A. Deraemaeker, Automated data-based damage localization under ambient vibration using local modal filters and dynamic strain measurements: Experimental applications, J. Sound Vib. 333 (2014) 7364–7385. [34] G. Tondreau, A. Deraemaeker, Local modal filters for automated data-based damage localization using ambient vibrations, Mech. Syst. Signal Process. 39 (2013) 162–180. [35] G. Tondreau, A. Deraemaeker, Experimental Localization of Small Damages using Modal Filters, Special Topics in Structural Dynamics, Springer, Berlin, Germany (2013), pp. , 2013, 585–591. [36] A. Katunin, Nondestructive damage assessment of composite structures based on wavelet analysis of modal curvatures: state-of-the-art review and description of wavelet-based damage assessment benchmark, Shock Vib. 501 (2015) 735219. [37] A. Katunin, K. Dragan, M. Dziendzikowski, Damage identification in aircraft composite structures: a case study using various non-destructive testing techniques, Compos. Struct. 127 (2015) 1–9. [38] A. Katunin, Diagnostics of Composite Structures using Wavelets, Radom, Gliwice, Poland, 2015.
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp
Fourier spectral-based modal curvature analysis and its application to damage detection in beams Zhi-Bo Yang a,b,n, Maciej Radzienski b, Pawel Kudela b, Wieslaw Ostachowicz b,c a b c
School of Mechanical Engineering, Xi’an Jiaotong University, 710049 Xi’an, PR China Institute of Fluid Flow Machinery, Polish Academy of Science, 80-231 Gdansk, Poland Faculty of Automotive and Construction Machinery, Warsaw University of Technology, 02-524 Warsaw, Poland
a r t i c l e i n f o
abstract
Article history: Received 15 February 2016 Received in revised form 10 May 2016 Accepted 2 July 2016 Available online 18 July 2016
In this paper, a simple Fourier spectral-based method is proposed to calculate the modal curvature (MC) of beams instead of the traditional central difference method. Based on the present method, damages in beam-like structures are localized. The present method provides an alternative selection to estimate MC in damage detection. There are two advantages of the present method. Firstly, the spectral calculation of spatial derivatives is conducted globally, which provides the suppression for noise. In addition, signal processing in the wavenumber domain provides an alternative choice for spatial filtering for mode shapes. Secondly, the proposed method provides a precise estimation of the MC which is related to original definition. With the absence of numerical derivative, the estimated results can be more stable and robust. Statistical analysis is conducted to show the effectiveness and noise immunity of the proposed method. In order to obtain the better identification, the MC calculated by the proposed method is employed as the input of continuous wavelet transform, and then the hybrid method is generated. The validations of the present method and comparison with the traditional central difference method are numerically and experimentally demonstrated. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Fourier transform Modal curvature Damage detection Beam Wavenumber filtering
1. Introduction Structural damage detection has been a research focus with increasing interest in mechanical, civil, aerospace and some other fields during the last few decades [1–9]. The modal curvature (MC) has been investigated in a number of studies [10– 14]. The MC is the curvature of mode shapes of structures. In beam-like structures, it can be expressed as:
w″( x) =
d2w ( x) M = − dx2 EI ( x)
(1)
where w′′(x) is the MC, w(x) is the displacement mode shape of a beam, M is the bending moment and EI(x) is the bending stiffness. A reduction of the bending stiffness induced by damage causes alteration in the MC. This change is directly related to damage location. Therefore, the MC is a theoretically feasible damage index in beam-like structures. However, only few sensor types can provide satisfactory direct measurement of the MC. Instead, some local estimation methods, such as the n
Corresponding author at: School of Mechanical Engineering, Xi’an Jiaotong University, 710049 Xi’an, PR China. E-mail addresses: [email protected], [email protected], [email protected] (Z.-B. Yang).
http://dx.doi.org/10.1016/j.ymssp.2016.07.005 0888-3270/& 2016 Elsevier Ltd. All rights reserved.
764
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Nomenclature
ε λ ψ Δ a ai cr e h i f ( ∙) k t wd w ( x) ^ ( k) w wx
the convergence error the random parameter the signal to noise ratio parameter the height of rectangular impulse the depth of crack the constant about Nuttall window the constant about damage the constant about convergence error the sampling interval imaginary number the dimensionless local compliance function wavenumber time the damaged mode shape the displacement mode shape of beam the Fourier transform of w ( x ) the discrete form of w (x )
^ w k w′( x ) w″( x ) wx′ wx″ yn(t) yn∘ ( t ) B E EI ( x ) H Kc L Lc Ll M N WNuttall ‖∙‖∞
^ (x ) the discrete form of w the first-order derivative of w (x ) the modal curvature the discrete form of w′( x ) the discrete form of w″( x ) the noisy response of measuring points n the noise-free response of measuring points n the width of beam the Young's modulus the bending stiffness the thickness of beam the stiffness of crack the length of beam crack location the width of rectangular impulse. the bending moment the number of grid point the Nuttall window the infinite norm
numerical difference techniques, are usually used in the MC approximations. Among these numerical difference techniques, the central difference approximation is usually employed in the MC estimation. For simplicity, the MC calculated by the central difference approximation is denoted as the CD-MC. The formulation of the CD-MC is:
w″( x) =
w ( x − h) − 2w ( x) + w ( x + h) h2
(2)
where h is the sampling interval. The application of the MC in damage detection can be traced back to Pandey's work in 1991 [10], in which the abrupt change of the MC is seen as the indication of damage. Following this work, some investigations have been conducted, and the MC is also used as a baseline-free damage index as presented in [12,15–17]. Although the effectiveness of Eq. (2) has been verified in these works, we are willing to point some possible deficiencies and validate them in this paper: (1) Ill-posed condition. False peaks in the CD-MC can be easily affected by noise. The problem of numerical difference (containing the central difference) is known to be ill-posed in the sense that a small measurement error/noise induces a large error in the computed derivatives. (2) Approximation order. Mode shape can be represented by trigonometric function. However, local estimation of mode shape is often done by a polynomial approximation. Thus, the approximation error is relatively large, which influences the accuracy of numerical difference further. (3) Noise suppression. Due to the local estimation, the capacity of the CD-MC for noise suppression is relatively weak. To solve these problems, some investigations have been done [13,16,18–21]. However, these studies approximate the MC by using the central difference. Therefore, the aim of this research work is to provide the MC estimation without the central difference approximation. Paying attention to problems (1–3), one can find that the Fourier spectral method is the natural choice for solving them by aid of the following properties: (1) Special process in derivative calculation. The difference of functions can be easily conducted via multiplication in wavenumber domain of the Fourier transform (FT), and reconstruction via the inverse Fourier transform (IFT). (2) The approximation function in Fourier spectral method is trigonometric function. (3) Noise suppression via averaging in the mapped domain. Fourier spectral-based methods are global, thus the parts of the noise can be suppressed in the wavenumber domain via special filters. Compared with the frequency-based method [22–24], modal shape-based method (containing MC) performs better in damage localization. But the modal shape-based method [25] is also argued by its sensitivity to small damage. Modal test often shows the low-frequency domain information compared with the high-frequency domain information provided by wave propagation-based method [26,27]. The identification of them can be complementary to each other, especially for large structures. The aim of this paper is to provide an alternative and simple method to estimate the MC used in damage detection. The following parts of this paper are organized as follows: Section 2 presents the formulas and implementation of the Fourier
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
765
spectral method for the MC estimation. Section 3 mainly focuses on the numerical comparison between the proposed method and the CD-MC method. In addition, some auxiliary methods are also introduced to improve the performance of the proposed method. Section 4 validates the proposed method via benchmark test. Finally, some discussions and conclusions are presented in Section 5.
2. Fourier spectral method for MC estimation 2.1. Theoretical basis Fourier spectral method is one of the most important techniques for numerical solution of partial differential equations. For the sake of completeness theoretical basis are explained in this section. The FT is used in signal analysis to map signals from the time domain to the frequency domain. Let's consider the mode shape w(x) as the function to be transformed by the FT and the IFT in the spatial and wavenumber domain: +∞
^ ( k) = w
∫−∞
w ( x) =
1 2π
e−ikxw ( x) dx +∞
∫−∞
(3)
^ ( k ) dk eikxw
(4)
where x ∈ and k ∈ . k is called wavenumber with the specific physical meaning. Paying attention to the distinguished characteristic of the FT:
dw ( x) 1 i = dx 2π
w′( x) =
+∞
∫−∞
^ ( k ) dk eikxkw
(5)
one can get the analytical form of the MC instead of Eq. (2):
w″( x) =
dw′( x) 1 = − dx 2π
+∞
∫−∞
^ ( k ) dk eikxk2w
(6)
However, Eq. (6) cannot be employed to estimate the MC directly because the analytical form of the mode shape w(x) is unknown in application. What is worse, the integration in Eq. (6) is defined on the interval [ −∞ , + ∞]. Given that the mode shape w(x) is discrete and bounded in numerical estimation and experiment, the discrete Fourier transform (DFT) is the natural choice for the MC calculation. Without the loss of generality, it can be assumed that w(x) is defined on the interval [0, 2π], and the interval is further divided by N 1 equal subintervals between N grid points. Let's define the length of each subinterval as h, then the related DFT synthesis can be written as:
^ = h∑ e−ikxw w k x x
wx =
1 ∑ eikxw^ k 2π k
(7)
(8)
^ are the discrete forms of w ( x ) and w ^ ( x ) respectively. Corresponding to (Eqs. (7) and 8), (Eqs. (5) and 6) are where wx and w k rewritten in the DFT frame as:
wx′ =
⎛ ⎞ h i∑ keikx ⎜⎜ ∑ e−ikxwx ⎟⎟ 2π k ⎝ x ⎠
wx″ = −
⎛ ⎞ h ∑ k2eikx ⎜⎜ ∑ e−ikxwx ⎟⎟ 2π k ⎝ x ⎠
(9)
(10)
In (Eqs. (9) and 10), the x and k are all bounded and discrete, without any technical restrictions. Based on Eq. (10), the MC can be calculated easily. This process can be seen as a trigonometric interpolation, thus Eq. (10) provides a semi-analytical method to estimate the MC. According to the characteristic of trigonometric functions, the problems (1, 2) are solved. Problem (3) will be solved by filtering in the wavenumber domain. 2.2. Damage models Damage in structures can be manifested by the singularity in mode shapes. Thus, it is necessary to investigate and classify singularities, based on which generally applicable wavenumber domain filters are provided. Before the discussion about wavenumber domain filter, the procedure for establishing damage models of beam is presented in this part of the
766
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 1. The influences of linear spring crack on the 1st–3rd mode shapes: crack features in spatial domain (a–c) and in wavenumber domain (d–f).
filter selection. In this paper, the numerical simulations are conducted using the cantilever Euler–Bernoulli beam model with the rectangular cross-section. The structure used in numerical test is a beam with the length L ¼535 mm, the width B ¼20 mm and the thickness H¼12 mm. The selection of material can be arbitrary, thus it is not highlighted here. Given that transverse vibration of the Euler–Bernoulli beam concerns only the bending vibrations of beam, the rotational spring constant is assumed to be dominant in the local flexibility matrix [1,3], and the stiffness of rotational spring is given by [9,28–30]:
Kc =
BH 2E 2
72π ( a/H ) f ( a/H )
(11)
where E is the Young's modulus, B is the width of beam, H is the depth of beam, and f(a/H) is the dimensionless local compliance function calculated as:
Fig. 2. The influences of rectangular model on the mode shape: (a) the damaged 1st mode shape, (b) damage feature in spatial domain, (c) damage feature in wavenumber domain.
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
⎛ a⎞ ⎛ a ⎞2 ⎛ a ⎞3 ⎛ a ⎞4 ⎛ a ⎞5 ⎛ a ⎞6 a f ⎜ ⎟ = 0.6384 − 1.035 + 3.7201 ⎜ ⎟ − 5.1773 ⎜ ⎟ + 7.553 ⎜ ⎟ − 7.332 ⎜ ⎟ + 2.4909 ⎜ ⎟ ⎝ H⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ H⎠ H H H H H
767
(12)
Based on (Eqs. (11) and 12), the influence of a crack on mode shape can be calculated easily by finite element method or analytical model. To illustrate effects caused by the crack model, a severe crack (a¼4.8 mm) is arranged at the location Lc ¼321 mm. The differential waveforms between cracked and intact mode shapes and the corresponding wavenumber domain representations are displayed in Fig. 1. Although crack features of 1st–3rd mode shapes are totally different in the spatial domain, the wavenumber domain representations of them are similar (belonging to the low wavenumber domain). Namely, one can conclude that the crack related features are concentrated in the low wavenumber domain via DFT. Second model, the rectangular model, is simpler compared with the linear spring model. By adding a constant component cr on the specific interval of the intact mode shape w, the damaged mode shape wd can be obtained. The formulation of this model is:
⎧ ⎪ wx + cr wd = ⎨ ⎪ ⎩ wx
( x = Lc ) ( for the others)
(13)
Fig. 2 shows the damage features with the parameter cr ¼0.02. As the damage features in every order mode shapes are the same, thus only the first mode shape is presented. Comparing the spatial domain representation with the wavenumber domain representation, one can find the feature of this damage model mainly concentrates in low wavenumber domain, similar to the linear spring model. The third model (zero-frequency component model) is formulated as:
⎧ x ∈ ( L c , L c + h, ... , L ) ⎪ wx + cr wd = ⎨ ⎪ x ∈ ( 0, L/N , ... , L c − h) ⎩ wx
(14)
In this model, the effect induced by damage is considered as a leap at the damage location Lc, and the area from damage location Lc to the free end L is also affected as show in Fig. 3. Similarly, damage features concentrates in low wavenumber domain. 2.3. Wavenumber domain filter Fig. 4 presents the representations of noise in spatial domain and wavenumber domain. The waveform shown in Fig. 4 (a) consists of the pseudo-random value generated by ‘rand’ in Matlab to simulate a white noise. We can see the noise energy distributes evenly in spatial domain and wavenumber domain. It is different with the features of the singularities mentioned above. Thus, a low-pass filter can be used in the wavenumber domain to suppress noise. In this paper, a Nuttall window WNuttall is applied in the wavenumber domain directly, which is given by:
WNuttall = a 0 − a1 cos
π ( 2hk + L ) L
+ a2 cos
2π ( 2hk + L ) L
− a3 cos
3π ( 2hk + L ) L
−
L L ≤k≤ 2h 2h
(15)
where a0 ¼ 0.3635819, a1 ¼0.4891775, a2 ¼0.1365995 and a3 ¼ 0.0106411. The comparison of the Nuttall window, the Blackman window and the Hann window is displayed in Fig. 5, in which one can observe that the Nuttall window attenuates faster than the other two windows. This is an appropriate property for noise suppression in high wavenumber domain referring to Figs. (1–4). The procedure of wavenumber domain filtering is
Fig. 3. The influences of zero-frequency component model on the mode shape: (a) the damaged 1st mode shape, (b) damage feature in spatial domain, (c) damage feature in wavenumber domain.
768
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 4. The features of noise: (a) the contaminated 1st mode shape, (b) noise in the spatial domain, (c) noise in the wavenumber domain.
Fig. 5. Comparison of window functions.
formulated as:
^ i+1 = W ^i w Nuttall *wx x
(16)
The superscript depicts the step in iteration. Given that the required times of the wavenumber domain filtering is un-
^1 known, the iteration shown in Fig. 6 is required. In Fig. 6, ε depicts the convergence error. We suggest to choose ε = e w x
‖ ‖, ∞
where e is set as 1% in this paper (the choice is flexible). To present the effectiveness of the iteration filtering, the intact mode shape, which is contaminated by the same noise level (shown in Fig. 4), before and after 5 times filtering processes is presented via spatial and wavenumber representations in Fig. 7. Noise is considerably reduced by the proposed procedure in the wavenumber domain, and the reconstructed mode shape is smoother. Further investigations and validations about the iteration will be presented in the following section.
3. Method verification 3.1. Damage scenarios The performance of the proposed method (denoted as the DFT-MC) for damage identification is investigated for a wide spectrum of damage scenarios (Table 1), and the CD-MC is employed as the numerical benchmark of the DFT-MC. In a progressive manner in terms of damage identification for different scenarios and the noise immunity tests, comprehensive comparisons are conducted between the DFT-MC and the CD-MC. In Table 1, the spectrum of damage scenarios are provided following damage models mentioned in the above section. For the linear spring model, the damage parameter vector is defined as ( a, L c ). For the rectangular model, the damage parameter vector is defined as ( Δ, Ll, L c ), where Δ is the height of rectangular impulse and Ll is the width of rectangular impulse.
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
769
w
Fig. 6. The flow chart of the wavenumber domain filtering.
L
k
Fig. 7. Illustration for the iteration result: (a) the contaminated 1st mode shape before and after filtering, (b) the noise before and after filtering.
The damage parameter vector for the zero-frequency component model is defined as ( Δ, L c ), where Δ is the height of step change. Damage scenarios are organized as follow: (1) Scenarios 1–3 are used to verify the identification capacity of DFT-MC for the three models, respectively. (2) Scenarios 4 and 5 are used to verify the capability of the present method for mixed damage identification with different damage severities. Scenarios 5 and 6 are used to verify the capability of the present method for mixed damage identification with different damage locations. (3) In scenarios 1–6, the 1st–3rd mode shapes are used in validations. The 1st, 5th and 10th mode shapes are used in scenario 7 for the purpose of comparison in low, middle and high mode shapes, respectively. (4) The further comparisons between CD-MC and DFT-MC are based on scenario 6 via noise immunity tests. (5) Scenario 8 is used for further comparisons shown in Section 3.4. It should be mentioned that cyclic continuations have been used in the following examples to avoid boundary effects, either for CD-MC or DFT-MC.
770
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Table 1 Damage scenarios used in validations. Damage scenario
Mode shape order
Linear spring model
( Δ,
( a,
[mm, mm, mm]
[mm, mm]
– (0.02, 15, 321) – (0.02, 15, 321) (0.01, 15, 321) (0.01, 15, 421) (0.01, 15, 421) –
– – (0.02, 321) (0.02, 421) (0.01, 421) (0.01, 221) (0.01, 221) –
Lc )
Rectangular model
Ll , L c )
Zero-frequency component model ( Δ, L c )
[mm, mm] 1 2 3 4 5 6 7 8
1–3 1–3 1–3 1–3 1–3 1–3 1, 5, 10 1–3
(2.4, 321) – – (2.4, 221) (1.2, 221) (1.2, 321) (1.2, 321) (0.6, 321)
Fig. 8. Scenarios 1–3: (a–c) the 1st–3rd DFT-MCs of scenarios 1, (d–f) the 1st–3rd DFT-MCs of scenarios 2, (g–i) the 1st–3rd DFT-MCs of scenarios 3. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
771
Fig. 9. Scenarios 4 and 5: (a–c) the 1st–3rd DFT-MCs of scenarios 4, (d–f) the 1st–3rd DFT-MCs of scenarios 5.
3.2. Noise-free test Fig. 8 presents the identifications related to scenarios 1–3 via DFT-MC. The present method identifies all the three kinds of damage clearly. The damage localization also agrees well with the actual locations (depicted by red triangles) The damage features displayed by modal curvature are different for these three damage models. This provides an evidence for distinguishing different type of damage. This group of simulations verifies the capacity of DFT-MC for the identification of different damages. Scenarios 4 and 5 combine three models in one mode shape, and the results are displayed in Fig. 9. The results agree with the actual conditions. It is also observed that singularity induced by the linear spring model increases with the order of mode shape. Either in scenarios 4 or 5, three singularities are all identified clearly. In higher modal curvatures, the related amplitudes are also distinct. This group of simulations verifies that DFT-MC is capable of detecting mixed damage. In scenario 6, the locations of damage are changed compared with scenario 5. However, one can see that the identification accuracy is not affected by this change (see in Fig. 10). Thus, the adaptability of the present method is validated to some extent. Fig. 11 presents a preliminary comparison between DFT-MC and CD-MC via scenario 7. This scenario is designed to verify and compare the identification capabilities of DFT-MC and CD-MC via higher mode shapes and damage mixture. These two
Fig. 10. Scenario 6: (a–c) the 1st–3rd DFT-MCs.
772
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 11. Scenario 7: (a–c) the 1st–3rd DFT-MCs, (d–f) the 1th–3rd CD-MCs.
methods obtain accurate identification for three damage models. The detected singularity via CD-MC in Fig. 11(d) at 0.321 m induced by the linear spring model is weaker than the one detected by DFT-MC in Fig. 11(a) at the same location. However, the conclusion that DFT-MC is better than CD-MC can be hardly obtained just from the above performance. Thus, further investigations are displayed via noise immunity test in the following part. 3.3. Noise immunity test Before the comparison in noise immunity test, a fact should be noticed: In the following comparisons, filtering is just applied in DFT-MC in wavenumber domain. Here, the filtering is emphasized in spatial domain, the filtering in the time domain during measurement is not included. That is to say: the results given by CD-MC is not filtered, thus comparison seems to be unfair for CD-MC in this view. However, the implementation of filtering is easier to perform in the wavenumber domain than in the spatial domain. The wavenumber domain filtering is one of the advantages of the present method. Of course, the spatial or modal filtering methods are available in the literatures [31–35], which provide the possibilities for the filtering of CD-MC. This problem is not investigated further in this paper. In this paper, noise is defined as the white Gaussian noise (WGN). The simulated WGN is added on the waveform directly
Fig. 12. (a–c) the 1st–3rd CD-MCs based on noisy mode shapes, ψ ¼0.02.
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
773
Fig. 13. DFT-MCs based on noisy mode shapes, (a) the 1st DFT-MC with ψ ¼0.02, (b) the 2nd DFT-MC with ψ ¼0.02, (c) the 3rd DFT-MC with ψ ¼0.04, (d) the 1st DFT-MC with ψ ¼0.03, (e) the 2nd DFT-MC with ψ ¼ 0.04, (f) the 3rd DFT-MC with ψ ¼0.05.
referring to [32]:
yn ( t ) = yn0 ( t ) + ψλ max ⎡⎣ yn0 ( t ) ⎤⎦
(15)
where yn(t) and ( t ) are the noisy and noise-free responses of measuring points n, λ is the random parameter, with its continuous distribution following a Gaussian distribution with zero mean and unitary standard deviation, ψ is the parameter which controls the signal to noise ratio. Fig. 12 shows the CD-MC results for scenario 6 (ψ ¼0.02). The noise amplified by the central difference severely impairs the profile of the modal curvature. From these figures, useful damage features can be hardly extracted. This example illustrates the susceptibility of CD-MC to noise. The identifications given by DFT-MC are presented in Fig. 13. The noise levels are assigned as ψ ¼ 0.02 and 0.03 for the 1st order DFT-MC, ψ ¼0.02 and 0.04 for the 2nd order DFT-MC, and ψ ¼0.04 and 0.05 for the 3rd order DFT-MC. According to the flow chart shown in Fig. 6, most of the locations of the three defects are represented clearly in each noise level. The identified damage locations agree well with the actual damage locations. It is observed in Fig. 13 that the higher order of mode shape causes lower noise susceptibility of the DFT-MCs. In addition, some false peaks (the high-wavenumber disturbances) also appear in Fig. 13, which disturb the identification of damage locations. Thus, the misjudgment and location error may emerge. In the following part, we will discuss the influences and the inhibition methods for the low/high-wavenumber disturbances by aid of statistical analysis and Monte Carlo simulation.
yn∘
3.4. Statistical analysis and post-processing method Disturbances are divided into two groups, a low- and a high-wavenumber group. The low-wavenumber disturbances are related to the waveforms of higher mode shapes, whose fluctuation will be expressed in DFT-MC too. Some post-processing damage detection methods are suggested to restrain this part of disturbance. 3.4.1. Statistical analysis Appearance of a noise may cause an error in identification and a probability of misjudgment. Influence of the noise can be hardly eliminated, but statistical analysis can help to minimize undesirable conditions. There is a requirement for the damage index, which helps to indicate the peak related to damage location among the results given by curvature methods. As most points in the measured system are related to the intact parts of structures, the
774
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 14. The illustrations of typical samples: (a) correct identification, (b) identification with locating error, (c) false dismissal, (d) failure. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
damage points become the outliers compared with the intact points distributed in a certain field. Selection of quantile to determine outlier is important: If the quantile is chosen too small, risk of erroneous judgment will sharply increase, and vice versa. An appropriate choice of quantile can be determined by the balance of risk and false dismissal probability. In this paper, quantile is selected as 3s [30]. In this part, the noise-polluted scenarios 1 and 8 are employed as the basic cases to evaluate the error-ratio/accuracy of the DFT-MC method according to Monte Carlo simulation. Some typical samples are displayed in Fig. 14, in which the filled triangles on the bottom of these figures indicate the actual location of damage, and the unfilled triangles plotted on the curves indicate the identified location (if any). The red lines parallel to the horizontal axes are the upper and lower bounds determined by 3s. The following terms should be highlighted before the simulation: (1) It can be observed in Fig. 14(a, b, and d) that the peak beyond the bounds will be deemed as damage location. If there is no point beyond the bounds, the structure will be deemed intact as show in Fig. 14(c). (2) It is noticed in Fig. 14(a) that the points beyond the bounds determine a damaged interval instead of a specific location. In this paper, we employ the peak location in this interval as the damage location instead of the center of the interval. (3) The identification as shown in Fig. 14(b) is defined as damage localization error because the difference between identified and actual location is smaller than 1% of the total length of the beam. Other identification will be seen as incorrect damage localization which means that the damage identification procedure failed. This case is indicated within the paper as failure. (4) Other complex conditions can be deemed as the combinations of the typical cases shown in Fig. 14. The Monte Carlo simulation consists of some sub-simulations: scenarios 1 with different noise levels and scenarios 8 with different noise levels. Each sub-simulation contains 1000 samples. According to the Chebyshev law of large number, the capacity as 1000 in simulation is enough to fulfill the requirement that sample distribution should cover the entire solution field. By choosing the different noise levels, some typical evaluations are listed in Table 2 and shown in Fig. 15. It should be mentioned that different noise levels are employed for the 1st–3rd order DFT-MCs to obtain similar accuracy ratio, false dismissal ratio and incorrect damage localization ratio. Among the 1st–3rd DFT-MCs, we find the lower order modes are more sensitive to noise. It is observed that the upper limitations of noise for the 1st–3rd DFT-MCs are respectively 0.03, 0.04 and 0.06 to obtain the acceptable results (larger than 70% accuracy ratio). As the increase of noise blurs the characteristic of crack, the ratio of false dismissal increases sharply while noise level rises. The incorrect damage localization (failure) ratio will ascend too, but with an approximate linear Table 2 Comparisons of different orders of DFT-MC Monte Carlo simulation results, 1000 samples, outlier identified by 3s, ( a, L c ) ¼ (2.4 mm, 321 mm). Modal order
Noise level ψ
Accuracy
False dismissal
Failure
1st
0.02 0.03 0.04
98.5% 76.9% 47.7%
0.8% 17.5% 42.1%
0.7% 5.6% 10.2%
2nd
0.02 0.04 0.06
99.7% 83.9% 47.9%
0.0% 10.0% 37.2%
0.3% 6.1% 14.9%
3rd
0.04 0.06 0.08
99.6% 87.8% 61.5%
0.2% 9.5% 29.0%
0.2% 2.7% 9.5%
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
(a)
(b)
775
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Fig. 15. The DFT-MC “Accuracy” and “False dismissal” samples for the crack case (a, Lc) ¼ (2.4 mm, 321 mm): (a–c) The 1st order DFT-MCs with ψ ¼0.02, 0.04, 0.06, (d–f) The 2nd order DFT-MCs with ψ ¼ 0.02, 0.03, 0.04, (g–i) The 3rd order DFT-MCs with ψ ¼ 0.04, 0.06, 0.08.
growth rate. As the supplement and numerical benchmark, the classical CD-MC method is also tested by the Monte Carlo simulation, and the related results are displayed in Table 3. One can see that the results given by the CD-MC method make no sense when noise exists. To obtain an acceptable accuracy rate, noise levels should be extremely small (1 10 5, 2 10 4 and 4 10 4 for 1st–3rd mode shapes, respectively). Further, we will discuss detection of smaller cracks. The results related to change of crack parameter to (a, Lc)¼ (0.6 mm, 321 mm) are shown in Table 4. The “Accuracy” and “False dismissal” samples are plotted in Fig. 16. In this case, the DFT-MC fails, the accuracy rate is 0.4–13.4%. The failure can be interpreted by the noise-free samples as shown in Fig. 16(d–f): Although crack is highlighted by spikes, the amplitudes of these spikes are not big enough to exceed the barriers given by Table 3 Comparisons of different orders of CD-MC Monte Carlo simulation results, 1000 samples, outlier identified by 3s, ( a, L c ) ¼(2.4 mm, 321 mm). Modal order
Noise level ψ
Accuracy
False dismissal
Failure
1st
1 10 5 0.02
75.3% 0.0%
24.3% 100.0%
0.4% 0.0%
2nd
2 10 4 0.02
72.7% 0.0%
27.0% 100.0%
0.3% 0.0%
3rd
4 10 4 0.02
79.3% 0.0%
20.7% 100.0%
0.0% 0.0%
776
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Table 4 Comparisons of different orders of DFT-MC Monte Carlo simulation results, 1000 samples, outlier identified by 3s, ( a, L c ) ¼ (0.6 mm, 321 mm). Modal order
Noise level ψ
Accuracy
False dismissal
Failure
1st 2nd 3rd 5th 7th 9th 10th
0.02 0.02 0.02 0.02 0.02 0.02 0.02
0.4% 13.3% 1.4% 82.2% 83.8% 0.6% 93.2%
68.9% 63.0% 93.5% 16.9% 16.2% 99.4% 6.8%
30.7% 23.7% 5.1% 0.9% 0.0% 0.0% 0.0%
3s. For the higher order DFT-MCs, one can observe the obvious improvements of damage identification results (as shown in Table 4 and Fig. 17). However, not all higher modes perform well which is indicated by the 9th order DFT-MC presented in Fig. 17(c). In application, one can hardly know which higher DFT-MC performs well, therefore, the identification obtained from different order DFT-MCs should be compared to avoid failure of damage identification. 3.4.2. Post-processing methods The characteristic of crack has been revealed in Fig. 16(d–f), but it is not strong enough to be classified as outlier. Hence, some post-processing methods could be added to intensify these characteristics. The curvature modes provided by the DFTMCs are employed as the inputs for the post-processing methods, which have been proved effective in combination with curvature modes. Then the new curves will be tested by statistical analysis to identify the outliers. In this paper, we will employ the continuous wavelet transform (CWT) as the post-processing method for DFT-MCs. Damage parameters are selected as ( a, L c ) ¼(0.6 mm, 321 mm), and the 5th, 7th, 9th and 10th orders mode shapes or curvature mode shapes are investigated due to the sensibility of higher mode shapes. Biorthogonal wavelet ‘bior6.8’ is selected in the CWT and the scale upper boundary is settled as 64. Identifications given by the CWT, the CD-MC and the present hybrid method (the DFT-MC and the CWT) for the same noise level ψ ¼0.02 are displayed in Fig. 18. It is observed in Fig. 18(a–d) that only boundary effects and some undulations induced by intact parts of mode shape are presented when the CWT is employed to decompose mode shape directly. In Fig. 18(e–h) one only can see some peaks relating to noise in low scale. However, the damage location can be easily found in Fig. 18(i–l) when the present hybrid method is used. The comparison validates the effectiveness of the hybrid method. For more stable and reasonable identification, statistical analysis is required. We propose to superpose the curves of all the scales to get a lumped curve which reflects all the information in the CWT, and enhances the characteristics of damage. Then outliers will be studied via the statistical analysis method mentioned above as shown in Fig. 19. It is seen that damage positions are clearly identified.
Fig. 16. The DFT-MCs “Accuracy” and “False dismissal” samples for the crack (a, Lc) ¼ (0.6 mm, 321 mm): (a) The 1st order DFT-MC with ψ ¼ 0.02, (b) The 2nd order DFT-MC with ψ ¼0.02, (c) The 3rd order DFT-MC with ψ ¼ 0.02, (d–f) The noise-free DFT-MC curves for the case shown in (a–c).
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
777
Fig. 17. The higher order DFT-MCs for the crack case (a, Lc) ¼ (0.6 mm, 321 mm) with ψ ¼ 0.02: (a) the 5th order, (b) the 7th order, (c) the 9th order, and (d) the 10th order.
L
L
L
L
L
L
L
L
L
L
Fig. 18. The CWT results for the crack case (a, Lc) ¼(0.6 mm, 321 mm) with ψ ¼ 0.02: (a–d) the results calculated via the 5th, 7th, 9th and 10th orders directly, (e–h) the results calculated via the 5th, 7th, 9th and 10th CD-MCs, (i–l) the results calculated via the 5th, 7th, 9th and 10th DFT-MCs.
To evaluate the accuracy, false dismissal and failure ratios of the hybrid method, the Monte Carlo simulation is employed again and the results are displayed in Table 5. By aid of the CWT, the present method gets a higher accuracy ratio and higher robustness. It should be noticed that not all the higher modes perform better than the lower ones (even with the aid of the CWT), for example the 9th order (shown in Table 5). The failure of the method is mainly caused by the boundary effect.
778
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 19. The superposition of all the scale CWT results and the outliers: (a–d) the results of 5th, 7th, 9th and 10th orders.
Table 5 Evaluation of the hybrid method, ( a, L c ) ¼ (0.6 mm, 321 mm). Case
Modal order
Noise level ψ
Accuracy
False dismissal
Failure
1 2 3 4 5 6 7 8
5th 5th 7th 7th 9th 9th 10th 10th
0.02 0.04 0.06 0.07 0.06 0.07 0.15 0.18
98.6% 70.6% 92.1% 75.4% 52.2% 45.9% 98.5% 77.4%
0.0% 0.0% 0.0% 0.0% 0.0% 0.2% 0.0% 0.0%
1.4% 29.4% 7.9% 24.6% 47.8% 53.9% 1.5% 22.6%
4. Experimental investigation The experimental investigation is based on the benchmark given by Katunin [36–38]. The benchmark can be obtained from http://ipkm.polsl.pl/wavstructdamas/, and the data we choose for the comparison are mode shapes measured for the case of beam with single crack and the baseline data (the intact structure). In this case, specimens (beam, a strip from a composite plate) are made by 24-layered glass-fiber reinforced plastic (GFRP) with unidirectional impregnated fibers and the [0°/60°/ 60°] lay-up. The dimensions of the specimens are as follows: length L ¼250 mm, width B ¼25 mm, and thickness H ¼5.28 mm. The crack is simulated via notch with the depth of 1 mm locating near 150 mm from the clamped end. The clamped length is about 20 mm. Based on the scanning LDV Polytec PSV-400, 44 sampling points are measured on each specimen to get the 1st–4th order mode shapes with a 215 mm effective length. Readers are suggested to refer to [36– 38] for more details. Firstly, the results obtained by the DFT-MC for intact (baseline) and damaged beams are presented in Fig. 20. Compared with the baseline, the identified results for the damaged beam show some singularities near 150 mm but not clear enough, especially in the 3rd and 4th order modes. Thus, the hybrid approaches are required. The damage assessments using the CD-MC, the CWT and the present hybrid method are all displayed in Figs. 21 and 22 for intact case and damaged case, respectively. The upper boundary of scale is settled as 6 considering the length of data. In assessments for intact case (Fig. 21), all three methods perform well, without any false alarm. It is observed that some small peaks and the undulations related to the waveform of mode shape in the results given by the CD-MC. The undulations related to the waveform can also be seen in the identifications generated by the DFT-MC (Fig. 21(i–l)), but not deemed as damage location. In the assessments obtained by the CWT directly (Fig. 21(e–h)), only some undulations related to the waveform and boundary effects exist. The upper and the lower bounds generated by statistical analysis in the present method provide a quantified criterion to classify the peaks and undulations, without which the normal characteristics in damage indexes may be deemed as the evidence of damage. In the damaged case, the CD-MC generates some obvious characteristics near damage location (Fig. 22(a–d)). However, these characteristics may be blurred by the irregular peaks due to numerical difference. In Fig. 22(e–h), the performances of the CWT (biorthogonal wavelet ‘bior6.8’) are better than the CD-MC, some ridges related to the damage can be easily found. However, the maximum values of these ridges appear near 160–170 mm instead of 150 mm. This localization error can be interpreted by the inappropriate choice of wavelet and parameters. It should be emphasized that the dependence of selection of wavelet type and parameters may restrict the application of the CWT in damage detection. In Fig. 22(i–l), the last assessments via present hybrid method is performed. Damage locations can be easily identified from these figures and extracted by aid of the indications of outliers.
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
779
Fig. 20. The 1st–4th order DFT-MCs of intact and cracked structures (a–d) for the intact structure, (e–h) for the cracked structure.
0.1
0.01
0.05
0.02
0.05
0
0
-0.01
-0.05
50 100 150 200
50 100 150 200
-0.1
50 100 150 200
L [mm]
L [mm] 6
5
5
4
4
4
4
50
100
150
2
1
200
3
50
100
150
1
200
50
L [mm]
L [mm] 4
100
150
1
200
-4
-10 50
100 150
L [mm]
200
w"
w"
-5
150
200
L [mm]
5
0
-2
100
10
5
0
50
L [mm]
10
2
3 2
5
w"
1
3 2
2
Scale
6
5
Scale
6
5
3
50 100 150 200
L [mm]
6
Scale
Scale
L [mm]
w"
0 -0.05
-0.02
-0.02
w"
w"
w"
w"
0
0 -5
50
100 150 200
L [mm]
0 -5
50
100 150
L [mm]
200
50
100 150 200
L [mm]
Fig. 21. Identifications using the intact structure benchmark data: (a–d) the 1st–4th order CD-MCs, (e–h) the 1st–4th order CWTs, (i–l) the results of the present hybrid method.
780
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
Fig. 22. Identifications using cracked structure benchmark data: (a–d) the 1st–4th order CD-MCs, (e–h) the 1st–4th order CWTs, (i–l) the results of the present hybrid method.
5. Conclusions Modal curvature is an important dynamic feature in structural health monitoring. However, the susceptibility to noise induced by numerical derivative severely impairs the applicability in damage detection. To tackle this problem, the DFT-MC is proposed and some auxiliary methods are introduced in this paper. By aid of the characteristics of the FT, the numerical derivative in the CD-MC is accomplished by multiplication and the DFT synthesis, which enhances the numerical stability of modal curvature estimation. In addition, the wavenumber domain in the DFT provides an effective method for space filtering. Based on the iteration, statistical analysis and the post-processing method, a robust hybrid algorithm is generated. It should be mentioned that the present method is just a modification for estimating MC. Restricted by the fundamental properties of MC, the capacity of present method for small damage detection may be weak in applications. The singularity related to damage presence can be easily covered by noise. In that case, it appears few methods can conduct the satisfactory noise suppression and damage identification meanwhile based on limited time-domain measurement of mode shapes. Acknowledgments Special thanks should be given to Prof. A. Katunin, who shares the benchmark data in the website. This work is supported by the National Natural Science Foundation of China (No. 51405369), the China Scholarship Council, and the China Postdoctoral Science Foundation (No. 2014M560766).
References [1] S.P. Lele, S.K. Maiti, Modelling of transverse vibration of short beams for crack detection and measurement of crack extension, J. Sound Vib. 257 (2002)
Z.-B. Yang et al. / Mechanical Systems and Signal Processing 84 (2017) 763–781
781
559–583. [2] S.M. Murigendrappa, S.K. Maiti, H.R. Srirangarajan, Frequency-based experimental and theoretical identification of multiple cracks in straight pipes filled with fluid, Ndt E Int. 37 (2004) 431–438. [3] M.R. Naniwadekar, S.S. Naik, S.K. Maiti, On prediction of crack in different orientations in pipe using frequency based approach, Mech. Syst. Signal Process. 22 (2008) 693–708. [4] Z. Su, L. Ye, Y. Lu, Guided Lamb waves for identification of damage in composite structures: a review, J. Sound Vib. 295 (2006) 753–780. [5] Z. Su, L. Cheng, X. Wang, L. Yu, C. Zhou, Predicting delamination of composite laminates using an imaging approach, Smart Mater. Struct. 18 (2009) 074002. [6] M. Cao, L. Ye, L. Zhou, Z. Su, R. Bai, Sensitivity of fundamental mode shape and static deflection for damage identification in cantilever beams, Mech. Syst. Signal Process. 25 (2011) 630–643. [7] X. Qian, M. Cao, Z. Su, J. Chen, A hybrid particle swarm optimization (PSO)-simplex algorithm for damage identification of delaminated beams, Math. Probl. Eng. 2012 (2012) 607418. [8] Z.-B. Yang, X.-F. Chen, Y. Xie, H. Zuo, H.-H. Miao, X.-W. Zhang, Wave motion analysis and modeling of membrane structures using the wavelet finite element method, Appl. Math. Model. 40 (2016) 2407–2420. [9] Z. Yang, X. Chen, X. Li, Y. Jiang, H. Miao, Z. He, Wave motion analysis in arch structures via wavelet finite element method, J. Sound Vib. 333 (2014) 446–469. [10] A.K. Pandey, M. Biswas, M.M. Samman, Damage detection from changes in curvature mode shapes, J. Sound Vib. 145 (1991) 321–332. [11] Z.B. Yang, X.F. Chen, S.H. Tian, Z.J. He, Multiple damages detection in beam based approximate waveform capacity dimension, Struct. Eng. Mech. 41 (2012) 663–673. [12] M. Cao, W. Xu, W. Ostachowicz, Z. Su, Damage identification for beams in noisy conditions based on Teager energy operator-wavelet transform modal curvature, J. Sound Vib. 333 (2014) 1543–1553. [13] M. Cao, M. Radzieński, W. Xu, W. Ostachowicz, Identification of multiple damage in beams based on robust curvature mode shapes, Mech. Syst. Signal Process. 46 (2014) 468–480. [14] Z.-B. Yang, M. Radzienski, P. Kudela, W. Ostachowicz, Two-dimensional modal curvature estimation via Fourier spectral method for damage detection, Compos. Struct. 148 (2016) 155–167. [15] W. Fan, P.Z. Qiao, A strain energy-based damage severity correction factor method for damage identification in plate-type structures, Mech. Syst. Signal Process. 28 (2012) 660–678. [16] M.S. Cao, P.Z. Qiao, Novel Laplacian scheme and multiresolution modal curvatures for structural damage identification, Mech. Syst. Signal Process. 23 (2009) 1223–1242. [17] J. Xiang, T. Matsumoto, Y. Wang, Z. Jiang, Detect damages in conical shells using curvature mode shape and wavelet finite element method, Int. J. Mech. Sci. 66 (2013) 83–93. [18] R.B. Bai, M.S. Cao, Z.Q. Su, W. Ostachowicz, H. Xu, Fractal dimension analysis of higher-order mode shapes for damage identification of beam structures, Math. Probl. Eng. 2012 (2012) 454568. [19] M. Radzienski, M. Krawczuk, M. Palacz, Improvement of damage detection methods based on experimental modal parameters, Mech. Syst. Signal Process. 25 (2011) 2169–2190. [20] M. Radzienski, M. Krawczuk, W. Ostachowicz, Experimental verification and comparison of mode shape-based damage detection methods, Damage Assess. Struct. VIII 413-414 (2009) 699–706. [21] A. Katunin, P. Przystałka, Automated wavelet-based damage identification in sandwich structures using modal curvatures, J. Vibroeng. 17 (2015). [22] Z. Yang, X. Chen, J. Yu, R. Liu, Z. Liu, Z. He, A damage identification approach for plate structures based on frequency measurements, Nondestruct. Test. Eval. (2013) 1–21. [23] X. Zhang, R.X. Gao, R. Yan, X. Chen, C. Sun, Z. Yang, Multivariable wavelet finite element-based vibration model for quantitative crack identification by using particle swarm optimization, J Sound Vib. 375 (2016) 200–216. [24] J. Xiang, M. Liang, Y. He, Experimental investigation of frequency-based multi-damage detection for beams using support vector regression, Eng. Fract. Mech. 131 (2014) 257–268. [25] Z. Yang, X. Chen, Y. Jiang, Z. He, Generalised local entropy analysis for crack detection in beam-like structures, Nondestruct. Test. Eval. 29 (2014) 133–153. [26] M. Hong, Z. Su, Y. Lu, H. Sohn, X. Qing, Locating fatigue damage using temporal signal features of nonlinear Lamb waves, Mech. Syst. Signal Process. 60 (2015) 182–197. [27] Z. Su, C. Zhou, M. Hong, L. Cheng, Q. Wang, X. Qing, Acousto-ultrasonics-based fatigue damage characterization: linear versus nonlinear signal features, Mech. Syst. Signal Process. 45 (2014) 225–239. [28] B.P. Nandwana, S.K. Maiti, Detection of the location and size of a crack in stepped cantilever beams based on measurements of natural frequencies, J. Sound. Vib. 203 (1997) 435–446. [29] Z.B. Yang, X.F. Chen, Y. Xie, X.W. Zhang, The hybrid multivariate analysis method for damage detection, Struct. Control Health Monit. 23 (2016) 123–143. [30] Z.B. Yang, X.F. Chen, Y. Xie, H.H. Miao, J.J. Gao, K.Z. Qi, Hybrid two-step method of damage detection for plate-like structures, Struct. Control Health Monit. 23 (2016) 267–285. [31] A. Deraemaeker, A. Preumont, Vibration based damage detection using large array sensors and spatial filters, Mech. Syst. Signal Process. 20 (2006) 1615–1630. [32] G. Tondreau, A. Deraemaeker, Numerical and experimental analysis of uncertainty on modal parameters estimated with the stochastic subspace method, J. Sound Vib. 333 (2014) 4376–4401. [33] G. Tondreau, A. Deraemaeker, Automated data-based damage localization under ambient vibration using local modal filters and dynamic strain measurements: Experimental applications, J. Sound Vib. 333 (2014) 7364–7385. [34] G. Tondreau, A. Deraemaeker, Local modal filters for automated data-based damage localization using ambient vibrations, Mech. Syst. Signal Process. 39 (2013) 162–180. [35] G. Tondreau, A. Deraemaeker, Experimental Localization of Small Damages using Modal Filters, Special Topics in Structural Dynamics, Springer, Berlin, Germany (2013), pp. , 2013, 585–591. [36] A. Katunin, Nondestructive damage assessment of composite structures based on wavelet analysis of modal curvatures: state-of-the-art review and description of wavelet-based damage assessment benchmark, Shock Vib. 501 (2015) 735219. [37] A. Katunin, K. Dragan, M. Dziendzikowski, Damage identification in aircraft composite structures: a case study using various non-destructive testing techniques, Compos. Struct. 127 (2015) 1–9. [38] A. Katunin, Diagnostics of Composite Structures using Wavelets, Radom, Gliwice, Poland, 2015.