单层MoS2力学性能的有限元分析

杨 璐, 齐 博, 杨宏旭, 王天韵

(沈阳工业大学 建筑与土木工程学院, 沈阳 110870)

单层MoS2是近些年发展起来的新型材料,属于过渡金属硫化物,是一种类石墨烯结构和性能的二维材料,主要采用微机械力剥离、离子插层和液相剥离法等为主的“自上而下”的剥离法,以及以化学气相沉积法、水热/溶剂热合成法等为主的“自下而上”的合成制备方法.该材料具有优异的电学、光学以及力学性能,具有广阔的市场前景,可用作柔性光电器件、光电探测器、类脑逻辑器件和催化剂等[1-3].

单层MoS2的半导体性质允许其克服石墨烯的零带隙,同时分享石墨烯在电子应用中的许多优势[4].为了确保基于单层MoS2的传感器和纳米材料在制作和使用期间保持其结构完整性,防止其出现机械失效,有必要对单层MoS2的力学性能进行更全面的了解.对于MoS2力学性能的研究大多数人采用了原子力显微镜(AFM)方法、密度泛函理论(DFT)、分子动力学方法(MD)得出其弹性模量、泊松比、破坏强度等力学参数[5-10].采用AFM试验方法对所需仪器的要求较高,采用DFT和MD方法的模拟计算代价太大,经常用来模拟少量分子和原子的系统,以及相对短寿的现象.

部分学者还在致力于寻找其他高效准确的方法来研究纳米材料的力学性能,Li等[11]提出在三维空间中用空间框架结构来等效石墨烯,采用有限元方法得出了石墨烯的弹性模量以及剪切模量,为其他共价键合材料的模拟提供了一个新思路;韦琦等[12]计算了石墨烯的弹性模量、泊松比、剪切模量等力学参数;张续等[13]利用石墨烯的Morse势能函数作为C—C键的本构方程,得出了含缺陷石墨烯的抗拉强度.

同为二维材料的MoS2与石墨烯有相似的结构,为了探究出更加高效的方法,采用有限元法模拟MoS2力学性能并探究其正确性.本文计算了MoS2的弹性模量,并得出单层MoS2弹性模量具有尺寸效应的结论,通过ABAQUS软件模拟出了MoS2的破坏强度,为采用有限元法研究MoS2提供理论借鉴.然而采用有限元法的计算成本较低,只要采用足够精确的化学键拉伸、弯曲、扭转常数等参数,在对MoS2的屈曲、振动、拉伸断裂等线性或非线性范围内的力学特征进行模拟时也能通过有限元法得到较为精确的结果.

1 有限元模型建立

1.1 MoS2模型建立

单层MoS2中间一层由Mo原子的一个平面组成,两边连着两个S的原子平面,具体的原子结构如图1所示.在有限元计算中采用杆系结构来模拟Mo原子和S原子之间的共价键,杆系结构之间采用刚性连接.本文采用AutoCAD完成三维建模,再导入到大型商用有限元模型ABAQUS中进行计算,模型如图2所示,键长为0.241 nm[14].

图1 MoS2 球棍模型

图2 MoS2有限元模型

1.2 线性材料参数确定

在有限元模拟中采用圆形截面梁单元模拟Mo和S之间的共价键,共价键键长r0=0.241 nm.由于有限元方法需要材料属性以及单元截面参数,所以需要将微观分子计算与宏观结构力学建立联系.在分子力学中分子总能量等于分子势能与分子动能相加,总能量可表示为

U=∑Ur+∑Uθ+∑Uτ+∑Uω+∑Uvdw

(1)

式中:Ur为化学键拉伸能;Uθ为化学键弯曲能;Uτ为化学键扭转能;Uω为平面外扭转能;Uvdw为非键范德华相互作用.

一般来说,位能的作用主要来自前四项,为了方便研究,将平面内扭转能与平面外扭转能合并成一项,采用谐振子势表示各势能函数,可得

(2)

(3)

(4)

式中:kr为化学键拉伸长度;kθ为化学键弯曲常数;kτ为化学键扭转常数;Δr、Δθ和Δφ分别为化学键拉伸增量、键角转动增量和键扭转增量.

在结构力学中,框架结构中的梁也可以产生类似的拉伸变形和弯曲变形.在弹性小变形范围内,拉伸应变能、弯曲应变能、扭转应变能分别表示为

(5)

(6)

(7)

式中:A为截面面积,A=πd2/4,d为截面直径;I为截面惯性矩;E为弹性模量;ΔL、Δα和Δβ分别为梁轴向拉伸变形量、梁端点的转角、梁端点的扭转角;L为梁长度;F、M、T分别为梁单元所受的拉伸作用力、弯矩、扭矩;G为剪切模量;J为截面极惯性矩.

将分子势能与应变能等效,可得

(8)

(9)

(10)

进而得出

Li等[11]通过计算kτ得到石墨烯的泊松比不在合理范围内,且kτ对弹性模量的影响很小.假设圆截面梁单元的泊松比为0,Nicolini等[14]将不同谐振子势对比,认为谐振子势参数在预测MoS2的结构和振动性能方面有较好的应用前景,其参数为

计算得到截面直径d=0.317 nm,弹性模量E=452.859 GPa.

1.3 非线性材料参数确定

本文对理想单层MoS2的弹性性能进行了详细描述,但在弹性范围内研究MoS2纳米带在大变形下的拉伸断裂等性能是不准确的.因此,采用结构力学的方法来研究MoS2的非线性特征,需要将MoS2的非线性本构等效成圆形截面梁的本构方程.

与研究单层MoS2纳米带弹性性能类似,建立同样的空间框架结构.MoS2纳米带非线性模型的建立主要体现在Mo—S键的非线性上,采用Morse势参数作为等效圆截面梁的本构方程[14].Morse势的总方程可以表示为

E=Estretch+Eangle

(11)

Estretch=De[1-e-α(r-r0)]2

(12)

F=2Deα(1-e-αΔr)e-αΔr

(13)

F=2Deα(1-e-αr0)e-αr0

(14)

图3 圆形截面梁应力应变曲线

张续[13]和王少培等[15]在采用Morse势研究石墨烯时分别采用在应力最大值的应变以及应变为3%前的曲线斜率作为弹性模量值.经过计算参数的取值均合理,表明材料屈服时应变的取值对结果影响不大.本文将应变为4%视为材料屈服.在应变达到16%时,此时Mo—S键间应力值达到最大,承受的力最大为2.273 nN.当原子间的距离增大时,原子之间的吸引力逐渐变大,2.273 nN即为原子间最大的吸引力.随后原子间的作用力逐渐减小,在Mo—S键承受的力最大时,此时的Mo—S键可以认为发生断裂,等效为Mo—S键的圆形截面梁发生断裂.计算得到非线性阶段的屈服应力和应变,通过Mo—S键的Morse势能函数推导出结构力学中Mo—S键等效成的圆形截面梁的弹塑性本构.通过在ABAQUS中输入圆形截面梁的弹塑性本构,计算其拉伸破坏时的破坏应力以及应变.

1.4 荷载施加

扶手型和锯齿型的区别如图4所示.沿着x轴方向的边界形状似锯齿状,所以研究MoS2在x轴方向上的性质表述为锯齿型MoS2;沿着y轴方向的边界形状似扶手状,所以研究MoS2在y轴方向上的性质表述为扶手型MoS2.

图4 MoS2手性的定义

采用线性本构关系时,该方法是在弹性范围内应变能与分子势能等效,在施加力时只能施加小变形才能保证计算相对准确.采用非线性本构关系时,同样的约束方式,在另一端可施加较大的位移荷载.

荷载施加示意图如图5所示.采用一端约束,另一端单向拉伸的方法.对于扶手型MoS2,在最右侧原子上沿x方向施加0.001 nm的微小形变,在另一侧施加x方向以及z方向上的约束.对于锯齿型MoS2,在最上边每个原子上沿y方向施加0.001 nm的微小形变,在另一侧施加y方向以及z方向上的约束.

图5 单层MoS2加载示意图

2 线性本构下弹性模量的计算

2.1 模型尺寸确立

采用8种不同尺寸的单层MoS2,分别计算其在扶手型和锯齿型方向下的弹性模量大小.由于x方向与y方向的长度不是严格意义上的相等,所以在保证结构完整的情况下,尽可能接近矩形大小,基于此8个工况的具体尺寸如表1所示.

表1 各工况MoS2尺寸以及单元信息

2.2 弹性模量计算

在单向拉伸试验中,弹性模量是应力与应变的比值.对于扶手型MoS2,通过施加小变形,计算所对应方向支座反力Fx,并通过式(15)计算弹性模量大小.对于锯齿型MoS2,通过施加小变形,计算所对应方向的支座反力Fy,并通过式(16)计算弹性模量大小.

(15)

(16)

式中:t为MoS2薄片厚度;Lx为x方向长度;Ly为y方向长度;ΔLx为施加在x方向上的小变形;ΔLy为施加在y方向上的小变形;A0和A1分别为横向和纵向MoS2薄片横截面面积大小,采用Lxt和Lyt计算.无论是对石墨烯还是MoS2的研究中,其厚度的选取存在着较大争议,厚度的合理选取也是有限元模拟结果是否准确的关键,本文选取模型的厚度为0.61 nm[16].

2.3 数值模拟计算结果与分析

基于荷载施加方法,以及式(15)、(16)分别计算不同尺寸下扶手型和锯齿型MoS2的弹性模量.为了更直观地看出其变化,计算相邻两尺寸下MoS2纳米带弹性模量的相对误差,其表达式为

(17)

表2 单层MoS2弹性模量

图6 不同尺寸MoS2弹性模量变化曲线

由表2和图6可以看出,扶手型MoS2和锯齿型MoS2的弹性模量大小都随尺寸的增加逐渐减小,并且相邻两工况之间的相对误差也在逐渐减小.工况1两种类型的MoS2弹性模量均达到了最大值,分别为202.78 GPa和173.07 GPa.工况8两种MoS2弹性模量值均为8种工况下的最小值,分别为145.03 GPa和141.07 GPa.

扶手型和锯齿型MoS2的相对误差值都是逐渐减小的,最小值分别达到了1.08%和0.56%,可以看出,随着尺寸的增大,尺寸效应逐渐弱化.在工况8中无论是扶手型MoS2还是锯齿型MoS2的弹性模量均已趋于平稳.

从整体来看,在每种工况下,扶手型MoS2均比锯齿型MoS2的弹性模量大,这表明扶手型MoS2抵抗变形的刚度较大,具有更好的抗拉性能.锯齿型MoS2在每两种工况之间的相对误差比扶手型MoS2小.扶手型和锯齿型MoS2在工况1与工况2之间的相对误差达到最大,分别为14.42%和9.03%.扶手型和锯齿型MoS2在工况7和工况8之间的相对误差达到最小,最小值分别为1.08%和0.56%,表明锯齿型MoS2的尺寸效应没有扶手型MoS2明显.

采用单向拉伸不同工况下MoS2纳米带的弹性模量可以清晰地看到,无论是扶手型还是锯齿型MoS2的弹性模量都存在明显的尺寸效应.

3 非线性本构下破坏强度的计算

3.1 破坏应力应变的计算

拉伸过程中一端固定,另一端施加较大的位移荷载,大到Mo—S键达到其破坏应力且发生断裂,整个MoS2薄片无法再继续承受拉力,通过模拟结果计算MoS2薄片的应力和应变.MoS2薄片的应力和应变计算公式为

(18)

(19)

3.2 数值模拟结果与分析

图7 MoS2拉伸破坏过程应力应变曲线

扶手型MoS2的破坏应力达到18.9 GPa,此时对应的破坏应变为22.5%.锯齿型MoS2的破坏应力达到16.0 GPa,对应的破坏应变为21.2%.

4 结 论

本文采用有限元方法,利用杆系结构模拟共价键,分别采用线性本构关系以及非线性本构关系进行单轴拉伸,得出了MoS2纳米带的弹性模量以及破坏应力应变.

1) 采用线性本构关系,单轴拉伸得到MoS2的弹性模量具有尺寸效应,随着尺寸的增大,弹性模量不断减小.扶手型MoS2的弹性模量为145.03～202.78 GPa,锯齿型MoS2的弹性模量为141.07～173.07 GPa.扶手型MoS2的尺寸效应更为明显.

2) 采用非线性本构关系,单轴拉伸得到扶手型MoS2的破坏应力为18.9 GPa,破坏应变为22.5%.锯齿型MoS2的破坏应力为16.0 GPa,破坏应变为21.2%.扶手型MoS2的破坏应力和应变较锯齿型MoS2更大,力学性能更好.