Dynamic Instability of Functionally Graded Graphene Platelet-Reinforced Porous Beams on an Elastic Foundation in a Thermal Environment

Abstract

Under thermal environment and axial forces, the dynamic instability of functionally graded graphene platelet (GPLs)-reinforced porous beams on an elastic foundation is investigated. Three modes of porosity distributions and GPL patterns are considered. The governing equations are given by the Hamilton principle. On the basis of the differential quadrature method (DQM), the governing equations are changed into Mathieu–Hill equations, and the main unstable regions of the porous composite beams are studied by the Bolotin method. Thermal buckling and thermo-mechanical vibration problems are also studied. The effects of porosity coefficients and GPL weight fraction, dispersion pattern, initial thermal loading, slenderness ratio, geometry and size, boundary conditions, and foundation stiffness are discussed. The conclusions show that an elastic foundation has an obvious enhancement effect on thermal buckling, free vibration, and dynamic instability, which improves the stiffness of the beam.

1. Introduction

Many natural porous materials have been widely used for thousands of years. Compared with continuous medium materials, porous materials have excellent impact resistance, electrical conductivity, energy absorption, and thermal management properties [1,2,3,4,5]. Porous materials are often used in biological tissue, sound insulation materials, and new photoelectric elements [6]. As a kind of porous material, metal foam has high strength and stiffness [7]. At present, many scholars have probed into the mechanical behavior of porous materials and the influence of various factors on the materials.

Nguyen et al. [8] investigated the bucking, bending, and vibration of functionally graded porous (FGP) beams by the Ritz method. Akbaş [9,10,11] examined bucking and vibration by the finite element method. By the differential transform method (DTM), Ebrahimi and Mokhtari [12] presented the vibration of rotating FGP beams. Rjoub and Hamad [13] reported on the vibration of FGP beams by the Transfer Matrix Method. Wattanasakulpong and Chaikittiratana [14] found that a uniform distribution of porosities has an obvious effect on natural frequencies. Chen et al. [15,16,17] studied the buckling and bending of shear-deformable FGP beams by the Ritz method. Hoa et al. [18] studied nonlinear buckling and post-buckling of cylindrical shells by the three-terms solution and Galerkin’s method. Under various boundary conditions, Chan et al. [19] discussed nonlinear buckling and post-buckling of imperfect FG porous sandwich cylindrical panels by Galerkin’s method.

In order to meet the high efficiency of civil engineering structures and the high precision development of aerospace engineering devices, beams need to be thinner and thinner, and at the same time, it is necessary to improve the material strength of its structure so as to increase the effective space and load. Nanomaterials have good mechanical, thermal, optical, and electrical properties [20], so carbon nanomaterials are often regarded as nanofillers to heighten matrix materials’ properties. These include graphene platelets [21], carbon nanotubes [22], and fullerenes. In 2004, British scientists were the first to peel graphene sheets from graphite with an extremely high tensile, Young modulus, and surface area [23]. On account of the porous structure of metal foam, the stiffness and strength are weakened compared with that of dense metal. By filling the carbon nanomaterials into the matrix materials, the properties of porous materials are able to be efficiently improved. At present, there have been many studies on graphene-reinforced porous materials.

Many papers have been published on graphene-reinforced porous composite beams, plates, and shells. Kitipornchai et al. [24,25] investigated the static and dynamic mechanical behavior of graphene-reinforced FGP beams. Yas et al. [26,27] presented the buckling and vibration of graphene-reinforced FGP beams in thermal environments. Yang et al. [28] studied the buckling and vibration of graphene-reinforced FGP plates. Teng and Wang [29] explained the nonlinear forced vibration of simply supported graphene-reinforced FGP plates. Dong et al. [30] researched the buckling of spinning graphene-reinforced FGP shells. By Galerkin’s method, Zhou et al. [31] revealed the nonlinear buckling of graphene-reinforced FGP cylindrical shells. Under impulsive loading, Yang et al. [32] studied nonlinear forced vibration and the dynamic buckling of graphene-reinforced FGP arches.

Graphene-reinforced composite porous beams and plates are easily affected by the thermal environment, resulting in a decrease in their structural stiffness. Therefore, it is of great significance to study their thermal buckling, free vibration, and dynamic instability for engineering practices. To the best of the authors’ knowledge, no relevant literature has studied the dynamic instability of graphene-enhanced porous materials based on elastic foundations, thermal environments, and axial forces. The present paper mainly investigates the dynamic instability, thermal buckling, and free vibration of functionally graded graphene platelet-reinforced porous beams on an elastic foundation under a thermal environment and axial forces. Three modes of GPL patterns and porosity distributions are considered. Based on the theory of the Timoshenko beam, the governing equation is obtained by the Hamilton principle. On the basis of the differential quadrature method (DQM), the governing equations are changed into Mathieu–Hill equations and the main unstable regions of porous composite beams are studied by the Bolotin method. Moreover, we also use the two-step perturbation method (TSPM) to calculate the thermal buckling and free vibration. The effects of porosity coefficients and GPL’s weight fraction, initial thermal loading, slenderness ratio, geometry and size, boundary conditions, foundation stiffness, and dispersion pattern are discussed.

2. Model Construction

Under axial force $N_{x_0}$ and uniform temperature change $\Delta T=T-T_0$, we consider a FGP multilayer beam that rests on a two-parameter elastic foundation in an initial stress-free state at the reference temperature $T_0$.

As seen in Figure 1, L, b and h, respectively, represent the length, width, and thickness of the beam, and $k_w$ and $k_s$ are the Winkler stiffness and shearing layer stiffness. Among others, the thickness h is divided into n layers, each of which is $▵h=h/n$.

Figure 2 considers three porosity distributions and GPL patterns. Because of disparate dispersion patterns, GPL patterns can be divided as A, B, and C, and the GPLs volume content $V_{GPL}$ is smoothly on the $\mathrm{z}$-axis. According to different porosity distributions, $V_{GPL}$’s peak values can be denoted as $s_{ij}(i,j=1,2,3)$. Assuming three GPL patterns have the same total amount of GPLs will result in $s_{1i}\ne s_{2i}\ne s_{3i}$.

$E_1^{\prime }$ and $E_2^{\prime }$ are denoted as the maximum and minimum elastic moduli of the non-uniform porous beams without GPLs, respectively. In addition, $E^{\prime }$ represents the elastic moduli of the uniform porosity distribution beams.

The relationships of elastic moduli $E\left(z\right)$, mass density $\rho \left(z\right)$, and thermal expansion coefficient $\alpha \left(z\right)$ of FGP beams for three porosity distributions are given by the following formulas [25,33]

$$\left\{\begin{array}{c}E\left(z\right)=E_1[1-e_0\lambda \left(z\right)],\hfill \\ \rho \left(z\right)={\rho }_1[1-e_m\lambda \left(z\right)],\hfill \\ G\left(z\right)=E\left(z\right)/2[1+\nu (z\left)\right],\hfill \\ \alpha \left(z\right)={\alpha }_1\hfill \end{array}\right.$$

(1)

where

$$\lambda \left(z\right)=\left\{\begin{array}{cc}cos\left(\pi \beta \right)\hfill & Porositydistribution1\hfill \\ cos\left(\pi \beta /2+\pi /4\right)\hfill & Porositydistribution2\hfill \\ {\lambda }^{*}\hfill & Porositydistribution3\hfill \end{array}\right.$$

(2)

in which $\beta =z/h$. Further, $E_1$, ${\rho }_1$, and ${\alpha }_1$ are the maximum values of $E\left(z\right)$, $\rho \left(z\right)$, and $\alpha \left(z\right)$, respectively. The porosity coefficient $e_0$ is referred to as

$$e_0=1-\frac{E_2^{\prime }}{E_1^{\prime }},$$

(3)

Based on the Gaussian Random Field (GRF) scheme, the mechanical property of closed-cell cellular solids is denoted by [34]

$$\frac{E\left(z\right)}{E_1}={\left(\frac{\rho \left(z\right)/{\rho }_1+0.121}{1.121}\right)}^{2.3},(0.15<\frac{\rho \left(z\right)}{{\rho }_1}<1)$$

(4)

Using Equation (4), the coefficient of mass density $e_m$ is given by the following formula

$$e_m=\frac{1.121(1-\sqrt[2.3]{1-e_0\lambda \left(z\right)})}{\lambda \left(z\right)},$$

(5)

Similarly, using the closed-cell GRF scheme, Poisson’s ratio $\nu \left(z\right)$ is defined as [35]

$$\nu \left(z\right)=0.221p^{\prime }+{\nu }_1({0.342p\prime }^2-1.21p^{\prime }+1),$$

(6)

where ${\nu }_1$ is Poisson’s ratio of pure non-porous matrix materials and

$$p^{\prime }=1-\frac{\rho \left(z\right)}{{\rho }_1}=1.121(1-\sqrt[2.3]{1-e_0\lambda (z)}),$$

(7)

Due to the total masses being the same for the three porosity distributions, ${\lambda }^{*}$ in Equation (2) can be defined as

$${\lambda }^{*}=\frac{1}{e_0}-\frac{1}{e_0}{\left(\frac{M/{\rho }_{1}h+0.121}{1.121}\right)}^{2.3},$$

(8)

in which M represents all porosity distributions, as shown in the following equation

$$M={\int }_{-h/2}^{h/2}{\rho }_1(1-p^{\prime })\mathrm{d}z,$$

(9)

According to the distribution patterns, the volume fraction of GPLs $V_{GPL}$ is denoted by

$$V_{GPL}=\left\{\begin{array}{cc}{s}_{i1}[1-cos\left(\pi \beta \left(z\right)\right)]\hfill & GPLpatternA\hfill \\ s_{i2}[1-cos(\pi \beta \left(z\right)/2+\pi /4)]\hfill & GPLpatternB\hfill \\ s_{i3}\hfill & GPLpatternC\hfill \end{array}\right.$$

(10)

in which $i=1,2,3$.

The relationship between the weight fraction of GPLs ${\Lambda }_{GPL}$ and the volume fraction of GPLs $V_{GPL}$ is given by

$$\frac{{\Lambda }_{GPL}}{{\Lambda }_{GPL}+\frac{{\rho }_{GPL}}{{\rho }_M}-\frac{{\rho }_{GPL}}{{\rho }_M}{\Lambda }_{GPL}}{\int }_{-h/2}^{h/2}[1-e_m\lambda \left(z\right)]\mathrm{d}z={\int }_{-h/2}^{h/2}{V}_{GPL}[1-e_m\lambda \left(z\right)]\mathrm{d}z$$

(11)

Based on Halpin–Tsai micromechanics model [36,37,38,39], the elastic moduli $E_1$ of the nanocomposites is defined as

$$E_1=\frac{3}{8}\left(\frac{1+{\xi }_L^{GPL}{\eta }_L^{GPL}{V}_{GPL}}{1-{\eta }_L^{GPL}{V}_{GPL}}\right)E_M+\frac{5}{8}\left(\frac{1+{\xi }_W^{GPL}{\eta }_W^{GPL}{V}_{GPL}}{1-{\eta }_W^{GPL}{V}_{GPL}}\right)E_M$$

(12)

where

$$\begin{array}{c}\hfill {\xi }_L^{GPL}=2\frac{a_{GPL}}{t_{GPL}},{\xi }_W^{GPL}=2\frac{b_{GPL}}{t_{GPL}}\\ \hfill {\eta }_L^{GPL}=\frac{(E_{GPL}/E_m)-1}{(E_{GPL}/E_M)+{\xi }_L^{GPL}},\\ \hfill {\eta }_W^{GPL}=\frac{(E_{GPL}/E_m)-1}{(E_{GPL}/E_M)+{\xi }_W^{GPL}}.\end{array}$$

(13)

where $a_{GPL}$, $b_{GPL}$, and $t_{GPL}$ are the average length, width, and thickness of GPLs. $E_M$ and $E_{GPL}$ represent the elastic moduli of the metal and GPLs.

By the following mixture rule, the mass density ${\rho }_1$, Poisson’s ratio ${\upsilon }_1$, and thermal expansion coefficient ${\alpha }_1$ of the metal matrix reinforced by GPLs can be obtained as

$$\begin{array}{cc}\hfill {\rho }_1& ={\rho }_{GPL}{V}_{GPL}+{\rho }_M{V}_M,\hfill \\ \hfill {\upsilon }_1& ={\upsilon }_{GPL}{V}_{GPL}+{\upsilon }_M{V}_M,\hfill \\ \hfill {\alpha }_1& ={\alpha }_{GPL}{V}_{GPL}+{\alpha }_M{V}_M.\hfill \end{array}$$

(14)

in which ${\rho }_M$, ${\upsilon }_M$, ${\alpha }_M$, and $V_M=1-V_{GPL}$ are the mass density, Poisson’s ratio, thermal expansion coefficient, and volume fraction of the metals. Furthermore, ${\rho }_{GPL}$, ${\upsilon }_{GPL}$, ${\alpha }_{GPL}$, and $V_{GPL}$ are the corresponding properties of GPLs.

3. Formulations

3.1. Equations of Governing

Based on the Timoshenko beam theory, the displacement components are expressed as

$$\left\{\begin{array}{c}\overline{U}(x,z,t)=U(x,t)+z\Psi (x,t),\hfill \\ \overline{W}(x,z,t)=W(x,t).\hfill \end{array}\right.$$

(15)

in which U and W are the displacements of the x and z-axes, $\Psi$ is expressed as the normal transverse rotation of the y-axis, and t represents time. On the basis of the linear stress-displacement relationships

$$\left\{\begin{array}{c}{\epsilon }_{xx}=\frac{\partial U}{\partial x}+z\frac{\partial \Psi }{\partial x},\hfill \\ {\gamma }_{xz}=\frac{\partial W}{\partial x}+\Psi .\hfill \end{array}\right.$$

(16)

The linear stress-strain constitutive relationships are as follows

$$\left\{\begin{array}{c}{\sigma }_{xx}=Q_{11}\left(z\right)[{\epsilon }_{xx}-\alpha \left(z\right)\Delta T],\hfill \\ {\sigma }_{xz}=Q_{55}\left(z\right){\gamma }_{xz}.\hfill \end{array}\right.$$

(17)

where elastic elements $Q_{11}\left(z\right)$ and $Q_{55}\left(z\right)$ are denoted by

$$\left\{\begin{array}{c}{Q}_{11}\left(z\right)=\frac{E\left(z\right)}{1-{\upsilon }^2\left(z\right)},\hfill \\ Q_{55}\left(z\right)=G\left(z\right).\hfill \end{array}\right.$$

(18)

Using Hamilton’s principle, the considered problems can be expressed as

$${\int }_0^t\delta (T+V-\Pi )\mathrm{d}t=0$$

(19)

in which

$$\begin{array}{cc}\hfill \delta T& ={\int }_0^L{\int }_{-h/2}^{h/2}\rho \left(z\right)\left[{\left(\frac{\partial \overline{U}}{\partial t}\right)}^2+{\left(\frac{\partial \overline{W}}{\partial t}\right)}^2\right]\mathrm{d}z\mathrm{d}x,\hfill \\ \hfill \delta V& ={\int }_0^L\left[(N_{x0}+N_x^T){\left(\frac{\partial W}{\partial x}\right)}^2\right]\mathrm{d}x,\hfill \\ \hfill \delta \Pi & ={\int }_0^L{\int }_{-h/2}^{h/2}({\sigma }_{xx}{\epsilon }_{xx}+{\sigma }_{xz}{\epsilon }_{xz})\mathrm{d}z\mathrm{d}x+{\int }_0^L\left[k_w{W}^2+k_s{\left(\frac{\partial W}{\partial x}\right)}^2\right]\mathrm{d}x.\hfill \end{array}$$

(20)

$\delta$ stands for the variational symbol, T is the kinetic energy of the beam, V is made up of the axial force $N_{x0}$ and thermally axial force $N_x^T$ due to uniform temperature change $\Delta T$, and II consists of the strain energy of the beam and the elastic potential energy of the foundation.

The governing equations are given by applying Equation (20) to Hamilton’s principle in Equation (19), integrating through the beam thickness

$$\begin{array}{cc}\hfill 0=& {\int }_0^t{\int }_0^L\left[\left(\frac{\partial N_x}{\partial x}-I_1\frac{{\partial }^{2}U}{\partial t^2}-I_2\frac{{\partial }^2\Psi }{\partial t^2}\right)\delta U\right]\mathrm{d}x\mathrm{d}t+{\int }_0^t{\int }_0^L\left[\left(\frac{\partial Q_x}{\partial x}-k_{w}W+k_s\frac{{\partial }^{2}W}{\partial x^2}-(N_{x0}+N_x^T)\frac{{\partial }^{2}W}{\partial x^2}-I_1\frac{{\partial }^{2}W}{\partial t^2}\right)\delta W\right]\mathrm{d}x\mathrm{d}t\hfill \\ \hfill & +{\int }_0^t{\int }_0^L\left[\left(\frac{\partial M_x}{\partial x}-Q_x-I_2\frac{{\partial }^{2}U}{\partial t^2}-I_3\frac{{\partial }^2\Psi }{\partial t^2}\right)\delta \Psi \right]\mathrm{d}x\mathrm{d}t-{\int }_0^t\left[\left(N_x\delta U\right){|}_0^L+(M_x\delta \Psi ){|}_0^L-\left(Q_x-N_{x0}\frac{\partial W}{\partial x}-N_x^T\frac{\partial W}{\partial x}+k_s\frac{\partial W}{\partial x}\right)\delta W{|}_0^L\right]\mathrm{d}t\hfill \end{array}$$

(21)

Letting coefficients $\delta U$, $\delta W$, and $\delta \Psi$ from Equation (21) go to zero separately,

$$\begin{array}{cc}\hfill & \frac{\partial N_x}{\partial x}=I_1\frac{{\partial }^{2}U}{\partial t^2}+I_2\frac{{\partial }^2\Psi }{\partial t^2},\hfill \\ \hfill & \frac{\partial Q_x}{\partial x}-k_{w}W+k_s\frac{{\partial }^{2}W}{\partial x^2}-(N_{x0}+N_x^T)\frac{{\partial }^{2}W}{\partial x^2}=I_1\frac{{\partial }^{2}W}{\partial t^2},\hfill \\ \hfill & \frac{\partial M_x}{\partial x}-Q_x=I_2\frac{{\partial }^{2}U}{\partial t^2}+I_3\frac{{\partial }^2\Psi }{\partial t^2}.\hfill \end{array}$$

(22)

the force and moment are referred to

$$\left\{\begin{array}{c}{N}_x\hfill \\ M_x\hfill \\ Q_x\hfill \end{array}\right\}={\int }_{-h/2}^{h/2}\left\{\begin{array}{c}{\sigma }_{xx}\hfill \\ z{\sigma }_{xx}\hfill \\ {\sigma }_{xz}\hfill \end{array}\right\}\mathrm{d}z$$

(23)

Applying Equations (17) and (18) to Equation (23),

$$\begin{array}{cc}\hfill N_x& =A_{11}\frac{\partial U}{\partial x}+B_{11}\frac{\partial \Psi }{\partial x}-N_x^T,\hfill \\ \hfill M_x& =B_{11}\frac{\partial U}{\partial x}+D_{11}\frac{\partial \Psi }{\partial x}-M_x^T,\hfill \\ \hfill Q_x& =\kappa A_{55}\left(\frac{\partial W}{\partial x}+\Psi \right).\hfill \end{array}$$

(24)

where $\kappa =5/6$ denotes the shear correction factor. The stiffness components, inertia terms, and thermally induced force and moment are defined as

$$\begin{array}{cc}\hfill & (A_{11},B_{11},D_{11})={\int }_{-h/2}^{h/2}{Q}_{11}\left(z\right)(1,z,z^2)\mathrm{d}z=\frac{h}{n}\sum _{i=1}^n\frac{E_1\left(h{\beta }_i\right)[1-e_0\lambda \left(h{\beta }_i\right)]}{1-{\nu }^2\left(h{\beta }_i\right)}(1,h{\beta }_i,{\left(h{\beta }_i\right)}^2),\hfill \\ \hfill & A_{55}={\int }_{-h/2}^{h/2}{Q}_{55}\left(z\right)\mathrm{d}z=\frac{h}{n}\sum _{i=1}^n\kappa \frac{E_1\left(h{\beta }_i\right)[1-e_0\lambda \left(h{\beta }_i\right)]}{2[1+\nu \left(h{\beta }_i\right)]},\hfill \\ \hfill & (I_1,I_2,I_3)={\int }_{-h/2}^{h/2}\rho \left(z\right)(1,z,z^2)\mathrm{d}z=\frac{h}{n}\sum _{i=1}^n{\rho }_1\left(h{\beta }_i\right)[1-e_m\lambda \left(h{\beta }_i\right)](1,h{\beta }_i,{\left(h{\beta }_i\right)}^2),\hfill \\ \hfill & (N_x^T,M_x^T)={\int }_{-h/2}^{h/2}{Q}_{11}\left(z\right)\alpha \Delta T(1,z)\mathrm{d}z=\frac{h}{n}\sum _{i=1}^n\frac{E_1\left(h{\beta }_i\right)[1-e_0\lambda \left(h{\beta }_i\right)]}{1-{\nu }^2\left(h{\beta }_i\right)}\alpha \left(h{\beta }_i\right)\Delta T(1,h{\beta }_i).\hfill \end{array}$$

(25)

in which ${\beta }_i=\frac{1}{2}+\frac{1}{2n}-\frac{i}{n}(i=1,2,3,\dots ,n)$.

The governing equations and related boundary conditions are expressed through Equations (22)–(25)

$$\begin{array}{cc}\hfill & A_{11}\frac{{\partial }^{2}U}{\partial x^2}+B_{11}\frac{{\partial }^2\Psi }{\partial x^2}=I_1\frac{{\partial }^{2}U}{\partial t^2}+I_2\frac{{\partial }^2\Psi }{\partial t^2},\hfill \\ \hfill & \kappa A_{55}\left(\frac{{\partial }^{2}W}{\partial x^2}+\frac{\partial \Psi }{\partial x}\right)-k_{w}W+k_s\frac{{\partial }^{2}W}{\partial x^2}-(N_{x0}+N_x^T)\frac{{\partial }^{2}W}{\partial x^2}=I_1\frac{{\partial }^{2}W}{\partial t^2},\hfill \\ \hfill & B_{11}\frac{{\partial }^{2}U}{\partial x^2}+D_{11}\frac{{\partial }^2\Psi }{\partial x^2}-\kappa A_{55}\left(\frac{\partial W}{\partial x}+\Psi \right)=I_2\frac{{\partial }^{2}U}{\partial t^2}+I_3\frac{{\partial }^2\Psi }{\partial t^2}.\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & Clamped\left(C\right):U=W=\Psi =0,\hfill \\ \hfill & Hinged\left(H\right):U=W=M_x=0.\hfill \end{array}$$

(27)

By introducing dimensionless quantities

$$\begin{array}{cc}\hfill & \xi =\frac{x}{L},(u,w)=\frac{(U,W)}{h},\eta =\frac{L}{h},\psi =\Psi ,(P,P^T,M^T)=(N_{x0},N_x^T,M_x^T/h)/A_{110},\hfill \\ \hfill & ({\overline{I}}_1,{\overline{I}}_2,{\overline{I}}_3)=\left(\frac{I_1}{I_{10}},\frac{I_2}{I_{10}h},\frac{I_3}{I_{10}{h}^2}\right),(a_{11},a_{55},b_{11},d_{11})=\left(\frac{A_{11}}{A_{110}},\frac{A_{55}}{A_{110}},\frac{B_{11}}{A_{110}h},\frac{D_{11}}{A_{110}{h}^2}\right),\hfill \\ \hfill & K_w=\frac{k_w{L}^2}{A_{110}},K_s=\frac{k_s}{A_{110}},\omega =\Omega L\sqrt{\frac{I_{10}}{A_{110}}},\tau =\frac{t}{L}\sqrt{\frac{A_{110}}{I_{10}}}.\hfill \end{array}$$

(28)

in which $A_{110}$ and $I_{10}$ are the value of $A_{11}$ and $I_1$ of the pure metal beams without any pores and nanofillers.

Then Equations (26) and (27) are rewritten into the following dimensionless form

$$\begin{array}{cc}\hfill & a_{11}\frac{{\partial }^{2}u}{\partial {\xi }^2}+b_{11}\frac{{\partial }^2\psi }{\partial {\xi }^2}={\overline{I}}_1\frac{{\partial }^{2}u}{\partial {\tau }^2}+{\overline{I}}_2\frac{{\partial }^2\psi }{\partial {\tau }^2},\hfill \\ \hfill & \kappa a_{55}\left(\frac{{\partial }^{2}w}{\partial {\xi }^2}+\eta \frac{\partial \psi }{\partial \xi }\right)-K_{w}w+K_s\frac{{\partial }^{2}w}{\partial {\xi }^2}-(P+P^T)\frac{{\partial }^{2}w}{\partial {\xi }^2}={\overline{I}}_1\frac{{\partial }^2\overline{w}}{\partial {\tau }^2},\hfill \\ \hfill & b_{11}\frac{{\partial }^{2}u}{\partial {\xi }^2}+d_{11}\frac{{\partial }^2\psi }{\partial {\xi }^2}-\kappa \eta a_{55}\left(\frac{\partial w}{\partial \xi }+\eta \psi \right)={\overline{I}}_2\frac{{\partial }^{2}u}{\partial {\tau }^2}+{\overline{I}}_3\frac{{\partial }^2\psi }{\partial {\tau }^2}.\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & Clamped\left(C\right):u=w=\psi =0,\hfill \\ \hfill & Hinged\left(H\right):u=w=b_{11}\frac{\partial u}{\partial \xi }+d_{11}\frac{\partial \psi }{\partial \xi }-M^T=0.\hfill \end{array}$$

(30)

3.2. Solution Method

Based on the DQM rule, the displacement components $u,w$, and $\psi$, and the rth-order partial derivative is estimated in the following form [40,41]

$$\begin{array}{cc}\hfill & (u,w,\psi )=\sum _{m=1}^N{l}_m\left(\xi \right)(u_m,w_m,{\psi }_m),\hfill \\ \hfill & \frac{{\partial }^j}{\partial {\xi }^j}(u,w,\psi )=\sum _{m=1}^N{c}_{im}^{\left(j\right)}(u_m,w_m,{\psi }_m),i=1,\dots ,N,j=1,\dots ,N-1.\hfill \end{array}$$

(31)

in which $(u_m,w_m,{\psi }_m)$ are the values of $(u,w,\psi )$, and $l_m\left(\xi \right)$ is the Lagrange interpolating polynomial; $c_{im}^{\left(j\right)}$ is the weighting coefficient of the jth-order derivative [42]. N is the total number of grid points along the $\xi$-axis. The distribution of the $\xi$-axis is defined by a cosine pattern

$${\xi }_i=\frac{1}{2}\left[1-cos\left(\frac{i-1}{N-1}\pi \right)\right],i=1,2,\dots ,N.$$

(32)

By taking Equations (31) and (32) into Equations (29) and (30), the governing equations and boundary conditions are written as

$$\begin{array}{cc}\hfill & a_{11}\sum _{m=1}^N{C}_{im}^{\left(2\right)}{u}_m+b_{11}\sum _{m=1}^N{C}_{im}^{\left(2\right)}{\psi }_m={\overline{I}}_1\ddot{u_i}+{\overline{I}}_2\ddot{{\psi }_i},\hfill \\ \hfill & \kappa a_{55}\left(\sum _{m=1}^N{C}_{im}^{\left(2\right)}{w}_m+\eta \sum _{m=1}^N{C}_{im}^{\left(1\right)}{\psi }_m\right)-K_w{w}_i+K_s\sum _{m=1}^N{C}_{im}^{\left(2\right)}{w}_m-(P+P^T)\sum _{m=1}^N{C}_{im}^{\left(2\right)}{w}_m={\overline{I}}_1\ddot{w_i},\hfill \\ \hfill & b_{11}\sum _{m=1}^N{C}_{im}^{\left(2\right)}{u}_m+d_{11}\sum _{m=1}^N{C}_{im}^{\left(2\right)}{\psi }_m-\kappa \eta a_{55}\left(\sum _{m=1}^N{C}_{im}^{\left(1\right)}{w}_m+\eta {\psi }_i\right)={\overline{I}}_2\ddot{u_i}+{\overline{I}}_3\ddot{{\psi }_i}.\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill & u_1=w_1={\psi }_1=0,\hfill \\ \hfill & u_N=w_N={\psi }_N=0.\hfill \end{array}$$

(34)

clamped at both ends of $\xi =0,1$.

$$\begin{array}{cc}\hfill & u_1=w_1=b_{11}\sum _{m=1}^N{C}_{1m}^{\left(1\right)}{u}_m+d_{11}\sum _{m=1}^N{C}_{1m}^{\left(1\right)}{\psi }_m-\eta M^T{|}_{\xi ={\xi }_1}=0,\hfill \\ \hfill & u_N=w_N=b_{11}\sum _{m=1}^N{C}_{Nm}^{\left(1\right)}{u}_m+d_{11}\sum _{m=1}^N{C}_{Nm}^{\left(1\right)}{\psi }_m-\eta M^T{|}_{\xi ={\xi }_N}=0.\hfill \end{array}$$

(35)

hinged at both ends of $\xi =0,1$, where $C_{im}^{\left(1\right)}$ and $C_{im}^{\left(2\right)}$ represent the first- and second-order weighting coefficients.

By combining the discretized governing equation, Equation (33), and boundary conditions, Equations (34) and (35), a series of dimensionless algebraic formulas has been obtained as

$$\mathbf{M}\ddot{\mathbf{d}}+\left({\mathbf{K}}_{\mathrm{L}}-\Delta T{\mathbf{K}}_{\mathrm{T}}-P{\mathbf{K}}_{\mathrm{p}}\right)\mathbf{d}=0$$

(36)

where $\mathbf{d}={\left\{\left\{u_1,u_2,\dots ,u_N\right\},\left\{w_1,w_2,\dots ,w_N\right\},\left\{{\psi }_1,{\psi }_2,\dots ,{\psi }_N\right\}\right\}}^T$ is the unknown coefficient vector, ${\mathbf{K}}_{\mathrm{L}}$ and $\mathbf{M}$ represent the stiffness matrix and the mass matrix, and ${\mathbf{K}}_{\mathrm{T}}$ and ${\mathbf{K}}_{\mathrm{p}}$ represent the geometric stiffness matrix.

For the beam under a time-varying axial excitation, the dimensionless axial force P is defined as

$$P=P_s+P_{d}cos\theta \tau$$

(37)

in which $P_s$ and $P_d$ represent the static and dynamic force components. By putting Equation (37) into Equation (36), we have

$$\mathbf{M}\ddot{\mathbf{d}}+\left[{\mathbf{K}}_{\mathrm{L}}-\Delta T{\mathbf{K}}_{\mathrm{T}}-(P_s+P_{d}cos\theta \tau ){\mathbf{K}}_{\mathrm{p}}\right]\mathbf{d}=0$$

(38)

Under axial force and initial thermal loading, Equation (38) is a Mathieu–Hill-type equation, which is used to solve the problems of dynamic instability of the FGP beams. The boundary of the unstable region is obtained by Bolotin’s method [43]. According to the previous study, the solution with period 2 $T_{\theta }$ ($T_{\theta }=2\pi /\theta$) has a larger principal unstable region than the solution with period $T_{\theta }$, which is closer to the practical engineering significance. The solution to Equation (38) with period 2 $T_{\theta }$ uses the trigonometric form

$$\mathbf{d}=\sum _{k=1,3,..}^{\infty }{\mathbf{a}}_{k}sin\left(\frac{k\theta \tau }{2}\right)+{\mathbf{b}}_{k}cos\left(\frac{k\theta \tau }{2}\right)$$

(39)

where $a_k$ and $b_k$ are arbitrary constant vectors. Bolotin verified that the first-order approximation of $k=1$ accurately describes the boundary of the unstable region [43], so a homogeneous linear system of equations represented by $a_1$ and $b_1$ can be obtained by Bolotin’s method

$$\begin{array}{cc}\hfill & \left[{\mathbf{K}}_{\mathrm{L}}-\Delta T{\mathbf{K}}_{\mathrm{T}}-(P_s-\frac{P_d}{2}){\mathbf{K}}_{\mathrm{p}}-\frac{{\theta }^2}{4}\mathbf{M}\right]{\mathbf{a}}_1=0,\hfill \\ \hfill & \left[{\mathbf{K}}_{\mathrm{L}}-\Delta T{\mathbf{K}}_{\mathrm{T}}-(P_s+\frac{P_d}{2}){\mathbf{K}}_{\mathrm{p}}-\frac{{\theta }^2}{4}\mathbf{M}\right]{\mathbf{b}}_1=0.\hfill \end{array}$$

(40)

For a given axial force, Equation (40) gives two critical excitation frequencies. The two curves in the figure of $\theta$ and $P_d$ are used to describe the principle unstable regions. When $P_d=0$, it represents the origin of the principle unstable region, and $\theta$ represents the doubled fundamental frequency of the beam.

As for the thermal buckling problem, we form the equation by neglecting the inertia terms and making $P_s=P_d=0$ from Equation (38). Thus the critical buckling temperature rise can be obtained by solving the minimum positive eigenvalue of Equation (38)

$$\left({\mathbf{K}}_{\mathrm{L}}-\Delta T{\mathbf{K}}_{\mathrm{T}}\right)\mathbf{d}=0$$

(41)

Like thermal buckling, by setting $P_d=0$ and letting $\mathbf{d}={\mathbf{d}}^{*}{e}^{i\omega t}$, the free frequency of the beam makes from the following formula

$$\left({\mathbf{K}}_{\mathrm{L}}-\Delta T{\mathbf{K}}_{\mathrm{T}}-P_s{\mathbf{K}}_{\mathrm{p}}-{\omega }^2\mathbf{M}\right){\mathbf{d}}^{*}=0$$

(42)

4. Discussion

The effects of various factors on thermal buckling, thermo-mechanical vibration, and dynamic instability of FGP beams are discussed. Copper is often chosen as a matrix material, the material parameters of which are $E_M$ = 130 GPa, ${\rho }_M$ = 8960 kg/m${}^3$, $v_M$ = 0.34, and ${\alpha }_M$ = 17 $\times {10}^{-6}$${\mathrm{K}}^{-1}$. The material and size parameters of GPLs as reinforced materials are $w_{GPL}$ = 1.5 $\mathsf{\mu }$m, $l_{GPL}$ = 2.5 $\mathsf{\mu }$m, $t_{GPL}$ = 1.5 nm, $E_{GPL}$ = 1.01 TPa, ${\rho }_{GPL}$ = 1062.5 kg/m${}^3$, $v_{GPL}$ = 0.186, and ${\alpha }_{GPL}$ = −3.75 $\times {10}^{-6}$${\mathrm{K}}^{-1}$ [22,44,45,46].

4.1. Validation and Convergence Study

First, validation analysis is conducted. We adopt the degenerate forms of $K_w=0$, $K_s=0$, and $\Delta T=0$ to compare and validate with references [25,47].

Table 1. The dimensionless fundamental frequency of C-C beams (porosity 1, GPL A, $L/h=20$).

${\mathsf{\Lambda }}_{\mathit{GPL}}$	n	${\mathit{e}}_0=0$		${\mathit{e}}_0=0.2$		${\mathit{e}}_0=0.4$		${\mathit{e}}_0=0.6$
		Ref. [25]	Present	Ref. [25]	Present	Ref. [25]	Present	Ref. [25]	Present
1 wt.%	2	0.3376	0.3376	0.3217	0.3217	0.3042	0.3042	0.2845	0.2845
	6	0.4390	0.4390	0.4336	0.4336	0.4289	0.4289	0.4259	0.4259
	10	0.4464	0.4464	0.4421	0.4421	0.4388	0.4388	0.4372	0.4372
	14	0.4484	0.4484	0.4444	0.4444	0.4415	0.4415	0.4403	0.4403
	18	0.4492	0.4492	0.4454	0.4454	0.4426	0.4426	0.4416	0.4416
	10,000	0.4505	0.4505	0.4468	0.4468	0.4442	0.4442	0.4436	0.4436
0 wt.%	14	0.3159	0.3159	0.3134	0.3134	0.3121	0.3121	0.3128	0.3128
	10,000	0.3167	0.3167	0.3144	0.3144	0.3132	0.3132	0.3142	0.3142

and

Table 2. The dimensionless critical buckling load of C-C beams (porosity 1, GPL A, $L/h=20$).

${\mathsf{\Lambda }}_{\mathit{GPL}}$	n	${\mathit{e}}_0=0$		${\mathit{e}}_0=0.2$		${\mathit{e}}_0=0.4$		${\mathit{e}}_0=0.6$
		Ref. [25]	Present	Ref. [25]	Present	Ref. [25]	Present	Ref. [25]	Present
1 wt.%	2	0.008550	0.008550	0.007218	0.007219	0.005917	0.005918	0.004647	0.004647
	6	0.014323	0.014324	0.013057	0.013058	0.011784	0.011784	0.010486	0.010486
	10	0.014798	0.014798	0.013572	0.013573	0.012333	0.012334	0.011063	0.011063
	14	0.014929	0.014929	0.013714	0.013715	0.012486	0.012486	0.011224	0.011224
	18	0.014982	0.014983	0.013773	0.013774	0.012549	0.012549	0.011290	0.011290
	10,000	0.015065	0.015066	0.013863	0.013864	0.012645	0.012646	0.011392	0.011392
0 wt.%	14	0.007946	0.007947	0.007316	0.007316	0.006693	0.006693	0.006076	0.006076
	10,000	0.007986	0.007986	0.007362	0.007362	0.006745	0.006746	0.006135	0.006135

compare the fundamental frequency and critical buckling load with the calculation results in reference [25]. Figure 3 verifies the dynamic instability of Wu et al. [47]. In general, our results are consistent with the existing results.

Figure 4a,b separately show the convergence results of the critical buckling temperature rise and the dimensionless natural frequency under various conditions. When $N=9$, their values gradually approach a certain amount.

Table 3. The critical buckling temperature rise of H-H beams (porosity 1, GPL A, $L/h=20$, $K_w=K_s=0$).

${\mathsf{\Lambda }}_{\mathit{GPL}}$	n	${\mathit{e}}_0=0$			${\mathit{e}}_0=0.2$			${\mathit{e}}_0=0.4$			${\mathit{e}}_0=0.6$
		DQM	TSPM	Error	DQM	TSPM	Error	DQM	TSPM	Error	DQM	TSPM	Error
1 wt.%	2	97.7638	97.7646	0.0008%	97.5424	97.5433	0.0009%	97.2608	97.2616	0.0008%	96.8839	96.8847	0.0008%
	6	162.4627	162.4641	0.0009%	170.1854	170.1868	0.0008%	179.9406	179.9421	0.0008%	192.7597	192.7613	0.0008%
	10	167.7434	167.7448	0.0008%	176.4998	176.5012	0.0008%	187.5116	187.5131	0.0008%	201.8757	201.8774	0.0008%
	14	169.1999	169.2013	0.0008%	178.2459	178.2474	0.0008%	189.6084	189.6100	0.0008%	204.4014	204.4031	0.0008%
	18	169.7995	169.8009	0.0008%	178.9652	178.9667	0.0008%	190.4726	190.4742	0.0008%	205.4423	205.4440	0.0008%
	1000	170.7175	170.7189	0.0008%	180.0671	180.0686	0.0008%	191.7968	191.7984	0.0008%	207.0375	207.0392	0.0008%
0 wt.%	14	119.4444	119.4454	0.0008%	127.2565	127.2576	0.0009%	137.5587	137.5598	0.0008%	151.8699	151.8712	0.0008%
	1000	120.0523	120.0533	0.0008%	128.0345	128.0356	0.0009%	138.5566	138.5578	0.0009%	153.1614	153.1627	0.0008%

and

Table 4. The dimensionless fundamental frequency of H-H beams (porosity 1, GPL A, $L/h=20$, $K_w=K_s=0$, $\Delta T=50\mathrm{K}$).

${\mathsf{\Lambda }}_{\mathit{GPL}}$	n	${\mathit{e}}_0=0$			${\mathit{e}}_0=0.2$			${\mathit{e}}_0=0.4$			${\mathit{e}}_0=0.6$
		DQM	TSPM	Error	DQM	TSPM	Error	DQM	TSPM	Error	DQM	TSPM	Error
1 wt.%	2	0.1053	0.1051	0.1899%	0.1002	0.1000	0.1996%	0.0946	0.0944	0.2114%	0.0882	0.0881	0.1134%
	6	0.1639	0.1636	0.1830%	0.1636	0.1633	0.1834%	0.1638	0.1635	0.1832%	0.1649	0.1645	0.2426%
	10	0.1679	0.1676	0.1787%	0.1682	0.1678	0.2378%	0.1690	0.1686	0.2367%	0.1707	0.1704	0.1757%
	14	0.1690	0.1687	0.1775%	0.1694	0.1691	0.1771%	0.1704	0.1700	0.2347%	0.1723	0.1720	0.1741%
	18	0.1695	0.1691	0.2360%	0.1699	0.1696	0.1766%	0.1710	0.1706	0.2339%	0.1730	0.1726	0.2312%
	1000	0.1702	0.1698	0.2350%	0.1707	0.1703	0.2343%	0.1718	0.1715	0.1746%	0.1740	0.1736	0.2299%
0 wt.%	14	0.1078	0.1076	0.1855%	0.1094	0.1092	0.1828%	0.1117	0.1114	0.2686%	0.1151	0.1148	0.2606%
	1000	0.1083	0.1081	0.1847%	0.1100	0.1097	0.2727%	0.1123	0.1121	0.1781%	0.1158	0.1156	0.1727%

present the effect of the porosity coefficient $e_0$ and the total number of layers n on the critical buckling temperature rise and dimensionless natural frequency by DQM and the two-step perturbation method (TSPM). It turns out that the error between them is within $0.1\%$, and the accuracy and efficiency of the calculation results are verified again. Suppose n = 1000 is a continuous beam; it is found that a relative difference between n = 14 and n = 1000 is less than 1.5%. Considering the manufacturing process and manufacturing costs, n = 14 and $N=9$ are used in the following calculation. In addition, when n = 14 and $N=9$, the results show the natural frequency and critical buckling temperature rise are both increasing as $e_0$ increases.

4.2. Thermal Buckling

Figure 5 examines the effect of the foundation’s stiffness in critical buckling temperature rise. Where $(K_w,K_s)=(0,0)$ stands for no foundation, $(K_w,K_s)=(0.1,0)$ stands for Winkler foundation, and $(K_w,K_s)=(0.1,0.02)$ stands for Pasternak foundation. As observed, the critical buckling temperature increases as the foundation stiffness increases. The shearing layer stiffness $K_s$ contributes to more enhancement than the Winkler foundation stiffness $K_w$.

Figure 6a,b show the critical buckling temperature rise and its percentage increment at GPL weight fraction ${\Lambda }_{GPL}$. The results identify that symmetric GPL A with porosity 1 provides the best reinforcement, which takes the largest critical buckling temperature rise of the nine models. In addition, GPL C plays no role in the critical buckling temperature rise for the three porosity distributions.

Table 5. Effect of boundary conditions and slenderness ratio $L/h$ on critical buckling temperature rise. (Porosity 1, GPL A, $e_0=0.5$, ${\Lambda }_{GPL}=1.0\mathrm{wt.\%}$, $K_w=K_s=0$).

BC	$\mathit{L}/\mathit{h}=10$	$\mathit{L}/\mathit{h}=15$	$\mathit{L}/\mathit{h}=20$	$\mathit{L}/\mathit{h}=25$	$\mathit{L}/\mathit{h}=30$
C-C	2730.8457	1316.1778	762.8928	495.2309	346.6018
C-H	1487.9949	693.6980	397.0135	256.1591	178.6797
H-H	762.8528	346.5832	196.4819	126.2065	87.8172

illustrates the effect of boundary conditions and slenderness ratio $L/h$ on critical buckling temperature rise. As expected, the C-C beam with a smaller slenderness ratio has a maximum critical buckling temperature rise. With the increment of $L/h$, the result shows a downward trend.

Figure 7 depicts the effect of geometry and size $a_{GPL}/b_{GPL}$ and $b_{GPL}/t_{GPL}$ on the critical buckling temperature rise. Larger $a_{GPL}/b_{GPL}$ and $b_{GPL}/t_{GPL}$ efficiently enhance the critical buckling temperature rise. Moreover, when $b_{GPL}/t_{GPL}$ reaches ${10}^3$, the critical buckling temperature rise reaches a certain level, it will no longer increase any more, and as $a_{GPL}/b_{GPL}$ increases, the change in critical buckling temperature rise becomes less and less obvious.

4.3. Thermo-Mechanical Vibration

Figure 8 presents the effect of GPL weight fractions ${\Lambda }_{GPL}$ and normalized static axial force $P_s/P_{cr}$ on the dimensionless fundamental frequency. The positive and negative values of $P_s/P_{cr}$ indicate the compressive force and tensile forces. $P_{cr}$ represents the critical buckling load at $\Delta T=0$ K. As observed, increasing ${\Lambda }_{GPL}$ leads to better mechanical behavior.

Table 6. Effect of normalized static axial force $P_s/P_{cr}$ on the dimensionless fundamental frequency under porosity distributions and GPL patterns ($e_0=0.5$, ${\Lambda }_{GPL}$ = 1.0 wt.%, C-C, $L/h=20$, $\Delta T=0K$, $K_w=K_s=0$).

Multilayer Beam	1A	2A	3A	1B	2B	3B	1C	2C	3C
$P_s/P_{cr}=0.00$	0.4264	0.3863	0.3863	0.3748	0.3318	0.3389	0.3829	0.3452	0.3441
$P_s/P_{cr}=0.25$	0.3668	0.3314	0.3314	0.3213	0.2830	0.2896	0.3283	0.2949	0.2941
$P_s/P_{cr}=0.50$	0.2939	0.2641	0.2641	0.2556	0.2227	0.2287	0.2613	0.2329	0.2323

examines the effect of a normalized static axial force $P_s/P_{cr}$ on the dimensionless fundamental frequency under porosity distributions and GPL patterns. Due to compression forces compressing the beam, the free vibration frequency decreases when the compression force increases. Like thermal buckling, GPL A with porosity 1 has the highest free-vibration frequency.

The effects of normalized static axial force $P_s/P_{cr}$ on the dimensionless fundamental frequency under initial thermal loading $\Delta T$ are shown in Figure 9. With the increment of initial thermal loading, the overall trend is downward. The changes in the dimensionless fundamental frequency were more apparent at larger values of $P_s/P_{cr}$.

Figure 10 reveals the effects of normalized static axial force $P_s/P_{cr}$ on the dimensionless fundamental frequency under various foundation stiffness. The dimensionless fundamental frequency increases as the foundation stiffness increases. Compared with compressive force, tensile force enhances the dimensionless fundamental frequency, which strengthens the stiffness of the beam.

4.4. Dynamic Instability

Figure 11 illustrates the effect of porosity distributions and GPL patterns on dynamic instability. Like thermal buckling and thermo-mechanical vibration, GPL A with porosity 1 has the largest origin and the narrowest unstable region among the nine models. Then, it is the best-enhanced model because the beam has a smaller pore distribution and more GPL in the top and bottom layers, which is where the normal bending stress is highest.

Figure 12 depicts the effect of GPL weight fractions ${\Lambda }_{GPL}$ on dynamic instability. It is found that the origin becomes larger and the unstable region becomes narrower with the increment of ${\Lambda }_{GPL}$, which indicates that the addition of GPL nanofillers effectively raises the stiffness of the beam.

Figure 13 investigates the effect of the porosity coefficient $e_0$ on dynamic instability. As observed, the increment in $e_0$, which means that the beam has larger pores and a denser distribution of pores, causes a reduction in beam stiffness, the origin becomes lower, and the unstable region becomes wider.

The effect of static axial compressive force $P_s/P_{cr}$ and initial thermal loading $\Delta T$ on dynamic instability are shown in Figure 14 and Figure 15. Due to the change in $P_s/P_{cr}$ and $\Delta T$, the beam produces compression force, thus reducing the beam stiffness. As $P_s/P_{cr}$ and $\Delta T$ decrease, the origin increases and the unstable region narrows. $P_s/P_{cr}$ is able to achieve better structural stiffness, even more significantly than $\Delta T$ at the dynamic instability.

The effect of the slenderness ratio $L/h$, foundation stiffness, and boundary conditions on dynamic instability is demonstrated in Figure 16, Figure 17 and Figure 18. The results depict that the C-C beam with a smaller slenderness ratio on a Pasternak elastic foundation has a bigger origin and narrower unstable region, and the slenderness ratio and foundation stiffness have more obvious effects than the boundary conditions.

Figure 19 examines the effect of the geometry and size $a_{GPL}/b_{GPL}$ and $b_{GPL}/t_{GPL}$ on dynamic instability. For a given $b_{GPL}/t_{GPL}={10}^2$, as $a_{GPL}/b_{GPL}$ increases, the change in the unstable region is insignificant, and the effect is tiny. For a given $a_{GPL}/b_{GPL}=4$, the change in $b_{GPL}/t_{GPL}$ from 10 to ${10}^2$ has an obvious effect on the dynamic instability, but the effect of $b_{GPL}/t_{GPL}$ over ${10}^2$ is negligible. The results show that when the GPL contains less than a single graphene layer, it is better to reinforce the stiffness of the beam and enhance the mechanical behavior.

5. Conclusions

The effect of GPL nanofillers on FGP composite beams under thermal environments, thermal buckling, thermo-mechanical vibration, and dynamic instability are investigated. Among them, weight fraction, normalized static axial force, porosity coefficient, dispersion pattern, boundary conditions, initial thermal loading, geometry and size, foundation stiffness, and slenderness ratio are studied. The following conclusions are obtained:

•

Porosity 1 reinforced by GPL A of the beam has the biggest value of critical buckling, temperature rise, dimensionless fundamental frequency, and the origin of dynamic instability. The non-uniform, symmetric porosity distribution and GPL pattern have the strongest enhancement.

•

The porosity coefficient has an important influence on thermal buckling, thermo-mechanical vibration, and dynamic instability. When the porosity coefficient grows, the origin of the dynamic instability shows a decreasing trend, but the dimensionless fundamental frequency and critical buckling temperature rise both increase.

•

The addition of GPL nanofillers can enhance the beam stiffness significantly, and the mechanical performance is enhanced with ${\Lambda }_{GPL}$ increases.

•

The values of thermo-mechanical vibration and dynamic instability decrease with normalized static axial force and initial thermal loading increase.

•

Winkler and Pasternak foundations both strengthen the stiffness of the beam. It is noted that shearing layer stiffness has a better enhancement effect than Winkler foundation stiffness.