1 引 言

在细观尺度上，高聚物粘结炸药（Polymer Bonded Explosive, PBX）可以看成由形状不规则、尺寸满足一定级配关系、位置分布随机的高能炸药颗粒，高聚物粘结剂基体，以及炸药颗粒‑基体界面组成的多相复合材料，其炸药颗粒的体积分数都比较高，常高达90%以上。PBX的力学性能对其使用寿命、安全性和作用可靠性具有重要影响。数值模拟是研究PBX力学性能的重要手段，其中有效的PBX细观结构模型的建立，是在细观尺度上开展PBX力学行为数值模拟研究时首先面临的一个关键问题。

从文献来看，韦兴文^[1]、贾宪振等^[2]采用蒙特卡洛方法和随机投放法，分别生成了PBX随机圆形颗粒细观模型，由于圆形颗粒的限制，其得到的细观模型的颗粒体积分数均难达到85%以上；戴开达^[3]、傅华^[4]等人采用规则多边形颗粒替代圆形颗粒，在一定程度上提高了所生成PBX细观结构模型的颗粒体积分数，但其没有考虑颗粒间的级配关系，模型中多边形颗粒的尺寸和位置分布接近均匀，不符合真实PBX的细观结构特征；Seidel^[5]、Wu^[6]、Guo等^[7]采用Voronoi多边形颗粒建立PBX细观结构模型，与采用规则多边形颗粒相比，他们的工作保证了颗粒位置分布的随机性，并同时得到了较高的体积分数，但生成的Voronoi多边形颗粒尺寸均匀，因此也不符合真实PBX的细观结构特征；Barua和Zhou^[8]基于PBX的显微图像和高斯分布，分别建立了与PBX细观结构高度类似的二维细观结构模型，虽然其中的颗粒形状、尺寸及位置分布随机性很强，但从三维细观结构转换为二维细观模型所带来的误差，使其所得PBX细观结构模型的颗粒体积分数最高只达到了82%，与真实PBX中的颗粒体积分数相差甚远。

为了建立颗粒体积分数较高，并且颗粒形状、尺寸和位置分布特征均与真实PBX较为接近的PBX细观结构计算机模型，本研究在利用Voronoi方法^[9]生成多边形颗粒时，引入蒙特卡洛方法^[10,11]中的级配概念，形成了一种可以考虑级配分布的Voronoi多边形颗粒生成方法，进而通过颗粒修正、收缩和平滑等技术手段，形成了一种建立PBX细观结构模型的新方法。通过该方法建立的PBX细观结构模型的体积分数高达94%及以上，且同时还可以保证模型在边界处具有周期连续性。在此基础上，采用数值流形方法(Numerical Manifold Method，NMM)^[12,13]对不同体积分数和不同颗粒级配条件下的PBX细观模型的单轴压缩等效模量进行了计算和分析。

2 PBX细观结构计算模型生成算法

2.1 PBX细观结构计算模型生成流程

PBX细观结构计算模型的生成过程如图1所示。首先由蒙特卡洛方法生成给定级配的随机圆^[10]，然后以圆心作为Voronoi多边形的中心点生成Voronoi多边形，并将每个Voronoi多边形的顶点按逆时针顺序排序，进而对Voronoi多边形进行修正以保持级配关系并消除细长失真颗粒。在此基础上，根据指定的颗粒体积分数对Voronoi多边形进行收缩，生成颗粒间的粘结层，并对颗粒进行切角和光滑处理等操作，得到PBX细观结构几何模型。最后，将细观几何模型形成ANSYS命令流，通过ANSYS进行网格划分，最终得到用于数值模拟的PBX细观结构计算模型。

图1 PBX细观结构计算模型生成过程

Fig.1 The generation procedure of meso‑scale structure computational model for PBX

2.2 考虑级配的Voronoi多边形颗粒模型的生成

蒙特卡洛方法生成圆形颗粒时，需指定圆形颗粒的尺寸范围、级配个数和各级配颗粒尺寸等特征参数，采用转化后的二维Fuller公式^[11]计算出各级颗粒在投放区域所占比例，进而计算出各级颗粒的个数。最后，按均匀分布的方式，确定圆形颗粒的位置，完成圆形颗粒的蒙特卡洛方法生成过程。

根据文献[8]所采用的颗粒尺度，将颗粒尺寸从小到大分为3个等级，即3个级配，投放区域的大小为1.5 mm×1.5 mm，表1给出了三种不同级配类型下的颗粒尺寸范围及各级配颗粒尺寸。图2所示为三种级配类型下各生成的一组随机圆形颗粒。级配类型的不同，得到的圆形颗粒所占的体积分数也不同，随着级配中大颗粒尺寸的增加，颗粒的体积分数随之增加。图2中，颗粒体积分数最高仅达到62%，远低于90%，这也是不直接采用这种模型作为PBX细观模型的主要原因。

图2 蒙特卡洛方法生成的圆形颗粒模型

Fig.2 The circular particle models generated by the Monte Carlo method

为了提高颗粒体积分数，对于以上生成的圆形颗粒模型，将其中的圆心视为投放区域内的一组随机分布点，并以这些随机分布点为Voronoi多边形的中心，进一步生成相应的Voronoi多边形模型。为了使得到的Voronoi多边形模型具有周期性边界条件，在投放区域内，每生成一个随机点，同时也在周边8个邻域分区内相对应的另外8个位置额外生成8个点。以图中的模型为例，生成的周期性点的分布情况如图3a所示，相应生成的Voronoi多边形模型如图3b所示。

图3 周期性点及投放区域内的Voronoi多边形模型

Fig.3 The periodically distributed points and the Voronoi polygon model in launch area

polygon model

2.3 Voronoi多边形颗粒模型的修正

从图3b可以看到，对应于图2a中的随机圆形颗粒，以3个大圆的圆心（A、B、C）为中心的Voronoi多边形颗粒，周围布满了细长的Voronoi多边形。根据Voronoi多边形的生成算法，介于大圆圆心和小圆圆心之间的Voronoi多边形的边界，是这两个圆心连线的垂直平分线，所以这些细长的多边形是由大圆和小圆之间的尺寸突变引起的。为了消除这种细长颗粒，且同时使得模型中颗粒间保持指定的级配关系，需要对生成的Voronoi多边形进行修正。如图3b所示，筛选出相应面积变化较大的Voronoi多边形颗粒作为修正对象，比如圆心A所在的Voronoi多边形颗粒，以点A为圆心，用点A所对应的图2a中的圆截切与之相交的Voronoi多边形，删除圆内的线段，形成新的Voronoi多边形颗粒。按照这种修正方法，图3b所示的Voronoi多边形颗粒模型经修正后如图4所示。一般，修正后的颗粒模型在边界处仍然保持有周期连续性。

图4 修正后的Voronoi多边形颗粒模型

Fig.4 The modified Voronoi polygon particle model

图5所示为修正后的Voronoi多边形颗粒模型与相应的蒙特卡洛随机圆形颗粒模型中的颗粒尺度统计规律对比，其中横轴为模型中颗粒按尺度由小到大的排列序号。图5中的横坐标显示，蒙特卡洛随机圆形颗粒有180个，修正后得到的Voronoi多边形颗粒有220个，颗粒数目明显增多，增多的颗粒填充了原随机圆形颗粒间的空隙。从图5中还可以看出，在修正后的Voronoi多边形颗粒模型中，粒径接近0.07 mm的颗粒数量有所增多，且其粒径有缓慢增大的趋势，而中等粒径及大粒径的颗粒的尺寸却有所减小，整个模型中颗粒的级配关系仍较为明显，与原随机圆形颗粒的级配相对较一致。

图5 修正后的Voronoi多边形颗粒模型与相应蒙特卡洛随机圆形颗粒模型间的颗粒级配比较

Fig.5 Particle gradation comparison between the modified Voronoi polygon particle model and the corresponding Monte Carlo circular particle model

2.4 PBX细观结构计算模型的建立

修正后的Voronoi多边形颗粒细观模型的颗粒体积分数实际为100%。为了得到指定的体积分数，同时形成颗粒间的粘结剂层，需要对多边形颗粒进行收缩处理。个别不规则的颗粒在收缩后还可能产生尖角，为了避免其对随后网格划分质量的影响，对小于给定角度阀值的角进行切削，并规定切削后形成的短边的边长不能小于多边形颗粒收缩的距离。图6所示为模型经颗粒收缩、切角后，得到的一组具有不同体积分数的PBX细观结构模型，其中颗粒间的空白区域为粘结剂层。

图6 一组不同颗粒体积分数的PBX细观结构模型

Fig. 6 A group of meso‑scale structure models for PBX with different particle volume fractions

对于图6所示的细观结构模型，给出了模型的颗粒体积分数随粘结剂厚度的变化关系，如图7所示。从图7可以看出，随粘结剂厚度的增加，细观模型的体积分数按接近线性的规律减小。一般地，在生成粘结层时，模型中的个别细小颗粒经过收缩和切角后，可能会出现极其细长的形状，为了避免这些颗粒对网格划分和模型计算性能的影响，对其予以删除。当删除发生在模型边界附近时，这会导致模型的周期连续性变弱。此外，在收缩、切角后，模型中多边形颗粒的某些边，常会出现边长极小的情况，这也会影响模型的网格划分。为此，收缩、切角后，进一步在多边形颗粒的每条边上插入若干点，通过插入点对多边形进行样条曲线平滑处理。将颗粒平滑后的PBX细观结构模型形成ANSYS命令流，通过ANSYS软件建模并划分网格，建立PBX的细观结构计算模型。图8所示为通过以上方法最终建立的具有不同级配类型的三个PBX细观结构计算模型示例。

图7 颗粒体积分数随粘结剂厚度的变化关系

Fig.7 Relationship between the particle volume fraction and the matrix thickness

图8 PBX细观结构计算模型

Fig.8 The computational models of PBX meso‑scale structure

3 细观模型等效弹性模量的NMM模拟

数值流形方法(NMM)^[12,13]通过采用双层覆盖技术，将有限元、非连续变形分析（DDA）和解析法统一在同一个理论框架下，能够方便的实现带裂纹的结构、多物体连续与非连续大变形与大位移，以及多物理场的耦合等复杂问题的求解。自20世纪90年代诞生以来，NMM在理论研究和不连续变形及裂纹扩展分析应用等方面得到了较快的发展^[14]。

本研究采用NMM方法计算单轴压缩条件下PBX细观模型的等效弹性模量，以验证本文所提出的PBX细观结构计算模型建立方法的有效性，并初步探索颗粒体积分数及颗粒级配对PBX等效模量的影响。模拟中炸药颗粒采用各向同性弹性本构，粘结剂采用22单元的Prony级数粘弹性本构^[6,8]，炸药颗粒与粘结剂之间的界面采用双线性接触损伤模型^[15,16]。各组分的材料参数取值与文献[3]相同，具体见表2。

表2 PBX各组分的材料参数

Table 2 Material parameters of each component in PBX

parameters	elastic modulus / GPa	Poisson′s ratio	density / kg·m-3	tensile strength / MPa	friction angle / (°)	cohesive strength / MPa
particle	8.86	0.322	1910	10	30	30
binder	0.397	0.495	1280	6	25	6
interface	-	-	-	1.66	15	1.66

根据表1中的三种不同级配参数，建立三组不同级配类型的PBX细观结构计算模型。在每一组模型的生成过程中，通过改变Voronoi多边形颗粒收缩时的缩进量，分别得到十个不同颗粒体积分数的计算模型。其中不同级配类型的三个体积分数接近85%的PBX细观结构模型如图8所示。三组细观模型的总体尺寸均为1.5 mm×1.5 mm。模型底部设置一组固定点，约束底部Y方向的位移，而其沿X方向自由。模型上端设置一组位移加载点，沿Y方向向下加载，加载速度4 mm·s^-1。模型左端外侧设置一固定板，通过固定板限制细观模型左端沿X方向上的位移，而细观模型左端沿Y方向自由。为保证细观模型右边界在X方向上的位移始终相等，在模型右边界设置一耦合运动约束板，约束板在X方向自由，在Y方向不允许发生位移，细观模型右边界与约束板进行X方向的粘结，允许细观模型右边界与约束板间发生Y方向的相对位移。

在NMM模拟中，在每一个加载步计算得到单元的应力应变值后，按照文献[3,17-18]中的方法计算当前时刻的等效弹性模量。根据该方法，基于平面应变假设，细观模型上的平均应力和平均应变可以表示为：

式中,i表示应力/应变张量的方向，A表示二维PBX细观结构模型的面积。根据平面弹性应力‑应变关系，可以得到PBX在平面应变假设下的等效弹性模量和泊松比为：

再根据文献[19]中平面应变和平面应力之间的转换关系，得到PBX的等效泊松比和等效弹性模量最终表达式如下：

在单轴压缩模拟中，根据Y方向的加载位移换算得到模型的即时Y方向应变，结合式（3）中得到的当前等效弹性模量，就可以得到模型在Y方向上的总体等效应力‑应变关系。

图9所示为NMM模拟所得三组细观模型在Y方向的等效应力‑应变关系曲线。在三种不同的颗粒级配条件下，当Y方向上的应变值一定（即加载位移相同）时，随着颗粒体积分数的增大，Y方向上的应力都逐渐增大。这是由于当体积分数较小时，颗粒间较厚的粘结剂层相对于炸药颗粒起到了更多的变形作用，模型整体体现出粘结剂的“较软”性质，应力随应变的变化较缓慢；而当体积分数较大，如大于90%时，颗粒间的粘结层很薄，在压载作用下炸药颗粒起到了主要的支撑作用，所以模型整体体现出颗粒的“较硬”性质，此时应力应变曲线陡峭。

a. type Ⅰ

b. type Ⅱ

c. type Ⅲ

图9 不同颗粒体积分数下的应力‑应变曲线

Fig.9 Stress‑strain curves under different particle volume fractions

图10给出了三种颗粒级配条件下细观模型的等效模量随体积分数的变化情况，其中等效模量是图9中应力‑应变曲线初始线弹性段的斜率。在单轴压缩加载过程中，随着加载位移的增大，颗粒和基体界面上会出现损伤，从而导致等效应力‑应变关系逐渐呈现非线性。从图10可以看到，当体积分数较小时，等效弹性模量的增速较为缓慢，此时粘结剂起到较大的变形作用；而当体积分数较大时，等效弹性模量的增速变大，这表明炸药颗粒起到的支撑作用越趋明显。这一现象和文献[3]中的结果吻合较好。对于Ⅰ，Ⅱ，Ⅲ三组不同级配的模型，当体积分数低于30%时，三种颗粒级配条件下得到的等效模量基本相等，这是因为此时粘结剂对等效模量的确定起到了支配作用；与此相似，当体积分数高于85%时，三种颗粒级配条件下得到的等效模量也基本相等，但这时炸药颗粒对等效模量的确定起到了支配作用；而当体积分数为30%~85%时，颗粒和粘结剂的支配作用都相对较弱，因此炸药颗粒的级配对等效模量产生了较明显的影响。

图10 单轴压缩下等效模量随颗粒体积分数的变化关系

Fig.10 Relationships between the equivalent modulus and the particle volume fraction under uniaxial compression

在文献[3]中，作者在对其数值模拟结果进行验证时，列出了PBX材料在体积分数分别为90%、93%、95%和97.5%时的一组弹性模量实验测量值，见表3。表3也给出了在这4个体积分数条件下本研究的NMM模拟结果，这些结果由图10中的曲线插值得到。以实验值为基础，NMM模拟结果与实验结果间的百分比偏差也列于表3。从表3中可以看出，NMM模拟结果和实验结果之间存在较大的偏差。例如，在体积分数为90%时，三种不同级配条件下，NMM模拟结果都较实验值偏大，最小偏差是9.21%。造成模拟结果和实验结果之间出现如此大偏差的因素有很多，比如模拟中所使用的各组分材料参数与实验中所使用材料的真实参数之间本身可能就存在较大差别，实验中一般难以精准测量出炸药颗粒的真实弹性模量。另外实验是针对三维模型，而模拟采用的是平面模型，二维与三维之间的转换，也可能会在很大程度上引起偏差。其中一个显著而重要的因素是，PBX材料在合成、成型、机械加工等过程中，会产生例如粘结剂包覆不全、炸药晶体内部含有孔穴、气泡及微裂纹等缺陷，这些缺陷会使得PBX宏观等效模量降低^[20]，而本文NMM模拟中采用的是一种理想化的PBX细观结构，没有考虑缺陷的存在，因此NMM模拟结果比实验值偏大，在一定程度上可以认为是合理的。在下一步的工作中，可以通过模拟来研究缺陷对PBX等效弹性模量的影响。

从表3还可看出，随着体积分数从90%逐渐增大到97.5%，在三种不同级配条件下，NMM模拟结果与实验值间的偏差都随体积分数的增大而增大，这也可以认为是由于真实的PBX材料中的缺陷是随着体积分数的增大而增加的原因造成的。

4 结 论

将蒙特卡洛方法中级配的概念，与Voronoi多边形颗粒生成算法相结合，形成了一种可以考虑级配分布的Voronoi多边形颗粒生成方法，通过颗粒修正、收缩、切角和平滑等处理，建立考虑级配和满足体积分数要求的PBX细观结构模型。得到的PBX细观模型的颗粒体积分数可达到94%及以上，同时还较好的保证了模型在边界处的周期连续性。进而，采用NMM方法对不同颗粒级配条件下的三组PBX细观模型进行了单轴压缩模拟。结果表明：

（1） PBX的等效弹性模量随颗粒体积分数的增加而增大，当体积分数较小时，等效弹性模量增速较为缓慢，而当体积分数较大时，等效弹性模量的增速较为显著；当体积分数低于30%及高于85%时，颗粒的级配对等效模量的影响较小，当体积分数为30%~85%时，颗粒级配对等效模量影响显著。以上现象归功于粘结剂或炸药颗粒在PBX变形过程中的支配作用程度，当两者的支配作用都不明显时，等效模量受炸药颗粒级配的影响就得以突显。

（2） 与实验结果相比，考虑到PBX各组分材料参数的精确性、二维与三维之间的转换，以及缺陷等因素，可以认为本文模拟所得到的PBX等效弹性模量是合理的，从而也验证了本文所提出的PBX细观结构模型建立方法的合理性。在下一步研究中，将考虑通过数值模拟来研究缺陷对PBX等效弹性模量的影响等问题。

本研究可为进一步开展PBX的细观力学，乃至热力学和起爆性能研究奠定基础。

参考文献