# Arbitrary branched and intersecting cracks with the eXtended Finite Element Method

Christophe Daux, Nicolas Moës, John Dolbow, Natarajan Sukumar, Ted Belytschko

► To cite this version:

Christophe Daux, Nicolas Moës, John Dolbow, Natarajan Sukumar, Ted Belytschko. Arbitrary branched and intersecting cracks with the eXtended Finite Element Method. International Journal for Numerical Methods in Engineering, 2000, 48 (12), pp.1741-1760. 10.1002/1097-0207(20000830)48:123.0.CO;2-L. hal-01005274

HAL Id: hal-01005274

https://hal.science/hal-01005274

Submitted on 24 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Arbitrary branched and intersecting cracks with
the extended finite element method

Christophe Daux<sup>‡</sup>, Nicolas Moës<sup>*,†,§</sup>, John Dolbow<sup>¶</sup>,
Natarajan Sukumar<sup>||</sup> and Ted Belytschko<sup>**</sup>

Department of Mechanical and Civil Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, U.S.A.

Extensions of a new technique for the finite element modelling of cracks with multiple branches, multiple holes and cracks emanating from holes are presented. This extended finite element method (X-FEM) allows the representation of crack discontinuities and voids independently of the mesh. A standard displacement-based approximation is enriched by incorporating discontinuous fields through a partition of unity method. A methodology that constructs the enriched approximation based on the interaction of the discontinuous geometric features with the mesh is developed. Computation of the stress intensity factors (SIF) in different examples involving branched and intersecting cracks as well as cracks emanating from holes are presented to demonstrate the accuracy and the robustness of the proposed technique.

KEY WORDS: partition of unity; extended FEM; branched cracks; voids

# 1. INTRODUCTION

Solving crack problems in fracture mechanics is imperative to quantify and predict the behaviour
of cracked structures under service conditions. To this end, the accurate evaluation of fracture

<footnote>*Correspondence to: Nicolas Moës, Department of Mechanical Engineering, Northwestern University, 215 Sheridan Road, Evanston, IL 60208, U.S.A.</footnote>

<footnote>†E-mail: moes@taml.mech.nwu.edu</footnote>

<footnote>‡Research Assistant, Ecole Normale Superieure de Cachan, France</footnote>

<footnote>§Research Associate, Department of Mechanical Engineering</footnote>

<footnote>¶Assistant Professor, Department of Civil & Environmental Engineering, Duke University</footnote>

<footnote>||Post-Doctoral Research Fellow, Department of Civil Engineering</footnote>

<footnote>**Walter P. Murphy, Professor of Mechanical Engineering</footnote>

<footnote>Contract/grant sponsor: Office of Naval Research</footnote>

<footnote>Contract/grant sponsor: Army Research Office</footnote>

<footnote>Contract/grant sponsor: Natural Science Foundation; contract/grant number: CMS-9732319</footnote>
parameters such as the stress intensity factors (SIF) is required for simulation-based life-cycle
design analysis.

Over the past few decades, many numerical methods have been proposed to model crack prob-
lems. Finite element methods with non-singular and singular elements [1, 2] enable the accurate
computation of stress intensity factors. However, these methods require the finite element edges to
coincide with the crack. This often complicates mesh generation, since both the regular geometric
features and the crack must be considered. Some of the other prominent numerical methods for
cracks analysis are the boundary elements method [3], the boundary collocation method [4], the
body force method [5] and the integral equation method [6]. The dislocation method is also often
used for cracks with multiple branches ([7–9]). Recently, meshless methods, and in particular the
element-free Galerkin method, have been applied to two-dimensional crack problems [10].

The extended finite element method (X-FEM) allows for the modelling of arbitrary geometric
features independently of the finite element mesh. The crack modelling technique was presented in
References [11, 12]; see also Reference [13]. This method allows the modelling of crack growth
without remeshing. In this paper, we apply it to voids and to more complex geometries such
as cracks with multiple branches. This facilitates the accurate modelling of interactions between
cracks and holes, or systems with cracks emanating from holes.

By using the notion of partition of unity [14], we enrich the standard approximation with addi-
tional functions. For branched cracks, we use the near-tip asymptotic fields [11], the discontinuous
function for a single crack [12], and a discontinuous function which models the branching. For
voids, we use another discontinuous function which can be multiplied with the previous ones if
the crack is near a hole or emanates from it.

The use of the partition of unity as a means to enrich a finite element space may also be found
in the work of Strouboulis *et al.* [15] for modelling voids and the work of Strouboulis *et al.* [16]
to model reentrant corners and cracks. Also, Duarte *et al.* [17] applied the partition of unity for
the simulation of dynamic crack propagation in three-dimensions. Finally, in the work of Oden
*et al.* [18] a low order set of finite element shape function is enriched with hierarchical polynomials
to create a Cartesian *p* finite element which performs well under distortion.

In the next section, we recall the governing equations for linear elastostatics in their strong and
weak forms. In Section 3, we describe the ‘branched’ or ‘junction’ function and the algorithm
used to enrich the classical displacement field to model a multiple branched crack. In Section
4, the modelling of holes which are not part of the mesh is described. Numerical results are
presented in Section 5 to illustrate the accuracy and the robustness of the X-FEM. Finally, in
Section 6, some concluding remarks on the potential of the extended finite element method are
indicated.

2. FORMULATION

In this section, we briefly review the governing equations for elastostatics. Specifically, we consider
the case when internal boundaries are present.

2.1. *Governing equations*

Consider the domain Ω bounded by Γ. The boundary Γ is composed of Γ_u, Γ_t, Γ_c and Γ_h such that Γ = Γ_u ∪ Γ_t ∪ Γ_c ∪ Γ_h as shown in Figure 1. Prescribed displacements are imposed on Γ_u, while
<figure><img src="image_2.png" /><figcaption>Figure 1. Notations for a body with traction-free cracks and holes subjected to loads and imposed displacements.</figcaption></figure>

prescribed tractions are imposed on $\Gamma_t$. The boundary $\Gamma_c$ consists of the boundaries of all cracks (both faces for each crack) and $\Gamma_h$ gathers all boundaries of holes. The hole and crack faces are assumed to be traction-free. The equilibrium equations and the boundary conditions are

$$ \nabla \cdot \boldsymbol{\sigma} + \mathbf{b} = 0 \quad \text{in } \Omega \qquad (1a) $$

$$ \mathbf{u} = \bar{\mathbf{u}} \quad \text{on } \Gamma_u \qquad (1b) $$

$$ \boldsymbol{\sigma} \cdot \mathbf{n} = \bar{\mathbf{t}} \quad \text{on } \Gamma_t \qquad (1c) $$

$$ \boldsymbol{\sigma} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_c \qquad (1d) $$

$$ \boldsymbol{\sigma} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_h \qquad (1e) $$

where $\mathbf{n}$ is the unit outward normal, $\boldsymbol{\sigma}$ the Cauchy stress, and $\mathbf{b}$ is the body force per unit volume.

We consider small strains and displacements. The kinematics equations therefore consist of the strain-displacement relation

$$ \varepsilon = \varepsilon(\mathbf{u}) = \nabla_s \mathbf{u} \qquad (2) $$

where $\nabla_s$ is the symmetric part of the gradient operator. The constitutive relation is given by Hooke's law

$$ \boldsymbol{\sigma} = C\varepsilon \qquad (3) $$

where $C$ is Hooke's tensor.

## 2.2. Weak form

The space of trial functions is defined by

$$ \mathcal{U} = \{ \mathbf{u} \in \mathcal{V} : \mathbf{u} = \bar{\mathbf{u}} \text{ on } \Gamma_u, \mathbf{u} \text{ discontinuous on } \Gamma_c \cup \Gamma_h \} \qquad (4) $$
where the space $\mathcal{V}$ is related to the regularity of the solution. The test function space is defined similarly as

$$ \mathcal{U}_0 = \{ \mathbf{v} \in \mathcal{V} : \mathbf{v} = 0 \text{ on } \Gamma_u, \mathbf{v} \text{ discontinuous on } \Gamma_c \cup \Gamma_h \} \quad (5) $$

The weak form of the equilibrium equations is given by

$$ \int_{\Omega} \boldsymbol{\sigma}(\mathbf{u}):\boldsymbol{\varepsilon}(\mathbf{v}) \, d\Omega = \int_{\Omega} \mathbf{b} \cdot \mathbf{v} \, d\Omega + \int_{\Gamma_t} \bar{\mathbf{t}} \cdot \mathbf{v} \, d\Gamma \quad \forall \mathbf{v} \in \mathcal{U}_0 \quad (6) $$

Using the constitutive relation, Equation (3), and the strain-displacement relation, Equation (2), the weak form can be stated as: find $\mathbf{u} \in \mathcal{U}$ such that

$$ \int_{\Omega} \boldsymbol{\varepsilon}(\mathbf{u}) : \mathbf{C} : \boldsymbol{\varepsilon}(\mathbf{v}) \, d\Omega = \int_{\Omega} \mathbf{b} \cdot \mathbf{v} \, d\Omega + \int_{\Gamma_t} \bar{\mathbf{t}} \cdot \mathbf{v} \, d\Gamma \quad \forall \mathbf{v} \in \mathcal{U}_0 \quad (7) $$

It is shown in Reference [11] that the above is equivalent to the strong form (1), including the traction-free conditions on any lines of discontinuity: the two crack faces and the holes, i.e. $\Gamma_c \cup \Gamma_h$.

### 3. MODELLING A BRANCHED CRACK WITH THE X-FEM

It has been shown in Reference [12] that a crack can be modelled independently of the mesh by enriching the approximation by step functions and asymptotic near-tip fields. The finite element approximation for a single crack in a two-dimensional body can be written as

$$ \begin{aligned} \mathbf{u}^h(\mathbf{x}) = & \sum_{i \in I} \mathbf{u}_i \phi_i(\mathbf{x}) + \sum_{i \in L} \mathbf{a}_i \phi_i(\mathbf{x}) H(\mathbf{x}) \\ & + \sum_{i \in K_1} \phi_i(\mathbf{x}) \left( \sum_{l=1}^4 \mathbf{b}_{i,1}^l F_1^l(\mathbf{x}) \right) + \sum_{i \in K_2} \phi_i(\mathbf{x}) \left( \sum_{l=1}^4 \mathbf{b}_{i,2}^l F_2^l(\mathbf{x}) \right) \end{aligned} \quad (8) $$

where $I$ is the set of all nodes in the mesh; $\mathbf{u}_i$ is the classical (vectorial) degree of freedom at node $i$; $\phi_i$ is the shape function associated with node $i$. Each shape function $\phi_i$ has compact support $\omega_i$ given by the union of the elements connected to node $i$; $L \subset I$ is the subset of nodes that are enriched for the crack discontinuity and $\mathbf{a}_i$ are the corresponding additional degrees of freedom; the nodes in $L$ are such that their support (we mean the support of the nodal shape function) intersects the crack but do not contain any of its crack tips; and $K_1 \subset I$ and $K_2 \subset I$ are the subset of nodes that are enriched for the first and second crack tip, respectively. The corresponding additional degrees of freedom are $\mathbf{b}_{i,1}^l$ and $\mathbf{b}_{i,2}^l$, $l = 1, \dots, 4$, for the first and second crack tip, respectively; the nodes in $K_1$ ($K_2$) are such that their support contain the first (second) crack tip.
<figure><img src="image_3.png" /><figcaption>Figure 2. Normal and tangential co-ordinates for a smooth crack, (a), and a kinked crack, (b). $\mathbf{x}^*$ is the closest point to $\mathbf{x}$ on the crack. The value of the function $H$ is -1 at $\mathbf{x}$ for (a) and (b).</figcaption></figure>

The near-tip functions $F_1^l(\mathbf{x})$, $l=1,...,4$, are given by

$$ \{F_1^l(\mathbf{x})\} = \left\{ \sqrt{r} \sin\left(\frac{\theta}{2}\right), \sqrt{r} \cos\left(\frac{\theta}{2}\right), \sqrt{r} \sin\left(\frac{\theta}{2}\right) \sin(\theta), \sqrt{r} \cos\left(\frac{\theta}{2}\right) \sin(\theta) \right\} \quad (9) $$

where $(r, \theta)$ are the local polar co-ordinates at the first crack tip with $\theta=0$ coinciding with the tangent to the crack at the tip. Similarly, the near-tip functions $F_2^l(\mathbf{x})$ are also given by (9) but the local polar co-ordinates being now defined at the second crack tip.

The function $H(\mathbf{x})$ is a discontinuous function across the crack surface and is constant on each side of the crack: +1 on one side of the crack and -1 on the other. More precisely, the function $H(\mathbf{x})$ is defined as follows. The crack is considered to be a curve parametrized by the curvilinear co-ordinate $s$, as in Figure 2. The origin of the curve is taken to coincide with one of the crack tips. Given a point $\mathbf{x}$ in the domain, we denote by $\mathbf{x}^*$ the closest point on the crack to $\mathbf{x}$. At $\mathbf{x}^*$, we construct the tangent and normal vector to the curve, $\mathbf{e}_s$ and $\mathbf{e}_n$, with the orientation of $\mathbf{e}_n$ taken such that $\mathbf{e}_s \times \mathbf{e}_n = \mathbf{e}_z$ where the unit vector $\mathbf{e}_z$ points out of the page. The function $H(\mathbf{x})$ is then given by the sign of the scalar product $(\mathbf{x} - \mathbf{x}^*) \cdot \mathbf{e}_n$. In the case of a kinked crack as shown in Figure 2(b), where no unique normal exists but a cone of normals is defined at $\mathbf{x}^*$, $H(\mathbf{x})=1$ if the vector $(\mathbf{x} - \mathbf{x}^*)$ belongs to the cone of normals at $\mathbf{x}^*$ and -1 otherwise.

### 3.1. Branched cracks

In this section, we introduce the discontinuous functions for modelling branched cracks. We first consider a simple case of a branched crack modelled by four elements as shown in Figure 3. We wish to illustrate how an equivalent discrete space can be constructed with the mesh shown in Figure 4 by enriching with the discontinuous function $H(\mathbf{x})$ and another function which takes into account the junction between the three branches. The finite element approximation associated with the mesh in Figure 3 is

$$ \mathbf{u}^h(\mathbf{x}) = \sum_{i \in I} \mathbf{u}_i \phi_i(\mathbf{x}), \quad I = \{1, \dots, 17\} \qquad (10) $$
<figure><img src="image_4.png" /><figcaption>Figure 3. FEM mesh.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. X-FEM mesh.</figcaption></figure>

where $\mathbf{u}_i$ is the displacement at node $i$ and $\phi_i$ is the bilinear shape function associated with node $i$.
Let us introduce the nodal vectorial variables $\alpha_8, \beta_8, \alpha_{15}, \beta_{15}, \alpha_{18}, \beta_{18}, \alpha_{12}, \beta_{12}$ and $\gamma_{12}$ such that

$$ 
\begin{aligned}
\alpha_8 &= \frac{\mathbf{u}_6 + \mathbf{u}_7}{2}, & \beta_8 &= \frac{\mathbf{u}_6 - \mathbf{u}_7}{2} \\
\alpha_{15} &= \frac{\mathbf{u}_{13} + \mathbf{u}_{14}}{2}, & \beta_{15} &= \frac{\mathbf{u}_{14} - \mathbf{u}_{13}}{2} \\
\alpha_{18} &= \frac{\mathbf{u}_{16} + \mathbf{u}_{17}}{2}, & \beta_{18} &= \frac{\mathbf{u}_{17} - \mathbf{u}_{16}}{2} \\
\alpha_{12} &= \frac{\mathbf{u}_{11} + (\mathbf{u}_9 + \mathbf{u}_{10})}{2}, & \gamma_{12} &= \frac{\mathbf{u}_9 - \mathbf{u}_{10}}{2} \\
\beta_{12} &= \frac{\mathbf{u}_{11} - (\mathbf{u}_9 + \mathbf{u}_{10})}{2},
\end{aligned} 
\quad (11) $$

Inverting the previous system, we obtain the following equations:

$$ 
\begin{aligned}
\mathbf{u}_6 &= \alpha_8 + \beta_8, & \mathbf{u}_7 &= \alpha_8 - \beta_8 \\
\mathbf{u}_9 &= \alpha_{12} - \beta_{12} + \gamma_{12}, & \mathbf{u}_{10} &= \alpha_{12} - \beta_{12} - \gamma_{12} \\
\mathbf{u}_{11} &= \alpha_{12} + \beta_{12}, & \mathbf{u}_{13} &= \alpha_{15} - \beta_{15} \\
\mathbf{u}_{14} &= \alpha_{15} + \beta_{15}, & \mathbf{u}_{16} &= \alpha_{18} - \beta_{18} \\
\mathbf{u}_{17} &= \alpha_{18} + \beta_{18}
\end{aligned} 
\quad (12) $$

Now by expressing $\mathbf{u}_i$ in terms of $\alpha_j, \beta_j, \gamma_j$ in Equation 10, we obtain

$$ 
\begin{aligned}
\mathbf{u}^h = & \sum_{i=1}^{5} \mathbf{u}_i \phi_i + \alpha_8 (\phi_6 + \phi_7) + \beta_8 H_{II}(\mathbf{x})(\phi_6 + \phi_7) \\
& + \alpha_{12} (\phi_9 + \phi_{10} + \phi_{11}) + \beta_{12} H_1(\mathbf{x})(\phi_9 + \phi_{10} + \phi_{11}) \\
& + \gamma_{12} J(\mathbf{x})(\phi_9 + \phi_{10}) \\
& + \alpha_{15} (\phi_{13} + \phi_{14}) + \beta_{15} H_1(\mathbf{x})(\phi_{13} + \phi_{14}) \\
& + \alpha_{18} (\phi_{16} + \phi_{17}) + \beta_{18} H_1(\mathbf{x})(\phi_{16} + \phi_{17})
\end{aligned} 
\quad (13) $$

where the functions $H_1(\mathbf{x})$ and $H_{II}(\mathbf{x})$ and $J(\mathbf{x})$ are defined in Figure 5.
<figure><img src="image_6.png" /><figcaption>Figure 5. Definition of the discontinuous enrichment functions $H_I(\mathbf{x})$, $H_{II}(\mathbf{x})$ and $J(\mathbf{x})$.</figcaption></figure>

We next relate the functions $H_I(\mathbf{x})$ and $H_{II}(\mathbf{x})$ to the function $H(\mathbf{x})$ defined in the previous section. Consider the branched crack as consisting of two cracks denoted by cracks I and II in Figure 5. Deciding that crack I is parametrized so that we traverse it from nodes 15 to 18 and that crack II is parametrized so that we traverse it from nodes 8 to 12, the functions $H_I(\mathbf{x})$ and $H_{II}(\mathbf{x})$ are then the function $H(\mathbf{x})$ associated with cracks I and II, respectively.

The function $J(\mathbf{x})$ is referred to as a discontinuous 'junction' function. It may be expressed in terms of the functions $H_I(\mathbf{x})$ and $H_{II}(\mathbf{x})$ as

$$J(\mathbf{x}) = \begin{cases} H_{II}(\mathbf{x}) & \text{for } H_I(\mathbf{x}) < 0 \\ 0 & \text{for } H_I(\mathbf{x}) > 0 \end{cases}$$

If we now consider the mesh in Figure 4, $\phi_6 + \phi_7$, $\phi_9 + \phi_{10} + \phi_{11}$, $\phi_{13} + \phi_{14}$, $\phi_{16} + \phi_{17}$ can be replaced by $\phi_8$, $\phi_{12}$, $\phi_{15}$, $\phi_{18}$, respectively, and $\alpha_8$, $\alpha_{12}$, $\alpha_{15}$, $\alpha_{18}$ by $\mathbf{u}_8$, $\mathbf{u}_{12}$, $\mathbf{u}_{15}$, $\mathbf{u}_{18}$. The finite element approximation now reads

$$\begin{aligned} \mathbf{u}^h = & \sum_{i \in I'} \mathbf{u}_i \phi_i \\ & + \beta_{12} \phi_{12} H_I(\mathbf{x}) + \beta_{15} \phi_{15} H_I(\mathbf{x}) + \beta_{18} \phi_{18} H_I(\mathbf{x}) \\ & + \beta_8 \phi_8 H_{II}(\mathbf{x}) \\ & + \gamma_{12} \phi_{12} J(\mathbf{x}) \end{aligned} \quad (14)$$
<figure><img src="image_7.png" /><figcaption>Figure 6. Example of an enrichment for a branched crack.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Crossing cracks, — = main crack 1, --- = main crack 2.</figcaption></figure>

where $I' = \{1, 2, 3, 4, 5, 8, 12, 15, 18\}$. The first term on the right-hand side represents the classical finite element approximation associated with the mesh and the next four represent the discontinuous enrichment for cracks I and II taken separately, finally, the last term takes into account the junction. This shows that the discrete spaces obtained by FEM and X-FEM are equivalent for a branched crack aligned with the mesh.

We now extend the previous derivation to the case where the branched crack geometry is represented independently of the mesh. Consider the branched crack shown in Figure 6. Let us view this crack as a main crack (solid line) and a branched secondary crack (dashed line). The enrichment proceeds as follows. The main crack is enriched as if the secondary crack were absent. This enrichment involves the tip enrichments for the crack tips B and C, denoted by circles and the discontinuous interior enrichment denoted by squares. The secondary crack is enriched almost as if the main crack was absent: the tip D is enriched (plus) and all nodes for which the support intersects the secondary crack but do not contain points A or D are part of the discontinuous interior enrichment (diamonds). Finally, we need to enrich for the junction: the junction function is added to all nodes for which the support of the finite element shape functions contain the junction point A (crosses).

## 3.2. Cracks with multiple branches

Having introduced the junction function $J(x)$, we now describe how we can extend it to model multiple branched cracks. In the case of crossing cracks, the first idea that comes to mind is to consider the two cracks as independent cracks as shown in Figure 7. But it is easy to show that this approach is incorrect. For instance, consider a cross aligned with a mesh. If we do not use the function $J(x)$, we obtain only $3n_D$ degrees of freedom to represent the displacement at the centre node where $n_D = 2$ is the number of displacement components (x- and y-). Indeed, we have the $n_D$ standard finite element degrees of freedom and $n_D$ degrees of freedom per crack for the jump function leading to a total of $3n_D$ degrees of freedom. The finite element description involves however $4n_D$ degrees of freedom at the intersection of the two cracks. Consider now a
<figure><img src="image_9.png" /><figcaption>Figure 8. A crack with multiple branches, --- = main crack, ---- = secondary 1, --- = secondary 2.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Zoom on four finite elements near a branching point to show the subpolygons for the quadrature of the weak form.</figcaption></figure>

main crack and two cracks which join the main crack as shown in Figure 8. In this case, we obtain for the centre node $4n_D$ degrees of freedom: the $n_D$ classical degree of freedom, the $n_D$ degrees of freedom associated to the main crack and $n_D$ degrees of freedom for each junction leading to the correct $4n_D$ number of degrees of freedom.

In essence, using the X-FEM, we recast the independent unknowns in the FEM discretization in a different form. For multiple branched crack, we must have the same number of degrees of freedom per node as twice ($n_D = 2$) the number of pieces its support are cut into. Thus, the procedure for modelling a multiple branched crack consists of discretizing the branched crack into a main crack and several cracks joining the main crack. This approach leads to the right number of degrees of freedom per node: $n_D$ per crack which joins the main crack, $n_D$ associated with the main crack, and the classical $n_D$ degrees of freedom.

Consider now a main crack to which several secondary cracks are connected. We assume that the main and the secondary cracks do not have any branches, although it is easy to generalize to this case. Let $N_c$ be the number of cracks (the main crack plus the number of secondary cracks), $N_t$ be the number of crack tips and $N_x$ be the number of junctions ($N_x = N_c - 1$). The displacement approximation is

$$ \mathbf{u}^h(\mathbf{x}) = \sum_{i \in I} \mathbf{u}_i \phi_i(\mathbf{x}) + \sum_{j=1}^{N_c} \sum_{i \in L_j} \mathbf{a}_{i,j} \phi_i(\mathbf{x}) H_j(\mathbf{x}) + \sum_{j=1}^{N_t} \sum_{i \in K_j} \phi_i(\mathbf{x}) \left( \sum_{l=1}^{4} \mathbf{b}_{i,j}^l F_j^l(\mathbf{x}) \right) \\ + \sum_{j=1}^{N_x} \sum_{i \in J_j} \mathbf{c}_{i,j} \phi_i(\mathbf{x}) J_j(\mathbf{x}) \qquad (15) $$
where the notations not defined so far are:

(1) $L_j \subset I$ is the subset of nodes to enrich for the $j$th crack discontinuity and $a_{i,j}$ are the corresponding additional degrees of freedom; the nodes in $L_j$ are such that their support intersects the $j$th crack but do not contain any of its two extremities;

(2) $K_j \subset I$ is the subset of nodes to enrich for the $j$th crack tip and $b_{i,j}^l$, $l=1, \dots, 4$, are the corresponding additional degrees of freedom; the nodes in $K_j$ are such that their support contain the $j$th crack tip;

(3) $J_j \subset I$ is the subset of nodes to enrich for the $j$th junction and $c_{i,j}$ are the corresponding additional degrees of freedom; the nodes in $J_j$ are such that their support contain the $j$th junction.

### 3.3. Numerical integration of the weak form

For the elements cut by the crack, it is necessary to modify the element quadrature to accurately evaluate the contribution to the weak form for both sides of the discontinuity. As the crack is allowed to be arbitrarily oriented in an element, the use of standard Gauss quadrature may not adequately integrate the discontinuous field. If the integration of the discontinuous enrichment is indistinguishable from that of a constant function, the system of equations may be rank deficient. In this section, we present the modifications made to the numerical integration scheme for elements cut by a crack which are extensions of Moës et al. [12].

The discrete weak form is normally constructed with a loop over all elements, as the domain is approximated by

$$ \Omega = \cup \Omega_e \qquad (16) $$

where $\Omega_e$ is the element subdomain. For elements cut by a crack, we divide the element into triangular subdomains $\Omega_s$ with boundaries aligned with the crack geometry

$$ \Omega_e = \cup \Omega_s \qquad (17) $$

It is emphasized that the subdomains $\Omega_s$ are only necessary for quadrature; no additional degrees of freedom are associated with them. An example of subdomains creation near a branching point is shown in Figure 9.

Nodes are removed from the sets $L$ or $J$ if their support is almost entirely on one side of the crack. This is implemented by

$$ \text{node belongs to set } L \text{ if } \frac{\min(A_{-1}, A_{+1})}{A_{\omega}} > \varepsilon \qquad (18) $$

$$ \text{node belongs to set } J \text{ if } \frac{\min(A_0, A_{-1}, A_{+1})}{A_{\omega}} > \varepsilon \qquad (19) $$

where $A_0, A_{+1}$ and $A_{-1}$ are shown in Figure 10, $A_\omega$ is the area of the support and, $\varepsilon$ a specified tolerance; $10^{-4}$ is typically used. This prevents rank deficiency in the discrete system of equations.
<figure><img src="image_11.png" /><figcaption>Figure 10. Definition of the areas associated with the support of a nodal-shape function.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Enrichment for two holes.</figcaption></figure>

## 4. MODELLING HOLES WITH THE X-FEM

We show here how to model holes which are not defined by the mesh. The idea is to define an enrichment function $V(\mathbf{x})$ which is zero in the holes and one in the body

$$V(\mathbf{x}) = \begin{cases} 1 & \text{if } \mathbf{x} \in \Omega \\ 0 & \text{if } \mathbf{x} \notin \Omega \end{cases} \qquad (20)$$

where $\Omega$ is the domain occupied by the body. If the support of a nodal shape function intersects a hole, the nodal shape function is multiplied by the $V(\mathbf{x})$ function so that the support size is reduced to its material fraction. Also, the nodal degrees of freedom for which the support are totally in the void are eliminated (or set to zero depending on the implementation). Figure 11 illustrates the enrichment strategy for a body containing two holes.

To model a crack emanating from a hole, we only have to replace the classical shape functions $\phi_i(\mathbf{x})$ by $\psi_i(\mathbf{x}) = V(\mathbf{x})\phi_i(\mathbf{x})$ in Equation (16) and to set to zero all the degrees of freedom attached to a node whose support is entirely in the hole. In practice, we do not replace the classical shape
<figure><img src="image_13.png" /><figcaption>Figure 12. Symmetrical branched crack in a plate under uniaxial traction.</figcaption></figure>

functions $\phi_i(\mathbf{x})$ by $\psi_i(\mathbf{x}) = V(\mathbf{x})\phi_i(\mathbf{x})$, we simply omit the integral in the holes. Furthermore, in Equations (19) and (20), we replace the area $A_\omega$ by $A_\omega - A_{\text{hole}}$ where $A_{\text{hole}}$ is the area of the nodal support in the hole. Note that since the function is discontinuous on the boundary of the holes, zero tractions are automatically implied by the weak form.

## 5. NUMERICAL EXAMPLES

In this section, we apply the X-FEM to various crack problems. A branched crack in a finite plate is studied. We also study the influence of the relative size of the plate on a cross and a star-shaped crack, and we consider the case of two cracks emanating from a hole in a finite plate. The stress intensity factors are computed in all cases, and the accuracy and the robustness of the method are demonstrated.

### 5.1. Branched crack in an infinite plate

We examine the example of a symmetric branched crack in a plate of width $2w$ and height $2H$ subjected to uniaxial tension perpendicular to the main crack. In all computations, $w=20$, $H=16$ and the length of the main crack, $a$, is unity. These dimensions are chosen to ‘model’ an infinite plate.

We use the notations shown in Figure 12. The normalized stress intensity factors for an infinite plate at tips A and B

$$F_{I}^{A} = K_{I}^{A} / \sigma \sqrt{\pi c}, \quad F_{I}^{B} = K_{I}^{B} / \sigma \sqrt{\pi c}, \quad F_{II}^{B} = K_{II}^{B} / \sigma \sqrt{\pi c}$$

are given in Reference [8] and used as a benchmark.

We initially consider a geometry where $b/a=1$ and $\theta = \pi/4$. The average size of the finite elements in the vicinity of the crack is $h$. The stress intensity factors are computed for several ratios $h/a$ and compared to those given by Chen and Hasebe [8]. They are computed using domain
<table><caption>Table I. Normalized $K_i^A$, $K_i^B$, $K_{II}^B$ values for various ratios $h/a$.</caption><thead><tr><th>$h/a$</th><th>0.40</th><th>0.30</th><th>0.22</th><th>0.18</th><th>0.16</th><th>0.14</th><th>0.12</th><th>0.10</th><th>0.05</th><th>*</th></tr></thead><tbody><tr><th>$F_i^A$</th><td>0.963</td><td>1.009</td><td>1.027</td><td>1.056</td><td>1.038</td><td>1.042</td><td>1.045</td><td>1.045</td><td>1.044</td><td>1.044</td></tr><tr><th>$F_i^B$</th><td>0.460</td><td>0.468</td><td>0.498</td><td>0.493</td><td>0.493</td><td>0.494</td><td>0.493</td><td>0.495</td><td>0.496</td><td>0.495</td></tr><tr><th>$F_{II}^B$</th><td>0.458</td><td>0.464</td><td>0.501</td><td>0.506</td><td>0.503</td><td>0.505</td><td>0.504</td><td>0.507</td><td>0.508</td><td>0.506</td></tr></tbody></table>

forms of the interaction integrals [19, 20]. The X-FEM results, shown in Table I, are in good agreement with the reference solution as the mesh is refined.

The SIFs obtained by the X-FEM for different values of $b/a$ and $\theta$ are shown in Table II. All results are computed using the mesh shown in Figure 13 in which the crack is not explicitly model. Our results are comparable to those found by Kitagawa and Yuuki [21] and Chen and Hasebe [8]. The variation of $F_i^A$ is plotted in Figure 14. For the ratio $b/a=0.2$ and $\theta=15^\circ$, the mesh was not fine enough to compute the stress intensity factors: the J-integral domain contained several crack tips.

We note a negative value of $K_1$ in Table II when $b/a=2.0$ and $\theta=75^\circ$ suggesting the need of taking contact into account in this case. In Reference [22], it is mentioned that at the branching point a stress field singularity exists which is weaker than the crack tips singularity. This singularity has not been modelled here but can easily be included.

## 5.2. Branched crack in a finite plate

We now analyse the influence of the size of the plate on the stress intensity factors by varying the ratio $a/w$. The crack geometry, Figure 12, is given by $\theta=45^\circ$ and $b/a=1$. The stress intensity factors obtained are given in Table III and plotted in Figure 15 to show their deviations from the case of an infinite plate. The results presented by Theocaris et al. [23] do not agree with the present formulation as the ratio $a/w$ tends to zero.

## 5.3. Cross crack in a finite plate

In the problem of a cross crack shown in Figure 16, a square plate with bi-axial loading is analysed. The SIFs have been computed for the ratio $a/w$ varying from 0.1 to 0.9. The mesh used was obtained by the mesh generator gmsh [24]. It is an unstructured mesh which is refined in the vicinity of the cracks (but not conforming to it). In Table IV, we give the normalized SIF $F_1$, where $F_1=K_1/\sigma\sqrt{\pi a}$. The SIFs obtained by the X-FEM are close to those computed by Cheung et al. [25] provided the ratio $a/w$ is less than 0.6. Furthermore, the results we have obtained when $a/w$ is rather small can be compared with those obtained by an analytical method for a cross crack in an infinite plate [6], $F_1=0.8636$. The $F_1$ value obtained for $a/w=0.1$ is 0.8653 which is very close to the above reference value.

## 5.4. Star-shaped crack in a finite plate

We model a star-shaped crack in a finite plate under bi-axial loading as shown in Figure 17 and we analyse the influence of the size of the plate on the stress intensity factors. The normalized

<footnote>*Reference solution by Chen and Hasebe [8].</footnote>
<table><caption>Table II. Normalized SIFs for various ratios b/a and angles θ.</caption><thead><tr><th rowspan="2">b/a</th><th rowspan="2">θ</th><th colspan="2">15°</th><th colspan="2">30°</th><th colspan="2">45°</th><th colspan="2">60°</th><th colspan="2">75°</th></tr><tr><th>X-FEM</th><th>*</th><th>X-FEM</th><th>*</th><th>X-FEM</th><th>*</th><th>X-FEM</th><th>*</th><th>X-FEM</th><th>*</th></tr></thead><tbody><tr><td rowspan="3">2.0</td><td>F<sup>A</sup><sub>I</sub></td><td>1.023</td><td>—</td><td>1.044</td><td>—</td><td>1.078</td><td>—</td><td>1.135</td><td>—</td><td>1.133</td><td>—</td></tr><tr><td>F<sup>B</sup><sub>I</sub></td><td>0.749</td><td>—</td><td>0.663</td><td>—</td><td>0.494</td><td>—</td><td>0.282</td><td>—</td><td>-0.036</td><td>—</td></tr><tr><td>F<sup>B</sup><sub>II</sub></td><td>0.138</td><td>—</td><td>0.364</td><td>—</td><td>0.516</td><td>—</td><td>0.562</td><td>—</td><td>0.203</td><td>—</td></tr><tr><td rowspan="3">1.5</td><td>F<sup>A</sup><sub>I</sub></td><td>1.020</td><td>—</td><td>1.035</td><td>—</td><td>1.059</td><td>—</td><td>1.103</td><td>—</td><td>1.180</td><td>—</td></tr><tr><td>F<sup>B</sup><sub>I</sub></td><td>0.744</td><td>—</td><td>0.663</td><td>—</td><td>0.495</td><td>—</td><td>0.283</td><td>—</td><td>0.075</td><td>—</td></tr><tr><td>F<sup>B</sup><sub>II</sub></td><td>0.130</td><td>—</td><td>0.357</td><td>—</td><td>0.515</td><td>—</td><td>0.568</td><td>—</td><td>0.500</td><td>—</td></tr><tr><td rowspan="3">1.0</td><td>F<sup>A</sup><sub>I</sub></td><td>1.016</td><td>1.018</td><td>1.027</td><td>1.030</td><td>1.045</td><td>1.044</td><td>1.066</td><td>1.069</td><td>1.118</td><td>1.117</td></tr><tr><td>F<sup>B</sup><sub>I</sub></td><td>0.750</td><td>0.737</td><td>0.659</td><td>0.658</td><td>0.493</td><td>0.495</td><td>0.281</td><td>0.281</td><td>0.061</td><td>0.061</td></tr><tr><td>F<sup>B</sup><sub>II</sub></td><td>0.123</td><td>0.114</td><td>0.344</td><td>0.343</td><td>0.504</td><td>0.506</td><td>0.576</td><td>0.577</td><td>0.535</td><td>0.541</td></tr><tr><td rowspan="3">0.8</td><td>F<sup>A</sup><sub>I</sub></td><td>1.015</td><td>1.016</td><td>1.028</td><td>1.026</td><td>1.036</td><td>1.036</td><td>1.053</td><td>1.053</td><td>1.086</td><td>1.087</td></tr><tr><td>F<sup>B</sup><sub>I</sub></td><td>0.736</td><td>0.735</td><td>0.660</td><td>0.658</td><td>0.494</td><td>0.495</td><td>0.281</td><td>0.281</td><td>0.057</td><td>0.056</td></tr><tr><td>F<sup>B</sup><sub>II</sub></td><td>0.107</td><td>0.107</td><td>0.336</td><td>0.333</td><td>0.497</td><td>0.498</td><td>0.573</td><td>0.567</td><td>0.546</td><td>0.551</td></tr><tr><td rowspan="3">0.6</td><td>F<sup>A</sup><sub>I</sub></td><td>1.014</td><td>1.014</td><td>1.022</td><td>1.023</td><td>1.030</td><td>1.029</td><td>1.038</td><td>1.039</td><td>1.057</td><td>1.060</td></tr><tr><td>F<sup>B</sup><sub>I</sub></td><td>0.735</td><td>0.732</td><td>0.659</td><td>0.657</td><td>0.498</td><td>0.497</td><td>0.284</td><td>0.284</td><td>0.058</td><td>0.054</td></tr><tr><td>F<sup>B</sup><sub>II</sub></td><td>0.097</td><td>0.095</td><td>0.322</td><td>0.319</td><td>0.486</td><td>0.485</td><td>0.567</td><td>0.568</td><td>0.542</td><td>0.554</td></tr><tr><td rowspan="3">0.4</td><td>F<sup>A</sup><sub>I</sub></td><td>1.012</td><td>1.011</td><td>1.020</td><td>1.018</td><td>1.022</td><td>1.023</td><td>1.028</td><td>1.028</td><td>1.035</td><td>1.037</td></tr><tr><td>F<sup>B</sup><sub>I</sub></td><td>0.752</td><td>0.729</td><td>0.659</td><td>0.659</td><td>0.502</td><td>0.504</td><td>0.295</td><td>0.295</td><td>0.067</td><td>0.066</td></tr><tr><td>F<sup>B</sup><sub>II</sub></td><td>0.096</td><td>0.078</td><td>0.299</td><td>0.295</td><td>0.459</td><td>0.460</td><td>0.546</td><td>0.546</td><td>0.542</td><td>0.542</td></tr><tr><td rowspan="3">0.2</td><td>F<sup>A</sup><sub>I</sub></td><td>—</td><td>—</td><td>1.014</td><td>1.012</td><td>1.019</td><td>1.015</td><td>1.018</td><td>1.018</td><td>1.022</td><td>1.021</td></tr><tr><td>F<sup>B</sup><sub>I</sub></td><td>—</td><td>—</td><td>0.674</td><td>0.667</td><td>0.532</td><td>0.528</td><td>0.339</td><td>0.335</td><td>0.115</td><td>0.122</td></tr><tr><td>F<sup>B</sup><sub>II</sub></td><td>—</td><td>—</td><td>0.254</td><td>0.247</td><td>0.412</td><td>0.405</td><td>0.493</td><td>0.492</td><td>0.498</td><td>0.497</td></tr></tbody></table>

SIFs are defined as

$$F_{I}^{A} = K_{I}^{A} / \sigma \sqrt{\pi a}, \quad F_{I}^{B} = K_{I}^{B} / \sigma \sqrt{\pi a}, \quad F_{II}^{B} = K_{II}^{B} / \sigma \sqrt{\pi a}$$

We have computed the normalized SIFs for ratios of a/w from 0.1 to 0.9 using two different meshes obtained by the mesh generator gmsh. One is an unstructured mesh which is refined in the vicinity of the cracks (but not conforming to it), and the second is a uniform mesh (40 × 40 quadrilateral elements). The two meshes have about the same number of nodes. The results are given in Table V and demonstrate the accuracy and robustness of the X-FEM: the stress intensity

<footnote>*Reference solution by Chen and Hasebe [8].</footnote>
<figure><img src="image_14.png" /><figcaption>Figure 13. (a) Mesh used to compute the SIFs in an ‘infinite’ plate (1218 nodes, $h/a = 0.12$), (b) zoom in the region of the branched crack.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14. Normalized $K_i^A$ values for a branched crack in a finite plate for various ratios $b/a$ and angles $\theta$.</figcaption></figure>

factors computed with the uniform mesh are comparable to those obtained with the unstructured mesh and both are in good agreement with the values determined by Cheung *et al.* [25]. Furthermore, when the ratio $a/w$ tends to zero, the computed results converge to the analytical value for an infinite plate calculated by Ouchterlong [26] ($F_1 = 0.7454$).
<table><caption>Table III. Normalized $K_1^A$, $K_1^B$, $K_{II}^B$ values for a branched crack in a finite plate for various ratios $a/w$.</caption><thead><tr><th>a/w</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th><th>0.6</th><th>0.7</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><th>F<sub>I</sub><sup>A</sup></th><td>1.052</td><td>1.091</td><td>1.157</td><td>1.258</td><td>1.368</td><td>1.514</td><td>1.706</td><td>1.993</td><td>2.612</td></tr><tr><th>F<sub>I</sub><sup>B</sup></th><td>0.504</td><td>0.511</td><td>0.534</td><td>0.574</td><td>0.634</td><td>0.727</td><td>0.855</td><td>1.039</td><td>1.290</td></tr><tr><th>F<sub>II</sub><sup>B</sup></th><td>0.514</td><td>0.542</td><td>0.586</td><td>0.648</td><td>0.723</td><td>0.811</td><td>0.914</td><td>1.028</td><td>1.158</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>Figure 15. Ratios between the SIFs for a branched crack in a finite plate and infinite plate.</figcaption></figure>

<table><caption>Table IV. Normalized $K_1$ values for a cross in a finite plate for various ratio $a/w$.</caption><thead><tr><th>a/w</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th><th>0.6</th><th>0.7</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><th>X-FEM</th><td>0.8653</td><td>0.8844</td><td>0.9147</td><td>0.9572</td><td>1.0253</td><td>1.1348</td><td>1.3170</td><td>1.6388</td><td>2.4395</td></tr><tr><th>*</th><td>0.8641</td><td>0.8800</td><td>0.9092</td><td>0.9537</td><td>1.0223</td><td>1.1300</td><td>1.2866</td><td>1.4857</td><td>—</td></tr></tbody></table>

### 5.5. Two cracks emanating from a hole in a finite plate

Finally, we apply the X-FEM to two cracks emanating from a hole in a finite plate (Figure 18). To be able to compare our results with Newman [4], Bowie [27] and Woo *et al.* [28], the following dimensions for the half-height $H$ and the radius $r$ are taken: $H = 2w$, $r = 0.25w$. The normalized SIFs are defined as follows: $F_I = K_I/\sigma\sqrt{\pi a}$, $F_{II} = K_{II}/\sigma\sqrt{\pi a}$.

<footnote>*Reference solution by Cheung *et al.* [25].</footnote>
<figure><img src="image_17.png" /><figcaption>Figure 16. Cross in a finite square plate under bi-axial tension.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17. A star-shaped crack in a finite square under biaxial traction.</figcaption></figure>

<table><caption>Table V. Normalized $K_I$ and $K_{II}$ values for a star-shaped crack in a finite plate for various ratios $a/w$ and for two different meshes with about the same number of nodes: an unstructured mesh refined in the vicinity of the crack (mesh a) and a uniform mesh (mesh b).</caption><thead><tr><th rowspan="2">a/w</th><th colspan="3">F<sub>I</sub><sup>A</sup></th><th colspan="3">F<sub>I</sub><sup>B</sup></th><th colspan="3">F<sub>II</sub><sup>B</sup></th></tr><tr><th>Mesh a</th><th>Mesh b</th><th>*</th><th>Mesh a</th><th>Mesh b</th><th>*</th><th>Mesh a</th><th>Mesh b</th><th>*</th></tr></thead><tbody><tr><td>0.1</td><td>0.7511</td><td>†</td><td>0.7408</td><td>0.7690</td><td>†</td><td>0.7408</td><td>0.0001</td><td>†</td><td>0.0000</td></tr><tr><td>0.2</td><td>0.7670</td><td>0.7746</td><td>0.7570</td><td>0.7683</td><td>†</td><td>0.7578</td><td>0.0005</td><td>†</td><td>0.0004</td></tr><tr><td>0.3</td><td>0.7931</td><td>0.7942</td><td>0.7846</td><td>0.7983</td><td>0.7973</td><td>0.7884</td><td>0.0021</td><td>0.0021</td><td>0.0022</td></tr><tr><td>0.4</td><td>0.8287</td><td>0.8332</td><td>0.8255</td><td>0.8466</td><td>0.8466</td><td>0.8365</td><td>0.0080</td><td>0.0064</td><td>0.0070</td></tr><tr><td>0.5</td><td>0.8864</td><td>0.8928</td><td>0.8815</td><td>0.9255</td><td>0.9208</td><td>0.9087</td><td>0.0184</td><td>0.0168</td><td>0.0168</td></tr><tr><td>0.6</td><td>0.9673</td><td>0.9760</td><td>0.9758</td><td>1.0445</td><td>1.0401</td><td>1.0182</td><td>0.0364</td><td>0.0350</td><td>0.0338</td></tr><tr><td>0.7</td><td>1.0971</td><td>1.1120</td><td>1.1142</td><td>1.2367</td><td>1.2369</td><td>1.1936</td><td>0.0593</td><td>0.0614</td><td>0.0529</td></tr><tr><td>0.8</td><td>1.3423</td><td>1.3581</td><td>—</td><td>1.5624</td><td>1.5593</td><td>—</td><td>0.0864</td><td>0.0826</td><td>—</td></tr><tr><td>0.9</td><td>1.9037</td><td>1.9578</td><td>—</td><td>2.1927</td><td>2.1659</td><td>—</td><td>0.0868</td><td>0.0880</td><td>—</td></tr></tbody></table>

†The SIF was not computed since the J-integral domain contained several crack tips.

*Reference solution by Cheung *et al.* [25].

All computations are performed with a uniform 28 × 56 quadrilateral element mesh. As shown in Table VI, the X-FEM results are in good agreement with the reference solutions ($K_I$).

Table VII shows the computed SIFs $F_I$ and $F_{II}$ for several ratios of $a/w$ and different angles $\theta$. The plots in Figure 19 compare the X-FEM results with those obtained by Woo *et al.* [28]. We deduce from these results that the mesh is not sufficiently refined for ratios of $a/w$ less than 0.5.

## 6. CONCLUSIONS

An extended finite element method has been developed for modelling arbitrary geometries such as multiple branched cracks, voids and cracks emanating from holes without the need for the
<figure><img src="image_19.png" /><figcaption>Figure 18. Two cracks emanating from a hole in a finite plate.</figcaption></figure>

<table><caption>Table VI. Normalized $K_I$ values for different $a/w$.</caption><thead><tr><th>a/w</th><th>0.3</th><th>0.4</th><th>0.5</th><th>0.6</th><th>0.7</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><td>X-FEM</td><td>1.082</td><td>1.207</td><td>1.286</td><td>1.389</td><td>1.558</td><td>1.857</td><td>2.611</td></tr><tr><td>*</td><td>1.08</td><td>1.22</td><td>1.28</td><td>1.40</td><td>1.57</td><td>1.89</td><td>2.58</td></tr><tr><td>†</td><td>1.078</td><td>1.216</td><td>1.285</td><td>1.397</td><td>1.580</td><td>1.904</td><td>2.625</td></tr><tr><td>‡</td><td>1.078</td><td>1.216</td><td>1.285</td><td>1.396</td><td>1.574</td><td>1.892</td><td>2.498</td></tr></tbody></table>

<table><caption>Table VII. Normalized $K_I$ and $K_{II}$ values for two cracks emanating from a hole for various ratio $a/w$ obtained with the X-FEM.</caption><thead><tr><th>θ<br>a/w</th><th>0°<br>F<sub>I</sub></th><th>15°<br>F<sub>I</sub></th><th>30°<br>F<sub>I</sub></th><th>45°<br>F<sub>I</sub></th><th>60°<br>F<sub>I</sub></th></tr></thead><tbody><tr><td>0.4</td><td>1.207</td><td>1.120</td><td>0.864</td><td>0.517</td><td>0.155</td></tr><tr><td>0.5</td><td>1.286</td><td>1.178</td><td>0.943</td><td>0.586</td><td>0.212</td></tr><tr><td>0.6</td><td>1.389</td><td>1.288</td><td>1.031</td><td>0.654</td><td>0.260</td></tr><tr><td>0.7</td><td>1.558</td><td>1.459</td><td>1.145</td><td>0.726</td><td>0.301</td></tr><tr><td>0.8</td><td>1.857</td><td>1.702</td><td>1.297</td><td>0.812</td><td>0.348</td></tr><tr><td>0.9</td><td>2.611</td><td>2.178</td><td>1.531</td><td>0.916</td><td>0.388</td></tr></tbody></table>

<footnote>*Reference solution by Bowie [27].</footnote>

<footnote>†Reference solution by Newman [4].</footnote>

<footnote>‡Reference solution by Woo [28].</footnote>
<figure><img src="image_20.png" /><figcaption>Figure 19. Normalized $K_I$ (a) and $K_{II}$ (b) for two cracks emanating from a hole compared to the reference values $K_I^{ref}$ (a) and $K_{II}^{ref}$ obtained by Woo et al. [28].</figcaption></figure>

geometric entities to be meshed. The classical FEM approximation is enriched by the asymptotic near-tip field used by Belytschto and Black [11], the discontinuous function presented by Moës et al. [12] and a new discontinuous function to account for the multiple branched crack. The enrichment functions are introduced by employing the regular finite element shape functions as a partition of unity. The partition of unity has the desirable feature that the discrete equations retain the sparsity properties of the original mesh, albeit with an increase in the bandwidth. The voids are modelled by applying a quadrature scheme which partitions elements along the boundary of the voids and only consider the integrals within the body.

Unlike the boundary element methods, this method based on finite elements is readily generalized to problems involving non-linear constitutive laws.

The SIF-results using the X-FEM for various problems involving branched and intersecting cracks as well as cracks emanating from a hole were found to be in good agreement with reference solutions and analytical results.

In the extended finite element method, the crack is treated as a separate geometric entity and the only interaction with the mesh occurs through the selection of the enriched nodes and the numerical integration of the weak form. We have reported here quadrature schemes that use a subdivision of the elements cut by the cracks, with numerical integration performed separately on either side of the crack discontinuity.

The need of a variable number of degrees of freedom is a factor that contributes to greater computational cost of the new method. The relative ease of crack growth simulation without remeshing and the treatment of various flaws within a unified framework are, however, overriding advantages of the X-FEM over existing FE-based technology for modelling cracks.

## ACKNOWLEDGEMENTS

The support of the Office of Naval Research and Army Research Office, and the National Science Foundation through contract CMS-9732319 to Northwestern University, is gratefully acknowledged. The authors are grateful for the support provided by the DOE Computational Science Graduate Fellowship program, to John Dolbow.
REFERENCES

1. Barsoum R. Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements. *International Journal for Numerical Methods in Engineering* 1977; **11**:85-98.

2. Kwon YW, Akin JE. Development of a derivative singular element for application to crack propagation problems. *Computers and Structures* 1989; **31**(3):467-471.

3. Cruse T. *Boundary Element Analysis in Computational Fracture Mechanics*. Kluwer: Dordrecht, 1988.

4. Newman J. An improved method of collocation for the stress analysis of cracked plates with various shaped boundaries. *Technical Report TN D-6376, NASA, 1971.*

5. Nisitani H. Body force method for determination of the stress intensity factors. *Journal Aeronautical Society of India* (Special Issue on Fracture Mechanics) 1985; **37**:21-41.

6. Sneddon I. *Integral transform methods. In: Methods of Analysis and Solutions of Crack Problems*. Nordhoff International: Leyden, 1973.

7. Vitek V. Plane strain stress intensity factors for branched cracks. *International Journal of Fracture* 1977; **13**:481-501.

8. Chen Y, Hasebe N. New integration scheme for the branch crack problem. *Engineering Fracture Mechanics* 1995; **52**(5):791-801.

9. Obata M, Nemat-Nasser S, Goto Y. Branched cracks in anistropic elastic solids. *Journal of Applied Mechanics and Transactions ASME* 1989; **56**(4):858-864.

10. Fleming M, Chu YA, Moran B, Belytschko T. Enriched element-free Galerkin methods for crack tip fields. *International Journal for Numerical Methods in Engineering* 1997; **40**(8):1483-1504.

11. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. *International Journal for Numerical Methods in Engineering* 1999; **45**(5):601-620.

12. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. *International Journal for Numerical Methods in Engineering* 1999; **46**:131-150.

13. Dolbow J. An extended finite element method with discontinuous enrichment for applied mechanics. *Ph.D. Thesis*, Theoretical and applied mechanics, Northwestern University, 1999.

14. Melenk JM, Babuška I. The partition of unity finite element method: Basic theory and applications. *Computer Methods in Applied Mechanics and Engineering* 1996; **39**:289-314.

15. Strouboulis T, Copps K, Babuška I. The generalized finite element method: an example of its implementation and illustration of its performance. *International Journal for Numerical Methods in Engineering* 2000; **47**(8):1401-1417.

16. Strouboulis T, Babuška I, Copps K. The design and analysis of the generalized finite element method. *Computer Methods in Applied Mechanics and Engineering* 2000; **181**:43-71.

17. Duarte CA, Hamzeh ON, Liszka TJ, Tworzydło WW. The element partition method for the simulation of three-dimensional dynamic crack propagation. *Computer Methods in Applied Mechanics and Engineering* 2000, submitted for publication.

18. Oden J, Duarte C, Zienkiewicz O. New cloud-based hp finite element method. *Computer Methods in Applied Mechanics and Engineering* 1998; **153**(1):117-126.

19. Yau J, Wang S, Corten H. A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity. *Journal of Applied Mechanics* 1980; **47**:335-341.

20. Shih C, Asaro R. Elastic-plastic analysis of cracks on bimaterial interfaces: Part I – small scale yielding. *Journal of Applied Mechanics* 1988; **55**:299-316.

21. Kitagawa H, Yuuki R. Stress intensity factors for branched cracks in a two-dimensional stress state. *Transactions of the Japan Society of Mechanical Engineering* 1975; **41**:1641-1649.

22. Lo KK. Analysis of branched cracks. *Journal of Applied Mechanics* 1978; **45**:797-802.

23. Theocaris P, Tsamasphyros G, Theototokoglou E. An alternating coupling of finite elements and singular integral equations for the solution of branched cracks in finite sheets. *Engineering Fracture Mechanics* 1984; **20**:583-589.

24. Remacle J-F, Geuzaine C. *Gmsh finite element grid generator*, 1998, available at http://scorec.rpi.edu/~remacle/Gmsh_Eng.html.

25. Cheung Y, Wang Y, Woo C. A general method for multiple crack problems in a finite plate. *Computational Mechanics* 1992; **10**:335-343.

26. Ouchterlong F. Stress intensity factors for the expansion loaded star crack. *Engineering Fracture Mechanics* 1976; **8**:447-448.

27. Bowie O. Analysis of an infinite plate containing radial cracks originating at the boundary of an infinite circular hole. *Journal of Mathematical Physics* 1956; **35**:60-71.

28. Woo C, Wang Y, Cheung Y. The mixed mode problems for the cracks emanating from a circular hole in a finite plate. *Engineering Fracture Mechanics* 1989; **32**:279-288.