王 艺，黄西成

（中国工程物理研究院总体工程研究所，四川 绵阳 621999）

1 引言

高聚物黏结炸药（Polymer Bonded Explosive，PBX），在细观上可以看作由形状不规则、位置分布随机、尺寸满足一定级配关系的颗粒性主体炸药，高聚物黏结剂，以及炸药颗粒/黏结剂界面组成的多相非均质材料或混合物［1-2］，属高填充复合材料（highly⁃filled polymer composites）［3］。本研究选取PBX⁃9502 作为研究对象，其含有质量分数95%的三氨基三硝基苯（TATB）晶粒和5%的黏结剂［4］。

在实际使用过程中，由于昼夜、季节的变化，PBX材料会经受一系列的热循环加载，产生不可逆的体积增长，这种不可逆的体积增长现象被称作热棘轮增长［5］。在热作用下，炸药热胀冷缩产生宏观上的尺寸变化和各向异性体积膨胀，产生的药壳间隙则使得武器战斗部装药在振动、冲击、热加载环境中容易与壳体发生摩擦、碰撞，影响系统的安全使用［6］。国内外学者对PBX 材料的热棘轮增长现象开展了相关研究，发现热循环加载后PBX 炸药感度增加［7-9］，其主要原因是热膨胀过程中产生了新的孔隙或已有的孔隙发生变化［10］，例如Willey 等［11］发现Φ10×0.8 mm 的PBX⁃9502圆柱体经过热循环后空腔数量和尺寸增加，孔隙体积分数从1.6%增长至2.5%。

Thompson 等［12-13］通过实验探索PBX⁃9502 在热加载过程中出现的不可逆体积增长现象，发现其体积增长的幅度会随着周期的增加而逐渐减小，在发生体积增长的样品中，其拉伸断裂应力和弹性模量均较低，包含了更细的孔隙，孔隙的分布也更广，并且Mcgrane 等［14］发现PBX⁃9502 在不可逆变形过程中产生的孔隙主要是晶间裂纹。Schwarz 等［15］猜测在PBX⁃9502 中，不可逆变形主要是TATB 晶粒的各向异性所致，但是学术界对于影响PBX⁃9502 不可逆变形和孔隙变化的主要因素还缺乏清晰的认识，对机理的研究仍显薄弱。目前，在数值计算中很少考虑到材料的热学性能和TATB 晶粒的各向异性，缺乏明确的模型来预测PBX⁃9502 中的热棘轮增长现象。

认识和模拟多晶TATB 基复合材料的热膨胀是高能炸药力热模型研究中的一个重要方面，这种模拟可用于指导设计和操作决策制定，尤其是在武器系统部件的安全性和完整性方面。几十年来，TATB 复合材料在热膨胀方面的若干问题，一直是研究的重点和难点［16-17］。因此，开展PBX⁃9502 在热循环载荷下的不可逆变形机理研究，有助于深入认识TATB 基炸药的环境适应性，进而提高系统的安全性和可靠性。

本研究针对PBX⁃9502 的不可逆变形试验，根据材料中TATB 晶粒的随机分布特征，建立三相微结构模型，分别采用各向异性热弹性、热弹塑性、双线性本构描述TATB 晶粒变形、黏结剂损伤以及晶粒/黏结剂界面的破坏，利用扩展有限元（XFEM）法计算裂纹扩展。通过分析PBX⁃9502 在热循环加载过程中的黏结剂损伤和晶粒/黏结剂界面的破坏等对不可逆变形的影响，研究PBX⁃9502 不可逆变形的细观机理。

2 数值模拟

2.1 热棘轮增长实验的分析

Thompson 等［12］对PBX⁃9502 的热棘轮增长现象开展了相关测量，试件为Ф5X5 mm 的PBX⁃9502 圆柱体，在两个极端温度-43 ℃和113 ℃之间以交替加热和冷却的方式进行循环加载，每次循环的初始点和终点温度为室温（23 ℃），即一次热循环的加载历程为：从23 ℃开始，先冷却到-43 ℃，再加热到113 ℃，最后重新回到23 ℃，温升/降率为1 ℃·min-1。图1 为该试验测得的前10 个热循环过程中，每一次循环加载后的不可逆应变ε随循环次数n增加而变化的历程。

图1 PBX⁃9502 试件在-43 ℃至113 ℃之间热棘轮增长［12］Fig.1 Dimensional changes and irreversible growth ofPBX⁃9502 specimens under cyclic heating / cooling between-43 ℃and 113 ℃［12］

从图1 中的试验数据可知：不可逆应变ε与循环次数或周期n的关系近似于对数曲线：

随着循环次数n的增加，不可逆应变的增长趋于平缓。

2.2 PBX⁃9502 细观建模

研究PBX⁃9502 炸药的热棘轮增长现象，采用微结构表征技术并与其宏观力学性能相关联是一种非常有效的，也是非常重要的研究手段［18-20］。然而，PBX⁃9502 的微结构建模仍是一个挑战［21］，这种挑战在于如何精确表征晶粒级配，相材料性能等微结构特征。本研究利用Voronoi 多边形颗粒建立材料的微结构模型［2］，同时考虑了PBX⁃9502 中各向异性TATB 晶粒、黏结剂以及晶粒/黏结剂界面的热力学性能。

2.2.1 PBX⁃9502 微结构建模

PBX⁃9502 材料的主要成分为TATB 晶粒和黏结剂，在细观尺度上，可将之看作三相非均质复合材料：TATB 晶粒、黏结剂、晶粒/黏结剂界面，如图2所示［22］。

图2 PBX 的SEM 图像［22］Fig.2 SEM photograph of PBX［22］

根据材料的微结构特点，在蒙特卡洛方法中引入级配的概念，与Voronoi 方法结合来生成多边形颗粒，建立高颗粒体积分数的PBX⁃9502 微结构模型［2］。

由蒙特卡洛级配思想［23］生成给定级配的随机圆，然后以圆心作为Voronoi 多边形的中心点生成Vor⁃onoi 多边形，接着对Voronoi 多边形进行收缩，在颗粒间生成3 μm 厚的黏结层，最后对颗粒进行切角和光滑处理。设置的级配参数见表1。

表1 Voronoi 模型级配参数Table 1 Grading parameters of voronoi model

由该级配得到炸药颗粒体积占比为87% 的PBX⁃9502 微结构几何模型，如图3 所示。

图3 Voronoi 模型Fig.3 Voronoi model

由于TATB 晶粒具有高度的各向异性，在Voronoi模型中，以随机取向的方式给每个颗粒赋予局部的材料方向。

2.2.2 微结构材料本构模型

在PBX⁃9502 的细观建模过程中，还需设置各组分材料的材料性能。本研究考虑了PBX⁃9502 中各向异性TATB 晶粒、黏结剂以及晶粒/黏结剂界面的热力学性能，还考虑了黏结剂的拉伸破坏和晶粒/黏结剂界面的脱黏［14，24］。

2.2.2.1 TATB 材料各向异性弹性和各向异性热膨胀

张丘等［25］的研究表明，单个完好的TATB 晶体经热循环后不会出现不可逆应变，因此，可将TATB 晶体视为各向异性的线弹性体。各向异性线弹性体本构方程的一般形式为［26］：

式中，σij为应力张量，GPa；εkl为应变张量，为无量纲参数；Cijkl为四阶弹性张量，GPa；αkl为热膨胀系数张量，℃-1；ΔT为温度变化量，℃。

由弹性模量的局部对称及模量的主对称性，Cijkl=Cklij=Cjikl=Cijlk，四阶弹性张量Cijkl可通过Voigt 标记法简化为含21 个分量的二阶弹性张量Eij，GPa。同理，热膨胀系数张量αkl可通过Voigt标记法简化为含6个分量的向量βj，℃-1。采用Schlenker 方法［27］，通过晶格参数［28-29］即可求得TATB的弹性张量Eij和热膨胀系数βj。

2.2.2.2 界面模型

引入界面双线性本构模型，即线弹性应力⁃分离准则［30］来描述晶粒/黏结剂界面的脱黏，图4 为该模型在一维时的示意图。

图4 界面双线性本构模型［30］Fig.4 Interface bilinear constitutive model［30］

图4 中d 代表界面张开位移，m；σa为界面破坏应力，MPa；OA 边的斜率K，MPa·m-1，是界面初始刚度，三角形OAB 的面积为断裂能G，J·m-2。当界面受到外力加载时，加载路径为OA，若加载的应力首次达到σa，界面开始脱黏；继续加载，加载路径变为AB，界面刚度逐渐降低，且不可恢复，最后达到B 点处时刚度变为0，界面完全破坏。

在三维空间中，破坏应力σa包含了σ11，σ22，σ33三个分量，分别代表法向和两个切向的最大破坏应力；初始界面刚度K包含了K11，K22，K33三个分量，分别代表法向和两个切向的初始界面刚度。

实验测得［31］：界面的法向破坏应力σ11= 0.759 MPa，切向破坏应力σ22=σ33= 0.97 MPa，设置晶粒/黏结剂界面的材料参数如表2 所示。

表2 界面材料参数［31］Table 2 Material parameters of interface［31］

2.2.2.3 黏结剂热弹塑性模型

PBX⁃9502 的黏结剂为Kel⁃F 800，材料性质受温度影响较大，当达到特定温度时，出现玻璃化转化现象，材料属性急剧变化，在建模时采用温度相关的弹塑性模型，其杨氏模量与泊松比随温度的变化如表3所示［32］。

表3 黏结剂的杨氏模量与泊松比随温度的变化［32］Table 3 Young′s modulus and poisson′s ratio of binder with temperature［32］

黏结剂的塑性参数和拉伸强度见表4［32］。

表4 黏结剂的塑性参数和拉伸强度［32］Table 4 Plastic parameters and tensile strength of binder［32］

黏结剂的热膨胀系数也是一个很重要的材料参数，这会影响黏结剂与TATB 之间的互相作用以至于影响材料内应力的分布，设置黏结剂的热膨胀系数如表5 所示［33］。

表5 黏结剂的热膨胀系数［33］Table 5 Thermal expansion coefficient of binder［33］

2.3 仿真计算

设定仿真计算的条件：设定加载的起始点和终点温度为23 ℃，从23 ℃开始，先加热到103 ℃，再冷却到-57 ℃，最后重新回到23 ℃，以此作为一个周期的加载路径，温升/降率1 ℃·min-1，如图5，加载10 个周期；温度加载采用场变量形式，黏结剂的破坏采用基于最大主应力的启裂准则以及基于断裂能的裂纹张开准则，取黏结剂的断裂能为100 J·m-2［34］。

图5 一个周期的温度加载历程Fig.5 Temperature loading history for one cycle

计算中的试件假定为平面应力模型，竖直方向为轴向。利用2.2.1 中的方法建立了规模为0.5×0.5 mm 的Voronoi 模型，图3。TATB 晶粒的弹性模量和热膨胀系数分别采用公式3 中的Eij和βj；晶粒/黏结剂的界面设置见表2；黏结剂的弹性、塑性参数，拉伸强度和热膨胀系数见表3 至表5，利用扩展有限元法（XFEM）来计算裂纹的启裂及扩展；其他计算参数设置如表6所示［33］。

表6 黏结剂和TATB 的热力学参数［33］Table 6 Thermodynamic parameters of binder and TATB［33］

3 结果与分析

加载过程中，PBX⁃9502 试件（图3）中的应力云图如图6。加载结束后，PBX⁃9502 试件（图3）的损伤如图7，其中图7a 为加载结束后的模型；区域A 发生了黏结剂的开裂，见图7b；图7c～7e 分别为B 区域晶粒/黏结剂界面脱黏形成的裂缝、裂缝的法向间隙和切向滑移。

由图6、图7 可见，在温度加载过程中，PBX⁃9502试件内部产生了不均匀的应力场，当局部应力大于材料强度时，产生黏结剂的破坏和晶粒/黏结剂界面的脱黏。推测其原因为：TATB 晶粒的各向异性弹性以及各向异性热膨胀、TATB 晶粒与黏结剂热力学性能的差异，导致TATB 晶粒变形与黏结剂变形不协调。

图6 两个不同时刻PBX⁃9502 试件中的应力云图Fig.6 Stress field at two certain time of PBX⁃9502

图7 加载结束后PBX⁃9502 试件的裂纹和界面脱黏Fig.7 Cracking and debonding after loading of PBX⁃9502

图8 和图9 分别为加载过程中当温度为103 ℃和-57 ℃时PBX⁃9502 试件（图3）晶粒/黏结剂界面的法向间隙和切向滑移。由图8 和图9 可知：在加热到103 ℃时，试件中的界面几乎不存在法向间隙，仅存在切向的滑移。在冷却到-57 ℃时，大颗粒附近出现界面的法向间隙，切向的滑移距离也更大。结果表明：PBX⁃9502 在冷却阶段更容易出现界面的脱黏，受到较大的拉应力，这也意味着在冷却阶段更容易产生黏结剂开裂。

图8 晶粒/黏结剂界面的法向间隙和切向滑移（103 ℃）Fig.8 Normal gap and tangential slip of interface（103 ℃）

图9 晶粒/黏结剂界面的法向间隙和切向滑移（-57 ℃）Fig.9 Normal gap and tangential slip of interface（-57 ℃）

图10 为计算过程中颗粒尖端A 点处裂纹的萌生及演化过程。

裂纹起始点位于图10a 中的A 点所示颗粒尖端处；图10b 中裂纹起裂时温度为⁃15 ℃，此时PBX⁃9502 处于冷却阶段；裂纹的起点和终点都位于颗粒边界上。

图10 A 点处的裂纹萌生及演化Fig.10 Crack initiation and evolution at a certain point at A

数值计算中，每一次循环加载后PBX⁃9502 的轴向应变ε随循环次数n增加而变化的历程如图11 所示。由图11 可知：当黏结剂取弹性本构时，最大轴向应变接近0.17%，小于使用塑性本构时的0.2%，而且周期性的应变增长也更缓慢。这意味着黏结剂的塑性是PBX⁃9502 应变增长的来源之一。

图11 循环加载后PBX⁃9502 的轴向应变Fig.11 The axial strain at the end of cyclic loading of PBX⁃9502

Bennett 等［35］通过数值模拟也对PBX⁃9502 的不可逆变形进行了研究，该研究中，PBX⁃9502 在第10 个加载周期的不可逆应变约为1.20%，远大于实验测得的0.45%［12］和0.25%［9］。与之相比，本研究通过计算得到的不可逆应变为0.2%，更贴近实验结果。

目前为止，学术界对PBX⁃9502 的热循环不可逆变形现象的相关研究以宏观试验居多，积累了一定的理论基础，但是对机理的认识仍显薄弱，由于PBX⁃9502 微结构复杂，在细观尺度上的研究容易受到试验仪器和建模方法的限制。本研究引入带有级配的Voronoi 模型，采用各向异性热力学细观建模方法，通过扩展有限元计算PBX⁃9502 在温度循环加载过程中的应力不协调状态，将PBX⁃9502 的细观损伤与宏观不可逆变形结合起来，对于PBX⁃9502 热循环不可逆变形现象的机理研究，具有一定的指导性。然而，通过数值方法计算PBX⁃9502 的细观损伤是一种前沿且较为困难的问题，还处于理论研究阶段，许多问题尚未解决，有待进一步的探索和完善。

4 结论

（1）在温度场作用下，TATB 晶粒的各向异性弹性以及各向异性热膨胀、TATB 晶粒与黏结剂热力学性能的差异，导致晶粒变形与黏结剂变形不协调，产生非均匀应力场。当局部应力超过材料的拉伸强度时，PBX⁃9502 内部萌生裂纹，随着温度的循环，裂纹进一步扩展，最终导致PBX⁃9502 的不可逆变形。

（2）仿真计算显示：黏结剂的塑性也是PBX⁃9502不可逆变形的来源之一，而且对应变的周期性变化趋势会造成一定的影响。

（3）在冷却过程中，PBX⁃9502 更容易产生晶粒/黏结剂界面的脱黏。