Article

Analytical model for determining
effective stiffness and mechanicalbehavior of polymer matrix
composite laminates using
continuum damage mechanics
Sota Onodera1

International Journal of Damage
Mechanics
2020, Vol. 29(10) 1512–1542! The Author(s) 2020
Article reuse guidelines:
sagepub.com/journals-permissions
DOI: 10.1177/1056789520939624
journals.sagepub.com/home/ijd

and Tomonaga Okabe1,2Abstract
The present paper proposes a new analytical model for predicting the effective stiffness of composite
laminates with fiber breaks and transverse cracks. The model is based on continuum damage mechanics
and the classical laminate theory. We derived damage variables describing stiffness reduction due to fiber
breaks and its maximum value during ultimate tensile failure from the global load-sharing model.
Furthermore, a simplified analytical model is presented for obtaining two damage variables for a cracked
ply subjected to transverse tensile loading or in-plane shear loading. This model was developed assumingthat the displacement field of the longitudinal direction can be expressed in the form of a quadric function
by loosening the boundary condition for the governing differential equation. For verifying the developedmodel, the elastic constants of damaged composite laminates were predicted for cross-ply and angle-ply
laminates and compared with the finite element analysis results. As for the appropriate expression of theeffective elastic stiffness matrix of the damaged ply, we verified four types of effective compliance/stiffness
matrices including the Murakami, Yoshimura, Li, and Maim!ı models. We found the Maim!ı model to be themost appropriate among these four models. Moreover, we successfully simplified the expressions fordamage variables in the complicated infinite series obtained in our previous study. We also proved that
this could contribute toward improving the accuracy of our analysis. After verifying the present model,
the stress–strain response and failure strength of carbon- or glass-fiber-reinforced plastic cross-ply
laminates were predicted using Maim!ı’s compliance model and the simplified damage variables.
Keywords
Transverse crack, continuum damage mechanics, damage variable, composite laminate, polymer matrix
composite

1
2

Department of Aerospace Engineering, Tohoku University, Miyagi, Japan
Department of Materials Science and Engineering, University of Washington, Seattle, WA, USA

Corresponding author:
Sota Onodera, Department of Aerospace Engineering, Tohoku University, 6-6-01, Aramaki Aza Aoba, Aoba-ku, Sendai, Miyagi 9808579, Japan.
Email: sonodera@plum.mech.tohoku.ac.jp

Onodera and Okabe

1513

Introduction
The application scope of advanced composite materials is increasing due to their superior mechanical properties like high specific modulus and strength. The specific modulus plays an important role
in real structures, but the specific strength of these materials has not been effectively utilized in realapplications owing to the incapability of failure analysis. The reason why effective designs considering the high specific strength are not yet developed is that structural designs are developed using
the strain criterion based on initial ply failure. When load is applied to composite structures, matrix
cracks are generated in the transverse ply. However, if the applied strain level is low, these cracks do
not cause secondary damage and affect the reduction of stiffness. Furthermore, if the applied strainreaches a certain level, the damage migrates from the ply cracks, leading to catastrophic damage
phenomena like delamination and fiber breakage. Thus, similar to metals, if designers consider the
initial damage and users can predict the progress of the damage, the structural weight can be
reduced. Therefore, damage tolerance designing is the next step for expanding the application
scope of composite structures.
Designers need to analyze the stiffness reduction in composite structures owing to transverse
cracks accurately. Transverse cracks occurring in cross-ply laminates have been studied since the1970s (Bailey and Parvizi, 1981; Garrett and Bailey, 1977a, 1977b; Manders et al., 1983; Parvizi andBailey, 1978; Parvizi et al., 1978). Initially, there were many stress analysis studies performed using
the shear-lag model, which was developed to analyze the stress profile around discontinuous fibers(Cox, 1952). Using the shear-lag model, the damage progress could be predicted using the stresscriterion and/or energy criterion. Specifically, the onset of transverse cracks from the viewpoint of
micromechanics was analyzed. In the 1980s, the shear-lag model was improved (Han et al., 1988;
Lee and Daniel, 1990). Finite element analysis (FEA) (Wang et al., 1984) and the minimum energy
principle model (Hashin, 1986, 1987; Lee et al., 1989) utilizing variation methods were proposed toanalyze the stress field around the transverse cracks. These models enable us to calculate the stressfield around the transverse cracks in cross-ply laminates very accurately.
In addition to the above approaches, continuous damage mechanics (CDM) models (Allen et al.,1987a, 1987b; Talreja, 1985) were developed simultaneously. CDM models are divided into the following two types: the first one involves the effective elastic stiffness or compliance matrix utilizing thedamage tensor derived from the thermomechanical relationship, while the other one calculates theinelastic strain with internal state variables so that the average stress of the laminate can be estimated.
Both models are fitted with the classical laminate theory (CLT) because the detailed stress distribution
in the laminate is unnecessary, and the damaged ply is homogenized. Therefore, arbitrary laminates
other than cross-ply laminates can be analyzed easily. On the contrary, problems regarding theobjective determination of damage tensor and internal state variables must be addressed.
In 1992, McCartney presented an accurate anisotropic elastic solution for cross-ply laminates,
satisfying both the equilibrium and compatible equation. Pagano et al. (1998) proposed an accurate
solution based on cylindrical bending. Moreover, several researchers analyzed the delamination
migrating from the matrix crack from the view point of energy balance (Nairn and Hu, 1992;
Noda et al., 2004). Xia et al. (1993) used the McCartney model to analyze damage progress inceramic–matrix cross-ply laminates. Furthermore, using CDM, Gudmundson and Zang (1993)
extended the two-dimensional model of Dvorak et al. (1985, 1987) which is applicable to an isotropic
body to obtain a three-dimensional model with the generalized laminate theory for analyzing threedimensional laminates. As stated above, CDM can be applied to angle-ply laminates. Gudmundson
and Zang (1993) applied the model to study stiffness reduction in angle-ply laminates. Kobayashi
et al. (2000) used the Gudmundson and Zang (1993) model to predict the damage progress in1514International Journal of Damage Mechanics 29(10)quasi-isotropic carbon-fiber-reinforced plastic (CFRP) laminates. After the 2000s, several researchers(Hajikazemi and McCartney, 2018; McCartney, 2017; Vinogradov and Hashin, 2010) continued todevelop models for analyzing arbitrary laminates other than cross-ply laminates. However, thederived formulations are extremely complicated for ordinary people to use them easily.
Recently, we derived the analytical form of the damage tensor with transverse crack density bysolving the Laplace equation for the displacement fields under the generalized plane strain condition(Okabe et al., 2018). The laminate stiffness can be estimated with the obtained damage tensor
through the CLT and CDM model proposed by Murakami and his colleagues (2012) and
Kanagawa et al. (1996). We have also proved that the estimated stiffness agrees with the finite
element calculations and experimental results (Talreja, 1985). This approach is extended to a threedimensional model (Onodera and Okabe, 2019) using a sophisticated model proposed by Li et al.
(2019). However, the effective compliance matrix proposed by Murakami and his colleagues (2012)is asymmetric. This asymmetry of the compliance matrix was also observed in the study done by Lee
et al. (1989). Furthermore, because the damage tensor is derived by solving the Laplace equation, it
is written as a complicated infinite series (Okabe et al., 2018).
The fiber breakage generally causes nonlinearity in the stress–strain curves when the applied
strain is large. It is well known that Rosen (1994) originally studied the stress transfer length aroundfiber breaks as well as the unidirectional composite failure strength and applied to glass-epoxy
monolayer specimens. Analytical models considering fiber breakage have been mainly developedin the field of ceramic matrix composites (CMCs) since the late-1980s. This is because a stressconcentration due to a fiber break can almost be neglected, so that the global load sharing
model (GLS), wherein the load shed by a broken fiber is shared by all intact fibers, can be directly
applied to the experiments. Curtin (1991) presents a famous theory to predict the ultimate tensile
strength of the unidirectional composites assuming the GLS model. Neumeister (1993) developed an
approximate solution of the stress–strain curve and strength of the composite and compared to anumerical result based on the fiber fragmentation theory formulated by Curtin (1991). Hui et al.
(1995) derived an explicit closed-form solution of the evolution of fiber fragments and applied toestimate the strength of unidirectional composite materials. In general, the GLS models give a larger
strength than the local load-sharing (LLS) model that considers local stress concentrations
caused by fiber breaks (Ibnabdeljalil and Curtin, 1997). However, the GLS models are theoreticallyeasy to handle. Therefore, it is suitable to utilize the GLS models to analytically formulate thedamage tensor.
This study developed an appropriate expression for the effective stiffness matrix of damaged
composite laminates and derived simplified expressions for the damage tensor with the fiber breakage and transverse crack density. Damage variable for fiber breakage was formulated using the GLS
model for unidirectional fiber-reinforced composites. An approximate simplified expression for thedamage tensor with transverse crack density was obtained from the Laplace equation, assuming thequadratic function. For validating the present model, the obtained damage tensor was substituted
into CDM and CLT models, and the elastic constants of damaged composite laminates were thenpredicted for cross-ply and angle-ply laminates. As for the appropriate expression for the effective
compliance/stiffness matrix of a damaged ply, we verified four types of damage matrices including
the Murakami (2012) (Kanagawa et al., 1996), Yoshimura (2011), Li et al. (2019) and Maim!ı models(Maim!ı et al., 2007). We determined that the Maim!ı model was the most appropriate among thesefour models. Moreover, we succeeded in simplifying the expression for damage variables withcomplicated infinite series obtained in our previous study. We also proved that this could contributetoward the improving the analysis accuracy. Further, we determined the nonlinear stress–strain
response and ultimate tensile strength of CFRP and glass-fiber-reinforced-plastic (GFRP) cross-ply

Onodera and Okabe

1515

laminates. We employed the Maim!ı model and simplified damage variables proposed in this study
when calculating the stress–strain response. The ultimate tensile strength of CFRP cross-ply laminates was predicted based on a strain-based fiber failure criterion (McCartney, 2005) because fiber
fracture in CFRP laminates often occurs at the point of laminate failure. For GFRP laminates, thestiffness reduction due to multiple-fiber fragmentation in 0! plies and the subsequent final laminate
failure were determined using the GLS model.

Model formulation
Effective elastic compliance/stiffness matrix of damaged ply
Four models for the elastic compliance matrix of a damaged composite ply under the plane stresscondition are introduced here, and their characteristics are explained.

Model 1: Murakami model
Murakami published a sophisticated book (Murakami, 2012) introducing CDM and proposed acompliance matrix for damaged composite plies with his colleagues (Kanagawa et al., 1996). They
used the strain equivalent principle to obtain the following stress–strain relationship
0

1
0
1 B
B ð1 $ d11 ÞE1
e11
B
! 12
B
C B
@ e22 A ¼ B
B $ ð1 $ d11 ÞE1
B
2e12
B
@
0

! 12
$
ð1 $ d22 ÞE1
1
ð1 $ d22 ÞE2
0

1

0
0
!

"

1
1
1
1
þ
2 1 $ d11 1 $ d22 G12

C0
1
C
C r11
CB
C@ r22 C
A
C
C
C r12
A

(1)

Ei is the Young’s modulus, Gij is the shear modulus, ! ij is the Poisson’s ratio, eij is the strain tensor,
rij is the stress tensor, and dij is the damage variable. In equation (1), subscripts 1 and 2 denote thefiber direction and the in-plane transverse direction, respectively. In our previous study, we foundthat this matrix is effective for obtaining the elastic constants of damaged composite laminates if theappropriate damage tensor is used. However, the compliance matrix in equation (1) is obviously
asymmetric, and this violates the existence of the strain energy function. Moreover, this form cannot
be directly entered into the numerical simulation code as in FEA. As stated above, this asymmetryof the compliance matrix was also observed in the work of Lee et al. (1989).

Model 2: Yoshimura model
Yoshimura et al. (2011) proposed the concept of the stress equivalent principle to overcome thedifficulty of the asymmetry of the compliance matrix. The stress–strain relation can be written as
0

r11

1

0

B11

B12

B
C B
B22
@ r22 A ¼ @
r12
sym:

0

10

e11

1

CB
C
0 A@ e22 A
2e12
B33

(2)

1516

International Journal of Damage Mechanics 29(10)

where
B11 ¼

B12 ¼

#

ð1 $ d11 ÞE1 d211 ! 12 ! 21 ð1 þ ! 23 ÞE1$
1 $ ! 12 ! 21
4 Dð1 $ ! 12 ! 21 Þ$
22
! 12 E2
1 $ d11 þd2
1 $ ! 12 ! 21

B22 ¼

$

d11 d22 ! 12 ð! 23 þ ! 12 ! 21 ÞE2
4
Dð1 $ ! 12 ! 21 Þ

ð1 $ d22 ÞE2 d222 ð! 23 þ ! 12 ! 21 Þ2 E2$
1 $ ! 12 ! 21
4 Dð1 þ ! 23 Þð1 $ ! 12 ! 21 Þ
B33 ¼

!

"
d11 þ d22
G12
1$
2

(3)

(4)

(5)

(6)

D ¼ 1 $ ! 23 $ 2! 12 ! 21(7)

In this case, the stiffness matrix is obviously symmetric. This effective stiffness matrix has alreadybeen inserted into the ABAQUS user subroutine to simulate the high-velocity impact problem. The
simulation results roughly reproduced the damage progress observed experimentally, but the simulation underestimated the extent of damage compared to the experiments. Moreover, the degradation coefficient of shear modulus does not have any physical meaning. These points are discussed
later when comparing the four models explained in this section.

Model 3: Li model
Li et al. (2019) formulated the three-dimensional effective stiffness matrix for a ply with transverse
cracking assuming a small damage. The two-dimensional effective stiffness matrix can be written as
0

C213
C
$
11
B
r11
C33
B
B
C B
@ r22 A ¼ B
B
@
r12
sym:
0

1

C13 C23
C12 $
C33
C2
C22 $ 23
C33

1

0 C0 e 1
11
CCCB
C@ e22 A
C
0 A
2e12
C66

(8)

where
1 $ ! 23
E1 ! 12 ! 21
E1 $ x
D
D2

(9)

! 12
E2 ! 12 ð1 $ ! 12 ! 21 Þ
E2 $ xD
ð1 þ ! 23 ÞD2

(10)

C11 ¼
C12 ¼

Onodera and Okabe

1517

! 12
E2 ! 12 ð! 23 þ ! 12 ! 21 Þ
E2 $ xD
ð1 þ ! 23 ÞD2(11)

1 $ ! 12 ! 21
E2 ð1 $ ! 12 ! 21 Þ2
E2 $ xð1 þ ! 23 ÞD
ð1 þ ! 23 Þ2 D2(12)! 23 þ ! 12 ! 21
E2 ð1 $ ! 12 ! 21 Þð! 23 þ ! 12 ! 21 Þ
E2 $ xð1 þ ! 23 ÞD
ð1 þ ! 23 Þ2 D2(13)

1 $ ! 12 ! 21
E2 ð! 23 þ ! 12 ! 21 Þ2
E2 $ xð1 þ ! 23 ÞD
ð1 þ ! 23 Þ2 D2

(14)

C13 ¼

C22 ¼
C23 ¼

C33 ¼

C66 ¼ G12 $ nxG12(15)Here, x is the damage parameter related to tensile damage, and n is the constant damage parameter related to shear damage. x and n can be expressed as follows (Onodera and Okabe, 2019)
x¼1$

Ed2
¼ d22
E2

(16)

!
"
1
Gd12
1
1$
'
n¼
x
k
G12

(17)

Here, superscript d denotes the property of the damaged ply and k is the material constantdescribed later. Onodera and Okabe (2019) utilized this effective stiffness for predicting the stiffness
reduction and steady-state cracking of laminates with ply cracking. The proposed CDM model
agrees with the experimental and numerical results (Onodera and Okabe, 2019). However, thedamage parameter x and transverse crack density have upper limits because a small damage wasassumed.

Model 4: Maim!ı model
Maim!ı et al. (2007) proposed the complementary free energy function and thermodynamicallyderived the effective compliance matrix as follows
0

1
0
1 B
B ð1 $ d11 ÞE1
e11
B
B
C B
@ e22 A ¼ B
B
B
2e12
B
@
sym:

! 12
$
E1
1
ð1 $ d22 ÞE2001
ð1 $ d12 ÞG12

1

C0
1
C
C r11
CB
C@ r22 C
A
C
C
C r12
A

(18)

1518

International Journal of Damage Mechanics 29(10)The effective compliance matrix in this case is also explicitly symmetric. However, this case needs
an additional damage variable d12 for evaluating the shear stiffness reduction due to transverse
cracks or fiber breaks. For the ply with transverse cracks, this parameter is derived from the displacement field of the damaged ply in this study. For the ply with fiber breaks, d12 can be approximated as follows (Murakami, 2012)
!
"
1
1
1
¼
þ1
1 $ d12 2 1 $ d11(19)

The approximate relation of equation (19) was used when considering the ply with fiber breaks.

Damage variables
Damage variable d11 for multi-fiber fragmentation. An analytical expression for d11 was developed by
employing Curtin’s GLS model (1991) considering the ultimate tensile strength of unidirectional
fiber-reinforced composites. Now, we consider the mechanical equilibrium on an arbitrary crosssection in a unidirectional fiber-reinforced composite with broken fibers. Utilizing the shear-lag
model, the slip length lf can be described as
lf ¼

rT
2s

(20)

where T is the far-field axial fiber stress, s is the fiber/matrix interfacial sliding stress, and r
is the fiber radius. If fiber breaks occur equally within the range (lf , then the average
stress carried by the broken fibers in the selected plane is half of that carried by the
unbroken fibers. The stress carried by the broken, and unbroken fibers in the selected plane canbe described as%& T" f ¼ T 1 $ Pf ðT; 2lf Þ þ Pf ðT; 2lf Þ
r
2

(21)

Here, the first term on the right-hand side represents the load carried by the unbroken fiber withthe fraction 1 $ Pf ðT; 2lf Þ, and the second term is the average load carried by the broken fibers,
whose fraction is Pf ðT; 2lf Þ. It is well known that the probability Pf ðT; 2lf Þ of the fiber failure withapplied stress T and length 2lf can be simplified as! "! "mþ1
2lf T mT
Pf ðT; 2lf Þ '
¼
L0 r0
rc

if

! "
2lf T m
)1
L0 r0

(22)

where m is the Weibull modulus and r0 is the mean strength at gauge length L0. rc is the characteristic strength defined as follows
rc ¼!

r0 sL0r"1=ðmþ1Þ

(23)

Onodera and Okabe

1519

In the selected cross section, the following mechanical equilibrium condition must be satisfied
rf
r1 ¼ f"(24)

r1 is the stress applied on the unidirectional composite along the fiber direction and f is the fiber
volume fraction. Using the Voigt model for representing the Young’s modulus of composite, theunidirectional axial composite stress r1, axial fiber stress T, and unidirectional Yong’s modulus E1
are expressed as
r1 ¼ ð1 $ d11 ÞE1 e1(25)

T ¼ Ef e1(26)

E1 ¼ fEf þ ð1 $ fÞEm ' fEf(27)

Here, e1 is the unidirectional composite strain, Ef is the fiber axial modulus, and Em is the matrix
modulus. We neglected the term ð1 $ fÞEm in equation (27) because Ef * Em is approximately
satisfied in the case of polymer matrix composites. Using equations (21), (22), and (25) to (27) inequation (24), the damage variable was simply derived as
d11 ¼

! "
1 T mþ1
2 rc(28)

The remaining unknown variable is the far-field fiber axial stress T. Substituting equations (21)and (22) into equation (24), the following nonlinear equation for T was obtained! " !
1 T mþ1
r1 ¼ fT 1 $
2 rc

(29)

Stress T at applied load r1 was numerically obtained by solving the above nonlinear equation. When!
!
a cross-ply laminate having 0 plies with volume fraction V1 and the 90 plies with volume fraction V2 is" L can be described using the rule of mixture as follows
considered, the applied axial laminate stress r" L ¼ V1 r1 þ V2 r2 ' V1 r1
r
!

(30)

!

r1 is the 0 ply stress and r2 is the 90 ply stress. We decided to neglect the second term on the!
right-hand side in equation (30). This is because fiber breaks in the 0 ply may start after the!
saturation of transverse cracks, and the 90 plies have a very small load-bearing capacity during
saturation of transverse cracks. Using equation (30) in equation (29), the nonlinear expression for
stress T is rewritten as! " !
"L
1 T mþ1r¼ fT 1 $
V1
2 rc

(31)

1520

International Journal of Damage Mechanics 29(10)The ultimate tensile strength of cross-ply laminates can be derived using equation (31). The
ultimate failure condition of composite laminates is described by the following equation
d"
rL
¼0
dT

(32)

" ULT
The maximum far-field fiber stress Tmax and ultimate tensile strength rof cross-ply laminatesL
were simply obtained by substituting equations (31) into (32)Tmax ¼!

" ULT¼ V1 frcr
L"1=ðmþ1Þ
2
rc
mþ2

(33)

!

(34)

"1=ðmþ1Þ
2
mþ1
mþ2
mþ2

The maximum value of damage variable d11 was obtained using equation (33) in equation (28)
dmax
11 ¼

1
mþ2

(35)

d11 increased to its maximum value in equation (35) under the ultimate tensile failure condition. In
this study, the far-field stress T was calculated using the Newton method to solve the nonlinear" L until the ultimate tensile failure conequation (equation (31)) at a given laminate tensile stress r
dition in equation (32) was satisfied. d11 was obtained using the calculated stress T in equation (28).
Damage variable d22 associated with transverse cracking. In any of the effective stiffness/compliancemodels stated above, d22 is necessary to evaluate the compliance/stiffness matrix of the damaged
ply. The governing equation (i.e. Laplace equation) of the corresponding problem was derived. The
simplified solution of the equation obtained assuming that the displacement field in the longitudinal
direction can be expressed as a quadric function. Further, the boundary condition for the governing
equation was loosened.
As shown in Figure 1, the solved problem is the displacement field for the damaged transverse ply
including periodic matrix cracks having an interval of 2l. These matrix cracks fully propagate in thez-direction and reach the interface. The longitudinal strain of the neighboring ply is assumed to beuniform and unaffected by these cracks. Considering the symmetry of the problem, only the first
quadrant in the xy plane was addressed. The plane strain condition was assumed, and the constitutive relations are given byexx ðx; yÞ ¼@uðx; yÞ
¼ C1 rxx ðx; yÞ $ C2 ryy ðx; yÞ
@x

(36)

eyy ðx; yÞ ¼@vðx; yÞ
¼ C1 ryy ðx; yÞ $ C2 rxx ðx; yÞ
@y

(37)

rxy ðx; yÞ @vðx; yÞ @uðx; yÞ @vðx; yÞ
þ
'
¼
@x
@y
@x
G23

(38)

cxy ðx; yÞ ¼

Onodera and Okabe

1521

3
2
1

Transverse crack

x
y
t
t

z
l

l

Figure 1. Representative volume element in a cracked ply with transverse crack spacing 2l and ply thickness 2t.
(1, 2, 3) is the material coordinate system, and (x, y, z) is the representative volume element coordinate system.

Here, u and v are the displacements along the x and y directions, respectively. C1 and C2 aregiven by
C1 ¼

1 $ ! 12 ! 21
;
E2

C2 ¼

! 23 þ ! 12 ! 21
E2

(39)

We assumed that @u=@y was negligible in equation (38) so that the compatible condition of thestrain tensor was not satisfied. To obtain the governing equation for displacement v, the following
relationship was introducedexx ðx; yÞ ¼ aeyy ðx; yÞ(40)

where a is the effective Poisson’s ratio derived from the equilibrium equation. Substituting equation(40) into equations (36) and (37), we obtained the stress with displacement v as follows
rxx ðx; yÞ ¼

aC1 þ C2 @vðx; yÞ@yC21 $ C22(41)

ryy ðx; yÞ ¼

C1 þ aC2 @vðx; yÞ@yC21 $ C22

(42)

These stresses satisfy the following two-dimensional equilibrium equations@rxx @rxy
þ
¼0
@x
@y

(43)

@rxy @ryy
þ
¼0
@x
@y

(44)

1522

International Journal of Damage Mechanics 29(10)Substituting equations (38) and (41) into equation (43) and removing the trivial case (i.e.
@ 2 v=@x@y 6¼ 0), a was calculated as
C2 þ G23 ðC21 $ C22 Þ
C1
E2 ð! 23 þ ! 12 ! 21 Þ þ G23 fð1 $ ! 12 ! 21 Þ2 $ ð! 23 þ ! 12 ! 21 Þ2 g¼$
E2 ð1 $ ! 12 ! 21 Þ

a¼$

(45)

Furthermore, substituting equations (38) and (42) into equation (44), the following governing
equation of the displacement field was obtained
2
@2v
2@ v
þ
k
¼0
@x2
@y2

(46)

sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
C1 þ aC2
E2 $ G23 ð! 23 þ ! 12 ! 21 Þ
k¼¼
G23 ð1 $ ! 12 ! 21 Þ
G23 ðC21 $ C22 Þ

(47)

where

The boundary condition of the governing equation is written as
rxy ð0; yÞ ¼ 0 ()@v
ð0; yÞ ¼ 0
@x

(48)

vðt; yÞ ¼ e1
yy

(49)

vðx; 0Þ ¼ 0(50)

ryy ðx; lÞ ¼ 0 ()@v
ðx; lÞ ¼ 0
@y

(51)

where e1
yy is the applied strain of the ply parallel to the y-axis and 2t is the ply thickness. The
displacement field can be solved with equation (46). In our previous study (Okabe et al., 2018), thispartial differential equation was solved by the variable separation method. However, the obtained
solution was an infinite series and very complicated. Therefore, this study obtained a simplified formof the displacement field.
To obtain the simplified solution, the displacement form of the following type was assumed
vðx; yÞ ¼ AðyÞx2 þ BðyÞx þ CðyÞ(52)

Using the displacement in equation (52) in equations (48) and (49) yieldedBðyÞ ¼ 0(53)

Onodera and Okabe1523

2CðyÞ ¼ e1yy y $ AðyÞt

(54)

Thus, equation (52) can be written as
vðx; yÞ ¼ AðyÞðx2 $ t2 Þ þ e1
yy y

(55)

Next, our main task was to explicitly determine A(y). The boundary condition of equation (50)givesAð0Þ ¼ 0(56)

However, it is impossible to satisfy the boundary condition of equation (51) exactly as it is. Thus,
according to McCartney (1992), the following averaged boundary condition was used
1
t

Z t

@vðx; lÞdx ¼ 0@y0(57)Hence, using equation (55) into equation (57) results indAðlÞ 3e1
yy
¼ 2
2t
dy

(58)

Similar to equation (51), the government equation, i.e. equation (46), cannot be solved directly.
Instead of equation (46), the following averaged government equation was here solved
1
t

Z t
0

!
2
@2v
2@ v
þk
dx ¼ 0
@x2
@y2

(59)

Using equation (55) in equation (59) yielded the following ordinary differential equationd2 AðyÞ
$ b2 AðyÞ ¼ 0
dy2

(60)

where
b¼

pﬃﬃﬃ
3
kt

(61)

The solution of equation (60) with the boundary conditions in equations (56) and (58) isAðyÞ ¼

3e1
yy sinhby
2t2 b coshbl

(62)

1524

International Journal of Damage Mechanics 29(10)Thus, it can be shown that
vðx; yÞ ¼

3e1
yy sinhby
2t2 b coshbl

ðx2 $ t2 Þ þ e1
yy y

(63)

As stated in our previous paper (Okabe et al., 2018), the damage variable is the ratio of theinelastic strain due to crack opening to the total strain and is finally given by

d22 ¼

1
e1
yy $ ltZ tvðx; lÞdx

0
e1
yy

¼

tanhbl 2qb¼ tanh
bl
b
2q

(64)

where q ¼ 1=ð2lÞ is the transverse crack density in the damaged ply. As stated above, we derived thefollowing equation (Okabe et al., 2018)d22 ¼ 1 $1
8 X
1tanh½ð2n $ 1Þpktq,
p3 kt n¼1 ð2n $ 1Þ3
q(65)

Equation (64) is obviously simpler than equation (65), and loosening the boundary condition
equation (51) increases crack opening so that equation (64) gives better predictions than equation(65) in some cases. We will describe this section later through verifications and validations.
Damage variable d12 associated with transverse cracking. When the effective compliance matrix shown inequation (18) is used, the analytical form for damage variable d12 is required. The schematic shownin Figure 1 is also considered here. The neighboring ply deforms under a constant shear strain c1yz in
the yz plane, and the displacement field is assumed as
uðx; y; zÞ ¼ vðx; y; zÞ ¼ 0(66)

w ¼ wðx; yÞ(67)Thus, these equations give rx ¼ ry ¼ rz ¼ rxy ¼ 0 and ex ¼ ey ¼ ez ¼ cxy ¼ 0 so that the out ofplane shear stresses are nonzero and given asrxz ðx; yÞ ¼ G12 cxz ðx; yÞ ¼ G12@wðx; yÞ@x(68)

ryz ðx; yÞ ¼ G12 cyz ðx; yÞ ¼ G12@wðx; yÞ@y(69)

The stresses should satisfy the following three-dimensional equilibrium equations@rxx @rxy @rxz
þ
þ
¼0
@x
@y
@z

(70)

Onodera and Okabe

1525@rxy @ryy @ryz
þ
þ
¼0
@x
@y
@z

(71)

@rxz @ryz @rzz
þ
þ
¼0
@x
@y
@z

(72)

These stresses given in equations (68) and (69) automatically satisfy equations (70) and (71) andlead to the following governing equation derived from equation (72)
@2w @2w
þ¼0@x2 @y2(73)

As mentioned in the previous section, the following boundary conditions should be satisfied:rxz ð0; yÞ ¼ 0 ()@w
ð0; yÞ ¼ 0
@x

(74)

wðt; yÞ ¼ c1
yz y

(75)

wðx; 0Þ ¼ 0(76)

ryz ðx; lÞ ¼ 0 ()@w
ðx; lÞ ¼ 0
@y

(77)

This boundary problem defined by equations (73) to (77) is almost the same as the one defined byequations (46) to (51); therefore, the same derivation procedure can be adopted for this boundary
problem. Thus, the displacement field is given by
wðx; yÞ ¼

3c1
yz sinhb12 y
2t2 b12 coshb12 l

ðx2 $ t2 Þ þ c1
yz y

pﬃﬃﬃ
3
b12 ¼
t

(78)

(79)

The shear damage variable d12 can be defined as
d12 ¼ 1 $

Gd12
G12

(80)

where Gd12 is the shear stiffness of the ply with transverse cracking. The average ply shear strain caveyz
not considering the crack-opening displacement and applied shear strain c1
yz are described as follows
cave
yz ¼ 1 $

ryz
G12

(81)

1526

International Journal of Damage Mechanics 29(10)c1
yz ¼ 1 $

ryz
Gd12

(82)

Substituting equations (81) and (82) into equation (80), d12 can be rewritten as
d12 ¼ 1 $

cave
yz
c1
yz

(83)

The average shear strain is defined by
cave
yz ¼

1
lt

Z t

wðx; lÞdx

(84)

0

We obtained d12 from equations (78), (83), and (84)

d12 ¼

1
c1
yz $ ltZ twðx; lÞdx

0

c1
yz

¼

tanhb12 l 2qb¼
tanh 12
2q
b12 l
b12

(85)

The reduction in shear modulus was initially determined by Hashin (1986), by using the variation
approach. Recently, Onodera and Okabe (2019) also obtained the shear damage variable with an
infinite seriesd12 ¼ 1 $1
8 X
1tanh½ð2n $ 1Þptq,3
3
p t n¼1 ð2n $ 1Þq(86)

However, equation (85) is obviously simpler than the infinite series of equation (86).

Effective elastic constants for the damaged composite laminates
The substitution of damage variables d22 and d12 into equations (1), (2), (8), and (18) gives theeffective compliance/stiffness matrix of a ply with transverse cracks as a function of crack density q,
and fiber breaks can be considered by substituting d11 into equations (1), (2), and (18). The effective
stiffness of composite laminates with transverse cracks can be estimated through the CLT using theeffective compliance matrix of a damaged ply.
Let us consider the coordinate systems shown in Figure 2. The effective compliance matrix andthermal expansion coefficient of the kth ply having transverse crack density q in material coordinate
system (1, 2, 3) are defined as Ck ðqÞ and ak , respectively. Similarly, those of laminate coordinate" k ðqÞ and "ak , respectively. When the fiber angle is given by hk, C" k ðqÞ
system (X, Y, Z) are defined as Cand "ak are given as" k ðqÞ ¼ Rk Ck ðqÞT$1
C
k

(87)

"ak ¼ Rak

(88)

Onodera and Okabe

1527

1st
2nd
3rd

Y

2

1
θk

kth

Z

X

(N-1)thNthtL

kth plyY

X
Figure 2. Relationship between the laminate coordinate system (X, Y, Z) and material coordinate system (1, 2, 3).

where(
ak ¼ a112

6
Rk ¼ 64

a22

0cos2 hksin2 hksin2 hkcos2 hk2sinhk coshk$2sinhk coshk2

6
Tk ¼ 64

cos2 hksin2 hksin2 hksinhk coshk

cos2 hk$sinhk coshk

)T

(89)

3
$sinhk coshk
7
sinhk coshk 7
5
cos2 hk $ sin2 hk$2sinhk coshk

3

7
2sinhk coshk 7
5
cos2 hk $ sin2 hk

(90)

(91)

where a is the thermal expansion coefficient. Therefore, according to the CLT, the effective laminate
strain "e L is expressed as" Lr"e L ¼ C"L þ "aL DT(92)" L ðqÞ and thermal expansion coefficient "aL ðqÞ of the cracked
where the effective stiffness matrix C
laminate can be written as()" L ðqÞ ¼ C" Lij ¼
C("aL ðqÞ ¼ "a Li

)

"

N
1X
" $1 ðqÞ
tk C
ktL k¼1

"

#$1

(93)

#

(94)

N
X
"L 1
" $1 ðqÞ"ak
¼C
tk C
k
tL k¼1"L is the applied axial laminate stress, tL is the thickness of the corresponding composite lamr
inates, and DT ð¼ T $ Tsf Þ is the difference between the testing temperature T and stress-free1528

International Journal of Damage Mechanics 29(10)temperature Tsf . Therefore, the effective elastic constants and effective thermal expansion coefficients of the laminates are given as1;
EX ¼ "C 11 ðqÞ

C" 12 ðqÞ! XY ¼ $ ";C 11 ðqÞ

1
EY ¼ ";C 22 ðqÞ
aX ¼ "a L1 ðqÞ;1
GXY ¼ "C 66 ðqÞ

aY ¼ "a L2 ðqÞ

(95)
(96)

We can derive the effective elastic constants analytically with the crack density and fiber breaks
through the CLT. The stress–strain relationship of damaged laminates can be obtained using
equation (92).

Results and discussion
Verification of damage models
The stiffness reduction in laminates with transverse cracking was studied using five approaches forverifying the effective compliance/stiffness matrices and proposed damage variables d22 and d12.
Murakami_simplified shows the Murakami model using d22 expressed in a simplified form in equation (64). Yoshimura_simplified is the Yoshimura model with d22 expressed in a simplified form,
identical to Murakami_simplified. Li_simplified is the Li model including the damage parameterx in equations (16) and (64) and constant damage parameter n in equations (17) and (47).
Maim!ı_simplified gives the Maim!ı model with simplified forms of d22 in equation (64) and d12
in equation (85). Maim!ı_infinite_series describes the result of the Maim!ı model with
infinite forms of d22 in equation (65) and d12 in equation (86). In Li_simplified, the upper limit of
the transverse crack density resulting from the assumption of a small damage postulated in Li’s
effective stiffness matrix was determined in the same manner as in our previous study (Li et al.,2019). Table 1 presents the upper limit of the transverse crack density in GFRP laminates usedin Li_simplified.
For verification, comparisons to results obtained using FEA are presented for ½90=0,s GFRP
cross-ply laminates. The material properties listed in Table 2 were used. Transverse cracks are!
assumed to exist in the 90 plies and face the outer surface. Gudmundson and Zang (1993) calculated the reduction in Poisson’s ratio and Young’s modulus due to the increase in the transverse
crack density. When the matrix crack faces the surface, the crack length should be doubled toinclude the surficial effect and satisfy the boundary condition. Therefore, we obtaind22 ¼ 1 $

d22 ¼

4q
b
tanh
b
4q

(97)

d12 ¼

4q
b
tanh 12
4q
b12

(98)

1
4 X
1
tanh½2ð2n $ 1Þpktq,
33
p kt n¼1 ð2n $ 1Þq(99)

Onodera and Okabe

1529

Table 1. Upper limit of transverse crack density used in Li_simplified.
LaminteUpper limit of transverse crack density (cracks/mm)

GFRP ½90=0,s
GFRP ½55=$55,N

2.01
4.02

Table 2. Material constants and thickness of GFRP ply used in stiffness reduction calculation Gudmundson andZhang (1993).
E1 (GPa)
41.7

E2 (GPa)
13

! 12
0.3

! 23
0.42

G12 (GPa)3.4

d12 ¼ 1 $G23 (GPa)
4.58

a1 (! C)6.72-10a2 (! C)$6

29.3-101
4 X
1tanh½2ð2n $ 1Þptq,p3 t n¼1 ð2n $ 1Þ3q

Ply thickness (mm)
$6

0.203

(100)

The comparison results are shown in Figures 3 and 4. As for Young’s modulus, all approaches
reproduce the tendency of the results observed in FEA. When comparing Maim!ı_simplified and
Maim!ı_infinite_series, it is shown that the simple forms of equations (64) and (85) yielded results
closer to the FEA results than the infinite series forms of equations (65) and (86). Thus, equations(65) and (86) gave smaller crack openings than the FEA results, and loosening the boundary condition may improve the tendency of the crack opening displacement to match with those obtainedusing FEA. Murakami_simplified, Li_simplified, Maim!ı_simplified, and Maim!ı_infinite_series captured the degradation tendency of the Poisson’s ratio; however, Yoshimura_simplified obviously
underestimated the degradation. This causes the significant errors in the case of angle-ply laminates
because the non-diagonal component directly affects the Young’s modulus of the laminate.
Gudmundson and Zhang (1993) also calculated the stiffness reduction in several types of angleply laminates using FEA. The present model was verified by comparison with the results for ½(55,Nangle-ply GFRP laminates. The material properties listed in Table 2 were used. Figures 5 to 10 show
the changes in the Young’s moduli, shear modulus, Poisson’s ratio, and thermal expansion coefficients of ½(55,N laminate as a function of the transverse crack density. Clearly, Maim!ı_simplified
shows the best agreement among these five approaches for any of the elastic constants, whereas
Yoshimura_simplified cannot capture the trends of Young’s modulus and Poisson’s ratio at all. As
stated above, Poisson’s ratio is a critical parameter for evaluating the properties of angle-ply
laminates as it reflects the behavior of non-diagonal components in compliance/stiffness matrix.
The trends of Yoshimura_simplified are substantially different from those observed in FEA.
Therefore, we do not recommend using Yoshimura_simplified for evaluating angle-ply laminates.
Although the results obtained using all the models agree well with the FEA results,
Maim!ı_simplified shows the best agreement. In fact, we proved that Maim!ı_simplified is the mostappropriate approach among five approaches considering in this study.

Nonlinear stress–strain response of CFRP cross-ply laminates
The nonlinear stress–strain response of CFRP cross-ply laminates resulting from transverse cracking was calculated using the present model and compared to the experimental results obtained inthis study. We also applied our model to the results reported by Fikry et al. (2018). The stress–strain

1530

International Journal of Damage Mechanics 29(10)Young’s modulus EX (GPa)28

FEA result
Murakami_simplified
Yoshimura_simplified
Li_simplified
Maimí_simplified
Maimí_infinite_series

26
24
22
20

0

2

4
6
8
Crack density ρ (/mm)

10Figure 3. Comparison of Young’s modulus reduction in GFRP ½90=0,s laminate due to transverse crack predictedusing present models and FEA (Gudmundson and Zhang (1993)).

0.16

FEA result
Murakami_simplified
Yoshimura_simplified
Li_simplified
Maimí_simplified
Maimí_infinite_series

Poisson’s ratio νXY

0.14
0.12
0.1
0.08
0.06

0

2

4
6
8
Crack density ρ (/mm)

10Figure 4. Comparison of Poisson’s ratio reduction in GFRP ½90=0,s due to transverse cracking predicted usingpresent models and FEA (Gudmundson and Zhang (1993)).

Young’s modulus EX (GPa)12

FEA result
Murakami_simplified
Yoshimura_simplified
Li_simplified
Maimí_simplified
Maimí_infinite_series

8
4
0

0

2
4
6
8
Crack density ρ (/mm)

10Figure 5. Comparison of Young’s moduli of GFRP ½(55,N axial laminates due to transverse cracking predicted bypresent models and FEA (Gudmundson and Zhang (1993)).

Onodera and Okabe

1531

Young’s modulus EY (GPa)

18
16
14
12
10

FEA result
Murakami_simplified
Yoshimura_simplified
Li_simplified
Maimí_simplified
Maimí_infinite_series

8
6
4

0

2
4
6
8
Crack density (/mm)

10Figure 6. Comparison of Young’s moduli reduction in GFRP ½(55,N transverse laminate due to transverse cracking
predicted by present models and FEA (Gudmundson and Zhang (1993)).

Shear modulus GXY (GPa)12

FEA result
Murakami_simplified
Yoshimura_simplified
Li_simplified
Maimí_simplified
Maimí_infinite_series

11
10
9
8

0

2

4

6

Crack density (/mm)

810Figuer 7. Comparison of shear modulus reduction in GFRP ½(55,N laminate due to transverse cracking predicted bypresent models and FEA (Gudmundson and Zhang (1993)).

response of cross-ply laminates with transverse cracking was calculated using equations (92) to (94).
From the verification result described in the previous section, the effective stiffness model ofMaim!ı_simplified composed of the Maim!ı model with simplified forms of d22 in equation (64)and d12 in equation (85) was utilized to predict the stiffness reduction due to transverse cracks in!
90 plies. In CFRP laminates subjected to quasi-static monotonic tensile loading, fiber fractureoften occurs at the ultimate failure stress of the laminate. Therefore, the stiffness reduction
owing to fiber breaks was not considered here, and the ultimate tensile strength of CFRP crossply laminates was predicted based on a strain-based fiber failure criterion presented by McCartney(2005) to determine the fracture point in the nonlinear stress–strain curve. In this criterion,
the calculated effective laminate axial strain using equation (92) is considered equal to the fiber!
axial strain in 0 plies. When the effective laminate axial strain is equal to the critical strain, efc ,during ultimate fiber failure, ultimate failure of the laminate occurs. The ultimate fiber failure
strain, efc , is derived from the spring element model (SEM) simulation (Okabe et al., 2005),
which considers the surficial stress concentration obtained from multiple-fiber fragmentation
tests (Yamamoto et al., 2019). The bimodal Weibull distribution was applied to reproduce
the statistically distributed strength in fibers. A stress concentration factor of 1.9 was used.

1532

International Journal of Damage Mechanics 29(10)Poisson’s ratio XY

0.5
0.45
FEA result
Murakami_simplified
Yoshimura_simplified
Li_simplified
Maimí_simplified
Maimí_infinite_series

0.4
0.35

0

2
4
6
8
Crack density (/mm)

10Thermal expansion coefficient αX( 10-6 /℃)

Figure 8. Comparison of Poisson’s ratio change in GFRP ½(55,N laminate due to transverse cracking predicted bypresent models and FEA (Gudmundson and Zhang (1993)).
25
20

15
FEA result
Murakami_simplified
Yoshimura_simplified
Li_simplified
Maimí_simplified
Maimí_infinite_series

10
5

0

2

4

6

8

10

Crack density ρ (/mm)

Thermal expansion coefficient( 10-6 /℃)

αYFigure 9. Comparison of axial laminate thermal expansion coefficient reduction in GFRP ½(55,N due to transverse
cracking predicted by present models and FEA (Gudmundson and Zhang (1993)).
8

6
4
FEA result
Murakami_simplified
Yoshimura_simplified
Li_simplified
Maimí_simplified
Maimí_infinite_series

2
0

0

2

4

6

8

Crack density ρ (/mm)

10Figure 10. Comparison of transverse laminate thermal expansion coefficient reduction in GFRP ½(55,N due totransverse cracking predicted by present models and FEA (Gudmundson and Zhang (1993)).

The average ultimate fiber failure strain obtained by 100 SEM simulations scaled with each specimen size was regarded as the ultimate fiber failure strain efc . Table 3 lists the material properties of
the T700S and resin, and Table 4 shows the bimodal Weibull parameters of the T700S fiber used in
SEM simulations. The evolution of transverse crack density q as a function of the applied laminate
stress is needed to predict the nonlinear stress–strain curves of the laminates. In this study, toOnodera and Okabe

1533

Table 3. Properties of T700S fiber and resin used in SEM simulations(Okabe et al., 2015).
Fiber axial Young’s modulus Ef (GPa)Fiber radius r (mm)Matrix Young’s modulus Em (GPa)Matrix Poisson’s ratio ! m (–)

230
3.5
3.4
0.31

Table 4. Bimodal Weibull parameters of the T700S fiber usedin SEM simulations ( Yamamoto et al., 2019).
r01 (GPa)m1

r02 (GPa)

m2

5.2

4.8

6.1

12.0

validate the proposed stiffness reduction model, we assumed the following Weibull form for transverse crack density evolution
"
( !
" )#
" L $ rL0 mTrq ¼ qs 1 $ exp $
(101)
rLc
where qs is the saturated transverse crack density, rL0 is the initiation stress of the transverse crack,rLc is the Weibull scale parameter, and mT is the Weibull modulus. These parameters were determined by fitting with the experimental results.
We performed the experiments by ourselves to validate our model presented in this study. Here,
we briefly explain the experimental procedures followed to determine the stiffness reduction, transverse crack density, and stress–strain curves of the cross-ply laminates. Loading and unloading tests
were conducted to obtain the nonlinear stress–strain response of cross-ply laminates, transverse!
crack density in 90 plies, and stiffness reduction. The ½0=908 =0, cross-ply laminate was composed of
T700S/2592 (Toray Industries) unidirectional plies. The laminate specimen had a length, width, andthickness of 249, 25, and 1.42 mm, respectively. A strain gauge (Kyowa Electronic Instruments Co.)was used to measure the strain in the laminate. The free edge of the laminate was polished bya diamond grinding disk (Buehler Ltd.) to observe the transverse crack density. Loading andunloading tests were performed using a servo-hydraulic testing machine MTS810 (MTS Systems
Corporation) and involved the following steps: (1) the specimen was loaded with a crosshead speed
1.0 mm/min; (2) when the load reached a specific value, the specimen was held at this load and thetransverse cracks were observed at the free edge of the specimen using RepliSet (Struers); (3) the
specimen was then unloaded with a crosshead speed 1.0 mm/min until the load reached zero. Thisloading/unloading cycle comprising steps (1) to (3) was repeated until the specimen failed.
The stress–strain relationship for the T700S/2592 ½0=908 =0, laminate predicted from the present
model was compared to the obtained experiment results. We denote T700S/2592 as CFRP-I. The
material properties of the ply, parameters of transverse crack evolution, and ultimate fiber failure
strain, efc , derived from SEM simulations are listed in Tables 5 to 7, respectively. The out-of-plane
shear modulus, G23, of CFRP-I ply given in Table 5 was derived as follows
G23 ¼

E2
2ð1 þ ! 23 Þ

(102)

1534

International Journal of Damage Mechanics 29(10)Table 5. Mechanical properties of plies used in the calculation of stress-strain behavior of cross-ply laminates.
CFRP-I (Kumagai et al., 2020) CFRP-II (Fikry et al., 2018)
115
Longitudinal Young’s modulus E1 (GPa)Transverse Young’s modulus E2 (GPa)7.74
In-plane Poisson’s ratio ! 12 (–)
0.315
0.394
Out-of-plane Poisson’s ratio ! 23 (–)
In-plane shear modulus G12 (GPa)3.54
Out-of-plane shear modulus G23 (GPa)
3.78a
$6 !
Longitudinal thermal expansion coefficient a1 (10 / C) 0.4bTransverse thermal expansion coefficient a2 (10$6 /! C) 36b
Temperature difference DT (! C)–100
Fiber volume fraction f (–)
0.51
Ply thickness (mm)
0.142

123
8.68
0.33
0.394c
3.92
2.91a
0.4b
36b
–110
0.528d0.15

a

G23 was estimated using equation (102).
These values are quoted from Okabe et al. (2015).
cThis value is quoted from Kumagai et al. (Kumagai et al., 2020).
dThis was estimated from equation (103). Ef and Em given in Table 3 were used.
b

Table 6. Simulation results of ultimate fiber failure strain at stressconcentration factor of 1.9.

Ultimate fiber failure strain efc (%)

CFRP-I

CFRP-II

1.77

1.83

Table 7. Parameters of transverse crack evolution in
CFRP laminate used in equation (102).

rL0 (MPa)rLc (MPa)mT
qs (mm)

CFRP-I

CFRP-II

180
74
1.5
0.76

180
110
1.7
0.9

Parameters of transverse crack density evolution were determined by fitting the calculated resultsof the present model with the experimental data as shown in Figure 11(a). Figure 11(b) presents acomparison of the predicted results for stiffness reduction as a function of transverse crack density
determined by using the present model with experimental results. The results predicted using thepresent model agree quantitatively with the experimental results obtained through loading/unloading tests. The predicted stress–strain responses of the T700S/2592 ½0=908 =0, laminate were comparedto the experimental results, as shown in Figure 11(c). The results of the present model obtained
considering stiffness reduction due to transverse cracking show good agreement with the experimental data. According to these results, the nonlinearity of the stress–strain responses observed in
the 200–350 MPa range is mainly due to transverse cracking. The ultimate failure stress andstrain determined from the critical-strain-based criterion also agree with the experimental data.

Onodera and Okabe

1535

(a)

(b)
Normalized Young's modulus (GPa)Crack density (cracks/mm)

0.8

Experiment

0.7

Present model

0.6
0.5
0.4
0.3
0.2

CFRP T700S/2592[0/908/0]

0.1
0

0

100

200

300

400

500

1.05

Experiment1

Present model

0.95
0.9
0.85
0.8

Stress (MPa)

CFRP T700S/2592[0/908/0]
0

0.2

0.4

0.6

Crack density (cracks/mm)

0.8

(c)

Experiment
Present model

Figure 11. Comparison of results predicted using present model for (a) transverse crack density evolution, (b)stiffness reduction, and (c) stress–strain response of T700S/2592 ½0=908 =0, laminate with experimental results.
!

This result implies that the ultimate failure strength of CFRP laminates with 0 plies can be determined using the strain-based fiber criterion, and the ultimate fiber failure strain used in this criterion
can be predicted through SEM simulations.
The nonlinear stress–strain behavior of the T700SC/2592 ½0=908 =0, laminate was predicted using
the present model, and the results were compared to the experimental result obtained by Fikry et al.
(2018). We denote T700SC/2592 as CFRP-II. The material properties of the ply, parameters oftransverse crack evolution, and ultimate fiber failure strain, efc , derived from SEM simulations arealso listed in Tables 5 to 7, respectively. The out-of-plane shear modulus G23 of CFRP-II ply
given in Table 5 was estimated using equation (102). The fiber volume fraction, f, of CFRP-II isestimated as
f¼

E1 $ Em
Ef $ Em(103)

where Ef is the fiber axial Young’s modulus and Em is the matrix Young’s modulus given in Table 3.
The crack evolution parameters used in equation (101) were derived by fitting the predicted results

1536

(b)

1

Normalized Young's modulus

Crack density (cracks/mm)(a)International Journal of Damage Mechanics 29(10)Experiment

0.8

Present model

0.6
0.4

CFRP T700SC/2592[0/908/0]

0.2
0

0

100

200

300

400

500

600

1
0.95
0.9
0.85
0.8

Experiment
Present model

CFRP T700SC/2592[0/908/0]

0.2

0.6

0.75
0.7

0

Stress (MPa)
(c)

0.4

0.8

1

Crack density (cracks/mm)

600

Stress (MPa)

500
400
300
200

CFRP T700SC/2592[0/908/0]

100
0

-0.5

0

0.5

1

1.5

Experiment(Longitudinal strain)
Experiment(Transverse strain)
Present model(Longitudinal strain)
Present model(Transverse strain)

2

Strain (%)
Figure 12. Comparison of (a) transverse crack density, (b) stiffness reduction, and (c) stress-strain response of
T700SC/2592 ½0=908 =0, predicted using the present model with experimental results (Fikry et al. (2018)).

with the experiment results, as shown in Figure 12(a). The predicted results for stiffness reduction as
a function of the transverse crack density show excellent agreement with the experimental results, as
shown in Figure 12(b). Figure 12(c) presents the nonlinear stress as a function of the longitudinal ortransverse strain obtained using the present model and through experiments. The present model
considering the stiffness reduction due to transverse cracks can predict the nonlinearity of thestress–strain response up to the ultimate failure strength of the laminate estimated using thecritical-strain-based criterion.

Nonlinear stress–strain response of GFRP cross-ply laminates
The nonlinear stress–strain behavior of E-glass/TX24235 GFRP ½0=90n =0,s ðn ¼ 2; 4; 6Þ laminates
was predicted to validate our proposed damage variables d11, d22, and d12. The stress–strain
responses of cross-ply laminates with transverse cracking and fiber breaks were calculated using
equations (92) to (94). The effective stiffness model Maim!ı_simplified comprising of the Maim!ı
model with simplified forms of d22 in equation (64) and d12 in equation (85) was utilized to predict!
the stiffness reduction due to transverse cracks in 90 plies. In the case of GFRP laminates, the!
stiffness reduction due to multiple-fiber fragmentation in 0 plies and subsequent final laminate
failure was calculated using damage variable d11 in equation (28).

Onodera and Okabe

1537

The GFRP ply properties were estimated using the analytical model developed by Kravchenko
et al. (2017). According to them, the properties of unidirectional-fiber-reinforced ply composed of
the transversely isotropic fibers and the isotropic matrix are determined as follows"#ðfÞ4ð! ðmÞ $ ! 12 Þ2 kðmÞ kðfÞ GðmÞðfÞ ðfÞðmÞ ðmÞ
(104)f ðfÞ f ðmÞE1 ¼ E1 f þ E f þðkðfÞ þ GðmÞ ÞkðmÞ þ ðkðfÞ $ kðmÞ ÞGðmÞ f ðfÞ
*

1
1
!2
þ
þ 12
4k2 4G23 E1

E2 ¼
ðfÞ
! 12 ¼ ! 12 f ðfÞ þ ! ðmÞ f ðmÞ þ"

G12 ¼ GðmÞ

G23 ¼ G

"

"

E2
E2
$ 2! 212
2k2
E1

ðfÞ

ðfÞ

ðfÞ

ðfÞ

ðG12 þ GðmÞ Þ $ ðG12 $ GðmÞ Þf ðfÞ
ðfÞ

#

ðfÞ

ðG23 þ GðmÞ ÞkðmÞ þ 2G23 GðmÞ þ ðG23 $ GðmÞ ÞkðmÞ f ðfÞ
ðfÞ

ðfÞ

ðG23 þ GðmÞ ÞkðmÞ þ 2G23 GðmÞ $ ðG23 $ GðmÞ ÞðkðmÞ þ 2GðmÞ Þf ðfÞ
a1 ¼

ðfÞ

ðfÞ

(106)

(107)

ðG12 þ GðmÞ Þ þ ðG12 $ GðmÞ Þf ðfÞ

ðfÞ

ðfÞ

(105)

#
ðfÞ
ð! ðmÞ $ ! 12 ÞðkðmÞ $ kðfÞ ÞGðmÞ
f ðfÞ f ðmÞðkðfÞ þ GðmÞ ÞkðmÞ þ ðkðfÞ $ f ðmÞ ÞGðmÞ f ðfÞ! 23 ¼ 1 $
ðmÞ

+$1

ðfÞ

a1 E1 þ aðmÞ EðmÞ f ðmÞ

(108)
#

(109)

(110)

ðfÞ

E1 f ðfÞ þ EðmÞ f ðmÞðfÞ ðfÞ

ðfÞ

(111)

# $2
ðfÞ
ðfÞ
ðfÞ
ðfÞ
2ð1 $ ! 23 ÞE1 $ 4 ! 12 E2(112)

a2 ¼ ða2 þ ! 12 a1 Þf ðfÞ þ ð1 þ ! ðmÞ ÞaðmÞ f ðmÞ $ ð! 12 f ðfÞ þ ! ðmÞ f ðmÞ Þa1
wherekðfÞ ¼kðmÞ ¼

k2 ¼

ðfÞ

ðfÞ

E1 E2

EðmÞ
2ð1 $ ! ðmÞ Þ $ 4ð! ðmÞ Þ

2ðkðfÞ þ GðmÞ ÞkðmÞ þ ðkðfÞ $ kðmÞ ÞGðmÞ fðfÞðkðfÞ þ GðmÞ Þ $ ðkðfÞ $ kðmÞ ÞfðfÞ

(113)

(114)

Here, superscripts (f) and (m) denote the fiber and matrix property, respectively. Subscripts 1, 2,
and 3 represent the material coordinate axis of the ply shown in Figure 1. E, G, !, a, k, and f are1538

International Journal of Damage Mechanics 29(10)Table 8. Material properties of E-glass fiber and resin (Malik, 2011).
Fiber Young’s modulus Eðf Þ (GPa)Fiber Poisson’s ratio ! ðf Þ
Fiber shear modulus Gðf Þ (GPa)Fiber thermal expansion coefficient aðf Þ (10$5 /! C)
Fiber volume fraction f
Matrix Young’s modulus EðmÞ (GPa)Matrix Poisson’s ratio ! ðmÞ
Matrix shear modulus GðmÞ (GPa)Matrix thermal expansion coefficient aðmÞ (10$5 /! C)

72.4
0.2
30.17
0.5
0.5
3.4
0.35
1.259
11.0

Table 9. Material properties used for calculating d11 (Zhao et al.,1974).
Fiber radius r (mm)Gauge length L0 (mm)Weibull modulus m
Fiber strength r0 based on L0 (MPa)Fiber/matrix interfacial shear stress s (MPa)

6.5
24
6.34
1550
25

Table 10. Parameters of transverse crack evolution in E-glass/
TX24235 GFRP ½0=90n =0,s (n ¼ 2, 4, and 6) laminates.
rL0 (MPa)rLc (MPa)mT
qs (mm)½0=902 =0,s

½0=904 =0,s

½0=906 =0,s

140
100
3.1
21.5

70
100
3.1
16.0

55
100
3.1
10.5

the Young’s modulus, shear modulus, Poisson’s ratio, thermal expansion coefficient, bulk modulus,
and volume fraction, respectively. k2 is the transverse plane bulk modulus.
ðfÞ
ðfÞ
ðfÞ
ðfÞ
ðfÞ
ðfÞ
ðfÞ
ðfÞ
E1 ¼ E2 ¼ EðfÞ ; G12 ¼ G23 ¼ GðfÞ ; ! 12 ¼ ! 23 ¼ ! ðfÞ , and a1 ¼ a2 ¼ aðfÞ are assumed because
the E-glass fiber is isotropic. Table 8 lists the material properties of the E-glass fiber and resin
used in equations (104) to (114). Tables 9 and 10 present the material properties used for calculating
d11 and the parameters of transverse crack evolution in equation (101), respectively. The crack
evolution parameters were derived by fitting the results predicted by present model with the experimental results, as shown in Figure 13(a).
Figure 13(b) shows the comparison of the stress–strain responses of E-glass/TX24235 GFRP
½0=90n =0,s ðn ¼ 2; 4; 6Þ laminates with those obtained experimentally by Okabe et al. (2004). The
results predicted by the present model for the ½0=902 =0,s laminate does not agree with the experiment results. The present model overestimated the laminate failure strength, despite considering thestiffness reduction due to transverse cracking and fiber breaks. The reason for this overestimation
may be the stiffness reduction, due to additional damage modes such as delamination resulting fromstress concentration at the transverse crack tip. The predicted results for ½0=904 =0,s and ½0=906 =0,slaminates reasonably agree with the experimental results. Furthermore, the laminate failure strength
derived from the GLS model can be used to predict the results. From the results described in theprevious and current sections, we concluded that the present model can reproduce the stress–strain
responses of cross-ply laminates with transverse cracks and fiber breaks using the Maim!ı model withdamage variables d11, d22, and d12 proposed in this study.

1539

(a)

25

Transverse crack density(cracks/cm)

Onodera and Okabe

20

n=2
n=4

15

n=6

10
5
0

E-glass/TX24235[0/90n/0]0(b) 500

1
1.5
Laminate strain (%)

2.5

n=2
n=4

300

n=6

200
100
0

2

Experiment
Present model

400
Stress (MPa)

0.5E-glass/TX24235[0/90n/0]01Strain (%)

2

3

Figure 13. Comparison of (a) transverse cracking progression and (b) stress-strain responses of E-glass/TX24235
GFRP ½0=90n =0,s (n ¼ 2, 4, and 6) laminates predicted using the present model with experimental results (Okabe et al.
(2004)).

Conclusions
This study involved deriving a new analytical model for predicting the effective compliance ofcomposite laminates and simplified expressions for the damage tensor. We evaluated four models
that yielded asymmetric/symmetric compliance matrices of damaged composite plies. Damage
variable d11 describing stiffness reduction due to fiber breaks and its upper bound dmax
11 given by
1=ðm þ 2Þ were formulated using the GLS model. Moreover, the simplified analytical forms ofdamage variables d22 and d12 were developed, assuming that the displacement field along the longitudinal direction could be expressed as a quadric function by loosening the boundary condition
for the governing differential equation. We also proved that the effective laminate properties could
be obtained by substituting the damage variables into the CLT- and CDM-based models.
Comparisons of the effective stiffness of laminates with transverse cracking to the results obtained
via finite element calculations were performed for model verification. We found that the Maim!ı
model (Maim!ı et al., 2007) in combination with the proposed formulation for the damage variables,referred to as Maim!ı_simplified in this study, was the most appropriate model for predicting or
estimating damage in composite laminates with transverse cracking.
The nonlinear stress–strain behaviors of CFRP cross-ply laminates resulting from transverse
cracking were predicted using the present model, and the predicted results are compared with theexperimental results. The effective compliance model with simplified damage variable d22 and d12

1540International Journal of Damage Mechanics 29(10)was utilized to calculate the stiffness reduction of the CFRP laminates. The ultimate tensile strengthof CFRP cross-ply laminates was predicted based on a strain-based fiber failure criterion because
fiber fracture in CFRP laminate occurs at the ultimate laminate failure stress. We assumed that thelaminate failed when the effective axial laminate strain became equal to the ultimate fiber axial!
failure strain in 0 plies. The ultimate fiber axial failure strain was derived from SEM simulations.
The present model quantitatively reproduced the experimental results up to the ultimate failure
strength of the laminate.
We predicted the nonlinear stress–strain behavior of ½0=90n =0,s ðn ¼ 2; 4; 6Þ GFRP cross-ply!
!
laminates with transverse cracks in 90 plies and fiber breaks in 0 plies to validate our proposed
damage tensor d11, d22, and d12. d11 associated with fiber breaks and its upper bound describing theultimate failure condition were analytically derived using the GLS model (Curtin, 1991). The effective stiffness model with simple damage variables d11, d22, and d12 was utilized to evaluate thestiffness reduction due to transverse cracking, fiber breaks, and ultimate failure of the laminate.
In ½0=904 =0,s and ½0=906 =0,s laminates, the results predicted using the developed model showedreasonable agreement with the experimental results when the transverse cracking progression was
fitted with the experimental results. Thus, if transverse cracking progression is formulated properly,
the stress–strain diagram can be reproduced better.
Acknowledgements
We would like to thank S Date (Tohoku University) for his support in conducting the SEM simulations. We
would also like to express our appreciation to S Kochii (Tohoku University) for his assistance while performing
experiments. The authors would like to thank Editage (www.editage.com) for English language editing.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors disclosed receipt of the following financial support for the research, authorship, and/or publication: Council for Science, Technology and Innovation (CSTI), Cross-ministerial Strategic
Innovation Promotion Program (SIP), “Materials Integration” for revolutionary design system of structural
materials (Funding agency: JST); University of Washington-Tohoku University: Academic Open Space (UWTU: AOS); Toray Industries, Inc.; JSPS KAKENHI Grant No. JP 18J20899.

ORCID iD
Sota Onodera

https://orcid.org/0000-0002-1313-2921

References
Allen D, Harris C and Groves S (1987a) A thermomechanical constitutive theory for elastic composites with
distributed damage-I. Theoretical development. International Journal of Solids and Structures 23(9):1301–1318.
Allen D, Harris C and Groves S (1987b) A thermomechanical constitutive theory for elastic composites with
distributed damage-II. Application to matrix cracking in laminated composites. International Journal of
Solids and Structures 23(9): 1319–1338.

Onodera and Okabe

1541

Bailey JE and Parvizi A (1981) On fibre debonding effects and the mechanism of transverse-ply failure in crossply laminates of glass fibre/thermoset composites. Journal of Materials Science 16(3): 649–659.
Cox HL (1952) The elasticity and strength of paper and other fibrous materials. British Journal of Applied
Physics 3(3): 72–79.
Curtin WA (1991a) Exact theory of fibre fragmentation in a single-filament composite. Journal of Materials
Science 26(19): 5239–5253.
Curtin WA (1991b) Theory of mechanical properties of ceramic-matrix composites. Journal of the American
Ceramic Society 74(11): 2837–2845.
Dvorak GJ and Laws N (1987) Analysis of progressive matrix cracking in composite laminates II. First ply
failure. Journal of Composite Materials 21(4): 309–329.
Dvorak GJ, Laws N and Hejazi M (1985) Analysis of progressive matrix cracking in composite laminates I.
Thermoelastic properties of a ply with cracks. Journal of Composite Materials 19(3): 216–234.
Fikry MM, Ogihara S and Vinogradov V (2018) The effect of matrix cracking on mechanical properties in FRP
laminates. Mechanics of Advanced Materials and Modern Processes 4(1): 3.
Garrett KW and Bailey JE (1977a) Multiple transverse fracture in 90˚ cross-ply laminates of a glass fibrereinforced polyester. Journal of Materials Science 12(1): 157–168.
Garrett KW and Bailey JE (1977b) The effect of resin failure strain on the tensile properties of glass fibrereinforced polyester cross-ply laminates. Journal of Materials Science 12(11): 2189–2194.
Gudmundson P and Zang W (1993) An analytic model for thermoelastic properties of composite laminates
containing transverse matrix cracks. International Journal of Solids and Structures 30(23): 3211–3231.
Hajikazemi M and McCartney LN (2018) Comparison of variational and generalized plane strain approaches
for matrix cracking in general symmetric laminates. International Journal of Damage Mechanics 27(4):507–540.
Han Y, Hahn H and Croman R (1988) A simplified analysis of transverse ply cracking in cross-ply laminates.
Composites Science and Technology 31(3): 165–177.
Hashin Z (1986) Analysis of stiffness reduction of cracked cross-ply laminates. Engineering Fracture Mechanics
25(5–6): 771–778.
Hashin Z (1987) Analysis of orthogonally cracked laminates under tension. Journal of Applied Mechanics 54(4):872–879.
Hui CY, Phoenix S, Ibnabdeljalil M, et al. (1995) An exact closed form solution for fragmentation of Weibull
fibers in a single filament composite with applications to fiber-reinforced ceramics. Journal of the Mechanics
and Physics of Solids 43(10): 1551–1585.
Ibnabdeljalil M and Curtin W (1997) Strength and reliability of fiber-reinforced composites: Localized loadsharing and associated size effects. International Journal of Solids and Structures 34(21): 2649–2668.
Kanagawa Y, Murakami S, Liu Y, et al. (1996) Description of internal damage in FRP laminates by continuum
damage mechanics. Journal of the Society of Materials Science, Japan 45(2): 206–211.
Kobayashi S, Ogihara S and Takeda N (2000) Damage mechanics analysis for predicting mechanicalbehavior of general composite laminates containing transverse cracks. Advanced Composite Materials
9(4): 363–375.
Kravchenko OG, Kravchenko SG and Pipes RB (2017) Cure history dependence of residual deformation in athermosetting laminate. Composites Part A: Applied Science and Manufacturing 99: 186–197.
Kumagai Y, Onodera S, Salviato M, et al. (2020) Multiscale analysis and experimental validation of crack
initiation in quasi-isotropic laminates. International Journal of Solids and Structures 193–194: 172–191.
Lee JW, Allen D and Harris C (1989) Internal state variable approach for predicting stiffness reductions infibrous laminated composites with matrix cracks. Journal of Composite Materials 23(12): 1273–1291.
Lee JW and Daniel I (1990) Progressive transverse cracking of crossply composite laminates. Journal of
Composite Materials 24(11): 1225–1243.
Li S, Wang M, Jeanmeure L, et al. (2019) Damage related material constants in continuum damage mechanics
for unidirectional composites with matrix cracks. International Journal of Damage Mechanics 28(5):690–707.

1542

International Journal of Damage Mechanics 29(10)McCartney LN (1992) Theory of stress transfer in a 0˚-90˚-0˚ cross-ply laminate containing a parallel array oftransverse cracks. Journal of the Mechanics and Physics of Solids 40(1): 27–68.
McCartney L (2005) Energy-based prediction of progressive ply cracking and strength of general symmetric laminates using an homogenisation method. Composites Part A: Applied Science and Manufacturing 36(2): 119–128.
Maimı́ P, Camanho P, Mayugo J, et al. (2007) A continuum damage model for composite laminates: Part I –Constitutive model. Mechanics of Materials 39(10): 897–908.
Malik SA (2011) Damage detection using self-sensing composites. PhD Thesis, University of Birmingham, UK.
Manders PW, Chou TW, Jones FR, et al. (1983) Statistical analysis of multiple fracture in 0˚/90˚/0˚ glass fibre/
epoxy resin laminates. Journal of Materials Science 18(10): 2876–2889.
Murakami S (2012) Continuum Damage Mechanics, Solid Mechanics and Its Applications. vol. 185. Dordrecht,
The Netherlands: Springer.
Nairn JA and Hu S (1992) The initiation and growth of delaminations induced by matrix microcracks inlaminated composites. International Journal of Fracture 57(1): 1–24.
Neumeister JM (1993) A constitutive law for continuous fiber reinforced brittle matrix composites with fiber
fragmentation and stress recovery. Journal of the Mechanics and Physics of Solids 41(8): 1383–1404.
Noda J, Okabe T, Takeda N, et al. (2004) Damage process of GFRP cross-ply laminates. Transactions of theJapan Society of Mechanical Engineers Series A 70(698): 1364–1369.
Okabe T, Imamura H, Sato Y, et al. (2015) Experimental and numerical studies of initial cracking in CFRP
cross-ply laminates. Composites Part A: Applied Science and Manufacturing 68: 81–89.
Okabe T, Onodera S, Kumagai Y, et al. (2018) Continuum damage mechanics modeling of composite laminates including transverse cracks. International Journal of Damage Mechanics 2018; 27(6): 877–895.
Okabe T, Sekine H, Ishii K, et al. (2005) Numerical method for failure simulation of unidirectional fiberreinforced composites with spring element model. Composites Science and Technology 65(6): 921–933.
Okabe T, Sekine H, Noda J, et al. (2004) Characterization of tensile damage and strength in GFRP cross-ply
laminates. Materials Science and Engineering A 383(2): 381–389.
Onodera S and Okabe T (2019) Three-dimensional analytical model for effective elastic constants of transversely isotropic plates with multiple cracks: Application to stiffness reduction and steady-state cracking ofcomposite laminates. Engineering Fracture Mechanics 219: 106595.
Pagano N, Schoeppner G, Kim R, et al. (1998) Steady-state cracking and edge effects in thermo-mechanicaltransverse cracking of cross-ply laminates. Composites Science and Technology 58(11): 1811–1825.
Parvizi A and Bailey JE (1978) On multiple transverse cracking in glass fibre epoxy cross-ply laminates. Journalof Materials Science 13(10): 2131–2136.
Parvizi A, Garrett KW and Bailey JE (1978) Constrained cracking in glass fibre-reinforced epoxy cross-ply
laminates. Journal of Materials Science 13(1): 195–201.
Rosen BW (1964) Tensile failure of fibrous composites. AIAA Journal 2(11): 1985–1991.
Talreja R (1985) A continuum mechanics characterization of damage in composite materials. Proceedings of theRoyal Society A: Mathematical, Physical and Engineering Sciences 399: 195–216.
Vinogradov V and Hashin Z (2010) Variational analysis of cracked angle-ply laminates. Composite Science
Technology 70(4): 638–646.
Wang A, Chou P and Lei S (1984) A stochastic model for the growth of matrix cracks in composite laminates.
Journal of Composite Materials 18(3): 239–254.
Xia Z, Carr R and Hutchinson J (1993) Transverse cracking in fiber-reinforced brittle matrix, cross-ply
laminates. Acta Metallurgica Materials 41(8): 2365–2376.
Yamamoto G, Onodera M, Koizumi K, et al. (2019) Considering the stress concentration of fiber surfaces in
the prediction of the tensile strength of unidirectional carbon fiber-reinforced plastic composites.
Composites Part A: Applied Science and Manufacture 121: 499–509.
Yoshimura A, Okabe T, Yamada M, et al. (2011) Damage simulation of CFRP laminates under high velocity
projectile impact. In: 18th International conference on composites materials, Jeju, South Korea, 21–26 August.
Zhao FM, Okabe T and Takeda N (2000) The estimation of statistical fiber strength by fragmentation tests ofsingle-fiber composites. Composites Science Technology 60(10): 1965–1974.