# Vibration Analysis of Bilayered Graphene Sheets for Building Materials in Thermal Environments Based on the Element-Free Method

Authors: Zhengtian Wu; Yang Zhang; Fuyuan Hu; Qing Gao; Xinyin Xu; Ronghao Zheng
Journal: Journal of Nanomaterials
 (2018)

---

## Abstract

Graphene sheets are widely applied due to their unique and highly valuable properties, such as energy conservation materials in buildings. In this paper, we investigate the vibration behavior of double layer graphene sheets (DLGSs) in thermal environments which helps probe into the mechanism of energy conservation of graphene sheets in building materials. The nonlocal elastic theory and classical plate theory (CLPT) are used to derive the governing equations. The element-free method is employed to analyze the vibration behaviors of DLGSs embedded in an elastic medium. The accuracy of the solutions in this study is demonstrated by comparison with published results in the literature. Furthermore, a number of numerical results are presented to investigate the effects of various parameters (boundary conditions, nonlocal parameter, aspect ratio, side length, elastic foundation parameter, and temperature) on the frequency of DLGSs.


---

## Body

## 1. Introduction

In 2004, Novoselov and Geim in [1], physicists at the University of Manchester, successfully stripped off graphene sheets from the graphite. Graphene sheets (GSs) have carbon atomic structure of the monolayer or multilayer atomic thickness. DLGSs in which the interlayer width is 0.34 nm [2] have specific chemical, mechanical, thermal, and electronic properties [3–10] and have been applied as constituent in nanoelectromechanical systems (NEMS) and microelectromechanical systems (MEMS). DLGSs used in these applications generally rest on the elastic medium rather than being suspended. Due to the increasing applications of graphene sheets, understanding the vibration of DLGSs embedded in the elastic media is essential to optimize the design and manufacture in engineering.Applications of DLGSs depend on clear knowledge of their mechanical behavior. In recent years, world-wide researchers have engaged in the study of GSs using experimental, numerical, and theoretical tools. However, it is of difficulty to carry out experimental research due to the very small dimensions of the nanostructure elements [11]. Thus, theoretical research is becoming increasingly important for studying nanostructures. The theoretical tools can be divided into three categories [12]: atomistic modeling, continuum mechanics (CM), and the hybrid method. Generally, the atomistic modeling consists of classic tight-binding molecular dynamics (TBMD) and molecular dynamics (MD). Mirparizi and Aski [13] studied the interlayer shear effect on vibrational behavior for DLGSs using the MD simulation. Gajbhiye and Singh [14] studied the nonlinear frequency response of SLGSs via the atomistic finite element method (FEM). The MD simulation is still the most sophisticated and reliable approach which provides deep insight into nanomaterials and their properties. However, MD simulation is limited to the analysis of GSs with few atoms and time scales because of high computational cost. Thus, unlike MD simulations, the CM has a lower computational cost [15] and has been adopted by numerous researchers to deal with nanoscale structures. Therefore, some relatively large nanostructures can be analyzed by this method [16, 17]. The classical continuum approach is however a kind of scale-free theory which cannot account for the small scale effect arising in the nanostructures. Therefore, the mechanical behavior of nanostructures cannot be successfully captured by the classical continuum approach. It behooves us to modify the classical continuum model which can describe precise mechanical behaviors of the nanostructures. Some improved continuum methods such as couple stress theory [18], surface elasticity theory [19], nonlocal theory, and strain gradient theory [20] were used in the simulation of the size-dependent materials. It is found that the nonlocal theory has provided good matching results compared with lattice dynamics. Recent literature indicates that the theory of nonlocal elasticity is being increasingly used [21–23]. Peddieson et al. [24] were the first ones who adopted the nonlocal elasticity theory to analyze the static deformation behavior of Euler-Bernoulli nanobeams. Duan and Wang [25] firstly formulated the nonlocal continuum plate model to study the small scale effect in the bending regarding circular nanoplates. Dastjerdi and Jabbarzadeh [26] studied the multilayer orthotropic circular GSs nonlinear bending using nonlocal elasticity theory. Liew et al. [27] investigated the effect of polymer matrix on the vibration response of GSs. Zhang et al. [12] studied the free vibration behaviors of DLGSs employing nonlocal elasticity theory. Zhang et al. [28] analyzed the vibration of quadrilateral graphene sheets subjected to an in-plane magnet using nonlocal elasticity theory. Pradhan and Phadikar [29] conducted the vibration analysis of double-layered GSs resting on a polymer matrix utilizing the nonlocal continuum mechanics. Malekzadeh et al. [21] investigated the small scale effect on the thermal buckling characteristic of orthotropic arbitrary straight-sided quadrilateral nanoplates embedded in an elastic medium. Asemi et al. [30] presented an axisymmetric buckling analysis of circular SLGS using Eringen’s theory. Additionally, Asemi et al. [31] used nonlocal theory of Eringen to simulate postbuckling response of orthotropic SLGS.Numerical techniques are important for theoretical analysis. Various numerical tools have been utilized to analyze the mechanical behaviors of nanostructure [32, 33]. Ansari et al. [34] proposed a nonlocal FEM model to simulate the vibration behaviors of multilayered GSs based on the Mindlin plate theory. As a numerical technique, the meshless method has been successfully applied in solving engineering problems. Because of its advantage of not relying on elements, it has good adaptability [35–39].In this work, we study the vibration behavior of DLGSs resting on a Pasternak elastic medium employing the nonlocal theory and established the meshless framework to obtain the solution. The effects of nonlocal parameter, temperature, aspect ratio, side length, stacking types, and foundation parameters on the vibration frequency of DLGSs are examined. Additionally, the accuracy of the solutions in this study is demonstrated through comparison with published results in the literature.


## 2. Nonlocal Elasticity Theory and Model of Elastic Medium

### 2.1. Nonlocal Elasticity Theory

Eringen [36] proposed the nonlocal theory through experiment and the theoretical analysis in which the stress depends on the strain of all the points in the domain. Therefore, it can reflect the small scale effect. The most widely applied constitutive relation of the differential form is written as(1)1-e0a2∇2σijx=Cijklεklx,in which e0a is the nonlocal parameter and ∇2 is the Laplacian operator defined by ∇2=∂2/∂x2+∂2/∂y2.


### 2.2. Model of Elastic Medium

In nanoscale devices, vibration behaviors can be observed, which affect the performance significantly. As mentioned above, graphene based products are usually resting on an elastic medium. In the present study, the Pasternak-type foundation [40] is adopted. Pasternak-type foundation model is modeled with two parameters. The mathematical expressions of Pasternak-type foundation are written as (2)qpasternak=kww-kg∇2w,where kw is the Winker modulus and kg is the shear modulus.


### 2.3. Temperature Field Distribution Type

In this paper, the temperature field can be divided into three types as linear, nonlinear, and uniform. The temperature changes along the thickness of the GSs. The temperature change obeys the following law [22]:(3)ΔTz=z+h/2zφT,-h2≤z≤h2,where φ is the temperature distribution coefficient. When φ>1, there is a nonlinear temperature distribution along the thickness of the GSs; for φ=1 we have a linear temperature distribution, for φ=0 the temperature distribution is uniform, and the temperature change is T=Tb-Tr, where Tb is the temperature value when the buckling occurs, while Tr is the referent temperature. The thermal load caused by temperature change satisfies the following formulation:(4)NT=-αE1-υ∫-h/2h/2z+h/2zφTdz=-αhE1-υ1+φT.


## 2.1. Nonlocal Elasticity Theory

Eringen [36] proposed the nonlocal theory through experiment and the theoretical analysis in which the stress depends on the strain of all the points in the domain. Therefore, it can reflect the small scale effect. The most widely applied constitutive relation of the differential form is written as(1)1-e0a2∇2σijx=Cijklεklx,in which e0a is the nonlocal parameter and ∇2 is the Laplacian operator defined by ∇2=∂2/∂x2+∂2/∂y2.


## 2.2. Model of Elastic Medium

In nanoscale devices, vibration behaviors can be observed, which affect the performance significantly. As mentioned above, graphene based products are usually resting on an elastic medium. In the present study, the Pasternak-type foundation [40] is adopted. Pasternak-type foundation model is modeled with two parameters. The mathematical expressions of Pasternak-type foundation are written as (2)qpasternak=kww-kg∇2w,where kw is the Winker modulus and kg is the shear modulus.


## 2.3. Temperature Field Distribution Type

In this paper, the temperature field can be divided into three types as linear, nonlinear, and uniform. The temperature changes along the thickness of the GSs. The temperature change obeys the following law [22]:(3)ΔTz=z+h/2zφT,-h2≤z≤h2,where φ is the temperature distribution coefficient. When φ>1, there is a nonlinear temperature distribution along the thickness of the GSs; for φ=1 we have a linear temperature distribution, for φ=0 the temperature distribution is uniform, and the temperature change is T=Tb-Tr, where Tb is the temperature value when the buckling occurs, while Tr is the referent temperature. The thermal load caused by temperature change satisfies the following formulation:(4)NT=-αE1-υ∫-h/2h/2z+h/2zφTdz=-αhE1-υ1+φT.


## 3. Theoretical Approach

### 3.1. The Van der Waals (vdW) Force Expression

In DLGSs, the interaction force between layers is the vdW force which is a nonbonded interaction. By neglecting the nonlinear terms, the vdW’s interaction can be expressed as [41](5)qij=cijwi-wj,i,j=1,2in which cij represents the vdW coefficient. Generally, DLGSs consist of two types in terms of the relative positions of the upper and lower layers. The one with the upper graphene layer stacking directly on lower layer is called AA-type whereas the other is called AB-type [41]. Figure 1 depicts these two categories of DLGSs. According to [41], the vdW coefficient is(6)cij=-27.3546Gpa/nmforAA-type,cij=-146.8247Gpa/nmforAB-type.If there is no special description, the vdW coefficient is taken as AB-type in the following simulations.Figure 1
Two stacking types of DLGSs.


### 3.2. Governing Equations

The DLGSs resting on elastic foundation are demonstrated in Figure2 in which they are in the thermal condition. The external foundation is modeled as Pasternak-type. Lx and Ly, respectively, denote the length and width of GSs. The thickness of graphene is h. Transverse displacements of the top-GSs and bottom-GSs for the buckling analysis are represented by w1(x,y,t) and w2(x,y,t), respectively.Figure 2
DLGSs resting on Pasternak-type elastic foundation under thermal environment.The displacement field of DLGSs according to CLPT can be written as(7)ux,y,z;t=-z∂w1x,y;t∂xvx,y,z;t=-z∂w1x,y;t∂ywx,y,z;t=w1x,y;tin which (u,v,w), respectively, denote the displacements along the axes (x,y,z). w1 is the displacement along z-axis of a mid-plane material point. According to the strain-displacement relations(8)ε1xx=-z∂2w1x,y;t∂x2(9)ε1yy=-z∂2w1x,y;t∂y2(10)γ1xy=-2z∂2w1x,y;t∂x∂y(11)ε1zz=γ1xz=γ1yz=0the constitutive relationship is given by(12)1-e0a2∇2σ1xxσ1yyσ1xy=E1-υ2υE1-υ20υE1-υ2E1-υ2000Gε1xx-αΔTε1yy-αΔTε1xy,where G=E/21+υ.According to (10), the nonlocal moments can be written as(13)1-e0a2∇2M1xx=-D∂2w1x,y;t∂x2+υ∂2w1x,y;t∂y2,(14)1-e0a2∇2M1yy=-D∂2w1x,y;t∂y2+υ∂2w1x,y;t∂x2,(15)1-e0a2∇2M1xy=-D1-υ∂2w1x,y;t∂x∂y,where (M1xx,M1yy,Mxy)=∫(σ1xx,σ1yy,σ1xy)zdz. D is the bending stiffness of the plate and D=Eh3/12(1-υ2).With the principle of virtual work, the equilibrium equation can be derived as(16)-M1xx,xx-M1yy,yy-2M1xy,xy-q12+kww1-kg∇2w1-NT∇2w1-I2w¨1,xx+w¨1,yy+I0w¨1=0.From (11)–(13), we can get(17)-D∇4w1=1-e0a2∇2M1xx,xx+M1yy,yy+2M1xy,xy.Substituting (17) into (16), the governing equation can be written as(18)D∇4w1+1-e0a2∇2-q12+kww1-kg∇2w1-NT∇2w1-I2w¨1,xx+w¨1,yy+I0w¨1.Similarly, we can obtain the governing equation of bottom-GSs(19)D∇4w2+1-e0a2∇2-q21+kww2-kg∇2w2-NT∇2w2-I2w¨2,xx+w¨2,yy+I0w¨2.


## 3.1. The Van der Waals (vdW) Force Expression

In DLGSs, the interaction force between layers is the vdW force which is a nonbonded interaction. By neglecting the nonlinear terms, the vdW’s interaction can be expressed as [41](5)qij=cijwi-wj,i,j=1,2in which cij represents the vdW coefficient. Generally, DLGSs consist of two types in terms of the relative positions of the upper and lower layers. The one with the upper graphene layer stacking directly on lower layer is called AA-type whereas the other is called AB-type [41]. Figure 1 depicts these two categories of DLGSs. According to [41], the vdW coefficient is(6)cij=-27.3546Gpa/nmforAA-type,cij=-146.8247Gpa/nmforAB-type.If there is no special description, the vdW coefficient is taken as AB-type in the following simulations.Figure 1
Two stacking types of DLGSs.


## 3.2. Governing Equations

The DLGSs resting on elastic foundation are demonstrated in Figure2 in which they are in the thermal condition. The external foundation is modeled as Pasternak-type. Lx and Ly, respectively, denote the length and width of GSs. The thickness of graphene is h. Transverse displacements of the top-GSs and bottom-GSs for the buckling analysis are represented by w1(x,y,t) and w2(x,y,t), respectively.Figure 2
DLGSs resting on Pasternak-type elastic foundation under thermal environment.The displacement field of DLGSs according to CLPT can be written as(7)ux,y,z;t=-z∂w1x,y;t∂xvx,y,z;t=-z∂w1x,y;t∂ywx,y,z;t=w1x,y;tin which (u,v,w), respectively, denote the displacements along the axes (x,y,z). w1 is the displacement along z-axis of a mid-plane material point. According to the strain-displacement relations(8)ε1xx=-z∂2w1x,y;t∂x2(9)ε1yy=-z∂2w1x,y;t∂y2(10)γ1xy=-2z∂2w1x,y;t∂x∂y(11)ε1zz=γ1xz=γ1yz=0the constitutive relationship is given by(12)1-e0a2∇2σ1xxσ1yyσ1xy=E1-υ2υE1-υ20υE1-υ2E1-υ2000Gε1xx-αΔTε1yy-αΔTε1xy,where G=E/21+υ.According to (10), the nonlocal moments can be written as(13)1-e0a2∇2M1xx=-D∂2w1x,y;t∂x2+υ∂2w1x,y;t∂y2,(14)1-e0a2∇2M1yy=-D∂2w1x,y;t∂y2+υ∂2w1x,y;t∂x2,(15)1-e0a2∇2M1xy=-D1-υ∂2w1x,y;t∂x∂y,where (M1xx,M1yy,Mxy)=∫(σ1xx,σ1yy,σ1xy)zdz. D is the bending stiffness of the plate and D=Eh3/12(1-υ2).With the principle of virtual work, the equilibrium equation can be derived as(16)-M1xx,xx-M1yy,yy-2M1xy,xy-q12+kww1-kg∇2w1-NT∇2w1-I2w¨1,xx+w¨1,yy+I0w¨1=0.From (11)–(13), we can get(17)-D∇4w1=1-e0a2∇2M1xx,xx+M1yy,yy+2M1xy,xy.Substituting (17) into (16), the governing equation can be written as(18)D∇4w1+1-e0a2∇2-q12+kww1-kg∇2w1-NT∇2w1-I2w¨1,xx+w¨1,yy+I0w¨1.Similarly, we can obtain the governing equation of bottom-GSs(19)D∇4w2+1-e0a2∇2-q21+kww2-kg∇2w2-NT∇2w2-I2w¨2,xx+w¨2,yy+I0w¨2.


## 4. Discretized System Equations

### 4.1. Weak Form

In order to get the weak form, the equilibrium equation is multiplied by the virtual displacementδw. Therefore, we obtain(20)∫Ωδw1D∇4w1+1-e0a2∇2-q12+kww1-NT∇2w1-kg∇2w1+I0w¨1-I2w¨1,xxtt+w¨1,yyttdxdy=0.Equation (20) can be integrated by parts as(21)
D
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
-
c
∫
Ω
w
1
δ
w
1
-
w
2
δ
w
1
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
k
w
∫
Ω
w
1
δ
w
1
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
e
0
a
2
k
w
∫
Ω
w
1
,
x
δ
w
1
,
x
+
w
1
,
y
δ
w
1
,
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
N
T
∫
Ω
w
1
,
x
δ
w
1
,
x
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
-
c
e
0
a
2
∫
Ω
w
1
,
x
δ
w
1
,
x
+
w
1
,
y
δ
w
1
,
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
e
0
a
2
N
T
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
x
y
δ
w
1
,
x
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
c
e
0
a
2
∫
Ω
w
2
,
x
δ
w
1
,
x
+
w
2
,
y
δ
w
1
,
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
k
g
∫
Ω
w
1
,
x
δ
w
1
,
x
+
w
1
,
y
δ
w
1
,
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
e
0
a
2
∫
Ω
I
0
δ
w
1
,
x
w
1
,
x
t
t
+
δ
w
1
,
y
w
1
,
y
t
t
+
I
2
δ
w
1
,
x
x
w
1
,
x
x
t
t
+
2
δ
w
1
,
x
y
w
1
,
x
y
t
t
+
δ
w
1
,
y
y
w
1
,
y
y
t
t
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
e
0
a
2
N
T
∫
Ω
w
1
,
x
y
δ
w
1
,
x
y
+
w
1
,
y
y
δ
w
1
,
y
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
N
T
∫
Ω
w
1
,
y
δ
w
1
,
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
k
g
e
0
a
2
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
∫
Ω
I
0
δ
w
1
w
1
,
t
t
+
I
2
δ
w
1
,
x
w
1
,
x
t
t
+
δ
w
1
,
y
w
1
,
y
t
t
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
-
e
0
a
2
k
w
∫
Γ
w
1
,
x
δ
w
1
n
x
-
w
1
,
y
δ
w
1
n
y
d
s
-
N
T
∫
Γ
w
1
,
x
δ
w
1
n
x
d
s
+
N
T
∫
Γ
w
1
,
y
δ
w
1
n
y
d
s
+
e
0
a
2
N
T
∫
Γ
w
1
,
x
x
x
δ
w
1
n
x
-
w
1
,
x
x
δ
w
1
,
x
n
x
+
w
1
,
x
y
y
δ
w
1
n
x
+
w
1
,
x
y
δ
w
1
,
x
n
y
d
s
+
e
0
a
2
N
T
∫
Γ
w
1
,
x
y
y
δ
w
1
n
x
+
w
1
,
x
y
δ
w
1
,
x
n
y
-
w
1
,
y
y
y
δ
w
1
n
y
+
w
1
,
y
y
δ
w
1
,
y
n
y
d
s
+
k
g
e
0
a
2
∫
Γ
w
1
,
x
x
x
δ
w
1
n
x
-
w
1
,
x
x
δ
w
1
,
x
n
x
-
w
1
,
y
y
y
δ
w
1
n
y
+
w
1
,
y
y
δ
w
1
,
y
n
y
+
2
w
1
,
x
y
y
δ
w
1
n
x
+
2
w
1
,
x
y
δ
w
1
,
x
n
y
d
s
+
D
∫
Γ
w
1
,
x
x
x
δ
w
1
n
x
-
w
1
,
x
x
δ
w
1
,
x
n
x
-
w
1
,
y
y
y
δ
w
1
n
y
+
w
1
,
y
y
δ
w
1
,
y
n
y
+
2
w
1
,
x
y
y
δ
w
1
n
x
+
2
w
1
,
x
y
δ
w
1
,
x
n
y
d
s
+
c
e
0
a
2
∫
Γ
w
1
,
y
δ
w
1
n
y
-
w
1
,
x
δ
w
1
n
x
+
w
2
,
x
δ
w
1
n
x
-
w
2
,
y
δ
w
1
n
y
d
s
+
k
g
∫
Γ
w
1
,
y
δ
w
1
n
y
-
w
1
,
x
δ
w
1
n
x
d
s
-
e
0
a
2
I
0
∫
Γ
w
1
,
x
t
t
δ
w
1
n
x
-
w
1
,
y
t
t
δ
w
1
n
y
d
s
·
I
2
∫
Γ
w
1
,
y
t
t
δ
w
1
n
y
-
w
1
,
x
t
t
δ
w
1
n
y
d
s
+
e
0
a
2
I
2
∫
Γ
w
1
,
x
x
x
t
t
δ
w
1
n
x
-
w
1
,
x
x
t
t
δ
w
1
,
x
n
x
-
w
1
,
y
y
y
t
t
δ
w
1
n
y
+
w
1
,
y
y
t
t
δ
w
1
,
y
n
y
+
2
w
1
,
x
y
y
t
t
δ
w
1
n
x
+
2
w
1
,
x
y
t
t
δ
w
1
,
x
n
y
d
s
=
0
.Equation (19) has the following weak form:(22)∫Ωδw2D∇4w2+1-e0a2∇2-q21+kww2-NT∇2w2-kg∇2w2+I0w¨2-I2w¨2,xxtt+w¨2,yyttdxdy=0Equation (22) can be dealt with similarly.


### 4.2. Element-Free Discretization Technique

UsingM values of the scattered particles (x1,x2,…,xM) to express the displacement field of SLGSs based on the kernel particle Ritz method [42] results in(23)w1w2=∑I=1MψIx,yw1Iw2I.where m is the number of nodes. wI and ψI, respectively, denote the displacement value and the shape function regarding node I. We write the shape function as(24)ψI=Cx;x-xIΦx-xI.Herein, Cx;x-xI is the coefficient function and Φx;x-xI represents the kernel function. Cx;x-xI can be expressed as(25)Cx;x-xI=HTx-xIBx,HTx-xI=1,x-xI,y-yI,x-xIy-yI,x-xI2,y-yI2T,Bx=k0x,y,k1x,y,k2x,y,k3x,y,k4x,y,k5x,yT,where H denotes the vector of the quadratic basis and B represents the function of x. Thus, we have(26)ΨIx;x-xI=BTxHx-xIΦx-xI,satisfying the shape function to reproduction conditions resulting in(27)∑I=1MΨIxxIpyIqforp+q=0,1,2.Coefficient B(x) has the following form:(28)Bx=M-1xH0,where(29)Mx=∑I=1MHx-xIΦx-xI,H0=1,0,0,0,0,0T.We have the kernel function(30)Φx-xI=Φx;xI·Φy;yI,where(31)Φx;xI=ϕx-xIaI.The cubic spline function is used in the present study(32)ϕs1=23-4s12+4s13for0≤s1≤1243-4s1+4s12-43s13for12≤s1≤10otherwise,where s1=x-xI/aI,aI=a(xI) are the dilation parameters. aI can be expressed as aI=dmaxcI, where dmax denotes the scaling factor ranging between 2 and 4. To ensure that the support area contains numerous enough nodes to prevent the matrix M from singularity, the factor should be selected appropriately. cI denotes the larger distances.


### 4.3. Discrete System Equations

WhereNT=-T(αhE)/[(1-υ)(1+φ)] and substituting (21) into (16), the discrete equation can be expressed as(33)M11w¨1+k11+kTw1+k12w2=0,where(34)M11=M11l+M11nl,(35)M11l=∫ΩI0ψIψJ+I2N1ITN1JdΩ,(36)M11nl=∫Ωe0a2I0N1ITN1J+N2ITIN2JdΩ,(37)k11=k11l+k11nl+kmedium,(38)k11l=∫ΩN2ITLDN2J-cψIψJdΩ,(39)k11nl=∫Ω-ce0a2N1ITN1JdΩ,(40)kmedium=∫ΩN4ITLGN4J+N3ITLwN3JdΩ,(41)k12=k12l+k12nl,(42)k12l=∫ΩcψIψJdΩ,(43)k12nl=∫Ωce0a2N1ITN1JdΩ,(44)kT=∫ΩN4ITLTN4JdΩ.Similarly, we can obtain(45)M22w¨2+k21w1+k22+kTw2=0(46)M11=M22,(47)k21=k21l+k21nl,(48)k21l=∫ΩcψIψJdΩ,(49)k21nl=∫Ωce0a2N1ITN1JdΩ,(50)k22=k22l+k22nl+kmedium,(51)k22l=∫ΩN2ITLDN2J-cψIψJdΩ,(52)k22nl=∫Ω-ce0a2N1ITN1JdΩ,(53)LT=-αhET1-υ1+φ100000e0a22e0a21e0a2,(54)I=I20002I2000I2,LD=D0002D000D,Lw=kw000e0a2kw000e0a2kw,(55)LG=kg00000kge0a2000002kge0a200000kg00000kge0a2,(56)N1I=ψI,xψI,y,N2I=ψI,xxψI,xyψI,yy,N3I=ψIψI,xψI,y,N4I=ψI,xψI,xxψI,xyψI,yψI,yy.We rearrange matrix equations (31) and (43) as(57)M1100M22w¨1w¨2+k11+kTk12k21k22+kT=0.We assume that(58)w1w2=w-1w-2eiωt.We rewrite (55) as(59)-ω2M1100M22+k11+kTk12k21k22+kTw-1w-2=0.By solving eigenvalue problem of (57), the buckling properties of DLCGSs can be presented.


## 4.1. Weak Form

In order to get the weak form, the equilibrium equation is multiplied by the virtual displacementδw. Therefore, we obtain(20)∫Ωδw1D∇4w1+1-e0a2∇2-q12+kww1-NT∇2w1-kg∇2w1+I0w¨1-I2w¨1,xxtt+w¨1,yyttdxdy=0.Equation (20) can be integrated by parts as(21)
D
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
-
c
∫
Ω
w
1
δ
w
1
-
w
2
δ
w
1
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
k
w
∫
Ω
w
1
δ
w
1
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
e
0
a
2
k
w
∫
Ω
w
1
,
x
δ
w
1
,
x
+
w
1
,
y
δ
w
1
,
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
N
T
∫
Ω
w
1
,
x
δ
w
1
,
x
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
-
c
e
0
a
2
∫
Ω
w
1
,
x
δ
w
1
,
x
+
w
1
,
y
δ
w
1
,
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
e
0
a
2
N
T
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
x
y
δ
w
1
,
x
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
c
e
0
a
2
∫
Ω
w
2
,
x
δ
w
1
,
x
+
w
2
,
y
δ
w
1
,
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
k
g
∫
Ω
w
1
,
x
δ
w
1
,
x
+
w
1
,
y
δ
w
1
,
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
e
0
a
2
∫
Ω
I
0
δ
w
1
,
x
w
1
,
x
t
t
+
δ
w
1
,
y
w
1
,
y
t
t
+
I
2
δ
w
1
,
x
x
w
1
,
x
x
t
t
+
2
δ
w
1
,
x
y
w
1
,
x
y
t
t
+
δ
w
1
,
y
y
w
1
,
y
y
t
t
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
e
0
a
2
N
T
∫
Ω
w
1
,
x
y
δ
w
1
,
x
y
+
w
1
,
y
y
δ
w
1
,
y
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
N
T
∫
Ω
w
1
,
y
δ
w
1
,
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
k
g
e
0
a
2
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
+
∫
Ω
I
0
δ
w
1
w
1
,
t
t
+
I
2
δ
w
1
,
x
w
1
,
x
t
t
+
δ
w
1
,
y
w
1
,
y
t
t
∫
Ω
w
1
,
x
x
δ
w
1
,
x
x
+
w
1
,
y
y
δ
w
1
,
y
y
+
2
w
1
,
x
y
δ
w
1
,
x
y
d
x
d
y
-
e
0
a
2
k
w
∫
Γ
w
1
,
x
δ
w
1
n
x
-
w
1
,
y
δ
w
1
n
y
d
s
-
N
T
∫
Γ
w
1
,
x
δ
w
1
n
x
d
s
+
N
T
∫
Γ
w
1
,
y
δ
w
1
n
y
d
s
+
e
0
a
2
N
T
∫
Γ
w
1
,
x
x
x
δ
w
1
n
x
-
w
1
,
x
x
δ
w
1
,
x
n
x
+
w
1
,
x
y
y
δ
w
1
n
x
+
w
1
,
x
y
δ
w
1
,
x
n
y
d
s
+
e
0
a
2
N
T
∫
Γ
w
1
,
x
y
y
δ
w
1
n
x
+
w
1
,
x
y
δ
w
1
,
x
n
y
-
w
1
,
y
y
y
δ
w
1
n
y
+
w
1
,
y
y
δ
w
1
,
y
n
y
d
s
+
k
g
e
0
a
2
∫
Γ
w
1
,
x
x
x
δ
w
1
n
x
-
w
1
,
x
x
δ
w
1
,
x
n
x
-
w
1
,
y
y
y
δ
w
1
n
y
+
w
1
,
y
y
δ
w
1
,
y
n
y
+
2
w
1
,
x
y
y
δ
w
1
n
x
+
2
w
1
,
x
y
δ
w
1
,
x
n
y
d
s
+
D
∫
Γ
w
1
,
x
x
x
δ
w
1
n
x
-
w
1
,
x
x
δ
w
1
,
x
n
x
-
w
1
,
y
y
y
δ
w
1
n
y
+
w
1
,
y
y
δ
w
1
,
y
n
y
+
2
w
1
,
x
y
y
δ
w
1
n
x
+
2
w
1
,
x
y
δ
w
1
,
x
n
y
d
s
+
c
e
0
a
2
∫
Γ
w
1
,
y
δ
w
1
n
y
-
w
1
,
x
δ
w
1
n
x
+
w
2
,
x
δ
w
1
n
x
-
w
2
,
y
δ
w
1
n
y
d
s
+
k
g
∫
Γ
w
1
,
y
δ
w
1
n
y
-
w
1
,
x
δ
w
1
n
x
d
s
-
e
0
a
2
I
0
∫
Γ
w
1
,
x
t
t
δ
w
1
n
x
-
w
1
,
y
t
t
δ
w
1
n
y
d
s
·
I
2
∫
Γ
w
1
,
y
t
t
δ
w
1
n
y
-
w
1
,
x
t
t
δ
w
1
n
y
d
s
+
e
0
a
2
I
2
∫
Γ
w
1
,
x
x
x
t
t
δ
w
1
n
x
-
w
1
,
x
x
t
t
δ
w
1
,
x
n
x
-
w
1
,
y
y
y
t
t
δ
w
1
n
y
+
w
1
,
y
y
t
t
δ
w
1
,
y
n
y
+
2
w
1
,
x
y
y
t
t
δ
w
1
n
x
+
2
w
1
,
x
y
t
t
δ
w
1
,
x
n
y
d
s
=
0
.Equation (19) has the following weak form:(22)∫Ωδw2D∇4w2+1-e0a2∇2-q21+kww2-NT∇2w2-kg∇2w2+I0w¨2-I2w¨2,xxtt+w¨2,yyttdxdy=0Equation (22) can be dealt with similarly.


## 4.2. Element-Free Discretization Technique

UsingM values of the scattered particles (x1,x2,…,xM) to express the displacement field of SLGSs based on the kernel particle Ritz method [42] results in(23)w1w2=∑I=1MψIx,yw1Iw2I.where m is the number of nodes. wI and ψI, respectively, denote the displacement value and the shape function regarding node I. We write the shape function as(24)ψI=Cx;x-xIΦx-xI.Herein, Cx;x-xI is the coefficient function and Φx;x-xI represents the kernel function. Cx;x-xI can be expressed as(25)Cx;x-xI=HTx-xIBx,HTx-xI=1,x-xI,y-yI,x-xIy-yI,x-xI2,y-yI2T,Bx=k0x,y,k1x,y,k2x,y,k3x,y,k4x,y,k5x,yT,where H denotes the vector of the quadratic basis and B represents the function of x. Thus, we have(26)ΨIx;x-xI=BTxHx-xIΦx-xI,satisfying the shape function to reproduction conditions resulting in(27)∑I=1MΨIxxIpyIqforp+q=0,1,2.Coefficient B(x) has the following form:(28)Bx=M-1xH0,where(29)Mx=∑I=1MHx-xIΦx-xI,H0=1,0,0,0,0,0T.We have the kernel function(30)Φx-xI=Φx;xI·Φy;yI,where(31)Φx;xI=ϕx-xIaI.The cubic spline function is used in the present study(32)ϕs1=23-4s12+4s13for0≤s1≤1243-4s1+4s12-43s13for12≤s1≤10otherwise,where s1=x-xI/aI,aI=a(xI) are the dilation parameters. aI can be expressed as aI=dmaxcI, where dmax denotes the scaling factor ranging between 2 and 4. To ensure that the support area contains numerous enough nodes to prevent the matrix M from singularity, the factor should be selected appropriately. cI denotes the larger distances.


## 4.3. Discrete System Equations

WhereNT=-T(αhE)/[(1-υ)(1+φ)] and substituting (21) into (16), the discrete equation can be expressed as(33)M11w¨1+k11+kTw1+k12w2=0,where(34)M11=M11l+M11nl,(35)M11l=∫ΩI0ψIψJ+I2N1ITN1JdΩ,(36)M11nl=∫Ωe0a2I0N1ITN1J+N2ITIN2JdΩ,(37)k11=k11l+k11nl+kmedium,(38)k11l=∫ΩN2ITLDN2J-cψIψJdΩ,(39)k11nl=∫Ω-ce0a2N1ITN1JdΩ,(40)kmedium=∫ΩN4ITLGN4J+N3ITLwN3JdΩ,(41)k12=k12l+k12nl,(42)k12l=∫ΩcψIψJdΩ,(43)k12nl=∫Ωce0a2N1ITN1JdΩ,(44)kT=∫ΩN4ITLTN4JdΩ.Similarly, we can obtain(45)M22w¨2+k21w1+k22+kTw2=0(46)M11=M22,(47)k21=k21l+k21nl,(48)k21l=∫ΩcψIψJdΩ,(49)k21nl=∫Ωce0a2N1ITN1JdΩ,(50)k22=k22l+k22nl+kmedium,(51)k22l=∫ΩN2ITLDN2J-cψIψJdΩ,(52)k22nl=∫Ω-ce0a2N1ITN1JdΩ,(53)LT=-αhET1-υ1+φ100000e0a22e0a21e0a2,(54)I=I20002I2000I2,LD=D0002D000D,Lw=kw000e0a2kw000e0a2kw,(55)LG=kg00000kge0a2000002kge0a200000kg00000kge0a2,(56)N1I=ψI,xψI,y,N2I=ψI,xxψI,xyψI,yy,N3I=ψIψI,xψI,y,N4I=ψI,xψI,xxψI,xyψI,yψI,yy.We rearrange matrix equations (31) and (43) as(57)M1100M22w¨1w¨2+k11+kTk12k21k22+kT=0.We assume that(58)w1w2=w-1w-2eiωt.We rewrite (55) as(59)-ω2M1100M22+k11+kTk12k21k22+kTw-1w-2=0.By solving eigenvalue problem of (57), the buckling properties of DLCGSs can be presented.


## 5. Numerical Simulations, Results, and Discussion

In this paper, the nonlocal elastic theory and element-free method are developed for the thermal vibration behavior of DLGSs. The mechanical properties of GSs are shown in Table1. To verify the effectiveness of the present study, the comparison of the results is given for circular nature frequency for DLGs with the side length of 10 nm for all edges simply supported by the result obtained by [43]. ω(m,n)i in this paper represents the ith interlayer vibration mode frequency. From Table 2, the results of the present study agree well with those of He et al. [43].Table 1
Mechanical properties.
E h ρ υ 1.06 TPa 0.34 nm 2250 Kg/m3 0.25Table 2
Comparison of the present study with the published ones for square DLGSs with widthb = 10 nm.
ω (THz) ω ( 1,1 ) 1 ω(1,2)1 ω(2,1)1 ω(2,2)1 ω(1,1)2 ω(1,2)2 ω(2,2)2 He et al. 0.069 0.173 0.173 0.276 2.683 2.688 2.697 Present 0.0692 0.175 0.175 0.2787 2.680 2.684 2.693Firstly, the influence of the side length and the mode (n,m) on the circular natural frequency of the square DLGSs under SSSS boundary condition is plotted in Figure 3. It can be studied that, for all of the first vibration mode, the circular frequencies decrease with the increasing side length whereas, for the second vibration mode, the circular frequencies are not affected by the side length.Figure 3
Frequency response with the side length for square DLGSs (T=0,kg=kw=0,SSSS,e0a=0).Secondly, we conduct a case study on the square DLGSs with 10 nm width to demonstrate the influence of nonlocal parameter on the circular natural frequency. From Figure4, we also can observe that, for all of the first vibration mode, the circular frequencies decrease with the increasing nonlocal parameter. However, for the second vibration mode, the circular natural frequencies are not influenced by the nonlocal parameters.Figure 4
Variation of frequency with the nonlocal parameter for square DLGSs (T=0,kg=kw=0,SSSS,Lb=10nm).Then, in Figure5, the influence of the aspect ratio on the frequency of DLGSs with the length of 10 nm under various boundary conditions is considered. It can be seen in Figure 5 that, under various boundary conditions, all the first vibration mode circular natural frequencies decrease as the aspect ratio increases. We also can study that, under the same condition, the first interlayer circular natural frequency under CCCC boundary condition is higher than the other boundary condition. However, the second interlayer circular natural frequencies under various boundary conditions are not affected by the aspect ratio and the frequencies are the same.Figure 5
Frequency response with the aspect ratio for DLGSs (T=0,kg=kw=0,Lb=10nm,e0a=0).The influence of stacking types of DLGSs on the vibration is also studied. Figure6 shows the frequency change with the side length under various boundary conditions and two stacking types of DLGSs. Figure 6 indicates that the first interlayer mode circular natural frequency is not affected by the stacking types of DLGSs. However, for second interlayer mode, the circular natural frequency for the AA stacking types is smaller than the AB stacking types.Figure 6
Frequency response with the stacking types for square DLGSs (T=0,kg=kw=0,e0a=0).To study the effects of Pasternak elastic foundation on the interlayer mode circular natural frequency of DLGSs, a case is simulated on square DLGSs with side length of 10 nm under various boundary conditions. The nondimensional Winker modulus parameterk-W and shear modulus parameter k-G are adopted. The nondimensional parameters are defined as (60)k-W=kwLx4D,k-G=kgLx2D.As shown in Figure 7, the Winker modulus parameter k-W=0. Figure 7 indicates that, for the first vibration mode, the circular frequencies increase with the increasing of k-G. The first natural frequency under CCCC boundary condition is larger than that under other boundary conditions. In contrast, the second natural frequency is not affected by the shear modulus parameter. Figure 8 shows the relationship between the interlayer circular natural frequency and Winker modulus parameter. It is evident that, for the first mode, the natural frequencies increase with the increasing of Winker modulus parameter. However, for the second mode, the circular natural frequency is not influenced by the Winker modulus parameter.Figure 7
Variation of frequency with the shear modulus parameter for square DLGSs (T=0,kw=0,Lb=10nm,e0a=0).Figure 8
Variation of frequency with the Winker modulus parameter for square DLGSs (T=0,kg=0,Lb=10nm,e0a=0).Finally, to comprehensively understand the effect of the temperature and temperature distribution coefficient on the interlayer mode circular natural frequency, we simulate the vibration behavior of square DLGSs with 10 nm width under SSSS boundary condition. Figure9 presents the frequency responses under low temperature environment; as shown in Figure 9, under low temperature environment, the interlayer mode circular natural frequency is not affected by the temperature and temperature distribution types. However, as is evident from Figure 10, under high temperature environment, the first interlayer mode circular natural frequency decreases with increasing of the temperature under various temperature distribution types. However, for the second mode, the natural frequency is not affected by either the temperature or the temperature distribution types.Figure 9
Variation of frequencies with the low temperature for square DLGSs (kw=kg=0,Lb=10nm,e0a=0,SSSS).Figure 10
Frequency change with the high temperature for square DLGSs (kw=kg=0,Lb=10nm,e0a=0,SSSS).


## 6. Conclusion

In this paper, the nonlocal elasticity theory and element-free method are utilized to reveal the vibration properties of DLGDs resting on a Winker’s elastic medium in thermal environment. Influence of side length, nonlocal parameters, boundary conditions, elastic medium parameter, stacking types of graphene, aspect ratio, and temperature distribution types on the interlayer mode circular natural frequency is examined. The results reveal that, for the first mode, the natural frequency is related to the nonlocal parameter, side length, boundary condition, elastic medium parameter, and aspect ratio, but is unrelated to the vdW force. On the contrary, for the second mode, the natural frequency is related to the vdW force, but is not affected by the side length, boundary condition, nonlocal parameter, aspect ratio, elastic medium parameter, and temperature distribution types.



---

*Source: 6508061-2018-04-11.xml*