Structural Engineering and Mechanics, Vol. 41, No. 1 (2012) 25-41

25

Crack identification based on Kriging surrogate modelHai-yang Gao*, Xing-lin Guoa and Xiao-fei Hub
State Key Laboratory of Structural Analysis for Industrial Equipment,
Dalian University of Technology, Dalian, P.R. China(Received March 20, 2011, Revised July 25, 2011, Accepted November 7, 2011)
Abstract. Kriging surrogate model provides explicit functions to represent the relationships between the
inputs and outputs of a linear or nonlinear system, which is a desirable advantage for response estimationand parameter identification in structural design and model updating problem. However, little research has
been carried out in applying Kriging model to crack identification. In this work, a scheme for crack
identification based on a Kriging surrogate model is proposed. A modified rectangular grid (MRG) is
introduced to move some sample points lying on the boundary into the internal design region, which will
provide more useful information for the construction of Kriging model. The initial Kriging model is thenconstructed by samples of varying crack parameters (locations and sizes) and their corresponding modal
frequencies. For identifying crack parameters, a robust stochastic particle swarm optimization (SPSO)
algorithm is used to find the global optimal solution beyond the constructed Kriging model. To improve
the accuracy of surrogate model, the finite element (FE) analysis soft ANSYS is employed to deal withthe re-meshing problem during surrogate model updating. Specially, a simple method for crack number
identification is proposed by finding the maximum probability factor. Finally, numerical simulations andexperimental research are performed to assess the effectiveness and noise immunity of this proposed
scheme.
Keywords: Kriging surrogate model; crack identification; stochastic particle swarm optimization; probability factor

1. Introduction
Fatigue damage commonly occurs in engineering structures under working conditions of highrotational velocity and impact loads, which may induce intensive stresses and dangerous damages.
A significant number of methods for health monitoring and damage detection based on finite
element model updating technique have been proposed and developed during recent years. A
comprehensive review was given by Doebling et al. (1998), of most related works, as well as future
research requirements.
Previous studies on damage identification can be categorized into two branches: the element-based
method and the DOF-based method. For the element-based method, the model updating strategy iscommonly employed, which is considered as some changes in physical properties of a structure*Corresponding author, Ph.D. Candidate, E-mail: eric415@126.comaProfessor, E-mail: xlguo@dlut.edu.cnb
Ph.D. Candidate, E-mail: beyondhxf@163.com

26

Hai-yang Gao, Xing-lin Guo and Xiao-fei Husuch as Young’s modulus, the section inertia, density and Poisson ratio of elements. Once the FEmodel has been constructed, the position of elements and nodes would never be changed during theidentification. For the DOF-based method, damage is often described as added massless export node
which connects two neighboring elements (e.g., crack). After that, the identification of crack
locations and sizes is carried out. Differently from element-based method, the construction of FE
model depends on damage locations and sizes which remain to be confirmed. Therefore, re-meshingtechniques should be considered.
Due to its particularity and uncertainty, the crack identification problem can hardly employ theelement-based method. In the early stages of the development of crack detection, the frequency
contour plot method (Nandwana and Maiti 1997, Chaudhari and Maiti 2000, Lele and Maiti 2002)
was one of the most favorite techniques to identify a single crack using the first few modal
frequencies since these parameters can be determined by measuring at only one point of thestructure. Gudmundson (1982, 1983) discussed the effect of geometrical imperfection on the
eigenvalue by means of perturbation analysis. Liang et al. (1991) proposed a massless rotational
spring model for the crack based on the Euler-Bernoulli beam theory, in which the location and sizeof a crack can be identified via finding the intersection point of a few frequency contour lines.
Ostachowicz and Krawczuk (1991) studied the forced vibrations of the beam and the effects of thecrack locations and sizes on the vibration behavior of the structure. Nandwana and Maiti (1997)
developed the frequency contour plot method from constant section beam model to stepped beamand truncated wedged beam. Lele and Maiti (2002) extended the frequency contour plot method inbeams based on Timoshenko beam theory. However, the frequency contour plot method suffersfrom the drawback that the curves of frequency contour plot might not intersect because of
inaccuracies in the modeling with respect to the measured results in many cases.
Many efforts have been devoted to single crack identification in beams. However, there are few
works on multiple cracks identification. In response, Hu and Liang (1993) first used the continuous
damage model to identify the discrete elements of a structure. Patil and Maiti (2003) adopted themathematical approach of Hu and Liang (1993) and applied the transfer matrix of vibration
modeling to multi-crack identification. Later, Shifrin and Ruotolo (1999) presented a strategy usingn + 2 characteristic equations for a beam with n cracks. The component mode synthesis method(Kisa and Brandon 2000) and modied Fourier series method (Zheng and Fan 2001a, b) were also
proposed for multi-crack identification. Unfortunately, a serious limitation is that the number ofcracks present in a beam is usually not known as a priori in a practical damaged engineering
structure.
Crack identification can be attributed to an inverse problem. Besides the direct least square
method, a large volume of iterative algorithms was reported. A simple sensitivity-based method wasapplied by Lee (2009a, b), which relies on the minimization of certain objective function based on
the residual errors between the measured and predicted frequencies; the Newton-Raphson iterative
procedure was carried out to identify n cracks in a beam where 2n natural frequencies of thecracked beam were strictly required. In his work, the mesh of beam model was regenerated at each
iterative step. Nevertheless, it is pointed out by Fang and Perera (2011) that proper selection of theinitial guesses is important for the convergence of the solution. Additionally, ill-conditioning andnon-uniqueness in the solution of inverse problem may appear as inevitable difficulties for a largescale structure.
To avoid the shortcoming of sensitivity-based method, artificial intelligence techniques such asneural networks and generic algorithms have been increasingly utilized owing to their excellent

Crack identification based on Kriging surrogate model

27

pattern recognition capability and the global convergence property (Atalla and Inman 1998, Zitzler
and Thiele 1999, Lee et al. 2005, Perera and Ruiz 2008). However, the required number of training
samples would exponentially increase, resulting in considerable computational efforts.
As a combination of mathematical and statistical techniques, surrogate model might constitute a
good alternative for parameter identification, which relaxes the strict requirement on thequantitative relation between inputs and outputs compared to sensitivity-based method. And it ispotential to build any mathematical relationship between input parameters and output responses inboth linear and nonlinear systems. Response surface model (RSM) for damage detection has
drawn much attention. In the thesis by Cundy (2002), a four-step process based on RSM wasdeveloped for simple physical systems and damage identification was performed successfully,
given the limited amount of “training” data used. Faravelli and Casciati (2004) and Casciati(2010) employed acceleration time histories collected under different loading conditions for
identifying the presence of damage and locating cracks by a comparison of the sum of thesquared errors (SSE) histograms. Moreover, a derivative of SSE, defined as ADM, was proposed
and used for detecting distributed cracks in an actual masonry (Casciati 2010). Huang et al.
(2011) pointed out that since RSM requires an understanding of the qualitative tendency of theentire design space, it is likely awkward when used for describing non-linearity commonly
appeared in complex problem.
Compared to RSM, Kriging surrogate model (KSM) has more flexibility to model response datawith multiple local extrema, which is already popular in the chemical and industrial design (Shyyet al. 2001, Gao and Wang 2008, Gao et al. 2008, Forrester and Keane 2009). However, to date,little research has been carried out in applying KSM to crack identification. In this work, we
explore an effective scheme for multi-crack identification based on KSM. The initial KSM isconstructed by the samples of crack parameters and their corresponding frequencies. To estimate
the actual crack parameters, then the SPSO algorithm is employed, which is demonstrated to beeffective, accurate and robust for searching for a global optimal solution. After that, the optimal
solution will be used in FE analysis and inserted into the initial sample set to update the initial
KSM, until the surrogate model is sufficiently accurate and the optimization process converges.
For a practical damaged engineering structure, a simple procedure for crack number identification
is proposed by finding the maximum probability factor. Numerical examples and experimental
test have been performed on cantilever beams. Finally, some appealing merits of the proposed
scheme were also discussed.

2. Finite element model
The Finite element model of a cantilever beam with double cracks is shown in Fig. 1(a).
Parameters αi = ai/h and βi = si/L (i = 1.2…) represent the normalized size and location of the ith
crack. In previous literatures, cracks are modeled as massless rotational springs or finite element
equation of a beam segment based on the Euler-Bernoulli theory or Timoshenko theory. Also some
empirical equations given by Ostachowicz and Krawczuk (1991) and Dimarogonas and Paipetis(1983) are employed to deal with the crack stiffness matrix for rectangular beam. However, it issomewhat discommodious when the FE models need to be re-meshed. In this study, a twodimensional plane elastic rectangular beam is modeled by using ANSYS code to avoid this trouble
instead of Euler-Bernoulli or Timoshenko beam theory, as shown in Fig. 1(b).

28Hai-yang Gao, Xing-lin Guo and Xiao-fei Hu

Fig. 1 Schematic: (a) cantilever beam with double cracks (b) two-dimensional mesh at crack tip

3 Basic theory of Kriging model
3.1 Sampling method
The basis of building a Kriging model is sample information. Here, a modified rectangular grid
(MRG) approach presented by Gao et al. (2008) is adopted to ensure its uniformity within design
region. Moreover, it can move some points lying on the boundary into the internal design region,
which will provide more useful information for the Kriging model, and it can ensure that the points
have less replicated coordinate values. A brief review is given in Fig. 2.
The advantage of this method is shown in Fig. 3 given by Gao et al. (2008). It can be seen that
MRG can avoid the case that the sample points are spaced close to each other, which may occur
using Latin Hypercube Sampling (LHS).
3.2 Construction of Kriging surrogate model
Kriging surrogate model is a statistics-based interpolated method (Sacks et al. 1989, Sakata et al.
2007). It maps the input parameters to the corresponding responses by a model function, which canbe written as

Crack identification based on Kriging surrogate model

29

Fig. 2 Flowchart of modified rectangular grid approach

Fig. 3 Three different methods of sample distributions: (a) RG (b) MRG (c) LHS

i

T

i

i

ŷ ( x ) = f ( x )β + z ( x )
i

i

i

i

i = 1, 2, …, n(1)

where x = [ x1, x2, …, xm ] is the ith sample points expressed as a m-dimension variable vector,iiŷ ( x ) is an approximate function corresponding to a linear or nonlinear function f ( x ) . β denotesithe matrix of regression coefficients to be estimated from the samples, and z ( x ) represents a model2
of Gaussian and stationary stochastic process with mean of zero and variance of σ . The covariancematrix of z ( x ) is

30

Hai-yang Gao, Xing-lin Guo and Xiao-fei Hu
i

j

2

i

j

Cov [ z ( x ), z ( x ) ] = σ R ( θ, x , x )
i

(2)

j

The n × n correlation function matrix between x and x can be formed by Gaussian correlation
function with only a single correlation parameter θ in this paper. It provides a smooth and infinitely
differentiable surface, given by
i

m

j

i

j 2

R ( θ, x , x ) = exp –θ ∑ ( xl – xl )

(3)

l=1

2

The unknown correlation parameters β and σ can be estimated by maximizing the log-likelihood
function (Jones et al. 1998)
2
2
1- ( Y – f T β )T R–1 ( Y – fT β ) + C--- n ln σ + ln R + ---ln ( β, σ , θ ) = – 1
2
2
σ
1

2

(4)

n

where Y = [ ŷ ( x ), ŷ ( x ), …, ŷ ( x ) ] , denotes the response value of samples, and C denotes a2
constant. By differentiating the log-likelihood function with respect to b and σ , respectively, andletting them be equal to zero. We can then obtain
T –1
2
T
T
--- ( Y – F β ) R ( Y – F β )σ̂ = 1
n

(5)

and
T

–1

–1

T

–1

β̂ = ( F R F ) ( F R Y )
1

n

(6)

T

where F = [ f ( x ), …, f ( x ) ] . This model leads to a best linear unbiased predictor and an*
associated mean squared error. The function value ŷ ( x ) at new points x* can be approximately
predicted by
*

* T

* T

*

2

–1

ŷ ( x ) = f ( x ) β̂ + r ( x ) R ( Y – Fβ̂ )(7)in which
*

1

*

n

*

r ( x ) = [ R ( θ , x , x ) , R ( θ, x , x ) , … , R ( θ , x , x ) ]T(8)Kriging model can not only predict responses by the new input parameters, but also estimate
optimal parameters satisfied with equality or inequality response constraints in both linear andnonlinear system (Simpson et al. 2001). When the original surrogate model is constructed, thequality of the model can be assessed according to the accuracy of predictions. The squared multiple
correlation coefficient (R2) and the empirical integrated squared error (EISE) criterion (Ren andChen 2010) can be used to select additional sample points for updating the initial KSM.

Crack identification based on Kriging surrogate modelN2∑ [ ŷj – yj ]
j=1

31

2

R = 1 – -------------------------N2∑ [ yj – ỹ ]

(9)

j=1

2
1 N
EISE = ---- ∑ [ ŷ j – y j ]Nj = 1(10)

where ŷj and yj are the jth component of the response vector of surrogate model and the true value
calculated via FE analysis respectively, ỹ is the mean of all true values, N is the length of vector ŷ .
4. Optimization procedure
The identification of crack locations and sizes is actually an inverse problem. The optimal
parameters can be determined by minimizing the discrepancies between the predicted responses
based on KSM and measured ones from the practical structure. To construct KSM, the samples ofvarying crack locations and sizes are used as the input parameters and their corresponding model
frequencies are adopted as the output responses. The single-objective function can be constructed byusing first few frequencies, which can be expressed asN*minΠ = ∑ W j ( λ̂j ( x ) – λ j )j=1

0 < Wj < 1*

lb < x < ub(11)where λ̂ is a vector of frequencies predicted by KSM at sample point x*, and λ represents themeasured frequencies caused by actual crack locations and sizes in a practical structure; Wj is
weighting factor to impose to the different order of natural frequencies. In this study, all its element
are set to 1. lb and up denote the lower and upper bounds.
To avoid the local optimal solution in conventional sensitivity-based method due to improper
selection of the initial guesses, a SPSO algorithm (Qi et al. 2007) is adopted for crack identification
based on the constructed Kriging model, which has been demonstrated to be effective, accurate androbust for searching a global optimal solution. When weight parameter is set to 0.0, the simple
updating procedure of the PSO can be expressed asXi ( t + 1 ) = Xi ( t ) + c1 r1 [ Pi ( t ) – Xi ( t ) ] + c2 r2 [ Pg ( t ) – Xi ( t ) ](12)where c1 and c2 are two positive constants called acceleration coefficients, Pi and Pg are local andglobal best locations, respectively, r1 and r2 are random numbers in the interval (0,1). If Pg = Pi,
Xi (t + 1) = Xi (t). Hence, the particle at the global best position will stop evolution. To improve theglobal searching ability, this algorithm randomly generates an extra particle labeled j with position
Xj to continue evolution in the search domain. It means that at least one particle is generated in the
searching domain randomly to improve the global searching quality of PSO algorithm.

32

Hai-yang Gao, Xing-lin Guo and Xiao-fei Hu

Fig. 4 Flowchart of crack identification based on Kriging surrogate model

5. Computation procedures
The proposed scheme can be organized as a surrogate model updating and optimal estimation
procedure as illustrated in Fig. 4. The identification of crack locations and sizes consists of thefollowing steps:
Step 1 Generate initial sample matrix X (crack locations and sizes) using MRG method and run
the simulation program ANSYS to obtain the corresponding output response Y (say modal
frequencies).
Step 2 Construct the initial Kriging surrogate model with sample parameters and output responses
obtained in Step 1,
Step 3 Find the optimal parameters x*k by minimizing discrepancies between response measured
in practical structures and that calculated by means of SPSO algorithm based on the initial
Kriging surrogate model, and set the iterative index k = 1,
Step 4 Check criterions: If the current response y*k predicted by ANSYS and the one based on
Kriging surrogate model satisfy the given criterions (both R2 > 0.98 and EISE < 0.01), then stop
Kriging surrogate model updating and predict the parameters with this surrogate model; else, go
to Step 5,
Step 5 Add x*k and y*k behind the initial sample points generated in Step 1, then update the
Kriging model with setting the iterative index k = k + 1,
Step 6 Loop to Step 3 and repeat the process till the criterions are satisfied.
6. Numerical examples
6.1 Single crack in a beamFirst of all, four cases of single crack in a beam are performed to assess the proposed method.
The followings are the properties adopted: Young’s modulus E = 210 GPa, Poisson’s ratio µ = 0.3,
density ρ = 7850 Kg/m3, length L = 5 m, height h = 0.02 m. Finite element analysis model is built byusing ANSYS code. The crack parameters (locations and sizes) are identified using SPSO algorithm

Crack identification based on Kriging surrogate model

33

Table 1 Comparison of actual and identified crack parameters after 100 generationsαβThe number of Kriging surrogate
model updating step (k)ActualIdentified

ActualIdentified

A

1

0.1

0.09920.2

0.1996B

1

0.2

0.19980.50.5004C

1

0.3

0.30060.6

0.5998D

1

0.4

0.40050.80.8000Case

Fig. 5 Kriging model for each frequency as function of depth (size) and location after 100 generations inCase B

based on Kriging surrogate model. To establish the KSM, 30 groups of initial samples in theinterval of (0~1) and the corresponding first four calculated model frequencies are selected as the
input and output, respectively. The c1 and c2 of SPSO are set to 1.0. The Swarm Size and Maximum
Swarm are set to 20 and 100, respectively. The inverse problem of identifying single crack in a
cantilever beam is solved for four cases (A, B, C and D) as shown in Table 1.
From Table 1, it can be found that the identified crack parameters are very close to the assumed
actual ones for all cases. It means that the proposed scheme has a good accuracy for single crack
detection in a beam. Fortunately, the initial Kriging surrogate model (k = 1) is precise enough topredict the crack parameters. And more importantly, there is no limitation for initial guesses of each

34

Hai-yang Gao, Xing-lin Guo and Xiao-fei Hu

Fig. 6 The sum modal frequency residuals of each mode during SPSO for Case A to Case D
Table 2 Parameters setting of initial sample and SPSO algorithm
Case E

Parameters

Case F

Swarm Size

30

50

Maximum Swarm

300

500

c1 and c2

1.0

1.0

Frequency used

First 5 modes

First 6 modes

The number of sample points

200

550

Table 3 Comparison of actual and identified crack parameters in Case E and F
Case

k

α1

α2

α3

α4

α5

β1

β2

β3

β4

β5

0.10

0.15

0.18

-

-

0.35

0.55

0.75

-

-

0.093

0.143

0.186

-

-

0.347

0.552

0.742

-

-

0.08

0.10

0.15

0.18

0.20

0.20

0.30

0.40

0.50

0.60

0.075

0.106

0.142

0.169

0.211

0.198

0.310

0.392

0.496

0.611

E
ActualIdentified

2

F
ActualIdentified

3

parameter which is a common problem in sensitivity-based approach. The relations between crack
parameters and first four modal frequencies in Case B can be extracted from Kriging model asshown in Fig. 5. Fig. 6 illustrates the convergence of cost function residue during SPSO. The
horizontal axis is the number of generation while the vertical axis is the sum modal frequency
residuals of each mode.

Crack identification based on Kriging surrogate model

35

Fig. 7 The sum modal frequency residuals during SPSO algorithm and Kriging model updating for Case E and F

6.2 Multiple cracks in a beam
The presented scheme can easily be extended to a beam with multiple cracks. In this section, twocases (E with three cracks; F with five cracks) are considered to validate the effectiveness of the
proposed scheme. The identification procedure is carried out with parameters setting in Table 2, andthe actual and identified crack parameters are given in Table 3. The cost function values in Case E
and F of each generation are plotted in Fig. 7.
It was found that the proposed scheme gave satisfactory predictions for multi-crack cases. The
errors of the estimated crack locations and sizes are within 6.25 percent of the actual ones. At thesame time it is interesting to remark that the identified crack locations seem to be more accurate
than the identified crack sizes. One reason could be that the location parameters are more sensitive
than size parameters to the model frequencies of the structure. Particularly, this method can beextended to detect any number of cracks without the strict requirement on the quantitative relation
between inputs and outputs compared to sensitivity-based method. In this sense, Kriging surrogate
model provides us a better choice for estimating with a moderate number of variables (less than 50)(Simpson et al. 2001).
To further study the effect of referring mode number on the identification result, Case F wasrecalculated by using 4, 5 and 7 referring modal frequencies, respectively. The results are listed inTable 4, in which the required number of Kriging model updating steps is also given. It was foundthat more updating steps should be made to describe a relatively accurate surrogate model when asmall number of referring modes is used. All results are uniquely given by SPSO algorithm.
Remarkably, the identifications with 6 and 7 modes yield very similar results due to a more prefectsurrogate model compared to the ones with 4 or 5 modes.

36

Hai-yang Gao, Xing-lin Guo and Xiao-fei Hu

Table 4 Identification results using different number of referring modes
Number of referring modes

k

α1

4

7

0.071 0.114 0.152 0.162 0.245 0.192 0.322 0.444 0.475 0.632

5

4

0.075 0.104 0.162 0.178 0.188 0.202 0.292 0.395 0.502 0.610

63

0.075 0.106 0.142 0.169 0.211 0.198 0.310 0.392 0.496 0.611

7

3

0.076 0.106 0.145 0.170 0.210 0.198 0.311 0.392 0.496 0.611

α2

α3

α4

α5

β1

β2

β3

β4

β5

Fig. 8 (a) Dimensional drawing; (b) a steel beam with two cracks

7. Experimental study
7.1 Problem description
The validity of this scheme has been verified through simulated examples. For experimental
measurements, it is expected that there would be some deviation due to noise originating fromenvironment as well as electronic devices. In this section, a steel beam with two cracks cut by wire
EDM is used to assess the feasibility and anti-noise ability of this proposed scheme, as shown in
Fig. 8, which has the elastic modulus of E = 209.2 GPa, Poisson's ratio of µ = 0.25, mass density ofρ = 7830 Kg/m3, length of L = 1.013m, width of S = 0.088 m, thickness of t = 0.009 m. A twodimensional finite element numerical analysis model of the beam is built by using ANSYS program.
As we know, frequency discrepancy between modal predictions and experimental results is usually
generated due to inevitable difficulties or oversights in the design, material property and boundary
condition. Therefore, model updating strategy should be employed for minimizing this discrepancy
before crack identification. An undamaged beam with same material property is prepared for
updating the properties of the FE model.
During the experiment, a 6-channel B&K 3050 data acquisition system with Pulse 13.1 software,
an impact hammer with an accelerometer is employed as shown in Fig. 9. The acquired frequency
response functions are then inputted into the modal analysis software ME’scope 4.0 to identify the
first six frequencies of undamaged and cracked beams in free boundary condition, respectively, asCrack identification based on Kriging surrogate model

37

Fig. 9 Experimental setup for crack identification

Table 5 First six theoretical and measured frequencies
Mode

1

2

3

4

5

6

Measured frequency (Hz)439.661158.62152.32538.43339.14656.3

Theoretical frequency after model updating(Hz)440.591159.72156.32537.13341.64660.7

Measured frequency (Hz)425.611111.32138.52528.73385.44536.9Theoretical frequency after model updating (Hz)423.951115.72135.42535.83389.34543.8

Undamaged beam

Cracked beamshown in Table 5. It can be seen that relative accurate FE model can be constructed after model
updating, which lays the foundation for crack identification based on the Kriging model.
7.2 An additional procedure for crack number identification
Like other methodologies, the number of cracks present in a beam is usually not known as apriori in a practical damaged engineering structure, which is a serious limitation in its application.
Some strategies (Hu and Liang 1993, Lee 2009, Lam and Ng 2008) were presented to overcome
this difficulty.
Here a simple method related to SPSO algorithm is presented for crack number identification.
Different from Newton-Raphson iteration method, SPSO algorithm can provide a converged
solution in most case. It means that this searching algorithm ensures the existence of solutions. The
value of cost functions should be a constant in this work. When the assumed number of cracks isequal to the actual one, the value of cost function should be minimal. That is, we can find the mostprobable number of cracks by looking for the number which yields the minimum value of costfunction. Fig.10 shows the procedure for crack number identification. Firstly, a single crack is38

Hai-yang Gao, Xing-lin Guo and Xiao-fei Hu

Fig. 10 Flowchart of crack number identification

Table 6 Comparison of actual and identified crack parameters in experimental case
Case
G

α1

k
2

β1

α2

β2

Actualidentified

Actualidentified

Actualidentified

Actualidentified

0.28410.29020.29620.30050.17050.16830.54290.5473assumed and crack identification is carried out, then we calculate the cost function (CF). The values
under conditions of several different numbers of cracks in the beam are calculated in the samemanner. The probability factor is defined as–1

CF ( M )P ( M ) = ----------------------------------Pm

(13)

–1

∑ ( CF( M ) )M=1

where M denotes the crack number we assumed to calculate the CF, Pm is the maximum number ofcracks to be considered. More assumptive calculations should be carried out when the largest valueof P(M) appears. Finally, the number of cracks can be then determined by comparing with the P(M)values. The maximum value of P(M) indicates the most probable number of cracks. Numerical casewill be discussed in the next section.
7.3 Results and discussion
To construct the Kriging model and identify the location and size of cracks, 20 groups of sample
parameters corresponding to their first four theoretical model frequencies are selected as the input
and output, respectively. The bounds of crack sizes are (0~0.5), and that of the crack locations are(0~1). The c1 and c2 of SPSO are set to 1.0. The Swarm Size and Maximum Swarm are set to 20
and 400, respectively. The identified results are shown in Table 6. It can be found that the twocracks can be located with the errors of 1.45% and 0.81% and the observed errors of the crack sizes
are 2.15% and 1.29%. The validity of Kriging surrogate model seems to be insensitive to random

Crack identification based on Kriging surrogate model

39

Fig. 11 Probability factor P(M) under conditions of five different assumed numbers of cracks

perturbations and noise information (Gao 2008). To check the crack number, P(M) (M = 1.2…5) is
calculated under conditions of five different numbers of cracks in the beam, as shown in Fig. 11.
The maximum value of P(M) indicates a correct assumption of the crack number.
Although the proposed scheme shows satisfactory predictions for multi-crack identification, it
should be mentioned that the natural frequency is the characteristics of the whole dynamic system.
For complex structure, frequencies are not adequate since cracks result more in local features of thestructure. Therefore, more responses such as model shapes, vibration amplitudes and frequency
response functions might be taken into account for multiobjective SPSO. In addition, a simple
method for adding sample point as illustrated in Section 5 was used to updating the Kriging modelin this work. To improve the efficiency of Kriging model updating, the use of some advanced point
addition criterions should be encouraged. Sakata et al. (2003, 2004) introduced an empiricalsemivariogram based semivariogram fitting approach for large-scale sampling problems on the
optimization of stiffened straight cylinder and wing structure. Such an approach may serve as an
excellent module for present method to improve the efficiency and accuracy of the surrogate modelreconstruction.

8. Conclusions
The main aim of the present work is to provide an efficient scheme for crack identification basedon Kriging surrogate model. The identification results in numerical and experimental cases show
good performance of this scheme. Some appealing merits can be summarized as follow:• Kriging surrogate model is applied to provide a simple relationship between crack parametersand corresponding modal frequencies for avoiding the expensive FE analysis at every iterative
step of optimization.
• A modified rectangular grid approach was adopted to move inward the sample points on theboundary and ensure the uniformity within design region, which makes the initial Kriging model

40

Hai-yang Gao, Xing-lin Guo and Xiao-fei Hu

more accurate inside the design space.
• A robust SPSO algorithm was used to solving a global optimization problem. The initial guessesof the crack parameters could be avoided, which is commonly difficult in Newton-Raphson
iteration method.
• For a practical damaged engineering structure, the number of cracks is usually not known. In
order to solve this issue, a simple and effective procedure is presented related to SPSO
algorithm, in which the most probable number of cracks can be obtained by finding the
maximum of P(M) factor.
To assess the performance of the proposed scheme, Numerical studies of a cantilever beam withsingle and multiple cracks were carried out. The results show that the identified crack locations andsizes are in good agreement with the actual ones. Moreover, experimental research is also
considered by using a cracked steel beam structure under free-free boundary condition. The
proposed scheme is then proved to be effective and promising for crack identification in engineering.

References
Atalla, M.J. and Inman, D.J. (1998), “On model updating using neural networks”, Mech. Syst. Signal Pr., 12(1),135-161.
Casciati, S. (2010), “Response surface models to detect and localize distributed cracks in a complex continuum”,
J. Eng. Mech. ASCE, 136(9), 1131-1142.
Casciati, S. (2010), “Statistical approach to a SHM benchmark problem”, Smart Struct. Syst., 6(1), 17-27.
Chaudhari, T.D. and Maiti, S.K. (2000), “A study of vibration of geometrically segmented beams with andwithout crack”, Int. J. Solids Struct., 37(5), 761-779.
Cundy, A.L. (2002), “Use of response surface metamodels in damage identification of dynamic structures”,
Master Thesis, Virgina Polytechnic Institute and State University, Virginia, US.
Dimarogonas, A.D. and Paipetis, S.A. (1983), Analytical Methods in rotor dynamics, Elsevier Applied Science,
London.
Doebling, S.W., Farrar, C.R. and Prime, M.B. (1998), “A summary review of vibration-based damage
identification methods”, Shock. Vib., 30(2), 91-105.
Fang, S.E. and Perera, R. (2011), “Damage identification by response surface based model updating using Doptimal design”, Mech. Syst. Signal Pr., 25(2), 717-733.
Faravelli, L. and Casciati, S. (2004), “Structural damage detection and localization by response change
diagnosis”, Prog. Struct. Eng. Mater., 6(2), 104-115.
Forrester, A.I.J. and Keane, A.J. (2009), “Recent advances in surrogate-based optimization”, Prog. Aerosp. Sci.,
45(1-3), 50-79.
Gao, Y.H. (2008), “Optimization methods based on Kriging surrogate model and their application in injection
molding”, Ph. D. Thesis, Dalian University of Technology, Dalian, PRC.
GaoY.H., Turng, L.S. and Wang, X.C. (2008), “Adaptive geometry and process optimization for injection
molding using the Kriging surrogate model trained by numerical simulation”, Adv. Poly. Tech., 27(1), 1-16.
Gao, Y.H. and Wang, X.C. (2008), “An effective warpage optimization method in injection molding based on the
Kriging model”, Int. J. Adv. Manuf. Tech., 37(9-10), 953-960.
Gudmundson, P. (1982), “Eigenfrequency changes of structures due to cracks, notches or other geometricalchanges”, J. Mech. Phys. Solids., 30(5), 339-353.
Gudmundson, P. (1983), “The dynamic behaviour of slender structures with cross-sectional cracks”, J. Mech.
Phys. Solids., 31(4), 329-345.
Hu, J. and Liang, R.Y. (1993), “An integrated approach to detection of cracks using vibration characteristics”, J.
Frankin Institute, 330(5), 841-853.
Huang, Z., Wang, C., Chen, J. and Tian, H. (2011), “Optimal design of aeroengine turbine disc based on Kriging

Crack identification based on Kriging surrogate model

41

surrogate models”, J. Comput. Struct., 89(1-2), 27-37
Jones, D.R., Schonlau, M. and Welch, W.J. (1998), “Efficient global optimization of expensive black-box
functions”, J. Glob. Optim., 13(4), 455-492.
Kisa, M. and Brandon, J.A (2000), “Free vibration analysis of multiple open-edge cracked beams by component
mode synthesis”, Struct. Eng. Mech., 10(l), 81-92.
Lam, H.F. and Ng, C.T. (2008), “A probabilistic method for the detection of obstructed cracks of beam-type
structures using spatial wavelet transform”, Probab. Eng. Mech., 23(2-3), 237-245.
Lee, J. (2009a), “Identification of multiple cracks in a beam using vibration amplitudes”, J. Sound Vib., 326(1-2),
205-212.
Lee, J. (2009b), “Identification of multiple cracks using natural frequencies”, J. Sound Vib., 320(3), 482-490.
Lee, J.J., Lee, J.W., Yi, J.H., Yun, C.B. and Jung, H.Y. (2005), “Neural networks-based damage detection forbridges considering errors in baseline finite element models”, J. Sound Vib., 280(3-5), 555-578.
Lele, S.P. and Maiti, S.K. (2002), “Modelling of transverse vibration of short beams for crack detection andmeasurement of crack extension”, J. Sound Vib., 257(3), 559-583.
Liang, R.Y., Choy, F.K. and Hu, J. (1991), “Detection of cracks in beam structures using measurements ofnatural frequencies”, J. Frankin Institute, 328(4), 505-518.
Nandwana, B.P. and Maiti, S.K. (1997), “Detection of the location and size of a crack in stepped cantilever
beams based on measurements of natural frequencies”, J. Sound Vib., 203(3), 435-446.
Ostachowicz, W.M. and Krawczuk, M. (1991), “Analysis of the effect of cracks on the natural frequencies of a
cantilever beam”, J. Sound Vib., 150(2), 191-201.
Patil, D.P. and Maiti, S.K. (2003), “Detection of multiple cracks using frequency measurements”, Eng. Fract.
Mech., 70(12), 1553-1572.
Perera, R. and Ruiz, A. (2008), “A multistage FE updating procedure for damage identification in large-scale
structures based on multiobjective evolutionary optimization”, Mech. Syst. Signal Pr., 22(4), 970-991.
Qi, H., Ruan, L.M., Zhang, H.C., Wang, Y.M. and Tan, H.P. (2007), “Inverse radiation analysis of a onedimensional participating slab by stochastic particle swarm optimizer algorithm”, Int. J. Therm. Sci., 46(7),649-661.
Ren, W.X. and Chen, H.B. (2010), “Finite element model updating in structural dynamics by using the response
surface method”, Eng. Struct., 32(8), 2455-2465.
Sacks, J., Welch, W.J., Mitchell, T.J. and Wynn, H.P. (1989), “Design and analysis of computer experiments”,
Stat. Sci., 4(4), 409-435.
Sakata, S., Ashida, F. and Zako, M. (2003), “Structural optimization using Kriging approximation”, Comput.
Methods Appl. Mech. Eng., 192(7-8), 923-939.
Sakata, S., Ashida, F. and Zako, M. (2004), “An efficient algorithm for Kriging approximation and optimization
with large-scale sampling data”, Comput. Methods Appl. Mech. Eng., 193(3-5), 385-404.
Sakata, S., Ashida, F. and Zako, M. (2007), “On applying Kriging-based approximate optimization to inaccuratedata”, Comput. Meth. Appl. Mech. Eng., 196(13-16), 2055-2069.
Shifrin, E.I. and Ruotolo, R. (1999), “Natural frequencies of a beam with an arbitrary number of cracks”, J.
Sound Vib., 222(3), 409-423.
Shyy, W., Papila, N., Vaidyanathan, R. and Tucker, K. (2001), “Global design optimization for aerodynamics androcket propulsion components”, Prog. Aerosp. Sci., 37(1), 59-118.
Simpson, T.W., Poplinski, J.D., Koch, P.N. and Allen, J.K. (2001), “Meta-models for computer-based engineering
design: survey and recommendations”, Eng. Comput., 17(2), 129-150.
Zheng, D.Y. and Fan, S.C. (2001a), “Natural frequencies of a non-uniform beam with multiple cracks via
modied Fourier series”, J. Sound Vib., 242(4), 701-717.
Zheng, D.Y. and Fan, S.C. (2001b), “Natural frequency changes of a cracked Timoshenko beam by modied
Fourier series”, J. Sound Vib., 246(2), 297-317.
Zitzler, E. and Thiele, L. (1999), “Multiobjective evolutionary algorithms: a comparative case study and thestrength Pareto approach”, IEEE T. Evolut. Comput., 3(4), 257-271.