六角形二维材料的超弹性建模：一种基于互补能量的公式化方法，该方法考虑了更高阶的柔顺性特性

《Mechanics of Materials》：Hyperelastic modeling of hexagonal two-dimensional materials: A complementary-energy formulation with higher-order compliances

　　摘要
二维（2D）材料，如石墨烯，在纳米压痕和拉伸测试中可以发生超出线性弹性范围的变形，这需要使用非线性本构模型。现有的超弹性公式将应变能量表达为应变的函数，因此能够得到弹性刚度。相比之下，顺应性公式将应变表示为应力的函数，更适合于应力已知或由平衡状态确定的问题。本工作提出了一种针对六角形2D材料的超弹性本构模型，其中补充应变能量密度在第二个Piola–Kirchhoff应力项中展开到五阶，并将Green–St. Venant应变作为其工作共轭导数来恢复。六角对称性将五阶顺应性展开简化为十四个独立常数。这些常数是根据使用第一性原理计算和经典原子间势能沿着结合了双向拉伸和剪切的变形路径计算的应变能量和应力数据来确定的。另一种方法是首先拟合一个刚度展开式，然后将其反转来获得顺应性，在低阶时是准确的，但在高阶时条件较差。拟合得到的顺应性可以给出应力依赖的杨氏模量、泊松比和剪切模量的闭合形式表达式。作为应用示例，推导出了预拉伸膜条在压痕过程中的载荷-变形关系的闭合形式，这些关系分别达到了二阶和三阶。通过密度泛函理论得到的石墨烯和二硫化钼的顺应性与AIREBO-M和Stillinger–Weber势能得到的顺应性进行了比较，并且分析预测结果通过直接的分子静态压痕模拟得到了验证。补充能量展开式直接给出了更高阶的弹性顺应性，并为六角形2D材料的非线性弹性建模提供了本构框架。

引言
二维（2D）材料是原子级薄的晶体固体，只有几埃厚，具有很强的面内结合（Tan等人，2017年；Shanmugam等人，2022年；Das等人，2014年）。自从石墨烯被发现以来（Novoselov等人，2004年），这类材料迅速扩展到包括六角形氮化硼（Novoselov等人，2005年）、过渡金属硫属化合物（如二硫化钼MoS2）（Mak等人，2010年）、磷烯（Li等人，2014年）、硼烯（Feng等人，2016年）以及IV族单硫属化合物（如硒化锡SnSe）（Gomes和Carvalho，2015年），高通量计算筛选预测了更多候选的单层材料（Mounet等人，2018年）。尽管只有几原子厚，它们却属于已知最硬和最强的材料之一（Lee等人，2008年；Akinwande等人，2017年）。它们还能在较大应变下保持弹性，例如石墨烯在断裂前仍然保持弹性（Lee等人，2008年），远远超出了线性弹性范围。在许多与柔性电子器件（Kim等人，2015年）和应变工程（Deng等人，2018年）相关的变形范围内，它们的机械响应是非线性的，其中故意施加有限的晶格应变以调节电子和光学性质。因此，准确的本构描述对于评估2D材料在大变形下的机械性能以及设计应变调谐器件至关重要。
现有的二维材料非线性弹性响应的超弹性本构模型大致分为两类。第一类将应变能量密度表示为Green–St. Venant应变分量的多项式。Cadelano等人（2009年）通过紧束缚计算从三次应变能量展开式中得到了石墨烯的非线性弹性常数，而Wei等人（2009年）通过第一性原理计算表明，需要五阶项来捕捉其大应变响应，采用了Murnaghan展开式（Murnaghan，1937年；Birch，1947年；Hiki，1981年）。此后，针对不同对称性的2D材料（如六角形氮化硼（Peng等人，2012年）、二硫化钼（Cooper等人，2013年）、碳化硅（Peng，2020年）、硼烯（Faghihnasiri等人，2019年）和磷烯（Maalouf和Vel，2022年）提出了这种多项式本构模型，独立常数的数量由对称性决定，并随着展开阶数的增加而急剧增加。第二类公式将应变能量密度表示为应变不变量的函数，这种方法直接在函数形式中编码了材料对称性，并通常需要较少的常数。Xu等人（2012年）使用Green–St. Venant应变的各向同性不变量对石墨烯进行了建模，并通过切换函数来捕捉大应变各向异性。Kumar和Parks（2015年）基于对数（Hencky）应变张量的三个对称不变量构建了一个超弹性模型（Hencky，1928年；Ogden，1997年），H?ller等人（2020年）将该框架适应于石墨烯的Green应变不变量。Maalouf和Vel（2025a）引入了物理启发的Green应变张量不变量，以分离六角形2D材料中的膨胀、偏应力和各向异性贡献，并随后将该方法推广到Seth–Hill应变测量家族（Vel和Maalouf，2026年）。
在多项式和不变量公式中，应变能量密度都被写为应变的函数，因此得到的弹性常数是刚度，它们将施加的应变与产生的应力相关联。当应变是已知的时候，这样的模型非常方便。然而，在许多应用中，已知的是应力，无论是直接施加的还是从平衡状态确定的。一个常见的例子是沿轴向受单轴应力的2D材料条带。然后横向应变是未知的，要从刚度关系中恢复它需要在横向应力为零的条件下求解非线性应力-应变方程。相比之下，基于顺应性的公式将应变张量表示为应力张量的函数，未知的应变直接从已知的应力中得出。
通常报告的2D材料的工程模量也都是在对称应力下定义的。在线性极限下，杨氏模量是轴向顺应度的倒数，泊松比是顺应度的比值。它们的有限应力行为来自于更高阶的顺应性。这些更高阶的顺应性本身就是弹性常数，但很少直接确定，通常是通过反转拟合的刚度得到的。在线性阶数下，这种反转仅仅是刚度的矩阵逆。超过线性阶数，Barsch（1968年）和Kube（2022年）分别将三阶和四阶顺应性表示为相应刚度的收缩，类似的关联随着阶数的增加而变得更加复杂。这些关系在每个阶数上都是代数精确的，但它们放大了拟合刚度的不确定性，特别是对于分辨率较低的更高阶常数。因此，反转得到的更高阶顺应性的条件较差，这促使人们直接确定它们而不是从刚度中恢复它们。
在本工作中，直接确定了2D材料的更高阶顺应性，而不需要反转。该公式基于Green–St. Venant应变和第二个Piola–Kirchhoff应力形成的有限应变工作共轭对。然后，补充应变能量密度（W?）作为应变能量（W）的勒让德变换得到，并作为第二个Piola–Kirchhoff应力的多项式展开，其中顺应性作为展开式中的系数出现。六角对称性将二阶模型的独立常数数量减少到两个，五阶模型的独立常数数量减少到十四个，并且依赖的顺应性关系被明确推导出来。为了确定顺应性，沿着结合了双向拉伸和剪切的路径采样变形状态，并使用密度泛函理论（DFT）和经典原子间势能计算应变能量和应力。然后通过勒让德变换结合采样的应变能量、应力和应变来获得补充能量，进而通过对称性简化的线性最小二乘回归得到顺应性。这个过程与用于刚度的应变能量拟合类似（Maalouf和Vel，2022年；Maalouf和Vel，2023年），使顺应性和刚度处于同一基础。相同的采样状态适用于这两种拟合，无需应力控制计算即可获得顺应性。拟合得到的顺应性以闭合形式给出了应力依赖的工程模量，包括杨氏模量、泊松比和剪切模量的正割和切线形式。
报告了两种六角形单层材料——石墨烯和二硫化钼的二阶到五阶顺应性及其相应的工程模量，并与经典原子间势能AIREBO-M（O’Connor等人，2015年）（对于石墨烯）和Stillinger–Weber（Stillinger和Weber，1985年；Wen等人，2017年）（对于二硫化钼）得到的常数进行了比较。这些势能也为条带压痕问题提供了原子级参考，因为在所需系统规模下直接进行第一性原理模拟在计算上是不可行的。刚度-顺应性反转关系是通过五阶推导出来的，并将反转得到的顺应性和直接拟合得到的顺应性与DFT数据进行了比较。两种方法在低阶时一致，在高阶时直接拟合得到的顺应性更能准确地再现数据。
该公式应用于预拉伸的2D材料条带的横向压痕，这是一个完全非线性的边界值问题，直接检验了顺应性形式的有效性。压痕是测量悬浮2D材料面内弹性响应的标准方法（Cao和Gao，2019年），已用于测量石墨烯（Lee等人，2008年）和二硫化钼（Bertolazzi等人，2011年）的刚度，这两种材料也是本研究的对象。这种悬浮膜的弹性响应也可以通过静电驱动的鼓的非线性动态响应来提取（Davidovikj等人，2017年）。大多数研究采用圆形膜，而这里检查的双端固定条带遵循经典的夹紧梁配置（Pruessner等人，2003年）。同样的配置也被用于表征悬浮在光刻沟槽上的石墨烯（Frank等人，2007年）和由无尖端原子力显微镜悬臂加载的网状碳纳米管条带（Wang等人，2016年）。条带沿其长度预拉伸，其中间的横向载荷使其变形为对称的V形配置。V形运动学通过变形固定了轴向应变，平衡关系将压头载荷与膜张力联系起来，它们之间的本构联系由单轴应力顺应性多项式提供。在二阶和三阶，载荷-变形关系以闭合形式获得，其中顺应性多项式分别简化为应力的线性和二次方程，并在高阶时进行数值评估。这些预测与使用相应经典势能计算得到的石墨烯和二硫化钼条带的直接分子静态模拟进行了比较。由于条带模拟不是用来拟合顺应性的，因此比较测试了拟合的顺应性关系是否适用于独立的边界值问题。
本文的其余部分组织如下。第2节介绍了补充能量公式、对称性简化、正割和切线工程模量以及通过刚度反转恢复顺应性的方法及其在高阶时的局限性。第3节推导了预拉伸条带压痕的解析载荷-变形关系，提供了后来用于验证顺应性公式的非线性边界值问题。第4节描述了DFT和经典原子间势能计算以及通过射线采样提取顺应性的方法。第5节介绍了石墨烯和二硫化钼的顺应性，并将条带解与分子静态模拟进行了验证。最后，第6节提出了总结和结论。

超弹性公式
图1展示了一种2D材料A及其等效的连续体C，它们具有相同的面内尺寸和机械响应。在这种等价假设下，连续体中存储的应变能量密度W与2D材料在相同变形状态下的应变能量密度WA相等，即WA=WC。图1中显示的正交坐标基向量e1、e2和e3定义了参考坐标系，未变形的基态为

预拉伸2D材料条带在压痕下的非线性载荷-变形响应
本工作中提出的基于顺应性的本构模型通过将其应用于一个完全非线性的边界值问题——预拉伸的六角形2D材料条带的压痕进行了验证。该条带的参考长度为2L0，宽度为W0，如图2所示，在两个锚点处被固定。压痕之前，锚点向外移动以施加均匀的轴向拉伸λ0，使锚点之间的距离扩大到2λ0L0。条带处于单轴应力状态

计算方法
顺应性是根据在一定范围内均匀变形中采样的应变、应力和能量数据确定的。这些数据可以基于第一性原理或原子间势获得，在本研究中两种方法都被采用了：第4.1节描述了基于第一性原理的计算，而第4.2节则描述了使用经典原子间势的计算方法。**结果与讨论**：本节展示了两种代表性的六角形二维材料——石墨烯（graphene）和二硫化钼（MoS2）的本构建模框架。石墨烯是已知最坚硬且研究最广泛的二维材料之一，提供了大量可用于验证的数据；而MoS2则是一种更具柔性的过渡金属硫属化合物，具有三层原子结构。**总结与结论**：本文为六角形二维材料在有限变形条件下的行为开发了一个超弹性本构模型，在该模型中，主要弹性常数是弹性柔顺性而非刚性。相关的应变能密度W?被展开为第二Piola–Kirchhoff应力的多项式形式，通过对应力求导可以得到Green–St. Venant应变。六角对称性使得五阶展开简化为十四个独立的柔顺性参数。**利益冲突声明**：作者声明不存在可能影响本文研究结果的已知财务利益冲突或个人关系。**致谢**：本工作得到了缅因大学工程学院Arthur O. Willey教授职位的支持。作者感谢Serge Maalouf对早期草稿提出的建设性意见。
Senthil S. Vel
Arthur O. Willey教授
机械工程系
缅因大学
奥罗诺，美国