点缺陷对单晶硅力学性能影响的分子动力学研究

万 伟， 唐昌新， 项永康， 邱 安

(南昌大学 光伏研究院， 南昌 330031)

1 引 言

硅的丰度在地球上仅次于氧，储量丰富，易于被人类提取与利用. 随着科技的发展，单晶硅已在大规模集成电路和太阳能电池等领域获得了广泛的应用，例如：单晶硅晶圆作为太阳电池等产品的原料，是由线锯工艺从硅锭上切割而成[1]. 然而，在硅锭的生产加工过程中会产生崩边、隐裂和碎片等不良品[2]，降低了原材料的利用率，提高了经济成本，不利于单晶硅产业的高速发展，更会对其他以硅材料为基础的产业产生负面影响. 因此，有必要对单晶硅的力学性能展开研究，探究产生切割不良品的原因.

从单晶硅力学性能面临的问题入手，Liu等人[3]采用压痕法等实验手段对单晶硅加工过程中产生的变形与表面损伤等问题进行研究，并对比结合分子动力学等手段所给出的模拟值来验证晶体亚结构的变形与破坏理论. 由于分子动力学模拟方法可以提供原子微观结构演变的详细信息，为理解材料微观性能提供了一个有力的研究工具. 在单晶硅产业发展的数十年历程中，通过多学科综合研究等一系列方法，我们对氮[4]、氧[5]等杂质原子以及位错[6]等因素对硅力学性能的强化效应及其影响机理有了充分的认识. 另外在硅片的生产过程中，由于热应力、快速结晶等因素，单晶硅中会形成点缺陷、点缺陷团簇甚至原子规模的小孔洞，这些因素都会对产品的质量造成不可忽略的影响. 因此，本文将点缺陷纳入研究范围，对单晶硅力学性能的影响因素加以讨论和完善. 目前在点缺陷对金属[7, 8]、化合物[9-11]和碳基材料[12, 13]的力学性能研究中，已经通过理论分析、实验测试和分子动力学模拟等方式开展了较多高水平的工作. 作为单晶硅主要竞争对象的铸造多晶硅，点缺陷对其力学性能的影响亦有过系统性研究[14]. 但点缺陷对单晶硅的力学性能的影响如何，未见相关报道.

另一方面，单晶硅因其具有的高理论比容量和清洁环保等优点，被广泛应用在以硅基负极为首的新型电极材料上. 而硅在作为负极材料发生的嵌锂过程中会发生巨大的体积膨胀现象，造成电池首次库仑效率低和容量衰减等问题[15]. 因此，抑制硅负极电池充放电过程中的体积膨胀效应，探究单晶硅在循环应力作用下的力学性能变化情况，已经成为硅负极研究的关键问题之一. 为解决困扰硅负极电池的体积膨胀效应，Darbaniyan等人[16]研究了应变率和锂质量分数对硅晶体的力学性能影响. Kim等人[17]以预锂化的方式成功抑制了体积的变化，还有一些研究者以负极结构设计为着力点，提出了多种改善硅负极力学性能的优化方式. 例如，天津大学杨全红教授研究组在硅负极的构建中提出了一种通过可控收缩结合碳笼网络稳定碳/硅界面的力学缓冲增强策略[18]，有效地改善了嵌锂过程中的体积膨胀效应. Xie等人[19]将纳米硅颗粒包覆到富氮多孔炭基体中. 所制备的蜂窝状多孔复合材料具有显著的循环稳定性和优异的倍率性能，并且能够很好地适应嵌锂过程中的膨胀效应. 但尚未有研究者从原子级别的微观尺度，以分子动力学手段讨论在较大原子尺度的孔洞或点缺陷团簇作用下的单晶硅晶体力学性能的变化情况.

本文运用分子动力学方法，选择硅的典型势函数，在恒定应变率作用下模拟单晶硅的拉伸破坏行为，通过在晶体中构造不同原子大小的点缺陷，以此来分析单个点缺陷，点缺陷团簇乃至大型孔洞对单晶硅力学性能的影响. 同时利用数值模拟给出的数据验证了米塞斯强度理论，讨论了单晶硅中的应力集中效应与点缺陷之间的联系，阐明了点缺陷对单晶硅力学性能的影响机理. 期望能为单晶硅力学性能的同类研究问题提供参考.

2 模型与方法

2.1 模型构造

单晶硅具有金刚石型晶体构型，晶格常数为0.543 nm. 本文采用LAMMPS软件构建理想的单晶硅超晶胞模型，坐标轴X、Y、Z方向分别对应单晶硅的[100]、[010]、[001]晶向，计算区域取40×40×40个晶胞，包含512000个硅原子，模拟盒子几何大小为217×217×217 nm3. 模型在三个方向上均运用周期性边界条件，因为周期性边界条件可以消除表面效应，使模拟更加符合实际情况[20]. 构建完成后的理想晶体如图1所示. 体系建立后再在理想晶体的几何中心位置划分出不同半径的球形区域，以删除该区域内晶硅原子的方式构建不同大小的点缺陷. 点缺陷的大小则使用其生成时删除的晶硅原子数度量.

图1 单晶硅晶胞模型和6个含点缺陷的样品

2.2 势函数

分子动力学模拟是基于势函数描绘的相互作用进行的，因此选择一个能较好地描绘硅晶体性质的势函数对于模拟的准确性而言至关重要. 硅的原子间势函数发展至今已有三十多个版本[21]，经过一系列比较与参考后，本文决定使用Tersoff势的第二个版本进行分子动力学模拟. Tersoff 势是一种键级势, 首次报道于1986年, 后分别于1988年和1989 年先后经过两次修改, 形成了 T1[22], T2[23]和 T3[24]三个版本, 该势函数除了能较好地描述金刚石结构外, 对硅的非正四面体构型也能描述, 如团簇、晶向和液态等结构. 其系统的总势能为：

(1)

Vij=fc(rij)[aijfR(rij)+bijfA(rij)]

(2)

式中：fR与fA分别表示排斥作用函数和吸引作用函数，rij为原子i，j之间的距离，fc为光滑截断函数，作用是使两原子在距离较远时的相互作用光滑地趋向于零，aij和bij分别表示排斥和吸引作用的系数. Zhou等人[25]对比了Tersoff，SW和MEAM等势函数对硅晶体熔化特性的描述，结果表明Tersoff 势相对更适合描述硅的熔化和凝固过程.

2.3 模拟方法

单晶硅拉伸过程采用恒定应变率加载方式，应变率加载方向为X轴上的[100]晶向，应变率为1×108/S-1. 模拟过程分为弛豫与加载两个阶段. 由于初始构造的原子体系是不平衡的，为使体系达到平衡状态，将体系在NPT系综(等温等压恒定原子数目)弛豫300 ps，使其温度达到300 K. 弛豫完毕后开始加载应变：该过程采用NVE系综，利用Berendsen[26]恒温器控制温度，Berendsen恒压器控制Y轴与Z轴上的压力趋于零以消除边界压力对模拟结果的影响. 本模型时间步长设置为1 fs；使用Verlet算法计算原子的运动轨迹，结果的可视化分析则通过Ovito[27]软件完成.

表1 模拟参数设置

3 结果与讨论

3.1 点缺陷对力学性能的影响现象

材料的力学性能受晶体缺陷的影响最为严重[28]，而硅晶体又易受加工时产生的热应力等因素的影响而产生常见的点缺陷. 为分析点缺陷对单晶硅力学性能的影响现象，将所有单晶硅样本沿X轴方向施以恒定应变率加载的应力—应变曲线绘制为图2.

图2 所有样品的应力—应变曲线

理想晶体的应力-应变曲线表明，单晶硅在恒定应变率的拉伸过程包括弹性阶段和屈服阶段. 在弹性阶段中应变不断增加，应力也呈线性增加的趋势，此阶段中的变形为可恢复的弹性变形. 当应力超过屈服强度之后，晶体发生破坏，应力迅速下降至下屈服点，转入屈服阶段. 对比不同大小点缺陷的应力-应变曲线则表明，相较于不含点缺陷的理想晶体，点缺陷使单晶硅的屈服强度、极限应变显著降低，下屈服点升高.

进一步讨论点缺陷对屈服强度的影响时发现，点缺陷大小与屈服强度之间存在着函数关系. 图3的拟合结果表明，点缺陷大小与屈服强度σs的变化服从公式(3)所述的指数函数关系.

图3 屈服强度与点缺陷大小的指数拟合曲线

σs(c)=σ0+A×Exp(R×c)

(3)

式中，c为点缺陷的大小，σ0为屈服强度下限值，约为12.05±0.46 GPa，A=5.41±0.64，R=-0.0051±0.0015.A与R均是与点缺陷自身性质有关的参数. 至于上述两个参数的影响因素，还有待更深入的研究和讨论.

3.2 从应力变化的角度讨论点缺陷的影响机理

Bullegas等人[29]研究了应力集中在复合材料破坏时的作用，从中归纳出内应力的聚集与释放是材料破坏的主要原因. 而材料的破坏通常离不开应力的影响. 所以为更深入地讨论单晶硅的破坏机理，就有必要先观察和统计其在破坏前一刻的应力分布情况.

图4和图5利用0.5 ps内的应力张量数据将σx在晶体中的分布情况绘制为可视化图像. 从中可见，在单个点缺陷和点缺陷团簇中，σx的分布大致是相同的. 图5(a)说明应力的集中现象总是分布在垂直于外加应力方向的危险截面上. 在图5(b)，(c)和(d)中，随着应变的增加，与[100]方向成一定夹角的(111)晶面上的应力水平也逐步增加. 图4和图5已经从原子尺度上展示了σx的集中现象和点缺陷之间的联系，但仅凭应力云图还无法获取达成破坏的所需的极限内应力范围，故还需要对应力进行一些统计工作.

图6与图7定量地给出了所有单晶硅样本中破坏前1 ps内的最大应力与点缺陷大小的变化关系，图6中σx的变化情况表明，晶体内部的最大正应力是始终大于应力—应变图像上的屈服强度的，两者的变化趋势相同，说明单晶硅的破坏可能取决于该点的σx，但不能说是绝对. 另外σy和σz虽然被控制为接近0，但随着点缺陷大小的增大，却有着上升的趋势. 图7中的剪应力τxy，τyz和τxz则随着样品中点缺陷大小的增大而增大. 这一变化情况则表明，含有点缺陷的单晶硅晶体在单向拉伸时其内部存在着复杂应力状态，仅单纯考虑σx的变化情况与屈服强度的关系过于片面，晶体的破坏应该是各向应力共同作用的结果. 为讨论正应力和剪应力变化情况之间的联系以及此类复杂应力状态的决定因素，还需从强度理论着手；使用分子动力学模拟得到的应力张量数据，对晶体内每个原子的受力情况按各强度理论的公式进行校核，确定一种能够较好地描述单晶硅屈服强度降低现象的强度理论，以供后续的分析与讨论.

图6 破坏前1 ps内所有晶体样本内最大正应力与点缺陷大小的关系

图7 破坏前1 ps内所有晶体样本内最大剪应力与点缺陷大小的关系

图8 按四类强度理论给出的理论值与单晶硅实际屈服强度的校核

(4)

式中，σs为米塞斯强度理论给出的理论屈服强度，σx，σy和σz分别指某点X，Y和Z轴上的正应力，τxy，τyz和τxz分别对应某点XY，YZ和XZ平面上的切应力. 当材料内部某一点的正应力和切应力满足米塞斯等效应力公式给出的理论屈服强度时，该点即进入塑性状态. 因为已经发现米塞斯强度理论可以很好地描述单晶硅的屈服强度的变化，所以为进一步讨论点缺陷对单晶硅屈服强度的影响机理，就有必要对晶体内部的米塞斯应力水平进行分析.

图9 样本3中面在发生破坏时的裂纹和米塞斯应力的演化情况，(a) 时间步2042500的裂纹快照；(b) 时间步2044500的裂纹快照；(c) 时间步2046000的裂纹快照；(d) 时间步2047500的裂纹快照；(a1) 时间步2042500的米塞斯应力快照；(b1) 时间步2044500的米塞斯应力快照；(c1) 时间步2046000的米塞斯应力快照；(d1) 时间步2047500的米塞斯应力快照

在了解了点缺陷对单晶硅力学性能影响的一系列现象后，将点缺陷对屈服强度的影响机理概括为以下结论：点缺陷会改变晶体内部的应力分布，在其周围区域产生应力集中效应，使这些区域的破坏条件符合米塞斯屈服准则. 因为应力集中区域的米塞斯应力总是最容易达到理论屈服强度，所以实际屈服强度会按米塞斯强度理论的预计趋势变化. 在讨论单晶硅的理论屈服极限时，可按米塞斯强度理论所给出的参考值来预测其力学性能.

3.3 破坏后的单晶硅微结构

3.2节中提到，小规模的破坏是在<111>晶面族上开始的，经过图9所示的延展后在晶体中产生宏观裂纹. 通过Ovito软件观察破坏后的结构时，发现晶体发生破坏的结构又表现出一定的规律性. 本文通过去除晶硅原子的方式，将破坏后的微结构几何构造绘制成图10.

图10(a)与图10(c)中的破坏微结构说明，于[100]方向施加恒定应变率时，裂纹在(100)面与(110)面上表现出明显的周期性重复. 而从图10(b)中的[110]晶向视角来看，发生破坏后的结构类似于一种非正交二维网格. 进一步观察后发现该网格是沿着解理面分布的. 另一方面，图10(d)中所给出的单晶胞破坏结构表明，二维网格在单个样品中的基本组成单元是两个夹角为70.53 °，且同属于<111>晶面族的晶面. 上述现象说明破坏不完全由拉力引起，剪力也参与了破坏，侧面验证了米塞斯强度理论对单晶硅力学性能描述的准确性.

图10 单晶硅晶体受[100]晶向拉力破坏后的微结构，(a) [100]晶向视角；(b) [110]晶向视角；(c) [010]晶向视角；(d) 单样品中的网格结构

表2 破坏后结构的面分布情况

4 结 论

本文通过分子动力学模拟研究了点缺陷对单晶硅力学性能的影响和机理，结果归纳为以下几点：

(1) 点缺陷对单晶硅弹性模量的影响很小，但会降低其屈服强度. 屈服强度σs与点缺陷大小c符合公式(3)中所述的指数函数关系.

(2) 通过对硅晶体中每个原子的应力进行计算和可视化处理，结果发现在点缺陷附近区域存在着应力集中现象，该现象说明点缺陷能够改变其附近区域的应力状态. 接着使用原子的六个应力张量对四大强度理论进行了校核，给出了理论屈服强度的变化情况. 其中，第四强度理论所给出的理论屈服强度最符合单晶硅实际屈服强度的变化情况，其描述的米塞斯应力的演化情况也能较好地解释单晶硅裂纹的扩展现象. 所以本文认为点缺陷对单晶硅屈服强度的影响机理是：点缺陷通过改变晶体内部的应力分布，使这些区域的发生破坏的应力条件符合米塞斯屈服准则. 因为米塞斯应力的最大值总是在应力的集中区域出现，所以实际屈服强度会按第四强度理论的理论预计趋势变化.

(3)通过直观地观察发生破坏后的微结构，发现破坏后的单晶硅结构类似于一种二维网格. 产生的微结构沿<111>解理晶面族分布，但在<111>晶面族中的4个可能发生破坏的晶面中，每次仅会有两个夹角为70.53°的晶面发生破坏，通过统计各晶面上的应力推测该现象的产生原因是由于破坏总倾向于在应力满足了破坏条件的晶面上产生.