3D crack identification using the Nelder-Mead Simplex algorithm combined with a random generation of crack positions Method;

A BSTRACT . In this paper, we present a scheme for cracks identification in three-dimensional linear elastic mechanical components. The scheme uses a boundary element method for solving the forward problem and the Nelder-Mead simplex numerical optimization algorithm coupled with a low discrepancy sequence in order to identify an embedded crack. The crack detection process is achieved through minimizing an objective function defined as the difference between measured strains and computed ones, at some specific sensors on the domain boundaries. Through the optimization procedure, the crack surface is modelled by geometrical parameters, which serve as identity variables. Numerical simulations are conducted to determine the identity parameters of an embedded elliptical crack, with measures randomly perturbed and the residual norm regularized in order to provide an efficient and numerically stable solution to measurement noise. The accuracy of this method is investigated in the identification of cracks over two examples. Through the treated examples, we showed that the method exhibits good stability with respect to measurement noise and convergent results could be achieved without restrictions on the selected initial values of the crack parameters.


INTRODUCTION
he identification of cracks in mechanical components remains as one of the challenges of inverse problems in mechanics. To identify an internal crack in a structure (aeronautic, nuclear, transport…), in order to make a decision about its structural integrity, is an issue that faces engineers. The detection of such flaws, in order to evaluate the damage level, is an important step in the process of health monitoring of those components. In the experimental methods, engineers process the response of the component to an applied field (static or dynamic), and give an approximate evaluation of the crack's identity (size and position). These experimental methods have their specific advantages and limitations depending on the application and the setup of the measurement system employed [1][2][3]. If this response is used in a numerical technique, a closer solution to the crack's identity could be made. Processing experimental data to complete the information needed to analyze a mechanical system leads to solve an inverse problem [4][5][6]. For example, a cracked component subjected to known loads and displacements on the boundaries and unknown crack's size and position. The given boundary conditions and geometry are insufficient to analyze the system and an iterative procedure is necessary, starting by a random initial crack geometry. Several previous works have already dealt with problems of crack identification by 2D or 3D inverse analysis using various numerical technics [7][8][9][10][11][12]. In these earlier works, the component undergoes a quasi-static load or potential and the functional to minimize is the amplitude of the differences between theoretical and computed values of responses at sensors on the boundaries. To solve the direct problem, several numerical techniques could be used to compute the responses at fixed nodes. However, for iterative methods, the computation time of the direct problem is a determining parameter in the choice of the method and in this context, the extended finite element method (XFEM) [13,14]and the Dual Boundary Element Method (DBEM) [15,16] are the most used ones, since they do not need remeshing and domain subdivision. Amoura et al [10] developed a crack identification algorithm for 2-D and axisymmetric structures using a coupled DBEM and Nelder-Mead function minimization method (The simplex method) [17] to set a stable procedure. In their work, they used a lowdiscrepancy sequence (LDS) to produce the initial crack's identity for the simplex launch, which considerably reduced the computing time. In this work, an extension of the procedure to 3D problems is done with application to the identification of elliptical embedded cracks.

CRACK IDENTIFICATION PROCEDURE
or the solution of the direct problem, we consider a 3D convex domain with an elliptic internal discontinuity as a crack model. Modelling the crack geometry in the form of an ellipse makes it possible to approximate roughly a variety of real shapes with a reduced number of parameters. The crack consists of two identical surfaces of opposite normals. Moreover, we assume that the strains are measurable at certain points (sensors) on the skin of the domain. The first step is to perform a direct calculation with an initial configuration of the crack's identity (geometry and position), considered as the real configuration. The direct simulation by the Dual Boundary Element Method (DBEM) allows evaluation of the strain field at fixed points. The latter is used in an iterative optimization procedure in order to find the real position and size of the crack (actual identity). The optimization procedure uses the Nelder-Mead Simplex Algorithm (NMSA). This algorithm requires a starting identity, which is close to the actual one, in order to avoid convergence to local minima. Therefore, and to give the simplex its optimal starting solution, we will use a Low Discrepancy Sequence (LDS) to generate a sample of identities (solutions) distributed over the whole domain. The main advantage of quasi-random sequence in comparison to real-random sequence is it distributes uniformly hence there is no larger gaps and no cluster formation, this leads to spread the sample of identities regularly over the entire domain. The process of finding the optimal solution consists in minimizing the following objective function (OF): where: ( ) e x : The objective function (or error function) produced by the crack's identity vector x  ex : Strain's field at sensor points for experimental data T F  : Strain's field at sensor points for guessed position of the crack 0 e : Value used to normalize the functional. The functional (1) to be minimized is expressed as the norm of the gaps vector in the two sets of readings at sensors: guessed crack shape and experimental crack shape. The crack surface is parameterized by an identity vector, which groups shape and position parameters. For an elliptical crack ( Fig. 1), identity parameters may be: the coordinates of its center point C (x c , y c , z c ) the lengths of the two radii a and b, and the angles of rotation of its normal about the three axes θx, θy and θz (Euler angles). In the present work, restrictions are made on the values of the parameters, assuring that the crack remains inside the domain of interest (Embedded crack). In order to stabilize the problem and select a useful and stable solution, the inverse problem requires regularization [18,19]. The latter introduces additional information about the desired solution. Though many types of supplementary information on the solution x are possible, the prevalent technique to regularization of ill-posed problem is to require that the norm of the solution be bounded. An initial approximation x* of the solution is included in the side constraint, which is defined as: We can then control the smoothness of the regularized solution by means of the side constraint Ω(x). When the side constraint Eqn. 2 is added, we must waive the requirement that    * in the least square problem Eqn. 1 and instead search for a solution that offers a fair balance between minimizing Ω(x) and minimizing the residual norm e(x). It means that the regularized solution, with small norm and residual error norm is not too far, from the needed unknown solution to the perturbed problem underlying the real problem. Tikhonov's regularization is one of the most used technique and the one, which well adapts with our problem [20]. This technique defines the regularized solution x α as the minimizer of the following weighted combination of the residual norm and the side constraint: The parameter α is added to balance the perturbation error and the regularization error in the regularized solution. The discrepancy principle [21,22], due to Morozov is the method adopted to choose α, which equals to fixing the regularization parameter such that the residual norm for the regularized solution satisfies: In order to validate the stability of the solution with respect to measurement noise, the values of the strains measured at the sensor points are disturbed with the standard deviation of noise defined by a Gaussian law.

DUAL BOUNDARY ELEMENT METHOD FOR 3D ELASTOSTATIC PROBLEMS
he dual boundary element method is based on dual equations, which are the displacement and traction boundary integral equations [23,24]. Considering a linear elastic solid Ω with boundary Γ (Fig. 2), the Somigliana's identity relating the displacement ui at an internal point P to the displacement uj and traction tj on the surface is given by Where i, j= 1, 2, 3;    Differentiating the displacement Eqn. (5) with respect to P and using the Hooke's law yields the Somigliana's identity for stresses at an interior point P T P Q respectively. Considering the limiting process as an internal point P goes to the boundary, the displacement boundary equation can be written as where  denotes Cauchy principal value integral and     ij c P is a jump term arising from taking the first integral of the boundary, dependent on the geometry of the body. For a smooth boundary    1 2 ij c P The stress boundary integral equation is obtained by taking the limiting form of the interior stress Eqn. 6, as an internal point P goes to the boundary. It can be written as The dual boundary element method eliminates need of domain sub-division and reduces both the computational effort and the perturbations due to mesh adjustments, since only domain boundaries are meshed. The finite elements used for meshing are of two types: regular elements and singular ones. The regular elements are isoparametric triangular elements with six nodes employed to mesh the domain boundaries. The singular elements, used to mesh the two crack surfaces, have their nodes shifted towards the center of the element in order to satisfy, smoothness at the boundary nodes, continuity of the displacement derivatives and boundary curvatures at these points; these conditions are required for the evaluation of the improper integrals inherent in a dual boundary equation formulation [24]. The dual equations Eqn. 7 and 9 are applied to the elements of each side of the crack to ensure a non-singular system to be solved. After evaluation of integrals over each element, a discretized system of equations is obtained and solved to get the unknown tensions and displacements. Afterwards, the strains at the selected sensor points are post-processed. The numerical process of 3-D DBEM have been implemented in the KSP software, developed by H. Kebir from the University of Technology of Compiegne. The software uses C++ programming language to solve the direct problem consisting of computing stresses and strains on the domain frontiers for a cracked structure subjected to a mixed field of boundary conditions (displacements and tractions). In this present study, we have implemented both the Low-Discrepancy Sequence (LDS) algorithm and the Nelder-Mead Simplex Algorithm (NMSA) to solve the inverse problem of crack identification. The inverse problem uses the computed (or measured) strains at some sensor points from a real case of cracked structure and minimizes the regularized objective function Eqn. 3.

NELDER-MEAD OPTIMIZATION ALGORITHM
ased on the works of Box [25] and Spendley [26], the Nelder-Mead [17] version of the Simplex Algorithm (NMSA) is an optimization method applied to construct a non-local linear approximation of a n-dimensional constrained problem from a set of points at sufficient intervals. The method finds a locally optimal solution to an objective function (OF) in  n when it varies smoothly. At each iteration, the NMSA produces a new test position using four operations named: reflection, expansion, contraction and the similarity transformation. Through these operations, the algorithm extrapolates the behavior of the OF measured at each test point and arranged as a simplex. The algorithm then rejects the position which gives the highest value of the OF with the new test point and so the algorithm progresses [27]. The first step is to change the worst point with a point reflected through the centroid of the opposite side of the Simplex. If this point is better than the best current one, then we can attempt extending exponentially out along this line. Now, if this new point do not give a much better value for the OF than the previous value, then we contract the simplex towards the best point (Fig. 3). Additional constraints must be considered on domain limits to ensure embedded cracks, on the first guess for the crack's identity to start the minimization process, which must converge to the actual solution as well as a convergence test to stop the progression of the algorithm. The first 3D example treated is relating to a cylindrical shaft in axial traction with an embedded elliptical crack. The dimensions used for the shaft are L/D=2; L and D, are the length and the diameter of the shaft respectively; a/b=2 and D/a=10, a and b the crack radii. The material has an elasticity modulus E= 200 GPa and a Poisson's ratio ν=0.3. The NMSA was applied starting with an initial crack position close to the actual one to carry out a fast convergence and avoiding local minima. Convergence to the actual position of the crack (Fig. 4) was achieved in less than 1100 iterations with a normalized OF less than 10 -5 . Figs. 5-7 show the convergence of the normalized OF and the crack's identity parameters to B the actual position of the crack. The fluctuating shape of the curve is a consequence of the application of the NMSA steps described in Fig. 3.  When the structure under analysis is a non-convex domain or includes cavities and other geometry irregularities, the NMSA does not converge toward the required solution and often converges to a local minimum. Indeed, the convergence of the algorithm strongly depends on the choice of the initial identity of the crack [10]. The NMSA is an algorithm with a slow convergence, then more the initial identity is close to the actual one, better is the probability to avoid local minimum. Since we do not know a priori where to look, the idea for selecting the initial crack identity to start the search algorithm for the actual identity is to breed a random sequence of identities for the crack and to hold the one with the smallest value of the objective function. Selected identity parameters are used as an initial crack identity to start the NMSA algorithm. Of course, the random crack generation procedure must uniformly cover the entire area of interest. In order to achieve this goal, we use a quasi-random sequence that has a high speed of convergence even if the stochastic independence of draws of the same element is not observed. The independence of the elements is the main condition to respect.

LOW-DISCREPANCY SEQUENCE (LDS)
low-discrepancy sequence still called a quasi-random sequence of cracks' identities fills the geometry more uniformly than uncorrelated random identities. Although the ordinary uniform random numbers and quasi-random sequences both produce uniformly distributed sequences, there is a major difference between the two. A uniform random generator will produce outputs so that each trial has the same probability of generating an identity on equal subspaces. Thus, it is possible for n trials to coincidentally all lie in the first half of the geometry, while the (n+1)st crack still falls within the other half. This is not the case with the quasi-random sequences, in which the outputs are constrained by a low-discrepancy requirement that has a net effect of identities being generated in a highly correlated manner (i.e., the next generated crack "knows" where the previous cracks are). In order to offer the NMSA its best starting solution, a random sample of cracks' identities is produced using a quasirandom sequence. Afterwards, the identity which gives the smallest value of the OF, given by Eqn. 1 is selected and used as initial solution to start the minimization process.  The random generation of crack identities may produce cracks that lie out of the domain's boundaries. In order to retain only embedded cracks restrictions are made to exclude invalid locations within the domain of interest. The sequence is generated by the series [28] k is a succession of integers and p is the kernel of the LDS. There is as much distinct p than random variables. In this sequence, p is a prime number and the brackets indicate the integer part of the product. Fig. 7 shows the scatter plot of the first 100 points with the kernel p taking the values of the first four prime numbers. The same example of the cylindrical shaft in Fig. 4 is taken again. The quasi-random sequence is used to generate a sample of 50 crack identities and we can note the regularity of the distribution witch cover all regions of the domain (Fig. 9). Strains are calculated and the OF is evaluated at the end of each iteration. In the second step, NMSA is applied starting the perturbation by the best random approximation of the crack's identity. Convergence is carried out with less than 1100 iterations with normalized OF less than 10 -5 (Fig. 10). Figs.11-13 show respectively the convergence of the OF and the crack's identity parameters to the actual position. As can be seen in Fig. 11 the error on the OF is less than 10 -4 after only 600 iterations. This represents about half of the total number of iterations and in the second half, the graphs show small fluctuations around the actual solution.      . In order to assess the stability of the regularized crack identification algorithm to measurement errors, a Gaussian law is used to generate a fixed standard deviation noise level, which is added to inputs as perturbations and the deviations for the solutions are compared to the noise level. The same crack configuration and loading conditions as in the above example were used. Tests performed with noisy data show that deviations obtained for the cracks' identities confirmed the robustness of the coupled LDS-NMSA algorithm in 3D crack identification (Fig. 14). The second example, illustrated by Fig. 15, deals with the identification of a crack in a spur gear tooth of a large girth gear used in driving heavy industrial equipment like tube mills and rotary kilns. The Fig. 15 gives a drawing of a gear tooth of 40 mm module with dimensions, loading and an inner crack to be identified. The results are obtained for carbon-steel gear material with an elasticity modulus E=210 GPA and a Poisson's ration ν=0.3.
The optimal random crack's identity is obtained with a sequence of 50 iterations (Fig. 16). The final position was reached after 1200 NMSA iterations with a normalized OF less than 10 -5 (Fig.17). Figs. 18-20 show respectively the convergence of the OF and the identity crack parameters to the actual position for the spur gear tooth example. Figure 15: Elliptical crack in a spur gear tooth. Figure 16: LDS generation of cracks' identities (Spur gear tooth example).
In Tab.1, the first line of values gives the eight real identity parameters used to identify an elliptical crack in a spur gear tooth. The second line gives the optimal crack identity (with the lowest value of the OF) among fifty random LDS iterations. From the third line and down, results for the convergence process by the NMSA algorithm are given, and we can see that for a value of the objective function fixed at 10 -4 , convergence is reached after 1100 iterations.

CONCLUSIONS
three-dimensional crack identification method in the context of linear fracture mechanics has been developed. The inverse problem of crack identification has been defined as a minimization of a least squares functional given as an objective function. The direct problem was solved by the Dual boundary Element Method to compute strains at selected sensor points on the boundaries, since strains allow to better capture the effect of crack opening on the domain boundaries, by a judicious fixture of elastic strain gauges. The crack surface is parameterized by an identity vector, which groups shape and position parameters. The direct search method used to minimize the objective function is the Nelder-Mead Simplex Algorithm, which is a heuristic search method that can converge to non-stationary solutions. To overcome this shortcoming, a low-discrepancy sequence has been used to generate a sample of crack identities, uniformly distributed over the entire domaine of interest. From the generated sample, the crack identity with the smallest value of the objective function is held as a starting solution for the NMSA, which refines the search process to get the nearest solution to the actual one.
Two numerical examples of identification of elliptical cracks were treated. The first one dealt with a cylindrical shaft in traction. In this example, the stability of the regularized crack identification algorithm to measurements errors is assessed with introduction of noise via a Gaussian law with fixed standard deviation level. The perturbations are added to input strains, and the deviations for the solutions are compared to the noise level. Tests performed with noisy data has exhibited a stable convergence with deviations in cracks identities of the same order as those used for noises. The second example concerned a spur gear tooth of a large girth gear. Despite the relatively important number of iterations related to the number of identity parameters to be identified, the handled examples has shown that the coupling of the two methods provides reliable identification results.