# Investigation of Through-Thickness Stresses in Composite Laminates Using Layerwise Theory

Authors: Hamidreza Yazdani Sarvestani; Ali Naghashpour
Journal: International Journal of Engineering Mathematics
 (2013)

---

## Abstract

In this study, an analytical method is developed to exactly obtain the interlaminar stresses near the free edges of laminated composite plates under the bending moment based on the reduced form of elasticity displacement field for a long laminate. The analytical and numerical studies were performed based on the Reddy’s layerwise theory for the boundary layer stresses within cross-ply, symmetric, angle-ply, and general composite laminates. Finally, a variety of numerical results are presented for the interlaminar normal and shear stresses along the interfaces and through thickness of laminates near the free edges. The results showed high stress gradient of interlaminar normal and shear stresses near the edges of laminates.


---

## Body

## 1. Introduction

Due to high specific strength and stiffness of fiber reinforced polymer composite, laminated composites have found augmented use in many industrial applications. In the boundary-layer regions, due to geometry and material discontinuities, the interlaminar stresses exhibit much higher values than those predicted by the classical lamination theory (CLT). These highly concentrated stresses cause delaminate failure in the laminates. However, no exact solution is found for elasticity equations because of inherent complexities involved in the problem of finding exact stress values in the edges. Hence, different analytical and numerical methods for finding the interlaminar stresses are used to describe the interlaminar stresses at the free edges of composite laminates. Complete literature surveys on this subject are available in review articles of Kant and Swaminathan [1] which obviously show the detailed path of development of methods. The first approximate analysis of interlaminar stresses was presented by Puppo and Evensen [2]. They studied interlaminar shear stresses in an idealized laminate consisting of orthotropic layers separated by isotropic shear layers with interlaminar normal stress being neglected through the laminate. Other approximate analytical methods utilized to examine the problem consist of the use of the higher-order plate theory by Pagano [3], the perturbation technique by Hsu and Herakovich [4], the boundary layer theory by Tang and Levy [5], and the approximate elasticity solutions by Pipes and Pagano [6]. An approximate theory is also employed by Pagano [7] based on the assumed inplane stresses and the use of Reissner’s variational principle. Wang and Choi [8] utilized Lekhnistskii’s stress potential and the theory of anisotropic elasticity for examining the free edge singularities. A variational approach concerning Lekhnitskii’s stress functions is used by Yin [9] for the evaluation of free-edge stresses in laminates under uniaxial tension, bending, and torsion. Wang and Crossman [10] developed a quasi-three-dimensional finite element solution to determine the free-edge stresses in a symmetric balanced composite laminate under uniaxial tension and uniform thermal loading. Whitcomb et al. [11] studied the differences in numerical results for interlaminar stresses obtained by various methods (finite difference methods, finite element methods, and perturbation techniques). Boundary element method and the integral equation theory were used by Davì [12] to study the stresses in a general laminate under uniform axial strain. Carrera and Demasi [13, 14] studied the accuracy of the finite-element mixed layerwise solutions by using the Reissner mixed variational theorem (RMVT). They compared the numerical results for interlaminar stresses in several finite-element models and elasticity theory within composite laminates and sandwich plates. Nguyen and Caron [15] employed a multiparticle finite element method to study the interlaminar stresses near the free edges of general composite laminates under mechanical and thermal loading. Robbins and Reddy [16] used a displacement-based variable kinematic global-local finite element method. Mittelstedt and Becker [17] utilized Reddy’s layerwise laminate plate theory to find the closed-form analysis of free-edge effects in layered plates of arbitrary nonorthotropic layups. The approach consists of the subdivision of the physical laminate layers into an arbitrary number of mathematical layers through the plate thickness. Na [18] used a finite element model based on the layerwise theory. He employed the von Kármán type nonlinear strains to analyze damage in laminated composite beams. In his formulation, the Heaviside step function is employed to express the discontinuous interlaminar displacement field at the delaminated interfaces. Plagianakos and Saravanos [19] presented a higher-order layerwise theoretical framework, which enables prediction of the static response of thick composite and sandwich composite plates. Ullah et al. [20] carried out some experimental tests to characterize the behavior of a woven CFRP material under large-deflection bending. Two-dimensional finite element (FE) models were implemented in the commercial code Abaqus. They performed series of simulations to study the deformation behavior and damage in CFRP for cases of high-deflection bending. Hélénon et al. [21] presented an experimental and numerical investigation into failure of T-shaped laminated composite structures. Three out-of-plane bending cases are studied. They found that very high free-edge maximum principal transverse tensile stresses perpendicular to the fiber direction occur at the failure locations. Thai et al. [22] indicated an isogeometric finite element formulation for static, free vibration and buckling analysis of laminated composite and sandwich plates. Their method allows removing shear correction factors and improves the accuracy of transverse shear stresses. Thai et al. [23] investigated the behavior of laminated composites using several high order or layerwise finite element calculations. A layerwise model and its dedicated C° eight-node finite element were specifically developed for interlaminar stresses analysis in free edge problem. Malekzadeh [24] developed a high accuracy and rapid convergence hybrid approach for two-dimensional static analyses of circular arches with different boundary conditions. The method essentially consists of a layerwise theory used for the thickness direction in conjunction with differential quadrature method in the axial direction.There have been very limited works to study the interlaminar stresses subjected to the bending moment. In the present paper, by the use of Reddy’s LWT, an analytical solution is presented to evaluate interlaminar stresses in cross-ply, symmetric, angle-ply, and general composite laminates under the bending moment. To commence with, based on physical arguments regarding the deformations of a long generally and other laminated composite plates, an appropriate reduced elasticity displacement field is established. The boundary-layer stresses within the laminate are obtained analytically based on Reddy’s LWT.


## 2. Problem Formulation

### 2.1. Elasticity Displacement Field

AnNth-layered composite plate (with arbitrary lamination) under the bending moment is considered as shown in Figure 1. The coordinate system (x,y,z) is located at the middle plane of the laminate, that is of thickness h, width 2b, and length 2a and is assumed to be long in the x direction so that the strains away from the ends (x=±a) of the laminate are functions of only y and z.Figure 1
Laminate geometry and coordinate system.The integrations of the three-dimensional elasticity strain-displacement relations [25] inside the kth layer of the laminate will generate the most general form of displacement field which can be shown to be [26]
(1)u1(k)(x,y,z)=B4(k)xy+B6(k)xz+B2(k)x+u(k)(y,z)u2(k)(x,y,z)=-B1(k)xz+B3(k)x-12B4(k)x2+v(k)(y,z)u3(k)(x,y,z)=B1(k)xy+B5(k)x-12B6(k)x2+w(k)(y,z),
where u1(k), u2(k), and u3(k) represent the displacement components of the material point (x,y,z) in the x, y, and z directions, respectively. To satisfy the displacement continuity conditions at the interfaces of the adjoining layers, the integration constants in (1) must be the same for all layers within the laminate. Therefore, the most general displacement field in the kth layer is shown:
(2)u1(k)(x,y,z)=B4xy+B6xz+B2x+u(k)(y,z),u2(k)(x,y,z)=-B1xz+B3x-12B4x2+v(k)(y,z),u3(k)(x,y,z)=B1xy+B5x-12B6x2+w(k)(y,z).
If loading conditions at x=-a and a are similar, based on the physical grounds, the following conditions hold:
(3)u1(k)(x,y,z)=-u1(k)(-x,-y,z),u2(k)(x,y,z)=-u2(k)(-x,-y,z),u3(k)(x,y,z)=u3(k)(-x,-y,z).
Upon imposing these restrictions on (1) it is readily seen that the constants B4 and B5 must vanish and the displacement field in (1) is reduced to the following equations:(4a)u1(k)(x,y,z)=B2x+B6xz+u(k)(y,z),(4b)u2(k)(x,y,z)=-B1xz+B3x+v(k)(y,z),(4c)u3(k)(x,y,z)=B1xy-12B6x2+w(k)(y,z).Also, by replacing u(k)(y,z) by -B3y+u(k)(y,z) in (4a), it becomes apparent that terms involving B3 in (4a)–(4c) can be neglected since no strains are generated by such terms (these terms will correspond to an infinitesimal rigid-body rotation of the laminate about the z-axis in Figure 1). Thus, the most general form of the displacement field of an arbitrary laminate is given by(5a)u1(k)(x,y,z)=B2x+B6xz+u(k)(y,z),(5b)u2(k)(x,y,z)=-B1xz+v(k)(y,z),(5c)u3(k)(x,y,z)=B1xy-12B6x2+w(k)(y,z).The displacement field in (5a), (5b), and (5c) may be used, in principle, for computing the stress field in any composite laminate subjected to arbitrary combinations of self-equilibrating mechanical and uniform hygrothermal loads. The terms involving B1, B2, and B6 in (5a), (5b), and (5c) demonstrate certain global deformations that occurred in the laminate. On the other hand, the unknown functions appearing in (5a), (5b), and (5c) illustrate the local deformations. They will be evaluated by LWT of Reddy that happened in a laminate.For cross-ply laminates subjected to mechanical loadings, based on physical grounds, the following restrictions hold (see Figure1):(6a)u1(k)(x,y,z)=-u1(k)(-x,y,z),(6b)u2(k)(x,y,z)=u2(k)(-x,y,z).Upon imposing these conditions on (5a), it is readily seen that u(k)(y,z)=0. Also it is concluded from (6b) and (5b) that B1=0. Thus, for cross-ply laminates the most general form of the displacement field is given:
(7)u1(k)(x,y,z)=B6xz+B2x,u2(k)(x,y,z)=v(k)(y,z),u3(k)(x,y,z)=-12B6x2+w(k)(y,z).Also, for symmetric laminates based on physical grounds the following conditions must hold (see Figure1): (8a)u1(k)(x,y,z)=u1(N+1-k)(x,y,-z)(8b)u2(k)(x,y,z)=u2(N+1-k)(x,y,-z)(8c)u3(k)(x,y,z)=-u3(N+1-k)(x,y,-z)(k=1,2,3,…,N2).Upon imposing these conditions on (5a), (5b), and (5c), it is found that B1=B6=0. Thus, for symmetric laminates the most general form of displacement field is expressed by
(9)u1(k)(x,y,z)=B2x+u(k)(y,z),u2(k)(x,y,z)=v(k)(y,z),u3(k)(x,y,z)=w(k)(y,z).Finally, for antisymmetric angle-ply laminates withN layers subjected to mechanical loadings, the following condition must hold:
(10)u1(k)(x,y,z)=u1(N+1-k)(x,-y,-z)(k=1,2,3,…,N2).
Concluding from (10) and (5a), B6=0 and, therefore, the most general displacement field for such laminates are written as
(11)u1(k)(x,y,z)=B2x+u(k)(y,z),u2(k)(x,y,z)=-B1xz+v(k)(y,z),u3(k)(x,y,z)=B1xy+w(k)(y,z).The displacement fields in (7), (9), and (11) can be represented in one place as
(12)u1(k)(x,y,z)=B2x+δCB6xz+δSu(k)(y,z),u2(k)(x,y,z)=-δAB1xz+v(k)(y,z),u3(k)(x,y,z)=δAB1xy-12δCB6x2+w(k)(y,z),
where δA=δS=0 and δC=1 for general cross-ply laminates, δA=δC=0 and δS=1 for symmetric laminates, δA=δS=1 and δC=0 for angle-ply laminates, and δA=δS=δc=1 for general laminates.


## 2.1. Elasticity Displacement Field

AnNth-layered composite plate (with arbitrary lamination) under the bending moment is considered as shown in Figure 1. The coordinate system (x,y,z) is located at the middle plane of the laminate, that is of thickness h, width 2b, and length 2a and is assumed to be long in the x direction so that the strains away from the ends (x=±a) of the laminate are functions of only y and z.Figure 1
Laminate geometry and coordinate system.The integrations of the three-dimensional elasticity strain-displacement relations [25] inside the kth layer of the laminate will generate the most general form of displacement field which can be shown to be [26]
(1)u1(k)(x,y,z)=B4(k)xy+B6(k)xz+B2(k)x+u(k)(y,z)u2(k)(x,y,z)=-B1(k)xz+B3(k)x-12B4(k)x2+v(k)(y,z)u3(k)(x,y,z)=B1(k)xy+B5(k)x-12B6(k)x2+w(k)(y,z),
where u1(k), u2(k), and u3(k) represent the displacement components of the material point (x,y,z) in the x, y, and z directions, respectively. To satisfy the displacement continuity conditions at the interfaces of the adjoining layers, the integration constants in (1) must be the same for all layers within the laminate. Therefore, the most general displacement field in the kth layer is shown:
(2)u1(k)(x,y,z)=B4xy+B6xz+B2x+u(k)(y,z),u2(k)(x,y,z)=-B1xz+B3x-12B4x2+v(k)(y,z),u3(k)(x,y,z)=B1xy+B5x-12B6x2+w(k)(y,z).
If loading conditions at x=-a and a are similar, based on the physical grounds, the following conditions hold:
(3)u1(k)(x,y,z)=-u1(k)(-x,-y,z),u2(k)(x,y,z)=-u2(k)(-x,-y,z),u3(k)(x,y,z)=u3(k)(-x,-y,z).
Upon imposing these restrictions on (1) it is readily seen that the constants B4 and B5 must vanish and the displacement field in (1) is reduced to the following equations:(4a)u1(k)(x,y,z)=B2x+B6xz+u(k)(y,z),(4b)u2(k)(x,y,z)=-B1xz+B3x+v(k)(y,z),(4c)u3(k)(x,y,z)=B1xy-12B6x2+w(k)(y,z).Also, by replacing u(k)(y,z) by -B3y+u(k)(y,z) in (4a), it becomes apparent that terms involving B3 in (4a)–(4c) can be neglected since no strains are generated by such terms (these terms will correspond to an infinitesimal rigid-body rotation of the laminate about the z-axis in Figure 1). Thus, the most general form of the displacement field of an arbitrary laminate is given by(5a)u1(k)(x,y,z)=B2x+B6xz+u(k)(y,z),(5b)u2(k)(x,y,z)=-B1xz+v(k)(y,z),(5c)u3(k)(x,y,z)=B1xy-12B6x2+w(k)(y,z).The displacement field in (5a), (5b), and (5c) may be used, in principle, for computing the stress field in any composite laminate subjected to arbitrary combinations of self-equilibrating mechanical and uniform hygrothermal loads. The terms involving B1, B2, and B6 in (5a), (5b), and (5c) demonstrate certain global deformations that occurred in the laminate. On the other hand, the unknown functions appearing in (5a), (5b), and (5c) illustrate the local deformations. They will be evaluated by LWT of Reddy that happened in a laminate.For cross-ply laminates subjected to mechanical loadings, based on physical grounds, the following restrictions hold (see Figure1):(6a)u1(k)(x,y,z)=-u1(k)(-x,y,z),(6b)u2(k)(x,y,z)=u2(k)(-x,y,z).Upon imposing these conditions on (5a), it is readily seen that u(k)(y,z)=0. Also it is concluded from (6b) and (5b) that B1=0. Thus, for cross-ply laminates the most general form of the displacement field is given:
(7)u1(k)(x,y,z)=B6xz+B2x,u2(k)(x,y,z)=v(k)(y,z),u3(k)(x,y,z)=-12B6x2+w(k)(y,z).Also, for symmetric laminates based on physical grounds the following conditions must hold (see Figure1): (8a)u1(k)(x,y,z)=u1(N+1-k)(x,y,-z)(8b)u2(k)(x,y,z)=u2(N+1-k)(x,y,-z)(8c)u3(k)(x,y,z)=-u3(N+1-k)(x,y,-z)(k=1,2,3,…,N2).Upon imposing these conditions on (5a), (5b), and (5c), it is found that B1=B6=0. Thus, for symmetric laminates the most general form of displacement field is expressed by
(9)u1(k)(x,y,z)=B2x+u(k)(y,z),u2(k)(x,y,z)=v(k)(y,z),u3(k)(x,y,z)=w(k)(y,z).Finally, for antisymmetric angle-ply laminates withN layers subjected to mechanical loadings, the following condition must hold:
(10)u1(k)(x,y,z)=u1(N+1-k)(x,-y,-z)(k=1,2,3,…,N2).
Concluding from (10) and (5a), B6=0 and, therefore, the most general displacement field for such laminates are written as
(11)u1(k)(x,y,z)=B2x+u(k)(y,z),u2(k)(x,y,z)=-B1xz+v(k)(y,z),u3(k)(x,y,z)=B1xy+w(k)(y,z).The displacement fields in (7), (9), and (11) can be represented in one place as
(12)u1(k)(x,y,z)=B2x+δCB6xz+δSu(k)(y,z),u2(k)(x,y,z)=-δAB1xz+v(k)(y,z),u3(k)(x,y,z)=δAB1xy-12δCB6x2+w(k)(y,z),
where δA=δS=0 and δC=1 for general cross-ply laminates, δA=δC=0 and δS=1 for symmetric laminates, δA=δS=1 and δC=0 for angle-ply laminates, and δA=δS=δc=1 for general laminates.


## 3. Layerwise Laminated Plate Theory of Reddy

In this section, LWT is used for calculation of interlaminar stresses in arbitrary laminated composite plates with free edges. In LWT, it is possible to replace the actual physical layers to many mathematical layers [27]. This can be done by subdividing each physical layer through the thickness, to get the proper accuracy. The theory assumes that the displacement components of a generic point in the laminate are given by [27]
(13)u1(x,y,z)=uk(x,y)ϕk(z)u2(x,y,z)=vk(x,y)ϕk(z)u3(x,y,z)=wk(x,y)ϕk(z)k=1,2,…,N+1.
It should be emphasized that a repeated index indicates summation from 1 to N+1. The functions uk(x,y), vk(x,y), and wk(x,y) represent the displacements of the points initially located at the kth plane of the laminate in the x, y, and z directions, respectively. The variable N in (13) is the total number of numerical layers introduced in any laminate. Also the functions ϕk’s are the global approximation functions of the thickness coordinate which are assumed here to be linear [27]. This function expressed as(14a)ϕk(z)={0z≤zk-1ψk-12(z)zk-1≤z≤zkψk1(z)zk≤z≤zk+10z≥zk+1iiiiiii(k=1,2,…,N+1),
where the local Lagrangian interpolation functions ψkj(z)(j=1,2) related with the kth surface in the laminate are defined as
(14b)ψk1(z)=1hk(zk+1-z),ψk2(z)=1hk(z-zk),where hk is the thickness of the kth mathematical layer. Based on the displacement field in (12), the displacement field of LWT in (13) takes the following form:
(15)u1(x,y,z)=δCB6xz+B2x+δSUk(y)ϕk(z)u2(x,y,z)=-δAB1xz+Vk(y)ϕk(z)u3(x,y,z)=δAB1xy-12δCB6x2+Wk(y)ϕk(z)iiiiiiiiiiiiiiiiiiiiiiik=1,2,…,N+1.
It is noted that, by the means of through-the-thickness linear interpolation functions, the continuity of displacement components through the thickness of laminate is identically satisfied. On the other hand, the transverse strain components remain discontinuous at the interfaces. They will ultimately emphasize the feasibility of having continuous interlaminar stresses at the interfaces of adjoining layers by increasing the number of numerical layers through the physical laminate. Substituting (15) into the linear strain-displacement relations of elasticity [25] yields
(16)εx=δCB6z+B2,εy=Vk′ϕk,εz=Wkϕk′,γyz=Vkϕk′+Wk′ϕk,γxz=δSUkϕk′+δAB1y,γxy=δSUk′ϕk-δAB1z,
where a prime in (16) denotes ordinary differentiation with respect to a suitable independent variable (i.e., either y or z). The three-dimensional constitutive law within the kth layer (with arbitrary fiber orientation) of the laminate may also be shown as [28]
(17)(σxσyσzσyzσxzσxy)(k)=[C-11C-12C-1300C-16C-12C-22C-2300C-26C-13C-23C-3300C-36000C-44C-450000C-45C-550C-16C-26C-3600C-66](k)×(εxεyεzγyzγxzγxy)(k),
where [C-](k) is the transformed (i.e., off axis) stiffness matrix of the kth layer. By using the strain-displacement relations (16) in the principle of the minimum total potential energy [25], the equilibrium equations within LWT are found as(18a)δUk:δS(Qxk-dMxykdy)=0k=1,2,…,N+1,(18b)δVk:Qyk-dMykdy=0k=1,2,…,N+1,(18c)δWk:Nzk-dRykdy=0k=1,2,…,N+1,(19a)
δ
B
1
:
δ
A
∫
-
b
b
(
Q
x
y
-
M
x
y
)
d
y
=
0
,
(19b)
δ
B
2
:
∫
-
b
b
N
x
d
y
=
0
,
(19c)
δ
B
6
:
δ
C
∫
-
b
b
M
x
d
y
=
M
0
.
Also the traction-free boundary conditions at the free edges of the laminate are given as
(20)Myk=Ryk=δSMxyk=0aty=±b(k=1,2,…,N+1).
Equations (18a), (18b), and (18c) present 3(N+1) local equilibrium equations associated with N+1 surface in the laminate within LWT. Also, (19a), (19b), and (19c) represent the global equilibrium equation of a laminate. The generalized stress and moment resultants in (18a), (18b), (18c), (19a), (19b), (19c), and (20) are defined as(21a)(Myk,Mxyk,Nzk)=∫-h/2h/2(σyϕk,σxyϕk,σzϕ′)dz,(21b)(Qxk,Qyk,Ryk)=∫-h/2h/2(σxzϕk′,σyzϕk′,σyzϕk)dz,(21c)(Mx,Nx,Mxy,Qx)=∫-h/2h/2(σxz,σx,σxyz,σxz)dz.Substitution (16) into (17) and subsequently into (21a), (21b), and (21c) generates the following relations:(22a)(Myk,Mxyk,Nzk)=δS(D26kj,D66kj,B36kj)Uj′+(D22kj,D26kj,B23jk)Vj′+(B23kj,B36kj,A33kj)Wj-δA(D26k,D66k,B-36k)B1+δC(D12k,D16k,B-13k)B6+(B12k,B16k,A13k)B2,(22b)(Qxk,Qyk,Ryk)=δS(A55kj,A45kj,B45kj)Uj+(A45kj,A44kj,B44kj)Vj+(B45jk,B44jk,D44kj)Wj′-δA(A55k,A45k,B45k)B1y,(22c)(My,Mxy,Nx)=δS(D16k,D66k,B16k)Uj′+(D12k,D26k,B12k)Vj′+(B13k,B-36k,A13k)Wj-δA(D16,D66,B16)B1+δC(B11,B16,A11)B6+(D11,D16,B11)B2,(22d)Qx=δSA55kUj+A45kVj+B45kWj′+δAA55B1y.The explicit expressions for the laminate rigidities appearing in (22a), (22b), (22c), and (22d) in terms of C-i’s are presented in the Appendix. Direct substitution of (22a), (22b), (22c), and (22d) into (18a), (18b), and (18c) provides the local displacement equilibrium equations within LWT. The results are expressed as(23a)δUk:δS(D66kjUj′′-A55kjUj+D26kjVj′′000000000-A45kjVj+(B36kj-B45jk)Wj′)=δAA55kB1y0000000000000000000000000000k=1,2,…,N+1.(23b)δVk:δSD26kjUj′′-δSA45kjUj+D22kjVj′′-A44kjVj+(B23kj-B44jk)Wj′=δAA45kB1y0000000000000000000k=1,2,…,N+1,(23c)δWk:δS(B45kj-B36jk)Uj′-(B44kj-B23jk)Vj′+D44kjWj′′-A33kjWj=-δA(B45k+B-36k)B1+A13kB2+δCB-13kB60000000000000000000k=1,2,…,N+1.In the same way, upon substitution of (16) into (17) and the subsequent results into (19a), (19b), and (19c), the global equations, (19a), (19b), and (19c), are found in terms of displacement functions as follows: (24a)δB1:δA∫-bb((δSA55kUj+A45kVj+B45kWj′)y00000000000-(δSD66kUj′+D26kVj′+B-36kWj)00000000000+δA(A55B1y2+D66)B1-B16B200000000000-δCD16B6(δSA55kUj+A45kVj+B45kWj′))dy=0,(24b)δB2:∫-bb(δSB16kUj′+B12kVj′+A13kWj00000000-δAB16B1+A11B2+δCB11B6B16kUj′)dy=0,(24c)δB6:δC∫-bb(B-13kδSD16kUj′+D12kVj′+B-13kWj0000000000-δAD16B1+B11B2+δCD11B6B-13k1)dy=M0,where the extra laminate rigidities appearing in (24a), (24b), and (24c) are defined in the Appendix. The system of equations in (23a), (23b), and (23c) shows 3(N+1) coupled ordinary differential equations with constant coefficients which may be indicated in a matrix form as
(25)[M]{η′′}+[K]{η}=[L]{B}·y,
where
(26){η}={{U}T,{V}T,{W-}T}T{U}={U1,U2,…,UN+1}T{V}={V1,V2,…,VN+1}T{W-}={W-1,W-2,…,W-N+1}T(27){B}={B1,B2,B6}T,W-j=∫yWjdy.
The coefficient matrices [M], [K], and [L] in (25) are displayed in the Appendix. Also from the conditions in (3) and the displacement field in (23a), (23b), and (23c) it is observed that the functions Uj, Vj, and W-j are all odd functions of the independent variable y. It can be confirmed that the general solution of (25) may be presented as
(28){η}=[ψ][sinh(λy)]{H}+[K]-1[L]{B}·y,
where [sinh(λy)] is a 3(N+1)×3(N+1) diagonal matrix, and
(29)[sinh(λy)]=diag(sinh(λ1y),sinh(λ2y),…,sinh(λ3(N+1)y)).Also [ψ] and (λ12,λ22,…,λ3(N+1)2) are the model matrix and eigenvalues of (-[M]-1[K]), respectively. {K} is an unknown vector representing 3(N+1) integration constants. The constants Bj(j=1,2,6) must be calculated within LWT analysis. So, the boundary conditions in (20) are first imposed to found the vector {K} in terms of the unknown parameters Bj(j=1,2,6). These constants are next calculated in terms of the specific bending moment M0 by satisfaction of the global equilibrium conditions in (19a), (19b), and (19c).


## 4. Numerical Results and Discussions

In this part, several numerical examples are shown for the distribution of interlaminar stress. All physical layers are assumed to be of equal thickness (=0.5 mm). The mechanical properties of unidirectional graphite/epoxy T300/5208 used in this study are shown in Table1 [28].Table 1
Engineering properties of unidirectional graphite/opoxy [28].
E
1  (GPa)
E
2
=
E
3 (GPa)
G
12
=
G
13 (GPa)
G
23 (GPa)
υ
12
=
υ
13
υ
23
132
10.8
5.65
3.38
0.24
0.59Also, each physical ply is modeled as being made up of 12 numerical layers within LWT (i.e.,p=12). In addition, the width thickness ratio (i.e., 2b/h) is assumed to be equal to 10. Furthermore, the stress components are normalized as
(30)σ-ij=σijσ0,
where σ0=M0/bh2.Figure2 shows the distribution of interlaminar normal and shear stresses along the width of [45°/−45°]s laminate. Good correspondence was found between the present results and the results obtained by Goodsell and Pipes [29] which confirms the accuracy of the present solution. Figure 3 shows the distribution of interlaminar stresses σz and σyz along the −45°/0° and σxz along the 0°/90° interfaces of the general [45°/−45°/0°/90°/30°/−30°] laminate. It is seen that the interlaminar stresses show high stress gradient near the free edge. It is observed that the interlaminar stresses σz and σxz grow quickly near the free edge, while being zero in the inner region of the laminate. Also the magnitude of the transverse shear stress σxz is greater than that of transverse normal stress σz. On the other hand, σyz rises toward the free edge and decreases suddenly to zero at free edge.Figure 2
Interlaminar stresses along the 45°/−45° interface of [45°/−45°]s laminate.Figure 3
Distribution of interlaminar stressesσz and σyz along the −45°/0° and σxz along the 0°/90° interfaces of the [45°/−45°/0°/90°/30°/−30°] laminate.The interlaminar normal stressσz along the 90°/0° interface of the unsymmetric cross-ply [90°/0°/90°/0°/90°/0°] laminate for various width-to-thickness ratios is demonstrated in Figure 4. It is observed that by decreasing the width-to-thickness ratio, the boundary-layer region can be expanded towards the internal region of the laminate with its width being almost equal to the thickness of the laminate. It is seen that the magnitude of the interlaminar stress at the free edge is not varied as the width-to-thickness ratio of the laminate is changed. In addition, the highly localized nature of interlaminar stresses near and exactly at the free edges of the laminate is apparent.Figure 4
Distribution of the interlaminar normal stressσz along the 90°/0° interface of the laminate for various width-to-thickness ratios.The distribution of interlaminar stresses along the 0°/90° and 90°/0° interfaces of the unsymmetric cross-ply [0°/90°/0°/90°] and [90°/0°/90°/0°] laminates, respectively, are compared in Figure 5. It is seen that the maximum interfacial value σz at the 0°/90° and 90°/0° interfaces take place at the free edge. Figure 6 shows the variation of interlaminar normal stress σz on through-the-thickness and near the free edge of the angle-ply [45°/−45°/45°/−45°] laminate as the free edge is approached. It is observed that the maximum negative and positive values of σz occur within the top -80° layer and the bottom 80° layer at the free edge (i.e., y=b), respectively. It is also noticed that σz diminishes away from the free edge as the interior region of the laminate is approached. Moreover, σz reduces by moving slightly away from the free edge and its distribution becomes smoother.Figure 5
Distribution of interlaminar stresses along the 0°/90° and 90°/0° interfaces of the [0°/90°/0°/90°] and [90°/0°/90°/0°] laminates, respectively.Figure 6
Variations of interlaminar normal stressσz through the thickness of the [45°/−45°/45°/−45°] laminate for various width-to-thickness ratios.The variations of the interlaminar shear stressσyz at y=b through the thickness of the [0°/60°/−60°]s  laminate are displayed in Figure 7. By increasing the number of numerical layers in each lamina, σyz becomes slightly closer to zero but it never becomes zero. It should be noted that increasing the number of subdivisions results in no convergence for σz and σxz at interface-edge junction of two dissimilar layers. The numerical values of these components (σz and σxz) continue to grow as the number of sublayers increased. Since the generalized stress resultants Ryk instead of σyz are forced to disappear at the free edge in LWT, the numerical value of σyz may never become zero at the interface-edge intersection (even by increasing the number of sublayers in each physical layer).Figure 7
Variations of interlaminar shear stressσyz on through-thickness of the [0°/60°/−60°]s laminate as a function of layer subdivision number p.The distribution of interlaminar normal stressσz on through thickness of the general [−15°/90°/90°/−15°] laminate is displayed in Figure 8. It is noted that increasing the number of numerical layers (p) in each lamina has no significant effect on the values of the interlaminar stress σz within the boundary-layer region of the laminate, expect exactly at the free edge (i.e., y=b). Furthermore, the interlaminar normal stress σz grows rapidly in the vicinity of the free edge, while being zero in the interior region of the laminate.Figure 8
Variations of interlaminar normal stressσz on through-thickness of the [−15°/90°/90°/−15°] laminate as a function of layer subdivision number p.


## 5. Conclusions

In the present paper, analytical solutions are developed to analyze the interlaminar stresses for different lay-up configurations in laminated composite plates subjected to the bending moment. The solutions were obtained based on the reduced elasticity displacement field for long laminated plates. It is found that the components of the displacements field are composed of two distinct parts, signifying the global and local deformations within laminates. The local displacement functions as well as the interlaminar stresses through the layers of the laminate were found by Reddy’s layerwise theory (LWT). The numerical results based on LWT were shown for the free-edge interlaminar stresses through the thickness and across the interfaces of various cross-ply, symmetric, angle-ply, and general composite laminates.



---

*Source: 676743-2013-12-12.xml*