Dynamic Stability of Temperature-Dependent Graphene Sheet Embedded in an Elastomeric Medium

Abstract

This work applies the first-order shear deformation theory (FSDT) to study the dynamic stability of orthotropic temperature-dependent single-layered graphene sheet (SLGS) embedded in a temperature-dependent elastomeric medium and subjected to a biaxial oscillating loading in a thermal environment. Possible thermal effects are considered in the size-dependent governing equations of the problem. These last ones are derived by means of the Hamilton’s variational principle combined with the Eringen’s differential constitutive model. Navier’s solution as well as Bolotin’s approach are applied to obtain the dynamic instability region (DIR) of the graphene sheet. Thus, a parametric study is carried out to explore the sensitivity of the DIR of the graphene sheet to the temperature variation, the static load factor, the aspect ratio, the foundation type, and the nonlocal parameter (NP). Results indicate that the dimensionless pulsation frequency reduces for increasing values of temperature and NP, whereas the size effect becomes even more pronounced for increasing temperatures. In addition, the adoption of temperature-dependent mechanical properties, rather than independent ones, yields a global shift of the DIR to smaller pulsating frequencies. This proves the relevance of the temperature-dependent mechanical properties to obtain reliable results, in a physical sense.

1. Introduction

Among nanostructures, graphene sheets have increased the interest of the scientific community, owing to their astonishing thermal, chemical, electrical, and mechanical properties. Due to the distinguishing features, theses nanostructures have been adopted as important components in various high technology industries such as nanoactuators, nanoresonators, nanosensors, and nanocomposites.

In order to investigate the mechanical characteristics of micro/nano-scale structures accurately, the small-scale effects should be considered. It has been recognized that the main properties of materials and structures at micro/nano-scale are size-dependent and differ significantly from their behavior at larger scales. Hence, non-classical continuum mechanics such as couple stress theory (CST), nonlocal strain gradient theory (NSGT), and nonlocal elasticity theory (NET) have been applied by the scientific community to analyze the mechanical behavior of micro/nanomaterials due to the difficulties in experimental investigations at the nanoscale and the high computational costs of molecular dynamics (MD) simulations.

Among these theories, the NET introduced by Eringen [1] has been broadly employed in order to include the size effects on the static and dynamic analyses of nanomaterials. For example, Murmu et al. [2] investigated the in-plane magnetic field effect on the vibration of a magnetically sensitive single-layered graphene sheet (SLGS) resting on an elastic foundation via the NET in conjunction with Kirchhoff plate theory. Wang et al. [3] presented an analytical solution for the nonlinear vibration behavior of viscoelastic double layered nanoplates. The problem was tackled through the Kirchhoff plate theory and the multiple scales method. Karličić et al. [4] performed the nonlinear vibration and dynamic instability of single-walled carbon nanotubes (SWCNTs) resting on a Kelvin–Voigt viscoelastic medium under magnetic field using NET and the Euler–Bernoulli beam theory (EBT). Arani and Jalaei [5] employed the Fourier–Laplace transformation method as well as the sinusoidal shear deformation theory (SSDT) to examine the longitudinal magnetic field effect on the transient behavior of a viscoelastic simply-supported orthotropic nanoplate embedded on a visco-Pasternak medium. Kolahchi et al. [6] studied the dynamic buckling of embedded viscoelastic laminated nanoplates on the basis of NET in conjunction with refined zigzag theory. Jalaei et al. [7] studied the dynamic stability of a viscoelastic graphene sheet embedded on viscose matrix with various boundary conditions according to the FSDT and NET by means of the Ritz technique and the Bolotin’s approach. Most recently, Huang et al. [8] applied an analytical solution to perform the dynamic stability of nanobeams exposed to an axially oscillating loading on the basis of the NET and EBT. Additionally, several works in literature have utilized the NET for the study of buckling, bending, vibration, and wave propagation problems of nanobeams/rods [9,10,11,12], nanoshells [13,14,15], and nanoplates [16,17,18,19,20,21,22].

It is recognized that in practical applications, nanostructures or nano reinforced composites often work in environmental conditions. Hence, the mechanical response of these structures at different environmental conditions has been a great subject of investigation among the scientific community. In this regard, Malekzadeh et al. [23] investigated the thermal buckling behavior of arbitrary straight-sided quadrilateral orthotropic nanoplates resting on Pasternak foundation by means of the nonlocal Kirchhoff plate theory and differential quadrature method (DQM). Ansari et al. [24] applied the NET within the framework of the EBT and the Timoshenko beam theory (TBT) to examine the dynamic stability of embedded SWCNTs in a thermal environment. In line with the previous works, Karličić et al. [25] analyzed the thermal vibration and buckling behavior of the embedded multi-layered graphene sheets (MLGSs) according to the NET and Kirchhoff–Love plate theory. Sobhy [26] applied state-space concept to carry out the static bending of Levy-type nanoplates subjected to hygrothermal and mechanical loads based on the NET in conjunction with two-variable plate theory. A DQM numerical approach was applied by Ebrahimi and Hosseini [27] to investigate the nonlinear vibration behavior of bilayer viscoelastic nanoplates in a thermal environment. Ansari and Gholami [28] studied the dynamic instability of multi-walled carbon nanotubes (MWCNTs) embedded on a Winkler medium including thermal effects, while applying a nonlocal TBT. The governing equations were solved via generalized DQM in conjunction with Bolotin’s approach. Wu et al. [29] employed DQM and Bolotin’s approach to explore the dynamic instability behavior of functionally graded graphene-reinforced nanocomposite plates under an axially oscillating loading and uniform temperature variation according to the FSDT. Jouneghani et al. [30] employed NET to study the bending characteristics of porous functionally-graded (FG) Timoshenko nanobeams exposed to hygro-thermo-mechanical loads. Further valuable research of nanostructures in a thermal environment can be found in [31,32,33,34].

Despite many works in literature have focused on the dynamic stability of nanobeams subjected to a thermal loading, there is a general lack of works focusing on the dynamic stability of nanoplates such as graphene sheets in a thermal environment. In most cases, the study of the thermal sensitivity of mechanical behavior of a SLGS is associated with temperature-independent material properties, whereas the experimental outcomes and MD-based simulations have proved the dependence of nanomaterial properties to temperature.

Hence, we propose a numerical study on the dynamic stability of a temperature-dependent SLGS embedded in a temperature-dependent elastomeric medium and subjected to a thermo-mechanical loading. This represents the key point of the present investigation, where the governing equations of motion are derived according to the FSDT in conjunction with the NET. These equations are then converted into a Mathieu–Hill type equation via the Navier’s method, while obtaining the principle unstable region through the application of the Bolotin’s method. A large systematic investigation evaluates the influences of the nonlocal parameter (NP), aspect ratio, temperature variation, foundation type, and static load factor on the dynamic instability region (DIR) of an orthotropic SLGS, as useful in design of nano electro-mechanical systems and micro electro-mechanical systems (NEMS/MEMS).

2. Basic Formulation

A temperature-dependent SLGS is here studied, with length a, width b, and thickness h subjected to a biaxial periodic compressive loading in a thermal environment (Figure 1). The graphene sheet rests on an elastomeric temperature-dependent foundation, which is described by means of the spring constant k_w and the shear layer parameter k_p. Among many possibilities of selecting the plate theory and solution methodology, in this work we select a FSDT, while considering possible shear effects on the structural response of the SLGS despite their atomic dimensions, as already demonstrated in [35]. The selected FSDT is here combined to the NET, whose basic notions are provided in the following overview.

2.1. Overview of NET

The NET introduced by Eringen [1], unlike the classical elasticity theory, states the stress tensor as a nonlocal function of the strain for the whole points in the domain. Hence, for a linear homogenous body, the nonlocal stress tensor at a point x is given as:

$${{\sigma }_{ij}}^{nl}(x)={\displaystyle {\int }_v\alpha (\left|x-x^{\prime }\right|,\tau ){{\sigma }_{ij}}^{l}dV(x^{\prime })},\forall x\in V$$

(1)

where ${\sigma }_{ij}^{nl}$ is the nonlocal stress, ${\sigma }_{ij}^l$ is the local stress tensors, $\alpha (\left|x-x^{\prime }\right|,\tau )$ stands for the nonlocal kernel function. Additionally, $\tau$ is defined as $\tau =e_{0}a/l$, where $a$ and l represent the internal and external characteristic lengths, respectively, and the constant $e_0$ is introduced for adjusting the model according to experimental data or MD simulations. The differential constitutive form of Equation (1) can be written as follows [1,25], including the thermal effects, namely:

$$\left(1-\mu {\nabla }^2\right){\sigma }_{ij}^{nl}={\sigma }_{ij}^l=C_{ijkl}\left({\epsilon }_{kl}-{\alpha }_{kl}\Delta T\right),{\nabla }^2={\partial }^2/\partial x^2+{\partial }^2/\partial y^2.$$

(2)

In Equation (2), $\mu ={(e_{0}l)}^2$ denotes the NP, ${\alpha }_{ij}$ is the thermal expansion coefficient, $\Delta T=T-T_0$ is the temperature variation with respect to the reference temperature $T_0$. In addition, ${\epsilon }_{kl}$ indicates the strain tensor and $C_{ijkl}$ describes the fourth-order stiffness tensor.

Based on Equation (2), the constitutive relations based on the NET for the orthotropic nanoplate in thermal environment reads as follows:

$$\left\{\begin{array}{c}{\sigma }_{xx}^{nl}\\ {\sigma }_{yy}^{nl}\\ {\sigma }_{yz}^{nl}\\ \begin{array}{l}{\sigma }_{xz}^{nl}\\ {\sigma }_{xy}^{nl}\end{array}\end{array}\right\}-\mu {\nabla }^2\left\{\begin{array}{c}{\sigma }_{xx}^{nl}\\ {\sigma }_{yy}^{nl}\\ {\sigma }_{yz}^{nl}\\ \begin{array}{l}{\sigma }_{xz}^{nl}\\ {\sigma }_{xy}^{nl}\end{array}\end{array}\right\}=\left[\begin{array}{ccccc}{C}_{11}& C_{12}& 0& 0& 0\\ C_{21}& C_{22}& 0& 0& 0\\ 0& 0& C_{44}& 0& 0\\ 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& C_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}-{\alpha }_{xx}\Delta T\\ {\epsilon }_{yy}-{\alpha }_{yy}\Delta T\\ {\gamma }_{yz}\\ \begin{array}{l}{\gamma }_{xz}\\ {\gamma }_{xy}\end{array}\end{array}\right\}.$$

(3)

in which the elastic constants of orthotropic materials can be determined by

$$C_{11}=\frac{E_1\left(T\right)}{1-{\nu }_{12}{\nu }_{21}},C_{12}=\frac{{\nu }_{12}{E}_2\left(T\right)}{1-{\nu }_{12}{\nu }_{21}},C_{22}=\frac{E_2\left(T\right)}{1-{\nu }_{12}{\nu }_{21}},C_{44}=G_{23}\left(T\right),C_{55}=G_{13}\left(T\right),C_{66}=G_{12}\left(T\right).$$

(4)

2.2. Strain Displacement Relationships

In what follows we apply the FSDT to formulate the governing equations of the problem, while including the thickness shear deformations and rotary effects. Thus, the displacement components associated to mid-surface $(u_0,v_0,w_0)$, transverse normal rotations $({\varphi }_x,{\varphi }_y)$, and the displacement of an arbitrary point $(u,v,w)$ are related as [36]

$$\begin{array}{l}u\left(x,y,z,t\right)=u_0\left(x,y,t\right)+z{\varphi }_x\left(x,y,t\right),\\ v\left(x,y,z,t\right)=v_0\left(x,y,t\right)+z{\varphi }_y\left(x,y,t\right),\\ w\left(x,y,z,t\right)=w_0\left(x,y,t\right),\end{array}$$

(5)

in which t was the time variable. The nonzero normal and transverse shear strains of the FSDT can be derived by

$$\left\{\begin{array}{l}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\gamma }_{xy}\end{array}\right\}=\left\{\begin{array}{l}{\epsilon }_{xx}^{\left(0\right)}\\ {\epsilon }_{yy}^{\left(0\right)}\\ {\gamma }_{xy}^{\left(0\right)}\end{array}\right\}+z\left\{\begin{array}{l}{\epsilon }_{xx}^{\left(1\right)}\\ {\epsilon }_{yy}^{\left(1\right)}\\ {\gamma }_{xy}^{\left(1\right)}\end{array}\right\},\left\{\begin{array}{l}{\gamma }_{yz}\\ {\gamma }_{xz}\end{array}\right\}=\left\{\begin{array}{l}{\gamma }_{yz}^{\left(0\right)}\\ {\gamma }_{xz}^{\left(0\right)}\end{array}\right\},$$

(6)

where

$$\left\{\begin{array}{l}{\epsilon }_{xx}^{\left(0\right)}\\ {\epsilon }_{yy}^{\left(0\right)}\\ {\gamma }_{xy}^{\left(0\right)}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial u_0}{\partial x}\\ \frac{\partial v_0}{\partial y}\\ \frac{\partial u_0}{\partial y}+\frac{\partial v_0}{\partial x}\end{array}\right\},\left\{\begin{array}{l}{\epsilon }_{xx}^{\left(1\right)}\\ {\epsilon }_{yy}^{\left(1\right)}\\ {\gamma }_{xy}^{\left(1\right)}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial {\varphi }_x}{\partial x}\\ \frac{\partial {\varphi }_y}{\partial y}\\ \frac{\partial {\varphi }_x}{\partial y}+\frac{\partial {\varphi }_y}{\partial x}\end{array}\right\},\left\{\begin{array}{l}{\gamma }_{yz}^{\left(0\right)}\\ {\gamma }_{xz}^{\left(0\right)}\end{array}\right\}=\left\{\begin{array}{c}{\varphi }_y+\frac{\partial w_0}{\partial y}\\ {\varphi }_x+\frac{\partial w_0}{\partial x}\end{array}\right\}.$$

(7)

3. Energy Method

3.1. Kinetic Energy

The kinetic energy for the nanoplate can be computed as

$$\begin{array}{ll}\hfill K& =\frac{\rho }{2}{\displaystyle {\int }_A}{\displaystyle {\int }_{-h/2}^{h/2}\left({\left(\frac{\partial u}{\partial t}\right)}^2+{\left(\frac{\partial v}{\partial t}\right)}^2+{\left(\frac{\partial w}{\partial t}\right)}^2\right)\mathrm{d}A\mathrm{d}z}\hfill \\ & =\frac{1}{2}{\displaystyle {\int }_A\{I_0\left({\left(\frac{\partial u_0}{\partial t}\right)}^2+{\left(\frac{\partial v_0}{\partial t}\right)}^2+{\left(\frac{\partial w_0}{\partial t}\right)}^2\right)}+2I_1\left(\left(\frac{\partial u_0}{\partial t}\frac{\partial {\varphi }_x}{\partial t}\right)+\left(\frac{\partial v_0}{\partial t}\frac{\partial {\varphi }_y}{\partial t}\right)\right)+I_2\left({\left(\frac{\partial {\varphi }_x}{\partial t}\right)}^2+{\left(\frac{\partial {\varphi }_y}{\partial t}\right)}^2\right)\}\mathrm{d}A\end{array}$$

(8)

where $\rho$ denotes the density of the orthotropic graphene sheet and $I_i$ refers to the mass inertias defined by

$$I_i=\rho {\displaystyle {\int }_{-h/2}^{h/2}{z}^i\mathrm{d}z(i=0,1,2)}.$$

(9)

3.2. Strain Energy

The strain energy for the nanoplate could be obtained as

$$U=\frac{1}{2}{\displaystyle {\int }_A}{\displaystyle {\int }_{-h/2}^{h/2}({\sigma }_{xx}{\epsilon }_{xx}+{\sigma }_{yy}{\epsilon }_{yy}+K_s{\sigma }_{xz}{\gamma }_{xz}+K_s{\sigma }_{yz}{\gamma }_{yz}+{\sigma }_{xy}{\gamma }_{xy})\mathrm{d}A\mathrm{d}z}.$$

(10)

In the above equation, K_s stands for the shear correction factor that is here set to 5/6, as typically assumed in the literature for nanoplates and single-layer graphene sheets [20].

The combination of Equations (6), (7), and (10) yields

$$\begin{array}{ll}\hfill U=\frac{1}{2}\underset{A}{\int }& {\displaystyle (N_{xx}\frac{\partial u_0}{\partial x}+M_{xx}\frac{\partial {\varphi }_x}{\partial x}+N_{yy}\frac{\partial v_0}{\partial y}+M_{yy}\frac{\partial {\varphi }_y}{\partial y}+N_{xy}\left(\frac{\partial u_0}{\partial y}+\frac{\partial v_0}{\partial x}\right)}\hfill \\ & +M_{xy}\left(\frac{\partial {\varphi }_x}{\partial y}+\frac{\partial {\varphi }_y}{\partial x}\right)+Q_{xz}\left(\frac{\partial w_0}{\partial x}+{\varphi }_x\right)+Q_{yz}\left(\frac{\partial w_0}{\partial y}+{\varphi }_y\right))\mathrm{d}A\end{array}$$

(11)

where

$$\left\{\begin{array}{l}{N}_{\alpha \beta }\\ M_{\alpha \beta }\end{array}\right\}={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_{\alpha \beta }^{nl}\left\{\begin{array}{l}1\\ z\end{array}\right\}}\mathrm{d}z\alpha =x,y\beta =x,y$$

(12)

$$Q_{\alpha z}=K_s{\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_{\alpha z}^{nl}}\mathrm{d}z\alpha =x,y.$$

(13)

$N_{\alpha \beta }$ denotes the resultant force, $M_{\alpha \beta }$ describes the resultant moment, and $Q_{\alpha z}$ indicates the resultant shear force.

By substitution of Equations (3), (6), and (7) into Equations (12) and (13), the stress resultant-displacement relations are determined as

$$\begin{array}{l}{N}_{xx}-\mu {\nabla }^2\left(N_{xx}\right)=A_{11}\frac{\partial u_0}{\partial x}+A_{12}\frac{\partial v_0}{\partial y}-N_{xx}^T,\\ N_{yy}-\mu {\nabla }^2\left(N_{yy}\right)=A_{12}\frac{\partial u_0}{\partial x}+A_{22}\frac{\partial v_0}{\partial y}-N_{yy}^T,\\ N_{xy}-\mu {\nabla }^2\left(N_{xy}\right)=A_{66}\left(\frac{\partial u_0}{\partial y}+\frac{\partial v_0}{\partial x}\right),\end{array}$$

(14)

$$\begin{array}{l}{M}_{xx}-\mu {\nabla }^2\left(M_{xx}\right)=D_{11}\frac{\partial {\varphi }_x}{\partial x}+D_{12}\frac{\partial {\varphi }_y}{\partial y}-M_{xx}^T,\\ M_{yy}-\mu {\nabla }^2\left(M_{yy}\right)=D_{12}\frac{\partial {\varphi }_x}{\partial x}+D_{22}\frac{\partial {\varphi }_y}{\partial y}-M_{yy}^T,\\ M_{xy}-\mu {\nabla }^2\left(M_{xy}\right)=D_{66}\left(\frac{\partial {\varphi }_x}{\partial y}+\frac{\partial {\varphi }_y}{\partial x}\right),\end{array}$$

(15)

$$\begin{array}{l}{Q}_{yz}-\mu {\nabla }^2\left(Q_{yz}\right)=J_{44}\left({\varphi }_y+\frac{\partial w_0}{\partial y}\right),\\ Q_{xz}-\mu {\nabla }^2\left(Q_{xz}\right)=J_{55}\left({\varphi }_x+\frac{\partial w_0}{\partial x}\right).\end{array}$$

(16)

More specifically, the extensional stiffness $A_{ij}$, the bending stiffness $D_{ij}$, and the shear stiffness $J_{ii}$ of the SLGS, are given as

$$\left(A_{ij},D_{ij}\right)={\displaystyle {\int }_{-h/2}^{h/2}{C}_{ij}}\left(1,z^2\right)\mathrm{d}z(i,j=1,2,6),$$

(17)

$$J_{ii}=K_s{\displaystyle {\int }_{-h/2}^{h/2}{C}_{ii}}\mathrm{d}z=K_{s}h{C}_{ii}(i=4,5).$$

(18)

Additionally, $\left(N_{xx}^T,N_{yy}^T\right)$ and $\left(M_{xx}^T,M_{yy}^T\right)$ denote the thermally induced forces and moments due to a uniform temperature variation, respectively, determined by

$$\left\{\begin{array}{l}{N}_{xx}^T\\ N_{yy}^T\end{array}\right\}={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}\left\{\begin{array}{l}{A}_{11}{\alpha }_{11}+A_{12}{\alpha }_{22}\\ A_{12}{\alpha }_{11}+A_{22}{\alpha }_{22}\end{array}\right\}\Delta T}\mathrm{d}z,$$

(19)

$$\left\{\begin{array}{l}{M}_{xx}^T\\ M_{yy}^T\end{array}\right\}={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}\left\{\begin{array}{l}{A}_{11}{\alpha }_{11}+A_{12}{\alpha }_{22}\\ A_{12}{\alpha }_{11}+A_{22}{\alpha }_{22}\end{array}\right\}\Delta T}z\mathrm{d}z.$$

(20)

3.3. External Work

The external work done due to the temperature-dependent elastomeric medium and axially oscillating loading could be expressed as

$$V=\frac{1}{2}{\displaystyle \underset{A}{\int }\left(-P_{xx}\frac{{\partial }^2{w}_0}{\partial x^2}-P_{yy}\frac{{\partial }^2{w}_0}{\partial y^2}-k_w{w}_0+k_p\left(\frac{{\partial }^2{w}_0}{\partial x^2}+\frac{{\partial }^2{w}_0}{\partial y^2}\right)\right)}{w}_0\mathrm{d}A,$$

(21)

where $P_{xx}=P_{yy}=P$ denotes the axial dynamic forces. Additionally, the elastomeric foundation stiffness $k_w$ could be obtained as [37]

$$k_w=\frac{E_0}{4a\left(1-{\nu }_0^2\right){\left(2-c_1\right)}^2}\left[5-\left(2{\gamma }_1^2+6{\gamma }_1+5\right)\exp \left(-2{\gamma }_1\right)\right],$$

(22)

where

$$c_1=\left({\gamma }_1+2\right)\exp \left(-{\gamma }_1\right),{\gamma }_1=\frac{H_s}{a},E_0=\frac{E_s}{\left(1-{\nu }_s^2\right)},{\nu }_0=\frac{{\nu }_s}{1-{\nu }_s},$$

(23)

while $E_s,{\nu }_s$ and $H_s$ stand for the Young’s modulus, Poisson’s ratio, and depth of the foundation, respectively. In this research, $E_s$ is supposed to be temperature-dependent, whereas ${\nu }_s$ is kept constant.

4. Motion Equations

Hamilton’s variational principle has been utilized in order to determine the equilibrium equations of nanoplate. Accordingly, one may write

$${\displaystyle \underset{0}{\overset{t}{\int }}}\left(\delta U+\delta V-\delta K\right)dt=0.$$

(24)

Inserting Equations (8), (11) and (21) into Equation (24) yields the following equations of motion

$$\delta u_0:\frac{\partial N_{xx}}{\partial x}+\frac{\partial N_{xy}}{\partial y}=I_0\frac{{\partial }^2{u}_0}{\partial t^2}+I_1\frac{{\partial }^2{\varphi }_x}{\partial t^2},$$

(25)

$$\delta v_0:\frac{\partial N_{xy}}{\partial x}+\frac{\partial N_{yy}}{\partial y}=I_0\frac{{\partial }^2{v}_0}{\partial t^2}+I_1\frac{{\partial }^2{\varphi }_y}{\partial t^2},$$

(26)

$$\delta w_0:\frac{\partial Q_{xz}}{\partial x}+\frac{\partial Q_{yz}}{\partial y}-\left(P_{xx}+N_{xx}^T\right)\frac{{\partial }^2{w}_0}{\partial x^2}-\left(P_{yy}+N_{yy}^T\right)\frac{{\partial }^2{w}_0}{\partial y^2}-k_w{w}_0+k_p\left(\frac{{\partial }^2{w}_0}{\partial x^2}+\frac{{\partial }^2{w}_0}{\partial y^2}\right)=I_0\frac{{\partial }^2{w}_0}{\partial t^2},$$

(27)

$$\delta {\varphi }_x:\frac{\partial M_{xx}}{\partial x}+\frac{\partial M_{xy}}{\partial y}-Q_{xz}=I_1\frac{{\partial }^2{u}_0}{\partial t^2}+I_2\frac{{\partial }^2{\varphi }_x}{\partial t^2},$$

(28)

$$\delta {\varphi }_y:\frac{\partial M_{xy}}{\partial x}+\frac{\partial M_{yy}}{\partial y}-Q_{yz}=I_1\frac{{\partial }^2{v}_0}{\partial t^2}+I_2\frac{{\partial }^2{\varphi }_y}{\partial t^2}.$$

(29)

Finally, by substituting Equations (14)–(16) into Equations (25)–(29), the nonlocal governing equations of orthotropic graphene sheet resting on an elastomeric medium including the thermal effect can be determined as follows:

$$A_{11}\left(\frac{{\partial }^2{u}_0}{\partial x^2}\right)+\left(A_{12}+A_{66}\right)\left(\frac{{\partial }^2{v}_0}{\partial x\partial y}\right)+A_{66}\left(\frac{{\partial }^2{u}_0}{\partial y^2}\right)+\left(1-\mu {\nabla }^2\right)\left[-I_0\frac{{\partial }^2{u}_0}{\partial t^2}-I_1\frac{{\partial }^2{\varphi }_x}{\partial t^2}\right]=0,$$

(30)

$$\left(A_{12}+A_{66}\right)\left(\frac{{\partial }^2{u}_0}{\partial x\partial y}\right)+A_{66}\left(\frac{{\partial }^2{v}_0}{\partial x^2}\right)+A_{22}\left(\frac{{\partial }^2{v}_0}{\partial y^2}\right)+\left(1-\mu {\nabla }^2\right)\left[-I_0\frac{{\partial }^2{v}_0}{\partial t^2}-I_1\frac{{\partial }^2{\varphi }_y}{\partial t^2\partial y}\right]=0,$$

(31)

$$\begin{array}{cc}\hfill J_{55}\left(\frac{{\partial }^2{w}_0}{\partial x^2}+\frac{\partial {\varphi }_x}{\partial x}\right)& +J_{44}\left(\frac{{\partial }^2{w}_0}{\partial y^2}+\frac{\partial {\varphi }_y}{\partial y}\right)\hfill \\ & +\left(1-\mu {\nabla }^2\right)\left(\begin{array}{l}-k_w{w}_0+k_p\left(\frac{{\partial }^2{w}_0}{\partial x^2}+\frac{{\partial }^2{w}_0}{\partial y^2}\right)-\left(P_{xx}+N_{xx}^T\right)\frac{{\partial }^2{w}_0}{\partial x^2}\\ -\left(P_{yy}+N_{yy}^T\right)\frac{{\partial }^2{w}_0}{\partial y^2}-I_0\frac{{\partial }^2{w}_0}{\partial t^2}\end{array}\right)=0,\hfill \end{array}$$

(32)

$$\begin{array}{l}{D}_{11}\left(\frac{{\partial }^2{\varphi }_x}{\partial x^2}\right)+D_{12}\left(\frac{{\partial }^2{\varphi }_y}{\partial x\partial y}\right)+D_{66}\left(\frac{{\partial }^2{\varphi }_x}{\partial y^2}+\frac{{\partial }^2{\varphi }_y}{\partial x\partial y}\right)-J_{55}\left(\frac{\partial w_0}{\partial x}+{\varphi }_x\right)\\ +\left(1-\mu {\nabla }^2\right)\left(-I_1\frac{{\partial }^2{u}_0}{\partial t^2}-I_2\frac{{\partial }^2{\varphi }_x}{\partial t^2}\right)=0,\end{array}$$

(33)

$$\begin{array}{l}{D}_{12}\left(\frac{{\partial }^2{\varphi }_x}{\partial x\partial y}\right)+D_{22}\left(\frac{{\partial }^2{\varphi }_y}{\partial y^2}\right)+D_{66}\left(\frac{{\partial }^2{\varphi }_x}{\partial x\partial y}+\frac{{\partial }^2{\varphi }_y}{\partial x^2}\right)-J_{44}\left(\frac{\partial w_0}{\partial y}+{\varphi }_y\right)\\ +\left(1-\mu {\nabla }^2\right)\left(-I_1\frac{{\partial }^2{v}_0}{\partial t^2}-I_2\frac{{\partial }^2{\varphi }_y}{\partial t^2}\right)=0.\end{array}$$

(34)

5. Solution Procedure

The governing differential equations are solved analytically via Navier’s approach for simply-supported graphene sheet in this research. According to the selected approach, the mid-plane displacement and rotation field are solved by means of the Fourier series as

$$u_0(x,y,t)={\displaystyle \sum _{m=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }{U}_{mn}(t)\cos (\alpha x)\sin (\beta y)}},$$

(35)

$$v_0\left(x,y,t\right)={\displaystyle \sum _{m=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }{V}_{mn}}}(t)\sin \left(\alpha x\right)\cos \left(\beta y\right),$$

(36)

$$w_0\left(x,y,t\right)={\displaystyle \sum _{m=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }{W}_{mn}}}(t)\sin \left(\alpha x\right)\sin \left(\beta y\right),$$

(37)

$${\varphi }_x\left(x,y,t\right)={\displaystyle \sum _{m=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }{X}_{mn}}}(t)\cos \left(\alpha x\right)\sin \left(\beta y\right),$$

(38)

$${\varphi }_y\left(x,y,t\right)={\displaystyle \sum _{m=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }{Y}_{mn}}}(t)\sin \left(\alpha x\right)\cos \left(\beta y\right),$$

(39)

where $\alpha =\frac{m\pi }{a},\beta =\frac{n\pi }{b},$ and m is the half wave number in the x direction and n is the half wave number in the y direction.

Finally, substituting Equations (35)–(39) into Equations (30)–(34), leads to the following matrix form, namely

$$M\ddot{d}+\left(K_L-P\left(t\right)K_G\right)d=0,$$

(40)

in which $d={\left\{UVWXY\right\}}^T$ denotes the displacement vector. Additionally, M refers to the mass, K_L describes the stiffness, and K_G stands for the geometric stiffness matrices.

We consider a periodic compressive load which includes a static and dynamic part, as follows:

$$P(t)=P_s+P_d\cos \left(\omega t\right)=\alpha P_{cr}+\beta P_{cr}\cos \left(\omega t\right).$$

(41)

Here, $P_s$ refers to the static force component and $P_d$ indicates the dynamic force component, $P_{cr}$ shows the static buckling load and $\omega$ describes the pulsation frequency. Additionally, $\alpha$ is the static load factor and $\beta$ is the dynamic load factor. The substitution of Equation (41) into Equation (40) read

$$M\ddot{d}+\left(K_L-\left(\alpha +\beta \cos \left(\omega t\right)\right)P_{cr}{K}_G\right)d=0.$$

(42)

Equation (42) is a Mathieu–Hill type equation discussing the dynamic instability behavior of graphene sheet under a periodic compressive load. The Bolotin’s method [38] is employed to determine the boundaries of the DIR. Based on this method, the displacement vector of d can be written in the following series with period 2T as

$$d={\displaystyle \sum _{k=1,3,\dots }^{\infty }\left[a_k\sin \frac{k\omega t}{2}+b_k\cos \frac{k\omega t}{2}\right]},$$

(43)

in which $a_k$ and $b_k$ denote arbitrary constant vectors. On the basis of Bolotin’s approach [38], the first approximation with $k=1$ could be used to compute the instability boundary. In this case, by setting Equation (43) into Equation (42) and considering a null value for the coefficients of the harmonic functions and the sum of the constant terms to zero, we get to the following relation:

$$\left|\left(K_L-P_{cr}\left(\alpha \pm \frac{\beta }{2}\right)K_G-M\frac{{\omega }^2}{4}\right)\right|=0.$$

(44)

At the end, we can compute the critical static buckling load, assuming $\beta =\omega =0$. On the basis of the eigenvalue problem, the variation of $\omega$ versus $\beta$ can be plotted as instability regions of the nanoplate to the periodic load including thermal effect.

6. Numerical Examples

This section is devoted to the numerical study on the dynamic instability of orthotropic temperature-dependent SLGS resting on temperature-dependent elastomeric foundation undergoing biaxial oscillating loading in a thermal environment. Here, two types of zigzag SLGSs are considered with two different values of aspect ratio. The material properties of the SLGS are highly dependent to the temperature. Shen et al. [39] performed a MD simulation to evaluate these properties at three different temperatures, as listed in

Table 1. Thermo-mechanical properties of the material [39].

(K)	E1 (TPa)	E2 (TPa)	G12 (TPa)	${\mathit{\nu }}_{\mathbf{12}}$	${\mathit{\alpha }}_{11}({10}^{-6}/\mathbf{K})$	${\mathit{\alpha }}_{22}({10}^{-6}/\mathbf{K})$
Zigzag sheet I: a = 4.855 nm, b = 4.888 nm, h = 0.154 nm
300	1.987	1.974	0.857	0.205	2.1	1.9
500	1.974	1.968	0.870	0.205	2.3	2.0
700	1.961	1.948	0.870	0.205	2.4	2.1
Zigzag sheet II: a = 9.496 nm, b = 4.877 nm, h = 0.145 nm
300	2.145	2.097	0.938	0.223	1.7	1.5
500	2.103	2.055	0.959	0.223	2.0	1.7
700	2.069	2.014	0.959	0.223	2.1	2.0

. In addition, a polydimethylsiloxane (PDMS) is here selected for the elastomeric foundation, while assuming the following material properties: a Poisson’s ratio ${\nu }_s=0.48$ and a Young’s modulus $E_s=\left(3.22-0.0034T\right)GPa$, where $T=T_0+\Delta T$ and $T_0=300K$ refers to the room temperature [37].

6.1. Validation of Results

A preliminary validation of the proposed mathematical formulation is performed against some results available in literature, see

Table 2. Dimensionless natural frequencies of the square graphene sheet for different nonlocal and foundation parameters (a/h = 10).

μ (nm2)	KW	KP	Sobhy [40]	Ebrahimi and Barati [41]	Present	Error (%)
0	0	0	1.93861	1.93861	1.93851	0.0051
	100	0	2.18396	2.18396	2.18386	0.0045
	100	20	2.96017	-	2.96010	0.0023
1	0	0	1.17816	1.17816	1.17809	0.0059
	100	0	1.54903	1.54903	1.54898	0.0032
	100	20	2.52831	-	2.52827	0.0015
2	0	0	0.92261	0.92261	0.92257	0.0043
	100	0	1.36479	1.36479	1.36475	0.0029
	100	20	2.41979	-	2.41976	0.0012
3	0	0	0.78347	0.78347	0.78343	0.0051
	100	0	1.27485	1.27485	1.27482	0.0023
	100	20	2.37020	-	2.37018	0.0008

and

Table 3. Comparative evlauation of the critical buckling load (nN/nm) for isotropic square graphene sheet under biaxial compression loads.

a (nm)	MD Results [42]	DQM [43]	FSM [19]	Present
4.99	1.0837	1.0749	-	1.0711
8.08	0.6536	0.6523	0.6592	0.6514
10.77	0.4331	0.4356	0.4387	0.4353
14.65	0.2609	0.2645	0.2656	0.2643
18.51	0.1714	0.1751	0.1757	0.1751
26.22	0.0889	0.0917	0.0918	0.0917
33.85	0.0554	0.0561	0.0562	0.0561
37.81	0.0449	0.0453	0.0453	0.0452
45.66	0.0315	0.0313	0.0313	0.0313

. More specifically, Table 2 summarizes the results in terms of dimensionless natural frequencies for a simply-supported square graphene sheet, as computed by the present method, against predictions by Sobhy [40] based on SSDT, and or those ones by Ebrahimi and Barati [41] according to the higher-order refined plate theory and the DQM approach. A parametric investigation analyzes the variation of the NPs and foundation constants. The outcomes in Table 2 clearly show that the natural frequencies gradually decrease for an increased value of μ and increased for increasing values of K_W and K_P. The very good agreement between our results and predictions from literature, clearly confirms the reliability of the present approach to capture the response of the problem, where few negligible differences are only related to the different selected solution methodology and plate theory.

As another attempt for validation, we compute the critical buckling loads for isotropic square graphene sheet under biaxial load with different lengths, whose results are summarized and compared in Table 3 with the ones based on an MD approach [42], a DQM approach [43], or a refined strip method (RFSM) [19]. The gradual increase of the geometrical dimension a yields a meaningful reduction of the critical loading, whereas the very good agreement between our results and those reported in literature, confirms once again the accuracy of the proposed method.

6.2. Parametric Studies and Discussions

In the following examples, the DIR is traced as a dimensionless excitation frequency $\mathrm{\Omega }=\omega b^2\sqrt{\rho h/D_{11}}\mathrm{in}\mathrm{which}{D}_{11}=E_1{h}^3/12(1-{\nu }_{12}{\nu }_{21})$ against the dynamic load factor $\beta$. It was supposed, unless otherwise explained, that $a/b=1$(zigzag sheet I) and $\mu =0.5{\mathrm{nm}}^2$.

Figure 2 investigates the influence of the temperature change on the dynamic instability of orthotropic graphene sheet without elastomeric foundation. As clearly visible in Figure 2, an increase in temperature change reduces the excitation frequencies and shifts the origins of the DIR to the left side and also decreases gradually its width. The temperature rise leads to an increased compressive prestress in the graphene sheet, thus weakening the nanoplate stiffness.

As mentioned before, in this study we analyze the effect of the temperature on the mechanical properties of both the nanostructure and elastic foundation. Thus, Figure 3a,b plot the DIR of the SLGS resting on the elastomeric foundation under the twofold assumption of temperature-independent and temperature-dependent material properties. A comparative evaluation of the results is illustrated in Figure 3a,b, for two different values of temperature, i.e., $T=500\mathrm{K}$ and $T=700\mathrm{K}$, respectively. As visible in Figure 3a,b, when temperature-dependent properties are used, the origins of the instability regions occur at smaller excitation frequencies. The difference between the two cases becomes even more pronounced as the temperature is enhanced. This means that it is significant to account for the effect of temperature on the material properties in order to obtain more reliable results.

The surrounding elastic foundation effect on the DIR of SLGS is displayed in Figure 4a,b, for T = 300 K and T = 700 K, respectively. As visible in Figure 4a,b, the presence of an elastic foundation increases the excitation frequency and hence, the origin of unstable region moves to the right side. The existence of an elastomeric foundation magnifies the nanoplate stiffness which yields an increased stability of the nanoplate. In addition, by considering an elastic foundation, the instability region of the graphene sheet tends to become wider, especially at increased levels of temperature. Moreover, it can be observed that the effect of the elastomeric Pasternak medium $({\mathrm{K}}_{\mathrm{W}}\ne 0\mathrm{and}{\mathrm{K}}_{\mathrm{P}}\ne 0)$ is higher than the elastomeric Winkler $({\mathrm{K}}_{\mathrm{W}}\ne 0\mathrm{and}{\mathrm{K}}_{\mathrm{P}}=0)$ one, on the DIR of the SLGS. The elastomeric Winkler foundation is able to describe only the normal loading condition of the elastomeric foundation, whereas the Pasternak model considers normal stresses as well as transverse shear loading conditions. From these figures, it is also worth observing that the foundation effects become even more significant for increased temperatures.

Figure 5a,b illustrates the aspect ratio $a/b$ effect on the unstable region of graphene sheets embedded in the elastomeric foundation for T = 300 K and T = 700 K, respectively. In this example, two different values of aspect ratio (i.e., $a/b\simeq 1$ (zigzag sheet I) and $a/b\simeq 2$ (zigzag sheet II)) are considered [39]. This figure indicates that the excitation frequency reduces and the dynamic instability moves to the left side as $a/b$ increases. This trend is almost predictable, since the nanoplate becomes softer for increasing aspect ratios, with a global reduction of the structural stiffness. Moreover, the width of the instability region increases with a decrease in the slenderness ratio.

The NP effect on the DIR of embedded SLGS is displayed in Figure 6a,b, for T = 300 K and T = 700 K, respectively. The dimensionless pulsation frequency seems to decrease as the NP becomes higher, since the DIR of graphene sheet forms at lower excitation frequencies. Such phenomenon is owing to the fact that NP indicates the softening effect on stiffness which leads to lower pulsation frequency. Furthermore, comparing these figures exhibits that an increase in temperature change increases the effect of NP on the unstable region of the nanoplate.

Finally, Figure 7a,b demonstrates the static load factor effect on the unstable region of the SLGS resting on the elastomeric foundation for T = 300 K and T = 500 K, respectively. As can be concluded, with an increase in the static load factor the excitation frequency decreases and therefore, dynamic instability shifted to the left side due to the reduction in the stiffness of the nanoplate. Furthermore, the results indicate that the DIR becomes wider for a higher static load factor. Any static force component is present, in the particular case when $\alpha =0$.

7. Conclusions

This study investigated the dynamic instability of temperature-dependent SLGS embedded in an elastic medium in thermal environment. This was modelled by means of a temperature-dependent elastomeric foundation. Based on FSDT within the framework of Eringen’s differential constitutive model, the equations of motion were here derived from the Hamilton’s variational principle as well as the energy method. The Navier’s solution in conjunction with the Bolotin’s approach were here applied to calculate the unstable region for the graphene sheet. A large numerical investigation examined the influence of different parameters such as the NP, temperature change, aspect ratio, static load factor, and foundation type on the dynamic stability behavior. Based on the numerical results, it was found that the elastic foundation increases the width of the instability region and the origin of the unstable region shift to the right side which leads to higher excitation frequencies. An increased static load factor and/or temperature weakens the nanoplate stiffness and leads to a lower pulsating frequency, while moving the DIR to the left side. It seemed also that the DIR occurs at smaller values of excitation frequency, when temperature-dependent mechanical properties were assumed for both the SLGS and the foundation, compared to the assumption of constant properties.