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

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.


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,

Model Construction
Under axial force N x 0 and uniform temperature change ∆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.  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 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 = s 2i = s 3i . E 1 and E 2 are denoted as the maximum and minimum elastic moduli of the nonuniform porous beams without GPLs, respectively. In addition, E represents the elastic moduli of the uniform porosity distribution beams. The relationships of elastic moduli E(z), mass density ρ(z), and thermal expansion coefficient α(z) of FGP beams for three porosity distributions are given by the following formulas [25,33]  where cos(πβ/2 + π/4) Porosity distribution 2 λ * Porosity distribution 3 (2) in which β = z/h. Further, E 1 , ρ 1 , and α 1 are the maximum values of E(z), ρ(z), and α(z), respectively. The porosity coefficient e 0 is referred to as Based on the Gaussian Random Field (GRF) scheme, the mechanical property of closed-cell cellular solids is denoted by [34] Using Equation (4), the coefficient of mass density e m is given by the following formula Similarly, using the closed-cell GRF scheme, Poisson's ratio ν(z) is defined as [35] ν where ν 1 is Poisson's ratio of pure non-porous matrix materials and Due to the total masses being the same for the three porosity distributions, λ * in Equation (2) can be defined as in which M represents all porosity distributions, as shown in the following equation According to the distribution patterns, the volume fraction of GPLs V GPL is denoted by in which i = 1, 2, 3. The relationship between the weight fraction of GPLs Λ GPL and the volume fraction of GPLs V GPL is given by Based on Halpin-Tsai micromechanics model [36][37][38][39], the elastic moduli E 1 of the nanocomposites is defined as where 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 ρ 1 , Poisson's ratio υ 1 , and thermal expansion coefficient α 1 of the metal matrix reinforced by GPLs can be obtained as in which ρ M , υ M , α M , and V M = 1 − V GPL are the mass density, Poisson's ratio, thermal expansion coefficient, and volume fraction of the metals. Furthermore, ρ GPL , υ GPL , α GPL , and V GPL are the corresponding properties of GPLs.

Equations of Governing
Based on the Timoshenko beam theory, the displacement components are expressed as in which U and W are the displacements of the x and z-axes, Ψ is expressed as the normal transverse rotation of the y-axis, and t represents time. On the basis of the linear stressdisplacement relationships The linear stress-strain constitutive relationships are as follows where elastic elements Q 11 (z) and Q 55 (z) are denoted by Using Hamilton's principle, the considered problems can be expressed as in which δ 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 T x due to uniform temperature change ∆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 Letting coefficients δU, δW, and δΨ from Equation (21) go to zero separately, the force and moment are referred to Applying Equations (17) and (18) to Equation (23), where κ = 5/6 denotes the shear correction factor. The stiffness components, inertia terms, and thermally induced force and moment are defined as in which β i = 1 2 + 1 2n − i n (i = 1, 2, 3, . . . , n). The governing equations and related boundary conditions are expressed through Equations (22)-(25) By introducing dimensionless quantities 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

Solution Method
Based on the DQM rule, the displacement components u, w, and ψ, and the rth-order partial derivative is estimated in the following form [40,41] in which (u m , w m , ψ m ) are the values of (u, w, ψ), and l m (ξ) is the Lagrange interpolating polynomial; c (j) im is the weighting coefficient of the jth-order derivative [42]. N is the total number of grid points along the ξ-axis. The distribution of the ξ-axis is defined by a cosine pattern By taking Equations (31) and (32) into Equations (29) and (30), the governing equations and boundary conditions are written as clamped at both ends of ξ = 0, 1.
hinged at both ends of ξ = 0, 1, where C (1) im and C (2) im 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 . . , ψ N }} T is the unknown coefficient vector, K L and M represent the stiffness matrix and the mass matrix, and K T and K 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 θτ (37) in which P s and P d represent the static and dynamic force components. By putting Equation (37) into Equation (36), we have 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 θ (T θ = 2π/θ) has a larger principal unstable region than the solution with period T θ , which is closer to the practical engineering significance. The solution to Equation (38) with period 2 T θ uses the trigonometric form 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 For a given axial force, Equation (40) gives two critical excitation frequencies. The two curves in the figure of θ 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 θ 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) Like thermal buckling, by setting P d = 0 and letting d = d * e iωt , the free frequency of the beam makes from the following formula

Discussion
The

Validation and Convergence Study
First, validation analysis is conducted. We adopt the degenerate forms of K w = 0, K s = 0, and ∆T = 0 to compare and validate with references [25,47]. Tables 1 and 2 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.     Tables 3 and 4 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.    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 Λ 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 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. Table 5. Effect of boundary conditions and slenderness ratio L/h on critical buckling temperature rise. (Porosity 1, GPL A, e 0 = 0.5, Λ GPL = 1.0 wt.%, K w = K s = 0).  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.  Figure 8 presents the effect of GPL weight fractions Λ 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 ∆T = 0 K. As observed, increasing Λ GPL leads to better mechanical behavior. Table 6 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. Figure 8. Effect of GPL weight fractions Λ GPL and normalized static axial force P s /P cr on the dimensionless fundamental frequency. 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 , Λ GPL = 1.0 wt.% , C-C, L/h = 20, ∆T = 0 K, K w = K s = 0). The effects of normalized static axial force P s /P cr on the dimensionless fundamental frequency under initial thermal loading ∆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 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 11. Effect of porosity distributions and GPL patterns on dynamic instability. Figure 12 depicts the effect of GPL weight fractions Λ GPL on dynamic instability. It is found that the origin becomes larger and the unstable region becomes narrower with the increment of Λ 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.

Dynamic Instability
The effect of static axial compressive force P s /P cr and initial thermal loading ∆T on dynamic instability are shown in Figures 14 and 15. Due to the change in P s /P cr and ∆T, the beam produces compression force, thus reducing the beam stiffness. As P s /P cr and ∆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 ∆T at the dynamic instability.   The effect of the slenderness ratio L/h, foundation stiffness, and boundary conditions on dynamic instability is demonstrated in Figures 16-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.

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, thermomechanical 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 Λ 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.