Nonlinear thermo-mechanical buckling and postbuckling behavior of FG-GPLRC panels with auxetic core and piezoelectric layers

Abstract

In this paper, the nonlinear large deflection and higher-order shear deformation shell theories are combined to analyze the nonlinear thermo-mechanical buckling and postbuckling behaviors of functionally graded graphene platelet reinforced composite (FG-GPLRC) complexly curved panels with porous core and piezoelectric layers. Different distributions of graphene platelet (GPL), and different types of curvatures, including parabola, sinusoidal, and cylindrical curvatures, are investigated. Four types of simply supported boundary conditions are taken into account, and the potential energy expression is formulated. The nonlinear three-parameter elastic foundation model and Ritz energy method are used to obtain the critical buckling and postbuckling responses of panels. Numerical examples present different effects of the GPL distribution and volume fraction, curvature of panels, piezoelectric layers, foundation, and uniformly distributed temperature change on the nonlinear buckling responses of the considered panels.

Introduction

Plate and panel structures typically have flat or slightly curved shapes, often rectangular or circular, and are designed to withstand and distribute loads across their surfaces. Their notable features include good load-bearing capacity and lightweight, which help reduce material usage and construction costs. Common applications of these structures include construction (walls, floors), industry (machine casings), and transportation infrastructure (bridge decks, tunnels).

Functionally Graded Materials (FGMs) are materials with gradually changing composition and structure, leading to continuous variations in physical and mechanical properties. FGMs optimize properties such as strength, heat resistance, and wear resistance, and are widely used in aerospace, energy, and mechanical applications. Linear mechanical and thermal buckling, nonlinear thermomechanical bending, and nonlinear dynamic analysis of FGM plates were investigated using the first-order shear deformation theory (FSDT) with element-free method, ¹ using the FSDT with meshless local Petrov–Galerkin, ² using the higher-order shear deformation theory (HSDT) with perturbation method. ³ By improving the smeared stiffener technique, the Galerkin method combined with FSDT ⁴ and HSDT ⁵ was used to study the nonlinear thermo-mechanical buckling behavior of FGM plates with oblique stiffeners. The thermal stresses and panel-foundation interaction were used to study the nonlinear thermal bending analysis of FGM cylindrical panels, taking into account the HSDT and perturbation method. ⁶ The physical neutral surface was considered, and HSDT was applied to the problem of nonlinear free vibration behavior of FGM cylindrical panels with large longitudinal edge lengths. ⁷ Three-dimensional elasticity theory and state space method were employed to investigate the thermo-elastic behavior of FGM cylindrical panels.^8,9 Nonlinear buckling analysis of FGM doubly curved panels was performed using the HSDT with the perturbation method. ¹⁰ Navier’s solution for linear dynamic responses of FGM doubly curved panels was mentioned using the FSDT with Sander’s curvature assumption. ¹¹ By using the Galerkin and Runge-Kutta methods with Donnell shell theory, the nonlinear vibration and dynamic buckling responses of FGM doubly curved panels were studied. ¹² The bending response of FGM doubly curved shells was performed by Kolapkar and Sayyad ¹³ using the Navier solution technique.

Nanocomposite materials incorporate matrix material and nanomaterials to enhance mechanical, thermal, and electrical properties. Widely used in aerospace, automotive, electronics, and biomedical fields, they significantly improve the strength and functionality of the matrix material. Functionally Graded Carbon Nanotube-Reinforced Composites (FG-CNTRCs) are advanced materials that incorporate carbon nanotubes into a composite matrix with a gradient distribution to enhance mechanical properties. These materials offer superior strength, stiffness, and thermal stability, making them ideal for aerospace, automotive, and structural applications. The HSDT with perturbation method was considered to investigate the postbuckling behavior, ¹⁴ thermal buckling and postbuckling behavior, ¹⁵ and thermally postbuckled vibration behavior ¹⁶ of FG-CNTRC plates. The shear buckling and dynamic responses of FG-CNTRC plates and cylindrical panels were performed using the Ritz formulations with FSDT.^17,18 The dynamic stability responses and wave propagation behavior of FG-CNTRC with the viscous damping effect of FG-CNTRC plates were mentioned with an orthotropic elastomeric medium ¹⁹ and with a sensor and actuator. ²⁰ By using the FSDT and the Ritz method, FG-CNTRC thick quadrilateral plates were considered in linear vibration responses. ²¹ The tangential-edge restraints were considered to evaluate the boundary condition effects on the nonlinear thermal buckling behavior of FG-CNTRC plates and cylindrical panels using the Galerkin method.^22,23 The static analysis of FG-CNTRC cylinder shells with various boundary restrictions was investigated using singular trigonometric series and the Laplace transform. ²⁴ The FG-CNTRC cylindrical shells with stiffeners ²⁵ and toroidal shell segments with corrugated core ²⁶ were considered in global torsional, compressed, and tensile buckling problems using Donnell shell theory with von Karman nonlinearities.

Functionally Graded Graphene Platelet-Reinforced Composites (FG-GPLRCs) are advanced materials that utilize graphene platelets distributed with a gradient to improve mechanical and thermal properties. GPLs can be incorporated into metal matrices, such as copper ²⁷ and aluminum, ²⁸ thermoplastic polymer matrices, such as polyurethane, ²⁹ polycarbonate, ³⁰ and Poly(methyl methacrylate) (PMMA), ³¹ as well as thermosetting polymer matrices, such as epoxy. ³² These composites are known for their exceptional strength, electrical conductivity, and lightweight nature, making them suitable for high-performance engineering applications. By using various methods and theories, the vibration, large deflection bending, nonlinear buckling, and stress analysis problems of plates and panels were performed in many works.^32–39 By using the HSDT and Ritz procedure, the nonlinear buckling and vibration problems of FG-GPLRC circular plates and spherical shells with and without porous core were investigated.^40–44

Piezoelectric materials are substances capable of generating or responding to electricity when placed in an electric field. They are used in sensors, automatic controls, and the electronics industry. The mechanical behavior of composite plates and shells with piezoelectric layers was mentioned in many works.^45–47 Auxetic materials are unique materials that exhibit a negative Poisson’s ratio, meaning they become thicker perpendicular to an applied force when stretched. This counterintuitive property provides enhanced energy absorption, indentation resistance, and mechanical robustness, making them valuable for various engineering applications. The auxetic material can be designed as a separate layer or as a core for advanced composite surface layers, and its special effects have also been evaluated.^48–53 The FG-GPLRC plates and cylindrical panels with auxetic core were mentioned in vibration and static problems, taking into account the hygrothermal effects.^54,55

Although extensive studies were reported on FG-GPLRC structures, auxetic cores, and piezoelectric layers individually, the combination of these features into a unified structure remains challenging, especially under piezo-thermo-mechanical conditions and complex curvature. The development of analytical models capable of analyzing such multi-functional interactions offers opportunities for flexible and efficient engineering design. In addition, the combination of FG-GPLRC layers, auxetic core, and piezoelectric layer allows the complementary advantages of these materials to be effectively exploited within a unified structure. The nonlinear HSDT is used to analyze the nonlinear electro-thermo-mechanical buckling and postbuckling behaviors of FG-GPLRC complexly curved panels with porous core and piezoelectric layers. Different distributions of graphene platelet (GPL), and different types of curvatures, including parabola, sinusoidal, and cylindrical curvatures, are investigated. The mechanical, thermal, and electrical components are combined simultaneously in Hooke’s law equation and the potential energy expression. To overcome the difficulty in obtaining exact stress functions arising from complex curvature, an approximate stress function technique is applied. The Ritz energy method is used to obtain the critical buckling and postbuckling responses of panels. Numerical examples present different effects of the GPL distribution and volume fraction, curvature of panels, piezoelectric layers, foundation, porosity coefficient, and uniformly distributed temperature change on the buckling responses of the considered panels.

Configuration properties and governing formulations

Consider the panels with complex curvature with the thickness $h$ , lengths $a$ and $b$ , and mid-plane rise $H_{ri}$ . Three configuration types are included parabola, sinusoidal, and cylindrical panels. The Cartesian coordinate system $Oxyz$ can be approximately obtained as in Figure 1, based on the assumption of shallow curvature of the panels.

Figure 1.

Configuration properties of panels, auxetic cell, and GPL distributions.

The geometrical equation and radius equation of the curved edges in the $y$ -direction for the parabola are presented as

$z_{pp}=\frac{4H_{ri}}{b^2}{y}^2-\frac{4H_{ri}}{b}y,R_{pp}=\frac{{\left[b^4+16H_{ri}^2{\left(b-2y\right)}^2\right]}^{3/2}}{8b^4{H}_{ri}}.$

(1)

Similarly, the geometrical equation and radius equation of the curved edges in the $y$ -direction for the sinusoidal panels can be obtained by

$z_{sp}=-H_{ri}\hspace{0.17em}\sin \left(\mathrm{\pi }y/b\right),R_{sp}=\frac{{\left[H_{ri}^2{\mathrm{\pi }}^2\hspace{0.17em}{\cos }^2\left(y\mathrm{\pi }/b\right)+b^2\right]}^{3/2}}{b{H}_{ri}{\mathrm{\pi }}^2\hspace{0.17em}\sin \left(y\mathrm{\pi }/b\right)}.$

(2)

The extended Halpin-Tsai model can be used to anticipate the elastic modulus of FG-GPLRC layers, as

$E_{gp}=\frac{3+3{\mathrm{\kappa }}_1^{gp}{\mathrm{\tau }}_1^{gp}{V}_{GPL}^{*}}{8-8{\mathrm{\tau }}_1^{gp}{V}_{GPL}^{*}}{E}_m^{*}+\frac{5+5{\mathrm{\kappa }}_2^{gp}{\mathrm{\tau }}_2^{gp}{V}_{GPL}^{*}}{8-8{\mathrm{\tau }}_2^{gp}{V}_{GPL}^{*}}{E}_m^{*},$

(3)

where

${\mathrm{\tau }}_1^{gp}=\frac{E_{GPL}^{*}/E_m^{*}-1}{E_{GPL}^{*}/E_m^{*}+{\mathrm{\kappa }}_1^{gp}},{\mathrm{\tau }}_2^{gp}=\frac{E_{GPL}^{*}/E_m^{*}-1}{E_{GPL}^{*}/E_m^{*}+{\mathrm{\kappa }}_2^{gp}},{\mathrm{\kappa }}_1^{gp}=2a_{GPL}/t_{GPL},{\mathrm{\kappa }}_2^{gp}=2b_{GPL}/t_{GPL},$

with $E_m^{*}$ and $E_{GPL}^{*}$ are the elastic modulus of the isotropic matrix and GPL, respectively. The length, width, and thickness of the GPL are denoted by $a_{GPL},b_{GPL},$ and $t_{GPL}$ , respectively. The GPLs volume fraction $V_{GPL}^{*}$ is calculated as

$V_{GPL}^{*}=\frac{W_{GPL}^{*}}{W_{GPL}^{*}+\left({\mathrm{\rho }}_{GPL}/{\mathrm{\rho }}_M\right)\left(1-W_{GPL}^{*}\right)},$

(4)

where ${\mathrm{\rho }}_M$ and ${\mathrm{\rho }}_{GPL}$ are the mass densities of the matrix and GPLs, respectively.

The GPLs mass fraction $W_{GPL}^{*}$ respected to the thickness direction of the FG-GPLRC layers by the following functions

+) The upper layer $\left(\frac{-h}{2}+h_{pe}\le z\le \frac{-h_{au}}{2}\right)$

$-\mathrm{X}\hspace{0.17em}\text{distribution}:W_{GPL}^{*}=\left(\left|\frac{8z-4h_{pe}+2h+2h_{au}}{h-h_{au}-2h_{pe}}\right|\right){\overline{W}}_{GPL}$

(5)

$-\mathrm{O}\hspace{0.17em}\text{distribution}:W_{GPL}^{*}=\left(2-\left|\frac{8z-4h_{pe}+2h+2h_{au}}{h-h_{au}-2h_{pe}}\right|\right){\overline{W}}_{GPL}$

(6)

$-\mathrm{V}\hspace{0.17em}\text{distribution}:W_{GPL}^{*}=\left(\frac{-4z-2h_{au}}{h-h_{au}-2h_{pe}}\right){\overline{W}}_{GPL}$

(7)

$-\mathrm{\Lambda }\hspace{0.17em}\text{distribution}:W_{GPL}^{*}=\left(\frac{4z+2h-4h_{pe}}{h-h_{au}-2h_{pe}}\right){\overline{W}}_{GPL}$

(8)

+) The lower layer $\left(\frac{h_{au}}{2}\le z\le -h_{pe}+\frac{h}{2}\right)$

$-\mathrm{X}\hspace{0.17em}\text{distribution}:W_{GPL}^{*}=\left(\left|\frac{8z+4h_{pe}-2h-2h_{au}}{h-h_{au}-2h_{pe}}\right|\right){\overline{W}}_{GPL}$

(9)

$-\mathrm{O}\hspace{0.17em}\text{distribution}:W_{GPL}^{*}=\left(2-\left|\frac{8z+4h_{pe}-2h-2h_{au}}{h-h_{au}-2h_{pe}}\right|\right){\overline{W}}_{GPL}$

(10)

$-\mathrm{V}\hspace{0.17em}\text{distribution}:W_{GPL}^{*}=\left(\frac{4z-2h_{au}}{h-h_{au}-2h_{pe}}\right){\overline{W}}_{GPL}$

(11)

$-\mathrm{\Lambda }\hspace{0.17em}\text{distribution}:W_{GPL}^{*}=\left(\frac{-4z+2h-4h_{pe}}{h-h_{au}-2h_{pe}}\right){\overline{W}}_{GPL}$

(12)

The expressions of the Poisson ratio, coefficient of thermal expansion, and density of FG-GPLRC layers are anticipated using the mixture rule as

$\begin{array}{l}{\mathrm{\nu }}_{gp}={\mathrm{\nu }}_m\left(1-V_{GPL}^{*}\right)+{\mathrm{\nu }}_{GPL}{V}_{GPL}^{*},{\mathrm{\alpha }}_{gp}={\mathrm{\alpha }}_m\left(1-V_{GPL}^{*}\right)+{\mathrm{\alpha }}_{GPL}{V}_{GPL}^{*},\\ {\mathrm{\rho }}_{gp}={\mathrm{\rho }}_m\left(1-V_{GPL}^{*}\right)+{\mathrm{\rho }}_{GPL}{V}_{GPL}^{*}.\end{array}$

(13)

The piezoelectric layers are assumed to be PZT-5A with $E_{pe}=63$ GPa, $G_{pe}=24.2$ GPa, ${\mathrm{\nu }}_{pe}=0.3$ , ${\mathrm{\alpha }}_{pe}=0.9\times {10}^{-6}$ /K, ${\mathrm{\rho }}_{pe}=7600$ kg/m³, $d_{31}^{pe}=d_{32}^{pe}=2.54\times {10}^{-10}$ m/V. ⁵⁶

The auxetic layers are made of UD-GPLRC or isotropic material by ensuring that the materials of the auxetic layer and FG-GPLRC layers at the contact surfaces are the same. The elastic constants of auxetic material can be determined by

$E_{gp}^0={E_{gp}|}_{z=\pm h_{au}/2},{\mathrm{\nu }}_{gp}^0={{\mathrm{\nu }}_{gp}|}_{z=\pm h_{au}/2},{\mathrm{\alpha }}_{gp}^0={{\mathrm{\alpha }}_{gp}|}_{z=\pm h_{au}/2}.$

(14)

The orthotropic constants of auxetic structures are presented as

$\begin{array}{l}{E}_{11}^{au}=\frac{{\mathrm{\delta }}_2^3\left({\mathrm{\delta }}_1-\sin \hspace{0.17em}\mathrm{\beta }\right)}{{\cos }^3\hspace{0.17em}\mathrm{\beta }\left[1+\left({\tan }^2\hspace{0.17em}\mathrm{\beta }+{\mathrm{\delta }}_1\hspace{0.17em}{\sec }^2\hspace{0.17em}\mathrm{\beta }\right){\mathrm{\delta }}_2^2\right]}{E}_{gp}^0,E_{22}^{au}=\frac{{\mathrm{\delta }}_2^3}{\cos \hspace{0.17em}\mathrm{\beta }\left({\mathrm{\delta }}_1-\sin \hspace{0.17em}\mathrm{\beta }\right)\left({\tan }^2\hspace{0.17em}\mathrm{\beta }+{\mathrm{\delta }}_2^2\right)}{E}_{gp}^0,\\ G_{12}^{au}=\frac{{\mathrm{\delta }}_2^3}{{\mathrm{\delta }}_1\left(1+2{\mathrm{\delta }}_1\right)\cos \hspace{0.17em}\mathrm{\beta }}{E}_{gp}^0,G_{23}^{au}=\frac{{\mathrm{\delta }}_2\hspace{0.17em}\cos \hspace{0.17em}\mathrm{\beta }}{{\mathrm{\delta }}_1-\sin \hspace{0.17em}\mathrm{\beta }}\frac{E_{gp}^0}{2\left(1+{\mathrm{\nu }}_{gp}^0\right)},\\ G_{13}^{au}=\frac{{\mathrm{\delta }}_2}{4\left(1+{\mathrm{\nu }}_{gp}^0\right)\cos \hspace{0.17em}\mathrm{\beta }}\left[\frac{{\mathrm{\delta }}_1-\sin \hspace{0.17em}\mathrm{\beta }}{1+2{\mathrm{\delta }}_1}+\frac{{\mathrm{\delta }}_1+2\hspace{0.17em}{\sin }^2\hspace{0.17em}\mathrm{\beta }}{2\left({\mathrm{\delta }}_1-\sin \hspace{0.17em}\mathrm{\beta }\right)}\right]E_{gp}^0,\\ {\mathrm{\nu }}_{12}^{au}=-\frac{\sin \hspace{0.17em}\mathrm{\beta }\left(1-{\mathrm{\delta }}_2^2\right)\left({\mathrm{\delta }}_1-\sin \hspace{0.17em}\mathrm{\beta }\right)}{\left[1+\left({\mathrm{\delta }}_1\hspace{0.17em}{\sec }^2\hspace{0.17em}\mathrm{\beta }+{\tan }^2\hspace{0.17em}\mathrm{\beta }\right){\mathrm{\delta }}_2^2\right]{\cos }^2\hspace{0.17em}\mathrm{\beta }},{\mathrm{\nu }}_{21}^{au}=-\frac{\sin \hspace{0.17em}\mathrm{\beta }\left(1-{\mathrm{\delta }}_2^2\right)}{\left({\tan }^2\hspace{0.17em}\mathrm{\beta }+{\mathrm{\delta }}_2^2\right)\left({\mathrm{\delta }}_1-\sin \hspace{0.17em}\mathrm{\beta }\right)},{\mathrm{\alpha }}_{au}={\mathrm{\alpha }}_{gp}^0.\end{array}$

(15)

where ${\mathrm{\delta }}_1=s_2/s_1,{\mathrm{\delta }}_2=t_1/s_1.$

For the FG-GPLRC layers, the X, O, V, and Λ types of GPL distributions are considered. Piezoelectric layers, FG-GPLRC layers, and auxetic-core are combined into four models as PE/X/AU/X/PE, PE/O/AU/O/PE, PE/V/AU/Λ/PE, and PE/Λ/AU/V/PE.

Hooke’s law is applied to the FG-GPLRC layers and auxetic layer, taking into account the thermal strains, presented as

$\left[\begin{array}{c}{\mathrm{\sigma }}_x\\ {\mathrm{\sigma }}_y\\ {\mathrm{\sigma }}_{xy}\\ {\mathrm{\sigma }}_{xz}\\ {\mathrm{\sigma }}_{yz}\end{array}\right]=\left[\begin{array}{ccccc}{\overline{Q}}_{11}& {\overline{Q}}_{12}& 0& 0& 0\\ {\overline{Q}}_{12}& {\overline{Q}}_{22}& 0& 0& 0\\ 0& 0& {\overline{Q}}_{66}& 0& 0\\ 0& 0& 0& {\overline{Q}}_{44}& 0\\ 0& 0& 0& 0& {\overline{Q}}_{55}\end{array}\right]\left[\begin{array}{c}{\mathrm{\epsilon }}_x-{\mathrm{\alpha }}_{gp}\mathrm{\Delta }T\\ {\mathrm{\epsilon }}_y-{\mathrm{\alpha }}_{gp}\mathrm{\Delta }T\\ {\mathrm{\gamma }}_{xy}\\ {\mathrm{\gamma }}_{xz}\\ {\mathrm{\gamma }}_{yz}\end{array}\right],$

(16)

and for the piezoelectric layers, as

$\left[\begin{array}{c}{\mathrm{\sigma }}_x\\ {\mathrm{\sigma }}_y\\ {\mathrm{\sigma }}_{xy}\\ {\mathrm{\sigma }}_{xz}\\ {\mathrm{\sigma }}_{yz}\end{array}\right]=\left[\begin{array}{ccccc}{Q}_{11}^{pe}& Q_{12}^{pe}& 0& 0& 0\\ Q_{12}^{pe}& Q_{22}^{pe}& 0& 0& 0\\ 0& 0& Q_{66}^{pe}& 0& 0\\ 0& 0& 0& Q_{44}^{pe}& 0\\ 0& 0& 0& 0& Q_{55}^{pe}\end{array}\right]\left\{\left[\begin{array}{c}{\mathrm{\epsilon }}_x\\ {\mathrm{\epsilon }}_y\\ {\mathrm{\gamma }}_{xy}\\ {\mathrm{\gamma }}_{xz}\\ {\mathrm{\gamma }}_{yz}\end{array}\right]-\left[\begin{array}{c}{\mathrm{\alpha }}_{pe}\\ {\mathrm{\alpha }}_{pe}\\ 0\\ 0\\ 0\end{array}\right]\mathrm{\Delta }T-\left[\begin{array}{c}{d}_{31}^{pe}/h_{pe}\\ d_{32}^{pe}/h_{pe}\\ 0\\ 0\\ 0\end{array}\right]V_0\right\},$

(17)

where $\mathrm{\Delta }T$ is the temperature difference from the initial state of the free thermal stress shells to the final state of deformed panels, and the reduced stiffnesses can be determined for GPLRC layers as

${\overline{Q}}_{11}^{gp}={\overline{Q}}_{22}^{gp}=\frac{E_{gp}}{1-{\mathrm{\nu }}_{gp}^2},{\overline{Q}}_{12}^{gp}=\frac{E_{gp}{\mathrm{\nu }}_{gp}}{1-{\mathrm{\nu }}_{gp}^2},{\overline{Q}}_{44}^{gp}={\overline{Q}}_{55}^{gp}={\overline{Q}}_{66}^{gp}=\frac{E_{gp}}{2+2{\mathrm{\nu }}_{gp}},$

(18)

for the auxetic-core

$\begin{array}{l}{\overline{Q}}_{11}^{au}=\frac{E_{11}^{au}}{1-{\mathrm{\nu }}_{12}^{au}{\mathrm{\nu }}_{21}^{au}},{\overline{Q}}_{22}^{au}=\frac{E_{22}^{au}}{1-{\mathrm{\nu }}_{12}^{au}{\mathrm{\nu }}_{21}^{au}},{\overline{Q}}_{12}^{au}=\frac{E_{11}^{au}{\mathrm{\nu }}_{12}^{au}}{1-{\mathrm{\nu }}_{12}^{au}{\mathrm{\nu }}_{21}^{au}},\\ {\overline{Q}}_{44}^{au}=G_{13}^{au},{\overline{Q}}_{55}^{au}=G_{23}^{au},{\overline{Q}}_{66}^{au}=G_{12}^{au}\end{array}$

(19)

and for the piezoelectric layers

$Q_{11}^{pe}=Q_{22}^{pe}=\frac{E_{pe}}{1-{\mathrm{\nu }}_{pe}^2},Q_{12}^{pe}=\frac{E_{pe}{\mathrm{\nu }}_{pe}}{1-{\mathrm{\nu }}_{pe}^2},Q_{44}^{pe}=Q_{55}^{pe}=Q_{66}^{pe}=G_{pe}.$

(20)

In this paper, a uniform electric voltage $V_0$ is applied across the thickness of the piezoelectric layer, inducing in-plane strains through the piezoelectric coefficients $d_{31}^{pe}$ and $d_{32}^{pe}$ , providing electro-mechanical coupling that controls the pre-buckling stress state and consequently the buckling and postbuckling behavior of the panel.

The nonlinear large deflection theory and HSDT are combined to formula the governing equations of the considered structures. The strain-displacement relations are presented by

$\begin{array}{l}\left\{\begin{array}{l}{\mathrm{\epsilon }}_x\\ {\mathrm{\epsilon }}_y\\ {\mathrm{\gamma }}_{xy}\end{array}\right\}=\left\{\begin{array}{l}{\mathrm{\epsilon }}_{0x}^{*}\\ {\mathrm{\epsilon }}_{0y}^{*}\\ {\mathrm{\gamma }}_{0xy}^{*}\end{array}\right\}+z\left\{\begin{array}{l}{\mathrm{\epsilon }}_{1x}^{*}\\ {\mathrm{\epsilon }}_{1y}^{*}\\ {\mathrm{\gamma }}_{1xy}^{*}\end{array}\right\}+z^3\left\{\begin{array}{l}{\mathrm{\epsilon }}_{3x}^{*}\\ {\mathrm{\epsilon }}_{3y}^{*}\\ {\mathrm{\gamma }}_{3xy}^{*}\end{array}\right\},\\ \left\{\begin{array}{l}{\mathrm{\gamma }}_{xz}\\ {\mathrm{\gamma }}_{yz}\end{array}\right\}=\left\{\begin{array}{l}{\mathrm{\gamma }}_{0xz}^{*}\\ {\mathrm{\gamma }}_{0yz}^{*}\end{array}\right\}+z^2\left\{\begin{array}{l}{\mathrm{\gamma }}_{2xz}^{*}\\ {\mathrm{\gamma }}_{2yz}^{*}\end{array}\right\},\end{array}$

(21)

where

$\begin{array}{l}\left\{\begin{array}{l}{\mathrm{\epsilon }}_{0x}^{*}\\ {\mathrm{\epsilon }}_{0y}^{*}\\ {\mathrm{\gamma }}_{0xy}^{*}\end{array}\right\}=\left\{\begin{array}{l}{u}_{,x}+w_{,x}^2/2+w_{,x}{\overline{w}}_{,x}\\ v_{,y}+w_{,y}^2/2+w_{,y}{\overline{w}}_{,y}-w/R\left(y\right)\\ v_{,x}+u_{,y}+w_{,x}{w}_{,y}+w_{,y}{\overline{w}}_{,x}+w_{,x}{\overline{w}}_{,y}\end{array}\right\},\left\{\begin{array}{l}{\mathrm{\epsilon }}_{1x}^{*}\\ {\mathrm{\epsilon }}_{1y}^{*}\\ {\mathrm{\gamma }}_{1xy}^{*}\end{array}\right\}=\left\{\begin{array}{c}{\mathrm{\varphi }}_{x,x}\\ {\mathrm{\varphi }}_{y,y}\\ {\mathrm{\varphi }}_{y,x}+{\mathrm{\varphi }}_{x,y}\end{array}\right\},\overline{\mathrm{\lambda }}=\frac{4}{3h^2},\\ \left\{\begin{array}{l}{\mathrm{\epsilon }}_{3x}^{*}\\ {\mathrm{\epsilon }}_{3y}^{*}\\ {\mathrm{\gamma }}_{3xy}^{*}\end{array}\right\}=-\overline{\mathrm{\lambda }}\left\{\begin{array}{c}{w}_{,xx}+{\mathrm{\varphi }}_{x,x}\\ w_{,yy}+{\mathrm{\varphi }}_{y,y}\\ {\mathrm{\varphi }}_{y,x}+2w_{,xy}+{\mathrm{\varphi }}_{x,y}\end{array}\right\},\left\{\begin{array}{l}{\mathrm{\gamma }}_{0xz}^{*}\\ {\mathrm{\gamma }}_{0yz}^{*}\end{array}\right\}=\left\{\begin{array}{l}{w}_{,x}+{\mathrm{\varphi }}_x\\ w_{,y}+{\mathrm{\varphi }}_y\end{array}\right\},\left\{\begin{array}{l}{\mathrm{\gamma }}_{2xz}^{*}\\ {\mathrm{\gamma }}_{2yz}^{*}\end{array}\right\}=-3\overline{\mathrm{\lambda }}\left\{\begin{array}{l}{w}_{,x}+{\mathrm{\varphi }}_x\\ w_{,y}+{\mathrm{\varphi }}_y\end{array}\right\}.\end{array}$

(22)

By integrating Hooke’s law through the panel thickness, the force and moment expressions are established as

$\begin{array}{l}{N}_x=L_{11}{\mathrm{\epsilon }}_{0x}^{*}+L_{12}{\mathrm{\epsilon }}_{0y}^{*}-{\mathrm{\Phi }}_{1x}^T\mathrm{\Delta }T-{\mathrm{\Phi }}_{1x}^{pe}{V}_0,N_y=L_{12}{\mathrm{\epsilon }}_{0x}^{*}+L_{22}{\mathrm{\epsilon }}_{0y}^{*}-{\mathrm{\Phi }}_{1y}^T\mathrm{\Delta }T-{\mathrm{\Phi }}_{1y}^{pe}{V}_0,\\ N_{xy}=L_{66}{\mathrm{\gamma }}_{0xy}^{*},M_x=F_{11}{\mathrm{\epsilon }}_{1x}^{*}+F_{12}{\mathrm{\epsilon }}_{1y}^{*}+A_{11}{\mathrm{\epsilon }}_{3x}^{*}+A_{12}{\mathrm{\epsilon }}_{3y}^{*},\\ M_y=F_{12}{\mathrm{\epsilon }}_{1x}^{*}+F_{22}{\mathrm{\epsilon }}_{1y}^{*}+A_{12}{\mathrm{\epsilon }}_{3x}^{*}+A_{22}{\mathrm{\epsilon }}_{3y}^{*},M_{xy}=F_{66}{\mathrm{\gamma }}_{1xy}^{*}+A_{66}{\mathrm{\gamma }}_{3xy}^{*},\\ T_x=A_{11}{\mathrm{\epsilon }}_{1x}^{*}+A_{12}{\mathrm{\epsilon }}_{1y}^{*}+D_{11}{\mathrm{\epsilon }}_{3x}^{*}+D_{12}{\mathrm{\epsilon }}_{3y}^{*},T_y=A_{12}{\mathrm{\epsilon }}_{1x}^{*}+A_{22}{\mathrm{\epsilon }}_{1y}^{*}+D_{12}{\mathrm{\epsilon }}_{3x}^{*}+D_{22}{\mathrm{\epsilon }}_{3y}^{*},\\ T_{xy}=A_{66}{\mathrm{\gamma }}_{1xy}^{*}+D_{66}{\mathrm{\gamma }}_{3xy}^{*},\end{array}$

(23)

The shear forces $Q_x,Q_y$ and the higher-order shear forces $S_x,S_y$ are expressed by

$\begin{array}{l}{Q}_x=C_{44}{\mathrm{\gamma }}_{0xz}^{*}+B_{44}{\mathrm{\gamma }}_{2xz}^{*},Q_y=C_{55}{\mathrm{\gamma }}_{0yz}^{*}+B_{55}{\mathrm{\gamma }}_{2yz}^{*}\\ S_x=B_{44}{\mathrm{\gamma }}_{0xz}^{*}+E_{44}{\mathrm{\gamma }}_{2xz}^{*},S_y=B_{55}{\mathrm{\gamma }}_{0yz}^{*}+E_{55}{\mathrm{\gamma }}_{2yz}^{*}\end{array}$

(24)

In equations (23) and (24), the stiffnesses, thermal and electric forces are presented by

$\begin{array}{l}\left(L_{ij},F_{ij},A_{ij},D_{ij}\right)={{\displaystyle \int }}_{-h/2}^{-h/2+h_{pe}}{Q}_{ij}^{pe}\left(1,z^2,z^4,z^6\right)\mathrm{d}z+{{\displaystyle \int }}_{-h/2+h_{pe}}^{-h_{au}/2}{\overline{Q}}_{ij}^{gp}\left(1,z^2,z^4,z^6\right)\mathrm{d}z\\ +{{\displaystyle \int }}_{-h_{au}/2}^{h_{au}/2}{\overline{Q}}_{ij}^{au}\left(1,z^2,z^4,z^6\right)\mathrm{d}z+{{\displaystyle \int }}_{h_{au}/2}^{-h_{pe}+h/2}{\overline{Q}}_{ij}^{gp}\left(1,z^2,z^4,z^6\right)\mathrm{d}z\\ +{{\displaystyle \int }}_{-h_{pe}+h/2}^{h/2}{Q}_{ij}^{pe}\left(1,z^2,z^4,z^6\right)\mathrm{d}z,\left(\hspace{0.17em}ij=1,2,6\hspace{0.17em}\right),\end{array}$

$\begin{array}{l}\left(C_{ij},B_{ij},E_{ij}\right)={\int }_{-h/2}^{-h/2+h_{pe}}{Q}_{ij}^{pe}\left(1,z^2,z^4\right)\mathrm{d}z+{\int }_{-h/2+h_{pe}}^{-h_{au}/2}{\overline{Q}}_{ij}^{gp}\left(1,z^2,z^4\right)\mathrm{d}z+{\int }_{-h_{au}/2}^{h_{au}/2}{\overline{Q}}_{ij}^{au}\left(1,z^2,z^4\right)\mathrm{d}z\\ +{\int }_{h_{au}/2}^{-h_{pe}+h/2}{\overline{Q}}_{ij}^{gp}\left(1,z^2,z^4\right)\mathrm{d}z+{\int }_{-h_{pe}+h/2}^{h/2}{Q}_{ij}^{pe}\left(1,z^2,z^4\right)\mathrm{d}z,\left(\hspace{0.17em}ij=4,5\right),\end{array}$

$\begin{array}{l}{\mathrm{\Phi }}_{1x}^T={\int }_{-h/2}^{-h/2+h_{pe}}\left(Q_{11}^{pe}+Q_{12}^{pe}\right){\mathrm{\alpha }}_{pe}\mathrm{d}z+{\int }_{-h/2+h_{pe}}^{-h_{au}/2}\left({\overline{Q}}_{11}^{gp}+{\overline{Q}}_{12}^{gp}\right){\mathrm{\alpha }}_{gp}\mathrm{d}z\\ +{\int }_{-h_{au}/2}^{h_{au}/2}\left({\overline{Q}}_{11}^{au}+{\overline{Q}}_{12}^{au}\right){\mathrm{\alpha }}_{au}\mathrm{d}z+{\int }_{h_{au}/2}^{-h_{pe}+h/2}\left({\overline{Q}}_{11}^{gp}+{\overline{Q}}_{12}^{gp}\right){\mathrm{\alpha }}_{gp}\mathrm{d}z+{\int }_{-h_{pe}+h/2}^{h/2}\left(Q_{11}^{pe}+Q_{12}^{pe}\right){\mathrm{\alpha }}_{pe}\mathrm{d}z,\end{array}$

$\begin{array}{l}{\mathrm{\Phi }}_{1y}^T={\int }_{-h/2}^{-h/2+h_{pe}}\left(Q_{12}^{pe}+Q_{22}^{pe}\right){\mathrm{\alpha }}_{pe}\mathrm{d}z+{\int }_{-h/2+h_{pe}}^{-h_{au}/2}\left({\overline{Q}}_{12}^{gp}+{\overline{Q}}_{22}^{gp}\right){\mathrm{\alpha }}_{gp}\mathrm{d}z\\ +{\int }_{-h_{au}/2}^{h_{au}/2}\left({\overline{Q}}_{12}^{au}+{\overline{Q}}_{22}^{au}\right){\mathrm{\alpha }}_{au}\mathrm{d}z+{\int }_{h_{au}/2}^{-h_{pe}+h/2}\left({\overline{Q}}_{12}^{gp}+{\overline{Q}}_{22}^{gp}\right){\mathrm{\alpha }}_{gp}\mathrm{d}z+{\int }_{-h_{pe}+h/2}^{h/2}\left(Q_{12}^{pe}+Q_{22}^{pe}\right){\mathrm{\alpha }}_{pe}\mathrm{d}z,\end{array}\begin{array}{l}{\mathrm{\Phi }}_{1x}^{pe}={\int }_{-h/2}^{-h/2+h_{pe}}\frac{Q_{11}^{pe}{d}_{31}^{pe}+Q_{12}^{pe}{d}_{32}^{pe}}{h_{pe}}\mathrm{d}z+{\int }_{-h_{pe}+h/2}^{h/2}\frac{Q_{11}^{pe}{d}_{31}^{pe}+Q_{12}^{pe}{d}_{32}^{pe}}{h_{pe}}\mathrm{d}z,\\ {\mathrm{\Phi }}_{1y}^{pe}={\int }_{-h/2}^{-h/2+h_{pe}}\frac{Q_{12}^{pe}{d}_{31}^{pe}+Q_{22}^{pe}{d}_{32}^{pe}}{h_{pe}}\mathrm{d}z+{\int }_{-h_{pe}+h/2}^{h/2}\frac{Q_{12}^{pe}{d}_{31}^{pe}+Q_{22}^{pe}{d}_{32}^{pe}}{h_{pe}}\mathrm{d}z.\end{array}$

By using equation (22), the nonlinear deformation compatibility equation can be established by

${\mathrm{\epsilon }}_{0x,yy}^{*}+{\mathrm{\epsilon }}_{0y,xx}^{*}-{\mathrm{\gamma }}_{0xy,xy}^{*}=w_{,xy}^2-\frac{w_{,xx}}{R\left(y\right)}-w_{,xx}{w}_{,yy}-w_{,xx}{\overline{w}}_{,yy}+2w_{,xy}{\overline{w}}_{,xy}-{\overline{w}}_{,xx}{w}_{,yy}.$

(25)

Introducing the stress function $\mathrm{\mu }\left(x,y\right)$ which satisfies

$N_x={\mathrm{\mu }}_{,yy},N_{xy}=-{\mathrm{\mu }}_{,xy},N_y={\mathrm{\mu }}_{,xx}.$

(26)

The nonlinear deformation compatibility equation (25) can be reestablished considering equations (26) and (23), as follows

$\begin{array}{l}\mathrm{\Upsilon }\equiv L_{22}^{*}{\mathrm{\mu }}_{,xxxx}+\left(L_{66}^{*}+2L_{12}^{*}\right){\mathrm{\mu }}_{,xxyy}+L_{11}^{*}{\mathrm{\mu }}_{,yyyy}+\frac{w_{,xx}}{R\left(y\right)}-w_{,xy}^2\\ +w_{,xx}{w}_{,yy}+w_{,xx}{\overline{w}}_{,yy}-2w_{,xy}{\overline{w}}_{,xy}+{\overline{w}}_{,xx}{w}_{,yy}=0,\end{array}$

(27)

where

$\mathrm{\Pi }=L_{11}{L}_{22}-L_{12}^2,L_{11}^{*}=L_{22}/\mathrm{\Pi },L_{12}^{*}=-L_{12}/\mathrm{\Pi },L_{22}^{*}=L_{11}/\mathrm{\Pi },L_{66}^{*}=1/L_{66}.$

Boundary conditions and solving problem

In this study, three forms of panels are considered with four boundary conditions.

Firstly, the four edges of the panels are freely movable and simply supported (4F), as

$\begin{array}{l}{w|}_{x=0,a}=0,N_x=N_{0x}=-h{P}_x,{T_x|}_{x=0,a}=0,{N_{xy}|}_{x=0,a}=0,{{\mathrm{\varphi }}_y|}_{x=0,a}=0,{M_x|}_{x=0,a}=0,\\ {w|}_{y=0,b}=0,N_y=N_{0y}=0,{T_y|}_{y=0,b}=0,{N_{xy}|}_{y=0,b}=0,{{\mathrm{\varphi }}_x|}_{y=0,b}=0,{M_y|}_{y=0,b}=0,\end{array}$

(28)

Secondly, all four edges are simply supported, where two edges $x=0,x=a$ are allowed to be freely movable and to be immovable for two edges $y=0,y=b$ (2F2I), as follows

$\begin{array}{l}{N}_x=-h{P}_x=N_{0x},{M_x|}_{x=0,a}=0,{T_x|}_{x=0,a}=0,{{\hspace{0.17em}{\mathrm{\varphi }}_y|}_{x=0,a}=0,N_{xy}|}_{x=0,a}=0,{w|}_{x=0,a}=0,\\ N_y=N_{0y},{M_y|}_{y=0,b}=0,{T_y|}_{y=0,b}=0,{{\mathrm{\varphi }}_x|}_{y=0,b}=0,{w|}_{y=0,b}=0,v=0,\end{array}$

(29)

Thirdly, all four edges are simply supported, where two edges $y=0,y=b$ are allowed to be freely movable and to be immovable for two edges $x=0,x=a$ (2I2F), as follows

$\begin{array}{l}{N}_x=N_{0x},{T_x|}_{x=0,a}=0,{M_x|}_{x=0,a}=0,{N_{xy}|}_{x=0,a}=0,{{\mathrm{\varphi }}_y|}_{x=0,a}=0,{w|}_{x=0,a}=0,\\ N_y=N_{0y}=0,{T_y|}_{y=0,b}=0,{M_y|}_{y=0,b}=0,u=0,{{\mathrm{\varphi }}_x|}_{y=0,b}=0,{w|}_{y=0,b}=0.\end{array}$

(30)

Finally, all edges of the panels are simply supported and immovable (4I), as

$\begin{array}{l}{{\mathrm{\varphi }}_y|}_{x=0,a}=0,N_x=N_{0x},{T_x|}_{x=0,a}=0,{M_x|}_{x=0,a}=0,u=0,{w|}_{x=0,a}=0,\\ {{\mathrm{\varphi }}_x|}_{y=0,b}=0,N_y=N_{0y},{T_y|}_{y=0,b}=0,{M_y|}_{y=0,b}=0,v=0,{w|}_{y=0,b}=0.\end{array}$

(31)

To satisfy the four boundary conditions (28-31), the approximate solutions for the deflection, rotations, and imperfect deflection of three forms of the panel are presented as

$\begin{array}{l}w=W\hspace{0.17em}\sin \hspace{0.17em}{a}_{1}x\hspace{0.17em}\sin \hspace{0.17em}{b}_{1}y,\overline{w}=\mathrm{\iota }h\hspace{0.17em}\sin \hspace{0.17em}{a}_{1}x\hspace{0.17em}\sin \hspace{0.17em}{b}_{1}y,\\ {\mathrm{\varphi }}_x={\mathrm{\Phi }}_x\hspace{0.17em}\cos \hspace{0.17em}{a}_{1}x\hspace{0.17em}\sin \hspace{0.17em}{b}_{1}y,{\mathrm{\varphi }}_y={\mathrm{\Phi }}_y\hspace{0.17em}\sin \hspace{0.17em}{a}_{1}x\hspace{0.17em}\cos \hspace{0.17em}{b}_{1}y,\end{array}$

(32)

where $a_1=m\mathrm{\pi }/a,b_1=n\mathrm{\pi }/b$ ; the numbers of half waves in the $x$ and $y$ directions are denoted by $m$ and $n$ , respectively, and $\mathrm{\iota }$ is the dimensionless imperfection of the panels.

In this paper, the imperfect deflection of panels $\overline{w}$ is assumed to be the same form as the deflection of panels $w$ .

The approximate stress function can be chosen in nonlinear form, as

$\mathrm{\mu }={\mathrm{\mu }}_1\hspace{0.17em}\cos \hspace{0.17em}2a_{1}x+{\mathrm{\mu }}_2\hspace{0.17em}\cos \hspace{0.17em}2b_{1}y+{\mathrm{\mu }}_3\hspace{0.17em}\sin \hspace{0.17em}{a}_{1}x\hspace{0.17em}\sin \hspace{0.17em}{b}_{1}y+\frac{1}{2}{N}_{0y}{x}^2+\frac{1}{2}{N}_{0x}{y}^2.$

(33)

Substituting equations (32) and (33) into equation (27), then, using the Galerkin method, the coefficients of stress function are obtained as

${\mathrm{\mu }}_1={\overline{V}}_{11}\left(W+2\mathrm{\iota }h\right)W+{\overline{V}}_{12}W,{\mathrm{\mu }}_2={\overline{V}}_{21}\left(W+2\mathrm{\iota }h\right)W+{\overline{V}}_{22}W,{\mathrm{\mu }}_3={\overline{V}}_{33}W,$

(34)

where

${\overline{V}}_{11}=\frac{b_1^2}{32a_1^2{L}_{22}^{*}},{\overline{V}}_{21}=\frac{a_1^2}{32b_1^2{L}_{11}^{*}},a_m={\left(-1\right)}^m-1,b_n={\left(-1\right)}^n-1,$

$\begin{array}{l}{\overline{V}}_{12}=\frac{-12ab{a}_1{b}_1{a}_m{b}_n{r}_{1\left(cp,pp,sp\right)}+\left(8a_m^2{b}_n^2-9a^2{b}^2{a}_1^2{b}_1^2\right)r_{2\left(cp,pp,sp\right)}-8a_m^2{b}_n^2{r}_{3\left(cp,pp,sp\right)}}{L_{22}^{*}ab{a}_1^4\left(72a^2{b}^2{a}_1^2{b}_1^2-128a_m^2{b}_n^2\right)},\\ {\overline{V}}_{22}=\frac{-12ab{a}_1{b}_1{a}_m{b}_n{r}_{1\left(cp,pp,sp\right)}-8a_m^2{b}_n^2{r}_{2\left(cp,pp,sp\right)}+\left(8a_m^2{b}_n^2-9a^2{b}^2{a}_1^2{b}_1^2\right)r_{3\left(cp,pp,sp\right)}}{L_{11}^{*}ab{b}_1^4\left(72a^2{b}^2{a}_1^2{b}_1^2-128a_m^2{b}_n^2\right)}\end{array}$

${\overline{V}}_{33}=-\frac{a_1{b}_1\left\{36ab{a}_1{b}_1{r}_{1\left(cp,pp,sp\right)}+24a_m{b}_n\left[r_{2\left(cp,pp,sp\right)}+r_{3\left(cp,pp,sp\right)}\right]\right\}}{\left[a_1^4{L}_{22}^{*}+b_1^2\left(2L_{12}^{*}+L_{66}^{*}\right)a_1^2+b_1^4{L}_{11}^{*}\right]\left(9a^2{b}^2{a}_1^2{b}_1^2-16a_m^2{b}_n^2\right)},$

$r_{1cp}=\frac{-a_1^{2}ab{H}_{ri}}{2\left[H_{ri}^2+{\left(b/2\right)}^2\right]},r_{2cp}=r_{3cp}=\frac{2H_{ri}{a}_1{a}_m{b}_n}{3b_1\left[H_{ri}^2+{\left(b/2\right)}^2\right]},r_{1\left(pp,sp\right)}=-\frac{1}{2}{a}_1^{2}a\underset{0}{\overset{b}{{\displaystyle \int }}}\frac{{\sin }^2\hspace{0.17em}{b}_{1}y}{R_{\left(pp,sp\right)}}\hspace{0.17em}\mathrm{d}y,$

$r_{2\left(pp,sp\right)}=-\frac{1}{3}{a}_1{a}_m\underset{0}{\overset{b}{{\displaystyle \int }}}\frac{\sin \hspace{0.17em}{b}_{1}y}{R_{\left(pp,sp\right)}}\hspace{0.17em}\mathrm{d}y,r_{3\left(pp,sp\right)}=a_1{a}_m\underset{0}{\overset{b}{{\displaystyle \int }}}\frac{\sin \hspace{0.17em}{b}_{1}y\hspace{0.17em}\cos \hspace{0.17em}2b_{1}y}{R_{\left(pp,sp\right)}}\hspace{0.17em}\mathrm{d}y.$

The immovable condition can be approximated by the following expressions, as

$\underset{0}{\overset{b}{{\displaystyle \int }}}\underset{0}{\overset{a}{{\displaystyle \int }}}{u}_{,x}\mathrm{d}x\mathrm{d}y=j_{41}W+j_{44}W\left(2\mathrm{\iota }h+W\right)+j_{45}{N}_{0x}+j_{46}{N}_{0y}+j_{47}\mathrm{\Delta }T+j_{48}{V}_0=0,$

(35)

$\underset{0}{\overset{b}{{\displaystyle \int }}}\underset{0}{\overset{a}{{\displaystyle \int }}}{v}_{,y}\mathrm{d}x\mathrm{d}y=j_{51}W+j_{54}W\left(2\mathrm{\iota }h+W\right)+j_{55}{N}_{0x}+j_{56}{N}_{0y}+j_{57}\mathrm{\Delta }T+j_{58}{V}_0=0,$

(36)

where

$j_{41}=-{\overline{V}}_{33}\left(a_1^2{L}_{12}^{*}+b_1^2{L}_{11}^{*}\right)a_m{b}_n/\left(a_1{b}_1\right),j_{44}=-ab{a}_1^2/8,j_{45}=ba{L}_{11}^{*},j_{46}=ba{L}_{12}^{*},$

$j_{47}=ab\left(L_{11}^{*}{\mathrm{\Phi }}_{1x}^T+L_{12}^{*}{\mathrm{\Phi }}_{1y}^T\right),j_{48}=ab\left(L_{11}^{*}{\mathrm{\Phi }}_{1x}^{pe}+L_{12}^{*}{\mathrm{\Phi }}_{1y}^{pe}\right),$

$j_{51}=-{\overline{V}}_{33}\left(a_1^2{L}_{22}^{*}+b_1^2{L}_{12}^{*}\right)a_m{b}_n/\left(a_1{b}_1\right)+r_{4\left(cp,pp,sp\right)},j_{54}=-ab{b}_1^2/8,j_{55}=ab{L}_{12}^{*},j_{56}=ab{L}_{22}^{*},$

$j_{57}=ab\left(L_{12}^{*}{\mathrm{\Phi }}_{1x}^T+L_{22}^{*}{\mathrm{\Phi }}_{1y}^T\right),j_{58}=ab\left(L_{12}^{*}{\mathrm{\Phi }}_{1x}^{pe}+L_{22}^{*}{\mathrm{\Phi }}_{1y}^{pe}\right),$

$r_{4cp}=\frac{2H_{ri}{a}_m\hspace{0.17em}{b}_n}{\left[H_{ri}^2+{\left(b/2\right)}^2\right]a_1\hspace{1em}{b}_1},r_{4\left(pp,sp\right)}=\frac{a_m}{a_1}\underset{0}{\overset{b}{{\displaystyle \int }}}\frac{\sin \hspace{0.17em}{b}_{1}y}{R_{\left(pp,sp\right)}}\hspace{0.17em}\mathrm{d}y.$

The total potential energy is expressed by

$U_{Total}^{*}=U_{in}^{*}-U_{ext}^{*}.$

(37)

where the strain energy of the panels can be established by

$U_{in}^{*}=\frac{1}{2}\underset{-h/2}{\overset{h/2}{{\displaystyle \int }}}\underset{0}{\overset{b}{{\displaystyle \int }}}\underset{0}{\overset{a}{{\displaystyle \int }}}\left[\begin{array}{l}{\mathrm{\sigma }}_{xy}{\mathrm{\gamma }}_{xy}+{\mathrm{\sigma }}_{xz}{\mathrm{\gamma }}_{xz}+{\mathrm{\sigma }}_{yz}{\mathrm{\gamma }}_{yz}+{\mathrm{\sigma }}_x{\mathrm{\epsilon }}_x+{\mathrm{\sigma }}_y{\mathrm{\epsilon }}_y\\ -\alpha \mathrm{\Delta }T\left({\mathrm{\sigma }}_x+{\mathrm{\sigma }}_y\right)-V_0\left(d_{31}^{pe}{\mathrm{\sigma }}_x+d_{32}^{pe}{\mathrm{\sigma }}_y\right)/h_{pe}\end{array}\right]\mathrm{d}x\mathrm{d}y\mathrm{d}z.$

(38)

and the work done by the external loads and foundation interaction is presented as

$\begin{array}{l}{U}_{ext}^{*}=-\underset{0}{\overset{b}{{\displaystyle \int }}}\underset{0}{\overset{a}{{\displaystyle \int }}}\left\{\frac{1}{2}w\left[K_{1}w-K_2\left(w_{,xx}+w_{,yy}\right)+\frac{K_3{w}^3}{2}\right]\right\}\mathrm{d}x\mathrm{d}y\\ +N_{0x}\underset{0}{\overset{b}{{\displaystyle \int }}}\underset{0}{\overset{a}{{\displaystyle \int }}}\left({\mathrm{\epsilon }}_{0x}^{*}-\frac{1}{2}{w}_{,x}^2-w_{,x}{\overline{w}}_{,x}\right)\mathrm{d}x\mathrm{d}y.\end{array}$

(39)

with $K_1$ (N/m³) and $K_2$ (N/m) are the Winkler and Pasternak parameters of the foundation. The nonlinear parameter of the foundation $K_3$ (N/m⁵) may be positive or negative values.

Applying the Ritz energy method as

$\frac{\partial U_{Total}^{*}}{\partial W}=\frac{\partial U_{Total}^{*}}{\partial {\mathrm{\Phi }}_x}=\frac{\partial U_{Total}^{*}}{\partial {\mathrm{\Phi }}_y}=0,$

(40)

leads to

$\begin{array}{l}{j}_{11}W+j_{12}{\mathrm{\Phi }}_x+j_{13}{\mathrm{\Phi }}_y+j_{16}W\left(W+4/3\mathrm{\iota }h\right)+j_{17}W\left(W+2\mathrm{\iota }h\right)\left(W+\mathrm{\iota }h\right)\\ +j_{19}{N}_{0x}\left(W+\mathrm{\iota }h\right)+j_{110}{N}_{0y}+j_{111}{K}_3{W}^3=0,\end{array}$

(41)

$j_{12}W+j_{22}{\mathrm{\Phi }}_x+j_{23}{\mathrm{\Phi }}_y=0,$

(42)

$j_{13}W+j_{23}{\mathrm{\Phi }}_x+j_{33}{\mathrm{\Phi }}_y=0,$

(43)

where

$j_{11}=\frac{ab}{4}\left\{\begin{array}{l}\left(32\hspace{0.17em}{L}_{22}^{*}{\overline{V}}_{12}^2+L_{22}^{*}{\overline{V}}_{33}^2+{\overline{\mathrm{\lambda }}}^2{D}_{11}\right)a_1^4+\left[\left(2D_{12}+4D_{66}\right){\overline{\mathrm{\lambda }}}^2+{\overline{V}}_{33}^2\left(2L_{12}^{*}+L_{66}^{*}\right)\right]a_1^2{b}_1^2\\ +\left(9\hspace{0.17em}{E}_{44}{\overline{\mathrm{\lambda }}}^2-6B_{44}\overline{\mathrm{\lambda }}+C_{44}\right)a_1^2+\left(32L_{11}^{*}{\overline{V}}_{22}^2+L_{11}^{*}{\overline{V}}_{33}^2+{\overline{\mathrm{\lambda }}}^2{D}_{22}\right)b_1^4\\ +\left(9\hspace{0.17em}{E}_{55}{\overline{\mathrm{\lambda }}}^2-6B_{55}\overline{\mathrm{\lambda }}+C_{55}\right)b_1^2+K_1+K_2\left(a_1^2+b_1^2\right)\\ -32\left[{\overline{V}}_{12}{L}_{22}^{*}{a}_1^4+b_1^2{L}_{12}^{*}\left({\overline{V}}_{12}+{\overline{V}}_{22}\right)a_1^2+b_1^4{\overline{V}}_{22}{L}_{11}^{*}\right]a_m{b}_n{\overline{V}}_{33}/\left(3ab{a}_1{b}_1\right)\end{array}\right\},$

$j_{12}=\frac{ab{a}_1}{4}\left\{\overline{\mathrm{\lambda }}\left(\overline{\mathrm{\lambda }}{D}_{11}-A_{11}\right)a_1^2+\left[\left(D_{12}+2D_{66}\right){\overline{\mathrm{\lambda }}}^2-\left(A_{12}+2A_{66}\right)\overline{\mathrm{\lambda }}\right]b_1^2+\left(9{\overline{\mathrm{\lambda }}}^2{E}_{44}-6B_{44}\overline{\mathrm{\lambda }}+C_{44}\right)\right\},$

$j_{13}=\frac{ab{b}_1}{4}\left\{\overline{\mathrm{\lambda }}\left(\overline{\mathrm{\lambda }}{D}_{22}-A_{22}\right)b_1^2+\left[\left(D_{12}+2D_{66}\right){\overline{\mathrm{\lambda }}}^2-\left(A_{12}+2A_{66}\right)\overline{\mathrm{\lambda }}\right]a_1^2+\left(9{\overline{\mathrm{\lambda }}}^2{E}_{55}-6\overline{\mathrm{\lambda }}{B}_{55}+C_{44}\right)\right\},$

$j_{17}=16ba\left({\overline{V}}_{11}^2{a}_1^4{L}_{22}^{*}+{\overline{V}}_{21}^2{b}_1^4{L}_{11}^{*}\right),j_{19}=a_1^{2}ab/4,$

$j_{16}=24ab\left({\overline{V}}_{11}{\overline{V}}_{12}{a}_1^4{L}_{22}^{*}+{\overline{V}}_{21}{\overline{V}}_{22}{b}_1^4{L}_{11}^{*}\right)-4\left[{\overline{V}}_{11}{a}_1^4{L}_{22}^{*}+a_1^2{b}_1^2\left({\overline{V}}_{11}+{\overline{V}}_{21}\right)L_{12}^{*}+b_1^4{\overline{V}}_{21}{L}_{11}^{*}\right]a_m{b}_n{\overline{V}}_{33}/\left(a_1{b}_1\right),$

$j_{110}=-{\overline{V}}_{33}\left(a_1^2{L}_{22}^{*}+b_1^2{L}_{12}^{*}\right)a_m{b}_n/\left(a_1{b}_1\right),j_{111}=9ab/64,$

$j_{22}=\frac{ab}{4}\left[\left({\overline{\mathrm{\lambda }}}^2{D}_{11}-2A_{11}\overline{\mathrm{\lambda }}+F_{11}\right)a_1^2+\left({\overline{\mathrm{\lambda }}}^2{D}_{66}-2A_{66}\overline{\mathrm{\lambda }}+F_{66}\right)b_1^2+9E_{44}{\overline{\mathrm{\lambda }}}^2-6B_{44}\overline{\mathrm{\lambda }}+C_{44}\right],$

$j_{23}=\frac{ab{a}_1{b}_1}{4}\left[\left(D_{12}+D_{66}\right){\overline{\mathrm{\lambda }}}^2-2\left(A_{12}+A_{66}\right)\overline{\mathrm{\lambda }}+F_{12}+F_{66}\right],$

$j_{33}=\frac{ab}{4}\left[\left({\overline{\mathrm{\lambda }}}^2{D}_{22}-2A_{22}\overline{\mathrm{\lambda }}+F_{22}\right)b_1^2+\left({\overline{\mathrm{\lambda }}}^2{D}_{66}-2A_{66}\overline{\mathrm{\lambda }}+F_{66}\right)a_1^2+9E_{55}{\overline{\mathrm{\lambda }}}^2-6B_{55}\overline{\mathrm{\lambda }}+C_{55}\right].$

By combining equations (41)–(43) with the conditions $N_{0x}=-P_{x}h$ , $N_{0y}=0$ , the 4F panels are obtained, meanwhile, the expressions of 2F2I panels are obtained by combining equations (41)–(43), equation (36), and the condition $N_{0x}=-P_{x}h$ is taken into account. For 2I2F panels, the compressive load is not considered, the expression of panels is obtained by combining equations (41)–(43), equation (35) with $N_{0y}=0$ . Equations (41)–(43), equations (35) and (36) are combined for the 4I panels.

The axially compressive postbuckling expression is obtained as

$\begin{array}{l}{P}_x=[\left({\mathrm{\eta }}_1{T}_{11}+{\mathrm{\eta }}_2{T}_{21}\right)W_1+j_{16}{W}_{1}h\left(W_1+4\mathrm{\iota }/3\right)+{\mathrm{\eta }}_2{T}_{24}{W}_{1}h\left(2\mathrm{\iota }+W_1\right)+\frac{{\mathrm{\eta }}_2}{h}\left(T_{28}\mathrm{\Delta }T+T_{29}{V}_0\right)\\ +\left({\mathrm{\eta }}_1+{\mathrm{\eta }}_2\right)j_{17}{h}^2{W}_1\left(2\mathrm{\iota }+W_1\right)\left(W_1+\mathrm{\iota }\right)+j_{111}{W_1}^3{h}^2{K}_3]/\left[{\mathrm{\eta }}_2{T}_{26}+\left({\mathrm{\eta }}_1+{\mathrm{\eta }}_2\right)j_{19}h\left(W_1+\mathrm{\iota }\right)\right]\end{array}$

(44)

where, ${\mathrm{\eta }}_1=1,{\mathrm{\eta }}_2=0$ for the case of the 4F condition, ${\mathrm{\eta }}_1=0,{\mathrm{\eta }}_2=1$ for the case of 2F2I.

The expression of the thermal postbuckling curve of the panels is achieved as

$\begin{array}{l}\mathrm{\Delta }T=-[\left({\mathrm{\eta }}_2{T}_{21}+{\mathrm{\eta }}_3{T}_{11}+{\mathrm{\eta }}_4{T}_{41}\right)W_1+\left({\mathrm{\eta }}_3{T}_{38}+{\mathrm{\eta }}_4{T}_{48}\right)V_0\left(\mathrm{\iota }+W_1\right)+j_{16}{W}_{1}h\left(W_1+4\mathrm{\iota }/3\right)\\ +\left({\mathrm{\eta }}_3{T}_{33}+{\mathrm{\eta }}_4{T}_{43}\right)W_{1}h\left(\mathrm{\iota }+W_1\right)+\left({\mathrm{\eta }}_2{j}_{17}+{\mathrm{\eta }}_3{T}_{35}+{\mathrm{\eta }}_4{T}_{45}\right)h^2{W}_1\left(2\mathrm{\iota }+W_1\right)\left(W_1+\mathrm{\iota }\right)\\ +\left({\mathrm{\eta }}_2{T}_{24}+{\mathrm{\eta }}_4{T}_{44}\right)W_{1}h\left(2\mathrm{\iota }+W_1\right)-{\mathrm{\eta }}_2{T}_{26}{P}_x+\left({\mathrm{\eta }}_2{T}_{29}+{\mathrm{\eta }}_4{T}_{47}\right)\frac{V_0}{h}\\ -{\mathrm{\eta }}_2{j}_{19}{P}_{x}h\left(W_1+\mathrm{\iota }\right)+j_{111}{W_1}^3{h}^2{K}_3]/\left[\frac{{\mathrm{\eta }}_2{T}_{28}}{h}+\frac{{\mathrm{\eta }}_4{T}_{46}}{h}+\left({\mathrm{\eta }}_3{T}_{32}+{\mathrm{\eta }}_4{T}_{42}\right)\left(\mathrm{\iota }+W_1\right)\right]\end{array}$

(45)

where, ${\mathrm{\eta }}_2=1,{\mathrm{\eta }}_3=0,{\mathrm{\eta }}_4=0$ for the case of 2F2I, ${\mathrm{\eta }}_2=0,{\mathrm{\eta }}_3=1,{\mathrm{\eta }}_4=0$ for the case of 2I2F, ${\mathrm{\eta }}_2=0,{\mathrm{\eta }}_3=0,{\mathrm{\eta }}_4=1$ for the case of 4I, and

$\begin{array}{l}{W}_1=W/h,T_{11}=j_{12}{I}_{14}+j_{13}{I}_{24}+j_{11},T_{21}=T_{11}-\frac{j_{110}{j}_{51}}{j_{56}},T_{24}=-\frac{j_{110}{j}_{54}}{j_{56}},T_{26}=-\frac{j_{110}{j}_{55}}{j_{56}},\\ T_{28}=-\frac{j_{110}{j}_{57}}{j_{56}},T_{29}=-\frac{j_{110}{j}_{58}}{j_{56}},T_{32}=-\frac{j_{19}{j}_{47}}{j_{45}},T_{38}=-\frac{j_{19}{j}_{48}}{j_{45}},T_{33}=-\frac{j_{19}{j}_{41}}{j_{45}},T_{35}=j_{17}-\frac{j_{19}{j}_{44}}{j_{45}},\\ T_{41}=T_{11}-\frac{j_{110}\left(I_{33}{I}_{42}+I_{32}\right)}{I_{33}{I}_{43}-1},T_{42}=-\frac{j_{19}\left(I_{34}{I}_{43}+I_{44}\right)}{I_{33}{I}_{43}-1},T_{48}=-\frac{j_{19}\left(I_{35}{I}_{43}+I_{45}\right)}{I_{33}{I}_{43}-1},\\ T_{43}=-\frac{j_{19}\left(I_{32}{I}_{43}+I_{42}\right)}{I_{33}{I}_{43}-1},T_{44}=-\frac{j_{110}\left(I_{33}{I}_{41}+I_{31}\right)}{I_{33}{I}_{43}-1},T_{45}=j_{17}-\frac{j_{19}\left(I_{31}{I}_{43}+I_{41}\right)}{I_{33}{I}_{43}-1},\\ T_{46}=-\frac{j_{110}\left(I_{33}{I}_{44}+I_{34}\right)}{I_{33}{I}_{43}-1},T_{47}=-\frac{j_{110}\left(I_{33}{I}_{45}+I_{35}\right)}{I_{33}{I}_{43}-1},\\ I_{14}=-\frac{j_{13}{j}_{23}-j_{33}{j}_{12}}{\hspace{1em}{j}_{23}^2-j_{33}{j}_{22}},I_{24}=\frac{j_{13}{j}_{22}-j_{23}{j}_{12}}{j_{23}^2-j_{33}{j}_{22}},I_{31}=-\frac{j_{54}}{j_{56}},I_{32}=-\frac{j_{51}}{j_{56}},I_{33}=-\frac{j_{55}}{j_{56}},\\ I_{34}=-\frac{j_{57}}{j_{56}},I_{35}=-\frac{j_{58}}{j_{56}},I_{41}=-\frac{j_{44}}{j_{45}},I_{42}=-\frac{j_{41}}{j_{45}},I_{43}=-\frac{j_{46}}{j_{45}},I_{44}=-\frac{j_{47}}{j_{46}},I_{45}=-\frac{j_{48}}{j_{46}}.\end{array}$

The critical buckling compressive load is obtained for 4F panels from equation (44), as

$P_x^{cr}=T_{11}/\left(j_{19}h\right).$

(46)

The critical temperature expression can be obtained from equation (45), as

$\mathrm{\Delta }{T}_{cr}=-\left(T_{11}+T_{38}{V}_0\right)/T_{32}.$

(47)

Results and discussion

To validate the temperature and mechanical buckling loads of the present study, the present results are compared with those of previous works. The critical buckling temperatures of FG-CNTRC plates are compared with those of previous results of Shen and Zhang ¹⁵ and Tung ²² in Table 1. As can be observed, the validations confirm the accuracy of the present work.

Table 1.

Comparisons of critical buckling temperatures $T_{cr}=T_0+\mathrm{\Delta }{T}_{cr}$ (in K) of for perfect CNTRC plates under uniform temperature rise (4I, $h$ = 2 mm, $b/h=20$ , $m=1$ , $n=1$ , $K_1$ = 0 N/m³, $K_2$ = 0 N/m, $T_0=300K$ ).

$a/b$	Source	$V_{CNT}^{*}$ = 0.12		$V_{CNT}^{*}$ = 0.17		$V_{CNT}^{*}$ = 0.28
		UD	X	UD	X	UD	X
1	Shen and Zhang 15	340.74	355.85	344.54	362.79	345.88	366.21
	Tung 22	352.89	380.31	355.39	386.64	362.56	404.75
	Present	336.92	350.34	340.65	357.00	341.39	359.75
1.5	Shen and Zhang 25	323.47	332.67	325.59	336.57	326.62	340.41
	Tung 22	326.31	338.43	326.90	340.88	329.65	348.93
	Present	322.31	330.76	324.67	334.95	324.97	337.97
2	Shen and Zhang 15	317.52	323.59	319.34	326.77	319.51	329.69
	Tung 22	318.53	325.50	318.23	326.42	319.25	330.98
	Present	317.74	323.60	320.05	327.46	319.28	329.49

The numerical investigations on the thermo-mechanical buckling and postbuckling responses of cylindrical, parabola, and sinusoidal panels are mentioned with various material and geometrical properties in this section. The graphene platelet, copper, and PMMA matrix are chosen, and their properties are used by Wang et al. ⁵⁷ and Shen and Zhang, ¹⁵ with $a_{GPL}=2.5\times {10}^{-6}$ m, $b_{GPL}=1.5\times {10}^{-6}$ m, $t_{GPL}=1.5\times {10}^{-9}$ m, $E_{GPL}^{*}=1010\times {10}^9$ N/m², ${\nu }_{GPL}=0.186$ , ${\alpha }_{GPL}=5\times {10}^{-6}$ K⁻¹, ${\rho }_{GPL}=1062.5$ kg/m³ for GPL, $E_m^{*}=130\times {10}^9$ N/m², ${\nu }_m=0.34$ , ${\alpha }_m=17\times {10}^{-6}$ K⁻¹, ${\rho }_m=8960$ kg/m³ for copper matrix, and $E_m^{*}=2.5\times {10}^9$ N/m², ${\nu }_m=0.34$ , ${\alpha }_m=51.75\times {10}^{-6}$ K⁻¹, ${\rho }_m=1150$ kg/m³ for PMMA matrix.

The critical buckling compressive loads $P_x^{cr}$ of the cylindrical, parabola, and sinusoidal 4F panels are considered in Table 2 with four distribution laws of mass fraction of copper- and PMMA-based FG-GPLRC panels. The small differences in compressive critical buckling loads of the three considered panel types can be observed at the same panel rise, and the small advantage of the sinusoidal panel over the other two panel types is unclearly expressed. Due to the assumption of shallow panels, the geometrical investigations show the small differences between panel geometries, leading to the small differences in critical buckling loads. The four distribution laws of mass fraction of FG-GPLRC layers largely affect the critical buckling loads of the panels. The critical buckling loads of PE/Λ/AU/V/PE panels are the smallest and those of PE/V/AU/Λ/PE panels are the largest, respectively. Although in current designs the auxetic core materials are different with different distribution laws of the FG-GPLRC layer, this observation still has many important implications for the designs of these structures. It is evident that the elastic modulus of copper is substantially higher than that of PMMA. As a result, the effective elastic modulus of FG-GPLRC with a copper matrix is considerably greater than that with a PMMA matrix, which in turn leads to a markedly higher critical buckling load for panels fabricated from copper-based FG-GPLRC compared to those with a PMMA matrix. Although the critical buckling loads of PMMA-based panels are significantly lower than those of copper-based panels, the variation trends with respect to different GPL distribution patterns and panel configurations remain consistent.

Table 2.

The critical buckling compressive load $P_x^{cr}$ (MPa) of FG-GPLRC with piezoelectric layers and auxetic core (4F, $h$ = 0.02 m, $h_{au}=0.56h$ , $h_{gp}=0.21h$ , $h_{pe}=0.01h$ , $a=b=30h$ , $H_{ri}=h$ , ${\overline{W}}_{GPL}$ = 0.5%, $V_0$ = 0V, $\mathrm{\iota }$ = 0, $\mathrm{\Delta }T$ = 0K, $\left(m,n\right)$ =(1,1), ${\mathrm{\delta }}_1$ = 1.5, ${\mathrm{\delta }}_2$ = 0.1, $\mathrm{\beta }$ = 35°, $K_1$ = 20 MN/m³, $K_2$ = 0.2 MN/m).

	Cylindrical panel	Parabola panel	Sinusoidal panel
Copper-based FG-GPLRC
PE/X/AU/X/PE	710.05	710.11	714.27
PE/Λ/AU/V/PE	690.93	690.99	695.15
PE/O/AU/O/PE	705.13	705.19	709.26
PE/V/AU/Λ/PE	723.78	723.84	727.91
PMMA-based FG-GPLRC
PE/X/AU/X/PE	208.89	208.88	210.10
PE/Λ/AU/V/PE	190.69	190.68	191.82
PE/O/AU/O/PE	203.80	203.79	204.96
PE/V/AU/Λ/PE	219.02	219.01	220.23

Table 3 presents the critical buckling temperatures $\mathrm{\Delta }{T}_{cr}$ of copper-based FG-GPLRC 2I2F panels with piezoelectric layers and auxetic core. Similar to the compressive load case, the small advantage of the sinusoidal panel over the other two panel types can also be obtained in the case of critical buckling temperatures. When the panel types are compared at an identical panel rise of the middle surface, their geometric variations are relatively minor; consequently, only slight differences in the critical compressive buckling loads are observed among the panel types. The critical temperature load of the PE/V/AU/Λ/PE panel reaches the outstanding value; on the contrary, the PE/Λ/AU/V/PE panel reaches the smallest value in both cases of auxetic core panel and solid core panel with equal material amount. It is clear that with the same amount of core material, the critical buckling temperatures of the auxetic core panels are significantly larger than those of the corresponding solid core panels. Due to the hollow nature of the auxetic core, the thickness of the auxetic core is significantly larger than that of the solid core for the same amount of material. This results in significantly larger eccentricities of the face layers relative to the middle surface of the auxetic core panel than those of the corresponding solid core panel, resulting in the overall stiffnesses of the auxetic core panels being superior to those of the corresponding solid core panels.

Table 3.

The critical buckling temperatures $\mathrm{\Delta }{T}_{cr}$ (K) of copper-based FG-GPLRC panels with piezoelectric layers and auxetic core (2I2F, $h$ = 0.02 m, $h_{au}=0.56h$ , $h_{gp}=0.21h$ , $h_{pe}=0.01h$ , $a=b=30h$ , $H_{ri}=h$ , ${\overline{W}}_{GPL}$ = 0.5%, $V_0$ = 0V, $\mathrm{\iota }$ = 0, $P_x$ = 0 MPa, $\left(m,n\right)$ =(1,1), ${\mathrm{\delta }}_1$ = 1.5, ${\mathrm{\delta }}_2$ = 0.1, $\mathrm{\beta }$ = 35°, $K_1$ = 20 MN/m³, $K_2$ = 0.2 MN/m).

		Cylindrical panel	Parabola panel	Sinusoidal panel
Auxetic core	PE/X/AU/X/PE	581.62	581.67	585.08
	PE/Λ/AU/V/PE	565.97	566.01	569.42
Solid core	PE/X/AU/X/PE	442.70	442.77	447.45
	PE/Λ/AU/V/PE	435.91	435.99	440.67
Auxetic core	PE/O/AU/O/PE	589.78	589.83	593.23
	PE/V/AU/Λ/PE	605.38	605.43	608.83
Solid core	PE/O/AU/O/PE	461.80	461.87	466.23
	PE/V/AU/Λ/PE	470.09	470.16	474.52

Effects of the geometrical properties of the auxetic core on the critical buckling compressive loads and temperatures of the sinusoidal panels with piezoelectric layers and FG-GPLRC auxetic core are presented in Table 4. Different values of incline angles $\mathrm{\beta }$ and geometrical parameters ${\mathrm{\delta }}_1$ of auxetic cores are considered in this investigation. While the incline angle $\mathrm{\beta }$ increases from 15° to 55°, the geometric parameter ${\mathrm{\delta }}_1$ increases from 1.5 to 3. The decrease in the critical buckling compressive loads and temperatures with increasing incline angle $\mathrm{\beta }$ and geometrical parameter ${\mathrm{\delta }}_1$ can be clearly observed in the calculation examples. In addition, the change can be observed more clearly with the critical temperature than with the critical compressive load.

Table 4.

Effects of the geometrical properties of auxetic core on the critical buckling compressive loads and temperatures of the sinusoidal panels with piezoelectric layers and copper-based FG-GPLRC auxetic core (PE/V/AU/Λ/PE, $h$ = 0.02 m, $h_{au}=0.56h$ , $h_{gp}=0.21h$ , $h_{pe}=0.01h$ , $a=b=30h$ , $H_{ri}=h$ , ${\overline{W}}_{GPL}$ = 0.5%, $V_0$ = 0V, $\mathrm{\iota }$ = 0, $\left(m,n\right)$ =(1,1), ${\mathrm{\delta }}_2$ = 0.1, $K_1$ = 20 MN/m³, $K_2$ = 0.2 MN/m).

	$\mathrm{\beta }$	${\mathrm{\delta }}_1$
		1.5	2	3
$P_x^{cr}$ (MPa)(4F)	15°	728.04	726.65	725.09
	35°	727.91	725.18	720.66
	55°	726.52	718.74	699.88
$\mathrm{\Delta }{T}_{cr}$ (K)(2I2F)	15°	624.17	620.99	614.24
	35°	608.83	594.22	567.31
	55°	586.28	548.58	511.75

Table 5 presents the effects of the auxetic core thickness of the core $h_{au}$ on the critical buckling compressive loads and temperatures of the panels with piezoelectric layers and FG-GPLRC auxetic core. Note that when the auxetic core thickness increases, the piezoelectric layer thickness remains the same and only the thickness of the FG-GPLRC layers decreases. As can be observed in this table, the critical buckling compressive loads decrease sharply, however, the critical buckling temperatures sharply increase as the auxetic core thickness increases. These trends can be observed for all three considered panel types. It can be observed that as the thickness of the auxetic core increases and of the FG-GPLRC layers decreases, the overall stiffnesses of the panel decreases, and the compressive critical buckling compressive load decreases sharply, however, this causes the overall internal thermal forces of the panel to decrease, leading to an increase in the critical buckling temperature. Unlike mechanical buckling, which is primarily governed by the stiffnesses, critical buckling temperature is controlled by the interaction between stiffnesses and the effective thermal expansion response of the layers with two opposite effects when the auxetic core thickness increases.

Table 5.

Effects of the auxetic core thickness of auxetic core on the critical buckling compressive loads and temperatures of the panels with piezoelectric layers and copper-based FG-GPLRC auxetic core (PE/X/AU/X/PE, $h$ = 0.02 m, $h_{gp}=\left(h-h_{au}-2h_{pe}\right)/2$ , $h_{pe}=0.01h$ , $a=b=30h$ , $H_{ri}=h$ , ${\overline{W}}_{GPL}$ = 0.5%, $V_0$ = 0V, $\mathrm{\iota }$ = 0, $\left(m,n\right)$ =(1,1), ${\mathrm{\delta }}_1$ = 1.5, ${\mathrm{\delta }}_2$ = 0.1, $\mathrm{\beta }$ = 35°, $K_1$ = 20 MN/m³, $K_2$ = 0.2 MN/m).

		$h_{au}$
		0.5 $h$	0.6 $h$	0.7 $h$
$P_x^{cr}$ (MPa)(4F)	Cylindrical panel	781.38	681.32	560.08
	Parabola panel	781.44	681.37	560.12
	Sinusoidal panel	786.04	685.09	562.95
$\mathrm{\Delta }{T}_{cr}$ (K)(2I2F)	Cylindrical panel	577.04	625.18	679.04
	Parabola panel	577.09	625.23	679.09
	Sinusoidal panel	580.49	628.64	682.53

Figure 2 presents the thermal and mechanical postbuckling responses of panels with different panel types, panel rises, and edge lengths. As can be seen in Figure 2(a) and (b), the compressive and thermal postbuckling curves of the cylindrical panel and the parabola panel are not significantly different for the same rise between the panels for both 4F and 2I2F boundary conditions and in both small and large deflection regions. Additionally, significant superiority of the postbuckling load-carrying capacity of the sinusoidal panels compared with those of the two remaining panel types can be observed. A similar tendency is also observed in Tables 2 and 3 for the corresponding critical buckling loads and temperatures. The slight differences reported for panel types with the same rise are attributed to their marginal geometric variations, which have a limited influence on the overall stiffnesses and buckling behavior. The snap-through can be clearly observed in the positive deflection region of the compressive postbuckling curve for both three-panel types (Figure 2(a)). Meanwhile, the snap-through can only be observed in the positive region of the deflection of the thermal postbuckling curve for the sinusoidal panel (Figure 2(b)). For the increase in the panel rises in Figure 2(c), it can be seen that the postbuckling curves of the panels increase obviously. The geometrical imperfections are also examined in this figure. With positive amplitude of imperfection, the postbuckling curve is obtained in the positive region of deflection, and the bifurcation point does not appear as in the cases of perfect panels. The panel thickness is kept the same and the edge lengths are changed in Figure 2(d), resulting in decreases in the mechanical and thermal postbuckling load-carrying capacity.

Figure 2.

Thermal and mechanical postbuckling responses of panels with different panel types, panel rises, and edge lengths. (a) Mechanical postbuckling with different panel types; (b) Thermal postbuckling with different panel types; (c) Mechanical postbuckling with different panel rises; (d) Thermal postbuckling with different edge lengths.

Figure 3 presents the thermal and mechanical postbuckling responses of panels with different boundary conditions, mass fraction of GPL, and auxetic core widths. As can be seen that the trends of mechanical and thermal postbuckling of panels with different boundary conditions are different in all investigated cases. The 4 F and 2F2I panels are considered for mechanical postbuckling behavior and the 4I and 2I2F panels are considered for thermal postbuckling behavior. The mass fraction of GPL greatly affects the postbuckling behavior of the panels as in Figure 3(a). A small amount of GPL added can significantly increase the postbuckling load-carrying capacity of the panels. The postbuckling curves of GPLRC auxetic core panels and GPLRC solid core panels are compared in Figure 3(b). Although the postbuckling load-carrying capacity of auxetic core panels is significantly larger than that of the corresponding solid core, the upward trends of the curves do not appear to be much different. The increased auxetic core thickness rapidly increases the mechanical and thermal postbuckling load-carrying capacity of the panel as shown in Figure 3(c)–(f), respectively. It is observed that the occurrence of snap-through strongly depends on the degree of in-plane restraint imposed by the boundary conditions (4F for mechanical load and 2I2F for thermal load). When the unloaded edges permit translational movement, the structure is able to develop the membrane state before snap-through. Conversely, when the unloaded edges are restrained, the combined effect of compressive membrane forces and the initial shell curvature induces deflection immediately upon loading; therefore, the bending develops from the early loading stage.

Figure 3.

Thermal and mechanical postbuckling responses of panels with different boundary conditions, mass fraction of GPL, and auxetic core widths. (a) Mechanical postbuckling with different mass fractions of GPL; (b) Thermal postbuckling with auxetic and solid cores; (c) Mechanical postbuckling with different core thicknesses (sinusoidal panel); (d) Mechanical postbuckling with different core thicknesses (parabola panel); (e) Thermal postbuckling with different core thicknesses (sinusoidal panel); (f) Thermal postbuckling with different core thicknesses (parabola panel).

Mechanical and thermal postbuckling responses of panels with different core geometry parameters are considered in Figure 4. Effects of incline angle $\mathrm{\beta }$ of auxetic core on the mechanical postbuckling behavior of panels are presented in Figure 4(a) and (b). When the incline angle increases from 15° to 55°, the load-carrying capacity decreases slightly. Effects of the geometrical parameter ${\mathrm{\delta }}_1$ on the thermal postbuckling behavior of panels are also investigated in Figure 4(c) and (d). As can be seen, the load-carrying capacity decreases significantly for thermal postbuckling cases. Although the postbuckling load-carrying capacities of the panels in these investigations varied significantly, the trends of the curves did not change much.

Figure 4.

Thermal and mechanical postbuckling responses of panels with different core geometry parameters. (a) Mechanical postbuckling with different incline angles (parabola panel); (b) Mechanical postbuckling with different incline angles (sinusoidal panel); (c) Thermal postbuckling with different geometrical parameters δ₁ (sinusoidal panel); (d) Thermal postbuckling with different geometrical parameters δ₁ (cylindrical panel).

Figure 5 investigates the thermal and mechanical postbuckling responses of panels with different pre-loads, voltages, and nonlinear foundation parameters. The pre-temperature applied to the panel creates significant negative pre-deflection, resulting in a marked reduction in the load-carrying capacity of the panels as investigated in Figure 5(a) and (b). Similarly, in Figure 5(c), the pre-compressive load also causes negative pre-deflection and reduces the thermal load capacity of the panel. The hardening foundation (positive nonlinear parameter) and softening foundation (negative nonlinear parameter) are considered in Figure 5(d). It is clear that the nonlinear elastic foundation parameters do not significantly affect the postbuckling behavior in the small deflection region, but strongly affect the large deflection region. The electrical voltage largely affects the postbuckling behavior of panels in Figure 5(e) and (f). Positive and negative voltages increase and decrease the load-carrying capacity of the panel, respectively. Positive or negative voltages may induce compressive or tensile membrane forces that combine with compressive thermal loads, thereby decreasing or increasing the effective stiffness and stabilizing the postbuckling response. This demonstrates that the postbuckling characteristics of the panel can be actively tailored via electrical control in the applied structures.

Figure 5.

Thermal and mechanical postbuckling responses of panels with different pre-loads, voltages, and nonlinear foundation parameters. (a) Mechanical postbuckling with different pre-temperatures (parabola panel); (b) Mechanical postbuckling with different pre-temperatures (cylindrical panel); (c) Thermal postbuckling with different pre-compressive loads; (d) Thermal postbuckling with different nonlinear parameters of foundation; (e) Thermal postbuckling with different voltages (parabola panel); (f) Thermal postbuckling with different voltages (sinusoidal panel).

Concluding remarks

The nonlinear large deflection assumption and the HSDT are used to analyze the buckling and postbuckling behaviors of FG-GPLRC complexly curved panels with porous core and piezoelectric layers. Four types of simply supported boundary conditions are considered, and the approximate stress function is determined using the Galerkin method. Additionally, the nonlinear three-parameter elastic foundation model and Ritz energy method are used to obtain the critical buckling and postbuckling responses of panels. The remarkable observations can be made as follows:

− The small differences in compressive critical buckling loads and postbuckling curves of the three considered panel types can be observed at the same panel rise, and the small advantage of the sinusoidal panel over the other two panel types is unclearly expressed.

− The decrease in the critical buckling compressive loads and temperatures with increasing the incline angle and geometrical parameter can be clearly observed in the calculation examples.

− The increased auxetic core thickness rapidly increases the mechanical and thermal critical buckling loads and postbuckling load-carrying capacity of the panel, respectively.

The present formulation allows for the analysis of a wide range of conditions, including mechanical and thermal loading conditions, auxetic and solid cores, different panel geometries, various elastic foundation models, and, notably, piezoelectric actuation for active control. Consequently, it provides flexible and comprehensive capacities for predicting and tailoring the buckling and postbuckling behavior of advanced composite panels in engineering applications.