冲击波加载下RDX初始反应机理及分解产物光谱特征

牛振威，曾淑琼，唐 梅，葛妮娜，马兆侠

(1.西南科技大学 国防科技学院，四川 绵阳 621010；2.西南科技大学 核废物与环境安全国防重点学科实验室，四川 绵阳 621010；3.西南科技大学 环境友好能源材料国家重点实验室，四川 绵阳 621010；4.中国空气动力研究与发展中心超高速碰撞研究中心，四川 绵阳 621000)

引 言

RDX是性能优异的高能炸药，由于爆速高，常用于导弹战斗部装药、核武器的起爆装药以及火箭推进剂的添加剂。理解RDX在外界高速撞击作用下的微观起爆机理和化学反应过程，对于提高含能材料的毁伤效率、设计安全高效的新型含能材料、加深理解爆轰产物状态方程这一武器设计中的关键参数至关重要。

对于RDX分解过程的研究得到了国内外的广泛关注，但到目前为止研究主要集中在温度条件下的热分解[1-5]。Wight和Botche等[2]通过对固态RDX快速淬火和红外激光加热后发现固态RDX的热初始分解始于N-N键断裂；Irikura[5]运用密度泛函理论模拟结果同样表明固相RDX初始分解始于N-N键的断裂。相比于温度条件下的热分解，由于温度和压强存在强烈的耦合效应，冲击加载下RDX初始分解过程则更为复杂且难于模拟。2003和2005年，加州理工Goddard[6-7]课题组利用ReaxFF反应力场方法对RDX在冲击波作用下的化学反应机制进行了研究，从微观角度分析了含能材料RDX的冲击分解反应机理，指出低速冲击下RDX主要产物为NO2。2016年，陈芳和程新路[8]同样利用ReaxFF反应力场方法在恒温下对RDX了进行冲击模拟，结果表明RDX 初始分解主要为 N-NO2键的断裂。尽管这些研究给出了冲击作用下RDX可能的初始分解反应机理，但由于势函数参数选取对模拟结果影响较大，ReaxFF反应力场方法对于极端条件下材料性能的预测仍不够准确。2017年，袁姣楠等[9-11]采用半经验量子力学方法模拟了高速冲击加载下RDX的初始分解过程，结果表明初始分解是由于C-N键断裂导致的开环反应，而N-N键则受到了抑制。但是袁姣楠等[9-11]在研究中对冲击反应过程模拟时间较短(最长仅3.5ps)，且仅模拟了8、9和10km/s三个冲击波速，这不足以确切给出初始分解产物分布以及不同产物随冲击波速的变化趋势。

综上所述，虽然研究人员已经利用多种理论和实验手段，从多方面对RDX分解机理进行了详细研究并取得了巨大进步，但是其在高速撞击下初始分解产物的分布及其与温度和压强的关系仍然不明确。特别是高速冲击下含能材料分解过程不仅会形成高温高压环境，且反应时间极短，实验上准确测定反应过程中间产物的物理状态、组分、辐射系数并理解其演化过程是极其困难的。基于上述原因，本研究采用自洽电荷紧束缚密度泛函理论(SCC-DFTB)[12]结合多尺度冲击模拟技术(MSST)[13]分析RDX在较宽冲击波速度范围内初始分解机理及产物组分变化，并从理论上给出冲击Hugoniot关系，通过含时密度泛函理论(TDDFT)[14]获取主要反应产物的发射谱及辐射常数等光学特征。

1 计算方法

冲击模拟主要基于自洽电荷紧束缚密度泛函理论[12]及多尺度冲击技术[13]方法，冲击模拟过程在分子动力学程序包CP2K[15]中实现。初始结构是具有64个RDX分子的晶胞，晶格长度分别为：a=26.364nm，b=23.148nm，c=20.158nm，空间对称群为PBCA。对初始结构进行优化后，300K下晶格参数为：a=26.414Å，b=23.12Å，c=21.396Å，且α=β=γ=90.000°，密度为1.806g/cm3，这与Choi等[16]的实验结果相一致。采用一系列冲击波速度沿a轴方向进行冲击模拟，冲击波速度范围在6～13km/s。对冲击波速度低于11km/s的模拟，所选时间步长为0.1fs，而对于冲击波速度大于11km/s的模拟所选时间步长为0.05fs，总的模拟时间均在6ps以上。波函数的收敛标准为10-5a.u.。光谱计算主要基于分子光谱理论和含时密度泛函理论，其中发射光谱精细结构的计算采用FC(Franck-Condon)近似。光谱模拟在Gaussian 16程序中[17]实现。

2 结果与讨论

2.1 RDX在冲击加载下体积、温度以及应力随时间的变化

由于航天器中含能材料在超高速撞击产生的高温高压环境中会发生剧烈的化学反应，而通常超高速撞击速度约为3～7km/s。通过飞片撞击法对不同碰撞速度下材料内部所产生的冲击波速度进行预估，如图1所示。

图1 不同碰撞速度与材料内部冲击波速关系预估Fig.1 Relationship prediction between different collision velocities and the shock wave velocities inside the material

由图1可以发现，在碰撞速度为3～7km/s范围内，阻抗不同的飞片材料在靶材料(RDX)内部产生的冲击波速度有较大差异。密度高的铁作为飞片碰撞靶材料时比密度低的铝更容易在靶材料内部产生高速的冲击波，尤其是飞片速度较高时。因此，考虑到飞片材料的差异以及预估可能存在的偏差，在分子动力学模拟中将冲击波速度定义在一个较为宽泛的范围内(6～13km/s)。

RDX分解反应过程错综复杂，针对不同冲击速度下分解反应产物的动态变化进行分析，有助于进一步理解含能材料的化学反应机理。不同冲击波速度沿晶格矢量a方向加载下体系的应力、体积以及温度随时间的变化关系模拟结果如图2所示。从图2中可以发现，在冲击波扫过体系后，系统应力和温度升高，体积则被大幅压缩。从图2(a)应力变化来看，应力随着冲击波速度的增加而增加。当冲击波速度增加到9km/s时，应力约为61GPa，超过了凝聚相RDX的实验临界爆轰压力17GPa以及8.35km/s的爆轰速度，体系发生完全分解反应。这与袁姣楠等[9-11]计算结果相一致。从图2(b)体积变化来看，在模拟所采用的冲击波速度范围内(6～15km/s)，体积被压缩为原来的41%～75%。整体来说，体积压缩存在一个较长的时间过程，大致持续1.5ps，然后才趋于平衡。在达到平衡后体积随时间变化仍有一定的起伏，这与应力结果相类似，是由于在平衡时间内反应仍在不断进行。从图2(c)温度变化来看，首先，冲击波速度在6～8km/s范围内，体系温度上升缓慢，基本在2～3ps后才达到平衡状态，这是由于体系内RDX分子存在一个缓慢持续分解释放能量的过程，因此温度达到平衡所用时间较长。另外，计算所得冲击波速度为8km/s时，体系平衡温度比袁姣楠等[9-11]计算结果高300K左右，这可能是由于结构模型和时间尺度不同造成的。其中主要原因可能是，相对于袁姣楠等[9-11]模拟中仅采用24个RDX分子，在模拟中则采用了包含64个RDX分子的超胞，从而避免了由于体系过小对结果准确性的影响。另外，模拟时长均超过6ps，较长的模拟时间可以充分避免体系分解演化所引起温度变化对统计误差的影响。当冲击波速度增加到9km/s时，体系温度随时间快速增加，在1ps内即达到平衡状态，这是由于此时冲击波速度超过了凝聚相RDX的爆轰速度，系统内RDX分子结构变形以及环上化学键快速断裂，从而使得体系在较短时间内就能达到平衡状态。

图2 高速冲击加载下体系的应力、体积以及温度随模拟时间的变化Fig.2 Evolutions of stress, volume and temperature of the system with simulation time under high-speed shock loading

2.2 RDX在冲击加载下的初始反应路径

在以往的研究中，RDX的热分解主要是通过3种分解路径[1-11]：(1)N—NO2键的断裂；(2)C—H键的断裂；(3)C—N键的断裂。其中多数研究认为N—NO2键的断裂是引起RDX热分解的活性键。本研究模拟结果显示键的断裂与冲击波速度有关。通过对比小分子和长链结构产物在不同时刻的分布变化，发现整个初始分解过程由两类3种反应路径组成，如图3所示。

图3 冲击加载下固相RDX主要初始分解反应路径Fig.3 The main initial decomposition reaction paths of solid-phase RDX under shock loading

在低速冲击下(≤8km/s)，N—NO2键首先断裂，这与之前热分解结论[1-5]以及基于ReaxFF反应力场方法的低速冲击结果相符[6-7]。其中存在少量RDX分子中C—H键先断裂的情况。随着冲击波速度提高到9km/s时，N—NO2键受到高压抑制，C—H键先于N—NO2键断裂，大量H离子转移形成OH而脱离RDX分子。当冲击波速度超过10km/s后，N—NO2键受到极大抑制，C—H和C—N键则几乎同时断裂，这与袁姣楠等[9-11]所得结果基本一致。但是，本研究的计算显示，当RDX分子主环断裂解体时，相邻N—NO2和C—H键会优先形成五元环，然后断裂生成HONO分子，随后HONO分子进一步分解为更小的NO和OH分子。分析可以发现，相比于C—H和C—N键，N—NO2键在高压下更易被压缩而成为强健，从而变得更为稳定。这一结果解释了为什么低冲击速度(≤8km/s)下RDX初始分解反应和热分解反应类似，而高冲击速度(≥10 km/s)下RDX初始分解反应则完全不同，这可为研制开发更高性能的炸药提供新思路。

2.3 RDX在冲击加载下的产物组分变化

高速冲击下，伴随着高温高压环境，RDX分解过程错综复杂。不同冲击波速下，生成的中间产物和小分子不尽相同。了解小分子产物随冲击波速的变化有助于深入理解含能材料爆轰产物状态。图4显示了不同冲击作用下RDX分解产物分子数随时间的变化。

图4 不同冲击波速度下RDX初始反应小分子产物数量随时间演化曲线Fig.4 Time evolution curves of the number of small molecule products formed in RDX initial reaction under different shock wave velocities

根据图4(a)产生的小分子数目可以发现，当冲击波速度为6～8km/s时，RDX并未完全引爆，约10ps后产物随时间变化趋于稳定，产物以NO、NO2和H2O为主。此时NO分子主要来自于中间产物分子HONO的进一步分解，而NO2则主要来自于N—NO2键的断裂。因此，反应初期NO2分子含量会率先增加，这和Strachan等[6-7]的结果相一致。如图4(b)所示，当冲击波速度为9km/s时，RDX完全引爆，约5ps后产物随时间变化趋于稳定，在压缩过程中有大量NO2、NO、H2O和N2分子形成，在演化阶段，NO、H2O和N2分子数目随时间有明显的增加，最终趋于平衡，而NO2分子数目随时间的增加快速下降，最终也趋于平衡，这与袁姣楠等[9-11]给出的结果相类似。出现这种现象的原因可能是，此时体系温度在1900K以上， NO2分子不稳定并大量分解生成NO分子和游离O原子，游离O原子与游离H原子反应生成H2O。如图4(c)和(d)所示，当冲击波速度为10～13km/s时，体系温度在2200～6000K之间，约3ps后产物随时间变化趋于稳定，平衡阶段反应产物主要为N2、CO、H2O和OH，同时也伴随着大量游离态H和O原子的产生，而NO数目出现显著下降。因此，随着温度和压强的增加，氮化合物变化趋势是n(N2)>n(NO)>n(NO2)，这个趋势表明了冲击反应中温度热点对RDX分解演化的重要性，这与陈芳和程新路[8]的结论相一致。此外，特别要注意的是，游离态H和O原子的数量比接近2∶1，在体系完全反应降至常温后这些游离态的H和O原子大部分会相互反应生成H2O分子。

表1给出了不同冲击波速度下RDX分解产物的摩尔分数，从表1中可以观察到小分子产物随冲击波速度的变化趋势。其中NO2和NO含量随冲击波速度的增加呈快速下降趋势，而N2含量随冲击波速度的增加呈现出逐渐增加的趋势。相较于氮化合物，CO的含量随冲击波速度的增加而明显增加，但摩尔分数基本不超过10%。对于OH分子，除了冲击波速度为9km/s与10km/s的结果略高，其余冲击波速度下OH摩尔分数基本维持在5%左右。冲击波速度为9km/s与10km/s的结果略高是因为此时RDX完全引爆，C—N键断裂生成大量中间产物碎片。如图5所示，由于主环上C—N键断裂，生成了多种包含C—N长链和短链，这些中间产物会继续分解生成小分子，比如CN、C2N、N2O2、N2O等。

表1 不同温度和压强下RDX冲击分解产物的摩尔分数Table 1 The mole fraction of RDX shock decomposition products at different temperatures and pressures

2.4 冲击Hugoniot关系

根据表1中模拟所得密度、应力等数据，结合冲击Hugoniot关系计算冲击后粒子速度(u)，并给出了冲击后波速(Us)—粒子速度(u)的拟合曲线，结果见图6。因为计算所得数据点存在拐折，因此仅利用前5个数据点进行线性拟合。拟合结果显示初始密度为1.806g/cm3的单晶RDX声速C0约为3.98km/s。图中黑色线为拟合结果，蓝色虚线为LASL数据库中一种以RDX为主体炸药(含2.5% B方腊及2.5%树脂，平均密度为1.711g/cm3)的冲击Hugoniot实验数据及其线性延伸[18]。由于单晶RDX密度较大，计算所得声速C0要略微高于LASL数据给出的声速C0的值(3.11km/s)，同时给出的斜率则略低于LASL数据。其中差异的来源多是由于实验中并不是纯单晶RDX，包含少量其他物质且密度低。虽然结果与LASL数据库中的结果[18]有一定差异，但整体上保持一致。因此在一定程度上是可靠的。最后，由于不同冲击速度下的数据点出现了明显的拐折，因此高压和低压段的数据拟合应分开进行。

图6 RDX冲击波速度和粒子速度关系Fig.6 The relationship of shock wave velocities and particle velocities

2.5 典型产物的荧光发射谱及其精细结构

采用目前应用最为广泛的含时密度泛函理论[14]结合def2-TZVP基组模拟了主要小分子稳定产物的荧光发射谱及其精细结构。根据与实验结果对比以及相关计算经验[19]可以发现，wB97XD密度泛函计算电子发射光谱的结果较为理想，因此采用此泛函进行光谱模拟。

由于RDX是传统CHON类含能材料，它的分解会产生CN、OH、NO、NO2和CO等常见的小分子。表2列出了RDX分解产生的主要稳定产物荧光发射谱的相关辐射常数数据，对应于电子态之间的跃迁。但实际上每个电子态还对应诸多振动模式，电子可以从某一电子态的振动基态跃迁到另一电子态的各种振动态上。因此，一个电子态跃迁的峰，如果将发射光谱分辨率增加来获得精细结构，就会看到它是由许多与振动相关峰构成的，称之为振动分辨的发射光谱。将理论得到的跃迁数据用高斯函数展宽成峰形即是光谱图，半峰半宽取135.00cm-1。考虑核振动波函数获得了各种小分子产物对应的振动分辨的荧光发射谱，其计算结果如图7所示。

表2 RDX在极端条件下主要小分子产物的发射光谱信息Table 2 Emission spectrum information of main small molecule products of RDX under extreme conditions

结合表2和图7可以看到，对于OH自由基，电子从第二激发态(E=-2056.8eV)向基态(E=-2061.0eV)跃迁，发射谱线波长为297.5nm，摩尔吸收系数为36520932dm3/mol/cm。若将光谱分辨率增加获得精细结构，OH自由基振动分辨的荧光发射谱显示，两个特征峰分别位于298nm和335nm，与实验测得的带状分子光谱(300～330nm)[20]符合得较好。对于CN2分子，电子从第二激发态(E=-4007.3eV)向基态(E=-4010.8eV)跃迁，发射谱线波长为354.3nm，摩尔吸收系数为7413934801dm3/(mol·cm)；电子从第四激发态(E=-4004.8eV)向基态(E=-4010.4eV)跃迁，发射谱线波长为222.6nm，摩尔吸收系数为10629487dm3/(mol·cm)；电子分别从第一激发态、第三激发态向基态跃迁，振子强度为0，因此实验中得到的是“暗”态。CN2分子振动分辨的荧光发射谱显示，电子分别从第二、第四及第六激发态跃迁，发射光谱均由多个与振动相关峰构成。

图7 冲击加载下RDX主要小分子产物不同激发态下振动分辨的荧光发射谱Fig.7 Vibrational resolved fluorescence emission spectra of main small molecules products of RDX at different excited states under shock loading

因为振动分辨光谱的原理是不同温度下各个振动模式之间的玻尔兹曼分布发生改变，而OH自由基等双原子分子只有一个振动模式，所以不存在玻尔兹曼分布问题，因此，这些双原子分子0 K下振动分辨的发射光谱与其他温度下的光谱一致。CN2与NO3等多原子分子具有多个振动模式，以这些分子为例，其振动分辨的荧光发射谱的特征谱线线型随温度的变化关系如图8所示。研究结果表明，随着温度的升高，发射光谱的谱线线型发生显著变化，特征峰的波长基本不变，但谱线强度逐渐减弱。此外，随着温度的升高，光谱强度分布对温度的敏感性逐渐降低，不同温度下的特征谱线线型逐渐接近。这些发射谱特征谱线模拟结果有望为超高速碰撞条件下的试验诊断提供新思路。

图8 冲击加载下RDX小分子产物CN2(S2)和NO3(S3)在有限温度下振动分辨的荧光发射谱Fig.8 Vibrational resolved fluorescence emission spectra of small molecule products CN2(S2) and NO3(S3) of RDX at finite temperature under shock loading

3 结 论

(1)冲击加载下固相RDX的初始分解路径与冲击波速度有关，低速冲击下(≤8km/s)分解反应源于N—NO2键的断裂，高速冲击下(≥9km/s) ，N—NO2键受到高压抑制，初始分解反应源于C—H键和主环上C—N键的断裂。

(2)低速冲击下(≤8km/s)反应最终稳定产物以NO2、NO和H2O为主，高速冲击下(≥9km/s) 反应最终稳定产物以N2、CO和H2O为主。其中，NO含量随冲击波速度的增加呈快速下降趋势，而N2含量随冲击波速度的增加呈现出逐渐增加的趋势，冲击波速超过10km/s后，NO摩尔分数不超过5%，而N2摩尔分数在12%以上。

(3)基于模拟所得密度、应力等数据，低压段冲击波速(Us)和粒子速度(u)拟合关系为：Us=1.33u+3.98。

(4)OH自由基振动分辨的荧光发射谱显示，两个特征峰分别位于298nm和335nm，与实验测得的带状分子光谱(300～330nm)相符。多种稳定小分子产物发射谱的光谱特性将为实验光谱诊断提供参考。