Homotopy Iteration Algorithm for Crack Parameters Identification with Composite Element Method

An approach based on homotopy iteration algorithm is proposed to identify the crack parameters in beam structures. In the forward problem, a fully open crackmodel with the composite elementmethod is employed for the vibration analysis.The dynamic responses of the cracked beam in time domain are obtained from theNewmark direct integrationmethod. 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 proposedmethod. And the effects of the influencing parameters, such asmeasurement time duration,measurement points, division of the homotopy parameter and measurement noise, are studied.


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][6][7][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][10][11][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.

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 where  is Young's modulus and  0 =  3 /12 and   = ( − ℎ  ) 3 /12 are 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.

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  CEM is 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 N() and the nodal displacement vector (), where 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 () 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 where    is the estimation of the th circular frequency with -terms in the CT.
is the difference of the th circular frequency obtained with the -terms and ( − 1)terms.Tol is the tolerance value.
The stiffness matrix of the cracked beam can be obtained from the following equation: 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 K 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: 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 In this study, Rayleigh damping model [14] is adopted; that is, C =  1 M +  2 K, where  1 and  2 are two constants obtained from two different undamped natural frequencies,   and   , and their associated modal damping ratios,   and   , with the expression of and ).The vector of generalized acceleration Q, velocity Q, and displacement Q of the structure can be obtained from ( 9) by direct integration.Further, the vector of physical acceleration d (), velocity ḋ (), and displacement d() can be obtained from The generalized force vector is G() = S(  )  (), where () is an external force acting at the location of   from the left support.

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 d * = ( d 1 , d 2 , . . ., d  ).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 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 where  d (, )/  and  d (, )/ℎ  are the acceleration response sensitivities with respect to the crack location and depth, respectively.

Homotopy Iteration Algorithm for Crack Identification.
The basic idea of homotopy iteration algorithm is to introduce a parameter  into the nonlinear equation () = 0 to construct a family of homotopy mapping (, ), where  ∈ [0, 1] is named the homotopy parameter.The homotopy equation is assumed to be (, ) = () + (1 − ) 0 () in this paper.When  = 1,  is the mapping of () as (12) shows, and, when  = 0,  is assumed to be the mapping of  0 ().And  0 () is set to be similar to () but, related to the initial point, can be expressed as  0 () = ( d (, ) − d ( 0 , )) d (, )/, where  0 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: Equations ( 14) and (15) show that as the homotopy parameter  increases from zero to one, the homotopy equation (, ) varies from  0 () 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 2. Trace the homotopy path from the initial value ( 0 ,  0 ) and give an increment  to the homotopy parameter  as  1 =  0 + .
Step 3. Using Newton iterative method to calculate an updated vector of crack parameter  1 = [ 1 , ℎ 1 ] in the homotopy path, we have where (,  1 )/ is the homotopy equation with respect to the crack parameters.The value of ( 0 ,  1 ) and (,  1 )/ will be introduced in next section.
Step 4. Check whether the updated crack parameter vector ] 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: where Tol1 is the first tolerance value for convergence and is taken as 1.0 × 10 −3 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 ( 1 ,  1 ) meets the following criterion: where Tol2 is 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  1 would then be the required set of identification results.

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 where  Q/  ,  Q/ℎ  ,  Q/  ,  Q/ℎ  , Q/  , and Q/ℎ  are the vectors of generalized acceleration, velocity, and displacement sensitivities with respect to the crack location and depth, respectively.K/  and K/ℎ  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 K can be obtained from (7), the derivatives K/  and K/ℎ  can be obtained directly.As the dynamic responses Q and Q 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 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 where  2 K/ 2  ,  2 K/  ℎ  ,  2 K/ℎ 2  , and  2 K/ℎ    are the second partial derivative of the stiffness matrix with respect to the crack parameters, and  2 K/  ℎ  =  2 K/ ℎ    .Since the response sensitivities Q/ℎ  ,  Q/  , and  Q/ℎ  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 Furthermore, the first derivative of the homotopy mapping in (13) with respect to the unknown vector (, )/ can be obtained as where The first derivatives  0 ()/  and  0 ()/ℎ  can also be obtained in a similar way.(25)

Numerical Simulations
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 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.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.

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.

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.

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.

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 where ̂d is the vectors of measured structural acceleration response,   is the noise level,  oise is standard normal distribution vector with zero mean and unit standard deviation, and var(⋅) is the variance of the time history.
where  id and  true are 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.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.

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.79GPa, mass density  = 2600 kg/m 3 , length  = 1000 mm, width  = 50 mm and height ℎ = 25 mm.when calculating the dynamic response, a sinusoidal force () = 40 sin(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.2shows.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.

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  + 1 homotopy 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.

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.

Figure 1 :
Figure 1: A simply supported beam with cracks.

Figure 5 :
Figure 5: Effect of division of homotopy parameter.

Table 1 :
Comparison of the convergence property with different number of terms in CT.

Table 2 :
Effect of the initial value on crack identification.

Table 3 :
Effect of measurement time.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

Table 4 :
Effect of measurement points.

Table 5 :
Effect of division of homotopy parameter.

Table 6 :
Effect of measurement noise.

Table 7 :
Identification of two cracks in a simply supported beam.

Table 8 :
Crack identification in a cantilever beam.