# A step-by-step global crack-tracking approach in E-FEM simulations of quasi-brittle materials

Francesco Riccardi, Ejona Kishta, Benjamin Richard

► To cite this version:

Francesco Riccardi, Ejona Kishta, Benjamin Richard. A step-by-step global crack-tracking approach in E-FEM simulations of quasi-brittle materials. Engineering Fracture Mechanics, 2017, 170, pp.44-58. 10.1016/j.engfracmech.2016.11.032. cea-02381488

HAL Id: cea-02381488

https://cea.hal.science/cea-02381488

Submitted on 29 Jan 2020

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
# A step-by-step global crack-tracking approach in E-FEM simulations of quasi-brittle materials

Francesco Riccardi<sup>a</sup>, Ejona Kishta<sup>a,b</sup>, Benjamin Richard<sup>a,*</sup>

The numerical simulation of two-dimensional fracture processes of quasi-brittle materials by means of the Embedded Finite Element Method is dealt with. Attention is paid to the coupling with the global crack-tracking strategy, which has been proposed in the literature in the form of a heat conduction-like problem. It turns out that the stiffness-like matrix associated with this formulation is singular and a numerical perturbation has to be intro-duced in order to overcome the ill-posedness of the problem. The dependence of the solution on this parameter may represent a limitation for the global tracking approach. Furthermore, it is found that if the root of each discontinuity is not updated during an incremental analysis, a loss of continuity of the crack path may appear when principal stress directions rotate. This paper aims to provide a solution to the aforementioned issues. An alternative mathematical formulation of the problem is thus given in terms of Navier-Stokes equations, linking the diffusive contribution to a characteristic mesh length. Additionally, a modified crack-tracking algorithm, considering the evolution of the root for the identification of the crack path, is proposed. The numerical assessment of the pro-posed tracking strategy is reported by means of benchmark tests at the structural level.

## 1. Introduction

In the last decade, the Embedded Finite Element Method (E-FEM) has gained wide popularity for the description of cracking phenomena [1–6]. Due to the local nature of the kinematic enrichment, this approach presents some advantages concerning the computational effort with respect, for example, to the Extended Finite Element Method (X-FEM) [7–11]. Indeed the additional degrees of freedom can be statically condensed, since they are expressed at the element level in terms of crack opening displacement components [5], thus leaving unchanged the dimension of the global system of equations. As a counterpart, since no information is available in the vicinity of cracked elements, the continuity of the crack path is not intrinsically guaranteed and a loss in objectivity may be encountered in numerical simulations [12,13]. For this reason, the E-FEM is usually associated with the adoption of tracking algorithms leading to C⁰-continuous crack paths [14–17]. Different crack-tracking strategies can be found in the literature. In general, two families may be distinguished: local tracking and global tracking algorithms [14,5,18,19].

Local tracking relies on geometry-based or energy-based schemes, which are applied for each element able to crack. The main difference between the two strategies stands in the evaluation of the direction of the propagating discontinuity. In case

<footnote><sup>a</sup> CEA, DEN, DANS, DM2S, SEMT, Laboratoire d'Études de Mécanique Sismique, F-91191 Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup> LMT, ENS Cachan, CNRS, Université Paris-Saclay, 61 avenue du Président Wilson, F-94235 Cachan, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: benjamin.richard@cea.fr (B. Richard).</footnote>
Nomenclature

<table><tr><td>Ω</td><td>bi-dimensional domain</td></tr><tr><td>T</td><td>tangent vector field</td></tr><tr><td>N</td><td>normal vector field</td></tr><tr><td>θ</td><td>scalar temperature-like field</td></tr><tr><td>Γ</td><td>isovalue of the temperature-like field θ</td></tr><tr><td>q</td><td>conduction-like flux vector</td></tr><tr><td>K</td><td>conductivity-like tensor</td></tr><tr><td>ε</td><td>numerical conductivity-like parameter</td></tr><tr><td>α</td><td>scalar diffusion-like parameter</td></tr><tr><td>J</td><td>normal diffusion-like flux</td></tr><tr><td>v</td><td>unit outward normal vector</td></tr><tr><td>h<sub>T</sub></td><td>mesh characteristic length</td></tr><tr><td>E</td><td>Young's modulus</td></tr><tr><td>ν</td><td>Poisson's ratio</td></tr><tr><td>G<sub>f</sub></td><td>tensile fracture energy</td></tr><tr><td>f<sub>f</sub></td><td>uniaxial tensile strength</td></tr><tr><td>r</td><td>root of the propagating discontinuity</td></tr><tr><td>σ</td><td>Cauchy's stress tensor</td></tr><tr><td>I<sub>f</sub></td><td>crack input point</td></tr><tr><td>s</td><td>curvilinear abscissa</td></tr><tr><td>[u]</td><td>displacement jump</td></tr><tr><td>SEN</td><td>single-edge notched</td></tr><tr><td>CMOD</td><td>Crack Mouth Opening Displacement</td></tr><tr><td>CMOD</td><td>Crack Mouth Sliding Displacement</td></tr></table>

of a geometrical approach, the crack orientation is given by the assumed failure criterion, e.g. accordingly to Rankine it is supposed to be perpendicular to the maximum tensile principal direction [20,5,21], whereas in case of an energetic approach, it is computed from linear-elastic fracture mechanics (LEFM) by minimization of the mechanical energy [18,22]. From the knowledge of the root of each discontinuity, i.e. the material point experiencing failure and not associated to any pre-existing crack path, it is possible to make the input point of the crack inside an element match the output point of the crack present in an adjacent element. Local tracking techniques are able to reproduce continuous crack paths in a robust manner by exploiting the informations of nearby elements. Nevertheless, their implementation in the case of multiple crack problems may be cumbersome and this strategy can lose much of its robustness.

Global tracking has been introduced to overcome the limits of local strategies when dealing with multi-cracking modeling [14,5]. The main idea is to trace the envelopes of the tangent vector field to the discontinuities as the isovalues of the temperature field of a heat conduction-like problem in the case of steady state conditions and no internal sources. Dirichlet and Neumann boundary conditions must be prescribed on the respective portions of the boundary, in terms of fixed temperature values and (null) heat flux. Once the assumed failure criterion is fulfilled for the first time in a certain material point, the latter becomes the root of a new discontinuity line, which can be traced as the isovalue passing through that point. Numerical applications can be found in [23,18,24]. The main advantage of this approach is the fact that no information from the neighbourhood of the cracked elements is required to perform the analysis: indeed, since the isovalues are available at every point of the domain, only root element coordinates shall be provided in order to ensure a continuous crack path.

The finite element formulation of the heat conduction-like problem is straightforward. However, the stiffness-like matrix deduced from the anisotropic conductivity tensor reveals to be singular and a user-defined perturbation (isotropic algorithmic conductivity) is introduced in numerical simulations to circumvent this drawback. The dependence of the solution on this parameter may then represents a limitation for the application of global tracking, since its value changes for each specific structural problem. In addition, it is found out that the capability of the thermal-like isovalues to envelop the vector tangent field is reduced as soon as cracking occurs. This fact, due to the rotation of the principal stress directions outside the region crossed by the discontinuity, may lead to a loss of continuity of the crack path and, in the worst case, to a wrong evaluation of the enriched shape functions whenever the element domain is not decomposed properly [12,25]. Such a circumstance becomes increasingly critical as the evolution of the principal stress field is important, see Fig. 1.

The objective of this paper is to provide an alternative formulation of the global crack-tracking strategy able to improve the performance of this technique with respect to the aforementioned issues. This paper is organized as follows. In Section 2 a new physical interpretation of the problem is given in terms of Navier-Stokes equations, where the concept of numerical diffusion is introduced in order provide a stable and consistent solution of the initially ill-posed discrete problem. A revised algorithm for ensuring continuous crack paths in case of step-by-step analysis is presented in Section 3, considering the evolution of the root for the choice of the isovalue enveloping the propagating discontinuity. Section 4 investigates the
<figure><img src="image_2.png" /><figcaption>Fig. 1. Evolution of the thermal isovalue in a double-edge notched specimen under tension.</figcaption></figure>

application of the proposed model to the E-FEM by means of benchmark tests. Section 5 then concludes with a critical comparison between the proposed approach and the original formulation.

## 2. Global equations

The problem of tracing the envelopes of a vector field $T(x)$ in a domain $\Omega$ has been treated in [14,5]. For the sake of simplicity, let us focus upon the bi-dimensional case. If we indicate with $N(x)$ the normal vector field to $T(x)$, we can consider a scalar function $\theta(x)$ whose gradient is parallel to $N(x)$, i.e. such that:

$$ N(x) = \frac{\nabla \theta(x)}{\|\nabla \theta(x)\|}, \quad x \in \Omega \qquad (1) $$

Hence, the following partial differential equation will hold in the domain $\Omega$:

$$ T(x) \cdot \nabla \theta(x) = 0 \qquad (2) $$

Since the level contours of the function $\theta(x)$ are orthogonal to the gradient, the envelope of the vector field $T(x)$ passing through a generic point $P$ can be defined as:

$$ \Gamma_P = \{x \in \Omega \mid \theta(x) = \theta_P\} \qquad (3) $$

The envelopes of the vector field $T$ thus provide $C^0$-continuous curves, which are well-suited to model crack-paths within the framework of the E-FEM. From now on the dependence of all the quantities on $x$ will be omitted.

### 2.1. Heat conduction-like problem

Eq. (2) can be manipulated by multiplying it by the vector field $T$. After some analytical computations and using the same notations as in [14], the previous problem can be reformulated as the following boundary value problem for the unknown function $\theta$:

$$ \nabla \cdot q = 0 \quad \forall x \in \Omega \qquad (4a) $$

$$ q = -K \nabla \theta \quad \forall x \in \Omega \qquad (4b) $$

$$ q \cdot v = 0 \quad \forall x \in \partial_q \Omega \qquad (4c) $$

$$ \theta = \theta^* \quad \forall x \in \partial_0 \Omega \qquad (4d) $$

with:

$$ K := T \otimes T \qquad (5) $$

where $\otimes$ denotes the tensor product. The boundary value problem (4) defines a heat conduction-like problem in the domain $\Omega$, where $\theta$ is the temperature field and where $q$ is the conduction flux vector. Boundary conditions are prescribed on the boundary $\partial\Omega = \partial_q\Omega \cup \partial_0\Omega$ such that $\partial_q\Omega \cap \partial_0\Omega = \emptyset$. More precisely, Eq. (4c) expresses the Neumann condition of a null heat flux on the set $\partial_q\Omega$, whereas Eq. (4d) represents the Dirichlet condition of fixed temperature values on the boundary $\partial_0\Omega$.

From expression (5), it turns out that the conductivity tensor $K$ is singular. In order to avoid the ill-posedness of the conduction problem, the following isotropic perturbation is introduced [14,5]:

$$ [K]_e = \begin{bmatrix} T_x^2 & T_x T_y \\ T_x T_y & T_y^2 \end{bmatrix} + \epsilon \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \qquad (6) $$
where $T_x$ and $T_y$ are the Cartesian components of the vector field $\mathbf{T}$ and where $\epsilon$ is a user-defined numerical parameter. No rigorous criterion is available for its choice. However, it should be as small as possible in order to fulfil Eq. (4a), but sufficiently large to break down the singularity of $\mathbb{K}$. The dependence of the results on this numerical parameter may then limit the applicability of global tracking, eventually leading to numerical instability issues.

## 2.2. Heat convective-diffusion-like problem

With the aim of overcoming the aforementioned limitation, a new interpretation of the original problem defined by Eq. (2) is here proposed. To start with, let us consider a convection-diffusion-like problem, which consists in finding a temperature field $\theta$ such that:

$$ \underbrace{\mathbf{T} \cdot \nabla \theta}_{\text{Convective term}} - \underbrace{\mathbf{div}(\alpha \nabla \theta)}_{\text{Diffusive term}} = 0 \quad \forall \mathbf{x} \in \Omega \qquad (7a) $$

$$ J = -\alpha \nabla \theta \cdot \mathbf{v} = 0 \quad \forall \mathbf{x} \in \partial_j \Omega \qquad (7b) $$

$$ \theta = \theta^* \quad \forall \mathbf{x} \in \partial_0 \Omega \qquad (7c) $$

where $\alpha$ is a diffusion coefficient and $\mathbf{v}$ the unit outward normal vector to the boundary. Problem (7) can be derived from Navier-Stokes equations in case of incompressible fluids. In this context, the vector field $\mathbf{T}$ takes the physical meaning of a fluid velocity. Therefore Eq. (7a) describes a heat transfer process where two contributions can be distinguished: the first one of convective nature, the second one of diffusive nature. Boundary condition (7b) expresses the Neumann condition on the diffusive flux $J$, while Eq. (7c) assumes the same meaning as in problem (4). However, since $\mathbf{T}$ is now a velocity field, Dirichlet boundary conditions should be prescribed only on the portions of the boundary where the former is directed inward the domain $\Omega$.

In order to better understand the preceding formulation, let us focus upon the one-dimensional case. Assuming a constant diffusion coefficient and considering, on the one hand, a first order upwind scheme for the convective term and, on the other hand, a second order centered scheme for the diffusive term, Eq. (7a) is discretized as:

$$ T \frac{\theta_i - \theta_{i-1}}{\Delta x} - \alpha \frac{\theta_{i+1} - 2\theta_i + \theta_{i-1}}{\Delta x^2} = 0 \qquad (8) $$

where $\Delta x$ is the one-dimensional spatial discretization step. If we consider now a centered discretization of the convective term, it turns out that the difference between the upwind discretization and the aforementioned one is given by:

$$ \begin{aligned} \left[ T \frac{\partial \theta}{\partial x} \right]_{\text{Upwind}} - \left[ T \frac{\partial \theta}{\partial x} \right]_{\text{Centered}} &= T \frac{\theta_i - \theta_{i-1}}{\Delta x} - T \frac{\theta_{i+1} - \theta_{i-1}}{2\Delta x} \\ &= -\frac{T \Delta x}{2} \frac{\theta_{i+1} - 2\theta_i + \theta_{i-1}}{\Delta x^2} \\ &= \left[ \frac{T \Delta x}{2} \frac{\partial^2 \theta}{\partial x^2} \right]_{\text{Centered}} \end{aligned} \qquad (9) $$

Eq. (9) shows that the discretized expression of the convective term by means of the first order upwind scheme is equal to the second order centered scheme discretization of the same term plus an additional diffusive contribution. By comparing Eqs. (8) and (9), we notice that a numerical diffusive term comes out naturally and it is characterized by a diffusion coefficient $T \frac{\Delta x}{2}$ which is function of a mesh characteristic length - $\Delta x$ in the case of one-dimensional problems. This term is similar to the numerical conductivity coefficient $\epsilon$ introduced in Eq. (6) but it is no more user-defined and it tends towards zero as $\Delta x \to 0$, which means that the centered and the upwind discretization schemes are consistent. This observation constitutes the fundamentals of upwind discretization methods, classically used in fluid mechanics [26]. Higher dimensional extensions of the concept of numerical diffusion have been well-established in fluid mechanics in the case of finite element discretizations. More precisely, two formulations are considered in this study. Given the vector field $\mathbf{T}(\mathbf{x})$, the following bi-dimensional extension of the concept of numerical diffusion has been proposed [26] (Streamline Upwind method - SU):

$$ \mathbf{T} \cdot \nabla \theta - \mathbf{div} \left( \frac{h_T \| \mathbf{T} \|}{2} \nabla \theta \right) = 0 \quad \forall \mathbf{x} \in \Omega \qquad (10a) $$

$$ J = - \left( \frac{h_T \| \mathbf{T} \|}{2} \cdot \nabla \theta \right) \cdot \mathbf{v} = 0 \quad \forall \mathbf{x} \in \partial_j \Omega \qquad (10b) $$

$$ \theta = \theta^* \quad \forall \mathbf{x} \in \partial_0 \Omega \qquad (10c) $$

where $h_T$ is a mesh characteristic length. Problem (10) describes the numerical diffusion as an isotropic process and this may lead to imprecise results if coarse meshes are adopted. As a matter of fact, since in this case a large amount of diffusion is introduced in all directions, the isovalues of the temperature-like field do not longer envelop the tangent vector field. Another formulation, which considers an anisotropic numerical diffusion only in the direction of the vector field $\mathbf{T}$ (Streamline Upwind Petrov Galerkin method - SUPG) has also been proposed:
$$ \mathbf{T} \cdot \nabla \theta - \mathbf{div} \left( \frac{h_T \|\mathbf{T}\|}{2} \frac{\mathbf{T} \otimes \mathbf{T}}{\|\mathbf{T}\|^2} \nabla \theta \right) = 0 \quad \forall \mathbf{x} \in \Omega \qquad (11a) $$

$$ J = - \frac{h_T \|\mathbf{T}\|}{2} \frac{\mathbf{T} \otimes \mathbf{T}}{\|\mathbf{T}\|^2} \nabla \theta \cdot \mathbf{v} = 0 \quad \forall \mathbf{x} \in \partial_j \Omega \qquad (11b) $$

$$ \theta = \theta^* \quad \forall \mathbf{x} \in \partial_0 \Omega \qquad (11c) $$

It can be demonstrated that the tensor $\frac{\mathbf{T} \otimes \mathbf{T}}{\|\mathbf{T}\|^2}$ has only one non-zero eigenvalue associated to $\mathbf{T}$ as eigenvector. The product $\frac{\mathbf{T} \otimes \mathbf{T}}{\|\mathbf{T}\|^2} \nabla \theta$ then retains only the part of $\nabla \theta$ parallel to $\mathbf{T}$.

## 2.3. Preliminary results

Concrete structures often experience mixed mode fracture mechanisms, due to the strong interaction between shear and bending. This occurrence can translate into curved crack paths developing from the zones where maximum tensile principal stresses concentrate. Therefore, global tracking must allow to reproduce such behavior. In particular, the capability of the three formulations presented in Sections 2.1 and 2.2 to envelop the principal stress field is now investigated. These are recalled in Table 1.

The well-known experiment performed by Schlangen [27] has been simulated. The specimen is a single-edge notched (SEN) plain concrete beam subjected to four-point shear (see Fig. 2).

The isovalue of the thermal-like field obtained by solving problems (4), (10) and (11) are shown respectively in Figs. 3–5 for three different meshes of 3-node triangular finite elements (879, 3672, 13370 elements). The isovalue crossing the first material point experiencing cracking is also traced in bold. In all simulations, Dirichlet boundary conditions have been imposed on the lower edge of the beam, such that a linearly increasing temperature profile is established onto this boundary.

In case of the HC formulation, a numerical perturbation $\epsilon = 10^{-5}$ has been considered (see [14,28] for typical values of this parameter). As it can be observed in Fig. 3b, a loss of uniqueness of the solution is encountered near the notch for the intermediate discretization.

The isovalue obtained by the HCD-SU formulation, implemented in the finite element software Cast3M-CEA [29], are depicted in Fig. 4. It can be noticed that the solution is now more diffused, which means that the prescribed temperature values are more easily transported all over the domain. This feature gives stability to the solution, but at the same time it provides less accurate results if coarse meshes are adopted. However, for finer discretizations, the results are closed to those obtained for the HC problem.

The HCD-SUPG formulation seems to be a good candidate to reduce the dependence on the mesh refinement. The results are shown in Fig. 5. In fact, the solution reveals to be unstable, even though less diffused than in the previous case. Therefore, the introduction of a numerical diffusion only in the direction of the tangent field $\mathbf{T}$ is not sufficient to solve out our problem. As a consequence of the above considerations, the most appropriate convection-diffusion-like formulation of the global tracking problem is the one given by Eqs. (10).

## 3. Crack-tracking algorithm

In Section 2.2 the problem of tracing the envelopes of a vector field $\mathbf{T}(\mathbf{x}, t)$ has been formulated in terms of Navier-Stokes equations for incompressible fluids. As for the heat-conduction-like formulation, the scalar function $\theta(\mathbf{x}, t)$ represents the temperature field whose isovalue describe all the possible discontinuity lines in the domain $\Omega$.

The choice of the right isovalue stands on the stress distribution at time $t$. In particular, global tracking associates to each discontinuity line $\Gamma_i$ a root $r_i$, i.e. the material point (or the element) at time $t_0$ not belonging to any crack path and satisfying for the first time the activation condition [14,5]. The discontinuity is thus represented as follows:

$$ \Gamma_i(t) = \{ \mathbf{x} \in \Omega \mid \theta(\mathbf{x}, t) = \theta(\mathbf{x}_{r_i}, t) \} \qquad (12a) $$

$$ \mathbf{x}_{r_i} \in \Omega \mid \| \underline{\sigma}(\mathbf{x}_{r_i}, t_0) \| \ge f_t \qquad (12b) $$

where the activation condition is expressed in terms of a certain tensorial norm $\|\cdot\|$ and the material strength $f_t$. The reference isovalue is considered to pass through the centroid of the root element of the discontinuity [14,5].

<table><caption>**Table 1**

Possible strategies to envelop the principal stress field.</caption><thead><tr><th></th><th>Formulation</th></tr></thead><tbody><tr><td>HC</td><td>Heat conduction</td></tr><tr><td>HCD-SU</td><td>Heat convection-diffusion SU</td></tr><tr><td>HCD-SUPG</td><td>Heat convection-diffusion SUPG</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Four point shear of a single-edge notched beam according to Schlangen [27] – dimensions in mm.</figcaption></figure>

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

<figure><img src="image_5.png" /><figcaption>Fig. 4. Isovalues of the HCD-SU problem. Isotropic diffusion ranges: (a) $[1.66 \times 10^{-3}, 6.12 \times 10^{-3}]$, (b) $[9.89 \times 10^{-4}, 3.48 \times 10^{-4}]$, (c) $[3.95 \times 10^{-4}, 2.39 \times 10^{-4}]$.</figcaption></figure>

The previous definition implicitly assumes that the portion of isovalue associated to the *active part* of the crack, i.e. with points characterized by $[\mathbf{u}]_t \neq 0$, does not change any more. In reality, due to the approximative nature of the finite element solution, this is not generally guaranteed: indeed, since principal stress directions are free to rotate where the material is linear elastic, the isovalues of the thermal field may evolve also inside the elements already exhibiting a crack. As a consequence, the nodes of the element domain $\Omega_e$ traversed by the crack may be incorrectly separated, with a loss of continuity of the crack-path and stress-locking effects taking place. A possible solution would be to freeze the nodal temperature of the cracked element by imposing additional Dirichlet boundary conditions [6,28]. However, the initial boundary value problem (4) or (10) would be contradicted with respect to the initial boundary conditions and local techniques should be adopted in order to perform the analysis at each time step [25].
<figure><img src="image_6.png" /><figcaption>Fig. 5. Isovalues of the HCD-SUPG problem. Anisotropic diffusion ranges: (a) [1.66 × 10⁻³, 6.12 × 10⁻³], (b) [9.89 × 10⁻⁴, 3.48 × 10⁻⁴], (c) [3.95 × 10⁻⁴, 2.39 × 10⁻⁴].</figcaption></figure>

This drawback can be found if the root $r_i$ is not updated for $t \ge t_0$. Thus, the tracking procedure should explicitly take into account both the active part $\Gamma_{ia}$ and the potential part $\Gamma_{ip}$ of each discontinuity $\Gamma_i$ and, at the same time, impose the match between the root $r_i$, used to trace the prosecution of the crack, and the crack-tip. In these circumstances, the reference point for tracing the isovalue defining $\Gamma_{ip}$ does not coincide with the root element centroid but instead with a point belonging to one of the root element edges, i.e. the crack-tip. If we consider a time discretization $\{t_0, \dots, t_k, \dots, t_n\}$, the total discontinuity $\Gamma_i$ at time $t_n$ can then be represented as:

$$ \Gamma_i = \Gamma_{ia} \cup \Gamma_{ip} = \bigcup_{k=0}^{n-1} \tilde{\Gamma}_i^{(k)} \cup \tilde{\Gamma}_i^{(n)} \quad (13a) $$

$$ \tilde{\Gamma}_i^{(k)} = \{ \mathbf{x} \in \Omega \mid \theta(\mathbf{x}, t_k) = \theta(\mathbf{x}_{r_i}^{(k)}, t_k), \|\mathbf{u}(\mathbf{x}, t_k)\| \neq 0 \} \quad (13b) $$

$$ \tilde{\Gamma}_i^{(n)} = \{ \mathbf{x} \in \Omega \mid \theta(\mathbf{x}, t_n) = \theta(\mathbf{x}_{r_i}^{(n)}, t_n), \|\mathbf{u}(\mathbf{x}, t_n)\| = 0 \} \quad (13c) $$

$$ \mathbf{x}_{r_i}^{(k)} \in \Omega \mid \| \underline{\sigma}(\mathbf{x}_{r_i}^{(k)}, t_k) \| \ge f_t \quad (13d) $$

From the previous definition, the potential part of the discontinuity $\Gamma_{ip}$ has been defined as the portion of the isovalue triggered off the root $r_i$ at time $t_n$ and characterized by linear elastic behavior. Consequently, only the portion of the isovalue that does not cross the element associated to $r_i$ should be taken into account. Since root $r_i$ divides the potential line into two parts, in order to avoid ambiguity, it seems convenient to orient the curve by setting the origin at the root itself and assume as positive the direction of the propagating discontinuity. This can be done by means of a curvilinear abscissa $s_i$, whose origin is set to coincide with root $r_i$ (see Fig. 6).

Thus, the prosecution of the crack path $\Gamma_i$ at time $t_n$ will be associated to the positive values of $s_i$, with origin at the root $r_i$. This procedure can be translated by the following steps:

1. Make the input point $I_{\Gamma_i}$ match the crack-tip.

2. Trace of the isovalue passing trough $I_{\Gamma_i}$.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Crack path tracking in a bi-dimensional domain by means of the variable-root algorithm.</figcaption></figure>
3. Create the curvilinear abscissa $s_i$ with origin in $I_{r_i}$ and with positive values in the sense of the propagating discontinuity.

4. Choose the potential continuation of discontinuity $\Gamma_i$ as the part of isovalue associated to $s_i > 0$.

The algorithms derived by problem (12) and problem (13) are now compared. The test reported in [30] has been simulated. The geometry is depicted in Fig. 7 and consists in a concrete specimen under tension, characterized by two notches with an offset in the direction of the load. An imposed displacement $\delta$ is applied horizontally on the right side of the structure, while the left side is fixed in this direction. A hinged support is introduced at the upper left corner in order to forbid any rigid body motion.

In Fig. 8a the case of fixed roots $r_1$ and $r_2$ for the crack paths $\Gamma_1$ and $\Gamma_2$ is shown. It appears that the evolution of the reference isovvalues due to the rotation of the principal stress directions inside linear elastic elements may lead to an imprecise evaluation of the propagating discontinuities. In particular, even if the principal stress directions are freezed inside cracked elements, the isovvalues may no longer be able to envelop the tangents to the active part of the discontinuities. A change in the decomposition mode of the elemental domain $\Omega_e$ may then occur whenever the nodes shared by adjacent elements do not belong to the same sub-domain $\Omega_e^+$ or $\Omega_e^-$.

If the root of each discontinuity is updated, it is possible to separate its active part from its potential prosecution. This strategy allows to enforce the continuity of the crack path even if the isovalue distribution evolves during the analysis. As depicted in Fig. 8b, the requirement of a domain of unique decomposition is attained for all the cracked elements.

# 4. Numerical examples

In this section, the performances of the global crack-tracking technique applied to the E-FEM are investigated. A comparison between the version presented in [14,5] and the version proposed in this paper is made by means of bi-dimensional benchmark tests, presented in order of increasing complexity. Both the mathematical formulation and the tracking scheme are here discussed.

## 4.1. Double-edge notched specimen under tension

The first example is the same as the double-edge notched specimen under tension that has been preliminarily studied in Section 3. The formation of two principal crack paths has been observed experimentally in [30]. A correct evaluation of the propagating discontinuities is essential in order to properly simulate the structural response. Firstly, the fixed-root algorithm deriving from Eqs. (12) is applied in the case of the HCD-SU formulation. Three discretizations counting 496, 1691 and 6762 3-node triangular elements respectively have been considered. The corresponding load-displacement curves are plotted in Fig. 9a. The localization zones at the end of each simulation are depicted in Fig. 9b.

It can be noticed that the crack paths start developing correctly from the notches, but then, as the distance from the respective root element increases, a loss of continuity occurs in all the three simulations. As it can be observed in Fig. 9b, the localizations zones are not always defined by simple bands of elements, which means that stress-locking and spurious cracking take place. Consequently the knowledge of the crack-tip position is lost and the constitutive response is not well evaluated. In addition, numerical issues are encountered already at early stages in particular when fine meshes are adopted. These drawbacks have already been reported in [25] and may constitute a limitation to the applicability of the global tracking scheme.

If the variable-root algorithm derived from Eqs. (13) is adopted, the position of the crack-tip is always available during the analysis. This information allows to impose the continuity of each crack path by separating it into an active part and a potential part without the use of any further strategy. As shown in Fig. 10a, the numerical simulations fit pretty well with the experimental result in terms of load-displacement curves, denoting a good mesh-size independence. The localization zones, depicted in Fig. 9b, consist in a fully developed crack-band initiating at the upper notch and a partially developed crack-band starting from the lower notch.

The mathematical formulation is now discussed. A comparison between the heat conduction-like problem and the heat convection-diffusion-like problem is drawn in case of the variable-root global tracking. The load-displacement curves and

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Tension test on double-edge notched specimen according to [30] – dimensions in mm.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Crack propagation scheme.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Results of the fixed-root global tracking algorithm for the Shi's test.</figcaption></figure>

the localizations zones are shown in Figs. 11 and 12 for the intermediate mesh containing 1691 elements. The global struc-
tural response shows minor differences between the two simulations, in particular it seems that the HC approach allows a
better energy dissipation in the final stage of the analysis. By observing the crack trajectories, it appears that an asymmet-
rical propagation takes place for the HCD-SU problem, whereas two almost identical crack paths are found for the HC prob-
lem, which is coherent with the fact that the stress distribution is symmetrical with respect to the vertical axis. This
discrepancy may be due to the different physical meaning of the mathematical formulation. In particular, for the HCD-SU
problem the tangent vector field represents a fluid velocity, deduced from the principal stress directions. Therefore, the solu-
tion is affected by the sense given to the velocity field. Such operation may lead to less accurate results with respect to the HC
formulation if not properly done, especially in presence of singular points or high gradients in the stress distribution.
<figure><img src="image_11.png" /><figcaption>Fig. 10. Results of the variable-root global tracking algorithm for the Shi's test.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Load-displacement curve comparison of the HC and HCD-SU formulations for the Shi's test.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Crack path comparison of the HC and HCD-SU formulations for the Shi's test.</figcaption></figure>
## 4.2. Four-point shear test on a single-edge notched beam

The second example is the numerical simulation of four-point shear of a single-edge notched beam. In literature this problem has been addressed by several authors, see e.g. [31–33,22]. At first, the same SEN beam configuration studied in Section 2.3 is considered. An imposed vertical displacement is applied centrally to the loading platen. A curved crack path is expected to initiate at the right lower corner of the notch.

The fixed-root algorithm is compared to the variable-root algorithm when the HCD-SU is adopted. Fig. 13a and b depict the total applied load F versus the Crack Mouth Opening Displacement (CMOD) for the two cases. In Fig. 14, the same comparison is made in terms of load-Crack Mouth Sliding Displacement (CMSD). It can be noticed that in the case of the SEN beam the global structural response is less sensitive to the evolution of the thermal-like field. This can be explained by considering the fact that the structure undergoes a sudden failure mechanism as soon as the peak-load is reached. As a consequence of this brittle behavior, a fully developed localization zone is observed once equilibrium is re-obtained. Thus, the reference isovalue at the end of the simulation does not differ substantially from the initial one. Nevertheless, in the case of fine meshes, the difference between the two approaches is more important: indeed, spurious cracking arises for the fixed-root algorithm, as it can be seen in Fig. 15a. The band of elements that undergo failure thus does not define a strictly

<figure><img src="image_14.png" /><figcaption>Fig. 13. CMOD curves for the SEN beam.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. CMSD curves for the SEN beam.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 15. Localization zones for the SEN beam.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. CMOD and CMSD comparison of the HC and HCD-SU formulations for the SEN beam.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 17. Crack path comparison of the HC and HCD-SU formulations for the SEN beam.</figcaption></figure>

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

(a) loss of uniqueness of the solution
in the case of the HC formulation

(b) HCD-SU result

<caption>Fig. 18. Problem in the tracking procedure.</caption>
continuous crack path, although it reproduces pretty well the shape of the experimental failure surface. On the contrary, the
continuity requirement is recovered by the variable-root algorithm (see Fig. 15b).

The capability to well reproduce the failure mode is now analyzed once again for the two reference formulations consid-
ered here, i.e. the HC and HCD-SU problems. This is done by means of the variable-root algorithm applied to the intermediate
discretization containing 3672 elements. The global structural response is plotted in Fig. 16 in terms of load-CMOD and load-
CMSD curves, while the resulting cracks paths are depicted in Fig. 17. A good agreement between the simulations can be
observed, although in the case of the heat conduction-like problem a loss of uniqueness of the thermal solution has been
encountered in the proximity of the peak-load. This numerical issue has obliged the authors to stop the tracking procedure
in order to complete the analysis (see Fig. 18).

**5. Conclusions and perspectives**

In this paper a modified global crack-tracking strategy has been applied to the E-FEM simulation of quasi-brittle mate-
rials. Firstly, the problem of tracing the envelopes of a vector field has been addressed by comparing different mathematical
formulations. The Streamline Upwind discretization of the heat convection-diffusion-like problem (HCD-SU) has been cho-
sen as the privileged candidate for this task. Secondly, a revised algorithm for the representation of the propagating discon-
tinuity has been presented. The main ingredient is the updating of the material point used to trace the potential part of the
crack, instead of assuming a fixed-root for its representation. These issues have been discussed in detail by means of repre-
sentative numerical examples and the following conclusions can be drawn:

* Independently from the mathematical formulation of the thermal-like problem, a fixed-root scheme does not prevent the loss of continuity of the crack path due to the rotation of the isovalue. As a consequence, spurious cracking appears and both stress-locking effects and numerical instability issues occur;

* The variable-root scheme is able to provide a $C^0$-continuous crack path at no additional computational cost and without integrating any further technique;

* A good agreement with the experiments is found when coupling the variable-root algorithm to the HCD-SU formulation, although in the case of coarse meshes its precision is reduced with respect to the heat conduction-like approach (HC). However, from the authors' experience, good results have been obtained for all the spatial discretizations adopted in the analysis;

* Attention must be paid to the sense of the tangent vector field, especially in the case of strong gradients of the stress dis-
  tribution. This consideration is important since Dirichlet boundary conditions can be applied only on the portions of the
  boundary where the velocity field is directed inward the domain;

* The HCD-SU formulation provides better stability performances with respect to the HC formulation. The latter is strongly dependent on the numerical conductivity parameter, which may be not sufficient to guarantee stable solutions and therefore to perform the global tracking procedure.

The applicability of the variable-root HCD-SU tracking strategy to three-dimensional problems has not been investigated yet, in particular the variable-root algorithm could be object of further studies. The extension of the proposed approach to branching scenarios should be also evaluated, such as the possibility to handle intersecting cracks.

**Acknowledgements**

The authors wish to express their most grateful thanks to CEA/DEN for its financial support. The work carried out under the SINAPS@ project has benefited from French funding managed by the National Research Agency under the program Future Investments (SINAPS@ reference No. ANR-11-RSNR-0022). This work has also been supported by the SEISM Institute (http://www.institut-seism.fr) and the SIMEX project, managed by Jean-Luc Fayard. The authors thank also Stephane Gounand for its fruitful discussions.

**References**

[1] Simo JC, Oliver J, Armero F. An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. Comput Mech 1993;12(5):277-96.

[2] Simo J, Oliver J. A new approach to the analysis and simulation of strain softening in solids. Fract Damage Quasibrittle Struct 1994;25-39.

[3] Armero F, Garikipati K. An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids. Int J Solids Struct 1996;33(20):2863-85.

[4] Oliver J, Cervera M, Manzoli O. Strong discontinuities and continuum plasticity models: the strong discontinuity approach. Int J Plast 1999;15 (3):319-51.

[5] Oliver J, Huespe A, Samaniego E, Chaves E. Continuum approach to the numerical simulation of material failure in concrete. Int J Numer Anal Methods Geomech 2004;28(7-8):609-32.

[6] Oliver J, Huespe A. Continuum approach to material failure in strong discontinuity settings. Comput Methods Appl Mech Eng 2004;193(30):3195-220.

[7] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Eng 1999;45(5):601-20.

[8] Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Eng 1999;46(1):131-50.

[9] Dolbow J, Moës N, Belytschko T. Discontinuous enrichment in finite elements with a partition of unity method. Finite Elem Anal Des 2000;36 (3):235-60.
[10] Belytschko T, Moës N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. Int J Numer Methods Eng 2001;50(4):993-1013.

[11] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Eng Fract Mech 2002;69(7):813-33.

[12] Samaniego E. Contributions to the continuum modelling of strong discontinuities in two-dimensional solids (PhD thesis); 2003.

[13] Feist C, Hofstetter G. Three-dimensional fracture simulations based on the SDA. Int J Numer Anal Methods Geomech 2007;31(2):189-212.

[14] Oliver J, Huespe A, Samaniego E, Chaves E. On strategies for tracking strong discontinuities in computational failure mechanics. In: Fifth world congress on computational mechanics.

[15] Ventura G, Budyn E, Belytschko T. Vector level sets for description of propagating cracks in finite elements. Int J Numer Methods Eng 2003;58 (10):1571-92.

[16] Mosler J, Meschke G. Embedded crack vs. smeared crack models: a comparison of elementwise discontinuous crack path approaches with emphasis on mesh bias. Comput Methods Appl Mech Eng 2004;193(30):3351-75.

[17] Meschke G, Dumstorff P. Energy-based modeling of cohesive and cohesionless cracks via X-FEM. Comput Methods Appl Mech Eng 2007;196 (21):2338-57.

[18] Dumstorff P, Meschke G. Crack propagation criteria in the framework of X-FEM-based structural analyses. Int J Numer Anal Methods Geomech 2007;31 (2):239.

[19] Jäger P, Steinmann P, Kuhl E. Modeling three-dimensional crack propagation comparison of crack path tracking strategies. Int J Numer Methods Eng 2008;76(9):1328-52.

[20] Alfaiate J, Wells G, Sluys L. On the use of embedded discontinuity elements with crack path continuity for mode-I and mixed-mode fracture. Eng Fract Mech 2002;69(6):661-86.

[21] Cervera M, Pelà L, Clemente R, Roca P. A crack-tracking technique for localized damage in quasi-brittle materials. Eng Fract Mech 2010;77 (13):2431-50.

[22] Zhang Y, Lackner R, Zeiml M, Mang HA. Strong discontinuity embedded approach with standard sos formulation: element formulation, energy-based crack-tracking strategy, and validations. Comput Methods Appl Mech Eng 2015;287:335-66.

[23] Chaves EW. Tracking 3d crack path. In: Proceeding of the international conference on mathematical and statistical modeling in honor of Enrique Castillo ICMSM. Ciudad Real, Spain.

[24] Dufour F, Pijaudier-Cabot G, Choinska M, Huerta A. Extraction of a crack opening from a continuous approach using regularized damage models. Comput Concr 2008;5(4):375-88.

[25] Feist C. A numerical model for cracking of plain concrete based on the strong discontinuity approach (PhD thesis); 2004.

[26] Hughes TJ. Recent progress in the development and understanding of SUPG methods with special reference to the compressible euler and Navier-Stokes equations. Int J Numer Methods Fluids 1987;7(11):1261-75.

[27] Schlangen H. Experimental and numerical analysis of fracture processes in concrete. HERON 1993;38(2).

[28] Cervera M, Chiumenti M. Mesh objective tensile cracking via a local continuum damage model and a crack tracking technique. Comput Methods Appl Mech Eng 2006;196(1):304-20.

[29] Verpeaux P, Millard A, Charras T. Cast3m-cea; 2016. <http://www-cast3m.cea.fr>.

[30] Shi C, Van Dam AG, Van Mier JG, Sluys B. Crack interaction in concrete. Mater Build Struct 2005;6:125-31.

[31] Geers M, De Borst R, Peerlings R. Damage and crack modeling in single-edge and double-edge notched concrete beams. Eng Fract Mech 2000;65 (2):247-61.

[32] Oliver J, Huespe A, Pulido M, Chaves E. From continuum mechanics to fracture mechanics: the strong discontinuity approach. Eng Fract Mech 2002;69 (2):113-36.

[33] Bocca P, Carpinteri A, Valente S. Mixed mode fracture of concrete. Int J Solids Struct 1991;27(9):1139-53.