Study on shock compression phase transition of single crystal siliconbased on molecular dynamics simulation
-
摘要: 晶体硅具有复杂的相变机制，在相图研究中受到广泛关注，其在动载荷下的变形机制是当前研究热点。为揭示晶体硅在强动加载下的变形和相变行为特征，基于分子动力学方法，采用平板冲击加载方式，模拟研究了单晶硅在初始环境温度为300 K时分别沿[001]、[110]和[111]晶向的不同强度下的冲击压缩行为，冲击粒子速度为0.3～3.2 km/s。研究发现，随着冲击粒子速度的增加，单晶硅剪切应力在逐渐增加后由于结构相变发生急剧下降，相变阈值和相变机制均呈现各向异性。其中，沿[001]晶向冲击压缩下观察到多种固-固相变以及固-液相变，并观察到与最新文献的实验高度一致的固-液共存现象。研究结果可为动加载下晶体硅的相变研究提供纳米尺度的结果支撑。Abstract: Crystalline silicon has a complicated phase transition mechanism, which has received extensive attention in the research field of phase diagram, and the deformation mechanism of silicon crystals under dynamic loading is the current research hotspot. In order to reveal its deformation and phase transition behaviors under intensive dynamic loading, molecular dynamics method was used to simulate the shock compression behavior of single crystal silicon along the crystal directions [001], [110] and [111] at an initial ambient temperature of 300 K, respectively. All simulations were carried out basing on the classical open-source codes LAMMPS and a Tersoff interatomic potential was adopted to describe the material responses of silicon under dynamic compression. Before shock loading, periodic boundary conditions were applied along the three independent directions, and an NPT ensemble was used to equilibrate the systems; then shock compression was applied by using the piston method, where a virtual piston wall impinges the sample such that the particle velocity in the sample is the same as the piston speed after the shock reaches a steady state. The shock particle velocities varied from 0.3 km/s to 3.2 km/s, and a timestep of 0.001 ps was adopted. During the stress wave formation and propagation, the simulation system was in the NVE ensemble with the absence of temperature control. The loading method and effect are similar to typical plane impact experiments. The results show that with the increase of shock particle velocity, the shear stress of single crystal silicon increases gradually and then decreases sharply due to the structural phase change. Both the phase transition threshold and the phase transition mechanism are anisotropic. Among them, a variety of solid-solid phase transitions and solid-liquid phase transitions are observed under shock compression along the [001] crystal direction. The phenomenon of solid-liquid coexistence is highly consistent with the recent international experiments. The research results provides new nano-scale results to support the study of phase transition of crystalline silicon under dynamic loading.
-
晶体硅具有纯度高、晶体取向良好等特点，在电子等领域应用广泛[1]，同时由于其复杂的结构相变机制，在相图研究中也被广泛关注[2-9]。
在准静态荷载下，单晶硅表现出显著的多态性[2-9]，目前已观察到13种相结构。Mujica 等[10]指出，单晶硅在10～11 GPa的载荷下发生了从立方金刚石（cubic diamond, CD）结构到β-白锡（β-tin）结构的相变，并通过第一性原理计算出存在中间暂稳态 Imma相。而增加压力时，晶体硅继续发生向简单六角形（simple hexagonal，SH）结构的转变；再进一步增加压力，变为Cmca空间群结构，而在41 GPa的高压下，晶体硅又从Cmca结构转变为六方密堆积结构（hexagonal close pack，HCP）。在动态荷载实验中，Gilev等[11]研究了在冲击压缩下单晶硅的金属化现象，发现在冲击压力达到23 GPa时，晶体结构的缺陷导致电导率发生了偏离。Loveridge-Smith等[12]用X射线衍射技术测量了冲击加载下单晶硅的正交平面晶格参数，发现沿压缩方向单晶硅的晶格间距减少了近11%，而与压缩方向正交平面的晶格间距没有显著变化。Turneaure 等[13]研究了单晶硅[100]、 [111]晶向在15.9、21.7 GPa冲击压力下的动力学响应，发现其弹性波和非弹性波晶向依赖性很大，相变波晶向依赖性很小。Zhao等[14]对单晶硅[001]晶向进行了激光冲击压缩和回收实验，观察到两个塑性响应：靠近冲击表面的大块非晶层和沿{111}平面产生的滑移带，该结果揭示了单晶硅在一维应变条件下的塑性响应。Smith 等[15]对单晶硅[001]、[110]和[111]晶向进行了冲击实验，揭示了单晶硅材料冲击响应的各向异性，发现了其在冲击压力高于13 GPa时转变为简单六角形结构。Liu等[16]研究了单晶硅沿[100]晶向冲击压缩的塑性变形，发现其塑性变形能力随着温度增加而增加。Kishimura等[17]对在38 GPa冲击压力下的单晶硅进行回收并进行了X射线分析和拉曼光谱分析。XRD分析显示，单晶硅在11 GPa冲击压力下有少量亚稳态相产生，在38 GPa冲击压力下有Cu3Si这一物质生成，而在其余加载压力下并没有明显的高压相和亚稳相存在。拉曼光谱的分析显示，其峰值的整体偏移是由晶体尺寸减少引起的，该结果说明了单晶硅冲击压缩下的晶粒化过程。Turneaure等[18]用实时XRD检测了单晶硅在冲击压力低于54 GPa下的晶体结构，发现单晶硅在冲击压力为31～33 GPa时产生冲击熔化，并从熔融边界再结晶到密堆积六角形结构。由于硅和锗属于同族元素，具有相似的物理力学特性，因此Renganathan 等[19]用原位XRD检测了锗在[100]晶向的相变，发现锗在冲击压力为15.7 GPa时转变为白锡结构，在冲击压力为31.5 GPa时转变为熔融态；应力卸载后，又转变为立方金刚石结构，该结果表明锗从白锡结构到立方金刚石结构的相变是可逆的。McBride等[20]用实时X射线检测技术揭示了多晶硅的冲击相变行为，观察到不同压力下的多种固-固相变以及固-液相变。Paul 等[21]则利用密度泛函理论、进化算法和晶格动力学等多种方法，构建了4 TPa和6000 K下硅的相图，为极高压相图的研究提供了参考。
尽管以上实验与计算研究结果丰富，但实验之间也存在较显著的差异。一方面，实验直接获取的物理量过少，难以直接表征结果；另一方面，各个实验设备和实验过程也存在差异，变量如应变率、加载时间等不一致，而这些因素都会对单晶硅的相变有显著影响[22]。同时，单晶硅作为一种脆性材料，不仅加工困难，在冲击压缩时样品也很难回收，而且单晶硅的相图比较复杂，存在稳态相和非稳态中间相等多个相，因此它在冲击压缩下的物理和力学特性尚无完全统一的结论。此外，单晶硅冲击实验中加载时间较短，很难实时测量其相变过程中的晶格结构。所以，尽管检测微观结构变化的原位实时诊断技术取得了很大的进展，但对于可逆相变往往还是求助于把宏观力学量的实验测试结果与静高压数据进行比较，以确认新相的结构、相变的类型及性质，这也会导致研究结果的差异。
由于分子动力学具有超高的时间和空间分辨率，在微纳米尺度研究材料的物理和力学特性时具有独特的优势。而在当前的超算能力下，分子动力学模拟规模可以达到千万甚至数十亿原子，并能同时获取各个原子运动的轨迹，进而从原子尺度到微纳米尺度揭示材料的变形与力学参量的关系，使得其能够对理论计算和实验研究提供重要的补充。在分子动力学模拟单晶硅方面也有诸多发现。Swift 等[23]用第一性原理计算得到单晶硅的物态方程，这个物态方程可以反应单晶硅从金刚石结构到体心四方结构的转变。Demkowicz 等[24]用SW势函数模拟非晶硅在施加恒定压力下的变形情况，发现压缩过程中存在类晶体和类液体部分，后者更容易发生塑性流动。Kumagai等[25]改进了Tersoff 势函数，使之既能模拟弹性常数又能准确估测熔点。Higginbotham[26]等通过分子动力学模拟，表明压力诱导相变是弹性波异常的原因。Mogni 等[27]则基于Tersoff-like 势函数模拟[001]晶向单晶硅在激光冲击下的相变，发现了新的Imma相。Zhao 等[14]研究了基于分子动力学模拟的冲击压缩下单晶硅的力学响应，发现硅在冲击压力低于10 GPa时为弹性响应，在稍高压力下，由于位错和缺陷沿着特定的晶体学方向扩展，会进一步产生新的缺陷，此外在更高压力下会产生交叉的非晶带。
尽管基于分子动力学的单晶硅变形与相变的模拟研究成果丰富，但仍有若干关键问题尚未明确，包括单晶硅的动态压缩各向异性，以及单晶硅复杂相变机制的物理描述和不同因素对其相变的影响规律和机制等。深入研究这些关键问题有利于更加全面地认清单晶硅在动态压缩下的变形和相变的物理和力学特性。本文中将基于大规模分子动力学方法模拟研究单晶硅的冲击压缩响应，从微纳米尺度认识晶体硅的冲击压缩特性，揭示其在动加载下的各向异性和结构相变等特征。
1. 研究方法
1.1 模型的构建与冲击压缩模拟方法
基于广泛使用的LAMMPS[28]平台，采用分子动力学方法模拟研究单晶硅的冲击压缩响应，主要对其非弹性变形和相变行为展开分析。晶体硅通常情况下具有立方金刚石结构，在300 K温度下晶格常数为0.543 nm。模型的建立是由单个晶胞沿x、y、z晶向分别复制30、30和390个晶胞得到的超胞结构。因此，模型尺寸约为16.3 nm×16.3 nm×212 nm，见表1。
表 1 单晶硅计算模型详细参数Table 1. Parameters of single crystal Si sample for MD simulation加载晶向 x轴 y轴 z轴 模型原子数 晶向 模型尺寸/nm 晶向 模型尺寸/nm 晶向 模型尺寸/nm [001] [100] 16.3 [010] 16.2 [001] 217.0 ～2.84×106 [110] [ˉ110] 16.3 [001] 16.2 [110] 217.0 ～2.84×106 [111] [ˉ1ˉ12] 16.0 [1ˉ10] 16.2 [111] 214.8 ～2.76×106 晶体硅中Si-Si原子间的相互作用采用改进的半经验Tersoff势函数描述，并使用Erhart等[29]优化的参数。与改进前的势函数相比，该势函数可以成功模拟晶体硅的多种性质，包括弹性常数和一般热力学性质，以及在静水压缩和冲击压缩条件下立方金刚石相到β-白锡相的结构特征[27]。
在施加冲击压缩前，先对构建的初始模型进行弛豫，弛豫过程中保持三维方向均为周期性边界条件，使整个晶体模型在NPT系综下弛豫了20 ps。从能量和残余应力的曲线上确认了材料模型为能量趋于稳定且残余应力在误差允许范围内趋于零的状态，此时获得了合理稳定的初始结构。在冲击压缩加载阶段，为了消除边界对应力波的影响，在沿z轴施加自由边界条件的同时，沿x、y方向施加周期性边界条件。模拟平板冲击采用的是活塞法，活塞对自由表面产生冲击波。单晶硅冲击加载的示意图如图1所示。冲击压缩模拟过程中系统处于NVE系综，并且不施加温度控制，时间步长为1 fs。这样的加载方式和效果，与典型的平面冲击实验相似，均为一维应变问题。为揭示单晶硅冲击响应的各向异性，尤其是冲击结构的演化，分别对[001]、[110]和[111]晶向进行系列冲击压缩研究，冲击粒子速度从0.3 km/s逐渐增加至3.2 km/s。
图 1 分子动力学模拟中单晶硅的冲击压缩示意图Figure 1. Schematic of shock compression of single crystal Si in molecular dynamic simulations1.2 主要物理量的计算和可视化分析方法
应力的计算采用Virial定律[30]，对原子i，其应力张量分量
σij(i,j∈{x,y,z}) 为：σij=−1Ωk[mkvk,ivk,j+12Np∑n=1(r1iF1j+r2iF2j)+13Na∑n=1(r1iF1j+r2iF2j+r3iF3j)] (1) 式中：
Ωk 、mk 和vk 分别为原子k的体积、质量和速度。选取应力张量
σ33 （或σzz ）作为Hugoniot应力，同时定义剪切应力如下：τ=12(σ33−12(σ11+σ22)) (2) 为了深入分析单晶硅的冲击压缩响应，利用OVITO[31]进行可视化分析，并采用配位数、中心对称参数、键角分析等多种分析技术，对晶体硅的变形和结构变化进行识别和分析。
2. 结果与讨论
2.1 冲击响应与各向异性
分子动力学模拟计算结果与不同文献中实验所得的晶体硅冲击Hugoniot应力和体积变化的关系如图2所示。从图2中可以看出，在体积压缩较小时，所有结果高度吻合；压缩较大时，实验的应力-体积关系呈现两大分歧：一是以Goto等[32]和Gust等[33]的结果为代表，这两个独立研究团队的结果比较接近，二是以Turmeaure等[13]和McBride等[20]的实验结果为代表。因此，对晶体硅不同冲击实验研究的差异值得关注，这可能与实验手段、技术和材料样品差异有关。一方面，Gust等[33]早期采用爆炸驱动的冲击波压缩实验，而Goto等[32]和Turmeaure等[13]采用轻气炮驱动的平板冲击装置对不同晶向的单晶硅进行冲击实验，McBride等[20]则采用激光驱动的应力脉冲加载对多晶硅进行冲击压缩实验；另一方面，尽管所有实验都是运用Rankine-Hugoniot关系推算试样的内部压力，但压力推算中所需的粒子速度或相变波速度测定的差异显著[15]，例如Gust等[33]提及其结果与相关历史数据的对比时，对应于HEL和相变点的临界体积（或密度）是很吻合的，但压力数值却存在差异，他们将此归因于粒子速度的测定差异。值得注意的是，两大存在分歧的压力-体积关系均有相互独立的研究团队的结果支撑，可能预示着深层次的原因在于速度测定的技术差异与硅材料响应有关。Goto等[32]的轻气炮驱动的平面冲击实验中采用了与Gust等[33]相似的测量技术（斜镜与条纹相机）对硅试件的自由面测速，其Hugoniot压力推算值也与Gust等[33]非常吻合；然而，Turneaure等[13]的轻气炮驱动平面冲击实验与McBride等[20]的激光冲击实验则采用了VISAR测试系统对硅/LiF窗口或者硅自由面速度测定，其结果相接近。此外，Smith等[15]指出硅自由表面端存在的应力松弛和初始的速度溅射下的自由表面测量结果可能会影响应力和体积的关系测定。因此，这些结果在压力数值上的差异深刻反映了实验技术差异和数据处理方法对实验结果的重要影响。
与实验结果相比，本文计算结果与Goto等[32]和Gust等[33]的结果更接近，但也在一定程度上高于实验值。尽管如此，基于该原子间势函数模拟计算的应力-体积整体变化趋势与实验结果相似。例如，随着体积压缩的增加，可以明显看出在体积比低于0.8时单晶硅已经处于非弹性阶段；而进一步压缩至体积比为约0.67时，材料压力急剧增加显示出典型的高压相变特征。因此，尽管该原子间势函数对单晶硅的Hugoniot压力存在高估的不足，但其可以成功反映晶体硅在冲击压缩下的若干重要变形响应特征，对本文研究动高压下硅的变形和相变仍有价值。
在不同冲击压缩下，单晶硅的剪切应力随粒子速度变化如图3所示。在[001]、[110]和[111]晶向中，其剪切应力在一定范围内均随着粒子速度的增加而逐渐增加，达到一定临界粒子速度后开始急剧下降。但是，该临界粒子速度在不同晶向差异显著，在[001]晶向最小，为1.4 km/s，在[110]晶向最大，为1.8 km/s，而在[111]晶向时介于两者之间，为1.5 km/s。在这3个晶向中，最大剪切应力分别为5.0、7.8和6.2 GPa。剪切应力的大幅度下降意味着单晶硅在这些冲击加载下发生了重要的变形（如塑性或相变等），使得剪切应力释放。可以观察到 [001]和[111]晶向中剪切应力经历了由正值转为负值的阶段：分别在1.8、2.4 km/s达到最大负剪切应力，而后逐渐减小。负剪切应力现象存在于单晶碳化硅[34-35]和单晶金属铁中[36]，学者们普遍认为这是由于材料过度弛豫释放的结果，并且具有晶向依赖性，多发生在[001]晶向。从剪切应力的计算方法角度分析，负剪应力意味着在冲击加载中材料内部的横向应力高于纵向应力。最终，[001]、[110]和[111]晶向的剪切应力分别在2.0、1.9、2.6 km/s之后趋近零，因此材料失去剪切强度。此外，结合图1中应力和体积压缩的关系，单晶硅在[110]晶向的奇异表现还体现在弹塑性或相变转变的临界状态，相比于[001]和[111]晶向的在约0.85体积压缩比时的“弹塑性”转变，[110]晶向则表现出奇异的弹性特征。这种冲击压缩下奇异的弹性响应现象在金刚石[37]等材料的分子动力学模拟研究也有过发现，这些材料呈现出非单调的剪应力-应变关系。由于剪切应力是塑性等变形的驱动力，因此，单晶硅冲击压缩下奇异的弹性响应应该与其剪切应力的非单调性相关[38]。
图 3 冲击Hugoniot状态的剪切应力(τ)与粒子速度(up) 曲线Figure 3. Shock Hugoniot shear stress as a function of particle velocity为了进一步探究单晶硅冲击加载下的材料特性及其各向异性，图4中给出了沿[001]、[110]和[111]晶向加载时单晶硅的纵向应力、剪切应力和密度波剖面。根据单晶硅的冲击Hugoniot曲线，单晶硅3个晶向的弹塑性或相变阈值是不同的，[110]是3个晶向相变阈值最大的，阈值粒子速度为1.8 km/s，其次是[111]晶向，阈值粒子速度为1.5 km/s，[001]晶向最小，为1.4 km/s。因此，波剖面分析中选取了各晶向中具有代表性的3个粒子速度状态。如图4(a)所示，在[001]晶向中，粒子速度分别为1.3、1.4、1.5 km/s的波剖面体现了单晶硅从单一弹性波到弹性-塑性/相变波的冲击波结构转变过程。在纵向应力波剖面中，当粒子速度为1.3 km/s时，波剖面平整，且波阵面十分狭窄陡峭。当粒子速度增加至1.4 km/s时，弹性前驱波幅值增加，但同时材料出现软化现象，导致弹性波之后出现应力幅值下降的现象，同时由于塑性或相变的激发，在靠近活塞端的波剖面区域呈现比软化区域更高幅值的塑性或相变波。当粒子速度为1.5 km/s时，类似的现象依然存在，但由于冲击强度增加，弹性前驱波幅值增加，这可能是由于单晶硅冲击Hugoniot弹性极限对加载速率敏感性很强导致的；同时，更高的冲击强度使单晶硅的软化越发明显，使得弹性前驱波之后的波幅值下降更快；同样，塑性或相变波部分幅值增加，对应的波速度也在增加。从[001]晶向的剪切应力波剖面可以进一步揭示该现象的变化。剪切应力对材料的响应非常灵敏，单晶硅在冲击压缩下呈现复杂的材料响应，包括弹性波及其软化和塑性或相变波。剪切应力随着弹性波的出现急剧上升，然后随着软化作用而略有降低，并在塑性或相变波区域迅速下降，更强的塑性或相变波引起更大幅值的剪切应力释放。密度波剖面的变化对应了包括剪切应力波剖面在内的各物理量波剖面情况，即密度随着弹性波抵达而快速上升，同时由于软化作用略有回落，并在相变波区域由于材料的进一步压缩而显著上升。本文中研究获得的单晶硅的波剖面特征在Higginbotham等[26]和Stubley等[39]的研究工作也得到了体现。
图 4 单晶硅中分别沿[001]、[110]和[111]晶向的冲击波剖面Figure 4. Shock profiles in single crystal Si along [001], [110] and [111] crystal orientation在[110]和[111]晶向中，复杂的多波结构依旧存在，但软化作用更加微弱。同时，与[001]晶向不同的是，在[110]和[111]晶向中波剖面弹性波幅值几乎保持不变，显示出与加载速率无关，塑性或相变波更加复杂。
基于波剖面的分析，采用可视化分析技术进一步探究冲击引起的单晶硅非弹性变形和相变等材料响应。图5中分别给出了沿着[001]、[110]和[111]晶向冲击压缩时的可视化结果，粒子速度分别为1.4、1.8、1.5 km/s。结合上文分析可知，这些粒子速度冲击压缩下单晶硅3个晶向中的剪切应力均开始产生急剧的下降，意味着单晶硅在这些加载区域发生了重要变形响应。中心对称参数和配位数是反映晶体结构特征的重要指标，可以体现单晶硅的发生与结构相关的变化。图5(a)中单晶硅根据中心对称参数分析着色，而图5(b)则根据原子最小配位数着色。在中心对称参数分析中，3个晶向的单晶硅均发生了显著的中心对称参数变化，但呈现出不同的带状纹理。未干扰区域具有比压缩区域更高的中心对称参数，而在带状条纹区域该参数则更低，并且在交叉区域该参数呈现出趋于零的变化。与此同时，在原子最小配位数分析中也可以发现3个晶向的差异显著。[001]单晶硅存在配位数为5和6的带状和交叉区域；而[110]单晶硅中在靠近加载端一段区域几乎全部变成具有更大配位数的原子；在[111]单晶硅中则仅有极少数原子发生最小配位数的变化。这说明在3个晶向中单晶硅的晶体结构均发生了结构性的变化，但是呈现明显的晶体各向异性。从不同的配位数和中心对称参数看，单晶硅冲击压缩下发生了多种复杂结构相变，这体现在某些结构相变可引起中心参数的变化但不会引起配位数的变化。
图 5 单晶硅沿[001]、[110]、[111]晶向分别在1.4、1.8、1.5 km/s冲击压缩下的可视化Figure 5. Visualizations of single crystal Si in [001], [110] and [111] crystal orientations with shock particle velocities of 1.4, 1.8, 1.5 km/s2.2 沿[001]晶向冲击的变形与相变
为进一步探究单晶硅的冲击结构相变等特征，重点分析[001]晶向单晶硅的情况。图6中给出了不同粒子速度下的单晶硅冲击响应可视化分析，通过晶格原子的偏转角可以显著地区分出冲击压缩下单晶硅变形和相变的局部变化，蓝色区域原子的偏转角非常小，黄色区域原子较初始结构偏转约54°，绿色区域原子较初始结构偏转45°。图中[001]晶向的相变带状区域延伸方向与加载方向呈57°夹角，与偏转角接近，这样的少许差异与材料处于压缩状态相关。这表明相变波的相变带遵循特定的晶体学取向，是沿着靠近{111}晶面向前传播，其中{111}晶面与{001}晶面夹角为54°。随着粒子速度的增加，相变带状区域宽度有略微增加的趋势。当粒子速度为2.0 km/s时，靠近冲击加载端的区域出现无序的状态，这个无序区域随着粒子速度的增加而扩展。当粒子速度为2.4 km/s时，一个显著的无序状态与晶体状态同时存在于相同加载方向位置。考虑到强冲击可能引起剧烈温升和更高的压力状态，因此可能引起单晶硅的熔化，这意味着在此冲击条件下单晶硅处于固-液共存的区域。值得注意的是，本文中所发现的冲击波阵面后方固-固相变和固-液相变共存的现象也出现在其他两个晶向，该发现与McBride等[20]的激光冲击实验诊断结果一致，本文中冲击压缩下固-液共存甚至完全熔化所对应的临界体积压缩比，也与McBride等[20]的实验高度吻合。如图6所示，当粒子速度高达2.8 km/s时熔化区域已经扩展很大，并在3.2 km/s时波阵面后几乎完全熔化。
图 6 不同冲击粒子速度下[001]单晶硅的结构演化Figure 6. Structural evolutions of single crystal Si in [001] orientation at various shock particle velocities基于图6的代表性局部区域，进一步通过键角分析和晶向分布函数分析冲击压缩下单晶硅的结构相变。图7(a)为几个代表性的局部区域的键角分析。区域 Ⅰ代表冲击波阵面处于带状区域之前的部分，区域Ⅱ代表黄色带状区域，区域Ⅲ代表绿色的交叉区域，区域Ⅳ代表中间合围区域，区域Ⅴ代表无定型区域。初始结构中，键角呈现109°，区域Ⅰ主要键角峰值在105°，区域Ⅱ和Ⅲ则依次主要呈现96.5°和接近90°，而区域Ⅳ则与初始结构非常相似。这说明除了区域Ⅳ，其余区域均出现了与结构变化相关的键角变化。图7(b)的径向分布函数进一步揭示了可能的结构变化。从图7可知，与初始结构相比，相变区域Ⅱ～Ⅲ、Ⅴ变化明显。其中区域Ⅱ～Ⅲ中的峰值位置和数量与初始结构差异显著，体现了晶体结构的固-固转变，而区域Ⅴ中峰值之后趋于平滑的变化则说明硅材料向无定型状态的变化。
图 7 对代表性局部区域的键角分析和径向分布函数分析Figure 7. Bond-angle and radial distribution functions analysis for representative local regions2.3 沿[110]和[111]晶向冲击的变形与相变
利用相似的后处理技术，对单晶硅沿[110]和[111]晶向冲击压缩的变形和相变进行分析。图8～9分别为单晶硅[110]晶向中不同粒子速度下的结构演化和代表性局部的键角统计和径向分布函数分析。当粒子速度为1.8 km/s时，在靠近左端一定范围内出现了两种变化：垂直于加载方向呈带状条纹区域的晶格原子有很小的晶向偏转，紧随其后的是一定程度的原子无序状态区域。带状区域中区域Ⅰ的键角分布显示了两个主峰值，分别对应约90°和120°的键角，而区域Ⅱ的主峰值仍为109°。区域Ⅲ的径向分布函数在第1个峰值后趋于平缓且键角分布没有独立峰值，该区域为无序状态，这可能是冲击引起的熔化。对比图4的波剖面可知，区域Ⅰ和区域Ⅲ这两种变化与波阵面后方的相变波对应，而且是分离的两列波，对应着材料内部不同程度的剪应力释放和密度上升。随着冲击粒子速度的增加，两列波波速均在增加。并且当粒子速度高于2.4 km/s时，固-液相变区域逐渐占据主导。
图 8 不同冲击粒子速度下[110]单晶硅的结构演化Figure 8. Structural evolutions of single crystal Si in [110] orientation at various shock particle velocities
图 9 对图8代表性局部区域的键角分析和径向分布函数分析Figure 9. Bond-angle and radial distribution functions analysis for representative local regions图10～11为单晶硅[111]晶向中不同冲击粒子速度下的结构演化和代表性局部的键角统计与径向分布函数分析。当粒子速度为1.6 km/s时，在波阵面后一定区域内也出现了带状条纹以及靠近端部的无序状态区域。与[001]和[110]晶向不同的是，[111]晶向中带状条纹出现后，随着冲击粒子速度的增加，条纹宽度逐渐减小，并在粒子速度为2.4 km/s时逐渐消失。无序状态区域随着粒子速度的增加也不断扩展，并在3.2 km/s时扩展速度急剧增加。这一现象可能是冲击加载下温升越发显著引起的软化作用和压力引起的结构变化共同作用的结果。在相对低的冲击粒子速度下软化作用相对较弱，而压力诱发的晶格原子偏转更显著，当冲击强度增加时，温升会更显著，软化作用影响增加，晶格偏转受到一定的削弱。类似于压缩下的变形孪晶，高温会抑制变形孪晶的形成。由于硅晶体材料的变形和相变对温度、压力、应变率等因素比较敏感，因此其动态物理和力学行为会受到这些竞争性作用因素的影响。当粒子速度高达3.2 km/s时，固-液相变波波速急剧增大，这是软化区域大面积激活转化为熔化状态导致的。
图 10 不同冲击粒子速度下[111]单晶硅的结构演化Figure 10. Structural evolutions of single crystal Si in [111] orientation at various shock particle velocities
图 11 对代表性局部区域的键角分析和径向分布函数分析Figure 11. Bond-angle and radial distribution functions analysis for representative local regions2.4 讨论
从早期实验与模拟研究以及本文的模拟工作中可以看到晶体硅的复杂相变。在单轴加载条件下，沿着[001]、[110]、[111]晶向的加载过程都会发生相变，但是相变特征和机制具有各向异性。实验中已探明晶体硅在压力逐渐增加作用下发生一些列相变，如cd相→β-tin相→Imma相→sh相等。其中cd相是四配位结构，sh相为六配位结构。β-tin相是介于四配位和六配位之间的一种体心四角结构，而Imma相又是介于β-tin相和sh相之间的正交相。因此，本文模拟中涉及的固-固相变可能包括β-tin相和Imma相。Mogni等[27]基于相同的原子间势函数的模拟分析指出，[001]晶向单晶硅的冲击压缩固-固相变结构是Imma相，而且他们认为仅[001]晶向中带状区域为Imma相。通过本文中对[001]单晶硅冲击压缩下的交叉区域的偏转角和配位数以及键角的分析，可以看出交叉区域的结构应该与Imma相有所不同。该区域的配位数显示可达5甚至6，与sh相具有相同配位，但本文中尚无足够完整的结构信息用于确定具体结构。从图5中可以发现，单晶硅[001]、[110]和[111]晶向具有明显各向异性，如沿[001]晶向加载时配位数发生变化，但在[110]和[111]晶向时配位数没有变化，该现象是由于不同加载方式的原子堆积方式不同导致的。此外，不同晶向的变形与相变受到滑移系的影响，位错阻力小的滑移系更容易滑移。对于单晶硅这种金刚石结构，密排面为(111)晶面，密排方向为<110>方向，所以位错更容易沿着滑移系{111}<110>产生。所以相较于[001]晶向，[110]晶向偏转角非常小，发生结构变化的平面几乎平行于(110)平面。
晶体硅的冲击压缩研究揭示了不同的物理和力学行为，但不同实验之间的确存在较大差异。在已经报道出的晶体硅变形中除了熔化还包括位错滑移、无定型化等。McBride等[20]指出非原位探测很难识别出单晶硅冲击压缩阶段的熔化，一些中间相变也很难通过回收试样的方法准确观察到，因为晶体硅在高压卸载时可能产生其他相结构。对分子动力学模拟工作，尽管本文定性地描述了单晶硅的冲击压缩变形和相变，观测到弹性、非弹性、固-固相变以及固-液相变和固-液共存等物理特征，对比了晶体的各向异性，但是也应指出，本文所采用的势函数仍然不能准确包含单晶硅的全部变形和相变特征。因此，应对不同原子间势函数进行系统评估和对相关势函数进行修改。
总之，本文已经初步明确单晶硅冲击压缩下可依次呈现弹性、非弹性到固-固相变和固-液相变的变化趋势及其各向异性特征，对进一步理解单晶硅的冲击压缩变形和相变提供了一定的支持。针对动态加载下单晶硅相变的复杂性，后续仍需对单晶硅在动态条件下的变形和相变研究提供更深入全面的纳观尺度信息。
3. 结 论
采用大规模分子动力学方法，模拟研究了单晶硅的冲击压缩响应，获得了其冲击压力和剪切应力的变化，分析了其冲击波传播特性各向异性，并重点研究了不同冲击强度下单晶硅变形的非弹性转变和相变行为，得到主要研究结论如下：
（1）单晶硅在冲击压缩下呈现弹性到相变的变化趋势，其剪切应力随着冲击粒子速度的增加先增加，达到临界点后由于结构的变化导致应力急剧下降，其剪应力和相应的临界点存在明显的各向异性；
（2）基于偏转角和键角以及径向分布函数的分析发现，沿[001]晶向冲击压缩下可发生多种固-固相变以及固-液相变，并观察到与McBride等[20]的实验高度一致的固-液共存的重要现象。
本文研究工作可为动加载下晶体硅的相变研究提供新的纳米尺度的结果支撑，为相关实验提供新的补充。
-
图 1 分子动力学模拟中单晶硅的冲击压缩示意图
Figure 1. Schematic of shock compression of single crystal Si in molecular dynamic simulations
图 3 冲击Hugoniot状态的剪切应力(τ)与粒子速度(up) 曲线
Figure 3. Shock Hugoniot shear stress as a function of particle velocity
图 4 单晶硅中分别沿[001]、[110]和[111]晶向的冲击波剖面
Figure 4. Shock profiles in single crystal Si along [001], [110] and [111] crystal orientation
图 5 单晶硅沿[001]、[110]、[111]晶向分别在1.4、1.8、1.5 km/s冲击压缩下的可视化
Figure 5. Visualizations of single crystal Si in [001], [110] and [111] crystal orientations with shock particle velocities of 1.4, 1.8, 1.5 km/s
图 6 不同冲击粒子速度下[001]单晶硅的结构演化
Figure 6. Structural evolutions of single crystal Si in [001] orientation at various shock particle velocities
图 7 对代表性局部区域的键角分析和径向分布函数分析
Figure 7. Bond-angle and radial distribution functions analysis for representative local regions
图 8 不同冲击粒子速度下[110]单晶硅的结构演化
Figure 8. Structural evolutions of single crystal Si in [110] orientation at various shock particle velocities
图 9 对图8代表性局部区域的键角分析和径向分布函数分析
Figure 9. Bond-angle and radial distribution functions analysis for representative local regions
图 10 不同冲击粒子速度下[111]单晶硅的结构演化
Figure 10. Structural evolutions of single crystal Si in [111] orientation at various shock particle velocities
图 11 对代表性局部区域的键角分析和径向分布函数分析
Figure 11. Bond-angle and radial distribution functions analysis for representative local regions
表 1 单晶硅计算模型详细参数
Table 1. Parameters of single crystal Si sample for MD simulation
加载晶向 x轴 y轴 z轴 模型原子数 晶向 模型尺寸/nm 晶向 模型尺寸/nm 晶向 模型尺寸/nm [001] [100] 16.3 [010] 16.2 [001] 217.0 ～2.84×106 [110] [ˉ110] 16.3 [001] 16.2 [110] 217.0 ～2.84×106 [111] [ˉ1ˉ12] 16.0 [1ˉ10] 16.2 [111] 214.8 ～2.76×106 
下载: 导出CSV
-
[1] EL-KAREH B, HUTTER L N. Review of single-crystal silicon properties [M] // EL-KAREH B, HUTTER L N. Silicon Analog Components. Cham: Springer, 2020: 25–63. DOI: 10.1007/978-3-030-15085-3_2. [2] JAMIESON J C. Crystal structures at high pressures of metallic modifications of silicon and germanium [J]. Science, 1963, 139(3556): 762–764. DOI: 10.1126/science.139.3556.762. [3] MCMAHON M I, NELMES R J. New high-pressure phase of Si [J]. Physical Review B, 1993, 47(13): 8337–8340. DOI: 10.1103/PhysRevB.47.8337. [4] MCMAHON M I, NELMES R J, WRIGHT N G, et al. Pressure dependence of the Imma phase of silicon [J]. Physical Review B, 1994, 50(2): 739–743. DOI: 10.1103/PhysRevB.50.739. [5] OLIJNYK H, SIKKA S K, HOLZAPFEL W B. Structural phase transitions in Si and Ge under pressures up to 50 GPa [J]. Physics Letters A, 1984, 103(3): 137–140. DOI: 10.1016/0375-9601(84)90219-6. [6] DUCLOS S J, VOHRA Y K, RUOFF A L. Hcp to fcc transition in silicon at 78 GPa and studies to 100 GPa [J]. Physical Review Letters, 1987, 58(8): 775–777. DOI: 10.1103/PhysRevLett.58.775. [7] HANFLAND M, SCHWARZ U, SYASSEN K, et al. Crystal structure of the high-pressure phase silicon Ⅵ [J]. Physical Review Letters, 1999, 82(6): 1197–1200. DOI: 10.1103/PhysRevLett.82.1197. [8] WENTORF R H JR, KASPER J S. Two new forms of silicon [J]. Science, 1963, 139(3552): 338–339. DOI: 10.1126/science.139.3552.338-a. [9] PILTZ R O, MACLEAN J R, CLARK S J, et al. Structure and properties of silicon Ⅻ: a complex tetrahedrally bonded phase [J]. Physical Review B, 1995, 52(6): 4072–4085. DOI: 10.1103/PhysRevB.52.4072. [10] MUJICA A, RUBIO A, MUÑOZ A, et al. High-pressure phases of group-Ⅳ, Ⅲ-Ⅴ, and Ⅱ-Ⅵ compounds [J]. Reviews of Modern Physics, 2003, 75(3): 863–912. DOI: 10.1103/RevModPhys.75.863. [11] GILEV S D, TRUBACHEV A M. Metallization of monocrystalline silicon under shock compression [J]. Physica Status Solidi (B), 1999, 211(1): 379–383. DOI: 10.1002/(SICI)1521-3951(199901)211:1<379::AID-PSSB379>3.0.CO;2-4. [12] LOVERIDGE-SMITH A, ALLEN A, BELAK J, et al. Anomalous elastic response of silicon to uniaxial shock compression on nanosecond time scales [J]. Physical Review Letters, 2001, 86(11): 2349–2352. DOI: 10.1103/PhysRevLett.86.2349. [13] TURNEAURE S J, GUPTA Y M. Inelastic deformation and phase transformation of shock compressed silicon single crystals [J]. Applied Physics Letters, 2007, 91(20): 201913. DOI: 10.1063/1.2814067. [14] ZHAO S, HAHN E N, KAD B, et al. Amorphization and nanocrystallization of silicon under shock compression [J]. Acta Materialia, 2016, 103: 519–533. DOI: 10.1016/j.actamat.2015.09.022. [15] SMITH R F, BOLME C A, ERSKINE D J, et al. Heterogeneous flow and brittle failure in shock-compressed silicon [J]. Journal of Applied Physics, 2013, 114(13): 133504. DOI: 10.1063/1.4820927. [16] LIU Y X, WU X Q, WANG X, et al. Deformation behavior of single crystal silicon induced by laser shock peening [C]// Proceedings of SPIE 8796, 2nd International Symposium on Laser Interaction with Matter (LIMIS 2012). Xi’an, 2013: 87962M. DOI: 10.1117/12.2011314. [17] KISHIMURA H, MATSUMOTO H. Effect of phase transition in shock-recovered silicon [J]. Journal of Applied Physics, 2008, 103(2): 023505. DOI: 10.1063/1.2830805. [18] TURNEAURE S J, SHARMA S M, GUPTA Y M. Nanosecond melting and recrystallization in shock-compressed silicon [J]. Physical Review Letters, 2018, 121(13): 135701. DOI: 10.1103/PhysRevLett.121.135701. [19] RENGANATHAN P, TURNEAURE S J, SHARMA S M, et al. Structural transformations including melting and recrystallization during shock compression and release of germanium up to 45 GPa [J]. Physical Review B, 2019, 99(13): 134101. DOI: 10.1103/PhysRevB.99.134101. [20] MCBRIDE E E, KRYGIER A, EHNES A, et al. Phase transition lowering in dynamically compressed silicon [J]. Nature Physics, 2019, 15(1): 89–94. DOI: 10.1038/s41567-018-0290-x. [21] PAUL R, HU S X, KARASIEV V V. Anharmonic and anomalous trends in the high-pressure phase diagram of silicon [J]. Physical Review Letters, 2019, 122(12): 125701. DOI: 10.1103/PhysRevLett.122.125701. [22] DOMNICH V V. Phase transformations in silicon induced by contact loading [D]. Chicago: University of Illinois at Chicago, 2002. [23] SWIFT D C, ACKLAND G J, HAUER A, et al. First-principles equations of state for simulations of shock waves in silicon [J]. Physical Review B, 2001, 64(21): 214107. DOI: 10.1103/PhysRevB.64.214107. [24] DEMKOWICZ M J, ARGON A S. Liquidlike atomic environments act as plasticity carriers in amorphous silicon [J]. Physical Review B, 2005, 72(24): 245205. DOI: 10.1103/PhysRevB.72.245205. [25] KUMAGAI T, IZUMI S, HARA S, et al. Development of bond-order potentials that can reproduce the elastic constants and melting point of silicon for classical molecular dynamics simulation [J]. Computational Materials Science, 2007, 39(2): 457–464. DOI: 10.1016/j.commatsci.2006.07.013. [26] HIGGINBOTHAM A, STUBLEY P G, COMLEY A J, et al. Inelastic response of silicon to shock compression [J]. Scientific Reports, 2016, 6: 24211. DOI: 10.1038/srep24211. [27] MOGNI G, HIGGINBOTHAM A, GAÁL-NAGY K, et al. Molecular dynamics simulations of shock-compressed single-crystal silicon [J]. Physical Review B, 2014, 89(6): 064104. DOI: 10.1103/PhysRevB.89.064104. [28] PLIMPTON S. Fast parallel algorithms for short-range molecular dynamics [J]. Journal of Computational Physics, 1995, 117(1): 1–19. DOI: 10.1006/jcph.1995.1039. [29] ERHART P, ALBE K. Analytical potential for atomistic simulations of silicon, carbon, and silicon carbide [J]. Physical Review B, 2005, 71(3): 035211. DOI: 10.1103/PhysRevB.71.035211. [30] THOMPSON A P, PLIMPTON S J, MATTSON W. General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions [J]. The Journal of Chemical Physics, 2009, 131(15): 154107. DOI: 10.1063/1.3245303. [31] STUKOWSKI A. Visualization and analysis of atomistic simulation data with OVITO-the Open Visualization Tool [J]. Modelling and Simulation in Materials Science and Engineering, 2009, 18(1): 015012. DOI: 10.1088/0965-0393/18/1/015012. [32] GOTO T, SATO T, SYONO Y. Reduction of shear strength and phase-transition in shock-loaded silicon [J]. Japanese Journal of Applied Physics, 1982, 21(6A): L369–L371. DOI: 10.1143/jjap.21.l369. [33] GUST W H, ROYCE E B. Axial yield strengths and two successive phase transition stresses for crystalline silicon [J]. Journal of Applied Physics, 1971, 42(5): 1897–1905. DOI: 10.1063/1.1660465. [34] LI W H, HAHN E N, YAO X H, et al. Shock induced damage and fracture in SiC at elevated temperature and high strain rate [J]. Acta Materialia, 2019, 167: 51–70. DOI: 10.1016/j.actamat.2018.12.035. [35] LI W H, HAHN E N, BRANICIO P S, et al. Rate dependence and anisotropy of SiC response to ramp and wave-free quasi-isentropic compression [J]. International Journal of Plasticity, 2021, 138: 102923. DOI: 10.1016/j.ijplas.2020.102923. [36] KADAU K, GERMANN T C, LOMDAHL P S, et al. Atomistic simulations of shock-induced phase transitions [J]. AIP Conference Proceedings, 2004, 706(1): 229–234. DOI: 10.1063/1.1780223. [37] LANG JR J M, GUPTA Y M. Strength and elastic deformation of natural and synthetic diamond crystals shock compressed along [100] [J]. Journal of Applied Physics, 2010, 107(11): 113538. DOI: 10.1063/1.3448027. [38] OLEYNIK I I, ZYBIN S V, ELERT M L, et al. Nanoscale molecular dynamics simulaton of shock compression of silicon [J]. AIP Conference Proceedings, 2006, 845(1): 413–416. DOI: 10.1063/1.2263349. [39] STUBLEY P G, HIGGINBOTHAM A, WARK J S. Simulations of the inelastic response of silicon to shock compression [J]. Computational Materials Science, 2017, 128: 121–126. DOI: 10.1016/j.commatsci.2016.11.006. 期刊类型引用(1)
1. 唐晓畅，邓杰仁，莫泳晖，孟令怡，姚小虎. 动态加载下非晶合金的塑性事件演化及非近邻相互作用. 中国科学:物理学 力学 天文学. 2024(05): 6-24 . 
百度学术其他类型引用(0)
-
下载:
百度学术
计量
- 文章访问数: 588
- HTML全文浏览量: 283
- PDF下载量: 112
- 被引次数: 1