A nonlocal finite element model for progressive failure analysis of composite laminates

Accepted Manuscript A nonlocal finite element model for progressive failure analysis of composite laminates P.F. Liu, Z.P. Gu, Y.H. Yang, X.Q. Peng PII:

S1359-8368(15)00625-3

DOI:

10.1016/j.compositesb.2015.09.061

Reference:

JCOMB 3838

To appear in:

Composites Part B

Received Date: 2 July 2015 Revised Date:

19 September 2015

Accepted Date: 30 September 2015

Please cite this article as: Liu PF, Gu ZP, Yang YH, Peng XQ, A nonlocal finite element model for progressive failure analysis of composite laminates, Composites Part B (2015), doi: 10.1016/ j.compositesb.2015.09.061. 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 nonlocal finite element model for progressive failure analysis of composite laminates

RI PT

P.F. Liu a, *, Z.P. Gu a, Y.H.Yang a, X.Q. Peng b a. Institute of Chemical Machinery and Process Equipment, Zhejiang University, Hangzhou 310027, China b. School of Materials Science and Engineering, Shanghai Jiaotong University, Shanghai 200030, China

SC

Abstract Localized damage and failure of composites often lead to spurious mesh dependence and zero energy dissipation due to the lack of a length scale by using finite element analysis.

M AN U

Based on the variational form of the stress boundary value problem with an interface discontinuity, this paper originally proposed a finite element model on the nonlocal intralaminar damage and interlaminar delamination of composite laminates. First, strain-based initial failure criteria and damage evolution criteria are used for the fiber and matrix, and the delamination is modeled by bilinear cohesive model. It is shown two length scales are introduced into the damage evolution

TE D

behaviors of the fiber and matrix to solve the localization problem of composites. Second, the global nonlocal stiffness equation is derived including the nonlocal intralaminar consistent tangent stiffness and the cohesive interface stiffness, and the node forces for bulk and interface elements.

EP

Third, numerical algorithm and finite element codes are developed by combining ABAQUS-UMAT, UEL, USDFLD subroutines and UEXTERNALDB utility subroutine to

AC C

implement the proposed model, where nonlocal averaging is performed on the derivatives of the equivalent strains with respect to the strain components during each load increment which are included in the nonlocal consistent tangent stiffness. Finally, damage evolution behaviors and load responses for two [0°/90°]4s T700/8911 composite specimens with different central hole sizes are studied by nonlocal FEA. Numerical results demonstrate the ability of the proposed model to address localized deformation, strain softening and strength prediction of composites effectively. Keywords A. Laminate, Polymer-matrix composites (PMCs) ; B. Strength, Mechanical properties, Delamination; C. Finite element analysis (FEA) *Corresponding author. E-mail: [email protected] and [email protected] (P.F.Liu)

1

ACCEPTED MANUSCRIPT 1.

Introduction Currently, laminated composite structures are increasingly applied to aerospace and aircraft

fields due to high strength and stiffness, low density and excellent designability. Lightweight

RI PT

design of composites as a fundamental work relates physical and mechanical properties of materials to geometry structures. However, failure mechanisms and damage evolution behaviors

prediction of load-bearing abilities of composites.

SC

of composite structures are often very complicated, which pose a large challenge to the accurate

M AN U

Beyond the initial failure of composites after the stress levels exceed the material strengths, mechanical properties of composites which are represented by progressive failure and continuous stiffness degradation become highly nonlinear. It is important to understand how the damage evolves after the initial failure as most of the energy dissipation takes place there. In general,

TE D

failure modes of composites include fiber pullout and breakage or kinking, matrix cracking, fiber/matrix shear failure and delamination. Degradation mechanisms and damage evolution behaviors depend on the ply’s orientation, the stacking sequence, the geometry size as well as the

EP

load and environment.

AC C

Currently, progressive failure analysis (PFA) based on continuum damage mechanics using finite element analysis (FEA) has grown up to be a popular method to predict stiffness degradation behaviors and load-bearing abilities of composite structures [1-13]. Associative intralaminar damage and interlaminar delamination of notched composite laminates by FEA are studied, where an important mechanism is interaction between matrix cracking and delamination [14-19]. However, progressive failure of composites inevitably meets the strain localization problem for geometrically discontinuous composite structures with holes, notches and cut-outs [9,10].

2

ACCEPTED MANUSCRIPT Localization means smoothly varying deformation field suddenly gives rise to one with highly localized field where high-gradient strain field, severe damage and energy dissipation concentrate on the local region, leading to ill-posed boundary value problems with the loss of ellipticity of the

RI PT

governing equation for quasi-static problems. No matter by experimental or numerical approaches, more efforts are needed to address localized deformation and damage evolution behaviors of composites with instabilities. Geers et al. [20] revealed the process zone where damage is

SC

accumulated by monitoring the strain localization behaviors of composites. Ibnabdeljalil and

M AN U

Curtin [21] and Martínez-Esnaola et al. [22] studied the localized stress concentration and stress redistribution around the fiber breakage of composites by micromechanical damage modeling. Kelly and Hallström [23] found the dominating localized matrix cracking under transverse load induces the delamination by experimental and numerical methods. Hallett et al. [24] divided the

TE D

damage evolution process of [45°/90°/-45°/0°]4s composite laminates with a central hole under tension into four stages by some experimental approaches (such as X-ray and C-scanning): 1 isolated matrix cracking at the hole and free edge and 2 interconnected damage at the hole (inner

EP

delamination regions) and localized damage at the free edge from matrix cracking and 3 damage

AC C

across the width of the specimen in an “influence zone” of the hole and 4 final catastrophic failure. The second stage shows that the strain localization with softening behavior starts to appear after initial matrix cracking. Mechraoui et al. [25] showed by acoustic emission test that the localization signals come from different failure mechanisms of composites and the extent of localization is represented by a real wave velocity. Bažant and Jirásek [26] showed the model needs to be enriched so as to capture the real process more adequately if the deformation field is expected to have important components with wave lengths below the resolution level of material model, i.e.,

3

ACCEPTED MANUSCRIPT the minimum size of the region into which the strain can localize. This shows there exists an intrinsic physical characteristic length scale that transmits the deformation information from the continuum body to the microstructure when the localization occurs. Unfortunately, FEA based on

RI PT

continuum mechanics lacks a length scale to describe the deformation field of microstructure, leading to spurious mesh sensitivity. In fact, zero characteristic length scale is approached when the mesh is gradually refined within the localized region. In this case, the ellipticity of quasi-static

SC

governing equation transforms into a hyperbolic system even though the boundary conditions are

M AN U

still prescribed in an elliptical system. It is beyond doubt that the reproduction of the localized damage zone of composites in their entirety, e.g., the delamination cracks and various types of transverse macrocracks is an important work. For example, Tay et al. [27] studied the localized damage of composites under transverse loads by FEA. Flatscher et al. [28] explored the

TE D

intralaminar plasticity and localized damage mechanisms of open-hole composites by FEA. Moure et al. [29] investigated the damage localization around material defects in composites using a discrete damage model. Crivelli et al. [30] monitored the localization behaviors for composite

EP

panels with fatigue matrix cracks and delamination using acoustic emission. However, these work

AC C

did’nt attempt to introduce any length scale to solve the localization problem. In recent decades, efforts to introduce a length scale as the localization limiter to regularize

ill-posed boundary value problems have led to a resurgence of research. The length scale is generally considered as the mean spacing of neighbouring representative microstructures. A common approach to deal with spurious mesh dependency by FEA is to refer to the smeared crack band model proposed by Bažant and Oh [31] who imposed a minimum element size on the fracture toughness, which was later adopted by such as Lapczyk and Hurtado [9] and Maimí et al.

4

ACCEPTED MANUSCRIPT [10] in intralaminar damage models for composites as built-in ABAQUS modules. However, Jirásek and Zimmermann [32] showed FEA by the crack band model is still suffered from stress locking with spurious stress transfer across widely opening cracks, mesh-induced directional bias

RI PT

and possible instability at late stages of the loading process. By comparison, a fundamental and general way to deal with the localization problem is to introduce the nonlocal integral theory developed mainly by Eingern and Edelen [33] and Pijaudier-Cabot and Bažant [34], in which the

SC

nonlocal internal variable is assumed as a spatial convolution of the local corresponding internal

M AN U

variable and the kernel function with a length scale. Pijaudier-Cabot and Bažant [34] pointed out that a key idea of the nonlocal integral theory is to perform nonlocal averaging on those physically motivated internal variables that control strain softening and to treat the strain, stress and other variables as local definitions. Bažant and Jirásek [26] further pointed out the motivation of

TE D

nonlocal averaging is that the damage model and FEA with strain softening lead to spurious localization of damage into a zone of zero volume, which causes unobjective numerical solutions with a vanishing energy dissipation during structural failure. This showed nonlocal averaging

EP

plays an important role in solving the localization problem with strain softening. Although they

AC C

demonstrated only one-dimensional (1D) problem, Bažant and Jirásek [26] concluded objective numerical results can also be achieved in multiple spatial dimensions. Recently, Li et al. [35] proposed non-local fracture model for composite laminates and performed numerical calculations using Fast Fourier Transforms (FFTs). Patel and Gupta [36] proposed a nonlocal continuum damage model for composite laminated panels, but numerical implementation using FEA was not introduced. Forghani et al. [37] and Liu et al. [38] performed nonlocal FEA on the localized intralaminar damage and interlaminar delamination of composites using LS-DYNA and ABAQUS

5

ACCEPTED MANUSCRIPT respectively. Germain et al. [39], Patel and Gupta [40] and Wu et al. [41] used the implicit gradient theory as an approximation to the nonlocal integral theory to solve the localization problem of composites. It is worth pointing out nonlocal averaging by FEA should ultimately

RI PT

contributes to stiffness equations instead of only local internal variables in order to achieve quadratic convergence. However, there are few nonlocal finite element models including the nonlocal consistent tangent stiffness for composites.

SC

The purpose of this paper is to solve the localization problem of composites by introducing

M AN U

the nonlocal integral theory developed by Pijaudier-Cabot and Bažant [34] into the intralaminar damage model of composites. An important originality is to derive the nonlocal implicit stiffness equation including the nonlocal intralaminar consistent tangent stiffness and the cohesive interface stiffness based on the variational form of the stress boundary value problem with an interface

TE D

discontinuity. Two length scales are introduced into the damage evolution behaviors of the fiber and matrix, indicating long-range interactions among different material points. Different from the work of Forghani et al. [37] and Liu et al. [38] who performed nonlocal averaging on the local

EP

equivalent strains for the fiber and matrix, a remarkable feature in this research is that the

AC C

derivatives of the equivalent strains with respect to the strain components are taken as the nonlocal variables, which are included in the nonlocal consistent tangent stiffness that is assembled with direct integration contribution from neighbouring elements, leading to a larger band width in the stiffness matrix. Therefore, the advantage of proposed nonlocal finite element model over the work by Forghani et al. [37] and Liu et al. [38] is to guarantee quadratic convergence of nonlinear implicit solutions. The proposed nonlocal model is implemented by combining ABAQUS-UEL, UMAT and USDFLD subroutines and UEXTERNALDB utility subroutine. The update of

6

ACCEPTED MANUSCRIPT nonlocal data is semi-implicit because nonlocal averaging is performed during each load increment in UEXTERNALDB. 3D FEA on the intralaminar damage and interlaminar delamination of two [0°/90°]4s T700/8911 composite specimens with different hole diameters is

properties and strengths.

SC

2. Nonlocal progressive failure model of composite laminates

RI PT

performed using the developed numerical technique to predict the nonlocal progressive failure

2.1 Intralaminar damage model

M AN U

The intralaminar anisotropic damaged stress-strain relationship is written as

σ = Cd : ε

(1)

where C d (1 − d i (i = 1, 2, 3), C ) is the damaged anisotropic four-order elastic tensor. C is the

TE D

undamaged four-order elastic tensor. di (i = 1, 2, 3) represent the fiber breakage, matrix cracking and fiber/matrix interface debonding, respectively. The damage configuration in the fiber principle coordinate system proposed by Wei et al. [42] is given by

 (1 − d1 ) E1  C d = MCM T =  1 − d1 1 − d 2 v12 E2 1 − v12 v21  0 

  0 , (1 − d2 ) E2 0 (1 − d3 )(1 − v12v21 ) G12   1 − d1  0    (2)  0 , = 1 − d , M   2    (1 − v12v21 ) G12  1 − d3    1 − d1 1 − d 2 v21 E1

0

EP

1

v21 E1

AC C

 E1  v12 E2 C=  1 − v12 v21  0 1

E2 0

d 2t , ε 22 + ε 33 > 0 d1t , ε11 > 0 d3 = 1 − (1 − d1t )(1 − d1c )(1 − d 2t )(1 − d 2 c ) , d1 =  and d 2 =  d 2 c , ε 22 + ε 33 < 0 d1c , ε11 < 0 where E1 and E2 are the longitudinal and transverse Young’s moduli, and v12 and v21 are the Poisson’s ratios, and G12 is the shear modulus. All damage variables change from zero (no damage) to one (complete damage). The matrix C d is symmetric due to v21 E1 = v12 E2 . d1t ,

7

ACCEPTED MANUSCRIPT d1c , d 2t and d 2c are damage variables for the fiber tension, fiber compression, matrix tension and matrix compression, respectively.

ε ij (i , j = 1, 2, 3) are the strain components.

Hashin [43] proposed the stress-based failure criteria for the fiber and matrix which can

RI PT

identify detailed failure modes under tensile and compressive states. Huang and Lee [44] showed by experiments that the strain may be a more favorable approach to assess the failure of composites because the strain is more continuous and smoother than the stress during the damage

SC

evolution of composites. Later, Linde et al. [45] proposed the strain-based failure criteria, which

M AN U

are yet limited to plane cases. Recently, Forghani et al. [37] and Liu et al. [38] proposed simple expressions for the equivalent strains as the failure criteria of fiber and matrix. However, tensile and compressive failure cannot be differentiated. In this research, we adopt the modified 3D strain-based Hashin failure criteria for the fiber and matrix proposed by Huang and Lee [44], but

TE D

the forms are changed to those proposed by Linde et al. [45] who introduced the equivalent strains

ε 1eq and ε 2eq for the damage evolution laws. The failure criteria for the fiber and matrix under tensile and compressive states are written as

ε11 > 0

EP

(a) Fiber tensile failure

2

2

ε =

(ε11 ) + (ε12 ) 2

AC C

eq 1

2

0  ε10t  2  ε 1t  0  0  + ( ε13 )  0  ≥ ε1t  ε12   ε13 

(3)

ε11 < 0

(b) Fiber compressive failure

ε1eq = ε11 ≥ ε10c

(c) Matrix tensile failure

(4)

ε 22 + ε 33 > 0 2

ε = eq 2

(ε 22 + ε 33 )

2

2

2

0 0  ε 20t   2 ε 22ε 33  2  ε 2t  2  ε 2t  +  0   ε 23 − 2 E2 E3  + ( ε12 )  0  + ( ε13 )  0  ≥ ε 20t (5) G23  ε12   ε 23     ε13 

(d) Matrix compressive failure in the through-thickness direction

8

ε 22 + ε 33 < 0

ACCEPTED MANUSCRIPT

 2 E2 E3  ε  ε  ε 23 − 2 ε 22ε 33  + ε 20c  120  + ε 20c  130  +  2 G23   ε12   ε13  (ε 230 ) 

2

≥ ε 20c

(6)

ε 10t = X T / E1 and ε 10c = X C / E1 are the fiber tensile and compressive strains,

respectively.

RI PT

where

ε

2

0 2c

ε 20t = YT / E2 and ε 20c = YC / E2 are the matrix tensile and compressive strains,

respectively. X T , X C , YT

and YC are the tensile and compressive strengths along and

perpendicular to the fiber direction, respectively.

ε120 = S12 / ( 2G12 ) , ε130 = S13 / ( 2G13 ) and

SC

ε 2eq =

2  E ε 0  2  0  ε 22 E2 + ε 33 E3  2 2c − 1 + ε 2 c  ( ε 22 + ε 33 )   + 0  0  2G12ε12    2G12ε 12 

M AN U

ε 230 = S23 / ( 2G23 ) are the shear strain strengths in the 1-2, 1-3 and 2-3 plane, respectively. S12 , S13 and S 23 are the shear strengths. E1 and E2 are the elastic moduli. G12 , G23 and G13 are the shear moduli.

Similar to the work of Lapczyk and Hurtado [9], Maimí et al. [10], Forghani et al. [37] and

TE D

Wei et al. [42], the damage variables for fiber breakage and matrix cracking are given by

EP

 ε1ft ,1c ( κ1 − ε10t ,1c ) d1t ,1c = , fiber tensile or compressive failure, κ1 ( ε1ft ,1c − ε10t ,1c )   ε 2ft ,2c (κ 2 − ε 20t ,2c )  , matrix tensile or compressive failure d 2t ,2 c = κ 2 ( ε 2ft ,2 c − ε 20t ,2 c ) 

(6)

AC C

where the damage load functions F1,2 = ε1,2 − κ1,2 and the Kuhn-Tucker loading-unloading conditions variables.

eq

F1,2 ≤ 0, κ&1,2 ≥ 0, κ&1,2 F1,2 = 0 are introduced. κ1 and κ 2 are the internal

ε1ft ,1c and ε 2ft ,2c are the final fiber breakage and matrix cracking strains which are

determined by experiments.

ε1ft ,1c > ε10t ,1c and ε 2ft ,2 c > ε 20t ,2c are required.

2.2 Bilinear cohesive models for predicting delamination Camanho et al. [46] and Turon et al. [47] proposed bilinear cohesive models. The single-mode traction-displacement jump relationships are given by

9

ACCEPTED MANUSCRIPT Ti = K ∆i , ∆imax ≤ ∆i0  s 0 max f Ti = (1 − di ) K ∆i , ∆i ≤ ∆i < ∆i ,  ∆imax ≥ ∆i f Ti = 0,

and

∆ if ( ∆ imax − ∆ i0 ) ∆

max i

∆ = [[u]]

(∆

f i

−∆

0 i

)

(0 ≤ d

s i

)

≤ 1, i = 1, 2,3) are damage variables.

u

(7)

is the displacement

c is the displacement jump. ∆3f = 2G1c / σ max , ∆1,2f = 2G2,3 / τ max .

σ max and

RI PT

where d = s i

(

 ∆imax = max ∆imax , ∆i , i = 1, 2, mode II and III   max max max  ∆3 = max ( ∆3 , ∆3 ) with ∆3 ≥ 0, mode I

τ max are the maximum tractions in the normal and tangential directions, respectively.

(

)

SC

∆10 = ∆20 = τ max / K , ∆30 = σ max / K . Gic i = 1, 2,3 are mode-I, II and III interlaminar fracture toughnesses, respectively. ∆imax ( i = 1, 2, 3) are the maximum displacement jumps in the load

M AN U

history indicating the irreversibility. In order to avoid the interpenetration of the completely failed interface under compression, the penalty stiffness T3 = K ∆3 ( ∆3 ≤ 0 ) is introduced in the normal direction, where the penalty stiffness K = 1E6 N/mm3 was suggested by Camanho et al. [46] and Turon et al. [47].

TE D

The mixed-mode traction-displacement jump relationships are written as

EP

Ti = K ∆i , ∆imax ≤ ∆i0 ,  s 0 max f Ti = (1 − d ) K ∆i , ∆i ≤ ∆i < ∆i ,  ∆imax ≥ ∆i f Ti = 0,

where the scalar mixed-mode damage variable d

∆mf ( ∆mmax − ∆m0 )

, ∆m = ∆12 + ∆22 + ∆3

AC C

ds =

∆

max m

(∆

f m

−∆

0 m

)

2

(8)

s

(0 ≤ d

s

≤ 1) is written as

, ∆mmax = max ( ∆mmax , ∆m ) ,

η  0 0 2  2  1+ β 2   c c c  β ∆ , > 0 ∆3 ∆1   G + G − G 3 ( 2 1 )  1 + β 2   , ∆3 > 0 2 2 1 0   ( ∆10 ) + ( β∆30 )    ∆m0 =  , ∆mf =  K ∆ m    2 2 0 2 0 2 ∆3 ≤ 0 ∆3 ≤ 0  ( ∆10 ) + ( ∆20 ) ,  ( ∆1 ) + ( ∆2 ) , 

where ⋅

denotes the MacAuley bracket. β =

(∆

2 1

(9)

+ ∆22 ) / ∆32 ( ∆3 > 0 ) is the mode mixity ratio

and η is the parameter defined in the B-K law [48]. The cohesive model introduces a length scale parameter that is generally called cohesive

10

ACCEPTED MANUSCRIPT zone length due to cohesive softening behavior, which is defined as the distance from the crack tip to the position where the maximum cohesive traction is attained [49,50]. Turon et al. [49] suggested three elements in the cohesive zone are sufficient to predict the delamination growth,

RI PT

which is adopted in this research. 2.3 Nonlocal integral theory

Strain localization as a real physical phenomenon appears widely in composites. Forghani et

SC

al. [37] showed the microscopic crack propagation of composites is essentially a structural

M AN U

problem instead of a material problem, which means the stress and strain at a material point are nonlocal. In order to solve the localization problem, a length scale is required to transmit the deformation information between the localized region and the entire continuum solid. Eringen [33] and Pijaudier-Cabot and Bažant [34] developed the nonlocal integral theory which assumes the

TE D

motion of one point depends not only on this point itself, but also on all points around this point. This actually introduces the nonlocal concept with long-range interactions among different material points. Pijaudier-Cabot and Bažant [34] showed only those internal variables that control

EP

the strain softening can be taken as the localization limiter. Forghani et al. [37] further showed the

AC C

damage variable is not a good choice for nonlocal averaging due to locking problem. Because nonlocal averaging directly contributes to the consistent tangent stiffness in this research, the derivatives of the equivalent strains

ε1eq and ε 2eq for the fiber and matrix in Eqs.(3)-(6) with

respect to the strain ε are used for nonlocal averaging, which are included in the nonlocal tangent stiffness as shown in Appendix-A eq ∂ε 1,2 (X p )

∂ε

=

1

α ( X p ) ∫Ω

α0 ( X p , X q )

eq ∂ε1,2 (Xq )

∂ε

dX q

(10)

where X q is the neighbouring material point around the point X p . The weight function 11

ACCEPTED MANUSCRIPT α 0 ( X p , X q ) and the normalization factor α ( X p ) are given by  X −X q − p α0 ( X p , X q ) = exp 32 3 2  2l ( 2π ) l  1

⋅

  , α ( X p ) = α 0 ( X p , X q ) dX q ∫Ω  

(11)

denotes the norm. Here, the length scale l is l( f ) and l( m) for the fiber and

matrix respectively. The nonlocal variables eq ∂ε 1,2 (X p)

eq ∂ε 1,2 (X p )

∂ε

at l( f ) → 0 and l(m) → 0 , respectively.

are reduced to the local variables

SC

∂ε

RI PT

where

2

Bažant and Jirásek [26] showed there is a subtle difference between the notions of the

M AN U

intrinsic material length and the length scale l . The latter is a special case of the former because nonlocal averaging is only one of the possible enrichments that introduce a length parameter into the continuum theory. This demonstrates the length scale l is introduced from not only physical but also mathematical considerations. Bažant and Jirásek [26] also showed the incremental

TE D

nonlocal variables are separated into the average part and interactive part and the latter represents the statistical continuum smearing of discrete cracks. This indicates the motivation of nonlocal

EP

averaging due to microcrack interactions contains the smeared crack band concept [31]. A controversial issue using the nonlocal integral theory is the treatment of the boundary. The

AC C

reason is that the regions which contribute to the nonlocal average near the boundary are smaller than those far from the boundary. Thus, damage evolution near the boundary is delayed and these regions are stronger than the internal parts of the body. Jirásek et al. [51] showed the kernel function

( ∫ α ( X , X ) dX ) / α ( X ) = 1 Ω

0

p

q

q

p

in Eq.(10) is a feasible strategy to eliminate the boundary

effect during nonlocal averaging. 2.4 Unified finite element model of nonlocal intralaminar damage and interlaminar delamination of composite laminates 12

ACCEPTED MANUSCRIPT The virtual work equation for the body Ω containing the discontinuous interface Γ d is given by the following weak form

Ω+ U Ω−

∇ X δ u : σ dV + ∫ δ ∆ ⋅ T dA = ∫ δ u ⋅ t b dA Γd

where δ u is the virtual displacement.

(12)

∂Ωt

ub is the Dirichlet boundary condition and tb is the

RI PT

∫

Neumann boundary condition. It is shown the virtual work done by external load is divided into the virtual strain energy (including the elastic deformation energy and damage dissipation energy)

SC

and the virtual cohesive interface energy.

implicit stiffness equation is written as nbulk + ncoh

A ( K

e =1

where

A

L

+ K c )  ∆u = Fext %

M AN U

After performing finite element discretization and linearization on Eq.(12), the nonlocal

denotes the assembly of the bulk and interface element stiffnesses.

(13)

nbulk and ncoh

TE D

represent the number of bulk elements and cohesive elements. ∆u is the increment of the node

%

displacement. The derivations of the nonlocal element stiffness K L and node force FL for

EP

bulk elements see Appendix-A. The derivations of the stiffness cohesive elements see Appendix-B. The external node force

K c and node force Fc for

Fext is taken as the incremental form

AC C

according to the Newton-Raphson algorithm. Numerical calculation ends until the residual node force Fext − FL − Fc reduces to the tolerance R . The implicit FEA often meets the numerical convergence problem due to softening behavior.

In this research, artificial viscous effect is introduced to improve the convergence [52]

 Fext − FL − Fc − Fv = R ,  *  Fv = cM v , v = ∆u / ∆t %

(14)

where M is the mass matrix calculated with unity density, v is the node velocity, ∆u is the *

%

13

ACCEPTED MANUSCRIPT node displacement increment, Fv is the artificial viscous force, c is a constant damping factor and ∆ t is the time increment during nonlinear iterations. It is noted introduced viscous forces should be large enough to regularize cohesive softening behaviors but small enough not to affect

3.

RI PT

numerical accuracy.

Numerical results and discussion

SC

In this research, FEA is used to implement the proposed nonlocal model. Although

M AN U

Pijaudier-Cabot and Bažant [34] and Jirásek et al. [51] demonstrated the role of the length scale by nonlocal averaging for 1D and 2D localization problems, little work is done for 3D problems. Nonlocal averaging contributes to the consistent tangent stiffness according to Appendix–A, in sharp contrast with that performed only on the equivalent strains ε1eq and ε 2eq in Eq.(10) and the

TE D

damage variables d1 and d 2 in Eq.(6) according to Forghani et al. [37] and Liu et al. [38]. Thus, numerical implementation in this work can be truly called as the nonlocal FEA. Fig.1 shows the flow chart of nonlocal FEA, in which ABAQUS-UEL (User element subroutine), UMAT (User

EP

material subroutine), USDFLD (User defined field subroutine) and UEXTERNALDB (User

AC C

external data utility subroutine) are used. The nonlocal intralaminar damage model is implemented by combining ABAQUS-UMAT, USDFLD and UEXTERNALDB, and the cohesive model for delamination is implemented by ABAQUS-UEL, where zero-thickness 3D eight-node cohesive elements are established. After the initial failure of fiber or matrix elements are detected, the nonlocal bulk element stiffness K L and cohesive element stiffness K c are assembled and the node forces FL and Fc are calculated at each iteration according to Eq.(13). Nonlocal averaging includes two key points: 1 the coordinate and volume at each integration point are

14

ACCEPTED MANUSCRIPT needed in ABAQUS-UEXTERNALDB, which are initially obtained by ABAQUS-USDFLD subroutine and assigned to those field variables stored in the COMMON BLOCK data module, which are shared by both ABAQUS-UEXTERNALDB and ABAQUS-USDFLD and 2 the

RI PT

derivatives of the local equivalent strains with respect to the strain components are also placed in a defined COMMON BLOCK data module shared by both ABAQUS-UMAT and ABAQUSUEXTERNALDB. Nonlocal averaging is performed on the derivatives of the equivalent strains

SC

with respect to the strain components by ABAQUS-UEXTERNALDB, which are returned to

M AN U

ABAQUS-UMAT by COMMON BLOCK. It is noted the update of nonlocal data is semi-implicit because nonlocal averaging is performed during each load increment instead of each iteration. After performing nonlinear calculations on the stiffness equation, all variables including the strains, equivalent strains, damage variables and stresses become nonlocal.

TE D

Because nonlocal FEA requires data exchanges among different integration points, parallel calculations are needed to improve computational efficiency. Besides, numerical technique by dynamically allocating and releasing arrays is also developed in UEXTERNALDB in order to

EP

increase computational capacity. Jirásek et al. [51] pointed out the convergence and equilibrium

AC C

iteration in the global level are usually more regular for the nonlocal model than for the local one. The reason is that the transition of strain fields from a diffuse damage mode to a localized one is relatively smooth, similar to the conclusion by Huang and Lee [44]. Therefore, required iterations that achieve equilibrium are reduced, which partly compensates increased computational cost for nonlocal search for neighbouring integration points. The geometry structure of two [0°/90°]4s composite specimens is shown in Fig.2. The thickness of each ply is 0.1125mm. Liu et al. [38] studied the effect of different layup patterns on

15

ACCEPTED MANUSCRIPT the damage evolution behaviors of composites by nonlocal FEA. Here, we further study the effects of two hole diameters: 5mm and 10mm. Fig.3 shows two symmetric mesh models including 1080 and 1440 bulk elements and 540 and 720 cohesive interface elements, respectively. Material

RI PT

parameters for T700/8911 composites are listed in Table 1 [53]. Turon et al. [47] obtained

η = 2.284 for AS4/PEEK composite by experiments. As η in Eq.(9) changes from 1.5 to 2.5 for T700/8911 composite, delamination results are not affected distinctly and η =2.0 is adopted in

SC

this research.

M AN U

The length scale controls the spread of the kernel function, but cannot be directly measured. In general, these parameters are explained as the mean half spacing between the representative microstructures [26]. However, the fiber/matrix microstructures of composites often take on random arrangements. Bažant and Jirásek [26] pointed out there are two length scale parameters:

TE D

one of them characterizing the length of process zone and the other describing the width of process zone, and any identification procedure through experimental data extracts from only a single length scale parameter when the fixed ratio of these two length scale parameters are

EP

assumed. However, this assumption can hardly be considered as a general rule and a favorable

AC C

method is to extract two length scale parameters. In the proposed model, two length scales l(f) and l(m) for longitudinal fiber breakage and transverse matrix cracking are introduced. Currently, there are some advanced experimental approaches for the identification of the length scale of composite structure, such as the strain field method proposed by Geers et al. [54], the elastic wave dispersion method proposed by Dontsov et al. [55] and the micro-cantilever test method proposed by Dehrouyeh-Semnani and Nikkhah-Bahrami [56]. Yet, it is almost impossible to determine the length scales rigorously for complex structures [26]. It is more believed the length scale is

16

ACCEPTED MANUSCRIPT introduced from not only the physical concept but also the numerical consideration. In this research, l(f) and l(m) are assumed as adjustable parameters comparable to the minimum element size suggested by Rousseller et al. [57]. Ahad et al. [58] further showed the length scale must be

RI PT

greater than (or at least equal to) the minimum element size in the regions where the damage may occur so that the regularization by nonlocal averaging is variable. In the work of Liu et al. [38], same values for l(f)=l(m)=1mm and 2mm are adopted for simplification. Yet, Forghani et al. [37]

SC

showed employing a relatively small transverse nonlocal radius by orthotropic averaging is helpful

M AN U

to capture a narrow matrix cracking pattern, and took l(f)=2mm and l(m) =0.5mm, 1mm and 2mm respectively. In this research, we attempt to study the effect of different length scales l(f), l(m)=0.5mm, 1mm and 2mm on the load-displacement curves for two mesh models. Furthermore, the length scales are approximately determined by the load stress-strain curves, similar to the work

TE D

by Ahad et al. [58] and Korsunsky and Kim [59] who determined the length scale parameter for ductile metal materials by notched tensile tests.

The tensile stress-strain curves for two composite specimens using different length scales

EP

and two mesh models are shown in Fig.4. Table 2 lists the tensile strengths and fracture strains by

AC C

FEA. After the initial failure of fiber or matrix, the stiffness of structure degrades which leads to nonlinear load response. Intralaminar matrix cracking and delamination are main failure modes before fiber breakage, where matrix cracking induces a little delamination growth near the hole. From the nonlocal FEA, the decreasing parts of tensile curves tend to be consistent for two mesh models, which show objective description of the damage evolution and energy dissipation after performing nonlocal averaging. By comparison, the local FEA that lacks a length scale in the governing equation, leads to some mesh dependence. When the length scales l(f) and l(m) vary from

17

ACCEPTED MANUSCRIPT 0.5mm to 2.0mm, the load curves do not change too much. By referring to the suggestion by Rousseller et al. [57] that the length scale should be larger than the minimum element side length in Fig.3, we take l(f) =l(m)=0.5mm in the following analysis.

ε 1eq and ε 2eq and the damage variables d1

RI PT

The contours of the nonlocal equivalent strains

and d 2 of the fiber and matrix for the specimen with a 5mm-diameter hole by nonlocal FEA are shown in Figs.5-8. The corresponding results for the specimen with a 10mm-diameter hole by

SC

nonlocal FEA are shown in Figs.9-12. The predicted damage evolution behaviors are close using

M AN U

two mesh models for two specimens. The initial failure for both the fiber and matrix appears near the hole. For the specimen with a 5mm-diameter hole, the matrix starts to fail at the strain 0.342% and the equivalent matrix strain 0.85%, and the fiber fails initially at the strain 0.702% and equivalent fiber strain 2.654%. For the specimen with a 10mm-diameter hole, the matrix starts to

TE D

fail at the strain 0.229% and the equivalent matrix strain 0.637%, and the fiber fails initially at the strain 0.472% and equivalent fiber strain 2.367%. The load-bearing ability decreases when the diameter of hole adds. After the initial failure, matrix cracks propagate to the whole specimen, but

EP

fiber breakage propagates at about 30º angle to the loading direction, which are consistent with the

AC C

results obtained by Lee et al. [18] who adopted the Puck’s failure criteria [60] and damage mechanics method to study local failure of the fiber and matrix of [(0º/90º)6]s composite laminates under a particular stress combination. Chang and Chang [1] proposed a damage model to study in-plane damage mechanisms of

[0º/±45º/90º]2s composites although the delamination was not considered. They showed the damage initiates through matrix cracking in the 90º plies before propagating across the 45º plies, and fiber breakage follows as the cracks reach the 0º plies. Kortschot and Beaumont [61]

18

ACCEPTED MANUSCRIPT performed experimental research on the observed subcritical notch tip cracking patterns in [00n / 90n0 ]ns carbon fiber/epoxy laminates and established a qualitative relationship between the

terminal damage and notched strength. They divided the visible failure modes into the transverse

RI PT

cracks in the 90º plies, the axial splitting in the 0º plies and the triangular delamination zone at the 0º/90º interface. Kongshavn and Poursartip [62] developed a strain softening approach to predict the failure of [-45º/90º/45º/0º/]2s notched composites and showed the damage zone propagates

SC

from the notch-tip in a stable, self-similar manner which includes matrix cracks and delamination

M AN U

in the off-axis plies. Green et al. [63] performed a series of experiments on open-hole

[450m / 900m / −450m / 00m ]ns composite laminates by employing three different scaling methods: in-plane scaling, sublaminate-scaling and ply-scaling. Failure modes for the sublaminate-scaling specimens (m=1, n=1,2,4,8) are the fiber pull-out or brittle failure, and load-drop is caused by

TE D

fiber breakage in the 0º plies. By comparison, ply-scaled specimens (m=1,2,4,8, n=1) all fail in the delamination mode, and load-drop is caused by delamination at the -45º/0º interface prior to fiber breakage in the 0º plies. Hallett et al. [24] divided the damage evolution of sublaminate-scaling

EP

[45°/90°/-45°/0°]4s laminates into four stages mentioned above. For the sublaminate-scaling

AC C

[0°/90°]4s laminates in this research, transverse matrix cracking appears at the 90º plies and then induces a little delamination near the hole. With progressive delamination, a little splitting (longitudinal matrix cracking) appears in the 0º plies, accompanied by a little fiber breakage. Suddenly much fiber breakage appears indicating the final collapse of structures. Either fiber breakage or matrix cracking occurs across the laminate where the weakest sections of these plies are located. Figs.13-16 compare the results by local and nonlocal FEA for two specimens using fine mesh model. After performing nonlocal averaging, the damage evolution of fiber and matrix

19

ACCEPTED MANUSCRIPT is alleviated to some extent, which is represented by the reduction in the intensity of damage within the localized area. In addition, selected nonlocal variables, i.e., equivalent strains for the fiber and matrix which control the damage evolution and strain softening are also demonstrated to

RI PT

be suitable for nonlocal averaging. The effect of different viscous constants c in Eq.(14) on the load responses for two specimens is shown in Fig.17. In general, c=(1E-3)-(5E-3) is acceptable for achieving excellent convergence

SC

and numerical accuracy. Figs.18-19 show the delamination results at the final failure strain using

M AN U

different cohesive strengths for two specimens. Here, delamination along the load direction is detected by the dislocation of pair nodes at the same position. The delamination initiation and growth concentrate mainly on the area near the hole. Besides, predicted delamination using the cohesive model becomes more difficult when the normal and tangential cohesive strengths

τ max increase from 5MPa to 20MPa. When σ max and τ max exceed 20MPa, no delamination

TE D

and

σ max

appears. In addition, it is found the load-displacement curves are the same for the same specimen and mesh model when

σ max and τ max exceed 1MPa. This shows the proposed nonlocal model

σ max and τ max increase, similar to the conclusions obtained by Turon et al. [49]

AC C

arises when

EP

can eliminate the mesh dependency without considering delamination. The convergence difficulty

and Liu et al. [38]. Compared with the specimen with a 10mm-diameter hole, the specimen with a 5mm-diameter hole leads to more distinct delamination because of more severe deformation mismatch between the 0° plies and 90° plies near the hole, similar to the conclusion obtained by Forghani et al. [37]. In general, fiber breakage plays a dominating role in governing the tensile strength of open-hole composite laminates although matrix cracking and delamination decrease the stiffness and integrity of structures largely [64].

20

ACCEPTED MANUSCRIPT 4.

Concluding remarks Strain localization is commonly considered as a challenging problem in composite structures.

Local FEA cannot predict the damage evolution behaviors and load-bearing abilities of composites

RI PT

accurately due to the lack of a length scale. This paper introduces the nonlocal integral theory developed by Pijaudier-Cabot and Bažant [34] into the damage model to solve the localization problem of composites. Based on the variational form of the stress boundary value problem with

SC

an interface discontinuity, this paper originally develops a finite element model on the nonlocal

M AN U

intralaminar damage and interlaminar delamination of composite laminates, where two length scales for the fiber and matrix are introduced. It should be emphasized nonlocal averaging on the derivatives of equivalent strains with respect to the strain components for the fiber and matrix in this research ultimately contributes to the implicit nonlocal consistent tangent stiffness. The

TE D

proposed model is implemented by combining ABAQUS-UMAT, UEL, USDFLD subroutines and UEXTERNALD utility subroutine. For two [0°/90°]4s T700/8911 composite specimens with 5mmand 10mm-diameter central holes respectively, effects of different length scales on the load

EP

responses are comparatively studied. By comparing the local and nonlocal results in terms of

AC C

load-displacement curves and damage evolution behaviors, introduced length scales play an important role in predicting localized damage evolution behaviors of composites accurately after eliminating spurious mesh dependence. The developed finite element model and numerical codes by ABAQUS as robust and reliable computational technique will be further applied to aerospace structures with a variety of geometry discontinuities including holes, notches and cut-outs. However, perhaps the localization problem for composite structures that is expected to fundamentally solved using FEA would have to resort to high-performance computation.

21

ACCEPTED MANUSCRIPT Acknowledgements Dr. Pengfei Liu would sincerely like to thank the support of the National Natural Science Funding of China (No.51375435), the Natural Science Funding of Zhejiang Province of China (No.

RI PT

LY13E050002) and Aerospace Science and Technology Innovation Funding.

Appendix-A: Nonlocal consistent tangent stiffness tensor K L for bulk elements

SC

Jirásek and Patzák [65] derived an isotropic consistent tangent stiffness for damaged concrete

M AN U

materials. In the following, we further derive an anisotropic nonlocal consistent tangent stiffness for damaged composites. The discrete incremental form of stress tensor ∆σ is given by

∑ ∆σ ( X

p

)=

p

∑ C ∆ε = ∑ C d

p

p

(

)

d ( p)

B ( X p ) ∆u %

(A.1)

where C d ( p ) = C d C ,1 − di ( p ) (i = 1, 2) is defined and the strain matrix B ( X p ) is simplified

TE D

as B p . ω p is the weight of the integration point X p . The element stiffness tensor K L is derived by

 T  ∂C d B C d B + BT B ∆u  d V  Ω UΩ ∂∆u % e =1  % nbulk   ∂C d = A  ∑ ω p BTp C d ( p ) B p + ∑ ω p BTp B p ∆u  ∂∆u % e =1  p p % nbulk

A∫

+

−

AC C

EP

KL =

(A.2)

where it is noted the last term in Eq.(A.2) includes two parts: fiber and matrix. From Eq.(2), differentiating the damaged elastic tensor C d with respect to the node

displacement increment ∆u leads to

%

22

ACCEPTED MANUSCRIPT

where

∂κ1

∂ε

eq 1

∂κ 2

= 1 and

∂ε

eq 2

M AN U

SC

RI PT

 eq eq  ∂C d = ∂C d ∂d1 ∂κ1 ∂ε1 + ∂C d ∂d 2 ∂κ 2 ∂ε 2 ,  ∂∆u ∂d1 ∂κ1 ∂ε1eq ∂∆u ∂d 2 ∂κ 2 ∂ε 2eq ∂∆u % %  % ε1ft ε1ct  ∂d1t  ∂κ = ε f − ε 0 κ 2 , Fiber tensile failure ε11 > 0, 1t 1t ( 1 )  1 f  ∂d ε1c ε1cc 1c = f , Fiber compressive failure ε11 < 0,  2 0  ∂κ1 ε1c − ε1c (κ1 )  ε 2ft ε 2ct  ∂d 2t = , Matrix tensile failure ε 22 + ε 33 > 0,  2 0 f ε ε ∂ κ − κ ( ) 2t 2t  2 2  f ε 2cc  ∂d 2 c = ε 2 c , Matrix compressive failure ε 22 + ε 33 < 0,  ∂κ 2 ε 2fc − ε 20c ( κ 2 )2   ∂ε1eq  ∂ε1eq  ∂ε q  ∂ε1eq  = ω α = ω α ∑q q pq  ∂ε  ∂∆u ∑q q pq  ∂ε  Bq ,  ∂∆u  q  q  % %  eq eq ∂ε q  ∂ε 2eq   ∂ε 2 = ∑ ω α  ∂ε 2  = ∑ ωqα pq   Bq q pq  ∂∆u ∂ε  q ∂∆u ∂ε  q q q    % % ∂κ1

= 1 represent loading, and

∂ε

eq 1

the strain

TE D

unloading. The nonlocal derivatives of the local equivalent strains

= 0 and

(A.3) ∂κ 2

∂ε 2eq

ε 1eq and ε 2eq with respect to

ε are derived by Eq.(10)

EP

eq eq eq  ∂ε1eq ∂ε1 ∂ε 2 ∂ε 2 ( X p ) = ∑ ω qα pq ( X q ), ( X p ) = ∑ ωqα pq ( X q ),  ∂ε ∂ε ∂ε q q  ∂ε  α0 X p − X q α pq =  ∑r wrα 0 X p − X r 

(

)

AC C

(

to the components of strain tensor (a) Fiber tensile failure

∂ε ∂ε11

(A.4)

)

Based on Eqs.(3-6), the differentiation of the equivalent strains

eq 1

= 0 represent

ε 1eq and ε 2eq with respect

ε leads to the following expressions

ε11 > 0

2 2 0 0   2 2  ε 1t  2  ε 1t  = ε11 ( ε11 ) + ( ε12 )  0  + ( ε13 )  0     ε12   ε13  

23

−

1 2

(A.5)

ACCEPTED MANUSCRIPT ∂ε ∂ε12

eq 1

2 2 0 0  ε   2 2  ε 1t  2  ε 1t   ε ε ε = ε 12  + + ( ) ( ) ( )  0     0  11 12 13  ε    ε12   ε13   0 1t 0 12

2

 ε   ∂ε 2 2ε  2ε  = ε13   ( ε11 ) + ( ε12 )   + ( ε13 )    ∂ε13  ε   ε   ε   2

0 1t 0 12

(b) Fiber compressive failure

2

0 1t 0 13

ε11 < 0

∂ε1eq =1 ∂ε11

−

1 2

(A.7)

SC

ε 22 + ε 33 > 0

2 2 2 0 0   ε 20t   2 ε 22ε 33  ε   2 2  ε 2t  2  ε 2t   ε ε ε ε ε = ε12  + + − E E + + ( 13 )  0    ( 22 33 )  0   23 2 3  ( 12 )  0  2 ε ε ε G      23     12   ε 13   23 0 2t 0 12

2

M AN U

∂ε ∂ε12

(A.6)

(A.8)

(c) Matrix tensile failure

eq 2

2

1 2

RI PT

0 1t 0 13

eq 1

−

2

TE D

0 2t 0 13

(A.10)

EP

 E E  ε 0 2   ∂ε 2eq  = ε 22 + 1 − 2 2 3  20t   ε 33  × ∂ε 22  G23  ε 23       

AC C

2 2 2 0 0    ε 20t   2 ε 22ε 33  2 2  ε 2t  2  ε 2t  ( ε 22 + ε 33 ) +  0   ε 23 − 2 E2 E3  + ( ε12 )  0  + ( ε13 )  0   G23   ε12   ε 23     ε 13  

∂ε ∂ε 23 eq 2

−

1 2

(A.11)

2 2 0 2 0  ε    ε 20t   2 ε 22ε 33  2 2  ε 2t  2  ε 2t   = ε 23  + + − E E + + ε ε ε ε ε ( 13 )  0    ( 22 33 )  0   23 2 3  ( 12 )  0  G232  ε12   ε    ε 23     ε13   0 2t 0 23

2

(A.12)

24

1 2

(A.9)

2 2 2 0 0  ε    ε 20t   2 ε 22ε 33  ∂ε 2 2  ε 2t  2  ε 2t  =ε13   ( ε 22 + ε 33 ) +  0   ε 23 − 2 E2 E3  + ( ε12 )  0  + ( ε13 )  0   ∂ε13 G23  ε    ε 23     ε12   ε13   eq 2

−

−

1 2

−

1 2

ACCEPTED MANUSCRIPT  E E  ε 0 2   ∂ε 2eq  = ε 33 + 1 − 2 2 3  20t   ε 22  × ∂ε 33  G23  ε 23        −

1 2

RI PT

2 2 2 0 0    ε 20t   2 ε 22ε 33  2 2  ε 2t  2  ε 2t  ( ε 22 + ε 33 ) +  0   ε 23 − 2 E2 E3  + ( ε12 )  0  + ( ε13 )  0   G23   ε12   ε 23     ε 13  

(A.13)

(d) Matrix compressive failure in the through-thickness direction

SC

2ε 0 ε ∂ε 2eq = 2c 12 2 ∂ε12 (ε120 )

ε 22 + ε 33 < 0

M AN U

0 2ε 2c ε13 ∂ε 2eq = 2 0 ∂ε13 (ε13 )

(A.14)

(A.15)

ε 0 E (ε E + ε E ) ε 0 ε E E  E ε 0  ∂ε 2eq = 2c 2 22 2 2 33 3 − 2c 33 2 23 +  2 2c0  − 1 0 ∂ε 22 2 ( G12ε12 ) ( G23ε 230 )  2G12ε12 

(A.16)

0 2ε 2c ε 23 ∂ε 2eq = 2 ∂ε 23 (ε 230 )

(A.17)

TE D

2

ε 0 E (ε E + ε E ) ε 0 ε E E  E ε 0  ∂ε 2eq = 2c 3 22 2 2 33 3 − 2c 22 2 23 +  2 2c0  − 1 0 ∂ε 33 2 ( G12ε12 ) ( G23ε 230 )  2G12ε12 

EP

2

(A.18)

AC C

Substituting Eq.(A.3) into Eq.(A.2) leads to the weighted nonlocal consistent tangent

stiffness K L

  ∂C d   ∂d1   ∂κ1   ∂ε 1eq  T T K L = A ∑ ω p B p C d ( p ) B p + ∑ ω pωqα pq B p B p      eq    Bq ∆u % e =1  p p,q  ∂d1  p  ∂κ1  p  ∂ε1  p  ∂ε q nbulk

 ∂C d   ∂d 2   ∂κ 2   ∂ε 2eq   + ∑ ω pωqα pq B B p      eq    Bq ∆u  % p ,q  ∂d 2  p  ∂κ 2  p  ∂ε 2  p  ∂ε  q  T p

(A.19) where it is noted K L is not symmetric. Last two terms represent interactions among different 25

ACCEPTED MANUSCRIPT integration points. Here, nonlocal averaging is only performed over those integration points whose distance is smaller than the nonlocal radius R, preserving the sparsity of stiffness matrix.

RI PT

Appendix-B: Stiffness matrix K c and internal force Fc for cohesive elements Segurado and LLorca [66] proposed a finite element formulation for implementing the cohesive model by adopting the middle-plane numerical technique, as shown in Fig.20. Similar to

(ξ ,η , ζ ) and the global coordinate

SC

standard isoparametric elements, the local coordinate system

M AN U

system (1, 2,3) are defined. The displacement jump in the local coordinate system [[ u]] for the cohesive element is calculated as

[[u]] = ϕ T [[u(ξ ,η )]], ϕ = ϕ (n, t%1 , t%2 ) where

(B.1)

ϕ (n, t%1 , t%2 ) is a 3×3 transformation matrix of the local coordinate system to the global one

TE D

and three orthotropic directions are given by

∂x R ∂ξ , ∂x R ∂ξ

EP

 ∂x R ∂x R  ×   ∂ξ ∂η   n= , t%1 = ∂x R ∂x R × ∂ξ ∂η

t%2 = n × t%1

where the coordinate of middle-plane point x

1 H ( ξ ,η ) ( I12×12 2

AC C

xR =

R

(B.2)

is written as

I12×12 ) ( x + u ) %

(B.3)

where x and u are the 24×1 node coordinate matrix and the 24×1 node displacement matrix

%

for 3D eight-node elements respectively. H (ξ ,η ) is a 3×12 matrix for 3D eight-node elements including the shape function and I12×12 is the identity matrix. n is the normal direction and t%1 and t%2 are tangential directions. The relative displacements at the point

(ξ ,η ) between element sides are interpolated as a

function of relative displacements between paired nodes 26

ACCEPTED MANUSCRIPT [[u1 (ξ ,η )]]    [[u (ξ ,η )]] = [[u2 (ξ ,η )]] = H ( ξ ,η ) [[uN ]] = H (ξ ,η ) φ u = Lu % % [[u (ξ ,η )]] 3  

L is a 3×24 matrix and −1 ≤ ξ , η ≤ 1 are the local coordinates of elements.

(

φ = I12×12

I12×12 ) is a 12×24 matrix.

The node force vector

RI PT

where

(B.4)

Fc ( 24 × 1 matrix ) and the stiffness matrix K c ( 24 × 24 matrix ) are

written as

i

SC

Fc = ∫ ∫ LT ϕ T TdA = ∫ ∫ M T TdA = ∑∑ ωiω j M T T J , j

∂Fc ∂T ∂T ∂[[u]] ∂T ∂[[u]] dA = ∫ ∫ M T dA = ∫ ∫ M T dA = ∫ ∫ MT ∂u ∂u ∂[[u]] ∂u ∂[[u]] ∂u % % % % ∂T T T M dA== ∑∑ ωiω j M KM J = ∫∫M ∂[[u]] i j

(B.5)

M AN U

Kc =

where J is the Jacobian of isoparametric transformation and

ω is the weight. M is a 3×24

shape function matrix for the cohesive element and K is a 3 × 3 cohesive stiffness matrix in

References

TE D

which diagonal values are K in Eq.(7).

EP

[1] Chang FK, Chang KY. A progressive damage model for laminated composites containing stress concentrations. J Compos Mater 1987; 21(9): 834-855. [2] Allix O, Ladeveze P, Gilletta D, Ohayon R. A damage prediction method for composite

AC C

structures. Int J Numer Meth Eng 1989; 27(2): 271-283. [3] Tan SC. A progressive failure model for composite laminates containing openings. J Compos Mater 1991; 25(5): 556-577. [4] Ladeveze P, LeDantec E. Damage modeling of the elementary ply for laminated composites. Compos Sci Technol 1992; 43(3): 257-267. [5] Matzenmiller A, Lubliner J, Taylor RL. A constitutive model for anisotropic damage in fiber-composites. Mech Mater 1995; 20(2):125-152. [6] Lubineau G, Ladeveze P. Construction of a micromechanics-based intralaminar mesomodel, and illustrations in ABAQUS/Standard. Comput Mater Sci 2008; 43(1): 137-145. [7] Reddy YSN, Moorthy CMD, Reddy JN. Non-linear progressive failure analysis of laminated composite plates. Int J Non-Linear Mech 1995; 30(5): 629-649. 27

ACCEPTED MANUSCRIPT [8] Barbero EJ, Vivo LD. A constitutive model for elastic damage in fiber-reinforced PMC Laminae. Int J Damage Mech 2001; 10(1): 73-93. [9] Lapczyk I, Hurtado JA. Progressive damage modeling in fiber-reinforced materials. Compos Part A 2007; 38(11): 2333-2341. [10] Maimí P, Camanho PP, Mayugo JA. A continuum damage model for composite laminates: Part I-constitutive model. Mech Mater 2007; 39(10): 897-908.

RI PT

[11] Tay TE, Liu G, Tan VBC, Sun XS, Pham DC. Progressive failure analysis of composites. J Compos Mater 2008; 42(18): 1921-1966.

[12] Liu PF, Zheng JY. Progressive failure analysis of carbon fiber/epoxy composite laminates using continuum damage mechanics. Mater Sci Eng A 2008; 485(1-2): 711-717.

SC

[13] Van Der Meer FP, Sluys LJ. Continuum Models for the Analysis of Progressive Failure in Composite Laminates. J Compos Mater 2009; 43(20): 2131-2156.

[14] Goyal VK, Jaunky NR, Johnson ER, Ambur DR. Intralaminar and interlaminar progressive

M AN U

failure analysis of composite panels with circular cutouts. Compos Struct 2004; 64(1): 91-105. [15] Kashtalyan M, Soutis C. Analysis of composite laminates with intra- and interlaminar damage. Prog Aero Sci 2005; 41(2):152-173.

[16] Hallett SR, Jiang WG, Khan B, Wisnom MR. Modelling the interaction between matrix cracks and delamination damage in scaled quasi-isotropic specimens. Compos Sci Technol 2008; 68(1): 80-89.

TE D

[17] Chen JF, Morozov EV, Shankar K. Simulating progressive failure of composite laminates including in-ply and delamination damage effects. Compos Part A 2014; 61(7): 185-200. [18] Lee CS, Kim JH, Kim SK, Ryu DM, Lee JM. Initial and progressive failure analyses for composite laminates using Puck failure criterion and damage-coupled finite element method.

EP

Compos Struct 2015;121: 406-419.

[19] Mohammadi B, Olia H, Hosseini-Toudeshky H. Intra and damage analysis of laminated composites using coupled continuum damage mechanics with cohesive interface layer. Compos

AC C

Struct 2015; 20: 519-530.

[20] Geers MGD, Peijs T, Brekelmans WAM, de Borst R. Experimental monitoring of strain localization and failure behaviour of composite materials. Compos Sci Technol 1996; 56(11):1283-1290.

[21] Ibnabdeljalil M, Curtin WA. Strength and reliability of fiber-reinforced composites: Localized load-sharing and associated size effects. Int J Solids Struct 1997; 34(21): 2649-2668. [22] Martínez-Esnaola JM, Martín-Meizoso A, Daniel AM, Sánchez JM, Elizalde MR, Puente I, Fuentes M. Localized stress redistribution after fibre fracture in brittle matrix composites. Compos Part A 1997; 28(4): 347-353 [23] Kelly G, Hallström S. Strength and failure mechanisms of composite laminates subject to localised transverse loading. Compos Struct 2005; 69(3): 301-314. 28

ACCEPTED MANUSCRIPT [24] Hallett SR, Green BG, Jiang WG, Wisnom MR. An experimental and numerical investigation into the damage mechanisms in notched composites. Compos Part A 2009; 40(5): 613-624. [25] Mechraoui SE, Laksimi A, Benmedakhene S. Reliability of damage mechanism localisation by acoustic emission on glass/epoxy composite material plate. Compos Struct 2012; 94(5): 1483-1494. [26] Bažant ZP, Jirásek M. Nonlocal integral formulations of plasticity and damage: survey of

RI PT

progress. ASCE J Eng Mech 2002; 128(11): 1119-1149.

[27] Tay TE, Lam KY, Cen Z. Analysis of composite structures with distributed and localized damage by the finite-element method. Compos Struct 1997; 37(2): 135-143.

[28] Flatscher TH, Wolfahrt M, Pinter G, Pettermann HE. Simulations and experiments of open

SC

hole tension tests-assessment of intra-ply plasticity, damage, and localization. Compos Sci Technol 2012; 72(10): 1090-1095.

[29] Moure MM, Sanchez-Saez S, Barbero E, Barbero EJ. Analysis of damage localization in

M AN U

composite laminates using a discrete damage model. Compos Part B 2014;66: 224-232. [30] Crivelli D, Guagliano M, Eaton M, Pearson M, Al-Jumaili S, Holford K, Pulli R. Localisation and identification of fatigue matrix cracking and delamination in a carbon fibre panel by acoustic emission. Compos Part B 2015; 74(1): 1-12.

[31] Bažant ZP, Oh BH. Crack band theory for fracture of concrete. Mater Struct 1983; 16(3):155-177.

TE D

[32] Jirásek M, Zimmermann T. Rotating crack model with transition to scalar damage. ASCE J Eng Mech 1998; 124(3): 277-284.

[33] Eringen AC, Edelen BL. Nonlocal elasticity. Int J Eng Sci 1972; 10(3): 233-248.

1512-1533.

EP

[34] Pijaudier-Cabot GU, Bažant ZP. Nonlocal damage theory. ASCE J Eng Mech 1987; 113(10):

[35] Li J, Meng SH, Tian XX, Song F, Jiang CP. A non-local fracture model for composite laminates and numerical simulations by using the FFT method. Compos Part B 2012; 43(3):

AC C

961-971.

[36] Patel BP, Gupta AK. An investigation on nonlocal continuum damage models for composite laminated panels. Compos Part B 2014; 60: 485-494. [37] Forghani A, Zobeiry N, Poursartip A, Vaziri R. A structural modelling framework for prediction of damage development and failure of composite laminates. J Compos Mater 2013; 47: 20-21. [38] Liu PF, Yang YH, Gu ZP, Zheng JY. Finite element analysis of progressive failure and strain localization of composite laminates. Appl Compos Mater (In Press) [39] Germain N, Besson J, Feyel F. Composite layered materials: Anisotropic nonlocal damage models. Comput Meth Appl Mech Eng 2007; 196(41-44): 4272-4282. [40] Patel BP, Gupta AK. An investigation on nonlocal continuum damage models for composite 29

ACCEPTED MANUSCRIPT laminated panels.Compos Part B 2014; 60: 485-494. [41] Wu L, Sket F, Molina-Aldareguia JM, Makradi A, Adam L, Doghri I, Noels L.A study of composite laminates failure using an anisotropic gradient-enhanced damage mean-field homogenization model. Compos Struct 2015; 126: 246-264. [42] Wei XD, de Vaucorbeil A, Tran P, Espinosa HD. A new rate-dependent unidirectional

61(6): 1305-1318.

RI PT

composite model-application to panels subjected to underwater blast. J Mech Phys Solids 2013;

[43] Hashin Z. Failure criteria for unidirectional fiber composites. J Appl Mech 1980; 47: 329-334.

[44] Huang CH, Lee YJ. Experiments and simulation of the static contact crush of composite

SC

laminated plates. Compos Struct 2003; 61(3): 265-270.

[45] Linde P, Pleitner J, DeBoer H, Carmone C. Modelling and simulation of fiber metal laminates. In: ABAQUS Users' Conference, 2004.

M AN U

[46] Camanho PP, Davila CG, de Moura MF. Numerical simulation of mixed-mode progressive delamination in the composite materials. J Compos Mater 2003; 37(16): 1415-1438. [47] Turon A, Camanho PP, Costaa J, Dávila CG. A damage model for the simulation of delamination in advanced composites under variable-mode loading. Mech Mater 2006; 38(11): 1072-1089.

[48] Benzeggagh ML, Kenane M. Measurement of mixed-mode delamination fracture toughness

TE D

of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos Sci Technol 2006; 56(4): 439-449.

[49] Turon A, Dávila CG, Camanho PP, Costa J. An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Eng Fract Mech 2007; 74(10):

EP

1665-1682.

[50] Harper PW, Hallett SR. Cohesive zone length in numerical simulations of composite delamination. Eng Fract Mech 2008; 75(16): 4774-4792.

AC C

[51] Jirásek M, Rolshoven S and Grassl P. Size effect on fracture energy induced by non-locality. Int J Numer Analy Meth Geomech 2004; 28(7): 653-670. [52] Abaqus Version 6.12 Documentation-Abaqus Analysis Users Manual. [53] Liu PF, Chu JK, Liu YL, Zheng JY. A study on the failure mechanisms of carbon fiber/epoxy composite laminates using acoustic emission. Mater Des 2012; 37: 228-235. [54] Geers MGD, de Borst R, Brekelmans WAM, Peerlings RHJ. Validation and internal length scale determination for a gradient damage model: application to short glass-fibre-reinforced polypropylene. Int J Solids Struct 1999; 36(17): 2557-2583. [55] Dontsov EV, Tokmashev RD, Guzina BB. A physical perspective of the length scales in gradient elasticity through the prism of wave dispersion. Int J Solids Struct 2013; 50(22-23): 3674-3684. 30

ACCEPTED MANUSCRIPT [56] Dehrouyeh-Semnani AM, Nikkhah-Bahrami M.A discussion on evaluation of material length scale parameter based on micro-cantilever test. Compos Struct 2015; 122: 425-429. [57] Rousselier G. In: Nemat-Nasser (Ed.), Finite deformation constitutive relations and ductile fracture. North-Holland, Amsterdam, pp. 201-217. [58] Ahad FR, Enakoutsa K, Solanki KN, Bammann DJ. Nonlocal modeling in high-velocity impact failure of 6061-T6 aluminum. Int J Plast 2014; 55(4): 108-132.

RI PT

[59] Korsunsky AM, Kim K. Determination of essential work of necking and tearing from a single tensile test. Int J Fract 2005; 132: 137-144.

[60] Puck A, Schürmann H. Failure analysis of FRP laminates by means of physically based phenomenological models. Compos Sci Technol 1998; 58(7): 1045-1067.

SC

[61] Kortschot MT, Beaumont PWR. Damage mechanics of composite materials: I-measurements of damage and strength. Compos Sci Technol 1990; 39(4): 289-301.

[62] Kongshavn I, Poursartip A. Experimental investigation of a strain-softening approach to

M AN U

predicting failure in notched fibre-reinforced composite laminates. Compos Sci Technol 1999; 59(1): 29-40.

[63] Green BG, Wisnom MR, Hallett SR. An experimental investigation into the tensile strength scaling of notched composites. Compos Part A 2007; 38(3): 867-878.

[64] Talreja R. Stiffness properties of composite laminates with matrix cracking and interior delamination. Eng Fract Mech 1986; 25(5-6): 751-762.

TE D

[65] Jirásek M, Patzák B. Consistent tangent stiffness for nonlocal damage models. Comput Struct 2002; 80(14-15): 1279-1293.

[66] Segurado J, LLorca J. A new three-dimensional interface finite element to simulate fracture in

EP

composites. Int J Solids Struct 2004; 41(11-12): 2977-2993.

Figure captions

AC C

Fig.1 Flow chart of the nonlocal FEA of progressive failure of composite laminates Fig.2 Geometry structures and sizes for two composite specimens with a central hole Fig.3 Two mesh models for two composite specimens Fig.4 Stress-strain curves for [0º/90º]4s specimens with a (a) 5mm-diameter hole and (b) a 10mm-diameter hole using different length scales and mesh models Fig.5 Contours of the equivalent strain

ε1eq and damage variable d1 of the fiber for the

[0º/90º]4s laminate with a 5mm-diameter hole using fine mesh model at the strains (a) 0.702%, (b) 1.195% and (c) 1.557%, respectively Fig.6 Contours of the equivalent strain

ε 2eq and damage variable d 2 of the matrix for the

[0º/90º]4s laminate with a 5mm-diameter hole using fine mesh model at the strains (a) 0.342%, (b) 1.195% and (c) 1.557%, respectively 31

ACCEPTED MANUSCRIPT Fig.7 Contours of the equivalent strain

ε1eq and damage variable d1 of the fiber for the

[0º/90º]4s laminate with a 5mm-diameter hole using coarse mesh model at the strains (a) 0.794%, (b) 1.244% and (c) 1.557%, respectively Fig.8 Contours of the equivalent strain

ε 2eq and damage variable d 2 of the matrix for the

[0º/90º]4s laminate with a 5mm-diameter hole using coarse mesh model at the strains (a) 0.342%, (b) 1.244% and (c) 1.557%, respectively

ε1eq and damage variable d1 of the fiber for the

RI PT

Fig.9 Contours of the equivalent strain

[0º/90º]4s laminate with a 10mm-diameter hole using fine mesh model at the strains (a) 0.472%, (b) 1.092% and (c) 1.445%, respectively Fig.10 Contours of the equivalent strain

ε 2eq and damage variable d 2 of the matrix for the

SC

[0º/90º]4s laminate with a 10mm-diameter hole using fine mesh model at the strains (a) 0.229%, (b) 1.092% and (c) 1.445%, respectively Fig.11 Contours of the equivalent strain

ε1eq and damage variable d1 of the fiber for the

(b) 1.053% and (c) 1.445%, respectively

M AN U

[0º/90º]4s laminate with a 10mm-diameter hole using coarse mesh model at the strains (a) 0.635%,

Fig.12 Contours of the equivalent strain

ε 2eq and damage variable d 2 of the matrix for the

[0º/90º]4s laminate with a 10mm-diameter hole using coarse mesh model at the strains (a) 0.310%, (b) 1.053% and (c) 1.445%, respectively

Fig.13 Contours of the equivalent strain

ε1eq and damage variable d1 of the fiber for the

TE D

[0º/90º]4s laminate with a 5mm-diameter hole using fine mesh model at the strain 1.557% Fig.14 Contour of the equivalent strain

ε 2eq and damage variable d 2 of the matrix for the

[0º/90º]4s laminate with a 5mm-diameter hole using fine mesh model at the strain 1.557% Fig.15 Contours of the equivalent strain

ε1eq and damage variable d1 of the fiber for the

EP

[0º/90º]4s laminate with a 10mm-diameter hole using fine mesh model at the strain 1.445% Fig.16 Contours of the equivalent strain

ε 2eq and damage variable d 2 of the matrix for the

AC C

[0º/90º]4s laminate with a 10mm-diameter hole using fine mesh model at the strain 1.445% Fig.17 Stress-strain curves of two [0º/90º]4s specimens using different viscous constants using fine mesh model

Fig.18 Delamination of the [0º/90º]4s laminate with a 5mm-diameter hole using (1) coarse mesh model and (2) fine mesh model at the strain 1.557% and cohesive strengths (a)

σ max = τ max =5MPa (b) σ max = τ max =10MPa and (c) σ max = τ max =20MPa, respectively Fig.19 Delamination of the [0º/90º]4s laminate with a 10mm-diameter hole using (1) coarse mesh model and (2) fine mesh model at the strain 1.445% and cohesive strengths (a)

σ max = τ max =5MPa (b) σ max = τ max =10MPa and (c) σ max = τ max =20MPa, respectively Fig.20 Eight-node cohesive element using middle-plane interpolation technique

32

ACCEPTED MANUSCRIPT Tables Ply longitudinal modulus

E1

135 GPa

Ply transverse modulus

E2

11.41 GPa

Out-of-plane modulus

E3

11.41 GPa

Inplane shear modulus

G12

7.92 GPa

Poisson’s ratio

G13

7.92 GPa

G23

3.79 GPa

v12

0.33

v13

0.33

RI PT

Out-of-plane shear modulus

0.49

XT

2600 MPa

Longitudinal compressive strength

XC

1422 MPa

Transverse tensile strength

YT

603 MPa

SC

v23 Longitudinal tensile strength

YC

241 MPa

S12

94 Mpa

Shear strength

S13

94 Mpa

Shear strength

S23

40 Mpa

M AN U

Transverse compressive strength Shear strength

Interlaminar mode-I fracture toughness

c

0.252 N/mm

G1 c

Interlaminar mode-II/III fracture toughness G2

= G3c 0.665 N/mm

Initial fiber tensile breakage strain

ε10t

= X T / E1

0.019

Initial matrix tensile cracking strain

ε 20t

= YT / E2

0.005

ε1ft

0.021

ε 2ft

0.035

TE D

Critical fiber tensile breakage strain

Critical matrix tensile cracking strain

AC C

EP

Table 1 Mechanical parameters of T700/8911 composites [53]

33

ACCEPTED MANUSCRIPT

Tensile strengths (MPa)

AC C

Composite laminate with a 10mm-diameter hole

EP

TE D

(1)Test results (2)Coarse mesh and nonlocal model at l(f) = l(m) = 0.5mm (3)Coarse mesh and nonlocal model at l(f) = l(m) = 1.0mm (4)Coarse mesh and nonlocal model at l(f) = l(m) = 2.0mm (5)Coarse mesh and nonlocal model at l(f) = 2.0mm and l(m) = 0.5mm (6)Coarse mesh and nonlocal model at l(f) = 2.0mm and l(m) = 1.0mm (7) Fine mesh and nonlocal model at l(f) = l(m) = 0.5mm (8) Fine mesh and nonlocal model at l(f) = l(m) = 1.0mm (9) Fine mesh and nonlocal model at l(f) = l(m) = 2.0mm (10) Fine mesh and nonlocal model at l(f) = 2.0mm and l(m) = 0.5mm (11) Fine mesh and nonlocal model at l(f) = 2.0mm and l(m) = 1.0mm (12) Coarse mesh and local model (13) Fine mesh and local model

Fracture strains (%) 1.560 1.575

795

1.585

800

1.595

791 793

SC

781

RI PT

809 790

1.577 1.580

1.557

785

1.566

784

1.564

M AN U

(1) Test results (2)Coarse mesh and nonlocal model at l(f) = l(m) = 0.5mm (3)Coarse mesh and nonlocal model at l(f) =l(m) = 1.0mm (4)Coarse mesh and nonlocal model at l(f) = l(m) = 2.0mm (5)Coarse mesh and nonlocal model at l(f) = 2.0mm and l(m) = 0.5mm (6)Coarse mesh and nonlocal model at l(f) = 2.0mm and l(m) = 1.0mm Composite laminate with a (7)Fine mesh and nonlocal model at 5mm-diameter l(f) = l(m) = 0.5mm (8)Fine mesh and nonlocal model at hole l(f) = l(m) = 1.0mm (9)Fine mesh and nonlocal model at l(f) =l(m) = 2.0mm F (10) Fine mesh and nonlocal model at l(f) = 2.0mm and l(m) = 0.5mm (11) Fine mesh and nonlocal model at l(f) = 2.0mm and l(m) = 1.0mm (12) Coarse mesh and local model (13) Fine mesh and local model

782

1.559

784

1.564

756 776

1.503 1.549

714

1.524

723

1.480

726

1.489

729

1.493

722

1.479

725

1.486

706

1.445

705

1.443

713

1.460

710

1.453

711

1.456

677 711

1.379 1.456

Table 2 Tensile strengths and fracture strains using nonlocal and local FEA for two mesh models

34

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT