Document Type: Research Paper
Authors
^{1} Department of Mechanical Engineering, Semnan University, Semnan, Iran
^{2} School of new Technology, Iran University of Science and Technology, Tehran, Iran
Abstract
Keywords

Mechanics of Advanced Composite Structures 6 (2019) 131  138


Semnan University 
Mechanics of Advanced Composite Structures journal homepage: http://MACS.journals.semnan.ac.ir 
M. Mehrabani ^{a}, M.R. Ashory ^{a,}^{*}, M.M. Khatibi ^{a}, S. Sadeghzadeh ^{b}
^{a }Department of Mechanical Engineering, Semnan University, Semnan, 3513119111, Iran
^{b }School of New Technologies, Iran University of Science and Technology (IUST), Tehran, 1684613114, Iran
Paper INFO 

ABSTRACT 
Paper history: Received 20181116 Received in revised form 20190418 Accepted 20190426 
Graphene sheets are combined of Honeycombs lattice carboncarbon bonds which have high natural frequencies, high strength, and high conductivity. Due to important applications of the graphene sheets particularly at higher frequencies, the study of their dynamic behavior is important in this frequency range. From Molecular Dynamics (MD) point of view as the dimensions of graphene sheet incline, the number of atoms increases, and as a result, its modeling becomes more timeconsuming. Besides the experimental methods in small dimensions are difficult to conduct and not economical. In this research Finite Element Method (FEM) is used for frequency analysis of graphene sheets in various dimensions in order to study the capability of FEM in simulating the dynamic behavior of graphene sheets at small scales. In this research, the objective function is to find the minimum size of the sheet in which both methods have good convergence. Also, the timeconsuming for the simulation is investigated. The timeconsuming for analysis in the Finite Element Method is less than other methods, including Molecular Dynamics (MD), Generalized Differential Quadrature (GDQ), etc. Also, The results indicated that for circular singlelayer graphene sheets simulation, using Finite Element Method (FEM) is in good agreement with the results obtained from the Molecular Dynamics (MD) simulation, in the radius more than 100 nm. In this research, the ABAQUS has been used for Finite Element Method (FEM) simulation. 



Keywords: Graphene Finite Element Method Frequency Analysis 


© 2019 Published by Semnan University Press. All rights reserved. 
Graphene, as one of the main allotropes of carbon in the singlelayer state, is made of honeycombs lattice carboncarbon bonds. This lattice has a very high natural frequency and also high strength and conductivity and can have higher properties in an electric field. The recent applications of graphene have been in different fields such as large scale measurement devices, sensors, transparent electrodes, solar cells, energy storage devices, polymer composites, and nanocomposites. Regarding the widespread applications of graphene at higher frequencies, the study of its different features particularly its dynamic behavior at this frequency range is important [1]. Recently, several studies have been conducted to clarify the dynamic behavior of graphene. In some of these investigations, the analytical methods have been used to study the vibrations of graphene sheets while in some others Molecular Dynamics simulations have been applied to derive the dynamic behavior of graphene. Leissa and Narita [2] studied the effect of Poisson's coefficient on the natural frequencies of a simply supported circular plate. Yongqiang and Jian [3] considered the effect of the different boundary conditions and different degrees of freedom on the natural frequencies of graphene. Aghababaei and Reddy [4] studied the vibrational properties of a single layer graphene sheet using thirdorder nonlocal shear deformation theory. Arash and Wang [5] used both nonlocal continuum mechanics theory and molecular dynamics simulation to study the free vibrations of single layer and bilayer graphene sheet and compared the results. Murmu and Pradhan [6] performed nonlocal elasticity theory to study the vibrational properties of singlelayer graphene sheets and determined the effect of elastic circumstances on the first mode frequency of singlelayer graphene sheets. NeekAmal and Peeters [7] investigated the effect of radial loads on buckling and stiffness properties of circular graphene singlelayer sheet using Molecular Dynamics model. Mahmoudinezhad and Ansari [8] modeled the vibrational properties of circular and square singlelayer graphene sheet using Finite Element Method (FEM). Asemi and Farajpour [9] studied the effects of thermal changes on vibrational behavior of circular graphene sheets. Gong et al. [10] modeled the circular graphene sheets and studied the vibrations of mass sensors. They also compared the natural frequencies of these sheets with the results of the Rayleigh method and investigated the effects of temperature and stress on the natural frequencies. Natsuki et al. [11] obtained the natural frequencies of bilayer graphene sheets. Mohammadi et al. [12] modeled the circular and annular graphene sheets with various boundary conditions by deriving their analytical equations. Mortazavi et al. [13] investigated the thermal conductivity of graphene epoxy nanocomposites using a hybrid method of FEM and molecular dynamics. Ansari et al. [14] used a nonlocal model in small dimensions to determine the characteristics of Multilayer graphene sheet vibrational properties in the elastic area for different boundary conditions. In this simulation, FEM was used to investigate the effect of length and modulus of elasticity on the vibrational behavior of graphene laminated sheets. Rouhi and Ansari [15] designed an atomic model of singlelayer graphene sheet using FEM in order to study its vibrational behavior. Based on their model, the natural frequencies of the graphene sheet was obtained with different boundary conditions and various dimensions. This study indicated that the vibrational behavior of graphene sheet could be modeled by using Molecular Dynamics simulation, numerical solution of the governing equations or analytical methods. The implementation of these methods requires the various stages, such as: extracting governing equations, a computational method and analyzing the vibrational signals. Molecular Dynamics (MD) method is commonly used to simulate the vibrational behavior of nanoscale sheets. Based on the literature [12], the accuracy of MD method is verified. When the size of the sheet is enlarged (for example, a circular sheet with a radius of more than 10 nanometer), the CPU time is usually prolonged, due to the increase in the number of particles. Based on the other researches [16], it is possible to simulate nanosheets in nanoscale by a Finite Element (FE) method (without considering nonlocal coefficients in the FE simulation). Although, the main problem is to find the minimum dimension of nanosheet, which the results of Finite Element method would have good agreement with Molecular Dynamics method. The main objective of this paper is to find the minimum size of the circular sheet in which both methods have good convergence. Furthermore, the solving time is presented in both finite element and molecular dynamics method. For this purpose, At the First, the theory of this problem is explained, then the effect of mesh resizing on the convergence of the solutions for the natural frequency values is checked and the best mesh size is suggested. Finally, the simulation results are verified. The result of FE modeling in various dimensions with those of the Molecular Dynamics modeling and the results obtained by other researchers would be compared. Furthermore, the simulation CPU time is studied.
2.1. Free vibrations of a singlelayer circular sheet based on the nonlocal theory
The graphene singlelayer sheet displacement equation based on the nonlocal method is given by Eq. (1) [12]:
(1) 
In this equation, is the transverse displacement, is the nonlocal coefficient, and are the Winkler modulus, and is the laplacian operator. It is worth to mention that if , Eq. (1) transforms into a classic one. Also, is the flexural strength calculated according to Eq. (2) [12];
(2) 
where E is the modulus of elasticity, h is the sheet thickness and is the Poisson coefficient. The displacement of the sheet can be calculated from Eq. (3) [12]:
(3) 
where and indicates the natural frequency.
By inserting Eq. (2) into Eq. (1) and using Eq. (3), the following would be derived [12]:
(4) 
In which and are obtained as follows [12],
(5) 

(6) 
By solving the eigenvalue problem given by Eq. (4), the values of the dimensionless natural frequencies are obtained as Eq. (7) [12]:
(7) 
2.2. Finite Element Model
FEM is one of the numerical solution methods of Partial Differential Equations (PDEs) and also integral equations. In this method, the numerical solution of equations is given by converting the partial differential equations into Ordinary Differential Equations (ODEs). The aim of solving partial differential equations is to achieve a simple and stable equation so that it does not lead to inaccurate and unreasonable results. This method solves the problem by dividing a continuous domain into subdomains called elements. Researchers have invented various types of element especially triangular and quadrilateral elements, some of them are very complex. According to the Kirchhoff law, shear deformation is neglected in thin plates. Using equilibrium relations and assuming that the flexural strength, D, and the force, are constant, the differential equation of the fourth order is obtained as [17]:
(8) 
If the mechanical behavior of the material is homogeneous, the flexural stiffness can be obtained from Eq. (2) [17].
In this research, square elements with 4 nodes are used. This element has 12 degrees of freedom in each node; therefore in order to approximate the freemotion of , a polynomial with 12 parameters are implemented (Eq. (9)) [17].
(9) 
whereα coefficients are constant and are obtained considering the boundary conditions. Also, considering the deflection, in each of four nodes of the element, the shape function is determined as follows (10):
(10) 
In this Equation, a and b are constant coefficients and depend on the type of element [17]. Also, and are the local coordinates of each element node, and and are local coordinate axes and are defined as Eqs. (11) and (12).
(11) 

(12) 
Considering the above equations, the stiffness matrix of each element can be calculated as Eq. (13) [17]:
(13) 
B Matrix is defined as follows:
(14) 
In this equation, L is the Lagrangian and N is a shape function of the element [17].
At first, according to the mechanical properties of graphene, the sheet is modeled by ABAQUS. In the present work, the vibrational properties of the circular graphene sheet are examined in different dimensions. This result is compared with those of other studies. Another simulation is implemented using molecular dynamics (MD) model. Taking into account intermolecular forces, MD model has a high degree of accuracy in calculating the behavior of nanoscale sheets. In this research, MD and FE methods are compared. As a result, the best simulation method is represented, based on accuracy and CPU time. In order to model in terms of Finite Element (FE), square element with 4 nodes has been selected in this research. In this type of element, each node has 3 degrees of freedom. Therefore, the square element has 12 degrees of freedom. It is worth to mention that the strain matrix is not fixed in linear square elements; therefore this type of element, gives more accurate results for the stress and strain matrices. Also, the governing equations are simpler due to their regular and symmetrical geometry. Therefore, the CPU time of analysis decreases.
3.1. Mechanical Properties of the Circular Sheet
In order to calculate the natural frequencies and mode shapes of the circular graphene sheet, ABAQUS is implemented. The properties of this sheet are given in Table 1. In this table, refers to mass density, E is elasticity modulus, h is the thickness of graphene sheet and is Poisson’s ratio.
3.2. Mesh Convergence
The circular sheet is simulated using ABAQUS, according to the properties listed in Table 1. In order to achieve the desired accuracy, the optimum number of elements is selected. After choosing the best type of element, the simulation is done. The meshed graphene sheet is shown in Fig. 1. This figure is a sample. The optimum number of elements should be selected.
Table 1. Mechanical properties of the singlelayer graphene sheet [12].
h(nm) 
E(TPa) 

0.33 
0.34 
1.06 
2300 
Fig. 1. Graphene Elemental Model. 
In order to ensure that the results of the Finite Element (FE) Model are accurate, the convergence analysis is implemented, and the effect of change in number of elements on natural frequencies in the simulated model was checked. Obviously, as the number of elements increases, the size of the elements becomes smaller. This study has been carried out for various sheet size. For example; the variations diagram of the third and sixth natural frequencies of the circular sheet with a radius of 20 nm have been presented in Figs. 2 and 3. The most suitable mesh size of 8400 elements has been suggested for the computation of natural frequencies by taking into account the lowest error value.
Fig. 2. Study the Effect of number of elements on the third natural frequency. 
Fig. 3. Study the Effect of number of elements on the sixth natural frequency. 
4.1. Results of Finite Element Model
In this section, after achieving desire accuracy, the results of a circular graphene sheet analysis, with the properties given in Table 1, are presented. The natural frequencies and mode shapes of this sheet are shown.
For validation, the comparison of the results of this work with those of another study is indicated in Table 2. The value of the nonlocal coefficient is set according to the desired frequencies. This comparison is executed for sheets with the radius of 6, 7, 8, 10, 15 and 20 nm.
Table 2. Comparison of natural frequencies (GHz) of the circular sheet with simply support boundary conditions.
error (%) 
Finite Element (present study) 
Mohammadi et al. [12] 
Number of modes 
radius (nm) 

Frequency 
a 

0.292 
49.0878 
49.232 
0 
1 
6 
0.433 
136.828 
136.2383 
0.5 
2 

0.433 
136.828 
136.2383 
0.5 
3 

0.821 
249.741 
247.7065 
0.5 
4 

1.764 
251.037 
255.5455 
0 
5 

1.635 
291.667 
296.5175 
0 
6 

0.165 
36.1107 
36.1705 
0 
1 
7 
0.768 
100.863 
100.0935 
0.5 
2 

0.768 
100.863 
100.0935 
0.5 
3 

1.393 
184.524 
181.9884 
0.5 
4 

1.114 
185.655 
187.7477 
0 
5 

1.037 
215.590 
217.8496 
0 
6 

0.33 
27.7846 
27.693 
0 
1 
8 
0.092 
78.0700 
77.9978 
0 
2 

0.395 
78.3064 
77.9978 
0 
3 

0.211 
143.440 
143.7443 
0 
4 

1.608 
146.056 
143.7443 
0 
5 

1.485 
169.268 
166.7911 
0 
6 

0.041 
17.7308 
17.7235 
0 
1 
10 
0.398 
49.7199 
49.9186 
0 
2 

0.279 
49.7792 
49.9186 
0 
3 

0.723 
91.3311 
91.9964 
0 
4 

0.206 
92.1862 
91.9964 
0 
5 

0.262 
107.027 
106.7463 
0 
6 

0.856 
7.85775 
7.791 
0 
1 
15 
0.049 
21.9545 
21.9436 
0 
2 

0.049 
21.9545 
21.9436 
0 
3 

0.286 
40.3246 
40.4404 
0 
4 

0.286 
40.3246 
40.4404 
0 
5 

0.168 
46.8451 
46.9243 
0 
6 

0.894 
4.4217 
4.3825 
0 
1 
20 
0.182 
12.3658 
12.3433 
0 
2 

0.183 
12.3660 
12.3433 
0 
3 

0.062 
22.7334 
22.7477 
0 
4 

0.021 
22.7447 
22.7477 
0 
5 

0.050 
26.4083 
26.3949 
0 
6 
Table 3. Comparison of natural frequencies (GHz) of the simply supported sheet with radius 20 (nm) with other references.
Finite Element (present study) 
Reference [18] 
Reference [2] 
Reference [12] 
Number of modes 
4.4217 
4.38308 
4.38308 
4.3825 
1 
12.3658 
12.3432 
12.3432 
12.3433 
2 
12.3660 
12.3432 
12.3432 
12.3433 
3 
22.7334 
22.7493 
22.7476 
22.7477 
4 
22.7427 
22.7493 
22.7476 
22.7477 
5 
26.4083 
26.3942 
26.3951 
26.3949 
6 
c mode 3 
b mode 2 
a mode 1 
f mode 6 
e mode 5 
d mode 4 
Fig. 4. Mode shape of a circular sheet with the radius of 20 nm from FEM model: a) Mode 1, b) Mode 2, c) Mode3, d) Mode 4, e) Mode 5 and f) Mode 6.
In order to ensure the accuracy of the results, the natural frequencies of the circular sheet with the radius of 20 nm are validated against several previous works which are summarized in Table 3.
In addition to the natural frequencies, the mode shapes of the singlelayer graphene sheet are also obtained, which are shown in Fig. 4, for a sheet with the radius of 20 nm.
4.2. Results of Molecular Dynamics Model
The aim of this section is to find the lowest radius of the circular singlelayer graphene sheet, which corresponds to the natural frequency calculated by FEM compared to the natural frequency obtained from the MD method. Molecular Dynamics (MD) has been suggested in studies that are costly and timeconsuming. This method has often been applied in very small dimensions by researchers. According to the physical laws, simulation using Molecular Dynamics (MD) has been conducted by determining the type, position, properties of atoms and the force between them. In the following, the Molecular Dynamics (MD) model of the circular graphene sheet is also carried out, and its results are compared with those of FEM. A circular graphene sheet that is modeled by Molecular Dynamics method is shown in Fig. 5.
For simulating the graphene sheet in the Molecular Dynamics (MD) model, CarbonCarbon bonds are simulated by Tersoff potential model. The simulation is executed in room temperature. As the initial condition, a pulse is applied to graphene sheet by steptime of 0.001 ps. In this method, the particle displacement is considered as an output [19]. In the frequency decomposition method that is presented by Brincker et al. [20], the natural frequencies and mode shapes are extracted after estimation of spectral density matrix and applying the decomposition method for the singular values on it. For this purpose, after ABAQUS simulation, the first natural frequency of the sheet with different radius is compared with MD Model and the results are summarized in Table 4. Also, the CPU time of simulation using both methods is presented in this table. The configuration of the system that is utilized for the simulation is Intel^{®} Xeon^{®} Processor X5675 (12M Cache, 3.06 GHz, 6.40 GT/s Intel® QPI) FCLGA10) with Random Access Memory as 64 GB.
Fig. 5. A circular Graphene sheet with the radius of 8 nm that is modeled by Molecular Dynamic Method in LAMMPS software. 
Table 4. Comparison of the first mode frequency (GHz) of the sheet with simply supported boundary condition in two Finite Element (FE) Models and Molecular Dynamics (MD) Models.
ref. [14] 
MD 
FEM 
Radius (nm) 

Timeconsuming (min) 
Frequency (GHz) 
Timeconsuming (min) 
Frequency (GHz) 

27.693 
372 
39.0625 
0.11 
27.7846 
8 
4.3825 
1254 
13.4277 
0.23 
4.4263 
20 
1.8234 
4856 
8.5449 
0.43 
1.9674 
30 
0.9576 
6249 
7.3242 
0.64 
1.1062 
40 
0.4365 
7203 
4.8828 
0.82 
0.49161 
60 
As a result of the comparison, the first natural frequency of FE and MD models versus radius of the graphene sheet is plotted in Fig. 6. According to this figure, the predicted radius that results of FE and MD models are matched is 100 nm. In low radius, the nonlocal effects are high. The Finite Element Model does not consider these effects. Therefore, the difference between MD and FE models in low radius is high. This difference in the radius greater than 100 nm is reduced to a desirable level. Also, according to table 4, the CPU time for MD simulation is more than FE simulation. Simulation time for graphene sheets with radius of higher than 20 nm, take long more than one day, although this time for FE simulation is less than one minute. According to table 4 in comparison with table 2, the nonlocal coefficient by increasing the radius of the sheet tends to zero. In other words, by increasing the sheet size, the effects of intermolecular forces decrease as well as the difference between natural frequencies of MD and FE simulation decreases.
According to Table 2 and 3, the second and third modes, as well as the fourth and fifth modes, are doubled modes. As shown in Fig. 6, whatever the radius of the sheet increased, the FE and MD results became closer together. It’s due to a reduction in the nonlocal effects in Molecular Dynamics (MD) model, Also the model is approached to the classic model. Therefore, due to the timeconsuming simulation by Molecular Dynamics (MD) model, results larger than 60 nm radius is predicted by extrapolation. Consequently, in large dimensions, graphene singlelayer circular sheet, due to a large number of atoms, simulation by Molecular Dynamics (MD) is very time consuming, and ABAQUS can be used for this modeling.
In this article, the dynamic behavior of circular graphene sheets has been investigated using FE and MD methods. In order to guarantee the accuracy of the results of the FE method, the natural frequencies of the nanosheet are verified by those of other studies [12]. Also, the effect of sheet size on the accuracy of simulation by Finite Element is studied. This research aimed to find the minimum size of nanosheet, which the results of finite element method would have good agreement with molecular dynamics method. The main objective of this paper is to find the minimum size of the circular sheet in which both methods have good convergence. Furthermore, the solving time is presented in both finite element and molecular dynamics method. Generally, the results of this study can be summarized as follows:
References
[1] Singh V, Joung D, Zhai L, Das S, Khondaker SI, Seal S. Graphene based materials: past, present and future. Progress in materials science 2011; 56(8): 1178271.
[2] Leissa AW, Narita Y. Natural frequencies of simply supported circular plates. Journal of Sound and Vibration 1980; 70(2): 2219.
[3] Yongqiang L, Jian L. Free vibration analysis of circular and annular sectorial thin plates using curve strip Fourier pelement. Journal of Sound and Vibration 2007; 305(3): 45766.
[4] Aghababaei R, Reddy JN. Nonlocal thirdorder shear deformation plate theory with application to bending and vibration of plates. Journal of Sound and Vibration 2009; 326(1): 27789.
[5] Arash B, Wang Q. Vibration of Single and DoubleLayered Graphene Sheets. Journal of Nanotechnology in Engineering and Medicine 2011; 2(1): 0110127.
1.
(a) 
(b) 
(c) 
(d) 
(e) 
(f) 
Fig. 6. Compare the frequencies of the two methods and predict their matching: (a) Mode 1, (b) Mode 2, (c) Mode 3, (d) Mode 4, (e) Mode 5 and (f) Mode 6. 
[6] Murmu T, Pradhan SC. Vibration analysis of nanosinglelayered graphene sheets embedded in elastic medium based on nonlocal elasticity theory. Journal of Applied Physics 2009; 105(6): 064319.
[7] NeekAmal M, Peeters FM. Buckled circular monolayer graphene: a graphene nanobowl. Journal of Physics: Condensed Matter 2011; 23(4): 045002.
[8] Mahmoudinezhad E, Ansari R. Vibration analysis of circular and square singlelayered graphene sheets: An accurate spring mass model. Physica E: Lowdimensional Systems and Nanostructures 2013; 47: 126.
[9] Asemi SR, Farajpour A. Decoupling the nonlocal elasticity equations for thermomechanical vibration of circular graphene sheets including surface effects. Physica E: Lowdimensional Systems and Nanostructures 2014; 60: 8090.
[10] Gong X, Jiang S, Wang X, Liu S, Wang S. Vibration analysis of nanomechanical mass sensor based on circular graphene sheets2014.
[11] Natsuki T, Shi JX, Ni QQ. Vibration analysis of circular doublelayered graphene sheets. Journal of Applied Physics 2012; 111(4): 044310.
[12] Mohammadi M, Ghayour M, Farajpour A. Free transverse vibration analysis of circular and annular graphene sheets with various boundary conditions using the nonlocal continuum plate model. Composites Part B: Engineering 2013; 45(1): 3242.
[13] Mortazavi B, Benzerara O, Meyer H, Bardon J, Ahzi S. Combined molecular dynamicsfinite element multiscale modeling of thermal conduction in graphene epoxy nanocomposites. Carbon 2013; 60: 35665.
[14] Ansari R, Rajabiehfard R, Arash B. Nonlocal finite element model for vibrations of embedded multilayered graphene sheets. Computational Materials Science 2010; 49(4): 8318.
[15] Rouhi S, Ansari R. Atomistic finite element model for axial buckling and vibration analysis of singlelayered graphene sheets. Physica E: Lowdimensional Systems and Nanostructures 2012; 44(4): 76472.
[16] Tserpes K, Papanikos P. Finite element modeling of singlewalled carbon nanotubes. Composites Part B: Engineering 2005; 36(5): 46877.
[17] Mirzaei M, Rouzegar J. Crack simulation in shells and plates using Extended Finite Element Method. Modares Civil Engineering 1389; 10(4): 112.
[18] Zhou ZH, Wong KW, Xu XS, Leung AYT. Natural vibration of circular and annular thin plates by Hamiltonian approach. Journal of Sound and Vibration 2011; 330(5): 100517.
[19] Rapaport DC. The Art of Molecular Dynamics Simulation: Cambridge University Press; 1995.
[20] Rune B, Lingmi Z, Palle A. Modal identification of outputonly systems using frequency domain decomposition. Smart Materials and Structures 2001; 10(3): 441.