# A damage criterion based on energy balance for isotropic cohesive zone model

André Chrysochoos, Loïc Daridon, Mathieu Renouf

► To cite this version:

André Chrysochoos, Loïc Daridon, Mathieu Renouf. A damage criterion based on energy balance for isotropic cohesive zone model. Journal of Theoretical, Computational and Applied Mechanics, 2022, Volume 2, pp.1-20. 10.46298/jtcam.7056. hal-03098095v5

HAL Id: hal-03098095

https://hal.science/hal-03098095v5

Submitted on 11 Mar 2022

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

DOI 10.46298/jtcam.7056

OAI hal-03098095v5

**History**

Received Jan 7, 2021

Accepted Nov 10, 2021

Published Mar 7, 2022

**Associate Editor**

Laurence BRASSART

**Reviewers**

Anonymous

Jean-Jacques MARIGO

Anonymous

**Open Review**

OAI hal-03482933

**Supplementary Material**

Information on last page

**Licence**

CC BY 4.0

©The Authors

# A damage criterion based on energy balance for an isotropic cohesive zone model

André CHRYSOCHOOS<fnref n="1" />,<fnref n="2" />, Loic DARIDON<fnref n="1" />,<fnref n="2" />, and Mathieu RENOUF<fnref n="1" />,<fnref n="2" />

The objective of this paper is to present an energy damage criterion for cohesive zone models within the framework of the non-linear thermodynamics of irreversible processes. An isotropic elastic damageable material is considered for isothermal transformations. Damage is then the only irreversible effect accompanying the deformation process and this mechanism is assumed to be fully dissipative. Once a separation law and a damage state variable have been chosen, it is demonstrated that the damage evolution law can be automatically derived from the energy balance. From this observation, a cohesive zone model is derived for a given choice of traction-separation law and damage state variable and the quality of its numerical predictions is analyzed using an experimental benchmark bending test. Damage, elastic and dissipated energy fields around the crack path are shown during this rupture test. Finally, a numerical simulation of a Brazilian test is proposed where no pre-crack is present in the specimen. Then, as before, the evolution of the dissipated energy fields are plotted during the loading until the total failure of the specimen. In the case of the bending test, the crack path obtained with our model is in the experimental envelope until the end of the simulation. Concerning the force vs. CMOD curve, it fits perfectly in the section corresponding to the linear increase and estimates very well the maximum value and the initial decay phase. In the case of the Brazilian test, with a heterogeneous material, we obtain a macro-crack zigzagging through the all microstructure which starts near the loading plates, as observed in experimental test. In this numerical simulation, a contact law is considered between the lips of the crack when the interface represented by our model is completely ruptured.

**Keywords:** cohesive zone, damage, fracture, thermodynamics of irreversible processes, energy balance, finite element analysis, Brazilian test

## 1 Introduction

In many industrial situations, the management of damage and failure of mechanical structures is crucial. This is the reason why many academic and industrial laboratories have intensively studied and still study this problem from an experimental, theoretical, and numerical point of view. Behavioral models taking into account the damage, cracking and failure of structures have followed roughly two distinct paths. Since the pioneering work of Dugdale (1960) and Barrenblatt (1962), surface approaches were proposed. Their objective was to describe in a practical way the material behavior during its rupture, more precisely during the onset and the propagation of crack. The concept of traction-separation curve associated with the crack tip was introduced to depict the gradual separation of material elements. This type of approach has led to the so-called cohesive zone models. The other path can be characterized by continuum damage approaches that have gradually developed since the works of Rabotnov (1969) and Kachanov (1986). Volume descriptions have often used a scalar or tensorial damage variable to describe the progressive degradation of the material. These variables are still often linked to the loss of material elastic stiffness (Lemaitre 1996).

During these last thirty years, whether mechanical approaches are surface or volume, some have been progressively presented within a thermodynamics of irreversible processes (TIP) framework (Costanzo and Allen 1995; Chandrakanth and Pandey 1995). Thermodynamics

<footnote n="1">¹ LMGC, Université de Montpellier, CNRS, Montpellier, France</footnote>

<footnote n="2">² MIST, Université de Montpellier, IRSN, CNRS, France</footnote>
provides indeed a consistent framework both for discussing the admissibility of the constitutive equations that account for the irreversibility of damage mechanisms, and for characterizing the energy properties of damage and crack growth phenomena. A significant number of works dealing with volume approaches used the TIP with internal variables including state variables related to the damage. The behavioral constitutive equations are then divided into two groups: the state equations derived from the thermodynamic potential, characterizing the properties of equilibrium states of the material, and the complementary (or evolution) equations derived either from threshold functions generally defined in the space of thermodynamic forces associated with the model, or from a dissipation potential in the framework of Generalized Standard Materials (Kondo et al. 2007). Particular attention has been paid to the form of the evolution equations such that the model predictions conform to the second principle of thermodynamics.

For cohesive zone models, the introduction of thermodynamics has been much the same as for the bulk/continuum approaches. As noted by Costanzo and Allen (1995), Gurtin (1979) was probably the first to propose a thermodynamic framework for cohesive zones in fracture. He proposed to consider the crack surface as a two-dimensional thermodynamic system endowed with a potential, (e.g., free energy), dependent on the crack temperature and the crack tip opening displacement. Regarding the evolution laws, and particularly those related to the damage variables, some, like Costanzo and Allen (1995) advocate the use of a dissipation potential in the GSM framework. Others, probably more numerous, proposed forms of differential (kinetic) equations. Although it is impossible to mention all the works related to how damage kinetics were constructed, some references spanning the last twenty years include (Ortiz and Pandolfi 1999; Roe and Siegmund 2003; Evangelista et al. 2013; Serpieri et al. 2015; Kuna and Roth 2015), and more recently (Shu and Stanciuulescu 2020). Here again, when the evolution law is not derived from a dissipation potential, it is necessary to check that the irreversible evolution of the system is in accordance with the second principle of thermodynamics, often formulated via the Clausius-Duhem inequality.

In this paper, we focus on the formulation of a CZ model for an isotropic elastic damageable material. Its main objective is to show that when the damage is the only irreversible mechanism and when this latter is fully dissipative, it is no longer necessary to formulate any hypothesis concerning the damage evolution law, this one being entirely fixed by the energy balance. The potential interest of such an observation is that from now on, experimental techniques dealing with thermal and kinematic full-field measurements allow one to evaluate local energy balances whose results will help to identify the CZ model (Richefeu et al. 2012).

The following sections are devoted to the construction of an energy damage criterion derived from the energy balance for an isotropic elastic damageable material within the TIP framework. Contrary to what is classically done, the damage kinetics is not derived, in this work, from an ad-hoc threshold criterion or a dissipation potential but is based on the premise that the damage progress is linked to a prescribed evolution in the maximum elastic energy that can be stored within the material for a given damage state. From an energy standpoint, it must be noted that the damage mechanisms are considered as the only microstructural irreversible effects accompanying the deformation process and these mechanisms are fully dissipative (no energy storage induced by the material degradation). These restrictive hypotheses are however often implicitly present in the literature we have previously mentioned. To illustrate this statement, let us consider for example the paper proposed by Bouvard et al. (2009). In this paper, the fact that the damage is the only dissipative mechanism is described by a traction force defined via a state law (i.e. no irreversibility is associated with the displacement jump). The fact that the damage mechanisms are considered as fully dissipative can be established once looking at the dissipation form: the dissipation is the product of the thermodynamic force associated with the damage multiplied by the damage rate. In such a case, no energy storage (or release) accompanying the damage progress should appear in the model simulations. Naturally, damage dissipation may induce self-heating leading to non-isothermal deformation processes that are consequently irreversible due to heat diffusion. However, for sake of simplicity, only isothermal transformations are considered in the following sections, and the chosen state variables are the displacement jump vector **u** and a scalar damage variable, denoted by $u_d$. This damage variable $u_d$ is the maximum of the equivalent opening displacement as used in many papers such as (Ortiz and Pandolfi 1999;
van den Bosch et al. 2006; Park et al. 2009; Daridon et al. 2011; Blal et al. 2011) for example.

The layout of the paper is as follows. The energy criterion of the damageable elastic cohesive zone model is presented in Section 2 through a one-dimensional scenario within the TIP framework. In Section 3, a vectorial extension of the cohesive zone law is proposed for an isotropic damage evolution. In Section 4, the capability of the model is investigated using an experimental benchmark test (i.e. a single-edge notch-bending specimen for fracture toughness testing) (Moës et al. 2011; Wojtacki et al. 2017; Galvez et al. 1996). Mechanical and energy responses are shown and discussed. Several damage, elastic and dissipated energy fields around the fracture paths are plotted during the crack propagation. As already mentioned, the computation of the dissipated energy fields is of special interest since they can be compared with the ones derived from quantitative infrared techniques (Chrysochoos 2012; Benaarbia and Chrysochoos 2017). With this perspective, in the final section, a numerical simulation of a Brazilian disc test is proposed. This type of test is well adapted to infrared imaging since the flat surface of the specimen remains perpendicular to the optical axis of the camera until the crack occurs.

## 2 One-dimensional scenario

The objective of the following section is to briefly review the mechanical concepts classically introduced with cohesive zone models (CZM) in the case of a one-dimensional monotonic traction and to embed them into the TIP framework to derive, through an energy criterion, a damage evolution law.

### 2.1 Mechanical aspects

In the literature, the mechanical response of the cohesive zone is described by the correspondence between the normal traction force $f$ supported by the interface and its normal opening displacement $u$, often called separation, during a monotonic opening. Depending on the chosen form of the traction-separation diagram, the relationships are called bilinear, polynomial or exponential cohesive laws. In Figure 1, a polynomial form is chosen to illustrate the most common features of these curves. We find the cohesive strength $f_0$ corresponding to the maximum

<figure><img src="image_2.png" /><figcaption>Figure 1 Traction-separation diagram. Monotonic envelope (solid line), elastic unload or reload (dashed line). An arbitrary polynomial cohesive law is chosen.</figcaption></figure>

of the traction-separation curve or its associated opening displacement $u_0$, the maximum value of separation $u_c$ corresponding to the crack opening. An energy parameter is also often mentioned (Ortiz and Pandolfi 1999): this is the fracture energy (work of separation)

$$A_c = \int_{0}^{u_c} f(u) du \qquad (1)$$

which is the area below the traction-separation curve.

This traction-separation curve is considered as a threshold over which the damage develops irreversibly. This threshold is an intrinsic characteristic of the cohesive zone behavior. When
unloading is considered, it is supposed to be purely elastic, assuming that the damage progress stops as soon as the loading point is below the threshold curve. For convenience, the elastic unloading paths are often directed towards the origin of the traction-separation diagram, see Figure 1. This implies that the elasticity remains linear and that there is no residual opening at the end of the unloading.

The progress of the damage can be depicted by a continuous decrease of the secant stiffness $K = f/u$ towards zero until rupture at $u_c$. A classical scalar definition of the damage variable is

$$D_k = \frac{K_0 - K}{K_0} \qquad (2)$$

where $K_0$ is the initial stiffness of the cohesive zone. The variable $D_k$ progressively increases from 0 to 1 when the opening displacement increases from 0 to $u_c$ (or from $u_e$ to $u_c$ when a pure elastic domain, $[0, u_e]$, is introduced in the traction-separation curve, see Figure 1).

A second possibility is to consider a normalized deformation energy definition of the damage (Ortiz and Pandolfi 1999):

$$D_A = \frac{A}{A_c} \quad \text{where} \quad A = \int_{0}^{u} f(v) \, dv. \qquad (3)$$

Here again, this last definition slightly changes when an elastic domain limited by the point ($u_e, f_e$) is introduced. In such a case, Equation (3) requires a renormalization:

$$D_A^* = \frac{A^*}{A_c^*} \quad \text{where} \quad A^* = \int_{u_e}^{u} f(v) \, dv \quad \text{and} \quad A_c^* = \int_{u_e}^{u_c} f(v) \, dv. \qquad (4)$$

Then, by construction, $D_A$ and $D_A^*$ belong to $[0, 1]$. In fact, there are many ways to define damage. The damage process being assumed irreversible, the damage variable rate is often chosen to be non-negative whatever the loading history, to depict its monotonic evolution. Damage develops when the mechanical state ($u, f$) corresponds to a point of the cohesive threshold curve. In what follows we have chosen a kinematic definition of the damage variable. Like previously done by numerous authors, e.g. (Serpieri et al. 2015), we have chosen the maximum value of the separation $u_d$ ever reached by the cohesive zone until instant $t$. This damage state variable is then defined at instant $t$ by

$$u_d = \max\{u(\tau), \forall \tau \le t\}. \qquad (5)$$

This variable monotonically increases during the damage progress from 0 to $u_c$ whatever the loading path, see Figure 1.

## 2.2 Energy aspects

Usually during a load cycle, the deformation energy $w_{def}$, which corresponds to the area surrounded by the loading curve, see Equation (8), is transformed into dissipated energy, denoted by $w_d$, and stored energy, denoted by $w_s$, due to the irreversible microstructural transformations accompanying the deformation process. Part of $w_{def}$ can also involve strong thermomechanical coupling energy (heat) $w_{thm}$ (Chrysochoos 2012). An illustrative example of such coupling effects in the mechanical response is the famous thermoelastic damping introduced by Zener (1938). The general form of the energy balance over a loading cycle can then be written as

$$w_{def} = w_d + w_s + w_{thm}. \qquad (6)$$

For any other loading, the elastic energy $w_e$ has to be added so that

$$w_{def} = w_e + w_d + w_s + w_{thm}, \qquad (7)$$

where $w_e$ vanishes, by construction, over a loading cycle. In the present situation, we only consider isothermal transformations with no thermomechanical coupling. Moreover, we assume that damage is a pure dissipative mechanism and that, consequently, no energy storage or
release of stored energy, due to microstructural changes, occurs during the loading. These assumptions imply $w_s = w_{\text{thm}} = 0$. For any kind of separation-controlled loading {$u(\tau), \forall \tau \le t$}, the deformation energy at instant $t$ is here defined by

$$w_{\text{def}}(t) = \int_{0}^{t} f(\tau) \dot{u}(\tau) d\tau. \qquad (8)$$

For monotonic loadings, the mechanical state follows the traction-separation curve. The deformation energy then represents the mechanical energy required to reach the damage state $u_d = u(t)$. This cost in deformation energy can be defined by

$$w_{\text{def}}^{\text{d}}(u_{\text{d}}) = \int_{0}^{u_{\text{d}}} f(v) dv. \qquad (9)$$

Another important mechanical energy term is the elastic energy $w_e(u, u_d)$ in the cohesive zone at a given state of damage $u_d$. It is defined by

$$w_e(u, u_d) = \frac{1}{2} K(u_d) u^2. \qquad (10)$$

Note that this energy is mechanically recoverable during the unloading. This is why it does not appear in the general form of the energy balance in Equation (6) for a complete loading cycle.

As previously done for the deformation energy during monotonic loading, we can define the elastic energy by

$$w_e^{\text{d}}(u_d) = \frac{1}{2} K(u_d) u_d^2 = w_e(u_d, u_d) \qquad (11)$$

which represents the maximum elastic energy mechanically recoverable for a given damage state $u_d$.

As previously assumed (no thermomechanical coupling energy, no energy storage), the difference between $w_{\text{def}}^{\text{d}}(u_d)$ and $w_e^{\text{d}}(u_d)$ is attributed to the energy dissipation accompanying the irreversibility of damage mechanisms. The dissipated energy $w_d^{\text{d}}(u_d)$ is then defined by

$$w_d^{\text{d}}(u_d) = w_{\text{def}}^{\text{d}}(u_d) - w_e^{\text{d}}(u_d), \qquad (12)$$

all quantities being illustrated in Figure 2. Based on the mechanical response chosen in Figure 1,

<figure><img src="image_3.png" /><figcaption>**Figure 2** Energy version of the traction-separation diagram. Monotonic envelope (solid line) and elastic unload (dashed line).</figcaption></figure>

the evolutions of the three different energies associated with a loading-unloading tensile test are pictured in Figure 3 showing the deformation, elastic, and dissipated energies. The deformation energy $w_{\text{def}}$ is naturally the sum of the dissipated and elastic energies, respectively $w_d$ and $w_e$ in Equation (7), since the damage is supposed to be the only microstructural transformation which is fully dissipative during loading, see Equation (12) (no energy storage is induced by the microstructural transformations). Figure 3(a) illustrates that during the elastic unloading, $w_e^{\text{d}}$
<figure><img src="image_4.png" /><figcaption>**Figure 3** Energy balance evolution during a load-unload-reload process. Solid lines are associated with the monotonic envelope, dashed lines correspond to the elastic unload and reload. (a) Loading up to $u = u_d$ and unloading. (b) Reloading until rupture.</figcaption></figure>

remains constant (no evolution of damage) while $w_e$ returns to zero. In parallel, the deformation energy $w_{def}$ also decreases and tends towards the energy previously dissipated during the first loading cycle, $w_d^d$. In Figure 3(b), the elastic reloading while $u \le u_d$ is shown and extended by a monotonic loading until rupture for $u_d = u_c$.

Under these restrictive assumptions, the area under the traction-separation curve, Figure 1 (equivalent to a monotonic traction rupture) is completely dissipated when the cohesive zone vanishes. In the next sub-section once the thermodynamic working framework has been specified, this important property is discussed. Then, another point to underline is that if the traction-separation curve is classically considered as the constituent element of the behavior of the cohesion zone, it is thus the same for the evolutions of $w_{def}^d(u_d)$ and $w_e^d(u_d)$. Therefore, instead of using the tension-separation curve to describe the damage progress, associated with the loss of stiffness, it is also possible to use the evolution of the allowable maximum elastic energy $w_e^d$ to define the threshold function associated with the damage rate.

## 2.3 Thermodynamics aspects

In this sub-section, the above results and comments are integrated into the TIP framework.

### 2.3.1 Cohesive zone potential and state laws

In the case of isothermal transformations, the chosen state variables are $(u, u_d)$. The thermomechanical approach starts with the assumption of the existence of a potential $\psi(u, u_d)$ capable of gathering all the state laws. Here, we identify this potential to the elastic energy $w_e$ defined in Equation (11):

$$ \psi(u, u_d) = \frac{1}{2} K(u_d) u^2. \qquad (13) $$

The state laws are, by construction, the partial derivatives of the potential with respect to the state variables. We then define the conjugate variable $f^r$, associated with $u$ which represents the reversible part of the traction force, and $A_d$, the conjugate variable associated with $u_d$:

$$ f^r = \frac{\partial \psi}{\partial u} = K(u_d)u \quad \text{and} \quad A_d = \frac{\partial \psi}{\partial u_d} = \frac{1}{2} K'(u_d)u^2. \qquad (14) $$

### 2.3.2 Clausius-Duhem inequality

The irreversibility of the mechanisms accompanying the opening of the cohesive zone is depicted by the Clausius-Duhem inequality which enables the definition of the intrinsic dissipation $w_d^o$. In the present framework, it can be written as

$$ w_d^o = w_{def}^o - \dot{\psi} = f\dot{u} - \frac{\partial \psi}{\partial u}\dot{u} - \frac{\partial \psi}{\partial u_d}\dot{u}_d = f^{ir}\dot{u} + X_d\dot{u}_d \ge 0. \qquad (15) $$
The terms $w_d^o$ and $w_{def}^o$ determine the dissipated and deformation energy rates, respectively. The symbol ●° is introduced to underline that $w_d$ and $w_{def}$ are not a priori state functions and are then path dependent. Equation (15) also introduces the irreversible part of the traction force $f^{ir} = f - f^r$, and the thermodynamic force $X_d$ associated with $\dot{u}_d$. Note that during an irreversible transformation $\dot{u}_d > 0$, we get $X_d = -A_d$. If damage is the only irreversible process, no dissipation has to be associated with $\dot{u}$. In such a case the irreversible traction force vanishes $f^{ir} = 0$. The traction force $f$ can then be directly defined via the state law

$$f = f^r = K(u_d)u. \qquad (16)$$

Moreover, with Equation (14) and Equation (15), the intrinsic dissipation becomes

$$w_d^o = X_d \dot{u}_d = -\frac{1}{2} K'(u_d) u^2 \dot{u}_d \ge 0. \qquad (17)$$

The fact that $\dot{u}_d \ge 0$ implies $K'(u_d) \le 0$ is physically consistent. The irreversible nature of damage leads to a degradation of the secant stiffness.

### 2.3.3 Threshold function and damage evolution law

In the TIP framework the thermodynamic forces are supposed to be function of the state variable rates. In the case of the linear TIP proposed by Onsager (1931), the correspondence between thermodynamic forces and state variable fluxes is linear. The Onsager matrix is supposed to be symmetric positive definite in order to verify the Clausius-Duhem inequality (positive dissipation) whatever the thermodynamic process. Extension to non-linear theory exists as for example the formalism of Generalized Standard Materials (Halphen and Quoc-Son 1975). Based on the hypothesis of normal dissipation, the thermodynamic forces derive from a convex dissipation potential or equivalently, state variables rates derived from a dual dissipation potential, function of the thermodynamic forces. This dissipation potential can also involve the state variables of the model as parameters (Lemaitre 1996). A common approach is then

* to define a threshold function depending on the thermodynamic forces (and possibly state variables),

* to write that irreversibility occurs and develops if the thermodynamic state is on the threshold and remains on it during a time increment.

Note that once the state laws (derived from the thermodynamic potential) and complementary laws (derived from the dissipation potential) have been derived, it is possible to deduce the evolution of the energy balance associated with the transformation.

In what follows, the existence of the threshold function will not be associated with the normal dissipation hypothesis as the energy balance form imposes by construction non-negative dissipation. Indeed, the current elastic domain is characterized by $w_e^d(u_d)$, the maximum elastic energy available for a given damage state which also corresponds to the energy required to further damage the material. Then, the damage energy criterion based on the energy balance is defined by

$$w_e(u, u_d) \le w_e^d(u_d). \qquad (18)$$

The evolution law for $u_d$ is then derived from the fact that for the damage to occur, the maximum elastic energy allowable in the material has to be and remain on the threshold during the loading step, i.e.

$$w_e(u, u_d) = w_e^d(u_d) \qquad (19a)$$

$$\dot{w}_e(u, u_d) = \dot{w}_e^d(u_d). \qquad (19b)$$

The first equality implies $u_d = u$. The second equality leads to a proposal of evolution equation for the damage:

$$\dot{u}_d = \begin{cases} \dot{u} & \text{if } u = u_d \text{ and } \dot{u} \ge 0 \\ 0 & \text{if } u < u_d \text{ or } \dot{u} \le 0, \end{cases} \qquad (20)$$
which is consistent if we remind the definition of the damage state variable Equation (5) and the fact that the damage increases irreversibly, $\dot{u}_d \ge 0$.

To be fully compatible with non-linear TIP framework, the final step is to introduce a threshold function that accounts for the thermodynamic force $X_d$. As previously stated, we consider a derivative form of the energy balance to get this threshold function Equation (19b). Using Equation (14) and Equation (17), we get

$$ \dot{w}_e = -X_d \dot{u}_d + K(u_d) u \dot{u}. \qquad (21) $$

On the threshold, Equation (21) becomes

$$ (w_e^d)' \dot{u}_d = (-X_d + K(u_d)u_d) \dot{u}_d. \qquad (22) $$

Then, a threshold function $F$ involving the thermodynamic force $X_d$ and the state variables can be considered in the form

$$ F(X_d; u, u_d) = K(u_d)u - X_d - (w_e^d)' \le 0. \qquad (23) $$

To be consistent with the incremental form of the energy balance, the equality $F(X_d; u, u_d) = 0$ implies, once again, $u_d = u$ while the consistency condition $dF = 0$ leads to $du = du_d$, or equivalently to Equation (20). To be precise, the full calculation of $dF = 0$ at $u = u_d$ leads to

$$ (K(u_d) + 2K'(u_d)u_d)(du - du_d) = 0, \qquad (24) $$

and then, $du = du_d$, except possibly when $u_d = -K(u_d)/(2K'(u_d))$.

### 2.3.4 Comments on the damage evolution equations

To depict the evolution of damage, in addition to the traction-separation curve data, the literature often proposes a specific evolution equation in the form of $\dot{D} = \dot{D}(f, D, \dot{u})$ whatever the definition of the damage variable $D$ (Roe and Siegmund 2003; Bouvard et al. 2009; Kuna and Roth 2015).

In the foregoing, because of the hypotheses explicitly made on the energy balance (i.e. damage is the only dissipative mechanism and it is totally dissipative), the damage evolution law is fixed by the definition of the damage variable itself and by the explicit form of the energy balance. Note that the damage evolution law (20) deduced from the energy criterion (19) is perfectly compatible with the definition of the damage variable itself given in Equation (5). We can also note that this evolution law is an extremely simple form of the general equation proposed by Roe and Siegmund (2003), but here, this law is totally imposed by the shape of the traction-separation curve, or equivalently, by the threshold $w_e^d(u_d)$ of Equation (11).

To set ideas, let's consider the following simple case: let $f(u)$ be a one-dimensional traction-separation law. We suppose that the elastic energy is, as often, written as $2\psi(u, D) = (1 - D)K_0 u^2$, where $D$ is the isotropic damage variable, $u$, the displacement jump, and $K_0$, the elastic stiffness of the virgin cohesive zone. We consider a monotonic loading. The deformation energy rate is, by definition,

$$ \dot{w}_{\text{def}} = f(u) \dot{u} \qquad (25) $$

where $f(u)$ follows the traction-separation curve. The elastic energy rate can be split into two parts:

$$ \dot{w}_e = (1 - D)K_0 u \dot{u} - K_0 \frac{u^2}{2} \dot{D}. \qquad (26) $$

If we assume now that the damage is the only irreversible mechanism, then the traction force is the conjugate variable of the displacement jump, where

$$ f(u) = \frac{\partial \psi(u, D)}{\partial u} = (1 - D)K_0 u \quad \text{and} \quad \dot{w}_e = \dot{w}_{\text{def}} - K_0 \frac{u^2}{2} \dot{D}. \qquad (27) $$
If the damage is supposed to be exclusively dissipative (no internal stored energy), then the dissipation is given by

$$ \dot{w}_d = K_0 \frac{u^2}{2} \dot{D} = \dot{w}_{\text{def}} - \dot{w}_e. \qquad (28) $$

Following the traction-separation curve, the damage evolution has to satisfy

$$ \dot{D} = 2 \frac{\dot{w}_{\text{def}} - \dot{w}_e}{\frac{1}{2} K_0 u^2}. \qquad (29) $$

Noting that for each current point (u, f(u)) of the traction separation curve, we have $2w_e = f(u)u$, the time differentiation, following the curve, reads

$$ 2\dot{w}_e = f(u)\dot{u} + \dot{f}(u)u \qquad (30) $$

and, as a consequence,

$$ \dot{D} = \frac{f(u)\dot{u} - \dot{f}(u)u}{\frac{1}{2}K_0 u^2}. \qquad (31) $$

The right-hand member of this equation is fully determined by the traction-separation curve. Any form of damage evolution law, incompatible with this previous equation, would lead to an energy balance form incompatible with the initial energy assumptions (i.e. form of the free energy, damage unique and exclusive dissipative mechanism). The consequences could be the emergence of energy storage mechanisms, i.e. $\dot{w}_d^d < \dot{w}_{\text{def}}^d - \dot{w}_e^d$, or internal energy transformation into dissipated energy (release of stored energy), i.e. $\dot{w}_d^d > \dot{w}_{\text{def}}^d - \dot{w}_e^d$. Taking into account this stored energy, variations should lead to the introduction of new internal state variables and/or to a change of the deformation energy rate definition (Frémond 2002).

# 3 Three-dimensional cohesive zone model

In this section, we propose an extension to a three-dimensional vectorial version of the CZM where the isotropic damage is controlled by the evolution of the maximum recoverable elastic energy, $w_d^e(u_d)$. Isotropic damage means here that a scalar state variable is solely used to describe the damage evolution. This generalization has been made by following the same approach as the one previously proposed, namely define a damage variable and a energy balance where the damage is the only dissipative phenomenon.

## 3.1 Mechanical variables

Regarding the mechanical description of the cohesive zone, the traction force and the separation become now vectors. Let us introduce a frame of reference where directions 1 and 2 correspond to the tangent plane of the cohesive zone while direction 3, is the normal direction. The traction vector, f, with components ($f_{t_1}$, $f_{t_2}$, $f_n$) and the separation vector, u, with components ($u_{t_1}$, $u_{t_2}$, $u_n$), are introduced. As conventionally admitted in CZM, the normal move jump denoted by $u_n$ is positive or null. This unilateral condition is taken into account by a Signorini type relationship in the numerical simulations using the open source software IMGc90 (Dubois and Jean 2005) and exposed at the end of this article.

## 3.2 Cohesive zone potential and state equations

A set of state variables has first to be chosen. Here, we selected the components ($u_{t_1}$, $u_{t_2}$, $u_n$) of the separation vector and a scalar damage variable denoted by $u_d$. Then, to generalize the form of the cohesive zone potential proposed in Equation (13), the form

$$ \begin{aligned} 2\psi(\mathbf{u}, u_d) &= 2w_e(\mathbf{u}, u_d) \\ &= K_n(u_d)u_n^2 + K_t(u_d)u_{t_1}^2 + K_t(u_d)u_{t_2}^2 \\ &= K_n(u_d)(u_n^2 + \alpha u_{t_1}^2 + \alpha u_{t_2}^2) \\ &= K_n(u_d)u_{\text{eq}}^2, \quad \text{with } u_{\text{eq}} = \sqrt{u_n^2 + \alpha u_{t_1}^2 + \alpha u_{t_2}^2}. \end{aligned} \qquad (32) $$
inspired by Bouvard et al. (2009) is adopted. The parameter $\alpha$ is the ratio between $K_t(u_d)$ and $K_n(u_d)$, the tangential and normal secant stiffnesses, respectively, at a given $u_d$. In the case of isotropic damage, $\alpha$ is a constant.

A three-dimensional formulation of the scalar depicting the isotropic damage is

$$u_d = \max \{u_{eq}(\tau), \forall \tau \le t\}. \qquad (33)$$

By construction, $u_d$ takes the three-dimensional aspect of the separation vector $\mathbf{u}$ into account and then $\dot{u}_d$ is non-negative and de facto respects the irreversibility of the damage progress.

By definition, the state laws are the partial derivatives of the cohesive zone potential of Equation (32). They introduce the components of the reversible traction vector $\mathbf{f}^r$ and the conjugate variable $A_d$ associated with $(u_{t_1}, u_{t_2}, u_n)$ and $u_d$, respectively:

$$\left\{
\begin{aligned}
f_n^r &= \frac{\partial \psi}{\partial u_n} = K_n(u_d)u_n, \\
f_{t_1}^r &= \frac{\partial \psi}{\partial u_{t_1}} = \alpha K_n(u_d)u_{t_1} = K_t(u_d)u_{t_1}, \\
f_{t_2}^r &= \frac{\partial \psi}{\partial u_{t_2}} = \alpha K_n(u_d)u_{t_2} = K_t(u_d)u_{t_2}, \\
A_d &= \frac{\partial \psi}{\partial u_d} = \frac{1}{2} K'(u_d)u_{eq}^2.
\end{aligned}
\right.
\qquad (34)$$

Because only damage induces irreversibility, no dissipation has to be associated with the component of the separation vector. The reversible part $\mathbf{f}^r$ of the separation can therefore be identified with $\mathbf{f}$, that is $\mathbf{f} = \mathbf{f}^r$.

### 3.2.1 Energy definition of the damage threshold

To extend the damage energy criterion to a three-dimensional isotropic damageable CZM, it is possible to choose the damage variable, $u_d$, whose evolution is directly related to that of the elastic energy $w_e^d$. In the three-dimensional case, this maximum elastic energy, for a given damage state $u_d$, describes, in the displacement space, a half spheroid of radii

$$r_n(u_d) = \left( \frac{2w_e^d(u_d)}{K_n(u_d)} \right)^{\frac{1}{2}} \quad \text{and} \quad r_t(u_d) = \left( \frac{2w_e^d(u_d)}{\alpha K_n(u_d)} \right)^{\frac{1}{2}} = \frac{r_n(u_d)}{\sqrt{\alpha}} \qquad (35)$$

as shown in Figure 4. As the normal jump $u_n$ is, by definition, positive or null, only half of the

<figure><img src="image_5.png" /><figcaption>**Figure 4** Reachable separation states for a given damage state $u_d$. The color represents the value of $u_n$.</figcaption></figure>

spheroid is reachable for any separation states.

As long as the further separation states, $\mathbf{u}$, respect the damage energy criterion (i.e. $w_e(\mathbf{u}, u_d) < w_e^d(u_d)$), the behavior remains elastic. Then, for a given opening such that $u_{eq} = u_d$, the elastic energy reaches the maximal value associated with this damage state (i.e. $w_e(\mathbf{u}, u_d) = w_e^d(u_d)$). Once the surface of the spheroid is reached:

* either the separation increment $\delta\mathbf{u}$ is directed towards the inside of the spheroid, and an elastic unloading at constant damage can be observed,

* or $\delta\mathbf{u}$ is directed towards the outside of the spheroid, and then the damage develops and defines a new elastic limit surface.
For isotropic damage, a single evolution equation for $u_d$ is required. We have already underlined that for threshold behavior law, the yield function depends on the thermodynamic forces and possibly on the states variables themselves, acting as parameters. In the present case, the thermodynamic force of the model, associated with the damage variable rate, is $X_d$. A generalized form of the yield criterion (23) is chosen where the role of $u$ used in the one-dimensional scenario is played by $u_{eq}$. Accordingly, the proposed yield energy criterion (18) may be rewritten using the thermodynamic force $X_d$ as

$$F(X_d; u_{eq}, u_d) = K_n(u_d)u_{eq} - X_d - (w_e^d)' \le 0. \quad (36)$$

Damage develops if the threshold is reached, that is $F(X_d; u_{eq}, u_d) = 0$, and if the consistency condition $\dot{F}(X_d; u_{eq}, u_d) = 0$ is satisfied. For the same reasons as the ones shown for the one-dimensional model, the evolution law of the parameter $u_d$ is expressed as

$$\dot{u}_d = \dot{u}_{eq} \quad \text{if} \quad u_d = u_{eq} \quad \text{and} \quad \dot{u}_{eq} \ge 0. \quad (37)$$

This evolution law, as already underlined, is induced by the very definition of the damage variable.

An illustration of the energy criteria is given in Figure 5. Following a monotonic loading (i.e. remaining on the $w_e^d(u_d)$ curve), $K_n(u_d)u_d$ is the slope of the deformation energy $w_{def}'$, $X_d$ is the slope of the dissipated energy $w_d^d$ and $(w_e^d)'$ is naturally the slope of the maximal allowable elastic energy $w_e^d$.

<figure><img src="image_6.png" /><figcaption>Figure 5 Damage energy criterion $u_{eq} = u_d$ where $K_n(u_d)u_d = (w_e^d)'u_d + X_d$.</figcaption></figure>

To conclude the present section, we would like to stress once again the fact that the damage evolution law is not here a matter of choice. It is imposed by the form of the energy balance and by the definition of the damage state variable.

The state (34) and the evolution (37) are, in what follows, implemented in a home-made finite element code. The material parameters of the constitutive equations will be specified. In order to show the capabilities of such a CZM, two types of simulations are undertaken. The first one is a bending test whose numerical results are compared with experimental ones. The second one is a Brazilian disc test whose material is made of heterogeneous elastic grains.

## 4 Numerical implementation

To illustrate the potential of the proposed model, simulations reproducing a common benchmark extracted from the literature (Galvez et al. 1996) were carried out. It is important to notice that the objective of this practical comparison is simply to demonstrate the use of the model and not to optimize its parameters in order to fit the benchmark. The previous model is then numerically implemented in the code IMC90 based on the Non-Smooth Contact Dynamics (NSCD) method (Moreau 1988; Jean 1999; Jean et al. 2001). The NSCD method is dedicated to solving problems related to dynamic systems with unilateral constraints. It is therefore particularly suitable for unilateral contact problems with friction. It offers a non-smooth treatment (no
compliance, no penalty) of the conditions of contact (Jean 1999), which is explicit in the definition of $u_n$. The way adhesion is taken into account in this method makes it possible to consider each point of contact as a cohesive zone. The mechanical behavior of the cohesive zones may then vary at any point of the spatial discretization of the problem. This relevant modeling framework was adopted to numerically simulate crack propagation with cohesive zone (Champagne et al. 2014).

## 4.1 Bending test

To compare the proposed model with a benchmark solution (Galvez et al. 1996), the form of the maximum storable elastic energy $w_e^d(u_d)$ must be specified in order to be implemented in IMGC90 (Dubois et al. 2022), the open source platform used to carry out the simulations (Dubois et al. 2011). This benchmark, illustrated in Figure 6, traces the evolution of a crack in mixed mode to be followed. In the context of this feasibility study, a simple quadratic form of $w_e^d(u_d)$ is proposed.

<figure><img src="image_7.png" /><figcaption>**Figure 6** Characteristics of the benchmark issued from Galvez et al. (1996) used for simulation</figcaption></figure>

In what follows, we also assumed the existence of a pure elastic domain and thus the existence of a threshold equivalent elastic displacement $u_{eq}^e$, denoted by $u_e$. The maximum storable elastic energy as a function of the damage parameter $u_d$ simply reads

$$w_e^d(u_d) = A(u_d - u_c)^2 + B(u_d - u_c) \quad \text{if } u_e \le u_d \le u_c \qquad (38)$$

where $u_c$ is the critical equivalent displacement corresponding to the crack onset. Parameters A and B are two constants chosen to ensure the $C^1$ continuity of the maximum storable elastic energy $w_e^d(u_d)$ at the threshold equivalent elastic displacement $w_e^d(u_e) = \frac{1}{2}K_n^0 u_e^2$. They are defined by

$$A = -\frac{1}{2}K_n^0 u_e \frac{2u_c - u_e}{(u_c - u_e)^2} \quad \text{and} \quad B = -K_n^0 \frac{u_c u_e}{u_c - u_e} \qquad (39)$$

where $K_n^0$ is the initial normal stiffness of the CZM.

The thickness $H$ of the sample is equal to 0.3 m while its length is equal to 1.2 m. A 0.15 m pre-crack is located in the middle of the lower edge. The point $B$ is fixed in both $x$ and $y$ directions whereas the point $A$ is only fixed in the $y$ direction. A displacement is imposed on point $A$ to load the structure. The mesh is composed of three parts: two coarse meshes, the left and the right parts of the structure composed respectively of 958 and 2063 T3 elements, where no interface elements have been introduced between the different meshes and a finer mesh, guarantying the continuity of the structure (domain $\Omega_1$ in Figure 6), composed of 6723 T3 elements where the crack path is supposed to appear and where interface elements are therefore introduced between each element. The interactions between elements of $\Omega_1$ are governed by the proposed cohesive zone model where the initial secant elastic stiffnesses, $K_n^0$ and $K_t^0$, are chosen to satisfy the criterion proposed in (Blal et al. 2011) to limit the reduction of stiffness due to the presence of CZM. It is important to highlight that the objective of this practical comparison is simply to show the functional features of the model and not to optimize its parameters in order to fit the benchmark. The values of the CZM parameters are listed in Table 1.
<table><caption>**Table 1** Parameter values of the CZM for the bending test.</caption><thead><tr><th>K<sup>0</sup><sub>n</sub> (N m<sup>-1</sup>)</th><th>α</th><th>u<sub>e</sub> (m)</th><th>u<sub>c</sub> (m)</th></tr></thead><tbody><tr><td>2.48 10<sup>9</sup></td><td>0.5</td><td>0.5 10<sup>-6</sup></td><td>1.5 10<sup>-6</sup></td></tr></tbody></table>

Figure 7 shows the evolution, at various simulation times, of different characteristic quantities associated with the model: the damage variable, the elastic energy $w_e$ and the dissipated energy $w_d$. In order to present a quantity varying from 0 to 1, the damage ratio, as a function of the damage variable, is introduced and defined by $\frac{\langle u_d - u_e \rangle^+}{u_c - u_e}$. To improve the visibility of these different quantities supported by the interfaces, they are projected on adjacent elements.

<figure><img src="image_8.png" /><figcaption>**Figure 7** Damage ratio (top), elastic energy (center) and dissipated energy (bottom) during crack propagation. Labels (a), (b) and (c) apply column-wise.</figcaption></figure>

Figure 7(a), corresponding to a pre-cracking state, shows a concentration of the elastic energy at the outset of the crack tip. However, the damage criterion has not been reached within the cohesive zone so that no damage or dissipation has yet occurred, see Equation (18). The maps corresponding to $\frac{\langle u_d - u_e \rangle^+}{u_c - u_e}$ and $w_d$ are then uniformly equal to 0. As expected, Figure 7(b) and Figure 7(c), corresponding to two post-cracking steps, highlight the correlation between the evolution of the dissipated energy and the damage ratio. The elastic energy is still concentrated ahead of the crack tip, then returns to zero along the crack lips. In contrast, the dissipated energy related to the damage evolution can be exhibited all along the crack path. Similarly, the damage field allows the cracking path to be tracked.

To exhibit the capability of our CZM where only the shape of the cohesive energy associated with a simple energy balance is needed, see Equation (38), various quantities, numerically obtained, are compared with experimental measurements from the literature (Cendón et al. 2000). For such comparisons, Figure 8 presents both the classical crack path monitoring and the load vs. CMOD curve (Crack Mouth Opening Displacement).

In Figure 8(a) and Figure 8(b), dots correspond to the simulation while the dashed lines represent the crack envelop obtained experimentally by Cendón et al. (2000). In Figure 8(a), the crack obtained numerically corresponds closely to the experimental envelope. The starting angle is strongly related to the discretization around the initiation point, explaining the slight difference at the beginning of the initiation. Then, the path is corrected and repositioned in the experimental envelope until the end of the simulation.

Concerning the force vs. CMOD curves, they fit perfectly in the section corresponding to the linear increase. This highlights that the introduction of a two-dimensional interface element,
<figure><img src="image_9.png" /><figcaption>**Figure 8** Comparison of numerical macroscopic measurements associated to the crack evolution with experimental results (Cendón et al. 2000). (a) Crack path and (b) load vs. CMOD curve.</figcaption></figure>

where the values of $K_n^0$ and $K_t^0$ satisfy the criterion proposed in (Blal et al. 2011) between each elements of $\Omega_1$ do not affect the global stiffness of the sample. The maximum force obtained is also in good agreement with that obtained in the experiment, as well as the beginning of the non-linear decreasing part of the CMOD curve occurring at the initiation of cracking. In the last part, the curves diverge. This difference is partly explained by the fact that the numerical simulation is two-dimensional while the experiments are three-dimensional. Indeed, not all deformation modes are taken into account (especially out-of-plane modes), which explains this different behavior at the end of the simulation. Moreover, we have arbitrarily chosen a second degree polynomial to characterize the damage of the cohesive zone model, see Equation (38). This choice could be fine-tuned in order to better account for experiences by taking a Needleman-type damage (van den Bosch et al. 2006; Needleman 1990).

## 4.2 Sensitivity analysis

Finally, in order to see the impact of a variation in the parameters $u_e$ and $u_c$ on the overall behavior of the system and more particularly on the evolution of the force vs. CMOD curves, a sensitivity study is proposed. The influence of these parameters on the crack path is not presented because it is not very significant. The influence of these parameters on the energy available to be dissipated in the model is pointed out in Figure 9. The parametric study is carried out relative to the reference point (0, 0) corresponding to the results presented in Figure 8(b) with the parameters defined in Table 1. With the chosen law, a variation of $u_c$ has almost the same consequence as a variation of $u_e$ in terms of energy available to be dissipated. The map in Figure 9 is quasi-symmetric with respect to the circle-triangle diagonal. Concerning the parametric studies, the color code used for the curves is defined in Figure 9.

<figure><img src="image_10.png" /><figcaption>**Figure 9** Map of the normalized dissipated energy variations as a function of the variations of $u_e$ and $u_c$. Symbols at the four corners of the map identify the curves in Figure 10.</figcaption></figure>
Figure 10 presents the normalized plots of $w_e^d(u_d)$ for different values of $u_e$ and $u_c$ and the corresponding force vs. CMOD curves. The normalization parameters are $u_{e,0} = 0.5 \ 10^{-2}$ and $w_{e,0} = w_{e,0}^d(u_{e,0})$ using the values in Table 1. Even if the shape of the curves is significantly different in Figure 10(a), the energies available to be dissipated for the case represented by a cross and the one represented by a square are of same order of magnitude. The maxima order observed at the scale of the CZM models in Figure 10(a) is preserved at the scale of the structure in Figure 10(b).

<figure><img src="image_11.png" /><figcaption>**Figure 10** (a) Normalized $w_e^d(u_d)$ for different values of $u_e$ and $u_c$ with $u_{e,0} = 0.5 \ 10^{-2}$ and $w_{e,0} = w_{e,0}^d(u_{e,0})$ using values in Table 1. (b) Corresponding force vs. CMOD curves.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11** Comparison between the reference load vs. CMOD curve and the ones related to the variation of (a) $u_c$ of ±20% while $u_e = u_{e,0}$ and of (b) $u_e$ of ±20% while $u_c = u_{c,0}$.</figcaption></figure>

In Figure 11 to Figure 12, we observe respectively the influence of $u_c$ and $u_e$ on the force vs. CMOD curves. These figures show that the influence of the variation of $u_c$ is less than that of $u_e$. Indeed, where we observe a variation of less than 10% on the critical values of the curve ($F_{max}$ and CMOD$_{max}$) for a variation of $u_c$, we observe a variation of more than 20% for an equivalent variation of $u_e$. Nevertheless, in both cases, an increase of the damage energy $w_d$ induces an increase of the CMOD and loading maxima in the Load vs. CMOD curve. In this model, where an elastic domain is assumed, $u_e$ is the threshold where the damage begins to occur. This value determines the outset of the non-linear response of the structure. This is exhibited in Figure 11(b) where an increase of $u_e$ at the local scale induces an increase of the maximal force at the macroscopic scale and a delay of the occurrence of the nonlinear response of the curves.

Figure 12 demonstrates that the non-linear region of the curve is also governed by the shape of the energy curve, see Figure 10. Although the dissipative energy in this parametric study is almost constant, we observe a variation of about 10% on the characteristic values of the response curve. Accordingly, by combining the effects of $u_e$, $u_c$ and the shape of the local curve, see Figure 10, it is possible to obtain a better optimal result to fit experiments. The experimental characterization of this type of local curve depicting the micro-structural phenomenon linked to
<figure><img src="image_13.png" /><figcaption>**Figure 12** Comparison between the reference Load vs. CMOD curve (dash line) and the ones related to the variation of $u_c$ and $u_e$ of ±20 % while $w_d^o$ is constant.</figcaption></figure>

fracture is relevant and is still an on-going problem.

## 4.3 Brazilian test

As a complement to the previous numerical simulation and as an opening to the continuation of the present work, the developed CZM law is used in the simulation of a Brazilian test, with CZM parameters listed in Table 2. This test consists of the compression of a circular sample

<table><caption>**Table 2** Parameter values of the CZM for the Brazilian test.</caption><thead><tr><th>K<sub>n</sub><sup>0</sup> (N m<sup>-1</sup>)</th><th>α</th><th>u<sub>e</sub> (m)</th><th>u<sub>c</sub> (m)</th></tr></thead><tbody><tr><td>2.48 10<sup>15</sup></td><td>0.4</td><td>0.8 10<sup>-6</sup></td><td>10<sup>-6</sup></td></tr></tbody></table>

located between two rigid plates. In contrast to the previous case, no pre-crack is introduced in the numerical model. The microstructure of the sample used is presented on the left side of Figure 13. This microstructure has been generated using the open-source Software Neper (Quey

<figure><img src="image_14.png" /><figcaption>**Figure 13** Meshed microstructure used in the simulation of the Brazilian disc test.</figcaption></figure>

et al. 2011). It is composed of 1000 elastic grains following a normal size distribution to make the microstructure heterogeneous (Ma and Huang 2018). In these case, the cohesive zones are only introduced at the grain boundaries. The mesh size used for meshing is identical for all grains and calibrated so that the smallest grains have at least two elements on their smallest side. The total number of elements is 98378. A zoom of the mesh is shown on the right hand side of Figure 13. For the sake of simplicity, the diameter of the sample is unitary. A vertical velocity is imposed on both walls to compress the sample. The simulation is carried out within the framework of large deformations in order to manage possible strain localization and grain rotations. The elastic constitutive equations are those of the linear elasticity where stresses and strains are respectively represented by the second Piola-Kirchhoff stress tensor and the Green-Lagrange strain tensor.

For an homogeneous material during Brazilian test, a tensile state is induced in the center of the disc perpendicular to the load direction. Increasing the load leads to an increase in tensile stress until a crack appears in the center of the disc. Under the effect of the load, the crack develops until the disc eventually separates into, at least, two parts. For a heterogeneous material,
the damage occurs near the rigid plates and then develops along the loading axis until the disk breaks (Na et al. 2017). Figure 14 shows a visual of the sample at the end of the simulation.

<figure><img src="image_15.png" /><figcaption>**Figure 14** Final state of the simulation. (a) Crack through the microstructure. (b) Displacement field norm normalized by the diameter *D* within the sample.</figcaption></figure>

Figure 14(a) shows the macro crack zigzagging through the microstructure. This extends from the contact between the sample and the rigid plates away from the centre of the sample. Figure 14(b) shows the norm of the nondimensionalised displacement field within the sample. The discontinuities highlight the multiple cracking paths generated during compression. Numerous disjointed fragments can be particularly seen in the volume near the top wall.

To complete these observations, the evolution of the dissipated energy fields is presented in Figure 15. Figure 15(a) corresponds to the initiation of the crack while Figure 15(f) corresponds

<figure><img src="image_16.png" /><figcaption>**Figure 15** Evolution of the dissipated energy field from the crack initiation (a) to the end of simulation (f).</figcaption></figure>

to the end of the simulation. The other images are captured at intermediate times. Through Figure 15(a) to Figure 15(f), the damage evolution is exhibited where branching is observed until the coalescence of the macro crack. The next step is to experimentally perform the same type of test using an experimental setup coupling kinematic and thermal full-field measurements. The kinematic measurements will allow us to locate zones of strain localization and even discontinuities of the displacement fields while the thermal measurements will be used to determine the zones where the dissipation is localized. The confrontation of these two pieces of
information should help us check the relevance of this energy approach of cohesive zones.

## 5 Conclusion

In this paper, we present an energy criterion for cohesive zone models where the damage progress is assessed along with the ability of the material to store energy elastically. The damage state variable used is $u_d(t) = \sup\{u_{eq}(\tau), \tau \le t\}$ where $u_{eq}$ is an equivalent displacement compatible with the isotropic evolution of the damage progress. The paper shows that if damage is the unique and exclusive dissipative mechanism, the damage evolution law is automatically set by the evolution of the maximum storable elastic energy $w_e^d(u_d)$. We have also underlined that the data of this energy is equivalent, in a one-dimensional formalism, to the one of a traction-separation law. The interest of this energy approach is its immediate generalization to three-dimensional cohesive zone models. In order to check the operational character of this type of approach, the isotropic damage model has been implemented in the open source software LMGc90 based on Non-Smooth Contact Dynamics and used to perform numerical simulations in the case of bending and Brazilian tests. The results obtained for this plain stress modeling are encouraging. Using a simple quadratic function $w_e^d(u_d)$ for the interface, we obtained a close correlation between the simulations and the experimental observations of the crack path for the bending test along with realistic multicrack propagations in the case of the Brazilian test. A parametric study of the macroscopic response of the structure naturally demonstrates the importance of the shape of the function $w_e^d(u_d)$ which characterizes the interface behavior between two elements. It is indeed this quantity that should be identified experimentally. In subsequent theoretical developments, we will consider an extension to a non-isotropic degradation of the material elastic properties. From an experimental stand point, Brazilian tests will be performed using full field techniques during monotonic loadings, the goal being to extract from the experimental data valuable information on the form of the energy balance and particularly on $w_e^d(u_d)$.

## References

* Barrenblatt, G. I. (1962). The mathematical theory of equilibrium of cracks in brittle fracture. *Advances in Applied Mechanics* 7:55-129. [DOI], [HAL].

* Benaarbia, A. and A. Chrysochoos (2017). Proper orthogonal decomposition preprocessing of infrared images to rapidly assess stress-induced heat source fields. *Quantitative InfraRed Thermography Journal* 14(1):132-152. [DOI], [HAL].

* Blal, N., L. Daridon, Y. Monerie, and S. Pagano (2011). Criteria on the artificial compliance inherent to the intrinsic cohesive zone. *Comptes Rendus Mécanique* 339(12):789-795. [DOI], [HAL].

* van den Bosch, M., P. Schreurs, and M. Geers (2006). An improved description of the exponential Xu and Needleman cohesive zone law for mixed-mode decohesion. *Engineering Fracture Mechanics* 73(9):1220-1234. [DOI].

* Bouvard, J.-L., J. L. Chaboche, F. Feyel, and F. Gallerneau (2009). A cohesive zone model for fatigue and creep-fatigue crack growth in single crystal superalloys. *International Journal of Fatigue* 31(5):868-879. [DOI], [HAL].

* Cendón, D., J. Gálvez, M. Elices, and J. Planas (2000). Modelling the fracture of concrete under mixed loading. *International Journal of Fracture* 103(3):293-310. [DOI], [OA].

* Champagne, M., M. Renouf, and Y. Berthier (2014). Modeling wear for heterogeneous bi-phasic materials using discrete elements approach. *Journal of Tribology* 136(2). [DOI], [HAL].

* Chandrakanth, S. and P. Pandey (1995). An isotropic damage model for ductile material. *Engineering Fracture Mechanics* 50(4):457-465. [DOI].

* Chrysochoos, A. (2012). Infrared thermography applied to the analysis of material behavior: a brief overview. *Quantitative InfraRed Thermography Journal* 9(2):193-208. [DOI], [HAL].

* Costanzo, F. and D. Allen (1995). A continuum thermodynamic analysis of cohesive zone models. *International Journal of Engineering Science* 33(15):2197-2219. [DOI].

* Daridon, L., B. Wattrisse, A. Chrysochoos, and M. Potier-Ferry (2011). Solving fracture problems using an asymptotic numerical method. *Computers & Structures* 89(5):476-484. [DOI], [HAL].
Dubois, F. and M. Jean (2005). The non smooth contact dynamic method: recent LMGc90 software developments and application. *4th Contact Mechanics International Symposium* (Hannover, Germany, July 4, 2005–July 5, 2005). [DOI], [HAL].

Dubois, F., M. Jean, R. Mozul, and M. Renouf (2022). *LMGc90*. Version r2022. [URL].

Dubois, F., M. Jean, M. Renouf, R. Mozul, A. Martin, and M. Bagnéris (2011). *LMGc90: 10e colloque national en calcul des structures* (Giens, France, May 9, 2011–May 13, 2011). [HAL].

Dugdale, D. (1960). Yielding of steel sheets containing slits. *Journal of the Mechanics and Physics of Solids* 8(2):100–104. [DOI].

Evangelista, F., J. Roesler, and S. Proença (2013). Three-dimensional cohesive zone model for fracture of cementitious materials based on the thermodynamics of irreversible processes. *Engineering Fracture Mechanics* 97:261–280. [DOI].

Frémond, M. (2002). *Non-Smooth Thermomechanics*. Springer-Verlag. [DOI].

Galvez, J., M. Elices, G. V. Guinea, and J. Planas (1996). Crack trajectories under mixed mode and non-proportional loading. *International Journal of Fracture* 81(2):171–193. [DOI], [OA].

Gurtin, M. (1979). Thermodynamics and the cohesive zone in fracture. *Zeitschrift für angewandte Mathematik und Physik ZAMP* 30(6):991–1003. [DOI].

Halphen, B. and N. Quoc-Son (1975). Sur les matériaux standards généralisés. *Journal de Mécanique* 14:39–63. [HAL].

Jean, M. (1999). The non-smooth contact dynamics method. *Computer Methods in Applied Mechanics and Engineering* 177(3):235–257. [DOI], [HAL].

Jean, M., V. Acary, and Y. Monerie (2001). Non-smooth contact dynamics approach of cohesive materials. *Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences* 359(1789):2497–2518. [DOI], [HAL].

Kachanov, L. (1986). *Introduction to continuum damage mechanics*. Springer. [DOI].

Kondo, D., H. Welemane, and F. Cormery (2007). Basic concepts and models in continuum damage mechanics. *Revue Européenne de Génie Civil* 11(7-8):927–943. [DOI], [HAL].

Kuna, M. and S. Roth (2015). General remarks on cyclic cohesive zone models. *International Journal of Fracture* 196(1):147–167. [DOI].

Lemaitre, J. (1996). *A Course on Damage Mechanics*. Springer. [DOI].

Ma, Y. and H. Huang (2018). DEM analysis of failure mechanisms in the intact Brazilian test. *International Journal of Rock Mechanics and Mining Sciences* 102:109–119. [DOI].

Moreau, J. J. (1988). Unilateral Contact and Dry Friction in Finite Freedom Dynamics. *Nonsmooth Mechanics and Applications*. Ed. by J. J. Moreau and P. D. Panagiotopoulos. Springer, pp 1–82. [DOI], [HAL].

Moës, N., C. Stolz, P.-E. Bernard, and N. Chevaugeon (2011). A level set based model for damage growth: The thick level set approach. *International Journal for Numerical Methods in Engineering* 86(3):358–380. [DOI], [HAL].

Na, S., W. Sun, M. Ingraham, and H. Yoon (2017). Effects of spatial heterogeneity and material anisotropy on the fracture pattern and macroscopic effective toughness of Mancos Shale in Brazilian tests. *Journal of Geophysical Research: Solid Earth* 122(8):6202–6230. [DOI], [OA].

Needleman, A. (1990). An analysis of tensile decohesion along an interface. *Journal of the Mechanics and Physics of Solids* 38(3):289–324. [DOI].

Onsager, L. (1931). Reciprocal relations in irreversible processes. I. *Physical Review* 37 (4):405–426. [DOI], [OA].

Ortiz, M. and A. Pandolfi (1999). Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. *International Journal of Numerical Methods in Engineering* 44(9):1267–1282. [DOI], [HAL].

Park, K., G. Paulino, and J. Roesler (2009). A unified potential-based cohesive model of mixed-mode fracture. *Journal of the Mechanics and Physics of Solids* 57(6):891–908. [DOI].

Quey, R., P. Dawson, and F. Barbe (2011). Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing. *Computer Methods in Applied Mechanics and Engineering* 200(17):1729–1745. [DOI], [HAL].

Rabotnov, Y. (1969). *Creep Problems In Structural Members*. Series in Applied Mathematics and Mechanics. North-Holland.

Richefeu, V., A. Chrysochoos, V. Huon, Y. Monerie, R. Peyroux, and B. Wattrisse (2012). Toward
local identification of cohesive zone models using digital image correlation. *European Journal of Mechanics-A/Solids* 34:38-51. [DOI], [HAL].

Roe, K. and T. Siegmund (2003). An irreversible cohesive zone model for interface fatigue crack growth simulation. *Engineering Fracture Mechanics* 70(2):209-232. [DOI].

Serpieri, R., E. Sacco, and G. Alfano (2015). A thermodynamically consistent derivation of a frictional-damage cohesive-zone model with different mode I and mode II fracture energies. *European Journal of Mechanics-A/Solids* 49:13-25. [DOI], [OA].

Shu, W. and I. Stanciuşescu (2020). Fully coupled thermo-mechanical cohesive zone model with thermal softening: Application to nanocomposites. *International Journal of Solids and Structures* 188-189:1-11. [DOI], [OA].

Wojtacki, K., L. Daridon, F. Dubois, N. Moës, and Y. Monerie (2017). Analyse comparative de trois méthodes performantes de simulation numérique de la fissuration. *13e colloque national en calcul des structures* (Giens, France, May 15, 2017-May 19, 2017). [HAL].

Zener, C. (1938). Internal friction in solids II. General theory of thermoelastic internal friction. *Physical Review* 53 (1):90-99. [DOI], [HAL].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** AC and LD conceived the presented ideas. MR performed the numerical computations. All the authors wrote the paper and discussed the results.

**Supplementary Material** Animation of the crack propagation during Brazilian test and r2022 IMGc90 source files for the proposed simulations.

**Acknowledgements** All simulations are performed with the OpenSource IMGc90 platform. The authors sincerely acknowledge Frédéric DUBOIS and the IMGc90 team from Laboratoire de Mécanique et de Génie Civil de l'Université de Montpellier.

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.