全文HTML
--> --> -->国内外仅有少数研究团队开展了含能材料热力学性质的理论计算研究, 且多集中于单质含能材料. Gee等[9]和Bedrov等[10]分别开发了针对TATB分子动力学(MD)计算的力场参数; Stevens等[11]和Rykounov[12]研究了TATB的状态方程; Kroonblawd和Sewell[13,14]利用分子动力学方法系统研究了TATB的各向异性热导; Fan等[15]采用德拜声子理论方法研究了TATB热导与压力和温度的关系; 蒋文灿等[16]、Wu等[17]和Liu等[18]也基于第一性原理方法研究了TATB的声子振动模式. 目前, 单质TATB的状态方程、热容、热膨胀、热导率等性质目前已较为完善, 但是关于PBX等复合含能材料的热导率、界面热阻的研究则鲜有报道. Long和Chen[19]利用分子动力学方法最早研究了奥克托今(HMX)基PBX的界面热导率. 该方法基于线性响应理论和德拜声子理论, 认为声子在界面处的散射率等同于弹性波的散射率, 以此得到了HMX和氟聚物(F2312, F2313), TATB等界面处波的散射率. 界面热导在低频区正比于振动频率的平方, 在高频区因分子振动模的耦合而变得无规律. Long和Chen[20]采用第一性原理方法, 通过计算TATB和石墨的声子谱, 得到TATB和石墨界面的声子折射率和界面热导率, 并界定出五个主导TATB/石墨界面热导率的入射声子模. 总体来看, 目前针对含能材料的热力学性质的理论计算研究已有初步进展, 开展了部分含能材料的声子谱研究, 这对于今后从晶格动力学角度开展热导率计算有着很好的助益.
含能材料的热力学理论计算主要有第一性原理和分子动力学计算两种途径. 但是, 受制于含能材料结构的复杂性, 目前用于炸药分子的力场参数的可靠性尚未得到有效检验. 第一性原理计算是目前能准确针对复杂体系含能材料进行声子研究的有效方法, 但需要注意的是该方法仍在计算能力方面存在着一定的制约. 本文对于TATB和PVDF的声子色散问题将采用第一性原理计算方法进行. 而在热传导问题上, 分子动力学方法通过牛顿力学方程描述体系中原子的相互作用, 相较于第一性原理的计算更为简单和直观, 且计算的结果与实验符合较好. 本文用分子动力学方法进行在TATB热导率的计算. 第一性原理和分子动力学方法可以从微观尺度上解释其结构与性能之间的关系. 这不仅有助于含能材料的配方设计, 还可缩短研制周期并提高研制过程中的安全性.
本文的研究基于晶格动力学理论[21], 结合第一性原理计算和分子动力学方法系统探讨了TATB基PBX炸药的热力学性质. 首先采用冻声子法计算得到TATB和PVDF的声子色散曲线和热容等热力学性质; 根据散射错配模型, 详细讨论了TATB与高聚物界面的声子透射率和界面热导随温度的变化情况; 最后, 结合平衡分子动力学获得的热导率数据, 对TATB基PBX炸药的热导率做了简要分析. 研究结果为PBX混合炸药的热力学参数研究方法的构建, 以及对理解PBX中热传导的本质可以提供理论上的解释.
2.1.第一性原理计算模型
第一性原理计算的TATB和PVDF结构见图1. 其中TATB(C6H6N6O6)晶体属三斜晶系, 单胞中包含两个TATB分子, 共48个原子. TATB分子的构型与石墨层状结构类似. PVDF是一种高分子聚合物材料, 目前实验已发现α, β, γ, δ四种晶相的九种构型形式[22,23], 其中β相单胞中有两条全反式分子链, 分子链中相邻的CH2和CF2具有相同朝向. 目前针对β相的研究最为广泛, 为便于结果比较, 在此我们选取β相构型为代表构型, 进行本文的声子谱计算和界面热导率研究. 针对上述两种晶体, 构建了2 × 2 × 2的TATB超胞和3 × 3 × 3的PVDF超胞分别用于声子谱的计算.![](https://wulixb.iphy.ac.cn/fileWLXB/journal/article/wlxb/2019/10/PIC/10-20190075-1_mini.jpg)
Figure1. The crystal structure of TATB (a), and PVDF (b).
TATB和PVDF的声子谱由冻声子法[24]计算得到, 结构优化和力的计算利用VASP[25]软件进行. 计算中, 我们选取投影缀加平面波(PAW)赝势[26,27]结合Perdew-Burke-Ernzerhof (PBE)交换关联泛函[28,29]的方法. 其平面波截断能为550 eV, 能量的收敛标准为10–6 eV, 结构优化过程中力的收敛标准为10–4 eV/?, K点的设置为4 × 4 × 4. 二阶力常数和声子谱的处理均通过phonopy[30]软件进行.
2
2.2.简谐晶格动力学
声子是晶格振动的简正模式的量子化描述, 在晶格的动力学行为和热力学性质中起着重要的作用, 晶体中的声子色散曲线可以依据简谐晶格动力学理论[31,32]计算. 晶格动力学理论认为, 原子能够在格点附近发生振动. 对于包含N个原子的晶体体系, 其势能函数可以在平衡位置附近展开成泰勒级数的形式:![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M1.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M2.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M3.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M4.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M5.png)
2
2.3.界面热传导模型
虽然基于傅里叶定律的传热理论宏观上可以使用, 但界面的研究多涉及微观的原子或分子尺度, 此时界面声子的散射影响加大, 必须从微观角度重新考虑界面传热问题. 首先定义界面处的净热流密度和界面热阻[33], 其计算公式分别为![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M6.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M7.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M9.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M8.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M9.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M8.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M10.png)
2
2.4.平衡态分子动力学
平衡态分子动力学(EMD)计算热导率的工作, 主要是基于格林-久保公式进行[34]. 格林-久保公式建立了非平衡过程的输运系数与平衡态中相应物理量的涨落之间的联系. 输运系数等于自相关函数对于时间的积分. 因此热导率的计算公式[35,36]为![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M11.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M12.png)
本文的计算使用由Gee等[9]提出的TATB全原子力场进行分子动力学计算. 上述力场已在研究中广为应用[37,38], 其模拟结果与实验匹配较好. 力场详细参数见文献[9]. 我们构建了3 × 4 × 5 的超胞用于分子动力学模拟. 平衡态分子动力学的计算主要利用LAMMPS(large-scale atomic molecular massively parallel simulator)模拟程序包[39]. 在分子动力学的计算中, 计算时间步长设定为0.5 fs, 力的收敛标准为10–6 kcal/mol·?. 首先利用分子力学对建立的模型进行构型优化, 确保体系能量最小化. 随后在300 K温度及大气压条件下进行等压等温(NPT)分子动力学模拟, 使体系得到充分的弛豫, 模拟时长为100 ps. 最后使体系在微正则(NVE)系综下继续模拟, 总时长为1000 ps, 关联时长为100 ps.
3.1.声子谱计算
声子色散关系是固体的一个重要物理量, 可以反映出晶格振动模式与动量、能量的关系. 声子色散关系中, Γ点处的晶格振动与材料的红外和Raman特性有关. 基于冻声子法, 我们得到了TATB, PVDF的声子色散关系和声子振动模式结果(见图2(a)和图2(b)以及图3(a)和图(b)). TATB晶体的单胞有48个原子, 共有144个振动模式, 其中包含3个声学支和141个光学支.![](https://wulixb.iphy.ac.cn/fileWLXB/journal/article/wlxb/2019/10/PIC/10-20190075-2_mini.jpg)
Figure2. The phonon dispersion relation of TATB (a), and PVDF (b); (c) the heat capacity of TATB as a function of temperature; the Brillouin zone and high symmetry points of TATB (d), and PVDF (e).
![](https://wulixb.iphy.ac.cn/fileWLXB/journal/article/wlxb/2019/10/PIC/10-20190075-3_mini.jpg)
Figure3. Decomposition of the gamma-point eigenmodes of TATB (a), and PVDF (b); the weighted phonon density of states of TATB (c), and PVDF (d).
对比表1和图2(c)中理论计算结果和实验数据, 本文采用广义梯度近似(GGA)交换关联泛函的计算结果与其他研究组采用不同交换关联泛函的计算结果符合较好. 交换关联泛函及考虑范德瓦耳斯修正等因素对计算结果的影响不大. 在高频区, 尤其是大于3000 cm–1时, GGA的结果会略小于Jiang(局域梯度近似(LDA))的计算结果, 更接近于实验值. TATB晶格比热的计算数值与其他计算结果相符. 其中在温度293 K时计算的比热(498.3 J/mol·K)与实验值(495.7 J/mol·K)[41]的误差为0.52%. 结合图3和表1分析TATB的晶格振动情况, 可以发现在低频区( < 300 cm–1)尤其是频率低于100 cm–1的区域, TATB晶体中硝基(–NO2)的振动是主要的振动模式. 随着振动频率提高, 在300—1600 cm–1区间, 苯环的振动贡献有着明显的提升, 同时硝基(–NO2)和氨基(–NH2)的振动也有着一定程度的贡献. 在高频区(> 3000 cm–1), TATB晶格振动的贡献绝大部分来自于氨基(-NH2)的振动. 伴随着晶格的振动, 固体内的传热问题随之而生. 在此, 我们基于振动模式的结果, 简要分析材料热导率这一表征传热能力的重要的物理参数. 热导率是声子弛豫时间、声子群速度和声子比热的函数. 声子群速度和声子比热可由简谐晶格动力学获得, 而声子弛豫时间的获得, 则需要通过非简谐晶格动力学方法求得三阶等高阶力常数. TATB等复杂材料的高阶力常数计算仍较为困难, 目前仍暂未获得. 我们在此仅对TATB和PVDF的声子比热进行分析, 以初步探究TATB和PVDF的导热情况. 晶格动力学中晶格热容与声子态密度
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M16.png)
声子模 | 本文工作 | 蒋文灿等[16] | Liu等[18]/cm–1 | Exp.[40]/cm–1 | 不可约表示 | |||
波数/cm–1 | 偏差/% | 波数/cm–1 | 偏差/% | |||||
Q27 | 283 | –4.30 | 288 | –2.70 | 292 | 296 | E′ | |
Q30 | 284 | –4.32 | 289 | –2.69 | 295 | 297 | E′ | |
Q32 | 330 | –0.45 | 331 | –0.30 | 312 | 332 | E″ | |
Q33 | 331 | –1.01 | 332 | –0.59 | 318 | 334 | E″ | |
Q36 | 352 | –4.53 | 359 | –2.71 | 370 | 369 | E′ | |
Q38 | 354 | –4.55 | 362 | –2.42 | 371 | 371 | E′ | |
Q42 | 380 | –2.72 | 382 | –2.30 | 391 | 391 | ${\rm{A}}_2'$ | |
Q44 | 431 | –3.19 | 440 | –1.12 | 436 | 445 | E′ | |
Q46 | 432 | –3.88 | 441 | –1.78 | 438 | 449 | E′ | |
Q50 | 507 | –3.27 | 518 | –1.14 | 520 | 524 | E′ | |
Q64 | 699 | –0.72 | 680 | –3.41 | 704 | 704 | E′ | |
Q65 | 700 | –0.72 | 696 | –1.27 | 708 | 705 | E′ | |
Q88 | 869 | –0.28 | 846 | –2.87 | 870 | 871 | E′ | |
Q89 | 872 | –0.39 | 851 | –2.74 | 874 | 875 | E′ | |
Q91 | 1027 | –0.06 | 996 | –3.11 | 1026 | 1028 | E′ | |
Q94 | 1029 | –0.23 | 1001 | –2.91 | 1032 | 1031 | E′ | |
Q105 | 1217 | 0.14 | 1154 | –5.02 | 1215 | 1215 | E′ | |
Q107 | 1221 | 0.19 | 1162 | –4.67 | 1219 | 1219 | E′ | |
Q109 | 1308 | –0.34 | 1244 | –5.18 | 1320 | 1312 | E′ | |
Q111 | 1313 | –0.35 | 1250 | –5.16 | 1327 | 1318 | E′ | |
Q119 | 1442 | –0.33 | 1407 | –2.76 | 1446 | 1447 | E′ | |
Q127 | 1551 | –0.93 | 1542 | –1.53 | 1575 | 1566 | ${\rm{A}}_1'$ | |
Q129 | 1560 | –2.07 | 1548 | –2.82 | 1586 | 1593 | E′ | |
Q130 | 1564 | –2.33 | 1549 | –3.25 | 1596 | 1601 | E′ | |
Q134 | 3281 | 2.39 | 3325 | 3.78 | 3313 | 3204 | E′ | |
Q138 | 3298 | 2.66 | 3351 | 4.29 | 3334 | 3213 | ${\rm{A}}_1'$ | |
Q142 | 3399 | 2.91 | 3439 | 4.12 | 3436 | 3303 | E′ |
表1TATB部分Raman活性模计算数值与文献结果的比较
Table1.Comparison of the Raman modes of TATB crystal obtained in the present and previous calculations with experimental results.
另外需要注意到的是, 我们的计算结果中TATB在部分高对称点附近存在着一定的虚频, 该现象在蒋文灿等[16]和Wu等[17]的结果同样也有出现. 虚频的产生, 有计算中未考虑非简谐效应和温度效应的原因, 冻声子计算超胞的结构较小等. 本文计算中, TATB的2 × 2 × 2超胞有384个原子. 受制于结构的复杂性, 目前计算资源无法高效针对更大超胞做出计算. 与已有实验和理论数据对比, 本文的声子谱计算结果具有可靠性, 可用于后续的界面声子散射研究.
2
3.2.界面热导
界面热导率的计算结果如图4所示. 从图4(a)可以看出, 界面热导随温度的升高而相应提高, 并在高温时趋于平稳. 在低温区间, 界面两侧的晶格中只有少部分声子模式被激发, 随着温度的升高, 更多的声子模式随之被激发, 进而参与界面处的传热, 当温度接近或达到德拜温度时, 被激发的声子模式数量接近饱和, 声子模对界面热导的贡献不再增加. 进一步分析, 根据声子透射率公式(7)式, 可以发现界面的热传导主要发生在界面两侧振动频率重叠的地方, 且与重叠程度大小有关. 图4(b)展示了声子透射率和累积界面热导率随频率的变化关系, TATB-PVDF界面的热传导集中在声子频率小于45 THz的范围内. 其中在0—5 THz区间, 累积界面热导率已经达到总热导率的50%, 随后的声子贡献程度逐渐下降, 在大于45 THz的高频区间, 声子对界面热传导几乎无贡献. 界面热导在声子态密度重叠程度最大处达到最大值. 另外, 声子从TATB侧到PVDF侧的透射率, 明显大于反向的情况. 从图4(b)中可以看出, 声子透射率在2, 10和20 THz附近出现了拐点, 这与图3(c)TATB的加权声子态密度的变化趋势基本一致. DMM模型认为, 界面的声子透射率与入射声子的声子模、声子群速度和入射角无关, 声子通过界面的概率与该声子对应的态密度成正比, 因此我们并没有分析各声子模对声子透射率的影响. TATB声子主导着TATB-PVDF界面的声子透射, 也进一步主导着界面的热传导.![](https://wulixb.iphy.ac.cn/fileWLXB/journal/article/wlxb/2019/10/PIC/10-20190075-4_mini.jpg)
Figure4. (a) The interface thermal conductance of TATB-PVDF as a function of temperature; (b) phonon transmission of TATB-PVDF interface as a function of frequency.
以上的研究表明, PBX炸药各组分的声子谱计算, 可以作为指导炸药配方设计的一种手段. 通过计算各组分声子谱, 我们可以根据振动频率匹配程度, 初步筛选出界面热阻小的配方. 另外, 以上研究也预示着通过高分子功能化、接枝手段, 改变界面两侧振动频率的情况, 可以获得更高界面热导率的配方. 目前在实验方面, 我们尚无有效测量手段获知炸药界面的热导数据, 而且PBX炸药界面的理论研究也不多见, 我们较完整地研究了PBX炸药的界面热传导问题.
2
3.3.热导率
我们通过EMD方法获得TATB的热导率为0.58 W/mK, 结果见图5(a), 计算细节参考2.4节. EMD的交换关联时长为100 ps, 本文TATB的计算共重复10组, 并对最后的热导率结果做平均. 可以看到交换关联时间在25 ps逐渐收敛. 由于TATB的低热导率, 其声子平均自由程也仅有0.54 nm, 远小于硅、锗等中高热导率材料的自由程. 因此, 我们用于EMD模拟的参数, 包括模型尺寸和交换关联时间均较合理. 其300 K热导率计算结果也与实验值0.54 W/mK[42]符合较好. 此外我们选取PVDF的热导率实验值0.13 W/mK, 和3.2节中获得的界面热导率进行PBX热导率计算.![](https://wulixb.iphy.ac.cn/fileWLXB/journal/article/wlxb/2019/10/PIC/10-20190075-5_mini.jpg)
Figure5. (a) EMD calculation of TATB thermal conductivity; (b) thermal conductivity of TATB-based PBXs as a function of particle size.
对于复合材料的热导率, 通常我们考虑采用串并联的基本模型预测热导. 结合3.2节中DMM模型所获得的界面热导率, 我们基于串联模型计算得到了TATB-PVDF高聚物粘结炸药的热导率. 该体系的热导率由三部分贡献组成, 分别是主炸药TATB、黏结剂PVDF和TATB-PVDF界面的热导率. 当我们从某一方向对尺寸为L的炸药体系施加一个热流q, 则该体系的温度变化
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M17.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M18.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M19.png)
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M20.png)
PBX炸药的热导由(14)式得出,
![](https://wulixb.iphy.ac.cn/cn/article/doi/10.7498/10-20190075_M21.png)
1)计算了TATB和PVDF的声子色散关系, 计算结果与实验符合较好. TATB在低频区主要以硝基(-NO2)的振动为主, 中频区苯环的振动贡献有着明显的提升, 而在高频区TATB晶格振动的贡献绝大部分来自于氨基(-NH2)的振动. 低频声子在导热中有着明显的主导贡献.
2)利用DMM模型和声子色散关系, 研究了TATB与PVDF界面的热传导. 热导率随温度升高而上升, 并且在高温情况下接近不变. 通过提高界面两侧的振动频率的重合程度, 可以改善界面热导. DMM模型仅考虑了界面二声子的弹性散射情况, 多声子的非弹性散射修正的界面热传导将在接下来的研究进行.
3)基于串联模型, 分析了界面热导率和颗粒尺寸对PBX热导率的关系. 当颗粒尺寸大于100 nm时, 界面热阻对于PBX热导率的影响有限. 现实情况中, 由于PBX炸药中主炸药颗粒并非以单晶状态存在, 其粘结剂包覆的炸药颗粒是一个诸多炸药小颗粒的集合体, 颗粒间的接触热阻、以及颗粒尺寸都会影响PBX的热导.
鉴于炸药分子的复杂结构, PBX的大规模第一性原理计算目前仍难以实现. 目前计算结果的准确程度仍有提高的空间. 本文构建了PBX炸药界面热导的研究方法, 对炸药配方研究、PBX设计有着指导作用.