An investigation of matrix cracking damage evolution in composite laminates – Development of an advanced numerical tool

Composite Structures 108 (2014) 937–950

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

An investigation of matrix cracking damage evolution in composite laminates – Development of an advanced numerical tool G. Sadeghi a, H. Hosseini-Toudeshky a,⇑, B. Mohammadi b a b

Department of Aerospace Engineering, Amirkabir University of Technology, 424 Hafez Avenue, Tehran, Iran School of Mechanical Engineering, Iran University of Science and Technology, Narmak, Tehran, Iran

a r t i c l e

i n f o

Article history: Available online 16 October 2013 Keywords: Matrix cracking Progressive damage evolution General symmetric laminates Stress transfer Finite element method

a b s t r a c t In this study, progressive matrix cracking damage analyses are performed for general symmetric laminate with any number of plies under in-plane tensile/shear loading condition. To predict initiation and evolution of matrix cracking in each layer of such laminate, an advanced numerical tool is developed here. The tool consists of, (i) a micromechanic-based damage model to find stiffness reduction of the laminate and damage parameters due to the presence of matrix cracking, (ii) an energy-based evolution law to predict initiation and progression of the damage and (iii) a finite element (FE) basement to implement the constitutive damage law analyzing laminates under complex boundary conditions. The stress transfer method is used to find the displacement and stress fields in a unit cell confined between two consecutive matrix cracks. Then, the degradation of elastic constants and damage parameters are calculated for each damaged layer of the laminate. New formulation is also developed to implement the evolution law at integration points to eliminate the length dependency of the energy-based damage criterion. The crack density is used as the only state variable which controls the damage state. An eight node element is developed using the full layerwise plate theory and the damage constitutive law is implemented in the formulation of a user defined element routine in Abaqus commercial software. The developed procedure is validated favorably with published experimental and numerical results for matrix cracking damage mechanism in general symmetric laminates. Concurrent progressive damage analyses in different plies of various laminates are also investigated and some results which are rarely available in the literature are presented in this paper. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction The application of composite materials has been promoted from secondary to primary structures in different industries particularly in aerospace industry. Despite of the desirable mechanical properties, composites are susceptible to damage from a number of sources, both during manufacturing and in service. The integrity of structural components during their service life is an important design requirement when composite materials are used in loadbearing conditions. But damage occurrence can never be entirely avoided, so composite structures should be designed for safe operation despite the presence of damages. The first form of damage in composite laminates in the early stages of loading is usually matrix cracks, which are intra-laminar cracks that traverse the thickness of the ply and run parallel to the fibers in that ply. Although this type of damage is not the final failure mode of the laminates, it can degrade their thermo-mechanical properties and trigger the other damage mechanisms such as induced delamination. Usually matrix cracking happens when the ⇑ Corresponding author. Tel.: +98 21 66405032; fax: +98 21 66959020. E-mail address: [email protected] (H. Hosseini-Toudeshky). 0263-8223/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruct.2013.10.007

tensile load is applied in a direction different from the fiber direction, such as a cross-ply laminate in which the load is perpendicular to 90° plies and in off-axis plies in balanced/unbalanced laminates. In complex in-plane loading condition, the cracks may appear in every laminate ply. Most of the previously performed studies in progressive matrix cracking analyses have been restricted to special layup configurations in which the damage evolution in a laminate layer was investigated [1,2]. Also in most of the previously performed researches, the loading condition was limited to the axial tensile load applied to un-notched laminates that experiencing uniform in-plane stress field. Fan and Zhang [3] developed the equivalent constraint model (ECM) based on the shear-lag stress transfer method to predict the progressive matrix cracking in 90° plies of [±h/90n]s laminates [4]. The ECM was used to calculate the stiffness reduction of the laminate and the energy-based criterion employed to control the damage initiation and growth. McCartney [5] developed a methodology similar to ECM using a more accurate stress transfer method. He used a homogenization technique to replace the thermomechanical effects of fully-developed discrete cracks in one ply by the effective homogenized properties such that the homogenized laminate contains exactly the same properties as the

938

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

laminate with discrete ply cracks. This procedure enables the prediction of progressive matrix cracking formation, based on the energy methods, in any ply of a general symmetric laminate. Hosseini-Toudeshky et al. [6–8] developed a damage model on ply level and implements that in a FE based code to find the structural response of the general symmetric laminates having arbitrary layup configuration. The major drawback of such damage analyses procedure is that the used rectangular unit cell in the micromechanic model cannot be quite compatible with the finite element discretization around a discontinuity such as a notch or a circular hole. For [/m/un]s laminates containing matrix cracks, Yakozeki and Aoki [9] developed a 2D shear-lag model in an oblique coordinate system along the matrix cracks to find the effective in-plane thermo-elastic properties of the laminates. By fixing the crack density, k1, in one ply, the stiffness reduction of the laminate was calculated as a function of the crack density, k2, in the neighboring ply. No results were reported to predict progressive matrix cracking in both bunches of the laminate layers in this study. Using Yakozeki’s progressive damage concept and based on the shear-lag approach, Barbero et al. [10–13] developed a damage constitutive law to find the laminate stiffness reduction. Using FEM based formulation of a shell element; the model was implemented at integration points to predict the matrix crack evolution in general symmetric laminates. No results were reported for simultaneous progression of damage in multi plies of the laminate in this study. For cross ply laminates having any number of plies and arbitrary stacking sequence, the authors have been developed a model to perform progressive matrix crack analysis [14]. In the present study, a progressive damage analysis procedure is developed to predict matrix cracking in symmetric laminates with arbitrary stacking sequence. In this hypothesis each laminate ply is susceptible to damage and therefore matrix cracking evolution has to be considered for all layers. In this way, transverse and shear damage parameters are defined to consider the degradation of elastic constants of each ply. A micromechanical damage model based on the stress transfer method is used to find the stiffness degradation of the unit cell. The method may leads to more accurate results compared to other shear-lag approaches as it benefits from a three dimensional stress–displacement field of a cracked unit cell [15]. The unit cell which is defined between two consecutive matrix cracks is extracted for each ply with a proper coordinate system transformation. Classical lamination theory is used to find the ply properties in the laminate and transverse and shear damage parameters are calculated for each cracked ply. The finite fracture mechanics and energy based failure criterion are used to predict the initiation and evolution of the damage in all laminate layers. The strain energy, which is calculated using the degraded properties of the laminate at integration points, is decomposed to the related terms according to the normal and shear loading condition. The effect of mixed mode condition on matrix cracking is also considered. To implement the model at an integration point, a fictious plate is considered there and the strain energy release rate is evaluated for such plate as a function of the unit cell degraded properties and crack density. User element definition capability of Abaqus commercial software is used to implement the damage constitutive law. The capability of the standard code to entitle Fortran routines helps us to construct and solve the governing differential equations arising from the micromechanic damage model. This is a key point in the performed progressive damage analysis procedure in the present study. The developed tool can be used to predict the progressive matrix damage of laminates with complex geometry and boundary conditions. It is employed to predict the progressive matrix crack damage for several cross-ply, balanced and general symmetric CFRP and GFRP laminates and the obtained results are compared with the available experimental and numerical results.

2. Micromechanic damage model The stress transfer model, similar to the other micromechanicbased models, is formulated based on a unit cell which is a representative continuum of a laminate limited between the two consecutive matrix cracks. The goal of such model is to predict stiffness deterioration of cracked laminate using stress–displacement fields of the unit cell. The stress–displacement field of such unit cell is calculated using the constitutive law, equilibrium equation, interface and edge boundary conditions and simplifying assumptions. The main assumptions are piecewise linear form of the out of plane perturbed shear stress components through the thickness and the generalized plain strain condition of cracked laminate due to the complexity of the problem. Also some of the boundary conditions are not met exactly and must be averaged through the thickness of layers. There is no restriction on the number and orientation of plies and by the refinement of layers, the accuracy of the solution can be improved. When the stress and displacement fields of the problem are obtained, the stiffness reduction of the unit cell can be calculated as a function of crack density. Fig. 1 shows the extracted unit cell from cracked laminate with N + 1 layer at the structural scale. The matrix cracks may evolve in every layer in accordance with the boundary and loading conditions, but the crack evolution is investigated in only one layer at each calculation step. Suppose that the jth layer is under consideration and damage evolution occurs in that ply only. For more clarity and better understanding of the modeling, the exaggerated cracked layer is shown in Fig. 1. Three coordinate systems are also introduced as depicted in this figure. The structure or global coordinate system (XYZ) which is a fixed system during the analysis is used to define the layers orientations in the laminate. The angles of the plies are measured with respect to the X direction of the global coordinate system. The unit cell coordinate system (xyz) is defined according to the extracted unit cell in which the x direction is normal to the crack face, y direction is in the laminate plane and along the crack direction and the z direction is normal to the laminate plane. It is clear that the unit cell coordinate system varies during the analysis and must be defined in proper rotation with respect to the global coordinate system for jth layer. The last coordinate system is the material coordinate system (123) which is defined for each layer, as direction 1 is along the fibers, direction 2 is in the layer plane transverse to the fibers and direction 3 is normal to the lamina plane. The unit cell is limited to the region in which jxj 6 a,jyj 6 b of the laminate and due to the symmetry condition, half of the laminate, 0 6 z 6 H is modeled only. The laminate has N + 1 perfectly bonded layers with arbitrary stacking sequence. The locations of the N interfaces are specified by z = zk , k = 1, . . . , N. The laminate mid-plane is specified by z = z0 = 0 and the external surface by z = zN+1 = H where 2H is the laminate total thickness. The thickness of the kth layer is denoted by hk = zk  zk1. The constitutive law for the kth layer (k = 1, 2, . . . , N+1) in the material coordinate system (123) is defined as: 8 9k e1 > > > > > > > > > e2 > > > > > > < > =

e3

> e4 > > > > > > > > > e5 > > > > > : > ;

e6

21 6 6 6 6 6 6 ¼6 6 6 6 6 4

E1

t12 E1 1 E2 ð1d2 Þ

symm:

t12 E1 t23 E2

0

0

0

0

0

0

1 E2

0

0

0

1 G23

0

0

1 G12

0 1 G12 ð1d6 Þ

3k 7 7 7 7 7 7 7 7 7 7 7 5

8 9k r1 > > > > > > > > > r2 > > > > > > < > =

r3

> r4 > > > > > > > > > r5 > > > > > : > ;

r6

þ

8 9k a1 > > > > > > > > > a2 ð1  da Þ > > > > > > > < =

a2

> > > > > > > > :

0 0 0

> > > > > > > > ;

DT

ð1Þ

By transverse isotropy hypothesis, seven thermo-mechanical constants are only needed to describe the stress–strain relations. The constants are longitudinal and transverse young’s modulus E1 and

939

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

Extracted Unit Cell Matrix Cracks

Uncracked Constraint Layers

Craked Layer

Symm. Plane

3

z

2b (N+1)th ply 2a th

. . .

j ply H st

2

. . .

1 ply

1

Symm. Plane Cracked Layer Fig. 1. Typical unit cell extracted from damaged laminate (the top figure is schematically shown the extraction of the unit cell and its orientation with respect to loading direction).

E2, in-plane shear modulus G12, poison’s ratios m12 and m23 and longitudinal and transverse thermal expansion coefficients, a1 and a2. The isotropy condition relates the out-of plane shear modulus, G23, to the transverse young modulus and m23 by G23 = E2/2(1 + m23). Such assumptions are used for virgin or un-cracked laminate material properties. The presence of the damage in a layer will change the state to an orthotropic behavior and such condition is dictated by the damage parameters. Three damage parameters, d2, d6 and da are introduced according to the matrix cracking effects on the degradation of individual material constants. The matrix cracking has no effect on the ply properties in fiber and thickness directions, and therefore no damage parameters are defined for related constants. But this damage type causes the degradation of thermo-mechanical constants in transverse direction, and the damage parameters are considered for related terms. These parameters are related to the crack density as an internal state variable which defines the damage state. They are all zero-valued for virgin material and by the damage evolution will grow to the maximum value of unity which defines the fully loss of the material stiffness. As an example the damage parameter, d2 has the scale of degradation of transverse elastic modulus by the evolution of matrix cracking. In compact notation in material coordinate system, Eq. (1) cab be written as,

fegkm ¼ ½Sk frgkm þ fagkm DT

ð2Þ

where, the index m refers to the material coordinate system, {e} = {e1 e2 e3 e4 e5 e6}T is the strain vector, {r} = {r1 r2 r3 r4 r5 r6}T is the stress vector, {a} = {a1a2(1  da)a2 0 0 0}T is the thermal expansion coefficients vector, [S] is the compliance matrix of kth layer in the material coordinate system and DT is the temperature change with respect to the free stress state. The superscript T denotes the transpose of the vector. Using transformation matrix in the unit cell coordinate system,

^ g k DT fegkUC ¼ ½Sk frgkUC þ fa

ð3Þ

where

fegkUC ¼ h exx

eyy ezz 2exz 2eyz 2exy ik

ð4Þ

frgkUC ¼ h rxx

ryy rzz rxz ryz rxy ik

ð5Þ

½Sk ¼ ½RTk ½Sk ½Rk

ð6Þ

^ gk ¼ ½RTk fagk fa

ð7Þ

The index UC refers to the unit cell coordinate system and the transformation matrix is [16],

2

c2 6 s2 6 6 6 0 ½Rk ¼ 6 6 0 6 6 4 0 cs

s2 c2

0 0 0 0

0

1 0

0

0

0

c

s

0

0

s

c

cs

0 0

0

0 0

3k 2cs 2cs 7 7 7 0 7 7 0 7 7 7 0 5

ð8Þ

c2  s2

hk Þ; s ¼ sinð^ hk Þ and ^ hk ¼ hk  hRj . hk is the angle bewhere c ¼ cosð^ tween X direction of the global system and the fiber direction of kth layer (measured in clockwise direction). hRj is the rotation angle of the unit cell for current jth cracked layer under investigation and can be calculated as hRj ¼ hj  90  SGNðhj Þ, in which SGN is the sign function. By this transformation, the constitutive law is formulated in the unit cell coordinate system in which cracked faces of jth layer are normal to the x direction. When the matrix cracks form in the planes normal to the x direction, away from the edge of laminate, it can be assumed that the transverse strains of all layers are the same. This is the generalized in-plane strain condition which releases the complexity of the

940

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

3D constitutive law. In this way, the stress field is independent from y direction and the displacement field is defined as,

u ¼ uðx; zÞ;

v ¼ f ðx; zÞ þ ec y;

w ¼ wðx; zÞ

ð9Þ

where u, v and w are deformation in the x, y, and z direction respectively, is the transverse strain of the cracked laminate and with equal value in all layers. By substitution of this definition with ekyy component in (3) yield to:

ek2 ¼ ekyy ¼

  @v k ^ 2 DT ¼ ec ¼ Sk12 rkxx þ Sk22 rkyy þ Sk23 rkzz þ Sk26 rkxy þ a @y ð10Þ

and the stress component in y direction is:

rkyy ¼

1  Sk22

ec  Sk12 rkxx  Sk23 rkzz  Sk26 rkxy þ a^ 2 DT

By eliminating of



ð11Þ

rkyy from (3) using Eq. (11) yields to,

 gk DT fegk ¼ ½Fk frgk þ fGgk ec þ fa

ð12Þ

where {e}k and {r}k are the strain and stress vectors in the unit cell coordinate system and the non-zero components of matrix [F] and  g are presented in Appendix A. Therefore, the vectors {G} and fa 3D stress–strain field is simplified to the generalized plane strain hypothesis, and the equilibrium equations are also simplified to,

8 k @ rk @ rxx > þ @zxz ¼ 0 > > < @x @ rkxy > @x

> > : @rkxz @x

þ

@ rkyz @z

¼0

þ

@ rkzz @z

¼0

ð13Þ

/0k ðaÞ ¼ 0;

k ¼ 1; 2; . . . ; N

ð17Þ

For cracked layer, the normal rxx and shear stress rxy components are zero and from Eqs. (B2) and (B3) in Appendix B,

/k ðxÞ  /k1 ðxÞ ¼ hk rk

ð18Þ

wk ðxÞ  wk1 ðxÞ ¼ hk sk

ð19Þ

rk and sk are normal and shear stress components of kth layer in undamaged laminate before occurrence of any matrix cracking. For un-cracked layers, the edge boundary conditions are not satisfied exactly because of dependency on the z (through the thickness) direction, so the average conditions are imposed on the x = a plane, namely,  k ðaÞ ¼ ec a u

ð20Þ

v k ðaÞ ¼ cc a þ ec y

ð21Þ

Using Eqs. (B21) and (B22) from Appendix B,

ec ¼

cc ¼

k    F k11  ~ ~ k ðaÞ þ F 15 w ~ k1 ðaÞ  w ~ k ðaÞ /k1 ðaÞ  / hk a hk a    k1 DT þ F k11 rk þ F k15 sk þ Gk1 ec þ a

ð22Þ

k    F k15  ~ ~ k ðaÞ þ F 55 w ~ k1 ðaÞ  w ~ k ðaÞ /k1 ðaÞ  / hk a hk a    k5 DT þ F k15 rk þ F k25 sk þ Gk5 ec þ a

ð23Þ Rx

The basic assumption in the stress transfer mechanics approach is that the out of plane shear stresses varies in linear form along the thickness direction and they are unknown functions of x along the length of unit cell. For the kth layer,

(

the x = a plane, thus from (14) considering the zero value for / and w on the mid-plane and free surface of the unit cell yields to,





rkxz ¼ /0k1 ðxÞ þ /0k ðxÞ  /0k1 ðxÞ zzhkk1 ; k ¼ 1; 2; . . . ; N þ 1   rkyz ¼ w0k1 ðxÞ þ w0k ðxÞ  w0k1 ðxÞ zzhkk1

ð14Þ

where / and w are unknown functions of x at the interface of layers and appear due to the presence of cracks which perturb the stress field. These functions vanish at mid-plane and free surface of the laminate. By obtaining of these functions, the stress and displacement field of the unit cell and the stiffness reduction of laminate due to the presence of damage can be calculated. In the above equations, the prime sign shows the derivative of functions with respect to x variable. Using equilibrium equations, definition of unknown functions in (14) and interface and edge boundary conditions, the coupled system of governing ordinary differential equation (ODE) is obtained (B25). The details of formulations leading to ODE system of the problem are outlined in Appendix B. Solution of such system reveals the stress and displacement field of the unit cell. It is worth to note that due to the symmetry condition with respect to the x = 0 plane, only 3N boundary conditions are required to find a unique solution for the system of differential equations (N is the number of interfaces between the layers). The edge conditions at x = a for the cracked layer are defined as:

rkxx ða; y; zÞ ¼ rkxz ða; y; zÞ ¼ rkxy ða; y; zÞ ¼ 0

ð15Þ

and for un-cracked layer(s),

rkxz ða; y; zÞ ¼ 0; uk ða; y; zÞ ¼ ec a; v k ða; y; zÞ ¼ cc a þ ec y

ð16Þ

in which ec, ec and cc are the effective applied axial, transverse and shear strain components. The shear stress, rxz is zero everywhere on

k= 1, . . . , N + 1 for un-cracked layers only and g~k ðxÞ ¼ 0 g k ðsÞds. But the parameters, ec and cc are not yet known and must be eliminated from these boundary conditions. Let l be the label of the nearest uncracked layer to the center plane of the laminate, it then follows by subtraction that the required boundary conditions for un-cracked layers derived from (22) and (23) are in the form of: l    F k11  ~ ~ k ðaÞ  F 11 / ~ l1 ðaÞ  / ~ l ðaÞ /k1 ðaÞ  / hk a hl a l    F k15  ~ ~ k ðaÞ  F 15 w ~ l1 ðaÞ  w ~ l ðaÞ ¼ 0 wk1 ðaÞ  w þ hk a hl a l    F k15  ~ ~ k ðaÞ  F 15 / ~ l1 ðaÞ  / ~ l ðaÞ /k1 ðaÞ  / hk a hl a l    F k55  ~ ~ k ðaÞ  F 55 w ~ l1 ðaÞ  w ~ l ðaÞ ¼ 0 wk1 ðaÞ  w þ hk a hl a

ð24Þ

ð25Þ

where k and l take values appropriate for un-cracked layers such that k – l. The conditions in (17)–(19) and (24) and (25) provide the 3N boundary conditions to ensure that the system of governing differential equation of the problem has a unique solution. 3. Mechanical properties degradation and damage parameters The solution of differential equations in (B25), using 3N boundary conditions in Eqs. (17)–(19) and (24) and (25), reveals the stress and displacement fields of the unit cell as a function of crack density. Such fields are used to find the degradation of mechanical properties of the cell due to the presence of cracks. Then the damage parameters are calculated according to the degradation of inplane thermo-mechanical constants. In this way, the in-plane stress–strain relations for the undamaged laminate is extracted from Eq. (12), and by performing some algebraic manipulation, the stress components for kth layer are obtained as,

941

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

rkxx ¼

h     i k k k k k k k   k1  F k15 a  k5 DT  2 F 55 e  G1 S66  G5 S16 e  F 15 c  F 55 a k  F 15 1

F k11 F k55

ð26Þ 2  3 8 9 1 > >  Gk1 Gk1 Sk66  Gk5 Sk16   > > Sk22 > > 6 7 k k k k > > > F 15 G5  F 55 G1 e þ 4   5e > > > > > k k k k k > > < = G G S  G S 5 1 16 5 11 1 2   3 rkyy ¼  2 > k k k > >  k  F 15 a  k5  G F a  >  > F k11 F k55  F k15 > > 6 1 55 1 7 > > > >  þ F k15 Gk1  F k11 Gk5 c þ 4  DT > 5 > > > > ^ k : ;  k1  F k11 a  k5  ak2 G5 F k15 a S 22

ð27Þ

rkxy ¼

1



F k11 F k55  F k15

h     i k k k k k k k   k1  F k11 a  k5 DT 2 F 15 e þ G1 S16  G5 S11 e þ F 11 c þ F 15 a

1 2a

Z

a

8 9 > =

r > : s



8 9k 8 9 8 9 0 1 > > > Nþ1 Nþ1 < rxx > = = = 1X 1X k B k C  ¼ hk ryy ¼ hk @½Q  e þ fRg DT A ¼ ½Q  e þ fRgDT > > > > H k¼1 ; H k¼1 > : ; : > ; : > ;

rxy

c

ð30Þ

cc ¼ s21 r þ s22 s þ s23 ec þ s24 DT

ð31Þ

r ¼ s31 r þ s32 s þ s33 ec þ s34 DT

ð32Þ

By solving (32) for



ec ¼ s11  s13

8 92 8 93 8 94 3 0> 0> 0> > > > > > > > > > > > > > 7 <1> <0> <0= = = 7 7 ; ; 7 > > > > > > 0 1 0 > > > > > > > 5 > > > > > : ; : ; : ; 0 0 1

28 9 r 1 >1 > > 6> > < = 6 0> s 6 ¼ ; 6  > e > > > 4> >0> > > > : c > ; : > ; DT 0 8 > > > <

9 > > > =

ð33Þ

r ¼

Nþ1 1 X hk 2aH k¼1

Z

zk

zk1

Z

a a

rkyy dxdz ¼

Nþ1 1 X hk 2aH k¼1

Z

a

a

r kyy dx

ð34Þ

r kyy is calculated by averaging the transverse stress component through the thickness. Using the edge boundary conditions at x = ±a, the averaged stress components are as below. 1 2a 1 2a

Z Z

r kzz ðxÞdx ¼ 0

a

a



cc ¼ s21  s23

r kxx ðxÞdx ¼

 1 ~ ~ k ðaÞ þ rk /k1 ðaÞ  / hk a

s31 s33

ð38Þ

ec and substituting in Eqs. (30) and (31),

rþ

 s13  s r þ s12  s13 32 s s33 s33

ð39Þ

s31 1 s s r þ r  32 s  34 DT s33 s33 s33 s33 

rþ

ð40Þ

 s23  s r þ s22  s23 32 s s33 s33

 s34 þ s24  s23 DT s33

ð41Þ

and by inversing the stress–strain relations in (39)–(41) the constitutive law for cracked laminate is obtained. In compact notation it yields to,

8 > <

8

9

9

r> > < ec > = = r ¼ ½Q Deg ec þ fRgDeg DT > > : > : > cc ; s; where ½Q Deg and fRgDeg are the degraded stiffness matrix and thermal expansion coefficient vectors for cracked unit cell. The last stage is the calculation of transverse, shear and thermal damage parameters for the cracked layer as a function of degraded mechanical properties of the laminate. The degraded stiffness matrix and thermal expansion coefficient vectors of the cracked laminate can be decomposed to the related components of jth cracked layer and all remaining layers as,

H ½Q ¼ hj j

fRg ¼

(

H hj

½Q Deg

j1 Nþ1 X k k 1 X  hk ½Q þ hk ½Q  H k¼1 k¼jþ1

( fRgDeg 

!)

j1 Nþ1 X 1 X k k hk fRg þ hk fRg H k¼1 k¼jþ1

ð43Þ !) ð44Þ

j

½Q  and fRg are the stiffness matrix and thermal expansion coefficients vectors for the jth cracked layer. The right hand side of the Eqs. (43) and (44) is known from the reduced stiffness of laminate in Eq. (42) and the thermo-mechanical properties of all layers except the cracked layer. These layers may contain no damage, so the virgin properties are used to calculate their stiffness matrix and thermal expansion coefficients vector. If they contain a converged damage state, so the stored damage parameters are used to calculate the properties of the layers. Therefore, the stiffness matrix and thermal expansion coefficients vector of the cracked layer is available and the degraded thermo-elastic constants can be calculated. The damage parameters of the cracked layer are defined as, j

a

a

ec ¼ 

j

each set of sij (compliance matrix components) is calculated. The axial and shear strain components of the damaged unit cell, ec and cc are calculated from equations (22) and (23). The effective in-plane transverse stress, r⁄ is calculated by summation of corresponding stress components at each layer as,

s31 s33



 s34 þ s14  s13 DT s33

j

By applying the following set of load cases to the presented model,

ð37Þ

   Nþ1  X hk 1 1 ~  k ~ k ðaÞ þ rk / / e  S ðaÞ  k1 12 hk a H Sk22 c k¼1     1 ~ ~ k ðaÞ þ sk  a ^ 2 DT wk1 ðaÞ  w Sk26 hk a

c

ec ¼ s11 r þ s12 s þ s13 ec þ s14 DT

 1 ~ ~ k ðaÞ þ sk wk1 ðaÞ  w hk a

r ¼

ð29Þ

where e⁄ is the transverse normal strain of the undamaged laminate. ½Q  and fRg are the stiffness matrix and thermal expansion coefficient vector of the laminate and the non-zero components of them are presented in Appendix C. The tensile and shear stress components in each layer of the un-cracked laminate, rk and sk in Eqs. (B2) and (B3) are obtained from (29). To find the stiffness reduction of the laminate due to the presence of matrix cracks, the stress and displacement fields of the unit cell are used. Because of the linear elasticity assumption, it follows that the effective strains and stress, ec, cc and r⁄ are linearly related to the r, s, ec and DT parameters as below.

r kxy ðxÞdx ¼

By inserting (35)–(37) to the average of (11), the expression for r⁄ is obtained as:

ð28Þ

For the undamaged laminate, the strain components of all layers are equal and the effective applied axial, transverse and shear stress components can be written as:

a

ð35Þ

d2 ¼ 1 

ð36Þ

d6 ¼ 1 

j

Ej2 E02 Gj12 G012

ð45Þ

ð46Þ

942

j

da ¼ 1 

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

aj2 a02

ð47Þ

where superscript 0 refers to the virgin material constants. Thus for each cracked layer under investigation, with proper transformation, the loss of stiffness of the layer as a function of crack density can be calculated from micromechanical (unit cell) model. It is worth to note that the procedure explained in this section is independent from finite element data like stress or strain tensors of the element. The damage criterion, which will be explained in the next section, bridges between the FEM (macromechanics) and the micromechanic damage model. 4. Energy based criterion and damage surface A composite laminate under various loading conditions experiences initiation and propagation of multiple microcracks in each ply at early stages of loading. By increasing the load and coalescence of neighboring microcracks, fully developed through the thickness cracks appear within a ply. These cracks will be increased in number and trigger other types of damage until the final failure of a laminate. The conventional fracture mechanics does not answer the problem of predicting the extent or number of matrix cracks. Hashin [17] extended the fracture mechanics method and proposed the finite fracture mechanics. In the finite fracture mechanics, the damage is assumed to occur when the finite energy released by the generation of new crack faces exceeds the toughness of material. Using the finite fracture principle, several models have been developed to predict the initiation and growth of the matrix cracking damage. Nairn [18,19] obtained the stress field of the unit cell using the variational approach. In this model, the strain energy of the unit cell is calculated by doubling the number of matrix cracks at each step and the strain energy release rate is evaluated according to the new generated crack area. Lim and Li [20] also used the

Model, Discretize and Apply Load and BCs.

variational approach to find the stress field of the [S/90n]s laminates under uni-axial tensile loading condition and to obtain the energy release rate by the incremental and continuous variation of crack density. For cross-ply laminates containing transverse and longitudinal cracks, Rebie‘re et al. [21,22] proposed an energetic criterion using the decomposition of strain energy to predict the initiation and growth of the damage. The strain energy release rate was calculated as a continuous function of unit cell length and width for different damage mechanisms. These studies need the stress field of the unit cell which is sometimes a cumbersome task for laminates with general layup configuration under arbitrary in-plane loading condition as stated in the present study. The strain energy of the unit cell can be written as,

U cell ¼

1 2

Z

heifrgdV

ð48Þ

V

According to this equation, the strain energy of the unit cell can be calculated in two different ways. In the first method, the energy relation is expanded in terms of the stress tensor and compliance matrix. The stress field of the unit cell under certain loading condition is obtained from the micromechanic damage models such as variational approach [18,20]. The compliance matrix consists of the virgin material properties and perturbed the stress filed of the unit cell as a function of the crack density and controls the damage criterion. The method is restricted to the special laminates and loading condition due to the difficulty in calculating the stress field of more general problems. To overcome such limitations, a method is developed and used in the present study, where the energy equation is expanded in terms of the stiffness matrix and the strain vector as,

U cell ¼

1 2

Z

heifrgdV ¼

V

b 2

Z A

hei½Q Deg fegdA ¼ U cell þ U cell I II

Coordinate Transformation to Layer j (j=1..N+1)

Define Material Properties and Layup Configuration

State variable λ

Load Increment

Construct and Solve Governing ODE Systems

Update Strain, State Variables and Damage Parameters

Calculate Laminate Stiffness Degradation

Construct Element Stiffness Matrix and Solve [K]{u} = {F}

λ = λ ± Δλ

Calculate Strain Energy, Decompose it and Find its Release Rate Damage Surface

Residual Force Calculation

N

Layer Convergence N

Y

Convergence

Find Properties Deterioration of Layer and Calculate Damage Parameters

Y N

Final Load Y

Y

Output Results Finite Element Procedure

Laminate Convergence

N

Micromechanic Damage Model

Fig. 2. Progressive damage analysis algorithm.

ð49Þ

943

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

where ½Q Deg is the degraded stiffness matrix of the unit cell and is calculated from the micromechanic damage model. hei is the transpose of strain vector that is used from FEM results at the integration point. V is the volume of the unit cell, (V = 2abH) and H is the thickness of laminate. So the crack density and degraded stiffness matrix controls the strain energy of the unit cell. For matrix cracking damage mechanism, it is assumed that the in-plane normal and shear strain components cause the opening and sliding mode of the crack respectively. The tensile in-plane strain components, exx and eyy drives the opening mode and shear strain components, and exy drives the sliding mode of the matrix cracking. The strain energy for opening (I) and sliding (II) modes are,

U cell ¼ I

Fig. 4a. Normalized modulus versus applied strain for [15/15/904]s laminate.

1 1 hexx iþ ðQ 11 exx þ Q 12 eyy þ Q 16 exy Þ þ eyy þ ðQ 12 exx 2 2 þ Q 22 eyy þ Q 26 exy Þ

ð50Þ

Table 1 Thermo-mechanical properties of the laminates of study.

E1 (MPa) E2 (MPa) G12 (MPa) G23 (MPa)

m12 a1 (106/°C) a2 (106/°C) Ply thickness (mm) GIC (K J/m2) GIIC (K J/m2) DT (°C) ⁄

Glass/epoxy [32]

IM600/#133 [9]

44,730 12,760 5800 4490⁄ 0.297 8.6 22.1 0.144 0.175 [26] 1.561 [26] 105

147,000 8310 4700 2866⁄ 0.352 0.55 22.5 0.14 0.04⁄ 1.5⁄ 105

Fig. 4b. Normalized shear modulus and poison’s ratio versus applied strain for [15/ 15/904]s laminate.

Assumed value.

1 1

EX/EX

0

0.8

0.8

0.6

0.4 0.6 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

0

0.5

2

ε 0 (%)

1.5

2

Fig. 5a. Normalized modulus versus applied strain for [30/30/904]s laminate.

GXY/GXY0

Fig. 3a. Normalized modulus versus applied strain for [02/904]s laminate.

1

ε0 (%)

1 0.9 0.8

ν

ν

ν XY/ν XY0

0.7

0

0.5

1

1.5

2

ε Fig. 3b. Normalized shear modulus and degraded poison’s ratio versus applied strain for [02/904]s laminate.

0.6 0.5

0

0.5

1

1.5

2

ε0 (%) Fig. 5b. Normalized shear modulus and poison’s ratio versus applied strain for [30/ 30/904]s laminate.

944

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

EX/EX

0

and tcr is the thickness of the cracked layer. Using the derivative of Eqs. (52) and (53) with respect to k, the strain energy release rate associated with the transverse cracking damage is,

GI;II

0

0.5

1

1.5

2

ε0 (%) Fig. 6a. Normalized modulus versus applied strain for [40/40/904]s laminate.

@U cell @U cell 1 I;II I;II @k ¼ ¼ ¼ 2tcr @Acr @k @Acr

U cell I;II

@U cell I;II þk @k

! ð54Þ

It’s clear that the L and b (length and width of the considered plate) have been in the energy release rate equation and G is the function of the crack density only. U cell and U cell I II are calculated from Eqs. (50) and (51) using the reduced stiffness parameter from the damage model and the strain tensor from the FEM results. The critical energy release rate for microcracking is the criterion to predict the damage evolution. A number of mixed mode fracture criterion have been proposed for describing the interlaminar fracture of composites [23–25]. These damage surfaces are formulated on the basis of energy release rate of the normal and shear damage modes. In the present study, the fracture criteria proposed in [26] is used for the damage surface as follows,

2 GI GII þ 1P0 GIC GIIC

ð55Þ

where GIc and GIIc are the critical values of the energy release rate for normal and shear damage modes which are used from experimentally determined data in the literature. 5. Numerical procedure The developed damage model is implemented in the FEM formulation of a layerwise element. The full layerwise plate theory with three degrees of freedom for each node is used to approximate the displacement field. In contrast with the 3-D solid elements, the layerwise elements provide the capability of discretizing the

λ (Cr/mm)

Damage Parameters

Fig. 6b. Normalized shear modulus and poison’s ratio versus applied strain for [40/ 40/904]s laminate.

ε Fig. 6c. Damage parameters for 90° layers of [40/40/904]s laminate versus applied strain at an integration point of a single element.

1 exy ðQ 16 exx þ Q 26 eyy þ Q 66 exy Þ 2

Fig. 7. Crack density versus applied stress for [±h/904]s glass/epoxy laminates.

ð51Þ

For compressive normal stresses, the cracks are closed, hence no evolution occurs. Thus the notation in (50) is defined as hai+ = a if a P 0 otherwise hai+ = 0. To calculate the strain energy release rate in mixed mode conditions a plate with the length of L and width of b is considered at the integration point. In the presence of M matrix cracks with uniform distance of 2a, the plate is divided to M unit cells. The strain energy of the entire plate is obtained from the summation of all energy terms of the unit cells, cell U I;II ¼ MU cell I;II ¼ LkU I;II

ð52Þ

where M = L/2a = Lk and k is the crack density. The cracks surface area is,

Acr ¼ 2Mbtcr ¼ 2Lbktcr

crack density (cr/mm)

U cell II ¼

Stress (MPa)

ð53Þ

ε0 (%) Fig. 8a. Crack density versus applied strain in h layer of [0/554/554/00.5]s laminate.

EX/EX0

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

λ (cr/mm) Fig. 8b. Normalized elastic modulus versus crack density for [0/554/554/00.5]s laminate.

structure by a 2-D mesh. Therefore, the amount of input data is reduced and the number of operations for constructing the element stiffness matrix is decreased [27]. The element is a 2-D eight node serendipity quadrilateral element in conjunction with one linear 1-D Lagrange element through the thickness [16]. It is developed as a user-defined element (UEL) in the Abaqus commercial finite element code [28]. The standard code contains an interface that allows the users to define the element and it is capable to link the Fortran language based subroutines. The governing differential equation in (B25) is constructed and solved numerically using the methodology presented in [29] for any number of layers during the analysis. This is a key point for progressive damage analysis in the present study as the damage parameters are updated in an iterative manner during the analysis. Fig. 2 shows the analysis algorithm consists of two main parts distinguished in the dashed boxes. The left box, finite element procedure, is the general procedure to solve the nonlinear problems including the progressive damage analysis. It uses the capability of standard code for modeling, discretization, applying the initial conditions, assembly of element stiffness matrices, calculating of

945

the displacement filed for the nonlinear problem iteratively and finally post-processing of the results. The second part of the algorithm is the developed damage model used to calculate the deteriorated thermo-mechanical properties of the layers of the laminate. At each integration point of each layer, with proper coordinate transformation, the typically presented unit cell in Fig. 1 is extracted from the laminate and the stiffness reduction is calculated. The initial value of crack density, k, is 0.01 (cracks/mm) for each layer and the same value is used as an incremental crack density, Dk to calculate the energy release rate in Eq. (54) numerically. Using the damage surface, the damage initiation and progression is investigated for the layer. While the crack density is converged for the layer, the damage parameters are calculated and stored to be used during the progressive damage analysis of remaining layers. The procedure is continued to the last layer of the laminate. Since the solution in a layer depends on the stiffness of the remaining layers, a converged iteration for a lamina does not guarantee the full convergence of the problem. The damages in the remaining layers are updated, unless the remaining layers remain undamaged. Therefore, the laminate convergence loop is repeated until the convergence of the k value for every layer. The updated damage parameters, represent the new properties of the damaged layer, are used to construct the stiffness matrix of the laminate. 6. Results To verify the developed procedure, the predicted stiffness reduction and crack density evolution for different layups are compared with the available experimental/numerical results. In the previously performed investigations, the stiffness reduction as a function of crack density or applied strain and the damage evolution including the onset and growth of k in a layer of the laminate have been performed. These results are obtained at an integration point of an element with unit length under uni-axial tensile displacement and boundary condition as depicted in Fig. 3a. The stress field of such model is uniform far from the edge and the

Fig. 9. Normalized elastic modulus versus crack density evolution for [h/h]s laminates.

946

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

predicted results are comparable with the available results for unnotched laminates. Most of the studies in this scope are limited to simple layup configurations and prediction of k in single lamina and there are only a few resources to find the results for more complex problems. To evaluate the accuracy of the proposed model, Varna et al. [30–32] experimental results for glass/epoxy laminate, Barbero’s model results [11,12] for the same material and Yakozeki’s [9] numerical model results for IM600/#133 laminates are used. The virgin material properties are presented in Table 1. In Figs. 3–6 variations of the deteriorated mechanical properties including the normal modulus, shear modulus and poison’s ratio of [h/h/904]s glass/epoxy laminates for h = 0°, 15°, 30° and 40° are presented as a function of the applied strain. For h = 0° the governing differential equation of the unit cell was decoupled and the solution algorithm is somehow different. Superscript 0 refers to the undamaged laminate properties. The degraded normal modulus of the laminate are obtained for different ply refinements to exhibit the effect of such technique on the accuracy of the results. Increasing the number of sub-layers in the micromechanic damage model improves the accuracy of the results. The degraded shear modulus and poison’s ratio of the laminate are obtained for ply refinement of np = 5 only. No shear modulus data has been reported in the literature and for normal modulus and poison’s ratio; the results of the present study in these figures follow the available results in the literature with acceptable agreement. Variations of the damage parameters d2, d6 and da versus the applied strain are presented in Fig. 6c. As explained in Section 1, these parameters control the deterioration of the transverse modulus, shear modulus and thermal expansion coefficient in the 90° layers. The figure shows that the damage parameters are zero-valued until the applied strain reaches about 0.6%. By increasing the applied strain, the matrix cracking is detected and the damage parameters grow. Theses parameters control the damage state at the integration point. The evolution of crack density, k, for [h/h/904]s laminates are compared with the available experimental results in Fig. 7. As the stiffness of the constraint layer is increased, the onset of damage occurs in the higher load levels. The initiation and progress of the matrix cracking in the early stage of loading are in good agreement with the experimental results. But, by increasing the load, the present model predicts larger crack density values with respect to the experimental values especially for [40/40/904]s laminate. For this laminate, no more matrix cracks detected in the test and the saturated crack density happens at k = 0.4. One reason is the effect of other damage mechanisms like the induced delamination or fiber fracture that are not included in the modeling. Fig. 8a shows the variation of crack density evolution versus applied strain in h layers of [0/554/554/00.5]s glass/epoxy laminate in which the cracks are subjected to both mode I and mode II fractures. The prediction of initiation and evolution of matrix cracks are in good agreement with the experimental findings. Barbero et al. [12] predicted the onset of damage in smaller applied strain and shown the smooth behavior of the damage evolution, while the experiential results show a rapid increase of crack density during the loading implementation. The present study predicts the onset of damage at the strain value of about e = 1.3% for np = 5 and shows the rapid evolution of matrix cracking which is more compatible with the experimental results. The degraded elastic modulus of the laminate is depicted in Fig. 8b and shows a good agreement with the experimental results. This may be due to the considering of simultaneous damage evolution in both layers. The IM600/#133 [h/h]s laminates with h = 15°, 30°, 45°, 60° and 75° are considered for progressive damage analyses. The predicted results by Yakozeki et al. [9] using 2D shear-lag method in oblique coordinate system, are used to evaluate the results with

the developed procedure in this study. The unit cell in this model is confined between matrix cracks in both bunches of layers. They obtained the degradation of mechanical properties based on the evolution of cracks in h layers while the crack density in h layer is k = 1 (crack/mm) without occurrence of simultaneous crack evolution. Due to the limitations in their analytical model, it’s difficult to consider the evolution of damage in both layers simultaneously. However, using the developed model in the present study, there are no limitations in the simultaneous evolution of matrix cracks in any number of layers. Fig. 9 shows the variation of normalized elastic modulus versus crack density in the layers of [±h]s laminates with different h angles. The damage onset occurs in h layers firstly and the elastic modulus starts to degrade from virgin property. By increasing the loading value, the value of k in h layers and the degradation of EX are increased. This is also clear in this figure for k(h) curves in laminates

Fig. 10. Normalized shear modulus versus crack density evolution for [h/h]s laminates.

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

with h 6 6°. By increasing the crack density values the predicted elastic modulus becomes smaller than those by shear-lag approach (Yakozeki et al.) which is due to the consideration of possible simultaneous damage progress in both bunches of layers. Variations of normalized shear modulus versus crack density are also shown in Fig. 10. The same conclusion as explained above for variation of elastic modulus may be interpreted here. When the matrix cracks initiates at h layers, the GXY property is deteriorated rapidly because of the progression of the damage in both bunches of the layers. Fig. 11 shows the evolution of matrix cracking in both layers of [±45]s glass/epoxy laminate with a central hole. Two snapshots in different loading stages are presented to show the onset and evolution of damage in the layers of laminate. Due to the presence of hole, the stress field of the laminate is non-uniform and the previously developed models in the literature could not properly predict the damage onset and growth of such component. At the early stage of loading, the damage grows in stress concentration region near the hole. While the stiffness of the laminate in this region is deteriorated due to the presence of the damage, the load shedding is occurred and the neighboring area carries the load. The damage progresses in both layers simultaneously and the interaction of

947

stiffness reduction of the layers controls the damage state of the laminate at the higher load levels. These predicted damage patterns in various layers of the laminate with central geometry discontinuity demonstrate the capability of the developed numerical procedure in the present study to find the damage evolution in the laminates with complex geometry and boundary conditions. 7. Conclusion Micromechanic damage model based on the stress transfer model was successfully implemented in a finite element based tool to predict the matrix cracking onset and evolution in general symmetric laminates. Energy-based criterion and damage surface was used as the evolution law and the unit cell length dependency removed from the formulations to be able to implement the model at integration points. A full layer-wise element was developed in the commercial finite element software and the damage model was implemented inside the formulations of the stiffness matrix of the element. The obtained results for various layup configurations were validated favorably using the available numerical/experimental results from the literature. Also, progressive damage analysis

Fig. 11. Matrix crack evolution in layers of IM600/#133 [±45]s laminate.

948

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

results for laminates with various layups which have been rarely obtained by the researchers were presented here.

rkxx ¼

/k1 ðxÞ  /k ðxÞ þ rk hk

ðB2Þ

Appendix A

rkxy ¼

wk1 ðxÞ  wk ðxÞ þ sk hk

ðB3Þ

 g are, Non-zero components of the matrices [F], {G} and fa

F k11

¼

Sk11



F k12 ¼ Sk13 

F k15 ¼ Sk16 

F k22 ¼ Sk33 

F k25

¼

Gk2 ¼

Gk5 ¼

Sk12 Sk23 Sk22 Sk12 Sk26  2 Sk23

F k44 ¼ Sk55

 2 Sk26 Sk22

ðA7Þ

ðA9Þ

Sk22 Sk26

ðB5Þ

Integration with respect to x caused to appear through the thickness displacement Wk(x) which can be found from boundary conditions. The expression for uk is calculated from definition of exz and Eq. (B5) as, @uk @wk 1 2exz ¼ þ ¼ F k33 rkxz þ F k34 rkyz ! uk ¼ 2hk @z @x 8 9   > > < 2hk ðz  zk Þ F k33 /0k ðxÞ þ F k34 w0k ðxÞ þ = h  i    2 k k 0 0 0 0 > : ðz  zk Þ F 33 /k ðxÞ  /k1 ðxÞ þ F 34 wk ðxÞ  wk1 ðxÞ > ;    ðz  zk Þ3  000 000 ðz  zk Þ /000 k1 ðxÞ  /k ðxÞ  4hk /k ðxÞ 24hk    i ðz  zk Þ2 h  hk F k22 Y0k ðxÞ þ F k12 /0k1 ðxÞ  /0k ðxÞ þ F k25 w0k1 ðxÞ  w0k ðxÞ 2hk  ðz  zk ÞW 0k ðxÞ þ U k ðxÞ; zk1 6 z 6 zk ðB6Þ

Sk12 Sk22 Sk23 Sk22 Sk26 Sk22

Uk(x) is the axial displacement on the kth interface at z = zk arising from integration with respect to z. Similarly the expression for vk,

a^ 2

ðA11Þ

a^ 2

ðA12Þ

a^ 2

ðA13Þ

Stress and displacement fields of the unit cell are calculated by constructing and solving of the fourth order governing differential equation of the problem. Using equilibrium equations (Eq. (10)) and definition of function / (Eq. (11)),

r

   ðz  zk Þ2  ðz  zk Þ /00k1 ðxÞ  /00k ðxÞ  3hk /00k ðxÞ 6hk þW k ðxÞ; zk1 6 z 6 zk þF k22

ðA10Þ

Sk22

Appendix B

k zz

@wk  k2 DT ¼ F k12 rkxx þ F k22 rkzz þ F k25 rkxy þ Gk2 ec þ a ðB4Þ @x h   k ðxÞ W k ðxÞ ¼ ðz  zk Þ F k12 /k1 ðxÞ/ þ rk hk

  w ðxÞ  wk ðxÞ  k2 DT þ F k22 Yk ðxÞ þF k25 k1 þ sk þ Gk2 ec þ a hk

 F k22

Sk23

a k5 ¼ a^ 6 

ðA6Þ

ðA8Þ

Sk22

a k2 ¼ a^ 3 

ekzz ¼

ðA5Þ

Sk12

a k1 ¼ a^ 1 

ðA2Þ

ðA4Þ

Sk22

F k34 ¼ Sk45 ;

ðA1Þ

The crack will perturb the stress field and in the absence of transverse cracks, / and w are vanished and only the remote applied stresses on undamaged laminate, rk and sk remain. Strain and displacement in the thickness direction are:

ðA3Þ

Sk22

Sk22

¼ Sk66 

Gk1 ¼

Sk22

Sk23 Sk26

F k33 ¼ Sk44 ;

F k55

 2 Sk12

 ðz  zk Þ2  00 ¼ /k1 ðxÞ  /00k ðxÞ  ðz  zk Þ/00k ðxÞ þ Yk ðxÞ; zk1 6 z 6 zk 2hk

@ v k @wk 1 2eyz ¼ þ ¼ F k34 rkxz þ F k44 rkyz ! v k ¼ 2hk @z @y 8 9   > > < 2hk ðz  zk Þ F k34 /0k ðxÞ þ F k44 w0k ðxÞ þ = h  i    > : ðz  zk Þ2 F k34 /0k ðxÞ  /0k1 ðxÞ þ F k44 w0k ðxÞ  w0k1 ðxÞ > ; þ V k ðxÞ þ ec y; zk1 6 z 6 zk

The latest term in the above equation is added to satisfy the generalized plane strain condition of Eq. (6). The boundary conditions of the problem are categorized to two sets. The interface conditions which implying the continuity of stress and displacement at the layer interfaces and edge conditions. The first set is used to find recurrence relation for unknown terms which are arisen from integration in the previous expressions. On any interface z = zk, mid-plane z = z0 and free surface of the laminate z = zN+1,

(

ðB1Þ

Yk(x) is appeared due to the integration with respect to z and will be find later using the interface boundary conditions. Using the equilibrium equations and definition of / and w functions,

ðB7Þ

on z ¼ zk !

k kþ1 k kþ1 rkzz ¼ rkþ1 zz ; rxz ¼ rxz ; ryz ¼ ryz

uk ¼ ukþ1 ; v k ¼ v kþ1 ; wk ¼ wkþ1

;

k ¼ 1; 2; . . . ; N ðB8Þ

on z ¼ z0 ! r1xz ¼ r1yz ¼ u1 ¼ 0

ðB9Þ

949

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

on z ¼ zNþ1 ! rNþ1 ¼ rNþ1 ¼ rNþ1 ¼0 zz xz yz

ðB10Þ

The edge boundary conditions are used to solve the governing differential equations of the problem and will be discussed later. Using B8,B9 and B10 the recurrence relation and Wk can be calculated as,

 hk  00 Yk1 ðxÞ ¼ / ðxÞ þ /00k ðxÞ þ Yk ðxÞ; 2 k1 k ¼ N þ 1; N; . . . ; 1; YNþ1 ðxÞ ¼ 0

ðB11Þ

h   k ðxÞ W k ðxÞ ¼ hk F k12 /k1 ðxÞ/ þ rk hk

  w ðxÞ  wk ðxÞ  k2 DT þ F k22 Yk ðxÞ þF k25 k1 þ sk þ Gk2 ec þ a hk 2

þF k22

 hk  00 /k1 ðxÞ þ 2/00k ðxÞ þ W k1 ðxÞ 6

ðB12Þ

Also by definition of DUk(x)  Uk(x)  UN+1(x) and DVk(x)  Vk(x)  VN+1(x) using B8,B9 and B10,   i hk h k  0 F / ðxÞ þ /0k1 ðxÞ þ F k34 w0k ðxÞ þ w0k1 ðxÞ 2 33 k "  # 3  hk hk F k22 Y0k ðxÞ þ F k12 /0k1 ðxÞ  /0k ðxÞ h  000  F k22 k /000   k1 ðxÞ þ 3/k ðxÞ  24 2 þF k w0 ðxÞ  w0 ðxÞ

DU k1 ðxÞ ¼ 

25

k1

k

  i hk h k  0 F 34 2/k ðxÞ þ /0k1 ðxÞ þ F k44 2w0k ðxÞ þ w0k1 ðxÞ 6 þ V k ðxÞ þ ec y

v k ¼ 

Inserting the (B18)–(B21) and (B22) to (B16) and (B17) and by some algebraic manipulation, 2  i h h  ð4Þ ð4Þ U 0k ðxÞ ¼ F k22 k hk /k1 ðxÞ þ 4/k ðxÞ þ 20F k22 Y00k ðxÞ 120 2 3    k k k k 00 00 þ F ðxÞ þ 2F  F ðxÞþ 2F / w k k 33 12 34 25 hk 6 7    þ 4 5 6 F k33 þ 2F k12 /00k1 ðxÞ þ F k34 þ F k25 w00k1 ðxÞ i 1h k  F 11 ð/k1 ðxÞ  /k ðxÞÞ þ F k15 ð2wk1 ðxÞ  wk ðxÞÞ þ hk F k12 Yk ðxÞ hk   h  k1 DT  k W 00k ðxÞ þ F k11 rk þ F k15 sk þ Gk1 ec þ a 2 ðB23Þ

   i hk h k V 0k ðxÞ ¼ F 34 þ F k25 2/00k ðxÞ þ /00k1 ðxÞ þ F k44 2w00k ðxÞ þ w00k1 ðxÞ 6 i 1h k þ F 15 ð/k1 ðxÞ  /k ðxÞÞ þ F k55 ð2wk1 ðxÞ  wk ðxÞÞ þ hk F k25 Yk ðxÞ h k   k5 DT þ F k15 rk þ F k55 sk þ Gk5 ec þ a

þ hk W 0k ðxÞ þ DU k ðxÞ; DU Nþ1 ðxÞ ¼ 0

ðB24Þ ðB13Þ

  i hk h k  0 DV k1 ðxÞ ¼  F 34 /k ðxÞ þ /0k1 ðxÞ þ F k44 w0k ðxÞ þ w0k1 ðxÞ 2 þ DV k ðxÞ; DV Nþ1 ðxÞ ¼ 0 ðB14Þ Now the stress–displacement relations in (B1)–(B4) and (B5) are functions of unknowns, / and w and their derivatives. As some of the boundary conditions are not satisfied exactly, an attempt is made to satisfy them in an average manner. The average of any quantity through the thickness is denoted by,

1 gk ðx; yÞ ¼ hk

Z

zk1

hk ¼ zk  zk1

 k  k k k k exx ¼ @ uk ¼ F k r k k 11 xx þ F 12 rzz þ F 15 rxy þ G1 ec þ a1 DT @x

ðB16Þ

 k @ v k @u k  k k  k þ Fk r k k 2exy ¼ þ ¼ F k15 r xx 25 zz þ F 55 rxy þ G5 ec þ a5 DT @y @x

ðB17Þ

and corresponding expression for stress components rzz, rxx and rxy in B1,B2 and B3,

/k1 ðxÞ  /k ðxÞ þ rk hk

w ðxÞ  wk ðxÞ r kxy ¼ rkxy ¼ k1 þ sk hk

r kzz ¼

k¼1

 hk  00 / ðxÞ þ 2/00k ðxÞ þ Yk ðxÞ 6 k1

The average of displacement terms uk and are:

ðB25Þ

ðB15Þ

hk is the thickness of kth layer. By averaging exx and exy,

r kxx ¼ rkxx ¼

using (B23) and (B24) in conjunction with the derivative of (B13) and (B14) with respect to x and insert the expressions for the definition of Uk(x)  UN+1(x)  DUk(x) = 0 and Vk(x)  VN+1 (x)  DVk(x) = 0, the following fourth order coupled ordinary differential equations are obtained: 8 N >X ð4Þ > > Akj /k ðxÞ þ Bkj /00k ðxÞ þ C kj /k ðxÞ þ bkj w00k ðxÞ þ ckj wk ðxÞ ¼ 0 > < k¼1 ; j ¼ 1;2;...;N N > X > 00 00 > > E / ðxÞ þ F / ðxÞ þ e w ðxÞ þ f w ðxÞ ¼ 0 : kj k kj k kj k kj k

zk

g k ðx; y; zÞdz;

ðB22Þ

The coefficients of the system are the matrices with order of N and are functions of material constants, ply thicknesses and damage parameters. As stated earlier, some of the boundary conditions satisfied in the average manner and to increase the accuracy of the method, ply refinement technique is used. Using such technique increases the rank of the coefficient of matrices in (B23) and also computation time. Appendix C The components of matrix [Q] and vector fRg are given as,

Q 11

0 1 Nþ1 1 @X F k55 A ¼ hk H k¼1 F k F k  ðF k Þ2 11 55 15

ðC1Þ

Q 12

0 1 Nþ1 1 @X F k15 Gk5  F k55 Gk1 A ¼ hk H k¼1 F k F k  ðF k Þ2

ðC2Þ

0 1 Nþ1 1 @X F k15 A ¼ hk H k¼1 F k F k  ðF k Þ2 11 55 15

ðC3Þ

ðB18Þ

ðB19Þ

ðB20Þ

vk in (B6) and (B7)

3        k ¼  hk F k /0 ðxÞ þ F k w0 ðxÞ  hk F k /0 ðxÞ þ F k w0 ðxÞ  F k hk /000 ðxÞ þ 4/000 ðxÞ u 34 k 34 k1 22 k 3 33 k 6 33 k1 120 k1    i hk hk h hk F k22 Y0k ðxÞ þ F k12 /0k1 ðxÞ  /0k ðxÞ þ F k25 w0k1 ðxÞ  w0k ðxÞ  W 0k ðxÞ þ U k ðxÞ  6 2 ðB21Þ

11 55

Q 13

15

0 " #1 N þ1 1 X 1 1 k k k k k k k k k k A Q 22 ¼ @ hk  G ðG S  G S Þ  G ðG S  G S Þ 1 1 66 5 16 5 1 16 5 11 H k¼1 F k F k  ðF k Þ2 Sk22 11 55 15 ðC4Þ

950

G. Sadeghi et al. / Composite Structures 108 (2014) 937–950

0

Q 23 ¼

1

Nþ1 1 @X F k Gk  F k11 Gk5 A hk 15 1 H k¼1 F k F k  ðF k Þ2 11 55

Q 33

ðC5Þ

15

0 1 Nþ1 1 @X F k11 A ¼ hk H k¼1 F k F k  ðF k Þ2 11 55 15

0 1 Nþ1  k1  F k15 a  k5 1 @X F k55 a A R1 ¼  hk H k¼1 F k F k  ðF k Þ2 11 55

ðC6Þ

ðC7Þ

15

0 " #1 N þ1 ^2 1 @X 1 a k k k k k k k k k k A   F 15 a  5 Þ  G5 ðF 15 a  1  F 11 a 5 Þ  R2 ¼ hk G ðF a H k¼1 F k F k  ðF k Þ2 1 55 1 Sk22 11 55 15 ðC8Þ

0 1 N þ1  k1  F k11 a  k5 1 @X F k15 a A R3 ¼ hk H k¼1 F k F k  ðF k Þ2 11 55

ðC9Þ

15

References [1] Liu S, Nairn JA. The formation and propagation of matrix microcracks in crossplylaminates during static loading. J Reinf Plast Compos 1992;11:158–73. [2] Tan SC, Nuismer RJ. A theory for progressive matrix cracking in composite laminates. J Compos Mater 1989;23:1029–47. [3] Fan J, Zhang J. In-situ damage evolution an micro/macrotransition for laminated composites. Comput Sci Technol 1993;47:107–18. [4] Zhang J, Fan J, Soutis C. Analysis of multiple matrix cracking in [±hm/90n]s composite laminates, part 2: development of transverse ply cracks. Composites 1992;23:299–304. [5] McCartney LN. Energy-based prediction of failure in general symmetric laminates. Eng Fract Mech 2005;72:909–30. [6] Hosseini Toudeshky H, Farrokhabadi A, Mohammadi B. Consideration of concurrent transverse cracking and induced delamination propagation using a generalized micro–meso approach and experimental validation. Fatigue Fract Eng Mater Struct 2012;35:885–901. [7] Farrokhabadi A, Hosseini Toudeshky H, Mohhammadi B. A generalized micromechanical approach for the analysis of transverse crack and induced delamination in composite laminates. Compos Struct 2011;93:443–55. [8] Farrokhabdi A, Hosseini Toudeshky H, Mohammadi B. Damage analysis of laminated composites using a new coupled micro–meso approach. Fatigue Fract Eng Mater Struct 2010;33:420–35. [9] Yokozeki T, Aoki T. Overall thermoelastic properties of symmetric laminates containing obliquely crossed matrix cracks. Compos Sci Technol 2005;65:1647–54. [10] Cortes DH, Barbero E. Stiffness reduction and fracture evolution of obliquely matrix cracks in composite laminates. Ann Solid Struct Mech 2010;1:29–40.

[11] Barbero EJ, Cortes DH. A mechanistic model for transverse damage initiation, evolution, and stiffness reduction in laminated composites. Composites Part B 2010;41:124–32. [12] Barbero E, Sgambitterra G, Adumitroaie A, Martinez X. A discrete constitutive model for transverse and shear damage of symmetric laminates with arbitrary stacking sequence. Compos Struct 2011;93:1021–30. [13] Sgambitterra G, Adumitroaie A, Barbero E, Tessler A. A robust three-node shell element for laminated composites with matrix damage. Composites: Part B 2011;42:41–50. [14] Sadeghi G, Hosseini Toudeshky H, Mohammadi B. In-plane progressive matrix cracking analysis of symmetric cross-ply laminates with holes. Fatigue Fract Eng Mater Struct 2013. http://dx.doi.org/10.1111/ffe.12113. [15] McCartney LN. Model to predict effects of triaxial loading on ply cracking in general symmetric laminates. Comput Sci Technol 2000;60:2255–79. [16] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. 2nd ed. CRC Press; 2004. [17] Hashin Z. Finite thermoelastic fracture criterion with application to laminate cracking analysis. J Mech Phys Solids 1996;44:1129–45. [18] Nairn JA. The strain energy release rate of composite microcracking: a variational approach. J Compos Mater 1989;23(11):1106–29. [19] Nairn JA, Hu S. The initiation and growth of delaminations induced by matrix micro-cracks in laminated composites. Int J Fract 1992;57:1–24. [20] Lim S, Li S. Energy release rates for transverse cracking and delaminations induced by transverse cracks in laminated composites. Composites: Part A 2005;36:1467–76. [21] Rebie‘re JL, Gamby D. Initiation and growth of transverse and longitudinal cracks in composite cross-ply laminates. Compos Struct 2001;53:173–87. [22] Rebie‘re JL, Gamby D. A decomposition of the strain energy release rate associated with the initiation of transverse cracking, longitudinal cracking and delamination in cross-ply laminates. Compos Struct 2008;84:186–97. [23] Di Leonardo G. Fracture toughness characterization of materials under multiaxial loading. Int J Fract 1979;15:537–52. [24] Charalambides M, Kinloch AJ, Wang Y, Williams JG. On the analysis of mixedmode failure. Int J Fract 1992;54:269–91. [25] Hahn HT. A mixed-mode fracture criterion for composite materials. J Compos Technol Res 1983;5:26–9. [26] Rikards R, Buchholz FG, Wang H, Bledzki AK, Korjakin A, Richard HA. Investigation of mixed mode I/II interlaminar fracture toughness of laminated composites by using a CTS type specimen. Eng Fract Mech 1998;61:325–42. [27] Van Hoa S, Feng W. Hybrid finite element method for stress analysis of laminated composites. Dordrecht, The Netherlands: Kluwer Academic Publisher; 1998. [28] Inc A. Abaqus 6.5 user’s manual; 2005. [29] Hannaby SA. The solution of ordinary differential equations arising from stress transfer mechanics of general symmetric laminates. NPL Report CISE 13/97; 1997. [30] Varna J, Joffe R, Akshantala NV, Talreja R. Damage in composite laminates with off-axis plies. Compos Sci Technol 1999;59:2139–47. [31] Varna J, Joffe R, Talreja R. A synergistic damage-mechanics analysis of transverse cracking in [±h/904]s laminates. Compos Sci Technol 2001;61:657–65. [32] Varna J, Joffe R, Talreja R. Mixed micromechanics and continuum damage mechanics approach to transverse cracking in [S, 90n] laminates. Mech Compos Mater 2001;37(2):115–26.