# An extended energy-based cohesive zone model to depict anisotropic degradation of elastic-brittle interface

André Chrysochoos, Loïc Daridon, Mathieu Renouf

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

## To cite this version:

André Chrysochoos, Loïc Daridon, Mathieu Renouf. An extended energy-based cohesive zone model to depict anisotropic degradation of elastic-brittle interface. 2026. hal-05503492

HAL Id: hal-05503492

https://hal.science/hal-05503492v1

Preprint submitted on 10 Feb 2026

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An extended energy-based cohesive zone model to depict anisotropic degradation of elastic-brittle interface

Chrysochoos Andréª, Loïc Daridonª,*, Renouf Mathieuª

ªLMGC, University of Montpellier, CNRS, Montpellier, France

## Abstract

This paper presents a cohesive zone model (CZM) in which a scalar state variable characterizing damage evolution is used to depict an anisotropic degradation of elastic properties. The proposed model extends a previous modelling work where the damage criterion was based on the energy balance form associated with cohesive zone degradation and particularly on the evolution of the maximum elastic energy that can be stored for a given damage state. This extended model allows the description of various interactions between normal and tangential stiffness degradation. This model relies on a single scalar damage variable combined with an energy-based formulation to capture various types of coupling between normal and tangential damage mechanisms. The model's capabilities and limitations to consider different normal and tangential separation curve are then illustrated through a few simulations. Following in the footsteps of Rice and Willis, the paper then compares the local CZ approach to damage and the mechanics of brittle fracture, which introduces the concept of energy release rate proposed by Griffith. For this purpose, various numerical bending tests inspired by conventional

*Corresponding author: loic.daridon@umontpellier.fr

Preprint submitted to European Journal of Mechanics A

February 5, 2026
tests (DCB, ENF, MMB) are performed by adjusting several parameters related to the normal and tangential characteristics of the CZ. When a stable propagation of the crack is observed, the energy released rate, estimated macroscopically via an area method, is compared with the dissipated energy developed by the CZs.

***Keywords:*** Anisotropic damage effects, brittle fracture, cohesive zone, thermodynamics of irreversible processes, energy balance, finite element analysis.

# 1. Foreword

The following lines aim to position the present contribution within the broader legacy of Professor Jean Lemaitre's seminal work. Jean Lemaitre (1948–2019) was a pioneer in the development of continuum damage mechanics applied to materials. His work profoundly influenced the modeling of ductile materials—particularly metals—under mechanical and thermal loading. By embedding his theory within the thermodynamics of irreversible processes and introducing internal state variables, he created a coherent framework linking damage to material free energy and its evolution to the second law of thermodynamics.

One of his major achievements was the coupling of damage with viscoplasticity, making material degradation sensitive to the accumulation of plastic strain—a key aspect for predicting failure in ductile materials. His models have since been widely implemented in finite element software to simulate damage in metallic structures such as turbines, reactors, and pipelines. Several of his landmark publications are cited in the bibliography at the end of
17 this paper.

18 Beyond his scientific achievements, it is Jean's remarkable human quali-
ties that we would like to emphasize. Together with his friend and colleague
19 Jean-Louis Chaboche, he co-led the renowned Greco 47 initiative "Grandes
20 Déformations en Endommagement," which later evolved into the Mecamat
21 association. The first author of this paper had the privilege of witnessing
22 the formative years of the French "Mechanics of Materials" community—an
23 era that saw the convergence of mechanicians, materials physicists, applied
24 mathematicians, computational scientists, academics, and industry partners.
25 Jean was an exceptional communicator and unifier, playing a central role in
26 shaping the French solid mechanics community from the early 1980s onward.
27

28 This work builds on the foundation laid by Jean Lemaître, focusing on the
damage to brittle elastic materials. Although we adopt a different modelling
29 strategy (the cohesive zone approach), the damage processes are still treated
30 as irreversible phenomena within a thermodynamic framework.
31

32 Thanks, Jean.

## 2. Introduction

34 In an age of energy and raw materials savings, the durability of structures
35 is becoming a major economic and ecological issue. Whether the irreversibil-
36 ity of material behavior is induced by mechanical fatigue or physicochemical
37 ageing, the progressive degradation of loaded structures inevitably leads to
38 the formation of cracks which, as they propagate, cause their ruin.

39 Over the past decades, a variety of modeling frameworks have been de-
40 veloped to describe damage and rupture, reflecting the intrinsic complex-
ity, irreversibility, and multiscale nature of these phenomena. Among them, 41 Continuum Damage Mechanics (CDM), cohesive zone models (CZM), and 42 variational phase-field approaches constitute the most influential paradigms. 43

Continuum Damage Mechanics originated from the pioneering work of Kachanov [1],[2], who introduced a continuous damage variable to represent the loss of effective load-bearing area. This concept was subsequently formalized within a thermodynamically consistent framework by Lemaitre [3],[4], who established CDM as a constitutive theory based on internal variables, free energy potentials, and dissipation inequalities. Extensions by Chaboche and others [5] enabled the coupling of damage with plasticity and viscoplasticity, making CDM particularly effective for modeling diffuse degradation processes in ductile materials. 51 52

However, classical CDM formulations rely on local evolution laws, which leads to severe limitations when stress softening occurs. In quasi-brittle materials, damage localizes into narrow zones, resulting in ill-posed boundary value problems, mesh dependency, and an inability to reproduce experimentally observed size effects. These issues motivated the introduction of spatial regularization mechanisms. A major advance was achieved through non-local damage theories, notably developed by Bažant and Pijaudier-Cabot [6]. By replacing local damage-driving variables with spatially averaged quantities, non-local CDM introduces an intrinsic material length scale that regularizes localization and restores well-posedness. Closely related gradient-enhanced formulations were later proposed, providing an alternative but equivalent way to embed non-local interactions within the continuum framework [7]. 53 54 55 56 57 58 59 60 61 62 63

Beyond localization issues, the modeling of damage raises fundamental ques- 65
tions related to irreversibility and dissipation. In this respect, the work of Frémond [8] represents a decisive contribution. Through the framework of non-smooth thermomechanics, damage is described as an irreversible process governed by energy balances, dissipation potentials, and variational inequalities, providing a rigorous thermodynamic foundation for models involving unilateral constraints and non-differentiable energies.

A further paradigm shift was introduced by Marigo with the variational re-
formulation of fracture. In contrast to classical approaches based on local
criteria, fracture is treated as a global energy minimization problem sub-
ject to irreversibility constraints. This perspective was formalized in the
seminal work of Francfort and Marigo [9], which revisited Griffith's theory
[10] within a variational framework. Crack initiation and propagation then
emerge naturally from energy competition, without ad hoc propagation laws.
This variational viewpoint led to the development of phase-field models, in
which cracks are represented by a continuous scalar field regularized over a
finite length scale. These models can be interpreted as gradient-enhanced
damage formulations with a strong energetic basis and have proven highly
effective for simulating complex fracture phenomena [11].

Today, CDM, CZM, and phase-field approaches should be viewed as com-
plementary descriptions operating at different levels of localization. In this
general context, the present work focuses on a new type of CZM. Unlike con-
ventional cohesive zone models (CZMs), which rely on traction-separation
laws and stiffness degradation, the present model is formulated within an
energy-based framework. It is based on the evolution of the cohesive zone's
capacity to store elastic energy, which directly reflects its damage state. This
91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

damage state is quantified through a damage variable defined by an implicit
constitutive equation.

Proposed by Dugdale [12], Barenblatt [13], and then Rice [14] and Willis
[15], CZMs considered initially the normal forces exerted on the lips of a
crack [12], then the tangential forces linked to degradation [16], and later
those linked to friction of lips once the crack has formed, considering the
unilateral aspects inherent in non-penetration [17], [18]. CZMs are widely
used in composite materials to capture the R-curve effect, which arises from
fiber bridging between crack tips. Crack propagation can then be simulated
by combining the global fracture response with local energy dissipation within
the cohesive zone [19].

From the 80s onwards, these CZMs were sometimes presented in a thermo-
dynamic framework. Gurtin [20] was probably the first to consider cohesive
zones as a thermodynamic system endowed with a potential, (e.g., free en-
ergy), dependent on the crack temperature and the crack tip opening dis-
placement. Naturally, in this framework, the status of state variable was
given to a scalar (or tensorial) parameter depicting the damage progress.
Regarding the evolution laws, and particularly those related to the damage
variable, some, like Costenzo and Allen [21] advocate the use of a dissipation
potential as usually done in the Generalized Standard Materials (GSM) the-
oretical framework [22]. Others, probably more numerous, proposed forms
of differential (kinetic) equations compatible with the 2nd principle of ther-
modynamics. Although it is impossible to mention all the works related to
how damage kinetics were constructed, a previous paper [23] reminded some
references spanning the last twenty years: e.g. [24, 25, 26, 27, 28, 29], and
more recently [30].

In our previous paper [23], it was shown that the damage evolution law
is automatically defined as soon as this damage is assumed to be the only
irreversible mechanism and is fully dissipative. These assumptions are par-
ticularly well suited to describing elastic brittle media whose cracking, at
the CZ level, does not induce internal energy storage. The main physical
ingredient of this model is the maximum elastic energy that can be stored
for a given level of damage. Once the (constant) ratio of tangent stiffness to
normal stiffness is known, the normal and tangential separation curves can
be reconstructed. For such a constant ratio, a damageable isotropic elastic
CZM was then considered for isothermal transformations and applied to two
standard tests: the bending of a pre-cracked beam and the Brazilian test.
For monotonous loadings, the evolution of the cracked zones and the form of
the energy balances around the cracks were briefly shown.

The current work presents in a first part an extension of this previous
model. The basic ingredient is still the maximum storable elastic energy (for
a given damage level), but this time the ratio of tangential to normal stiffness
becomes an implicit function of the damage state variable. This degree of
freedom then allows depicting anisotropic degradations of elastic properties.
Several numerical exercises then illustrate the potential of this approach. For
the sake of simplicity, we restrict ourselves to examples who consider trans-
verse isotropy, the isotropy plane being the plane of the CZ. By adjusting
both the evolution of the maximum storable elastic energy and that of the
stiffness ratio as a function of damage, these examples show the wide variety
of normal and tangential separation curves that can be obtained, depict-
ing different interaction modes between tangential and normal degradation
effects.

The second part of the article examines some links between fracture me-
chanics (global scale) and the cohesive zone approach (local scale) from a
energy standpoint. In the past, famous works have been carried out on this
subject (e.g. [13, 14, 31]). We also want to compare both approaches in terms
of the dissipated energy induced by the advance of the crack(s) and the pro-
gressive degradation of cohesive zones close to the crack tip, for elastic brittle
materials. Indeed the dissipation in fracture mechanics is associated with ge-
ometrical state variables depicting the crack grows while the cohesive zones
may locally dissipate energy without showing at a more macroscopic scale
a crack onset or progress. Furthermore, how do these elastic and dissipated
energies vary according to the normal and tangential properties of the CZ
itself? What are the dynamic effects (i.e. dynamic crises) when a crack
suddenly expands, even under quasi-static loading, and what is their energy
cost? Finally, can we systematically define a macroscopic energy quantity
and the scale at which this latter is defined, that can be sufficiently repre-
sentative of the crack progress and be considered as the strain energy release
rate?

For this purpose, numerical bending tests close to standard experimental
ones are performed : crack propagation in mode I (double cantilever beam,
DCB), mode II (bending with end-notch, ENF), or modes I and II (mixed
mode bending mode, MMB) [32]. In this particular situation, a pre-crack
exists from the very start of the test, allowing Griffith's theory to be used without any initiation problem. Tests are made by varying several parameters related to the normal and tangential characteristics of the CZs. In each situation, global (i.e. beam scale) and local (i.e. CZs scale) energy balance are performed. Computations are developed within the framework of plane strain brittle elasticity.

The layout of this paper is as follows. In Section 3, the CZM and its underlying physical assumptions are presented, using the framework of the thermodynamics of irreversible processes (TIP). The damage rate and its associated thermodynamic force are then defined on the basis of the maximum storable elastic energy criterion defined for a given damage state. Several numerical examples show the responses of CZ for different interactions between normal and tangential stiffnesses. In Section 4, the energy concepts related to fracture mechanics are briefly reviewed. Section 5 presents the computational tools used. We essentially use LMGc90, an open-source numerical platform for the simulation of interacting multi-body systems, with a particular emphasis on non-smooth contact dynamics involving friction and impacts [33, 34]. Its implicit formulation ensures robust and stable treatment of complex contact interactions without resorting to penalty-based regularization. In Section 6, the energy balances associated with the various bending tests are systematically presented and discussed. The paper concludes by reviewing the contributions and limitations of this energy-based approach to formulating cohesive zone models and proposes several ways to improve future works.
### 3. Cohezive zone model

#### 3.1. Mechanical variables

Regarding the standard mechanical description of the CZ, the traction force **R** and the separation **u** naturally have 3 components. For the sake of simplicity, we introduce in what follows only 2 components, the normal and the tangential ones denoted by $(R_n, R_t)$ for the traction vector and $(u_n, u_t)$ for the separation vector. This is equivalent to considering the CZ as having a transverse isotropic behavior, with the isotropy plane corresponding to that of the CZ. Moreover, as is conventionally admitted for CZM, the normal displacement jump $u_n$ is non negative. This unilateral condition is taken into account by a Signori type relationship in the numerical simulations performed with LMGc90 [35] in the final part of this work. Although the code handles contact and potential friction of cracks, the latter is not taken into account to limit dissipation to damage mechanisms.

#### 3.2. Thermodynamic potential and state equations

Working in a thermomechanical framework implies first to define a set of state variables. Again, for simplicity's sake, we consider only isothermal processes. We then select the components $(u_n, u_t)$ of the separation vector and a scalar damage variable denoted by $u_d$. This variable, through its evolution, has to be able to depict the irreversible nature of the damage progress. As explained in [23], the Helmholtz free energy is chosen as thermodynamic potential, and in case of isothermal processes, this potential is classically identified to the elastic energy :
$$
\begin{align}
\psi(\mathbf{u}, u_d) &= w_e(\mathbf{u}, u_d), \nonumber \\
&= \frac{1}{2}(K_n(u_d)u_n^2 + K_t(u_d)u_t^2), \nonumber \\
&= \frac{1}{2}K_n(u_d)(u_n^2 + \alpha(u_d)u_t^2), \nonumber \\
&= \frac{1}{2}K_n(u_d)u_{eq}^2, \tag{1}
\end{align}
$$

212 where $K_n$ and $K_t$ are the normal and tangential secant stiffnesses and where
213 $u_{eq}$ is defined by :

$$
u_{eq} = (u_n^2 + \alpha(u_d)u_t^2)^{\frac{1}{2}} . \qquad (2)
$$

214 In the last writing of the free energy Eq. (1), an equivalent displacement
215 jump $u_{eq}$ appears, leading us to define the state variable $u_d$ representative of
216 the damage progress as:

$$
u_d = \max \{u_{eq}(\tau), \forall \tau \leq t\} . \quad (3)
$$

Indeed, the goal is to quantify the damage progress via the CZ aptitude to store the elastic energy. The state variable $u_d$ defined in Eq. (3) takes the 2D aspect of the separation vector $\mathbf{u}$ into account. By construction, $\dot{u}_d$ is non-negative and de facto respects the irreversibility of the damage progress. Note also that under a monotonous separation-controlled loading such that $\dot{u}_n \ge 0$ and $\dot{u}_t \ge 0$ (or $\dot{u}_t \le 0$), and whatever the loading direction, $u_d = u_{eq}$. The ratio $\alpha = K_t/K_n$, introduced in Eq. (1), is defined as a function of the damage state variable $u_d$ and is, by construction, strictly positive ($\alpha \in \mathbb{R}^+$). In the case of isotropic damage, $\alpha$ remains constant [23]. Allowing $\alpha$ to vary with $u_d$ enables the model to capture different levels of interaction between
normal and tangential degradation mechanisms. For instance, a significant
decrease in $\alpha$ results in a rapid reduction of $K_t$ even under purely normal
loading ($\dot{u}_n \ge 0$, $\dot{u}_t = 0$), illustrating a strong coupling between the normal
and tangential responses.

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

$$
\left\{
\begin{aligned}
R_n^r &= \frac{\partial \psi}{\partial u_n} = K_n(u_d)u_n \\
R_t^r &= \frac{\partial \psi}{\partial u_t} = K_t(u_d)u_t = \alpha(u_d)K_n(u_d)u_t \\
A_d &= \frac{\partial \psi}{\partial u_d} = \frac{1}{2}K'_n(u_d)u_n^2 + \frac{1}{2}K'_t(u_d)u_t^2 \\
    &= \frac{1}{2}K'_n(u_d)u_{eq}^2 + \frac{1}{2}\frac{\alpha'(u_d)}{\alpha(u_d)}K_t(u_d)u_t^2
\end{aligned}
\right.
\qquad (4)
$$

where $K'_n(u_d) = \frac{dK_n}{du_d}$ and $K'_t(u_d) = \frac{dK_t}{du_d}$.

3.3. *Taking irreversibility into account*

The irreversibility of mechanisms accompanying the opening of the co-
hesive 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 :

$$
\begin{align}
w_d^o &= w_{def}^o - \dot{\psi} = \mathbf{R} \cdot \dot{\mathbf{u}} - \frac{\partial \psi}{\partial \mathbf{u}} \cdot \dot{\mathbf{u}} - \frac{\partial \psi}{\partial u_d} \dot{u}_d \tag{5} \\
&= \mathbf{R}^{ir} \cdot \dot{\mathbf{u}} + X_d \dot{u}_d \ge 0 \nonumber
\end{align}
$$

The terms $w_d^o$ and $w_{def}^o$ determine the dissipated and deformation energy
rates, respectively. The symbol $(-)^o$ underlines that $w_d$ and $w_{def}$ are not a
priori state functions and are then path-dependent. Eq. (5) also introduce the
irreversible part of the traction force, $\mathbf{R}^{ir} = \mathbf{R} - \mathbf{R}^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 $\mathbf{R}^{ir} = \mathbf{0}$. The traction force $\mathbf{R}$ is then derived from the potential $\psi$:

$$ \mathbf{R} = \frac{\partial \psi}{\partial \mathbf{u}} = \begin{pmatrix} K_n(u_d)u_n \\ K_t(u_d)u_t \end{pmatrix}. \qquad (6) $$

Consequently, the intrinsic dissipation is reduced to:

$$ \begin{aligned} w_d^o &= X_d \dot{u}_d = -\frac{1}{2}(K'_n u_n^2 + K'_t u_t^2) \dot{u}_d \\ &= -\frac{1}{2}(K'_n u_{eq}^2 + \frac{\alpha'}{\alpha} K_t u_t^2) \dot{u}_d \quad \geq 0. \end{aligned} \qquad (7) $$

Since the rates $\dot{u}_d$ and $w_d^o$ are non-negative for any loading path (i.e. for all $(u_n, u_t)$), it follows that $K'_n \leq 0$ and $K'_t \leq 0$, which is physically consistent. The irreversible nature of the damage leads to a degradation of the secant stiffnesses. The derivative $\alpha'$ is therefore also expected to be negative or null.

In general, at this stage of model presentation, the complementary laws depicting the evolution of the thermodynamic system have to be defined. As already mentioned, within the GSM framework, a dissipation potential is introduced with the right convexity properties, so as to ensure the positivity of the Clausius-Duhem inequality whatever the system evolution. In [23], we showed that this approach is no longer necessary. The damage evolution law is governed by energy-based hypotheses. Accordingly, the current elasticity domain is defined through the maximum storable elastic energy $w_e^d$,
independently of the loading path.

$$w_e(\mathbf{u}, u_d) \le w_e^d(u_d) = w_e(\mathbf{u}^*, u_d), \text{ when } u_{eq}^* = u_d. \quad (8)$$

The evolution of $u_d$ stems from the fact that, for damage to occur, the maximum elastic energy that can be stored in the CZ must remain on the threshold during the loading step, i.e.

$$\left\{ \begin{array}{l} w_e(\mathbf{u}, u_d) = w_e^d(u_d) \quad (a) \\ \dot{w}_e(\mathbf{u}, u_d) = \dot{w}_e^d(u_d) \quad (b) \end{array} \right. . \qquad (9)$$

Equations Eq. (9)(a,b) can be rewritten in the following compact form, exhibiting a formal structure similar to that of the Signorini conditions in unilateral contact problems:

$$(w_e(\mathbf{u}, u_d) - w_e^d(u_d)) \dot{u}_d = 0. \qquad (10)$$

Let us also note that Eq. (9)(a) provides a non-classical implicit relation for the damage variable $u_d$. For a given pair $(u_n, u_t)$ such that $u_{eq}(u_n, u_t) = u_d$, we obtain:

$$u_d = (u_n^2 + \alpha(u_d)u_t^2)^{\frac{1}{2}}. \qquad (11)$$

This definition of $u_d$ is naturally employed in the numerical examples presented hereafter to illustrate both the capabilities and the limitations of this type of CZM.
### 3.4. Numerical illustrations of damage-induced normal-tangential coupling

First, let us recall once more the energy hypothesis underlying the CZMs considered here. The maximum elastic energy that can be stored by the CZ for a given damage level is a predefined function $w_e^d$ of the damage variable $u_d$, Eq. (8). This means that, regardless of the loading path supported by the CZ during its lifetime, it will pass through the same levels of stored elastic energy.

For the sake of simplicity, the following form of $w_e^d(u_d)$ has been chosen as:

$$w_e^d(u_d) = \frac{1}{2} K_n^0 \left( 1 - \left( \frac{u_d}{u_n^c} \right)^{r_n} \right) u_d^2 \quad \text{with } u_d \le u_n^c, \qquad (12)$$

where $K_n^0$, $u_n^c$ and $r_n$ are material constants. Of course, other analytical expressions or experimental fitting could have been chosen. However, the form given in Eq. (12) is particularly convenient: it allows us to describe an initial phase of growth of the maximum storable elastic energy, followed by a phase of decay until all elastic energy vanishes. The parameter $K_n^0$ controls the elastic energy magnitude, $u_n^c$ limits the displacement jump, and the exponent $r_n$ advances or delays the degradation mechanism. The $n$ index has been chosen because for a monotonous normal separation, $u_d = u_{eq} = u_n$, with $u_t \equiv 0$. The normal secant stiffness can also be written as:

$$K_n(u_d) = K_n^0 \left( 1 - \left( \frac{u_d}{u_n^c} \right)^{r_n} \right) . \qquad (13)$$

Instead of introducing an a priori $\alpha(u_d)$, it seemed more intuitive to introduce the elastic energy stored under a monotonous tangential separation-
controlled loading $w_e^t(u_t)$, considering an analogous form to those in Eq. (12):

$$w_e^t(u_t) = \frac{1}{2} K_t^0 \left( 1 - \left( \frac{u_t}{u_t^c} \right)^{r_t} \right) u_t^2 . \qquad (14)$$

The constants $K_t^0$, $u_t^c$ and $r_t$ must then be set, taking care to ensure that the range of elastic energy $w_e^d$ remains the same whatever the type of loading considered in Eq. (12). In particular, if we compare a monotonous pure normal separation-controlled loading with a pure tangential one, we must have:

$$\max(K_n(u_d)u_n^2) = \max(K_t(u_d)u_t^2) \quad \text{with } u_n \in [0, u_n^c] \text{ and } |u_t| \in [0, u_t^c]. \qquad (15)$$

This property ensures a one-to-one relationship between the normal and tangential separation components during the increasing and then decreasing phase of $w_e^d(u_d)$. The correspondence is made in terms of elastic energy and for each $u_t$ it is then possible to find a unique $u_d$ value such that :

$$\overbrace{K_t^0 \left(1 - \left(\frac{u_t}{u_t^c}\right)^{r_t}\right)}^{K_t(u_d)} u_t^2 = K_n(u_d) u_d^2, \qquad (16)$$

and to derive the function $K_t(u_d)$. For the particular expressions of the secant stiffnesses (Eq. (13), Eq. (16)), the crucial constraint for ensuring a bijective relationship, between $u_t$ and $u_d$ under monotonous loadings, is as follows:

$$K_n^0 \frac{r_n}{r_n + 2} \left( \frac{2}{2 + r_n} \right)^{\frac{2}{r_n}} (u_n^c)^2 = K_t^0 \frac{r_t}{r_t + 2} \left( \frac{2}{2 + r_t} \right)^{\frac{2}{r_t}} (u_t^c)^2 . \qquad (17)$$
The relation Eq. (16), which implies Eq. (15), formalizes the assumption of a unique function $w_e^d(u_d)$ defining the elastic domain of the cohesive zone, independently of the loading path. This formulation establishes a consistent link between the normal and tangential damage evolution, relying on a single scalar damage variable. If two of the three parameters defining the stiffness have been selected, the third one is determined by Eq. (17), which nevertheless allows the description of a wide range of relevant cases. In the following examples, the maximum tangential separations, $u_n^c$ and $u_t^c$, have been arbitrarily fixed. The initial stiffnesses, $K_n^0$ and $K_t^0$, and the exponents, $r_n$ and $r_t$, were then adjusted using Eq. (17) to obtain, in each case, a maximum elastic energy of $1 \text{ J m}^{-2}$.

### 3.4.1. CZM parameters and characteristic curves

Let's take a look at three examples illustrating, for a given loading, different forms of interaction between normal and tangential degradation.

Table 1 lists the parameters used in the three numerical tests presented hereafter. The first cohesive zone model (CZM1) describe an isotropic degradation of the elastic properties. The ratio $\alpha(u_d) = K_n^0/K_n$ remains then constant during the loading. For the second and third cohesive zone models, we simply swapped the exponents $r_n$ and $r_t$, modifying the initial stiffnesses $K_n^0$ and $K_t^0$ accordingly. The aim here is to see the effects of different interactions between normal and tangential degradation effects. For CZM2, tangential damage is triggered first ($\alpha(u_d)$ decreasing) while for CZM3, normal damage is the first to occur ($\alpha(u_d)$ increasing).
<table><caption>Table 1: Parameter values of CZM</caption><thead><tr><th>CZM</th><th>K<sub>n</sub><sup>0</sup> (N/m<sup>3</sup>)</th><th>u<sub>n</sub><sup>c</sup> (m)</th><th>r<sub>n</sub></th><th>K<sub>t</sub><sup>0</sup> (N/m<sup>3</sup>)</th><th>u<sub>t</sub><sup>c</sup> (m)</th><th>r<sub>t</sub></th></tr></thead><tbody><tr><td>1</td><td>2.3 × 10<sup>12</sup></td><td>1.5 × 10<sup>-6</sup></td><td>4</td><td>0.83 × 10<sup>12</sup></td><td>2.5 × 10<sup>-6</sup></td><td>4</td></tr><tr><td>2</td><td>2.3 × 10<sup>12</sup></td><td>1.5 × 10<sup>-6</sup></td><td>4</td><td>2.9 × 10<sup>12</sup></td><td>2.5 × 10<sup>-6</sup></td><td>0.7</td></tr><tr><td>3</td><td>8.1 × 10<sup>12</sup></td><td>1.5 × 10<sup>-6</sup></td><td>0.7</td><td>0.83 × 10<sup>12</sup></td><td>2.5 × 10<sup>-6</sup></td><td>4</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>(a) Isotropic degradation</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Tangential degradation priority</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Normal degradation priority</figcaption></figure>

<caption>Figure 1: Normal and tangential elastic energies for the three CZM</caption>

<footnote n="334">334 FIGURE 1 shows the elastic energies corresponding to the three cohesive zones.</footnote>

<footnote n="335">335 The separations at crack initiation are arbitrarily the same in all cases, the tan-</footnote>
gential crack separation $u_i^c$ being chosen larger than the normal one. Note that, according to Eq. (15), the maximum stored elastic energy (indicated by a circle) in FIGURE 1) is the same whether the loading is normal or tangential.

<figure><img src="image_6.png" /><figcaption>Figure 2: Normal and tangential separation curves for the three CZM</figcaption></figure>

Knowing the evolution of elastic energy for monotonous normal or tangential loading, it is possible to derive the corresponding separation curves (see FIGURE 2). Dotted straight lines were superimposed, with slopes corresponding to the initial normal (red) and tangential (blue) stiffnesses. The dissociation of these virtual elastic responses with the separation curves indicates the value of the displacement
ment jump ($u_n$ or $u_t$) at which damage begins to have a significant effect on the
mechanical response.

<figure><img src="image_7.png" /><figcaption>Figure 3: Normal and tangential stiffnesses for the three CZM</figcaption></figure>

FIGURE 3 shows the evolution of the stiffnesses, $K_n(u_n)$ and $K_t(u_t)$, for the three CZMs, clearly demonstrating the consequences of the chosen model parameters, *Table 1*. For CZM1, the parameters were chosen to produce an isotropic degradation of the elastic properties. In contrast, CZM2 exhibits a rapid decrease in tangential stiffness at low separations while the normal stiffness remains nearly constant. The opposite behavior is observed for CZM3, where the tangential
stiffness remains constant at low separations, and the normal stiffness decreases
rapidly.

3.4.2. Damage threshold evolution under loading

The separation-controlled loading chosen to illustrate the CZMs’ responses is
the following:

• stage 1: equibiaxial loading $0 \le u_n \le u_n^{(1)}$ and $0 \le u_t \le u_t^{(1)}$, with $\dot{u}_n = \dot{u}_t = \dot{u}_0 \ge 0$.

• stage 2: tangential unloading and reloading at constant normal separation
$u_n = u_n^{(1)}$ and $\dot{u}_t = -\dot{u}_0$ until the elastic domain boundary is reached again
in $(u_n^{(1)}, u_t^{(2)})$, i.e. $u_t^{(2)} \le u_t \le u_t^{(1)}$

• stage 3 continues stage 2 (i.e. $u_n = u_n^{(1)}$ and $\dot{u}_t = -\dot{u}_0$, ) until the crack onset of the CZ. In fact, in FIGURE 4 and in FIGURE 5, the elasticity domain (dotted blue 1/2 ellipse) is plotted a few iterations before the crack onset.

Indeed, when $u_d = u_n^c$, the elasticity domain of the CZM is reduced to zero in
the $(R_t, R_n)$ plane (see FIGURE 5); by contrast, it is maximal in the $(u_n, u_t)$ plane,
as illustrated in FIGURE 4.
<figure><img src="image_8.png" /><figcaption>Figure 4: Evolution of the elasticity domain in $(u_t, u_n)$ plane</figcaption></figure>

FIGURE 4 shows the evolution of the elasticity domain under the separation-controlled loading path. The continuous red straight line corresponds to stage 1 (i.e. $\dot{u}_n = \dot{u}_t = \dot{u}_0 \ge 0$). At the end of this stage, the ellipsoid boundary of the elasticity domain is plotted (continuous blue line). During stages 2 and 3 ($\dot{u}_n = 0, \dot{u}_t = -\dot{u}_0$), the red dashed line represents the path where the normal separation remains constant. Markers (black circles) indicate the beginning and the end of stage 2. A black triangle indicates the end of the loading. The numerical test has been ended before the crack onset. The corresponding boundary of elasticity domain was plotted (blue dashed line).
<figure><img src="image_9.png" /><figcaption>Figure 5: Evolution of the elasticity domain in $(R_t, R_n)$ plane</figcaption></figure>

FIGURE 5 display the same scenarios as in FIGURE 4, but this time in terms of traction forces $(R_t, R_n)$. During stage 1, a difference in curvature induced by the damage kinetics is visible: in FIGURE 5a, damage effects develop isotropically, so the normal and tangential forces $R_t, R_n$ remain proportional (no curvature). In FIGURE 5b, damage in the tangential direction appears first (positive curvature), while in FIGURE 5c, the opposite occurs (negative curvature). Once the maximum elastic energy has been reached, (see Eq. (16) or the blue or red circles at the top of the curves in FIGURE 1), the loading path lies outside the demi-ellipse represent-
ing the elasticity domain. During subsequent damaging loading (end of stage 3),
the forces progressively decrease, thereby reducing the size of the elasticity domain.

In FIGURE 4a and FIGURE 5a, the evolution of the surface of the elastic domain is homothetic, since these cases correspond to isotropic damage. In contrast, FIGURE 4b, FIGURE 4c and FIGURE 5b, FIGURE 5c show a non-homothetic evolution of the elastic domain, illustrating the coupling between normal and tangential degradations.

<figure><img src="image_10.png" /><figcaption>Figure 6: Energy balance associated with the separation-controlled loading path shown in FIGURE 4</figcaption></figure>
The energy balance associated with the loading previously described are shown for the three CZs in FIGURE 6. As already mentioned in Eq. (5), the deformation energy $w_{def}$ (green) is the sum of the stored elastic energy $w_e$ (blue) and the dissipated energy $w_d$ (red). The corresponding energy courses have been plotted as a function of the damage state variable $u_d$ and the 3 loading stages have been indicated. Note that stage 2 is carried out at constant damage ($\dot{u}_d = 0$) and thus results in a vertical line running downwards under tangential elastic unloading, then upwards when elastic loading is resumed (with $u_t \le 0$). At the beginning of loading, all CZMs tend to store elastically the deformation energy ($w_{def}^0 \approx w_e^0$, $w_d^0 \approx 0$), whereas at the end of the test, the deformation energy rate and the rate of released elastic energy (previously stored) are dissipated ($w_d^0 \approx w_{def}^0 - w_e^0$). Finally note that for this complex loading, graphs show that the elastic energy reaches in each case a maximum of $1 J m^{-2}$, as was arbitrarily proposed. However, the dissipated and deformation energies at failure differ from one model to another.

## 4. Brief reminders of fracture mechanics

Griffith approaches the problem of brittle fracture of cracked bodies from an energy standpoint [10]. He considers the solid to be elastic and when cracking occurs its geometry changes. In the case of isothermal processes, cracking is the only irreversible mechanisms and it is fully dissipative. This energy framework is therefore very similar to the one used to introduce CZMs in the previous sections. The main difference is that energy is dissipated if and only if the crack or cracks progress. Griffith then introduces the concept of energy release rate $G$. To define $G$, let us consider the kinetic energy theorem:

$$W_{ext}^o - W_{def}^o = \frac{dK}{dt} \qquad (18)$$
where we used capital letters to designate energies at the scale of the structure and where $W_{ext}$, $W_{def}$ and $K$ stand for the work of external forces applied on the body, the deformation energy and the kinetic energy, respectively. Neglecting the volume density of external forces, these energy terms are defined by:

$$ 
\begin{aligned}
W_{ext}(t) &= \int_{0}^{t} \int_{\partial\Omega(\tau)} t_i v_i d\partial\Omega d\tau \\
W_{def}(t) &= \int_{0}^{t} \int_{\Omega(\tau)} \sigma_{ij}\dot{\varepsilon}_{ij} d\Omega d\tau \\
K(t) &= \int_{\Omega(t)} \frac{1}{2}\rho v_i v_i d\Omega
\end{aligned} 
\quad (19) $$

• where $t_i$, ($i = 1, 2, 3$) are the components of the stress vector at the boundary $\partial\Omega(\tau)$ of the structure, $v_i$ are the component of the Galilean velocity vector,

• where $\sigma_{ij}$ and $\varepsilon_{ij}$ are usually the components of the stress and strain tensors,

• where $\rho$ is the mass density.

The deformation energy rate can be divided into dissipated and elastic parts, so that:

$$ W_{def}^{o} = W_{d}^{o} + W_{el}^{o} \quad (20) $$

So, whether the crack progression is stable (i.e. $\frac{dK}{dt} \approx 0$) or not, the associated dissipation can be expressed as the difference between the deformation energy rate and the elastic energy rate:

$$ W_{d}^{o} = W_{def}^{o} - W_{el}^{o} \quad (21) $$

If we now take into account the variation in energy dissipated in relation to the advance $\delta S$ of the cracked surface, we obtain $G$, the famous energy release rate defined by:
$$G = \frac{\delta}{\delta S}(W_{ext} - W_{el} - K) = \frac{\delta W_d}{\delta S} \qquad (22)$$

How the rate of dissipated energy evolves in relation to crack advancement, in comparison to the dissipation associated with all cohesive zones, will be addressed in Section 6, following a brief overview of the numerical tools used. Moreover, although quasi-static loading is considered in the following sections, we will estimate and compare the kinetic energy associated with dynamic crises, wherein cracks suddenly progress, with other terms in the energy balance.

## 5. Numerical Implementation

To demonstrate the potential of the proposed CZ models, simulations were carried out based on common experiments extracted from the literature (e.g. [32]). It is important to note that the objective of these simulations is merely to demonstrate the functionality of the cohesive zone models, rather than to optimise parameters in order to fit any benchmark. The numerical implementation of the previous models is then done in the code IMGe90 based on Non-Smooth Contact Dynamics (NSCD) [33, 34]. The NSCD method is dedicated to solving problems related to dynamic systems with unilateral constraints. It is therefore particularly suitable for contact friction problems. It proposes a non-smooth treatment (no compliance, no penalty) of the conditions of contact [34].

In this method, the treatment of adhesion enables the definition of a cohesive zone at each point of contact. Then the mechanical behaviour of the cohesive zones may vary at any point of the spatial discretization of the problem. This relevant modeling framework was then adopted to numerically simulate crack propagation with cohesive zone [36]. While the IMGe90 software is capable of modeling friction at crack edges, this effect was deliberately neglected in the present work to facilitate
itate the energy analysis, particularly in the case of ENF tests. Accordingly, the dissipation predicted by the numerical simulations originates exclusively from the degradation of the cohesive zones.

A two-dimensional finite element model is developed to investigate interfacial damage and debonding mechanisms in adhesively bonded beam assemblies. The structure consists of two elastic beams bonded along a common interface modeled using a cohesive zone formulation. The left end of the assembly is fully clamped, while various loading conditions are applied at the right end in order to activate different fracture modes.

The discretization is performed using a conforming unstructured mesh composed of 4325 linear triangular (T3) elements (FIGURE 7(a)). A cohesive interface is introduced along the bonding line between the two beams, governed by a traction-separation law derived from the proposed CZMs.

Three classical fracture configurations are considered to selectively activate the main interfacial failure modes: the Double Cantilever Beam (DCB) test for mode I (opening), the End-Notched Flexure (ENF) test for mode II (shear), and the Mixed-Mode Bending (MMB) test for combined mode I/II loading. These configurations allow a systematic assessment of the cohesive model response under pure and mixed-mode conditions, as well as a consistent comparison of the associated energy dissipation mechanisms. For the three loading configurations and the three CZMs considered (FIGURE 7), crack propagation is necessarily confined between the two beams, as the path is governed by cohesive zones introduced along the edges of selected mesh elements. The case of a pre-cracked homogeneous beam subjected to flexural loading, with cohesive zones implemented along all element boundaries, is not addressed hereafter and will be widely discussed in the Perspectives section.
<figure><img src="image_11.png" /><figcaption>Figure 7: Delamination tests and amplified deformed configurations; white-filled red circles indicate the displacement application points, while the cohesive interface is shown in red; (a) finite element mesh; (b) modified DCB test with symmetric opening displacements; (c) modified ENF test with identical vertical displacements applied to both beams; (d) modified MMB test with displacement applied only to the upper beam.</figcaption></figure>

# 6. Analysis of results and discussion

This section presents, in turn, the bending test specimens considered, the various loading characteristics and the kinematic boundary conditions imposed in order to eliminate any rigid body motion under loading. Then, it systematically presents and comments on the energy results obtained via simulations carried out on the three cohesive zone models and the three bending modes.

## 6.1. Material parameters, bending tests and loadings

The volume elastic properties of the material are 200 GPa for the Young's modulus and 0.3 for the Poisson's ratio. The mass density is $2.5 \times 10^3 \text{ kg m}^{-3}$. The samples are 18 mm long and 2 mm height.

Basic sketch of the bending specimens are shown in FIGURE 8. The name of a bending test includes the flexural mode and the type of cohesive zone selected. For example, ENF2 denotes an end-notched flexural test using the CZM parameters listed in Table 1 (line 2), referred to as CZM2. For these initial exploratory tests, the cohesive zones were drastically limited on a line placed halfway along
the longitudinal axis of the flexural test specimen. The finite element mesh is cre-
ated accordingly. The pre-crack is 0.5 mm for BCD tests. This was extended to
1.5 mm for ENF tests and 2.5 mm for MMB tests to avoid total instability of crack
propagation.

<figure><img src="image_12.png" /><figcaption>Figure 8: Basic sketch of bending specimens and tests. The black line symbolizes the locations where CZs are placed.</figcaption></figure>

Tests are displacement-controlled and the displacement rate is about $4 \times 10^{-5} \text{ mm s}^{-1}$.

The beams are embedded at their left ends to prevent any rigid body displacement.

6.2. DCB tests

FIGURE 9 shows the time evolution of the load $F_{kn}$, the process zone length $a_{pz}$, and the crack tip position $a_{ct}$ for the three CZMs. From the very beginning of loading, the process zone $a_{pz}$ advances steadily as the crack tip remains stationary. Then, the process zone and crack tip move forward together, with the size of the process zone ($a_{pz} - a_{ct}$) remaining almost constant. As usual, for elastic brittle materials, the force increases at the start of loading as long as the crack does not
propagate and then gradually decreases as the crack progresses. The evolution curves are similar for strength, damage area and cracks, regardless of the cohesive zone model. For CZ2, the maximum strength is slightly higher, which is logical given the parameters used to model normal behaviour. We will come back to this point a little later.

<figure><img src="image_13.png" /><figcaption>Figure 9: Load, process zone and crack tip time courses</figcaption></figure>

FIGURE 10 shows the evolution of the energy balance for each of the three cohesive zone families. These balances are plotted as a function of the regular process zone advance, unlike the irregular crack tip progression. The work developed by the cohesive forces ($w_{def}$), the stored elastic energy ($w_e$) and the dissipated energy ($w_d$) are shown, this latter being broken down into normal ($w_{dn}$) and tangential ($w_{dt}$) components according to the following form of Eq. (7):

$$w_d^o = -\frac{1}{2} K'_n u_n^2 u_d - \frac{1}{2} K'_t u_t^2 u_d = w_{dn}^o + w_{dt}^o. \quad (23)$$

As previously, the results associated with the three types of cohesive CZs are similar. At the very beginning of the tests, the first cohesive zones subjected to loading respond quasi-elastically, resulting in an equilibrium between the work developed by the cohesive forces and the stored elastic energy. Then, an equilibrium state is progressively established between the CZs that are reaching rupture and
those that are beginning to be elastically stressed. These particular CZs correspond to the process zone whose size remains constant throughout the test. This leads to stagnation of the elastic energy associated with all CZs and to an almost parallel increase in dissipation $w_d$ with the work of cohesive forces $w_{def}$. Finally, it should be noted that in all three cases, the tangential dissipation $w_{dt}$ remains zero, which means that only the normal forces of the cohesive zones dissipate, corresponding to crack propagation in mode I, a well-known result for DCB-type tests [32].

<figure><img src="image_14.png" /><figcaption>Figure 10: Local energy balance</figcaption></figure>

FIGURE 11 shows the global energy balances. 'Global' means that the energies are associated with the entire structure. As already mentioned, IMGc90 code considers dynamic situations. The kinetic energy $K$ is therefore computed. In the case of a quasi-static loaded DCB test, the kinetic energy remains low, though not zero. It is multiplied by 1000 to observe its evolution during the test. The other terms of the energy balance are the work of external forces $W_{ext}$, the dissipated energy $W_d$ and the elastic energy $W_{el}$. Insofar as the volumetric behavior of the structure is purely elastic, only the stressed CZs dissipate energy, i.e. $W_d = w_d$. Furthermore, the total elastic energy $W_{el}$ includes the volumetric deformation energy of the structure plus the elastic energy of the CZs. As can be seen from the local energy balance, crack propagation is in mode I. Given the shape of the mesh
used, it is primarily the 'normal' characteristics of the cohesive zones that come into play. Therefore, it is only natural that the CZM2 shows greater resistance to propagation than the other two models.

<figure><img src="image_15.png" /><figcaption>Figure 11: Global energy balance</figcaption></figure>

Given the levels of kinetic energy involved, it can be concluded that crack propagation is stable, regardless of the CZ model. Furthermore, as the tangential dissipative effects remain negligible, it is then possible to estimate $G_{Ic}$, an energy release rate in mode I, using a rule of areas (global estimate) based on a finite difference of Eq. (22) computed between times $t_i$ and $t_{i+1}$.

$$G_{Ic}^{(g)} \approx \frac{F_{kn_i} U_{l_{i+1}} - F_{kn_{i+1}} U_{l_i}}{a_{ct_{i+1}} - a_{ct_i}}, \quad (24)$$

where $U_l$ is the load displacement. This global estimate can then be compared to the energy dissipated by a CZ when $u_n = u_c$. This dissipated energy is equal to the area under the separation curve. Using the analytical expression for secant stiffness (see Eq. (13)), the energy dissipated at the rupture in mode I of the cohesive zone can be written as follows (local estimate):

$$G_{Ic}^{(l)} = \frac{K_n^0 r_n}{2(r_n + 2)} \quad (25)$$
Table 2 gives the different values involved in these local and global estimates of $G_{Ic}$ for the 3 CZMs.

<table><caption>Table 2: Local and global estimates of $G_{Ic}$</caption><thead><tr><th>CZM</th><th>t<sub>i</sub> (µs)</th><th>t<sub>i+1</sub> (µs)</th><th>a<sub>cr<sub>i</sub></sub> (µm)</th><th>a<sub>cr<sub>i+1</sub></sub> (µm)</th><th>F<sub>kn<sub>i</sub></sub> (N m<sup>-1</sup>)</th><th>F<sub>kt<sub>i</sub></sub> (N m<sup>-1</sup>)</th><th>U<sub>l<sub>i</sub></sub> (nm)</th><th>U<sub>l<sub>i+1</sub></sub> (nm)</th><th>G<sup>(g)</sup><sub>Ic</sub> (N m<sup>-2</sup>)</th><th>G<sup>(l)</sup><sub>Ic</sub> (N m<sup>-2</sup>)</th></tr></thead><tbody><tr><td>1</td><td>0.4</td><td>1.0</td><td>5.9</td><td>12.2</td><td>1504</td><td>578</td><td>1.6</td><td>8.2</td><td>1.8</td><td>1.7</td></tr><tr><td>2</td><td>0.4</td><td>1.0</td><td>5.9</td><td>11.9</td><td>1685</td><td>642</td><td>1.7</td><td>7.4</td><td>1.9</td><td>1.7</td></tr><tr><td>3</td><td>0.4</td><td>1.0</td><td>6.1</td><td>12.1</td><td>1497</td><td>578</td><td>1.6</td><td>10</td><td>2.3</td><td>2.4</td></tr></tbody></table>

Since Eq. (24) is a finite difference approximation of Eq. (22) when kinetic energy is neglected, the comparison of local and global estimates of $G_{Ic}$ is deemed to be satisfactory.

## 6.3. ENF tests

The results of the ENF tests for the three cohesive zones are gathered in FIGURE 12. The first thing that can be seen in the upper curves (a-c) is the progression of the process zone, even though the crack itself is not progressing. In the middle of the test, the crack suddenly progresses, causing a sharp decrease in force, followed by rapid oscillations around a slowly increasing average force value. During these oscillations, the process zone and the crack now progress together. These progressions are likely facilitated by the kinetic energy produced in the medium during the previous sudden advance of the crack.

The energy balances associated with cohesive zones correspond to the curves in the middle (FIGURE 12(d-f)). In contrast to the findings of the DCB tests, the crack does not initiate and the strained CZs, which are linked to the development of the process zone, multiply and become elastically loaded. Throughout this process, the elastic energy of the CZs remains very close to the work of the cohesive forces. This elastic energy reaches a maximum when the most strained CZs begin to dissipate, thereby reducing their elastic energy. This trend contin-
<figure><img src="image_16.png" /><figcaption>Figure 12: ENF test results for the 3 CZMs.</figcaption></figure>
ues as the crack progresses rapidly. At the end of the test, small oscillations in
deformation energy corresponding to structural vibrations can be seen for the 3
CZMs. As expected for this type of test, 'normal' dissipation is zero. The 'total'
and 'tangential' dissipation curves merge as soon as the process zone is triggered.

The bottom curves (FIGURE 12(g-i)) are associated with overall energy bal-
ances. This time, kinetic energy K becomes apparent as soon as the crack prop-
agates in an unstable manner. This quasi-instantaneous propagation can be seen
on the curve representing the evolution of the work of external forces Wext in the
form of a plateau. Initially, the structure loads elastically, the work of external
forces equals the elastic energy of the structure; when the process zone progresses,
some CZs start dissipating and the elastic energy decreases. At the end of the test,
the structure's very strict embedding conditions cause the external forces and the
structure's elastic energy to increase again. Meanwhile, in terms of cohesive zones,
the elastic energy gradually approaches zero, which indicates that all strained Czs
have failed completely.

6.4. MMB tests

As previously mentioned, the pre-crack in the MBB tests was lengthened from
0.5 mm to 2.5 mm. This was done to achieve relatively stable crack propagation.
Tests carried out with the pre-crack of 0.5 mm resulted in a sudden fracture for
the three CZMs, with the process zone almost exactly corresponding to the length
of the sample.

The results of the MMB tests are shown in FIGURE 13. For the three CZMs, the
upper figures (see FIGURE 13(a-c)) show that the process zone develops without
the crack propagating, as seen in ENF tests. The crack then progresses more
steadily. The load $F_{kn}$ increases at the start of the test, then decreases as the
structure becomes damaged. At the end of the tests, the load increased again. It is hypothesised that this is due to the embedding conditions at the left end of the structure when it is almost fully delaminated.

The energy balance at the cohesive zones (Fig.13(d-f)) is now standard. At the start of the test, during propagation of the process zone and without cracking, the cohesive zones are elastically loaded. The work of the cohesive forces is therefore equal to the stored elastic energy. Then, certain cohesive zones begin to dissipate energy and deform until they reach their displacement limit, which induces crack progression.

The uniqueness of MMB tests is linked to the appearance of modes I and II of decohesion. The graphs (Fig.13(d-f)) show the 'normal' and 'tangential' dissipations associated with these two modes, respectively.
<figure><img src="image_17.png" /><figcaption>Figure 13: MMB tests' results for the 3 CZMs. Subfigures (a)–(i) correspond to different specimens.</figcaption></figure>

Initially, delamination occurs predominantly in mode II, whereby the process zone progresses without cracking. Once the crack progresses, mode I begins to reach the same dissipated energy levels.

The overall energy balances shown in graphs (Fig.13(g-i)) demonstrate that cracking is proceeding stably, with the kinetic energy becoming negligible again
by the end of the test. However, a small amount of kinetic energy appears at the moment of the sudden advance of the crack. But this kinetic energy is quickly reabsorbed, probably consumed by the more regular advance of the crack and the process zone [37]. At the beginning of the tests, the structure loads elastically ($W_{el} \approx W_{ext}$). This elastic energy decreases slightly as the process zone grows. However, as soon as the crack advances due to the boundary conditions and loading mode, the structure's elastic energy begins to increase again, despite the CZs' regular decrease in elastic energy.

## 7. Concluding comments

This work presents a cohesive zone model based on the concept of available elastic energy, whose evolution characterizes the damage state through a scalar variable, independently of the loading path, along with its implementation within the LMGc90 code. The coupling of damage modes is formulated based on the following energy hypothesis: the maximum elastic energy that a cohesive zone can store for a given damage state is a function, $w_e^d$, of its damage level. One key advantage of this approach is that it allows anisotropic degradation of the elastic properties of the cohesive zone to be represented using a single scalar state variable, $u_d$, while accommodating various interactions between failure modes.

The modelling work is carried out within a thermodynamic framework that incorporates the irreversible nature of damage into the definition of the damage state variable (cf. Eq. (3)). Since damage is assumed to be the only irreversible, totally dissipative mechanism, no dissipation potential is necessary.

A wide range of possible interactions between degradation modes was revealed by small-scale numerical tests carried out on a single cohesive zone. One numerical challenge is managing the implicit definition of the variable $u_d$ (see Eq. (11)). Its implicit nature necessitates numerical resolution schemes that can significantly
increase computation times. For this first attempt, we used a basic dichotomy method. Although this method is not very efficient, it ensures convergence once the root has been initially separated. However, it is clear that other methods will need to be considered in the near future, particularly if cohesive areas are spread across the entire mesh. This is the reason why the final section of the paper focuses on using cohesive zones for flexural delamination tests on structures where the CZs are distributed solely at the interface between two elastic layers. These numerical tests are based on traditional experimental flexural tests used to study delamination in modes I (DCB), II (ENF) and mixed modes I/II (MMB). Although using an embedding instead of standard simple roller supports probably distorts the results, particularly at the end of the test when the crack has propagated significantly, the numerical results nevertheless align with those from classic mode I (DCB), mode II (ENF) and mixed mode (MMB).

DCB tests remain stable, and estimates of energy release rates obtained using the area rule correspond well with the energy dissipated by the cohesive zone up to failure. Although the loading remains quasi-static, ENF tests show that cracks can propagate dynamically. This process involves transforming deformation energy into dissipation and kinetic energy. The latter probably accelerates the progression of cracks [37]. To ensure stable propagation of the crack over the end of the test and to make modes I and/or II observable, the pre-crack is significantly increased for the ENF and MMB tests.

In the present study, the crack trajectory is prescribed by introducing cohesive zones along an initially straight interface and by using a mesh aligned with this interface. In this configuration, the crack propagation was verified to be independent of the time discretization and the mesh refinement. The next step consists
in performing simulations in which cohesive zones are distributed throughout the mesh. How would different mesh topologies and time integration schemes affect the results? How would the crack trajectory evolve under such conditions? Would the crack remain single or could branching occur? Could dynamic instabilities still develop? Finally, how should an equivalent decohesion surface be defined within the framework of Griffith's criterion? Such questions reflect the limitations generally associated with the CZ damage approach. Indeed, if they are effective in representing localised fracture processes, cohesive zone models are generally considered to have several limits. They often require prior anticipation of crack paths, since cohesive interfaces must be inserted a priori or introduced through adaptive strategies, which may restrict their predictive capability when crack trajectories are not known in advance. Regarding the CZMs proposed, the damage variable is kinematic in nature, with its evolution linked to an elastic energy. The deformation process remains stable, as observed in DCB tests, even though softening behaviour is present. Whether the reorganization of the deformation and stress fields induced by the inception of a crack will be sufficient to impose realistic crack paths remains to be seen.

Furthermore, cohesive zone models are often considered poorly suited to the description of diffuse damage, as they do not naturally account for the progressive degradation of materials that may occur prior to crack initiation. Their predictions may also be sensitive to mesh discretization and interface placement, which can lead to biased results if cohesive surfaces are not appropriately positioned.

This limitation partly arises from the fact that many CZMs introduce their own elastic (or even rigid) domain. In the CZMs considered in the present work, damage initiates as soon as deformation begins. The ENF and MMB simulations reveal that the process zone can develop spatially without any visible crack propagation.

From an experimental standpoint, infrared thermography studies [38] have
704 shown that, even during nominally “elastic” fatigue cycles [39], dissipation fields
705 are highly heterogeneous from the very first loading cycles. Moreover, regions of
706 high dissipation are not necessarily located where macroscopic stress levels are
707 expected to be maximal [40]. Such heterogeneous dissipation patterns can be at-
708 tributed to an initial microstructural state that is far more complex than commonly
709 assumed, challenging the notion of a spatially uniform initial material state.

710 Part of the material response therefore appears to be predetermined from the
711 onset of loading, which may explain the occurrence of diffuse damage. In this
712 context, the forthcoming numerical investigations, based on distributing cohesive
713 zones throughout the mesh, are expected to provide further insight into the ques-
714 tions raised in this study.

8. Acknowledgements

716 The authors would like to sincerely thank F. Dubois and R. Mozul, from the
717 COMPEX service, for their assistance in getting started with the IMc90 software
718 and for their role in ensuring the development and the dissemination of the code.

719 References

[1] L. Kachanov, Rupture time under creep conditions, Izv. Akad. Nauk SSSR 8 (1958) 26-31.

[2] L. Kachanov, Introduction to continuum damage mechanics, Vol. 10, Springer Science & Business Media, 2013.

[3] J. Lemaitre, A continuous damage mechanics model for ductile fracture, Journal of engineering materials and technology 107 (1) (1985) 83-89.
726

727

[4] J. Lemaitre, A course on damage mechanics, Springer science & business media, 2012.

728

729

[5] J.-L. Chaboche, Continuum damage mechanics: Part i—general concepts, Journal of Applied Mechanics 55 (1988) 59–64.

730

731

[6] G. Pijaudier-Cabot, Z. P. Bažant, Nonlocal damage theory, Journal of engineering mechanics 113 (10) (1987) 1512–1533.

732

733

734

735

[7] R. H. Peerlings, R. de Borst, W. Brekelmans, M. G. Geers, Gradient-enhanced damage modelling of concrete fracture, Mechanics of Cohesive-frictional Materials: An International Journal on Experiments, Modelling and Computation of Materials and Structures 3 (4) (1998) 323–342.

736

737

[8] M. Frémond, M. Shitikova, Non-smooth thermomechanics, Appl. Mech. Rev. 55 (5) (2002) B99–B100.

738

739

740

[9] G. A. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, Journal of the Mechanics and Physics of Solids 46 (8) (1998) 1319–1342.

741

742

743

744

[10] A. A. Griffith, Vi. the phenomena of rupture and flow in solids, Philosophical Transactions of the Royal Society of London, Series A: Containing Papers of a Mathematical or Physical Character 221 (582-593) (1921) 163–198.
doi:10.1098/rsta.1921.0006.

745

746

[11] B. Bourdin, G. A. Francfort, J.-J. Marigo, The variational approach to fracture, Journal of elasticity 91 (1) (2008) 5–148.

747

748

749

[12] D. Dugdale, Yielding of steel sheets containing slits, Journal of the Mechanics and Physics of Solids 8 (2) (1960) 100 – 104.
doi:https://doi.org/10.1016/0022-5096(60)90013-2.
[13] G. I. Barenblatt, The mathematical theory of equilibrium cracks in brittle fracture, in: H. L. Dryden, T. von Kármán, G. Kuerti, F. H. van den Dungen, L. Howarth (Eds.), Advances in Applied Mechanics, Vol. 7, Elsevier, 1962, pp. 55-129. doi:10.1016/S0065-2156(08)70121-2.

[14] J. R. Rice, A path independent integral and the approximate analysis of strain concentration by notches and cracks, Journal of Applied Mechanics 35 (2) (1968) 379-386. doi:10.1115/1.3601206.

[15] J. Willis, Fracture mechanics of interfacial cracks, Journal of the Mechanics and Physics of Solids 19 (6) (1971) 353-368. doi:https://doi.org/10.1016/0022-5096(71)90004-4.
URL https://www.sciencedirect.com/science/article/pii/0022509671900044

[16] A. C. Palmer, J. R. Rice, The growth of slip surfaces in the progressive failure of over-consolidated clay, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 332 (1591) (1973) 527-548. doi:10.1098/rspa.1973.0040.

[17] M. Raous, L. Cangémi, M. Cocu, A consistent model coupling adhesion, friction, and unilateral contact, Computer Methods in Applied Mechanics and Engineering 177 (3) (1999) 383-399. doi:https://doi.org/10.1016/S0045-7825(98)00389-2.

[18] Y. Monerie, V. Acary, Formulation dynamique d'un modèle de zone cohésive tridimensionnel couplant endommagement et frottement, Revue européenne des éléments finis 10 (2-4) (2001) 489-503.

[19] L. Daridon, B. Cochelin, M. Potier-Ferry, Delamination and fiber bridging
773 modelling in composite samples, Journal of Composite Materials 31 (9) (1997) 874-888.
774

[20] M. E. Gurtin, Thermodynamics and the cohesive zone in fracture, Journal of Applied Mathematics and Physics (ZAMP) 30 (6) (1979) 991-1003.
775
776

[21] F. Costanzo, D. Allen, A continuum thermodynamic analysis of cohesive zone models, International Journal of Engineering Science 33 (15) (1995) 2197-2219, the Edelen Symposium. doi:https://doi.org/10.1016/0020-7225(95)00066-7.
777
778
779
780

[22] B. Halphen, Q.-S. Nguyen, Sur les matériaux standard generalisés, Journal de mécanique 14 (1975) 39,63.
781
782

[23] A. Chrysochoos, L. Daridon, M. Renouf, A damage criterion based on energy balance for isotropic cohesive zone model, Journal of Theoretical, Computational and Applied Mechanics (Mar 2022). doi:10.46298/jtcam.7056.
783
784
785

[24] M. Ortiz, A. Pandolfi, Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis, International journal for numerical methods in engineering 44 (9) (1999) 1267-1282.
786
787
788

[25] K. Roe, T. Siegmund, An irreversible cohesive zone model for interface fatigue crack growth simulation, Engineering fracture mechanics 70 (2) (2003) 209-232.
789
790
791

[26] F. Evangelista Jr, J. R. Roesler, S. P. Proença, Three-dimensional cohesive zone model for fracture of cementitious materials based on the thermodynamics of irreversible processes, Engineering Fracture Mechanics 97 (2013) 261-280.
792
793
794
795
796

797

798

[27] R. Serpieri, E. Sacco, G. Alfano, 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 (2015) 13–25.

799

800

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

801

802

803

[29] L. Daridon, B. Wattrisse, A. Chrysochoos, M. Potier-Ferry, Solving fracture problems using an asymptotic numerical method, Computers & Structures 89 (5) (2011) 476 – 484.

804

805

806

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

807

808

809

[31] C. Shet, N. Chandra, Analysis of energy balance when using cohesive zone models to simulate fracture processes, J. Eng. Mater. Technol. 124 (4) (2002) 440–450.

810

811

[32] J. R. Reeder, J. H. Crews, Mixed-mode bending method for delamination testing, AIAA Journal 28 (7) (1990) 1270–1276. doi:10.2514/3.25204.

812

813

814

[33] J. J. Moreau, Unilateral contact and dry friction in finite freedom dynamics, in: J. J. Moreau, P. D. Panagiotopoulos (Eds.), Nonsmooth Mechanics and Applications, Springer, Vienna, 1988, pp. 1–82.

815

816

817

[34] M. Jean, The non-smooth contact dynamics method, Computer Methods in Applied Mechanics and Engineering 177 (3) (1999) 235 – 257. doi:https://doi.org/10.1016/S0045-7825(98)00383-1.

818

819

[35] F. Dubois, V. Acary, M. Jean, The contact dynamics method: A nonsmooth story, Comptes Rendus Mécanique 346 (3) (2018) 247–262.
[36] M. Champagne, M. Renouf, Y. Berthier, Modeling Wear for Heterogeneous Bi-Phasic Materials Using Discrete Elements Approach, Journal of Tribology 136 (2), 021603 (01 2014). doi:10.1115/1.4026053.

[37] H. BUI, Brittle fracture mechanics(french book), Paris, Masson, Editeur, 1978. 231 p, In French (1978).

[38] N. Ranc, A. Blanche, D. Ryckelynck, A. Chrysochoos, Pod preprocessing of ir thermal data to assess heat source distributions, Experimental Mechanics 55 (4) (2015) 725-739.

[39] B. Berthel, B. Wattrisse, A. Chrysochoos, A. Galtier, Thermographic analysis of fatigue dissipation properties of steel sheets, Strain 43 (3) (2007) 273-279.

[40] A. Blanche, A. Chrysochoos, N. Ranc, V. Favier, Dissipation assessments during dynamic very high cycle fatigue tests, experimental mechanics 55 (4) (2015) 699-709.