A methodology for identification of damage in beams

Abstract

In this work, a parameter based on correlation of signals using mathematical morphology and other two parameters based on energy of a vibration signal using wavelet transform and bispectrum theory are used for approximate location of damage in a steel clamped-free beam. The experimental data are obtained through accelerometers placed along the sample. The system is excited using impact hammer. A preliminary location of the damage is found using the position of the accelerometer that has a higher energy change when it is compared to the signal of the system with and without damage. The mathematical models are obtained using 2D elasticity theory and the finite element method. The numerical and experimental data are approximated using the particle swarm optimization method, and this way, it is possible to adjust the location and severity of the damage. The procedure is also applied to a free–free steel beam with two damages and a sandwich beam.

Keywords:

• damage
• beam
• PSO
• bispectrum
• wavelet transform
• mathematical morphology
• FEM

AMS Subject Classifications:

• 70J10
• 70J30
• 74H45
• 74R99

1. Introduction

An area of great interest in the field of engineering is the study and identification of structural damage. There are a lot of mechanical tests and techniques used to detect cracks and possible failures in mechanical structures. One of the major problems for the industry lies on acquiring reliable and non-destructible results.

Most identification techniques are based on modal parameters measured using only a few modes of vibration and/or modal frequencies of the structure that can easily be obtained by dynamic testing. The results are obtained with dynamic measurements at different times, corresponding to two moments of the structure lifetime (usually with the model in good condition and damaged).

The main methods based on vibration signals are as follows: mode shape curvature method,[1,2] change in flexibility method,[1,2] change in flexibility curvature method,[1,2] strain energy method,[1,3] statistical methods,[4–11] wavelet analysis [12–16] and response surface methodology.[17] Defects in components of machinery and structures can be detected by monitoring vibration.[6]

The bispectrum, a third-order statistic and kurtosis, a fourth-order moment, help to identify faults in mechanical components. The bispectrum technique relates one set of mixing waves through the spectral coupling.[4–10]

Mathematical morphology (MM) is a nonlinear analysis method, which has been developed and widely applied to various fields of image processing and analysis. When MM is used in signal processing, the information of local morphological features of the signal is the only determinative factor. By morphological signal decomposition, a complex signal can be separated from the background and decomposed into various components reserving morphological features of the signal. Literature review reveals that the research of morphology analysis in one-dimensional (1-D) signals, especially in machine fault diagnosis, is still limited.[18,19]

Identification, localization and quantification of structural damage can be performed through a model updating procedure. Model updating methods require a baseline finite element model of the undamaged structure, which imposes a restriction on their applicability and can become very problematic especially for large and complex structures.[20]

The correlation techniques are a mixture of visual and numerical means to identify the differences between measurements and predictions. Whereas numerical correlation techniques return a numerical value, visual means of correlation are subjective and of qualitative nature.[21] Some of the basic correlation tools include simple tabulation or plotting the measured and predicted eigenvalues. A more strict correlation is the use of the so-called modal assurance criteria (MAC).[20–22]

When materials in a structure suffer damage, the mechanical response of the structure is modified, particularly in the vicinity of the damaged zones. Such local perturbations in the structural response are generally very small, but they can be detected using wavelet transform techniques.[16]

Recently, wavelet analysis has shown its inherent merits in damage detection over traditional methods due to its ability to closely examine the signal with multiple scales to provide various levels of detail and approximations. The use of wavelet transform to identify damage from mode shape has been one of the most popular techniques. These methods treat mode shape data as a signal in spatial domain, and they use spatial wavelet transform technique to detect the signal irregularity caused by damage. A damage detection index called wavelet packet energy rate index (WPERI) was proposed by [15] for the damage detection of beam structures. The measured dynamic signals are decomposed into the wavelet packet components, and the wavelet energy rate index is computed to indicate the structural damage.

The aim of this work was to analyse structural changes through the inclusion of damage (crack) with known formats and positions in the samples. The damage in steel clamped-free beam is analysed considering samples with and without structural modifications. The damage location is investigated by varying the damage position in the mathematical model and comparing the numerical and experimental data. The mathematical models are obtained using the finite element method (FEM). The experimental data are obtained using the impact hammer and an accelerometer placed in different positions along the sample. The statistics procedures are used in qualitative analysis. In order to approximate the experimental and numeric frequency response functions (FRF) data, the particle swarm optimization (PSO) [23–25] method is used. The first approximations of the damage position are obtained using three different parameters based on the energy of the signal and using morphological analysis, bispectrum and wavelet analysis. Using this procedure and making a sweep of the modified finite element position along the sample, it is possible to identify, with great precision, the location and severity of the damage through the comparison of the experimental and numeric modal parameters. The procedure is also applied to a free–free steel beam with two damages and a sandwich beam.

2. Mathematical models

2.1. Two-dimensional elastic model

The two-dimensional (2-D) elastic analysis of the sandwich beam was performed by assuming that the materials are homogeneous and present a linear behaviour and that the materials are perfectly bonded in the interfaces. The continuity of the displacements at the interfaces of the materials that make up the sandwich beam (continuous and fully glued interfaces) allows the finite elements to assemble in the conventional way, even with different materials.

In this work, cubic Lagrangian Finite Element (16 nodes) is used for the analyses and the elastic components of the displacement u = (u, v) can be interpolated as follows:(1) $$\left\{\begin{array}{c}u(\mathbf{x})\\ v(\mathbf{x})\end{array}\right\}\cong \left[\begin{array}{ccccc}{\mathit{\varphi }}_1(\mathbf{x})& 0& \dots & {\mathit{\varphi }}_{16}(\mathbf{x})& 0\\ 0& {\mathit{\varphi }}_1(\mathbf{x})& \dots & 0& {\mathit{\varphi }}_{16}(\mathbf{x})\end{array}\right]\{\mathbf{q}\}=[\mathrm{\Phi }(\mathbf{x})]\{\mathbf{q}\}$$(1)

where $\{\mathbf{q}\}={\left\{\begin{array}{ccccc}{u}_1& v_1& .....& u_{16}& .\end{array}{v}_{16}\right\}}^T$ are the pairs $\left(\begin{array}{cc}{u}_i& v_i\end{array}\right)$ that represent the displacement components at node i of the element and ϕ_i(x) represents the Lagrangian interpolation functions (well known in the classical books of finite elements).

The linear strains vector is related to the generalised displacement vector in the following manner:(2) $$\{\mathit{\epsilon }\}=\left\{\begin{array}{c}{\mathit{\epsilon }}_{xx}\\ {\mathit{\epsilon }}_{yy}\\ {\mathit{\gamma }}_{xy}\end{array}\right\}=\left\{\begin{array}{c}\mathit{\partial }u(\mathbf{x})/\mathit{\partial }x\\ \mathit{\partial }v(\mathbf{x})/\mathit{\partial }y\\ \mathit{\partial }u(\mathbf{x})/\mathit{\partial }y+\mathit{\partial }v(\mathbf{x})/\mathit{\partial }x\end{array}\right\}\cong [B]\{\mathbf{q}\}$$(2)

where [B] is a (3 × 32) matrix relating the three strains to the generalised displacements and is given by:(3) $$[B]=\left[\begin{array}{ccccc}\frac{\mathit{\partial }{\mathit{\varphi }}_1}{\mathit{\partial }x}& 0& ...& \frac{\mathit{\partial }{\mathit{\varphi }}_{16}}{\mathit{\partial }x}& 0\\ 0& \frac{\mathit{\partial }{\mathit{\varphi }}_1}{\mathit{\partial }y}& ...& 0& \frac{\mathit{\partial }{\mathit{\varphi }}_{16}}{\mathit{\partial }y}\\ \frac{\mathit{\partial }{\mathit{\varphi }}_1}{\mathit{\partial }y}& \frac{\mathit{\partial }{\mathit{\varphi }}_1}{\mathit{\partial }x}& ...& \frac{\mathit{\partial }{\mathit{\varphi }}_{16}}{\mathit{\partial }y}& \frac{\mathit{\partial }{\mathit{\varphi }}_{16}}{\mathit{\partial }x}\end{array}\right]$$(3)

The element stiffness matrix, [K] is obtained by integration and its expression is found in the conventional finite element classical books:(4) $$[K]=\underset{A}{\int }{[B]}^t[D][B]t\text{d}A$$(4)

where t denotes the thickness of the element, A denotes the area of the element and [D] is the (3 × 3) matrix, which relates three stresses and three strains,(5) $$\{\mathit{\sigma }\}=\left\{\begin{array}{c}{\mathit{\sigma }}_{xx}\\ {\mathit{\sigma }}_{yy}\\ {\mathit{\tau }}_{xy}\end{array}\right\}=[D]\left\{\begin{array}{c}{\mathit{\epsilon }}_{xx}\\ {\mathit{\epsilon }}_{yy}\\ {\mathit{\gamma }}_{xy}\end{array}\right\}=[D]\{\mathit{\epsilon }\}$$(5)

The [D] matrix is written as a function of the elastic material parameters, and the adjustment of these parameters as a function of frequency was the main objective of this work. The mass matrix is obtained using the following integral:(6) $$[M]={\int }_A\mathit{\rho }{\left[\mathrm{\Phi }\right]}^t\left[\mathrm{\Phi }\right]t\text{d}A$$(6)

where ρ is the material density and [Φ] is a matrix of the order (2×32) of Equation (1(1) $$\left\{\begin{array}{c}u(\mathbf{x})\\ v(\mathbf{x})\end{array}\right\}\cong \left[\begin{array}{ccccc}{\mathit{\varphi }}_1(\mathbf{x})& 0& \dots & {\mathit{\varphi }}_{16}(\mathbf{x})& 0\\ 0& {\mathit{\varphi }}_1(\mathbf{x})& \dots & 0& {\mathit{\varphi }}_{16}(\mathbf{x})\end{array}\right]\{\mathbf{q}\}=[\mathrm{\Phi }(\mathbf{x})]\{\mathbf{q}\}$$(1) ).

For punctual external forces as f(x, t) = f₀ sin (ωt) acting on node i and the balance of the element in relation to the harmonic analysis can be written as:(7) $$[K_e(E,\mathit{\omega })]\{\mathbf{q}\}=\{\mathbf{f}\}$$(7)

where(8) $$[K_e(E,\mathit{\omega })]=[K(E)]-{\mathit{\omega }}^2[M]$$(8)

and(9) $$\{f\}={\left\{\begin{array}{ccc}0...& f_0& ...0\end{array}\right\}}^t$$(9)

and the position is not equal to zero at 2i when the force is applied in the y direction and 2(i − 1) when the force is applied in the x direction.

Figure 1 shows the detail of how the boundary condition of type clamped end was imposed on the sandwich beam. As can be seen, this condition is much better represented with 2-D elements, because null displacements are imposed only on the nodes in contact with the rigid base (unlike the beam element).

Figure 1. Detail of the boundary conditions for the clamped end.

Figure 2 shows the pattern construction of the homogeneous finite element mesh (it shows only a 3 mm thickness band along the height of the beam). Note that elements with a length of 3 mm were used to keep the aspect ratio limited to 5 for the elements of the steel sheet. This formulation is also used to simulate the dynamical behaviour of a steel clamped beam. In this case, the model is simplified.

Figure 2. Construction pattern mesh along the thickness of the beam (rotated 90°) and cubic Lagrangian element used in the analysis.

2.2. Bispectrum theory

The statistical properties of a stationary random process are completely described by its mean value μ (first order moment) and variance σ² (second-order central moment) if and only if its probability density function has a Gaussian distribution. Higher order statistics (HOS) measures, such as higher order moments, are extensions of second-order measures to higher orders and cannot provide additional information about a signal if it is Gaussian. However, many signals found in the industrial field are non-Gaussian, e.g. vibrational signals of rotating machines. Therefore, HOS may be used to extract information about signals and systems which cannot be obtained from conventional statistics. The most used traditional signal processing measure is the power spectrum, that is the decomposition over frequency of the signal power and is therefore related to the signal variance σ². The bispectrum (third-order spectrum) can be viewed as a decomposition of the third moment (skewness) of a signal over frequency and as such can detect non-symmetric nonlinearities. For a stationary random process, the discrete bispectrum B(k, l) can be defined in terms of the signal’s discrete Fourier transform X(k) as:

(10) $$B(k,l)=E⌊X(k)X(l)X{(k+l)}^{\ast }⌋$$(10)

where $E\lfloor \rfloor$ denotes the expectation operator. It should be noted that the bispectrum is complex-valued (it contains phase information) and that it is a function of two independent frequencies, k and l. Furthermore, it is not necessary to compute B(k, l) for all (k, l) pairs, due to several symmetries existing in the (k, l) plane. If ω_k, ω_l and ω_k + l are independent, each one will have an independent random phase (relative to each other). Therefore, the bispectrum is a statistical is of great importance in the study of those nonlinear vibrations where the relationships between three spectral components are in question, such as generation of combinational resonance modes or quadratic mode couplings [4].

2.3. Wavelet theory

The continuous wavelet transform is defined as follows [15]:(11) $$C(a,b)=\underset{-\infty }{\overset{+\infty }{\int }}f(t){\mathit{\psi }}_{a,b}(t)\text{d}t$$(11)

where(12) $${\mathit{\psi }}_{a,b}(t)=a^{1/2}\mathit{\psi }\left(\frac{t-b}{a}\right)$$(12)

is a window function called the mother wavelet, where a is a scale and b is a translation. The term wavelet means a small wave. The smallness refers to the condition in which this (window) function is in the finite length (compactly supported). The wave refers to the oscillatory condition of this function. The term mother implies that the functions with different support regions that are used in the transformation process are derived from one main function, or the mother wavelet. In other words, the mother wavelet is a prototype for generating the other window functions.

Wavelet packets consist of a set of linearly combined usual wavelet functions. The wavelet packets inherit the properties such as orthonormality and time–frequency localization from their corresponding wavelet functions. A wavelet packet ${\mathit{\psi }}_{j,k}^i(t)$ is a function with three indices where integers i, j and k are the modulation, scale and translation parameters, respectively,

(13) $${\mathit{\psi }}_{j,k}^i(t)=2^{j/2}{\mathit{\psi }}^j(2^{j}t-k),i=1,2,3,\dots$$(13)

The wavelet packet component signal $f_j^i(t)$ can be represented by a linear combination of wavelet packet functions ${\mathit{\psi }}_{j,k}^i(t)$ as follows:

where the wavelet packet coefficients $c_{j,k}^i(t)$ can be obtained from(15) $$c_{j,k}^i(t)=\underset{-\infty }{\overset{+\infty }{\int }}f(t){\mathit{\psi }}_{j,k}^i(t)\text{d}t$$(15)

and(16) $$f(t)=\sum _{i=1}^{2j}{f}_j^i(t)$$(16)

The wavelet packet energy index is used to identify the initial locations of damage. In this case, the signal energy $E_{f_j}$ at j level is first defined as:(17) $$E_{f_j}=\underset{-\infty }{\overset{+\infty }{\int }}{f}^2(t)\text{d}t=\sum _{m=1}^{2^j}\sum _{n=1}^{2^j}\underset{-\infty }{\overset{+\infty }{\int }}{f}_j^m(t)f_j^n(t)\text{d}t$$(17)

After manipulations, it is possible to find the wavelet packet component $E_{f_j^i}$ as the stored energy in the component signal $f_j^i(t)$ as:(18) $$E_{f_j^i}=\underset{-\infty }{\overset{+\infty }{\int }}{f}_j^i{(t)}^2\text{d}t$$(18)

The WPERI is used to indicate the localization of the structural damage. The rate of signal wavelet packet energy $\mathrm{\Delta }(E_{f_j})$ at j level is defined as(19) $$\mathrm{\Delta }(E_{f_j})=\sum _{i=1}^{2^j}\frac{\left|{(E_{f_j^i})}_b-{(E_{f_j^i})}_a\right|}{{(E_{f_j^i})}_a}$$(19)

where ${(E_{f_j^i})}_a$ is the component signal energy $E_{f_j^i}$ at j level without damage, and ${(E_{f_j^i})}_b$ is the component signal energy $E_{f_j^i}$ with some damage.

2.4. Mathematical morphology

The vibration signal dealt with in this study is a discrete 1-D signal, and the multivalued morphological transformation for this type of signal is addressed in following paragraphs.[18,19]

If f(n) is the original 1-D signal, which is the discrete function over a domain F = (0; 1; 2; … ; N − 1) and g(n) is the SE (flat structuring element), which is the discrete function over a domain G = (0; 1; 2; …; M − 1), two basic morphological operators, the erosion and the dilation, can be defined as:(20) $$\left(f\ominus g\right)\left(n\right)=\text{min}\left[f\left(n+m\right)-g\left(m\right)\right],m\in 0,1,2,\dots ,M-1$$(20) (21) $$(f\oplus g)\left(n\right)=\text{max}\left[f\left(n+m\right)-g\left(m\right)\right],m\in 0,1,2,\dots ,M-1$$(21)

where $\ominus$ denotes the erosion operator and ⊕ denotes the dilation operator.

Based on the dilation and erosion, two other basic morphological operators, the opening and the closing, can be further defined as:(22) $$(f\circ g)\left(n\right)=(f\ominus g\oplus g)\left(n\right)$$(22) (23) $$(f\bullet g)\left(n\right)=(f\oplus g\ominus g)\left(n\right)$$(23)

where $\circ$ stands for the opening operator and $\bullet$ for the closing operator.

These four morphological operators can all be used to extract morphological features of a signal, but different operators fit different morphological features.

Multi-scale MM refers to a morphology analysis with SEs at different scales. The SE scale, especially the length scale, is important for the multi-scale morphology analysis of 1-D signals. The multi-scale morphological operations have also opening and closing operations. Table 1 shows the double-dot structuring element used in the analysis. The correlation coefficient of two pattern spectra is expressed as(24) $$\mathit{\rho }=\frac{\text{Cov}\left[P_1,P_2\right]}{\sqrt{\text{Var}\left[P_1\right]\text{Var}\left[P_2\right]}}$$(24)

Table 1. Multi-scale structuring elements.

Scale	Double-dot structuring element
1	{1 0 1}
2	(1 0 0 1}
3	{1 0 0 0 1}
4	(1 0 0 0 0 1}
n	…

where P₁ and P₂ represent two different pattern spectra, and ρ is their correlation coefficient which measures the similarity of two signals (with and without damage).

2.5. PSO method

In this work, the physical parameters of a steel clamped-free beam are estimated using measured and numeric FRFs. The mathematical models are obtained using the FEM. In order to approximate the experimental and numeric FRF data, the PSO method is used. The PSO algorithm is a biologically inspired algorithm motivated by a social analogy. Sometimes, it is related to the evolutionary computation techniques, basically with the genetic algorithms and evolutionary strategies, but there are significant differences among those techniques.[23–25]

The PSO algorithm is population based: a set of potential solutions evolves to approach a convenient solution (or set of solutions) for a problem. Being an optimization method, the aim is finding the global optimum of a real-valued function (fitness function) defined in a given space (search space).

In the PSO algorithm, each individual is called a ‘particle’ and is subject to a movement in a multidimensional space that represents the belief space. Particles have memory, thus retaining part of their previous state. There is no restriction for particles to share the same point in belief space, but in any case, their individuality is preserved. Each particle’s movement is the composition of an initial random velocity and two randomly weighted influences: individuality, the tendency to return to the particle’s best previous position, and sociality, the tendency to move towards the neighbourhood’s best previous position.

The ‘continuous’ version uses a real-valued multidimensional space as belief space and evolves the position of each particle in that space using the following equations:(25) $$v_i(k+1)=v_i(k)+{\mathit{\gamma }}_{1i}(p_i-x_i(k))+{\mathit{\gamma }}_{2i}(G-x_i(k))$$(25) (26) $$x_i(k+1)=x_i(k)+v_i(k+1)$$(26)

where i is the particle index; k is the discrete time index; v is the velocity on ith particle; x is the position of the ith particle; p is the best position found by ith particle (personal best); G is the best position found by swarm (global best, best of personal bests); γ_1,2 are random numbers on the interval [0,1] applied to ith particle

The common PSO algorithm is as follows:(27) $$v_i(k+1)=\mathit{\varphi }(k)v_i(k)+{\mathit{\alpha }}_1\left[{\mathit{\gamma }}_{1i}(p_i-x_i(k))\right]+{\mathit{\alpha }}_2\left[{\mathit{\gamma }}_{2i}(G-x_i(k))\right]$$(27)

where ϕ is the inertia function, and α₁ and α₂ are the acceleration constants.

2.6. Modal parameters

Methods that consider the modal changes are also used in the identification of damage (structural changes). A method of setting modal (MAC) is used to compare pairs of modes. The matrix coefficients are obtained by [26]:(28) $${\text{MAC}}_{i,j}=\frac{{\left({\left\{{\mathit{\varphi }}_A\right\}}_i^T{\left\{{\mathit{\varphi }}_A\right\}}_j^{\ast }\right)}^2}{{\left\{{\mathit{\varphi }}_A\right\}}_i^T{\left\{{\mathit{\varphi }}_A\right\}}_i^{\ast }{\left\{{\mathit{\varphi }}_X\right\}}_j^T{\left\{{\mathit{\varphi }}_X\right\}}_j^{\ast }}$$(28)

where $\left\{{\mathit{\varphi }}_A\right\}$ and $\left\{{\mathit{\varphi }}_X\right\}$ denote the numeric and experimental modes; the superscript symbol T denotes the transpose of a vector and the subscript superscript symbol * denotes the complex conjugate vector.

A MAC value close to 1 suggests that the two modes are well correlated. An overall mode shape indicator may be calculated from:(29) $${\mathit{\epsilon }}_{\mathit{\varphi }}=\left[1-\frac{1}{n}\sqrt{\sum _{i=1}^n{(\text{MAC})}_i^2}\right]\times 100$$(29)

where n is the number of modes.

3. Results

Figure 3 shows the flowchart of the damage identification method with four steps.

Figure 3. Flowchart of the damage identification method.

3.1. Experiment

The geometric dimensions of the experimental clamped-free beam steel sample are as follows: length = 0.84 m; thickness = 0.0127 m and width = 0.0254 m. The effective length for the clamped free is 0.825 m.

The experimental data are obtained using the impact hammer, and one accelerometer displaced along the sample in the positions 0.125, 0.225, 0.325, 0.425, 0.525, 0.625, 0.725 and 0.825 m. The excitation force was applied in the same position (between first and second accelerometers) for all tests. Figure 4 shows schematically the basic components of the vibration testing system. The system consists of a rigid structure (1), accelerometers PCB 333B (2), impact hammer PCB 086C03 (3), beam sample (4), rigid termination (5), signal analyzer HP 3567 and a microcomputer (7).

Figure 4. Schematic view of the testing system for beam sample.

The structural modifications (Figure 5) were introduced in positions near 0.075 m (position 1) and 0.245 m (position 2). For these two cases, the data were collected for the system with and without damage with three levels of damage (damage 1, damage 2 and damage 3) whose geometrical parameters for both cases are as follows: length of 0.001 m and heights of 0.001, 0.0025 and 0.0045 m. A1–A8 represents the positions of the accelerometers.

Figure 5. Schematic view beam sample with damages.

The first step was to obtain data in the frequency domain and time domain. The data in the frequency domain are used to obtain the modal parameters, and the data in the time domain are used for the statistical methods application. The inverse of FRF was used aiming for the signals to have the same excitation force amplitude.

3.2. Visualisation

In an attempt to identify the damage, we applied the bispectrum technique in the signal in time domain. Figures 6 and 7 show the bispectrum of the signal of the accelerometer placed at position 0.425 m for the system with and without damage. The curves present differences, but it is visually difficult to establish a reliable parameter for analysis.

Figure 6. Bispectrum (system without damage).

Figure 7. Bispectrum (system with damage).

3.3. Damage position

It was defined a relative parameter that contained the sum of all values of the bispectrum (spectral energy) defined by (30(30) $$IndB=\sum _{k=1}^n\sum _{l=1}^{n}B(k,l)$$(30) ) and (31(31) $$B_{\text{rel}}=\frac{\left|Ind{B}_b-Ind{B}_a\right|}{\left|Ind{B}_a\right|}$$(31) ):(30) $$IndB=\sum _{k=1}^n\sum _{l=1}^{n}B(k,l)$$(30) (31) $$B_{\text{rel}}=\frac{\left|Ind{B}_b-Ind{B}_a\right|}{\left|Ind{B}_a\right|}$$(31)

where IndB_b is the bispectrum parameter of the system with damage; IndB_a is the bispectrum parameter of the system without damage; n is the number of points in the frequency range (in this case n = 400); k and l are the frequency points varying from 1 to 400 (increment frequency of 1 Hz).

Figures 8 and 9 show the values of the parameter B_rel for all positions of the accelerometer. Figure 8 shows the curves for the system with damage in position 1 and Figure 9 for the system with damage in position 2. It can be noticed that the great values for all cases of damage are found for the position 1 of the accelerometer (Figure 8) and position 2 of the accelerometer (Figure 9). In all situations, the great values correspond to the real near position of the damages.

Figure 8. Bispectrum index (damage in position 1).

Figure 9. Bispectrum index (damage in position 2).

Similar results were found using the energy index based on the wavelet transform signal. Figures 10 and 11 show the energy parameter values obtained by (19).

Figure 10. Energy index (damage in position 1).

Figure 11. Energy index (damage in position 2).

Finally, a correlation coefficient of signals was used taking into account the pattern spectrum of two signals, with and without damage. The results obtained with this coefficient were more stable (lower fluctuations) and, in this case, the lower values were near the damage position (Figures 12 and 13).

Figure 12. Average correlation coefficient (damage in position 1).

Figure 13. Average correlation coefficient (damage in position 2).

3.4. Damage parameters

After verifying the approximated position of the damage, a mathematical beam model was obtained through the FEM. For all analyses were written subroutines in Fortran and Matlab languages. The influence of the severity degree of damage was also evaluated. The MAC number (28) and the error parameter (29) were used to validate the numeric and experimental data. Tables 2 and 3 show the values of these parameters for the beam system with three different damage levels. The maximum error is of the order of 0.08% for the sample with damage in position 1 and 0.16% for the sample with damage in position 2.

Table 2. MAC number (damage in position 1).

System	1st mode shape	2nd mode shape	3rd mode shape	Error
No damage	0.9997	0.9987	0.9990	0.0855
Damage 1	0.9996	0.9987	0.9994	0.0767
Damage 2	0.9999	0.9993	0.9991	0.0571
Damage 3	0.9998	0.9985	0.9989	0.0938

Table 3. MAC number (damage in position 2).

System	1st mode shape	2nd mode shape	3rd mode shape	Error
No damage	0.9999	0.9996	0.9993	0.0413
Damage 1	0.9998	0.9994	0.9960	0.1621
Damage 2	0.9999	0.9997	0.9977	0.0911
Damage 3	0.9993	0.9982	0.9989	0.1198

Finally, we tried to adjust the possible location of the damage based on the natural frequencies variations and a performance index defined by Equation (32(32) $${\mathit{\epsilon }}_{\mathit{\omega }}={\left[\frac{{\sum }_{i=1}^L{\left(\mathrm{\Delta }{\mathit{\omega }}_{A_i}-\mathrm{\Delta }{\mathit{\omega }}_{X_i}\right)}^2}{{\sum }_{i=1}^L\mathrm{\Delta }{\mathit{\omega }}_{A_i}^2}\right]}^{1/2}$$(32) ). To adjust the damage position, it was made a sweep of the damage position along the sample. The position analysed is in agreement with the position of the accelerometers showed in Figures 7–12. As the greater levels of energy are near the first accelerometer (damage in position 1) and second accelerometer (damage in position 2), it was considered in the mathematical model (FEM model) the damage position by varying from 0 to 0.275 m with increment of 1 mm (the length of the damage). The position 0.275 m (position in the middle of the second and third accelerometers) was used for identification of the damage in position 2. To identify the damage in position 1, it was necessary to vary the position of the damage from 0 to 0.175 m. In all cases, it was considered damage with length of 1 mm and the height varying from 1 to 5 mm with increment of 1 mm. Figures 14–16 show the variation of the first, second and third natural frequencies.

Figure 14. First natural frequency (clamped-free beam).

Figure 15. Second natural frequency (clamped-free beam).

Figure 16. Third natural frequency (clamped-free beam).

The numerical and experimental data are approximated using the PSO method, and this way, it is possible to adjust the localization and the severity of the damage. The first three natural frequencies were considered, and the following parameter error checking was used:

(32) $${\mathit{\epsilon }}_{\mathit{\omega }}={\left[\frac{{\sum }_{i=1}^L{\left(\mathrm{\Delta }{\mathit{\omega }}_{A_i}-\mathrm{\Delta }{\mathit{\omega }}_{X_i}\right)}^2}{{\sum }_{i=1}^L\mathrm{\Delta }{\mathit{\omega }}_{A_i}^2}\right]}^{1/2}$$(32)

where $\mathrm{\Delta }{\mathit{\omega }}_{A_i}$ and $\mathrm{\Delta }{\mathit{\omega }}_{X_i}$ represent the variation of the experimental and numeric frequencies. Three input variables were used in the optimization process: position, length and height of the damage.

Figures 17 and 18 show the values of the error parameter defined in Equation (32(32) $${\mathit{\epsilon }}_{\mathit{\omega }}={\left[\frac{{\sum }_{i=1}^L{\left(\mathrm{\Delta }{\mathit{\omega }}_{A_i}-\mathrm{\Delta }{\mathit{\omega }}_{X_i}\right)}^2}{{\sum }_{i=1}^L\mathrm{\Delta }{\mathit{\omega }}_{A_i}^2}\right]}^{1/2}$$(32) ) for two different samples with different damage location. It can be noticed in Figure 17, that the region of minimal values (intersection of curves) for the three natural frequencies is near 75 mm (near position 1 of damage), and in Figure 18, the region is near 250 mm (near position 2 of damage).

Figure 17. Variation of the error parameter (clamped-free beam).

Figure 18. Variation of the error parameter (clamped-free beam).

The optimised parameters for the three systems are shown in Tables 4 and 5, and the real parameters (reference values) present good agreement. The PSO parameters used are as follows: particles number: 25; velocity 2; acceleration constant: 2.1; inertia weights: 0.9 and 0.6 and number of iterations: 400. The convergence demands great computational time because a finite element mesh with quadratic elements of 16 nodes and 1 mm of size was used. The second sample with damage located at 250 mm from the clamped end presented similar results.

Table 4. Optimised damage parameters (damage in position 1).

	Position (m)	Length (m)	Height (m)
Damage 1	0.0745	0.0010	0.00095
Damage 2	0.0751	0.0011	0.00249
Damage 3	0.0750	0.0011	0.00457
Reference values	0.0750	0.0010	0.00100
			0.00250
			0.00450

Table 5. Optimised damage parameters (damage in position 2).

	Position (m)	Length (m)	Height (m)
Damage 1	0.2435	0.0009	0.00096
Damage 2	0.2455	0.0012	0.00252
Damage 3	0.2460	0.0012	0.00462
Reference values	0.2450	0.0010	0.00100
			0.00250
			0.00450

Other analyses were carried out using a steel free–free beam whose dimensions are as follows: length = 0.84 m; thickness = 0.0127 m and width = 0.0254 m. In this case, was used only one accelerometer placed in position 35 cm and the excitation force was applied using an impact hammer in 16 different positions (5–80 cm with step of 5 cm). Two structural modifications were introduced in positions near 0.08 m (length of 0.001 m and height of 0.002 m) and 0.38 m (length of 0.001 m and height of 0.004 m). Figure 19 shows the energy parameter values obtained with bispectrum theory and wavelet transform and the correlation coefficient (the values of this parameter are multiplied by 3.5 to better visualisation). It can be seen two peak values in energy curves and two regions of minimum values in the range correlation curve, near to the positions 2 and 8 of the excitation force. These positions are close to the damage positions. To verify the correct positions and geometric data of the damages, it was used the error parameter defined by (32) and the PSO. As Figure 19 shows two different peaks of energy the geometric constraints imposed for the solution of the optimization problems were as follows: number of damage = 2; positions of the damages: 5 < x₁ <15 cm and 25 < x₂ <35 cm; length of the damages: 0 < length < 2 mm and height of the damages: 0 < height < 10 mm. The adjusted parameters values were as follows: position of damage 1 x₁ = 8 mm, length of 0.975 mm and the height of 2.05 mm; position of damage 2 x₂ = 0.278 m, length of 1.005 mm and the height of 3.95 mm. The values are close to the reference values.

Figure 19. Energy index (free–free beam).

In order to test the efficiency of the method in a more complex system, the parameters of a sandwich beam with damage were adjusted. The sandwich beam experimental sample made with hot-rolled steel associated with polyurethane rigid foam and high impact polystyrene is shown in Figure 20. The steel thickness is 0.6 mm; the polyurethane is 38.25 mm and the Polystyrene is 1.25 mm and the beam width is 39.18 mm. The experimental data are obtained using the impact hammer and four accelerometers displaced along the sample (A1–A4). The damage was introduced into the system by removing 2 mm from the polystyrene face at the position of 75 mm from the rigid base.

Figure 20. Sandwich beam.

Figure 21 shows the energy parameter values obtained with bispectrum theory and wavelet transform and the correlation coefficient. It can be noticed that the great energy values and the lower values of the correlation coefficient are in the first accelerometer (the accelerometer nearest the damage). Using the error parameter defined in Equation (32(32) $${\mathit{\epsilon }}_{\mathit{\omega }}={\left[\frac{{\sum }_{i=1}^L{\left(\mathrm{\Delta }{\mathit{\omega }}_{A_i}-\mathrm{\Delta }{\mathit{\omega }}_{X_i}\right)}^2}{{\sum }_{i=1}^L\mathrm{\Delta }{\mathit{\omega }}_{A_i}^2}\right]}^{1/2}$$(32) ) and the PSO whose constraints were the damage parameters for the finite elements of the meshes (Figure 2) of the polyurethane rigid foam or high impact polystyrene and the position of the damage. The geometric constraints imposed for the solution of the optimization problems were as follows: number of damage =1, positions of the damages: 0 < x < 0.225 m, length of the damages: 0 < length < 4 mm and height of the damages: 0 < height < 2 mm. The adjusted parameters values were as follows: position = 75.5 mm (from the rigid base), length of 1.975 mm (polystyrene face) and the height of 1.25 mm. The values are close to the reference values.

Figure 21. Energy index (sandwich beam).

4. Conclusion

It was noticed that the application of statistical methods for assessing damage to structures can be a potential tool. The specific application for cantilever beams presents difficulties due to variations in the vibration modes of the system. Large variations in values for the energy index using bispectrum and wavelet transform were observed, but, in all cases, the two parameters presented a good estimate of the damage position. The results obtained using wavelet transform became evident. This is because this technique differs from conventional power spectrum using parameters that vary increment of frequency in the analysis. Furthermore, the bispectrum is highly recommended for nonlinear systems.

Another parameter, correlation coefficient, using MM and taking into account the correlation of two signals, was used. The results presented minor fluctuations converging to the expected, i.e. lower values for the accelerometer near the damage. Morphological analysis appears to be a potential tool for analysing damage to structures. One of the difficulties lies in the choice of the structuring element.

The correct damage position in a clamped-free steel beam was obtained considering a performance index minimization procedure based on the correlation of numerical and experimental variations of the natural frequencies. To adjust the position of the damage, a sweep of the damage position was made along the sample using the PSO method. The numeric damage position and the optimal geometric damage parameters are close to the real values.

The results are also obtained for a free–free steel beam with two different damages. The damages parameters and positions were adjusted using only the energy parameters to find approximate location of damages and the PSO method to approximate the natural frequencies obtained experimentally and numerically.

The same procedure was applied to sandwich beam.

Disclosure statement

No potential conflict of interest was reported by the authors.