# Homotopy Iteration Algorithm for Crack Parameters Identification with Composite Element Method

Authors: Ling Huang; Zhongrong Lv; Weihuan Chen; Jike Liu
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

An approach based on homotopy iteration algorithm is proposed to identify the crack parameters in beam structures. In the forward problem, a fully open crack model with the composite element method is employed for the vibration analysis. The dynamic responses of the cracked beam in time domain are obtained from the Newmark direct integration method. In the inverse analysis, an identification approach based on homotopy iteration algorithm is studied to identify the location and the depth of a cracked beam. The identification equation is derived by minimizing the error between the calculated acceleration response and the simulated measured one. Newton iterative method with the homotopy equation is employed to track the correct path and improve the convergence of the crack parameters. Two numerical examples are conducted to illustrate the correctness and efficiency of the proposed method. And the effects of the influencing parameters, such as measurement time duration, measurement points, division of the homotopy parameter and measurement noise, are studied.


---

## Body

## 1. Introduction

Identification of crack parameters using the dynamic responses with finite element method (FEM) has been studied extensively for many years. In order to improve the degree of accuracy of vibration analysis, most of the approaches [1, 2] are based on densely refining finite element mesh. But these methods may lead to a large amount of computations. Unlike the classical FEM, Zeng [3] proposed a composite element method (CEM) combining the conventional FEM and classical theory (CT) to solve the structural dynamical problems. And Lu and Law [4] further improve the CEM using certain special boundary conditions of the beam. It is known that the CEM has two available approaches, theh-version and thec-version. As the former by increasing the number of element mesh is just like that ofp-version of the FEM [3], the latter by increasing the number of analytical functions will reduce the number of degrees of freedom in the FEM. And the latter can obtain the fine approximate solution with less computational effort in structural dynamics. Thus, a super convergence equation of motion of the structure can be established by using thec-version of CEM.There exist many reports on damage detection using structural dynamic response. Some of the methods are based on the time domain [5–8]. Cattarius and Inman [5] used the time histories of vibration response of the structure to identify damage in smart structures. Lu and Liu [6] made use of the dynamic responses of bridge under a moving vehicle to identify both the local damage of bridge and vehicular parameters. Lu and Law [7] proposed an approach for structural damage identification based on response sensitivity analysis. However, the response sensitivity approach is a local convergence approach, and it requires the initial values of the unknowns to be close to the true values. Then, Lu et al. [8] proposed a two-step approach based on mode shape curvature and response sensitivity analysis for crack identification.As both the crack location and depth are unknowns in the identification, a method with large convergence range should be sought. The homotopy method [9–12] which is based on the concept of homotopy is a widely convergent method for solving system of equations. It has been extended for other algorithm and applied in many fields; for example, He [10] proposed a homotopy perturbation method, which combines the homotopy in topology and traditional perturbation method, to provide an approximate solution to a wide range of linear and nonlinear problems. Liao [11] proposed a homotopy analysis method for solving highly nonlinear problems in science, finance, and engineering areas. Alexander and Yorke [12] employed a homotopy continuation method to solve the fixed points and singularities of vector fields and bifurcation problems.In this paper, an approach for identifying the parameters of cracks in a cracked beam based on homotopy iteration algorithm is presented. A fully open crack model with composite element model is adopted to establish the dynamics equation of the Euler-Bernoulli beam system. In inverse problem, the identification equation is derived by the minimization of the error between the calculated acceleration response and the simulated measured one. The equation is solved iteratively to identify the crack parameters. Meanwhile, Newton’s method with the homotopy equation is used to track the correct path and improve the convergence of the crack parameters. A simply supported beam and a cantilever beam are studied to illustrate the accuracy and efficiency of the proposed method. Both single and multiple cracks in the beam can be identified successfully using several measured acceleration responses. And the effects of measurement time duration, measurement points, division of the homotopy parameter, and measurement noise on the identified results are investigated.


## 2. Forward Problem

### 2.1. Crack Model

Figure1 shows a simple beam with multiple fully opened cracks along its length. It is assumed that the cracks have uniform depth across the width of the beam and that the mass distribution along the beam does not change. The stiffness, for the simplicity, is assumed to be varied linearly with a triangular reduction in a local region as Sinha et al. [13] proposed. That is, the crack only leads to local stiffness reduction in a specified length adjacent to the crack. The effective length of the stiffness reduction for the crack islc=1.5d. As shown in Figure 2, the flexural rigidityEI(x)adjacent to theith crack can be written as
(1)EI(x)={EI0-E(I0-Ici)(x-xi1)(xi-xi1),xi1≤x≤xci,EI0-E(I0-Ici)(xi2-x)(xi2-xi),xci≤x≤xi2,
whereEis Young’s modulus andI0=wd3/12andIci=w(d-hci)3/12are the moments of inertia of the intact and cracked cross sections, respectively.wanddare the width and height of the beam, respectively.xi1=xci-lcandxi2=xci+lcare the positions where the reduction of the flexural rigidity begins and finishes, respectively. If the beam has multiple cracks, the same procedure can be followed to calculate the flexural rigidity of other cracks.Figure 1
A simply supported beam with cracks.Figure 2
Variation on flexural rigidityEI(x)adjacent to a crack.


### 2.2. Dynamic Response of the Structure with CEM

As shown in Figure1, the beam is discretized into one element together with several terms eigenfunction of classical theory. Lu and Law [4, 7] have proved that only using one finite element is effective in thec-version technique of CEM. And it can reduce the total number of degrees of freedom in the finite element model. For the crack identification, what is important is that we do not need to judge whether the cracks affect one or more elements stiffness as the crack(s) is always in one finite element for the entire beam. Otherwise, it needs to calculate the stiffness distribution due to the crack(s) along the beam in each iteration this one-beam-one-element strategy simplifies greatly the process in crack identification.In the CEM, the displacement fielddCEMis combined with the FEM part and the CT part. The FEM part of the displacement field should satisfy the nodal boundary conditions, and it can be expressed as the product of the shape function matrixN(x)and the nodal displacement vectorq(t),
(2)dFEM(x,t)=N(x)q(t),
where
(3)N(x)=[1-3x2L2+2x3L3,x-2x2L+x3L2,h3x2L2-2x3L3,x3L2-x2L],
whereLis the length of the beam.The second part can be obtained by the linear combination of the multiplication of analytical shape functionϕ(x)with a vector of coefficientc(t)(4)dCT(x,t)=∑r=1nϕr(x)cr(t),(r=1,2,…,n),
whereϕr(x)is a special value which can be obtained according to the boundary conditions of the beam andnis the number of the mode functions used from the CT. Such that, for a simply supported beam,ϕr(x)=sin(rπx/L)can be selected as Lu and Law [4] proposed. And the coefficientcr(t)denotes the contribution of theϕr(x)in the total displacement field. For a uniform beam, less terms of the shapes function in CT also can be obtained with accuracy results [4]. And the number of termsncan be determined by a frequency convergence test as the previous study by Lu and Law [7]. The frequency convergence criterion is defined as
(5)maxr=1,2,…,10|Δωrnωrn|<Tol,
whereωrnis the estimation of therth circular frequency withn-terms in the CT.Δωrn=ωrn-ωrn-1is the difference of therth circular frequency obtained with then-terms and(n-1)-terms.Tolis the tolerance value.The displacement field of the CEM for a uniform Euler-Bernoulli beam element can be written as(6)dCEM(x,t)=dFEM(x,t)+dCT(x,t)=S(x)Q(t),
whereS(x)=[N(x),ϕ1(x),ϕ2(x),…,ϕn(x)]is the generalized shape function of CEM.The stiffness matrix of the cracked beam can be obtained from the following equation:(7)K=∫0Ld2S(x)Tdx2EI(x)d2S(x)dx2dx.For a beam with multiple cracks, as the whole crack beam is represented as a single finite element, the whole length of the beam can be divided according to the varied flexural rigidity. ForEI(x)adjacent to theith crack can be obtained from (1); the global stiffness matrixKcan be obtained by superposition.As the effect of crack on the mass of the beam is neglected and the consistent mass matrix of the cracked beam can be obtained from the following equation:(8)M=∫0LS(x)TρAS(x)dx,
whereρis the mass densityAis the area of the cross section.After introducing the boundary condition, the equation of motion of the forced vibration of the cracked Euler-Bernoulli beam is expressed as(9)MQ¨+CQ˙+KQ=G(t).In this study, Rayleigh damping model [14] is adopted; that is,C=a1M+a2K, wherea1anda2are two constants obtained from two different undamped natural frequencies,ωiandωj, and their associated modal damping ratios,ξiandξj, with the expression ofa1=2ωjωi(ωjξi-ωiξj)/(ωj2-ωi2)anda2=2(ωjξj-ωiξi)/(ωj2-ωi2). The vector of generalized accelerationQ¨, velocityQ˙, and displacementQof the structure can be obtained from (9) by direct integration. Further, the vector of physical accelerationd¨(x), velocityd˙(x), and displacementd(x)can be obtained from(10a)d¨(x)=S(x)Q¨,(10b)d˙(x)=S(x)Q˙,(10c)d(x)=S(x)Q.The generalized force vector isG(t)=S(xl)Tg(t), whereg(t)is an external force acting at the location ofxlfrom the left support.


## 2.1. Crack Model

Figure1 shows a simple beam with multiple fully opened cracks along its length. It is assumed that the cracks have uniform depth across the width of the beam and that the mass distribution along the beam does not change. The stiffness, for the simplicity, is assumed to be varied linearly with a triangular reduction in a local region as Sinha et al. [13] proposed. That is, the crack only leads to local stiffness reduction in a specified length adjacent to the crack. The effective length of the stiffness reduction for the crack islc=1.5d. As shown in Figure 2, the flexural rigidityEI(x)adjacent to theith crack can be written as
(1)EI(x)={EI0-E(I0-Ici)(x-xi1)(xi-xi1),xi1≤x≤xci,EI0-E(I0-Ici)(xi2-x)(xi2-xi),xci≤x≤xi2,
whereEis Young’s modulus andI0=wd3/12andIci=w(d-hci)3/12are the moments of inertia of the intact and cracked cross sections, respectively.wanddare the width and height of the beam, respectively.xi1=xci-lcandxi2=xci+lcare the positions where the reduction of the flexural rigidity begins and finishes, respectively. If the beam has multiple cracks, the same procedure can be followed to calculate the flexural rigidity of other cracks.Figure 1
A simply supported beam with cracks.Figure 2
Variation on flexural rigidityEI(x)adjacent to a crack.


## 2.2. Dynamic Response of the Structure with CEM

As shown in Figure1, the beam is discretized into one element together with several terms eigenfunction of classical theory. Lu and Law [4, 7] have proved that only using one finite element is effective in thec-version technique of CEM. And it can reduce the total number of degrees of freedom in the finite element model. For the crack identification, what is important is that we do not need to judge whether the cracks affect one or more elements stiffness as the crack(s) is always in one finite element for the entire beam. Otherwise, it needs to calculate the stiffness distribution due to the crack(s) along the beam in each iteration this one-beam-one-element strategy simplifies greatly the process in crack identification.In the CEM, the displacement fielddCEMis combined with the FEM part and the CT part. The FEM part of the displacement field should satisfy the nodal boundary conditions, and it can be expressed as the product of the shape function matrixN(x)and the nodal displacement vectorq(t),
(2)dFEM(x,t)=N(x)q(t),
where
(3)N(x)=[1-3x2L2+2x3L3,x-2x2L+x3L2,h3x2L2-2x3L3,x3L2-x2L],
whereLis the length of the beam.The second part can be obtained by the linear combination of the multiplication of analytical shape functionϕ(x)with a vector of coefficientc(t)(4)dCT(x,t)=∑r=1nϕr(x)cr(t),(r=1,2,…,n),
whereϕr(x)is a special value which can be obtained according to the boundary conditions of the beam andnis the number of the mode functions used from the CT. Such that, for a simply supported beam,ϕr(x)=sin(rπx/L)can be selected as Lu and Law [4] proposed. And the coefficientcr(t)denotes the contribution of theϕr(x)in the total displacement field. For a uniform beam, less terms of the shapes function in CT also can be obtained with accuracy results [4]. And the number of termsncan be determined by a frequency convergence test as the previous study by Lu and Law [7]. The frequency convergence criterion is defined as
(5)maxr=1,2,…,10|Δωrnωrn|<Tol,
whereωrnis the estimation of therth circular frequency withn-terms in the CT.Δωrn=ωrn-ωrn-1is the difference of therth circular frequency obtained with then-terms and(n-1)-terms.Tolis the tolerance value.The displacement field of the CEM for a uniform Euler-Bernoulli beam element can be written as(6)dCEM(x,t)=dFEM(x,t)+dCT(x,t)=S(x)Q(t),
whereS(x)=[N(x),ϕ1(x),ϕ2(x),…,ϕn(x)]is the generalized shape function of CEM.The stiffness matrix of the cracked beam can be obtained from the following equation:(7)K=∫0Ld2S(x)Tdx2EI(x)d2S(x)dx2dx.For a beam with multiple cracks, as the whole crack beam is represented as a single finite element, the whole length of the beam can be divided according to the varied flexural rigidity. ForEI(x)adjacent to theith crack can be obtained from (1); the global stiffness matrixKcan be obtained by superposition.As the effect of crack on the mass of the beam is neglected and the consistent mass matrix of the cracked beam can be obtained from the following equation:(8)M=∫0LS(x)TρAS(x)dx,
whereρis the mass densityAis the area of the cross section.After introducing the boundary condition, the equation of motion of the forced vibration of the cracked Euler-Bernoulli beam is expressed as(9)MQ¨+CQ˙+KQ=G(t).In this study, Rayleigh damping model [14] is adopted; that is,C=a1M+a2K, wherea1anda2are two constants obtained from two different undamped natural frequencies,ωiandωj, and their associated modal damping ratios,ξiandξj, with the expression ofa1=2ωjωi(ωjξi-ωiξj)/(ωj2-ωi2)anda2=2(ωjξj-ωiξi)/(ωj2-ωi2). The vector of generalized accelerationQ¨, velocityQ˙, and displacementQof the structure can be obtained from (9) by direct integration. Further, the vector of physical accelerationd¨(x), velocityd˙(x), and displacementd(x)can be obtained from(10a)d¨(x)=S(x)Q¨,(10b)d˙(x)=S(x)Q˙,(10c)d(x)=S(x)Q.The generalized force vector isG(t)=S(xl)Tg(t), whereg(t)is an external force acting at the location ofxlfrom the left support.


## 3. Inverse Problem

### 3.1. Objective Function

The problem of crack identification can be treated as an optimization problem. In the identification, the unknowns consist of the crack locationxcand crack depthhc; that is,α=[xc,hc]. If there aremmeasured points in the beam structure, the measured acceleration response can be expressed asd¨*=(d¨1,d¨2,…,d¨m). The objective function for the optimization problem can be derived by minimizing the error between the calculated acceleration response and the simulated measured one as
(11)J=12∑i=1nt∑j=1m(d¨(α,t)-d¨*(α,t))2,i=1,2,…,nt,j=1,2,…,m,
wheretdenotes the the measured time instants,nt=t/dtis the number of time step, anddtis the increment of time step.The inverse problem can be solved by minimizing the objective function. Taking partial derivative with respect to the identification vector in (11), we have
(12)f(α)={∑i=1nt∑j=1m(d¨(α,t)-d¨*(α,t))∂d¨(α,t)∂xc=0,∑i=1nt∑j=1m(d¨(α,t)-d¨*(α,t))∂d¨(α,t)∂hc=0,
where∂d¨(α,t)/∂xcand∂d¨(α,t)/∂hcare the acceleration response sensitivities with respect to the crack location and depth, respectively.


### 3.2. Homotopy Iteration Algorithm for Crack Identification

The basic idea of homotopy iteration algorithm is to introduce a parameterλinto the nonlinear equationf(α)=0to construct a family of homotopy mappingH(α,λ), whereλ∈[0,1]is named the homotopy parameter. The homotopy equation is assumed to beH(α,λ)=λf(α)+(1-λ)f0(α)in this paper. Whenλ=1,His the mapping off(α)as (12) shows, and, whenλ=0,His assumed to be the mapping off0(α). Andf0(α)is set to be similar tof(α)but, related to the initial point, can be expressed asf0(α)=(d¨(α,t)-d¨(α0,t))∂d¨(α,t)/∂α, whereα0is the arbitrarily selected initial value vector. In this way, a family of homotopy mappingHcan be constructed instead of a single mappingf. Thus, the homotopy equation can meet the following equation:
(13)H(α,λ)=λf(α)+(1-λ)f0(α)=0,(14)H(α,0)=f0(α)=0,(15)H(α,1)=f(α)=0.
Equations (14) and (15) show that as the homotopy parameterλincreases from zero to one, the homotopy equationH(α,λ)varies fromf0(α)tof(α). Hence, the parameters of the crack can be obtained by following the variations of parameterλ, and, whenH(α,λ)=0, the identification vectorαis in the homotopy path.After obtaining the homotopy equation, the Newton iterative method is used to track the path in the process of identification. The solving procedure is explained as follows.Step 1.
Guess an initial value vector of crack parametersα0=[xc0,hc0], divide the range of homotopy parameterλ∈[0,1]intonequal parts, and then obtaindλ=1/n, letλ0=0.Step 2.
Trace the homotopy path from the initial value(α0,λ0)and give an incrementdλto the homotopy parameterλasλ1=λ0+dλ.Step 3.
Using Newton iterative method to calculate an updated vector of crack parameterα1=[xc1,hc1]in the homotopy path, we have
(16)α1=α0-(∂H(α,λ1)∂α|α=α0)-1H(α0,λ1),
where∂H(α,λ1)/∂αis the homotopy equation with respect to the crack parameters. The value ofH(α0,λ1)and∂H(α,λ1)/∂αwill be introduced in next section.Step 4.
Check whether the updated crack parameter vectorα1=[xc1,hc1]is physically meaningful or not. If not, the result is considered diverged, then stop iteration and move to Step 1 to start with a new initial value or division of homotopy parameter. If otherwise, go to the next step.Step 5.
Letα0=α1; repeat from Steps 3 to 4, until the obtained(α1,λ1)in homotopy path satisfies the following convergence conditions:
(17a)∥α1-α0α1∥≤Tol1,
whereTol1is the first tolerance value for convergence and is taken as1.0×10-3in this study.Step 6.
Repeat the Newton iterative again, until the following convergence criterion is satisfied. That is, the identification parameters are assumed to in the homotopy path whenH(α1,λ1)meets the following criterion:
(17b)∥H(α1,λ1)∥≤Tol2,
whereTol2is the second tolerance value and is taken to be 0.1 for all the study cases.Step 7.
The pair of unknowns(α1,λ1)is taken as the new set of initial value(α0,λ0)and repeat Steps 2 to 5 until the homotopy parameterλreaches 1.0. The finalα1would then be the required set of identification results.


### 3.3. Homotopy Equation and Its First Derivative

In order to obtain the homotopy equation, the acceleration response with respect to the crack parameter should be obtained first. We take the derivative on both sides of (9) with respect to the crack location and crack depth; then we have
(18)M∂Q¨∂xc+C∂Q˙∂xc+K∂Q∂xc=-a2∂K∂xcQ˙-∂K∂xcQ,M∂Q¨∂hc+C∂Q˙∂hc+K∂Q∂hc=-a2∂K∂hcQ˙-∂K∂hcQ,
where∂Q¨/∂xc,∂Q¨/∂hc,∂Q˙/∂xc,∂Q˙/∂hc,∂Q/∂xc, and∂Q/∂hcare the vectors of generalized acceleration, velocity, and displacement sensitivities with respect to the crack location and depth, respectively.∂K/∂xcand∂K/∂hcare the first partial derivatives of the stiffness matrix with respect to the crack parameters. Since the flexural rigidityEI(x)involves the crack location and depth and the global stiffness matrixKcan be obtained from (7), the derivatives∂K/∂xcand∂K/∂hccan be obtained directly. As the dynamic responsesQandQ˙have been calculated from (9), the right-hand side of (18) can be considered as a form of external force. Thus, the dynamic response sensitivity (i.e., the generalized acceleration response sensitivity, velocity response sensitivity, and displacement response sensitivity) can also be calculated from the Newmark integration method. Then, the physical acceleration response sensitivity can be obtained as(19a)∂d¨∂xc=S(x)∂Q¨∂xc,(19b)∂d¨∂hc=S(x)∂Q¨∂hc.Furthermore, the objective function in (12) and homotopy equation in (13) can then be obtained.To obtain the first derivative of the homotopy equation with respect to the crack parameters, we should apply the first derivative to both sides of (18) with respect to the crack parameters as
(20)M∂2Q¨∂xc2+C∂2Q˙∂xc2+K∂2Q∂xc2=-a2∂2K∂xc2Q˙-∂2K∂xc2Q-2(a2∂K∂xc∂Q˙∂xc+∂K∂xc∂Q∂xc),M∂2Q¨∂xc∂hc+C∂2Q˙∂xc∂hc+K∂2Q∂xc∂hc=-a2∂2K∂xc∂hcQ˙-∂2K∂xc∂hcQ-a2∂K∂hc∂Q˙∂xc-∂K∂hc∂Q∂xc-a2∂K∂xc∂Q˙∂hc-∂K∂xc∂Q∂hc,M∂2Q¨∂hc2+C∂2Q˙∂hc2+K∂2Q∂hc2=-a2∂2K∂hc2Q˙-∂2K∂hc2Q-2(a2∂K∂hc∂Q˙∂hc+∂K∂hc∂Q∂hc),M∂2Q¨∂hc∂xc+C∂2Q˙∂hc∂xc+K∂2Q∂hc∂xc=-a2∂2K∂hc∂xcQ˙-∂2K∂hc∂xcQ-a2∂K∂xc∂Q˙∂hc-∂K∂xc∂Q∂hc-a2∂K∂hc∂Q˙∂xc-∂K∂hc∂Q∂xc,
where∂2K/∂xc2,∂2K/∂xc∂hc,∂2K/∂hc2, and∂2K/∂hc∂xcare the second partial derivative of the stiffness matrix with respect to the crack parameters, and∂2K/∂xc∂hc=∂2K/∂hc∂xc. Since the response sensitivities∂Q/∂hc,∂Q˙/∂xc, and∂Q˙/∂hchave been calculated from (18), the second partial derivative of generalized acceleration response with respect to the unknown variables(xc,hc)can be calculated by the Newmark integration with (20). Then, the physical acceleration response sensitivity can be obtained as(21a)∂2d¨∂xc2=S(x)∂2Q¨∂xc2,(21b)∂2d¨∂hc2=S(x)∂2Q¨∂hc2,(21c)∂2d¨∂xc∂hc=∂2d¨∂hc∂xc=S(x)∂2Q¨∂xc∂hc.Furthermore, the first derivative of the homotopy mapping in (13) with respect to the unknown vector∂H(α,λ)/∂αcan be obtained as
(22)∂H(α,λ)∂α=[λ∂f(α)∂xc+(1-λ)∂f0(α)∂xcλ∂f(α)∂hc+(1-λ)∂f0(α)∂hc],
where
(23)∂f(α)∂xc=[∑i=1nt∑j=1m((d¨-d¨*)∂2d¨∂xc2+(∂d¨∂xc)2),hhh∑i=1nt∑j=1m((d¨-d¨*)∂2d¨∂hc∂xc+∂d¨∂xc∂d¨∂hc)],(24)∂f(α)∂hc=[∑i=1nt∑j=1m((d¨-d¨*)∂2d¨∂xc∂hc+∂d¨∂hc∂d¨∂xc),hhh∑i=1nt∑j=1m((d¨-d¨*)∂2d¨∂hc2+(∂d¨∂hc)2)].
The first derivatives∂f0(α)/∂xcand∂f0(α)/∂hccan also be obtained in a similar way.


## 3.1. Objective Function

The problem of crack identification can be treated as an optimization problem. In the identification, the unknowns consist of the crack locationxcand crack depthhc; that is,α=[xc,hc]. If there aremmeasured points in the beam structure, the measured acceleration response can be expressed asd¨*=(d¨1,d¨2,…,d¨m). The objective function for the optimization problem can be derived by minimizing the error between the calculated acceleration response and the simulated measured one as
(11)J=12∑i=1nt∑j=1m(d¨(α,t)-d¨*(α,t))2,i=1,2,…,nt,j=1,2,…,m,
wheretdenotes the the measured time instants,nt=t/dtis the number of time step, anddtis the increment of time step.The inverse problem can be solved by minimizing the objective function. Taking partial derivative with respect to the identification vector in (11), we have
(12)f(α)={∑i=1nt∑j=1m(d¨(α,t)-d¨*(α,t))∂d¨(α,t)∂xc=0,∑i=1nt∑j=1m(d¨(α,t)-d¨*(α,t))∂d¨(α,t)∂hc=0,
where∂d¨(α,t)/∂xcand∂d¨(α,t)/∂hcare the acceleration response sensitivities with respect to the crack location and depth, respectively.


## 3.2. Homotopy Iteration Algorithm for Crack Identification

The basic idea of homotopy iteration algorithm is to introduce a parameterλinto the nonlinear equationf(α)=0to construct a family of homotopy mappingH(α,λ), whereλ∈[0,1]is named the homotopy parameter. The homotopy equation is assumed to beH(α,λ)=λf(α)+(1-λ)f0(α)in this paper. Whenλ=1,His the mapping off(α)as (12) shows, and, whenλ=0,His assumed to be the mapping off0(α). Andf0(α)is set to be similar tof(α)but, related to the initial point, can be expressed asf0(α)=(d¨(α,t)-d¨(α0,t))∂d¨(α,t)/∂α, whereα0is the arbitrarily selected initial value vector. In this way, a family of homotopy mappingHcan be constructed instead of a single mappingf. Thus, the homotopy equation can meet the following equation:
(13)H(α,λ)=λf(α)+(1-λ)f0(α)=0,(14)H(α,0)=f0(α)=0,(15)H(α,1)=f(α)=0.
Equations (14) and (15) show that as the homotopy parameterλincreases from zero to one, the homotopy equationH(α,λ)varies fromf0(α)tof(α). Hence, the parameters of the crack can be obtained by following the variations of parameterλ, and, whenH(α,λ)=0, the identification vectorαis in the homotopy path.After obtaining the homotopy equation, the Newton iterative method is used to track the path in the process of identification. The solving procedure is explained as follows.Step 1.
Guess an initial value vector of crack parametersα0=[xc0,hc0], divide the range of homotopy parameterλ∈[0,1]intonequal parts, and then obtaindλ=1/n, letλ0=0.Step 2.
Trace the homotopy path from the initial value(α0,λ0)and give an incrementdλto the homotopy parameterλasλ1=λ0+dλ.Step 3.
Using Newton iterative method to calculate an updated vector of crack parameterα1=[xc1,hc1]in the homotopy path, we have
(16)α1=α0-(∂H(α,λ1)∂α|α=α0)-1H(α0,λ1),
where∂H(α,λ1)/∂αis the homotopy equation with respect to the crack parameters. The value ofH(α0,λ1)and∂H(α,λ1)/∂αwill be introduced in next section.Step 4.
Check whether the updated crack parameter vectorα1=[xc1,hc1]is physically meaningful or not. If not, the result is considered diverged, then stop iteration and move to Step 1 to start with a new initial value or division of homotopy parameter. If otherwise, go to the next step.Step 5.
Letα0=α1; repeat from Steps 3 to 4, until the obtained(α1,λ1)in homotopy path satisfies the following convergence conditions:
(17a)∥α1-α0α1∥≤Tol1,
whereTol1is the first tolerance value for convergence and is taken as1.0×10-3in this study.Step 6.
Repeat the Newton iterative again, until the following convergence criterion is satisfied. That is, the identification parameters are assumed to in the homotopy path whenH(α1,λ1)meets the following criterion:
(17b)∥H(α1,λ1)∥≤Tol2,
whereTol2is the second tolerance value and is taken to be 0.1 for all the study cases.Step 7.
The pair of unknowns(α1,λ1)is taken as the new set of initial value(α0,λ0)and repeat Steps 2 to 5 until the homotopy parameterλreaches 1.0. The finalα1would then be the required set of identification results.


## 3.3. Homotopy Equation and Its First Derivative

In order to obtain the homotopy equation, the acceleration response with respect to the crack parameter should be obtained first. We take the derivative on both sides of (9) with respect to the crack location and crack depth; then we have
(18)M∂Q¨∂xc+C∂Q˙∂xc+K∂Q∂xc=-a2∂K∂xcQ˙-∂K∂xcQ,M∂Q¨∂hc+C∂Q˙∂hc+K∂Q∂hc=-a2∂K∂hcQ˙-∂K∂hcQ,
where∂Q¨/∂xc,∂Q¨/∂hc,∂Q˙/∂xc,∂Q˙/∂hc,∂Q/∂xc, and∂Q/∂hcare the vectors of generalized acceleration, velocity, and displacement sensitivities with respect to the crack location and depth, respectively.∂K/∂xcand∂K/∂hcare the first partial derivatives of the stiffness matrix with respect to the crack parameters. Since the flexural rigidityEI(x)involves the crack location and depth and the global stiffness matrixKcan be obtained from (7), the derivatives∂K/∂xcand∂K/∂hccan be obtained directly. As the dynamic responsesQandQ˙have been calculated from (9), the right-hand side of (18) can be considered as a form of external force. Thus, the dynamic response sensitivity (i.e., the generalized acceleration response sensitivity, velocity response sensitivity, and displacement response sensitivity) can also be calculated from the Newmark integration method. Then, the physical acceleration response sensitivity can be obtained as(19a)∂d¨∂xc=S(x)∂Q¨∂xc,(19b)∂d¨∂hc=S(x)∂Q¨∂hc.Furthermore, the objective function in (12) and homotopy equation in (13) can then be obtained.To obtain the first derivative of the homotopy equation with respect to the crack parameters, we should apply the first derivative to both sides of (18) with respect to the crack parameters as
(20)M∂2Q¨∂xc2+C∂2Q˙∂xc2+K∂2Q∂xc2=-a2∂2K∂xc2Q˙-∂2K∂xc2Q-2(a2∂K∂xc∂Q˙∂xc+∂K∂xc∂Q∂xc),M∂2Q¨∂xc∂hc+C∂2Q˙∂xc∂hc+K∂2Q∂xc∂hc=-a2∂2K∂xc∂hcQ˙-∂2K∂xc∂hcQ-a2∂K∂hc∂Q˙∂xc-∂K∂hc∂Q∂xc-a2∂K∂xc∂Q˙∂hc-∂K∂xc∂Q∂hc,M∂2Q¨∂hc2+C∂2Q˙∂hc2+K∂2Q∂hc2=-a2∂2K∂hc2Q˙-∂2K∂hc2Q-2(a2∂K∂hc∂Q˙∂hc+∂K∂hc∂Q∂hc),M∂2Q¨∂hc∂xc+C∂2Q˙∂hc∂xc+K∂2Q∂hc∂xc=-a2∂2K∂hc∂xcQ˙-∂2K∂hc∂xcQ-a2∂K∂xc∂Q˙∂hc-∂K∂xc∂Q∂hc-a2∂K∂hc∂Q˙∂xc-∂K∂hc∂Q∂xc,
where∂2K/∂xc2,∂2K/∂xc∂hc,∂2K/∂hc2, and∂2K/∂hc∂xcare the second partial derivative of the stiffness matrix with respect to the crack parameters, and∂2K/∂xc∂hc=∂2K/∂hc∂xc. Since the response sensitivities∂Q/∂hc,∂Q˙/∂xc, and∂Q˙/∂hchave been calculated from (18), the second partial derivative of generalized acceleration response with respect to the unknown variables(xc,hc)can be calculated by the Newmark integration with (20). Then, the physical acceleration response sensitivity can be obtained as(21a)∂2d¨∂xc2=S(x)∂2Q¨∂xc2,(21b)∂2d¨∂hc2=S(x)∂2Q¨∂hc2,(21c)∂2d¨∂xc∂hc=∂2d¨∂hc∂xc=S(x)∂2Q¨∂xc∂hc.Furthermore, the first derivative of the homotopy mapping in (13) with respect to the unknown vector∂H(α,λ)/∂αcan be obtained as
(22)∂H(α,λ)∂α=[λ∂f(α)∂xc+(1-λ)∂f0(α)∂xcλ∂f(α)∂hc+(1-λ)∂f0(α)∂hc],
where
(23)∂f(α)∂xc=[∑i=1nt∑j=1m((d¨-d¨*)∂2d¨∂xc2+(∂d¨∂xc)2),hhh∑i=1nt∑j=1m((d¨-d¨*)∂2d¨∂hc∂xc+∂d¨∂xc∂d¨∂hc)],(24)∂f(α)∂hc=[∑i=1nt∑j=1m((d¨-d¨*)∂2d¨∂xc∂hc+∂d¨∂hc∂d¨∂xc),hhh∑i=1nt∑j=1m((d¨-d¨*)∂2d¨∂hc2+(∂d¨∂hc)2)].
The first derivatives∂f0(α)/∂xcand∂f0(α)/∂hccan also be obtained in a similar way.


## 4. Numerical Simulations

### 4.1. A Simply Supported Beam

A uniform cross section simply supported beam with cracks is studied as shown in Figure1. The parameters of the beam are: Young’s modulusE=69.79 GPa, mass densityρ=2600 kg/m3, lengthL=2000 mm, widthw=50 mm, and heighth=25 mm. An impulsive force acts at the 700 mm from the left support in the negativey-direction fromt=0 s tot=0.04 s with
(25)g(t)={10000tN,(0s≤t≤0.02s),10000(0.04-t)N,(0.02s<t≤0.04s).
In calculating the dynamic response, the time step is 0.001 s. Artificial measurement noise with different levels is added to the calculated responses to simulate the “measured” structural responses. All the time history data of 1.0 second are used in the identification unless otherwise specified. The two damping coefficients used for calculating Rayleigh damping matrix are both assumed to be 0.01.Using the traditional FEM, the first eight frequencies of the beam are 14.65, 58.73, 131.85, 234.95, 366.31, 528.84, 718.45, and 941.15 HZ. Comparison with Table1, one can find that when the number of term in CT isn, only the firstnfrequency is convergence for the uniform cross section beam. When only 10 modes are used, the first 8 natural frequencies can be obtained with good accuracy (the max error for the 8th natural frequency is only 0.17%). Indeed, when more terms of the harmonic functions are used, the accuracy can further be improved, but more computational time is needed in the calculation. When dynamic loads act on the beam, the dynamic responses of the beam will only include the first few lower modes; generally speaking, the contribution of the higher modes on the dynamic responses can be neglected. Thus, the number of the mode functions used from the CT is set to be 10 in the following study.Table 1
Comparison of the convergence property with different number of terms in CT.
Number of term in CT
Natural frequencies in Hz (xc=1000 mm, hc=1 mm)
1
2
3
4
5
6
7
8
4
14.65
58.73
131.86
234.92
592.11
890.33
53510.81
58925.59
5
14.65
58.73
131.86
234.92
366.28
890.33
1252.84
58925.59
6
14.65
58.73
131.86
234.92
366.28
528.55
1252.83
1680.70
9
14.65
58.73
131.86
234.92
366.28
528.56
717.93
939.60
10
14.65
58.73
131.86
234.92
366.28
528.56
717.93
939.60


#### 4.1.1. Effect of Initial Value on Crack Identification

In this case, single crack identification is conducted to illustrate the proposed method. Three measurement points located at 600 mm, 1100 mm, and 1600 mm from the left support are used in the crack identification. Totally, four Scenarios are studied as listed in Table2. One can find that the crack parameters have been identified successfully with high accuracy in all Scenarios. It should be pointed out as the initial crack location is far away from the true crack location, the crack can also be identified successfully but the number of iteration will increase accordingly. And in general, finer division in the homotopy parameterλshould be considered as to track the correct path. Scenarios 3 and 4 show that the method is less affected by the initial values.Table 2
Effect of the initial value on crack identification.
Scenarios
Initial value (mm)
True value (mm)
Identified results (mm)
Division ofλ
Measurement points (mm)
No. of iteration
x
c
0
h
c
0
x
c
h
c
x
c
h
c
1
1500
2
1720
10
1720.000000
10.000000
5
600, 1100, 1600
31
2
1500
2
1000
10
999.999990
10.000000
5
600, 1100, 1600
38
3
1500
2
1000
1
1000.000234
1.000000
3
600, 1100, 1600
26
4
260
3
1000
1
999.999997
1.000000
21
600, 1100, 1600
146


#### 4.1.2. Effect of Measurement Time Duration

In this case, the effect of measurement time duration in crack identification is studied. The identification of Scenario 3 in Table2 is reexamined; except the measurement time duration is taken to be 1.0 s, 2.0 s, and 4.0 s, respectively. The crack parameters are identified with satisfactory accuracy as shown in Table 3. It also can be notes that longer time duration has little effect on the degree of accuracy but will increase the number of iteration. And the process of iteration for this study is shown in Figure 3; one can see that longer measurement time will obtain more stable convergence pattern when tracking the approximate values in the homotopy path.Table 3
Effect of measurement time.
Scenarios
Identified results (mm)
Division ofλ
Measurement time (s)
No. of iteration
x
c
h
c
1
1000.000234
1.000000
3
1
26
2
1000.000237
1.000000
3
2
30
3
1000.000187
1.000000
3
4
32Effect of measurement time duration.
(a)
(b)


#### 4.1.3. Effect of Measurement Points

In this case, how the measurement points affect the accuracy and the iterative process is studied. Three, four, and six measurement points are used for crack identification to give a comparison. And the other parameters are the same as the last case. The identification results are listed in Table4. It can be found that increasing measurement points can improve the accuracy of the identification results. Figure 4 gives a comparison on the identified process of iteration. It can be noted that different number of measurement points has the similar convergence process. And the proposed method does not need a large number of measurement points.Table 4
Effect of measurement points.
Scenarios
Identified results (mm)
Measurement points (mm)
No. of iteration
x
c
h
c
1
1000.000234
1.000000
600, 1100, 1600
26
2
1000.000005
1.000000
400, 800, 1200, 1600
25
3
999.999999
1.000000
300, 600, 900, 1200, 1500, 1800
25Effect of measurement points.
(a)
(b)


#### 4.1.4. Effect of Division of Homotopy Parameter

In this study, the effect of division of homotopy parameter on the identified results is discussed. Other parameters are the same as those in the last study except that the homotopy parameter is divided into 3, 4, and 5 parts. And the identified results for the study are listed in Table5. We can find that increase of division parts of the homotopy parameter, the number of iteration will increase accordingly in general but it has little effect on the accuracy of identified results. Figure 5 shows the evolution of the crack parameters in the process of iteration for different division parts of the homotopy parameter; one can see that the convergence of Newton’s method to track the homotopy path is similar.Table 5
Effect of division of homotopy parameter.
Scenarios
Identified results (mm)
Division ofλ
Measurement time (s)
No. of iteration
x
c
h
c
1
1000.000234
1.000000
3
1
26
2
1000.000003
1.000000
4
1
28
3
999.9999346
1.000000
5
1
30Effect of division of homotopy parameter.
(a)
(b)


#### 4.1.5. Effect of Measurement Noise

In this section, the effect of measurement noise on the accuracy of identified results is taken into account. Again, the identification of Scenario 3 in Table2 is studied. The effect of measurement noise is simulated as a normally distributed random error with zero mean and a unit standard deviation is added to the calculated acceleration as
(26)d¨^=d¨+Ep×Noise×var(d¨),
whered¨^is the vectors of measured structural acceleration response,EPis the noise level,Noiseis standard normal distribution vector with zero mean and unit standard deviation, andvar(·)is the variance of the time history.The relative errors of crack location and depth identification are defined as(27)Relativeerror=∥αid-αtrueαtrue∥×100%,
whereαidandαtrueare the identified and true values, respectively. Table 6 gives a comparison of the identified results for 0%, 1%, 5%, and 15% noise level, respectively. This study shows that the crack parameters have been identified successfully even with 15% measurement noise. With the increase of the measurement noise level, identification errors will become larger and the max identification error is 8.93% in the location of crack and 5.24% in the crack depth.Table 6
Effect of measurement noise.
Scenarios
Noise
Identified results (mm)/rel. (%)
Division ofλ
No. of iteration
x
c
h
c
1
Nil
1000.000234
/
0.00
1.000000
/
0.00
3
26
2
1%
1021.498303
/
2.15
1.003232
/
0.32
3
22
3
5%
1070.401667
/
7.04
1.022658
/
2.27
3
23
4
15%
1089.318837
/
8.93
1.052358
/
5.24
3
19
Note: “Rel. (%)” denotes relative error.


#### 4.1.6. Identification of Multiple Cracks

The same simply supported beam is studied with two cracks in this Section. The cracks are assumed to be located at 450 mm and 1000 mm from the left support of the beam, respectively. The depths of the crack are 5 mm and 1 mm, respectively. In the identification, the initial values of locations for these two cracks are set to be 500 mm and 1500 mm, respectively, and depths are both set to be 2 mm. Six measurement points as the Scenario 3 in Table4 are used in this case. And the setting of other parameters is shown in Table 7. The identification results and the relative errors for the cases with 0% and 5% noise levels are considered. It can be noted that the identified results for the two cracks are as good as those for a single crack when the measurement noise is free. Comparing the identified results with measurement noise and the free noise one, we can find that the effect of measurement noise on crack identification is noticeable, but the identified results are still satisfactory. For 5% noise level, the maximum relative errors are 1.82% and 0.28% in the crack location and depth for the first crack, respectively and 6.49% and 7.22% for the second crack, respectively. These results further illustrate the effectiveness of the proposed method.Table 7
Identification of two cracks in a simply supported beam.
Scenarios
Noise
Identified results (mm)/rel. (%) in crack 1
Identified results (mm)/rel. (%) in crack 2
Division ofλ
Measurement time (s)
No. of iteration
x
c
1
h
c
1
x
c
2
h
c
2
1
Nil
450.000000
/
0.00
5.000000
/
0.00
999.999997
/
0.00
1.000000
/
0.00
5
1
78
2
5%
441.804646
/
1.82
4.986154
/
0.28
1064.867252
/
6.49
1.072199
/
7.22
7
1
74


### 4.2. A Cantilever Beam

A uniform cross section cantilever beam with cracks is studied as another example. As Figure6 shows, the physical parameters of the beam are: elastic modulus of materialE=69.79 GPa, mass densityρ=2600 kg/m3, lengthL=1000 mm, widthw=50 mm and heighth=25 mm. when calculating the dynamic response, a sinusoidal forceg(t)=40sin(15πt)Nacts at the free end of the beam. Four measurement points located at 300 mm, 450 mm, 600 mm, and 750 mm from the clamped end are used to record the acceleration data. The two damping coefficients used for calculating Rayleigh damping matrix are both assumed to be 0.01. The time step is 0.005 s and the time duration is taken to be 2.0 s and 4.0 s, respectively. The number of the mode functions used from the CT is also set to be 10, and the way to choose the number of term is the same as the simply supported beam. Both single crack and two cracks are simulated in this case. Good identification results are listed in Table 8. For a single crack identification, the measurement of time duration may have little effect on the degree of accuracy as Section 4.1.2 shows. But for two cracks, it can improve the accuracy of the identified results in general as the unknown parameters increases in the identification. This case further illustrates the effectiveness of the proposed method.Table 8
Crack identification in a cantilever beam.
Initial
True
Identified results
2s
4s
Single crack
Crack locationxc (mm)
750
250
250.000162
250.000259
Crack depthhc (mm)
3
1
1.000001
1.000001
Division ofλ
—
—
26
20
No. of iteration
—
—
92
96
Two cracks
Crack locationsxc1/xc2 (mm)
360
/
700
250
/
600
249.603356
/
599.789017
250.009037
/
600.005717
Crack depthshc1/hc2 (mm)
2
/
7.8
1
/
6.5
0.997240
/
6.498047
1.000067
/
6.500058
Division ofλ
—
—
16
15
No. of iteration
—
—
137
137Figure 6
A cantilever beam with cracks.


## 4.1. A Simply Supported Beam

A uniform cross section simply supported beam with cracks is studied as shown in Figure1. The parameters of the beam are: Young’s modulusE=69.79 GPa, mass densityρ=2600 kg/m3, lengthL=2000 mm, widthw=50 mm, and heighth=25 mm. An impulsive force acts at the 700 mm from the left support in the negativey-direction fromt=0 s tot=0.04 s with
(25)g(t)={10000tN,(0s≤t≤0.02s),10000(0.04-t)N,(0.02s<t≤0.04s).
In calculating the dynamic response, the time step is 0.001 s. Artificial measurement noise with different levels is added to the calculated responses to simulate the “measured” structural responses. All the time history data of 1.0 second are used in the identification unless otherwise specified. The two damping coefficients used for calculating Rayleigh damping matrix are both assumed to be 0.01.Using the traditional FEM, the first eight frequencies of the beam are 14.65, 58.73, 131.85, 234.95, 366.31, 528.84, 718.45, and 941.15 HZ. Comparison with Table1, one can find that when the number of term in CT isn, only the firstnfrequency is convergence for the uniform cross section beam. When only 10 modes are used, the first 8 natural frequencies can be obtained with good accuracy (the max error for the 8th natural frequency is only 0.17%). Indeed, when more terms of the harmonic functions are used, the accuracy can further be improved, but more computational time is needed in the calculation. When dynamic loads act on the beam, the dynamic responses of the beam will only include the first few lower modes; generally speaking, the contribution of the higher modes on the dynamic responses can be neglected. Thus, the number of the mode functions used from the CT is set to be 10 in the following study.Table 1
Comparison of the convergence property with different number of terms in CT.
Number of term in CT
Natural frequencies in Hz (xc=1000 mm, hc=1 mm)
1
2
3
4
5
6
7
8
4
14.65
58.73
131.86
234.92
592.11
890.33
53510.81
58925.59
5
14.65
58.73
131.86
234.92
366.28
890.33
1252.84
58925.59
6
14.65
58.73
131.86
234.92
366.28
528.55
1252.83
1680.70
9
14.65
58.73
131.86
234.92
366.28
528.56
717.93
939.60
10
14.65
58.73
131.86
234.92
366.28
528.56
717.93
939.60


### 4.1.1. Effect of Initial Value on Crack Identification

In this case, single crack identification is conducted to illustrate the proposed method. Three measurement points located at 600 mm, 1100 mm, and 1600 mm from the left support are used in the crack identification. Totally, four Scenarios are studied as listed in Table2. One can find that the crack parameters have been identified successfully with high accuracy in all Scenarios. It should be pointed out as the initial crack location is far away from the true crack location, the crack can also be identified successfully but the number of iteration will increase accordingly. And in general, finer division in the homotopy parameterλshould be considered as to track the correct path. Scenarios 3 and 4 show that the method is less affected by the initial values.Table 2
Effect of the initial value on crack identification.
Scenarios
Initial value (mm)
True value (mm)
Identified results (mm)
Division ofλ
Measurement points (mm)
No. of iteration
x
c
0
h
c
0
x
c
h
c
x
c
h
c
1
1500
2
1720
10
1720.000000
10.000000
5
600, 1100, 1600
31
2
1500
2
1000
10
999.999990
10.000000
5
600, 1100, 1600
38
3
1500
2
1000
1
1000.000234
1.000000
3
600, 1100, 1600
26
4
260
3
1000
1
999.999997
1.000000
21
600, 1100, 1600
146


### 4.1.2. Effect of Measurement Time Duration

In this case, the effect of measurement time duration in crack identification is studied. The identification of Scenario 3 in Table2 is reexamined; except the measurement time duration is taken to be 1.0 s, 2.0 s, and 4.0 s, respectively. The crack parameters are identified with satisfactory accuracy as shown in Table 3. It also can be notes that longer time duration has little effect on the degree of accuracy but will increase the number of iteration. And the process of iteration for this study is shown in Figure 3; one can see that longer measurement time will obtain more stable convergence pattern when tracking the approximate values in the homotopy path.Table 3
Effect of measurement time.
Scenarios
Identified results (mm)
Division ofλ
Measurement time (s)
No. of iteration
x
c
h
c
1
1000.000234
1.000000
3
1
26
2
1000.000237
1.000000
3
2
30
3
1000.000187
1.000000
3
4
32Effect of measurement time duration.
(a)
(b)


### 4.1.3. Effect of Measurement Points

In this case, how the measurement points affect the accuracy and the iterative process is studied. Three, four, and six measurement points are used for crack identification to give a comparison. And the other parameters are the same as the last case. The identification results are listed in Table4. It can be found that increasing measurement points can improve the accuracy of the identification results. Figure 4 gives a comparison on the identified process of iteration. It can be noted that different number of measurement points has the similar convergence process. And the proposed method does not need a large number of measurement points.Table 4
Effect of measurement points.
Scenarios
Identified results (mm)
Measurement points (mm)
No. of iteration
x
c
h
c
1
1000.000234
1.000000
600, 1100, 1600
26
2
1000.000005
1.000000
400, 800, 1200, 1600
25
3
999.999999
1.000000
300, 600, 900, 1200, 1500, 1800
25Effect of measurement points.
(a)
(b)


### 4.1.4. Effect of Division of Homotopy Parameter

In this study, the effect of division of homotopy parameter on the identified results is discussed. Other parameters are the same as those in the last study except that the homotopy parameter is divided into 3, 4, and 5 parts. And the identified results for the study are listed in Table5. We can find that increase of division parts of the homotopy parameter, the number of iteration will increase accordingly in general but it has little effect on the accuracy of identified results. Figure 5 shows the evolution of the crack parameters in the process of iteration for different division parts of the homotopy parameter; one can see that the convergence of Newton’s method to track the homotopy path is similar.Table 5
Effect of division of homotopy parameter.
Scenarios
Identified results (mm)
Division ofλ
Measurement time (s)
No. of iteration
x
c
h
c
1
1000.000234
1.000000
3
1
26
2
1000.000003
1.000000
4
1
28
3
999.9999346
1.000000
5
1
30Effect of division of homotopy parameter.
(a)
(b)


### 4.1.5. Effect of Measurement Noise

In this section, the effect of measurement noise on the accuracy of identified results is taken into account. Again, the identification of Scenario 3 in Table2 is studied. The effect of measurement noise is simulated as a normally distributed random error with zero mean and a unit standard deviation is added to the calculated acceleration as
(26)d¨^=d¨+Ep×Noise×var(d¨),
whered¨^is the vectors of measured structural acceleration response,EPis the noise level,Noiseis standard normal distribution vector with zero mean and unit standard deviation, andvar(·)is the variance of the time history.The relative errors of crack location and depth identification are defined as(27)Relativeerror=∥αid-αtrueαtrue∥×100%,
whereαidandαtrueare the identified and true values, respectively. Table 6 gives a comparison of the identified results for 0%, 1%, 5%, and 15% noise level, respectively. This study shows that the crack parameters have been identified successfully even with 15% measurement noise. With the increase of the measurement noise level, identification errors will become larger and the max identification error is 8.93% in the location of crack and 5.24% in the crack depth.Table 6
Effect of measurement noise.
Scenarios
Noise
Identified results (mm)/rel. (%)
Division ofλ
No. of iteration
x
c
h
c
1
Nil
1000.000234
/
0.00
1.000000
/
0.00
3
26
2
1%
1021.498303
/
2.15
1.003232
/
0.32
3
22
3
5%
1070.401667
/
7.04
1.022658
/
2.27
3
23
4
15%
1089.318837
/
8.93
1.052358
/
5.24
3
19
Note: “Rel. (%)” denotes relative error.


### 4.1.6. Identification of Multiple Cracks

The same simply supported beam is studied with two cracks in this Section. The cracks are assumed to be located at 450 mm and 1000 mm from the left support of the beam, respectively. The depths of the crack are 5 mm and 1 mm, respectively. In the identification, the initial values of locations for these two cracks are set to be 500 mm and 1500 mm, respectively, and depths are both set to be 2 mm. Six measurement points as the Scenario 3 in Table4 are used in this case. And the setting of other parameters is shown in Table 7. The identification results and the relative errors for the cases with 0% and 5% noise levels are considered. It can be noted that the identified results for the two cracks are as good as those for a single crack when the measurement noise is free. Comparing the identified results with measurement noise and the free noise one, we can find that the effect of measurement noise on crack identification is noticeable, but the identified results are still satisfactory. For 5% noise level, the maximum relative errors are 1.82% and 0.28% in the crack location and depth for the first crack, respectively and 6.49% and 7.22% for the second crack, respectively. These results further illustrate the effectiveness of the proposed method.Table 7
Identification of two cracks in a simply supported beam.
Scenarios
Noise
Identified results (mm)/rel. (%) in crack 1
Identified results (mm)/rel. (%) in crack 2
Division ofλ
Measurement time (s)
No. of iteration
x
c
1
h
c
1
x
c
2
h
c
2
1
Nil
450.000000
/
0.00
5.000000
/
0.00
999.999997
/
0.00
1.000000
/
0.00
5
1
78
2
5%
441.804646
/
1.82
4.986154
/
0.28
1064.867252
/
6.49
1.072199
/
7.22
7
1
74


## 4.1.1. Effect of Initial Value on Crack Identification

In this case, single crack identification is conducted to illustrate the proposed method. Three measurement points located at 600 mm, 1100 mm, and 1600 mm from the left support are used in the crack identification. Totally, four Scenarios are studied as listed in Table2. One can find that the crack parameters have been identified successfully with high accuracy in all Scenarios. It should be pointed out as the initial crack location is far away from the true crack location, the crack can also be identified successfully but the number of iteration will increase accordingly. And in general, finer division in the homotopy parameterλshould be considered as to track the correct path. Scenarios 3 and 4 show that the method is less affected by the initial values.Table 2
Effect of the initial value on crack identification.
Scenarios
Initial value (mm)
True value (mm)
Identified results (mm)
Division ofλ
Measurement points (mm)
No. of iteration
x
c
0
h
c
0
x
c
h
c
x
c
h
c
1
1500
2
1720
10
1720.000000
10.000000
5
600, 1100, 1600
31
2
1500
2
1000
10
999.999990
10.000000
5
600, 1100, 1600
38
3
1500
2
1000
1
1000.000234
1.000000
3
600, 1100, 1600
26
4
260
3
1000
1
999.999997
1.000000
21
600, 1100, 1600
146


## 4.1.2. Effect of Measurement Time Duration

In this case, the effect of measurement time duration in crack identification is studied. The identification of Scenario 3 in Table2 is reexamined; except the measurement time duration is taken to be 1.0 s, 2.0 s, and 4.0 s, respectively. The crack parameters are identified with satisfactory accuracy as shown in Table 3. It also can be notes that longer time duration has little effect on the degree of accuracy but will increase the number of iteration. And the process of iteration for this study is shown in Figure 3; one can see that longer measurement time will obtain more stable convergence pattern when tracking the approximate values in the homotopy path.Table 3
Effect of measurement time.
Scenarios
Identified results (mm)
Division ofλ
Measurement time (s)
No. of iteration
x
c
h
c
1
1000.000234
1.000000
3
1
26
2
1000.000237
1.000000
3
2
30
3
1000.000187
1.000000
3
4
32Effect of measurement time duration.
(a)
(b)


## 4.1.3. Effect of Measurement Points

In this case, how the measurement points affect the accuracy and the iterative process is studied. Three, four, and six measurement points are used for crack identification to give a comparison. And the other parameters are the same as the last case. The identification results are listed in Table4. It can be found that increasing measurement points can improve the accuracy of the identification results. Figure 4 gives a comparison on the identified process of iteration. It can be noted that different number of measurement points has the similar convergence process. And the proposed method does not need a large number of measurement points.Table 4
Effect of measurement points.
Scenarios
Identified results (mm)
Measurement points (mm)
No. of iteration
x
c
h
c
1
1000.000234
1.000000
600, 1100, 1600
26
2
1000.000005
1.000000
400, 800, 1200, 1600
25
3
999.999999
1.000000
300, 600, 900, 1200, 1500, 1800
25Effect of measurement points.
(a)
(b)


## 4.1.4. Effect of Division of Homotopy Parameter

In this study, the effect of division of homotopy parameter on the identified results is discussed. Other parameters are the same as those in the last study except that the homotopy parameter is divided into 3, 4, and 5 parts. And the identified results for the study are listed in Table5. We can find that increase of division parts of the homotopy parameter, the number of iteration will increase accordingly in general but it has little effect on the accuracy of identified results. Figure 5 shows the evolution of the crack parameters in the process of iteration for different division parts of the homotopy parameter; one can see that the convergence of Newton’s method to track the homotopy path is similar.Table 5
Effect of division of homotopy parameter.
Scenarios
Identified results (mm)
Division ofλ
Measurement time (s)
No. of iteration
x
c
h
c
1
1000.000234
1.000000
3
1
26
2
1000.000003
1.000000
4
1
28
3
999.9999346
1.000000
5
1
30Effect of division of homotopy parameter.
(a)
(b)


## 4.1.5. Effect of Measurement Noise

In this section, the effect of measurement noise on the accuracy of identified results is taken into account. Again, the identification of Scenario 3 in Table2 is studied. The effect of measurement noise is simulated as a normally distributed random error with zero mean and a unit standard deviation is added to the calculated acceleration as
(26)d¨^=d¨+Ep×Noise×var(d¨),
whered¨^is the vectors of measured structural acceleration response,EPis the noise level,Noiseis standard normal distribution vector with zero mean and unit standard deviation, andvar(·)is the variance of the time history.The relative errors of crack location and depth identification are defined as(27)Relativeerror=∥αid-αtrueαtrue∥×100%,
whereαidandαtrueare the identified and true values, respectively. Table 6 gives a comparison of the identified results for 0%, 1%, 5%, and 15% noise level, respectively. This study shows that the crack parameters have been identified successfully even with 15% measurement noise. With the increase of the measurement noise level, identification errors will become larger and the max identification error is 8.93% in the location of crack and 5.24% in the crack depth.Table 6
Effect of measurement noise.
Scenarios
Noise
Identified results (mm)/rel. (%)
Division ofλ
No. of iteration
x
c
h
c
1
Nil
1000.000234
/
0.00
1.000000
/
0.00
3
26
2
1%
1021.498303
/
2.15
1.003232
/
0.32
3
22
3
5%
1070.401667
/
7.04
1.022658
/
2.27
3
23
4
15%
1089.318837
/
8.93
1.052358
/
5.24
3
19
Note: “Rel. (%)” denotes relative error.


## 4.1.6. Identification of Multiple Cracks

The same simply supported beam is studied with two cracks in this Section. The cracks are assumed to be located at 450 mm and 1000 mm from the left support of the beam, respectively. The depths of the crack are 5 mm and 1 mm, respectively. In the identification, the initial values of locations for these two cracks are set to be 500 mm and 1500 mm, respectively, and depths are both set to be 2 mm. Six measurement points as the Scenario 3 in Table4 are used in this case. And the setting of other parameters is shown in Table 7. The identification results and the relative errors for the cases with 0% and 5% noise levels are considered. It can be noted that the identified results for the two cracks are as good as those for a single crack when the measurement noise is free. Comparing the identified results with measurement noise and the free noise one, we can find that the effect of measurement noise on crack identification is noticeable, but the identified results are still satisfactory. For 5% noise level, the maximum relative errors are 1.82% and 0.28% in the crack location and depth for the first crack, respectively and 6.49% and 7.22% for the second crack, respectively. These results further illustrate the effectiveness of the proposed method.Table 7
Identification of two cracks in a simply supported beam.
Scenarios
Noise
Identified results (mm)/rel. (%) in crack 1
Identified results (mm)/rel. (%) in crack 2
Division ofλ
Measurement time (s)
No. of iteration
x
c
1
h
c
1
x
c
2
h
c
2
1
Nil
450.000000
/
0.00
5.000000
/
0.00
999.999997
/
0.00
1.000000
/
0.00
5
1
78
2
5%
441.804646
/
1.82
4.986154
/
0.28
1064.867252
/
6.49
1.072199
/
7.22
7
1
74


## 4.2. A Cantilever Beam

A uniform cross section cantilever beam with cracks is studied as another example. As Figure6 shows, the physical parameters of the beam are: elastic modulus of materialE=69.79 GPa, mass densityρ=2600 kg/m3, lengthL=1000 mm, widthw=50 mm and heighth=25 mm. when calculating the dynamic response, a sinusoidal forceg(t)=40sin(15πt)Nacts at the free end of the beam. Four measurement points located at 300 mm, 450 mm, 600 mm, and 750 mm from the clamped end are used to record the acceleration data. The two damping coefficients used for calculating Rayleigh damping matrix are both assumed to be 0.01. The time step is 0.005 s and the time duration is taken to be 2.0 s and 4.0 s, respectively. The number of the mode functions used from the CT is also set to be 10, and the way to choose the number of term is the same as the simply supported beam. Both single crack and two cracks are simulated in this case. Good identification results are listed in Table 8. For a single crack identification, the measurement of time duration may have little effect on the degree of accuracy as Section 4.1.2 shows. But for two cracks, it can improve the accuracy of the identified results in general as the unknown parameters increases in the identification. This case further illustrates the effectiveness of the proposed method.Table 8
Crack identification in a cantilever beam.
Initial
True
Identified results
2s
4s
Single crack
Crack locationxc (mm)
750
250
250.000162
250.000259
Crack depthhc (mm)
3
1
1.000001
1.000001
Division ofλ
—
—
26
20
No. of iteration
—
—
92
96
Two cracks
Crack locationsxc1/xc2 (mm)
360
/
700
250
/
600
249.603356
/
599.789017
250.009037
/
600.005717
Crack depthshc1/hc2 (mm)
2
/
7.8
1
/
6.5
0.997240
/
6.498047
1.000067
/
6.500058
Division ofλ
—
—
16
15
No. of iteration
—
—
137
137Figure 6
A cantilever beam with cracks.


## 5. Discussions

The homotopy iteration algorithm is to divide the range of homotopy parameterλ∈[0,1]into multiple subintervals; then the Newton algorithm is used to search the best convergence solution in each subinterval to obtain the updated crack parameters and another homotopy equation. Assuming that the homotopy parameter is divided intoNparts, we will haveN+1homotopy equations asα=0,1/N,…,1. When the division of homotopy parameter is only one, the homotopy iteration algorithm decreases to a simple Newton algorithm. But the Newton algorithm is with local convergence and is easy to divergent. Thus, the division of homotopy parameter is an integer greater than or equal to 2. And with the range between the initial parameters and the true parameters increasing, the finer division of homotopy parameter should be considered to improve the convergence. But in some subdivision, the Newton algorithm to track the path may be divergent to the crack depth very close to zero or the height of the beam. Because in these two conditions, the homotopy equation and its first derivative are so small that it leads to misjudgment. Thus, the new division or initial parameters should be considered. But according to authors’ calculation experience, when the subdivision increases to some value (greater than 30), the misjudgment will still occur as the above two conditions exist. That is why increasing the subdivision over some larger value does not modify the convergence properties of the algorithm.


## 6. Conclusions

In this paper, an open crack model with the CEM for Euler-Bernoulli beam system is adopted to establish the dynamics equation in the forward problem. In the inverse problem, an identification approach based on homotopy iteration algorithm is presented to identify the parameters of cracks. Numerical simulation shows that the proposed method is less affected by the initial values and is effective and accurate for crack identification when the measurement noise is free. With the increase of the measurement noise, the identification errors will become larger, but the identified results are still satisfactory. The study also shows that more measurement points can obtain more accuracy of the identification results in general but has the similar convergence process.



---

*Source: 524562-2013-11-26.xml*