(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

Contents lists available at ScienceDirect

Engineering Fracture Mechanics
journal homepage: www.elsevier.com/locate/engfracmech

Three-dimensional analytical model for eﬀective elastic constantsof transversely isotropic plates with multiple cracks: Application to
stiﬀness reduction and steady-state cracking of composite
laminates

7⁎

Sota Onoderaa, , Tomonaga Okabea,babDepartment of Aerospace Engineering, Tohoku University, 6-6-01 Aramaki-Aza-Aoba, Aoba-ku, Sendai 980-8579, Japan
Department of Materials Science and Engineering, University of Washington, Seattle, WA 98195, USA

A R T IC LE I N F O

ABS TRA CT

Keywords:
Fibre reinforced materials
Polymer matrix composites
Damage mechanics
Micromechanics and/or materials
Crack initiation

An analytical three-dimensional eﬀective elastic constant of transversely isotropic plates that
include ply cracks is proposed using a continuum damage mechanics approach. Two damage
parameters associated with tensile and shear damage are formulated as functions of ply crack
density using local stress ﬁelds that satisfy the equilibrium equations. Three-dimensional laminate theory is then employed to formulate the eﬀective compliance of the laminate using theeﬀective compliance of a damaged ply, and an analytical steady-state cracking model is established. The proposed model reproduces the thermomechanical properties and the crack initiationstress of laminates.

1. Introduction
In recent years, polymer matrix composites that have high speciﬁc strength and rigidity have been widely used in aerospacematerials. Unidirectional ﬁber-reinforced polymer matrix composites are highly anisotropic; therefore, practical applications generally require various types of laminates made by stacking unidirectional ﬁber-reinforced ply, the fracture processes of which must be
known to ensure safe design. The ﬁrst form of damage in ﬁber-reinforced laminates is typically a ply crack [1], which grows totraverse the thickness of the ply and penetrate in the width direction parallel to the ﬁbers in that ply. Although ply cracks are notcritical to the ﬁnal failure of composite laminates, these cracks signiﬁcantly reduce the laminate’s stiﬀness. The initiation of ply
cracking results in stress concentration at the crack tip; therefore more severe damage, such as delamination and/or ﬁber breakage,occurs at the crack tip [2,3]. Hence, the thermomechanical properties of laminates that include ply cracks should be appropriately
formulated to clarify the fracture mechanism of such laminates.
Several models for the stiﬀness reduction of laminates have been developed for cross-ply laminates [4,5], angle-ply laminates[6–8], and general symmetric laminates [9–11]. However, a model for general laminates should also be developed to clarify thefracture mechanism. Also, models for the stiﬀness reduction of laminates having arbitrary conﬁgurations have seldom been formulated.
Continuum damage mechanics (CDM) [12–17] is widely used to formulate the eﬀective stiﬀness of a ply (or laminate) including
ply cracks, and CDM has exhibited favorable compatibility with the laminate theory, which can address arbitrary conﬁgurations.

⁎

Corresponding author.
E-mail address: sonodera@plum.mech.tohoku.ac.jp (S. Onodera).

https://doi.org/10.1016/j.engfracmech.2019.106595
Received 5 April 2019; Received in revised form 17 July 2019; Accepted 30 July 2019
$YDLODEOHRQOLQH$XJXVW
(OVHYLHU/WG$OOULJKWVUHVHUYHG

(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

Nomenclature

λ1, λ2Latin characters

ν
ξ
ρ
ρN
ρkc
σ
σ

proportionality constants between εx , εy , and εza, bA
conversion matrix of stress and strain
C
stiﬀness matrix
shear stiﬀness reduction
Ds
E
Young’s modulus
f (x , y ), g (x , y ) functions used in local stresses of damaged
plyG
Shear modulushhalf of ply widthlhalf of crack spacingNtotal number of plies
elementary matrix that is interchanged in two
P (i, j )rows (or two columns) i and jR
coordinate conversion matrix of the strain
S
compliance matrix
S
compliance matrix for ε and σTcoordinate conversion matrix of the stressthalf of ply thickness
ply thickness of k-th ply
tktL
laminate thickness
Ttesting temperaturestress-free temperature
TsfU
strain energy
displacements in the x-, y-, and z-directionsu, v, w
vCOD
crack opening displacement
coordinate system of the representative volume
x , y, zelement of damaged ply
X , Y , Z coordinate system of the laminate

σ th
σ L,app
σcL, kω

Sub/superscripts•0•1•1 , •2 , •3

related to undamaged property
related to damaged property
related to the longitudinal, in-plane, and out-ofplane transverse directions of ply•23 , •31 , •12 related to the 2–3, 3–1, and 1–2 plane•app
related to applied ply stress or strain of ply in O123
•ave
related to average ply stress or strain of ply in O123•k
components of the k-th ply•Lcomponent of the laminate•max
related to maximum value•I , •O
related to in-plane parts and out-of-plane parts•II , •IO , •OO 3 × 3 compliance submatrix of the 6 × 6 compliance matrix•x , •y , •z related to the 3-, 2-, and 1-directions•xy , •yz , •zx related to the x- y, y -z, and z-x plane•X , •Y , •Z related to the longitudinal, in-plane, and out-ofplane transverse directions of laminate•XY , •YZ , •ZX related to the X-Y , Y -Z, and Z-X plane

Greek characters

α
α
γ
Γk

Γc
∆
∆T
ε
ε
θk

material constants of the Laplace equation for
displacement v
Poisson’s ratio
damage parameter related to shear loading
ply crack density
normalized ply crack density
critical ply crack density
stress vector
stress vector divided into in-plane and out-of-planeparts
stress to cancel out the thermal residual stress of
the k-th ply
uniaxial monotonic tensile stress of the laminate
critical applied laminate stress when the cracks
propagate in k-th ply
damage parameter related to tensile loading

thermal expansion coeﬃcient vector
thermal expansion coeﬃcient vector divided into
in-plane and out-of-plane parts
engineering shear strainenergy release rate associated with ply crack of kth ply
critical energy release ratematerial constant used in the stiﬀness matrix
temperature change
strain vector
strain vector divided into in-plane and out-ofplane partsﬁber angle of k-th ply between the 1- and X- axis

Abbreviations
CDM
CFRP
FEA
GFRP
LRAM
NPL

continuum damage mechanics
carbon ﬁber-reinforced plasticﬁnite-element analysis
glass ﬁber-reinforced plastic
Large Radius Axisymmetric Damage Model
National Physical Laboratory

Kachanov ﬁrst proposed CDM [18]. Later, Allen et al. [12] and Talreja [16] applied CDM to ﬁber-reinforced composites. In CDM
models, damage parameters are used to represent the extent of damage in a material. The damage parameters can be formulatedusing the average crack opening displacement based on linear fracture mechanics [19] or the local stress ﬁeld model in a ply
including ply cracks [14,15]. Okabe et al. [15] formulated a two-dimensional eﬀective compliance matrix of laminates including ply
cracks based on the damage tensor given by Murakami [13]. However, this compliance matrix is not symmetric, and does not allow
all properties of a damaged laminate including out-of-plane properties to be determined. A three-dimensional eﬀective compliance(or eﬀective stiﬀness) matrix is absolutely essential for completely determining the thermoelastic properties of damaged laminate.
Talreja [16] characterizes eight material constants in the three-dimensional eﬀective stiﬀness matrix of transversely isotropic
composite materials with ply cracks by formulating the Helmholtz free energy for isothermal small deformation and small damage
based on vector damage variables. Li et al. [20] associated these material constants derived by Talreja with elastic constants oftransversely isotropic materials and two damage parameters. They conducted a parametric study of crack geometry and distribution(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabein a ﬁnite element-based study to determine the sensitivity of the eﬀective stiﬀness matrix of a unidirectional composite plate with anelliptical crack on the damage parameters; however, the formulation of these damage parameters was not explicitly derived.
Therefore, this study will attempt to determine formulations of these damage parameters analytically to accomplish the formulationof a three-dimensional symmetric eﬀective stiﬀness matrix of transversely isotropic material with ply cracks.
In this study, we developed a stiﬀness reduction model of transversely isotropic composite plates with ply cracks based on the
CDM approach formulated by Li et al. The two damage parameters were formulated as a function of ply crack density using local
stress ﬁeld models that satisfy the equilibrium equations subjected to tensile loading normal to the ﬁber and in-plane shear loading.
These two damage parameters are implemented into the eﬀective stiﬀness matrix of ply derived by Li et al. The three-dimensionallaminate theory is then employed to describe the stiﬀness reduction of composite laminate with arbitrary lay-up conﬁgurations andply cracks, and we validated this stiﬀness reduction model by comparing its results to experiment results and results of ﬁnite-element
analysis. Finally, the energy release rate associated with ply cracking is formulated using the eﬀective compliance of composite
laminates. An energy-based steady-state cracking model is proposed, and the steady-state cracking stress of cross-ply laminates iscalculated for comparison with previous analytical models.
This paper is organized as follows. Section 2 presents the theoretical model proposed in this study. In Section 3, the present model
is compared with experiments and previous models for stiﬀness reduction and steady-state cracking of composite laminates. Section 4
presents the conclusions of this study. Furthermore, to help readers derive the formulation more easily, the Appendix describes thecoordinate conversion to the laminate coordinate system, the formulation of conversion matrix A for the stress and strain used in thethree-dimensional laminate theory, and the formulation of the three-dimensional laminate theory.

2. Theory
2.1. Three-dimensional eﬀective stiﬀness matrix for a ply with ply cracks
Under the condition of small deformation and small damage, a stiﬀness matrix of unidirectional ply with multiple ply cracks wasderived using a CDM-based model [20]. As indicated in Fig. 1, the 1-axis is the ﬁber direction, the 2-axis is the transverse direction,and the 3-axis is the thickness direction.
The ply crack planes are assumed to be parallel to the 1–3 plane, and these cracks penetrate in the ﬁber direction and thethickness direction. According to Li et al., the relationship between the applied average ply stress σ app and the applied ply strain ε app
not considering the thermal residual strain in the damaged ply in Fig. 1 iswhere

σ app = Cε app,

(1)

app
app
app T
σ app = [ σ1app σ2app σ3app σ23
] ,
σ13
σ12

(2)

app
app
app
γ13
γ12
ε app = [ ε1app ε2app ε3app γ23
] ,
T

(3)

and the stiﬀness matrix of damaged ply C is the sum of the stiﬀness matrix C 0 of undamaged ply and the damage-related stiﬀness
matrix C1, expressed as

C = [Cij] = C 0 + C1,(4)

Fig. 1. Schematic of multiple ply cracks parallel to the ﬁber direction.
(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe⎡ 1 − 0ν23 E10
⎢ ∆
⎢
⎢
⎢
0
C 0 = [Cij ] = ⎢
⎢
⎢
⎢
⎢
⎢
⎣ sym.
0

∆0

0
ν12

E20

0 ν0
1 − ν12210 )∆0(1 + ν23

∆0

0

0 0

E20

0

0

E20

0

0

E20

0

0

0
G23

0

0 + ν0 ν0
ν23
12 21
0 )∆0
(1 + ν23E20

0 ν0
1 − ν12210 )∆0(1 + ν23

ν21 E2 ν12 (1 − ν12 ν21)
⎡ E1 ν120 )(∆0)2(1 + ν23⎢ (∆0)2⎢0
0 ν 0 )2E2 (1 − ν1221⎢0 )2 (∆0)2
(1
+
ν
1
23
C1 = [Cij ] = −ω ⎢
⎢
⎢
⎢
⎢
⎢
⎣ sym.
0
0 0 ⎤
0
0 0 ⎥
0
0 0 ⎥
0
⎥
G23
0 0 ⎥,
0)2(1 + ν23
⎥
0 0 ⎥
0⎥
ξG12 ⎦
0 0

where

0
ν12

0

0

0
G12

0 (ν 0 + ν 0 ν 0 )E20 ν12
23
12 21

0 ⎤
⎥
⎥
0 ⎥
⎥
0 ⎥,
⎥
0 ⎥
⎥
0 ⎥
0⎥
G12
⎦

0 )(∆0)2(1 + ν230 ν 0 )(ν 0 + ν 0 ν 0 )
E20 (1 − ν12
21 23
12 21
0 )2 (∆0)2(1 + ν23
0 + ν 0 ν 0 )2
E20 (ν23
12 21
0 )2 (∆0)2(1 + ν23

(6)

0
0 0
∆0 = 1 − ν23− 2ν12
ν21,

0
ν21
=

0
G23
=

(5)

(7)

E20 0
ν12,
E10

(8)

E20
.
0
2(1 + ν23)

(9)

Here, E is Young’s modulus, G is the shear modulus, and ν is Poisson’s ratio. The superscript 0 indicates the undamaged ply, and thesubscripts indicate the coordinate axes of the ply. The parameter ω is associated with tensile damage, and ξ is associated with shear
damage. These damage parameters are expressed as

ω=1−
ξ=

E2
,
E20

(10)

1⎛
G12 ⎞
,
⎜1 −
0 ⎟
ω⎝
G12
⎠

(11)

where E2 is the Young’s modulus of the damaged ply along the 2-axis, and G12 is the shear modulus of damaged ply in the 1–2 plane.
We formulate the damage parameters with the local stress ﬁeld models using the proportion of the average ply strain, which
ignores the crack opening displacement, to the applied ply strain considering crack opening displacement. When average stress σ2app isapplied to the ply along the 2-axis, the average ply strain ε2ave and the applied ply strain ε2app along the 2-axis can be described as

ε2ave =

ε2app =

σ2app
,
E20

(12)

σ2app
.
E2

(13)

Using Eqs. (12) and (13), Eq. (10) is rewritten as

ω=1−

ε2ave
.
ε2app

(14)

app
in the 1–2 plane, the average ply
As for damage parameter ξ , which is a constant, when the ply is subjected to the shear stress σ12
ave
app
γ
γ
engineering shear strain 12 and the applied ply engineering shear strain 12 in the 1–2 plane are written as(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe
ave
γ12
=

app
γ12
=

app
σ12
,
0
G12

(15)

app
σ12
.
G12

(16)

Using Eqs. (15) and (16), Eq. (11) can be reformulated as

ξ=

ave
γ12
1⎛
⎞.
⎜1 −
app ⎟
ω⎝
γ12 ⎠

(17)

Eqs. (14) and (17) can be formulated as functions of ply crack density (i.e., the number of ply cracks per unit length normal to thecrack plane) using the local stress ﬁeld models for ply including cracks subjected to tensile loading along the 2-axis and shear loadingin the 1–2 plane. The eﬀective stiﬀness matrix C is valid when damage parameter ω is very small (i.e., ω2 ≪ 1). The maximum valueof damage parameter ω when the model works eﬀectively will discussed in Section 3.1. The next subsection describes the formulationof the three-dimensional local stress ﬁeld model with the damage parameters.
2.2. Damage parameters
A three-dimensional local stress ﬁeld model of a ply that includes ply cracks was formulated to evaluate stiﬀness reduction. Whenformulating the local stress ﬁeld model, the ply was assumed to be thin, and the damage was assumed to be caused mainly by ply
cracks. Based on these assumptions, the damage due to delamination and ﬁber breakage were ignored.
A three-dimensional local stress ﬁeld model is ﬁrst formulated to subject the stress σ2app along the ply’s 2-axis. Fig. 2 presents a
representative volume element (RVE) that includes ply cracking on both sides of a part of the ply illustrated in Fig. 1.
The coordinate system in Fig. 2 diﬀers from that in Fig. 1: the coordinates x , y , and z in Fig. 2 correspond to the 3-, 2-, and 1-axes
in Fig. 1. The ply cracks are assumed to have tunnel-like crack surfaces that are symmetrical about the y-axis and are postulated not to
propagate into the neighboring ply. The crack spacing is 2l , the thickness of the ply is tk = 2t , and the width of the ply is 2h . The
displacements in the x-, y-, and z-directions are deﬁned as u, v , and w. The ply strain ε2app along the y-axis (or 2-axis) is applied to theRVE (or ply). Because RVE is symmetric, the region of interest in this deﬁnition is limited to 0 ⩽ x ⩽ t , 0 ⩽ y ⩽ l, 0 ⩽ z ⩽ h . The
strain-stress relationship in the RVE is expressed as

εx =

εy =

εz =

ν0
ν0
∂u
σ
= x0 − 230 σy − 210 σz ,
∂x
E2
E2
E2

(18)

σy
ν0
ν0
∂v
= − 230 σx + 0 − 210 σz ,
∂y
E2
E2
E2

(19)

σz
ν0
ν0
∂w
= − 210 σx − 210 σy + 0 ,
∂z
E2
E2
E1

(20)

Fig. 2. Representative volume element subjected to tensile loading.
(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

γxy =

γyz =
γzx =

σxy

0
G23

σyz

0
G12

∂v
∂u
∂v
+
≈
,
∂x
∂y
∂x

=

(21)

∂v
∂w
∂v
+
≈
,
∂z
∂y
∂z

=

(22)

σzx
∂u
∂w
=
+
≈ 0,
0
∂z
∂x
G12

(23)

where σi (i = x , y, z ) is the stress along the i-axis, εi is the strain along the i- axis, σij ( j = x , y, z;j ≠ i ) is the shear stress in the i-j plane,and γij is the engineering shear strain in the i-j plane. For engineering shear strains γxy and γyz , the partial diﬀerential coeﬃcients ∂u/ ∂yand ∂w / ∂y are assumed to be quite small, and the engineering shear strain γzx is approximated as zero. Following Okabe et al. [15], therelationships between the strains εi are assumed as

εx = aεy,(24)

εz = bεy,(25)

where a and b are proportionality constants deﬁned to satisfy the equilibrium equations of stress. The physical meaning of a and b isthe average Poisson’s ratio in the RVE. Using Eqs. (18)–(20), (24) and (25), σx , σy , and σz are expressed as
0 0
0 0
0
0
ν 0 + ν12
ν21 + (1 − ν12ν21) a + ν12(1 + ν23) b ∂vσx = E20 23
,
0
0
0 0
∂y
(1 + ν23)(1 − ν23
− 2ν12ν21)

σy = E20

(26)

0 0
0
0 0
0
0
1 − ν12
ν21 + (ν23+ ν12ν21) a + ν12(1 + ν23) b ∂v,
000 0∂y(1 + ν23)(1 − ν23 − 2ν12 ν21)

(27)

0
0
ν 0 + ν21
a + (1 − ν23) b ∂vσz = E10 21
.
0
0 0
∂y
1 − ν23 − 2ν12 ν21(28)

In addition, using Eqs. (21)–(23), the shear stresses σxy, σyz , and σzx are described as0σxy = G230σyz = G12σzx = 0.

∂v
,
∂x

(29)

∂v
,
∂z

(30)
(31)

The equilibrium equations of stress are expressed as

∂σxy
∂σzx
∂σx
+
= 0,
+
∂y
∂z
∂x

∂σxy
∂x

+

∂σy
∂y

+

∂σyz
∂z

(32)

= 0,

(33)

∂σyz
∂σz
∂σzx
+
= 0.
+
∂y
∂z
∂x

(34)

Assuming that ∂2v / ∂x ∂y ≠ 0 and ∂2v / ∂y ∂z ≠ 0 , and substituting Eqs. (26)–(31) into Eqs. (32)–(34), the proportionality constants aand b and the Laplace equation for displacement v are obtained as

G0
G0 0 ⎛⎡ 00 2
0 ⎞⎤+ 230 {1 − (ν23a = −⎢ν23) } − 120 ν12⎜1 + ν12 ⎟ ,
⎥
E
E1
2
⎝
⎠⎦
⎣

(35)

0
G0 0 ⎛
G12⎡ 0
0⎞⎛⎜1 − ν 0 ν 0 ⎞⎟⎤,
b = −⎢ν21− 230 ν21
⎜1 + ν23⎟ +12 21 ⎥
E2
E10 ⎝
⎝
⎠
⎠⎦
⎣

(36)

∂ 2v
∂ 2v
∂ 2v
+ λ12 2 + λ 22 2 = 0,
∂x 2
∂y
∂z

(37)

where constants λ1 and λ2 are deﬁned as

λ1 =

0 0
0
0 0
0
0
+ ν12
ν21 + (ν23ν21) a + ν12(1 + ν23)b
E20 1 − ν12
,
0
0
0
0 0
G23
(1 + ν23)(1 − ν23 − 2ν12 ν21)

(38)


(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

λ2 =

0
G12
.
0
G23

(39)

To determine displacement v, the boundary conditions of the Laplace equation in Eq. (37) are given by

v = 0 on y = 0,∂v
= 0 on y = l,
∂y

(40)

∂v
= 0 on z = h,
∂z

(42)

v = ε2app y

(44)

∂v
= 0 on x = 0,
∂x

(41)

∂v
= 0 on z = 0,
∂z

(43)

on x = t .

(45)

At y = 0 , no displacement v is considered in Eq. (40). Eq. (41) was determined from Eqs. (26)–(28) considering the stress condition
σx = σy = σz = 0 at the crack surface ( y = l ). Eq. (42) was deﬁned based on Eq. (29), considering the shear stress condition σxy = 0 on
the center plane at x = 0 of the ply. Eq. (43) means that shear stress σyz expressed by Eq. (30) is zero on the surface at z = h , and Eq.
(44) expresses that the shear stress σxy is zero on the center plane at z = 0 of the ply. In addition, the displacement distribution in theinterface, presented as Eq. (45), is assumed based on the previous study [14]. Therefore, the neighboring ply is thought to be
deformed uniformly by mechanical loading ε2app , regardless of the ply cracks. Separating variables and assigning the boundary
conditions presented above to Eq. (37), the solution v (x , y, z ) of the Laplace equation that satisﬁes the boundary conditions can be
expressed as(−1)n + 1 cosh[(2n − 1) πλ1 x /(2l)]8l ⎛2n − 1⎞⎛⎞v ⎜x , y, z ⎟ = 2 ⎜ ∑sin ⎡πy⎤ ⎟ ε2app.
π n = 1 (2n − 1)2 cosh[(2n − 1) πλ1 t /(2l)]
2l
⎣
⎦
⎠
⎝
⎠
⎝
∞

(46)

When Eq. (46) is substituted into Eqs. (26)–(31), the local stress distribution under tensile loading along the 2-axis in the RVE isderived as
0
0 0
0 0
0
0
+ ν12
ν21 + (1 − ν12ν21) a + ν12(1 + ν23)b ⎛4E20 ν23
σx ⎛⎜x , y, z ⎟⎞ =f ⎜x , y ⎟⎞ ε2app,
0
0
0 0
π
+
−
−
(1
ν
)(1
ν
2
ν
ν
)
23
23
12 21
⎝
⎠
⎝
⎠
0 0
0
0 0
0
0
+ ν12
ν21 + (ν23ν21) a + ν12(1 + ν23)b ⎛4E20 1 − ν12σy ⎛⎜x , y, z ⎟⎞ =f ⎜x , y ⎟⎞ ε2app,
0
0
0 0
π
+
−
−
(1
ν
)(1
ν
2
ν
ν
)
23
23
12 21
⎝
⎠
⎝
⎠
0
0
0
+ ν21
a + (1 − ν23)b ⎛4E10 ν21σz ⎛⎜x , y, z ⎟⎞ =f ⎜x , y ⎞⎟ ε2app,
0
0 0
π
−
−
1
ν
2
ν
ν
23
12 21
⎝
⎠
⎝
⎠

(48)

(49)

4 0
σxy ⎛⎜x , y, z ⎞⎟ = G23
λ1 g ⎛⎜x , y ⎟⎞ ε2app,
π
⎝
⎠
⎝
⎠

(50)

σyz (x , y, z ) = 0,

where

(47)

(51)

σzx (x , y, z ) = 0,(52)(−1)n + 1 cosh[(2n − 1) πλ1 x /(2l)]2n − 1⎞⎛f ⎜x , y ⎟ = ∑
cos ⎡
πy⎤,
2n − 1 cosh[(2n − 1) πλ1 t /(2l)]
⎣ 2l
⎦
⎠ n=1
⎝
∞

(−1)n + 1 sinh[(2n − 1) πλ1 x /(2l)]2n − 1⎞⎛g ⎜x , y ⎟ = ∑
sin ⎡
πy⎤.
−
−
2
n
1
cosh[(2
n
1)
πλ
t
/(2
l
)]
2l
⎣
⎦
1
⎠ n=1
⎝
∞

(53)

(54)

Thus, the local stress distribution in a ply that includes ply cracking is formulated. As indicated in Eqs. (46)–(54), the displacement v
and the local stress distribution in a ply are constant with respect to the z-axis, and the shear stresses σyz and σzx are always zero. Whenthe boundary condition given in Eqs. (40)–(45) is utilized, this local stress ﬁeld model is close to the generalized plane strain state.
The compatibility condition for the strain is not satisﬁed because of the assumption used in Eqs. (24) and (25). For carbon ﬁber
reinforced plastic (CFRP), the local stress ﬁeld model almost satisﬁes the plane strain condition because a is much larger than b.
However, a is the same order as b when considering glass ﬁber reinforced plastic (GFRP). Therefore, parameters a and b that are(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabeaverage Poisson’s ratio are needed for GFRP ply. Strictly, a and b are functions of x , y , and z. However a and b are approximated byconstants that satisfy the equilibrium equation of stress in the present model. Therefore, the proposed stress distributions are analytical solutions, but not exact solutions. Damage parameter ω can be formulated using the three-dimensional stress ﬁeld model under
tensile loading along the 2-axis. The average ply strain ε2ave is deﬁned as

ε2ave =

1lth

∫0 ⎛⎜∫0 v ⎛x, l, z ⎞ dx ⎞⎟ dz.
h

t

⎝

⎝

⎠

⎜

⎟

⎠

(55)

By substituting Eq. (46) into Eq. (55), ε2ave/ ε2app is derived as

ε2avetanh[(2n − 1) πλ1 tk ρ/2]
16
1
= 3
,
∑
ε2app
π λ1 tk n = 1 (2n − 1)3
ρ
∞

(56)

where ρ = 1/(2l) is ply crack density, and the thickness of the ply is deﬁned as tk = 2t . Using the above equation, Eq. (14) isreformulated as a function of ply crack density.

ω=1−tanh[(2n − 1) πλ1 tk ρ/2]
161∑
π 3λ1 tk n = 1 (2n − 1)3
ρ
∞

(57)

As illustrated in Fig. 3 (a), a crack with a semi-tunnel-like surface may occur at the laminate surface.
In a thin laminate, surface cracks seriously aﬀect the laminate’s mechanical properties. Fig. 3 (b) depicts an RVE of a ply that has
semi-tunnel-like cracks on the surface. Here, the area delineated by 0 ⩽ x ⩽ 2t , 0 ⩽ y ⩽ l , and 0 ⩽ z ⩽ h should be considered, because the RVE in Fig. 3 (b) is symmetrical with respect to the z-x and x-y planes. The mechanical behavior of a ply with a surface crack
is assumed to be equal to the mechanical behavior of a ply with a ply crack that has twice the length of a surface crack. Based on thisassumption, the surface crack can be considered by replacing ply thickness tk with a thickness of 2tk . Thus, damage parameter ω canbe represented to consider the surface crack by replacing tk in Eq. (57) with 2tk as follows.

ω=1−

8

∑ (2n − 1)3∞π 3λ1 tk n = 1

1tanh[(2n − 1) πλ1 tk ρ]ρ(58)

When the cracks in a ply are considered, Eq. (57) is appropriate; Eq. (58) is appropriate to consider a ply with a surface crack.
Next, damage parameter ξ is derived to consider the damage due to ply cracking parallel to the ﬁber by formulating the local
stress distribution of a ply with ply cracking subjected to shear loading. Fig. 4 depicts an RVE that includes a ply crack on both sides,
which is a part of the ply illustrated in Fig. 1.
app
The x-, y-, and z-coordinates in Fig. 4 correspond to the 3-, 2- and 1-axes in Fig. 1. The in-plane ply shear stress σ12
is applied tothe corresponding ply. The crack distance is 2l , and the ply crack is assumed to have a tunnel-like surface that is symmetrical about
the y-axis. Because the model is symmetric, the range is limited to 0 ⩽ x ⩽ t , 0 ⩽ y ⩽ l, 0 ⩽ z ⩽ h . The deformations are assumed asu (x , y, z ) = v (x , y, z ) = 0,(59)

Fig. 3. (a) Ply including surface cracks. (b) Representative volume element.
(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

Fig. 4. Representative volume element subjected to shear loading.

w = w (x , y ).

(60)

Under these assumptions, the strains of the cracked ply are expressed as

εx =
εy =
εz =

γxy =

γyz =
γzx =

ν0
ν0
∂u
σ
= x0 − 230 σy − 210 σz = 0,
∂x
E2
E2
E2

(61)

σy
ν0
ν0
∂v
= − 230 σx + 0 − 210 σz = 0,
∂y
E2
E2
E2

(62)

σz
ν0
ν0
∂w
= − 210 σx − 210 σy + 0 = 0,
∂z
E2
E2
E1
σxy

0
G23

σyz

0
G12

=

=

(63)

∂v
∂u
+
= 0,
∂x
∂y

(64)

∂v
∂w
∂w
,
+
=
∂z
∂y
∂y

(65)

σzx
∂u
∂w
∂w
=
+
=
.
0
∂z
∂x
∂x
G12

(66)

From Eqs. (59)–(66), the stresses of the cracked ply are expressed asσx = σy = σz = σxy = 0,0σyz = G120σzx = G12

(67)

∂w ⎛
⎞
⎜x , y ⎟ ,
∂y ⎝
⎠

(68)

∂w ⎛
⎞
⎜x , y ⎟.
∂x ⎝
⎠

(69)

When the stresses expressed by Eqs. (67)–(69) are substituted into equilibrium Eqs. (32)–(34), the Laplace equation for displacement
w is obtained as

∂ 2w
∂ 2w
+
= 0.
2
∂x
∂y 2

(70)

To determine the displacement w, the boundary conditions of Eq. (70) are given by

w = 0 on y = 0,

(71)

∂w
= 0 on y = l,
∂y

(72)


(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe∂w= 0 on x = 0,∂xappw = γ12y(73)

on x = t ,

(74)

app
where γ12is the applied ply shear strain. At y = 0 , displacement w is not considered in Eq. (71). Eq. (72) is determined from Eq. (68),considering the shear stress condition σyz = 0 at the crack surface ( y = l ). Eq. (73) is deﬁned based on Eq. (69) considering the shear
stress condition σzx = 0 on the center plane ( x = 0 ) of the ply. In addition, in the displacement distribution at the interface ( x = t ),apppresented as Eq. (74), the neighboring ply is assumed to deform uniformly by γ12, regardless of the ply cracks. Separating variables
and assigning the boundary conditions presented above to Eq. (70), the displacement w (x , y ) of the Laplace equation can be expressed as(−1)n + 1 cosh[(2n − 1) πx /(2l)]8l2n − 1⎞⎛appw ⎜x , y ⎟ = 2 ∑sin ⎡πy⎤ γ12.
π n = 1 (2n − 1)2 cosh[(2n − 1) πt /(2l)]
2l
⎣
⎦
⎠
⎝
∞

(75)

When Eq. (75) is substituted into Eqs. (68) and (69), the local stress distribution under shear loading in the 1–2 plane in the RVE isderived as04G12(−1)n + 1 cosh[(2n − 1) πx /(2l)]2n − 1appcos ⎡πy⎤ γ12,∑π n = 1 2n − 1 cosh[(2n − 1) πt /(2l)]
⎣ 2l
⎦
∞

σyz =

0
4G12
(−1)n + 1 sinh[(2n − 1) πx /(2l)]2n − 1app
sin ⎡πy⎤ γ12.
∑π n = 1 2n − 1 cosh[(2n − 1) πt /(2l)] ⎣ 2l
⎦
∞

σzx =

(76)

(77)

Thus, the local stress distribution is formulated for a ply that includes ply cracking under shear loading. As indicated in Eq. (67) andEqs. (75)–(77), displacement w and the local stress distribution in a ply are constant with respect to the z-axis. With shear loading, thecompatibility condition for the strain is satisﬁed; therefore, the proposed stress distribution is an exact solution. Damage parameter ξ
can be formulated using the three-dimensional stress ﬁeld model for shear loading in the 1–2 plane. The average ply shear strain isdeﬁned as
ave
γ12
=

1
lth

∫0 ⎛⎜∫0 w ⎛x, l⎞ dx ⎞⎟ dz.
h

⎝

t

⎝
⎜

⎠

⎟

⎠

(78)

By substituting Eq. (75) into Eq. (78), the proportion of the average ply shear strain to the applied ply shear strain is derived as
ave
γ12

app
γ12

=

tanh[(2n − 1) πtk ρ /2]
16
1.
∑
π 3tk n = 1 (2n − 1)3
ρ
∞

(79)

Using Eq. (79), shear stiﬀness reduction Ds (ρ) associated with ply cracking subjected to shear loading can be deﬁned as [21]Ds (ρ) ≡ 1 −

ave
γ12

app
γ12

=1−

tanh[(2n − 1) πtk ρ /2]
16
1.
∑
π 3tk n = 1 (2n − 1)3
ρ
∞

(80)

Under the condition of low crack density (or the range in which damage parameter ω is small), we found that the damage value Ds in
Eq. (80) can be written with λ1 given by Eq. (38) as

Ds = 1 −

ave
γ12

app
γ12

≈

1
ω
λ1

(81)

Substituting Eq. (81) into Eq. (17), damage parameter ξ , which is a constant value, can be expressed using material parameter λ1.

ξ=

1
λ1

(82)

When considering the surface crack, damage parameter ξ is also written as Eq. (82). According to Li et al. [21], the constant shear
damage parameter ξ depends on the material of ply and laminate lay-up conﬁguration. Our derived shear damage parameter in Eq.
(82) depends on the material of ply only because our local stress ﬁeld model did not consider the stress transfer at the interface of thecracked ply.
When considering ply cracking, the eﬀective stiﬀness matrix C can be calculated using Eqs. (4), (57), and (82). The eﬀective
stiﬀness matrix C for surface cracking can be calculated using Eqs. (4), (58), and (82). Because damage parameter ω is formulated asfunctions of ply crack density using the local stress ﬁeld model, the stiﬀness matrix of the damaged ply is derived as a function of
crack density. In the next subsection, the stiﬀness matrix of damaged ply was employed in the three-dimensional laminate theory to
formulate the eﬀective compliance and the eﬀective thermal expansion coeﬃcient of laminates with arbitrary conﬁgurations as afunction of ply crack density.
(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

2.3. Eﬀective thermo-elastic properties of three-dimensional laminates with ply cracks
The eﬀective thermo-elastic properties of three-dimensional laminates with ply cracks were formulated with the help of the threedimensional laminate theory [19,22,23]. To determine the thermo-elastic properties of laminate utilizing the three-dimensionallaminate theory, two-stage coordinate conversion is applied to the eﬀective compliance and the thermal expansion coeﬃcient of thek-th ply (k = 1, 2, ⋯, N ;N is the total number of plies in the laminate) in the laminate. First, the compliance and the thermal
expansion coeﬃcient of the ply in principal axis O-123 is converted into those of the k-th ply in the coordinate system of laminate OXYZ. We assume that the X-Y plane is parallel to the 1–2 plane of the ply and the Z-axis is in the same direction as the 3-axis. The
direction of a ﬁber is tilted at an angle θ k between the 1- and X-axes, as indicated in Fig. 5.
The constitutive law of the k-th ply in the laminate coordinate system O-XYZ is described as

ε k = Skσ k + α k ∆T ,(83)

Sk = Rk (θ k ) C −1T k (−θ k ),(84)

where ∆T is the temperature change T-Tsf from the stress-free temperature Tsf to the testing temperature T. The eﬀective compliance
Sk and the thermal expansion coeﬃcient α k of the k-th ply with cracks in the laminate coordinate system O-XYZ are formulated as(see Eqs. (106) and (108)α k = Rk (θ k ) α 0,(85)where superscript k denotes the components of the k-th ply, and the thermal expansion coeﬃcient α 0 of a ply in the coordinate systemof the principal axis of ply O-123 is deﬁned as

α 0 = [ α10 α 20 α 20 0 0 0] .
T

(86)Here, α10 is the thermal expansion coeﬃcient of an undamaged ply along the 1-axis, and α 20 is that along the 2-axis. T k (θ k ) and Rk (θ k )are the coordinate conversion matrices of the stress and strain indicated in Eqs. (A.4) and (A.5). Sk is calculated by substituting Eqs.

(4), (A.4), and (A.5) into Eq. (84). Second, the stress, strain, and thermal expansion coeﬃcient of the k-th ply are converted into thosethat are divided into in-plane parts and out-of-plane parts. The components of stress σ k , strain ε k , and thermal expansion coeﬃcient
α k in the constitutive law of the k-th ply in Eq. (83) are deﬁned as
k
k
k ] ,
σ k = [ σXk σYk σZk σYZ
σZX
σXY
T

(87)

k
k
k
] ,
ε k = [ εXk εYk εZk γYZ
γZX
γXY
T

(88)

k ] .
α k = [ αXk αYk αZk 0 0 2αXY
T

(89)

The constitutive law of the stress σ k and strain ε k of the k-th ply that are divided into in-plane parts and out-of-plane parts can be
described aswhere

ε k = S kσ k + α k ∆T ,⎡ σI ⎤
σ k = ⎢ k ⎥,⎣ σO ⎦
k

(90)

⎡ε ⎤
ε k = ⎢ Ik ⎥,⎣ εO ⎦k⎡α ⎤α k = ⎢ Ik ⎥,⎣ αO ⎦
k(91)

Fig. 5. In-plane coordinate conversion to the laminate coordinate system O − XYZ .
(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabek⎡ σX ⎤k⎢σI = σYk ⎥,
⎢
⎥
k ⎥
⎢ σXY
⎣
⎦

⎡ σZ ⎤k ⎥,σOk = ⎢ σZX
⎢
⎥
k ⎥
⎢ σYZ
⎣
⎦
k

k
⎡ εX ⎤⎢kεI = ⎢ εYk ⎥⎥,⎢γk ⎥⎣ XY ⎦⎡ αX ⎤
,
αYk ⎥
⎥
k ⎥
⎢ 2αXY⎣⎦αIk = ⎢⎢⎡ εZ ⎤⎢ k ⎥εOk = ⎢ γZX⎥,⎢ γk ⎥⎣ YZ ⎦

k(92)⎡ αZ ⎤k ⎥.
αOk = ⎢ 2αZX⎢⎥k ⎥⎢ 2αYZ
⎣
⎦

k

k

(93)

Here, subscript I denotes the in-plane component, and O denotes the out-of-plane component. The relationship between Eqs. (87) and(88) and Eq. (91) can be formulated as

σ k = Aσ k ,(94)

ε k = Aε k ,(95)

where A is the conversion matrix in Eq. (B.3) (see Appendix B). Substituting Eqs. (94) and (95) into Eq. (83), the eﬀective complianceand thermal expansion coeﬃcient for the constitutive law between the strain and stress in Eq. (90) can be expressed using conversion
matrix A .
k
k
S IO
⎤
⎡ S II= ASk ATSk = ⎢ kk ⎥T⎣ (S IO ) SOO ⎦

(96)

⎡α ⎤
α k = ⎢ Ik ⎥ = Aα k
⎣ αO ⎦
k(97)

Finally, using the three-dimensional laminate theory (see Appendix C), the constitutive law of laminate with arbitrary lay-ups isformulated as
L
L
L
L
L
S IO
⎡ S II⎤ ⎡ σI ⎤ ⎡ αI ⎤⎡ε ⎤ε L = ⎢ IL ⎥ = S Lσ L + α L∆T = ⎢ L+ ⎢ L ⎥ ∆T ,L ⎥ ⎢ L⎥T
⎣ εO ⎦⎣ (S IO ) SOO ⎦ ⎣ σO ⎦ ⎣ αO ⎦(98)

where the stress and strain components are deﬁned asσI ⎤σL = ⎡⎢ σ L ⎥,⎣ O⎦
L

L
⎡ σX ⎤
L⎢
σI = σYL ⎥,
⎥
⎢
L
⎢
⎦
⎣ σXY ⎥L
⎡ σZ ⎤LL⎢
σO = σZX ⎥,
⎥
⎢
L
⎢
⎦
⎣ σYZ ⎥⎡ε ⎤ε L = ⎢ IL ⎥,⎣ εO ⎦LL

L⎡ εX ⎤
L⎢
εI = εYL ⎥,
⎥
⎢
L ⎥
⎢ γXY
⎦
⎣

L
⎡ εZ ⎤L ⎥⎢L
εO = ⎢ γZX ⎥,⎢γL ⎥⎣ YZ ⎦and the eﬀective compliancefollows.

SL⎡α ⎤α L = ⎢ IL ⎥,⎣ αO ⎦(99)⎡ αX ⎤
α L ⎥,⎢ Y ⎥L⎢⎦⎣ 2αXY ⎥αIL = ⎢

L

(100)

L
⎡ αZ ⎤LL ⎥,⎢αO = 2αZX⎥⎢
L⎢⎦⎣ 2αYZ ⎥(101)

and eﬀective thermal expansion coeﬃcient

L
L
S IO
L
⎡ S II⎤S L = [Sij ] = ⎢ L
,
L ⎥
T
(
S
)
S
IO
OO
⎣
⎦

αI ⎤
α L = [αiL] = ⎡⎢ α L ⎥,⎣ O⎦αLof laminate with arbitrary lay-ups can be obtained as

(102)

L

where

(103)

t
⎤
⎡
S IIL = ⎢ ∑ k (S IIk )−1⎥ ,t⎦⎣n=1 L
N

−1

(104)

t
⎡
L
k⎤
S IO
= S IIL ⎢ ∑ k (S IIk )−1S IO⎥,
t⎦⎣n=1 LN(105)


(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe
L
L T
L
SOO
= (S IO) (S IIL)−1S IO
+ ∑

tk k
k T
k
[SOO − (S IO) (S IIk )−1S IO
],
t
L
n=1
N

(106)

t
⎡
⎤
αIL = S IIL ⎢ ∑ k (S IIk )−1αIk⎥,
tL
n
=
1
⎣
⎦
N

(107)

L T
αOL = (S IO) (S IIL)−1αIL + ∑

tk k
k T
[αO − (S IO) (S IIk )−1αIk ],
t
n=1 L
N

where the 6 × 6 eﬀective compliance matrices S k and S L are divided into 3 × 3 submatrices(108)

Smk and SmL (m = II, IO, OO), as indicatedin Eqs. (96) and (102). Superscript L denotes the laminate components, and tL is laminate thickness. From the eﬀective complianceLS L = [Sij ] in Eq. (102) and the eﬀective thermal expansion coeﬃcient α L = [αiL] in Eq. (103) of laminate with ply cracks, the eﬀective
thermo-elastic constants of the damaged laminate can be calculated as

EXL =

1

L

S11

,

S
L
νXY
= − 12L ,
S11
L

L
GXY
=

1

L

S33

αXL = α1L,

L
αXY
= α3L,

,

EYL =

1

L

S22

EZL =

,

S
L
νXZ
= − 14L ,
S11
L

L
=
GXZ

αYL = α 2L,

1

L

S55

L
αXZ
= α5L,

,

1

,

L

S44

(109)

L
νYZ
=−

L
=
GYZ

αZL = α4L,

1

L

S66

L

S24
L

S22

,

(110)

,

(111)
(112)

L
αYZ
= α6L.

(113)

Using Eqs. (109)–(113), changes in the thermo-elastic properties of the laminate due to change in ply crack density are obtained
analytically by determining the thermo-elastic properties of a ply and the laminated constitution, with no ﬁtting parameter. The
proposed model for predicting the thermo-elastic properties of the laminate is an analytical model of ply cracking in general composite laminates and exhibits favorable compatibility with laminate theory because it handles only damaged plies. The present model
is not limited to symmetric laminates, including free surface; it formulates damage parameters ω and ξ analytically and therefore
incurs little computational cost.
2.4. Steady-state cracking analysis
This subsection considers the energy-based model for steady-state ply cracking. Steady-state cracking is the fracture mode that aply crack will propagate over the full-width of the specimen under constant thermomechanical loading. It is also assumed that a newply crack propagates between two pre-existing cracks only in the k-th ply under constant applied laminate tensile stress, as indicatedin Fig. 6, when the ply cracks are equally spaced. Ply crack density is deﬁned as ρ = 1/(2l) .
With these assumptions, the energy release rate Γk associated with ply cracking in a k-th ply under constant applied laminate
tensile stress σ L isFig. 6. Formation of a new crack between two pre-existing cracks in the k-th ply.
(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. OkabeΓk (ρ) = −U k (ρ /2) − 2U k (ρ)
,
tk

(114)

where U k (ρ) is strain energy in a laminate with length 2l and laminate thickness tL . The strain energy U k (ρ) is given by

U k (ρ) =tL L(σ − σ th ) TS L ⎛⎜σ L − σ th⎟⎞,
2ρ
⎝
⎠

(115)

is the laminate stress to cancel out the thermal residual stress of the k-th ply. When the thermomechanical loadingof thewhere
T
T T
laminate is equal to the stress σ th , all stress components of k-th ply are zero. The stress σ th can be written as σ th = [(σIth ) (σOth ) ] .
Under this loading condition, σIth and σOth are calculated using the ﬁrst row of Eqs. (C.1) and (C.2) and Eq. (C.3).

σ th

σL

k
L
σIth = (S IIL)−1 (S IIk σIk + (S IO− S IO) σOth − (αIL − αIk )∆T ),(116)

σOth = σOk.

(117)

Substituting σIk = σOk = 0 into Eqs. (116) and (117), the stress σ th is formulated as
L −1
Lk⎡ σI ⎤σ th = ⎢ th ⎥ = ⎡− (S II ) (αI − αI )∆T ⎤.
⎢⎥0⎦⎣ σO ⎦ ⎣th(118)

The energy release rate Γk can be calculated substituting Eqs. (115) and (118) into Eq. (114). In this study, the uniaxial monotonic
tensile stress σ L,app of the laminate is considered as

σ L = [ σ L,app 0 0 0 0 0]T(119)

Cracking analysis of the energy-based model for initiation of steady-state ply cracking in laminate is conducted as follows. It is
assumed that ply cracks form in just one ply. Using Eq. (114), the energy release rate Γk of each ply is calculated as a function of ply
crack density ρ at applied average laminate stress σ L,app . The critical ply crack density ρkc and the critical applied laminate stress σcL, kwhen the maximum value of the energy release rate max ρ (Γk) associated with ply cracking in a k-th ply is equal to the critical energy
release rate Γ c are calculated. The ply having minimum cracking laminate stress mink (σcL, k ) is determined as steady-state cracking
stress.
3. Results and discussion
The damage parameter ω has upper limit because the eﬀective stiﬀness C in Eq. (4) is formulated by assuming small damage.
Small damage refers to the condition in which the damage parameter ω is very small (i.e., ω2 ≪ 1). It also means relatively low crack
density because the damage parameter ω is a monotone increasing function of crack density. First, the upper limit of ω was investigated in the case of various CFRP and GFRP plies. Then, the proposed model was validated by calculating the eﬀective thermoelastic properties and steady-state cracking stress in comparison with experiment and ﬁnite-element analysis (FEA) results in previous
works on GFRP and CFRP laminates. The material properties of GFRP and CFRP unidirectional plies are listed in Table 1.
3.1. Maximum values of damage parameter and ply crack density
This subsection describes the upper limits of damage parameter ω and ply crack density. The eﬀective stiﬀness C is valid while theapproximate expression of shear stiﬀness reduction Ds in Eq. (81) is equal to that calculated using Eqs. (57) and (80). These upper
limits can be obtained by comparing shear stiﬀness reduction Ds - damage parameter ω curve calculated using Eqs. (57) and (80) withthat estimated using Eq. (81) for low damage conditions. Fig. 7 plots shear stiﬀness reduction Ds as a function of damage parameter ω
for CFRP unidirectional ply listed in Table 1 calculated using Eqs. (57) and (80) or Eq. (81).
In the range of small ω , shear stiﬀness reduction Ds estimated using Eq. (81) is in good agreement with that calculated using Eqs.
(57) and (80). A similar trend was observed for GFRP unidirectional plies. We determine that the damage parameter ω when therelative error of stiﬀness reduction Ds estimated using Eq. (81) to that calculated using Eqs. (57) and (80) is equal to 10 % is the upper
limit of damage parameter ω . The upper limit of the damage parameter is converted into the maximum value of normalized crack
Table 1
Material properties of GFRP and CFRP unidirectional plies.
E10
(GPa)

E20
(GPa)

0
ν12

0
ν23

0
G12
(GPa)

0
G23
(GPa)

α10
(/°C)

α20
(/°C)

Ply thickness

GFRP-1 [19]

41.7

13

0.3

0.42

3.4

4.58

46
44.7
142

13
12.7
9.85

0.3
0.297
0.3

0.42
0.3
0.46

5
5.8
4.48

4.6
4.885
3.37

× 10−6
–
–
–

0.203

GFRP-2 [24]
GFRP-3 [7]
CFRP [19]× 10−6
–
–
–

Type



6.72

29.3

(mm)

0.5
0.125
0.127

(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

Fig. 7. Shear stiﬀness reduction Ds as a function of damage parameter ω for CFRP ply.

density, which is the product of crack density and ply thickness using Eq. (57). The maximum values of damage parameter ωmax andnormalized crack density ρN max for CFRP and GFRP plies in Table 1 are listed in Table 2.
The maximum value of damage parameter ω is approximately 0.57. For laminates that are used to calculate stiﬀness reduction,
maximum crack densities ρmax ( ρN max divided by ply thickness) are listed in Table 3.
The stiﬀness reduction of laminates with ply cracking was investigated in the range of maximum ply crack density ρmax listed inTable 3.

3.2. Cross-ply laminate
This subsection discusses the elastic modulus and the crack opening displacement of a cross-ply laminate with damage due to ply
cracks or surface cracks in the 90° plies. Here [90/0]s GFRP laminate with surface cracks in the 90° plies was considered, and thematerial properties of GFRP-1 listed in Table 1 were used. Fig. 8(a) and (b) compare the results for the Young’s modulus and Poisson’s
ratio of the laminate as a function of ply crack density as determined with the present model and previously published FEA results[19].
The results of the present model agree well with the FEA results. The normalized elastic moduli as a function of ply crack densityof [0/90]s and [02/902]s CFRP laminates with ply cracks in the 90° plies are plotted in Figs. 9 and 10.
The experiment results [25] are also presented in Figs. 9 and 10. Here, EXL0 is the elastic modulus in the X-direction of theundamaged laminate. The results obtained from the present model are in good agreement with the experiment results. Fig. 11 illustrates the normalized Young’s modulus and normalized Poisson’s ratio as a function of the ply crack density of [0/908/01/2]s GFRP
laminate with ply cracks in the 90° plies as calculated by the present model, using the material properties of GFRP-3 listed in Table 1.
L0
is Poisson’s ratio in the X − Y plane of the undamaged laminate. The experiment results [7] are also presented in
Here, νXY
Fig. 11. The results predicted by the present model were in good agreement with the experiment results. For cross-ply laminate, ourmodel does not have an advantage over the simplistic shear-lag model. However, damage parameters derived in this paper aretheoretical solutions under the assumption that the ply next to the cracked ply is rigid, and therefore the present model exhibits
favorable compatibility with laminate theory that can be used to calculate the eﬀective thermo-elastic constants of laminates with
arbitrary lay-ups. Regarding this point, our model is superior to the shear-lag model that has diﬃculty solving the eﬀective stiﬀness ofdamaged laminates for arbitrary lay-ups.
The crack opening displacement for [0/90]s CFRP laminate is compared with the semi-analytical results of the National PhysicalLaboratory (NPL) model calculated by McCartney et al. [26]. Using Eq. (47), the crack opening displacement vCOD can be calculated
asTable 2
Maximum value of damage parameter ω and normalized crack density ρN .
Material

ωmax

ρN max

CFRP
GFRP-1GFRP-1 (surface crack)
GFRP-2GFRP-3

0.568
0.569
0.569
0.570
0.577

0.781
0.788
0.394
0.794
0.832



(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

Table 3
Maximum crack density of laminates used to calculate stiﬀness reduction.
Lay-upMaterial[90/0]s[0/90]s[02/902]s[55/ −55]N[67.5/ −67.5]N[0/90/ −45/ +45]s[0/908/01/2]s[0/ ±704 /01/2]s

GFRP-1 (surface crack)
CFRP
CFRP
GFRP-1GFRP-1GFRP-2GFRP-3
GFRP-3ρmax (/mm)
1.94
3.07
1.54
3.88
3.88
0.794
0.832
1.66

Fig. 8. (a) Young’s modulus and (b) Poisson’s ratio as a function of ply crack density for [90/0]s GFRP laminate with surface cracks.

Fig. 9. Normalized Young’s modulus as a function of ply crack density for [0/90]s CFRP laminate with ply cracks.

vCOD=2[lε2app − v (x , l, z )]cosh[(2n − 1) πλ1 x / (2l)] ⎤ app
8l
1⎡=2 ⎢1 − 2 ∑lε .
(2n − 1)2 cosh[(2n − 1) πλ1 t / (2l)] ⎥ 2
π
n=1
⎦
⎣
∞

(120)

The following material properties are used in the calculation of the crack opening displacement: E1 = 144.78 GPa, E2 = 9.58 GPa,
G12 = 4.785 GPa, ν12 = 0.31, ν23 = 0.55, and 2t = 0.25 mm. Crack opening displacement in 90° ply for CFRP [0/90]s laminate at crack
spacing 2l = 4.0 mm and applied laminate stress σ L,app = 0.2 GPa is plotted in Fig. 12.
The crack opening displacement calculated by the present model is less than that computed by the NPL model because our local
stress ﬁeld models or damage parameters did not consider stress transfer at the interface. It is assumed that the plies adjacent to thecracked ply are deformed uniformly in our model. Based on this assumption, the crack opening displacement in our model is less thanthat in the NPL model, which considers stress transfer at the interface. Our model takes into account stress transfer by reconsidering
the boundary conditions of Laplace equations for v and w in Eqs. (46) and (75). This improvement of our model is our main future
work.
There are many stiﬀness reduction models for cross-ply laminates. Lee et al. [14] derived stiﬀness reductions in cross-ply(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

Fig. 10. Normalized Young’s modulus as a function of ply crack density for [02/902]s CFRP laminate with ply cracks.

Fig. 11. (a) Normalized Young’s modulus and (b) normalized Poisson’s ratio as a function of ply crack density for [0/908/01/2]s GFRP laminate withply cracks.

Fig. 12. Crack opening displacement in 90° ply for CFRP [0/90]s laminate at crack spacing 2l = 4.0 mm and applied laminate stress σ L,app = 0.2 GPa.

laminates with ply cracks using an internal state variable approach. They formulated an internal state variable that describes thedamage state in composite materials using a stress ﬁeld model of cross-ply laminates including ply cracks. However, this stress ﬁeld
model does not satisfy the stress-free condition at the crack plane. Therefore, their model overestimated the stiﬀness of laminates. In
contrast, our local stress ﬁeld models fulﬁlled the stress-free condition on the crack plane. One- or two-dimensional shear-lag analyses[27] have been very useful for approximating the stress transfer at the interface by shear stress. However, shear-lag models must be
conﬁned to symmetric laminates and symmetric damage, because they neglect the bending eﬀect in their formulations. For example,(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

Nairn and Hu [28] formulated [90m/0n]s laminate including staggered ply cracks using variational mechanics analysis of the stresses.
For [90m/0n]s laminates, these shear-lag analyses were inadequate because antisymmetric or staggered ply cracks were observed inexperiments. Although other variational models [4,28–30] (e.g., the Nairn and Hu model) can treat laminates including staggered ply
cracks, these models are limited to simpliﬁed laminate conﬁgurations. In contrast, our model can handle the staggered ply crack
pattern and can consider laminates with any conﬁguration. Gudmundson and Zang [19] derived a general three-dimensional laminate model with ply cracks based on crack opening displacements, and their model can predict the average stress in each ply.
However, the local stress ﬁeld in a ply is not included in their model, whereas our model can analytically predict the local stress ﬁeld
model in 90° plies, as indicated in Eqs. (47)–(52).
3.3. Angle-ply laminate

This subsection discusses all the eﬀective thermoelastic properties of angle-ply GFRP laminates of [± 55]N and [± 67.5]N that
include ply cracks in each ply. The material properties of GFRP-1 listed in Table 1 were used in the calculations, and Figs. 13 and 14
present the Young’s moduli, shear moduli, Poisson’s ratios, and thermal expansion coeﬃcients for [± 55]N and [± 67.5]N laminates asfunctions of ply crack density estimated using the present model or FEA results [19].
For [± 55]N and [± 67.5]N laminates, the behavior predicted by the present model is approximately equal to that of the FEA results.
Fig. 15 indicates the normalized Young’s modulus and normalized Poisson’s ratio as a function of crack density of [0/ ± 704 /01/2]s GFRP
laminate with ply cracking in ± 70° plies as calculated by the present model, using the material properties of GFRP-3 listed in Table 1.
The experiment results by Varna et al. [7] are also presented in Fig. 15. Their experiment results are much smaller than the resultsof our model when crack density exceeds 0.2/mm. The stiﬀness of the experiment results declines drastically when crack density
increases from 0.2/mm to 0.4/mm. The formation of local delamination is cited as a possible cause of the diﬀerence between thepresent model and the experiment results. Experiment observation [31] indicates that matrix crack-induced delamination or edge
delamination occurs in the angle-ply laminates. Furthermore, as pointed out in [7,8], stiﬀness reduction for the [0/ ±θ4 /01/2]s GFRP
angle-ply laminate is due to ply cracking as well as shear-induced damage depending on the angle θ .
For angle-ply laminates, Vinogradov and Hashin [6] modeled [θm(1) / θn(2) ] laminates containing ply cracks in the middle ply using the

Fig. 13. (a) Young’s modulus, (b) shear modulus, (c) Poisson’s ratio, and (d) thermal expansion coeﬃcient as a function of ply crack density for[55/ −55]N GFRP laminate with ply cracks.
(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

Fig. 14. (a) Young’s modulus, (b) shear modulus, (c) Poisson’s ratio, and (d) thermal expansion coeﬃcient as a function of ply crack density for[67.5/ −67.5]N GFRP laminate with ply cracks.

Fig. 15. (a) Normalized Young’s modulus and (b) Normalized Poisson’s ratio as a function of ply crack density for [0/ ± 704 /01/2]s GFRP laminate.

principle of minimum complementary energy. Their model underestimated the stiﬀness of [0m/ θn]s angle-ply laminates in comparison
with experiment results [32], because the principle of minimum complementary energy provided a lower bound for the stiﬀness of acracked angle-ply. Gudmundson and Zang [19] and Lundmark and Varna [33–35] formulated the eﬀective properties of laminateswith ply cracks using crack opening displacement-based methods, which describe the elastic response changes caused by ply cracks ina medium by considering the crack opening displacements of individual cracks. In Gudmundson and Zang’s model, crack opening
displacement was calculated by ﬁtting the homogeneous isotropic medium with cracks under the action of uniform tractions on crack(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabesurfaces given by Benthem and Koiter [36] and Tada et al. [37]. In the model formulated by Lundmark and Varna, crack opening
displacement was computed by ﬁtting empirical results with ﬁnite-element calculations. In contrast, our stress ﬁeld model analytically derived the crack opening displacement in the transversely isotropic elastic body from the displacement.
3.4. Quasi-isotropic laminate

The quasi-isotropic [0/90/−45/+45]s GFRP laminate was analyzed to compare the results of our model with the experiment results
obtained by Tong et al. [38], using the material properties of GFRP-2 listed in Table 1. We assume that all plies except 0° plies havethe same damage, because the crack density in ± 45° plies is not available from the experiment of Tong et al. With this assumption,
the eﬀective stiﬀness calculated by the present model is the lower bound because cracks in ± 45° plies do not propagate perfectly
through the laminate width, and cracks in ± 45° initiate later than in 90° and contribute little to total stiﬀness reduction. Fig. 16 plots
the normalized Young’s modulus and normalized Poisson’s ratio for the laminate obtained from the proposed model in comparison
with the experiment results.
Although the contribution of ± 45° cracks to laminate stiﬀness reduction is small, the eﬀective Young’s modulus of laminate
obtained by our model is in good agreement with those of the experimental results. Fig. 16 is also plotted the result of present model
without ± 45° cracks. The present model without ± 45° cracks is slightly higher than the experimental results at the range of crack
density from 0/mm to 0.7/mm although the ± 45° cracks did not cause at this range [38]. This reason is possibly due to local
delamination at the crack tips.
A stiﬀness reduction model of quasi-isotropic laminates must be formulated, because most laminates used in real structures are
quasi-isotropic. However, few analytical models enable such predictions [10,11,17]. Tay and Lim [17] suggested a stiﬀness reduction
model of general laminates including ply cracks, using internal state variables that were ﬁtted by ﬁnite-element calculation.
McCartney [11] derived the stiﬀness reduction and progressive ply crack formation of general symmetric laminates with ply cracks
using a generalized plane strain analysis and the homogenization technique, and considering the stress transfer at the interfaces ofneighboring plies. However, this model could not be applied to asymmetric laminates, whereas the proposed model can predict thestiﬀness reduction of laminates with arbitrary conﬁgurations.
3.5. Steady-state cracking stress of cross-ply laminates
Steady-state cracking stress of cross-ply laminate estimated by the energy-based model described in Section 2.4 is compared withanalytical results in previous studies. Here, [0/90m/0] (m = 1, 2, 3, 4 ) cross-ply laminate of IM7/5250 is analyzed; properties of theply are listed in Table 4.
Fig. 17 plots steady-state cracking stress assuming steady-state cracking as a function of thickness of 90° plies per thickness of 0°plies.
For comparison, the analytical results using Gudmundson and Zang’s [19] model’s eﬀective compliance matrix and the energybased model (Eq. (114)) in this study are also plotted in Fig. 17. Furthermore, Fig. 17 includes the semi-analytical results of the NPL
model [39] and the Large Radius Axisymmetric Damage Model (LRAM) [39], as well as the analytical results of Dvorak and Laws[40]. The present model is in quantitatively good agreement with the analytical results of the LRAM and NPL models [26,41]. Bothsemi-analytical models consider the stress transfer of neighboring ply caused by ply cracking. The LRAM model was developed by
Schoeppner and Pagano [42] to approximate the thermoelastic stress ﬁeld and energy release rate in ﬂat laminates with ply cracks
and delaminations. When the radius-to-laminate thickness ratio is equal to 100,000, the large radius axisymmetric damage model isreasonably similar to the stress ﬁelds and stiﬀness reduction of ﬂat laminates estimated by the NPL model, which is considered to be
the most accurate generalized plane strain model. Our damage mechanics model is formulated explicitly and is highly simpliﬁed,

Fig. 16. (a) Normalized Young’s modulus and (b) normalized Poisson’s ratio as a function of ply crack density for [0/90/ −45/ +45]s GFRP laminate
with ply cracks in 90° and + 45° plies.
(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

Table 4
Material properties of IM7/5250-4 [39].
Longitudinal Young’s modulus E10

165.475 GPa0In-plane Poisson’s ratio ν12

0.31

Transverse Young’s modulus E20

10.342 GPa0Out-of-plane Poisson’s ratio ν23

0.56

0In-plane shear modulus G125.7922 GPa0Out-of-plane shear modulus G23

Longitudinal thermal expansion coeﬃcient α100.45 × 10−6 /°C

24.66 × 10−6 /°C
180 °C

Ambient testing temperature T
Ply thickness
Critical energy release Rate Γc

24 °C
0.127 mm3.3147 GPa

Transverse thermal expansion coeﬃcient α20
Stress-free temperature Tsf

225 J/m2

Fig. 17. Laminate crack initiation stress assuming steady-state cracking as a function of thickness of 90° plies per thickness of 0° plies.

whereas the LRAM and NPL models should solve the simultaneous diﬀerential equation numerically. The steady-state analyticalresults of Gudmundson and Zang’s model and those of Dvorak and Laws’s model are lower than the results of the present model, theLRAM, and the NPL model. Gudmundson and Zang’s model and Dvorak and Laws’s model assume that a cracked ply is identical to theinﬁnite cracked medium. This assumption does not work well, because cracked ply in laminate is signiﬁcantly constrained by adjacent plies. Therefore, the crack opening displacements of those models are higher than those of the other three models described
here. In this study, we validated the present steady-state cracking analysis model by comparing it to semi-analytical and analyticalresults for cross-ply laminates, because there is an abundance of models for cross-ply laminate. The present model can be used to
calculate crack density evolution for laminates with any lay-up. This is our future work, and we plan to publish the cracking analysisof laminate with any lay-up in our next paper.

4. Conclusions
In this study, the eﬀective compliance and elastic constants of laminates were formulated based on a CDM approach and laminate
theory, in an eﬀort to predict the thermoelastic properties of laminates of arbitrary conﬁgurations as a function of ply crack density.
The damage parameters ω and ξ were expressed as functions of ply crack density based on analytically formulated local stress ﬁeld
models subjected to tensile loading and shear loading. The model proposed in this study considers the eﬀect of damage due to ply
cracks (or surface cracks) based on only the thermomechanical properties of the ply and the laminate constitution. This model
quantitatively reproduced FEA and experiment results for the thermomechanical properties of cross-ply, angle-ply, and quasi-isotropic laminates including ply cracks. Following the approach of energy-based steady-state cracking analysis, the laminate crack
initiation stress for cross-ply is calculated and compared with analytical results in previous studies. Our model is quantitatively ingood agreement with the semi-analytical results of the large radius axisymmetric damage model by Pagano and the NPL model by
McCartney.

(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

Acknowledgements
This work was supported by 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). The
authors would like to acknowledge the vitally important encouragement and support made through the University of WashingtonTohoku University: Academic Open Space (UW-TU: AOS). S. O. appreciates the support of JSPS KAKENHI Grant No. JP 18J20899.
Appendix A. Coordinate conversion to the laminate coordinate system
The constitutive law of the principal axis of ply O-123 depicted in Fig. 1 is converted to that of the principal axis of laminate OXYZ in Fig. 5. The constitutive law of the ply with cracks in the principal axis O-123 is described as

ε app = C −1σ app + α 0∆T ,(A.1)

σ k = T k (θ k ) σ app,

(A.2)

where ∆T is the temperature change T-Tsf from the stress-free temperature Tsf to the testing temperature T. The relationship between
the stress σ app and strain ε app in the coordinate O-123 and the stress σ k and strain ε k in the coordinate O-XYZ can be expressed as

ε k = Rk (θ k ) ε app,

(A.3)

where T k (θ k ) and Rk (θ k ) are the coordinate conversion matrices of the stress and strain, described as
2 ksin2 θ k⎡ cos θ⎢ sin2 θ kcos2 θ k⎢00
T k (θ k ) = ⎢
0
0
⎢
⎢
0
0
⎢
k
k
k
k
⎣ sinθ cosθ − sinθ cosθ− sin2θ k ⎤
0
0
0
0
0
0
sin2θ k ⎥
⎥
1
0
0
0
⎥,
0 cosθ k sinθ k
0⎥⎥
0 − sinθ k cosθ k
0
⎥
0
0
0
cos2θ k ⎦000
− sinθ k cosθ k ⎤000sinθ k cosθ k ⎥
⎥
1
0
0
0
⎥.
0 cosθ k sinθ k
0
⎥
k
k
⎥
0 − sinθ cosθ0000
cos2θ k ⎥⎦2 ksin2 θ k⎡ cos θ⎢ sin2 θ k cos2 θ k⎢ 00
Rk (θ k ) = ⎢
0
⎢ 0
⎢ 0
0
⎢ sin2θ k − sin2θ k
⎣

(A.4)

(A.5)

The inverse matrices of T k (θ k ) and Rk (θ k ) are simply given by{T k (θ k )}−1 = T k (−θ k ),(A.6)

{Rk (θ k )}−1 = Rk ( −θ k ).

(A.7)

Substituting Eqs. (A.2), (A.3), and (A.6) into Eq. (A.1), the constitutive law of the k-th ply with ply cracks in the principal axis oflaminate O-XYZ is given by

ε k = Skσ k + α k ∆T ,(A.8)

where the eﬀective compliance matrix Sk and the thermal expansion coeﬃcient α k of the k-th ply in the coordinate system O-XYZ are

indicated in Eqs. (84) and (85).

Appendix B. Formulation of the conversion matrix A
The stress and strain of the k-th ply in Eqs. (87) and (88) are converted into thoses in Eq. (91) that are divided into in-plane partsand out-of-plane parts using the conversion matrix A , described as

σ k = Aσ k ,(B.1)

ε k = Aε k.

(B.2)

The conversion matrix A can be formulated as

⎞ ⎡1
⎞ ⎛
⎛
⎟ ⎢0
⎟ ⎜
⎜
0
A = P ⎜4, 6⎟ P ⎜3, 6⎟ = ⎢
⎢0
⎟ ⎢
⎟ ⎜
⎜
0
⎟ ⎢
⎟ ⎜
⎜
⎠ ⎣0
⎠ ⎝
⎝

0
1
0
0
0
0

0
0
0
1
0
0

0
0
0
0
0
1

0
0
0
0
1
0

0⎤
0⎥
1⎥
,
0⎥
0⎥
0⎥
⎦

(B.3)

where 6 × 6 matrix P (i, j ) = [Pkl] is the elementary matrix that is interchanged in two rows (or two columns) i and j. Pkl is deﬁned as(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

follows:⎧1 (k = l, k ≠ i, k ≠ j )⎪1 (k = i , l = j )Pkl =.
⎨1 (k = j , l = i)
⎪ 0 (other)
⎩

(B.4)

The nature of elementary matrix P (i, j ) isP −1 (i, j ) = P (i, j ),(B.5)

P T (i, j ) = P (i, j ).

(B.6)

Using Eqs. (B.5) and (B.6), the inverse matrix of A can be calculated as

A−1 = AT .

(B.7)

Substituting Eqs. (B.1), (B.2), and (B.7) into Eq. (A.8), the constitutive law between the stress and strain in Eq. (91) is expressed as

ε k = S kσ k + α k ∆T ,(B.8)

where the eﬀective compliance matrix S k and the thermal expansion coeﬃcient α k of the k-th ply for σ k and ε k are indicated in Eqs.
(96) and (97). The compliance matrix Sk of the orthotropic plate is described as

⎡ S11
⎢S k
⎢ 12
k
⎢ S13
Sk = ⎢
⎢ 0
⎢
⎢ 0
⎢S k
⎣ 16
k

k
k
S12
S13

k
S22
k
S23

k
S23
k
S33

0

0

0

0

k
k
S26
S36

0

0

0

0

0

0

k
S16
⎤
k⎥
S26 ⎥
k⎥
S36
⎥.
0 ⎥
⎥
0 ⎥
k⎥
S66
⎦

k
k
S44
S45

k
k
S45
S55

0

0

(B.9)

Substituting Eq. (B.9) into Eq. (96), the compliance matrix S k can be expressed as

⎡ S11
⎢S k
⎢ 12
k
⎢ S16
k
S =⎢ k
⎢ S13
⎢
⎢ 0
⎢ 0
⎣
k

k
k
k
S12
S16
S13

k
k
k
S22
S26
S23

0
0

k
k
k
S26
S66
S36

0

0

0

0

k
S55

0

0

0

k
k
k
S23
S36
S33

0
k
S45

0⎤
0⎥
⎥
0⎥
⎥.
0⎥
k⎥
S45
⎥
k ⎥
S44 ⎦

(B.10)

Appendix C. Formulation of the three-dimensional laminate theory
The three-dimensional laminate theory is utilized to formulate the thermoelastic properties of the composite laminate. The
constitutive law of the laminate for average stress and strain in Eq. (99) is expressed as
L
L
L
L
L
S IO
⎡ S II⎤ ⎡ σI ⎤ ⎡ αI ⎤⎡ε ⎤ε L = ⎢ IL ⎥ = S Lσ L + α L∆T = ⎢ L+
L ⎥ ⎢ L⎥⎢ L ⎥ ∆T ,T
⎣ εO ⎦⎣ (S IO ) SOO ⎦ ⎣ σO ⎦ ⎣ αO ⎦(C.1)

and the constitutive law of the k-th ply for average stress and strain in Eq. (91) is given by
k
k
k
k
k
S IO
⎡ S II⎤ ⎡ σI ⎤ ⎡ αI ⎤
⎡ε ⎤
∆T ,
+
ε k = ⎢ Ik ⎥ = S kσ k + α k ∆T = ⎢ k
⎥
⎢
⎥
⎢
k
k⎥
k
T
⎣ εO ⎦⎣ (S IO ) SOO ⎦ ⎣ σO ⎦ ⎣ αO ⎦(C.2)

where superscript L denotes the laminate component and k denotes the k-th ply component. From the compatibility and equilibrium
conditions in the laminate, the following relationships must be satisﬁed.

εIk = εIL,σOk = σOL(C.3)

The laminate average stresses and strains are deﬁned as

σIL ⎤
tk k
tk ⎡ σIk ⎤σL = ⎡⎢ σ L ⎥ = ∑ tL σ = ∑ tL ⎢ σ k ⎥,k=1
⎣ O ⎦ k=1
⎣ O⎦
N

N

(C.4)



(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe

where

N
N
L
t
t ⎡ε k ⎤⎡ε ⎤ε L = ⎢ IL ⎥ = ∑ k ε k = ∑ k ⎢ Ik ⎥,t
k = 1 L ⎣ εO ⎦⎣ εO ⎦ k = 1 tL

(C.5)

tL = ∑ tk.
N

k=1

(C.6)

Here, tL is laminate thickness, and tk is k-th ply thickness. From the ﬁrst row of Eqs. (C.2) and (C.3), in-plane average stress of k-th ply
σIk can be obtained ask
σIk = (S IIk )−1 (εIL − S IO
σOL − αIk ∆T ).

(C.7)

Substituting Eq. (C.7) into the ﬁrst row of Eq. (C.4), the in-plane average laminate strain εIL is expressed as
t
t
t
⎤ ⎡
⎡
⎤
⎡
k ⎤ L
εIL=⎢ ∑ tk (S IIk )−1⎥ σIL + ⎢ ∑ tk (S IIk )−1⎥ ⎢ ∑ tk (S IIk )−1S IO
⎥ σO
L
L
L
⎦
⎦ ⎣ k=1
⎣ k=1
⎦
⎣ k=1
−1

N

N

−1

N

t
t
⎡
⎤ ⎡
⎤
+ ⎢ ∑ tk (S IIk )−1⎥ ⎢ ∑ tk (S IIk )−1αIk ⎥ ∆T .
L
L
⎣ k=1
⎦ ⎣ k=1
⎦
N

−1

N

(C.8)

L
By comparing Eq. (C.8) and the ﬁrst row of Eq. (C.1), the compliance submatrices of the laminate S IIL and S IO
and in-plane thermal
expansion coeﬃcient of laminate αIL are obtained as Eqs. (104), (105), and (107). In terms of the out-of-plane components, from thesecond row of Eqs. (C.2) and (C.3), the out-of-plane average strain of k-th ply εOL is formulated ask T kkεOk = (S IO) σI + SOO
σOL + αOk ∆T .

(C.9)

Inserting Eq. (C.9) into the second row of Eq. (C.5) and using Eq. (C.7) and the ﬁrst row of Eq. (C.1), the out-of-plane average strain oflaminate can be expressed ask TεOL=∑ tk [(S IO) (S IIk )−1S IIL σIL
N

k=1

t

L

k T
L
k
k
+ {(S IO
) (S IIk )−1 (S IO− SIO) + SOO} σOL
k T+ {(S IO
) (S IIk )−1 (αIL − αIk ) + αOk }∆T ].

(C.10)

The following relationship can be obtained using Eq. (105).

∑ tk (SIOk )T (SIIk )−1 = (SIOL )T (SIIL)−1
N

k=1

t

(C.11)L

Substituting Eq. (C.11) into Eq. (C.10), the out-of-plane average strain of laminate is rewritten as
L T L
εOL=(S IO) σIt⎡ L T L −1 L
k
k T
k ⎤ L
− (S IO+ ⎢ (S IO) (S II ) S IO + ∑ tk (SOO
) (S IIk )−1S IO) ⎥ σO
L
1
k
=
⎦
⎣
N

t
⎤
⎡ L T L −1 Lk T
+ ⎢ (S IO) (S II ) αI + ∑ tk (αOk − (S IO) (S IIk )−1αIk ) ⎥ ∆T .
L
k
1
=
⎦
⎣
N

(C.12)L
By comparing Eq. (C.12) and the second row of Eq. (C.1), the compliance submatrix of the laminate SOO
and out-of-plane thermalL
expansion coeﬃcient of laminate αO are obtained as Eqs. (106) and (108).

References[1] Nairn J, Hu S. Micromechanics of damage: a case study of matrix microcracking. In: Talreja R, editor. Damage mechanics of composite materials. Amsterdam:
Elsevier; 1994. p. 187–243.
[2] Takeda N, Ogihara S. In situ observation and probabilistic prediction of microscopic failure processes in CFRP cross-ply laminates. Compos Sci Technol
1994;52:183–95. https://doi.org/10.1016/0266-3538(94)90204-6.
[3] Talreja R, Singh C. Damage and failure of composite materials. New York: Chambridge University Press; 2012. https://doi.org/10.1017/CBO9781139016063.
[4] Hashin Z. Analysis of cracked laminates: a variational approach. Mech Mater 1985;4:121–36. https://doi.org/10.1016/0167-6636(85)90011-0.
[5] McCartney L. Theory of stress transfer in a 0° -90° -0° cross-ply laminate containing a parallel array of transverse cracks. J Mech Phys Solids 1992;40:27–68.
https://doi.org/10.1016/0022-5096(92)90226-R.
[6] Vinogradov V, Hashin Z. Variational analysis of cracked angle-ply laminates. Compos Sci Technol 2010;70:638–46. https://doi.org/10.1016/j.compscitech.
2009.12.018.
[7] Varna J, Joﬀe R, Akshantala N, Talreja R. Damage in composite laminates with oﬀ-axis plies. Compos Sci Technol 1999;59:2139–47. https://doi.org/10.1016/
S0266-3538(99)00070-6.

(QJLQHHULQJ)UDFWXUH0HFKDQLFV  

S. Onodera and T. Okabe[8] Singh CV, Talreja R. Analysis of multiple oﬀ-axis ply cracks in composite laminates. Int J Solids Struct 2008;45:4574–89. https://doi.org/10.1016/j.ijsolstr.2008.
04.004.
[9] Hajikazemi M, McCartney L, Van Paepegem W, Sadr M. Theory of variational stress transfer in general symmetric composite laminates containing non-uniformly
spaced ply cracks. Compos Part A: Appl Sci Manuf 2018;107:374–86. https://doi.org/10.1016/j.compositesa.2018.01.021.
[10] Hajikazemi M, McCartney L. Comparison of Variational and Generalized Plane Strain approaches for matrix cracking in general symmetric laminates. Int J
Damage Mech 2018;27:507–40. https://doi.org/10.1177/1056789516685381.
[11] McCartney L. Model to predict eﬀects of triaxial loading on ply cracking in general symmetric laminates. Compos Sci Technol 2000;60:2255–79. https://doi.org/
10.1016/S0266-3538(00)00086-5.
[12] Allen D, Harris C, Groves S. A thermomechanical constitutive theory for elastic composites with distributed damage-I. Theoretical development. Int J Solids
Struct 1987;23:1301–18. https://doi.org/10.1016/0020-7683(87)90107-7.
[13] Murakami S. Continuum damage mechanics. Dordrecht: Springer; 2012. https://doi.org/10.1007/978-94-007-2666-6.
[14] Lee J-W, Allen D, Harris C. Internal state variable approach for predicting stiﬀness reductions in ﬁbrous laminated composites with matrix cracks. J Compos
Mater 1989;23:1273–91. https://doi.org/10.1177/002199838902301205.
[15] Okabe T, Onodera S, Kumagai Y, Nagumo Y. Continuum damage mechanics modeling of composite laminates including transverse cracks. Int J Damage Mech
2018;27:877–95. https://doi.org/10.1177/1056789517711238.
[16] Talreja R. A continuum mechanics characterization of damage in composite materials. Proc Roy Soc Lond, Ser A: Math Phys Sci 1985;399:195–216. https://doi.
org/10.1098/rspa.1985.0055.
[17] Tay T, Lim E. Analysis of composite laminates with transverse cracks. Compos Struct 1996;34:419–26. https://doi.org/10.1016/0263-8223(96)00010-4.
[18] Kachanov L. Time of the rupture process under creep conditions. Izvestiia Akademii Nauk SSSR, Otdelenie Teckhnicheskikh Nauk 1958;8:26–31.
[19] Gudmundson P, Zang W. An analytic model for thermoelastic properties of composite laminates containing transverse matrix cracks. Int J Solids Struct
1993;30:3211–31. https://doi.org/10.1016/0020-7683(93)90110-S.
[20] Li S, Wang M, Jeanmeure L, Sitnikova E, Yu F, Pan Q, et al. Damage related material constants in continuum damage mechanics for unidirectional composites
with matrix cracks. Int J Damage Mech 2019;28:690–707. https://doi.org/10.1177/1056789518783239.
[21] Li S, Reid SR, Soden PD. A continuum damage model for transverse matrix cracking in laminated ﬁbre-reinforced composites. Philos Trans Roy Soc Lond Ser A:
Math Phys Eng Sci 1998;356:2379–412. https://doi.org/10.1098/rsta.1998.0278.
[22] Pagano NJ. Exact moduli of anisotropic laminates. In: Reddy J, editor. Mech compos mater Dordrecht: Springer; 1994. p. 210–31. https://doi.org/10.1007/97894-017-2233-9_18.
[23] Sun C, Li S. Three-dimensional eﬀective elastic constants for thick laminates. J Compos Mater 1988;22:629–39. https://doi.org/10.1177/002199838802200703.
[24] Tong J, Guild FJ, Ogin SL, Smith PA. On matrix crack growth in quasi-isotropic laminates – II. Finite element analysis. Compos Sci Technol 1997. https://doi.
org/10.1016/S0266-3538(97)00083-3.
[25] Groves S, Harris C, Highsmith A, Allen D, Norvell R. An experimental and analytical treatment of matrix cracking in cross-ply laminates. Exp Mech
1987;27:73–9. https://doi.org/10.1007/BF02318867.
[26] McCartney LN, Schoeppner GA, Becker W. Comparison of models for transverse ply cracks in composite laminates. Compos Sci Technol 2000;60:2347–59.
https://doi.org/10.1016/S0266-3538(00)00030-0.
[27] Berthelot J-M. Transverse cracking and delamination in cross-ply glass-ﬁber and carbon-ﬁber reinforced plastic laminates: Static and fatigue loading. Appl Mech
Rev 2003;56:111–47. https://doi.org/10.1115/1.1519557.
[28] Nairn JA, Hu S. The formation and eﬀect of outer-ply microcracks in cross-ply laminates: a variational approach. Eng Fract Mech 1992;41:203–21. https://doi.
org/10.1016/0013-7944(92)90181-D.
[29] Li S, Reid SR, Soden PD. A ﬁnite strip analysis of cracked laminates. Mech Mater 1994;18:289–311. https://doi.org/10.1016/0167-6636(94)90041-8.
[30] Li S, Hafeez F. Variation-based cracked laminate analysis revisited and fundamentally extended. Int J Solids Struct 2009;46:3505–15. https://doi.org/10.1016/j.
ijsolstr.2009.03.031.
[31] Johnson P, Chang FK. Characterization of matrix crack-induced laminate failure – Part I: Experiments. J Compos Mater 2001;35:2009–35. https://doi.org/10.
1106/7RN1-PFBN-XQR9-3KDK.
[32] Fikry M, Kitamura R, Ogihara S. Eﬀect of matrix cracking on mechanical properties in FRP angle-ply laminates. M&M 2017 Zairyorikigaku conference The Japan
Society of Mechanical Engineers Materials & Mechanics Division; 2017. p. 818–21. https://doi.org/10.1299/jsmemm.2017.OS1014.
[33] Lundmark P, Varna J. Constitutive relationships for laminates with ply cracks in in-plane loading. Int J Damage Mech 2005;14:235–59. https://doi.org/10.1177/
1056789505050355.
[34] Lundmark P, Varna J. Crack face sliding eﬀect on stiﬀness of laminates with ply cracks. Compos Sci Technol 2006;66:1444–54. https://doi.org/10.1016/j.
compscitech.2005.08.016.
[35] Varna J. Physical interpretation of parameters in synergistic continuum damage mechanics model for laminates. Compos Sci Technol 2008;68:2592–600.
https://doi.org/10.1016/j.compscitech.2008.04.037.
[36] Benthem J, Koiter WT. Asymptotic approximations to crack problems. In: Sih G, editor. Mechanics of fracture 1: methods of analysis and solutions of crack
problems The Netherlands: Leyden; 1973. p. 131–78. https://doi.org/10.1007/978-94-017-2260-5_3.
[37] Tada H, Paris P, Irwin G. The stress analysis of cracks handbook. 3rd ed. New York: The American Society of Mechanical Engineers; 2000. https://doi.org/10.
1115/1.801535.
[38] Tong J, Guild FJ, Ogin SL, Smith PA. On matrix crack growth in quasi-isotropic laminates-I. Experimental investigation. Compos Sci Technol 1997. https://doi.
org/10.1016/S0266-3538(97)00080-8.
[39] Pagano N, Schoeppner G, Kim R, Abrams F. Steady-state cracking and edge eﬀects in thermo-mechanical transverse cracking of cross-ply laminates. Compos Sci
Technol 1998;58:1811–25. https://doi.org/10.1016/S0266-3538(98)00047-5.
[40] Dvorak G, Laws N. Analysis of ﬁrst ply failure in composite laminates. Eng Fract Mech 1986;25:763–70. https://doi.org/10.1016/0013-7944(86)90039-1.
[41] McCartney LN, Schoeppner GA. Predicting the eﬀect of non-uniform ply cracking on the thermoelastic properties of cross-ply laminates. Compos Sci Technol
2002;62:1841–56. https://doi.org/10.1016/S0266-3538(02)00091-X.
[42] Schoeppner GA, Pagano NJ. Stress ﬁelds and energy release rates in cross-ply laminates. Int J Solids Struct 1998;35:1025–55. https://doi.org/10.1016/S00207683(97)00107-8.

