Composite Structures 121 (2015) 406–419

Contents lists available at ScienceDirect

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

Initial and progressive failure analyses for composite laminates using

Puck failure criterion and damage-coupled finite element method
Chi-Seung Lee, Jeong-Hyeon Kim, Seul-kee Kim, Dong-Man Ryu, Jae-Myung Lee ⇑
Department of Naval Architecture and Ocean Engineering, Pusan National University, Busan 609-735, Republic of Korea

a r t i c l e i n f o a b s t r a c t

Article history: In the present study, an evaluation method for the initial and progressive failure of composite laminates
Available online 18 November 2014 was proposed based on the Puck failure criterion and damage mechanics, respectively. In other words, the
initial failure (crack initiation in the fiber and/or matrix) and progressive failure (crack growth in the fiber
Keywords: and/or matrix) were evaluated using the Puck failure criterion, and fiber- and matrix-dependent damage
Composite laminates variables, respectively. In addition, the ABAQUS user-defined subroutine UMAT was developed based on
Puck failure criterion coupling theories for the failure criterion and damage mechanics in order to efficiently analyze the
Damage mechanics
progressive failure phenomenon in glass/carbon fiber-reinforced composite laminates. The developed
Progressive failure
ABAQUS user-defined subroutine
subroutine was applied to the failure of industrial composite laminates, and the analysis results were
compared to the experimental/numerical results previously reported in the literature. This comparative
study confirmed that the simulation results were in good agreement with the reported composite failure
results.
Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction composite laminates are postulated as single orthotropic materials.

For this reason, the local failure in a region of the fiber and/or
Owing to the many advantages of glass/carbon fiber-reinforced matrix under a particular stress combination cannot be calculated
composite laminates, including their light weight, robust specific [4,5].
stiffness and strength, and superior vibration, noise, and In order to overcome the aforementioned problem, Hashin
electromagnetic wave damping capacities, the application of such proposed a separate failure criterion for fiber and matrix collapse
composite laminates to industrial structures has been rapidly by introducing fracture plane-dependent stress components [4].
increasing during the last few decades. One application example And then, Puck and Schürmann enhanced Hashin’s failure criterion
for composite laminates involves subsea umbilicals, risers, and by implementing the angle of fracture plane and proposed three
flowlines (SURF). The reeling performance is the most crucial kinds of fracture modes: inter-fiber failures under tensile stress
capability of SURF. In other words, a low material stiffness and high on the plane perpendicular to the fiber direction and in-plane shear
tensile strength are in great demand during the installation and stress (mode A), failure due to a small compressive stress on the
operation of SURF. Therefore, composite laminates are more plane perpendicular to the fiber direction and large in-plane shear
suitable materials than metal-class materials [1]. In addition, stress (mode B), and failure due to a large compressive stress on
composite marine propellers for naval ships and submarines have the plane perpendicular to the fiber direction and small in-plane
recently been fabricated from carbon-reinforced composite shear stress (mode C) [5]. In particular, the Puck criterion has been
materials rather than Ni–Al–Bronze or Mn–Ni–Al–Bronze because verified by many researchers and engineers through experiments
they provide high fuel efficiency and low underwater radiated and numerical analyses [6–9].
noise [2,3]. However, it is quite arduous to apply the Puck criterion to
In these recent applications of composite laminates, only classic industrial fields because the theory is very complex, and there
failure criteria such as the maximum stress/strain, Tsai–Hill, and are no appropriate commercial finite element analysis (FEA) codes
Tsai–Wu criteria have been considered during the static analyses that contain the Puck criterion as a default option. Accordingly,
of composite structures. According to above failure criteria, the there have been few reports containing application examples of
the Puck criterion in the fields of naval architectural engineering
and aerospace engineering. Nevertheless, both the fiber failure
⇑ Corresponding author. Tel.: +82 51 510 2342; fax: +82 51 512 8836.
and inter-fiber failure should be estimated prior to the design
E-mail address: jaemlee@pusan.ac.kr (J.-M. Lee).

http://dx.doi.org/10.1016/j.compstruct.2014.11.011
0263-8223/Ó 2014 Elsevier Ltd. All rights reserved.
 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419 407

and fabrication of particular industrial structures such as wind tur- cracking, and fiber/matrix interface shear failure), along with their
bine blades [10]. associated damage evolution laws, which incorporated different
On the other hand, the initial fiber failure is the main design damage variables and conjugate forces. The numerical analysis
criterion to ensure the structural safety of a composite structure. results were compared to the amounts of internal pressure and
In other words, the failure of a single fiber represents the total radial displacement in a burst test of a composite cylindrical
failure of the structure. However, this design criterion is too storage vessel.
redundant and conservative, and it is laborious to describe the In addition, there have been a few similar studies on simulation
progressive failure of the fiber and/or matrix under increasing methods and damage models for the progressive failure of
static loads or cyclic fatigue loads. composite laminates, e.g., a failure model based on a generalized
In order to remove the aforementioned obstacle, tensile, laminate plate theory [14], an LaRC04 failure criterion-orthotropic
compressive, shear, and fatigue tests of composite laminates are damage model that took into account the crack closure effect
commonly carried out to investigate the failure features from the [15,16], a fully 3D continuum damage model [17], a micro-
initial to final failure. Although such material tests are considered mechanical model based on a unit cell [18], and a Hashin failure
to be reliable and straightforward, a huge investment in time and criterion-material properties degradation model [19].
cost are needed to establish the experimental facilities. In these studies, the onset and growth of composite internal
As an alternative approach, a damage mechanics-coupled FEA defects were successfully analyzed, but there were two kinds of
method is widely adopted nowadays. Based on damage mechanics, challenging factors for industrial utilization. First, the assessment
it is feasible to quantitatively describe the relationship between method was conducted in the framework of in-house FEA code
the increase in a material’s internal defects and the decrease in rather than commercial FEA code such as ABAQUS. Second, the
its stiffness. Although it requires a somewhat substantial amount initial failure was estimated using only the criterion of Hashin
of computational time, it is possible to simulate the progressive rather than that of Puck. However, many classification societies
damage to a composite structure using the simple material proper- such as Germanischer Lloyd specify the use of the Puck criterion
ties of a composite laminate.
There is a large body of literature on the progressive failure
analysis of composites, much of which exploits the damage Table 3
mechanics and material stiffness degradation methods. Chang EWM conditions for each failure mode.
and Chang [11] carried out two kinds of composite laminate Failure mode Condition for validity Damage value
analyses. They assessed the stress and strain on the basis of the
Tensile fiber mode (1) dft = 1
classical lamination theory and calculated the damage accumula- Compressive fiber mode (2) dfc = 1
tions using a material degradation-implemented failure criterion Tensile matrix mode (3) dmt = 1
that they proposed. Their research confirmed that the results of a Compressive matrix mode (4) and (5) dmc = 1
numerical analysis of the initial and progressive failures of cen-
ter-notched composite laminates under tensile loads coincided
well with those of experiments. Tay et al. [12] proposed an element
failure method to effectively express the damage growth in a
framework of finite elements. Progressive failures were detected
using various types of failure models such as the Tsai–Wu and
micromechanics-based failure criteria, and the results of numerical
analyses were compared to the tensile test results for double-
notched composite laminates. Liu and Zheng [13] developed an
energy-based damage model to predict the damage/flaw propaga-
tion in three-dimensional (3D) composite cylindrical laminates,
which were fabricated using a filament winding method. Three
kinds of failure modes were proposed (i.e., fiber breakage, matrix

Table 1
Application conditions for Puck failure criterion.

Equation Condition for validity

(1) e1 þ mEff121 mrf r22 P 0
(2) e1 þ mEff121 mrf r22 < 0
(3) r22 P 0
 
(4) r22 < 0 and 0 6 rs2122  6 jsR21c
A
vv
j
 
(5) r22 < 0 and 0 6 rs2122  6 jsR21c
A
j
vv

Table 2
Recommended values for inclination parameters.

Parameter Value
GFRP/Epoxy CFRP/Epoxy

pþ
vp 0.30 0.35
p
vp 0.25 0.30

pvv 0.20–0.25 0.25–0.30
Fig. 1. Algorithm for ABAQUS user-defined subroutine.
408 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419

Table 4 composite structure failure analysis. Finally, in order to validate

Material properties of E-glass 21xK43/LY556 and AS4/3501-6 test specimens. the proposed analysis method, the analysis results were compared
Property E-glass 21xK43/LY556 AS4/3501-6 to the failure test results for composite laminates reported in the
E11 (MPa) 53,480 126,000 literature.
E22 (MPa) 17,700 11,000
G12 (MPa) 5830 6600 2. Theories for initial and progressive failure of composite
m12 0.278 0.28
laminates
Ef1 (MPa) 80,000 225,000
mf12 0.2 0.2
e1T (%) 2.132 1.380 2.1. Initial failure theory based on Puck failure criterion
e1C (%) 1.065 1.175
XT (MPa) 1140 1950
Hashin pointed out the limitation of classical failure criteria
XC (MPa) 570 1480
YT (MPa) 35 48
such as the maximum stress/strain, Tsai–Hill, and Tsai–Wu criteria,
YC (MPa) 114 200 and proposed four kinds of unique failure functions [4]. Unlike the
S12 (MPa) 72 79 classical failure criteria, in the Hashin failure function, the
equations are divided into the tensile and compressive parts, and
more complex terms are introduced for matrix failure under
Table 5 compressive loads.
Thicknesses of E-glass 21xK43/LY556 and AS4/3501-6 test specimens. After two decades, Puck and Schürmann modified and enhanced
the Hashin failure criterion, i.e., the fiber failure was calculated
E-glass 21xK43/LY556 AS4/3501-6
based on the material properties of the fiber itself rather than those
Ply Ply thickness Ply orientation Ply thickness
of the ply, and the matrix failure was divided into three
orientation (°) (mm) (°) (mm)
different stress states, including the transverse tension (mode A),
90 0.172 90 0.1375
moderate transverse compression (mode B), and large transverse
+30 0.414 +45 0.1375
30 0.414 45 0.1375
compression (mode C). In addition, an equation was proposed for
– – 0 0.1375 determining the angle of the new fracture plane that was
Total thickness (mm) 2.0 1.1 generated during the fracture of a composite laminate.
The Puck failure criterion can be written as follows.

for the identification of the initial failure, as discussed above. For

(1) Fiber failure in tension
this reason, it is extremely burdensome to employ the aforemen-  
tioned methodologies in industrial applications. 1 mf 12
e1 þ mrf r22 ¼1 ð1Þ
Hence, in the present study, the initiation and progression of e1T Ef 1
the internal damage to a composite laminate were practically ana-
lyzed using the Puck failure criterion, as well as a damage mechan-
(2) Fiber failure in compression
ics-based material properties degradation model. In addition, the  
introduced model was implemented in the ABAQUS user-defined 1  mf 12 
 þ ð10c Þ2 ¼ 1
 e1 þ mr f r 22  21 ð2Þ
subroutine UMAT to facilitate the accessible utilization of a e1C E f1

(a) Ply stacking sequence [90/±30]S (b) 90° fiber angle direction at 1st ply

(c) +30° fiber angle direction at 2nd ply (d) -30° fiber angle direction at 3rd ply

Fig. 2. Ply stacking sequence and fiber angle direction for E-glass/LY556 composite laminate.
 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419 409

(3) Matrix failure in transverse tension The application conditions for Eqs. (1)–(5) are listed in Table 1.
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
 2  2  2 In this table, RAvv is the fracture resistance of the action plane
s21 YT r22 r22 r11 against its fracture due to transverse/transverse shear stressing,
þ 1  pþv p þ pþv p þ ¼1 ð3Þ
S21 S21 YT S21 r11D and s21c is the shear stress at the ‘‘turning point’’ of the (r22, s21)
fracture curve.
On the other hand, the main parameters for the Puck failure
(4) Matrix failure in moderate transverse compression A
criterion such as pþ  
v p , pv p , pvv , Rvv , and s21c can be obtained by using
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 
1 r11 the following formulae.
s221 þ ðpv p r22 Þ2 þ pv p r22 þ ¼1 ð4Þ  
S21 r11D ds21
pþv p ¼  of ðr22 ; s21 Þcur v e; r22 P 0 ð6Þ
dr22 r2 ¼0
(5) Matrix failure in large transverse compression
 
" 2  2 # ds21
s21 r22 YC r11 pv p ¼  of ðr22 ; s21 Þcur v e; r22 6 0 ð7Þ
þ þ ¼1 ð5Þ dr22 r2 ¼0
2ð1 þ pvv ÞS21 YC ðr22 Þ r11D
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !
where e1T and e1C are the tensile and compressive failure A Y S21 YC
strains of a unidirectional layer in the x1 direction, respec- Rvv ¼ ¼ 1 þ 2pv p 1 ð8Þ
2ð1 þ pvv Þ 2pv p S21
tively; e1 is the normal strain of a unidirectional layer; mf12
and Ef1 are Poisson’s ratio and Young’s modulus for the fiber
RA
in the x1 direction, respectively; mrf is the mean stress mag- pvv ¼ pv p vv ð9Þ
nification factor for the fibers in the x2 direction; r11 and r22 S21
are the normal stresses in a unidirectional layer; c21 and s21 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
are the shear strain and stress of a unidirectional layer in the s21c ¼ S21 1 þ 2pvv ð10Þ
elastic symmetry direction, respectively; S21 is the shear
strength of a unidirectional layer transverse and parallel to However, it is arduous to obtain some parameters such as pþ 
v p , pv p ,
the fiber direction; pþ   and p because the (r , s ) curve should be obtained using a ser-
v p , pv p , and pvv are the fracture plane vv 22 21

angle-dependent parameters; and r11D is the stress value ies of experiments. Accordingly, Puck et al. [20] recommended
for linear degradation. using the values listed in Table 2.

2.2. Progressive failure theory based on damage mechanics

Table 6
Loading cases for initial failure analysis. From the phenomenological viewpoint, the defect initiation/prop-
agation phenomenon for a material is strongly related to the degrada-
Case Laminate Material Loading case
layup
tion of the material’s capacity, and the most visible aspect is the
deterioration of the material’s stiffness. Therefore, it is possible to
A [90/±30]S E-glass/LY556 Biaxial failure stress envelope (ry vs rx)
describe the progressive failure of a composite laminate by numeri-
B Biaxial failure stress envelope (rx vs sxy)
C [90/±45/0]S AS4/3501-6 Biaxial failure stress envelope (ry vs rx) cally representing the correlation between an increase in the mate-
rial’s internal damage and a decrease in the material’s stiffness.
There are numerous damage models for identifying the internal
damage for various materials. Toi and Lee [21] and Lee et al. [22]

Table 7
Dimensions and ply orientations of test specimen.

Laminate layup D (mm) W (mm) T (mm) L (mm)

[(0/90)6]S 6.35 25.4 3.175 203.2
[(0/±45/90)3]S
[(±45)6]S

Fig. 3. Example of loading and boundary conditions for case B.

Table 8
Material properties of T300/1304-C test specimen.

Property T300/1304-C
E11 (MPa) 146,858
E22 (MPa) 11,376
G12 (MPa) 6185
m12 0.30
Ef1 (MPa) 230,000
mf12 0.2
e1T (%) 1.807
e1C (%) 0.652
XT (MPa) 1731
XC (MPa) 1379
YT (MPa) 67
YC (MPa) 268
S12 (MPa) 134
Fig. 4. Configurations of T300/1304-C test specimen.
410 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419

(a) Ply stacking sequence [(0/±45/90)3]S

(b) 0° fiber angle direction at 1st ply

(c) +45° fiber angle direction at 2nd ply

(d) -45° fiber angle direction at 3rd ply

(e) 90° fiber angle direction at 4th ply

Fig. 5. Ply stacking sequence and fiber angle direction for T300/1034-C composite laminate.

Fig. 6. Loading and boundary conditions for T300/1034-C composite laminate.

adopted the Lemaitre damage model to solve the thermal-induced the debonding failure problem of adhesively bonded joints. Kim
failure and fatigue failure problems of hot-dip galvanized struc- et al. [24,25], Lee et al. [26] and Kim et al. [27] used the isotropic
tural members and welded structures, respectively. Lee et al. [23] Bodner–Chan damage model to describe the crack propagation
employed the Gurson–Tvergaard void growth model to identify problems as well as the damage-induced viscoplastic behavior of
 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419 411

austenitic stainless steel structural members. Lee and Lee [28] C 012 ¼ E011 ðm21 þ m31 m23 ÞD ð26Þ
applied the anisotropic Bodner–Chan damage model to calculate
the progressive failure of a glass-fiber-reinforced polyurethane C 023 ¼ E022 ðm32 þ m12 m31 ÞD ð27Þ
foam-based liquefied natural gas insulation structure under a cyc-
lic impact load, as well as a unidirectional static load. C 013 ¼ E011 ðm31 þ m21 m32 ÞD ð28Þ
Although these damage models are well defined and widely
adopted, both the numerical derivation and implementation D ¼ 1=ð1  m12 m21  m23 m32  m31 m13  2m21 m32 m13 Þ ð29Þ
procedure for finite elements are considerably complicated. There
are two reasons. (1) Most of the damage models were designed for where Cij and C 0ij
are the initial and damaged material stiffness ten-
isotropic ductile materials. Thus, it is difficult to apply a general sors, respectively; df and dm are the total damage variables for the
damage model to anisotropic brittle materials. (2) Much effort is fiber and matrix, respectively; dft, dfc, dmt, and dmc are the fiber
required to identify the material parameters of a damage model, and matrix damage variables in relation to the tensile and compres-
i.e., an enormous number of material test results should be sive stress states, respectively; mij is the Poisson ratio tensor; and smt
obtained to define the exact material parameters. and smc are the loss control factors for the shear stiffness caused by
Consequently, in the present study, stress-based Puck failure the matrix tensile and compressive failures, respectively. In the
criterion were adopted to predict the evolution of damage. Namely, present study, the loss control factors were postulated as smt = 0.9
the damage initiation and growth were detected based on Puck’s and smt = 0.5.
criteria and material stiffness degradation method, respectively.
This approach is relatively simple and extensively used by compos- 3. Computational analysis procedures
ite researchers for strength estimation, as well as for progressive
failure analysis, and there have been a few reports of good agree- 3.1. Initial/progressive failure analysis procedure using ABAQUS user-
ment with the results of experiments, such as McCarthy et al. defined subroutine
[29] and Tay et al. [12].
The elastic stress–strain constitutive relation, fiber and matrix The stresses and strains for each element, as well as the mate-
damage variable with respect to the tensile and compressive stress rial properties, were implemented into the Puck failure criterion
states, and material constants for the initial and damage-coupled in order to predict the initial failures of composite laminates under
material stiffness for orthotropic composite laminates can be arbitrary loads.
written as follows. Once the failure was disclosed, the material stiffness of the
8
>
>
r11 9
>
>
2
C 11 C 12 C 13 0 0 0
38
>
>
9
e11 >
>
initially failed element was replaced by zero in compliance with
>
> r22 >
> 6C >
> e22 >
7> > Eqs. (14)–(22). This method is the so-called ‘‘element weakening
>
> >
> 6 12 C 22 C 23 0 0 0 7> >
>
>
> >
> 6 7>
> >
> method (EWM)’’ in the finite element (FE) procedure [26,28]. The
< r33 = 6 C 13 C 23 C 33 0 0 0 7< e33 = application conditions for the EWM for the tensile and compressive
¼6
6 0
7
7> e > ð11Þ
>
>
> r12 >>
> 6 0 0 2G12 0 0 7>
> 12 >>
fiber and matrix failure modes are listed in Table 3.
>
> >
> 6 7> >
> r > 4 0 5>
> >
> The aforementioned theory for the Puck failure criterion and the
>
>
:
23 >
>
;
0 0 0 2G23 0 > e23 >
>
: >
; material stiffness degradation model were implemented in the
r31 0 0 0 0 0 2G31 e31
commercial FEA code (ABAQUS) user-defined subroutine UMAT
in the present study. The computational algorithm that includes
df ¼ 1  ð1  dft Þð1  dfc Þ ð12Þ
the EWM is shown in Fig. 1.
dm ¼ 1  ð1  dmt Þð1  dmc Þ ð13Þ
Table 9
C 11 ¼ ð1  df ÞC 011 ð14Þ Analysis cases for initial failure analyses.

Case Loading condition Initial failure

C 22 ¼ ð1  df Þð1  dm ÞC 022 ð15Þ A-1 rx = 25 MPa, ry = 25 MPa Xa
A-2 rx = 100 MPa, ry = 100 MPa Ob
C 33 ¼ ð1  df Þð1  dm ÞC 033 ð16Þ A-3 rx = 150 MPa, ry = 150 MPa X
A-4 rx = 250 MPa, ry = 250 MPa O
A-5 rx = 10 MPa, ry = 50 MPa Cc
C 12 ¼ ð1  df Þð1  dm ÞC 012 ð17Þ A-6 rx = 100 MPa, ry = 50 MPa O
A-7 rx = 200 MPa, ry = 10 MPa X
A-8 rx = 250 MPa, ry = 10 MPa O
C 23 ¼ ð1  df Þð1  dm ÞC 023 ð18Þ
B-1 rx = 10 MPa, sxy = 30 MPa X
B-2 rx = 10 MPa, sxy = 60 MPa O
C 13 ¼ ð1  df Þð1  dm ÞC 013 ð19Þ B-3 rx = 50 MPa, sxy = 40 MPa X
B-4 rx = 50 MPa, sxy = 80 MPa O
B-5 rx = 150 MPa, sxy = 40 MPa X
G12 ¼ ð1  df Þð1  smt dmt Þð1  smc dmc ÞG012 ð20Þ B-6 rx = 150 MPa, sxy = 80 MPa O
B-7 rx = 200 MPa, sxy = 40 MPa C
G23 ¼ ð1  df Þð1  smt dmt Þð1  smc dmc ÞG023 ð21Þ B-8 rx = 200 MPa, sxy = 80 MPa O
C-1 rx = 100 MPa, ry = 100 MPa X
C-2 rx = 300 MPa, ry = 300 MPa O
G31 ¼ ð1  df Þð1  smt dmt Þð1  smc dmc ÞG031 ð22Þ C-3 rx = 600 MPa, ry = 600 MPa X
C-4 rx = 900 MPa, ry = 900 MPa O
C-5 rx = 100 MPa, ry = 200 MPa X
C 011 ¼ E011 ð1  m23 m32 ÞD ð23Þ C-6 rx = 300 MPa, ry = 200 MPa O
C-7 rx = 200 MPa, ry = 100 MPa X
C 022 ¼ E022 ð1  m13 m31 ÞD ð24Þ C-8 rx = 600 MPa, ry = 100 MPa O
A
Not failed.
C 033 ¼ E033 ð1  m12 m21 ÞD
B
ð25Þ Failed.
c
Almost failed.
412 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419

(a) Case A (b) Case B

(c) Case C
Fig. 7. Initial failure stress envelopes and analysis cases.

As shown in this figure, two types of analyses were carried out: There are numerous experimental results available for compos-
an initial failure analysis and a complete failure analysis. During ite laminate failures under arbitrary loads. However, in the present
the former analysis, the stresses and strains of each element were study, the results presented by Soden’s research team [6–8], and
calculated using the constitutive relation (Eq. (11)) and its associ- Chang’s research team [11,30] were selected for validations of fail-
ated engineering constants (Eqs. (23)–(29)). In addition, the initial ure stress envelope and damage propagation, respectively.
failure was estimated based on the Puck failure criterion (Eqs. (1)–
(5)). If a failure was detected in any single element, the EWM was 3.2.1. Composite laminates for failure stress envelope validation
applied to the damaged element. Then, the material stiffness was The composite laminates for validation of failure stress were
modified using Eqs. (12)–(22), and the stresses and strains were fabricated from E-glass 21xK43 glass fiber and an LY556 epoxy
calculated repeatedly until the regions of interest completely matrix, and AS4 carbon fiber and a 3501-6 epoxy matrix. The lam-
failed. Finally, it was possible to obtain the stress, strain, and dam- inate layups for the glass-fiber-reinforced polymer (GFRP) and car-
age contours of the composite laminate finite elements, and we bon-fiber-reinforced polymer (CFRP) were [90/±30]S and [90/±45/
could predict the structural failure under given loading and bound- 0]S, respectively. The material properties and ply thickness for each
ary conditions. material are listed in Tables 4 and 5, respectively.
Although a parameter used in Eqs. (5) and (6), mrf (i.e., the mag-
3.2. Target materials for computational analyses nification factor of the mean stress for fiber in the x2 material
direction), could be determined by experiment, simple values for
In the present study, two types of comparison approaches were this parameter (mrf = 1.3 for GFRP and mrf = 1.1 for CFRP) were
introduced to verify the proposed computational analysis method, recommended by Puck and Schürmann [5].
i.e., validation through (1) failure stress envelope, and (2) damage Fig. 2 illustrates the total stacking sequence and fiber direction
propagation simulation. Hence, two different kinds of target mate- at each ply of the E-glass/LY556 laminate in the ABAQUS FE model.
rials were used for both validations. The hexahedral continuum shell element, i.e., SC8R in ABAQUS,
 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419 413

(a) Fiber failurecontour (case A-2) (b) Matrix failurecontour (case A-2)

(c) Fiber failure contour (case B-2) (d) Matrix failure contour (case B-2)

(e) Fiber failure contour (case C-2) (f) Matrix failure contour (case C-2)
Fig. 8. Fiber/matrix failure contour based on Puck failure criterion for case A-2, B-2 and C-2.

Table 10 Table 11
Initial failure analysis results and agreement with experimental results. Analysis cases for progressive failure analyses.

Case PFFIa PMFIb Failure Agreementc Case Loading condition Failure state
A-1 0.0371 0.4363 X O B-1 rx = 10 MPa, sxy = 30 MPa No failure
A-2 0.1483 1.7450 O O B-2 rx = 10 MPa, sxy = 60 MPa Initiation
A-3 0.4454 0.9012 X O B-2–1 rx = 10 MPa, sxy = 150 MPa Progression
A-4 0.7424 1.5020 O O B-2–2 rx = 10 MPa, sxy = 320 MPa Final
A-5 0.2222 1.5900 C O B-7 rx = 200 MPa, sxy = 40 MPa No failure
A-6 0.3327 1.5740 O O B-8 rx = 200 MPa, sxy = 80 MPa Initiation
A-7 0.4113 0.9386 X O B-8–1 rx = 200 MPa, sxy = 125 MPa Progression
A-8 0.5116 1.1720 O O B-8–2 rx = 200 MPa, sxy = 200 MPa Final
B-1 0.0661 0.7311 X O
B-2 0.1384 1.8240 O O
B-3 0.1674 3.4950 X X
B-4 0.2333 4.3120 O O
B-5 0.3727 0.8922 X O and progressive failures of a test specimen under an in-plane stress
B-6 0.4379 9.7880 O O state [6,7].
B-7 0.4767 1.0430 C O
The loading cases for the initial/progressive failure analyses are
B-8 0.5415 1.6640 O O
C-1 0.1018 0.4098 X O
listed in Table 6. During the analysis, the in-plane stress states (ry
C-2 0.3054 1.2290 O O vs rx and rx vs sxy) were considered because of the limitation of the
C-3 0.7174 0.5901 X O experimental results of the research by Soden et al. [6,7].
C-4 1.0760 0.8852 O O In addition, an example of the loading and boundary conditions
C-5 0.3804 0.5977 X O
for case B is illustrated in Fig. 3. As shown in this figure, two edges
C-6 0.5059 1.2550 O O
C-7 0.3804 0.5977 X O of the FE model are simply supported, namely, UX = 0 and UY = 0,
C-8 1.0470 1.3040 O O and the other edges are incrementally loaded by normal and shear
Accuracy of analysis results (%) 95.8 stresses.
a
Puck fiber failure index.
b
Puck matrix failure index. 3.2.2. Composite laminates for damage propagation validation
c
Agreement between Puck failure analysis results and experimental results. The composite laminates for validation of damage propagation
were made of T300/1034-C graphite/epoxy prepreg tapes. The con-
was adopted for the FE modeling, and the length and width of the figurations and the dimensions/ply orientations of the test speci-
test specimen were postulated to each be 10 mm because these men are shown in Fig. 4 and Table 7, respectively, and the
two dimensions were not important when investigating the initial material properties of T300/1034-C are listed in Table 8.
414 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419

Fig. 5 demonstrates the total stacking sequence and fiber direc-

tion at each ply of the T300/1034-C composite laminate in the
ABAQUS FE model. The element type is equal to the case of failure
stress envelope validation, namely, the hexahedral continuum
shell element, i.e., SC8R in ABAQUS, was selected for the FE model-
ing. In order to reduce the computational time, 1/4 model was
adopted during analysis.
The boundary and loading conditions are indicated on the FE
model of composite laminate as shown in Fig. 6. As shown in this
figure, left and bottom edge of the FE model are X- and Y-symmetry
conditions, respectively, and right edge is displacement-controlled
to X-direction.

4. Computational analysis results

4.1. Outline of computational analysis

As aforementioned in Section 3, two types of comparison

approaches were adopted in order to validate the developed
analysis method. Namely, initial and final failures of composite
laminate plates under complex loads were evaluated by comparing
with the failure stress envelope for test results of composite
Fig. 9. Initial and final failure stress envelopes and analysis cases. laminates. In addition, damage growth of laminated composites

(a) Fiber failure contour (case B-1) (b) Fiber failure contour (case B-2)

(c) Fiber failure contour (case B-2-1) (d) Fiber failure contour (case B-2-2)

(e) Matrix failure contour (case B-1) (f) Matrix failure contour (case B-2)

(g) Matrix failure contour (case B-2-1) (h) Matrix failure contour (case B-2-2)
Fig. 10. Progressive fiber/matrix failure contours based on Puck failure criterion for case B-1/2/2–1/2–2.
 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419 415

(a) Fiber failure contour (case B-7) (b) Fiber failure contour (case B-8)

(c) Fiber failure contour (case B-8-1) (d) Fiber failure contour (case B-8-2)

(e) Matrix failure contour (case B-7) (f) Matrix failure contour (case B-8)

(g) Matrix failure contour (case B-8-1) (h) Matrix failure contour (case B-8-2)
Fig. 11. Progressive fiber/matrix failure contours based on Puck failure criterion for case B-7/8/8–1/8–2.

containing stress concentration under axial loads were predicted

by comparing with the tensile test results of composite plies con- Table 12
Failed elements from no-failure to final-failure analyses.
taining center hole.
Case No. of failed elements
4.2. Validation for failure stress envelope PFFa PMFb
B-1 0 0
4.2.1. Analysis scenario for initial failure B-2 0 10
In order to verify the proposed method for prediction of initial B-2–1 0 110
failure, the failure stress envelope for E-glass/LY556 and AS4/ B-2–2 0 200
B-7 0 32
3501-6 composite laminates that was already reported by Soden’s
B-8 0 74
research group [6–8] was introduced in this study. B-8–1 0 94
Table 9 lists the analysis scenarios for initial failure analyses B-8–2 0 176
based on three types of loading cases. In order to validate the pre- a
Puck fiber failure.
diction capability of the proposed computational method for com- b
Puck matrix failure.
posite laminate damage, failure and no failure cases were equally
selected. Fig. 7 represents the initial failure stress envelopes and
analysis cases for the A, B, and C analysis scenarios.
failures were separately predicted using the proposed ABAQUS
4.2.2. Analysis results and discussions for initial failure UMAT.
Fig. 8 illustrates the initial failure analysis results for cases A-2, Although no fiber failure was observed, matrix failures were
B-2, and C-2 based on the Puck failure criterion. In this figure, the detected under the 100 MPa (case A-2) and 300 MPa (case C-2)
gray domain means that the failure index exceeds 1.0 at each biaxial tension stress states, as shown in Fig. 8(b) and (f),
element. As shown in this figure, the fiber- and matrix-dependent respectively. Similarly, a matrix failure occurred at the edge of
416 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419

(a) Case B-1/B-2 (b) Case B-7/B-8

Fig. 12. Relationship between equivalent (von Mises) stress and number of damaged elements for case B-1/B-2 and Case B-7/B-8.

(a) Fiber damage at P = (b) Fiber damage at P = (c) Fiber damage at P =

17,000N 26,000N 32,000N

(d) Matrix damage at P = (e) Matrix damage at P = (f) Matrix damage at P =

17,000N 26,000N 32,000N

Fig. 13. Damage propagation of [(0/90)6]S laminate under each loading increment.

the boundary condition of the FE model, but no fiber failure is design tool for the detection of the initial failure of a composite
shown in Fig. 8(d). laminate because the analysis results are conservative from the
In particular, a matrix failure index cannot be estimated using viewpoint of the safety factor.
the classical failure criteria such as the maximum stress/strain,
Tsai–Hill, and Tsai–Wu criteria. Therefore, it is essential to 4.2.3. Analysis scenario for progressive failure
adopt the advanced failure criterion, e.g., the Puck criterion, for a Similar to the initial failure prediction analysis, the failure stress
precise estimation of the composite laminate-based structural envelope from initial to final failures that was already reported by
failure. Soden’s research group [6–8] was also adopted to validate the
Table 10 lists the two kinds of initial failure analysis results computational method in the present study. In particular, a pro-
using the proposed ABAQUS UMAT: the Puck fiber failure index gressive failure analysis was performed using case B (E-glass/
(PFFI) and Puck matrix failure index (PMFI). In addition, the agree- LY556 under in-plane normal and shear stresses) because a non-
ment and disagreement between the analysis and experiment are uniform distribution of the in-plane stress was observed in case
listed. B compared to the other analysis cases.
As shown in this table, the analysis results agree with the exper- Table 11 lists the analysis scenarios for the progressive failure
imental results with an accuracy of 95.8%. Although there is one analyses, and Fig. 9 illustrates the initial and final failure stress
overestimation case such as case B-3 in PMFI under the rx and envelopes, analysis cases, and failure states under the various rx
sxy stress state, the proposed analysis method might be a robust and sxy stress conditions.
 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419 417

(a) Fiber damage at P = (b) Fiber damage at P = (c) Fiber damage at P =

13,000N 17,000N 21,000N

(d) Matrix damage at P = (e) Matrix damage at P = (f) Matrix damage at P =

13,000N 17,000N 21,000N

Fig. 14. Damage propagation of [(0/±45/90)3]S laminate under each loading increment.

(a) Fiber damage at P = (b) Fiber damage at P = (c) Fiber damage at P =

9,800N 10,000N 10,200N

(d) Matrix damage at P = (e) Matrix damage at P = (f) Matrix damage at P =

9,800N 10,000N 10,200N

Fig. 15. Damage propagation of [(±45)6]S laminate under each loading increment.

Table 13
Comparison of failure stresses of T300/1304-C laminates.

Stacking sequence Chang and Chang Present study Error (%)

Initial failure stress (MPa) Final failure stress (MPa) Initial failure stress (MPa) Final failure stress (MPa) Initial failure stress Final failure stress
[(0/90)6]S 232 414 322 397 9.01 4.08
[(0/±45/90)3]S 177 281 161 260 8.67 7.43
[(±45)6]S 121 127 109 114 10.08 10.08
Average 9.25 7.20

In this figure, cases B-1 and B-7 show the no failure state, cases 4.2.4. Analysis results and discussions for progressive failure
B-2 and B-8 illustrate the initial failure state, cases B-2–1 and B-8– Figs. 10 and 11 represent the numerical analysis results for pro-
1 demonstrate the progressive failure state, and cases B-2–2 and B- gressive fiber and matrix failures based on the Puck criterion for
8–2 show the final failure state of the composite laminates. Hence, cases B-1, B-2, B-2–1, and B-2–2, and cases B-7, B-8, B-8–1, and
it was confirmed that the safe, failure onset, failed–not failed B-8–2, respectively. As shown in these figures, as the magnitude
mixed, and totally failed regions can be observed in the FE model of the applied stress increases, the gray domain that represents
during the computational analysis. the failed region increases.
418 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419

It is interesting that all of the elements in the matrix failure Fig. 13 shows the growth of fiber and matrix damage of [(0/
mode change to gray, but not all of the elements in the fiber failure 90)6]S laminate at the 17,000 N, 26,000 N and 32,000 N loading
mode do. In particular, none of the elements in the Puck fiber states. The initial matrix damage was observed in the region of
failure mode become gray. In other words, no progressive fiber fail- center hole under approximately 17,000 N. As the load was
ure could be observed in the Puck fiber failure during the calcula- increased, the matrix damage was propagated through the lami-
tion because the material damage of the fiber was not affected by nates. In addition, the fiber damage was also investigated at the
the aforementioned stress state. domain of center hole under about 26,000 N. The failed region
On the other hand, the characteristics of the progressive failure induced by matrix crack is much wider than that induced by fiber
rate were investigated in this study. Table 12 lists the failed ele- breakage. Finally, the test specimen was totally failed when the
ments from the no-failure to final-failure analyses in cases B-2 load was approached approximately 32,000 N.
and B-8, and Fig. 12 illustrates the relationship between the equiv- Fig. 14 demonstrates the growth of fiber and matrix damage of
alent (von Mises) stress and the number of damaged elements. The [(0/±45/90)3]S laminate at the 13,000 N, 17,000 N and 21,000 N
total number of elements in case B was 200. loading states. Similar to the [(0/90)6]S laminate, the matrix crack-
As shown in this table and figure, as the equivalent stress ing was initiated on the center hole region under approximately
increases, the growth rate of the Puck matrix failure rapidly 13,000 N. As the load is increased from 17,000 N to 21,000 N, the
increases. However, no Puck fiber failure was observed during cracking was propagated to the perpendicular direction from the
the analysis. Hence, it could be confirmed that the in-plane com- loading direction, especially, it was concentrated in the region of
plex stress could lead to a matrix failure rather than a fiber failure center hole. However, the mixed failure of fiber breakage and
of a composite laminate. matrix cracking was represented to the 45° from the loading direc-
Furthermore, the slope of the damage accumulation under a tion when the load was reached approximately 21,000 N.
large amount of compressive stress in the fiber direction (case B Fig. 15 represents the growth of fiber and matrix damage of
-7/B-8) is larger than that with a small amount of compressive [(±45)6]S laminate at the 8,800 N, 9,000 N and 9,200 N loading
stress (case B-1/B-2). Therefore, it could be theorized that a matrix states. This laminate show slightly different failure aspect among
failure expeditiously propagates under a compressive stress in the other laminates. Namely, the fiber breakage was occurred on the
fiber direction. This result means that if the combined stresses such center hole under about 9,800 N. And then, as the load is margin-
as the in-plane shear stress and compressive stress in the fiber ally increased, the small amount of fiber damage was propagated
direction are applied to a composite structure, a catastrophic fail- along the 90° from the loading direction. At last, the test specimen
ure can occur. was entirely fractured under approximately 10,200 N. One of the
Although the progressive failure rate could not be specifically distinguished characteristic of this test specimen is that the failed
validated because of the limitation of the experimental results, region of matrix cracking as well as fiber breakage is relatively
the progressive failure magnitude could be verified using the intro- sharp.
duced computational analysis method and ABAQUS UMAT, and the Although the damaged zone of each test specimen is not exactly
analysis results were in good agreement with the results reported coincided with the literature [11] due to the difference of failure
in the literature [6,7]. Hence, it was confirmed that the developed criterion, the damage propagation tendency of matrix cracking
method could be a good analysis tool for initial and progressive and fiber breakage under each load increment is corresponded
failure evaluations of composite laminate-based structures. with the literature as shown in Figs. 13–15.
On the other hand, Table 13 represents the comparison of stres-
4.3. Validation for damage propagation ses at initial and final failure states between literature and present
study. As shown in this table, the simulation results underesti-
4.3.1. Analysis scenario mated the ultimate stress approximately 7–10% compared to the
On the other hand, in order to validate the developed analysis calculation results of literature. Although there is somewhat of
method for simulation of damage propagation, the fiber and matrix an error in present analysis results, it was confirmed that the pro-
damage growth data of center-holed T300/1304-C test specimen posed analysis method could be considered as a practical design
during the uniaxial tensile test that was previously addressed by tool for the detection of the progressive failure of a composite lam-
Chang’s research team [11,30] was applied in this study. inate since the amount of error is reasonable and the analysis
In the experimental results carried by Chang’s research group, results are conservative from the viewpoint of the safety factor.
the ultimate tensile strengths as well as the damage propagations
were investigated with respect to seven kinds of stacking
sequences, namely, [(0/90)6]S, [(0/±45/90)3]S, [(902/±60/±30)2]S, 5. Concluding remarks
[(±45)6]S. [0/±45/907]S, [0/(±45)2/907]S and [0/(±45)3/903]S, and
three kinds of width ratio (W/D), namely, 4.0, 6.0 and 8.0 for In the present study, the onset and growth of cracks in a CFRP/
T300/1304-C composite laminates. However, in the present study, GFRP composite laminate were successfully estimated and
load–displacement and damage growth are simulated with respect predicted using the Puck failure criterion and a damage mecha-
to three kinds of stacking sequence and W/D = 4.0 width ratio, as nism-based computational method. Furthermore, the numerical
aforementioned in Table 7. Furthermore, the right edge is displace- analysis results were compared to the failure test results
ment-controlled to the fiber direction (X-direction) during the previously reported in the literature.
analysis. The primary research results and further study topics are listed
below.
4.3.2. Analysis results and discussions
Figs. 13–15 demonstrates the computational simulation results  A new computational analysis method for the practical evalua-
for fiber and matrix damage growth regarding to three kinds of tion of the initial and progressive failures of glass/carbon fiber-
T300/1304-C test specimens, namely, [(0/90)6]S, [(0/±45/90)3]S reinforced composite laminates using advanced assessment
and [(±45)6]S laminates based on Puck failure criterion. The specific function such as the Puck failure criterion was developed based
failure characteristics of each simulation results are discussed on the user-defined subroutine UMAT in the commercial FEA
below. code ABAQUS.
 C.-S. Lee et al. / Composite Structures 121 (2015) 406–419 419

 Through a comparative study, i.e., a comparison to the initial/ [8] Hinton MJ, Kaddour AS, Soden PD. A comparison of the predictive capabilities
of current failure theories for composite laminates, judged against
final failure stress envelope as well as the damage propagation
experimental evidence. Compos Sci Technol 2002;62:1725–97.
of various types of test specimens, it was confirmed that the [9] Kaddour AS, Hinton MJ, Soden PD. A comparison of the predictive capabilities
simulation results agreed well with the reported composite fail- of current failure theories for composite laminates: additional contributions.
ure results. Compos Sci Technol 2004;64:449–76.
[10] Germanischer Lloyd, Guideline for the Certification of Wind Turbines,
 Based on the proposed subroutine and methodology, it might be Germanischer Lloyd, 2010.
possible to conduct the structural failure analyses of large and [11] Chang FK, Chang KY. A progressive damage model for laminated composites
complex composite structures. Most significantly, a matrix containing stress concentrations. J Compos Mater 1987;21:834–55.
[12] Tay TE, Liu G, Tan VBC, Sun XS, Pham DC. Progressive failure analysis of
failure, which is the most difficult failure mode to predict, composites. J Compos Mater 2008;42:1921–46.
could be effectively estimated using the developed analysis [13] Liu PF, Zheng JY. Progressive failure analysis of carbon fiber/epoxy composite
tool. laminates using continuum damage mechanics. Mater Sci Eng A 2008;485:
711–7.
 The proposed assessment method could be a robust analysis/ [14] Reddy YSN, Dakshina Moorthy CM, Reddy JN. Non-linear progressive failure
design tool for the structural safety evaluation of composite analysis of laminated composite plates. Int J Non Linear Mech 1995;30:
laminate-based structures such as wind turbine blades and 629–49.
[15] Maimi P, Camanho PP, Mayugo JA, Davila CG. A continuum damage model for
marine composite propellers under static loads because the composite laminates: Part I – constitutive model. Mech Mater 2007;39:
method satisfies the evaluation methodology, which is regular- 897–908.
ized by classification society rules [10]. [16] Maimi P, Camanho PP, Mayugo JA, Davila CG. A continuum damage model for
composite laminates: Part II – computational implementation and validation.
 In a future study, the failure characteristics of composite
Mech Mater 2007;39:909–19.
laminates under fatigue and impact loads will be investigated, [17] Maimi P, Mayugo JA, Camanho PP. A three-dimensional damage model for
and an advanced estimation method will be proposed in order transversely isotropic composite laminates. J Compos Mater 2008;42:
to predict the structural performance of a composite laminate. 2717–45.
[18] Zhang BM, Zhao L. Progressive damage and failure modeling in fiber-
reinforced laminated composites containing a hole. Int J Damage Mech
2012;21:893–911.
Acknowledgements [19] Ubaid J, Kashfuddoja M, Ramji M. Strength prediction and progressive failure
analysis of carbon fiber reinforced polymer laminate with multiple interacting
holes involving three dimensional finite element analysis and digital image
This research was financially supported by the Ministry of correlation. Int J Damage Mech 2014;23:609–35.
Education (MOE) and National Research Foundation of Korea [20] Puck A, Kopp J, Knops M. Guidelines for the determination of the parameters in
Puck’s action plane strength criterion. Compos Sci Technol 2002;62:371–8.
(NRF-2013H1B8A2078346) through the Human Resource Training [21] Toi Y, Lee JM. Thermal elasto-viscoplastic damage behavior of structural
Project for Regional Innovation (2013 Pilot Project). In addition, members in hot-dip galvanization. Int J Damage Mech 2002;11:171–85.
this research was supported by Basic Science Research Program [22] Lee CS, Kim MH, Mahendran M, Lee JM. Computational study on the fatigue
behavior of welded structures. Int J Damage Mech 2011;20:423–63.
through the National Research Foundation of Korea (NRF) funded [23] Lee CS, Chun MS, Kim MH, Lee JM. Numerical evaluation for debonding failure
by the Ministry of Education (No. 2014R1A1A2006102). phenomenon of adhesively bonded joints at cryogenic temperatures. Compos
Sci Technol 2011;71:1921–9.
[24] Kim JH, Lee CS, Kim MH, Lee JM. Prestrain-dependent viscoplastic damage
References
model for austenitic stainless steel and implementation to ABAQUS user-
defined material subroutine. Comput Mater Sci 2013;67:273–81.
[1] Pavlou D. Composite materials in piping applications: design, analysis and [25] Kim SK, Lee CS, Kim JH, Kim MH, Lee JM. Computational evaluation of
optimization of subsea and onshore pipelines from FRP materials. DEStech resistance of fracture capacity for SUS304L of liquefied natural gas insulation
Publications; 2013. system under cryogenic temperatures using ABAQUS user-defined material
[2] Young YL. Fluid–structure interaction analysis of flexible composite marine subroutine. Mater Des 2013;50:522–32.
propeller. J Fluids Struct 2008;24:799–818. [26] Lee CS, Yoo BM, Kim MH, Lee JM. Viscoplastic damage model for austenitic
[3] Motley MR, Liu Z, Young YL. Utilizing fluid–structure interactions to improve stainless steel and its application to the crack propagation problem at
energy efficiency of composite marine propellers in spatially varying wake. cryogenic temperatures. Int J Damage Mech 2013;22:95–115.
Compos Struct 2009;90:304–13. [27] Kim JH, Kim SK, Kim MH, Lee JM. Numerical model to predict deformation of
[4] Hashin Z. Failure criteria for unidirectional fiber composites. J Appl Mech corrugated austenitic stainless steel sheet under cryogenic temperatures for
1980;47:329–34. design of liquefied natural gas insulation system. Mater Des 2014;57:26–39.
[5] Puck A, Schürmann H. Failure analysis of FRP laminates by means of [28] Lee CS, Lee JM. Failure analysis of reinforced polyurethane foam-based LNG
physically based phenomenological models. Compos Sci Technol 1998;58: insulation structure using damage-coupled finite element analysis. Compos
1045–67. Struct 2014;107:231–45.
[6] Soden PD, Hinton MJ, Kaddour AS. Lamina properties, lay-up configurations [29] McCarthy CT, McCarthy MA, Lawlor VP. Progressive damage analysis of multi-
and loading conditions for a range of fibre-reinforced composite laminates. bolt composite joints with variable bolt-hole clearances. Compos B Eng
Compos Sci Technol 1998;58:1011–22. 2005;36:290–305.
[7] Soden PD, Hinton MJ, Kaddour AS. A comparison of the predictive capabilities [30] Chang FK, Scott RA, Springer GS. Failure of composite laminates containing pin
of current failure theories for composite laminates. Compos Sci Technol loaded holes – method of solution. J Compos Mater 1984;18:255–78.
1998;58:1225–54.