引 言
高聚物黏结炸药(PBX)是由炸药晶粒及高分子黏结剂等经过一定工艺压制加工而成的高能炸药[1]．它在运输、存贮、使役中常会遇到以热力变化为主的复杂环境，使炸药出现损伤等物理力学性能发生改变的现象，进而影响其爆轰性能与稳定性．由于使用实验复现实际工况操作复杂且耗时、PBX或其代用材料试样制作繁琐，因此采用较为简便的数值模拟方法对PBX在不同热力环境中损伤行为进行研究也是十分重要的手段．
在现有对PBX的仿真实验中，主要存在两种几何建模方法．一种方法是将炸药视为均质材料: 如Wang等[2]将PBX视作均质材料，数值模拟了HMX基聚合物黏结炸药的燃点反应，研究了热点和潜在热点的位置．Hu等[3]在模拟了经历温度循环的均质炸药后认为，炸药晶粒的相变在一定程度上会使炸药内部产生裂纹．这种视炸药整体为均质的建模方法能突出材料整体在进行温度循环后发生的变化，但无法体现出炸药晶粒与黏结剂由于各自热力学性能不同而产生的界面损伤．另一种是分别为炸药晶粒和黏结剂建模，通过指定炸药晶粒的形状或者运用图像识别等技术建立一种可以分别表示炸药晶粒和黏结剂几何结构的几何模型：Wu等[1]根据Voronoi理论[4-5]建立的PBX细观模型保证了晶粒分布的随机性且有较高的体积分数，但是每个炸药晶粒的大小随机，与实际相的晶粒级配不一致．Barua等[6]和Ambos等[7-8]基于炸药内部的电镜图像，通过图像识别、Fourier变换等手段，生成了与PBX实际细观结构较为一致的二维几何模型，但这类以电镜图像为基础的炸药模型只是对于PBX特定区域的复现．
已有的研究表明，以PBX 9502[9]为代表的PBX由温度变化引起的损伤主要表现为炸药晶粒与高分子黏结剂之间界面的脱黏，以及不可逆增长、热分解等．Pi等[10]通过实验和数值模拟探究了在不同温度条件下炸药对冲击的敏感程度，认为炸药的敏感性与温度在一定程度上有线性关系．林聪妹等[11]在实验的基础上综合分析了TATB基PBX的力学性能和热物理性能，认为加入高玻璃化转变温度和高力学强度的黏结剂有助于降低PBX在热冲击下的损伤程度．韦兴文等[12]、Willey等[13]、张伟斌等[14]国内外学者分别采用实验与数值模拟的方法研究了热冲击对炸药产生损伤的影响．从这些研究可以发现，在热载荷作用下，尤其是在降温过程中，PBX更易产生损伤，但由于相关研究均采用均质材料代表PBX，因此在温度变化时晶粒-黏结剂界面的损伤机理仍有待研究．
为了探究晶粒-黏结剂界面的损伤机理，本文建立了带有微结构的PBX二维几何模型．相较于三维模型，二维模型能更直观地反映在温度变化时的晶粒-黏结剂界面损伤的变化，且建模计算更为简便．本文根据Voronoi模型的原理，结合Monte Carlo的级配思想，提出了一种炸药晶粒带有级配的、高分子黏结剂层厚度非均匀的二维几何模型生成方法．该方法能根据实际需求对这些晶粒的棱角进行平滑和优化，更好地控制炸药晶粒的形状以及晶粒的级配，可以得到体积分数超过90%的PBX 9502几何模型．采用弹塑性和黏弹塑性分别描述炸药晶粒和黏结剂的力学性质，引入内聚力模型描述晶粒与黏结剂界面的力学性能．在此基础上，对含有微结构以及不同体积分数的炸药进行了热力学行为的数值分析，研究了升温和降温过程对炸药以及晶粒与黏结剂界面力学行为的影响．
1 高聚物黏结炸药(PBX)的几何建模
PBX 9502的细观结构复杂，从组成成分来看主要有炸药晶粒和高分子黏结剂．为了反映炸药的真实结构，本文提出了一种几何模型生成方法．
1.1 PBX二维几何模型的生成过程
PBX的二维几何建模的整体流程如图1所示．其中基于Voronoi理论的几何建模以MATLAB编程为主，对模型的进一步优化则在AutoCAD中实现．在建模时将用MATLAB编程所得的模型转化为命令流文件导入AutoCAD中进行进一步优化．该方法既能保证炸药晶粒形状的随机性，又能通过人工优化保证后续网格划分的质量．
1.2 几何建模的编程过程
Voronoi图形又叫Thiessen多边形[4-5]，它在几何建模中运用广泛，在MATLAB中主要是通过定义基点坐标的分布来控制整体Voronoi图形的生成．本文根据Monte Carlo的级配思想[15]，采用了三种半径不同的圆，并通过编程的方法控制这些圆的位置，以这些圆心作为Voronoi图形的基点最终生成相应的Voronoi图形．具体圆的尺寸与相应的个数如表1所示．
表1 各级级配圆的半径与数量
Table 1 The number of circles with different sizes
circle radius r/μmnumber of circles N5212263851153
图1 PBX二维几何建模流程
Fig. 1 The procedure of 2D geometrical modeling of PBX
文本采用的圆共有450个，以此为基础生成相应级配的炸药晶粒．现有含能晶粒的制备技术已十分完备，可控制其粒径在300 μm以内，因此文中采用的级配关系具有一定代表性．将这些圆随机投放到1 mm×1 mm的区域中，通过限制各圆之间的距离，使每个圆不能相交或相切，同时保证这些圆在投放区域中保持完整．最终所得的圆的分布图像如图2(a)所示．
(a) 级配圆的分布 (b) MATLAB中生成的Voronoi图像
(a) Distribution of graded circles (b) The Voronoi diagram generated in MATLAB
(c) 上色的Voronoi图像 (d) 级配圆与Voronoi图像的重合
(c) The Voronoi diagram filled with colors (d) The coincidence between graded circles and the Voronoi diagram
图2 编程生成的图像
Fig. 2 Diagram generated in programming
注 为了解释图中的颜色，读者可以参考本文的电子网页版本，后同．
图2(b)是以这些圆的圆心为基点的Voronoi图形，图2(c)为上色后的Voronoi图形．将圆的分布图与Voronoi图形重叠(如图2(d))，可发现Voronoi图形各个基点均为相应圆的圆心，可证明该建模方法能与Voronoi建模理论相对应．
将每个多边形进行一定范围的随机缩放以生成厚度不均匀的黏结剂层，并将所得的二维模型转化成命令流文件导入AutoCAD中，如图3(a)所示．由于生成的模型颗粒数量众多，为控制后续的数值模拟计算规模，在生成的模型中随机框选出0.5 mm×0.5 mm的区域，如图3(b)．将此框选出的模型作为PBX的基础模型．采用Voronoi理论进行几何建模会带有一定的缺陷，即多边形的每条边均为对应相邻两个基点连线的垂直平分线，因此在实际建模时会将图3(c)中代表大晶粒的红色多边形与其四周中小尺寸的多边形的大小平均化．在后续对模型的优化时，通过人工干预等手段突出大晶粒的尺寸并将各多边形进行倒角处理，以保证后续进行数值分析时网格划分的可靠性．经过优化后的最终二维几何模型如图3(d)所示．
(a) 导入命令流文件生成的几何模型 (b) PBX的基础模型
(a) The geometric model generated with (b) The basic model for PBX
the imported command stream file
(c) 大晶体在PBX的基础模型的位置 (d) 优化后的PBX模型
(c) The positions of large crystals in (d) The optimized PBX model
the basic model for PBX
图3 AutoCAD中的模型优化
Fig. 3 Modifying in AutoCAD
图4(a)～4(e)为根据黏结剂层不同厚度范围并采用上述方法生成的二维几何模型．通过图像识别等手段，得出PBX二维模型的炸药晶粒体积分数分别为72.4%，76.0%，81.7%，85.3%和91.6%，它们对应的黏结剂层厚度范围分别为4～10 μm， 3～9 μm， 2～7 μm， 1.5～5.5 μm和1～3 μm．
图4(f)为每个模型的黏结剂层平均厚度与炸药晶粒体积分数的关系．从图中可以发现，随着黏结剂平均厚度的减少，炸药晶粒的体积分数不断增加．图4(g)为体积分数91.6%的模型所含炸药晶粒的级配图，其晶粒分布近似单峰正态分布，平均粒径约为43.5 μm．图4(h)为体积分数91.6%的模型的平均粒径与实际晶粒大小的对比，其中横坐标是该模型中炸药晶粒由小到大排列的序号，纵坐标是各个晶粒的粒径， 中间横线表示的是平均粒径43.5 μm．通过图4(g)和图4(h)可以看出, 该模型的大部分晶粒粒径在40～60 μm之间．本文所采用的级配关系与实际晶粒的粒径分布较为一致， 可认为能反映PBX内部炸药晶粒的真实情况．
(a) 晶粒体积分数(72.4%) (b) 晶粒体积分数(76.0%)
(a) The crystal volume fraction(72.4%)(b) The crystal volume fraction(76.0%)
(c) 晶粒体积分数(81.7%) (d) 晶粒体积分数(85.3%)
(c) The crystal volume fraction(81.7%)(d) The crystal volume fraction(85.3%)
(e) 晶粒体积分数(91.6%)(f) 黏结剂平均厚度与晶粒体积分数的关系
(e) The crystal volume fraction(91.6%)
(f) The relationship between the average thickness of the binder layer and the crystal volume fraction
(g) 晶粒的级配 (h) 粒径与平均粒径的对比
(g) The gradation of crystals (h) Comparison between the average crystal size and the actual crystal size
图4 不同体积分数的PBX二维几何模型
Fig. 4 The 2D geometric models for PBX with different crystal volume fractions
2 升温和降温过程中PBX及其界面的力学响应数值模拟
2.1 各组分材料参数
PBX 9502主要成分是TATB炸药晶粒[16-17]和高分子黏结剂F2314[18]．其中，TATB(三氨基三硝基苯)是一种具有各向异性的热力学属性的炸药颗粒．根据温茂萍等[19]和Sun等[20]的实验研究成果，认为TATB为弹塑性材料，其正交各向异性热膨胀系数等相关参数如表2所示．
表2中ρ代表密度，E为弹性模量，μ为Poisson比，σ为屈服极限，λ为导热系数，a11,a22,a33是正交各向异性的热膨胀系数，c为TATB的比热容．
表2 TATB热力学参数
Table 2 Thermodynamic parameters of TATB
ρ/(kg·m-3)E/MPaμσ/MPaλ/(W·m-1·K-1)1 9408 8600.350.544a11/K-1a22/K-1a33/K-1c/(J·kg-1·K-1)1.56×10-42.07×10-57.38×10-51 170
黏结剂F2314的材料参数如表3和图5所示．黏结剂为黏弹塑性材料，表3中A,n,m,f 为黏弹塑性材料参数，为拟合PBX空心半球顶点变形实验结果[21]所得．由黏结剂相关实验研究[19]得到了图5所示材料参数为黏结剂弹性模量与热膨胀系数随温度的变化情况．如图所示，黏结剂的弹性模量与热膨胀系数均随温度升高而不断降低．
表3 F2314热力学参数
Table 3 Thermodynamic parameters of F2314
ρ/(kg·m-3)μσ/MPaλ/(W·m-1·K-1)c/(J·kg-1·K-1)2 0300.30.10.15470Anmf1×10-80.435-0.580.82
(a) 弹性模量 (b) 热膨胀系数
(a) Young’s modulus (b) The thermal expansion rate
图5 F2314的弹性模量与热膨胀系数
Fig. 5 Young’s modulus and the thermal expansion rate of F2314
为探究炸药晶粒与黏结剂之间界面的脱黏行为，引入界面双线性本构模型，即应力-分离模型[22]，如图6所示．图中的d是界面张开位移，η是界面上包括法向与两个切向的应力，K是应力加载初始阶段的界面初始刚度，整个三角形的面积为断裂能．如图6所示，在AB段的应力加载时，界面的刚度保持不变，界面张开位移随应力的增大而增大．当界面所受应力不断增大，同时张开位移达到d0后，界面刚度开始下降，界面开始发生损伤，可以认为此处界面开始脱黏．随着界面张开位移继续增加，应力会沿着BC段下降．
图6 界面双线性本构模型
Fig. 6 The bilinear cohesive contact relationship model
对于界面刚度下降系数，采用如下算法[22]：
(1)
Kt=(1-D)K,
(2)
D为刚度下降系数，Kt为界面残余刚度值．当界面张开位移达到或超过dt时，刚度下降系数等于1，此时可认为界面完全损坏，其残余刚度为0．
通过PBX的拉剪复合实验等[23-24]得到了炸药晶粒-黏结剂界面相关参数，如表4所示．表中下标11，22，33分别表示法向与两个切向方向，G为界面断裂能．
表4 界面材料参数
Table 4 Parameters of the crystal-binder interface
K11/(MPa·m-1)K22/(MPa·m-1)K33/(MPa·m-1)T11/MPaT22/MPaT33/MPaG/(J·m-2)1061061060.7590.970.97100
2.2 体积分数72.4%的PBX在升温时的界面损伤
本文所用数值模拟方法是将PBX模型的左侧与下侧进行位移约束并由四周向内传热，用10 min将其从219 K(-54 ℃)升温到347 K(74 ℃)．温度随时间变化的曲线如图7(a)所示．由于模型尺寸较小，在升温过程中模型中心温度与表面温度温差较小，可以忽略内外温差．
(a) 升温时的温度曲线 (b) 界面刚度下降系数
(a) The temperature curve in a heating process (b) The interfacial stiffness reduction rate
(c) 界面法向应力云图 (d) 界面切向应力云图
(c) The normal stress nephogram of interface(d) The tangential stress nephogram of interface
(e) 界面最大应力与损伤初始应力的比值
(e) The ratio of the maximum interface stress to the initial damage stress
图7 升温过程中的界面状态
Fig. 7 Condition of PBX interface during heating process
图7(b)是在升温过程中PBX界面刚度下降系数云图．从图中可以看出，炸药晶粒与黏结剂之间的界面刚度下降系数为零，即在升温过程中界面刚度没有下降，晶粒与黏结剂的界面未发生损伤．图7(c)为界面法向拉应力最大时的整体云图．从图中可以看出，界面整体颜色为黑色与蓝色．其中黑色对应的是法向压应力，界面在受压应力时不会产生损伤，而彩色代表的拉应力最大值2.3×102 kPa(0.23 MPa)也远远达不到实验[23]所得的0.79 MPa，所以无法使界面脱黏．可以认为在升温中，黏结剂阻碍了TATB晶粒的膨胀，因此在大部分晶粒与黏结剂的界面上存在压应力．
图7(d)表示的是界面受到的剪应力最大时的整体云图．从该图中可看出，界面剪应力最大值在为PBX的右边界处，内部界面的剪应力不高，同样无法造成界面脱黏．图7(e)为整个升温过程中，界面所受最大应力与表4中设置的最大应力之比，可以看到，比值的最大值0.5出现在模型右边界上，对应图7(d)．根据图7(c)和图7(d)，同样说明在整体界面上均未出现脱黏现象，而界面所受的最大剪应力比最大拉应力更大，因此若界面在升温时发生损伤，可能由剪应力引起．
2.3 体积分数72.4%的PBX在降温时的界面损伤
对PBX模型的降温模拟与升温模拟相似，将炸药模型向四周散热，从347 K(74 ℃)经历10 min降温至219 K(-54 ℃)．图8(a)为模型降温时的整体温度曲线．
(a) 降温时的温度曲线 (b) 界面刚度下降系数
(a) The temperature curve in a cooling process (b) The interfacial stiffness reduction rate
(c) X方向位移云图 (d) Y方向位移云图
(c) The nephogram of the X-direction displacement (d) The nephogram of the Y-direction displacement
(e) 界面法向应力云图 (f) 界面切向应力云图
(e) The normal stress nephogram of interface(f) The tangential stress nephogram of interface
图8 降温过程中的界面状态
Fig. 8 Condition of PBX interface during cooling process
图8(b)是PBX在降温时的刚度下降系数达到最大时的云图．从图中可以看到，界面刚度下降的主要区域为大晶粒周围与黏结剂之间的界面．相较于升温时的界面刚度没有发生下降，可以认为在降温阶段，PBX的界面更容易发生脱黏的损伤，这与PBX热冲击实验[12]的现象相符合．图8(c)、8(d)是相应时刻PBX在X和Y方向的位移云图，可以看出，PBX整体在横向的位移在炸药大晶粒与左右的黏结剂层是不连续的，也能从一定程度反映出晶粒与黏结剂之间的界面脱黏情况．图8(e)是界面的刚度下降系数达到最大时的界面法向应力云图，其中负值为拉应力．从图中可以看出，拉应力最大，即颜色为红色的区域大部分集中在TATB大晶粒周围．结合界面刚度下降系数达到最大时的界面切向应力云图(图8(f)与图8(b))，可以认为界面的刚度下降和脱黏主要是由法向拉应力引起的．
2.4 体积分数对PBX降温时界面损伤的影响
实际PBX中的TATB含量可随着工况而发生改变，因此探究不同体积分数的PBX在降温时的界面损伤情况也十分重要．
图9(a)～9(d)为体积分数76.0%，81.7%，85.3%和91.6%的PBX模型在降温过程中界面刚度下降系数最大时的界面法向应力云图．图10(a)～10(d)为对应体积分数在降温过程中界面刚度下降系数最大时的界面切向应力云图．
(a) 76.0%的界面法向应力分布 (b) 81.7%的界面法向应力分布
(a) The interfacial normal stress of 76.0%(b) The interfacial normal stress of 81.7%
(c) 85.3%的界面法向应力分布 (d) 91.6%的界面法向应力分布
(c) The interfacial normal stress of 85.3%(d) The interfacial normal stress of 91.6%
图9 不同体积分数的PBX降温时的界面法向应力分布
Fig. 9 Interfacial normal traction of models with different crystal volume fraction during cooling process
对比图9与图8(e)可以发现，随着炸药晶粒的体积分数不断提高，界面上达到最大拉应力的区域逐渐减少，即由拉应力产生的界面脱黏现象不断减少．同时，对比图10与图8(f)可以发现，随着体积分数的提高，界面剪应力的最大值也有下降的趋势，从界面剪应力分布来看，大部分区域的剪应力不高，无法使界面产生脱黏的情况．
图11(a)为五种体积分数的PBX对应的整体界面上出现的最大刚度下降系数．根据图11(a)与式(2)可算得PBX界面最小残余刚度．图11(b)为不同体积分数对应的PBX界面上最小的残余刚度．结合图8～11，可以看出随着炸药晶粒体积分数的提高，界面上的最大刚度下降系数不断上升，整体界面的脱黏情况不断减少．
(a) 76.0%的界面切向应力分布 (b) 81.7%的界面切向应力分布
(a) The interfacial tangential stress of 76.0% (b) The interfacial tangential stress of 81.7%
(c) 85.3%的界面切向应力分布 (d) 91.6%的界面切向应力分布
(c) The interfacial tangential stress of 85.3% (d) The interfacial tangential stress of 91.6%
图10 不同体积分数的PBX降温时的界面剪应力分布
Fig. 10 Interfacial tangential traction of models with different crystal volume fraction during cooling process
(a) 最大刚度下降系数与体积分数的关系 (b) 界面最小残余刚度与体积分数的关系
(a) The relationship between the maximum stiffness (b) The relationship between the minimum residual stiffness reduction rate and the crystal volume fraction of interface and the crystal volume fraction
图11 降温后界面刚度与晶粒体积分数的关系
Fig. 11 Interface stiffnessess of different crystal volume fractions after a cooling process
2.5 炸药晶粒的级配对PBX降温时界面损伤的影响
在制备TATB晶粒时，可按需求给定不同的晶粒级配．本文针对几种体积分数相近但级配不同的PBX，探究其在降温时的晶粒-黏结剂界面的损伤情况．
图12～14为体积分数相近，各炸药晶粒级配不同的PBX在降温时界面产生损伤的数值分析．
(a) 几何模型 (b) 晶粒的级配
(a) The geometric model (b) The gradation of crystals
(c) 界面刚度下降系数 (d) X方向位移云图
(c) The interfacial stiffness reduction rate(d) The nephogram of the X-direction displacement
(e) 界面法向应力云图 (f) 界面切向应力云图
(e) The normal stress nephogram of interface (f) The tangential stress nephogram of interface
图12 各晶粒粒径相近的体积分数73%的PBX在降温过程中的界面状态
Fig. 12 Conditions of the PBX interface with a 73% crystal volume fraction and similar crystal sizes during a cooling process
(a) 几何模型 (b) 晶粒的级配
(a) The geometric model (b) The gradation of crystals
(c) 界面刚度下降系数 (d) X方向位移云图
(c) The interfacial stiffness reduction rate(d) The nephogram of the X-direction displacement
(e) 界面法向应力云图 (f) 界面切向应力云图
(e) The normal stress nephogram of interface (f) The tangential stress nephogram of interface
图13 体积分数71.7%的有一较大晶粒的PBX在降温过程中的界面状态
Fig. 13 Conditions of the PBX interface with a 71.7% crystal volume fraction and a large crystal during a cooling process
对比图12与图13可发现，当炸药晶粒的粒径相近时，整体在横向的位移更为连续(图12(d))，晶粒-黏结剂界面上的刚度下降系数更小(图12(c))，即界面损伤更小．
将图13的较大晶粒进一步放大，即得到图14(a)的PBX二维模型．从级配图(图14(b))可看出，该模型的最大晶粒的粒径约120 μm．以此为PBX模型进行降温数值分析，通过对比图13和图14可明显看出，当模型中最大晶粒更大时，整体在横向的位移在该晶粒处更不连续(图14(d))，对应的界面最大刚度下降系数也更高(图14(c))，在降温时界面上出现的损伤更严重．
(a) 几何模型 (b) 晶粒的级配
(a) The geometric model (b) Gradation of crystals
(c) 界面刚度下降系数 (d) X方向位移云图
(c) The interfacial stiffness reduction rate(d) The nephogram of the X-direction displacement
(e) 界面法向应力云图 (f) 界面切向应力云图
(e) The normal stress nephogram of interface (f) The tangential stress nephogram of interface
图14 体积分数71.9%的有一更大晶粒的PBX在降温过程中的界面状态
Fig. 14 Conditions of the PBX interface with a 71.9% crystal volume fraction and a larger crystal during a cooling process
结合图12～14可以看出，炸药晶粒的级配对PBX在降温过程中产生的损伤有很大的影响．从数值分析中可看出，在晶粒体积分数越相近时，各晶粒的粒径越趋于一致，整体的界面损伤会更小，同时那些粒径较大的晶粒周围也更易产生损伤．
3 结 论
本文对含能晶粒和黏结剂构成的含微结构的PBX进行了数值分析，重点研究了升温和降温过程对炸药晶粒与黏结剂界面脱黏的影响，主要结论如下：
1) 以Voronoi理论为基础的PBX几何建模有助于表达炸药晶粒与黏结剂构成的微结构，通过控制各晶粒的形状与级配，能生成与PBX真实结构更为相近的二维几何模型．
2) 升温时PBX界面会产生较大剪切应力，更易导致界面脱黏．在降温过程中，界面上的法向拉应力较剪应力更大，是界面脱黏的主导因素．相同体积分数下的升温与降温的数值结果检验了PBX在降温过程中更容易发生界面脱黏．
3) 随着炸药晶粒体积分数的提高，降温过程中界面的残余刚度就越高，表明提高晶粒体积分数有助于减少界面脱黏．
4) 在相近的体积分数下，炸药晶粒的级配对降温时的界面损伤有较大影响，晶粒粒径越相近，则在降温时界面损伤越小．
[1] WU Y, HUANG F. A micromechanical model for predicting combined damage of particles and interface debonding in PBX explosives[J]. Mechanics of Materials, 2009, 41(1): 27-47.
[2] WANG G, WANG Y, WEN Q. Thermal-mechanical analysis for confined HMX-based polymer-bonded explosives[J]. Journal of Thermal Stresses, 2019, 42(8): 1011-1034.
[3] HU W, WU Y, HUANG F, et al. Numerical simulation analyses of β↔δ phase transition for a finite-sized HMX single crystal subjected to thermal loading[J]. RSC Advances, 2018, 8(44): 24873-24882.
[4] KLEIN R. Voronoi Diagrams and Delaunay Triangulations[M]. New York: Springer, 1975.
[5] VORONOI G. Nouvelles applications des paramètres continus à la théorie des formes quadratiques, deuxième mémoire: recherches sur les parallélloèdres primitifs[J]. Journal für die Reine und Angewandte Mathematik (Crelles Journal), 1907, 134: 97-178.
[6] BARUA A, KIM S, HORIE Y, et al. Ignition criterion for heterogeneous energetic materials based on hotspot size-temperature threshold[J]. Journal of Applied Physics, 2013, 113(6): 64906.
[7] AMBOS A, WILLOT F, JEULIN D, et al. Numerical modeling of the thermal expansion of an energetic material[J]. International Journal of Solids and Structures, 2015, 60/61(4): 125-139.
[8] AMBOS A, TRUMEL H, WILLOT F, et al. A fast Fourier transform micromechanical upscaling method for the study of the thermal expansion of a TATB-based pressed explosive[C]//The 15th International Detonation Symposium. San Francisco, USA, 2014.
[9] MCGRANE S D, ASLAM T D, PIERCE T H, et al. Temperature of shocked plastic bonded explosive PBX 9502 measured with spontaneous Stokes/anti-Stokes Raman[J]. Journal of Applied Physics, 2018, 123(4): 045902.
[10] PI Z, LANG C, WU J. Temperature-dependent shock initiation of CL-20 based high explosives[J]. Central European Journal of Energetic Materials, 2017, 14(2): 361-374.
[11] 林聪妹, 刘佳辉, 曾贵玉, 等. 苯乙烯共聚物改性TATB基PBX的抗热冲击性能[J]. 含能材料, 2016, 24(2): 149-154.(LIN Congmei, LIU Jiahui, ZENG Guiyu, et al. Thermal shock resistance of styrene copolymer modified TATB-based polymer bonded explosive[J]. Chinese Journal of Energetic Materials, 2016, 24(2): 149-154.(in Chinese))
[12] 韦兴文, 吴束力, 唐兴. HMX基PBX炸药热损伤的数值计算与实验研究[J]. 火炸药学报, 2014, 37(4): 9-13.(WEI Xingwen, WU Shuli, TANG Xing. Numerical calculation and experimental study on thermal damage of HMX based polymer bonded explosive[J]. Chinese Journal of Explosives & Propellants, 2014, 37(4): 9-13.(in Chinese))
[13] WILLEY T M, LAUDERBACH L, GAGLIARDI F, et al. Comprehensive characterization of voids and microstructure in TATB-based explosives from 10 nm to 1 cm: effects of temperature cycling and compressive creep[C]//The 14th International Detonation Symposium. Coeur d’Alene, USA, 2010.
[14] 张伟斌, 田勇, 温茂萍, 等. JOB-9003炸药热冲击损伤的超声波检测[J]. 含能材料, 2004, 12(2): 85-88.(ZHANG Weibin, TIAN Yong, WEN Maoping, et al. Experimental study on the thermal shock damage of explosive by ultrasonic testing[J]. Chinese Journal of Energetic Materials, 2004, 12(2): 85-88.(in Chinese))
[15] 柏巍, 彭刚. 蒙特卡洛法生成混凝土随机骨料模型的ANSYS实现[J]. 石河子大学学报(自然科学版), 2007, 25(4): 504-507.(BAI Wei, PENG Gang. ANSYS implementation of Monte Carlo method for generating random concrete aggregate model[J]. Journal of Shihezi University (Natural Science), 2007, 25(4): 504-507.(in Chinese))
[16] BUECHLER M A, MILLER N A, LUSCHER D J, et al. Modeling the effects of texture on thermal expansion in pressed PBX 9502 components[C]//ASME 2016 International Mechanical Engineering Congress and Exposition. Phoenix, USA, 2016.
[17] THOMPSON D G, BROWN G W, OLINGER B, et al. The effects of TATB ratchet growth on PBX 9502[J]. Propellants Explosives Pyrotechnics, 2010, 35(6): 507-513.
[18] XU X, XIAO J, HUI H, et al. Molecular dynamic simulations on the structures and properties of ε-CL-20(0 0 1)/F2314 PBX[J]. Journal of Hazardous Materials, 2010, 175(1): 423-428.
[19] 温茂萍, 唐维, 董平, 等. 粘结剂含量对热压TATB基PBX残余应力的影响[J]. 含能材料, 2017, 25(8): 661-666.(WEN Maoping, TANG Wei, DONG Ping, et al. Effect of binder content on residual stress of thermally compacted TATB based PBX[J]. Chinese Journal of Energetic Materials, 2017, 25(8): 661-666.(in Chinese))
[20] SUN J, KANG B, ZHANG H, et al. Investigation on irreversible expansion of 1, 3, 5-triamino-2, 4, 6-trinitrobenzene cylinder[J]. Central European Journal of Energetic Materials, 2011, 8(1): 69-79.
[21] 唐维, 李明, 张丘, 等. PBX部件机械加工过程中的夹持变形预测[J]. 含能材料, 2008, 16(6): 703-707.(TANG Wei, LI Ming, ZHANG Qiu, et al. Prediction for clamping deformation of PBX parts on machining process[J]. Chinese Journal of Energetic Materials, 2008, 16(6): 703-707.(in Chinese))
[22] Dassault Systems. ABAQUS 6.14 documentation[DB/CD]. Providence, Rhode Island, USA, 2014.
[23] 颜熹琳, 唐明峰, 甘海啸, 等. 拉剪复合试验测试炸药晶体/粘结剂界面力学特性[J]. 含能材料, 2016, 24(6): 587-591.(YAN Xilin, TANG Mingfeng, GAN Haixiao, et al. Mechanical properties of explosive crystal/binder interface based on tension-shear test[J]. Chinese Journal of Energetic Materials, 2016, 24(6): 587-591.(in Chinese))
[24] 黄西成, 李尚昆, 魏强, 等. 基于XFEM与Cohesive模型分析PBX裂纹产生与扩展[J]. 含能材料, 2017, 25(8): 694-700.(HUANG Xicheng, LI Shangkun, WEI Qiang, et al. Analysis of crack initiation and growth in PBX energetic material using XFEM-based Cohesive method[J]. Chinese Journal of Energetic Materials, 2017, 25(8): 694-700.(in Chinese))