康歌1，陈鹏万1，曾乙伦2，甯尤军2，3
（1.北京理工大学机电学院，北京100081；2.西南科技大学制造科学与工程学院，四川绵阳621010；3.西南科技大学工程材料与结构冲击振动四川省重点实验室，四川绵阳621010）
摘要：为了从细观尺度上对高聚物黏结（PBX）的力学性能进行模拟，将蒙特卡洛方法中级配的概念引入高体积分数的PBX 细观模型，与Voronoi方法相结合，形成了一种可以考虑级配分布的Voronoi多边形颗粒生成方法，通过颗粒修正、收缩、切角和平滑等处理，建立了考虑级配和体积分数达94%及以上的PBX细观结构模型，模型在边界处还可保持周期连续性。采用数值流形（NMM）方法对不同颗粒级配条件下的三组PBX细观模型进行了单轴压缩模拟。将NMM模拟所得等效弹性模量与文献实验结果进行了比较。分析了模拟结果与实验结果之间出现偏差的原因，讨论了PBX细观结构模型生成方法的合理性。结果表明，PBX的等效弹性模量随体积分数的提高而增大，并且体积分数越高，增速越快；同时，当体积分数低于30%及高于85%时，颗粒的级配对等效模量的影响较小，而当30%≤体积分数≤85%时，颗粒级配对等效模量影响显著。
关键词：复合材料；高聚物粘结（PBX）；蒙特卡洛方法；Voronoi方法；颗粒级配；等效模量
环境模拟舱中图分类号：TJ55；TB34文献标志码：A DOI：10.11943/CJEM2018023
1引言
在细观尺度上，高聚物粘结（Polymer Bond‑ed 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多边形
文章编号：1006⁃9941（2018）09⁃0772⁃07
收稿日期：2018⁃01⁃22；修回日期：2018‑04⁃26
网络出版日期：2018‑06‑14
基金项目：国家自然科学基金资助（51204137）；四川省科技计划项
目资助（2017JY0128）
作者简介：康歌（1989-），男，在读博士，主要从事含能材料力学性
塑料水嘴能研究。e‑mail：1208339811@qq
通信联系人：甯尤军（1981-），男，研究员，主要从事计算方法和非
连续介质力学行为研究。e‑mail：cnningyj@foxmail
引用本文：康歌，陈鹏万，曾乙伦，等.一种高颗粒体积分数PBX细观模型的生成方法［J］.含能材料，2018，26（9）：772-778.
KANG Ge，CHEN Peng⁃wan，ZENG Yi⁃lun，et al.A Method of Generating Mesoscopic Models for PBXs with High Particle Volume Fraction［J］.Chinese Journal of Energetic Materials（Hanneng Cailiao），2018，26（9）：772-778.
颗粒生成方法，进而通过颗粒修正、收缩和平滑等技术手段，形成了一种建立PBX细观结构模型的新方法。通过该方法建立的PBX细观结构模型的体积分数高达94%及以上，且同时还可以保证模型在边界处具有周期连续性。在此基础上，采用数值流形方法（Numerical Manifold Method，NMM）［12-13］对不同体积分数和不同颗粒级配条件下的PBX细观模型的单轴压缩等效模量进行了计算和分析。
2PBX细观结构计算模型生成算法
2.1PBX细观结构计算模型生成流程
PBX细观结构计算模型的生成过程如图1所示。首先由蒙特卡洛方法生成给定级配的随机圆［10］，然后以圆心作为Voronoi多边形的中心点生成Voronoi多边形，并将每个Voronoi多边形的顶点按逆时针
顺序排序，进而对Voronoi多边形进行修正以保持级配关系并消除细长失真颗粒。在此基础上，根据指定的颗粒体积分数对Voronoi多边形进行收缩，生成颗粒间的粘结层，并对颗粒进行切角和光滑处理等操作，得到PBX细观结构几何模型。最后，将细观几何模型形成ANSYS命令流，通过ANSYS进行网格划分，最终得到用于数值模拟的PBX细观结构计算模型。
2.2考虑级配的Voronoi多边形颗粒模型的生成
蒙特卡洛方法生成圆形颗粒时，需指定圆形颗粒的尺寸范围、级配个数和各级配颗粒尺寸等特征参数，采用转化后的二维Fuller公式［11］计算出各级颗粒在投放区域所占比例，进而计算出各级颗粒的个数。最后，按均匀分布的方式，确定圆形颗粒的位置，完成圆形颗粒的蒙特卡洛方法生成过程。
根据文献［8］所采用的颗粒尺度，将颗粒尺寸从小到大分为3个等级，即3个级配，投放区域的大小为1.5mm×1.5mm，表1给出了三种不同级配类型下的颗粒尺寸范围及各级配颗粒尺寸。图2所示为三种级配类型下各生成的一组随机圆形颗粒。级配类型的不同，得到的圆形颗粒所占的体积分数也不同，随着级配中大颗粒尺寸的增加，颗粒的体积分数随之增加。图2中，颗粒体积分数最高仅达到62%，远低于90%，这也是不直接采用这种模型作为PBX细观模型的主要原因。
为了提高颗粒体积分数，对于以上生成的圆形颗粒模型，将其中的圆心视为投放区域内的一组随机分布点，并以这些随机分布点为Voronoi多边形的中心，进一步生成相应的Voronoi多边形模型。为了使
得到的Voronoi多边形模型具有周期性边界条件，在投放区域内，每生成一个随机点，同时也在周边8个邻域分区内相对应的另外8个位置额外生成8个点。以图中的模型为例，生成的周期性点的分布情况如图3a所示，相应生成的Voronoi多边形模型如图3b所示。
2.3Voronoi多边形颗粒模型的修正
从图3b可以看到，对应于图2a中的随机圆形
颗a.periodically distributed points b.generated Voronoi
polygon model
图3周期性点及投放区域内的Voronoi多边形模型
Fig.3The periodically distributed points and the Voronoi polygon model in launch area
表1三种不同级配类型的给定参数
Table1Giving parameters for the three different gradation types
gradation type
Ⅰ
Ⅱ
硬质合金密封环Ⅲ
range of particle
diameter/mm
0.07-0.35
0.07-0.40
0.07-0.50
particle diameter of
each gradation/mm
0.07，0.18，0.35
0.07，0.16，0.40
0.07，0.18，
0.50
图2蒙特卡洛方法生成的圆形颗粒模型
Fig.2The circular particle models generated by the Monte Carlo
method
图1PBX细观结构计算模型生成过程
Fig.1The generation procedure of meso‑scale structure
computational model for PBX
粒，以3个大圆的圆心（A、B、C）为中心的Voronoi多边形颗粒，周围布满了细长的Voronoi多边形。根据Voronoi多边形的生成算法，介于大圆圆心和小圆圆心之间的Voronoi多边形的边界，是这两个圆
心连线的垂直平分线，所以这些细长的多边形是由大圆和小圆之间的尺寸突变引起的。为了消除这种细长颗粒，且同时使得模型中颗粒间保持指定的级配关系，需要对生成的Voronoi多边形进行修正。如图3b所示，筛选出相应面积变化较大的Voronoi多边形颗粒作为修正对象，比如圆心A所在的Voronoi多边形颗粒，以点A为圆心，用点A所对应的图2a中的圆截切与之相交的Voronoi多边形，删除圆内的线段，形成新的Voronoi多边形颗粒。按照这种修正方法，图3b所示的Voronoi多边形颗粒模型经修正后如图4所示。一般，修正后的颗粒模型在边界处仍然保持有周期连续性。
图5所示为修正后的Voronoi多边形颗粒模型与相应的蒙特卡洛随机圆形颗粒模型中的颗粒尺度统计规律对比，其中横轴为模型中颗粒按尺度由小到大的排列序号。图5中的横坐标显示，蒙特卡洛随机圆形颗粒有180个，修正后得到的Voronoi多边形颗粒有220个，颗粒数目明显增多，增多的颗粒填充了原随机圆形颗粒间的空隙。从图5中还可以看出，在修正后的Voronoi多边形颗粒模型中，粒径接近0.07mm的颗粒数量有所增多，且其粒径有缓慢增大的趋势，而中等粒径及大粒径的颗粒的尺寸却有所减小，整个模型中颗粒的级配关系仍较为明显，与原随机圆形颗粒的级配相对较一致。
2.4PBX细观结构计算模型的建立
修正后的Voronoi多边形颗粒细观模型的颗粒体积分数实际为100%。为了得到指定的体积分数，同时形成颗粒间的粘结剂层，需要对多边形颗粒进行收缩处理。个别不规则的颗粒在收缩后还可能产生尖
角，为了避免其对随后网格划分质量的影响，对小于给定角度阀值的角进行切削，并规定切削后形成的短边的边长不能小于多边形颗粒收缩的距离。图6
所示为模型经颗粒
图5修正后的Voronoi多边形颗粒模型与相应蒙特卡洛随机
圆形颗粒模型间的颗粒级配比较
Fig.5Particle gradation comparison between the modified
Voronoi polygon particle model and the corresponding Mon‑
te Carlo circular particle
model
图4修正后的Voronoi多边形颗粒模型
Fig.4The modified Voronoi polygon particle
model
图6一组不同颗粒体积分数的PBX细观结构模型
Fig.6A group of meso‑scale structure models for PBX with
different particle volume fractions
收缩、切角后，得到的一组具有不同体积分数的PBX 细观结构模型，其中颗粒间的空白区域为粘结剂层。
对于图6所示的细观结构模型，给出了模型的颗粒体积分数随粘结剂厚度的变化关系，如图7所示。从图7可以看出，随粘结剂厚度的增加，细观模型的体积分数按接近线性的规律减小。一般地，在生成粘结层时，模型中的个别细小颗粒经过收缩和切角后，可能会出现极其细长的形状，为了避免这些颗粒对网格划分和模型计算性能的影响，对其予以删除。当删除发生在模型边界附近时，这会导致模型的周期连续性变弱。此外，在收缩、切角后，模型中多边形颗粒的某些边，常会出现边长极小的情况，这也会影响模型的网格划分。为此，收缩、切角后，进一步在多边形颗粒的每条边上插入若干点，通过插入点对多边形进行样条曲线平滑处理。将颗粒平滑后的PBX细观结构模型形成ANSYS 命令流，通过ANSYS软件建模并划分网格，建立PBX的细观结构计算模型。图8所示为通过以上方法最终建立的具有不同级配类型的三个PBX细观结构计算模型示例。
3细观模型等效弹性模量的NMM模拟
数值流形方法（NMM）［12-13］通过采用双层覆盖技术，将有限元、非连续变形分析（DDA）和解析法统一在同一个理论框架下，能够方便的实现带裂纹的结构、多物体连续与非连续大变形与大位移，以及多物理场的耦合等复杂问题的求解。自20世纪90年代诞生以来，NMM在理论研究和不连续变形及裂纹扩展分析应用等方面得到了较快的发展［14］。
本研究采用NMM方法计算单轴压缩条件下PBX 细观模型的等效弹性模量，以验证本文所提出的PBX 细观结构计算模型建立方法的有效性，并初步探索颗粒体积分数及颗粒级配对PBX等效模量的影响。模拟中颗粒采用各向同性弹性本构，粘结剂采用22单元的Prony级数粘弹性本构［6，8］，颗粒与粘结剂之间的界面采用双线性接触损伤模型［15-16］。各组分的材料参数取值与文献［3］相同，具体见表2。
根据表1中的三种不同级配参数，建立三组不同级配类型的PBX细观结构计算模型。在每一组模型的生成过程中，通过改变Voronoi多边形颗粒收缩时的缩进量，分别得到十个不同颗粒体积分数的计算模型。其中不同级配类型的三个体积分数接近85%的PBX细观结构模型如图8所示。三组细观模型的总体尺寸均为1.5mm×1.5mm。模型底部设置一组固定点，约束底部Y方向的位移，而其沿X方向自由。模型上端设置一组位移加载点，沿Y方向向下加载，加载速度4mm·s-1。模型左端外侧设置一固定板，
通过固定板限制细观模型左端沿X方向上的位移，而细观模型左端沿Y方向自由。为保证细观模型右边界在X方向上的位移始终相等，在模型右边界设置一耦合运动约束板，约束板在X方向自由，在Y方向不允许发生位移，细观模型右边界与约束板进行X方向的粘结，允许细观模型右边界与约束板间发生Y方向的相对位移。
在NMM模拟中，在每一个加载步计算得到单元的应力应变值后，按照文献［3，17-18］中的方法计算当前时刻的等效弹性模量。根据该方法，基于平面应变假设，细观模型上的平均应力和平均应变可以表示为：-σ
i=
1
A∫Aσi d A，-εi=1A∫Aεi d A（1）式中，i表示应力/应变张量的方向，A表示二维PBX细表2PBX各组分的材料参数
Table2Material parameters of each component in PBX parameters
particle
binder
interface
elastic
modulus
/GPa
8.86
0.397
-
Poisson′s
ratio
0.322
0.495
-
density
/kg·m-3
1910
1280
-
tensile
strength
/MPa
10
6
1.66
friction
angle
/（°）
30
25
管道防爬刺15
cohesive
彩印业务strength
/MPa
30
6
1.66
图8PBX细观结构计算模型
Fig.8The computational models of PBX meso‑scale
structure
图7颗粒体积分数随粘结剂厚度的变化关系
Fig.7Relationship between the particle volume fraction
and the matrix thickness
搪瓷标牌观结构模型的面积。根据平面弹性应力‑应变关系，可以得到PBX 在平面应变假设下的等效弹性模量和泊松比为：
---E eff =-σy 2--σx 2-σy -εy --σx -εx ，----μeff =--εx -εy
（2）
再根据文献［19］中平面应变和平面应力之间的转换关系，得到PBX 的等效泊松比和等效弹性模量最终表达式如下：μeff =----μeff 1+----μeff
，E eff =---E eff (1-μ2
eff )
（3）
在单轴压缩模拟中，根据Y 方向的加载位移换算得到模型的即时Y 方向应变，结合式（3）中得到的当前等效弹性模量，就可以得到模型在Y 方向上的总体等效应力‑应变关系。
图9所示为NMM 模拟所得三组细观模型在Y 方向的等效应力‑应变关系曲线。在三种不同的颗粒级配条件下，当Y 方向上的应变值一定（即加载位移相同）时，随着颗粒体积分数的增大，Y 方向上的应力都逐渐增大。这是由于当体积分数较小时，颗粒间较厚的粘结剂层相对于颗粒起到了更多的变形作用，模型整体体现出粘结剂的“较软”性质，应力随应变的变化较缓慢；而当体积分数较大，如大于90%时，颗粒间的粘结层很薄，在压载作用下颗粒起到了主要的支撑作用，所以模型整体体现出颗粒的“较硬”性质，此时应力应变曲线陡峭。
图10给出了三种颗粒级配条件下细观模型的等效模量随体积分数的变化情况，其中等效模量是图9中应力‑应变曲线初始线弹性段的斜率。在单轴压缩加载过程中，随着加载位移的增大，颗粒和基体界面上会出现损伤，从而导致等效应力‑应变关系逐渐呈现非线性。从图10可以看到，当体积分数较小时，等效弹性模量的增速较为缓慢，此时粘结剂起到较大的变形作用；而当体积分数较大时，等效弹性模量的增速变大，这表明颗粒起到的支撑作用越趋明显。这一现象和文献［3］中的结果吻合较好。对于Ⅰ，Ⅱ，Ⅲ三组不同级配的模型，当体积分数低于30%时，三种颗粒级配条件下得到的等效模量基本相等，这是因为此时粘结剂对等效模量的确定起到了支配作用；与此相似，当体积分数高于85%时，三种颗粒级配条件下得到的等效模量也基本相等，但这时颗粒对等效模量的确定起到了支配作用；而当体积分数为30%~85%时，颗粒和粘结剂的支配作用都相对较弱，因此
颗粒的级配对等效模量产生了较明显的影响。
图10单轴压缩下等效模量随颗粒体积分数的变化关系Fig.10
Relationships between the equivalent modulus and
the particle volume fraction under uniaxial
compression
Ⅰ
Ⅱ
c.
type Ⅲ
图9不同颗粒体积分数下的应力‑应变曲线
Fig.9Stress‑strain curves under different particle volume
fractions
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议