A nonlinear finite thickness cohesive interface element for modeling delamination in fibre-reinforced composite laminates

Accepted Manuscript A nonlinear finite thickness cohesive interface element for modeling delamination in fibre-reinforced composite laminates J. Reinoso, M. Paggi, A. Blázquez PII:

S1359-8368(16)31256-2

DOI:

10.1016/j.compositesb.2016.10.042

Reference:

JCOMB 4640

To appear in:

Composites Part B

Received Date: 8 July 2016 Revised Date:

5 October 2016

Accepted Date: 16 October 2016

Please cite this article as: Reinoso J, Paggi M, Blázquez A, A nonlinear finite thickness cohesive interface element for modeling delamination in fibre-reinforced composite laminates, Composites Part B (2016), doi: 10.1016/j.compositesb.2016.10.042. 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.

ACCEPTED MANUSCRIPT

A nonlinear finite thickness cohesive interface element for modeling delamination in fibre-reinforced composite laminates J. Reinosoa,∗, M. Paggib , A. Bl´azqueza and Strength of Materials Group, School of Engineering, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092, Seville, Spain b IMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100 Lucca, Italy

RI PT

a Elasticity

Abstract

M AN U

SC

Delamination events are major issues which notably affect the integrity of composite structures. To minimize the experimental efforts, there is an increasing demand for developing reliable numerical tools that can accurately simulate delamination initiation and propagation under mixed-mode loading conditions. The current investigation is concerned with the formulation and the finite element (FE) implementation of a new nonlinear finite thickness cohesive interface model for delamination analysis of fibre-reinforced composite laminates relying on the solid shell concept. The incorporation of geometrically nonlinear effects into the proposed interface formulation is motivated by the recent trend of producing composite structures that can experience large displacements prior to failure, as is the case of postbuckling in stiffened panels. The inelastic material behavior of the interface is modeled using two standard nonlinear decohesion laws: (i) an exponential-based, and (ii) a polynomial-based interface laws. Finally, the performance of the proposed interface element is demonstrated by means of several examples focusing on double cantilever beam (DCB) and rib-stiffened specimens. A excellent level of accuracy is achieved when comparing the numerical predictions and the available experimental data.

1. Introduction

TE

D

Keywords: A. Finite thickness interface. B. Cohesive zone modelling. C. Delamination. D. Finite Element Method. E. Fracture

AC C

EP

The setting up of advanced numerical models for the simulation of decohesion and/or delamination events has been an area of extensive research over the last two decades. These failure mechanisms play a crucial role in the structural integrity in many engineering applications involving layered composites, especially in aeronautical and aerospace components and more recently in automotive parts, among many others.. From the numerical perspective, these inelastic processes have been widely modeled using two popular FE-based strategies: (i) the so-called virtual crack closure technique (VCCT) [1], and (ii) Cohesive Zone Models (CZMs) relying on the corresponding interface model [2, 3]. In massive simulations, CZMs have been successfully adopted due to their high versatility and relative simplicity to be incorporated into research and commercial FE codes, especially in applications where the crack path can be clearly identified a priori. In contrast to linear elastic fracture mechanics (LEFM) strategies, CZMs enable the reproduction of the fracture process zone (FPZ) ahead the crack tip, wherein the nonlinear interfacial softening response plays an important role. A general interpretation of CZMs can be conceived by referring the inelastic traction-separation law (TSL), which characterizes the degradation response, to a reference surface. This reference surface is typically identified as the middle surface between the two adjacent bodies (bulks) joined through the common interface. ∗ Corresponding

authors Email address: [email protected] (J. Reinoso)

Preprint submitted to Composites Part B

October 17, 2016

ACCEPTED MANUSCRIPT

AC C

EP

TE

D

M AN U

SC

RI PT

Finite elements based on this methodology are usually denominated as zero-thickness or surface-like CZ elements. In the related literature, a great number of CZMs using this zero-thickness (surface-like) element topology have been proposed in conjunction with different descriptions of the TSL [4]. In this regard, the vast majority of the proposed CZMs are based on phenomenological arguments or simple mathematical relations to describe the TSL that governs the interface such as polynomial [5], bi-linear [6, 7], exponential [8, 9, 10] or following continuum-based formulations [11, 12], to quote a few of them. A possible categorization of these cohesive relations regards potential-derived and non-potential derived inelastic laws [13, 14, 15]. An attempt to define the CZM relation from damage evolution taking place within a finite thickness interface has been performed in [16, 17, 18]. An alternative interface approach to characterize the response is the so-called spring type interface, which can be idealized as a continuous distribution of linear-elastic springs with the appropriate mechanical parameters in terms of stiffness and cut-off strength [19, 20]. This interface conception is usually denominated as linear elastic brittle interface model (LEBIM) [21], which is also equipped with the corresponding strength and/or fracture criterion to account for the interface failure under mixed mode conditions [21, 22]. One of the principal ingredients of LEBIM regards its simplicity and robustness, leading to very accurate predictions within the context of the Boundary Element Method (BEM) [23] and FEM [24]. More recently, further contributions with regard to CZMs have dealt with the incorporation of geometrically nonlinear effects in the description of the interface failure, see [25, 26] and the references therein given. The thermodynamic foundations of such interface models have been thoroughly revisited in [27], where a novel interface model based on Helmholtz free-energy definition complying with the Clausius-Duhem inequality is advocated. However, from the operative point of view in numerical simulations, one of the most critical limitations of interface elements regards the determination of their location in the corresponding discretization. In this concern, according to the clasification established by Paulino and coauthors (see [13] and the references therein given), CZMs can be clasified in terms of the way they are inserted in the computations. Thus, on the one hand, intrinsic cohesive models correspond to those which are introduced between the bulk with elements before running the computations [6, 18]. On the other hand, extrinsic cohesive models are inserted in the corresponding mesh after damage initiation is predicted to take place [9, 10]. Note that though the extrinsic approach prevents to introduce an artificial compliance in the system, which can alter the accuracy of the computations, it requires the development of a systematic procedure for their incorporation. Moreover, in 3D applications and parallel computations, the implementation of the extrinsic framework is not straightforward, suffering from different algorithmic difficulties. These latter aspects are difficult to be performed in most of the standard FE codes, especially in commercial packages. In this study, differing from the previous zero-thickness formulations, the interface between two adjoining bodies is considered to have a finite thickness, see [28, 29] and the references therein given. This approximation allows the description of the interface using a pseudo-continuum approach to be accomplished by means of assuming that the nonlinear failure phenomena take place within a very narrow and confined area [30, 31]. This is the key concept of the formulation proposed by Balzani and Wagner [32], in which the use of a degenerated 3D isoparametric continuum element formulation to model decohesion events is proposed, see also [33] for finite strain applications. Remarkable contributions in this area concern with the hierarchical formulation developed in [16], and more recent investigations accounting for the interface dilation [34] or incorporating the in-plane deformation of heterogeneous interfaces using multiscale methods [35]. Based on this finite thickness interface idealization, this paper is concerned with the development of a novel interface element for large deformation analysis based on the solid shell concept [36, 37, 38]. This kinematic description has been recently invoked in [39], and also used in contact applications under the denomination of contact layer formulation [40]. The geometric nonlinear effects are principally motivated by the current trend in most of composite structures (especially in the aeronautical and aerospace industries) for producing lightweight components. This is the case for instance of composite stiffened panels whose optimized designs are intended to be used for composite fuselage components undergoing postbuckling deformation patterns prior to collapsing point, see [41] and the references therein given. Consequently, due to the characteristics of the interface model herein envisaged and to assess the capabilities of the present formulation, the focus of the present study is twofold. First, with regard to the geometric and kinematic 2

ACCEPTED MANUSCRIPT

2. Finite thickness interface model

M AN U

SC

RI PT

description and differing from previous approaches, the current element formulation is completely defined in the curvilinear setting, in which the combined effects of material and geometrical nonlinearities are specifically taken into account. The adoption of the solid shell concept exploits the R6 -parametrization of the shell deformation [37] which accounts for displacements of two material points on the top and on the bottom surfaces of the body to describe the kinematic field. Second, regarding the constitutive definition of the interface, two types of standard irreversible interface laws are integrated into the element definition using a continuum-based approach in order to preserve the consistency with the interface formulation: (i) an exponential-based, and (ii) a polynomial-based laws. The manuscript is organized as follows. Section 2 outlines fundamentals of the proposed interface formulation. Particularly, Section 2.1 gives the geometric and kinematic description of the interface, whereas the corresponding FE formulation and implementation aspects are covered in Section 2.2. This novel interface finite element formulation accounting for both interface laws aforementioned is integrated into the general purposes FE codes FEAP [42] and ABAQUS [43] by means of user-defined routines. The description of both irreversible TSLs that define the inelastic response of the interface considered in this investigation is given in Section 3. Section 4 demonstrates the applicability of the developed formulation by means of several numerical applications. Numerical predictions are also correlated with available experimental data. Finally, the main conclusions of this investigation are provided in Section 5.

This section presents the description of the geometrical and kinematic aspects of the degenerated solid shell interface element for finite strain analysis and its corresponding finite element formulation. 2.1. Interface kinematics in the curvilinear setting (i)

AC C

EP

TE

D

Let consider two adjacent bulk regions B0 with i = 1, 2 in the reference configuration, corresponding (i) to the bulk regions Bt in the current configuration, respectively, see Figure 1. Between these two bulks, we assume the existence of an interface of finite thickness, which constitutes the body of interest (see Figure 1), whose reference configuration is denoted by B¯0 . The nonlinear deformation map is identified by ϕ(X) : B¯0 × [0, t] → R3 , where [0, t] is the time step interval, that maps the reference material points (X ∈ B¯0 ) onto the current material points (x ∈ B¯t ) 1 . The interface is idealized under the assumption stating that one of its geometric dimensions (the one identified with the thickness direction) is significantly smaller than the other two dimensions (the in-plane dimensions), see [29, 32, 40]. Specifically, Balzani and Wagner [32] proposed that the initial thickness of the interface element H0 should be about 1/100 the characteristic thickness of the adjoining bulks (dealing with thin walled structures). This adoption aims at minimizing the effects of the potential bending and torsional moments due to the eccentricities of the corresponding nodal forces in the subsequent finite element discretization. Consequently, a constant strain state along the thickness direction can be postulated. The local curvilinear co-variant convected basis in the reference and in the current configurations of the interface are defined as follows: Ga =

∂x(ξ) ∂X(ξ) ; ga = , a = s, t, n a ∂ξ ∂ξ a

(1)

where the subscripts s and t denote the in-plane directions, and the subscript n stands for the out-of-plane direction. In this setting, the interface is parametrized as follows: A(ξ s , ξ t , ξ n ) = A(ξ) = M ×F ⊂ R3 denotes the manifold of the body in which M (ξ s , ξ t ) ⊂ R2 and identifies the reference midsurface of the interface and F (ξ n ) ⊂ R denotes the thickness space, which are defined in the bi-unit cube  = [−1, 1] × [−1, 1] × [−1, 1]. 1 In the following, magnitudes in capital and small letter denote quantities in the reference and current configurations, respectively.

3

ACCEPTED MANUSCRIPT

Current configuration

Reference configuration Bulk-1

RI PT

Bulk-1

Bulk-2

SC

Bulk-2

Interface element: FE discretization 8

3

4

M AN U

5

2

n

1

1

t

s

6

7

3

2

Element node

D

Figure 1: Geometric definition of the interface between the two bulk regions in the reference and current configurations. Discretization of the three-dimensional finite thickness interface element.

TE

The dual contra-variant base vectors satisfy the standard relations: Ga · Gb = δab and ga · gb = δab , with a, b = s, t, n. The metric tensors in the curvilinear setting associated with the reference surface takes the form: G = Gab Ga ⊗ Gb , g = gab ga ⊗ gb . The vectors perpendicular to the interface midsurface in the reference and in the current configurations (see Figure 2) are defined as:

EP

¯n ¯n = G , G ¯ nn G

AC C

¯n = g

¯n g , g¯nn

¯ n = Gs × Gt ; G ¯n ·G ¯ n. ¯ nn = G where G ¯ n = gs × gt ; g¯nn = g ¯n · g ¯n. where g

(2) (3)

With regard to the kinematic description, the so-called solid shell concept is herein exploited [37, 38]. According to this parametrization, any material point of the interface is expressed in terms of a pair of material points located on the top and on the bottom surfaces of the interface. These top and bottom position vectors (being identified by the subscripts t and b, respectively) are denoted as {Xt , Xb } in the reference configuration and {xt , xb } in the current configuration. Relying on this approximation, the positioning vectors in the reference and in the current configurations with respect to the curvilinear setting read: X(ξ) =

1 1 (1 + ξ n ) Xt (ξ s , ξ t ) + (1 − ξ n ) Xb (ξ s , ξ t ) 2 2

(4)

1 1 (1 + ξ n ) xt (ξ s , ξ t ) + (1 − ξ n ) xb (ξ s , ξ t ). (5) 2 2 The current configuration is obtained through use of the the displacement field u as follows: x(ξ) := X(ξ) + u(ξ) (see Figure 1). Following [44], the kinematic field can be easily reformulated in terms of the x(ξ) =

4

ACCEPTED MANUSCRIPT

8 7 t

n

6

5

4

RI PT

s

3

1 2 3

Element node

SC

2

1

M AN U

Figure 2: Definition of the local basis of the degenerated solid shell interface model.

displacements of the reference surface of the interface v(ξ s , ξ t ) and the director vector w(ξ s , ξ t ): u(ξ) = v(ξ s , ξ t ) + ξ n w(ξ s , ξ t ). With the previous definitions at hand, the Jacobi matrices referred to the transformations between the parametric space in the reference, J(ξ), and in the current, j(ξ), configurations are defined as follows: T

J(ξ) = [Gs , Gt , Gn ] ,

T

j(ξ) = [gs , gt , gn ]

(6)

∂x ∂x ∂ξ a = a = ga ⊗ Ga ∂X ∂ξ ∂X

(7)

D

Correspondingly, the deformation gradient reads:

TE

F := ∇X ϕ =

where ∇X denotes the material gradient operator. The Green-Lagrange strain tensor is defined as:  1 T F F − I2 ; 2

Eij =

1 [gi · gj − Gi · Gj ] , 2

(8)

EP

E :=

AC C

where I2 stands for the material metric matrix. ¯ associated with the interlaminar strain The reduced version of the Green-Lagrange strain tensor, E, ¯ (corresponding to the second components, and their corresponding energetically conjugated stress tensor, S Piola-Kirchhoff stress tensor), are expanded using a vector notation as follows:        sn  ¯sn E gs · gn − Gs · Gn gsn − Gsn S¯ ¯ ¯ ¯        E = Etn = gt · gn − Gt · Gn = gtn − Gtn , S = S¯tn  (9) ¯ Enn gn · gn − Gn · Gn gnn − Gnn S¯nn From the kinematic definition of the interface herein adopted, Eqs.(4)-(5), we assume that any arbitrary transverse material fiber that connects the bottom and the top surfaces of the interface is not required to remain neither normal to the reference surface nor inextensible along the deformation process. This assumption has a direct consequence on the contribution of the interface deformation to the different fracture modes. In this respect, the convective basis vectors aligned with the directions s and t are associated with the interlaminar fractures modes II and III respectively (shear modes), whilst the normal vector to the n ¯– ¯ n ) regards the interlaminar mode I of failure (out-of-plane opening of the interface), which plane (vector G is perpendicular to the midsurface surface of the interface. Note that only in the case of planar structures ¯ n in the reference configuration (plates) and under very specific loading conditions, the vectors Gn and G 5

ACCEPTED MANUSCRIPT

RI PT

¯n in the current configuration are coincident. Therefore, in a general framework, a modification and gn and g of the metric components associated with such failure modes, i.e. gsn , gtn and gnn need to be accounted for, with special emphasis in curved interfaces. It is also worth mentioning that in contrast to alternative formulations, the current interface model does not assume any simplifying hypothesis with regard to the kinematics field, properly taking into account the consideration of the deflection of the interface along the deformation process. A more comprehensive discussion regarding the role of such effects has been addressed in [29]. Based on the previous arguments, a modification of the metric components related to the interlaminar actions is required. This modification results from the projection of the vectors Gn and gn onto the directions associated with each of the failure modes in both configurations, respectively. These alternative convected vectors in the reference and current configurations can be expressed as: (10)

¯s = gs + Pgs (gn ); g ¯t = gt + Pgt (gn ); g ¯n = Pg¯n (gn ); g

(11)

SC

¯ s = Gs + PG (Gn ); G ¯ t = Gt + PG (Gn ); G ¯ n = PG G ¯ n (Gn ); s t

where the operators Pr (q) denote the projection of the vector q onto the direction r. These operators rely on the standard vector projection relationship, which reads: q·r r. |r|2

M AN U Pr (q) =

(12)

D

The modified Green-Lagrange strain tensor based on the alternative definitions given in Eqs.(10)-(11) renders:      gn · gs Gn · Gs ¯     ¯ g + g G · g − G + · G s s n s n ¯s ·G ¯n ¯sn  s ¯s · g ¯n − G g E |gs |2 |Gs |2      ¯ ¯ ¯ ¯  .     Gn · Gt gn · gt ¯t · g ¯n − Gt · Gn =  E = Etn = g (13) ¯  ¯ g G · g − G + · G g + t t n t n t ¯ ¯ ¯  2 2 ¯n · g ¯n − Gn · Gn g Enn |gt | |Gt | ¯n ·G ¯n ¯n · g ¯n − G g

TE

Note that one of the most appealing aspects of the proposed interface model regards its complete kinematic compatibility with the adjoining bulk elements in simulations where solid shell or standard solid elements can be employed.

AC C

EP

2.2. Variational basis and finite element formulation In this section, the principal aspects of the variational formulation and the corresponding FE discretization of the proposed interface element are briefly outlined. The weak form of the interface contribution to the total mechanical potential energy of the system composed by the interface and the adjoining bulk regions is discussed in Section 2.2.1. Subsequently, the description of the finite element discretization and the consistent linearization procedure are addressed in Section 2.2.2, whereas the specification of the suitable operators that are required for the finite element implementation of the current formulation are introduced in Section 2.2.3. 2.2.1. Weak form of equilibrium In order to derive the equilibrium equations of the interface model, assuming isothermal conditions, the total potential energy of the system , ΠT , is composed by the internal contribution due to the continuum adjacent bulks, Πint,bulk , the internal contribution of the interface, Πint,cohe , and the external applied actions, Πext : ΠT (u) = Πint,bulk (u) + Πint,cohe (u) + Πext (u). (14) Invoking the principle of virtual work, the variational form of Eq.(14) can be expressed as: δΠT (u, δu) = δΠint,bulk (u, δu) + δΠint,cohe (u, δu) + δΠext (u, δu), 6

(15)

ACCEPTED MANUSCRIPT

where δu identify admissible virtual displacements. In the reference configuration, the variation of the term corresponding to the internal interface energy is given by: Z ¯ : δE ¯ dΩ, δΠint,cohe (u, δu) = S (16) B¯0

RI PT

¯ denotes the virtual reduced strain vector, whilst S ¯ stands for the reduced stress vector. The where δ E ¯ ¯ relationship E − S is governed by a suitable nonlinear interface law as is described in Section 3. The virtual variation of the modified Green-Lagrange strain tensor concerning the interlaminar components yields   ¯n + g ¯s · δ¯ δ¯ gs · g gn ¯ =  δ¯ ¯n + g ¯t · δ¯ gt · g gn  . δE (17) ¯n + g ¯n · δ¯ δ¯ gn · g gn

M AN U

SC

The variation of the metric components are computed through the variation of the displacements associated with the current interface formulation. In addition, the prevention of locking phenomena can be precluded by means of the combined use of the Enhanced Assumed Strain (EAS) and the Assumed Natural Strain (ANS) methods as was comprehensively described in [44], see also the corresponding procedure for their computational implementation outlined in [41]. Nevertheless, since the interface displacements is driven by the adjoining bulk elements, the use of numerical techniques to prevent locking pathologies are not considered in the current formulation. 2.2.2. Finite element discretization According to an iso-parametric map, the corresponding finite element formulation introduces the discretization of the interface through ne ∈ N solid shell-like interface elements. The standard tri-linear shape functions NA (ξ) take the form:


1 1 2 3 1 + ξ 1 ξA 1 + ξ 2 ξA 1 + ξ 3 ξA , A = 1, . . . , nn , 8

(18)

D

NA (ξ) =

TE

with ξ denoting the isoparametric domain and nn = 8 identifying the number of nodes of the element. The position vectors of any material point in the reference and current configurations can be respectively expressed as: nn nn X X X ≈ Xh = NA (ξ)XA ; x ≈ xh = NA (ξ)xA , (19) A=1

A=1

AC C

EP

where XA and xA stand for the nodal position vectors in the reference and current configurations, respectively. Note that in the sequel, discrete magnitudes are denoted by the superscript h. The real, virtual and incremental displacement vectors are interpolated using the previous shape functions according to the following scheme: u ≈ uh =

nn X

NA (ξ)dA = Nd; δu ≈ δuh =

A=1

nn X

NA (ξ)δdA = Nδd; ∆u ≈ ∆uh =

nn X

NA (ξ)∆dA = N∆d,

A=1

A=1

(20) where N and d are the suitable operators that collect the interpolation functions and the nodal displacements of the element; dA , δdA and ∆dA identify the real, virtual and incremental displacement vectors at the node A. The virtual and incremental strain vectors are accordingly defined as: ¯ ≈ δE ¯h = δE

nn X

¯ ≈ ∆E ¯h = BA δdA = Bδd; ∆E

A=1

nn X

BA ∆dA = B∆d,

(21)

A=1

where B denotes the so-called compatibility strain-displacement operator, which in the current formulation is function of the nodal displacements. 7

ACCEPTED MANUSCRIPT

The insertion of the expressions Eq.(19)-(21) into the weak form of the term associated with the internal energy of the interface, Eq.(16), yields to: Z  h T T¯ δΠint,cohe (d, δd) = δd B S dΩ = δdT R, (22) B¯0

RI PT

where R is the residual vector at the element level. The consistent linearization of Eq.(22) through the use of the directional derivative concept [45] in the direction ∆d is given by: # "Z  T Z ¯ ∂ S ∂B T h T ¯ dΩ + S dΩ ∆d = δdT [K] ∆d, B (23) ∆δΠint,cohe (d, δd, ∆d) = δd ∂d ∂d ¯ ¯ B0 B0

M AN U

SC

with K identifying the element stiffness matrix. The linearization of the stress vector can be accomplished using the material tangent operator Cc of the interface constitutive law (see Section 3): Z Z ¯ ∂S ¯ dΩ = (24) BT ∆S BT Cc B dΩ, with Cc = ¯ ∂E B¯0 B¯0 Through the integration of Eq.(24) into the definition of the element stiffness matrix, Eq.(23), one can identify the classical material Kmat and geometrical Kgeom contributions, which can be expressed as: Z Kmat =

B¯0

T

Z

B Cc B dΩ; Kgeom =

B¯0



∂B ∂d

T

¯ dΩ. S

(25)

TE

D

Note that in case of geometrically linear applications, the reference and the current configurations are coincident, and therefore the geometrical contribution to the element stiffness matrix vanishes. Finally, it is worth mentioning that the standard Gauss integration rule with 2×2×2 integration points is used to perform the numerical integration of the element matrices previously defined.

AC C

EP

2.2.3. Implementation aspects This section briefly describes several operators that are required for the numerical implementation of the current solid shell-like interface element. The discrete version of the curvilinear basis vectors in the current configuration renders:     nn nn XA dAx X X ∂x ∂u NA,i  YA  + NA,i dAy  , with i = s, t, n. (26) gi = i = Gi + i ≈ ∂ξ ∂ξ ZA dAz A=1 A=1 In Eq.(26), NA,i identifies the partial derivative of the shape function NA with respect to the natural coordinate ξ i . The partial derivative of the vectors gi with respect to the displacement vector dk of the node k takes the form:   ∂gix   ∂dk   Nk,i 0 0  ∂g  ∂gi  iy   0 Nk,i 0  , with i = s, t, n and k = 1, . . . , nn . = (27) =  ∂dk  ∂dk 0 0 Nk,i  ∂g  iz

∂dk ¯n , that is perpendicular to the midsurface of the Similarly, the partial derivative of the modified vector g interface in the deformed configuration, with respect the displacement vector, dk , corresponding to the node 8

ACCEPTED MANUSCRIPT

k is given by:   ∂¯ gnx  ∂dk   0  ∂¯  ¯n ∂g  gny   =  = Nk,t gsz − Nk,s gtz  ∂dk  ∂dk Nk,s gty − Nk,t gsy  ∂¯ gnz  ∂dk

 Nk,t gsy − Nk,s gty Nk,s gtx − Nk,t gsx  , with k = 1, . . . , nn . (28) 0

RI PT

Nk,s gtz − Nk,t gsz 0 Nk,t gsx − Nk,s gtx

SC

With these definition at hand, the strain-displacement operator B(d) in Eq.(22) can be recast as:  ∂ g¯  ¯n g ¯sT ∂∂d ¯nT ∂dks + g g k ¯ ∂E  T ∂ g¯t ¯n  g ¯n ∂dk + g ¯tT ∂∂d B= (29) = g . k ∂dk ¯n ∂g T ¯n ∂dk g

M AN U

The geometric contribution to the stiffness matrix incorporates the derivative of the B-operator with respect to the kinematic field, see Eq. (25)2 . This derivative is computed as: h  h iT iT 2 2 ¯s ¯s ¯s ¯n ¯n ¯n ∂g ∂g ∂g ∂g T ∂ g T ∂ g ¯ ¯ + + g + g n ∂dk ∂dm s ∂dk ∂dm  ∂d ∂d h ∂dm i ∂dk h m iT k ¯   T ∂B ∂2E 2 2 ¯ ¯ ¯ ¯ ¯ ¯ n  , with k, m = 1, . . . , n . ∂ g ∂ g ∂ g ∂ g ∂ g g  T n t t t n = =  ∂d ¯n ∂dk ∂dm + ∂dm ¯tT ∂d∂k ∂d +g +g n  ∂d ∂d m m k k ∂dm ∂dk ∂dm   h iT 2 ¯n ¯n ∂g ∂g ¯ T ∂ g¯n ∂dm ∂dk + gn ∂dk ∂dm

3. Interface constitutive law

D

(30) The term expressed in Eq.(30) yields to a third-order tensor that needs to be consistently multiplied by the interlaminar stress vector in order to obtain the final form of the geometric term of the element stiffness matrix.

EP

TE

The simulation of delamination effects at the interface requires the incorporation of a suitable inelastic constitutive law. Without loss of generality, this section describes two types of nonlinear decohesion laws, which are accordingly modified for their integration into the current interface formulation: (i) the exponential-based constitutive law proposed by Ortiz and Pandolfi [10], whose modification for finite thickness interface formulations was derived in [32], and (ii) the polynomial law proposed by Tvergaard [5]. For a clear identification of the contribution of the intelaminar components of the Green–Lagrange strain tensor and their corresponding second Piola-Kirchhoff stress tensor to the fracture modes I, and II, and III, the following notation is adopted:

AC C

    ¯ = E ¯ = S¯sn , S¯tn , S¯nn T = [SII , SIII , SI ]T . ¯sn , E ¯tn , E ¯nn T = [EII , EIII , EI ]T ; S E

(31)

Note however, that the stress components introduced in Eq.(31) should be pushed forward to the current configuration in order to obtain the Cauchy stress tensor σ: σ=

1 FSFT , det[F]

(32)

where det[F] identifies the determinant of the deformation gradient. 3.1. Exponential-based interface law The exponential-based interface model relies on the universal binding relation by Rose et al. [47], which has been successfully applied to a wide range of engineering applications such as crack propagation in metals [9, 10] and delamination in composite structures [32], among many others. We postulate the existence 9

ACCEPTED MANUSCRIPT

¯ of a cohesive free energy as function of the interlaminar components of the strain tensor, ψc = ψc (E). Following the previous formulation proposed in [32], this interface free energy can be expressed in terms of the equivalent strain E that is defined as: q 2 2 + E 2 ), E := hEI i + β (EII (33) III

RI PT

where the parameter β accounts for different weights for peeling fracture mode (mode I) and shear modes (modes II and III), and the operator h•i denotes the negative Macaulay brackets. Based on the previous definition, the cohesive free energy that describes the decohesion at the interface reads:     E −E/E 0 0 0 , (34) ψc (E) = eS E 1 − 1 + 0 e E

M AN U

SC

where e = exp[1] denotes the Euler number, S 0 is critical equivalent stress and E 0 identifies the critical equivalent strain for the initiation of decohesion. Figure 3 portrays the response for the single modes I and II, according to the decohesion law given in Eq.(34). Note that the area under the stress-strain curve is set equal to the fracture toughness of the corresponding mode divided by the interface thickness H0 . This also holds for mixed-mode conditions, so that Gc = lim ψc (E). E→∞ H0

(35)

−

EP

TE

D

Gc∗ =

−

AC C

Figure 3: Interface response corresponding to the fracture Modes I (left) and II (right) according to the exponential-based model [32].

We also introduce the definition of the internal variable Emax that identifies the maximum decohesion strain along the loading history. Therefore, the loading/unloading relations can be expressed as: E = Emax E < Emax

E˙ > 0 for loading E˙ ≤ 0 for unloading/reloading

(36)

where E˙ denotes the time derivative of the equivalent strain. The interlaminar stress components are ¯ with respect to the interlaminar strains as follows: computed by taking partial derivatives of ψc (E) ¯ = ∂ψc (E) ∂E for loading. S ¯ ∂E ∂ E 10

(37)

ACCEPTED MANUSCRIPT

¯ = ∂ψc (Emax ) ∂E for unloading/reloading. S ¯ ∂E ∂E The second Piola-Kirchhoff stress tensor can be accordingly expressed in vector notation as:

(39)

RI PT

T  ¯ = [SII , SIII , SI ]T = ∂ψc (E) ∂E , ∂ψc (E) ∂E , ∂ψc (E) ∂E . S ∂E ∂EII ∂E ∂EIII ∂E ∂EI

(38)

∂EII

∂EIII

SC

The linearization of the interlaminar stress components yields to the definition of the tangent constitutive operator of the interface Cc :   ∂SII ∂SII ∂SII  ∂EII ∂EIII ∂EI     ∂SIII ∂SIII ∂SIII   . (40) Cc =  ∂EI   ∂EII ∂EIII   ∂SI ∂SI ∂SI  ∂EI

M AN U

We also define a damage variable to trigger the progression of the delamination along the loading process, which is defined as: ψc (Emax ) d= (41) Gc∗ This damage variable ranges from 0 to 1, identifying intact and fully damaged interface states, respectively. The interface model is further equipped by introducing a penalty contact condition to avoid interpenetration, whose stiffness k reads [32]: k=

eS 0 E0

(42)

TE

D

Under mixed mode fracture conditions, the interface failure is attained using the 3D version of the Benzeggah-Kenane failure criterion [46]:  η GII + GIII Gc = GIc + (GIIc − GIc ) , (43) GI + GII + GIII

being

AC C

EP

where η is a fitting parameter with respect to the single mode experiments, namely double cantilever beam (DCB) and end notch flexure (ENF) for mode I and mode II conditions, respectively, and mixed mode bending (MMB) for mixed mode scenarios. However, note that any other fracture criterion can be easily integrated into the present formulation. In the above criterion, Eq.(43), due to the lack of experimental data that accurately determine the fracture toughness of mode III, identical fracture toughness for modes II and III is assumed, i.e. GIIc = GIIIc . In line with [32], the BK criterion can be expressed as follows:

Gshear β¯ = = GT

Gc = GIc + (GIIc − GIc ) β¯η 

Eshear E

2 ,

Eshear =

q 2 + E 2 ), β (EII III

(44)

(45)

where Gshear = GII + GIII and GT = GI + GII + GIII identify the shear and total energy release rates, respectively.

11

ACCEPTED MANUSCRIPT

RI PT

3.2. Polynomial-based interface law The second interface law under consideration in this study regards the modified version of the polynomial law proposed by Tvergaard [5]. In order to take into account the interaction between the different fracture modes, the effective dimensionless strain λ and the factor P (λ) are defined as follows: s 2  2  2 hEI i EII EIII + + , (46a) λ= EIc EIIc EIIIc (
27 1 − 2λ + λ2 , for 0 ≤ λ ≤ 1 , (46b) P (λ) = 4 0, otherwise

SC

where EIIc , EIIIc and EIc are the critical relative strains corresponding to the fracture Modes II, III and I, respectively. ¯ = [SII , SIII , SI ]T is a The cohesive traction vector associated with the interlaminar components S nonlinear polynomial function of the corresponding strain components at the interface: EII P (λ) EIIc EIII SIII =SIII,max P (λ) EIIIc hEI i SI =SI,max P (λ) EIc

M AN U

SII =SII,max

(47a) (47b) (47c)

TE

D

where SII,max , SIII,max and SI,max identify the maximum cohesive traction components corresponding to the fracture Modes II, III and I, respectively. Figure 4 depicts the normal traction SI vs. EI . Similar curves hold for the two other fracture modes, whose representations are omitted here for conciseness reasons.

/

,

AC C

EP

1

1

/

Figure 4: Schematic representation for interface response for Mode I according to the polynomial law [5].

12

ACCEPTED MANUSCRIPT

SC

RI PT

The constitutive tangent derived from the stress vector given in Eq.(47) reads:    P (λ) EII ∂P (λ) ∂P (λ) ∂P (λ) + SII,max SII,max  SII,max E E ∂E ∂E ∂EI IIc IIc II III     P (λ) EIII ∂P (λ) ∂P (λ) ∂P (λ)  Cc =  SIII,max + SIII,max SIII,max  ∂EII EIIIc EIIIc ∂EIII ∂EI     ∂P (λ) ∂P (λ) P (λ) EI ∂P (λ) SI,max SI,max SI,max + ∂EII ∂EIII EIc EIc ∂EI (48) The detailed form of the terms of the tangent operator Cc in Eq.(48) can be found in [26] adopting a zero-thickness interface model. Finally, the current interface law also accounts for a penalty condition for compressive Mode I in order to avoid interpenetration between the two flanks of the interface, see Figure 4. 4. Representative applications

M AN U

The practicability of the proposed interface formulation is demonstrate by means of several examples. The current model has been implemented into the general purpose FE codes FEAP [42] and ABAQUS [43] through user-defined elements. Fiber-reinforced composite entities are discretized using SC8R solid shell elements of ABAQUS (any other solid shell formulation can be easily incorporated [48]), therefore a complete kinematic compatibility in terms of nodal degrees of freedom between the bulk and interface elements is accomplished. This fact prevents the use of additional numerical constraints of any alternative numerical scheme to interrelate the bulk and interface nodal degrees of freedom. The simulations include both geometrical and material nonlinearities using a displacement-control solution scheme to solve the nonlinear governing equations. The current model is generated using the complete geometric description of each specimen under consideration.

AC C

EP

TE

D

4.1. Double cantilever beam (DCB) test The first application to evaluate the performance of the proposed element formulation is the standard double cantilever beam (DCB) test [39]. This test has been widely used for the determination of the mode I fracture toughness, see [6, 32] and the references therein given. The geometry, load and boundary conditions are shown in Figure 5, replicating the data reported in [32] whereby quasi-static loading conditions are assumed for the corresponding computations. The geometric dimensions are: (i) length 210 mm, (ii) width 20 mm, and (iii) laminate thickness 2h = 1.5 mm (2 × 12 composite layers, being the individual layer thickness equal to 0.125 mm and oriented along the reference direction, see Figure 5). An initial pre-crack with 45 mm in length (Teflon layer of 40 mm and an additional pre-crack of 5 mm in length) is introduced. The reference material orientation is aligned with the longitudinal direction (x-direction) of the specimen. The specimen is manufactured from the unidirectional carbon fiber reinforced polymer (CFRP) composite material IM7/8552, whose mechanical properties are reported in Tables 1 and 2. A series of numerical computations using single fracture mode models are performed to adjust the fracture behavior of the exponential cohesive law (Section 3.2) to the fracture parameters given in Table 2. E11 (GPa) 144

E22 (GPa) 7.7

G12 (GPa) 5.9

G13 (GPa) 5.9

ν12 0.3

thickness (mm) 0.125

Table 1: Mechanical properties of IM7/8552.

σ c (MPa) 51.7

τ c (MPa) 115.7

GIC (J/m2) 270

GIIC (J/m2) 687

η 2.3

Table 2: Interlaminar fracture properties for IM7/8552.

13

H0 (mm) 0.00125

        

ACCEPTED MANUSCRIPT

a

z

y

x

0º 165 u

b

45

RI PT

clamped side

M AN U

SC

mesh detail

Figure 5: Double cantilever beam (DCB) application. (a) Geometric and loading description. (b) FE mesh.

AC C

EP

TE

D

The beam is loaded at the free edge prescribing a uniform displacement along the z-direction according to the frame depicted in Figure 5.a, whilst the displacements are restrained at the opposite longitudinal side. Both interface laws (Sections 3.1 and 3.2) are used to simulate the current application. Preliminary simulations are performed to asses the mesh dependency leading to the use of a graded mesh along the longitudinal direction based on the following guidelines: region 1 is discretized with 15 elements, while region 2 consists of 240 elements. The width of the beam is discretized using 20 elements. The final characteristics of the resulting mesh are: (i) 7400 SC8R elements of ABAQUS for the bulks, and (ii) 3000 interface elements. Note that a single element over the laminate thickness is considered, see the detail in Figure 5.b. The developed interface element are inserted between the two portions of the composite beam in region 2. The initial thickness of the interface is set equal to H0 = 0.00125 mm. Computations are performed by the incremental-iterative Newton–Raphson method with an initial and maximum load step of 1% of the final displacement applied. The final load is reached without exhibiting remarkable difficulties in achieving equilibrium solutions. Figure 6.a shows the load-deflection curves of the current FE simulations using the exponential- (labelled as Exp Law) and the polynomial-based (labelled as Polyn Law) interface laws herein considered and the builtin cohesive interface formulation complying with the proposed finite thickness formulation. The cohesive implementation of ABAQUS relies on the formulation developed by Camanho and coauthors [6] using specific interface elements COH3D8, which has been extensively used by several authors for the validation purposes for different cohesive formulations, see [39] and the references therein given. In particular, a linear softening for the decohesion evolution is employed for the definition of the ABAQUS interface model (labelled as ABAQUS Bilinear). The required penalty parameters for the initial slope of the bilinear TSL law replicate those reported in [24]. In this graph it can be seen that almost identical strength values are predicted for both interface models, which are in excellent agreement with the experimental data (see the shading area) and the computations 14

ACCEPTED MANUSCRIPT

60

a

50

Exp. data [50]

30

20

Num Exp Law Num Polyn Law

10

ABAQUS Bilinear 0 0

2.5

5

7.5

10

crack opening [mm]

15

17.5

20

SC

b

12.5

RI PT

Force [N]

40

M AN U

direction of propagation

TE

D

Figure 6: Double cantilever beam (DCB) application. (a) Load-deflection evolution curves: correlation between the experimental and the numerical data using both proposed interface laws and the built-in ABAQUS cohesive formulation obeying a bilinear TSL. (b) Contour plot of the damage variable corresponding to the exponential decohesion law for a crack opening equal to 15 mm.

AC C

EP

reported in [32, 33]. Only slight differences are observed at the initial stage of the post-peak evolution, the polynomial law slightly underestimating the peak point of the exponential interface model. With regard to the interface model of ABAQUS, a sensitive dependency of the initial stage of the evolution (before the softening branch) upon the values for the initial stiffness is observed. Nevertheless, the corresponding data lie within the experimental range and are in close agreement with the numerical estimation using the formulation herein proposed. Finally, a contour plot corresponding to the damage variable defined in Section 3.1 is shown in Figure 6.b for an intermediate state. This internal variable predicts a curved delamination front through the width of the beam as was previously discussed in [32, 33]. 4.2. Rib-sttiffened composite test The second example under consideration is concerned with the analysis of the skin-stringer separation of a composite rib-stiffened specimen, see Figure 7. The structural set consists of: (i) a stringer with Ttransverse section, with flange width and web height are equal to 60 mm and 20 mm, respectively, and (ii) a curved skin with inner radius equal to 938 mm and 100 mm in width. The stringer and the skin have a longitudinal length (x-direction) of 250 mm. Similarly to the standard DCB test, an initial pre-crack of 43 mm in length was inserted via a Teflon layer. The experimental setup also includes two aluminium block, which are tied to the specimen in order to apply the external loading. Further details with regard to the experimental investigation of the current application can be found in [49, 50]. 15

ACCEPTED MANUSCRIPT

z x

y

SC

0º

RI PT

The stringer and the skin are manufactured from the unidirectional CFRP material IM7/8552, whose mechanical properties are reported in Tables 1 and 2. The reference material direction is coincident with the longitudinal direction. The laminated disposal of the skin and the stringer are [0/45/-45/90]S and [45/-45/0/0]3S , respectively. Due to the lack of symmetry, this applications is used in order to asses the mixed–mode behavior of the proposed interface formulation.

R938

60 100

14

40

M AN U

250

20

R2

auxiliary aluminum blocks

x

z Lengths in mm

TE

D

F

F

3

EP

Figure 7: Rib-stiffened application: geometric and loading description.

AC C

In line with the previous example, a graded mesh is generated, which is defined in three different areas along the x-direction of the specimen (Figure 8.a): (i) region 1 (pre-crack) is discretized using 15 elements, (ii) 240 elements are used in region 2, and (iii) region 3 is meshed using 10 elements. The skin and the stringer are meshed using a solid shell-based mesh with 38425 SC8R elements of ABAQUS. The interface between the skin and the stringer in the region 2 is discretized using 8160 elements, whereas the interface in the region 3, which does not include any damage modeling capability, is meshed using 340 elements C3D8 of ABAQUS. Note that again a complete kinematic compatibility with regard to the nodal degrees of freedom between the bulk and the interface elements is achieved. The aluminum blocks mentioned above are not considered in the current simulations, therefore their effects are taken into consideration via the definition of the corresponding boundary conditions. The applied boundary conditions aim at replicating the experimental setup in [50], according to [49] the following restrictions are defined: • Nodal displacements along the x-direction for the skin and the stringer are restrained for x = 243 mm. • The external loading is imposed as a prescribed displacement onto the external surface of the skin in the region 1 where the aluminum blocks are located, see Figure 4.2. 16

ACCEPTED MANUSCRIPT

• Nodal displacements along the z-direction for the stringer flanges at the regions of the aluminum blocks are constrained.

M AN U

SC

RI PT

The initial and maximum time step of the analysis is set equal to 1% of the final displacement. Differing from [49] and complying with the developed model, geometric nonlinear effects are considered in the current simulations. Computations concluded by achieving the final imposed displacement with 158 increments and 1495 iterations in total. Figure 8.b depicts the load–crack opening evolution curve corresponding to the example under investigation (the crack opening is identified with the displacement of the skin at x = 243 mm). In this graph, a satisfactory agreement between the experimental and numerical data can be observed. Indeed, a closer agreement with respect to that reported in [49] (labelled as WG[47] in the graph) between the proposed interface model and the experimental data is achieved along the initial elastic path before the delamination propagation. This fact can be attributed to two main issues: (i) the role of considering geometric nonlinear effects into the performed simulations and (ii) the complete kinematic compatibility of the bulk and interface elements. Continuing the analysis of the simulated response, the post-peak evolution (degradation) paths are in excellent agreement with the experimental data [49, 50]. Similar estimations are obtained with the polynomial-based interface model herein addressed (Section 3.2) and using the cohesive formulation of ABAQUS introduced in the previous section. The results corresponding to these computations are not included in the present graph for the sake of clarity. Figure 9.a shows deformed configurations of the specimen for different values of the crack opening, where an increasing skin-stringer separation can be observed. Finally, the estimation of the corresponding damage variable is depicted in Figure 9.b. Note again that a curved delamination is predicted due to the different mixed mode conditions along the delamination front in the present application. 5. Concluding remarks

AC C

EP

TE

D

A new nonlinear finite thickness cohesive interface element for modeling delamination in fibre-reinforced composite laminates and similar damage events has been developed. The proposed formulation relied on a degenerated solid shell concept, which is especially suitable for the analysis of thin-walled applications. One of the most appealing aspects of the model regarded the incorporation of geometric and material nonlinear effects into the analysis. Interlaminar damage has been considered by means of accommodating a modified version of two popular interface formulations, particularly the exponential-based model derived in [10] and the polynomial formulation addressed in [5]. Note however that any other interface law can be easily integrated into the proposed framework. The relative simplicity of the developed model allowed its implementation into general purpose nonlinear FE codes such as ABAQUS [43] and FEAP [42], among others, to be carried out. The accuracy of the current interface element has been assessed by means of two exemplary applications regarding unidirectional CFRP composites: (i) the standard double cantilever beam (DCB) test, and (ii) a rib-stiffened specimen. Numerical predictions were in close agreement with the experimental data. Simulations were computed without exhibiting remarkable difficulties in achieving equilibrium states, highlighting the robustness of the proposed formulation. The effectiveness of the proposed element lies in the following aspect: (i) the high versatility to integrate any interfacial law, (ii) the possibility of considering both geometric and material nonlinear effects in a consistent manner, (iii) the complete kinematic compatibility between the adjacent bulk (solid shell) elements, which is especially relevant for the analysis of thin-walled composite structures that are discretized using efficient solid shell formulations. Acknowledgments JR and AB gratefully acknowledge the financial support of the Andalusian Government (Project of Excellence No. TEP-7093) and the Spanish Ministry of Economy and Competitiveness (MAT2015-71036-P). MP 17

ACCEPTED MANUSCRIPT

a

RI PT

Stringer

b

SC

Skin

300

[49] WB [47]

250

M AN U

WB [47] [49]

Force [N]

200

Num Exp Law

150

100

Exp. data [50]

50

0 10

20

30

D

0

40

50

60

70

80

TE

crack opening [mm]

EP

Figure 8: Rib-stiffened application: (a) Mesh detail. (b) Load-crack opening evolution curve: comparison between the experimental data and numerical predictions.

AC C

would like to thank the European Research Council for supporting the ERC Starting Grant “Multi-field and multi-scale Computational Approach to Design and Durability of PhotoVoltaic Modules” - CA2PVM, under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 306622. References [1] [2] [3] [4] [5] [6] [7] [8]

Krueger R (2002) The virtual crack closure technique: history, approach and applications. NASA/CR-2002-211628. Dugdale D (1960) Yielding of steel sheets containing slits. J Mech Phys Solids 8(2):100–104. Barenblatt G (1962) The mathematical theory of equilibrium cracks in brittle fracture. Adv. Appl. Mech. 7, 55-129, 1962. Alfano G (2006) On the influence of the shape of the interface law on the application of cohesive-zone models. Compos Sci Techn, 66(6):723–730. Tvergaard V (1990) Effect of fiber debonding in a whisker-reinforced metal, Mat Sci Eng A, 107:23–40. Camanho PP, D´ avila CG, de Moura MF (2003). Numerical simulation of mixed-mode progressive delamination in composite materials. J. Compos. Mater. 37, 1415–1438. Harper PW, Hallett SR (2008) Cohesive zone length in numerical simulations of composite delamination. Eng Fract Mech, 75(16):4774–4792. Allix O, Corigliano A, 1999. Geometrical and interfacial non-linearities in the analysis of delamination in composites. Int J Solids Struct 36, 2189–2216.

18

ACCEPTED MANUSCRIPT

a

crack opening: 14 mm

crack opening: 0 mm

crack opening: 20 mm

RI PT

crack opening: 50 mm

crack opening 14 mm

M AN U

SC

b

crack opening 20 mm

crack opening 50 mm

D

Figure 9: Rib-stiffened application: (a)Deformed shape at different crack opening values. (b) Contour plot of the damage variable for different crack opening values.

AC C

EP

TE

[9] de-Andr´ es A, P´ erez JL, Ortiz M (1999) Elastoplastic finite element analysis of three-dimensional fatigue crack growth in aluminium shafts subjected to axial loading. Int J Solids Struct 36:2231–2258. [10] Ortiz M, Pandolfi A (1999) Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J Num Meth Eng 44:1267–1282. [11] de Borst R, Schipperen JHA (2002)Computational methods for delamination and fracture in composites. In: Allix O, Hild F, editors. Continuum damage mechanics of materials and structures. Elsevier Science Ltd. p. 325–352. [12] Sprenger W, Gruttmann F, Wagner W (2000) Delamination growth analysis in laminated structures with continuum based 3D-shell elements and a viscoplstic softeniung model. Comput Methods Appl Mech Eng 185:123–139. [13] Park K, Paulino GH, and Roesler JR (2009) A unified potential-based cohesive model of mixed-mode fracture. J Mech Phys Solids 57 (6), 891–908. [14] Xu XP, Needleman A (1993) Void nucleation by inclusion debonding in a crystal matrix. Modell Simul Mater Sci Eng 1(2):111–132. [15] van den Bosch MJ, Schreurs PJG, Geers MGD (2006) An improved description of the exponential Xu and Needleman cohesive zone law for mixed-mode decohesion. Eng Fract Mech 73(9):1220–34. [16] Paggi M, Wriggers P (2012) Stiffness and strength of hierarchical polycrystalline materials with imperfect interfaces. J Mech Phys Solids, 60:557–572. [17] Paggi M, Wriggers P. A nonlocal cohesive zone model for finite thickness interfaces part I: mathematical formulation and validation with molecular dynamics. Comput Mater Sci 2011;50(5):1625–33. [18] Paggi M, Wriggers P. A nonlocal cohesive zone model for finite thickness interfaces part II: FE implementation and application to polycrystalline materials. Comput Mater Sci 2011;50(5):1634–1643. [19] Prandtl L (1933) Ein Gedankenmodell f¨ ur den Zerreivorgang spr¨ oder K¨ orper (A thought model for the fracture of brittle solids). Zeitschrift f¨r Angenwandte Mathematik und Mechanik 13(2):129133. [20] Bruno D, Greco F, Lonetti P (2003) A coupled interface-multilayer approach for mixed mode delamination and contact analysis in laminated composites. Int J Solids Struct 40:7245–7268. [21] T´ avara L, Mantic V, Graciani E, Par´ıs F (2011) BEM analysis of crack onset and propagation along fiber-matrix interface under transverse tension using a linear elastic-brittle interface model. Eng Anals Boundary Elements 35:207–222. [22] Cornetti P, Mantic V, Carpinteri A (2012) Finite Fracture Mechanics at elastic interfaces. Int J Solids Struct 49:1022–1032.

19

ACCEPTED MANUSCRIPT

AC C

EP

TE

D

M AN U

SC

RI PT

[23] T´ avara L, Mantiˇ c V, Graciani E, Ca˜ nas J, Par´ıs F. (2010) Analysis of a crack in a thin adhesive layer between orthotropic materials. An application to composite interlaminar fracture toughness test. Comp Model in Eng Sc 58 (3) 247–270. [24] Reinoso J, Bl´ azquez A, T´ avara L, Par´ıs F, Arellano C (2016) Damage tolerance of composite runout panels under tensile loading. Comp Part B: Engineering 96, 79–93. [25] van den Bosch MJ, Schreurs PJG, Geers MGD (2007) A cohesive zone model with a large displacement formulation accounting for interfacial fibrilation. European J Mech A/Solids 26:1–19. [26] Reinoso J, Paggi M (2014) A consistent interface element formulation for geometrical and material nonlinearities, Comput Mech, 54:1569–1581. [27] Mosler J, Scheider I (2011). A thermodynamically and variationally consistent class of damage-type cohesive models. J Mech Phys Solids 59 (8),1647–1668. [28] Giambanco G, Scimeni GF, Spada A (2012) The interphase finite element. Comp Mech 50:353–366. [29] Freddi F, Sacco E (2014) An interface damage model accounting for in–plane effects. Int J Solids Struct 51:4230–4244. [30] Sarrado C, Leone F, Tur´ on A (2016) Finite-thickness cohesive elements for modeling thick adhesives. Engng Fract Mech, http://dx.doi.org/10.1016/j.engfracmech.2016.03.020 [31] Corrado M, Paggi M (2015) Nonlinear fracture dynamics of laminates with finite thickness adhesives. Mech of Mat 80, 183–192. [32] Balzani C, Wagner W (2008) An interface element for the simulation of delamination in unidirectional fiber-reinforced composite laminates. Eng Fract Mech 75:2597–615. [33] Areais P, Rabczuk T, Camanho PP (2016) Finite–strain laminates: Bending–enhanced hexaedron and delamination. Compos Struct 139:277–290. [34] Sørensen BF, Goutianos S (2014) Mixed mode cohesive law with interface dilatation. Mech Mater, 70(0):76–93. [35] Arag´ on AM, Soghrati S, Geubelle PH (2013) . Effect of in-plane deformation on the cohesive failure of heterogeneous adhesives. J Mech Phys Solids, 61(7): 1600 1611. [36] Parisch H (1995) A continuum-based shell theory for non-linear applications. Int J Numer Methods Eng 38:1855–1883. [37] Miehe C (1998) A theoretical and computational model for isotropic elastoplastic stress analysis in shells at large strains. Comp Meth Appl Mech Eng 155, 193–233. [38] Klinkel S, Wagner W (1997) A geometrical nonlinear brick element based on the EAS method, Int J Num Meth Eng, 40:4529–4545. [39] Simon J-W, Hoewer D, Stier B, Reese S (2015) Meso-mechanically motivated modeling of layered fiber reinforced composites accounting for delamination. Compos Struct, 122:477–487. [40] Weißenfels C, Wriggers P. (2015) A contact layer element for large deformations. Comp Mech http://dx.doi.org/10.1007/s00466-015-1140-7. [41] Reinoso J, Bl´ azquez A (2016) Application and finite element implementation of 7-parameter shell element for geometrically nonlinear analysis of layered CFRP composites. Compos Struct 139:263–276. [42] O.C. Zienkiewicz, R.L. Taylor (2000) The Finite Element Method. Butterworth-Heinemann, Woburn, MA, 5th Edition, Vol. I. ISBN: 0750650494. [43] ABAQUS (2011) ABAQUS Documentation, Dassault Syst` emes, Providence, RI, USA. [44] Bischoff M, Ramm E (1997) Shear Deformable Shell Elements for Large Strains and Rotations. Int J. Numer Meth. Engng, 40:4427–4449. [45] Bonet J, Wood RD (1997) Nonlinear continuum mechanics for finite element analysis. Cambridge University Press. [46] Benzeggagh ML, Kenane M. (1996) Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixedmode bending apparatus. Comp Sci Tech 56(4):43949. [47] Rose JH, Ferrante J, Smith JR, (1981) Universal binding energy curves formetals and bimettallic interfaces. Phys Rev Lett 47:675–678. [48] Reinoso, J., Paggi, M., Rolfes, R. (2016) A computational framework for the interplay between delamination and wrinkling in functionally graded thermal barrier coatings, Computational Materials Science 116:82–95. [49] Wagner W, Balzani C (2008) Simulation of delamination in stringer stiffened fiber-reinforced composite shells, Computers and Structures 86(9), 930–939. [50] Korjakins A, Ozolinsh O, Rikards R. (2005) COCOMAT Technical Report, WP 2, task 2.2: investigation of degradation by tests and development of degradation models. Final Report. Riga Technical University, Institute of Materials and Structures, Latvia.

20