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

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.

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 𝐶^∘ 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

An Nth-layered composite plate (with arbitrary lamination) under the bending moment is considered as shown in Figure 1. The coordinate system (𝑥,𝑦,𝑧) is located at the middle plane of the laminate, that is of thickness ℎ, width 2𝑏, and length 2𝑎 and is assumed to be long in the 𝑥 direction so that the strains away from the ends (𝑥=±𝑎) of the laminate are functions of only 𝑦 and 𝑧.

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] 𝑢₁^(𝑘)𝑥,𝑦,𝑧=𝐵₄^(𝑘)𝑥𝑦+𝐵₆^(𝑘)𝑥𝑧+𝐵₂^(𝑘)𝑥+𝑢^(𝑘)𝑢𝑦,𝑧_2(𝑘)𝑥,𝑦,𝑧=−𝐵₁^(𝑘)𝑥𝑧+𝐵₃^(𝑘)1𝑥−2𝐵₄^(𝑘)𝑥²+𝑣^(𝑘)𝑢𝑦,𝑧_3(𝑘)𝑥,𝑦,𝑧=𝐵₁^(𝑘)𝑥𝑦+𝐵₅^(𝑘)1𝑥−2𝐵₆^(𝑘)𝑥²+𝑤^(𝑘)𝑦,𝑧, (1) where 𝑢₁^(𝑘), 𝑢₂^(𝑘), and 𝑢₃^(𝑘) represent the displacement components of the material point (𝑥,𝑦,𝑧) in the 𝑥, 𝑦, and 𝑧 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(𝑘)𝑥,𝑦,𝑧=−𝐵₁𝑥𝑧+𝐵₃1𝑥−2𝐵₄𝑥²+𝑣^(𝑘),𝑢𝑦,𝑧_3(𝑘)𝑥,𝑦,𝑧=𝐵₁𝑥𝑦+𝐵₅1𝑥−2𝐵₆𝑥²+𝑤^(𝑘)𝑦,𝑧. (2) If loading conditions at 𝑥=−𝑎 and 𝑎 are similar, based on the physical grounds, the following conditions hold: 𝑢₁^(𝑘)𝑥,𝑦,𝑧=−𝑢₁^(𝑘),𝑢−𝑥,−𝑦,𝑧_2(𝑘)𝑥,𝑦,𝑧=−𝑢₂^(𝑘),𝑢−𝑥,−𝑦,𝑧_3(𝑘)𝑥,𝑦,𝑧=𝑢₃^(𝑘)−𝑥,−𝑦,𝑧. (3) Upon imposing these restrictions on (1) it is readily seen that the constants 𝐵₄ and 𝐵₅ must vanish and the displacement field in (1) is reduced to the following equations:𝑢₁^(𝑘)𝑥,𝑦,𝑧=𝐵₂𝑥+𝐵₆𝑥𝑧+𝑢^(𝑘),𝑢𝑦,𝑧_2(𝑘)𝑥,𝑦,𝑧=−𝐵₁𝑥𝑧+𝐵₃𝑥+𝑣^(𝑘),𝑢𝑦,𝑧_3(𝑘)𝑥,𝑦,𝑧=𝐵₁1𝑥𝑦−2𝐵₆𝑥²+𝑤^(𝑘)𝑦,𝑧. (4a)(4b)(4c)Also, by replacing 𝑢^(𝑘)(𝑦,𝑧) by −𝐵₃𝑦+𝑢^(𝑘)(𝑦,𝑧) in (4a), it becomes apparent that terms involving 𝐵₃ 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𝑢₁^(𝑘)𝑥,𝑦,𝑧=𝐵₂𝑥+𝐵₆𝑥𝑧+𝑢^(𝑘),𝑢𝑦,𝑧_2(𝑘)𝑥,𝑦,𝑧=−𝐵₁𝑥𝑧+𝑣^(𝑘),𝑢𝑦,𝑧_3(𝑘)𝑥,𝑦,𝑧=𝐵₁1𝑥𝑦−2𝐵₆𝑥²+𝑤^(𝑘)𝑦,𝑧. (5a)(5b)(5c)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 𝐵₁, 𝐵₂, and 𝐵₆ 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 Figure 1):𝑢₁^(𝑘)𝑥,𝑦,𝑧=−𝑢₁^(𝑘),𝑢−𝑥,𝑦,𝑧_2(𝑘)𝑥,𝑦,𝑧=𝑢₂^(𝑘)−𝑥,𝑦,𝑧. (6a)(6b)Upon imposing these conditions on (5a), it is readily seen that 𝑢^(𝑘)(𝑦,𝑧)=0. Also it is concluded from (6b) and (5b) that 𝐵₁=0. Thus, for cross-ply laminates the most general form of the displacement field is given: 𝑢₁^(𝑘)𝑥,𝑦,𝑧=𝐵₆𝑥𝑧+𝐵₂𝑢𝑥,_2(𝑘)𝑥,𝑦,𝑧=𝑣^(𝑘),𝑢𝑦,𝑧_3(𝑘)1𝑥,𝑦,𝑧=−2𝐵₆𝑥²+𝑤^(𝑘)𝑦,𝑧. (7)

Also, for symmetric laminates based on physical grounds the following conditions must hold (see Figure 1): 𝑢₁^(𝑘)𝑥,𝑦,𝑧=𝑢₁^{(𝑁+1−𝑘)}𝑢𝑥,𝑦,−𝑧_2(𝑘)𝑥,𝑦,𝑧=𝑢₂^{(𝑁+1−𝑘)}𝑢𝑥,𝑦,−𝑧_3(𝑘)𝑥,𝑦,𝑧=−𝑢₃^{(𝑁+1−𝑘)}𝑁𝑥,𝑦,−𝑧𝑘=1,2,3,…,2. (8a)(8b)(8c)Upon imposing these conditions on (5a), (5b), and (5c), it is found that 𝐵₁=𝐵₆=0. Thus, for symmetric laminates the most general form of displacement field is expressed by 𝑢₁^(𝑘)𝑥,𝑦,𝑧=𝐵₂𝑥+𝑢^(𝑘),𝑢𝑦,𝑧_2(𝑘)𝑥,𝑦,𝑧=𝑣^(𝑘),𝑢𝑦,𝑧_3(𝑘)𝑥,𝑦,𝑧=𝑤^(𝑘)𝑦,𝑧. (9)

Finally, for antisymmetric angle-ply laminates with 𝑁 layers subjected to mechanical loadings, the following condition must hold: 𝑢₁^(𝑘)𝑥,𝑦,𝑧=𝑢₁^{(𝑁+1−𝑘)}𝑁𝑥,−𝑦,−𝑧𝑘=1,2,3,…,2. (10) Concluding from (10) and (5a), 𝐵₆=0 and, therefore, the most general displacement field for such laminates are written as 𝑢₁^(𝑘)𝑥,𝑦,𝑧=𝐵₂𝑥+𝑢^(𝑘),𝑢𝑦,𝑧_2(𝑘)𝑥,𝑦,𝑧=−𝐵₁𝑥𝑧+𝑣^(𝑘),𝑢𝑦,𝑧_3(𝑘)𝑥,𝑦,𝑧=𝐵₁𝑥𝑦+𝑤^(𝑘)𝑦,𝑧. (11)

The displacement fields in (7), (9), and (11) can be represented in one place as 𝑢₁^(𝑘)𝑥,𝑦,𝑧=𝐵₂𝑥+𝛿_𝐶𝐵₆𝑥𝑧+𝛿_𝑆𝑢^(𝑘),𝑢𝑦,𝑧_2(𝑘)𝑥,𝑦,𝑧=−𝛿_𝐴𝐵₁𝑥𝑧+𝑣^(𝑘),𝑢𝑦,𝑧_3(𝑘)𝑥,𝑦,𝑧=𝛿_𝐴𝐵₁1𝑥𝑦−2𝛿_𝐶𝐵₆𝑥²+𝑤^(𝑘)𝑦,𝑧, (12) where 𝛿_𝐴=𝛿_𝑆=0 and 𝛿_𝐶=1 for general cross-ply laminates, 𝛿_𝐴=𝛿_𝐶=0 and 𝛿_𝑆=1 for symmetric laminates, 𝛿_𝐴=𝛿_𝑆=1 and 𝛿_𝐶=0 for angle-ply laminates, and 𝛿_𝐴=𝛿_𝑆=𝛿_𝑐=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] 𝑢₁𝑥,𝑦,𝑧=𝑢_𝑘𝑥,𝑦𝜙_𝑘(𝑧)𝑢₂𝑥,𝑦,𝑧=𝑣_𝑘𝑥,𝑦𝜙_𝑘(𝑧)𝑢₃𝑥,𝑦,𝑧=𝑤_𝑘𝑥,𝑦𝜙_𝑘(𝑧)𝑘=1,2,…,𝑁+1. (13) It should be emphasized that a repeated index indicates summation from 1 to 𝑁+1. The functions 𝑢_𝑘(𝑥,𝑦), 𝑣_𝑘(𝑥,𝑦), and 𝑤_𝑘(𝑥,𝑦) represent the displacements of the points initially located at the kth plane of the laminate in the 𝑥, 𝑦, and 𝑧 directions, respectively. The variable 𝑁 in (13) is the total number of numerical layers introduced in any laminate. Also the functions 𝜙_𝑘’s are the global approximation functions of the thickness coordinate which are assumed here to be linear [27]. This function expressed as𝜙_𝑘(𝑧)=0𝑧≤𝑧_𝑘−1𝜓²_𝑘−1(𝑧)𝑧_𝑘−1≤𝑧≤𝑧_𝑘𝜓¹_𝑘(𝑧)𝑧_𝑘≤𝑧≤𝑧_𝑘+10𝑧≥𝑧_𝑘+1(𝑘=1,2,…,𝑁+1), (14a) where the local Lagrangian interpolation functions 𝜓^𝑗_𝑘(𝑧)(𝑗=1,2) related with the kth surface in the laminate are defined as 𝜓¹_𝑘(𝑧)=1ℎ_𝑘𝑧_𝑘+1−𝑧,𝜓²_𝑘(𝑧)=1ℎ_𝑘𝑧−𝑧_𝑘, (14b)where ℎ_𝑘 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: 𝑢₁𝑥,𝑦,𝑧=𝛿_𝐶𝐵₆𝑥𝑧+𝐵₂𝑥+𝛿_𝑆𝑈_𝑘𝑦𝜙_𝑘(𝑧)𝑢₂𝑥,𝑦,𝑧=−𝛿_𝐴𝐵₁𝑥𝑧+𝑉_𝑘𝑦𝜙_𝑘(𝑧)𝑢₃𝑥,𝑦,𝑧=𝛿_𝐴𝐵₁1𝑥𝑦−2𝛿_𝐶𝐵₆𝑥²+𝑊_𝑘𝑦𝜙_𝑘(𝑧)𝑘=1,2,…,𝑁+1. (15) 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) where a prime in (16) denotes ordinary differentiation with respect to a suitable independent variable (i.e., either 𝑦 or 𝑧). The three-dimensional constitutive law within the kth layer (with arbitrary fiber orientation) of the laminate may also be shown as [28] 𝜎_𝑥𝜎_𝑦𝜎_𝑧𝜎_𝑦𝑧𝜎_𝑥𝑧𝜎_𝑥𝑦^(𝑘)=𝐶₁₁𝐶₁₂𝐶₁₃00𝐶₁₆𝐶₁₂𝐶₂₂𝐶₂₃00𝐶₂₆𝐶₁₃𝐶₂₃𝐶₃₃00𝐶₃₆000𝐶₄₄𝐶₄₅0000𝐶₄₅𝐶₅₅0𝐶₁₆𝐶₂₆𝐶₃₆00𝐶₆₆^(𝑘)×𝜀_𝑥𝜀_𝑦𝜀_𝑧𝛾_𝑦𝑧𝛾_𝑥𝑧𝛾_𝑥𝑦^(𝑘), (17) where [𝐶]^(𝑘) 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𝛿𝑈_𝑘:𝛿_𝑆𝑄^𝑘_𝑥−d𝑀^𝑘_𝑥𝑦d𝑦=0𝑘=1,2,…,𝑁+1,𝛿𝑉_𝑘:𝑄^𝑘_𝑦−d𝑀^𝑘_𝑦d𝑦=0𝑘=1,2,…,𝑁+1,𝛿𝑊_𝑘:𝑁^𝑘_𝑧−d𝑅^𝑘_𝑦d𝑦=0𝑘=1,2,…,𝑁+1, (18a)(18b)(18c)

𝛿𝐵₁:𝛿_𝐴^𝑏_−𝑏𝑄_𝑥𝑦−𝑀_𝑥𝑦d𝑦=0,𝛿𝐵₂:^𝑏_−𝑏𝑁_𝑥d𝑦=0,𝛿𝐵₆:𝛿_𝐶^𝑏_−𝑏𝑀_𝑥d𝑦=𝑀₀. (19a)(19b)(19c) Also the traction-free boundary conditions at the free edges of the laminate are given as 𝑀^𝑘_𝑦=𝑅^𝑘_𝑦=𝛿_𝑆𝑀^𝑘_𝑥𝑦=0at𝑦=±𝑏(𝑘=1,2,…,𝑁+1). (20) Equations (18a), (18b), and (18c) present 3(𝑁+1) local equilibrium equations associated with 𝑁+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𝑀^𝑘_𝑦,𝑀^𝑘_𝑥𝑦,𝑁^𝑘_𝑧=^ℎ/2_−ℎ/2𝜎_𝑦𝜙_𝑘,𝜎_𝑥𝑦𝜙_𝑘,𝜎_𝑧𝜙^′𝑄d𝑧,_𝑘𝑥,𝑄^𝑘_𝑦,𝑅^𝑘_𝑦=^ℎ/2_−ℎ/2𝜎_𝑥𝑧𝜙^′_𝑘,𝜎_𝑦𝑧𝜙^′_𝑘,𝜎_𝑦𝑧𝜙_𝑘𝑀d𝑧,_𝑥,𝑁_𝑥,𝑀_𝑥𝑦,𝑄_𝑥=^ℎ/2_−ℎ/2𝜎_𝑥𝑧,𝜎_𝑥,𝜎_𝑥𝑦𝑧,𝜎_𝑥𝑧d𝑧. (21a)(21b)(21c)Substitution (16) into (17) and subsequently into (21a), (21b), and (21c) generates the following relations:𝑀^𝑘_𝑦,𝑀^𝑘_𝑥𝑦,𝑁^𝑘_𝑧=𝛿_𝑆𝐷^𝑘𝑗₂₆,𝐷^𝑘𝑗₆₆,𝐵^𝑘𝑗₃₆𝑈^′_𝑗+𝐷^𝑘𝑗₂₂,𝐷^𝑘𝑗₂₆,𝐵^𝑗𝑘₂₃𝑉^′_𝑗+𝐵^𝑘𝑗₂₃,𝐵^𝑘𝑗₃₆,𝐴^𝑘𝑗₃₃𝑊_𝑗−𝛿_𝐴𝐷^𝑘₂₆,𝐷^𝑘₆₆,𝐵^𝑘₃₆𝐵₁+𝛿_𝐶𝐷^𝑘₁₂,𝐷^𝑘₁₆,𝐵^𝑘₁₃𝐵₆+𝐵^𝑘₁₂,𝐵^𝑘₁₆,𝐴^𝑘₁₃𝐵₂,𝑄^𝑘_𝑥,𝑄^𝑘_𝑦,𝑅^𝑘_𝑦=𝛿_𝑆𝐴^𝑘𝑗₅₅,𝐴^𝑘𝑗₄₅,𝐵^𝑘𝑗₄₅𝑈_𝑗+𝐴^𝑘𝑗₄₅,𝐴^𝑘𝑗₄₄,𝐵^𝑘𝑗₄₄𝑉_𝑗+𝐵^𝑗𝑘₄₅,𝐵^𝑗𝑘₄₄,𝐷^𝑘𝑗₄₄𝑊^′_𝑗−𝛿_𝐴𝐴^𝑘₅₅,𝐴^𝑘₄₅,𝐵^𝑘₄₅𝐵₁𝑀𝑦,_𝑦,𝑀_𝑥𝑦,𝑁_𝑥=𝛿_𝑆𝐷^𝑘₁₆,𝐷^𝑘₆₆,𝐵^𝑘₁₆𝑈^′_𝑗+𝐷^𝑘₁₂,𝐷^𝑘₂₆,𝐵^𝑘₁₂𝑉^′_𝑗+𝐵^𝑘₁₃,𝐵^𝑘₃₆,𝐴^𝑘₁₃𝑊_𝑗−𝛿_𝐴𝐷₁₆,𝐷₆₆,𝐵₁₆𝐵₁+𝛿_𝐶𝐵₁₁,𝐵₁₆,𝐴₁₁𝐵₆+𝐷₁₁,𝐷₁₆,𝐵₁₁𝐵₂,𝑄_𝑥=𝛿_𝑆𝐴^𝑘₅₅𝑈_𝑗+𝐴^𝑘₄₅𝑉_𝑗+𝐵^𝑘₄₅𝑊^′_𝑗+𝛿_𝐴𝐴₅₅𝐵₁𝑦. (22a)(22b)(22c)(22d)The explicit expressions for the laminate rigidities appearing in (22a), (22b), (22c), and (22d) in terms of 𝐶_𝑖’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𝛿𝑈_𝑘:𝛿_𝑆𝐷^𝑘𝑗₆₆𝑈^′′_𝑗−𝐴^𝑘𝑗₅₅𝑈_𝑗+𝐷^𝑘𝑗₂₆𝑉^′′_𝑗−𝐴^𝑘𝑗₄₅𝑉_𝑗+𝐵^𝑘𝑗₃₆−𝐵^𝑗𝑘₄₅𝑊^′_𝑗=𝛿_𝐴𝐴^𝑘₅₅𝐵₁𝑦𝑘=1,2,…,𝑁+1.𝛿𝑉_𝑘:𝛿_𝑆𝐷^𝑘𝑗₂₆𝑈^′′_𝑗−𝛿_𝑆𝐴^𝑘𝑗₄₅𝑈_𝑗+𝐷^𝑘𝑗₂₂𝑉^′′_𝑗−𝐴^𝑘𝑗₄₄𝑉_𝑗+𝐵^𝑘𝑗₂₃−𝐵^𝑗𝑘₄₄𝑊^′_𝑗=𝛿_𝐴𝐴^𝑘₄₅𝐵₁𝑦𝑘=1,2,…,𝑁+1,𝛿𝑊_𝑘:𝛿_𝑆𝐵^𝑘𝑗₄₅−𝐵^𝑗𝑘₃₆𝑈^′_𝑗−𝐵^𝑘𝑗₄₄−𝐵^𝑗𝑘₂₃𝑉^′_𝑗+𝐷^𝑘𝑗₄₄𝑊^′′_𝑗−𝐴^𝑘𝑗₃₃𝑊_𝑗=−𝛿_𝐴𝐵^𝑘₄₅+𝐵^𝑘₃₆𝐵₁+𝐴^𝑘₁₃𝐵₂+𝛿_𝐶𝐵^𝑘₁₃𝐵₆𝑘=1,2,…,𝑁+1. (23a)(23b)(23c)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: 𝛿𝐵₁:𝛿_𝐴^𝑏_−𝑏𝛿_𝑆𝐴^𝑘₅₅𝑈_𝑗+𝐴^𝑘₄₅𝑉_𝑗+𝐵^𝑘₄₅𝑊^′_𝑗𝑦−𝛿_𝑆𝐷^𝑘₆₆𝑈^′_𝑗+𝐷^𝑘₂₆𝑉^′_𝑗+𝐵^𝑘₃₆𝑊_𝑗+𝛿_𝐴𝐴₅₅𝐵₁𝑦²+𝐷₆₆𝐵₁−𝐵₁₆𝐵₂−𝛿_𝐶𝐷₁₆𝐵₆d𝑦=0,𝛿𝐵₂:^𝑏_−𝑏𝛿_𝑆𝐵^𝑘₁₆𝑈^′_𝑗+𝐵^𝑘₁₂𝑉^′_𝑗+𝐴^𝑘₁₃𝑊_𝑗−𝛿_𝐴𝐵₁₆𝐵₁+𝐴₁₁𝐵₂+𝛿_𝐶𝐵₁₁𝐵₆d𝑦=0,𝛿𝐵₆:𝛿_𝐶^𝑏_−𝑏𝛿_𝑆𝐷^𝑘₁₆𝑈^′_𝑗+𝐷^𝑘₁₂𝑉^′_𝑗+𝐵^𝑘₁₃𝑊_𝑗−𝛿_𝐴𝐷₁₆𝐵₁+𝐵₁₁𝐵₂+𝛿_𝐶𝐷₁₁𝐵₆d𝑦=𝑀₀, (24a)(24b)(24c)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(𝑁+1) coupled ordinary differential equations with constant coefficients which may be indicated in a matrix form as [𝑀]𝜂^′′+[𝐾]𝜂=[𝐿]{𝐵}⋅𝑦, (25) where 𝜂={𝑈}^𝑇,{𝑉}^𝑇,𝑊^𝑇𝑇{𝑈}=𝑈₁,𝑈₂,…,𝑈_𝑁+1^𝑇{𝑉}=𝑉₁,𝑉₂,…,𝑉_𝑁+1^𝑇𝑊=𝑊₁,𝑊₂,…,𝑊_𝑁+1^𝑇{𝐵}=𝐵₁,𝐵₂,𝐵₆^𝑇,𝑊_𝑗=^𝑦𝑊_𝑗d𝑦. (26)(27) The coefficient matrices [𝑀], [𝐾], and [𝐿] 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 𝑈_𝑗, 𝑉_𝑗, and 𝑊_𝑗 are all odd functions of the independent variable 𝑦. It can be confirmed that the general solution of (25) may be presented as 𝜂=𝜓sinh𝜆𝑦{𝐻}+[𝐾]⁻¹[𝐿]{𝐵}⋅𝑦, (28) where [sinh(𝜆𝑦)] is a 3(𝑁+1)×3(𝑁+1) diagonal matrix, and 𝜆sinh𝜆𝑦=diagsinh₁𝑦,sinh𝜆₂𝑦,…,sinh𝜆_3(𝑁+1)𝑦. (29)

Also [𝜓] and (𝜆²₁,𝜆²₂,…,𝜆²_3(𝑁+1)) are the model matrix and eigenvalues of (−[𝑀]⁻¹[𝐾]), respectively. {𝐾} is an unknown vector representing 3(𝑁+1) integration constants. The constants 𝐵_𝑗(𝑗=1,2,6) must be calculated within LWT analysis. So, the boundary conditions in (20) are first imposed to found the vector {𝐾} in terms of the unknown parameters 𝐵_𝑗(𝑗=1,2,6). These constants are next calculated in terms of the specific bending moment 𝑀₀ 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 Table 1 [28].

Table 1 

Engineering properties of unidirectional graphite/opoxy [28].

Also, each physical ply is modeled as being made up of 12 numerical layers within LWT (i.e., 𝑝=12). In addition, the width thickness ratio (i.e., 2𝑏/ℎ) is assumed to be equal to 10. Furthermore, the stress components are normalized as 𝜎_𝑖𝑗=𝜎_𝑖𝑗𝜎₀, (30) where 𝜎₀=𝑀₀/𝑏ℎ².

Figure 2 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 𝜎_𝑧 and 𝜎_𝑦𝑧 along the −45°/0° and 𝜎_𝑥𝑧 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 𝜎_𝑧 and 𝜎_𝑥𝑧 grow quickly near the free edge, while being zero in the inner region of the laminate. Also the magnitude of the transverse shear stress 𝜎_𝑥𝑧 is greater than that of transverse normal stress 𝜎_𝑧. On the other hand, 𝜎_𝑦𝑧 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 𝜎𝑧 and 𝜎𝑦𝑧 along the −45°/0° and 𝜎𝑥𝑧 along the 0°/90° interfaces of the [45°/−45°/0°/90°/30°/−30°] laminate.

The interlaminar normal stress 𝜎_𝑧 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 𝜎𝑧 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 𝜎_𝑧 at the 0^∘/90^∘ and 90^∘/0^∘ interfaces take place at the free edge. Figure 6 shows the variation of interlaminar normal stress 𝜎_𝑧 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 𝜎_𝑧 occur within the top −80^∘ layer and the bottom 80^∘ layer at the free edge (i.e., 𝑦=𝑏), respectively. It is also noticed that 𝜎_𝑧 diminishes away from the free edge as the interior region of the laminate is approached. Moreover, 𝜎_𝑧 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 𝜎𝑧 through the thickness of the [45°/−45°/45°/−45°] laminate for various width-to-thickness ratios.

The variations of the interlaminar shear stress 𝜎_𝑦𝑧 at 𝑦=𝑏 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, 𝜎_𝑦𝑧 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 𝜎_𝑧 and 𝜎_𝑥𝑧 at interface-edge junction of two dissimilar layers. The numerical values of these components (𝜎_𝑧 and 𝜎_𝑥𝑧) continue to grow as the number of sublayers increased. Since the generalized stress resultants 𝑅^𝑘_𝑦 instead of 𝜎_𝑦𝑧 are forced to disappear at the free edge in LWT, the numerical value of 𝜎_𝑦𝑧 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 𝜎𝑦𝑧 on through-thickness of the [0°/60°/−60°]s laminate as a function of layer subdivision number 𝑝.

The distribution of interlaminar normal stress 𝜎_𝑧 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 (𝑝) in each lamina has no significant effect on the values of the interlaminar stress 𝜎_𝑧 within the boundary-layer region of the laminate, expect exactly at the free edge (i.e., 𝑦=𝑏). Furthermore, the interlaminar normal stress 𝜎_𝑧 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 𝜎𝑧 on through-thickness of the [−15°/90°/90°/−15°] laminate as a function of layer subdivision number 𝑝.

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.

Appendix

The laminate rigidities in (22a), (22b), (22c), (22d), (24a), (24b), and (24c) are defined as 𝐴^𝑘𝑗_𝑝𝑞,𝐵^𝑘𝑗_𝑝𝑞,𝐷^𝑘𝑗_𝑝𝑞=^𝑁_𝑖=1^ℎ/2_−ℎ/2𝐶^(𝑖)_𝑝𝑞𝜙^′_𝑘𝜙^′_𝑗,𝜙_𝑘𝜙^′_𝑗,𝜙_𝑘𝜙_𝑗𝐴d𝑧_𝑘𝑝𝑞,𝐵^𝑘_𝑝𝑞,𝐵^𝑘_𝑝𝑞,𝐷^𝑘_𝑝𝑞=^𝑁_𝑖=1^ℎ/2_−ℎ/2𝐶^(𝑖)_𝑝𝑞𝜙^′_𝑘,𝜙_𝑘,𝜙^′_𝑘𝑧,𝜙_𝑘𝑧𝐴d𝑧𝑘,𝑗=1,2,…,𝑁+1_𝑝𝑞,𝐵_𝑝𝑞,𝐷_𝑝𝑞=^𝑁_𝑖=1^𝑧𝑖+1_𝑧𝑖𝐶^(𝑖)_𝑝𝑞1,𝑧,𝑧²d𝑧 (A.1) which, upon integration, are presented in the following form:𝐴^𝑘𝑗_𝑝𝑞,𝐵^𝑘𝑗_𝑝𝑞,𝐷^𝑘𝑗_𝑝𝑞=−𝐶^(𝑘−1)_𝑝𝑞ℎ_𝑘−1,−𝐶^(𝑘−1)_𝑝𝑞2,ℎ_𝑘−1𝐶^(𝑘−1)_𝑝𝑞6if𝑗=𝑘−1𝐶^(𝑘−1)_𝑝𝑞ℎ_𝑘−1+𝐶^(𝑘)_𝑝𝑞ℎ_𝑘,𝐶^(𝑘−1)_𝑝𝑞2−𝐶^(𝑘)_𝑝𝑞2,ℎ_𝑘−1𝐶^(𝑘−1)_𝑝𝑞3+ℎ_𝑘𝐶^(𝑘)_𝑝𝑞3−if𝑗=𝑘𝐶^(𝑘)_𝑝𝑞ℎ_𝑘,𝐶^(𝑘)_𝑝𝑞2,ℎ_𝑘𝐶^(𝑘)_𝑝𝑞6()𝐴if𝑗=𝑘+10,0,0if𝑗<𝑘−1or𝑗>𝑘+1,_𝑘𝑝𝑞,𝐵^𝑘_𝑝𝑞,𝐵^𝑘_𝑝𝑞,𝐷^𝑘_𝑝𝑞=−𝐶⁽¹⁾_𝑝𝑞,ℎ₁𝐶⁽¹⁾_𝑝𝑞2,𝐶⁽¹⁾_𝑝𝑞𝑧²₁−𝑧²₂2ℎ₁,𝐶⁽¹⁾_𝑝𝑞ℎ₁𝑧³₁−𝑧³₂3−𝑧₂𝑧²₁−𝑧²₂2−if𝑘=1𝐶^(𝑘−1)_𝑝𝑞,ℎ_𝑘−1𝐶^(𝑘−1)_𝑝𝑞2,𝐶^(𝑘−1)_𝑝𝑞𝑧²_𝑘−𝑧²_𝑘−12ℎ_𝑘−1,𝐶^(𝑘−1)_𝑝𝑞ℎ_𝑘−1𝑧³_𝑘−𝑧³_𝑘−13−𝑧_𝑘−1𝑧²_𝑘−𝑧²_𝑘−12if𝑘=𝑁+1𝐶^(𝑘−1)_𝑝𝑞−𝐶^(𝑘)_𝑝𝑞,ℎ_𝑘−1𝐶^(𝑘−1)_𝑝𝑞2+ℎ_𝑘𝐶^(𝑘)_𝑝𝑞2,𝐶^(𝑘−1)_𝑝𝑞𝑧²_𝑘−𝑧²_𝑘−12ℎ_𝑘−1+𝐶^(𝑘)_𝑝𝑞𝑧²_𝑘−𝑧²_𝑘−12ℎ_𝑘,𝐶^(𝑘−1)_𝑝𝑞ℎ_𝑘−1𝑧³_𝑘−𝑧³_𝑘−13−𝑧_𝑘−1𝑧²_𝑘−𝑧²_𝑘−12+𝐶^(𝑘)_𝑝𝑞ℎ_𝑘𝑧³_𝑘−𝑧³_𝑘+13−𝑧_𝑘−1𝑧²_𝑘−𝑧²_𝑘+12if1<𝑘<𝑁+1. (A.2)Also 𝐴_𝑝𝑞,𝐵_𝑝𝑞,𝐷_𝑝𝑞=^𝑁_𝑖=1𝐶^(𝑖)_𝑝𝑞𝑧_𝑖+1−𝑧_𝑖,𝑧²_𝑖+1−𝑧²_𝑖2,𝑧³_𝑖+1−𝑧³_𝑖3. (A.3)

The coefficient matrices [𝑀],[𝐾], and [𝐿] appearing in (25) are given as: [𝑀]=𝛿_𝑆𝐷₆₆𝛿_𝑆𝐷₂₆𝛿_𝑆𝐵₃₆−𝐵₄₅^𝑇𝛿_𝑆𝐷₂₆𝐷₂₂𝐵₂₃−𝐵₄₄^𝑇[0][0]𝐷₄₄,[𝐾]=−𝛿_𝑆𝐴₅₅+[𝛼]−𝛿_𝑆𝐴₄₅[0]−𝛿_𝑆𝐴₄₅−𝐴₄₄+[𝛼][0]𝛿_𝑆𝐵₄₅−𝐵₃₆^𝑇𝐵₄₄−𝐵₂₃^𝑇−𝐴₃₃+[𝛼],[𝐿]=𝛿_𝐴𝐴₅₅00}𝛿{}{_𝐴𝐴₄₅{0}{0}−𝛿_𝐴𝐵₄₅+𝐵₃₆𝐴₁₃𝛿_𝐶𝐵₁₃, (A.4) where [𝐴_𝑝𝑞], [𝐵_𝑝𝑞], and [𝐷_𝑝𝑞] are (𝑁+1)×(𝑁+1) square matrices containing 𝐴^𝑘𝑗_𝑝𝑞, 𝐵^𝑘𝑗_𝑝𝑞, and 𝐷^𝑘𝑗_𝑝𝑞, respectively, and the vectors {𝐴_𝑝𝑞}, {𝐵_𝑝𝑞}, and {𝐵_𝑝𝑞} are (𝑁+1)×1 column matrices containing 𝐴^𝑘_𝑝𝑞, 𝐵^𝑘_𝑝𝑞, and 𝐵^𝑘_𝑝𝑞, respectively. Also, [0] is (𝑁+1)×(𝑁+1) square zero and {0} is a zero vector with 𝑁+1 rows. The artificial matrix [𝛼] is also an (𝑁+1)×(𝑁+1) square matrix whose elements are given by 𝛼^𝑘𝑗=𝛼^ℎ/2_−ℎ/2𝜙_𝑘𝜙_𝑗d𝑧 (A.5) with 𝛼 being a relatively small parameter in comparison with the rigidity constants 𝐴^𝑘𝑗_𝑝𝑞(𝑝𝑞=33,44,55). It is noted that the insertion of [𝛼] in the matrix [𝐾] makes the eigenvalues of the matrix (−[𝑀]⁻¹[𝐾]) distinct.