An atomistic-based F\"oppl-von K\'arm\'an model for graphene

TL;DR

This paper develops a non-linear continuum F"oppl-von Kármán model for graphene based on atomistic interactions, capturing out-of-plane displacements and lattice effects with derived constitutive coefficients.

Contribution

It introduces a novel atomistic-to-continuum derivation of a F"oppl-von Kármán model for graphene, accounting for lattice shifts and self-stress effects.

Findings

Derived continuum model depends on macroscopic displacement and lattice shift.

Constitutive coefficients explicitly linked to atomistic Brenner's REBO potential.

Formal minimization yields a F"oppl-von Kármán type model for graphene.

Abstract

We deduce a non-linear continuum model of graphene for the case of finite out-of-plane displacements and small in-plane deformations. On assuming that the lattice interactions are governed by the Brenner's REBO potential of 2nd generation and that self-stress is present, we introduce discrete strain measures accounting for up-to-the-third neighbor interactions. The continuum limit turns out to depend on an average (macroscopic) displacement field and a relative shift displacement of the two Bravais lattices that give rise to the hexagonal periodicity. On minimizing the energy with respect to the shift variable, we formally determine a continuum model of F\"oppl-von K\'arm\'a type, whose constitutive coefficients are given in terms of the atomistic interactions.

Equations

$$\begin{array}[]{l}L_{1}(\ell)=\{\mathbf{x}\in\mathbb{R}^{2}:\mathbf{x}=n^{1}\ell\mathbf{d}_{1}+n^{2}\ell\mathbf{d}_{2}\quad\mbox{with}\quad(n^{1},n^{2})\in\mathbb{Z}^{2}\},\\ L_{2}(\ell)=\ell\mathbf{p}+L_{1}(\ell),\end{array}$$

$${\mathbf{d}}_1=\sqrt{3}{\mathbf{e}}_1,{\mathbf{d}}_2=\frac{\sqrt{3}}{2}{\mathbf{e}}_1+\frac{3}{2}{\mathbf{e}}_2\text{mbox}and\mathbf{p}=\frac{\sqrt{3}}{2}{\mathbf{e}}_1+\frac{1}{2}{\mathbf{e}}_2.$$

$$\mathbf{p}_{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 267\relax$}}{\mbox{\boldmath$\textstyle\mathchar 267\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 267\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 267\relax$}}}=\mathbf{d}_{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 267\relax$}}{\mbox{\boldmath$\textstyle\mathchar 267\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 267\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 267\relax$}}}-\mathbf{p}\ \ ({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 267\relax$}}{\mbox{\boldmath$\textstyle\mathchar 267\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 267\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 267\relax$}}}=1,2)\quad\mbox{and}\quad\mathbf{p}_{3}=-\,\mathbf{p}.$$

$${\mathbf{d}}_3={\mathbf{d}}_2-{\mathbf{d}}_1.$$

$${\mathbf{x}}^{\ell }=n^1\ell {\mathbf{d}}_1+n^2\ell {\mathbf{d}}_2+m\ell \mathbf{p},(n^1,n^2,m)\in {\mathbb{Z}}^2\times \{0,1\}$$

$$V=\sum_{i}\sum_{j<i}V_{ij}\,,\quad V_{ij}=V_{R}(l_{ij})+b_{ij}({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 291\relax$}}{\mbox{\boldmath$\textstyle\mathchar 291\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 291\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 291\relax$}}}_{hij},\Theta_{hijk})V_{A}(l_{ij}),$$

$$\mathbf{u}^{\ell}(\mathbf{x}^{\ell})={{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}\operatorname{\boldsymbol{\mathfrak{v}}^{\ell}}(\mathbf{x}^{\ell})+{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}^{1/2}\operatorname{\mathfrak{w}^{\ell}}(\mathbf{x}^{\ell})\mathbf{e}_{3},$$

$$\displaystyle{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 270\relax$}}{\mbox{\boldmath$\textstyle\mathchar 270\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 270\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 270\relax$}}}l_{i}(\mathbf{x}^{\ell})$$

$$=\ell |{\mathbf{p}}_i+\frac{{\mathbf{u}}^{\ell }({\mathbf{x}}^{\ell }+\ell {\mathbf{p}}_i)-{\mathbf{u}}^{\ell }({\mathbf{x}}^{\ell })}{\ell }|-\ell .$$

$$D_{{\mathbf{p}}_i}^{\ell }f(\mathbf{x}):=\frac{f(\mathbf{x}+\ell {\mathbf{p}}_i)-f(\mathbf{x})}{\ell },$$

$$\frac{{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 270\relax$}}{\mbox{\boldmath$\textstyle\mathchar 270\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 270\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 270\relax$}}}\ell_{i}}{\ell_{i}}=|\mathbf{p}_{i}+D^{\ell}_{\mathbf{p}_{i}}\mathbf{u}^{\ell}(\mathbf{x}^{\ell})|-1=|\mathbf{p}_{i}+{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}D^{\ell}_{\mathbf{p}_{i}}\operatorname{\boldsymbol{\mathfrak{v}}^{\ell}}(\mathbf{x}^{\ell})+{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}^{1/2}D^{\ell}_{\mathbf{p}_{i}}\operatorname{\mathfrak{w}^{\ell}}(\mathbf{x}^{\ell})\mathbf{e}_{3}|-1,$$

$$\frac{{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 270\relax$}}{\mbox{\boldmath$\textstyle\mathchar 270\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 270\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 270\relax$}}}\ell_{i}}{\ell_{i}}={{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}\left(D^{\ell}_{\mathbf{p}_{i}}\operatorname{\boldsymbol{\mathfrak{v}}^{\ell}}(\mathbf{x}^{\ell})\cdot\mathbf{p}_{i}+\frac{1}{2}|D^{\ell}_{\mathbf{p}_{i}}\operatorname{\mathfrak{w}^{\ell}}(\mathbf{x}^{\ell})|^{2}\right)+o({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}),$$

$${{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 290\relax$}}{\mbox{\boldmath$\textstyle\mathchar 290\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 290\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 290\relax$}}}_{i}(\mathbf{x}^{\ell}):=D^{\ell}_{\mathbf{p}_{i}}\mathbf{v}^{\ell}(\mathbf{x}^{\ell})\cdot\mathbf{p}_{i}+\frac{1}{2}|D^{\ell}_{\mathbf{p}_{i}}w^{\ell}(\mathbf{x}^{\ell})|^{2}\qquad i=1,2,3,$$

$${\mathbf{u}}_i^{\ell }:={\mathbf{u}}^{\ell }({\mathbf{x}}^{\ell }+\ell {\mathbf{p}}_i),\text{mbox}and{\mathbf{u}}_0^{\ell }:={\mathbf{u}}^{\ell }({\mathbf{x}}^{\ell }).$$

$$\mathbf{q}_{i+1}=\ell\left(\mathbf{p}_{i+1}+{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}D^{\ell}_{\mathbf{p}_{i+1}}\operatorname{\boldsymbol{\mathfrak{v}}^{\ell}}(\mathbf{x}^{\ell})+{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}^{1/2}D^{\ell}_{\mathbf{p}_{i+1}}\operatorname{\mathfrak{w}^{\ell}}(\mathbf{x}^{\ell})\mathbf{e}_{3}\right)$$

$$\mathbf{q}_{i+2}=\ell\left(\mathbf{p}_{i+2}+{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}D^{\ell}_{\mathbf{p}_{i+2}}\operatorname{\boldsymbol{\mathfrak{v}}^{\ell}}(\mathbf{x}^{\ell})+{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}^{1/2}D^{\ell}_{\mathbf{p}_{i+2}}\operatorname{\mathfrak{w}^{\ell}}(\mathbf{x}^{\ell})\mathbf{e}_{3}\right)$$

$$\operatorname{\overset{(2)}{{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 291\relax$}}{\mbox{\boldmath$\textstyle\mathchar 291\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 291\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 291\relax$}}}}}_{i}({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}})=\arccos\left(\frac{{\mathbf{q}_{i+1}}\cdot{\mathbf{q}_{i+2}}}{|{\mathbf{q}_{i+1}}||{\mathbf{q}_{i+2}}|}\right).$$

$$\frac{1}{\ell^{2}}\mathbf{q}_{i+1}\cdot\mathbf{q}_{i+2}=\mathbf{p}_{i+1}\cdot\mathbf{p}_{i+2}+{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}\Big{(}\mathbf{p}_{i+1}\cdot D^{\ell}_{\mathbf{p}_{i+2}}\operatorname{\boldsymbol{\mathfrak{v}}^{\ell}}(\mathbf{x}^{\ell})+\mathbf{p}_{i+2}\cdot D^{\ell}_{\mathbf{p}_{i+1}}\operatorname{\boldsymbol{\mathfrak{v}}^{\ell}}(\mathbf{x}^{\ell})+D^{\ell}_{\mathbf{p}_{i+1}}\operatorname{\mathfrak{w}^{\ell}}(\mathbf{x}^{\ell})D^{\ell}_{\mathbf{p}_{i+2}}\operatorname{\mathfrak{w}^{\ell}}(\mathbf{x}^{\ell})\Big{)}+o({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}),$$

$$\displaystyle\frac{1}{\ell^{2}}|\mathbf{q}_{i+1}|^{2}=1+{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}\big{(}|D^{\ell}_{\mathbf{p}_{i+1}}\operatorname{\mathfrak{w}^{\ell}}(\mathbf{x}^{\ell})|^{2}+2\mathbf{p}_{i+1}\cdot D^{\ell}_{\mathbf{p}_{i+1}}\operatorname{\boldsymbol{\mathfrak{v}}^{\ell}}(\mathbf{x}^{\ell})\big{)}+o({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}),$$

$$\displaystyle\frac{1}{\ell^{2}}|\mathbf{q}_{i+2}|^{2}=1+{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}\big{(}|D^{\ell}_{\mathbf{p}_{i+2}}\operatorname{\mathfrak{w}^{\ell}}(\mathbf{x}^{\ell})|^{2}+2\mathbf{p}_{i+2}\cdot D^{\ell}_{\mathbf{p}_{i+2}}\operatorname{\boldsymbol{\mathfrak{v}}^{\ell}}(\mathbf{x}^{\ell})\big{)}+o({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}),$$

$$\operatorname{\overset{(2)}{{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 291\relax$}}{\mbox{\boldmath$\textstyle\mathchar 291\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 291\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 291\relax$}}}}}_{i}({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}})=\arccos\left(\frac{\mathbf{p}_{i+1}\cdot\mathbf{p}_{i+2}+a{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}+o({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}})}{\sqrt{1+b{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}+o({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}})}\,\sqrt{1+c{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}+o({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}})}}\right)=\frac{2}{3}{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 281\relax$}}{\mbox{\boldmath$\textstyle\mathchar 281\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 281\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 281\relax$}}}-\frac{1}{2\sqrt{3}}(4a+b+c)\,{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}+o({{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 280\relax$}}{\mbox{\boldmath$\textstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 280\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 280\relax$}}}),$$

$$a:={\mathbf{p}}_{i+1}\cdot D_{{\mathbf{p}}_{i+2}}^{\ell }{\mathfrak{v}}^{\mathrm{\ell }}({\mathbf{x}}^{\ell })+{\mathbf{p}}_{i+2}\cdot D_{{\mathbf{p}}_{i+1}}^{\ell }{\mathfrak{v}}^{\mathrm{\ell }}({\mathbf{x}}^{\ell })+D_{{\mathbf{p}}_{i+1}}^{\ell }{\mathfrak{w}}^{\mathrm{\ell }}({\mathbf{x}}^{\ell })D_{{\mathbf{p}}_{i+2}}^{\ell }{\mathfrak{w}}^{\mathrm{\ell }}({\mathbf{x}}^{\ell }),$$

$$b:=|D_{{\mathbf{p}}_{i+1}}^{\ell }{\mathfrak{w}}^{\mathrm{\ell }}({\mathbf{x}}^{\ell }){|}^2+2{\mathbf{p}}_{i+1}\cdot D_{{\mathbf{p}}_{i+1}}^{\ell }{\mathfrak{v}}^{\mathrm{\ell }}({\mathbf{x}}^{\ell }),$$

$$c:=|D_{{\mathbf{p}}_{i+2}}^{\ell }{\mathfrak{w}}^{\mathrm{\ell }}({\mathbf{x}}^{\ell }){|}^2+2{\mathbf{p}}_{i+2}\cdot D_{{\mathbf{p}}_{i+2}}^{\ell }{\mathfrak{v}}^{\mathrm{\ell }}({\mathbf{x}}^{\ell }).$$

$$\operatorname{\overset{(2)}{{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 288\relax$}}{\mbox{\boldmath$\textstyle\mathchar 288\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 288\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 288\relax$}}}}}_{i}(\mathbf{x}^{\ell})=\operatorname{\overset{(2)}{{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 291\relax$}}{\mbox{\boldmath$\textstyle\mathchar 291\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 291\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 291\relax$}}}}}_{i}-\frac{2}{3}{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 281\relax$}}{\mbox{\boldmath$\textstyle\mathchar 281\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 281\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 281\relax$}}}=\\ -\frac{1}{\sqrt{3}}\Bigg{(}D^{\ell}_{\mathbf{p}_{i+1}}\mathbf{v}^{\ell}(\mathbf{x}^{\ell})\cdot\mathbf{p}_{i+1}+D^{\ell}_{\mathbf{p}_{i+2}}\mathbf{v}^{\ell}(\mathbf{x}^{\ell})\cdot\mathbf{p}_{i+2}+2(D^{\ell}_{\mathbf{p}_{i+1}}\mathbf{v}^{\ell}(\mathbf{x}^{\ell})\cdot\mathbf{p}_{i+2}+D^{\ell}_{\mathbf{p}_{i+2}}\mathbf{v}^{\ell}(\mathbf{x}^{\ell})\cdot\mathbf{p}_{i+1})+\\ \frac{1}{2}\left(|D^{\ell}_{\mathbf{p}_{i+1}}w^{\ell}(\mathbf{x}^{\ell})|^{2}+|D^{\ell}_{\mathbf{p}_{i+2}}w^{\ell}(\mathbf{x}^{\ell})|^{2}\right)+2D^{\ell}_{\mathbf{p}_{i+1}}w^{\ell}(\mathbf{x}^{\ell})\,D^{\ell}_{\mathbf{p}_{i+2}}w^{\ell}(\mathbf{x}^{\ell})\Bigg{)}\,,$$

$$\sum_{i=1}^{3}\operatorname{\overset{(2)}{{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 288\relax$}}{\mbox{\boldmath$\textstyle\mathchar 288\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 288\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 288\relax$}}}}}_{i}(\mathbf{x}^{\ell})=-3\sqrt{3}\Bigg{(}\frac{1}{3}\sum_{i=1}^{3}D^{\ell}_{\mathbf{p}_{i}}w^{\ell}(\mathbf{x}^{\ell})\Bigg{)}^{2}.$$

$$\overset{(1)}{{{}\mathchoice{\mbox{\boldmath$\displaystyle\mathchar 288\relax$}}{\mbox{\boldmath$\textstyle\mathchar 288\relax$}}{\mbox{\boldmath$\scriptstyle\mathchar 288\relax$}}{\mbox{\boldmath$\scriptscriptstyle\mathchar 288\relax$}}}}_{i}(\mathbf{x}^{\ell})=\frac{1}{\sqrt{3}}\Bigg{(}D^{\ell}_{-\mathbf{p}_{i+1}}\mathbf{v}^{\ell}(\mathbf{x}^{\ell})\cdot\mathbf{p}_{i+1}+D^{\ell}_{-\mathbf{p}_{i+2}}\mathbf{v}^{\ell}(\mathbf{x}^{\ell})\cdot\mathbf{p}_{i+2}+2(D^{\ell}_{-\mathbf{p}_{i+1}}\mathbf{v}^{\ell}(\mathbf{x}^{\ell})\cdot\mathbf{p}_{i+2}+D^{\ell}_{-\mathbf{p}_{i+2}}\mathbf{v}^{\ell}(\mathbf{x}^{\ell})\cdot\mathbf{p}_{i+1})+\\ -\frac{1}{2}\left(|D^{\ell}_{-\mathbf{p}_{i+1}}w^{\ell}(\mathbf{x}^{\ell})|^{2}+|D^{\ell}_{-\mathbf{p}_{i+2}}w^{\ell}(\mathbf{x}^{\ell})|^{2}\right)-2D^{\ell}_{-\mathbf{p}_{i+1}}w^{\ell}(\mathbf{x}^{\ell})\,D^{\ell}_{-\mathbf{p}_{i+2}}w^{\ell}(\mathbf{x}^{\ell})\Bigg{)}\,.$$

$$\cos {\stackrel{(\mathrm{c})}{\mathrm{\Theta }}}_{{\mathbf{p}}_i^{+}}=\frac{({\mathbf{q}}_i\times {\mathbf{q}}_{i+1})\cdot ({\mathbf{q}}_i\times {\mathbf{q}}_{i^{+}})}{|{\mathbf{q}}_i\times {\mathbf{q}}_{i+1}||{\mathbf{q}}_i\times {\mathbf{q}}_{i^{+}}|},$$

$$\cos {\stackrel{(\mathrm{c})}{\mathrm{\Theta }}}_{{\mathbf{p}}_i^{-}}=\frac{({\mathbf{q}}_{i+2}\times {\mathbf{q}}_i)\cdot ({\mathbf{q}}_{i^{-}}\times {\mathbf{q}}_i)}{|{\mathbf{q}}_{i+2}\times {\mathbf{q}}_i||{\mathbf{q}}_{i^{-}}\times {\mathbf{q}}_i)|},$$

$$\cos {\stackrel{(\mathrm{z})}{\mathrm{\Theta }}}_{{\mathbf{p}}_i{\mathbf{p}}_{i+1}}=\frac{({\mathbf{q}}_i\times {\mathbf{q}}_{i+1})\cdot ({\mathbf{q}}_{i^{-}}\times {\mathbf{q}}_i)}{|{\mathbf{q}}_i\times {\mathbf{q}}_{i+1}||{\mathbf{q}}_{i^{-}}\times {\mathbf{q}}_i|},$$

$$\cos {\stackrel{(\mathrm{z})}{\mathrm{\Theta }}}_{{\mathbf{p}}_i{\mathbf{p}}_{i+2}}=\frac{({\mathbf{q}}_{i+2}\times {\mathbf{q}}_i)\cdot ({\mathbf{q}}_i\times {\mathbf{q}}_{i^{+}})}{|{\mathbf{q}}_{i+2}\times {\mathbf{q}}_i||{\mathbf{q}}_i\times {\mathbf{q}}_{i^{+}}|},$$

$${\mathbf{q}}_{i^{+}}=$$

$$=$$

$${\mathbf{q}}_{i^{-}}=$$

$$=$$

Full text

**An atomistic-based

Föppl–von Kármán model for graphene **

Cesare Davini1

  

Antonino Favata2

  

Roberto Paroni3

1 Via Parenzo 17, 33100 Udine

cesare.davini@uniud.it

2 Department of Structural and Geotechnical Engineering

Sapienza University of Rome, Rome, Italy

antonino.favata@uniroma1.it

3 Dipartimento di Ingegneria Civile e Industriale

Università di Pisa, Pisa, Italy

roberto.paroni@unipi.it

Abstract

We deduce a non-linear continuum model of graphene for the case of finite out-of-plane displacements and small in-plane deformations. On assuming that the lattice interactions are governed by the Brenner’s REBO potential of 2nd generation and that self-stress is present, we introduce discrete strain measures accounting for up-to-the-third neighbor interactions. The continuum limit turns out to depend on an average (macroscopic) displacement field and a relative shift displacement of the two Bravais lattices that give rise to the hexagonal periodicity. On minimizing the energy with respect to the shift variable, we formally determine a continuum model of Föppl–von Kármán type, whose constitutive coefficients are given in terms of the atomistic interactions.

Keywords: Graphene, Continuum Modeling, Föppl–von Kármán equations.

Contents

1. 1 Introduction

2. 2 Kinematics and energetics

3. 3 Approximated strain measures

4. 3.1 Change of the edge lengths

5. 3.2 Change of the wedge angles

6. 3.3 Change of the dihedral angles

7. 4 Membrane and bending energy

8. 5 The continuum limit of the membrane energy

9. 6 Conclusions

1 Introduction

For its extraordinary mechanical, electrical, and thermal properties graphene is one of the most studied materials of the last two decades. Its discovery gave the occasion to renew a classical debate on the stability of 2D materials in nature [21, 25], and opened the way to the discussion of many theoretical issues. From the onset it was clear that the applications in various fields of technology might have been revolutionary. Though it’s fair to say that graphene’s potentialities are far from being fully explored and exploited, and remain the object of an intensive study, as much as it would be difficult to give an account of the huge literature on the subject. For a general picture of the subject and the new technological applications see the review by Ferrari et al. [13].

The availability of macroscopic models is crucial to design applications and experiments. The simplest models of structural mechanics such as membranes, plates, and shells have been often adopted in the past; in some cases, they have been assumed as a priori models and the relevant constitutive constants have been estimated from ab initio or Molecular Dynamics simulations. Huang et al. [18] highlight how this has resorted to a bit of a stretch in some cases and led to paradoxical conclusions. More recent contributions start from atomistic analyses based on appropriate constitutive assumptions on the interatomic potentials to obtain the continuum models of structural mechanics.

Thus, Lu and Huang [23] estimated the elastic modulus and bending stiffness of a graphene sheet from Molecular Mechanics calculations by considering the one-dimensional stretching and the rolling on cylinders of various radii of a rectangular piece of graphene, and assuming that the interatomic forces are ruled by Brenner’s REBO potential of the second type. In this way they obtained values of the elastic constant and bending stiffness that closely agree with those found by Kudin et al. [20] from ab initio calculations. These values are twice as much as those found by Arroyo and Belitschko [2] in a paper where the dihedral contribution in Brenner’s potential is taken into account; Arroyo and Belitschko also gave an atomistic-based membrane model for single layer crystalline films, [1]. Finally, Davini [4] deduced a 2D continuum model for the in-plane deformations of a graphene sheet within the framework of $\Gamma$-convergence.

The out-of-plane deformations have been considered by various authors. In particular, by exploiting a formal analysis we deduced a continuum model of a graphene sheet [5], and provided explicit expressions for the bending and Gaussian stiffness by starting from the study of the lattice kinematics and assuming the reactive empirical bond-order potential (REBO) of 2nd generation by Brenner et al. [3]. The approach takes into account the role of self-stress and provides a quantitative estimate of the self-stress contribution to the overall bending and Gaussian stiffness. Indeed, the continuum model turns out to be the $\Gamma$-limit of the discrete graphene sheet, as proven in [8].

To understand the bending behavior of graphene is of the essence for several technological applications. The bending behavior controls the ripple formation and the performance of graphene nano-electro-mechanical devices [18, 24, 37, 22, 31, 19, 30, 17, 16, 35, 28, 11, 12, 36, 14], and it is regarded as crucial in order to produce efficient hydrogen-storage devices [32, 15, 33]; moreover it can be instrumental to get inspiration for designing new metamaterials [7]. Indeed, the intrinsic ripples are believed to be essential for the structural stability of the 2D graphene lattice and may have major impacts on the electronic and mechanical properties of graphene [23].

In a recent review on Materials Today, Deng and Berry [9] give an overview of the hot problem of wrinkling, rippling, and crumpling, highlighting both formation mechanism and applications. The formation of these corrugations may have various explanations, see [10, 27, 26]. Basically, the out-of-plane deformations (wrinkles and ripples) can significantly reduce the magnitude of in-plane stresses generated, for instance, by defects, [29, 34, 38]. Zhang et al. [38] adopted a generalized Föppl–von Kármán equation for a flexible solid membrane to describe ripples near defects such as disclinations (heptagons or pentagons) and dislocations (heptagon-pentagon dipoles) on graphene, and predicted the large scale graphene configurations under specific defect distributions. The paper closely follows a study of Seung and Nelson [29]. Comparison with atomistic simulations indicates that the proposed model is capable to predict the atomic scale wrinkles near disclination/dislocation cores. The analysis shows that considering the buckling into 3-dimensional deformations is energetically more favorable than restricting to the in-plane ones. Similar defect-guided ripples in graphene were also simulated and discussed in the work of Wang et al. [34].

With an eye toward wrinkling and ripple formation, here we deduce a continuum model of graphene for the case of finite out-of-plane displacements and small in-plane deformations. We consider an array of C-atoms sitting at the nodes of a hexagonal lattice, and assume that the lattice interactions are governed by the Brenner’s REBO potential of 2nd generation and that self-stress is present. Thus, the starting point is the same as in [5], but the changes of edge lengths, wedge angles and dihedral angles are calculated by keeping the quadratic term in the out-of-plane displacements, according to the form of the in-plane Green-Lagrange strain used in Föppl–von Kármán plate theory.

The computation of the approximated measures of strain is done in Section 3. Unlike Zhang et al., [38], that assume a triangular lattice for the continuum analysis as done by Seung and Nelson in [29], here we use the real geometry of a hexagonal lattice. With due modifications, we adopt a harmonic approximation of the interatomic potential, which yields a splitting of the energy into membrane and bending parts, see Section 4. It follows that the bending part keeps the form already discussed in [5], while the membrane part turns out to be affected by the non-linearity of the assumed in-plane Green-Lagrange strain. The continuum limit of the membrane energy is computed in Section 5 according to the formal approach followed in [5]. The founding assumption is that, to within a remainder tending to zero with the lattice size, the nodal displacements can be described by an average (macroscopic) displacement and a relative shift displacement of the two Bravais lattices that give rise to the hexagonal periodicity. On minimizing the energy with respect to the shift variable, we formally determine a continuum model of Föppl–von Kármán type, whose constitutive coefficients are given in terms of the atomistic interactions. A full validation of the obtained continuum limit within the scheme of $\Gamma$-convergence is left for future work.

2 Kinematics and energetics

As reference configuration we use the $2$–lattice generated by two simple Bravais lattices

[TABLE]

simply shifted with respect to one another by $\ell\mathbf{p}$, see Fig. 1.

In (1), $\ell$ denotes the lattice size (the reference interatomic distance), while \ell\mathbf{d}_{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 267\relax}}{\mbox{\boldmath\textstyle\mathchar 267\relax}}{\mbox{\boldmath\scriptstyle\mathchar 267\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 267\relax}}}} and $\ell\mathbf{p}$ respectively are the lattice vectors and the shift vector, with

[TABLE]

The sides of the hexagonal cells in Figure 1 stand for the bonds between pairs of next nearest neighbor atoms and are represented by the vectors

[TABLE]

For convenience we also set

[TABLE]

In what follows we denote by

[TABLE]

the lattice points and label them by the triplets $(n^{1},n^{2},m)$: the points with $m=0$ belong to $L_{1}(\ell)$, while those in $L_{2}(\ell)$ correspond to $m=1$.

Graphene energetics depends on the description chosen to mimic atomic interactions. Our model is based on the 2nd-generation Brenner potential [3], which is one of the most used in molecular dynamics simulations of graphene. Accordingly, the binding energy $V$ of an atomic aggregate is given as a sum over nearest neighbors:

[TABLE]

where the individual effects of the repulsion and attraction functions $V_{R}(l_{ij})$ and $V_{A}(l_{ij})$, which model pair-wise interactions of the atoms $i$ and $j$ depending on their distance $l_{ij}$, are modulated by the bond-order function $b_{ij}$; for a given bond chain $h,i,j,k$ the function $b_{ij}$ depends in a complex manner on the angle between the edges $hi$ and $ij$ and on the dihedral angle between the planes spanned by $(hi,ij)$ and $(ij,jk)$. This potential reveals that, in order to properly account for the mechanical behavior of graphene, it is necessary to consider three types of energetic contributions:

binary interactions between next nearest atoms (edge bonds), 2. 2.

three-body interactions between consecutive pairs of next nearest atoms (wedge bonds), 3. 3.

and four-body interactions between three consecutive pairs of next nearest atoms (dihedral bonds). There are two types of relevant dihedral bonds: the Z-dihedra, in which the edges connecting the four atoms form a Z-shape, and the C-dihedra, in which the edges form a C-shape (see Fig. 2).

Moreover, it is possible to show [11] that the angle at ease between consecutive edges is greater than \frac{2}{3}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 281\relax}}{\mbox{\boldmath\textstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 281\relax}}}: this means that in the flat reference configuration the graphene sheet is not stress-free, and we will proper account for this feature.

3 Approximated strain measures

In this section we calculate the strain measures associated to a change of configuration described by a displacement field $\mathbf{u}^{\ell}:(L_{1}(\ell)\cup L_{2}(\ell))\cap\Omega\to\mathbb{R}^{3}$. Since we have in mind to deduce a model with the same non-linearities as the Föppl–von Kármán one, we write the displacements of the nodes in the form

[TABLE]

where $\textstyle\mathchar 280\relax$ is a positive scalar measuring smallness, \operatorname{\boldsymbol{\mathfrak{v}}^{\ell}}:={{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}^{-1}\mathbf{v}^{\ell}:={{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}^{-1}(\mathbf{e}_{1}\otimes\mathbf{e}_{1}+\mathbf{e}_{2}\otimes\mathbf{e}_{2})\mathbf{u}^{\ell} and \operatorname{\mathfrak{w}^{\ell}}:={{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}^{-1/2}:=w^{\ell}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}^{-1/2}\mathbf{u}^{\ell}\cdot\mathbf{e}_{3} stand for the in-plane and out-of-plane normalized displacements, respectively.

3.1 Change of the edge lengths

With {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}l_{i}(\mathbf{x}^{\ell}) we denote the change in length of the edge parallel to $\mathbf{p}_{i}$ and starting from the lattice point $\mathbf{x}^{\ell}\in(L_{1}(\ell)\cup L_{2}(\ell))\cap\Omega$. We fix our attention to lattice points in $L_{2}(\ell)$. Thus,

[TABLE]

On introducing the notation

[TABLE]

the axial strain measure can be recast as

[TABLE]

where we have made use of (6). The expansion up to the first order in $\textstyle\mathchar 280\relax$ of the non-linear strain measure (8) is

[TABLE]

which allows to define the edge strain measure, once rescaled-back by $\textstyle\mathchar 280\relax$:

[TABLE]

and {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}\ell_{i}={{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 290\relax}}{\mbox{\boldmath\textstyle\mathchar 290\relax}}{\mbox{\boldmath\scriptstyle\mathchar 290\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 290\relax}}}_{i}\ell_{i}.

3.2 Change of the wedge angles

For each fixed node $\mathbf{x}^{\ell}\in L_{2}(\ell))\cap\Omega$ we denote by \operatorname{\overset{(2)}{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}}}_{i}(\mathbf{x}^{\ell}) the angle of the wedge delimited by the edges $\mathbf{p}_{i+1}$ and $\mathbf{p}_{i+2}$; that is, the wedge angle opposite to the $i$-th edge (see Fig. 3). Here, $i,i+1$, and $i+2$ take values in $\{1,2,3\}$ and the sums should be interpreted mod 3: for instance, if $i=2$ then $i+1=3$ and $i+2=1$.

To keep the notation compact, we set

[TABLE]

Let

[TABLE]

and

[TABLE]

be the images of the edges parallel to $\mathbf{p}_{i+1}$ and $\mathbf{p}_{i+2}$ and starting at $\mathbf{x}^{\ell}$. Then, the angle \operatorname{\overset{(2)}{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}}}_{i}({{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}})= is given by

[TABLE]

Now, in the light of (11) and (12), we get

[TABLE]

and

[TABLE]

so that the expansion up to the first order in $\textstyle\mathchar 280\relax$ of \operatorname{\overset{(2)}{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}}}_{i}({{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}) yields

[TABLE]

with

[TABLE]

This allows to define the wedge strain measure, once the displacement is rescaled-back by $\textstyle\mathchar 280\relax$:

[TABLE]

for $\mathbf{x}^{\ell}\in L_{2}(\ell))\cap\Omega$. In particular, algebraic manipulations allow to conclude that

[TABLE]

If we consider a lattice point belonging to $L_{1}(\ell)$, it is not difficult to see that

[TABLE]

3.3 Change of the dihedral angles

For each fixed node $\mathbf{x}^{\ell}\in L_{2}(\ell)\cap\Omega$ and for each edge parallel to $\mathbf{p}_{i}$ and starting at $\mathbf{x}^{\ell}$ we need to define four types of dihedral angles $\operatorname{\overset{(c)}{\Theta}}_{\mathbf{p}_{i}^{+}}(\mathbf{x}^{\ell}),\operatorname{\overset{(c)}{\Theta}}_{\mathbf{p}_{i}^{-}}(\mathbf{x}^{\ell}),\operatorname{\overset{(z)}{\Theta}}_{\mathbf{p}_{i}\mathbf{p}_{i+1}}(\mathbf{x}^{\ell})$ and $\operatorname{\overset{(z)}{\Theta}}_{\mathbf{p}_{i}\mathbf{p}_{i+2}}(\mathbf{x}^{\ell})$:

[TABLE]

where

[TABLE]

are the images of vectors $\ell\mathbf{p}_{i^{+}}$ and $\ell\mathbf{p}_{i^{-}}$ (see Fig. 4, for $i=1$), parallel to $\mathbf{p}_{i+2}$ and $\mathbf{p}_{i+1}$ and starting at the image of the point $\mathbf{x}^{\ell}+\ell\mathbf{p}_{i}$.

Also here, $i,i+1$, and $i+2$ take values in $\{1,2,3\}$ and the sums should be interpreted mod 3: for instance, if $i=3$ then $i+1=1$ and $i+2=2$.

The C-dihedral angle $\operatorname{\overset{(c)}{\Theta}}_{\mathbf{p}_{i}^{+}}(\mathbf{x}^{\ell})$ is the angle corresponding to the C-dihedron with middle edge $\ell\mathbf{p}_{i}$ and oriented as $\mathbf{p}_{i}^{\perp}$, while $\operatorname{\overset{(c)}{\Theta}}_{\mathbf{p}_{i}^{-}}(\mathbf{x}^{\ell})$ is the angle corresponding to the C-dihedron oriented opposite to $\mathbf{p}_{i}^{\perp}$ (see Fig. 4 for $i=1$). The Z-dihedral angle $\operatorname{\overset{(z)}{\Theta}}_{\mathbf{p}_{i}\mathbf{p}_{i+1}}(\mathbf{x}^{\ell})$ corresponds to the Z-dihedron with middle edge $\ell\mathbf{p}_{i}$ and the other two edges parallel to $\mathbf{p}_{i+1}$ (see Fig. 4 for $i=1$).

To fix the ideas, we focus on the dihedral angle $\operatorname{\overset{(z)}{\Theta}}_{\mathbf{p}_{1}\mathbf{p}_{2}}$, sketched in Fig. 4; the other strains can be obtained in analogous manner.

The first order approximation of the dihedral angle is all we need to evaluate the corresponding energy contribution.

Let us introduce the vector

[TABLE]

image of $\mathbf{p}_{5}$ under the deformation. We have that

[TABLE]

Cumbersome computations yield:

[TABLE]

The expansion up to the first order in $\textstyle\mathchar 280\relax$ of $\operatorname{\overset{(z)}{\Theta}}_{\mathbf{p}_{1}\mathbf{p}_{2}}$ yields

[TABLE]

This allows to define the Z-dihedron strain measure, once the displacement be rescaled-back by $\textstyle\mathchar 280\relax$:

[TABLE]

For a generic Z-dihedral angle centered in $\mathbf{p}_{i}$, we get

[TABLE]

Analogous computations allows to determine the C-dihedron strain measure

[TABLE]

4 Membrane and bending energy

The above calculations show that $\textstyle\mathchar 290\relax$ and $\textstyle\mathchar 288\relax$ depend upon both the in-plane and the out-of-plane components of $\mathbf{u}^{\ell}$, while $\operatorname{\overset{(c)}{\Psi}}$, and $\operatorname{\overset{(z)}{\Psi}}$ depend upon the out-of-plane component of $\mathbf{u}^{\ell}$; moreover, (19) shows that the sum of all $\textstyle\mathchar 288\relax$ depends on the out-of-plane component of $\mathbf{u}^{\ell}$. We introduce the following splitting of the energy into membrane and bending parts:

[TABLE]

defined by

[TABLE]

where $\mathcal{U}_{\ell}^{(e)}$, $\mathcal{U}_{\ell}^{(w)}$, $\mathcal{U}_{\ell}^{(s)}$ and $\mathcal{U}_{\ell}^{(d)}$ are the edge, wedge, self-, and dihedral energy, respectively.

The self-stress term is the outcome of the fact that the angle at ease between consecutive edges is greater than \frac{2}{3}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 281\relax}}{\mbox{\boldmath\textstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 281\relax}}}: this means that in the flat reference configuration the graphene sheet is not stress-free, and {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 284\relax}}{\mbox{\boldmath\textstyle\mathchar 284\relax}}{\mbox{\boldmath\scriptstyle\mathchar 284\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 284\relax}}}_{0} represents the pre-stress couple (see [11]). The constants $k^{l}$, k^{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}, $k^{\mathcal{Z}}$, and $k^{\mathcal{C}}$ can be deduced by means of the 2nd-generation Brenner potential. In the wedge energy we may interpret the presence of $\ell^{2}$ as a scaling of the constant k^{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}, introduced to keep the energy finite as the lattice size $\ell$ goes to zero. We notice that the multiplication factor $\ell^{2}$ appears in all contributions to the membrane energy.

With the notation introduced in Section 3 we now write the energy more explicitly. The edge energy can be written as

[TABLE]

while the wedge energy reads:

[TABLE]

Analogously, the self-energy becomes

[TABLE]

We further split the dihedral energy $\mathcal{U}_{\ell}^{(d)}$ in

[TABLE]

where

[TABLE]

is the contribution of the Z-dihedra, and

[TABLE]

is the contribution of the C-dihedra.

We notice that the bending energy here obtained is the same as [5, 8]. In the following, we then focus on the membrane energy.

5 The continuum limit of the membrane energy

In the previous sections the discrete energy $\mathcal{U}_{\ell}$ was defined over the lattice $(L_{1}(\ell)\cup L_{2}(\ell))\cap\Omega$. By letting the lattice size $\ell$ go to zero the discrete set $(L_{1}(\ell)\cup L_{2}(\ell))\cap\Omega$ invades the domain $\Omega$, and the displacement functions $w^{\ell}$ and $\mathbf{v}^{\ell}$ will approach two functions $w$ and $\mathbf{v}$ defined over $\Omega$. To derive a continuous energy, defined over the domain $\Omega$, from the discrete energy $\mathcal{U}_{\ell}$ we need to specify the relations between $w^{\ell}$ and $w$ and between $\mathbf{v}^{\ell}$ and $\mathbf{v}$.

We assume $w:\Omega\to\mathbb{R}$ to be twice continuously differentiable and to, almost, coincide with ${w}^{\ell}$ over the lattice $(L_{1}(\ell)\cup L_{2}(\ell))\cap\Omega$. More precisely, we assume that

[TABLE]

If we think of $w$ and ${w}^{\ell}$ as the macroscopic and microscopic displacements, respectively, then (35) can be thought as a Cauchy–Born rule. The assumption (35) is motivated and essentially justified in [8].

For the in-plane displacement we assume $\mathbf{v}:\Omega\to\mathbb{R}^{2}$ to be continuously differentiable and

[TABLE]

Thus, over the lattice $L_{1}(\ell)$ we make the Cauchy–Born assumption for the in-plane displacement, while over the lattice $L_{2}(\ell)$ this assumption is relaxed by introducing a “shift displacement” $\mathbf{s}:\Omega\to\mathbb{R}^{2}$ that we assume to be continuously differentiable. Clearly, if the shift displacement is set equal to zero we have the Cauchy–Born rule over both lattices, but energetically it might be convenient to have a shift displacement $\mathbf{s}$ different from zero. The minus sign in front of the term containing $\mathbf{s}$ is introduced simply for later convenience. The assumption (36) is motivated and essentially justified in [4].

We note that for $\mathbf{x}^{\ell}\in L_{2}(\ell)$ we have that $\mathbf{x}^{\ell}+\ell\mathbf{p}_{i}\in L_{1}(\ell)$, for $i=1,2,3$, and

[TABLE]

while for $\mathbf{x}^{\ell}\in L_{1}(\ell)$ we have that $\mathbf{x}^{\ell}-\ell\mathbf{p}_{i}\in L_{1}(\ell)$, for $i=1,2,3$, and

[TABLE]

Thus,

[TABLE]

and similarly for $\mathbf{x}^{\ell}\in L_{2}(\ell)$ we have that

[TABLE]

For $i=1,2,3$ and for $\mathbf{x}^{\ell}\in L_{2}(\ell)$ the edge strain measure defined in (45) writes as

[TABLE]

where

[TABLE]

are the linearized and the von Kármán strain tensors, respectively. Similarly, the wedge strain measure defined in (18) rewrites as

[TABLE]

where we set

[TABLE]

and where we used the fact that $\mathbf{p}_{i+1}+\mathbf{p}_{i+2}=-\mathbf{p}_{i}$. A similar computation shows that, see (20),

[TABLE]

We now compute the energies. The edge energy (30) becomes

[TABLE]

where the second equality is obtained by noticing that the number of points $\mathbf{x}^{\ell}$ in $L_{2}(\ell)\cap\Omega$ is of order $1/\ell^{2}$. Let $|E^{\ell}(\mathbf{x}^{\ell})|=\ell^{2}3\sqrt{3}/2$ be the area of the hexagon $E^{\ell}(\mathbf{x}^{\ell})$ of side $\ell$ centred at $\mathbf{x}^{\ell}$, see Figure 6, and let {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 287\relax}}{\mbox{\boldmath\textstyle\mathchar 287\relax}}{\mbox{\boldmath\scriptstyle\mathchar 287\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 287\relax}}}_{E^{\ell}(\mathbf{x}^{\ell})}(\mathbf{x}) be the characteristic function of $E^{\ell}(\mathbf{x}^{\ell})$, i.e., the function equal to $1$ if $\mathbf{x}\in E^{\ell}(\mathbf{x}^{\ell})$ and [math] otherwise. The energy $\mathcal{U}_{\ell}^{(e)}$ may be rewritten as

[TABLE]

and since the function \sum_{\mathbf{x}^{\ell}\in L_{2}(\ell)\cap\Omega}\big{(}\mathbf{E}^{\rm vk}[\mathbf{v},w](\mathbf{x}^{\ell})\cdot\mathbf{p}_{i}\otimes\mathbf{p}_{i}+\mathbf{s}(\mathbf{x}^{\ell})\cdot\mathbf{p}_{i}\big{)}^{2}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 287\relax}}{\mbox{\boldmath\textstyle\mathchar 287\relax}}{\mbox{\boldmath\scriptstyle\mathchar 287\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 287\relax}}}_{E^{\ell}(\mathbf{x}^{\ell})}(\mathbf{x}) converges, as $\ell$ goes to zero, to \big{(}\mathbf{E}^{\rm vk}[\mathbf{v},w](\mathbf{x})\cdot\mathbf{p}_{i}\otimes\mathbf{p}_{i}+\mathbf{s}(\mathbf{x})\cdot\mathbf{p}_{i}\big{)}^{2} we have that

[TABLE]

By taking into account that the expression for \overset{(1)}{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 288\relax}}{\mbox{\boldmath\textstyle\mathchar 288\relax}}{\mbox{\boldmath\scriptstyle\mathchar 288\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 288\relax}}}}_{i} and \overset{(2)}{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 288\relax}}{\mbox{\boldmath\textstyle\mathchar 288\relax}}{\mbox{\boldmath\scriptstyle\mathchar 288\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 288\relax}}}}_{i} are identical, the wedge energy (31) rewrites:

[TABLE]

By introducing the triangles $T^{\ell}(\mathbf{x}^{\ell})$ centered at $\mathbf{x}^{\ell}$ of area $3\sqrt{3}\ell^{2}/4$ as depicted in Figure 6 and proceeding as above, we deduce that

[TABLE]

The limit of the membrane energy is

[TABLE]

By using the relations

[TABLE]

we find that

[TABLE]

where (\mathbf{E}^{\rm vk})_{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 267\relax}}{\mbox{\boldmath\textstyle\mathchar 267\relax}}{\mbox{\boldmath\scriptstyle\mathchar 267\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 267\relax}}}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 268\relax}}{\mbox{\boldmath\textstyle\mathchar 268\relax}}{\mbox{\boldmath\scriptstyle\mathchar 268\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 268\relax}}}} are the components of $\mathbf{E}^{\rm vk}=\mathbf{E}^{\rm vk}[\mathbf{v},w]$ with respect to the basis $\{\mathbf{e}_{i}\}$. A simple computation then shows that

[TABLE]

Similarly we find

[TABLE]

With these identities it follows that

[TABLE]

The shift displacement $\mathbf{s}$ enters into the membranal energy without any derivatives and therefore it can be minimized once and for all. We readily see that the shift displacement that minimizes the energy is

[TABLE]

and that

[TABLE]

By expanding the squares and reorganizing the terms we find:

[TABLE]

which is clearly an isotropic energy. By means of the relation \det\mathbf{A}=\big{(}(\operatorname{\textrm{tr}}\mathbf{A})^{2}-|\mathbf{A}|^{2})/2, which holds for every two by two matrix $\mathbf{A}$, we may write

[TABLE]

with

[TABLE]

6 Conclusions

We have proposed a non-linear continuum model for the mechanical behavior of graphene inferred from Molecular Dynamics potentials. Starting from a harmonic approximation of the energy as depicted by the 2nd-generation Brenner potential, we have found a discrete energy, depending on the displacement of each atom, which is sensitive to change of (i) the distance between two atoms (edge energy), (ii) the angle spanned by three subsequent atoms (wedge energy), (iii) two types of dihedral angles generated by the plane spanned by four subsequent atoms (C- and Z-dihedral energy). Thus, up-to-third neighbors interactions have been considered. Moreover, we have taken into account the presence of the self-stress, as predicted by the 2nd-generation Brenner potential (self-energy).

We have introduced a different scaling for in-plane and out-of-plane components, and this feature has produced a coupling in the strain measures at the discrete level. In particular, while the dihedral and the self-stress energies depend just on the out-of-plane components, the edge and the wedges energies depend on both. These two latter contributions determine the membrane energy, while the former two are part of the bending energy. With the scaling here adopted, the discrete bending energy turns out to be the same as that already considered in [5, 6, 8]; for this reason, we focused on the membrane energy.

The deduced discrete energy is defined over the two Bravais lattices generating the graphene sheet. By letting the size of the lattices to zero, they invade a continuum domain $\Omega$ and the discrete displacement functions approach two continuous functions $w$ and $\mathbf{v}$, representing the in-plane and the out-of-plane continuum displacements, defined over $\Omega$. To obtain the continuum energy, it has been necessary to specify the relation between the discrete displacements and the corresponding continuum functions. To this end, motivated by [4], we have made the Cauchy–Born assumption for the in-plane displacement defined over one of the Bravais lattices, while we have relaxed this assumption for the second lattice, by introducing a shift displacement $\mathbf{s}$.

With these assumptions, we have found a continuum membrane energy depending on $w$, $\mathbf{v}$ and $\mathbf{s}$. Since the shift displacement $\mathbf{s}$ enters into the membranal energy without any derivatives, we have minimized it once and for all; this has lead to a membrane energy depending just on $w$ and $\mathbf{v}$.

Thus, on considering the bending energy already deduced in [5, 6, 8] and the membrane energy here found, we can state that the total continuum energy of graphene reads:

[TABLE]

where the membrane energy is given by

[TABLE]

and the bending energy is

[TABLE]

We find that

[TABLE]

where $k^{l}$, k^{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}}, $k^{\mathcal{C}}$, $k^{\mathcal{Z}}$ and {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 284\relax}}{\mbox{\boldmath\textstyle\mathchar 284\relax}}{\mbox{\boldmath\scriptstyle\mathchar 284\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 284\relax}}}_{0} are the constants entering the edge, wedge, C-dihedral, Z-dihedral, and self-, atomistic energies.

Within the limits of this formal deduction, we have found that graphene can be modeled as a classical Föppl–von Kármán plate, where the constitutive constants depend on the atomistic interactions, as described by the 2nd-generation Brenner potential.

Acknowledgments

A.F. acknowledges support from Sapienza University of Rome through the projects RP116154C92AF8A4 “Multiscale Mechanics of 2D Materials: Modeling and Applications’ and RM11715C7F61C3E8 “Shape morphing. From advanced differential geometry to applications in engineering and architecture”.

R.P.  acknowledges support from the Università di Pisa through the project PRA_2018_61 “Modellazione multi-scala in ingegneria strutturale”.