Direct Evaluation of the Stress Intensity Factors for the Single and Multiple Crack Problems Using the P-Version Finite Element Method and Contour Integral Method

: Stress intensity factor (SIF) is one of three important parameters in classical linear elastic fracture mechanics (LEFM). The evaluation of SIFs is of great signiﬁcance in the ﬁeld of engineering structural and material damage assessment, such as aerospace engineering and automobile industry, etc. In this paper, the SIFs of a central straight crack plate, a slanted single-edge cracked plate under end shearing, the offset double-edge cracks rectangular plate, a branched crack in an inﬁnite plate and a cruciﬁx crack in a square plate under bi-axial tension are extracted by using the p-version ﬁnite element method (P-FEM) and contour integral method (CIM). The above single-and multiple-crack problems were investigated, numerical results were compared and analyzed with results using other numerical methods in the literature such as the numerical manifold method (NMM), improved approach using the ﬁnite element method, particular weight function method and exponential matrix method (EMM). The effectiveness and accuracy of the present method are veriﬁed.


Introduction
In the field of fracture mechanics, the stress field near the crack tip is usually described by the SIFs. For brittle materials, it can also be used to describe the non-singular stresses on the crack surface. The SIF is of great significance for the research of material fracture and crack propagation in practical engineering problems. At present, there have been many numerical methods which were applied to solve the SIFs, such as the traditional finite element method (FEM) and boundary element method (BEM), extended finite element method (XFEM), meshless method (MLM), EMM, the improvement of the FEM, the particular weight function method and NMM, etc. After many efforts of numerous researchers, many methods were also used to extract the SIFs, such as the extrapolation method, interaction integral method and J integral method, etc.
In the past twenty years, the XFEM has mostly been used to analyze fracture mechanics problems. The XFEM is based on the partition of the unity finite element method, adding the step function and the asymptotic field of crack tip displacement to describe the crack, and has achieved fruitful results in the analysis and research of discontinuous problems such as crack and crack propagation [1][2][3][4][5]. The EMM [6] converts the linear, elastic, two-dimensional and steady-state equations of the crack problem into a Hamiltonian system. By changing the angular stress to the radial stress, the eigensolutions satisfying the adjoint symplectic orthogonality are obtained, and then the desired solution of the problem is obtained by using the linear combination of these eigensolutions. The improved

Development of P-FEM
In past three decades, the P-FEM has been developed rapidly and its mathematical theory has been established completely. The convergence of numerical solutions of P-FEM has also been strictly proved, and its error estimation has also been obtained [12][13][14]. Compared with the traditional FEM, the P-FEM has the advantages of less preprocess, faster convergence rate and higher accuracy.
At present, the P-FEM has been widely used in various engineering practice fields, and has obtained fruitful research results. The P-FEM has been applied to the unsteady temperature field problem [15], beam vibration problem [16][17][18] and micro rectangular plate thermal buckling analysis [19]. In the three-dimensional elastic viscoelastic problems in the hydraulic structure [20], seepage analysis [21], the mechanical response of the arteries in the field of biomechanics of bone problems [22] and sandwich and Kirchhoff plate deformation analysis [23,24], the P-FEM has been also successfully applied. The application of P-FEM in elastic fracture mechanics (including the calculation of SIF and predicting the growth of irregular crack fronts) can be found in references [25][26][27][28][29].

Shape Functions of the P-FEM
For two-dimensional problems, there are usually two types of shape functions (interpolation shape function and hierarchic shape function) defined in standard quadrilateral element S = (−1, 1) × (−1, 1) and triangular element T: Taking a two-dimensional standard triangle element (as shown in Figure 1) as an example, the high-order hierarchic shape function of the P-FEM on the standard triangle element is constructed as follows [30]:

Nodal Modes Shape Functions for p =1
The nodal modes shape function is the shape function related to nodes. The nodal modes shape function for p = 1 is the same as the conventional Lagrangian interpolation function. Taking the standard triangle element as an example, the nodal modes shape function are as follows: 2.2.2. Side Modes Shape Functions for p ≥ 2 The side modes shape functions are the shape functions related to the edges. For example, the side modes shape function of the edge N i (k) is the shape function related to the corresponding side Γ k . The side modes shape function associated with the edge Γ 1 = {η = 0, −1 ≤ ξ ≤ 1} is defined as follows: where P n (t) is Legendre polynomial of n ≥ 0: Similarly, the side modes shape functions associated with the side Γ k (2 ≤ k ≤ 3) are defined as: N When p = 4, the new internal modes shape function is added: Similarly, when p = 8, the new internal modes shape function is added: where P i is the Lagrange interpolation polynomial, and The number of hierarchic shape functions is determined by the degree p. In this paper, the triangular element and high-order hierarchic shape functions are used to calculate the stresses and strains near the crack tip, and the computation is carried out under the condition of p = 1, 2, . . . , 8.

Contour Integral Method
The CIM is a super-convergent numerical method based on the Betti's theorem of reciprocal works. In recent years, the CIM has been widely used in various fields to solve various practical or theoretical problems, and has yielded many excellent results. For example, Pereira et al. [31] analyzed the loading crack based on the CIM and the generalized finite element method. Leblanc et al. [32] solved the acoustic nonlinear eigenvalue problem based on the CIM. Garzon et al. [33] studied the three-dimensional crack problems based on the CIM and truncation function method. Feng et al. [34] analyzed the reflection radiation error between turbine blades and blades in computational fluid dynamics based on the CIM. Cui et al. [9] calculated the SIFs based on the CIM and displacement discontinuity method.
In this paper, the CIM proposed by Szabó et al. [35] and Pereira et al. [31] is adopted to extract the SIFs of the single-and multiple-crack problems. The detailed discussions about the CIM were reported in [31,35]. According to Pereira et al. [31], the SIFs can be obtained by calculating the follow integrals: are traction vectors respectively calculated by v −I i and v −II i , which are the export functions respectively corresponding to modes I and II SIFs, and their definitions are as follows: The SIFs K I and K II calculated by Equations (11) and (12) are exact. Definitions of symbols in Equations (11)- (14) refer to [31]. Due to the reason that displacement field u and the corresponding traction vector T (u) are unknown, based on the P-FEM, the approximate solutionū of u and the approximate solution of T (u) is obtained, and the numerical solutions of the SIFs K I and K II are finally extracted. The selection of contour Γ 2 in the computation is flexible and important. Generally, we choose the contour Γ 2 surrounding the innermost elements near the crack tip.

A Central Straight Crack Plate
In order to demonstrate the validation of the present approach, the first benchmark example is a central straight crack plate subjected to a uniform uniaxial tension σ = 3 MPa on the upper and bottom edges shown in Figure 2. For this example, the length of the plate is 2L (L = 2 m), the width is 2W (W = 1 m), the ratio of crack length to width a/W = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 and 0.7 are considered, and the plane stress condition is assumed. The material parameters are elastic modulus E = 1000 MPa, Poisson's ratio ν = 0.3. The analytical solution of the SIF K I for this model is given in [36] as follows: For this simulation, 252 elements are used for discretization (as shown in Figure 2b). The computation results of the SIF K I are summarized in Tables 1 and 2 Table 1; here we only list the results for the case a/W = 0.4, the results for other cases are similar. It can be seen from Table 2, compared with the high-order NMM [37], that an obvious improvement in accuracy is observed.      The stress nephogram of the example is shown in Figure 3. For the sake of simplicity, only the stress nephogram for the first example is plotted, the others are similar.

A Slanted Single-Edge Crack Plate
The rectangular plate under end shearing with oblique crack is shown in Figure 4. The bottom of the model is fixed, and the top is imposed to mixed-mode shear load τ = 1MPa, where W = L = 10m, the initial crack length a = 6m, and the initial crack direction β =30 • . The material parameters of the model are elastic modulus, E = 210 GPa, Poisson's ratio ν = 0.25. The model is computed under plane stress condition. The present computation results are shown in Table 3, compared with results in the literature [6,38], which are in good agreement and verify effectiveness of the present method. The results K I = 20.485 MPa·m 1/2 , K II = 9.8188 MPa·m 1/2 , T = 5.7802 MPa in [6] are referred as the reference solutions, the present computation results are K I = 20.490 MPa·m 1/2 , K II = 9.821 MPa·m 1/2 , T = 5.781 MPa. The relative errors of K I , K II , and T are 0.024%, 0.022% and 0.0138% respectively. The errors of energy norm of the P-FEM solutions and relative errors of normalized SIFs comparison with Abaqus solutions [6] for the example are listed in Table 4.

The Offset Double Edge Cracks Plate
The rectangular plate with the offset double-edge cracks is shown in Figure 5. The top and bottom of the model are subjected to uniformly distributed tensile load σ = 1 MPa, the geometric size of the model is W × L, where W = 4 m, and the initial crack length on both sides is a = 2 m. The material parameters of the model plate are elastic modulus E = 1 MPa, Poisson's ratio ν = 0.25, and the plane stress condition is assumed.
The SIFs under various crack spacing conditions are computed, and the influence of the ratio of length to width of different plates is also considered, as shown in Tables 5-7. The present results are in good agreement with those in the literature [7,8], as shown in Figures 6-8.      The errors of energy norm of the P-FEM solutions and relative errors of SIFs comparison with Chen et al. [8] for the case L/W = 3, c/L = 0.3 are listed in Table 8, the others are similar. Table 8. Errors of energy norm of the P-FEM solutions and relative errors of SIFs (L/W = 3, c/L = 0.3) comparison with Chen et al. [8].
p DOF Error of Energy Norm/(%)

A Branched Crack in an Infinite Plate
The example is a branched crack in an infinite plate, as shown in Figure 9. A branched cracks embedded in an infinite plate are simulated by the p-version FEM and CIM to demonstrate the accuracy and effectiveness when dealing with complex cracks. Uniform tension σ = 1 MPa unit is applied to the boundary of the plate. The dimensions of the model are W = 5 m, L = 4 m, a/W = 0.05, b/a = 0.9 and θ = 45 • , and the plane stress condition is assumed. The material parameters are considered as Young's modulus E = 1000 MPa and Poisson's ratio ν = 0.3. The normalized SIFs for crack tips A and C are as follows: In this example, the normalized SIFs reported in the literature [39] are used as reference solutions. The numerical results obtained by the present method are shown in Table 9, which are compared with the results using the NMM in the literature [40,41] and the reference solutions. It is easy to see that relative error is no more than 0.21% for all the cases, which is less than the relative error in Cai [40,41].

A Crucifix Crack in a Square Plate under Bi-Axial Tension
Last example is a crucifix crack in a square plate under bi-axial tension, as shown in Figure 10. The parameters for this example are taken the same as [42,43]: W = 0.5 m, σ = 1 MPa, E = 1.0 × 10 3 MPa, ν = 0.3, and the plane strain condition is assumed. Due to symmetry, the SIF at each crack tip is theoretically the same; therefore, only the results at crack tip A are considered. The normalized SIFs F I = K I /σ √ πa and present numerical results comparison with references [42,43] are shown in Table 10.   Table 11. In the computations of the first model, the total number of discretizing meshes is 252, and the maximum DOF is 15,678; in the second to the fifth models, the total number of discretizing meshes is 39, 142, 329 and 580, respectively, and the maximum DOF is 2524, 9064, 20,722 and 36,206, respectively. The above examples show that coarse meshes are enough to obtain quite accurate SIFs using the present method. It can be seen that the convergences are very fast with the increasing of degree p from 1 to 8, the accuracy of the SIFs obtained by the present method is also satisfactory, and the SIFs have better stability.

Conclusions
The SIFs of several classical LEFM models were evaluated based on the P-FEM and CIM, and compared with the results in the literature to verify the effectiveness and accuracy of the present method. The following conclusions can be drawn:

•
The computation based on the P-FEM and CIM only needs a fewer number of meshes (less preprocessing), has a faster convergence rate, which can not only be used to analyze the single-crack problems, but also can be used to calculate the SIFs of multiplecrack problems. Furthermore, the present method can also be used to analyze threedimensional fracture problems and simulate the crack growth path.

•
The CIM is a super-convergent method based on the Betti's theorem of reciprocal works. The SIFs are derived by using the displacements and stresses of higher accuracy on the integral path. Therefore, the CIM can obtain results of relatively higher precision.

•
The P-FEM does not need fine meshes, has a faster convergence rate and higher accuracy. However, in the study of crack propagation behavior, it is inevitable to meet the problem of remeshing. On the other hand, the XFEM need not consider the internal geometry of structures and materials; it employs the partition of unity (PU) concept to simulate fracture and crack propagation. The XFEM uses the enrichment functions to reflect the local characteristics of the crack tip and avoids remeshing of meshes. The high-order extended finite element method (high-order XFEM), a combination of the P-FEM with XFEM, will probably have the advantages of both: avoiding remeshing as well as obtaining higher accuracy. The study of the high-order XFEM will be an interesting issue.  Data Availability Statement: Data is contained within the present article. Other data presented in this research are available in [6][7][8][37][38][39][40][41][42][43].