基于分子动力学模拟的单晶硅冲击压缩相变研究*

2022-02-11刘梦婷李旺辉奉兰西张晓晴姚小虎

爆炸与冲击 2022年1期

刘梦婷，李旺辉，奉兰西，张晓晴，姚小虎

（华南理工大学土木与交通学院, 广东广州 510641）

晶体硅具有纯度高、晶体取向良好等特点，在电子等领域应用广泛，同时由于其复杂的结构相变机制，在相图研究中也被广泛关注。

在准静态荷载下，单晶硅表现出显著的多态性，目前已观察到13 种相结构。Mujica 等指出，单晶硅在10～11 GPa 的载荷下发生了从立方金刚石（cubic diamond, CD）结构到β-白锡（β-tin）结构的相变，并通过第一性原理计算出存在中间暂稳态 Imma 相。而增加压力时，晶体硅继续发生向简单六角形（simple hexagonal，SH）结构的转变；再进一步增加压力，变为Cmca 空间群结构，而在41 GPa 的高压下，晶体硅又从Cmca 结构转变为六方密堆积结构（hexagonal close pack，HCP）。在动态荷载实验中，Gilev 等研究了在冲击压缩下单晶硅的金属化现象，发现在冲击压力达到23 GPa 时，晶体结构的缺陷导致电导率发生了偏离。Loveridge-Smith 等用X 射线衍射技术测量了冲击加载下单晶硅的正交平面晶格参数，发现沿压缩方向单晶硅的晶格间距减少了近11%，而与压缩方向正交平面的晶格间距没有显著变化。Turneaure 等研究了单晶硅[100]、 [111]晶向在15.9、21.7 GPa 冲击压力下的动力学响应，发现其弹性波和非弹性波晶向依赖性很大，相变波晶向依赖性很小。Zhao 等对单晶硅[001]晶向进行了激光冲击压缩和回收实验，观察到两个塑性响应：靠近冲击表面的大块非晶层和沿{111}平面产生的滑移带，该结果揭示了单晶硅在一维应变条件下的塑性响应。Smith 等对单晶硅[001]、[110]和[111]晶向进行了冲击实验，揭示了单晶硅材料冲击响应的各向异性，发现了其在冲击压力高于13 GPa 时转变为简单六角形结构。Liu 等研究了单晶硅沿[100]晶向冲击压缩的塑性变形，发现其塑性变形能力随着温度增加而增加。Kishimura 等对在38 GPa 冲击压力下的单晶硅进行回收并进行了X 射线分析和拉曼光谱分析。XRD 分析显示，单晶硅在11 GPa 冲击压力下有少量亚稳态相产生，在38 GPa 冲击压力下有CuSi 这一物质生成，而在其余加载压力下并没有明显的高压相和亚稳相存在。拉曼光谱的分析显示，其峰值的整体偏移是由晶体尺寸减少引起的，该结果说明了单晶硅冲击压缩下的晶粒化过程。Turneaure 等用实时XRD 检测了单晶硅在冲击压力低于54 GPa 下的晶体结构，发现单晶硅在冲击压力为31～33 GPa 时产生冲击熔化，并从熔融边界再结晶到密堆积六角形结构。由于硅和锗属于同族元素，具有相似的物理力学特性，因此Renganathan 等用原位XRD 检测了锗在[100]晶向的相变，发现锗在冲击压力为15.7 GPa 时转变为白锡结构，在冲击压力为31.5 GPa 时转变为熔融态；应力卸载后，又转变为立方金刚石结构，该结果表明锗从白锡结构到立方金刚石结构的相变是可逆的。McBride 等用实时X 射线检测技术揭示了多晶硅的冲击相变行为，观察到不同压力下的多种固-固相变以及固-液相变。Paul 等则利用密度泛函理论、进化算法和晶格动力学等多种方法，构建了4 TPa 和6 000 K 下硅的相图，为极高压相图的研究提供了参考。

尽管以上实验与计算研究结果丰富，但实验之间也存在较显著的差异。一方面，实验直接获取的物理量过少，难以直接表征结果；另一方面，各个实验设备和实验过程也存在差异，变量如应变率、加载时间等不一致，而这些因素都会对单晶硅的相变有显著影响。同时，单晶硅作为一种脆性材料，不仅加工困难，在冲击压缩时样品也很难回收，而且单晶硅的相图比较复杂，存在稳态相和非稳态中间相等多个相，因此它在冲击压缩下的物理和力学特性尚无完全统一的结论。此外，单晶硅冲击实验中加载时间较短，很难实时测量其相变过程中的晶格结构。所以，尽管检测微观结构变化的原位实时诊断技术取得了很大的进展，但对于可逆相变往往还是求助于把宏观力学量的实验测试结果与静高压数据进行比较，以确认新相的结构、相变的类型及性质，这也会导致研究结果的差异。

由于分子动力学具有超高的时间和空间分辨率，在微纳米尺度研究材料的物理和力学特性时具有独特的优势。而在当前的超算能力下，分子动力学模拟规模可以达到千万甚至数十亿原子，并能同时获取各个原子运动的轨迹，进而从原子尺度到微纳米尺度揭示材料的变形与力学参量的关系，使得其能够对理论计算和实验研究提供重要的补充。在分子动力学模拟单晶硅方面也有诸多发现。Swift 等用第一性原理计算得到单晶硅的物态方程，这个物态方程可以反应单晶硅从金刚石结构到体心四方结构的转变。Demkowicz 等用SW 势函数模拟非晶硅在施加恒定压力下的变形情况，发现压缩过程中存在类晶体和类液体部分，后者更容易发生塑性流动。Kumagai 等改进了Tersoff 势函数，使之既能模拟弹性常数又能准确估测熔点。Higginbotham等通过分子动力学模拟，表明压力诱导相变是弹性波异常的原因。Mogni 等则基于Tersoff-like 势函数模拟[001]晶向单晶硅在激光冲击下的相变，发现了新的Imma 相。Zhao 等研究了基于分子动力学模拟的冲击压缩下单晶硅的力学响应，发现硅在冲击压力低于10 GPa 时为弹性响应，在稍高压力下，由于位错和缺陷沿着特定的晶体学方向扩展，会进一步产生新的缺陷，此外在更高压力下会产生交叉的非晶带。

尽管基于分子动力学的单晶硅变形与相变的模拟研究成果丰富，但仍有若干关键问题尚未明确，包括单晶硅的动态压缩各向异性，以及单晶硅复杂相变机制的物理描述和不同因素对其相变的影响规律和机制等。深入研究这些关键问题有利于更加全面地认清单晶硅在动态压缩下的变形和相变的物理和力学特性。本文中将基于大规模分子动力学方法模拟研究单晶硅的冲击压缩响应，从微纳米尺度认识晶体硅的冲击压缩特性，揭示其在动加载下的各向异性和结构相变等特征。

1 研究方法

1.1 模型的构建与冲击压缩模拟方法

基于广泛使用的LAMMPS平台，采用分子动力学方法模拟研究单晶硅的冲击压缩响应，主要对其非弹性变形和相变行为展开分析。晶体硅通常情况下具有立方金刚石结构，在300 K 温度下晶格常数为0.543 nm。模型的建立是由单个晶胞沿、、晶向分别复制30、30 和390 个晶胞得到的超胞结构。因此，模型尺寸约为16.3 nm×16.3 nm×212 nm，见表1。

表1 单晶硅计算模型详细参数Table 1 Parameters of single crystal Si sample for MD simulation

晶体硅中Si-Si 原子间的相互作用采用改进的半经验Tersoff 势函数描述，并使用Erhart 等优化的参数。与改进前的势函数相比，该势函数可以成功模拟晶体硅的多种性质，包括弹性常数和一般热力学性质，以及在静水压缩和冲击压缩条件下立方金刚石相到β-白锡相的结构特征。

在施加冲击压缩前，先对构建的初始模型进行弛豫，弛豫过程中保持三维方向均为周期性边界条件，使整个晶体模型在NPT 系综下弛豫了20 ps。从能量和残余应力的曲线上确认了材料模型为能量趋于稳定且残余应力在误差允许范围内趋于零的状态，此时获得了合理稳定的初始结构。在冲击压缩加载阶段，为了消除边界对应力波的影响，在沿轴施加自由边界条件的同时，沿、方向施加周期性边界条件。模拟平板冲击采用的是活塞法，活塞对自由表面产生冲击波。单晶硅冲击加载的示意图如图1所示。冲击压缩模拟过程中系统处于NVE 系综，并且不施加温度控制，时间步长为1 fs。这样的加载方式和效果，与典型的平面冲击实验相似，均为一维应变问题。为揭示单晶硅冲击响应的各向异性，尤其是冲击结构的演化，分别对[001]、[110]和[111]晶向进行系列冲击压缩研究，冲击粒子速度从0.3 km/s逐渐增加至3.2 km/s。

图1 分子动力学模拟中单晶硅的冲击压缩示意图Fig. 1 Schematic of shock compression of single crystal Si in molecular dynamic simulations

1.2 主要物理量的计算和可视化分析方法

为了深入分析单晶硅的冲击压缩响应，利用OVITO进行可视化分析，并采用配位数、中心对称参数、键角分析等多种分析技术，对晶体硅的变形和结构变化进行识别和分析。

2 结果与讨论

2.1 冲击响应与各向异性

分子动力学模拟计算结果与不同文献中实验所得的晶体硅冲击Hugoniot 应力和体积变化的关系如图2 所示。从图2 中可以看出，在体积压缩较小时，所有结果高度吻合；压缩较大时，实验的应力-体积关系呈现两大分歧：一是以Goto 等和Gust 等的结果为代表，这两个独立研究团队的结果比较接近，二是以Turmeaure 等和McBride 等的实验结果为代表。因此，对晶体硅不同冲击实验研究的差异值得关注，这可能与实验手段、技术和材料样品差异有关。一方面，Gust 等早期采用爆炸驱动的冲击波压缩实验，而Goto 等和Turmeaure 等采用轻气炮驱动的平板冲击装置对不同晶向的单晶硅进行冲击实验，McBride 等则采用激光驱动的应力脉冲加载对多晶硅进行冲击压缩实验；另一方面，尽管所有实验都是运用Rankine-Hugoniot 关系推算试样的内部压力，但压力推算中所需的粒子速度或相变波速度测定的差异显著，例如Gust 等提及其结果与相关历史数据的对比时，对应于HEL 和相变点的临界体积（或密度）是很吻合的，但压力数值却存在差异，他们将此归因于粒子速度的测定差异。值得注意的是，两大存在分歧的压力-体积关系均有相互独立的研究团队的结果支撑，可能预示着深层次的原因在于速度测定的技术差异与硅材料响应有关。Goto 等的轻气炮驱动的平面冲击实验中采用了与Gust 等相似的测量技术（斜镜与条纹相机）对硅试件的自由面测速，其Hugoniot 压力推算值也与Gust 等非常吻合；然而，Turneaure 等的轻气炮驱动平面冲击实验与McBride 等的激光冲击实验则采用了VISAR 测试系统对硅/LiF 窗口或者硅自由面速度测定，其结果相接近。此外，Smith 等指出硅自由表面端存在的应力松弛和初始的速度溅射下的自由表面测量结果可能会影响应力和体积的关系测定。因此，这些结果在压力数值上的差异深刻反映了实验技术差异和数据处理方法对实验结果的重要影响。

图2 不同实验[13,20,32-33]与模拟中晶体硅的冲击压缩Hugoniot 应力和体积变化的关系Fig. 2 Shock Hugoniot stress as a function of volume change in various experiments[13,20,32-33] and simulations of Si crystals

与实验结果相比，本文计算结果与Goto 等和Gust 等的结果更接近，但也在一定程度上高于实验值。尽管如此，基于该原子间势函数模拟计算的应力-体积整体变化趋势与实验结果相似。例如，随着体积压缩的增加，可以明显看出在体积比低于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 达到最大负剪切应力，而后逐渐减小。负剪切应力现象存在于单晶碳化硅和单晶金属铁中，学者们普遍认为这是由于材料过度弛豫释放的结果，并且具有晶向依赖性，多发生在[001]晶向。从剪切应力的计算方法角度分析，负剪应力意味着在冲击加载中材料内部的横向应力高于纵向应力。最终，[001]、[110]和[111]晶向的剪切应力分别在2.0、1.9、2.6 km/s 之后趋近零，因此材料失去剪切强度。此外，结合图1 中应力和体积压缩的关系，单晶硅在[110]晶向的奇异表现还体现在弹塑性或相变转变的临界状态，相比于[001]和[111]晶向的在约0.85 体积压缩比时的“弹塑性”转变，[110]晶向则表现出奇异的弹性特征。这种冲击压缩下奇异的弹性响应现象在金刚石等材料的分子动力学模拟研究也有过发现，这些材料呈现出非单调的剪应力-应变关系。由于剪切应力是塑性等变形的驱动力，因此，单晶硅冲击压缩下奇异的弹性响应应该与其剪切应力的非单调性相关。

图3 冲击Hugoniot 状态的剪切应力(τ)与粒子速度(up) 曲线Fig. 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 等和Stubley等的研究工作也得到了体现。

图4 单晶硅中分别沿[001]、[110]和[111]晶向的冲击波剖面Fig. 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 冲击压缩下的可视化Fig. 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

2.2 沿[001]晶向冲击的变形与相变

为进一步探究单晶硅的冲击结构相变等特征，重点分析[001]晶向单晶硅的情况。图6 中给出了不同粒子速度下的单晶硅冲击响应可视化分析，通过晶格原子的偏转角可以显著地区分出冲击压缩下单晶硅变形和相变的局部变化，蓝色区域原子的偏转角非常小，黄色区域原子较初始结构偏转约54°，绿色区域原子较初始结构偏转45°。图中[001]晶向的相变带状区域延伸方向与加载方向呈57°夹角，与偏转角接近，这样的少许差异与材料处于压缩状态相关。这表明相变波的相变带遵循特定的晶体学取向，是沿着靠近{111}晶面向前传播，其中{111}晶面与{001}晶面夹角为54°。随着粒子速度的增加，相变带状区域宽度有略微增加的趋势。当粒子速度为2.0 km/s 时，靠近冲击加载端的区域出现无序的状态，这个无序区域随着粒子速度的增加而扩展。当粒子速度为2.4 km/s 时，一个显著的无序状态与晶体状态同时存在于相同加载方向位置。考虑到强冲击可能引起剧烈温升和更高的压力状态，因此可能引起单晶硅的熔化，这意味着在此冲击条件下单晶硅处于固-液共存的区域。值得注意的是，本文中所发现的冲击波阵面后方固-固相变和固-液相变共存的现象也出现在其他两个晶向，该发现与McBride 等的激光冲击实验诊断结果一致，本文中冲击压缩下固-液共存甚至完全熔化所对应的临界体积压缩比，也与McBride 等的实验高度吻合。如图6 所示，当粒子速度高达2.8 km/s 时熔化区域已经扩展很大，并在3.2 km/s 时波阵面后几乎完全熔化。

图6 不同冲击粒子速度下[001]单晶硅的结构演化Fig. 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 对代表性局部区域的键角分析和径向分布函数分析Fig. 7 Bond-angle and radial distribution functions analysis for representative local regions

2.3 沿[110]和[111]晶向冲击的变形与相变

利用相似的后处理技术，对单晶硅沿[110]和[111]晶向冲击压缩的变形和相变进行分析。图8～9分别为单晶硅[110]晶向中不同粒子速度下的结构演化和代表性局部的键角统计和径向分布函数分析。当粒子速度为1.8 km/s 时，在靠近左端一定范围内出现了两种变化：垂直于加载方向呈带状条纹区域的晶格原子有很小的晶向偏转，紧随其后的是一定程度的原子无序状态区域。带状区域中区域Ⅰ的键角分布显示了两个主峰值，分别对应约90°和120°的键角，而区域Ⅱ的主峰值仍为109°。区域Ⅲ的径向分布函数在第1 个峰值后趋于平缓且键角分布没有独立峰值，该区域为无序状态，这可能是冲击引起的熔化。对比图4 的波剖面可知，区域Ⅰ和区域Ⅲ这两种变化与波阵面后方的相变波对应，而且是分离的两列波，对应着材料内部不同程度的剪应力释放和密度上升。随着冲击粒子速度的增加，两列波波速均在增加。并且当粒子速度高于2.4 km/s 时，固-液相变区域逐渐占据主导。

图8 不同冲击粒子速度下[110]单晶硅的结构演化Fig. 8 Structural evolutions of single crystal Si in [110]orientation at various shock particle velocities

图9 对图8 代表性局部区域的键角分析和径向分布函数分析Fig. 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]单晶硅的结构演化Fig. 10 Structural evolutions of single crystal Si in [111]orientation at various shock particle velocities

图11 对代表性局部区域的键角分析和径向分布函数分析Fig. 11 Bond-angle and radial distribution functions analysis for representative local regions

2.4 讨论

从早期实验与模拟研究以及本文的模拟工作中可以看到晶体硅的复杂相变。在单轴加载条件下，沿着[001]、[110]、[111]晶向的加载过程都会发生相变，但是相变特征和机制具有各向异性。实验中已探明晶体硅在压力逐渐增加作用下发生一些列相变，如cd 相→β-tin 相→Imma 相→sh 相等。其中cd 相是四配位结构，sh 相为六配位结构。β-tin 相是介于四配位和六配位之间的一种体心四角结构，而Imma 相又是介于β-tin 相和sh 相之间的正交相。因此，本文模拟中涉及的固-固相变可能包括β-tin 相和Imma 相。Mogni 等基于相同的原子间势函数的模拟分析指出，[001]晶向单晶硅的冲击压缩固-固相变结构是Imma 相，而且他们认为仅[001]晶向中带状区域为Imma 相。通过本文中对[001]单晶硅冲击压缩下的交叉区域的偏转角和配位数以及键角的分析，可以看出交叉区域的结构应该与Imma 相有所不同。该区域的配位数显示可达5 甚至6，与sh 相具有相同配位，但本文中尚无足够完整的结构信息用于确定具体结构。从图5 中可以发现，单晶硅[001]、[110]和[111]晶向具有明显各向异性，如沿[001]晶向加载时配位数发生变化，但在[110]和[111]晶向时配位数没有变化，该现象是由于不同加载方式的原子堆积方式不同导致的。此外，不同晶向的变形与相变受到滑移系的影响，位错阻力小的滑移系更容易滑移。对于单晶硅这种金刚石结构，密排面为(111)晶面，密排方向为<110>方向，所以位错更容易沿着滑移系{111}<110>产生。所以相较于[001]晶向，[110]晶向偏转角非常小，发生结构变化的平面几乎平行于(110)平面。

晶体硅的冲击压缩研究揭示了不同的物理和力学行为，但不同实验之间的确存在较大差异。在已经报道出的晶体硅变形中除了熔化还包括位错滑移、无定型化等。McBride 等指出非原位探测很难识别出单晶硅冲击压缩阶段的熔化，一些中间相变也很难通过回收试样的方法准确观察到，因为晶体硅在高压卸载时可能产生其他相结构。对分子动力学模拟工作，尽管本文定性地描述了单晶硅的冲击压缩变形和相变，观测到弹性、非弹性、固-固相变以及固-液相变和固-液共存等物理特征，对比了晶体的各向异性，但是也应指出，本文所采用的势函数仍然不能准确包含单晶硅的全部变形和相变特征。因此，应对不同原子间势函数进行系统评估和对相关势函数进行修改。

总之，本文已经初步明确单晶硅冲击压缩下可依次呈现弹性、非弹性到固-固相变和固-液相变的变化趋势及其各向异性特征，对进一步理解单晶硅的冲击压缩变形和相变提供了一定的支持。针对动态加载下单晶硅相变的复杂性，后续仍需对单晶硅在动态条件下的变形和相变研究提供更深入全面的纳观尺度信息。