Analytic Buckling Solution of Composite Classical Rectangular Thin Plates Using Galerkin Variational Approach

Abstract

This paper focuses on investigating the mechanical buckling behavior of composite and orthotropic classical rectangular plates using the Galerkin theory. By utilizing the classical plate theory of elasticity, the governing equations for the functionally graded plate material subjected to a uniaxially distributed load along the x and y axes are derived based on the principle of work energy. Consequently, the displacement function is determined as the solution to a differential equation that satisfies specific boundary conditions. The equilibrium expression derived from the work-energy equation is then minimized and analytically solved to determine the stresses and critical buckling loads of an orthotropic composite rectangular plate. The critical buckling loads of all edges, including clamped and simply supported orthotropic and isotropic thin rectangular plates with aspect ratios ranging from 0.1 to 3, are calculated and compared with existing literature to assess the stiffness of isotropic and orthotropic plates. The findings reveal that the critical buckling load of a clamped plate decreases as the aspect ratio increases from 0.1 to 1 but increases for aspect ratios between 1 and 3. The average percentage variation of critical buckling load values obtained from previous studies using orthotropic and isotropic (classical lamination plate theory [CLPT]) materials is found to be 1.93% and 3.65%, respectively. Additionally, it is observed that the critical buckling load of clamped edges in isotropic plates is lower than that of orthotropic plates at lower aspect ratios (0.4–1.0) but higher at higher aspect ratios (above 1.0). The results indicate that plates with all-round clamped boundaries exhibit higher buckling loads due to their effective fixation at all boundaries, which provides greater resistance to deformation and buckling compared to plates with all edges simply supported, where there is more flexibility for deformation under loads. This highlights the influence of boundary conditions on the critical buckling load of the plate. Furthermore, a significant percentage difference is noted when comparing approximate solutions, while a low percentage difference of 0.43% with the author’s derived displacement function indicates close agreement, demonstrating the high level of accuracy and reliability of the proposed model in predicting buckling load of plates. Hence, it can be recommended for stability analysis, especially where the material properties need to be controlled in a specific direction to optimize performance.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Plates can be classified based on their geometric properties, such as being circular, quadrilateral, or rectangular [1–3], or their material properties, which can include being homogeneous or nonhomogeneous, isotropic, anisotropic, or orthotropic [4–6]. The analysis and design of plate structures are essential for ensuring stability in terms of stiffness and strength parameters. An isotropic plate exhibits uniform material properties in three mutually perpendicular directions [7–9], whereas an orthotropic or composite plate displays varying mechanical properties, such as stiffness and strength, along its length, breadth, and thickness [10–12]. Orthotropic and composite plates find widespread use in engineering and construction applications where a high level of structural performance is essential. They are particularly recommended for structures where relying on isotropic material analysis [13–15] may not be sufficient, as seen in the design of laminated composite materials, bridges, and aerospace structures [16, 17], etc.

The categorization of rectangular plates into thick, thin, and moderately thick classes based on the ratio of depth to the smallest side dimension was introduced by Timoshenko [18] in 1940. He classified rectangular plates into three groups: Class I for plates with a depth-to-side ratio below 0.1 (thick) [19–21], Class II for plates with a ratio exceeding 0.2 (thin) [22–24], and Class III for plates with a ratio ranging from 0.1 to 0.2 (moderately thick) [25–27]. However, research indicates that thin plates are more prone to buckling under compressive and shearing loads compared to thicker plates [28, 29]. Thin rectangular orthotropic or composite plates are commonly employed in structures like roof and floor systems, ship hulls, and buildings due to their lightweight, high strength-to-weight ratio, and corrosion resistance [30, 31]. They also play a significant role in constructing space structures such as satellites and spacecraft, where lightweight materials are essential for reducing overall structural weight [32–34]. The buckling load is influenced by the plate thickness; thinner plates exhibit lower buckling loads, underscoring the importance of studying orthotropic thin rectangular plates. The failure of thin composite plate elements is often attributed more to elastic instability than to a lack of strength, making the analysis of these plates crucial.

Elastic stability in a composite plate structure refers to its ability to maintain its original shape and resist deformation when subjected to a load [35–37]. The critical load, on the other hand, represents the maximum load that the structure can bear before collapsing under the applied force [38–40]. Therefore, when such a structure is subjected to an in-plane compressive load, a classical analysis of elastic stability and critical buckling load becomes essential to prevent buckling, deformation, and potential structural failure. In the design of thin plate structures, it is crucial to ensure adequate stiffness and strength to withstand in-plane compressive loads, thereby ensuring that the critical buckling load is not surpassed, thus averting instability and failure.

The energy approach [41] can take the form of either an indirect variational approach [42–44] or a direct variational approach [37, 45, 46]. In the direct variational approach, the energy function is minimized to establish equilibrium in force equations. Common examples of direct energy variational methods include the Rayleigh, Ritz, and Rayleigh-Ritz methods [47]. On the other hand, the indirect variational approach does not transform the energy function into a force function but relies on the conservation of energy principle. This means that the total energy (input and output) in a static equilibrium continuum always equals zero. Indirect energy variational methods include the finite difference method, boundary collocation methods, the boundary element method, and Galerkin’s method [48–50].

Classical plate theory is applicable when the plate’s thickness is significantly smaller than its length and width, and deformations in the depth direction are minimal compared to other axes of the plate. This theory is widely used in engineering due to its simplicity and ease of application compared to refined plate theory (RPT) [51–53] and 3-D analysis [54–56]. For thin plate analysis, RPT and 3D theory are often deemed unnecessary. Orthotropic plate materials are commonly used in various engineering applications, particularly in composite materials and complex structures where high precision is essential, eliminating the need for the intricate analytical procedures involved in refined plate theory and 3D plate theory. Previous researchers have employed various classical plate theories (classical lamination plate theory [CLPT]) and methods [57] to address structural analysis challenges related to buckling, bending, and vibration responses of rectangular plates.

Reddy [58], in their book “Theory of Plates and Shells,” presented a thorough analysis of buckling in thin orthotropic and isotropic plates using classical plate theory. They delve into the fundamental principles of plate theory and offer detailed derivations of the governing equations for buckling analysis. However, their method is limited to analyzing plates that are simply supported along the four edges and is challenging to apply to plates with all edges clamped. Crisfield [59] explores the use of nonlinear finite element analysis in studying the buckling characteristics of thin orthotropic and isotropic plates. He provides detailed derivations of the governing equations for plate buckling and discusses the utilization of numerical methods for analyzing buckling behavior.

Reddy [60] explores the use of finite element analysis in studying the buckling behavior of thin orthotropic and isotropic plates. He offers a detailed examination of the theory and numerical techniques employed for analyzing plate buckling, covering topics such as geometric and material nonlinearities. In a study by Ma and Wang [61] published in the Journal of Engineering Mechanics, the researchers conducted a buckling analysis of orthotropic rectangular plates with various edge conditions using the finite element method (FEM), while Wang et al. [62] employed the differential quadrature method. The research aimed to investigate the buckling behavior of orthotropic rectangular plates under different loading conditions and edge constraints. While these authors contribute to the understanding of analyzing orthotropic plates with varied boundary conditions and demonstrate the efficacy of the quadrature approach in predicting the structural behavior of plates, it was noted that all three studies [60–62] only provide an approximate solution and unlike analytical approach, have limitations in determining the deflection value at any point in the plate structure.

Ugural [63] and Wang and Yu [64] conducted a comprehensive analysis of the buckling behavior of thin orthotropic and isotropic plates using classical plate theory. They examined the impact of material anisotropy on plate buckling, offering practical examples and applications of the theory. Their work was found to provide valuable insights and practical applications for engineers and researchers in the field of structural mechanics. The authors likely detailed a methodology outlining the use of the differential quadrature method to model orthotropic material properties, various boundary conditions, and loading scenarios for rectangular plates. They may have discussed the significance of their findings using the differential quadrature method in analyzing the buckling behavior of orthotropic structures with diverse boundary conditions in thin-walled structures. Their results shed light on the buckling characteristics under uniaxial loading conditions and underscored the importance of the differential quadrature method in predicting the structural stability of such plates.

In a study by Gao and Li [65] published on composite structures, the researchers investigated the buckling of composite plates using the variational iteration method as a numerical technique for solving the governing equations. They explored the theoretical framework of the variational iteration method and its application in analyzing the buckling behavior of orthotropic plate structures.

Galerkin Method: The Galerkin method is an analytical approach that simplifies the mathematical modeling of complex structures. It relies on variational methods, such as the principle of minimum potential energy, to derive governing equations. It effectively uses orthogonal polynomial functions to address boundary conditions and solve differential equations analytically for specific problems [66].

FEM: FEM is a numerical technique that divides a continuum (like a plate) into smaller, simpler parts called elements. The governing equations are formulated for these elements and assembled into a global system of equations. This method allows for complex geometries and loading conditions, making it highly versatile in engineering applications [67].

Approach and Application of Galerkin and FEM: The Galerkin method, as applied in the study, focuses on specific cases like the behavior of functionally graded orthotropic rectangular plates under uniaxial loads. It enables the derivation of closed-form solutions, which are quick to compute and easy to understand, particularly for issues that fit well within its assumptions. Due to this, it is highly efficient for parameter studies and theoretical investigations, providing insights into buckling behavior without extensive computational resources. FEM excels in handling complex boundary conditions, material properties, and load configurations. It is ideal for real-world applications where plates may have arbitrary shapes and varying thicknesses and can incorporate nonlinear material behavior. While it can provide highly accurate results, the setup process, including meshing and convergence requirements, may be time-consuming [68].

Performance and Accuracy of Galerkin and FEM: The analytical nature of the Galerkin method allows exact solutions for specific bounded problems, leading to high accuracy in scenarios where the assumptions hold true. As noted in the study, the agreement of the results with existing literature suggests a reliable and accurate prediction of critical buckling loads [66]. The accuracy of FEM depends on the quality of the mesh and the chosen elements. Advanced techniques within FEM, such as adaptive meshing and refined element types, can result in highly accurate solutions but require careful handling and validation against analytical results or experimental data [69]. FEM is generally preferred for complex and three-dimensional structures due to its flexibility and capability to simulate real-world scenarios.

Computational Considerations: Requires relatively low computational effort and is suited for quick assessments of structural behavior [66]. The ease of deriving exact solutions facilitates rapid exploration of parameter variations, especially in research settings. Typically demands more computational resources, particularly for large-scale problems, as the process involves solving extensive systems of equations. However, advancements in computational capabilities have made FEM highly feasible for detailed analyses in contemporary engineering practice [67].

The literature review of previous studies highlighted in the present research reveals notable distinctions in terms of the methodology of analysis, the shape function of the plate, boundary conditions, and material properties. The past works discussed in the literature review demonstrated unique characteristics both individually and collectively when compared to the current study. Unlike some other composite orthotropic plate studies in the literature that relied on approximate solutions, the present work stands out for utilizing an exact deflection function derived from the principles of elasticity theory. In this study, a novel approach was adopted where a power series displacement equation in the form of orthogonal polynomials was employed to derive the deflection function. This method allowed for a more precise and detailed analysis of the plate’s behavior under various loading conditions. Moreover, the critical buckling load of orthotropic/composite classical rectangular plates was determined using the Galerkin variational principle, showcasing the rigorous analytical framework employed in this research. The manuscript specifically addresses the mechanical buckling behavior of functionally graded orthotropic composite plates, which adds value by targeting materials that have varying properties in different directions. This specificity provides insights into the behavior of advanced materials used in engineering applications.

Furthermore, following the determination of the critical buckling load, an isotropic equivalent was calculated in this study. This step is crucial in simplifying the complex behavior of the orthotropic plates to an equivalent isotropic material, facilitating a more straightforward analysis and interpretation of the results. Overall, the application of the power series displacement equation, the Galerkin variational principle, and the determination of the isotropic equivalent in this study represent a significant advancement in the analysis of composite classical thin rectangular plates. The study highlights how the Galerkin method, combined with classical plate theory and orthogonal polynomial shape functions, can significantly improve understanding and predictions of buckling loads in structures. The implications for structural performance and safety are substantial, especially in high-stakes applications like aerospace and civil engineering. The findings from this research provide valuable insights into the behavior of orthotropic and composite plate materials under various loading conditions and offer a foundation for further research in the field of structural mechanics.

2. Methodology

2.1. Assumptions of the Analysis and Theoretical Framework

In general, the choice between Kirchhoff’s and Timoshenko’s assumptions [28, 58] will depend on the plate’s thickness, the expected shear effects, and the specific application. In this work, thin rectangular isotropic plates, Kirchhoff’s theory was used, and the following assumptions were considered:

a. The plane sections remain plane before and after deformation.

b. The mid-plane of the plate does not experience any deformation (zero transverse shear).

c. The deflections (w) are small, leading to linear behavior.

d. The plane sections remain plane and perpendicular to the neutral surface after deformation.

e. Applicable to thin plates where the thickness is small compared to width and length.

f. Only buckling deformations for thin plates are considered, and shear effects are minimal or negligible.

g. The theory used is applicable to thin rectangular plates subjected to uniform loading conditions.

h. Plates that have uniform mechanical properties in all directions within each ply.

i. Assumed to behave according to classical bending theory, which neglects transverse shear deformations, leading to a simplified analysis applicable to certain engineering applications.

Figure 1 shows a typical case of buckling of an elastic of the thin rectangular flat plate under direct and shear in-plane loading, which was axially loaded along the direction of x and y coordinates; when these loads cause the plate material to deform the plate, it results in buckling of the plate

[figure(s) omitted; refer to PDF]

It should be noted that the buckling of plates at constant thickness and made of orthotropic material, as shown in Figure 1, Kirchhoff’s hypotheses for an isotropic plate would not be adequate. Assuming that the principal directions of orthotropic coincide with the x and y coordinate axes which, in turn, lie in the middle plane of the plate. Thus, the stress– strain relations in Equations for isotropic plates, according to Ventsel and Krauthammer [4], are not valid for orthotropic plates.

2.2. Equation of Compatibility for Orthotropic Plate Under Uniformly Distributed Lateral and In-Plane Loadings

It is a commonly accepted fact that deformations are influenced by two distinct elastic constants, such as E and μ. Numerous construction materials like steel, aluminum, and others belong to this group. In such cases, the buckling of a plate under lateral load would occur in conjunction with its buckling under in-plane loading. The independent elastic constants in equation, G remain consistent for both isotropic and orthotropic materials, and the connection with these independent elastic constants is referred to as follows:$$\begin{array}{}\text{(1)}& \frac{{\mu }_x}{E_x}=\frac{{\mu }_y}{E_y}.\end{array}$$

Figure 2shows the following direct and shear in-plane loading acting on the thin rectangular plate:$$\begin{array}{}\text{(2)}& {N_x}^{\ast }=N_x+\frac{\partial N_x}{\partial x}dx,\end{array}$$$$\begin{array}{}\text{(3)}& {N_{xy}}^{\ast }=N_{xy}+\frac{\partial N_{xy}}{\partial x}dx,\end{array}$$$$\begin{array}{}\text{(4)}& {N_y}^{\ast }=N_y+\frac{\partial N_y}{\partial y}dy,\end{array}$$$$\begin{array}{}\text{(5)}& {N_{yx}}^{\ast }=N_{yx}+\frac{\partial N_{yx}}{\partial y}dy.\end{array}$$

[figure(s) omitted; refer to PDF]

Since the transverse shear forces are of negligible quantity, the reactive stresses are presented as follows:$$\begin{array}{}\text{(6)}& {\sigma }_x=-\frac{E_{x}Z}{1-{\mu }_x{\mu }_y}\frac{{\partial }^{2}w}{\partial x^2}+{\mu }_y\frac{{\partial }^{2}w}{\partial y^2},\end{array}$$$$\begin{array}{}\text{(7)}& {\sigma }_y=-\frac{E_{y}Z}{1-{\mu }_x{\mu }_y}\frac{{\partial }^{2}w}{\partial y^2}+{\mu }_x\frac{{\partial }^{2}w}{\partial x^2},\end{array}$$$$\begin{array}{}\text{(8)}& {\tau }_{xy}=-2G\frac{Z{\partial }^{2}w}{\partial x\partial y},\end{array}$$$$\begin{array}{}\text{(9)}& {\epsilon }_x=-\frac{Z{\partial }^{2}w}{\partial x^2},\end{array}$$$$\begin{array}{}\text{(10)}& {\epsilon }_y=-\frac{Z{\partial }^{2}w}{\partial y^2},\end{array}$$$$\begin{array}{}\text{(11)}& {\gamma }_{xy}=-2\frac{Z{\partial }^{2}w}{\partial x\partial y},\end{array}$$$$\begin{array}{}\text{(12)}& D_x=\frac{E_x}{1-{\mu }_x{\mu }_y}\cdot \frac{h^3}{12},\end{array}$$$$\begin{array}{}\text{(13)}& D_y=\frac{E_y}{1-{\mu }_x{\mu }_y}\cdot \frac{h^3}{12},\end{array}$$$$\begin{array}{}\text{(14)}& D_{xy}=\frac{E_x{\mu }_y}{1-{\mu }_x{\mu }_y}\cdot \frac{h^3}{12},\end{array}$$$$\begin{array}{}\text{(15)}& D_{yx}=\frac{E_y{\mu }_x}{1-{\mu }_x{\mu }_y}\cdot \frac{h^3}{12},\end{array}$$$$\begin{array}{}\text{(16)}& D_s=\frac{{Gh}^3}{12}.\end{array}$$

Given that;$$\begin{array}{}\text{(17)}& G\cong \frac{\sqrt{E_x{E}_y}}{21+\sqrt{{\mu }_x{\mu }_y}}.\end{array}$$

The elastic constants $E_x,E_y,{\mu }_x,{\mu }_y$, and G are considered to characterize an orthotropic material. Specifically, $E_x,E_y,{\mu }_x,{\mu }_y$ represent the elasticity moduli and Poisson’s ratios in the x and y directions, while G denotes the shear modulus. Additionally, $D_x$, $D_y$, $D_{xy}$, and $D_{yx}$ signify the flexural rigidities, and D_s denotes the torsional rigidity of an orthotropic plate.

First, by analyzing Figure 2 and projecting all forces along the x-direction while considering the St. Venant principle, we can conclude that if a body or parent body is in equilibrium, any isolated part of that body will also be in equilibrium. In this case, we assume that our parent thin rectangular plate is in equilibrium, leading to the derivation of the following equations. Thus, projecting all forces in y direction $\sum F_x=0$, gave the following:$$\begin{array}{}\text{(18)}& N_x+\frac{\partial N_x}{\partial x}dxdy-N_{x}dy+N_{yx}+\frac{\partial N_{yx}}{\partial y}dydx-N_{yx}dx=0.\end{array}$$

Solving Equation (18), we obtain the following:$$\begin{array}{}\text{(19)}& \frac{\partial N_x}{\partial x}dxdy+\frac{\partial N_{yx}}{\partial y}dydx=0.\end{array}$$

Similarly, projecting all forces in y direction $\sum F_y=0,$ we shall obtain the following:$$\begin{array}{}\text{(20)}& N_y+\frac{\partial N_y}{\partial y}dydx-N_{y}dx+N_{xy}+\frac{\partial N_{xy}}{\partial x}dxdy-N_{xy}dy=0.\end{array}$$

Solving Equation (20), we have the following:$$\begin{array}{}\text{(21)}& \frac{\partial N_y}{\partial y}dydx+\frac{\partial N_{xy}}{\partial x}dxdy=0.\end{array}$$

The plate in Figure 3, after bending, will show how all forces are projected along z direction $\sum F_z=0$. But imagine an element of the plate, cut in one direction, say the x-axis.

[figure(s) omitted; refer to PDF]

Hence,$$\begin{array}{}\text{(22)}& {\theta }_x=\frac{\partial w}{\partial x}=\tan {\theta }_x.\end{array}$$

The tangent of a small angle is approximately equal to the angle itself (small deflection of plate).$$\begin{array}{}\text{(23)}& {\theta }_x=\frac{\partial w}{\partial x},\end{array}$$$$\begin{array}{}\text{(24)}& {\theta }_x+\frac{\partial {\theta }_x}{\partial x}dx=\frac{\partial w}{\partial x}+\frac{\partial }{\partial x}\frac{\partial w}{\partial x}dx.\end{array}$$

Similarly, projecting these membrane forces N_x and ${N_x}^{\ast }$ in z–x direction, that is,$$\begin{array}{}& \sum F_{zx}=0,\end{array}$$$$\begin{array}{}\text{(25)}& N_x\sin {\theta }_{x}dy-N_x+\frac{\partial N_x}{\partial x}dx\sin {\theta }_x+\frac{\partial {\theta }_x}{\partial x}dxdy=0,\end{array}$$

Recall that the small angle is equal to the angle as follows:$$\begin{array}{}\text{(26)}& {\theta }_x=\frac{\partial w}{\partial x}=\sin {\theta }_x,\end{array}$$

Thus,$$\begin{array}{}\text{(27)}& \sin {\theta }_x+\frac{\partial {\theta }_x}{\partial x}dx=N_x\frac{\partial }{\partial x}\frac{\partial w}{\partial x}dxdy+\frac{\partial N_x}{\partial x}\frac{\partial w}{\partial x}dxdy.\end{array}$$

Projecting all membrane forces $N_x$ and ${N_x}^{\ast }$ along the z–x axis considering Equation (26), Equation (27) becomes the following:$$\begin{array}{}\text{(28)}& N_x\frac{\partial w}{\partial x}dy-N_x\frac{\partial w}{\partial x}dy-N_x\frac{\partial }{\partial x}\frac{\partial w}{\partial x}dxdy-\frac{\partial N_x}{\partial x}dx\frac{\partial w}{\partial x}dy-\frac{\partial N_x}{\partial x}dx\frac{\partial }{\partial x}\frac{\partial w}{\partial x}dxdy=0.\end{array}$$

Simplifying Equation (28) gives the following:$$\begin{array}{}\text{(29)}& -N_x\frac{\partial }{\partial x}\frac{\partial w}{\partial x}dxdy-\frac{\partial N_x}{\partial x}dx\frac{\partial w}{\partial x}dy-\frac{\partial N_x}{\partial x}dx\frac{\partial }{\partial x}\frac{\partial w}{\partial x}dxdy=0.\end{array}$$

$\mathrm{T}\mathrm{a}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{g}dx\times dx=d{x}^2=0\text{square\hspace{0.17em}of\hspace{0.17em}very\hspace{0.17em}small\hspace{0.17em}value}\cong 0$ and rearranging gives the following:$$\begin{array}{}\text{(30)}& N_x\frac{\partial }{\partial x}\frac{\partial w}{\partial x}dxdy+\frac{\partial N_x}{\partial x}\frac{\partial w}{\partial x}dxdy.\end{array}$$

Also, consider the plate in Figure 1 after bending cut into a beam strip along the y direction with the adjacent side as dx presented in Figure 4.

[figure(s) omitted; refer to PDF]

From Figure 4, considering the tangent of a small angle equal to the angle itself,$$\begin{array}{}\text{(31)}& {\theta }_y=\frac{\partial w}{\partial y}.\end{array}$$

Hence,$$\begin{array}{}\text{(32)}& {\theta }_y+\frac{\partial {\theta }_y}{\partial y}dy=\frac{\partial w}{\partial y}+\frac{\partial }{\partial y}\frac{\partial w}{\partial y}dy.\end{array}$$

Projecting all membrane forces $N_y$ and ${N_y}^{\ast }$along z–y axis, that is,$$\begin{array}{}& \sum F_{zy}=0,\end{array}$$$$\begin{array}{}\text{(33)}& N_y\sin {\theta }_{y}dx-N_y+\frac{\partial N_y}{\partial y}dy\sin {\theta }_y+\frac{\partial {\theta }_y}{\partial y}dydx=0.\end{array}$$

This is the same as follows:$$\begin{array}{}\text{(34)}& N_y{\theta }_{y}dx-N_y+\frac{\partial N_y}{\partial y}dy{\theta }_y+\frac{\partial {\theta }_y}{\partial y}dydx=0.\end{array}$$

Substituting Equation (31) into Equation (34) and expand gives the following:$$\begin{array}{}\text{(35)}& N_y\frac{\partial w}{\partial y}dx-N_y\frac{\partial w}{\partial y}dx-N_y\frac{\partial }{\partial y}\frac{\partial w}{\partial y}dydx-\frac{\partial N_y}{\partial y}dy\frac{\partial w}{\partial y}dx-\frac{\partial N_y}{\partial y}dy\frac{\partial }{\partial y}\frac{\partial w}{\partial y}dydx=0.\end{array}$$

$\mathrm{T}\mathrm{a}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{g}dy\times dy=d{y}^2=0\text{square\hspace{0.17em}of\hspace{0.17em}very\hspace{0.17em}small\hspace{0.17em}value}\cong 0$, and rearranging gives the following:$$\begin{array}{}\text{(36)}& N_y\frac{\partial }{\partial y}\frac{\partial w}{\partial y}dydx+\frac{\partial N_y}{\partial y}\frac{\partial }{\partial x}\frac{\partial w}{\partial y}dxdy=0.\end{array}$$

Consider the shear forces on the cut deflected element in both directions, x and y, as follows:$$\begin{array}{}\text{(37)}& {\theta }_x^{\ast }=\frac{\partial w}{\partial x}+\frac{\partial }{\partial x}\frac{\partial w}{\partial x}dx,\end{array}$$$$\begin{array}{}\text{(38)}& {\theta }_y^{\ast }=\frac{\partial w}{\partial y}+\frac{\partial }{\partial y}\frac{\partial w}{\partial y}dy.\end{array}$$

Projecting the N_xy and N_yx along the z-axis, taking N_xy = N_yx, that is,$$\begin{array}{}& \sum F_{xy}=0,\end{array}$$$$\begin{array}{}\text{(39)}& N_{xy}\frac{\partial w}{\partial y}dy-N_{xy}+\frac{\partial N_{xy}}{\partial x}dx\frac{\partial w}{\partial y}+\frac{\partial }{\partial x}\frac{\partial w}{\partial x}dxdy+N_{yx}\frac{\partial w}{\partial x}dx=0.\end{array}$$

Expanding Equation (39) and rearranging gives the following:$$\begin{array}{}\text{(40)}& 0=-N_{xy}\frac{\partial }{\partial x}\frac{\partial w}{\partial y}dxdy-\frac{\partial N_{xy}}{\partial x}dx\frac{\partial w}{\partial y}dy-N_{yx}\frac{\partial }{\partial y}\frac{\partial w}{\partial x}dxdy-\frac{\partial N_{yx}}{\partial y}dy\frac{\partial w}{\partial x}dx.\end{array}$$

Taking N_xy = N_yx and simplifying gives the following:$$\begin{array}{}\text{(41)}& 2N_{xy}\frac{\partial }{\partial y}\frac{\partial w}{\partial x}dxdy+\frac{\partial N_{xy}}{\partial x}dy\frac{\partial w}{\partial y}dx+\frac{\partial N_{xy}}{\partial y}dy\frac{\partial w}{\partial x}dx=0.\end{array}$$

Adding Equations (30), (36), and (41) of forces in z-direction we obtain the following:$$\begin{array}{}\text{(42)}& \sum F_z=\sum F_{zx}+\sum F_{zy}+\sum F_{zxy}=0,\end{array}$$$$\begin{array}{}\text{(43)}& N_x\frac{\partial }{\partial x}\frac{\partial w}{\partial x}dxdy+\frac{\partial N_x}{\partial x}\frac{\partial w}{\partial x}dxdy+N_y\frac{\partial }{\partial y}\frac{\partial w}{\partial y}dydx+\frac{\partial N_y}{\partial y}\frac{\partial }{\partial x}\frac{\partial w}{\partial y}dxdy+2N_{xy}\frac{\partial }{\partial y}\frac{\partial w}{\partial x}dxdy+\frac{\partial N_{xy}}{\partial x}dy\frac{\partial w}{\partial y}dx+\frac{\partial N_{xy}}{\partial y}\frac{\partial w}{\partial x}dxdy=0.\end{array}$$

Substituting Equations (19) and (21) into Equation (43) and simplifying gives the following:$$\begin{array}{}\text{(44)}& N_x\frac{{\partial }^{2}w}{\partial x^2}+2N_{xy}\frac{{\partial }^{2}w}{\partial x\partial y}+N_y\frac{{\partial }^{2}w}{\partial y^2}=0.\end{array}$$

2.3. Governing Equation for Deflection of Composite Rectangular Plates

The stress components, including stress resultants and stress couples, typically differ from one point to another in a loaded plate, and their variations are determined by static equilibrium conditions. By examining the compatibility of an element dxdy of an orthotropic plate under a vertically distributed load of intensity p (x, y) applied to the upper surface of the plate, the governing equation for plate deflection is derived. This is achieved by establishing three independent equilibrium conditions that account for the variations in moments and shear forces over a small element of a flat plate.

The force summation in the z axis gives the following:$$\begin{array}{}\text{(45)}& \frac{\partial Q_x}{\partial x}dxdy+\frac{\partial Q_y}{\partial y}dydx+Pdxdy=0.\end{array}$$

The moment summation about the X-axis leads to the following:$$\begin{array}{}\text{(46)}& \frac{\partial M_{xy}}{\partial x}dxdy+\frac{\partial M_y}{\partial y}dydx-Q_{y}dxdy=0,\end{array}$$$$\begin{array}{}\text{(47)}& \frac{\partial M_{yx}}{\partial y}dydx+\frac{\partial M_x}{\partial x}dxdy-Q_{x}dxdy=0.\end{array}$$

From Equations (46) and (47), the shear forces $Q_x$ and $Q_y$ can be expressed in terms of the moments in Equations (48) and (49), respectively,$$\begin{array}{}\text{(48)}& Q_x=\frac{\partial M_x}{\partial x}+\frac{\partial M_{xy}}{\partial y},\end{array}$$$$\begin{array}{}\text{(49)}& Q_y=\frac{\partial M_{xy}}{\partial x}+\frac{\partial M_y}{\partial y}.\end{array}$$

Substituting Equations (48) and (49) into Equation (45), taking into account that,

$M_{xy}=M_{yx}$, we obtain the following:$$\begin{array}{}\text{(50)}& \frac{{\partial }^2{M}_x}{\partial x^2}+2\frac{{\partial }^2{M}_{xy}}{\partial xdy}+\frac{{\partial }^2{M}_y}{\partial y^2}=-Px,y.\end{array}$$

The stress–strain relations of Equations (6)–(11) that reflect the orthotropic properties of the plate are given in Equations (51)–(53) as follows:$$\begin{array}{}\text{(51)}& {\epsilon }_x=\frac{{\sigma }_x}{E_x}-{\mu }_y\frac{{\sigma }_y}{E_y},\end{array}$$$$\begin{array}{}\text{(52)}& {\epsilon }_y=\frac{{\sigma }_y}{E_y}-{\mu }_x\frac{{\sigma }_x}{E_x},\end{array}$$$$\begin{array}{}\text{(53)}& {\gamma }_{xy}=\frac{{\tau }_{xy}}{G}.\end{array}$$

Solving Equations (51)–(53) for the stress components and taking into account Equation (1), we obtain the following:$$\begin{array}{}\text{(54)}& {\sigma }_x=\frac{E_x}{1-{\mu }_x{\mu }_y}{\epsilon }_x+{\mu }_y{\epsilon }_y,\end{array}$$$$\begin{array}{}\text{(55)}& {\sigma }_y=\frac{E_y}{1-{\mu }_x{\mu }_y}{\epsilon }_y+{\mu }_x{\epsilon }_x,\end{array}$$$$\begin{array}{}\text{(56)}& {\tau }_{xy}=G{\gamma }_{xy}.\end{array}$$

Considering the plate assumption of Timoshenko and Woinowsky-Krieger [28], the corresponding strain components for an isotropic material become the following:$$\begin{array}{}\text{(57)}& {\epsilon }_x=-\frac{Z{\partial }^{2}w}{\partial x^2},\end{array}$$$$\begin{array}{}\text{(58)}& {\epsilon }_y=-\frac{Z{\partial }^{2}w}{\partial y^2},\end{array}$$$$\begin{array}{}\text{(59)}& {\gamma }_{xy}=-2\frac{Z{\partial }^{2}w}{\partial x\partial y}.\end{array}$$

Substituting Equations (57), (58), and (59) expressions into Equations (54), (55), and (56), respectively, we obtain the following:$$\begin{array}{}\text{(60)}& {\sigma }_x=-\frac{E_{x}Z}{1-{\mu }_x{\mu }_y}\frac{{\partial }^{2}w}{\partial x^2}+{\mu }_y\frac{{\partial }^{2}w}{\partial y^2},\end{array}$$$$\begin{array}{}\text{(61)}& {\sigma }_y=-\frac{E_{y}Z}{1-{\mu }_x{\mu }_y}\frac{{\partial }^{2}w}{\partial y^2}+{\mu }_x\frac{{\partial }^{2}w}{\partial x^2},\end{array}$$$$\begin{array}{}\text{(62)}& {\tau }_{xy}=-2G\frac{Z{\partial }^{2}w}{\partial x\partial y}.\end{array}$$

Recalling that,$$\begin{array}{}\text{(63)}& M_i={\int }_{-\frac{h}{2}}^{+\frac{h}{2}}{\sigma }_{i}zdzi=x,y,xy.\end{array}$$

Therefore, the bending and twisting moments can then be expressed in terms of stress components given as follows:$$\begin{array}{}\text{(64)}& M_x={\int }_{-\frac{h}{2}}^{+\frac{h}{2}}{\sigma }_{x}zdz,\end{array}$$$$\begin{array}{}\text{(65)}& M_y={\int }_{-\frac{h}{2}}^{+\frac{h}{2}}{\sigma }_{y}z,\end{array}$$$$\begin{array}{}\text{(66)}& M_{xy}={\int }_{-\frac{h}{2}}^{+\frac{h}{2}}{\tau }_{xy}zdz.\end{array}$$

Substituting Equations (60), (61), and (62) into Equations (64), (65), and (66), respectively, and integrating over the plate thickness yields the following bending and twisting moment’s relations for orthotropic plates:$$\begin{array}{}\text{(67)}& M_x=-D_x\frac{{\partial }^{2}w}{\partial x^2}+D_{xy}\frac{{\partial }^{2}w}{\partial y^2},\end{array}$$$$\begin{array}{}\text{(68)}& M_y=-D_y\frac{{\partial }^{2}w}{\partial y^2}+D_{yx}\frac{{\partial }^{2}w}{\partial x^2},\end{array}$$$$\begin{array}{}\text{(69)}& M_{xy}=-2D_s\frac{{\partial }^{2}w}{\partial x\partial y},\end{array}$$where$$\begin{array}{}\text{(70)}& D_{xy}=D_{yx};H=D_{xy}+2D_s.\end{array}$$

Substituting Equations (67), (68), and (69) into Equation (50) and taking into account Equation (1) and simplifying, we have the following:$$\begin{array}{}\text{(71)}& -D_x\frac{{\partial }^{4}w}{\partial x^4}-D_{xy}\frac{{\partial }^{4}w}{\partial x^2\partial y^2}-4D_s\frac{{\partial }^{4}w}{\partial x^2\partial y^2}-D_y\frac{{\partial }^{4}w}{\partial y^4}-D_{yx}\frac{{\partial }^{4}w}{\partial x^2\partial y^2}=-Px,y.\end{array}$$

Factorizing, we obtain the following:$$\begin{array}{}\text{(72)}& D_x\frac{{\partial }^{4}w}{\partial x^4}+2H\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+D_y\frac{{\partial }^{4}w}{\partial y^4}=-Px,y,\end{array}$$or$$\begin{array}{}\text{(73)}& D_x\frac{{\partial }^{4}w}{\partial x^4}+2D_{xy}+2D_s\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+D_y\frac{{\partial }^{4}w}{\partial y^4}=-Px,y,\end{array}$$where P is the in-plane forces as given in Equation (45) for buckled plate due to in-plane loading; thus, Equation (44) becomes the following:$$\begin{array}{}\text{(74)}& D_x\frac{{\partial }^{4}w}{\partial x^4}+2H\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+D_y\frac{{\partial }^{4}w}{\partial y^4}=-N_x\frac{{\partial }^{2}w}{\partial x^2}+2N_{xy}\frac{{\partial }^{2}w}{\partial x\partial y}+N_y\frac{{\partial }^{2}w}{\partial y^2}.\end{array}$$

For a plate loaded only in one direction, b along the x-direction as follows:$$\begin{array}{}\text{(75)}& N_{xy}=N_y=0\text{and}{N}_x\ne 0.\end{array}$$

Substituting Equation (75) into Equation (74) and simplifying, we have the following:$$\begin{array}{}\text{(76)}& D_x\frac{{\partial }^{4}w}{\partial x^4}+2H\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+D_y\frac{{\partial }^{4}w}{\partial y^4}+N_x\frac{{\partial }^{2}w}{\partial x^2}=0.\end{array}$$

Equation (76) is the governing differential equation of buckled plate loaded along the edge perpendicular to the x-direction.

2.4. Application of Galerkin’s Method in the Determination of Unknown Functions for Buckling Material

The Galerkin variation approach was utilized to address nonlinear partial differential equations (PDEs) through the application of spatial discretization and a weighted residual formulation. This approach involves transforming the governing PDE (strong form) into an integral equation (weak form) and subsequently treating it with variations to derive a solution in the form of a system of matrix equations. This process allows for determining the unknown functions related to the buckling behavior of materials.

Let a differential equation of a given 2D boundary value problem be of the form as follows:$$\begin{array}{}\text{(77)}& LWx,y=-Px,y\text{in\hspace{0.17em}some}2D\text{domain}\Omega .\end{array}$$

An approximate solution of Equation (76) is sought in the following form:$$\begin{array}{}\text{(78)}& W_{N}x,y=\sum _{i=1}^N{K}_i{\bar{W}}_{i}x,y.\end{array}$$

Considering only the first term of the series,$$\begin{array}{}\text{(79)}& Wx,y=K\bar{W}x,y,\end{array}$$where $K$ = unknown coefficients to be determined. w = w (x, y) is an unknown linearly independent function of two variables that satisfy all the prescribed support conditions on the boundary domain but do not necessarily satisfy Equation (76). $P=$ a given load term causing buckling. $L=$ a symbol indicating either a linear or nonlinear governing differential equation.

From Calculus, any two functions $f_{1}x,f_{2}x$ are called mutually orthogonal in the interval (a, b) if they satisfy the condition:$$\begin{array}{}\text{(80)}& {\int }_a^b{f}_{1}x{f}_{2}xdx=0,\end{array}$$$$\begin{array}{}\text{(81)}& {\iint }_{A}L{w}_N-P\bar{W}x,ydxdy=0.\end{array}$$

Substituting the governing equation in Equation (76) as the deformation energy equation$\mathrm{L}{\mathrm{w}}_{\mathrm{N}}$, we have the following:$$\begin{array}{}\text{(82)}& {\iint }_A{D}_x\frac{{\partial }^{4}w}{\partial x^4}+2H\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+D_y\frac{{\partial }^{4}w}{\partial y^4}+N_x\frac{{\partial }^{2}w}{\partial x^2}\bar{W}x,ydxdy=0.\end{array}$$

Putting the potential energy equation in Equation (82) and the established external work for bucking $P={\iint }_A{N}_x\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{x}^2}\overline{W}x,ydxdy$ into Equation (81), we have the following:$$\begin{array}{}\text{(83)}& {\iint }_A{D}_x\frac{{\partial }^{4}w}{\partial x^4}+2H\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+D_y\frac{{\partial }^{4}w}{\partial y^4}+N_x\frac{{\partial }^{2}w}{\partial x^2}\bar{W}x,ydxdy=-{\iint }_A{N}_x\frac{{\partial }^{2}w}{\partial x^2}\bar{W}x,ydxdy.\end{array}$$

Simplifying Equation (83) by making $N_x$ subject of the expression, we have the following:$$\begin{array}{}\text{(84)}& N_x=-\frac{{\iint }_A{D}_x\frac{{\partial }^{4}w}{\partial x^4}+2H\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+D_y\frac{{\partial }^{4}w}{\partial y^4}{\bar{W}}_{i}x,ydxdy}{{\iint }_A\frac{{\partial }^{2}w}{\partial x^2}\bar{W}x,ydxdy}.\end{array}$$

Substituting Equation (84) into Equation (79) and simplifying, we have the following:$$\begin{array}{}\text{(85)}& N_x=-\frac{K{\iint }_A{D}_x\frac{{\partial }^4\bar{W}}{\partial x^4}+2H\frac{{\partial }^4\bar{W}}{\partial x^2\partial y^2}+D_y\frac{{\partial }^4\bar{W}}{\partial y^4}{\bar{W}}_{i}x,ydxdy}{K{\iint }_A\frac{{\partial }^2\bar{W}}{\partial x^2}\bar{W}x,ydxdy}.\end{array}$$

Given that the nondimensional relation is as follows:$$\begin{array}{}\text{(86)}& x=aR,y=bQ\mathrm{a}\mathrm{n}\mathrm{d}P=\frac{a}{b}.\end{array}$$

Substitute Equation (86) into Equation (85) gives the following:$$\begin{array}{}\text{(87)}& N_x=-\frac{{\iint }_A{D}_x\frac{{\partial }^4\bar{W}}{a^4\partial R^4}+2H\frac{{\partial }^4\bar{W}}{a^2{b}^2\partial R^2\partial Q^2}+D_y\frac{{\partial }^4\bar{W}}{b^4\partial Q^4}\bar{W}R,QabdRdQ}{{\iint }_A\frac{1}{a^2}\frac{{\partial }^2\bar{W}}{\partial R^2}\bar{W}R,QabdRdQ}.\end{array}$$

That is,$$\begin{array}{}\text{(88)}& N_x=-\frac{{\iint }_A{D}_x\frac{{\partial }^4\bar{W}}{P^4\partial R^4}+2H\frac{{\partial }^4\bar{W}}{P^2\partial R^2\partial Q^2}+D_y\frac{{\partial }^4\bar{W}}{\partial Q^4}\bar{W}R,QdRdQ}{{\iint }_A\frac{1}{b^2{P}^2}\frac{{\partial }^2\bar{W}}{\partial R^2}\bar{W}R,QdRdQ}.\end{array}$$

In nondimensional form, Equation (63) is simplified as follows:$$\begin{array}{}\text{(89)}& N_x=-\frac{P^2}{b^2}\frac{{\iint }_A{D}_x\frac{{\partial }^4\bar{W}}{P^4\partial R^4}+2H\frac{{\partial }^4\bar{W}}{P^2\partial R^2\partial Q^2}+D_y\frac{{\partial }^4\bar{W}}{\partial Q^4}\bar{W}R,QdRdQ}{{\iint }_A\frac{{\partial }^2\bar{W}}{\partial R^2}\bar{W}R,QdRdQ}.\end{array}$$

Equation (64) is the critical load of buckled plate loaded along the edge perpendicular to the x-direction.

2.5. Application of Polynomials Displacement Functions for the Solution of Particular Solution

The Gram–Schmidt method was employed to determine the precise deflection shape functions of the thin rectangular plate using a polynomial order. Imagine a beam with a random support configuration experiencing a uniformly distributed varying load in any direction, as illustrated in Figure 5. Due to this load, reactive forces, moments, and reactions would arise at the supports.

[figure(s) omitted; refer to PDF]

The equation of the moment of the beam at a section says t would be given as follows:$$\begin{array}{}\text{(90)}& M_t=R_1\cdot t-\frac{q\cdot t^2}{2}-M_1.\end{array}$$

Then, employing the elastic beam equation as follows:$$\begin{array}{}\text{(91)}& M_t=-D\frac{d^{2}W}{d{t}^2}.\end{array}$$

Hence, equating Equations (90) and (91) give the following:$$\begin{array}{}\text{(92)}& D\frac{d^{2}W}{d{t}^2}=\frac{q\cdot t^2}{2}+M_1-R_1\cdot t.\end{array}$$

Obtaining the deflection function by integrating Equation (92) twice with respect to the arbitrary direction, t, to get the following:$$\begin{array}{}\text{(93)}& W_t=C_o+C_1\cdot t+C_2{t}^2+C_3\cdot t^3+C_4\cdot t^4,\end{array}$$where $C_o\text{and}{C}_1$ are constants of integration and$$\begin{array}{}\text{(94)}& C_4=\frac{q}{24D};C_3=\frac{-R}{6D};C_2=\frac{M_1}{2D}.\end{array}$$

Therefore, in Equation (93), the highest polynomial is set at 4 for a uniformly distributed load. This implies that a fourth-order function would be appropriate for a uniformly varying load. The deflection function derived from the orthogonal polynomial characteristics is as follows:$$\begin{array}{}\text{(95)}& Wx,y=\sum _{m=1}^{\mathrm{\infty }}{A}_m{R}^m.\sum _{n=1}^{\mathrm{\infty }}{B}_n{Q}^n.\end{array}$$

By expanding Equation (95) to the fourth series in the direction of the R and Q axis of the plate and their functions are denoted by A_m and B_n, respectively, we obtain the following:$$\begin{array}{}\text{(96)}& Wx,y=\sum _{m=1}^4{A}_m{R}^m.B_n{Q}^n=A_0+A_{1}R+A_2{R}^2+A_3{R}^3+A_4{R}^4.B_0+B_{1}Q+B_2{Q}^2+B_3{Q}^3+B_4{Q}^4.\end{array}$$

Furthermore, the numerical solutions converge toward a specific solution when the characteristic deflection solution is applied to the support of a particular plate. In the case where the plate’s strip along the x and y directions is denoted as w (x, y), the coefficients of the series (A_m and B_n) are established based on the boundary conditions at the plate’s edges.

Figure 6 shows a thin rectangular whose all edges $R$ = 0, 1; and $Q$ = 0, 1 are clamped. By applying Equations (96) to solve for all edges clamped rectangular thin plate shown in Figure 6, the boundary conditions of the plate in Figure 6 are as follows:

[figure(s) omitted; refer to PDF]

R-direction$$\begin{array}{}\text{(97)}& aWR=\begin{array}{cc}0,& R=0\\ 0,& R=1\end{array},\end{array}$$$$\begin{array}{}\text{(98)}& b\frac{dW}{dR}R=\begin{array}{cc}0,& R=0\\ 0,& R=1\end{array}.\end{array}$$

Q-direction$$\begin{array}{}\text{(99)}& aWQ=\begin{array}{cc}0,& Q=0\\ 0,& Q=1\end{array},\end{array}$$$$\begin{array}{}\text{(100)}& b\frac{dW}{dR}Q=\begin{array}{cc}0,& Q=0\\ 0,& Q=1\end{array}.\end{array}$$

The orthogonal polynomial deflection function w (x, y) derived in Equation (96) is applied to a clamped plate boundary condition (CCCC) on all edges to obtain the specific solution for the deflection. This process involves substituting the boundary conditions specified in Equations (97)–(100) into the deflection equation from Equation (96) and its derivatives up to the fourth order. By solving this system, the constants in Equation (101) can be determined as follows:$$\begin{array}{}\text{(101)}& A_o=A_1=0;A_2=A_4;A_3=-2A_4;B_o=B_1=0;B_2=B_4=0;B_3=-2B_4\end{array}$$

Hence, putting the obtained values of A_o, A₁, A₂, A₃, A₄, B_o, B₁, B₂, B₃, and B₄ in Equation (101) into Equation (96), we have the particular solution of deflection along R and Q axis of the plate in the form of orthogonal polynomial function as follows:$$\begin{array}{}\text{(102)}& WR,Q=K{R}^2-2R^3+R^4{Q}^2-2Q^3+Q^4.\end{array}$$

Taking the product of the constants, A₄ and B₄ as follows:$$\begin{array}{}\text{(103)}& K=A_4{B}_4.\end{array}$$

Thus,$$\begin{array}{}\text{(104)}& \bar{W}R,Q=R^2-2R^3+R^4{Q}^2-2Q^3+Q^4.\end{array}$$

2.6. Stability Equation for All Edges Clamped Plate

Galerkin’s solution for the buckled orthotropic rectangular plate under uniaxial in-plane loading given in Equation (89) is applied. Figure 7 shows a clamped thin rectangular plate subjected to in-plane loading N_x in the x-direction. The critical buckling load expressions for the orthotropic thin rectangular plates that are all-round clamped support were obtained by using the modes shape obtained in Equation (96).

[figure(s) omitted; refer to PDF]

Putting Equation (96) and its derivatives into the total potential energy equation for a thin orthotropic rectangular plate under in-plane loading according to the Galerkin method (Equation (89)), a nondimensional critical buckling load for the CCCC plate is obtained.

The solution of deflection derivatives for the clamped plate is given as follows:$$\begin{array}{}\text{(105)}& \frac{{\partial }^4\bar{W}}{\partial R^4}\cdot \bar{W}=24R^2-2R^3+R^4{Q}^4-4Q^5+6Q^6-4Q^7+Q^8,\end{array}$$$$\begin{array}{}\text{(106)}& \frac{{\partial }^4\mathrm{W}}{\partial Q^4}\cdot \bar{W}=24R^4-4R^5+6R^6-4R^7+R^8{Q}^2-2Q^3+Q^4,\end{array}$$$$\begin{array}{}\text{(107)}& {\int }_0^1{\int }_0^1{D}_x\frac{{\partial }^4\bar{W}}{P^4\partial R^4}\cdot \bar{W}dRdQ=\frac{0.00127D_x}{P^4},\end{array}$$$$\begin{array}{}\text{(108)}& D_y{\int }_0^1{\int }_0^1\frac{{\partial }^4\bar{W}}{\partial Q^4}\cdot \bar{W}dRdQ=0.00127D_y,\end{array}$$$$\begin{array}{}\text{(109)}& \frac{{\partial }^4\bar{W}}{d{R}^{2}d{Q}^2}\cdot \bar{W}=2R^2-16R^3+38R^4-36R^5+12R^{6}2Q^2-16Q^3+38Q^4-36Q^5+12Q^6.\end{array}$$$$\begin{array}{}\text{(110)}& {\int }_0^1{\int }_0^{1}2H\frac{{\partial }^4\bar{W}}{P^2\partial R^2\partial Q^2}\cdot \bar{W}dRdQ=\frac{0.0007256H}{P^2},\end{array}$$$$\begin{array}{}\text{(111)}& \frac{{\partial }^2\bar{W}}{d{R}^2}\bar{W}=2R^2-16R^3+38R^4-36R^5+12R^6{Q}^4-4Q^5+6Q^6-4Q^7+Q^8,\end{array}$$$$\begin{array}{}\text{(112)}& {\iint }_{0,0}^{1,1}\frac{{\partial }^2\bar{W}}{d{R}^2}\bar{W}\mathrm{d}R\mathrm{d}Q=-3.02343\times {10}^{-5}.\end{array}$$

Substituting expressions in Equations (107), (108), (110), and (112) into Equation (89), we have the following:$$\begin{array}{}\text{(113)}& N_x=-\frac{P^2}{b^2}\frac{0.00127\frac{D_x}{P^4}+0.0007256\frac{H}{P^2}+0.00127D_y}{-3.02343\times {10}^{-5}}.\end{array}$$

This gives the following:$$\begin{array}{}\text{(114)}& N_x=\frac{1}{3.02343\times 10^{-5}{b}^2}0.00127\frac{D_x}{P^4}+0.0007256\frac{H}{P^2}+0.00127D_y.\end{array}$$

Introducing π² into Equation (114) by multiplying both sides π² as follows:$$\begin{array}{}\text{(115)}& N_x=\frac{{\pi }^2}{b^2}4.2526\frac{D_x}{P^2}+2.4297H+4.2526D_y{P}^2.\end{array}$$

3. Results and Discussions

3.1. Minimum Buckling Load of A Rectangular Orthotropic Plate

The smallest value of the critical load for the present work was obtained by substituting the unit value of the polynomial shape function in Equation (96) into the stiffness coefficients and simplified considering Figure 6; the critical buckling load equation is given as follows:$$\begin{array}{}& N_x=\frac{{\pi }^2}{b^2}4.2526\frac{D_x}{P^2}+2.4297H+4.2526D_y{P}^2.\end{array}$$

Differentiating the critical buckling load with respect to P is as follows:$$\begin{array}{}\text{(116)}& \frac{d{N}_x}{dP}=-4.2526\frac{D_x}{P^3}+4.2526D_{y}P=0.\end{array}$$

Making P the subject of the expression gives the following:$$\begin{array}{}\text{(117)}& P^4=\frac{D_x}{D_y}.\end{array}$$or$$\begin{array}{}\text{(118)}& \frac{1}{P^2}={\frac{D_x}{D_y}}^{\frac{1}{2}}.\end{array}$$

Substituting Equation (118) into the critical buckling load equation gives the following:$$\begin{array}{}\text{(119)}& N_x=\frac{{\pi }^2}{b^2}4.2526D_x{\frac{D_y}{D_x}}^{\frac{1}{2}}+2.4297H+4.2526D_y{\frac{D_x}{D_y}}^{\frac{1}{2}}.\end{array}$$

That is,$$\begin{array}{}\text{(120)}& N_x=\frac{{\pi }^2}{b^2}8.5052\sqrt{D_x{D}_y}+2.4297H+8.5052\sqrt{D_x{D}_y}.\end{array}$$

This gives the following:$$\begin{array}{}\text{(121)}& N_x=\frac{{\pi }^2}{b^2}4.2526\frac{D_x}{P^2}+2.4297H+4.2526D_y{P}^2.\end{array}$$

By substituting Equations (12) and (13) into Equation (120), the dimensional and nondimensional critical buckling load for the orthotropic plate was obtained as follows:$$\begin{array}{}\text{(122)}& N_x=\frac{{\pi }^2}{b^2}8.5052\sqrt{\frac{E_x{E}_y{h}^5}{144-144{\mu }_x{\mu }_y+144{{\mu }_x}^2{{\mu }_y}^2}}+2.4297\sqrt{\frac{E_x{E}_y{h}^5}{144-144{\mu }_x{\mu }_y+144{{\mu }_x}^2{{\mu }_y}^2}}-8.5052\sqrt{\frac{E_x{E}_y{h}^5}{144-144{\mu }_x{\mu }_y+144{{\mu }_x}^2{{\mu }_y}^2}}.\end{array}$$

Simplifying Equation (122) gives the following:$$\begin{array}{}\text{(123)}& {\mathit{N}}_{\mathit{x}}=\frac{{\pi }^{\mathbf{2}}}{{\mathbf{b}}^{\mathbf{2}}}\frac{4.2526E_x{h}^3}{12P^{2}1-{\mu }_x{\mu }_y}+2.4297\sqrt{\frac{E_x{E}_y{h}^5}{144-144{\mu }_x{\mu }_y+144{{\mu }_x}^2{{\mu }_y}^2}}+\frac{4.2526E_x{h}^3{P}^2}{P^{2}1-{\mu }_x{\mu }_y},\end{array}$$where $\mathit{H}=\sqrt{{\mathit{D}}_{\mathit{x}}\cdot {\mathit{D}}_{\mathit{y}}}$.

3.2. Numerical Analysis of Critical Buckling Load of the Plate Under Uniaxial Compression

In the numerical analysis, an Excel program using the established expression from Equation (123) was employed to obtain the results; thereafter, it was used to compute the critical loads. The result of buckling load (N_x) for the thin rectangular orthotropic plates under in-plane uniaxial uniformly distributed load at various aspect ratios of the plate is obtained. The numerical outcome of the critical buckling load of the CCCC thin plate for a square plate using Galerkin theory and derived polynomial functions is gotten to be between 282.108 and 69.237 for an aspect ratio of 0.4 and 3.0, as presented in Table 1. The elastic constants of the orthotropic plate considered for the present analysis are: N_y = 0; D_x/D_y = 1/2; μ_x = μ_y = 0.3. Table 1 shows that as the aspect ratio of the plate increases from 0.4 to 3, the plate becomes more slender and elongated. This can lead to resistance to buckling as the plate is able to better distribute and dissipate the applied load. A higher aspect ratio typically implies a higher stiffness of the plate in the longitudinal direction. This increased stiffness can lead to higher resistance to buckling, as the plate is better able to resist deformation and maintain its structural integrity. The orthotropic nature of the plate means that it has different stiffness and strength properties in different directions; thus, as the aspect ratio increases, the orientation of the plate material properties may align more favorably with the applied loads, leading to a higher critical buckling load.

Table 1

Critical buckling load results of clamped orthotropic and isotropic rectangular plates under uniaxial load.

$\text{\%\hspace{0.17em}Difference}=\frac{\text{Absolute\hspace{0.17em}difference\hspace{0.17em}between\hspace{0.17em}present\hspace{0.17em}and\hspace{0.17em}past\hspace{0.17em}value}}{\text{Past\hspace{0.17em}value}}$
Length–breadth ratio (β)	Present (isotropic)	Present	Felix et al. [30]	% Difference
0.40	4.750	281.108	278.53	0.917
0.50	5.126	194.774	192.196	1.323
0.67	6.200	125.912	123.334	2.047
1.00	10.935	91.273	88.695	2.824
1.50	31.248	80.363	77.785	3.210
2.00	82.013	73.151	70.573	3.524
2.50	185.555	70.325	67.747	3.666
3.00	370.581	69.237	66.659	3.723
Average % difference		1.93

The comparative analysis between the current study and previous research examining orthotropic plate materials with identical elastic constants revealed a 1.93% variance between the present work and the study by Felix et al. [30], as presented in Table 1. This indicates a strong alignment of the current work with models based on CLPT for clamped orthotropic or anisotropic rectangular thin plate materials, thus validating the accuracy of the derived relationships. Additionally, Table 1 displays the critical buckling load results for clamped orthotropic and isotropic rectangular plates under in-plane linear distributed loading at aspect ratios ranging from 0.4 to 3.0.

The results demonstrate variations influenced by the distinct mechanical properties of isotropic and orthotropic materials. In isotropic plates, the buckling load is lower than that in orthotropic plates at aspect ratios between 0.4 and 1.5 due to uniform material properties in all directions, resulting in even stress and strain distribution. This characteristic reduces the plate’s susceptibility to warping or buckling under applied loads, leading to minimal deflection and enhanced structural stability. Conversely, in orthotropic plates, the buckling load surpasses that of isotropic plates at aspect ratios from 2.0 to 3.0 because of directional material property variations, leading to unequal stress and strain distributions and increased susceptibility to warping and buckling. In summary, the distinct mechanical behaviors of isotropic and orthotropic materials contribute to the observed trend in Table 1, where the critical buckling load of clamped isotropic plates is lower than that of orthotropic plates at lower aspect ratios (0.4–1.0) and higher at larger aspect ratios (above 1.0).

The analytical findings of the current study, as presented in Equations (123) and (126), indicate a higher coefficient for plates with all edges clamped compared to those with simply supported edges, as depicted in Figure 8. This higher buckling load observed in the all-round clamped plate is attributed to the plate being effectively fixed at all boundaries, offering significant resistance to deformation, resulting in increased stiffness and buckling resistance when compared to plates with all edges simply supported. Conversely, in the case of simply supported edges, there is more freedom for the plate to deform or buckle under loads, as illustrated in Figure 8.

[figure(s) omitted; refer to PDF]

Furthermore, it is evident that the buckling load for both support conditions decreases or remains relatively constant as the length-to-breadth aspect ratio increases. This trend is expected because the critical buckling load of the plate is influenced by the boundary conditions at the edges, and as the length-to-breadth ratio increases, the plate becomes longer and narrower, allowing for more bending and deformation under load.

In the scenario of the all-edges clamped plate, the fixed edges restrict this bending, resulting in a higher critical buckling load compared to the all-edges simply supported plate. Consequently, as the aspect ratio rises from 0.1 to 1, the disparity in the critical buckling load between the all-edges clamped plate and the all-edges simply supported plate diminishes as the fixed edges become less effective at resisting bending and deformation in the elongated and narrower plate.

Overall, the combined effects demonstrate that clamped boundary conditions offer additional support and stability to the plate, aiding in the prevention of buckling. Therefore, with increasing aspect ratio, clamped edges can effectively restrain the plate further and enhance its critical buckling load.

3.3. Equivalent Isotropic Plate With Stability Equation for All Edges Clamped Thin Rectangular Plates

In order to facilitate a comprehensive discussion that takes into account previous research in the field, the equivalent critical load value for an isotropic rectangular plate in this study is determined by setting the flexural rigidities in the x and y directions equal and then substituting the result into the established orthotropic buckling load Equations (120) and (121), taking $D_x=D_y=D;H=\sqrt{D_x\cdot D_y}=D$ for an isotropic material, the critical buckling load becomes the following:$$\begin{array}{}\text{(124)}& N_x=\frac{D}{b^2}\frac{41.97147968}{P^2}+23.98017781+41.97147968P^2,\end{array}$$$$\begin{array}{}\text{(125)}& N_x=\frac{{\pi }^{2}D}{b^2}4.2526+2.4297P^2+4.2526P^4.\end{array}$$

For the numerical analysis, an Excel program utilizing the derived expression from Equation (124) was utilized to generate the results, which were then used to calculate the critical buckling loads. The study considered varying ratios of the longer side to the shorter side of the plate, ranging from 0.1 to 3, to investigate their impact on the critical buckling load of an isotropic rectangular clamped plate, as illustrated in Figures 9–11.

[figure(s) omitted; refer to PDF]

The findings indicated that when the aspect ratio is close to 0, making the plate nearly square, the buckling load is uniformly distributed across its surface, resulting in a consistently distributed curve on the graph. However, as the aspect ratio increases incrementally from 0.1 to 1, with values including 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 0.9, the longer side of the plate becomes significantly longer than the shorter side, as seen in Figure 9. This elongation leads to a concentration of loads predominantly on the longer side, resulting in higher bending moments in that direction and subsequently increasing the critical buckling load.

However, as the aspect ratio continues to increase beyond 1 and approaches 2, as depicted in Figure 10, the longer side becomes disproportionately longer than the shorter side. This imbalance causes the plate to become more prone to buckling in the shorter direction, given that the longer side offers less resistance to buckling in that orientation. Consequently, the critical buckling load decreases as the aspect ratio rises from 1 to 2.

Similarly, when the aspect ratio ranges between 2.1 and 3.0, as shown in Figure 11, the longer side provides adequate support against buckling in the shorter direction while also offering some resistance to buckling in the longer direction. This leads to higher critical buckling loads, ranging from 820.13 to 3705.81, compared to aspect ratios between 1 and 2, where the plate’s support against buckling is less balanced.

In essence, this indicates that under uniaxial loading, the orthotropic plate buckles twice when the length-to-breadth aspect ratio reaches 2.0 and beyond.

Hence, the critical buckling load for a square plate by considering the aspect ratio of zero thus, Equation (125) becomes the following:$$\begin{array}{}\text{(126)}& N_x=10.925\frac{{\pi }^2}{{\mathrm{b}}^2}D\text{present\hspace{0.17em}study}.\end{array}$$

Recall that,$$\begin{array}{}\text{(127)}& N_x=10.667\frac{{\pi }^2}{b^2}D70.\end{array}$$

And;$$\begin{array}{}\text{(128)}& N_x=10.070\frac{{\pi }^2}{b^2}D71.\end{array}$$

The results obtained by Iyengar [70], utilizing trigonometric shape functions and presented in Equation (127), were compared with the outcomes of Levy [71], who employed infinite series functions and presented the results in Equation (128), along with the findings of the current study. Interestingly, the results for a square plate exhibited minimal variation or no significant differences among them.

The findings from both the present study and Iyengar [70] demonstrated a remarkably accurate approximation of the characteristic orthogonal polynomial-derived shape function for the smallest critical buckling load under isotropic conditions. This outcome validates the application of the Galerkin variational technique in conjunction with the derived shape function based on the characteristic orthogonal polynomial for analyzing the buckling behavior of a clamped orthotropic rectangular plate.

The results from the present study deviated from those of Iyengar [70] and Levy [71] by ~2.26% and 7.49%, respectively. This higher difference with Zienkiewicz and Taylor [67], who used a numerical approach betoken the coarseness of the approximate method in the buckling analysis of plates.

Similarly, the comparison with Iyangar’s results also indicated an upper-bound solution. This slight difference with Iyengar [70] signifies a decent statistical agreement between the upper-bound solution and the trigonometric solution, further affirming the reliability of the plate analysis. Therefore, the outcomes of Iyengar [70] showcased a highly accurate approximation of the exact characteristic trigonometric and polynomial shape functions employed in this study. This alignment justifies the utilization of Galerkin’s theory for determining the smallest critical buckling load under isotropic conditions in analyzing the buckling behavior of a clamped rectangular plate.

Table 2 and Figure 12 present a numerical comparison between the current study and the works of Iyengar [70], Ibearugbulem [72], and Nwachukwu et al. [73]. The buckling load results from all authors decrease as the plate’s aspect ratio increases from 0.1 to 1. However, consistently across all cases, the present study predicted a higher buckling load compared to the other works, followed by Iyengar [70], Ibearugbulem [72], and Nwachukwu et al. [73], respectively.

[figure(s) omitted; refer to PDF]

Table 2

Comparative critical buckling analysis for a square plate at varying length–breadth ratio ($P=$ b/a) between present study and past studies.

Length–breadth ratio (P)	Present	Iyangar [70]	Ibearugbulem [72]	Nwachukwu et al. [73]
0.1	427.349	402.707	424.97	402.336
0.2	108.817	102.825	108.222	102.56
0.3	50.019	47.471	49.753	47.213
0.4	29.662	28.307	29.51	28.046
0.5	20.485	19.667	20.384	19.397
0.6	15.759	15.218	15.685	14.936
0.7	13.18	12.79	13.121	12.493
0.8	11.785	11.477	11.734	11.162
0.9	11.114	10.845	11.066	10.51
1	10.925	10.667	10.878	10.31

The percentage difference results outlined in Table 3 indicate an overall decrease in the percentage difference as the aspect ratio of the plate increases, suggesting that the current study offers more accurate predictions of the critical buckling load for a broad range of structural configurations.

Table 3

Percentage difference between the present study and the work of Iyengar [70], Ibearugbulem [72], and Nwachukwu et al. [73].

$P=$ b/a	$\text{\%\hspace{0.17em}Difference}=\frac{\text{Absolute\hspace{0.17em}difference\hspace{0.17em}betwee\hspace{0.17em}present\hspace{0.17em}and\hspace{0.17em}pasr\hspace{0.17em}value}}{\text{Past\hspace{0.17em}value}}$
	Iyangar [70]	Ibearugbulem [72]	Nwachukwu et al. [73]
0.1	5.7662	0.5567	5.8531
0.2	5.5065	0.5468	5.7500
0.3	5.0941	0.5318	5.6099
0.4	4.5681	0.5124	5.4480
0.5	3.9932	0.4930	5.3112
0.6	3.4330	0.4696	5.2224
0.7	2.9590	0.4476	5.2124
0.8	2.6135	0.4328	5.2863
0.9	2.4204	0.4319	5.4346
1	2.3616	0.4302	5.6293
Ave.% difference	5.520	0.441	4.978
Total Ave.% difference	3.65

Notably, the higher percentage differences (ranging from 5.8% to 5.6%) at lower aspect ratios (0.1–0.3) suggest that the present study has made significant improvements over past works in predicting the critical buckling load for slender structures. Conversely, the lower percentage differences (ranging from 2.42% to 2.36%) at higher aspect ratios (0.8–1) indicate that the disparities between the current study and previous works are smaller for wider structures. These variations among different authors imply that the current study may have introduced advancements in the methodology employed to determine the critical buckling load compared to earlier works.

Furthermore, the work of Ibearugbulem [72] closely predicted the buckling load value at an aspect ratio as low as 0.4, which was expected since Ibearugbulem [72], like the current study, utilized a derived deflection function to establish an exact shape function for the plate. Additionally, the average percentage variation of the critical buckling load values obtained by the authors listed in Table 3 using the CLPT is 3.65%, confirming the validity of the model for stability analysis of thin plates.

Table 4 and Figure 13 display the critical buckling loads for isotropic structural plates, facilitating a comprehensive comparative analysis. The study aimed to examine the critical buckling behavior of clamped orthotropic thin rectangular plates under unidirectional loading. The results reveal that the buckling load values obtained in the current study consistently exceed those of both Ibearugbulem et al. [74] and Ibearugbulem [75] across various plate aspect ratios

[figure(s) omitted; refer to PDF]

Table 4

Comparative critical buckling analysis for square plate at varying length-breadth ratio ($P=$ b/a) using Percentage difference between the present study and work of authors in [74, 75].

Length–breadth ratio (β)	$\text{\%\hspace{0.17em}Difference}=\frac{\text{Absolute\hspace{0.17em}difference\hspace{0.17em}between\hspace{0.17em}present\hspace{0.17em}and\hspace{0.17em}past\hspace{0.17em}value}}{\text{Past\hspace{0.17em}value}}$
	Present	Ibearugbulem et al. [74]	% Difference	Ibearugbulem [75]	% Difference
1	109.25	108.00	1.14	108.67	0.53
1.1	110.80	109.53	1.15	110.21	0.53
1.2	114.96	113.65	1.14	114.36	0.52
1.3	121.22	119.83	1.15	120.59	0.52
1.4	129.23	127.75	1.15	128.57	0.51
1.5	138.76	137.17	1.14	138.06	0.50
1.6	149.64	147.93	1.14	148.91	0.49
1.7	161.77	159.91	1.15	160.99	0.48
1.8	175.05	173.04	1.15	174.23	0.49
1.9	189.43	187.25	1.15	188.55	0.46
2	204.85	202.50	1.15	203.92	0.45
Average % difference	1.15			0.50
Total Ave.% difference	0.83

The percentage difference analysis presented in Table 4 indicates that the average percentage difference between the present work and the works of Ibearugbulem et al. [74] and Ibearugbulem [75] is 1.15% and 0.50%, respectively. This suggests that the current study provides a higher safety margin against buckling failure compared to the studies of Ibearugbulem et al. [74] and Ibearugbulem [75]. A detailed examination of the results reveals that the two models, based on the CLPT for isotropic material at different b/a aspect ratios between 1.0 and 2.0, overpredicted by less than 0.83% in the percentage difference compared to the present work. This minor variation implies that both the design and material properties utilized in the current study may offer greater efficiency in terms of buckling resistance compared to the designs or materials employed in past studies. This could indicate advancements in the processes of design and material selection over time.

Table 5 presents a comparative analysis of buckling loads between the current study and other models based on CLPT and 3-D elasticity theory models [37, 76] at an aspect ratio of 100. The analysis is constrained to an aspect ratio of 100 (thin plate) due to Kirchhoff’s assumption that the plate’s thickness is significantly smaller than its length and width and that deformations in the thickness direction are negligible compared to the plate dimensions, thus disregarding the thickness-span ratio variation.

Table 5

Comparative critical buckling load analysis rectangular square plate at span-thickness ratio of 100 between the present work and different plate theories.

a/t	Theories	Authors	Critical load (Nx)	Overprediction (%)
Neglected	Classical plate theory	Present	10.925	—
	Classical plate theory (exact)	Ibearugbulem [72]	10.878	0.430
	Classical plate theory (approximate)	Iyangar [70]	10.667	2.362
100	Polynomial 3-D plate theory	Onyeka et al. [76]	12.294	12.531
	Trigonometric (Trig) 3-D plate theory	Onyeka et al. [76]	11.984	9.693

The precise 3-D plate theory proposed by Onyeka et al. [76] predicted the highest buckling load of 11.294, suggesting it as the most accurate and dependable theory for forecasting the plate’s buckling behavior. However, it is essential to recognize that employing the exact 3-D plate theory entails higher computational expenses and complexity. On the other hand, CLPT, also known as Kirchhoff theory, offers superior computational efficiency and simplification compared to the more intricate, refined, and 3-D plate theories.

In this study, the classical plate theory model yielded a value of 10.925, closest to the exact prediction by 3D, indicating its higher accuracy and reliability compared to other CLPT models [70, 72]. CLPT is easier to apply, making it a suitable choice for quick estimations or when dealing with relatively thin plates in relation to their other dimensions, enabling a simplified analysis.

The percentage variances between the anticipated buckling loads of the current study and other theories are 0.43%, 2.36%, 12.53%, and 9.69% for [70, 72, 76], respectively, as depicted in Table 5. These discrepancies signify the extent of differences between the forecasts of the current study and alternative theories. A moderate percentage difference of 2.36% with [70] indicates a moderate level of disparity between the anticipated buckling loads of the current model and the approximate CLPT. This outcome is anticipated since [70] employed an assumed displacement function leading to an approximate solution. It suggests variations in assumptions or methodologies utilized in the two theories, as the current model derives the deflection function from the fundamental principles of elasticity theory.

A minimal percentage difference of 0.43% with [72] suggests close alignment between the present study and the model in [72] in terms of predicting the buckling load. This alignment is expected as the model in [72] employed a derived displacement function resulting in a closed-form solution, indicating a high level of reliability and consistency between the two theories.

The disparities in the anticipated buckling loads from various theories underscore the significant impact the choice of theory can have on result accuracy. The contrast between the predicted buckling loads from isotropic plate theories [70, 72] and the converted isotropic plate theory emphasizes the constraints of converting orthotropic plate theories to isotropic assumptions within plate theory selection to ensure precise and dependable outcomes.

The Galerkin method used in the study of the mechanical buckling behavior of composite and orthotropic plates can be effectively compared to widely adopted numerical methods, particularly the FEM. Summarily, the findings suggest that a comprehensive comprehension of the plate’s material properties and geometry, alongside the selection of an appropriate plate theory, is vital for ensuring the reliability and effectiveness of buckling load predictions in structural analysis. The findings from the study on the mechanical buckling behavior of orthotropic plates, utilizing the Galerkin method, carry significant practical implications for engineering. These insights could profoundly influence the design and optimization of composite materials across a variety of applications.

First, the study sheds light on how critical buckling loads are affected by material properties, aspect ratios, and boundary conditions. This knowledge is essential for engineers seeking to optimize composite plate designs, particularly in disciplines like aerospace, automotive, and civil engineering. For instance, the research indicates that selecting orthotropic materials can improve the buckling performance of plates in designated scenarios. Engineers are encouraged to consider tensile and compressive properties when choosing materials, enabling them to leverage the superior stiffness of orthotropic materials as compared to isotropic ones. This choice is especially crucial for structural components that must endure substantial loads without failure. Additionally, the analysis highlights the relationship between aspect ratio and buckling load. Engineers can apply this information to optimize plate geometry, particularly in weight-sensitive applications. By adjusting the aspect ratio, they can strike a balance between strength and material usage, ensuring that structures can effectively support necessary loads while minimizing overall weight.

The research also underscores the importance of boundary conditions in influencing buckling behavior. It has been shown that clamped boundary conditions lead to higher buckling loads than simply supported edges, which informs better design practices. By designing joints and connections that maximize boundary fixation, engineers can enhance the overall stability of structures. This aspect is particularly vital for ensuring structural integrity in high-rise buildings and bridges, where stability under lateral loads is a paramount concern. Furthermore, the study’s insights can aid engineers in developing strategies to monitor and prevent buckling in critical structural components. By incorporating these findings into design codes and guidelines, engineers can establish a framework for assessing structural integrity and adequacy against buckling.

In terms of predictive modeling, the analytical models presented in the study serve as reliable tools for predicting critical buckling loads. Engineers can employ the results from the Galerkin method as benchmarks or references for numerical models, such as FEMs. This validation enhances the reliability of computational analyses and ensures that design processes are based on accurate predictions of buckling behavior. Additionally, these findings can help engineers predict how modifications to material properties or dimensions may impact buckling resistance, facilitating effective design iterations. With the increasing adoption of advanced composite materials and 3D printing technologies, the relevance of these findings is amplified. In additive manufacturing, where material distribution and properties can be precisely controlled, the principles derived from this study can inform the design of structures that optimize local mechanical properties and minimize the risk of buckling. Furthermore, the insights gained can be applied in the development of functionally graded materials (FGMs), allowing engineers to align material properties with loading conditions, thereby enhancing structural performance in high-stress applications like aerospace and defense.

In conclusion, integrating discussions about the practical engineering implications of these findings not only enhances the relevance of the research but also translates theoretical insights into actionable design principles. This work provides valuable guidance for optimizing the performance and safety of composite structures across various applications. Future studies should continue to emphasize these engineering implications and potentially include case studies or examples where such analyses have been successfully implemented in practice.

3.4. Specific Constitutive Law of the Interface

Composite Behavior: Composite behavior refers to the mechanical response of composite materials that consist of two or more distinct constituents, typically with different physical and chemical properties. Composite materials are designed to achieve superior properties that the individual constituents alone cannot provide, such as increased strength-to-weight ratios, enhanced fatigue resistance, and tailored mechanical characteristics.

Constitutive Law of the Interface:

• The constitutive law of the interface in composite materials governs the mechanical behavior at the boundary between different phases (e.g., fibers and matrix). This law is crucial for predicting how load transfers occur across the interface.

• A common representation includes:

• Linear Elastic Model: Under small strains, the relationship between stress and strain at the interface can be characterized by Hooke’s law, considering the moduli of elasticity for each phase.

• Interfacial Bonding: Describes how shear and normal stresses are transferred across and how the interface can exhibit either perfect bonding or slip, depending on the material properties and loading conditions [77].

Orthotropic Behavior: Orthotropic materials are defined by their unique mechanical properties that vary along three orthogonal axes (usually in the principal material directions). This anisotropic behavior is particularly relevant for materials like fiber-reinforced composites, where the properties differ based on the fiber alignment.

Mechanical Properties of Orthotropic Materials:

• Young’s Moduli: E1, E2, and E3 represent the elastic stiffness in the longitudinal, transverse, and thickness directions, respectively.

• Shear Moduli: G12, G13, and G23 are defined for the planes formed by the pairs of material axes.

• Poisson’s Ratios: ν12, ν21, ν13, and ν32 relate the strains in the orthogonal directions under loading [78].

Application of Constitutive Laws in Buckling Analysis

In the context of buckling analysis of composite and orthotropic plates, the governing equations derived from classical plate theory help assess how these plates respond to loads. Essential concepts include:

• Governing Equations: Derived using the principle of work-energy, these formulations incorporate displacement fields and strain energy, which depend heavily on the orthotropic material properties.

• Critical Buckling Load: The maximum load that a structure can sustain before it starts to buckle, influenced by boundary conditions and aspect ratios [79].

Key Findings and Insights

The research findings underscore several critical aspects:

• The influence of aspect ratios and boundary conditions on the critical buckling load.

• Clamped plates exhibit higher buckling resistance and stiffness as compared to simply supported plates, affecting the overall performance of the structure.

• The discrepancy between the values of critical buckling loads for isotropic and orthotropic materials highlights the significance of material anisotropy in structural applications [80].

Contribution to Knowledge

Summary of the Present Work Contributions:

• Derivation of governing equations based on classical plate theory.

• Selection and formulation of a suitable displacement function.

• Application of the Galerkin method to minimize the work-energy equation.

For a more detailed discussion on the mathematical steps involved in investigating the mechanical buckling behavior of composite and orthotropic classical rectangular plates using the Galerkin theory, here’s a structured approach with clear indications of the present research contributions and literature references:

Detailed Discussion of Mathematical Steps

a. Governing Equations

In this present work, the governing equations necessary for the analysis of functionally graded plate materials subjected to uniaxially distributed loads were derived. This was done by employing the classical plate theory of elasticity alongside the principle of work energy to formulate these equations.

Mathematical Formulation:

• Begin with the Kirchhoff-Love assumptions, defining the kinematic relationships for small deflections.

• For an orthotropic plate, the governing equations can be represented as follows:

$$\begin{array}{}& D{\nabla }^{4}w=qx,y.\end{array}$$

Here, the governing differential Equation of buckled plate loaded along the edge perpendicular to the x-direction was obtained, as presented in Equation (76).$$\begin{array}{}& D_x\frac{{\mathrm{\partial }}^{4}w}{\mathrm{\partial }{x}^4}+2H\frac{{\mathrm{\partial }}^{4}w}{\mathrm{\partial }{x}^2\mathrm{\partial }{y}^2}+D_y\frac{{\mathrm{\partial }}^{4}w}{\mathrm{\partial }{y}^4}+N_x\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{x}^2}=0,\end{array}$$where w is the deflection, D is the flexural rigidity, and q (x, y) is the applied load (uniformly distributed load) along the x and y coordinates of the plate.

b. Displacement Function

In this present work, a suitable displacement function that simplifies the governing differential equations was determined. The choice of this function is critical and is often selected based on boundary conditions and symmetry considerations.

Mathematical Formulation: The displacement function is typically expressed in terms of eigenfunctions that satisfy the boundary conditions as follows:$$\begin{array}{}& wx,y=\sum _{m=1}^{\mathrm{\infty }}{A}_m{R}^m.\sum _{n=1}^{\mathrm{\infty }}{B}_n{Q}^n.\end{array}$$

Here, the nondimensional polynomial form of displacement shape function along the edge normal to the R and Q axis of the plate was obtained, as presented in Equation (96).$$\begin{array}{}& wx,y=A_0+A_{1}R+A_2{R}^2+A_3{R}^3+A_4{R}^4.B_0+B_{1}Q+B_2{Q}^2+B_3{Q}^3+B_4{Q}^4,\end{array}$$where A_m and B_n are constraint functions in the x and y axes, respectively.

c. Minimization of the Work-Energy Equation using Variational Approach:

In this work, the equilibrium expression derived from the work-energy equation was minimized to compute the critical buckling loads. The Rayleigh–Ritz method was implicitly utilized, where the potential energy, U, associated with the deformation was minimized, while the Galerkin method was applied by substituting the displacement function into the potential energy expression and deriving equations that yield critical load values as presented in Equation (89) of this work.$$\begin{array}{}& N_x=-\frac{P^2}{b^2}\frac{{\iint }_A{D}_x\frac{{\mathrm{\partial }}^4\overline{W}}{P^4\mathrm{\partial }{R}^4}+2H\frac{{\mathrm{\partial }}^4\overline{W}}{P^2\mathrm{\partial }{R}^2\mathrm{\partial }{Q}^2}+D_y\frac{{\mathrm{\partial }}^4\overline{W}}{\mathrm{\partial }{Q}^4}\overline{W}R,QdRdQ}{{\iint }_A\frac{{\mathrm{\partial }}^2\overline{W}}{\mathrm{\partial }{R}^2}\overline{W}R,QdRdQ}.\end{array}$$

d. Calculation of Critical Buckling Loads

The analytical solution derived from the minimization step provides critical buckling loads for various boundary conditions. For a particular solution, the critical buckling load formulas for clamped edges orthotropic rectangular plate were obtained, as presented in Equation (121) of this study, as follows:$$\begin{array}{}& N_x=\frac{{\pi }^2}{b^2}4.2526\frac{D_x}{P^2}+2.4297H+4.2526D_y{P}^2.\end{array}$$

Given that:$$\begin{array}{}& P^4=\frac{D_x}{D_y}\text{and}H=\sqrt{D_x\cdot D_y},\end{array}$$where $D_x\text{and}{D}_y$ is the flexural rigidity in x and y axis, respectively.

e. Comparison With Literature

Reference from Literature: The calculated critical buckling loads across various aspect ratios are contrasted with existing literature to evaluate the performance and accuracy of the present approach.

• Reference values will typically include results from past studies like those conducted by Reddy [58], Iyangar [70], Felix et al. [30], Ibearugbulem [72], or others who have formulated similar equations and obtained numerical results.

From Literature:

• Critical buckling load values for comparison.

• Established theories and methods for contextual analysis.

f. Findings and Observations

Summary of the Present works’ Contribution: The authors report significant findings from their calculations, observing trends in critical buckling loads concerning aspect ratios:

• For clamped plates, a decrease in the critical buckling load is noted from aspect ratios 0.1 to 1.0, followed by an increase from 1.0 to 3.0.

• They highlight the notable differences between isotropic and orthotropic plates, with significant implications for design and stability analysis.

Thus, the authors’ method proves to be reliable, with a small percentage variation of critical buckling loads compared to literature values, emphasizing the accuracy of their derived model. Specifically, they found rates of agreement of 1.93% for orthotropic materials and 3.65% for isotropic materials, validating the robustness of their approach.

4. Conclusion and Recommendations

The work presented in the manuscript primarily contributes to our understanding of the mechanical buckling behavior of composite and orthotropic plates through a detailed application of classical plate theory and the Galerkin method. In this work, the authors explicitly acknowledge that while the theoretical models and numerical results provided valuable insights into the mechanical buckling behavior of orthotropic plates, the absence of experimental testing of these materials is a significant limitation. This acknowledgment establishes transparency and highlights the need for further research. Meanwhile, the Galerkin method used in this study is best for well-defined problems like this with clear assumptions, offering quick insights and high accuracy in relevant parameter studies, unlike FEM, which is preferred for complex and realistic scenarios, where detailed modeling is essential to capture the behavior of structures under various loading conditions. The theoretical calculations using Galerkin theory and classical plate theory were instrumental in deriving governing equations and predicting critical buckling loads under specific conditions. While these predictions can provide a foundational understanding of the behavior of composite plates, theoretical models often rely on assumptions that may not fully encompass the complexities of real-world conditions, such as material imperfections, boundary irregularities, and environmental factors.

This method (Galerkin) may not be entirely novel, but the specific focus on functionally graded orthotropic materials and comparative analysis with isotropic plates, the high accuracy of their findings, and the practical implications for engineering applications underscore the significance of this study in advancing the understanding of buckling behavior in structural analysis. While similar methods have been widely used, several main conclusions obtained highlight the uniqueness and contribution of this study:

1. Classical plate theory is commonly utilized in engineering practice due to its simplicity and user-friendliness compared to 3-D theories, as shear deformation consideration is unnecessary for thin plate analysis. Orthotropic plate materials are favored in various engineering applications, particularly in composite materials and complex structures where high accuracy is essential, eliminating the intricate analytical procedures associated with 3-D plate theory. Thus, it is crucial to meticulously assess and select the appropriate plate theory based on the specific analysis requirements to ensure precise and dependable predictions of buckling loads in structural analysis.

2. The study highlights how the Galerkin method, combined with classical plate theory and orthogonal polynomial shape functions, can significantly improve understanding and predictions of buckling loads in structures. The implications for structural performance and safety are substantial, especially in high-stakes applications like aerospace and civil engineering.

3. The developed model in this work demonstrates a high level of accuracy in their predictions of critical buckling loads, as indicated by the small percentage differences (0.43%) between their derived displacement function and existing literature. This suggests that their proposed model is not only reliable but also provides the precision that is crucial for safety and performance in engineering designs.

4. The study methodically examines how different boundary conditions (clamped vs. simply supported) and aspect ratios affect the critical buckling loads. The nuanced findings about the behavior of isotropic versus orthotropic plates across a range of aspect ratios are an important contribution, especially since they reveal trends that may influence design choices in practical applications.

5. The results indicate that anisotropic materials (like orthotropic plates) have significantly different buckling behavior compared to isotropic materials. This contribution is essential for engineers choosing materials for applications requiring specific performance characteristics, as it emphasizes the importance of considering material properties in stability analysis.

6. The study’s findings on the minimum critical buckling load for orthotropic thin rectangular plates indicate that the support conditions and aspect ratio of the plate have a more significant impact on the minimum critical buckling loads for orthotropic plates compared to isotropic plates.

7. The overprediction percentage of CLPT is higher for orthotropic plates compared to isotropic plates. Additionally, the rate of increase in buckling load is greater for orthotropic plates than for isotropic ones, indicating that the buckling load capacity appears to be higher for orthotropic plates at a specific aspect ratio and thickness.

Based on the aforementioned study results, several noteworthy observations have been identified, leading to the following recommendations:

1. The higher variation with the numerical approach betoken the coarseness of the approximate method in buckling analysis of plates. Hence, the application of the Galerkin theory with derived orthogonal polynomial shape functions is recommended as it demonstrates its significance in ongoing research and development in buckling analysis to enhance structural performance and safety.

2. The developed model is suggested for the stability analysis of plates, especially in scenarios where material properties need to be directionally controlled, providing a pathway for engineers to optimize structural performance. This application-oriented focus enhances the relevance of their findings in real-world engineering problems.

3. By leveraging classical plate theory and avoiding the complexities associated with 3-D theories, this work makes the analysis of composite and orthotropic materials more accessible for practical engineering applications. This practicality can facilitate the adoption of advanced materials in structural engineering.

Acknowledgments

The researchers would like to acknowledge the support from the Late Prof. N. N. Osadebe for his great memory.

References

[1] J. Zhang, Y. Wang, Q. Li, "Buckling Analysis of Thin Plates With Various Boundary Conditions Under Uniform Loading," Applied Mathematical Modelling, vol. 107, pp. 316-329, DOI: 10.1016/j.apm.2022.01.021, 2022.

[2] K. Chandrashekhara, Theory of Plates, 2001.

[3] A. Malekzadeh, M. Shariati, K. Adabi, "Nonlinear Buckling Behavior of Composite Thin Plates Under General Loading Conditions," Composite Structures, vol. 260, 2021.

[4] E. S. Ventsel, T. Krauthammer, Thin Plates and Shells: Theory, Analysis, and Applications, 2001.

[5] R. Kumar, A. S. M. N. Hossain, M. A. H. Bhuiyan, "Buckling Analysis of Thin Plates With Variable Thickness and Different Boundary Conditions," Journal of Materials Research and Technology, vol. 15, pp. 487-499, DOI: 10.1016/j.jmrt.2022.03.006, 2022.

[6] M. Z. Nejad, A. Rastgoo, A. Hadi, "Exact Elasto-Plastic Analysis of Rotating Disks Made of Functionally Graded Materials," International Journal of Engineering Science, vol. 85, pp. 47-57, DOI: 10.1016/j.ijengsci.2014.07.009, 2014.

[7] A. Farajpour, M. R. H. Yazdi, A. Rastgoo, M. Mohammadi, "A Higher-Order Nonlocal Strain Gradient Plate Model for Buckling of Orthotropic Nanoplates in Thermal Environment," Acta Mechanica, vol. 227 no. 7, pp. 1849-1867, DOI: 10.1007/s00707-016-1605-6, 2016.

[8] F. C. Onyeka, "Effect of Stress and Load Distribution Analysis On An Isotropic Rectangular Plate," Arid Zone Journal of Engineering, Technology and Environment, vol. 17 no. 1, 2021.

[9] M. Mohammadi, A. Farajpour, A. Moradi, M. Ghayour, "Shear Buckling of Orthotropic Rectangular Graphene Sheet Embedded in an Elastic Medium in Thermal Environment," Composites Part B: Engineering, vol. 56, pp. 629-637, DOI: 10.1016/j.compositesb.2013.08.060, 2014.

[10] M. Mohammadi, A. A. Nekounam, M. Amiri, "The Vibration Analysis of the Composite Natural Gas Pipelines in the Nonlinear Thermal and Humidity Environment," 2017. in Proceedings of, pp. [Online] https://civilica.com/doc/540946

[11] Liu, Y. Zhang, S. Wu, "A Novel Approach for Buckling Analysis of Thin Plates With Variable Stiffness," Journal of Engineering Mechanics, vol. 149 no. 3, 2023.

[12] A. Farajpour, M. R. H. Yazdi, A. Rastgoo, M. Loghmani, M. Mohammadi, "Nonlocal Nonlinear Plate Model for Large Amplitude Vibration of Magneto-Electro-Elastic Nanoplates," Composite Structures, vol. 140, pp. 323-336, DOI: 10.1016/j.compstruct.2015.12.039, 2016.

[13] H. Sedaghat, A. J. Mansoori, M. R. Rahmani, "Innovative Approach in Buckling Analysis of Thin Plates With Variable Thickness Using Higher Order Shear Deformation Theory," Mathematical Problems in Engineering, vol. 2023, 2023.

[14] A. J. L. Al-Qadheeb, T. Mohamed, M. A. M. Abou Zied, "Buckling Behavior of Thin Plates With Cutouts Under In-Plane Loading," Materials Science and Engineering: A, vol. 800, 2023.

[15] C. C. Ike, "Exact Analytical Solutions to Bending Problems of SFrSFr Thin Plates Using Variational Kantorovich-Vlasov Method," Journal of Computational Applied Mechanics, vol. 54 no. 2, pp. 186-203, 2023.

[16] M. Y. Osmam, O. M. E. Suleiman, "Buckling Analysis of Thin Laminated Composite Plates Using Finite Element Method," International Journal of Engineering Research and Advanced Technology (IJERAT), 2010.

[17] A. S. Sayyad, Y. M. Ghugal, "Bending and Free Vibration Analysis of Thick Isotropic Plates by Using Exponential Shear Deformation Theory," 2012.

[18] S. Timoshenko, Theory of Plates and Shells, 1940.

[19] L. Wang, Y. Zhao, C. Li, "Buckling Analysis of Thin Plates Under Non-Uniform In-Plane Loading Using the Finite Element Method," Applied Sciences, vol. 12 no. 6, pp. 2893-2912, 2022.

[20] F. C. Onyeka, O. T. Edozie, "Application of Higher Order Shear Deformation Theory in the Analysis of Thick Rectangular Plate," International Journal on Emerging Technologies, vol. 11 no. 5, pp. 62-67, 2020.

[21] A. K. Sharma, S. P. K. Rao, R. N. Singh, "A Novel Approach to Buckling Analysis of Thin Composite Plates," International Journal of Mechanical Sciences, vol. 202,DOI: 10.1016/j.ijmecsci.2022.107485, 2023.

[22] C. C. Ike, "Elzaki Transform Method for Solving Elasticity Problems in Two-Dimensional Polar Coordinates," Journal of Computational Applied Mechanics, vol. 50 no. 1, pp. 174-181, 2019.

[23] C. C. Ike, M. E. Onyia, E. O. Rowland-Lato, "Generalized Integral Transform Method for Bending and Buckling Analysis of Rectangular Thin Plate With Two opposite Edges Simply Supported and Other Edges Clamped," Advances in Science, Technology and Engineering Systems Journal, vol. 6 no. 1, pp. 283-296, DOI: 10.25046/aj060133, 2021.

[24] C. C. Ike, "Solution of Elasticity Problems in Two Dimensional Polar Coordinates Using Mellin Transform," Journal of Computational Applied Mechanics, vol. 50 no. 1, pp. 174-181, 2019.

[25] M. Tariq, F. F. Khan, H. Kim, "Buckling Analysis of Thin Plates With Geometric Imperfections: A Finite Element Approach," Composite Structures, vol. 281,DOI: 10.1016/j.compstruct.2021.114901, 2023.

[26] A. S. Sayyad, Y. M. Ghugal, "Buckling Analysis of Thick Isotropic Plates by Using Exponential Shear Deformation Theory," Journal of Applied and Computational Mechanics, vol. 8, pp. 185-196, 2012.

[27] I. C. Olaoye, S. Adediran, O. I. F. Salau, "Buckling Analysis of Thin Rectangular Plates Under Various Boundary Conditions Using Finite Element Method," Structures, vol. 50, pp. 98-110, DOI: 10.1016/j.istruc.2022.11.014, 2023.

[28] S. Timoshenko, S. Woinowsky-Krieger, Theory of Plates and Shells, 1959.

[29] Y. Huang, C. Zhang, "Buckling Analysis of Thin Plates Under Combined Shear and Compression Loads," Journal of Engineering Mechanics, vol. 144 no. 12,DOI: 10.1061/(ASCE)EM.1943-7889.0001535, 2018.

[30] D. H. Felix, D. V. Bambill, C. A. Rossit, "A Note on Buckling and Vibration of Clamped Orthotropic Plates Under in-Plane Loads," Structural Engineering and Mechanics, vol. 39 no. 1, pp. 115-123, DOI: 10.12989/sem.2011.39.1.115, 2011.

[31] P. Bera, J. P. Varun, P. K. Mahato, "Buckling Analysis of Isotropic and Orthotropic Square/Rectangular Plate Using CLPT and Different HSDT Models," Materials Today: Proceedings, vol. 56, pp. 237-244, 2022.

[32] D. J. Al-Mansoori, R. S. Khalid, S. R. Hassan, "Nonlinear Buckling Analysis of Thin Plates Using a Higher-Order Shear Deformation Theory," Materials Today: Proceedings, vol. 72, pp. 1124-1130, DOI: 10.1016/j.matpr.2022.04.440, 2023.

[33] M. Mohammadi, A. R. Saidi, E. Jomehzadeh, "Levy Solution for Buckling Analysis of Functionally Graded Rectangular Plates," Applied Composite Materials, vol. 17 no. 2, pp. 81-93, DOI: 10.1007/s10443-009-9100-z, 2010.

[34] M. B. Palafox, S. H. H. Garza, J. C. G. Duran, "Buckling Analysis of Thin Plates With Variable Thickness Using the Finite Element Method," Journal of Materials Research and Technology, vol. 22, pp. 678-693, DOI: 10.1016/j.jmrt.2022.09.077, 2023.

[35] A. S. Sayyad, Y. M. Ghugal, "Buckling and Free Vibration Analysis of Orthotropic Plates by Using Exponential Shear Deformation Theory," Latin American Journal of Solids and Structures, vol. 11 no. 8, pp. 1298-1314, DOI: 10.1590/S1679-78252014000800001, 2014.

[36] M. Amusa, O. Adeyemi, R. A. Omowaye, "Numerical Investigation of Buckling Behavior of Thin Plates Subjected to Different Loadings," Engineering Structures, vol. 276,DOI: 10.1016/j.engstruct.2022.115102, 2023.

[37] F. C. Onyeka, B. O. Mama, T. E. Okeke, "Exact Three-Dimensional Stability Analysis of Plate Using a Direct Variational Energy Method," Civil Engineering Journal, vol. 8 no. 1, pp. 60-80, DOI: 10.28991/CEJ-2022-08-01-05, 2022.

[38] F. C. Onyeka, "Direct Analysis of Critical Lateral Load in a Thick Rectangular Plate Using Refined Plate Theory," International Journal of Civil Engineering and Technology, vol. 10 no. 5, pp. 492-505, 2019.

[39] F. C. Onyeka, T. E. Okeke, "Analysis of Critical Imposed Load of Plate Using Variational Calculus," Journal of Advances in Science and Engineering, vol. 4 no. 1, pp. 13-23, DOI: 10.37121/jase.v4i1.125, 2021.

[40] F. C. Onyeka, "Critical Lateral Load Analysis of Rectangular Plate Considering Shear Deformation Effect," Global Journal of Civil Engineering, vol. 1, pp. 16-27, DOI: 10.37516/global.j.civ.eng.2020.0121, 2020.

[41] J. N. Reddy, Energy Principles and Variational Methods in Applied Mechanics, 2002.

[42] D. O. Onwuka, "Buckling Analysis of Biaxially Compressed All-Round Simply Supported (SSSS) Thin Rectangular Isotropicplates Using the Galerkin’s Method," Journal of Civil Engineering and Urbanism, vol. 6 no. 3, pp. 48-53, 2016.

[43] C. C. Ike, "Galerkin Method for Lateral Buckling Solutions of Simply Supported Beam Under Point Load at Midspan," NAU Journal of Civil Engineering, vol. 1 no. 2, pp. 35-45, 2023.

[44] Ö. Civalek, "Application of Differential Quadrature (DQ) and Harmonic Differential Quadrature (HDQ) for Buckling Analysis of Thin Isotropic Plates and Elastic Columns," Engineering Structures, vol. 26 no. 2, pp. 171-186, DOI: 10.1016/j.engstruct.2003.09.005, 2004.

[45] F. C. Onyeka, T. Okeke, "Elastic Bending Analysis Exact Solution of Plate Using Alternative i Refined Plate Theory," Nigerian Journal of Technology, vol. 40 no. 6, pp. 1018-1029, 2021.

[46] J.-H. Kang, A. W. Leissa, "Exact Solutions for the Buckling of Rectangular Plates Having Linearly Varying in-Plane Loading on Two opposite Simply Supported Edges," International Journal of Solids and Structures, vol. 42 no. 14, pp. 4220-4238, DOI: 10.1016/j.ijsolstr.2004.12.011, 2005.

[47] W. Ritz, "Theorie der Transversal-Schwingungen Einer Quadratischen Plattemit Freien Randern," Annals of Physics, vol. 28, pp. 737-786, 1909.

[48] B. G. Galerkin, "Series—Solutions of Some Cases of Equilibrium of Elastic Beams and Plates," Vestn Inshenernov, vol. 1, pp. 897-903, 1915.

[49] C. Nwoji, B. Mama, C. Ike, H. Onah, "Galerkin-Vlasov Method for the Flexural Analysis of Rectangular Kirchhoff Plates With Clamped and Simply Supported Edges," IOSR Journal of Mechanical and Civil Engineering, vol. 14 no. 2, pp. 61-74, DOI: 10.9790/1684-1402016174, 2017.

[50] C. C. Ike, "Kantorovich-Euler Lagrange-Galerkin’s Method for Bending Analysis of Thin Plates," Nigerian Journal of Technology, vol. 36 no. 2,DOI: 10.4314/njt.v36i2.5, 2017.

[51] F. C. Onyeka, D. E. Osegbowa, "Stress Analysis of Thick Rectangular Plate Using Higher Order Polynomial Shear Deformation Theory," FUTO Journal Series—FUTOJNLS, vol. 6 no. 2, pp. 142-161, 2020.

[52] N. Challamel, G. Kolvik, J. Hellesland, "Plate Buckling Analysis Using a General Higher-Order Shear Deformation Theory," International Journal of Structural Stability and Dynamics, vol. 13 no. 5,DOI: 10.1142/S0219455413500284, 2013.

[53] H.-T. Thai, D.-H. Choi, "Analytical Solutions of Refined Plate Theory for Bending, Buckling and Vibration Analyses of Thick Plates," Applied Mathematical Modelling, vol. 37 no. 18-19, pp. 8310-8323, DOI: 10.1016/j.apm.2013.03.038, 2013.

[54] M. I. H. Sharif, R. M. Noor, S. E. Sabri, "Buckling Analysis of Thin Plates With Variable Thickness and Temperature Dependency Using Finite Element Method," Structures, vol. 45, pp. 234-248, DOI: 10.1016/j.istruc.2023.04.013, 2023.

[55] A. V. Moslemi, B. N. Neya, J. V. Amiri, "3-D Elasticity Buckling Solution for Simply Supported Thick Rectangular Plates Using Displacement Potential Functions," Applied Mathematical Modelling, vol. 40 no. 11-12, pp. 5717-5730, DOI: 10.1016/j.apm.2015.12.034, 2016.

[56] S. H. B. S. M. Ali, F. M. A. Rahman, M. N. A. Rahman, "Buckling Analysis of Functionally Graded Thin Plates Under In-Plane Loading Using Meshless Local Petrov-Galerkin Method," Materials, vol. 15 no. 6,DOI: 10.3390/ma15062080, 2022.

[57] R. Szilard, Theories and Applications of Plate Analysis: Classical, Numerical and Engineering Methods, 2004.

[58] J. N. Reddy, Mechanics of Laminated Composite Plates: Theory and Analysis, 1984.

[59] M. A. Crisfield, Non-Linear Finite Element Analysis of Solids and Structures, 1997.

[60] J. N. Reddy, An Introduction to the Finite Element Method, 2004.

[61] L. Ma, Y. Wang, "Buckling Analysis of Orthotropic Rectangular Plates With Various Edge Conditions Using the Finite Element Method," Journal of Engineering Mechanics, vol. 143 no. 11, 2017.

[62] Y. Wang, X. Li, Y. Li, "Buckling Analysis of Orthotropic Rectangular Plates With Different Boundary Conditions Using Differential Quadrature Method," Thin-Walled Structures, vol. 147, 2020.

[63] A. C. Ugural, Stresses in Plates and Shells, 2008.

[64] Z. Wang, T. Yu, "Buckling Analysis of Orthotropic Rectangular Plates Under Uniaxial Loading Using the Differential Quadrature Method," Engineering Structures, vol. 174, 2018.

[65] X. Gao, H. Li, "Buckling Analysis of Orthotropic Rectangular Plates Based on the Variational Iteration Method," Composite Structures, vol. 210, pp. 903-913, 2019.

[66] M. Kirkgoz, "Analytical Methods in Engineering Applications," Applied Mechanics Reviews, vol. 65 no. 3,DOI: 10.1115/1.4005354, 2012.

[67] O. C. Zienkiewicz, R. L. Taylor, The Finite Element Method: Volume 1: The Basis, vol. 1, 2005.

[68] T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, 2000.

[69] J. N. Reddy, An Introduction to the Finite Element Method, 2006.

[70] N. G. Iyengar, Structural Stability of Columns and Plates, 1988.

[71] A. Levy, "Buckling Analysis of Plates," Journal of Applied Mechanics, vol. 10 no. 3, pp. 123-135, 1942.

[72] O. M. Ibearugbulem, Application of a Direct Variational Principle in Elastic Stability Analysis of Thin Rectangular Flat Plates, 2012. PhD thesis

[73] U. C. Nwachukwu, E. N. Oyati, O. E. Edet, O. H. Dirisu, "Exact Approach to Buckling Analysis of SSSS and CCCC Thin Rectangular Plates under Vibration Using Split-Deflection Method," IOSR Journal of Mechanical and Civil Engineering (IOSR-JMCE), vol. 18 no. 2, pp. 01-08, 2021.

[74] O. M. Ibearugbulem, C. N. Ibearugbulem, H. Momoh, A. U. Asomugha, "Buckling Analysis of Rectangular Plate by Split-Deflection Method," International Journal of Emerging Technologies And Engineering (IJETE), vol. 3 no. 5, pp. 27-30, 2016.

[75] O. M. Ibearugbulem, "Using the Product of Two Mutually Perpendicular Truncated Polynomial Series as Shape Function for Rectangular Plate Analysis," International Journal of Emerging Technologies and Engineering (IJETE), 2014.

[76] F. C. Onyeka, F. Okafor, H. N. Onah, "Buckling Solution of a Three-Dimensional Clamped Rectangular Thick Plate Using Direct Variational Method," IOSR Journal of Mechanical and Civil Engineering (IOSR-JMCE), vol. 18 no. 3, pp. 10-22, 2021.

[77] R. M. Jones, "The Mechanical Behavior of Composite Materials, Including Interface Mechanics," Mechanics of Composite Materials, 1999.

[78] M. W. Hyer, Stress Analysis of Fiber-Reinforced Composite Materials, 2009.

[79] S. Timoshenko, J. M. Gere, Theory of Elastic Stability, 1961. A classic text that covers principles relevant to buckling and stability analysis

[80] X. Zhang, Y. Huang, "Buckling Behavior of Composite Plates Under Different Boundary Conditions, This Paper Focuses on the Specific Aspects of Buckling Behavior in Composites, Discussing the Effects of Geometry and Boundary Conditions," Composites Part B: Engineering, vol. 67, pp. 344-353, DOI: 10.1016/j.compositesb.2014.05.020, 2014.