Direct Computation of 3-D Stress Intensity Factors of Straight and Curved Planar Cracks with the P-Version Finite Element Method and Contour Integral Method

This paper presents direct computations of 3-D fracture parameters including stress intensity factors (SIFs) and T-stress for straight and curved planar cracks with the p-version finite element method (P-FEM) and contour integral method (CIM). No excessive singular element or enrichment function is required for the computation. To demonstrate the accuracy and efficiency of the proposed approaches, several benchmark numerical models of 3-D planar straight and curved cracks with single and mixed-mode fractures are considered and analyzed: a through thickness edge straight crack in a homogeneous material, a through thickness inclined straight crack, a penny-shaped crack embedded in a cube and a central ellipse shaped crack in a homogeneous cube. Numerical results are analyzed and compared with analytical solutions and those reported by the extended finite element method (XFEM) and the scaled boundary finite element method (SBFEM) in the available literature. Numerical experiments show the accuracy, robustness and effectiveness of the present method.


Introduction
Inherent flaws or cracks in complex engineering materials and structures are unavoidable. The existence of cracks will affect the safety and durability of structures. The analysis of 3-D crack behavior is of great significance in practical engineering. The stress intensity factors are important fracture parameters in linear elastic fracture analysis, which are often applied to characterize the displacements and stresses fields near the crack tip and further to predict crack propagations. In an analysis of three-dimensional fracture problems, the SIFs K I , K II and K III accurately describe the stress state at the crack tip. The accurate computation of 3-D fracture parameters (including SIFs and T-stress) for straight or curved planar cracks with single and mixed-mode loading is very important and has always been a challenging issue in linear elastic fracture analysis.
Classical problems in 3-D fracture mechanics include: a through thickness edge straight crack in a homogeneous specimen, a through thickness inclined crack, a pennyshaped crack in a homogeneous rectangular plate and an ellipse shaped crack embedded in a homogeneous large cube. A lot of techniques are proposed to analyze the aforementioned fracture problems. Numerical results of the SIFs for a through thickness edge straight crack in a homogeneous specimen have been reported by using different numerical methods in [1,13,14]. Wang et al. [13] evaluated the three-dimensional SIFs of an edge parallel crack specimen in a rectangular block using the adaptive multiscale XFEM and interaction integral method. Wang et al. [14] computed the SIFs of an edge straight crack in a cube and a central ellipse shaped crack in a large cube using local mesh refinement in combination with the XFEM and greatly improved the accuracy of SIFs computed by the classical extended finite element method. Wang et al. [14] also studied a mixed-mode edge inclined straight crack in a plate and considered the effects of various crack lengths for the SIFs. Saputra et al. [15] computed the SIFs and T-stress of a penny-shaped crack in a homogeneous cube using the scaled boundary finite element method (SBFEM).
Bhat et al. [16] analyzed the effect of specimen parameters on mixed-mode I/II SIFs for additive manufactured slant edge crack plate. The mixed-mode I/II stress intensity factors were studied for an inclined crack in a specimen. The authors observed that the specimen parameters, the crack length ratio, crack angle and layer-built angle significantly affect the stress intensity factors of slant edge crack plate. Khosravani et al. [17] experimentally investigated fracture characteristics and load-carrying capacity of 3-D-printed components (polycarbonate and nylon filaments) produced using fused deposition modeling, and their numerical simulation results are in very good agreement with reported experimental results. The above studies are important for the domain of manufacturing and engineering. We aim to generalize the present method to analyze similar industrial applications problems in our next manuscript.
The objective of this paper is to propose an effective numerical approach to accurately estimate the stress intensity factors (SIFs) and T-stress for 3-D planar straight and curved cracks. Especially, for the accurate computation of 3-D fracture parameters, which is a challenging task in 3-D cracked structures, the numerical approach is able to improve the accuracy of fracture parameters compared to other numerical methods. In this paper, the fracture parameters, including the stress intensity factors (SIFs) and T-stress for 3-D planar straight and curved cracks with single and mixed-mode of several benchmark problems, were computed directly using P-FEM in combination with CIM. The validity and accuracy of the present approaches are illustrated by numerical models, and numerical results were analyzed and compared with analytical solutions and published results obtained by the adaptive multiscale extended finite element method (XFEM), local mesh refinement XFEM (Lm-XFEM) and scaled boundary finite element method (SBFEM) [13][14][15]. Computation examples show that the SIFs computed by the present method have faster convergence speed, smaller relative error and better stability.

P-Version Finite Element Method in Three Dimensions
Approximation theory of the three-dimensional p-version of the finite element method has been established [18][19][20]. The p-version of the finite element method has turned out to be an efficient discretization technique for many engineering problems. The p-version FEM has an exponential rate of convergence in energy norm with a proper mesh design. Generally, the p-version FEM is superior to the traditional h-version FEM. The application of the p-version FEM in three dimensions, including in three -dimensional linear elastic fracture analysis, can be found in [21][22][23][24][25][26][27].

3-D Hierarchical Shape Functions for Standard Tetrahedral Elements
There are several element types for three-dimensional problems: the hexahedral, tetrahedral, pentahedral (wedge), and pyramid elements. The most common elements are the hexahedral and tetrahedral elements (shown in Figure 1).

3-D Hierarchical Shape Functions for Standard Tetrahedral Elements
There are several element types for three-dimensional problems: the hexahedral, tetrahedral, pentahedral (wedge), and pyramid elements. The most common elements are the hexahedral and tetrahedral elements (shown in Figure 1). There are four types of hierarchical shape functions of the three-dimensional p-version FEM on standard tetrahedral element T: (1) nodal modes shape functions, (2) edge modes shape functions, (3) face modes shape functions and (4) internal modes shape functions. Details about these shape functions defined on a standard tetrahedral element can be found in [28]. In this paper, we use tetrahedral elements to compute fracture parameters, and we find the results obtained from tetrahedral elements are better than those from hexahedral elements.

Contour Integral Method in Three Dimensions
Szabó et al. [29] first reported the contour integral method as a super-convergent technique to calculate stress intensity factors. The contour integral method can be utilized to calculate modes I, II and III SIFs. Garzon et al. [9] applied CIM in combination with the generalized finite element method to extract the SIFs of 3-D crack. The three-dimensional formulae used in [9] are as follows: where  [9]. The stress intensity factors K I , K II and K III in (1), (2) and (3) are exact solutions, but u and the corresponding traction vector T (u) are unknown. Any numerical method can be used to approximate the exact solutions u. In this study, numerical solutions calculated by the three-dimensional p-version finite element method are used to replace the exact solutions u, the corresponding traction vectors calculated are used to replace the exact traction T (u) and the approximate solutions of K I , K II and K III are finally obtained. There are four types of hierarchical shape functions of the three-dimensional p-version FEM on standard tetrahedral element T: (1) nodal modes shape functions, (2) edge modes shape functions, (3) face modes shape functions and (4) internal modes shape functions. Details about these shape functions defined on a standard tetrahedral element can be found in [28]. In this paper, we use tetrahedral elements to compute fracture parameters, and we find the results obtained from tetrahedral elements are better than those from hexahedral elements.

Contour Integral Method in Three Dimensions
Szabó et al. [29] first reported the contour integral method as a super-convergent technique to calculate stress intensity factors. The contour integral method can be utilized to calculate modes I, II and III SIFs. Garzon et al. [9] applied CIM in combination with the generalized finite element method to extract the SIFs of 3-D crack. The three-dimensional formulae used in [9] are as follows: where . A detailed analysis of the 3-D CIM can be found in [9]. The stress intensity factors K I , K II and K III in (1), (2) and (3) are exact solutions, but u and the corresponding traction vector T (u) are unknown. Any numerical method can be used to approximate the exact solutions u. In this study, numerical solutions calculated by the three-dimensional p-version finite element method are used to replace the exact solutions u, the corresponding traction vectors calculated are used to replace the exact traction T (u) and the approximate solutions of K I , K II and K III are finally obtained.

Numerical Examples and Discussion
In order to illustrate the effectiveness and the accuracy of the present approaches, the SIFs and T-stress of three-dimensional planar straight and curved cracks in 3-D objects were computed for the following four crack configurations: (1) a through thickness edge straight crack in a homogeneous material; (2) a mixed-mode edge inclined straight crack Based on the 3-D p-version FEM and CIM, numerical models are presented in two kinds of problems: planar straight cracks and planar curved cracks in 3-D objects with single and mixed-mode fractures. The first two models are used to investigate single and mixed-mode fractures in 3-D objects by studying an edge straight crack and an edge inclined straight crack and to verify the validity and accuracy of the present method by comparing numerical results with analytical solutions [30] and those obtained by the adaptive multiscale extended finite element method [13], local mesh refinement extended finite element method [14] and scaled boundary finite element method [15]. In the last two examples, the 3-D p-version finite element method and contour integral method are applied to investigate planar curved cracks in 3-D by studying a central penny-shaped crack and a central elliptical crack and to demonstrate the accuracy and effectiveness by comparing numerical results with reference solutions [31,32] and those obtained by the scaled boundary finite element method [15] and the local mesh refinement extended finite element method [14]. All examples can be compared with those published in the literature to evaluate the accuracy of the present method; these examples show the simplicity in computation of the SIFs without the need to use excessive singular elements or enrichment functions near the crack front.

An Edge Straight Crack
In this subsection, an example of an edge straight crack in a rectangular block is presented to demonstrate the efficiency and accuracy of the present approach. The same example has been investigated by Wang et al. [13] and Saputra et al. [15]. The SIFs K I are computed by the p-version FEM in combination with CIM and compared with the analytical solutions and results reported by the adaptive multiscale extended finite element method and the scaled boundary finite element method. This example is considered under the plane strain condition.

SIFs of an Edge Straight Crack in a Rectangular Block
Consider a rectangular specimen of size 10 m × 10 m × 20 m which contains an edge straight crack as shown in Figure 2a. The crack length is 5 m; the top surface and bottom surface of the specimen are subjected to a unidirectional uniform traction of σ = 1 Pa. The material parameters are considered for E = 206 GPa and ν = 0.3 (Young's modulus and Poisson's ratio). The analytical solution of the SIFs K I for the edge straight crack problem is given by [30] under the plane strain condition and d l ≤ 0.6: where d is the crack length, l is the length of the rectangular block along the crack direction The SIFs K I computed by the present approach and the analytical solutions are listed in Table 1. A comparison with the results obtained by the adaptive multiscale extended finite element method in [13] are shown in Figure 3. In this example, the normalized analytical SIFs K I = 2.826. Obviously, from Figure 3, the present SIFs K I are closer to the analytical solutions at most points of the crack front.  The SIFs K I computed by the present approach and the analytical solutions are listed in Table 1. A comparison with the results obtained by the adaptive multiscale extended finite element method in [13] are shown in Figure 3. In this example, the normalized analytical SIFs K I = 2.826. Obviously, from Figure 3, the present SIFs K I are closer to the analytical solutions at most points of the crack front.
A similar model has been considered by Saputra et al. [15]. At the same conditions, namely, when elastic material properties are the same as Young's modulus E and Poisson's ratio ν = 0.3, when uniform uniaxial tension is subjected at the top and bottom surfaces of the specimen and when the other boundaries are traction free, the SIFs and T-stress are to be calculated by the present method. The meshes for the present method are shown in Figure 2b,c. Comparisons with the results obtained by the scaled boundary finite element method in [15] are shown in  A similar model has been considered by Saputra et al. [15]. At the same conditions, namely, when elastic material properties are the same as Young's modulus E and Poisson's ratio 0.3 ν = , when uniform uniaxial tension is subjected at the top and bottom surfaces of the specimen and when the other boundaries are traction free, the SIFs and T-stress are to be calculated by the present method. The meshes for the present method are shown in Figure 2b,c. Comparisons with the results obtained by the scaled boundary finite element method in [15] are shown in Figure 4.     In the computation of the example, the number of discretizing meshes is 18,447, and the degrees of hierarchical shape functions are p = 1, · · · , 8. The results of K I along the crack front are presented in Figure 4. The present results are in very good agreement with the results of the SBFEM obtained with element order 9 and are closer to the analytical solution.  Figure 5b,c shows the P-FEM mesh of the edge inclined straight crack and the present local mesh refinements in the crack front.

An Edge Inclined Straight Crack
Normalized SIFs K I and K II are computed for different crack lengths ( a varies from 3 to 5) using the present method and are listed in Tables 2 and 3. The reference solutions are given by the Stress Intensity Factors Handbook [30]. Comparisons with the reference solution and results of the normalized SIFs K I and K II in [14] are shown in Figures 6 and 7. (a) (b) (c)    Normalized SIFs K I and K II are computed for different crack lengths (a varies from 3 to 5) using the present method and are listed in Tables 2 and 3. The reference solutions are given by the Stress Intensity Factors Handbook [30]. Comparisons with the reference solution and results of the normalized SIFs K I and K II in [14] are shown in Figures 6 and 7.   In the computation of the example, the number of discretizing meshes is 33,292, and the degrees of hierarchical shape functions are p from 1 to 8. It is easy to see from Figures 6 and 7 that the present results of normalized SIFs K I and K II are closer the reference solutions.    In the computation of the example, the number of discretizing meshes is 33,292, and the degrees of hierarchical shape functions are p from 1 to 8. It is easy to see from Figures 6 and 7 that the present results of normalized SIFs K I and K II are closer the reference solutions. In the computation of the example, the number of discretizing meshes is 33,292, and the degrees of hierarchical shape functions are p from 1 to 8. It is easy to see from Figures 6 and 7 that the present results of normalized SIFs K I and K II are closer the reference solutions.

A Penny-Shaped Crack
Accurate computation of 3-D fracture parameters (such as SIFs and T-stress) for curved cracks is a challenging problem in linear elastic fracture analysis. Many efforts have been made to improve the accuracy of 3-D fracture parameters in the literature in the past two decades. A central penny-shaped crack with radius a in a homogeneous material is investigated as shown in Figure 8a. The dimensions of the model are a/w = 0.08, w = 80. The material parameters are used for E and ν = 0.3 (Young's modulus and Poisson's ratio). Accurate computation of 3-D fracture parameters (such as SIFs and T-stress) for curved cracks is a challenging problem in linear elastic fracture analysis. Many efforts have been made to improve the accuracy of 3-D fracture parameters in the literature in the past two decades. A central penny-shaped crack with radius a in a homogeneous material is investigated as shown in Figure 8a. The dimensions of the model are  Figure 8b-d shows the present p-version FEM mesh of the penny-shaped crack and the local mesh discretization in the crack front. The analytical solutions of the stress intensity factor K I and the T-stress T in the penny-shaped crack in an unbounded domain are given in [31,32]: The analytical solutions (4) are used as the reference solutions in the model. The SIFs K I and T-stress T of the penny-shaped crack calculated by the present method, the scaled boundary finite element method [15] and analytical solutions are compared in Tables 4  and 5   A uniaxial uniform tension σ = 1 Pa is subjected to the top and bottom surfaces of the cube. In this example, the p-version FEM and CIM are applied to extract 3-D SIFs and T-stress for the planar curved cracks in three dimensions. Figure 8b-d shows the present p-version FEM mesh of the penny-shaped crack and the local mesh discretization in the crack front. The analytical solutions of the stress intensity factor K I and the T-stress T in the penny-shaped crack in an unbounded domain are given in [31,32]: The analytical solutions (4) are used as the reference solutions in the model. The SIFs K I and T-stress T of the penny-shaped crack calculated by the present method, the scaled boundary finite element method [15] and analytical solutions are compared in Tables 4 and 5 separately and shown in Figures 9 and 10.
In Tables 4 and 5, as the computation values of SIFs K I and T-stress T are not given in [15], only the differences of the computed values K I and T to the exact values are listed from [15]. In the computation of the example, the number of discretizing meshes is 33,478, and the degrees of hierarchical shape functions are p from 1 to 8. It is easy to see from Tables 4 and 5 that the relative errors of the SIFs K I computed by the present method are less than and equal to 2%, and the relative errors of the T-stress T computed by the present method are less than 2.15%. All of them are less than the results of the SBFEM.    In Tables 4 and 5, as the computation values of SIFs K I and T-stress T are not given in [15], only the differences of the computed values K I and T to the exact values are listed from [15]. In the computation of the example, the number of discretizing meshes is 33,478, and the degrees of hierarchical shape functions are p from 1 to 8. It is easy to see from Tables 4 and 5 that the relative errors of the SIFs K I computed by the present method are less than and equal to 2%, and the relative errors of the T-stress T computed by the present method are less than 2.15%. All of them are less than the results of the SBFEM. Table 6 shows that the errors of energy norm and relative errors of SIF K I for interpolation polynomial order p from 1 to 8. It is easy to see that the error of energy norm and relative error of SIF K I gradually decrease with the increasing of the polynomial degree p, and the errors gradually stabilize. Other examples are similar to this one; for simplicity, we list only Table 6 for reference.    Table 6 shows that the errors of energy norm and relative errors of SIF K I for interpolation polynomial order p from 1 to 8. It is easy to see that the error of energy norm and relative error of SIF K I gradually decrease with the increasing of the polynomial degree p, and the errors gradually stabilize. Other examples are similar to this one; for simplicity, we list only Table 6 for reference.

A Central Ellipse Shaped Crack
where θ is the elliptic angle and ( ) E k is the second kind of elliptic integral, The model has been analyzed by Wang et al. [14] using the local mesh refinement extended finite element method, and excellent results have been obtained in [14]. The SIFs K I computed by the present method and analytical solutions are given in Table 7, and a comparison with the analytical solutions and results of the SIFs K I extracted by the local mesh refinement extended finite element method in [14] are shown in Figure 12.  The exact SIFs for a central ellipse shaped planar crack in an infinite domain is listed in [32]: where θ is the elliptic angle and E(k) is the second kind of elliptic integral, The model has been analyzed by Wang et al. [14] using the local mesh refinement extended finite element method, and excellent results have been obtained in [14]. The SIFs K I computed by the present method and analytical solutions are given in Table 7, and a comparison with the analytical solutions and results of the SIFs K I extracted by the local mesh refinement extended finite element method in [14] are shown in Figure 12. Figure 12 shows clearly that the SIFs K I extracted by the present method approximate the analytical solutions and are in good agreement with those obtained by the local mesh refinement extended finite element method. In the computation of the example, the number of discretizing meshes is 36,247, and the degrees of hierarchical shape functions are p from 1 to 8. It is easy to see that the relative errors of the SIFs K I when 30 ≤ ≤ 60 are bigger than those when 0 ≤ ≤ 30 and 60 ≤ ≤ 90; the curve is more stable and there is no large jump.    Figure 12 shows clearly that the SIFs K I extracted by the present method approximate the analytical solutions and are in good agreement with those obtained by the local mesh refinement extended finite element method. In the computation of the example, the number of discretizing meshes is 36,247, and the degrees of hierarchical shape functions are p from 1 to 8. It is easy to see that the relative errors of the SIFs K I when 30

Conclusions
The 3-D p-version FEM and contour integral method are applied in order to analyze 3-D fracture problems. The hierarchical high-order shape functions are used to approximate 3-D fracture parameters. No excessive singular element or enrichment function near the crack front is required for the computation. Numerical experiments show that the p-version FEM and CIM are efficient and accurate for the direct computation of stress intensity factors and T-stress in three dimensions. It is effective not only for straight planar cracks but also for curved planar cracks in the analysis of three-dimensional fracture problems. Numerical examples show that the fracture parameters determined by

Conclusions
The 3-D p-version FEM and contour integral method are applied in order to analyze 3-D fracture problems. The hierarchical high-order shape functions are used to approximate 3-D fracture parameters. No excessive singular element or enrichment function near the crack front is required for the computation. Numerical experiments show that the p-version FEM and CIM are efficient and accurate for the direct computation of stress intensity factors and T-stress in three dimensions. It is effective not only for straight planar cracks but also for curved planar cracks in the analysis of three-dimensional fracture problems. Numerical examples show that the fracture parameters determined by the present method are in very good agreement with the reported results in the literature and have higher precision and better numerical stability in most cases.
Although the present method is effective and straightforward for the computation of 3-D fracture parameters, it inevitably needs to remesh meshes near the crack front. Local mesh discretization or refinement is required in order to obtain results with higher accuracy. To overcome the drawback, it would be interesting if the present method could be combined with the extended finite element method to avoid remeshing. Data Availability Statement: Data is contained within the present article. Other data presented in this research are available in [13][14][15] and [30][31][32].