基于分子动力学模拟的纳米多晶α-碳化硅变形机制

doi:10.11901/1005.3093.2019.600

基于分子动力学模拟的纳米多晶α-碳化硅变形机制

施渊吉¹, 陈显冰¹, 吴修娟¹, 王红军^,¹, 郭训忠², 黎军顽³

1 南京工业职业技术大学机械工程学院南京 210046

2 南京航天航空大学材料科学与技术学院南京 210016

3 上海大学材料科学与工程学院上海 200072

Deformation Mechanism of Nanoscale Polycrystalline α-Silicon Carbide Based on Molecular Dynamics Simulation

SHI Yuanji¹, CHEN Xianbing¹, WU Xiujuan¹, WANG Hongjun^,¹, GUO Xunzhong², LI Junwan³

1 School of Mechanical Engineering, Nanjing Vocational University of Industry Technology, Nanjing 210046, China

2 College of Material Science and Technology, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

3 School of Materials Science and Engineering, Shanghai University, Shanghai 200072, China

通讯作者: 王红军，教授，1991100119@niit.edu.cn，研究方向为精密加工

责任编辑: 黄青

收稿日期: 2019-12-25 修回日期: 2020-04-05 网络出版日期: 2020-08-25

基金资助:

江苏省高等学校自然科学研究面上项目.  19KJB430024
江苏省自然科学基金.  BK20181036
南工院自然科学培育基金.  YK190109

Corresponding authors: WANG Hongjun, Tel: (025)85864009, E-mail:1991100119@niit.edu.cn

Received: 2019-12-25 Revised: 2020-04-05 Online: 2020-08-25

Fund supported:

Natural Science Foundation of the Jiangsu Higher Education Institutions of China.  19KJB430024
Natural Science Foundation of Jiangsu Province.  BK20181036
National Science Foundation of NIIT.  YK190109

作者简介 About authors

施渊吉，男，1989年生，博士

摘要

在考虑晶界和温度效应影响的条件下，基于分子动力学法使用Vashishta势函数研究多晶α-碳化硅基体在纳米压痕作用下的塑性变形机制，分析载荷位移曲线并通过识别变形结构描述了变形区域中的原子破坏和迁移轨迹变化。在下压过程中，因接触载荷不断增大在接触区的晶粒内产生无定型化相变并不断向晶体内部扩展，扩展到晶界处被阻碍住。随着载荷的持续增大，晶界作为位错发射源在高应力水平下出现1/2〈110〉全位错滑移。同时，随着温度的升高α-碳化硅多晶的承载能力下降，特别是材料内部出现塑性变形，位错从晶界处形核长大并向晶体内部扩展，最后形成‘U型’位错环。

关键词： 无机非金属材料 ; 分子动力学 ; 多晶 ; α-碳化硅 ; 结构变形 ; 原子尺度

Abstract

Based on the molecular dynamics method, Vashishta potential function was used to study the plastic deformation mechanism of polycrystalline α-silicon carbide matrix under the action of nano indentation in terms of the effect of grain boundary and temperature. The load displacement curve was analyzed, and the change of atomic failure and migration path in the deformation area was described by identifying the deformation structure. As the contact load increased the amorphous phase transformation occurred in the contact zone and expanded to the crystal interior, which was blocked by the grain boundary. With the increase of the load, the grain boundary as the source of 1/2<110> perfect dislocation emission will slip at high stress level. In addition, with the increase of temperature the bearing capacity of SiC polycrystal decreases, especially the plastic deformation occurs inside the material, the dislocation grows from the grain boundary to the inside of the crystal, and finally forms 'U-shaped' dislocation ring.

Keywords： inorganic non-metallic materials ; molecular dynamics ; polycrystalline ; α-silicon carbide ; structural deformation ; atomic scale

PDF (16155KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

施渊吉, 陈显冰, 吴修娟, 王红军, 郭训忠, 黎军顽. 基于分子动力学模拟的纳米多晶α-碳化硅变形机制. 材料研究学报[J], 2020, 34(8): 628-634 DOI:10.11901/1005.3093.2019.600

SHI Yuanji, CHEN Xianbing, WU Xiujuan, WANG Hongjun, GUO Xunzhong, LI Junwan. Deformation Mechanism of Nanoscale Polycrystalline α-Silicon Carbide Based on Molecular Dynamics Simulation. Chinese Journal of Materials Research[J], 2020, 34(8): 628-634 DOI:10.11901/1005.3093.2019.600

硬度仅次于金刚石的α-碳化硅具有质轻、耐高温等优点，在摩擦器件、发动机涡轮叶片以及抗冲击涂层等领域得到了广泛的应用^[1]。α-碳化硅这些优异的宏观性能本质上是由其微观结构决定的^[2]。但是，微观结构影响宏观力学性能的机制尚不明确，从微观结构到宏观性能的跨尺度研究成为当前的热点和难点^[3]。

Bourne等^[4]和Millett等^[5]对三种等级的α-碳化硅进行了平板冲击和SHPB实验。他们发现，在平板冲击达到最大应力后受冲击的α-碳化硅出现延迟失效现象；而SHPB实验揭示了三种材料之间冲击响应的显著不同，即HEL和剪切强度随着制备工艺的不同而变化。Pickup等^[6]也对三种α-碳化硅样品进行了冲击试验，发现这些样品具有相似的密度和晶粒尺寸，但是材料制备工艺不同。Hopkinson压杆实验结果表明，制备方法不同不仅使材料抗冲击强度显著不同，还改变了材料失效的形态。他们还报道了这三种材料在平面冲击中的失效波现象。而实验法研究的局限性，在于很难揭示其内部微观结构和演化规律。相反，分子动力学法(MD)不仅能动态模拟和观察体系中的每个原子运动变状态的变化，还能观察材料内部的微细结构演变成为纳观接触学和实验法辅助研究中的有效计算法^[7,8]。

在过去的研究中，不少文献给出了单晶α-碳化硅陶瓷在高速载荷下的变形机理和力学性能。南加州大学的Branicio等^[9]用分子动力学模拟研究了刚性弹丸以15 km/s的超高速对α-碳化硅单晶板的撞击，发现撞击产生的巨大应力达到150 GPa，冲击波速最高可达24 km/s。模拟结果还发现，应力达到89 GPa以上可观察到α-碳化硅从闪锌矿结构到岩盐矿结构的相变。他们的模拟揭示了由单个位错核心参与的纳米延性诱导的裂纹成核机制。原子损伤机制，涉及从冲击引起的结构转变到塑性变形以及到脆性断裂的动态转变。压缩过程中的剪切应力高达45 GPa，压缩冲击波在自由表面反射后局部拉伸应力达到20 GPa。Zhang等^[10]采用分子动力学方法初步研究了α-碳化硅的 Hugoniot性质，发现α-碳化硅中冲击引起的塑性变形主要以变形孪晶的形式存在。Makeev等^[11]基于Tersoff势的分子动力学模拟，研究了在金刚石弹丸的高超音速速度冲击下非晶α-碳化硅靶的动态损伤响应。在他们的模拟中，在侵彻深度图中发现了四种不同的损伤状态:浅坑形成，深入靶板，深度穿透与靶板的局部熔化以及入射弹体的完全崩解。Makeev等^[12,13]也模拟了非晶α-碳化硅及其碳纳米管增强材料在冲击载荷下响应并进行了比较。上述文献主要研究了单晶材料，对多晶的微观结构尚不清楚。鉴于此，本文在考虑晶界和温度效应影响的条件下，基于分子动力学法使用Vashishta势函数研究多晶α-碳化硅基体在纳米压痕作用下的塑性变形机制，分析载荷位移曲线并通过识别变形结构描述变形区域的原子破坏和迁移轨迹变化。

1 分子动力学建模和计算

1.1 MD模型的建立

图1给出了多晶α-碳化硅的三维压痕模拟示意图。球体压头的半径为R，模拟设置盒子的几何尺寸为15 nm×15 nm×15 nm。所有的分子动力学压痕模拟都执行1 fs的恒定时间步长。用共轭梯度(CG)算法在压痕模拟前优化样本，使系统达到最小平衡能量从而处于稳定状态。在模拟中选择NVE系综，纳米压痕的加载速度设为40 m/s。使用Langevin恒温热浴控制系统的温度，分别在10 K、300 K、600 K和900 K进行压痕。为了在压痕过程中防止基底移动，样品底部的原子是固定的。z方向的边界条件设置固定边界条件，x、y方向的边界条件分别设置为周期边界条件。在压痕前将压头放置在上表面中心的上方0.2 nm处。在压入过程中压头向下移动压入样品。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 多晶α-碳化硅的纳米压痕分子动力学模拟模型

Fig.1 Simulation model of nanoindentation of polycrystalline α-silicon carbide (α-SiC)，the size of the simulation box is 15 nm×15 nm×15 nm

为了节省计算时间，将压头设为一个无原子的整体刚性球体/柱体。它排斥与其接触的所有原子，相互作用势为

V (r) = \{\begin{array}{l} k^{(R - r) 3}, r < R, \\ 0, r \geq R . \end{array}

(1)

其中r为压头中心与原子中心的距离，R为压头半径，k为压头刚度，设为1 eV·nm^-3[15]。

1.2 MD模拟计算

1.2.1 基本原理和计算方法

分子动力学模拟基于体系初始构型的位置和速度，按照牛顿第二定律计算调整其新的位置和速度。分子动力学模拟的系统通常为孤立的能量守恒系统。但是当有能量耗散到周围环境时，则采用非能量守恒法。

分子动力学模拟系统中原子的运动，遵循牛顿运动方程

M_{α} \frac{\partial^{2} R_{α}}{\partial t^{2}} = f_{i} i = 1, 2, \cdot \cdot \cdot, N

(2)

其中M_α为第α个原子的质量，R_α表示第α个原子的三维坐标，t为模拟时间，f_i为原子间的相互作用力。

确切描述纳尺度下的力学行为的关键，是确定原子间相互作用力f_i，即确定原子间相互作用势E(R^N)。因为

F_{α} = - \frac{E (R^{N})}{\partial R_{α}}

(3)

其中E(R^N)为系统的总势能，R^N为体系原子的三维坐标。

确定原子间的相互作用势后，需要求得原子的运动轨迹，可用有限差分法求解运动方程，如Verlet算法、Verlet算法的速度形式、Leap-frog算法、Beeman算法和Gear算法等。分子动力学模拟中常用的是蛙跳算法(Leapfrog Algorithm):

v_{i} (t + \frac{Δ t}{2}) = (1 - η) v_{i} (t - \frac{Δ t}{2}) + \frac{Δ t}{m_{i}} f_{i} (t)

(4)

r_{i} (t + Δ t) = r_{i} (t) + v_{i} (t + \frac{Δ t}{2}) Δ t

(5)

其中v_i表示第i个原子的速度，m_i表示第i个原子的质量，η表示原子运动时所受的阻尼大小。由原子的瞬时位置r^N(t)可求出原子间的作用力f^N(t)，进而求出系统原子的速度v^N(t+∆t/2)和更新位置r^N(t+∆t)。

在正则系综(NVT，即系统的原子数、体积和温度恒定不变)下，由于原子间的能量传输使系统温度升高，而能量必须传输到外界以维持整个体系的平衡。由此可见，控制系统的温度非常重要。温度也反应了系统中的耗散动力学机制和能量交换。

温度对系统的平衡着很大的影响，体系的动能也与温度有关，因此在进行分子动力学模拟时控温很重要。除了Berendsen热浴，Gaussian热浴和Nose-Hoover^[16]热浴也是分子动力学模拟中广泛应用的控温手段。分子动力学模拟结果表明，虽然控温方法不同但是其计算结果差异较小。

1.2.2 原子间的势函数

势函数描述系统中原子与原子之间的相互作用势，是决定计算所需时间和仿真结果准确性的核心因素。本文采用SW势的改进势Vashishta势。Vashishta势基于SW的改进势^[17]，已经成功地用于描述半导体和陶瓷材料的行为。本文用它描述Si-C系统的相互作用^[18]。选用这种势函数的主要原因，是其可预测SiC的各种相的性质。在SW势的基础上，Vashishta势将排斥力、屏蔽库仑力、屏蔽电荷偶极子、弥散相互作用与键角能相结合。其表达式为

V = \sum_{i < j} V_{i j}^{(2)} ({r_{i}}_{j}) + \sum_{i, j < k} V_{j i k}^{(3)} (r_{i j}, r_{i k})

(6)

V_{i j}^{(2)} (r) = \frac{H_{i j}}{r^{η_{i j}}} + \frac{Z_{i} Z_{j}}{r} e^{- r / λ} - \frac{D^{i j}}{r^{4}} e^{- r / ξ} - \frac{W_{i j}}{r^{6}}

(7)

式中包括二体势和三体势。H_ij和 $η_{i j}$ 为排斥力的强度和指数，W_ij和D_ij为范德瓦尔斯引力和电荷偶极子的强度，Z_i表示以 $|e|$ 为单位的有效电荷，ξ和λ分别代表电荷偶级和库伦作用的屏蔽长度，r=r_ij表示i和j原子之间的距离。

1.2.3 局部结构变形的识别方法

为了识别结构变形、缺陷的成核和演化，本文使用金刚石结构识别法(IDS)^[19]。IDS可用于识别六方或立方金刚石晶格中的原子排布。为了将中心原子类型分类，考虑了第二近邻原子来区分六方和立方金刚石结构。这种方法涉及到第二临近原子的几何排列特征，还可用来可视化识别堆垛层错(SFs)和位错。为了更清楚的演示缺陷区域，用不同的颜色表示不同的晶格结构，六方金刚石结构(第一邻近)用红色表示，六方金刚石结构(第二邻近)用橘色表示，立方金刚石结构用蓝色表示。

2 结果和分析

2.1 多晶α-碳化硅变形机制分析

不同温度下压痕的P-h曲线，如图2所示。在初始加载阶段弹性变形很明显，例如，10K的样本P-h曲线在h=1.5 nm前没出现明显的载荷跌落，整体表现为弹性变形。对于300K，600K和900K样本的P-h曲线，与塑性变形机制相关的初始跌落点在曲线上很明显，且落差较大。对比不同温度h=3.0 nm下的载荷可以发现，随着温度的升高载荷不断降低，多晶α-碳化硅材料表现出高温软化。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 多晶α-碳化硅的纳米压痕力-压痕深度曲线

Fig.2 The relationship between contact force (P) and displacement ( $δ$ ), (a) - (d) load-displacement curves of 10K, 300K, 600k and 900K are given respectively

为了更深入地研究多晶α-碳化硅压痕载荷下的变形机制，分别将原子在h=1.0、1.4、1.8、2.2、2.6和3.0 nm处的IDS分布提取出来，如图3所示。由图3可以看出，随着压痕深度的增大压痕变形的主要机制是非晶化变形。因接触载荷不断增加在接触区的晶粒内产生非晶化相变并不断向晶体内部扩展, 当增加到晶界处被晶界阻碍住。随着载荷的持续增加，晶界作为位错发射源在高应力水平下出现位错滑移，如图中红色箭头所示。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 压痕过程中基体的变形

Fig.3 Structural deformation during indentation (a)~(f) the distribution of IDS at h=1.0, 1.4, 1.8, 2.2, 2.6 and 3.0 nm

材料的结构变形，是在应力值超过一个临界值后发生的。因为结构失效破坏与von Mises应力相关，有必要计算von Mises应力以确定结构变形在压痕过程中的相关性。von Mises应力计算可表述为

σ_{M i s e s} = \sqrt[]{\frac{(σ_{x x} - σ_{y y})^{2} + (σ_{y y} - σ_{z z})^{2} + (σ_{z z} - σ_{x x})^{2} + 6 ({τ_{x y}}^{2} + {τ_{y z}}^{2} + {τ_{z x}}^{2})}{2}}

(8)

其中σ_xx、σ_yy、σ_zz、τ_xy、τ_yz和τ_zx表示应力张量的各分量，可根据维里定理计算。从图4可以看出，在晶界处和压头下方，von Mises应力值最大。这意味着，相变和位错等结构变形容易在压头和晶界处发生。但是，由于压头下方的应力集中区较大且随着压头下压应力值的持续增加，压头下方的区域更容易发生结构变形。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 在压痕过程中的von Mises应力分布

Fig.4 Distributions of von Mises stress of different depths under nanoindentation (a)~(f) the von Mises stress distribution at h =1.0, 1.4, 1.8, 2.2, 2.6 and 3.0 nm

2.2 温度效应对多晶α-碳化硅变形机制的影响

从图2可以看出，在h=3.0 nm处10 K，300 K，600 K和900 K对应的载荷分别为925 nN, 875 nN, 825 nN和800 nN。显然，温度的升高导致材料的承载能力的下降。为了更清楚地分析在压痕过程中的温度效应的影响，图5和图6给出了900 K下的变形过程。对比图5和图3给出的原子结构变形信息，可以发现在900K和10K下的变形机制相同。压痕变形的主要机制是非晶化变形，随着载荷增加出现位错结构。因接触载荷不断增加，在接触区的晶粒内产生非晶化相变并不断向晶体内部扩展，当增加到晶界处被晶界阻碍住。随着载荷的持续增加，晶界作为位错发射源在高应力水平下出现位错滑移，如图中红色箭头所示。将未变形结构删掉，从图6可见，在900K位错核心从压头下和晶界处萌生，并不断向晶体内部发射。压痕产生了钉扎在滑移面基面{111}面底部的“U形”半环，这种半环就是所谓的“U形”位错环。这些半环绝大部分由三个部分组成:一部分平行于基面，另外两部分垂直于基面，都沿着Burgers矢量b= $\frac{1}{2} < 110 >$ 方向上运动。从广义层错能曲线可以发现，在α-SiC中沿 $\frac{1}{2} < 110 >$ 的shuffle-set位错的滑移具有最小能量，这是多晶α-碳化硅中出现这种类型位错的原因^[20]。图6给出了位错环形成的“U形”位错环机制。可以看出，位错环受压痕应力区的影响向样本底部方向滑移。与10 K下的变形相对比可见，在900 K出现很多位错环，也是结构软化的原因。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 在900K压痕过程中基体的变形

Fig.5 Structural deformation during indentation at 900 K, (a)~(f) the distribution of IDS at h=1.0, 1.4, 1.8, 2.2, 2.6 and 3.0 nm

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 在900K的压痕过程中位错的演化

Fig.6 Dislocation evolution during indentation at 900K (a)~(f) the distribution of atomics structures at h=1.0, 1.4, 1.8, 2.2, 2.6 and 3.0 nm

3 结论

(1) 在下压过程中，因接触载荷不断增大α-碳化硅多晶接触区的晶粒内产生无定形化相变并不断向晶体内部扩展, 在晶界处被晶界阻碍住。随着载荷的持续增大，晶界作为位错发射源在高应力水平下出现位错滑移。

(2) 随着温度的升高α-碳化硅多晶的承载能力下降，特别是材料内部出现塑性变形。在高温下位错更容易从晶界处形核长大并逐渐向晶体内部扩展，最后形成“U形”位错环。这是高温软化的主要原因。

(3) 在压痕过程中，载荷位移曲线都呈现剧齿状趋势，波动现象随接触的作用而不断重复出现。这是器件接触表面间发生黏着接触失效的主要原因。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Padture

N P

Advanced structural ceramics in aerospace propulsion

[J]. Nature materials, 2016, 15(8): 804

URL PMID [本文引用: 1]

[2]

Pan

, Zhou

, Lu

, et al.

History-independent cyclic response of nanotwinned metals

[J]. Nature, 2017, 551(7679): 214

URL PMID [本文引用: 1]

[3]

Taloni

, Vodret

, Costantini

, et al.

Size effects on the fracture of microscale and nanoscale materials

[J]. Nature Reviews Materials, 2018, 3(7): 211