冲击加载下蓝宝石力热耦合响应的分子动力学模拟

周孟谦; 战金辉; 贺文; 操秀霞; 张伟; 刘晓星

doi:10.11858/gywlxb.20240749

冲击加载下蓝宝石力热耦合响应的分子动力学模拟

doi: 10.11858/gywlxb.20240749

周孟谦^{1, 2,},
战金辉²,
贺文²,
操秀霞³,
张伟^1,*, ,,
刘晓星^2,*, ,

1.
北京石油化工学院新材料与化工学院燃料清洁化及高效催化减排技术北京市重点实验室, 北京　102617
2.
中国科学院过程工程研究所介科学与工程全国重点实验室, 北京　100190
3.
中国工程物理研究院流体物理研究所冲击波物理与爆轰物理全国重点实验室, 四川绵阳　621999

基金项目: 国家自然科学基金（11872344）；中国工程物理研究院创新发展基金（CXXM20200700208）

详细信息

作者简介:
周孟谦（1997－），男，硕士，主要从事冲击分子动力学研究. E-mail：1010407831@qq.com

通讯作者:
张　伟（1986－），男，博士，副教授，主要从事油品脱氯技术和选择性加氢研究. E-mail：wzhang@bipt.edu.cn

刘晓星（1978－），男，博士，研究员，主要从事计算颗粒技术与应用以及多相流与多尺度数值模拟研究. E-mail：xxliu@ipe.ac.cn

中图分类号: O469; O521.2
计量
- 文章访问数: 140
- HTML全文浏览量: 36
- PDF下载量: 24
出版历程
- 收稿日期: 2024-03-07
- 修回日期: 2024-03-27
- 录用日期: 2024-04-22
- 网络出版日期: 2024-11-25
- 刊出日期: 2024-12-05

Force-Thermal Coupling Response of Sapphire under Impact Loading Based on Molecular Dynamics Simulation

ZHOU Mengqian^{1, 2
,},
ZHAN Jinhui²,
HE Wen²,
CAO Xiuxia³,
ZHANG Wei^{1,*
, ,},
LIU Xiaoxing^{2,*
, ,}

1.
Beijing Key Laboratory of Clean Fuel and High-Efficiency Catalytic Emission Reduction Technology, School of New Materials and Chemical Engineering, Beijing Institute of Petrochemical Technology, Beijing 102617, China
2.
State Key Laboratory of Mesoscience and Engineering, Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, China
3.
National Key Laboratory of Shock Wave and Detonation Physics, Institute of Fluid Physics, China Academy of Engineering Physics, Mianyang 621999, Sichuan, China

摘要

摘要: 蓝宝石因其出色的强度、硬度和光学透明度，常被选为冲击波实验中的观测窗口。深入了解蓝宝石在冲击载荷下的力学和热力学响应机制以及内部损伤原因，对准确评估其性能和稳定性至关重要。利用分子动力学模拟，从原子层面探讨蓝宝石单晶在沿(0001)晶面（C面）冲击作用下的力热响应行为。模拟结果表明，蓝宝石C面冲击作用下激活的滑移系为基于R面{ $0 \overline 1 12$ }的菱形面滑移。冲击速度为1～3 km/s时未出现滑移现象，冲击速度为4 km/s时出现菱形面滑移，冲击速度为5～6 km/s时试样出现以不规则条带为主的非均匀形变。研究表明，蓝宝石滑移系的激活不仅依赖其晶格结构，还需分剪切应力达到临界值。温度场的分析结果表明，局域温升与滑移之间存在对应关系，剪切应变集中区域的温度较高。
- 分子动力学 /
- 蓝宝石 /
- 冲击响应 /
- 损伤类型
Abstract: Sapphire is often chosen as the observation window in shock wave experiments due to its excellent strength, hardness and optical transparency. A deep understanding of the mechanical and thermodynamic response mechanisms of sapphire under impact loading and the causes of internal damage is crucial for accurately evaluating its performance and stability. In this work, molecular dynamics simulations were performed to explore the mechanical and thermal response of a sapphire single crystal under shock loading along the C-plane. The results indicate that the activated slip system after the impact loading is the rhombic plane slip based on the R-plane { $0 \overline 1 12$ }. When the impact velocity is in the range of 1−3 km/s, no slip occurs; when the impact velocity reaches 4 km/s, slip occurs. When the impact velocity reaches to the range of 5−6 km/s, the sample shows inhomogeneous deformation, mainly composed of irregular stripes. Such results suggest that the activation of the slip system in sapphire depends not only on its lattice structure, but also on the partial shear stress (which needs to reach a critical value). The analysis of the temperature field indicates that there is an intrinsic relation between the local slip and temperature increase, i.e., the formation of intense shear localization is accompanied by the higher temperature.
- molecular dynamics /
- sapphire /
- impact response /
- damage type

HTML全文

蓝宝石（Al₂O₃）因其优异的理化性能和力学强度，在航空航天、国防军事、信息技术等诸多领域得到了广泛应用。在动高压领域，蓝宝石由于具有极佳的光学透明性，常被用作试验中的光学测量窗口^[1]。在冲击作用下，蓝宝石测试窗口可能会出现发光、消光、升温、损伤、破碎等现象，从而导致冲击实验中所测数据严重失真^[2–3]。因此，深入认识冲击载荷下蓝宝石材料的力热耦合响应特性以及材料内部的力学损伤演化过程及其机制，对于理解和分析通过蓝宝石光学窗口测得的数据具有重要的指导意义。

目前，针对冲击加载下蓝宝石材料响应行为的研究主要集中在发光、消光、升温、损伤和破碎等现象的规律和机理。例如，张宁超等^[4]在轻气炮加载平台上研究了冲击加载条件下蓝宝石的发光辐射特性，结果显示，在低于60 GPa的加载压力区间，蓝宝石的辐射温度随加载压力的增加而上升，但在87 GPa压力附近出现了一个明显的降低拐点，通过对比实验结果与静高压蓝宝石高压结构相图后发现，该压力对应蓝宝石的液相、corundum相和Rh₂O₃(Ⅱ)相的三相点边界，即在此条件下蓝宝石发生了结构相变。Hare等^[2]研究了蓝宝石单晶在沿(0001)晶面(C面)、(1120)晶面(A面)和( $0\overline 1 12$ )晶面（R面）冲击下的发光辐射特性，高速摄像结果表明：蓝宝石中发光的热点并不是均匀地分布在整个冲击截面上，而是集中在某些局部区域；相同冲击压力下，沿R面的辐射发光强度要低于沿A面和C面的辐射发光强度，因此，推测蓝宝石材料中非均匀热点的形成可能与材料内部形成的剪切带有关。

受检测手段的限制，目前还难以通过实验在线精细刻画冲击加载下试样内部的损伤演化过程。随着计算能力的提升以及相关力学理论的完善，越来越多的学者采用数值模拟来探究和刻画冲击加载下材料的力热响应特性。在诸多模拟方法中，分子动力学（molecular dynamics, MD）模拟能够从原子/分子层面揭示材料的力-热-反应响应行为，得到了诸多学者的关注。Le等^[5]基于分子动力学模拟研究了蓝宝石在45 GPa压力条件下冷却过程中的相变过程，模拟结果表明：在高温条件下会出现γ-Al₂O₃和k-Al₂O₃；当温度降低至300 K时，γ-Al₂O₃会转变为δ-Al₂O₃，k-Al₂O₃转变为α-Al₂O₃。Zhang等^[6–8]采用包含高达5.4亿个原子的大规模分子动力学模拟研究了α-Al₂O₃在高速冲击载荷下的变形机制，分析了冲击加载过程中裂纹萌生过程、结构相变、撞击点附近高应变率变形模式等。

蓝宝石作为一种常用的冲击波窗口材料，在高压冲击作用下，因局部热点形成而导致的强烈冲击发光，以及因内部应力损伤导致光学散射进而造成的冲击消光都会严重干扰各类材料在冲击下的光学测量温度和压强。本研究旨在尝试从原子层面认识冲击加载下蓝宝石材料的力热响应特性，尤其是应力损伤的萌生和扩展以及加载过程中局部的温升特性，以期为后续解读蓝宝石光学窗口测得的数据提供理论支撑。

1. 分子动力学计算

1.1 模型构建

本研究工作基于Materials Studio中的Building Tools模块进行蓝宝石单胞的构建。单胞模型的晶格常数取值如下：a=b=0.4758 nm，c=1.2997 nm，α=β=90°，γ=120°。为方便后续建模及冲击模拟，将六方单胞模型通过正交变换转变为正交单胞模型，而后分别在x、y、z方向上进行扩胞。图1给出了扩胞后蓝宝石单晶的结构示意图，其中，c轴与ab面垂直，所以建模时只需简单将ab面正交即可。本研究中，构建的蓝宝石单晶试样包含约3.8×10⁶个原子，试样在x、y、z轴方向上的长度分别为16、19、104 nm， x、y、z对应的晶向指数分别为[ $10\overline 1 0$ ]、[ $\overline 1 2 \overline 1 0$ ]、[0001]。

图 1 蓝宝石单晶扩胞示意图（为显示方便，a、b、c轴上的长度均在单胞长度的基础上扩大了1倍）

Figure 1. Schematic diagrams of an enlarged model of sapphire single crystal (For visualization reason,the unit cell is replicated once in the a, b, and c axes.)

下载: 全尺寸图片幻灯片

1.2 力场参数的选择及模拟方法

分子动力学模拟中，采用能够准确描述晶格参数及结构相变的Vashishta势函数^[9]描述原子间的势函数。该势函数由双体相互作用势和三体相互作用势组成，涵盖了空间尺寸效应、库仑相互作用、电荷-偶极子相互作用和范德华相互作用。Xu等^[10]的研究结果表明，Vashishta势可以较准确地计算α-Al₂O₃的材料属性，计算结果与密度泛函理论（DFT）计算以及实验数据吻合良好。该势函数目前已被用于模拟冲击加载下多种陶瓷材料（包括AlN^[11–12]、Al₂O₃^[6–8]、SiC^[13–17]等）的力学行为。

Vashishta势函数的数学表达式为

$V={\displaystyle \sum _{i < j}{V}_{ij}^{(2)}({r}_{ij})}+{\displaystyle \sum _{i,j < k}{V}_{jik}^{(3)}}({r}_{ij},{r}_{ik})$

(1)

式中： $V$ 为系统的总势能， $V_{ij}^{(2)}({r_{ij}})$ 为双体相互作用势， $V_{ijk}^{(3)}({r_{ij}},{r_{ik}})$ 为三体相互作用势， ${r_{ij}} = | {{{\boldsymbol{r}}_i} - {{\boldsymbol{r}}_j}} |$ 为第 $i$ 个原子与第 $j$ 个原子间的距离。

双体相互作用势的表达式如下

$V_{ij}^{(2)}(r) = \frac{{{H_{ij}}}}{{{r^{{\eta _{ij}}}}}}+\frac{{{Z_i}{Z_j}}}{r}{\text{e}^{ - r/\lambda }} - \frac{{{D_{ij}}}}{{2{r^4}}}{\text{e}^{ - r/\xi }} - \frac{{{W_{ij}}}}{{{r^6}}}$

(2)

式中：H_ij为空间排斥的强度，η_ij为空间排斥力的指数，r为粒子之间的距离，Z_i、Z_j为以电子电荷e为单位的有效电荷，D_ij和W_ij分别为电荷偶极和范德华力的强度，λ为库伦项，ξ为电荷耦极项的屏蔽长度。

三体相互作用势为

$V_{ijk}^{(3)}({r_{ij}},{r_{ik}}) = {R^{(3)}}({r_{ij}},{r_{ik}}){P^{(3)}}({\theta _{jik}})$

(3)

${R^{(3)}}({r_{ij}},{r_{ik}}) = {B_{jik}}\exp \left( {\frac{\gamma }{{{r_{ij}}+{r_0}}}+\frac{\gamma }{{{r_{ik}} - {r_0}}}} \right)\Theta \left( {{r_0} - {r_{ij}}} \right)\Theta \left( {{r_0} - {r_{ik}}} \right)$

(4)

${P^{(3)}}({\theta _{jik}}) = \frac{{{{\left( {\cos {\theta _{jik}} - \cos {{\bar \theta }_{jik}}} \right)}^2}}}{{1+{C_{jik}}{{\left( {\cos {\theta _{jik}} - \cos {{\bar \theta }_{jik}}} \right)}^2}}}$

(5)

式中：B_jik为三体相互作用的强度， ${r_0}$ 为粒子间相互作用力最小时对应的距离， ${\theta _{jik}}$ 为r_ij与r_ik的夹角， ${\bar \theta _{jik}}$ 和C_jik为常量， $\Theta \left( {{r_0} - {r_{ij}}} \right)$ 为Heavisid阶跃函数。

本研究的所有数值模拟都是基于开源分子动力学模拟软件LAMMPS开展的。对于通过扩胞操作获得的蓝宝石试样，首先，基于共轭梯度法对试样进行能量最小化处理，并在模拟中沿x、y、z轴设置周期性边界条件；随后，在NPT系综下，以0.1 fs的时间步长进行总时长为40 ps的结构优化，以确保系统达到完全平衡；在冲击加载前，将z方向（冲击加载方向）的周期性边界条件更改为自由边界条件，随后，在NVE系综条件下，在试样左侧设置一个以恒定速度沿z方向运动的刚体活塞，活塞持续向右侧（沿z轴）冲击蓝宝石试样，如图2所示。本研究考察了6种不同的冲击加载工况，冲击速度在1～6 km/s范围内。

图 2 冲击加载示意图

Figure 2. Schematic diagram of the impact setting

下载: 全尺寸图片幻灯片

随着加载的进行，当冲击波抵达试样右侧面，形成一个向后传播的稀疏波。在二维数值模拟中，学者们一般采用增加试样长度的方式来避免后界面稀疏波的影响^[18]。受分子动力学三维数值模拟计算量的限制，本研究很难进一步增加试样的长度。为了弱化稀疏波的影响，借鉴Yu等^[19–20]在研究透明陶瓷冲击力学响应特性中的格点弹簧模型的模拟策略，在试样右侧也设置一个刚体活塞。当冲击波阵面传播至试样右侧面时，右侧活塞以与试样右侧自由表面相同的速度运动，进而削弱向后传播稀疏波的强度。图3对比了无、有右侧活塞情况下试样内部的应力状态。当试样右侧无刚体活塞时，稀疏波在试样内反向传播过程中形成一个明显的低压区，如图3(a)所示；当试样右侧存在一个紧贴自由面的刚体活塞时，尽管不能完全消除反向传播的稀疏波，但能明显降低稀疏波的强度，这是因为试样内部稀疏波传播导致的正应力和剪切应力的下降幅度显著降低，如图3(b)所示。

图 3 冲击速度为1 km/s时试样中心位置处沿z轴的正应力和剪切应力随时间的变化

Figure 3. Temporal variations of the normal and shear stresses at the central position of the sample along z axis at the shock velocity of 1 km/s

下载: 全尺寸图片幻灯片

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

为了分析试样的冲击响应及损伤演变过程，模拟中考察了冲击加载过程中应力、应变、速度和温度的变化。采用Virial应力定义方法^[21–22]计算应力分量。原子i的应力张量分量 ${\sigma }_{\alpha \beta }\left(\alpha ,\beta \in \left\{x,y,\textit{z}\right\}\right)$ 为

$\sigma_{\alpha\beta}=-\frac{1}{\mathit{\Omega_i}}\left[m_iv_{i,\alpha}v_{i,\beta}+\frac{1}{2}\sum\limits_{n=1}^{N_{\text{p}}}(r_{1\alpha}F_{1\beta}+r_{2\alpha}F_{2\beta})+\frac{1}{3}\sum\limits_{n=1}^{N_{\text{α}}}(r_{1\alpha}F_{1\beta}+r_{2\alpha}F_{2\beta}+r_{3\alpha}F_{3\beta})\right]$

(6)

式中：等号右侧的3项分别表示速度脉动、两体作用和三体作用对应力的贡献， $\mathit{\Omega}_i$ 为原子体积，v为速度，m为原子质量，F₁、F₂为一对原子上的力，r₁、r₂为原子在成对相互作用中的位置，α、β表示笛卡尔分量，N_p为成对近邻原子数，N_α为原子i的三体相互作用数。

本研究中，应力张量分量σ₃₃（σ_zz）为冲击应力，剪切应力 $\sigma _{\tau}$ 由下式计算得到

${\sigma _{\tau}} = \frac{1}{2}\left[ {{\sigma _{\textit{z}\textit{z}}} - \frac{1}{2}({\sigma _{xx}}+{\sigma _{yy}})} \right]$

(7)

式中： ${\sigma _{xx}}$ 、 ${\sigma _{yy}}$ 、 ${\sigma _{\textit{z}\textit{z}}}$ 为在x、y、z方向上的正应力。

剪切应变则采用由各应变张量分量 ${\eta _i}$ 计算得到的Mises应变

$\eta _i^{{\text{Mises}}} = \sqrt {\eta _{xy}^2+\eta _{y\textit{z}}^2+\eta _{x\textit{z}}^2+\frac{1}{6}\left[ {{{({\eta _{xx}} - {\eta _{yy}})}^2}+{{({\eta _{xx}} - {\eta _{\textit{z}\textit{z}}})}^2}+{{({\eta _{yy}} - {\eta _{\textit{z}\textit{z}}})}^2}} \right]}$

(8)

式中： $\eta _{xy}^{}$ 、 $\eta _{y\textit{z}}^{}$ 、 $\eta _{x\textit{z}}^{}$ 分别为xy、yz、xz平面上的应变张量， ${\eta _{xx}}$ 、 ${\eta _{yy}}$ 、 ${\eta _{\textit{z}\textit{z}}}$ 分别为x、y、z方向上的应变张量。

在分子动力学框架下，温度由原子/分子的脉动速度计算得到

$T=\frac{1}{3Nk_{\boldsymbol{\mathrm{B}}}}\sum\limits_{i=1}^Nm_i(v_{i,x}^2+v_{i,y}^2+v_{i,\textit{z}}^2)$

(9)

式中： ${k_{\text{B}}}$ 为玻尔兹曼常数， ${m_i}$ 为原子 $i$ 的原子质量， ${v_{i,x}}$ 、 ${v_{i,y}}$ 、 ${v_{i,\textit{z}}}$ 为原子/分子的脉动速度分量，N为原子/分子数量。

蓝宝石单晶中氧原子呈密排六方结构。为了刻画冲击加载过程中蓝宝石试样的结构变化，本研究采用OVITO软件自带的共近邻分析法考察了氧原子排列的演化。

2. 结果与讨论

2.1 冲击波的传播与力学响应

图4给出了冲击加载5.0 ps后试样内部的压强剖面以及沿z轴（冲击加载方向）氧原子排列方式的变化。图4(a)表明，在所考察的6个加载速度条件下，冲击速度对压强剖面的线形影响不大，波阵面都十分陡峭。图4(b)中，红色为密排六方结构的氧原子，灰色为其他结构的氧原子。由图4(b)可以看到：在1 km/s冲击速度条件下，原子排列方式未发生破坏，氧原子间仍保持密排六方结构；随着冲击速度的增大（2、3 km/s时），波阵面后方部分氧原子的密排六方结构被破坏；当冲击速度增大至4 km/s时，波阵面后方氧原子排列结构呈现出2种状态，即在冲击波前端区域，部分氧原子仍保持密排六方结构，在冲击波后端区域，氧原子密排六方结构被完全破坏；当冲击速度达到5 km/s时，波阵面后氧原子的密排六方结构被完全破坏。

图 4 (a) t =5.0 ps时不同加载速度下的压强剖面，(b) 对应的氧原子结构示意图

Figure 4. (a) Pressure profiles at t=5.0 ps under different loading velocities; (b) the corresponding schematic diagrams of the oxygen atom structure

下载: 全尺寸图片幻灯片

图5给出了试样内部中心截面上剪切应变的空间分布。可以看到，当冲击速度为1～3 km/s时，试样内部的剪切应变场均匀，没有出现明显的结构特征。这与图4(b)所示的1～3 km/s冲击速度下氧原子的排列情况相吻合，在此条件下，波阵面后方绝大部分氧原子仍维持初始的密排六方结构。当冲击速度增加至4 km/s时，波阵面后方出现明显的剪切应变集中，形成了宽度约为0.5～1 nm、与冲击方向近似呈50°夹角的条带状结构。当冲击速度进一步增加至5 km/s时，试样内部的剪切应变在空间上再次呈现出均匀分布。而当冲击速度达到6 km/s时，应变分布情况再次发生变化。此时，剪切应变场中出现了宽度约为2.2 nm的不规则条带结构，这些条带相互交错连接，并与冲击方向大致成45°角。该结果与Ma等^[23]给出的冲击加载下金属玻璃的分子动力学模拟结果定性吻合。

图 5 t=20.0 ps时不同冲击速度下试样内剪切应变等高图

Figure 5. Contour maps of shear strain inside samples under different impact velocities at t=20.0 ps

下载: 全尺寸图片幻灯片

图5表明4 km/s冲击加载下试样内部出现了明显的剪切应变集中现象。为了分析其对应的滑移系，提取了剪切应变条带中的原子信息。由于冲击过程的压缩作用，试样沿z轴的长度随着加载的进行而变化，由初始的104 nm压缩至t=20.0 ps时的75 nm左右（见图5），因此，冲击后剪切滑移面与C面的夹角与冲击前剪切滑移面上原子所在平面与C面的夹角不同。图6显示了4 km/s冲击下滑移面上所有原子的提取过程。图6(a)和图6(b)显示了20.0 ps时滑移面上的所有原子，零时刻它们的位置如(c)所示。计算表明，零时刻这些原子构成的平面与C面（xy平面）的夹角为58°6′。蓝宝石R面与C面的夹角为57°35′^[24]，推测沿冲击加载时蓝宝石内激活的滑移系是基于R面的滑移系{ $0 \overline 1 12$ }。

图 6 滑移面上原子的提取

Figure 6. Extraction of the atoms located in the shear plane

下载: 全尺寸图片幻灯片

图7给出了图6所示的4 km/s冲击条件下滑移带的形成和演化过程，图中的颜色代表剪切应变。为便于三维可视化，图7中仅给出了剪切应变大于0.18的原子。可以看到，在约11.2 ps时试样中后部出现了一个局部高剪切应变区，见图7中红圈所示；随着加载的进行，该高剪切区域沿与冲击方向成一定夹角的平面（见图7中12.0～15.0 ps时的红色方框区域）逐渐向前扩展，最终形成完整的滑移带（见图7中20.0 ps时的情形）。

图 7 冲击速度为4 km/s时试样内部滑移面的形成和演化过程

Figure 7. Initiation and development of the shear plane inside the sample under an impact velocity of 4 km/s

下载: 全尺寸图片幻灯片

2.2 滑移的演变

为理解图7所示滑移面的形成机制，分析了垂直滑移面的截面上原子的相对位置变化，如图8(a)所示。可以看到，11.2 ps时，在截面与后期滑移面相交的位置（图8(a)中红色方框区域）出现了应变集中。图8(b)给出了该区域的局部放大图，可以看到，该截面形成了一个倾斜的剪切局域化区域。进一步分析表明，该区域内发生最大相对滑移的原子为图8(b)中黑圈所示的2个氧原子。20.0 ps时（图8(c)），整个截面被分为上、下两部分，对比图8(b)可以发现，两部分发生了明显的相对滑移，两者间形成了滑移带。在该截面上，冲击加载所导致的剪切变形主要集中在滑移带，滑移带中原子的相对位置发生了显著变化，例如图8(b)红圈中相邻的2个氧原子的相对距离在20.0 ps时显著增大。

图 8 (a) 11.2 ps时垂直于滑移面截面上的剪切应变分布，(b) (a)中红框的放大（黑圈内2个氧原子的剪切应变大于0.18），(c) 20.0 ps时红框内的应变分布（黑圈内的原子对应(b)中的2个氧原子）

Figure 8. (a) Distribution of shear strain in the cross-section vertical to the later shear plane at t =11.2 ps; (b) zoomed-in view of the red box region in (a) (The shear strains of the two oxygen atoms highlighted by black circles are greater than 0.18.); (c) distributionof shear strain at t=20.0 ps (The two atoms highlighted by black circles correspond to that shown in (b).)

下载: 全尺寸图片幻灯片

为探究11.2 ps时应变局域化现象是如何形成的，进一步追踪了冲击压缩过程中原子位置的变化，结果如图9所示。4.0 ps时，冲击波阵面抵达目标区域（图8中红色方框区域），挤压该区域中的原子，如图9(b)所示； 6.0 ps时，冲击波阵面已扫过目标区域，冲击加载作用使沿冲击方向的原子间距显著减小，但原子的排列形式并没有发生变化，如图9(c)所示；11.2 ps时，绿圈内2个氧原子的排列形式发生了变化，由初始时刻的水平排列变为上下排列，如(d)所示，说明该局部区域发生了结构变化。采用OVITO软件中的Wigner-Seitz缺陷分析，得到的结果如所示，其中：Occupancy 1表示未离位的原子，Occupancy 2表示间隙原子，Occupancy 3表示有2个间隙原子。可以看出，由于冲击压缩作用，4.0 ps时出现间隙原子；随着冲击加载的进行，6.0 ps时2个间隙原子出现的概率逐渐增多；之后，由于应力释放，2个间隙原子出现的概率逐渐降低；11.2 ps时，只有一个位置出现了2个间隙原子，而该位置上的一个原子恰好对应中绿圈内的2个氧原子之一，同时，样品内部沿[ $01 \overline 11$ ]方向发生了剪切作用，导致2个氧原子的相对位置发生改变。

图 9 图8中放大区域的冲击压缩过程

Figure 9. Impact compression process of the enlarged region in Fig.8

下载: 全尺寸图片幻灯片

图 10 Wigner-Seitz缺陷分析的结果

Figure 10. Wigner-Seitz defect analysis result of defects

下载: 全尺寸图片幻灯片

以6.0 ps时原子的位移作为参考，对11.2和20.0 ps时原子的位置进行位移矢量分析，结果如所示。以6.0 ps时刻作为参考，箭头代表对应原子的位移方向，箭头长短和原子颜色代表位移大小；大球为铝原子，小球为氧原子。11.2 ps时，在所考察区域内原子沿冲击方向的位移大小已呈现出较为明显的空间差异；20.0 ps时，在所考察区域内原子位移场可分为上、下两部分，两者之间存在位移突变，即发生了滑移，滑移方向为[ $01 \overline 1 1$ ]。

图 11 不同时刻的位移矢量场

Figure 11. Displacement vector fields at different moments

下载: 全尺寸图片幻灯片

图12给出了垂直于冲击方向且与滑移面相交的截面上的剪切应变空间分布。图12(c)为零时刻截取的平面上的原子在20.0 ps时的位置与剪切应变的关系，可以看到，在应变集中的区域，原子在滑移带上呈现出明显的上下位移，并且2个特征氧原子相互远离。从图12(b)和图12(c)可以看出，沿滑移方向，应变集中区域扩展，应变增大。

图 12 垂直于冲击方向截面的剪切应变分布：(a)截面位置，(b) 11.2 ps时剪切应变的分布，(c) 20.0 ps时剪切应变的分布

Figure 12. Shear strain distributions in the cross-section plane perpendicular to the impact direction: (a) schematic diagram of the location of the target plane; (b) shear strain distribution at t=11.2 ps; (c) shear strain distribution at t=20.0 ps

下载: 全尺寸图片幻灯片

进一步分析滑移面上的缺陷。由于11.2 ps时尚未出现滑移现象，而20.0 ps时已发生滑移面相交，因此，选择16.0 ps时的原子排列作为代表进行分析。(a)给出了滑移面( $0 \overline 1 12$ )上剪切应变大于0.28的原子，该区域存在4条位错线（见(b)中红色箭头所指的蓝色线段），长度共计4.2 nm。位错线对应的伯格斯矢量方向沿A面( $\overline 2 110$ )的对角线方向[ $01 \overline 1 1$ ]和[ $0\overline 1 1 \overline 1$ ]，与滑移方向< $01 \overline 1 1$ >一致，与位错线既不平行也不垂直，表明该位错为混合位错。

图 13 (a) t=16.0 ps时的缺陷分析结果（粉色箭头为伯格斯矢量），(b) 位错分析结果（蓝线为位错线，灰色区域为缺陷）

Figure 13. (a) Defect analysis result at t=16.0 ps (Pink arrows denote the Burgers vectors.); (b) the dislocation analysis result (The blue line represents the dislocation line, and the gray area represents the defect.)

下载: 全尺寸图片幻灯片

2.3 温升现象

蓝宝石材料在高速冲击条件下会形成非均匀空间分布的局部热点^[5]。为了探究局域高温区与应力损伤（尤其是局域滑移带）之间可能的联系，本研究进一步分析了冲击加载过程中试样内部温度的空间分布。因为在本工作所考察的参数范围内，当冲击速度为4 km/s时，试样前端出现了局域剪切滑移，所以，以下将对4 km/s冲击加载工况进行分析。

图14给出了t=5.6 ps时试样内部剪切应变的空间分布以及对应的剪应力和温度剖面。可以看到，在试样内部z=30 nm附近形成了明显的剪切应变集中，与之对应的是该处剪切应力较低，而温度较高。这种特征一直延伸至剪切应变逐渐平缓的z=65 nm处。随后，剪切应力增大，温度降低。该结果与文献[25]的结果定性一致，说明剪切应变、剪应力与温度之间存在一定的关联。

图 14 4 km/s冲击速度下5.6 ps时试样内剪切应变等高图以及剪切应力和温度沿z轴的分布

Figure 14. Contour of shear strain inside the sample and the distributions of shear stress and temperature along the z-axis at 5.6 ps under an impact velocity of 4 km/s

下载: 全尺寸图片幻灯片

在分子动力学模拟的框架下，局部温度可由分子/原子的脉动速度计算得到，见式(9)。为了得到冲击加载过程中试样内部的温度分布，将试样在xyz三维空间中均分成2.6×10⁵个网格，计算每个网格内原子的平均运动速度，得到每个原子相对于该空间平均速度的脉动速度，而后，基于式(9)计算得到试样内部的温度。图15(a)显示了计算得到的试样内部温度的空间分布。为便于比较，图15(b)和图15(c)给出了对应时刻剪应力和剪切应变的空间分布。可以看到，该时刻试样前端出现了明显的剪切应变集中，形成了滑移带。对比后发现，剪切应变集中区域的剪应力低于周边区域的剪应力，而剪切应变集中区域的温度高于周边区域温度。这一结果说明，局域高温区（热点）的形成区域与剪切局域化区域存在一定的对应关系。

图 15 4 km/s冲击速度下20.0 ps时剪切应力(a)、剪切应变(b)和温度(c)的空间分布云图

Figure 15. Spatial distributions of shear stress (a), shear strain (b), and temperature (c) at 20.0 ps under an impact velocity of 4 km/s

下载: 全尺寸图片幻灯片

3. 结　论

基于分子动力学模拟研究了蓝宝石单晶在受到沿C面的冲击加载后的力-热响应行为，分析了冲击过程中剪切局域化的形成和演化过程，以及试样内部的温度变化，获得的主要结论如下。

(1) 蓝宝石在受到沿C面的冲击后，激活的滑移系为基于R面{ $0 \overline 1 12$ }的菱形面滑移系。在所考察的参数范围内，仅在冲击速度为4 km/s时出现菱形面滑移，在冲击速度为1～3 km/s和5～6 km/s范围均未出现滑移现象，说明蓝宝石滑移系的激活不仅依赖其晶格结构，还依赖分剪切应力（需达到临界值）。在高速冲击下，试样出现以不规则条带为主的非均匀形变。

(2) 滑移面( $0 \overline 1 12$ )的缺陷分析结果表明，滑移面内部出现了位错，对应的伯格斯矢量沿A面( $\overline 2 110$ )的对角线方向[ $01 \overline 1 1$ ]和[ $0 \overline 1 1 \overline 1$ ]，与滑移方向< $01 \overline 1 1$ >一致，与位错线既不平行也不垂直，表明该位错为混合位错。

(3) 剪切局域化与局部温升存在对应关系，剪切应变集中区域的温度较高。

图 1 蓝宝石单晶扩胞示意图（为显示方便，a、b、c轴上的长度均在单胞长度的基础上扩大了1倍）

Figure 1. Schematic diagrams of an enlarged model of sapphire single crystal (For visualization reason,the unit cell is replicated once in the a, b, and c axes.)

下载: 全尺寸图片幻灯片

图 2 冲击加载示意图

Figure 2. Schematic diagram of the impact setting

下载: 全尺寸图片幻灯片

图 3 冲击速度为1 km/s时试样中心位置处沿z轴的正应力和剪切应力随时间的变化

Figure 3. Temporal variations of the normal and shear stresses at the central position of the sample along z axis at the shock velocity of 1 km/s

下载: 全尺寸图片幻灯片

图 4 (a) t =5.0 ps时不同加载速度下的压强剖面，(b) 对应的氧原子结构示意图

Figure 4. (a) Pressure profiles at t=5.0 ps under different loading velocities; (b) the corresponding schematic diagrams of the oxygen atom structure

下载: 全尺寸图片幻灯片

图 5 t=20.0 ps时不同冲击速度下试样内剪切应变等高图

Figure 5. Contour maps of shear strain inside samples under different impact velocities at t=20.0 ps

下载: 全尺寸图片幻灯片

图 6 滑移面上原子的提取

Figure 6. Extraction of the atoms located in the shear plane

下载: 全尺寸图片幻灯片

图 7 冲击速度为4 km/s时试样内部滑移面的形成和演化过程

Figure 7. Initiation and development of the shear plane inside the sample under an impact velocity of 4 km/s

下载: 全尺寸图片幻灯片

图 8 (a) 11.2 ps时垂直于滑移面截面上的剪切应变分布，(b) (a)中红框的放大（黑圈内2个氧原子的剪切应变大于0.18），(c) 20.0 ps时红框内的应变分布（黑圈内的原子对应(b)中的2个氧原子）

Figure 8. (a) Distribution of shear strain in the cross-section vertical to the later shear plane at t =11.2 ps; (b) zoomed-in view of the red box region in (a) (The shear strains of the two oxygen atoms highlighted by black circles are greater than 0.18.); (c) distributionof shear strain at t=20.0 ps (The two atoms highlighted by black circles correspond to that shown in (b).)

下载: 全尺寸图片幻灯片

图 9 图8中放大区域的冲击压缩过程

Figure 9. Impact compression process of the enlarged region in Fig.8

下载: 全尺寸图片幻灯片

图 10 Wigner-Seitz缺陷分析的结果

Figure 10. Wigner-Seitz defect analysis result of defects

下载: 全尺寸图片幻灯片

图 11 不同时刻的位移矢量场

Figure 11. Displacement vector fields at different moments

下载: 全尺寸图片幻灯片

图 12 垂直于冲击方向截面的剪切应变分布：(a)截面位置，(b) 11.2 ps时剪切应变的分布，(c) 20.0 ps时剪切应变的分布

Figure 12. Shear strain distributions in the cross-section plane perpendicular to the impact direction: (a) schematic diagram of the location of the target plane; (b) shear strain distribution at t=11.2 ps; (c) shear strain distribution at t=20.0 ps

下载: 全尺寸图片幻灯片

图 13 (a) t=16.0 ps时的缺陷分析结果（粉色箭头为伯格斯矢量），(b) 位错分析结果（蓝线为位错线，灰色区域为缺陷）

Figure 13. (a) Defect analysis result at t=16.0 ps (Pink arrows denote the Burgers vectors.); (b) the dislocation analysis result (The blue line represents the dislocation line, and the gray area represents the defect.)

下载: 全尺寸图片幻灯片

图 14 4 km/s冲击速度下5.6 ps时试样内剪切应变等高图以及剪切应力和温度沿z轴的分布

Figure 14. Contour of shear strain inside the sample and the distributions of shear stress and temperature along the z-axis at 5.6 ps under an impact velocity of 4 km/s

下载: 全尺寸图片幻灯片

图 15 4 km/s冲击速度下20.0 ps时剪切应力(a)、剪切应变(b)和温度(c)的空间分布云图

Figure 15. Spatial distributions of shear stress (a), shear strain (b), and temperature (c) at 20.0 ps under an impact velocity of 4 km/s

下载: 全尺寸图片幻灯片

参考文献(25)

[1]	OCKENFELS T, VEWINGER F, WEITZ M. Sapphire optical viewport for high pressure and temperature applications [J]. Review of Scientific Instruments, 2021, 92(6): 065109. doi: 10.1063/5.0047609
[2]	HARE D E, HOLMES N C, WEBB D J. Shock-wave-induced optical emission from sapphire in the stress range 12 to 45 GPa: images and spectra [J]. Physical Review B, 2002, 66(1): 014108. doi: 10.1103/PhysRevB.66.014108
[3]	ZHANG N C, LI D, LI Y Q, et al. The radiation temperature characteristics of sapphire under shock loading [J]. Crystals, 2022, 12(10): 1364. doi: 10.3390/cryst12101364
[4]	张宁超, 王鹏, 华翔, 等. 兆巴压力下蓝宝石发光辐射特性与结构相变 [J]. 光学学报, 2019, 39(7): 0730002. doi: 10.3788/AOS201939.0730002 ZHANG N C, WANG P, HUA X, et al. Optical radiation characteristics and structural phase transition of sapphire under Megabar pressure [J]. Acta Optica Sinica, 2019, 39(7): 0730002. doi: 10.3788/AOS201939.0730002
[5]	LE C T, NGUYEN T T, NGUYEN T T, et al. Molecular dynamics simulation of phase transformation and mechanical behavior in Al₂O₃ model [J]. Vacuum, 2019, 167: 175–181. doi: 10.1016/j.vacuum.2019.06.010
[6]	ZHANG C, KALIA R K, NAKANO A, et al. Hypervelocity impact induced deformation modes in α-alumina [J]. Applied Physics Letters, 2007, 91(7): 071906. doi: 10.1063/1.2753092
[7]	ZHANG C, KALIA R K, NAKANO A, et al. Deformation mechanisms and damage in α-alumina under hypervelocity impact loading [J]. Journal of Applied Physics, 2008, 103(8): 083508. doi: 10.1063/1.2891797
[8]	ZHANG C, KALIA R K, NAKANO A, et al. Fracture initiation mechanisms in α-alumina under hypervelocity impact [J]. Applied Physics Letters, 2007, 91(12): 121911. doi: 10.1063/1.2786865
[9]	VASHISHTA P, KALIA R K, NAKANO A, et al. Interaction potentials for alumina and molecular dynamics simulations of amorphous and liquid alumina [J]. Journal of Applied Physics, 2008, 103(8): 083504. doi: 10.1063/1.2901171
[10]	XU Q Q, SALLES N, CHEVALIER J, et al. Atomistic simulation and interatomic potential comparison in α-Al₂O₃: lattice, surface and extended-defects properties [J]. Modelling and Simulation in Materials Science and Engineering, 2022, 30(3): 035008. doi: 10.1088/1361-651X/ac4d76
[11]	BRANICIO P S, KALIA R K, NAKANO A, et al. Atomistic damage mechanisms during hypervelocity projectile impact on AlN: a large-scale parallel molecular dynamics simulation study [J]. Journal of the Mechanics and Physics of Solids, 2008, 56(5): 1955–1988. doi: 10.1016/j.jmps.2007.11.004
[12]	BRANICIO P S, NAKANO A, KALIA R K, et al. Shock loading on AlN ceramics: a large scale molecular dynamics study [J]. International Journal of Plasticity, 2013, 51: 122–131. doi: 10.1016/j.ijplas.2013.06.002
[13]	BRANICIO P S, KALIA R K, NAKANO A, et al. Nanoductility induced brittle fracture in shocked high performance ceramics [J]. Applied Physics Letters, 2010, 97(11): 111903. doi: 10.1063/1.3478003
[14]	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. doi: 10.1063/1.2894188
[15]	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–2455. doi: 10.1016/j.compscitech.2008.04.040
[16]	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. doi: 10.1063/1.3152587
[17]	FENG L X, LI W H, HAHN E N, et al. Structural phase transition and amorphization in hexagonal SiC subjected to dynamic loading [J]. Mechanics of Materials, 2022, 164: 104139. doi: 10.1016/j.mechmat.2021.104139
[18]	JIANG T L, YU Y, HE H L, et al. Macroscopic shock plasticity of brittle material through designed void patterns [J]. Journal of Applied Physics, 2016, 119(9): 095905. doi: 10.1063/1.4943227
[19]	YU Y, WANG W Q, CHEN K G, et al. Controllable fracture in shocked ceramics: shielding one region from severely fractured state with the sacrifice of another region [J]. International Journal of Solids and Structures, 2018, 135: 137–147. doi: 10.1016/j.ijsolstr.2017.11.016
[20]	喻寅, 王文强, 杨佳, 等. 多孔脆性介质冲击波压缩破坏的细观机理和图像 [J]. 物理学报, 2012, 61(4): 048103. doi: 10.7498/aps.61.048103 YU Y, WANG W Q, YANG J, et al. Mesoscopic picture of fracture in porous brittle material under shock wave compression [J]. Acta Physica Sinica, 2012, 61(4): 048103. doi: 10.7498/aps.61.048103
[21]	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
[22]	HAHN E N, GERMANN T C, RAVELO R, et al. On the ultimate tensile strength of tantalum [J]. Acta Materialia, 2017, 126: 313–328. doi: 10.1016/j.actamat.2016.12.033
[23]	MA C, WANG G X, YE C, et al. Shocking of metallic glass to induce microstructure heterogeneity: a molecular dynamics study [J]. Journal of Applied Physics, 2017, 122(9): 095102. doi: 10.1063/1.5000366
[24]	郭志越. 单晶蓝宝石湿法刻蚀机理及表面形貌研究 [D]. 南京: 东南大学, 2019: 11–13. GUO Z Y. The anisotropic mechanism and topography research of wet etching of single-crystal sapphire [D]. Nanjing: Southeast University, 2019: 11–13.
[25]	胡博, 郭亚洲, 魏秋明, 等. 绝热剪切变形中温升现象的研究进展 [J]. 高压物理学报, 2021, 35(4): 040106. doi: 10.11858/gywlxb.20210728 HU B, GUO Y Z, WEI Q M, et al. Temperature rise during adiabatic shear deformation [J]. Chinese Journal of High Pressure Physics, 2021, 35(4): 040106. doi: 10.11858/gywlxb.20210728

施引文献

资源附件(0)

访问统计

点击查看大图

图(15)

计量

文章访问数: 140
HTML全文浏览量: 36
PDF下载量: 24

1. 分子动力学计算
1.1 模型构建
1.2 力场参数的选择及模拟方法
1.3 主要物理量的计算和可视化分析方法
2. 结果与讨论
2.1 冲击波的传播与力学响应
2.2 滑移的演变
2.3 温升现象
3. 结　论

冲击加载下蓝宝石力热耦合响应的分子动力学模拟

doi: 10.11858/gywlxb.20240749

作者简介: 周孟谦（1997－），男，硕士，主要从事冲击分子动力学研究. E-mail：1010407831@qq.com

计量

出版历程