329-359

Computer methods
in applied
mechanics and
englneerlng
ELSEVIER Comput. Methods Appl. Mech. Engrg. 140 (1997) 329-359
A discrete formulation for elastic solids with damaging

interfaces’
Gabriella Bolzon, Albert0 Corigliano
Department of dbuctural Engineering, Politecnico of Milano, Piazza Leonardo da Vinci 32, Milano, I-20133, Italy
Received 19 September 1995; revised 2 May 1996
Abstract
An elastic solid with embedded interfaces, loci of possible displacement discontinuities, is considered here. Decohesion and
quasi-brittle fracture processes are simulated by making use of softening interface laws which relate tractions to displacements
jumps. A discrete formulation in terms of interface variables only, is obtained. The space discretization is carried out by means of
a mixed finite element approach in which all interface variables are modelled. Study of uniqueness of the rate problem formulated
in terms of interface variables and of stability of the equilibrium states is presented. Some examples are shown in order to clarify
the formulation.
1. Introduction
Many structural problems require the consideration of the appearance and the evolution of important
damage phenomena. These can cause local ruptures and material decohesion, and then develop in
macroscopic fractures and separation between parts. Due to the lack of homogeneity and/or anisotropy
of the structure, or due to particular loading and constraint conditions, concentrated damage and
macroscopic separation of parts is in many cases forced to appear in a limited number of fixed critical
zones, while the remaining part of the structure can still be considered in the elastic regime. Then, a
suitable schematization for this kind of global behaviour is that of an elastic medium with a discrete set
r of embedded interfaces c on which displacement discontinuities are possible. In this context,
interface r is assumed to be a set of zero measure, i.e. r consists of a set of surfaces in three
dimensional situations, lines in two dimensional cases and points for one dimensional problems.
In order to simulate the progressive damage which concentrates along r, a softening interface
constitutive law, described by a relationship between tractions and displacement discontinuities, is
attributed to the points of r. The cohesion between surfaces initially in contact along r gradually
vanishes, until a complete separation is reached.
The present paper focuses on the class of structural problems amenable to the above schematization,
and aims to obtain a unified formulation for finite element computations.
Representative examples of typical situations which can be conceived are those briefly described in
the following paragraphs.
Fracture of quasi-brittle materials when the fracture path is known a priori e.g. for symmetry
conditions and the presence of initial defects. As specific references, the cohesive-crack models for
‘Dedicated to Prof. Giulio Maier on the occasion of his 65”’ anniversary.
0045-7825/97/$17.00 10 1997 Elsevier Science S.A. All rights reserved

PII SOO45-7825(96)01101-2
330 G. Bolzon, A. Corigliano I Compur. Methods Appl. Mech. Engrg. 140 (1997) 329-359
concrete and concrete-like materials (see e.g. [l-6]), comprising the case of joints separating rocks in
geo-mechanics [7,8] can be quoted.
Related meaningful applications refer also to masonry structures (see e.g. [9-111). In this case, when
the damage inside the bricks is considered to be of minor importance, the set of interfaces is naturally
represented by the mortar joints where the structural degradation is thought to be concentrated.
Composite materials, particularly fibre reinforced composites, show different damage mechanisms.
Among these, particle-matrix and fibre-matrix debonding, layer-layer debonding or delamination (see
e.g. [12-191) are suited to the proposed schematization.
Other typical applications for which the scheme of an elastic continuum with damaging interfaces
naturally applies are those connected with the study of glued elements. In this case the interfaces consist
of the glue layers [20,21].
The same schematization can be applied in the study of post-localisation processes (see e.g. [22-24]),
with the interface represented by the thin localisation band in which displacement discontinuities can be
thought of as developing. This class of problems has been studied since the early seventies, with
reference to elastic-plastic frames with softening hinges (see e.g. [25-281). In this latter case, interfaces
coincide with the softening hinges which form in those critical sections where plastic strains strongly
localise. As natural extensions, elastic-plastic plates with linear softening hinges can be also conceived.
An unified formulation is presented here for the study of the problems listed above, and is used as a
basis for the derivation of a discrete model. This is done within the so-called generalized variables
approach specialised to interface variables.
The notion of generalised variables, introduced by Prager [29] with reference to beam problems, was
developed by Maier [30] and by Corradi and co-workers [31-331. Recently, Comi et al. [34] and Comi
and Perego [35] proposed generalised variable formulations for elastic-plastic problems, Comi and
Perego [36] and Comi and Corigliano [37] for gradient-dependent plasticity in static and dynamic cases
respectively.
The discrete formulation introduced here is reduced to interface variables only, thus obtaining a
discrete set of equations which are particularly advantageous both from a computational point of view
and for the study of uniqueness and stability of the discretized rate problem.
An outline of the paper is as follows.
In Section 2, the governing relations for the problem of an elastic body with embedded interfaces are
presented. The softening interface law is formulated in the class of the so-called standard dissipative
systems, and some examples are given.
Section 3 deals with a weak formulation for the problem under study. All the relations involving the
interface, included compatibility, are presented in a weak format, while elastic behaviour and
compatibility inside the elastic portions of the body are assumed to be a priori satisfied.
The discrete finite-element formulation is derived in Section 4. The governing equations of the
discrete (in space) problem are reduced to interface variables only.
In Section 5, uniqueness of the solution of the rate problem and stability of the equilibrium states are
discussed. This is done by following both the classical approach of Hill [38] and the approach (here
called intrinsic) which makes reference to the work by Maier [25], Maier and Drucker [39], and Nguyen
and co-workers [40-441. Sufficient conditions for uniqueness and stability are explicitly given for the
problem when it is reduced to interfaces.
Finally, Section 6 presents some examples aimed at clarifying the formulation.
2. Problem formulation
As a suitable schematization of the situations described in Section 1, an elastic body 0 is considered

here with the following assumptions (see Fig. 1):
l strain and displacements are small;

0 loads tand displacements U are assigned on the parts aL$ and 80, of boundary aa, body forces b are
assigned in 0;
G. Bolzon, A. Corigliano I Comput. Methods Appl. Mech. Engrg. I40 (I997) 329-359 331
Fig. 1. Elastic body with damaging interfaces.

Fig. 2. Notations for the interface I:
l a set of n interfaces r = U y=, c is given. Interfaces 4 are surfaces embedded in three dimensional
(3D) domains, lines in 2D domains and points in 1D domains.
l displacement diacontinuities w are possible on interfaces r.
Each interface 1: is thought to divide the body R into the portions 0 + and fi -, and the unit vector
normal to & is oriented toward 0 + , see Fig. 2; R’ is defined as 0’ =0-r; tractions and
displacements re:levant to surfaces r+ and r- are respectively denoted by t+ , u+ and t- , u-. In the
following text the distinction among r+, r- and r will be considered, r+ and r- being the sides of
r on 0 + and R - respectively, as defined in Fig. 2.
The problem to be solved is that of finding the equilibrium configuration(s) of the body 0 under the
assigned loads.
2.1. Governing relations
The problem introduced in the above section is governed by the following relationships:
0 compatibility
E=CU inn’ (1)

w=u+- U-=C+I 0nr (2)
u=U onfU&, (3)

l equilibrium
CTa+b=O in0 (4)
t=t-= -t+=e 0nr (5)
No = t on an,, (6)
l constitutive equation for the elastic body 0’
844 in n,
a=ds (7)
l constitutive equation for the interfaces r

w = we + wp (8)
t= wW, 4’. a= _ w(wey4

aa (9)
ad ’
vT(44 i. VT@,
4 li
kp _
_--
at 9
&=
aa (10)
332 G. Bolzon, A. Corigliano I Comput. Methods Appl. Mech. Engrg. 140 (1997) 329-359
A 20; f(4 4 s 0; f’i = 0. (11)
In the above equations:
l C represents the differential operator of linear compatibility in matrix form, therefore CT is the
matrix of equilibrium;
l vectors u and E respectively collect the six independent components of the Cauchy stress tensor and
of the small strain tensor;
l N is the matrix whose entries are the components of the outward normal to the interface r or to the
boundary &$.
An elastic behaviour is assumed for the body 0 through Eq. (7), where o is the strain energy density
function. In the following, a linear elastic behaviour is considered in 0, with w being a quadratic
function of E.
The set of relations (8)-(11) represents the interface constitutive behaviour, linking the displacement
discontinuity vector w to the traction vector t on r. In Eq. (8) the displacement jump w is given by the
sum of a reversible part we and an unrecoverable part wp. Variables wp and a set of further kinematic
internal variables (Y govern the irreversible behaviour of the interface. State Eq. (9) give tractions t and
static internal variables a as derivatives of a free energy function $(we, cu), defined per unit surface lY
Differentiability but not convexity is required for the free energy function +,, thus allowing us to
reproduce softening. The evolution of the kinematic internal variables wp and (Y is governed by
relations (10) and (ll), where f(t, a) is a set of activation functions and A a vector of plastic multipliers.
Relations (10) and (11) define an associative flow rule. Functions f(t, a) are required to be convex and
differentiable.
With the above assumptions, the interface constitutive law belongs to the category of standard
dissipative models [40,42]. The Hill’s maximum dissipation principle is then fulfilled, i.e.:
GtTip+SaTiy~O;V&Jia EE,, E,,={&=t* -t,&z=u* -ulf(t*,u*)sO}. (12)
The restriction to this category of constitutive laws allows us to derive results for stability and
bifurcation of the considered interface degradation processes as shown in Section 5.
2.2. Interface behuviour
Several widely employed categories of interface behaviours can be obtained by specialisation of Eqs.
(8-11) and many mechanically meaningful situations, ranging from perfectly continuous bodies to
fracture problems, can thus be modelled. In this section, some representative situations are discussed.
Perfectly bonded surfaces can be obtained by simply assuming w =0 on r (or on a portion of it).
Elastic (reversible) behuviour can be derived under the assumptions: Cc, = Jl(w’), wp =0 and by using
Eq. (9a) only, for the definition of the interface law. In this case, the free energy $ acquires the
meaning of strain energy per unit surface. As a particular case, $ =constant gives zero traction vector t
on r, thus representing the situation of disconnected bodies. The non-linear elastic behaviours which
can be obtained for the interface, by varying the choice of function $I, can be used in order to simulate
elastic adhesive joints (see e.g. [20,21]), fibre debonding in composites (see e.g. [12]), or progressive
crack propagation in brittle and quasi-brittle homogeneous solids (see e.g. [45]). In these cases, the
degradation of the interface is simulated through the introduction of an elastic, possibly non-linear
softening law but, obviously, local unloading or cyclic loading situations cannot be correctly taken into
account.
In simulating crack propagation, the hypotheses introduced above can be contemporary considered
by distinguishing in r those portions which are completely damaged (t = 0), damaged but still able to
transmit some traction and undamaged (w = 0).
In addition to the cases mentioned above, the use of an elastic interface model allows us to reproduce
more classical situations in which elastic bodies are connected one to each other by means of elastic
G. Bolzon, A. Corigliano I Comput. Methods Appl. Mech. Engrg. 140 (1997) 329-359 333
springs. Beams over elastic foundations and elastic joints in frames are practical cases which belong to
this category.
Another meaningful case which can be represented by means of a potential function 9 is that of
unilateral contact. As an example, the one-dimensional behaviour shown in Fig. 3a is obtained from the
continuous and differentiable function Jl(w) = K( w ) ‘_/2 represented in Fig. 3b, where the symbol ( 0) _
means the negative part of 0, i.e. (0) _ = l for 0~0; (0) _ =0 for l>O.
The cohesive model recently used e.g. in [3-61 to characterize the fracture behaviour of concrete-like
materials in mode I, can be represented by the piece-wise linear law drawn in Fig. 4. Below the critical
level t, of interface tractions, a perfect bonded interface is assumed with zero displacement dis-
continuities; after attaining the value t,, a (linear) softening branch is activated. In the softening phase,
the behaviour can be considered as rigid-plastic with fully unrecoverable displacement discontinuity W.
The softening regime ends when a critical value w0 is reached with zero traction. The subsequent phase
representing completely debonded interfaces, can be described by an elastic behaviour with constant $
as in the case of disconnected bodies.
Let TVrepresent the current time instant during an hypothetical load history, and define the following
sets r,, r, and r, of points x belonging to r:
(I3a)
rP = {x E Tl~!(x, 7) = t, for some 7 C T,}, (13b)
r,={xErllv>wo}, (13c)
An analytical description of the behaviour in Fig. 4 is given by the following assumptions which
specialise th$ interface model (8)-(11):
w=O on<!, (144
$=tj/(a)=f,,a+;hc?; f=t+a onr,, (14b,c)
*=O enc., (14d)
(4
Fig. 3. Unilateral contact: a) interface law; b) corresponding energy function.

Fig. 4. The interface law for the mode-1 cohesive model.
Eqs. (14a), (14d) represent respectively the conditions of perfectly bonded interfaces and completely
disconnected bodies already introduced in the previous discussion. Eqs. (14b,c) define the free energy
function and the activation function to be introduced in Eqs. (8-11) in order to recover the complete
single-mode interface model. In Eq. (14b) h = -t,lw,, represents the negative slope of the traction-
displacement discontinuity diagram of Fig. 4. It is worth noticing that, since a rigid-plastic-softening
model is considered for r,, we=0 and, as a consequence, Eq. (9a) is useless and the set of kinematic
internal variables cx reduces in this case to w = rvp=A.
A unified description for the behaviour on q and < can be obtained by replacing the definition (14b)
of function $, with the following one:
@=+(a)=t,n+;ha2
-;h(a-iv,>:, (15)
where the symbol (o)+ means the positive part of 0, i.e. ( l)+ =* for 030; (o)+ =0 for 0~0.
By replacing the quadratic term ha’/2 in Eqs. (14b) and (15) with a generic non-linear, non-convex
function of W, it is possible to obtain softening curves of different shapes which can better reproduce
some experimental results (see e.g. [2,46,47]).
A simple generalization to mixed-mode conditions of the above pure-mode cohesive-crack model
given by Eqs. (13)-(14) can be obtained by defining as kinematic internal variables cy= [cr, at] and by
introducing functions $ and f governing the behaviour on the cohesive process zone rp as follows:
(16)
(17)
where parameters t,_, t,“, h,, h, are defined in Fig. 5, while a,, and a, are the static internal variables
conjugated to cu, and a,.‘ihrough Eq. (9b). Portion c of the interface r where t, = t, = t, =O, is given by:
r, = {x E r/a,(x) > 9,; c+) > LytO}. (18)
The elastic domain defined by relations (llb) and (17) in the space ti is the cone represented in Fig.
6. Subscripts 1, 2, 3 denote respectively the direction 1 normal to the interface and directions 2, 3
tangent (parallel to the interface) as suggested by Fig. 2. For pure mode II situations the evolution of
the elastic domain is such that the cone progressively reduces to the interval t, < -a, on axis t, , keeping
(4
Fig. 5. The interface law for the mixed-mode cohesive model: a) opening mode; b) shearing mode.
Fig. 6. The elastic domain for the mixed-mode cohesive model.

fixed its vertex. In pure mode I, the cone vertex moves toward the origin of the reference frame
keeping constant the intercepts with the t,- t, plane.
The interface model described by Eqs. (8-11) supplemented by definitions (16)-(17) can be seen as a
prototype of more complete models, like those defined by Stankowski et al. [19].
An elastic-plastic interface constitutive law similar to the rigid-plastic cohesive crack model above,
but including a reversible (elastic) behaviour, can be formulated by defining (see [16-181)
(19)
(20)
In Eq. (19) K, , K2, K, are interface elastic stiffnesses with the dimension of a force over a length
cube; &(a) is a non-convex function of the unique kinematic internal variable CYgoverning the
degradation of the interface. A possible choice for t&(a), which gives a linear softening branch as
previously considered, is G2(~) = (Y+ 4 ha *, with h a negative non-dimensional parameter. The
introduction of the elastic part of the interface is physically justified when the interface is seen as the
limit of a thin layer for vanishing thickness.
The activation function f of Eq. (20) defines an elastic domain with the shape drawn in Fig. 7; the
role of parameters t,, and t,, is represented in the picture, while t,, and t,, are related to t,, as:
1 ltl, + 1/to, = 1lt,, . The evolution of the elastic domain is governed by the kinematic internal variable (Y
as in isotropic softening plasticity; the final elastic domain consists of the interval t, ~0 on axis t,. The
portion r, of the interface where t, = t, = t, =0 is now defined as follows:
r,= xEr]cX(x)> -+ . (21)

1 >
Interface dumqge models can also be recast in the class of interface models described by relations
(8)-(11). A particularly simple interface damage model is the one described by the functions:
(22)
(23)
In Eq. (22) kinematic internal variables (Y~,i= 1, 2, 3, with aiEIO, l] are damage variables which
Fig. 7. The elastic domain for the elastic-plastic interface model.

Fig. 8. Interface damage models: a) brittle; b) progressive.
affect the behaviour in opening and sliding modes; (Yeis a further kinematic internal variable. The
choice of function I,$((Y~)has a direct influence on the shape of the softening branch in the traction-
displacement jump diagram. As meaningful examples: the choice &(a,)=0 gives the brittle behaviour
of Fig. 8a; the choice (J12(@$)=q++& with h a positive non-dimensional parameter, gives the curve
of Fig. 8b (see [42]).
3. Weak formulation
A weak formulation of the problem described in Section 2.1 is here given. This allows for the
introduction of a spatial discretization through finite elements as discussed in Section 4. In the present
formulation, compatibility in 0’ and the linear elastic constitutive law in 0’ are locally satisfied.
Equilibrium, compatibility along the interfaces r and the interface constitutive law are given a weak
format. This choice is motivated by the kind of spatially discretized equations which will be obtained in
the following Section 4, in which interface variables only will appear explicitly.
l Compatibility conditions on r are given by:
(w-u+ +u-)T&ddT=O V6t, (24)
l Equilibrium is derived from the virtual work principle which in the present case reads:
I R’
uT(e(u))SE(U)do’ -
I r-
tT 6u dr +
I rf
t=Gu dT
(25)
- bT6udOn’- tT Su dS = 0, u E U, V 6u E U,,,
where:
(264
(26b)
In Eq. (25) it has been implicitly assumed that E depends on u through Eq. (1) and that t=t- = -t+
on r; regular in (26) implies continuity and differentiability up to the order required by the previous
relations, in particular E =E(u). Having assumed a linear elastic behaviour in 0’, one has: u(E(u))=
&Y(U), with E the matrix of elastic moduli.
Relation (25) gives the elastic solution of the problem in 0’ for assigned tractions t on r. In the
general case, tractions t are unknowns, and the solution can only be obtained after introduction of the
interface behaviour.
A meaningful interpretation of problem (24)-(26) can be obtained by considering the displacement
jumps w as assigned quantities. Having introduced a linear behaviour in a’, the complete solution in
terms of the unknowns: displacements u in 0’ and tractions t on r, can then be obtained by the
superposition of two sub-problems.
pl-Perfect body subjected to given loads 6 and L (superscript P)
Find u’(x) in a’, t’(x) on r which satisfy (24)-(26) for given external loads and zero displacement
discontinuities w on r.
p2-Body subjected to distortions w. (superscript w)
Find U”(X) in 0’, t”(x) on r which satisfy (24)-(26) for zero external loads and given displacement
discontinuities w on r.
The decomposition above will be used in Section 4 to derive a discrete formulation in which interface
variables only will appear.
l Constitutive Eq. (9)-(11) for the interface behaviour are here replaced by the following variational
statements.
SW”dT = 0, V 6~“; 8adT=O,V8a, (27)
Ir(*pT&+iyT&)d~c0,V&,6aEE,, E~“{6t=t*-t,sa=a*-a~f(t*,a*)~o}, (28)
where St=&(x), &z=&(x), xfZr.

Relations (27) are the weak form of state Eq. (9); relations (28) represent a weak form of the
maximum dissipation inequality (12). The last, in turn, is a restatement of the flow rules (10) and of the
loading-unloading conditions (11).
In order to obtain the weak form of the flow rule equations which will be directly used for the finite
element discretized version of the problem, relations (28) are rewritten as substationary conditions of
the following functional:
L(t, Q, i) = - Jr,
(tTi p+ aTir) dr + I, i’f(t, a) dr,
defined for i~kl={&IisO on r}.

A sub-stationa:ry point of L is defined by (t, a, A) such that:
(30)
(31)
4. Space discretimtion
4.1. Generalized variable formulation
The weak formulation presented in Section 3 is used here as a basis for space discretization by means
of finite elements. The relevant field variables, like t and w, a and a, conjugate in pairs, are replaced by
the interpolations of the corresponding nodal quantities in a way which guarantees the conservation of
scalar products. The discretization variables can then be considered as generalized variables in Prager’s
sense [29,31].
With reference to, e.g., the couple of conjugate fields t and w, denote with T and W the interpolation
parameters for the whole system and with M, and M, the matrices collecting the relevant interpolation
functions:
t(x) = M,(x)T; ~(4 = M&W, (34
where only the dependence on spatial coordinate x and not on time has been made explicit.
T and W can be interpreted as generalized variables if the following relation is satisfied
(33)
Eq. (33) implies that M, and M,,, are forced to satisfy the orthogonality condition
I MT(x)M,,(x)
r
dr = I.
From relations (33) and (34), T and W can be interpreted as weighted averages of the corresponding
local fields:
T = r M;(x)t(x) dT; W = r M;(x)w(x) dT. (35)

I I
The generalized variable formulation has been applied recently to elastic-plastic continua, see
[33-371; the orthogonality relation (34) representing the key-point of this kind of formulation can be
satisfied by different choices of the interpolations of conjugate fields, see e.g. [35].
For the purposes to be pursued in the present paper, it is also meaningful to consider the fields of the
static and kinematic internal variables a and (Y as given by the interpolation of the corresponding nodal
parameters A and (Y through the functions collected in matrices M, and M,.
According to the weak formulation of Section 3, the displacement field u is the only one to be
modelled in 0’. This means that in a’, here assumed to behave linearly, a standard displacement-based
formulation is considered. Hence:
U(X)= N, (x)U in 0 ‘, (36a)
t(x) = M,(x)T; w(x) = M,(x)W on r, (36bc)
a(x) = M, @)A ; a(x) = M,(x)cY on f. (36d,e)

The interpolation introduced by Eqs. (36) are thought to be expanded to the whole relevant domains
0’ or I’, hence assembling procedures have been implicitly assumed.
It is important to notice that displacement interpolation (36a) in 0’ has to be considered as the result
of continuous interpolations in the portions of 0’ separated by the interfaces. It is discontinuous on the
interfaces where nodes and interpolation for sides l-‘+ and r- do not necessarily coincide, moreover the
vector U contains distinct nodal parameters for sides r+ and r-.
4.2. Compatibility
The compatibility condition in a’, given by Eq. (l), has been assumed to be a priori satisfied in the
weak formulation of Section 3. Hence the modelled strains E(X) are related to the modelled
displacements U(X) through the usual relations:
E(X) = B,(x)U; B,(x) = W”(x) in 0’. (37)
Compatibility conditions on do,, Eq. (3), are here assumed to be homogeneous and to be satisfied by
the choice of the degrees of freedom collected by vector V. Non homogeneous conditions can be
introduced in the numerical model in a standard way; their introduction does not affect the generality of
the subsequent developments.
Compatibility on the interface r is given in weak form by Eq. (24). Then, substituting in Eq. (24) the
modelled fields t(x) and w(x) given by Eqs. (36b,c), one obtains:
6TT M~(x)(i14w(x)W-Nu(x)~r+U + N,,(x)l,-U) dT = 0, V 6T.
In Eq. (3% K($],-+, N,(X)],- represent the matrices of displacement shape functions computed on
sides r’ and f - respectively. Taking into account the orthogonality condition (34), the above relation
(38) is equivalent to:
Relation (39~) is the discretized version of compatibility condition (2) on r. It is worth noting that the
above weak statement of interface compatibility allows different discretizations of sides r+ and r- of
r, and also different displacement interpolations on parts fin+ and 0- of R adjacent to r (see Fig. 9a).
Moreover, in the general case, the choice of interpolation functions and sampling points for
displacement discontinuities w along r is independent from that of the displacements U.
Fig. 9. Discretization of the interface: a) general case; b) same discretization on r’ and r-. Circles: sample points for
displacements. Crosses: sample points for displacement discontinuities.
As a particular case, frequently adopted in the literature (see e.g. [23,48]), the same discretizations
and displacement interpolations are assumed on sides r+ and F- (see Fig. 9b). In this case the
compatibility rela.tion (39~) on r can be given the following expression:
w= [i, ikf~(x)&(x)l, d++ - C-)U = [I,- ~~(~)~Jx)~, d+Jr+ - u’-). (40)
In the above relation, tiU(x)I, denotes the matrix which collects the interpolation functions on r, which
are in this case the same for sides r + and r -, reduced to the dimension of degrees of freedom on I’
only. C+ and C‘- are the Boolean connectivity matrices which extract, from the vector U of global
degrees of freedom, those pertaining to the two sides of the interface r, denoted Ur+ and Ur-
respectively.
As a further particularisation, sampling points for displacement discontinuities are made to coincide
with those for displacements. Hence interpolation for w coincides with that for displacements,
computed on r, i.e. w=~$‘~(x)[, W. From the orthogonality condition (34) and relation (40) it is
therefore deduced that:
W= ur+ - ,IJr-*
(41)
4.3. Equilibrium
The discrete form of equilibrium is obtained by introducing Eqs. (36a,b), (37) in the virtual work
principle (25):
(42)
The above statement is equivalent to the solution of the following system of equations:
KU+C;T==P, (43)
where K is the elastic stiffness matrix, C, is the interface compatibility matrix defined by (39) and P is
the equivalent nodal loads vector.
Matrix K can lbe given different properties depending on the set of interfaces r. When r does not
divide the body LJ into completely disconnected parts, K results as being non-singular (provided that
rigid-body motions are excluded) and hence can be inverted as a whole. When the body R is split into
different parts by the set r, still K can be assembled in a block-diagonal shape, each block being
invertible.
Properties of matrix K suggest different ways of reducing the system of Eqs. (39) and (43) to the
degrees of freedom T and W on the interface r, as described in the following section.
340 G. Bokon, A. Corigliano I Comput. Methods Appl. Mech. Engrg. 140 (1997) 329-359
4.4. Reduction to interface degrees of freedom
Let us assume the invertibility of matrix K in Eq. (43). Analogously to what was done in Section 3,
Eqs. (39) and (43), which are the discretized counterparts of relations (24) and (25), can be solved for
assigned P and W by superposing the solution of the following two sub-problems:
Pl-Perfect body subjected to given loads P.
c&J p = 0; KVP+ C;T'=p. (440)

P2-Body subjected to distortions W.
crvw=w; KVW + C;TW = 0. (450)
By eliminating unknowns VP and VW from the problems Pl and P2, vector T can be given the
following expression
T=TP+TW, (46a)
TP = K&K-‘P; TW= -K,W, (Mb,c)
K, = (CJ-‘C;)? (46d)
The set of Eqs. (46) can be interpreted as the version of problem (39), (43) reduced to the interface
variables T and W only.
For the particular case already introduced in Section 4.2, which assumes the same discretizations and
displacement interpolations on sides f + and r-, it is interesting to give a more explicit form of
relations (46).
To this purpose, it is convenient to derive the expression of K, and TP according to the following
operative scheme.
i) A partition of degrees of freedom in vector V is introduced:
v~ = {C@ VrTj; UrT = {p+T Ur-Tj,

(47)
where Vr+ and Vr- collect the 12,.d.o.f.s of corresponding nodes on r+ and r-, see Fig. 9b, while the
remaining nR, d.o.f.s are collected in the vector Va’. The matrix K is accordingly partitioned.
The interface compatibility matrix C,, taking into account Eq. (39~) and (41) acquires the following
expression:
Z(“rxnr) - $IrXtlr)l’
c, = [O(nj-Xnn*) (48)
where 0 is the rectangular zero matrix and Z the identity matrix of the dimensions noted in (48).
ii) D.o.f.s Vn’ are eliminated through static condensation from Eq. (43), thus obtaining:
w= Q.vr; RVr+d;T=B, (490)
er E [Zcnrxnr) -z~+~n,jl; f = K,, - KrO,K,f,sK,.,; p = Pr - K,,.K,:o.Po’. (49cAe)
iii) As already done for the complete ptoblem,_Eqs. (46), the solution to problem (49) is interpreted
as the superposition of two sub-problems Pl and P2 which are defined analogously to Pl and P2 in Eqs.
(44) and (45):
vr = vrP + vrW,
[;::]=[;:]+[:;J* (500)
In the above equations, fi’_is the common value of displacements Vr+ =Vr- in the case of perfectly
bonded interfaces (problem Pl); this common value V can also be interpreted as the mean value
(Ur++Ur-)/2. The expression given to UrW in Eq. (50b), i.e.: Urw=~~W12, follows from the
introduction of if’ = (U + + Ur-)/2 into (5Oa), from the comparison with Eq. (49a) and from the
definition (49~) of matrix e,-.
Pre-multiplying Eq. (49b) by the matrix e,-, substituting the expression for Ur given by Eqs. (50) and
taking into account that @?f =2ZC,,+,,,), it follows that:
+2T=&fi. (51)
Hence the following relations, analogous to (46), are obtained:
T=TP+TH’, (524
TP = ;(er$ - tfr&J”); TW= -K,W, (52b,c)
The expression (5#2d) of the matrix K, is computationally advantageous with respect to that obtained in
the general case and given by Eq. (46d). It implies only the inversion of matrix Knrns, see Eq. (49d).
The reduced version of relations (39~) and (43) given by Eqs. (46) or by Eqs. (52), represents the
direct link betwee.n the interface variables T and Was a conseyence of elasticity of 0’ and equilibriym.
Matrix Kr reflects the elastic stiffness of the body 0’; vector T , solution of the sub-problem Pl (or Pl),
depends linearly on the given external loads P.
The complete solution to the discrete version of the problem governed by the equations introduced in
Section 2.1, can then be obtained after supplementing relations (46), or (52), with a discretized version
of the interface constitutive law.
4.5. Interface constitutive law
The discretized version of the interface constitutive law is obtained after introducing in Eqs. (27),
(30), (31), the interpolation of conjugate fields t, w and a, (Y given by relations (36) and the
interpolation of plastic multipliers i=M,A. The additivity relation (8) is preserved in terms of nodal
values, i.e. W=W”+W’, the same interpolation is used for w, we, wp.
8WeTlr M:($)(M,(x)T - -$(we(x), a(x))) dr = 0, V &We, (53)
SCYT/r M:(x) (&f&$4 + -$w~(x), a(x))) dT = 0, V 6a!, (54)
6TT M:(x) M,(x)~#~ -g@(x), a(x))M,(x)~ dr = 0, V 6T, (55)

Ir ( >
&AT MT(x) M,(x)& - g(<x), a(x))M,(x)A dF = 0, V &A, (56)

Ir ( >
sAT M;f(x)f(t(X), a(x)) dT G 0 V S/i EM, M={6A=A*-/iIA*>o}, (57a,b)

Ir
where the arguments in functions $(w”(x), a(x)) andf(t(x), a(x)) are meant as interpolated quantities.
Relations (54)-(57) can be rewritten in a way formally analogous to that governing the interface
constitutive law locally (relations (9)-( 11)). For this purpose let us first introduce the functions p and
F, defined as:
ty(We, (u) =s1, JI(M,(x)W’, M,(x)(Y) dF, (58)

With the above definitions and by making use of orthogonality relations between interpolations of
conjugate variables analogous to Eq. (34), it can be seen that (53) and (54) are the discrete version of
(9a) and (9b), respectively. Relations (55)-(57) are in turn equivalent to the discrete version of the
maximum dissipation principle, and therefore equivalent to a flow rule, for the discretized system,
formally analogous to the local flow rule given by Eqs. (10,ll). This can be seen by reasoning, as was
done at the end of Section 3.
On the basis of the above considerations, the discrete (in space) version of the interface constitutive
law can be written as follows:
w=we+wp, (60)
T = aF(we, a). A = _ fww”7 a).

da ’ (61)
awe ’
(62)
A 30; F(T, A) s 0; FT(T, A)i = 0. (63)
The above set of relations represents the interface behaviour in a global average sense; this can be
considered as a multi-mode constitutive law. Together with relations (46), relations (60)-(63) govern
the problem of an elastic continuum with damaging interfaces discretized in space and reduced to
interface unknowns only.
It is important to remark that, in order to have a global constitutive law which reflects the properties
of the local one, it must be guaranteed that the properties of convexity (or lack of convexity) of the
local energy function JI and that of convexity of the local activation functions f are transferred to the
global functions !P and F. Due to the definition of ?P, it is straightforward to show that the local
character of II, is transferred to ?P. As far as the definition of F is concerned, one can observe that
convexity is surely transferred to the global level if the interpolation of A is such as to give always
positive MA(x) for any x along r. If this is not the case, convexity of F should in any case be guaranteed
by taking care of interpolation for A. For more details on this and other subjects concerning the
generalized variable formulation (see e.g. [32-351).
5. Uniqueness of the rate problem discretized in space and stability of equilibrium states
It is well known that uniqueness and stability are not guaranteed in the presence of damaging
phenomena, such those attributed here to the interface behaviour which shows softening as a peculiar
feature [3-5,25-27,38-44,49-511. In this section uniqueness of the rate problem and stability of the
equilibrium states are discussed.
The reduced formulation in terms of interface variables, obtained in Section 4, is particularly
convenient for this study. This can be easily recognised since it has been assumed that all non-linearities
and irreversible phenomena, which can cause non-uniqueness and/or instability, are concentrated along
the interface.
For given increments of the external loads i, starting from a completely determined situation, the
rate problem in point is governed by the rate form of Eqs. (46):
iYK,ti’==‘; P’ = K,C,K-‘P, (64)

and the rate form of Eqs. (60-63).
For assigned 1”, tractions tP can be uniquely determined from Eq. (64b), hence in the rate problem
vector fp must be considered as a given quantity.
G. Bolzon, A. Corigliano I Comput. Methodr Appl. Mech. Engrg. 140 (1997) 329-359 343
In order to derive uniqueness and stability conditions, the classical path of reasoning proposed by Hill
[38], is first followed in Section 5.1. In the subsequent Section 5.2 uniqueness and stability are
addressed following an approach called here intrinsic which dates back to the works of Maier [25],
Maier and Drucker [39] and more recently of Nguyen and co-workers [40-441.
5.1. Uniqueness and stability by the classical approach
In order to derive a sufficient condition for uniqueness, let us suppose that two distinct solutions (I#‘i,
f,), (w2, f,) exist for the same f’. From Eq. (64) it then follows:
AW’A? + Alir’K,AI$‘= 0, (65)

where the symbol A denotes the difference of quantities pertaining to the two solutions 1 and 2. From
Eq. (65) it is deduced that a sufficient condition for uniqueness of the solution in terms of i‘ and w is
that
Aw=A? + Aw=K,AI+ > 0, V Arir = w, - I+$ # 0. (66)

When inequality (66) is satisfied, Eq. (65) implies that A*=O, i.e. w1 =I$‘*; from Eq. (64a) uniqueness
of f for given f” follows.
In the above inequality (66), the quadratic form AriVTK,.Aw is a positive term; this can be seen by
recalling the expression given to matrix K,, Eq. (46d), and by considering that a linear elastic behaviour
has been assumed1 for the portion 0’ of the body. The only term which can cause non-uniqueness in the
solution of the rate problem is then the scalar product ApTA? which in turn depends on the interface
behaviour.
A sufficient condition for the stability of an equilibrium state is the non-negativeness of the second
order work performed by an external agency which preserves equilibrium between internal and external
forces, in any kinematically admissible configuration change infinitesimally close to the initial one,
complying with the constitutive law.
In the present context, the second order work for the space-discretized system is given by (cp. Eq.
(43)):
2fj2Z = GUTK 6U + @CT, 6T = GUTK NJ -i-6WT 6T, (67)
where use has been made of the interface compatibility condition, Eq. (39).
It is now possible to take advantage of the decomposition introduced in Section 4.4 in order to derive
an alternative expression of the second order work which allows for a meaningful interpretation and
which gives some advantages for numerical applications. The kinematically admissible configuration
change Xl can be thought to consist of the sum of two addends NJ’ and SU w such that
6U = 6UP + 6UW, (68a)
SUP = K-l@- C;K,-C,K-‘) SP, (68b)
6U w = K -‘C;K, SW, (68c)

where matrices C, and K, have been defined previously by Eqs. (39~) and (46d) respectively.
The introduction of the above relations (68) in the expression of the second order work (67) gives:
2ij2Z = GUpTK 6Up= + 2 GUprK MJw + SWTKr 6W + SW= 6T. (69)
The second adadend in the r.h.s. of relation (69) can be shown to be zero by making use of relations
(68b,c) and of the definition (46d) of matrix Kr. The same conclusion can be obtained by applying the
principle of virtual work with static quantities pertaining to the sub-problem P2, relations (45), and
kinematic quantities pertaining to the sub-problem Pl, relations (44). The expression of s2Z’ thus
changes as follows:
2S2Z = SU PTK SU pT + SW’K,- SW + SW’ ST.
The above expression of S2Z shows completely de-coupled contributions of kinematic disturbances SUP
and SUw. The first addend in the r.h.s. of Eq. (70), g e rt aining to disturbances SUP, is found to be never
negative and to be equal to zero only for zero SU . Hence, as disturbances SUw, linked to non-zero
displacement discontinuities SW, are considered, it follows a sufficient condition for stability of the
discretized system in point in the form:
SWT ST + SW=K, SW > 0, V SW # 0, (71)

where disturbances SW are kinematically admissible, and quantities SW and ST comply with the
interface constitutive law.
By comparing the sufficient condition for uniqueness expressed by relation (66) with the above
sufficient condition for stability, one can observe that the uniqueness condition is stronger than the
stability one, i.e. uniqueness implies stability but the opposite is not true. This is classical when non
uniqueness and stability depend on incrementally non-linear constitutive behaviours, like in Shanley’s
column [49].
It is also worth noticing that the remark made with reference to the uniqueness condition (66),
concerning the positivity of the quadratic form based on matrix K,, still holds. This means that
instability is caused by constitutive softening only (the interface constitutive behaviour considered in the
present work is of an associate kind).
By examining conditions (66) and (71), it clearly appears that they require the evaluation of the
terms AWrAp and SWTST.
In the most general case of an incrementally non-linear interface behaviour, a concept similar to that
proposed by Hill of a linear comparison solid must be introduced. This allows us to obtain over-
sufficient uniqueness and stability conditions, which can be practically applied.
Let us define a linear comparison interface as an interface r which is governed by a linear incremental
relation based on matrix HL, i.e. f= HLw, which satisfies the following inequality:
AW’AP 2 AWTH,AW, (72)

for any W1, I$‘; solution of the rate problem, with ?, , f”, linked to them via the incrementally non-linear
interface constitutive law.
From the above definition it follows that uniqueness and stability conditions can be replaced by the
following over-sufficient conditions.
A sufficient condition for the uniqueness of the solution of the rate problem in terms of i’ and W is:
AW=(H,+K,)AW>O, VAW=riv,-W2#0, (73)
where W1 and W2 must be admissible displacement discontinuity rates for the incremental problem.
A sufficient condition for the stability of the equilibrium state is:
SW=(H, + K,) SW > 0, V SW # 0, (74)
where SW is kinematically admissible.

The above conditions are both guaranteed when matrix (HL + K,) is positive definite.
For the derivation of matrix HL, which must comply with property (72), it is important to notice that
in the present generalized variable context, matrix HL must be derived by examining directly the
incremental constitutive law which applies at the global level and from which the tangent matrix of the
non-linear incremental process can be derived. In other words, in general it is not possible to derive
inequality (72) from an analogous one which holds point-wise locally.
In the particular case in which interface interpolations are based on Gauss integration points, the
interface law de-couples and can be imposed independently at the Gauss point level, see e.g. [35]. This
allows the deduction of inequality (72) from an analogous one which holds locally, i.e. in this case the
linear comparison interface is defined in the same way as at the local level.
G. Bolzon, A. Corigliano I Comput. Methods Appl. Mech. Engrg. 140 (1997) 329-3.59 345
5.2. Uniqueness and stability by the intrinsic approach
As shown in Section 5.1, sufficient conditions for uniqueness and stability can be derived involving
the interface variables T and W only, after condensing the remaining variables pertaining to the solid 0’
outside the interfa.ce r.
A further reduction is possible: it consists of the elimination of those variables which describe the
reversible (elastic) part of the interface constitutive law. One can thus obtain uniqueness and stability
conditions which involve the kinematic internal variables Wp and CY only.
This kind of intrinsic approach has been followed by Maier and Drucker [39] with reference to
elastic-plastic structures and more recently by Nguyen and co-workers [40-441 in a series of papers
concerning a large class of irreversible phenomena governed by standard dissipative constitutive
equations. This approach is followed here in order to obtain an alternative, explicit expression of
uniqueness and st.ability conditions.
Let us first rewrite explicitly the relations governing the rate problem and recall some meaningful
properties of the :Row rule.
i-= -K@-ki-‘, (75)
(76)
f = a*W(W”,
a) *’ + a2qwe, a!) &
(77)
aWe ffWeT xr aweT ’
A = _ -~
a2qwe, Ly)& _ a2?qwe, cqw,
(78)
aa! ac.xT awedaT ’
wp _ aFT(Ty
AjA. iw = aFTPWA
aT ’ aA ’ (79)
A 30; F(T, A) s 0; FT(T, A)li = 0. (80)
It is worth recalling that, due to the particular class of standard dissipative interface constitutive laws
considered here, the flow rule (79)-(80) satisfies the maximum dissipation principle, i.e. :
(T*-T)TI+p+(A*-A)TiY~O; V(T*,A*)EE, E = {T, AjF(T, A) s O}. (81)
It can be shown (see e.g. [40]) that for this class of dissipative systems the flow rule in rate form
satisfies the following relations:
PTWP + A Tth! = 0, (82a)
Lix*1
*i.lp*
fT@* +i:r&f* 60, v E N,(T,A), (82b)
where N,(T, A) d,enotes the normal cone to the elastic domain E:
NC(T,A)={~p,~~(T*-T)T~p+(A*-A)T&~O;V(T*,A*)EE}. (83)
From relations (82) it is deduced that the incremental process satisfies the following variational
statement
W
E N,(,T, A), (84a)
[ du I
1
J,pyP*
(W- w*)=f + (iu - &!*)‘A 20, v E WT, 4.
[ ti!*
Another way to characterise the incremental process is to make direct use of relations (79)-(80). For
the active modes at the current time instant, i.e. for those k modes collected in vector #, such that F).(T,
A)=O, i=l.. . k, the flow rule can be written
By means of Eqs. (75-78) and of the incremental flow rule, alternatively and indifferently expressed
by Eq. (84) or by Eqs. (85), (86), it is now possible to reformulate the whole rate problem after the
elimination of the reversible variables rir”. This is done under the hypothesis that the following
invertibility condition holds
a2w
sw>o, vswzo. (87)
awe aweT>
For the present purpose, k’ is expressed as a function of wp, du and Pp by means of relations
(75)-(77) and the result is substituted in relations (75) and (78). This transformation allows us to
obtain the rate of static variables F and k as functions of the rate of the kinematic irreversible variables,
tip and &, and of the input quantity fp:
i’ = - (GWPWPI%”
+ GaWPiy) + BTfP, (88a)
k= -(G WP,wp + G,,iY) + Bjp, (88b)

where the following definitions have been introduced:
a2vf -’
G wpwp =K,- -& K, + K
aweaweT > r’
aZF -’ a2*
G cxwp = -K, K,+ = G&m, @%c)
awe aweT > acuaweT
(894
(894
BAp _ a2q 0-W

awe actrT
For the sake of brevity, let us introduce the following compact notation:
(90)
noticing that matrix G defined above is symmetric (cp. Eq. (89~)).

The introduction of relations (88) in the variational inequality (84) by making use of the compact
notation (90), transforms the rate problem in the following way:
i-f=X0, (91a)
(F-F*)'(ci--st')'~o, vi-* ~zv,(x). Plb)

Existence and uniqueness conditions for the solution(s) of (91) can be found in the literature, see e.g.
[401.
The existence o:E at least one solution is guaranteed, for arbitrary fp, if
FTGF>O, Vi,#OEN,(Z). (92)
The uniqueness of the solution in terms of i- is guaranteed when
Vi;#i;EN,(X).
(i;-i72)TG(,i;-i;)>0, (93)
Existence and uniqueness conditions can be further transformed by substituting the condition
peN,(Z,) with a Imore explicit one. This can be indirectly done by making use of relations (85), (86)
characterising the incremental process. The introduction of Eqs. (%a), (88b) (with the compact
notation) in relations (85), (86) transforms the rate problem into:
Jj= _ ($)T+!& +(i!!?)TBj’p<O;

$‘_&O. (94)
The above set of relations (94) represents a Linear Complementarity Problem, for which existence
and uniqueness conditions have been stated, see e.g. [25].
A unique solution to problem (94) exists in terms of A, for any fp, if and only if the symmetric
matrix
(g)‘G(!$),
being matrix G symmetric, is positive definite.
The number of possibly multiple solutions is finite for all f Pif and only if all the principal minors of
the matrix
are non zero.

The necessary and sufficient uniqueness condition in terms of R transforms in a sufficient uniqueness
condition in terms of f, because i- is uniquely defined, through (79), by /i.
As has been done above for uniqueness conditions, stability is now discussed by making use of the
alternative formulations of the rate problem given previously.
In Section 5.1 it has been shown that the second order work S2.Z * due to disturbances 6W can be
expressed by (71). Then, from (75) it also follows that
s2.Y* = W”j’ + WT&W = WTpP,

2- (95)
a27
where 7 denotes time.

Eliminating W’ as already done in the derivation of Eqs. (SSa), (88b), it is now possible to rewrite
(95) as
s22?*
2 x= i”+ pTBfp. (96)
Relations (88) allows us to express Bfp as a function of 2 and p; by exploiting this relation and Eqs.
(82a), Eq. (96) further transforms into:
(97)
The first addend in Eq. (97) is non-negative due to (87), hence a sufficient condition for stability,
based on the second order work criterion, can be reformulated as:
pTGp>O, VFWJX). (98)
Again, the above condition can be made more explicit by making use of the flow rule (85), (86).
Relation (98) thus transforms into:
(99)
which means co-positiveness of the matrix involved.
6. Examples
In this section some examples are presented and discussed. These examples are only intended to
better clarify what has been presented in the previous Sections. Numerical results concerning
applications of the proposed approach will be presented in a forthcoming paper.
In Section 6.1 two elastic beams connected through an interface r along their length are considered.
This is done in order to focus on some aspects of the generalized variable discretization and to explicitly
show an example of a possible interface interpolation. It is also shown how the global interface
behaviour results in being de-coupled with respect to the contribution of each field point if the
Newton-Cotes rule for numerical integration is adopted.
In Section 6.2 two meaningful cases of interface models, formulated in terms of generalized variables,
are briefly commented upon; uniqueness and stability conditions are derived for them on the basis of
the results presented in Section 5.
In Section 6.3 two examples concerning a frame and a beam structure with softening hinges are
considered. This is done in order to show the perhaps most classical case of elastic structures with
damaging interfaces, dealt with in the framework of generalized variables. This kind of structural model
dates back to the end of the sixties and has been extensively studied, see e.g. [25-28,49,50]. The two
examples are completely solved, and a discussion on possible bifurcation and instability phenomena is
presented which shows the complexity that the structural response can reach also for few degrees of
freedom.
6.1. Beam model of layered composites connected through a damaging interface
Two straight elastic beams connected along their sides through a damaging interface are considered
as a simple model of layered composites which can undergo delamination (Fig. 10a). Perfect adhesion
between the two beams before decohesion can be simulated through the elastic part of the interface
constitutive law. Delamination, i.e. the progressive separation of adjacent layers, is described through
the irreversible, softening behaviour of the interface [16,17].
The Timoshenko kinematics is chosen for the beam structural model, and a beam finite element with
linear interpolations for displacements and rotations is considered. It is worth noting that in the present
formulation the beam generalised displacements are those of the lower side for the upper beam (cp. F+
in Fig. 10a) and those of the upper side for the lower beam (cp. r- in Fig. 10a). This choice differs
from the classical one in which generalised displacements coincide with those of the beam axis, and
allows for an easy introduction of the interface constitutive law since in the numerical model the upper
X r+
_________________________-_____________
(4 I-- I-
I
1
I
+
(b) ~2_.;__________~_________~___-_-_
0
r-O
Fig. 10. Beam model of layered composite: a) schematic representation; b) discretization.
and lower beama are substituted by lines which coincide with interfaces r+ and r-. The same
interpolation is chosen for the upper and lower beam elements, with nodes having the same coordinate
along axis x (Fig,. lob). Moreover, sampling points for displacement discontinuities are chosen to
coincide with those for displacements on r+ and r-.
Due to the above hypotheses, the interpolation for displacement discontinuities w(x) coincides with
that for displacements u along the sides r’ and r-, i.e. M,(x)=NU(x), as shown in Section 4.2. The
interface compatibility matrix is simply given by Cr=(C+ -C), where C+ and C- are the Boolean
connectivity matrices which extract from the global displacement vector U those degrees of freedom
pertaining to the interface, excluding rotations.
The following expression for the modelling of displacement discontinuities can then be deduced.
1
N 1" ..Ni.. ..N,,, 0.. ..O.. ..O
wT=[...(u+-u-)i...
...(7J+-u-)i...]
M,=
i '0.. ..o.. ..O Nl.. ..Ni....N,,,
; (100)
In the above relations: i = 1, . . . nn is the index which denotes the discretization nodes along x; Ni
denotes the shape function relevant to node i expanded to the whole axis x as shown in Fig. 11.
As discussed in Section 4.1, in order to satisfy the requirements of the generalized variable
discretization, the interpolation for tractions t must be chosen such that to satisfy the orthogonality
condition (34). A. simple way for doing this is to chose M, as:
-1
M,(x) = M,(x) dr
> . (101)
With the above choice, and with analogous choices for internal variables (Y and a, the complete set of
equations governing the structural problem can be recovered.
In Section 4.5 it has been observed that the global interface law resulting from the present approach
is generally coupl’ed, i.e. it cannot be obtained by imposing point-wise the local interface law. However,
the de-coupling can be performed when the field points used in the interpolation of those variables
which are relevant to the constitutive law, coincide with the sampling points of the numerical
integration (see e.g. [32,35]).
In the present context, the complete de-coupling is obtained if the integration points coincide with
the nodes of the discretization on interface r. For the two-nodes elements here chosen, this is
equivalent to adopt the Newton-Cotes or Lobatto integration rules along each finite element (see [48]).
The integral along r of a function g(x) is then computed as:
i-l i i+l
Fig. 11. Expanded shape function for the i-th node.
(102)
where L, and L, are the length of elements before and after node i, and wi = 1, i = 1, 2, are the
corresponding weights.
The numerical integration of the r.h.s. of Eq. (101) through rule (102), gives matrix M, of the form
1
..NilOi.. ..N,,,IG,, 0..
.. .
..O.. ..O
(103)
0 0 N,/W,.. ..N;lWi.. ..N,,,IW,,
The introduction of interpolations (100) and (103) in Eq. (35), shows that variables collected by
vector T acquire the meaning of local nodal values oft multiplied by the corresponding nodal weight Gi,
while the variables collected by W are simply the local nodal values of w.
The numerical integration of the global energy function ?P and of the global activation functions F
defined by (58) and (59), gives
(104b)
The above expressions of ‘I’ and F show how the constitutive behaviour is completely decoupled node
by node; therefore, the properties of local stiffness matrices are directly transferred to the global level.
In particular, the matrix H,_ of the global linear comparison interface can be built by simply assembling
the stiffness matrices h, of the local linear comparison interface.
6.2. Uniqueness and stability for two particular generalized variable interface laws
As a first example, a non-linear elastic behaviour for the interface is considered. As shown in Section
2.2, this is defined at the local level by the function ?P(w) only, together with relation (9a). The
incremental behaviour is hence linear and is given by:
~“wu~ a2wV = a2*

p=
awawT ; aWawT- Ir
MT
WawawT
M dr
w ’
(105)
where use has been made of the definition (58) of the global energy function q(W).
In this elastic case, both uniqueness and stability conditions, Eqs. (66) (71), require the positive
definiteness of matrix (K, + (a2?P/dWaWT)).
It is worth noticing that, in the previous discussion of Section 5.2, the above matrix has been assumed
to be always positive definite in order to concentrate on the effects of the irreversible interface
behaviour.
As a second case, the cohesive-crack model is treated. For pure mode I, the constitutive behaviour is
defined by relations (8), (9b), (lo), (ll), (13), (14), and has the character of a single-mode
rigid-plastic model, with the unique kinematic irreversible variable w.
By their definitions, Eqs. (58), (59), the global functions ?P and F acquire the format
(106)
F = r M;(M,T + M,A) dT = (j-$M,dT)T + (lrM:M,dT)A=R:T +R;A. (107)

I
The above expressions (106), (107), with relations (60), (61b), (62) and (63), define the interface
behaviour. The structural problem in point is then completed by Eqs. (46), expressing the overall
equilibrium.
The above global constitutive law does not apply to the portions r, and r, of the interface r which
are either undamaged or completely damaged. In this framework of a multi-mode global constitutive
law, part r, of r is now defined by the set of inactive modes which satisfy T CT,,, while part c is
defined by the set of modes for which T,+H(Y =O. Notice that, due to the linearity of F with respect to
T and A, the convexity of the local activation function f is directly transferred to the global one.
Uniqueness and stability conditions can be stated by following unaltered the path of reasoning of
Section 5.2, starting from Eqs. (75) (76), (78)-(80), with !3’=ly(a~) in Eq. (78) in the present
rigid-plastic context.
Then a unique solution exists in terms of^A for any fp if and only if the matrix (R;K,R, + RiHRA)-
is positive definite, where the superscript means that the matrices involved must be built by taking
into account the active modes only; a sufficient condition for stability of the equilibrium state requires
matrix (RGK,R, + RiHR,)* to be co-positive.
If, as a particular case, it is assumed that interpolations M,, M,, M,, M,, M, satisfy the following
relations (notice that, locally, t, a, h, (Y, w are all scalars in this case),
W (4 = M, (4 ; M,(x) = M,(x) = M,(x);
matrices R, and R, reduce to the identity. The uniqueness and stability conditions for this last case
have been recently derived by Bolzon et al. [51].
6.3. Frame and bleam with softening hinges
As a first example, the framed structure of Fig. 12 is considered, with the geometrical characteristics
reported in the picture, where J denotes the cross section moment of inertia. The elastic solution for
given load, in terlms of the bending moment distribution, is also drawn in Fig. 12. Parameters r and k
are defined as: r==(3+2k)/3; k=.I,alJ,b; k>O, r>l.
Sections denoted by r,, r, and r, in Fig. 12 are considered as critical, in the sense that plustic-
softening hinges can there develop when the elastic limit M, (or -M,) for bending moment is reached.
The load multiplier q is parametrized on the value PE of the load at the elastic limit for the structure;
this is related to M, by: PE=(M,,Ib) 8rl(2r-1).
The structural problem represented in Fig. 12 belongs to the class of problems studied in this paper.
Critical sections j’, i = 1, 2, 3, represent the set r of damaging interfaces defined in Section 2. The
behaviour of each softening hinge (interface 4) is given in terms of the bending moment M and relative
rotation 0 which play the role of traction t and displacement discontinuity w, respectively. Notice that
mixed-mode interface laws should be considered if interaction between axial and flexural behaviour
Fig. 12. Framed structure and its elastic solution.

should be taken into account. Notice also that, due to the discrete number of critical sections, the global
vectors T and W are naturally defined as those collecting the relevant moments and relative rotations.
As far as the elastic behaviour is concerned, the application of the beam theory gives analytical
solutions to the two sub-problems Pl and P2, introduced in Section 4.4, which allow for the reduction
to the interface degrees of freedom only.
The interface law linking the bending moments to the anelastic relative rotation in the critical sections
q, is represented in Fig. 13. This kind of model corresponds, in the present context, to the cohesive
crack model of Section 2.2, depicted in Fig. 4; it can therefore be considered as a rigid-plastic softening
model until O< 10 1<MO / I/zI. Unlike the cohesive crack model for mode-1 crack propagation, where only
opening (w>O) causes interface degradation, in the present case both positive and negative relative
rotations must be considered.
An intrinsic sign-convention has been assumed for bending moments and anelastic relative rotations
0. In particular, it has been assumed that bending is positive when the inner fibres of the frame are
under tensile stress, while relative rotations are positive when the internal work is positive for positive
moment. Hence, by defining:
@=(@+-@-)with@+sO, O-20, (109)

a description of the behaviour along the softening branches in Fig. 13 can be given in the framework of
the class of models introduced in Section 2.1. For each softening hinge, the interface model can be
obtained from the following energy and activation functions:
1
Cc,= rn;fa +yThm; f=nTA4+a, (110)
where,
m,,=[tI]; -_[:I]; h=[: ho];I-[:‘]; u-[I:]; n=[l-11. (111)
The interface model is completely defined by supplementing the above definitions with Eqs. (8), (9b),
(lo), (11) formulated in the case of a two-modes plasticity model, taking care of the identifications
M+t; O-+w. The negative parameter h = -MO/O,, governs the softening behaviour.
The constitutive law in generalized variables at the global level can be simply derived due to the
independent behaviour of the interfaces 4. Any integral over r reduces in fact to a sum of quantities
pertaining to each critical section &, without any need of interpolating functions. Then, the global
energy and activation functions defined by (58) and (59) are given by:
3
P(a!)=C~;(ai)=M~a ++a, (112a)

i=l
(112b)
Fig. 13. Interface law for the softening hinge.

G. Bolzon, A. Corigliano I Compur. Methodr Appl. Mech. Engrg. 140 (1997) 329-359 353
wherethe vectors Ac6x1j,Mt3,,, collect the sub-vectors of each softening hinge

M0c6x1j,ac6xlj, Fc6x11,
(interface), while matrices HChX6)and NC3X61are block diagonal.
Functions defined by (112) supplemented by Eqs. (60), (61b), (62) and (63), completely define the
global interface behaviour.
In order to discuss the structural response of the frame in Fig. 12, the solution of sub-problems Pl
and P2 is derived first. This amount to compute the bending moments Mj in the critical sections c as
sum of the effects of the load qPE on the perfect structure and the effect of the relative rotations ei
(internal distortions) on the structure with zero loads. The resulting values are:
M=MS+M”=M”-&@, (113)
(114)
where E denotes the Young modulus of the material.

The evolution of the system can be traced by supplementing the above equilibrium relations with the
interface behaviour defined previously. The structural response will depend on the stiffness ratio r and,
to a greater extent, on the slope h. The whole overall behaviour of the structure is summarised in Fig.
14, where the deflection 6 in correspondence of section r,, where the external load is applied, is drawn
versus the load multiplier 77.The results have been obtained by taking into account the possible critical
situations for uniqueness and stability, and by following for each equilibrium path the evolution of the
structure at varying load multiplier.
Uniqueness ant1 stability of the structural response can be analyzed from the conditions derived in
Section 5. Uniqueness in terms of A and stability of the equilibrium state are-in this example based on
the positiveness and co-positiveness respectively of the matrix (NTKJV+H) where the superscript
means that active modes only must be considered. This condition can be derived as done in Section 6.2
with reference to the cohesive-crack model in generalized variables. By analysing the behaviour of the
structure, it emerges that the critical situations are defined by the following meaningful cases.
In the first phase, only section r, is activ_ated with positive relative rotation Ol, situation E of Fig.
14; the corresponding matrix (iVTK,-N+H) simply reduces to:
(ij’TKrA’+H)A= (z+h) =$(1-h*) (115)
where h* = (hlbrl,EJ,,, and uniqueness and stability of the overall response are guaranteed for O<h * < 1.
Stability is lost for h* > 1, and above some critical level (here h * = 32/ 11) the severe snap-back
instability is met.
In a second phase, as bending moments at sections r, and r, reach their ultimate value, -MO,
softening hinges can potentially activate in these sections asAwell, allowing for rotations 0; and 0;
beside 0:) situation P of Fig, 14. The matrix (NTKJV+H) now assumes the following expression:
(l-h*) 1 -1
(NTKrN+H)A=~ (l-h*) -1 . (116)

-1 (l-h*) I
The eigen-value analysis of the above matrix gives the eigenvalues:
-EJs/r$ .
PI +(3 -h*); &=p3=- ,.b 9 (117a)
with the corresponding eigenvectors

rl
Fig. 14. Overall behaviour of the framed structure and schematic representation of active (black) and inactive (white) softening
hinges.
(117b)
The two coincident eigen-values ~12,k are less than zero for any value of h* #O, and correspond to
unsymmetric rotation modes, revealing the possible bifurcation of the overall behaviour which, as usual
in softening, destroys the symmetry of the structural response; in fact, in the observed bifurcated path,
one hinge activates while the other stops rotating.
It is worth noticing that, because of the discontinuous slope of the moment-rotation diagram, the
eigen-values of the stiffness matrix can become negative without passing through the critical value zero.
Other meaningful situations arise when the inelastic rotation of hinge r, reaches its critical value O,,
beyond which it cannot transfer the bending moment any longer. From this stage, as far as 0: >O,,,
rotations have been assumed to be reversible, see Fig. 13.
This may happen after inelastic rotations of section c and r, have been already activated, situation B
or, as limit case, situation Q of Fig. 14. In this case that the matrix (NTK,-iV+H)A becomes:
(118)
which also has negative eigenvalues corresponding to unsymmetric inelastic modes (one hinge rotates
while the other unloads).
The condition of zero bending moments in section r, allows us to condense the degree of freedom
pertaining to the recoverable rotation 0;. The matrix in Eq. (118) can then be replaced by:
-;* _“h*
(lvTzQv+H)A=~ [ I )
(119)
whose negative eigenvalues are evident.

Some particular situations are further discussed here below.
When h” = h =O, the critical sections have perfectly plastic behaviour. The load multiplier can grow
above the elastic limit (n= 1) until the combined effect of increasing the load and of the relative
rotation 0, (positive) causes sections r, and r, reach their elastic limit; this situation represent the
so-called collapse configuration for the considered structure, and can be recovered by applying the
theory of limit analysis.
When O< h* s 1, it happens that inelastic deformations of section r,, as well as the structure
deflection, can grow until the bending moments of sections r, and r, progressively decrease up to reach
their elastic limit, -44, (notice that for h* =l this happens at constant load). From this situation,
labelled as P in Fig. 14, the evolution of the system could be either symmetric or unsymmetric.
Symmetric response corresponds to increasing inelastic deformations and decreasing bending moments
at section r, and r, simultaneously, while the unsymmetric path foresees one section deforming
inelastically while the other unloads, up to the critical situation indicated as F in the figure; at this stage,
both section r, and r, (or, respectively, r, and r,) have exhausted their load-carrying capacities:
O;=O,=M,llhl; M,=M,=O (or, respectively 0 : = 0 i = M, / 1h 1; M, = M, = 0). Thermodynamical
considerations suggest that this latter path will be followed by a real structure, nevertheless it can be
interesting to note that, if the symmetric path is followed, a further bifurcation point B is met after
section r, is completely unloaded, while sections r, and r, can still transfer some load.
For 1 <h * <2r/(2r - l), the same sequence of phenomena is repeated, but with the activation of
hinges in correspondence of sections r, and r, occurring at lower and lower load, and with steeper and
steeper softening branches as h * increases; h * = 2rl(2r - 1) represents the limiting situation beyond
which section r, completely exhausts its load carrying capacity before other hinges can eventually
activate and points P and B coalesce into point Q.
When h* >2rl(2r - l), section r, loses its capability of transferring bending moments before softening
hinges of sections r,, r, activate. This transforms the original structure in a new elastic one of the
reduced stiffness represented by line OAQ of Fig. 14. The followed path is then OEAQ. When point Q
is reached (still with M, = 0 but now for 0, > 0, = M, / Ih I), hinges form in sections r, and/or r, and
bifurcation is again possible. The final evolutions from this situation, up to exhausting the load carrying
capacity, is represented again by softening branches of slope increasing with h * .
As a second example, the simple clamped beam shown in Fig. 15 is studied, by assuming that inelastic
rotations can activate at sections r, , r, , r, , and r,. The same piece-wise linear interface constitutive law
assumed for the previous example has been chosen.
Stability conditions for this example have already been discussed by Maier [50] restricting the
discussion to the equilibrium paths which result from the imposition of symmetry conditions for the
whole beam.
The overall be:haviour of the structure, after removing the above restriction, is represented in Fig. 16:
the deflection of sections r, and r,, where the load is applied, is drawn as a function of the load
multiplier. For unsymmetric deformations, the mean value of these deflections has been considered.
The load has been parametrized again on the value at the elastic limit for the structure, attained here
for PE = 3M,/2C
Fig. 15. Clamped beam and its elastic solution.
_ Q _ _ F _ G
r t t-3 4 * 4
Fig. 16. Overall behaviour of the clamped beam and schematic representation of active (black) and inactive (white) softening
hinges.
The remarks already made for the previous example still hold; it is interesting to note how the
complexity of the response greatly increases by simply adding one more critical section. In particular,
stability is lost for h* >3 (here h * = 9(h(k’I2EJ) after the first softening hinge forms, and bifurcated
paths can start from the situation at the elastic limit (q= 1) for all values of h*>9 since, in the present
case, softening hinges can contemporary appear at sections r, and r,. In fact, in this situation matrix
(NTK,-N+H)A is given by
(120)
with eigenvalues
p1=$(3-h*); &=%(9-h*); (121a)
which respectively correspond to the eigenvectors
(121b)
More than one phase in which a positive stiffness is recovered, after some hinge has completely
exhausted its load carrying capacity, can be observed. In these phases, the overall behaviour is
analogous to that of a similar beam with reduced redundancies eliminated by introducing structural
hinges in their place, allowing for completely free rotations.
7. Conclusions
This paper focused on the study of elastic media with embedded softening interfaces. The main
assumptions and results are the following.
a) A priori known interfaces, embedded in an elastic body, have been considered as loci of possible
displacement discontinuities which result from fracture and decohesion phenomena.
b) Softening imerface constitutive laws, relating tractions to displacement discontinuities, have been
used. The considered interface constitutive laws belong to the class of standard dissipative systems.
Some examples have been given and discussed.
c) The spatially discretized model of the system has been derived from a weak formulation of the
governing equations. Spatial discretization has been carried out by modelling all interface fields in a
mixed finite element context, based on the notion of generalized variables, while a standard
displacement-based formulation has been adopted for the elastic body.
d) Compatibility has also been imposed in a weak form along the interfaces, thus allowing for
completely free discretizations of the parts of the body separated by the interfaces, and along the
interfaces themselves.
e) A systematic procedure, apt to reduce the unknowns of the discrete problem to the interface
variables only, has been proposed. The reduced problem can be solved with a limited computational
effort.
f) A discussion on uniqueness of the rate problem and on stability of the equilibrium states has been
presented. This has been done by following both the classical Hill’s approach, applied here to the
problem discretiz,ed in space, and the so-called intrinsic approach which allows for the explicit
derivation of sufficient conditions for uniqueness and stability.
g) Some examples have been presented, which help in clarifying the formulation and the related
results.
The proposed #discrete formulation can be used to model different phenomena, ranging from the
study of frames with softening hinges to the simulation of delamination in composite structures. The
introduced reduction to interface variables only is particularly advantageous, since uniqueness and
stability must be checked in a step-by-step numerical analysis.
The paper focused on general aspects of the formulation. A study of the algorithmic aspects and
numerical examples which show the potentiality of the proposed methodology are in progress, and will
be presented in a forthcoming paper.
Various related topics deserve further study:
l the choice of the interpolations for the interface variables;
l the choice of the numerical integration scheme on the interfaces;
l the development of criteria to choose among, and to follow, various equilibrium paths in the presence
of multiplicity of solutions, and the related algorithms;
l the case of moving interfaces;
l the extension of the present formulation to the case of interfaces embedded inside the finite elements.
Acknowledgments
This work has been carried out in the framework of a project financed by the Italian Ministry for
University and Scientific Research (MURST 40%).
References
[l] A. Hillerborg, M. Modeer and P.E. Petersson, Analysis of crack formation and crack growth in concrete by means of
fracture mechanics and finite elements, Cem. Concr. Research 6 (1976) 773-782.
[2] L. Cedolin, S. dei Poli and I. Iori, Tensile behaviour of concrete, ASCE J. Engrg. Mech. 113 (1987) 431-439.
[3] A. Carpinteri, Softening and snap-back instability of cohesive solids, Int. J. Numer. Methods Engrg. 29 (1989) 1521-1538.
[4] Z. Cen and G. Maier, Bifurcations and instabilities in fracture of cohesive-softening structures: a boundary element analysis,
Fatigue Fract. Engrg. Mat. Struct. 15 (1992) 911-928.
[5] Y.N. Li and R.Y. Liang, Stability theory of cohesive crack model, J. Engrg. Mech. ASCE 118 3 (1992) 587-603.
[6] W.H. Gerstle and Xie Ming, FEM modeling of fictitious crack propagation in concrete, ASCE J. Engrg. Mech. 118 (1992)
416-434.
[7] A.R. Ingraffea, Theory of crack initiation and propagation in rock, in: B. Atkinson, ed., Rock Fracture Mechanics, Chap. 3
(Academic Press, New York, 1987).
[8] E.E. Alonso, I. Carol and A. Gens, An interface element formulation for the analysis of soil-reinforcement interaction,
Comput. Geotechnics 7 (1988) 133-151.
[9] A.W. Page, Finite element model for masonry, ASCE J. Struct. Engrg. 104 (1978) 1267-1285.
[lo] J.G. Rots, Numerical simulation of cracking in structural masonry, Heron 36 (1991) 49-63.
[ll] H.R. Lotfi and P.B. Shing, Interface model applied to fracture of masonry structures, ASCE J. Struct. Engrg.120 (1994)
63-80.
[12] A. Needleman, A continuum model for void nucleation by inclusion debonding, J. Appl. Mech. 54 (1987) 525-531.
[13] A.J. Levy, The debonding of elastic inclusions and inhomogeneities, J. Mech. Phys. Solids 39 (1991) 477-505.
[14] J. Backlund, Fracture analysis of notched composites, Comput. Struct. 13 (1981) 145-154.
[15] P. Ladeveze, A damage computational method for composite structures, Comput. Struct. 46 (1992) 79-87.
[16] A. Corigliano, Formulation, identification and use of interface models in the numerical analysis of composite delamination,
Int. J. Solids Struct. 30 (1993) 2779-2811.
[17] 0. Allix, P. Lade&e and A. Corigliano, Damage analysis of interlaminar fracture specimens, Composites Struct. 31 (1995)
61-74.
[18] J.C.J. Schellekens and R. de Borst, Application of anisotropic softening plasticity to mixed mode delamination in
composites, Proc. COMPLAS III, Barcelona (1992) 1671-1682.
[19] T. Stankowski, K. Runesson and S. Sture, Fracture and slip of interface in cementitious composites, ASCE J. Engrg. Mech.
119 (1993) 292-327.
[20] U. Edlund and A. Klarbring, Analysis of elastic-plastic adhesive joints using a mathematical programming approach,
Comput. Methods Appl. Mech. Engrg. 78 (1990) 19-47.
[21] U. Edlund, Surface adhesive joint description with coupled elastic-plastic damage behaviour and numerical applications,
Comput. Methods Appl. Mech. Engrg. 115 (1994) 253-276.
[22] E.N. Dvorkin, A.M. Cuitiiio and G. Gioia, Finite element with displacement interpolated embedded localization lines
insensitive to mesh size and distortions, Int. J. Numer. Methods Engrg. 30 (1990) 541-564.
[23] R. Larsson, K. Runesson and N.S. Ottosen, Discontinuous displacement approximation for capturing plastic localization,
Int. J. Numer. Methods Engrg. 36 (1993) 2087-2105.
[24] J.C. Simo, J. Oliver and F. Armero, An analysis of strong discontinuities induced by strain-softening in rate-independent
inelastic solids, Comput. Mech. 12 (1993) 277-296.
[25] G. Maier, On structural instability due to strain softening, in Instability of Continuous Systems, Proc. IUTAM Symposium
Herrenalb 1969, Leipholz editor (1971) 411-417.
[26] G. Maier, A. Zavelani Rossi and J.C. Dotreppe, Equilibrium branching due to structural softening, ASCE J. Struct. Div. 99
(1973) 897-906.
[27] Z.P. Baiant, G. Pijaudier-Cabot and J. Pan, Ductility, snapback, size effect and redistribution in softening beams or frames,
ASCE J. Struct. Eng. 113 (1987) 2348-2364.
[28] A. Cipollina, A. Lopez-Inojosa and J. Florez-Lopez, A simplified damage mechanics approach to nonlinear analysis of
frames, Comput. Struct. 54 (1995) 1113-1126.
[29] W. Prager, The general theory of limit design, Proc. of the 8th Int. Conf. on Appl. Mech., Instanbul 2 (1952) 65-72.
[30] G. Maier, Quadratic programming and theory of elastic-perfectly plastic structures, Meccanica 3 (1968) l-9.
[31] L. Corradi, On compatible finite element models for elastic plastic analysis, Meccanica 13 (1978) 133-150.
[32] L. Corradi, On stress computation in displacement finite element models, Comput. Methods Appl. Mech. Engrg. 54 (1986)
325-339.
[33] A. Capsoni and L. Corradi, A mixed finite element model for plain-strain elastic-plastic analysis: part I-formulation and
assessment of the overall behaviour, Comput. Methods Appl. Mech. Engrg. (to appear).
[34] C. Comi, G. Maier and U. Perego, Generalized variable finite element modelling and extremum theorems in stepwise
holonomic elastoplasticity with internal variables, Comput. Methods Appl. Mech. Engrg. 96 (1992) 213-237.
[35] C. Comi and U. Perego, A unified approach for variationally consistent finite elements in elastoplasticity, Comput. Methods
Appl. Mech. Engrg. 121 (1995) 323-344.
[36] C. Comi and U. Perego, A regularization technique for elastoplastic softening analyses based on generalized variables, Proc.
COMPLAS IV, Barcelona (1995) 535-546.
G. Bolzon, A. Corigliano I Comput. Methoak Appl. Mech. Engrg. 140 (1997) 329-359 359
[37] C. Comi and A. Corigliano, On uniqueness of the dynamic finite-step problem in gradient-dependent softening plasticity,
Int. J. Solids Struct. 33 (1996) 3881-3902.
[38] R. Hill, A general theory of uniqueness and stability in elastic-plastic solids, J. Mech. Phys. Solids 6 (1958) 236-249.
[39] G. Maier and D.C. Drucker, Effects of geometry change on essential features of inelastic behaviour, ASCE J. Mech. Div. 99
(1973) 819-834.
[40] Q.S. Nguyen, Bifurcation et stabilitt des systtmes irrtversibles obtissant au principe de dissipation maximal, J. de
Micanique Thiorique et Appliquee 3 (1984) 41-61.
[41] Q.S. Nguyen, Bifurcation and post-bifurcation analysis in plasticity and brittle fracture, J. Mech. Phys. Solids 35 (1987)
303-324.
[42] Q.S. Nguyen, Bifurcation and stability of dissipative systems, CISM n 327, Udine (Springer-Verlag, 1993).
[43] Q.S. Nguyen, Some remarks on plastic bifurcation, Eur. J. Mech. A/Solids 9 (1994) 485-500.
[44] Q.S. Nguyen, C. Stolz and G. Debruyne, Energy methods in fracture mechanics: stability, bifurcation and second variations,
Eur. J. Mech. A/Solids 9 (1990) 157-173.
[45] X.P. Xu and A. Needleman, Numerical simulation of fast crack growth in brittle solids, J. Mech. Phys. Solids 42 (1994)
1397-1434.
[46] A.M. Alvaredo and R.J. Torrent, The effect of the shape of the strain-softening diagram on the bearing capacity of concrete
beams, Mat. Struct. 20 (1987) 448-454.
[47] G.V. Guinea, Modelling the fracture of concrete: the cohesive crack, Mat. Struct. 28 (1995) 187-194.
[48] J.C.J. Schellekens and R. de Borst, On the numerical integration of interface elements, Int. J. Numer. Methods Engrg. 36
(1993) 43-66.
[49] Z.P. Baiant and L. Cedolin, Stability of Structures (Oxford Univ. Press, New York, 1991).
[50] G. Maier, On softening flexural behaviour in elastic-plastic beams, Rend. 1st. Lombard0 Sci. Lett. 102 (1968) 648-677, (in
Italian).
[51] G. Bolzon, G. Maier and G. Novati, Some aspects of quasi-brittle fracture analysis as a linear complementarity problem, in:
Z.P. Bazant, Z. Bittnar, M. Jirasek and J. Mazars, eds., Fracture and Damage in Quasi-brittle Structures (E&IN Spon,
London, 1994) 159-174.

329-359

Uploaded by

Copyright:

Available Formats

329-359

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

329-359

Uploaded by

Copyright:

Available Formats

Computer methods

A discrete formulation for elastic solids with damaging

Received 19 September 1995; revised 2 May 1996

‘Dedicated to Prof. Giulio Maier on the occasion of his 65”’ anniversary.

0045-7825/97/$17.00 10 1997 Elsevier Science S.A. All rights reserved

As a suitable schematization of the situations described in Section 1, an elastic body 0 is considered

l strain and displacements are small;

Fig. 1. Elastic body with damaging interfaces.

2.1. Governing relations

E=CU inn’ (1)

u=U onfU&, (3)

CTa+b=O in0 (4)

t=t-= -t+=e 0nr (5)

l constitutive equation for the interfaces r

t= wW, 4’. a= _ w(wey4

A 20; f(4 4 s 0; f’i = 0. (11)

In the above equations:

GtTip+SaTiy~O;V&Jia EE,, E,,={&=t* -t,&z=u* -ulf(t*,u*)sO}. (12)

2.2. Interface behuviour

rP = {x E Tl~!(x, 7) = t, for some 7 C T,}, (13b)

$=tj/(a)=f,,a+;hc?; f=t+a onr,, (14b,c)

*=O enc., (14d)

Fig. 3. Unilateral contact: a) interface law; b) corresponding energy function.

Fig. 6. The elastic domain for the mixed-mode cohesive model.

r,= xEr]cX(x)> -+ . (21)

Fig. 7. The elastic domain for the elastic-plastic interface model.

(w-u+ +u-)T&ddT=O V6t, (24)

SW”dT = 0, V 6~“; 8adT=O,V8a, (27)

Ir(*pT&+iyT&)d~c0,V&,6aEE,, E~“{6t=t*-t,sa=a*-a~f(t*,a*)~o}, (28)

where St=&(x), &z=&(x), xfZr.

defined for i~kl={&IisO on r}.

4.1. Generalized variable formulation

T = r M;(x)t(x) dT; W = r M;(x)w(x) dT. (35)

t(x) = M,(x)T; w(x) = M,(x)W on r, (36bc)

a(x) = M, @)A ; a(x) = M,(x)cY on f. (36d,e)

6TT M~(x)(i14w(x)W-Nu(x)~r+U + N,,(x)l,-U) dT = 0, V 6T.

w= [i, ikf~(x)&(x)l, d++ - C-)U = [I,- ~~(~)~Jx)~, d+Jr+ - u’-). (40)

4.4. Reduction to interface degrees of freedom

c&J p = 0; KVP+ C;T'=p. (440)

crvw=w; KVW + C;TW = 0. (450)

TP = K&K-‘P; TW= -K,W, (Mb,c)

v~ = {C@ VrTj; UrT = {p+T Ur-Tj,

w= Q.vr; RVr+d;T=B, (490)

er E [Zcnrxnr) -z~+~n,jl; f = K,, - KrO,K,f,sK,.,; p = Pr - K,,.K,:o.Po’. (49cAe)

TP = ;(er$ - tfr&J”); TW= -K,W, (52b,c)

4.5. Interface constitutive law

8WeTlr M:($)(M,(x)T - -$(we(x), a(x))) dr = 0, V &We, (53)

SCYT/r M:(x) (&f&$4 + -$w~(x), a(x))) dT = 0, V 6a!, (54)

6TT M:(x) M,(x)~#~ -g@(x), a(x))M,(x)~ dr = 0, V 6T, (55)

&AT MT(x) M,(x)& - g(<x), a(x))M,(x)A dF = 0, V &A, (56)

sAT M;f(x)f(t(X), a(x)) dT G 0 V S/i EM, M={6A=A*-/iIA*>o}, (57a,b)

ty(We, (u) =s1, JI(M,(x)W’, M,(x)(Y) dF, (58)

T = aF(we, a). A = _ fww”7 a).

iYK,ti’==‘; P’ = K,C,K-‘P, (64)

5.1. Uniqueness and stability by the classical approach

AW’A? + Alir’K,AI$‘= 0, (65)

Aw=A? + Aw=K,AI+ > 0, V Arir = w, - I+$ # 0. (66)

SUP = K-l@- C;K,-C,K-‘) SP, (68b)

6U w = K -‘C;K, SW, (68c)

2ij2Z = GUpTK 6Up= + 2 GUprK MJw + SWTKr 6W + SW= 6T. (69)

GtTip+SaTiy~O;V&Jia EE,, E,,={&=t* -t,&z=u* -ulf(t,u)sO}. (12)

Ir(pT&+iyT&)d~c0,V&,6aEE,, E~“{6t=t-t,sa=a-a~f(t,a*)~o}, (28)

sAT M;f(x)f(t(X), a(x)) dT G 0 V S/i EM, M={6A=A-/iIA>o}, (57a,b)

(T-T)TI+p+(A-A)TiY~O; V(T,A)EE, E = {T, AjF(T, A) s O}. (81)

(F-F)'(ci--st')'~o, vi- ~zv,(x). Plb)

p1=$(3-h); &=%(9-h); (121a)