Three-Dimensional Buckling Analysis of Functionally Graded Saturated Porous Rectangular Plates under Combined Loading Conditions

Abstract

The present work studies the buckling behavior of functionally graded (FG) porous rectangular plates subjected to different loading conditions. Three different porosity distributions are assumed throughout the thickness, namely, a nonlinear symmetric, a nonlinear asymmetric and a uniform distribution. A novel approach is proposed here based on a combination of the generalized differential quadrature (GDQ) method and finite elements (FEs), labeled here as the FE-GDQ method, while assuming a Biot’s constitutive law in lieu of the classical elasticity relations. A parametric study is performed systematically to study the sensitivity of the buckling response of porous structures, to different input parameters, such as the aspect ratio, porosity and Skempton coefficients, along with different boundary conditions (BCs) and porosity distributions, with promising and useful conclusions for design purposes of many engineering structural porous members.

1. Introduction

In the last decades, an increased interest in porous materials has arisen among scientists and designers regarding engineering materials and structures due to their remarkable mechanical properties, electrical conductivity and high permeability. Besides, porous materials can be used in the aerospace industry and sea structures because of their very low density, but also in submarines, reformers and catalysts owing to their high specific surfaces. Thus, many investigations on the mechanical behavior of functionally graded (FG) porous plate and shell structures have been increasingly conducted in the literature from a theoretical, experimental and computational standpoint. Biot [1] was one of the pioneers who investigated the buckling response of a fluid-saturated porous slab under an axial compression, and checked for the sensitivity of the buckling load to pore compressibility. Similarly, Magnucki and Stasiewicz [2] suggested an analytical determination of the critical buckling load of a compressed porous beam based on a broken-line hypothesis and the principle of stationary action for the total potential energy. A shear deformation theory was applied in [3] for the buckling study of porous beams with varying material properties, and in [4] for the bending and buckling of rectangular plates made of a foam material with a nonlinear symmetric porosity distribution. In the further work by Chen et al. [5], the elastic buckling behavior of shear deformable FG porous beams was studied systematically to check for the effect of different porosity distributions on the mechanical response. A multiple analytical, numerical and experimental approach was proposed by Jasion et al. [6] for the buckling study of plates and beams, with a foam core and external layers of perfect material. In the last decades, different higher-order assumptions have been integrated with high-performance computational methods to treat buckling problems of perfect and/or porous composite structures. A nonlinear dynamic buckling of FG porous beams was performed by Kang Gao et al. [7]. The Galerkin method was applied by the authors to determine the governing equations of the problem, which was then solved numerically by means of a fourth-order Runge–Kutta method. A different approach based on a generalized differential quadrature (GDQ) method was applied by Tang et al. [8] to analyze the nonlinear and linear buckling behavior of FG porous Euler–Bernoulli beams. A refined theory was also proposed by Ebrahimi and Jafari [9] to treat the buckling problem of smart magneto-electro-elastic-FG porous plates, accounting for two different FG distributions. Hyperbolic higher-order shear deformation theory (HSDT) was combined with a mesh-free approach in [10] to investigate the buckling and free vibration behavior of porous FG plates resting on an elastic foundation. Among coupled problems, Cong et al. [11] focused on the nonlinear thermomechanical buckling and post-buckling of porous FG plates with two poro/nonlinear symmetric and non-symmetric distributions, by using the Reddy’s HSDT and Galerkin method. Further recent contributions on the buckling and free vibration response of perfect and porous FG plates applied first-order shear deformation theory (FSDT) combined with the Chebyshev Polynomials-Ritz method [12,13,14,15], even for graphene-reinforced nanocomposites. Tu et al. [16], instead, proposed a Galerkin-based solution for the nonlinear buckling and post-buckling study of imperfect porous plates subjected to different mechanical loads while applying classical shell theory–Von Karman nonlinearity. In addition, Sekkal et al. [17] focused on a novel quasi-3D HSDT to assess the buckling and vibration response of FG plates, whose solution was determined analytically. Shahsavari et al. [18] investigated the shear buckling of porous nanoplates with even, uneven and logarithmic-uneven distribution templates, by means of the Galerkin method, a novel size-dependent quasi-3D shear deformation theory and Eringen’s nonlocal elasticity. Another successful application of mesh-free methods can be found in [19] for the thermal buckling response of porous sandwich plates with CNT-reinforced nanocomposite layers. At the same time, Li et al. [20] studied the nonlinear vibration and dynamic buckling of a sandwich FG porous plate reinforced by graphene platelets and resting on a Winkler–Pasternak elastic foundation, where the Galerkin method was proposed together with the fourth-order Runge–Kutta approach as theoretical and numerical tools. A conventional FSDT approach was also employed by Shahgholian et al. [21] for the study of the buckling behavior of FG graphene-reinforced porous cylindrical shells combined with the Rayleigh–Ritz numerical method, whereas Zhao et al. [22] applied the classical Euler–Bernoulli theory and the Galerkin method to check for the dynamic instability of FG porous arches reinforced by graphene platelets.

Based on the current literature on the buckling of FG porous structures, however, most studies rely on the use of simple elastic Hooke’s laws, with limited attention to the effect of pore fluid pressures stemming from poroelastic constitutive Biot’s laws. In such a context, Jabbari et al. [23,24] proposed a closed-form solution for the axial buckling of FG-saturated, porous, rectangular, simply supported Kirchhoff plates, immersed in a piezoelectric [23] or thermal field [24], respectively. In [25], the same authors studied the axisymmetric buckling of a saturated circular porous-cellular plate as provided by FSDT. In the further work by Jabbari et al. [26,27], classical plate theory (CPT) or HSDT was implemented for the analysis of the buckling capacity of circular porous plates under a radial compressive load, and its sensitivity to some important poroelastic material properties. In another work, Jabbari et al. [28] performed a buckling study of thin circular FG plates made of saturated porous-soft ferromagnetic materials in transverse magnetic fields, whereas in [29,30], a FSDT closed solution was proposed to the buckling problem of transversely graded saturated porous plates with piezoelectric layers, and the axisymmetric post-buckling study of saturated porous circular plates under a uniform radial compression. Among moderately thick plates, Rezaei and Saidi [31] assessed the buckling behavior of fluid-infiltrated porous annular sector plates, as provided by Mindlin plate theory involving fluid-saturated and fluid-free conditions. Additional buckling studies for structural members made of metals or composite materials can be found in [32,33,34,35]. The available literature, however, shows the potential application of buckling issues in many porous structures, although at the present state, a proper study of the buckling response of saturated porous rectangular plates subjected to normal and shear loads is still lacking. Based on the literature overview, it seems that the analysis of porous structures is usually based on FSDT and HSDT. Moreover, in most studies, the Hooke’s law or drained condition is commonly assumed to model the porous behavior of structures. Based on the above-mentioned lacking aspects of the problem, in this work, the buckling behavior is investigated for FG-saturated porous rectangular plates subjected to a double normal and shear loads. To this end, 3D elasticity theory and Biot’s constitutive law are applied, while proposing a mixed FE-DQM based on a Rayleigh–Ritz energy formulation as an efficient computational tool to solve the problem. The application of Biot’s constitutive law in lieu of the simple Hooke’s law provides more realistic results and conclusions, even from a practical standpoint. Based on the fact that plate theories overestimate the buckling loads for thick plates, 3D elasticity is implemented here to account for the thickness stretching effects, for the sake of accuracy, together with a more efficient mixed FE-GDQ method rather than conventional FEs. Three different porosity distributions are selected here in the thickness direction, namely, a nonlinear symmetric, a nonlinear asymmetric and uniform distribution. The objective of the work is to check the effects of different porosity distributions, as well as the porosity and Skempton coefficients on the critical buckling load of undrained rectangular plates with different geometrical dimensions and BCs, as useful for many engineering applications.

The remainder of the paper is structured as follows. In Section 2, the geometrical and mechanical properties of porous rectangular plates are briefly described, together with the governing equations of the problem, as determined by means of the virtual work principle and Biot’s constitutive poroelastic law. Section 3 presents the main basics of the mixed FE-GDQ numerical formulation, as proposed here to solve the problem, whose numerical examples are tested and discussed in Section 4 among a large systematic investigation. Conclusions are finally drawn in Section 5.

2. Theoretical Definition of the Problem

2.1. Poroelastic Modeling of Plates

Let us consider a rectangular FG porous plate, with in-plane dimensions $a$ and $b$, and thickness $h$, as depicted in Figure 1, along with three different porosity patterns throughout the thickness direction $\left(0\le z\le h\right)$, namely, a non-symmetric nonlinear porous distribution (PNND), a symmetric nonlinear porous distribution (PNSD) and a uniform porous distribution (PUD). Except for uniform porosities, the mechanical properties of the material in terms of shear modulus, Young’s modulus and mass density, for a PNND and PNSD, are defined as in Equation (1) [36,37,38,39,40].

$$\begin{gathered} E=E_1\left[1-e_0\hspace{0.17em}Q\right] \\ G=G_1\left[1-e_0\hspace{0.17em}Q\right] \\ \rho ={\rho }_1\left[1-e_m\hspace{0.17em}Q\right] \end{gathered}$$

(1)

In which

$$Q(z)=\{\begin{array}{l}(a)PNND\hspace{1em}\hspace{1em}\hspace{1em}cos\left(\frac{\pi z}{2h}\right)\\ \left(\mathrm{b}\right)-\mathrm{PNSD}\hspace{1em}\hspace{1em} cos\left(\frac{\pi }{2}-\frac{\pi z}{h}\right)\end{array}$$

(2)

and $0\le e_0\le 1$ is the porosity coefficient. Moreover, $E_1$, $G_1$ and ${\rho }_1$ denote the Young’s modulus, the shear modulus and the mass density at $z=h$ (for a PNND) and at $z=0$ (for a PNSD), whereby $E_j=2G_j\left(1+\nu \right)$, $j=0,1$ and the Poisson’s ratio, $\nu$, is assumed to be constant in the z-direction. The constitutive equations of FG-saturated porous rectangular plates are derived from Biot’s theory, which accounts for the displacements field of the solid, the pore fluid movement as well as their interactions owing to the applied loads [41]. Based on Biot’s theory, the constitutive law is thus written as [42]:

$${\sigma }_{ij}=2G{\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}-p\alpha {\delta }_{ij}$$

(3)

where

$$\begin{gathered} p=\overline{M}\left(\Psi -\alpha {\epsilon }_{kk}\right) \\ \overline{M}=\frac{2G\left(v_u-v\right)}{{\alpha }^2\left(1-2v_u\right)\left(1-2v\right)} \\ {v}_u=\frac{v+\alpha \beta \left(1-2v\right)/3}{1-\alpha \beta \left(1-2v\right)/3} \\ {v}_{=}\frac{{\epsilon }_{jj}}{{\epsilon }_{ii}}{|}_{{\sigma }_{ii}=0},p=0,i\ne \mathrm{j} \\ {v}_u{}_{=}\frac{{\epsilon }_{jj}}{{\epsilon }_{ii}}{|}_{{\sigma }_{ii}=0},\Psi =0,i\ne \mathrm{j} \\ \alpha =1-\frac{K}{K_S} \\ K=\frac{2\left(1+v\right)}{3\left(1-2v\right)}G \\ {K}_u=\frac{2\left(1+v_u\right)}{3\left(1-2v_u\right)}G \end{gathered}$$

Note that $p$ is the pore fluid pressure, such that, for $p=0$, Biot’s law reverts to the classical Hooke’s law (or drained condition). In addition, $\lambda$ denotes the Lamè constant, ${\delta }_{ij}$ is the Kronecker delta and $\alpha$ is the Biot’s effective stress coefficient (with $0<\alpha <1$). This parameter accounts for the effect of porosity on the structural behavior and resistance of porous materials in the absence of an internal fluid. At the same time, $\overline{M},\hspace{0.17em}G,\hspace{0.17em}\hspace{0.17em}{\nu }_u,\hspace{0.17em}\hspace{0.17em}{\epsilon }_{kk},\hspace{0.17em}\Psi ,\hspace{0.17em}{K}_s$ and $\beta$ stand for the Biot’s modulus, shear modulus, undrained Poisson’s ratio $\left(\nu <{\nu }_u<0.5\right)$, volumetric strain, variation of fluid volume content, bulk modulus of a homogeneous material and the Skempton coefficient, which introduces the pore fluid property, respectively. This last coefficient, $\beta$, in particular, denotes a dimensionless parameter to include the impact of a fluid within cavities on the overall response of a porous material in undrained condition ($\Psi =0$), and it is described as the ratio of the cavity pressure to the total body stress, namely,

$$\beta ={\frac{dp}{d\sigma }|}_{\Psi =0}=\frac{1}{1+e_0\frac{C_P}{C_S}}=\frac{K_u-K}{\alpha K_u}$$

(4)

In which $K_u,\hspace{0.17em}K$ refer to the bulk modulus in undrained and drained conditions respectively, and $C_P$ and $C_S$ stand for the fluid and solid compressibility in pores. Thus, the Skempton coefficient defines the effect of fluid compressibility on the elastic modulus and compressibility of the whole porous material [42].

2.2. Governing Equations

The governing equations of the problem are derived here from the principle of virtual work, as follows:

$$\delta U-\delta V_g=0$$

(5)

where $U$ is the total strain potential energy of the plate defined on the domain $\Omega$ as:

$$U=\frac{1}{2}{{\displaystyle \int }}_{\Omega }{\sigma }_{ij}{\epsilon }_{ij}d\Omega$$

(6)

and $V_g$ is the potential energy related to geometry, which takes the following form:

$$V_g=\frac{1}{2}{{\displaystyle \int }}_{\Omega }\left\{\begin{array}{l}{P}_x\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial w}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial x}\right)}^2\right]+P_y\left[{\left(\frac{\partial u}{\partial y}\right)}^2+{\left(\frac{\partial w}{\partial y}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2\right]+\\ P_{xy}\left(\frac{\partial w}{\partial x}\frac{\partial w}{\partial y}+\frac{\partial v}{\partial x}\frac{\partial v}{\partial y}-\frac{\partial u}{\partial x}\frac{\partial u}{\partial y}\right)\end{array}\right\}d\Omega$$

(7)

By substitution of Equations (6) and (7) into Equation (5), the following relation is obtained:

$${\displaystyle {\int }_0^h{\displaystyle {\int }_0^a{\displaystyle {\int }_0^b\left[\begin{array}{l}\left({\sigma }_{xx}\delta {\epsilon }_{xx}+{\sigma }_{yy}\delta {\epsilon }_{yy}+{\sigma }_{zz}\delta {\epsilon }_{zz}+{\sigma }_{yz}\delta {\gamma }_{yz}+{\sigma }_{xz}\delta {\gamma }_{xz}+{\sigma }_{xy}\delta {\gamma }_{xy}\right)\\ -P_x\left(\frac{\partial u}{\partial x}\frac{\partial \delta u}{\partial x}+\frac{\partial v}{\partial x}\frac{\partial \delta v}{\partial x}+\frac{\partial w}{\partial x}\frac{\partial \delta w}{\partial x}\right)-P_y\left(\frac{\partial u}{\partial y}\frac{\partial \delta u}{\partial y}+\frac{\partial v}{\partial y}\frac{\partial \delta v}{\partial y}+\frac{\partial w}{\partial y}\frac{\partial \delta w}{\partial y}\right)\\ -P_{xy}\left(\frac{\partial \delta w}{\partial x}\frac{\partial w}{\partial y}+\frac{\partial \delta w}{\partial y}\frac{\partial w}{\partial x}+\frac{\partial \delta v}{\partial x}\frac{\partial v}{\partial y}+\frac{\partial \delta v}{\partial y}\frac{\partial v}{\partial x}-\frac{\partial \delta u}{\partial x}\frac{\partial u}{\partial y}-\frac{\partial \delta u}{\partial y}\frac{\partial u}{\partial x}\right)\end{array}\right]}}}\hspace{0.17em}dv=0$$

(8)

where the constitutive relations for FG-saturated porous plates in 3D poroelasticity can be defined according to Biot’s constitutive law, as $[{\sigma }_{ij}]=\left[C\right][{\epsilon }_{ij}]$. For FG-saturated porous rectangular plates, the elasticity matrix reads as follows:

$$\begin{array}{l}C=\frac{E\left(1-\nu \right)}{\left(1+\nu \right)\left(1-2\nu \right)}\left(\begin{array}{cccccc}1& \frac{\nu }{1-\nu }& \frac{\nu }{1-\nu }& 0& 0& 0\\ \frac{\nu }{1-\nu }& 1& \frac{\nu }{1-\nu }& 0& 0& 0\\ \frac{\nu }{1-\nu }& \frac{\nu }{1-\nu }& 1& 0& 0& 0\\ 0& 0& 0& \frac{1-2\nu }{2(1-\nu )}& 0& 0\\ 0& 0& 0& 0& \frac{1-2\nu }{2(1-\nu )}& 0\\ 0& 0& 0& 0& 0& \frac{1-2\nu }{2(1-\nu )}\end{array}\right)+\left(\begin{array}{cccccc}\overline{M}{\alpha }^2& 0& 0& 0& 0& 0\\ 0& \overline{M}{\alpha }^2& 0& 0& 0& 0\\ 0& 0& \overline{M}{\alpha }^2& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)=\\ =E \Lambda +\Upsilon =\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{12}& 0& 0& 0\\ C_{12}& C_{11}& C_{12}& 0& 0& 0\\ C_{12}& C_{12}& C_{11}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right]\end{array}$$

(9)

The elasticity modulus, E, is assumed to vary along the z-direction, whereas the Poisson’s ratio, $\nu$, remains constant.

3. Mixed FE-GDQ Numerical Formulation

A set of rectangular-quadratic elements, $N_e$, is considered to discretize the x–y plane of the domain. Each element should be differentiable at both in-plane and transverse displacements for solving the governing equations. Hence, the displacement components are approximated as:

$$\begin{array}{l}u\left(x,y,z,t\right)={\displaystyle \sum _{j=1}^N}{\phi }_j\left(x,y\right)U_j\left(z,t\right),\\ v\left(x,y,z,t\right)={\displaystyle \sum _{j=1}^N}{\phi }_j\left(x,y\right)V_j\left(z,t\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,2,\dots ,N\\ w\left(x,y,z,t\right)={\displaystyle \sum _{j=1}^N}{\phi }_j\left(x,y\right)W_j\left(z,t\right)\end{array}$$

(10)

where $N$ stands for the total number of nodes in the discretized x–y plane, and ${\phi }_j\left(x,y\right)$ denotes the global Lagrange interpolation functions. By combination of Equations (3) and (8), and integrating by parts in the z-direction, the following governing equations per each node $i$ ($i=1,2,\dots ,N$) are obtained:

$$\begin{array}{l}\delta U_i:{\displaystyle \sum _{j=1}^N{A}^{ij}{C}_{55}(z_j)\frac{{\partial }^2{U}_j}{\partial z^2}+\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{j=1}^N\frac{\partial C_{55}(z_j)}{\partial z}}\hspace{0.17em}\frac{\partial U_j}{\partial z}+{\displaystyle \sum _{j=1}^N\left(D^{ij}{C}_{55}(z_j)-D^{ji}{C}_{12}(z_j)\right)\frac{\partial W_j}{\partial z}}}\\ +{\displaystyle \sum _{j=1}^N\frac{\partial C_{44}}{\partial z}(z_j)\hspace{0.17em}{W}_j}-{\displaystyle \sum _{j=1}^N\left(E^{ji}{C}_{12}(z_j)+E{\hspace{0.17em}}^{ij}{C}_{66}(z_j)\right)}{U}_j\hspace{0.17em}-\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{j=1}^N\left(F^{ij}{C}_{11}(z_j)+H^{ji}{C}_{66}(z_j)\right)V_j}\\ -P_x{\displaystyle \sum _{j=1}^N{D}_b{}^{ij}{V}_j}-P_y{\displaystyle \sum _{j=1}^N{B}_b{}^{ij}}{V}_j-P_{xy}{\displaystyle \sum _{j=1}^N{S}_b{}^{ij}{V}_j-P_{xy}{\displaystyle \sum _{j=1}^N{K}_b{}^{ij}{V}_j=0}}\end{array}$$

(11)

$$\begin{array}{l}\delta V_i:{\displaystyle \sum _{j=1}^N{A}^{ij}\hspace{0.17em}\hspace{0.17em}{C}_{44}(z_j)\frac{{\partial }^2{V}_j}{\partial z^2}+\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{j=1}^N\frac{\partial C_{44}(z_j)}{\partial z}}\hspace{0.17em}\frac{\partial V_j}{\partial z}+{\displaystyle \sum _{j=1}^N\left(D^{ij}\hspace{0.17em}{C}_{44}(z_j)-D^{ji}\hspace{0.17em}{C}_{12}(z_j)\right)\frac{\partial W_j}{\partial z}}}\\ +{\displaystyle \sum _{j=1}^N\frac{\partial C_{44}}{\partial z}(z_j)\hspace{0.17em}{W}_j}-{\displaystyle \sum _{j=1}^N\left(E^{ji}\hspace{0.17em}\hspace{0.17em}{C}_{12}(z_j)+E{\hspace{0.17em}}^{ij}{C}_{66}(z_j)\right)}{U}_j\hspace{0.17em}-\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{j=1}^N\left(F^{ij}\hspace{0.17em}{C}_{11}(z_j)+H^{ji}\hspace{0.17em}{C}_{66}(z_j)\right)V_j}\\ -P_x{\displaystyle \sum _{j=1}^N{D}_b{}^{ij}{V}_j}-P_y{\displaystyle \sum _{j=1}^N{B}_b{}^{ij}}{V}_j-P_{xy}{\displaystyle \sum _{j=1}^N{S}_b{}^{ij}{V}_j-P_{xy}{\displaystyle \sum _{j=1}^N{K}_b{}^{ij}{V}_j=0}}\end{array}$$

(12)

$$\begin{array}{l}\delta W_i:{\displaystyle \sum _{j=1}^N{A}^{ij}\hspace{0.17em}{C}_{11}(z_j)\hspace{0.17em}\frac{{\partial }^2{W}_j}{\partial z^2}+{\displaystyle \sum _{j=1}^N\left(B^{ij}{C}_{12}(z_j)-B^{ij}{C}_{55}(z_j)\right)\frac{\partial U_j}{\partial z}}}\\ +{\displaystyle \sum _{j=1}^N\left(D^{ij}{C}_{12}(z_j)-D^{ji}\hspace{0.17em}{C}_{44}(z_j)\right)}\frac{\partial V_j}{\partial z}-{\displaystyle \sum _{j=1}^N\hspace{0.17em}\hspace{0.17em}{B}^{ij}\frac{\partial C_{12}(z)}{\partial z}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{U}_j\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ +{\displaystyle \sum _{j=1}^N\hspace{0.17em}\hspace{0.17em}{D}^{ij}\frac{\partial C_{12}(z)}{\partial z}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{V}_j+{\displaystyle \sum _{j=1}^N\left(F^{ij}\hspace{0.17em}{C}_{44}(z_j)+H^{ij}{C}_{55}(z_j)+A^{ij}\frac{\partial C_{11}(z)}{\partial z}\hspace{0.17em}\right)W_j}\\ -P_x{\displaystyle \sum _{j=1}^N{D}_b{}^{ij}{W}_j}-P_y{\displaystyle \sum _{j=1}^N{B}_b{}^{ij}}{W}_j-P_{xy}{\displaystyle \sum _{j=1}^N{S}_b{}^{ij}{W}_j-P_{xy}{\displaystyle \sum _{j=1}^N{K}_b{}^{ij}{W}_j=0}}\end{array}$$

(13)

where,

$$\begin{gathered} H^{ij}={\displaystyle {\displaystyle \underset{0}{\overset{a}{\int }}}{\displaystyle {\displaystyle \underset{0}{\overset{b}{\int }}}\frac{\partial {\phi }_i}{\partial x}}}\frac{\partial {\phi }_j}{\partial x}dxdy,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{F}^{ij}={\displaystyle {\displaystyle \underset{0}{\overset{a}{\int }}}{\displaystyle {\displaystyle \underset{0}{\overset{b}{\int }}}\frac{\partial {\phi }_i}{\partial y}}}\frac{\partial {\phi }_j}{\partial y}dxdy,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{E}^{ij}={\displaystyle {\displaystyle \underset{0}{\overset{a}{\int }}}{\displaystyle {\displaystyle \underset{0}{\overset{b}{\int }}}\frac{\partial {\phi }_i}{\partial x}}}\frac{\partial {\phi }_j}{\partial y}dxdy \\ {H}^{ij}={\displaystyle {\displaystyle \underset{0}{\overset{a}{\int }}}{\displaystyle {\displaystyle \underset{0}{\overset{b}{\int }}}\frac{\partial {\phi }_i}{\partial x}}}\frac{\partial {\phi }_j}{\partial x}dxdy,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{F}^{ij}={\displaystyle {\displaystyle \underset{0}{\overset{a}{\int }}}{\displaystyle {\displaystyle \underset{0}{\overset{b}{\int }}}\frac{\partial {\phi }_i}{\partial y}}}\frac{\partial {\phi }_j}{\partial y}dxdy,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{E}^{ij}={\displaystyle {\displaystyle \underset{0}{\overset{a}{\int }}}{\displaystyle {\displaystyle \underset{0}{\overset{b}{\int }}}\frac{\partial {\phi }_i}{\partial x}}}\frac{\partial {\phi }_j}{\partial y}dxdy \\ {D}_b{}^{ij}={\displaystyle {\displaystyle \underset{0}{\overset{a}{\int }}}{\displaystyle {\displaystyle \underset{0}{\overset{b}{\int }}}\frac{\partial {\phi }_i}{\partial x}}}\frac{\partial {\phi }_j}{\partial x}dydx,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{B}_b{}^{ij}={\displaystyle {\displaystyle \underset{0}{\overset{a}{\int }}}{\displaystyle {\displaystyle \underset{0}{\overset{b}{\int }}}\frac{\partial {\phi }_i}{\partial y}}}\frac{\partial {\phi }_j}{\partial y}dydx,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{S}_b{}^{ij}={\displaystyle {\displaystyle \underset{0}{\overset{a}{\int }}}{\displaystyle {\displaystyle \underset{0}{\overset{b}{\int }}}\frac{\partial {\phi }_i}{\partial x}}}\frac{\partial {\phi }_j}{\partial y}dydx,\hspace{0.17em}\hspace{0.17em}{K}_b{}^{ij}={\displaystyle {\displaystyle \underset{0}{\overset{a}{\int }}}{\displaystyle {\displaystyle \underset{0}{\overset{b}{\int }}}\frac{\partial {\phi }_i}{\partial y}}}\frac{\partial {\phi }_j}{\partial x}dydx \end{gathered}$$

(14)

The BCs at the lower and upper surfaces ($z=0$ and $z=h$) associated with Equations (11)–(13) are defined as:

$$\mathrm{Either} \delta U_i=0,\hspace{1em}\mathrm{or}\{\begin{array}{c}{\displaystyle \sum _{j=1}^N}{A}^{ij}{C}_{55}(z_j)\frac{\partial U_j}{\partial z}+{\displaystyle \sum _{j=1}^N}{B}^{ij}\hspace{0.17em}{C}_{55}(z_{ij})W_j=0 \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{at} \hspace{0.17em}\hspace{0.17em}z=0,\\ {\displaystyle \sum _{j=1}^N}{A}^{ij}{C}_{55}(z_j)\frac{\partial U_j}{\partial z}+{\displaystyle \sum _{j=1}^N}{B}^{ij}\hspace{0.17em}{C}_{55}(z_{ij})W_j=0 \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{at} \hspace{0.17em}\hspace{0.17em}z=h,\end{array}$$

(15)

$$\mathrm{Either}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\delta V_i=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{or}\{\begin{array}{l}{\displaystyle \sum _{j=1}^N{A}^{ij}{C}_{44}(z_j)\frac{\partial V_j}{\partial z}+{\displaystyle \sum _{j=1}^N{D}^{ij}\hspace{0.17em}{C}_{44}(z_j)W_j=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{at}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}z=0,}}\\ {\displaystyle \sum _{j=1}^N{A}^{ij}{C}_{44}(z_j)\frac{\partial V_j}{\partial z}+{\displaystyle \sum _{j=1}^N{D}^{ij}\hspace{0.17em}{C}_{44}(z_j)W_j=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{at}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}z=h,}}\end{array}$$

(16)

$$\mathrm{Either}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\delta W_i=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{or}\{\begin{array}{l}{\displaystyle \sum _{j=1}^N{A}^{ij}{C}_{11}(z_j)\frac{\partial W_j}{\partial z}+{\displaystyle \sum _{j=1}^N{D}^{ij}\hspace{0.17em}{C}_{12}(z_j)V_j+{\displaystyle \sum _{j=1}^N{B}^{ij}{C}_{12}(z_j)U_j=0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{at}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}z=0}}\\ {\displaystyle \sum _{j=1}^N{A}^{ij}{C}_{11}(z_j)\frac{\partial W_j}{\partial z}+{\displaystyle \sum _{j=1}^N{D}^{ij}\hspace{0.17em}{C}_{12}(z_j)V_j+{\displaystyle \sum _{j=1}^N{B}^{ij}{C}_{12}(z_j)U_j=0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{at}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}z=h}}\end{array}$$

(17)

As far as the GDQ method is concerned, this approach discretizes the spatial derivatives of a function $f(z,t)$ as a weighted linear sum of the functional values at all nodes in the solution domain, by means of some fixed weighting coefficients. Thus, the first- and second-order derivatives of a one-dimensional function read as follows:

$$\begin{array}{c}{\frac{\partial f\left(z,t\right)}{\partial z}|}_{z=z_i}={\displaystyle \sum _{j=1}^{N_z}}{A}_{ij}^{z}f\left(z_j,t\right)={\displaystyle \sum _{j=1}^{N_z}}{A}_{ij}^z{f}_j\left(t\right)\\ {\frac{{\partial }^{2}f\left(z,t\right)}{\partial z^2}|}_{z=z_i}={\displaystyle \sum _{j=1}^{N_z}}{B}_{ij}^{z}f\left(z_j,t\right)={\displaystyle \sum _{j=1}^{N_z}}{B}_{ij}^z{f}_j\left(t\right)\end{array}$$

(18)

where $A_{ij}^z$ and $B_{ij}^z$ are the weighted coefficients at the grid nodes of the solution domain. To derive the weighting coefficients, the following relations are employed:

$$A_{ij}^z=\{\begin{array}{l}\frac{M\left(z_i\right)}{\left(z_i-z_j\right)M\left(z_i\right)}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i\ne j\\ -{\displaystyle \sum _{k=1,k\ne i}^{N_z}}{A}_{ik}^z\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=j\end{array}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i,j=1,2,\dots ,N_z,$$

(19)

$$B_{ij}^z=\{\begin{array}{l}2\left[A_{ii}^z{A}_{ij}^z-\frac{A_{ij}^z}{z_i-z_j}\right]\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i\ne j,\\ -{\displaystyle \sum _{k=1,k\ne i}^{N_z}}{B}_{ik}^z\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=j\end{array}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i,j=1,2,\dots ,N_z,$$

(20)

being

$$M^{\left(1\right)}\left(z_i\right)={\displaystyle \prod }_{j=1,j\ne i}^N\left(z_i-z_j\right) \hspace{0.17em}\mathrm{for}\hspace{0.17em}i=1,2,\dots ,N$$

To obtain more accurate results, a Chebyshev–Gauss–Lobatto quadrature-mesh size is assumed here, in line with findings by Malik and Bert [43]. At the current stage, the GDQ method is employed to discretize the system of equations through the thickness direction (i.e., along the z-axis). A set of $N_z$ grid points is assumed to discretize the domain along the thickness direction for each quadratic grid point. This means that Equations (11)–(13) can be rewritten in the domain (i.e., for each node $k=2,3,\dots ,N_z-1$), as follows:

$$\begin{array}{l}\delta U_i:{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_z}{A}^{ij}}{C}_{55}(z_j)\hspace{0.17em}{B}_{km}{U}_{jm}+\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{j=1}^N\frac{\partial C_{55}(z_j)}{\partial z}}\hspace{0.17em}{A}_{km}^z{U}_{jm}+{\displaystyle \sum _{j=1}^N\hspace{0.17em}{\displaystyle \sum _{m=1}^{N_Z}\left(B^{ij}{C}_{55}(z_j)-B^{ij}\hspace{0.17em}{C}_{12}(z_j)\right)}{A}_{km}^z{W}_{jm}}}\\ +{\displaystyle \sum _{j=1}^N\frac{\partial C_{44}}{\partial z}(z_j)\hspace{0.17em}{W}_j}-{\displaystyle \sum _{j=1}^N\left(E^{ji}{C}_{12}(z_j)+E{\hspace{0.17em}}^{ij}{C}_{66}(z_j)\right)}{U}_j\hspace{0.17em}-\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{j=1}^N\left(F^{ij}{C}_{11}(z_j)+H^{ji}{C}_{66}(z_j)\right)V_j}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ -P_x{\displaystyle \sum _{j=1}^N{D}_b{}^{ij}{V}_j}-P_y{\displaystyle \sum _{j=1}^N{B}_b{}^{ij}}{V}_j-P_{xy}{\displaystyle \sum _{j=1}^N{S}_b{}^{ij}{V}_j-P_{xy}{\displaystyle \sum _{j=1}^N{K}_b{}^{ij}{V}_j=0}}\end{array}$$

(21)

$$\begin{array}{l}\delta V_i:{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_Z}{A}^{ij}}{C}_{44}(z_j)B_{km}{V}_{jm}+\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_Z}\frac{\partial C_{44}(z_j)}{\partial z}}{A}_{km}^z{V}_{km}}\hspace{0.17em}+{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_Z}\left(D^{ij}{C}_{44}(z_j)-D^{ji}{C}_{12}(z_j)\right)}{A}_{km}^z{W}_{km}}}\\ +{\displaystyle \sum _{j=1}^N\frac{\partial C_{44}}{\partial z}(z_j)\hspace{0.17em}{W}_j}-{\displaystyle \sum _{j=1}^N\left(E^{ji}{C}_{12}(z_j)+E^{ij}{C}_{66}(z_j)\right)}{U}_j\hspace{0.17em}-\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{j=1}^N\left(F^{ij}{C}_{11}(z_j)+H^{ji}{C}_{66}(z_j)\right)V_j}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ -P_x{\displaystyle \sum _{j=1}^N{D}_b{}^{ij}{V}_j}-P_y{\displaystyle \sum _{j=1}^N{B}_b{}^{ij}}{V}_j-P_{xy}{\displaystyle \sum _{j=1}^N{S}_b{}^{ij}{V}_j-P_{xy}{\displaystyle \sum _{j=1}^N{K}_b{}^{ij}{V}_j=0}}\end{array}$$

(22)

$$\begin{array}{l}\delta W_i:{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_Z}{A}^{ij}}{C}_{11}(z_j)\hspace{0.17em}{B}_{km}^z{W}_{jm}+{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_Z}\left(B^{ij}{C}_{12}(z_j)-B^{ij}{C}_{55}(z_j)\right)}\hspace{0.17em}{A}_{km}^z{U}_{jm}}}\\ +{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_z}\left(D^{ij}{C}_{12}(z_j)-D^{ji}{C}_{44}(z_j)\right)}}{A}_{km}^z{V}_{jm}-{\displaystyle \sum _{j=1}^N\hspace{0.17em}\hspace{0.17em}{B}^{ij}\frac{\partial C_{12}(z)}{\partial z}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{U}_j+{\displaystyle \sum _{j=1}^N\hspace{0.17em}\hspace{0.17em}{D}^{ij}\frac{\partial C_{12}(z)}{\partial z}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{V}_j\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ +{\displaystyle \sum _{j=1}^N\left(F^{ij}{C}_{44}(z_j)+H^{ij}{C}_{55}(z_j)+A^{ij}\frac{\partial C_{11}(z)}{\partial z}\hspace{0.17em}\right)W_j}-P_x{\displaystyle \sum _{j=1}^N{D}_b{}^{ij}{W}_j}-P_y{\displaystyle \sum _{j=1}^N{B}_b{}^{ij}}{W}_j\\ -P_{xy}{\displaystyle \sum _{j=1}^N{S}_b{}^{ij}{W}_j-P_{xy}{\displaystyle \sum _{j=1}^N{K}_b{}^{ij}{W}_j=0}}\end{array}$$

(23)

Likewise, the BCs in Equations (15)–(17) at the top and bottom sides of the structures take the following form:

$$\mathrm{Either}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{U}_{ik}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{or}\{\begin{array}{l}{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_z}{A}^{ij}\hspace{0.17em}{C}_{55}(z_j)A_{km}^z{U}_{jm}+}{\displaystyle \sum _{j=1}^N{B}^{ij}{C}_{55}(z_j)W_{jk}=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1}}\\ {\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_z}{A}^{ij}\hspace{0.17em}{C}_{55}(z_j)A_{km}^z{U}_{jm}+}{\displaystyle \sum _{j=1}^N{B}^{ij}{C}_{55}(z_j)W_{jk}=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=N_Z}}\end{array}$$

(24)

$$\mathrm{Either}\hspace{0.17em}\hspace{0.17em}{V}_{ik}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{or}\{\begin{array}{l}{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_z}{A}^{ij}{C}_{44}(z_j)A_{km}^z{V}_{jm}+}{\displaystyle \sum _{j=1}^N{D}^{ij}{C}_{44}(z_j)W_{jk}=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1}}\\ {\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_z}{A}^{ij}{C}_{44}(z_j)A_{km}^z{V}_{jm}+}{\displaystyle \sum _{j=1}^N{D}^{ij}{C}_{44}(z_j)W_{jk}=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=N_Z}}\end{array}$$

(25)

$$\mathrm{Either}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{W}_{ik}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{or}\{\begin{array}{l}{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_z}{A}^{ij}{C}_{11}(z_j)A_{km}^z{W}_{jm}+}{\displaystyle \sum _{j=1}^N{D}^{ij}\hspace{0.17em}{C}_{12}(z_j)V_{jk}+{\displaystyle \sum _{j=1}^N{B}^{ij}{C}_{12}(z_j)U_{jk}}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1}}\\ {\displaystyle \sum _{j=1}^N{\displaystyle \sum _{m=1}^{N_z}{A}^{ij}{C}_{11}(z_j)A_{km}^z{W}_{jm}+}{\displaystyle \sum _{j=1}^N{D}^{ij}\hspace{0.17em}{C}_{12}(z_j)V_{jk}+{\displaystyle \sum _{j=1}^N{B}^{ij}{C}_{12}(z_j)U_{jk}}=0,}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=N_Z\end{array}\hspace{0.17em}\hspace{0.17em}$$

(26)

For a unified treatment of the problem, the degrees of freedom (DOFs) can be divided into the domain- and boundary-type DOFs, as follows:

$$\begin{array}{l}{U}_d=\left\{\begin{array}{l}{U}_{12}\\ U_{13}\\ ⋮\\ U_{N\left(N_z-1\right)}\end{array}\right\},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{V}_d=\left\{\begin{array}{l}{V}_{12}\\ V_{13}\\ ⋮\\ V_{N\left(N_z-1\right)}\end{array}\right\},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{W}_d=\left\{\begin{array}{l}{W}_{12}\\ W_{13}\\ ⋮\\ W_{N\left(N_z-1\right)}\end{array}\right\},\\ U_b=\left\{\begin{array}{l}{U}_{11}\\ U_{21}\\ ⋮\\ U_{N{N}_z}\end{array}\right\},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{V}_b=\left\{\begin{array}{l}{V}_{11}\\ V_{21}\\ ⋮\\ V_{N{N}_z}\end{array}\right\},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{W}_b=\left\{\begin{array}{l}{W}_{11}\\ W_{21}\\ ⋮\\ W_{N{N}_z}\end{array}\right\},\end{array}$$

(27)

In which $W_{mn}=W_m\left(z_n,t\right),V_{mn}=V_m\left(z_n,t\right),U_{mn}=U_m\left(z_n,t\right).$

This means that Equations (21)–(23) can be rearranged in matrix form as:

$$\left[\begin{array}{l}{k}_{bb}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{k}_{bd}\\ k_{db}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{k}_{dd}\end{array}\right]\left[\begin{array}{l}{d}_b\\ d_d\end{array}\right]=P\left[\begin{array}{l}\mathbf{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathbf{0}\\ \mathbf{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}G\end{array}\right]\left[\begin{array}{l}{d}_b\\ d_d\end{array}\right]$$

(28)

where the stiffness quantities $k_{bb},$$k_{db},$$k_{bd},$$k_{dd}$ refer to the boundary, $b$, and domain, $d$, weighting coefficients of the plate respectively, ${\left[d_b\hspace{0.17em}\hspace{0.17em}{d}_d\right]}^T$ is the displacement vector, $P$ refers to the buckling load and $G$ is the stability matrix due to the in-plane stresses.

At the same time, from Equations (23)–(26), the boundary weight coefficients can be replaced by the domain weight coefficients as follows:

$$d_b=k_{bb}{}^{-1}{k}_{bd}{d}_d$$

(29)

By substitution of $d_b$ from Equation (29) into Equation (28), and considering the harmonic solution $\left[\begin{array}{l}{U}_d\\ V_d\\ W_d\end{array}\right]=\left[\begin{array}{l}{\overline{U}}_d\\ {\overline{V}}_d\\ {\overline{W}}_d\end{array}\right]e^{iwt}$, the governing equations of the problem can be rewritten in terms of the domain unknowns, as follows:

$$K\left[\begin{array}{c}{\overline{U}}_d\\ {\overline{V}}_d\\ {\overline{W}}_d\end{array}\right]-PG\left[\begin{array}{c}{\overline{U}}_d\\ {\overline{V}}_d\\ {\overline{W}}_d\end{array}\right]=0$$

(30)

where $K$ is the stiffness matrix, defined as:

$$K=K_{dd}-K_{db}{K}_{bb}{}^{-1}{K}_{bd}$$

(31)

The solution of Equation (30) corresponds to the critical buckling load of the structure under in-plane conditions, labeled hereafter as $\lambda$.

In what follows, three different BCs for the buckling analysis of the plate structure are considered:

(i)

Simply supported BCs at all edges (SSSS):

$$\begin{array}{l}w\left(0,y,z\right)=w\left(a,y,z\right)=w\left(x,0,z\right)=w\left(x,b,z\right)=0\\ u\left(0,y,z\right)=u\left(a,y,z\right)=v\left(x,0,z\right)=v\left(x,b,z\right)=0\end{array}$$

(32)

(ii)

Clamped BCs at edges parallel to the y-axis (i.e., at $x=0,\hspace{0.17em}\hspace{0.17em}a$) and free BCs at edges parallel to the x-axis (i.e., at $y=0,\hspace{0.17em}\hspace{0.17em}b$) (CFCF):

$$u,v,w\left(0,y,z\right)=u,v,w\left(a,y,z\right)=0$$

(33)

(iii)

Clamped BCs at edges parallel to the x-axis (i.e., at $y=0,\hspace{0.17em}\hspace{0.17em}b$) and free BCs at edges parallel to the y-axis (i.e., at $x=0,\hspace{0.17em}\hspace{0.17em}a$) (FCFC):

$$u,v,w\left(x,0,z\right)=u,v,w\left(x,b,z\right)=0$$

(34)

Note that the (i)-type BC is considered for both uniaxial and biaxial loading conditions.

4. Numerical Investigation

The numerical study starts with a preliminary validation of the proposed formulation against the classical FE results (Ansys Workbench). At the present stage, the porosity coefficient is assumed to be zero $\left(e_0=0\right)$, together with a null Skempton coefficient, a null Biot’s modulus $\overline{M}=0,$ a null pore fluid pressure $p=0$ and ${\nu }_u=\nu =1/3$, $E_1=210 \mathrm{GPa}$. The rectangular plate selected here for the analysis has $b=1 \mathrm{m},$ $a/b=2$ and $h=0.1 \mathrm{m}$.

Table 1. Comparative evaluation of the first three buckling loads (10 GPa), as provided by our formulation and from FEM for different boundary conditions.

			FE-GDQ			FE		Difference (%)
	BC	${\mathit{\lambda }}_1$	${\mathit{\lambda }}_2$	${\mathit{\lambda }}_3$	${\mathit{\lambda }}_1$	${\mathit{\lambda }}_2$	${\mathit{\lambda }}_3$	${\mathit{\lambda }}_1$	${\mathit{\lambda }}_2$	${\mathit{\lambda }}_3$
Shear load	CFCF	0.314	0.529	0.624	0.317	0.574	0.678	0.94	7.83	7.96
Shear load	FCFC	1.405	1.406	2.168	1.430	1.430	2.117	1.74	1.67	2.35
Shear load	SSSS	2.534	2.642	3.062	2.575	2.691	3.0628	1.59	1.82	0
Axial load (x-direction)	SSSS	1.197	1.231	2.381	1.206	1.240	2.401	0.74	0.72	6.7
Biaxial load	SSSS	0.259	0.489	0.790	0.264	0.494	0.794	1.89	1.01	0.506

summarizes the results for the first three buckling loads $\left({\lambda }_1,\hspace{0.17em}\hspace{0.17em}{\lambda }_2,\hspace{0.17em}\hspace{0.17em}{\lambda }_3\right)$, as provided by our proposed FE-GDQ formulation and FEs, while considering the three different BCs (32)–(34) alternatively, as well as a uniaxial or biaxial loading condition. To model the problem in Ansys Workbench, the most accurate 20-node hexahedral quadratic elements were chosen to mesh the plate. First, a linear static analysis for edge loads equal to 1 Pa was performed in a static structural environment, and then the solutions were transferred to an eigenvalue buckling environment. Based on the results from this table, the good correspondence among predictions from the two alternative computational strategies proves the reliability and accuracy of the proposed FE-GDQ method to handle the problem. This is also confirmed in terms of mode shapes, as visible in the contour plots of Figure 2 and Figure 3, at least for the SSSS rectangular plate under a shear and biaxial loading, respectively.

After this validation step, the numerical study aimed at computing the buckling load of FG-saturated porous plates in undrained conditions, accounting for the effects of different BCs, aspect ratios, Skempton coefficients, porosity distributions and porosity coefficients, while keeping the geometry and material properties of the structure fixed. In

Table 2. First four shear buckling loads for a SSSS FG-saturated porous plate with different aspect ratios, $a/b,$ porosity distributions and porosity coefficients, $e_0$, keeping $\beta =0.6$ fixed.

Aspect Ratio	Buckling Load (10 GPa)	PNND			PNSD			PUD
		${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$	${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$	${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$
a/b = 1	${\lambda }_1$	2.666	1.935	1.226	2.674	1.961	1.263	2.269	1.212	0.275
	${\lambda }_2$	2.666	1.936	1.226	2.674	1.961	1.263	2.269	1.212	0.275
	${\lambda }_3$	2.958	2.255	1.442	2.998	2.265	1.390	2.595	1.482	0.354
	${\lambda }_4$	2.984	2.261	1.447	3.016	2.266	1.390	2.607	1.483	0.354
a/b = 2	${\lambda }_1$	2.281	1.863	1.162	2.392	2.005	1.269	2.083	1.281	0.2952
	${\lambda }_2$	2.339	1.865	1.162	2.438	2.029	1.271	2.105	1.282	0.2954
	${\lambda }_3$	2.558	1.968	1.237	2.627	2.068	1.379	2.248	1.334	0.3464
	${\lambda }_4$	2.719	2.044	1.307	2.761	2.094	1.388	2.350	1.356	0.3467

,

Table 3. First four shear buckling loads for a FCFC FG-saturated porous plate with different aspect ratios, $a/b,$ porosity distributions and porosity coefficients, $e_0$, keeping $\beta =0.6$ fixed.

Aspect Ratio	Buckling Load (10 GPa)	PNND			PNSD			PUD
		${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$	${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$	${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$
a/b = 1	${\lambda }_1$	1.143	0.897	0.536	1.199	0.978	0.692	1.020	0.636	0.177
	${\lambda }_2$	1.150	0.900	0.543	1.208	0.983	0.693	1.020	0.645	0.182
	${\lambda }_3$	1.653	1.267	0.773	1.705	1.337	0.845	1.455	0.865	0.219
	${\lambda }_4$	1.664	1.275	0.778	1.717	1.343	0.849	1.464	0.868	0.219
a/b = 2	${\lambda }_1$	1.216	0.954	0.576	1.276	1.039	0.737	1.082	0.680	0.191
	${\lambda }_2$	1.217	0.955	0.576	1.277	1.040	0.738	1.083	0.680	0.191
	${\lambda }_3$	1.768	1.360	0.834	1.824	1.401	0.912	1.558	0.931	0.238
	${\lambda }_4$	1.768	1.360	0.834	1.824	1.402	0.912	1.558	0.931	0.238

and

Table 4. First four shear buckling loads for a CFCF FG-saturated porous plate with different aspect ratios, $a/b,$ porosity distributions and porosity coefficients, $e_0$, keeping $\beta =0.6$ fixed.

Aspect Ratio	Buckling Load (10 GPa)	PNND			PNSD			PUD
		${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$	${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$	${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$
a/b = 1	${\lambda }_1$	1.143	0.897	0.536	1.199	1.006	0.692	1.020	0.636	0.177
	${\lambda }_2$	1.150	0.900	0.543	1.208	1.007	0.693	1.020	0.645	0.182
	${\lambda }_3$	1.653	1.267	0.773	1.705	1.434	0.845	1.455	0.865	0.219
	${\lambda }_4$	1.664	1.275	0.778	1.717	1.350	0.849	1.464	0.868	0.219
a/b = 2	${\lambda }_1$	0.310	0.255	0.153	0.330	0.292	0.241	0.287	0.186	0.055
	${\lambda }_2$	0.497	0.393	0.229	0.538	0.470	0.375	0.449	0.281	0.078
	${\lambda }_3$	0.598	0.466	0.281	0.610	0.520	0.405	0.528	0.325	0.096
	${\lambda }_4$	0.805	0.622	0.374	0.817	0.682	0.504	0.708	0.430	0.116

, the first four shear buckling loads are summarized for a FG-saturated plate under the SSSS, FCFC and CFCF BCs respectively, while keeping the Skempton coefficient constant. The systematic study starts by considering a square plate with aspect ratio $a/b=1$, thus extending the analysis to a rectangular plate with $a/b=2$. In each case, the porosity coefficient, $e_0$, is gradually increased from 0.3 up to 0.6, by steps of 0.3, to check for the sensitivity of the buckling response to porosity. As expected, when the porosity coefficient increases, the stiffness of the structure decreases, and the buckling load decreases as well. The results demonstrate that the maximum and minimum values of the buckling load are associated with the symmetric (PNSD) and uniform (PUD) porosity distributions respectively, due to the highest and lowest stiffness reached in the structure. An intermediate buckling load level, instead, is always obtained for a PNND porosity distribution within the material. It also seems that the effect of the porosity coefficient, $e_0$, on the buckling load becomes more pronounced for a uniform porosity distribution than the other ones. Moreover, by increasing the aspect ratio $a/b$, the buckling load can decrease or increase, under the same assumptions of porosity coefficient and distribution, depending on the selected BC. Differently from the FCFC case, a clear reduction of the buckling load is noticed for a SSSS and CFCF porous plate with an increased aspect ratio. This confirms the strict dependence of the stiffness on the geometrical dimensions and BCs of the structural member.

The same parametric investigation was then repeated for a biaxial loading (

Table 5. First four biaxial buckling loads for a SSSS FG-saturated porous plate with different aspect ratios, $a/b,$ porosity distributions and porosity coefficients, $e_0$, keeping $\beta =0.6$ fixed.

Aspect Ratio	Buckling Load (10 GPa)	PNND			PNSD			PUD
		${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$	${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$	${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$
a/b = 1	${\lambda }_1$	0.392	0.303	0.166	0.406	0.362	0.309	0.346	0.216	0.069
	${\lambda }_2$	0.929	0.723	0.408	0.949	0.830	0.658	0.827	0.523	0.151
	${\lambda }_3$	0.929	0.723	0.408	0.949	0.830	0.658	0.827	0.523	0.151
	${\lambda }_4$	0.141	0.108	0.604	0.144	0.127	0.926	1.245	0.763	0.210
a/b = 2	${\lambda }_1$	0.233	0.186	0.106	0.240	0.221	0.193	0.209	0.138	0.0420
	${\lambda }_2$	0.434	0.336	0.188	0.405	0.400	0.329	0.384	0.239	0.0680
	${\lambda }_3$	0.698	0.564	0.333	0.703	0.634	0.519	0.634	0.426	0.120
	${\lambda }_4$	0.778	0.603	0.341	0.804	0.702	0.570	0.688	0.427	0.131

) and a uniaxial loading acting in the x-direction of a SSSS structure (

Table 6. First four uniaxial buckling loads (in the x-direction) for a SSSS FG-saturated porous plate with different aspect ratios, $a/b,$ porosity distributions and porosity coefficients, $e_0$, keeping $\beta =0.6$ fixed.

Aspect Ratio	Buckling Load (10 GPa)	PNND			PNSD			PUD
		${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$	${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$	${\mathit{e}}_0=0.3$	${\mathit{e}}_0=0.6$	${\mathit{e}}_0=0.9$
a/b = 1	${\lambda }_1$	0.781	0.603	0.330	0.808	0.721	0.614	0.690	0.430	0.122
	${\lambda }_2$	0.981	0.774	0.446	0.983	0.875	0.703	0.879	0.568	0.168
	${\lambda }_3$	1.554	1.225	0.724	1.526	1.320	0.981	1.394	0.895	0.261
	${\lambda }_4$	2.147	1.688	1.017	2.081	1.856	1.218	1.925	1.225	0.350
a/b = 2	${\lambda }_1$	1.053	0.816	0.462	1.073	0.943	0.771	0.930	0.574	0.159
	${\lambda }_2$	1.101	0.860	0.495	1.117	0.982	0.795	0.977	0.613	0.173
	${\lambda }_3$	2.097	1.553	0.880	2.082	1.725	1.314	1.801	1.020	0.260
	${\lambda }_4$	2.331	1.771	1.042	2.293	1.907	1.405	2.035	1.204	0.319

). Based on a comparative evaluation of results in Table 5 and Table 6 with Table 2, a uniaxial or biaxial loading condition clearly reduces the buckling load of the structure under the same values of $e_0$, $a/b$ and porosity distribution.

As expected, this reduction is much more pronounced for rectangular plates subjected to a biaxial loading condition, due to the overall decay of the structural stiffness. A further investigation considered the effect of the porosity distribution and Skempton coefficient on the first four buckling loads of FG-saturated plates, as listed in

Table 7. First four shear buckling loads for a SSSS FG-saturated porous plate with different aspect ratios, $a/b,$ porosity distributions and Skempton coefficients, $\beta$, while keeping $e_0=0.6$ fixed.

Aspect Ratio	Buckling Load (10 GPa)	PNND			PNSD			PUD
		β = 0.0	β = 0.6	β = 0.9	β = 0.0	β = 0.6	β = 0.9	β = 0.0	β = 0.6	β = 0.9
a/b = 1	${\lambda }_1$	1.839	1.935	1.941	1.954	1.961	1.982	1.196	1.212	1.303
	${\lambda }_2$	1.913	1.944	1.951	1.957	1.961	1.986	1.196	1.212	1.303
	${\lambda }_3$	2.216	2.255	2.264	2.236	2.265	2.274	1.478	1.480	1.486
	${\lambda }_4$	2.240	2.261	2.273	2.250	2.266	2.275	1.478	1.483	1.486
a/b = 2	${\lambda }_1$	1.689	1.863	1.963	1.942	2.005	2.012	1.190	1.281	1.326
	${\lambda }_2$	1.738	1.865	1.964	1.953	2.029	2.043	1.203	1.283	1.326
	${\lambda }_3$	1.913	1.968	2.018	2.050	2.068	2.094	1.284	1.334	1.341
	${\lambda }_4$	2.038	2.044	2.048	2.130	2.139	2.147	1.343	1.356	1.386

,

Table 8. First four shear buckling loads for a FCFC FG-saturated porous plate with different aspect ratios, $a/b,$ porosity distributions and Skempton coefficients, $\beta$, while keeping $e_0=0.6$ fixed.

Aspect Ratio	Buckling Load (10 GPa)	PNND			PNSD			PUD
		β = 0.0	β = 0.6	β = 0.9	β = 0.0	β = 0.6	β = 0.9	β = 0.0	β = 0.6	β = 0.9
a/b = 1	${\lambda }_1$	0.848	0.897	0.955	0.969	0.978	0.989	0.583	0.636	0.709
	${\lambda }_2$	0.853	0.900	0.967	0.978	0.983	0.993	0.583	0.645	0.728
	${\lambda }_3$	1.231	1.267	1.299	1.329	1.337	1.346	0.831	0.865	0.878
	${\lambda }_4$	1.240	1.275	1.306	1.338	1.343	1.371	0.836	0.868	0.879
a/b = 2	${\lambda }_1$	0.903	0.954	1.021	1.031	1.039	1.042	0.618	0.680	0.765
	${\lambda }_2$	0.904	0.955	1.026	1.032	1.040	1.068	0.619	0.680	0.765
	${\lambda }_3$	1.318	1.360	1.399	1.423	1.401	1.414	0.890	0.931	0.955
	${\lambda }_4$	1.318	1.360	1.399	1.423	1.402	1.416	0.890	0.931	0.955

and

Table 9. First four shear buckling loads for a CFCF FG-saturated porous plate with different aspect ratios, $a/b,$ porosity distributions and Skempton coefficients, $\beta$, while keeping $e_0=0.6$ fixed.

Aspect Ratio	Buckling Load (10 GPa)	PNND			PNSD			PUD
		β = 0.0	β = 0.6	β = 0.9	β = 0.0	β = 0.6	β = 0.9	β = 0.0	β = 0.6	β = 0.9
a/b = 1	${\lambda }_1$	0.848	0.897	0.955	0.969	0.978	0.989	0.583	0.636	0.709
	${\lambda }_2$	0.853	0.900	0.967	0.978	0.983	0.993	0.583	0.645	0.728
	${\lambda }_3$	1.231	1.267	1.299	1.329	1.337	1.346	0.831	0.865	0.878
	${\lambda }_4$	1.240	1.275	1.306	1.338	1.343	1.371	0.836	0.868	0.879
a/b = 2	${\lambda }_1$	0.237	0.255	0.280	0.284	0.292	0.296	0.164	0.186	0.222
	${\lambda }_2$	0.367	0.393	0.425	0.494	0.470	0.482	0.254	0.281	0.314
	${\lambda }_3$	0.458	0.466	0.500	0.520	0.524	0.529	0.302	0.325	0.387
	${\lambda }_4$	0.600	0.622	0.615	0.682	0.689	0.798	0.404	0.430	0.465

(for a shear loading condition), in

Table 10. First four biaxial buckling loads for a SSSS FG-saturated porous plate with different aspect ratios, $a/b$, porosity distributions and Skempton coefficients, $\beta$, while keeping $e_0=0.6$ fixed.

Aspect Ratio	Buckling Load (10 GPa)	PNND			PNSD			PUD
		β = 0.0	β = 0.6	β = 0.9	β = 0.0	β = 0.6	β = 0.9	β = 0.0	β = 0.6	β = 0.9
a/b = 1	${\lambda }_1$	0.288	0.303	0.325	0.361	0.362	0.389	0.198	0.216	0.246
	${\lambda }_2$	0.683	0.723	0.781	0.825	0.830	0.846	0.472	0.523	0.603
	${\lambda }_3$	0.683	0.723	0.781	0.825	0.830	0.846	0.472	0.523	0.603
	${\lambda }_4$	0.103	0.108	0.146	0.123	0.127	0.131	0.711	0.763	0.843
a/b = 2	${\lambda }_1$	0.171	0.186	0.207	0.213	0.221	0.229	0.119	0.138	0.168
	${\lambda }_2$	0.319	0.369	0.362	0.397	0.400	0.405	0.219	0.239	0.273
	${\lambda }_3$	0.513	0.564	0.635	0.616	0.634	0.638	0.362	0.426	0.480
	${\lambda }_4$	0.573	0.603	0.644	0.697	0.702	0.713	0.393	0.427	0.526

(for an axial loading condition) and in

Table 11. First four uniaxial buckling loads (x-direction) for a SSSS FG-saturated porous plate with different aspect ratios, $a/b$, porosity distributions and Skempton coefficients, $\beta$, while keeping $e_0=0.6$ fixed.

Aspect Ratio	Buckling Load (10 GPa)	PNND			PNSD			PUD
		β = 0.0	β = 0.6	β = 0.9	β = 0.0	β = 0.6	β = 0.9	β = 0.0	β = 0.6	β = 0.9
a/b = 1	${\lambda }_1$	0.572	0.603	0.648	0.719	0.721	0.737	0.394	0.430	0.491
	${\lambda }_2$	0.722	0.774	0.848	0.867	0.875	0.884	0.502	0.568	0.672
	${\lambda }_3$	1.147	1.225	1.331	1.318	1.320	1.384	0.796	0.895	1.046
	${\lambda }_4$	1.589	1.688	1.817	1.756	1.856	1.878	1.100	1.225	1.403
a/b = 2	${\lambda }_1$	0.777	0.816	0.868	0.938	0.943	0.941	0.531	0.574	0.636
	${\lambda }_2$	0.813	0.860	0.923	0.977	0.982	0.993	0.558	0.613	0.693
	${\lambda }_3$	1.560	1.553	1.550	1.720	1.725	1.749	1.009	1.020	1.042
	${\lambda }_4$	1.735	1.771	1.817	1.901	1.907	1.923	1.163	1.204	1.276

(for a biaxial loading condition), while keeping the porosity coefficient fixed at $e_0=0.6$. The same BCs from Table 1 are accounted here for the analyses. Based on a comparative estimation of results from these tables, it is confirmed that the maximum and minimum buckling loads are always associated with a symmetric (PNSD) and uniform (PUD) porosity distribution, respectively.

Results denote that in drained conditions (i.e., $\beta =0$), the plate always features the smallest buckling load, under a fixed aspect ratio, porosity coefficient and distribution. An increasing value of the Skempton coefficient, instead, enables an increased buckling load, because of a decreased compressibility of the fluid within pores. In other words, if the compressibility of the pore fluid becomes high ($\beta \to 0$), the mechanical response of the plate resembles that of a porous plate in drained conditions (i.e., in the absence of fluid). In this condition, the structural stiffness reaches its minimum value along with the lowest buckling load. Differently, when the compressibility of a pore fluid becomes small ($\beta \to 1$), the plate behaves as a rigid solid, thus reaching its highest load magnitude. Furthermore, the effect of the Skempton coefficient on the buckling load, for a uniform distribution, seems to be more pronounced than other porosity distributions. By comparing Table 2, Table 3, Table 4, Table 5 and Table 6 with Table 7, Table 8, Table 9, Table 10 and Table 11, it is worth observing the higher sensitivity of the buckling response to the porosity coefficient than the Skempton coefficient. The first four buckling mode shapes are finally plotted in Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8, for a rectangular plate with $a=2\hspace{0.17em}\mathrm{m},\hspace{0.17em}\hspace{0.17em}b=1\hspace{0.17em}\mathrm{m}$, under different loading and boundary conditions, and a fixed value of $e_0=\beta =0.6$. More specifically, in Figure 4, Figure 5 and Figure 6, the rectangular plate is subjected to a shear loading condition, with a clear compatibility among the displacement field and the selected BCs (i.e., CFCF, FCFC and SSSS, respectively). Figure 7 and Figure 8 plot the four mode shapes for the same plate under a uniaxial (Figure 7) and biaxial (Figure 8) loading, where the kinematic response clearly changes depending on the selected loading condition.

5. Conclusions

The present work focused on the buckling behavior of FG-saturated porous rectangular plates subjected to normal and shear loads, while adopting Biot’s constitutive law and proposing a combined FE-GDQ as an efficient computational tool to solve the problem. This means that the in-plane problem has been discretized in the x–y-directions by means of classical FEs, and follows a weak formulation. Along the thickness direction (the z-direction), instead, the problem is defined in a strong form based on a GDQ approximation. This mixed method deals with a three-dimensional theory of elasticity without any additional kinematic assumption for the plate deformability. Various numerical examples have been considered and solved systematically to check for the reliability of the proposed method against a pure FE response, as well as to study the sensitivity of the response to some input parameters, i.e., the geometrical aspect ratio, the Skempton coefficients, the porosity distribution and coefficient and the BCs. Based on the parametric analysis, the main conclusions can be summarized as follows:

• The porosity coefficient more significantly affects the buckling load than the Skempton coefficient. In detail, an increased porosity coefficient and a decreased Skempton coefficient yield an overall decrease of the buckling load.
• Among different boundary and loading conditions, the maximum and minimum values of the buckling load are reached for a FCFC plate under a shear loading and a SSSS plate under a biaxial loading condition, respectively.
• The influence of the porosity coefficient on the buckling load for a uniform distribution is larger than other types of non-uniform porosity distributions.
• The effect of the Skempton coefficient on the buckling load, for a uniform distribution, is larger than other types of porosity distributions.
• By increasing the ratio, the buckling load generally decreases, except for a FCFC plate under a shear load and a SSSS plate under a normal uniaxial load, because of the variability in stiffness of the overall structure.
• The proposed method is verified to be a reliable tool for the computational study of saturated porous materials and structures, even from a design standpoint.