An improved two-step method based on Kriging model for beam structures

Abstract

An improved two-step method for multiple cracks detection in beam structures is presented in this paper. First, the crack locations are identified by the mode shape curvatures. Second, samples of varying crack depths at the crack locations and corresponding natural frequencies are used to construct the initial Kriging model in place of the crack detection database. Finally, a particle swarm optimization algorithm is employed to detect the crack depths and update the Kriging model by an optimal point adding process. The numerical simulation and experimental verification of a cantilever beam with two cracks demonstrate that the proposed method is effective. The detection results indicate that the Kriging model can effectively reduce computation time with a similar accuracy. The detection results also indicate that the more natural frequencies used, the higher accuracy of crack depths detected.

Introduction

Beam structures appear in a wide range of structural systems such as civil, mechanical, and aeronautical engineering. The crack in structures can lead to significant loss of strength, and, even worse, resulting in catastrophic failure of structures. Therefore, how to monitor and evaluate crack should be considered seriously in both theoretical and experimental research.

A significant number of researches for health monitoring and crack detection have been proposed during recent years.^1–3 Some researchers focused on the characteristics of structures with initial crack, including fatigue crack propagation, residual strength, and residual life. Several kinds of nonprobabilistic interval conceptions are introduced by Qiu and Wang. ⁴ Wang et al. ⁵ established a new model of nonprobabilistic time-dependent reliability, which integrates interval process theory with the first-passage approach, to achieve the life-cycle safety estimation for cracked structures. A novel method of the time-dependent reliability, which combines the perturbation series expansions with the interval mathematics, was presented for life prediction of fatigue crack growth problems. ⁶ For crack detection problems, because modal shapes and natural frequencies can be obtained accurately and easily in practice, many researchers employed modal shapes and natural frequencies to identify crack locations and depths.^7–10

Generally, the methods using modal shape and natural frequencies for crack detection involve forward and inverse problem analysis. ¹¹ The purpose of forward problem analysis is to develop a numerical simulation model and construct a crack detection database. In the inverse problem analysis, the modal shape calculated via simulation model is used to detect the crack locations, and the measured natural frequencies are used as the inputs to an optimization algorithm to determine the crack depths from the crack detection database. However, the computational cost may be unacceptable because of the repeated analysis of computationally expensive finite element models during the process of constructing the crack detection database. So it is desirable to construct a simple relationship between crack depths and natural frequencies to replace the crack detection database. The surrogate model may be a good alternative to construct the relationship between crack depths and natural frequencies.

Recently, surrogate models were investigated by many researchers. Cho ¹² performed an investigation using response surface method to predict the accumulated cracks in concrete structures. And it was shown that the response surface method could be employed to predict the probability of failure for such cracks. Fang and Perera ¹³ proposed a response surface methodology based on model updating using D-optimal design for damage identification. Gao et al. ¹⁴ proposed an effective method based on a Kriging model for crack detection in a plate. The result shows that this method can effectively be used for identifying the crack parameters. Yang et al. ¹⁵ presented a new surrogate model based on model updating method using frequency response function for damage identification. The proposed method has good accuracy and robustness, which has been verified by a numerical simulation of a cantilever and experimental test data of a laboratory three-story structure. Ghadami et al. ¹⁶ presented an adaptive multiple cracks detection algorithm in beam-like structures. The main advantage of this method is that it can detect adaptively the unknown number of cracks.

Compared with other surrogate models, Kriging model has a higher accuracy, and it can provide both predictive location and prediction error, which is very suitable for global optimization algorithm. ¹⁷ Gao et al. ¹⁸ presented an idea for the application of Kriging model in crack detection and has a great interest of attentions but still the study of multiple cracks detection in beams with Kriging model is limited.

In this paper, we present an improved two-step method for detecting multiple cracks in beam structures. In the first step, the mode shape curvature is employed to detect the crack locations. The crack depths at the identified locations are then detected in the second step by using particle swarm optimization (PSO) algorithm. In the second step, the meshless local Petrov–Galerkin (MLPG) model is first constructed to simulate the multicracked beam. Then, the initial Kriging model, expressing the relationships between natural frequencies and crack depths, is constructed by samples of various crack depths and corresponding natural frequencies. The initial samples are generated by Latin hypercube sampling (LHS), and corresponding natural frequencies are calculated via MLPG model. The PSO algorithm is used to determine the crack depths by minimizing the objective function related to the measured natural frequencies. After that, the current optimal solution will be used in MLPG analysis and inserted into the initial sample set to update the initial Kriging model, until the surrogate model meets the accuracy requirements. Finally, numerical simulation and experimental investigation demonstrate the effectiveness and efficiency of the presented method.

Theory and formulation

In this section, the theory and formulation of Kriging model is presented in the “Construction of Kriging model” subsection. The initial samples of Kriging model are produced by LHS. “Construction of cracked beam” subsection describes how to construct the model of cracked beam. The model is used to calculate modal shape and corresponding natural frequencies of the cracked beam.

Construction of Kriging model

The Kriging model is an unbiased estimation model with a minimum estimation variance. Its basic concept was proposed by Krige and further deduced by Matheron. ¹⁹ The regression equation of Kriging model can be expressed as

$$\widehat{y}(x)=f{(x)}^{\text{T}}{\hspace{0.17em}}_{p\times 1}{\beta }_{p\times 1}+z(x)$$

(1)where $f{(x)}^{\text{T}}{\hspace{0.17em}}_{p\times 1}$ is the basis function for regression model, p is the number of basis function in regression model, ${\beta }_{p\times 1}$ is the regression coefficients, and $z(x)$ is a q-dimensional stochastic process with statistical properties as follows

$$\begin{array}{l}E(z(x))=0\\ E(z(\omega )z(x))={\sigma }^2\widehat{R}(\theta ,\omega ,x)\end{array}$$

(2)where ${\sigma }^2$ is the variance of the stochastic process and ${\sigma }^2\widehat{R}(\theta ,\omega ,x)$ is the correlation functions with parameters θ. The correlation functions include Gauss, spline, and so on. In this paper, the Gaussian correlation function is used to construct the Kriging model. Suppose there are m sampling points, the predictor value of an unknown point can be defined by

$$\widehat{y}(x)=\omega {(x)}^{\text{T}}Y$$

(3)where $\omega (x)={[{\omega }_1,{\omega }_2,...,{\omega }_m]}^{\text{T}}$ is the weight function. $Y={[y_1,y_2,...,y_m]}^{\text{T}}$ . Because the Kriging model is an unbiased estimation model and has a minimum estimation variance, ω can be expressed as

$$\begin{array}{l}\omega \left(x\right)=R^{\text{-}1}\left(r\text{-}F{(F^{\text{T}}{R}^{\text{-1}}F)}^{\text{-1}}\right(F^{\text{T}}{R}^{\text{-1}}r\text{-}f\left)\right)\\ R_{m\times m}=\left[\begin{array}{ccc}R\left(x_1,x_1\right)& \cdots & R\left(x_1,x_m\right)\\ ⋮& \ddots & ⋮\\ R\left(x_m,x_1\right)& \cdots & R\left(x_m,x_m\right)\end{array}\right]\\ F_{m\times p}=\left[\begin{array}{ccc}{f}_1\left(x_1\right)& \cdots & f_p\left(x_1\right)\\ ⋮& \ddots & ⋮\\ f_1\left(x_m\right)& \cdots & f_p\left(x_m\right)\end{array}\right]\end{array}$$

(4)

For the regression problem $F\beta \cong Y$ , we have

$$\widehat{\beta }={(F^{\text{T}}{R}^{\text{-1}}F)}^{\text{-1}}{F}^{\text{T}}{R}^{\text{-1}}Y$$

(5)

By substituting equations (4) and (5) into equation (3), the predictor can be expressed as

$$\widehat{y}(x)=f{(x)}^{\text{T}}\widehat{\beta }+r{(x)}^{\text{T}}{R}^{-1}(Y-F\widehat{\beta })$$

(6)where ${{\displaystyle r}}_{m\times 1}^{\text{T}}(x)=\left[R(\theta ,x,x^1),\dots ,R(\theta ,x,x_m)\right]$ .

When the initial Kriging model is constructed, its accuracy can be assessed by using the squared multiple correlations (SMCs) and the empirical integrated squared error (EISE) criteria^20,21 as follows

$$\begin{array}{l}S\mathit{MC}=1-\frac{{\displaystyle {\sum }_{j=1}^N{[{\widehat{y}}_j-y_j]}^2}}{{\displaystyle {\sum }_{j=1}^N{[y_j-{\overline{y}}_j]}^2}}\\ E\mathit{ISE}=\frac{1}{N}{\displaystyle \sum _{j=1}^N{[{\widehat{y}}_j\mathrm{-}{y}_j]}^2}\end{array}$$

(7)where ${\widehat{y}}_j$ and $y_j$ denote the jth component of response vector of Kriging model and the true value calculated via MLPG analysis, respectively. ${\overline{y}}_j$ is the mean of all true values from the real cracked beam.

Construction of cracked beam

The model of multicracked cantilever beam is shown in Figure 1. l_x is the length of cantilever beam, and h is the height of beam. The normalized parameters α = t_i/h, β = x_i/l_x are used to describe the crack locations and crack depths, respectively. Hence, the values of crack locations and depths are limited from 0 to 1. The free vibration equation of a beam with uniform cross-section can be defined by ²²

$$\mathit{EI}\frac{{\partial }^{4}u(x,t)}{\partial x^4}+\rho A\frac{{\partial }^{2}u(x,t)}{\partial t^2}=00<x<l_x$$

(8)where u(x,t) is the amplitude of the transverse displacement of the beam axis, ρ is the material density, A is the cross-sectional area, and EI is the bending stiffness. By assuming that the beam axis is under time-harmonic vibration, u(x,t) can be expressed as

$$u(x,t)=u(x)e^{\mathit{iwt}}=u{e}^{\mathit{iwt}}$$

(9)

Figure 1.

Multicracks cantilever beam model.

Substituting equation (9) back into equation (8), equation (8) can be simplified as

$$\mathit{EI}{u}^{(4)}={\omega }^2\rho \mathit{Au}0<x<l_x$$

(10)where ω is the natural frequency. The boundary conditions of the beam are given as

$$\left[\begin{array}{cc}u(x)=u_0\in {\Gamma }_u& \frac{\partial u(x)}{\partial x}={\theta }_0\in {\Gamma }_{\theta }\\ M=M_0\in {\Gamma }_M& V=V_0\in {\Gamma }_V\\ {\Gamma }_u\cap {\Gamma }_V=\varnothing & {\Gamma }_{\theta }\cap {\Gamma }_M=\varnothing \end{array}\right]$$

(11)where M is the moment and V is the shear force. ${\Gamma }_u,{\Gamma }_{\theta },{\Gamma }_M,{\Gamma }_V$ are the displacement, slope, moment, and shear force on the boundary, respectively.

For the cracked beam, the continuity conditions at the crack location $x=x_c$ are

$$\begin{array}{l}u(x_c^{-})=u(x_c^{+})\\ u^{(2)}(x_c^{-})=u^{(2)}(x_c^{+})\\ u^{(3)}(x_c^{-})=u^{(3)}(x_c^{+})\end{array}$$

(12)

And the discontinuity condition is ²³

$$u^{(1)}(x_c^{+})-u^{(1)}(x_c^{-})=\lambda u^{(2)}(x_c^{-})$$

(13)where λ is the flexibility coefficient of the cracks having a depth of α and it can be expressed as ²⁴

$$\lambda =5.346hg(\xi )$$

(14)where

$$\begin{array}{l}g(\xi )=1.8624{\xi }^2-3.95{\xi }^3+16.375{\xi }^4-37.226{\xi }^5+76.81{\xi }^6\\ -126.9{\xi }^7+172{\xi }^8-143.97{\xi }^9+66.56{\xi }^{10}\\ \xi \text{=}\alpha \text{/}h\end{array}$$

(15)

To apply local weighted-residual equation, ²⁵ equation (10) can be rewritten as

$${\displaystyle {\int }_{{\Omega }_q^{(i)}}{\varpi }_i(\mathit{EI}{u}^{(4)}-{\omega }^2\rho \mathit{Au})\text{d}\Omega }=0$$

(16)where ${\varpi }_i={\varpi }_i(x)$ is weight function. ${\Omega }_q^{(i)}$ is local sub-domain located entirely inside the global domain $\Omega$ .

Choosing a reasonable weight function, equation (16) can be rewritten as

$$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_q^{(i)}}{\varpi }_i{\hspace{0.17em}}^{(2)}EI{u}^{(2)}\text{d}\Omega }-{\displaystyle {\int }_{{\Omega }_q^{(i)}}{\varpi }_i{\omega }^2\rho \mathit{Au}\mathrm{d\Omega }}+[{\varpi }_{i}EI{u}^{(3)}{]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_u}-[{\varpi }_i{\hspace{0.17em}}^{(1)}EI{u}^{(2)}{]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_{\theta }}\\ +[{\varpi }_{i}EI{u}^{(3)}{]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_V}-[{\varpi }_i{\hspace{0.17em}}^{(1)}EI{u}^{(2)}{]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_M}=0\end{array}$$

(17)where ${\Gamma }_q^{(i)}$ is the boundary of ${\Omega }_q^{(i)}$ .

The moment M and shear force V related to the displacement u are

$$M=M_0=\mathit{EI}{u}^{(2)}V=V_0=-\mathit{EI}{u}^{(3)}$$

(18)

Substituting equation (18) into equation (17), we have

$$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_q^{(i)}}{\varpi }_i{\hspace{0.17em}}^{(2)}EI{u}^{(2)}\text{d}\Omega }-{\displaystyle {\int }_{{\Omega }_q^{(i)}}{\varpi }_i{\omega }^2\rho \mathit{Aud}\Omega }+{\left[{\varpi }_{i}EI{u}^{(3)}\right]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_u}-{\left[{\varpi }_i{\hspace{0.17em}}^{(1)}EI{u}^{(2)}\right]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_{\theta }}\\ -{\left[{\varpi }_i{V}_0\right]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_V}-{\left[{\varpi }_i{\hspace{0.17em}}^{(1)}{M}_0\right]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_M}=0\end{array}$$

(19)

When a crack occurs in ${\Omega }_q^{(i)}$ , equation (19) can be rewritten as

$$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_q^{(i)}}{\varpi }_i{\hspace{0.17em}}^{(2)}EI{u}^{(2)}\text{d}\Omega }-{\displaystyle {\int }_{{\Omega }_q^{(i)}}{\varpi }_i{\omega }^2\rho \mathit{Aud}\Omega }+{\left[{\varpi }_{i}EI{u}^{(3)}\right]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_u}-{\left[{\varpi }_i{\hspace{0.17em}}^{(1)}EI{u}^{(2)}\right]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_{\theta }}\\ +\lambda \mathit{EI}{\varpi }_i{\hspace{0.17em}}^{(2)}(x_c^{-})u^{(2)}(x_c^{-})-{\left[{\varpi }_i{V}_0\right]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_V}-{\left[{\varpi }_i{\hspace{0.17em}}^{(1)}{M}_0\right]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_M}=0\end{array}$$

(20)

The MLS approximation u^h(x) of dependent variable u (x) over the global domain $\Omega$ can be defined by

$$u^h(x)={\displaystyle \sum _{i=1}^n{\varphi }_i(x)u_i}$$

(21)where ${\varphi }_i(x)$ is the MLS shape function. Substituting equation (21) and weight function ${\varpi }_i={\varpi }_i(x)$ into equation (20), we can obtain the frequency equation

$$(K+K_c-{\omega }^{2}M)U=F$$

(22)where K is the stiffness matrix of the beam without a crack, K_c is the correction matrix for K due to the crack, M is the mass matrix, and F is the external load vector. The entries in K, K_c, and M are defined by

$$K_i={\displaystyle {\int }_{{\Omega }_q^{(i)}}{\varpi }_i{\hspace{0.17em}}^{(2)}EI{\varphi }_i{\hspace{0.17em}}^{(2)}\text{d}\Omega }+{\left[{\varpi }_{i}EI{\varphi }_i{\hspace{0.17em}}^{(3)}\right]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_u}-{\left[{\varpi }_i{\hspace{0.17em}}^{(1)}EI{\varphi }_i{\hspace{0.17em}}^{(2)}\right]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_{\theta }}$$

(23)

$$K_{c,i}=\lambda \mathit{EI}{\varpi }^{(2)}(x_c^{-}){\varphi }_i{\hspace{0.17em}}^{(2)}(x_c^{-})$$

(24)

$$M_i={\displaystyle {\int }_{{\Omega }_q^{(i)}}{\varpi }_i\rho A{\varphi }_i\text{d}\Omega }$$

(25)

$$F_i={\left[{\varpi }_{i}V\right]}_0{\hspace{0.17em}}_{{\Gamma }_q^{(i)}\cap {\Gamma }_V}+{\left[{\varpi }_i{\hspace{0.17em}}^{(1)}{M}_0\right]}_{{\Gamma }_q^{(i)}\cap {\Gamma }_M}$$

(26)

For an arbitrary number of cracks, equation (24) can be rewritten as

$$K_{c,i}={\displaystyle \sum _{n_c}{\lambda }_{n_c}}EI{\varpi }^{(2)}(x_{c,n_c}^{-}){\varphi }_i{\hspace{0.17em}}^{(2)}(x_{c,n_c}^{-})$$

(27)where n_c is the number of cracks. We can see clearly that the matrix K_C represents the cracks’ softening effect on the stiffness matrix K of the beam. By adding it into the matrix K, a general stiffness matrix for the cracked beam can be obtained as

$$\overline{K}=K+K_c$$

(28)

Substituting equation (28) into equation (22), we have

$$(\overline{K}-{\omega }^{2}M)U=F$$

(29)where ω (rad/s) is the natural circular frequency, $f=\omega /2\pi$ (Hz) is the natural frequency and U is corresponding mode of the cracked beam.

Crack detection method

As mentioned above, the crack locations and depths are detected separately in the improved two-step method. The flow chart of the proposed crack detection method is shown in Figure 2. First, the mode shape curvature is employed to detect the crack locations. Second, the initial Kriging model is constructed by samples of various crack depths and corresponding natural frequencies. Then an optimal point adding process for the Kriging model updating is carried out to improve the accuracy of Kriging model. Finally, the measured natural frequencies are employed as inputs to PSO algorithm to determine the crack depths. The details associated with the two steps are given in the “Crack locations detection” and “Crack depths evaluation” subsections.

Figure 2.

The flow chart of improved two-step method. MLPG: meshless local Petrov–Galerkin; PSO: particle swarm optimization.

Crack locations detection

Obtain an arbitrary modal shape of the cracked beam

In the simulation, the mth modal shape of cracked beam is calculated via the MLPG method. It should be pointed that the finite element model can also be used to obtain the modal shape.

2. Calculate the mode shape curvature

If the mth modal shape Y_m is measured, each point V_m,j of mth mode shape curvature V_m can be calculated from the mth modal shape using a central difference approximation as follows ²⁶

$$V_{m,j}=\frac{Y_{m,j+1}-2Y_{m,j}+Y_{m,j-1}}{h_c{\hspace{0.17em}}^2}$$

(30)where h_c is the length of two neighbor nodes; j is the node at different spatial locations; and $Y_{m,j-1}$ , $Y_{m,j}$ , and $Y_{m,j+1}$ are three continuous nodes of Y_m.

• 3. Detect n crack locations from mode shape curvature

The crack locations can be detected by the peaks or sudden changes in the mode shape curvature of cracked beam.

Crack depths evaluation

Construct MLPG model to simulate cracked beam

The MLPG model can be constructed following the procedure proposed in the “Construction of cracked beam” subsection. The correction stiffness matrix in the associated crack locations is used to simulate the cracked beam, whereas the global mass matrix remains unchanged.

2. Construct the initial Kriging model

As mentioned above, the Kriging model, expressing the relationships between natural frequencies and crack depths, is constructed by samples of various crack depths and corresponding natural frequencies. The initial samples are generated by LHS, and corresponding natural frequencies are calculated via MLPG model.

3. Update the Kriging model by an optimal point adding process

The optimal points are used to update the current Kriging model until the Kriging model meets the accuracy requirements. The SMC and EISE criteria are used to decide whether the model should be updated.

4. Evaluate crack depths using PSO algorithm based on the Kriging model

The crack depths detection is essentially an optimization problem. The PSO algorithm is employed to detect the crack depths of the n cracks, and the first q measured natural frequencies are used as inputs. We use the following objective function to search for the optimal crack depth based on Kriging model ²⁷

$$\begin{array}{l}\min {{\displaystyle \sum _{j=1}^q\Vert \frac{{\overline{f}}_j-f_j^{*}}{{\overline{f}}_j}\Vert }}_2\\ s.t.\hspace{0.5em}0.1<{\alpha }_i<0.9i=1,2,\dots ,n\end{array}$$

(31)where ${\Vert \cdot \Vert }_{\text{2}}$ is the Euclidean norm. ${\overline{f}}_j=F_j({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_n)$ is the discrete function of the n crack depths ${\alpha }_i(i=1,2,\dots ,n)$ , which is calculated by Kriging model. The constraint $0.1<{\alpha }_i<0.9$ limits the search space from 0.1 to 0.9. $f_j^{*}$ is the measured natural frequency, where $j=1,2,\dots ,q$ . It should be pointed out that q ≥ n is the necessary condition to obtain a robust solution of unknown ${\alpha }_i(i=1,2,\dots ,n)$ .

To make a clear description on how to evaluate crack depths, we summarized as follows:

Step 1: Generate the initial sampling matrix X (crack depths) by LHS, and calculate corresponding natural frequencies Y by MLPG analysis.

Step 2: The initial Kriging model is constructed with samples X and output Y obtained in step 1.

Step 3: Apply the PSO algorithm to search the optimal crack depth by minimizing discrepancies between the measured natural frequencies $f^{*}$ and the calculated natural frequencies $\overline{f}$ based on the initial Kriging model.

Step 4: Check criteria ( $\mathit{SMC}>0.995\text{ and }\mathit{EISE}<0.005$ ). If the Kriging model meets the criteria, then stop updating model and outputting n crack depths; else, go to step 5.

Step 5: Add x* and y^* into the X and Y generated in Step 1, respectively. Then update the current Kriging model.

Step 6: Go to step 3 and repeat the process till the criteria are satisfied, and output n crack depths.

Numerical simulation

In order to verify the efficiency of the presented method, a numerical study of a cantilever beam with two cracks is carried out. The parameters of the beam are length l_x = 0.5 m, width l_y = 0.012 m, height h = 0.02 m, Young’s modulus E = 210 GPa, Poisson’s ratio μ = 0.3, material density ρ = 7860 kg/m³. The normalized parameters α = t_i/h, β = x_i/l_x are used to describe the crack locations and the crack depths. Crack case is considered as β₁ = 0.2, β₂ = 0.4, and α₁ = α₂ = 0.1.

The second modal shape shown in Figure 3(a) is available on a 101 × 1 sample grid, which is dense enough to detect crack locations. The second mode shape curvature calculated by equation (30) is shown in Figure 3(b). The two peaks indicate the detected crack locations, i.e. β₁^* = 0.198, β₂^* = 0.396, and the detected results are round to the actual locations of cracks. Therefore, the crack locations can be detected by the mode shape curvature.

Figure 3.

Crack location detection using the second modal shape and the second mode shape curvature. (a) The second modal shape and (b) the second mode shape curvature.

Eight cases with different depth combinations (α₁, α₂) for the beam with two cracks are shown in Table 1. The first five noise-free frequencies associated with the eight cases calculated by MLPG method are also shown in Table 1.

Table 1.

The first five noise-free frequencies for crack cases.

Case	Crack depths		Noise-free frequencies
	α 1	α 2	f 1	f 2	f 3	f4	f5
1	0.1	0.2	66.20	413.95	1159.17	2347.34	3936.67
2	0.1	0.3	65.47	407.87	1142.56	2336.51	3847.29
3	0.2	0.1	65.57	415.92	1164.46	2332.47	3946.35
4	0.2	0.4	62.28	388.35	1096.96	2282.45	3616.63
5	0.3	0.1	63.38	414.50	1159.60	2288.87	3879.50
6	0.3	0.4	60.49	385.82	1095.05	2227.42	3559.20
7	0.4	0.2	56.48	407.49	1139.81	2166.55	3700.84
8	0.4	0.3	56.11	400.43	1126.62	2144.05	3620.61

To simulate the random noise, a random number generator, rand in MATLAB, is used to simulate ±e% noise, i.e.

$$f_j^N=(1+e\times (2\times \mathit{rand}-1))\times {\overline{f}}_j$$

(32)where ${\overline{f}}_j$ is the noise-free natural frequencies, $f_j^N$ is the noise-contaminated frequencies, rand is the random number, and e is the noise level.

Suppose ±2% (e = 2) random noise is added to the noise-free frequencies according to equation (32), the first five noise-contaminated frequencies are shown in Table 2.

Table 2.

The first five noise-contaminated frequencies for crack cases.

Case	Noise-contaminated frequencies (±2% noise)
	f 1	f 2	f 3	f4	f5
1	65.23	411.10	1143.02	2380.56	3890.68
2	66.11	413.73	1131.13	2382.72	3868.29
3	64.57	410.30	1176.41	2372.94	4002.46
4	63.12	394.43	1093.72	2240.25	3645.73
5	64.30	417.11	1138.92	2277.50	3901.42
6	61.41	386.87	1081.78	2224.30	3541.58
7	57.51	404.16	1139.13	2153.53	3761.95
8	56.12	397.43	1109.10	2149.77	3683.64

In order to verify the efficiency of Kriging model for crack depths detection, we use 20 samples and corresponding first five noise-contaminated natural frequencies to construct the initial Kriging model. The optimization is implemented with MATLAB using a PSO Toolbox coded by Birge. ²⁸ According to the advice of Birge, the population size is set to 25, and the learning factors c₁ and c₂ are set to 2. The evaluation results by using Kriging model are shown in Table 3. For the sake of comparison, the crack depths detection by using crack detection database is also applied. To obtain a more accurate crack detection database, α₁ and α₂ are varied from 0.1 to 0.9 with step length of 0.01. Therefore, there are 6561 (=81 × 81) data points in the crack detection database. The evaluation results by using crack detection database are also shown in Table 3. As shown in the table, we can see clearly that the Kriging model can significantly reduce computation time with a similar accuracy.

Table 3.

The crack depths evaluation results using first five noise-contaminated frequencies.

Case	Evaluation results using crack detection database				Evaluation results using Kriging model
	α 1 *	ε 1 a	α 2 *	ε 2 b	α 1 *	ε 1 a	α 2 *	ε 2 b	Kriging model updating step
1	0.10	0	0.27	7	0.1000	0	0.2692	6.92	5
2	0.10	0	0.27	3	0.1000	0	0.2718	2.82	5
3	0.19	1	0.10	0	0.1813	1.87	0.1000	0	6
4	0.25	5	0.38	2	0.2358	3.58	0.3847	1.53	4
5	0.27	3	0.16	6	0.2745	2.55	0.1549	5.49	3
6	0.26	4	0.41	1	0.2727	2.73	0.4058	0.58	4
7	0.39	1	0.21	1	0.3931	0.69	0.1960	0.4	5
8	0.40	0	0.29	1	0.3991	0.09	0.2967	0.33	4
Mean error		1.75		2.63		1.44		2.26

^aThe identification error, ε₁=|α₁−α₁^*|×100%.

^bThe identification error, ε₂=|α₂−α₂^*|×100%.

To examine the effect of the number of natural frequencies on the detection accuracy, we use 20 samples and corresponding first three natural frequencies to construct the initial Kriging model. The evaluation results of crack depths are shown in Table 4. Compared with Table 3, the detection results using first five natural frequencies (the mean errors of crack depths are 1.44 and 2.26%, respectively) are better than the results using first three natural frequencies (the mean errors of crack depths are 4.10 and 3.96%, respectively). Therefore, it is preferred to use as many natural frequencies as possible to improve the detection accuracy of crack depths.

Table 4.

The crack depths evaluation results using first three noise-contaminated frequencies.

Case	Evaluation results using first three noise-contaminated frequencies				Kriging model updating step
	α 1 *	ε 1 a	α 2 *	ε 2 b
1	0.1657	6.57	0.2738	7.38	4
2	0.1088	0.88	0.2769	2.31	4
3	0.2556	5.56	0.1000	0	4
4	0.1311	6.89	0.3931	0.69	6
5	0.2553	4.47	0.2094	10.94	4
6	0.2283	7.17	0.4084	0.84	6
7	0.3899	1.01	0.2573	5.73	4
8	0.3975	0.25	0.3381	3.81	4
Mean error		4.10		3.96

^aThe identification error, ε₁=|α₁−α₁^*|×100%.

^bThe identification error, ε₂=|α₂−α₂^*|×100%.

The above results demonstrate that the Kriging model is effective in evaluating crack depths when the data contain certain level of noise.

Experimental verification

In this section, experiments of cantilever beam with two cracks are constructed to verify the validity of the proposed method. The experimental setup is shown in Figure 4. The parameters of the beam are length l_x = 0.5 m, width l_y = 0.012 m, height h = 0.02 m, Young’s modulus E = 206 GPa, Poisson’s ratio μ = 0.3, material density ρ = 7860 kg/m³.

Figure 4.

Experimental setup.

The experiments are conducted using two beams. One is intact and the other has two cracks. The locations and depths of two cracks are β₁ = 0.16, β₂ = 0.76, α₁ = 0.42, and α₂ = 0.42, respectively. As we know, in a typical impact test, the accelerometer is attached to a single point on the beam (in the experiments, the accelerometer is fixed near the left end of the beam), and the hammer is used to impact it as many points to define its mode shape. In the experimental study, the sampling frequency f_s is 5000 Hz and 10,000 data points are collected at each impact points. The first three measured frequencies of cracked beam are f₁^* = 54.5 Hz, f₂^* = 338 Hz, f₃^* = 869 Hz.

The second modal shape is shown in Figure 5(a) and the second mode shape curvature is shown in Figure 5(b). The peaks in Figure 5(b) indicate the detected crack locations, i.e. β₁^* = 0.16, β₂^* = 0.76, which are the exact locations of the two cracks.

Figure 5.

Two crack locations detection using the second mode shape curvature. (a) The second modal shape and (b) the second mode shape curvature.

Sometimes, the large difference between the measured frequencies and the calculated frequencies will make detection results ineffective. To avoid this issue, the “zero-setting” method is used to reduce the difference. The “zero-setting” method is defined by ²⁹

$$\frac{E}{E_j^{*}}={\left(\frac{{\overline{f}}_j}{f_j^{*}}\right)}^2$$

(33)where $E_j^{*}(j=1,2,...,q)$ is the corrected value of Young’s modulus E, which can be obtained from equation (33) for each natural frequency. This method can greatly reduce the error between measured frequencies and calculated frequencies, which are caused by boundary conditions and material parameters.

The measured frequencies of the intact beam are f₁ = 58.5 Hz, f₂ = 338 Hz, f₃ = 869 Hz and the calculated frequencies of the intact beam are f₁ = 68.23 Hz, f₂ = 394.59 Hz, f₃ = 1101.48 Hz. Therefore, Young’s modulus calculated by equation (33) is modified as E₁^* = 151.44 GPa, E₂^* = 157.48 GPa, E₃^* = 139.37 GPa.

There are 20 samples and corresponding first three natural frequencies calculated by E_m^* used to construct the initial Kriging model. And the measured frequencies of the cracked beam are adopted as the inputs of PSO algorithm. The detection results are shown in Figure 6, i.e. α^*₁ = 0.3973, α₂^* = 0.3785. The relative errors between the estimated and the actual depths are ε₁ = |α₁−α₁^*|×100%=2.27%, ε₂ = |α₂−α₂^*|×100%=4.15%. Therefore, the presented method can obtain the results with reasonable accuracy.

Figure 6.

The PSO convergence progress and the optimal particle location represented two crack depths.

Conclusion

An improved two-step method is presented for detecting multiple cracks in beam structures, which is based on mode shape curvature and Kriging surrogate model. First, the crack locations are detected from the mode shape curvature. Then, the initial Kriging model is constructed by samples of various crack depths and corresponding natural frequencies. And an optimal point adding process is carried out to update the initial Kriging model until the Kriging model is sufficiently accurate. Finally, the PSO algorithm is used to evaluate crack depths based on the Kriging model. Both of the simulation and experimental results demonstrate that the presented method is effective to detect multiple cracks in beam structures. The detection results demonstrate that Kriging model can significantly reduce computation time with a similar accuracy. The detection results also indicate that the more natural frequencies used, the higher accuracy of crack depths detected.