# Damage Identification of Structures Based on Smooth Orthogonal Decomposition and Improved Beetle Antennae Search Algorithm

Authors: Zhixiang Hu; Peiguan Zhang
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

A novel damage identification method that utilizes the smooth orthogonal decomposition (SOD) combined with the improved beetle antennae search algorithm (BAS) presented by previous scholars is proposed. Firstly, the damage index which can track the curvature changing of mode shape identified by the SOD method is generated by an adaptive polynomial fit method. The locations of structure damages are determined according to the damage index. Thus, the number of possible damaged elements needed to be taken into account can be reduced when calculating the degree of damage. Then, the reduction in the stiffness at the damage location of the structure is calculated by the improved BAS in which the fitness function is constructed by calculated frequencies of the damaged structure in each iteration and the modal frequencies obtained by SOD. The BAS algorithm is improved through a fusion strategy of simulated annealing theory. Thus, the improved BAS algorithm is efficient and adaptive. The effect of this combined application in damage identification has been verified by numerical examples of a simply supported beam with single damage and a cantilever beam with double damage. The numerical results show that this combined algorithm exhibits high reliability in damage identification of beam-like structures.


---

## Body

## 1. Introduction

In recent years, damage identification based on modal parameter identification combined with intelligent optimization algorithms has attracted extensive attention and have been widely studied [1, 2]. And a key problem of structural damage identification is to search the characteristic parameters closely related to structural dynamic characteristics and sensitive to damage [3–5].Scholars have studied the inverse problem in beam-like structure damage identification. A common problem that still exists is the balance between accuracy and computational costs. Khatir et al. developed a two-stage approach to locate the potential of damaged elements followed by estimating damage and they verified that the proposed approach can be used to determine correctly the severity of damage in beam structures [6].Many structural characteristic deflection shapes (CDSs), such as mode shape, contain spatial information of a structure. Thus, these common modal parameters have great significance to the damage detection and localization [7, 8]. The curvature curves of CDSs can be obtained using the second-order central difference method. Local singularities caused by damage on curvature curves can be detected as damage location of structures. A damage identification method based on CDSs and adaptive gapped smoothing method (GSM) has been studied in [7]. However, this GSM-based damage identification method could not obtain the damage degree directly.To obtain the damage degree, scholars have investigated the application of many intelligent algorithms in the field of structural damage identification [9, 10]. The beetle antennae search algorithm (BAS) [11–13] is a novel biology-inspired intelligent optimization algorithm proposed by Jiang, which arises from the simulation of foraging behavior of beetle. BAS requires few parameters and could be implemented easily. It possesses better optimization ability in low-dimensional space problems; for example, the result is not unique or convergence to local extremum in multidimensional space may exist [14]. To avoid falling into local extremum in the early stages of the algorithm and improve the optimization ability of the algorithm in multidimensional space, an improved beetle antennae search algorithm (IBAS) that fuses simulated annealing has been proposed [15]. In the present work, the IBAS is firstly introduced to be applied in the field of structural damage identification. The identification of damage degree is essentially an optimization problem to determine the stiffness of all structural elements. However, in practical engineering, taking the stiffness of all elements as an unknown quantity will generate too many parameters to be processed, resulting in an oversized calculation scale of the optimization problem. Furthermore, the damage may be located at the symmetrical position of the actual damage location by the algorithm.In this paper, we proposed a novel damage identification method of beam structure based on SOD and IBAS algorithm. The SOD is used to obtain the modal frequencies and mode shapes of structures from the dynamic response excited by stochastic external forces, which is a new modal identification method that has a simple structure [16, 17]. In the field of vibration engineering, the SOD method is originally applied in multi-degree-of-freedom undamped free vibration systems [17]. Then, the SOD has been developed to identify modal parameters of systems excited by stochastic forces [18]. It can efficiently compute modal parameters from output-only responses of vibrating structures. The SOD method does not require uniform structure mass compared with POD method, which can directly identify structure frequency [16]. Specifically, after using SOD method to obtain the mode shape, the modal curvature is calculated by the second-order central difference method. Inspired by the adaptive GSM method, this paper employs the least square method to fit the modal curvatures of each mode, and the optimal order of the fitting polynomial is determined through cross-validation. Thus, the polynomial has a certain adaptive ability when fitting the modal curvatures. By comparing the modal curvature curves and fit curves, the damage index is constructed to locate damage through determining the damaged elements. Then, the IBAS algorithm and the structural frequencies identified by SOD method are used to calculate the damage degree of the structure. To improve the speed of convergence and the accuracy of the algorithm when calculating the degree of structural damage, the damage location results obtained by considering the damage index are utilized. The reduction in the stiffness of the identified damaged elements is termed as the unknown quantity in the subsequent IBAS algorithm. This operation significantly reduces the dimension of the optimization problem’s solution.The damage identification method of beam structures based on SOD and IBAS algorithm has been verified using numerical examples of beams under two damage cases in this paper. The results show that the SOD method can accurately identify the structural frequencies and mode shapes. The damaged elements can be accurately determined by the damage localization method using the results of SOD. The reduction in the stiffness of the damaged elements can be calculated accurately according to the structural frequencies identified by SOD method and IBAS algorithm.The remainder of this article is organized as follows. The damage localization method based on SOD is presented in Section2. The damage identification combined with IBAS algorithm is presented in Section 3. In Section 4, two numerical examples are carried out to verify the performance of the method. The conclusion is drawn in Section 5.


## 2. Damage Localization Method Based on SOD

### 2.1. Principle of SOD Method

Let matrixX denote an n×N ensemble of displacements obtained from N time samples at n lumped masses. The total sampling time T should contain several cycles of the first-mode vibration. Similarly, let V=X˙ and A=X¨ denote the ensembles of velocities and accelerations ensemble data, respectively. Next, the displacement covariance matrix Rxx=XXT/N−1 and velocity covariance matrix Rvv=VVT/N−1 are formed. Then, the smooth orthogonal decomposition can retrieve the modal parameters from the following generalized eigenvalue problem:(1)λRxxψ=Rvvψ,where the generalized eigenvalues λ are the smooth orthogonal values (SOVs), which approach the frequencies squared as N increases. The generalized eigenvectors ψ are termed as smooth orthogonal modes (SOMs) and the matrix Ψ=ψ1,ψ2,…,ψn can be formed, the inverse transpose of which approaches the linear modal matrix. The matrix form of equation (1) is(2)RxxΨΔ=RvvΨ,here, Δ=diagλ1,λ2,…,λn. Each diagonal element of Δ is an eigenvalue. It has been proved that the eigenvalue problem of SOD can reduce to the eigenvalue problem of the undamped vibration system [16]. For the undamped and slightly damped vibration systems, the modal matrix of the vibration system can be estimated by relationship Φ=Ψ−T, and the natural frequencies can be estimated by SOVs, ωi2=λi. In applications, the eigenvectors are sorted by the eigenvalues.To briefly present the concept of the SOD method and study the theory of modal identification, according to the modal superposition property, we have(3)X=ΦQ,V=ΦQ˙,where Q denotes the ensemble data of modal coordinate responses. For the undamped free vibration system, each row of the matrix Q is a discretely sampled sine waving with the natural frequency of the system. Therefore, Rqq=QQT/N−1 and Rq˙q˙=Q˙Q˙T/N−1 are diagonal matrices. It is easy to obtain that(4)ΨTRxxΨ=Rqq,ΨTRvvΨ=Rq˙q˙.Substituting this equation in equation (2), it can be deduced that the SOD directly produces estimates of the modal frequencies squared, i.e., Δ=Rqq−1Rq˙q˙. In other words, the SOM matrix Ψ can simultaneously diagonalize the matrices Rxx and Rvv.For detection and identification of the damage of the system, the modal identification results of the SOD are utilized. In practice, the lower modes are always excited and the higher modes are often not active. The number of available modesn1 needs to be selected [7]. There are a series of methods establishing n1, and two common methods are introduced here. The first method is to search the minimum integer n1 that satisfies(5)∑i=1n1λ^i∑i=1nλ^i≥be,where λ^i is the ith diagonal element of the matrix Rqq in which diagonal elements have been taken absolute values and sorted in descending order; be denotes the participation threshold value. Another method is to search λ^n1 that satisfies λ^n1≫λ^n1+1. For the latter method, the eigenvalue plot is commonly used to determine n1. The numerical examples in this paper adopt the first method to determine n1; be=0.96 is used. The eigenvector matrix Φ identified by SOD method is used to calculate the matrix Rqq=ΨTRxxΨ=Φ−1RxxΦ−T. The column vectors of Φ corresponding to the first n1 modes are taken as the mode shapes. The corresponding diagonal elements in Δ are used to calculate structural frequencies.


### 2.2. Damage Localization Index

In this paper, damage localization of structures is realized by curvature fitting and constructing damage index. To be specific, firstly, the modal curvatures are obtained by using the second-order central difference method to process the mode shapes identified by SOD method. Then, the smoothed modal curvatures are obtained by using the adaptive polynomial fit method to fit the modal curvatures. To localize the damage, a damage index is introduced and computed in the method.


#### 2.2.1. Modal Curvature

It has been introduced in [19] that modal curvature φ″ can be approximately expressed as(6)φ″=d2φdz2,where φ is displacement mode shape function and z denotes location coordinate.When the measurement points are evenly spaced, each element in the modal curvature matrix can be expressed as(7)φk,l″=φk,l−1−2φk,l+φk,l+1h2,where φk,l denotes the modal displacement at the lth measurement point of the kth mode and h denotes the distance between adjacent measurement points [20].The modal curvature values at the first and last measurement points cannot be directly calculated by equation (7). They can be estimated as follows (assuming there are a total of m measurement points):(8)φk,1″=2φk,2″−φk,3″,φk,m″=2φk,m−1″−φk,m−2″.The derivation and verification of these curvature computation formulas are described in detail in [19]. In this paper, φk″ denotes the kth modal curvature. It is shown that there is a significant difference between the modal curvature values and the smoothed modal curvature values at the measurement points near the damage [7]. Thus, it is crucial to estimate the smoothed modal curvature φ^k″ by using adaptive polynomial fit method and calculate the shape distortion Δφ^k″.


#### 2.2.2. Adaptive Polynomial Fit Method

To obtain smoothed CDS curvature such as smoothed modal curvature, Ratcliffe proposed a gapped smoothing method (GSM) [21], which is expressed as(9)φk″z,c=c0+c1z+c2z2+c3z3,where z denotes the location coordinate of the modal curvature and c=c0,c1,c2,c3 are the coefficients of the cubic polynomial. The two CDS curvature values φk,l−2″ and φk,l−1″ before measurement point l and the two CDS curvature φk,l+1″ and φk,l+2″ after measurement point l are used to calculate c and predict φk″ at measurement point l.The smoothed modal curvatures are calculated by using the least square method to conduct polynomial fitting of the modal curvatures in the presented work. To enhance the noise robustness and improve the accuracy of damage identification, an adaptive gapped smoothing method [7] is used in the present work:(10)φk″z,c,r=c0+c1z+c2z2+⋯+crzr,where c=c0,c1,c2,…,cr contains the polynomial coefficients and r is the adaptive polynomial order.Inspired by the adaptive gapped smoothing method, in order to obtain the optimal polynomial orderr for polynomial fit, K-fold cross-validation is adopted to determine r in the range of 0,rmax. If the number of measurement points m is large, it is better to separate measurement points into several segments and perform polynomial fit, respectively.The basic idea of K-fold cross-validation is to divide all measurement points intokc subsamples. One of the subsamples is taken as test data and the other subsamples are used for training. Each subsample should be selected as test data once, i.e., the K-fold cross-validation includes kc trainings. In this work, the number of subsets in each K-fold cross-validation is set to 3, i.e., kc=3. And the maximum polynomial order rmax=kc−1/kcm−1−1 (where ∗ is the floor operator). The optimal polynomial order r is calculated through rmax cycles. The K-fold cross-validation is performed at each cycle, and the number of the current cycle i is taken as the order of the fit polynomial in K-fold cross-validation. During each training in the K-fold cross-validation, the least square method is employed to solve the fit curve of the training subsamples. After obtaining the fit curve, the fit curve values at the coordinates of the test subsamples are calculated. Then, the differences between the modal curvature curves and fit curves at the coordinate of the test subsamples are calculated. After calculating the absolute values of the differences at all test subsamples, the sum of above absolute values in each training is calculated. And the mean value of the sums in kc trainings at each cycle is calculated later. The optimal polynomial order r is the order of fit polynomial with the minimum mean value in rmax cycles. Therefore, in this study, the polynomial has a certain adaptive ability when fitting the modal curvatures.


#### 2.2.3. Damage Index

The damage-induced shape distortions can be obtained using the polynomial smoothing method to process curvature of CDSs. Each CDS (such as mode shape of the first mode) is sensitive to damage at some locations while being less sensitive to damage at other locations [7]. Hence, a damage identification index including damage features of several mode shapes should be more robust and accurate. To reach the above objectives, this study refers to the damage identification index proposed in [7]. The distance between the original modal curvature Φd″ and smoothed modal curvature Φ^″ obtained by adaptive polynomial fit method at each measurement point is defined as the damage index (DI):(11)DIl=ρlγlδl,ρl=∏k=1Kmaxφ^k″−φkd″1/K,γlΦ^″Φd″=∑k=1Kwkφ^k,l″−φk,ld″,δlΦ^″Φd″=∑k=1Kwk2φ^k,l″−φk,ld″2,where n1 is obtained by equation (5); φk,ld″ denotes the value at the lth measurement point of the kth modal curvature calculated by second-order central difference method; φ^k,l″ denotes the value at the lth measurement point of the kth smoothed modal curvature obtained by using the adaptive polynomial fit method to process modal curvature; l denotes measurement point in the damaged structure; KK≤n1 denotes the number of the mode shapes used for damage identification; ρl provides a global magnitude measure of the damage-induced shape changes over several mode shapes; and the weighing coefficient wk=1/maxΔφk″ and wherein Δφk″=φ^k″−φkd″. Furthermore, γl and δl are the normalized first-order and second-order moments of shape distortions Δφ^k,l″, respectively, Δφk,l″=φ^k,l″−φk,ld″.Although it is relatively simple to calculate the value of modal curvatureφk,ld″ by using the second-order central difference method, the accuracy of this method is easily affected by the density of measurement points and noise level. In general, dense sensor deployment and accurate extraction of mode shapes are necessary to ensure good accuracy of mode curvature estimation. Nevertheless, the damaged elements are determined by surveying the mutation of damage index given by equation (11). Then, the results of damage localization and structural frequencies obtained by SOD method are utilized in IBAS algorithm to calculate the damage degrees of the structures.


## 2.1. Principle of SOD Method

Let matrixX denote an n×N ensemble of displacements obtained from N time samples at n lumped masses. The total sampling time T should contain several cycles of the first-mode vibration. Similarly, let V=X˙ and A=X¨ denote the ensembles of velocities and accelerations ensemble data, respectively. Next, the displacement covariance matrix Rxx=XXT/N−1 and velocity covariance matrix Rvv=VVT/N−1 are formed. Then, the smooth orthogonal decomposition can retrieve the modal parameters from the following generalized eigenvalue problem:(1)λRxxψ=Rvvψ,where the generalized eigenvalues λ are the smooth orthogonal values (SOVs), which approach the frequencies squared as N increases. The generalized eigenvectors ψ are termed as smooth orthogonal modes (SOMs) and the matrix Ψ=ψ1,ψ2,…,ψn can be formed, the inverse transpose of which approaches the linear modal matrix. The matrix form of equation (1) is(2)RxxΨΔ=RvvΨ,here, Δ=diagλ1,λ2,…,λn. Each diagonal element of Δ is an eigenvalue. It has been proved that the eigenvalue problem of SOD can reduce to the eigenvalue problem of the undamped vibration system [16]. For the undamped and slightly damped vibration systems, the modal matrix of the vibration system can be estimated by relationship Φ=Ψ−T, and the natural frequencies can be estimated by SOVs, ωi2=λi. In applications, the eigenvectors are sorted by the eigenvalues.To briefly present the concept of the SOD method and study the theory of modal identification, according to the modal superposition property, we have(3)X=ΦQ,V=ΦQ˙,where Q denotes the ensemble data of modal coordinate responses. For the undamped free vibration system, each row of the matrix Q is a discretely sampled sine waving with the natural frequency of the system. Therefore, Rqq=QQT/N−1 and Rq˙q˙=Q˙Q˙T/N−1 are diagonal matrices. It is easy to obtain that(4)ΨTRxxΨ=Rqq,ΨTRvvΨ=Rq˙q˙.Substituting this equation in equation (2), it can be deduced that the SOD directly produces estimates of the modal frequencies squared, i.e., Δ=Rqq−1Rq˙q˙. In other words, the SOM matrix Ψ can simultaneously diagonalize the matrices Rxx and Rvv.For detection and identification of the damage of the system, the modal identification results of the SOD are utilized. In practice, the lower modes are always excited and the higher modes are often not active. The number of available modesn1 needs to be selected [7]. There are a series of methods establishing n1, and two common methods are introduced here. The first method is to search the minimum integer n1 that satisfies(5)∑i=1n1λ^i∑i=1nλ^i≥be,where λ^i is the ith diagonal element of the matrix Rqq in which diagonal elements have been taken absolute values and sorted in descending order; be denotes the participation threshold value. Another method is to search λ^n1 that satisfies λ^n1≫λ^n1+1. For the latter method, the eigenvalue plot is commonly used to determine n1. The numerical examples in this paper adopt the first method to determine n1; be=0.96 is used. The eigenvector matrix Φ identified by SOD method is used to calculate the matrix Rqq=ΨTRxxΨ=Φ−1RxxΦ−T. The column vectors of Φ corresponding to the first n1 modes are taken as the mode shapes. The corresponding diagonal elements in Δ are used to calculate structural frequencies.


## 2.2. Damage Localization Index

In this paper, damage localization of structures is realized by curvature fitting and constructing damage index. To be specific, firstly, the modal curvatures are obtained by using the second-order central difference method to process the mode shapes identified by SOD method. Then, the smoothed modal curvatures are obtained by using the adaptive polynomial fit method to fit the modal curvatures. To localize the damage, a damage index is introduced and computed in the method.


### 2.2.1. Modal Curvature

It has been introduced in [19] that modal curvature φ″ can be approximately expressed as(6)φ″=d2φdz2,where φ is displacement mode shape function and z denotes location coordinate.When the measurement points are evenly spaced, each element in the modal curvature matrix can be expressed as(7)φk,l″=φk,l−1−2φk,l+φk,l+1h2,where φk,l denotes the modal displacement at the lth measurement point of the kth mode and h denotes the distance between adjacent measurement points [20].The modal curvature values at the first and last measurement points cannot be directly calculated by equation (7). They can be estimated as follows (assuming there are a total of m measurement points):(8)φk,1″=2φk,2″−φk,3″,φk,m″=2φk,m−1″−φk,m−2″.The derivation and verification of these curvature computation formulas are described in detail in [19]. In this paper, φk″ denotes the kth modal curvature. It is shown that there is a significant difference between the modal curvature values and the smoothed modal curvature values at the measurement points near the damage [7]. Thus, it is crucial to estimate the smoothed modal curvature φ^k″ by using adaptive polynomial fit method and calculate the shape distortion Δφ^k″.


### 2.2.2. Adaptive Polynomial Fit Method

To obtain smoothed CDS curvature such as smoothed modal curvature, Ratcliffe proposed a gapped smoothing method (GSM) [21], which is expressed as(9)φk″z,c=c0+c1z+c2z2+c3z3,where z denotes the location coordinate of the modal curvature and c=c0,c1,c2,c3 are the coefficients of the cubic polynomial. The two CDS curvature values φk,l−2″ and φk,l−1″ before measurement point l and the two CDS curvature φk,l+1″ and φk,l+2″ after measurement point l are used to calculate c and predict φk″ at measurement point l.The smoothed modal curvatures are calculated by using the least square method to conduct polynomial fitting of the modal curvatures in the presented work. To enhance the noise robustness and improve the accuracy of damage identification, an adaptive gapped smoothing method [7] is used in the present work:(10)φk″z,c,r=c0+c1z+c2z2+⋯+crzr,where c=c0,c1,c2,…,cr contains the polynomial coefficients and r is the adaptive polynomial order.Inspired by the adaptive gapped smoothing method, in order to obtain the optimal polynomial orderr for polynomial fit, K-fold cross-validation is adopted to determine r in the range of 0,rmax. If the number of measurement points m is large, it is better to separate measurement points into several segments and perform polynomial fit, respectively.The basic idea of K-fold cross-validation is to divide all measurement points intokc subsamples. One of the subsamples is taken as test data and the other subsamples are used for training. Each subsample should be selected as test data once, i.e., the K-fold cross-validation includes kc trainings. In this work, the number of subsets in each K-fold cross-validation is set to 3, i.e., kc=3. And the maximum polynomial order rmax=kc−1/kcm−1−1 (where ∗ is the floor operator). The optimal polynomial order r is calculated through rmax cycles. The K-fold cross-validation is performed at each cycle, and the number of the current cycle i is taken as the order of the fit polynomial in K-fold cross-validation. During each training in the K-fold cross-validation, the least square method is employed to solve the fit curve of the training subsamples. After obtaining the fit curve, the fit curve values at the coordinates of the test subsamples are calculated. Then, the differences between the modal curvature curves and fit curves at the coordinate of the test subsamples are calculated. After calculating the absolute values of the differences at all test subsamples, the sum of above absolute values in each training is calculated. And the mean value of the sums in kc trainings at each cycle is calculated later. The optimal polynomial order r is the order of fit polynomial with the minimum mean value in rmax cycles. Therefore, in this study, the polynomial has a certain adaptive ability when fitting the modal curvatures.


### 2.2.3. Damage Index

The damage-induced shape distortions can be obtained using the polynomial smoothing method to process curvature of CDSs. Each CDS (such as mode shape of the first mode) is sensitive to damage at some locations while being less sensitive to damage at other locations [7]. Hence, a damage identification index including damage features of several mode shapes should be more robust and accurate. To reach the above objectives, this study refers to the damage identification index proposed in [7]. The distance between the original modal curvature Φd″ and smoothed modal curvature Φ^″ obtained by adaptive polynomial fit method at each measurement point is defined as the damage index (DI):(11)DIl=ρlγlδl,ρl=∏k=1Kmaxφ^k″−φkd″1/K,γlΦ^″Φd″=∑k=1Kwkφ^k,l″−φk,ld″,δlΦ^″Φd″=∑k=1Kwk2φ^k,l″−φk,ld″2,where n1 is obtained by equation (5); φk,ld″ denotes the value at the lth measurement point of the kth modal curvature calculated by second-order central difference method; φ^k,l″ denotes the value at the lth measurement point of the kth smoothed modal curvature obtained by using the adaptive polynomial fit method to process modal curvature; l denotes measurement point in the damaged structure; KK≤n1 denotes the number of the mode shapes used for damage identification; ρl provides a global magnitude measure of the damage-induced shape changes over several mode shapes; and the weighing coefficient wk=1/maxΔφk″ and wherein Δφk″=φ^k″−φkd″. Furthermore, γl and δl are the normalized first-order and second-order moments of shape distortions Δφ^k,l″, respectively, Δφk,l″=φ^k,l″−φk,ld″.Although it is relatively simple to calculate the value of modal curvatureφk,ld″ by using the second-order central difference method, the accuracy of this method is easily affected by the density of measurement points and noise level. In general, dense sensor deployment and accurate extraction of mode shapes are necessary to ensure good accuracy of mode curvature estimation. Nevertheless, the damaged elements are determined by surveying the mutation of damage index given by equation (11). Then, the results of damage localization and structural frequencies obtained by SOD method are utilized in IBAS algorithm to calculate the damage degrees of the structures.


## 2.2.1. Modal Curvature

It has been introduced in [19] that modal curvature φ″ can be approximately expressed as(6)φ″=d2φdz2,where φ is displacement mode shape function and z denotes location coordinate.When the measurement points are evenly spaced, each element in the modal curvature matrix can be expressed as(7)φk,l″=φk,l−1−2φk,l+φk,l+1h2,where φk,l denotes the modal displacement at the lth measurement point of the kth mode and h denotes the distance between adjacent measurement points [20].The modal curvature values at the first and last measurement points cannot be directly calculated by equation (7). They can be estimated as follows (assuming there are a total of m measurement points):(8)φk,1″=2φk,2″−φk,3″,φk,m″=2φk,m−1″−φk,m−2″.The derivation and verification of these curvature computation formulas are described in detail in [19]. In this paper, φk″ denotes the kth modal curvature. It is shown that there is a significant difference between the modal curvature values and the smoothed modal curvature values at the measurement points near the damage [7]. Thus, it is crucial to estimate the smoothed modal curvature φ^k″ by using adaptive polynomial fit method and calculate the shape distortion Δφ^k″.


## 2.2.2. Adaptive Polynomial Fit Method

To obtain smoothed CDS curvature such as smoothed modal curvature, Ratcliffe proposed a gapped smoothing method (GSM) [21], which is expressed as(9)φk″z,c=c0+c1z+c2z2+c3z3,where z denotes the location coordinate of the modal curvature and c=c0,c1,c2,c3 are the coefficients of the cubic polynomial. The two CDS curvature values φk,l−2″ and φk,l−1″ before measurement point l and the two CDS curvature φk,l+1″ and φk,l+2″ after measurement point l are used to calculate c and predict φk″ at measurement point l.The smoothed modal curvatures are calculated by using the least square method to conduct polynomial fitting of the modal curvatures in the presented work. To enhance the noise robustness and improve the accuracy of damage identification, an adaptive gapped smoothing method [7] is used in the present work:(10)φk″z,c,r=c0+c1z+c2z2+⋯+crzr,where c=c0,c1,c2,…,cr contains the polynomial coefficients and r is the adaptive polynomial order.Inspired by the adaptive gapped smoothing method, in order to obtain the optimal polynomial orderr for polynomial fit, K-fold cross-validation is adopted to determine r in the range of 0,rmax. If the number of measurement points m is large, it is better to separate measurement points into several segments and perform polynomial fit, respectively.The basic idea of K-fold cross-validation is to divide all measurement points intokc subsamples. One of the subsamples is taken as test data and the other subsamples are used for training. Each subsample should be selected as test data once, i.e., the K-fold cross-validation includes kc trainings. In this work, the number of subsets in each K-fold cross-validation is set to 3, i.e., kc=3. And the maximum polynomial order rmax=kc−1/kcm−1−1 (where ∗ is the floor operator). The optimal polynomial order r is calculated through rmax cycles. The K-fold cross-validation is performed at each cycle, and the number of the current cycle i is taken as the order of the fit polynomial in K-fold cross-validation. During each training in the K-fold cross-validation, the least square method is employed to solve the fit curve of the training subsamples. After obtaining the fit curve, the fit curve values at the coordinates of the test subsamples are calculated. Then, the differences between the modal curvature curves and fit curves at the coordinate of the test subsamples are calculated. After calculating the absolute values of the differences at all test subsamples, the sum of above absolute values in each training is calculated. And the mean value of the sums in kc trainings at each cycle is calculated later. The optimal polynomial order r is the order of fit polynomial with the minimum mean value in rmax cycles. Therefore, in this study, the polynomial has a certain adaptive ability when fitting the modal curvatures.


## 2.2.3. Damage Index

The damage-induced shape distortions can be obtained using the polynomial smoothing method to process curvature of CDSs. Each CDS (such as mode shape of the first mode) is sensitive to damage at some locations while being less sensitive to damage at other locations [7]. Hence, a damage identification index including damage features of several mode shapes should be more robust and accurate. To reach the above objectives, this study refers to the damage identification index proposed in [7]. The distance between the original modal curvature Φd″ and smoothed modal curvature Φ^″ obtained by adaptive polynomial fit method at each measurement point is defined as the damage index (DI):(11)DIl=ρlγlδl,ρl=∏k=1Kmaxφ^k″−φkd″1/K,γlΦ^″Φd″=∑k=1Kwkφ^k,l″−φk,ld″,δlΦ^″Φd″=∑k=1Kwk2φ^k,l″−φk,ld″2,where n1 is obtained by equation (5); φk,ld″ denotes the value at the lth measurement point of the kth modal curvature calculated by second-order central difference method; φ^k,l″ denotes the value at the lth measurement point of the kth smoothed modal curvature obtained by using the adaptive polynomial fit method to process modal curvature; l denotes measurement point in the damaged structure; KK≤n1 denotes the number of the mode shapes used for damage identification; ρl provides a global magnitude measure of the damage-induced shape changes over several mode shapes; and the weighing coefficient wk=1/maxΔφk″ and wherein Δφk″=φ^k″−φkd″. Furthermore, γl and δl are the normalized first-order and second-order moments of shape distortions Δφ^k,l″, respectively, Δφk,l″=φ^k,l″−φk,ld″.Although it is relatively simple to calculate the value of modal curvatureφk,ld″ by using the second-order central difference method, the accuracy of this method is easily affected by the density of measurement points and noise level. In general, dense sensor deployment and accurate extraction of mode shapes are necessary to ensure good accuracy of mode curvature estimation. Nevertheless, the damaged elements are determined by surveying the mutation of damage index given by equation (11). Then, the results of damage localization and structural frequencies obtained by SOD method are utilized in IBAS algorithm to calculate the damage degrees of the structures.


## 3. Damage Identification with IBAS Algorithm

### 3.1. IBAS Algorithm

The beetle antennae search (BAS) algorithm is a novel biology-inspired intelligent optimization algorithm arising from the simulation of the foraging behavior of beetle [11–13]. When a beetle is foraging, its left and right antennae are used to sense the odor intensity of the food. If the left antennae receive stronger odor intensity compared with the right antennae, the beetle will fly to left in the next step and vice versa. In BAS algorithm, the location of beetle is the feasible solution for the problem; the odor intensity of food is the fitness function. The modeling process of basic BAS algorithm is shown as follows:Step 1. Initialize the position of individual beetle in the k-dimensional space.Step 2. Create the left and right beetle antennae’s space coordinates:(12)xl=xt−d⋅b,xr=xt+d⋅b,wheret is the current iteration number; d denotes the distance between the beetle’s centre of mass and beetle antennae; and b=rndk,1/rndk,1 is a random unit vector that represents the direction from the left antennae to the right antennae.Step 3. Update the coordinate of beetle:(13)xt+1=xt+St⋅b⋅signfxl−fxr,wherefx is the fitness function value at x; St denotes the step size of the tth iteration; and sign∗ is the sign function. In practical studies, St=α⋅St−1 is proposed, where α is the step size factor, which can scale down the search step size to improve the convergence speed of the algorithm [11].The step size factorα in the basic BAS algorithm is fixed. In order to enhance the search capability of the algorithm, a method with changing step size factor has been proposed in [15], i.e.,(14)β=α−0.2⋅i−15⋅N1+0.5,fi>fmin,β=α,fi≤fmin,where fi is the current fitness function value; fmin is the historical best-fitness function value; i and N1 denote current iteration number and total iteration number, respectively; α is the default step size factor, which is generally set to 0.95; and β is the current step size factor.In order to avoid falling into the local extremum in the early stage of the algorithm, the simulated annealing algorithm (SA) [15, 22] can be adopted, which is an intelligent algorithm simulating the annealing process of solid matter in physics. The SA can search the global optimal solution and has fast convergence speed. The feature of the SA is that it can jump out of partial optimization with a certain probability in the search process, which makes the algorithm effectively avoid falling into the local extremum. Specifically, the SA algorithm can also accept a bad solution with a certain probability except for accepting the good solution in the search process. The value of probability is reduced gradually as the temperature drops under the influence of the temperature parameters. This makes the algorithm have a large probability of jumping out of local extremum in the early stage while having high convergence speed in the later steps.In the SA algorithm, the Metropolis criterion is usually used as the standard for accepting the new solution [23]. The incremental ΔT=fxt+1−fxt is calculated firstly. If ΔT<0, xt+1 will be accepted as the new solution; otherwise, calculate p=exp−ΔT/T, which is the probability to accept xt+1 as the new solution, i.e., compare p with a random number between 0,1; if p is greater than the random number, xt+1 will be accepted as the new solution; otherwise, xt+1 will not be accepted, where exp∗ is the exponential function of the natural constant e; T is the current temperature.Therefore, there are two layers of loops in the IBAS algorithm. In the external cycle, the current temperatureT is used to control the initial step size S of each internal cycle in the basic BAS. In the internal cycle, the Metropolis criterion is used to update the next position of beetle. The temperature T will decrease as the iterative number increases; thus, the step size of the algorithm will also decrease. This relationship ensures the convergence of the algorithm in the later period.Refer to [24], the fitness function of the IBAS algorithm in this paper is(15)f=∑in1ωir−ωic2ωir2,where values of n1 are obtained by equation (5); ωir is the calculated damaged structural frequencies in each iteration of IBAS; ωic represents the natural frequency of the damaged structure extracted by SOD method. In the numerical simulations, the accurate damaged structural frequencies can be extracted from damaged finite element method (FEM) model with preset damage values, which can be used to verify the accuracy of modal frequencies obtained by SOD.


### 3.2. Damage Identification

The optimization ability of IBAS algorithm in multidimensional problems is proven to be better than the basic BAS algorithm [15]. However, in practical engineering, the structure is generally divided into a large number of elements. Taking the stiffness of all elements as unknown quantities in IBAS algorithm probably still causes problems such as slow convergence speed and low optimization accuracy. Based on the mode shapes obtained by using SOD, modal curvatures are calculated by the second-order central difference method. Smoothed modal curvatures are constructed by using the adaptive polynomial fit method to fit modal curvatures. By comparing the modal curvature curves with smoothed modal curvature curves, the damage index is constructed and the damaged elements can be determined. Therefore, it is feasible to reduce the dimension of the optimization problem by determining the damaged element in advance. Besides, when IBAS algorithm is used alone for damage identification of symmetric structures, the damage position could be located at the symmetrical position of the actual damage location. By prejudging the damage location, one can effectively avoid the occurrence of the above problem. To illustrate the above ideas, the application of IBAS algorithm in damage identification of beam-like structures is studied in the present work.The structural frequencies obtained by SOD method are used to solve the fitness function in IBAS algorithm. Damaged elements of the beam structure are identified by the abovementioned damage localization method firstly. The number of identified damaged elements is taken as the dimension of the solution in IBAS algorithm. The position coordinates of beetleX denote the stiffness reduction ratios of identified damaged elements.The basic process of the proposed algorithm to calculate the stiffness reduction ratios of damaged elements is shown in Figure1.Figure 1
Schematic diagram of damage identification.


## 3.1. IBAS Algorithm

The beetle antennae search (BAS) algorithm is a novel biology-inspired intelligent optimization algorithm arising from the simulation of the foraging behavior of beetle [11–13]. When a beetle is foraging, its left and right antennae are used to sense the odor intensity of the food. If the left antennae receive stronger odor intensity compared with the right antennae, the beetle will fly to left in the next step and vice versa. In BAS algorithm, the location of beetle is the feasible solution for the problem; the odor intensity of food is the fitness function. The modeling process of basic BAS algorithm is shown as follows:Step 1. Initialize the position of individual beetle in the k-dimensional space.Step 2. Create the left and right beetle antennae’s space coordinates:(12)xl=xt−d⋅b,xr=xt+d⋅b,wheret is the current iteration number; d denotes the distance between the beetle’s centre of mass and beetle antennae; and b=rndk,1/rndk,1 is a random unit vector that represents the direction from the left antennae to the right antennae.Step 3. Update the coordinate of beetle:(13)xt+1=xt+St⋅b⋅signfxl−fxr,wherefx is the fitness function value at x; St denotes the step size of the tth iteration; and sign∗ is the sign function. In practical studies, St=α⋅St−1 is proposed, where α is the step size factor, which can scale down the search step size to improve the convergence speed of the algorithm [11].The step size factorα in the basic BAS algorithm is fixed. In order to enhance the search capability of the algorithm, a method with changing step size factor has been proposed in [15], i.e.,(14)β=α−0.2⋅i−15⋅N1+0.5,fi>fmin,β=α,fi≤fmin,where fi is the current fitness function value; fmin is the historical best-fitness function value; i and N1 denote current iteration number and total iteration number, respectively; α is the default step size factor, which is generally set to 0.95; and β is the current step size factor.In order to avoid falling into the local extremum in the early stage of the algorithm, the simulated annealing algorithm (SA) [15, 22] can be adopted, which is an intelligent algorithm simulating the annealing process of solid matter in physics. The SA can search the global optimal solution and has fast convergence speed. The feature of the SA is that it can jump out of partial optimization with a certain probability in the search process, which makes the algorithm effectively avoid falling into the local extremum. Specifically, the SA algorithm can also accept a bad solution with a certain probability except for accepting the good solution in the search process. The value of probability is reduced gradually as the temperature drops under the influence of the temperature parameters. This makes the algorithm have a large probability of jumping out of local extremum in the early stage while having high convergence speed in the later steps.In the SA algorithm, the Metropolis criterion is usually used as the standard for accepting the new solution [23]. The incremental ΔT=fxt+1−fxt is calculated firstly. If ΔT<0, xt+1 will be accepted as the new solution; otherwise, calculate p=exp−ΔT/T, which is the probability to accept xt+1 as the new solution, i.e., compare p with a random number between 0,1; if p is greater than the random number, xt+1 will be accepted as the new solution; otherwise, xt+1 will not be accepted, where exp∗ is the exponential function of the natural constant e; T is the current temperature.Therefore, there are two layers of loops in the IBAS algorithm. In the external cycle, the current temperatureT is used to control the initial step size S of each internal cycle in the basic BAS. In the internal cycle, the Metropolis criterion is used to update the next position of beetle. The temperature T will decrease as the iterative number increases; thus, the step size of the algorithm will also decrease. This relationship ensures the convergence of the algorithm in the later period.Refer to [24], the fitness function of the IBAS algorithm in this paper is(15)f=∑in1ωir−ωic2ωir2,where values of n1 are obtained by equation (5); ωir is the calculated damaged structural frequencies in each iteration of IBAS; ωic represents the natural frequency of the damaged structure extracted by SOD method. In the numerical simulations, the accurate damaged structural frequencies can be extracted from damaged finite element method (FEM) model with preset damage values, which can be used to verify the accuracy of modal frequencies obtained by SOD.


## 3.2. Damage Identification

The optimization ability of IBAS algorithm in multidimensional problems is proven to be better than the basic BAS algorithm [15]. However, in practical engineering, the structure is generally divided into a large number of elements. Taking the stiffness of all elements as unknown quantities in IBAS algorithm probably still causes problems such as slow convergence speed and low optimization accuracy. Based on the mode shapes obtained by using SOD, modal curvatures are calculated by the second-order central difference method. Smoothed modal curvatures are constructed by using the adaptive polynomial fit method to fit modal curvatures. By comparing the modal curvature curves with smoothed modal curvature curves, the damage index is constructed and the damaged elements can be determined. Therefore, it is feasible to reduce the dimension of the optimization problem by determining the damaged element in advance. Besides, when IBAS algorithm is used alone for damage identification of symmetric structures, the damage position could be located at the symmetrical position of the actual damage location. By prejudging the damage location, one can effectively avoid the occurrence of the above problem. To illustrate the above ideas, the application of IBAS algorithm in damage identification of beam-like structures is studied in the present work.The structural frequencies obtained by SOD method are used to solve the fitness function in IBAS algorithm. Damaged elements of the beam structure are identified by the abovementioned damage localization method firstly. The number of identified damaged elements is taken as the dimension of the solution in IBAS algorithm. The position coordinates of beetleX denote the stiffness reduction ratios of identified damaged elements.The basic process of the proposed algorithm to calculate the stiffness reduction ratios of damaged elements is shown in Figure1.Figure 1
Schematic diagram of damage identification.


## 4. Numerical Examples

### 4.1. Calculation Model and Parameters

A simply supported beam model and a cantilever beam model both with a span of 0.75 m have been established in ANSYS software. The element type is the Beam3 Euler beam element, and the beam models are divided into 25 elements. The geometric and material parameters are given in Table1. Gaussian white noise arrays are imported to simulate random excitation. The modal analysis of the damaged structure is analyzed firstly to compute the theoretical modal parameters. And then, the undamped transient dynamic analysis is carried out using the full method.Table 1
Dimensions and material properties.
PropertiesMean valueWidth (m)0.025Thickness (m)0.006Young’s modulus (GPa)200Density (kg/m3)7800Poisson’s ratio0.3Numerical simulations of two damage cases in this study are shown in Table2. The damage is simulated by reducing the stiffness of the specific element model. In the process of damage identification, stiffness reduction is realized by reducing elasticity modulus of the damaged element as [25]:(16)Eid=1−aiE0,where Eid is the elasticity modulus of the ith damaged element, E0 denotes the pristine elasticity modulus, and ai denotes the reduction in stiffness of the ith damaged element, which is equal to the ratio between the stiffness reduction and the pristine stiffness, and it is also called damage degree in this paper.Table 2
Damage cases.
Damage caseDamaged element numberElement stiffness reductionSingle damage1030%Double damage8, 1830%, 30%The element stiffness matrix and element mass matrix of the middle element of the undamaged equal section beam model are used to assemble the global stiffness matrix and global mass matrix, respectively. And the global stiffness matrix and global mass matrix that are imported to compute the reduction in stiffness of the damaged element by iterating are used to calculate the value of damaged structural frequencyωir by finite element theory at each iteration of IBAS algorithm. The initial temperature T of IBAS is set to 1. The fixed ratio of the distance between the centre of mass and the antennae of beetle to step size is index, which is set to 2 in this work. The initial factor of step size α is set to the default value 0.95. The maximum number of iterations N1 is 200, the maximum number of annealing cycles L is 100, and the threshold of iteration termination is set to 10−20.In order to investigate noise robustness of the aforementioned modal identification and damage identification methods, 0.2% and 0.5% noise have been added into the dynamic time-history response data extracted. LetY∈ℝm×N denote vibration response such as displacement and velocity. Y is contaminated by white noise in terms of(17)Y^l=Yl+d1⋅nlevel⋅σYl,where d1∈ℝ1×N is a vector of normally distributed random values with a zero mean and a unit variance, nlevel denotes the white noise level, and σYl indicates the standard deviation of vibration response at lth measurement point [7].This proposed method was implemented on a PC with an Intel I7 Processor 2.20 GHz and 32 GB RAM.


### 4.2. Single Damage Case

The finite element model of a simply supported beam with single damage is shown in Figure2; the serial number is the measurement point number. In transient dynamic analysis, the load step is 0.001 s, and the duration of the excitation is 50 s. The random excitations are applied at the second to the 25th measurement points in the beam. The displacement response and velocity response in the random excitation process at the second to the 25th measurement points are extracted with the sampling frequency of 200 Hz. In Figure 3, the diagonal elements of the matrix Rqq obtained from the dynamic response with 0.5% noise interference are plotted. The diagonal elements of the matrix Rqq obtained from dynamic response data without noise interference and with 0.2% noise interference are similar to Figure 3; thus, they will not be reiterated here. From the distribution of diagonal elements of Rqq, it is shown that the first two modes are excited with the most energy. To detect the damaged elements, only the first two modes are used for damage identification.Figure 2
Schematic diagram of a simply supported beam with single damage.Figure 3
Diagonal element values of a simply supported beam with single damage (0.5% for the noise level).The theoretical frequencies, structural frequencies obtained by FEM modal, and frequencies identified by SOD method are shown in Table3. These frequencies of the simply supported beam with single damage show that the difference of the first two structural frequencies obtained by SOD method is very small when the dynamic response data are under different levels of noise disturbance. Therefore, SOD method can overcome the noise interference and identify accurately structural frequencies of the simply supported beam with single damage accurately.Table 3
Comparison of structural frequencies of the simply supported beam.
Damage caseNoise levelModeStructural frequencyTheoryFEMSODSingle damageWithout noise124.1824.1424.14297.9797.2097.200.2%124.1824.1424.14297.9797.2097.290.5%124.1824.1424.12297.9797.2097.15No damageWithout noise124.54——298.74——The MAC values [26] based on mode shapes identified by the SOD compared to mode shapes extracted by FEM are all close to 1. The mode shapes identified by the SOD method and mode shapes extracted by FEM modal analysis are shown in Figure 4. It is shown that SOD method can still identify the mode shapes of a single damage beam even after adding different levels of noise into the dynamic response data. Thus, the SOD method has certain noise robustness. Since the diagrams of mode shapes obtained from the dynamic response data without noise interference, with 0.2% noise interference, and 0.5% noise interference are similar, Figure 4 only shows the mode shapes obtained from the dynamic response data with 0.5% noise interference.Figure 4
Mode shapes of a simply supported beam with single damage (0.5% for the noise level). (a) First mode. (b) Second mode.
(a)(b)The modal curvature of each mode is obtained from mode shape by using a second-order central difference method. The smoothed modal curvature of each mode is obtained by using the adaptive polynomial fit method to fit the modal curvatures. Modal curvature values and smoothed modal curvature values at measurement points 2 to 25 in the simply supported beam are shown in Figure5. The figure shows that the first two modal curvatures and smoothed modal curvatures of the simply supported beam are very different at both ends of the damaged element, while there are negligible differences at other measurement points. Because the diagrams of modal curvatures and smoothed modal curvatures obtained from the dynamic response data without noise interference, with 0.2% noise interference, and 0.5% noise interference exhibit the same rule, only modal curvatures and smoothed modal curvatures obtained from the dynamic response data with 0.5% noise interference are shown in Figure 5.Figure 5
Mode shape curvatures of a simply supported beam with single damage (0.5% for the noise level). (a) First mode. (b) Second mode.
(a)(b)The damage index values at the second to the 25th measurement points are shown in Figure6. The results indicate that the damage index values at two measurement points near to the damaged element are obviously greater than those at other nearby measurement points. Therefore, the damage identification method based on SOD method can identify the damaged element accurately in the single damage case of a simply supported beam, which also has certain antinoise properties. The results also show that the method can locate damaged elements in the symmetrical structure. After determining that the 10th element is the damaged element, the subsequent optimization problem is reduced to a one-dimensional problem.Figure 6
DI values of a simply supported beam with single damage.The reduction in stiffness of the 10th element is set as the unknown quantity of IBAS algorithm, and the searching range is set as0,1. The first two frequencies have been applied in the fitness function. The results shown in Table 4 indicate that IBAS algorithm can calculate the stiffness loss of the damaged element accurately after reducing the dimension of the optimization problem. The iterative process of IBAS algorithm based on the results from dynamic response data with 0.5% noise interference is shown in Figure 7. The figure shows that the algorithm converges to the actual reduction in stiffness in the medium term.Table 4
Element stiffness reduction ratio.
Damage caseNoise levelDamaged element numberPreset value (%)Calculated value (%)Error (%)Single damageWithout noise103029.890.370.2%103030.110.370.5%103029.870.43Double damageWithout noise83030.351.17183029.950.170.2%83030.551.83183030.642.130.5%83031.324.40183031.434.77Figure 7
The reduction in stiffness calculated by IBAS (single damage case + 0.5% for the noise level).


### 4.3. Double Damage Case

The finite element model of a cantilever beam with double damage is shown in Figure8; the serial number is the measurement point number. In transient dynamic analysis, the load step is 0.001 s, and the duration of the excitation is 50 s. The random excitations are applied at the second to the 26th measurement points in the beam. The displacement response and velocity response in the random excitation process at the second to the 26th measurement points are extracted with the sampling frequency of 200 Hz. In Figure 9, the diagonal elements of the matrix Rqq obtained from the dynamic response with 0.5% noise interference are plotted. The diagonal elements of the matrix Rqq obtained from dynamic response data without noise interference and with 0.2% noise interference are similar to Figure 9; thus, they will not be reiterated here. From the distribution of diagonal elements of Rqq, it is shown that the first two modes are excited with the most energy. To detect the damaged elements, only the first two modes are used for damage identification.Figure 8
Schematic diagram of a cantilever beam with double damage.Figure 9
Diagonal element values of a cantilever beam with double damage (0.5% for the noise level).The theoretical frequencies, structural frequencies obtained by FEM modal, and frequencies identified by SOD method are shown in Table5. These frequencies of the cantilever beam with double damage show that the difference of the first two structural frequencies obtained by SOD method is very small when the dynamic response data are under different levels of noise disturbance. Therefore, SOD method can overcome the noise interference and identify accurately structural frequencies of the cantilever beam with double damage accurately.Table 5
Comparison of structural frequencies of the cantilever beam.
Damage caseNoise levelModeStructural frequencyTheoryFEMSODDouble damageWithout noise18.618.628.62254.0653.9953.980.2%18.618.628.61254.0653.9953.970.5%18.618.628.61254.0653.9953.94No damageWithout noise18.72——254.75——The MAC values based on mode shapes identified by the SOD compared to mode shapes extracted by FEM are all close to 1. The mode shapes identified by the SOD method and mode shapes extracted by FEM modal analysis are shown in Figure10. It is shown that SOD method can still identify the mode shapes of the beam with double damage even after adding different levels of noise into the dynamic response data. Thus, the SOD method has certain noise robustness. Since the diagrams of mode shapes obtained from the dynamic response data without noise interference, with 0.2% noise interference, and with 0.5% noise interference are similar, Figure 10 only shows the mode shapes obtained from the dynamic response data with 0.5% noise interference.Figure 10
Mode shapes of a cantilever beam with double damage (0.5% for the noise level). (a) First mode. (b) Second mode.
(a)(b)The modal curvature of each mode is obtained from mode shape by using a second-order central difference method. The smoothed modal curvature of each mode is obtained by using the adaptive polynomial fit method to fit the modal curvatures. Modal curvature values and smoothed modal curvature values at measurement points 2 to 26 in the cantilever beam are shown in Figure11. The figure shows that the first two modal curvatures and smoothed modal curvatures of the cantilever beam are very different at both ends of the damaged element, while there are negligible differences at other measurement points. Because the diagrams of modal curvatures and smoothed modal curvatures obtained from the dynamic response data without noise interference, with 0.2% noise interference, and with 0.5% noise interference exhibit the same rule, only modal curvatures and smoothed modal curvatures obtained from the dynamic response data with 0.5% noise interference are shown in Figure 11.Figure 11
Mode shape curvatures of a cantilever beam with double damage (0.5% for the noise level). (a) First mode. (b) Second mode.
(a)(b)The damage index values at the second to the 25th measurement points are shown in Figure12. The results indicate that the damage index values at two measurement points near to the damaged element are obviously greater than those at other nearby measurement points. Therefore, the damage identification method based on SOD method can identify the damaged element accurately in the double damage case of the cantilever beam, which also has certain antinoise properties. After determining that the 8th and 18th elements are damaged elements, the subsequent optimization problem is reduced to a two-dimensional problem.Figure 12
DI values of a cantilever beam with double damage.The reduction in stiffness of the 8th and 18th elements are set as the unknown quantities of IBAS algorithm, and the searching ranges are all set as0,1. The first two frequencies have been applied in the fitness function. The results shown in Table 4 indicate that IBAS algorithm can calculate the stiffness loss of the damaged elements accurately after reducing the dimension of the optimization problem. The iterative processes of IBAS algorithm based on the results from dynamic response data with 0.5% noise interference are shown in Figure 13. The red circle indicates that the fitness function value is less than the threshold value and the iteration ceases.Figure 13
The reduction in stiffness calculated by IBAS (double damage case + 0.5% for the noise level). (a) First damaged element. (b) Second damaged element.
(a)(b)


## 4.1. Calculation Model and Parameters

A simply supported beam model and a cantilever beam model both with a span of 0.75 m have been established in ANSYS software. The element type is the Beam3 Euler beam element, and the beam models are divided into 25 elements. The geometric and material parameters are given in Table1. Gaussian white noise arrays are imported to simulate random excitation. The modal analysis of the damaged structure is analyzed firstly to compute the theoretical modal parameters. And then, the undamped transient dynamic analysis is carried out using the full method.Table 1
Dimensions and material properties.
PropertiesMean valueWidth (m)0.025Thickness (m)0.006Young’s modulus (GPa)200Density (kg/m3)7800Poisson’s ratio0.3Numerical simulations of two damage cases in this study are shown in Table2. The damage is simulated by reducing the stiffness of the specific element model. In the process of damage identification, stiffness reduction is realized by reducing elasticity modulus of the damaged element as [25]:(16)Eid=1−aiE0,where Eid is the elasticity modulus of the ith damaged element, E0 denotes the pristine elasticity modulus, and ai denotes the reduction in stiffness of the ith damaged element, which is equal to the ratio between the stiffness reduction and the pristine stiffness, and it is also called damage degree in this paper.Table 2
Damage cases.
Damage caseDamaged element numberElement stiffness reductionSingle damage1030%Double damage8, 1830%, 30%The element stiffness matrix and element mass matrix of the middle element of the undamaged equal section beam model are used to assemble the global stiffness matrix and global mass matrix, respectively. And the global stiffness matrix and global mass matrix that are imported to compute the reduction in stiffness of the damaged element by iterating are used to calculate the value of damaged structural frequencyωir by finite element theory at each iteration of IBAS algorithm. The initial temperature T of IBAS is set to 1. The fixed ratio of the distance between the centre of mass and the antennae of beetle to step size is index, which is set to 2 in this work. The initial factor of step size α is set to the default value 0.95. The maximum number of iterations N1 is 200, the maximum number of annealing cycles L is 100, and the threshold of iteration termination is set to 10−20.In order to investigate noise robustness of the aforementioned modal identification and damage identification methods, 0.2% and 0.5% noise have been added into the dynamic time-history response data extracted. LetY∈ℝm×N denote vibration response such as displacement and velocity. Y is contaminated by white noise in terms of(17)Y^l=Yl+d1⋅nlevel⋅σYl,where d1∈ℝ1×N is a vector of normally distributed random values with a zero mean and a unit variance, nlevel denotes the white noise level, and σYl indicates the standard deviation of vibration response at lth measurement point [7].This proposed method was implemented on a PC with an Intel I7 Processor 2.20 GHz and 32 GB RAM.


## 4.2. Single Damage Case

The finite element model of a simply supported beam with single damage is shown in Figure2; the serial number is the measurement point number. In transient dynamic analysis, the load step is 0.001 s, and the duration of the excitation is 50 s. The random excitations are applied at the second to the 25th measurement points in the beam. The displacement response and velocity response in the random excitation process at the second to the 25th measurement points are extracted with the sampling frequency of 200 Hz. In Figure 3, the diagonal elements of the matrix Rqq obtained from the dynamic response with 0.5% noise interference are plotted. The diagonal elements of the matrix Rqq obtained from dynamic response data without noise interference and with 0.2% noise interference are similar to Figure 3; thus, they will not be reiterated here. From the distribution of diagonal elements of Rqq, it is shown that the first two modes are excited with the most energy. To detect the damaged elements, only the first two modes are used for damage identification.Figure 2
Schematic diagram of a simply supported beam with single damage.Figure 3
Diagonal element values of a simply supported beam with single damage (0.5% for the noise level).The theoretical frequencies, structural frequencies obtained by FEM modal, and frequencies identified by SOD method are shown in Table3. These frequencies of the simply supported beam with single damage show that the difference of the first two structural frequencies obtained by SOD method is very small when the dynamic response data are under different levels of noise disturbance. Therefore, SOD method can overcome the noise interference and identify accurately structural frequencies of the simply supported beam with single damage accurately.Table 3
Comparison of structural frequencies of the simply supported beam.
Damage caseNoise levelModeStructural frequencyTheoryFEMSODSingle damageWithout noise124.1824.1424.14297.9797.2097.200.2%124.1824.1424.14297.9797.2097.290.5%124.1824.1424.12297.9797.2097.15No damageWithout noise124.54——298.74——The MAC values [26] based on mode shapes identified by the SOD compared to mode shapes extracted by FEM are all close to 1. The mode shapes identified by the SOD method and mode shapes extracted by FEM modal analysis are shown in Figure 4. It is shown that SOD method can still identify the mode shapes of a single damage beam even after adding different levels of noise into the dynamic response data. Thus, the SOD method has certain noise robustness. Since the diagrams of mode shapes obtained from the dynamic response data without noise interference, with 0.2% noise interference, and 0.5% noise interference are similar, Figure 4 only shows the mode shapes obtained from the dynamic response data with 0.5% noise interference.Figure 4
Mode shapes of a simply supported beam with single damage (0.5% for the noise level). (a) First mode. (b) Second mode.
(a)(b)The modal curvature of each mode is obtained from mode shape by using a second-order central difference method. The smoothed modal curvature of each mode is obtained by using the adaptive polynomial fit method to fit the modal curvatures. Modal curvature values and smoothed modal curvature values at measurement points 2 to 25 in the simply supported beam are shown in Figure5. The figure shows that the first two modal curvatures and smoothed modal curvatures of the simply supported beam are very different at both ends of the damaged element, while there are negligible differences at other measurement points. Because the diagrams of modal curvatures and smoothed modal curvatures obtained from the dynamic response data without noise interference, with 0.2% noise interference, and 0.5% noise interference exhibit the same rule, only modal curvatures and smoothed modal curvatures obtained from the dynamic response data with 0.5% noise interference are shown in Figure 5.Figure 5
Mode shape curvatures of a simply supported beam with single damage (0.5% for the noise level). (a) First mode. (b) Second mode.
(a)(b)The damage index values at the second to the 25th measurement points are shown in Figure6. The results indicate that the damage index values at two measurement points near to the damaged element are obviously greater than those at other nearby measurement points. Therefore, the damage identification method based on SOD method can identify the damaged element accurately in the single damage case of a simply supported beam, which also has certain antinoise properties. The results also show that the method can locate damaged elements in the symmetrical structure. After determining that the 10th element is the damaged element, the subsequent optimization problem is reduced to a one-dimensional problem.Figure 6
DI values of a simply supported beam with single damage.The reduction in stiffness of the 10th element is set as the unknown quantity of IBAS algorithm, and the searching range is set as0,1. The first two frequencies have been applied in the fitness function. The results shown in Table 4 indicate that IBAS algorithm can calculate the stiffness loss of the damaged element accurately after reducing the dimension of the optimization problem. The iterative process of IBAS algorithm based on the results from dynamic response data with 0.5% noise interference is shown in Figure 7. The figure shows that the algorithm converges to the actual reduction in stiffness in the medium term.Table 4
Element stiffness reduction ratio.
Damage caseNoise levelDamaged element numberPreset value (%)Calculated value (%)Error (%)Single damageWithout noise103029.890.370.2%103030.110.370.5%103029.870.43Double damageWithout noise83030.351.17183029.950.170.2%83030.551.83183030.642.130.5%83031.324.40183031.434.77Figure 7
The reduction in stiffness calculated by IBAS (single damage case + 0.5% for the noise level).


## 4.3. Double Damage Case

The finite element model of a cantilever beam with double damage is shown in Figure8; the serial number is the measurement point number. In transient dynamic analysis, the load step is 0.001 s, and the duration of the excitation is 50 s. The random excitations are applied at the second to the 26th measurement points in the beam. The displacement response and velocity response in the random excitation process at the second to the 26th measurement points are extracted with the sampling frequency of 200 Hz. In Figure 9, the diagonal elements of the matrix Rqq obtained from the dynamic response with 0.5% noise interference are plotted. The diagonal elements of the matrix Rqq obtained from dynamic response data without noise interference and with 0.2% noise interference are similar to Figure 9; thus, they will not be reiterated here. From the distribution of diagonal elements of Rqq, it is shown that the first two modes are excited with the most energy. To detect the damaged elements, only the first two modes are used for damage identification.Figure 8
Schematic diagram of a cantilever beam with double damage.Figure 9
Diagonal element values of a cantilever beam with double damage (0.5% for the noise level).The theoretical frequencies, structural frequencies obtained by FEM modal, and frequencies identified by SOD method are shown in Table5. These frequencies of the cantilever beam with double damage show that the difference of the first two structural frequencies obtained by SOD method is very small when the dynamic response data are under different levels of noise disturbance. Therefore, SOD method can overcome the noise interference and identify accurately structural frequencies of the cantilever beam with double damage accurately.Table 5
Comparison of structural frequencies of the cantilever beam.
Damage caseNoise levelModeStructural frequencyTheoryFEMSODDouble damageWithout noise18.618.628.62254.0653.9953.980.2%18.618.628.61254.0653.9953.970.5%18.618.628.61254.0653.9953.94No damageWithout noise18.72——254.75——The MAC values based on mode shapes identified by the SOD compared to mode shapes extracted by FEM are all close to 1. The mode shapes identified by the SOD method and mode shapes extracted by FEM modal analysis are shown in Figure10. It is shown that SOD method can still identify the mode shapes of the beam with double damage even after adding different levels of noise into the dynamic response data. Thus, the SOD method has certain noise robustness. Since the diagrams of mode shapes obtained from the dynamic response data without noise interference, with 0.2% noise interference, and with 0.5% noise interference are similar, Figure 10 only shows the mode shapes obtained from the dynamic response data with 0.5% noise interference.Figure 10
Mode shapes of a cantilever beam with double damage (0.5% for the noise level). (a) First mode. (b) Second mode.
(a)(b)The modal curvature of each mode is obtained from mode shape by using a second-order central difference method. The smoothed modal curvature of each mode is obtained by using the adaptive polynomial fit method to fit the modal curvatures. Modal curvature values and smoothed modal curvature values at measurement points 2 to 26 in the cantilever beam are shown in Figure11. The figure shows that the first two modal curvatures and smoothed modal curvatures of the cantilever beam are very different at both ends of the damaged element, while there are negligible differences at other measurement points. Because the diagrams of modal curvatures and smoothed modal curvatures obtained from the dynamic response data without noise interference, with 0.2% noise interference, and with 0.5% noise interference exhibit the same rule, only modal curvatures and smoothed modal curvatures obtained from the dynamic response data with 0.5% noise interference are shown in Figure 11.Figure 11
Mode shape curvatures of a cantilever beam with double damage (0.5% for the noise level). (a) First mode. (b) Second mode.
(a)(b)The damage index values at the second to the 25th measurement points are shown in Figure12. The results indicate that the damage index values at two measurement points near to the damaged element are obviously greater than those at other nearby measurement points. Therefore, the damage identification method based on SOD method can identify the damaged element accurately in the double damage case of the cantilever beam, which also has certain antinoise properties. After determining that the 8th and 18th elements are damaged elements, the subsequent optimization problem is reduced to a two-dimensional problem.Figure 12
DI values of a cantilever beam with double damage.The reduction in stiffness of the 8th and 18th elements are set as the unknown quantities of IBAS algorithm, and the searching ranges are all set as0,1. The first two frequencies have been applied in the fitness function. The results shown in Table 4 indicate that IBAS algorithm can calculate the stiffness loss of the damaged elements accurately after reducing the dimension of the optimization problem. The iterative processes of IBAS algorithm based on the results from dynamic response data with 0.5% noise interference are shown in Figure 13. The red circle indicates that the fitness function value is less than the threshold value and the iteration ceases.Figure 13
The reduction in stiffness calculated by IBAS (double damage case + 0.5% for the noise level). (a) First damaged element. (b) Second damaged element.
(a)(b)


## 5. Conclusions

In this work, a two-stage damage identification method of structure that utilizes SOD combined with IBAS algorithm has been studied. In the first stage, the SOD method has been introduced to identify the damage localization firstly. Specifically, the curvatures of the mode shapes identified by SOD method are calculated by using the second-order central difference method. Then, the smoothed modal curvatures are obtained by using the adaptive polynomial fit method to fit modal curvatures. Finally, the damage index values are calculated to locate damages. In the second stage, after eliminating the healthy elements, the IBAS algorithm proposed by Zhou et al. is used to calculate the reduction in stiffness of the damaged elements based on the modal parameters obtained by SOD method and the results of damage localization. Numerical examples of two damage cases indicate that the damage identification method of beam structure based on SOD method and IBAS algorithm can locate the damaged elements and calculate the stiffness loss of the damaged elements. The damage identification method is accurate even for noise-contaminated data. IBAS makes BAS possible to apply for the optimization problems of large-scale structural models. Besides, the results of the single damage case show that this combined application also can locate damaged elements in the symmetrical structure. The calculation time of each numerical example is within 10 s, which indicates that the proposed method is numerically efficient. However, The accuracy of the proposed approach under strong noise interference needs to be further studied. And further research on the efficiency of the proposed approach applied for other complex structures should be conducted.



---

*Source: 8857356-2021-02-27.xml*