Fracture and Fatigue Analyses of Cracked Structures Using the Iterative Method

It is a quite challenging subject to efficiently perform fracture and fatigue analyses for complex structures with cracks in engineering. To precisely and efficiently study crack problems in practical engineering, an iterative method is developed in this study. The overall structure which contains no crack is analyzed by the traditional finite element method (FEM), and the crack itself is analyzed using analytical solution or other numerical solutions which are effective and efficient for solving crack problems. An iteration is carried out between the two abovementioned solutions, and the original crack problem could be solved based on the superposition principle. Several typical crack problems are studied using the present method, showing very high precision and efficiency of this method when making fracture and fatigue analyses of structures.


Introduction
The stress intensity factors (SIFs) characterize the singularity of the area near the crack tip, which are the most significant parameters for fracture mechanics. However, the lack of precise solutions of stress intensity factors has hindered the progress and application of fracture mechanics to the fracture and fatigue analyses for various kinds of structures. Once the SIFs are precisely and efficiently computed, the fatigue crack propagation rate of a cracked structure under cyclic load can be determined, and its fatigue life estimation can then be rationally made. So, the calculation of fracture mechanics parameters for arbitrary surface and embedded cracks in complex civil, mechanical, and aerospace structures remains an important task for the structural integrity assessment and damage tolerance analyses [1]. Since analytical solutions are usually very difficult, sometimes even impossible, to be obtained for crack problems in complex structures, it is of great importance to conduct numerical as well as experimental studies on crack-related problems in practical engineering.
Since the traditional finite element method (FEM) uses simple polynomial interpolations in numerical analysis, it is unsuitable and unwise to simulate cracks and their fatigue growth with FEM, which is largely attributed to its high inefficiency and labor costing of approximating stress and strain singularities using polynomial FEM shape functions. Then, embedded singularity elements [2,3] and singular quarterpoint elements [4,5] were proposed to overcome this difficulty, which are now incorporated in many commercial FEM software. However, the need for constant remeshing makes the automatic fatigue crack growth analyses with FEM extremely difficult, sometimes even impossible.
To overcome the inherent difficulty of FEM when solving crack problems, many other numerical methods have been proposed and studied during the past several decades. The traditional boundary element methods (BEM) and dual boundary element methods were developed in fracture analysis [6,7]. Compared with the traditional and dual boundary element methods, the symmetric Galerkin boundary element method (SGBEM) showed several advantages, which resulted in a symmetrical coefficient matrix of equation system, and avoided to specially treat sharp corner. Early derivations of SGBEMs involved regularization of hypersingular integrals [8][9][10][11]. Han and Atluri [12,13] proposed a systematic procedure to develop weakly singular symmetric Galerkin boundary integral equations. The derivation of this simple formulation involved only the nonhyper singular integral equations for tractions [14,15]. Since the coefficient matrix for SGBEM is usually fully populated, it is unwise to carry out large-scale simulation of complex structures with pure SGBEM. In order to explore the advantages of both FEM and SGBEM, they were coupled indirectly to analyze cracked 3D solids with surface flaws [16,17], which was called the SGBEM-FEM alternating method. Tian et al. [18][19][20] conducted three-dimensional fracture and fatigue analyses of many typical structure components in civil and mechanical engineering by employing the SGBEM-FEM alternating method, and obtained a series of numerical results with reference value for practical engineering. The extended finite element method (XFEM) has been widely used for fracture and fatigue analysis in engineering. The fatigue life of homogeneous plate containing multiple discontinuities was evaluated by XFEM under cyclic loading condition by Singh et al. [21]. Bergara et al. [22] presented the numerical simulation and validation of a fatigue propagation test of a semielliptical crack located at the side of the rectangular section of a beam subjected to four-point bending, and for the numerical simulations, XFEM implemented in the Abaqus (R) 2017 software has been used. Liu et al. [23] investigated the influence of the microdefects on the macrocrack propagation numerically, by using the multiscale method in conjugation with XFEM. This kind of the multiscale method can efficiently provide solution of accuracy on both scales without any restrictive assumptions as the analytical approximation method. To evaluate the fatigue life of cracked specimens in the presence of bimaterial interfaces, Jameel and Harmain [24] modified the standard finite element approximation by adding appropriate enrichment functions in order to consider the effect of these discontinuities present in the domain, and the level set method was used to track different discontinuities in the domain. Many other applications of XFEM to fracture and fatigue analyses can be found in the papers [25][26][27][28][29][30]. Besides these methods, the meshless method (or mesh-free method) and the numerical manifold method (NMM) were also used for analyzing crack problems. Rao and Rahman [31] presented an efficient meshless method for analyzing linearelastic cracked structures subject to single or mixed-mode loading conditions. The method involved an element-free Galerkin formulation in conjunction with an exact implementation of essential boundary conditions and a new weight function. Ching and Batra [32] used the meshless local Petrov-Galerkin (MLPG) method to determine the crack tip fields in linear elastostatics, by studying a single-edge cracked plate and a double-edge cracked plate with plate edges parallel to the crack axis loaded in tension. New enriched weight functions were introduced by Duflot and Nguyen-Dang [33] in the meshless method formulation to capture the singularity at the crack tip and simulate the fatigue growth of cracks in two-dimensional bodies. The efficiency of the classical truly meshless local Petrov-Galerkin method with linear test function approximation in the crack growth problems of complex configurations was investigated by Memari and Mohebalizadeh [34]. Rao and Rahman [35] presented a coupling technique for integrating the elementfree Galerkin method (EFGM) with the traditional finite element method (FEM) for analyzing linear-elastic cracked structures subject to mode I and mixed-mode loading conditions. The EFGM was used to model material behavior close to cracks and the FEM in areas away from cracks. Based on the contact technique of the numerical manifold method and the incorporation of the Mohr-Coulomb crack initiation criterion, the effects of the friction and cohesion on the crack growth from a closed flaw (crack) under compression were investigated by Wu and Wong [36]. By employing NMM, Ma et al. [37] studied complex crack problems such as multiple branched and intersecting cracks to exhibit the advantageous features of the NMM. Zhang et al. [38] made numerical analysis of 2D crack propagation problems using the NMM and focused on modeling complex crack propagation problems containing multiple or branched cracks. Zheng et al. [39] generated the mathematical cover with the influence domains of nodes in the moving least squares interpolation instead of commonly used finite element meshes, significantly simplifying the generation of the physical cover and the simulation of crack growth. A three-node triangular element fitted to the NMM with continuous nodal stress called Trig3-CNS (NMM) element for accurately modeling two-dimensional linear elastic fracture problems was presented by Yang and Zheng [40]. Many other applications of the meshless method and NMM to fracture and fatigue analyses can also be found in the literature [41][42][43][44][45]. Besides, many other researches [46][47][48][49][50][51][52][53] were conducted by numerical simulation or model test to study the fracture evolution and failure process of various kinds of engineering problems.
In this paper, an iterative method is employed to obtain the SIFs of structures with different crack configurations and examine the fatigue crack growth subjected to cyclic loading. Crack growth rates are determined by the Paris fatigue law. For damage tolerance evaluation, the crack growth paths and number of loading cycles are predicted. The validations of the current iterative method are illustrated by the comparison of numerical results with available analytical solutions as well as empirical observations. The finite element method is generally unsuitable for modeling cracks, because polynomial shape functions are inefficient for approximating singularities, and the shaped function-based method involves complicated remeshing procedures. On the other hand, the iterative method is a highly accurate and efficient method for solving fracture mechanics problems. The basic idea is to model the uncracked global structure using finite elements with very crude meshes, and model the crack with more suitable methods, such as the analytical methods, integral equations, and SGBEMs. By iterating between the solution for a crack in an infinite domain (or in a finite subdomain) and the solution by finite element for the uncracked global structure, the iterative methods obtain solution for the original problem.
Tractions at the crack's surface are obtained by FEM, and the stresses at field boundaries are corrected through the iteration process. No remeshing is involved, which is a significant advantage compared to the traditional FEM-based methods.

Principle of Superposition.
For linear elastic materials, stress intensity factors are additive like stress, strain, and displacements, as long as the mode of loading is consistent, which can be expressed as follows: So, tractions acting on the crack face can replace tractions that act on the boundary, as long as these two configurations result in the same SIFs, as shown in Figure 1.
2.3. Algorithm of the Iterative Procedure. Consider a finite body containing initial cracks. The crack surface with no traction is denoted as B c . Boundary of the finite domain is denoted as B, in which the boundary with prescribed tractions t 0 is B t , and the boundary with prescribed displacements u 0 is B u . Obviously, B = B u S B t . The iterative method involves two different solutions for solving the original problem. The first is denoted as P A (shown in Figure 2(c)), which is an infinite domain contain-ing cracks subjected to unknown crack surface loading T. The second one shown in Figure 2(b) is denoted as P F , which has the same geometrical properties as the original problem but without cracks. For P F , its boundary B u has the prescribed displacement u, and boundary B t has the prescribed traction t. Since cracks are ignored in P F , u and t are different from those for the original problem. The problems P F can be easily solved by the finite element method. Here, the original problem is shown in Figure 2(a), which is denoted as P O . By iteration of P A and P F , the solution of P O can be obtained.
For the problem of P F , the finite element method can be used to obtain the tractions at the location of crack with any given boundary conditions u and t, which can be expressed as follows: where K u 1 and K t 1 are the stiffness matrices, which are linear operators.
For the given crack surface load T in the problem of P A , its u a on boundary B u and t a on boundary B t can be obtained  3 Geofluids through analytical solutions [54], which is the same as that in P F . The solutions are expressed as follows: where K u 2 and K t 2 are also the stiffness matrices. Subtract the solution for P A from that one for P F , which results in zero tractions on the crack surface. To ensure the same boundary conditions on B, the following two equations should be satisfied.
Substituting Equations (4) and (5) into Equations (6) and (7) and then substituting the result into Equation (3), the following equation can be obtained: Collecting the terms containing T in one side: This equation can be used to solve for T. Also, substituting Equation (3) into Equations (4) and (5) and then substituting the results into Equations (6) and (7), the following equations can be obtained: Rearrange these two equations: These two equations can be rewritten as follows: where I is the identity operator and Equation (12) can be used to solve for u and t. Similarly, the following linear system can be obtained for the traction t a and displacement u a . where To directly solve these equations involves computation of K u 1 and K t 1 , which requires to obtain traction T for the uncracked domain with prescribed u and t. The problem of P F has to be solved for a larger number of times. So, it is very numerically expensive to directly solve Equation (14). Here, a fixed-point iteration scheme is used for solving this linear system, which is expressed as follows: The solution is as follows: Since the iterative scheme, Equation (16), can be rewritten as follows: Since this fixed-point iteration usually converges very quickly, the iterative method is efficient in solving crackrelated problems for practical engineering.

Fatigue Crack Propagation.
If the crack is subjected to constant amplitude cyclic loading, a plastic zone will form at the crack tip. Usually, the plastic zone is quite small and it is buried in an elastic singularity zone, so the conditions at the crack tip are defined by the value of stress intensity factors. The following equation is often used to express the functional relationship for fatigue crack growth: where ΔK = K max − K min , which is the range of the stress intensity factors, R = K min /K max represents the stress ratio, and da/dN denotes the crack growth during per loading cycle. Figure 3 illustrates the typical fatigue crack propagation behavior for metals. As can be seen from this figure, three distinct regions are contained in this figure. The curve is linear at intermediate ΔK values, but the crack growth rate shows nonlinear trend at both high and low levels. At the low level, da/dN reaches zero at a threshold, which means the crack will not grow below this threshold. On the other hand, the observed fatigue crack growth rate increases very rapidly at high values of ΔK for some materials.

Geofluids
The linear region in Figure 3 can be described by the following power law: where C and n are the material constants which are determined by the experiment. Based on Equation (22), the fatigue crack propagation rate depends only on ΔK, and da/dN is insensitive to the stress ratio R in the linear region. Equation (22) is widely known and called as the Paris fatigue law [55]. Studies over the past several decades have revealed that n ranges from 2 to 4 for most metals if no corrosion exists in the environment.
Based on the Paris fatigue law, the number of loading cycles required to make a crack grow from an initial length a o to the final length a f is determined by the following equation, which also represents the fatigue life of the cracked structure.
Indeed, plasticity-induced crack closure may account for the effects of plate thickness on crack growth rates. Here, the simple Paris law is used to predict the fatigue growth rate. However, other models which account for the effect of plate thickness can also be incorporated in the framework of the current iterative method, which will also be our future study. Moreover, some 3D effects of fracture mechanics are also neglected in this study, such as 3D corner singularity and mode II and III coupling, which are expected to have an insignificant effect for the damage tolerance of the cracked structures.  7 Geofluids elements (Q8), and the crack is meshed using 2-node bar elements (B2). The CPU time for computation will also be mentioned to show the efficiency of the current method in solving crack analysis problems. The Paris fatigue law is employed for the fatigue crack propagation analysis, and the error of stress intensity factor is defined as follows:

Numerical Examples
where K num is the computed SIF and K ana is the analytical or empirical SIF.

Shear Edge Crack.
A finite plate with an edge crack subjected to a shear load is studied in this example. The size of the plate is 7 mm × 16 mm, and the length of the crack is 3.5 mm. The crack tip of the right side is located at the center of the plate. The plate is clamped at the bottom, and a shear load of 1.0 N is applied to its top. The discretization of the plate and the crack is shown in Figure 4, including 800 Q8 elements and 20 B2 elements. After meshing the model in Patran, exporting it into a file and running the computational code, the computed stress intensity factors are listed as follows: The result is K I = 33:67, K II = 1:39 at the crack tip, computed in about 0.25 seconds by a personal computer.

Three Parallel Cracks.
This example studies the problem of three parallel cracks in a finite plate aligned normally to the tensile direction, which is shown in Figure 5. The size of the plate is 40 mm × 80 mm. In Figure 5, the related dimensions are a = 1:27 mm, d = 3:175 mm. The crack B is located in the middle of the plate. An evenly distributed tensile loading of 400 MPa is applied at the edges of the plate. The discretization of the model includes 32 Q8 elements for the plate and 60 B2 elements for the crack.
The computed SIF for K IA is 683.2, as compared to the reference value of 679.6 by Murakami [56], and the computed result for K IB is 604.2, as compared to the reference value of 599.2 by Murakami [56]. The time cost is about 0.19 seconds with a personal computer.

Fatigue Crack Growth of a Slanted Crack.
A finite plate with a slanted crack is studied in this example. The dimension and load are shown in Figure 6(a), with a = 10 mm, b = d = 50 mm, h = 100 mm, θ = 45°, and σ = 50 MPa. Figure 6(b) shows the mesh of the plate and the crack.
Here, fatigue crack growth is simulated by employing the Paris fatigue law, and the parameters for Equation (22) are C = 6:9 × 10 −12 and n = 3. The computed SIFs for the last fatigue step are K I = 39:17, K II = −0:12, and the crack propagation is shown in Figure 7. The time cost is about 1.34 seconds with a personal computer. Figure 7 shows plotting of the values of the crack tip coordinates to illustrate the crack propagation direction, in which the blue line represents initial crack, and the red line represents fatigue propagation of the initial crack. It is important to recognize that as the crack propagates, the values of the stress intensity factors increase, the number of load cycles per step required to produce propagation increment is decreasing since the structure is getting weaker, and that, even for a slanted crack, K I of mode I is dominating over K II of mode II.

Fatigue Growth of an Eccentric Crack in a Cantilever
Beam. In this example, a cantilever beam with the size of 10 mm × 2 mm is studied, as shown in Figure 8. The left side of the cantilever beam is clamped, while two concentrated forces with magnitude of 1000 N are applied to the upper-right and bottom-right corner of the beam, in two opposite directions. The initial length of the crack is 0.5 mm, which is 5% of the beam's length. The crack is off-aligned from the midline of the beam, with an eccentric distance of 0.1 mm, which is also 5% of the beam's height.
The discretization of the model is shown in Figure 8, including 240 Q8 elements for the plate and 15 B2 elements for the crack. The Paris equation is used with C = 6:9 × 10 −15 and n = 3. A maximum increment 1.2 mm is considered for each crack tip, in a total 30 crack growth steps. The predicted final crack shape after fatigue growth is shown in Figure 9, in which the blue line represents initial crack, and the red line represents fatigue propagation of the initial crack.
3.5. Mode I Cracks Growing from a Fastener Hole. In this example, a circular fastener hole with two edge cracks at its left and right side is considered. As shown in Figure 10(a), size of the plate is 50:8 mm × 101:6 mm. A circular hole is located at the center of the plate with radius R = 6:35 mm, and two cracks emanate from the hole. The initial length of each crack is 2.5 mm; thus, the summation of the radius and the crack length is a = 8:85 mm. An evenly distributed tensile loading of 82.7 MPa is applied to the edges of the plate.
The discretization of the model is shown in Figure 10(b), including 657 Q8 elements for the plate and 60 B2 elements for the crack. The Paris equation is used with C = 6:9 × 10 −15 and n = 3. A maximum increment 10 mm is considered for each crack tip, in a total 20 crack growth steps. The predicted final crack shape after fatigue growth of 4939 loading cycles is shown in Figure 11.
3.6. Mixed Mode Fatigue Growth of Multiple Slanted Cracks. Figure 12 shows the geometry of a plate with two holes and emanating cracks. The geometric parameters are as follows: H = 40 mm, w 1 = w 3 = 15 mm, and w 2 = 10 mm. The radius of the holes is R = 0:805 mm, the initial cracks' lengths are a 1 = a 2 = a 3 = a 4 = 0:3 mm, and the initial cracks' angles are θ = 30°. The fatigue-related parameters are C = 3 × 10 −13 and n = 3. Fatigue uniform stress loading is applied on the upper and lower edges, and the maximum stress is 12000 Pa and the stress ratio r = 0:1. Figure 13 shows the mesh of the plate surface. The loads are specified at the upper and lower edges of the plate, and the boundary conditions are defined properly as also shown in Figure 13.

Geofluids
The propagation of the cracks is shown in Figure 14. Total number of 80272 cycles is applied. The stress intensity factors are shown in Figure 15. The stress intensity factor of mode I, K I , is increasing while that of K II is constant. So, it is clear that crack mode I is dominating over mode II and hence the propagation goes horizontally as shown in the figures.

Conclusions
In this paper, based on superposition principle, the iterative method was developed and used to obtain the stress intensity factors and simulate fatigue crack propagation of various kinds of structures with different initial cracks. Several typical numerical examples were illustrated in the paper. The following conclusions can be drawn by carefully comparing the numerical results with analytical solutions as well as empirical observations available in the open literature.
The iterative method requires independent and very coarse meshes for both the uncracked global structure and the crack; it only requires very little computational burden and human labor cost for modeling the fatigue propagation of various kinds of cracks. The computed stress intensity factors for cracks using this method are in good agreement with analytical solutions or empirical solutions. The whole crack growth path up to failure of the structures can be easily and efficiently simulated using this method.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.