A two-step method for damage identification in beam structures based on influence line difference and acceleration data

Abstract

This article presents a two-step approach for structural damage identification in beam structure. Damages are located using the influence line difference before and after damage, the calculation of damage severity is accomplished by acceleration data and bird mating optimizer algorithm. Local damages are simulated as the reduction of both the elemental Young’s modulus and mass of the beam. The technique for damage localization based on displacement influence line difference and its derivatives for beam structure has been outlined. An objective function that comprises dynamic acceleration is utilized in bird mating optimizer. All data are originated from only a few measurement points. Two numerical examples, namely, a simply supported beam and a four-span continuous beam, are investigated in this article. Identification results from different objective functions are compared with results from objective function conventional modal assurance criterion, which shows the superiority of the proposed function. In addition, results of dynamic responses under different types of excitation are presented. The effect of measurement noise level on damage identification results is studied. Studies in the article indicate that the proposed method is efficient and robust for identifying damages in beam structures.

Introduction

It is well known that structures are suffering from all kinds of damage problems due to both the unpredictable internal and external factors that jeopardize the safety of the structure; thus, a reliable structural damage identification method is needed.

Generally, damages are quantified by the reduction of stiffness parameter,^1,2 which affects the dynamic property of the structure, such as natural frequencies and mode shapes. In the last few decades, many damage identification methods had been developed using the vibration data in frequency domain. Cawley and Adams ³ first localized the damage from natural frequency of structure. Pandey et al. ⁴ and Rizos et al. ⁵ identified the local damage by curvature of mode shapes. And flexibility of structure^6–9 is also used. Lee ¹⁰ applied frequencies in iteration to detect the location and severity of damages. Zhong and Oyadiji ¹¹ detected damages by combining wavelet and modal data. On the contrary, a great deal of methods based on data in time domain had also been developed. Majumder and Manohar ¹² proposed a time domain approach for beam structures using vibration data. Zhu and Law ¹³ applied wavelet to deflection of time history to identified damages. Law and Lu ¹⁴ proposed a method based on dynamic response and later presented response sensitivity–based model updating method to identify structural damage and parameter.^15,16 Recently, Khorram et al. ¹⁷ combined wavelet and factorial design to detect the crack from moving load responses. In addition, two-step approach was developed by Lu et al. ¹⁸ Based on curvature of mode shape, the crack location is predicted by a damage indicator, and then response sensitivity–based model updating is applied to identify the crack depth and location further. More recently, a two-step method based on modal strain energy and response sensitivity extended to identify the damages in isotropic plates by Fu et al. ¹⁹

In mathematical aspect, damage identification problem can be regarded as an optimization problem that aims to find the optimal solution of an objective function which manifests the difference between real structure and simulation. When adopting conventional methods to solve these problems, it is easy to reach a local optimum and are sensitive to initial value. More attention was paid in heuristic algorithm these years, Wu et al. ²⁰ is among the earliest one who detect damage using neural network, and the method was developed further later.^21–23 Majumdar et al. ²⁴ assessed damage of truss from changes in natural frequencies using ant colony optimization (ACO). Maity and Tripathy ²⁵ used genetic algorithm (GA) to identify the damage from changes in natural frequencies, and Ravanfar et al. ²⁶ utilized GA as the second step in his two-step method to assess the severity of the damage. Mohan et al. ²⁷ compared the particle swarm optimization (PSO) and GA in the efficiency and accuracy of identification. Seyedpoor ²⁸ developed a two-step method using PSO for damage identification. Sokhangou et al. ²⁹ made use of pseudo-residual force vector to determine the damage element, and then applied gravitational search algorithm (GSA) to get the damage extent. Chaotic artificial bee colony (CABC) algorithm derived from artificial bee colony (ABC) was used to detect damage by Xu et al. ³⁰

A new heuristic algorithm named bird mating optimizer (BMO) was proposed by Askarzadeh^31,32 and Askarzadeh and Rezazadeh^33,34 recently, which imitates the mating behavior of birds. The algorithm was simulated and compared with GA, PSO, and GSO ³¹ and had better accuracy, Askarzadeh further estimate parameters with the most influence on the electrochemical model, ³² modeled the fuel cell system, ³⁴ and extracted the maximum power point in solar cells. ³³ Li et al. ³⁵ utilized the method in damage to optimize the objective and optimum set of stiffness reduction parameters. Zhu et al. ³⁶ designed a disturbance procedure in BMO which enhances the precision of damage identification.

As it is mentioned above that some two-step approaches are applied in damage identification area, the localization of the damage greatly reduces the parameter dimension, which greatly reduces the amount of calculation; hence, the efficiency will improved. Seyedpoor ²⁸ localized the damages from modal strain energy–based index and assess them through PSO. Sokhangou et al. ²⁹ used pseudo-residual force vector to determine the damage element. Grandić et al. ³⁷ used the change in the rotation of displacement influence lines to indicate damage. Wang and colleagues^38,39 located and quantified damage by strain statistical moment in both beam and plate structure.

In this article, a two-step method for damage identification is proposed. Displacement influence line will be utilized to determine the damage location, which is conducted by difference quotient in finite element model. Then, the damage is quantified using BMO in the second step. An objective function is established by minimizing the discrepancies between the simulated “measured” dynamic acceleration responses and calculated ones using direct integration. As only a few number of measurement points are needed, and it has the advantages making of plenty of time domain data, it is convenient in practical applications. A simply supported beam and a four-span beam are studied to show the promising results of the method, and the influence of measurement noise on damage identification results is also investigated in the numerical simulations. In addition, comparison among the proposed objective function and other objective functions consisted of modal assurance criterion (MAC) and dynamic displacement is conducted to validate the superiority of the dynamic acceleration–based objective function.

Damage detection methodology

Damages are normally defined as the deduction of stiffness, beside they can be simulated in many forms. As it is shown in Figure 1, a rectangular portion of the cross section was hollowed out, and this damage will cause the stiffness to reduce and mass loss. Given this premise, we can quantify the local damage as the reduction of both stiffness and mass in some parts.

Figure 1.

A damage model with rectangular crack.

Damage model construction

In this article, analytical models of beam structures are established by finite element method (FEM). When damage occurred in a structure, the extent of reduction in stiffness and mass can be represented by elemental damages with parameters $\alpha$ and $\beta$ between 0 and 1. For a structure with $n$ elements, we set n-dimensional vectors $\overrightarrow{\alpha }$ and $\overrightarrow{\beta }$ , where ${\alpha }_i$ and ${\beta }_i$ $(i=1,2,\dots ,n)$ stand for the damage parameters of stiffness and mass of the ith element, respectively. Thus, the system stiffness and mass matrices after damage can be written as follows

$${\mathbf{K}}_{\mathbf{d}}=\sum _{i=1}^n{\alpha }_i{\mathbf{k}}_{\mathbf{i}}^{\mathbf{e}},{\mathbf{M}}_{\mathbf{d}}=\sum _{i=1}^n{\beta }_i{\mathbf{m}}_{\mathbf{i}}^{\mathbf{e}}$$

(1)

where ${\mathbf{k}}_{\mathbf{i}}^{\mathbf{e}}$ and ${\mathbf{m}}_{\mathbf{i}}^{\mathbf{e}}$ represent the elemental stiffness and mass matrices for element $i$ , and ${\mathbf{K}}_{\mathbf{d}}$ and ${\mathbf{M}}_{\mathbf{d}}$ as global stiffness and mass matrices, respectively. Normally, the damage extents for stiffness and mass may not be closely correlated and usually not the same. When ${\alpha }_i={\beta }_i=1$ , the corresponding element is intact.

Now, the free vibration equation of damage beam is acquired

$$\mathbf{M}\stackrel{··}{\mathbf{d}}+\mathbf{C}\stackrel{·}{\mathbf{d}}+\mathbf{Kd}=\mathbf{F}\left(t\right)$$

(2)

where $\mathbf{M}$ , $\mathbf{C}$ , and $\mathbf{K}$ are the system mass, damping, and stiffness matrices, respectively. Rayleigh damping is adopted that $\mathbf{C}={\eta }_1\mathbf{M}+{\eta }_2\mathbf{K}$ , and ${\eta }_1$ and ${\eta }_2$ are two constants determined from two given damping ratios that correspond to the first two modal frequencies of the structure. For the damaged structure in the study, the calculation is conducted by replacing the above matrices by ${\mathbf{M}}_{\mathbf{d}}$ , ${\mathbf{C}}_{\mathbf{d}}$ , and ${\mathbf{K}}_{\mathbf{d}}$ .

And the natural frequencies and the associated mode shapes can be obtained from the following eigen equation

$$(\mathbf{K}-{\omega }_j^2\mathbf{M}){\varphi }_j=0$$

(3)

where ${\omega }_j$ is the $j$ th natural frequency and ${\varphi }_j$ is the associated mode shape.

Localization of the damage by influence line difference

The reduction of stiffness and mass will cause deviation in responses of the structure. In a beam structure with length $l$ , if a force $P$ is acting on it, we can obtain the deflection of point located at $x_s$ by equation (4)

$$f(x_s)=\underset{0}{\overset{l}{\int }}\frac{M_P(x)\overline{M}(x_s,x)}{\mathrm{EI}(x)}\mathrm{dx}$$

(4)

where $M_P(x)$ is the bending moment through the beam when only one external force $P$ is acting on the beam, $\overline{M}(x_s,x)$ is the bending moment when only a unit force acting at the measuring point $x_s$ , respectively, and $f$ is the deflection function of the beam.

If we assume that there is a damage existing in the beam between $a$ and $b$ whose stiffness damage index is $\alpha$ , then the stiffness value between $a$ and $b$ is reduced to $\alpha \mathrm{EI}$ , where $\mathrm{EI}$ is the stiffness for undamaged part. Equation (4) can be divided to

$$f(x_s)={\displaystyle \underset{0}{\overset{a}{\int }}\frac{M_p(x)\overline{M}(x_s,x)}{EI}dx}+{\displaystyle \underset{a}{\overset{b}{\int }}\frac{M_p(x)\overline{M}(x_s,x)}{\alpha EI}dx}+{\displaystyle \underset{b}{\overset{l}{\int }}\frac{M_p(x)\overline{M}(x_s,x)}{EI}dx}$$

(5)

Subtracting equation (5) by equation (4), we obtain the following

$$\Delta f(x_s)=\underset{a}{\overset{b}{\int }}\frac{M_p(x)\overline{M}(x_s,x)}{\mathrm{EI}}(\frac{1}{\alpha }-1)\mathrm{dx}$$

(6)

If we set the external force $P$ as a single moving load with location $x_p$ , the displacement influence line value for point $x_s$ is formed as a two-dimensional function

$$\Delta f(x_p,x_s)=\underset{a}{\overset{b}{\int }}\frac{M_p(x_p,x)\overline{M}(x_s,x)}{\mathrm{EI}}(\frac{1}{\alpha }-1)\mathrm{dx}$$

(7)

Judging from the integral range of equation (7), the value of $M_p(x_p,x)$ within $[a,b]$ will change with external force moving from left side of $a$ to right side of $b$ , which is illustrated by Figure 2. In mechanical aspect, the change is caused directly by the change of shear force, noticing that there is no concentrated moment acting on the beam, so that $M_p(x_p,x)$ is continuous along $x$ direction. Taking the partial derivative on both sides of equation (7) with respective to $x_p$ , we obtain

$$\begin{array}{cc}\multicolumn{1}{c}{\Delta f_{x_p}(x_p,x_s)}& =\frac{\partial \underset{a}{\overset{b}{\int }}\frac{M_p(x_p,x)\overline{M}(x_s,x)}{\mathrm{EI}}(\frac{1}{\alpha }-1)\mathrm{dx}}{\partial x_p}\hfill \\ \multicolumn{1}{c}{}& =\underset{a}{\overset{b}{\int }}\frac{\partial M_p(x_p,x)\overline{M}(x_s,x)}{\mathrm{EI}\partial x_p}(\frac{1}{\alpha }-1)\mathrm{dx}\hfill \\ \multicolumn{1}{c}{}& =\underset{a}{\overset{b}{\int }}\frac{F_{\mathrm{sp}}(x_p,x)\overline{M}(x_s,x)}{\mathrm{EI}}(\frac{1}{\alpha }-1)\mathrm{dx}\hfill \end{array}$$

(8)

where $F_{\mathrm{sp}}(x_p,x)$ represents the shear force of the beam caused by moving load $P(x_p)$ . A sudden change of $\Delta f_{x_p}(x_s,x_p)$ will occur within $[a,b]$ , which can be explained as follows: when the moving force passes through a damaged area, the shear force within the area would change from positive to negative illustrated in Figure 3.

Figure 2.

The change of bending moment before and after a force is through the point.

Figure 3.

The change of shear force before and after a force is through the point.

For a structure with multiple damages, the values of $\Delta f_{x_p}(x_s,x_p)$ will be the linear sum of derivative of influence line difference (DILD) values of every damage individual in equation (8). For a system with $n$ elements with elemental length (1/n), the discrete influence line difference originated from equation (7) is as follows

$$\Delta f(x_p,x_s)=\sum _{i=1}^{\mathrm{ND}}\frac{M_p(x_p,i)\overline{M}(x_s,i)l}{\mathrm{nEI}}(\frac{1}{{\alpha }_i}-1)$$

(9)

where $\mathrm{ND}$ is the number of damages and ${\alpha }_i$ is the stiffness damage index of damage of node $i$ . Note that $\Delta f(x_p,x_s)$ is proportional to $x_s$ , to magnify the value of $\Delta f(x_p,x_s)$ and make the damage notable, $x_s$ , namely, the measurement location should not be adjacent to the supporting point.

To further calculating the DILD in discrete model, forward difference quotient scheme is applied, and DILD values are calculated by equation (10)

$$\Delta f_{x_p}(n_i)=\mathrm{DILD}(n_i)=\frac{\Delta f(n_{i+1})-\Delta f(n_i)}{l_e}$$

(10)

where $l_e$ is the distance between two adjacent nodes. By using forward difference quotient, we need to neglect the value of last node, which as the supporting point will not set to be damaged in damage detection.

One can further specify the change of DILD by taking derivative to DILD with respect to x_p again. In the analytical FE model, the second-order central difference quotient is used to calculate the second derivative. Difference quotient of DILD (DDILD) is calculated in equation (11)

$$DDILD(n_i)=\frac{DILD(n_i)-DILD(n_{i-1})}{l_e}=\frac{\Delta f(n_{i+1})-2\Delta f(n_i)+\Delta f(n_{i-1})}{l_e^2}$$

(11)

And it is obvious that the value of DDILD will not be zero if the corresponding point is damaged.

Note that in this study, a damage is assumed as an elemental damage associated with two nodes; on the contrary, the values of $\Delta f_{x_p}$ are located at each node, so for a single damage, the localization will be presented as two consecutive non-zero values in two nodes in the curve of DDILD.

Objective function for detection assessment

Since the localization of the damages by DDILD is based on the affection of the stiffness damage but not the mass damage, the DDILD values is not available on damage level assessment. Normally, the damage identification problem can be numerically transferred into an optimization, the goal is to minimize the difference of property between simulated and actual structure.

Damage always results in changes of property in structure, responses such as vibration frequency, mode or dynamic responses, are frequently used in damage identification.

MAC, as shown in equation (12), is usually the object function in damage assessment, where $\varphi$ is mode, the superscripts $C$ and $M$ indicate the data from calculation and measurement. In practice, only a few lower orders of mode shapes are available in identification

$$\mathrm{MAC}=\frac{{({\left[{\varphi }_j^C\right]}^T•\left[{\varphi }_j^M\right])}^2}{{\left|\left[{\varphi }_j^C\right]\right|}^2{\left|\left[{\varphi }_j^M\right]\right|}^2}$$

(12)

In addition to data from frequency domain, there will be influences on dynamic responses in time history to change, such as acceleration that contains numerous damage information. To node number $n$ , the acceleration responses of $m$ time steps can be expressed as a vector

$$\mathbf{ac}{\mathbf{c}}_n={\left[\begin{array}{cccc}{a}_n\left(t_1\right)\hfill & a_n\left(t_2\right)\hfill & \cdots \hfill & a_n\left(t_m\right)\hfill \end{array}\right]}^T$$

(13)

By measuring the deviance of acceleration responses between actual and calculation in the form of MAC, function accMAC is defined as follows

$$\mathrm{accMA}{C}_n=\frac{{({\left[{\mathbf{acc}}_n^C\right]}^T•\left[{\mathbf{acc}}_n^M\right])}^2}{{\left|\left[{\mathbf{acc}}_n^C\right]\right|}^2{\left|\left[{\mathbf{acc}}_n^M\right]\right|}^2}$$

(14)

where the superscripts $C$ and $M$ denote that the responses originated from calculation and measurement, respectively. When the calculation is in accordance with the measured ones, it indicates a correct assessment and the accMAC reaches 1, so the objective function is as follows

$$\mathrm{Objfunc}=\sum _{j=1}^{\mathrm{NA}}{w}_{\mathrm{aj}}^2[1-\frac{{({\left[{\mathbf{acc}}_j^C\right]}^T•\left[{\mathbf{acc}}_j^M\right])}^2}{{\left|\left[{\mathbf{acc}}_j^C\right]\right|}^2{\left|\left[{\mathbf{acc}}_j^M\right]\right|}^2}]$$

(15)

Where $w_{\varphi }$ is the weighting factor of the acceleration data and there are $\mathrm{NA}$ measuring points in total.

So the goal is to minimize equation (15). Compared to taking the least-squares of deviation between simulated and actual responses, accMAC prevents the inconformity of two acceleration vectors because acceleration with larger absolute values would weigh heavier in least-squares scheme.

BMO

BMO is an evolutionary algorithm proposed by Alireza Askarzadeh, which imitates the breeding behavior of bird spices, and its core idea is to pass on superior genes to next generations.

In this algorithm, a feasible solution is called a bird, and there are two kinds of birds in a society, male and female. The breeding process is equivalent to searching process in the algorithm. Individuals are categorized into five mating patterns judging by sex and gene quality, with two for female and three for male. Their traits to produce offspring are introduced as follows.

1. Monogamy

   A monogamy male bird $\overrightarrow{x}$ mates with his interested female $\overrightarrow{x_i}$ , a resultant brood is produced by

   $$\begin{array}{c}{\overrightarrow{x}}_b=\overrightarrow{x}+w\times \overrightarrow{r}.\times ({\overrightarrow{x}}_i-\overrightarrow{x})\hfill \\ c=\mathrm{a\; random\; number\; between}1\mathrm{and}n\hfill \\ \mathrm{if}{r}_1>\mathrm{mcf}\hfill \\ x_b(c)=l(c)-r_2\times (l(c)-u(c));\hfill \\ \mathrm{end}\hfill \end{array}$$

   (16)

   where $\overrightarrow{x_b}$ is the resultant brood, $w$ is a time-varying weighing factor, $\overrightarrow{r}$ is a vector whose elements is distributed randomly in [0, 1], the dimension of the bird individual is $n$ , $r_i$ is a random number between 0 and 1, $\mathrm{mcf}$ is a mutation control factor varying between 0 and 1, and $u$ and $l$ are the upper and lower bounds of the solution, respectively.

2. Polygyny

   A polygyny bird has several interested female birds and the resultant brood $\overrightarrow{x_b}$ is given by equation (17)

   $$\begin{array}{c}{\overrightarrow{x}}_b=\overrightarrow{x}+w\times \sum _{j=1}^{n_i}\overrightarrow{r}.\times ({\overrightarrow{x}}_{\mathrm{ij}}-\overrightarrow{x})\hfill \\ c=\mathrm{a\; random}\mathrm{number\; between}1\mathrm{and}n\hfill \\ \mathrm{if}{r}_1>\mathrm{mcf}\hfill \\ x_b(c)=l(c)-r_2\times (l(c)-u(c));\hfill \\ \mathrm{end}\hfill \end{array}$$

   (17)

   where $n_i$ is the number of interested birds, and $x_{\mathrm{ij}}$ represents the jth interested bird.

3. Parthenogenesis

   The best solutions play the role as parthenogenetic birds. Each parthenogenetic bird produces a brood by the following equation

   $$\begin{array}{c}\mathrm{for}i=1:n\hfill \\ \mathrm{if}{r}_1>\mathrm{mc}{f}_p\hfill \\ x_b(i)=x(i)+\mu \times (r_2-r_3)\times x(i);\hfill \\ \mathrm{else}\hfill \\ x_b(i)=x(i);\hfill \\ \mathrm{end}\hfill \\ \mathrm{end}\hfill \end{array}$$

   (18)

   where $\mathrm{mc}{f}_p$ is the parthenogenetic mutation control factor and $\mu$ is the step size.

4. Promiscuity

   Promiscuity birds mutated from birds with worst genes judging by the fitness function; this group is produced by chaotic sequence. And the way they breed is the same as that for monogamous birds, that is, equation (16).

5. Polyandry

   In this group, each female will mate with several males, so the resultant brood is produced by equation (17) too.

After every mating process, the fitness values of all the resultant broods are evaluated through a fitness function, the group with better individuals will substitute the worse one and prepare to breed the next generation.

Damage detection procedure

Based on BMO algorithm and in combination with the influence line difference, the damage identification procedures are shown as follows, also illustrated in Figure 4:

• Step 1. Measurement: The influence line before and after the damage (deflection data) and responses of the damaged structure are measured and recorded.

• Step 2. Localization: Using influence line difference data and calculated DDILD to localize the damage position, a few times of measurements can be applied till a stable trend is eminent in the curve, and the damaged elements are selected.

• Step 3. Initialization: A set of damage parameter vectors with dimension of damage number are randomly set in the searching space $[0,1]$ .

• Step 4. Ranking: Each parameter vector ranks according to the corresponding gene quality judging by the fitness value.

• Step 5. Classification: Vectors with the best fitness value are selected as female birds and the others are chosen as males and the ratio will be preset in the simulation. There is a ratio 50:30:10:5:5 ³¹ to monogamy, polygyny, parthenogenesis, polyandry, and promiscuity recommended by the author. Then females are equally divided into two groups that the better ones make parthenogenetic birds and the others make polyandrous group. The males are categorized into three groups. The males with better quality are selected as monogamous, the worst males are substituted by new individuals generated by chaotic sequence as the promiscuous, while the remaining ones are chosen as polygynous.

• Step 6. Breeding: Each bird breeds the resultant brood by its own way.

• Step 7. Replacement: Evaluate the quality of new and existed birds by fitness function, the number of a society is limited, and the better ones are remained to form a new society.

• Step 8. Searching: Loop steps 5–7 until a fitness value meets with the accuracy threshold, or the number of evolution reaches the maximum generation.

• Step 9. The bird with best quality is considered to be the identification result.

Figure 4.

Flow chart of the proposed method.

Numerical simulation

Parameter of BMO

The basic parameters of BMO are presented as follows:

• The society scale is 100, which means 100 damage parameter vectors are set in the assessment. The number of monogamous, polygynous, promiscuous, polyandrous, and parthenogenetic birds are 50, 30, 10, 5, and 5, respectively.

• Each polygynous or polyandrous bird has three interested mates.

• The mutation control factor $\mathrm{mcf}$ is 0.9, and the other mutation control factor used in parthenogenesis $\mathrm{mc}{f}_p$ is set as an increasing parameter that changes from 0.15 to 0.8 linearly.

• The value of parameter $\mu$ in equation (18) is set as $9.0\times {10}^{-3}$ .

• The time weighting factor w decreases linearly from 1.9 to 0.1.

• The maximum evolution step is 200.

A simply supported beam

A simply supported beam made of aluminum with 12 Euler–Bernoulli elements shown in Figure 5 is the first simulation example for damage identification. With length $l=10\mathrm{m}$ and cross-sectional width $b=0.6\mathrm{m}$ and height $h=0.4\mathrm{m}$ , the beam is discretized into 12 Euler–Bernoulli elements, Young’s modulus is $E=6.9\mathrm{GPa}$ and mass density $\rho =2700\mathrm{kg}/{\mathrm{m}}^3$ . There are three designed damages described in Table 1, elements 3, 8, and 9 with Young’s modulus reduction by 10%, 20%, and 15%, mass reduction by 5%, 10%, and 5%, respectively. The damping model is Rayleigh damping and the two Rayleigh coefficients are both assumed to be 0.01.

Figure 5.

A simply supported beam.

Table 1.

Damage status of the simply supported beam.

Damage degree (%)	Element 3	Element 8	Element 9
Stiffness	10	20	15
Mass	5	10	5

To ensure practicality and reliability of the identification, 3% random noise presented by equation (19) is added to the displacement values during localization

$$f_m=f[1+\mu (2·\mathrm{rand}-1)]$$

(19)

where $\mu$ is the noise level. To localize the damages precisely, two sensors are set in nodes 4 and 8 to get relatively larger DDILD values, to overcome the randomness of noise, the procedure is processed three times to reduce the affection of noise. The DILD and DDILD are shown in Figure 6. For the DDILD in Figure 6(b), nodes 3, 4, 8, 9, and 10 have stable changes in the chart of sensor set in node 4, and nodes 8, 9, and 10 have constant trend in the chart of the sensor in node 8. Hence, the corresponding elements 3, 8, and 9 are selected as suspected elements. The weighting factors $w_{\varphi }$ of the accMAC are all set as 1.

Figure 6.

Localization of the simply supported beam: (a) DILD and (b) difference quotient of DILD.

Case study 1: damage identification from noise-free measurements

In this numerical example, the objective function accMAC in equation (15) is used in BMO. The acceleration responses of the structure under a constant sinusoidal force acting at node 7 in the global $z$ direction $P(t)={10}^3\sin (10\pi t)\mathrm{N}$ is set as benchmark in objective function. The Newmark-beta integration method is used in calculating the dynamic response, the time increment is 0.01 s, and the time duration for the response calculation is 6.0 s and 600 time steps in total. The arrangements of acceleration measurements are two nodes locating at nodes 5 and 8.

The identification results are shown in Table 2 and Figure 7, the logarithmic best fitness value and the corresponding damage indices of the suspicious elements in a single run during the evolution process are shown in Figure 8. The six damage parameters in identification yield good results, and the maximum relative error is 0.9667% by elemental stiffness of the third element. It can be seen from Figure 8 that the fitness value converges after 80 generations, and after 70 generations, every damage index tends to convergence, the surplus iterations guarantee better identification results. The study case illustrates the availability and efficiency of the proposed method.

Table 2.

Identification results of the simply supported beam in multiple noise levels.

Results of stiffness
Noise level	${\alpha }_3$			${\alpha }_8$			${\alpha }_9$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.9	0.9087	0.9667	0.8	0.7983	0.2125	0.85	0.8497	0.0353
5%	0.9	0.9447	4.9667	0.8	0.8204	2.55	0.85	0.8351	1.7529
10%	0.9	0.9069	0.7667	0.8	0.8229	2.8625	0.85	0.8277	2.6235
Results of mass
Noise level	${\beta }_3$			${\beta }_8$			${\beta }_9$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.95	0.9559	0.6211	0.9	0.8993	0.0778	0.95	0.9461	0.4105
5%	0.95	0.9529	0.3053	0.9	0.9137	1.5222	0.95	0.9746	2.5894
10%	0.95	0.9999	5.2526	0.9	0.9068	0.7556	0.95	0.9275	2.3684

Figure 7.

Damage identification results (noise free): (a) stiffness indices and (b) mass indices.

Figure 8.

Evolution of the identification: (a) logarithmic best fitness and (b) damage indices.

Case study 2: damage identification in multiple noise levels

To illustrate the robustness and accuracy of the proposed method, damage identification is conducted under different noise levels. For acceleration responses $\mathbf{a}$ , a Gaussian distributed random noise ${\mathbf{N}}_{\mathrm{oise}}$ with zero mean and unit standard deviation in time history is added as equation (20), where $\mathrm{std}({\mathbf{a}}_{\mathrm{cal}})$ stands for the standard deviation of the calculated acceleration responses; in this study, noise levels ${\mathbf{E}}_{\mathrm{cal}}$ are set in 5% and 10%

$$\mathbf{a}={\mathbf{a}}_{\mathrm{cal}}+E_p·{\mathbf{N}}_{\mathrm{oise}}·\mathrm{std}({\mathbf{a}}_{\mathrm{cal}})$$

(20)

The identified results are shown in Table 2, with comparison to the result of noise free in Figure 9. In low noise level, the result is still satisfactory with maximum deviation less than 5%, and in high noise level the maximum relative error is 5.2526% which is still satisfactory compared to the result of 5% noise.

Figure 9.

Comparison of identification results under different noise levels: (a) stiffness indices and (b) mass indices.

In the best fitness value chart shown in Figure 10, it can be observed that with higher noise, the searching process reaches convergence faster; in the other words, the data were contaminated greater by the noise, thus the precision dropped and tend to converge easier.

Figure 10.

Comparison of logarithmic best fitness under different noise levels.

The results of the objective function based on conventional MAC are further demonstrated. The objective function that comprises MAC is shown in equation (21), where $\mathrm{NM}$ is the number of mode selected in assessment, and $w_{\varphi }$ is the weighting factor of the corresponding data.

Here, the weighting factors $w_{\varphi j}$ are all set to 1 and the first six orders of incomplete mode data at nodes 3, 5, 8, and 10 are used in the identification. The identification results are illustrated in Table 3 and Figure 11, the result without noise is very good; in two noise conditions, the maximum relative errors are similar to the results based on accMAC, but the results yield several deviations and the accuracy drops especially in the mass identification compared to the former one

$$f=\sum _{j=1}^{\mathrm{NM}}{w}_{\varphi j}^2[1-\frac{{({\left\{{\varphi }_j^C\right\}}^T•\left\{{\varphi }_j^M\right\})}^2}{{\left|\left\{{\varphi }_j^C\right\}\right|}^2{\left|\left\{{\varphi }_j^M\right\}\right|}^2}]$$

(21)

Table 3.

Identification results of the simply supported beam using MAC-based objective function.

Results of stiffness
Noise level	${\alpha }_3$			${\alpha }_8$			${\alpha }_9$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.9	0.9002	0.0222	0.8	0.7993	0.0875	0.85	0.8498	0.02353
5%	0.9	0.9299	3.3222	0.8	0.8175	2.1875	0.85	0.8284	2.5412
10%	0.9	0.9346	3.8444	0.8	0.8215	2.6875	0.85	0.8318	2.1412
Results of mass
Noise level	${\beta }_3$			${\beta }_8$			${\beta }_9$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.95	0.9497	0.0316	0.9	0.8994	0.0667	0.95	0.9494	0.0636
5%	0.95	0.9408	0.9684	0.9	0.8572	4.7556	0.95	0.9285	2.2632
10%	0.95	1	5.2632	0.9	0.8608	4.3556	0.95	0.9317	1.926

Figure 11.

Comparison of identification results between accMAC-based and MAC-based objective functions under different noise levels: (a) stiffness indices and (b) mass indices.

The damage model combines both stiffness and mass damages. More information is needed in the objective function should better results are to be acquired; the traditional MAC takes more measure points and it is hard to measure mode shape especially in high order. On the contrary, the proposed method can measure data from many time steps which provides more damage information.

A four-span beam

A four-span beam structure as shown in Figure 12 is the second numerical example in the study. The length of the beam is 6 m, with width $b=0.2\mathrm{m}$ and height $h=0.3\mathrm{m}$ , Young’s modulus is $E=20.5\mathrm{GPa}$ and mass density $\rho =7800\mathrm{kg}/{\mathrm{m}}^3$ . The beam is divided into 40 Euler–Bernoulli elements and the damage information is in Table 4, and the two Rayleigh damping coefficients are both 0.01.

Figure 12.

Four-span beam structure.

Table 4.

Damage status of the four-span beam.

Damage degree (%)	Element 3	Element 12	Element 25	Element 33
Stiffness	8	12	10	15
Mass	5	10	5	10

The deflection data of two sensors located at the midpoint of first and third spans, namely, nodes 6 and 26, are used in localization, and all the influence line data were contaminated by 3% random noise. As shown in Figure 13, nodes 3, 4, 12, and 13 show stable non-zero values in the chart of DDILD measured in node 6, and the bars of nodes 25, 26 and 33, 34 also show the potential of damage and is certified by the chart of DDILD measure in node 26. Hence, the corresponding elements 3, 12, 25, and 33 are the suspected damage elements of the four-span beam.

Figure 13.

Localization of the four-span beam: (a) DILD and (b) difference quotient of DILD.

Case study 3: damage identification in multiple noise levels

The acceleration responses were obtained by a sinusoidal force acting at the 16th node of the beam in the global z direction with $P(t)=200\sin (20\pi t)\mathrm{N}$ in calculating the dynamic response, the time increment is 0.01 s, and the time duration for the response calculation is 8.0 s. In this study case, the responses of two nodes locating at 5 and 25 are used in damage identification. The weighting factors $w_{\varphi }$ of the accMAC are all set as 1.

It takes almost 140 steps for the noise-free identification to reach a stable convergence; with the increase of noise level, the converge step decreases to 30 and 20 in 5% and 10% noise level, respectively. It can be implied from the convergence chart (Figure 14) that the contamination of noise influences the details of the data; thus, the algorithm tends to converge once an approximate optimal solution is found, though it converges faster, the accuracy drops.

Figure 14.

Comparison of logarithmic best fitness under different noise levels.

The final identification results of three noise levels are shown in Table 5 and Figure 15. The noise-free condition yields a good result with maximum relative error less than 0.5%; even in the 5% noise condition, only the mass damage index ${\beta }_{12}$ of the second suspicious element reaches the relative error of 2.2%; in 10% noise level, the data were greatly contaminated by the noise, both the stiffness and mass damage indices of the element 33 have relatively large identification errors 5.1294% and 7.9222%; the results of others elements are still good.

Table 5.

Identification results of the four-span beam in multiple noise levels.

Results of stiffness
Noise level	${\alpha }_3$			${\alpha }_{12}$			${\alpha }_{25}$			${\alpha }_{33}$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.92	0.9207	0.0761	0.88	0.8789	0.125	0.9	0.9008	0.0889	0.85	0.8446	0.6353
5%	0.92	0.914	0.6522	0.88	0.874	0.6818	0.9	0.8888	1.2444	0.85	0.8473	0.3176
10%	0.92	0.9363	1.7717	0.88	0.8623	2.0114	0.9	0.903	0.3333	0.85	0.8064	5.1294
Results of mass
Noise level	${\beta }_3$			${\beta }_{12}$			${\beta }_{25}$			${\beta }_{33}$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.95	0.9511	0.1158	0.9	0.8983	0.1889	0.95	0.9499	0.0105	0.9	0.8957	0.4778
5%	0.95	0.9474	0.2737	0.9	0.8802	2.2	0.95	0.9409	0.9579	0.9	0.8897	1.1444
10%	0.95	0.9724	2.3579	0.9	0.8796	2.2667	0.95	0.9583	0.8737	0.9	0.8287	7.9222

Figure 15.

Comparison of identification results under different noise levels: (a) stiffness indices and (b) mass indices.

An objective function originated from dynamic responses is further investigated. Here, the dynamic displacement of the same structure is acting as objective function, hence equation (22)

$$f=\sum _{j=1}^{\mathrm{NA}}{w}_{\mathrm{dj}}^2[1-\frac{{({\left\{{\mathbf{dis}}_j^C\right\}}^T•\left\{{\mathbf{dis}}_j^M\right\})}^2}{{\left|\left\{{\mathbf{dis}}_j^C\right\}\right|}^2{\left|\left\{{\mathbf{dis}}_j^M\right\}\right|}^2}]$$

(22)

where the superscripts $C$ and $M$ indicate the data from calculation and measurement separately, and the weighting factors $w_{\mathrm{dj}}$ are all set to 1. The displacement sensors are placed in two nodes locating at 5 and 25, too.

The results are shown in Figure 16. The results from noise-free data are satisfactory but in 5% noise condition there exists huge error between the identified results and the designed ones, and the identification of 10% is not shown here. The displacement cannot provide more dynamic information compared to acceleration, and the disturbance of noise makes the data even worse, thus in the noise condition it failed to identify damages. The comparison shows the superiority of the chosen data type for structural damage detection.

Figure 16.

Identification results of disMAC-based objective functions under different noise levels: (a) stiffness indices and (b) mass indices.

Identified results of MAC are also studied here, as shown in equation (21). The weighting factors $w_{\varphi j}$ are all set to 1, the first six orders of mode data at midpoint of each span, namely, nodes 5, 15, 25, and 35, are used in the identification. The identification results are illustrated in Table 6 and Figure 17, the result without noise is quite accurate, while in noise conditions there yields large relative errors. Unlike in simply supported beam, the maximum error reaches 14.1224% in 5% noise level, and 26.2942% in 10% condition, which are far worse than the corresponding results yielded by accMAC. The number of data in MAC is much less than from accMAC, which suggests that results from MAC are susceptible to noise, and the randomness from noises causes huge errors. In general, results from accMAC are better than ones from MAC, especially for rather complex structures in noise measurement conditions.

Table 6.

Identification results of the four-span beam using MAC-based objective function.

Results of stiffness
Noise level	${\alpha }_3$			${\alpha }_{12}$			${\alpha }_{25}$			${\alpha }_{33}$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.92	0.9195	0.0522	0.88	0.8796	0.05	0.9	0.9	0	0.85	0.8498	0.0224
5%	0.92	0.897	2.5033	0.88	0.8269	6.0341	0.9	0.8988	0.1367	0.85	0.7299	14.1224
10%	0.92	0.8636	6.1355	0.88	0.8206	6.751	0.9	0.9348	3.8614	0.85	0.8606	1.2528
Results of mass
Noise level	${\beta }_3$			${\beta }_{12}$			${\beta }_{25}$			${\beta }_{33}$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.95	0.949	0.1011	0.9	0.8996	0.04	0.95	0.9499	0.0158	0.9	0.8996	0.0489
5%	0.95	0.9564	0.6747	0.9	0.8109	9.8978	0.95	1	5.2632	0.9	0.7849	12.7859
10%	0.95	0.9442	0.6134	0.9	0.6634	26.2942	0.95	1	5.2631	0.9	0.9653	7.2596

Figure 17.

Identification results of MAC-based objective functions under different noise levels: (a) stiffness indices and (b) mass indices.

Case study 4: damage identification under impulse force

In the previous study cases, dynamic responses of sinusoidal excitation show good result in damage identification. This study case will further conduct damage identification using acceleration response from impulse excitation. An impulse is supposed to act at the 16th node of the beam in the global z direction with $\{\begin{array}{c}F(t)={10}^{4}t\mathrm{N}(0<t\le 0.03\mathrm{s})\hfill \\ \multicolumn{1}{c}{F(t)={10}^4(0.06-t)\mathrm{N}(0.03\mathrm{s}<t\le 0.06\mathrm{s})}\end{array}$ , the time increment is 0.001 s, and measurement lapsed for 1 s. Similarly, the acceleration responses of two nodes, 5 and 25, are used in damage identification.

Figure 18 and Table 7 show the results; one can find that the identification error becomes large especially in noise condition compared to results from sinusoidal load, the largest relative error is 2.9882% and 11.5529% in 5% and 10% noise level condition, respectively, but the identified results are still satisfactory. The time histories of the calculated acceleration and the simulated measured acceleration of 25th node are shown by Figure 19. As the calculated damage parameters converged to the true ones, the calculated dynamic responses match well with measured data, and it is obvious that the lower the noise level, the better the calculated acceleration data accord with the true. The reason is that the noise becomes significant when the acceleration responses ebb. This case study shows that the acceleration responses under the impulse excitation can also be used for damage identification by proposed method, though the results are worse than identification from sinusoidal force.

Figure 18.

Comparison of identification results using impulse excitation under different noise levels: (a) stiffness indices and (b) mass indices.

Table 7.

Identification results of the four-span beam using impulse excitation.

Results of stiffness
Noise level	${\alpha }_3$			${\alpha }_{12}$			${\alpha }_{25}$			${\alpha }_{33}$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.92	0.9159	0.4457	0.88	0.8725	0.8523	0.9	0.9001	0.0111	0.85	0.8504	0.0471
5%	0.92	0.9266	0.7174	0.88	0.8809	0.1023	0.9	0.897	0.3333	0.85	0.8246	2.988
10%	0.92	0.9039	1.75	0.88	0.8246	6.2955	0.9	0.8701	3.3222	0.85	0.9482	11.5529
Results of mass
Noise level	${\beta }_3$			${\beta }_{12}$			${\beta }_{25}$			${\beta }_{33}$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.95	0.9416	0.8842	0.9	0.8983	0.1889	0.95	0.9511	0.1158	0.9	0.899	0.1111
5%	0.95	0.9625	1.3158	0.9	0.8932	0.7556	0.95	0.9466	0.3579	0.9	0.8784	2.4
10%	0.95	0.9265	2.4737	0.9	0.8327	7.4778	0.95	0.9465	0.3684	0.9	0.9668	7.4222

Figure 19.

Comparison of the calculated acceleration responses and the measured responses.

Case study 5: effect of additional measure points to reduce the impact of noise

It can conclude from the above cases that in the conditions with high noise level, the accuracy of the damage identification dropped to some extent. Considering that data of only two nodes were used in previous cases, the result of more measured points is investigated in the study case.

The model of four-span beam with sinusoidal excitation $P(t)=200\sin (20\pi t)\mathrm{N}$ is re-examined. Two additional sensors are placed at nodes 15 and 35 as shown in scenario 2 in Table 8 so that there is a sensor distributed at each span at 5th, 15th, 25th, and 35th nodes. Results of three noise level at 0, 5%, and 10% are inspected.

Table 8.

Arrangement of the measure points.

Measurement scenarios	Number of measure points	Distribution of measure points (node)
1	2	5, 25
2	4	5, 15, 25, 35

Table 9 shows the identified results, compared with case 3; the results in noise conditions show better accuracy when the acceleration data of only two nodes are used, as shown in Figure 20. The maximum relative error dropped from 2.2% to 1.6235% in 5% noise condition; in high noise scenario with 10% noise, relative error dropped from 7.9222% to 4.7368%. The identification yields better results especially the mass damage indices. This case illustrates that with proper addition of the sensors along the beam, damage identification through the proposed method can yield better results in noise condition.

Table 9.

Identification results of the four-span beam with additional measure points.

Results of stiffness
Noise level	${\alpha }_3$			${\alpha }_{12}$			${\alpha }_{25}$			${\alpha }_{33}$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.92	0.9199	0.01087	0.88	0.8798	0.0227	0.9	0.9	0	0.85	0.85	0
5%	0.92	0.9248	0.5217	0.88	0.8809	0.1023	0.9	0.9033	0.3667	0.85	0.8638	1.6235
10%	0.92	0.8903	3.2283	0.88	0.8446	4.0227	0.9	0.9028	0.3111	0.85	0.8429	0.8353
Results of mass
Noise level	${\beta }_3$			${\beta }_{12}$			${\beta }_{25}$			${\beta }_{33}$
	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)	Real	Identified	Relative error (%)
0	0.95	0.9498	0.02105	0.9	0.9	0	0.95	0.9501	0.0105	0.9	0.9001	0.0111
5%	0.95	0.9541	0.4316	0.9	0.904	0.4444	0.95	0.9502	0.0211	0.9	0.9195	2.1667
10%	0.95	0.905	4.7368	0.9	0.8743	2.8556	0.95	0.956	0.6316	0.9	0.8927	0.8111

Figure 20.

Comparison of identification results using different measure arrangements under different noise levels: (a) stiffness indices and (b) mass indices.

Discussion

Through the numerical examples conducted in the simulation section, the proposed method presents several advantages in structural damage identification. However, there are several aspects of the proposed method cannot be covered in the studies above. Hence, some discussions are made as follows:

• In the study, damage model is parameterized in a simple way, while the mechanisms of damages in real structures are much more complicated. Therefore, complicated damage model, such as breathing crack or even crack with geometric parameters, should be considered to examine the proposed method in the future study.

• It requires both data from time domain and frequency domain in the proposed method. This will increase workload and may reduce the accuracy and efficiency of identification, and the measured data are always contaminated by noise which means there are noises exist in both deflection and acceleration data. Improvements should focus on identification derived from a single type of data using a newly constructed quantity ³⁹ or new method.

• It should consider conducting experiments in the study for the proposed method in further research, which may validate the practicality on real structures.

• Excitations are known in the identification process. However, the most common type of vibration in real structures is random vibration excited by ambient excitations. Further investigation should aim to apply the proposed method on structures under ambient excitation.

Conclusion

By making use of the influence line difference and acceleration data, a two-step damage detection method for beam structures is proposed in this study. A damage model combined with stiffness and mass reduction is applied. The influence line data of displacement is used and further investigated to localize the damage. An objective function based on acceleration and the form of MAC is established, and then BMO is used to minimize the discrepancies between the simulated measured data and the data from damaged structure. Two beam structures with several study cases reflect the effectiveness of the proposed method; comparison shows that the proposed function has better results than identification dynamic displacement or MAC. In addition, both the acceleration data from sinusoidal and impulse excitation proved successful in the identification. All numerical simulations are conducted in three different noise levels which demonstrate that the proposed method is insensitive to measurement noise, and by increasing the sampling number of acceleration measure points, the results in noise condition had improved.