Chinese Journal of Aeronautics, (2016), xxx(xx): xxx-xxx

Chinese Society of Aeronautics and Astronautics & Beihang University

Chinese Journal of Aeronautics

JOURNAL OF

AERONAUTICS

cja@buaa.edu.cn www.sciencedirect.com

Stress analysis and damage evolution in individual plies of notched composite laminates subjected to in-plane loads

Junshan Hua, Kaifu Zhang a*, Hui Chengb, Ping Liua, Peng Zoua, Danlong Song a

a The Ministry of Education Key Lab of Contemporary Design and Integrated Manufacturing Technology, Northwestern Polytechnical University, Xi'an 710072, China

b School of Mechanical Engineering, Northwestern Polytechnical University, Xi'an 710072, China Received 23 February 2016; revised 10 October 2016; accepted 10 October 2016

KEYWORDS

Damage evolution; Failure mode; Failure strength; Layer-by-layer stresses; User-defined subroutine

Abstract This work aims to investigate local stress distribution, damage evolution and failure of notched composite laminates under in-plane loads. An analytic method containing uniformed boundary equations using a complex variable approach is developed to present layer-by-layer stresses around the notch. The uniformed boundary equations established in series together with confor-mal mapping functions are capable of dealing with irregular boundary issues around the notch and at infinity. Stress results are employed to evaluate the damage initiation and propagation of notched composites by progressive damage analysis (PDA). A user-defined subroutine is developed in the finite element (FE) model based on coupling theories for mixed failure criteria and damage mechanics to efficiently investigate damage evolution as well as failure modes. Carbon/epoxy laminates with a stacking sequence of [45°/0°/—60°/90°]s are used to investigate surface strains, in-plane load capacity and microstructure of failure zones to provide analytic and FE methods with strong validation. Good agreement is observed between the analytic method, the FE model and experiments in terms of the stress (strain) distributions, damage evaluation and ultimate strength, and the layer-by-layer stress components vary according to a combination effect of fiber orientation and loading type, causing diverse failure modes in individuals.

© 2016 Production and hosting by Elsevier Ltd. on behalf of Chinese Society of Aeronautics and Astronautics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/

licenses/by-nc-nd/4.0/).

Corresponding author. E-mail address: zhangkf@nwpu.edu.cn (K. Zhang). Peer review under responsibility of Editorial Committee of CJA.

1. Introduction

Carbon fiber reinforced polymer (CFRP) composites are widely used in all fields of aerospace, automobile, electronic power, and mechanical engineering for several advantages offered over metals, ceramics, and plastics. These include low

http://dx.doi.org/10.1016/j.cja.2016.10.022

1000-9361 © 2016 Production and hosting by Elsevier Ltd. on behalf of Chinese Society of Aeronautics and Astronautics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

CJA 765

24 December 2016

2 J. Hu et al.

27 density, robust specific strength, low thermal expansion, corro-

28 sion resistance and designable characteristics for lightweight,

29 efficient structures. Cutouts, especially circular and elliptical

30 shapes, commonly appear in composite structural components

31 due to requirements for stability, maneuverability, low weight

32 optimization and accessibility of other systems. This is often

33 the case in aircraft where composite structures such as wing

34 spars and thin walls are drilled for electronic wires and hydrau-

35 lic pipes1 or to facilitate assembly operations.2 Other instances

36 occur not as part of the initial design of essential structures but

37 due to material defects or unexpected damage during a service

38 cycle.3 These are serious and dangerous cases since they are not

39 supposed to be there. Whatever their reasons for being, the

40 integrity and continuity of fiber and matrix in the composite

41 are destroyed; this makes the notch region the weakest part

42 of the structure and causes serious local stress concentrations

43 in the vicinity. As a consequence, structural capacity is reduced

44 and uncertain reactions to external loads may occur and lead

45 to unexpected fractures in service. In order to avoid potential

46 safety hazards, accurate analysis of local stress levels and dam-

47 age evolution of notched laminates are of great significance to

48 the utilization of the material and lay the foundation for engi-

49 neering applications.

50 The determination of stress fields in notched anisotropic

51 plates has been the focus of many scholars for a long time.

52 Muskhelishvili4 first introduced complex potential theory to

53 the isotropic elastic plates and successfully obtained an accu-

54 rate solution of stress distribution. Analytical solutions for

55 the stress distribution around holes of different shapes in ani-

56 sotropic plates were given by Lekhnitskii5 using series meth-

57 ods. Savin6 presented a much simpler approach by

58 conformal mapping of Cauchy integrals. Gao7 used a biaxial

59 loading factor together with an arbitrary orientation angle pre-

60 viously used to solve the problem of a plate with biaxial load-

61 ing at infinity, avoiding the superposition of a solution with

62 two uniaxial loading problems. The stress analysis research

63 above established the base for later study of notched compos-

64 ite laminates. Ukadgaonker et al. developed a general solution

65 for stress around oval holes,8 triangular holes9 with rounded

66 corners10 and irregular shaped holes11 in cross-ply and angle-

67 ply orthotropic plates under in-plane loading by superposition

68 of two stage solutions for boundary conditions. The stresses

69 were employed together with a series of failure criteria, which

70 were an extension of the Von Mises criterion for quadratic

71 interaction, to calculate the first ply failure (FPF) strength.

72 Similar solutions around a rectangular hole in an infinite iso-

73 tropic and anisotropic plate were given by Pan et al.12 and

74 Rao et al.,13 respectively, using the complex variable method.

75 Sharma14 suggested general solutions for determining the

76 stress distribution around polygonal holes and investigated

77 the effect of hole geometry and loading pattern on the stress

78 concentration factor. Batista15 used a modified solution to

79 solve problems of stress distribution around polygonal holes

80 of complex geometry in an infinite plate subjected to uniform

81 loading at infinity. Remeepazhand and Jafari16 also studied the

82 central polygonal hole problem in composite plates using two

83 simple equations with parameters k, c, n and w controlling the

84 size, shape and bluntness of corners. Toubal et al.17 experimen-

85 tally investigated the tensile strain field of composite plates in

86 the presence of stress concentration caused by geometrical cut-

87 outs consisting of circular holes by Electronic Speckle Pattern

Interferometer (ESPI). The stress obtained in experiments is 88

consistently lower than the analytical and numerical models. 89

This literature mostly concentrated on stress distribution in 90

anisotropic plates as well as laminated composites, which are 91

equivalent anisotropic plates, thus the stresses around cutouts 92

were simplified to be uniform throughout, regardless of thick- 93

ness. It should be noted that composite laminates contain sev- 94

eral plies which possess dissimilar properties due to different 95

fiber orientations, so the stress distribution in individual plies 96

should be thoroughly investigated. 97

The above stress research methodologies and conclusions 98

establish a theoretical basis for damage and failure of anisotro- 99

pic plates. As for failure analysis, the damage propagation pro- 100

cess is often depicted based on progressive damage analysis 101

(PDA) utilizing numerical and experimental methods.18 Lapc- 102

zyk and Hurtado19 proposed an anisotropic damage model 103

suitable for predicating failure and post-failure behavior in 104

fiber-reinforced materials. The plane stress formulation is used 105

and the response of the undamaged material is assumed to be 106

linearly elastic. The evaluation law is based on fracture energy 107

dissipation and implemented in a finite element code. Zahari 108

and EI-Zafrany20 developed a progressive analysis algorithm 109

based on Tsai-Hill failure to model the non-linear material 110

behavior and capture the compressive response of woven 111

glass/epoxy composite plates via non-linear finite element anal- 112

ysis. Rakesh et al.21 introduced a generic finite element model 113

to investigate the failure of unidirectional glass fiber reinforced 114

plastic (UD-GFRP) composite laminates with drilled holes 115

under tensile testing and compared the results with experimen- 116

tal work done earlier. Effects of joint geometry and stacking 117

sequence on the bearing strength and damage mode were 118

investigated by Ondurucu et al.22 and specimens were exam- 119

ined for failure modes using a scanning electron microscope 120

(SEM). Kim et al.23 employed rosette strain gauges to measure 121

the strain around the joint hole during insertion of stainless 122

steel pins into glass fiber reinforced plastic (GFRP) specimens, 123

and results were compared with the finite element method. Sat- 124

apathy et al.24 presented a modified fiber failure fatigue model 125

to characterize the behavior of laminated composites with a 126

central circular hole under in-plane fatigue loading. Martins 127

et al.25 studied the influence of diameter and thickness on fail- 128

ure pressure during tube burst tests and employed the progres- 129

sive failure analysis using a damage model by a user subroutine 130

(UMAT) implemented in ABAQUS software to understand 131

the behavior of composite tubes under internal pressure. Lee 132

et al.26 proposed an evaluation method for the progressive fail- 133

ure of composite laminates built upon Puck failure criterion by 134

implanted UMAT to efficiently analyze the progressive failure 135

phenomenon in glass/carbon fiber-reinforced composite lami- 136

nates. Compression tests were performed by Aljibori et al.27 137

on 16 fiber-glass laminated plates with and without circular 138

cut-outs to investigate the effects of varying the centrally 139

located circular cut-out sizes and fiber angle-ply orientations 140

on the ultimate load. Similar experiments were also conducted 141

by Abu et al.28 to investigate the influence of cut-outs on multi- 142

layer Kevlar-29/epoxy composite laminated plates. 143

These remarkable works contribute a lot to the failure anal- 144

ysis of notched anisotropic composites whose layups are quasi- 145

isotropic, cross-ply, and angle-ply. At the same time, numeri- 146

cal and advanced experimental methods for progressive failure 147

have greatly developed. However, the failure criteria employed 148

Individual plies of notched composite laminates

in the literature were linear expressions; whereas composite laminates in special application are not quasi-isotropic lay-up, composite materials accumulate damage in shear, leading to nonlinear stress-strain relations. The relationship between stress distribution and damage evolution thus needs further investigation.

In the present study, an analytical model for layer-by-layer stress around notches in individual plies of general composite laminates is investigated using a complex variable approach that contains uniformed boundary functions. The functions established in series together with conformal mapping functions are capable of dealing with irregular boundary issues both around the notch and at infinity. A stress solution is employed to evaluate the progressive damage and failure of notched laminates based on coupling theories for mixed failure criteria and damage mechanics. A user-defined subroutine is then developed in ABAQUS to simulate the damage evolution under tensile and compressive loading. Experiments for strain and failure are conducted using a 3D digital speckle measurement system and microscopy study to validate the proposed analysis method.

2. Analysis for local stress distribution

2.1. Problem configuration of notched laminates

A composite laminate with a traction free notch under inplane loading is investigated in Fig. 1. The inplane loading is in the form of stress resultants in an arbitrary coordinate system Xoy' rotated clockwise by an angle b from the global coordinate system xoy, as shown in Fig. 1(a). D is the diameter of the circle that encloses the notch. It can be written in terms of Gao's loading conditions7 using stress loading y as follows:

r1 = kr = 0

where r1, r1 and r1y are x' -direction, y'-direction and shear loading components, respectively, r is applied stress in y' direction, and k is biaxial loading factor defined by the ratio of r1 to ry1.

The notched plate is projected to a unit circle on f plane under normalized loading in Fig. 1(b) with an angle h to

denote hole position. The equivalent loading along xoy axes corresponding to original loading along x oy axes is obtained by transformation of axis:

~ = 2((k + 1) + (k - 1) cos2b)

«? = f ((k +1)

. = r ((k-

(k - 1) cos 2b) 1) sin 2b)

where r^0, r^0 and r^y are x-direction, y-direction and shear loading components along xoy axes, respectively. Through the variation of parameters k and b, the uniformed external loading form Eq. (2) covers several load conditions, namely uniaxial (k = 0), biaxial (k = 0, b—0) and shear stress (k = — 1, b = p/4 or 3p/4), as well as their combinations.

In states of generalized plane stress for laminates, according to Lekhnitskii,5 the generalized Hooke's law is valid when analysis is limited to the linear elastic range. Hooke's law in terms of mean stresses rx,y and mean strains exy is given by:

fix ~an a12 a16" <?x

£y = a12 an a26 ffy or Sx, y = aj ■ Ox, y (3)

Jxy. _a!6 a26 a66 rxy

where ex, sy and yxy are x-direction, y-direction and shear mean strain components, respectively. ax, ay and rxy are x-direction, y-direction and shear mean stress components, respectively. aij (i, j = 1, 2,6) are effective coefficients of the compliance matrix for symmetric multilayered composite plates.

In the absence of body force, the mean stress components <rx, y can be expressed in terms of Airy's stress function U(x, y) :

: = dy2

Sx2 d2

d2U dxdy

In order to ensure the integrity and continuity during deformation, the mean strains eX; y should meet the equilibrium equation:

d2ex d2Ëy

dy2 Sx2

By substituting Eqs. (3) and (4) into Eq. (7), a forth order differential equation is obtained:

ry = r

(a) A notched plate with arbitrary biaxial (b) The notched plate projected to a unite

loading on z plane circle on (plane with normalized loading

Fig. 1 Problem configuration and solution of a notched composite under in-plane loading.

J. Hu et al.

e4u „ a4u ,„ N eru

aii^J-2a2^ + (2an + a66)

1 ay4 =

dx4 dx3dy dx2dy2 dxdy3

The general solution of bi-harmonic Eq. (6) is determined by roots of the following characteristic equation:

«1114 - 2fli6l3 + (2fl!2 + «66)l2 - 2«26l + «22 — 0 (7)

In general, it can be proved that Eq. (7) has four distinct roots, which are always complex conjugates in pairs:

' 11 = ai + ib1

12 — a2 + i$2

Ju1 — a1 — ib1

12 — a2 — ib2

where a1, a2, b1, b2 are real numbers and b1 > 0, b2 > 0.11, i2 and their conjugates are complex parameters of a plane stress state. They are determined by effective compliance coefficients a] and can be used to describe the degree of material anisotropy.

Through the integral operation of Eq. (6), the stress function U(x, y) can be expressed by the real part of two complex functions F1(z1) and F2(z2):

U(x,y)=2Re(Fi(zi)+F2(z2)) zi,Z2 2 Sh (9)

where zp — x + ipy (p — 1,2) are anisotropic complex coordinates, and Sh represents the boundary conditions of the laminate in Fig. 1(a). So, the problem is reduced to finding two complex analytical stress functions that satisfy the boundary conditions at the notch edge as well as infinity.

2.2. Determination of analytic stress functions

Here, in the issue of plane stress state, the uniformed boundary equations u(z1) and W(z2) are developed in the form of derivatives of F1 and F2:

1 + 1 12 + 12 i( 12 - 12) B " ri "

(6) 1 1 0 B = 0.6r;°

.11 +11 12 + 12 i(12 - 12) . C si xy

U(z1)=f1 = C ln Z1 + Bz1 + £n=11 W(Z2) = f = C' lnz2 + (B + iC')z2 + E1=1 b

Z1, z2 2 Sh

where parameters y and y' are complex variables of the notch boundary. Parameters B, B and C are loading condition contents. Parameters an and bn are series coefficients determined by the boundary conditions Sh around the cutouts and at infinity.

2.2.1. Remote boundary conditions

At the infinite plate edge (|z| ! 1), the mean stress components reach the applied ^stress loading of the normalized boundary condition:

lim ffy = r°°

lim 1 = S1y

Through substituting Eqs. (9) and (11) into Eq. (4), the parameters B, B and C are obtained by the following equation group:

2.2.2. Boundary conditions at the notch edge At the notch edge, which is traction free, the notch area on the z plane is projected onto the exterior of a unit circle on f complex plane using the conformal mapping method (Fig. 1). The unit circle is defined by the Euler formula f — e'e in Fig. 1(b). This opens up the possibility of a simple general approach to determine the stress functions of the actual notch contour independently.

The mapping function for the elliptic notch in anisotropic cases from z plane to f plane is defined by:

= xp(0=§( (1 + m0( f +1) + i(1

- m1 )

1 - f K

P = 1; 2

where R — (a + b)/2 is a constant for size of the hole, and m1 — (a — b)/(a + b) is defined to control the ellipticity. Parameters a and b are major and minor axes of the elliptic hole, respectively. These are typically: (1) m1 — 0 or a — b: circular holes; (2) 0 < m1 < 1 or a—b: elliptical holes; (3) m1 ! 1 or a » b: blunted cracks.

By the projection process and taking Eq. (13) into the uniformed boundary Eq. (10) on f plane, the boundary equations in f plane are written as:

u(f) — u(z) — y ln(x1 (f)) + B(X1 (f)) + £ Akf—k

>(f) = W(z2)=C ln(x2(f)) + (B + iC )(«2(C))+E 1=1 Bk f-k

where parameters Ak, Bk (k — 1,2,...) are constants of stress functions. The resultant stress boundary conditions are determined by mean resultant forces XS in x-direction and YS in y-direction:

'2Re(i1u(f)+l2W(f))s = +fs*ds 2Re(u(f) + W(f))s = -Is Yds

It is seen that the boundary of the notch is free of traction, which indicates that the mean resultant forces on the notch are

X = 0 Y= 0

Parameters y and y are confirmed to be zero, and the logarithmic terms drop out. The rest of the parameters in series terms are obtained by substituting Eqs. (14) and (16) into Eq. (15):

_ (l2(K1+K2)-(K3+K4))+m1 (l2(K1+K2HK3+K4)) l1-l2

_ (11 (K1+K2)-(K3+K4))+m1 (11 (K1+K2)-(K3+K4 )) 11-12

,Ak = Bk = 0; k = 2; 3; 4 ...

where K1, K2, K3 and K4 are complex constants. They are determined by:

lim (Jx — r

x—>co

Individual plies of notched composite laminates

Ki — R (B(1 + ill) K2 — R (B(1 - ill) K3 = R (iiB(i + ill)

K4 = R(iiB(i - iii)

(B' + iC')(1 + il2)) (B' + iC')(1 - il2)) -l2(B'+ iC')(1 + il2)) -l2(B'+ iC')(1 - il2))

330 2.3. Local stresses around notches of individual plies

By substituting the analytic stress functions into Eq. (4), the mean stress components around the hole in global coordinates (Cartesian coordinates in Fig. 1) are obtained through the following equations:

rx — 2Re (1 u (f) + 1 w (f)

uy — 2Re

u'(f) .«1(f)

w'(f) «2©

sxy — -2Re ( 1 xTf 4

"12 W2®

According to the classical laminate theory29, laminate strains can be written in terms of the midplane strains s0xy and midplane curvatures Kxy as follows:

fix Kx

ey — 4 + z Ky

Jxy. / xy. Kxy

— ex ,y 4 zKx,y

where ex, ey and yxy are x-direction, y-direction and shear strain components in z thickness layer, and s0x, e0y and y0xy are x-direction, y-direction and shear strain components in the midplane. kx, Ky and Kxy are x-direction, y-direction and tor-sional curvatures in midplane, respectively.

For the symmetric laminates studied here, it can be proven that the coupling stiffness coefficients Bj (i,j — 1,2 ,6) are zero, which indicates the resultant and moment terms are uncoupled. Hence the curvature of midplane kxy is equal to zero. So, the strain at an arbitrary point of the notch edge throughout its thickness is the same as midplane strain and the mean strain:

8x, y — sx, y — ®x , y (21)

The stress components in the kth ply in principal direction are calculated by the following equation modified from Stress-Strain Equation and the stress transformation of axis:

rL 2 m2 n2 2mn Q11 Q12 Q16 7

iT — n2 2 m2 -2mn Q12 Q22 Q26

Tlt_ k - mn mn 22 m2 - n2 k .Q16 Q26 Q66_

a11 a12 a16 ffx

X a12 a22 a26 ffy

a16 a26 a66 sxy

where rL, rT and sLT are longitudinal, transverse and shear stress components in principal directions of the kth ply, respectively. m and n are parameters in stress transformed matrix Tk) decided by fiber orientation angle of the kth ply. Q(k) (i,j — 1,2 ,6) is the transformed reduced stiffness matrix for the kth layer. Eq. (22) can be written in contracted notation as:

3. Analysis for damage and failure 374

As mentioned before, one of the main purposes of determining 375

local stress distribution in individual plies around the notch is 376

to investigate the damage and failure of the perforated lami- 377

nate accurately. In progressive damage analysis (PDA) the- 378

ory,31 composites are materials characterized by complicated 379

failure phenomena that interact. The onset of failure in an 380

individual ply generally does not lead to the structure's col- 381

lapse, and this may still not occur until the failure has spread 382

to each of the multiple plies, indicating that some plies fail 383

before others. Since the notch region is the weakest part of a 384

composite plate, the damage initiates around the notch. 385

3.1. Failure criterion 386

2S2/G12 + 3aS4

2sLT/G12 + 3asLT 2Sl/G12 + 3aS4

2tLtI G12 + 3asLT lST/G^"30Sf

2tLt/G12 + 3asLT

IST/G^^Sf

Most failure criteria are extended from the von Mises criterion for quadratic interaction.32 They allow the failure of a whole

ply to be checked, but ignore the distinct failure modes of com- 389

posites. Given this, the mixed failure criteria, which combines 390

Hanshin criteria and maximum criteria, differs greatly from 391

them. It not only allows distinguishing between failure in ten- 392

sion compression as well as fiber and matrix, but also takes 393

nonlinear shear stress-strain behavior into consideration. The 394

mixed failure criteria33 are given by: emt = \jt) J

Lk)T — T<k). Qi

where emt, emc, eft, efc and efs are failure indices that denote 411

matrix tensile, matrix compression, fiber tensile, fiber compres- 412

sion and the fiber-matrix shear-out, respectively. The parame- 413

ter a is a material constant that is determined experimentally 414

according to Chang and Lessard34 and commonly takes the 415

value of 2.44 x 10~8 MPa~3 in the User's Manual of Abaqus 416

6.11. Once an index exceeds 1.0, the corresponding failure 417

occurs. The terms Xt, Xc are tensile and compressive strengths 418

in the fiber direction and Yt, Yc are tensile and compressive 419

strengths in the traverse direction. Sc is the in-plane shear 420

strength. 421

3.2. Material property degradation rules 422

From the phenomenological viewpoint, the defect initiation/ 423

propagation phenomenon for a material is strongly related 424

to the degradation of the material's capacity, and the most vis- 425

ible aspect is the degradation of the material's stiffness. There- 426

fore, it is possible to describe the progressive failure of a 427

CJA 765 24 December 2016 ARTICLE IN PRESS No. of Pages 14

6 J. Hu et al.

composite laminate by numerically representing the correlation between an increase in the material's internal damage and a decrease in the material's stiffness. The Camanho and Matthews degradation rules35 (Table 1) are adopted because they are closest to the real conditions according to the research of Chang.34 Parameters E1, E2, G12 and o12 are the longitudinal Young's modulus, transverse Young's modulus, shear modulus and Poison's ratios of the lamina. E1, E2, G'12 and o'12 are material properties after degradation.

3.3. Progressive failures of notched laminates

After using Eq. (23) to obtain the stress distribution around the notch in each individual ply, the mixed failure criteria is employed to check the failure state around the notch edge. On the occurrence of some sort of failure, the failing ply should be degraded according to Table ', which causes redistribution of stresses around the notch. Thus the failure state of the notch edge will be rechecked. The procedure of stress redistribution and failure recheck will continue until all plies fail or fiber failure occurs, causing sharp decline of carrying capacity. So, the location of damage initiation and final failure strength can be obtained during the analysis. Failure strength here refers to the maximum permissible remote stress loading applied.30

4. Numerical implementation

The stress distribution, damage and strength of notched specimens were also investigated using FE analysis implemented in the commercial software ABAQUS 6.11/standard. A composite laminate with a general stacking sequence of [45°/0°/—60° /90o]s was chosen as the investigated object, and the material properties are given in Table 2. The specimen configuration for the open hole tensile (OHT) test is illustrated in Fig. 2 according to ASTM D5766M-02a Standard.36 The open hole compressive (OHC) test followed ASTM D 6484/6484M-04.37 The specimen geometry was similar to that used for the OHC test. The aspect ratio met the condition W=D P 4 so that the external boundary effect was negligible.38

In the FE model, the 3D eight-node layered continuum shell element was employed for the laminate. The mesh was structured using linear elements with reduced-integration

(SC8R), and the notch region was highly refined by partition to capture the high near-hole stress gradients. A planar mesh size of 0.188 mm x 0.275 mm in (in x-y plane) in the vicinity of the hole was used while a coarser mesh was used away from the zone of the interest to computation cycles. In the thickness direction, a single layer of elements was used to model each ply. The approximate element size was 0.206 mm x 0.188 mm x 0.375 mm. Successive space discretization was compared, varying the element size surrounding the hole to verify that the final mesh size was enough to provide accurate results without expensive computational cost. A total of 92,345 eight-node brick elements were employed to model the specimen. The specimen was encastre on one side and applied tensile or compressive stress loading on the other side. Two lines, between which the length was L-2S, were marked to define the elongation corresponding to the part of the specimen that was not griped. The mesh and boundary conditions (BCs) are shown in Fig. 3.

The progressive model is implemented in ABAQUS by applying a separate subroutine called user-define field (USDFLD).39 In the analysis, the material properties depend on six field variables (FVs). The first five FVs from FV1 to FV5 represent failure indexes corresponding to emt, emc, eft, efc and efs. The last one, FV6, denotes the shear damage parameter. The algorithm is illustrated in Fig. 4. In the beginning, the FVs and the solution-dependent state variables (SDVs) are set to 0 and the material properties are equal to their initial values. Then the load is increased gradually; For each load increment, several iterations are necessary before the analysis converges to an equilibrium state. At the end of each increment, stresses, failure indices and shear damage parameter are computed at the integration points. The value of the shear damage parameter is directly assigned to the corresponding FVs and used to degrade the shear modulus. The values of the failure indices are stored as SDVs. Once an SDV exceeds 1.0, the corresponding FV is set to 1 and maintains this value until the end of the analysis. Meanwhile, the material properties are automatically reduced according to the degradation rules in Table 1. With each increasing applied load, the procedure is repeated until the program terminates due to excessive element distortion. Since equilibrium will not be reestablished after degradation of properties, small load increments should be set to obtain an accurate solution.40

5. Experimental investigation

The damage and failure behavior of notched composite laminates is also investigated on a Bairoe electronic universal testing machine with a load capacity of 100 kN. The tests were run in displacement control at a rate of 2 mm/min. Loading stopped after dropping 30% from the peak load. A non-contact measurement method, 3D digital speckle measurement system (DSMS), was used to investigate the strain field in the vicinity of the hole in the notched composite plate. The system is based on a 3D digital image correlation (DIC) system with binocular stereo vision technique that together perform full-field 3D measurements of surface strain and out-of-plane deformation of the joints, as shown in Fig. 5(a). The principle is to track an applied surface pattern during loading by continuously taking digital images of the surface.41 Two charge coupled device (CCD) cameras and two light emitting diode

Table 1 Property degradation Matthews. . rules of Camanho and

Failure mode Property degradation rules

Matrix tensile failure É2- = 0.2E2, m'12 ~ = 0, G'12 — 0.2G12

Matrix compression failure É2- = O.4E2, m'12 ~ = 0, G'X1 — 0.4G12

Fiber tensile failure Ei = = 0.07Ei , V'12 = 0

(4k) > 0)

Fiber compression failure Ei = = 0.14E , V'12 = 0

(4k) < 0)

Fiber-matrix shear-out failure m'12 " = 0; g'M = 0

(4k) < 0)

CJA 765 24 December 2016 ARTICLE IN PRESS No. of Pages 14

Individual plies of notched composite laminates 7

Table 2 Material properties of carbon-epoxy T300/3526.

E1 (GPa) E2 (GPa) »12 Gi2 (GPa) Xt (MPa) Xc (MPa) Yt (MPa) Yc (MPa) Sc (MPa)

133 10.4 0.29 4.14 1830 872 37.7 137 69.7

Fig. 2 Specimen configuration of composite for open hole test.

Fig. 3 Mesh and BCs of the FE model.

Fig. 4 Algorithm for ABAQUS user-defined subroutine.

526 (LED) sources are required for measurement. To ensure data

527 synchronization, DSMS was connected with the testing

528 machine by a synchronous data transmission device.

Fig. 5 Test machine with digital speckle measurement system and specimen processing.

The CFRP laminate was fabricated through hand lay-up 529

method using resin prepregs, and specimens were processed 530

with a water jet cutter to avoid defects caused by cutting heat. 531

The material properties and specimen geometry are the same 532

as the FE model in Section 5. Holes were drilled and reamed 533

CJA 765 24 December 2016 ARTICLE IN PRESS No. of Pages 14

8 J. Hu et al.

Fig. 6 Strain distribution around circular notch of composite laminates under tensile loading.

with carbide-tipped tools to attain fine surface quality, .then surfaces of the specimens were polished with sandpaper (first using a grain size of 800 then a grain size of 2000) to ensure a smooth surface. All the specimens were cleaned with alcohol to remove grease and dust before the test. A white background was painted on the testing surface of the specimens using spray paint, and later small black spray paint droplets were subtly applied with an appropriate density to make a high contrast-pattern, as shown in Fig. 5(b). Images of speckle pattern were recorded at a frequency of 5 Hz. The VIC 3D software calculated the accumulated strain by comparing the subsequent images with the first image, which was taken as a reference.

the opposite distribution law due to the tensile and compres-sive deformation state of specimen in either the x or y axial direction. For strains in the x axial direction, maximum strain occurs at the 90° point where by pass load goes through net tension plane and causes larger strain than any other part of the specimen. As for strains in the y axial direction, extreme points occur at regions around 0° and 90° due to transverse contracting under tension (Fig. 6(b))or transverse stretching under compression (Fig. 7(b)). Extreme points at 45° and 135° are caused by shear effect.

6.2. Stress distribution around the hole in individual layers

6. Results and discussion

6.1. Experimental verification of the strain distribution around the notch

The strain distributions around the notch on the surface of the specimen under tensile and compressive loading are presented in Figs. 6 and 7, respectively. The test was repeated three times, and the strains in the vicinity of the open hole were extracted by DSMS (Fig. 5(b)) to verify the validation of the theory and FE model. Since the strains around the notch are center symmetric, only those whose angular position h ranges from 0° to 180° are revealed.

The experimental results show relatively good trending alongside the proposed analytical solution and FE method with acceptable amplitude errors. There are two possible reasons for the errors: Either the open hole is a singularity at the surface of the specimen, or the fragmentary grid on the hole boundary cannot be used to calculate the strain field.42 In fact, the strains obtained by DSMS are the same as for the ring belt, which is a small offset from the hole edge. The strains fade much more sharply along the radial direction at the points where the strain concentrations are the most serious at the edge of the hole, so the obtained strains from ring belt near these points are apparently smaller than those from the hole edge.

The small error between analytical and FE results is due to the algorithm in the FE method, which is a numerical and approximate solution as compared to an analytical solution. Actually, the strains under tension and compression possess

The layer-by-layer stress distributions under tensile and com-pressive loading are provided separately in Figs. 8 and 9 to give clear insight into the mechanical behavior of notched composite laminates. Good agreement is observed between theory and FEA, thus providing confidence in the accuracy of the present results. Stress components in individual plies are shown to be highly non-uniform throughout the laminate thickness. The normalized longitudinal stresses in an individual layer have a much larger magnitude than the transverse and tangential stresses, reflecting the specific characteristic of anisotropic material in which carbon fibers bear most of the load.

It is seen that the stress concentration level in a layer depends on the relative angle between the fiber orientation and loading direction. The same dependence exists for concentration regions. For longitudinal stress, the smaller the angle is, the higher concentration level found near the region where the fiber direction is tangential to the hole perimeter. This can explain why the 0° layer exhibit the highest stress peaks of ±8.65 near 91.8° and is the most effective in carrying the axial load while the 90° layer has the smallest concentration level of ±2.25 at 0°. The transverse stresses in each layer exhibit peaks at locations where the fiber direction deviates most from the loading direction or where the fiber direction is normal to the notch boundary. So the 90° layer shows the largest concentration of ±0.62 near 90° while the 0° layer has a level of ±0.18 near 0°. Actually, there is adverse distribution rule between longitudinal and transverse stresses. For tangential stress influenced by the shear effect, higher degrees of concentration are found in layers whose fiber orientation possesses an angle in the loading direction.

Individual plies of notched composite laminates

Fig. 7 Strain distribution around circular notch of composite laminates under compressive loading.

i0° ™ 45°

Fig. 8 Stress distribution around circular notch in individual layers of composite laminates under tensile loading.

(a) Longitudinal stress (b) Transverse stress

— Theory x FEA Layer M 90° M -60° H0° «45°

Fig. 9 Stress distribution around circular notch in individual layers of composite laminates under compressive loading.

615 6.3. Damage initiation and propagation of notched composite

616 laminates

617 Damage initiation and propagation was evaluated by compar-

618 ing screenshots of various failure modes in different increment

619 steps (IS). Some of the important damage initiation locations

620 were also predicted using the analytical method in tion 3.3.

621 A microscopy study that investigated micro-scale failures

622 occurring in the vicinity of the hole using scanning electron

623 microscope (SEM) is shown in Fig. 12 to validate the damage

624 mode simulated in the FE model.

Fig. 10 shows the failure mode (i.e., matrix tensile, matrix 625

compress, fiber tensile, fiber compress, and fiber-matrix shear) 626

in different ply orientations under tensile loading. Matrix ten- 627

sile failure, which onset first at 92.72° in the 90° layer, was the 628

dominant failure mode in early increment steps. This damage, 629

though possibly not affecting tensile strength much, was 630

caused by transverse stress concentration and spread to other 631

layers with increasing loading. Fiber tensile failure occurred 632

at 88.71° in the 0° layer and was accompanied by matrix com- 633

pressive failure in later increment steps. As loading increased, 634

fiber damage propagated across the net tension plane causing 635

CJA 765 24 December 2016 ARTICLE IN PRESS No. of Pages 14

10 J. Hu et al.

dramatic decrease in load capacity. Matrix-fiber shear failure was observed due to excessive transverse shrinkage in late loading increments. The failure mode under compressive loading where matrix compressive failure initiation was observed at 88.4° in the 90° layer and propagates in early increment steps is presented in Fig. 11. Fiber compressive failure was found along with matrix-fiber shear in the 0° layer before specimen buckling deformation. The initiation location of fiber compres-sive damage occurred where fiber direction was normal to the hole perimeter. Tensile failure rarely appeared in the process.

Microscopy investigation of damaged specimen form tests was conducted to back up the progressive failure results from FE model. In the microscopic morphology, matrix crack (Fig. 12(a)) and fiber breakage (Fig. 12(b)), which corresponded to matrix tensile failure and fiber tensile failure from FE results in Fig. 10, were observed as common micro-scale failure modes in OHT specimens. Matrix crush (Fig. 12(c)), fiber buckling (Fig. 12(d)) and matrix-fiber shear (Fig. 12(e)) were dominant failure modes in compressive specimens, just as the FE model simulated in Fig. 11. Slight delamination is also found in the experiment because excessive load destroyed the compatibility of deformation among layers. This is beyond the elastic range and not considered in the analytical solution.

To further support the progressive failure results from the FE model, strain concentration and crack growth path on the specimen surface were captured by DSMS to compare with the damage of the 45° layer (surface layer) in the composites, and the results are listed in Figs. 13 and 14. For OHT testing, the initial damage points induced by matrix tensile failure (Fig. 13(c)) matched very well with the zones of strain concen-

tration (Fig. 13(a)), and the damage propagation path follows the crack profile. Fiber tensile failure (fiber tensile rupture) is also observed in both test results (Fig. 13(b)) and the FE results (Fig. 13(d)). As for OHC testing, the initial damage points of failure (Fig. 14(c)) were in accordance with the zones of strain concentration around the hole (Fig. 14(a)). The fibermatrix shear-out failure in the 45° layer (Fig. 14(d)) resulted in surface cracking of the specimen (Fig. 14(b)) because they possessed the same growth path. The above analysis provided sufficient credibility to the validation of the FE model for progressive failure analysis.

6.4. Failure strength

The load-displacement plots of experimental and FE simulation results for the notched composite laminates are shown in Figs. 15 and 16. The slope of the linearity in the experimental tests (i.e. elastic specimen stiffness) is seen to be very well captured by the numerical models, indicating that the composite can be regarded as linear-elastic material despite its brittle-ness. The stiffness of the specimen (slope of the curve) from the test is smaller than that from FE results due to the defects of the composite laminates and the degradation rules for the material properties. The OHC test shows slight nonlinear behavior at the start of the curve because the specimen slightly buckled under compressive loading. Since the OHC specimens were not ruptured, compressive failure is much more moderate than in the OHT test. The nonlinear behavior observed near the peak load is associated with severe damage in the specimens leading to fracture.

Fig. 10 Damage initiation and propagation of notched composite laminates under tensile loading.

Individual plies of notched composite laminates 11

Ply Matrix tensile Matrix compressive Fiber tensile Fiber compressive Shear-out

45° • • • • •

■ IS=68 ■ IS=30 ■ IS=68 «îsfllM is=6o sa^Viihafe is—68

0° A A A Vj1' 1 A

w w w w

IS—68 K Jl IS—30 1S—68 ■■i^H is—60 ■.'•.:■;■:!<!■ m is=&8

-60° • • • • •

■IS=68 UÊÊM 1S=30 ■ IS=68 ■Kfl is=60 x&smÊ^m is=68

90° • • • •

■ IS=68 ■ IS=30 ■ IS=68 ■ 1S=60 ■ IS=68 |

Fig. 11 Damage initiation and propagation of notched composite laminates under compressive loading.

(c) Fiber tensile (d) Fiber compression (e) Fiber shear

Fig. 12 Scanning electron microscope images of damage types in the vicinity of the notch.

694 In the tensile test, there is a protuberant point at a load of

695 2 47 MPa, which is mainly caused by the initiation of fiber ten-

696 sile damage around the notch in the 0° layer along with serious

697 matrix tensile failure in all plies. With each load increment, the

698 fiber tensile damage spreads across the tension plane at 0°

699 through thickness into 90°, —60°, + 45° layers in the vicinity

of the hole, resulting in drastic fluctuations in the curve. The 700

experimental ultimate strength was 484.6 MPa with an error 701

of 7.7% compared with the FE result of 524.9 MPa. In the 702

compressive test, the main damage mode in an earlier state 703

was matrix compressive failure. Fiber buckling failure occurs, 704

especially in the 0° layer in the net tension plan, causing a 705

CJA 765

24 December 2016

J. Hu et al.

Fiber tensile rupture

(a) Strain field of test

Fig. 13 Comparison of failure results between experimental test and FE model under tensile loading.

(a) Strain field of test (b) Crack growth patch of test (c) Matrix compressive (d) Fiber-matrix shear-out

failure of FEA tensile failure of FEA

Fig. 14 Comparison of failure results between experimental test and FE model under compressive loading.

" -»-FEA

Experiment

Mixture of severe fiber

- and matrix tensile failure

Damage initiation of fiber tensile failure

and propagation of matrix tensile failure

Axial displacemnt (mm)

Fig. 15 Load-displacement curve of notched composite laminates under tensile test.

706 quick shock in the load capacity curve. As the load increased,

707 buckling failure extended beyond the surroundings of the

708 notch and caused bending of the specimen. Since the specimens

709 were not ruptured, the compressive failure is much more mod-

710 erate than in the tensile test. There was a difference of 5.3%

711 between the test strength of 206.1 MPa and FE result of

712 217.5 MPa.

7. Conclusions 713

This research investigates the layer-by layer stress components, 714 and the initiation and propagation of damage and failure in 715 notched composite laminates by analytic method and FE 716

model, both of which are validated by experimental tests. 717 There is good agreement between the analytic method and 718

Individual plies of notched composite laminates

Fig. 16 Load-displacement curve of notched composite laminates under compressive test.

the FE model in terms of the stress (strain) distributions, damage initiation and location of failure modes, and experiments on strain by DSMS and microscopy study by SEM provide them with strong suitability and validation. The proposed analytic fundamentals present a dependable and manageable alternative to evaluate stress level and damage initiation quickly, while the FE model simulates the damage evolution process thoroughly. Several conclusions can be drawn from the research:

(1) The stress components are all non-uniform and periodic in different plies and caused by the anisotropy of the laminates. The magnitude of longitudinal stress is much larger than that of transverse stress and tangential stress due to the high elastic modulus of the fiber contrasting with the matrix.

(2) The layer-by-layer stress distributions around cutouts are affected by fiber orientation and loading direction. The high stress level of rL mainly occurs where the fiber direction is tangential to the hole perimeter. The concentration of rT is observed where the fiber direction is normal to the notch boundary. The sLT is more influenced by shear effect and possesses more extreme points than rL and rT.

(3) Matrix dominant failure is common in early states of inplane loading. It initiates at the hole edge in the ply whose fiber orientation is perpendicular to the load direction and then propagates to other plies. Fiber damage occurs and builds until matrix damage causes ruptures. While fiber buckling occurs in most plies, fiber breakage is found in the ply where the fiber orientation is parallel with loading direction.

(4) Load-displacement plots exhibit a long linear state, indicating that the composites can be regarded as linear-elastic material despite of their brittleness. Matrix damage evolution does not affect tensile strength much while fiber damage initiation may cause sudden shocks in plots. CFRP laminates possess better tensile load capacity than compressive ones.

Acknowledgements

This work was sponsored by the National Natural Science Foundation of China under Grant No. 51275410, No. 51305349 and No. 51305352. The authors would like to acknowledge the editors and the anonymous reviewers for their insightful comments.

References

1. Ostergaard MG, Ibbotson AR, Roux OL, Prior AM. Virtual testing of aircraft structures. CEAS Aeronaut J 2011;1(1-4):83-103.

2. Singh I, Bhatnagar N. Drilling-induced damage in uni-directional glass fiber reinforced plastic (UD-GFRP) composite laminate. Int J Adv Manuf Tech 2006;27(9-10):877-82.

3. Senthil K, Arockiarajan A, Palaninathan R, Santhosh B, Usha KM. Defects in composite structures: Its effects and prediction methods - a comprehensive review. Compos Struct 2013;106:139-49.

4. Muskhelishvili NI. Some basic problems of the mathematical theory of elasticity; fundamental equations, plane theory of elasticity, torsion, and bending, 1963. p. 123-9.

5. Lekhnitskii SG. Anisotropic plates. New York: Gordon and Breach; 1968. p. 41-5.

6. Savin GN. Stress concentration around holes. New York: Pergamon Press; 1961. p. 169-74.

7. Gao XL. A general solution of an infinite elastic plate with an elliptic hole under biaxial loading. Int J Pres Ves Pip 1996;67 (1):95-104.

8. Ukadgaonker VG, Awasare PJ, Gade SV. Stress analysis of door and window of a passenger aircraft. Aeronaut J 2000;104 (1039):409-14.

9. Ukadgaonker VG, Rao DKN. Stress distribution around triangular holes in anisotropic plates. Compos Struct 1999;45 (3):171-83.

10. Ukadgaonker VG, Rao DKN. A general solution for stresses around holes in symmetric laminates under inplane loading. Compos Struct 2000;49(3):339-54.

11. Ukadgaonker VG, Kakhandki V. Stress analysis for an ortho-tropic plate with an irregular shaped hole for different in-plane loading conditions - Part 1. Compos Struct 2005;70(3):225-74.

CJA 765 24 December 2016 ARTICLE IN PRESS No. of Pages 14

14 J. Hu et al.

12. Pan ZX, Cheng YS, Liu J. Stress analysis of a finite plate with a rectangular hole subjected to uniaxial tension using modified stress functions. Int J Mech Sci 2013;75:265-77.

13. Rao DKN, Babu MR, Reddy KRN, Sunil D. Stress around square and rectangular cutouts in symmetric laminates. Compos Struct 2010;92(12):2845-59.

14. Sharma DS. Stress distribution around polygonal holes. Int J Mech Sci 2012;65(1):115-24.

15. Batista M. On the stress concentration around a hole in an infinite plate subject to a uniform load at infinity. Int J Mech Sci 2011;53 (4):254-61.

16. Remeepazhand J, Jafari M. Stress analysis of perforated composite plates. Compos Struct 2005;71(3-4):463-8.

17. Toubal L, Karama M, Lorrain B. Stress concentration in a circular hole in composite plate. Compos Struct 2005;68(1):31-6.

18. Liu PF, Zheng JY. Recent developments on damage modeling and finite element analysis for composite laminates: a review. Mater Des 2010;31(8):3825-34.

19. Lapczyk I, Hurtado JA. Progressive damage modeling in fiber-reinforced materials. Compos Part A - Appl Sci 2007;38 (11):2333-41.

20. Zahari R, EI-Zafrany A. Progressive failure analysis of composite laminated stiffened plates using the finite strip method. Compos Struct 2009;87(1):63-70.

21. Rakesh PK, Singh I, Kumar D. Failure prediction in glass fiber reinforced plastics laminates with drilled hole under uni-axial loading. Mater Des 2010;31(6):3002-7.

22. Ondurucu A, Esendemir U, Tunay RF. Progressive failure analysis of glass-epoxy laminated composite pinned-joints. Mater Des 2012;36:617-25.

23. Kim SY, He B, Shim CS, Kim D. An experimental and numerical study on the interference-fit pin installation process for cross-ply glass fiber reinforced plastics (GFRP). Compos Part B-Eng 2013;54:153-62.

24. Satapathy MR, Vinayak BG, Jayaprakash K, Naik NK. Fatigue behavior of laminated composites with a circular hole under inplane multiaxial loading. Mater Des 2013;51:347-56.

25. Martins LAL, Bastian FL, Netto TA. Reviewing some design issues for filament wound composite tubes. Mater Des 2014;55:242-9.

26. 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. Compos Struct 2015;121:406-19.

27. Aljibori HSS, Chong WP, Mahlia TMI, Chong WT, Edi P, Alqrimli H, et al. Load-displacement behavior of glass fiber/epoxy composite plates with circular cut-outs subjected to compressive load. Mater Des 2010;31(1):466-74.

28. Talib AR Abu, Ramadhan AA, Rafie ASM, Zahari R. Influence of cut-out hole on multi-layer Kevlar-29/epoxy composite laminated plates. Mater Des 2013;43:89-98.

29. Jones RM. Mechanics of composite materials. 2nd ed. Philadelphia: Taylor & Francis; 1999. p. 109-21.

30. Muc A, Rys J, Latas W. Limit load-carrying capacity for spherical laminated shells under external-pressure. Compos Struct 1993;25 (1-4):295-303.

31. Zhang JY, Liu FR, Zhao LB, Chen YL, Fei BJ. A progressive damage analysis based characteristic length method for multi-bolt composite joints. Compos Struct 2014;108:915-23.

32. Pietropaoli E. Progressive failure analysis of composite structures using a constitutive material model (USERMAT) developed and implemented in ANSYS©. Appl Compos Mater 2012;19(3-4):657-68.

33. Ribeiro ML, Vandepitte D, Tita V. Damage model and progressive failure analyses for filament wound composite laminates. Appl Compos Mater 2013;20(5):975-92.

34. Chang FK, Lessard LB. Damage tolerance of laminated composites containing an open hole and subjected to compressive loadings. 1. Analysis. J Compos Mater 1991;25(1):2-43.

35. Camanho PP, Matthews FL. A progressive damage model for mechanically fastened joints in composite laminates. J Compos Mater 1999;33(24):2248-80.

36. ASTM International. D5766/D5766M-02a Standard test method for open hole tensile strength of polymer-matrix composite laminates; 2007.

37. ASTM International. D6484/D6484M-04 Standard test method for open-hole compressive strength of polymer-matrix composite laminates; 2004.

38. Wu HC, Mu B. On stress concentrations for isotropic/orthotropic plates and cylinders with a circular hole. Compos Part B-Eng 2003;34(2):127-34.

39. Hibbit, Karlssonand Sorensen Inc. ABAQUS/Standard version 6.11, user's manual; 2011.

40. Huang ZM. Inelastic and failure analysis of laminate structures by ABAQUS incorporated with a general constitutive relationship. J ReinfPlast Comp 2007;26(11):1135-81.

41. Ekh J, Schon J, Melin LG. Secondary bending in multi fastener, composite-to-aluminium single shear lap joints. Compos Part B-Eng 2005;36(3):195-208.

42. Vieille B, Aucher J, Taleb L. Overstress accommodation in notched woven-ply thermoplastic laminates at high-temperature: numerical modeling and validation by digital image correlation. Compos Part B-Eng 2013;45(1):290-302.

Hu Junshan is a Ph.D. student at School of Mechanical Engineering, Northwestern Polytechnical University, Xi'an, China. His main research interests are stress analysis, damage evolution and failure strength of composites.

Zhang Kaifu is a professor at School of Mechanical Engineering, Northwestern Polytechnical University, Xi'an, China. He received his Ph.D. from the same university in 2006. His research interests are the mechanical behavior of composite materials, including composite drilling and bolted composite interference-fit joints.