Size Dependent Buckling Analysis of a FG-CNTRC Microplate of Variable Thickness under Non-Uniform Biaxial Compression

Abstract

This paper combines third-order shear deformation theory (TSDT) and modified couple stress theory (MCST) with the principle of total potential energy to analyze the size-dependent buckling behavior of a functionally graded carbon nanotube-reinforced composite (FG-CNTRC) rectangular microplate of variable thickness subject to non-uniform biaxial compression when resting on an elastic medium. To determine the thickness qualities of the material, the extended rule of mixture was applied. In the context of microplate buckling in the presence of small length scale effects, the three kinds of Carbon Nanotube (CNT) distribution—(a) UD, (b) FG-O, and (c) FG-X—were used and compared. The equations governing various combinations of simply supported or clamped boundary conditions have been solved using the differential quadrature method (DQM). The correctness and precision of the solutions have been compared to another study. A numerical study was conducted to examine the dependence of buckling load on several parameters, including percentage change of thickness, length scale parameter, nonuniform edge loads, boundary conditions, volume percentage of the CNTs, CNT distribution, and elastic medium parameter. The results of their effects are presented in this paper.

1. Introduction

Both uniform and non-uniform plates are used in a wide variety of engineering contexts. Designers can accomplish this by using structures of varying thicknesses. Such plates are optimal when minimizing mass is crucial, as is the case with space-based construction. Compressive pressures cause plates to buckle, and knowing how devices behave under stress is crucial to making good design decisions [1]. Most buckling research investigates evenly distributed uniaxial and biaxial in-plane loads and is based on several plate theories [2,3,4,5,6]. Hamedani and Ranji [7] investigated the buckling analysis of stiffened plates exposed to varied biaxial edge compressions under different boundary conditions using super and conventional elements. Ruocco and Reddy [8] proposed an innovative discrete differential geometry-based computational method for buckling and vibration studies of non-uniformly thick rectangular plates. The bending, buckling, stress, and failure strength in inter-ply hybrid laminated composite have been investigated [9] using an isogeometric analysis based on non-uniform rational B-spline (NURBS) basis functions.

Structures of micron and submicron scales are widely used in many modern contexts, particularly in microelectromechanical systems (MEMS). Due to their small size, low power requirements, and reliable batch manufacture, these systems offer a wide range of possible engineering applications [10,11]. Lacking an intrinsic material length scale characteristic, classical continuum models are unable to account for size effects. Various size-dependent theories have been established in recent years, including the couple stress theory (CST) [12], the nonlocal elasticity theory [13,14], the strain gradient theory (SGT) [15,16], the modified strain gradient theory (MSGT) [17], the nonlocal strain gradient theory (NSGT) [18,19,20,21,22,23], and the modified couple stress theory (MCST) [24]. Among these size-dependent theories, MCST requires only a material length scale parameter, making it easy to implement. By eliminating the dilatation gradient tensor and the deviatoric stretch gradient tensor from SGT, as shown by MSCT, the size effect is shown to be a function of an internal length scale parameter. A variety of size-dependent beam and plate models, including the Timoshenko beam model [25,26], Kirchhoff plate model [27,28], Mindlin plate model [29,30,31], size-dependent Euler–Bernoulli beam model [32,33], and higher order shear deformable plate model [34,35] have recently been developed using MCST. Based on MCST in conjunction with TSDT, Thanh et al. [36] developed a numerical model to analyze the thermal bending and thermal buckling of composite laminate microplates while taking into account the effects of scale. Lou et al. [37] applied MSCT and first shear deformation theory (FSDT) to evaluate the vibrational dynamics of a buckled FG microplate supported by a nonlinear elastic medium. Zhang et al. [38] examined the mechanical buckling of a microplate consisting of two piezoelectric face sheets and one honeycomb structure considering MCST and sinusoidal shear deformation theory (SSDT). The effective mechanical characteristics of the honeycomb core were calculated using the Minghui and Jiuren relationships. Afshari and Adab [39] performed size-dependent mechanical buckling and free vibration calculations on simply supported microplates reinforced by FG graphene nanoplatelets using MCST SSDT. It was assumed that the composite microplate was made up of epoxy reinforced by GNPs, which were then arranged in a variety of symmetric and non-symmetric patterns along the thickness axis. By combining the Halpin–Tsai model and the rule of mixtures, they were able to make an educated assessment of the effective material qualities of the composite microplate.

Carbon nanotubes (CNTs) are a class of unique materials that exhibit remarkable electro-thermo-mechanical properties, and they have been receiving more and more attention in recent years [40,41,42,43]. CNTs are an attractive option for composite matrix reinforcement. The CNTs can be employed in a variety of technical applications when combined with a polymer composite material to create carbon nanotube reinforced composites (CNTRCs). Because of the uniform or random distribution of reinforcements within the system, the mechanical, physical, and thermal properties of conventional composites do not change spatially at the visible level. Functionally graded materials (FGMs) [44,45,46] have been created to improve this property. They have two or more elements with spatially variable characteristics and a non-uniform reinforcement distribution. Vaishali et al. [47] aligned between first-order shear deformation theory and higher-order zigzag theory using data and machine learning. This enabled uncertainty to be measured with a high degree of accuracy for laminated composite and sandwich plates. The reanalysis algorithm is a quick and precise solution that Li [48] applied to the vibration problem and incorporated to increase the effectiveness of metaheuristic optimization of composite laminates instead of a surrogate model. Savran and Aydin [49] used stochastic optimization on a laminated composite made of graphite, flax, and epoxy to obtain the highest fundamental frequency and lowest cost. Engineers can use this method to create materials better than homogenous ones. Shen [50] demonstrated that CNTRCs with nonuniform dispersion of CNTs through the media can be created by a powder metallurgy production process. Functionally graded carbon nanotube-reinforced composite is the name given to this type of reinforced composite media (FG-CNTRC). In recent years, more research has been conducted on these engineered materials in a variety of combinations to better understand their behavior in static and dynamic settings. The elastic instability and free vibration of FG-CNTRC annular sector plates on an elastic medium were presented by Zhang et al. [51] using FSDT and a generalized differential quadrature (GDQ) approach. Thai et al. [52] presented a size-dependent computational technique for static bending and free vibration studies of FG-CNTRC microplates based on MSGT and higher-order shear deformation theory. Cheshmeh et al. [53] investigated the buckling and vibration analysis of an FG-CNTRC rectangular plate under different temperature distributions and boundary conditions using the 12-uncertainty higher order shear deformation theory (HSDT). According to NSGT and isogeometric analysis (IGA), Phung-Van and Thai [54] predicted the natural frequencies of FG-CNTRC nanoplates employing HSDT.

Despite the extensive research on the uniaxial and biaxial buckling characteristics of composite microplates, the literature presented here shows that the buckling behavior of an FG-CNTRC microplate of varying thickness subjected to non-uniform biaxial compression resting on an elastic medium based on the TSDT in conjunction with MCST has not been studied theoretically or numerically. The main goal of the current research is to fill this gap by providing a novel combination of many existing methodologies. The three different types of CNT distribution—(a) UD, (b) FG-O, and (c) FG-X—are used and compared in the context of microplate buckling in the presence of small length scale effects. Total potential energy is used to obtain the governing equations, which are then solved by DQM for various combinations of simply supported or clamped boundary conditions.

2. Problem Description and Modeling

As can be seen in Figure 1, an FG-CNTRC microplate ($a\times b$) is subjected to biaxial non-uniform edge compression as its thickness varies in the y-direction (y-axis) while resting on an elastic medium with the stiffness coefficient $k_w$. Furthermore, $h_1$ indicates the thickness of the FG-CNTRC microplate’s left end, while $h_0$ represents the thickness of the right side of the FG-CNTRC microplate. In this illustration, UD represents the uniform distribution. FG-O and FG-X are two distinct FG distributions in the plate thickness direction.

2.1. MCST

Yang et al. [24] were the first to present MCST. Here, the strain energy density is analyzed as a function of the stress tensor, the curvature tensor, and the couple stress tensor. In an isotropic linear elastic body that occupies region $\mathsf{\Lambda }$, the strain energy ${\mathsf{\Pi }}_s$ may then be expressed as

$${\mathsf{\Pi }}_s=\frac{1}{2}{{\displaystyle \int }}_{\mathsf{\Lambda }}^{}\left(\sigma :\epsilon +m:\chi \right)\mathsf{\Lambda }$$

(1)

Here, $\epsilon$ signifies the strain tensor, $\epsilon$ represents the Cauchy stress tensor, $\chi$ represents the symmetric curvature tensor, and m represents the deviatoric component of the couple stress tensor. Specifically, these tensors are characterized by [24]

$$\begin{array}{c}\epsilon =\frac{1}{2}\left[\nabla \mathrm{u}+{\left(\nabla \mathrm{u}\right)}^{\mathrm{T}}\right],\chi =\frac{1}{2}\left[\nabla \mathsf{\phi }+{\left(\nabla \mathsf{\phi }\right)}^{\mathrm{T}}\right],\hfill \\ \sigma =\lambda \mathrm{tr}\left(\epsilon \right)\mathrm{I}+2\mu \epsilon ,m=2l^2\mu \chi ,\hfill \end{array}$$

(2)

where $l$ is a material length scale parameter that is considered to be a material feature describing the influence of the coupling stress, u is the displacement vector, $\lambda$ and $\mu$ are Lame’s constants, and $\mathsf{\phi }=\frac{1}{2}\mathrm{curl}\left(u\right)$ is the rotation vector. In this model, in addition to the standard Lame’s constants, only one length scale parameter $l$ is required. It should be noted that the additional parameter $l$ is technically the square of the modulus of curvature/modulus of the shear ratio and physically a property measuring the influence of couple stress.

2.2. Structural Model

Both shear deformation and rotating inertia in the transverse direction are taken into account by TSDT. Total lateral displacement (U,V,W) in Cartesian coordinates for microplate can be shown as [55,56]

$$\begin{array}{c}U\left(x,y,z\right)=u\left(x,y\right)+z{\theta }_x\left(x,y\right)-C_1{z}^3\left({\theta }_x\left(x,y\right)+\frac{\partial w\left(x,y\right)}{\partial x}\right),\hfill \\ V\left(x,y,z\right)=v\left(x,y\right)+z{\theta }_y\left(x,y\right)-C_1{z}^3\left({\theta }_y\left(x,y\right)+\frac{\partial w\left(x,y\right)}{\partial y}\right),\hfill \\ W\left(x,y,z\right)=w\left(x,y\right)\hfill \end{array}$$

(3)

in which $C_1=\frac{4}{3h^2\left(y\right),}$. In addition, $w$ is the transverse displacement component; ${\theta }_x$ and ${\theta }_y$ represent rotations of the middle plane around the x- and y-axes, respectively; and u and v are the membrane displacements along the x- and y-axes, respectively. However, using the linear strain-displacement connection, the specific components of the normal strains $\left({\epsilon }_{xx},{\epsilon }_{yy}\right)$ and shear strains $\left({\gamma }_{xz},{\gamma }_{yz},{\gamma }_{xy}\right)$ are stated as follows:

$$\begin{array}{c}{\epsilon }_{xx}=\frac{\partial u}{\partial x}+z\frac{\partial {\theta }_x}{\partial x}-c_1{z}^3\left(\frac{\partial {\theta }_x}{\partial x}+\frac{{\partial }^{2}w}{\partial x^2}\right),\hfill \\ {\epsilon }_{yy}=\frac{\partial v}{\partial y}+z\frac{\partial {\theta }_y}{\partial y}-c_1{z}^3\left(\frac{\partial {\theta }_y}{\partial y}+\frac{{\partial }^{2}w}{\partial y^2}\right),\hfill \\ {\gamma }_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}+z\left(\frac{\partial {\theta }_x}{\partial y}+\frac{\partial {\theta }_y}{\partial x}\right)-c_1{z}^3\left(\frac{\partial {\theta }_x}{\partial y}+\frac{\partial {\theta }_y}{\partial x}+2\frac{{\partial }^{2}w}{\partial x\partial y}\right),\hfill \\ {\gamma }_{xz}=\left(1-3c_1{z}^2\right)\left({\theta }_x+\frac{\partial w}{\partial x}\right),\hfill \\ {\gamma }_{yz}=\left(1-3c_1{z}^2\right)\left({\theta }_y+\frac{\partial w}{\partial y}\right)\hfill \end{array}$$

(4)

The rotation vector is calculated by substituting the results of Equations (1) and (2) into Equation (2) as

$$\begin{array}{c}{\phi }_x=\frac{1}{2}\left(\frac{\partial W}{\partial y}-\frac{\partial V}{\partial z}\right)=\frac{1}{2}\left(\left[\frac{\partial w}{\partial y}-{\theta }_y\right]+3C_1{z}^2\left[{\theta }_y+\frac{\partial w}{\partial y}\right]\right),\hfill \\ {\phi }_y=\frac{1}{2}\left(\frac{\partial U}{\partial z}-\frac{\partial W}{\partial x}\right)=-\frac{1}{2}\left(\left[\frac{\partial w}{\partial x}-{\theta }_x\right]+3C_1{z}^2\left[{\theta }_x+\frac{\partial w}{\partial x}\right]\right),\hfill \\ {\phi }_z=\frac{1}{2}\left(\frac{\partial V}{\partial x}-\frac{\partial U}{\partial y}\right)=\frac{1}{2}\left(\left[\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}\right]+\left(z-C_1{z}^3\right)\left[\frac{\partial {\theta }_y}{\partial x}-\frac{\partial {\theta }_x}{\partial y}\right]\right).\hfill \\ {\chi }_x=\frac{1}{2}\left(\frac{{\partial }^{2}w}{\partial x\partial y}-\frac{\partial {\theta }_y}{\partial x}\right)+\frac{3C_1{z}^2}{2}\left[\frac{\partial {\theta }_y}{\partial x}+\frac{{\partial }^{2}w}{\partial x\partial y}\right],\hfill \\ {\chi }_y=-\frac{1}{2}\left(\frac{{\partial }^{2}w}{\partial x\partial y}-\frac{\partial {\theta }_x}{\partial y}\right)-\frac{3C_1{z}^2}{2}\left(\frac{\partial {\theta }_x}{\partial y}+\frac{{\partial }^{2}w}{\partial x\partial y}\right),\hfill \\ {\chi }_z=\frac{1-3C_1{z}^2}{2}\left(\frac{\partial {\theta }_y}{\partial x}-\frac{\partial {\theta }_x}{\partial y}\right),\hfill \\ {\chi }_{xy}=\frac{1+3C_1{z}^2}{4}\left(\frac{{\partial }^{2}w}{\partial y^2}-\frac{{\partial }^{2}w}{\partial x^2}\right)+\frac{1-3C_1{z}^2}{4}\left(\frac{\partial {\theta }_x}{\partial x}-\frac{\partial {\theta }_y}{\partial y}\right),\hfill \\ {\chi }_{xz}=\frac{1}{4}\left(\frac{{\partial }^{2}v}{\partial x^2}-\frac{{\partial }^{2}u}{\partial x\partial y}\right)+\frac{z-C_1{z}^3}{4}\left(\frac{{\partial }^2{\theta }_y}{\partial x^2}-\frac{{\partial }^2{\theta }_x}{\partial x\partial y}\right)+\frac{3C_{1}z}{2}\left({\theta }_y+\frac{\partial \mathrm{w}}{\partial y}\right),\hfill \\ {\chi }_{yz}=\frac{1}{4}\left(-\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}v}{\partial x\partial y}\right)+\frac{z-C_1{z}^3}{4}\left(\frac{{\partial }^2{\theta }_y}{\partial x\partial y}-\frac{{\partial }^2{\theta }_x}{\partial y^2}\right)-\frac{3C_{1}z}{2}\left({\theta }_x+\frac{\partial \mathrm{w}}{\partial x}\right).\hfill \end{array}$$

(5)

Based on Equation (2), the strain–stress relationship for the FG-CNTRC microplate in the plain stress state can be written as:

$$\begin{array}{c}\left[\begin{array}{c}\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\\ {\sigma }_{yz}\\ {\sigma }_{xz}\end{array}\right]=\left\{\begin{array}{ccccc}{Q}_{11}& Q_{12}& 0& 0& 0\\ Q_{21}& Q_{22}& 0& 0& 0\\ 0& 0& Q_{66}& 0& 0\\ 0& 0& 0& Q_{44}& 0\\ 0& 0& 0& 0& Q_{55}\end{array}\right\}\left[\begin{array}{c}\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\gamma }_{xy}\end{array}\\ {\gamma }_{yz}\\ {\gamma }_{xz}\end{array}\right],\hfill \\ Q_{11}=\frac{E_{11}}{1-{\vartheta }_{12}{\vartheta }_{21}},Q_{22}=\frac{E_{22}}{1-{\vartheta }_{12}{\vartheta }_{21}},Q_{12}=\frac{{\vartheta }_{12}{E}_{22}}{1-{\vartheta }_{12}{\vartheta }_{21}},Q_{66}=G_{12},Q_{44}=G_{23},Q_{55}=G_{13}\hfill \end{array}$$

(6)

Here, $E_{11}$ and $E_{22}$ express the effective Young’s modulus of the FG-CNTRC microplate, ${\vartheta }_{12}$ and ${\vartheta }_{21}$ stand for Poisson’s ratio, while $G_{12}$, $G_{23}$, and $G_{13}$ denote the shear modulus. Applying the rule of mixture, we can characterize the material properties of the FG-CNTRC plate as [57]

$$\begin{array}{c}{E}_{11}={\zeta }_1{V}_{CNT}{E}_{11}^{CNT}+V_m{E}_m,\hfill \\ \frac{{\zeta }_2}{E_{22}}=\frac{V_{CNT}}{E_{22}^{CNT}}+\frac{V_m}{E_m},\hfill \\ \frac{{\zeta }_3}{G_{12}}=\frac{V_{CNT}}{G_{12}^{CNT}}+\frac{V_m}{G_m},\hfill \end{array}$$

(7)

where ($E_m$, $G_m$) represent the isotropic matrix moduli, and $\left(E_{11}^{CNT},E_{22}^{CNT},G_{12}^{CNT}\right)$ represent the material characteristics of CNTs. Furthermore, the CNT efficiency characteristics are denoted by ${\zeta }_1$, ${\zeta }_2$, and ${\zeta }_3$, and

$$\begin{array}{c}{V}_m+V_{CNT}=1,\hfill \\ {\vartheta }_{12}=V_{CNT}^{*}{\vartheta }_{12}^{CNT}+V_m{\vartheta }_m,\hfill \end{array}$$

(8)

where $V_{CNT}^{*}$ and $V_{CNT}$ denote the volume percentages of each FG-CNTRC layer, which can be expressed for three types of CNT distribution as UD, FG-O, and FG-X [57]:

$$V_{CNT}=\{\begin{array}{c}{V}_{CNT}^{*} \left(\mathrm{UD}\right)\\ \begin{array}{c}2\left(1-\frac{2\left|z\right|}{h_c}\right)V_{CNT}^{*}, \left(\mathrm{FG}-\mathrm{O}\right)\\ 4\frac{\left|z\right|}{h_c}{V}_{CNT}^{*}, \left(\mathrm{FG}-\mathrm{X}\right)\end{array}\end{array}$$

(9)

The total potential energy concept is applied to derive the FG-CNTRC microplate equilibrium equations as

$$\delta \left({\mathsf{\Pi }}_s+{\mathsf{\Pi }}_f\right)=0$$

(10)

where $\delta$ refers to a variation with regard to both x and y. The total energy of deformations is denoted by ${\mathsf{\Pi }}_s$, and the virtual work done by external forces is denoted by ${\mathsf{\Pi }}_f$. The variation of the total strain energy based on Equation (1) can be expressed as

$$\begin{array}{cc}\delta {\mathsf{\Pi }}_s={{\displaystyle \int }}_A^{}\hfill & {{\displaystyle \int }}_{-\frac{h\left(y\right)}{2}}^{\frac{h\left(y\right)}{2}}({\sigma }_{xx}\delta {\epsilon }_{xx}+{\sigma }_{yy}\delta {\epsilon }_{yy}+{\sigma }_{xy}\delta {\gamma }_{xy}+{\sigma }_{yz}\delta {\gamma }_{yz}+{\sigma }_{xz}\delta {\gamma }_{xz}+m_x\delta {\chi }_x+m_y\delta {\chi }_y\hfill \\ \hfill & +m_z\delta {\chi }_z+2m_{xy}\delta {\chi }_{xy}+2m_{xz}\delta {\chi }_{xz}+2m_{yz}\delta {\chi }_{yz})dzdA\hfill \\ \hfill & ={{\displaystyle \int }}_S^{}\{N_x\delta \frac{\partial u}{\partial x}+\left(M_x-C_1{P}_x\right)\delta \frac{\partial {\theta }_x}{\partial x}-C_1{P}_x\delta \frac{{\partial }^{2}w}{\partial x^2}+N_y\delta \frac{\partial v}{\partial y}\hfill \\ \hfill & +\left(M_y-C_1{P}_y\right)\delta \frac{\partial {\theta }_y}{\partial y}-C_1{P}_y\delta \frac{{\partial }^{2}w}{\partial y^2}+N_{xy}\left(\frac{\delta \partial u}{\partial y}+\frac{\delta \partial v}{\partial x}\right)\hfill \\ \hfill & +\left(M_{xy}-C_1{P}_{xy}\right)\left(\frac{\delta \partial {\theta }_x}{\partial y}+\frac{\delta \partial {\theta }_y}{\partial x}\right)-2C_1{P}_{xy}\frac{\delta {\partial }^{2}w}{\partial x\partial y}\hfill \\ \hfill & +\left(Q_x-3C_1{R}_x\right)\left(\delta {\theta }_x+\frac{\delta \partial w}{\partial x}\right)+\left(Q_y-3C_1{R}_y\right)\left(\delta {\theta }_y+\frac{\delta \partial w}{\partial y}\right)\hfill \\ \hfill & +\frac{X_x}{2}\left(\frac{\delta {\partial }^{2}w}{\partial x\partial y}-\frac{\delta \partial {\theta }_y}{\partial x}\right)+\frac{3C_1{Z}_x}{2}\left(\frac{\delta {\partial }^{2}w}{\partial x\partial y}+\frac{\delta \partial {\theta }_y}{\partial x}\right)-\frac{X_y}{2}\left(\frac{\delta {\partial }^{2}w}{\partial x\partial y}-\frac{\delta \partial {\theta }_x}{\partial y}\right)\hfill \\ \hfill & -\frac{3C_1{Z}_y}{2}\left(\frac{\delta {\partial }^{2}w}{\partial x\partial y}+\frac{\delta \partial {\theta }_x}{\partial y}\right)+\frac{X_z-3C_1{Z}_{xy}}{2}\left(\frac{\delta \partial {\theta }_y}{\partial x}-\frac{\delta \partial {\theta }_x}{\partial y}\right)\hfill \\ \hfill & +\frac{X_{xy}+3C_1{Z}_{xy}}{2}\left(\frac{\delta {\partial }^{2}w}{\partial y^2}-\frac{\delta {\partial }^{2}w}{\partial x^2}\right)+\frac{X_{xy}-3C_1{Z}_{xy}}{2}\left(\frac{\delta \partial {\theta }_x}{\partial x}-\frac{\delta \partial {\theta }_y}{\partial y}\right)\hfill \\ \hfill & +\frac{X_{xz}}{2}\left(\frac{\delta {\partial }^{2}v}{\partial x^2}-\frac{\delta {\partial }^{2}u}{\partial x\partial y}\right)+\frac{Y_{xz}-C_1{W}_{xz}}{2}\left(\frac{\delta {\partial }^2{\theta }_y}{\partial x^2}-\frac{\delta {\partial }^2{\theta }_x}{\partial x\partial y}\right)\hfill \\ \hfill & +3C_1{Y}_{xz}\left(\delta {\theta }_y+\frac{\delta \partial w}{\partial y}\right)+\frac{X_{yz}}{2}\left(\frac{\delta {\partial }^{2}v}{\partial x\partial y}-\frac{\delta {\partial }^{2}u}{\partial y^2}\right)\hfill \\ \hfill & +\frac{Y_{yz}-C_1{W}_{yz}}{2}\left(\frac{\delta {\partial }^2{\theta }_y}{\partial x\partial y}-\frac{\delta {\partial }^2{\theta }_x}{\partial y^2}\right)-3C_1{Y}_{yz}\left(\delta {\theta }_x+\frac{\delta \partial w}{\partial x}\right)\}dA\hfill \end{array}$$

(11)

where $A$ is the of cross-sectional area, and

$$\begin{array}{c}\left\{N_i,M_i,P_i\right\}={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}\left\{1,z,z^3\right\}{\sigma }_{i}dz, i=x,y,xy,\hfill \\ \left\{Q_i,R_i\right\}={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}\left\{1,z^2\right\}{\sigma }_{iz}dz, i=x,y,\hfill \\ \left\{X_i,Y_i,Z_i,W_i\right\}={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}\left\{1,z,z^2,z^3\right\}{m}_{i}dz,\hfill \end{array}$$

(12)

The potential energy variation caused by non-uniform forces and the elastic medium can be expressed as

$$\delta {\mathsf{\Pi }}_f={\displaystyle \int }\left(-f_x\left(\frac{{\partial }^{2}w}{\partial x^2}\right)-f_y\left(\frac{{\partial }^{2}w}{\partial y^2}\right)-2f_{xy}\left(\frac{{\partial }^{2}w}{\partial x\partial y}\right)+K_{w}w\right)\delta wdA$$

(13)

in which the compressive forces in the x and y directions are denoted by $f_x$ and $f_y$, respectively, while the shear force is denoted by $f_{xy}$. The final form of the governing equations is determined by inserting Equations (11) and (13) into Equation (10) and performing various manipulations while keeping the coefficients $\delta u, \delta v, \delta w, \delta {\theta }_x$, and $\delta {\theta }_y$ equal to zero.

$$\begin{array}{c}\delta u: \frac{\partial N_x}{\partial x}+\frac{\partial N_{xy}}{\partial y}+\frac{1}{2}\left(\frac{{\partial }^2 X_{xz}}{\partial x\partial y}+\frac{{\partial }^2 X_{yz}}{\partial y^2}\right)=0,\hfill \\ \delta v: \frac{\partial N_y}{\partial y}+\frac{\partial N_{xy}}{\partial x}-\frac{1}{2}\left(\frac{{\partial }^2 X_{xz}}{\partial x^2}+\frac{{\partial }^2 X_{yz}}{\partial x\partial y}\right)=0,\hfill \\ \delta w_0: \frac{\partial Q_x}{\partial x}-3C_1\frac{\partial R_x}{\partial x}+\frac{\partial Q_y}{\partial y}-3C_1\frac{\partial R_y}{\partial y}+C_1\left(\frac{{\partial }^2{P}_x}{\partial x^2}+2\frac{{\partial }^2{P}_{xy}}{\partial x\partial y}+\frac{{\partial }^2{P}_y}{\partial y^2}\right)+3C_1\left(\frac{\partial Y_{xz}}{\partial y}-\frac{\partial Y_{yz}}{\partial x}\right)+\frac{1}{2}(\frac{{\partial }^2{X}_{xy}}{\partial x^2}+\hfill \\ 3C_1\frac{{\partial }^2{Z}_{xy}}{\partial x^2}-\frac{{\partial }^2{X}_{xy}}{\partial y^2}-3C_1\frac{{\partial }^2{Z}_{xy}}{\partial y^2}+\frac{{\partial }^2{X}_y}{\partial x\partial y}+3C_1\frac{{\partial }^2{Z}_y}{\partial x\partial y}-\frac{{\partial }^2{X}_x}{\partial x\partial y}-3C_1\frac{{\partial }^2{Z}_x}{\partial x\partial y})-f_x\left(\frac{{\partial }^{2}w}{\partial x^2}\right)-f_y\left(\frac{{\partial }^{2}w}{\partial y^2}\right)-\hfill \\ 2f_{xy}\left(\frac{{\partial }^{2}w}{\partial x\partial y}\right)+K_{w}w=0,\hfill \\ \delta {\theta }_x: \frac{\partial M_x}{\partial x}-C_1\frac{\partial P_x}{\partial x}+\frac{\partial M_{xy}}{\partial y}-C_1\frac{\partial P_{xy}}{\partial y}-\left(Q_x-3C_1{R}_x\right)+\frac{1}{2}(\frac{\partial X_{xy}}{\partial x}-3C_1\frac{\partial Z_{xy}}{\partial x}+\frac{\partial X_y}{\partial y}-3C_1\frac{\partial Z_y}{\partial y}-\frac{\partial X_z}{\partial y}+\hfill \\ 3C_1\frac{\partial Z_z}{\partial y}+\frac{{\partial }^2{Y}_{xz}}{\partial x\partial y}-3C_1\frac{{\partial }^2{W}_{xz}}{\partial x\partial y}+\frac{{\partial }^2{Y}_{yz}}{\partial y^2}-3C_1\frac{{\partial }^2{W}_{yz}}{\partial y^2})=0,\hfill \\ \delta {\theta }_y: \frac{\partial M_{xy}}{\partial x}-C_1\frac{\partial P_{xy}}{\partial x}+\frac{\partial M_y}{\partial y}-C_1\frac{\partial P_y}{\partial y}-\left(Q_y-3C_1{R}_y\right)-\frac{1}{2}(\frac{\partial X_x}{\partial x}-3C_1\frac{\partial Z_x}{\partial x}+\frac{\partial X_{xy}}{\partial y}-3C_1\frac{\partial Z_{xy}}{\partial y}-\frac{\partial X_z}{\partial x}+3C_1\frac{\partial Z_z}{\partial x}+\frac{{\partial }^2{Y}_{xz}}{\partial x^2}-\hfill \\ 3C_1\frac{{\partial }^2{W}_{xz}}{\partial x^2}+\frac{{\partial }^2{Y}_{yz}}{\partial x\partial y}-3C_1\frac{{\partial }^2{W}_{yz}}{\partial x\partial y})=0,\hfill \end{array}$$

(14)

Finally, the governing equations in terms of the FG-CNTRC microplate displacement are obtained and presented in Appendix A by substituting Equation (12) into Equation (14).

The corresponding boundary conditions for clamped and simply supported boundary conditions are as follows:

Clamped (C) edge:

$$\begin{array}{c}u=v=w={\theta }_x={\theta }_y=\frac{\partial w}{\partial x}=0 \to x=0,a,\hfill \\ u=v=w={\theta }_x={\theta }_y=\frac{\partial w}{\partial y}=0 \to y=0,b,\hfill \end{array}$$

(15)

Simply supported (S) edge:

$$\begin{array}{c}v=w={\theta }_y=N_x=M_x-C_1{P}_x+\frac{1}{2}\left(X_{xy}-3C_1{Z}_{xy}\right)=P_x=0 \to x=0,a,\hfill \\ u=w={\theta }_x=N_y=M_y-C_1{P}_y-\frac{1}{2}\left(X_{xy}-3C_1{Z}_{xy}\right)=P_y=0 \to y=0,b,\hfill \end{array}$$

(16)

3. Solution Method

DQM [57,58] is used to discretize the governing Equations (A1)–(A5) and their corresponding boundary conditions (15)–(16) to solve them. By using suitable weighting coefficients, this technique transforms the differential equations into an algebraic equation of the first order. The weighting coefficients are independent of the nature of the problem and only depend on the grid resolution. In other words, at a discontinuous point in a defined domain ($0<x<a$ and $0<y<b$), the partial derivatives of a function (let’s call it $w$) are approximated as a set of linear weighting coefficients. However, partial derivatives of the function $f\left(x_i,y_j\right)$ at the position $\left(x_i,y_j\right)$ are written as in the DQM [58,59]:

$$\begin{array}{c}{f}_x^n\left(x_i,y_j\right)={\displaystyle \sum }_{K=1}^N{g}_{iK}^{\left(n\right)}f\left(x_K,y_j\right)\left\{\begin{array}{c}n=1,2,\dots ,N-1\\ i=1,2,\dots ,N\end{array}\right\},\hfill \\ f_y^m\left(x_i,y_j\right)={\displaystyle \sum }_{l=1}^M{d}_{jl}^{\left(m\right)}f\left(x_i,y_j\right)\left\{\begin{array}{c}m=1,2,\dots ,M-1\\ j=1,2,\dots ,M\end{array}\right\},\hfill \\ f_{xy}^{\left(n+m\right)}\left(x_i,y_j\right)={\displaystyle \sum }_{K=1}^N{\displaystyle \sum }_{l=1}^M{g}_{iK}^{\left(n\right)}{d}_{jl}^{\left(m\right)}f\left(x_K,y_l\right)\left\{\begin{array}{c}n+m=2,..,N-1\\ i=1,2,\dots ,N\\ j=1,2,\dots ,M\end{array}\right\}\hfill \end{array}$$

(17)

in which $N$ and $M$ denote the number of grid points in the x and y directions, respectively. The weighting coefficients for the $n$ th-order derivative in the x and y directions are denoted by $g_{ij}^{\left(n\right)}$ and $d_{ij}^{\left(n\right)}$, respectively, which can be defined as follows:

For the first-order derivative, $n=1$:

$$g_{iK}^1=\frac{M\left(X_i\right)}{\left(X_i-X_K\right)M\left(X_K\right)}i,K=1,2,\dots ,N, K\ne i$$

(18)

where

$$M\left(X_i\right)={\displaystyle \prod }_{\begin{array}{c}k=1\\ k\ne i\end{array}}^N\left(X_i-X_K\right)$$

(19)

The weighting coefficients for larger derivatives ($n>1$) can be calculated using the following recursion relationships

$$g_{iK}^{\left(n\right)}=n\left(g_{ii}^{\left(n-1\right)}{g}_{iK}^{\left(1\right)}-\frac{g_{iK}^{\left(n-1\right)}}{X_i-X_K}\right)\left\{\begin{array}{c}i,K=1,2,\dots ,NK\ne i\\ n=2,3,\dots ,N-1\end{array}\right\}$$

(20)

and

$$g_{ii}^{\left(n\right)}=-{{\displaystyle \sum }}_{K=1}^N{g}_{iK}^{\left(n\right)}\left\{\begin{array}{c}i=1,2,\dots ,NK\ne i\\ n=1,2,\dots ,N-1\end{array}\right\}$$

(21)

The DQM can be implemented with any number of discrete grid points and any distribution of grid points. It has been demonstrated, however, that the grid point distribution, which is based on the widely established Gauss–Chebyshev–Lobatto points [58,59], produces findings that are sufficiently precise. The coordinates of the grid points are as follows based on the Gauss–Chebyshev–Lobatto grid point distribution:

$$x_i=\frac{a}{2}\left\{1-\cos \left[\frac{\left(i-1\right)\pi }{\left(N-1\right)}\right]\right\},y_j=\frac{b}{2}\left\{1-\cos \left[\frac{\left(j-1\right)\pi }{\left(M-1\right)}\right]\right\}$$

(22)

By inserting Equation (17) into the governing Equations (A1)–(A5) and their corresponding boundary conditions (15)–(16) and making various simplifications, the following matrix form is obtained as

$$\left(\left[K_{tot}\right]-F^{*}\left[I\right]\right)\left\{\mathsf{\Gamma }\right\}=0$$

(23)

in which $I$ signifies the identity matrix and $F^{*}$ denotes the buckling load. A classic eigenvalue solver can solve Equation (23). The number of discrete mesh points employed in the current DQM (i.e., $N=M$) is the same in both the $x$ and $y$ directions.

4. Results and Discussion

4.1. Model Validation and Convergence Checking

The validation and convergence analysis of the created model are covered here. Note that the physical and geometric criteria used in this section and subsequent sections are as follows:

$$E_m=2.1GPa, E_{11}^{CNT}=5.6466TPa, E_{22}^{CNT}=7.08TPa, G_{12}^{CNT}=1.9445TPa, {\vartheta }_{12}^{CNT}=0.175, {\vartheta }_m=0.24$$

(24)

In addition,

Table 1. Efficiency coefficients for different CNT volume fractions.

${\mathit{V}}_{\mathit{C}\mathit{N}\mathit{T}}^{*}$	Efficiency Parameters
	${\mathit{\zeta }}_1$	${\mathit{\zeta }}_2$	${\mathit{\zeta }}_3$
0.11	0.149	0.934	0.934
0.14	0.15	0.941	0.941
0.17	0.14	1.381	1.381

lists the CNT efficiency metrics for different volume fractions [57]. After eliminating the length scale parameter, the current model is used to calculate the non-dimensional buckling load ${\overline{N}}_{cr}=\frac{F^{*}{a}^2}{D_{11}}$ of an FG-CNTRC plate subjected to uniform uniaxial compressive pressure in fully clamped and simply supported boundary conditions, and the results are compared with those from [60] in

Table 2. Convergence test and comparison of dimensionless buckling load of an FG-CNTRC plate subjected to uniform uniaxial compressive pressure.

B.C	Types	Element Mesh					Ref. [60]
		$5\times 5$	$7\times 7$	$9\times 9$	$11\times 11$	$13\times 13$
SSSS	UD	13.5834	14.1105	14.1091	14.1078	14.1078	14.1073
	FG-O	8.2569	9.6432	9.8357	9.8314	9.8314	9.8306
	FG-X	15.8670	16.8605	17.0809	17.0642	17.0642	17.0631
CCCC	UD	24.2967	25.4071	25.6825	25.7338	25.7338	25.7329
	FG-O	19.6568	20.7710	21.0564	21.1230	21.1230	21.1221
	FG-X	25.9276	27.5414	27.8009	27.8893	27.8893	27.8882

. As the DQM procedure’s findings are dependent on the number of grid points, a convergence test is also conducted in Table 2. Good agreement can be seen between the predicted buckling values and the data reported in the later reference.

In another example,

Table 3. Comparison of the non-dimensional critical buckling loads of an elastic plate with variable thickness subjected to uniform uniaxial compressive pressure.

$\frac{\mathit{a}}{\mathit{b}}$	$\mathit{\alpha }$	Present	Ref. [61]	Ref. [62]
0.5	0.125	7.45	7.46	7.46
	0.5	11.56	11.61	11.61
1.5	0.125	5.19	5.19	5.19
	0.5	8.22	8.23	8.23

compares the non-dimensional critical buckling loads ${\overline{N}}_{cr}=\frac{F^{*}{b}^2}{{\pi }^2{D}_{11}}$ of thin, isotropic plates subjected to uniform uniaxial compressive pressure under simply supported boundary conditions found with the current solution to those reported by [61,62] for different aspect ratios. The plate was compressed uniaxially in the x direction. Using the law $h\left(y\right)=h_0\left(1+\alpha y\right), \alpha =\frac{h_1-h_0}{h_0}$, the thickness varies linearly along the $y$ axis. For different aspect ratios, our results show excellent agreement with those reported in the literature. Please refer to the cited references for further information on the effect of the aspect ratio on buckling load performance.

Finally, as seen in Figure 2, using the current approach, the buckling load for various values of the length scale parameter have been computed and compared with the details available in [63] for a square microplate with uniform thickness subjected to uniform uniaxial compressive pressure. For a number of length scale parameters, there is very good agreement between our results and those of [63], which confirms the findings of our approach.

4.2. Parametric Study

In this section, numerical examples are presented to examine the influences of the percentage change of thickness, the material length scale parameter to thickness ratio, nonuniform edge loads, boundary conditions, volume percentage of the CNTs, CNT distribution, and elastic medium parameter on the dimensionless buckling load ${\overline{N}}_{cr}=\frac{F^{*}{a}^2}{D_0}, D_0=\frac{E_m{h}_0{}^3}{12\left(1-{\vartheta }_m{}^2\right)}$ of the FG-CNTRC microplate. The thickness of the microplate is expected to vary in the $y$ direction according to the law $h=h_0\left(1+\alpha y^{\beta }\right)$, where $\beta$ indicates a non-uniform characteristic for a microplate with variable thickness. Depending on the value of the non-uniform parameter, the y-axis thickness will display a range of patterns. If $\beta =0$, for instance, the plate’s height remains the same in both the x and y axes. If we set $\beta =1$, we see a linear relationship between thickness ($h_0$) at $y=0$ and thickness ($h_1$) at $y=b$. At the value of two for the non-uniform parameter ($\beta =2$), the thickness ($h$) varies parabolically along the y axis.

Figure 3, Figure 4 and Figure 5 illustrate the variations of the dimensionless biaxial buckling load of the UD-CNTRC microplate against the material length scale parameter to thickness ratio $\left(\frac{l}{h_0}\right)$ for different boundary conditions (SSSS, CCSS, CCCC) and different values of the elastic medium stiffness parameter under uniform (Figure 3), sinusoidal (Figure 4), and reverse sinusoidal (Figure 5) edge compressive loads. It is assumed that $h_0=\frac{a}{10},\mathrm{and} a=b, V_{CNT}^{*}=0.11$. The boundary conditions have a significant impact on the variability of the dimensionless biaxial buckling load. Additionally, as anticipated, the buckling load for the CCCC boundary condition has the largest value, which results from the high stiffness effect of this kind of boundary condition on the structure. On the other hand, the structure has the highest buckling load when subjected to reverse sinusoidal load and the lowest buckling load when subjected to sinusoidal load. This may be because, in the case of inverse sinusoidal load, most of the load is distributed towards the edges of the FG-CNTRC microplate, where the microplate’s stiffness is substantially greater. As demonstrated by these figures, the critical buckling load increases as the material length scale parameter to thickness ratio increases. Since the couple stress effect adds bending stiffness that grows with $\left(\frac{l}{h_0}\right)$, the total bending stiffness and critical buckling load also grow with $\left(\frac{l}{h_0}\right)$. It can be shown that the dimensionless buckling load decreases with a rise in the non-uniform parameter ($\beta$). This is because the average stiffness of the microplate decreases as the non-uniform parameter increases. The buckling load is found to be proportional to the thickness difference between the microplate’s ends and the elastic medium stiffness parameter. As can be observed, the influence of the non-uniform parameter on buckling reduces as the material length scale parameter to thickness ratio rises; however, the influence of the thickness difference increases.

Figure 6, Figure 7 and Figure 8 present the effects of the CNT distributions on the uniaxial and biaxial buckling loads of the FG-CNTRC microplate against the material length scale parameter to thickness ratio $\left(\frac{l}{h_0}\right)$ for different boundary conditions (SSSS, CCSS, CCCC) under uniform (Figure 6), sinusoidal (Figure 7), and reverse sinusoidal (Figure 8) edge compressive loads. It is assumed that $h_0=\frac{a}{10},a=b,V_{CNT}^{*}=0.11,k_w=0,\alpha =0$. These figures show that the critical buckling load for biaxial compression is significantly lower than for uniaxial compression because stiffness is reduced when loads are applied in both the horizontal and vertical planes. Based on these findings, FG-X microplates have higher buckling load values than UD microplates for all possible CNT distributions, while FG-O microplates have lower values than UD microplates. The physical causes of this behavior are as follows: For the FG-X type, the number of CNTs on the bottom and top of the FG-CNTRC microplate is made as high as possible. This makes the structure of the system more rigid and strong. The results indicate that FG-X > UD > FG-O is the stiffness order for the FG-CNTRC microplates.

Table 4. The effects of the CNT distributions, CNT volume fractions, and length to thickness ratio on the variations of the dimensionless biaxial buckling load of the FG-CNTRC microplate under uniform and nonuniform edge loads (SSSS).

Load Cases	$V_{CNT}^{*}$	$a/h_0$	UD	FG-O	FG-X
Uniform load	0.11	10	76.3982	62.4456	85.6194
		20	108.6706	78.6910	134.0614
	0.14	10	82.8363	67.4566	92.5827
		20	123.6956	87.7549	147.3007
Sine load	0.11	10	59.4015	48.5749	66.5321
		20	84.6411	61.3355	102.4814
	0.14	10	64.3922	52.4675	71.9237
		20	96.2908	68.3911	112.2957
Inverse sine load	0.11	10	109.2032	124.1952	132.7104
		20	192.4701	142.7652	229.9330
	0.14	10	137.7109	116.9394	146.7540
		20	227.1229	164.4443	269.8629

,

Table 5. The effects of the CNT distributions, CNT volume fractions, and length to thickness ratio on the variations of the dimensionless biaxial buckling load of the FG-CNTRC microplate under uniform and nonuniform edge loads (CCSS).

Load Cases	$V_{CNT}^{*}$	$a/h_0$	UD	FG-O	FG-X
Uniform load	0.11	10	92.7860	79.2604	99.2072
		20	149.6673	116.2542	170.1263
	0.14	10	98.7655	84.5839	106.6632
		20	163.0214	128.1712	185.1515
Sine load	0.11	10	70.7530	62.8014	76.2244
		20	113.5681	89.0498	128.6020
	0.14	10	75.2025	64.7786	81.0391
		20	123.3760	97.9230	139.7811
Inverse sine load	0.11	10	134.3347	118.1731	141.0803
		20	248.9409	196.2968	281.2835
	0.14	10	141.6780	126.8329	148.2357
		20	279.8703	224.8863	310.8059

and

Table 6. The effects of the CNT distributions, CNT volume fractions, and length to thickness ratio on the variations of the dimensionless biaxial buckling load of the FG-CNTRC microplate under uniform and nonuniform edge loads (CCCC).

Load Cases	$V_{CNT}^{*}$	$a/h_0$	UD	FG-O	FG-X
Uniform load	0.11	10	116.5609	104.5692	124.5312
		20	200.8065	167.8815	221.3473
	0.14	10	123.0866	108.5070	130.5892
		20	215.0219	181.3804	237.6610
Sine load	0.11	10	84.9358	75.2791	89.8126
		20	145.5112	121.6703	160.8139
	0.14	10	89.5611	79.5182	94.8094
		20	156.0347	130.9422	172.8222
Inverse sine load	0.11	10	160.4343	148.0554	164.2181
		20	403.4661	370.3635	450.0608
	0.14	10	194.8202	168.3681	202.8523
		20	447.2782	391.7958	476.9338

show the effects of the CNT distributions, CNT volume fractions, and length to thickness ratio on the variations of the dimensionless biaxial buckling load of the FG-CNTRC microplate under uniform and nonuniform edge loads for different boundary conditions when $\frac{l}{h_0}=0.1,a=b,k_w=0,\alpha =0$. Thus, the buckling load of the FG-CNTRC microplates has a higher value when the volume fraction of CNT is greater. This is because the stiffness of the CNTRC microplates increases when the value of the CNT volume fraction rises. Furthermore, as shown in these tables, increasing the length to thickness ratio increases the total stiffness of the structure and, as a result, the dimensionless buckling load. Based on the tables, the effects of CNT distributions on buckling changes are amplified when the microplate is subjected to reverse sinusoidal loading.

5. Conclusions

This paper reports our examination of the small-scale effect on the stability of a rectangular FG-CNTRC microplate with varying thickness under non-uniform biaxial compression resting on an elastic medium using MCST in combination with TSDT. The extended rule of mixture is used to evaluate the characteristics of the material in the thickness. The DQM was used to solve the derived equations for various combinations of simply supported or clamped boundary conditions. By comparing the established model to data from previous studies, its accuracy and correctness were assessed. A comprehensive numerical analysis was performed to investigate the impact of thickness variations, length scale parameter, nonuniform edge loads, boundary conditions, CNT volume fractions, CNT distributions, and elastic medium parameters on the buckling load. From our findings, we can draw the following conclusions:

The structure has the highest buckling load when it is subjected to reverse sinusoidal load and the lowest buckling load when it is subjected to sinusoidal load. This may be because, in the case of inverse sinusoidal load, most of the load is distributed towards the edges of the FG-CNTRC microplate, where the microplate’s stiffness is substantially greater.

The dimensionless buckling load decreases with a rise in the non-uniform parameter ($\beta$). This is because the average stiffness of the microplate decreases as the non-uniform parameter increases.

The buckling load is found to be proportional to the thickness difference between the microplate’s ends.

FG-X microplates have higher buckling load values than UD microplates for all possible CNT distributions, while FG-O microplates have lower values than UD microplates.

The effects of CNT distributions on buckling changes are amplified when the microplate is subjected to reverse sinusoidal loading.