A damage model for a finite thickness composite interface accounting for in-plane deformation

Engineering Fracture Mechanics 163 (2016) 396–415

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

A damage model for a finite thickness composite interface accounting for in-plane deformation Francesco Freddi a, Elio Sacco b,⇑ a b

Department of Civil Environmental Engineering and Architecture, University of Parma, Parco Area delle Scienze 181/A, I 43124 Parma, Italy Department of Civil and Mechanical Engineering, University of Cassino and SL, Via G. Di Biasio 43, I 03043 Cassino, Italy

a r t i c l e

i n f o

Article history: Received 27 October 2015 Received in revised form 24 May 2016 Accepted 3 June 2016 Available online 17 June 2016 Keywords: Cohesive interface Damage mechanics Enriched kinematics Detachment

a b s t r a c t The present paper deals with the formulation of an innovative interface model characterized by nonlinear response. The interface is assumed to be composed of several thin layers, one of which is characterized by nonlinear response. In this layer, denoted throughout the paper as detachment layer, the nonlinearity is modeled adopting a damage law and the decohesion process occurs. The kinematics of the interface is enriched accounting for the in-plane strain component; in such a way, the presence of the stress component arising in the plane of the interface can be computed and taken into account for the evaluation of the nonlinear response of the detachment layer. In fact, when in-plane compressive stresses are present a confinement effect occurs, slowing down the damage growth. On the contrary, when tensile in-plane stress is present in the detachment layer, the damage evolves more rapidly. The interface model is implemented in a finite element computer code and a numerical procedure is proposed. Two numerical applications are presented to assess the ability of the proposed interface model and of the implemented numerical procedure to reproduce specific decohesion problems. The first application deals with the debonding of a FRP reinforcement from a concrete support; comparison of the numerical results with available experimental evidences are given. The second application concerns the mixed mode decohesion in an edge-notch four-point bending specimen. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction One of the most evident damage mechanism of material is the formation of well defined regions where inelastic deformation localises and subsequently propagates. When the size of these regions is small enough (compared to the structural scale), it is computationally advantageous to idealize the damage zone as an interface, while the remaining parts of the structure can still be considered in elastic regime. Interface can be defined as the link between the two surfaces containing a small portion of material with high deformation gradient; it is characterized by zero thickness and it is able to transmit stresses from one surface to another. In this way, the full continuum constitutive law is replaced by a formulation relating the interfacial tractions to the separation between the two surfaces. Moreover, the interface stresses are assumed to be linear or nonlinear functions of relative displacements. Interface adhesion and cohesion are crucial in many technological areas ranging from the scales of kilometers in geosciences to micro- and nano-meters in biological systems, in microelectronics and nano-device applications. Moreover, ⇑ Corresponding author. E-mail addresses: [email protected] (F. Freddi), [email protected] (E. Sacco). http://dx.doi.org/10.1016/j.engfracmech.2016.06.001 0013-7944/Ó 2016 Elsevier Ltd. All rights reserved.

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

397

Nomenclature t t uþ u s E

m C K D

e

c

r ed r
d r
eq j ju

U N1 N2

interface total thickness detachment layer thickness displacement vector of the top interface surface displacement vector of the bottom interface surface relative displacement vector elastic modulus of the interface Poisson ratio of the interface elastic isotropic constitutive matrix of the interface elastic interface stiffness matrix damage variable strain in the interface control displacement vector stress in the interface strain in the detachment layer effective stress in the detachment layer Drucker–Prager equivalent effective stress limit threshold equivalent effective stress ultimate value of the equivalent effective stress element nodal displacement vector interpolation (shape) finite element functions

the continuous development of adhesive bonding technology has contributed to the diffusion of structural adhesives in civil and mechanical engineering applications. The interfaces are crucial regions governing the strength and stability of structures. In fact, local rupture and material de-cohesion develop in macroscopic fractures and separations between parts. Consequently, an accurate understanding of fracture initiation and propagation has become ever more important from the serviceability and safety standpoints for structures. Huge amount of literature is devoted to this topic. Just to cite few, studies have been conducted to describe the detachment mechanism in composite laminates [2,4,11,12,23,24,35]; to analyze the damaging process of composite (heterogeneous) materials from the micro-mechanical point of view; to model the adhesion between different type of solids, e.g. [25,30,41,42,45]; to determine the surfaces of potential crack in fracture mechanics and to describe the response of cohesive materials [1,22], even implemented in XFEM procedures [6]; to model the response of softening materials subjected to strain and damage localization [36,46] and to material instabilities due to the formation of shear bands [20,28,29]; to simulate the response of structural parts joined by adhesive [19,34,43]; to reproduce the detachment and the interaction between solids in other applications [3,8,16,21,26,27,33,38,39]. Structural collapse in composite structures is often caused by the appearance and evolution of different damage phenomena in a narrow region near the gluing surfaces. Usually the interface should reproduce the mechanical response of a material layer with a specific thickness where rupture occurs. The failure zone presents a size that is not necessarily equivalent to the thickness of the whole interface but mainly depends on the type of rupture and consequently on the involved materials (brittle, quasi-brittle or ductile). Failure may occur within a single material, an adherent as in the case of debonding or in the adhesive in the so called cohesion failure. Moreover, delamination appears interfacially between the adhesive and one the adherent [40]. At this level, the peculiar characteristic of each material is fundamental. For example, it is well known that in quasi-brittle materials the rupture is due to micro-cracks which form first from the undamaged material in small-width bands usually referred to as the process zone, and eventually coalesce to form a main crack [17]. Thus a material intrinsic length-scale, that characterizes the width of the process-zone can be defined; for example in conglomerates it equals 2
3 times the characteristic size of the constituent grains. Starting from this idea we introduce for the interface a similar intrinsic length-scale called rupture thickness that characterizes the size of the failure zone and depends directly on the type of failure and on the intrinsic material properties. Moreover, it is possible to properly define the mechanical response of the whole interface as a function of the failure region. A clear example is the case of a FRP bonded to the cohesive substrate by means of epoxidic adhesives which are characterized by high strength. The response of the FRP-substrate bond is often (and almost successfully) reproduced adopting interface models [15,44]. When applied, the adhesive results almost fluid so that it is able to penetrate the pores of the substrate for a significant thickness, improving the mechanical characteristics of a skin-deep layer of the support material. The detachment process of the FRP from the substrate occurs inside a thin layer of the cohesive support located below the skindeep layer [7], as schematically illustrated in Fig. 1.

398

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

Fig. 1. Interface of FRP-substrate.

Thin layer cohesive failure

Adhesive failure

Cohesive failure

Light-tear failure

Fig. 2. Failure modes of adhesive joints subjected to shear.

Another clear example is the failure mode of adhesive joints in shear tests. Ruptures are extremely vary and are illustrated in the guidelines of the ASTM D5573 [13]. The main rupture modes suggested by the standard that are usually modeled with cohesive interface models are shown in Fig. 2. The key factor which determines the failure mode of the specimen is the difference between the cohesive and adhesive resistances and the strength of the substrate. The failure firstly occurs at the weakest point. In such a way, the interface is composed by two or three ideal layers of which only one is characterized by nonlinear behavior. For example, for the cases denoted as thin-layer cohesive failure and cohesive failure two external layers surrounds a central zone where rupture occurs. The other two remaining cases can be idealized by two adjacent layers characterized by linear elastic behavior. In the present paper, a new interface model is proposed; its response is deduced assuming that the interface is composed by three different layers: two interface layers are characterized by a linear elastic response, while the third interface layer is responsible for the detachment process. In particular, the failure of the detachment layer is modeled considering a cohesive damage law, which accounts for the confinement effect according to the modeling approach recently proposed in [18]. In fact, an enriched kinematics is considered for cohesive interface, defined by the relative displacement occurring between the two surfaces of the interface and, even, by the strain arising in the plane of the interface. As in-plane stresses are determined according to the enriched model, the Drucker–Prager failure criterion can be adopted for reproducing the nonlinear response of the detachment interface, which takes into account the in-plane (confinement) stress component. The overall behavior of the interface system is recovered by a homogenization process accounting for the mechanical response of the linear elastic layer and of the nonlinear detachment layer. The proposed model is implemented in a finite element code able to solve nonlinear evolution problem by a predictor–corrector algorithm implemented in the backward Euler time-step procedure. Two numerical applications are presented to assess the ability of the proposed interface model and of the implemented numerical procedure to reproduce specific decohesion problems. The first application deals with the debonding of a FRP reinforcement from a concrete support; comparison of the numerical results with available experimental evidences are given. The second application concerns the mixed mode decohesion in an edge-notch four-point bending specimen. The main novelties of the present paper can be summarized in the following items:  a new interface model characterized by a finite thickness is proposed; the model is obtained assembling different layers one of which is responsible for the detachment process;  a simple homogenization procedure is proposed to derive the overall response of the composite interface;  a specific numerical procedure with the determination of the fully consistent tangent stiffness matrix for the interface element is derived;  two original numerical simulations have been presented, the first one concerning a comparison with available experimental results. The paper is organized as follows. In Section 2, the proposed interface model is presented with its main features; in Section 3 the numerical strategy is illustrated; Section 4 presents two interesting numerical simulations. Finally, Section 5 gives some conclusions of the present study.

399

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

2. Interface model In this section, the interface model is presented in the framework of the small displacement and strain regime, limiting the analysis to two-dimensional structural problems. Herein, the case of large displacement, necessary when the underformed configuration cannot any longer approximate the deformed one, is not considered. The presence of large displacements can be accounted for formulating the problem in finite strain framework or adopting a corotational approach [31,32,37]. In the following the Voigt notation is adopted, so that strain and stress are represented in vector form. Typically, interface models assume zero thickness and are based on the kinematics described by the relative displacements occurring between the two surfaces defining the thin detachment layer where the failure takes place. In such a way, the interface response is governed by the traction (normal and shear stresses) evaluated on the plane of the interface. Of course, this kinematics is approximated even if it results very effective for most of the engineering problems. In the following, it is assumed that the bond interface is composed by three layers, with the total thickness denoted as t  – 0. The detachment layer responsible for the damage and decohesion is positioned between two elastic layers and it is characterized by a thickness t. The distinction between t and t is important in the next modeling developments. Moreover, an enriched kinematics is assumed for the detachment layer, characterized not only by the normal and shear interface relative displacement components, but also by the presence of the in-plane strains. The use of a richer kinematics can be important in same interesting application, i.e. when the in-plane strains and, consequently, stresses cannot be neglected, playing a significant role in the overall response of the detachment layer. Because of the nonzero thickness of the interface and of a rich kinematics, a complete stress state can be evaluated inside the interface itself and, in particular, in the detachment layer, where nonlinear effects take place; in fact, the in-plane stress component arises, playing an important role in the nonlinear response of the interface [14,18]. The presence of compressive in-plane normal stress leads to a confinement effect, which can improve the behavior of the interface. On the contrary, tensile in-plane normal stress can induce a reduction of the interface strength. In Fig. 3 the interface geometry, the reference system ðO; xt ; xn Þ and the sketch of the stress state in the thin detachment layer are illustrated. þ T The displacement components of the top and bottom surfaces of the bond interface are denoted as uþ ¼ fuþ t un g and  T u ¼ fu t un g , respectively. The relative displacement vector, i.e. the jump of displacement arising between the two surfaces of the interface, is defined as s ¼ uþ  u resulting in components: 

sn ¼ uþn  un ;

st ¼ uþt  ut :

In order to define the kinematic state of the detachment layer, a very simple analysis is accomplished. It is assumed that the three layers composing the interface are characterized by the same elastic modulus E; the detachment layer, subjected to damage, has an instantaneous (secant) modulus equal to ð1  DÞ E, with D the damage variable. The strain in normal direction inside the two elastic layers and inside the detachment layer, eun and edn , respectively, are determined solving the compatibility and equilibrium equations:

eun ðt  tÞ þ edn t ¼ sn ð1  DÞE edn ¼ E eun

ð1Þ

Then, the strain in the detachment layer and the homogenized elastic modulus of the bond interface are:

edn ¼ q

sn ; t

E ¼ b E:

ð2Þ

where

q¼

1 ; hð1  DÞ þ D

ð3Þ

b ¼ q h ð1  DÞ;

ð4Þ

with h ¼ t  =t.

xn

t*

σn

t

σt

τ nt

τ nt

σt

xt

σn

Fig. 3. Interface geometry, reference system and sketch of the stress state.

400

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.2

0.4

0.6

0.8

1

Fig. 4. Overall adhesion parameter b vs damage D for different values of ratio h.

The same type of (simplified) analysis can be performed for the shear effect, leading to:

s t

G ¼ b G:

cdnt ¼ q t ;

ð5Þ

Note that the quantity b represents the overall adhesion parameter of the whole interface. In Fig. 4, the plot of the adhesion parameter versus the damage of the detachment layer is given for different values of the ratio h. It can be noted that when h ¼ 1, i.e. t  ¼ t, the linear law b ¼ 1  D is recovered; on the contrary, for h > 1 a nonlinear relationship between b and D is obtained. Denoting with r ¼ frt rn snt gT and e ¼ fet sn =t  st =t gT the stress and strain vectors in the interface, the overall constitutive law of the bond interface can be written in the form:

r ¼ b C e ¼ b K c;

ð6Þ

with

2 1 m E 6 m 1 C¼ 4 1  m2 0 0

0

3

0 7 5 1m 2

the elastic isotropic constitutive matrix of the interface material under the plane stress assumption, m denoting the Poisson ratio; in Eq. (6), K ¼ C=t  is the elastic interface stiffness matrix and c ¼ e t the control displacement vector, et representing the longitudinal strain in the interface which is set as:

et ¼

 x   1 þ e þ et þ n eþt  et ; 2 t t

ð7Þ

 where eþ t and et are the longitudinal strains occurring in the upper and lower surfaces of the interface. The strain state arising in the thin detachment layer is evaluated taking into account formulas (2)1, (5)1 and (7) for the computation of the normal, shear and longitudinal strains, respectively, resulting:

8 > <

9

8

9

edt > = > < et > = d e ¼ edn ¼ q stn ; > : d > ; > : st > qt; cnt

ð8Þ

Certainly, the position xn has to be specified for a proper evaluation of the longitudinal strain and, as consequence, of the stress state in the detachment layer. Particular cases can be obtained by setting xn ¼ 1=2: only the deformation of the upper/lower surface is taken into account. Contrarily, by choosing xn ¼ 0 the average longitudinal strain of the overall interface is considered.

401

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415


dt r
dn s
dnt g , i.e. the stress
d ¼ fr In the framework of the Damage Mechanics, the effective stress in the detachment layer r in the undamaged part of the material, is obtained by the relationship: T

r
d ¼ C ed ¼ Kd cd ;

ð9Þ

where Kd ¼ C=t is the detachment layer stiffness matrix and cd ¼ ed t is the detachment layer control displacement. Setting cde ¼ tfet 0 0gT and cds ¼ f0 sn st gT , Eq. (9) can be rewritten in the form:





r
d ¼ Kd cde þ q cds ;

ð10Þ

As remarked above, the detachment layer is considered to be subjected not only to the classical interface traction stresses rn and snt , but also to the in-plane normal stress rt . The presence of rt could play an important role in the damage evolution and in the failure of the detachment layer. For this reason, a pressure dependent failure criterion, that accounts for the triaxiality and confinement effect, is considered for the damage evolution of the thin detachment layer. In particular, the Drucker–Prager criterion is considered, as it is simple and generally effective for brittle and cohesive materials. Moreover, it results to be almost flexible, as a proper setting of the governing material parameters permits to recover as special cases the Mises as well as the regularized form of the Mohr-Coulomb criteria. Of course, alternative even more complex failure criteria can be adopted for describing the detachment process. In this framework, the first invariant of the effective stress, sm , and the second invariant of the deviatoric part of the effective stress, sd , are introduced as:

 1 d m T d
þr
dn
¼ r sm ¼ ðj Þ r 3 t  0:5 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  1 d T d d 1 2
Þ J r
ðr ¼ r^ þ r^ 2n þ r^ 2b þ 2ðs^2nt Þ sd ¼ 2 2 t

; ð11Þ ;

where it is set:

8 9 > <1> = 1 m 1 ; j ¼ 3> : > ; 0

2 3 2 1 0 1 6 7 Jd ¼ 4 1 2 0 5; 3 0 0 6

ð12Þ

In Eq. (11), the components of the deviatoric effective stress result:

 1 d
t  r
dn 2r 3  1 r^ n ¼ 2 r
dn  r
dt 3  1 r^ b ¼ r
dt  r
dn 3 s^nt ¼ s
dnt ;

r^ t ¼

ð13Þ

with sb the non-zero component evaluated in the direction orthogonal to the xt  xn plane. The equivalent effective stress is introduced as:

r
eq ¼ sd þ a sm ;

ð14Þ

where a is the material parameter ruling the influence of the hydrostatic pressure ðsm Þ on the failure. According to the Drucker–Prager criterion, failure occurs when it occurs:


eq  j ¼ 0; f ¼r

ð15Þ

with j the limit threshold equivalent effective stress, which allows the starting of the damage evolution.
t  r
n , setting In Fig. 5, the elastic domain defined by the Drucker–Prager criterion is illustrated in the plane r j ¼ 7:7 MPa: In particular, in Fig. 5(a) the elastic domains obtained for a ¼ 1; 1:65; 2:5 with s
nt ¼ 0 are reported; the considered values for j and a correspond to materials characterized by uni-axial tensile and compressive failure stresses equal to rc0 ¼ 31:56 MPa and rt0 ¼ 8:46 MPa, rc0 ¼ 281:53 MPa and rt0 ¼ 6:83 MPa, rc0 ¼ 1 and rt0 ¼ 5:46 MPa, in the last
nt ¼ 0; 5; 10 with a ¼ 1:65 are case no compressive uni-axial failure can occur. In Fig. 5(b) the elastic domains obtained for s
n  s
nt , setting j ¼ 7:7 MPa. Fig. 6(a) illustrates the given. Analogously, in Fig. 6 the elastic domain is plotted in the plane r
nt ¼ 0, while Fig. 6(b) the domains for r
n ¼ 0; 5; 10 and a ¼ 1:65. domains for a ¼ 1; 1:65; 2:5 and s A damage evolution law, characterized by a linear softening branch in the stress–strain relation, is considered; the damage variable is given by the relationship:

e D ¼ maxð0; minð1; DÞÞ

with D_ P 0;

ð16Þ

402

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415


t  r
n obtained (a) for a ¼ 1; 1:65; 2:5 with s
nt ¼ 0 or (b) for s
nt ¼ 0; 5; 10 and Fig. 5. Elastic domain defined by the Drucker–Prager in the plane r a ¼ 1:65.

with


eq Þ e ¼ ðj  r D ; r
eq ðj  ju Þ

ð17Þ

being ju the ultimate value of the equivalent effective stress. In Fig. 7 the shear stress snt versus the slip st is plotted for different values of the ratio h; computations are performed setting E ¼ 28; 700 MPa; m ¼ 0:2; j ¼ 7:7 MPa; ju ¼ 90 MPa and t ¼ 1 mm. It can be noted that for any value of h a linear softening branch occurs after the threshold relative displacement; moreover, increasing h the softening branch becomes more and more steep until leading to a snap-back phase.

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

403


n  s
nt obtained (a) for a ¼ 1; 1:65; 2:5 and r
t ¼ 0 or (b) for r
t ¼ 0; 5; 10 and a ¼ 1:65. Fig. 6. Elastic domain defined by the Drucker–Prager in the plane r

3. Numerical procedure The proposed damage interface model above presented has been implemented into a four-node interface element [18], which has been introduced as user element in an ad-hoc version of the Open Source library deal.II [5]. In particular, the displacement fields uþ and u are interpolated using linear shape functions. In such a way, the relative displacement s is approximated by a linear function, while the longitudinal strain et results to be constant along the interface. The time integration of the nonlinear evolutive constitutive relationships is performed in the framework of a quasi-static incremental/iterative solution procedure adopting a backward-Euler implicit algorithm. The time interval of interest is subdivided in subincrements and the evolutive problem is solved into a typical interval ½t n ; t nþ1 , being t nþ1 > t n . For brevity of notation, the index n denotes the quantities evaluated at the time t n , while subscript is omitted for all quantities evaluated at

404

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

Fig. 7. Shear stress vs slip for different values of the ratio h.

Fig. 8. Interface finite element: coordinate system and nodes.

the time t nþ1: The time-discrete constitutive model for the interface is solved at the Gauss quadrature points of the FE discretization. T

With reference to Fig. 8, let U ¼ fu1t u1n u2t u2n u3t u3n u4t u4n g denote the element nodal displacement vector, where the superscript indicates the node number; the control displacement vector c, is obtained as:

2 t 2

c ¼ BU

N01

t 2

0

6 with B ¼ 4 0

N1

N02

t 2

0

N01

0

t 2

N02

0

3

N1

0

N2

0

N1

0

7 N2 5;

0

N2

0

N1

0

N2

0

ð18Þ

and

N1 ¼

1 ð1  nÞ 2

N2 ¼

1 ð1 þ nÞ 2

N01 ¼ N02 ¼ 

1 2j

ð19Þ

being j ¼ x0t ¼ L=2 with L the length of the interface element; the prime apex denotes the derivative with respect to the abscissa n defined in the parent element. Once the nodal displacements in the typical FE has been given, the detachment layer control displacement can be determined at the typical Gauss point as:



cd ¼ Bde þ q Bds U;

ð20Þ

where it is:

2 1 x  0       t 2  tn N1 0 t 12  xtn N02 0 t 12 þ xtn N 01 0 t 12 þ xtn N02 6 Bde ¼ 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 3 0 0 0 0 0 0 0 0 6 7 Bds ¼ 4 0 N 1 0 N2 0 N1 0 N2 5; N1

0

N2

0

N1

0

N2

0

0

3

7 05 0

ð21Þ

405

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

Then, the possible damage evolution has to be computed developing a predictor–corrector algorithm. Initially, a trial solution is determined, assuming the damage as frozen, i.e. setting D ¼ Dn . Using Eq. (3), the value of the parameter q is determined and, as consequence, the strain vector ed is evaluated by formula (8), allowing to get the detach
d is computed, so that the ment layer control displacement using formulas (20) and (21). Through Eq. (9) the effective stress r equivalent effective stress and the yield function are evaluated by (14) and (15), respectively. If f < 0, there is no evolution of the damage and the trial solution is the solution of the problem. On the contrary, if f P 0 the damage state evolves. In order to determine the damage evolution Eqs. (3), (4), (10), (14) and (15) have to be solved, taking into account of formulas (20) and (21). As these equations are coupled and nonlinear, they are rewritten in residual form:

Rb ¼ b ½hð1  DÞ þ D  h ð1  DÞ;

ð22Þ

Rq ¼ q ½hð1  DÞ þ D  1


d  Kd Bde þ q Bds U Rr
d ¼ r  0:5 1 d T d d m T d
eq 


Þ J r Req ¼ r ðr  a ðj Þ r 2
eq ðj  ju Þ  ju ðj  r
eq Þ RD ¼ D r

ð23Þ ð24Þ ð25Þ ð26Þ

To solve the system of Eqs. (22)–(26) a Newton–Raphson algorithm is adopted. Setting R ¼ fRb Rq Rr
d Req RD gT and
eq Dg , it results:
d r X ¼ fb q r T

 !1 @R X ¼ Xk  RðXk Þ @XXk

ð27Þ

where Xk is the solution at the k-th iteration. Explicitly, the components of the matrix @R=@X are: @Rb @b @Rb @q

¼ hð1  DÞ þ D ½1  1

@Rq @b

¼0

¼0

½1  1

@Rq @q

¼ hð1  DÞ þ D ½1  1

½1  3

@Rq T
d ¼ 0 @r @Rq
eq ¼ 0 @r

½1  3

@Rq @D

½1  1

@Rb T
d ¼ 0 @r @Rb
eq ¼ 0 @r @Rb @D

½1  1

¼ b ð1  hÞ þ h

½1  1

½1  1

½1  1

¼ q ð1  hÞ

@Rr
d @b

¼0

½3  1

@Rr
eq @b

¼0

½1  1

@Rr
d @q

¼ Kd Bds U ½3  1

@Rr
eq @q

¼0

½1  1

@Rr
d
d @r @Rr
d
eq @r

¼I

½3  3

@Rr
eq
d @r

¼0

½3  1

¼1

½1  1

@Rr
d @D

@Rr
eq
eq @r

¼0

½3  1

@Rr
eq @D

¼0

½1  1

@RD @b

¼0

½1  1

@RD @q

¼0

½1  1

T @RD
d ¼ 0 @r @RD
eq ¼ Dð @r @RD @D



T
d =sd þ a jm ¼  12 Jd r

½1  3

½1  3

j  ju Þ

½1  1


eq ðj  ju Þ þ ju ¼r

½1  1

where I is the identity matrix, 0 ¼ f0 0 0gT and the square brackets indicate the dimension of the arrays. Finally, it results:

2 6 6 @R 6 6 ¼6 @X 6 6 4

hð1  DÞ þ D

0

0T

0

0

hð1  DÞ þ D

0T

0

0

Kd Bds U

I

0 m T

0

0


d =sd þ a j Þ ðJ r

0

0

0T

d

b ð1  hÞ þ h

3

7 q ð1  hÞ 7 7 0

1

0

Dðj  ju Þ þ ju

r
eq ðj  ju Þ

7 7 7 7 5

ð28Þ

406

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

The typical finite element, schematically depicted in Fig. 8, is subjected to tractions pþ and p applied on the upper and lower surfaces of the interface due to the interaction with the solid adjacent elements. The variational form of the equilibrium equation becomes:

Z

1

1

dcT r j dn ¼

Z

1

h

i T ðduþ Þ pþ þ ðdu ÞT p j dn;

ð29Þ

1

Taking into account Eqs. (6) and (18) and the condition pþ ¼ p ¼ p, the equilibrium equation of the finite element, written in residual form, results:

Z

1

r¼ 1

where

Nþ ¼



Z b BT K B U j dn þ

1

T

ðNþ  N Þ p j dn ¼ 0;

ð30Þ

1

0 0 0 0 N1 0 0 0 0 0

0 N1

N2 0

0 N2

N ¼



N1 0

0 N1

N2 0

0 N2

0 0 0 0 : 0 0 0 0

ð31Þ

The element consistent tangent stiffness matrix, in the case of damage evolution, is computed by differentiating the residual defined by Eq. (30) with respect to the element nodal displacement vector, resulting:

Kt ¼

@r ¼ @U

Z

1

1

Z b BT K B j dn þ

1

 BT K B U

1

@b @U

T j dn:

ð32Þ

In order to compute the derivative of b with respect to the nodal displacement vector U, the differential of the residual R is determined, accounting also for the variation of the nodal displacement vector:

dR ¼

@R @R dX þ dU ¼ 0; @X @U

ð33Þ

which allows to determine:

 dX ¼ 

@R @X

1

@R dU; @U

ð34Þ

with the components of the derivative of R with respect to U being:

@Rb ½1  8 ¼ 0T @U @Rq ½1  8 ¼ 0T @U

@Rr
d ¼ Kd Bde þ q Bds @U @Req ½1  8 ¼ 0T @U @RD ½1  8 ¼ 0T @U

½3  8

where the dimensions of the arrays are explicitly reported between square brackets. Taking into account the definition of the vector X it results b ¼ X 1 ; then, Eq. (34) gives:

"  #  1 @b @R @R ¼ @U i @X @U



ð35Þ

1i

Formula (35) allows to compute the final term of Eq. (32) and, hence, the consistent element matrix Kt . 4. Numerical simulations and results In the following some numerical simulations will be presented in order to illustrate the features and the capabilities of the proposed interface model which is characterized by failure localized at different position in the interface thickness. Moreover, the numerical applications will stress the influence of the confinement effect on the interface response and, as consequence, on the overall response of the structural systems. In particular, two applications will be illustrated in the following subsections, dealing with:  debonding tests of FRP reinforcements from concrete supports, for which the comparison with available experimental data is given. In this simulation it has been set xn ¼ t =2 in Eq. (7), being the rupture within the substrate material and localized at the bottom of the interface

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

407

 detachment of a sandwich beam in four point bending test. Failure is assumed to occur at different positions: in the middle and at the two extremities of the adhesive layer so that it is set xn ¼ 0; t  =2 in Eq. (7). The problems were analyzed in two dimensions in plane stress condition. In the numerical examples the Cartesian coordinate system ðx1 ; x2 Þ is introduced so that the following stress components are adopted r ¼ fr11 r22 r12 gT ;
22 r
12 gT . Finally, for the discretization of the domains classical quadrilateral four node element with linear shape
11 r r
d ¼ fr functions has been adopted whereas for the interface the finite element previously presented has been used. At the interface level the size of the element has been chosen small enough to guarantee accurate solution; in both the examples a length of 0:5 mm has been assumed for the interface elements. 4.1. FRP debonding The debonding of FRP laminate glued on concrete element has been studied by numerical test according to the single lap joint shear test proposed in [9]. The experimental setup is depicted in Fig. 9. Left and bottom sides of the specimen are constrained in order to avoid displacements in the normal direction to the surface and to permit displacements tangent to it. Plate bond is withdrawn of 50 mm from the front side of the specimen to limit the boundary effect that strongly influences the rupture mechanisms [7]. The anchorage length is equal to 150 mm. The thickness of the FRP plate tf , the mechanical properties of the specimen such as the Young modulus of concrete Ec and reinforcement Ep and the interface parameters are listed in Table 1. Finally, the adopted values of the whole interface thickness t are mean values suggested in CNR DT200 [10] and has been assumed equal to 10 mm and whereas the thickness of the failure zone t has been taken equal to 3 mm. In this case the rupture thickness t  is directly correlated with the internal length of the concrete. The failure domains at the interface level are the ones previously presented in Figs. 5–7. Other details of the experimental test can be found in [9]. In this application failure occurs at the bottom of the interface within the substrate. So, in the numerical simulation only the elongation of the lower surface has been taken into account by setting xn ¼ t  =2. Load displacement curve obtained with the proposed model is reported in Fig. 10. The result shows good agreement with two experimental results given in [9]. Additionally, the numerical solution obtained by adopting for the interface a pure Mode II bilinear law with equivalent fracture energy value is added in Fig. 10. The pure Mode II interface law assumes a shear stress-slip relationship equivalent to the one suggested in [10]. The two numerical solutions present similar ascending and descending branches while the debonding with the present formulation gives higher values of transferred load. The proposed model is in agreement with the experimental results while the pure Mode II interface law gives lower value of the transmitted force. The investigation of the stress state at the interface level permits a clear explanation of the failure mechanism. In Fig. 11 the triplet of stress components r11 ; r22 ; r12 is plotted along the interface for different points of the equilibrium path: (a) 0:7 Pmax in the ascending branch, (b) Pmax and (c) umax , with Pmax and umax the maximum transmitted force and displacement respectively.

Fig. 9. (a) Test setup adopted in [9] and reproduced in the numerical simulations. (b) Elastic domains adopted for the interface.

408

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

Table 1 Geometrical and mechanical properties of specimens considered in numerical simulations. tf (mm)

Ep (MPa)

Ec (MPa)

t (mm)

t (mm)

j (MPa)

ju (MPa)

a

1.3

168,500

28,700

10

3

7.7

90

1.65

20

16

F (kN)

12

8

proposed model Mode II Experimental

4

0 0

0.1

0.2

0.3

0.4

u1 (mm) Fig. 10. Detachment curves per unit width.

At 0:7 Pmax the debonding process has already affected 15% of the interface with a stress profile for r12 that is classical of the detachment process. The normal stress r22 reveals extremely small values; only the right extremity of the glued surface shows a compression values. The axial stress r11 has non-negligible values almost the entire length so each point of the interface presents a different elastic domain. The plotted confinement stress is always in tension being the FRP plate always in elongation. The value of the axial stress reduces as debonding approached the left extremity of the reinforcement as clearly outlined by Fig. 11c.
11 are plotted Now, the constitutive relationship at the interface level is explored; stress components r11 ; r12 ; r22 and r as a function of the slip s1 at specific points of the interface. If x1 ¼ 0 is the left unloaded extremities of the reinforcement, points at coordinates x1 ¼ 5; 40; 75; 110; 140 mm will be considered. The Fig. 12 evidenced highly different interface law: these differences are due to variation of the axial stress in the debonding process. For comparison, the pure Mode II interface law is reported in same figures. Point at x1 ¼ 5 mm always
11 in the whole presents tensile axial stress whereas the right extremity at x1 ¼ 140 mm reveals compression value of r debonding process. The other points are characterized by confinement stress that moves from tension to compression. These results are in accordance with the experimental evidences of [9]. Finally, by integrating the previous curve r12  s1 it is possible to determine the value of equivalent fracture energy of pure Mode II interface law. The value of this fracture energy has been calculated and it is plotted along the anchorage in Fig. 13. In the same figure the constant value of pure Mode II is reported. At the beginning of the reinforcement an important decrement of the fracture energy occurs due to presence of peeling stresses combined with traction value of the axial stress. The other portion of interface presents increased value of the fracture energy with respect to the pure Mode II due to the
11 that has non negligible effect on the failure mechanism thus leading to an increment confinement stress component r of the transmissible force. 4.2. Edge-notch four-point bending specimens Adhesion of deposited thin films has become a key integration issue in semiconductor industry. The four-point bending test is emerging as the preferred method for the study of adhesion of thin films in semiconductor devices. In four-point bending geometry, the ‘‘sandwich” specimen is compressed by a set of four offset pins, producing a constant bending moment in between the inner pins, as schematically illustrated in Fig. 14a. A notch is cut through the lower beam to provide a point for crack initiation. Under the influence of the bending forces, the crack propagates outward along the interface towards the loading pins on either side. An advantage of this technique is that the energy release rate does not vary with crack length,

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

409

12 σ11 σ12

σ (MPa)

8

σ22

4

0 0

50

100

150

100

150

x1 (mm) -4 10

8 σ11 σ12

σ (MPa)

6

σ22

4

2

0 0

50

x1 (mm)

-2 8

σ11 σ12

6

σ (MPa)

σ22

4

2

0 0

50

100

150

x1 (mm) -2 Fig. 11. Stress distribution along the interface for (a) 0:7 P max of the ascending branch, (b) P max and (c) umax .

simplifying data acquisition and test analysis. The fracture energy of the interface can be calculated from the critical load required to propagate an interface crack in a steady-state manner. The proposed model seems to be an effective tool to model this kind of failure process because:  the thin adhesive layer can be simply identified by the interface,  rupture occurs in a thin zone with limited thickness (different from the thickness of the interface) in the adhesive layer.

410

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

8

10

8 σ11

6

proposed model pure mode II

proposed model pure mode II

4

σ 12 (MPa)

σ12 (MPa)

6

4

2

2 0 0

σ11

0.04

0.08

-2

0 0

0.04

0.08

0.12

0.12

0.16

s1 (mm)

0.16

s1 (mm)

-4 10

10 8

8

proposed model pure mode II

6

proposed model pure mode II

6

4

σ 12 (MPa)

σ 12 (MPa)

4 2 0 0

0.04

0.12

0.16

0

s1 (mm)

-2 -4

0.08

2

0

0.04

0.08

0.16

0.2

σ11

-4

-6

0.12

s1 (mm)

-2

σ11

-6

-8 10

8 proposed model pure mode II

σ 12 (MPa)

6

4

2

0 0

0.1

-4

0.2

0.3

s1 (mm)

-2

σ11

Fig. 12. Local interface stresses vs slip s1 located at different positions in interface (a) x ¼ 5 mm, (b) x ¼ 40 mm, (c) x ¼ 75 mm, (d) x ¼ 110 mm and (e) x ¼ 140 mm.

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

411

0.7

Fracture energy (N/mm)

0.65

0.6

0.55

proposed model pure mode II

0.5

0.45 0

25

50

75

100

125

150

x1 (mm) Fig. 13. Fracture energy of pure Mode II calculated along the interface.

For the case at hand, two aluminum rectangular bar glued together along the horizontal line have been considered. Geometric and mechanical properties of the material and interface are reported in Tables 2 and 3, respectively. The detachment layer has been assumed to be positioned at a middle, at the top and at the bottom of the interface, i.e. it has been set xn ¼ 0; t  =2. Moreover, the thickness of the detachment layer is taken one half of the interface thickness, i.e. h ¼ 4. Computations are performed adopting for the four points bending test a two-dimensional plane stress finite element model characterized by a quite fine mesh: the entire interface is discretized with 100 finite elements. Three types of numerical experiments are developed considering:  the complete interface proposed model with m ¼ 0:2 and xn ¼ 0; t =2;  the complete interface proposed model with m ¼ 0 and xn ¼ 0, i.e. the Poisson effect is disregarded;  the simplified interface model obtained neglecting the axial strain and xn ¼ 0, i.e. the confinement effect is ignored.

Fig. 14. Four point bending test (a) Test setup adopted for the numerical simulations. (b) Load vs vertical displacement for different interface parameters for xn ¼ 0.

412

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

Table 2 Geometrical and mechanical properties of the two specimens considered in the four point bending test. h (mm)

l (mm)

s (mm)

E (MPa)

m

10

5

170

70,000

0.23

Table 3 Geometrical and mechanical properties of the interface considered in the four point bending test. t (mm)

t (mm)

a (mm)

j (MPa)

ju

E (MPa)

m

1

0.5

1.65

5

50

1000

0, 0.2

In Fig. 14b, the value of the applied forces F versus the maximum deflection of the beam ðu2 Þ, evaluated at the middle of its length, is plotted for the three types of computations setting in all cases xn ¼ 0. This assumption permits to reproduce the so called cohesive failure. Initially, the system behaves linearly. As the damage begins at the interface level the system softens and failure along the detachment layer evolves at a constant load. The load plateau is associated with a stable propagation of the rupture along the interface. Once that detachment has reached the inner pins the system regains a linear response with a reduced stiffness. The load of the plateau varies for the three simulations: the Poisson effect reduces the value of the load whereas the classical interface model without axial effect permits higher values of transmissible force. The differences of the simulation during the damaging process are investigated in Figs. 15 and 16, where the stress along half of the interface (from the left extremity to the center) for the cases with or without the Poisson effect. The components of stress r11 ; r22 ; r12 are plotted for four values of the vertical displacement of the beam 0.5, 1, 1,5, 2.5 mm. It is clear how delamination propagates as the displacement increases. The Poisson effect reduces the transmissible shear stress, thus leading to a reduction of the transmitted load. The axial stress r11 reaches non negligible values so that failure occurs at smaller load level. Moreover, the damage distribution along the interface is plotted in Fig. 17 for different values of the prescribed displacement u2 ¼ 0:5; 1; 1:5; 2:5 mm for the two cases at hand. It clearly emerges that the Poisson effect induces faster crack propagation. Finally, the influence of the position of the detachment layer has been investigated in the case of adhesive detachment. The solution has been obtained for xn ¼ t  =2; failure layer at the top and at the bottom of the interface. In Fig. 18a the value of the applied forces F versus the maximum deflection of the beam ðu2 Þ, evaluated at the middle of its length, is plotted for the two positions of the detachment layers. The value of the force is slightly higher for the case xn ¼ t =2. This fact can be
11 along the interface. understood by investigating the distribution of the effective stress component r

4

4

σ 11

σ 11

σ 12

σ 12 σ 22

0 0

50

σ 22

2

100

σ (MPa)

σ (MPa)

2

0 0

x (mm)

50

100

x (mm)

-2

-2

-4

-4

Fig. 15. Four point bending test for xn ¼ 0: Stress distribution along the interface for (a) u2 ¼ 0:5 mm and (b) u2 ¼ 1:0 mm. Red lines refer to m ¼ 0 case, whereas black lines represent the case m ¼ 0:2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

413

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

4

4

σ11

σ11

σ12

σ12

σ22

σ22

2

0 0

50

100

σ (MPa)

σ (MPa)

2

x (mm)

0 0

50

100

x (mm) -2

-2 -4

-4

-6

Fig. 16. Four point bending test for xn ¼ 0: Stress distribution along the interface for (a) u2 ¼ 1:5 mm and (b) u2 ¼ 2:5 mm. Red lines refer to m ¼ 0 case, whereas black lines represent the case m ¼ 0:2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

1

0.8

0.6

d 2.5

0.4

1.5

1

0.5

0.2

0 0

25

50

75

100

x (mm) Fig. 17. Four point bending test for xn ¼ 0: damage distribution along the interface for u2 ¼ 0:5; 1; 1:5; 2:5 mm. Red lines refer to m ¼ 0 case, whereas black lines represent the case m ¼ 0:2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)


11 are plotted along half of the interface for two values of the vertical displacement u2 ¼ 0:5, In Fig. 18a the values of r
11 always assumes very limited values contrarily for xn ¼ t =2, due to the specific position 1.5 mm. For xn ¼ t  =2, the stress r of the neutral axis that is non perfectly in the center of the beam for the influence of boundary conditions, non negligible traction values of the axial stress induces premature failure of the detachment layer. 5. Conclusions In solid mechanics it is common practice to adopt interface to model the cohesive rupture which concentrates in narrow and well-defined region or to describe the behavior of solids connected by adhesive material. This interface model should be able to reproduce the mechanical behavior of a layer of material with a small thickness. It is not unusual that the interface incorporates the behavior of a narrow thickness that is composed by different materials with dissimilar intrinsic characteristics. At the same time the constitutive law of the interface must permit a nonlinear response. Moreover, in the classical approach information concerning the position of rupture is completely lost and the possibility to define a parameter of the internal length of the material is forbidden.

414

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

50

2.5

xn=-t*/2

40

xn=-t*/2, u2=0.5 mm

2

xn=t*/2

xn=-t*/2, u2=1.5 mm xn=t*/2, u2=1.5 mm

1.5

30

σ (MPa)

F (kN)

xn=t*/2, u2=0.5 mm

1

11

20

0.5 10

0 0

0 0

1

2

u2 (mm)

3

50

100

x (mm)

4

-0.5

Fig. 18. Four point bending test for xn ¼ t =2: (a) Load vs vertical displacement for different non linear layer locations. (b) Distribution of effective stress
11 along the interface for u2 ¼ 0:5; 1:5 mm. component r

The proposed approach overcome these limitations. In fact, the rupture thickness that has been defined represents an internal length parameter that acts to be correlated with the material in which the rupture occurs, thus, defining a process zone for the cohesive interface. The proposed interface scheme is composed by different interface layers: portions which remain in linear elastic regime and a layer responsible of the detachment process modeled considering a cohesive damage law which accounts for the axial stress. The global behavior of the interface is determined by a homogenization procedure accounting for the response of the different layers. The numerical examples clearly illustrate these peculiarities. The interface model is governed by four mechanical parameters: the two classic constants of elasticity such as the Young modulus E and Poisson coefficient m and two parameters that define the elastic limit j and the ultimate limit ju states. Again, two are the interface geometric dimensions that have to be specified: the entire thickness t  and the detachment layer thickness t. Moreover, it is necessary to define within the thickness of the interface t  the position of the rupture zone. Classical interface model can be obtained by neglecting the axial stress and by setting t  ¼ t. Finally, this innovative interface model can be easily implemented as a generalization of the classic interface element. Numerical simulations show the ability of the composite interface model in reproducing the detachment mechanism. Even if the differences of the results obtained using the interface model presented in [18] and the new composite interface model for simulation of the FRP detachment test are limited, the proposed model is preferable as it better describes the physics of the detachment phenomenon. In fact, it considers a micro-structured interface allowing for the presence of the several layers in the interface. Further numerical simulations of experimental evidences would be performed to deeply assess the performances and the abilities of the proposed composite interface model in reproducing the detachment phenomenon of micro-structured interfaces. For a proper validation ad-hoc experimental campaign should be performed. As an example, the edge-notch four-point bending test should be equipped at the interface level of a series of strain gages properly included at the interface level thus permitting to evaluate the axial elongation. Acknowledgement This research has been possible thanks to the financial support of the project ReLUIS from the Italian Ministry of the Civil Protection. References [1] Achenbach JD, Baz˘ant ZP, Khetan RP. Elastodynamic near-tip fields for a rapidly propagating interface crack. Int J Eng Sci 1976;14(9):797–809. [2] Alfano G, Crisfield MA. Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. Int J Numer Meth Eng 2001;50(7):1701–36. [3] Alfano G, Sacco E. Combining interface damage and friction in a cohesive-zone model. Int J Numer Meth Eng 2006;68(5):542–82. [4] 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. [5] Bangerth W, Heister T, Kanschat G, et al. deal.I Differential equations analysis library. Technical reference. . [6] Benvenuti E. A regularized {XFEM} framework for embedded cohesive interfaces. Comput Method Appl Mech Eng 2008;197(49–50):4367–78.

F. Freddi, E. Sacco / Engineering Fracture Mechanics 163 (2016) 396–415

415

[7] Benzarti K, Freddi F, Frémond M. A damage model to predict the durability of bonded assemblies. Part I: Debonding behavior of FRP strengthened concrete structures. Constr Build Mater 2011;25(2):547–55. [8] Bonetti E, Bonfanti G, Rossi R. Modeling via the internal energy balance and analysis of adhesive contact with friction in thermoviscoelasticity. Nonlinear Anal: Real World Appl 2015;22:473–507. [9] Carrara P, Ferretti D, Freddi F, Rosati G. Shear tests of carbon fiber plates bonded to concrete with control of snap-back. Eng Fract Mech 2011;78 (15):2663–78. [10] CNR. DT-200/2013 – Istruzioni per la progettazione, l’esecuzione ed il controllo di interventi di consolidamento statico mediante l’utilizzo di compositi fibrorinforzati, Consiglio Nazionale delle Ricerche. Rome, Italy, 2013. [in Italian]. [11] 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. [12] Cottone A, Turetta T, Giambanco G. Delamination study of through-thickness reinforced composite laminates via two-phase interface model. Compos Part A: Appl Sci Manuf 2007;38(9):1985–95. [13] Active Standard ASTM D5573. Standard Practice for Classifying Failure Modes in Fiber-Reinforced-Plastic (FRP) Joints. 2012. [14] de Borst R, Remmers JJC, Verhoosel CV. Evolving discontinuities and cohesive fracture. In: Mechanics for the world: proceedings of the 23rd international congress of theoretical and applied mechanics, {ICTAM2012}. Proc {IUTAM}, 10(0). p. 125–37. [15] Ferracuti B, Savoia M, Mazzotti C. A numerical model for FRP-concrete delamination. Compos Part B - Eng 2006;37(4–5):356–64. [16] Freddi F, Frémond M. Damage in domains and interfaces: a coupled predictive theory. J Mech Mater Struct 2006;1:1205–34. [17] Freddi F, Royer-Carfagni G. Regularized variational theories of fracture: a unified approach. J Mech Phys Solids 2010;58(8):1154–74. [18] Freddi F, Sacco E. An interface damage model accounting for in-plane effects. Int J Solids Struct 2014;51(25–26):4230–44. [19] Frémond M. Non-smooth thermomechanics. Heidelberg: Springer-Verlag; 2001. [20] Gajo A, Bigoni D, Wood DM. Multiple shear band development and related instabilities in granular materials. J Mech Phys Solids 2004;52 (12):2683–724. [21] Guiamatsia I, Nguyen GD. A thermodynamics-based cohesive model for interface debonding and friction. Int J Solids Struct 2014;51(3-4):647–59. [22] Hillerborg A, Modéer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6:773–82. [23] Hosseini-Toudeshky H, Jahanmardi M, Goodarzi MS. Progressive debonding analysis of composite blade root joint of wind turbines under fatigue loading. Compos Struct 2015;120:417–27. [24] Hosseini-Toudeshky H, Jasemzadeh A, Mohammadi B. Fatigue debonding analysis of repaired aluminium panels by composite patch using interface elements. Appl Compos Mater 2011;18(6):571–84. [25] Kushch VI. Fibrous composite with interface cracks. In: Kushch VI, editor. Micromechanics of composites. Boston: Butterworth-Heinemann; 2013. p. 331–66 [chapter 10]. [26] Lebon F. Modeling the interfaces in masonry structures. Mechanics of masonry structures, vol. 551. p. 213–40. [27] Needleman A. An analysis of tensile decohesion along an interface. J Mech Phys Solids 1990;38(3):289–324. [28] Nguyen GD, Einav I, Korsunsky AM. How to connect two scales of behaviour in constitutive modelling of geomaterials. Geotech Lett 2012;2(7– 9):129–34. [29] Nguyen GD, Korsunsky AM, Einav I. A constitutive modelling framework featuring two scales of behaviour: fundamentals and applications to quasibrittle failure. Eng Fract Mech 2014;115:221–40. [30] Oliveira DV, Lourenço PB. Implementation and validation of a constitutive model for the cyclic behaviour of interface elements. Comput Struct 2004;82 (17–19):1451–61. [31] Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J Numer Methods Eng 1999;44(9):1267–82. [32] Paggi M, Reinoso J. An anisotropic large displacement cohesive zone model for fibrillar and crazing interfaces. Int J Solids Struct 2015;69–70:106–20. [33] Palmieri V, De Lorenzis L. Multiscale modeling of concrete and of the FRP – concrete interface. Eng Fract Mech 2014;131:150–75. [34] Parrinello F, Failla B, Borino G. Cohesive-frictional interface constitutive model. Int J Solids Struct 2009;46(13):2680–92. [35] Parrinello F, Marannano G, Borino G. Mixed mode delamination analysis by a thermodynamically consistent cohesive interface model with independent mode I and mode II fracture energies. Proc Eng 2015;109:327–37 [XXIII Italian Group of Fracture Meeting, IGFXXIII]. [36] Pietruszczak S, Xu G. Brittle response of concrete as a localization problem. Int J Solids Struct 1995;32(11):1517–33. [37] Qiu Y, Crisfield MA, Alfano G. An interface element formulation for the simulation of delamination with buckling. Eng Fract Mech 2001;68 (16):1755–76. [38] Ragueneau F, Dominguez N, Ibrahimbegovic A. Thermodynamic-based interface model for cohesive brittle materials: application to bond slip in {RC} structures. Comput Method Appl Mech Eng 2006;195(52):7249–63. [39] Raous M. Interface models coupling adhesion and friction. Comptes Rendus Mécanique 2011;339(7-8):491–501. [40] Ruocci G, Argoul P, Benzarti K, Freddi F. An improved damage modelling to deal with the variability of fracture mechanisms in FRP reinforced concrete structures. Int J Adhes Adhes 2013;45:7–20. [41] Sacco E, Lebon F. A damage-friction interface model derived from micromechanical approach. Int J Solids Struct 2012;49(26):3666–80. [42] Sahlaoui R, Sab K, Heck JV. Yield strength of masonry-like structures containing thin adhesive joints: 3d or 2d-interface model for the joints? C R Mécanique 2011;339(6):432–8. [43] Spada A, Giambanco G, Rizzo P. Damage and plasticity at the interfaces in composite materials and structures. Comput Methods Appl Mech Eng 2009;198(49–52):3884–901. [44] Toti J, Marfia S, Sacco E. Coupled body-interface nonlocal damage model for FRP detachment. Comput Method App Mech Eng 2013;260:1–23. [45] Wriggers P, Zavarise G, Zohdi TI. A computational study of interfacial debonding damage in fibrous composite materials. Comput Mater Sci 1998;12 (1):39–56. [46] Xu G, Pietruszczak S. Numerical analysis of concrete fracture based on a homogenization technique. Comput Struct 1997;63(3):497–509.