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

硬度仅次于金刚石的α-碳化硅具有质轻、耐高温等优点，在摩擦器件、发动机涡轮叶片以及抗冲击涂层等领域得到了广泛的应用[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



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

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

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

V(r)=k(R-r)3,r<R,0,r≥R.

(1)

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

1.2 MD模拟计算1.2.1 基本原理和计算方法

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

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

Mα?2Rα?t2=fii=1,2,???,N

(2)

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

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

Fα=-E(RN)?Rα

(3)

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

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

vit+Δt2=1-ηvit-Δt2+Δtmifit

(4)

rit+Δt=ri(t)+vit+Δt2Δt

(5)

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

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

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

1.2.2 原子间的势函数

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

V=∑i<jVij(2)(rij)+∑i,j<kVjik(3)(rij,rik)

(6)

Vij(2)(r)=Hijrηij+ZiZjre-r/λ-Dijr4e-r/ξ-Wijr6

(7)

式中包括二体势和三体势 Hij和ηij为排斥力的强度和指数，Wij和Dij为范德瓦尔斯引力和电荷偶极子的强度，Zi表示以e为单位的有效电荷，ξ和λ分别代表电荷偶级和库伦作用的屏蔽长度，r=rij表示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



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

Fig.2The 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



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

Fig.3Structural 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应力计算可表述为

σMises=(σxx-σyy)2+(σyy-σzz)2+(σzz-σxx)2+6(τxy2+τyz2+τzx2)2

(8)

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

图4



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

Fig.4Distributions 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=12<110>方向上运动 从广义层错能曲线可以发现，在α-SiC中沿12<110>的shuffle-set位错的滑移具有最小能量，这是多晶α-碳化硅中出现这种类型位错的原因[20] 图6给出了位错环形成的“U形”位错环机制 可以看出，位错环受压痕应力区的影响向样本底部方向滑移 与10 K下的变形相对比可见，在900 K出现很多位错环，也是结构软化的原因

图5



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

Fig.5Structural 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



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

Fig.6Dislocation 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) 在压痕过程中，载荷位移曲线都呈现剧齿状趋势，波动现象随接触的作用而不断重复出现 这是器件接触表面间发生黏着接触失效的主要原因

参考文献

View Option 原文顺序文献年度倒序文中引用次数倒序被引期刊影响因子

[1]

Padture N P.

Advanced structural ceramics in aerospace propulsion

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

URLPMID [本文引用: 1]

[2]

Pan Q, Zhou H, Lu Q, et al.

History-independent cyclic response of nanotwinned metals

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

URLPMID [本文引用: 1]

[3]

Taloni A, Vodret M, Costantini G, et al.

Size effects on the fracture of microscale and nanoscale materials

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

[本文引用: 1]

[4]

Bourne N, Millett J, Pickup I.

Delayed failure in shocked α-silicon carbide

[j]. journal of applied physics, 1997, 81(9): 6019

[本文引用: 1]

[5]

Millett J, Bourne N, Dandekar D.

Delayed failure in a shock-loaded α-silicon carbide

[J]. Journal of applied physics, 2005, 97(11): 113513

[本文引用: 1]

[6]

Pickup I, Barker A.

Damage kinetics in α-silicon carbide; proceedings of the AIP Conference Proceedings, F, 1998

[C]. AIP

[本文引用: 1]

[7]

Plimpton S.

Fast Parallel Algorithms for Short-Range Molecular-Dynamics

[J]. J Comput Phys, 1995, 117(1): 1

[本文引用: 1]

[8]

Kelchner C L, J P S, C H J.

Dislocation nucleation and defect structure during surface indentation

[J]. Phys Rev B, 1998, 5811085

[本文引用: 1]

[9]

Branicio P S, Zhang J, Rino J P, et al.

Shock-induced microstructural response of mono-and nanocrystalline SiC ceramics

[J]. Journal of Applied Physics, 2018, 123(14): 145902

[本文引用: 1]

[10]

Zhang J, Branicio P S.

Molecular dynamics simulations of plane shock loading in SiC

[J]. Procedia Engineering, 2014, 75150

[本文引用: 1]

[11]

Makeev M A, Srivastava D.

Hypersonic velocity impact on a-SiC target: A diagram of damage characteristics via molecular dynamics simulations

[J]. Applied Physics Letters, 2008, 92(15): 151909

[本文引用: 1]

[12]

Makeev M A, Srivastava D.

Molecular dynamics simulations of hypersonic velocity impact protection properties of CNT/a-SiC composites

[J]. Composites Science and Technology, 2008, 68(12): 2451

[本文引用: 1]

[13]

Makeev M A, Sundaresh S, Srivastava D.

Shock-wave propagation through pristine a-SiC and carbon-nanotube-reinforced a-SiC matrix composites

[J]. Journal of Applied Physics, 2009, 106(1): 014311

[本文引用: 1]

[14]

Hirel P.

Atomsk: a tool for manipulating and converting atomic data files

[J]. Comput Phys Commun, 2015, 197212

DOIURLPMID class="outline_tb" " />

Dyson-Schwinger equations are important tools for non-perturbative analyses of quantum field theories. For example, they are very useful for investigations in quantum chromodynamics and related theories. However, sometimes progress is impeded by the complexity of the equations. Thus automating parts of the calculations will certainly be helpful in future investigations. In this article we present a framework for such an automation based on a C++ code that can deal with a large number of Green functions. Since also the creation of the expressions for the integrals of the Dyson-Schwinger equations needs to be automated, we defer this task to a Mathematica notebook. We illustrate the complete workflow with an example from Yang-Mills theory coupled to a fundamental scalar field that has been investigated recently. As a second example we calculate the propagators of pure Yang-Mills theory. Our code can serve as a basis for many further investigations where the equations are too complicated to tackle by hand. It also can easily be combined with DoFun, a program for the derivation of Dyson-Schwinger equations. PROGRAM SUMMARY: Program title: CrasyDSE Catalogue identifier: AEMY _v1_0 Program summary : http://cpc.cs.qub.ac.uk/summaries/AEMY_v1_0.html Program obtainable from: CPC Program Library, Queen's University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 49030 No. of bytes in distributed program, including test data, etc.: 303958 Distribution format: tar.gz Programming language: Mathematica 8 and higher, C++. Computer: All on which Mathematica and C++ are available. Operating system: All on which Mathematica and C++ are available (Windows, Unix, Mac OS). Classification: 11.1, 11.4, 11.5, 11.6. Nature of problem: Solve (large) systems of Dyson-Schwinger equations numerically. Solution method: Create C++ functions in Mathematica to be used for the numeric code in C++. This code uses structures to handle large numbers of Green functions. Unusual features: Provides a tool to convert Mathematica expressions into C++ expressions including conversion of function names. Running time: Depending on the complexity of the investigated system solving the equations numerically can take seconds on a desktop PC to hours on a cluster.

[15]

Xiang H, Li H, Fu T, et al.

Formation of prismatic loops in AlN and GaN under nanoindentation

[J]. Acta Materialia, 2017, 138131-9

PMID

[16]

Hoover W G.

Constant-pressure equations of motion

[J]. Physical Review A, 1986, 34(3): 2499

[17]

Vashishta P, Kalia R K, Rino J P, et al.

Interaction potential for SiO 2: a molecular-dynamics study of structural correlations

[J]. Phys Rev B, 1990, 41(17): 12197

[18]

Vashishta P, Kalia R K, Nakano A, et al.

Interaction potential for α-silicon carbide: a molecular dynamics study of elastic constants and vibrational density of states for crystalline and amorphous α-silicon carbide

[J]. Journal of applied physics, 2007, 101(10): 103515

[19]

Maras E, Trushin O, Stukowski A, et al.

Global transition path search for dislocation formation in Ge on Si (001)

[J]. Comput Phys Commun, 2016, 20513

[20]

Sun S, Peng X, Xiang H, et al.

Molecular dynamics simulation in single crystal 3C-SiC under nanoindentation: Formation of prismatic loops

[J]. Ceramics International, 2017, 43(18): 16313

Advanced structural ceramics in aerospace propulsion

1

2016