Optimisation Method for Determination of Crack Tip Position Based on Gauss-Newton Iterative Technique

In the digital image correlation research of fatigue crack growth rate, the accuracy of the crack tip position determines the accuracy of the calculation of the stress intensity factor, thereby affecting the life prediction. This paper proposes a Gauss-Newton iteration method for solving the crack tip position. The conventional linear fitting method provides an iterative initial solution for this method, and the preconditioned conjugate gradient method is used to solve the ill-conditioned matrix. A noise-added artificial displacement field is used to verify the feasibility of the method, which shows that all parameters can be solved with satisfactory results. The actual stress intensity factor solution case shows that the stress intensity factor value obtained by the method in this paper is very close to the finite element result, and the relative error between the two is only − 0.621%; The Williams coefficient obtained by this method can also better define the contour of the plastic zone at the crack tip, and the maximum relative error with the test plastic zone area is − 11.29%. The relative error between the contour of the plastic zone defined by the conventional method and the area of the experimental plastic zone reached a maximum of 26.05%. The crack tip coordinates, stress intensity factors, and plastic zone contour changes in the loading and unloading phases are explored. The results show that the crack tip change during the loading process is faster than the change during the unloading process; the stress intensity factor during the unloading process under the same load condition is larger than that during the loading process; under the same load, the theoretical plastic zone during the unloading process is higher than that during the loading process.


Introduction
Fracture failure of a component originates from the initiation and propagation of cracks, and the arrest or growth of cracks during their development is governed by the stress field near the crack tip. The stress intensity factor (SIF) K, which characterises the stress field strength near the crack tip, can be used to estimate whether the crack will become unstable [1]. There are many ways to measure the SIF experimentally, such as strain-based techniques [2,3] and photoelastic methods [4]. However, the photoelastic method is usually applicable only to transparent materials. Since digital image correlation (DIC) technology was introduced in the 1980s, it has been widely used to measure the crack tip displacement field because it is a highly accurate non-contact method that can measure the full field in the near-tip region. In the crack tip field function, the crack tip is used as the coordinate origin. When the displacement field obtained by DIC technology is used to measure the SIF, it is important to accurately determine the crack tip position. Avtaev et al. [5] pointed out that it is difficult to directly identify the position of the crack tip from digital images, and scholars have proposed many methods for this problem. In the area near the crack tip, the growth rate of elastic and plastic strain will increase, and will reach the maximum at the crack tip. For I model cracks, the distribution law of the e yy strain value on the extension line of the crack can be used to identify the location of the crack tip, Avtaev and Yakovlev [6] use this method to determine the location of the crack tip. The displacement value is the first-hand information obtained by using DIC technology, and it is a common method to obtain the position of the crack tip by using the displacement field data. One method to achieve this is to observe whether the vertical displacement field near the crack tip is continuous and use the critical point at which the displacement suddenly changes as the crack tip position [7]. Vasco-Olmo et al. [8][9][10] used the projection of the vertical displacement field on the VOY coordinate plane to find the coordinates of the crack tip in the Y direction. The projection of the vertical displacement field on the VOX is combined with the displacement value corresponding to the coordinate of the crack tip in the Y direction to determine the coordinate value of the crack tip in the X direction. On the other hand, some scholars directly use the method of image recognition to obtain the position of the crack tip from the digital image. Gao et al. [11] locate the sub-image containing the crack from the crack tip image according to the image grayscale standard deviation (GSD) distribution, and then use the maximum entropy threshold segmentation method to extract the crack contour from the sub-image. Then, the crack skeleton image is extracted from the crack contour image by the morphological thinning algorithm to determine the position of the crack tip. This method is more complicated. Mokhtarishirazabad et al. [12] use a simpler method. It directly enlarges the digital image of the crack tip area to a certain multiple to determine the position of the crack tip. Cao et al. [13] used a similar method to determine the position of the crack tip of the rock sample. When Yue et al. [14] is studying the dynamic fracture parameters of marble, it directly uses two high-speed cameras to record the crack propagation process to determine the position of the crack tip. Tong [15] and Zhu et al. [16] used an optical microscope to determine the coordinates of the crack tip in the x direction. After the coordinates of the crack tip in the x direction are determined, the average value of the vertical displacement field is used to determine the coordinates of the crack tip in the y direction. This method is simple in mathematics, but requires equipment other than DIC (optical microscope). The artificial selection of the crack tip position with the help of sophisticated equipment may lead to differences in the crack tip positions selected by different researchers, and the use of mathematical methods combined with displacement field data can avoid this difference. Yoneyama et al. [17] used the Newton-Raphson iteration method with the crack tip position as an unknown parameter in the solution. Zanganeh et al. [18] suggested that the iteration matrix of this method might be ill-conditioned and that the initial iteration value pair would greatly affect the convergence. In order to avoid solving the nonlinear problem encountered in obtaining the crack tip position, Zhang et al. [19] set up a number of points in the crack tip neighborhood, and used these points as the crack tip coordinates to calculate the objective function value. The point corresponding to the minimum value of the objective function is the final coordinate of the crack tip. This paper proposes a crack tip location optimisation method based on the Gauss-Newton nonlinear iterative method. When this method is combined with observations of the discontinuity of the vertical displacement field near the crack tip, the crack tip position can be determined more accurately. This provides a solution for the iterative initial value selection of the nonlinear iterative solution of the crack tip position and the parameters of the Williams equation. At the same time, the preprocessing conjugate gradient method is used to solve the ill-conditioned iterative matrix encountered in the iterative process. This method provides more options for the subsequent description of the crack tip plastic zone and calculation of the SIF. Similar to the common use of nonlinear optimisation to solve the parameter methods such as crack tip position, this article assumes that the actual displacement field measured each time is independent, and the accuracy of the solution to the displacement field is independent each time. But Fernández [20] believes that this condition is not necessarily satisfied, and he suggests using the generalized least squares method to solve the above parameters. This provides guidance for the next step of research.

Crack Tip Field Function
In a paper published by Williams in 1957, the crack tip displacement field and stress field in a plane stress or plane strain state were expressed as an infinite series solution [21]. Although this solution is not a closedform solution, it has a broad range of applications. This series expansion is widely used to calculate the SIF of many types of samples commonly used to study fracture mechanics [22]. Double edge notched test (DENT) samples were used in this study and the crack can be treated as a mode I crack. The displacement and stress fields are expressed as follows [23]: where a 1 = K I / √ 2π , a 2 = σ 0x 4 , G is the shear modulus, κ = (3 − υ) (1 + υ) for the plane stress condition, and κ = 3 − 4υ for the plane strain condition, where υ is Poisson's ratio. The sample is rotated and translated during loading. The displacement field measured directly by DIC technology can be expressed as follows: where T x , T y , and R represent translation in the x direction, translation in the y direction, and rotation of the sample, respectively. Because Eq. (3) uses the crack tip position as the coordinate origin, the following (1a) u = ∞ n=1 r n/2 2G a n (κ + n 2 transformations must be made in order to use the results obtained by the DIC measurement system: where (x i , y i ) are the coordinates of a point in the DIC measurement system, and (x o , y o ) are the coordinates of the crack tip in the DIC measurement coordinate system.

Crack Tip Plastic Zone
When the equivalent stress is equal to the yield strength, metals experience plastic deformation. In addition, the SIF is a key parameter in linear elastic fracture mechanics. Therefore, if the displacement field is used to calculate the SIF in a test, the plastic region must be excised from the data by defining an inner radius for data collection that excludes the plastic zone. Parnas et al. [24] defined this inner radius of the fitting area as R inner > T 2 , where T is the sample thickness. However, the range of the plastic zone sometimes exceeds this inner diameter and hence, in either plane stress or plane strain, Dehnavi et al. [25] suggested y ] to determine the inner radius of the fitting area. However, the SIF range under given experimental conditions must be obtained before this formula can be applied, which makes this method not always applicable.
The displacement field measured by DIC can be used to obtain the strain field through numerical calculations. By combining the strain field with appropriate material parameters, the stress field and crack tip plastic region can be obtained. The strain field must be as accurate as possible and Pan et al. [26] proposed a full-field strain measurement method based on local least-squares fitting of the displacement field; because this method has strong noise reducing capability, it is much more accurate than the difference method. To obtain the strain field from the displacement field, Vasco-Olmo et al. [27] used the Green-Lagrange strain tensor, which is expressed as follows: where u p and v p are the local horizontal and local vertical displacement fields obtained using local least-squares fitting, respectively. The stress field can then be obtained from the strain field as follows: To determine the extent of the plastic zone, the von Mises stress or Tresca stress is used in combination with the yield strength of the material. In this study, the von Mises stress is used to identify the plastic zone.

Gauss-Newton Iterative Method
The inner radius of the region used for the iterative process can be calculated according to Section 2.2. The outer diameter of the fitting area and the number of Williams expansion terms used during the fitting process are determined according to the conditions. Let the horizontal and vertical displacements corresponding to the coordinates (x i , y i ) of the displacement field in the fitting area be U i and V i , respectively. The corresponding objective function is: where Z = a 1 , a 2 , · · · , a l , x o , y o , T x , T y , R , a i is the Williams expansion coefficient, and the Z value at which ψ(Z) reaches the minimum is the optimal solution of the objective function. Both the horizontal and vertical displacement fields contribute to the fitting: The iterative form of the Gauss-Newton method is [28]: is a Jacobi matrix, and the mathematical expression is a first-order partial derivative matrix.
T , where l is the number of expansion terms in the Williams equation, and z i corresponds to the above parameters. DF (Z) can be written as follows according to the derivation rule: Note that because of the nonlinear nature of the matrix components, sometimes the components are correlated, which causes DF (Z) to become degenerate or ill-conditioned, i.e., more sensitive to small errors, and DF (Z k ) T DF (Z k ) is also ill-conditioned [28]. In addition, the DIC technique measures the crack tip displacement field in units of pixels, with a single pixel generally equating to a distance between several microns and several tens of microns. Thus, DF (Z k ) T DF (Z k ) is always ill-conditioned in the actual calculation. There are two ways to solve this problem: one is to optimise the Gauss-Newton method itself, or the more commonly used damped least-squares method (the Levenberg-Marquardt method) and its variants [28,29]; the other is to improve the accuracy with which ill-conditioned equations are solved. The preconditioned conjugate gradient (PCG) method can be used to solve a system of equations whose coefficient matrix is an ill-conditioned sparse symmetric matrix.
is a non-sparse symmetric matrix, and the sparsity of the matrix is not a necessary condition for applying the PCG method in MATLAB. The crack tip optimisation iterative algorithm can be written as shown in Figure 1. Figure 1 can be run according to the following steps: a) Determine the initial crack tip position (x o , y o ) according to the discontinuity of the crack tip displacement filed. b) Using (x o , y o ) as the center, select a fitting area radius r. Take the fitting point in the displacement filed after removing the plastic area. c) Using the linear fitting method and (x o , y o ) , Z ′ = (a 1 , a 2 , · · · , a n , T x , T y , R) can be obtained.
f ) Judge whether || Z|| 2 satisfies the convergence condition, if so, output Z. g) If || Z|| 2 does not meet the convergence condition, perform Z = Z + Z , extract new crack tip coordinates, and update the fitting area. h) Repeat steps e), f ), g) until the convergence condition is meet, and Z is output.

Numerical Verification
The feasibility of the proposed crack tip optimisation algorithm can be confirmed using a simulated displacement field that includes translational and rotational components, and noise, whose parameters are given in Table 1.
The elastic modulus E is set to 210000 MPa, Poisson's ratio υ is 0.28, and the yield strength σ s is 780 MPa.  Figure 2 shows the size and shape of the plastic zone obtained by local least-squares fitting   according to Pan et al. [26] and the effect of increasing noise is clearly seen in Figure 1b. As shown in Table 2, the proposed algorithm can calculate the parameters of an artificial displacement field with noise.

Materials and Testing
The experimental work was performed on U71MnG rail steel. Its mechanical properties and chemical composition are given in Tables 3 and 4, respectively. DENT samples were used with thickness of 1 mm and the geometry shown in Figure 3. Fatigue crack growth testing was performed on an ElectroPuls E3000 dynamic testing machine. The test parameters are summarized in Table 5. During the experiment, the crack tip displacement field was captured using a Dantec Q-400 3D DIC system. As an example, Figure 4 shows cloud diagrams of the displacement fields of the sample under the maximum load for a crack size of 18 mm. Both the horizontal and vertical displacement fields are symmetric with respect to the crack surface as a whole, but the symmetry is accompanied by translation and rotation which is considered in Eq. (3) by adding translation and rotation terms.

Description of Displacement Fields and Plastic Regions
The calculation of the SIF from the displacement field also requires knowledge of the values of two important parameters: the outer radius of the fitting region and the number of expansion terms used in the Williams equation. The outer radius of the fitting area is usually a function of the number of fitting terms. For a larger fitting region, more terms are required and the calculation effort increases [23]. Mokhtarishirazabad et al.    [23] used 20 mm×20 mm squares as the fitting regions and noted that when the number of expansion terms is five, a stable SIF solution can be obtained. Ju et al. [30] determined the length of the crack as the outer radius of the fitting area. In this case, a stable solution can be obtained when the number of fitting terms is eight. A circular fitting region was used in the work described in this paper with an outer radius of R outer = 2.5 mm, and a stable convergence solution was also obtained when the number of fitting terms was five. The relative error of each point in the fitting area is defined as η = (ω ′ i − ω i )/ω i where ω ′ i is u ′ i or v ′ i and refers to the x i , y i displacement of the point obtained using the Williams equation; ω i is u i or v i and refers to the displacement of the point x i , y i measured using the DIC technique. To simplify the comparison, in the present paper, the method of using the discontinuity of the displacement field to determine the crack tip position is called the conventional method. Now, this paper uses the crack length of 18 mm and the DENT sample loaded to 2400 N as a typical working condition. Table 6 shows the results of statistical analysis of the fitting results for the conventional and proposed methods, specifically, the relative errors at corresponding points in the respective fitting regions. The mean and variance of the relative errors both indicate that the proposed method results in a better fit than the conventional method, and that the resulting displacement field is closer to the true values measured by the DIC technique.
A finite element technique was combined with displacement extrapolation method and used to solve the numerical solution for the SIF. The finite element model was loaded using a coupled loading method that reduced the relative error and achieved a lower calculation cost [31]. Now, we use the above conditions as a case to compare the results of SIF. Table 7 shows that the SIF obtained after the crack tip position was optimised is closer to the finite element solution with an relative error of only −0.621% than that obtained using the conventional method which gives an relative error of 7.21%. This further demonstrates that crack tip optimisation is necessary and effective.
From the Williams coefficients and crack tip position parameters obtained by the proposed and conventional methods, the contour lines of the plastic zone can be drawn and compared with that of the plastic zone   calculated directly using the displacement field ( Figure 5). The crack lengths corresponding to the three working conditions shown in Figure 5 are all 18 mm. Table 8 lists the relative error values of the plastic zone area obtained by the two methods and the experimentally calculated plastic zone area. Obviously, during loading, at the peak load, and during unloading, the contour line of the Williams equation obtained by the proposed method reproduces the crack tip plastic area and its changes during the test better than the conventional method does. It can be seen from Figure 5 and Table 8 that under the three working conditions, the area of the plastic zone calculated by the method proposed in this paper is smaller than the experimental result, while the result of the plastic zone obtained by the conventional method is larger than the experimental result. Table 9 lists the Williams parameters and the crack tip coordinates for the three working conditions. The x 0 coordinate value obtained by the proposed method in this paper is smaller than the x 0 value of the conventional method. In fact, the value of x 0 can directly reflect the length of the crack. As far as this article is concerned, the smaller the value of x 0 corresponds to the longer the crack. Among all the coefficients of Williams equation, |a 1 | and |a 2 | significantly affect the size of the theoretical plastic zone, and the two correspond to K and T respectively. From the data in Table 9, it can be found that under the same working conditions, the |a 1 | and |a 2 | obtained by the proposed method are significantly larger than the |a 1 | and |a 2 | obtained by the conventional method. Under the same displacement field data, the longer the crack corresponds to the smaller K and T, so in this paper, the smaller the value of x 0 , the smaller |a 1 | and |a 2 |. The conventional method cannot consider the influence of crack tip dullness when selecting the crack tip coordinates. At this time, the crack length calculated using the selected crack tip coordinates will be too small. According to the above analysis, it can be seen that at this time |a 1 | and |a 2 | will be too large, resulting in a larger theoretical plastic zone. Due to crack tip dullness and stress relaxation, the crack tip coordinates obtained by the method proposed in this paper are used to calculate the crack length will be long. According to the above analysis, the smaller |a 1 | and |a 2 | will be obtained at this time, so the plastic zone is also smaller.

Application
The analysis presented in Section 4.2 shows that the proposed method outperforms the conventional methods of calculating important parameters, such as the crack tip position, Williams coefficients, and the translational and rotational components of the actual displacement field. The proposed method was used to calculate the crack tip position and corresponding SIF during loading and unloading of a DENT sample of U71MnG rail steel for a certain load cycle and to characterize the crack propagation during this cycle. To obtain clear results, a load cycle with a long crack (18 mm) near the end of the test is selected. The loads are shown in Table 10, where the coordinates of the crack tip under a 1320 N load are used as the origin. Figure 6 shows the crack propagation path, while Figure 7 shows the SIF during loading and unloading in the chosen load cycle.
In Figure 6, the calculated crack tip position appears to move forward during the loading and unloading, and the crack tip moves forward in a zigzag manner. During the loading phase, from 1320 N to 2400 N, the crack tip moved 11.82 pixels forward. During the unloading stage from 2400 N to 1320 N, the crack tip moved 5.3 pixels forward. Generally speaking, the growth rate of fatigue cracks is much lower than the rate calculated in this paper. When the load increases, the plastic zone will increase, leading to more pronounced stress relaxation. The result is that the crack tip position is calculated from the displacement field, and it is estimated that the change of the crack length will show an increase in the crack length. This is the reason why the crack growth rate calculated in this paper is higher than the conventional fatigue crack growth rate. Due to the stress relaxation of the crack tip and the actual expansion of the crack, the SIF when unloaded under the same load level is slightly larger than the SIF during the loading process. Figure 8 shows the evolution of the theoretical plastic zone based on the Williams equation for the above load cycle. As the load increases in each loading cycle, the plastic area gradually increases, reaching its maximum value at the maximum load of 2400 N and then gradually decreasing with decreasing load. In addition, owing to the increases in the crack size and SIF during loading at these high crack growth rates, the plastic area during unloading becomes slightly larger than that found during loading at the same load level.
The characteristic horizontal size of the plastic region can be defined as the intersection of the maximum contour line with the x axis [22] and the maximum extent of the plastic region is then given by r x = (K I /σ s ) 2 /(2π) . The vertical extent is defined as the maximum vertical distance of the contour line of the plastic area from the x axis. When the Williams series expansion uses only the first term, the contour line of the plastic zone is expressed as follows [22]: In this case, the lateral feature size is r x = (K I /σ s 2 (2π) . When the expression: (15) f (θ) = sin θ cos 2 θ 2 1 + 3 sin 2 θ 2   has a maximum value, the vertical feature size r y can be obtained; the value of this function is near the maximum value at θ = 1.3984 and, r y = 0.20598(K I /σ s ) 2 . The characteristic size of the plastic region in this case clearly depends only on the SIF; that is, it depends only on the first coefficient of the Williams series. The characteristic size of the plastic zone is called the first-order feature size. When the number of terms n is expanded (n > 1), the characteristic size of the plastic region is called the nth order feature size. Generally, numerical solutions can only be obtained for the nth order feature size. Figure 9 shows the lateral feature size for loading and unloading under the above loads. During loading and unloading, the fifth-order lateral feature size has a quadratic relationship with the load. In addition, the fifth-order lateral feature size is smaller than the first-order lateral feature size, and the difference between the two increases with increasing load. Take loading process for example, the fifth-order lateral feature size is 9.018% smaller than that of the first-order lateral feature size at 1320 N, while the  difference reaches 18.75% at 2400 N. Figure 10 shows the corresponding vertical feature sizes for loading and unloading. Different from the lateral feature size, the difference between the first-order and fifth-order vertical feature sizes is rather small. Both the first-order and fifth-order vertical feature size have a quadratic relationship with the load. In general, the higher-order Williams coefficient strongly affects the lateral feature size, and this effect increases with increasing load. However, the higher-order coefficient has a weaker effect on the vertical feature size, which is essentially independent of the load. For present study, the SIF calculated at 2400 N using the first-order Williams equation is 56.71 MPa·m 0.5 . As described in Section 4.2, when using the expanded fifth-order Williams equation to calculate SIF, the value is 58.46 MPa·m 0.5 . Obviously, the SIF value obtained by using the first-order Williams equation is slightly smaller.

Conclusions
(1) This paper proposes an optimisation method for crack tip position based on Gauss-Newton iterative method. The linear fitting method is used to provide the initial value of the nonlinear iteration, and the preconditioning conjugate gradient method is used to solve the ill-conditioned equations.
(2) The relative error between the SIF value obtained by the proposed method and that by the finite element method is only -0.612%, which is much smaller than the error, 7.21%, generated between the conventional method and the finite element method. Moreover, the difference between the crack tip plastic zone obtained by the method proposed in this paper and the experimental plastic zone is even smaller, and the maximum is only -11.29%, which proves that the method proposed in this paper is effective.
(3) The crack propagation path, SIF, and plastic zone evolution during loading and unloading in a load cycle have been determined using this method. The results show that the crack propagates faster during loading than during unloading, and the crack tip moves forward in a zigzag pattern. The calculated crack growth rate includes the result of crack tip blunt and stress relaxation. As the crack continues to grow, during the un-loading process with the same load as the loading process, the SIF and plastic zone at the crack tip become larger, which provides a basis for further understanding of crack propagation behavior.
(4) The higher-order coefficients of the Williams equation strongly affect the lateral feature size of the plastic zone, and this effect increases with increasing load. However, the higher-order coefficients have a weaker effect on the vertical feature size. These characteristics are one of the reasons why different Williams equation orders are used to obtain different SIF values when solving the stress intensity factor using the least square method. Yang