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

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.

💡 Research Summary

This paper presents a fully nonlinear, anisotropic hyperelastic membrane model for single‑layer blue phosphorus (β‑P), also known as blue phosphorene, directly calibrated from density functional theory (DFT) calculations. The authors begin by motivating the need for a continuum description of 2‑D materials that can capture large‑strain behavior and lattice symmetry, noting that existing models for black phosphorus lack periodicity and fail for arbitrary loading directions.

The kinematic framework adopts curvilinear surface coordinates and defines the surface deformation gradient F, its polar decomposition (R U), and the logarithmic surface strain E^(0) = ln U = ½ ln C. E^(0) is split into an area‑changing part (E^(0)_area) and a deviatoric part (E^(0)_dev). For the hexagonal lattice of blue phosphorus, the authors construct structural tensors consistent with a six‑fold rotational symmetry (C₆ᵥ) and derive three independent invariants:

• J₁ = ln J, the logarithmic area change,
• J₂ = ½ E_dev : E_dev = (ln λ)², measuring isotropic shear,
• J₃ = (ln λ)³ cos 6θ, capturing the six‑fold anisotropy, where λ is the principal stretch ratio and θ is the angle between the principal stretch direction and the armchair axis.

The strain‑energy density per unit reference area is proposed as
W = f₁(J₁) + f₂(J₁) J₂ + f₃(J₁) J₂² + f₄(J₁) J₃,
where each coefficient f_i is a fourth‑ or fifth‑order polynomial in J₁. This formulation introduces 20 material parameters (n_i, μ_ij, η_i) that are identified by minimizing a least‑squares cost function that measures the discrepancy between the model predictions and DFT data for three deformation modes: pure dilatation, uniaxial stretch along the armchair direction, and uniaxial stretch along the zigzag direction. A total of 520 DFT data points are used, yielding the parameter sets listed in Tables 1–4.

Calibration results show that under pure dilatation the energy density increases monotonically with J₁, while the surface tension reaches a maximum and then declines, indicating the onset of lattice instability. In uniaxial tests, the model reproduces the anisotropic stress‑strain curves observed in DFT: the armchair direction exhibits lower failure strain and higher σ₁₁, whereas the zigzag direction shows higher σ₂₂. The agreement is excellent up to the material’s instability point.

To account for out‑of‑plane thickness change, the authors introduce an isotropic function λ₃(J₁, J₂, J₃) calibrated from the same DFT dataset, which behaves almost isotropically for both dilatation and uniaxial loading.

For finite‑element implementation, the logarithmic stress S^(0) = ∂W/∂E^(0) and the associated elasticity tensor C^(0) are derived analytically. Because transforming these quantities to the Cauchy stress σ and Kirchhoff stress τ involves complex tensor operations, the authors adopt a finite‑difference scheme to compute the required push‑forward operations efficiently, as detailed in the appendices.

The calibrated model is then embedded in a rotation‑free isogeometric shell formulation (based on Sauer et al.) and applied to two benchmark problems. First, an indentation simulation of a micrometer‑scale blue‑phosphorene sheet reproduces the force‑depth curve obtained from atomistic DFT indentation, confirming that the continuum model captures the correct stiffness and failure behavior. Second, a free‑vibration analysis of a suspended sheet yields natural frequencies and mode shapes that match analytical predictions and prior molecular dynamics studies. These examples demonstrate that the model can bridge the gap between quantum‑scale data and macroscopic engineering analyses.

In conclusion, the paper contributes (1) a symmetry‑consistent invariant formulation for hexagonal 2‑D crystals, (2) a DFT‑calibrated high‑order hyperelastic energy function capable of large‑strain anisotropic response, (3) an efficient computational strategy for stress and tangent transformations, and (4) validation through large‑scale shell simulations. The methodology is readily extensible to other 2‑D materials (e.g., graphene, black phosphorus) and provides a foundation for incorporating temperature effects and multi‑physics couplings in future work.