Comparing quantum, molecular and continuum models for graphene at large deformations

TL;DR

This study compares the accuracy of molecular and continuum models for graphene under large deformations, validating them against DFT data and analyzing their predictive capabilities for mechanical and vibrational properties.

Contribution

It evaluates three interatomic potentials and a continuum shell model, providing insights into their validity and limitations for modeling graphene's large deformation behavior.

Findings

Continuum model aligns well with DFT results.

MM3 potential is accurate up to material instability.

Tersoff potential predicts auxetic behavior.

Abstract

In this paper, the validity and accuracy of three interatomic potentials and the continuum shell model of Ghaffari and Sauer [1] are investigated. The mechanical behavior of single-layered graphene sheets (SLGSs) under uniaxial stretching, biaxial stretching and pure bending is studied for this comparison. The validity of the molecular and continuum models is assessed by direct comparison with density functional theory (DFT) data available in the literature. The molecular simulations are carried out employing the MM3, Tersoff and REBO+LJ potentials. The continuum formulation uses an anisotropic hyperelastic material model in the framework of the geometrically exact Kirchhoff-Love shell theory and isogeometric finite elements. Results from the continuum model are in good agreement with those from DFT. The results from the MM3 potential agree well up to the point of material instability, whereas those from the REBO+LJ and Tersoff potentials agree only for small deformations. Only the Tersoff potential is found to yield auxetic response in SLGSs under uniaxial stretch. Additionally, the transverse vibration frequencies of a pre-stretched graphene sheet and a carbon nanocone are obtained using the continuum model and molecular simulations with the MM3 potential. The variations of the frequencies from these approaches agree within an error of 5%.

Figures

Tables

Table 1. Table 1 : Elastic properties of graphene reported in the literature. NA = Not available.

Ref.	Year	Method/potential	$E_{\mathrm{AC}}$ [TPa]	$E_{\mathrm{ZZ}}$ [TPa]	Thickness [nm]
[17]	2000	Ab-initio	1.100	1.100	0.340
[18]	2006	Molecular structural mechanics	1.096–1.125	1.106–1.201	0.340
		(Tersoff-Brenner)
[19]	2006	DFT	1.250	1.250	0.340
[20]	2007	Ab-initio	1.050	1.050	0.334
[3]	2008	Experimental (nano indentation)	1 $\pm$ 0.100	1 $\pm$ 0.100	0.335
[21]	2009	DFT	0.964	0.964	0.340
[22]	2009	Molecular structural mechanics	1.040	1.042	0.340
[23]	2009	Quantum molecular dynamics	1.100	0.600	0.335
[24]	2009	Orthogonal tight-binding and MD	1.01 $\pm$ 0.030	1.01 $\pm$ 0.030	0.335
[25]	2009	Truss-type analytical models
		with AMBER	1.378	1.303	NA
		with MORSE	1.379	1.957	NA
[26]	2010	MD (AIREBO)	0.890	0.830	NA
[27]	2010	Molecular structural mechanics	0.721	0.737	0.340
[28]	2010	MD (Tersoff)	1.130	1.050	0.335
[29]	2010	Molecular statics (MM3)	3.380	3.400	0.100
[30]	2012	MD (AIREBO)	1.097	0.961	0.335
[31]	2013	Molecular structural mechanics and MD	1.070	1.070	0.335
[32]	2014	Molecular structural mechanics	1.100	1.100	0.34
		(Modified MORSE)
[33]	2017	Space frame approach
		with AMBER	0.780	0.819	NA
		with MORSE	0.890	0.938	NA
[34]	2018	Multiscale model (MM3)	0.927	0.927	0.335

Table 2. Table 2 : Membrane behavior: Hyperelastic material constants determined by fitting to DFT calculations based on generalized gradient approximation (GGA) according to Kumar and Parks [ 70 ] and Shirazian et al. [ 55 ] .

	$\hat{\alpha }$	$\epsilon \left[\text{N/m}\right]$	${\mu }_0\left[\text{N/m}\right]$	${\mu }_1\left[\text{N/m}\right]$	$\hat{\beta }$	${\eta }_0\left[\text{N/m}\right]$	${\eta }_1\left[\text{N/m}\right]$
Kumar and Parks [70]	1.53	93.84	172.18	27.03	5.16	94.65	4393.26
Shirazian et al. [55]	1.435	103.9	182.6	34.94	4.665	83.46	3932

Table 3. Table 3 : Bending stiffness, c b subscript 𝑐 b c_{\text{b}} in [nN ⋅ ⋅ \cdot nm] (in [eV]), according to various atomistic models. FGBP = first generation Brenner potential; SGBP = second generation Brenner potential; QM = quantum mechanics.

FGBP Lu et al. [71]	SGBP Lu et al. [71]	QM Kudin et al. [72]	Shirazian et al. [55]
0.133 (0.83)	0.225 (1.40)	0.238 (1.49)	0.238 (1.49)

Table 4. Table 4 : The mean and standard deviation of the bond lengths for the reference configuration of the SLGS.

Potential	Mean [nm]	Standard deviation [nm]
MM3	0.1345	0.00020
REBO+LJ	0.1395	0.00036
Tersoff	0.1460	0.00028
modified Tersoff	0.1474	0.00026

Table 5. Table 5 : Bending stiffness ( c b subscript 𝑐 b c_{\text{b}} ) obtained from an equivalent plate model. Here K 𝐾 K , ν 𝜈 \nu and h ℎ h are the basal plane stiffness, Poisson ratio and the effective thickness of a graphene sheet, respectively. Here, first three rows of the table are calculated from Eq.  27 and remaining are from Eq.  28 .

Potential/model	$K$ [N/m]	$\nu$	$h$ [nm]	$c_{\text{b}}$ [eV]
MM3 [29]	340.0	0.210	0.100	1.86
REBO-I [80]	235.0	0.412	0.062	0.56
REBO-II [80]	243.0	0.397	0.057	0.49
MM3 (Current work)	-	-	-	2.11
REBO+LJ (Current work)	-	-	-	2.27
Tersoff (Current work)	-	-	-	2.08
modified Tersoff (Current work)	-	-	-	0.72

Table 6. Table 6 : Bending stiffness c b subscript 𝑐 b c_{\text{b}} in [eV] from an equivalent CNT model.

Potential/model	$c_{\text{b}}$ (armchair)	$c_{\text{b}}$ (zigzag)
MM3 (Current work)	3.271	3.146
REBO+LJ (Current work)	2.184	2.235
Tersoff (Current work)	2.078	2.010
modified Tersoff (Current work)	0.915	0.969

Table 7. Table 7 : The potential parameters of MM3 [ 35 ] .

Parameter	value
$K_{\text{s}}$	4.49 mdyne/Å
$K_{\theta }$	0.67 mdyne-Å/rad2
$V_1$	0.185 Kcal/mol
$V_2$	0.170 Kcal/mol
$V_3$	0.520 Kcal/mol
${ϵ}_{\text{e}}$	0.027 Kcal/mol
$r_{\text{v}}$	2.04 Å
$K_{\text{s}\theta }$	0.130 mdyne/rad
$K_{\varphi s}$	0.059 mdyne/rad
$K_{\theta {\theta }^{{}^{\prime }}}$	0.240 mdyne-Å/rad2
$r^{\text{e}}$	1.5247 Å
${\theta }^{\text{e}}$	110.2°

Table 8. Table 8 : The REBO+LJ potential parameters for carbon [ 39 ] .

Parameter	value
$Q$ (Å)	0.313460
$\alpha$ (Å-1)	4.746539
$A$ (eV)	10953.544
$B^1$ (eV)	12388.792
$B^2$ (eV)	17.567065
$B^3$ (eV)	30.714932
${\beta }^1$ (Å-1)	4.720452
${\beta }^2$ (Å-1)	1.433213
${\beta }^3$ (Å-1)	1.382691
$ϵ$ (eV)	0.002840
$\sigma$ (Å)	3.4
$r^{\text{min}}$ (Å)	1.7
$r^{\text{max}}$ (Å)	2.0

Table 9. Table 9 : The Tersoff and modified Tersoff potential parameters for carbon.

Parameter	Tersoff [36]	Erhart and Albe [56]
$A$ (eV)	$1393.600$	2019.800
$B$ (eV)	$346.700$	175.400
$\lambda$ (Å-1)	3.487	4.184
$\mu$ (Å-1)	2.2119	1.931
$\beta$	$1.5724\times {10}^{-7}$	1.000
$n$	0.727	1.000
$c$	38049.00	181.910
$d$	$4.384$	6.284
$h$	-0.57058	-0.5556
$r^{\text{min}}$ (Å)	1.80	1.85
$r^{\text{max}}$ (Å)	2.10	2.15
$\chi$	1.0	1.0
$\omega$	1.0	1.0

Equations

$$\begin{array}[]{lll}\boldsymbol{F}_{\!I}=m_{I}\,\ddot{\boldsymbol{r}}_{\!I}~{},\end{array}$$

$$\begin{array}[]{lll}\boldsymbol{F}_{I}=-\displaystyle\frac{\partial{U}}{\partial{\boldsymbol{r}_{I}}}~{},\end{array}$$

$$\begin{array}[]{lll}\bar{\sigma}_{I}^{ij}\!\!\!&:=&\!\!\!\displaystyle\frac{1}{V_{I}}\displaystyle\sum_{J=1}^{N_{p}}f_{IJ}^{i}\,r_{IJ}^{j}~{},\end{array}$$

$$\begin{array}[]{lll}\sigma_{I}^{ij}\!\!\!&:=&\!\!\!\displaystyle\displaystyle\frac{V_{I}}{A_{I}}\bar{\sigma}_{I}^{ij}~{},\end{array}$$

$$\begin{array}[]{lll}U_{\text{MM3}}=\displaystyle\sum_{I}\sum_{J}\big{(}U_{\text{s}}+U_{\theta}+U_{\phi}+U_{\text{s}\theta}+U_{\theta\theta^{{}^{\prime}}}+U_{\phi\text{s}}\big{)}+\sum_{I}\sum_{K}U_{\text{vdW}}~{},\end{array}$$

$$\begin{array}[]{lll}U_{\text{REBO}}=\displaystyle\sum_{I}\sum_{J=I+1}\big{[}E_{\text{R}}(r_{IJ})+b_{IJ}\,E_{\text{A}}(r_{IJ})\big{]}~{},\end{array}$$

$$\begin{array}[]{lll}U_{\text{LJ}}=\displaystyle 4\epsilon\left[\left(\frac{\sigma}{r_{IJ}}\right)^{12}-\left(\frac{\sigma}{r_{IJ}}\right)^{6}\right]~{},\end{array}$$

$$\begin{array}[]{lll}U_{\text{Tersoff}}=\displaystyle\frac{1}{2}\displaystyle\sum_{I}\sum_{J\neq I}f_{\text{c}}(r_{IJ})\big{[}f_{\text{R}}(r_{IJ})+b_{IJ}\,f_{\text{A}}(r_{IJ})\big{]}~{},\end{array}$$

$$\begin{array}[]{l}W(\mathcal{J}_{1},\mathcal{J}_{2},\mathcal{J}_{3})=W_{\text{m}}^{\mathrm{dil}}(\mathcal{J}_{1})+W_{\text{m}}^{\mathrm{dev}}(\mathcal{J}_{2},\,\mathcal{J}_{3};\mathcal{J}_{1})+W_{\text{b}}(\kappa_{1},\kappa_{2};\mathcal{J}_{1})~{},\end{array}$$

$$\begin{array}[]{lll}W^{\mathrm{dil}}_{\mathrm{m}}\!\!\!&:=&\!\!\!\varepsilon\big{[}1-(1+\hat{\alpha}\,\mathcal{J}_{1})\,\exp(-\hat{\alpha}\,\mathcal{J}_{1})\big{]}~{},\end{array}$$

$$\begin{array}[]{lll}W^{\mathrm{dev}}_{\mathrm{m}}\!\!\!&:=&\!\!\!2\,\mu(\mathcal{J}_{1})\,\mathcal{J}_{2}+\eta(\mathcal{J}_{1})\,\mathcal{J}_{3}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle W_{\mathrm{b}}\!\!\!&:=&\!\!\!\displaystyle J\frac{c_{\text{b}}}{2}\left(\kappa_{1}^{2}+\kappa_{2}^{2}\right)~{},\end{array}$$

$$\begin{array}[]{lll}\mu(\mathcal{J}_{1})\!\!\!&:=&\!\!\!\mu_{0}-\mu_{1}\,e^{\hat{\beta}\,\mathcal{J}_{1}}~{},\\[5.69054pt] \eta(\mathcal{J}_{1})\!\!\!&:=&\!\!\!\eta_{0}-\eta_{1}\,\mathcal{J}_{1}^{2}~{}.\end{array}$$

$$\begin{array}[]{lll}\mathcal{J}_{1}\,\!\!\!&:=&\!\!\!\,\ln J=\ln(\lambda_{1}\lambda_{2})~{},\\[5.69054pt] \displaystyle\mathcal{J}_{2}\!\!\!&:=&\!\!\!\displaystyle\,\frac{1}{4}\left(\frac{\lambda^{2}_{1}}{\lambda^{2}_{2}}+\frac{\lambda^{2}_{2}}{\lambda^{2}_{1}}-2\right)~{},\\[5.69054pt] \displaystyle\mathcal{J}_{3}\!\!\!&:=&\!\!\!\displaystyle\,\frac{1}{8}\left(\frac{\lambda_{1}}{\lambda_{2}}-\frac{\lambda_{2}}{\lambda_{1}}\right)^{3}\,\cos(6\theta)~{}.\end{array}$$

$$\begin{array}[]{lll}\theta\,\!\!\!&:=&\!\!\!\,\arccos{(\boldsymbol{Y}_{\!1}\cdot\hat{\boldsymbol{x}})}~{},\end{array}$$

$$\begin{array}[]{lll}\mathbf{M}\,\ddot{\mathbf{u}}+\mathbf{f}_{\text{int}}(\mathbf{u})=\mathbf{f}_{\text{ext}}(\mathbf{u})~{},\end{array}$$

$$\begin{array}[]{lll}\mathbf{f}_{\text{int}}^{e}\!\!\!&:=&\!\!\!\displaystyle\int\limits_{\Omega_{0}^{e}}{\tau^{\alpha\beta}\,\mathbf{N}_{,\alpha}^{\text{T}}\,\boldsymbol{a}_{\beta}~{}\mathrm{d}A}+\int\limits_{\Omega_{0}^{e}}{M_{0}^{\alpha\beta}\,{\mathbf{N}}_{;\alpha\beta}^{\text{T}}\,\boldsymbol{n}~{}\mathrm{d}A}~{},\end{array}$$

$$\begin{array}[]{lll}\mathbf{f}_{\text{ext}}^{e}\!\!\!&=&\!\!\!\displaystyle\int\limits_{\partial_{t}\Omega^{e}}{\mathbf{N}^{\text{T}}\,\boldsymbol{t}~{}\mathrm{d}s}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle\tau^{\alpha\beta}\,\!\!\!&=&\!\!\!\displaystyle\,\frac{\partial{W}}{\partial{a_{\alpha\beta}}}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle M^{\alpha\beta}_{0}\,\!\!\!&=&\!\!\!\displaystyle\,\frac{\partial{W}}{\partial{b_{\alpha\beta}}}~{},\end{array}$$

$$\begin{array}[]{lll}\sigma^{\alpha\beta}=\displaystyle\frac{1}{J}\tau^{\alpha\beta}~{}.\end{array}$$

$$\begin{array}[]{lll}\mathbf{M}\,\Delta\ddot{\mathbf{u}}+\mathbf{K}\,\Delta\mathbf{u}=-(\mathbf{f}_{\text{int}}+\mathbf{M}\,\ddot{\widehat{\mathbf{u}}}-\mathbf{f}_{\text{ext}})~{},\end{array}$$

$$\begin{array}[]{lll}\Delta\mathbf{u}\,\!\!\!&:=&\!\!\!\,\Delta\bar{\mathbf{u}}\,e^{-i\omega t}~{},\end{array}$$

$$\begin{array}[]{lll}\mathbf{K}\,\Delta\bar{\mathbf{u}}=\omega^{2}\,\mathbf{M}\,\Delta\bar{\mathbf{u}}~{}.\end{array}$$

$$\begin{array}[]{lll}\left[\begin{matrix}x\\ y\\ \end{matrix}\right]=\left[\begin{matrix}\lambda_{1}&0\\ 0&\lambda_{2}\\ \end{matrix}\right]\left[\begin{matrix}X\\ Y\\ \end{matrix}\right],\end{array}$$

$$\begin{array}[]{lll}E(X)=\displaystyle\frac{X-X_{\text{DFT}}}{\max(X_{\text{DFT}})}~{},\end{array}$$

$$\begin{array}[]{lll}c_{\text{b}}=\displaystyle\frac{Eh^{3}}{12(1-\nu^{2})}~{},\end{array}$$

$$\begin{array}[]{lll}f_{mn}=\displaystyle\frac{\lambda_{mn}^{2}}{2\pi a^{2}}\sqrt{\frac{c_{\mathrm{b}}}{\rho_{\mathrm{s}}}}~{},\end{array}$$

$$\begin{array}[]{lll}U_{\text{s}}\!\!\!&=&\!\!\!\displaystyle 71.94K_{\text{s}}(r_{IJ}-{r^{\text{e}}})^{2}\left[1-2.55(r_{IJ}-r^{\text{e}})+\left(\frac{7}{12}\right)2.55(r_{IJ}-r^{\text{e}})^{2}\right]~{},\\[8.53581pt] U_{\theta}\!\!\!&=&\!\!\!\displaystyle 0.021914K_{\theta}(\theta_{IJK}-{\theta^{\text{e}}})^{2}\biggl{[}1-0.014(\theta_{IJK}-{\theta^{\text{e}}})+5.6(10)^{-5}(\theta_{IJK}-{\theta^{\text{e}}})^{2}\\[5.69054pt] \!\!\!&-&\!\!\!\displaystyle 7.0(10)^{-7}(\theta_{IJK}-{\theta^{\text{e}}})^{3}+9.0(10)^{-10}(\theta_{IJK}-{\theta^{\text{e}}})^{4}\biggr{]}~{},\\[8.53581pt] U_{\phi}\!\!\!&=&\!\!\!\displaystyle\frac{V_{1}}{2}(1+\cos\phi_{IJKL})+\frac{V_{2}}{2}(1-\cos 2\phi_{IJKL})+\frac{V_{3}}{2}(1+\cos 3\phi_{IJKL})~{},\\[8.53581pt] U_{\text{s}\theta}\!\!\!&=&\!\!\!\displaystyle 2.51118K_{\text{s}\theta}\left[(r_{IJ}-{r^{\text{e}}})+(r_{JK}-{r^{\text{e}}})\right](\theta_{IJK}-{\theta^{\text{e}}})~{},\\[8.53581pt] \end{array}$$

$$\begin{array}[]{lll}U_{\phi\text{s}}\!\!\!&=&\!\!\!\displaystyle 11.995\frac{K_{\phi\text{s}}}{2}(r_{IJ}-{r^{\text{e}}})(1+\cos(3\phi_{IJKL}))~{},\\[8.53581pt] U_{\theta\theta^{{}^{\prime}}}\!\!\!&=&\!\!\!\displaystyle-0.021914K_{\theta\theta^{{}^{\prime}}}(\theta_{IJK}-{\theta^{\text{e}}})(\theta_{IKL}-{\theta^{\text{e}}})~{},\\[8.53581pt] U_{\text{vdW}}\!\!\!&=&\!\!\!\displaystyle\epsilon_{\text{e}}\left[-2.25\left(\frac{r_{\text{v}}}{r_{IJ}}\right)^{6}+1.84(10)^{5}\text{exp}\left(-12.0\frac{r_{IJ}}{r_{\text{v}}}\right)\right]~{},\end{array}$$

Full text

\UseRawInputEncoding

**Comparing quantum, molecular and continuum models for graphene at large deformations

**

Aningi Mokhalingam*†111Email: aningi@iitk.ac.in, Reza Ghaffari§222Email: ghaffari@aices.rwth-aachen.de, Roger A. Sauer§333Email: sauer@aices.rwth-aachen.de and Shakti S. Gupta†*444Corresponding author, email: ssgupta@iitk.ac.in

§*Aachen Institute for Advanced Study in Computational Engineering Science (AICES),

RWTH Aachen University, Templergraben 55, 52056 Aachen, Germany

† Department of Mechanical Engineering, IIT Kanpur, Kanpur - 208016, India*

Published555This pdf is the personal version of an article whose journal version is available at https:/​/sciencedirect.com in Carbon, DOI: 10.1016/j.carbon.2019.12.014

Submitted on 2 October 2019; Revised on 6 December 2019; Accepted on 7 December 2019

Abstract

In this paper, the validity and accuracy of three interatomic potentials and the continuum shell model of Ghaffari and Sauer [1] are investigated. The mechanical behavior of single-layered graphene sheets (SLGSs) under uniaxial stretching, biaxial stretching and pure bending is studied for this comparison. The validity of the molecular and continuum models is assessed by direct comparison with density functional theory (DFT) data available in the literature. The molecular simulations are carried out employing the MM3, Tersoff and REBO+LJ potentials. The continuum formulation uses an anisotropic hyperelastic material model in the framework of the geometrically exact Kirchhoff-Love shell theory and isogeometric finite elements. Results from the continuum model are in good agreement with those from DFT. The results from the MM3 potential agree well up to the point of material instability, whereas those from the REBO+LJ and Tersoff potentials agree only for small deformations. Only the Tersoff potential is found to yield auxetic response in SLGSs under uniaxial stretch. Additionally, the transverse vibration frequencies of a pre-stretched graphene sheet and a carbon nanocone are obtained using the continuum model and molecular simulations with the MM3 potential. The variations of the frequencies from these approaches agree within an error of $\approx 5\%$.

Keywords: Anisotropic hyperelasticity; carbon nanocone; continuum mechanics; graphene; Kirchhoff-Love shells; molecular dynamics.

1 Introduction

Graphene is an atom-thick two-dimensional (2D) hexagonal lattice of covalently bonded carbon atoms, which can be exfoliated from bulk graphite [2]. Due to its excellent mechanical [3], electrical and thermal properties [4], it has many industrial applications in the fields of nanocomposites [5], nano-electromechanical systems [6] and electronic devices [7]. The electronic properties, in particular the band gap of a single-layered graphene sheet (SLGS) can be altered by applying uniaxial or biaxial strains [8, 9, 10] on it. This method is popularly called as strain engineering. Also, uniaxial or biaxial strain in an SLGS of a bilayered graphene sheets (BLGSs) leads to change in the stacking dislocations that commence changing electronic or optical properties [11]. However, successful operation of such strain-tunable devices will hinge on the satisfactory mechanical response of SLGS under applied strain field, small or large.

Another carbon based nanostructure obtained due to pentagonal and heptagonal defects in graphene is carbon nanocone (CNC) [12, 13]. CNCs have potential applications in the scanning probe microscopy [14], field emission electron source [15] and molecular pumps [16].

The mechanical response, static or dynamic, of a SLGS/CNC under imposed strains can be studied comprehensively using quantum-mechanical or molecular statics/dynamics (MS/MD) simulations. Alternatively, one can develop a suitable equivalent continuum model of the SLGS/CNC which reproduces identical response as obtained from these simulations for a large range of applied strain field. We note that the quantum-mechanical simulations of CNCs will be computationally inefficient, however, a continuum membrane model calibrated using quantum-mechanical or MS/MD simulations based data can be employed to study deformations/dynamics of CNCs.

In the following, we begin with surveying the literature reporting linear response of SLGSs and the material constants derived from it. In Table 1, the elastic moduli of a SLGS along the armchair, $E_{\text{AC}}$, and zigzag, $E_{\text{ZZ}}$, directions obtained from these methods are reported. As seen, some authors report different stiffness in the armchair and zigzag directions. This table also reveals ambiguity in the literature regarding the isotropic behavior of SLGS, even in the small deformation regime. Subsequently, in this paper, relevant discussion on this aspect from our findings are reported.

In MD and MS simulations, the accuracy of the results depends on the potential defining atomic interactions. For carbon structures, popular potentials are: MM3 [35], Tersoff [36], the first and second generation reactive empirical bond order (REBO) [37, 38], adaptive intermolecular reactive empirical bond order (AIREBO) [39], and ReaxFF [40]. Recently, Lebedeva et al. [41] investigated the applicability of these potentials (except MM3) up to 3$\%$ uniaxial elongation and reported that the considered potentials fail to reproduce precisely the experimental and ab-initio in-plane and out-of-plane deformations of a SLGS. Employing the MM3 potential in MS simulations, Gupta and Batra [29] reported significant increase of the transverse vibration frequencies of SLGS under pre-stretch. Similarly, Liao et al. [42] studied the influence of temperature, cone height, and cone angles on the mechanical behavior of CNCs under uniaxial tensile and compression employing the Tersoff potential.

Some equivalent continuum structures have been proposed in the linear regime to study vibrations of SLGSs [43, 44, 22, 27, 45, 46, 47] and CNCs [48, 49]. Using a continuum plate model, Jiang et al. [44] investigated the effect of size, shape and boundary conditions on vibration of SLGSs and multi-layered graphene sheets (MLGSs). Apart from the elastic properties, the vibrational behavior of SLGSs has also been studied using molecular structural mechanics [50, 51]. Sadeghi and Naghdabadi [51] studied the nonlinear vibrations of SLGS and reported that the fundamental frequency depends on the amplitude of vibration and length of SLGS. Using MS based on the universal force field (UFF) model, Chowdhury et al. [45] reported that the natural frequencies of SLGSs are insensitive to the chirality. Singh and Patel [52] used a multiscale method to study the effect of pre-tension on the nonlinear static and dynamic response of SLGSs. They have reported that the pre-tension significantly increases natural frequencies. Singh and Patel [53] also studied the nonlinear elastic response of SLGSs and reported that SLGSs show softening behavior at small strains and hardening behavior at large bending. Lu and Huang [54] studied the uniaxial in-plane and bending deformations of SLGS using molecular simulations employing REBO potential and equivalent continuum model. They have reported that SLGS exhibits nonlinear and anisotropic response under finite axial stretch.

There are some studies on vibration of nanocones which are in their relaxed state. For example, Fakhrabadi et al. [48] studied the vibrational properties of CNCs of different heights and cone angles using MS simulations and reported that the transverse vibrational frequencies reduce with increase in cone angle and height. Hu et al. [49] modeled CNCs as tapered beams using the Euler-Bernoulli and Timoshenko beam theory to study transverse vibrations.

In this paper, the nonlinear mechanical response of a square graphene sheet under uniaxial and biaxial in-plane stretch is studied using MS/MD, DFT and a hyperelastic continuum model based on DFT data [55]. The second generation REBO+LJ and Tersoff potentials are used in the MD simulations, and the MM3 potential is used in the MS simulations. For the Tersoff potential, two sets of parameterization are considered: the original parameters proposed by Tersoff [36] and the modified parameters proposed by Erhart and Albe [56]. These two parameterizations are denoted “Tersoff” and “modified Tersoff”, respectively, henceforth. The validity of the interatomic potentials is investigated in the linear and nonlinear deformation regime by comparing the results from MS/MD simulations with those obtained from DFT simulations [55]. The frequencies of the transverse vibrations of the SLGS and a CNC under stretch, obtained from MS simulations, are compared with those obtained from the continuum model.

The remainder of this paper is arranged as follows: Section 2 describes the molecular simulation methods and the interatomic potentials used in the present study. A brief introduction to the continuum model is given in Section 3. Numerical results are then presented in Section 4 followed by conclusions in Section 5.

2 Molecular simulations

This section presents details about molecular simulations and the mathematical expressions of the interatomic potentials considered.

Molecular simulations are carried out by solving the Newtonian equations of motion for the atoms,

[TABLE]

where $\boldsymbol{r}_{I}$ and $m_{I}$ are the position and mass of the atom $I$, and $\boldsymbol{F}_{I}$ is the interatomic force on atom $I$. $\boldsymbol{F}_{I}$ can be determined from the potential function $U\!\left(\text{{r}}^{\textit{N}}\right)$ as

[TABLE]

where $\mathbf{r}^{N}\mathrel{\mathop{:}}=\{\boldsymbol{r}_{1},\boldsymbol{r}_{2},\boldsymbol{r}_{3},...,\boldsymbol{r}_{N}\}$ are the positions of all the atoms.

The response of a graphene sheet under uniaxial and biaxial stretching, and bending can be computed using MS simulations for all the potentials discussed here. However, the computation of the bending stiffness based on plate vibration requires modal analysis. Whereas Tinker (coded with the MM3 potential) provides modal analysis, LAMMPS (coded with the REBO+LJ and Tersoff potentials) does not. Therefore, in Tinker the simulations are performed at 0 K, while in LAMMPS we perform MD at extremely low temperatures (0.1 K), that do not introduce entropic effects, and use fast Fourier transform to extract the vibration frequencies from the transverse time response of the atoms.

At 0 K, the average velocity fluctuations of the atoms are zero. These systems are thus static problems that can be solved by finding their minimum energy configurations. Here, this is done using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [57]. Before applying any loads, the system under consideration (here SLGS/CNC) should itself be brought into the minimum energy configuration called the reference configuration. We do this initial relaxation step by bringing the root mean square gradient of $U$ to 0.001 Kcal/mol/Å. Subsequently, we apply the desired boundary conditions and compute the eigenvalues (frequencies) and eigenvectors (mode shapes) of the mass-weighted Hessian of the SLGS/CNCs at various stretch values using the Vibrate subroutine of Tinker [58].

In MD simulations, the specification of a finite temperature leads to the fluctuating velocity in the system. MD systems thus require transient solution approaches. Further, they need to be thermally equilibrated in order to reach quasi-static states at macro-scales. However, first the relaxed or minimum energy configuration of the system should be obtained. This is done here with the Polak-Ribiere’s conjugate gradient method [59] in a quasi-static approach with absolute zero temperature and zero velocity. The energy666LAMMPS uses normalized energy units for the energy minimization criterion, defined as the relative change of the energy between two successive iterations divided by the energy magnitude at the current iteration. and force tolerances for terminating the energy minimization are set to $10^{-8}$ and $10^{-8}$ eV/Å,  respectively. Subsequently, the system is thermally equilibrated at constant volume and temperature of 0.1 K for 50 ps with a timestep of 1 fs to obtain the reference configuration. The temperature is maintained at 0.1 K employing a Nose-Hoover thermostat [60] with three chains during deformations. The MD simulations are performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [61].

In molecular simulations at approximately 0 K, the virial stress is defined by [62, 63]

[TABLE]

where $V_{I}$ is the volume occupied by atom $I$, and $i$ and $j$ denote the Cartesian components along the x, y and z directions, respectively. $N_{p}$ is the number of neighboring atoms of atom $I$, $f_{IJ}^{i}$ is the force on atom $I$ due to atom $J$, and $r_{IJ}$ is the distance between atom $I$ and $J$. For 2D materials, such as graphene, it is natural to introduce the stress as force per length. Given (3), this stress follows as

[TABLE]

where $A_{I}$ is the sheet area attributed to atom $I$. Through definition Eq. (4), the usage of a thickness for SLGSs is avoided.

The interatomic potentials are described in the following subsections.

2.1 MM3 Potential

The MM3 potential consists of the first and higher order expansions of bond stretching, angle bending, and torsion. This potential also incorporates the cross terms777The cross term between the bending and torsion is not considered. among the mentioned contributions and between angle bending and out-of-plane bending [35]. The expression of the MM3 potential is given by [35]

[TABLE]

where $U_{\text{s}}$, $U_{\theta}$ and $U_{\phi}$ are the energy contributions corresponding to changes in bond length, bond angle and dihedral angle, respectively. $U_{\text{s}\theta}$, $U_{\theta\theta^{{}^{\prime}}}$ and $U_{\phi\text{s}}$ account for energies of the cross-term interactions between bond stretch and angle bending, angle bending and out-of-plane-bending, and bond stretch and dihedral angle, respectively. $U_{\text{vdW}}$ defines van der Waals (vdW) attraction and steric repulsion in the form $({r_{\text{c}}}/r_{IJ})^{6}$ and $\exp(-12r_{IJ}/r_{\text{c}})$, where $r_{c}$ is the cut-off distance. Further details of these terms are given in Appendix A.1.

2.2 REBO+LJ Potential

The REBO+LJ potential consists of two parts. The covalent bonds between carbon atoms are modeled using the second generation REBO potential, which is widely used for the formation and breaking of bonds in carbon structures. The REBO part of the potential is [39]

[TABLE]

where $r_{IJ}$ is the distance between a pair of atoms and $b_{IJ}$ is an empirical bond-order term. $E_{\text{R}}$ and $E_{\text{A}}$ are, respectively, the repulsive and attractive terms and are given in Appendix A.2. The vdW attraction and steric repulsion are modeled by the standard 12-6 Lennard-Jones (LJ) potential [64]

[TABLE]

where $\sigma$ and $\epsilon$ are the LJ parameters. The vdW energy is only included when the covalent bond energy from the REBO potential become zero, i.e., after breakage of the covalent bond.

2.3 Tersoff Potential

The Tersoff potential is a pair-like potential in which the strength of a bond depends on the local environment, i.e., an atom with fewer neighboring atoms forms a stronger bond than the atom with more neighboring atoms. In the literature, the mechancial behavior of carbon structures has been studied successfully using the Tersoff [65, 66] and modified Tersoff [67] potentials. Both are described by [36, 56]

[TABLE]

where $r_{IJ}$ is the distance between the atom pairs $I$ and $J$, $b_{IJ}$ is the bond order term, $f_{\text{R}}$ and $f_{\text{A}}$ are the repulsive and attractive terms, respectively, and $f_{\text{c}}$ is a smooth cutoff function. Details of the potential are given in Appendix A.3.

3 Continuum model

As a homogenized structure, the SLGS and CNC are modeled based on the shell formulation of Duong et al. [68] and the anisotropic hyperelastic material model of Ghaffari and Sauer [1]. Ghaffari and Sauer [1] formulated the strain energy density, per unit area of the initial configuration, based on a set of invariants $\mathcal{J}_{i}$ , i.e. [69, 1]

[TABLE]

where $W_{\text{m}}^{\mathrm{dil}}$ and $W_{\text{m}}^{\mathrm{dev}}$ are the pure dilatational and deviatoric parts of the membrane strain energy density, respectively, and $W_{\text{b}}$ is the bending strain energy density. These terms are defined as

[TABLE]

[TABLE]

[TABLE]

where $\mu$ and $\eta$ are defined as

[TABLE]

The material constants $\varepsilon$, $\hat{\alpha}$, $\mu_{0}$, $\mu_{1}$, $\hat{\beta}$, $\eta_{0}$, $\eta_{1}$, and $c_{\text{b}}$ from Shirazian et al. [55] are used in the continuum model, see Tables 2 and 3.

$J$ is the surface area change, and $\kappa_{1}$ and $\kappa_{2}$ are the principal surface curvatures [73]. $\mathcal{J}_{1}$ and $\mathcal{J}_{2}$ capture isotropic dilatation and shear deformation, respectively, while $\mathcal{J}_{3}$ captures anisotropic shear deformation. $\mathcal{J}_{i}$ are given by

[TABLE]

$\lambda_{1}$ and $\lambda_{2}$ are the principal surface stretches with $\lambda_{1}\geq\lambda_{2}$. $\theta$ is the maximum stretch angle relative to the armchair direction and defined as

[TABLE]

where $\boldsymbol{Y}_{\!1}$ is the direction of the maximum stretch, see Ghaffari et al. [69] and Ghaffari and Sauer [1] for details. The material has isotropic behavior under pure dilatation, and anisotropic behavior only appears under large shear deformation. Material model (9) is implemented in the nonlinear finite shell element formulation of Duong et al. [68]. The discretized weak form can be written as [68]

[TABLE]

where $\mathbf{M}$ is the mass matrix (see Ghaffari and Sauer [74] for the mass matrix of graphene), $\mathbf{u}$ is the displacement vector and $\mathbf{f}_{\text{int}}$ and $\mathbf{f}_{\text{ext}}$ are the internal and external force vectors, respectively. They are assembled from the elemental contributions $\mathbf{f}_{\text{int}}^{e}$ and $\mathbf{f}_{\text{ext}}^{e}$, respectively, using the standard finite element assembly procedure. $\mathbf{f}_{\text{int}}^{e}$ and $\mathbf{f}_{\text{ext}}^{e}$ are defined by

[TABLE]

[TABLE]

for the special case of zero body forces and zero boundary moments. Here $\Omega_{0}^{e}$ and $\Omega^{e}$ denote element domains in the reference and current configuration, respectively, and along $\partial_{\boldsymbol{t}}\Omega$, the boundary traction $\boldsymbol{t}$ is applied. $\boldsymbol{a}_{\alpha}$ are the covariant tangent vectors and $\boldsymbol{n}$ is the normal vector to the shell surface. $\mathbf{N}$, $\mathbf{N}_{,\alpha}$ and $\mathbf{N}_{;\alpha\beta}$ denote the shape function matrix and its first parametric derivative and second covariant derivative [68], and $(\bullet)^{\text{T}}$ is the transpose operator. The contra-variant components of the surface Kirchhoff stress tensor and the moment tensor are given by [68]

[TABLE]

[TABLE]

where $a_{\alpha\beta}$ and $b_{\alpha\beta}$ are the surface metric and curvature tensor components. The Cauchy and Kirchhoff stress components are connected by

[TABLE]

Using a Taylor expansion about $\widehat{\mathbf{u}}$ such that $\mathbf{u}=\widehat{\mathbf{u}}+\Delta\mathbf{u}$, the linearized equations of motion become

[TABLE]

where $\mathbf{K}:=\partial(\mathbf{f}_{\text{int}}-\mathbf{f}_{\text{ext}})/\partial\mathbf{u}$ is the tangent stiffness matrix (see Ghaffari and Sauer [74] and Duong et al. [68] for details). For harmonic vibrations we have

[TABLE]

where $\Delta\bar{\mathbf{u}}$ are $\omega$ are the mode shape and frequency of the structure, respectively. Using (23), Eq. (22) can be transformed to the standard eigenvalue problem

[TABLE]

At each load increment, we solve Eq. (16) iteratively using the Newton-Raphson method. Subsequently eigenvalue problem (24) is solved at each load increment in order to obtain the variation of the frequencies with the loading.

4 Numerical results

This section presents numerical results pertaining to the deformation and vibrations of a square SLGS and a CNC. The three molecular models of Sec. 2 are compared with DFT data from Shirazian et al. [55] and the DFT-based continuum model of Sec. 3. First, the minimum energy configuration of the SLGS is given in Subsection 4.1. Then the procedure for the in-plane stretch of the SLGS is described in Subsection 4.2.1, and the variation of strain energy and stresses is studied in Subsections 4.2.2 and 4.2.3. The out-of-plane bending of SLGS is then discussed in Subsection 4.3, followed by the modal analysis of SLGS and CNC in Subsection 4.4. The transverse vibration frequencies of SLGS are calculated up to the instability point, where the ellipticity of the elasticity tensor is lost [70]. A pointwise summary of the main findings is finally given in Subsection 4.5.

4.1 Minimum energy configuration of SLGS

A SLGS with dimensions 10 nm $\times$ 10 nm and an initial bond length of 0.142 nm is considered. The mean and standard deviation of the bond length for the reference configuration employing the three potentials are shown in Table 4. In all cases, the equilibrium (i.e. mean) bond length is slightly different than the experimental value 0.1422 nm [75]. The difference can be attributed to the various functional forms of the potentials and the constants used therein. The deviations from the mean are due to the numerical errors that depend on the convergence tolerances considered for the energy minimization. Further, from Table 4, we note that after the minimization, the SLGS shrinks more when the MM3 potential is employed compared to REBO+LJ. On the other hand, the SLGS dilates after relaxation when the Tersoff potential is used.

4.2 In-plane stretching of SLGS

This section examines the elastic response of the square SLGS under uniaxial and biaxial stretch.

4.2.1 Procedure

The MS simulations (with MM3) of the SLGS are performed with and without periodic boundary conditions (PBCs) to investigate the size effect. In the absence of PBCs, the position in the current configuration ($x$, $y$) of the edge atoms is given by

[TABLE]

where $X$ and $Y$ are the initial position of the edge atoms in the reference configuration. For uniaxial stretch, $\lambda_{1}=\lambda$ and $\lambda_{2}=1$ or $\lambda_{1}=1$ and $\lambda_{2}=\lambda$ (i.e. the lateral direction is kept fixed), while for pure dilatation, $\lambda_{1}=\lambda_{2}=\lambda$ (see Fig. 1). When PBCs are used, the periodic simulation box is deformed along $X$ and/or $Y$ with a stretch increment of 0.1 Å. At each increment, the edge atoms are kept fixed and the system is relaxed to compute the potential energy and virial stress.

In the MD simulations (with REBO+LJ and Tersoff), PBCs are employed. The SLGS is stretched uniaxially or biaxially using the deformation control method. A side of the periodic box ($l$) at any instant of time ($t$) is deformed according to $l(t)=l_{0}(1+e_{rate}t)$, where $l_{0}=10\ \text{nm}$, is the initial length of the box and the applied strain rate is $e_{rate}=0.001$/ps888This strain rate translates to the velocity 10 m/s for a SLGS with 10 nm length, which is much much less than the speed of sound wave propagation in SLGS, $\approx$ 21,000 m/s, and hence eliminates any transient effect. [61]. These parameters values are optimum to study the elastic behavior of the SLGS according to Zhao et al. [24].

In the continuum simulations, (25) is applied quasi-statically without using PBCs.

4.2.2 Strain energy

Figure 2 shows the potential energy variation versus the uniaxial and biaxial stretch computed from molecular simulations, the DFT-based continuum model and DFT [55]. The relative error in the strain energy $W$ is shown in Fig. 2d, 2e and 2f. The relative error in the quantity $X$ with respect to the DFT results is defined as

[TABLE]

where $X$ results from the continuum and/or molecular simulations and $X_{\text{DFT}}$ is the DFT result. The maximum value is used to avoid division by zero. In MS simulations, it is found that the variation in the energy and stress (discussed subsequently) are almost identical with and without PBCs. These results confirm that the size effect on the elastic response is negligible when the diagonal length of the SLGS is over 10.0 nm as reported by Zhao et al. [24].

The strain energy results from the MM3 potential agree with the DFT results up to the stretch $\lambda_{1}\approx 1.15$ for uniaxial loading and $\mathcal{J}_{1}\approx 0.3$ for biaxial loading within $\approx 5\%$ error. The results from the REBO+LJ and Tersoff potentials agree within an error of $\approx$ 20$\%$. The continuum results are in excellent agreement with those from DFT simulations for the whole range under study (within $\approx 0.05\%$ error). This is due to the fact, that the continuum model has been calibrated directly from DFT data [55].

4.2.3 Stresses

Next, the stresses from the different approaches are compared. For uniaxial stretch, the SLGS is stretched either along the armchair or zigzag direction. Three different stresses are examined in the following.

1. *Stress along the stretch direction.

*The variation in $\sigma_{11}$ – the stress along the stretch direction – and the corresponding error are shown in Fig. 3 for stretch along the armchair direction, and Fig. 4, for stretch along the zigzag direction. Fig. 3 shows that the molecular simulation results are in good agreement with the DFT results up to $\lambda_{1}=1.1$ within an error of $\approx 5\%$.

Beyond this stretch, $\sigma_{11}$ computed from the MM3 potential shows gradual hardening. $\sigma_{11}$ computed from the other two potentials follow the DFT results up to $\lambda_{1}=1.19$ and then there is a sudden rise, which is discusssed subsequently.

Similarly Fig. 4 shows that the molecular simulation results follow the DFT results up to $\lambda_{1}=1.13$. Beyond this stretch, the MM3 results exhibit much higher stresses than the other cases. This may be due to the absence of cutoff function in MM3 as is present in the other two potentials, which initiates bond breaking in the structure. The results from the other two potentials continue to follow the DFT results up to $\lambda_{1}=1.25$.

As Fig. 3 shows a sudden rise in $\sigma_{11}$ is noticed in the results with the REBO+LJ and Tersoff potentials in the armchair direction. In order to explain this behavior, Figs. 5a and 5b show the distribution of the bond lengths at $\lambda_{1}=1.25$ for stretch along the armchair and zigzag directions, respectively, employing the REBO+LJ potential. It is found that around 36$\%$ of bonds are stretched to a bond length of more than 0.17 nm when the stretch is along the armchair direction, which activates the cutoff function of Eq. (32). Due to the discontinuity of the second derivative of the cutoff function a sudden rise in the stress is recorded999REBO+LJ is sensitive to the choice of cutoff distances. Other values have shown to lead to bond breaking (and hence a sudden stress drop) at much lower strains.. In the case of the Tersoff potential, the anomalous response in $\sigma_{11}$ in the armchair direction around $\lambda_{1}=1.2$ is due to the elongation of bonds beyond the cutoff radii 0.18 and 0.185 nm for the Tersoff and modified Tersoff potentials, respectively, see Eq. (37).

2. *Stress perpendicular to the stretch direction.

*The variation in the lateral stress $\sigma_{22}$ – the stress in the perpendicular direction to the stretch – and its error according to Eq. (26) as a function of $\lambda_{1}$ are shown in Figs. 6 and 7. The figures show that the results from the MM3 potential agree well with the DFT results up to $\lambda_{1}=1.15$ and $\lambda_{1}=1.1$ in the armchair and zigzag directions, respectively. As Fig. 6 shows for stretch in the armchair direction, the lateral stresses from the REBO+LJ potential agree well with the DFT results within the small deformation regime (upto $\lambda_{1}=1.025$) and then show gradual softening. However, for stretch in the zigzag direction (see Fig. 7), $\sigma_{22}$ from REBO+LJ matches well with $\sigma_{22}$ from DFT up to $\lambda_{1}=1.25$. Contrary to the results from the MM3 and REBO+LJ potentials, both Tersoff potentials produce negative lateral stress for uniaxial stretch.

To explain this exceptional behavior of the Tersoff potential, we have performed simulations without the constraints on the lateral edge atoms and calculated the Poisson’s ratio $\nu=-{\varepsilon_{\text{lat}}}/{\varepsilon_{\text{lon}}},$ where $\varepsilon_{\text{lat}}$ and $\varepsilon_{\text{lon}}$ are the lateral and longitudinal strains, respectively. The strains are computed by taking the ratio of the change in periodic box dimensions with the initial box dimensions. For stretch along the zigzag and armchair directions, the SLGS exhibits a negative Poisson’s ratio for all the stretch ratios according to the Tersoff potential. This has also been reported in the literature [76, 77, 41]. Due to this behavior, the Tersoff potential stands in sharp contrast to all the other methods.

It is noticed that the stress response computed from the continuum model for both the loading cases are in good agreement with those from DFT. Both the molecular and continuum stresses are equal at smaller stretches, and it is evident that the SLGS shows anisotropy at higher stretches. Additionally, the anisotropy predicted from both the approaches is different, i.e., molecular simulations predict that the armchair direction is stiffer whereas DFT and continuum predict otherwise.

To elucidate the anisotropy of the SLGS, we have plotted the two dominating energy terms101010Van der Waals interactions play an insignificant role here. Including them reveals that they are below 0.06 N/m, which is less than 1$\%$ of the bond-stretching term. of the MM3 potential in Fig. 8. The figure shows that the energy contribution from bond-stretching (term $U_{s}$ in Eq. (5)) is the same for stretch in both directions.

The contribution of angle-bending on the other hand (term $U_{\theta}$) to the total energy for both cases is equal up to $\lambda_{1}=1.17$, but deviates beyond this stretch. It increases more strongly when the SLGS is stretched in the armchair direction. The remaining energy contributions in Eq. (5) are not significant for this deformation state. Thus, we conclude that $U_{\theta}$ is responsible for the anisotory in SLGS at large deformations.

3. *Surface tension.

*Finally, the variation in the surface tension $\gamma=(\sigma_{11}+\sigma_{22})/2$, under pure dilatation and its corresponding percentage error is shown in Fig. 9.

As seen, the results from molecular simulation employing the MM3 potential agree with the results from DFT up to $\mathcal{J}_{1}=0.28$ and then show gradual hardening, whereas the results from the REBO+LJ and Tersoff potentials agree only up to $\mathcal{J}_{1}=0.08$ and then deviate. Additionally, the results from the later two potentials exhibit a sharp rise at $\mathcal{J}_{1}=0.39$ and 0.43, respectively. At those values of $\mathcal{J}_{1}$, $\approx 96\%$ of bonds are stretched to more than 0.17 nm and 0.18 nm ($r_{IJ}^{\mathrm{min}}$ of respective potentials) using the REBO+LJ and Tersoff potentials, respectively. On the other hand, the results from the continuum model and the modified Tersoff potential agree within an error of $\approx 5\%$.

4.3 Out-of-plane bending of SLGS

In molecular simulations, the bending stiffness of the SLGS can be calculated in two different ways. They both assume linear elastic bending behavior.

In the first approach, the SLGS is considered as a thin, linearly elastic, homogeneous and isotropic plate with thickness $h$, whose bending stiffness is given by [78]

[TABLE]

where $E$ and $\nu$ are the Young’s modulus and Poisson’s ratio, respectively. To avoid introducing a thickness of the SLGS in the first approach, the first bending frequency of a square SLGS is computed instead in our work. In MS simulations, the Vibrate module in Tinker is used to find the natural frequencies. In MD, on the other hand, the equilibrated SLGS is deformed into the mode shape corresponding to the first natural frequency and then allowed to vibrate freely in the NVE ensemble. The time history of the atoms closest to the center of the SLGS is then evaluated to determine the frequencies of the transverse vibration using fast Fourier transform (FFT). The obtained frequencies are then compared with those determined from an equivalent plate model. The formula for the transverse frequencies of a linearly elastic, homogeneous and isotropic square plate of side $a$ is given by [79]

[TABLE]

where $m$ and $n$ are the half wave numbers along the x and y directions, $c_{\mathrm{b}}$ and $\rho_{\mathrm{s}}$ are the bending stiffness and areal density, respectively. For $m=1$ and $n=1$, the constant $\lambda_{11}^{2}$ is 35.99 [79].

The bending stiffness determined from the MM3, REBO+LJ and Tersoff potentials, according to Eqs. (27) and (28) (see Table 5) is larger than 1.49 eV, which is the value from DFT [55, 72]. The discrepancy may be associated with the presence of pre-stress in the molecular simulations, which also causes higher frequencies at zero stretch. This is discussed further below. On the other hand, the bending stiffness determined from the modified Tersoff potential is 51 $\%$ lower than the DFT value.

In the second approach, the bending energy of the SLGS is obtained by computing the potential energy of relaxed carbon nanotubes of different radii with respect to the ground state energy of SLGS. The potential energy is fitted by a quadratic curve. The second derivative of this curve then corresponds to the bending stiffness. The obtained values are listed in Table 6. We note that the second approach is problematic, as relaxed CNTs usually change radius and therefore also contain in-plane strain energy that is usually not accounted for in the stiffness calculation. As a consequence the bending stiffness may be overestimated.

The variation of the bending energy with the curvature for armchair and zigzag carbon nanotubes is shown in Fig. 10a and 10b. The bending energy from the DFT and DFT-equivalent continuum model in (12) is given by $W_{\text{b}}=\displaystyle{c_{\text{b}}}\kappa_{1}^{2}/{2}~{},$ where $c_{\text{b}}$ is the bending stiffness and $\kappa_{1}$ is the curvature radius. The curvature radius of a CNT is the reciprocal of the radius $R=({\sqrt{3}a_{cc}}/{2\pi})\sqrt{n^{2}+m^{2}+nm}$ , where $a_{cc}$ is the equilibrium bond length of C-C, and $n$ and $m$ are the chirality indices.

The bending energy calculated from the MM3 potential differs more than the other two potentials considered. The difference may be attributed to the presence of higher order and cross-interaction terms of the potential. In their study, Gupta et al. [81] computed the bending stiffness of different radii and chiralities of CNTs by equating the frequencies from MS simulations employing the MM3 potential with that of a shell model. The bending energy computed from these calculations is shown in Fig. 10(b).

4.4 Modal analysis

In the following we discuss the effect of incremental prestretch in a SLGS and a CNC on their first few modes of vibrations.

4.4.1 Square SLGS

After having established the validity and accuracy of the MM3 potential for stretches up to $\approx 1.1$ in the previous section, we now study the transverse modes of vibrations of a SLGS using this potential. The vibration response is then compared with that obtained from the DFT-based continuum model at different stretch states. In the continuum model, the stretch is applied up to the loss of ellipticity of the elasticity tensor. This limiting value of stretch is also used in the molecular simulations.

The mode shapes and frequency variation of the transverse vibrations of the first three modes of the SLGS with increasing stretch/dilatation are shown in Fig. 11 and Fig. 12, respectively. The molecular simulations are performed without applying PBCs. As Fig. 12 shows, the frequency vs. stretch curves obtained from the two approaches agree within $\approx 95\%$. It is observed that for all the three loading cases, the frequencies increase monotonically with the stretch. The rate of increase in the frequency is higher at small stretches than at higher stretches. Singh and Patel [52] have reported a similar behaviour for a rectangular SLGS under uniaxial stretch. It is also observed that at zero strain, the frequencies from the molecular simulations are $\approx 15\%$ higher than those computed from the continuum model. This difference may be due to residual stresses in the relaxed configuration of SLGS, which are not accounted for in the present continuum model. For uniaxial stretch along the armchair direction, the frequencies from the molecular simulations and those from the continuum model agree up to $\lambda_{1}=1.15$ after which the MS simulation results show a stiffer behavior. This is consistent with the variation in the stress with stretch shown in Figs. 3 and 4. At higher stretches, the variation in the frequency is higher when the SLGS is stretched along the armchair direction than in the zigzag direction. This mild anisotropic response is attributed to the angle-bending energy, which contributes significantly to the total potential energy when the SLGS is stretched along the armchair direction. Contrary to the molecular simulation results, the continuum model shows the zigzag direction to be stiffer111111I.e. the mode shapes with a higher number of sinusoidal waves along the stiffer direction have a higher frequencies than the one with higher waves along the softer direction..

4.4.2 Carbon nanocone

The modal analysis of a pre-stretched simply-supported CNC is carried out using MS simulations employing the MM3 potential. The modal frequencies obtained from the MS simulations are compared with those obtained from the continuum model at each stretch. For this study, a truncated cone is selected with initial apex angle $\alpha=14.2^{\circ}$, height $H=12.4$ nm, and tip radius $R=1.017$ nm.

Before computing the modal frequencies in the MS simulations, the minimum energy configuration of the structure is found at each stretch level. The stretch is obtained by fixing all the degrees-of-freedom of the edge atoms of the larger radius and incrementally moving the atoms of the smaller radius by 0.1 Å  in the axial direction. We note that this leads to an inhomogeneous stretch state in the CNC, due to its tapered geometry. The following results are therefore taken at the average stretch $\langle\lambda_{1}\rangle=h/H$, where $h$ is the current height. The variations in the potential energy and the modal frequencies computed from the MM3 potential and the continuum model are shown in Fig. 13. The corresponding mode shapes computed from the MM3 potential are shown in Fig. 14.

For moderate values of the stretch, the potential energy from both the approaches are in good agreement. As in the case of SLGS, the modal frequencies from both the approaches differ slightly in the unstretched configuration, possibly due to the absence of residual stresses in the continuum model. The frequency of the lower order modes $\omega_{01}$ and $\omega_{21}$ obtained from the MM3 potential and the continuum model are in good agreement. Based on continuum shell theory, we conjecture that this agreement is due to the dominant membrane/stretching deformation in the lower modes for which the present continuum model is calibrated. This is also the reason for a good match in the strain energy since for moderate value of $\langle\lambda_{1}\rangle$ only the membrane deformations are dominating. At higher modes bending deformation dominates and hence there is a disagreement between the frequencies computed from the two approaches. Finally, we note that $\omega_{01}$ remains unaffected for the range of $\langle\lambda_{1}\rangle$ considered here.

4.5 Summary of main findings

While comparing the elastic response of SLGS obtained from the three potentials and the continuum model with that from DFT, we find the following:

1. The MM3 potential results agree up to a stretch of 1.1 (10$\%$ strain).

2. The REBO+LJ potential results agree only within the range of small deformation.

3. The two different parameterizations of the Tersoff potential fail to capture the behavior correctly. Further, they predict negative Poisson ratio.

4. The agreement of the energies and stresses for uniaxial and biaxial stretch between the continuum model and the DFT model are within an error of $\approx 0.05\%$.

While examining the vibrations of stretched SLGS and CNC studied using MM3 potential and DFT based continuum model, we find that:

5. The transverse vibrational frequencies increase monotonically with the pre-stretch.

6. The frequencies found in MS simulations differ from those obtained using the continuum model at zero stretch.

7. The bending frequency $\omega_{01}$ of CNC is nearly unaffected by the range of stretch studied in the paper.

5 Conclusions

The nonlinear mechanical response of square SLGS under large uniaxial stretch and pure dilatation is investigated with continuum and atomistic simulations employing the MM3, REBO+LJ, and Tersoff potentials. The obtained results are compared with DFT results available in the literature for the two stretch states. The elastic response obtained from the continuum model is in good agreement with DFT results. Among the potentials used for molecular simulations, the MM3 is found to be the most accurate in this study. The transverse vibrational frequencies of a square SLGS and a CNC found from the continuum model and molecular simulations employing the MM3 potential are found to be in good agreement, within an error of about $~{}5\%$. In future work, the present study will be extended to finite temperatures.

Acknowledgement

Financial support from the German Research Foundation (DFG) through grants SA 1822/8-1 and GSC 111 is gratefully acknowledged. The authors thank Farzad Shirazian for his comments on the atomistic simulations.

Appendix A Description of the molecular potentials

In this section, the mathematical expressions of the used potentials are given.

A.1 MM3 Potential

The terms in Eq. (5) are given for atoms $I$, $J$, $K$ and $L$ by

[TABLE]

[TABLE]

where $r_{IJ}$, $\theta_{IJK}$ and $\phi_{IJKL}$ are the bond length, the angle between the bonds and the torsion angle, respectively. The parameters with superscript e define the equilibrium values at which the total potential energy is minimum. $K_{\text{s}}$, $K_{\theta}$, $V_{1}$, $V_{2}$, $V_{3}$, $\epsilon_{\text{e}}$, $r_{\text{v}}$, $K_{\text{s}\theta}$, $K_{\phi s}$, $K_{\theta\theta^{{}^{\prime}}}$, $r^{\text{e}}$ and $\theta^{\text{e}}$ are the potential parameters and are listed in Table 7.

A.2 REBO+LJ Potential

The repulsive ($E_{\text{R}}$) and attractive ($E_{\text{A}}$) terms in Eq. (6) are

[TABLE]

[TABLE]

where $Q$, $A$, $\alpha$, $B^{n}$ and $\beta^{n}$ are the potential parameters that depend on the type of atoms $I$ and $J$ and are given in Table 8. $w$ is the bond-weighting factor that depends on the cutoff distances (${r^{\mathrm{min}}}$ and ${r^{\mathrm{max}}}$) and varies from zero to one. The REBO interactions are gradually turned off, when the bond length is in the range ${r^{\mathrm{min}}}<r_{IJ}<{r^{\mathrm{max}}}$. This is achieved through the bond-weighting factor

[TABLE]

where the switching and scaling functions $S^{{}^{\prime}}$ and $t_{c}$

[TABLE]

[TABLE]

Here $\Theta(t)$ is the Heaviside step function.

A.3 Tersoff Potential

The terms in the Tersoff potential of Eq. (8) are given by

[TABLE]

[TABLE]

and

[TABLE]

where $A$, $B$, $\lambda$ and $\mu$ are the parameters that depend on paired atom types.

The cutoff function ($f_{\text{c}}$) describes a gradual decrease of the bond strength between ${r^{\mathrm{min}}}$ and ${r^{\mathrm{max}}}$. In Eq. (8), $b_{IJ}$ is the bond order term. It depends on neighboring atoms and it is defined by

[TABLE]

[TABLE]

[TABLE]

where $\theta_{IJK}$ is the angle between bond $IJ$ and $IK$, while $A$, $B$, $\lambda$, $\mu$, $\beta^{n}$, $n$, $c$, $d$, $h$, ${r^{\text{min}}}$, ${r^{\text{max}}}$, $\chi$ and $\omega$ are the potential parameters that for carbon are given in Table 9.