# 3D crack propagation with cohesive elements in the extended finite element method

Guilhem Ferté, Patrick Massin, Nicolas Moes

► To cite this version:

Guilhem Ferté, Patrick Massin, Nicolas Moes. 3D crack propagation with cohesive elements in the extended finite element method. Computer Methods in Applied Mechanics and Engineering, 2016, 300, pp.347-374. 10.1016/j.cma.2015.11.018. hal-01302550

HAL Id: hal-01302550

https://hal.science/hal-01302550

Submitted on 26 Jul 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
3D crack propagation with cohesive elements in the extended finite
element method

G. Ferté$^{a,*}$, P. Massin$^a$, N. Moës$^b$

$^a$IMSIA, UMR 9219 EDF-CNRS-CEA-ENSTA ParisTech, Université Paris Saclay, 828 Boulevard des Maréchaux, 91762 Palaiseau Cedex, France

$^b$GeM UMR Ecole Centrale de Nantes - Université de Nantes - CNRS, Ecole Centrale de Nantes, 1 rue de la Noë, BP 92101, 44321 Nantes, France

Abstract

A model is presented that accurately describes brittle failure in the presence of cohesive forces, with a particular focus on the prediction of non planar crack paths. In comparison with earlier literature, the originality of the procedure lies in the a posteriori computation of the crack advance from the equilibrium, instead of a most common determination beforehand from the stress state ahead of the front. To this aim, a robust way of introducing brittle non-smooth cohesive laws in the X-FEM is presented. Then the *a posteriori* update algorithm of the crack front is detailed. The crack deflection angle is computed from cohesive quantities exclusively, by introducing equivalent stress intensity factors. The procedure shows good accordance with experiments from the literature.

Keywords: X-FEM; Quasi-brittle fracture; Cohesive zone models; Size effect

**1. Introduction**

Four families of models are generally admitted to compute the brittle failure of components or structures : non local
approaches relying firstly on regularization techniques [1,2] with limits described in [3], first order [4–6] and higher
order [7,8] gradient approaches sharing similarities with damage phase-field models [9,10], cohesive zone models
(CZM) [11] and the linear elastic fracture mechanics (LEFM) [12,13]. While the first two types of models introduce
continuous deterioration of the mechanical properties in the bulk, strain localization is modelled by surfaces of
discontinuity in the others, interpreted as cracks propagating in the structure in the case of brittle failure. Recently, the
thick level-set model [14–16] brought together the two approaches into a *united formalism*: surfaces of discontinuity
automatically appear in some continuous damaged zones.

So that possibly non-planar surfaces of discontinuity may be inserted in a model regardless of the underlying mesh,
an extended finite element method (X-FEM) was developed by [17]. It accounts for discontinuities *within* elements,

<footnote>* Corresponding author.</footnote>

<footnote>*E-mail addresses:* guilhem.ferte@edf.fr (G. Ferté), patrick.massin@edf.fr (P. Massin), nicolas.moes@ec-nantes.fr (N. Moës).</footnote>
by means of a local enrichment with discontinuous functions, based on the partition of unity method [18]. Most of the
time, cracks are implicitly represented by so-called level-set functions (see [19–21]) and the X-FEM is coupled to a
level-set update algorithm and crack deflection criterion to model the propagation of a crack.

The level-set update consists of solving Hamilton–Jacobi equations, which can be done with implicit finite
difference schemes [19], explicit finite differences [22–24] – more efficiently used on an auxiliary grid [25,26] –
or Fast Marching Methods [27–29]. However, since the existing crack surface is preserved, the level-set update is
more readily done by simple geometrical equations in a fully explicit way, by taking advantage of the existence of a
moving crack front from which distances can be measured, as proposed by [30–32] with classical level-sets, by [33]
with vector level-sets or by [34] with an implicit–explicit crack description.

Though the LEFM is the most mature engineering tool to study crack propagation, it sometimes shows limitations
for some advanced applications:

* If the size of the process zone approaches any relevant lengthscale of the model (structure typical size, initial crack length, distance between crack tips), the response of the structure is subjected to a *size effect*, that LEFM fails to predict.

* In 3D, it is very difficult to track the new crack front as an iso-G curve, as Griffith theory would suggest. Although some recent work based on the use of configurational mechanics tackle this issue in an implicit way [35], a rule of thumb is generally used to compute an explicit crack advance instead: for instance, the crack front is moved forward proportionally to G (see [34] for discussion).

* The crack may not initiate from a sound structure: a pre-crack has to be put in. Thus, the sensitivity of the response to the shape and orientation of the pre-crack has to be studied, and the failure load cannot always be accurately reproduced.

* The calculation of stress intensity factors (SIF) requires the construction of adequate contours surrounding the front [36–38], or specific elements to achieve optimal quality of the results [39,40], both techniques being generally coupled, which requires additional monitoring.

* The notion of stress intensity factor is originally intended to linear isotropic homogeneous elasticity. Extensions to anisotropic materials [41–43] and large deformations [44] do exist under some assumptions. Elastic–plastic materials may also be treated through asymptotic expansions near the crack tip for power-hardening laws by Hutchinson, Rice and Rosengren [45–47], but related solutions are not valid for (even local) unloading.

In order to address some of these issues, cohesive zone models were originally proposed by Barenblatt [11]. They consist of restraining forces which vanish as the crack opens, thus obeying a softening traction–displacement relation. This behaviour was introduced in the X-FEM and coupled with crack propagation algorithms as illustrated by the literature cited below.

In a first family of approaches (Fig. 1(a)), the crack path is fully extended before computing the equilibrium. While conserving the crack surface, the load increment is applied. The extension of the crack path is then determined from the stress state ahead of the front, where elements are going to exceed a critical stress. The cohesive calculation is performed with this new path, and the procedure is repeated iteratively (Fig. 1(a)). The crack advance increment is often monitored, as in [48,49], so as to prevent spurious oscillations or too rough results.

The crack path may be extended by segments through the determination of a deflection angle and crack advance, as in [48–51]. In [50,51] for instance, a bidimensional crack is propagated perpendicular to the direction of an averaged maximal principal stress in all cut elements, as long as a critical value is exceeded in the element. Alternatively, this crack advance may be determined by setting the overall energy variation (including cohesive energy) to zero [49].

However, the crack path is often determined as a whole, by so-called crack tracking algorithms: Jäger [52] offers a nice review of these methods. Areias and Belytschko [53] first proposed a local version: for each element adjacent to the current crack front, the crack surface is iteratively extended based on the adjacent preexisting crack front. Gasser and Holzapfel [54,55] proposed a non local version where the crack is iteratively extended based on principal directions and the pre-existing surface in a neighbourhood of the element to be cut. As this smoothes the crack front, it gives better results for sharply curved cracks. In a similar manner, de Borst et al. [56,57] proposed a cohesive segment approach that would initiate and propagate cracks in heterogeneous materials based on non-local stresses. Finally, Oliver et al. [58,59] proposed a global crack tracking algorithm where the crack path is determined from the field of principal directions and previous crack front by a global resolution of an anisotropic heat equation.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Existing approaches in the literature to handle cohesive propagation.</figcaption></figure>

In a second family of approaches (Fig. 1(b)), the crack advance is controlled at each step: from a given crack path, a deflection angle is determined. The crack is extended by a controlled length in this direction and the load is adapted so that the cohesive zone be coincident with the new crack tip (see [60–62]), this coincidence being ensured by setting the overall stress intensity factor (SIF) to zero.

As for the crack deflection angle, criteria rely either on stress intensity factors (SIF), local stress, or variational formulations.

Based on the SIF are the maximal hoop stress criterion of Erdogan and Sih [63], the minimum strain energy density criterion [64], the maximum energy release rate criterion by [65,66], the maximum KI or minimum KII criterion (see [67]). Now, in the LEFM framework, these SIF may be computed by integral methods [68,48] or related to the stress around the tip [34,69]. Equivalent stress intensity factors may also be computed straight from a cohesive result with integrals surrounding the process zone [60].

As for criteria based on local stresses, for brittle failure the direction of the maximal principal stress is most commonly considered as being normal to the propagation direction. This was converted into an angle expression in 2D and 3D by [69,70]. Some papers [62] recommend not to take axial stress into account when computing the principal directions. Finally, in [71] the crack angle and directions are included as new unknowns in the energy of the 2D structure to be minimized.

In all these techniques, the surface of discontinuity has to stop within the structure. The X-FEM enrichment strategies either have it stopped at the edges (2D) or faces of the mesh, as in [50,51] or as in the PU-FEM methods [55–57]. Alternatively, a crack tip may be located within an element, which would be enriched with regular functions as in [60,48,49,72]. An alternative is proposed by [61] to locate a crack tip within an element with no other enrichment than the Heaviside.

In this paper a new procedure is proposed where the crack advance is implicit: it is naturally given as a side result from cohesive fields at equilibrium. In this way, the crack advance is kind of embedded in the unknowns of the problems rather than postulated in advance with some criterion that would have to be discussed. It would be even better to also determine the crack deflection angle in an analogous way, but this is complex in 3D and will not be addressed in this paper. First, in Section 2, a robust implementation of non-smooth cohesive laws in the X-FEM is presented, that is capable of handling initial perfect adherence, even if it occurs over broad areas. This turns out to be necessary when the next location of the crack front is not assumed to be known beforehand. Then, the crack path prediction procedure is given in Section 3. Finally, numerical results are presented in Section 4.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Cohesive zone model.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Notations of the problem.</figcaption></figure>

## 2. XFEM cohesive zone models with large adherent zone

The cohesive law is thus inserted on a broad potential crack surface (see Fig. 3) and the equilibrium will naturally separate the cracked surface – which comprises a traction-free zone and a cohesive zone – from the adherent zone (see Fig. 2). A robust implementation of non-smooth cohesive laws in the X-FEM is thus presented here, that is capable of describing broad adherent zones.

In quasi-brittle materials such as concrete, the cohesive zone models the area where micro-cracking has occurred, but micro-cracks have not yet coalesced into a single macro-crack. The boundary between traction-free and cohesive zone thus locates the physical crack front of the observable macro-crack. On the contrary, the boundary between cohesive and adherent zone only locates the onset of micro-cracking: it is the front of a fictitious crack representing the zone of incipient material debonding ahead of the physical crack. In a purposeful misuse of language, the “fictitious crack front” is simply referred to as “crack front” from here onward.

### 2.1. Equilibrium of cracked bodies with cohesive forces

Let $\Omega \subset \mathbb{R}^d$ be a body that is cut by a potential crack surface $\Gamma$. Its external boundary $\partial\Omega$ is decomposed into non-overlapping parts $\Gamma_u$ and $\Gamma_g$, with prescribed displacements and forces $g$, respectively (see Fig. 3). If $u^i$ are the displacement fields on each of the crack lip $\Gamma^i$ ($\Gamma^+$ is the upper crack lip, and $\Gamma^-$ the lower crack lip), the displacement jump is $[u](x) = u^+(x) - u^-(x)$. Let $n$ be the outward normal vector to $\Gamma^-$ and $t_c$ be the restraining (cohesive) force that $\Gamma^+$ applies over $\Gamma^-$ (Fig. 3). Therefore, the stress tensor $\sigma$ satisfies $\sigma^+ \cdot n = \sigma^- \cdot n = t_c$ on $\Gamma$.
<figure><img src="image_5.png" /><figcaption>Fig. 4. 3D description of the crack and the covariant basis attached to the crack front.</figcaption></figure>

## 2.2. Representation of the crack by means of level-set functions

The crack has an implicit representation. A normal level-set function $\phi_n$ defines the position of the cracked surface and a $\langle \text{potential} \rangle$ extension of it, through the isosurface $\phi_n = 0$, which we call potential crack surface in this paper. A tangential level-set $\phi_t$ allows to discriminate between the cracked surface and the adherent zone (see Fig. 2), so that:

$$ \Gamma = \{ \mathbf{x}, \phi_n(\mathbf{x}) = 0 \text{ and } \phi_t(\mathbf{x}) \le 0 \}. \qquad (1) $$

The crack front is then represented by $T: (\phi_n = 0) \cap (\phi_t = 0)$ (see Fig. 4). This allows to define the covariant basis along the front (see Fig. 4) as $\mathbf{n} := \nabla\phi_n, \mathbf{t} := \nabla\phi_t$ and $\mathbf{b} := \nabla\phi_n \times \nabla\phi_t$.

## 2.3. Discrete displacement

The elements crossed by the potential crack are all enriched with a Heaviside-like function $H$, so that the discrete displacement space is:

$$ V_h := \left\{ \sum_{i \in \mathcal{N}} \mathbf{a}_i N_i(\mathbf{x}) + \sum_{i \in \mathcal{K}} \mathbf{b}_i N_i(\mathbf{x}) H(\mathbf{x}), \mathbf{a}_i \in \mathbb{R}^d, \mathbf{b}_i \in \mathbb{R}^d \right\} \qquad (2) $$

where $\mathcal{N}$ represents the set of nodes of the mesh, $\mathcal{K}$ is the set of enriched nodes (see Fig. 5), $N_i$ are the linear shape functions (linear elements are considered throughout this paper) and:

$$ H(\mathbf{x}) := \begin{cases} +1 & \text{if } \phi_n(\mathbf{x}) \ge 0 \\ -1 & \text{if } \phi_n(\mathbf{x}) < 0. \end{cases} \qquad (3) $$

## 2.4. Why a penalized law is not efficient to describe a large adherent zone

As a trial run, the cohesive traction $t_c$ is written as an explicit function of $[u]$ (see Fig. 6). Hence, an equivalent jump is defined as $[u]_{eq} = \sqrt{\langle [u]_n \rangle_+^2 + \beta^2 [u]_s^2}$, where:

* $[u]_n = [u] \cdot \mathbf{n}$ is the normal jump;

* $\langle [u]_n \rangle_+ = \max \{[u]_n, 0\}$ and $\langle [u]_n \rangle_- = \min \{[u]_n, 0\}$ are its positive and negative parts;

* $[u]_s = (\mathbf{1} - \mathbf{n} \otimes \mathbf{n}) \cdot [u]$ is the tangential jump;

* $\beta$ is a material parameter that quantifies the ratio between tensile and shear strengths.

A dimensionless internal variable $\alpha$ may then be defined at time $t$ as:

$$ \alpha(t) = \max_{t' \in [0, t]} \frac{[u]_{eq}(t')}{w_c} \qquad (4) $$
<figure><img src="image_6.png" /><figcaption>Fig. 5. Potential crack surface not matching the mesh edges, and subsequent enriched nodes.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Penalized linear-softening mode-coupling cohesive law.</figcaption></figure>

where $w_c$ corresponds to the first opening for which cohesive forces are reduced to zero. The internal variable has non-zero initial value $\alpha_0$: $\alpha = \alpha_0$ indicates adherence with the initial stiffness, $\alpha \in ]\alpha_0, 1[$ describes a damaged material, and $\alpha = 1$ characterizes full debonding. For loading conditions (that is, if $\alpha = [u]_{eq}/w_c$), the surface energy reads:

$$ \Pi([u]) = 2G_c\alpha \left(1 - \frac{\alpha}{2}\right) + \frac{\sigma_c}{2\alpha_0 w_c} \langle [u]_n \rangle^2. \quad (5) $$

The cohesive traction is deduced from (5) by $t_c([u]) = \frac{\partial\Pi}{\partial[u]}$. Its expression may be summarized as follows: defining an equivalent traction as $t_{c,eq} = \sqrt{\langle t_{c,n} \rangle_+^2 + \beta^{-2}t_{c,s}^2}$, it holds:

$$ t_{c,eq} = \sigma_c (1 - \alpha) \frac{[u]_{eq}}{\alpha w_c}. \quad (6) $$

The normal traction is then defined by $t_{c,n} = \frac{t_{c,eq}}{[u]_{eq}} \langle [u]_n \rangle_+ + \sigma_c \frac{\langle [u]_n \rangle_-}{\alpha_0 w_c}$ and the shear component is $t_{c,s} = \beta^2 \frac{t_{c,eq}}{[u]_{eq}} [u]_s$.

To see how such a formulation performs with large adherent zones, an inclusion debonding test is carried out, whose geometry and loading conditions are represented in Fig. 7 (all dimensions are millimetres). It consists of a plate in plane strain under tension, for which linear isotropic elasticity is assumed with coefficients $E = 36.56$ GPa and $\nu = 0.2$. A circular inclusion is inserted into the plate, which is prone to debonding and therefore subjected to the above cohesive law, with tensile strength $\sigma_c = 2.7$ MPa, fracture energy $G_c = 0.095$ N mm$^{-1}$ and mixed-mode parameter $\beta = 1$. With strong mode-coupling, a curved crack surface and some contact at the upper and lower extremities of the inclusion, this problem is a good challenger for the accuracy and robustness of a formulation.

Normal opening and tractions – these being computed straight from the displacement by (6) – have been plotted on Figs. 8–9 at load $u_0 = 0.04$ mm, for crack-surface-matching and non matching meshes, and with various penalization parameters $\alpha_0^{-1}$. It can be seen on Figs. 8–9 that such classical penalized cohesive laws $t_c([u])$ raise three issues:
<figure><img src="image_8.png" /><figcaption>Fig. 7. Inclusion debonding test.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Normal opening with the inclusion debonding test.</figcaption></figure>

1. *non physical opening* or interpenetration is observed, impairing the accuracy of the results, when penalization parameter $\alpha_0^{-1}$ is too low (tag 1 on Fig. 8): $\alpha_0^{-1} = 10^4$ is the minimal value which accurately resolves the adherent phase;

2. *spurious oscillations* of the traction are observed in the (almost) adherent zone, when $\alpha_0^{-1}$ is on the contrary too high (tag 2 on Fig. 9). This numerical issue has been widely reported and studied by X-FEM literature (see [73–75]). In one word, it comes from the fact that the discretization space for traction gets too rich in comparison with that of the displacement, when the stiffness of the interfacial law gets so high that it is akin to enforcing a Dirichlet condition;

3. *misevaluation of the regime* (adherence or debonding) may occur as a collateral damage of these oscillations, with some Gauss points in the adherent zone being mistakenly considered as dissipative (tag 3 on Fig. 9).

To conclude, classical penalized laws fail to achieve good resolution of large adherent zones, since a too low $\alpha_0^{-1}$ implies large opening in the adherent zone leading to a wrong solution in terms of accuracy, while a higher value of $\alpha_0^{-1}$ arises spurious oscillations due to stability issues, as soon as the penalization is about to get sufficiently stiff to correctly describe adherence.

Changes to the formulation are proposed as remedial actions to each of the aforementioned drawbacks:

1. a *mixed law with initial perfect adherence* is proposed to avoid non physical opening;

2. a *stable (mortar) formulation* is introduced to handle perfect adherence and the sudden switch to a Neumann-like condition when debonding starts;

3. a discretization with *blockwise diagonal operators* is suggested to avoid regime misevaluations.

## 2.5. The use of Lagrange multipliers in the X-FEM

The ability to prescribe Dirichlet conditions in the X-FEM is a prerequisite to introduce mixed laws. Literature on the topic can be divided into *stable* methods, which works out discrete spaces by enriching the displacement
<figure><img src="image_10.png" /><figcaption>Fig. 9. Normal tractions with the inclusion debonding test.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Mesh not matching the crack surface and reduced multiplier space.</figcaption></figure>

space [76] or reducing the 〈 traction 〉 multiplier vector space [74,75,77], and stabilized methods, which consists of adding stabilizing terms to the formulation, either by having the stress trace on the interface play the role of Lagrange multiplier, as in Nitsche's methods [78], or with stabilization term on the multiplier/stress discrepancy [79]. In this paper, the multiplier reduced vector space by Géniaut [77] is used.

In this space, multipliers components are supported by the nodes in $\mathcal{K}$ (see Fig. 5). This initial set is reduced into a fewer number $N_\lambda$ of degrees of freedom by prescribing equality relations between the components, supported by some intersected edges $\mathcal{V}$: a truly independent Lagrange degree of freedom $I$ is shared by a group of nodes $i$ in $\mathcal{K}$ (see Fig. 10), hence making a non-local shape function $\psi_I := \sum_{i \in I} N_i$ (see Fig. 10)—it is called non-local because its support is made up of non-adjacent elements (see Fig. 10). The multiplier field is then interpolated over cut elements and the discrete multiplier is finally the trace of this field on the crack surface, so that the multiplier reduced vector space is:

$$ M_h := \left\{ \sum_I \mu_I \psi_I |_{\Gamma_I}, \mu_I \in \mathbb{R}^d \right\}. $$

## 2.6. A mixed law, in the augmented lagrangian formalism

For the sake of conciseness, let us call $w$ the displacement jump. In the augmented Lagrangian formalism, a general expression of the surface density of energy for orthotropic potential laws is:

$$ \Pi(w, \lambda) = \phi(\lambda_n + r w_n, \lambda_s + r_s w_s) - \frac{\lambda_n^2}{2r} - \frac{\lambda_s \cdot \lambda_s}{2r_s} \quad (7) $$

where $\phi$ is a derivable function, $r$ and $r_s$ are the normal and tangential augmentation parameters. Here, tangential augmentation is taken as $r_s = \beta^2 r$.
<figure><img src="image_12.png" /><figcaption>Fig. 11. Mode-coupling mixed cohesive law.</figcaption></figure>

While the cohesive traction is still given by $t_c = \frac{\partial \Pi}{\partial w}$, dual equations are needed to know where Dirichlet conditions take place, and enforce them there. This additional interfacial law is determined by stating that $\frac{\partial \Pi}{\partial \lambda} = 0$ in the augmented Lagrangian formalism, which ensures that there is no need to verify any further inequality constraint throughout computation: all constraints are embedded in that (equality) interfacial law, which makes it eligible for resolution with a mere Newton-Raphson procedure.

The cohesive traction thus reads:

$$t_{c,n}(\lambda_n + r w_n, \lambda_s + r_s w_s) = \frac{\partial \Pi}{\partial w_n} = r \frac{\partial \phi}{\partial(\lambda_n + r w_n)} \quad (8)$$

$$t_{c,s}(\lambda_n + r w_n, \lambda_s + r \beta^2 w_s) = \frac{\partial \Pi}{\partial w_s} = r \beta^2 \frac{\partial \phi}{\partial(\lambda_s + r \beta^2 w_s)}. \quad (9)$$

The dependence of $t_c$ on $\lambda_n + r w_n$ and $\lambda_s + r_s w_s$ will be omitted from here onward to alleviate notations. The interfacial law is simply $\lambda = t_c$ (see Fig. 11).

Let $(\lambda + rw)_{eq} := \sqrt{(\lambda_n + r w_n)_+^2 + \beta^{-2} (\lambda_s + \beta^2 r w_s)^2}$ be an equivalent augmented traction. A threshold function $\varphi$ is then introduced as $\varphi((\lambda + rw)_{eq}) := \frac{(\lambda + rw)_{eq} - \sigma_c}{r w_c - \sigma_c}$, so that a scalar dimensionless internal variable $\alpha$ is defined as verifying:

$$\varphi((\lambda + rw)_{eq}) - \alpha \le 0 \quad (10)$$

$$\dot{\alpha} \ge 0 \quad (11)$$

$$\dot{\alpha} [\varphi ((\lambda + rw)_{eq}) - \alpha] = 0. \quad (12)$$

It holds $\alpha \le 0$ for an uncracked material, and $\alpha \ge 1$ for a fully cracked material. For loading conditions, that is to say if $\alpha = \varphi((\lambda + rw)_{eq})$, function $\varphi$ is defined by:

$$\phi(\lambda_n + r w_n, \lambda_s + r_s w_s) = 2G_c \left(1 - \frac{\sigma_c}{r w_c}\right) \alpha \left(1 - \frac{\alpha}{2}\right) + \frac{1}{2r} \langle \lambda_n + r w_n \rangle_c^2. \quad (13)$$

For contact-free situations, and not considering the related term in (13), the surface energy (7) depends upon $\alpha$ and $\lambda_{eq}$ only, through $\Pi(\alpha, \lambda_{eq}) = \phi(\alpha) - \frac{\lambda_{eq}^2}{2r}$. When dissipation starts, it holds $\alpha = 0$ and $\lambda_{eq} = \sigma_c$. When it ends, it holds $\alpha = 1$ and $\lambda_{eq} = 0$. Then, we have $\pi(\alpha = 1, \lambda_{eq} = 0) - \pi(\alpha = 0, \lambda_{eq} = \sigma_c) = G_c$ since $\phi(\alpha = 1) - \phi(\alpha = 0) = G_c - \frac{\sigma_c^2}{2r}$, which ensures that an energy $G_c$ shall be provided to fully open a unit surface of crack.

The resulting traction deduced from (13) by (8)-(9) may be synthesized as follows: still defining equivalent tractions as $t_{c,eq} = \sqrt{\langle t_{c,n} \rangle_+^2 + \beta^{-2} t_{c,s}^2}$, the traction-separation law is expressed in terms of equivalent quantities
as $t_{c,eq} = (1 - T_d) (\lambda + rw)_{eq}$, where $T_d$ is the damage tensor. Its expression for linear softening is:

$$T_d = \frac{\alpha}{\left(1 - \frac{\sigma_c}{r w_c}\right) \alpha + \frac{\sigma_c}{r w_c}}. \qquad (14)$$

We may readily check that $T_d = 0$ if $\alpha = 0$ (perfect adherence), and $T_d = 1$ if $\alpha = 1$ (full debonding). The vector traction-separation law is then defined in terms of normal component $t_{c,n} = (1 - T_d) (\lambda_n + rw_n)_+ + (\lambda_n + rw_n)_-$ (the latter term being added to account for unilateral contact) and shear component $t_{c,s} = (1 - T_d) (\lambda_s + r\beta^2 w_s)$.

## 2.7. Stable ≪ mortar ≫ formulation for inserting an interfacial law in XFEM

As mentioned before, interface quantities (tractions $t_c$ and energy density $\Pi$) should be defined over the reduced space $M_h$ to avoid spurious oscillations during adherence phases. This still applies to the displacement jump $w$, because it actually comes into play during adherence phases as well (through the augmentation). Hence, $w$ is introduced as a new unknown of the problem, to be discretized over a different space from that of $[u]: M_h$. The total energy of the problem reads:

$$E(u, \lambda, w) = \frac{1}{2} \int_{\Omega} \epsilon(u) : C : \epsilon(u) d\Omega - \int_{\Gamma_g} g \cdot u d\Gamma_g + \int_{\Gamma} \Pi(w, \lambda) d\Gamma. \qquad (15)$$

Now, finding the solution of the continuous problem would imply to find $(u, w, \lambda) = \text{argmin}_{w^*=[u^*]} E(u^*, \lambda^*, w^*)$, hence the Lagrangian of the problem:

$$\mathcal{L}(u, w, \lambda, \mu) = \frac{1}{2} \int_{\Omega} \epsilon(u) : C : \epsilon(u) d\Omega - \int_{\Gamma_g} g \cdot u d\Gamma_g + \int_{\Gamma} \Pi(w, \lambda) d\Gamma + \int_{\Gamma} \mu \cdot ([u] - w) d\Gamma. \qquad (16)$$

The optimality conditions of Lagrangian (16) give the following discrete weak form:

$$\forall u^* \in V_h, \quad \int_{\Omega} \sigma(u) : \epsilon(u^*) d\Omega - \int_{\Gamma_g} g \cdot u^* d\Gamma_g + \int_{\Gamma} \mu \cdot [u^*] d\Gamma = 0 \qquad (17)$$

$$\forall \mu^* \in M_h, \quad \int_{\Gamma} ([u] - w) \cdot \mu^* d\Gamma = 0 \qquad (18)$$

$$\forall w^* \in M_h, \quad -\int_{\Gamma} [\mu - t_c] \cdot w^* d\Gamma = 0 \qquad (19)$$

$$\forall \lambda_n^* \in M_h, \quad -\int_{\Gamma} \frac{[\lambda_n - t_{c,n}]}{r} \cdot \lambda_n^* d\Gamma = 0 \qquad (20)$$

$$\forall \lambda_s^* \in M_h, \quad -\int_{\Gamma} \frac{[\lambda_s - t_{c,s}]}{r\beta^2} \cdot \lambda_s^* d\Gamma = 0. \qquad (21)$$

## 2.8. Blockwise diagonal discrete operators at the interface

Brittle cohesive laws exhibit sudden changes of tangent behaviour matrix, and as such belong to the family of non-smooth interface laws (as unilateral contact for instance). It is thus better, to prevent cyclic ≪ flip-flop ≫ behaviour of a Newton-Raphson algorithm, to limit the number of points where these changes may occur to the number $N_\lambda$ of Lagrange degrees of freedom, which leads to the discretization adopted below.

The components of unknown vectors $u$ and $\mu$ are defined in a fixed basis $(e_X, e_Y, e_Z)$, while the components of unknown vectors $w$ and $\lambda$ are defined in the covariant basis $(n, t, b)$ to the crack surface $\Gamma$ at each point $x \in \Gamma$ (see Fig. 4), so that:

$$w(x) = \sum_{I=1}^{N_{\lambda}} \psi_I(x) (w_{I,n}n(x) + w_{I,t}t(x) + w_{I,b}b(x)). \qquad (22)$$
A similar equation to (22) holds for $\lambda$. At a Lagrange degree of freedom $I \in \{1..N_\lambda\}$, it is possible to compute cohesive components $t_{c,n}^I, t_{c,t}^I, t_{c,b}^I$ from components $(w_{I,n}, w_{I,t}, w_{I,b})$ and $(\lambda_{I,n}, \lambda_{I,t}, \lambda_{I,b})$ with the aforementioned cohesive law. These cohesive components are not meant to be associated with particular directions around degree of freedom $I$, but they are meant to be used in a weak sense. The discretization procedure is as follows:

* for (17), classical Gaussian integration is applied;

* for (18), classical Gaussian integration is performed in the fixed basis $\forall I \in \{1..N_\lambda\}, \int_\Gamma ([u_\chi](x) - w(x) \cdot e_\chi) \psi_I(x)d\Gamma = 0$, and similarly for directions $e_Y, e_Z$;

* for (19), a lump is applied to the cohesive traction only, so that, (19) becomes in the covariant basis $\forall I \in \{1..N_\lambda\}, t_{c,n}^I \int_\Gamma \psi_I(x)d\Gamma - \int_\Gamma \mu(x) \cdot n(x) \psi_I(x)d\Gamma = 0$ and similarly for directions $t, b$;

* as for (20)–(21), which represents the interfacial law, it is lumped everywhere, so that $\forall I \in \{1..N_\lambda\}, \frac{\int_\Gamma \psi_I(x)d\Gamma}{r}(t_{c,n}^I - \lambda_{I,n}) = 0$ and similarly for directions $t, b$. In this expression, the multiplying factor has been kept to have a symmetric tangent stiffness matrix.

This highly non-linear set of equations is solved with a Newton-Raphson algorithm. The tangent stiffness matrix of the problem reads:

$$ K := \begin{bmatrix} k^{uu} & (k^{\mu u})^T & 0 & 0 \\ k^{\mu u} & 0 & (-k^{w \mu})^T & 0 \\ 0 & -k^{w \mu} & d^{ww} & (d^{\lambda w})^T \\ 0 & 0 & d^{\lambda w} & d^{\lambda \lambda} \end{bmatrix} \qquad (23) $$

where:

* $k^{uu}$ is the bulk stiffness matrix;

* $k^{\mu u}$ and $k^{w \mu}$ are classical $\langle\langle$ mass $\rangle\rangle$ matrices stemming from the discretization of $\langle\langle$ mortar $\rangle\rangle$ operators, the latter also handling the change from fixed to covariant basis;

* matrices $d$ are all blockwise diagonal: for $I$ and $J$ distinct Lagrange degrees of freedom, it holds $d_{IJ} = 0$.

Owing to (8)–(9) we have $\frac{\partial t_{c,s}}{\partial(\lambda_n + r w_n)} = \beta^2 \left(\frac{\partial t_{c,n}}{\partial(\lambda_s + r \beta^2 w_s)}\right)^T$ so that a symmetric matrix derived from the cohesive law is introduced as:

$$ k^c := \begin{bmatrix} \frac{\partial t_{c,n}}{\partial (\lambda_n + r w_n)} & \beta^2 \frac{\partial t_{c,n}}{\partial (\lambda_s + r \beta^2 w_s)} \\ \frac{\partial t_{c,s}}{\partial (\lambda_n + r w_n)} & \beta^2 \frac{\partial t_{c,s}}{\partial (\lambda_s + r \beta^2 w_s)} \end{bmatrix} = \begin{bmatrix} k_{nn}^c & k_{ns}^c \\ k_{sn}^c & k_{ss}^c \end{bmatrix}. \qquad (24) $$

Diagonal blocks are then determined by:

$$ d_{II}^{ww} = rk^c \qquad (25) $$

$$ d_{II}^{\lambda w} = \begin{bmatrix} k_{nn}^c & k_{ns}^c \\ \beta^{-2}k_{ns}^c & \beta^{-2}k_{ss}^c \end{bmatrix} \qquad (26) $$

$$ d_{II}^{\lambda \lambda} = \begin{bmatrix} \frac{k_{nn}^c - 1}{r} & \text{sym.} \\ \frac{k_{ns}^c}{r\beta^2} & \frac{\beta^{-2}k_{ss}^c - 1}{r\beta^2} \end{bmatrix}. \qquad (27) $$

## 2.9. Numerical validation with the inclusion debonding test

The above three changes (mixed law, stable $\langle\langle$ mortar $\rangle\rangle$ formulation, blockwise diagonal operators) may be applied independently from one another. To assess their individual numerical effect, they were tested incrementally in intermediate formulations, which are all summarized on Table 1.
<table><caption>Table 1
Tested formulations and total number of required Newton iterations. The load is initially prescribed in three time steps, which may be further subdivided whenever needed to converge.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Classical formulation</th><th colspan="2">Stable «mortar» formulation</th></tr><tr><th>Consistent operators</th><th>Blockwise diagonal operators</th></tr></thead><tbody><tr><td>Penalized law</td><td>t<sub>c</sub>([u])</td><td>t<sub>c</sub>(w)</td><td>t<sub>c</sub>(w)</td></tr><tr><td>Newton iterations</td><td>729</td><td>311</td><td>147</td></tr><tr><td>Mixed law</td><td>t<sub>c</sub>([u], λ)</td><td>t<sub>c</sub>(w, λ)</td><td>t<sub>c</sub>(w, λ)</td></tr><tr><td>Newton iterations</td><td>63</td><td>43</td><td>15</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Total number of Newton iterations to solve the problem in 3 time steps.</figcaption></figure>

The inclusion debonding is considered (Fig. 7), with a load $u_0 = 0.2$ mm being originally applied in 3 time steps. At a given step, if the Newton-Raphson algorithm fails to converge, computation restarts with a load increment that is half its previous value. As an indicator for robustness, we take the total number of iterations needed to apply the full load: it is obtained by summing up the Newton iterations of all converged load increments (there can be much more than three, if several subdivisions had to be applied to reach convergence), and reported in Table 1 and Fig. 12. For penalized laws, the smaller penalization parameter which still ensures adequate adherence (non physical opening was not observed) is $\alpha_0^{-1} = 10^4$. This value was chosen in Table 1, and only higher values are tested on Fig. 12. The total number of Newton iterations is then observed to decrease significantly for each of the changes, illustrating the ability of each of them to bring additional robustness (see Table 1 and Fig. 12).

Moreover, normal cohesive tractions have been plotted along the inclusion in Fig. 13, for $u_0 = 0.04$ mm. It can be deduced from this plot that:

* the stable multiplier space, used along with a mixed law allows to remedy the aforementioned spurious oscillations issue;

* the implementation of mixed cohesive law is checked, since its results coincide with those of the linear-softening penalized law;

* neither the stable «mortar» formulation, nor blockwise diagonal operators alter the results, since they coincide with other laws.

### 3. 3D cohesive crack propagation with implicit crack advance

The procedure that we propose to study crack propagation is as follows (Fig. 14):

* given a potential crack surface, the equilibrium state is computed (Fig. 14(b)),

* update of the crack front based on the computed cohesive state (Fig. 14(c)),
<figure><img src="image_14.png" /><figcaption>Fig. 13. Cohesive stress along the interface for intermediate formulations.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Overview of the procedure.</figcaption></figure>

* determination of deflection angles along the front (Fig. 14(d)),

* update of the potential crack surface accordingly (Fig. 14(e)),

* correction of the cohesive internal variables (Fig. 14(a)).

### 3.1. Update of the crack front

This is the main novelty of our propagation procedure. When used on indefinite paths, the cohesive crack advance is generally determined beforehand. On the contrary, when used on predefinite paths, cohesive laws do not require a priori knowledge of the crack front: that information is naturally embedded in the model, which is an advantage of cohesive laws in comparison with LEFM. Hence, we would like to benefit from this feature for an indefinite path as well.

We shall then set up a detection phase, in which the tangential level-set – and consequently the crack front – should be updated a posteriori, based on the new computed cohesive state and on the crack front from the previous load
<figure><img src="image_16.png" /></figure>

Upcoming views :
representation of the crack
surface intersecting the mesh

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

(b) Computation of intersected
elements

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

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

(a) Nodal values of the internal
variables, negative values for
undamaged material

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

(c) Intersection points with the
mesh faces

<caption>Fig. 15. Computation of a rough crack front.</caption>

increment. A requirement for that detection is that it should produce a smooth $\phi_t$: it has to be suitable for further use by level-set update algorithms during the next propagation steps. To do so, a first rough front detection is performed using the cohesive internal variable state only (Fig. 15): that field is similar to a tangential level-set, but it is too rough for a straight use in the level-set update procedure. A smooth crack advance field will therefore be recreated a posteriori, from the former crack front (Fig. 16(a)–(b)), and put into a level-set update algorithm which will produce a smoother tangential level-set (Fig. 16(c)–(d)).

The detection process is decomposed as follows:

1. Rough crack front detection (Fig. 15)

(a) The nodal field of internal variables $\alpha$ is computed, which is allowed to be negative for the adherent part (Fig. 15(a)).

(b) Elements which are intersected by the isozero of that field are computed (Fig. 15(b)).

(c) Intersection points of the rough crack front with the faces of the mesh are computed, as in [25], creating a cloud of points (Fig. 15(c)).

2. Reconstruction of a smooth tangential level-set (Fig. 16)

(a) Construction of a crack advance field at each point *P* of the previous front, as the distance from *P* to the projection of the cloud of points in the propagation plane (*P*, *n*<sub>*P*</sub>, *t*<sub>*P*</sub>) (Fig. 16(a)).

(b) Smoothing of that advance field as a scalar function (Fig. 16(b)).

(c) Application of a level-set update algorithm (to the tangential level-set only) with that advance field and the deflection angle computed at earlier steps (Fig. 16(c)–(d)).
<figure><img src="image_21.png" /><figcaption>Fig. 16. Reconstruction of a smooth tangential level-set.</figcaption></figure>

## 3.2. Deflection angle

The use of SIF-based criteria in cohesive models has been suggested by [80,81]. Authors explain that this is a reasonable assumption for concrete, since experiments showed that the crack path, contrary to the load-deflection curve, is almost insensitive to the size of the cohesive zone, to such a point that the crack path may still reasonably be predicted with LEFM assumptions.

The stress intensity factors may be post-processed from a separate LEFM computation – cohesive forces are removed from the model – as in [82,62]. However, this adds complexity to the procedure by introducing a LEFM stage in addition to the cohesive one.

For this reason, Moës and Belytschko [60] used the SIF-based criteria by computing SIF straight from the cohesive results. The reason why this can be done is brought by Planas and Elices [83], who explained that mechanical fields asymptotically tend to those of an equivalent free crack far enough from the process zone. Hence, equivalent stress intensity factors may be defined by contour integrals away from the process zone.

A virtual crack extension field $\theta$ is defined as tangent to the cohesive zone, directed towards $t$ on the front and so that $|\theta| = 1$ (Fig. 17). Let us consider a contour $C$ surrounding the cohesive zone, and denote by $\Gamma_C^+$ and $\Gamma_C^-$ the crack lip segments between the extremities of $C$ and the tip (see Fig. 17), so that the whole defines a domain $D_C$ (see Fig. 17). A generalized Rice integral is defined as:

$$J := - \int_{\partial D_C} \theta \cdot E \cdot n \, d\Gamma \qquad (28)$$
<figure><img src="image_22.png" /><figcaption>Fig. 17. Defining invariant integrals in the presence of cohesive forces.</figcaption></figure>

where $E = \nabla u^T \cdot \sigma - \frac{1}{2} (\sigma : \epsilon) \mathbf{1}$ is Eshelby's tensor and $\mathbf{n}$ is the outward normal to the contour. Due to the fact that asymptotic cohesive near-fields are not singular, we have $J = 0$.

By $\partial D_C = C \cup \Gamma_C^+ \cup \Gamma_C^-$, $J$ may be split up into $J = J_{\text{coh}} + J_{\text{ext}}$, with $J_{\text{ext}} := -\int_C \theta \cdot E \cdot \mathbf{n} \, d\Gamma$ and $J_{\text{coh}} = \int_{\Gamma_C} \theta \cdot [E] \cdot \mathbf{n} \, d\Gamma$ (see Rice [68]) with $[E] = E^+ - E^-$. Given the definition for $E$, the fact that $\theta \cdot \mathbf{n} = 0$ and the fact that $\sigma \cdot \mathbf{n} = \mathbf{t}_c$ on $\Gamma$, it comes:

$$J_{\text{coh}} = \int_{\Gamma_C} \mathbf{t}_c \cdot \nabla[\mathbf{u}] \cdot \boldsymbol{\theta} \, d\Gamma. \qquad (29)$$

It is noticeable that the integrand of (29) vanishes outside the cohesive zone, since $\mathbf{t}_c = 0$ on the traction-free zone and $[\mathbf{u}] = 0$ on the adherent zone. Hence $J_{\text{coh}}$ is independent of contour $C$ as long as it surrounds the cohesive zone (see Fig. 17), and may be written as:

$$J_{\text{ext}} = -J_{\text{coh}} = -\int_{\Gamma} \mathbf{t}_c \cdot \nabla[\mathbf{u}] \cdot \boldsymbol{\theta} \, d\Gamma. \qquad (30)$$

Consequently, $J_{\text{ext}}$ does not depend upon $C$ for such contours. It may then be interpreted as an intensity measure for an equivalent singularity, when looking far fields.

From Irwin's formula, it comes that $J_{\text{ext}} = -J_{\text{coh}} = \frac{1-\nu^2}{E} (K_{I,eq}^2 + K_{II,eq}^2) + \frac{1}{2\mu} K_{III,eq}^2$. Besides, the cohesive traction is decomposed as $\mathbf{t}_c = t_{c,n}\mathbf{n} + t_{c,t}\mathbf{t} + t_{c,b}\mathbf{b}$. Making use of the notation $[\nabla\mathbf{u}] \cdot \boldsymbol{\theta} = \frac{\partial[\mathbf{u}]_n}{\partial\theta}\mathbf{n} + \frac{\partial[\mathbf{u}]_t}{\partial\theta}\mathbf{t} + \frac{\partial[\mathbf{u}]_b}{\partial\theta}\mathbf{b}$, the equivalent SIF may alternatively be computed by:

$$K_{I,eq}^2 = -\frac{E}{1-\nu^2} \int_{\Gamma} \frac{\partial[\mathbf{u}]_n}{\partial\theta} t_{c,n} d\Gamma \qquad (31)$$

$$K_{II,eq}^2 = -\frac{E}{1-\nu^2} \int_{\Gamma} \frac{\partial[\mathbf{u}_t]}{\partial\theta} t_{c,t} d\Gamma \qquad (32)$$

$$K_{III,eq}^2 = -2\mu \int_{\Gamma} \frac{\partial[\mathbf{u}_b]}{\partial\theta} t_{c,b} d\Gamma. \qquad (33)$$

Physically, expressions (31)–(33) quantify the dissipated energy in the three fracture modes for a homothetic propagation of the crack in the direction $\boldsymbol{\theta}$, provided the cohesive process zone is small compared to the sample size. There is therefore no need to construct contours surrounding the front and no need to compute auxiliary fields with this method.

The adopted deflection angle is then the classical maximum hoop stress criterion defined by Erdogan and Sih [63]:

$$\beta = 2 \arctan \left[ \frac{1}{4} \left( K_{I,eq} / K_{II,eq} - \operatorname{sign}(K_{II,eq}) \sqrt{(K_{I,eq} / K_{II,eq})^2 + 8} \right) \right]. \qquad (34)$$
a

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

b

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

c

<figure><img src="image_25.png" /><figcaption>Fig. 18. Geometric update algorithm, as introduced by Colombo [32].</figcaption></figure>

### 3.3. Update of the potential crack surface

At this stage, the potential crack surface should be extended from the crack front after the deflection angle determined at the previous step. In other terms, the cracked zone – which contains the cohesive zone and the traction-free zone – is kept unchanged while the adherent zone is $\langle\langle \text{forgotten} \rangle\rangle$ and repositioned according to the deflection angle. In practice, this means applying a level-set update algorithm to the normal level-set (the tangential level-set having already been updated during the detection phase of Section 3.1).

To this aim, simple geometrical algorithms have emerged in the literature that take advantage of the fact that the previous cracked surface is frozen, to such a point that the use of such procedures is now prevailing for crack propagation problems. One of them was adopted here, named geometrical level-set update by [32]. New level-set values are computed from the previous crack front and advance vector field as follows:

For every node $M$ of the mesh:

* node $M$ is projected onto the previous crack front (see Fig. 18(a));

* its projection $Q$ onto the new front is deduced from the advance vector field (see Fig. 18(b));

* the updated normal (for this stage) or tangential level-sets (for the crack front detection of Section 3.1) is directly computed, respectively as the normal and tangential components of vector $QM$. In the case where we had $\phi_t(M) \le 0$ before the update of the crack front, $\phi_n(M)$ is not updated, so that the previous cracked surface is kept unchanged (see Fig. 18(c)).

### 3.4. Extension of the multiplier space and initial internal variables

As illustrated by Fig. 19(a)–(b), the update of the crack surface implies that it has turned ahead of the front. Consequently, the edges which are intersected by the new crack surface are not the same as those which were intersected by the former crack surface, especially ahead of the crack front (see Fig. 19(a)–(b)).

In Fig. 19(d), $Q$ denotes the set of intersected edges by the new crack surface, $V$ the subset of edges with an equality relation, and $K$ denotes the set of enriched nodes. We denote $Q_0$, $V_0$ and $K_0$ these (different) sets if the former crack surface configuration is considered (Fig. 19(a)). Hence, the construction of the set of equality relations $V$ – defining the new reduced multiplier space $M_h$ – non longer relies upon $Q_0$: it has to be formed again starting from $Q$.

Now, a new $V$ started from scratch would possibly lead to completely different combinations of nodes sharing a Lagrange degree of freedom. As the internal variables were defined on the older $M_h$, they would have to be projected onto a quite different space, which would involve serious energetic issues. It is therefore wiser to start the new space abiding by preexisting combinations behind the crack front, and extend it with new groups on the formerly non intersected area, located ahead of the crack front (see Fig. 19(a)–(c)).

In practice, edges of $Q_0 \setminus V_0$ whose both vertices belong to $K$ are removed from the set $Q$ provided as an input to the restriction algorithm (see Fig. 19(b)). The restriction algorithm by Géniaut et al. [77] is then performed, thus naturally resulting in a $V$ which preserves preexisting combinations (see Fig. 19(c)).

The new restricted space involves nodes in $K \setminus K_0$ which had no attributed value for the internal variables (white nodes in Fig. 19(c)). For such a node $n$, an initial value is determined as follows (see Fig. 19(c)–(d)):

* if there exists an edge in $V$ connecting $n$ to a node $m \in K_0$, the value of $m$ is attributed to $n$;

* otherwise, the initial internal variable is set to 0 (uncracked material).
<figure><img src="image_26.png" /><figcaption>Fig. 19. Extension of the restricted multiplier space and update of the internal variables.</figcaption></figure>

Note that such update for internal variables leads to only slight equilibrium discrepancy, since implied nodes are very close to the crack front or ahead of it, and hence have almost zero values (see Fig. 19(a)).

There are also a few nodes in $K_0$ whose value must still be changed, since the related connecting edge has changed between $\mathcal{V}$ and $\mathcal{V}_0$ (grey nodes in Fig. 19(c)). Once more if for such a node, there exists an edge in $\mathcal{V}$ connecting it to a node $m \in K_0$, the value of $m$ is attributed to $n$.

# 4. Numerical tests

Several tests of increased complexity were carried out, which all have been numerically and experimentally investigated in earlier literature. Whenever possible, these benchmarks have been carried out with the same data.

## 4.1. An extruded test: the L-shaped panel

Based on the experimental results by Wrinkler [84] on plain concrete, this L-shaped panel test was numerically reproduced by [48,62,71]. The geometry and loading are summarized in Fig. 21. The documented material data of the
<figure><img src="image_27.png" /><figcaption>Fig. 20. Bilinear-softening law.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 21. L-shaped panel: geometry and loading.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 22. Comparison of computed and experimental crack paths for the L-shaped panel.</figcaption></figure>

experimented concrete are given in Table 2, as well as the values reported by [48] as best representing the results, also used here for computation. Instead of linear softening, the bilinear-softening law of Fig. 20 is used, which is closer to the actual fracture properties of concrete.

As shown by Fig. 22, the computed crack path appears to be within the experimental range. As for the load-deflection curve, represented in Fig. 23, the peak is rather accurately reproduced, but the computed post-peak
<table><caption>Table 2
Material data for the plain concrete of the L-shaped panel.</caption><thead><tr><th></th><th>Experimented concrete</th><th>Used for computation (see [48])</th></tr></thead><tbody><tr><td>Modulus of elasticity</td><td><i>E</i> = 25.85 GPa</td><td><i>E</i> = 25 GPa</td></tr><tr><td>Poisson's ratio</td><td><i>ν</i> = 0.18</td><td>ibid</td></tr><tr><td>Tensile strength</td><td><i>σ</i><sub>c</sub> = 2.7 MPa</td><td><i>σ</i><sub>c</sub> = 2.5 MPa</td></tr><tr><td>Fracture energy</td><td><i>G</i><sub>c</sub> = 0.095 N mm<sup>-1</sup></td><td><i>G</i><sub>c</sub> = 0.13 N mm<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_30.png" /><figcaption>Fig. 23. Computed and experimental load-deflection curves.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 24. L-shaped panel: deformed shape and stress modulus.</figcaption></figure>

behaviour exhibits a steeper slope than the experiments. This is probably due to the use of bilinear instead of exponential softening (see [38]), the latter being closer to the actual softening behaviour of concrete [55]. As expected from the values of Table 2, the process zone is of the same order than the characteristic dimension of the specimen, as can be seen on the map of normal cohesive tractions on the crack surface (see Fig. 25). The stress modulus map is represented in Fig. 24.

## 4.2. Three-point bending test with an initial skew crack

This test consists of an experimental and numerical study of the propagation of a fatigue crack in PMMA (Plexiglas®) by [85,86], and was used as a numerical benchmark by [32]. Geometry and loading conditions are
<figure><img src="image_32.png" /><figcaption>Fig. 25. L-shaped panel: profile of the normal cohesive traction on the crack surface.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Fig. 26. Three-point-bending test with an initial skew crack.</figcaption></figure>

<table><caption>Table 3
Material data for the PMMA sample.</caption><thead><tr><td></td><td>Physical values</td><td>Used for computation</td></tr></thead><tbody><tr><td>Modulus of elasticity</td><td><i>E</i> = 2.8 GPa</td><td>ibid</td></tr><tr><td>Poisson's ratio</td><td><i>ν</i> = 0.38</td><td>ibid</td></tr><tr><td>Fracture energy</td><td><i>G</i><sub>c</sub> = 0.5 N mm<sup>-1</sup></td><td>ibid</td></tr><tr><td>Tensile strength</td><td><i>σ</i><sub>c</sub> = 40 MPa</td><td><i>σ</i><sub>c</sub> = 15 MPa</td></tr></tbody></table>

summarized on Fig. 26. The actual material parameters are given in Table 3. Here, the applied load is monotonic instead of cyclic, and focus is made on an accurate prediction of the crack path. The typical length of the cohesive zone is then expected to be (see [87,88]) $l_c = \frac{EG_c}{(1-ν^2)σ_c^2} = 0.3$ mm. As such the process zone is much smaller than a characteristic dimension of the sample. The response of the structure should tend to that given by LEFM, and a larger size of the process zone may even be used for computation without significant change of results (see Table 3) as long as it remains small compared to the specimen size. Consequently, as the crack path is controlled by the stress intensity factors ratio, it is expected to be identical for monotonic and cyclic loading: this is why our computed crack path in the monotonic case should compare well to fatigue experiments for this specific case.

The computed crack path (see Figs. 27–28) is in good accordance with the experiment (Fig. 29), with an initial twist of the crack path to recover mode I loading conditions (see Figs. 27–28): the qualitative evolution of the crack front is similar to that of the experiments (see Fig. 30). Fig. 31 shows the comparison of the positions of the endpoint of the crack in the cut plane $x = -5$ mm.

### 4.3. Brokenshire's torsion test

Again, a prismatic specimen with an initial skew crack is considered, and submitted to torsion (see Fig. 32). Experiments (see Figs. 36 and 33) were carried out by Brokenshire [89], and the test was numerically reproduced
<figure><img src="image_34.png" /><figcaption>Fig. 27. Perspective view of the crack surface and field of cohesive tractions.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Fig. 28. Top view of the crack surface.</figcaption></figure>

by [34,54], as a benchmark to test the ability of several crack propagation algorithms to capture non-planar crack paths. Geometry and loading are represented in Fig. 32. The material parameters of plain concrete of Table 4 and the bilinear-softening law of Fig. 20 are adopted.

The experimental «S-shaped» crack path is accurately reproduced for this test, as can be seen in Fig. 34. Plotting the applied load as a function of the crack mouth opening displacement (CMOD) at the centre top of the notch (see Fig. 32), the peak is accurately reproduced while the post-peak behaviour appears to be too stiff again (see Fig. 33),
<figure><img src="image_36.png" /><figcaption>Fig. 29. Experimental crack path from [85,86].</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Fig. 30. 3-point bending test with an initial skew crack: deformed shape and map of stress modulus.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Fig. 31. Crack path in the cut plane ($X = -5$).</figcaption></figure>

which is still due to the use of bilinear instead of exponential softening (see [38]). This time, the size of the process zone is large compared to the size of the specimen, as can be seen on the profile of normal cohesive tractions on the
<figure><img src="image_39.png" /><figcaption>Fig. 32. Geometry and loading for Brokenshire's torsion test.</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Fig. 33. Load-deflection curves for Brokenshire's torsion test.</figcaption></figure>

<table><caption>Table 4</caption><caption>Material data for the plain concrete of Brokenshire's test.</caption><tbody><tr><td>Modulus of elasticity</td><td><i>E</i> = 35 GPa</td></tr><tr><td>Poisson's ratio</td><td><i>ν</i> = 0.2</td></tr><tr><td>Fracture energy</td><td><i>G</i><sub>c</sub> = 0.0825 N mm<sup>-1</sup></td></tr><tr><td>Tensile strength</td><td><i>σ</i><sub>c</sub> = 2.3 MPa</td></tr></tbody></table>

crack surface (see Fig. 34). In spite of this, the crack path is still rather well predicted with LEFM assumptions. The stress modulus map is represented on Fig. 35.

## 5. Conclusion

In this paper, an original method is presented to predict crack paths and load-deflection curves of brittle or quasi-brittle materials, with a combined use of the X-FEM and cohesive zone models. Instead of specifying the location of the crack front for the next load increment in advance, the cohesive law is defined over a broad interface, so that the next location of the crack front is naturally embedded in the equilibrium. This demands a robust insertion of cohesive laws with an exact enforcement of initial adherence in the X-FEM. The key ideas to achieve this are:
<figure><img src="image_41.png" /><figcaption>Fig. 34. Shape of the computed crack path and map of normal cohesive tractions for Brokenshire’s test.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>Fig. 35. Deformed shape and map of stress modulus for Brokenshire’s test.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>Fig. 36. Experimental crack path from [89].</figcaption></figure>

* an expression of the cohesive law within the augmented Lagrangian formalism,

* the use of a reduced space for the discrete cohesive traction and displacement jump,

* a cohesive law that is written componentwise on the reduced space.
We think that the methodology could be extended to most non-smooth interfacial behaviour (unilateral contact for instance) with on-going work in that direction.

The crack kink angle is determined by classical criteria deduced from equivalent stress intensity factors, computed creatively by surface integrals over the cohesive zone. Three numerical tests were carried out, for which there exist experimental data and computation with other methods. The method is efficient to track the crack path and predict the effect of the size of the sample over the load-deflection curve.

## References

[1] C. Comi, U. Perego, Numerical aspects of non-local damage analyses, Rev. Eur. Elem. Finis 10 (2001) 227-242.

[2] G. Pijaudier-Cabot, Z. Bazant, Nonlocal damage theory, J. Eng. Mech. 113 (1987) 1512-1533.

[3] E. Lorentz, S. Andrieux, Analysis of non-local models through energetic formulations, Int. J. Solids Struct. 40 (2003) 2905-2936.

[4] Q. Nguyen, S. Andrieux, The non-local generalized standard approach: a consistent gradient theory, C. R. Méc. 33 (2005) 139-145.

[5] R. Bargellini, J. Besson, E. Lorentz, S. Michel-Ponnelle, A non local finite element based on volumetric strain gradient: application to ductile fracture, Comput. Mater. Sci. 33 (2009) 762-767.

[6] K. Pham, J.-J. Marigo, Approche variationnelle de l'endommagement: ii. les modèles à gradient, C. R. Méc. 338 (2010) 199-206.

[7] R. Fernandes, C. Chavant, R. Chambon, A simplified second gradient model for dilatant materials: theory and numerical implementation, Int. J. Solids Struct. 45 (2008) 5289-5307.

[8] R. Plassart, R. Fernandes, A. Giraud, D. Hoxha, F. Laigle, Hydromechanical modelling of an excavation in an underground research laboratory with an elasto-visco-plastic behaviour law and regularization by second gradient of dilation, Int. J. Rock Mech. Min. Sci. 58 (2013) 23-33.

[9] V. Hakim, A. Karma, Laws of crack motion and phase-field models of fracture, J. Mech. Phys. Solids 57 (2009) 342-368.

[10] C. Miehe, F. Welschinger, Hofhacker, Thermodynamically consistent phase field models of fracture: variational principles and multifields fe implementations, Internat. J. Numer. Methods Engrg. 83 (2010) 1273-1311.

[11] G. Barenblatt, The mathematical theory of equilibrium of cracks in brittle failure, Adv. Appl. Mech. 7 (1962) 55-129.

[12] H. Westergaard, Exact solution for a cracked elastic bulk: bearing pressures and cracks, J. Appl. Mech. 6 (1939) 49-53.

[13] S. Géniaut, E. Galenne, A simple method for crack growth in mixed mode with x-fem, Int. J. Solids Struct. 49 (2012) 2094-2106.

[14] N. Moës, C. Stolz, P. Bernard, N. Chevaugeon, A level set based model for damaged growth: the thick level-set approach, Internat. J. Numer. Methods Engrg. 86 (2011) 358-380.

[15] P. Bernard, N. Moës, N. Chevaugeon, Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings, Comput. Methods Appl. Mech. Engrg. 233-236 (2012) 11-27.

[16] C. Stolz, N. Moës, A new model of damage: a moving thick layer approach, Int. J. Fract. 174 (2012) 49-60.

[17] N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Internat. J. Numer. Methods Engrg. 46 (1999) 135-150.

[18] J.M. Melenk, I. Babuška, The partition of unity finite element method: basic theory and applications, Comput. Methods Appl. Mech. Engrg. 39 (1996) 289-314.

[19] S. Osher, J. Sethian, Fronts propagating with curvature-dependent speed: algorithm based on Hamiton-Jacobi formulations, J. Comput. Phys. 79 (1988) 12-49.

[20] N. Sukumar, D. Chopp, N. Moës, T. Belytschko, Modeling holes and inclusions by level-sets in the extended finite element method, Comput. Methods Appl. Mech. Engrg. 190 (2001) 6183-6200.

[21] M. Stolarska, D. Chopp, N. Moës, T. Belytschko, Modelling crack growth by level sets in the extended finite element method, Internat. J. Numer. Methods Engrg. 51 (2001) 943-960.

[22] D. Peng, B. Merriman, S. Osher, H. Zhao, M. Kang, A pde-based fast local level-set method, J. Comput. Phys. 155 (1999) 410-438.

[23] N. Moës, A. Gravouil, T. Belytschko, Non-planar 3d crack growth by the extended finite element and level-sets-part ii : level-set update, Internat. J. Numer. Methods Engrg. 53 (2002) 2569-2586.

[24] B. Prabel, Modélisation avec la méthode X-FEM de la propagation dynamique et de l'arrêt de fissure de clivage dans un acier de cuve REP (Ph.D. thesis), INSA de Lyon, 2007.

[25] D. Colombo, P. Massin, Fast and robust level-set update for 3d non-planar x-fem crack propagation modelling, Comput. Methods Appl. Mech. Engrg. 200 (2011) 2160-2180.

[26] B. Prabel, A. Combescure, A. Gravouil, S. Marie, Level set x-fem non-matching meshes: application to dynamic crack propagation in elastic-plastic media, Internat. J. Numer. Methods Engrg. 155 (2006) 410-438.

[27] T. Barth, J. Sethian, Numerical schemes for the Hamilton-Jacobi and level-set equations on triangulated domains, J. Comput. Phys. 145 (1998) 1-40.

[28] N. Sukumar, D. Chopp, B. Moran, Extended finite element method and fast marching method for three dimensional fatigue crack propagation, Eng. Fract. Mech. 70 (2003) 29-48.

[29] N. Sukumar, D. Chopp, E. Béchet, N. Moës, Three-dimensional non-planar crack growth by a coupled extended finite element and fast marching method, Internat. J. Numer. Methods Engrg. 76 (2008) 727-748.

[30] M. Duflot, A study of the representation of cracks with level-sets, Internat. J. Numer. Methods Engrg. 70 (2007) 1261-1302.

[31] T. Menouillard, J. Réthoré, N. Moës, A. Combescure, H. Bung, Mass lumping strategies for x-fem explicit dynamics: Application to crack propagation, Internat. J. Numer. Methods Engrg. 74 (2008) 447-474.

[32] D. Colombo, An implicit geometrical approach to level-set update for 3d non-planar x-fem crack propagation, Comput. Methods Appl. Mech. Engrg. 240 (2012) 39-50.
[33] G. Ventura, E. Budyн, T. Belytschko, Vector level sets for description of propagating cracks in finite elements, Internat. J. Numer. Methods Engrg. 58 (2003) 1571–1592.

[34] T. Fries, M. Baydoun, Crack propagation with the extended finite element and a hybrid explicit-implicit crack description, Internat. J. Numer. Methods Engrg. 89 (2012) 1527–1558.

[35] L. Kaczmarczyk, M.M. Nezhad, C. Pearse, Three-dimensional brittle fracture: configurational-force-driven crack propagation, Internat. J. Numer. Methods Engrg. 97 (2014) 531–550.

[36] P. Destyunder, M. Djaoua, S. Lescure, Quelques remarques sur la mécanique de la rupture élastique, J. Méc. Théor. Appl. 2 (1983) 113–135.

[37] T. Hellen, On the method of virtual crack extentions, Internat. J. Numer. Methods Engrg. 9 (1975) 187–207.

[38] G. Alfano, On the influence of the shape of the cohesive law on the application of cohesive zone models, Compos. Sci. Technol. 66 (2006) 723–730.

[39] P. Bouchard, F. Bay, Y. Chastel, I. Tovena, Crack propagation modelling using an advanced remeshing technique, Comput. Methods Appl. Mech. Engrg. 189 (2000) 723–742.

[40] R. Barsoum, Application of quadratic isoparametric finite elements in linear fracture mechanics, Int. J. Fract. 10 (1974) 603–605.

[41] L. Banks-Sills, I. Hershkovitz, P.A. Wawrzynek, R. Eliasi, A.R. Ingraffea, Methods for calculating stress intensity factors in anisotropic materials: Part i—z = 0 is a symmetric plane, Eng. Fract. Mech. 72 (2005) 2328–2358.

[42] L. Banks-Sills, P.A. Wawrzynek, B. Carter, A.R. Ingraffea, I. Hershkovitz, Methods for calculating stress intensity factors in anisotropic materials: Part ii—arbitrary geometry, Eng. Fract. Mech. 74 (2007) 1293–1307.

[43] Z. Wen-bin, Z. Xue-xia, C. Xiao-chao, Y. Wei-yang, Analysis of stress intensity factor in orthotropic bi-material mixed interface crack, Appl. Math. Mech. 35 (2014) 1271–1292.

[44] R.S. Alwar, S. Thiagarajan, Combined effect of shear and large deformation on stress intensity factors, Internat. J. Numer. Methods Engrg. 28 (1989) 1951–1964.

[45] J. Hutchinson, Singular behaviour at the end of a tensile crack in a hardening material, J. Mech. Phys. Solids 16 (1968) 13–31.

[46] J. Rice, G. Rosengren, Plane strain deformation near a crack tip in a power-law hardening material, J. Mech. Phys. Solids 16 (1968) 1–12.

[47] J. Hutchinson, Plastic stress and strain fields at a crack tip, J. Mech. Phys. Solids 16 (1968) 337–342.

[48] J. Unger, S. Eckardt, C. Könke, Modelling of cohesive crack growth in concrete structures with the extended finite element method, Comput. Methods Appl. Mech. Engrg. 196 (2007) 4087–4100.

[49] P. Dumstorff, G. Meschke, Crack propagation criteria in the framework of x-fem-based structural analyses, Int. J. Numer. Anal. Methods Geomech. 31 (2007) 239–259.

[50] G. Wells, L. Sluys, A new method for modelling cohesive cracks using finite elements, Internat. J. Numer. Methods Engrg. 50 (2001) 2667–2682.

[51] J. Merghem, E. Kuhl, P. Steinmann, A finite element method for the computational modelling of cohesive cracks, Internat. J. Numer. Methods Engrg. 63 (2005) 276–289.

[52] P. Jäger, P. Steinmann, E. Kuhl, Modelling three dimensional crack propagation—a comparison of crack path tracking strategies, Internat. J. Numer. Methods Engrg. 76 (2008) 1328–1352.

[53] P. Areias, T. Belytschko, Analysis of three-dimensional crack initiation and propagation using the extended finite element method, Internat. J. Numer. Methods Engrg. 63 (2005) 760–788.

[54] T. Gasser, G. Holzapfel, 3d crack propagation in unreinforced concrete. A two-step algorithm for tracking 3d crack paths, Comput. Methods Appl. Mech. Engrg. 195 (2006) 5198–5219.

[55] T. Gasser, G. Holzapfel, Modelling 3d crack propagation in unreinforced concrete using pufem, Comput. Methods Appl. Mech. Engrg. 194 (2005) 2859–2896.

[56] J. Remmers, R. de Borst, A. Needleman, A cohesive segments method for the simulation of crack growth, Comput. Mech. 31 (2003) 69–77.

[57] R. de Borst, J. Remmers, A. Needleman, Mesh-independent discrete numerical representations of cohesive zone models, Eng. Fract. Mech. 73 (2006) 160–177.

[58] J. Oliver, A. Huespe, On strategies for tracking strong discontinuities in computational failure mechanics, in: *Proceeding of the 5th WCCM, Vienna, Austria, 2002*.

[59] J. Oliver, A. Huespe, E. Samaniego, E. Chaves, Continuum approach to the numerical simulation of material failure in concrete, Int. J. Numer. Anal. Methods Geomech. 28 (2004) 609–632.

[60] N. Moës, T. Belytschko, Extended finite element method for cohesive crack growth, Eng. Fract. Mech. 69 (2002) 813–833.

[61] G. Zi, T. Belytschko, New crack-tip elements for x-fem and application to cohesive cracks, Internat. J. Numer. Methods Engrg. 57 (2003) 2221–2240.

[62] A. Zamani, R. Gracie, M. Eslami, Cohesive and non-cohesive fracture by higher-order enrichment of x-fem, Internat. J. Numer. Methods Engrg. 90 (2012) 452–483.

[63] F. Erdogan, G. Sih, On the crack extension in plane loading and transverse shear, J. Basic Engng. 85 (1963) 519–527.

[64] G. Sih, Mechanics of Fracture Initiation and Propagation, vol. 85, Kluwer Academic Publisher, Boston, 1991.

[65] M. Hussain, S. Pu, J. Underwood, Strain energy release rate for a crack under combined mode i and mode ii, Fract. Anal. ASTM STP 560 (1974) 2–28.

[66] J. Chang, J. Xu, Y. Mutoh, A general mixed-mode brittle fracture criterion for cracked materials, EFM 73 (2006) 1246–1263.

[67] N. Hasebe, S.N. Nasser, L. Keer, Stress analysis of a kinked crack initiating from a rigid line inclusion. Part II: direction of propagation, Mech. Mater. 3 (1984) 147–156.

[68] J. Rice, Elastic fracture mechanics concepts for interfacial cracks, J. Appl. Mech. 55 (1988) 98–103.

[69] D. Haboussa, Modélisation de la transition traction cisaillement des métaux sous choc par la X-FEM (Ph.D. thesis), INSA de Lyon, 2012.

[70] M. Schöllman, H. Richard, G. Kullmer, M. Fulland, A new criterion for the prediction of crack development in multiaxially loaded structures, Int. J. Fract. 117 (2002) 129–141.
[71] G. Meschke, P. Dumstorff, Energy based modelling of cohesive and cohesionless cracks via x-fem, Comput. Methods Appl. Mech. Engrg. 196 (2007) 2338-2357.

[72] Q. Xiao, B. Karihaloo, Asymptotic fields at frictionless and frictional cohesive crack tips in quasibrittle materials, J. Mech. Mater. Struct. 1 (2006) 881-910.

[73] H. Ji, J. Dolbow, On strategies for enforcing interfacial constraints and evaluating jump conditions with the extended finite element method, Internat. J. Numer. Methods Engrg. 61 (2004) 2508-2535.

[74] E. Béchet, N. Moës, B. Wohlmuth, A stable lagrange multiplier space for stiff interface conditions within the extended finite element method, Internat. J. Numer. Methods Engrg. 78 (2009) 931-954.

[75] N. Moës, E. Béchet, M. Tourbier, Imposing dirichlet boundary conditions in the extended finite element method, Internat. J. Numer. Methods Engrg. 67 (2006) 1641-1669.

[76] J. Dolbow, L. Franca, Residual-free bubbles for embedded dirichlet problems, Comput. Methods Appl. Mech. Engrg. 197 (2008) 3751-3759.

[77] S. Géniaut, P. Massin, N. Moës, A stable 3d contact formulation for cracks using x-fem, Rev. Eur. Méc. Numér. 75 (2007) 259-276.

[78] M. Hautefeuille, C. Annavarapu, J. Dolbow, Robust imposition of dirichlet boundary conditions on embedded surfaces, Internat. J. Numer. Methods Engrg. 90 (2012) 40-64.

[79] S. Amdouni, P. Hild, V. Lleras, M. Moakher, Y. Renard, A stabilized lagrange multiplier method for the enriched finite-element approximation of contact problems of cracked elastic bodies, Math. Modelling Numer. Anal. 49 (2012) 813-839.

[80] P. Bocca, A. Carpitieri, S. Valente, Mixed-mode fracture of concrete, Int. J. Solids Struct. 27 (1991) 1139-1153.

[81] D. Cendon, J. Galvez, M. Elices, J. Planas, Modeling the fracture of concrete under mixed loading, Int. J. Fract. 103 (2000) 293-310.

[82] G. Meschke, P. Dumstorff, W. Fleming, S. Jox, Computational failure analysis of concrete structures using the extended finite element method, Neue Bauweisen-Trends in Statik und Dynamik 2 (2006) 395-408.

[83] J. Planas, M. Elices, Asymptotic analysis of a cohesive crack: 1. Theoretical background, Int. J. Fract. 55 (1992) 153-177.

[84] B. Wrinkler, Traglastuntersuchungen von Unbewehrten und Bewehrten Betonstrukturen auf der Grundlage eines objectiven Werkstoffgesetzes für Beton (Ph.D. thesis), Universität Innsbruck, 2001.

[85] R. Citarella, F. Buchholz, Comparison of crack growth simulation by dbem and fem for sen-specimen undergoing torsion or bending loading, Eng. Fract. Mech. 75 (2008) 489-509.

[86] V. Lazarus, F.-G. Buchholz, M. Fulland, J. Wiebesiek, Comparison of predictions by mode ii or mode iii criteria on crack front twisting in three or four point bending experiments, Int. J. Fract. 153 (2008) 141-151.

[87] E. Lorentz, V. Godard, Gradient damage models: Toward full-scale computations, Comput. Methods Appl. Mech. Engrg. 200 (2011) 1927-1944.

[88] J. Rice, The mechanics of earthquake rupture, Phys. Earth Inter. (1980) 555-649.

[89] D. Brokenshire, B. Barr, A study of torsion fracture tests (Ph.D. thesis), Cardiff's University, 1996.