Convergence Investigation of XFEM Enrichment Schemes for Modeling Cohesive Cracks

: When simulating cohesive cracks in the XFEM framework, speciﬁc enrichment schemes are designed for the non-singular near-tip ﬁeld and an iteration procedure is used to solve the nonlinearity problem. This paper focuses on convergence and accuracy analysis of XFEM enrichment schemes for cohesive cracks. Four different kinds of enrichment schemes were manufactured based on the development of XFEM. A double-cantilever beam specimen under an opening load was simulated by Matlab programming, assuming both linear and exponential constitutive models. The displacement and load factors were solved simultaneously by the Newton–Raphson iterative procedure. Finally, based on a linear or an exponential constitutive law, the inﬂuences of variations in these enrichment schemes, including (i) specialized tip branch functions and (ii) corrected approximations for blending elements, were determined and some conclusions were drawn.


Introduction
In quasi-brittle materials, such as geomaterials and concrete, the fracture behavior is quite different from that of brittle materials. A fracture process zone (FPZ) of negligible size develops at the crack front due to plasticity or micro-cracking [1]. The assumption of linear elastic fracture mechanics (LEFM) is quite restrictive for certain types of failure, where the nonlinear zone ahead of the crack tip is negligible in comparison with the dimension of the crack. Employing LEFM may produce dangerous results for fracture propagation in quasi-brittle materials [2]. Therefore, cohesive crack models have been developed to analyze metal materials. Hillerborg et al. [3] introduced fracture energy into the cohesive crack model and proposed various traction-separation relations for concrete. The cohesive crack models have been extensively used in studies on the FPZ and nonlinear failure in engineering structures.
Within the FPZ ahead of the crack tip, although damage develops to a certain degree, cohesive stress can still be transferred. In the cohesive model, the nonlinear FPZ, where degradation or damage mechanisms occur as a result of micro-cracking or micro-voids, ahead of the crack tip is lumped into a discrete line [4,5]. This stress-softening type of behavior in the FPZ is represented by a cohesive constitutive relation [6]. The FPZ is characterized by two crack tips: a mathematical (or fictitious) crack tip and a physical one. As shown in Figure 1, at the mathematical tip, the crack opening is zero and the cohesive traction equals the tensile strength of the material, while, at the real crack tip, the crack opening equals the critical crack opening and the cohesive traction is zero. The crack smoothly closes from the physical tip to the fictitious tip, and the drawback of infinite stress due to the LEFM theory is avoided [7]. Over the last two decades, the methodological development of t element method (XFEM) has led to a phenomenal increase in applicat framework, localized phenomena are modeled by incorporating a prior the solution into the FEM approximation using a partition of the unity p ture propagation can be handled even on a structured mesh by dynami pre-selected local approximation spaces in the problem domain. To i approximation space, specialized enrichment functions and correspo freedom are added onto local existing FE nodes. In contrast to FEM, it straints such as mesh conformance to physical discontinuities, mesh r the crack tip, and burdensome adaptive remeshing whenever the cra enrichment schemes have been specialized to apply the XFEM in mode problems, such as bi-material [8,9], three-dimensional crack [10,11], i [12], microcrack [13,14], two-phase flow [15,16], and frictional contact These applications have reached a high degree of robustness and are n rated into general software such as LS-DYNA and ABAQUS.
For linear elastic fracture simulation, two types of enrichment func [19,20]: (i) heaviside functions, which model the jump in the displacem crack surface; and (ii) tip branch functions, which are derived from the t of the displacement field in the neighborhood of the crack tip, to captu The XFEM provides higher numerical accuracy than FEM without an refinement. However, the convergence rate with respect to the mesh improved, because the enriched zone becomes smaller as the mesh is r other hand, researchers [22] have found that the parasitic terms present elements may drastically reduce the convergence rate of XFEM appr ever, the influence of the parasitic terms cannot easily be predicted [23 ments or problems, such as bi-material structures, they may reduce the while for others, such as strong discontinuities, they may only increase while keeping the convergence rate unchanged. In order to reach the op rate, special treatments to blending elements have been proposed, su riched and standard regions [24], hierarchical shape functions [25], th Over the last two decades, the methodological development of the extended finite element method (XFEM) has led to a phenomenal increase in applications. In the XFEM framework, localized phenomena are modeled by incorporating a priori knowledge about the solution into the FEM approximation using a partition of the unity property. The fracture propagation can be handled even on a structured mesh by dynamically adjusting the pre-selected local approximation spaces in the problem domain. To incorporate a local approximation space, specialized enrichment functions and corresponding degrees of freedom are added onto local existing FE nodes. In contrast to FEM, it relaxes mesh constraints such as mesh conformance to physical discontinuities, mesh refinement around the crack tip, and burdensome adaptive remeshing whenever the crack grows. Various enrichment schemes have been specialized to apply the XFEM in modeling discontinuity problems, such as bi-material [8,9], three-dimensional crack [10,11], inclusion and void [12], microcrack [13,14], two-phase flow [15,16], and frictional contact [17,18] problems. These applications have reached a high degree of robustness and are now being incorporated into general software such as LS-DYNA and ABAQUS.
For linear elastic fracture simulation, two types of enrichment functions are required [19,20]: (i) heaviside functions, which model the jump in the displacement field across the crack surface; and (ii) tip branch functions, which are derived from the theoretical solution of the displacement field in the neighborhood of the crack tip, to capture the singularity. The XFEM provides higher numerical accuracy than FEM without any significant mesh refinement. However, the convergence rate with respect to the mesh parameter h is not improved, because the enriched zone becomes smaller as the mesh is refined [21]. On the other hand, researchers [22] have found that the parasitic terms presented in the blending elements may drastically reduce the convergence rate of XFEM approximations. However, the influence of the parasitic terms cannot easily be predicted [23]; for some enrichments or problems, such as bi-material structures, they may reduce the convergence rate, while for others, such as strong discontinuities, they may only increase the absolute error while keeping the convergence rate unchanged. In order to reach the optimal convergence rate, special treatments to blending elements have been proposed, such as coupling enriched and standard regions [24], hierarchical shape functions [25], the enhanced strain technique [26], and the corrected XFEM [23]. Among these techniques, the corrected XFEM may be the easiest to implement while producing the optimal result.
For cohesive crack simulation, enrichment schemes used for traction-free cracks are no longer suitable. This is because, in the cohesive crack model, the singularity of the crack tip field vanishes. The XFEM enriches the standard local FE approximations with prior known information about the problem. The cohesive crack model abandons the singularity of the crack tip stress field, which is an unrealistic assumption of LEFM. Therefore, new enrichment functions have to be designed in the XFEM framework to simulate the true asymptotic field for the cohesive crack model [2]. To date, various enrichment schemes have been developed for modeling cohesive cracks. Originally, only a heaviside function Mathematics 2022, 10, 383 3 of 17 was employed. Because the singularity vanishes in the near-tip field, the heaviside function can be suitable for the entire crack, including the crack tip. This approach is used by Wells and Sluys [27]. However, if only the heaviside function is applied to all nodes, the crack is restricted to ending at the element edges to ensure that the jump in the displacement field at the mathematical crack tip equals zero. The approaches given in Duarte et al. [28] and Zi and Belytschko [29] overcome this deficiency by modification of the shape functions within the tip element, so that the crack tip can lie within the element. Mariani and Perego [30] proposed enrichment functions as a product of the heaviside function and ramp functions. Some references provide special tip branch functions for cohesive cracks. Moës and Belytschko [31] suggest the following tip branch function: φ(r, θ) = r k sin θ 2 , with k being either 1, 1.5, or 2. Other enrichment functions based on analytical considerations are given by Cox [32]. Meschke and Dumstorff [33] use four tip branch functions similar to those for traction-free cracks, only replacing √ r with r, e.g., φ α (r, θ) = {r sin θ 2 , r cos θ 2 , r sin θ 2 sin θ, r cos θ 2 sin θ}. With the employment of tip branch functions, the crack can end arbitrarily within the element. However, a loss of the partition of unity in the blending elements may lead to a reduction in accuracy. Convergence and accuracy studies of these enrichment schemes are needed for a suitable choice.
As far as convergence rates are concerned, when numerically simulating traction-free crack by the XFEM, the factors that influence the convergence rate include the enrichment zone size [21], the shape function polynomial order [24], the special treatment of the blending elements, and the choice of enrichment functions. V. Gupta et al. [34] studied the influence of enrichment zone size on convergence rate and found that, for traction-free crack simulation, the convergence rate is controlled by the stress gradient outside the enrichment zone and the error is caused by the blending element. When it comes to the cohesive crack problem, the smoother stress gradient and the nonlinearity of the governing equation make the accuracy and convergence properties new problems that require study.
In this paper, we focus on investigating the accuracy and convergence properties of different enrichment schemes for cohesive crack simulation. A numerical simulation was conducted on a double-cantilever beam specimen, assuming a linear or an exponential constitutive law, in order to provide useful information for the choice of enrichment scheme for cohesive crack simulation. The enrichment schemes we considered can briefly be stated as follows.
(i) XFEM-h. Only the heaviside function is used, with a small modification of the shape functions in the tip element. (ii) XFEM-s. Both the heaviside function and the tip branch function φ(r, θ) = r sin θ 2 are used. (iii) XFEM-c1. Only φ(r, θ) = r sin θ 2 is used as the tip branch function, and a corrected approximation for φ(r, θ) = r sin θ 2 is used in the blending elements. (iv) XFEM-c2. Both φ(r, θ) = r sin θ 2 and φ(r, θ) = r cos θ 2 are used as tip branch functions in addition to a corrected approximation for these two tip branch functions in blending elements.
The remainder of this paper is organized as follows. A description of the cohesive crack problem domain and the XFEM formulation for the cohesive crack problem are provided in Section 2. Information on these four enrichment schemes is provided in Section 3. Section 4 presents the Newton iterative algorithm for solving the nonlinear problem. In Section 5, numerical results of convergence and accuracy studies on these enrichment schemes for different cohesive constitutive models are presented. The effect of variations in these enrichment schemes is investigated. Finally, our main conclusions are summarized as Section 5.

Model Problem Definition
Consider a two-dimensional domain Ω crossed by a cohesive discontinuity, as shown in Figure 2. The strong form of the equilibrium equation of this body can be expressed as where ∇ is the gradient operator, σ is the Cauchy stress, and b is the body force. The behavior of the bulk material is assumed to be linearly elastic, and the constitutive relation is defined as σ = D·ε. The essential and natural boundary conditions are presented as follows where n Γ is the outward unit normal vector to the external boundary Γ, t is the prescribed load vector on the boundary Γt, u is the the prescribed displacement on the boundary Γ u , and f c is the cohesive traction transferred across the Γc, which is related to the displacement gap ω across the discontinuity according to the stress softening model.

Model Problem Definition
Consider a two-dimensional domain Ω crossed by a cohe in Figure 2. The strong form of the equilibrium equation of th

Discretization of Governing Equations
In the XFEM, the displacement discontinuity can be directly embedded by introducing additional degrees of freedom onto existing nodes whose supports are intersected by discontinuities. Comprehensive overviews of the XFEM have been given by numerous studies [35][36][37].
The generalized form of the XFEM approximation of the displacement field can be written as In the above function, the standard FE approximation ∑ A N i (x)·u i represents the continuous part of the displacement field, while the second term represents the discontinuous part, where u i and a j are standard and enriched DOFs, respectively, φ j (x) are the enrichment functions, which take different forms for specific kinds of discontinuity problems, A denotes the set of all nodes, and J denotes the pre-selected set of local nodes associated with discontinuities. The weak form of the governing partial differential equation can be derived from the principle of virtual work and the Galerkin procedure. When the cohesive crack model is assumed, the cohesive traction fc that transferred is a function of the crack opening ω. The weak form of the equilibrium equation can be expressed as: Discretization of Equation (5) in the XFEM framework results in: where q is the generalized nodal displacement vector, q e = [ u e i a e i ] for each element, and λ is the load factor.
The stiffness matrix K is composed of With where T c is the tangential modulus matrix of the cohesive crack determined by the cohesive crack behavior and is obtained from the relation The external nodal force f ext and the cohesive nodal force f coh can be obtained as where the crack opening displacement ω can be given by It can be observed from the Equations (8) and (9) that the cohesive behavior has a direct effect on both the stiffness matrix K and the nodal force vector f coh . The relation between the cohesive force and the crack opening makes the problem nonlinear.
The four enrichment schemes we examined are designed to consider the effect of their variations, including the employment of tip branch functions and a corrected approximation in blending elements. These four schemes are denoted XFEM-h, XFEM-s, XFEM-c1, and XFEM-c2, and detailed as follows.

XFEM-h
Because the singularity of the displacement field around the crack tip vanishes, a heaviside function is suitable for the entire crack, including the crack tip. In this scheme, the approximation of the displacement field can be written as where H(x) is the heaviside function, which takes +1 on one side of the crack and −1 on the other side, and J is the set of nodes whose supports are fully cut by the crack, which is depicted in Figure 3a.
Mathematics 2022, 10, 383 where H(x) is the heaviside function, which takes +1 on one side of the crack and −1 on the other side, and J is the set of nodes whose supports are fully cut by the crack, which is depicted in Figure 3a. In order to make the displacement gap vanish to zero at the crack tip within the tip element, we extended the method proposed by Zi and Belytschko [29] to quadrilateral elements. Specifically, for the tip element, the modified shape function ( ) j N x was used instead of the standard shape function Figure 4, if the crack intersects with boundary 14 within the tip element, we make a straight line through the crack tip point and intersect the element boundary at points 5 and 6. Then, the shape function ( ) j N x used for the tip element is actually the standard shape function of virtual element 1564. Since nodes 1 and 4 are enriched, the discontinuous part of the displacement can be written as where x* are the coordinates of virtual element 1564. Since this scheme treats the entire domain with the heaviside function only, the blending with the unenriched subdomain does not occur, which implies that the PU holds in the entire domain. In order to make the displacement gap vanish to zero at the crack tip within the tip element, we extended the method proposed by Zi and Belytschko [29] to quadrilateral elements. Specifically, for the tip element, the modified shape function N j (x) was used instead of the standard shape function N j (x). As shown in Figure 4, if the crack intersects with boundary 14 within the tip element, we make a straight line through the crack tip point and intersect the element boundary at points 5 and 6. Then, the shape function N j (x) used for the tip element is actually the standard shape function of virtual element 1564. Since nodes 1 and 4 are enriched, the discontinuous part of the displacement can be written as where x* are the coordinates of virtual element 1564.

XFEM-s
Another way to allow the crack tip to be located arbitrarily is to employ branch functions. For traction-free cracks, the branch functions are chosen based on the analytical solution of the displacement field in the vicinity of the crack tip, that is However, the combination of these functions does not produce the non-singular stress field at the tip of the cohesive crack. On the basis of the analytical solution of the cohesive crack problem, some researchers proposed  Since this scheme treats the entire domain with the heaviside function only, the blending with the unenriched subdomain does not occur, which implies that the PU holds in the entire domain.

XFEM-s
Another way to allow the crack tip to be located arbitrarily is to employ branch functions. For traction-free cracks, the branch functions are chosen based on the analytical solution of the displacement field in the vicinity of the crack tip, that is However, the combination of these functions does not produce the non-singular stress field at the tip of the cohesive crack. On the basis of the analytical solution of the cohesive crack problem, some researchers proposed that only one non-singular enrichment function be used for the two-dimensional cohesive crack tip, which takes the following form.
Others proposed that four branch functions be used to enrich the tip element, which is In this enrichment scheme, r sin θ 2 is used as a branch function, which is presented in Figure 5a. It is obvious that this branch function is suitable for capturing the displacement gap at the crack tip.

XFEM-c1 and XFEM-c2
When branch functions are used in conjunction with a heaviside function, the partition of the unity property does not hold in the blending elements. As shown in Figure 3, in blending elements only some of the nodes are enriched, which means that In addition, the branch function is not a piecewise constant function like the heaviside function, so the parasitic terms resulting from do not vanish at the edges of the tip element. The presentation of parasitic terms in blending elements can reduce the convergence rate and accuracy [38]. Fries T.P. [23] proposed a corrected approximation in which all nodes in blending elements are enriched and the enrichment functions are modified, solving the problem most efficiently. The approximation of the displacement field of the corrected XFEM can be written as: where L is the set of second-layer nodes around the tip element, as marked with triangles in Figure 3c. In this scheme, the set of nodes L is also enriched with branch functions and is a ramp function, which is defined as follows and depicted in Figure 6.
It can be seen in Figure 6, within the tip element, that we have R(x) = 1, while within the blending element, the ramp function varies continuously and deceases to zero at the In this scheme, the XFEM approximation of the displacement field can be expressed as As marked in Figure 3b, J is the set of nodes whose supports are intersected with the crack, and M is the set of nodes whose supports contain the crack tip. If a node simultaneously belongs to J and M, then it belongs to M.

XFEM-c1 and XFEM-c2
When branch functions are used in conjunction with a heaviside function, the partition of the unity property does not hold in the blending elements. As shown in Figure 3, in blending elements only some of the nodes are enriched, which means that ∑ k N k (x) = 1.
In addition, the branch function is not a piecewise constant function like the heaviside function, so the parasitic terms resulting from ∑ k N k (x)φ(x) do not vanish at the edges of the tip element. The presentation of parasitic terms in blending elements can reduce the convergence rate and accuracy [38]. Fries T.P. [23] proposed a corrected approximation in which all nodes in blending elements are enriched and the enrichment functions are modified, solving the problem most efficiently. The approximation of the displacement field of the corrected XFEM can be written as: where L is the set of second-layer nodes around the tip element, as marked with triangles in Figure 3c. In this scheme, the set of nodes L is also enriched with branch functions and additional DOFs b α k . R(x) is a ramp function, which is defined as follows and depicted in Figure 6.
Mathematics 2022, 10, x FOR PEER REVIEW Figure 6. The value of R(x) on a discretized mesh.
In order to investigate the effect of an additional singular branch function in Equ 14, only 2 sin r  was employed in XFEM-c1, while both 2 sin r  and 2 cos r  wer ployed in XFEM-c2. It can be seen from their plots in Figure 5 that the branch fun 2 cos r  can help to capture the stress gradient at the rear of the crack tip. In both sch a corrected approximation in blending elements is used to eliminate the negative ence of parasitic terms.

Nonlinear Algorithm
In order to guarantee the smooth closing of the crack as required by the definit the cohesive crack model, one more condition is required. This condition is usua ferred to as the zero stress intensity factor condition. It is assumed that the crack p gates under the mode Ⅰ loading condition, so only the mode Ⅰ stress intensity facto is taken into account, i.e., A SIF at the crack tip equal to zero implies that the stress component normal crack tip is finite [39]. Alternatively, smooth closure can also be guaranteed by a condition, where the stress projection in the normal direction n  of the crack is eq the tensile strength of the material, i.e., where tip  is the stress at the crack tip and t f is the tensile strength of the materia zero SIF condition and the stress condition can be used interchangeably, but the condition is simpler to implement; therefore, it was adopted in this paper. The discretized form of the stress condition can be written as T Figure 6. The value of R(x) on a discretized mesh.
It can be seen in Figure 6, within the tip element, that we have R(x) = 1, while within the blending element, the ramp function varies continuously and deceases to zero at the element edges. After this modification, the PU holds everywhere in the domain. Improved convergence rates were verified in applications to bi-material problems, while in other applications only the error level was reduced.
In order to investigate the effect of an additional singular branch function in Equation (14), only r sin θ 2 was employed in XFEM-c1, while both r sin θ 2 and r cos θ 2 were employed in XFEM-c2. It can be seen from their plots in Figure 5 that the branch function r cos θ 2 can help to capture the stress gradient at the rear of the crack tip. In both schemes, a corrected approximation in blending elements is used to eliminate the negative influence of parasitic terms.

Nonlinear Algorithm
In order to guarantee the smooth closing of the crack as required by the definition of the cohesive crack model, one more condition is required. This condition is usually referred to as the zero stress intensity factor condition. It is assumed that the crack propagates under the mode I loading condition, so only the mode I stress intensity factor (SIF) is taken into account, i.e., where K Itip is the mode I SIF calculated at the crack tip. In FEM-based methods, the SIF can be obtained by means of the domain integration method. A SIF at the crack tip equal to zero implies that the stress component normal to the crack tip is finite [39]. Alternatively, smooth closure can also be guaranteed by a stress condition, where the stress projection in the normal direction n Γ of the crack is equal to the tensile strength of the material, i.e., where σ tip is the stress at the crack tip and f t is the tensile strength of the material. The zero SIF condition and the stress condition can be used interchangeably, but the stress condition is simpler to implement; therefore, it was adopted in this paper. The discretized form of the stress condition can be written as where S = M T · D · B is the operator by which the stress at the crack tip is calculated, and M = n Γ ⊗ n Γ is the Voigt notation. It is obvious that in the equilibrium condition (Equation (6)), as the cohesive force depends on the crack opening ω, the problem is nonlinear. The scheme recommended in [29] was employed, we combined Equation (6) and Equation (20), and q and λ were solved simultaneously by the Newton-Raphson iterative procedure. The residual vector of the governing equation is given by where the independent unknowns are q and λ. The Jacobian matrix is is the additional stiffness term effective on the crack surface in the FPZ, which can be obtained by At the ith iteration, the increments in independent variables obtained from Equations (21) and (22) are ∆q ∆λ

Numerical Study
In this study, a double-cantilever-beam (DCB) specimen containing a level cohesive crack was numerically simulated by the above four different enrichment schemes in order to examine their accuracy and convergence performance. This configuration was taken from the literature [29]. The boundary conditions and dimensions of this case study are provided as a sketch in Figure 7. Uniformly distributed forces were applied on the left side of the beam, and the plane stress condition was assumed to hold. The Young's modulus was 36.5 GPa and the Poisson's ratio was 0.18, which are the material properties of common concrete. ext coh ( ) where the independent unknowns are q and λ. The Jacobian matrix is is the additional stiffness term effective on the crack surface in t which can be obtained by

Numerical Study
In this study, a double-cantilever-beam (DCB) specimen containing a level c crack was numerically simulated by the above four different enrichment schemes to examine their accuracy and convergence performance. This configuration w from the literature [29]. The boundary conditions and dimensions of this case st provided as a sketch in Figure 7. Uniformly distributed forces were applied on side of the beam, and the plane stress condition was assumed to hold. The Young ulus was 36.5 GPa and the Poisson's ratio was 0.18, which are the material prop common concrete.

Linear Softening Model
For the cases with a linear softening model, the cohesive force can be expressed as Six finite element meshes were used in the convergence studies (17 × 9, 31 × 15, 61 × 31, 121 × 61, 301 × 151, and 601 × 301 grids of quadrilateral elements). The meshes were created by making the element length in the x and y directions approximately equal, with the number of elements being odd, such that the crack tip lies within an element. The mesh size h was represented by the square root of the area of an element.

Linear Softening Model
For the cases with a linear softening model, the cohesive force can be expressed as In this paper, the material properties of common concrete were used, where the tensile strength f t = 3.18 MPa and the critical crack opening ω c = 0.0314 mm. The fracture energy was E f = 0.5 f t ω c = 50 N/m. The plot of this softening model is provided in Figure 8.
In this paper, the material properties of common concrete we sile strength 3.18MPa  Due to the nonlinearity of the problem, an explicit analytical so ment field around the crack tip is not available. In order to evaluate for different mesh densities, we took the results obtained by the f enrichment scheme as reference exact solutions. The h-version c finite element method was quantified by means of the relative erro was calculated by the following equations. In Figure 9, the deformed geometry of the cohesive crack pro that of the traction-free crack problem when the same load factor hesive force exists between crack faces, the crack closes smoothly f the fictitious tip. The Contour plots of σyy for the cohesive crack pr free crack problem are provided in Figure 10. It can be seen that model, a stress concentration appears ahead of the crack tip, which zone (FPZ), rather than at the back of the crack tip, which is the ca Due to the nonlinearity of the problem, an explicit analytical solution for the displacement field around the crack tip is not available. In order to evaluate the relative error level for different mesh densities, we took the results obtained by the finest meshes with each enrichment scheme as reference exact solutions. The h-version convergence rate of the finite element method was quantified by means of the relative error in the L2-norm, which was calculated by the following equations.
where the superscript ref denotes a reference solution.
In Figure 9, the deformed geometry of the cohesive crack problem is compared with that of the traction-free crack problem when the same load factor is applied. When a cohesive force exists between crack faces, the crack closes smoothly from the physical tip to the fictitious tip. The Contour plots of σ yy for the cohesive crack problem and the tractionfree crack problem are provided in Figure 10. It can be seen that, in the cohesive crack model, a stress concentration appears ahead of the crack tip, which is the fracture process zone (FPZ), rather than at the back of the crack tip, which is the case for the traction-free crack problem. The stress gradient at the crack tip is much smaller compared with the case with the traction-free crack. The stress at the crack tip in Figure 10b is finite, and equal to material tensile strength. This means that the cohesive crack model abandons the unrealistic assumption in the LEFM that the stress at the crack tip is infinite, which can be seen in Figure 10a.  The stress σyy along the axis of symmetry is plotted in Figure 11. The stress profiles obtained by different enrichment schemes are difficult to distinguish from one another. They also show an obvious FPZ ahead of the crack tip, and the stress σyy is equal to the tensile strength at the crack tip. In contrast, a stress singularity appears around the crack tip in the traction-free crack problem. Figure 11. Stress profiles of σyy obtained by different enrichment schemes for the linear softening model. Figure 12 shows the relative error in the L2-norm plotted against the inverse of the element size h, which is taken as the square root of the area of the element. The rates of convergence were obtained by means of polynomial curve fitting of those data points. It is interesting that, as the linear softening model is considered, with the employment of the branch function  The stress σyy along the axis of symmetry is plotted in Figure 11. The stress profiles obtained by different enrichment schemes are difficult to distinguish from one another. They also show an obvious FPZ ahead of the crack tip, and the stress σyy is equal to the tensile strength at the crack tip. In contrast, a stress singularity appears around the crack tip in the traction-free crack problem. Figure 11. Stress profiles of σyy obtained by different enrichment schemes for the linear softening model. Figure 12 shows the relative error in the L2-norm plotted against the inverse of the element size h, which is taken as the square root of the area of the element. The rates of convergence were obtained by means of polynomial curve fitting of those data points. It is interesting that, as the linear softening model is considered, with the employment of the branch function The stress σ yy along the axis of symmetry is plotted in Figure 11. The stress profiles obtained by different enrichment schemes are difficult to distinguish from one another. They also show an obvious FPZ ahead of the crack tip, and the stress σ yy is equal to the tensile strength at the crack tip. In contrast, a stress singularity appears around the crack tip in the traction-free crack problem. The stress σyy along the axis of symmetry is plotted in Figure 11. The stress obtained by different enrichment schemes are difficult to distinguish from one a They also show an obvious FPZ ahead of the crack tip, and the stress σyy is equa tensile strength at the crack tip. In contrast, a stress singularity appears around th tip in the traction-free crack problem.  Figure 12 shows the relative error in the L2-norm plotted against the invers element size h, which is taken as the square root of the area of the element. The convergence were obtained by means of polynomial curve fitting of those data p is interesting that, as the linear softening model is considered, with the employmen branch function (b) Figure 11. Stress profiles of σ yy obtained by different enrichment schemes for the linear softening model. Figure 12 shows the relative error in the L2-norm plotted against the inverse of the element size h, which is taken as the square root of the area of the element. The rates of convergence were obtained by means of polynomial curve fitting of those data points. It is interesting that, as the linear softening model is considered, with the employment of the branch function r sin θ 2 , XFEM-s, XFEM-c1, and XFEM-c2 achieve a better convergence rate of more than 1, compared with the 0.7 obtained by XFEM-h. The employment of the additional branch function r cos θ 2 results in a similar convergence rate but a higher accuracy.
additional branch function 2 cos r  results in a similar convergence ra curacy. As far as accuracy is concerned, the XFEM-s scheme provides less sive crack problems, especially when coarse meshes are used. It can above stress contour plot in Figure 10 that the singularity vanishes a only a finite stress gradient exists. If the parasitic terms resulting from are not corrected, it can reduce the accuracy severely for the case of co the corrected approximation for blending elements as in XFEM-c1 and level is improved by around 2 times, while the convergence rate remai which is similar to the case of strong discontinuities.

Exponential Softening Model
Because of its simplicity, the linear softening model is frequently certain brittle materials, a nonlinear softening model may be more a cohesive traction-displacement relation changes, the stress gradient a differs, and that may affect the convergence rates of enrichment schem specimen was used for the numerical study of convergence rates. The displacement gap relation can be expressed as exp( ) The tensile strength was the same as ft = 3.18 MPa, and the fractur to be smaller than Ef = 12.5 N/m to increase the gradient of the cohesi and exponential softening models used in this study are depicted in seen from the profile that they will result in a similar FPZ length, but dients.
The stress σyy along the symmetric line produced by different enr provided in Figure 13, as well as a comparison of the two cohesive c

1/h Relative Error
XFEM-h, m=0.71 XFEM-s, m=1.14 XFEM-c1, m=1.06 XFEM-c2, m=1.07 As far as accuracy is concerned, the XFEM-s scheme provides less accuracy for cohesive crack problems, especially when coarse meshes are used. It can be seen from the above stress contour plot in Figure 10 that the singularity vanishes at the crack tip, and only a finite stress gradient exists. If the parasitic terms resulting from the branch function are not corrected, it can reduce the accuracy severely for the case of cohesive cracks. With the corrected approximation for blending elements as in XFEM-c1 and XFEM-c2, the error level is improved by around 2 times, while the convergence rate remains almost the same, which is similar to the case of strong discontinuities.

Exponential Softening Model
Because of its simplicity, the linear softening model is frequently used; however, for certain brittle materials, a nonlinear softening model may be more accurate. When the cohesive traction-displacement relation changes, the stress gradient around the crack tip differs, and that may affect the convergence rates of enrichment schemes. The same DCB specimen was used for the numerical study of convergence rates. The cohesive force and displacement gap relation can be expressed as The tensile strength was the same as f t = 3.18 MPa, and the fracture energy was made to be smaller than E f = 12.5 N/m to increase the gradient of the cohesive force. The linear and exponential softening models used in this study are depicted in Figure 8. It can be seen from the profile that they will result in a similar FPZ length, but different stress gradients.
The stress σ yy along the symmetric line produced by different enrichment schemes is provided in Figure 13, as well as a comparison of the two cohesive constitutive models. Likewise, the stress profiles produced by these enrichment schemes are hard to distinguish from one another. It can be seen from Figure 13 b that, when the traction-separation relation changes, although the cohesive force remains equal to the tensile strength at the crack tip, the stress gradient differs in the FPZ, and does not make much difference at the back of the crack tip.  Figure 14 shows the convergence rates of these enrichment schemes when the exponential softening model is inserted. Likewise, the enrichment schemes with tip branch functions exhibit a higher convergence rate. The employment of the tip branch function 2 cos r  increases both the convergence rate and accuracy substantially. However, especially when coarse meshes are used, these enrichment schemes achieve lower accuracy than XFEM-h in terms of error level. In comparison with the cases of the linear constitutive law, the difference between the convergence rates of these enrichment schemes is more pronounced for the cases of the exponential constitutive law.

Mixed-Mode Crack Problem
In this case, a plate with an inclined cohesive crack was analyzed using all four enrichment schemes in order to investigate their convergence properites in depth. The boundary conditions are shown in the sketch in Figure 15. The dimensions of the plate are 200 by 400 mm, with a thickness of 20 mm. The inclined crack is located at [0 150; 100 200]. A uniformly distributed tensile force fext = 1 Mpa was applied on the top edge with the plane stress condition. All the material properties and softening laws were the same as in previous cases.  Figure 14 shows the convergence rates of these enrichment schemes when the exponential softening model is inserted. Likewise, the enrichment schemes with tip branch functions exhibit a higher convergence rate. The employment of the tip branch function r cos θ 2 increases both the convergence rate and accuracy substantially. However, especially when coarse meshes are used, these enrichment schemes achieve lower accuracy than XFEM-h in terms of error level. In comparison with the cases of the linear constitutive law, the difference between the convergence rates of these enrichment schemes is more pronounced for the cases of the exponential constitutive law.  Figure 14 shows the convergence rates of these enrichment s nential softening model is inserted. Likewise, the enrichment sc functions exhibit a higher convergence rate. The employment of 2 cos r  increases both the convergence rate and accuracy substa cially when coarse meshes are used, these enrichment schemes than XFEM-h in terms of error level. In comparison with the cases law, the difference between the convergence rates of these enric pronounced for the cases of the exponential constitutive law.

Mixed-Mode Crack Problem
In this case, a plate with an inclined cohesive crack was ana richment schemes in order to investigate their convergence pr boundary conditions are shown in the sketch in Figure 15. The dim 200 by 400 mm, with a thickness of 20 mm. The inclined crack is lo A uniformly distributed tensile force fext = 1 Mpa was applied on plane stress condition. All the material properties and softening la previous cases.

Mixed-Mode Crack Problem
In this case, a plate with an inclined cohesive crack was analyzed using all four enrichment schemes in order to investigate their convergence properites in depth. The boundary conditions are shown in the sketch in Figure 15 The convergence rates of these enrichment schemes are provided inFig follow similar tendencies. The enrichment schemes with tip branch functions convergence properties, while the corrected approximation and additiona functions can increase the accuracy.  The convergence rates of these enrichment schemes are provided in Figure 16. They follow similar tendencies. The enrichment schemes with tip branch functions have similar convergence properties, while the corrected approximation and additional tip branch functions can increase the accuracy. The convergence rates of these enrichment schemes are provided inFigure 16. They follow similar tendencies. The enrichment schemes with tip branch functions have similar convergence properties, while the corrected approximation and additional tip branch functions can increase the accuracy. Figure 16. Plots of the convergence rate for the mixed-mode crack problem with (a) the linear softening law and (b) the exponential softening law (m is the convergence rate).

Conclusions
The present work focuses on investigating the convergence properties and accuracy of different enrichment schemes in the XFEM for modeling the cohesive crack problem. Four kinds of enrichment schemes for cohesive cracks were manufactured to examine the

Conclusions
The present work focuses on investigating the convergence properties and accuracy of different enrichment schemes in the XFEM for modeling the cohesive crack problem. Four kinds of enrichment schemes for cohesive cracks were manufactured to examine the influences of their variations on convergence performance and accuracy level. The convergence study was conducted on a double-cantilever beam specimen with both the pure mode I problem and the mixed-mode problem, and cases of linear and exponential constitutive laws were considered. On the basis of the simulation results, our main conclusions are as follows.

1.
When both linear and exponential constitutive laws are assumed, the enrichment schemes with tip branch functions achieve a higher convergence rate than XFEM-h; however, they have lower accuracy in terms of the absolute error value.

2.
In the case of the cohesive crack simulation, the corrected approximation for blending elements did not change the convergence rate much, but the error level improved substantially, which is similar to the case of traction-free cracks. The enrichment schemes with tip branch functions have similar convergence properties, while the corrected approximation and additional tip branch functions can increase the accuracy.

3.
As far as accuracy is concerned, the enrichment schemes with tip branch functions perform worse than XFEM-h when coarse meshes are used. If the parasitic terms resulting from the branch function are not corrected, it can reduce the accuracy severely for the simulation of cohesive cracks. 4.
The traction-displacement relation can also affect the convergence properties of these enrichment schemes. In the case of the exponential constitutive law, the difference between the convergence rates of these enrichment schemes is more pronounced.