Computational Methods for FractureView this Special Issue
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 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.
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  proposed a composite element method (CEM) combining the conventional FEM and classical theory (CT) to solve the structural dynamical problems. And Lu and Law  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 , 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  used the time histories of vibration response of the structure to identify damage in smart structures. Lu and Liu  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  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.  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  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  proposed a homotopy analysis method for solving highly nonlinear problems in science, finance, and engineering areas. Alexander and Yorke  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.  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. As shown in Figure 2, the flexural rigidityadjacent to theth crack can be written as whereis Young’s modulus andandare the moments of inertia of the intact and cracked cross sections, respectively.andare the width and height of the beam, respectively.andare 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.
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 fieldis 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 matrixand the nodal displacement vector, where whereis the length of the beam.
The second part can be obtained by the linear combination of the multiplication of analytical shape functionwith a vector of coefficient whereis a special value which can be obtained according to the boundary conditions of the beam andis the number of the mode functions used from the CT. Such that, for a simply supported beam,can be selected as Lu and Law  proposed. And the coefficientdenotes the contribution of thein the total displacement field. For a uniform beam, less terms of the shapes function in CT also can be obtained with accuracy results . And the number of termscan be determined by a frequency convergence test as the previous study by Lu and Law . The frequency convergence criterion is defined as whereis the estimation of theth circular frequency with-terms in the CT.is the difference of theth circular frequency obtained with the-terms and-terms.is the tolerance value.
The displacement field of the CEM for a uniform Euler-Bernoulli beam element can be written as whereis the generalized shape function of CEM.
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. Foradjacent to theth crack can be obtained from (1); the global stiffness matrixcan 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: whereis the mass densityis 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  is adopted; that is,, whereandare two constants obtained from two different undamped natural frequencies,and, and their associated modal damping ratios,and, with the expression ofand. The vector of generalized acceleration, velocity, and displacementof the structure can be obtained from (9) by direct integration. Further, the vector of physical acceleration, velocity, and displacementcan be obtained from
The generalized force vector is, whereis an external force acting at the location offrom 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 locationand crack depth; that is,. If there aremeasured 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 wheredenotes the the measured time instants,is the number of time step, andis 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 whereandare 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 parameterinto the nonlinear equationto construct a family of homotopy mapping, whereis named the homotopy parameter. The homotopy equation is assumed to bein this paper. When,is the mapping ofas (12) shows, and, when,is assumed to be the mapping of. Andis set to be similar tobut, related to the initial point, can be expressed as, whereis the arbitrarily selected initial value vector. In this way, a family of homotopy mappingcan 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 parameterincreases from zero to one, the homotopy equationvaries fromto. Hence, the parameters of the crack can be obtained by following the variations of parameter, and, when, the identification vectoris 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, divide the range of homotopy parameterintoequal parts, and then obtain, let.
Step 2. Trace the homotopy path from the initial valueand give an incrementto the homotopy parameteras.
Step 3. Using Newton iterative method to calculate an updated vector of crack parameterin the homotopy path, we have whereis the homotopy equation with respect to the crack parameters. The value ofandwill be introduced in next section.
Step 4. Check whether the updated crack parameter vectoris 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 obtainedin homotopy path satisfies the following convergence conditions: whereis the first tolerance value for convergence and is taken asin 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 whenmeets the following criterion: whereis the second tolerance value and is taken to be 0.1 for all the study cases.
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 where,,,,, andare the vectors of generalized acceleration, velocity, and displacement sensitivities with respect to the crack location and depth, respectively.andare the first partial derivatives of the stiffness matrix with respect to the crack parameters. Since the flexural rigidityinvolves the crack location and depth and the global stiffness matrixcan be obtained from (7), the derivativesandcan be obtained directly. As the dynamic responsesandhave 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 asFurthermore, 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,,, andare the second partial derivative of the stiffness matrix with respect to the crack parameters, and. Since the response sensitivities,, andhave been calculated from (18), the second partial derivative of generalized acceleration response with respect to the unknown variablescan 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 vectorcan be obtained as where The first derivativesandcan 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 GPa, mass density kg/m3, length mm, width mm, and height mm. An impulsive force acts at the 700 mm from the left support in the negative-direction from s to s with 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 firstfrequency 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.
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 parametershould be considered as to track the correct path. Scenarios 3 and 4 show that the method is less affected by the initial values.
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.
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.
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.
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 whereis the vectors of measured structural acceleration response,is the noise level,is standard normal distribution vector with zero mean and unit standard deviation, andis the variance of the time history.
The relative errors of crack location and depth identification are defined as whereandare 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.
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 GPa, mass density kg/m3, length mm, width mm and height mm. when calculating the dynamic response, a sinusoidal forceacts 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.
The homotopy iteration algorithm is to divide the range of homotopy parameterinto 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 intoparts, we will havehomotopy equations as. 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.
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.
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.
M. I. Friswell and J. E. Mottershead, Finite Element Model Updating in Structural Dynamics, Kluwer Academic, Dodrecht, The Netherlands, 1995.
A. G. Peano, “Hierarchies of conforming finite elements for plane elasticity and plate bending,” Computers and Mathematics with Applications, vol. 2, no. 3-4, pp. 211–224, 1976.View at: Google Scholar
P. Zeng, “Composite element method for vibration analysis of structure,” Journal of Sound and Vibration, vol. 218, no. 4, pp. 619–696, 1998.View at: Google Scholar
J. Cattarius and D. J. Inman, “Time domain analysis for damage detection in smart structures,” Mechanical Systems and Signal Processing, vol. 11, no. 3, pp. 409–423, 1997.View at: Google Scholar
X. B. Lu, J. K. Liu, and Z. R. Lu, “A two-step approach for crack identification in beam,” Journal of Sound and Vibration, vol. 332, no. 2, pp. 282–293, 2013.View at: Google Scholar
B. C. Eaves, F. J. Gould, H. O. Peitgen, and M. J. Todd, Homotopy Methods and Global Convergence, Plenum Press, New York, NY, USA, 1983.
K. J. Bathe, Finite Element Procedures in Engineering Analysis, Prentice Hall, Upper Saddle River, NJ, USA, 1982.