郭支明，秦栋泽，贾宪振

(1.中北大学，太原 030051；2.西安近代化学研究所，西安 710065)

0 引言

PBX9501炸药(高聚物粘结炸药)在外载荷作用下内部微缺陷会形核长大，并形成微裂纹进行演化，形成内部损伤，使其力学性能劣化。炸药的损伤断裂性能实验主要有直接拉伸实验、巴西实验、三点弯曲实验等。巴西实验的试样制备简单，而且方便快捷[1]。由于PBX9501炸药属于粘弹性-准脆性材料，准确描述其力学响应及断裂特性比较困难。国内一些学者对此方向进行了研究，赵四海[2]等利用VISCO-SCRAM 模型，计算了巴西压裂实验中PBX9501炸药圆盘试样的力学响应，并和试验数据进行了对比。陈鹏万等[3]采用巴西实验研究了PBX9501炸药的断裂性能，并通过显微镜观察了炸药颗粒断裂、界面脱粘、粘结剂撕裂等破坏形式。庞海燕[4]等采用标准巴西实验、圆弧巴西实验和橡胶垫巴西实验测量得到PBX炸药的拉伸强度。在数值计算方面，崔云霄[5]等基于EFG方法和线性内聚力模型，模拟了准静态加载圆弧巴西实验中PBX9501炸药的断裂行为。戴开达等[6]采用扩展有限元方法，对不同种类巴西实验中PBX炸药的裂纹起裂、扩展和断裂行为进行了数值模拟。颜学坚[7]基于液压致裂原理，自主搭建了测试平台,发展了拉伸强度液压致裂测试方法，实现了PBX炸药的拉伸强度的准确测试。袁洪魏等[8]基于Boltzmann本构模型，建立了考虑静水压力影响的TATB基PBX准静态Boltzmann-P非线性弹性本构模型。李尚昆等[9]从材料的力学行为特性、实验方法、本构模型和强度理论四个方面，对高聚物粘结炸药(PBX)的力学性能特征进行了归纳和评述。蒙君煚等[10]用分子动力学方法，模拟了功能助剂对界面结合能的影响规律，试验结果与数值模拟结果吻合。原曾年等[11]采用CT技术和数值仿真的方式，研究了PBX材料样品内的缺陷和损伤行为，模拟结果与试验结果吻合较好。 性能测试方面，温茂萍等[12]采用圆弧压头与差动变压器引伸计相结合的巴西试验方法，测试了脆性炸药力学性能。陈科全等[13]以压装HMX基PBX炸药为对象,采用改进后的圆弧巴西试验,研究了不同预制缺陷条件下炸药的破坏过程,并与无缺陷炸药的试验结果进行对比。刘晨等[14]设计了一种带中心圆孔的平板试样进行实验，基于数字图像相关方法(DICM)，对实验进行应变场分析。

综上所述，炸药巴西压裂过程的数值计算方法采用有限元或EFG方法(有背景网格)等基于控制方程弱形式的离散格式，这些基于网格的方法在计算固体断裂等强间断问题时，受到网格连通性的限制。本文使用基于控制方程强形式的离散格式——CRKSPH(守恒型再生核粒子方法)，对炸药巴西压裂过程进行模拟。这种纯拉格朗日方法不再依赖于网格，在求解断裂和裂纹扩展问题上有独特的优势。

1 基本理论

CRKSPH方法是SPH(Smoothed particle hydrodynamics，光滑粒子流体动力学)方法的一种改进。SPH是一种拉格朗日方法，它基于粒子近似对控制介质运动的积分或偏微分方程进行离散，连续场变量的值通过将离散粒子的性质(如质量，动量和能量)与内插核(通常用W表示)进行卷积来表示，是一种完全拉格朗日无网格数值方法，物质界面清晰，材料的变形不依赖于网格，而通过粒子的运动来描述。因此，SPH既具有拉氏计算中描述物质界面准确的优势，又兼备无网格方法的长处，适宜计算带有流体大变形及运动边界的各类问题，如材料损伤断裂、爆炸爆轰、冲击侵彻等[15-22]。SPH方法最早由LUCY[23]和GINGOLD[24]等于1977年提出，并应用于天体物理学，SWEGLE 等[25]将SPH方法用于模拟爆炸问题。

尽管SPH方法已成功应用于许多领域，但传统的SPH方法仍存在严重的缺陷。其中，最严重的是SPH缺乏零阶一致性，即所谓的“E0级误差”[26-27]。SPH的另一个常见问题是精确捕捉冲击波阵面所需的人工粘性公式。 MONAGHAN和GINGOLD[28]的标准粘性计算公式因阻尼过大，使结果严重失真。目前为止，人工粘性计算仍是一个尚未解决的问题。

解决零阶精度问题一个有意义的尝试是再生核(RK)方法[18,29]， RK方法以附加的自由度对SPH插值内核进行了改进，可以精确地再生任意阶的多项式。这完全消除了SPH的零阶误差，但这种方法使粒子之间的内核不对称。本文使用CRKSPH(守恒型再生核粒子方法)[29]，其使用线性再生核，基于RK插值重新推导了动量守恒方程。这种方法可以精确地保持线性动量守恒，以达到一定精度，可有效消除零阶插值误差和过度使用人工粘度。

CRKSPH保留了传统SPH方法的许多优点，同时对SPH的缺点进行了改进，特别是人工粘性过高和零阶精度误差。为此，CRKSPH方法对传统SPH方法进行三个变化：线性再生核函数、守恒型的控制方程、精确的人工粘性算子。

1.1 线性再生核函数

为解决SPH无法重构所需阶数的多项式，在传统SPH插值核函数中添加一些修正项，以精确再生常量、线性或更高阶数的项。这种插值方法，称为再生核方法(RPKM)，再生核公式可以推广到任意阶。 RK插值公式和变量场的梯度表示为

(1)

(2)

1.2 守恒型控制方程

关于再生核插值和梯度，在建立守恒型控制方程时，CRKSPH方法利用了MLSPH构造中的一些派生方法。

假设插值使用通用的核函数ψ(x)，其中守恒密度U和通量F定义为

(3)

通过插值函数和体积积分：

(4)

可得到动力动量和能量方程：

(5)

其中的力对是反对称的。因此，这样推导出来的动量方程，再加上再生核形式，可以保证精确的线性动量守恒。

1.3 人工粘性算子

这里使用经典的van Leer限制器[30-31]作为新算子的基础，这个算子是基于网格方法的计算流体力学理论。修正后的算子如下：

(6)

(7)

(8)

这里使用普通的RK梯度算子来计算该速度梯度：

(9)

等式(6)的第一部分，就是普通的van Leer限制器；当ηij<ηcrit时激活，当点靠近时将限制器强制为零。

2 损伤演化模型

PBX9501材料在拉伸过程中，随着变形的增加,损伤演化可表示为应变的函数。对于PBX9501材料，可建立如下损伤演化方程[32]，损伤参数D(ε)演化根据公式：

(10)

式中B和n为材料常数；εc为材料的断裂应变；ε0为损伤开始发展时的临界应变。

对式(10)积分，有

(11)

式中D0为材料的初始损伤。

PBX9501材料的初始损伤D0可通过材料的实际密度与理论密度的偏差来初步估计。根据文献[30-31]实验结果，可确定损伤演化方程中的材料常数，如表1所示。

表1 PBX9501材料参数

从而可建立材料的损伤本构关系：

σ=(1-D)Eε

(12)

式中σ为材料应力张量；E为弹性模量矩阵；ε为应变张量。

3 数值模型和仿真结果

计算时，将三维巴西实验问题简化为平面应力问题，建立的CRKSPH模型如图1左图所示，圆盘试样直径为20 mm，共划分14 400个粒子。为了模拟真实的实验加载，圆盘下部放置在刚性面上，上部施加位移载荷0.05 mm/min。图1右图为标准巴西压裂试验基本破坏形式示意图。

图1 CRKSPH计算粒子模型(左)和标准巴西压裂试验基本破坏形式 (右)

3.1 仿真结果

以上文所述加载工况(加载速率0.05 mm/min)和粒子模型进行巴西圆盘仿真实验，位移演化过程如表2所示，损伤及裂纹演化情况如表3所示。

表2 PBX9501试件不同时刻位移云图

在整个仿真试验过程中，试件沿载荷方向处于压缩状态，沿水平方向处于拉伸状态；在载荷作用点(接触点)附件一定区域内会发生剪切破坏，最终在试件内部沿着载荷作用方向会形成一条贯穿裂纹，将试件一分为二，计算结果可充分地展现这一过程。

从表2和表3中的结果可知：

(1)t=0.36 s圆盘中心位置微缺陷发展为局部损伤，在随后的加载过程中逐步形成裂纹，并沿着纵向向上和向下扩展；

(2)在圆盘的加载点两侧，由于局部剪切作用产生局部破坏；

(3)在试件的底部附件区域由于剪切和拉伸综合作用也发展出了多处不同程度的损伤裂纹，损伤情况基本以圆盘纵向对称轴呈对称分布。

将表2中位移云图和图1右图中基本破坏形式示意图做对比，可看出仿真结果和巴西压裂基本破坏模式保持一致。图2所示为0.5 s时微观缺陷分布、损伤参数D和文献中试验结果[6]。

图2 0.5 s时微观缺陷分布(左)、损伤参数D和文献中试验结果(中和右)

从图2断裂情况和文献中试验结果[6]进行对比可知，试件损伤、裂纹的产生和扩展方向均和试验基本吻合，而且基于CRKSPH方法仿真，除了可追踪明显的裂纹外，还可对产生损伤的位置进行预测和追踪，并且可有效处理多裂纹的产生和合并。

从图3所示试件顶部裂纹分布和试验数据对比情况可看出，仿真结果中顶部3条主裂纹分布和试验结果保持一致[4]。图4所示为基于本构关系(8)计算得到的试件中心位置应力-应变曲线。

图3 试件顶部裂纹分布情况对比 (左为试验结果和右为仿真结果)

图4 计算得到的应力位移曲线

从图4中可看出，应变在0～0.002 6这个阶段，应力基本呈线性增大，材料处于弹性变形阶段；当应变大于0.002 6时，应力达到最大值5.87 MPa，此时试件中心的拉应力超过抗拉强度，从而产生初始裂纹；裂纹继续扩展，导致试件中应力减小，当应力减小到2.857 MPa时，形成贯穿裂纹。

3.2 不同加载速率对试件断裂形式的影响

由于炸药在运输和使用过程中会遇到各种载荷情况，有必要考虑加载速率对断裂形式的影响。本节基于上文的计算方法和模型，考虑不同加载速率下试件的断裂形态和规律。

工况和计算结果如表4所示。这里主要考虑的加载速率为0.03、0.05、0.07 mm/min。

表4 三种加载速率下相同中心位置位移对应的断裂形态

从计算结果可看出：

(1)加载速率不同，试件中裂纹起始位置不同(如表4中第一行各图中的黄色圆圈标记处)。加载速率为0.03 mm/min时，裂纹起始位置位于试件中心点下方2.4 mm处，加载速率为0.05 mm/min时，裂纹起始位置位于中心点，加载速率为0.07 mm/min时，裂纹起始位置有2个，一个位于中心点上方2.5 mm处，另一个位于中心点下方2.5 mm处。

(2)试件中心点位移相同的情况下，各加载速率对应的断裂形态不同。加载速率越大，压裂所产生的裂纹数量越多，而且在试件外表面，裂纹数量从底部往上沿外边缘快速增加。

4 结论

(1)采用 CRKSPH方法，对PBX9501炸药标准巴西圆盘压裂过程计算结果中，试件破坏方式和裂纹走向与文献中实验结果基本保持一致。

(2)不同加载速率对试件中裂纹的产生和扩展都有一定的影响，随着加载速率的增加，试件中裂纹数量迅速增加，碎裂程度增大。

(3)计算结果表明，CRKSPH和材料的损伤本构相结合，能够有效求解固体炸药损伤、断裂和裂纹扩展问题，可有效处理弹脆性材料中多裂纹产生、合并等物理现象。计算方法和结果对炸药材料工程设计有一定的指导意义。

本文主要以探索新算法在炸药材料破坏仿真方面的应用，试验数据和图片均取自相关文献，下一步研究中，将会结合试验数据对计算模型进行进一步修正。