Analytic solution of angle-ply laminated plates under extension,bending,and torsion

Abstract

This paper develops a Hamiltonian state space approach for analytic determination of deformation and stress fields in multilayered monoclinic angle-ply laminates under the combined action of extension, bending, and torsion. The present solution satisfies the equations of anisotropic elasticity, the end conditions, the traction-free boundary conditions on the four edge surfaces of the rectangular section, and the interfacial continuity conditions in multilayered laminates. The proposed method only requires the solutions of matrix and eigen equations, regardless of the number or lamination of the layers. The finite element analyses are used to validate the accuracy of the analysis. The analytical solution and the numerical solutions are in excellent agreement.

Keywords

Composite laminates finite element analysis free edges generalized plane strain monoclinic composite plate state space formalism symplectic orthogonality

Introduction

This paper determines the deformation and stress fields of multilayered monoclinic angle-ply composite laminates theoretically. Although this class of problems has been studied for many years, there are still few exact solutions in the literature, because it is difficult to satisfy both the interfacial continuity conditions and the edge boundary conditions. For the traction-free boundary conditions along the edges, Saint-Venant's principle can be used such that the stress resultants across the thickness of the free edge are equal to zero, and the solutions then obtained are considered to be valid away from the edge surfaces.^1,2 The composite problems with simply supported boundary conditions have been solved in the literature. Pagano³ presented the well-known analytic solutions on bending for cross-ply laminates. Ren⁴ presented exact solutions in terms of Fourier series for laminated cylindrical shells. Savoia and Reddy⁵ solved the three-dimensional (3D) problem by minimizing the total potential energy functional for laminated plates. Heyliger⁶ presented exact solutions for the piezoelectric laminate subjected to transverse loadings. Pan⁷ derived exact solutions for 3D magneto-electro-elastic plates under static loadings. Demasi⁸ developed a Navier-type method to obtain an exact 3D solution for simply supported rectangular plates. Yang et al.⁹ used pseudo-Stroh formalism to derive an exact solution for simply supported plates under surface loadings.

There are a number of analytical solutions for the plates with free edge boundary conditions, and detailed reviews of these can be found by Mittelstedt and Becker.¹⁰ Cho and Kim¹¹ used an iterative method to analyze the free-edge stresses of composite laminates subjected to extension, bending, twisting, and thermal loads using the extended Kantorovich method. Nosier and Maleki¹² used an improved first-order shear deformation theory to study free-edge stresses in composite laminates subjected to extension loads. Sarvestani and Sarvestani¹³ solved the interlaminar stress in composite laminates subjected to extension, torsion, and bending moment based on the reduced form of elasticity displacement field. Saeedi et al.¹⁴ analyzed multilayered plates subjected to invariant loading along the longitudinal direction using the layerwise stress model (LS1), and excellent agreement with regard to the free-edge interlaminar stresses were found between the LS1 and finite element models. Tarn and Chang¹⁵ formulated the basic equations of anisotropic elasticity and piezoelasticity into the state space framework using a Hamiltonian variation formulation. Andakhshideh and Tahani¹⁶ developed a 3D elasticity-based solution using the iterative analytical method for determining the free edge stress field in laminates under uniaxial extension, bending, twisting, and thermal loading. Dhanesh et al.¹⁷ extended this iterative analytical method to solve free-edge laminate problems, while this approach satisfies all the boundary and interfacial continuity conditions. Wu et al.¹⁸ presented a hybrid analysis of the interfacial stresses near the free edges of piezoelectric laminated plates under uniform extension using a state space equation for cross-ply asymmetric piezoelectric laminates. Alipour¹⁹ presented an analytical approach for the static analysis of sandwich circular plates based on a layerwise theory linked with the 3D theory of elasticity. Liang et al.²⁰ developed a Hamiltonian state space approach for the analytic determination of deformation and stress fields in multilayered cross-ply laminates under extension and bending, while the laminates should be arranged symmetrically in 0 ° and 90 °. Schnabel et al.²¹ presented an analysis method for the determination of displacement, strain, and stress fields in cylindrically curved cross-ply laminates under bending load, and their method was compared to the results of finite element simulations in good agreement. Hajikazemi and Van Paepegem²² developed a variational model to determine stress and displacement fields at free edges of symmetric composite laminate strips under thermomechanical loads, and the method is also applicable to thin-ply laminates. Vukasović et al.²³ presented an analytic solution based on the classical Vlasov's beam theory for torsion of thin-walled laminated composite beams of symmetrical open cross sections, and this analytical model can be efficiently used as a tool in the early design stages of thin-walled laminated composite beams. Yan et al.²⁴ developed an exact solution based on refined beam theory which uses Lagrange-polynomials expansions for the multiscale analysis of composite laminates and sandwich beams, and the stress/displacement fields at different scales can be successfully detected by increasing the order of Lagrange polynomials opportunely. Wackerfuß and Kroker²⁵ presented an efficient semi-analytical simulation framework that allows a highly systematic analysis of the mechanical behavior of laminated prismatic thin-walled beams, and this approach is suitable for extensive case studies. Cas et al.²⁶ presented a new strain-based finite element formulation that consider non-linear material and contact model of steel-concrete beam. The computational efficiency and accuracy of the formulation was proved by comparing the numerical results with the experimental results of full-scale laboratory test.

The literature review presented above indicates there are very few analytical solutions for finite monoclinic or angle-ply laminates plates, especially for the unsymmetrical arrangement of the laminate layers. Based on the Hamiltonian state space formalism,²⁷ the present study thus derives the analytic solutions for finite monoclinic laminates under the combination of extension, bending, and torsion, where the traction-free boundary conditions on the four edge surfaces of the rectangular section are well satisfied, as are the interfacial continuity conditions in multilayered laminates. The major advantage of the current method is that the symmetric arrangement of laminates is not necessary, and this advantage is rarely found in the literature. Moreover, this is the first investigation using the Hamiltonian state space approach for the above problem in the literature.

Illustration of the analyzed monoclinic composite plate

For monoclinic materials, the stress and strain relationship is as follows

\begin{matrix} [\begin{matrix} σ_{11} \\ σ_{22} \\ σ_{33} \\ σ_{23} \\ σ_{31} \\ σ_{12} \end{matrix}] = [\begin{matrix} C_{11} & C_{12} & C_{13} & 0 & C_{15} & 0 \\ C_{12} & C_{22} & C_{23} & 0 & C_{25} & 0 \\ C_{13} & C_{23} & C_{33} & 0 & C_{35} & 0 \\ 0 & 0 & 0 & C_{44} & 0 & C_{46} \\ C_{15} & C_{25} & C_{35} & 0 & C_{55} & 0 \\ 0 & 0 & 0 & C_{46} & 0 & C_{66} \end{matrix}] [\begin{matrix} ɛ_{11} \\ ɛ_{22} \\ ɛ_{33} \\ 2 ɛ_{23} \\ 2 ɛ_{31} \\ 2 ɛ_{12} \end{matrix}] \end{matrix}

(1)

where

σ_{ij}

is the stress,

ɛ_{ij}

is the strain, and

C_{ij}

is the elasticity matrix, which may be calculated from a cross-ply material with a rotation on its plane.

The analyzed problem is a composite plate composed of m elastic monoclinic layers of the rectangular section subjected to a bending moment, a tension force, and a torque at the end sections, as shown in Figure 1, where the only restriction of the problem is that the curvature of bending in the x₂ axis should be zero. Using Cartesian coordinates $(x_{1}, x_{2}, x_{3})$ with the origin at the center of the middle plane, the $x_{1}$ , $x_{2}$ , and $x_{3}$ axes are pointing in the width, thickness, and longitudinal directions, respectively, while the end sections are defined by $x_{2} = \pm h / 2$ , $x_{1} = \pm a / 2$ , and $x_{3} = 0$ to l.

Figure 1.

The free-edge problem for a rectangular angle-ply laminate plate subjected to an axial force P, a bending moment $M_{1}$ and a torque $M_{t}$ at the end sections.

The general expression of the displacement field is^27,28

\begin{matrix} [\begin{matrix} u_{1} \\ u_{2} \\ u_{3} \end{matrix}] = [\begin{matrix} u (x_{1}, x_{2}) \\ v (x_{1}, x_{2}) \\ w (x_{1}, x_{2}) \end{matrix}] + ɛ [\begin{matrix} 0 \\ 0 \\ x_{3} \end{matrix}] + b_{1} [\begin{matrix} 0 \\ - x_{3}^{2} / 2 \\ x_{2} x_{3} \end{matrix}] \\ + ϑ [\begin{matrix} - x_{2} x_{3} \\ x_{1} x_{3} \\ 0 \end{matrix}] \end{matrix}

(2)

where the rigid body displacement is excluded, u, v, and w are unknown functions, and ɛ, b₁, and ϑ are the uniform tension strain, the curvature of bending, and the curvature of twisting. On the basis of the state space formalism,²⁷ equations (1) and (2) can be formulated as the following two equations.

\begin{matrix} \frac{\partial}{\partial x_{2}} [\begin{matrix} u \\ v \\ w \\ σ_{12} \\ σ_{22} \\ σ_{23} \end{matrix}] = [\begin{matrix} 0 & - \partial 1 & 0 & s_{66} & 0 & s_{46} \\ - {\overset{\land}{c}}_{12} \partial 1 & 0 & - {\overset{\land}{c}}_{25} \partial 1 & 0 & c_{22}^{- 1} & 0 \\ 0 & 0 & 0 & s_{46} & 0 & s_{44} \\ - {\tilde{c}}_{11} \partial 11 & 0 & - {\tilde{c}}_{15} \partial 11 & 0 & - {\overset{\land}{c}}_{12} \partial 1 & 0 \\ 0 & 0 & 0 & - \partial 1 & 0 & 0 \\ - {\tilde{c}}_{15} \partial 11 & 0 & - {\tilde{c}}_{55} \partial 11 & 0 & - {\overset{\land}{c}}_{25} \partial 1 & 0 \end{matrix}] [\begin{matrix} u \\ v \\ w \\ σ_{12} \\ σ_{22} \\ σ_{23} \end{matrix}] \\ - (ɛ + b_{1} x_{2}) [\begin{matrix} 0 \\ {\overset{\land}{c}}_{23} \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}] - ϑ x_{1} [\begin{matrix} 0 \\ 0 \\ 1 \\ 0 \\ 0 \\ 0 \end{matrix}] + ϑ x_{2} [\begin{matrix} 0 \\ {\overset{\land}{c}}_{25} \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}] \end{matrix}

(3)

\begin{matrix} [\begin{matrix} σ_{11} \\ σ_{13} \\ σ_{33} \end{matrix}] = [\begin{matrix} {\tilde{c}}_{11} \partial 1 & {\tilde{c}}_{15} \partial 1 & {\overset{\land}{c}}_{12} \\ {\tilde{c}}_{15} \partial 1 & {\tilde{c}}_{55} \partial 1 & {\overset{\land}{c}}_{25} \\ {\tilde{c}}_{13} \partial 1 & {\tilde{c}}_{35} \partial 1 & {\overset{\land}{c}}_{23} \end{matrix}] [\begin{matrix} u \\ w \\ σ_{22} \end{matrix}] \\ + (ɛ + b_{1} x_{2}) [\begin{matrix} {\tilde{c}}_{13} \\ {\tilde{c}}_{35} \\ {\tilde{c}}_{33} \end{matrix}] - ϑ x_{2} [\begin{matrix} {\tilde{c}}_{15} \\ {\tilde{c}}_{55} \\ {\tilde{c}}_{35} \end{matrix}] \end{matrix}

(4)

where

\begin{matrix} {\overset{\land}{c}}_{ij} = c_{ij} c_{22}^{- 1}, {\tilde{c}}_{ij} = c_{ij} - c_{i 2} c_{22}^{- 1} c_{j 2}, \\ [\begin{matrix} s_{66} & s_{46} \\ s_{46} & s_{44} \end{matrix}] = \frac{1}{c_{44} c_{66} - c_{46}^{2}} [\begin{matrix} c_{44} & - c_{46} \\ - c_{46} & c_{66} \end{matrix}] \end{matrix}

(5)

The traction-free lateral boundary conditions are

\begin{matrix} [\begin{matrix} σ_{11} & σ_{12} & σ_{13} \end{matrix}]^{T} = 0 on x_{1} = \pm a / 2, and \\ [\begin{matrix} σ_{12} & σ_{22} & σ_{23} \end{matrix}]^{T} = 0 on x_{2} = \pm h / 2 \end{matrix}

(6)

The conditions of interfacial continuity require

[\begin{matrix} u & v & w & σ_{2} \end{matrix}]_{k + 1} = [\begin{matrix} u & v & w & σ_{2} \end{matrix}]_{k} on x_{2} = y_{(k)}

(7)

for

k = 1, 2, \dots, m - 1

, where m is the total number of layers,

y_{(k)}

denotes the thickness coordinate at the interface between the k^th layer and the (k + 1)^th layer. The end conditions require that the stress resultants over the end surfaces reduce to an axial force P, bending moments M₁, and a torque M_t at

x_{3} = 0

and l such that

\sum_{k = 1}^{m} \int_{y_{(k - 1)}}^{y_{(k)}} \int_{- a / 2}^{a / 2} σ_{33} d x_{1} d x_{2} = P

(8)

\sum_{k = 1}^{m} \int_{y_{(k - 1)}}^{y_{(k)}} \int_{- a / 2}^{a / 2} σ_{33} x_{2} d x_{1} d x_{2} = M_{1}

(9)

\sum_{k = 1}^{m} \int_{y_{(k - 1)}}^{y_{(k)}} \int_{- a / 2}^{a / 2} (x_{1} σ_{23} - x_{2} σ_{13}) k d x_{1} d x_{2} = M_{t}

(10)

To illustrate this problem, we use finite element analysis to determine the displacement field of the angle-ply laminates with stacking sequences of [–45°/+45°/–45°/+45°] under a uniform tension force.

The elastic constants of the rotation angle of 0° are $c_{11} = c_{22} = 15, 300$ $N / m m^{2},$ $c_{33} = 140, 000$ $N / m m^{2},$ $c_{44} = c_{55} = 5, 900$ $N / m m^{2},$ $c_{23} = c_{13} = 3, 900$ $N / m m^{2},$ and $c_{12} = 3300$ $N / m m^{2}$ . The dimensions of the laminate are taken to be $a / h = 10$ , where $a = 10$ mm and $h = 1$ mm are the width and thickness of the laminate, respectively. Each layer of the laminate has an equal thickness of h/4. Finite element analysis with 20-node high-order isoparametric elements was used, while the high-order element can avoid the shear locking due to the first-order elements using linear interpolation function for displacement. It is noted that shear locking is a serious problem in the finite element analysis.²⁹ Figure 2 shows the 3D finite element mesh in the x₁–x₂ plane, in which the fine mesh near the two edges of the x₁ axis is generated in order to obtain an accurate stress distribution near the plate boundary. In the x₃ direction, there are 140 element layers with the space of 4h, while the nodes are fixed along the right boundary of the x₃ axis, and a master node at the center of the x₁–x₂ plane is used to control the nodes along the left boundary of the x₃ axis. The applied force or moment is then applied to this master node. Figure 3 shows the deformed plate of the finite element analysis, which represents that if a uniform tension force is applied to the plate then it will produce not only the tension strain but also the rotation in the x₃ direction. It is noted that a very fine finite element mesh is required to solve this problem with acceptable accuracy, and a compressed computer time is necessary. Generally, a theoretical solution is the best way to understand the mechanical behavior of composite plates in detail.

Figure 2.

Finite element mesh in the x₁–x₂ plane for the analyzed composite plate with the thickness of h and width of a = 10h (There are 2928 20-node solid elements in one layer of the Z direction, and 140 layers are set (total number of elements = 409,920, and total number of degrees of freedom = 5,124,267).)

Figure 3.

Deformed plate with the displacement and rotation in the x₃ direction due to a uniform tension force on the angle-ply laminates with stacking sequences of [–45°/+45°/–45°/+45°].

The state space formalism for rectangular plates with monoclinic materials

Eigensolution

We use an eigensolution, as shown in equation (11), to satisfy the homogeneous equations, as shown in equations (3) and (4), in each layer.

\begin{matrix} [\begin{matrix} u \\ v \\ w \\ σ_{12} \\ σ_{22} \\ σ_{23} \end{matrix}] = \sum_{j = 1}^{6} C_{j} [\begin{matrix} U \\ V \\ W \\ S_{12} \\ S_{22} \\ S_{23} \end{matrix}] e^{χ_{j} μ y} e^{μ x}, \\ [\begin{matrix} σ_{11} \\ σ_{13} \\ σ_{33} \end{matrix}] = \sum_{j = 1}^{6} C_{j} [\begin{matrix} S_{11} \\ S_{13} \\ S_{33} \end{matrix}] e^{χ_{j} μ y} e^{μ x} \end{matrix}

(11)

One thus obtains the displacement and stress field of each layer as follows

\begin{matrix} [\begin{matrix} U \\ V \\ W \\ S_{12} \\ S_{22} \\ S_{23} \end{matrix}] = [\begin{matrix} α_{12} α_{23} - α_{13} α_{22} \\ μ^{- 1} (- χ μ U - s_{66} χ S_{22} + s_{46} S_{23}) \\ α_{13} α_{21} - α_{11} α_{23} \\ - χ S_{22} \\ α_{11} α_{22} - α_{12} α_{21} \\ s_{44}^{- 1} χ (μ W + s_{46} S_{22}) \end{matrix}] \\ [\begin{matrix} S_{11} \\ S_{13} \\ S_{33} \end{matrix}] = [\begin{matrix} {\tilde{c}}_{11} μ & {\tilde{c}}_{15} μ & {\overset{\land}{c}}_{12} \\ {\tilde{c}}_{15} μ & {\tilde{c}}_{55} μ & {\overset{\land}{c}}_{25} \\ {\tilde{c}}_{13} μ & {\tilde{c}}_{35} μ & {\overset{\land}{c}}_{23} \end{matrix}] [\begin{matrix} U \\ W \\ S_{22} \end{matrix}] \end{matrix}

(12)

where

\begin{matrix} [\begin{matrix} α_{11} & α_{12} & α_{13} \\ α_{21} & α_{22} & α_{23} \end{matrix}] = \\ \times [\begin{matrix} ({\overset{\land}{c}}_{12} - χ^{2}) μ & ({\overset{\land}{c}}_{25} + s_{44}^{- 1} s_{46} χ^{2}) μ & (s_{44}^{- 1} s_{46}^{2} - s_{66}) χ^{2} - c_{22}^{- 1} \\ {\tilde{c}}_{11} μ & {\tilde{c}}_{15} μ & {\overset{\land}{c}}_{12} - χ^{2} \end{matrix}] \end{matrix}

(13)

χ^{6} + η_{1} χ^{4} + η_{2} χ^{2} + η_{3} = 0

(14)

in which

\begin{matrix} η_{1} = s_{66} {\tilde{c}}_{11} + s_{44} {\tilde{c}}_{55} - 2 ({\overset{\land}{c}}_{12} - s_{46} {\tilde{c}}_{15}), \\ η_{3} = s_{44} c_{22}^{- 1} (c_{11} c_{55} - c_{15}^{2}) \\ η_{2} = {\overset{\land}{c}}_{11} + (s_{44} s_{66} - s_{46}^{2}) ({\tilde{c}}_{11} {\tilde{c}}_{55} - {\tilde{c}}_{15}^{2}) \\ + 2 s_{44} ({\overset{\land}{c}}_{25} {\tilde{c}}_{15} - {\overset{\land}{c}}_{12} {\tilde{c}}_{55}) + 2 s_{46} ({\overset{\land}{c}}_{25} {\tilde{c}}_{11} - {\overset{\land}{c}}_{12} {\tilde{c}}_{15}) \end{matrix}

Using the transfer matrix method with interfacial continuity conditions and the traction-free boundary conditions at $x_{2} = - h / 2, h / 2,$ one obtains

\begin{matrix} [\begin{matrix} U (h / 2) \\ V (h / 2) \\ W (h / 2) \\ 0 \\ 0 \\ 0 \end{matrix}] = [\begin{matrix} P_{11} & P_{12} & P_{13} & P_{14} & P_{15} & P_{16} \\ P_{21} & P_{22} & P_{23} & P_{24} & P_{25} & P_{26} \\ P_{31} & P_{32} & P_{33} & P_{34} & P_{35} & P_{36} \\ P_{41} & P_{42} & P_{43} & P_{44} & P_{45} & P_{46} \\ P_{51} & P_{52} & P_{53} & P_{54} & P_{55} & P_{56} \\ P_{61} & P_{62} & P_{63} & P_{64} & P_{65} & P_{66} \end{matrix}] \\ \times [\begin{matrix} U (- h / 2) \\ V (- h / 2) \\ W (- h / 2) \\ 0 \\ 0 \\ 0 \end{matrix}] \end{matrix}

(15)

Thus,

\begin{matrix} | \begin{matrix} P_{41} & P_{42} & P_{43} \\ P_{51} & P_{52} & P_{53} \\ P_{61} & P_{62} & P_{63} \end{matrix} | = 0 \end{matrix}

(16)

Muller's method³⁰ can thus be used to solve equation (16) to obtain a number of eigenvalues $μ_{k}$ , where the above [P] matrix is

[P] = Π_{k = N}^{1} [T_{k}] = [T_{N}] [T_{N - 1}] [T_{N - 2}] \dots [T_{1}]

(17)

in which

[T_{k}] = [M_{k}] [E_{k}] (y_{k + 1} - y_{k}) [M_{k}]^{- 1}

(18)

\begin{matrix} [M]_{(k)} = [\begin{matrix} U_{χ_{1}} & U_{χ_{2}} & \dots . & U_{χ_{6}} \\ V_{χ_{1}} & V_{χ_{2}} & \dots . & V_{χ_{6}} \\ W_{χ_{1}} & W_{χ_{2}} & \dots . & W_{χ_{6}} \\ {S_{12}}_{χ_{1}} & {S_{12}}_{χ_{2}} & \dots . & {S_{12}}_{χ_{6}} \\ {S_{22}}_{χ_{1}} & {S_{22}}_{χ_{2}} & \dots . & {S_{22}}_{χ_{6}} \\ {S_{23}}_{χ_{1}} & {S_{23}}_{χ_{2}} & \dots . & {S_{23}}_{χ_{6}} \end{matrix}]_{(k)} \end{matrix}

(19a)

\begin{matrix} [E]_{(k)} (y) = [\begin{matrix} e^{χ_{1} μ (y)} & 0 & 0 & 0 & 0 & 0 \\ 0 & e^{χ_{2} μ (y)} & 0 & 0 & 0 & 0 \\ 0 & 0 & e^{χ_{3} μ (y)} & 0 & 0 & 0 \\ 0 & 0 & 0 & e^{χ_{4} μ (y)} & 0 & 0 \\ 0 & 0 & 0 & 0 & e^{χ_{5} μ (y)} & 0 \\ 0 & 0 & 0 & 0 & 0 & e^{χ_{6} μ (y)} \end{matrix}]_{(k)} \end{matrix}

(19b)

Equation (19) can be obtained using equation (12) with the material properties and $χ = χ_{j}$ of each layer. The parameter (C_j) $({C} (k) = [\begin{matrix} C_{1} & C_{2} & C_{3} & C_{4} & C_{5} & C_{6} \end{matrix}]_{(k)} T)$ in equation (11) of each layer should be calculated. First, the eigenvector of the bottom layer is generated as follows

\begin{matrix} {Z_{1} (- h / 2)} = [\begin{matrix} U (- h / 2) & V (- h / 2) & W (- h / 2) & 0 & 0 & 0 \end{matrix}]^{T} \end{matrix}

(20)

and the eigenvector of each layer is

{Z_{k} (y_{k})} = [T_{k - 1}] {Z_{k - 1}}

(21)

where

\begin{matrix} U (- h / 2) = | \begin{matrix} P_{42} & P_{43} \\ P_{52} & P_{53} \end{matrix} |, V (- h / 2) = | \begin{matrix} P_{43} & P_{41} \\ P_{53} & P_{51} \end{matrix} |, \\ W (- h / 2) = | \begin{matrix} P_{41} & P_{42} \\ P_{51} & P_{52} \end{matrix} | \end{matrix}

(22)

From ${Z_{k} (y)} = [M_{k}] [E_{k} (y)] {C_{k}}$ , one obtains

{C_{k}} = [M_{k}]^{- 1} [E_{k} (y)]^{- 1} {Z_{k} (y)}

(23)

Finally, the eigensolution is

\begin{matrix} [\begin{matrix} u \\ v \\ w \\ σ 12 \\ \begin{matrix} . \\ . \\ . \end{matrix} \\ σ 33 \end{matrix}] = \sum_{n = 1}^{\infty} D_{n} [\begin{matrix} \sum_{j = 1}^{6} C_{j} U e^{χ_{j} μ_{n} y} \\ \sum_{j = 1}^{6} C_{j} V e^{χ_{j} μ_{n} y} \\ \sum_{j = 1}^{6} C_{j} W e^{χ_{j} μ_{n} y} \\ \sum_{j = 1}^{6} C_{j} S_{12} e^{χ_{j} μ_{n} y} \\ \begin{matrix} . \\ . \\ . \end{matrix} \\ \sum_{j = 1}^{6} C_{j} S_{33} e^{χ_{j} μ_{n} y} \end{matrix}] e^{μ_{n} x} \end{matrix}

(24)

Particular solutions

This section determines a particular solution of equations (3) and (4), which satisfies the interfacial continuity conditions, and the boundary conditions at x₂ = −h/2, h/2. Equations (3) and (4) indicate that the particular solution can be the function only dependent on x₂ as follows

\begin{matrix} [\begin{matrix} u \\ v \\ w \\ σ_{12} \\ σ_{22} \\ σ_{23} \end{matrix}] = ɛ [\begin{matrix} 0 \\ ψ_{p 1} (y) \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}] + b_{1} [\begin{matrix} 0 \\ ψ_{p 2} (y) \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}] + ϑ [\begin{matrix} 0 \\ ψ_{p 3} (y) \\ - xy \\ 0 \\ 0 \\ 0 \end{matrix}] \end{matrix}

(25)

\begin{matrix} [\begin{matrix} σ_{11} \\ σ_{13} \\ σ_{33} \end{matrix}] = (ɛ + b_{1} y) [\begin{matrix} {\overset{\land}{c}}_{13} \\ {\overset{\land}{c}}_{35} \\ {\overset{\land}{c}}_{33} \end{matrix}] - 2 ϑ y [\begin{matrix} {\overset{\land}{c}}_{15} \\ {\overset{\land}{c}}_{55} \\ {\overset{\land}{c}}_{35} \end{matrix}] \end{matrix}

(26)

where

\begin{matrix} ψ_{p 1} (y) = - \sum_{n = 1}^{k - 1} [{\overset{\land}{c}}_{23 (n)} (y_{n} - y_{n - 1})] - {\overset{\land}{c}}_{23 (k)} (y - y_{k - 1}), \\ ψ_{p 2} (y) = - \sum_{n = 1}^{k - 1} [{\overset{\land}{c}}_{23 (n)} (\frac{y_{n}^{2} - y_{n - 1}^{2}}{2})] - {\overset{\land}{c}}_{23 (k)} (\frac{y_{n}^{2} - y_{n - 1}^{2}}{2}), \\ ψ_{p 3} (y) = \sum_{n = 1}^{k - 1} [{\overset{\land}{c}}_{25 (n)} (y_{n}^{2} - y_{n - 1}^{2})] + {\overset{\land}{c}}_{25 (k)} (y_{n}^{2} - y_{n - 1}^{2}) \end{matrix}

Zero eigenvalue solutions

In equation (24), $μ = 0$ is a repeated root, which can be determined by considering the Jordan-chain solution given as follows

\begin{matrix} [\begin{matrix} u \\ v \\ w \\ σ_{21} \\ σ_{22} \\ σ_{23} \end{matrix}]_{μ_{0}} = A ([\begin{matrix} x \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}] + [\begin{matrix} 0 \\ ψ_{01} (y) \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}]) \\ + C ([\begin{matrix} 0 \\ 0 \\ x \\ 0 \\ 0 \\ 0 \end{matrix}] + [\begin{matrix} 0 \\ ψ_{03} (y) \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}]) \\ + B ([\begin{matrix} 0 \\ \frac{x^{2}}{2} \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}] + [\begin{matrix} - xy \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}] + [\begin{matrix} 0 \\ ψ_{02} (y) \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}]) \end{matrix}

(27)

\begin{matrix} [\begin{matrix} σ_{11} \\ σ_{13} \\ σ_{33} \end{matrix}]_{(k)} = A [\begin{matrix} - {\tilde{c}}_{11} \\ - {\tilde{c}}_{15} \\ - {\tilde{c}}_{13} \end{matrix}]_{(k)} + B [\begin{matrix} - {\tilde{c}}_{11} y \\ - {\tilde{c}}_{15} y \\ - {\tilde{c}}_{13} y \end{matrix}]_{(k)} + C [\begin{matrix} {\tilde{c}}_{15} \\ {\tilde{c}}_{55} \\ {\tilde{c}}_{35} \end{matrix}]_{(k)} \end{matrix}

(28)

where

\begin{matrix} ψ_{01} (y) = - \sum_{n = 1}^{k - 1} [{\overset{\land}{c}}_{12 (n)} (y_{n} - y_{n - 1})] - {\overset{\land}{c}}_{12 (k)} (y - y_{k - 1}), \\ ψ_{02} (y) = \sum_{n = 1}^{k - 1} [{\overset{\land}{c}}_{12 (n)} (\frac{y_{n}^{2} - y_{n - 1}^{2}}{2})] + {\overset{\land}{c}}_{12 (k)} (\frac{y^{2} - y_{k - 1}^{2}}{2}), \end{matrix}

\begin{matrix} ψ_{03} (y) = - \sum_{n = 1}^{k - 1} [{\overset{\land}{c}}_{25 (n)} (y_{n} - y_{n - 1})] - {\overset{\land}{c}}_{25 (k)} (y - y_{k - 1}) \end{matrix}

There are three parameters, A, B, and C, in the zero solution, and we use the zero total traction forces and the zero total moment along the x₁ axis at the two edges (x = ±a/2) to find them, as follows

\begin{matrix} {\begin{matrix} A \\ B \\ C \end{matrix}} = [\begin{matrix} \frac{1}{2} \sum_{k = 1}^{N} (y_{k}^{2} - y_{k - 1}^{2}) {\tilde{c}}_{11} & \frac{- 1}{3} \sum_{k = 1}^{N} (y_{k}^{3} - y_{k - 1}^{3}) {\tilde{c}}_{11} & \frac{1}{2} \sum_{k = 1}^{N} (y_{k}^{2} - y_{k - 1}^{2}) {\tilde{c}}_{15} \\ \sum_{k = 1}^{N} (y_{k} - y_{k - 1}) {\tilde{c}}_{11} & \frac{- 1}{2} \sum_{k = 1}^{N} (y_{k}^{2} - y_{k - 1}^{2}) {\tilde{c}}_{11} & \sum_{k = 1}^{N} (y_{k} - y_{k - 1}) {\tilde{c}}_{15} \\ \sum_{k = 1}^{N} (y_{k} - y_{k - 1}) {\tilde{c}}_{15} & \frac{- 1}{2} \sum_{k = 1}^{N} (y_{k}^{2} - y_{k - 1}^{2}) {\tilde{c}}_{15} & \sum_{k = 1}^{N} (y_{k} - y_{k - 1}) {\tilde{c}}_{55} \end{matrix}]^{- 1} {F} \end{matrix}

(29)

\begin{matrix} {F} = - ɛ {\begin{matrix} \frac{1}{2} \sum_{k = 1}^{N} (y_{k}^{2} - y_{k - 1}^{2}) {\tilde{c}}_{13} \\ \sum_{k = 1}^{N} (y_{k} - y_{k - 1}) {\tilde{c}}_{13} \\ \sum_{k = 1}^{N} (y_{k} - y_{k - 1}) {\tilde{c}}_{35} \end{matrix}} - b_{1} {\begin{matrix} \frac{1}{3} \sum_{k = 1}^{N} (y_{k}^{3} - y_{k - 1}^{3}) {\tilde{c}}_{13} \\ \frac{1}{2} \sum_{k = 1}^{N} (y_{k}^{2} - y_{k - 1}^{2}) {\tilde{c}}_{13} \\ \frac{1}{2} \sum_{k = 1}^{N} (y_{k}^{2} - y_{k - 1}^{2}) {\tilde{c}}_{35} \end{matrix}} + 2 ϑ {\begin{matrix} \frac{1}{3} \sum_{k = 1}^{N} (y_{k}^{3} - y_{k - 1}^{3}) {\tilde{c}}_{15} \\ \frac{1}{2} \sum_{k = 1}^{N} (y_{k}^{2} - y_{k - 1}^{2}) {\tilde{c}}_{15} \\ \frac{1}{2} \sum_{k = 1}^{N} (y_{k}^{2} - y_{k - 1}^{2}) {\tilde{c}}_{55} \end{matrix}} \end{matrix}

(30)

Finally, adding the zero solutions (27) and (28) and the particular solutions (25) and (26) to the linear combination of the eigensolution (24) produces the complete solution, where only the parameters D_n are unknown, and they will be calculated using the symplectic orthogonality theory.

Satisfaction of the edge conditions

Symplectic orthogonality in the Hamiltonian state space holds for the eigenvector derived from the state vector, which is composed of the displacement vector and the associated stress vector. With the symplectic orthogonality, given in Zhong,³¹ we are able to determine the unknown constants D_n in the complete solution using symplectic orthogonality, as follows

\sum_{n = 1}^{N} K_{mn} D_{n} = L_{m} (m = 1, 2, 3 \dots, N)

(31)

where

K_{mn} = \sum_{(k) = 1}^{N} \int_{y_{k - 1}}^{y_{k}} (S_{m}^{T} U_{n})_{(k)} dy

(32)

L_{m} = \sum_{(k) = 1}^{N} \int_{y_{k - 1}}^{y_{k}} (U_{m}^{T} σ_{a / 2})_{(k)} dy

(33)

\begin{matrix} σ_{a / 2} (y) = - [\begin{matrix} P_{11} \\ P_{12} \\ P_{13} \end{matrix}] = - [\begin{matrix} σ_{11} \\ σ_{12} \\ σ_{13} \end{matrix}]_{Particularsolution} \\ - [\begin{matrix} σ_{11} \\ σ_{12} \\ σ_{13} \end{matrix}]_{Zerosolution} \end{matrix}

(34)

\begin{matrix} U_{m} = [\begin{matrix} U_{(χ_{1} m)} & U_{(χ_{2} m)} & \dots \\ V_{(χ_{1} m)} & V_{(χ_{2} m)} & \dots \\ W_{(χ_{1} m)} & W_{(χ_{2} m)} & \dots \end{matrix}] \times [\begin{matrix} e^{χ_{1} μ_{m} (y)} & 0 & 0 & 0 & 0 & 0 \\ 0 & e^{χ_{2} μ_{m} (y)} & 0 & 0 & 0 & 0 \\ 0 & 0 & e^{χ_{3} μ_{m} (y)} & 0 & 0 & 0 \\ 0 & 0 & 0 & e^{χ_{4} μ_{m} (y)} & 0 & 0 \\ 0 & 0 & 0 & 0 & e^{χ_{5} μ_{m} (y)} & 0 \\ 0 & 0 & 0 & 0 & 0 & e^{χ_{6} μ_{m} (y)} \end{matrix}] [\begin{matrix} C_{1} \\ C_{2} \\ C_{3} \\ C_{4} \\ C_{5} \\ C_{6} \end{matrix}]_{m} \end{matrix}

(35)

\begin{matrix} S_{m} = [\begin{matrix} {S_{1}}_{(χ_{1} m)} & {S_{1}}_{(χ_{2} m)} & \dots \\ {S_{2}}_{(χ_{1} m)} & {S_{2}}_{(χ_{2} m)} & \dots \\ {S_{3}}_{(χ_{1} m)} & {S_{3}}_{(χ_{2} m)} & \dots \end{matrix}] [\begin{matrix} e^{χ_{1} μ_{m} (y)} & 0 & 0 & 0 & 0 & 0 \\ 0 & e^{χ_{2} μ_{m} (y)} & 0 & 0 & 0 & 0 \\ 0 & 0 & e^{χ_{3} μ_{m} (y)} & 0 & 0 & 0 \\ 0 & 0 & 0 & e^{χ_{4} μ_{m} (y)} & 0 & 0 \\ 0 & 0 & 0 & 0 & e^{χ_{5} μ_{m} (y)} & 0 \\ 0 & 0 & 0 & 0 & 0 & e^{χ_{6} μ_{m} (y)} \end{matrix}] [\begin{matrix} C_{1} \\ C_{2} \\ C_{3} \\ C_{4} \\ C_{5} \\ C_{6} \end{matrix}]_{m} \end{matrix}

(36)

Solution procedures

In this study, we only assume that the bending curvature in the x₂ axis is zero, which can cover most cases of composite plates. For some special cases, this assumption is still appropriate in practice, since the curvature in the plate normal direction is often very small due to the large in-plane stiffness. As such, the symmetric arrangement of the monoclinic laminates along the x₁ axis is not necessary in the proposed method. The applied pure tension force (P), bending moment or torsion may cause the coupling of the tension strain (ɛ), bending curvature (b₁), and twisting curvature (ϑ). Therefore, to find the stress field due to the applied tension force (P), bending moment (M₁), and torsion (M_t), as shown in equations (8) to (10), the following relationships should be used.

\begin{matrix} {\begin{matrix} P \\ M_{1} \\ M_{t} \end{matrix}} = [\begin{matrix} K_{pp} & K_{pb} & K_{pt} \\ K_{bp} & K_{bb} & K_{bt} \\ K_{tp} & K_{tb} & K_{tt} \end{matrix}] {\begin{matrix} ɛ \\ b_{1} \\ ϑ \end{matrix}} or \\ {\begin{matrix} ɛ \\ b_{1} \\ ϑ \end{matrix}} = [\begin{matrix} K_{pp} & K_{pb} & K_{pt} \\ K_{bp} & K_{bb} & K_{bt} \\ K_{tp} & K_{tb} & K_{tt} \end{matrix}]^{- 1} {\begin{matrix} P \\ M_{1} \\ M_{t} \end{matrix}} \end{matrix}

(37)

One can set ɛ, b₁, or ϑ to one unit, respectively, to obtain the above stiffness values using equations (8) to (10).

The known data of the proposed equations are the plate thickness (h), plate width (a), layer thickness (y_k−y_k−1), material properties, uniform tension strain (ɛ), the curvature of bending (b₁), and the curvature of twisting (ϑ). First, equation (16) is used to obtain the eigenvalues ( $(μ_{n}), n = 1, N_{μ}$ ), equations (25) to (30) are then used to find the particular and zero solutions, and finally equation (31) is used to find the parameters (D_n). The displacement and stress fields are obtained using the combination of equations (24) to (28). The integration of equations (32) and (33) can be determined explicitly, so the above procedures only require the solutions of matrix and eigen equations. Two Fortran programs using the above procedures, one solving the eigen-problem and the other generating stress fields, can be found on the web site http://myweb.ncku.edu.tw/∼juju.

Examples and validation

Symmetric angle-ply laminates with tacking sequences of [+45°/–45°/–45°/+45°]

The first example is the symmetric angle-ply laminates with stacking sequences of [45°/–45°/–45°/45°] under applied bending. Since there is no analytic solution for this problem, we use the finite element analysis to validate the proposed theoretical solution. The plate dimensions, material properties, and finite element analysis are the same as those mentioned in the section ‘Illustration of the analyzed monoclinic composite plate’. Four hundred eigenvalues solved by the Muller's method are used in equation (16), where Table 1 shows the first 20 non-dimensional eigenvalues

μ_{i} h

. At the first step, the ɛ, b₁, or ϑ are set to one unit in order to obtain the stiffness of equation (37), as follows

\begin{matrix} {\begin{matrix} P \\ M_{1} \\ M_{t} \end{matrix}} = c_{33} h^{2} [\begin{matrix} 0.68678 & 0 & 0 \\ 0 & 8.837 h^{2} & 1.0859 h^{2} \\ 0 & 1.0859 h^{2} & 2.0827 h^{2} \end{matrix}] {\begin{matrix} ɛ \\ b_{1} \\ ϑ \end{matrix}} \end{matrix}

Table 1.

Eigenvalues of the three examples in the section ‘Examples and validation’.

Number	Example 1 [45°/–45°/–45°/45°]	Example 2 [–45°/45°/–45°/45°]	Example 3 [63.43°/–63.43°/26.56°/–26.56°]
1	0.14215E+01	0.12237E+01	0.14229E+01
2	0.31236E+01	0.30876E+01	0.29734E+01
3	0.37528E+01	0.35716E+01	0.31468E+01
4	0.40468E+01	0.37579E+01	0.37125E+01
5	0.41318E+01	0.41775E+01	0.48305E+01
6	0.44434E+01 + 0.12528E+01i	0.47331E+01 + 0.11245E+01i	0.50818E+01 + 0.14195E+01i
7	0.44434E+01–0.12528E+01i	0.47331E+01–0.11245E+01i	0.50818E+01–0.14195E+01i
8	0.53537E+01	0.52132E+01	0.55659E+01
9	0.62625E+01	0.67787E+01	0.61905E+01
10	0.71318E+01	0.68956E+01	0.74105E+01
11	0.74264E+01	0.74337E+01	0.81959E+01
12	0.77328E+01	0.77318E+01	0.89919E+01
13	0.93796E+01	0.94573E+01	0.93747E+01 + 0.78307E+00i
14	0.97247E+01 + 0.20842E+01i	0.96788E+01 + 0.19694E+01i	0.93747E+01–0.78307E+00i
15	0.97247E+01–0.20842E+01i	0.96788E+01–0.19694E+01i	0.99562E+01
16	0.10953E+02	0.10987E+02	0.11306E+02
17	0.11291E+02	0.11217E+02	0.11651E+02
18	0.11792E+02	0.11935E+02	0.11856E+02
19	0.12628E+02	0.12409E+02	0.12814E+02
20	0.13721E+02 + 0.32518E+01i	0.13762E+02	0.13462E+02 + 0.31555E+01i

The above equation indicates that the tension force produces only the tension strain, but the bending or torsional moment will produce the deformation coupled with the bending and twisting curvatures, although the plate is the symmetric arrangement of angle-ply laminates. Figure 4 shows the non-dimensional stress fields along the x₁ axis at the center of the fourth layer (x₂ = 3h/8) using the proposed theoretical method and the finite element method under M₁ = 1 Nmm. Figure 5 shows the non-dimensional stress fields along the x₂ axis for x₁ = 0.4625a using the above two methods. It can be seen that they are in excellent agreement, even though the stress curves are rather uneven along the x₁ or x₂ axis.

Figure 4.

Comparison of the stress fields obtained from the theoretical and finite element results along the x₁ axis (x₂=3h/8) for the laminates with stacking sequences of [45°/–45°/–45°/45°] under the applied bending of 1 Nmm.

Figure 5.

Comparison of the stress fields obtained from the theoretical and finite element results along the x₂ axis (x₁=0.4625a) for the laminates with stacking sequences of [45°/–45°/–45°/45°] under the applied bending of 1 Nmm.

Angle-ply laminates with tacking sequences of [–45°/+45°/–45°/+45°]

The second example is another type of symmetric angle-ply laminate with stacking sequences of [–45°/+45°/–45°/+45°], and other conditions are the same as those of the first example. Four hundred eigenvalues solved using Muller's method are used in equation (16), where Table 1 shows the first 20 non-dimensional eigenvalues $μ_{i} h$ . At the first step, the ɛ, b₁, or ϑ are set to one unit to obtain the stiffness of equation (37), as follows

\begin{matrix} {\begin{matrix} P \\ M_{1} \\ M_{t} \end{matrix}} = c_{33} h^{2} [\begin{matrix} 0.69727 & 0 & 0.13347 h \\ 0 & 8.38975 h^{2} & 0 \\ 0.13347 h & 0 & 1.5504 h^{2} \end{matrix}] {\begin{matrix} ɛ \\ b_{1} \\ ϑ \end{matrix}} \end{matrix}

The above equation indicates that the applied tension or torque produces the deformation coupled with the tension strain and twisting curvature, but the bending moment will only produce the bending curvature deformation. Figure 6 shows the non-dimensional stress fields along the x₁ axis at the center of the fourth layer (x₂ = 3h/8) using the proposed theoretical method and the finite element method under M_t = 1 Nmm. Figure 7 shows the non-dimensional stress fields along the x₂ axis for x₁ = 0.925a using the above two methods. Similar to the first example, the two results are also in excellent agreement.

Figure 6.

Comparison of the stress fields obtained from the theoretical and finite element results along the x₁ axis (x₂=3h/8) for the laminates with stacking sequences of [–45°/+45°/–45°/+45°] under the applied torque of 1 Nmm.

Figure 7.

Comparison of the stress fields obtained from the theoretical and finite element results along the x₂ axis (x₁=0.4625a) for the laminates with stacking sequences of [–45°/+45°/–45°/+45°] under the applied torque of 1 Nmm.

Angle-ply laminates with stacking sequences of [+63.43°/–63.43°/+26.56°/–26.56°]

The third example is an unsymmetrical arrangement of angle-ply laminates with stacking sequences of [+63.43°/–63.43°/+26.56°/–26.56°], and the other conditions are the same as those in the first example. Four hundred eigenvalues solved by Muller's method are used in equation (16), where Table 1 shows the first 20 non-dimensional eigenvalues $μ_{i} h$ . At the first step, the ɛ, b₁, or ϑ are set to one unit to obtain the stiffness of equation (37), as follows

\begin{matrix} {\begin{matrix} P \\ M_{1} \\ M_{t} \end{matrix}} = c_{33} h^{2} [\begin{matrix} 0.50897 & 1.0309 h & 0.012627 h \\ 1.0309 h & 7.0798 h^{2} & - 1.4834 h^{2} \\ 0.012627 h & - 1.4834 h^{2} & 2.5812 h^{2} \end{matrix}] {\begin{matrix} ɛ \\ b_{1} \\ ϑ \end{matrix}} \end{matrix}

The above equation indicates that the applied tension, bending moment or torque produces the deformation coupled with the tension strain, bending curvature, and twisting curvature. Figure 8 shows the non-dimensional stress fields along the x₁ axis at the center of the fourth layer (x₂ = 3h/8) using the proposed theoretical method and the finite element method under P = 1 N. Figure 9 shows the non-dimensional stress fields along the x₂ axis for x₁ = 0.925a using the above two methods. Similar to the first example, the two results are also in excellent agreement.

Figure 8.

Comparison of the stress fields obtained from the theoretical and finite element results along the x₁ axis (x₂=3h/8) for the laminates with stacking sequences of [+63.43°/–63.43°/+26.56°/–26.56°] under the applied tension of 1 N.

Figure 9.

Comparison of the stress fields obtained from the theoretical and finite element results along the x₂ axis (x₁=0.4625a) for the laminates with stacking sequences of [+63.43^o/–63.43^o/+26.56^o/–26.56^o] under the applied tension of 1 N.

Conclusion

Based on the Hamiltonian state space formalism, the stress fields in composite plates with monoclinic angle-ply laminates under extension, bending, and torsion are investigated. The Hamiltonian characteristics enables us to use the eigenfunction expansion to analyze the problem of multi-layered composite laminates efficiently, regardless of the number of layers, in which the stress-free boundary condition along the plate's upper and bottom surface is fitted. The symplectic orthogonality is then used systematically to fit the stress-free boundary condition along the left and right sides of the plate. Finally, a matrix equation is obtained to find the applied tension strain, bending curvature, and torsional curvature from the applied tension, bending, and torque. One can thus solve the stress field of composite plates using the force or displacement control. The integration of the proposed method can be determined explicitly, and the numerical procedures only require the solutions of matrix and eigen equations. The material properties of the laminate layers do not necessarily need to be symmetrically arranged, and the only limitation is that the out-of-plane curvature needs to be zero in the proposed method. The finite element analyses are used to validate the accuracy of the theoretical solution, and three examples including symmetric and unsymmetrical arrangements of the composite layers show that the solutions of the theoretical and finite element methods are in excellent agreements.

The laminates problems solved by the finite element method require a very fine mesh, because each laminates layer is quite thin compared to the cross section. Moreover, high-order elements should be used to avoid excessive errors for the bending effect. For example, each finite element analysis of the three examples in this paper required 12 hours to obtain the solution using an Intel I7-4790 computer. The solution scheme is an efficient parallel conjugate gradient method with the SSOR pre-conditioning. For the proposed analytical method, most of the computer time is used to solve the eigen-problem, where the efficient Muller's method was used. When the eigenvalues are obtained, one can solve the problem for different tensions, bending moments, and torsions without having to solve the eigenvalues again. For the three examples in this paper, we solved 1000 complex eigenvalues using the computer time of 1.3 hours. Thus, for laminates problems, the proposed analytic method can be much fast than the finite element method. It is worth noting that analytical solutions are often milestone solutions that can be used to verify the accuracy of other methods, such as finite element and experimental methods, while the computer time of the analytical method is not a serious issue. Nevertheless, the computational efficiency of the proposed method is still superior to that of the finite element method.

Footnotes

Declaration of Conflicting Interests

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

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Rogers

Watson

Spencer

AJM

. An exact three-dimensional solution for normal loading of inhomogeneous and laminated anisotropic elastic plates of moderate thickness. Proc R Soc Lond A 1992; 437: 199–213.

Tarn

Huang

. Saint-Venant end effects in multilayered piezoelectric laminates. Int J Solids Struct 2002; 39: 4979–4998.

Pagano

. Exact solutions for composite laminates in cylindrical bending. J Compos Mater 1969; 3: 398–411.

Ren

. Exact solutions for laminated cylindrical shells in cylindrical bending. Compos Sci Technol 1987; 29: 169–187.

Savoia

Reddy

. A variational approach to three-dimensional elasticity solutions of laminated composite plates. J Appl Mech 1992; 59: S166–S175.

Heyliger

. Exact solutions for simply supported laminated piezoelectric plates. ASME J Appl Mech 1997; 64: 299–306.

Pan

. Exact solution for simply supported and multilayered magneto-elactro-elastic plates. J Appl Mech 2001; 68: 608–618.

Demasi

. Three-dimensional closed form solutions and theories for orthotropic plates. Mech Adv Mater Struc 2010; 17: 20–39.

Yang

Gao

Pan

, et al. An exact closed-form solution for a multilayered one-dimensional orthorhombic quasicrystal plate. Acta Mech 2015; 226: 3611–3621.

10.

Mittelstedt

Becker

. Free-edge effects in composite laminates. ASME Appl Mech Rev 2007; 60: 217–245.

11.

Cho

Kim

. Iterative free-edge stress analysis of composite laminates under extension, bending, twisting and thermal loadings. Int J Solids Struct 2000; 37: 435–459.

12.

Nosier

Maleki

. Free-edge stresses in general composite laminates. Int J Mech Sci 2008; 50: 1435–1447.

13.

Sarvestani

. Free-edge stress analysis of general composite laminates under extension, torsion and bending. Appl Math Model 2012; 36: 1570–1588.

14.

Saeedi

Sab

Caron

. Stress analysis of long multilayered plates subjected to invariant loading: analytical solutions by a layerwise stress model. Compos Struct 2013; 100: 307–322.

15.

Tarn

Chang

. A Hamiltonian state space approach to anisotropic elasticity and piezoelasticity. Acta Mech 2013; 224: 1271–1284.

16.

Andakhshideh

Tahani

. Free-edge stress analysis of general rectangular composite laminates under bending, torsion and thermal loads. Eur J Mech A-Solid 2013; 42: 229–240.

17.

Dhanesh

Kapuria

Achary

. Accurate prediction of three-dimensional free edge stress field in composite laminates using mixed-field multi-term extended Kantorovich method. Acta Mech 2017; 228: 2895–2919.

18.

Han

Niu

. A 3D exact analysis of the boundary layer effect of asymmetric piezoelectric laminates with electromechanical coupling. Int J Solids Struct 2015; 72: 118–129.

19.

Alipour

. A novel economical analytical method for bending and stress analysis of functionally graded sandwich circular plates with general elastic edge conditions, subjected to various loads. Compos Part B – Eng 2016; 95: 48–63.

20.

Liang

Tarn

. Analytic determination of stress fields in cross-ply symmetric composite laminates. Int J Solids Struct 2016; 94–95: 87–99.

21.

Schnabel

Yousfi

Mittelstedt

. Free-edge stress fields in cylindrically curved symmetric and unsymmetric cross-ply laminates under bending load. Compos Struct 2017; 180: 862–875.

22.

Hajikazemi

Van Paepegem

. A variational model for free-edge interlaminar stress analysis in general symmetric and thin-ply composite laminates. Compos Struct 2017; 184: 443–451.

23.

Vukasović

Pavazza

Vlak

. Analytic solution for torsion of thin-walled laminated composite beams of symmetrical open cross sections with influence of shear. Arch Appl Mech 2017; 87: 1371–1384.

24.

Yan

Pagani

Carrera

, et al. Exact solutions for the macro-, meso- and micro-scale analysis of composite laminates and sandwich structures. J Compos Mater 2018; 52: 3109–3124.

25.

Wackerfuß

Kroker

. An efficient semi-analytical simulation framework to analyse laminated prismatic thin-walled beams. Comput Struct 2018; 208: 32–50.

26.

Cas

Bratina

Saje

, et al. Non-linear analysis of composite steel-concrete beams with incomplete interaction. Steel Compos Struct 2004; 4: 489–507.

27.

Tarn

. A state space formalism for anisotropic elasticity Part I: rectilinear anisotropy. Int J Solids Struct 2002; 39: 5143–5155.

28.

Lekhnitskii

. Theory of elasticity of an anisotropic body, Moscow: Mir, 1981.

29.

Schnabl

Saje

Turk

, et al. Locking-free two-layer Timoshenko beam element with interlayer slip. Finite Elem Anal Des 2007; 43: 705–714.

30.

Press

Flannery

Teukolsky

, et al. Numerical recipes, the art of scientific computing, New York: Cambridge University Press, 1986.

31.

Zhong WX. A new systematic methodology for theory of elasticity. Dalian (in Chinese): Dalian University of Technology Press, 1995.