Nonlinear Dynamics Suppression of Bending Deflection in a Composite Laminated Plate Using a Beam Stiffener Due to Critical Buckling Load and Shear Force

Abstract

In this paper, the nonlinear dynamics behavior of the bending deflection of a stiffened composite laminated plate is suppressed using beam stiffeners at different fiber volume fractions and different aspect ratios. The non-periodic motion and chaos in a stiffened composite laminated plate is detected using the largest Lyapunov exponent parameter and power density function of a fast Fourier transform (FFT). The critical buckling load is calculated at different thickness ratios, numbers of stiffeners, lamination angles and stiffener–depth ratios based on different boundary conditions. The nonlinear response of the bending deflection is analyzed analytically, numerically and experimentally. The analytic solution has been derived using Levy and Navier solutions of classical laminate plate theory at different boundary conditions (CLPT). The numerical simulation was conducted using the ANSYS program while the experiment test was carried out using a strain gauge through a strain meter device. Experimentally, a Southwell plot is used to investigate the value of the critical buckling load. The combined loading are the in-plane compression mechanical load and shear force. All the values of the largest Lyapunov exponent are positive, which gives indication to non-periodic motion and chaos. The nonlinear dynamics behavior of the bending deflection is decreased with the increasing of number of stiffeners in which the value of largest Lyapunov exponent has been decreased. The nonlinear dynamics behavior is increased with the increasing of aspect ratios and fiber volume fractions. The system with an aspect ratio (2.5) and fiber volume fraction (${\mathsf{\upsilon }}_{\mathrm{f}}$ = 80%) for an un-stiffened plate is more chaotic than the other systems.

1. Introduction

The stiffened composite laminated plate is widely used in steam boilers with longitudinal tubes, bodies’ ships and submarines’ space crafts by taking into consideration the high stiffness and the ratio between the high strength and low weight. The nonlinear dynamic phenomenon in a stiffened plate due to the bending deflection through the plate thickness could be either a periodic motion or non-periodic motion and chaos. Keshav and Patel studied the non-linear dynamic buckling behavior of laminated composite curved panels when the plate is subjected to a rectangular pulse load at various amplitude and durations. They used finite element analysis to study the effects of aspect ratio, the radius of curvature and the thickness ratio on the nonlinear dynamics phenomenon of composite plate [1]. Mondal and Ramachandra used B-spline function in finite element analysis to analyze the unidirectional in-plane compression pulse loading [2]. They studied the effect of nonlinear dynamic pulse buckling such as sinusoidal, exponential and rectangular on the shock spectrum using the Tsai–Wu quadratic interaction criterion. A nonlinear finite element method is presented by Taczała et al. to investigate the nonlinear stability of stiffened functionally graded materials (FGM) plates [3]. They studied the influence of material properties, geometrical properties of stiffeners and initial deflections on the buckling and post-buckling responses of the stiffened plates. Zhang et al. applied Hamilton’s principle to present the solution of the dynamic buckling of stiffened plates under in-plane impact loading [4]. They used the Galerkin method with the aid of Fourier series to determine the deflection of the plate. Yousuf studied the nonlinear dynamics behavior of the bending deflection for the composite laminated plate using the conception of the largest Lyapunov exponent. A power spectrum analysis has been conducted using the amplitude of Fast Fourier Transform (FFT) to detect the non-periodic motion of the bending deflection [5,6]. Hegaze presented the dynamic analysis of a stiffened and unstiffened composite-laminated plate using a high-order finite element [7]. The shear effect had been taken into consideration in the derivation of the deflection to avoid the problem of un-symmetry in the stiffness matrix. Less and Abramovich solved the nonlinear dynamics equation numerically for different loading amplitudes and various loading durations to obtain the buckling response [8]. Aboudi and Gilat performed the stability of infinitely wide composite plates under the effect of suddenly applied thermal and mechanical loading. They showed an efficient tool of the dynamic buckling in composite plates under various type of loading [9]. Chai et al. studied the nonlinear vibration of composite laminated plate against time by solving the dynamic equation of motion using von Kármán nonlinear plate theory. They analyzed the amplitude–frequency curve of the displacement against time at different ply angles of the composite laminated plate [10]. Borkowski applied dynamic tools such as phase portraits, Poincaré maps, a fast Fourier transform (FFT) analysis and largest Lyapunov exponents on the isotropic plate. He used the vibration theory of the dynamic to describe the theory of bifurcation and chaos in terms of its dynamic stability and instability regions [11]. Han-Gyu studied the nonlinear dynamic behavior of post-buckled composite plates experimentally and theoretically. Experimentally he applied harmonic loads by using an electro-dynamic shaker to measure the dynamic response of the plate specimens based on the single-point displacement-sensing laser approach. The chaotic and periodic motions behavior of the plate response were investigated using classical laminate plate theory and a finite element model [12]. In this paper, the composite plate is subjected to an in-plane compression mechanical load and shear force. The nonlinear dynamics phenomenon is suppressed using beam stiffeners at different fiber volume fractions and aspect ratios. The nonlinear dynamics behavior of the bending deflection through the plate thickness for the stiffened composite laminated plate is detected using the largest Lyapunov exponent parameter and power density function of fast Fourier transform (FFT). The manuscript is organized as in the follower steps: (a) Section 2 discusses the experiment test to track the bending deflection in the direction of plate thickness using a strain gauge through a strain meter device. The critical buckling load is calculated using a Southwell plot. (b) Section 3 discusses the analytical solution of the bending deflection and critical buckling load for the (S-F-S-F) boundary condition using the Levy solution. (c) Section 4 discusses the analytic solution of the bending deflection and critical buckling load for the (S-S-S-S) boundary condition using the Navier solution. (d) Section 5 discusses the calculation of the constants that are needed in the solution of the bending deflection and critical buckling load at different boundary conditions. (e) Section 6 discusses the numerical simulation of the bending deflection and critical buckling load using the ANSYS program for an un-stiffened plate and plate with stiffeners. (f) Section 7 discusses the steps how to find the bending deflection and critical buckling load using the ANSYS program. (g) Section 8 and Section 9 discuss how to detect the non-periodic motion of the bending deflection in un-stiffened and stiffened composite laminated plate using the largest Lyapunov exponent parameter and power density function of the fast Fourier transform (FFT).

2. Experiment Test

In this study, the fabrication of composite material is consisted of plane E-glass fiber plies in a thermosetting polyester resin. The composite laminate plate is manufactured by hand in the lay-up technique. The fabrication of composite laminates is conducted from six layers of E-glass fabric sheets. The fabric sheets are placed in a (24 cm × 24 cm) aluminum mold. Two X-ray photo sheets are used to avoid abrasion in a specimen surface. The fiber volume fraction is determined experimentally for the composites plate from the following relationship [13]:

(1)${\mathsf{\upsilon }}_{\mathrm{f}}=\frac{1}{1+\frac{1-\mathsf{\phi }}{\mathsf{\phi }}\ast \frac{{\mathsf{\rho }}_{\mathrm{f}}}{{\mathsf{\rho }}_{\mathrm{m}}}}$

where:

φ: is the fiber weight fraction.

And

$\mathsf{\phi }=\frac{{\mathrm{m}}_{\mathrm{f}}}{{\mathrm{m}}_{\mathrm{C}}}$

When the composite with the dimensions (24 × 24) gets dried and gets out from the mold, we can weigh the composite plate to determine (${\mathrm{m}}_{\mathrm{C}}$). Before fabrication we can weigh the E-glass fiber to find the value of (${\mathrm{m}}_{\mathrm{f}}$). Divide (${\mathrm{m}}_{\mathrm{f}}$) on (${\mathrm{m}}_{\mathrm{C}}$) to calculate (ψ). By substituting all these values in Equation (1) we can find the value of (${\mathsf{\upsilon }}_{\mathrm{f}}$) as shown in Figure 1.

The composite laminated plate with length (200 mm), width (100 mm) and thickness (4 mm) is selected in the experiment test to calculate critical buckling load. One and two stiffeners of width (6 mm) and depth (4 mm) are selected along the length of the main composite plate. The lamination angle for all stiffeners is equal to (0) degrees, while the lamination angle of composite plate is assumed to be cross ply. The dimensions of the specimen that were used in the buckling test has been taken from ASTM E1876 [14]. A tensile test machine was used in a buckling test with an axial compression load of (200 KN) in the vertical direction, as illustrated in Figure 2.

The boundary conditions of the specimen is simply supported by two ends and free at the other two ends (S-F-S-F). The simply supported boundary condition is simulated along the top and bottom edges of the compression test device. The specimen is compressed slowly with a constant cross-head speed of (5 mm/min) until buckling. Figure 3, Figure 4 and Figure 5 show the buckling test for the un-stiffened plate, stiffened plate with one stiffener and stiffened plate with two stiffeners.

For the buckling test loading, the test specimen is placed between two extremely stiff machine heads. The lower head is fixed during the test, whereas the upper one is moved downwards by a servo hydraulic cylinder. The strain gauge is placed in the middle point of the composite laminated plate through the z-axis at the connection point between the composite plate and the stiffener. The strain meter is used to catch the voltage of the bending deflection through the plate thickness due to the in-plane compression mechanical load and the shear force. When the upper head of the servo hydraulic cylinder moves downward, the composite laminated plate will buckle under the effect of the compression load. The data of the bending deflection through the plate thickness against time is stocked in an Excel spreadsheet, and these data were processed using MatLab software. The deflection through the plate thickness is calculated using the following equation:

${\mathsf{\delta }}_{\mathrm{zz}}={\mathsf{\epsilon }}_{\mathrm{zz}}\ast \mathrm{h}$

where:

${\mathsf{\epsilon }}_{\mathrm{zz}}$: is the experiment strain through the composite plate thickness.

${\mathsf{\delta }}_{\mathrm{zz}}$: is the experiment bending deflection through the composite plate thickness.

h: is the composite plate thickness.

The Southwell plot is a histogram between the bending deflection (${\mathsf{\delta }}_{\mathrm{zz}}$) and the in-plane compression mechanical load [15]. The in-plane compression mechanical load is readable from the buckling test machine screen since it has the values 100 N, 150 N, 200 N, 250 N, 300 N, 350 N, 400 N, 450 N and 500 N. The relationship between the bending deflection (${\mathsf{\delta }}_{\mathrm{zz}})$ and in-plane compression mechanical load is a straight line in which the slop represents the critical buckling load. The Southwell plot is shown in Figure 6 and Figure 7. It can be noticed from Figure 7 that the bending deflection of the un-stiffened and stiffened plates is increased with the increasing of the in-plane compression mechanical load. Moreover, the bending deflection is decreased with the use of beam stiffeners.

3. Analytic Solution of Bending Deflection

The Levy solution is used in the derivation of bending deflection through plate thickness in which the Levy solution applies to any type of boundary conditions under the effect of in-plane compression mechanical load and shear force. Based on the classical laminate plate theory, the stress–strain relationship is:

(2)${\left\{\begin{array}{c}{\mathsf{\sigma }}_{\mathrm{xx}}\\ {\mathsf{\sigma }}_{\mathrm{yy}}\\ {\mathsf{\sigma }}_{\mathrm{xy}}\end{array}\right\}}^{\mathrm{k}}=\left[\begin{array}{ccc}{\overline{\mathrm{Q}}}_{11}& {\overline{\mathrm{Q}}}_{12}& {\overline{\mathrm{Q}}}_{16}\\ {\overline{\mathrm{Q}}}_{12}& {\overline{\mathrm{Q}}}_{22}& {\overline{\mathrm{Q}}}_{26}\\ {\overline{\mathrm{Q}}}_{16}& {\overline{\mathrm{Q}}}_{26}& {\overline{\mathrm{Q}}}_{66}\end{array}\right]\ast \left\{\begin{array}{c}{\mathsf{\epsilon }}_{\mathrm{xx}}\\ {\mathsf{\epsilon }}_{\mathrm{yy}}\\ {\mathsf{\gamma }}_{\mathrm{xy}}\end{array}\right\}$

where:

$\left\{\begin{array}{c}{\mathsf{\epsilon }}_{\mathrm{xx}}\\ {\mathsf{\epsilon }}_{\mathrm{yy}}\\ {\mathsf{\gamma }}_{\mathrm{xy}}\end{array}\right\}=\left\{\begin{array}{c}{\mathsf{\epsilon }}_{\mathrm{xx}}^{\left(0\right)}\\ {\mathsf{\epsilon }}_{\mathrm{yy}}^{\left(0\right)}\\ {\mathsf{\gamma }}_{\mathrm{xy}}^{\left(0\right)}\end{array}\right\}+\mathrm{z}\ast \left\{\begin{array}{c}{\mathsf{\epsilon }}_{\mathrm{xx}}^{\left(1\right)}\\ {\mathsf{\epsilon }}_{\mathrm{yy}}^{\left(1\right)}\\ {\mathsf{\gamma }}_{\mathrm{xy}}^{\left(1\right)}\end{array}\right\}$

Additionally,

$\left\{{\mathsf{\epsilon }}^{\left(0\right)}\right\}=\left\{\begin{array}{c}{\mathsf{\epsilon }}_{\mathrm{xx}}^{\left(0\right)}\\ {\mathsf{\epsilon }}_{\mathrm{yy}}^{\left(0\right)}\\ {\mathsf{\gamma }}_{\mathrm{xy}}^{\left(0\right)}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial {\mathrm{u}}_{\mathrm{o}}}{\partial \mathrm{x}}\\ \frac{\partial {\mathrm{v}}_{\mathrm{o}}}{\partial \mathrm{y}}\\ \frac{\partial {\mathrm{u}}_{\mathrm{o}}}{\partial \mathrm{y}}+\frac{\partial {\mathrm{v}}_{\mathrm{o}}}{\partial \mathrm{x}}\end{array}\right\}$

$\left\{{\mathsf{\epsilon }}^{\left(1\right)}\right\}=\left\{\begin{array}{c}{\mathsf{\epsilon }}_{\mathrm{xx}}^{\left(1\right)}\\ {\mathsf{\epsilon }}_{\mathrm{yy}}^{\left(1\right)}\\ {\mathsf{\gamma }}_{\mathrm{xy}}^{\left(1\right)}\end{array}\right\}=\left\{\begin{array}{c}-\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2}\\ -\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^2}\\ -2\ast \frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial \mathrm{x}\partial \mathrm{y}}\end{array}\right\}$

The plate equation in terms of bending moments is

(3)$\frac{{\partial }^2{\mathrm{M}}_{\mathrm{xx}}}{\partial {\mathrm{x}}^2}+2\ast \frac{{\partial }^2{\mathrm{M}}_{\mathrm{xy}}}{\partial \mathrm{x}\partial \mathrm{y}}+\frac{{\partial }^2{\mathrm{M}}_{\mathrm{yy}}}{\partial {\mathrm{y}}^2}={\mathrm{I}}_0\ast \frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{t}}^2}$

The bending moments of the rectangular plate are

(4)$\left\{\begin{array}{c}{\mathrm{M}}_{\mathrm{xx}}\\ {\mathrm{M}}_{\mathrm{yy}}\\ {\mathrm{M}}_{\mathrm{xy}}\end{array}\right\}= \left[\begin{array}{ccc}{\mathrm{B}}_{11}& {\mathrm{B}}_{12}& 0\\ {\mathrm{B}}_{12}& {\mathrm{B}}_{22}& 0\\ 0& 0& {\mathrm{B}}_{66}\end{array}\right]\ast \left\{\begin{array}{c}\frac{\partial {\mathrm{u}}_{\mathrm{o}}}{\partial \mathrm{x}}\\ \frac{\partial {\mathrm{v}}_{\mathrm{o}}}{\partial \mathrm{y}}\\ \frac{\partial {\mathrm{u}}_{\mathrm{o}}}{\partial \mathrm{y}}+\frac{\partial {\mathrm{v}}_{\mathrm{o}}}{\partial \mathrm{x}}\end{array}\right\}-\left[\begin{array}{ccc}{\mathrm{D}}_{11}& {\mathrm{D}}_{12}& 0\\ {\mathrm{D}}_{12}& {\mathrm{D}}_{22}& 0\\ 0& 0& {\mathrm{D}}_{66}\end{array}\right]\ast \left\{\begin{array}{c}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2}\\ \frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^2}\\ 2\ast \frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial \mathrm{x}\partial \mathrm{y}}\end{array}\right\}$

where

${\mathrm{B}}_{\mathrm{ij}}=\frac{1}{2}\ast {\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{N}}{\overline{\mathrm{Q}}}_{\mathrm{ij}}{}^{\left(\mathrm{k}\right)}\ast ({\mathrm{z}}_{\mathrm{k}+1}^2-{\mathrm{z}}_{\mathrm{k}}^2)=0$, since the stacking sequences of plate orientation are the same above and below the z-axis, and this case is called special orthotropic.

${\mathrm{D}}_{\mathrm{ij}}=\frac{1}{3}\ast {\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{N}}{\overline{\mathrm{Q}}}_{\mathrm{ij}}{}^{\left(\mathrm{k}\right)}\ast \left({\mathrm{z}}_{\mathrm{k}+1}^3-{\mathrm{z}}_{\mathrm{k}}^3\right)$

$\left\{{\mathrm{I}}_{\mathrm{o}}\right\}={{\displaystyle \int }}_{-\frac{\mathrm{h}}{2}}^{\frac{\mathrm{h}}{2}}{\mathsf{\rho }}_{\mathrm{o} }\mathrm{dz}$

Substitute Equation (4) into Equation (3) to obtain the plate equation in terms of the bending stiffness:

(5)${\mathrm{D}}_{11}\ast \left(\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^4}\right)+2\ast \left({\mathrm{D}}_{12}+2{\ast \mathrm{D}}_{66}\right)\ast \left(\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2\partial {\mathrm{y}}^2}\right)+{\mathrm{D}}_{22}\ast \left(\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^4}\right)+{\mathsf{\rho }}_{\mathrm{o} }\ast \frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{t}}^2}+{\mathrm{N}}_{\mathrm{xx}}\ast \frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2}=0$

The composite laminated plate with stiffeners is assumed to be orthotropic. The distribution of the beam stiffeners the along (x and y) axis with the dimension is shown in Figure 8.

The bending stiffness of the composite laminated plate in the presence of the beam stiffeners are:

When the stiffeners are distributed along y-axis,

${ \mathrm{D}}_{11}=\frac{{\mathrm{E}}_{11}{\mathrm{h}}^3}{12\left(1-{\mathsf{\nu }}_{12}{\mathsf{\nu }}_{21}\right)}+\frac{{\mathrm{E}}_{11}}{{\mathrm{b}}_{\mathrm{x}}}\left(\frac{{\mathrm{t}}_{\mathrm{x}}\ast {\mathrm{h}}_{\mathrm{x}}^3}{12}\right)$

When the stiffeners are distributed along x-axis,

${ \mathrm{D}}_{22}=\frac{{\mathrm{E}}_{22}{\mathrm{h}}^3}{12\left(1-{\mathsf{\nu }}_{12}{\mathsf{\nu }}_{21}\right)}+\frac{{\mathrm{E}}_{22}}{{\mathrm{b}}_{\mathrm{y}}}\left(\frac{{\mathrm{t}}_{\mathrm{y}}\ast {\mathrm{h}}_{\mathrm{y}}^3}{12}\right)$

${ \mathrm{D}}_{12}={\mathsf{\nu }}_{12}{\mathrm{D}}_{22}={\mathsf{\nu }}_{21}{\mathrm{D}}_{11}$

${\mathrm{D}}_{66}=\frac{\left(1-\sqrt{{\mathsf{\nu }}_{12}{\mathsf{\nu }}_{21}}\right)}{2}\sqrt{{\mathrm{D}}_{11}{\mathrm{D}}_{22}}$

The solution of Equation (5) is

(6)${\mathrm{w}}_{\mathrm{o}}{\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)}_{\mathrm{H}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }{\mathrm{Y}}_{\mathrm{m}}\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

where

${\mathsf{\alpha }}_{\mathrm{m}}=\frac{\mathrm{m}\ast \mathsf{\pi }}{\mathrm{a}}$

Substitute Equation (6) into Equation (5) to obtain the general solution of the plate equation:

(7)${\mathrm{Y}}_{\mathrm{m}}^{″″}-\frac{2\ast \left({\mathrm{D}}_{12}+2{\ast \mathrm{D}}_{66}\right)}{{\mathrm{D}}_{22}}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2{\ast \mathrm{Y}}_{\mathrm{m}}^{″}+\left(\frac{{\mathrm{D}}_{11}}{{\mathrm{D}}_{22}}{\ast \mathsf{\alpha }}_{\mathrm{m}}^4-\frac{{\mathsf{\rho }}_{\mathrm{o} }{\ast \mathsf{\omega }}_{\mathrm{mn}}^2}{{\mathrm{D}}_{22}}-\frac{{\mathrm{N}}_{\mathrm{xx}}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2}{{\mathrm{D}}_{22}}\right){\ast \mathrm{Y}}_{\mathrm{m}}=0$

The roots of Equation (7) are

$\mathsf{\alpha }=\mp \sqrt{\frac{\left({\mathrm{D}}_{12}+2{\ast \mathrm{D}}_{66}\right){\ast \mathsf{\alpha }}_{\mathrm{m}}^2}{{\mathrm{D}}_{22}}+\sqrt{\frac{{\left({\mathrm{D}}_{12}+2{\ast \mathrm{D}}_{66}\right)}^2{\ast \mathsf{\alpha }}_{\mathrm{m}}^4}{{\mathrm{D}}_{22}^2}-\frac{{\mathrm{D}}_{11}{\ast \mathsf{\alpha }}_{\mathrm{m}}^4}{{\mathrm{D}}_{22}}+\frac{{\mathsf{\rho }}_{\mathrm{o} }{\ast \mathsf{\omega }}_{\mathrm{mn}}^2}{{\mathrm{D}}_{22}}+\frac{{\mathrm{N}}_{\mathrm{xx}}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2}{{\mathrm{D}}_{22}}}}$

$\mathsf{\beta }=\mp \mathrm{i}\sqrt{\frac{-\left({\mathrm{D}}_{12}+2{\ast \mathrm{D}}_{66}\right){\ast \mathsf{\alpha }}_{\mathrm{m}}^2}{{\mathrm{D}}_{22}}+\sqrt{\frac{{\left({\mathrm{D}}_{12}+2{\ast \mathrm{D}}_{66}\right)}^2{\ast \mathsf{\alpha }}_{\mathrm{m}}^4}{{\mathrm{D}}_{22}^2}-\frac{{\mathrm{D}}_{11}{\ast \mathsf{\alpha }}_{\mathrm{m}}^4}{{\mathrm{D}}_{22}}+\frac{{\mathsf{\rho }}_{\mathrm{o} }{\ast \mathsf{\omega }}_{\mathrm{mn}}^2}{{\mathrm{D}}_{22}}+\frac{{\mathrm{N}}_{\mathrm{xx}}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2}{{\mathrm{D}}_{22}}}}$

Then

(8)${\mathrm{Y}}_{\mathrm{m}}\left(\mathrm{y}\right)={\mathrm{N}}_1\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{y}\right)+{\mathrm{N}}_2\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{y}\right)+{\mathrm{N}}_3\ast \cos \left(\mathsf{\beta }\ast \mathrm{y}\right)+{\mathrm{N}}_4\ast \sin \left(\mathsf{\beta }\ast \mathrm{y}\right)$

Additionally,

(9)${\mathrm{w}}_{\mathrm{o}}{\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)}_{\mathrm{H}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }[{\mathrm{N}}_1\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{y}\right)+{\mathrm{N}}_2\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{y}\right)+{\mathrm{N}}_3\ast \cos \left(\mathsf{\beta }\ast \mathrm{y}\right)+{\mathrm{N}}_4\ast \sin \left(\mathsf{\beta }\ast \mathrm{y}\right)]\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

Figure 8 shows the loading and the boundary condition. The Levy solution is assumed that the two sides of the plate are simply-supported boundary conditions along the y-axis, while the x-axis is located in the middle of the plate [5] as shown in Figure 9.

The bending moments and shear force in terms of the bending stiffness of orthotropic plate are:

(10a)${\mathrm{M}}_{\mathrm{xx}}=-{\mathrm{D}}_{11}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2}-{\mathrm{D}}_{12}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^2}$

(10b)${\mathrm{M}}_{\mathrm{yy}}=-{\mathrm{D}}_{22}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^2}-{\mathrm{D}}_{12}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2}$

(10c)${\mathrm{M}}_{\mathrm{xy}}=2{\mathrm{D}}_{66}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial \mathrm{x}\partial \mathrm{y}}$

(11a)${\mathrm{Q}}_{\mathrm{x}}=-\frac{\partial }{\partial \mathrm{x}}\left({\mathrm{D}}_{11}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2}+{\mathrm{D}}_{12}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^2}+2{\mathrm{D}}_{66}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^2}\right)$

(11b)${\mathrm{Q}}_{\mathrm{y}}=-\frac{\partial }{\partial \mathrm{y}}\left({\mathrm{D}}_{12}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2}+2{\mathrm{D}}_{66}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2}+{\mathrm{D}}_{22}\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^2}\right)$

As mentioned earlier in experimental section, the boundary conditions are simply supported along the y-axis and free supported along the x-axis. The application of the boundary condition and the solution of bending deflection and critical buckling load can be shown in Appendix A.

4. Special Case Using the Navier Solution

The special case assumes that the boundary conditions from all edges are simply supported. The analytic solution of the bending deflection through the plate thickness is completed using the Navier solution of classical laminate plate theory. Recall Equation (5):

(12)${\mathrm{D}}_{11}\ast \left(\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^4}\right)+2\ast \left({\mathrm{D}}_{12}+2{\ast \mathrm{D}}_{66}\right)\ast \left(\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2\partial {\mathrm{y}}^2}\right)+{\mathrm{D}}_{22}\ast \left(\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^4}\right)+{\mathsf{\rho }}_{\mathrm{o} }\ast \frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{t}}^2}+{\mathrm{N}}_{\mathrm{xx}}\ast \frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2}=0$

The solution of the above equation is below:

(13)${\mathrm{w}}_{\mathrm{o}}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }{\mathrm{w}}_{\mathrm{mn}}\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\beta }}_{\mathrm{n}}\ast \mathrm{y}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

where

$\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^4}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }{\mathsf{\alpha }}_{\mathrm{m}}^4{\ast \mathrm{w}}_{\mathrm{mn}}\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\beta }}_{\mathrm{n}}\ast \mathrm{y}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

$\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^4}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }{\mathsf{\beta }}_{\mathrm{n}}^4{\ast \mathrm{w}}_{\mathrm{mn}}\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\beta }}_{\mathrm{n}}\ast \mathrm{y}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

$\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2\partial {\mathrm{y}}^2}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }{\mathsf{\alpha }}_{\mathrm{m}}^2{\ast \mathsf{\beta }}_{\mathrm{n}}^2{\ast \mathrm{w}}_{\mathrm{mn}}\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\beta }}_{\mathrm{n}}\ast \mathrm{y}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

$\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{t}}^2}=-{\mathsf{\omega }}_{\mathrm{mn}}^2{\displaystyle \sum }_{\mathrm{m}=1}^{\infty }{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }{\mathrm{w}}_{\mathrm{mn}}\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\beta }}_{\mathrm{n}}\ast \mathrm{y}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

and

${\mathsf{\alpha }}_{\mathrm{m}}=\frac{\mathrm{m}\ast \mathsf{\pi }}{\mathrm{a}}$

${\mathsf{\beta }}_{\mathrm{n}}=\frac{\mathrm{n}\ast \mathsf{\pi }}{\mathrm{b}}$

Substitute the expressions ($\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^4}$, $\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^4}$ and $\frac{{\partial }^4{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2\partial {\mathrm{y}}^2}$) into Equation (12) and solve the plate equation for (${\mathrm{N}}_{\mathrm{xx}}$):

(14)${\mathrm{N}}_{\mathrm{xx}}=\frac{{\mathsf{\pi }}^2{\ast \mathrm{a}}^2}{{\mathrm{m}}^2}\left({\mathrm{D}}_{11}\ast \frac{{\mathrm{m}}^4}{{\mathrm{a}}^4}+2\ast \left({\mathrm{D}}_{12}+2{\ast \mathrm{D}}_{66}\right)\ast \frac{{\mathrm{m}}^2{\ast \mathrm{n}}^2}{{\mathrm{a}}^2{\ast \mathrm{b}}^2}+{\mathrm{D}}_{22}\ast \frac{{\mathrm{n}}^4}{{\mathrm{b}}^4}-\frac{{\mathsf{\rho }}_{\mathrm{o} }{\ast \mathsf{\omega }}_{\mathrm{mn}}^2}{{\mathsf{\pi }}^4}\right)$

By putting the buckling mode (n = 1) while Equation (14) is derived with respect to the buckling mode (m), to calculate the critical buckling load is as below:

(15)$\frac{{\mathrm{dN}}_{\mathrm{xx}}}{\mathrm{dm}}=0$

The critical buckling load is:

(16)${\mathrm{N}}_{\mathrm{xx}}=\frac{{\mathsf{\pi }}^2}{{\mathrm{b}}^2}\sqrt{{\mathrm{D}}_{11}{\mathrm{D}}_{22}}\left(\sqrt{\frac{{\mathrm{D}}_{22}-{\mathsf{\rho }}_{\mathrm{o} }{\ast \mathsf{\omega }}_{\mathrm{mn}}^2{\ast \mathrm{b}}^2}{{\mathrm{D}}_{22}}}+\frac{2\ast \left({\mathrm{D}}_{12}+2{\ast \mathrm{D}}_{66}\right)}{\sqrt{{\mathrm{D}}_{11}{\mathrm{D}}_{22}}}+\sqrt{\frac{{\mathrm{D}}_{22}}{{\mathrm{D}}_{22}-{\mathsf{\rho }}_{\mathrm{o} }{\ast \mathsf{\omega }}_{\mathrm{mn}}^2{\ast \mathrm{b}}^2}}-\frac{{\mathsf{\rho }}_{\mathrm{o} }{\ast \mathsf{\omega }}_{\mathrm{mn}}^2{\ast \mathrm{b}}^2}{{\mathsf{\pi }}^4}\ast \sqrt{\frac{1}{{\mathrm{D}}_{22}^2-{\mathrm{D}}_{22}{\ast \mathsf{\rho }}_{\mathrm{o} }{\ast \mathsf{\omega }}_{\mathrm{mn}}^2{\ast \mathrm{b}}^2}}\right)$

and

(17)$\frac{\mathrm{a}}{\mathrm{b}}=\mathrm{m}\ast {\left(\frac{{\mathrm{D}}_{11}}{{\mathrm{D}}_{22}-{\mathsf{\rho }}_{\mathrm{o} }{\ast \mathsf{\omega }}_{\mathrm{mn}}^2{\ast \mathrm{b}}^2}\right)}^{\frac{1}{4}}$

The buckling factor is

(18)${\mathrm{k}}_1=\frac{{\mathrm{m}}^2}{{\mathrm{c}}^2}\sqrt{\frac{{\mathrm{D}}_{11}}{{\mathrm{D}}_{22}}}+\frac{2\ast \left({\mathrm{D}}_{12}+2{\ast \mathrm{D}}_{66}\right)}{\sqrt{{\mathrm{D}}_{11}{\mathrm{D}}_{22}}}+\frac{{\mathrm{c}}^2}{{\mathrm{m}}^2}\sqrt{\frac{{\mathrm{D}}_{22}}{{\mathrm{D}}_{11}}}-\frac{{\mathsf{\rho }}_{\mathrm{o} }{\ast \mathsf{\omega }}_{\mathrm{mn}}^2{\ast \mathrm{c}}^2{\ast \mathrm{b}}^2}{{\mathrm{m}}^2{\ast \mathsf{\pi }}^4}\ast \frac{1}{\sqrt{{\mathrm{D}}_{11}{\mathrm{D}}_{22}}}$

5. Different Boundary Conditions

(1) The free–free boundary conditions along the x-axis are assumed to be clamped–clamped while the top and bottom edges of the compression test device stayed simply-supported (S-C-S-C).

B.C. (a) At y = 0, $\frac{\partial \mathrm{wo}}{\partial \mathrm{y}}=0$, ${\mathrm{Q}}_{\mathrm{y}}={\mathrm{V}}_{\mathrm{y}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }\frac{\tau }{\mathrm{a}\ast \mathrm{h}}\sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

B.C. (b) At y = b, $\frac{\partial \mathrm{wo}}{\partial \mathrm{y}}=0$, ${\mathrm{Q}}_{\mathrm{y}}={\mathrm{V}}_{\mathrm{y}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }\frac{\tau }{\mathrm{a}\ast \mathrm{h}}\sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

(2) The free–free boundary conditions along the x-axis are assumed to be clamped–simply supported while the top and bottom edges of the compression test device stayed simply supported (S-C-S-S).

B.C. (a) At y = 0, $\frac{\partial \mathrm{wo}}{\partial \mathrm{y}}=0$, ${\mathrm{Q}}_{\mathrm{y}}={\mathrm{V}}_{\mathrm{y}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }\frac{\tau }{\mathrm{a}\ast \mathrm{h}}\sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

B.C. (b) At y = b, ${\mathrm{w}}_{\mathrm{o}}=0$, ${\mathrm{Q}}_{\mathrm{y}}={\mathrm{V}}_{\mathrm{y}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }\frac{\tau }{\mathrm{a}\ast \mathrm{h}}\sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

(3) The free–free boundary conditions along x-axis are assumed to be free-clamped while the top and bottom edges of the compression test device stayed simply-supported (S-F-S-C).

B.C. (a) At y = 0, ${\mathrm{M}}_{\mathrm{yy}}=0$, ${\mathrm{Q}}_{\mathrm{y}}={\mathrm{V}}_{\mathrm{y}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }\frac{\tau }{\mathrm{a}\ast \mathrm{h}}\sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

B.C. (b) At y = b, $\frac{\partial \mathrm{wo}}{\partial \mathrm{y}}=0$, ${\mathrm{Q}}_{\mathrm{y}}={\mathrm{V}}_{\mathrm{y}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }\frac{\tau }{\mathrm{a}\ast \mathrm{h}}\sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

(4) The free–free boundary conditions along x-axis are assumed to be free-simply supported while the top and bottom edges of the compression test device stayed simply-supported (S-F-S-S).

B.C. (a) At y = 0, ${\mathrm{M}}_{\mathrm{yy}}=0$, ${\mathrm{Q}}_{\mathrm{y}}={\mathrm{V}}_{\mathrm{y}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }\frac{\tau }{\mathrm{a}\ast \mathrm{h}}\sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

B.C. (b) At y = b, ${\mathrm{w}}_{\mathrm{o}}=0$, ${\mathrm{Q}}_{\mathrm{y}}={\mathrm{V}}_{\mathrm{y}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }\frac{\tau }{\mathrm{a}\ast \mathrm{h}}\sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

The constants (A, B, C, D, E, F, G, H, I, J, K, L and M) and $({\mathrm{N}}_1, {\mathrm{N}}_2$, ${\mathrm{N}}_3$ and ${\mathrm{N}}_4$) are shown in Appendix B.

6. Numerical Simulation

In this paper, the finite element analysis is carried out by using ANSYS [16]. An element (SHELL 99) is selected for the main plate in which this element works with thin and thick plate-shell structures. The element (SHELL 99) has eight nodes with six degrees of freedom at each node; translations in the (x, y and z) axes and rotations about the (x, y and z) axes. This kind of elements used for layered applications for modeling composite plates and shells includes already the effect of transverse shear deformation. For the stiffened plate, the element (SOLSH 190) is used in the simulation, which is having eight nodes connectivity with three degrees of freedom at each node: translations in the nodal (x, y and z) directions. Another three degrees of freedom (rotation) has been added at each node in the nodal (x, y and z) directions. A first order shear deformation theory with the aid of Mindlin–Reissner of shell theory is used by ANSYS software. The global (x) coordinate is directed along the width of the plate, while the global (y) coordinate is directed along the length. The global (z) coordinate corresponds to the thickness direction, which is taken to be outward normal of the plate surface, as shown in Figure 8. For the contact analysis between the main plate and the stiffener, the element (TARGE 170) and the element (CONTA 174) have been used to create the mesh between the target and the contact surfaces. The stiffener represents the target, while the main plate is the source. The mesh generation for un-stiffened plates, stiffened plates with one stiffener and stiffened plates with two stiffeners are illustrated in Figure 10, Figure 11 and Figure 12. The convergence test is completed using a high number of elements to determine the exact size of the element when the value of bending deflection is settled down, as indicated in Figure 13 and Figure 14.

7. Largest Lyapunov Exponent Parameter

In the design, the largest Lyapunov exponent parameter is one of the indicators to detect the non-periodic motion and chaos of the bending deflection of the stiffened composite laminated plate. The higher the value of the Lyapunov exponent parameter, the more non-periodic the motion and chaos. The value of the Lyapunov exponent of the bending deflection being positive means that the motion of the bending deflection is non-periodic. The negative largest Lyapunov exponent value of the bending deflection indicates periodic motion. The Wolff algorithm code based on MATLABB software is used to extract the values of the largest Lyapunov exponent by monitoring the orbital divergence for the contact force [17]. Equations (19) and (20) were used to build the Wolff algorithm code of the dynamic tool [18].

(19)$\mathrm{d}\left(\mathrm{t}\right)={\mathrm{De}}^{\mathsf{\lambda }\mathrm{t}}$

(20)$\mathrm{y}\left(\mathrm{i}\right)=\frac{1}{\mathsf{\Delta }\mathrm{t}}\left[{\mathrm{lnd}}_{\mathrm{j}}\left(\mathrm{i}\right)\right]$

Figure 15 shows the numerical value of the largest Lyapunov exponent against the number of samples at a fiber volume fraction (80%) and aspect ratio (2.5) for an un-stiffened composite laminated plate. The value of the Lyapunov exponent parameter is taken at its respective equilibrium point. One column of the bending deflection with the values of time delay and embedding dimension are used in the Wolff algorithm to extract the numerical value of the largest Lyapunov exponent parameter. One column of the bending deflection has been taken from the ANSYS program, in which it used the dynamic tool of the Wolff algorithm code.

Time delay and embedding dimensions are numeric values in which it is needed in the Wolff algorithm code in the estimation of the local Lyapunov exponent parameter. The algorithm of the dynamic code of the global false nearest neighbors (GFNN) of the bending deflection is used to obtain the value of the embedding dimensions, in which it is estimated when the trend of global false nearest neighbors approaches zero [19]. In this paper, a local and optimal time delay was estimated in which it is recommended that the local time delay should be chosen to be dependent on the embedding dimensions [20]. The algorithm code of the average mutual information (AMI) is used to obtain the value of the time delay, and the first minimum time in the average mutual information trend represents the value of the time delay. The optimal time delay is calculated from the following equation:

(21)${\mathsf{\tau }}_{\mathrm{w}}=\left(\mathrm{P}-1\right){\mathsf{\tau }}^{\ast }$

where

${\mathsf{\tau }}^{\ast }$: is the local time delay.

${\mathsf{\tau }}_{\mathrm{w}}$: is the optimal for independence of time series.

P: is the embedding dimensions.

Figure 16 and Figure 17 show the comparison of time delay, while Figure 18 and Figure 19 show the comparison of embedding dimensions.

There is another way for estimating the largest Lyapunov exponent of the bending deflection that is called the average logarithmic divergence. The straight line represents the slope of the average logarithmic divergence, which reflects the value of the Lyapunov exponent parameter. The curve represents the logarithm function against time of the bending deflection. As known, the logarithm function treats any set of data by a straight line, which gives the value of the largest Lyapunov exponent directly. Figure 20 shows the comparison of the average logarithmic divergence of bending deflection against time at a fiber volume fraction (80%) and aspect ratio (2.5) for a plate with two stiffeners. The analytic largest Lyapunov exponent of the bending deflection is calculated after applying Equation (9) on one column of the bending deflection using the average logarithmic divergence approach, while the numerical simulation of the largest Lyapunov exponent of the bending deflection is determined using the ANSYS program. The experiment value of the largest Lyapunov exponent has been obtained after tracking the bending deflection in the direction of the plate thickness using the strain gauge through the strain meter.

8. Power Density Function of Fast Fourier Transfer (FFT)

Power density function is used to detect the non-periodic motion and chaos of the bending deflection of the stiffened composite laminated plate since it gives six frequencies peaks’ alongside the fundamental frequency [21]. The signal power is the most important factor for signal quality, in which the noise can be measured by using the indicator SNR ratio. SNR is the ratio between the wanted signal and the unwanted noise. When the SNR of the bending deflection signal has a maximum value, it means that there is no error or noise in the bending deflection signal. The SNR power of the bending deflection signal is measured by the dB scale since it could be either positive or negative. A negative SNR means that the signal power is lower than the noise power. When the power of the amplitude peak of the fundamental frequency and the other frequencies are decreased with the increasing of number of samples, it means that the motion of the bending deflection is non-periodic motion and chaos. When the frequencies’ peaks have not been started disappearing from (FFT) diagram, the motion of the bending deflection is quasi-periodic. Figure 21 and Figure 22 show the comparison of the power density function of the fast Fourier transform (FFT) of the bending deflection at the fiber volume fraction (25%) and aspect ratio (2.5) for the un-stiffened plate. The SNR of the bending deflection signal is positive and equal to (17.24 dB).

9. Results and Discussion

Table 1 shows the comparison of the critical buckling load of (S-F-S-F) boundary conditions for different orientations when the plate is stiffened with one and two stiffeners. The critical buckling load is increased with the increasing of orientations and number of stiffeners. The critical buckling load results at the orientations (0/90/0/90) and (0/90/90/0) and are the same when the main plate is stiffened by one stiffener. When the main plate is stiffened by two stiffeners, the results of critical buckling load are not the same at the orientations (0/90/0/90) and (0/90/90/0). The numerical results of the critical buckling load were done using the ANSYS program while the experiment results were carried out using strain gauge through a strain meter device. The analytic results of the critical buckling load were completed using the Levy solution of classical laminate plate theory.

Figure 23, Figure 24 and Figure 25 show the deflection contour distribution for the un-stiffened plate, stiffened plate with one stiffener and stiffened plate with two stiffeners of fiber orientation (0/90). The frequency was 12.31 Hz, 20.96 Hz and 23.94 Hz for the un-stiffened plate, stiffened plate with one stiffener and stiffened plate with two stiffeners, respectively. The central point in the middle of the main plate was considered in the calculation of the bending deflection contour. The maximum bending deflection decreased with the increasing of number of stiffeners. The ANSYS program was used in the calculation of the bending deflection contour.

Figure 26 shows the critical buckling load against the aspect ratio of (S-C-S-C and S-F-S-F) boundary conditions for plate with one stiffener in the presence of in-plane compression mechanical load and shear force. The critical buckling load increased with the increasing of the aspect ratio. The critical buckling load value for (S-F-S-F) boundary conditions is greater than the value of the critical buckling load for (S-C-S-C) boundary conditions. The two clamped boundary conditions in (S-C-S-C) generate shear stress in the opposite direction, which reduces the value of critical buckling load. The ANSYS program is used in the calculation of the critical buckling load.

Figure 27 shows the critical buckling load against the thickness ratio for (S-C-S-C and S-F-S-F) boundary conditions for plate with one stiffener in the presence of the in-plane compression mechanical load and shear force. The critical buckling load is varied sinusoidal with a thickness ratio in which it reaches a maximum value at thickness ratio (a/h = 20). The ANSYS program is used in the calculation of the critical buckling load.

Figure 28 shows the critical buckling load against the lamination angle for the plate with one stiffener for (S-C-S-C and S-F-S-F) different boundary conditions in the presence of the in-plane compression mechanical load and shear force. The critical buckling load is varied sinusoidal for (S-C-S-C) boundary conditions, while the critical buckling load is varied parabolic for (S-F-S-F) boundary conditions. The critical buckling load is the same at the lamination angle (40 and 60) degrees for (S-F-S-F) boundary conditions. The ANSYS program is used in the calculation of the critical buckling load.

Figure 29 shows the critical buckling load against a number of stiffeners at (S-F-S-F, S-S-S-S and S-C-S-C) boundary conditions in the presence of in-plane compression mechanical load and shear force. The critical buckling load is increased with the increasing of number of stiffeners, while the critical buckling load has a greater value at (S-C-S-C) boundary conditions. The ANSYS program is used in the calculation of the critical buckling load.

Figure 30 shows the critical buckling load against the stiffener–depth ratio for (S-S-S-S) boundary conditions in the presence of the in-plane compression mechanical load and shear force for the plate with one and two stiffeners. The buckling load is varied sinusoidal with the increasing stiffener–depth ratio. As mentioned, the critical buckling load is increased with the increasing of number of stiffeners. The ANSYS program is used in the calculation of the critical buckling load.

Figure 31 shows the buckling factor against the aspect ratio of (S-S-S-S) boundary conditions for the plate with one stiffener in the presence of in-plane compression mechanical load and shear force. The buckling factor of the buckling mode is decreased with the increasing aspect ratio until it settles down at buckling factor (20). The analytical results of the critical buckling load of the buckling mode have been completed using the Navier solution of classical laminate plate theory. Figure 32 shows the buckling factor against the aspect ratio for the plate with two stiffeners of (S-S-S-S) boundary conditions.

Table 2 shows the comparison of the critical buckling load when the composite laminated plate is supported with one and two stiffeners at different boundary conditions. The value of the critical buckling load is increased with the increasing of number of stiffeners for different boundary conditions. The numerical results of the critical buckling load were calculated using the ANSYS program, while the analytic results were calculated using the Levy solution of (S-C-S-C, S-C-S-S, S-F-S-C, S-F-S-S and S-F-S-F) boundary conditions. The analytic results of the (S-S-S-S) boundary condition were determined using the Navier solution of classical laminate plate theory.

Figure 33, Figure 34 and Figure 35 show the largest Lyapunov exponent against aspect ratios for the un-stiffened plate and plate with one and two stiffeners at different fiber volume fractions. As mentioned, that Lyapunov exponent is the indicator to the non-periodic motion and chaos (nonlinear dynamic phenomena) for the bending deflection of stiffened and un-stiffened composite laminated plate. The nonlinear dynamic behavior of the bending deflection is increased with the increasing of the aspect ratios and fiber volume fractions. All the values of the largest Lyapunov exponent are positive and above zero, which means that the motion of the bending deflection is non-periodic. The nonlinear dynamic phenomena are decreased with the use of beam stiffeners. The system with aspect ratio (2.5) and fiber volume fraction (${\mathsf{\upsilon }}_{\mathrm{f}}$ = 80%) for the un-stiffened plate is more chaotic than the other systems, as indicated in Figure 33. The nonlinear dynamic phenomena is very low for the system with aspect ratio (2.5) and fiber volume fraction (${\mathsf{\upsilon }}_{\mathrm{f}}$ = 80%) after using plate with two stiffeners, as shown in Figure 35. The ANSYS program is used in the calculation of bending deflection.

Figure 36 and Figure 37 show the comparison of the amplitude of fast Fourier transform (FFT) against frequency at aspect ratio (a/b = 1) and fiber volume fraction (${\mathsf{\upsilon }}_{\mathrm{f}}$ = 60%) for the plate with two stiffeners. The bending deflection is used in the calculation of the power density function of the fast Fourier transform (FFT). The power of the frequencies’ peaks is disappeared from (FFT) diagram, which gives indication to the chaotic system. The experiment results of the bending deflection were carried out using a strain gauge through a strain meter device, while the numerical results of the bending deflection were completed using the ANSYS program.

Figure 38 and Figure 39 show the experiment mapping of the amplitude of the fast Fourier transform (FFT) against frequency at different fiber volume fractions and different aspect ratios, respectively. The power density function of the fast Fourier transform (FFT) is decreased with the increasing of fiber volume fractions and aspect ratios. All the systems in Figure 38 and Figure 39 have non-periodic motion, since the frequencies’ peaks start disappearing from the (FFT) diagram for the unstiffened plate and plate with two stiffeners, respectively. The experiment results of the bending deflection were carried out using a strain gauge through a strain meter device.

Figure 40 shows the comparison of average logarithmic divergence against time at aspect ratio (a/b = 1.5) and fiber volume fraction (${\mathsf{\upsilon }}_{\mathrm{f}}$ = 60%). The nonlinear curve represents the average logarithmic divergence of the bending deflection motion for a plate with one stiffener in the presence of an in-plane compression mechanical load and shear force. The straight line represents the slope of average logarithmic divergence, which reflects the value of largest Lyapunov exponent. Moreover, Figure 41 shows the comparison of the average logarithmic divergence against the time at aspect ratio (a/b = 1) and fiber volume fraction (${\mathsf{\upsilon }}_{\mathrm{f}}$ = 60%) for an un-stiffened plate. The experiment results of the bending deflection that were used in the average logarithmic divergence were carried out using a strain gauge through a strain meter device. The numerical results of the bending deflection were found using the ANSYS program.

10. Conclusions

This article discussed how to suppress the nonlinear dynamic phenomena in a composite laminated plate using a beam stiffener in the directions of the applied load. The nonlinear dynamic behavior of the bending deflection was detected using the largest Lyapunov exponent parameter and power density function of the fast Fourier transform (FFT). The nonlinear dynamic behavior of the bending deflection increased with the increase of the aspect ratios and fiber volume fractions. All the values of the largest Lyapunov exponent are positive and above zero, which means that the motion of the bending deflection is non-periodic. The power density function of the fast Fourier transform (FFT) is decreased with the increasing of the fiber volume fractions and aspect ratios. All the systems have non-periodic motion, since the frequencies’ peaks started disappearing from the (FFT) diagram for the unstiffened plate and plate with two stiffeners, respectively. The system with aspect ratio (2.5) and fiber volume fraction (${\mathsf{\upsilon }}_{\mathrm{f}}$ = 80%) for the un-stiffened plate was more chaotic than the other systems. The nonlinear dynamic phenomena was very low for the system with aspect ratio (2.5) and fiber volume fraction (${\mathsf{\upsilon }}_{\mathrm{f}}$ = 80%) after using the plate with two stiffeners.

Funding

This research received no external funding.

Data Availability Statement

The author declares that the manuscript data will be available to the public upon request.

Acknowledgments

I want to first thank the reviewers and the editor for reviewing my paper. Additionally, I want to acknowledge the Department of Mechanical Engineering at the University of Baghdad for their total support for finishing this research.

Conflicts of Interest

The author declares that he has no conflict of interest.

Nomenclature

a	Length of large span of rectangular plate, mm
b	Length of small span of rectangular plate, mm
h	Thickness of the composite plate, mm
${\mathsf{\sigma }}_{\mathrm{xx}}$, ${\mathsf{\sigma }}_{\mathrm{yy}}, {\mathsf{\sigma }}_{\mathrm{xy}}$	Lamina stresses field vector, N/mm2
${\mathsf{\epsilon }}_{\mathrm{xx}}^{\left(0\right)}, {\mathsf{\epsilon }}_{\mathrm{yy}}^{\left(0\right)}$, ${\mathsf{\gamma }}_{\mathrm{xy}}^{\left(0\right)}$	Extensional strain vector
${\mathsf{\epsilon }}_{\mathrm{xx}}^{\left(1\right)}, {\mathsf{\epsilon }}_{\mathrm{yy}}^{\left(1\right)}$, ${\mathsf{\gamma }}_{\mathrm{xy}}^{\left(1\right)}$	Flexural (bending) strain vector
${\mathrm{u}}_{\mathrm{o}}$,${ \mathrm{v}}_{\mathrm{o}}$	Displacement components in the 2-D coordinate system, mm
${\mathrm{Q}}_{\mathrm{ij}}$	Reduced stiffness elements, N/mm2
${\mathrm{w}}_{\mathrm{o}}$	Mid plane deflection along z-direction, mm
${\mathsf{\upsilon }}_{\mathrm{f}}$	Volume fraction for the fiber
$\mathsf{\phi }$	Fiber weight fraction
${\mathrm{m}}_{\mathrm{f}}$	Mass density for the fiber, kg
${\mathrm{m}}_{\mathrm{C}}$	Mass density for the composite, kg
${\mathrm{E}}_{11}$, ${\mathrm{E}}_{22}$	Modulus of elasticity along (x and y) directions, respectively, N/mm2
${\mathsf{\nu }}_{12}$, ${\mathsf{\nu }}_{21}$	Poisson’s ratio in plane 1-2 and the perpendicular plane 2-1, respectively
${\mathrm{t}}_{\mathrm{x}}, {\mathrm{h}}_{\mathrm{x}}, {\mathrm{b}}_{\mathrm{x}}$	Thickness, depth, distance between stiffeners when the stiffener placed along x-axis, mm
${\mathrm{t}}_{\mathrm{y}}, {\mathrm{h}}_{\mathrm{y}}, {\mathrm{b}}_{\mathrm{y}}$	Thickness, depth, distance between stiffeners when the stiffener placed along y-axis, mm
${\mathrm{N}}_{\begin{array}{c}cr\\ \end{array}}$	Critical buckling load, N
${\mathsf{\rho }}_{\mathrm{o} }$	Mass per unit area of the composite plate, kg/m2
m, n	Double trigonometric of Fourier series
${\mathsf{\omega }}_{\mathrm{mn}}$	Natural frequency of the composite laminated plate
${\mathrm{k}}_1$	Buckling factor
GFNN	Global false nearest neighbors for extracting embedding dimensions
AMI	Average mutual information algorithm for quantifying time delay
D	Distance between trajectories
d(t)	Percentage of variation in the distance between trajectories
${\mathrm{d}}_{\mathrm{j}}\left(\mathrm{i}\right)$	Distance between the pair (j-th) at nearest neighbors (i) of the trajectories
t	Time strides
y(i)	Data of follower displacement after curve fitting
Δt	Various time strides
λ	Largest Lyapunov exponent parameter
S-S-S-S	Simply-simply-simply-simply supported boundary conditions
S-F-S-F	Simply-free-simply-free supported boundary conditions
S-C-S-C	Simply-clamped-simply-clamped supported boundary conditions
S-F-S-S	Simply-free-simply-simply supported boundary conditions
S-F-S-C	Simply-free-simply-clamped supported boundary conditions
S-C-S-S	Simply-clamped-simply-simply supported boundary conditions

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Manufacturing of composite laminated plate.

Figure 2. Buckling test.

Figure 3. Buckling test of composite plate without stiffener.

Figure 4. Buckling test of composite plate with one stiffener.

Figure 5. Buckling test of composite plate with two stiffeners.

Figure 6. Bending deflection against in-plane compression mechanical load for plate with one and two stiffeners.

Figure 7. Bending deflection against in-plane compression mechanical load for un-stiffened plate and plate with one and two stiffeners.

Figure 8. Distribution of beam stiffeners along (x and y) axis.

Figure 9. Loading and the boundary condition of Levy solution.

Figure 10. Mesh contour of composite laminated plate without stiffener.

Figure 11. Mesh contour for composite laminated plate with one stiffener.

Figure 12. Mesh contour for composite laminated plate with two stiffeners.

Figure 13. Convergence test of bending deflection against number of elements for plate with one and two stiffeners.

Figure 14. Convergence test of bending deflection against numbers of elements for plate with and without stiffeners.

Figure 15. Local Lyapunov exponent against number of samples for ([Forumla omitted. See PDF.] = 80%) and (a/b = 2.5).

Figure 16. Time delay of numerical simulation.

Figure 17. Time delay of experiment test.

Figure 18. Embedding dimensions of numerical simulation.

Figure 19. Embedding dimensions of experiment test.

Figure 20. Average logarithmic divergence against the time for ([Forumla omitted. See PDF.] = 80%) and (a/b = 2.5).

Figure 21. Power spectrum analysis of the experiment test of the bending deflection for ([Forumla omitted. See PDF.] = 25%) and (a/b = 2.5).

Figure 22. Power spectrum analysis of the numerical simulation of the bending deflection for ([Forumla omitted. See PDF.] = 25%) and (a/b = 2.5).

Figure 23. Bending deflection contour for un-stiffened plate.

Figure 24. Bending deflection contour for plate with one stiffener.

Figure 25. Bending deflection contour for plate with two stiffeners.

Figure 26. Critical buckling load against aspect ratio for (S-C-S-C) and (S-F-S-F) boundary conditions.

Figure 27. Critical buckling load against thickness ratio for (S-C-S-C) and (S-F-S-F) boundary conditions.

Figure 28. Critical buckling load against lamination angle for (S-C-S-C) and (S-F-S-F) boundary conditions.

Figure 29. Critical buckling load against number of stiffeners at (S-F-S-F, S-S-S-S and S-C-S-C) boundary conditions.

Figure 30. Critical buckling load against stiffener–depth ratio for (S-S-S-S) boundary conditions.

Figure 31. Buckling factor against aspect ratio for (S-S-S-S) simply supported boundary conditions for plate with one stiffener.

Figure 32. Buckling factor against aspect ratio for (S-S-S-S) simply supported boundary conditions for plate with two stiffeners.

Figure 33. Largest Lyapunov exponent against aspect ratio for un-stiffened plate at different fiber volume fractions.

Figure 34. Largest Lyapunov exponent against aspect ratio for two stiffeners at different fiber volume fractions.

Figure 35. Largest Lyapunov exponent against aspect ratio for one stiffener at different fiber volume fractions.

Figure 36. Power spectrum analysis of experiment test.

Figure 37. Power spectrum analysis of numerical simulation.

Figure 38. Mapping of power density function at different fiber volume fractions for un-stiffened plate.

Figure 39. Mapping of power density function at different aspect ratios for plate with two stiffeners.

Figure 40. Comparison of average logarithmic divergence against time for stiffened plate with one stiffener.

Figure 41. Comparison of average logarithmic divergence against time for un-stiffened plate.

Table 1

Comparison of critical buckling load for (S-F-S-F) boundary conditions against fiber orientation for one and two stiffeners.

Stiffeners	Orientation	Experiment	Numerical	Analytic
One Stiffener	0/90	12.07	11.645	11.923
	0/90/0	13.274	12.249	13.03
	0/90/0/90	29.94	28.944	29.53
	0/90/90/0	30.387	29.26	29.856
Two Stiffeners	0/90	14.016	13.334	13.738
	0/90/0	15.446	14.665	15.127
	0/90/0/90	34.046	33.337	33.691
	0/90/90/0	44.366	43.603	43.752

Table 2

Comparison critical buckling load at different boundary conditions for one and two stiffeners.

The Support	One Stiffener		Two Stiffeners
	Numerical	Analytic	Numerical	Analytic
S-C-S-C	247.5	245.23	264.25	260.18
S-C-S-S	245.1	242.52	279.74	276.24
S-S-S-S	220.6	218.24	228.33	225.79
S-F-S-C	100.8	97.1	131.25	126.34
S-F-S-S	118.2	114.97	133.55	130.25
S-F-S-F	7.43	7.25	14.544	13.176

Appendix A

B.C. (1) At y = 0, ${\mathrm{M}}_{\mathrm{yy}}=0$, ${\mathrm{Q}}_{\mathrm{y}}={\mathrm{V}}_{\mathrm{y}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }\frac{\tau }{\mathrm{a}\ast \mathrm{h}}\sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

B.C. (2) At y = b, ${\mathrm{M}}_{\mathrm{yy}}=0$, ${\mathrm{Q}}_{\mathrm{y}}={\mathrm{V}}_{\mathrm{y}}={\displaystyle \sum }_{\mathrm{m}=1}^{\infty }\frac{\tau }{\mathrm{a}\ast \mathrm{h}}\sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)$

The expressions ($\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2}$, $\frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^2}$, $\frac{{\partial }^3{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2\partial \mathrm{y}}$, $\frac{{\partial }^3{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^3}$) can be prepared before I substitute for the boundary conditions:$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2}=-{\mathsf{\alpha }}_{\mathrm{m}}^2[{{\displaystyle \sum }}_{\mathrm{m}=1}^{\infty }[{\mathrm{N}}_1\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{y}\right)+{\mathrm{N}}_2\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{y}\right)+{\mathrm{N}}_3\ast \cos \left(\mathsf{\beta }\ast \mathrm{y}\right)+{\mathrm{N}}_4\ast \sin \left(\mathsf{\beta }\ast \mathrm{y}\right)]\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)]}\\ {\displaystyle \frac{\partial {\mathrm{w}}_{\mathrm{o}}}{\partial \mathrm{y}}={{\displaystyle \sum }}_{\mathrm{m}=1}^{\infty }[{\mathrm{N}}_1\ast \mathsf{\alpha }\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{y}\right)+{\mathrm{N}}_2\ast \mathsf{\alpha }\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{y}\right)-{\mathrm{N}}_3\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{y}\right)+{\mathrm{N}}_4\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{y}\right)]\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)}\\ {\displaystyle \frac{{\partial }^2{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^2}={{\displaystyle \sum }}_{\mathrm{m}=1}^{\infty }[{\mathrm{N}}_1{\ast \mathsf{\alpha }}^2\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{y}\right)+{\mathrm{N}}_2{\ast \mathsf{\alpha }}^2\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{y}\right)-{\mathrm{N}}_3{\ast \mathsf{\beta }}^2\ast \cos \left(\mathsf{\beta }\ast \mathrm{y}\right)-{\mathrm{N}}_4{\ast \mathsf{\beta }}^2\ast \sin \left(\mathsf{\beta }\ast \mathrm{y}\right)]\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)}\\ {\displaystyle {\mathrm{Q}}_{\mathrm{y}}=-\left({\mathrm{D}}_{12}\frac{{\partial }^3{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2\partial \mathrm{y}}+2{\mathrm{D}}_{66}\frac{{\partial }^3{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2\partial \mathrm{y}}+{\mathrm{D}}_{22}\frac{{\partial }^3{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^3}\right)}\\ {\displaystyle \frac{{\partial }^3{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{x}}^2\partial \mathrm{y}}=-{\mathsf{\alpha }}_{\mathrm{m}}^2[{{\displaystyle \sum }}_{\mathrm{m}=1}^{\infty }[{\mathrm{N}}_1\ast \mathsf{\alpha }\ast \sin \mathrm{h}\left(\mathsf{\alpha }\ast \mathrm{y}\right)+{\mathrm{N}}_2\ast \mathsf{\alpha }\ast \cos \mathrm{h}\left(\mathsf{\alpha }\ast \mathrm{y}\right)-{\mathrm{N}}_3\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{y}\right)+{\mathrm{N}}_4\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{y}\right)]\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)]}\\ {\displaystyle \frac{{\partial }^3{\mathrm{w}}_{\mathrm{o}}}{\partial {\mathrm{y}}^3}={{\displaystyle \sum }}_{\mathrm{m}=1}^{\infty }[{\mathrm{N}}_1{\ast \mathsf{\alpha }}^3\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{y}\right)+{\mathrm{N}}_2{\ast \mathsf{\alpha }}^3\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{y}\right)+{\mathrm{N}}_3{\ast \mathsf{\beta }}^3\ast \sin \left(\mathsf{\beta }\ast \mathrm{y}\right)-{\mathrm{N}}_4{\ast \mathsf{\beta }}^3\ast \cos \left(\mathsf{\beta }\ast \mathrm{y}\right)]\ast \sin \left({\mathsf{\alpha }}_{\mathrm{m}}\ast \mathrm{x}\right)\ast \sin \left({\mathsf{\omega }}_{\mathrm{mn}}\ast \mathrm{t}\right)}\end{array}$$

Apply the boundary conditions listed above on Equation (9) and after simplification to obtain (A1)$${\mathrm{N}}_1\ast \left(-{\mathrm{D}}_{22}\ast {\mathsf{\alpha }}^2-{\mathrm{D}}_{12}\right)+{\mathrm{N}}_3\ast \left({\mathrm{D}}_{22}\ast {\mathsf{\beta }}^2-{\mathrm{D}}_{12}\right)=0$$ (A2)$${\mathrm{N}}_2\ast \left({\mathrm{D}}_{12}\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }+2{\mathrm{D}}_{66}\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }-{\mathrm{D}}_{22}\ast {\mathsf{\alpha }}^3\right)+{\mathrm{N}}_4\ast \left({\mathrm{D}}_{12}\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+2{\mathrm{D}}_{66}{\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+{\mathrm{D}}_{22}\ast {\mathsf{\beta }}^3\right)=\frac{\tau }{\mathrm{a}\ast \mathrm{h}}$$ (A3)$$\begin{array}{l}{\displaystyle {\mathrm{N}}_1\ast \left(-{\mathrm{D}}_{22}\ast {\alpha }^2\ast \cosh \left(\alpha \ast \mathrm{b}\right)+{\mathrm{D}}_{12}\ast {\alpha }_{\mathrm{m}}^2\ast \cosh \left(\alpha \ast \mathrm{b}\right)\right)+{\mathrm{N}}_2\ast \left(-{\mathrm{D}}_{22}\ast {\alpha }^2\ast \sinh \left(\alpha \ast \mathrm{b}\right)+{\mathrm{D}}_{12}\ast {\alpha }_{\mathrm{m}}^2\ast \sinh \left(\alpha \ast \mathrm{b}\right)\right)+{\mathrm{N}}_3}\\ {\displaystyle \ast \left({\mathrm{D}}_{22}\ast {\mathsf{\beta }}^2\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)+{\mathrm{D}}_{12}\ast {\alpha }_{\mathrm{m}}^2\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)+{\mathrm{N}}_4\ast \left({\mathrm{D}}_{22}\ast {\mathsf{\beta }}^2\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)+{\mathrm{D}}_{12}\ast {\alpha }_{\mathrm{m}}^2\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)=0}\end{array}$$ (A4)$$\begin{array}{l}{\displaystyle {\mathrm{N}}_1\ast \left({\mathrm{D}}_{12}\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}\ast {\mathsf{\alpha }}^3\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)+{\mathrm{N}}_2}\\ {\displaystyle \ast \left({\mathrm{D}}_{12}\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}\ast {\mathsf{\alpha }}^3\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)+{\mathrm{N}}_3}\\ {\displaystyle \ast \left(-{\mathrm{D}}_{12}\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)-2{\mathrm{D}}_{66}\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}\ast {\mathsf{\beta }}^3\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)+{\mathrm{N}}_4}\\ {\displaystyle \ast \left({\mathrm{D}}_{12}\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)-{\mathrm{D}}_{22}\ast {\mathsf{\beta }}^3\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)=\frac{\mathsf{\tau }}{\mathrm{a}\ast \mathrm{h}}}\end{array}$$

After simplification, Equations (A1)–(A4) will be (A5)$${\mathrm{N}}_1\ast \mathrm{A}+{\mathrm{N}}_3\ast \mathrm{B}=0$$ (A6)$${\mathrm{N}}_2\ast \mathrm{C}+{\mathrm{N}}_4\ast \mathrm{D}=\mathrm{M}$$ (A7)$${\mathrm{N}}_1\ast \mathrm{E}+{\mathrm{N}}_2\ast \mathrm{F}+{\mathrm{N}}_3\ast \mathrm{G}+{\mathrm{N}}_4\ast \mathrm{H}=0$$ (A8)$${\mathrm{N}}_1\ast \mathrm{I}+{\mathrm{N}}_2\ast \mathrm{J}+{\mathrm{N}}_3\ast \mathrm{K}+{\mathrm{N}}_4\ast \mathrm{L}=\mathrm{M}$$ where $$\begin{array}{l}{\displaystyle \mathrm{A}=\left(-{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^2-{\mathrm{D}}_{12}\right)}\\ {\displaystyle \mathrm{B}=\left(-{\mathrm{D}}_{22}\ast {\mathsf{\alpha }}^2-{\mathrm{D}}_{12}\right)}\\ {\displaystyle \mathrm{C}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }-{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\right)}\\ {\displaystyle \mathrm{D}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+2{\mathrm{D}}_{66}{\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\right)}\\ {\displaystyle \mathrm{E}=\left(-{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^2\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{F}=\left(-{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^2\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{G}=\left({\mathrm{D}}_{22}{\ast \mathsf{\beta }}^2\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)+{\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{H}=\left({\mathrm{D}}_{22}{\ast \mathsf{\beta }}^2\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)+{\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{I}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{J}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{K}=\left(-{\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)-2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{L}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)-{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{M}=\frac{\tau }{\mathrm{a}\ast \mathrm{h}}}\end{array}$$

And;

A, B, C, D, E, F, G, H, I, J, K, L and M are constants.

Maple software is used to solve for the constants $({\mathrm{N}}_1, {\mathrm{N}}_2$, ${\mathrm{N}}_3$, ${\mathrm{N}}_4$):(A9)$$\begin{array}{l}{\displaystyle {\mathrm{J}}_1=\frac{\mathrm{B}\ast \left(\mathrm{F}\ast \mathrm{L} - \mathrm{H}\ast \mathrm{J}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J}-\mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_2=\frac{\mathrm{B}\ast \left(\mathrm{C}\ast \mathrm{H} - \mathrm{D}\ast \mathrm{F}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_1={\mathrm{J}}_1+{\mathrm{J}}_2}\end{array}$$ (A10)$$\begin{array}{l}{\displaystyle {\mathrm{J}}_3=\frac{\left(\mathrm{A}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{H}\ast \mathrm{K} - \mathrm{B}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{H}\ast \mathrm{I}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_4=\frac{\mathrm{D}\ast \left(\mathrm{A}\ast \mathrm{G} - \mathrm{B}\ast \mathrm{E}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_2={\mathrm{J}}_3-{\mathrm{J}}_4}\end{array}$$ (A11)$$\begin{array}{l}{\displaystyle {\mathrm{J}}_5=\frac{\mathrm{A}\ast \left(\mathrm{F}\ast \mathrm{L} - \mathrm{H}\ast \mathrm{J}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_6=\frac{\left(\mathrm{C}\ast \mathrm{H} - \mathrm{D}\ast \mathrm{F}\right)\ast \mathrm{A}\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_3=-{\mathrm{J}}_5-{\mathrm{J}}_6}\end{array}$$ (A12)$$\begin{array}{l}{\displaystyle {\mathrm{J}}_7=\frac{\left(\mathrm{A}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{G}\ast \mathrm{J} + \mathrm{B}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{F}\ast \mathrm{I}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_8=\frac{\left(\mathrm{A}\ast \mathrm{G} - \mathrm{B}\ast \mathrm{E}\right)\ast \mathrm{C}\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_4={\mathrm{J}}_7+{\mathrm{J}}_8}\end{array}$$

Substitute the constants $({\mathrm{N}}_1, {\mathrm{N}}_2$, ${\mathrm{N}}_3$, ${\mathrm{N}}_4$) into Equation (9) to obtain the general equation of the plate. Solve for the plate equation to determine the in-plane compression mechanical load and derive the output equation with respect to the buckling mode (m) along the x-axis to calculate the critical buckling load. From Equation (A1), the critical buckling load can be calculated as below in the following equation:(A13)$${\mathrm{N}}_{\mathrm{xx}}=\frac{{\left({\mathrm{N}}_1{\ast \mathrm{D}}_{12}+{\mathrm{N}}_3{\ast \mathrm{D}}_{12}+{\mathrm{N}}_1{\ast \mathrm{D}}_{22}{\ast \mathrm{O}}_1{\ast \mathsf{\alpha }}_{\mathrm{m}}^2-{\mathrm{N}}_3{\ast \mathrm{D}}_{22}{\ast \mathrm{O}}_1{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\right)}^2}{{\left({\mathrm{N}}_1+{\mathrm{N}}_3\right)}^2\ast {\left({\mathrm{D}}_{22}\right)}^2\ast {\mathsf{\alpha }}_{\mathrm{m}}^2\ast {\mathrm{O}}_4}-\frac{{\mathrm{O}}_1^2\ast {\mathsf{\alpha }}_{\mathrm{m}}^2}{{\mathrm{O}}_4}+\frac{{\mathrm{O}}_2\ast {\mathsf{\alpha }}_{\mathrm{m}}^2}{{\mathrm{O}}_4}-\frac{{\mathrm{O}}_3}{{\mathsf{\alpha }}_{\mathrm{m}}^2\ast {\mathrm{O}}_4}$$ where $${\mathrm{O}}_1=\frac{\left({\mathrm{D}}_{12}+2{\ast \mathrm{D}}_{66}\right)}{{\mathrm{D}}_{22}}, {\mathrm{O}}_2=\frac{{\mathrm{D}}_{11}}{{\mathrm{D}}_{22}}, {\mathrm{O}}_3=\frac{{\mathsf{\rho }}_{\mathrm{o} }{\ast \mathsf{\omega }}_{\mathrm{mn}}^2}{{\mathrm{D}}_{22}}, {\mathrm{O}}_4=\frac{1}{{\mathrm{D}}_{22}}$$

Derive Equation (A13) with respect to (${\mathsf{\alpha }}_{\mathrm{m}}$):(A14)$$\frac{\mathrm{d}\left({\mathrm{N}}_{\mathrm{xx}}\right)}{\mathrm{d}({\mathsf{\alpha }}_{\mathrm{m}})}=\left(\left(2\ast {\mathrm{O}}_2\ast {\mathsf{\alpha }}_{\mathrm{m}}^4+2\ast {\mathrm{O}}_3\right)\ast {\mathrm{N}}_1^2-8\ast \left(\left({\mathrm{O}}_1^2-\frac{{\mathrm{O}}_2}{2}\right)\ast {\mathsf{\alpha }}_{\mathrm{m}}^4-\frac{{\mathrm{O}}_3}{2}\right)\ast {\mathrm{N}}_1\ast {\mathrm{N}}_3+2\ast {\mathrm{N}}_3^2\ast \left({\mathrm{O}}_2\ast {\mathsf{\alpha }}_{\mathrm{m}}^4+{\mathrm{O}}_3\right)\right)\ast {({\mathrm{D}}_{22})}^2-2\ast {({\mathrm{D}}_{12})}^2\ast {\left({\mathrm{N}}_1+{\mathrm{N}}_3\right)}^2=0$$

Solve Equation (A14) to find the buckling mode:(A15)$$\mathrm{m}=\frac{\mathrm{a}}{\mathsf{\pi }}\ast {\left(\frac{{\left({\mathrm{N}}_1+{\mathrm{N}}_3\right)}^2\ast \left(-{\mathrm{D}}_{22}{}^2\ast {\mathrm{O}}_3\right)+{\mathrm{D}}_{12}{}^2}{{\mathrm{D}}_{22}{}^2\ast \left(-4\ast {\mathrm{N}}_1\ast {\mathrm{N}}_3\ast {\mathrm{O}}_1^2+{\mathrm{N}}_1^2\ast {\mathrm{O}}_2+2\ast {\mathrm{N}}_1\ast {\mathrm{N}}_3\ast {\mathrm{O}}_2+{\mathrm{N}}_3^2\ast {\mathrm{O}}_2\right)}\right)}^{1/4}$$

Substitute Equation (A15) into Equation (A13) to determine the critical buckling load:(A16)$${\mathrm{N}}_{\mathrm{xx}}\left(\mathrm{critical}\right)=\frac{{\left({\mathrm{N}}_1\ast {\mathrm{D}}_{12}+{\mathrm{N}}_3\ast {\mathrm{D}}_{12}+{\mathrm{N}}_1\ast {\mathrm{D}}_{22}\ast {\mathrm{O}}_1\ast {\mathrm{O}}_5-{\mathrm{N}}_3\ast {\mathrm{D}}_{22}\ast {\mathrm{O}}_1\ast {\mathrm{O}}_5\right)}^2}{{\left({\mathrm{N}}_1+{\mathrm{N}}_3\right)}^2\ast {\mathrm{D}}_{22}{}^2\ast {\mathrm{O}}_4\ast {\mathrm{O}}_5-\frac{{\mathrm{O}}_1^2\ast {\mathrm{O}}_5}{{\mathrm{O}}_4}+\frac{{\mathrm{O}}_2\ast {\mathrm{O}}_5}{{\mathrm{O}}_4}-\frac{{\mathrm{O}}_3}{{\mathrm{O}}_4\ast {\mathrm{O}}_5}}$$

Appendix B

The constants (A, B, C, D, E, F, G, H, I, J, K, L and M) and $({\mathrm{N}}_1, {\mathrm{N}}_2$, ${\mathrm{N}}_3$, ${\mathrm{N}}_4$) using different boundary conditions are as below:

• (1)

  (S-C-S-C) Boundary Conditions:

$$\begin{array}{l}{\displaystyle \mathrm{A}=\mathsf{\alpha }}\\ {\displaystyle \mathrm{B}=\mathsf{\beta }}\\ {\displaystyle \mathrm{C}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }-{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\right)}\\ {\displaystyle \mathrm{D}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+2{\mathrm{D}}_{66}{\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\right)}\\ {\displaystyle \mathrm{E}= \mathsf{\alpha }\ast \sin \mathrm{h}\left(\mathsf{\alpha }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{F}= \mathsf{\alpha }\ast \cos \mathrm{h}\left(\mathsf{\alpha }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{G}=-\mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{H}= \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{I}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{J}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{K}=\left(-{\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)-2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{L}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)-{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{M}=\frac{\tau }{\mathrm{a}\ast \mathrm{h}}}\\ {\displaystyle {\mathrm{N}}_1=\frac{\left(\mathrm{A}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{B}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{B}\ast \mathrm{G}\ast \mathrm{J}\right)\ast \mathrm{M}}{\left(\mathrm{E}\ast \mathrm{K} - \mathrm{G}\ast \mathrm{I}\right)\ast \left(\mathrm{A}\ast \mathrm{D} - \mathrm{B}\ast \mathrm{C}\right)}-\frac{\mathrm{G}\ast \mathrm{M}}{\mathrm{E}\ast \mathrm{K} - \mathrm{G}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_2=-\frac{\mathrm{B}\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{D} - \mathrm{B}\ast \mathrm{C}}}\\ {\displaystyle {\mathrm{N}}_3=-\frac{\left(\mathrm{A}\ast \mathrm{E}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{H}\ast \mathrm{I} - \mathrm{B}\ast \mathrm{E}\ast \mathrm{J} + \mathrm{B}\ast \mathrm{F}\ast \mathrm{I}\right)\ast \mathrm{M}}{\left(\mathrm{E}\ast \mathrm{K} - \mathrm{G}\ast \mathrm{I}\right)\ast \left(\mathrm{A}\ast \mathrm{D} - \mathrm{B}\ast \mathrm{C}\right)}+\frac{\mathrm{E}\ast \mathrm{M}}{\mathrm{E}\ast \mathrm{K} - \mathrm{G}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_4=\frac{\mathrm{A}\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{D} - \mathrm{B}\ast \mathrm{C}}}\end{array}$$

• (2)

  (S-C-S-S) Boundary Conditions:

$$\begin{array}{l}{\displaystyle \mathrm{A}=\mathsf{\alpha }}\\ {\displaystyle \mathrm{B}=\mathsf{\beta }}\\ {\displaystyle \mathrm{C}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }-{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\right)}\\ {\displaystyle \mathrm{D}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+2{\mathrm{D}}_{66}{\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\right)}\\ {\displaystyle \mathrm{E}= \cos \mathrm{h}\left(\mathsf{\alpha }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{F}= \sin \mathrm{h}\left(\mathsf{\alpha }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{G}= \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{H}= \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)}\\ {\displaystyle {\mathrm{N}}_1=\frac{\left(\mathrm{A}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{B}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{B}\ast \mathrm{G}\ast \mathrm{J}\right)\ast \mathrm{M}}{\left(\mathrm{E}\ast \mathrm{K} - \mathrm{G}\ast \mathrm{I}\right)\ast \left(\mathrm{A}\ast \mathrm{D} - \mathrm{B}\ast \mathrm{C}\right)}-\frac{\mathrm{G}\ast \mathrm{M}}{\mathrm{E}\ast \mathrm{K} - \mathrm{G}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_2=-\frac{\mathrm{B}\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{D} - \mathrm{B}\ast \mathrm{C}}}\\ {\displaystyle {\mathrm{N}}_3=-\frac{\left(\mathrm{A}\ast \mathrm{E}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{H}\ast \mathrm{I} - \mathrm{B}\ast \mathrm{E}\ast \mathrm{J} + \mathrm{B}\ast \mathrm{F}\ast \mathrm{I}\right)\ast \mathrm{M}}{\left(\mathrm{E}\ast \mathrm{K} - \mathrm{G}\ast \mathrm{I}\right)\ast \left(\mathrm{A}\ast \mathrm{D} - \mathrm{B}\ast \mathrm{C}\right)}+\frac{\mathrm{E}\ast \mathrm{M}}{\mathrm{E}\ast \mathrm{K} - \mathrm{G}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_4=\frac{\mathrm{A}\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{D} - \mathrm{B}\ast \mathrm{C}}}\end{array}$$

• (3)

  (S-F-S-C) Boundary Conditions:

$$\begin{array}{l}{\displaystyle \mathrm{A}=\left(-{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^2-{\mathrm{D}}_{12}\right)}\\ {\displaystyle \mathrm{B}=\left(-{\mathrm{D}}_{22}\ast {\mathsf{\alpha }}^2-{\mathrm{D}}_{12}\right)}\\ {\displaystyle \mathrm{C}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }-{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\right)}\\ {\displaystyle \mathrm{D}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+2{\mathrm{D}}_{66}{\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\right)}\\ {\displaystyle \mathrm{E}= \mathsf{\alpha }\ast \sin \mathrm{h}\left(\mathsf{\alpha }\ast \mathrm{b}\right). }\\ {\displaystyle \mathrm{F}= \mathsf{\alpha }\ast \cos \mathrm{h}\left(\mathsf{\alpha }\ast \mathrm{b}\right). }\\ {\displaystyle \mathrm{G}=-\mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right). }\\ {\displaystyle \mathrm{H}= \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{I}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{J}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{K}=\left(-{\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)-2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{L}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)-{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{M}=\frac{\tau }{\mathrm{a}\ast \mathrm{h}}}\\ {\displaystyle {\mathrm{J}}_1=\frac{\mathrm{B}\ast \left(\mathrm{F}\ast \mathrm{L} - \mathrm{H}\ast \mathrm{J}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J}-\mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_2=\frac{\mathrm{B}\ast \left(\mathrm{C}\ast \mathrm{H} - \mathrm{D}\ast \mathrm{F}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_1={\mathrm{J}}_1+{\mathrm{J}}_2}\\ {\displaystyle {\mathrm{J}}_3=\frac{\left(\mathrm{A}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{H}\ast \mathrm{K} - \mathrm{B}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{H}\ast \mathrm{I}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_4=\frac{\mathrm{D}\ast \left(\mathrm{A}\ast \mathrm{G} - \mathrm{B}\ast \mathrm{E}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_2={\mathrm{J}}_3-{\mathrm{J}}_4}\\ {\displaystyle {\mathrm{J}}_5=\frac{\mathrm{A}\ast \left(\mathrm{F}\ast \mathrm{L} - \mathrm{H}\ast \mathrm{J}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_6=\frac{\left(\mathrm{C}\ast \mathrm{H} - \mathrm{D}\ast \mathrm{F}\right)\ast \mathrm{A}\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_3=-{\mathrm{J}}_5-{\mathrm{J}}_6}\\ {\displaystyle {\mathrm{J}}_7=\frac{\left(\mathrm{A}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{G}\ast \mathrm{J} + \mathrm{B}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{F}\ast \mathrm{I}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_8=\frac{\left(\mathrm{A}\ast \mathrm{G} - \mathrm{B}\ast \mathrm{E}\right)\ast \mathrm{C}\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_4={\mathrm{J}}_7+{\mathrm{J}}_8}\end{array}$$

• (4)

  (S-F-S-S) Boundary Conditions:

$$\begin{array}{l}{\displaystyle \mathrm{A}=\left(-{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^2-{\mathrm{D}}_{12}\right)}\\ {\displaystyle \mathrm{B}=\left(-{\mathrm{D}}_{22}\ast {\mathsf{\alpha }}^2-{\mathrm{D}}_{12}\right)}\\ {\displaystyle \mathrm{C}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }-{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\right)}\\ {\displaystyle \mathrm{D}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+2{\mathrm{D}}_{66}{\mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }+{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\right)}\\ {\displaystyle \mathrm{E}= \cos \mathrm{h}\left(\mathsf{\alpha }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{F}= \sin \mathrm{h}\left(\mathsf{\alpha }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{G}= \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{H}= \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)}\\ {\displaystyle \mathrm{I}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\ast \sinh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{J}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\alpha }\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\alpha }}^3\ast \cosh \left(\mathsf{\alpha }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{K}=\left(-{\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)-2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)+{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\ast \sin \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{L}=\left({\mathrm{D}}_{12}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)+2{\mathrm{D}}_{66}{\ast \mathsf{\alpha }}_{\mathrm{m}}^2\ast \mathsf{\beta }\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)-{\mathrm{D}}_{22}{\ast \mathsf{\beta }}^3\ast \cos \left(\mathsf{\beta }\ast \mathrm{b}\right)\right)}\\ {\displaystyle \mathrm{M}=\frac{\tau }{\mathrm{a}\ast \mathrm{h}}}\\ {\displaystyle {\mathrm{J}}_1=\frac{\mathrm{B}\ast \left(\mathrm{F}\ast \mathrm{L} - \mathrm{H}\ast \mathrm{J}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J}-\mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_2=\frac{\mathrm{B}\ast \left(\mathrm{C}\ast \mathrm{H} - \mathrm{D}\ast \mathrm{F}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_1={\mathrm{J}}_1+{\mathrm{J}}_2}\\ {\displaystyle {\mathrm{J}}_3=\frac{\left(\mathrm{A}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{H}\ast \mathrm{K} - \mathrm{B}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{H}\ast \mathrm{I}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_4=\frac{\mathrm{D}\ast \left(\mathrm{A}\ast \mathrm{G} - \mathrm{B}\ast \mathrm{E}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_2={\mathrm{J}}_3-{\mathrm{J}}_4}\\ {\displaystyle {\mathrm{J}}_5=\frac{\mathrm{A}\ast \left(\mathrm{F}\ast \mathrm{L} - \mathrm{H}\ast \mathrm{J}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_6=\frac{\left(\mathrm{C}\ast \mathrm{H} - \mathrm{D}\ast \mathrm{F}\right)\ast \mathrm{A}\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_3=-{\mathrm{J}}_5-{\mathrm{J}}_6}\\ {\displaystyle {\mathrm{J}}_7=\frac{\left(\mathrm{A}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{G}\ast \mathrm{J} + \mathrm{B}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{F}\ast \mathrm{I}\right)\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{J}}_8=\frac{\left(\mathrm{A}\ast \mathrm{G} - \mathrm{B}\ast \mathrm{E}\right)\ast \mathrm{C}\ast \mathrm{M}}{\mathrm{A}\ast \mathrm{C}\ast \mathrm{G}\ast \mathrm{L} - \mathrm{A}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{K} + \mathrm{A}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{K} - \mathrm{A}\ast \mathrm{D}\ast \mathrm{G}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{C}\ast \mathrm{E}\ast \mathrm{L} + \mathrm{B}\ast \mathrm{C}\ast \mathrm{H}\ast \mathrm{I} + \mathrm{B}\ast \mathrm{D}\ast \mathrm{E}\ast \mathrm{J} - \mathrm{B}\ast \mathrm{D}\ast \mathrm{F}\ast \mathrm{I}}}\\ {\displaystyle {\mathrm{N}}_4={\mathrm{J}}_7+{\mathrm{J}}_8}\end{array}$$

References

1. Keshav, V.; Patel, S.N. Non-Linear dynamic pulse buckling of laminated composite curved panels. J. Struct. Eng. Mech.; 2020; 73, pp. 181-190.

2. Mondal, S.; Ramachandra, L.S. Nonlinear Dynamic Pulse Buckling of Imperfect Laminated Composite Plate with Delamination. Int. J. Solids Struct.; 2020; 198, pp. 170-182.

3. Taczała, M.; Buczkowski, R.; Kleiber, M. Nonlinear buckling and post-buckling response of stiffened FGM plates in thermal environments. J. Compos. Part B Eng.; 2017; 109, pp. 238-247.

4. Tao, Z.; Tu-guang, L.; Yao, Z.; Luo, J. Nonlinear dynamic buckling of stiffened plates under in-plane impact load. J. Zhejiang Univ. Sci. A; 2004; 5, pp. 609-617.

5. Yousuf, L.S. Nonlinear dynamics investigation of flexural stiffness of composite laminated plate under the effect of temperature and combined loading using Lyapunov exponent parameter. J. Compos. Part B Eng.; 2021; 219, 108926. [DOI: https://dx.doi.org/10.1016/j.compositesb.2021.108926]

6. Yousuf, L.S. Nonlinear dynamics simulation of bending deflection for composite laminated plate under varied temperature using Lyapunov exponent parameter. Proceedings of the ASME 2021 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference 2021; Online, 17–19 August 2021.

7. Hegaze, M.H. Nonlinear Dynamic Analysis of Stiffened and Unstiffened Laminated Composite Plates Using a High-order Element. J. Compos. Mater.; 2009; 44, pp. 327-346. [DOI: https://dx.doi.org/10.1177/0021998309345356]

8. Less, H.; Abramovich, H. Dynamic buckling of a laminated composite stringer–stiffened cylindrical panel. J. Compos. Part B Eng.; 2012; 43, pp. 2348-2358. [DOI: https://dx.doi.org/10.1016/j.compositesb.2011.11.070]

9. Aboudi, J.; Gilat, R. The Lyapunov exponents as a quantitative criterion for the dynamic buckling of composite plates. Int. J. Solids Struct.; 2002; 39, pp. 467-481.

10. Chai, Y.Y.; Li, F.M.; Song, Z.G. Nonlinear Vibration Behaviors of Composite Laminated Plates with Time-Dependent Base Excitation and Boundary Conditions. Int. J. Nonlinear Sci. Numer. Simul.; 2017; 18, pp. 145-161. [DOI: https://dx.doi.org/10.1515/ijnsns-2016-0138]

11. Borkowski, L. Numerical analysis of dynamic stability of an isotropic plate by applying tools used in dynamics. Dyn. Syst. Theory Appl.; 2017; 248, pp. 63-73.

12. Kim, H.G. Experimental nonlinear dynamics and snap-through of post-buckled composite plates. Nonlinear Dyn.; 2017; 1, pp. 21-35.

13. Pan, N. Theoretical determination of the optimal fiber volume fraction and fiber-matrix property compatibility of short fiber composites. Polym. Compos.; 1993; 14, pp. 85-93. [DOI: https://dx.doi.org/10.1002/pc.750140202]

14. Gopalan, V.; Suthenthiraveerappa, V.; David, J.S.; Subramanian, J.; Annamalai, A.R.; Jen, C.P. Experimental and Numerical Analyses on the Buckling Characteristics of Woven Flax/Epoxy Laminated Composite Plate under Axial Compression. Polymers; 2021; 13, 995. [DOI: https://dx.doi.org/10.3390/polym13070995] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33805018]

15. Sun, C.T. Mechanics of Aircraft Structures; 2nd ed. Wiley: Hoboken, NJ, USA, 2021.

16. Thompsom, M.K.; Thompson, J.M. ANSYS Mechanical APDL for Finite Element Analysis; Butterworth-Heinemann: Oxford, UK, 2017.

17. Parlitz, U. Estimating lyapunov exponents from time series. Chaos Detection and Predictability; Springer: Berlin/Heidelberg, Germany, 2016; Volume 915, pp. 1-34.

18. Terrier, P.; Fabienne, R. Maximum Lyapunov exponent revisited: Long-term attractor divergence of gait dynamics is highly sensitive to the noise structure of stride intervals. Gait Posture; 2018; 66, pp. 236-241. [DOI: https://dx.doi.org/10.1016/j.gaitpost.2018.08.010] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30212783]

19. Hussain, V.S.; Spano, M.L.; Lockhart, T.E. Effect of data length on time delay and embedding dimension for calculating the Lyapunov exponent in walking. J. R. Soc. Interface; 2020; 17, 20200311. [DOI: https://dx.doi.org/10.1098/rsif.2020.0311] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32674711]

20. Garcia, M.M.; Morales, I.; Rodriguez, J.M.; Marin, M.R. Selection of embedding dimension and delay time in phase space reconstruction via symbolic dynamics. Entropy; 2021; 23, 221. [DOI: https://dx.doi.org/10.3390/e23020221] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33670103]

21. Stoica, P.; Moses, R.L. Spectral analysis of signals. Digital Signal Processing; Rao, K.D.; Swamy, M.N.S. Springer: Berlin/Heidelberg, Germany, 2005; pp. 721-751.