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

1School of Electronic and Information Engineering, Suzhou University of Science and Technology, Suzhou 215009, China 2School of Science, Nanjing University of Science and Technology, Nanjing 210094, China 3Institute for Automatic Control and Complex Systems (AKS), Faculty of Engineering, University of Duisburg-Essen, Duisburg, Germany 4College of Electrical Engineering, Zhejiang University, Hangzhou 310027, China


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][4][5][6][7][8][9][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 scalefree 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][22][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][36][37][38][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.

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 in which  0  is the nonlocal parameter and ∇ 2 is the Laplacian operator defined by

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.Pasternaktype foundation model is modeled with two parameters.The mathematical expressions of Pasternak-type foundation are written as where   is the Winker modulus and   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]: 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  =   −   , where   is the temperature value when the buckling occurs, while   is the referent temperature.The thermal load caused by temperature change satisfies the following formulation: AB stacking AA stacking

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] in which   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 3546 Gpa/nm (for AA-type) ,   = −146.8247Gpa/nm (for AB-type) . ( If there is no special description, the vdW coefficient is taken as AB-type in the following simulations.

Governing Equations
in which (, V, ), respectively, denote the displacements along the axes (, , ). 1 is the displacement along axis of a mid-plane material point.According to the straindisplacement relations the constitutive relationship is given by where  = /2(1 + ).
According to (10), the nonlocal moments can be written as where ( 1 ,  1 ,   ) = ∫( 1 ,  1 ,  1 ) . is the bending stiffness of the plate and  = ℎ 3 With the principle of virtual work, the equilibrium equation can be derived as From ( 11)-( 13), we can get Substituting ( 17) into ( 16), the governing equation can be written as Similarly, we can obtain the governing equation of bottom-GSs

Discretized System Equations
4.1.Weak Form.In order to get the weak form, the equilibrium equation is multiplied by the virtual displacement .Therefore, we obtain Equation ( 20) can be integrated by parts as Equation ( 19) has the following weak form: Equation ( 22) can be dealt with similarly.

Element-Free Discretization Technique.
Using  values of the scattered particles (x 1 , x 2 , . . ., x  ) to express the displacement field of SLGSs based on the kernel particle Ritz method [42] results in where  is the number of nodes.  and   , respectively, denote the displacement value and the shape function regarding node .We write the shape function as where H denotes the vector of the quadratic basis and B represents the function of x.Thus, we have satisfying the shape function to reproduction conditions resulting in Coefficient B(x) has the following form: where We have the kernel function where The cubic spline function is used in the present study where  1 = |−  |/  ,   = (  ) are the dilation parameters.  can be expressed as   =  max   , where  max 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.  denotes the larger distances.

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 Table 1.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].  (,) in this paper represents the th interlayer vibration mode frequency.From Table 2, the results of the present study agree well with those of He et al. [43].
Firstly, the influence of the side length and the mode (, ) 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.
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 Figure 4, 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.
Then, in Figure 5, 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 Side length (nm) (1,1)     natural frequencies under various boundary conditions are not affected by the aspect ratio and the frequencies are the same.
The influence of stacking types of DLGSs on the vibration is also studied.Figure 6 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.
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 parameter   and shear modulus parameter   are adopted.The nondimensional parameters are defined as As shown in Figure 7, the Winker modulus parameter   = 0. Figure 7 indicates that, for the first vibration mode, the circular frequencies increase with the increasing of   .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.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.Figure 9 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.

Conclusion
In this paper, the nonlocal elasticity theory and elementfree 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.

Table 2 :
Comparison of the present study with the published ones for square DLGSs with width  = 10 nm.