A nonlinear hyperelasticity model for single layer blue phosphorus based on ab-initio calculations

TL;DR

This paper introduces a fully nonlinear hyperelasticity model for single layer blue phosphorus, calibrated from ab-initio DFT calculations, capturing anisotropic large-strain behavior and validated through simulations.

Contribution

The paper develops the first fully nonlinear hyperelastic model for blue phosphorus based on ab-initio data, enabling accurate large-strain simulations.

Findings

Model shows good agreement with DFT results.

Can be integrated into shell finite element methods.

Elasticity constants at small deformations are provided.

Abstract

A new hyperelastic membrane material model is proposed for single layer blue phosphorus ($\beta\text{-P}$), also known as blue phosphorene. The model is fully nonlinear and captures the anisotropy of $\beta\text{-P}$ at large strains. The material model is calibrated from density functional theory (DFT) calculations considering a set of elementary deformation states. Those are pure dilatation and uniaxial stretching along the armchair and zigzag directions. The material model is compared and validated with additional DFT results and existing DFT results from the literature, and the comparison shows good agreement. The new material model can be directly used within computational shell formulations that are for example based on rotation-free isogeometric finite elements. This is demonstrated by simulations of the indentation and vibration of single layer blue phosphorus sheets at micrometer scales. The elasticity constants at small deformations are also reported.

Tables

Table 1. Table 1 : Membrane material constants: Pure dilatation part.

	$n_2$	$n_3$	$n_4$	$n_5$
[N/m]	$24.49$	$-66.07$	$276.19$	$-444$

Table 2. Table 2 : Membrane material constants: Isotropic shear part.

	${\mu }_{10}$	${\mu }_{12}$	${\mu }_{13}$	${\mu }_{14}$	${\mu }_{20}$	${\mu }_{22}$	${\mu }_{23}$	${\mu }_{24}$
[N/m]	$61.88$	$-346.44$	$-670.88$	$559.77$	$-1029.84$	$487.59$	$-1076.74$	$-708.98$

Table 3. Table 3 : Membrane material constants: Anisotropic shear part.

	${\eta }_0$	${\eta }_2$	${\eta }_3$	${\eta }_4$
[N/m]	$62.22$	$-838.95$	$-161.41$	$-427.24$

Table 4. Table 4 : Constants for the thickness variation.

$p_1$	$p_2$	$s_2$	$s_3$	$s_3$
$-0.22$	$-27.48$	$-1245$	$119$	$-3.77$

Table 5. Table 5 : Linear elasticity constants: AC = armchair and ZZ = zigzag.

	$E$ [N/m]	$G$ [N/m]	$K$ [N/m]	$\nu$
Current work	$75.85$ (in AC and ZZ)	$30.94$	$48.98$	$0.226$
Sun et al. [14]	-	-	42.66	-
Peng et al. [27]	69.36 (in AC) and 66.24 (in ZZ)	-	-	-

Table 6. Table 1 : Lattice structure of blue phosphorus.

Lattice parameter (Å)	Thickness (Å)	Bond length (Å)
$3.28$	$1.23$	$2.26$

Equations

$$\begin{array}[]{lll}\boldsymbol{X}\!\!\!&=&\!\!\!\boldsymbol{X}(\xi^{\alpha})~{},\end{array}$$

$$\begin{array}[]{lll}\boldsymbol{x}\!\!\!&=&\!\!\!\boldsymbol{x}(\xi^{\alpha},t)~{},\end{array}$$

$$\begin{array}[]{lll}\boldsymbol{A}_{\alpha}=\boldsymbol{X}_{\!,\alpha}~{},\end{array}$$

$$\begin{array}[]{lll}\boldsymbol{a}_{\alpha}=\boldsymbol{x}_{,\alpha}~{},\end{array}$$

$$\begin{array}[]{lll}\boldsymbol{A}_{\alpha}\cdot\boldsymbol{A}^{\beta}\!\!\!&=&\!\!\!\boldsymbol{a}_{\alpha}\cdot\boldsymbol{a}^{\beta}=\delta_{\alpha}^{\beta}~{},\end{array}$$

$$\begin{array}[]{lll}A_{\alpha\beta}\!\!\!&:=&\!\!\!\boldsymbol{A}_{\alpha}\cdot\boldsymbol{A}_{\beta}~{},\end{array}$$

$$\begin{array}[]{lll}A^{\alpha\beta}\!\!\!&:=&\!\!\!\boldsymbol{A}^{\alpha}\cdot\boldsymbol{A}^{\beta}~{}.\end{array}$$

$$\begin{array}[]{lll}\boldsymbol{N}\!\!\!&=&\!\!\!\boldsymbol{A}_{1}\times\boldsymbol{A}_{2}/\|\boldsymbol{A}_{1}\times\boldsymbol{A}_{2}\|~{},\end{array}$$

$$\begin{array}[]{lll}\boldsymbol{n}\!\!\!&=&\!\!\!\boldsymbol{a}_{1}\times\boldsymbol{a}_{2}/\|\boldsymbol{a}_{1}\times\boldsymbol{a}_{2}\|~{}.\end{array}$$

$$\begin{array}[]{lll}\displaystyle\mathrm{d}\boldsymbol{x}\!\!\!&=&\!\!\!\displaystyle\boldsymbol{F}\,\mathrm{d}\boldsymbol{X}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle\boldsymbol{F}\!\!\!&=&\!\!\!\displaystyle\displaystyle\frac{\partial{\boldsymbol{x}}}{\partial{\boldsymbol{X}}}~{}.\end{array}$$

$$\begin{array}[]{lll}\boldsymbol{F}\!\!\!&=&\!\!\!\boldsymbol{R}\,\boldsymbol{U}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle\boldsymbol{U}\!\!\!&=&\!\!\!\displaystyle\sum_{\alpha=1,2}{\lambda_{\alpha}\,\boldsymbol{Y}_{\!\!\alpha}\otimes\boldsymbol{Y}_{\!\!\alpha}}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle\boldsymbol{C}\!\!\!&=&\!\!\!\displaystyle\boldsymbol{F}^{\text{T}}\,\boldsymbol{F}=\boldsymbol{U}^{2}=\sum_{\alpha=1,2}{\Lambda_{\alpha}\,\boldsymbol{Y}_{\!\!\alpha}\otimes\boldsymbol{Y}_{\!\!\alpha}}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle\boldsymbol{E}^{(0)}\!\!\!&:=&\!\!\!\displaystyle\ln\,\boldsymbol{U}=\frac{1}{2}\ln\,\boldsymbol{C}=\sum_{\alpha=1,2}{\ln(\lambda_{\alpha})\,\boldsymbol{Y}_{\!\!\alpha}\otimes\boldsymbol{Y}_{\!\!\alpha}}~{}.\end{array}$$

$$\begin{array}[]{lll}\displaystyle\boldsymbol{E}^{(0)}\!\!\!&=&\!\!\!\displaystyle\boldsymbol{E}^{(0)}_{\text{area}}+\boldsymbol{E}^{(0)}_{\text{dev}}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle\boldsymbol{E}^{(0)}_{\text{area}}\!\!\!&=&\!\!\!\displaystyle\frac{1}{2}\mathrm{tr}\,\left(\boldsymbol{E}^{(0)}\right)\boldsymbol{I}\end{array}$$

$$\begin{array}[]{lll}\displaystyle\boldsymbol{E}^{(0)}_{\text{dev}}\!\!\!&=&\!\!\!\displaystyle\boldsymbol{E}^{(0)}-\boldsymbol{E}^{(0)}_{\text{area}}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle\mathbb{H}_{n}\!\!\!&:=&\!\!\!\Re\left[\left(\widehat{\boldsymbol{x}}+i\widehat{\boldsymbol{y}}\right)^{(n)}\right]=\displaystyle\left\{{\begin{array}[]{*{20}{ll}}\Re\left[\left(\widehat{\boldsymbol{M}}+i\widehat{\boldsymbol{N}}\right)^{(m)}\right];&n=2m~{},\\[8.53581pt] \Re\left[\left(\widehat{\boldsymbol{x}}+i\widehat{\boldsymbol{y}}\right)\otimes\left(\widehat{\boldsymbol{M}}+i\widehat{\boldsymbol{N}}\right)^{(m)}\right];&n=2m+1~{},\end{array}}\right.\end{array}$$

$$\begin{array}[]{lll}\widehat{\boldsymbol{M}}\!\!\!&:=&\!\!\!\widehat{\boldsymbol{x}}\,\otimes\,\widehat{\boldsymbol{x}}-\widehat{\boldsymbol{y}}\,\otimes\,\widehat{\boldsymbol{y}}~{},\\[8.53581pt] \widehat{\boldsymbol{N}}\!\!\!&:=&\!\!\!\widehat{\boldsymbol{x}}\,\otimes\,\widehat{\boldsymbol{y}}+\widehat{\boldsymbol{y}}\,\otimes\,\widehat{\boldsymbol{x}}~{}.\end{array}$$

$$\begin{array}[]{lll}\displaystyle\mathcal{J}_{1}\!\!\!&=&\!\!\!\displaystyle\epsilon_{\text{a}}=\ln J~{},\\[8.53581pt] \displaystyle\mathcal{J}_{2}\!\!\!&=&\!\!\!\displaystyle\frac{1}{2}\boldsymbol{E}^{(0)}_{\text{dev}}:\boldsymbol{E}^{(0)}_{\text{dev}}=(\ln\lambda)^{2}~{},\\[8.53581pt] \displaystyle\mathcal{J}_{3}\!\!\!&=&\!\!\!\displaystyle(\ln\lambda)^{3}\,\cos(6\theta)~{},\end{array}$$

$$\begin{array}[]{lll}\mbox{\boldmath$\sigma$}=\displaystyle\frac{1}{\det\boldsymbol{F}}\,\boldsymbol{F}\,\boldsymbol{S}\,\boldsymbol{F}^{\text{T}}~{},\end{array}$$

$$\begin{array}[]{lll}\boldsymbol{S}=\displaystyle 2\frac{\partial{W}}{\partial{\boldsymbol{C}}}{}\end{array}$$

$$\begin{array}[]{lll}\displaystyle\chi\!\!\!&=&\!\!\!\displaystyle\sum_{I=1,N_{\text{QM}}}{\|W^{I}-W_{\text{QM}}^{I}\|_{2}+\|\mbox{\boldmath$\sigma$}^{I}-\mbox{\boldmath$\sigma$}_{\text{QM}}^{I}\|_{2}}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle W\!\!\!&=&\!\!\!\displaystyle f_{1}(\mathcal{J}_{1})+f_{2}(\mathcal{J}_{1})\,\mathcal{J}_{2}+f_{3}(\mathcal{J}_{1})\,\mathcal{J}_{2}^{2}+f_{4}(\mathcal{J}_{1})\,\mathcal{J}_{3}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle f_{1}\!\!\!&=&\!\!\!\displaystyle n_{2}\,\mathcal{J}_{1}^{2}+n_{3}\,\mathcal{J}_{1}^{3}+n_{4}\,\mathcal{J}_{1}^{4}+n_{5}\mathcal{J}_{1}^{5}~{},\\[8.53581pt] \displaystyle f_{2}\!\!\!&=&\!\!\!\displaystyle\mu_{10}+\mu_{12}\,\mathcal{J}_{1}^{2}+\mu_{13}\,\mathcal{J}_{1}^{3}+\mu_{14}\,\mathcal{J}_{1}^{4}~{},\\[8.53581pt] \displaystyle f_{3}\!\!\!&=&\!\!\!\displaystyle\mu_{20}+\mu_{22}\,\mathcal{J}_{1}^{2}+\mu_{23}\,\mathcal{J}_{1}^{3}+\mu_{24}\,\mathcal{J}_{1}^{4}~{},\\[8.53581pt] \displaystyle f_{4}\!\!\!&=&\!\!\!\displaystyle\eta_{0}+\eta_{2}\,\mathcal{J}_{1}^{2}+\eta_{3}\,\mathcal{J}_{1}^{3}+\eta_{4}\,\mathcal{J}_{1}^{4}~{}.\end{array}$$

$$\begin{array}[]{lll}\displaystyle\boldsymbol{S}^{(0)}\!\!\!&:=&\!\!\!\displaystyle\frac{\partial{W}}{\partial{\boldsymbol{E}^{(0)}}}=\left[f_{1}^{\prime}+f_{2}^{\prime}\,\mathcal{J}_{2}+f_{3}^{\prime}\,\mathcal{J}_{2}^{2}+f_{4}^{\prime}\,\mathcal{J}_{3}\right]\,\boldsymbol{I}+f_{2}\,\boldsymbol{E}^{(0)}_{\text{dev}}+2f_{3}\,\mathcal{J}_{2}\,\boldsymbol{E}^{(0)}_{\text{dev}}+f_{4}\,\frac{\partial{\mathcal{J}_{3}}}{\partial{\boldsymbol{E}^{(0)}}}~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle f_{1}^{\prime}\!\!\!&=&\!\!\!\displaystyle 2n_{2}\,\mathcal{J}_{1}+3n_{3}\,\mathcal{J}_{1}^{2}+4n_{4}\,\mathcal{J}_{1}^{3}+5n_{5}\,\mathcal{J}_{1}^{4}~{},\\[8.53581pt] \displaystyle f_{2}^{\prime}\!\!\!&=&\!\!\!\displaystyle 2\mu_{12}\,\mathcal{J}_{1}+3\mu_{13}\,\mathcal{J}_{1}^{2}+4\mu_{14}\,\mathcal{J}_{1}^{3}~{},\\[8.53581pt] \displaystyle f_{3}^{\prime}\!\!\!&=&\!\!\!\displaystyle 2\mu_{22}\,\mathcal{J}_{1}+3\mu_{23}\,\mathcal{J}_{1}^{2}+4\mu_{24}\,\mathcal{J}_{1}^{3}~{},\\[8.53581pt] \displaystyle f_{4}^{\prime}\!\!\!&=&\!\!\!\displaystyle 2\eta_{2}\,\mathcal{J}_{1}+3\eta_{3}\,\mathcal{J}_{1}^{2}+4\eta_{4}\,\mathcal{J}_{1}^{3}~{}.\end{array}$$

$$\begin{array}[]{lll}\displaystyle\mathbb{C}^{(0)}:=\displaystyle\frac{\partial^{2}{W}}{\partial{\boldsymbol{E}^{(0)}}\,\partial{\boldsymbol{E}^{(0)}}}\!\!\!&=&\!\!\!\displaystyle\left[f_{1}^{\prime\prime}+f_{2}^{\prime\prime}\,\mathcal{J}_{2}+f_{3}^{\prime\prime}\,\mathcal{J}_{2}^{2}+f_{4}^{\prime\prime}\,\mathcal{J}_{3}\right]\boldsymbol{I}\otimes\boldsymbol{I}+2f_{3}\,\boldsymbol{E}^{(0)}_{\text{dev}}\otimes\boldsymbol{E}^{(0)}_{\text{dev}}\\[8.53581pt] \!\!\!&+&\!\!\!\displaystyle(f_{2}^{\prime}+2f_{3}^{\prime}\mathcal{J}_{2})\,\left[\boldsymbol{I}\otimes\boldsymbol{E}^{(0)}_{\text{dev}}+\boldsymbol{E}^{(0)}_{\text{dev}}\otimes\boldsymbol{I}\right]+f_{4}\,\frac{\partial^{2}{\mathcal{J}_{3}}}{\partial{\boldsymbol{E}^{(0)}}\,\partial{\boldsymbol{E}^{(0)}}}\\[8.53581pt] \!\!\!&+&\!\!\!\displaystyle f_{4}^{\prime}\,\left[\boldsymbol{I}\otimes\frac{\partial{\mathcal{J}_{3}}}{\partial{\boldsymbol{E}^{(0)}}}+\frac{\partial{\mathcal{J}_{3}}}{\partial{\boldsymbol{E}^{(0)}}}\otimes\boldsymbol{I}\right]+(f_{2}+2f_{3}\,\mathcal{J}_{2})\left[\mathcal{I}-\frac{1}{2}\boldsymbol{I}\otimes\boldsymbol{I}\right]~{},\end{array}$$

$$\begin{array}[]{lll}\displaystyle f_{1}^{\prime\prime}\!\!\!&=&\!\!\!\displaystyle 2n_{2}+6n_{3}\,\mathcal{J}_{1}+12n_{4}\,\mathcal{J}_{1}^{2}+20n_{5}\,\mathcal{J}_{1}^{3}~{},\\[8.53581pt] \displaystyle f_{2}^{\prime\prime}\!\!\!&=&\!\!\!\displaystyle 2\mu_{12}+6\mu_{13}\,\mathcal{J}_{1}+12\mu_{14}\,\mathcal{J}_{1}^{2}~{},\\[8.53581pt] \displaystyle f_{3}^{\prime\prime}\!\!\!&=&\!\!\!\displaystyle 2\mu_{22}+6\mu_{23}\,\mathcal{J}_{1}+12\mu_{24}\,\mathcal{J}_{1}^{2}~{},\\[8.53581pt] \displaystyle f_{4}^{\prime\prime}\!\!\!&=&\!\!\!\displaystyle 2\eta_{2}+6\eta_{3}\,\mathcal{J}_{1}+12\eta_{4}\,\mathcal{J}_{1}^{2}~{}.\end{array}$$

Full text

**A nonlinear hyperelasticity model for single layer blue phosphorus based on ab-initio calculations

**

Reza Ghaffari§111Email: ghaffari@aices.rwth-aachen.de, Farzad Shirazian§222Email: shirazian@aices.rwth-aachen.de, Ming Hu*†*333Email: hu@sc.edu and Roger A. Sauer§444Corresponding author, email: sauer@aices.rwth-aachen.de

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

RWTH Aachen University, Templergraben 55, 52056 Aachen, Germany

†Department of Mechanical Engineering, University of South Carolina,

541 Main Street, Columbia, SC 29208, USA *

Abstract: A new hyperelastic membrane material model is proposed for single layer blue phosphorus ($\beta$-P), also known as blue phosphorene. The model is fully nonlinear and captures the anisotropy of $\beta\text{-P}$ at large strains. The material model is calibrated from density functional theory (DFT) calculations considering a set of elementary deformation states. Those are pure dilatation and uniaxial stretching along the armchair and zigzag directions. The material model is compared and validated with additional DFT results and existing DFT results from the literature, and the comparison shows good agreement. The new material model can be directly used within computational shell formulations that are for example based on rotation-free isogeometric finite elements. This is demonstrated by simulations of the indentation and vibration of single layer blue phosphorus sheets at micrometer scales. The elasticity constants at small deformations are also reported.

Keywords: Anisotropic hyperelasticity; curvilinear membrane formulation; density functional theory; nonlinear finite element methods; single layer blue phosphorus (phosphorene).

1 Introduction

Two-dimensional materials are a fascinating group of materials that can have very different mechanical, electronic, chemical, and optical properties than their bulk form [1]. Since the discovery of graphene, scientists have looked for other interesting 2D materials, and phosphorus allotropes are very promising candidates due to their unique properties. Up to now, many stable two dimensional structures have been predicted for phosphorus [2, 3]. Single layer black phosphorus (also know as phosphorene), as the first of this family to be discovered, has a large range of applications in field effect transistors [4], optoelectronic devices [5], sensors [6], batteries [7], and energy storage [8, 9]. Zhu and Tománek [10] predict the structure of single layer blue phosphorus. They use ab-initio calculations to investigate its structure and predict a higher band gap for single layer blue phosphorus than for black phosphorus. For simplicity, black and blue phosphorus refer to two dimensional structures of phosphorus in the rest of this paper. Zhang et al. [11] grow a single layer of blue phosphorus through epitaxial growth on Au(111) by using black phosphorous as precursor. While black phosphorus has a band gap of 0.3-0.4 eV, blue phosphorus has a much higher band gap of 2-3 eV [10, 12], indicating its potential application in field effect transistors and optoelectronic devices.

Mogulkoc et al. [13] study the electronic and optical properties of phosphide/blue phosphorus heterostructures and suggest their applicability in new generation optoelectric devices duo to their transparency to visible light and great absorption over the UV range. Sun et al. [14] use density functional theory (DFT) to study thermomechanical properties of black and blue phosphorus. Li et al. [4] study the application of single-layer and double-layer black and blue phosphorus in Li-Ion batteries and predict both of them can be good electrode materials due to their small diffusion energy barriers and ability to maintain their layered structures during lithiation and delithiation processes. Liu et al. [15] investigate the effect of external strain on the electronic properties of blue phosphorus. Xiao et al. [12] predict blue phosphorus to be an indirect p-type semiconductor with anisotropic properties. They calculate the carrier mobility of a monolayer blue phosphorus under uniaxial and biaxial strains and show strain engineering can tune the properties of blue phosphorus.

Cooper et al. [16] use Taylor expansion of the strain energy in order to describe molybdenum disulfide. Setoodeh and Farahmand [17] propose an anisotropic continuum model for black phosphorus. They give two parameter sets for uniaxial stretch along the armchair and zigzag directions. This model does not satisfy periodicity of the black phosphorus lattice and does not predict the behaviour of the structure for loading along an arbitrary direction. So it cannot be used in continuum models such as the membrane model of Sauer et al. [18] and the shell model of Duong et al. [19]. Kumar and Parks [20] propose an anisotropic, hyperelastic membrane material model for graphene. This model is nonlinear and directly calibrated from DFT data. It therefore does not require an interatomic interaction potential, and thus avoids any inaccuracies resulting from it. Kumar and Parks [20] use their model to simulate nano indentation of a micro scale graphene sheet. Ghaffari et al. [21] extend their model to a Kirchhoff-Love shell model and implement it within the finite element formulation of Duong et al. [19]. Indentation and peeling of graphene sheets, and torsion and bending of carbon nanotubes (CNTs) are simulated with this shell model, and it is used for the nonlinear modal analysis of graphene sheets and CNTs by Ghaffari and Sauer [22]. Ghaffari and Sauer [23] propose a new efficient computational shell model for graphene and apply it to carbon nanocones (CNCs). Shirazian et al. [24] propose a new set of material constants for the membrane material model of Kumar and Parks [20], Ghaffari et al. [21] and Ghaffari and Sauer [23].

In the current work, a new hyperelastic membrane material model is proposed for blue phosphorus. It is based on a set of invariants that are obtained from the symmetry of the lattice. This material model is nonlinear, anisotropic, and calibrated with DFT data. In summary, the novelties of the current work are:

A new hyperelastic continuum membrane material model is proposed for blue phosphorus that is directly based on ab-initio data. It thus circumvents the use of atomistic potentials.

• 2.

It is fully nonlinear and can capture anisotropic behaviour of the material.

• 3.

Due to the inherent efficiency of continuum formulations, the new model can be used to simulate systems a large length scales.

• 4.

This is demonstrated by simulating the indentation and vibration of phosphorene at the micrometer scales.

• 5.

The new model admits extension to other 2D materials, and it can be extended to finite temperatures based on the new formulation of Ghaffari and Sauer [25].

The remainder of this paper is organised as follows: In Sec. 2 the kinematics of deforming surfaces is summarised. In Sec. 3, a suitable structural tensor and a set of invariants are introduced for blue phosphorus. Based on these, a new hyperelastic membrane material model is proposed in Sec. 4. In Sec. 5, the model is calibrated and verified with the DFT results considering various test cases. Sec. 6 presents numerical indentation and modal analysis examples. The paper is concluded in Sec. 7.

2 Kinematics

In this section, the kinematics of deforming surface is summarised.

It will be used in the next sections to propose an anisotropic, hyperelastic continuum material model for blue phosphorus.

2.1 Surface description

In order to describe surfaces in very general terms, curvilinear coordinates are used. Accordingly, a point on the surface in the reference configuration $\mathcal{S}_{0}$ is indicated by

[TABLE]

and in the current configuration $\mathcal{S}$ by

[TABLE]

see Fig. 1. Here $\xi^{\alpha}$, for $\alpha=1,2$, are the curvilinear coordinates and $t$ is time. The tangent vectors of $\mathcal{S}_{0}$ then follow as

[TABLE]

and the tangent vectors of $\mathcal{S}$ as

[TABLE]

where ${\bullet,\alpha}=\partial\bullet/\partial\xi^{\alpha}$. The dual vectors $\boldsymbol{A}^{\alpha}$ and $\boldsymbol{a}^{\alpha}$ are defined such that

[TABLE]

where $\delta_{\alpha}^{\beta}$ is Kronecker delta defined by $[\delta_{\alpha}^{\beta}]=[1~{}0;0~{}1]$. The covariant surface metric of $\mathcal{S}_{0}$ is defined by

[TABLE]

while the contra-variant surface metric is

[TABLE]

The covariant, $a_{\alpha\beta}$, and contra-variant, $a^{\alpha\beta}$, metric of $\mathcal{S}$ are defined analogously. The unit normal vector of $\mathcal{S}_{0}$ is

[TABLE]

while the unit normal vector of $\mathcal{S}$ is

[TABLE]

2.2 Kinematics of deformation

The material model will be developed based on the logarithmic surface strain. The logarithmic surface strain can be either defined from the surface stretch tensor or the right Cauchy-Green surface deformation tensor, which both follow from the surface deformation gradient. Differential line elements in the reference and current configuration, denoted $\mathrm{d}\boldsymbol{X}$ and $\mathrm{d}\boldsymbol{x}$, are connected by the rank-two surface deformation gradient $\boldsymbol{F}=\boldsymbol{a}_{\alpha}\otimes\boldsymbol{A}^{\alpha}$ as

[TABLE]

where $\boldsymbol{F}$ is

[TABLE]

The polar decomposition of $\boldsymbol{F}$ can be written as

[TABLE]

where $\boldsymbol{R}$ and $\boldsymbol{U}$ are the surface rotation and right surface stretch tensor. $\boldsymbol{R}$ is proper orthogonal, i.e. $\boldsymbol{R}^{-1}=\boldsymbol{R}^{\text{T}}$ and $\det\boldsymbol{R}=1$, while $\boldsymbol{U}$ is symmetric. The spectral decomposition of $\boldsymbol{U}$ can be written as

[TABLE]

where $\lambda_{\alpha}$ and $\boldsymbol{Y}_{\!\!\alpha}$ are the eigenvalues and eigenvectors of $\boldsymbol{U}$. The right Cauchy-Green surface deformation tensor is

[TABLE]

where $\Lambda_{\alpha}=\lambda^{2}_{\alpha}$, and the eigenvectors of $\boldsymbol{C}$ and $\boldsymbol{U}$ are the same. $\boldsymbol{U}$ or $\boldsymbol{C}$ can be used to define the logarithmic surface strain $\boldsymbol{E}^{(0)}$ as

[TABLE]

$\boldsymbol{E}^{(0)}$ can be additively decomposed into the area-changing part $\boldsymbol{E}^{(0)}_{\text{area}}$ and area-preserving part $\boldsymbol{E}^{(0)}_{\text{dev}}$ as

[TABLE]

with

[TABLE]

and

[TABLE]

where $\boldsymbol{I}$ is the surface identity tensor on $\mathcal{S}_{0}$.

3 Structural tensor and invariants of blue phosphorus

In this section, the structural tensor and a set of invariants for hexagonal structures such as blue phosphorus are given. They are needed in order to model the anisotropic behaviour of the material. The structural tensors of a lattice with a symmetry group of $n$-fold rotational symmetry and reflection plane $C_{nv}$ can be written as [26]

[TABLE]

where $i$ is the unit imaginary number, $(\bullet)^{(n)}=(\bullet)\otimes(\bullet)...(\bullet)$ is tensor product of taking ($\bullet$) $n$ times, $\Re$ indicates the real part, $\widehat{\boldsymbol{x}}$ and $\widehat{\boldsymbol{y}}$ are two orthonormal vectors, where at least one of them is in the symmetry plane of the crystal (see Fig. 2), and $\widehat{\boldsymbol{M}}$ and $\widehat{\boldsymbol{N}}$ are defined as

[TABLE]

For blue phosphorus, $\widehat{\boldsymbol{x}}$ is assumed to be in the armchair direction. The blue phosphorus lattice has a rotational symmetry of $2\pi/3$. But, as will be shown later, its constitutive law has the higher symmetry of $\pi/3$, i.e. $n=6$. In this case, the material invariants based on the logarithmic surface strain can be obtained as [20, 23]

[TABLE]

with $\lambda=\sqrt{\lambda_{1}/\lambda_{2}}$, where $\lambda_{1}$ is the larger eigenvalue ($\lambda_{1}\geqslant\lambda_{2}$), and $\theta=\arccos(\widehat{\boldsymbol{x}}\cdot\boldsymbol{Y}_{\!\!1})$, where $\boldsymbol{Y}_{\!\!1}$ is the eigenvector corresponding to the largest eigenvalue. $\theta$ is the maximum stretch angle relative to the armchair direction. $\mathcal{J}_{1}$ and $\mathcal{J}_{2}$ capture isotropic behaviour and $\mathcal{J}_{3}$ captures anisotropic behaviour. $\mathcal{J}_{1}$ is characterising area-changing deformations, while $\mathcal{J}_{2}$ is characterising shape-changing deformations.

4 Material model

In this section, the previously obtained invariants are used to propose a functional for the strain energy density $W$. From $W$, the Cauchy surface555with units of [N/m] and in plane-stress format. stress tensor follows as

[TABLE]

where

[TABLE]

is the second Piola-Kirchhoff surface stress tensor. In general, $W$ contains several material constants that need to be calibrated from experimental, atomistic or quantum data. In the current work, $W$ is calibrated from DFT results of the strain energy density and Cauchy surface stress tensor $\sigma$. Details on the DFT simulations are given in A. A pure dilatation test, and two uniaxial stretching tests (one along the armchair and one along the zigzag direction) are used for the calibration. The material parameters are obtained such that the cost function

[TABLE]

is minimised. Here, “QM” indicates the quantum data, $N_{\text{QM}}$ is the number of the data points and $\|\bullet\|_{2}$ is the L2-norm. In this work, the strain energy functional, per unit reference area, is proposed as

[TABLE]

where $f_{i}$ are polynomial functions of $\mathcal{J}_{1}$ defined as

[TABLE]

Here $n_{i}$, $\mu_{ij}$, $\eta_{i}$ are material parameters. The three test cases used for their determination are presented in Sec. 5. The result of this calibration step are the material parameters given in Tab. 1, 2 and 3.

Using the partial derivatives of the three invariants $\mathcal{J}_{i}$ with respect to $\boldsymbol{E}^{(0)}$, see B, the logarithmic surface stress tensor, $\boldsymbol{S}^{(0)}$, that is conjugate to $\boldsymbol{E}^{(0)}$, can be computed as

[TABLE]

where $f_{i}^{\prime}:=\partial f_{i}/\partial\mathcal{J}_{1}$ are

[TABLE]

Using the second partial derivatives of $\mathcal{J}_{i}$ with respect to $\boldsymbol{E}^{(0)}$, see B, the corresponding elasticity tensor can be written as

[TABLE]

where $f_{i}^{\prime\prime}$ are

[TABLE]

For a finite element implementation, e.g. within the model of Sauer et al. [18], the contra-variant components of the Kirchhoff surface stress tensor $\tau^{\alpha\beta}$ and its corresponding elasticity tensor are needed. Therefore, the logarithmic surface stress and elasticity tensor need be transformed to the Cauchy surface stress tensor $\sigma$. $\tau^{\alpha\beta}$ and the contra-variant components of the Cauchy surface stress tensor, $\sigma^{\alpha\beta}$, can be connected as

[TABLE]

$\tau^{\alpha\beta}$ can be calculated from the second Piola-Kirchhoff surface stress tensor $\boldsymbol{S}$ as

[TABLE]

where $\boldsymbol{S}$ can be obtained from $W$ using the chain-rule as

[TABLE]

where $\partial W/\partial\boldsymbol{E}^{(0)}$ is given in (27) and $\partial\boldsymbol{E}^{(0)}/\partial\boldsymbol{C}$ can be found in Kumar and Parks [20] and Ghaffari et al. [21]. The transformation of (33) and the corresponding transformation for the elasticity tensor are provided by Kumar and Parks [20] and Ghaffari et al. [21] for graphene. Ghaffari and Sauer [23] pointed out that this transformation is computationally expensive and rather complicated. Therefore a finite difference formulation is considered here as discussed in C.

Next, the thickness variation $\lambda_{3}$ is discussed. $\lambda_{3}$ is computed from the DFT simulations which are used to propose the following isotropic function based on the first and second invariants

[TABLE]

where $p_{i}$ and $s_{i}$ are given in Tab. 4. A cost function similar to (24) is used to determine these constants from the three test cases presented in the following section.

5 Model calibration and validation

In this section, first the three calibration tests are described.. Then, the continuum results from the proposed blue phosphorus material model are compared and verified with DFT results of the current work and the literature. Finally the linear elastic constants are reported and compared with the literature.

5.1 Calibration

The calibration is conducted by performing one pure dilatation test and two uniaxial stretching tests. Altogether, 520 DFT data sets are used for this calibration, leading to the material parameters listed in Tab. 1, 2 and 3.

For pure dilatation, the hexagonal structure does not change and the material behaviour is isotropic. Fig. 3 shows the variation of the energy density and surface tension666\gamma=1/2\,\mathrm{tr}\,\mbox{\boldmath\sigma}. with respect to the area change invariant $\mathcal{J}_{1}$.

The energy density increases monotonically. The surface tension reaches a maximum and then begins to decrease. The DFT results are only considered up to this point since the lattice collapses beyond that point. This maximum occurs sooner in the continuum model than in the DFT results.

Next the lattice is stretched in the armchair or zigzag direction, while being fixed in the perpendicular direction in order to produce uniaxial stretch. The energy density and Cauchy surface stress components are shown in Fig. 4.

The material behaves anisotropic under uniaxial stretch and it fails sooner if stretched in the armchair direction rather than the zigzag direction. The strain energy density and stress $\sigma_{11}$ for stretch along the armchair direction are higher than for stretch along the zigzag direction, while it is vice versa for $\sigma_{22}$. The continuum and DFT results are in good agreement for all tests. The results are reported up to the maximum stress point beyond which the material becomes unstable.

Fig. 5 shows the calibration of the out-of-plane stretch $\lambda_{3}$ under pure dilatation and uniaxial stretch. The same 520 DFT data sets are used for this calibration, leading to the parameters listed in Tab. 4.

The material shows an approximate isotropic response for both cases. In both cases $\lambda_{3}$ decreases monotonically up to the failure point.

5.2 Validation

For uniaxial stretch along the armchair and zigzag direction, the current continuum model and DFT results are compared and validated with the DFT results of Peng et al. [27] in Fig. 4. As seen they are a good agreement.

For further validation of the model first, the Young’s modulus $E$, shear modulus $G$, bulk modulus $K$ and Poisson’s ratio $\nu$ are obtained in the small deformation regime and compared with reported values from the literature in Tab. 5. Given the plane-stress material model in Sec. 4, those material constant turn out to be $E=4G\,K/(G+K)$, $G=\mu_{10}/2$, $K=2n_{2}$ and $\nu=(K-G)/(K+G)$. Compared to other results, the proposed continuum model behaves stiffer by about 10-15%.

Next, the material is stretched in different directions and the variation of the stresses with respect to the stretch direction is reported in Fig. 6. The continuum and DFT results are in very good agreement. The DFT results for $\theta=10,20,40,50^{\circ}$ have not been used for the material calibration and so they can be used for the validation of the model. The material has a rotational periodicity of $\pi/3$ since a rotation of $\pi/3$ results in a mirror transformation of the structure with respect to its plane, which does not affect the mechanical properties of the blue phosphorus monolayer.

6 Simulation results at large length scales

The model has been validated in the previous section. In this section, it will be used to model a large micro-meter phosphorene specimen under indentation and vibration. The specimen consists of more than 5 million atoms. Molecular dynamics simulations would be too expensive for this length scale. The continuum model, on the other hand, has only 9000 nodes for the finest mesh. Thus the continuum model has more than 500 times less degrees of freedom than the atomistic system. Beyond that, there presently exists no suitable atomistic potential for blue phosphorus.

6.1 Phosphorene indentation

A circular blue phosphorus sheet is indented with a spherical indentor. The boundary is fixed (but free to rotate) as shown in Fig. 7a.

The membrane is pre-strained by different magnitudes in order to stabilize it. For simplicity, only one quarter of the specimen is modeled using the FE formulation of Sauer et al. [18]. The phosphorene specimen and indentor radii are 500 nm and 16.5 nm, respectively. Quadratic Lagrangian meshes with 2536, 2670, 3366, 4462 and 6506 finite elements are used for the convergence study (see Fig. 7b for the coarsest mesh). The reaction force only changes by about 0.05% between the finest and second finest mesh. The deformed geometry and force-displacement graph are shown in Figs. 8 and 9.

The force-displacement curve is linear for high pre-strains but shows pronounced nonlinearity otherwise. The curves terminate at a certain load level due to material instability (i.e. failure).

6.2 Phophorene vibration

Next, the variation of the frequencies for a square phosphorene sheet is obtained analytically. Since the model is validated in Sec. 5.2, we expect the results presented here to be correct. The frequencies of a pre-strained rectangular sheet with edge lengths $a$ and $b$ are [28, 23]

[TABLE]

where $\omega_{(m,n)}^{2}$ are the square of the frequencies and $\hat{\omega}^{2}_{(m,n)}$ are related to the bending stiffness, which can be neglected for large sheets ($\hat{\omega}_{(m,n)}\approx 0$) [23], $\rho=1.101\times 10^{-6}$ kg/m2 is the surface mass density, $a$ and $b$ are the half length and half width of the sheet and $\sigma_{x}$ and $\sigma_{y}$ are the normal surface stress components along the edges $a$ and $b$ [28]. $m$ and $n$ are the number of half waves in the mode shapes along the edges $a$ and $b$.

Under pure dilatation, $\sigma_{x}$ and $\sigma_{y}$ can be written in closed form as

[TABLE]

Based on this expression, the variation of the frequencies for a square sheet with $a=b=250$ nm follows as is shown in Fig. 10a.

Under uniaxial stretch along the armchair and zigzag directions, $\sigma_{x}$ and $\sigma_{y}$ can be computed from (22)777For uniaxial stretch, $\sigma_{x}=\sigma_{11}$ and $\sigma_{y}=\sigma_{22}$, where $\sigma_{11}$ and $\sigma_{22}$ are the Cartesian components of $\sigma$, which are computed from (22).. Based on this, the variation of the frequencies against $\lambda_{1}$ follows as is shown for a square sheet in Fig. 10b. The material behavior is anisotropic for uniaxial stretch and the frequencies are increasing faster if the sheet is stretched along the zigzag direction instead of the armchair direction. The frequencies increase monotonically up to the instability point, beyond which the material is unstable.

7 Conclusions

A new hyperelastic material model is proposed for blue phosphorus ($\beta$-P). The model is fully nonlinear and captures anisotropic behaviour of the material. It is written based on a set of invariants that are obtained from the symmetry group of the lattice. The continuum model is calibrated with DFT data. The model is in good agreement with further data sets from DFT results. The model is implemented in the curvilinear finite element membrane formulation of Sauer et al. [18] and thus used for the simulation of nano-indentation. To the best of our knowledge, there is no atomistic potential for blue phosphorus atoms and the proposed model is the only way to simulate specimens at micro-scale. The proposed material model can be extended to finite temperatures based on the new anisotropic thermoelastic shell formulation of Ghaffari and Sauer [25].

Authors’ contributions

RG and RAS derived the continuum material model and designed the continuum examples. RG implemented the finite element formulation, ran the numerical continuum examples, produced all continuum results, and designed the quantum experiments. FS conducted all density functional theory simulations. RAS and MH supervised and directed the research together. The paper has been written by all authors.

Competing Interests

The authors have no conflict of interests.

Funding

Financial support from the German Research Foundation (DFG) through grant GSC 111 is gratefully acknowledged.

Acknowledgement

The authors would like to thank Mr. Mohammad Sarkari Khorrami for helpful discussions.

Appendix A DFT simulations

Monolayer blue phosphorus is a two-dimensional material with hexagonal structure and lattice parameter of 3.28 Å. The thickness of the monolayer, due to out-of-plane distance of neighbouring phosphor atoms, is 1.23 Å. The atomic structure of blue phosphorus is shown in Fig. 2 and Tab. 1. DFT calculations of the mechanical properties of blue phosphorus are conducted using the Quantum ESPRESSO package [29, 30]. An ultrasoft pseudopotential is used to approximate the effect of non-valence electrons and the exchange correlation energy is approximated using the Perdew, Burke, and Ernzerhof (PBE) exchange-correlation functional [31]. The periodic unit cell contains 2 phosphor atoms and – in order to eliminate the interaction with replicates in $z$ direction – an interlayer spacing of 25 Å is chosen. The Brillouin zone integration is performed within a Monkhorst-Pack method with a $15\times 15\times 1$ k-mesh. The kinetic energy cutoff is set to 45 Ry and 360 Ry for wavefunctions and charge density, respectively. The convergence threshold for self-consistency is set to $10^{-8}$ Ry. Geometries have been optimized by the BFGS algorithm until the forces on each atom are less than $10^{-5}$ $\text{R}_{\text{y}}/\text{r}_{\text{Bohr}}$.

Appendix B Derivative of the invariants with respect to logarithmic strain

The first and second derivatives of the invariants $\mathcal{J}_{i}$ with respect to the logarithmic strain $\boldsymbol{E}^{(0)}$ are needed in the derivation of the stress and elasticity tensors. The first derivatives are [20]

[TABLE]

[TABLE]

[TABLE]

while the second derivatives are

[TABLE]

[TABLE]

[TABLE]

where $\mathcal{O}$ is the fourth order zero tensor and $\mathcal{I}$ is the fourth order identity tensor (see Kumar and Parks [20] for $\mathcal{O}$ and $\mathcal{I}$).

Appendix C Finite difference computation of the stress and elasticity tensor

The contra-variant components of the Kirchhoff surface stress tensor $\tau^{\alpha\beta}$ and the corresponding elasticity tensor $c^{\alpha\beta\gamma\delta}$ are needed for the curvilinear membrane FE formulation of Sauer et al. [18]. Here a central difference scheme is proposed for their calculation. $W$ is a function of the logarithmic surface strain tensor $\boldsymbol{E}^{(0)}$. $\boldsymbol{E}^{(0)}$ can be written as $\boldsymbol{E}^{(0)}=1/2\ln\boldsymbol{C}$, where $\boldsymbol{C}=a_{\alpha\beta}\,\boldsymbol{A}^{\alpha}\otimes\boldsymbol{A}^{\beta}$. Thus $W(a_{\alpha\beta})$ and accordingly, $\tau^{\alpha\beta}$ can be computed as

[TABLE]

where $\left[a_{\alpha\beta}\right]$ is the matrix form of the current surface metric and $[\Delta a_{\alpha\beta}^{+}]$ are

[TABLE]

and $[\Delta a_{\alpha\beta}^{-}]=-[\Delta a_{\alpha\beta}^{+}]$. The corresponding elasticity tensor is

[TABLE]

where $\tau^{\alpha\beta}(\Delta a_{\gamma\delta}^{+})$ and $\tau^{\alpha\beta}(\Delta a_{\gamma\delta}^{-})$ are defined by

[TABLE]

and

[TABLE]

In these formulae, $a_{\alpha\beta}$ and $\Delta a_{\alpha\beta}$ can be replaced by the Cartesian components of $\boldsymbol{C}$ and its increment $\Delta\boldsymbol{C}$, in order to compute the Cartesian components of the second Piola-Kirchhoff surface stress tensor $\boldsymbol{S}$ and its conjugate elasticity tensor $\mathbb{C}$ since $\boldsymbol{S}=2\partial W/\partial\boldsymbol{C}$ and $\mathbb{C}=\displaystyle 4\partial^{2}W/\partial\boldsymbol{C}\partial\boldsymbol{C}$. In the computation of $\tau^{\alpha\beta}(\Delta a_{\gamma\delta}^{+})$, the components of $[a_{\alpha\beta}]+[\Delta a_{\alpha\beta}^{+}]+[\Delta a_{\gamma\delta}^{+}]$ are summed. For example as

[TABLE]