炸药装药损伤行为数值模拟研究进展

张萌昭，屈可朋，沈 飞，吴翰林，周 涛

(西安近代化学研究所， 西安 710065)

炸药是武器战斗部的核心能源，是为武器提供杀伤性、破坏性、动力性的关键材料，其安全性也一直备受关注[1]。然而，在生产、运输、发射过程中，由于加载冲击[2]、应力状态、环境温度[3]、生产工艺等各方面复杂因素的影响，炸药会出现各种形式的损伤，其力学性能、感度、破坏性等特性也会因此产生变化，甚至导致炸药提前引爆，产生不可预估的后果[4]。为了保障战斗部能够安全稳定地实现打击效果，必须对炸药的损伤行为进行系统而深入的研究。

目前国内外对炸药损伤行为的研究主要依靠两种方法，一种主要依赖于实验获得相关数据，如目前常用的Hopkinson压杆实验、轻气炮实验、巴西实验[5]、落锤试验[6]等方式。实验方式能够真实、准确地得到不同炸药在不同条件下的应力应变关系及损伤情况，但是由于炸药的不稳定性，这种研究方法具有一定的危险性。此外试验成本、实验条件也往往会限制试验过程的开展，得到的实验结果一般不具有广泛性和预测性，对于炸药细微观损伤的实时数据采集也很难实现。另一种则是数值模拟方法，国内外许多学者对炸药的损伤特性进行了模拟研究，基于LS-DYNA等软件，可以同时得到空间和时间尺度上的仿真数据，进而能够观测到材料细微观尺度的内部变化过程。但是这种方法的准确性和适用范围常受到材料的本构模型和模拟计算方法的影响。为了对目前炸药损伤问题中常见的本构模型及算法进行详细分析，列举了国内外炸药损伤本构模型及炸药装药力学行为的数值计算模型，对比分析其优缺点及适用范围，以期为研究人员选择本构模型和仿真算法提供参考。

1 炸药含损伤本构模型研究现状

材料的本构关系能够通过数学表达式来描述材料在外界因素作用下的力学行为及变形规律。通过给定相关参数，通过本构模型即可得到细致、准确的力学状态。下面介绍几种损伤本构模型。

1.1 粘弹性统计裂纹模型

(1)

其中:K为含裂纹扩展的等效力强度因子;K0为材料的断裂韧性；vmax为最大扩展速度；m为速度系数。

(2)

其中：c为裂纹体的平均半径；G为弹性元的剪切模量；τ为松弛时间；a为初始缺陷尺寸；n和N分别代表第n个弹性元、体元数量。

Vioco-SCRAM模型与SCRAM模型相比结构简单、参数少，因此应用更加广泛。除了可以计算裂纹的生长过程外，还可用于模拟炸药的非冲击点火过程，如孙宝平就基于该模型计算了侵彻过程中装药的整体升温和点火机制[11]。但是该模型考虑的因素较多，各参量的确定虽然与SCRAM模型相比较为简单，但仍然是一项繁琐的工作。此外，该模型得到的计算结果与实验结果相比，在某些情况下仍然有较大差距，而且对于炸药细微观的损伤现象表述仍不够成熟和完善，需要进一步发展和改进。

针对Visco-SCRAM模型的不足，一些研究人员改进了模型参数的选取方式及参数取值范围，提高计算精度的同时简化了参数设置。如张延耿等[12]，基于币形微裂纹面的对称特性，利用半球体而非球体来描述裂纹的取向，将传统Visco-SCRAM中微裂纹方向矢量的取值范围缩小为原来的1/2。为了适应本构参数的变化，同时对初始缺陷尺寸表达式进行修正，得到了新的偏量本构关系率的表达形式，以不同的初始缺陷尺寸参数对应不同的应力状态。此外，考虑到拉伸条件下线性的体量本构关系不能够描述体积应变增量，因此结合文献[10]中给出的微裂纹扩展速度方程对体量本构关系进行改进，最终得到改进后的本构关系表达式。对该模型利用LS-DYNA程序进行仿真，通过对比PBX炸药的平板撞击试验数值模拟的计算结果，发现改进后的模型能够更加合理、细致地反应炸药在低速冲击作用下的细微观损伤破坏机制，但是该模型对一些参数有较强依赖性，数值模拟结果仍存在一定差异。

1.2 含裂纹扩展损伤的粘塑性本构方程

现在国内外的微裂纹统计损伤模型，研究的重点较多基于炸药的粘弹性行为，但一些力学性能试验表明，炸药在外力的作用下也会表现出明显的塑性特性。基于这种现象，成丽蓉等[13]建立了能够反映PBX材料的粘塑性特性的含损伤本构方程。该模型同样假设初始状态的微裂纹各向同性分布，并将一个阻尼器和塑性单元并联以后与弹性体串联以构建弹粘塑性模型，如图1。

图1 弹粘塑性模型示意图

得到的偏应力张量表达式如下：

(3)

为了验证该模型的准确性，成丽蓉等人对PBX炸药试件进行了静动态单轴压缩力学实验及断裂性能试验，得到的试验数据与仿真结果进行比较，发现虽然该模型能够描述PBX炸药的粘塑性及脆性损伤力学的响应过程，但是模拟曲线在应力强度超过最大值以后裂纹扩展速率迅速增大，模拟曲线与实验数据并未完全吻合，尤其在应变率较高的情况下拟合效果有较大偏差。

1.3 含各向异性损伤的粘弹性统计微裂纹本构模型

除了上文提出的基于初始材料各向同性假设的本构模型外，也有一些模型考虑了微裂纹法向造成的影响，提出了含各向异性损伤的粘弹性统计微裂纹本构模型(Aniso-Visco-SCRAM)。如图2。

图2 Aniso-Visco-SCRAM模型示意图

张延耿等[14]改变了原本Visco-SCRAM模型中利用积分形式计算微裂纹扩展产生的应变的方法，利用叠加法将不同法向上的微裂纹应变分量求和，得到宏观上的微裂纹扩展产生的附加应变。与传统Visco-SCRAM模型得到的裂纹在各个角度扩展情况相同相比，Aniso-Visco-SCRAM更好地描述了PBX炸药在低速冲击载荷下的损伤机理和变形特征。但是该模型所需的部分参数无法通过试验确定，而且这些参数的取值对模型的准确性有较大的影响，所以需要更进一步的定量研究以更好地描述损伤破坏过程。Marwen Chatt等[15]也对含各向异性损伤的本构模型进行了研究。他所提出的模型是一种新的粘弹性微平面模型，使用Maxwell链模型和微平面方法并与含各向异性损伤结合来模拟准脆性PBX材料的行为。该模型基于Carol等人提出的V-D-T分解方法，在微平面上将应变张量沿微平面方向投影，之后采用内部微平面法，在每个方向都定义一个Maxwell链模型，该模型由10个粘弹性体和1个弹性体并联而成。V-D-T分解方法如图3所示。

图3 V-D-T分解示意图

定义微平面内部的行为需要确定每个微平面方向和每个元素的弹性应变和弹性应力。并且该模型认为损伤的演化仅与弹性元件的相互作用力有关而与粘弹性元件无关。与实验数据相比，当考虑包络曲线、横向应变和松弛步阶时，数值预测与实验数据拟合较好。但在描述压缩的模拟卸载过程时，所提出的模型仍然缺乏准确性，没有考虑材料的塑性和不可逆性的损伤。同时该模型所需的参数较多，需要高性能的计算机提供支持。

1.4 内聚裂纹模型

由于炸药在某些情况下与混凝土等有相似的力学特性，也有一些由基于混凝土材料的本构模型发展而来的用于炸药的本构模型，例如内聚裂纹模型。内聚裂纹模型可以计算材料任意位置的裂纹的产生或扩展，该模型所需的材料参数少，物理意义明确。该模型的特点在于能够避免处理裂纹尖端的奇异的特点，通过单元函数的扩展而直接作用于网格内部。内聚裂纹模型的基本方程如下：

(4)

崔云霄[16]等人基于内聚裂纹模型对PBX模拟材料进行了动态巴西试验的数值模拟，利用LS-DYNA软件得到变形场及应变测试结果，得到拉伸应变的峰值误差小于5%，对于裂纹张开角度的模拟误差也能控制在1°之内，因此内聚裂纹模型能够较好地模拟PBX-M材料在动态巴西实验过程中的裂纹扩展过程。吴艳青[17]等人也采用了内聚力模型来研究PBX的损伤，该模型的特点在于分别考虑了粘弹性和脆性损伤的来源，认为炸药的粘弹性主要来自粘结剂，因此基于HMX的弹塑性模型将粘结剂的功能通过内聚力模型表现，进而直接置于HMX晶体的外界面。

1.5 Z-W-T模型

我国国内使用较多的本构模型由朱兆祥、王礼立和唐志平等人，基于Green-Rivlin本构理论提出的Z-W-T模型。该模型将一个非线性弹簧和两个Maxwell体并联进而描述高聚物的非线性粘弹性力学行为。其方程形式为：

(5)

其中:E0、α、β、E1、E2、θ1、θ2均为材料常数。但是这种模型难以用于描述PBX炸药的损伤情况。因此孙文旭等[18]将损伤行为引入Z-W-T模型，建立了含损伤的非线性粘弹性本构模型。他们对一种新型的抗过载浇铸PBX炸药进行准静态力学实验和分离式霍普金森压杆实验，通过遗传算法拟合，对高应变率下PBX炸药在到达破坏前的动态力学行为进行了较好的描述。屈可朋等[19]采用含颗粒相结构应变率效应的Z-W-T模型对RDX基含铝炸药的冲击特性进行研究，与通过分离式霍普金森压杆实验得到的实验数据拟合较好，表明该模型可以较好地描述RDX基含铝炸药在高应变率下的损伤行为。

2 数值模拟方法研究现状

随着计算机技术的发展，炸药损伤的数值模拟方法得到了极大的发展，一定程度上弥补了传统试验方法的限制和不足，得到了广泛应用。数值模拟方法能够全方位且直观精确地反应炸药在不同尺度上的力学响应特性。下面将对一些数值模拟方法进行介绍。

2.1 有限元方法

有限元方法是一种常用的研究PBX力学行为的方法，这种方法基于连续介质力学，能够在宏观尺度上研究材料的力学性能，其主要思想在于通过变分方法求解出偏微分方程的近似解，求解时将物体离散化，并对其中每个微小单元进行求解，并通过单元与单元之间的关系进行约束，具有计算速度快，而且相关的计算软件较为成熟的优点。石啸海等[20]就利用有限元方法，针对战斗部在侵彻过程中产生的装药损伤问题进行研究，分别对两组模型进行了仿真模拟，计算了PBX的装药损伤演化过程，并且分析了裂纹宽度、过载、轴向应力和损伤程度等因素，比较了有机玻璃、聚四氟乙烯、酚醛树脂、尼龙四种材料的缓冲性能，验证了内聚力本构模型的准确性。魏强等[21]则基于 Karagozian & Case模型，利用有限元方法对PBX炸药在撞击作用下的损伤情况，引入了特征长度以保证计算精度。但是有限元方法也存在不足，在有限元方法中，网格的边界与模型的几何边界一致，因而在变形过程中会出现网格的畸变、缠绕等影响计算精度的问题，同时对于裂纹的生长过程的研究也有一定限制，很难对结构的细微观变化进行准确描述。

2.2 扩展有限元方法

为了弥补有限元方法的不足，发展出了扩展有限元方法(XFEM)。扩展有限元方法由Belytschko和Black[22]首先提出，该方法改变了传统有限元方法的划分方法，并不对结构内部的裂纹和不连续面等问题进行网格划分，而是将增强形函数加至传统有限元的位移模式中来代表计算区域内的间断，图4表示裂纹平板有限元模型，其中浅黑色部分为裂纹，图中方框节点表示被完全切断的集合，而圆圈节点则表示围绕裂尖的节点。当网格密度足够时，可以更加方便地模拟材料中的裂纹、空洞等情况[23]，具有计算精度高、无需重新划分网格等优点。

图4 裂纹平板有限元模型示意图

戴开达等[24]用扩展有限元方法对PBX炸药在多种巴西实验中的变形破坏过程进行了研究，当单元被裂纹完全穿过时使用的增强型函数为

ψJ(x)=NJ(x)H(f(x))

(6)

其中NJ(x)为通用增强型函数，H(x)在(x-x*)·n≥0取值为1，否则为-1。f(x)的计算式如下：

(7)

当节点处于裂尖周围时，则增强形函数为

(8)

基于上述2个增强形函数，可用下式来表示含裂纹二维板的位移场：

(9)

其中:S为所有节点的集合；Sh为被裂纹完全切断的单元的节点集合；Sc为围绕裂尖的单元的节点；aJ、bK为节点的附加自由度。

利用上述公式，基于扩展有限元法模拟了巴西实验中裂纹起裂、扩展和断裂等行为，并且得到了在不同加载形式下圆盘式样的不同应力应变分布。

Huang X C等[25]则基于扩展有限元模型和内聚模型，对含孔洞平板的力学性能进行研究，发现对该平板施加整体压力时，孔洞周围的局部拉伸力会导致局部裂纹产生。并且通过数值模拟得到了裂纹在整个试验过程中的整体走势、拐点，包括起裂时刻、裂纹初期展速等数值，证明了有限元方法和内聚模型的准确性。但是扩展有限元方法在解决三微裂纹扩展和群裂纹等方面仍然存在缺陷。

2.3 物质点法

上面介绍的方法计算过程依赖于网格的划分，存在一定的局限性，因此近些年来无网格法被提出并且逐渐发展。物质点法(MPM)就是常用的无网格法之一。MPM与网格法本质上的区别是物质点法是一种基于粒子的计算方法[26]，MPM首先将连续体通过某种方法离散成一个个具有集中质量的物质点，在计算过程中物质点的质量保持不变，之后通过建立背景网格节点上的有限元函数实现质点和背景的固连关系。他们之间的映射关系表示为：

(10)

其中：ψp为物质点的变量；ψi为网格节点变量；nu为单元节点数；xp为物质点坐标。在此映射关系下，动量方程的虚功形式可以表示为离散形式：

(11)

(12)

通过上式和材料本身的本构关系方程即可计算得到物质点的应力，并且物质点法的背景网格随时间不断更新，避免了有限元法的缺陷，同时在计算过程中MPM也满足了质量守恒定律和动量守恒定律[27]。唐红[28]等人基于物质点法，对奥克托今基高聚物粘结炸药的力学性质进行模拟，并且建立了炸药颗粒随机分布的力学模型，对PBX炸药压制过程中的二维数值模拟进行模拟，得到了炸药在炸药过程中颗粒和缝隙的变化过程，同时还检测了炸药体系温度的变化情况。MPM与其他方法相比，有着较高的效率、精度和稳定性。但是由于质点积分存在误差，对于小变形问题，MPM的精度要低于拉格朗日有限元方法，因此该方法更加适用于大变形，

2.4 无网格伽辽金方法

无网格伽辽金方法(EFG)是另外一种使用较为广泛的无网格方法。这种方法是由Belytschko T[29]提出，对形函数的导数和边界条件的施加方法进行改进，提高了计算的收敛速度和精度，弥补了传统有限元方法中容易出现的体积锁死现象的问题。其具体方法是通过最小二乘法构建基于一定已知点构造出区域内的近似函数，即：

(13)

其中：u(X)为计算域内函数；p(X)为基函数；a(X)为J(X)的极小值，与X相关，其中J(X)为：

(14)

其中权函数wi(X)有一个影响范围，其值非负，当超出该范围时其值取0。

崔云霄[30]等人基于EFG对PBX炸药的圆弧巴西实验进行了数值模拟，得到了x方向的应变峰值约为6%，与实验结果较一致，较为准确地得到了准静态压缩下裂纹的产生即扩展过程。

2.5 离散元方法

离散元方法(DEM)是一种较为成熟的数值模拟方法。这种方法是由Cundall[31]最早提出，基于牛顿第二定律，将块体单元设定为理想刚体，并且基于块体单元的几何形状和变形特征得到不同的形状、性质离散元，用显式差分的计算方法，先得到各个元素的运动方程，之后通过进一步计算得到物质的整体运动形态，适用于非连续截止问题即非连续体破坏问题。在炸药的冲击过程中，对细观模拟能够取得不错的计算结果和精度。金博[32]在研究炸药的反应速率与裂纹的扩展和生长之间关系时，采用DEM方法进行数据模拟，得到了不同发射速度下的裂纹分布、成长状况，进而利用裂纹成核理论描述反应过程，进一步得到了较为完整的反应速率模型。但是这种方法的精确度对主要参数的依赖性较高，计算效率较低。

2.6 近场动力学方法

此外近些年也有学者基于近场动力学来研究炸药损伤问题。近场动力学方法(PD)最早由美国的Silling[33]提出，这种方法和MPM相似，首先将固体离散为一系列带质量的物质点，但是PD考虑了物质点和其相邻的一定范围内物质点之间的相互作用力，如图5，这种相互作用力包含材料的物性信息，从而避免了传统的应力-应变形式，进而不需要连续性假设。

图5 物质点的运动及相互之间作用力示意图

但是以键为基础的动力学过于简化，在使用过程中受到限制，故PD理论又进一步发展成为以状态为基础的近场动力学理论Silling[34]。李潘[35]基于状态PD理论对裂纹在PBX炸药裂纹扩展中的应用进行了研究，使用了显式的计算方法，通过巴西圆盘实验、带长方孔平板拉伸实验、含裂纹板动态断裂试验等多种试验研究PBX炸药的力学行为，建立了近场动力学蠕变模拟方法。秦洪远等[36]则通过非局部键形近场动力学建模，对含初始裂纹的巴西圆盘劈裂过程中裂纹的扩展及破坏形式进行研究，形成了较为完整的近场动力学数值体系，得到了不同初始裂纹长度、不同初始中心裂纹倾角下的圆盘裂纹生长情况。PD方法不需要借助外部准则，避免了裂纹尖端奇异性问题，同时提高了计算效率。

对以上几种数值模拟方法进行总结，具体情况如表1所示。

表1 几种数值模拟方法

3 结论

针对炸药的损伤，列举了一些较为常用的含损伤本构模型，介绍了一些数值模拟方法。由于炸药较为复杂的力学特征，很难得到一个广泛适用的本构模型，虽然所介绍的本构模型都可以用来描述炸药的损伤行为，但是其侧重点和适用范围及验证方法都不尽相同。含损伤的本构模型在细微观层面的发展还有欠缺。

一般受限于实验条件，对炸药的损伤进行全面系统的受力情况、力学性能分析是比较困难的，因此数值模拟方法在实际工程中相当重要。但是现有文献介绍的数值模拟方法在多尺度分析方面有很大不足，而且对于数值模拟的理论分析较为欠缺。

为了更好地研究炸药装药损伤，需要在细微观层面进行更加深入的探讨和研究，构建更加合理、准确的含损伤细微观本构模型，提升数值模拟精度和计算效率。