A 3D two-scale multiplane cohesive-zone model for mixed-mode fracture with finite dilation

Accepted Manuscript A 3D two-scale multiplane cohesive-zone model for mixed-mode fracture with finite dilation R. Serpieri, M. Albarella, E. Sacco PII: DOI: Reference:

S0045-7825(16)31354-8 http://dx.doi.org/10.1016/j.cma.2016.10.021 CMA 11182

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date: 7 February 2016 Revised date: 14 September 2016 Accepted date: 8 October 2016 Please cite this article as: R. Serpieri, M. Albarella, E. Sacco, A 3D two-scale multiplane cohesive-zone model for mixed-mode fracture with finite dilation, Comput. Methods Appl. Mech. Engrg. (2016), http://dx.doi.org/10.1016/j.cma.2016.10.021 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript Click here to download Manuscript: Revised_Manuscript.pdf

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Click here to view linked Referenc

A 3D two-scale multiplane cohesive-zone model for mixed-mode fracture with finite dilation R. Serpieri1 , M. Albarella Dipartimento di Ingegneria, Universit` a degli Studi del Sannio, Piazza Roma, 21 - I. 82100, Benevento, Italy.

E. Sacco Universit` a di Cassino e del Lazio Meridionale, Dipartimento di Ingegneria Civile e Meccanica, Cassino(FR), Italy.

Abstract A 3D multi-scale cohesive-zone model (CZM) combining friction and finite dilation by a multi-plane approach (M-CZM), based on the concept of Representative Multiplane Element (RME), is developed within the mechanics of generalized continua for the analysis of mixed-mode fracture. The proposed M-CZM formulation captures the increase of measured fracture energy in mode II as a natural effect of multi-scale coupling between cohesion, friction and interlocking, employing a reduced set of micromechanical parameters characterized by a well-defined micromechanical interpretation. This permits to devise clear calibration and identification procedures for 3D fracture problems. Upon assessing the retrieval, by a regular 5-plane RME, of a quasi-isotropic response for fracture resistance and for dilation the M-CZM is employed in FEM simulations of Double-Cantilever Beam (DCB) tests to obtain predictions of mixed mode I-II and mixed mode I-III fracture resistance. The DCB analyses show the key role of the characteristic height of asperities in determining the macroscopic fracture resistance in both mixed mode I-II and I-III interactions. Numerical results also show the independence of the mode-I fracture resistance on the geometry of the beam section and a marked dependence of the measured mixed-mode fracture resistance on the section aspect ratio. Keywords: three-dimensional cohesive-zone model; interlocking; finite dilation; mixed-mode delamination; damage-friction coupling.

1

Corresponding author

Preprint submitted to Computer Methods in Applied Mechanics and Engineering September 14, 2016

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1. Introduction Structural interfaces are mechanical models adopted in a wide range of engineering applications to reproduce the physical response of thin layers subjected to high strain gradients or to simulate the interaction between bodies in contact. They are generally constructed considering two surfaces at zero distance, whose kinematics is governed by the relative displacement field. Nonlinear localized phenomena, such as formation of shear bands, de-cohesion, detachment and delamination, fracture and contact, are successfully reproduced with the aid of interfaces, also for three-dimensional problems. In particular, among several proposed 3D cohesive interface models, Ortiz and Pandolfi [1] developed a finite-deformation cohesive element within a class of irreversible cohesive laws for the analysis of the dynamic crack growth. On this basis, Pandolfi et al. [2] developed numerical simulations to study dynamic fracture in C300 maraging steel. A 12-noded triangular interface element has been implemented in [3] in the framework of a large-displacement formulation accounting for cohesive crack propagation in mode I and in mixed modes; the effectiveness of the developed finite element has been validated by studying damage due to particle fracture and interface decohesion in composites. Finite element simulations have been performed using 3D cohesive interface elements to investigate the effect of triaxiality of the stress state occurring at the 3D crack tip of a fractured body in [4]. Non-planar 3D interfaces have been also embedded into linear elastic finite elements, adopting the unit partition technique to simulate quasi-static crack growth and propagation [5]. In several cohesive interfaces the damaging response is often accompanied by significant frictional effects influencing damage evolution and the overall nonlinear mechanical behavior. In these problems, a proper account of the presence of interface roughness is demanded in order to retrieve reliable predictions of damage evolution. Examples of structural applications in which coupling of cohesion and friction appears to be relevant are masonry panels [6, 7, 8] and metal-ceramic composites [9]. Also in cross-ply laminates, rough delamination patterns can be observed, characterized by an oscillating microstructure with a typical length-scale [10] (as shown in Figure 1). These patterns are also observed in Central Cut Ply (CCP) tests [11] on fiber-reinforced laminates, and non planarity of fracture surfaces has been recognized to play a role in mixed-mode resistance measured by Double Cantilever Beam (DCB) tests with Uneven Bending Moments (UBM) [12]. The interplay between these roughness-induced features and the cohesive response raises several research issues. From the point of view of computational mechanics, cohesion-roughness interaction poses the problem of devising affordable and mechanically reliable modeling strategies capable of capturing the complex interplay among friction, damage and the geometry of asperities at different 2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 1: Typical oscillating delamination path observed in cross-ply laminates. (left): view of fracture surface; (right): side view. Photo reprinted from Brunner et al., (2008) [13], with permission from Elsevier.

length-scales. From the point of view of material testing, this interaction introduces a difficulty in the definition of standards for mixed-mode delamination resistance characterization, where a well-known debated issue is the meaningfulness of mode I and mode II fracture energies as real material properties of a given structural interface [13], considered that friction strongly influences the quantity of dissipated energy before collapse. In particular, concerning cohesive-frictional interface models, Raous et al. [14] proposed a thermodynamically consistent interface model accounting for unilateral contact, nonassociative Coulomb friction and adhesion. A micromechanical formulation coupling damage and friction has been presented by Alfano and Sacco [15] considering an ideally flat microstructure where unilateral contact and friction are originated in the micro-cracked part of a representative area element. This idea has been extended to the case of hydraulic crack growth in [16] and, in a more complex context, in [17]. Del Piero and Raous [18] also derived an interface model able to account for decreasing of adhesion intensity, friction and viscosity, combining the first two laws of thermodynamics by phenomenological assumptions. More recently, the behavior of interfaces subjected to shear has been also modeled employing associative plasticity [19]. Dilation and interlocking phenomena due to the presence of asperities can also significantly influence the mechanical response of the interface. Several models accounting for dilation have been proposed in the literature, based on phenomenological or micro-mechanical considerations. Structural interface models for geomaterials and concrete accounting for dilation of the crack surface have been proposed in [20], and in [21] considering the presence of sawtooth asperities. An interface model accounting for unilateral contact, friction and dilation based on a phenomenological decomposition of the friction angle in components responsible for dissipation and for dilation is proposed in [22] to analyze the stability of masonry elements interconnected by dry joints. Dilation is modeled by Mr´ oz and Giambanco [23] assuming an interaction among spherical asperities with variation of the contact area and a characteristic distinction between sliding strains and 3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

micro-slip strains. Models of friction and micro-dilation for fiber-matrix degradation of composite materials have been presented by Stupkiewicz [24] while the inclusion of de-cohesion evolution and damage have been subsequently examined by Stupkiewicz and Mr´oz in [25] with special reference to cyclic loading. More recently, dilation has been phenomenologically included by a non-associative rule and governed by two characteristic mechanical angles [26] to reproduce the response of mortar in masonry assemblages, while a micromechanical constitutive model for simulating the behavior of fully formed cracks in geomaterials has been proposed by Mihai and Jefferson in [27]. The description of friction and interlocking mechanics by purely-macroscale interface models, which are typically computationally efficient, is based on a high idealization and simplification of the geometry of contacting surfaces involved at the meso- and micro-scale of the fracture process zone and requires often the addition of heuristic corrections to find an improved correlation with experimental data [28]. By a different perspective, multiscale approaches to modeling of microstructured interfaces pursue the description of the mechanics of fracture, friction and contact resorting to a higher degree of geometric detail of the interface mesostructure, considering differentiated material models for matrix, inclusions (i.e., aggregates or fibers) and for the relevant interfacial transition zone, with the objective of obtaining homogenized macroscopic stress-relative displacement relationships. The related computational methods are typically based on algorithms generating regular or random aggregate meso-structures, see, e.g., [29], upon introducing some degree of upfront idealization of the mesostructure geometry and of the distribution of aggregates. In particular, Palmieri and De Lorenzis [30], following an approach analogous to Matouˇs et al. [31], Alfaro et al. [32] and Verhoosel et al. [33], employ Hill’s principle and a periodic 2D cell to retrieve homogenized stress-relative displacement curves for a concrete layer idealized assuming plane-stress conditions with aggregates modeled as randomly shaped polygons, using the 7-parameter isotropic damage law by Mazar [34] for cement paste and for the interfacial transition zone. In the above quoted papers, 2D models are considered, and friction and finite dilation are not included in the analyses. Account of 3D kinematics and friction is considered instead in [35] where, although not directly focused on the determination of macroscopic cohesive laws, a multiscale approach is deployed for the 3D meso-structural analysis of concrete specimens, albeit without considering finite dilation. As reported, the consideration of 3D representative elements involves a high computational demand and major efforts of mesh generation, code efficiency and post-processing. It emerges from the above quoted works that obtaining realistic simulations of the macroscopic 3D cohesive-frictional response of structural interfaces by multi4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

scale analyses employing such a high degree of geometric detail of the mesostructure appears to be not indicated for ordinary analyses of 3D delamination and debonding problems in presence of macroscopic structures having a conspicuous number of integration points. The use of these approaches appears even more problematic when finite dilation is considered, due to the challenging computational demand implied by the use of contact-mechanics models at the mesoscale in view of the large number of potentially contacting surfaces. It can be concluded that realistic multiscale modelling of finite dilation in 3D interface problems appears to still represent a challenge for the computational analysis techniques and resources currently available, at least for analyses involving large macroscale structure with a large number of integration points. From the above surveyed literature it emerges that, while the formulation of damage, frictional and dilating behaviors in cohesive zone models has been broadly investigated, and a large library of computational methods is also available for three-dimensional cohesive zone modeling, these research problems have been mainly examined in separate form. The formulation of numerical methods implementing macroscale 3D cohesive zone models for microstructured interfaces comprehensively combining damage, friction and, in particular, finite dilation, seems to have been by far less investigated. A suitable ground for the investigation of the feasibility of the inclusion of all the above listed mechanical features into a comprehensive 3D cohesive-zone formulation by a simplified multiscale approach may be represented by the CohesiveZone Model (CZM) framework proposed in [36, 37, 38], hereby abbreviated in Multi-Plane Cohesive Zone Model (M-CZM). M-CZM formulations are based on the description of the small-scale geometry by means of elementary plane patterns (RME) in which the response of each elementary plane is defined exploiting a further decomposition into an undamaged and a fully damaged part with individual free energy functions and inelastic variables. The use of a simplified RME description of the interface geometry makes the M-CZM approach an intermediate multiscale option between macroscale CZMs and full multiscale approaches, where the number of elementary planes can be determined as a compromise between the sought level of detail in representing damage distribution over the different material space orientations and the computational effort. Previous investigations [36, 37, 38, 39, 40] have shown that, although the geometric resolution is very limited compared to full multiscale approaches, M-CZM captures the essential dissipation processes during crack initiation and propagation in quasi-brittle interfaces, by adding a reduced number of degrees of freedom associated with each integration point, compared to multiscale techniques. M-CZM shares with multiscale approaches the employment of material parameters having a precise physical meaning, which simplifies their identification based on the macroscale or microscale available data. In particu5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

lar, in [39] a sensitivity analysis with respect to the key model parameters defining the RME geometry is carried out, and numerical-experimental assessments are reported concerning the simulation of tests by Lee et al. [41] on granite joints subjected to cyclic shear loading, and of pull-out tests of a steel bar from a concrete block, by Shima et al. [42]. These assessments have indicated that, while M-CZM is of course not able to capture the entire complex damage processes occurring at the mesoscale and microscale, the essential qualitative aspects of the dissipative processes simulated can be captured with a quantitative accuracy that appears satisfactory, compared to the overall uncertainty that affects part of the experimental measures. Moreover, the reported numerical-experimental comparisons are based on upfront declared calibration procedures deliberately devised to leave the recourse to heuristic assumptions or curve-fitting to minimum, thus permitting an easier independent validation of the predicting capabilities of this approach in relation to the specific engineering problems at hand. The above referenced works have separately addressed finite dilation in planestrain 2D kinematics [39] and have presented a preliminary assessment of the capability of a 3D M-CZM to reproduce in-plane isotropy of fracture energy at the single Gauss point level [40], although without accounting for finite dilation. Inclusion of finite dilation into a 3D M-CZM cohesive-zone model is not trivial, and raises specific computational and theoretical issues. In particular, the identification of an optimal strategy for retrieving different fracture energies in mode I and in mode II by a thermodynamically consistent approach suitable for mixed-mode cyclic analyses turns out to be a still debated research problem, see, e.g., [43]. Specific issues posed by the three dimensionality of the formulation are a proper description of in-plane isotropy of dilation and the assessment of the predicted response for both mode II and mode III . The objective of the present study is the formulation and the development of suitable computational methods for the implementation in FEM analyses of a class of M-CZMs capable of effectively combining friction and finite dilation into a 3D framework. This study covers the assessment of the predicted response in structural FEM simulations of mixed-mode decohesion tests throughout the entire mode-mixity spectra I-II and I-III. The work is organized as follows. Section 2 describes the finite-dilation 3D M-CZM, its general hypotheses rooted in the thermo-mechanics of generalized continua with internal variables, placing emphasis on the concepts of rupture energy and total measured fracture energy. Details on the finite-step relationships required for the numerical integration of the M-CZM and on the associated algorithmic tangent operator are reported in Appendix A. The assessment of the predicted mechanical response is carried out at the single Gauss-point level in Section 3. Finally, in Section 4, structural simulations of mixed mode I-II and mixed mode I-III delamination in DCB-UBM tests are presented. 6

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

2. Finite-dilation 3D M-CZM formulation The M-CZM formulation is a class of CZMs for structural interfaces herein described as proposed in previous works [36, 37, 38, 39, 40]. The objective of the M-CZM formulation is to provide a simple multiscale mathematical representation of the coupled inelastic phenomena occuring during decohesion processes due to the interaction between the multiplicity of material surfaces which can be observed at different length-scales in the fracture process zones of quasi brittle solids, and which typically have different space orientations. Definition of Representative M-CZM Element The description of a Representative M-CZM Element (RME) consisting of a finite number of elementary planes and state variables is detailed. The geometry of a RME associated with a macroscopic point of a structural interface is specified by the unit normal vector N defining the average surface of the interface, plus a set of Np elementary planes, each having an individual orientation in space defined by the unit vector n(k) , and effective area A(k) , with k = 1, . . . , Np . The projected area of the k-th elementary plane on the average surface is A¯(k) = A(k) n(k) · N = A(k) cos θa(k) ,

(1)

(k)

where θa is the relevant inclination angle of the elementary plane above the ¯ and the average surface. Accordingly, the total projected area of the RME, A, total effective area, A, of the RME are given by: A¯ =

Np X

A¯(k)

A=

,

Np X

A(k) .

(2)

k=1

k=1

In the formulas above and henceforth the notation convention of marking all quantities related to projected areas by an overbar is adopted. Figure 2 shows an example of a RME composed of a pattern of Np = 5 microplanes. The symbols TX and TY are used to denote two orthogonal unit vectors associated with the average surface such that TX , TY and N form a left-handed reference-frame. The symbols h(k) and t(k) denote two orthogonal unit normal vectors lying on the elementary plane k, which form together with n(k) a left-handed local reference-frame of such elementary plane. Excluded the case of plane k oriented as the average surface (i.e., n(k) =N), h(k) and t(k) are chosen by the following rule: h(k) =

N × n(k) k N × n(k) k

,

t(k) = n(k) × h(k) .

(3)

Conversely, in the trivial case n(k) =N, local axes are set coincident with global ones, i.e., t(k) =TX and h(k) =TY . 7

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Geometrical and mechanical quantities associated with the k-th elementary plane are denoted by a round bracketed superscript (.)(k) . A sample RME is shown in Figure 2 whose first elementary plane has the same orientation of the TX -TY plane while the remaining elementary planes (2) to (5) are oriented so as to generate a truncated pyramid shape with square basis. (5)

n N

t

(5)

h

Ty

(5)

(3)

Tx t

(2)

n

(1)

n

(2)

n

(3)

h

(1)

h

(4)

n

(2)

h

t

(4)

h t

t

(1)

(3)

(4)

Figure 2: Example of RME geometry composed of 5 elementary planes

On the basis of the representation of the RME geometry above introduced, the free energy of the whole RME, F (RM E) , equals the sum of free energy of each elementary plane composing the RME, viz: F

(RM E)

=

Np X

(4)

F (k) .

k=1

Relation (4) is written in terms of densities of free energy as follows ψ¯A¯ =

Np X

A¯(k) ψ¯(k) =

k=1

Np X

A(k) ψ (k) ,

(5)

k=1

where ψ¯ = F (RM E) /A¯ is the density of overall RME energy per unit interface projected area, whereas ψ¯(k) = F (k) /A¯(k) and ψ (k) = F (k) /A(k) are the densities of free energy of the k-th elementary planes per unit projected area and per unit effective area, respectively. Introducing the dimensionless area fraction coefficients A¯(k) A(k) A(k) A γ¯ (k) = ¯ , γ (k) = , γ˜ (k) = ¯ = γ (k) ¯ , (6) A A A A on account of (1) and (2), one has trivially Np X k=1

γ¯ (k) = 1,

Np X k=1

γ (k) = 1, 8

Np X k=1

γ˜ (k) ≥ 1,

(7)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

PNp (k) γ˜ = 1 only in the case of a perfectly flat RME having n(k) · N = 1 with k=1 for any k. Recalling equations (1) and (6), the balance (5) can be equivalently also written in terms of area fractions in any of the following two forms ψ¯ =

Np X

γ¯

(k)

ψ¯(k) =

k=1

ψ¯ =

Np X

Np X

γ¯ (k)

k=1

γ˜

(k)

ψ

(k)

k=1

=

Np X k=1

1 (k)

cos θa

ψ (k)

A γ (k) ¯ ψ (k) A

(8)

(9)

General hypotheses for constitutive modelling In the M-CZM approach, the free energy balance, expressed in any of the forms shown in (8) and (9), constitutes the first basic state relation for combining different constitutive responses defined at the level of the elementary planes into a macroscopic constitutive response of the whole RME. A second relation of kinematic type employed in the M-CZM formulation hereby considered is the assumption that the small scale mechanical response associated with a macroscopic point x of the cohesive surface depends only on the rigid translational part of the small scale field of the relative displacements between the joined surfaces, which equals the macroscopic displacement ¯s attained in x, denoted by an overbar. Accordingly, the kinematic state of the RME is assumed to be completely defined, for constitutive purposes, by the macroscopic relative displacement vector ¯s which equals the relative displacement s(k) of each elementary plane k, viz.: s(k) = ¯s. (10) Relation (8), (or (9)) and (10) together constitute the basic relations coupling the constitutive responses of the elementary planes. In the following, the global normal component of ¯s is denoted by symbol s¯N , i.e. s¯N = ¯s · N and the local (k) normal component of s(k) relative to the elementary plane k is denoted by sn , (k) i.e. sn = s(k) · n(k) (with lowercase subscript). Concerning the constitutive hypotheses of the M-CZM formulation, it is admitted the possibility of identifying for each elementary plane k a set of constitutive internal variables specifically pertaining to such k-th plane alone. Adopting the constitutive scheme described above, a cohesive behavior is introduced employing a set of scalar dimensionless variables α(k) , k = 1 . . . , Np with α(k) representing the damage of plane k. This quantity admits the micromechanical interpretation, largely exploited in damage mechanics [44, 45], of being (k) the ratio between the damaged area Ad of the subset of the k-th elementary plane where bonds are fully broken over its total effective area, viz.: (k)

α

(k)

Ad = (k) , A 9

(11)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(k)

with 0 ≤ Ad ≤ A(k) . The range of these variables is 0 ≤ α(k) ≤ 1 where the attainment of the lower and upper thresholds corresponds to the base condition of null damage in plane k and to the achievement of complete decohesion in plane k, respectively. Full damage in the RME corresponds to the achievement of condition α(k) = 1 for all k. An hypothesis of local inelastic evolution on the elementary planes is applied according to which, denoting by I (k) the set of all remaining internal inelastic variables of the RME (other than α(k) ) governing the response of the damaged part of elementary plane k, the constitutive response of each plane can be more generally defined in terms of free energy state functions having the special form:     (12) ψ¯(k) = ψ¯(k) I (k) , α(k) , s(k) , ψ (k) = ψ (k) I (k) , α(k) , s(k) ,

and in terms of complementary evolution laws for the inelastic variables of plane k having the form:   (13) I˙ (k) = f I (k) , s(k) .

Specific hypotheses for the cohesive response Damage is formulated in the framework of the Thermo-Mechanics of Generalized Continua with internal variables (TMGC) [46]. In particular, for damage evolution, general hypotheses analogous to those considered in [37] are also herein introduced. These are hereby recalled with reference to plane k: • Coupling of damaged and undamaged response by a convex combination For the specific constitutive response of each individual elementary plane the strategy proposed by Alfano and Sacco [15] is applied for partitioning the elementary plane into an undamaged part with elastic behavior and a damaged inelastic part by using a convex combination in α(k) :    
(k) (k) (k) (k) ψ (k) α(k) , I (k) , s(k) = ψu α(k) , Iu , s(k) + ψd α(k) , Id , s(k) =

  
(k)  (k) (k) (k) 1 − α(k) ψus Iu , s(k) + α(k) ψds Id , s(k) (14) (k) (k) where ψu and ψd are, respectively, the free energy density contributed to ψ (k) by the undamaged and the damaged part of the elementary plane (k) (k) k, and ψus and ψds are their respective effective counterparts. The set of (k) (k) inelastic variables, I (k) , is divided into two subsets, Iu , Id , respectively associated with the mechanical response over the damaged and undamaged parts.

• Normal law for damage evolution and uncoupling with the remaining inelastic variables 10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

An uncoupled evolution law between damage and the other inelastic variables is considered and is formulated in agreement with the TMGC framework. Accordinlgy, denoting by X (k) the damage-driving force of plane k (k) associated with the free energy of the undamaged zone ψu (k)

X (k) = −

∂ψu (k) = ψus ∂α(k)

(15)

the evolution of α(k) is stated with the aid of a complementary evolution law for X (k) which, expressed in dual form in terms of a force potential ϕ(k) , admits the representation: α˙ (k) ∈

∂ϕ(k) . ∂X (k)

(16)

According to this hypothesis the potential ϕ(k) is only function of the
thermodynamic force X (k) and damage α(k) , viz.: ϕ(k) = ϕ(k) X (k) , α(k) .

• Rate-independent damage behavior A rate independent behavior is considered for ϕ(k) . Such behavior is described choosing for α(k) a characteristic
(k) (indicator) function IP of the interval Pα(k) = [0, X0 α(k) ], [47]:    0 if X (k) ∈ Pα(k) (k) (k) . (17) = X ϕ +∞ otherwise Concerning the assumption of the previous bullet point, the threshold value (k) (k) for the damage stress X0 is a function only of α(k) , this function, X0 =
(k) X0 α(k) , is defined so as to modulate a selected softening response for the elementary plane.

The rupture energy density per unit effective area of the elementary plane k, is defined as the free energy spent by decohesion forces in bringing α(k) from 0 to 1 (complete damage): Z (k) Gc = X (k) dα(k) (18) (k) Gc ,

Γ (k)

with Γ (k) defining a curve in the s(k) − α(k) plane that starts from an undamaged configuration P0 = (s(k) = s0 , α(k) = 0) and ends into a full damaged configuration P1 = (s = se , α(k) = 1). The curve Γ (k) can be partitioned into the union (ki) of curves ΓU characterized by a stationary damage condition, with α˙ (k) = 0, (ki) (or unloading branches) and branches ΓL characterized by a growing damage condition, α˙ (k) > 0, (loading branches) where, owing to (16) and (17) (k)

X (k) (s(k) ) = X0 (α(k) ). 11

(19)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

On account of this decomposition, Equation (18) yields Z Z Z (k) (k) (k) (k) (k) X dα = S X dα + S X (k) dα(k) = Gc = (ki) (ki) Γ i ΓU i ΓL Z 1 XZ XZ (k) X0 dα(k) . X (k) dα(k) = X (k) dα(k) + = i

(ki)

ΓU

(ki)

0

ΓL

i

(20)

(k)

The last equality in Equation (20) highlights the univocal character of Gc as a
(k) definite integral only depending on function X0 α(k) . An expressive representation is given to this integral applying to the integrand in the last term of (20) the following identity  i0    0 h   (k) 0 (k) (k) (k) (21) + X0 1 − α(k) 1 − α(k) X0 = −X0 1 − α(k) = − X0

with prime indicating differentiation with respect to α(k) . Accordingly, one computes Z 1  dX (k) (k) (k) 0 (22) 1 − α(k) Gc = X0 (0) + dα(k) . dα 0 (k) (k) Since X0 (0) equals the free energy at the onset of damage, ψu (k) , it is α

(k)

=0

recognized that if X0 is constant with respect to α(k) , which corresponds to a (k) completely brittle behavior, then Gc simply coincides with the free energy at (k) (k) the onset of damage, viz.: Gc = ψu . α=0 Averages of rupture energy per unit projected area are introduced by considerations analogous to those previously applied for free energy potentials. In ¯ (k) particular, the rupture energy density per unit projected area of plane k, G c , is defined as: (k) ¯ (k) = G(k) A . (23) G c c A¯(k) Recalling (1) and (2) one has ¯ (k) G c =

1 (k) cos θa

(24)

G(k) c

and relations analogous to (8) and (9) hold between the rupture energies of ele¯ c: mentary planes and the overall rupture energy density G ¯c = G

Np X k=1

¯ (k) = γ¯ (k) G c

Np X k=1

12

γ¯ (k)

1 (k) cos θa

G(k) c ,

(25)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

or, even ¯c = G

Np X

γ˜ (k) G(k) c =

k=1

Np X k=1

A γ (k) ¯ Gc(k) . A

(26)

(k)

(k)

For simplicity, it is hereby assumed that Gc are uniform, i.e., that Gc In such a case, recalling (7), relation (26) yields ¯ c = A Gc . G A¯

= Gc .

(27)

Moreover, since the following identity holds γ˜ (k) G(k) c =

A A(k) (k) A(k) A (k) Gc = Gc = γ (k) ¯ G(k) , ¯ ¯ A A A A c

(28)

(k)

a further consequence of the uniformity of Gc , inferred from (27) and (28) is that ¯ c. γ˜ (k) Gc = γ (k) G (29) 2.1. Specific constitutive assumptions for elementary planes The specific constitutive response of the individual elementary planes is now described by detailing the   expressions  considered  for the associated free energy (k) (k) (k) (k) (k) (k) and ψds Id , s of a given plane k. Modelling functions ψus Iu , s choices are analogous to those applied in previous works [37, 39, 40] accounting for unilateral contact, friction and interlocking, under the assumption that the individual elementary planes are perfectly flat. The unilateral character of the response is analytically described with the aid of the unit step function (or Heaviside function) u(·) and of the positive and negative parts of x, hxi+ , hxi− , defined as:

u(x) =



1 x>0 , 0 x<0

hxi+ =

1 (x + |x|) , 2

hxi− =

1 (x − |x|) , 2

(30)

and with the aid of the associated positive and negative ’unilateral projectors’ ˆ (+) ˆ (−) P e1 , Pe1 associated with a direction of unit vector e1 . These are vector to vector applications such that applied to a test vector v return the vector which has the following components in a orthonormal frame e1 , e2 , e3 : h i  t ˆ (+) P (v) = hv · e1 i+ , v · e2 , v · e3 (31) e1 h

i  t ˆ (−) P (v) = hv · e1 i− , v · e2 , v · e3 e1 13

(32)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Response over the undamaged part Over the undamaged parts of the elementary planes a linear elastic behavior is (k) considered. Accordingly, dependence on inelastic variables is ruled out, Iu = ∅, (k) and ψus is defined as follows   1 (k) (33) ψus s(k) = Ks(k) · s(k) , 2

where K is a stiffness operator defined by a constant 3 × 3 matrix [K] of stiffness coefficients. For simplicity this operator is the same for all the planes. Also on account of the hypothesis of perfectly flat elementary planes the related matrix in the local reference frame of plane k is assumed diagonal: [K] = diag(Kn , Kt , Kt ),

(34)

where Kn and Kt are normal and tangential penalty stiffness coefficients. Damage evolution The remaining constitutive assumptions for damage evolution are introduced
(k) by specifying K in (33) and function X0 α(k) entering the TMGC framework previously stated (Equations (16) and (17)). Herein modelling choices analogous to those in [37] are followed, with the purpose of defining a simplest possible cohesive response under the assumption that elementary planes are perfectly flat. Accordingly, it is set K(k) = KI with K = Kn = Kt being the value of a penalty modulus taken equal for tangential and normal stiffness. This choice corresponds to a direct dependence of the strain energy over the undamaged part upon the

2 displacement norm s(k) , viz.:

2 1

(k) ψus = K s(k) . (35) 2 (k)

Function X0 is defined in order to retrieve a simple linear softening branch. Its analytical expression is accordingly equal to   2 s0 sf    1 K α<1 (k) X0 α(k) = (36) 2 sf − α(k) (sf − s0 )  ∞ α=1

where s0 is the displacement at the onset of the damage and sf is the displacement corresponding to complete rupture. Relations (35) and (36) yield that, under monotonic decohesion paths, the same bilinear displacement-stress constitutive diagrams, shown in Figure 3, are retrieved in local mode I and mode II over the individual elementary plane, with the following relationship between stiffness and fracture energy: s0 2Gc , η =1− , (37) K= 2 sf (1 − η)sf 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where η is a ductility parameter such that a unit value reproduces a rigid-cohesive behavior with infinite initial stiffness.

ʍ ƌĞĂс'Đ ʍϬ

Ɛ ƐϬϭ

ƐĨϭ

Figure 3: Mode I bilinear stress-displacement law of damage response.

It is important to recall that, as shown in [37], the existence of a relativedisplacement energy of the form (35) implies the existence of a displace density

ment scalar norm, s(k) α(k) driving damage evolution, which defines iso-damage curves in the relative-displacement space, s(k) . In particular, the evolution law of variable α(k) , determined by Equations (15), (16),(17) and (36), coincides with the condition for α(k) to achieve the following maximum value

oo n n

(k) (k) α = max 0, min 1, s (k) , (38) history

α

where the norm recovering (36) turns out to be a 3D generalization of the norm considered in [36, 40]: !

(k) 1 β

(k) , (39)

s (k) = η 1 + β (k) α

with β (k) defined as

β

(k)

=

ˆ (+) kP s(k) k n(k) s0

− 1,

(40)

It is also worth to recall also from [37], that, in this equivalent restatement of (15), (16),(17) the unilateral character of the loading-unloading conditions, represented by the inequalities entering the usual Karush-Kuhn-Tucker conditions [48, 49] (KKT), is determined by the completely equivalent, yet computationally simpler, optimality condition (38). 15

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

ˆ (k) s(k) in (40) accounts for the adSimilar to [36, 40] the projected term P n ditional assumption that damage is insensitive to normal compressive displacement components. The existence of such a norm is particularly convenient in displacement-driven integration algorithms since it allows updating of damage variables with a direct computation from incremental displacements, [36, 40]. Response over the damaged part The response over the damaged part implements unilateral contact, friction and interlocking. The frictional law is described exploiting a standard formulation (k) of nonassociative plasticity [50], defined by a free energy function ψf : (+)

(k)

ψf

    1 1 (k) (k) (k) (k) = H(k) sde · sde = H(k) s(k) − sdi · s(k) − sdi 2 2

(41)

ˆ (−) where H(k) = P K is a stiffness operator that accounts for (abrupt) loss of n(k) frictional energy as a result of contact loss over the individual elementary plane k. (k) Evolution of sdi is described in the framework of nonassociated elastoplasticity [15, 50] by introducing a threshold function of Coloumb type: (k)

(k)

(k)

φ(σ f ) = µhσf n i + k τ f k, −

(42)

and a nonassociated plastic potential function, which determines sliding in the tangential direction of the elementary plane k (k)

(k)

g(σ f ) =k τ f k, (k)

(k)

(k)

(k)

(43)

(k)

where σf n = σ f · n(k) and τ f = σ f − σf n n(k) are the normal and the tangential component, respectively, of the frictional stress contributed by the damaged part of the elementary plane k, (k)

(k)

σf =

∂ψf

. (44) ∂s(k) The plastic flow rate equation is expressed with the aid of a nonnegative plastic multiplier λ(k) ∂g s˙di (k) = λ˙ (k) (k) (45) ∂σ f and by introducing Kuhn-Tucker conditions [49] to address the irreversible nature of frictional sliding: λ˙ (k) ≥ 0

,

(k)

φ(σ f ) ≤ 0

,

(k) λ˙ (k) φ(σ f ) = 0.

(46)

When the finite height of asperities is not considered, the stress acting over the (k) (k) damaged surface of the k-th elementary plane, σ d = ∂ψd /∂s(k) , consists only 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(k)

(k)

of the frictional contribution, σ d = σ f , and the formulation so far described achieves the 3D extension of the one examined in [37]. Since an indefinitely dilating behavior is acceptable only when relative displacements are small compared to the characteristic height of asperities, herein, the effect of the finite length of the asperities is implemented employing the same strategy proposed in [39]. Accordingly, a dimensionless projection function Ag is introduced which accounts, at a global level, for interlocking loss outside of an ¯ N (which is interaction zone set equal to the characteristic height of asperities H a global parameter). This function has the following analytical expression:       s¯N s¯N Ag ¯ = 1− ¯ . (47) HN HN + + (k)

Following [39] a modified expression for σ d is considered   s¯N (k) (k) σf . σ d = Ag ¯ HN

(48)

Macroscopic Fracture energy and Rupture energy The inclusion of friction adds a further mechanical process which dissipates free energy. Considering a relative-displacement history Γ (k) such that the initial configuration is undamaged and the final one completely damaged, the total spent (k) mechanical work density over the k-th elementary plane, GT , along Γ (k) Z (k) GT = σ (k) · ds(k) (49) Γ (k)

is thus balanced by the (path independent) rupture energy, plus the additional (k) (k) dissipation due to friction and the change in stored elastic energy (ψ1 − ψ0 ), over this elementary plane, viz.: Z (k) (k) (k) (k) (k) (k) GT = Gc + α(k) σ d · dsdi + ψ1 − ψ0 . (50) Γ (k)

Considering the typical condition in complete decohesion tests such that, for any plane (k), Γ (k) starts and ends with zero stored free energy in the interface (k) (k) (ψ0 = 0, ψ1 = 0), then the total spent mechanical work density over the entire ¯ T , along the global deformation path Γ (which RME, per unit projected area, G (k) due to (10) is such that Γ = Γ ) turns out to be equal to ¯T G

=

Np X k=1

(k)

γ˜ (k) GT =

Np X

γ˜ (k) G(k) c +

k=1

Np X k=1

17

γ˜ (k)

Z

Γ (k)

(k)

(k)

α(k) σ d · dsdi .

(51)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(k)

When Gc

(k)

is uniform (Gc Np X

= Gc ), then due to (7) and (29),

γ˜ (k) G(k) c =

k=1

Np X

γ˜ (k) Gc =

k=1

Np X

¯c = G ¯ c, γ (k) G

(52)

k=1

¯ T is recognized to be equal to the rupture energy, G ¯ c , plus the contribuand G tion of parasitic forces, possibly originated during decohesion by the interplay of friction and interlocking, viz.: ¯T G

¯c + = G

Np X

γ˜

(k)

Z

Γ (k)

k=1

(k)

(k)

α(k) σ d · dsdi .

(53)

¯ T is termed total measured fracture energy. This new denominaIn this case G ¯ c, tion is deliberately used to remark the difference with the rupture energy, G ¯ T may be taken as a macroscopic measure of the and remark that, although G ¯c energy required to produce decohesion, it is a path-dependent quantity while G is instead path independent. 2.2. Summary of material parameters and stress relations Stress relations The stress relations resulting from the formulation detailed in the previous sections are hereby summarized in view of their numerical integration carried out in the next section. Denoting by σ (k) = ∂ψ (k) /∂¯s the stress contributed by elementary plane (k) ¯ s the overall stress, the latter can be computed on account of and by σ = ∂ ψ/∂¯ either (8) or (9). In particular use of (9) yields for the overall interface stress the expression: Np X σ= γ˜ (k) σ (k) , (54) k=1

while decomposition (14), on account of the constitutive assumptions of Section 2.1, yields   ∂ψ (k) s¯N (k) (k) σ (k) = = (1 − α(k) )σ (k) + α A (55) g u ¯ N σf , ∂¯s H

with and

(k) σ (k) u = Ks

(56)

  (k) (k) σ f = H(k) s(k) − sdi .

(57)

18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

To achieve simpler relations in the computations carried out in the next section the following normal/tangential decompositions of stress quantities are also considered (k) (k) (k) σun = σ u · n(k) , τ (k) (58) u = σ u − σun · n , (k)

τ d = σ d − σdn · n(k) ,

(k)

τ f = σ f − σf n · n(k) ,

σdn = σ d · n(k) , σf n = σ f · n(k) .

(k)

(k)

(59)

(k)

(k)

(60)

As already pointed out in previous works in relation to the 2D formulation[37, 39], it is important to remark that also for the 3D model, although as previously observed the equations of Section 2.1 determine the same bilinear law on the individual elementary planes, this condition, which stems as a thermodynamic requirement, does not limit the ability of the multiplane model to describe different macroscopic mode I and mode II responses. In particular, as shown in these references, coincidence of fracture energies in modes I and II, over the individual elementary plane, is the thermodynamic requirement implied by the use of an associated type of damage evolution law with a single damage variable, combined with the choice of having a threshold damage function only depending on α(k) , and combined with the employment of an equivalent-displacement norm. As shown in the quoted references and by the numerical results in the next sections, difference between macroscopic mode-I and mode-II fracture energies emerges instead as the joint effect of a nonzero friction angle and of the presence inclined elementary planes. Similarly, coincidence of initial normal and tangential stiffnesses of the elementary planes entails no relevant limits to the macroscopic response describable by the present model given the penalty character of these stiffnesses. Material parameters Considering uniform properties over the elementary planes, the material parameters which have to be assigned in order to specify a particular interface response within the M-CZM formulation so far detailed are hereby summarized together with their physical meaning: • the RME geometry is specified by the sequence of n(k) and by any of the sequences of parameters γ¯ (k) , γ (k) or γ˜ (k) (with k = 1, . . . , Np ); • fracture energy is specified either microstructurally, by Gc , or by the overall ¯ c; global rupture energy density G • the definition of the cohesive response is completed by specifying any two of parameters s0 , sf , η, K, σ0 (with σ0 = Ks0 being the stress at damage onset); 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

• friction is specified by the tangent of the friction angle µ over the elementary planes; ¯N . • interlocking is specified by the global height of asperities H The previous bullet list summarizes the relatively small number of parameters which have to be identified to specify a given interface response in the M-CZM approach, under the hypothesis of uniform material parameters over the elementary planes. It should be remarked that while, on the one hand, the identification of the parameters entering the M-CZM formulation in relation to a given engineering problem is made easier as each of them has a very precise physical meaning, on the other hand, their number (and the associated computational demand) depends linearly from the number of elementary planes of the RME. Previous investigations have shown that in the M-CZM formulation the RME geometry primarily determines the qualitative and quantitative features of the macroscopic interface response [37, 39, 40]. As the choice of the RME is necessarily the result of a trade-off between the sought level of accuracy in the predicted mechanical response and the computational demand, these preparatory investigations have permitted to fix some elements for a proper RME selection and for the calibration of the associated material parameters in relation to some classes of problems of simulation of structural interfaces. These elements are summarized below: • The procedure for the calibration of the RME parameters depends on the type of identification problem at hand: specifically on the available measures of the mode-I and mode-II macro-scale response as well as on the information available on the small-scale properties and on the small-scale structure of the interface. In particular, macroscopic decohesion tests pro¯ c while the consideration of the vide data suitable for the identification of G (k) distribution of Gc over the elementary planes (or the common value Gc in presence of uniform material properties) requires a detailed knowledge of the microstructure. • Concerning the selection of the geometry of the RME for the retrieval of a macroscopic interface response with isotropy in the tangential plane, a dedicated study [40] has shown that a quasi-isotropic response can be conveniently approximated by employing an RME geometry corresponding to the lateral surface and upper base of a truncated pyramid surface with regular polygonal base, assuming equipartition of γ (k) (i.e., γ (k) = 1/Np , in agreement with (7) ) and uniform material properties over the elementary planes, since the computed degree of anisotropy becomes negligible as the number of sides of the polygon is increased. In this case the selection of the 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

geometry of the RME is dictated by the maximum level of anisotropy that is deemed acceptable for the application at hand. A 5-plane RME with square base was found to be suitable for several applications. • As it can be easily checked, when material parameters are specified using the ¯ c , recollecting stress relations (54)macroscopic fracture energy density G (57) and recalling (29), (35) and (37), one obtains the following most explicit macroscopic expression for global interface stresses in terms of parameters ¯ c: γ (k) and G " Np   ¯c X 2G (k) (k) σ= s(k) + α(k) Ag 1 − α γ (1 − η)s2f k=1

s¯N ¯ (k) H N

!

ˆ (−) P n(k)



(k) sd

−

(k) sdi



(61) ¯ In particular, parameters Gc , σ0 are obtained from measurements in pure mode-I experiments of the fracture energy and of the stress at damage on¯ N and other geometrical properties of the set, respectively. Parameter H RME, such as the inclination of the elementary planes are inferred from dilation diagrams in mode-II test. Parameter η has to be set to a sufficiently high value in order to not introduce a spurious deformability in the global structural model, yet not too high to avoid ill-conditioning. Parameter µ governs friction and is calibrated based on measurements of tangential and normal stresses in mode II tests. Concerning this parameter, it is important to point out the distinction between the microscopic friction µ of the model and the measured macroscopic friction µ ¯ which is defined as the ratio of global tangential stress to global normal stress after complete decohesion, the latter being the result of an interplay between friction and geometry of the asperities. In particular, simple equilibrium considerations on tangential and normal forces on a flat elementary plane (k) with inclination angle (k) θa yield the following relation between µ and µ ¯: (k)

µ ¯

(k)

=

(k)

µ cos θa + sin θa (k)

(k)

cos θa − µ sin θa

.

(62)

3. Single point assessment Results of numerical analyses of the mechanical response determined by the formulation described in Section 2 at the level of the single interface point are hereby illustrated. These computations were carried out employing the return mapping scheme and the tangent operator reported in Appendix A, and have confirmed robustness of the algorithm with quadratic convergence as well as a 21

#

.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

limited request of computational resources when a reasonable number of elementary planes is employed. Results are hereby reported with the objective of analyzing, under complete decohesion tests, the sensitivity of the interface response to material parameters, with a specific view to the assessment of such response along different space orientations and under different loading conditions. In particular, space orientation of loading and of the computed effects are specified in terms of azimuth angle θL (which is the angle contained in the plane TX -TY of a given space orientation measured from the TX -axis), and zenith angle, θN (measuring the inclination of the orientation off the global TX -TY plane), as shown in Figure 4.

Figure 4: Sketch of the RME geometry employed in single-point simulations and of angles defining test orientation; θa : inclination angle from plane TX -TY (common to all inclined planes); θN : zenith angle of test orientation; θL : azimuth angle of test orientation

Two types of analyses have been considered: analyses of complete displacementcontrolled decohesion, in which all three components of displacement are prescribed along a predetermined orientation in space, while the emerging three global stress components σN , τX , τX are the measured output data, and confinedslip decohesion analyses in which the two components of displacement are prescribed along a direction contained in the TX -TY plane, while the interface is subjected to normal confinement kept to a fixed value σN , and dilation s¯N in the N direction is unconstrained. For confined-slip tests only angle θL is specified. In all single-point simulations the 5-plane RME whose geometrical parameters are specified in Table 1 is used. This geometry, is defined by the inclination angle θa from plane TX -TY , which is the same for elementary planes 2-5, and is sketched by the truncated-pyramid shape with square basis in Figure 4. This RME geometry is defined based on a simplest possible equipartition hypothesis of area fractions and space orientations and is specifically selected since, as shown in [40], it provides a reasonable approximation of isotropy even with a limited number of elementary planes. When not otherwise specified, mechanical param22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

eters used in single-point analyses are those reported in Table 2.

(k) nN (k) nTX (k) nTY γ (k)

k=1

k=2

k=3

k=4

k=5

1.0

cos θa

cos θa

cos θa

cos θa

0.0

sin θa

0.0

0.0

0.0

0.0

− sin θa 0.0

sin θa

0.2

0.2

0.2

0.2

− sin θa 0.2

Table 1: 5-plane RME pattern (normals and area fractions)

Section 3.1 shows analyses of sensistivity of the stress-displacement response ¯ N when the interface undergoes tangential and to the height of the asperities H normal decohesion paths. Section 3.2 reports a sensitivity analysis of the mea¯ T to the friction angle µ, to θa and to H ¯ N . Section 3.3 sured fracture energy G finally shows an assessment on the degree of anisotropy of the model in the TX TY plane. σ0 [MPa] 3.0

¯ c [KJ/m2 ] G 0.3

η [–] 0.9

µ [–] 0.5

θa [◦ ] 40.0

Table 2: Material parameters employed in the single point simulations.

3.1. Sensitivity of the stress-displacement response to the height of the asperities A first set of confined-slip decohesion simulations were performed to investi¯ N . Figure 5a shows the stress-displacement gate the sensitivity of the model to H curves that are obtained by impressing a displacement-controlled decohesion in ¯ N in the interval 0.001 mm to 1.000 mm. The plot pure mode I, and spanning H shows that, the Mode I bilinear curve defined by equations (40), previously shown ¯ N and that the fracture energy is exin Figure 3, is recovered independent of H ¯ actly Gc . Figure 5b shows the tau-slip curves that are obtained by a confined-slip test in direction θL = 0 while applying a confinement pressure of σN = −1.0 MPa and spanning the height of the asperities in the interval 0.001 mm to 1.000 mm, ¯ N has an influence on the tau-slip response. The corresponding showing that H dilation curves are shown in Figure 6.

23

(a)

(b) 3.5

3

H = 1.000

H = 1.000

N

N

2.5

3

HN = 0.500 HN = 0.250

2

HN = 0.500 HN = 0.250

2.5

HN = 0.100

H = 0.100 HN = 0.050

1.5

τ [MPa]

σ [MPa]

N

H = 0.010 N

HN = 0.001

1

HN = 0.010

0

0.5

0.05

0.1

0.15

0.2

0.25

0.3

0

0.35

HN = 0.001

1.5

1

0

HN = 0.050

2

0.5

−0.5

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

st [mm]

s [mm] N

Figure 5: (a) Stress-dilation curves in a complete displacement controlled in the normal direction ¯ N spanned in the interval 0.001 mm to 1.000 mm. θN = 90◦ with parameters in Table 2 and H (b) Tau-slip curves in a confined-slip decohesion analysis with parameters in Table 2, direction ¯ N spanned in the interval 0.001 mm to 1.000 mm. of the displacement θL = 0 and H

(a)

(b)

1.2

0.3

HN = 1.000 1

HN = 1.000 0.25

HN = 0.500

sN [mm]

HN = 0.050 HN = 0.010 HN = 0.001

0.4

HN = 0.010

0

0

1

s [mm] t

1.5

2

HN = 0.050 HN = 0.001

0.1

0.05

0.5

HN = 0.100

0.15

0.2

0

HN = 0.250

0.2

HN = 0.100

0.6

−0.2

HN = 0.500

HN = 0.250

0.8

sN [mm]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

−0.05

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

st [mm]

Figure 6: (a) Dilation curves in a confined-slip decohesion analysis with parameters in Table 2, ¯ N spanned in the interval 0.001 mm to 1.000 direction of the displacement θL = 0 and with H mm. (b) Magnification of the diagrams contained in the dashed box in the figure (a)

¯ N of displacement controlled Mode I/Mode II interac3.2. Sensitivity to µ and H tion Complete displacement-controlled decohesion simulations were performed by applying displacements in a direction with a fixed θL and a variable θN and by ¯ T during decohesion. The interface paramcomputing the total fracture energy G eters are again those of Table 2. Polar plots are used to display the results of mode I-mode II interaction as function of θN . Sensitivity to friction and to asperity angle is illustrated by Figure 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

7a which shows the polar plots of the measured fracture energy as function of ¯ N = 10.0 mm. θN setting θL = 0◦ (hence the test is in the N − TX plane) and H Friction µ is spanned in the range [0, 1.4] and the angle of asperities θa is spanned in the range 0◦ − 30◦ . ¯ N is shown by Figure 7b which contains the polar plots of G ¯T Sensitivity to H ◦ ¯ as a function of θN with θL = 0 and with HN spanned in the range 0.01 − 10.0 mm. Figures 8a and 8b show the interaction diagrams obtained with the same modality and methods adopted for the analyses of Figures 7a and Figure 8a, but with θL = 45◦ , i.e., for a decohesion path which, projected onto the TX − TY plane is the bisector of the first quadrant. Plots in Figure 7a and 8a show that when µ is equal to 0 the polar diagrams collapse to circumferences. This is recognized to be a consequence of the fact (k) that, when frictional dissipation is canceled, frictional stresses σ f are null on all elementary planes as well as their dissipative contribution in (53), so that ¯T = G ¯ c and the work spent for decohesion is independent of the direction of G the impressed displacement. In contrast, when µ 6= 0, the polar curves overlap such circumference until θN > θa . This feature can be interpreted as an evidence of the absence of mechanical interference (and hence of friction phenomena) for ¯ T is an increasing θN > θa . When θN > θa the total measured fracture energy G ¯ function of µ, θa and HN and also increases with the mode II component. A comparison between Figures 7 and 8 also reveals a sensible, altough limited, deviation from isotropy. Such a deviation was already pointed out in [40] where a broader investigation on the anisotropy of the formulation was carried out. Therein, in agreement with [40] anisotropy was found to be an increasing function of µ and θa and could be reduced by using a RME with a higher number of elementary planes. It was shown in [40] that adoption of the 5-planes RME ¯ T , evaluated pure mode II, from shown in Figure 4 results in a deviation for G in-plane isotropy which is lower than 0.2%. 3.3. In-plane isotropy Confined-slip decohesion analyses with σN = 0.0 MPa (i.e., in absence of a confining pressure) and variable θL were performed to investigate the isotropy of the model in plane TX − TY . Interface parameters are again those in Table 2 ¯ N = 0.2 mm. Figure 9a and Figure 9b show the θL -polar graphs of the and H total measured fracture energy and of the maximum dilation s¯N , respectively, for a confined-slip decohesion setting σN = 0. Deviation from isotropy is confirmed to be very limited. Also, Figure 9b shows that the dilation is very close to the height of the asperities, which is an expected result.

25

(a)

(b) µ=0.0, θa=0 a

µ=1.4, θ =10

0.5

a

µ=0.7, θa=20

0.4

µ=1.4, θa=20 µ=0.7, θa=30

0.3

H =0.01

0.6

µ=0.7, θ =10

Mode I Specific Area Energy [kJ/m2]

Mode I Specific Area Energy [kJ/m2]

0.6

µ=1.4, θa=30 0.2

0.1

0

N

H =0.02 N

0.5

H =0.05 N

HN=0.1

0.4

HN=0.2 HN=0.5

0.3

HN=10.0 0.2

0.1

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0

Mode II Specific Area Energy [kJ/m2]

0.1

0.2

0.3

0.4

0.5

0.6

Mode II Specific Area Energy [kJ/m2]

¯ T in a complete displacement-controlled test as a function of θN with Figure 7: Polar graphs of G ¯ N = 10.0 mm, µ spanned in 0 − 1.4, θa θL = 0 using interface parameters in Table 2 and (a): H ¯ N spanned in the interval 0.2 − 2.0 mm spanned in 0◦ − 30◦ , (b) H

(a)

(b) µ=0.0, θa=0

0.6

µ=1.4, θa=10

0.5

µ=0.7, θa=20

0.4

µ=1.4, θa=20 µ=0.7, θa=30

0.3

HN=0.01

0.6

µ=0.7, θa=10

Mode I Specific Area Energy [kJ/m2]

Mode I Specific Area Energy [kJ/m2]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

µ=1.4, θa=30 0.2

0.1

0

HN=0.02 0.5

HN=0.05 HN=0.1

0.4

HN=0.2 HN=0.5

0.3

HN=10.0 0.2

0.1

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0

Mode II Specific Area Energy [kJ/m2]

0.1

0.2

0.3

0.4

0.5

0.6

Mode II Specific Area Energy [kJ/m2]

¯ T in a complete displacement-controlled test as a function of θN Figure 8: Polar graphs of G ¯ N spanned in with θL = 45◦ and (a): µ spanned in 0 − 1.4 and θa spanned in 0◦ − 30◦ , (b) H the interval 0.01 − 10.0 mm.

4. Structural examples: DCB-UBM Simulations To assess the capability of the model to reproduce the experimentally observed decohesion curves and to further investigate at a structural level the mixed mode interaction as predicted by the M-CZM formulation, 3D FEM simulations of decohesion tests with a DCB-UBM were performed using the commercial FEM software Abaqus, [51], implementing the M-CZM formulation via a user-defined 26

(a)

(b) 0.2

0.3

Y

0.15

Y

0.2 0.1

0.1

Dilation [mm]

Specific Area Energy [kJ/m2]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0

X −0.1

0.05 0

X −0.05 −0.1

−0.2

−0.15

−0.3 −0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

−0.2 −0.25 −0.2 −0.15 −0.1 −0.05

0.4

2

Specific Area Energy [kJ/m ]

0

0.05

0.1

0.15

0.2

0.25

Dilation [mm]

¯ T (solid line) and (b) dilation Figure 9: Polar graphs of (a) total measured fracture energy G sn as a function of θL in a confined test slip with σN = 0 confinement pressure (solid line). ¯ N = 0.2. The dotted curve in Figure b is a Interface parameters are those in Table 2 and H perfect circumference with radius equal to the value of the function for θL = 0, plotted as a reference to perfect isotropy.

Fortran subroutine coded with Intel Fortran compiler. 4.1. Simulated test setup The geometry of the structural model employed in the simulations is defined based on the experimental setup described in [12], and reproduces a double cantilever beam DCB with the geometry, constraints and loads shown in Figure 10. The DCB is subjected to quasi-static decohesion tests by slowly increasing the tip-bending moments, Mu , Ml , applied at the free ends of the upper (u) and lower (l) arms of the DCB. The direction and magnitude of the tip bending moments are controlled by three scalar multipliers M1 , M2 and MIII with ˆ the unit vectors of Mu = M1ˆi + MIIIˆj, and Ml = M2ˆi − MIIIˆj, being ˆi, ˆj, k the global reference system shown in Figure 10. The cross section of the arms is rectangular with height H and width B. On the initially unbonded region no interaction is considered, while, in the remaining part of the horizontal mid-plane, interaction between upper and lower arms is modeled with the CZM formulation detailed in Section 2. In all simulations the DCB length is L = 600 mm, the interface thickness is 0.015 mm and the length of the initially unbonded area is Lunb = 60 mm. To examine the sensitivity to the cross-section of the beam arms, three different sections, specified in Table 3, are employed. These sections are identified by the H/B ratios and in particular section H/B = 0.33 corresponds to the set-up actually employed in [12]. 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 10: Sketch of the model used in the simulations of subsection 4.6.1

H/B 0.33 1.00 3.33

B [mm] 30.00 30.00 30.00

H [mm] 9.00 30.00 100.00

Size of elements [mm] 5.00 5.00 10.00

# C3D8 1862 6517 1800

# COH3D8 840 840 210

Table 3: Sections of the beams used in the structural DCB-UBM models, and size and number of elements in the mesh for each section.

4.2. Employed material properties A linear elastic isotropic model is used for the arms of the beams employing the material parameters reported in [12] for the polyester E-glass specimens; these parameters are recalled in Table 4. E [N/mm2 ] 37 · 103

ν [−] 0.30

δ [kg/mm3 ] 2 · 10−6

Table 4: Material parameters used for the arms in the DCB simulations.

The CZM material parameters and the RME employed are selected according to the calibration carried out in [37] based on the experimental data in [12]. The parameters common to all simulations are collected in Table 5, which also reports the inclination angle θa of the inclined elementary planes. The value employed for η was found be sufficiently high so as to not induce in any simulation appreciable variations in the computed initial stiffness. Two RME geometries were employed in the simulations: 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

• the 3D isosceles trapezoidal RME (also examined in reference [40]), shown in Figure 11 with Np = 3, with equipartition of area fractions γ (k) = 1/3, k = 1, . . . , 3 and with inclination angle θm = 40◦ common to the two inclined planes, reproducing in 3D the geometry of the 2D RME of reference [37]; • the 5-plane truncated-pyramid shaped RME with square base of Figure 4 whose geometrical parameters are specified in Table 1, and which is selected again as a computationally inexpensive approximation of in-plane isotropy with a limited number of elementary planes [40], with equipartition of γ (k) employed as a simple, yet reasonable, approximation of an isotropic response.

N Ty

(2)

n

Tx t

(3)

(2)

(1)

h

(2)

h

n

(1)

n t

(3)

h t (3)

(1)

m

m

Figure 11: RME geometry used for the simulations in Section 4.6.1

¯ N , which completes the specification of the The height of the asperities H multi-plane cohesive response, has been object of sensitivity analyses and the employed values are specified next in relation to the particular numerical examples considered. σ0 [MPa] 1.00

¯ c [kJ/m2 ] G 0.97

η [–] 0.95

µ [–] 1.00

θa [◦ ] 40

Table 5: Material parameters employed in the structural DCB-UBM mixed mode I-II simulations.

Simulations consist of full dynamic analyses which are performed reproducing loading conditions generating quasi-static crack growth. For the beam density a value of 2000 kg/m3 is used, which is a rough intermediate between plastics and glass, and no viscosity is considered. No additional damping is considered in the analyses so that dissipation of the model is only originated at the interface as predicted by the M-CZM formulation. The applied loading consists of a slow 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

proportional increase of tip-moment multipliers, M1 , M2 , MIII , at established fixed ratios M1 /M2 , M1 /MIII , and using ramp time-laws. 4.3. Mesh details 3D 8-noded fully integrated solid brick stress-displacement elements (C3D8) are used for the domains of the beam arms while for the interface 3D 8-noded cohesive elements (COH3D8) are used. Structured meshes are employed both for the beam arms and for the interface. The size and number of elements of each mesh are reported in Table 3. 4.4. General features of the computed response The typical response obtained numerically corresponds, as expected, to the classical rate-independent delamination of fracture mechanics, consistently with the rate-independent character of the employed cohesive zone model, and this feature turns out to be not affected by the specific choice for the RME geometry and for the height of asperities. Below specific thresholds of the multipliers, which are found to be only dependent on mode mixity, constitutive parameters and cross section geometry, simulations show absence of delamination; beyond these thresholds, onset of delamination occurs and fracture propagation is observed at constant energy release rates. These thresholds were detected by a simple trial-and-error method. The values of the tip-moments triggering decohesion are found to be largely dependent on the cross section geometry and on mode mixity, and range from a lower bound of approximately 50 kN·m, for pure mode I on the H/B=0.33 section, to an upper bound of 12.8 · 103 kN·m for mixed mode I-II simulation with dominant mode II mixity, M1 /M2 = +0.87, achieved for the H/B = 3.33 section. The general property observed in the previous section is confirmed also by the structural DCB simulations: independently from the employed RME geometry and DCB geometry, when friction is set to zero or when a flat RME is employed (with zero inclination angle for all elementary planes), the plateaus of all computed mixed-mode G-δ curves overlap with the plateau of the mode I plot, recovering the independence of the fracture resistance from the mode mixity, characteristic of basic Linear Elastic Fracture Mechanics (LEFM) [52]. Dynamic oscillations with release of kinetic energy are observed, in particular, once complete delamination is attained. These kinetic effects are found to be increasing with the mode II and mode III components which show a more brittle character. As expected, oscillations increase with the rate of the applied tipmoment multipliers. Kinetic effects were minimized by employing sufficiently slow time laws for moments. In particular, the values of the maximum tip moments attained at the end of the ramps are selected just above the decohesion thresholds, 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

and employing a time range of an order of magnitude of approximately 60 s. In this way conditions of quasi-static crack growth are effectively reproduced. Examples of typical deformed meshes obtained in pure mode III debonding and in mixed mode I-II debonding are reported in Figure 12.

Figure 12: Deformed meshes of a simulation under pure mode III loading with H/B=1.00 (left) and of a simulation under mixed mode I-II loading with H/B=0.33 (right).

Inspection of the process zones resulting from simulations confirm self-similarity of stress, strain and damage fields during crack advancement and that the damaged area contour remains straight and orthogonal to the z-axis in all mixed mode tests, as shown, e.g., by Figures 13 and 14. 4.5. Choices for data postprocessing The above described features permit a meaningful analysis of the interface response as function of mode-mixity in terms of energy release rate G computed as J-integrals [53, 52]. The resulting expressions for G turn out to be particularly simple for the DCB since they are univocal functions of the tip bending moments and of parameters E, B and H. The structural response is examined similar to [12] in terms of G-δ curves with δ being the norm of the relative displacement at the initial crack tip. Almost all computed G -δ curves actually show large plateau regions with fracture propagating at constant G, with few exceptions only close to pure mode II and pure mode III delamination, see Figure 15a. The presence of a plateau is a further direct evidence of self similarity of the process zone during crack progression. In those cases close to pure mode II and pure mode III where the plateau could not be identified, this was found to be the consequence of the large extension of the fracture process zone which reaches the clamped end since the very beginning of the test, thus impeding self-similar process zone advancement characteristic of steady crack propagation in DCBs. The values of the energy release rate G were sampled in the plateau regions and polar plots of G as function of mode mixity angles were constructed. For instance, in Figure 15a the sampled points located at δ = 2.00 mm, are marked by crosses. 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

Figure 13: Examples of typical self-similar contour plots of damage variable α(1) originated over the cohesive zone during the plateau region for a cross section H/B=1.00, employing interface parameters in Table 5. (a): pure mode I decohesion; (b): pure mode III decohesion.

In particular, in mixed mode I-II simulations MIII is set to 0 and M1 and M2 are selected so as to obtain the desired mixed mode ratio M1 /M2 . The related energy release rate GI−II turns out to be:
21 M12 + M22 − 6M1 M2 , (63) GI−II = 4EH 3 B 2 while the mode I-II angle, θI−II , is defined as: θI−II

π = − arctan 4



M1 M2



(64)

With this last definition, in the mode I-II polar plots, pure mode I (M1 /M2 = −1.0) is represented by points on the vertical axes (θI−II = 90◦ ) while pure mode II (M1 /M2 = +1.0) is represented by points on the horizontal axes (θI−II = 0◦ ). For mixed mode I-III simulations, mode II is excluded setting M2 = −M1 and parameters M1 and MIII are selected so as to obtain the desired mixed mode ratio M1 /MIII . The corresponding expression of the energy release rate is: GI−III = 12

2 MI2 MIII + 12 , EB 2 H 3 EHB 4

32

(65)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

Figure 14: Examples of typical self-similar contour plots of damage variable α(1) originated over the cohesive zone during the plateau region in a mixed mode I-II simulation with M1 /M2 = 0.87, employing interface parameters in Table 5 . (a): cross section H/B=0.33 ; (b): cross section H/B=1.00.

and the mode I-III angle, θI−III , is defined as:   M1 θI−III = arctan MIII

(66)

so as to make pure mode I (M1 /MIII = +∞) correspond again to points on the vertical axis of the polar plot, θI−III = 90◦ , and pure mode III to points on the horizontal axis, θI−III = 0◦ . 4.6. Mixed mode I-II analyses Figure 15a shows a family of G-δ plots, with G defined as in (63), generated by spanning M1 /M2 in the interval −1 ≤ M1 /M2 ≤ 0.87 for the H/B=1.00 section employing the 5-plane RME of Figure 4. The polar plot of the plateau values of G as function of θI−II is reported in Figure 15b. As shown in Figure 15a the size of the plateau regions in the G-δ curves reduces as the mode mixity approaches pure mode II. This reduction was found to be strictly related to the enlargement of the process zones in the analyses with increasing mode II component and is examined in more detail in a next paragraph. For M1 /M 2 = 0.87 a plateau cannot be clearly recognized. For this reason the mixed mode I-II simulations were limited to mode ratios not exceeding this value. 33

(a)

(b) 5

Mode I energy release rate [kJ/m2]

4.5 5

Energy release rate [MPa mm]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

4.5

Sim. MI/MII = +0.87

4

Sim. MI/MII = +0.50 Sim. MI/MII = +0.25

3.5

Sim. MI/MII = −0.01

3

Sim. MI/MII = −0.52

2.5

Sim. MI/MII = −1.00 2 1.5 1

60o

4 3.5

45o

3 2.5

30o

2 1.5 1 0.5

0.5

0 0

0

5

10

displacement [mm]

15

20

0

1

2

3

4

5

Mode II energy release rate [kJ/m2]

Figure 15: For the mixed mode I-II model described in section 4.6.2 with the H/B=1.00 section (a) Plot of the Energy release rate as a function of the opening of the crack with the mode mixity spanned in the range -1.0 to +0.87 (b) Polar plot of the plateau value of the energy release rate (crosses in (a)) as a function of the mode mixity angle θI−II .

4.6.1. Mixed mode I-II numerical-experimental comparison A first set of simulations was carried out to compare and validate the 3D model against the results of the 2D FEM model in [37] and against the experimental measurements of the crack growth resistance obtained with a DCB-UBM setup by Sørensen et al. [12]. To suitably compare with the 2D interface formulation, the RME of Figure 11 is employed, which corresponds to the geometry of the 2D RME of reference [37], and is used in combination with the same interface material parameters of this reference, recalled in Table 5. In order to suitably reproduce the feature of unbounded asperities of the CZM formulation in [37] also within the present ¯ N = 20.0 mm which is a value much larger than finite-asperity model, it is set H the range of relative displacements obtained in the analyses at the onset of the plateau regions. In Figure 16a the G-δ plots computed with the present 3D M-CZM formulation, and obtained for the values of mode mixity ratios M1 /M2 ranging in the interval -1 to +0.87, (for which experimental data inn [12] are reported) are plotted by continuous lines. These curves exactly match those in [39] confirming the consistent recovery of the 2D interface formulation by the present one. Figure 16a also reports the experimental measures for fracture toughness obtained by Sørensen et al. [12], which are plotted by markers, showing a reasonable numerical-experimental agreement. 34

(a)

(b) 5 Model

Mode I energy release rate [kJ/m2]

4.5 5 Exp. r = +0.87 Exp. r = +0.50 Exp. r = +0.25 Exp. r = −0.52 Exp. r = −1.00 Sim. r = +0.87 Sim. r = +0.50 Sim. r = +0.25 Sim. r = −0.52 Sim. r = −1.00

4.5 4

J Integral [MPa mm]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3.5 3 2.5 2 1.5 1

60o

4

Exp. Sorensen

3.5

45o

3 2.5

30o

2 1.5 1 0.5

0.5 0

0 0

1

2

3

δ [mm]

4

5

6

0

1

2

3

4

5

Mode II energy release rate [kJ/m2]

Figure 16: For simulations in Section 4.6.1: (a) Comparison between numerical values of G (continuous lines) and experimental values the J-Integral (markers) of as a function of δ with M1 /M2 spanned in the interval -1 to +0.87. (b) Polar plot of the plateau values in figure (a) as a function of the mode mixity angle θI−II .

Figure 17 alternatively shows the numerical results of Figure 16 relevant to pure mode I (M1 /M2 = −1) and to mixed mode M1 /M2 = +0.5 in terms of tipmoment vs. tip rotation diagrams for the upper (U) and lower (L) arms of the DCB. Figure 17 also reports the ideal bilateral curves corresponding to LEFM. The LEFM curves are computed analytically by setting the initial tip momentrotation stiffness to the value EBH 3 /(12Lunb ) corresponding to a Euler-Bernoulli beam having the same length Lunb of the initially unbonded region and with plateau moment set to the value corresponding to pure mode I (i.e., obtained by setting in (63) M1 = −M2 and setting the energy release rate to the value of the ¯ c of Table 5 employed in the simulations). rupture energy G The comparison shows that the numerical pure mode I curves essentially match with the LEFM curves, although the initial elastic slope of the numerical model is obviously lower than the one provided by the Euler-Bernoulli model since the latter considers a perfect clamp at the initial crack tip. A marked transition from the linear elastic regime to steady crack propagation regime can be observed also in Figure 17, as well as the onset of deviations from the plateau region at relatively large rotations. As previously observed, these deviations were found to be associated with the process zone reaching the clamped end of the DCB. Comparison of Figures 17 and 16a also shows that G-δ plots permit an easier reading of the output data for mixed-mode responses. 35

5

2

x 10

Ideal Euler−Bernoulli DCB LEFM M1/M2=−1 (L)

1.5

Ideal Euler−Bernoulli DCB LEFM M1/M2=−1 (U) M−CZM simulation M1/M2=+0.5 (L)

Tip Moment [N mm]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1

M−CZM simulation M1/M2=+0.5 (U) M−CZM simulation M1/M2=−1 (L) M−CZM simulation M1/M2=−1 (U)

0.5 0 −0.5 −1

−0.5

0

0.5

1

Tip Rotation [−]

Figure 17: Tip-moment vs. tip rotation diagrams for the upper (U) and lower (L) arms of the DCB corresponding to pure mode I (M1 /M2 = −1) and to mixed mode M1 /M2 = +0.5. The figure also shows the ideal bilateral curves corresponding to LEFM. See legend for details.

4.6.2. Mixed mode I-II sensitivity analyses A set of mixed mode I-II simulations were performed to analyze the sensitivity of the formulation to model parameters and to the DCB geometry. For these simulations the RME geometry shown in Figure 4 and the model parameters in ¯ N = 2.00 mm. Table 5 are used again, setting H Figure 15a shows the family of G-δ plots obtained by spanning M1 /M2 in the range -1.0 to 0.87 for the section H/B=1.00 with the corresponding polar plot of the plateau values of G shown in Figure 15b. Contour plots of damage in the cohesive horizontal plane were post-processed from the output numerical data to inspect the fracture process zone in deformed states corresponding to the plateau regions of Figure 15a. The process zone is identified as the region of this plane having partial damage state, 0 < α(k) < 1 for some k. In particular, Figure 13 shows samples of the self-similar contour plots of damage to the first elementary plane, α(1) , of the RME (which is oriented as the horizontal X-Z global plane, see Table 1 and Figure 2) for the arm section H/B=1.00, subjected to pure mode I decohesion. The blue regions correspond to absence of damage in this plane (α(1) = 0), and the red ones to full damage (α(1) = 1). Figure 14b shows the corresponding plot obtained for mode mixity +0.87. An inspection of the amplitude of the process zones in these two figures highlights the increase in their extent as mode II becomes dominant. Figure 14a shows the process zone determined under the same conditions of Figures 13a and 14b, upon setting M1 /M2 = +0.87 and employing the cross section H/B = 0.33. An amplitude which is intermediate between the ones of 36

Figures 13a and 14b is observed which shows that the length of the process zone is also affected by the beam stiffness. Further simulations were performed to investigate the influence of the height of the asperities and of the H/B ratio on the computed interface toughness. Figure 18a shows a polar plot of G as a function of M1 /M2 for three H/B ratios in Table 3. Figure 18b shows the family of polar plots of G for the H/B=1.0 section obtained ¯ N in the range 0.02 mm-2000 mm. These plots show a pronounced by spanning H effect on fracture toughness of both the height of the asperities as well as of the beam geometry, which increases in proximity of pure mode II mixity. Conversely, ¯ N and H/B on fracture toughness is observed to decrease as the influence of H mode I mixity is approached, and completely vanishes in pure mode I where the ¯ c is invariably recovered. coincidence G ≡ G (a)

(b)

5

5 H/B=0.33

(a)

H/B=1.00

60o

H/B=3.33

4 3.5

45o

3 2.5

o

30

2 1.5 1 0.5 0

HN = 2000 mm HN = 2.0 mm

4.5

Mode I energy release rate [kJ/m2]

4.5

Mode I energy release rate [kJ/m2]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

HN = 0.2 mm

4

60o

H = 0.02 mm N

3.5

45o

3 2.5

o

30

2 1.5 1 0.5

0

1

2

3

4

5 2

Mode II energy release rate [kJ/m ]

0

0

1

2

3

4

5

Mode II energy release rate [kJ/m2]

Figure 18: Polar plot of the plateau value of GI−II as function of θI−II , for the mixed mode I-II ¯ N = 2.00 mm (b) with fixed delaminations of Section 4.6.2. (a) with a variable section and H ¯ section H/B=0.33 and HN spanned in 0.2 mm-2000 mm

4.7. Mixed mode I-III simulations Results of mixed mode I-III simulations are also reported in terms of families of G-δ plots, with G now computed by (65). The 5-plane RME of Figure 2 and the model parameters of Table 5 are again employed. Figure 19 shows the family of G-δ plots obtained by spanning M1 /MIII in the interval 0.0-10.0, for section H/B=1.00. Figure 20 shows the corresponding polar plots of the plateau values of G, as a function of the mode mixity angle θI−III . The sensitivity of mixed mode I-III toughness to the cross section geometry is examined in Figure 21 which shows a polar plot of the plateau value of G for the 37

sections listed in Table 3. The polar plot shows that the plateau value of G from θI−III = 90◦ downwards is nearly constant until 40◦ ≤ θI−III ≤ 90◦ . For θI−III < 40◦ an abrupt increase of the energy release rate can be observed which can be explained as an effect of the interplay between friction and interlocking. Figure 21 also shows that different H/B ratios result in a very pronounced difference in the polar diagrams of the energy release rate. (a)

(b)

25

2

Energy release rate [MPa mm]

1.8

Energy release rate [MPa mm]

20

15

Sim. MI/MIII = +0.000 Sim. MI/MIII = +0.125 Sim. MI/MIII = +0.250

10

Sim. MI/MIII = +0.500 Sim. MI/MIII = +0.840 Exp. MI/MIII = +1.000

5

Exp. MI/MIII = +2.000 0

0

10

20

30

δ [mm]

40

50

1.6 1.4 1.2

Sim. MI/MIII = +0.000

1

Sim. MI/MIII = +0.125

0.8

Sim. MI/MIII = +0.250 Sim. MI/MIII = +0.500

0.6

Sim. MI/MIII = +0.840

0.4

Exp. MI/MIII = +1.000 Exp. MI/MIII = +2.000

0.2

Exp. MI/MIII = +10.00 60

0

Exp. MI/MIII = +10.00 0

1

2

3

δ [mm]

4

5

6

Figure 19: For the mixed mode I-III delaminations described in section 4.7 (a) GI−III as a function of the opening of the crack with the mode III to mode I ratio spanned in the range 0 − 10.0 (b) Magnification of the diagrams contained in the dashed box in the figure (a)

(a)

(b)

25

5 4.5

60o

20

Mode I energy release rate [kJ/m2]

Mode I energy release rate [kJ/m2]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

45o

15

30o

10

5

60o

4 3.5

45o

3 2.5

30o

2 1.5 1 0.5

0

0

5

10

15

20

25

Mode III energy release rate [kJ/m2]

0

0

1

2

3

4

5

Mode III energy release rate [kJ/m2]

Figure 20: For the mixed mode I-III delaminations described in section 4.7 (a) Polar plot of the plateau value of GI−III as function of the mode mixity angle θI−III with section H/B=1.00 and ¯ N = 2.00 mm (b) Magnification of the diagrams contained in the dashed box in figure (a) H

38

(a)

(b) 5

25

H/B=0.33

H/B=0.33 H/B=1.00 H/B=3.33

20

4.5

Mode I energy release rate [kJ/m2]

Mode I energy release rate [kJ/m2]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

o

60

45o

15

o

30

10

5

H/B=1.00

60o

H/B=3.33

4 3.5

45o

3 2.5

30o

2 1.5 1 0.5

0

0

5

10

15

20

Mode II energy release rate [kJ/m2]

25

0

0

1

2

3

4

5

Mode II energy release rate [kJ/m2]

Figure 21: For the mixed mode I-III delaminations described in section 4.7 (a) Polar plot of the plateau value of GI−III as function of the mode mixity angle θI−III with the sections in table 3 ¯ N = 2.00 mm (b) Magnification of the diagrams contained in the dashed box in figure (a) and H

5. Discussion and Conclusions A 3D multi-scale multi-plane cohesive-zone model has been developed for the FEM analysis of structural interfaces. The quasi-isotropic character of the local response of the interface has been assessed in terms of fracture resistance and dilation. The proposed 3D interface model has been also employed to simulate double-cantilever beam tests and to determine predictions of the mixed mode I-II and mixed mode I-III fracture resistance. The study has highlighted theoretical, numerical and predictive advantages of the M-CZM approach, which are hereby summarized. From the theoretical point of view, the M-CZM formulation is framed within the thermo-mechanics of generalized continua, and employs a reduced set of micromechanical parameters characterized by a well-defined physical meaning which permits to devise clear calibration and identification procedures also for 3D fracture problems. Moreover, increase of apparent total measured fracture energy from mode I to mode II naturally emerges as the effect of multi-scale coupling between cohesion, friction and interlocking. From the numerical viewpoint, a robust return-mapping algorithm is developed, in the wake of well-established procedures for computational inelasticity, taking also advantage of the derivation of an explicit formula for the tangent operator which ensures quadratic convergence in Newton-Raphson schemes. The reported results show that the adopted integration procedure of the M-CZM formulation permits, by a proper selection of the RME, an efficient multiscale analysis with a very low computational effort. 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

DCB simulations of mixed mode I-II and mixed mode I-III decohesion tests show that, by properly calibrating few material parameters, the 3D model can reproduce, in a wide range of mode mixity, the experimentally observed mixedmode I-II response in polyester E-glass laminates. The DCB simulations also provide a wide database of numerical predictions of mixed mode I-II and mixed mode I-III fracture resistance, open to future experimental quantitative validation, which confirms known trends of fracture mechanics, and allows to identify open lines of research for the fine-tuning of testing standards with the aid of numerical procedures. In particular these analyses have highlighted that: a) the size of the process zone of the interface grows larger as the beam becomes stiffer and with increasing mode II and mode III components; b) the height of the asperities plays an important role in determining the energy release rate both in mixed mode I-II and mixed mode I-III interaction; c) different DCB sections result in significantly different energy release rate in mixed mode I-mode III interaction, to the extent that it is not possible to analyze mode III neglecting the geometry of the beam. The shown dependence of the macroscopically measured mode-II and modeIII energies on the geometry of the DCB cross section points out the capability of the multiplane approach to identify, with a computational effort considerably lower than full multiscale approaches, an important factor of dependence of the crack growth resistances upon the test-setup geometry. Such a dependence appears to be relevant in all those debonding and delamination phenomena where friction and mesoscale irregularities are significant in the fracture process. It is worth being observed that these results seem to indicate that any macroscopic model postulating the existence of a predetermined mode II or mode III fracture energy without accounting for friction should not predict such dependence (in absence of heuristic corrections accounting for the geometry of the test setup). As a final consideration, it is worth to remark that even if the investigation has been limited herein to a rate-independent quasi brittle response and to the recovery of isotropy, the versatility of the 3D multi-plane description makes the proposed M-CZM approach prone to the description of the behavior of anisotropic interfaces (e.g., such as those considered in [54, 55, 56]), as well as of rate-dependent behaviors, by properly selecting the specific 3D RME geometry and the constitutive response of the elementary planes. Appendix A. Numerical integration In this appendix finite-step relationships for the integration of the state variables according to the formulation of Section 2 in a generic timestep [τ n , τ ] are 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

developed for numerical implementation purposes. The notation convention is used of denoting variables at the beginning of the timestep by the addition of a (·)n superscript and the value at the end of the timestep without any additional superscript. Section Appendix A.2 reports the derivation of the expression of the consistent tangent operator which allows the efficient implementation of the M-CZM formulation in nonlinear gradient-based integration algorithms, such as in Newton-Raphson schemes. Appendix A.1. Finite step relationships Finite step relationships for the integration of the mechanical state of plane (k) are obtained by replacing equations (55), (57), (40), (38) and the KuhnTucker conditions (46) with their discretized forms and by integrating equation (45) in the timestep [τ n , τ ] via a backward Euler scheme [50] and using a linear (k) time interpolation for s(k) and sdi :

σ (k) = (1 − α(k) )Ks(k) + α(k) Ag (k)

σf



(k)

= H(k) s(k) − sdi ˆ (k) s(k) k kP n (+)

β

(k)

α(k)

=

s0 (



s¯N ¯ (k) H N

!

(k)

(A.1)

σf

(A.2)

−1

(A.3) (

1 β (k) = max α(k),n , min 1, · η 1 + β (k)

))

(A.4)

(k)

(k) sdi

=

(k),n sdi

+ ∆λ

τf

(k)

(A.5)

(k)

k τf k

∆λ(k) ≥ 0

(A.6)

(k)

(k)

µhσf n i + k τ f k≤ 0

(A.7)

(k)

(A.8)

−

(k)

∆λ(k) (µhσf n i + k τ f k) = 0 −

(k)

(k)

Equations (A.1)-(A.8) allow to compute the value of σ f , sdi , α(k) at the end of the timestep given the values of the same parameters at the beginning of the timestep and the displacement s(k),n , s(k) at the beginning and at the end of the timestep. 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The solution scheme is briefly summarized. As a first step, the updated values of β (k) and α(k) are computed via (A.3) and (A.4). Solution of Equations (A.2)(A.8) is achieved by an elastic prediction and inelastic correction scheme [50]. Accordingly, an elastic stress trial value is computed first as:   (k) (k),n σ f e = H(k) s(k) − sdi . (A.9)   (k) If the corresponding trial value of the yield function φ σ f e is negative, then

equation (A.8) is trivially solved by ∆λ(k) = 0 and, consequently:  (k) (k),n  sdi = sdi    (k) (k) σ f = H(k) (s(k) − sdi )     s¯N (k) (k) (k) (k)   σ = (1 − α )Ks + α Ag ¯ (k) σ f (k) Conversely, if φ σ f e

(A.10)

HN


(k)

is positive, consistently with the adoption of a full im˙ (k) ) = 0 is plicit backward Euler scheme ∆λ(k) , the consistency condition φ(σ f enforced at the end of the time step, and, stated explicitly, reads (k) φ˙ (k) (σ f ) =

=

∂φ(k) (k)

∂σ f

∂φ(k) (k) ∂σ f

(k)

· σ˙ f = 

 ∂φ(k)  (k)  (k) (k) ˙ ˙ s − s = · H di (k) ∂σ f 



(A.11)

∂g · H(k) s˙ (k) − λ˙ (k) (k)  = 0. ∂σ f

Owing to the specific choices for the yield function and plastic potential, represented by (42) and (43), the vectors in (A.11) turn out to be: ∂g (k) (k)

∂σ f

(k)

=

τf

(k)

k τf k

,

∂φ(k) (k)

∂σ f

(k)

= −u(−σf n )µn

(k)

+

τf

(k)

k τf k

.

(A.12)

Due to linear time interpolation for state variables, one has (k)

(k) σf

so that

=

(k) σf e

− ∆λ

(k)

τf

(k)

k τf k

(k)

K

τf

(k)

k τf k

,

(A.13)

(k)

= 42

τ fe

(k)

k τ fe k

,

(A.14)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

and from (A.11) one computes  −1   (k) ∂φ ∂φ(k) ∂g (k) (k) (k),n ∆λ(k) =  (k) · H(k) (k)  · H s − s , (k) ∂σ f ∂σ f ∂σ f

(A.15)

Due to (A.14), the partial derivatives in (A.15) can be computed at the trial stress (k) values σ f e , and, according to the previous expressions, ∆λ and the updated frictional stresses are computed as follows:   1  (k) (k) (k)  µσ + k τ k ∆λ =  f en fe   K       (k) (k) σf n = σf en (A.16)     (k)   τ fe  (k) (k)    τ f = −µσf en (k) k τ fe k

Once ∆λ(k) is determined, Equations (A.1), (A.2) and (A.5) are solved for sdi (k) , σ d (k) and σ (k) and, finally, upon integrating the stress of each elementary plane, the overall stress σ is reassembled via (54). Appendix A.2. Tangent operator The algorithmic tangent operators are derived consistently with the finite step relationships introduced in section Appendix A.1. To this end the consistent (k) tangent operator of the RME and of the elementary planes (k), KT and KT , respectively defined as: KT =

∂σ ∂¯s

(k)

KT =

,

∂σ (k) ∂s(k)

(A.17)

and related via (10) and (54) by KT =

Np X

(k)

γ˜ (k) KT

(A.18)

k=1

can be evaluated. Differentiation of equation (55) yields:   ∂α(k) s¯N ∂α(k) ∂σ (k) (k) (k) (k) (k) = −K s ⊗ + (1 − α )K + Ag ¯ σf ⊗ + ∂¯s ∂¯s ∂¯s HN α

(k)

(k) σf

      (k) ∂ s¯N s¯N ⊗ Ag ¯ + Ag ¯ α(k) Cep . f ∂¯s HN HN 43

(A.19)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The term ∂α(k) /∂¯s in (A.19) is computed as follows:  ! (+) ¯ ˆ  s s P s 1 0 ∂α(k)  0 1− if α(k),n < <1 ˆ (+)¯s k3 ˆ (+) s(k) k = η kP η k P  ∂¯s  0 otherwise while ∂Ag /∂¯s is given by

being ∂Ag (x) = ∂x

(A.20)

∂Ag 1 ∂Ag = ¯N N ∂¯s ∂x H

(A.21)



(A.22)

−1 0 1

(k)  in (A.19) is the consistent algorithmic elastoplastic tangent The term Cep f operator [50] associated with the frictional behavior of plane (k), and is implicitly defined as the tensor satisying  (k) (k) ep ∆σ f = Cf ∆s(k) , (A.23)

in agreement with the discretization strategy outlined in Section Appendix A.1 for the updating of friction variables. Accordingly, when ∆λ(k) is equal to zero (k) one trivially has ∆σ f = H(k) ∆s(k) and the tangent operator coincides with the elastic one. Otherwise, relations (A.5) and (A.15) yield  −1   (k) (k) (k) ∂φ ∂g ∂g ∂φ  (k) (k) ∆sdi (k) =  (k) · H(k) (k)   (k) ⊗ H ∆s . (A.24) (k) ∂σ ∂σ ∂σ ∂σ f

f

f

The previous equation combined with (A.2) yields the following final expression of the algorithmic tangent operator, summarized below for unloading and loading conditions:  (k) ∆λ = 0  H  −1     (k)  (k) Cep = ∂φ(k)  (k) ∂φ(k) f (k) ∂g  (k)  ∂g (k)   H H ∆λ < 0 H − · H ⊗  (k) (k) (k) (k)  ∂σ ∂σ ∂σ ∂σ f

f

f

f

(A.25)

Appendix B. Acknowledgements This study was partially funded by Regione Campania through European funds (POR Campania FSE 2007-2013) for the Dottorato in azienda project (CUP F82I11001150002). 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

References [1] M. Ortiz, A. Pandolfi, Finite-deformation irreversible cohesive elements for threedimensional crack-propagation analysis, International Journal for Numerical Methods in Engineering 44 (9) (1999) 1267–1282. [2] A. Pandolfi, P. Guduru, M. Ortiz, A. Rosakis, Three dimensional cohesive-element analysis and experiments of dynamic fracture in c300 steel, International Journal of Solids and Structures 37 (27) (2000) 3733–3760. [3] J. Segurado, J. LLorca, A new three-dimensional interface finite element to simulate fracture in composites, International Journal of Solids and Structures 41 (11) (2004) 2977–2993. [4] C. Chen, O. Kolednik, J. Heerens, F. Fischer, Three-dimensional modeling of ductile crack growth: Cohesive zone parameters and crack tip triaxiality, Engineering Fracture Mechanics 72 (13) (2005) 2072–2094. [5] N. Sukumar, D. Chopp, E. B´echet, N. Mo¨es, Three-dimensional non-planar crack growth by a coupled extended finite element and fast marching method, International Journal for Numerical Methods in Engineering 76 (5) (2008) 727. [6] L. Macorini, B. A. Izzuddin, A non-linear interface element for 3d mesoscale analysis of brick-masonry structures, International Journal for numerical methods in Engineering 85 (12) (2011) 1584–1608. [7] R. R. Van Der Pluijm, Out-of-plane bending of masonry: behaviour and strength, Ph.D. thesis, Technische Universiteit Eindhoven (1999). [8] J. R. B. Popehn, A. E. Schultz, M. Lu, H. K. Stolarski, N. J. Ojard, Influence of transverse loading on the stability of slender unreinforced masonry walls, Engineering Structures 30 (10) (2008) 2830–2839. [9] A. G. Evans, J. W. Hutchinson, Effects of non-planarity on the mixed mode fracture resistance of bimaterial interfaces, Acta Metallurgica 37 (3) (1989) 909–916. [10] A. Brunner, P. Fl¨ ueler, Prospects in fracture mechanics of engineering laminates, Engineering fracture mechanics 72 (6) (2005) 899–908. [11] G. Allegri, M. Jones, M. Wisnom, S. Hallett, A new semi-empirical model for stress ratio effect on mode ii fatigue delamination growth, Composites Part A: Applied Science and Manufacturing 42 (7) (2011) 733–740. [12] B. F. Sørensen, T. K. Jacobsen, Characterizing delamination of fibre composites by mixed mode cohesive laws, Composites Science and Technology 69 (3) (2009) 445–456. [13] A. Brunner, B. Blackman, P. Davies, A status report on delamination resistance testing of polymer–matrix composites, Engineering Fracture Mechanics 75 (9) (2008) 2779–2794. [14] M. Raous, L. Cang´emi, M. Cocu, A consistent model coupling adhesion, friction, and unilateral contact, Computer Methods in Applied Mechanics and Engineering 177 (3-4) (1999) 383 – 399. doi:http://dx.doi.org/10.1016/S0045-7825(98)00389-2. URL http://www.sciencedirect.com/science/article/pii/S0045782598003892 [15] G. Alfano, E. Sacco, Combining interface damage and friction in a cohesive-zone model, International Journal for Numerical Methods in Engineering 68 (5) (2006) 542–582. doi:10.1002/nme.1728. URL http://dx.doi.org/10.1002/nme.1728 [16] G. Alfano, S. Marfia, E. Sacco, A cohesive damage-friction interface model accounting for water pressure on crack propagation, Computer Methods in Applied Mechanics and Engineering 196 (1-3) (2006) 192 – 209. doi:http://dx.doi.org/10.1016/j.cma.2006.03.001. URL http://www.sciencedirect.com/science/article/pii/S004578250600096X [17] E. Sacco, F. Lebon, A damage–friction interface model derived from micromechanical approach, International Journal of Solids and Structures 49 (26) (2012) 3666–3680. [18] G. Del Piero, M. Raous, A unified model for adhesive interfaces with damage, viscosity, and friction, European Journal of Mechanics-A/Solids 29 (4) (2010) 496–507.

45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[19] P. Carrara, L. De Lorenzis, A coupled damage-plasticity model for the cyclic behavior of shear-loaded interfaces, Journal of the Mechanics and Physics of Solids 85 (2015) 33–53. [20] D. Aubry, A. Modaressi, H. Modaressi, A constitutive model for cyclic behaviour of interfaces with variable dilatancy, Computers and Geotechnics 9 (1-2) (1990) 47 – 58, special Issue on Soil-Structure Interaction. doi:http://dx.doi.org/10.1016/0266-352X(90)90028-T. URL http://www.sciencedirect.com/science/article/pii/0266352X9090028T [21] M. F. Snyman, J. B. Martin, A consistent formulation of a dilatant interface element, International Journal for Numerical and Analytical Methods in Geomechanics 16 (7) (1992) 493–527. doi:10.1002/nag.1610160703. URL http://dx.doi.org/10.1002/nag.1610160703 [22] H. Parland, Stability of rigid body assemblages with dilatant interfacial contact sliding, International Journal of Solids and Structures 32 (2) (1995) 203 – 234. doi:http://dx.doi.org/10.1016/0020-7683(94)00100-B. URL http://www.sciencedirect.com/science/article/pii/002076839400100B [23] Z. Mroz, G. Giambanco, An interface model for analysis of deformation behaviour of discontinuities, International Journal for Numerical and Analytical Methods in Geomechanics 20 (1) (1996) 1–33. [24] S. Stupkiewicz, Fiber sliding model accounting for interfacial micro-dilatancy, Mechanics of Materials 22 (1) (1996) 65 – 84. doi:http://dx.doi.org/10.1016/0167-6636(95)00026-7. URL http://www.sciencedirect.com/science/article/pii/0167663695000267 [25] S. Stupkiewicz, Z. Mrz, Modelling of friction and dilatancy effects at brittle interfaces for monotonic and cyclic loading, Journal of Theoretical and Applied Mechanics 39 (3). URL http://www.ptmts.org.pl/jtam/index.php/jtam/article/view/v39n3p707 [26] A. Spada, G. Giambanco, P. Rizzo, Damage and plasticity at the interfaces in composite materials and structures, Computer Methods in Applied Mechanics and Engineering 198 (49-52) (2009) 3884 – 3901. doi:http://dx.doi.org/10.1016/j.cma.2009.08.024. URL http://www.sciencedirect.com/science/article/pii/S0045782509002692 [27] I. Mihai, A. Jefferson, A multi-asperity plastic-contact crack plane model for geomaterials, International Journal for Numerical and Analytical Methods in Geomechanics 37 (11) (2013) 1492–1509. [28] L. De Lorenzis, Some recent results and open issues on interface modeling in civil engineering structures, Materials and structures 45 (4) (2012) 477–503. [29] P. Wriggers, S. Moftah, Mesoscale models for concrete: Homogenisation and damage behaviour, Finite elements in analysis and design 42 (7) (2006) 623–636. [30] V. Palmieri, L. De Lorenzis, Multiscale modeling of concrete and of the frp–concrete interface, Engineering Fracture Mechanics 131 (2014) 150–175. [31] K. Matouˇs, M. G. Kulkarni, P. H. Geubelle, Multiscale cohesive failure modeling of heterogeneous adhesives, Journal of the Mechanics and Physics of Solids 56 (4) (2008) 1511–1533. [32] M. C. Alfaro, A. Suiker, C. Verhoosel, R. De Borst, Numerical homogenization of cracking processes in thin fibre-epoxy layers, European Journal of Mechanics-A/Solids 29 (2) (2010) 119–131. [33] C. V. Verhoosel, J. J. Remmers, M. A. Guti´errez, R. De Borst, Computational homogenization for adhesive and cohesive failure in quasi-brittle solids, International Journal for Numerical Methods in Engineering 83 (8-9) (2010) 1155–1179. [34] G. Pijaudier-Cabot, J. Mazars, Damage models for concrete, Handbook of materials behavior models 2 (2001) 500–512. [35] A. Caballero, C. L´ opez, I. Carol, 3d meso-structural analysis of concrete specimens under uniaxial tension, Computer Methods in Applied Mechanics and Engineering 195 (52) (2006) 7182–7195. [36] R. Serpieri, G. Alfano, Bond-slip analysis via a thermodynamically consistent interface model combining interlocking, damage and friction, International Journal for Numerical

46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[37]

[38]

[39]

[40]

[41]

[42] [43]

[44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56]

Methods in Engineering 85 (2) (2011) 164–186. doi:10.1002/nme.2961. URL http://dx.doi.org/10.1002/nme.2961 R. Serpieri, E. Sacco, G. Alfano, A thermodynamically consistent derivation of a frictional-damage cohesive-zone model with different mode i and mode ii fracture energies, European Journal of Mechanics - A/Solids 49 (2015) 13 – 25. doi:http://dx.doi.org/10.1016/j.euromechsol.2014.06.006. URL http://www.sciencedirect.com/science/article/pii/S0997753814000837 R. Serpieri, L. Varricchio, E. Sacco, G. Alfano, Bond-slip analysis via a cohesive-zone model simulating damage, friction and interlocking, Fracture and Structural Integrity (29) (2014) 284–292. R. Serpieri, G. Alfano, E. Sacco, A mixed-mode cohesive-zone model accounting for finite dilation and asperity degradation, International Journal of Solids and Structures (67–68) (2015) 102–115. M. Albarella, R. Serpieri, G. Alfano, E. Sacco, A 3d multiscale cohesive zone model for quasi-brittle materials accounting for friction, damage and interlocking, European Journal of Computational Mechanics 24 (4) (2015) 144–170. H. S. Lee, Y. J. Park, T. F. Cho, K. You, Influence of asperity degradation on the mechanical behavior of rough rock joints under cyclic shear loading, International Journal of Rock Mechanics and Mining Sciences 38 (7) (2001) 967–980. C. L.-L. O. H. Shima, Hiroshi, Bond characteristics in post-yield range of deformed bars, Concrete Library of jsce 10 (1987) 113–124. F. Parrinello, G. Marannano, G. Borino, A thermodynamically consistent cohesivefrictional interface model for mixed mode delamination, Engineering Fracture Mechanics 153 (2016) 61–79. J. Lemaitre, J.-L. Chaboche, Mechanics of solid materials, Cambridge University press, 1994. L. Kachanov, Introduction to continuum damage mechanics, Vol. 10, Springer Science & Business Media, 2013. P. Germain, Cours de m´ecanique des milieux continus, Vol. 1, Masson, 1973. R. T. Rockafellar, Convex analysis, Princeton University Press, 2015. W. Karush, Minima of functions of several variables with inequalities as side constraints, Ph.D. thesis, Master’s thesis, Dept. of Mathematics, Univ. of Chicago (1939). H. W. Kuhn, A. W. Tucker, Proceedings of 2nd berkeley symposium (1951). J. C. Simo, T. J. R. Hughes, Computational inelasticity, Springer, 2008. D. Simulia, Abaqus 6.14 analysis user’s manual, Abaqus 6.14 (2014) 22–2. Z. P. Bazant, J. Planas, Fracture and size effect in concrete and other quasibrittle materials, Vol. 16, CRC press, 1997. J. R. Rice, A path independent integral and the approximate analysis of strain concentration by notches and cracks, Journal of Applied Mechanics 35 (2) (1968) 379–386. A. Konyukhov, P. Vielsack, K. Schweizerhof, On coupled models of anisotropic contact surfaces and their experimental validation, Wear 264 (7) (2008) 579–588. A. Zmitrowicz, Models of kinematics dependent anisotropic and heterogeneous friction, International Journal of Solids and Structures 43 (14) (2006) 4407–4451. M. Paggi, J. Reinoso, An anisotropic large displacement cohesive zone model for fibrillar and crazing interfaces, International Journal of Solids and Structures 69 (2015) 106–120.

47