Analysis of Free Edge Stresses in Composite Laminates Using Higher Order Theories

Abstract

This paper presents the determination of the interlaminar stresses close to the free edges of general cross-ply composite laminates based on higher order equivalent single-layer theory (HESL). The laminates with finite dimensions were subjected to a bending moment, an axial force, and/or a torque for investigation. Full three-dimensional stresses in the interior and the boundary-layer regions were determined. The computed results were compared with those obtained from Reddy’s layerwise theory. It was found that HESL theory predicts precisely the interlaminar stresses near the free edges of laminates. Besides, high efficiency in terms of computational time is obtainable when HESL theory is used as compared with layerwise theory. Finally, various numerical results were presented for the cross-ply laminates. Also design guidelines were proposed to minimize the edge-effect problems in composite laminates.

1. Introduction

Laminated composite materials are being used in several industries due to their high strength-to-weight ratio and stiffness-to-weight ratio. However, they are susceptible to different types of damage such as delamination which occurs due to high stress concentration near the edge of composite laminates. These stresses are induced by mismatch in elastic properties between adjacent plies of composite laminates [1]. It has been shown that the state of stresses in the edge zone of the laminate is three-dimensional (3D) in nature. Many attempts have been made to compute these stresses next to laminate’s free edges [1–12]. However, because of intrinsic complexities involved in the problem, no exact solution is known for elasticity equations. Therefore, many approximate methods used to determine the interlaminar stresses are recorded in the survey paper by Kant and Swaminathan [3]. Based on a laminated model containing anisotropic layers, the first approximate solution of interlaminar shear stresses was presented by Puppo and Evensen [4]. Other approximate analytical methods used to examine the problem are the employment of the higher order plate theory proposed by Pagano [2], the perturbation technique by Hsu and Herakovich [5], the boundary-layer theory by Tang and Levy [6], and the approximate elasticity solutions by Pipes and Pagano [7]. An approximate theory is also utilized by Pagano [8, 9] based on assumed in-plane stresses and the use of Reissner’s variational principle. The principle of minimum complementary energy and the force balance method are used by Kassapoglou and Lagace [10] to study the symmetric laminates under uniaxial loading. A variational method involving Lekhnitskii’s stress function is utilized by Yin [11, 12] to determine the interlaminar stresses in a multilayer strip of a laminate subjected to combinations of mechanical loads. Lin et al. [13] improved the technique developed by Kassapoglou and Lagace [10] for symmetric laminates under uniaxial tension to evaluate the interlaminar stress distribution near the straight free edges of symmetric and unsymmetric laminates under different types of loading conditions. The first numerical method to solve the 2D governing elasticity equations is given by Pipes and Pagano [1]. They utilized a finite-difference technique to establish the interlaminar stresses in a long symmetric laminate under uniform axial strain. A layer reduction technique and a layerwise theory (LWT) are employed by Lee and Chen [14] in the analysis of a reduced laminate. They neglected the through-thickness stretching and solved a simply supported plate subjected to bidirectional sinusoidal distributed loading. A displacement-based variable kinematic global-local finite element method is offered by Robbins and Reddy [15]. Their displacement field hierarchy contains both a conventional plate expansion (2D) and a full layerwise (3D) expansion. Neves et al. [16, 17] developed a higher order theory that considers deformations in the thickness direction under Carrera’s unified formulation to predict the buckling behaviour of laminated plates and modeling functionally graded plates accounting for extensibility in the thickness direction. The obtained governing equations and boundary conditions are then interpolated by collocation with radial basis functions. Mantari et al. [18] developed a new shear deformation theory for sandwich and composite plates. The presented theory is relatively close to 3D elasticity bending solutions. The theory accounts for adequate distribution of the transverse shear strains through the plate thickness and tangential stress-free boundary conditions on the plate boundary surface; thus, a shear correction factor is not required. A new model is proposed by Rahmani et al. [19] based on the high-order sandwich panel theory to study the effect of external loads on the free vibration of circular cylindrical composite sandwich shells with transversely compliant core, including also the calculation of the buckling loads. Thai et al. [20] present a novel finite element formulation for static, free vibration and buckling analyses of laminated composite plates. The higher order shear deformation plate theory (HSDT) is introduced in the present method to remove the shear correction factors and improve the accuracy of transverse shear stresses. A new improved high-order theory is presented by Kheirikhah et al. [21] for biaxial buckling analysis of sandwich plates with soft orthotropic core. Third-order plate theory is used for face sheets, quadratic and cubic functions which are assumed for transverse and in-plane displacements of the core, respectively.

It was found from the literature that no work has been found to study interlaminar stresses by higher order equivalent single-layer theory (HESL). Here, HESL theory is used to analytically study the interlaminar stresses in both finite and long cross-ply laminates subjected to a bending moment, an axial force, and/or a torque. Then HESL results are compared with those calculated from LWT. The presented works (HESL) are either accurate enough (and general) or computationally efficient.

2. Higher Order Equivalent Single-Layer Theory (HESL)

General cross-ply laminates are subjected to the bending moment, the axial force, and/or the torque in order to accurately determine the interlaminar stresses. The geometry of the laminate is illustrated in Figure 1. The formulation is limited to linear elastic material behavior, small strain, and displacements. The coordinate system (𝑥,𝑦,𝑧) is located at the middle plane of the laminate. Thickness, width, and length of the laminate are ℎ,2𝑏, and 2𝑎, respectively.

Figure 1 

Laminate geometry and coordinate system.

2.1. Displacement Field and Strains

The integrations of the three-dimensional elasticity strain-displacement relations [22] within the 𝑘th layer of the laminate producing the most general form of displacement field are given by𝑢₁^(𝑘)𝑥,𝑦,𝑧=𝐵₂𝑥+𝐵₆𝑥𝑧+𝑢^(𝑘)𝑦,𝑧,𝑢₂^(𝑘)𝑥,𝑦,𝑧=−𝐵₁𝑥𝑧+𝑣^(𝑘)𝑦,𝑧,𝑢₃^(𝑘)𝑥,𝑦,𝑧=𝐵₁1𝑥𝑦−2𝐵₆𝑥²+𝑤^(𝑘)𝑦,𝑧, (1a)(1b)(1c)where 𝑢₁, 𝑢₂, and 𝑢₃ represent the displacement components in the 𝑥-, 𝑦-, and 𝑧-directions, respectively, of a material point initially located at (𝑥,𝑦,𝑧) in the undeformed laminate. The displacement field in (1a), (1b), and (1c) may be used, in principle, for obtaining the stress field in any composite laminate subjected to arbitrary combinations of self-equilibrating mechanical and uniform hygrothermal loads. In the present work, however, our attention is focused on symmetric and unsymmetric cross-ply laminates under the bending moment, the axial force, and/or the torque. General cross-ply laminates based on physical grounds can be established as (see Figure 1) 𝑢₁^(𝑘)𝑥,𝑦,𝑧=−𝑢₁^(𝑘)−𝑥,𝑦,𝑧. (2) Upon imposing this condition on (1a), (1b), and (1c), it is readily seen that 𝑢^(𝑘)(𝑦,𝑧)=0. Thus, for cross-ply laminates the most general form of the displacement field is given as 𝑢₁^(𝑘)𝑥,𝑦,𝑧=𝐵₆𝑥𝑧+𝐵₂𝑢𝑥,_2(𝑘)𝑥,𝑦,𝑧=−𝐵₁𝑥𝑧+𝑣^(𝑘),𝑢𝑦,𝑧_3(𝑘)𝑥,𝑦,𝑧=𝐵₁1𝑥𝑦−2𝐵₆𝑥²+𝑤^(𝑘)𝑦,𝑧. (3) The unknown constants, namely, 𝐵₁, 𝐵₂, and 𝐵₆ appearing in (3), are global response of the laminate. On the other hand, unknown functions 𝑣^(𝑘)(𝑦,𝑧) and 𝑤^(𝑘)(𝑦,𝑧) are local response of the laminate.

The purpose of the present section is to show the HESL theory with infinity number of terms which provides sufficiently accurate results for the interlaminar stresses in composite laminate. In HESL, the components of the displacement vector at a material point in cross-ply laminated composite are expressed as [23] 𝑢₁𝑥,𝑦,𝑧=𝑢𝑥,𝑦+𝑧𝜂₁𝑥,𝑦+𝑧²𝜂₂𝑢𝑥,𝑦+⋯₂𝑥,𝑦,𝑧=𝑣𝑥,𝑦+𝑧𝜓₁𝑥,𝑦+𝑧²𝜓₂𝑢𝑥,𝑦+⋯₃𝑥,𝑦,𝑧=𝑤𝑥,𝑦+𝑧𝜙₁𝑥,𝑦+𝑧²𝜙₂𝑥,𝑦+⋯. (4) Relations (4) may more conveniently be presented as 𝑢₁𝑥,𝑦,𝑧=𝑢𝑥,𝑦+𝑧^𝑘𝜂_𝑘𝑢𝑥,𝑦𝑘=1,2,…,𝑙,₂𝑥,𝑦,𝑧=𝑣𝑥,𝑦+𝑧^𝑖𝜓_𝑖𝑢𝑥,𝑦𝑖=1,2,…,𝑛,₃𝑥,𝑦,𝑧=𝑤𝑥,𝑦+𝑧^𝑗𝜙_𝑗𝑥,𝑦𝑗=1,2,…,𝑚, (5) where 𝑘, 𝑖, and 𝑗 are dummy indexes indicating the summation of terms from 𝑘=𝑖=𝑗=1 to desirable number (i.e., 𝑙 or 𝑛 or 𝑚). It is expected that the accuracy of such theories can, generally, be increased by taking more terms in (5). Also, 𝑢₁, 𝑢₂, and 𝑢₃ are the displacement components of any point within a laminate in the 𝑥-, 𝑦-, and 𝑧-directions, respectively. 𝑘+𝑖+𝑗+3 unknown displacement functions should be found for a single layer or more layers in the laminate. The displacement field in (3) within HESL is simplified by 𝑢₁𝑥,𝑦,𝑧=𝐵₆𝑥𝑧+𝐵₂𝑢𝑥,₂𝑥,𝑦,𝑧=−𝐵₁𝑥𝑧+𝑣𝑦+𝑧^𝑖𝜓_𝑖𝑦𝑢𝑖=1,2,…,𝑛,₃𝑥,𝑦,𝑧=𝐵₁1𝑥𝑦−2𝐵₆𝑥²+𝑤𝑦+𝑧^𝑗𝜙_𝑗𝑦𝑗=1,2,…,𝑚. (6) By using the displacement field (6) in the principle of minimum total potential energy [22] and considering the constants as unknown parameters, (5+𝑖+𝑗) equilibrium equations are expressed as𝛿𝑣∶𝑁^′_𝑦=0,𝛿𝑤∶𝑁^′_𝑦𝑧=0,𝛿𝜓_𝑖∶𝐵_𝑖𝑦𝑧−𝑑𝐴_𝑖𝑦𝑑𝑦=0𝑖=1,2,…,𝑛,𝛿𝜙_𝑗∶𝐵_𝑗𝑧−𝑑𝐴_𝑗𝑦𝑧𝑑𝑦=0𝑗=1,2,…,𝑚, (7a)(7b)(7c)(7d)𝛿𝐵₁∶^𝑏_−𝑏𝑄_𝑥𝑦−𝑀_𝑥𝑦𝑑𝑦=𝑇₀,𝛿𝐵₂∶^𝑏_−𝑏𝑁_𝑥𝑑𝑦=𝐹₀,𝛿𝐵₆∶^𝑏_−𝑏𝑀_𝑥𝑑𝑦=𝑀₀, (8a)(8b)(8c)where a prime indicates an ordinary differentiation with respect to variable 𝑦 and the stress and moment resultants appearing in (7a), (7b), (7c), and (7d) and (8a), (8b), and (8c) which are defined as 𝑀_𝑥,𝑀_𝑥𝑦=^ℎ/2_−ℎ/2𝜎_𝑥,𝜎_𝑥𝑦𝑁𝑧𝑑𝑧,_𝑥,𝑁_𝑦,𝑁_𝑦𝑧,𝑄_𝑥=^ℎ/2_−ℎ/2𝜎_𝑥,𝜎_𝑦,𝜎_𝑦𝑧,𝜎_𝑥𝑧𝐴𝑑𝑧,_𝑖𝑦,𝐴_𝑗𝑦𝑧=^ℎ/2_−ℎ/2𝜎_𝑦𝑧^𝑖,𝜎_𝑦𝑧𝑧^𝑗𝐵𝑑𝑧,_𝑗𝑧,𝐵_𝑖𝑦𝑧=^ℎ/2_−ℎ/2𝜎_𝑧𝑗𝑧^𝑗−1,𝜎_𝑦𝑧𝑖𝑧^𝑖−1𝑑𝑧. (9) By substituting the displacement field (6) into (9), through the linear strain-displacement relations of elasticity and the plane-stress constitutive law [24] of a lamina, the stress and moment resultants can be expressed as𝑀_𝑥𝑦,𝑀_𝑥,𝑁_𝑥,𝑁_𝑦,𝐴_𝑖𝑦,𝐵_𝑗𝑧=−𝐷₆₆,𝐷₁₆,𝐵₁₆,𝐵₂₆,0,0𝐵₁+𝐷₁₆,𝐷₁₁,𝐵₁₁,𝐵₁₂,𝐹_𝑖+112,𝑗𝐻_𝑗13𝐵₆+𝐵₁₂,𝐵₁₁,𝐴₁₁,𝐴₁₂,𝐹_𝑖12,𝑗𝐻_𝑗−113𝐵₂+𝐵₂₆,𝐵₁₂,𝐴₁₂,𝐴₂₂,𝐹_𝑖22,𝑗𝐻_𝑗−123𝑉^′+0,𝐹_𝑖+112,𝐹_𝑖12,𝐹_𝑖22,𝐹_2𝑖22,𝑗𝐿_{𝑖𝑗−123}𝜓^′_𝑖+0,𝑗𝐻_𝑗13,𝑗𝐻_𝑗−113,𝑗𝐻_𝑗−123,𝑗𝐿_{𝑖𝑗−123},𝑗²𝐻_2𝑗−233𝜙_𝑗,𝑁_𝑦𝑧,𝐴_𝑗𝑦𝑧,𝐵_𝑖𝑦𝑧=𝑖𝐹_𝑖−144,𝑖𝐿_{𝑖−1𝑗44},𝑖²𝐹_2𝑖−244𝜓_𝑖+𝐴₄₄,𝐻_𝑗44,𝑖𝐹_𝑖−144𝑤^′+𝐻_𝑗44,𝐻_2𝑗44,𝑖𝐿_{𝑖−1𝑗44}𝜙^′_𝑗,𝑄_𝑥=𝐴₅₅𝐵₁𝑦, (10a)(10b)(10c)where 𝐴_𝑝𝑞,𝐵_𝑝𝑞,𝐷_𝑝𝑞=^𝑁_𝑖=1^𝑧𝑖+1_𝑧𝑖𝐶^(𝑘)_𝑝𝑞1,𝑧,𝑧²𝐹𝑑𝑧𝑝,𝑞=1,2,6,_𝑖𝑝𝑞=^𝑁_𝑖=1^𝑧𝑖+1_𝑧𝑖𝐶^(𝑘)_𝑝𝑞𝑧^𝑖𝐻𝑑𝑧,_𝑗𝑝𝑞=^𝑁_𝑖=1^𝑧𝑖+1_𝑧𝑖𝐶^(𝑘)_𝑝𝑞𝑧^𝑗𝐿𝑑𝑧,_{𝑖𝑗𝑝𝑞}=^𝑁_𝑖=1^𝑧𝑖+1_𝑧𝑖𝐶^(𝑘)_𝑝𝑞𝑧^𝑖𝑧^𝑗𝑑𝑧, (11) where 𝐶_𝑝𝑞 are the transformed stiffness of an orthotropic lamina. For free edges of the laminate at 𝑦=±𝑏 according to the principle of minimum total potential energy the traction-free boundary conditions must be imposed as 𝑀_𝑥𝑦=𝑁_𝑦=𝑁_𝑦𝑧=𝐴_𝑖𝑦=𝐴_𝑗𝑦𝑧=0at𝑦=±𝑏. (12) Substituting (10a), (10b), (10c) into (7a), (7b), (7c), the equilibrium equations are obtained in terms of the displacement components which can be given as𝛿𝑣∶𝐵₁₂𝐵₆+𝐴₁₂𝐵₂−𝐵₂₆𝐵₁+𝐴₂₂𝑣^′+𝐹_𝑖22𝜓^′_𝑖+𝑗𝐻_𝑗−123𝜙_𝑗=0,𝛿𝑤∶𝐴₄₄𝑤^′+𝑖𝐹_𝑖−144𝜓_𝑖+𝐻_𝑗44𝜙^′_𝑗=0,𝛿𝜓_𝑖∶𝑖²𝐹_2𝑖−244𝜓_𝑖+𝑖𝐹_𝑛−144𝑤^′+𝑖𝐿_{𝑖−1𝑗44}𝜙^′_𝑗−𝐹_𝑖22𝑣^′′−𝐹_2𝑖22𝜓_𝑖^′′−𝑗𝐿_{𝑖𝑗−123}𝜙^′_𝑗=0,𝛿𝜙_𝑗∶𝑗𝐻_𝑗13𝐵₆+𝑗𝐻_𝑗−113𝐵₂+𝑗𝐻_𝑗−123𝑣^′+𝑗𝐿_{𝑖𝑗−123}𝜓^′_𝑖+𝑗²𝐻_2𝑗−233𝜙_𝑗−𝑖𝐿_{𝑖−1𝑗44}𝜓^′_𝑖−𝐻_𝑗44𝑤^′′−𝐻_2𝑗44𝜙_𝑗^′′=0. (13a)(13b)(13c)(13d)To determine the parameters 𝐵₁, 𝐵₂, and 𝐵₆ as well as the interlaminar stresses in (13a)–(13d), the general solution of the ordinary differential equations in (13a)–(13d) is first obtained in terms of 𝐵₁, 𝐵₂, and 𝐵₆. Then, by applying the boundary conditions in (12) and using the global equilibrium equation in (8a), (8b), and (8c), the constant parameters 𝐵₁, 𝐵₂, and 𝐵₆ will be found in terms of the bending moment 𝑀₀, axial force 𝐹₀, and torque 𝑇₀. For completeness, the details of the steps involved are displayed in Appendix A.

3. Layerwise Theory

3.1. Displacement Field and Strains

It should be emphasized that the solution proposed by Tahani and Nosier [25] was used to verify HESL theory for predicting stress in the edge of composite laminates. The proposed solution was used while different boundary conditions are considered. This offers to develop new analytical solution which is expressed in the present work. The displacement field in this theory may be represented as [25] 𝑢₁𝑥,𝑦,𝑧=𝑈_𝑘𝑥,𝑦Φ_𝑘(𝑧),𝑢₂𝑥,𝑦,𝑧=𝑉_𝑘𝑥,𝑦Φ_𝑘(𝑧),𝑢₃𝑥,𝑦,𝑧=𝑊_𝑘𝑥,𝑦Φ_𝑘(𝑧),(𝑘=1,2,…,𝑁+1), (14) where 𝑢₁, 𝑢₂, and 𝑢₃ represent the displacement components in the 𝑥-, 𝑦-, and 𝑧-directions, respectively [26]. Also 𝑈_𝑘(𝑥,𝑦), 𝑉_𝑘(𝑥,𝑦), and 𝑊_𝑘(𝑥,𝑦) show the displacement components of all points located on the 𝑘th plane in the undeformed laminate, and Φ_𝑘(𝑧) is continuous function of the thickness coordinate 𝑧. Moreover, 𝑁 indicates the total number of numerical layers considered in a laminate. It should be noted that a repeated index indicates summation over all values of that index. Substituting (14) into the linear strain-displacement relations of elasticity [27], the results are obtained as 𝜀_𝑥=𝜕𝑈_𝑘Φ𝜕𝑥^𝑘,𝛾_𝑦𝑧=𝑉_𝑘𝑑Φ_𝑘+𝑑𝑧𝜕𝑊_𝑘Φ𝜕𝑦^𝑘,𝜀_𝑦=𝜕𝑉_𝑘Φ𝜕𝑦^𝑘,𝛾_𝑥𝑧=𝑈_𝑘𝑑Φ_𝑘+𝑑𝑧𝜕𝑊_𝑘Φ𝜕𝑥^𝑘,𝜀_𝑧=𝑑Φ_𝑘𝑊𝑑𝑧^𝑘,𝛾_𝑥𝑦=𝜕𝑈_𝑘Φ𝜕𝑦^𝑘+𝜕𝑉_𝑘Φ𝜕𝑥^𝑘. (15)

3.2. Equilibrium Equations

Using the principle of minimum total potential energy [23], 3(𝑁+1) equations of equilibrium corresponding to 3(𝑁+1) unknowns 𝑈_𝑘, 𝑉_𝑘, and 𝑊_𝑘 can be shown by 𝜕𝑀^𝑘_𝑥+𝜕𝑥𝜕𝑀^𝑘_𝑥𝑦𝜕𝑦−𝑄^𝑘_𝑥=0,𝜕𝑀^𝑘_𝑥𝑦+𝜕𝑥𝜕𝑀^𝑘_𝑦𝜕𝑦−𝑄^𝑘_𝑦=0,𝜕𝑅^𝑘_𝑥+𝜕𝑥𝜕𝑅^𝑘_𝑦𝜕𝑦−𝑁^𝑘_𝑧=0, (16) where 𝑘=1,2,…,𝑁+1. In (16), the generalized stress and moment resultants are defined as 𝑁^𝑘_𝑧=^ℎ/2_−ℎ/2𝜎_𝑧𝑑Φ_𝑘𝑅𝑑𝑧𝑑𝑧,_𝑘𝑥,𝑅^𝑘_𝑦=^ℎ/2_−ℎ/2𝜎_𝑥𝑧,𝜎_𝑦𝑧Φ_𝑘𝑄𝑑𝑧,_𝑘𝑥,𝑄^𝑘_𝑦=^ℎ/2_−ℎ/2𝜎_𝑥𝑧,𝜎_𝑦𝑧𝑑Φ_𝑘𝑀𝑑𝑧𝑑𝑧,_𝑘𝑥,𝑀^𝑘_𝑦,𝑀^𝑘_𝑥𝑦=^ℎ/2_−ℎ/2𝜎_𝑥,𝜎_𝑦,𝜎_𝑥𝑦Φ_𝑘𝑑𝑧. (17) The boundary conditions for a laminated plate with a rectangular platform in the layerwise theory at an edge parallel to 𝑦-axis involve the specification of 𝑈_𝑘 or 𝑀^𝑘_𝑥, 𝑉_𝑘 or 𝑀^𝑘_𝑥𝑦, and 𝑊_𝑘 or 𝑅^𝑘_𝑥. Similarly, at an edge parallel to 𝑥-axis, the boundary conditions involve the specification of 𝑈_𝑘 or 𝑀^𝑘_𝑥𝑦, 𝑉_𝑘 or 𝑀^𝑘_𝑦, and 𝑊_𝑘 or 𝑅^𝑘_𝑦.

3.3. Displacement Equations of Equilibrium

The linear constitutive relations for the 𝑘th orthotropic lamina, with fiber orientations of 0° and 90° only, are given by [24] 𝜎_𝑥𝜎_𝑦𝜎_𝑧𝜎_𝑦𝑧𝜎_𝑥𝑧𝜎_𝑥𝑦^(𝑘)=𝐶₁₁𝐶₁₂𝐶₁₃000𝐶₂₂𝐶₂₃000𝐶₃₃000𝐶₄₄00𝐶₅₅0Sym.𝐶₆₆^(𝑘)×𝜕𝑈_𝑘Φ𝜕𝑥^𝑘𝜕𝑉_𝑘Φ𝜕𝑦^𝑘𝑑Φ_𝑘𝑊𝑑𝑧^𝑘𝑉_𝑘𝑑Φ_𝑘+𝑑𝑧𝜕𝑊_𝑘Φ𝜕𝑦^𝑘𝑈_𝑘𝑑Φ_𝑘+𝑑𝑧𝜕𝑊_𝑘Φ𝜕𝑥^𝑘𝜕𝑈_𝑘Φ𝜕𝑦^𝑘+𝜕𝑉_𝑘Φ𝜕𝑥^𝑘, (18) where 𝐶^(𝑘)_𝑖𝑗 represent the off-axis stiffnesses within the 𝑘th layer of the laminate. By substituting (18) into (17), the stress and moment resultants are written as 𝑀^𝑘_𝑥,𝑀^𝑘_𝑦=𝐷^𝑘𝑗₁₁,𝐷^𝑘𝑗₁₂𝜕𝑈_𝑗+𝐷𝜕𝑥^𝑘𝑗12,𝐷^𝑘𝑗₂₂𝜕𝑉_𝑗+𝐵𝜕𝑦_𝑘𝑗13,𝐵^𝑘𝑗₂₃𝑊_𝑗,𝑁^𝑘_𝑧=𝐵^𝑗𝑘₁₃𝜕𝑈_𝑗𝜕𝑥+𝐵^𝑗𝑘₂₃𝜕𝑉_𝑗𝜕𝑦+𝐴^𝑘𝑗₃₃𝑊_𝑗,𝑀^𝑘_𝑥𝑦=𝐷^𝑘𝑗₆₆𝜕𝑈_𝑗𝜕𝑦+𝐷^𝑘𝑗₆₆𝜕𝑉_𝑗,𝑄𝜕𝑥_𝑘𝑦,𝑅^𝑘_𝑦=𝐴^𝑘𝑗₄₄,𝐵^𝑘𝑗₄₄𝑉_𝑗+𝐵^𝑗𝑘₄₄,𝐷^𝑘𝑗₄₄𝜕𝑊_𝑗,𝑄𝜕𝑦_𝑘𝑥,𝑅^𝑘_𝑥=𝐴^𝑘𝑗₅₅,𝐵^𝑘𝑗₅₅𝑈_𝑗+𝐵^𝑗𝑘₅₅,𝐷^𝑘𝑗₅₅𝜕𝑊_𝑗,𝜕𝑥 (19) where the coefficients 𝐴^𝑘𝑗, 𝐵^𝑘𝑗, and 𝐷^𝑘𝑗 appearing in (19) are defined in Appendix B. Finally, substituting (19) into (16) yields the governing equations of equilibrium as 𝐷^𝑘𝑗₁₁𝑈_{𝑗,𝑥𝑥}+𝐷^𝑘𝑗₆₆𝑈_{𝑗,𝑦𝑦}−𝐴^𝑘𝑗₅₅𝑈_𝑗+𝐷^𝑘𝑗₁₂+𝐷^𝑘𝑗₆₆𝑉_{𝑗,𝑥𝑦}+𝐵^𝑘𝑗₁₃−𝐵^𝑗𝑘₅₅𝑊_𝑗,𝑥𝐷=0,_𝑘𝑗12+𝐷^𝑘𝑗₆₆𝑈_{𝑗,𝑥𝑦}+𝐷^𝑘𝑗₆₆𝑉_{𝑗,𝑥𝑥}+𝐷^𝑘𝑗₂₂𝑉_{𝑗,𝑦𝑦}−𝐴^𝑘𝑗₄₄𝑉_𝑗+𝐵^𝑘𝑗₂₃−𝐵^𝑗𝑘₄₄𝑊_𝑗,𝑦𝐵=0,_𝑘𝑗55−𝐵^𝑗𝑘₁₃𝑈_𝑗,𝑥+𝐵^𝑘𝑗₄₄−𝐵^𝑗𝑘₂₃𝑉_𝑗,𝑦+𝐷^𝑘𝑗₅₅𝑊_{𝑗,𝑥𝑥}+𝐷^𝑘𝑗₄₄𝑊_{𝑗,𝑦𝑦}−𝐴^𝑘𝑗₃₃𝑊_𝑗=0, (20) where a comma followed by a variable indicates differentiation with respect to that variable.

For completeness, the details of the steps involved in the analytical solutions are displayed in Appendix C.

4. Results

4.1. Numerical Results and Discussion

To verify the accuracy and efficiency of the present method, several numerical examples are presented for general cross-ply laminates subjected to the bending moment, axial force, and/or torque. The analyses were performed for 4-layer and 6-layer laminates made of graphite/epoxy. The used mechanical and physical properties of the layers are presented in Table 1 [24].

Table 1 

Engineering properties of graphite/epoxy.

In addition, the thickness of each physical ply is assumed to be 0.5 mm (i.e., ℎ_𝑘=0.5mm). HESL (𝑖) represents that 𝑖 number of (𝑗=𝑖) term is taken in (6). Clearly, as the number 𝑖 is increased, the accuracy of the results is also increased. In the numerical examples that follow the interlaminar stresses are determined by integrating the local equilibrium equations of elasticity. Also, the width-to-thickness ratio (i.e., 2𝑏/ℎ) is assumed to be, unless otherwise mentioned, equal to 10. Furthermore, the stress components are normalized as 𝜎_𝑖𝑗=𝜎_𝑖𝑗𝜎₀, (21) where 𝜎₀=(1/𝑏ℎ)(𝐹₀+𝑀₀/ℎ+𝑇₀/ℎ).

To study the convergence of the stresses near free edges, two symmetric laminates [90^∘/0^∘]_𝑠 and [0^∘/90^∘]_𝑠 subjected to the bending moment, the axial force, and the torque are considered. On account of the nature of the cross-ply laminates, both 𝜎_𝑥𝑦 and 𝜎_𝑥𝑧 are zero here as expected. Figures 2 and 3 demonstrate the convergence of the solution for 𝜎_𝑧 versus 𝑦 at the 𝑧=ℎ/4 for the [90^∘/0^∘]_𝑠 and [0^∘/90^∘]_𝑠 laminates. It is seen that, for both [90^∘/0^∘]_𝑠 and [0^∘/90^∘]_𝑠, 𝜎_𝑧 is seen to rise or fall suddenly near the free edge, while being zero in the inner region of the laminate also the peak stress steadily increases as the number of terms (𝑖) taken is increased. This is often attributed to a possible singularity at the 0°/90° interfaces.

Figure 2 

Convergence of interlaminar normal stress along the 0°/90° interfaces of [90∘/0∘/0∘/90∘] laminate.

Figure 3 

Convergence of interlaminar normal stress along the 0°/90° interfaces of [0∘/90∘/90∘/0∘] laminate.

Numerical results are obtained from equilibrium equations. The present results are directly compared with those obtained from LWT. Figure 4 shows the 𝜎_𝑧 distribution at the 0°/90° interfaces of [0°/90°/90°/0°] and [0°/90°/0°/90°] laminates under the bending moment and the axial force. An excellent agreement is found between the present solutions and those of LWT.

Figure 4 

Distribution of interlaminar normal stress along the interfaces of [90∘/0∘/0∘/90∘] and [0∘/90∘/90∘/0∘] laminates.

Figures 5 and 6 show the 𝜎_𝑧 and 𝜎_𝑦𝑧 distribution at the 0°/90° interfaces of [0°/90°/90°/90°/90°/0°] and [0°/0°/0°/90°] laminates under the bending moment and the torque. Good accordance is seen between the results of the two theories. It is noted that the accuracy of HESL theory can be improved by taking more terms. Also the layerwise theory needs to take many terms to approach accurately the results of the present theory and these terms cause that LWT is more complex and computationally more time consuming than the HESL theory, so, it is better that HESL theory is used to compute the local phenomena such as free-edge-effect problems and the distribution of interlaminar stresses more precisely and less computationally than the layerwise theory in laminate composites.

Figure 5 

Distribution of interlaminar stress along the 0∘/90∘ interface of [0∘/90∘/90∘/90∘/90∘/0∘] laminate.

Figure 6 

Distribution of interlaminar stress along the 0∘/90∘ interface of [0∘/0∘/0∘/90∘] laminate.

The effect of the laminate width-to-thickness ratio on the interlaminar stresses for the [90°/0°/90°/0°/90°/0°] laminate under the bending moment is investigated in Figure 7. By decreasing the width-to-thickness ratio, the boundary-layer region is 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 does not change while the width-to-thickness ratio of the laminate changes. On the other hand, this is evident that the highly localized nature of interlaminar stresses occurs near and exactly at the free edges of the laminate.

Figure 7 

Distribution of interlaminar normal stress along the first 90∘/0∘ of [90∘/0∘/90∘/0∘/90∘/0∘] laminate for various width-to-thickness ratios.

The variation of interlaminar shear stress 𝜎_𝑦𝑧 through the thickness and near the free edge of the [90°/0°/90°/0°] laminate under the bending moment and the torque is revealed in Figure 8. It is observed that the maximum negative and the maximum positive values of 𝜎_𝑦𝑧 occur within the bottom 90° layer and the top 0° layer at the 90°/0° interfaces of the laminate, respectively. The variations of interlaminar normal stress 𝜎_𝑧 through the thickness of the unsymmetrical cross-ply laminate [90°/0°/90°/0°] under the bending moment and the axial force are portrayed in Figure 9. The maximum negative and the maximum positive values of 𝜎_𝑧 occur within the bottom 90° layer and the top 0° layer both near the 90°/0° interfaces at the free edge (i.e., 𝑦=𝑏), respectively. From Figure 9, it is seen that 𝜎_𝑧 decreases away from the free edge as the inner region of the laminate is approached.

Figure 8 

Distribution of interlaminar shear stress through thickness of [90∘/0∘/90∘/0] laminate.

Figure 9 

Distribution of interlaminar normal stress through thickness of [90∘/0∘/90∘/0] laminate.

The distributions of the interlaminar stresses 𝜎_𝑧 and 𝜎_𝑦𝑧 along the upper and lower interfaces of the unsymmetrical cross-ply [90°/90°/0°/90°] laminate subjected to the bending moment, the axial force, and the torque are demonstrated in Figure 10. It is observed that the interlaminar stresses demonstrate high stress gradient near the free edge. Both stresses are seen to grow suddenly near the free edge, while being zero in the interior region of the laminate. It is also observed that the interlaminar shear stress 𝜎_𝑦𝑧 rises toward the free edge and decreases (or increases) rather abruptly to zero at the free edge.

Figure 10 

Distribution of interlaminar stress along the width in [90∘/90∘/0∘/90∘] laminate.

4.2. Preliminary Design Guidelines

A comprehensive structural analysis program for designing composite laminates is very extensive and complex. This generally involves several analysis phases such as laminate stress and strength analysis. Hence, there is still a need to provide some preliminary knowledge of the lay-up sequence of composite laminates. Usually, structural properties of composite laminates such as stiffness, strength, and dimensional stability are affected by the laminate stacking sequences. Because each property has different relations with a particular stacking sequence, the choice of stacking sequence suited for a particular application may require a compromise. The obtained preliminary design guidelines consist of the following.(1)The laminates stacking sequence (LSS) should be symmetric about the midplane to avoid extension-bending coupling. If this is not possible due to other requirements, locate the asymmetry or imbalance as near to the laminate midplane as possible. Avoid symmetric LSS that create high interlaminar tension stresses (𝑠) at free edges.(2)Avoid grouping of 90° plies and separate 90° plies by a 0° ply to minimize interlaminar shear and normal stress. Minimize groupings of plies with the same orientations to create a more homogeneous laminate and to minimize interlaminar stress. If plies must be grouped, avoid grouping more than four plies of the same orientation. Minimizing grouping helps to increase strength and minimize interlaminar shear and normal stresses and therefore minimize the tendency to delaminate.(3)Shield primary load carrying plies by positioning them inside of laminate to increase tensile strength and buckling resistance.(4)An LSS should have at least both distinct ply angles (e.g., 0°, 90°) with a minimum of 10% of the plies oriented at each angle. Ply angles should be selected such that fibers are oriented with principal load axes.(5)Locating 90° ply toward the exterior surfaces improves the buckling allowable in many cases.

5. Conclusions

In this research, analytical solutions were established within HESL theory for the edge-effect problem of general cross-ply composite laminates with finite dimensions under the bending moment, the axial force, and/or the torque. The edges of the laminates at 𝑦=±𝑏 were assumed to have traction-free boundary conditions. The accuracy and effectiveness of HESL theory in describing the localized three-dimensional effects were demonstrated by comparing the results of HESL theory with those calculated from the layerwise theory. Good agreement was observed between the results of HESL and LWT theories. Furthermore, the analysis using HESL was found to be more cost effective and accurate, so HESL was employed to assess the local phenomena instead of LWT theory. Here, several numerical results were shown for the different loading problems. The design guidelines were developed to provide cross-ply laminate stacking sequences with minimized the free edges effects.

Appendices

A.

According to the principle of minimum total potential energy [22] at the equilibrium configuration of a body the variation of the total potential energy Π of the body must vanish as 𝛿Π≡𝛿𝑈+𝛿𝑉=0, (A.1) where 𝛿𝑈 is the variation of total strain energy of the body; that is, 𝛿𝑈=_𝑉𝜎_𝑥𝛿𝜀_𝑥+𝜎_𝑦𝛿𝜀_𝑦+𝜎_𝑧𝛿𝜀_𝑧+𝜎_𝑦𝑧𝛿𝛾_𝑦𝑧+𝜎_𝑥𝑧𝛿𝛾_𝑥𝑧+𝜎_𝑥𝑦𝛿𝛾_𝑥𝑦𝑑𝑥𝑑𝑦𝑑𝑧, (A.2) and 𝑉 is negative of the work done on the body by the specified external forces. Here, 𝑉=−2𝑀₀𝑎𝐵₆−2𝑇₀𝑎𝐵₁−2𝐹₀𝑎𝐵₂ and therefore, 𝛿𝑉=−2𝑀₀𝑎𝛿𝐵₆−2𝑇₀𝑎𝛿𝐵₁−2𝐹₀𝑎𝛿𝐵₂. Also, the variations of strains in (A.2) are found as 𝛿𝜀_𝑥=𝑧𝛿𝐵₆+𝛿𝐵₂,𝛿𝛾_𝑦𝑧=𝑖𝑧^𝑖−1𝛿𝜓_𝑖+𝛿𝑤^′+𝑧^𝑗𝛿𝜙^′_𝑗,𝛿𝜀_𝑦=𝛿𝑣^′+𝑧^𝑖𝛿𝜓^′_𝑖,𝛿𝛾_𝑥𝑧=𝐵₁𝑦,𝛿𝜀_𝑧=𝑗𝑧^𝑗−1𝛿𝜙_𝑗,𝛿𝛾_𝑥𝑦=−𝐵₁𝑧. (A.3) Upon substituting (A.3) and 𝛿𝑉=−2𝑀₀𝑎𝛿𝐵₆−2𝑇₀𝑎𝛿𝐵₁−2𝐹₀𝑎𝛿𝐵₂ into (A.1), carrying out the necessary integrations, and employing the fundamental lemma of calculus of variations the equilibrium equations and the associated boundary conditions of a laminate under bending are obtained to be as in (7a), (7b), (7c), and (7d), (8a), (8b), and (8c), and (12), respectively. To solve the linear equations in (13a), (13b), (13c), and (13d), first we obtained 𝑣^′ and 𝑤^′ from (13a) and (13b); then substituting 𝑣^′ and 𝑤^′ into (13c) and (13d) and integrating (13c) and (13d) yield 𝛿𝜓_𝑖∶𝑆₁𝑖,𝑗𝜓_𝑖^′′+𝑆₂𝑖,𝑗𝜓_𝑖+𝑆₃𝑖,𝑗𝜙_𝑗^′′=0𝑖=1,2,…,𝑛,𝛿𝜙_𝑗∶𝑆₄𝑖,𝑗𝜓_𝑖+𝑆₅𝑖,𝑗𝜙_𝑗^′′+𝑆₆𝑖,𝑗𝜙_𝑗=𝑆₇𝑗𝐵₆𝑦+𝑆₈𝑗𝐵₂𝑦+𝑆₉𝑗𝐵₁𝑦𝑗=1,2,…,𝑚,𝜙_𝑗=^𝑦𝜙_𝑗𝑑𝑦, (A.4) where 𝑆₁=𝐹𝑖,𝑗^𝑖22𝐹_𝑗22𝐴₂₂−𝐹_{0.5(𝑖+𝑗)22}𝐹𝑖+𝑗=even_𝑖22𝐹_𝑗22𝐴₂₂−𝐹_𝑗22𝐹𝑖>𝑗_𝑖22𝐹_𝑗22𝐴₂₂−𝐹_𝑖22𝑆𝑗>𝑖,₂=𝑖𝑖,𝑗²𝐹_{𝑖+𝑗−244}−𝑖𝐹_𝑖−144𝐹_𝑗−144𝐴₄₄𝑗𝑖+𝑗=even₂𝐹_2𝑗−244−𝑖𝐹_𝑖−144𝐹_𝑗−144𝐴₄₄𝑖𝑖>𝑗₂𝐹_2𝑖−244−𝑖𝐹_𝑖−144𝐹_𝑗−144𝐴₄₄𝑆𝑗>𝑖,₃𝑖,𝑗=𝑖𝐿_{𝑖−1𝑗44}−𝑖𝐹_𝑖−144𝐻_𝑗44𝐴₄₄−𝑗𝐿_{𝑖𝑗−123}+𝑗𝐹_𝑖22𝐻_𝑗−123𝐴₂₂,𝑆₄=𝑖,𝑗𝑗𝐿_{𝑖𝑗−123}+𝑖𝐹_𝑖−144𝐻_𝑗44𝐴₄₄−𝑖𝐿_{𝑖−1𝑗44}−𝑗𝐹_𝑖22𝐻_𝑗−123𝐴₂₂𝑖+𝑗=even𝑖𝐿_{𝑖𝑗−123}+𝑗𝐹_𝑗−144𝐻_𝑖44𝐴₄₄−𝑗𝐿_{𝑖−1𝑗44}−𝑗𝐹_𝑖22𝐻_𝑗−123𝐴₂₂𝑆otherwise,₅=𝐻𝑖,𝑗^𝑖44𝐻_𝑗44𝐴₄₄−𝐻_{0.5(𝑖+𝑗)44}𝐻𝑖+𝑗=even_𝑖44𝐻_𝑗44𝐴₄₄−𝐻_𝑗44𝐻𝑖>𝑗_𝑖44𝐻_𝑗44𝐴₄₄−𝐻_𝑖44𝑆𝑗>𝑖,₆=𝑗𝑖,𝑗²𝐻_{𝑖+𝑗−233}−𝑗𝐻_𝑗−123𝐻_𝑖−123𝐴₂₂𝑗𝑖+𝑗=even₂𝐻_2𝑗−233−𝑗𝐻_𝑗−123𝐻_𝑖−123𝐴₂₂𝑖𝑖>𝑗₂𝐻_2𝑖−233−𝑗𝐻_𝑗−123𝐻_𝑖−123𝐴₂₂𝑆𝑗>𝑖,₇𝑗𝐵=𝑗¹²𝐴₂₂𝐻_𝑗−123−𝑗𝐻_𝑗13,𝑆₈𝑗𝐴=𝑗¹²𝐴₂₂𝐻_𝑗−123−𝑗𝐻_𝑗−113,𝑆₉𝑗𝐵=𝑗²⁶𝐴₂₂𝐻_𝑗−123𝑗=1,2,…,𝑚𝑖=1,2,…,𝑛, (A.5) where [𝑆_𝑃](𝑃=1,2,…,6) and the vectors {𝑆_𝑄}(𝑄=7,8,9) are 𝑛×𝑚 matrices containing 𝑆_𝑃(𝑖,𝑗) and 𝑚×1 column matrices containing 𝑆_𝑄(𝑗), respectively. We introduced (A.4) in a matrix form as follows: [𝑀]𝜇^′′+[𝐾]𝜇=[𝐿]{𝐵}⋅𝑦, (A.6) where 𝜇=𝜓₁,𝜙₁,𝜓₂,𝜙₂,…,𝜓_𝑛,𝜙_𝑚^𝑇,{𝐵}=𝐵₁,𝐵₂,𝐵₆^𝑇, (A.7) with[𝑀]=𝑆₁(1,1)𝑆₃(1,1)𝑆₁(1,2)𝑆₃(1,2)⋯𝑆₁(1,𝑚)𝑆₃(1,𝑚)0𝑆₅(1,1)0𝑆₅(1,2)⋯0𝑆₅()𝑆1,𝑚₁(2,1)𝑆₃(2,1)𝑆₁(2,2)𝑆₃(2,2)⋯𝑆₁(2,𝑚)𝑆₃(2,𝑚)0𝑆₅(2,1)0𝑆₅(2,2)⋯0𝑆₅()⋯𝑆2,𝑚₁(𝑛,1)𝑆₃(𝑛,1)𝑆₁(𝑛,2)𝑆₃(𝑛,2)⋯𝑆₁(𝑛,𝑚)𝑆₃(𝑛,𝑚)0𝑆₅(𝑛,1)0𝑆₅(𝑛,2)⋯0𝑆₅(),[𝐾]=𝑆𝑛,𝑚₂(1,1)0𝑆₂(1,2)0⋯𝑆₂()0𝑆1,𝑚₄(1,1)𝑆₆(1,1)𝑆₄(1,2)𝑆₆(1,2)⋯𝑆₄(1,𝑚)𝑆₆()𝑆1,𝑚₂(2,1)0𝑆₂(2,2)0⋯𝑆₂()0𝑆2,𝑚₄(2,1)𝑆₆(2,1)𝑆₄(2,2)𝑆₆(2,2)⋯𝑆₄(2,𝑚)𝑆₆()⋯𝑆2,𝑚₂(𝑛,1)0𝑆₂(𝑛,2)0⋯𝑆₂()0𝑆𝑛,𝑚₄(𝑛,1)𝑆₆(𝑛,1)𝑆₄(𝑛,2)𝑆₆(𝑛,2)⋯𝑆₄(𝑛,𝑚)𝑆₆(),[𝐿]=𝑆𝑛,𝑚000₉(1)𝑆₈(1)𝑆₇(1)𝑆000₉(2)𝑆₈(2)𝑆₇(2)𝑆⋯⋯⋯000₉(𝑚)𝑆₈(𝑚)𝑆₇(𝑚). (A.8)The general solution of (A.6) can be readily shown by 𝜇=𝜓sinh𝜆𝑦{𝐻}+[𝐾]⁻¹[𝐿]{𝐵}⋅𝑦, (A.9) where [sinh(𝜆𝑦)] is a 2𝑛×2𝑚 diagonal matrix. Also in (A.9), [𝜓] and (𝜆²₁,𝜆²₂,…,𝜆²_𝑛+𝑚) represent the modal matrix and eigenvalues of (−[𝑀]⁻¹[𝐾]), respectively. It is clear that [𝑀]and[𝐾] are (2𝑛)×(2𝑚) matrices containing 𝑆_𝑃(𝑖,𝑗)(𝑃=1,2,…,6) and [𝐿] is 2𝑚×3 matrix containing 𝑆_𝑄(𝑖,𝑗)(𝑄=7,8,9). Also {𝐻} is an unknown vector containing 𝑖+𝑗 constants of integration. Next, upon imposing the remaining boundary conditions (i.e., 𝐴_𝑖𝑦=𝐴_𝑗𝑦𝑧=0 at either 𝑦=𝑏 or −𝑏) in (12), the unknown vector {𝐻} is determined in terms of vector{𝐵}.

B.

The coefficients 𝐴^𝑘𝑗, 𝐵^𝑘𝑗, and 𝐷^𝑘𝑗 appearing in (19) are given as 𝐴^𝑘𝑗₃₃=^𝑁_𝑖=1^𝑧𝑖+1_𝑧𝑖𝐶^(𝑖)₃₃𝑑Φ_𝑘𝑑𝑧𝑑Φ_𝑗𝐴𝑑𝑧𝑑𝑧,_{𝑘𝑗𝑝𝑞}=^𝑁_𝑖=1^𝑧𝑖+1_𝑧𝑖𝐶^(𝑖)_𝑝𝑞𝑑Φ_𝑘𝑑𝑧𝑑Φ_𝑗𝐵𝑑𝑧𝑑𝑧𝑝𝑞=44,55,_{𝑘𝑗𝑝𝑞}=^𝑁_𝑖=1^𝑧𝑖+1_𝑧𝑖𝐶^(𝑖)_𝑝𝑞Φ_𝑘𝑑Φ_𝑗𝐵𝑑𝑧𝑑𝑧𝑝𝑞=13,23,_{𝑘𝑗𝑝𝑞}=^𝑁_𝑖=1^𝑧𝑖+1_𝑧𝑖𝐶^(𝑖)_𝑝𝑞Φ_𝑘𝑑Φ_𝑗𝐷𝑑𝑧𝑑𝑧𝑝𝑞=44,55,_{𝑘𝑗𝑝𝑞}=^𝑁_𝑖=1^𝑧𝑖+1_𝑧𝑖𝐶^(𝑖)_𝑝𝑞Φ_𝑘Φ_𝑗𝐷𝑑𝑧𝑝𝑞=11,12,22,66,_{𝑘𝑗𝑝𝑞}=^𝑁_𝑖=1^𝑧𝑖+1_𝑧𝑖𝐶^(𝑖)_𝑝𝑞Φ_𝑘Φ_𝑗𝑑𝑧𝑝𝑞=44,55. (B.1) The linear global interpolation function is defined as Φ_𝑘(𝑧)=0𝑧≤𝑧_𝑘−1𝜓²_𝑘−1(𝑧)𝑧_𝑘−1≤𝑧≤𝑧_𝑘(),𝜓𝑘=1,2,…,𝑁+1_1𝑘(𝑧)𝑧_𝑘≤𝑧≤𝑧_𝑘+10𝑧≥𝑧_𝑘+1, (B.2) where 𝜓^𝑗_𝑘(𝑗=1,2) are the local Lagrangian linear interpolation functions within the 𝑘th layer which are defined as 𝜓¹_𝑘(𝑧)=1ℎ_𝑘𝑧_𝑘+1−𝑧,𝜓²_𝑘(𝑧)=1ℎ_𝑘𝑧−𝑧_𝑘, (B.3) with ℎ_𝑘 being the thickness of the 𝑘th layer.

C.

C.1. Analytical Solutions

The investigation is performed when analytically a rectangular composite laminate is subjected to the bending moment, the axial force, and the torque at its two opposite ends (𝑥=−𝑎 and 𝑥=𝑎) (see Figure 1). The linear Lagrangian interpolation functions through the thickness are used as in [26]. It is assumed that the laminates have invariably the traction-free boundary conditions at 𝑦=𝑏 and 𝑦=−𝑏. The boundary conditions at these edges are considered as 𝑀^𝑘_𝑥𝑦=𝑀^𝑘_𝑦=𝑅^𝑘_𝑦=0at𝑦=±𝑏. (C.1) Here the edges of the laminate are subjected to the bending moment 𝑀₀, the axial force 𝐹₀, and the torque 𝑇₀ at 𝑥=±𝑎: 𝑅^𝑘_𝑥=−2𝑇₀^ℎ/2_−ℎ/2𝑑Φ_𝑘𝑀𝑑𝑧𝑑𝑧,_𝑘𝑥𝑦=−2𝑇₀^ℎ/2_−ℎ/2Φ_𝑘(𝑧)1𝑧𝑀−1𝑑𝑧,_𝑘𝑥=−2𝑀₀^ℎ/2_−ℎ/2Φ_𝑘(𝑧)𝑧𝐹𝑑𝑧−2⁰𝐴^ℎ/2_−ℎ/2Φ_𝑘(𝑧)𝑑𝑧at𝑥=±𝑎. (C.2) To find analytical solutions for (C.2), it is assumed that 𝑈_𝑗𝑥,𝑦=𝑢_𝑗𝑥,𝑦+𝑈_𝑗(𝑥),𝑉_𝑗𝑥,𝑦=𝑣_𝑗𝑥,𝑦+𝑉_𝑗(𝑥),𝑊_𝑗𝑥,𝑦=𝑤_𝑗𝑥,𝑦+𝑊_𝑗(𝑥),𝑗=1,2,…,𝑁+1. (C.3) Upon substitution of these expressions into (20), two sets of equations are obtained. The first set contains 𝑈_𝑗(𝑥), 𝑉_𝑗(𝑥), and 𝑊_𝑗(𝑥) which are expressed by 𝐷^𝑘𝑗₁₁𝑈_𝑗^′′−𝐴^𝑘𝑗₅₅𝑈_𝑗+𝐵^𝑘𝑗₁₃−𝐵^𝑗𝑘₅₅𝑊^′_𝑗𝐷=0,_𝑘𝑗66𝑉_𝑗^′′−𝐴^𝑘𝑗₄₄𝑉_𝑗𝐵=0,_𝑘𝑗55−𝐵^𝑗𝑘₁₃𝑈^′_𝑗+𝐷^𝑘𝑗₅₅𝑊_𝑗^′′−𝐴^𝑘𝑗₃₃𝑊_𝑗=0. (C.4) The second set of equations contains 𝑢_𝑗(𝑥,𝑦), 𝑣_𝑗(𝑥,𝑦), and 𝑤_𝑗(𝑥,𝑦) which are given by 𝐷^𝑘𝑗₁₁𝑢_{𝑗,𝑥𝑥}+𝐷^𝑘𝑗₆₆𝑢_{𝑗,𝑦𝑦}−𝐴^𝑘𝑗₅₅𝑢_𝑗+𝐷^𝑘𝑗₁₂+𝐷^𝑘𝑗₆₆𝑣_{𝑗,𝑥𝑦}+𝐵^𝑘𝑗₁₃−𝐵^𝑗𝑘₅₅𝑤_𝑗,𝑥𝐷=0,_𝑘𝑗12+𝐷^𝑘𝑗₆₆𝑢_{𝑗,𝑥𝑦}+𝐷^𝑘𝑗₆₆𝑣_{𝑗,𝑥𝑥}+𝐷^𝑘𝑗₂₂𝑣_{𝑗,𝑦𝑦}−𝐴^𝑘𝑗₄₄𝑣_𝑗+𝐵^𝑘𝑗₂₃−𝐵^𝑗𝑘₄₄𝑤_𝑗,𝑦𝐵=0,_𝑘𝑗55−𝐵^𝑗𝑘₁₃𝑢_𝑗,𝑥+𝐵^𝑘𝑗₄₄−𝐵^𝑗𝑘₂₃𝑣_𝑗,𝑦+𝐷^𝑘𝑗₅₅𝑤_{𝑗,𝑥𝑥}+𝐷^𝑘𝑗₄₄𝑤_{𝑗,𝑦𝑦}−𝐴^𝑘𝑗₃₃𝑤_𝑗=0. (C.5) It should be emphasized that a repeated index in (C.4) and (C.5) is a summation index from 1 to 𝑁+1. When the boundary conditions (C.4) are considered, substituting (C.3) into (C.2) yields𝐵^𝑘𝑗₅₅𝑈_𝑗+𝐷^𝑘𝑗₅₅𝑊^′_𝑗=−2𝑇₀𝑛^𝑘_ℎ,𝐷^𝑘𝑗₆₆𝑉^′_𝑗=−2𝑇₀𝑛^𝑘_𝑧−𝑛^𝑘,𝐷^𝑘𝑗₁₁𝑈^′_𝑗+𝐵^𝑘𝑗₁₃𝑊_𝑗=−2𝑀₀𝑛^𝑘_𝑧𝐹−2⁰𝐴𝑛^𝑘at𝑥=±𝑎, (C.6a)(C.6b)(C.6c)where 𝑛^𝑘_𝑧=^ℎ/2_−ℎ/2Φ_𝑘(𝑧)𝑧𝑑𝑧,𝑛^𝑘=^ℎ/2_−ℎ/2Φ_𝑘(𝑧)𝑛𝑑𝑧,_𝑘ℎ=^ℎ/2_−ℎ/2𝑑Φ_𝑘𝑑𝑧𝑑𝑧,𝐵^𝑘𝑗₅₅𝑢_𝑗+𝐷^𝑘𝑗₅₅𝑤_𝑗,𝑥=𝐷^𝑘𝑗₆₆𝑢_𝑗,𝑦+𝐷^𝑘𝑗₆₆𝑣_𝑗,𝑥𝐷=0,_𝑘𝑗11𝑢_𝑗,𝑥+𝐷^𝑘𝑗₁₂𝑣_𝑗,𝑦+𝐵^𝑘𝑗₁₃𝑤_𝑗=0,at𝑥=±𝑎. (C.7)(C.8) Equations (C.4) subjected to the boundary conditions (C.6a), (C.6b), and (C.6c) can be solved analytically. It is to be noted that there exist repeated zero roots (or eigenvalues) in the characteristic equation of the set of equations in (C.4). To improve the solution scheme of these equations, some small artificial terms are added to these equations so that the characteristic roots become all distinct. So, (C.4) are rewritten as 𝐷^𝑘𝑗₁₁𝑈_𝑗^′′−𝐴^𝑘𝑗₅₅𝑈_𝑗+𝐵^𝑘𝑗₁₃−𝐵^𝑗𝑘₅₅𝑊^′_𝑗=𝛼^𝑘𝑗𝑈_𝑗,𝐷^𝑘𝑗₆₆𝑉_𝑗^′′−𝐴^𝑘𝑗₄₄𝑉_𝑗=𝛼^𝑘𝑗𝑉_𝑗,𝐵^𝑘𝑗₅₅−𝐵^𝑗𝑘₁₃𝑈^′_𝑗+𝐷^𝑘𝑗₅₅𝑊_𝑗^′′−𝐴^𝑘𝑗₃₃𝑊_𝑗=𝛼^𝑘𝑗𝑊_𝑗, (C.9) where, for convenience, 𝛼^𝑘𝑗 is assumed to be 𝛼^𝑘𝑗=𝛼^ℎ/2_−ℎ/2Φ_𝑘Φ_𝑗𝑑𝑧, (C.10) where 𝛼 is a prescribed number such that 𝛼^𝑘𝑗’s in (C.10) are comparatively small compared to the numerical values of stiffnesses 𝐴^𝑘𝑗_𝑝𝑞 (𝑝𝑞=33, 44, 55). The system of equations appearing in (C.9) is 3(𝑁+1) coupled second-order ordinary differential equations with constant coefficients which may be introduced in a matrix form as [𝑀]𝜂^′′+[𝐾]𝜂={0}, (C.11) where𝜂={𝑈}^𝑇,{𝑉}^𝑇,𝑊^𝑇𝑇,{𝑈}=𝑈₁,𝑈₂,…,𝑈_𝑁+1^𝑇,{𝑉}=𝑉₁,𝑉₂,…,𝑉_𝑁+1^𝑇,𝑊=𝑊₁,𝑊₂,…,𝑊_𝑁+1^𝑇,𝑊_𝑗=^𝑦𝑊_𝑗𝑑𝑥. (C.12a)(C.12b) The coefficient matrices [𝑀] and [𝐾] appearing in (C.11) are defined as [𝑀]=𝐷₁₁[0]𝐵₁₃−𝐵₅₅^𝑇[0]𝐷₆₆0][00𝐷[][]₅₅,[𝐾]=−𝐴₅₅−[𝛼00][0]−𝐴][][₄₄−[𝛼0]𝐵][₅₅−𝐵₁₃^𝑇[0]−𝐴₃₃−[𝛼],[𝐻]=𝑍₁⁻¹𝑍₂𝑍⋅𝑎,₂=−2𝑇₀𝑛_ℎ−2𝑇₀𝑛_𝑧−{𝑛}−2𝑀₀𝑛_𝑧𝐹−2⁰𝐴{𝑛},𝑍₁=𝐻₁𝜓𝜆²sinh(𝜆𝑎)+𝐻₂𝜓(,𝐻[sinh𝜆𝑎)]₁=[0][0]𝐷₅₅[000][000],𝐻][][][][₂=𝐵₅₅00][0𝐷[][]₆₆0]𝐷₁₁[[0]𝐵₁₃. (C.12c)(C.12d)The general solution of (C.11) can be written as 𝜂=𝜓[sinh(𝜆𝑥)]{𝐻}, (C.13) where [sinh(𝜆𝑥)], [𝜓], and (𝜆²₁,𝜆²₂,…,𝜆²_3(𝑁+1)) are a 3(𝑁+1)×3(𝑁+1) diagonal matrix, the modal matrix, and eigenvalues of (−[𝑀]⁻¹[𝐾]), respectively. In addition, {𝐻} is an unknown vector containing 3(𝑁+1) integration constants. It should be noted that by satisfying the boundary condition in (C.6a), (C.6b), and (C.6c) the integration constants in {𝐻} are found; then the problem is solved entirely.

It remains to solve (C.5) subject to the boundary conditions in (C.1) and (C.8). To this end, it is emphasized that the boundary conditions in (C.8) will identically be satisfied if the following expressions for the displacement components in (C.5) are assumed: 𝑢_𝑘𝑥,𝑦=^∞_𝑚=1𝑓^𝑘_1𝑚𝑦sin𝛼_1𝑚𝑥+𝑓^𝑘_2𝑚𝑦cos𝛼_2𝑚𝑥+𝑓^𝑘₁₀𝑦,𝑣_𝑘𝑥,𝑦=^∞_𝑚=1𝑓^𝑘_3𝑚𝑦sin𝛼_2𝑚𝑥+𝑓^𝑘_4𝑚𝑦cos𝛼_1𝑚𝑥+𝑓^𝑘₂₀𝑦,𝑤_𝑘𝑥,𝑦=^∞_𝑚=1𝑓^𝑘_5𝑚𝑦sin𝛼_2𝑚𝑥+𝑓^𝑘_6𝑚𝑦cos𝛼_1𝑚𝑥+𝑓^𝑘₃₀𝑦,𝑓^𝑘′₃₀𝑦+𝑓^𝑘₂₀𝑦+(−1)^𝑚×𝐷^𝑘𝑗₁₁𝑓^𝑘_1𝑚𝛼_1𝑚+𝐷^𝑘𝑗₁₁𝑓^𝑘_2𝑚𝛼_2𝑚−𝐷^𝑘𝑗₁₂𝑓^𝑘′_3𝑚+𝐷^𝑘𝑗₁₂𝑓^𝑘′_4𝑚−𝐵^𝑘𝑗₁₃𝑓^𝑘_5𝑚+𝐵^𝑘𝑗₁₃𝑓^𝑘_6𝑚𝑓=0,_𝑘10𝑦=0, (C.14)(C.15) where 𝛼_1𝑚=𝑚𝜋/𝑎, 𝛼_2𝑚=𝑚𝜋/2𝑎 with 𝑚 are the Fourier integer. Upon substitution of (C.14) into (C.5), the following ordinary differential equations are obtained: 𝐷^𝑘𝑗₆₆𝑓^𝑗′′_1𝑚−𝐷^𝑘𝑗₁₁𝛼²_1𝑚+𝐴^𝑘𝑗₅₅𝑓^𝑗_1𝑚−𝛼_1𝑚𝐷^𝑘𝑗₁₂+𝐷^𝑘𝑗₆₆𝑓^𝑗′_4𝑚−𝛼_1𝑚𝐵^𝑘𝑗₁₃−𝐵^𝑗𝑘₅₅𝑓^𝑗_6𝑚𝛼=0,_1𝑚𝐷^𝑘𝑗₁₂+𝐷^𝑘𝑗₆₆𝑓^𝑗′_1𝑚+𝐷^𝑘𝑗₂₂𝑓^𝑗′′_4𝑚−𝐷^𝑘𝑗₆₆𝛼²_1𝑚+𝐴^𝑘𝑗₄₄𝑓^𝑗_4𝑚−𝐵^𝑘𝑗₂₃−𝐵^𝑗𝑘₄₄𝑓^𝑗′_6𝑚𝛼=0,_1𝑚𝐵^𝑘𝑗₅₅−𝐵^𝑗𝑘₁₃𝑓^𝑗_1𝑚+𝐵^𝑘𝑗₄₄−𝐵^𝑗𝑘₂₃𝑓^𝑗′_4𝑚+𝐷^𝑘𝑗₄₄𝑓^𝑗′′_6𝑚−𝐴^𝑘𝑗₃₃+𝐷^𝑘𝑗₅₅𝛼²_1𝑚𝑓^𝑗_6𝑚=0,𝐷^𝑘𝑗₆₆𝑓^𝑗′′_2𝑚−𝐷^𝑘𝑗₁₁𝛼²_2𝑚+𝐴^𝑘𝑗₅₅𝑓^𝑗_2𝑚−𝛼_2𝑚𝐷^𝑘𝑗₁₂+𝐷^𝑘𝑗₆₆𝑓^𝑗′_3𝑚−𝛼_2𝑚𝐵^𝑘𝑗₁₃−𝐵^𝑗𝑘₅₅𝑓^𝑗_5𝑚=0,−𝛼_2𝑚𝐷^𝑘𝑗₁₂+𝐷^𝑘𝑗₆₆𝑓^𝑗′_2𝑚+𝐷^𝑘𝑗₂₂𝑓^𝑗′′_3𝑚−𝐷^𝑘𝑗₆₆𝛼²_2𝑚+𝐴^𝑘𝑗₄₄𝑓^𝑗_3𝑚−𝐵^𝑘𝑗₂₃−𝐵^𝑗𝑘₄₄𝑓^𝑗′_5𝑚=0,−𝛼_2𝑚𝐵^𝑘𝑗₅₅−𝐵^𝑗𝑘₁₃𝑓^𝑗_2𝑚+𝐵^𝑘𝑗₄₄−𝐵^𝑗𝑘₂₃𝑓^𝑗′_3𝑚+𝐷^𝑘𝑗₄₄𝑓^𝑗′′_5𝑚−𝐴^𝑘𝑗₃₃+𝐷^𝑘𝑗₅₅𝛼²_2𝑚𝑓^𝑗_5𝑚=0, (C.16a)(C.16b)𝐷^𝑘𝑗₂₂𝑓^𝑗′′₂₀+𝐵^𝑘𝑗₂₃−𝐵^𝑗𝑘₄₄𝑓^𝑗′₃₀−𝐴^𝑘𝑗₄₄𝑓^𝑗₂₀𝐷=0,_𝑘𝑗44𝑓^𝑗′′₃₀+𝐵^𝑗𝑘₄₄−𝐵^𝑘𝑗₂₃𝑓^𝑗′₂₀−𝐴^𝑘𝑗₃₃𝑓^𝑗₃₀=0. (C.17) Considering (C.14), the boundary conditions in (C.1) at 𝑦=𝑏 and 𝑦=−𝑏 are given as 𝑀^𝑘_𝑥𝑦=^∞_𝑚=1𝐷^𝑘𝑗₆₆𝑓^𝑗′_1𝑚−𝐷^𝑘𝑗₆₆𝑓^𝑗_4𝑚𝛼_1𝑚sin𝛼_1𝑚𝑥+𝐷^𝑘𝑗₆₆𝑓^𝑗′_2𝑚+𝐷^𝑘𝑗₆₆𝑓^𝑗_3𝑚𝛼_2𝑚cos𝛼_2𝑚𝑥+𝐷^𝑘𝑗₆₆𝑉_𝑗^′𝑀=0,_𝑘𝑦=^∞_𝑚=1−𝐷^𝑘𝑗₁₂𝑓^𝑗_2𝑚𝛼_2𝑚+𝐷^𝑘𝑗₂₂𝑓^𝑗′_3𝑚+𝐵^𝑘𝑗₂₃𝑓^𝑗_5𝑚sin𝛼_2𝑚𝑥+𝐷^𝑘𝑗₁₂𝑓^𝑗_1𝑚𝛼_1𝑚+𝐷^𝑘𝑗₂₂𝑓^𝑗′_4𝑚+𝐵^𝑘𝑗₂₃𝑓^𝑗_6𝑚cos𝛼_1𝑚𝑥+𝐷^𝑘𝑗₂₂𝑓^𝑗′₂₀+𝐵^𝑘𝑗₂₃𝑓^𝑗₃₀+𝐷^𝑘𝑗₁₂𝑈_𝑗^′+𝐵^𝑘𝑗₂₃𝑊_𝑗𝑅=0,_𝑘𝑦=^∞_𝑚=1𝐵^𝑘𝑗₄₄𝑓^𝑗′_3𝑚+𝐷^𝑘𝑗₄₄𝑓^𝑗′_5𝑚sin𝛼_2𝑚𝑥+𝐵^𝑘𝑗₄₄𝑓^𝑗′_4𝑚+𝐷^𝑘𝑗₄₄𝑓^𝑗′_6𝑚cos𝛼_1𝑚𝑥+𝐵^𝑘𝑗₄₄𝑉_𝑗+𝐷^𝑘𝑗₄₄𝑓^𝑗′₃₀+𝐵^𝑘𝑗₄₄𝑓^𝑗₂₀=0. (C.18) By introducing Fourier sine and cosine expansions for the underlined terms in (C.18), the following boundary conditions at 𝑦=𝑏 and 𝑦=−𝑏 can readily be obtained as 𝐷^𝑘𝑗₂₂𝑓^𝑗′₂₀+𝐵^𝑘𝑗₂₃𝑓^𝑗₃₀𝐽=−^𝑘02,𝐷^𝑘𝑗₄₄𝑓^𝑗′₃₀+𝐵^𝑘𝑗₄₄𝑓^𝑗₂₀𝐾=−^𝑘02, (C.19)𝐷^𝑘𝑗₆₆𝑓^𝑗′_1𝑚−𝐷^𝑘𝑗₆₆𝑓^𝑗_4𝑚𝛼_1𝑚=−𝐼^𝑘_1𝑚,𝐷^𝑘𝑗₁₂𝑓^𝑗_1𝑚𝛼_1𝑚+𝐷^𝑘𝑗₂₂𝑓^𝑗′_4𝑚+𝐵^𝑘𝑗₂₃𝑓^𝑗_6𝑚=−𝐽^𝑘_1𝑚,𝐵^𝑘𝑗₄₄𝑓^𝑗′_4𝑚+𝐷^𝑘𝑗₄₄𝑓^𝑗′_6𝑚=−𝐾^𝑘_1𝑚,𝐷^𝑘𝑗₆₆𝑓^𝑗′_2𝑚+𝐷^𝑘𝑗₆₆𝑓^𝑗_3𝑚𝛼_2𝑚=−𝐼^𝑘_2𝑚,−𝐷^𝑘𝑗₁₂𝑓^𝑗_2𝑚𝛼_2𝑚+𝐷^𝑘𝑗₂₂𝑓^𝑗′_3𝑚+𝐵^𝑘𝑗₂₃𝑓^𝑗_5𝑚=−𝐽^𝑘_2𝑚,𝐵^𝑘𝑗₄₄𝑓^𝑗′_3𝑚+𝐷^𝑘𝑗₄₄𝑓^𝑗′_5𝑚=−𝐾^𝑘_2𝑚, (C.20a)(C.20b)where𝐷^𝑘𝑗₆₆𝑉^′_𝑗=𝐼^𝑘(𝑥)=^∞_𝑚=1𝐼^𝑘_1𝑚sin𝛼_1𝑚𝑥+𝐼^𝑘_2𝑚cos𝛼_2𝑚𝑥+𝐼^𝑘₀2𝐷^𝑘𝑗₁₂𝑈^′_𝑗+𝐵^𝑘𝑗₂₃𝑊_𝑗=𝐽^𝑘(𝑥)=^∞_𝑚=1𝐽^𝑘_2𝑚sin𝛼_2𝑚𝑥+𝐽^𝑘_1𝑚cos𝛼_1𝑚𝑥+𝐽^𝑘₀2𝐵^𝑘𝑗₄₄𝑉_𝑗=𝐾^𝑘(𝑥)=^∞_𝑚=1𝐾^𝑘_2𝑚sin𝛼_2𝑚𝑥+𝐾^𝑘_1𝑚cos𝛼_1𝑚𝑥+𝐾^𝑘₀2, (C.21a) where the coefficients 𝐼^𝑘, 𝐾^𝑘, and 𝐽^𝑘 appearing in (C.21a) are given as 𝐼^𝑘_1𝑚=1𝑎^𝑎_−𝑎𝐼^𝑘(𝑥)sin𝛼_1𝑚𝑥𝐾𝑑𝑥,_𝑘1𝑚=1𝑎^𝑎_−𝑎𝐾^𝑘(𝑥)cos𝛼_1𝑚𝑥𝐼𝑑𝑥,_𝑘2𝑚=1𝑎^𝑎_−𝑎𝐼^𝑘(𝑥)cos𝛼_2𝑚𝑥𝐾𝑑𝑥,_𝑘2𝑚=1𝑎^𝑎_−𝑎𝐾^𝑘(𝑥)sin𝛼_2𝑚𝑥𝐽𝑑𝑥,_𝑘1𝑚=1𝑎^𝑎_−𝑎𝐽^𝑘(𝑥)cos𝛼_1𝑚𝑥𝐽𝑑𝑥,_𝑘0=1𝑎^𝑎_−𝑎𝐽^𝑘(𝑥)𝐽𝑑𝑥,_𝑘2𝑚=1𝑎^𝑎_−𝑎𝐽^𝑘(𝑥)sin𝛼_2𝑚𝑥𝑑𝑥,𝐾^𝑘₀=1𝑎^𝑎_−𝑎𝐾^𝑘(𝑥)𝑑𝑥. (C.21b)(C.21c)In order to obtain the solution of (C.16a) and (C.16b), it is assumed that the boundary conditions at 𝑦=−𝑏 and 𝑦=𝑏 appearing in (C.20a) and (C.20b) are identical. To be able to impose the boundary conditions only at one edge, say, at 𝑦=𝑏 and, a significant advantage of this is to save some computational time, the system of equations appearing in (C.16a) are 3(𝑁+1) coupled second-order ordinary differential equations with constant coefficients which are introduced in a matrix form as [𝑆]𝜇^′′+[𝐺]𝜇={0}, (C.22) where𝜇=𝑓^𝑗𝑇_1𝑚,𝑓^𝑗𝑇_4𝑚,𝑓^𝑗𝑇_6𝑚^𝑇,𝑓^𝑗_1𝑚=𝑓¹_1𝑚,𝑓²_1𝑚,…,𝑓^𝑁+1_1𝑚^𝑇,𝑓^𝑗_4𝑚=𝑓¹_4𝑚,𝑓²_4𝑚,…,𝑓^𝑁+1_4𝑚^𝑇,𝑓^𝑗_6𝑚=𝑓¹_6𝑚,𝑓²_6𝑚,…,𝑓^𝑁+1_6𝑚^𝑇,𝑓^𝑗_4𝑚=^𝑦𝑓^𝑗_4𝑚𝑑𝑦. (C.23a)(C.23b)

The coefficient matrices [𝑆] and [𝐺] appearing in (C.22) for solving (C.16a) are given as follows:[𝑆]=𝐷₆₆−𝛼_1𝑚𝐷₁₂+𝐷₆₆[0][0]𝐷₂₂0][0𝐵[]₄₄−𝐵₂₃^𝑇𝐷₄₄,[𝐺]=−𝐴₅₅+𝛼²_1𝑚𝐷₁₁[0]−𝛼_1𝑚𝐵₁₃−𝐵₅₅^𝑇𝛼_1𝑚𝐷₁₂+𝐷₆₆−𝐴₄₄+𝛼²_1𝑚𝐷₆₆𝐵₂₃−𝐵₄₄^𝑇𝛼_1𝑚𝐵₅₅−𝐵₁₃^𝑇[0]−𝐴₃₃+𝛼²_1𝑚𝐷₅₅. (C.23c)

The coefficient matrices [𝑆] and [𝐺] appearing in (C.22) for solving (C.16b) are given as[𝑆]=𝐷₆₆−𝛼_2𝑚𝐷₁₂+𝐷₆₆[0][0]𝐷₂₂0][0𝐵[]₄₄−𝐵₂₃^𝑇𝐷₄₄,[𝐺]=−𝐴₅₅+𝛼²_2𝑚𝐷₁₁[0]−𝛼_2𝑚𝐵₁₃−𝐵₅₅^𝑇𝛼_2𝑚𝐷₁₂+𝐷₆₆−𝐴₄₄+𝛼²_2𝑚𝐷₆₆𝐵₂₃−𝐵₄₄^𝑇𝛼_2𝑚𝐵₅₅−𝐵₁₃^𝑇[0]−𝐴₃₃+𝛼²_2𝑚𝐷₅₅. (C.23d)The coefficient matrices [𝑆] and [𝐺] appearing in (C.22) for solving (C.17) are given as [𝑆]=𝐷₂₂𝐵₂₃−𝐵₄₄^𝑇𝐷^𝑘𝑗₄₄𝐵₄₄−𝐵₂₃,[𝐺]=−𝐴₄₄0][0]−𝐴[₃₃. (C.23e)The general solution of (C.22) can be written as 𝜇=𝜉sinh𝜆𝑦{𝐵}, (C.24) where [sinh(𝜆𝑦)], [𝜉], and (𝜆²₁,𝜆²₂,…,𝜆²_3(𝑁+1)) are a 3(𝑁+1)×3(𝑁+1) diagonal matrix, the modal matrix, and eigenvalues of (−[𝑆]⁻¹[𝐺]), respectively. Moreover, {𝐵} is an unknown vector containing 3(𝑁+1) integration constants. By satisfying the boundary condition in (C.20a) only at one edge, say, at 𝑦=𝑏, the integration constants in {𝐵} are obtained; then the problem is solved completely.

The solution procedure for (C.16b) with the boundary conditions in (C.20b) and (C.17) and with the boundary conditions in (C.19) is similar to the one discussed in (C.16b), and therefore, for the sake of brevity, it will not be taken up here. The boundary conditions in (C.15) will identically be satisfied.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.