**1. 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].

**Citation:** Liu, G.; Guo, J.; Bao, Y. Convergence Investigation of XFEM Enrichment Schemes for Modeling Cohesive Cracks. *Mathematics* **2022**, *10*, 383. https://doi.org/10.3390/ math10030383

Academic Editor: Marcin Kami ´nski

Received: 15 December 2021 Accepted: 18 January 2022 Published: 26 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** The stress profiles around the fracture process zone.

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 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<sup>k</sup>* sin *<sup>θ</sup>* <sup>2</sup> , 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 <sup>√</sup>*<sup>r</sup>* with *<sup>r</sup>*, e.g., *φα*(*r*, *<sup>θ</sup>*) = {*<sup>r</sup>* sin *<sup>θ</sup>* <sup>2</sup> , *<sup>r</sup>* cos *<sup>θ</sup>* 2 , *r* sin *<sup>θ</sup>* <sup>2</sup> sin *<sup>θ</sup>*, *<sup>r</sup>* cos *<sup>θ</sup>* <sup>2</sup> 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.


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.

#### **2. XFEM Formulation for Cohesive Crack Problems**

#### *2.1. 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

$$\nabla \cdot \sigma + b = 0 \text{in } \Omega \tag{1}$$

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

$$\begin{aligned} \mu &= \overline{u} \quad \text{on} \quad \Gamma\_u\\ \sigma \cdot n\_\Gamma &= \overline{t} \quad \text{on} \quad \Gamma\_t\\ \sigma \cdot n\_{\Gamma d} &= f^c \quad \text{on} \quad \Gamma\_c \end{aligned} \tag{2}$$

where *n*<sup>Γ</sup> 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 <sup>c</sup>* is the cohesive traction transferred across the Γ*c*, which is related to the displacement gap ω across the discontinuity according to the stress softening model.

**Figure 2.** A two-dimensional domain containing a cohesive discontinuity Γ*c*.
