Sensitivity analyses of resonant frequencies and modal strain energy of damaged beams by perturbation method

1. Introduction

Early detection and quantification of structural damages are of great importance to ensure the lifetime safety, prevent catastrophic events, and increase the service of structures. Currently available damage detection methods include magnetic field methods, ultrasonic methods, eddy-current methods, and thermal field methods [1], which are localized experimental methods. Most of these experimental methods require that the damage location is approximately known in advance and accessible for inspection. Due to these limitations, the need for global damage detection methods has led to the development of the static damage detection technology [2, 3] and dynamic damage detection technology [4]. The underlying idea behind the vibration-based technology is that, the changes in the physical properties (stiffness, mass, and damping) of a structure will cause detectable changes in its modal parameters (frequencies, modes, and strain energy) which are the functions of its physical properties [5, 6]. Hence, the changes in the modal parameters can be assumed as indictors for damage detection. The variations of the physical properties of a structure, caused by damages, include the reduction in the structural stiffness and modification of the structural geometry, etc. Vibration-based damage detection techniques have attracted much attention in past decades and many approaches have been developed, based on the vibration response of structures [7]. Early studies found that the resonant frequencies were less sensitive to low levels of damage [1, 8]. The frequencies, as global parameters, generally cannot provide spatial information about structural changes [1]. On the other hand, it has been experimentally proved that the mode shapes [7-9] and modal strain energy [10-13] are more sensitive to structural damages.

The damage detection methods using mode shapes and their derivatives as a basic feature are more sensitive to local damage as the mode shapes contain local information [7-9]. The traditional mode shape change method uses the mode shape changes between the intact and damaged structures as a basic feature for damage identification. The displacement mode shape itself is not very sensitive to a small damage. To enhance the sensitivity of mode shape data to the damage, the mode shape curvature has been proposed as a promising feature for damage identification [7, 8]. For beam-type or plate-type structures, the second derivatives of the mode shapes are highly sensitive to a damage and can be used to localize it [9]. The difference of the curvature mode shapes between the intact and damaged structures is introduced as a damage indicator [7]. The modal strain energy is directly related to the mode shape curvatures for beam-type or plate-type structures; as a special case of the mode shape curvature-based method, the modal strain energy has been widely used for damage identification [10-16]. A damage indictor was created [13, 17] based on the change in the modal strain energy between the intact and damaged structures. The selection of a damage indictor will affect the final results and accuracy of structural damage detection [4]. The effectiveness of a damage detection method is determined by its sensitivity to a structural damage.

The sensitivity of the modal strain energy-based methods is related to those of the eigen-pairs of the concerned system. A variety of methods for computing the derivatives of eigen-pairs have been developed by many researchers. When a damage occurs in a structure, it can be represented by a small perturbation of the original system [18, 19]. The variation of an eigen-function can be expressed as a linear combination of the eigen-functions of the undamaged structure [18-20]. This method employs the modal superposition idea and needs all the mode shapes to obtain an accurate sensitivity of each mode shape, which is almost an impossible task. Often only the lower order frequencies and associated mode shapes are calculated. Recently, the calculation of the derivative of an eigenvector was substantially simplified, it required only the information of the corresponding eigenvector to be differentiated [21]. Most of sensitivity analyses were based on the matrix theory of a discrete system. In matrix analyses, it is not easy to reveal the mechanism how the damage parameter influences the sensitivity of damage indictors [19], and the influential relationship are submerged by the matrix calculations. To the knowledge of the authors, there are not theoretical comparisons between the sensitivities of various damage detection methods based on a continuum model.

To introduce a damage indictor for structural damage detection, one essential issue is to effectively evaluate its sensitivity to the variation of the structural parameters. This paper aims to propose a perturbation method based on the Rayleigh quotient via a beam model to evaluate the sensitivities of the first frequency and the modal strain energy in a region. The perturbation of the eigen-function of the damaged beam is expressed as the Taylor series of the eigen-functions of the undamaged beam, and the sensitivities are obtained. The analytical relationships are established between the variation rates of the frequency and strain energy with the damage parameter. These relationships approximately illustrate the sensitivity of these damage indictors to a structural damage. The sensitivities of the frequency-based and modal strain energy-based methods are investigated and compared. Furthermore, the finite element (FE) simulations and experimental modal analyses are conducted to validate the results obtained in the perturbation analyses. It is expected that the proposed method would be able to approximately assess the quality of a damage indictor.

2. Method

Consider the bending vibration of an elastic beam (Fig. 1). The Young’s modulus and mass density of the beam are denoted by $E$ and $\rho$; the area and second moment of a cross section are denoted by $A$ and $I$. A damage zone with a length of $d$ exists at $z=s$. Assume that the length ($L$) of the beam is more than 10 times the characteristic length of its cross section, and it can be analyzed by the Euler beam theory [22].

2.1. Perturbation analysis of a beam

For the bending in a principal plane, the governing equation of free vibration of the beam is expressed in the variational form [23] as:

1

$\mathrm{\Pi }-\lambda \mathrm{\Psi }=0,$

where $\mathrm{\Pi }={\int }_0^{L}EI{w}^{\text{'}\text{'}2}dz$ is (2 times) the strain energy and $\mathrm{\Psi }={\int }_0^L\rho A{w}^{2}dz$ is (2 times) the amplitude of kinetic energy, and where $\lambda ={\omega }^2$ is the eigen-value with $\omega$ being the (angular) frequency. Define ${\mathrm{\Pi }}_{ij}={\int }_0^{L}EI\cdot w_i^{\text{'}\text{'}}\cdot w_j^{\text{'}\text{'}}dz$ and ${\mathrm{\Psi }}_{ij}={\int }_0^L\rho A\cdot w_i\cdot w_{j}dz$, where, ${\lambda }_i$, ${\lambda }_j$ and $w_i$, $w_j$ are the $i$th and $j$th eigen-values and eigen-functions of the undamaged beam, which satisfy the orthogonal conditions:

2

${\mathrm{\Psi }}_{ij}={\delta }_{ij},{\mathrm{\Pi }}_{ij}={\lambda }_i{\mathrm{\Psi }}_{ij},$

where ${\delta }_{ij}$ is the Kronecker delta. It should be noted that ${\mathrm{\Pi }}_{ij}$ and ${\mathrm{\Psi }}_{ij}$ are not tensors and a repeated index does not represent summation. The $i$th eigen-value ${\lambda }_i$ satisfies the Rayleigh quotient:

3

${\lambda }_i=\frac{{\mathrm{\Pi }}_{ii}}{{\mathrm{\Psi }}_{ii}}.$

Fig. 1Bending vibration of an elastic beam: a) location of a damage zone at the position z=s and the damage length is d, b) cut depth in the damage section

a)

b)

After a damage is introduced into the original intact beam, its $EI$ and $\rho A$ have perturbations of $\mathrm{\Delta }\left(EI\right)$and $\mathrm{\Delta }\left(\rho A\right)$ respectively; hence, for the damaged beam:

4

$\tilde{E}\tilde{I}=EI+\mathrm{\Delta }\left(EI\right),\tilde{\rho }\tilde{A}=\rho A+\mathrm{\Delta }\left(\rho A\right).$

The governing equation of the damaged beam in the variational form is:

5

$\tilde{\mathrm{\Pi }}-\tilde{\lambda }\tilde{\mathrm{\Psi }}=0,$

with $\tilde{\mathrm{\Pi }}={\int }_0^L\tilde{E}\tilde{I}{\tilde{w}}^{\text{'}\text{'}2}dz$ and $\tilde{\mathrm{\Psi }}={\int }_0^L\tilde{\rho }\tilde{A}{\tilde{w}}^{2}dz$. Similarly, we define the integrations:

6

${\tilde{\mathrm{\Pi }}}_{ij}={\int }_0^L\tilde{E}\tilde{I}\cdot {\tilde{w}}_i^{\text{'}\text{'}}\cdot {\tilde{w}}_j^{\text{'}\text{'}}dz,{\tilde{\mathrm{\Psi }}}_{ij}={\int }_0^L\tilde{\rho }\tilde{A}\cdot {\tilde{w}}_i\cdot {\tilde{w}}_{j}dz,$

where, ${\tilde{w}}_i$ and ${\tilde{w}}_j$ are the $i$th and $j$th eigen-functions of the damaged beam, and the corresponding eigen-values are ${\tilde{\lambda }}_i$ and ${\tilde{\lambda }}_j$. They satisfy the orthogonal conditions:

7

${\tilde{\mathrm{\Psi }}}_{ij}={\delta }_{ij},{\tilde{\mathrm{\Pi }}}_{ij}={\tilde{\lambda }}_i{\tilde{\mathrm{\Psi }}}_{ij}.$

It is assumed that the eigen-pair, ${\tilde{\lambda }}_i$ and ${\tilde{w}}_i$, of the damaged beam are close to those of the undamaged beam, ${\lambda }_i$ and $w_i$; hence they can be expressed as:

8

${\tilde{\lambda }}_i={\lambda }_i+∆{\lambda }_i,{\tilde{w}}_i=w_i+∆w_i,$

where $∆{\lambda }_i$ and $∆w_i$ are the perturbations of ${\lambda }_i$ and $w_i$, introduced by the damage.

By using Eqs. (4) and (8) and ignoring the terms higher than second order, Eq. (6) becomes:

9

${\tilde{\mathrm{\Pi }}}_{ij}={\mathrm{\Pi }}_{ij}+{\int }_0^L\left[∆\left(EI\right)\cdot w_i^{\text{'}\text{'}}\cdot w_j^{\text{'}\text{'}}+EI\cdot \left(w_i^{\text{'}\text{'}}\cdot ∆w_j^{\text{'}\text{'}}+w_j^{\text{'}\text{'}}\cdot ∆w_i^{\text{'}\text{'}}\right)\right]dz,$

and:

10

${\tilde{\mathrm{\Psi }}}_{ij}={\mathrm{\Psi }}_{ij}+{\int }_0^L\left[∆\left(\rho A\right)\cdot w_i\cdot w_j+\rho A\cdot \left(w_i\cdot ∆w_j+w_j\cdot ∆w_i\right)\right]dz,$

and then:

11

${{\tilde{\lambda }}_i\tilde{\mathrm{\Psi }}}_{ij}={\lambda }_i{\mathrm{\Psi }}_{ij}+∆{\lambda }_i{\mathrm{\Psi }}_{ij}+{\lambda }_i{\int }_0^L\left[∆\left(\rho A\right)\cdot w_i\cdot w_j+\rho A\cdot \left(w_i\cdot ∆w_j+w_j\cdot ∆w_i\right)\right]dz.$

Following [24] and [18], the perturbation of the $i$th eigen-function of the damaged beam can be expressed as the Taylor series of the eigen-functions of the undamaged beam:

12

$∆w_i=\sum _{n=1}^{\infty }{\zeta }_{in}{w}_n,$

where ${\zeta }_{in}$ is the coefficient of the $n$th eigen-function $w_n$. Substituting Eq. (12) into Eqs. (9) and (11) and utilizing the orthogonal conditions Eq. (2), we have:

13

${\tilde{\mathrm{\Pi }}}_{ij}={\mathrm{\Pi }}_{ij}+{\int }_0^L∆\left(EI\right)\cdot w_i^{\text{'}\text{'}}\cdot w_j^{\text{'}\text{'}}dz+{\mathrm{\Pi }}_{ii}{\zeta }_{ji}+{\mathrm{\Pi }}_{jj}{\zeta }_{ij},$

and:

14

${{\tilde{\lambda }}_i\tilde{\mathrm{\Psi }}}_{ij}={\lambda }_i{\mathrm{\Psi }}_{ij}+∆{\lambda }_i{\mathrm{\Psi }}_{ij}+{\lambda }_i{\int }_0^L∆\left(\rho A\right)\cdot w_i\cdot w_{j}dz+{\lambda }_i\left({\mathrm{\Psi }}_{ii}{\zeta }_{ji}+{\mathrm{\Psi }}_{jj}{\zeta }_{ij}\right).$

Substitution of Eqs. (13) and (14) into the second equality of Eq. (7) yields:

15

${\int }_0^L∆\left(EI\right)\cdot w_i^{\text{'}\text{'}}\cdot w_j^{\text{'}\text{'}}dz+\left({\lambda }_j-{\lambda }_i\right){\mathrm{\Psi }}_{jj}{\zeta }_{ij}-∆{\lambda }_i{\mathrm{\Psi }}_{ij}-{\lambda }_i{\int }_0^L∆\left(\rho A\right)\cdot w_i\cdot w_{j}dz=0.$

For $j=i$, Eq. (15) gives:

16

$∆{\lambda }_i=\left[{\int }_0^L∆\left(EI\right)\cdot w_i^{\text{'}\text{'}}\cdot w_i^{\text{'}\text{'}}dz-{\lambda }_i{\int }_0^L∆\left(\rho A\right)\cdot w_i\cdot w_{i}dz\right]/{\mathrm{\Psi }}_{ii}.$

By using Eqs. (3) and (16), the sensitivity of the $i$th eigen-value is:

17

$\frac{∆{\lambda }_i}{{\lambda }_i}=\left[{\int }_0^L∆\left(EI\right)\cdot w_i^{\text{'}\text{'}}\cdot w_i^{\text{'}\text{'}}dz-{\lambda }_i{\int }_0^L∆\left(\rho A\right)\cdot w_i\cdot w_{i}dz\right]/{\mathrm{\Pi }}_{ii}.$

For $j\ne i$, Eq. (15) gives:

18

${\zeta }_{ij}=\left[{\int }_0^L∆\left(EI\right)\cdot w_i^{\text{'}\text{'}}\cdot w_j^{\text{'}\text{'}}dz-{\lambda }_i{\int }_0^L∆\left(\rho A\right)\cdot w_i\cdot w_{j}dz\right]/\left[\left({\lambda }_i-{\lambda }_j\right){\mathrm{\Psi }}_{jj}\right].$

The quantity ${\zeta }_{ii}$ is found from Eq. (10) and the orthogonal conditions Eqs. (3) and (7), which leads to:

19

${\zeta }_{ii}=-\frac{\frac{1}{2}{\int }_0^L∆\left(\rho A\right)\cdot w_i\cdot w_{i}dz}{{\mathrm{\Psi }}_{ii}}.$

The modal strain energy stored in a region (say, the $k$th) of the $i$th mode is expressed as [19, 25, 26]:

20

${\mathrm{\Theta }}_{ii}^{\left(k\right)}={\int }_{Z_k}^{Z_{k+1}}EI\cdot w_i^{\mathrm{\text{'}}\mathrm{\text{'}}}\cdot w_i^{\mathrm{\text{'}}\mathrm{\text{'}}}dz,$

where [$Z_k,Z_{k+1}$] are the nodal coordinates of the $k$th region. The change of the modal strain energy (by ignoring the terms higher than second order) is:

21

$∆{\mathrm{\Theta }}_{ii}^{\left(k\right)}={\int }_{Z_k}^{Z_{k+1}}∆\left(EI\right)\cdot w_i^{\mathrm{\text{'}}\mathrm{\text{'}}}\cdot w_i^{\mathrm{\text{'}}\mathrm{\text{'}}}dz+{\int }_{Z_k}^{Z_{k+1}}2EI\cdot w_i^{\mathrm{\text{'}}\mathrm{\text{'}}}\cdot {∆w}_i^{\mathrm{\text{'}}\mathrm{\text{'}}}dz.$

Substitution of Eq. (12) into Eq. (21) yields:

22

$∆{\mathrm{\Theta }}_{ii}^{\left(k\right)}={\int }_{Z_k}^{Z_{k+1}}∆\left(EI\right)\cdot w_i^{\mathrm{\text{'}}\mathrm{\text{'}}}\cdot w_i^{\mathrm{\text{'}}\mathrm{\text{'}}}dz+2\sum _{n=1}^{\infty }{\zeta }_{in}{\int }_{Z_k}^{Z_{k+1}}EI\cdot w_i^{\mathrm{\text{'}}\mathrm{\text{'}}}\cdot w_n^{\mathrm{\text{'}}\mathrm{\text{'}}}d\mathrm{z}.$

The beam model (Fig. 1) has a length $L=$ 0.80 m and a rectangular cross section with a width $b=$ 0.03 m and a height $h=$ 0.06 m. The Young’s modulus, Poisson’s ratio, and mass density of the beam (steel) are $E=$ 210 GPa, $v=$ 0.3, and $\rho =$ 7800 kg/m³. Two damage levels (10 % and 30 %) are introduced into the damage zone ($d=$ 0.02 m), i.e., the portion of $\delta h$ is removed (here $\delta =$ 10 % and 30 % respectively) along the bottom of the damage zone (Fig. 2).

Table 1Reduction ratios of the area and its second moment induced by a cut in a cross section

$h_d$ / $h$	$A_d$ / $A$	$I_{dx}$ / $I_x$	$I_{dy}$ / $I_y$	$\delta$	$\alpha$	${\varphi }_x$	${\varphi }_y$
0.9	0.9	0.729	0.9	0.1	–0.1	–0.271	–0.1
0.7	0.7	0.343	0.7	0.3	–0.3	–0.657	–0.3

As the damage is localized in a small zone where the cross sectional area, $A$, of the beam and its second moment, $I$, are reduced; thus, $∆\left(\rho A\right)=\alpha \left(\rho A\right)$ and $∆\left(EI\right)=\phi \left(EI\right)$ in the damage zone, everywhere else $∆\left(\rho A\right)=$ 0 and $∆\left(EI\right)=$ 0. The reduction ratios of the cross sectional area and flexible rigidity ($\alpha$and $\phi$) are related to the parameter $\delta$ as listed in Table 1. As the second moment of a cross section area is different about $x$ and $y$-axis, the reduction ratio of the second moment is denoted by ${\phi }_x$ and ${\phi }_y$ for bending around $x$ and $y$-axis respectively. In the table, $h_d$ and $A_d$ are the height and area of the cross section in the damage zone and $I_{dx}$ and $I_{dy}$ are the corresponding second moment of the cross section $A_d$ about $x$ and $y$-axis; $I_x$ and $I_y$ are the second moment of the cross section $A$ about $x$ and $y$-axis. It is easy to obtain that: $\alpha =-\delta$, ${\phi }_y=-\mathrm{}\delta$, and ${\phi }_x={\left(1+\delta \right)}^3-1$.

For the beam concerned, Eq. (16) gives:

23

$∆{\lambda }_i=\left(\varphi {\int }_{s-d/2}^{s+d/2}EI\cdot w_i^{\text{'}\text{'}}\cdot w_i^{\text{'}\text{'}}dz-{\lambda }_i\alpha {\int }_{s-d/2}^{s+d/2}\rho A\cdot w_i\cdot w_{i}dz\right)/{\mathrm{\Psi }}_{ii}.$

From Eq. (17), the sensitivity of the $i$th eigen-value is simplified to:

24

$\frac{∆{\lambda }_i}{{\lambda }_i}=\left(\varphi {\int }_{s-d/2}^{s+d/2}EI\cdot w_i^{\text{'}\text{'}}\cdot w_i^{\text{'}\text{'}}dz-{\lambda }_i\alpha {\int }_{s-d/2}^{s+d/2}\rho A\cdot w_i\cdot w_{i}dz\right)/{\int }_0^{L}EI\cdot w_i^{\text{'}\text{'}}\cdot w_i^{\text{'}\text{'}}dz.$

Eqs. (18) and (19) give:

25

${\zeta }_{ij}=\left(\varphi {\int }_{s-d/2}^{s+d/2}EI\cdot w_i^{\text{'}\text{'}}\cdot w_j^{\text{'}\text{'}}dz-{\lambda }_i\alpha {\int }_{s-d/2}^{s+d/2}\rho A\cdot w_i\cdot w_{j}dz\right)/\left[\left({\lambda }_i-{\lambda }_j\right){\mathrm{\Psi }}_{jj}\right],j\ne i,$

and:

26

${\zeta }_{ii}=-\frac{\frac{1}{2}\alpha {\int }_{s-d/2}^{s+d/2}\rho A\cdot w_i\cdot w_{i}dz}{{\mathrm{\Psi }}_{ii}}.$

For the case that the $k$th region falls within the damage zone [$s-d/2,s+d/2$], i.e., $∆\left(EI\right)=\varphi \left(EI\right)$, from Eqs. (22) and (20), the change rate of the modal strain energy for the $i$th mode can be expressed as:

27

$\frac{∆{\mathrm{\Theta }}_{ii}^{\left(k\right)}}{{\mathrm{\Theta }}_{ii}^{\left(k\right)}}=\varphi +2\sum _{n=1}^{\infty }{\zeta }_{in}{\int }_{Z_k}^{Z_{k+1}}EI\cdot w_i^{\text{'}\text{'}}\cdot w_n^{\text{'}\text{'}}dz/{\int }_{Z_k}^{Z_{k+1}}EI\cdot w_i^{\text{'}\text{'}}\cdot w_i^{\text{'}\text{'}}dz.$

If the region concerned is outside the damage zone, $∆\left(EI\right)=$ 0, then Eqs. (22) and (20) lead to:

28

$\frac{∆{\mathrm{\Theta }}_{ii}^{\left(k\right)}}{{\mathrm{\Theta }}_{ii}^{\left(k\right)}}=2\sum _{n=1}^{\infty }{\zeta }_{in}{\int }_{Z_k}^{Z_{k+1}}EI\cdot w_i^{\mathrm{\text{'}}\mathrm{\text{'}}}\cdot w_n^{\mathrm{\text{'}}\mathrm{\text{'}}}dz/{\int }_{Z_k}^{Z_{k+1}}EI\cdot w_i^{\mathrm{\text{'}}\mathrm{\text{'}}}\cdot w_i^{\mathrm{\text{'}}\mathrm{\text{'}}}dz.$

The beam is divided into 80 regions, and the modal strain energy for the first modes and its change rates are calculated for the two damage levels. The integrations are carried out in a MATLAB program. The eigen-pairs of an intact beam under free-free boundary conditions are listed in the Appendix.

2.2. FE analyses

For the sake of validations, the FE analyses of the steel beam (Fig. 2) are conducted. Besides the intact beam model, 2 damage levels (10 % and 30 %) are considered in the damage zone ($d=$ 0.02 m). The software package ABAQUS (SIMULIA Inc., Providence, USA) is used to obtain the first 8 resonant frequencies, mode shapes, and the corresponding strain energy density in each element. Both ends are set free, and about 12400 continuum elements (C3D8R) are used.

Fig. 2FE mesh of the steel beam

2.3. Experimental modal analyses

The experimental set-up (Fig. 3) for modal analysis includes an excitation hammer, response measuring sensors (accelerometers) and an acquisition system (JM3840, Jing-Ming Technology Inc., Yangzhou, China). The specimen is laid on sponge (Fig. 3) to simulate free-free boundary condition. Totally 10 accelerometers with nominal sensitivity of 100 mV/g are attached on the top surface of the beam. The impulse force is provided by the micro-hammer instrumented with a force sensor to measure the input force. The time histories of the force and acceleration are collected by the acquisition instrument and analyzed in a computer. The frequency response functions (FRFs) and the modal parameters are gained by using the kit software associated with JM3840 acquisition instrument. The parameters of the experimental instruments are listed in Table 2.

Fig. 3Experimental set-up for modal analyses and testing specimens

a)

b)

Table 2Parameters of the experimental instruments

Name	Pattern	Specifications	Manufacturer
Excitation hammer/sensor	JML-03	Sensibility: 0.08 mv/N	Jing-Ming Technology Inc., Yangzhou, China
Accelerometer	YD2150	Sensibility: 4.99 mv/m/s2	Xi-Yuan Technology Inc., Yangzhou, China
Data acquisition system	JM3840	16 channels; sampling rate: 2 kHz	Jing-Ming Technology Inc., Yangzhou, China
Modal analysis software	JMTEST	Version 7.9	Jing-Ming Technology Inc., Yangzhou, China

Calculation of the strain energy in a region involves the displacement and rotation at the measurement points. As the rotation cannot be measured directly, it is usually estimated from the difference of the displacement values of two adjacent points [27], which may introduce substantial inaccuracy. Hence the strain energy density is not included into the experimental analyses in this paper.

3. Results and discussions

There are 2 vibration modes, i.e., bending around $x$-axis and $y$-axis. Fig, 4 shows the first mode shapes for the bending vibration around $x$-axis and $y$-axis.

The first resonant frequencies for bending vibration around $x$-axis and $y$-axis are obtained by using the 3 methods (perturbation, FE, and experimental) for 2 damage levels (Table 3). The first resonant frequencies decrease as the damage increases. The variation rates of the first resonant frequencies of the damaged beam from those of the intact beam are listed in the parentheses.

When the beam is cut off 10 % in the height in the damage zone, the bending stiffness is reduced 27 % around $x$-axis and 10 % around $y$-axis; which results in the frequency being reduced 2 % and 1 % respectively for the bending vibration around $x$-axis and $y$-axis. When the beam is cut off 30 % in the height, the bending stiffness is reduced 66 % around $x$-axis and 30 % around $y$-axis; which results in the frequency being reduced 5 % and 2 % respectively for the bending vibration around $x$-axis and $y$-axis.

Fig. 4The first mode shapes for bending vibration around x-axis and y-axis

a) First bending vibration mode around $x$-axis

b) a) First bending vibration mode around $y$-axis

The 3 methods obtain the similar variations of the first resonant frequencies. For the bending vibration around $x$-axis, the perturbation method predicts 2 % and 5 % reduction in the first resonant frequency for the damage of 10 %, and 30% respectively; while the FE method predicts 2 % and 12 % reduction and the experimental method predicts 2 % and 11 % reduction. For the bending vibration around y-axis, the perturbation method predicts 1 % and 2 % reduction in the first resonant frequencies for the damage of 10 % and 30 % respectively; while the FE method predicts 1 % and 3 % reduction and the experimental method predicts 1 % and 2 % reduction.

The low sensitivity of the first eigenvalue to a damage of a structure has been demonstrated analytically by the beam model. As shown in Eq. (24), the denominator, ${\int }_0^{L}EI\cdot w_1^{\text{'}\text{'}}\cdot w_1^{\text{'}\text{'}}dz$, actually represents the strain energy of the whole beam for the first mode; while the integrals, ${\int }_{s-d∕2}^{s+d∕2}EI\cdot w_1^{\text{'}\text{'}}\cdot w_1^{\text{'}\text{'}}dz$ and ${\lambda }_1{\int }_{s-d∕2}^{s+d∕2}\rho A\cdot w_1\cdot w_{1}dz$, represent the strain energy and kinetic energy stored at the damaged part for the first mode. The sensitivity of the first eigenvalue to the damage is proportional to the energy changes at the damage region and inversely related to the strain energy of the whole beam. As the damage region, [$s-d/2,s+d/2$], is much smaller than the beam length$\mathrm{}L$, so the strain energy stored in the damage zone usually is much smaller than the strain energy of the whole beam. Furthermore, the coefficients $\varphi$ is the same order as $\alpha$, both are smaller than 1. This explains why the change rate of the first eigen-value, $∆{\lambda }_1/{\lambda }_1$, is less sensitive to the changes of the stiffness and mass in the damage zone. Additionally, the resonant frequency is less sensitive than the corresponding eigenvalue to the damage. The change of the first eigenvalue of the damaged beam is ${\tilde{\lambda }}_1/{\lambda }_1=\left(1+∆{\lambda }_1/{\lambda }_1\right)$, and the change of the corresponding (angular) frequency is ${\tilde{\omega }}_1/{\omega }_1$=$\sqrt{1+∆{\lambda }_1/{\lambda }_1}$. By using Taylor expansion for a small variation $∆{\lambda }_1/{\lambda }_1$, its first-order approximation is ${\tilde{\omega }}_1/{\omega }_1\approx$ (1+0.5$∆{\lambda }_1/{\lambda }_1$). This result demonstrates that the change rate of the first frequency is only half that of the corresponding eigenvalue.

The perturbation method is able to predict relatively accurate results for the frequency variations for a small damage (see Table 3). For the case with 10 % damage, the three methods (perturbation, FE, and experimental) obtained similar results (–2.1 %, –2.0 %, and –1.7 %) for bending around $x$-axis and (–0.8 %, –0.5 %, and –0.8 %) for bending around y-axis. For the case with 30 % damage, the three methods still produced similar results (–2.3 %, –2.7 %, and –2.4 %) for bending around $y$-axis; however, for bending around $x$-axis there existed evident difference (–5.0 %, –11.6 %, and –11.4 %) among the results predicted by the three methods. Firstly, for the sake of simplicity, the perturbation analyses only consider the first order increment, so it is more suitable for a small perturbation. As shown in Table 1, the reduction in the flexible rigidity (for bending around $x$-axis) is more than 87 % for the case with 30 % damage, thus the perturbation analyses would have evident inaccuracy for this case (Table 3). Secondly, the damage is not symmetric about the neutral axis of the original beam and the stress distribution would be very complex in the damage location, hence the assumption of plane section in the Euler beam theory is not really suitable in the damage region, which will also induce the inaccuracy in this part.

Fig. 5Variation rate of strain energy density within elements: a) perturbation, b) FE

a) First mode of the bending vibration around $x$-axis

b) First mode of the bending vibration around $y$-axis

The change rates of the strain energy density stored in each region (compared with the value of the intact beam) are shown in Fig. 5 for the 2 methods (perturbation and FE) and the 2 damage levels. Both the perturbation and FE methods obtain the similar pattern of the change rates. The strain energy density varies dramatically in the damage zone and less obviously somewhere else. For the bending vibration around $x$-axis, the perturbation method predicts 27 % and 65 % increase in the damaged part for the damage levels of 10 % and 30 % respectively; while the FE method predicts 18 % and 75 % increase in the damaged part. For the bending vibration around $y$-axis, the perturbation method predicts 10 % and 28 % increase in the damaged part for the damage levels of 10 % and 30 % respectively; while the FE method predicts 20 % and 145 % increase in the damaged part.

The variation of the strain energy density indeed reflects the damage degree of the beam (Fig. 5). For the regions outside the damage zone, the change rate of the strain energy density is very small, as the coefficients ${\zeta }_{ij}$ are very small; for the regions within the damage zone, the change rate of the strain energy density is evident, as the coefficients $\varphi$ is the same order as $\alpha$. For the bending vibration around $x$-axis, the change rate of the modal strain energy density, which is mostly influenced by$\mathrm{}\varphi$ (i.e., ${\varphi }_x$), can be as high as 27 % and 65 % respectively for the 2 damage levels; for the bending vibration around $y$-axis, the change rate of the modal strain energy is 10 % and 28 % respectively for the 2 damage levels. These conclusions are also confirmed by the FE results (Table 3), both show same trend.

4. Conclusions

A perturbation method has been proposed for prediction of the changes in the resonant frequencies and modal strain energy due to a small damage of a structure. The sensitivities of the frequency and modal strain energy are expressed explicitly in analytical formulae. These formulations provides us with a clearer view of the mathematical framework to explain why the lowest resonant frequency, as a global mechanical parameter, is less sensitive to a local defect in a structure and why the modal strain energy density is more sensitive to a damage.

The proposed perturbation method is able to assess the quality of a damage indictor. It can be implemented more easily and faster than FE and experimental methods. The perturbation analyses have revealed that the change rate of the first resonant frequency is related to the energy change in the damage zone and inversely related to the strain energy of the whole beam; thus, it is less sensitive to damage in a small region. While the change rate of the strain energy density is closely related to damage in a small region; hence, it is more sensitive to a damage than the first resonant frequency. The conclusions obtained from the perturbation analyses are validated with the FE simulations of the damaged beam.