Engineering Fracture Mechanics 73 (2006) 2774–2801

[Link]/locate/engfracmech

A damage-mechanics-based approach for modelling

decohesion in adhesively bonded assemblies
a,* b
Nunziante Valoroso , Laurent Champaney
a
Istituto per le Tecnologie della Costruzione, Consiglio Nazionale delle Ricerche, Viale Palmiro Togliatti 1473,
Roma 00155, Italy
b
Laboratoire d’Etudes Mécaniques des Assemblages, Université de Versailles Saint-Quentin-en-Yvelines,
45 avenue des Etats-Unis, Versailles Cedex 78035, France

Received 31 March 2005; received in revised form 8 April 2006; accepted 28 April 2006
Available online 10 July 2006

Abstract

A cohesive interface model formulated within the framework of damage mechanics is presented for the simulation of
decohesion in adhesively bonded assemblies. Characteristic features of the model are: the introduction of a single energy-
based damage variable for describing the damage state of the interface; use of a decohesion propagation condition relying
upon the linear elastic fracture mechanics (LEFM) energy balance; a treatment for the mixed-mode situation based on the
definition of an equivalent energy release rate whose expression is consistently derived from the formulation. The compari-
sons between numerical and experimental responses obtained for typical test problems illustrate the capabilities of the pro-
posed approach.
Ó 2006 Elsevier Ltd. All rights reserved.

Keywords: Damage mechanics; Adhesive assemblies; Interfaces; Cohesive-zone models; Decohesion

1. Introduction

Adhesive bonding represents a fascinating technique to join members of a structure as it can offer enhanced
flexibility in design and lower maintenance costs compared to mechanical fastening; moreover, especially for
light-weight composite structures, the use of structural adhesives is highly desirable since unfavourable stress
concentrations due to the presence of bolted or riveted joints can be reduced. As for most structural compo-
nents consisting of the assembly of individual elements, failure of adhesive joints due to damage nucleation
and growth at bonded interfaces is a common occurrence; for its simulation cohesive interface modelling
can be usefully resorted to.
The cohesive-zone concept, originally introduced by Barenblatt [1,2] and Dugdale [3] to describe the
near-tip fracture process, and subsequently put forward by Hillerborg [4,5] in the fictitious crack model,

*
Corresponding author. Tel.: +39 06 7259 7014/40900 7345; fax: +39 06 7259 7005/40900 7367.
E-mail addresses: [Link]@[Link] (N. Valoroso), [Link]@[Link] (L. Champaney).

0013-7944/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved.
doi:10.1016/[Link].2006.04.029
 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2775

Nomenclature

sub displacement jump (single mode)

t interface traction (single mode)
k+, k interface stiffnesses in tension/compression (single mode)
D damage variable
Y damage-driving force (single mode)
Yf characteristic value of Y (single mode)
Y* critical damage-driving force (single mode)
Go elastic energy at damage onset (single mode)
Gc critical fracture energy (single mode)
F damage function (single mode)
FA, FB, FC damage functions for models A, B, C
sub displacement jump vector
sunb, susb displacement jump, normal/sliding components
t interface traction vector
tn, ts interface traction, normal/sliding components
kþn ; k 
n interface stiffness in tension/compression, normal component
ks interface stiffness, sliding (tangential) component
Ym damage-driving force (mixed-mode)
YI, YII damage-driving force, pure mode components
Ymf characteristic value of Ym
YmfA, YmfB, YmfC characteristic values of Ym for models A, B, C
Ymo damage-driving force at damage onset (mixed-mode)
Y m critical damage-driving force (mixed-mode)
GoI,GoII elastic energy at damage onset, pure mode components
GI, GII released energy, pure mode components
GT total released energy (=GI + GII)
GcI, GcII critical fracture energy, pure mode components
GTc total work of separation
u, b loading angle, mode-mixity parameter
d, td equivalent opening displacement, equivalent traction
Fm damage function (mixed-mode)
FmA, FmB, FmC mixed-mode damage functions for models A, B, C

has gained major popularity in recent years for simulating delamination, debonding, fracture and fragmen-
tation via finite element methods [6–31], see also [32] for a survey account. This approach is mainly
motivated by the consideration that prior to the development of macroscopic fractures there exists a
zone in a state of progressive damage located in front of the crack, the so-called cohesive process zone,
where an interaction across the two sides of the crack itself in terms of non-zero cohesive tractions is still
present.
In most current implementations, cohesive interfaces are used in conjunction with zero-thickness decohe-
sion or interface elements in which the process zone is lumped; these elements are located at sites where the
potential crack trajectories can develop, e.g. where elastic mismatches occur in the structure or a layered mate-
rial is present, and are equipped with a non-linear traction–relative displacement relationship describing the
process zone evolution and the formation of traction-free surfaces.
Differently from other methods that are directly based on fracture mechanics, cohesive-zone models can
quite naturally account for the variation of the experimentally measured fracture toughness as function of
mode mixity under general loading conditions, see e.g. [9,15,18,25,28,30]. More generally, they provide an
2776 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

extremely flexible tool for structural analysis, the only real limitation being represented by the need for a suf-
ficiently refined mesh in the vicinity of the process zone in order to follow the strongly non-linear structural
response, which might be difficult to track even by using sophisticated path-following techniques [33].
In this work we present an interface model accounting for progressive irreversible damage prior to complete
separation; the proposed formulation, that allows to recover several earlier models as special cases, relies upon
the introduction of a single energy-based damage variable whose evolution is related to that of a critical dam-
age-driving force. The construction of the different damage models is detailed for the one-dimensional (single-
mode) situation, in which case the propagation of decohesion is assumed to take place when the dissipated
energy equals the critical fracture energy. The mixed-mode case is dealt with by making reference to a
mode-mixity parameter and an equivalent energy release rate; unlike other models, this last one is not intro-
duced at the outset but is consistently derived from the formulation. Damage onset and decohesion propaga-
tion conditions are then obtained based upon two interaction criteria taken from the composites literature.
The outline of the paper is as follows. In Section 2.1 the one-dimensional model is first introduced. Its
extension to the mixed-mode case is discussed in Section 2.2 and the derivation of the material tangent oper-
ator is provided in Section 2.3. Numerical simulations and comparisons with experimental results relevant to
aluminum components bonded by epoxy resin are presented in Section 3. Summary and conclusions are finally
given in Section 4.

2. Interface model

As a model problem reference is made to an assembly consisting of two elastic bodies (adherends) joined by
a plane adhesive layer whose thickness is assumed to be negligible compared to both that of the joined bodies
and to its in-plane dimensions.
These features enable the adhesive layer to be conveniently schematized as an interface, i.e. as a zero-thick-
ness surface entity ensuring stress transfer between the adherends, see Fig. 1.
A constitutive relationship for the interface is given in terms of the displacement jump vector sub = u+  u-,
which has the meaning of interface strain, and the related traction vector t, which plays the role of the stress.
For the problem at hand the interface constitution will be considered as the only source of material non-
linearity; all expected dissipative phenomena are thus supposed to take place at this level only.

2.1. One-dimensional model

To begin with, we start by considering the one-dimensional (mode I) case, for which the stored energy func-
tion is introduced as

1 2 1 2
wðsut; DÞ ¼ ð1  DÞk þ hsutiþ þ k  hsuti ð1Þ
2 2

where D 2 [0,1] denotes a scalar damage variable in the usual sense [34], the symbols h Æ i+ and h Æ i stand for
the positive and negative part of the argument h Æ i, defined as hxi± = 1/2(x ± jxj), and k+ and k are the
undamaged interface stiffnesses in tension and compression, respectively, the latter representing a penalty stiff-
ness accounting for the impenetrability constraint.

Fig. 1. Interface schematization.

 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2777

In the classical approach, also the damageable stiffness k+ is usually taken as a penalty stiffness parameter
that is used to ensure that interface strains prior to damage onset are sufficiently small [25,30]. Having in mind
the analysis of an adhesive joint, in this study we shall however consider the damageable interface stiffnesses as
model parameters since they can be effectively estimated via guided waves propagating in the adhesive layer,
see e.g. [35–37] for a detailed discussion.
The constitutive equations for the model follow from standard continuum thermodynamics and from the
expression of mechanical dissipation as
ow
t¼ ¼ ð1  DÞk þ hsutiþ þ k  hsuti
osut
ð2Þ
ow 1 þ 2
Y ¼ ¼ k hsutiþ
oD 2
and identify the interface traction t and the damage energy release rate Y as the work-conjugates of the dis-
placement jump and of the damage variable, respectively.
As in classical internal variables theories, the damage-driving force is assumed to be bounded by a critical
value. This can be characterized by means of a criterion of the form:
/ ¼ Y  Y 6 0 ð3Þ
*
where Y denotes the critical damage-driving force, whose meaning is that of an instantaneous energy thresh-
old; prior to the application of any loading it is assumed Y* = Go P 0, that is not necessarily the elastic energy
corresponding to the peak of the local traction–separation relationship, but has to be intended in the more
general sense of the energy at the onset of damage.
The increase in the size of the elastic locus defined by (3) is associated with growth of damage and, as such,
is irreversible. This can be accounted for by specifying the evolution equations as
o/ oF
D_ ¼ c_ ; Y_  ¼ c_ ð4Þ
oY oD
where c_ is a Lagrange multiplier subject to the Karush–Kuhn–Tucker (KKT) conditions [38]:
/ 6 0; c_ P 0; c_ / ¼ 0 ð5Þ
and F is a monotonically increasing positive function of the argument D.
It is worth noting the interpretation of the above relationships for the situation of increasing damage
ð_c > 0Þ. In this case, using (5) along with the condition of damage persistency:
Y_  ¼ Y_ ð6Þ
one has
 

Y ðtÞ ¼ max Go ; maxfY ðsÞg ð7Þ
ðs6tÞ

t 2 [0, T] being the pseudo-time parameter.

The damage model can be completely characterized by the choice of the relationship between the energy
threshold Y* and the damage variable D. In this respect, a basic condition to be accounted for is the energy
criterion of fracture mechanics, i.e. the requirement that the energy dissipated during the decohesion process
equals the critical fracture energy Gc:
Z þ1
Y D_ dt ¼ Gc ð8Þ
0

The simplest non-smooth model, defined as

(
Gc if D 2 ½0; 1½

Y ¼ max Y ðsÞ if D ¼ 1 ð9Þ
s2½0;T 

trivially satisfies the criterion (8), see also Fig. 2 for a pictorial representation.
2778 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

traction t

Gc

displacement jump [u]

Fig. 2. Non-smooth model.

In essence, the cohesive-zone approach amounts to a regularization of Eq. (9), that are equivalent to Grif-
fith’s criterion of LEFM. When recast into a damage mechanics framework, a cohesive model can be given the
form:
8
> Go if D ¼ 0
>
>
<Rt 
Y ¼ 0
Y_ dt ¼ F ðDÞ if D 2 0; 1½ ð10Þ
>
>
>
: max Y ðsÞ if D ¼ 1
s2½0;T 

Typical forms of F are that of power laws or exponential functions, and their explicit expressions are con-
structed in a way to fulfill (8). The standard procedure, see e.g. [10,13,17,23,28,30] among others, amounts
to introducing, more or less explicitly, in the expression of the function F, a characteristic value Yf of the dam-
age-driving force that may well be the one corresponding to complete decohesion, and by computing it by
enforcing the condition (8).
In the following three different damage models, yielding a bilinear (A), polynomial (B) and exponential (C)
traction–separation relationship will be considered. The corresponding expressions obtained for the damage
function F are denoted by appending to it the relevant subscript.
As shown in [30], the widely used bilinear traction–relative displacement relationship, see Fig. 3, can be
obtained by taking:
8
>
> Go if D ¼ 0
>
< Y 2f Go
Y  ¼ ½Go DþY f ð1DÞ2 ¼ F A ðDÞ if D 2 0; 1½ ð11Þ
>
>
>
: max Y ðsÞ if D ¼ 1
s2½0;T 

whereby the characteristic value Yf of the damage energy release rate, that in this case does not correspond to
the completely damaged state, is found from (8) as
Y f ¼ Gc ð12Þ
so to obtain
 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2779

16

14

12

10
traction t

0
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
displacement jump [u]

Fig. 3. Traction–separation relationship for the bilinear model, see (11). k+ = 1e4, Gc = 0.125, g = 0.91.

ð1  gÞGc
F A ðDÞ ¼ ð13Þ
ð1  gDÞ2
where
Go
g¼1 ð14Þ
Gc
plays the role of a regularizing parameter.
Actually, for given (undamaged) elastic stiffness, as g ! 0 one has Go ! Gc, thus recovering a perfectly brit-
tle behaviour, whereas for g ! 1 it turns out to be Go ! 0, that corresponds to vanishing strength and infinite
displacements at complete decohesion. Obviously, the case g = 1 has to be ruled out when considering a linear
softening relationship, since in this case the regularized damage law (11) is not a workable one.
This limit situation can however be handled when considering more sophisticated damage models, in which
more than one regularizing parameter can be defined. This is the case of the Allix & Ladèveze polynomial
model, in particular the version implemented in [13,19,23,24,29], that can be obtained by defining the critical
damage-driving force as
8
>
> Go if D ¼ 0
< 1=N
Y  ¼ Go þ ðY f  Go ÞD ¼ F B ðDÞ if D 2 0; 1½ ð15Þ
>
>
: max Y ðsÞ if D ¼ 1
s2½0;T 

for N > 0 and not necessarily integer.

The characteristic value of the force Yf corresponding to complete decohesion is computed from the energy
balance (8) as
Z 1
N
Y  ðDÞ dD ¼ Go þ ðY f  Go Þ ¼ Gc ð16Þ
0 N þ1
Accordingly, one obtains:
 
N þ1
F B ðDÞ ¼ ð1  gÞGc þ gGc D1=N ð17Þ
N
2780 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

16 50

45
14

40
12
35

10
30
traction t

traction t
8 25
N
20
6

15
4
10

2
5

0 0
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
displacement jump [u] displacement jump [u]

(a) (b)

Fig. 4. Traction–separation relationships for the Allix & Ladevèze polynomial model, see (15). k+ = 1e4, Gc = 0.125, g = 1.0; N = 0.2 (a)
and variable N (b).

where, besides the parameter g defined by (14), also the exponent N acts as a regularizing parameter in that a
perfectly brittle behaviour can be recovered, independently from the value of g, in the limit as N ! +1. The
traction–separation behaviours corresponding to (17) are depicted in Fig. 4.
Following the same path of reasoning, other damage laws endowed with two regularizing parameters can
be defined. For instance, an exponential damage model can be obtained by defining the critical energy Y* as
8
>
> Go if D ¼ 0
< N

Y ¼ G o þ ðY f  G o Þ½ logð1  DÞ ¼ F C ðDÞ if D 2 0; 1½ ð18Þ
>
>
: max Y ðsÞ if D ¼ 1
s2½0;T 

For this model the energy dissipated in the decohesion process is obtained from (8) as
Z 1
Y  ðDÞdD ¼ Go þ ðY f  Go ÞCðN þ 1Þ ¼ Gc ð19Þ
0

where C is the Gamma function [39], defined by

Z þ1
CðN þ 1Þ ¼ xN ex dx ¼ N  CðN Þ ð20Þ
0

i.e. the function extending the definition of the factorial to the whole complex plane, except for negative inte-
gers. There are a variety of methods in use for calculating the function C(N) numerically; for positive argu-
ments one of the most accurate is the approximation derived by Lanczos [40].
For the above model the damage function is then computed as
1
F C ðDÞ ¼ ð1  gÞGc þ gGc ½ logð1  DÞN ð21Þ
CðN þ 1Þ
which shows that, analogous to the previous case, allows a non-smooth, perfectly brittle behaviour to be
recovered, independently from g, in the limit as N ! 0, see also Fig. 5.
Remark 2.1. In the present context the computation of the damage state does not require the time-
discretization of the evolution equation (4) and the damage update only amounts to function evaluation.
Actually, due to (5) and (6), subsequent to damage loading, characterized by a positive trial value of the
damage mode (3), i.e.:
 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2781

16 50

45
14

40
12
35

10
30
traction t

traction t
8 25
N
20
6

15
4
10

2
5

0 0
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
displacement jump [u] displacement jump [u]

(a) (b)
+
Fig. 5. Traction–separation relationships for the exponential model, see (18). k = 1e4, Gc = 0.125, g = 1.0; N = 1.7 (a) and variable
N (b).

/trnþ1 ¼ Y ðtnþ1 Þ  Y  ðtn Þ > 0 ) c>0 ð22Þ

one does not need to evaluate the discrete damage multiplier c since in this case the value of the instantaneous
energy threshold Y*, which is a given function of the damage variable D, has to equal that of the actual energy
release rate Y.
Hence, at each time t the damage variable can be computed from
 
1 
DðtÞ ¼ min 1; maxfF ðY ðsÞÞg ð23Þ
ðs6tÞ

because of the positiveness and monotonicity of the damage function F.

Remark 2.2. It is worth noting that, for certain choices of the model parameters g and N, the damage func-
tions FB and FC, see (17) and (21), can lead to traction–separation relationships that may possess substantially
different shapes with respect to those reported in Figs. 4 and 5.
In particular, for N < 1 (N > 1), since the derivative of the damage function FB (FC) with respect to the
damage parameter vanishes in the limit for D ! 0+, is not difficult to show that for g < 1 the material tangent
stiffness turns out to be infinite at the onset of damage. The resulting traction–separation relationships that are
obtained in this case are given in Fig. 6.
On the contrary, for g = 1 the shape of the cohesive relationships given thus far remain unchanged because
the derivative of the damage variable with respect to the displacement jump is zero at damage onset, since this
occurs for sub = 0.

Remark 2.3. Owing to the softening behaviour of the local traction–separation relationship, the structural
response will unavoidably suffer from some mesh-dependence. However, as opposite to continuum damage
models, for which mesh-dependence occurs due to strain localization, in a structural model containing dam-
aging interfaces the mesh-dependent character of the solution can manifest in the form of response bifurcation
and loss of overall stability. A simple though striking example of this pathologic behaviour is discussed in [41].
This problem has been reported by many analysts, see e.g. [30,42] among others, and the general
recommendation is that of using a sufficiently refined mesh around the decohesion front in order to allow
values of the peak stress of the local traction–separation relationship that are high enough to correctly predict
the decohesion.
2782 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

16 16

14 14

12 12

10 10
traction t

traction t
8 8

6 6

4 4

2 2

0 0
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
displacement jump [u] displacement jump [u]
(a) (b)
+
Fig. 6. Traction–separation relationships for the polynomial (a) and exponential (b) models. k = 1e4, Gc = 0.125, g = 0.91; N = 0.1 (a)
and N = 3.5 (b).

A possible remedy for alleviating the cited shortcoming is that of using a viscous-like regularization, so to
make the interface model able to react to a potentially sudden evolution of the damage state by an increase in
the cohesive tractions.
This may be accomplished in different ways. For example, following the Duvaut–Lions approach [43], one
can define the damage rate as
1
D_ ¼  ðD  DÞ ð24Þ
l
where l is a relaxation time parameter and D is the rate-independent solution, whose computation is carried
out in the standard way.
For given D the integration of (24) is straightforward and the solution for the rate-dependent damage
variable reads
Dðtnþ1 Þ ¼ #Dðtn Þ þ ð1  #ÞDðtnþ1 Þ ð25Þ
with # = exp(Dt/l).
The other classical way of introducing rate-dependency is a Perzyna-like penalization [44], see also [45].
Accordingly, one has the constitutive relation for the damage variable:
1
D_ ¼ Kð/ðY ÞÞ ð26Þ
l
K being the flow function. Other examples of viscous-like regularization are discussed in [19,29,31].
As already pointed out in [41], the introduction of rate effects requires a special care since for highly rate-
sensitive models the energy dissipated by viscosity may be several times greater than the one dissipated in the
damage process.

2.2. Two-dimensional model

As opposite to the one-dimensional (single-mode) case, where the criteria used for determining damage
onset and propagation up to complete failure only involve one single component of the energy release rate,
when considering mixed-mode conditions the total energy released during decohesion:
GT ¼ GI þ GII ð27Þ
 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2783

results from the interplay of the I and II pure-mode contributions that evolve together as a consequence of the
interaction between the traction–displacement jump relationships in the normal and tangential directions.
Basically, two ways appear to be suggested in the literature to describe the mixed-mode situation, either
based upon the definition of an equivalent energy release rate, see e.g. [18,23,28], or relying on the introduction
of an equivalent displacement parameter [20,27] or a non-dimensional measure of separation [11,15,30].
In the following we shall appeal to some equivalent quantities for the formulation of the mixed-mode
model. We underline that, unlike other models, the equivalent variables we shall make use of will not be intro-
duced at the outset but, rather, will be derived in the damage mechanics framework outlined in the previous
section for the one-dimensional interface model; in the authors’ opinion, this makes their meaning more
apparent.
Following these ideas, we start by considering the stored energy function
1 2 2 1  2
wðsut; DÞ ¼ ð1  DÞ½k þ n hsun tiþ þ k s sus t  þ k n hsun ti ð28Þ
2 2
where sunb and susb denote the normal and sliding components of the displacement jump vector sub, i.e.
sunb = sub Æ n, susb = sub Æ s, n and s being the outward unit normal and the unit tangent vector to the interface,
see also Fig. 1.
The constitutive equations for the interface traction vector and the work-conjugate of the damage variable
are obtained via the classical thermodynamic argument as
ow
t ¼ osut ¼ ð1  DÞ½k þ 
n hsun tiþ n þ k s sus ts þ k n hsun ti n
2 2
ð29Þ
ow
Y m ¼  oD ¼ 12 k þ 1
n hsun tiþ þ 2 k s sus t

In the above equation and in the remainder of the paper the subscript m is appended to the mixed-mode vari-
ables in order to emphasize the difference with the analogous unsuffixed variables, that refer to the single-mode
case.
Based on the above relationships, the equivalent, mixed-mode energy release rate Ym can be expressed as
1
Y m ¼ kþ d2 ð30Þ
2 n
where d is an equivalent opening displacement given by
2 2 1=2
d ¼ ðhsun tiþ þ a2 sus t Þ ð31Þ
being
sffiffiffiffiffi
ks
a¼ ð32Þ
kþn

the parameter giving different weights to the normal and sliding components depending on the relevant initial
elastic stiffnesses.
Introducing the loading angle:
 
sus t
u ¼ arctan 2 ½0; þp=2 ð33Þ
hsun tiþ
a mixed-mode parameter b can be defined as
b ¼ a tanðuÞ ð34Þ
whereby the expressions of the pure-mode contributions to (30) follow as
1
YI ¼ Ym
1 þ b2
ð35Þ
b2
Y II ¼ Ym
1 þ b2
2784 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

that respectively represent the mode I (opening) and mode II (sliding) shares. We note in passing that the
above definition of the mixed-mode parameter is closely related to that of the phase angle given by Suo
and Hutchinson [46].
By ruling out the penalty term accounting for non-interpenetration, the cohesive relationship can thus be
reformulated as
td ¼ ð1  DÞk þ
nd ð36Þ
where td is an equivalent scalar traction:
 1=2
1
td ¼ t2n þ 2 t2s ð37Þ
a
being
1 cosðuÞ
tn ¼ ht  niþ ¼ t
2 1=2 d
¼ td
ð1 þ b Þ cos2 ðuÞþ a2 sin2 ðuÞ
ð38Þ
ab a2 sinðuÞ
ts ¼ t  s ¼ t d ¼ t d
ð1 þ b2 Þ1=2 cos2 ðuÞ þ a2 sin2 ðuÞ
the normal and sliding components of the traction vector, whose direction is in general different from that of
the resulting separation.
Having identified the damage-driving force as the mixed-mode energy release rate (30), one can specify the
evolution equations as
o/ oF m
D_ ¼ c_ m ; Y_ m ¼ c_ ð39Þ
oY m oD
along with the KKT conditions:
/m 6 0; c_ P 0; c_ /m ¼ 0 ð40Þ
for the damage mode:
/m ¼ Y m  Y m ð41Þ
where, analogous to the one-dimensional case, Y m
denotes the mixed-mode instantaneous critical energy
release rate, whose evolution is governed by a monotonically increasing positive function Fm.
As opposite to the one-dimensional situation, where the damage onset is determined by comparing the
energy release rate with the initial pure-mode threshold Go, under mixed-mode loading damage can occur
before any single-mode component attains its initial allowable value. Accordingly, the definition of the critical
damage-driving force for a regularized mixed-mode model, that at least formally can be given as in the single-
mode case, i.e. as
8
>
> Y mo if D ¼ 0
<

Ym ¼ F m ðDÞ if D 2 0; 1½ ð42Þ
>
>
: max Y m ðsÞ if D ¼ 1
s2½0;T 

requires the determination of two parameters, Ymo and Ymf, based on suitable interaction criteria determining
the damage onset and propagation of decohesion.
Neglecting the effect of compressive normal tractions, it is assumed that the initiation of damage can be
predicted using a Hashin-type [47] criterion, that can be conveniently written as
 a1  a
YI Y II 2
þ 1¼0 ð43Þ
GoI GoII
where GoI and GoII are the initial pure-mode damage thresholds while a1, a2 are model parameters to be chosen
in accordance with experimental data, that are assumed to be both strictly positive and non-necessarily
integer.
 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2785

The initial mixed-mode threshold Ymo can be computed from the previous equation that, on account of
(35), can be expressed in the equivalent form:
a a
fo ðb; Y m Þ ¼ cI ðbÞðY m Þ 1 þ cII ðbÞðY m Þ 2  1 ¼ 0 ð44Þ
for
 a1  a2
1 b2
cI ðbÞ ¼ ; cII ðbÞ ¼ ð45Þ
ð1 þ b2 ÞGoI ð1 þ b2 ÞGoII
Simple inspection shows that the function fo has continuous partial derivatives and that it is a strictly increas-
ing function of the argument Ym. Therefore, Eq. (44) implicitly defines Ym as a function of the mode-mixity b
in the whole domain of interest; as a consequence, for a given b there is always a unique solution Ymo of (44)
and its computation only requires a local iterative scheme.
For a2 = a1 one has also the closed-form solution:
 
logðcI þ cII Þ ð1 þ b2 ÞGoI GoII
Y mo ¼ exp  ¼ ð46Þ
a1 a a 1=a1
ðGoII Þ 1 þ ðb2 GoI Þ 1
whence the pure-mode threshold energies GoI and GoII are recognized to be recovered in the limit as u ! 0+
and u ! p/2, respectively.
As for damage onset, under mixed-mode loading the dependence of the fracture thoughness on the mode
ratio has to be properly accounted for. In this respect the adopted failure condition, stemming from one of the
most widely used criteria to predict the propagation of delamination in composites, is a generalized ellipse-like
criterion [48]:
 b1  b
GI GII 2
þ 1¼0 ð47Þ
GcI GcII
where the exponents b1, b2 are strictly positive reals while the mode I and mode II released energies are given
by
Z þ1
Gi ¼ Y i D_ dt; i 2 fI; IIg ð48Þ
0

From (35) it follows that, for a radial loading path, the ratio between GII and GI is constant and equals b2.
Accordingly, from (47) one has the non-linear equation:
b b
fc ðb; GT Þ ¼ d I ðbÞðGT Þ 1 þ d II ðbÞðGT Þ 2  1 ¼ 0 ð49Þ
where GT is defined by (27) and dI, dII have the same expressions as in (45), with GoI, GoII, a1, a2 being replaced
by GcI, GcII, b1, b2. Clearly, all the considerations developed for (44) apply to (49) as well.
For the particular case of b2 = b1 the propagation of decohesion takes place for
ð1 þ b2 ÞGcI GcII
GTc ¼ ð50Þ
½ðGcII Þb1 þ ðb2 GcI Þb1 1=b1
where GTc is computed as the total work of separation:
Z þ1
GTc ¼ Y m D_ dt ð51Þ
0

whose expression depends upon that of the function Fm defining the critical damage-driving force in the range
D 2 ]0, 1[.
In particular, taking for Fm one of the expressions used in the one-dimensional case, i.e.

Y 2mf Y mo
F mA ðDÞ ¼ 2
ð52Þ
½Y mo D þ Y mf ð1  DÞ
2786 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

for the bilinear model, see (11)2,

F mB ðDÞ ¼ Y mo þ ðY mf  Y mo ÞD1=N ð53Þ

for the Allix & Ladevèze model, see (15)2, and

N
F mC ðDÞ ¼ Y mo þ ðY mf  Y mo Þ½ logð1  DÞ ð54Þ

for the exponential model, see (18)2, upon comparison of (50) and (51) the characteristic values of the mixed-
mode energy release rate Ymf are respectively obtained as

Y mfA ¼ GTc ð55Þ

1
Y mfB ¼ ½ðN þ 1ÞGTc  Y mo  ð56Þ
N
1
Y mfC ¼ Y mo þ ½GTc  Y mo  ð57Þ
CðN þ 1Þ

for the three models considered.

The complete traction–relative displacement relationships corresponding to (52),(52),(54) are shown in
Figs. 7–9, where the equivalent tractions (37) versus the two pure-mode relative-displacement components sunb
and susb are plotted.
The surfaces are constructed using radial proportional loading paths at varying angle u in the range [0, +p/2]
and the following material parameters: k þ n ¼ k s ¼ 10; 000; GcI ¼ 0:125; GcII ¼ 0:5; a1 ¼ a2 ¼ b1 ¼ b2 ¼ 2 with
GoI = 0.01125, GoII = 0.045 for the bilinear model, GoI = GoII = 0, N = 0.2 for the Allix & Ladèveze model
and GoI = GoII = 0, N = 1.7 for the exponential model.
Remark 2.4. It is worth noting that the proposed formulation allows a significant improvement in terms of
flexibility with respect to earlier mixed-mode models.
In particular, the use of a single scalar damage indicator instead of a damage variable of vectorial nature
avoids the introduction of further conditions on the evolution of damage, see e.g. [18,26] where, in order to
ensure that at complete debonding all the traction components vanish simultaneously, it assumed at the outset
that
Dn ¼ 1 ) Dt ¼ 1 ð58Þ
Dn and Dt denoting the damage variables in mode I and mode II, respectively.
As far as the decohesion propagation condition is concerned, in the cited works and in several others, see
e.g. [19,23,24,28], the mixed-mode failure locus that is analytically obtained for the case of isotropic damage is
given by
 a  a
GI GII
þ 1¼0 ð59Þ
GcI GcII

i.e. one has necessarily equal exponents for the two different single-mode contributions, a condition which
originates from the assumed form of the equivalent energy release rate. Accordingly, the failure locus can only
be either convex or concave and, as such, cannot account for mode interaction varying from negative to po-
sitive or vice-versa, as typically occurs for toughened epoxy composites, see also [48] and references therein.
The same limitation appears to be in the mixed-mode bilinear model discussed in [30], where the generalized
ellipse criterion can be fulfilled only with equal exponents for the two single-mode shares with, in addition, the
further condition:
GoI GoII
1 ¼1 ¼g ð60Þ
GcI GcII
that is introduced in order to characterize the damage variable.
 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2787

tδ
30.
0.0
2.0
4.0
6.0
8.0
15. 10.
[us ]
12.
0.04
14.
0.03
16.
18.
20.
22.

0.01 24.
26.
0.02
28.
0.03
30.
0.04 [un]
(a)

0.0
2.0
4.0
6.0
8.0
10.
12.
14.
16.
18.
20.
22.
24.
26.
28.
30.
(b)

Fig. 7. Traction–separation surface for the mixed-mode bilinear model. Side (a) and top (b) views.

2.3. Tangent operator

An essential ingredient for the solution of the discretized boundary value problem via Newton’s method is
the computation of the material tangent moduli for the interface model, i.e., the derivative of the traction vec-
tor with respect to the displacement jump vector:
Ktan ¼ d sut t ð61Þ
In order to work out a unique compact expression for the tangent operator Ktan, we introduce the two non-
linear complementary projection operators defined by
2788 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

tδ
30.
0.0
2.0
4.0
6.0
8.0

15. 10.
[us ] 12.
0.04
14.
0.03
16.
18.
20.
22.
24.
0.01
26.
0.02
28.
0.03 30.
0.04 [un]
(a)

0.0
2.0
4.0
6.0
8.0
10.
12.
14.
16.
18.
20.
22.
24.
26.
28.
30.
(b)

Fig. 8. Traction–separation surface for the mixed-mode Allix & Ladèveze model. Side (a) and top (b) views.

b
PðsutÞ ¼ hsun tiþ n þ sus ts
ð62Þ

PðsutÞ ¼ hsun ti n

so that the constitutive relationship for the traction vector (29)2 can be rephrased as
b PðsutÞ
t ¼ ð1  DÞ K b  PðsutÞ
þK  ð63Þ
where Kb and K are undamaged stiffness operators given by
 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2789

tδ
30.
0.0
2.0
4.0
6.0
8.0
15. 10.
[us ]
12.
0.04
14.
0.03
16.
18.
20.
22.

0.01 24.
26.
0.02
28.
0.03
30.
0.04 [un]
(a)

0.0
2.0
4.0
6.0
8.0
10.
12.
14.
16.
18.
20.
22.
24.
26.
28.
30.
(b)

Fig. 9. Traction–separation surface for the mixed-mode exponential model. Side (a) and top (b) views.

b ¼ k þ ðn  nÞ þ k s ðs  sÞ;
K  ¼ k  ðn  nÞ
K ð64Þ
n n

and the symbol  stands for the dyadic tensor product between vectors defined as
ða  bÞc ¼ ðb  cÞa 8a; b; c ð65Þ
whose component form is
ða  bÞij ¼ ai bj ð66Þ
The expression of the material tangent follows from the linearization of (63) with respect to the displacement
jump vector as
2790 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

b sut PðsutÞ
Ktan ¼ ð1  DÞ Kd b  sut PðsutÞ
þ Kd  b PðsutÞ
þ H D K b  d sut D ð67Þ
where the gradients of (62), that basically account for the sign of the normal displacement jump component,
are given by

b 1
d sut PðsutÞ ¼ ð1 þ sgnðsun tÞÞðn  nÞ þ ðs  sÞ
2 ð68Þ
 1
d sut PðsutÞ ¼ ð1  sgnðsun tÞÞðn  nÞ
2

while the term HD represents the Heaviside step function evaluated for the discrete damage multiplier c or,
equivalently, for the trial value of the damage mode (39), i.e.

1 if c > 0
H D ¼ H ðcÞ ¼ ð69Þ
0 otherwise

that is used to select between damage loading and elastic unloading.

In order to distinguish the terms accounting for the specific adopted form of the damage function Fm from
those depending upon the damage initiation criterion and decohesion propagation conditions, the damage
derivative dsubD is computed as

oD oD oD
d sut D ¼ d sut Y m þ d sut Y mo þ d sut GTc ð70Þ
oY m oY mo oGTc

In particular, the linearization of the mixed-mode energy release rate Ym yields:

b PðsutÞ
d sut Y m ¼ K b ð71Þ

while the gradients of the initial mixed-mode threshold Ymo and of the released energy at decohesion propa-
gation GTc are computed as
oY mo oGTc
d sut Y mo ¼ d sut b; d sut GTc ¼ d sut b ð72Þ
ob ob
Being:
ofo ofc
ðb; Y mo ðbÞÞ 6¼ 0; ðb; GTc ðbÞÞ 6¼ 0 ð73Þ
oY m oGT
see also (44) and (49), according to the implicit function theorem [49] the mixed-mode energies Ymo and GTc
are continuously differentiable with derivatives:
oY m ofo =ob 2Y mo ðg1 b2  g2 Þ
¼ ¼ ð74Þ
ob ofo =oY m Y m ¼Y mo ðbÞ bð1 þ b2 Þðg1 þ g2 Þ
oGTc ofc =ob 2GTc ðh1 b2  h2 Þ
¼ ¼ ð75Þ
ob ofc =oGT GT ¼GTc ðbÞ bð1 þ b2 Þðh1 þ h2 Þ
with
 a1  a2
Y mo Y mo b2
g1 ¼ a1 ; g2 ¼ a2 ð76Þ
ð1 þ b2 ÞGoI ð1 þ b2 ÞGoII
and
 b1  b2
GTc GTc b2
h1 ¼ b1 ; h2 ¼ b2 ð77Þ
ð1 þ b2 ÞGcI ð1 þ b2 ÞGcII
 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2791

The derivative of the mode-mixity parameter can be effectively computed starting from
rffiffiffiffiffiffi
Y II
b¼ ð78Þ
YI
and using (35) to get, after some algebra:
1 þ b2 1 þ b2 b
d sut b ¼ ½d sut Y m  ð1 þ b2 Þd sut Y I  ¼¼ b2 ðn  nÞ þ ðs  sÞ PðsutÞ ð79Þ
2bY m 2bY m

Collecting the previous results we thus obtain:

b PðsutÞ
d sut Y mo ¼ AðbÞ K b
ð80Þ
b PðsutÞ
d sut GTc ¼ BðbÞ K b

where A(b) and B(b) are rank-two tensors given by

AðbÞ ¼ a1 ðn  nÞ þ a2 ðs  sÞ; BðbÞ ¼ b1 ðn  nÞ þ b2 ðs  sÞ ð81Þ
with
 
Y mo g1 b2  g2 1
a1 ¼  ; a2 ¼  a1
Y m g1 þ g2 b2
  ð82Þ
GTc h1 b2  h2 1
b1 ¼  ; b2 ¼  2 b1
Y m h1 þ h2 b
In conclusion, the final form of the material tangent is arrived at as
b sut PðsutÞ
Ktan ¼ ð1  DÞ Kd b  sut PðsutÞ
þ Kd  þ H D Kb PðsutÞ
b
 
oD oD oD b PðsutÞ
b
 1þ AðbÞ þ BðbÞ K ð83Þ
oY m oY mo oGTc
1 being the rank-two identity tensor.
The above expression is well-posed for any mode-mixity ratio; actually, it is not difficult to show that one
has
½a1 ; a2 ; b1 ; b2  ! ½0; 1; 0; 1 ð84Þ
and
½a1 ; a2 ; b1 ; b2  ! ½1; 0; 1; 0 ð85Þ
+
in the limit for b ! 0 and for b ! +1, respectively.

3. Numerical examples

The interface model presented in the previous sections has been implemented as a part of the general-pur-
pose finite element code CAST3M [50], developed by the CEA (Commisariat à l’Energie Atomique, France).
For the numerical simulations described in the following, use has been made of 8-noded plane strain quadri-
laterals for the bulk material and 6-noded interface elements, for which Newton–Cotes integration is adopted
in order to reduce spurious oscillations of the traction profiles [51]. Computations have been carried out by
using a varying step size and a local-control-based arc-length algorithm as discussed in [33].
In order to compare the different cohesive laws considered, we have taken the same initial (undamaged)
interface stiffnesses for all of them, whose values have been identified via ultrasonic testing [37], and the same
interaction criteria for damage onset and decohesion propagation conditions for equal values of the charac-
teristic energies.
The critical energy release rates have been derived from linear beam theory using the classical Irwin–Kies
expression for the fracture energy [52]:
2792 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

P 2 oC
Gc ¼ ð86Þ
2B oa
where P is the load, B the specimen width, a the crack length and C the compliance. The characteristic energies
GoI, GoII have been respectively taken as a fraction of GcI, GcII in absence of any further information; however,
for the numerical examples considered hereafter, these values are found not to influence appreciably the com-
puted global response.
A pictorial representation of the pure-mode traction–relative displacement relationships used in computa-
tions is given in Fig. 10 and the full detail of the material parameters is reported in Table 1.

3.1. Example 1. ENF test

As a first example we consider the end-notched flexure (ENF) test, commonly used for obtaining pure
mode-II loading conditions. The samples tested have been prepared as shown in Fig. 11. They consist of
two 3 mm thick and 20 mm wide aluminum arms bonded with a layer of 0.5 mm resin adhesive made of diglyc-
idyl ether of bisphenol A (DGEBA) with diamino dyphenol methane (DDM) as curing agent. The test has
been performed via displacement control with a three-point bending fixture with span 2L = 120 mm and an
initial crack of length a0 = 17 mm, that was obtained by introducing a Teflon film.
The experimental load–deflection response is plotted in Fig. 14 along with the curves computed via FE
analysis, that refer to a regular mesh made of 1440 continuum elements and 240 interface elements, see also
Fig. 12.
The numerical predictions reproduce the experimental results with sufficient accuracy; moreover, no sub-
stantial difference appears between the responses obtained with the three different interface models, that yield
practically coincident results owing to the small stiffness of the sample and the very refined FE mesh used for
the simulation.
The very refined character of the adopted mesh can be recognized from Fig. 13, where the size of the evolv-
ing cohesive process zone for the three different damage models is shown. Actually, the present length of the
cohesive elements (0.5 mm) allows several of them to be included within the cohesive process zone; this results
in quite smooth responses computed via finite element analyses, that do not exhibit false instabilities caused by
coarse meshes, as reported by several analysts [30,33].

8
Bilinear model
Allix-Ladevèze model
7 Exponential model

5
traction t

0
0 0.005 0.01 0.015 0.02 0.025 0.03
displacement jump [u]

Fig. 10. Traction–separation relationships used in FE computations.

 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2793

Table 1
Material data set
Bulk material
E (GPa) 75
m 0.3
Interface
Bilinear model Allix & Ladevèze model Exponential model
kþ 3
n (N/mm ) 760 760 760
ks (N/mm3) 810 810 810
GoI (N/mm) .008 .008 .008
GoII (N/mm) .036 .036 .036
GcI (N/mm) .02 .02 .02
GcII (N/mm) .09 .09 .09
a1 2 2 2
a2 2 2 2
b1 2 2 2
b2 2 2 2
N – 1 1

3 mm

17 mm

20 mm 60 mm

Fig. 11. ENF test. Model problem.

Fig. 12. ENF test. Finite element mesh.

The unstable character of the response at decohesion propagation and the influence of friction at crack
sides due to through-the-thickness compression is evidenced in the experimental response reported in
Fig. 14, see also [54] and references therein for a detailed discussion on the effects of friction in mode-II dela-
mination tests. The fact that friction is not accounted for in the model explains the intersection of the last part
2794 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

16

14

Cohesive process zone length (mm)

12

10

2 Bilinear
Allix-Ladevèze
Exponential
0
10 20 30 40 50 60 70 80 90 100
Crack length (mm)

Fig. 13. ENF test. Cohesive zone versus crack length.

600

550

500

450

400

350
Load (N)

300

250

200

150

100 Bilinear
Allix-Ladevèze
50 Exponential
Experiment
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8
Deflection (mm)

Fig. 14. ENF test. Experimental and FE computed load–deflection curves.

of the experimental curve with those obtained via FE analysis, see also Fig. 15 for a detail of the crossing of
the experimental and numerically computed curves.
In Fig. 16 is plotted the experimentally measured load–deflection response in conjunction with the analyt-
ical solution obtained via linear beam theory. Since the specimen is relatively thin, transverse shear deforma-
tion is neglected.
The deflection of the loaded point at decohesion propagation, obtained assuming that the advancement of
the initial crack takes place at constant GII = GcII, is computed as
 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2795

600

550

500
Load (N)

450

400

350 Bilinear
Allix-Ladevèze
Exponential
Experiment
300
1.4 1.6 1.8 2 2.2 2.4 2.6 2.8
Deflection (mm)

Fig. 15. ENF test. Experimental and FE computed load–deflection curves. Enlarged view.

600

550

500

450

400

350
Load (N)

300

250

200

150
Initial loading line
100 propagation (a<L)
propagation (a>L)
50 splitted beams
Experiment
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8
Deflection (mm)

Fig. 16. ENF test. Experimental and theoretical load–deflection curves.

pffiffiffi
3P 3 L3 þ 256 3ðGcII BEJ Þ3=2
D¼ for a < L
144P 2 EJ ð87Þ
p ffiffiffi 3=2
3P 3 L3  64 3ðGcII BEJ Þ
D¼ for a > L
36P 2 EJ
where J denotes the second area moment of one arm of the structure. The initial and final loading lines are
obtained by setting respectively a = a0 and a = 2L in the following relationships:
2796 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

P ð2L3 þ 3a3 Þ
D¼
96EJ
3
ð88Þ
P ½8L3  3ð2L  aÞ 
D¼
96EJ
We emphasize that the above analytical solution is intended here only as an approximate solution for com-
parison purposes only; many possibilities exist for improving the agreement with FE analyses results by
including, for example, large deflections [53], crack length corrections [54], and allowing for the difference be-
tween plane stress and plane strain.

3.2. MMF test

The second example concerns the mixed-mode flexure (MMF) test originally proposed in [55]. The testing
apparatus and the dimensions of the specimen are the same as for the ENF test, see Fig. 17. Loading is sim-
ulated via displacement control and the initial flaw is a0 = 23 mm.
The MMF test is mode-I p dominated
ffiffiffi with a mode mixity ratio that, for equal thicknesses of the lower and
upper arms, is evaluated as 3=2 using linear beam theory.
The load–deflection curves computed via FE analysis are depicted in Fig. 20; these have been obtained
using a regular mesh consisting of 1323 continuum elements for the bulk material and 214 interface elements,
see also Fig. 18. As for the previous example, in Fig. 19 is reported the size of the evolving cohesive process
zone during decohesion propagation, which shows the adequacy of the adopted finite element mesh. The
responses computed for the different shapes of the traction–displacement jump relationships are in substantial
agreement and the comparison with the experimental response, that is also plotted in Fig. 20, shows a fairly
good correspondence up to the peak load, that is well captured.
On the contrary, a more pronounced difference is found during crack propagation where the experimentally
measured load drops down more rapidly with respect to the numerical simulations. This fact can be appreci-
ated also by comparing the experimental curve with the analytical ones evaluated via linear beam theory, see
Fig. 21.
The analytical non-linear propagation curves are obtained by first computing the midspan deflection as

P ð2L3 þ 7a3 Þ
D¼ for a < L
96EJ ð89Þ
P ð40L3 þ 84L2 a  42La2 þ 7a3 Þ
D¼ for a > L
96EJ
pffiffiffi
and then eliminating the crack length a using the fixed mode ratio b ¼ 3=2 and

a¼x for a < L

ð90Þ
a ¼ 2L  x for a > L

3 mm

23 mm

20 mm 60 mm

Fig. 17. MMF test. Model problem.

 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2797

Fig. 18. MMF test. Finite element mesh.

3.5

3
Cohesive process zone length (mm)

2.5

1.5

0.5 Bilinear
Allix-Ladevèze
Exponential
0
20 30 40 50 60 70 80 90
Crack length (mm)

Fig. 19. MMF test. Cohesive zone versus crack length.

for
1=2
8½ð1 þ b2 ÞBEJGcI GcII 
x ¼ x1 ¼ pffiffiffi ð91Þ
7P ðb2 GcI þ GcII Þ1=2
or
8½ð1 þ b2 ÞBEJGcI GcII 1=2
x ¼ x2 ¼ pffiffiffi 1=4
ð92Þ
7P ðb4 G2cI þ G2cII Þ
that respectively correspond to a linear (b2 = b1 = 1) and quadratic (b2 = b1 = 2) interaction according to the
decohesion propagation criterion (47).
2798 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

180

160

140

Load (N) 120

100

80

60

40
Bilinear
20
Allix-Ladevèze
Exponential
Experiment
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2
Deflection (mm)

Fig. 20. MMF test. Experimental and FE computed load–deflection curves.

180

160

140

120
Load (N)

100

80

60
Initial loading line
propagation (a<L) - lin
40 propagation (a>L) - lin
propagation (a<L) - qua
20
propagation (a>L) - qua
splitted beams
Experiment
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2
Deflection (mm)

Fig. 21. MMF test. Experimental and theoretical load–deflection curves.

It is noted that the experimental response in the decohesion propagation phase is very close to the loading
line of the completely splitted specimen, thus suggesting a sudden rupture of the ligament right after the max-
imum load has been attained.
This is actually what was recorded during the test, that has been conducted with a crosshead displacement
rate of approximately 0.005 mm/s. A similar behaviour for MMF bonded specimens, i.e. unstable crack prop-
agation at a displacement rate of 0.5 mm/min, has been also reported in [56].
 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2799

4. Summary and conclusions

A cohesive-zone model has been developed for the analysis of interface decohesion in adhesively bonded
assemblies. The proposed approach relies upon a quite general damage-mechanics formulation that provides
deeper insights as well as several improvements with respect to earlier interface models, two of which, the
bilinear one implemented by Crisfield and co-workers and the Allix & Ladevèze polynomial model, are shown
to be recovered as special cases.
The characteristic features of the present model are mainly put forward in the mixed-mode case, whose
treatment is based upon the use of a single scalar damage indicator and a consistently defined work-conjugate
energy release rate variable. This ensures that, on one hand, all the traction components vanish simultaneously
if complete debonding occurs and, on the other side, that the form of the failure locus can be taken indepen-
dently from the definition of the equivalent damage energy release rate.
The flexibility of the model is further evidenced in the derivation of the linearized form of the constitutive
equations, where the contributions to the material tangent originating from the specific form of the traction–
relative displacement relationship and those accounting for the interaction criteria for damage onset and dec-
ohesion propagation are clearly distinguished.
The model has been implemented in a finite element code using interface elements and a local-control-based
arc-length algorithm. In the numerical examples computations have been carried out for different cohesive
laws endowed with equal critical energy release rates; for the relatively flexible samples considered a substan-
tial agreement of the results have been recorded, thus confirming that fracture behaviour is in these cases truly
energy-driven.
Comparisons with available experimental data have also been considered, that have shown quite good pre-
dictive capabilities of the model.

Acknowledgements

This work has been partly carried out within the framework of the activities of the Laboratoire Lagrange,
an European research group gathering CNR, CNRS, Università di Roma ‘‘Tor Vergata’’, Université de
Montpellier II, ENPC, LCPC.
NV also acknowledges the Italian National Research Council (CNR) and the Italian Ministry of Educa-
tion, University and Research (MIUR) for the financial support through the Cofin 2003 research program
Materials with microstructure for technological innovation in civil structures coordinated by prof. Franco
Maceri.

References

[1] Barenblatt GI. The formation of equilibrium cracks during brittle fracture: general ideas and hypotheses, axially symmetric cracks.
Appl Math Mech 1959;23:622–36.
[2] Barenblatt GI. The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 1962;7:55–129.
[3] Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8:100–4.
[4] Hillerborg A, Modéer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics
and finite elements. Cement Concrete Res 1976;6:773–82.
[5] Hillerborg A. Application of fictitious crack model to different types of materials. Int J Fract 1991;51:95–102.
[6] Allix O. Délaminage. Approche par la mécanique de l’endommagement. In: Fouet JM, Ladevèze P, Ohayon R, editors. Calcul des
Structures et Intelligence Artificielle, Paris, 1987. p. 39–53, Pluralis.
[7] Needleman A. A continuum model for void nucleation by inclusion debonding. J Appl Mech 1987;54:525–31.
[8] Needleman A. An analysis of tensile decohesion along an interface. J Mech Phys Solids 1990;38(3):289–324.
[9] Tvergaard V. Effect of fibre debonding in a whisker-reinforced metal. Mater Sci Engng 1990;A125:203–13.
[10] Allix O, Ladevèze P. Interlaminar interface modelling for the prediction of delamination. Compos Struct 1992;22(4):235–42.
[11] Tvergaard V, Hutchinson JW. The relation between crack growth resistance and fracture process parameters in elastic–plastic solids.
J Mech Phys Solids 1992;40(6):1377–97.
[12] Corigliano A. Formulation, identification and use of interface models in the numerical analysis of composite delamination. Int J
Solids Struct 1993;30(20):2779–811.
[13] Daudeville L, Ladevèze P. A damage mechanics tool for laminate delamination. Comput Struct 1993;25:547–55.
2800 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801

[14] Schellekens JCJ, de Borst R. Free edge delamination in carbon epoxy laminates: a novel numerical/experimental approach. Compos
Struct 1993;28(14):357–73.
[15] Tvergaard V, Hutchinson JW. The influence of plasticity on mixed mode interface toughness. J Mech Phys Solids 1993;41(6):1119–35.
[16] Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42:1397–434.
[17] Allix O, Ladevèze P, Corigliano A. Damage analysis of interlaminar fracture specimens. Compos Struct 1995;31(1):61–74.
[18] Allix O, Corigliano A. Modeling and simulation of crack propagation in mixed-modes interlaminar fracture specimens. Int J Fract
1996;77:111–40.
[19] Allix O, Ladevèze P. Damage mechanics of interfacial media: basic aspects, identification and application to delamination. In:
Voyiadjis GZ, Allen DH, editors. Damage and interfacial debonding in composites, number 44 in studies in applied mechan-
ics. Amsterdam: Elsevier; 1996. p. 167–88.
[20] Camacho GT, Ortiz M. Computational modelling of impact damage in brittle materials. Int J Solids Struct 1996;33(20–22):2899–938.
[21] Point N, Sacco E. A delamination model for laminated composites. Int J Solids Struct 1996;33(4):483–509.
[22] Chaboche JL, Girard R, Schaff A. Numerical analysis of composite systems by using interphase/interface models. Comput Mech
1997;20:3–11.
[23] Allix O, Lévêque D, Perret L. Identification and forecast of delamination in composite laminates by an interlaminar interface model.
Compos Sci Technol 1998;58(5):671–8.
[24] Ladevèze P, Allix O, Gornet L, Lévêque D, Perret L. A computational damage mechanics approach for laminates: identification and
comparison with experimental results. In: Voyiadjis GZ, Ju JW, Chaboche JL, editors. Damage mechanics in engineering materials,
number 46 in studies in applied mechanics. Amsterdam: Elsevier; 1998. p. 481–500.
[25] Mi Y, Crisfield MA, Davies GAO, Hellweg HB. Progressive delamination using interface elements. J Compos Mater
1998;32(14):1246–72.
[26] Allix O, Corigliano A. Geometrical and interfacial non-linearities in the analysis of delamination in composites. Int J Solids Struct
1999;36(15):2189–216.
[27] Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J
Numer Meth Engng 1999;44(9):1267–82.
[28] Corigliano A, Allix O. Some aspects of interlaminar degradation in composites. Comput Meth Appl Mech Engng 2000;185(2–4):
203–24.
[29] Ladevèze P, Allix O, Deü JF, Lévêque D. A mesomodel for localisation damage computation in laminates. Comput Meth Appl Mech
Engng 2000;183:105–22.
[30] Alfano G, Crisfield MA. Finite element interface models for the delamination analysis of laminated composites: mechanical and
computational issues. Int J Numer Meth Engng 2001;50(7):1701–36.
[31] Corigliano A, Mariani S, Pandolfi A. Numerical modeling of rate-dependent debonding processes in composites. Compos Struct
2003;61(1–2):39–50.
[32] de Borst R. Numerical aspects of cohesive-zone models. Engng Fract Mech 2003;70(14):1743–57.
[33] Alfano G, Crisfield MA. Solution strategies for the delamination analysis based on a combination of local-control arc-length and line
searches. Int J Numer Meth Engng 2003;58:999–1048.
[34] Lemaitre J, Chaboche JL. Mechanics of solids materials. Cambridge: Cambridge University Press; 1990.
[35] Lavrentyev AI, Rokhlin SI. Ultrasonic spectroscopy of imperfect contact interfaces between a layer and two solids. J Acoust Soc Am
1998;103(2):657–64.
[36] Lowe MJS, Challis RE, Chan CW. The transmission of Lamb waves across adhesively bonded lap joints. J Acoust Soc Am
2000;107(3):1333–45.
[37] Vlasie V, Rousseau M. Acoustical validation of the rheological models for a structural bond. Wave Motion 2003;37(4):333–49.
[38] Bertsekas DP. Constrained optimization and Lagrange multiplier methods. New York: Academic Press; 1982.
[39] Andrews GE, Askey R, Roy R. Special functions. Cambridge: Cambridge University Press; 1999.
[40] Lanczos C. A precision approximation of the Gamma function. SIAM J Numer Anal 1964;1:86–96.
[41] Chaboche JL, Feyel F, Monerie Y. Interface debonding models: a viscous regularization with a limited rate dependency. Int J Solids
Struct 2001;38(18):3127–60.
[42] Crisfield MA, Alfano G. Adaptive hierarchical enrichment for delamination fracture using a decohesive zone model. Int J Numer
Meth Engng 2002;54(9):1369–90.
[43] Duvaut G, Lions JL. Les Inéquations en Mécanique et en Physique. Paris: Dunod; 1972.
[44] Perzyna P. Fundamental problems in viscoplasticity. Adv Appl Mech, 9. New York: Academic Press; 1966. p. 243–377.
[45] Corigliano A, Ricci M. Rate-dependent interface models: formulation and numerical applications. Int J Solids Struct 2001;38(4):
547–76.
[46] Suo Z, Hutchinson JW. Interface crack between two elastic layers. Int J Fract 1990;43:1–18.
[47] Hashin Z. Failure criteria for unidirectional composites. J Appl Mech 1980;47:329–34.
[48] Reeder JR. An evaluation of mixed-mode delamination failure criteria, Technical Memorandum 104210, NASA, 1992.
[49] Narasimhan R. Analysis on real and complex manifolds. 2nd edition. Amsterdam: North-Holland; 1985.
[50] CEA, Available from: [Link] CAST3M – User Manual, 2003.
[51] Schellekens JCJ, de Borst R. On the numerical integration of interface elements. Int J Numer Meth Engng 1993;36:43–66.
[52] Irwin GR, Kies JA. Critical energy release rate analysis of fracture strength. Weld J Res Suppl 1954;33:193–8.
[53] Wang Y, Williams JG. Corrections for mode-II fracture toughness specimens of composite materials. Compos Sci Technol
1992;43(3):251–6.
 N. Valoroso, L. Champaney / Engineering Fracture Mechanics 73 (2006) 2774–2801 2801

[54] Blackman BRK, Kinloch AJ, Paraschi M. The determination of the mode II adhesive fracture resistance, GIIc, of structural adhesive
joints: an effective crack length approach. Engng Fract Mech 2005;72(6):877–97.
[55] Russell AJ, Street KN. Moisture and temperature effects on the mixed-mode delamination fracture of unidirectional graphite/epoxy.
In Johnson WS, editor. Delamination and debonding of materials, ASTM STP 876, Philadelphia. American Society for Testing
Materials, 1985. p. 349–70.
[56] Loh WK, Crocombe AD, Abdel Wahab MM, Ashcroft IA. Environmental degradation of the interfacial fracture energy in an
adhesively bonded joint. Engng Fract Mech 2002;69(18):2113–28.