# Comparaison between thick level-set and cohesive zone models

Andres Parrilla-Gomez, Nicolas Moes, Claude Stolz

► To cite this version:

Andres Parrilla-Gomez, Nicolas Moes, Claude Stolz. Comparaison between thick level-set and cohesive zone models. Advanced Modeling and Simulation in Engineering Sciences, 2015, 25 (1), pp.18-22. 10.1186/s40323-015-0041-9 . hal-01191684

HAL Id: hal-01191684

https://hal.science/hal-01191684

Submitted on 30 May 2024

**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
<figure><img src="image_1.png" /></figure>

## RESEARCH ARTICLE

## Open Access

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

CrossMark

# Comparison between thick level set (TLS) and cohesive zone models

Andrés Parrilla Gómez¹, Nicolas Moës¹* and Claude Stolz¹,²

### Abstract

**Background:** Two main families of methods exist to model failure of quasi-brittle structures. The first one consists on crack based models, like cohesive zone models. The second one is a continuum damage approach that leads to a local loss of stiffness. Local damage models need some regularization in order to avoid spurious localization. A recent one, the thick level set damage model, bridges both families by using level-sets. Cohesive and TLS models are presented. The cohesive one represents quasi-brittle behaviors with good accuracy but requires extra equations to determine the crack path. The TLS has proved its capability to model complex crack paths while easily representing cracks (i.e. displacement jumps); contrary to most damage models.

**Methods:** A one-dimensional analytical relationship is exhibited between TLS and cohesive models. The local damage behavior needed to obtain the same global behavior of a bar than with cohesive model is derived. It depends on the choice of some TLS parameters, notably the characteristic length $l_c$. This local behavior is applied to bi-dimensional simulations of three point bending as well as mixed-mode single edge notched specimens are performed. Results are compared to cohesive simulations, regarding both crack paths and force-CMOD curves.

**Results:** Force-CMOD curves obtained are very similar with both models. Theoretical analysis in 1D and numerical results in 2D indicates that, as $l_c$ goes to zero, TLS results tend to CZM ones.

**Conclusions:** The TLS model yields very similar results to the cohesive one, without the need for extra equations to determine the crack path.

**Keywords:** Damage mechanics, Thick level set, Cohesive zone, X-FEM, Level-set, Three point bending, Single edge notch

### Background

Modeling the failure of quasi-brittle structures is an important aim of numerical simulations. Two main kind of models have been developed [1]: crack based ones, as the cohesive zone model, and damage based ones. The first ones deals with crack evolution in elastic materials whereas damage models consider a continuous transition from sound to totally damaged materials.

Quasi-brittle structures are characterized by the existence of a non-negligible fracture process zone. Linear elasticity fracture mechanics (LEFM) does not apply as it requires this zone to be negligible. In cohesive zone models (CZM), it is assumed that the process zone is concentrated over a line (2D problems) or a surface (3D), that defines the crack

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

<footnote>*Correspondence:
nicolas.moes@ec-nantes.fr
¹ GeM, UMR CNRS 6183,
Ecole Centrale de Nantes, 1
rue de la Noë, 44321 Nantes,
France
Full list of author information
is available at the end of the
article</footnote>
<figure><img src="image_2.png" /></figure>

## RESEARCH ARTICLE

## Open Access

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

CrossMark

# Comparison between thick level set (TLS) and cohesive zone models

Andrés Parrilla Gómez¹, Nicolas Moës¹* and Claude Stolz¹,²

*Correspondence:
nicolas.moes@ec-nantes.fr
¹GeM, UMR CNRS 6183,
Ecole Centrale de Nantes, 1
rue de la Noë, 44321 Nantes,
France
Full list of author information
is available at the end of the
article

### Abstract

**Background:** Two main families of methods exist to model failure of quasi-brittle structures. The first one consists on crack based models, like cohesive zone models. The second one is a continuum damage approach that leads to a local loss of stiffness. Local damage models need some regularization in order to avoid spurious localization. A recent one, the thick level set damage model, bridges both families by using level-sets. Cohesive and TLS models are presented. The cohesive one represents quasi-brittle behaviors with good accuracy but requires extra equations to determine the crack path. The TLS has proved its capability to model complex crack paths while easily representing cracks (i.e. displacement jumps); contrary to most damage models.

**Methods:** A one-dimensional analytical relationship is exhibited between TLS and cohesive models. The local damage behavior needed to obtain the same global behavior of a bar than with cohesive model is derived. It depends on the choice of some TLS parameters, notably the characteristic length $l_c$. This local behavior is applied to bi-dimensional simulations of three point bending as well as mixed-mode single edge notched specimens are performed. Results are compared to cohesive simulations, regarding both crack paths and force-CMOD curves.

**Results:** Force-CMOD curves obtained are very similar with both models. Theoretical analysis in 1D and numerical results in 2D indicates that, as $l_c$ goes to zero, TLS results tend to CZM ones.

**Conclusions:** The TLS model yields very similar results to the cohesive one, without the need for extra equations to determine the crack path.

**Keywords:** Damage mechanics, Thick level set, Cohesive zone, X-FEM, Level-set, Three point bending, Single edge notch

### Background

Modeling the failure of quasi-brittle structures is an important aim of numerical simulations. Two main kind of models have been developed [1]: crack based ones, as the cohesive zone model, and damage based ones. The first ones deals with crack evolution in elastic materials whereas damage models consider a continuous transition from sound to totally damaged materials.

Quasi-brittle structures are characterized by the existence of a non-negligible fracture process zone. Linear elasticity fracture mechanics (LEFM) does not apply as it requires this zone to be negligible. In cohesive zone models (CZM), it is assumed that the process zone is concentrated over a line (2D problems) or a surface (3D), that defines the crack

<figure><img src="image_4.png" /></figure>
path. Damage models make no strong assumption on its size and shape; damaged zone length and width are a priori not negligible compared to the size of the material domain.

The CZM definition needs to know the crack path. If it is unknown, LEFM based methods exist to predict it. Even if most cohesive models cannot handle complex crack paths, such as branching, some adaptations [2, 3] allow to do so. Regarding damage models, there exist different families. Indeed, purely local damage presents spurious localization and different regularization methods have been proposed: strain or damage averaging, higher-order gradient, phase-field, variational models... These models are able to deal with complex damaged zone shape evolution but present important difficulties to deal with jump in displacement. Moreover, as damage evolution law has to be enforced over the whole material domain, computational cost is more important than for cohesive models.

The crack band model [4] is an intermediate between the cohesive model and classical damage ones, as damage is represented as a non-linear behavior of some elements that defines the crack path. However, damage band width is only dependent on the chosen elements that have a non-linear behavior and the crack path has to be previously determined. More recently, a level-set based damage model, the thick level set (TLS) model, has been presented [5]. This damage model bridges both families of models as it is able to recover complex crack paths while introducing jumps in displacement in a natural way. Contrary to the crack band model, the width of the damage band is here a parameter that is not related to the mesh.

The goal of this paper is to establish a one-dimensional equivalence between the cohesive zone model and a damage model with TLS regularization. Note that similar analysis concerning cohesive and gradient damage models has been performed by Lorentz in [6]. Thus, based on this equivalence, a method to derive a local damage behavior to use in TLS from any cohesive behavior will be exhibited. Later, it will be used in two-dimensional simulations and results will be compared to classical cohesive simulations, that have proved their capability to reproduce experimental results.

We consider structures under quasi-static load. We restrict the study to small strain and displacement and no particular shear behavior is considered, unlike in [7]. The relationship between cohesive and TLS models is analyzed. A one-dimensional equivalence is analytically first established. Then, these results are used to compare two-dimensional behaviors.

Let us start with the one-dimensional comparison between the cohesive and the TLS models. We consider an elastic bar of length $2L: x \in [-L, L]$ and Young modulus $E$. Degradation in this bar will be modeled either with a cohesive zone located at $x = 0$ or with an evolving damage layer centered at $x = 0$. Due to symmetry, only half of the bar needs to be considered: $x \in [0, L]$. Let $u(x)$ denote the displacement along the bar and in particular $u(L) = u_L$ the displacement of the extremity. The stress $\sigma$ is here a scalar variable.

In the cohesive zone model, the displacement at $x = 0$ is discontinuous and the cohesive zone opening is given by $w = 2u(0^+)$. Regarding the TLS model, the extent of the damage layer over the half-bar is denoted by $[0, l]$ where $x = l$ corresponds to the current damage front location. An example of the displacement $u(x)$ for both cohesive and TLS models is shown in Fig. 1.
path. Damage models make no strong assumption on its size and shape; damaged zone length and width are a priori not negligible compared to the size of the material domain.

The CZM definition needs to know the crack path. If it is unknown, LEFM based methods exist to predict it. Even if most cohesive models cannot handle complex crack paths, such as branching, some adaptations [2, 3] allow to do so. Regarding damage models, there exist different families. Indeed, purely local damage presents spurious localization and different regularization methods have been proposed: strain or damage averaging, higher-order gradient, phase-field, variational models... These models are able to deal with complex damaged zone shape evolution but present important difficulties to deal with jump in displacement. Moreover, as damage evolution law has to be enforced over the whole material domain, computational cost is more important than for cohesive models.

The crack band model [4] is an intermediate between the cohesive model and classical damage ones, as damage is represented as a non-linear behavior of some elements that defines the crack path. However, damage band width is only dependent on the chosen elements that have a non-linear behavior and the crack path has to be previously determined. More recently, a level-set based damage model, the thick level set (TLS) model, has been presented [5]. This damage model bridges both families of models as it is able to recover complex crack paths while introducing jumps in displacement in a natural way. Contrary to the crack band model, the width of the damage band is here a parameter that is not related to the mesh.

The goal of this paper is to establish a one-dimensional equivalence between the cohesive zone model and a damage model with TLS regularization. Note that similar analysis concerning cohesive and gradient damage models has been performed by Lorentz in [6]. Thus, based on this equivalence, a method to derive a local damage behavior to use in TLS from any cohesive behavior will be exhibited. Later, it will be used in two-dimensional simulations and results will be compared to classical cohesive simulations, that have proved their capability to reproduce experimental results.

We consider structures under quasi-static load. We restrict the study to small strain and displacement and no particular shear behavior is considered, unlike in [7]. The relationship between cohesive and TLS models is analyzed. A one-dimensional equivalence is analytically first established. Then, these results are used to compare two-dimensional behaviors.

Let us start with the one-dimensional comparison between the cohesive and the TLS models. We consider an elastic bar of length $2L: x \in [-L, L]$ and Young modulus $E$. Degradation in this bar will be modeled either with a cohesive zone located at $x = 0$ or with an evolving damage layer centered at $x = 0$. Due to symmetry, only half of the bar needs to be considered: $x \in [0, L]$. Let $u(x)$ denote the displacement along the bar and in particular $u(L) = u_L$ the displacement of the extremity. The stress $\sigma$ is here a scalar variable.

In the cohesive zone model, the displacement at $x = 0$ is discontinuous and the cohesive zone opening is given by $w = 2u(0^+)$. Regarding the TLS model, the extent of the damage layer over the half-bar is denoted by $[0, l]$ where $x = l$ corresponds to the current damage front location. An example of the displacement $u(x)$ for both cohesive and TLS models is shown in Fig. 1.
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Displacement over the whole bar for CZM and TLS models for a partially damaged bar. L is the length of the bar and l the position of the damage front in TLS model.</figcaption></figure>

## Cohesive zone model

The cohesive zone model was first introduced by Dugdale in [8] for ductile materials and by Barenblatt in [9] for concrete. The model represents the progressive fracture process by condensing over a crack the effect of the whole fracture process zone. A cohesive behavior is imposed between the crack lips (see Fig. 2a). The bi-linear cohesive law, well-known for concrete behavior, is presented in Fig. 3a. This model is widely used [1] and is considered to have good capabilities to fit experimental results when the crack path is known. If it is not previously known, a complementary method, like the maximum tangential stress (MTS) criterion [10], has to be applied.

Let us consider the following free energy of the cohesive zone, located at $x = 0$.

$$ \Phi_{CZM}(w, \alpha) = \frac{1}{2} g_{CZM}(\alpha) kw^2 \quad (1) $$

where $k > 0$ and $g_{CZM}$ is a decreasing dimensionless function that characterizes the stiffness of the cohesive zone. We derive the dual quantities $t$ (tension) and $A$

$$ \sigma = \frac{\partial \Phi_{CZM}}{\partial w} = g_{CZM}(\alpha) kw \quad (2) $$

<figure><img src="image_5.png" /><figcaption>**Fig. 2** Representation of an open crack with a cohesive zone (left) and with a TLS damage zone (right). Dark (light) gray indicates (un) damaged.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 1 Displacement over the whole bar for CZM and TLS models for a partially damaged bar. L is the length of the bar and l the position of the damage front in TLS model.</figcaption></figure>

## Cohesive zone model

The cohesive zone model was first introduced by Dugdale in [8] for ductile materials and by Barenblatt in [9] for concrete. The model represents the progressive fracture process by condensing over a crack the effect of the whole fracture process zone. A cohesive behavior is imposed between the crack lips (see Fig. 2a). The bi-linear cohesive law, well-known for concrete behavior, is presented in Fig. 3a. This model is widely used [1] and is considered to have good capabilities to fit experimental results when the crack path is known. If it is not previously known, a complementary method, like the maximum tangential stress (MTS) criterion [10], has to be applied.

Let us consider the following free energy of the cohesive zone, located at $x = 0$.

$$ \Phi_{CZM}(w, \alpha) = \frac{1}{2} g_{CZM}(\alpha) kw^2 \quad (1) $$

where $k > 0$ and $g_{CZM}$ is a decreasing dimensionless function that characterizes the stiffness of the cohesive zone. We derive the dual quantities $t$ (tension) and $A$

$$ \sigma = \frac{\partial \Phi_{CZM}}{\partial w} = g_{CZM}(\alpha) kw \quad (2) $$

<figure><img src="image_6.png" /><figcaption>Fig. 2 Representation of an open crack with a cohesive zone (left) and with a TLS damage zone (right). Dark (light) gray indicates (un) damaged.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 3 Cohesive bi-linear law and TLS local behavior for different $l_c$.</figcaption></figure>

$$A = - \frac{\partial \Phi_{CZM}}{\partial \alpha} = - \frac{1}{2} g'_{CZM}(\alpha) kw^2 \quad (3)$$

The evolution of $\alpha$ is governed by the criterion $f_{CZM}(A, \alpha) = A - A_c(\alpha)$ with the associated law:

$$f_{CZM}(A, \alpha) \le 0, \quad \dot{\alpha} \ge 0 \quad \text{and} \quad f_{CZM}(A, \alpha) \dot{\alpha} = 0 \quad (4)$$

where $A_c(\alpha)$ is supposed to be finite for all $\alpha$.

The model depends on two functional parameters: $g_{CZM}(\alpha)$ and $A_c(\alpha)$. When they are known, and assuming the uniqueness of evolution of $\alpha$ for any opening history, we can derive a tension-opening curve $f$ that describes the tension of the cohesive zone under a monotonously increasing opening $w$. We will denote by $F$ the corresponding dimensionless function:

$$\sigma = f(w) \quad (5)$$

The relation between $g_{CZM}$, $A_c$ and $f$ writes from (2), (3), (4) and (5):

$$\forall \alpha, \quad kg_{CZM}(\alpha) \sqrt{\frac{-2A_c(\alpha)}{kg'_{CZM}(\alpha)}} = f\left(\sqrt{\frac{-2A_c(\alpha)}{kg'_{CZM}(\alpha)}}\right) \quad (6)$$

Here, we will only consider cohesive behaviors with infinite initial stiffness and null stiffness for $\alpha = 1$ (see Fig. 4), that is

$$\lim_{\alpha \to 0^+} g_{CZM}(\alpha) = +\infty \quad \text{and} \quad g_{CZM}(1) = 0 \quad (7)$$

We define two characteristic values: $\sigma_f$ and $w_f$

$$w_f = \lim_{\alpha \to 1^-} w(\alpha) = \lim_{\alpha \to 1^-} \sqrt{\frac{-2A_c(\alpha)}{kg'_{CZM}(\alpha)}} \quad (8)$$
<figure><img src="image_7.png" /><figcaption>Fig. 3 Cohesive bi-linear law and TLS local behavior for different $l_c$.</figcaption></figure>

$$A = - \frac{\partial \Phi_{CZM}}{\partial \alpha} = - \frac{1}{2} g'_{CZM}(\alpha) kw^2 \quad (3)$$

The evolution of $\alpha$ is governed by the criterion $f_{CZM}(A, \alpha) = A - A_c(\alpha)$ with the associated law:

$$f_{CZM}(A, \alpha) \le 0, \quad \dot{\alpha} \ge 0 \quad \text{and} \quad f_{CZM}(A, \alpha) \dot{\alpha} = 0 \quad (4)$$

where $A_c(\alpha)$ is supposed to be finite for all $\alpha$.

The model depends on two functional parameters: $g_{CZM}(\alpha)$ and $A_c(\alpha)$. When they are known, and assuming the uniqueness of evolution of $\alpha$ for any opening history, we can derive a tension-opening curve $f$ that describes the tension of the cohesive zone under a monotonously increasing opening $w$. We will denote by $F$ the corresponding dimensionless function:

$$\sigma = f(w) \quad (5)$$

The relation between $g_{CZM}$, $A_c$ and $f$ writes from (2), (3), (4) and (5):

$$\forall \alpha, \quad kg_{CZM}(\alpha) \sqrt{\frac{-2A_c(\alpha)}{kg'_{CZM}(\alpha)}} = f\left(\sqrt{\frac{-2A_c(\alpha)}{kg'_{CZM}(\alpha)}}\right) \quad (6)$$

Here, we will only consider cohesive behaviors with infinite initial stiffness and null stiffness for $\alpha = 1$ (see Fig. 4), that is

$$\lim_{\alpha \to 0^+} g_{CZM}(\alpha) = +\infty \quad \text{and} \quad g_{CZM}(1) = 0 \quad (7)$$

We define two characteristic values: $\sigma_f$ and $w_f$

$$w_f = \lim_{\alpha \to 1^-} w(\alpha) = \lim_{\alpha \to 1^-} \sqrt{\frac{-2A_c(\alpha)}{kg'_{CZM}(\alpha)}} \quad (8)$$
<figure><img src="image_7.png" /><figcaption>Fig. 4 Response of a typical cohesive zone model for an imposed increasing crack opening.</figcaption></figure>

$$ \sigma_f = \lim_{\alpha \to 0^+} \sigma(\alpha) = \lim_{\alpha \to 0^+} \left( g_{CZM}(\alpha) \sqrt{\frac{-2kA_c(\alpha)}{g'_{CZM}(\alpha)}} \right) \quad (9) $$

which are supposed to be finite. Thus, as $g_{CZM}(\alpha)$ is monotonous, (7) implies

$$ \lim_{\alpha \to 0^+} g'_{CZM}(\alpha) = -\infty \quad \text{and} \quad \lim_{\alpha \to 0^+} w(\alpha) = 0 \quad (10) $$

We denote by $F$ the dimensionless cohesive function depicted on Fig. 4 and associated to $f$ by

$$ \sigma = f(w) = \sigma_f F\left(\frac{w}{w_f}\right) \quad (11) $$

## Thick level set damage model

The bar is made of a material with the following free energy density

$$ \Phi_{\text{dam}}(\varepsilon, d) = \frac{1}{2}g_{\text{dam}}(d)E\varepsilon^2 \quad (12) $$

where $g_{\text{dam}}$ is a decreasing dimensionless function respecting $g_{\text{dam}}(0) = 1$ that characterizes the remaining stiffness for a given amount of damage $d$. We can derive the stress $\sigma$ and the local energy release rate $Y$

$$ Y = \frac{\partial \Phi_{\text{dam}}}{\partial \varepsilon} = -\frac{1}{2}g_{\text{dam}}(d)E\varepsilon \quad (13) $$

$$ Y = -\frac{\partial \Phi_{\text{dam}}}{\partial d} = -\frac{1}{2}g'_{\text{dam}}(d)E\varepsilon^2 \quad (14) $$

We can write a local damage evolution criterion $f_{\text{dam}}(Y, d) = Y - h(d)Y_c$ and an associated law:

$$ f_{\text{dam}}(Y, d) \le 0, \quad \dot{d} \ge 0 \quad \text{and} \quad f_{\text{dam}}(Y, d) \dot{d} = 0 \quad (15) $$
<figure><img src="image_8.png" /><figcaption>Fig. 4 Response of a typical cohesive zone model for an imposed increasing crack opening.</figcaption></figure>

$$ \sigma_f = \lim_{\alpha \to 0^+} \sigma(\alpha) = \lim_{\alpha \to 0^+} \left( g_{CZM}(\alpha) \sqrt{\frac{-2kA_c(\alpha)}{g'_{CZM}(\alpha)}} \right) \quad (9) $$

which are supposed to be finite. Thus, as $g_{CZM}(\alpha)$ is monotonous, (7) implies

$$ \lim_{\alpha \to 0^+} g'_{CZM}(\alpha) = -\infty \quad \text{and} \quad \lim_{\alpha \to 0^+} w(\alpha) = 0 \quad (10) $$

We denote by $F$ the dimensionless cohesive function depicted on Fig. 4 and associated to $f$ by

$$ \sigma = f(w) = \sigma_f F\left(\frac{w}{w_f}\right) \quad (11) $$

## Thick level set damage model

The bar is made of a material with the following free energy density

$$ \Phi_{\text{dam}}(\varepsilon, d) = \frac{1}{2}g_{\text{dam}}(d)E\varepsilon^2 \quad (12) $$

where $g_{\text{dam}}$ is a decreasing dimensionless function respecting $g_{\text{dam}}(0) = 1$ that characterizes the remaining stiffness for a given amount of damage $d$. We can derive the stress $\sigma$ and the local energy release rate $Y$

$$ Y = \frac{\partial \Phi_{\text{dam}}}{\partial \varepsilon} = -\frac{1}{2}g_{\text{dam}}(d)E\varepsilon \quad (13) $$

$$ Y = -\frac{\partial \Phi_{\text{dam}}}{\partial d} = -\frac{1}{2}g'_{\text{dam}}(d)E\varepsilon^2 \quad (14) $$

We can write a local damage evolution criterion $f_{\text{dam}}(Y, d) = Y - h(d)Y_c$ and an associated law:

$$ f_{\text{dam}}(Y, d) \le 0, \quad \dot{d} \ge 0 \quad \text{and} \quad f_{\text{dam}}(Y, d) \dot{d} = 0 \quad (15) $$
$Y_c$ is a constant. The dependence on $d$ of the critical energy release rate is concentrated on the dimensionless function $h(d)$. The interest of making the critical energy release rate depend on $d$ has been shown in [11, 12].

This local damage model is dependent on two functional parameters: $g_{dam}(d)$ and $h(d)$. When they are fixed, we can derive a strain-stress curve $\sigma = f(\varepsilon)$ that describes the local behavior under a monotonously increasing strain $\varepsilon$. An example of strain-stress curve is shown in Fig. 5a. The function $h(d)$ is chosen to be 1 for $d = 0$. So, damage begins when $Y$ is locally equal to $Y_c$. In order to ensure uniqueness of the history of damage for a given strain history, we must ensure that $\dot{d} > 0$ implies $\dot{\varepsilon} > 0$. This condition writes

$$g''_{dam}h - g'_{dam}h' \geq 0 \qquad (16)$$

This corresponds to curves without snap-back. An example of strain-stress curve for which this condition is not ensured, is illustrated in Fig. 5b.

Local damage models are well-known for their spurious localization due to softening. Therefore some regularization is required. The thick level set damage model [5, 13–15] is one of them. The main idea is to consider a damage front that drives damage evolution over the non-local damage zone. A level set function $\phi$ is defined over the domain as the signed distance to a curve, called the damage front. If $\phi \leq 0$ the material is undamaged and if $\phi > 0$, it is damaged (see Fig. 2b). The value of damage $d$ is chosen to depend only on the level set value:

$$d(\phi/\ell_c \leq 0) = 0, \quad d'(0 < \phi/\ell_c < 1) \geq 0 \quad \text{and} \quad d(\phi/\ell_c \geq 1) = 1 \qquad (17)$$

where $\ell_c$ is a given material length. The condition of irreversibility $\dot{d} \geq 0$ is thus equivalent to $\dot{\phi} \geq 0$. An example of damage law is given in Fig. 6.

As damage is not a free variable any more and $\phi$ is a signed distance function, damage evolution cannot occur independently at each point of the damage band. In the one-dimensional case of the bar, $\dot{\phi}$ is uniform over the whole domain and its value is $\dot{\varepsilon}$.

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

a. Strain-stress curve from an acceptable and continuous $h$ function.

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

b. Strain-stress curve from a non acceptable $h$ function. Red dashed line corresponds to the forbidden part of the curve.

<caption>Fig. 5 Strain-stress curves.</caption>
$Y_c$ is a constant. The dependence on $d$ of the critical energy release rate is concentrated on the dimensionless function $h(d)$. The interest of making the critical energy release rate depend on $d$ has been shown in [11, 12].

This local damage model is dependent on two functional parameters: $g_{dam}(d)$ and $h(d)$. When they are fixed, we can derive a strain-stress curve $\sigma = f(\varepsilon)$ that describes the local behavior under a monotonously increasing strain $\varepsilon$. An example of strain-stress curve is shown in Fig. 5a. The function $h(d)$ is chosen to be 1 for $d = 0$. So, damage begins when $Y$ is locally equal to $Y_c$. In order to ensure uniqueness of the history of damage for a given strain history, we must ensure that $\dot{d} > 0$ implies $\dot{\varepsilon} > 0$. This condition writes

$$g''_{dam}h - g'_{dam}h' \geq 0 \qquad (16)$$

This corresponds to curves without snap-back. An example of strain-stress curve for which this condition is not ensured, is illustrated in Fig. 5b.

Local damage models are well-known for their spurious localization due to softening. Therefore some regularization is required. The thick level set damage model [5, 13–15] is one of them. The main idea is to consider a damage front that drives damage evolution over the non-local damage zone. A level set function $\phi$ is defined over the domain as the signed distance to a curve, called the damage front. If $\phi \leq 0$ the material is undamaged and if $\phi > 0$, it is damaged (see Fig. 2b). The value of damage $d$ is chosen to depend only on the level set value:

$$d(\phi/\ell_c \leq 0) = 0, \quad d'(0 < \phi/\ell_c < 1) \geq 0 \quad \text{and} \quad d(\phi/\ell_c \geq 1) = 1 \qquad (17)$$

where $\ell_c$ is a given material length. The condition of irreversibility $\dot{d} \geq 0$ is thus equivalent to $\dot{\phi} \geq 0$. An example of damage law is given in Fig. 6.

As damage is not a free variable any more and $\phi$ is a signed distance function, damage evolution cannot occur independently at each point of the damage band. In the one-dimensional case of the bar, $\dot{\phi}$ is uniform over the whole domain and its value is $\dot{\varepsilon}$.

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

a. Strain-stress curve from an acceptable and continuous $h$ function.

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

b. Strain-stress curve from a non acceptable $h$ function. Red dashed line corresponds to the forbidden part of the curve.

<caption>Fig. 5 Strain-stress curves.</caption>
<figure><img src="image_10.png" /><figcaption>Fig. 6 TLS classical damage profile.</figcaption></figure>

Indeed, damage is simply driven by the position $l$ of the damage front. In 2D or 3D problems, $\dot{\phi}$ is only uniform over lines perpendicular to the damage front. The local damage evolution criterion is enforced in an averaged way:

$$f_{\text{TLS}} = \overline{f_{\text{dam}}} = \frac{\int_{\phi=0}^{\ell} f_{\text{dam}}(Y, d) d'(\phi/\ell_c) \, d\phi}{\int_{\phi=0}^{\ell} d'(\phi/\ell_c) \, d\phi}, \quad f_{\text{TLS}} \le 0, \quad \dot{l} \ge 0 \quad \text{and} \quad f_{\text{TLS}} \dot{l} = 0 \qquad (18)$$

## Methods

In this section, we establish relations that must be satisfied in order to have equivalent models for the particular problem of the bar. Then, the relations will be specialized for a particular choice of $g_{\text{dam}}(d)$, $d(\phi/\ell_c)$ and finally two examples of cohesive tension-opening functions $f$.

### One-dimensional equivalence

Imagine that, in a quasi-static analysis, some displacement $u_L$ is imposed. We say that a CZM and a TLS models are equivalent if, for a given displacement $u_L$, the same stress is applied and they have the same energy and dissipation [16–19].

The displacement, energy and dissipation of the bar are computed in the following sections, first for cohesive model and secondly for TLS one. The stress can be directly calculated from the constitutive models (2) and (13)

### Cohesive zone model

Equilibrium enforces uniformity of stress which value is $\sigma$. The displacement at the end of the bar writes

$$u_L^{\text{CZM}} = \int_{x=0^+}^{L} \varepsilon_{\text{CZM}} dx + \frac{w}{2} = \frac{L\sigma}{E} + \frac{\sigma}{2kg_{\text{CZM}}(\alpha)} \qquad (19)$$

where the second equality stems from the constitutive model. The energy of the half-bar has two terms: the contribution of the elastic energy of the bar and the energy of the cohesive zone.
<figure><img src="image_11.png" /><figcaption>Fig. 6 TLS classical damage profile.</figcaption></figure>

Indeed, damage is simply driven by the position *l* of the damage front. In 2D or 3D problems, *φ* is only uniform over lines perpendicular to the damage front. The local damage evolution criterion is enforced in an averaged way:

$$f_{\text{TLS}} = \overline{f_{\text{dam}}} = \frac{\int_{\phi=0}^{\ell} f_{\text{dam}}(Y, d) d'(\phi/\ell_c) \, d\phi}{\int_{\phi=0}^{\ell} d'(\phi/\ell_c) \, d\phi}, \quad f_{\text{TLS}} \le 0, \quad \dot{\ell} \ge 0 \quad \text{and} \quad f_{\text{TLS}} \dot{\ell} = 0 \qquad (18)$$

## Methods

In this section, we establish relations that must be satisfied in order to have equivalent models for the particular problem of the bar. Then, the relations will be specialized for a particular choice of $g_{\text{dam}}(d)$, $d(\phi/\ell_c)$ and finally two examples of cohesive tension-opening functions $f$.

### One-dimensional equivalence

Imagine that, in a quasi-static analysis, some displacement $u_L$ is imposed. We say that a CZM and a TLS models are equivalent if, for a given displacement $u_L$, the same stress is applied and they have the same energy and dissipation [16–19].

The displacement, energy and dissipation of the bar are computed in the following sections, first for cohesive model and secondly for TLS one. The stress can be directly calculated from the constitutive models (2) and (13).

### Cohesive zone model

Equilibrium enforces uniformity of stress which value is $\sigma$. The displacement at the end of the bar writes

$$u_L^{\text{CZM}} = \int_{x=0^+}^{L} \varepsilon_{\text{CZM}} dx + \frac{w}{2} = \frac{L\sigma}{E} + \frac{\sigma}{2kg_{\text{CZM}}(\alpha)} \qquad (19)$$

where the second equality stems from the constitutive model. The energy of the half-bar has two terms: the contribution of the elastic energy of the bar and the energy of the cohesive zone.
$$ \mathcal{E}_{\text{CZM}}(u, w, \alpha) = \int_{x=0^+}^{L} \frac{1}{2} E \varepsilon(u)^2 \, dx + \frac{1}{2} \Phi_{\text{CZM}}(w, \alpha) = \frac{\sigma^2 L}{2E} + \frac{\sigma^2}{4g_{\text{CZM}}(\alpha)k} \quad (20) $$

The dissipation writes

$$ \mathcal{D}_{\text{CZM}} = \frac{1}{2} A \dot{\alpha} = \frac{1}{2} A_c(\alpha) \dot{\alpha} \quad (21) $$

### Thick level set model

The same computations that have been done for the cohesive model are here performed for the TLS. Equilibrium leads to $\sigma(x) = \sigma$. The displacement at the end of the bar writes

$$ u_{L}^{\text{TLS}} = \int_{x=0}^{L} \varepsilon_{\text{TLS}} dx = \int_{0}^{L} \frac{\sigma}{E g_{\text{dam}}(d)} dx \quad (22) $$

The energy of the half-bar is

$$ \mathcal{E}_{\text{TLS}}(u, d) = \int_{x=0}^{L} \Phi_{\text{dam}}(u, d) \, dx = \int_{x=0}^{L} \frac{\sigma^2}{2g_{\text{dam}}(d)E} \, dx \quad (23) $$

The dissipation writes

$$ \mathcal{D}_{\text{TLS}}(Y, \dot{d}) = \int_{x=0}^{L} Y \dot{d} \, dx = Y_c \dot{l} \int_{\phi=0}^{\ell} h(d) \frac{d'(\phi/\ell_c)}{\ell_c} \, d\phi \quad (24) $$

Let us remark that in the one-dimensional case of the bar, the TLS variable $\phi$ is related to the coordinate of a point of the bar $x$ and the position of the damage front $l$ by

$$ \phi = \ell - x \quad (25) $$

Moreover, we remind for a better comprehension of the following developments that at the center of the bar ($x = 0$), the value of damage is $d(\ell/\ell_c)$.

### General relationships

We have determined the displacement $u_L$, the total energy of the half-bar $\mathcal{E}$ and the dissipation $\mathcal{D}$ linked to the evolution of the internal variable ($\alpha$ or $d$) for both cohesive and TLS models. The equality of each one of these three quantities for both models leads to (26) and (27). In order to close the system, we need to remind the normality rule for one of the models (28). The second one is redundant.

$$ \frac{1}{2kg_{\text{CZM}}(\alpha)} = \frac{1}{E} \int_{\phi=0}^{\ell} \left( \frac{1}{g_{\text{dam}}(d(\phi/\ell_c))} - 1 \right) d\phi \quad (26) $$

$$ \frac{1}{2} A_c(\alpha) \dot{\alpha} = Y_c \dot{l} \int_{\tilde{d}=0}^{d(\ell/\ell_c)} h(\tilde{d}) \, d\tilde{d} \quad (27) $$

$$ A_c(\alpha) \dot{\alpha} = - \frac{\sigma^2 g'_{\text{CZM}}(\alpha)}{2k g^2_{\text{CZM}}(\alpha)} \dot{\alpha} \quad (28) $$
$$ \mathcal{E}_{\text{CZM}}(u, w, \alpha) = \int_{x=0^+}^{L} \frac{1}{2} E \varepsilon(u)^2 \, dx + \frac{1}{2} \Phi_{\text{CZM}}(w, \alpha) = \frac{\sigma^2 L}{2E} + \frac{\sigma^2}{4g_{\text{CZM}}(\alpha)k} \quad (20) $$

The dissipation writes

$$ \mathcal{D}_{\text{CZM}} = \frac{1}{2} A \dot{\alpha} = \frac{1}{2} A_c(\alpha) \dot{\alpha} \quad (21) $$

### Thick level set model

The same computations that have been done for the cohesive model are here performed for the TLS. Equilibrium leads to $\sigma(x) = \sigma$. The displacement at the end of the bar writes

$$ u_{L}^{\text{TLS}} = \int_{x=0}^{L} \varepsilon_{\text{TLS}} dx = \int_{0}^{L} \frac{\sigma}{E g_{\text{dam}}(d)} dx \quad (22) $$

The energy of the half-bar is

$$ \mathcal{E}_{\text{TLS}}(u, d) = \int_{x=0}^{L} \Phi_{\text{dam}}(u, d) \, dx = \int_{x=0}^{L} \frac{\sigma^2}{2g_{\text{dam}}(d)E} \, dx \quad (23) $$

The dissipation writes

$$ \mathcal{D}_{\text{TLS}}(Y, \dot{d}) = \int_{x=0}^{L} Y \dot{d} \, dx = Y_c \dot{l} \int_{\phi=0}^{\ell} h(d) \frac{d'(\phi/\ell_c)}{\ell_c} \, d\phi \quad (24) $$

Let us remark that in the one-dimensional case of the bar, the TLS variable $\phi$ is related to the coordinate of a point of the bar $x$ and the position of the damage front $l$ by

$$ \phi = \ell - x \quad (25) $$

Moreover, we remind for a better comprehension of the following developments that at the center of the bar ($x = 0$), the value of damage is $d(\ell/\ell_c)$.

### General relationships

We have determined the displacement $u_L$, the total energy of the half-bar $\mathcal{E}$ and the dissipation $\mathcal{D}$ linked to the evolution of the internal variable ($\alpha$ or $d$) for both cohesive and TLS models. The equality of each one of these three quantities for both models leads to (26) and (27). In order to close the system, we need to remind the normality rule for one of the models (28). The second one is redundant.

$$ \frac{1}{2kg_{\text{CZM}}(\alpha)} = \frac{1}{E} \int_{\phi=0}^{\ell} \left( \frac{1}{g_{\text{dam}}(d(\phi/\ell_c))} - 1 \right) d\phi \quad (26) $$

$$ \frac{1}{2} A_c(\alpha) \dot{\alpha} = Y_c \dot{l} \int_{\tilde{d}=0}^{d(\ell/\ell_c)} h(\tilde{d}) \, d\tilde{d} \quad (27) $$

$$ A_c(\alpha) \dot{\alpha} = - \frac{\sigma^2 g'_{\text{CZM}}(\alpha)}{2k g^2_{\text{CZM}}(\alpha)} \dot{\alpha} \quad (28) $$
These equations summarize the relations between the two models. They link the cohesive parameters $\alpha$, $g_{CZM}$ and $A_c$ to the damage ones $\ell$, $g_{dam}$ and $h$ and the TLS regularization function $d(\phi/\ell_c)$. As these equations hold during the whole damage process, that is for all possible values of $\alpha$ or $\ell$, it is possible to derive a TLS damage model from a cohesive one (and vice versa) by fixing four of the parameters. For example, if $d(\phi/\ell_c)$, $g_{dam}(d)$, $g_{CZM}(\alpha)$ and $A_c(\alpha)$ are known, the function $h(d)$ can be deduced (Table 1 presents a summary of the main parameters).

In what follows, we look for $h$ by only using the stress-opening function $f$ and the opening $w$, instead of fully describing $g_{CZM}$ and $A_c$ and $\alpha$ [see (6)]. By differentiating Eq. (26) with respect to time and by combining Eqs. (27) and (28), one obtains

$$ -\frac{1}{2k} \frac{g'_{CZM}(\alpha)}{g_{CZM}^2(\alpha)} \dot{\alpha} = \frac{1}{E} \dot{l} \left( \frac{1}{g_{dam}(d(\ell/\ell_c))} - 1 \right) \quad (29) $$

$$ 2Y_c \dot{l} \int_{\tilde{d}=0}^{d(\ell/\ell_c)} h(\tilde{d}) \, d\tilde{d} = - \frac{\sigma^2}{2k} \frac{g'_{CZM}(\alpha)}{g_{CZM}^2(\alpha)} \dot{\alpha} \quad (30) $$

Moreover, as $\sigma = kg_{CZM}w$, the system of equations is rewritten as

$$ \sigma = f(w) \quad (31) $$

$$ w = \frac{2\sigma}{E} \int_{0}^{\ell} \left( \frac{1}{g_{\text{dam}}(d(\phi))} - 1 \right) d\phi \quad (32) $$

$$ \sigma = \sqrt{2EY_c \frac{g_{\text{dam}}(d(\ell/\ell_c))}{1 - g_{\text{dam}}(d(\ell/\ell_c))} \int_{\tilde{d}=0}^{d(\ell/\ell_c)} h(\tilde{d}) \, d\tilde{d}} \quad (33) $$

where the parameters are $w$ and $f$ for cohesive model and $\ell$, $g_{dam}$ and $h$ for damage model, $d$ for TLS regularization and $\sigma$ as a general variable.

Making $\ell$ tend to 0 leads to $\sigma = \sqrt{\frac{2EY_c}{-g'_{\text{dam}}(0)}}$, $w = 0$ and $\sigma = \sigma_f$. Thus

$$ \sigma_f = \sqrt{\frac{2EY_c}{-g'_{\text{dam}}(0)}} \quad (34) $$

By defining dimensionless quantities

$$ \hat{\sigma} = \frac{\sigma}{\sigma_f}, \quad \hat{w} = \frac{w}{w_f}, \quad \hat{\phi} = \frac{\phi}{\ell_c}, \quad \hat{\ell} = \frac{\ell}{\ell_c}, \quad \lambda = \frac{2\sigma_f \ell_c}{Ew_f} \quad (35) $$

<table><caption>Table 1 Summary of cohesive and damage functions</caption><thead><tr><th>Cohesive model</th><th>Damage model</th></tr></thead><tbody><tr><td>$g_{CZM}(\alpha)$</td><td>$g_{dam}(d)$</td></tr><tr><td>$A_c(\alpha)$</td><td>$h(d)$</td></tr><tr><td></td><td>TLS regularization</td></tr><tr><td></td><td>$d(\phi/\ell_c), \ell_c$</td></tr></tbody></table>
These equations summarize the relations between the two models. They link the cohesive parameters $\alpha$, $g_{CZM}$ and $A_c$ to the damage ones $\ell$, $g_{dam}$ and $h$ and the TLS regularization function $d(\phi/\ell_c)$. As these equations hold during the whole damage process, that is for all possible values of $\alpha$ or $\ell$, it is possible to derive a TLS damage model from a cohesive one (and vice versa) by fixing four of the parameters. For example, if $d(\phi/\ell_c)$, $g_{dam}(d)$, $g_{CZM}(\alpha)$ and $A_c(\alpha)$ are known, the function $h(d)$ can be deduced (Table 1 presents a summary of the main parameters).

In what follows, we look for $h$ by only using the stress-opening function $f$ and the opening $w$, instead of fully describing $g_{CZM}$ and $A_c$ and $\alpha$ [see (6)]. By differentiating Eq. (26) with respect to time and by combining Eqs. (27) and (28), one obtains

$$ -\frac{1}{2k} \frac{g'_{CZM}(\alpha)}{g_{CZM}^2(\alpha)} \dot{\alpha} = \frac{1}{E} \dot{l} \left( \frac{1}{g_{dam}(d(\ell/\ell_c))} - 1 \right) \quad (29) $$

$$ 2Y_c \dot{l} \int_{\tilde{d}=0}^{d(\ell/\ell_c)} h(\tilde{d}) \, d\tilde{d} = - \frac{\sigma^2}{2k} \frac{g'_{CZM}(\alpha)}{g_{CZM}^2(\alpha)} \dot{\alpha} \quad (30) $$

Moreover, as $\sigma = kg_{CZM}w$, the system of equations is rewritten as

$$ \sigma = f(w) \quad (31) $$

$$ w = \frac{2\sigma}{E} \int_{0}^{\ell} \left( \frac{1}{g_{\text{dam}}(d(\phi))} - 1 \right) d\phi \quad (32) $$

$$ \sigma = \sqrt{2EY_c \frac{g_{\text{dam}}(d(\ell/\ell_c))}{1 - g_{\text{dam}}(d(\ell/\ell_c))} \int_{\tilde{d}=0}^{d(\ell/\ell_c)} h(\tilde{d}) \, d\tilde{d}} \quad (33) $$

where the parameters are $w$ and $f$ for cohesive model and $\ell$, $g_{dam}$ and $h$ for damage model, $d$ for TLS regularization and $\sigma$ as a general variable.

Making $\ell$ tend to 0 leads to $\sigma = \sqrt{\frac{2EY_c}{-g'_{\text{dam}}(0)}}$, $w = 0$ and $\sigma = \sigma_f$. Thus

$$ \sigma_f = \sqrt{\frac{2EY_c}{-g'_{\text{dam}}(0)}} \quad (34) $$

By defining dimensionless quantities

$$ \hat{\sigma} = \frac{\sigma}{\sigma_f}, \quad \hat{w} = \frac{w}{w_f}, \quad \hat{\phi} = \frac{\phi}{\ell_c}, \quad \hat{\ell} = \frac{\ell}{\ell_c}, \quad \lambda = \frac{2\sigma_f \ell_c}{Ew_f} \quad (35) $$

<table><caption>**Table 1 Summary of cohesive and damage functions**</caption><thead><tr><th>Cohesive model</th><th>Damage model</th></tr></thead><tbody><tr><td>$g_{CZM}(\alpha)$</td><td>$g_{dam}(d)$</td></tr><tr><td>$A_c(\alpha)$</td><td>$h(d)$</td></tr><tr><td></td><td>TLS regularization</td></tr><tr><td></td><td>$d(\phi/\ell_c), \ell_c$</td></tr></tbody></table>
and

$$I(\hat{\ell}) = \int_{0}^{\hat{\ell}} \frac{1 - g_{\text{dam}}(d(\hat{\phi}))}{g_{\text{dam}}(d(\hat{\phi}))} d\hat{\phi}, \qquad (36)$$

the system (31), (32) and (33) can be written in a dimensionless and condensed form

$$\hat{\sigma} = F(\lambda I \hat{\sigma}) \quad \text{and} \quad \int_{\tilde{\phi}=0}^{\hat{\ell}} h(\tilde{d}) \, d\tilde{\phi} = -g'_{\text{dam}}(0) \frac{1 - g_{\text{dam}}(d(\hat{\ell}))}{g_{\text{dam}}(d(\hat{\ell}))} \hat{\sigma}^2 \qquad (37)$$

These equations hold for all dimensionless position of the damage front in the bar $\hat{\ell}$, and therefore for any value of $\hat{\phi}$. Using $\hat{\phi}$ instead of $\ell$ gives a more general formulation, not constrained to the bar case.

$$\hat{\sigma} = F(\lambda I \hat{\sigma}) \quad \text{and} \quad \int_{0}^{\hat{\phi}} h(\hat{\phi}) \, d\hat{\phi} = -g'_{\text{dam}}(0) \frac{1 - g_{\text{dam}}(d(\hat{\phi}))}{g_{\text{dam}}(d(\hat{\phi}))} \hat{\sigma}^2 \qquad (38)$$

### Meaning of $\lambda$

The dimensionless quantity $\lambda$ includes the main parameters of both models and in particular the TLS characteristic length $\ell_c$. Different lengths have been defined to characterize the size of the fracture process zone [20–22] in quasi-brittle materials. For the specific case of cohesive models, Smith estimated its length $\ell_{coh}$ in [23] for a crack evolving in an infinite body and for a large range of cohesive laws:

$$\ell_{\text{coh}} \approx \frac{3 E w_f}{2 \sigma_f} \qquad (39)$$

Thus, the parameter $\lambda$ is close to the ratio between damage band width and cohesive zone length (see Fig. 7).

$$\lambda \approx \frac{3 \ell_c}{2 \ell_{\text{coh}}} \qquad (40)$$

<figure><img src="image_11.png" /><figcaption>Fig. 7 Classical shape of a damaged zone and lengths definitions.</figcaption></figure>
and

$$I(\hat{\ell}) = \int_{0}^{\hat{\ell}} \frac{1 - g_{\text{dam}}(d(\hat{\phi}))}{g_{\text{dam}}(d(\hat{\phi}))} d\hat{\phi}, \qquad (36)$$

the system (31), (32) and (33) can be written in a dimensionless and condensed form

$$\hat{\sigma} = F(\lambda I \hat{\sigma}) \quad \text{and} \quad \int_{\tilde{\phi}=0}^{\hat{\ell}} h(\tilde{d}) \, d\tilde{\phi} = -g'_{\text{dam}}(0) \frac{1 - g_{\text{dam}}(d(\hat{\ell}))}{g_{\text{dam}}(d(\hat{\ell}))} \hat{\sigma}^2 \qquad (37)$$

These equations hold for all dimensionless position of the damage front in the bar $\hat{\ell}$, and therefore for any value of $\hat{\phi}$. Using $\hat{\phi}$ instead of $\ell$ gives a more general formulation, not constrained to the bar case.

$$\hat{\sigma} = F(\lambda I \hat{\sigma}) \quad \text{and} \quad \int_{0}^{\hat{\phi}} h(\hat{\phi}) \, d\hat{\phi} = -g'_{\text{dam}}(0) \frac{1 - g_{\text{dam}}(d(\hat{\phi}))}{g_{\text{dam}}(d(\hat{\phi}))} \hat{\sigma}^2 \qquad (38)$$

### Meaning of $\lambda$

The dimensionless quantity $\lambda$ includes the main parameters of both models and in particular the TLS characteristic length $\ell_c$. Different lengths have been defined to characterize the size of the fracture process zone [20–22] in quasi-brittle materials. For the specific case of cohesive models, Smith estimated its length $\ell_{coh}$ in [23] for a crack evolving in an infinite body and for a large range of cohesive laws:

$$\ell_{\text{coh}} \approx \frac{3 E w_f}{2 \sigma_f} \qquad (39)$$

Thus, the parameter $\lambda$ is close to the ratio between damage band width and cohesive zone length (see Fig. 7).

$$\lambda \approx \frac{3 \ell_c}{2 \ell_{\text{coh}}} \qquad (40)$$

<figure><img src="image_12.png" /><figcaption>Fig. 7 Classical shape of a damaged zone and lengths definitions.</figcaption></figure>
The value of $\lambda$ will be chosen to be less than $\frac{1}{2}$. Indeed, as it is explained in "Appendix 2", the condition of uniqueness of the solution enforces values of $\lambda$ that are close or smaller than $\frac{1}{2}$. It corresponds to the choice of a small enough damage band width in comparison to the fracture process length. Writing $\lambda \lesssim \frac{1}{2}$ and Eq. (40) leads to

$$ \ell_c \lesssim \frac{1}{3} \ell_{\text{coh}} \qquad (41) $$

### From general to particular relations

These general relations are now derived in the case of a particular choice of $g_{\text{dam}}(d)$. Later, a particular form of $d(\hat{\phi})$ is assumed. Finally, the equations are particularized for two cohesive softening functions F: the linear and the bi-linear ones.

### A choice of damage function $g_{\text{dam}}$

A classical choice is made:

$$ g_{\text{dam}}(d) = 1 - d, \quad \text{and then} \quad \Phi_{\text{dam}}(\varepsilon, d) = \frac{1}{2}(1 - d)E\varepsilon^2 \qquad (42) $$

The admissibility condition (16) becomes

$$ h'(d) \geq 0 \qquad (43) $$

and

$$ I(\hat{\phi}) = \int_{0}^{\hat{\phi}} \frac{d}{1-d} d\hat{\phi} \qquad (44) $$

The system of Eq. (38) can be rewritten as follows.

$$ \hat{\sigma} = F(\lambda I \hat{\sigma}) \quad \text{and} \quad \int_{0}^{\hat{\phi}} h(\hat{\phi}) \, d\hat{\phi} = \frac{d(\hat{\phi})}{1 - d(\hat{\phi})} \hat{\sigma}^2 \qquad (45) $$

Moreover,

$$ \sigma_f = \sqrt{2EY_c} \qquad (46) $$

### A choice of TLS damage profile $d(\hat{\phi})$

We consider a parabolic damage law (see Fig. 6)

$$ d(\hat{\phi}) = 2\hat{\phi} - \hat{\phi}^2 \qquad (47) $$

and

$$ I(\hat{\phi}) = \frac{\hat{\phi}^2}{1 - \hat{\phi}} \qquad (48) $$
The value of $\lambda$ will be chosen to be less than $\frac{1}{2}$. Indeed, as it is explained in "Appendix 2", the condition of uniqueness of the solution enforces values of $\lambda$ that are close or smaller than $\frac{1}{2}$. It corresponds to the choice of a small enough damage band width in comparison to the fracture process length. Writing $\lambda \lesssim \frac{1}{2}$ and Eq. (40) leads to

$$ \ell_c \lesssim \frac{1}{3} \ell_{\text{coh}} \qquad (41) $$

### From general to particular relations

These general relations are now derived in the case of a particular choice of $g_{\text{dam}}(d)$. Later, a particular form of $d(\hat{\phi})$ is assumed. Finally, the equations are particularized for two cohesive softening functions F: the linear and the bi-linear ones.

### A choice of damage function $g_{\text{dam}}$

A classical choice is made:

$$ g_{\text{dam}}(d) = 1 - d, \quad \text{and then} \quad \Phi_{\text{dam}}(\varepsilon, d) = \frac{1}{2}(1 - d)E\varepsilon^2 \qquad (42) $$

The admissibility condition (16) becomes

$$ h'(d) \geq 0 \qquad (43) $$

and

$$ I(\hat{\phi}) = \int_{0}^{\hat{\phi}} \frac{d}{1-d} d\hat{\phi} \qquad (44) $$

The system of Eq. (38) can be rewritten as follows.

$$ \hat{\sigma} = F(\lambda I \hat{\sigma}) \quad \text{and} \quad \int_{0}^{\hat{\phi}} h(\hat{\phi}) \, d\hat{\phi} = \frac{d(\hat{\phi})}{1 - d(\hat{\phi})} \hat{\sigma}^2 \qquad (45) $$

Moreover,

$$ \sigma_f = \sqrt{2EY_c} \qquad (46) $$

### A choice of TLS damage profile $d(\hat{\phi})$

We consider a parabolic damage law (see Fig. 6)

$$ d(\hat{\phi}) = 2\hat{\phi} - \hat{\phi}^2 \qquad (47) $$

and

$$ I(\hat{\phi}) = \frac{\hat{\phi}^2}{1 - \hat{\phi}} \qquad (48) $$
The system (45) writes now

$$ \hat{\sigma} = F\left(\frac{\lambda\hat{\phi}^2}{1-\hat{\phi}}\hat{\sigma}\right) \quad \text{and} \quad \int_0^{\hat{\phi}} h(\hat{\phi}) \, d\hat{\phi} = \frac{2\hat{\phi}-\hat{\phi}^2}{(\hat{\phi}-1)^2}\hat{\sigma}^2 \qquad (49) $$

This choice is motivated by a damage distribution obtained from non-local damage and fracture equivalence [18], damage profiles obtained by lattice model simulations [24] and their similarity with some acoustic emission profiles [25, 26].

## Two examples of cohesive laws

Two particular cases of cohesive stress-opening functions are derived from previous relationships: the linear and the bi-linear ones.

### Linear cohesive law

We consider first a linear cohesive law $F(\hat{w}) = 1 - \hat{w}$ (see Fig. 8). The softening function is obtained by (49).

$$ \hat{\sigma}(\hat{\phi}) = \frac{1 - \hat{\phi}}{1 - \hat{\phi} + \lambda\hat{\phi}^2} \quad \text{and} \quad \int_0^{\hat{\phi}} h(\hat{\phi}) \, d\hat{\phi} = \frac{2\hat{\phi} - \hat{\phi}^2}{(1 - \hat{\phi} + \lambda\hat{\phi}^2)^2} \qquad (50) $$

### Bi-linear cohesive law

The bi-linear cohesive law is considered as it is one of the most popular laws to describe concrete [27–32]. It is presented in Fig. 3a. The method to obtain $h$ is the same as previously. The result is a discontinuous and increasing $h$ function. Corresponding strain-stress curves are shown in Fig. 3b. All of them present a discontinuity that appears as a linear zone where strain and stress increase while damage value remains constant.

<figure><img src="image_12.png" /><figcaption>Fig. 8 Cohesive linear law and TLS local behavior.</figcaption></figure>
The system (45) writes now

$$ \hat{\sigma} = F\left(\frac{\lambda\hat{\phi}^2}{1-\hat{\phi}}\hat{\sigma}\right) \quad \text{and} \quad \int_0^{\hat{\phi}} h(\hat{\phi}) \, d\hat{\phi} = \frac{2\hat{\phi}-\hat{\phi}^2}{(\hat{\phi}-1)^2}\hat{\sigma}^2 \qquad (49) $$

This choice is motivated by a damage distribution obtained from non-local damage and fracture equivalence [18], damage profiles obtained by lattice model simulations [24] and their similarity with some acoustic emission profiles [25, 26].

## Two examples of cohesive laws

Two particular cases of cohesive stress-opening functions are derived from previous relationships: the linear and the bi-linear ones.

### Linear cohesive law

We consider first a linear cohesive law $F(\hat{w}) = 1 - \hat{w}$ (see Fig. 8). The softening function is obtained by (49).

$$ \hat{\sigma}(\hat{\phi}) = \frac{1 - \hat{\phi}}{1 - \hat{\phi} + \lambda\hat{\phi}^2} \quad \text{and} \quad \int_0^{\hat{\phi}} h(\hat{\phi}) \, d\hat{\phi} = \frac{2\hat{\phi} - \hat{\phi}^2}{(1 - \hat{\phi} + \lambda\hat{\phi}^2)^2} \qquad (50) $$

### Bi-linear cohesive law

The bi-linear cohesive law is considered as it is one of the most popular laws to describe concrete [27–32]. It is presented in Fig. 3a. The method to obtain $h$ is the same as previously. The result is a discontinuous and increasing $h$ function. Corresponding strain-stress curves are shown in Fig. 3b. All of them present a discontinuity that appears as a linear zone where strain and stress increase while damage value remains constant.

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

a Cohesive linear law.

b TLS equivalent local behavior for different $\ell_c$ values. Increasing values of $\ell_c$ are indicated by the arrow.

<caption>Fig. 8 Cohesive linear law and TLS local behavior.</caption>
This plateau of damage is caused by the discontinuity of slope of the cohesive behavior. Details of the calculations are given in the "Appendix 1". Some conditions on the choice of the cohesive and TLS parameters are analyzed in "Appendix 2".

## Results and discussion

The main numerical issues for implementing the TLS approach are

* the necessity of representing jumps in displacement, for which classical X-FEM enrichment is used [33];

* the calculation of non-local energy release rate $\bar{Y}$, performed by the resolution of a variational problem described in [13];

* the propagation of the damage front, that is performed by an explicit algorithm with prediction described in "Appendix 3".

* The only difference with what is described in [13] is this last point. Concerning damage initiation for an initial sane structure, it is not numerically possible to initiate damage at a single material point. So, it is chosen to introduce a very small damaged zone of about one element size. Moreover, let us precise that the potential given in (12) was symmetric in tension and compression for the sake of simplicity. Whereas, concrete is known to have a dissymmetric behavior. Thus, it is replaced by the following one.

$$ \Phi_{\text{dam}}(\varepsilon, d) = \mu(1 - \alpha_i d)\varepsilon_i^2 + \frac{\lambda}{2}(1 - \alpha d)\text{tr}(\varepsilon)^2 \quad (51) $$

where $\lambda$ and $\mu$ are the Lamé elasticity coefficients, $\varepsilon_i$ the eigenvalues of the strain tensor and

$$ \left\{ \begin{array}{ll} \alpha_i = 0 & \text{if } \varepsilon_i < 0 \\ = 1 & \text{if } \varepsilon_i \ge 0 \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{ll} \alpha = 0 & \text{if } \text{tr}(\varepsilon) < 0 \\ = 1 & \text{if } \text{tr}(\varepsilon) \ge 0 \end{array} \right. \quad (52) $$

It is important to note that the h function derived in the previous section from a one-dimensional equivalence between cohesive model and TLS is used here in bi-dimensional simulations. The results are compared to cohesive simulations found in the literature. The results are compared in terms of local behavior in the "Simulations and results—influence of $\ell_c$" and global behavior in "Simulations and results—comparison under mixed-mode load". In this case, the non-rectilinear crack path and the global force-CMOD curves are compared.

### Simulations and results—influence of $\ell_c$

Simulations concern a three point bending test (see Fig. 9) of concrete which Young modulus is 37 GPa and Poisson ratio is 0.21. The equivalent cohesive zone properties of the TLS simulations performed are $w_f = 10^{-4}$ m, $w_1 = 4 \cdot 10^{-5}$ m, $\sigma_f = 3.5$ GPa and $\sigma_k = 1.0$ GPa. An example of the results is plotted in Fig. 10. All simulations are performed with the same ratio $l_{\text{mesh}}/\ell_c = 1/20$ where $l_{\text{mesh}}$ is the characteristic length of the mesh over the center of the beam where the damage band should evolve. The local TLS behavior equivalent to those CZM parameters is plotted in Fig. 3b for four values of $\ell_c$: 15, 30, 60 and 100 mm.
This plateau of damage is caused by the discontinuity of slope of the cohesive behavior. Details of the calculations are given in the "Appendix 1". Some conditions on the choice of the cohesive and TLS parameters are analyzed in "Appendix 2".

## Results and discussion

The main numerical issues for implementing the TLS approach are

* the necessity of representing jumps in displacement, for which classical X-FEM enrichment is used [33];

* the calculation of non-local energy release rate $\bar{Y}$, performed by the resolution of a variational problem described in [13];

* the propagation of the damage front, that is performed by an explicit algorithm with prediction described in "Appendix 3".

* The only difference with what is described in [13] is this last point. Concerning damage initiation for an initial sane structure, it is not numerically possible to initiate damage at a single material point. So, it is chosen to introduce a very small damaged zone of about one element size. Moreover, let us precise that the potential given in (12) was symmetric in tension and compression for the sake of simplicity. Whereas, concrete is known to have a dissymmetric behavior. Thus, it is replaced by the following one.

$$ \Phi_{\text{dam}}(\varepsilon, d) = \mu(1 - \alpha_i d)\varepsilon_i^2 + \frac{\lambda}{2}(1 - \alpha d)\text{tr}(\varepsilon)^2 \quad (51) $$

where $\lambda$ and $\mu$ are the Lamé elasticity coefficients, $\varepsilon_i$ the eigenvalues of the strain tensor and

$$ \left\{ \begin{array}{ll} \alpha_i = 0 & \text{if } \varepsilon_i < 0 \\ = 1 & \text{if } \varepsilon_i \ge 0 \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{ll} \alpha = 0 & \text{if } \text{tr}(\varepsilon) < 0 \\ = 1 & \text{if } \text{tr}(\varepsilon) \ge 0 \end{array} \right. \quad (52) $$

It is important to note that the h function derived in the previous section from a one-dimensional equivalence between cohesive model and TLS is used here in bi-dimensional simulations. The results are compared to cohesive simulations found in the literature. The results are compared in terms of local behavior in the "Simulations and results—influence of $\ell_c$" and global behavior in "Simulations and results—comparison under mixed-mode load". In this case, the non-rectilinear crack path and the global force-CMOD curves are compared.

### Simulations and results—influence of $\ell_c$

Simulations concern a three point bending test (see Fig. 9) of concrete which Young modulus is 37 GPa and Poisson ratio is 0.21. The equivalent cohesive zone properties of the TLS simulations performed are $w_f = 10^{-4}$ m, $w_1 = 4 \cdot 10^{-5}$ m, $\sigma_f = 3.5$ GPa and $\sigma_k = 1.0$ GPa. An example of the results is plotted in Fig. 10. All simulations are performed with the same ratio $l_{\text{mesh}}/\ell_c = 1/20$ where $l_{\text{mesh}}$ is the characteristic length of the mesh over the center of the beam where the damage band should evolve. The local TLS behavior equivalent to those CZM parameters is plotted in Fig. 3b for four values of $\ell_c$: 15, 30, 60 and 100 mm.
<figure><img src="image_13.png" /><figcaption>**Fig. 9** Three point bending (TPB) test.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 10** Damage over the deformed specimen at a given step of the simulation. Blue is undamaged material. Red is damage close to 1. The white locates the crack.</figcaption></figure>

**Global behavior: three point bending tests** It is interesting to observe global CMOD-force curves shown in Fig. 11. The CMOD is here defined as the relative displacement of two points symmetrically located at the bottom face of the beam and separated by a distance equal to the beam depth. The global behavior is very close for the four values of $l_c$. The main difference is the value of the maximum load. However this difference is slight: the maximum gap to the mean value 24.6 kN is less than 3%.

**Damage zone shape at maximum load** It is interesting to analyze the shape of the damaged zone at the maximum load. For different $l_c$ values, the position of the damage front is drawn in Fig. 12. The width of that zone is wider as $l_c$ is bigger. It is consistent with the

<figure><img src="image_15.png" /><figcaption>**Fig. 11** Global behavior of TLS simulations of TPB test for different $l_c$. Raw data is drawn.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 9** Three point bending (TPB) test.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 10** Damage over the deformed specimen at a given step of the simulation. Blue is undamaged material. Red is damage close to 1. The white locates the crack.</figcaption></figure>

**Global behavior: three point bending tests** It is interesting to observe global CMOD-force curves shown in Fig. 11. The CMOD is here defined as the relative displacement of two points symmetrically located at the bottom face of the beam and separated by a distance equal to the beam depth. The global behavior is very close for the four values of $l_c$. The main difference is the value of the maximum load. However this difference is slight: the maximum gap to the mean value 24.6 kN is less than 3%.

**Damage zone shape at maximum load** It is interesting to analyze the shape of the damaged zone at the maximum load. For different $l_c$ values, the position of the damage front is drawn in Fig. 12. The width of that zone is wider as $l_c$ is bigger. It is consistent with the

<figure><img src="image_16.png" /><figcaption>**Fig. 11** Global behavior of TLS simulations of TPB test for different $l_c$. Raw data is drawn.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 12 Position of the damage front at the maximum load for different $\ell_c$ values.</figcaption></figure>

fact that in TLS, the width of the fracture process zone is driven by $\ell_c$. For example, the width of the damaged zone corresponding to a propagating crack is $2\ell_c$. The length of the damaged zone ahead of the crack tip is significantly the same for all $\ell_c$. The depth is in fact driven by the equivalent cohesive behavior. We can conclude in which concerns the fracture process zone that its length is probably the same that in cohesive simulations and its width, neglected in CZM, is now driven by TLS parameter $\ell_c$.

*Agreement of the damage model with local cohesive behavior* Cohesive theoretical local behavior is compared to TLS simulated one. As the cohesive zone model is concentrated over a line and has no width, it is necessary to determine which TLS quantities, denoted $\sigma^{\text{TLS}}$ and $w^{\text{TLS}}$, are compared to CZM ones, that is $\sigma_{xx}$ and $w$. All over the evolving zone of the damage band, the stress and the opening are computed for different positions $y$. Let $y_{\text{max}}$ be the maximum value of $y$ where $d > 0$. For $y < y_{\text{max}}$, about 15 positions are analyzed for the four $\ell_c$ values.

It is chosen to take the value of the stress $\sigma_{xx}$ on the damage front, for a given value of $y$. By denoting $l(y)$ the width of the damage band at a given $y$ and $S_y$ the points of the damage band where $d(x, y) > 0$ and $y$ is fixed (see Fig. 13),

$$ \sigma^{\text{TLS}}(y) = \sigma_{xx} \left( x = \frac{l(y)}{2}, y \right) \qquad (53) $$

Concerning the opening, it is necessary for TLS to consider the whole width of the damaged zone. So it is chosen to compare displacement, perpendicular to the crack lips, at the boundaries of $S_y$. Thus, an additional term has to be added to cohesive opening to have comparable quantities:

$$ w^{\text{CZM}}(y) + \int_{S_y} \varepsilon_{xx}^{\text{CZM}}(x, y) \, dx \quad \text{can be compared to} \quad \int_{S_y} \varepsilon_{xx}^{\text{TLS}}(x, y) \, dx \qquad (54) $$

<figure><img src="image_17.png" /><figcaption>Fig. 13 Definition of $S_y$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 12 Position of the damage front at the maximum load for different $\ell_c$ values.</figcaption></figure>

fact that in TLS, the width of the fracture process zone is driven by $\ell_c$. For example, the width of the damaged zone corresponding to a propagating crack is $2\ell_c$. The length of the damaged zone ahead of the crack tip is significantly the same for all $\ell_c$. The depth is in fact driven by the equivalent cohesive behavior. We can conclude in which concerns the fracture process zone that its length is probably the same that in cohesive simulations and its width, neglected in CZM, is now driven by TLS parameter $\ell_c$.

*Agreement of the damage model with local cohesive behavior* Cohesive theoretical local behavior is compared to TLS simulated one. As the cohesive zone model is concentrated over a line and has no width, it is necessary to determine which TLS quantities, denoted $\sigma^{\text{TLS}}$ and $w^{\text{TLS}}$, are compared to CZM ones, that is $\sigma_{xx}$ and $w$. All over the evolving zone of the damage band, the stress and the opening are computed for different positions $y$. Let $y_{\text{max}}$ be the maximum value of $y$ where $d > 0$. For $y < y_{\text{max}}$, about 15 positions are analyzed for the four $\ell_c$ values.

It is chosen to take the value of the stress $\sigma_{xx}$ on the damage front, for a given value of $y$. By denoting $l(y)$ the width of the damage band at a given $y$ and $S_y$ the points of the damage band where $d(x, y) > 0$ and $y$ is fixed (see Fig. 13),

$$ \sigma^{\text{TLS}}(y) = \sigma_{xx} \left( x = \frac{l(y)}{2}, y \right) \qquad (53) $$

Concerning the opening, it is necessary for TLS to consider the whole width of the damaged zone. So it is chosen to compare displacement, perpendicular to the crack lips, at the boundaries of $S_y$. Thus, an additional term has to be added to cohesive opening to have comparable quantities:

$$ w^{\text{CZM}}(y) + \int_{S_y} \varepsilon_{xx}^{\text{CZM}}(x, y) \, dx \quad \text{can be compared to} \quad \int_{S_y} \varepsilon_{xx}^{\text{TLS}}(x, y) \, dx \qquad (54) $$

<figure><img src="image_18.png" /><figcaption>Fig. 13 Definition of $S_y$.</figcaption></figure>
where $\int_{S_y} \varepsilon_{xx}^{CZM}(x, y) dx$ is the displacement caused by the deformation of the sane matter around the cohesive crack. We can write

$$ \varepsilon_{xx}^{CZM} = (1 - d)\varepsilon_{xx}^{TLS} \qquad (55) $$

Thus we define:

$$ w^{TLS}(y) = \int_{S_y} \varepsilon_{xx}(x, y) dx - \int_{S_y} (1-d)\varepsilon_{xx}(x, y) dx = \int_{S_y} d(x, y)\varepsilon_{xx}(x, y) dx \qquad (56) $$

The opening-stress curves obtained for different characteristic lengths $l_c$ are shown in Fig. 14. It can be observed that scatter plot is close to the equivalent cohesive law for small values of $l_c$. The TLS model behavior around the damaged zone tends to the cohesive one as $l_c$ gets smaller. Simulation for much smaller $l_c$ values is difficult. Indeed, as $l_c$ goes to zero, displacement and damage gradients become extremely high and are difficult to capture numerically.

## Simulations and results—comparison under mixed-mode load

In order to verify the consistency of global results for both models, a mixed-mode test is performed on single edge notched specimens. As the crack path is not a priori known, the cohesive model is coupled with classical fracture model under MTS criterion. Concerning TLS, no modification or enhancement is needed. The test, initially described in [34], is simulated with a cohesive model in [35]. The damage front evolution during the simulation is shown in Fig. 15. Note that the gradient of damage is not parallel to the boundary of the domain, whereas it is for phase-field models [36].

A first TLS simulation is performed with the smallest initial damaged zone considered here, that is of radius 0.10 $l_c$ (see Fig. 16). It can be observed that the cohesive crack path is close to the TLS crack lips even if some differences exist. Concerning the load-CMOD curve, there is an overestimation of the maximum load and of the end of the post-peak curve. The experimental envelope given by [35] of the tests described in [34] is also drawn.

<figure><img src="image_18.png" /><figcaption>Fig. 14 Stress distribution in the process zone for TLS damage model with different $l_c$. Abscissa is an equivalent cohesive displacement jump. The continuous black line is the cohesive bi-linear model.</figcaption></figure>
where $\int_{S_y} \varepsilon_{xx}^{CZM}(x, y) dx$ is the displacement caused by the deformation of the sane matter around the cohesive crack. We can write

$$ \varepsilon_{xx}^{CZM} = (1 - d)\varepsilon_{xx}^{TLS} \qquad (55) $$

Thus we define:

$$ w^{TLS}(y) = \int_{S_y} \varepsilon_{xx}(x, y) dx - \int_{S_y} (1-d)\varepsilon_{xx}(x, y) dx = \int_{S_y} d(x, y)\varepsilon_{xx}(x, y) dx \qquad (56) $$

The opening-stress curves obtained for different characteristic lengths $l_c$ are shown in Fig. 14. It can be observed that scatter plot is close to the equivalent cohesive law for small values of $l_c$. The TLS model behavior around the damaged zone tends to the cohesive one as $l_c$ gets smaller. Simulation for much smaller $l_c$ values is difficult. Indeed, as $l_c$ goes to zero, displacement and damage gradients become extremely high and are difficult to capture numerically.

## Simulations and results—comparison under mixed-mode load

In order to verify the consistency of global results for both models, a mixed-mode test is performed on single edge notched specimens. As the crack path is not a priori known, the cohesive model is coupled with classical fracture model under MTS criterion. Concerning TLS, no modification or enhancement is needed. The test, initially described in [34], is simulated with a cohesive model in [35]. The damage front evolution during the simulation is shown in Fig. 15. Note that the gradient of damage is not parallel to the boundary of the domain, whereas it is for phase-field models [36].

A first TLS simulation is performed with the smallest initial damaged zone considered here, that is of radius 0.10 $l_c$ (see Fig. 16). It can be observed that the cohesive crack path is close to the TLS crack lips even if some differences exist. Concerning the load-CMOD curve, there is an overestimation of the maximum load and of the end of the post-peak curve. The experimental envelope given by [35] of the tests described in [34] is also drawn.

<figure><img src="image_19.png" /><figcaption>**Fig. 14** Stress distribution in the process zone for TLS damage model with different $l_c$. Abscissa is an equivalent cohesive displacement jump. The continuous black line is the cohesive bi-linear model.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 15** Damage evolution in TLS simulation.</figcaption></figure>

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

a Zoom on the ligament of the specimen. TLS damage iso-values are drawn. Totally damaged zone is white. For comparison, CZM result from [35] is drawn.

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

b Load-CMOD curve. Plain curve is TLS (raw data), dashed one is CZM and dotted ones are the experimental envelope.

<caption>**Fig. 16** Comparison between TLS and CZM results. Initial damage radius in TLS simulation is equal to 0.10 $\ell_c$.</caption>

The cohesive model used in [35] forces the crack to be in the continuity of the notch. In the TLS model, no particular assumption is done on the damage zone path. This could explain the difference in the results presented in Fig. 16. To check this assumption, another TLS simulation is presented in Fig. 17, where initial damage is bigger than previously. Here, its radius is $0.75\ell_c$. The CZM crack path is here overlapped with the TLS crack lips. Furthermore, load-CMOD curves are coincident until about the first half of the post-peak curve. Below that point, the TLS curve remains in the experimental envelope whereas cohesive one does not. It can be concluded that TLS simulation without an important initial damage give results slightly different from the CZM ones, but whose load-CMOD curve is closer to experimental one. Forcing similar crack paths by introducing a more important initial damage leads to very close global results.
<figure><img src="image_20.png" /><figcaption>**Fig. 15** Damage evolution in TLS simulation.</figcaption></figure>

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

a Zoom on the ligament of the specimen. TLS damage iso-values are drawn. Totally damaged zone is white. For comparison, CZM result from [35] is drawn.

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

b Load-CMOD curve. Plain curve is TLS (raw data), dashed one is CZM and dotted ones are the experimental envelope.

<caption>**Fig. 16** Comparison between TLS and CZM results. Initial damage radius in TLS simulation is equal to 0.10 $\ell_c$.</caption>

The cohesive model used in [35] forces the crack to be in the continuity of the notch. In the TLS model, no particular assumption is done on the damage zone path. This could explain the difference in the results presented in Fig. 16. To check this assumption, another TLS simulation is presented in Fig. 17, where initial damage is bigger than previously. Here, its radius is $0.75\ell_c$. The CZM crack path is here overlapped with the TLS crack lips. Furthermore, load-CMOD curves are coincident until about the first half of the post-peak curve. Below that point, the TLS curve remains in the experimental envelope whereas cohesive one does not. It can be concluded that TLS simulation without an important initial damage give results slightly different from the CZM ones, but whose load-CMOD curve is closer to experimental one. Forcing similar crack paths by introducing a more important initial damage leads to very close global results.
<figure><img src="image_22.png" /></figure>

a Zoom on the ligament of the specimen. TLS damage iso-values are drawn. Totally damaged zone is white. For comparison, CZM result from [35] is drawn.

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

b Load-CMOD curve. Plain curve is TLS (raw data), dashed one is CZM and dotted ones are the experimental envelope.

<caption>**Fig. 17** Comparison between TLS and CZM results. Initial damage radius in TLS simulation is equal to 0.75 *l*<sub>c</sub>.</caption>

**Conclusions**

It has been shown that it is possible within the TLS framework to derive an equiva-
lent one-dimensional damage behavior from any cohesive model. The damage model
depends on the characteristic length l̅c. The damage behavior tends to the cohesive one
when l̅c is close to zero. Furthermore, global response of structures is independent on
l̅c for a given one-dimensional equivalent cohesive law. CMOD-force curves and crack
path are the same for both models.

The TLS model is able to provide results of the same quality as the cohesive zone
model, that is well-known for its accuracy. Besides, it has already proven its capability to
determine accurate crack paths without supplementary hypothesis or models. Moreover,
it is able to represent branching, coalescence and initiation of damaged zones [5, 13].
Even if this paper only deals with quasi-brittle materials, as damaged zone is concen-
trated over a damage band, recent developments couple local and non-local damage and
allows diffuse damage prior to localization [11].

**Abbreviations**

TLS: thick level set; CZM: cohesive zone model; LEFM: linear elasticity fracture mechanics; FPZ: fracture process zone; MTS: maximum Tangential Stress; TPB: three point bending; CMOD: crack mouth opening displacement.

**Authors' contributions**

NM and CS provided the main idea for a comparison between TLS and CZM models and performed one-dimen-
sional analysis. APG carried out the bi-dimensional simulations and analysis. All authors read and approved the final
manuscript.

**Author details**

**Acknowledgements**

We would like to thank Alexis Salzman and Nicolas Chevaugeon for their contributions in the implementation of the TLS.

**Compliance with ethical guidelines**

**Competing interests**

The authors declare that they have no competing interests.

<footnote n="1">¹ GeM, UMR CNRS 6183, Ecole Centrale de Nantes, 1 rue de la Noë, 44321 Nantes, France. ² IMSIA, UMR CNRS 9219, EdF, Av. Charles de Gaulle, 92141 Clamart, France.</footnote>
<figure><img src="image_23.png" /></figure>

a Zoom on the ligament of the specimen. TLS damage iso-values are drawn. Totally damaged zone is white. For comparison, CZM result from [35] is drawn.

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

b Load-CMOD curve. Plain curve is TLS (raw data), dashed one is CZM and dotted ones are the experimental envelope.

<caption>**Fig. 17** Comparison between TLS and CZM results. Initial damage radius in TLS simulation is equal to 0.75 *l*<sub>c</sub>.</caption>

**Conclusions**

It has been shown that it is possible within the TLS framework to derive an equiva-
lent one-dimensional damage behavior from any cohesive model. The damage model
depends on the characteristic length l̅c. The damage behavior tends to the cohesive one
when l̅c is close to zero. Furthermore, global response of structures is independent on
l̅c for a given one-dimensional equivalent cohesive law. CMOD-force curves and crack
path are the same for both models.

The TLS model is able to provide results of the same quality as the cohesive zone
model, that is well-known for its accuracy. Besides, it has already proven its capability to
determine accurate crack paths without supplementary hypothesis or models. Moreover,
it is able to represent branching, coalescence and initiation of damaged zones [5, 13].
Even if this paper only deals with quasi-brittle materials, as damaged zone is concen-
trated over a damage band, recent developments couple local and non-local damage and
allows diffuse damage prior to localization [11].

**Abbreviations**

TLS: thick level set; CZM: cohesive zone model; LEFM: linear elasticity fracture mechanics; FPZ: fracture process zone; MTS: maximum Tangential Stress; TPB: three point bending; CMOD: crack mouth opening displacement.

**Authors' contributions**

NM and CS provided the main idea for a comparison between TLS and CZM models and performed one-dimen-
sional analysis. APG carried out the bi-dimensional simulations and analysis. All authors read and approved the final
manuscript.

**Author details**

**Acknowledgements**

We would like to thank Alexis Salzman and Nicolas Chevaugeon for their contributions in the implementation of the TLS.

**Compliance with ethical guidelines**

**Competing interests**

The authors declare that they have no competing interests.

<footnote n="1">¹ GeM, UMR CNRS 6183, Ecole Centrale de Nantes, 1 rue de la Noë, 44321 Nantes, France. ² IMSIA, UMR CNRS 9219, EdF, Av. Charles de Gaulle, 92141 Clamart, France.</footnote>
# Appendix

## Appendix 1: Bi-linear cohesive law

We define the following dimensionless quantities

$$ \hat{w}_1 = w_1/w_f, \quad \hat{w}_k = w_k/w_f, \quad \hat{\sigma}_k = \sigma_k/\sigma_f \qquad (57) $$

where $w_k = f^{-1}(\sigma_k)$. By construction, we have the following relation

$$ \hat{\sigma}_k = 1 - \hat{w}_k/\hat{w}_1 \qquad (58) $$

The tension opening relation is

$$ \hat{\sigma} = f(\hat{w}) = 1 - \hat{w}/\hat{w}_1, \quad 0 \le \hat{w} \le \hat{w}_k \qquad (59) $$

$$ \hat{\sigma} = f(\hat{w}) = \frac{\hat{\sigma}_k}{1 - \hat{w}_k} (1 - \hat{w}), \quad \hat{w}_k \le \hat{w} \qquad (60) $$

For the first and second part of the law, we solve (49) for each part of the tension-opening function:

$$ \hat{\sigma} = 1 - \lambda \hat{\sigma} I / \hat{w}_1 = \frac{1}{1 + (\lambda / \hat{w}_1) I} \Rightarrow H = \frac{2 \hat{\phi} - \hat{\phi}^2}{(1 - \hat{\phi} + (\lambda / \hat{w}_1) \hat{\phi}^2)^2} \qquad (61) $$

$$ \hat{\sigma} = \frac{\hat{\sigma}_k}{1 - \hat{w}_k} (1 - \lambda \hat{\sigma} I) = \frac{\hat{\sigma}_k / (1 - \hat{w}_k)}{1 + \hat{\sigma}_k / (1 - \hat{w}_k) \lambda I} \Rightarrow H = \left( \frac{\hat{\sigma}_k}{1 - \hat{w}_k} \right)^2 \frac{2 \hat{\phi} - \hat{\phi}^2}{(1 - \hat{\phi} + \frac{\lambda \hat{\sigma}_k}{1 - \hat{w}_k} \hat{\phi}^2)^2} \qquad (62) $$

where

$$ H(\hat{\phi}) = \int_{0}^{\hat{\phi}} h(\hat{\phi}) \, d\hat{\phi} \qquad (63) $$

It remains to determine for which value $\hat{\phi}_k$ of $\hat{\phi}$ we move from the first to the second function. This value must ensure the continuity of $H$. Thus

$$ \hat{\phi}_k = \frac{a}{2} \left( \sqrt{1 + 4/a} - 1 \right) \qquad (64) $$

with $a = \frac{\hat{w}_k}{\lambda \hat{\sigma}_k}$ and $\hat{\phi}_k$ between 0 and 1.

## Appendix 2: Conditions on the parameters

It has been shown that $h$ must verify condition (43) to ensure uniqueness of the solution. In the particular case of a linear cohesive model, uniqueness condition writes

$$ \lambda \le \frac{1}{2}, \quad it \qquad (65) $$

If $\lambda > 1/2$, there exists a critical value $\hat{\phi}_k$ above which $h' < 0$ (corresponding to the red part of the curve in Fig. 5b). It is difficult to obtain an explicit form of its value.

Concerning bi-linear model, two parameters are defined
# Appendix

## Appendix 1: Bi-linear cohesive law

We define the following dimensionless quantities

$$ \hat{w}_1 = w_1/w_f, \quad \hat{w}_k = w_k/w_f, \quad \hat{\sigma}_k = \sigma_k/\sigma_f \qquad (57) $$

where $w_k = f^{-1}(\sigma_k)$. By construction, we have the following relation

$$ \hat{\sigma}_k = 1 - \hat{w}_k/\hat{w}_1 \qquad (58) $$

The tension opening relation is

$$ \hat{\sigma} = f(\hat{w}) = 1 - \hat{w}/\hat{w}_1, \quad 0 \le \hat{w} \le \hat{w}_k \qquad (59) $$

$$ \hat{\sigma} = f(\hat{w}) = \frac{\hat{\sigma}_k}{1 - \hat{w}_k} (1 - \hat{w}), \quad \hat{w}_k \le \hat{w} \qquad (60) $$

For the first and second part of the law, we solve (49) for each part of the tension-opening function:

$$ \hat{\sigma} = 1 - \lambda \hat{\sigma} I / \hat{w}_1 = \frac{1}{1 + (\lambda / \hat{w}_1) I} \Rightarrow H = \frac{2 \hat{\phi} - \hat{\phi}^2}{(1 - \hat{\phi} + (\lambda / \hat{w}_1) \hat{\phi}^2)^2} \qquad (61) $$

$$ \hat{\sigma} = \frac{\hat{\sigma}_k}{1 - \hat{w}_k} (1 - \lambda \hat{\sigma} I) = \frac{\hat{\sigma}_k / (1 - \hat{w}_k)}{1 + \hat{\sigma}_k / (1 - \hat{w}_k) \lambda I} \Rightarrow H = \left( \frac{\hat{\sigma}_k}{1 - \hat{w}_k} \right)^2 \frac{2 \hat{\phi} - \hat{\phi}^2}{(1 - \hat{\phi} + \frac{\lambda \hat{\sigma}_k}{1 - \hat{w}_k} \hat{\phi}^2)^2} \qquad (62) $$

where

$$ H(\hat{\phi}) = \int_{0}^{\hat{\phi}} h(\hat{\phi}) \, d\hat{\phi} \qquad (63) $$

It remains to determine for which value $\hat{\phi}_k$ of $\hat{\phi}$ we move from the first to the second function. This value must ensure the continuity of $H$. Thus

$$ \hat{\phi}_k = \frac{a}{2} \left( \sqrt{1 + 4/a} - 1 \right) \qquad (64) $$

with $a = \frac{\hat{w}_k}{\lambda \hat{\sigma}_k}$ and $\hat{\phi}_k$ between 0 and 1.

## Appendix 2: Conditions on the parameters

It has been shown that $h$ must verify condition (43) to ensure uniqueness of the solution. In the particular case of a linear cohesive model, uniqueness condition writes

$$ \lambda \le \frac{1}{2}, \quad it \qquad (65) $$

If $\lambda > 1/2$, there exists a critical value $\hat{\phi}_k$ above which $h' < 0$ (corresponding to the red part of the curve in Fig. 5b). It is difficult to obtain an explicit form of its value.

Concerning bi-linear model, two parameters are defined
$$ \lambda_1 = \frac{2\sigma_f l_c}{Ew_1}, \quad \lambda_2 = \frac{2\sigma_2 l_c}{Ew_f} = \frac{2\sigma_k l_c}{E(w_f - w_k)} \quad (66) $$

The first one, $\lambda_1$, characterizes the first slope of the cohesive law and $\lambda_2$ the second one. Let us consider $\lambda_1 > \frac{1}{2}$. Then, as seen above, there exists a critical value of $\hat{\phi}$ that must not be exceeded. This value corresponds to a point of the cohesive law that must not be reached. If the position of the slope-discontinuity point is located before the critical point, the condition $h' \ge 0$ is respected for the first part of the cohesive law. Finally, choosing $\lambda_2 \le \frac{1}{2}$ ensures the condition for the whole local behavior. Some examples of equivalent one-dimensional strain-stress curves for a bi-linear cohesive law are shown in Fig. 8b.

## Appendix 3: Damage evolution in TLS and load factor

We rewrite the evolution criterion (18) by defining the averaged quantities

$$ \bar{Y} = \frac{\int_{\phi=0}^{\ell} Y \frac{d'(\phi/\ell_c)}{\ell_c} \left(1 - \frac{\phi}{\rho(s)}\right) d\phi}{\int_{\phi=0}^{\ell} \frac{d'(\phi/\ell_c)}{\ell_c} \left(1 - \frac{\phi}{\rho(s)}\right) d\phi} \quad \text{and} \quad \bar{h}(d)Y_c = Y_c \frac{\int_{\phi=0}^{\ell} h(d) \frac{d'(\phi/\ell_c)}{\ell_c} \left(1 - \frac{\phi}{\rho(s)}\right) d\phi}{\int_{\phi=0}^{\ell} \frac{d'(\phi/\ell_c)}{\ell_c} \left(1 - \frac{\phi}{\rho(s)}\right) d\phi} \quad (67) $$

where $(s, \phi)$ is a local basis, $s$ a curvilinear abscissa of the front and $\rho(s)$ the curvature of the damage front. Let $f_n$ denote the evolution criterion (18) at a given time step $n$

$$ f_n = \mu_n^2 \bar{Y}_n - \bar{h}(d) Y_c \quad (68) $$

where $Y_n$ is the local release rate at the same iteration for a unitary load and $\mu_n$ is the load factor. The load factor at step $n$ is computed in order to have

$$ \max_{\Gamma}(f_n) = 0 \quad (69) $$

where $\Gamma$ is the damage front. Thus

$$ \mu_n = \sqrt{\frac{\overline{h(d)Y_c}}{\max_{\Gamma}(\overline{Y}_n)}} \quad (70) $$

According to (18), only a single point of the front would evolve. However, as moving a single point of the front is not possible for a time-discretised problem [13], it is chosen to predict $f_{n+1}$ and to apply a discretised version of the evolution law (18):

$$ \forall \text{ step } n, \quad \left\{ \begin{array}{l} f_{n+1}^{\text{pred}} \le 0 \\ \Delta\phi \ge 0 \\ f_{n+1}^{\text{pred}} \Delta\phi = 0 \end{array} \right. \quad (71) $$

and we impose

$$ \max \Delta\phi = \Delta\phi_{\max} \quad (72) $$

where $\Delta\phi_{\max}$ is chosen to be about 0.25 of the mesh characteristic length and

$$ f_{n+1}^{\text{pred}} = \overline{Y}_{n+1}^{\text{pred}} - \overline{h(d)Y_c}^{\text{pred}} \quad (73) $$
$$ \lambda_1 = \frac{2\sigma_f l_c}{Ew_1}, \quad \lambda_2 = \frac{2\sigma_2 l_c}{Ew_f} = \frac{2\sigma_k l_c}{E(w_f - w_k)} \qquad (66) $$

The first one, $\lambda_1$, characterizes the first slope of the cohesive law and $\lambda_2$ the second one. Let us consider $\lambda_1 > \frac{1}{2}$. Then, as seen above, there exists a critical value of $\hat{\phi}$ that must not be exceeded. This value corresponds to a point of the cohesive law that must not be reached. If the position of the slope-discontinuity point is located before the critical point, the condition $h' \ge 0$ is respected for the first part of the cohesive law. Finally, choosing $\lambda_2 \le \frac{1}{2}$ ensures the condition for the whole local behavior. Some examples of equivalent one-dimensional strain-stress curves for a bi-linear cohesive law are shown in Fig. 8b.

## Appendix 3: Damage evolution in TLS and load factor

We rewrite the evolution criterion (18) by defining the averaged quantities

$$ \bar{Y} = \frac{\int_{\phi=0}^{\ell} Y \frac{d'(\phi/\ell_c)}{\ell_c} \left(1 - \frac{\phi}{\rho(s)}\right) d\phi}{\int_{\phi=0}^{\ell} \frac{d'(\phi/\ell_c)}{\ell_c} \left(1 - \frac{\phi}{\rho(s)}\right) d\phi} \quad \text{and} \quad \bar{h}(d)Y_c = Y_c \frac{\int_{\phi=0}^{\ell} h(d) \frac{d'(\phi/\ell_c)}{\ell_c} \left(1 - \frac{\phi}{\rho(s)}\right) d\phi}{\int_{\phi=0}^{\ell} \frac{d'(\phi/\ell_c)}{\ell_c} \left(1 - \frac{\phi}{\rho(s)}\right) d\phi} \qquad (67) $$

where $(s, \phi)$ is a local basis, $s$ a curvilinear abscissa of the front and $\rho(s)$ the curvature of the damage front. Let $f_n$ denote the evolution criterion (18) at a given time step $n$

$$ f_n = \mu_n^2 \bar{Y}_n - \bar{h}(d) Y_c \qquad (68) $$

where $Y_n$ is the local release rate at the same iteration for a unitary load and $\mu_n$ is the load factor. The load factor at step $n$ is computed in order to have

$$ \max_{\Gamma}(f_n) = 0 \qquad (69) $$

where $\Gamma$ is the damage front. Thus

$$ \mu_n = \sqrt{\frac{\overline{h(d)Y_c}}{\max_{\Gamma}(\bar{Y}_n)}} \qquad (70) $$

According to (18), only a single point of the front would evolve. However, as moving a single point of the front is not possible for a time-discretised problem [13], it is chosen to predict $f_{n+1}$ and to apply a discretised version of the evolution law (18):

$$ \forall \text{ step } n, \quad \left\{ \begin{array}{l} f_{n+1}^{\text{pred}} \le 0 \\ \Delta\phi \ge 0 \\ f_{n+1}^{\text{pred}} \Delta\phi = 0 \end{array} \right. \qquad (71) $$

and we impose

$$ \max \Delta\phi = \Delta\phi_{\max} \qquad (72) $$

where $\Delta\phi_{\max}$ is chosen to be about 0.25 of the mesh characteristic length and

$$ f_{n+1}^{\text{pred}} = \bar{Y}_{n+1}^{\text{pred}} - \overline{h(d)Y_c}^{\text{pred}} \qquad (73) $$
In order to approximate $f_{n+1}^{\text{pred}}$, a first order development and two hypothesis are made:

* strain derivative following the front is only influenced by the variation of the load factor;

* the boundary variations are neglected, because of the numerical difficulties linked to their computation, so $\int_{a_{n+1}}^{b_{n+1}} X_{n+1} d\phi \approx \int_{a_n}^{b_n} (X_n + \Delta X_n) d\phi$.

Thus

$$f_{n+1}^{\text{pred}} \approx f_n + \alpha_n \Delta \mu_n^{\text{pred}} - \beta_n \Delta \phi_n \quad (74)$$

with

$$\left\{ \begin{aligned} \alpha_n &= 2\mu_n \frac{\int_{\phi=0}^{\ell} \frac{1}{2} \varepsilon_n : \mathcal{C} : \varepsilon_n \frac{d'(\phi/\ell_c)}{\ell_c} (1-\frac{\phi}{\rho(s)}) d\phi}{\int_{\phi=0}^{\ell} \frac{d'(\phi/\ell_c)}{\ell_c} (1-\frac{\phi}{\rho(s)}) d\phi} \\ \beta_n &= Y_c \frac{\int_{\phi=0}^{\ell} h'(d) \frac{d'(\phi/\ell_c)}{\ell_c} (1-\frac{\phi}{\rho(s)}) d\phi}{\int_{\phi=0}^{\ell} \frac{d'(\phi/\ell_c)}{\ell_c} (1-\frac{\phi}{\rho(s)}) d\phi} \end{aligned} \right. \quad (75)$$

that are computed all over $\Omega$ with the same variational formulation than $\bar{Y}$. Among all the quantities in (74), only $\Delta\mu$ has a unique value over the whole domain. So, as all $f_{n+1}$ are wanted to be negative and $0 < \Delta\phi_n < \Delta\phi_{\text{max}}$, it can be determined that

$$\Delta \mu_n^{\text{pred}} < \min_{\Gamma} \frac{f_n + \Delta \phi_{\text{max}} \beta_n}{\alpha_n} \quad (76)$$

This quantity is well defined since $\forall s$, $\alpha_n > 0$ which is true for $\mu_n > 0$. Equation (72) imposes

$$\Delta \mu_n^{\text{pred}} = \min_{\Gamma} \frac{f_n + \Delta \phi_{\text{max}} \beta_n}{\alpha_n} \quad (77)$$

As it is possible to have $\mu$ determined by a point of the front that is not advancing, it is chosen to replace its value by the point that has determined $\Delta\mu_n^{\text{pred}}$, that is

$$\mu^{\text{corr}} = \arg \min_{\Gamma} \frac{f_n + \Delta \phi_{\text{max}} \beta_n}{\alpha_n} \quad (78)$$

Finally, the front advance is chosen to be

$$\Delta\phi = \max \left( 0, \frac{f_n + \alpha_n^{\text{corr}} \Delta\mu_n^{\text{pred}}}{\beta_n} \right) \quad (79)$$

in order to have $f_{n+1}^{\text{pred}}$ as close to 0 as possible while verifying (71). $\alpha_n^{\text{corr}}$ denotes $\alpha$ value computed with $\mu^{\text{corr}}$ instead of $\mu$. Let us remark that $\Delta\mu_n^{\text{pred}}$ is just a “prediction” and is different from $\mu_{n+1} - \mu_n$. Moreover, its value is always positive even if the load factor is decreasing over steps. Once $\Delta\phi$ is known over the front, the level-set is updated [13].
In order to approximate $f_{n+1}^{\text{pred}}$, a first order development and two hypothesis are made:

* strain derivative following the front is only influenced by the variation of the load factor;

* the boundary variations are neglected, because of the numerical difficulties linked to their computation, so $\int_{a_{n+1}}^{b_{n+1}} X_{n+1} d\phi \approx \int_{a_n}^{b_n} (X_n + \Delta X_n) d\phi$.

Thus

$$f_{n+1}^{\text{pred}} \approx f_n + \alpha_n \Delta \mu_n^{\text{pred}} - \beta_n \Delta \phi_n \quad (74)$$

with

$$\left\{ \begin{aligned} \alpha_n &= 2\mu_n \frac{\int_{\phi=0}^{\ell} \frac{1}{2} \varepsilon_n \varepsilon_n : C : \varepsilon_n \frac{d'(\phi/\ell_c)}{\ell_c} (1-\frac{\phi}{\rho(s)}) d\phi}{\int_{\phi=0}^{\ell} \frac{d'(\phi/\ell_c)}{\ell_c} (1-\frac{\phi}{\rho(s)}) d\phi} \\ \beta_n &= Y_c \frac{\int_{\phi=0}^{\ell} h'(d) \frac{d'(\phi/\ell_c)}{\ell_c} (1-\frac{\phi}{\rho(s)}) d\phi}{\int_{\phi=0}^{\ell} \frac{d'(\phi/\ell_c)}{\ell_c} (1-\frac{\phi}{\rho(s)}) d\phi} \end{aligned} \right. \quad (75)$$

that are computed all over $\Omega$ with the same variational formulation than $\bar{Y}$. Among all the quantities in (74), only $\Delta\mu$ has a unique value over the whole domain. So, as all $f_{n+1}$ are wanted to be negative and $0 < \Delta\phi_n < \Delta\phi_{\text{max}}$, it can be determined that

$$\Delta \mu_n^{\text{pred}} < \min_{\Gamma} \frac{f_n + \Delta \phi_{\text{max}} \beta_n}{\alpha_n} \quad (76)$$

This quantity is well defined since $\forall s$, $\alpha_n > 0$ which is true for $\mu_n > 0$. Equation (72) imposes

$$\Delta \mu_n^{\text{pred}} = \min_{\Gamma} \frac{f_n + \Delta \phi_{\text{max}} \beta_n}{\alpha_n} \quad (77)$$

As it is possible to have $\mu$ determined by a point of the front that is not advancing, it is chosen to replace its value by the point that has determined $\Delta\mu_n^{\text{pred}}$, that is

$$\mu^{\text{corr}} = \arg \min_{\Gamma} \frac{f_n + \Delta \phi_{\text{max}} \beta_n}{\alpha_n} \quad (78)$$

Finally, the front advance is chosen to be

$$\Delta\phi = \max \left( 0, \frac{f_n + \alpha_n^{\text{corr}} \Delta\mu_n^{\text{pred}}}{\beta_n} \right) \quad (79)$$

in order to have $f_{n+1}^{\text{pred}}$ as close to 0 as possible while verifying (71). $\alpha_n^{\text{corr}}$ denotes $\alpha$ value computed with $\mu^{\text{corr}}$ instead of $\mu$. Let us remark that $\Delta\mu_n^{\text{pred}}$ is just a “prediction” and is different from $\mu_{n+1} - \mu_n$. Moreover, its value is always positive even if the load factor is decreasing over steps. Once $\Delta\phi$ is known over the front, the level-set is updated [13].
Received: 30 April 2015 Accepted: 23 July 2015

Published online: 31 July 2015

References

1. Bažant ZP, Planas J (1997) Fracture and size effect in concrete and other quasibrittle materials. vol 16. CRC Press, Boca Raton

2. Remmers JJC, de Borst R, Needleman A (2003) A cohesive segments method for the simulation of crack growth. Comp Mech 31:69–77

3. Alfaiate J, Wells GN, Sluys LJ (2002) On the use of embedded discontinuity elements with crack path continuity for mode-I and mixed-mode fracture. Eng Fract Mech 69(6):661–686

4. Bažant ZP, Oh BH (1983) Crack band theory for fracture of concrete. Matériaux et Constr 16(3):155–177

5. Moës N, Stolz C, Bernard P-E, Chevaugeon N (2011) A level set based model for damage growth: the thick level set approach. Int J Num Methods Eng 86:358–380

6. Lorentz E, Cuvilliez S, Kazymyrenko K (2011) Convergence of a gradient damage model toward a cohesive zone model. Comptes Rendus Mécanique 339(1):20–26

7. Van den Bosch MJ, Schreurs PJG, Geers MGD (2006) An improved description of the exponential Xu and Needleman cohesive zone law for mixed-mode decohesion. Eng Fract Mech 73(9):1220–1234

8. Dugdale DS (1960) Yielding of steel sheets containing slits. J Mech Phys Solids 8(2):100–104

9. Barenblatt GI (1961) The mathematical theory of equilibrium cracks formed in brittle fracture. Zhurnal Prikladnoy Mekhaniki i Tekhnicheskoy Fiziki 4:3–56

10. Erdogan F, Sih GC (1963) On the crack extension in plates under plane loading and transverse shear. J Basic Eng 85(4):519–527

11. Moës N, Stolz C, Chevaugeon N (2014) Coupling local and non-local damage evolutions with the thick level set model. Adv Model Simul Eng Sci 1(1):1–21

12. van der Meer FP, Sluys LJ (2015) The thick level set method: sliding deformations and damage initiation. Comp Methods Appl Mech Eng 285:64–82

13. Bernard P, Moës N, Chevaugeon N (2012) Damage growth modeling using the thick level set (TLS) approach: efficient discretization for quasi-static loadings. Comp Methods Appl Mech Eng 233–236:11–27

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

15. Stolz C, Moës N (2012) On the rate boundary value problem for damage modelization by thick level-set. ACOMME-2012, pp 205–220

16. Stolz C (2010) Thermodynamical description of running discontinuities: application to friction and wear. Entropy 12(6):1418–1439

17. Cazes F, Coret M, Combescure A, Gravouil A (2009) A thermodynamic method for the construction of a cohesive law from a nonlocal damage model. Int J Solids Struct 46:1476–1490

18. Mazars J, Pijaudier-Cabot G (1996) From damage to fracture mechanics and conversely: a combined approach. Int J Solids Struct 33(20):3327–3342

19. Mazars J (1984) Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture du béton de structure. PhD thesis

20. Irwin GR (1958) Fracture. Encyclopaedia of Physics, Vol. VI. Springer-Verlag, New York, p 168

21. Hillerborg A, Modéer M, Petersson P-E (1976) Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 6(6):773–781

22. Planas J, Elices M (1992) Shrinkage eigen stresses and structural size effect. In: Bazant ZP (ed) Fracture Mechanics of Concrete Structures. Elsevier, London, pp 939–950

23. Smith E (1975) The size of the cohesive zone at a crack tip. Eng Fract Mech 7(2):285–289

24. Grassl P, Grégoire D, Rojas-Solano B, Laura, Pijaudier-Cabot G (2012) Meso-scale modelling of the size effect on the fracture process zone of concrete. Int J Solids Struct 49:1818–1827

25. Alam SY, Saliba J, Loukili A (2013) Study of evolution of fracture process zone in concrete by simultaneous application of digital image correlation and acoustic emission. In: VIII International Conference on Fracture Mechanics of Concrete and Concrete Structures, pp 1–9

26. Haidar K, Pijaudier-Cabot G, Dubé JF, Loukili A (2005) Correlation between the internal length, the fracture process zone and size effect in model materials. Mater Struct 38:201–210

27. Petersson PE (1981) Crack growth and development of fracture zones in plain concrete and similar materials

28. Roelfstra PE, Wittmann FH (1986) Numerical method to link strain softening with failure of concrete. Fract Tough Fract Energy Concre 163–175

29. Park K, Paulino GH, Roesler JR (2008) Determination of the kink point in the bilinear softening model for concrete. Eng Fract Mech 75(13):3806–3818

30. Hoover CG, Bažant ZP (2014) Cohesive crack, size effect, crack band and work-of-fracture models compared to comprehensive concrete fracture tests. Int J Fract 187(1):133–143

31. Gustafsson PJ, Hillerborg A (1985) Improvements in concrete design achieved through the application of fracture mechanics. In: Application of fracture mechanics to cementitious composites, Springer, pp 667–680

32. Guinea GV, Planas J, Elices M (1994) A general bilinear fit for the softening curve of concrete. Mat Struct 27(2):99–105

33. Moës N, Dolbow J, Belytschko T (1999) A finite element method for crack growth without remeshing. Int J Num Methods Eng 46:131–150

34. Gálvez JC, Elices M, Guinea GV, Planas J (1998) Mixed mode fracture of concrete under proportional and nonproportional loading. Int J Fract 94(3):267

35. Cendón DA, Gálvez JC, Elices M, Planas J (2000) Modelling the fracture of concrete under mixed loading. Int J Fract 103(3):293–310

36. Cazes F, Moës N (2015) Comparison of a phase-field model and of a thick level set model for brittle and quasi-brittle fracture. Int J Num Methods Eng
Received: 30 April 2015 Accepted: 23 July 2015

Published online: 31 July 2015

## References

1. Bažant ZP, Planas J (1997) Fracture and size effect in concrete and other quasibrittle materials. vol 16. CRC Press, Boca Raton

2. Remmers JJC, de Borst R, Needleman A (2003) A cohesive segments method for the simulation of crack growth. Comp Mech 31:69–77

3. Alfaiate J, Wells GN, Sluys LJ (2002) On the use of embedded discontinuity elements with crack path continuity for mode-I and mixed-mode fracture. Eng Fract Mech 69(6):661–686

4. Bažant ZP, Oh BH (1983) Crack band theory for fracture of concrete. Matériaux et Constr 16(3):155–177

5. Moës N, Stolz C, Bernard P-E, Chevaugeon N (2011) A level set based model for damage growth: the thick level set approach. Int J Num Methods Eng 86:358–380

6. Lorentz E, Cuvilliez S, Kazymyrenko K (2011) Convergence of a gradient damage model toward a cohesive zone model. Comptes Rendus Mécanique 339(1):20–26

7. Van den Bosch MJ, Schreurs PJG, Geers MGD (2006) An improved description of the exponential Xu and Needleman cohesive zone law for mixed-mode decohesion. Eng Fract Mech 73(9):1220–1234

8. Dugdale DS (1960) Yielding of steel sheets containing slits. J Mech Phys Solids 8(2):100–104

9. Barenblatt GI (1961) The mathematical theory of equilibrium cracks formed in brittle fracture. Zhurnal Prikladnoy Mekhaniki i Tekhnicheskoy Fiziki 4:3–56

10. Erdogan F, Sih GC (1963) On the crack extension in plates under plane loading and transverse shear. J Basic Eng 85(4):519–527

11. Moës N, Stolz C, Chevaugeon N (2014) Coupling local and non-local damage evolutions with the thick level set model. Adv Model Simul Eng Sci 1(1):1–21

12. van der Meer FP, Sluys LJ (2015) The thick level set method: sliding deformations and damage initiation. Comp Methods Appl Mech Eng 285:64–82

13. Bernard P, Moës N, Chevaugeon N (2012) Damage growth modeling using the thick level set (TLS) approach: efficient discretization for quasi-static loadings. Comp Methods Appl Mech Eng 233–236:11–27

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

15. Stolz C, Moës N (2012) On the rate boundary value problem for damage modelization by thick level-set. ACOMME-2012, pp 205–220

16. Stolz C (2010) Thermodynamical description of running discontinuities: application to friction and wear. Entropy 12(6):1418–1439

17. Cazes F, Coret M, Combescure A, Gravouil A (2009) A thermodynamic method for the construction of a cohesive law from a nonlocal damage model. Int J Solids Struct 46:1476–1490

18. Mazars J, Pijaudier-Cabot G (1996) From damage to fracture mechanics and conversely: a combined approach. Int J Solids Struct 33(20):3327–3342

19. Mazars J (1984) Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture du béton de structure. PhD thesis

20. Irwin GR (1958) Fracture. Encyclopaedia of Physics, Vol. VI. Springer-Verlag, New York, p 168

21. Hillerborg A, Modéer M, Petersson P-E (1976) Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 6(6):773–781

22. Planas J, Elices M (1992) Shrinkage eigen stresses and structural size effect. In: Bazant ZP (ed) Fracture Mechanics of Concrete Structures. Elsevier, London, pp 939–950

23. Smith E (1975) The size of the cohesive zone at a crack tip. Eng Fract Mech 7(2):285–289

24. Grassl P, Grégoire D, Rojas-Solano B, Laura, Pijaudier-Cabot G (2012) Meso-scale modelling of the size effect on the fracture process zone of concrete. Int J Solids Struct 49:1818–1827

25. Alam SY, Saliba J, Loukili A (2013) Study of evolution of fracture process zone in concrete by simultaneous application of digital image correlation and acoustic emission. In: VIII International Conference on Fracture Mechanics of Concrete and Concrete Structures, pp 1–9

26. Haidar K, Pijaudier-Cabot G, Dubé JF, Loukili A (2005) Correlation between the internal length, the fracture process zone and size effect in model materials. Mater Struct 38:201–210

27. Petersson PE (1981) Crack growth and development of fracture zones in plain concrete and similar materials

28. Roelfstra PE, Wittmann FH (1986) Numerical method to link strain softening with failure of concrete. Fract Tough Fract Energy Concre 163–175

29. Park K, Paulino GH, Roesler JR (2008) Determination of the kink point in the bilinear softening model for concrete. Eng Fract Mech 75(13):3806–3818

30. Hoover CG, Bažant ZP (2014) Cohesive crack, size effect, crack band and work-of-fracture models compared to comprehensive concrete fracture tests. Int J Fract 187(1):133–143

31. Gustafsson PJ, Hillerborg A (1985) Improvements in concrete design achieved through the application of fracture mechanics. In: Application of fracture mechanics to cementitious composites, Springer, pp 667–680

32. Guinea GV, Planas J, Elices M (1994) A general bilinear fit for the softening curve of concrete. Mat Struct 27(2):99–105

33. Moës N, Dolbow J, Belytschko T (1999) A finite element method for crack growth without remeshing. Int J Num Methods Eng 46:131–150

34. Gálvez JC, Elices M, Guinea GV, Planas J (1998) Mixed mode fracture of concrete under proportional and nonproportional loading. Int J Fract 94(3):267

35. Cendón DA, Gálvez JC, Elices M, Planas J (2000) Modelling the fracture of concrete under mixed loading. Int J Fract 103(3):293–310

36. Cazes F, Moës N (2015) Comparison of a phase-field model and of a thick level set model for brittle and quasi-brittle fracture. Int J Num Methods Eng