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

<sup>√</sup>*<sup>r</sup>* cos *<sup>θ</sup>* 2 , <sup>√</sup>*<sup>r</sup>* sin *<sup>θ</sup>* <sup>2</sup> sin *θ*, <sup>√</sup>*<sup>r</sup>* cos *<sup>θ</sup>* <sup>2</sup> sin *θ*}. 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.

$$\phi(r,\theta) = r\sin\frac{\theta}{2}, \text{ or } r^{3/2}\sin\frac{\theta}{2}, \text{ or } r^2\sin\frac{\theta}{2} \tag{13}$$

Others proposed that four branch functions be used to enrich the tip element, which is

$$\{\phi\_a(r,\theta)\} = \begin{Bmatrix} r\sin\frac{\theta}{2} & r\cos\frac{\theta}{2} & r\sin\frac{\theta}{2}\sin\theta & r\cos\frac{\theta}{2}\sin\theta \end{Bmatrix} \tag{14}$$

In this enrichment scheme, *r* sin *<sup>θ</sup>* <sup>2</sup> 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.

**Figure 5.** Tip branch functions for cohesive cracks. (**a**) *φ*(*r*, *θ*) = *r* sin *<sup>θ</sup>* <sup>2</sup> and (**b**) *<sup>φ</sup>*(*r*, *<sup>θ</sup>*) <sup>=</sup> *<sup>r</sup>* cos *<sup>θ</sup>* 2 .

In this scheme, the XFEM approximation of the displacement field can be expressed as

$$\begin{aligned} u^h(\mathbf{x}) &= \sum\_{i \in I} N\_i(\mathbf{x}) \cdot u\_i + \sum\_{j \in J} N\_j(\mathbf{x}) \cdot \left[ H(\mathbf{x}) - H(\mathbf{x}\_j) \right] \cdot a\_j \\ &+ \sum\_{k \in M} N\_k(\mathbf{x}) \cdot \left[ \phi(\mathbf{x}) - \phi(\mathbf{x}\_k) \right] \cdot b\_k \end{aligned} \tag{15}$$

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

#### *2.5. 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 Nk*(*x*) = 1. In addition, the branch function is not a piecewise constant function like the heaviside function, so the parasitic terms resulting from ∑ *k Nk*(*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:

$$\begin{array}{ll} \mu^{h}(\mathbf{x}) = \sum\_{i \in I} N\_{i}(\mathbf{x}) \cdot u\_{i} + \sum\_{j \in J} N\_{j}(\mathbf{x}) \cdot \left[ H(\mathbf{x}) - H(\mathbf{x}\_{j}) \right] \cdot a\_{j} \\ \quad + \sum\_{k \in M \cup L} N\_{k}(\mathbf{x}) \cdot \mathcal{R}(\mathbf{x}) \cdot \sum\_{a} \left[ \phi\_{a}(\mathbf{x}) - \phi\_{a}(\mathbf{x}\_{k}) \right] \cdot b\_{k}^{a} \end{array} \tag{16}$$

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<sup>α</sup> <sup>k</sup>* . *R*(*x*) is a ramp function, which is defined as follows and depicted in Figure 6.

**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 *<sup>θ</sup>* <sup>2</sup> was employed in XFEM-c1, while both *<sup>r</sup>* sin *<sup>θ</sup>* <sup>2</sup> and *<sup>r</sup>* cos *<sup>θ</sup>* <sup>2</sup> were employed in XFEM-c2. It can be seen from their plots in Figure 5 that the branch function *r* cos *<sup>θ</sup>* 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.

#### **3. 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.,

$$K\_{\text{flip}} = 0\tag{18}$$

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.,

$$
\boldsymbol{m}\_{\Gamma}^{\mathsf{T}} \cdot \boldsymbol{\sigma}\_{\mathsf{tip}} \cdot \boldsymbol{n}\_{\Gamma} = \boldsymbol{f}\_{t} \tag{19}
$$

where *σ*tip is the stress at the crack tip and *ft* 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

$$\mathbf{M}^{\mathrm{T}} \cdot \mathbb{C} \cdot \mathcal{B} \cdot \boldsymbol{q} = \mathbb{S} \cdot \boldsymbol{q} = f\_t \tag{20}$$

where *<sup>S</sup>* = *<sup>M</sup>*<sup>T</sup> · *<sup>D</sup>* · *<sup>B</sup>* is the operator by which the stress at the crack tip is calculated, and *M* = *n*<sup>Γ</sup> ⊗ *n*<sup>Γ</sup> 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

$$r = \left\{ \begin{array}{c} K \cdot q - \lambda f^{\text{ext}} - f^{\text{coh}}(q) \\ f\_t - Sq \end{array} \right\} \tag{21}$$

where the independent unknowns are *q* and *λ*. The Jacobian matrix is

$$
\Lambda = \begin{bmatrix}
K - \frac{\partial f^{\text{coh}}(q)}{\partial q} & -f^{\text{ext}} \\
\end{bmatrix}
\tag{22}
$$

where *<sup>∂</sup> <sup>f</sup> coh*(*q*) *<sup>∂</sup><sup>q</sup>* is the additional stiffness term effective on the crack surface in the FPZ, which can be obtained by

$$\frac{\partial f^{\text{coh}}(q)}{\partial q} = -2 \int\_{\Gamma\_c} \frac{\partial \sigma\_{\mathcal{Y}}(\omega)}{\partial \omega} \mathbf{N}^T \cdot \boldsymbol{n} \cdot \boldsymbol{n}^T \cdot \mathbf{N} d\Gamma \tag{23}$$

At the *i*th iteration, the increments in independent variables obtained from Equations (21) and (22) are

$$\left\{ \begin{array}{c} \Delta q\\ \Delta \lambda \end{array} \right\}^{i} = -\left(\Lambda^{i-1}\right)^{-1} \cdot r^{i-1} \tag{24}$$
