Homotopy Iteration Algorithm for Crack Parameters Identification with Composite Element Method

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.

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, theℎ-version and the𝑐-version. As the former by increasing the number of element mesh is just like that of𝑝-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 the𝑐-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

Figure 1 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 is𝑙_𝑐=1.5𝑑. As shown in Figure 2, the flexural rigidity𝐸𝐼(𝑥)adjacent to the𝑖th crack can be written as 𝐸𝐼(𝑥)=𝐸𝐼₀−𝐸𝐼₀−𝐼_𝑐𝑖𝑥−𝑥_𝑖1𝑥_𝑖−𝑥_𝑖1,𝑥_𝑖1≤𝑥≤𝑥_𝑐𝑖,𝐸𝐼₀−𝐸𝐼₀−𝐼_𝑐𝑖𝑥_𝑖2−𝑥𝑥_𝑖2−𝑥_𝑖,𝑥_𝑐𝑖≤𝑥≤𝑥_𝑖2, (1) where𝐸is Young’s modulus and𝐼₀=𝑤𝑑³/12and𝐼_𝑐𝑖=𝑤(𝑑−ℎ_𝑐𝑖)³/12are the moments of inertia of the intact and cracked cross sections, respectively.𝑤and𝑑are the width and height of the beam, respectively.𝑥_𝑖1=𝑥_𝑐𝑖−𝑙_𝑐and𝑥_𝑖2=𝑥_𝑐𝑖+𝑙_𝑐are 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 rigidity𝐸𝐼(𝑥)adjacent to a crack.

2.2. Dynamic Response of the Structure with CEM

As shown in Figure 1, 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 the𝑐-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 field𝑑_CEMis 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 matrix𝐍(𝑥)and the nodal displacement vector𝑞(𝑡), 𝑑_FEM(𝑥,𝑡)=𝐍(𝑥)𝑞(𝑡), (2) where 𝐍(𝑥)=𝑥1−3²𝐿²𝑥+2³𝐿³𝑥,𝑥−2²𝐿+𝑥³𝐿²,3𝑥²𝐿²𝑥−2³𝐿³,𝑥³𝐿²−𝑥²𝐿, (3) where𝐿is the length of the beam.

The second part can be obtained by the linear combination of the multiplication of analytical shape function𝜙(𝑥)with a vector of coefficient𝑐(𝑡)𝑑_CT(𝑥,𝑡)=^𝑛_𝑟=1𝜙_𝑟(𝑥)𝑐_𝑟(𝑡),(𝑟=1,2,…,𝑛), (4) where𝜙_𝑟(𝑥)is a special value which can be obtained according to the boundary conditions of the beam and𝑛is the number of the mode functions used from the CT. Such that, for a simply supported beam,𝜙_𝑟(𝑥)=sin(𝑟𝜋𝑥/𝐿)can be selected as Lu and Law [4] proposed. And the coefficient𝑐_𝑟(𝑡)denotes the contribution of the𝜙_𝑟(𝑥)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 terms𝑛can be determined by a frequency convergence test as the previous study by Lu and Law [7]. The frequency convergence criterion is defined as max_{𝑟=1,2,…,10}Δ𝜔^𝑛_𝑟𝜔^𝑛_𝑟<Tol, (5) where𝜔^𝑛_𝑟is the estimation of the𝑟th circular frequency with𝑛-terms in the CT.Δ𝜔^𝑛_𝑟=𝜔^𝑛_𝑟−𝜔_𝑟^𝑛−1is the difference of the𝑟th circular frequency obtained with the𝑛-terms and(𝑛−1)-terms.Tolis the tolerance value.

The displacement field of the CEM for a uniform Euler-Bernoulli beam element can be written as 𝑑_CEM(𝑥,𝑡)=𝑑_FEM(𝑥,𝑡)+𝑑_CT(𝑥,𝑡)=𝐒(𝑥)𝐐(𝑡), (6) where𝐒(𝑥)=[𝐍(𝑥),𝜙₁(𝑥),𝜙₂(𝑥),…,𝜙_𝑛(𝑥)]is the generalized shape function of CEM.

The stiffness matrix of the cracked beam can be obtained from the following equation: 𝐊=^𝐿₀𝑑²𝐒(𝑥)^𝑇𝑑𝑥²(𝑥)𝑑𝐸𝐼²𝐒(𝑥)𝑑𝑥²𝑑𝑥. (7)

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. For𝐸𝐼(𝑥)adjacent to the𝑖th crack can be obtained from (1); the global stiffness matrix𝐊can 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) where𝜌is the mass density𝐴is 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)

In this study, Rayleigh damping model [14] is adopted; that is,𝐂=𝑎₁𝐌+𝑎₂𝐊, where𝑎₁and𝑎₂are two constants obtained from two different undamped natural frequencies,𝜔_𝑖and𝜔_𝑗, and their associated modal damping ratios,𝜉_𝑖and𝜉_𝑗, with the expression of𝑎₁=2𝜔_𝑗𝜔_𝑖(𝜔_𝑗𝜉_𝑖−𝜔_𝑖𝜉_𝑗)/(𝜔_𝑗²−𝜔_𝑖²)and𝑎₂=2(𝜔_𝑗𝜉_𝑗−𝜔_𝑖𝜉_𝑖)/(𝜔_𝑗²−𝜔_𝑖²). The vector of generalized acceleration𝐐̈, velocity𝐐̇, and displacement𝐐of the structure can be obtained from (9) by direct integration. Further, the vector of physical acceleration𝐝̈(𝑥), velocity𝐝̇(𝑥), and displacement𝐝(𝑥)can be obtained from𝐝̈(𝑥)(𝑥)̈𝐝̇(𝑥)(𝑥)̇𝐝(𝑥)(𝑥)=𝐒𝐐,=𝐒𝐐,=𝐒𝐐. (10a)(10b)(10c)

The generalized force vector is𝐆(𝑡)=𝐒(𝑥_𝑙)^𝑇𝑔(𝑡), where𝑔(𝑡)is an external force acting at the location of𝑥_𝑙from 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 location𝑥_𝑐and crack depthℎ_𝑐; that is,𝜶=[𝑥_𝑐,ℎ_𝑐]. If there are𝑚measured points in the beam structure, the measured acceleration response can be expressed as𝐝̈^∗=(𝐝̈₁,𝐝̈₂,…,𝐝̈_𝑚). 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 1𝐽=2^𝑛𝑡𝑚_𝑖=1𝑗=1𝐝̈(𝜶,𝑡)−𝐝̈^∗(𝜶,𝑡)²,𝑖=1,2,…,𝑛𝑡,𝑗=1,2,…,𝑚, (11) where𝑡denotes the the measured time instants,𝑛𝑡=𝑡/𝑑𝑡is the number of time step, and𝑑𝑡is 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 𝑓(𝜶)=^𝑛𝑡𝑚∑_𝑖=1∑_𝑗=1𝐝̈(𝜶,𝑡)−𝐝̈^∗()𝜕𝐝̈()𝜶,𝑡𝜶,𝑡𝜕𝑥_𝑐=0,_𝑛𝑡^𝑚∑_𝑖=1∑_𝑗=1𝐝̈(𝜶,𝑡)−𝐝̈^∗()𝜕𝐝̈()𝜶,𝑡𝜶,𝑡𝜕ℎ_𝑐=0, (12) where𝜕𝐝̈(𝜶,𝑡)/𝜕𝑥_𝑐and𝜕𝐝̈(𝜶,𝑡)/𝜕ℎ_𝑐are 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 equation𝑓(𝜶)=0to construct a family of homotopy mapping𝐻(𝜶,𝜆), where𝜆∈[0,1]is named the homotopy parameter. The homotopy equation is assumed to be𝐻(𝜶,𝜆)=𝜆𝑓(𝜶)+(1−𝜆)𝑓₀(𝜶)in this paper. When𝜆=1,𝐻is the mapping of𝑓(𝜶)as (12) shows, and, when𝜆=0,𝐻is assumed to be the mapping of𝑓₀(𝜶). And𝑓₀(𝜶)is set to be similar to𝑓(𝜶)but, related to the initial point, can be expressed as𝑓₀(𝜶)=(𝐝̈(𝜶,𝑡)−𝐝̈(𝜶₀,𝑡))𝜕𝐝̈(𝜶,𝑡)/𝜕𝜶, where𝜶₀is the arbitrarily selected initial value vector. In this way, a family of homotopy mapping𝐻can be constructed instead of a single mapping𝑓. Thus, the homotopy equation can meet the following equation: 𝐻(𝜶,𝜆)=𝜆𝑓(𝜶)+(1−𝜆)𝑓₀(𝜶)𝐻()=0,𝜶,0=𝑓₀(𝜶)𝐻()(𝜶)=0,𝜶,1=𝑓=0. (13)(14)(15) Equations (14) and (15) show that as the homotopy parameter𝜆increases from zero to one, the homotopy equation𝐻(𝜶,𝜆)varies from𝑓₀(𝜶)to𝑓(𝜶). Hence, the parameters of the crack can be obtained by following the variations of parameter𝜆, and, when𝐻(𝜶,𝜆)=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,ℎ_𝑐0], divide the range of homotopy parameter𝜆∈[0,1]into𝑛equal parts, and then obtain𝑑𝜆=1/𝑛, let𝜆₀=0.

Step 2. Trace the homotopy path from the initial value(𝜶₀,𝜆₀)and give an increment𝑑𝜆to the homotopy parameter𝜆as𝜆₁=𝜆₀+𝑑𝜆.

Step 3. Using Newton iterative method to calculate an updated vector of crack parameter𝜶₁=[𝑥_𝑐1,ℎ_𝑐1]in the homotopy path, we have 𝜶₁=𝜶₀−𝜕𝐻(𝜶,𝜆₁)𝜕𝜶_𝜶=𝜶0⁻¹𝐻𝜶₀,𝜆₁, (16) where𝜕𝐻(𝜶,𝜆₁)/𝜕𝜶is the homotopy equation with respect to the crack parameters. The value of𝐻(𝜶₀,𝜆₁)and𝜕𝐻(𝜶,𝜆₁)/𝜕𝜶will be introduced in next section.

Step 4. Check whether the updated crack parameter vector𝜶₁=[𝑥_𝑐1,ℎ_𝑐1]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𝜶₀=𝜶₁; repeat from Steps 3 to 4, until the obtained(𝜶₁,𝜆₁)in homotopy path satisfies the following convergence conditions: 𝜶₁−𝜶₀𝜶₁≤Tol1, (17a) whereTol1is the first tolerance value for convergence and is taken as1.0×10⁻³in 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 when𝐻(𝜶₁,𝜆₁)meets the following criterion: 𝐻𝜶₁,𝜆₁≤Tol2, (17b) whereTol2is the second tolerance value and is taken to be 0.1 for all the study cases.

Step 7. The pair of unknowns(𝜶₁,𝜆₁)is taken as the new set of initial value(𝜶₀,𝜆₀)and repeat Steps 2 to 5 until the homotopy parameter𝜆reaches 1.0. The final𝜶₁would 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) where𝜕𝐐̈/𝜕𝑥_𝑐,𝜕𝐐̈/𝜕ℎ_𝑐,𝜕𝐐̇/𝜕𝑥_𝑐,𝜕𝐐̇/𝜕ℎ_𝑐,𝜕𝐐/𝜕𝑥_𝑐, and𝜕𝐐/𝜕ℎ_𝑐are the vectors of generalized acceleration, velocity, and displacement sensitivities with respect to the crack location and depth, respectively.𝜕𝐊/𝜕𝑥_𝑐and𝜕𝐊/𝜕ℎ_𝑐are the first partial derivatives of the stiffness matrix with respect to the crack parameters. Since the flexural rigidity𝐸𝐼(𝑥)involves the crack location and depth and the global stiffness matrix𝐊can be obtained from (7), the derivatives𝜕𝐊/𝜕𝑥_𝑐and𝜕𝐊/𝜕ℎ_𝑐can be obtained directly. As the dynamic responses𝐐and𝐐̇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)(19b)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 𝐌𝜕²𝐐̈𝜕𝑥²_𝑐𝜕+𝐂²𝐐̇𝜕𝑥²_𝑐𝜕+𝐊²𝐐𝜕𝑥²_𝑐=−𝑎₂𝜕²𝐊𝜕𝑥²_𝑐𝐐̇−𝜕²𝐊𝜕𝑥²_𝑐𝐐−2𝑎₂𝜕𝐊𝜕𝑥_𝑐𝜕𝐐̇𝜕𝑥_𝑐+𝜕𝐊𝜕𝑥_𝑐𝜕𝐐𝜕𝑥_𝑐,𝐌𝜕²𝐐̈𝜕𝑥_𝑐𝜕ℎ_𝑐𝜕+𝐂²𝐐̇𝜕𝑥_𝑐𝜕ℎ_𝑐𝜕+𝐊²𝐐𝜕𝑥_𝑐𝜕ℎ_𝑐=−𝑎₂𝜕²𝐊𝜕𝑥_𝑐𝜕ℎ_𝑐𝐐̇−𝜕²𝐊𝜕𝑥_𝑐𝜕ℎ_𝑐𝐐−𝑎₂𝜕𝐊𝜕ℎ_𝑐𝜕𝐐̇𝜕𝑥_𝑐−𝜕𝐊𝜕ℎ_𝑐𝜕𝐐𝜕𝑥_𝑐−𝑎₂𝜕𝐊𝜕𝑥_𝑐𝜕𝐐̇𝜕ℎ_𝑐−𝜕𝐊𝜕𝑥_𝑐𝜕𝐐𝜕ℎ_𝑐,𝐌𝜕²𝐐̈𝜕ℎ²_𝑐𝜕+𝐂²𝐐̇𝜕ℎ²_𝑐𝜕+𝐊²𝐐𝜕ℎ²_𝑐=−𝑎₂𝜕²𝐊𝜕ℎ²_𝑐𝐐̇−𝜕²𝐊𝜕ℎ²_𝑐𝐐−2𝑎₂𝜕𝐊𝜕ℎ_𝑐𝜕𝐐̇𝜕ℎ_𝑐+𝜕𝐊𝜕ℎ_𝑐𝜕𝐐𝜕ℎ_𝑐,𝐌𝜕²𝐐̈𝜕ℎ_𝑐𝜕𝑥_𝑐𝜕+𝐂²𝐐̇𝜕ℎ_𝑐𝜕𝑥_𝑐𝜕+𝐊²𝐐𝜕ℎ_𝑐𝜕𝑥_𝑐=−𝑎₂𝜕²𝐊𝜕ℎ_𝑐𝜕𝑥_𝑐𝐐̇−𝜕²𝐊𝜕ℎ_𝑐𝜕𝑥_𝑐𝐐−𝑎₂𝜕𝐊𝜕𝑥_𝑐𝜕𝐐̇𝜕ℎ_𝑐−𝜕𝐊𝜕𝑥_𝑐𝜕𝐐𝜕ℎ_𝑐−𝑎₂𝜕𝐊𝜕ℎ_𝑐𝜕𝐐̇𝜕𝑥_𝑐−𝜕𝐊𝜕ℎ_𝑐𝜕𝐐𝜕𝑥_𝑐, (20) where𝜕²𝐊/𝜕𝑥²_𝑐,𝜕²𝐊/𝜕𝑥_𝑐𝜕ℎ_𝑐,𝜕²𝐊/𝜕ℎ²_𝑐, and𝜕²𝐊/𝜕ℎ_𝑐𝜕𝑥_𝑐are the second partial derivative of the stiffness matrix with respect to the crack parameters, and𝜕²𝐊/𝜕𝑥_𝑐𝜕ℎ_𝑐=𝜕²𝐊/𝜕ℎ_𝑐𝜕𝑥_𝑐. Since the response sensitivities𝜕𝐐/𝜕ℎ_𝑐,𝜕𝐐̇/𝜕𝑥_𝑐, and𝜕𝐐̇/𝜕ℎ_𝑐have been calculated from (18), the second partial derivative of generalized acceleration response with respect to the unknown variables(𝑥_𝑐,ℎ_𝑐)can be calculated by the Newmark integration with (20). Then, the physical acceleration response sensitivity can be obtained as𝜕²𝐝̈𝜕𝑥²_𝑐(𝑥)𝜕=𝐒²𝐐̈𝜕𝑥²_𝑐,𝜕²𝐝̈𝜕ℎ²_𝑐(𝑥)𝜕=𝐒²𝐐̈𝜕ℎ²_𝑐,𝜕²𝐝̈𝜕𝑥_𝑐𝜕ℎ_𝑐=𝜕²𝐝̈𝜕ℎ_𝑐𝜕𝑥_𝑐(𝑥)𝜕=𝐒²𝐐̈𝜕𝑥_𝑐𝜕ℎ_𝑐. (21a)(21b)(21c)

Furthermore, the first derivative of the homotopy mapping in (13) with respect to the unknown vector𝜕𝐻(𝜶,𝜆)/𝜕𝜶can be obtained as 𝜕𝐻(𝜶,𝜆)=𝜆(𝜶)𝜕𝜶𝜕𝑓𝜕𝑥_𝑐+(1−𝜆)𝜕𝑓₀(𝜶)𝜕𝑥_𝑐𝜆𝜕𝑓(𝜶)𝜕ℎ_𝑐+(1−𝜆)𝜕𝑓₀(𝜶)𝜕ℎ_𝑐, (22) where 𝜕𝑓(𝜶)𝜕𝑥_𝑐=^𝑛𝑡𝑚_𝑖=1𝑗=1𝐝̈−𝐝̈^∗𝜕²𝐝̈𝜕𝑥²_𝑐+𝜕𝐝̈𝜕𝑥_𝑐²_{𝑛𝑡𝑚𝑖=1𝑗=1},𝐝̈−𝐝̈^∗𝜕²𝐝̈𝜕ℎ_𝑐𝜕𝑥_𝑐+𝜕𝐝̈𝜕𝑥_𝑐𝜕𝐝̈𝜕ℎ_𝑐,𝜕𝑓(𝜶)𝜕ℎ_𝑐=^𝑛𝑡𝑚_𝑖=1𝑗=1𝐝̈−𝐝̈^∗𝜕²𝐝̈𝜕𝑥_𝑐𝜕ℎ_𝑐+𝜕𝐝̈𝜕ℎ_𝑐𝜕𝐝̈𝜕𝑥_{𝑐𝑛𝑡𝑚𝑖=1𝑗=1},𝐝̈−𝐝̈^∗𝜕²𝐝̈𝜕ℎ²_𝑐+𝜕𝐝̈𝜕ℎ_𝑐². (23)(24) The first derivatives𝜕𝑓₀(𝜶)/𝜕𝑥_𝑐and𝜕𝑓₀(𝜶)/𝜕ℎ_𝑐can 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 Figure 1. The parameters of the beam are: Young’s modulus𝐸=69.79 GPa, mass density𝜌=2600 kg/m³, length𝐿=2000 mm, width𝑤=50 mm, and heightℎ=25 mm. An impulsive force acts at the 700 mm from the left support in the negative𝑦-direction from𝑡=0 s to𝑡=0.04 s with 𝑔(𝑡)=(),()().10000𝑡𝑁,0s≤𝑡≤0.02s100000.04−𝑡𝑁,0.02s<𝑡≤0.04s (25) 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 Table 1, one can find that when the number of term in CT is𝑛, only the first𝑛frequency 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.

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 Table 2. 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.

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 Table 2 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.

(a)

(b)

Figure 3 

Effect of measurement time duration.

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 Table 4. 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.

(a)

(b)

Figure 4 

Effect of measurement points.

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 Table 5. 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.

(a)

(b)

Figure 5 

Effect of division of homotopy parameter.

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 Table 2 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 𝐝̈=𝐝̈+𝐸_𝑝×𝑁_oise×var𝐝̈, (26) where𝐝̈is the vectors of measured structural acceleration response,𝐸_𝑃is the noise level,𝑁_oiseis 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 𝜶Relativeerror=^id−𝜶_true𝜶_true×100%, (27) 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.

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 Table 4 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.

4.2. A Cantilever Beam

A uniform cross section cantilever beam with cracks is studied as another example. As Figure 6 shows, the physical parameters of the beam are: elastic modulus of material𝐸=69.79 GPa, mass density𝜌=2600 kg/m³, length𝐿=1000 mm, width𝑤=50 mm and heightℎ=25 mm. when calculating the dynamic response, a sinusoidal force𝑔(𝑡)=40sin(15𝜋𝑡)𝑁acts 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.

Figure 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 into𝑁parts, we will have𝑁+1homotopy equations as𝜶=0,1/𝑁,…,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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported in part by the NSFC (11172333 and 11272361), the Fundamental Research Funds for the Central Universities (13lgzd06), and the General Financial Grant from the China Postdoctoral Science Foundation (2013M531893). Such financial aids are gratefully acknowledged.