(PDF) A thermodynamics-based cohesive mo

International Journal of Solids and Structures 51 (2014) 647–659 Contents lists available at ScienceDirect International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr A thermodynamics-based cohesive model for interface debonding and friction Irene Guiamatsia a,⇑, Giang D. Nguyen b a School of Civil Engineering J05, University of Sydney, Sydney, NSW 2006, Australia b School of Civil, Environmental and Mining Engineering, University of Adelaide, Adelaide, SA 5005, Australia a r t i c l e i n f o a b s t r a c t Article history: A constitutive model for interface debonding is proposed which is able to account for mixed-mode cou- Received 1 May 2013 pled debonding and plasticity, as well as further coupling between debonding and friction including post- Received in revised form 15 October 2013 delamination friction. The work is an extension of a previous model which focuses on the coupling Available online 11 November 2013 between mixed-mode delamination and plasticity. By distinguishing the interface into two parts, a cracked one where friction can occur and an integral one where further damage takes place, the coupling Keywords: between frictional dissipation and energy loss through damage is seamlessly achieved. A simple frame- Interface work for coupled dissipative processes is utilised to derive a single yield function which accurately cap- Constitutive model Friction tures the evolution of interface strength with increasing damage, for both tensile and compressive Coupling regimes. The new material model is implemented as a user-deﬁned interface element in the commercial Delamination package ABAQUS and is used to predict delamination under compressive loads in several test cases. Thermodynamics Ó 2013 Elsevier Ltd. All rights reserved. 1. Introduction rectly capturing the frictional dissipation and inelastic deformation that accompany the damage process. Drastic improvements in pre- One of the most conspicuous mechanisms of material damage is dicting the permanent set observed experimentally for some mate- the formation of well deﬁned regions where inelastic deformation rials, especially polymeric and cementitious materials can be localises and subsequently propagates. When the size of these re- obtained by coupling damage and plasticity. While most of these gions is small enough (compared to the structural scale), it can models are based on plasticity theory with the loss of stiffness con- be computationally advantageous to idealise the localised damage sidered secondary to the plastic-like behaviour (Tvergaard and zone as an interface; this way, the full continuum constitutive law Hutchinson, 1993; Su et al., 2004; Matzenmiller et al., 2010; Kolluri, is replaced by a formulation relating the interfacial tractions to the 2011), it is our opinion that the reverse (loss of stiffness as the separation between the two surfaces. primary failure process) is more faithful to the physics of interface Examples are laminated composite materials where the excess damage. An energy-based formulation has hence been proposed by resin accumulated between consecutive plies usually fosters the the authors (Guiamatsia and Nguyen, 2012) to capture the initiation of delamination; in the Arctic sea ice, the fracture lines coupling between damage and friction, in a simple and elegant or leads separating blocks of ice that are hundreds of kilometres way rooting from the partition of the total dissipation into in size; in concrete, the width of the damage zone is usually of plastic/frictional and damage ones. However, like many other the order of only a few aggregates in size. In all these cases, the damage–plasticity models, such a model does not have the capabil- computational cost of the numerical model can be reduced by sev- ity to explain the apparent increase in strength and toughness of eral orders of magnitude by using an interface constitutive the material when the loading of the damaging area consists of formulation. the combination of friction and transverse compression. The problem with this simpliﬁcation is that traditional interface Such increase is well documented in a large body of experi- decohesive models (Hillerborg et al., 1976; Xu and Needleman, ments found in the literature. For instance, Wisnom and Jones 1994; Scheider, 2001; Allix and Corigliano, 1996; Camanho et al., (1996) pointed out and carried out experiments to investigate 2003; Davies et al., 2006) tend to focus on the loss of stiffness the role of friction in increasing mode II fracture energy of the (modelled with a damage parameter) and are not capable of cor- interface of laminated composite; de Teresa et al. (2004), within the context of composite delamination, tested cylindrical speci- mens made of laminated composite materials under combined ⇑ Corresponding author. Tel.: +61 9351 2134; fax: +61 2 9351 3343. shear and compression loading and observed a clear increase of E-mail addresses:

[email protected]

(I. Guiamatsia). shear strength with increasing pressure; Lugovy et al. (2005) 0020-7683/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijsolstr.2013.10.032 648 I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 measured increased values of fracture toughness for Si3N4 lami- traction/shear displacement plot measured for the interface under nates under the same loading conditions; ﬁnally, Hallett and Li compression, and (b) the Mohr–Coulomb coefﬁcient of friction de- (2005) performed numerical experiments of the impact of cross- noted l. It is noted that in the previous model being extended here ply carbon ﬁbre laminated composites and concluded that, to (Guiamatsia and Nguyen, 2012), the term ‘‘friction’’ and ‘‘plastic- obtain a realistic prediction, both the fracture toughness and the ity’’ were used interchangeably to denote the partition of the en- strength of delaminating interfaces under compression had to be ergy dissipation resulting in the non-reversible deformation of artiﬁcially increased in the numerical model. the interface. In the present paper, there is a distinction between The easiest way to handle the frictional interaction between that dissipative process linked to residual deformation of the inter- delaminating interfaces is to adopt the very popular two-step face and the sliding friction at microcracks which occurs only when approach consisting of delamination ﬁrst, then frictional contact. the interface is subjected to combined compression and shear. This is the approach taken by Tvergaard (1990) whose model fea- Therefore, ‘‘plasticity’’ is used to refer to the former process, while tures an unphysical unloading–reloading in the constitutive ‘‘friction’’ is reserved for only the second. response if the material is loaded to and past the failure point in The presentation begins with the thermo-mechanical descrip- compressive shear. This is somewhat similar to what can be tion of the model and derivation of the yield function and evolution obtained with standard tools in most commercial ﬁnite element of internal variables. This is followed by the description of the packages: by using element deletion for example, post-damage stress-return algorithm, as implemented in a user-deﬁned inter- frictional contact of the newly created surfaces can be activated; face element with the commercial package ABAQUS/Explicit but the technique is unable to capture the way in which friction (2010). Finally the model is applied to the analysis of a laminated actually affects the damage evolution itself. composite plate under low-velocity impact and the modelling of a A variety of solutions have been proposed by researchers over ﬁbre push out test. the years, two popular ones being (a) the enhancement of Tverg- aard model proposed by Chaboche et al. (1997), where the shear 2. Thermo-mechanical formulation load smoothly decreases to the Mohr–Coulomb limit when under compression, and (b) the (phenomenological) adjustment of the The following expression of the Helmholtz energy potential, W, yield function and/or the fracture toughness to include a normal is considered for the two-phase integral/cracked material: traction term (Matzenmiller et al., 2010; Li et al., 2008; Hou et al., 2001; Christensen and De Teresa, 2004). A comprehensive re- 1 21 2 W ¼ ð1 DÞK n un upn þ D 1 H un upn K n ðun upn Þ view can be found in Raous (2011). 2 2 In 2006, Alfano and Sacco (2006) proposed the representation of 1 2 1 2 þ ð1 DÞK s ðusi upsi Þ þ D 1 H un upn K fs ðusc ufsc Þ the interface as a two-phase material consisting of a damaged part 2 2 and an integral part, hence deﬁning the damage parameter as the ð3Þ proportion of fully damaged interface. Using subscripts c and i For each loading direction (n or s), there are two terms correspond- for ‘cracked’ and ‘integral’, A for the interface area, the damage var- ing to the sum of the cracked (D) and integral (1 D) contributions. iable D is the area ratio as follows: Here, D is a scalar variable representing the interface damage state; Ac u is the vector of interfacial separation, with normal and shear com- A ¼ Ac þ Ai ; with D ¼ ð1Þ A ponents, respectively represented by subscripts n and s; Kn is the elastic stiffness corresponding to the normal or transverse direc- Using the classical mixture theory, the constitutive law of this tion; Ks and K fs are elastic shear stiffness corresponding, respec- ‘composite’ material is obtain through appropriate compatibility tively, to the tensile and compressive loading regimes; ﬁnally, H() relations and mixed stress: is the Heaviside function, taking the value of unity if the argument r ¼ Drc þ ð1 DÞri ð2Þ () is positive, and zero otherwise. The superscripts p and f indicate plasticity and friction, respectively. In the above expression of the It hence becomes straightforward to derive a constitutive model energy potential, it is implicitly assumed that the normal stiffness that seamlessly introduces the effect of compressive friction1 by Kn, once completely lost in tension, is fully recovered upon com- using a traditional frictional contact model for the damaged part of pression. However it is not the case with the shear stiffness Ks: only the material rc. In their model, Alfano and Sacco (2006) use a a fraction can be recovered, e.g. K fs < Ks, depending on the roughness non-associative plasticity model with the Mohr–Coulomb cone as of the cracked surface. Further details on K fs will be elaborated later. the dissipative function. In our model, a fully coupled constitutive It is easy to visualise the compatibility relation between the formulation is entirely derived from standard thermodynamic prin- interfacial displacement jump at the integral part ui and that at ciples. The approach is similar to that used in our previous coupled the cracked part, uc (also cf. Alfano and Sacco (2006)): damage/plasticity interface model, but this time, the expression of the free energy is expanded to include the behaviour of the two- un ¼ uni ¼ unc u ¼ ui ¼ uc ) ð4Þ phase material described above. us ¼ usi ¼ usc The model is unifying in the sense that it is capable of account- Since the distinction has already been established between the ing for damage and plasticity in mixed-mode loading conditions, as inelastic deformation of the integral part being referred to as ‘plas- well as frictional effects on both strength and toughness under tic’ and that of the cracked part as ‘frictional’, the indices i and c can transverse compression. In Fig. 1, we graphically present various be dropped altogether, yielding: scenarios that the new model is capable of capturing. The new model is simple, requiring the calibration of only two 1 1 2 2 parameters in addition to the previous ‘tension-only’ mixed-mode W ¼ ð1 DÞK n un upn þ D 1 Hðun upn Þ K n ðun upn Þ 2 2 delamination model (Guiamatsia and Nguyen, 2012): (a) a 1 2 1 2 compressive (elastic) stiffness, K fs , being the slope of the shear þ ð1 DÞK s us ups þ D 1 Hðun upn Þ K fs us ufs 2 2 ð5Þ 1 The frictional dissipation between two fully damaged surfaces in compression is distinguished from that taking place at an integral interface, and which is modelled In the above expression, the classical additive decomposition of the through as plastic deformation. total jump into elastic and plastic components has been used, with I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 649 ( a ) Damage and ( b ) Da m a g e u n d e r ( c ) Po s t - d a m a g e P l a sticity under transverse friction tr a nsverse tension compression ts ts ts us us us Fig. 1. Range of constitutive responses covered by the proposed formulation. ( ue = u up being the elastic displacement jump. There are four inter- 1 K ðun 2 upn Þ ¼ 2ð1DÞ n t2 ; ðun upn Þ > 0 2 n 2 nal variables: the plastic normal jump upn , the plastic shear jump ups , vn ¼ K n ð13Þ the frictional shear jump ufs and the damage variable D. The corre- 0; ðun upn Þ 6 0 sponding generalised forces, interface tractions tn, tsi, tsc and damage dissipative force v are obtained as derivatives of the free energy 1 t 2si potential (5) with respect to the associated internal variables: vsi ¼ K s us ups 2 ¼ ð14Þ 2 2ð1 DÞ2 K s @W @W @W @W 8 tn ¼ ; t si ¼ ; t sc ¼ ; v¼ ð6Þ @uen @uesi @uesc @D < 0; ðun upn Þ > 0 vsc ¼ 2 t 2sc ð15Þ Expanding, : 12 K fs ðus ufs Þ ¼ f ; ðun upn Þ 6 0 2D2 K s @W @W tn ¼ ¼ ¼ ð1 DÞK n ðun upn Þ Therefore, three terms contribute to the energy driving the damage @uen @ ðun upn Þ process, including a negative contribution from the damaged inter- þ DK n H un upn un upn face in friction. This is consistent with physical observations, as the 1 h 2 2 i (compressive) friction should have the effect of slowing down the þ Ddðun upn Þ K n un upn þ K fs us ufs ð7Þ 2 damage evolution. However, that term must remain smaller or equal to the sum of the ﬁrst two terms, in such a way that the @W @W damage energy is always positive, so as to fulﬁll the irreversibility t si ¼ ¼ ¼ ð1 DÞK s us ups ð8Þ @uesi @ ðus ups Þ condition. For this reason, K fs should be much smaller than Ks (K fs Ks). This also makes general physical sense as the shear stiff- @W @W ness of a cracked interface upon compression must be much smaller t sc ¼ ¼ ¼ DK fs Hð un upn Þ us ufs ð9Þ e @usc @ us uf s than that of the intact material.2 Note that damage stops evolving once the critical value D = 1 corresponding to full loss of cohesion The last term in the expression of the normal traction involves the being reached. Thermodynamically, the damage energy as per Dirac delta, which is inﬁnity when the elastic normal displacement Eq. (12) ceases to exist at that moment. jump vanishes, i.e. uen ¼ un upn ¼ 0. Physically, as the loading re- Similar to the approach taken in Guiamatsia and Nguyen gime goes from tension to compression, the normal ‘elastic’ separa- (2012), a strong coupling between damage, plasticity and friction tion will tend to zero, as should the term us ufs . These terms must is assumed and formulated through the speciﬁcation of damage indeed vanish and remain zero throughout the tensile loading re- potentials as homogenous functions of the rates of internal gime. The result, only at uen ¼ un upn ¼ 0, is hence inﬁnity times variables, after the concept described in Einav et al. (2007). The dis- zero which is indeterminate, but for practical purpose, it will be as- sipation rate is assumed of the following quadratic form, using the sumed that the contribution can be neglected. Consequently, the indicial notation corresponding to the associated internal variable: tractions can be expressed in the reduced form: qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ( U¼ u2D þ up2 n þ us þ us p2 f2 ð16Þ ð1 DÞK n ðun upn Þ; uen >0 tn ¼ ð10Þ K n ðun upn Þ; uen 6 0 where different individual contributions uD , upn , ups and ufs are de- scribed in Table 1. We note the appearance of the additional term ( t si ¼ ð1 DÞK s ðus ups Þ; uen > 0 htni in these expressions to account for the increase of both ts ¼ ð11Þ strength and toughness due to compression. t si þ t sc ¼ ð1 DÞK s ðus ups Þ þ DK fs ðus ufs Þ; uen 6 0 The rate of frictional dissipation, The damage energy is: ufs ¼ ½lhtn i þ X dufs ð17Þ @W 1 2 1 2 v¼ ¼ K n ðun upn Þ 1 Hðun upn Þ K n ðun upn Þ @D 2 2 includes the Mohr–coulomb shear strength, l htni, because we are seeking a yield in the form of the Mohr–Coulomb model for the 1 p 2 1 2 þ K s ðus us Þ 1 Hðun upn Þ K fs us ufs cracked interface. There is also the traction-like term, X, which is 2 2 utilised to impose a condition that the yield function should be ¼ vn þ vsi vsc ð12Þ where 2 The realistic value for K fs is discussed again later in the manuscript. 650 I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 Z D¼1 Table 1 Comparison of the expressions for the dissipation rates between this model and a F n;s ðDÞdD ¼ GðI;IIÞc ð25Þ previous friction-free model. D¼0 This model Guiamatsia and Nguyen Therefore, the mode-mixity parameter b can be tuned with experi- (2012) – k is the mode ratio ments performed with the ratio k held constant, like Reeder and uD v pﬃﬃﬃﬃﬃﬃﬃﬃﬃ vFðDÞ Crews (1990) mixed-mode bending experiments, by using: vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ # dD ð18Þ pﬃﬃﬃﬃ dD u " rD u tr D v n þ 2vsi K s ð1DÞ2 Gc ðk; DÞ FðDÞ p ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 bðkÞ ¼ ð26Þ ð1DÞ 2FðDÞK s þlhtn i ½kGIIc þ ð1 kÞGIc In the model implementation, either a power-law interaction (Whitcomb, 1984) or the BK criterion for composite failure pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃ upn FðDÞ tn FðDÞ t n (Benzeggagh and Kenane, 1996) are used to determine Gc(k, D), pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dupn ð19Þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dupn vð1kÞð1rD Þ but any other function can be utilised for the mixed-mode vn ð1 rD Þ toughness, provided that it corresponds to constant mode-mixity experiments. Further details on the mixed mode interaction can " pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ # pﬃﬃﬃﬃﬃﬃﬃ be found in our earlier paper (Guiamatsia and Nguyen, 2012). ups ð1 DÞ 2FðDÞK s þ lhtn i FðDÞ t s pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dups pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ tsi dups ð20Þ vkð1rD Þ ð1 DÞ 2K s vsi ð1 r D Þ 3. Model behaviour 3.1. Strength u f s ½lht n i þ X dufs ð21Þ N/A With all the dissipative terms being homogeneous functions (in terms of the rates of change of the internal variables) of order 1, the resulting loading function (or yield curve) is obtained, in general- ised stress space, as the quadratic expression, Einav et al. (2007): !2 0 12 0 12 0 12 v t n tsi B tsc C continuous for tn = 0, at the transition of loading regimes, as y ¼ @ uD þ @ @ up A þ @ @ up A þ @ f A 1 6 0 ð27Þ n s @ us described in Appendix A. @dD @dupn @dups f @dus The other dissipation rates (damage, plastic normal, plastic shear) feature the parameter rD ¼ Gcf =Gc ; this is the partition of From the thermodynamic formulation described in Einav et al. the energy dissipation that is allocated to the creation of new (2007), y⁄ is obtained from the Legendre transformation of the dis- surfaces (pure fracture), Gc and Gcf being, respectively, the total sipation potential (16) and plays the role of the plastic potential in interface toughness and the total dissipation from pure fracture3; classical plasticity. It gives the evolution rules for the internal vari- F(D) is a monotonic function controlling the damage dissipation that ables. Expanding and simplifying the above using functions deﬁned veriﬁes: in Table 1 gives the yield condition, denoted here as y, in stress Z D¼1 space as: FðDÞdD ¼ Gc ð22Þ 2 v 2vsi K s ð1 DÞ2 t sc D¼0 y¼ n þ p ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ þ 1 6 0 ð28Þ FðDÞ ð1 DÞ 2FðDÞK þ lht i 2 lhtn i þ X For a given ratio of mode mixity k = vs/v deﬁned as the ratio be- s n tween energy for shear and the total damage energy (see Eqs. This can also be written separately for tensile (+) and compressive (12)–(15)), the (mixed mode) dissipative process controlled by () loading as follows, expanding the damage forces in terms of function F(D) is a combination of functions Fn(D) and Fs(D) in pure damage and respective tractions: modes I and II, respectively (see Guiamatsia and Nguyen, 2012 for 8 2 tn t2 details on the mixed mode formulation and behaviour). In addition, > > 1 þs > < 2ð1DÞ2 Kn Ks for the effects of matrix plasticity, a simple averaging of pure mode I ðþÞ : FðDÞ 1¼0 y¼ 2 ð29Þ (rDI) and pure mode II (rDII) parameters is also assumed: > > t 2si > : ðÞ : pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 þ lhttscn iþX 1 ¼ 0 r D ¼ krDII þ ð1 kÞrDI ð23Þ ð1DÞ 2FðDÞK s þlht n i Continuity is then imposed for tn = 0 (cf. Appendix A) to obtain the FðDÞ ¼ b½kF s ðDÞ þ ð1 kÞF n ðDÞ following expression for the term X that, interestingly, vanishes for where a damage value of either 0 or 1: 2 vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ðN;SÞ2 ðN;SÞ2 u 1 þ 2K u2FðDÞDð1 DÞ2 K s K f X¼t 2K n;s 2 s n;s r DðI;IIÞ GðI;IIÞc ðN;SÞ ð30Þ F n;s ðDÞ ¼ 2 ð24Þ 2ð1 DÞK s þ DK fs ðN;SÞ2 1Dþ ð2K n;s rDðI;IIÞ GðI;IIÞc ðN;SÞ2 Þ Fig. 2 shows isodamage yield curves for an interface with properties Here, classical notation is used for the interface fracture toughness as speciﬁed on the graph, obtained by imposing and keeping con- in normal (GIc) and shear (GIIc) modes, and the strength in normal stant the normal traction tn and loading in shear until the yield con- (N) and shear (S) modes. It can be readily veriﬁed that the integral dition is met. It can be veriﬁed that the curves are indeed of the proposed damage function is equal to the fracture toughness, continuous at the tension/compression transition. It is also noted that the curves in the compressive regime are straight lines, and are parallel to one another, suggesting a variation of the yield shear stress with pressure that is more or less linear: 3 This parameter is identiﬁed experimentally, for each pure mode, from a traction- separation plot as described in Guiamatsia and Nguyen (2012). t;yield s tþ;yield s þ lhtn i ð31Þ I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 651 vs t sc 100 N=60.; S=90. dufs ¼ dD ð40Þ K =K =6000. n s r D ½lht n i þ X 2 t 2sc =0.2;Kf =60. D=0 s 80 The yield equation in compression is used in Eq. (37). From the above, the rate of total energy loss becomes: D=0.5 ts 60 " # vs dD ½lht n i þ X 2 D=0.9 d/ ¼ ð41Þ 40 rD ½lht n i þ X 2 t 2sc In reference to our previous model without friction (Guiamatsia and 20 Nguyen, 2012), the calculated rate of energy dissipation was D=1 simply: 0 -60 -40 -20 0 20 40 60 vs dD FðDÞdD tn d/ ¼ ¼ ð42Þ rD rD Fig. 2. Yield locus for various damage levels. Therefore, the new expression clearly shows that the rate of dissipa- 2 lhtn iþX tion predicted by this model is higher, with the ratio ½l½ht 2 2 n iþX t sc being larger than 1. It is, however, not possible to simplify the 3.2. Toughness expression further, meaning that the rate of dissipation is also dependent, in this case, on the loading path. Considering the same The rate of energy dissipation is obtained through multiplying loading paths utilised to obtain the yield curves of Fig. 2, the total the generalised forces by the rate of change of the associated inter- energy dissipation is numerically integrated, assuming total inter- nal variables. face damage for D = 0.9999.4 For the purpose of comparison, we re- dU ¼ dUD þ dUpn þ dUps þ dUfs ¼ vdD þ t n dupn þ t si dups þ t sc dufs ð32Þ fer to the work of Li et al. (2008), who studied the effect of compressive delamination with traction–separation laws that were Because of the strong coupling used here, the consequence of enhanced with the transverse pressure. Two of their proposed mod- which is a single loading function, the increments of each internal els (A and B) are shown in Table 2, along with the corresponding for- variable are related to a single plastic-type multiplier dk through mula for the damage dissipation. The total energy loss, i.e. the energy the following ’’ﬂow rules’’ obtained from the loading function, dissipated to bring the compressed interface to complete failure, Eq. (27): D = 1, is calculated numerically for several values of transverse pres- 2 3 sure and compared with models A and B in Fig. 3. The sensitivities of @y v r D 4 vn 2vsi K s ð1 DÞ2 the predicted total dissipation with respect to (1) the coefﬁcient of dD ¼ dk @v ¼ 2dk 2 ¼ 2dk þ v FðDÞ ð1 DÞpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ 2 5 @/D 2K s FðDÞ þ lht n i friction and (2) the relative stiffness of the fully damaged partition @dD a ¼ K fs =K s are investigated; in this simulation rD = 1. ð33Þ As seen in Fig. 3(a), the total dissipation does not appear to de- @y tn ð1 rD Þvn pend on the stiffness ratio a, although the partitioning between dupn ¼ dk ¼ 2dk 2 ¼ 2dk ð34Þ damage dissipation and frictional dissipation does, as per @t n @/np FðDÞt n p @dun Fig. 3(b). This is physically reasonable, as smaller a, e.g. lower K fs , leads to lower tsc (see Eq. (9)) and hence delayed frictional sliding @y t sik ð1 r D Þ 2vsi K s ð1 DÞ2 for the same shear displacement under same normal stress. In such dupsk ¼ dk ¼ 2dk 2 ¼ 2dk pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 cases the elastic strain energy vsc will also decrease correspond- @t sik @/sp tsik ½ð1 DÞ 2K s FðDÞ þ lhtn i p @dus ingly, resulting in higher damage energy driving the debonding ð35Þ process (see Eq. (12)) and hence increasing damage dissipation. Therefore it can be said that the ratio a controls the energy parti- @y tsck tsck tion between the damage/plasticity dissipation on one side and dufsk ¼ dk ¼ 2dk 2 ¼ 2dk ð36Þ @tsck @/sp ½lht n i þ X 2 the frictional dissipation on the other, analogous to rD which con- @dup trols the partition between pure damage and plasticity. Therefore, sk the only parameter that controls the total energy dissipated is the In the tensile loading case, the expressions simplify, showing that coefﬁcient of friction l. In other words, if l can be properly mea- increments of all other internal variables are proportional to that sured, the model can predict the effects of compression on the in- of damage, yielding an analytical expression of the total energy dis- crease of both strength and toughness of the interface without the sipation in function of the interface properties. In compression, need of any unphysical tuning parameters. Practically, this friction however, the expressions can only be further reduced as far as parameter l can normally be calibrated experimentally based on the following: yield locii in the compressive section of the stress space such as 2 3 that obtained in the experiments of De Teresa et al. (2004). rD 4 2vsi K s ð1 DÞ2 Fig. 3(a) also shows that the interface toughness in compression dD ¼ 2dk v pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 5 predicted by the current model is higher than that predicted by s ð1 DÞ 2K s FðDÞ þ lht n i " # models A and B, however this increase remains within the same rD t 2sc order of magnitude. The key point is that our model relies only ¼ 2dk 1 ð37Þ on physical basis of debonding and friction upon compression, vs ½lht n i þ X 2 embedded in an energy consistent formulation. The increase in dupn ¼ 0 ð38Þ 4 With the current model, the damage variable appears to increase asymptotically towards the maximum value, and hence a threshold for total damage needs to be ð1 r D Þ vs dups ¼ dD ð39Þ applied. This is also shown in the subsequent section (one element test), with the rD t si asymptotic softening of the shear traction towards the Mohr–Coulomb limit. 652 I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 Table 2 3.3. One element tests Two compression-enhanced laws considered in Li et al. (2008). Model Shape of the compression-enhanced Gc The constitutive model was implemented in ABAQUS/Explicit traction separation law (2010) as a user-deﬁned interface element with the stress-return A Gc 1 l tSn algorithm provided in Appendix B. This is an 8-node directional element (4-node in two dimensions), with upper and lower faces S* clearly speciﬁed through node numbering. Nodal integration scheme was used, to be consistent with ABAQUS generic elements S Gc* for the purpose of comparison in the numerical examples. Single interface element tests with a fully constrained lower face and displacement loading on the upper face (cf. Fig. 4) were used for the validation of the constitutive model implementation. 3.3.1. Cyclic loading test 1 The interface used here had the following properties: 2 B Gc 1 2l tSn þ l tSn Kn = Ks = 0.5 1014 N/m3; GIc = 281 N m/m2; GIIc = 800 N m/m2; S = N = 5 107 Pa; rDI = rDII = 1.0; K fs = 0.5 1012 N/m3; l = 0.3. Fig. 4 S* also illustrates the loading proﬁle consisting of initial compression (u2) that was kept at a constant value u2min while it is being loaded/ S Gc* unloaded/reloaded in shear (u1). Three different levels of compres- sion were applied, namely u2min = 3 107 m, 2 107 m and 0 m, corresponding to normal tractions tn = 0.27 kN, 0.18 kN and 0 kN; the responses are reported in the traction–separation plots of Fig. 5(a). Although the partition of damage dissipation is one in this case, i.e. there is no plastic deformation ups linked to debonding, the curve will not return to the origin upon unloading, because the frictional deformation ufs is also inelastic, as shown. Under compression, the shear behaviour of the interface varies smoothly both strength and toughness is left to the model prediction, from ‘debonding’- dominated to friction – dominated. Under without having to impose any phenomenological rule on the increasing shear loads, the interface shear traction progressively strength and fracture properties of the material model. reduces and converges towards the critical shear traction 15 14 a=0.001, 0.1, 1.; μ=0.4 12 10 a= 1., μ=0.2 10 Model B, μ=0.4 ΦD 8 Φ a = 0.001, 0.1, 1. Model A, μ=0.4 5 6 4 0 2 -60 -40 -20 0 -60 -50 -40 -30 -20 -10 0 tn tn (a) (b) Fig. 3. Variation of the energy dissipation with transverse compression; (a) total dissipation, (b) damage dissipation for l = 0.4. u1(m) u2(m) 6.E-05 0 5.E-05 2 -1E-07 1 4.E-05 u1 u2 3 3.E-05 -2E-07 2.e-05m u2min 3m 2.E-05 00 -3E-07 0. 1.E-05 0.003 m 0.E+00 -4E-07 0.00 0.10 0.20 0.30 0.40 0.50 me(s) Fig. 4. Single element with load (left) and loading history (right). I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 653 be reached for all simulations. A further analysis of the time proﬁle 500 for the case tn = 0.27 kN is compared to the frictionless case, tn=0kN tn = 0 kN, in Fig. 6. 400 tn=-0.18kN Fig. 6 clearly shows that damage is the main source of energy loss at the interface during the initial stages of failure (little gap be- shear tracon (N) tn=-0.27kN tween damage and total curves), but as the damaged area evolves, 300 frictional effects become increasingly important and dominate by the end of the failure process. It is also interesting to note that 200 damage may continue to progress when the shear loading is inverted; this indeed happens for the compressive case shown in 100 Fig. 6, and the damage variable (and dissipation) also increases, al- 0.81kN beit, insigniﬁcantly in this occasion. This response seems sensible 0.54kN as one would expect further damage as a result of continued 0 friction at the incompletely delaminated interface; this is achieved 0.E+00 1.E-05 2.E-05 3.E-05 4.E-05 5.E-05 here thanks to the coupling between damage and friction. shear separaon (m) The effect of varying the compressive stiffness (a ¼ K fs =K s ) and Fig. 5. Shear traction against shear separation for different interfacial pressure partition of damage dissipation (r D ¼ Gcf =Gc was examined in more loads. detail, for the case tn = 0.18 kN, and reported in Fig. 7(a) and (b), respectively. As expected, increasing the plasticity (by reducing rD) causes an t sc ¼ lhtn i, as expected. What was less expected is the asymptotic increase of the permanent deformation. However, increasing the character of that convergence, observed in Fig. 5(a): in fact, whilst compressive stiffness, via increasing ratio a, has a much more under compression, the constitutive model predicts that complete signiﬁcant effect on increasing the permanent deformation (in this damage (D = 1) will actually never happen and this is linked to the case it is ufs ), while the overall stress-separation curves seem assumed form of the dissipative functions. However, since the invariable (as seen in Fig. 7 for values of a sufﬁciently small to com- main features of compressive damage are captured by the model, plete the analysis). This is consistent with the result presented in a practical solution to the ‘asymptotic’ damage evolution is setting Section 3.2 where the frictional dissipation increases for a higher a threshold value of the damage variable, for which the debonding value of the ratio a. Note must be made, however, that the param- is assumed to be complete. In the results presented here, this eter a must not be too high as per the basic model assumption threshold value is 0.9999. K fs Ks. As a rule of thumb, K fs should be at least two orders of Fig. 5 also shows that the energy dissipated (area under curve) magnitude smaller than Ks and in Fig. 7(a), the analyses corre- increases during the progression of coupled debonding/friction sponding to a P 0.1 were not completed due to numerical failure process, as expected, with the applied transverse pressure. difﬁculties. The energy dissipation was also integrated numerically and values The reason for a small K fs is, however, not merely numerical, but of damage dissipation and total dissipation are reported in Table 3. is underpinned by the physical signiﬁcance of our fundamental Note that, for this table and Fig. 6, only one stage of shear loading/ assumption, being the additive decomposition of the interface as unloading was applied at a maximum amplitude of separation of damaged and undamaged parts. It is straightforward to visualise 20 105 m instead of the maximum 5 105 m applied to obtain that the apparent shear stiffness K fs at a fully debonded interface the plot in Fig. 4; this was to allow the threshold damage value to is, in fact, dependent on the applied pressure. If that pressure is relatively small, the force needed to impose a relative sliding Table 3 displacement between the surfaces will also be relatively small, Energy dissipated up to interface debonding, for various levels of pressure (rD = 1). depending on the surface roughness; hence a small K fs . On the other hand, applying a very large pressure between the two sur- Interface Damage dissipation Total dissipation R R R faces crush surface asperities and result, in the limit, to remerging pressure tn (kN) /D ¼ vdD (N m) / ¼ vdD þ tsc dufs (N m) of the two materials, meaning an effective shear stiffness of the 0 799.9 799.9 interface that is of the same order of magnitude as the parent 0.18 2574.1 4941.5 0.27 3511.7 5591.1 material, i.e. K fs Ks (in Alfano and Sacco (2006), K fs = Ks). This latter scenario is, in our view, in conﬂict with the premise that damaged and undamaged zones at the interface can be distinguished, as illustrated in Fig. 8. Since the effects of loading on the magnitude of K fs are not taken into account in this paper, a simple rule 6000 Compression Loading Unloading proposed above is adopted to simplify the implementation. Since 5000 theoretically Ks is very high, we found that the current condition Ф, tn = - 0.27kN on the magnitude of K fs is not restrictive, at least for the examples friction 4000 becomes friction in this paper. during N.mm/mm^2 significant ФD, tn = - 0.27kN unloading 3000 3.3.2. Cyclic loading test 2 Alfano and Sacco (2006) cyclic test is now used to further exam- 2000 full damage ine the capabilities of the current model; in this test, the interface pressure is held constant at tn = 10 MPa and the interface proper- 1000 ФD, Ф, tn = 0 ties are: S = N = 3 MPa; Kn = Ks = 150 N/mm3; GIc = 0.3 N mm/mm2; GIIc = 0.3 N mm/mm2; l = 0.5. We have also used K fs = 0.15 N/mm3. 0 0 0.1 0.2 0.3 It is noted that Alfano and Sacco’s model only utilises an inelas- Time (s) tic displacement jump in the damage part of the interface, which corresponds to our frictional separation ufs . There is no equivalent Fig. 6. Total and damage partition of the energy dissipated for the loading scenario. for our ‘plastic’ separation in which plastic/frictional deformation 654 I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 700 600 600 500 500 rD=1 400 shear traction (N) shear traction (N) a=0.001 400 a=0.01 rD=0.5 a=0.1 300 300 a=1 200 200 100 100 0 0 0.E+00 2.E-05 4.E-05 0.E+00 2.E-05 4.E-05 -100 -100 -200 -200 shear separation (m) shear separation (m) (a) (b) Fig. 7. Shear traction against shear separation. straightforward and very similar to that of a damage/plasticity model, with the addition of the increased yield point. In both cases, it is clear that in our model, the damage component (loss of stiffness) is always present and accompany the other dissipative processes, whereas the loading/unloading curves in Alfano and Sacco’s prediction are more or less parallel to one another and to the elastic loading curve, suggesting that all energy dissipation is linked to the inelastic frictional deformation. 4. Simulation of impact driven delamination Fig. 8. Interface decomposition. Top: under low compressive stress, the interface can be partitioned into pristine (1 – D) and fully damaged (D); Bottom: high The new interface element is now utilised to predict the extent compressive stress cause crushing of asperities to the extent that no separate fully of delamination in laminated composite plates subjected to impact damaged partition can be isolated. loading. The experiment considered is that of Aymerich et al. (2008) to which the reader is referred to for more details on the is coupled with damage during the delamination of the interface, experiment that was performed for several energy levels. For the e.g. due to plastic deformation of the resin or frictional sliding of sake of completeness, key material properties and model parame- micro cracks during delamination (Guiamatsia and Nguyen, ters are presented in Table 4. 2012). Hence the coupling coefﬁcient is initially taken to be The case chosen for the present simulation was the 5.1 J impact, rD = 1, such that, accordingly, no plasticity is present in the inter- which corresponds to an initial velocity of 2108.5 mm/s of the face debonding process. With these settings, the response to the 2.3 kg impactor. The ﬁnite element model, shown in Fig. 10, con- cyclic loading is shown in Fig. 9(a) for two values of the friction sidered only a quarter of the plate, taking advantage of the problem coefﬁcient l. As expected, there is not much residual deformation symmetry; only the two 0/90 interfaces are modelled with as the compressive stiffness is very small. If however, plasticity interface elements, as experimental observations showed the during delamination is coupled into the model, in addition to fric- conﬁnement of delamination to these locations. These interfaces tion from delaminated part of the interface, the prediction is that are designated here by ‘upper’ and ‘lower’ according to their produced in Fig. 9(b), which shows signiﬁcantly more residual distance from the impactor. The upper interface is expected to be deformation. The prediction of our model in these two cases is subjected to compression while the lower one is loaded in tension Fig. 9. Shear traction against shear separation of a single interface. (a) Effect of the friction coefﬁcient, (b) effect of the damage/plasticity coupling. I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 655 Table 4 one) by specifying higher values of shear strength and toughness. Test parameters and material properties, from Aymerich et al. (2008). Such enhancement is purely artiﬁcial and the higher values utilised Interface Thickness: Kn = 8000 MPa/0.02 mm are more or less arbitrarily chosen with the sole purpose of ﬁtting 0.02 mm Ks = 4000 MPa/0.02 mm the experimental results available, hence the models are devoid of q = 1.6109 kg/mm3 any predictive capabilities. GIc = 0.52 N mm/mm2; GIIc = 0.97 N mm/mm2 Here, the following scenarios were examined: Power law failure, g = 1.0 CFRP layer Ply thickness: E11 = 93.7 GPa; E22 = E33 = 7.45 GPa (a) generic interface element utilised as is at both interfaces; 0.125 mm G12 = G23 = G13 = 3.97 GPa (b) generic interface element with higher shear strength m12 = m13 = m23 = 0.261 (100 MPa vs. 80 MPa) and toughness (2.5 N mm/mm2 vs. q = 1.6109 kg/mm3 0.97 N mm/mm2) at the upper interface, elements with nor- Parameters for Hashin yield mal properties at the lower interface; (ABAQUS, 2010): Y1t = 2400 MPa, Y1c = 2000 MPa (c) interface elements with the current constitutive model, Y2t = 100 MPa, Y2c = 300 MPa without friction activated (K fs = 0 MPa/mm, l = 0), used at Y12 = 300 MPa, Y23 = 300 MPa both interfaces; Progressive damage parameters (d) interface elements with the current constitutive model, with (ABAQUS, 2010): friction activated, used at both interfaces (K fs = 40 MPa/ Gc1t = 40 N/mm, Gc1c = 40 N/mm Gc2t = 2 N/mm, Gc2c = 3.5 N/mma 0.02 mm, l – 0). E11 = 93.7 GPa; E22 = E33 = 7.45 GPa Other test CFRP [03, 903]s, Plate dimensions: 45 mm 67.5 mm The predicted delaminated areas on the bottom and upper parameters Impactor: M = 2.3 kg, Diameter = 12.5 mm (spherical) interfaces are compared with the experimental ﬁndings (e) in a Fig. 11. When the upper interface is made artiﬁcially stronger with These values of fracture toughness for matrix in-plane damage are higher than the interface toughness, but were found to yield better correlation with experi- the ABAQUS generic cohesive element (b), the prediction of the mental results in our ﬁnite element analysis. delaminated area at the lower interface is also slightly smaller than the area predicted without such a ﬁx. The new interface model without friction activated (c) yields more or less the same delami- nated area as ABAQUS generic element, which is too large. By using the new model with friction activated (d), a smaller delaminated area is also predicted at both interfaces. If the coefﬁcient of friction is chosen to be l = 0.8, then an excellent match with the experi- mental ﬁndings is obtained, as seen in Fig. 11. It is noted that the partition of damage to total dissipation rD is set to 1 since there is no signiﬁcant plasticity reported in coupon testing also reported in Aymerich et al. (2008). It is worth mentioning, however, that the parameters needed for the frictional part of the model were not validated for the speciﬁc material submitted to the impact test: l = 0.5 is roughly Fig. 10. ABAQUS quarter ﬁnite element model for the impact test. within the range of measurements by Schon (2000), but it was necessary here to use l = 0.8 in order to obtain a good match with experimental results. For a truly predictive model, it is desirable from the bending of adjacent CFRP layers. A common practise that the coefﬁcient of friction be calibrated independently, for when using standard commercial tools to model this sort of impact instance, through combined compression–shear loading as problem is to enhance the interface under compression (the upper suggested in Section 2, although it may also vary in function of (a) (b) (c) (a) 9 (d) 8 (b) (e) Upper interface 7 6 (c) mm 5 (d) 4 3 (e) Lower interface 2 1 0 0 degree line 0 5 10 15 20 25 mm Fig. 11. Numerical predictions vs. experiment: delamination at the lower (plots and damage map) and the upper (damage map only) interface. 656 I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 p MATRIX FIBRE H Fiber rf rm rs Fig. 12. Fibre push out test: experimental setup (left) and axisymmetric ﬁnite element model (right). the transverse compressive traction, which varies during loading. the matrix which is assigned a coefﬁcient of thermal expansion of In addition, we emphasise that the aim of this demonstration is a = 103/°C, resulting in a compressive strain along most of the to achieve a better numerical modelling of the friction enhance- ﬁbre/matrix interface. Then the displacement Dp is applied directly ment to delamination. We do not aim for a high-ﬁdelity modelling onto the ﬁbre nodes in contact with the punch. A relatively ﬁne of damage propagation in composite laminates, in which the inter- mesh was utilised around the interfacial region in both the ﬁbre action between matrix cracking and delamination (Choi and Chang, and the matrix, resulting in a total of 9509 nodes and 9104 1992; de Moura and Goncalves, 2004) must be taken into account. elements. The quasi static analysis was performed with the current model with and without friction activated, as well as with the 5. Simulation of a ﬁbre push out experiment ABAQUS/Explicit solver.5 Fig. 13 shows the ﬁnal deformations ob- tained with both the ABAQUS generic interface element and the The experiment of pushing out a polyester ﬁbre embedded in an user-deﬁned interface element implementing the current constitu- epoxy resin was performed by Bechel and Sottos (1998) and stud- tive model. ied by other researchers (e.g. Hutchinson and Jensen, 1990; Lin The standard (generic) decohesive model is compared with et al., 2001). Fig. 12 illustrates the experimental setup and the ﬁ- the new uniﬁed constitutive model based on the predicted peak nite element model used in this analysis, and the geometric and load. Fig. 14 shows the various predictions for the total reac- material parameters are as follows: tion force recorded at the punch against its downward displacement. rf = 0.95 mm; rs = 1.025 mm; rm = 4.3 mm; H = 5.36 mm It can be noted that the prediction of the new model with Ef = 2500 MPa; mf = 0.35; Em = 4000 MPa; mm = 0.33 no friction activated is in good agreement the prediction of GIc = GIIc = 0.11 N/mm; Linear failure criterion g = 1.0 the ABAQUS generic interface element; this is a good veriﬁca- S = N = 22 MPa tion test. For these two analyses, the delamination is driven Kn = 8000/0.001 MPa/mm only by the relative shear displacement of the ﬁbre respective Ks = 2000/0.001 MPa/mm to the matrix, without consideration of friction. The introduc- K fs = 20/0.001 MPa/mm tion of friction through the coefﬁcient l results in an increase rD = 1.0 of the predicted peak load beyond the point of ﬁrst decohesion (approximately 210 N). In this instance, it is found that taking To simulate the reported initial compressive matrix strain of the value of l = 0.1 for the coefﬁcient of friction yields the best 0.0022, an initial temperature ﬁeld of 2.2 °C was imposed onto match to the experimental ﬁnding, with the peak load of about 400 N predicted to occur at a punch displacement of 0.15 mm. The coefﬁcient of friction was again tuned to obtain the best prediction, similar to the impact-driven delamination test, but should really be determined from independent experiments that were not available for the speciﬁc material system used here. In addition, as can be seen in Fig. 14, the effects of fric- tion on both the strength and toughness of the interface cannot be captured using ABAQUS generic elements with post failure friction. Our proposed model therefore provides a uniﬁed and consistent approach to treat both pre failure (increase of strength and toughness in compression) and post failure fric- tional effects. 5 Note that, only in this ﬁbre push-out example, for the simulation with the ABAQUS generic element to proceed in the post-peak stage, it is necessary to delete the failed cohesive elements, otherwise, severe distortion and hourglassing of the Fig. 13. Deformed ﬁnite element model; left: ABAQUS generic element, right: User- reduced integration continuum elements may cause the simulation to abort deﬁned element with new model. prematurely. I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 657 700 250 600 Generic, No Fricon New Model, No Fricon 200 500 New model, mu=0.1 New Model, mu=0.3 150 400 Exp 100 300 load (N) 50 load (N) 200 0 100 0 0.05 0.1 0.15 0.2 -50 0 Generic, no fricon -100 0 0.05 0.1 0.15 0.2 Generic, mu=0.1 -100 -150 Generic, mu=0.3 -200 -200 -300 displacement (mm) -250 displacement (mm) Fig. 14. Load against displacement for the ﬁbre push out test. Experimental results from Lin et al. (2001). 6. Conclusion In the examples presented, the coefﬁcient of friction was simply tuned until a good match with experiments was achieved. This is This paper presented an interface constitutive model coupling because the experimental results used were taken from the exist- damage, plasticity and friction within a consistent thermodynamic ing literature and no coefﬁcient of friction was provided for the framework. Two main ingredients constitute the basis of the mod- speciﬁc material systems. As the range for the coefﬁcient of friction el: (1) an expression for the free energy that uniﬁes the three can be rather wide and strongly depend on the interfacing materi- sources of dissipation and couples the evolution of friction to that als, an experimental programme needs to be developed for further of damage by separating a unit interface area into distinct validation and we are aiming to address this aspect in future work. damaged area undergoing friction and integral area undergoing We also acknowledge the strong assumption made in the damage; and (2) assumed functions of the damage dissipation formulation of this model, which is the lumping of the interface potentials that accounts for the effect of compression on the evolu- roughness into a single stiffness parameter K fs . While the demon- tion of all internal variables. These ingredients were introduced strations in this study show the usefulness and practicality of this into the existing framework developed by Einav et al. (2007), lead- assumption, further development to properly take into account the ing to a set of constitutive equations that are able to capture the effects of (cracked) surface roughness on both compressive shear observed behaviour of interfaces under a wide spectrum of loading stiffness and strength, e.g. (Serpieri and Alfano, 2011) is also scenarios, in particular, the increased yield point and total dissipa- planned for the next steps. tion observed when the interface is loaded in combined compres- sion and shear. Acknowledgments Three key parameters control the response of the constitutive model: The authors gratefully acknowledge the support of funding from the Australian Research Council discovery grant DP1093485. – The friction coefﬁcient, l, which directly determines the com- pressed interface shear strength. This total shear strength is Appendix A. Derivation of X(D) equal to the sum of the shear strength under zero compression and the equivalent critical shear strength of the Mohr–Coulomb Enforcing the continuity of the yield curve between tensile (+) cone (ltn). and compressive () side of the stress space gives: – the compressive shear stiffness K fs of the debonded part of the 2 interface, which controls the amount of inelastic frictional 2vsi K s ð1 DÞ2 tsc 2vþsi K s ð1 DÞ2 y¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ2 þ X 1¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ2 1 deformation and hence directly controls the partition of energy ð1 DÞ 2FðDÞK s ð1 DÞ 2FðDÞK s dissipation by friction. 2 v tsc vþ – The partition of fracture dissipation which is the ratio rD of pure ) si þ ¼ si damage (loss of stiffness) dissipation to the cumulative dissipa- FðDÞ X FðDÞ tion due to damage and plasticity under tensile loading only. In In tension, the maximum shear traction is found for zero normal this work, it was assumed that rD = 1 as the emphasis was traction: placed on illustrating the ability of the interface model to cap- 2 vþsi ðt þsi Þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ture the response under compressive loading. ¼ ¼ 1 ) t þs ¼ t þsi ¼ ð1 DÞ 2K s FðDÞ FðDÞ 2 2ð1 DÞ K s FðDÞ The constitutive model was implemented as a user-deﬁned That maximum must be identical to that in the compressive regime element in a commercial ﬁnite element package and utilised for at zero normal traction, i.e. t þ s ¼ t s ¼ t si þ t sc . the prediction of a number of test cases with readily available Rewriting the yield function in compression: experimental data. It was shown that toggling on the frictional 2 2 dissipation capability resulted in predictions that were more faith- vsi tsc ðt si Þ2 t sc t þ ¼1) 2 þ ¼ 1 ) sc ful to experimental results and expected trends. The key advanta- FðDÞ X tsi þ t sc X X ges of this model are (a) the thermodynamic framework used in its vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ u development, from which all the constitutive relations are consis- u ðt Þ2 tsc ¼ t1 si 2 ) X ¼ rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ tently derived with a single yield curve and (b) the meaningfulness t si þ tsc ðt Þ2 1 si 2 of the model parameters which have a clear physical signiﬁcance ðtsi þtsc Þ and can be straightforwardly calibrated from relatively simple t si pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ experiments. Using tsc ¼ ð1DÞ D Ks f and t si þ t sc ¼ ð1 DÞ 2K s FðDÞ gives Ks 658 I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 tsc t si þ tsc t si þ tsc where the individual components of increments of plastic and fric- X ¼ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sit þ t ¼ r ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ q ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sc 2t tion shear separation are, with j = 1, 2: ðt sc Þ 2t si þ tsc ð2tsi þtsc Þ t si þ1 t sc sc pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ t jsi tjsc ð1 DÞ 2K s FðDÞ ð1 DÞ 2K s FðDÞ dupj p s ¼ dus ; dufjs ¼ dufs ¼ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ) X t si tsc 2ð1DÞ 2ð1DÞK s þ1 f þ1 D DK s vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Step 2 is repeated iteratively until the value of the yield is within a u certain accuracy tolerance, |y| < tol. u2FðDÞDð1 DÞ2 K s K f ¼t s 2ð1 DÞK s þ DK fs References ABAQUS, Simulia Inc., 2010. Theory Manual, version 6.9. Appendix B. Stress-return algorithm Alfano, G., Sacco, E., 2006. Combining interface damage and friction in a cohesive zone model. International Journal for Numerical Methods in Engineering 68, The tractions tn, t1si , t 2si , t1sc , t 2sc , as well as the damage D are history 542–582. Allix, O., Corigliano, A., 1996. Modeling and simulation of crack-propagation in variables used to track the state of the interface and a ﬂag Comp- mixed-modes interlaminar fracture specimens. International Journal of Fracture Flag is used to determine whether the interface is in tension or 77, 111–140. compression. Note the two components of shear traction in three Aymerich, F., Dore, F., Priolo, P., 2008. Prediction of impact-induced delamination in cross-ply composite laminates using cohesive interface elements. Composites dimensions. Given an incremental displacement of the form Science and Technology 68, 2383–2390. d = (dn, d1s , d2s ), the stress is updated according to the following Bechel, V., Sottos, N.R., 1998. Application of debond length measurements to algorithm. examine the mechanics if ﬁber pushout. Journal of the Mechanics and Physics of Solids 46 (9), 1675–1697. Benzeggagh, M.L., Kenane, M., 1996. Measurement of mixed-mode delamination B.1. Step 1: Elastic predictor fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Composites Science and Technology 56 (4), 439–449. The normal (transverse) component is updated ﬁrst to assess Camanho, P.P., Davila, C.G., de Moura, M.F., 2003. Numerical simulation of mixed- mode progressive delamination in composite materials. Journal of Composite the loading regime: tensile or compressive. Materials 37, 1415–1438. Chaboche, J.L., Girard, R., Schaff, A., 1997. Numerical analysis of composite systems t n ¼ tn þ ð1 DÞK n dn by using interphase/interface models. Computational Mechanics 20, 3–11. Choi, H.Y., Chang, F.K., 1992. A model for predicting damage in graphite/epoxy Here, the prescript ‘‘’’ indicates the sate variable at the previous laminated composites resulting from low-velocity point impact. Journal of time step. The ﬂag C for compression is set (to 1) if the predicted nor- Composite Materials 26 (14), 2134–2169. mal traction is negative t n < 0, otherwise the ﬂag is zero. The predic- Christensen, R.M., De Teresa, S.J., 2004. Delamination failure investigation for out- of-plane loading in laminates. Journal of Composites Materials 38, 2231– tors for the three interface tractions (in 3D, j = 1, 2) are then further 2238. updated and the shear resultants for the integral and cracked parts Davies, G.A.O., Hitchings, D., Ankersen, J., 2006. Predicting delamination and are computed as the norm of their individual components. debonding in modern aerospace composite structures. Composites Science and Technology 66, 846–854. t n ¼ t n þ ½ð1 CÞð1 DÞ þ CK n dn De Moura, M.F.S.F., Goncalves, J.P.M., 2004. Modelling the interaction between matrix cracking and delamination in carbon-epoxy laminates under low t jsi ¼ t jsi þ ð1 DÞK s djs velocity impact. Composites Science and Technology 64 (7–8), 1021–1027. De Teresa, S.J., Freeman, D.C., Groves, S.E., 2004. The effects of through-thickness t jsc ¼ tjsc þ C DK fs djs compression on the interlaminar shear response of laminated ﬁber composites. Journal of Composite Materials 38, 681–697. and the damage energy is also calculated: Einav, I., Houlsby, G.T., Nguyen, G.D., 2007. Coupled damage and plasticity models " # " 2 2 # " 2 2 # derived from energy and dissipation potentials. International Journal of Solids ðt n Þ2 t 1si þ ðt 2si Þ t 1sc þ t 2sc and Structures 44 (7–8), 2487–2508. v¼ þ C Guiamatsia, I., Nguyen, G.D., 2012. A generic approach to constitutive modelling of ð1 DÞ2 K s ð1 DÞ2 K s ð DÞ2 K fs composite delamination under mixed-mode loading conditions. Composites Science and Technology 72, 269–277. Hallett, S.R., Li, X., 2005. A numerical investigation of dynamic transverse shear B.2. Step 2: Yield veriﬁcation and ﬂow rules failure. In: Proceedings of the 15th International Conference of Composite Materials, Durban, South Africa. Hillerborg, A., Modeer, M., Petersson, P.E., 1976. Analysis of crack formation and The value of the yield, y, is then calculated using Eq. (28). If it is crack growth in concrete by means of fracture mechanics and ﬁnite elements. violated, then the ﬂow rules are calculated and the state variables Cement and Concrete Research 6, 773–782. Hou, J.P., Petrinic, N., Ruiz, C., 2001. A delamination criterion for laminated updated using the increments of internal variables as per Eqs. composites under low-velocity impact. Composites Science Technology 61, (33)–(36) with the plastic multiplier given by: 2069–2074. y Hutchinson, J.W., Jensen, H.M., 1990. Models of ﬁber debonding and pullout in dk ¼ 8 9 brittle composites with friction. Mechanics of Materials 9, 139–163. > @y @tsi v > < @D þ @t@yn @tn @D þ @t@ysi @D þ @t@ysc @tsc @D ð@/D =@dDÞ2 þ = Kolluri, M., 2011. An in-situ experimental–numerical approach for interface delamination characterization. Ph.D. thesis. Eindhoven University of 2 @y tn @t si t si t sc > : @tn @t n p 2 þ @t@ysi p 2 þ @t@ysc @t sc f 2 > ; Technology. @us ð@/pn =@dupn Þ @us ð@/ps =@dups Þ @us ð@/fs =@dufs Þ Li, X., Hallett, S.R., Wisnom, M., 2008. Predicting the effect of through-thickness compressive stress on delamination using interface elements. Composites: Part Finally, the state variables are updated: A 39 (2008), 218–230. k Lin, G., Geubelle, P.H., Sottos, N.R., 2001. Simulation of ﬁber debonding with friction kþ1 tn in a model composite pushout test. International Journal of Solids and t n ¼ k t n ð1 CÞ dD þ ð1 DÞK n dupn Structures 38, 8547–8562. 1D Lugovy, M., Slyunyayev, V., Orlovskaya, N., Blugan, G., Kuebler, J., Lewis, M., 2005. ( ) k j Apparent fracture toughness of Si3N4-based laminates with residual kþ1 j tsi t si ¼ k t jsi dD þ ð1 DÞK s dupj s compressive or tensile stresses in surface layers. Acta Materialia 53, 289– 1D 296. k j Matzenmiller, A., Gerlach, S., Fiolka, M., 2010. A critical analysis of interface t sc constitutive models for the simulation of delamination in composites and kþ1 j t sc ¼ k tjsc þ dD DK fs dufjs failure of adhesive bonds. Journal of Mechanics of Materials and Structures 5 D (2), 185–211. kþ1 D ¼ k D þ dD Raous, M., 2011. Interface models coupling adhesion and friction. Compte Rendus Mecanique 339, 491–501. I. Guiamatsia, G.D. Nguyen / International Journal of Solids and Structures 51 (2014) 647–659 659 Reeder, J.R., Crews, J.H., 1990. Mixed-mode bending method for delamination Tvergaard, V., 1990. Effect of ﬁbre debonding in a whisker-reinforced metal. testing. AIAA Journal 28, 1270–1276. Materials Science and Engineering A 125, 203–213. Scheider, I., 2001. Cohesive model for crack propagation analyses of structures with Tvergaard, V., Hutchinson, J.W., 1993. The inﬂuence of plasticity on mixed mode elastic–plastic material behaviour; foundations and implementation. GKSS interface toughness. Journal of the Mechanics and Physics of Solids 41, 1119– research center Geesthacht, Report. 1135. Schon, J., 2000. Coefﬁcient of friction of composite delamination surfaces. Wear 237, Whitcomb, J.D., 1984. Analysis of instability-related growth of a through-width 77–89. delamination. NASA TM-86301, National Aeronautics and Space Serpieri, R., Alfano, G., 2011. Bond-slip analysis via a thermodynamically consistent Administration, Washington, DC. interface model combining interlocking, damage and friction. International Wisnom, M.R., Jones, M.I., 1996. Measurement of friction in mode II delamination Journal for Numerical Methods in Engineering 2011 (85), 164–186. with through-thickness compression. In: 7th European Conference on Su, C., Wei, Y.J., Anand, L., 2004. An elastic–plastic interface constitutive model: Composite Materials. London, UK. application to adhesive joints. International Journal of Plasticity 20 (12), 2063– Xu, X.P., Needleman, A., 1994. Numerical simulations of fast crack-growth in brittle 2081. solids. Journal of the Mechanics and Physics of Solids 42 (9), 1397–1434.

(PDF) A thermodynamics-based cohesive model for interface debonding and friction