Fracture analysis for materials by a stable generalized/extended finite element method

In this paper, a method is proposed for extracting fracture parameters in isotropic material cracking via a stable generalized/extended finite ele-mentmethod.Thenumericalresultsofthestressintensityfactorsandscaledconditionnumberofthesystemmatrixarepresentedandcompared with different enrichment schemes or those reported in related references. The good agreement and convergence of the results obtained by the developedmethodwiththoseobtainedbyothersolutionsorenrichmentschemesprovestheapplicabilityoftheproposedapproachandconfirmsitscapabilityofefficientlyextractingfractureparametersinisotropicmaterials.


INTRODUCTION
In recent years, the fracture of some advanced materials has an increasing research value with an in-depth investigation of various materials. There are many fracture analysis methods for materials, such as numerical methods. Among them, the finite element method (FEM) is the most widely used. However, there are some shortcomings in the FEM. For example, the mesh division is required to be consistent with the crack surface, and the remeshing is required during the crack propagation process. Therefore, the generalized/extended finite element method (G/XFEM) was proposed in 1999, inheriting all the advantages of the traditional FEM and possessing the high efficiency of fracture problem analysis. Particularly, its calculation grid does not need to be consistent with the crack surface, it can effectively analyze the discontinuity problem, and it does not need to conduct high-density grid division near the crack tip, contributing to greatly saving the calculation cost. Meanwhile, G/XFEM also has some disadvantages. For instance, it is difficult to deal with the blending element; the condition number of the corresponding system matrix is generally high, leading to ill-conditioned equations. This is due to linear dependence between the finite element functions and enrichment functions. Therefore, a stable generalized/extended finite element method (SG/XFEM) is developed to overcome the problem of the high condition number of the system matrix, and it can overcome the problem that the condition number of the system matrix is too high.
At present, there are various numerical methods to study material fractures, such as FEM, G/XFEM, boundary element method and explicit analytical method. Most of the analytical methods can only analyze the fracture problems of cracked com-ponents with simple geometry or simple loading conditions. For example, the numerical method is mainly used to analyze the material fracture because it is not suitable to solve the fracture problems of cracked components with complex boundaries or bearing complex loading conditions. Since 1990, the FEM has become one of the most efficient methods to investigate material fracture. It can be easily used in discontinuous fracture problems with arbitrary geometric boundary conditions. Chandra Kishen and Singh [1] employed the FEM to analyze the bending and bifurcating crack problems in concrete materials and obtained the stress intensity factors (SIFs) of fracture parameters. However, FEM requires the crack surface to be consistent with the element boundary when solving the crack problems; besides, it requires constant grid reconstruction with the crack growth. Thus, it is not suitable for analyzing the high stress concentration and three-dimensional crack problems. Moreover, there are also shortcomings in the analysis of some crack problems.
In 1999, Belytschko and Moës [2,3], American scholars, proposed G/XFEM, which is based on FEM and partition of unity finite element method [4] and inherits all the characteristics of FEM. Recently, it has become the most commonly used method to study material fracture problems. G/XFEM has the following advantages in analyzing fracture problems. First, the finite element mesh is independent of the crack surface. Second, there is no need to refine the mesh when dealing with the problem of high stress concentration near the crack tip. Third, its mesh only needs to be slightly densified near the crack surface and crack tip, achieving a great improvement. In 2003, Stefano and Umberto [5] analyzed quasi-brittle fracture by G/XFEM. Kumar et al. [6] investigated the bending crack problem by using the virtual joint G/XFEM in 2015. Besides, Jia and Li [7] used G/XFEM to study the crack problem of composite materials and analyzed the influence of grid density, strengthening function scheme, Gauss integral scheme and geometric enrichment radius on SIFs in 2015. Moreover, Liu et al. [8] applied multiscale G/XFEM to analyze the interaction between fatigue crack growth and microdefects in 2018.
However, the G/XFEM is difficult to deal with the blending element because the number of conditions of the generated system matrix is generally large, resulting in ill-conditioned equations [9,10]. Therefore, Babuška and Banerjee [11] first proposed the SG/XFEM in 2011. Afterward, Gupta et al. [12] analyzed the condition number and calculation accuracy of SG/XFEM in 2015 and employed it to solve the problem of three-dimensional fracture mechanics accurately. In 2019, de Oliveira et al. [13] addressed the problem of plane elastic crack with SG/XFEM and obtained the value of the SIF of the crack tip.
However, most of the above research works focused on the simple fracture problem with the single-mode crack in elastic materials. There is little research on the fracture problem of materials with mixed-mode cracks (such as I + II). Therefore, the research on SG/XFEM with mixed-mode cracks has certain social and economic value.
The novelty of this paper is performing fracture analysis via SG/XFEM in isotropic materials, and SIFs around crack tips and scaled condition number of the system matrix are obtained by domain formulation of the interaction integral method combined with SG/XFEM. On this basis, several numerical examples are given to validate the accuracy of results with different enrichment schemes.
The rest of the paper is organized as follows. In Section 2, some related theories and issues are introduced for G/XFEM and SG/XFEM and scaled condition matrix and mixed element. Then, the basic theory for fracture of SG/XFEM and some numerical examples are provided in Section 3. Some numerical examples are given in Section 4. Finally, a conclusion is drawn in Section 5.

REL ATED THEORIES AND ISSUES 2.1 Introduction to generalized/extended finite element and stable generalized/extended finite element
The generalized/extended finite element is used to add the discontinuity term into the displacement approximation of the finite element in the frame of the finite element. The steps of the fracture analysis of the generalized/extended finite element are provided as follows: (1) Mesh generation can be independent of component geometry.
(2) Generally, the level set method is used to determine the crack location and track the crack growth path. (3) Based on the partition of the unity approach, the element displacement interpolation functions near the crack surface and around the crack tip are improved, and the Lagrangian shape function is applied as the partition of unity function.
(4) The displacement field and temperature field are acquired by the FEM, the SIF at the crack tip is obtained by post-processing and then the crack growth is obtained by the crack growth conditions and fracture criteria.
A G/XFEM shape function φ ji depends on the node x j coordinates and the corresponding patch of elements ω j and is defined by the product of a set of linearly independent enrichment function sets {L ji (x)} q j i=1 and finite element shape functions N k , namely where forms a basis of local space χ j (ω j ), which is called patch approximation space. The test and trial spaces of G/XFEM are provided by Gupta et al. [12]: where G/XFEM also has weaknesses. Specifically, there is an illconditioned system matrix [11] and robustness is lacking. Therefore, Babuška and Banerjee proposed the SG/XFEM [11] in 2011. The biggest difference between SG/XFEM and G/XFEM is that the former is based on the latter, locally modifying the enrichment function to form a new patch of element approximation space. Besides, the value of the enrichment function at the patch of element node in SG/XFEM is equal to zero. The enrichment function [14] modified by SG/XFEM is expressed asL where I ω j denotes the piecewise finite element interpolation function ω j of the enrichment function L ji of the patch of element, which is determined by where x k represents the coordinate vector of nodes on the element e, n e is the number of nodes of the element e and L ji (x k ) is determined by Eq. (2).
If the global enrichment space associated with the patch of element approximation space isS ENR , the SG/XFEM trial space can be expressed as According to Gupta et al. [14], S FEM andS ENR are almost orthogonal about the inner product of their energy.
InS ENR , the shape function is composed of the product of the partition of unity form function and the modified enrichment function: φ ji where j has no summation.

Ill-conditioned matrix and mixed element
As mentioned earlier, G/XFEM still has its disadvantages: illconditioned system matrix problem and blending element problem, which are overcome by SG/XFEM. In SG/XFEM and G/XFEM, the condition number of the system matrix [12,14] is determined by the scaled condition number k, which is usually used for GFEM.
The condition number C of the system matrix K is defined as the ratio between its maximum (λ max ) and minimum eigenvalues (λ min ), as expressed in the following equation [15]: Introducing the scaled stiffness matrixK: where D represents a diagonal matrix, its matrix elements are The condition number of the scaled matrix is the condition number of scale (k) [14], which is defined as The number of scaled conditions k in G/XFEM is proportional to the negative fourth power of the mesh size of the finite element. The number of conditions generated by G/XFEM or SG/XFEM depends on the approximate linear correlation of the enrichment basis function used [12].
The mixed element is the element with some nodes enriched. The mixed element often leads to large discrete errors and poor convergence. Therefore, G/XFEM mixed element problems of both opening mode I fracture and blending fracture involve two parameters: element strain energy and SIF [16]. Among them, the first parameter is applied to measure the influence of SG/XFEM used in the blending element region on the reduction of discrete error; the second parameter SIF is an essential index in fracture mechanics. In this study, SIF is calculated according to the interactive integration method [3]; then, the blending element problem is studied.

FR ACTURE THEORY OF STABLE GENER ALIZED/E XTENDED FINITE ELE MENT 3.1 Problem model
A plane region = ∪ ∂ with linear cracks has the boundary ∂ = ∂ u ∪ ∂ σ , ∂ u ∩ ∂ σ = 0, of the empty set. The indexes u and σ represent Dirichlet and Neumann boundary conditions, respectively, as expressed in the following equations: where σ denotes Cauchy stress tensor, n represents the unit external normal vector of ∂ σ andt refers to specified external force on ∂ σ .
The equilibrium equation and the constitutive equation are provided as follows: where C represents Hooke's tensor and Є denotes small strain tensor.
An approximate solutionũ of the boundary value problem (13)- (15) can be expressed as whereχ 0 ( ) is the discretization of Hilbert space H 1 ( ).

Fracture theory and numerical simulation
It is necessary to enrich the vicinity of the crack surface to reflect the crack discontinuity in order to simulate the fracture behavior of cracked plates. In G/XFEM, the commonly used function to describe the crack discontinuity line is the Heaviside function, namely whereȳ represents the relative position between the point in the crack local coordinate system and the crack surface. Besides, the elements near the crack tip for mode I fracture require to be enriched to describe the singular stress behavior near the crack tip. The following enrichment functions are used, representing the first term of the exact displacement solution near the crack tip of linear elastic mode I crack [17]: where G denotes shear modulus, Q = 1/3, λ = 0.5, k = (3 − ν)/(1 + ν), and θ and r represent crack local coordinates. For mode II fracture, the following enrichment functions near the crack tip are applied [18]: In SG/XFEM, we use the linear Heaviside function [14] in another global coordinate system to reflect the displacement discontinuity of the crack surface; it is related to the element slice ω j , namely where H is determined by Eq. (18) and h j refers to the maximum size of the element around the node x j . The last two terms at the right end of the above formula are employed to reduce the error u − I ω j (u) of the element piece ω j [14]. The third term in Eq. (23) is redundant, and the first two terms are enough to reduce the patch of element errors. Therefore, a modified linear Heaviside function [13] is introduced in this study: The Heaviside function for G/XFEM enrichment is only enriched on the crack surface. However, the Heaviside function for SG/XFEM is enriched on both the crack surface nodes and a layer of nodes below the crack surface [14] if it is used. Simultaneously, the element piece ω j under the crack surface should be enriched by Eq. (23) to reduce the error of u − I ω j (u). In this paper, three kinds of Heaviside functions are used, including Standard (18), Linear (23) and Linear Mod (24).
For the enrichment near the crack tip, the same enrichment scheme is adopted for SG/XFEM and G/XFEM. The enrichment function given by Eqs (19) and (20) or Eqs (21) and (22) is only needed for the element node containing the crack tip. For the enrichment of mode I fracture tip, the geometric enrichment scheme is used in this study.
Thus, substituting Eqs (1), (2), (4)-(7), (18), (23), (24) and (19)-(22) into Eq. (8), the SG/XFEM displacement approximation can be obtained as where I is all the set of nodes in the finite element domain, N cr is the set of nodes whose support is crossed by the crack faces, and N tip is the set of nodes inside a fixed area around the crack tip (geometry enrichment). The selection of enriched node domain (crack surface enrichment or crack tip geometry enrichment) for two-dimensional crack problem is illustrated in Fig. 1. u i are usual finite element nodal displacements, and a n and b k are enriched nodal displacement freedoms. f cr (x) is the SG/XFEM crack surface enrichment function that is represented as Eq. (18)    In order to reduce the computational cost, the crack lip is aligned with the mesh; thus, the surface and tip enrichment schemes of mode I single edge crack obtained by SG/XFEM are illustrated in Fig. 1 [13].
The discontinuity of the crack is represented by duplicate nodes. All the following numerical simulations adopt rectangular elements and Gauss integration. Specifically, when the element node intersects with the crack, 12 Gauss integration points are used in each direction of the element, 4 Gauss points are used in each direction of the other enrichment areas and the minimum Gauss point integration scheme is used in the element area without enrichment points. SIF is an efficient and accurate interactive integration method [13]; some details on the SIF extraction process can be seen in [13].
In this paper, five enrichment schemes are applied. For the convenience of expression, G/XFEM (Standard) is adopted to represent the geometric enrichment scheme of extended finite element for crack tip enrichment, standard [Standard (18)] Heaviside function is selected for crack surface enrichment and G/XFEM (Linear) is employed to denote the

NUMERICAL E X A MPLES 4.1 Example 1 (mixed-mode crack under unilateral uniformly distributed pressure load)
A thin plate with an inclined crack is illustrated in Fig. 2; the same problem has been studied in the literature [13,17], in which only the loading mode is different from ours. The loading mode in the literature is uniformly distributed loading on the upper edge. The width of the plate is 1, the height is 4.5 and the plate contains a crack with an inclined angle of 135°and a length of 0.5. The upper edge of the plate bears the linear uniformly distributed pressure of 1. Assuming the thickness of the plate is 1, all units are unified and the plate is under plane stress condition, the material properties are isotropic: Poisson's ratio is 0.3 and Young's modulus is 1. There are three kinds of finite element mesh: mesh 1 is 7 × 11 mesh, mesh 2 is 14 × 22 mesh and mesh 3 is 28 × 33 mesh. Mesh 1 is exhibited in Fig. 3. The SIFs of modes I and II (K I ,K II ) at the crack tip under different enrichment combination schemes and different grid densities are listed in Tables 1 and 2. The reference values of SIFs [17] are selected, respectively. Table 3 displays each method's K II absolute error compared to the reference value.
As observed in Tables 1-3, the error between the SIF and the reference value gets smaller and smaller with the mesh refinement. The results in this paper are relatively accurate, and SG/XFEM (Linear) gives the best results. Second, G/XFEM (Standard) and G/XFEM (Linear) also provide the better results, while SG/XFEM (Standard) and SG/XFEM (Linear Mod) exhibit the worst results. These results verify that the stable generalized/extended finite element can reduce the numerical error of the mixed element (both enriched and unenriched nodes exist in the elements). Table 4 displays the scaled condition number (12), the total freedom of finite element and the time (T) taken to assemble the stiffness matrix of the finite element under different enrichment combination schemes corresponding to mesh 1. It can be revealed that the scaled condition number of the system matrix of SG/XFEM is lower than that of G/XFEM, verifying that the correctness of the method used in this paper and the numerical calculation is more robust. In terms of the time spent in assembling the stiffness matrix, SG/XFEM takes more time compared with G/XFEM; this is caused by the extra calculation of SG/XFEM enrichment functions (6) and (7). Figures 4 and 5 demonstrate the convergence of mesh refinement on the normalized mode I and mode II SIFs at the crack tip by different enrichment schemes, respectively. Three different mesh refinements are used: mesh 1 is 7 × 11 mesh, mesh 2 is 14 × 22 mesh and mesh 3 is 28 × 33 mesh. From Figs 4 and 5, it can be seen that the convergence of the results is very good and the SIFs gradually converge to the reference value with mesh refinement.

Example 2 (mixed-mode crack under shear load)
As illustrated in Fig. 6, a thin plate with an edge crack has been studied in the literature [2] by standard G/XFEM; besides, this problem has also been explored using the elastic theory. The width of the plate is 7, the height of the plate is 16, the length of the edge cracks is 3.5, the thickness of the plate is 1, the Poisson's ratio is 0.3 and the Young's modulus is 1. It is assumed that all units are unified, the plate is in the plane strain state, a uniformly distributed shear loading with the size of 1   is loaded on the top edge of the plate and the x, y direction displacement of all nodes at the lower edge of the plate is constrained to 0. There are two kinds of finite element grids: grid 1 contains 60 cells and 77 nodes; grid 2 contains 240 cells and 273 nodes. Two nested grids and enrichment schemes used in the calculation are presented in Fig. 7.
Tables 5 and 6 provide the regularized SIFs of modes I and II [2] at the crack tip under different enrichment com-bination schemes and different grid densities. The reference values of normalized SIFs (K I = 34.00, K II = 4.55) [19] are selected. Table 7 displays each method's absolute error for regularized SIFs of mode II compared to the reference value.
It can be revealed from Tables 5-7 that the error between the obtained SIF and the reference value gets smaller and smaller as the grid is densified. The result is similar to that of Example 1, which is still the best result provided by SG/XFEM  (Linear); meanwhile, the worst result is offered by SG/XFEM (Standard) and SXFEM (Linear Mod), as stated by Gupta et al. [14]. These results verify that the existing enrichment nodes in the mixed element (element) are reduced by the SG/XFEM. Table 8 exhibits the scaled condition number (12), the total freedom of finite element and the time T (in milliseconds) of assembling the stiffness matrix of the finite element under dif- ferent enrichment combination schemes corresponding to mesh 1. The proportion condition number of the system matrix of SG/XFEM is lower than that of G/XFEM, demonstrating that the correctness of the method used in this paper and the numerical calculation is more robust. In terms of the time spent in assembling the stiffness matrix, SG/XFEM takes significantly more time compared with G/XFEM; this is caused by the extra calculation of SG/XFEM enrichment functions (6) and (7).
It can be observed from Table 8 that the scaled condition number of the system matrix of SG/XFEM is greatly reduced compared with G/XFEM, suggesting that the numerical method used in this paper is more robust. From the perspective of the time spent in assembling the stiffness matrix, SG/XFEM takes significantly more time compared with G/XFEM due to the extra calculation of SG/XFEM enrichment functions (6) and (7).

CONCLUSIONS
This paper aimed to investigate the fracture problem of composite crack. First, the SG/XFEM was adopted to address the fracture problem of a cracked plate. Then, the developed crack strengthening function was applied to the problem of composite crack. Finally, the value of SIF at the tip of the component with mode I + II crack, the condition number of the system matrix and the numerical calculation time were obtained. The following conclusions can be drawn by the simulation: (1) The linear Heaviside function ensures the accuracy of the SG/XFEM. (2) The SG/XFEM is more robust and the condition number of the system matrix is greatly reduced compared with the XFEM; besides, the scaled condition number of the system matrix in SG/XFEM has the same order as that of the FEM. (3) It takes longer to assemble the stiffness matrix and calculate SIF compared with the XFEM. This is mainly due to the additional shape function calculation required by the SG/XFEM. The calculation efficiency of the SG/XFEM can be improved by encrypting the mesh. (4) SG/XFEM is robust because the conditioning of the system matrix does not worsen as the crack approaches the boundaries of elements.