烟幕初始云团爆炸分散模型建立及计算方法

陈 浩，高欣宝，李天鹏，张 倩，陈玉丹，杨 洋

(1.陆军工程大学石家庄校区， 石家庄 050003； 2.陆军北京军代局驻763厂军代室， 太原 030000)

为有效对抗制导武器的精确打击，提高部队战场生存能力，发展与之对抗的干扰理论与技术，受到世界各国军队的重视[1]。烟幕无源干扰凭借成本低、效费比高的优点已经成为干扰领域发展的重要研究方向。其中，分析烟幕初始云团的运动规律，计算烟幕初始云团的参数，是烟幕作战效能评估研究的热点之一[2-3]。爆炸分散型烟幕经中心扩爆管爆轰驱动，向四周迅速膨胀，当壳体破裂瞬间，形成了由爆轰产物与干扰剂小颗粒构成的气/固混合物云团，将其称为烟幕初始云团。从形成到烟幕粒子速度衰减为零的过程，称为烟幕初始云团的爆炸分散过程。将烟幕初始云团近似看成不断扩张的球体，其最大半径用于表征烟幕遮蔽面积，作为烟幕评估中的一项重要指标[4-6]。

为了探索烟幕初始云团运动规律，提高爆炸分散的理论模型精度和适用性，近年来，学者们做了大量工作。Chen等[7]研究并建立了真空环境中烟幕的膨胀模型；Zhu等[8]分析了赤磷烟幕初始云团膨胀过程，指出空气阻力对烟幕膨胀影响较大，利用常数变异法建立了基于膨胀力和空气阻力作用下的烟幕初始云团膨胀模型；Xu等[9]考虑忽略等熵过程中气动阻力作用的情况，建立烟幕初始云团膨胀模型，采用四阶龙格—库塔法对模型进行了计算。

上述研究针对烟幕初始云团爆炸分散形成过程建模与计算的不同方面进行了创新与改进，但仍存在以下问题：(1)根据数值计算结果分析，等熵膨胀阶段的气动阻力对模型精度有一定影响，建立理论模型时不能忽略；(2)理论模型采用四阶龙格—库塔法求解数值解的效率低，收敛性和平滑性存在不足。针对以上问题，为进一步提高模型精度及计算效率，本文对烟幕初始云团爆炸分散的理论模型和模型计算方法作进一步改进，将爆炸分散过程分为等熵膨胀和自由扩散两个阶段，考虑等熵膨胀阶段的气动阻力，推导并建立了烟幕初始云团爆炸分散理论模型，利用欧拉法求解模型的数值解，并与四阶龙格—库塔法进行对比分析。通过开展野外测试试验，对改进后的模型和所提算法进行了分析与验证。

1 理论部分

发烟装置起爆后，产生爆轰波，并迅速形成大量高温高压气体(中心黑索金药柱爆炸形成)向四周膨胀做功，驱动周围的碳基干扰剂向外挤压，壳体内部压力逐渐增大，当超过壳体承受极限后，发烟装置发生爆炸，形成一个近似球体且不断膨胀的高温高压云团。壳体膨胀到破裂过程时间极短，为使研究问题简化，假设爆轰前后爆轰产物和云团体积无变化，粒子速度为零，即认为爆轰瞬时完成。

爆炸分散过程实际为爆轰产物在空气介质中流动并驱动碳基干扰剂向外膨胀的过程，烟幕初始云团是爆轰产物与碳基干扰剂形成的混合物。根据等熵膨胀基本理论，同时，考虑便于计算爆轰产物压力，假设烟幕初始云团在内部压力衰减为零之前的过程视为爆轰波在空气中的绝热等熵膨胀过程，即爆炸分散过程的时间分界点为爆轰压力衰减为零时，分为等熵膨胀和自由膨胀两个阶段。

1.1 等熵膨胀阶段模型建立

取云团边缘一个小微元体，其质量为dm，受爆轰产物压力作用面积为dS(微元体体积远小于云团体积)，受到的作用力主要是爆轰产物压力和气动阻力，则有：

(1)

式(1)中：a为微元体膨胀过程的加速度，可表达为a=d2R/dt2；P为爆轰产物压强；k′为气动阻力系数；ρ′为标准条件下大气密度；v为微元体的速度。

微元体的状态方程可表示为：

p=A(S′)·ρk

(2)

式(2)中，A(S′)为常数。

为求得云团运动参数的数值解，在等熵膨胀假设前提下，得到以上微元体压力表达式，则可进一步得到微元体爆轰发生后任意时刻的状态与初始状态之间的关系如下：

(3)

式(3)中：P0为爆炸初始时刻形成的高温高压云团的压强；ρ0为爆炸初始时刻高温高压云团的密度；ρ为爆炸后等熵膨胀阶段任意时刻高温高压云团的密度；k为等熵指数。

假设爆炸后形成的云团为不断膨胀的球体，则云团密度可用体积和质量代替，并代入式(3)中可得：

(4)

式(4)中：m为烟幕初始云团质量；R为烟幕初始云团的膨胀半径。

将式(4)代入式(1)整理得到如下表达式：

(5)

将等熵膨胀阶段模型整理，得到如下微分方程组：

(6)

1.2 自由膨胀阶段模型建立

在等熵膨胀阶段，爆轰作用力在云团向外膨胀过程中不断衰减。当压力等于标准大气压时，进入自由扩散阶段。这一阶段，微元体主要受气动阻力作用，则：

(7)

由式(7)知，自由膨胀阶段，微元体加速度为负值，速度随时间不断衰减，直到速度衰减到零，此阶段结束。经整理化简得到烟幕初始云团爆炸分散自由膨胀阶段的微分方程组为：

(8)

2 实验

根据特种弹药战斗部设计原理，建立了一种外形为圆柱体的发烟装置装药模型，如图1所示。壳体材料为合金钢，厚度为2.50 mm，高为127.50 mm。中心装药为黑索金药柱，装填密度为1.70 g/cm3，装药半径为9.00 mm，雷管位于中心装药上端，用于起爆黑索金药柱。四周为碳基干扰剂，装填密度为2.30 g/cm3，装药半径为56.50 mm。

为有效获取烟幕初始云团爆炸分散过程中的图像数据，设计了发烟装置野外烟幕测试试验，图2为试验系统的示意图。试验主要利用高速摄像机对烟幕爆炸分散过程进行记录，并通过测距仪测量吊架两侧距离，基本参数值见表1。

图2 试验系统示意图

型号帧率/(帧/s)分辨率VEO4001 0001 280×800

利用获取的图像数据，采用图像边缘检测法，获取吊架两侧的像素点坐标和云团边界两侧的像素点坐标，通过比例关系求解云团径向膨胀的最大半径，其表达式为：

(9)

式(9)中：L为吊架两侧的距离；Δx1为图像中吊架两侧像素点距离；Δx2为图像中云团两侧像素点的最大距离。

3 结果与讨论

3.1 理论模型计算与分析

3.1.1基本原理

式(6)和式(8)为二阶非线性微分方程组，可采用数值积分法求得数值解。常用的数值积分法包括欧拉法和四阶龙格—库塔法等。一般针对低阶微分方程，为确保数值解的光滑性，采用欧拉法进行计算，计算公式如下：

(10)

为进一步计算模型的数值解，首先需确定初值。根据中心装药的体积和装药密度等参数，根据Liu等[10]提供的方法，通过计算不同氮当量计算爆压的公式为：

P′=1.092(ρD∑N)2-5.74

(11)

式(11)中：P′为炸药爆轰压力；ρD为炸药装药密度；∑N为炸药氮当量。

通过式(10)计算得到瞬时爆轰时刻，爆轰压力P′=27.76 GPa。圆柱壳体破裂的强度为1.20 GPa，因此，可确定求解理论模型爆炸瞬时的初值为：P0=1.20 GPa，R0=56.50 mm。

3.1.2数值计算与分析

根据初始参数，编写了欧拉法的MATLAB程序，综合考虑计算效率和精度，时间步长设为1×10-8s，根据图3～图5可知，第一阶段发生在0～4.62×10-5s时间范围内，由于破壳压力远小于中心炸药爆炸所产生的爆压，因此，压力在极短时间达到最大值，之后，随着爆轰能量的衰减，压力逐渐衰减到标准大气压。烟幕初始云团加速度迅速减小，当膨胀力小于气动阻力时，加速度成为负值，后缓慢增加；速度先增大后较小，最大半径膨胀到初始半径的2～3倍。在4.62×10-4s以后，烟幕初始云团进入爆炸分散的第二阶段。在第二阶段的短时间内，加速度和速度急剧衰减并趋于0为止。此时，膨胀半径增大到约3 m。

3.2 实验结果与讨论

采用VEO400高速摄影机进行拍摄，截取0～1 s烟幕初始云团膨胀的图像，如图6所示。以1 s这一时刻的图像为例，当背景容易区分情况下，采用图像的二值化，将云团和背景区分开，通过局部灰度拉伸，获取灰度突变的边缘灰度阀值，得到烟幕初始云团边缘形态图像，并求解得到烟幕初始云团最大半径。

图3 云团膨胀加速度随时间变化曲线

图4 云团膨胀速度随时间变化曲线

图5 云团膨胀半径随时间变化曲线

图6 发烟装置静爆试验图

通过上述算法处理得到的烟幕初始云团边缘曲线较为清晰，如图7所示，利用式(9)分别计算各时间节点烟幕初始云团最大半径，得到试验值。同时，分别计算了改进和未改进模型对应时刻云团最大半径并进行了比较与分析，如图8所示，改进模型求解得到的膨胀半径随时间变化曲线比未改进模型更接近试验结果，为进一步证明了改进模型的适用性，计算了二种模型的相对误差曲线，如图9所示。误差随时间逐步趋于稳定，未改进模型误差在10%左右，改进模型的误差约为5%，改进模型计算结果的相对误差较小，相比而言，相对误差降低了约5%～10%。烟幕初始云团膨胀半径，在爆炸发生短时间内，急剧增大。之后，逐渐趋于平缓，其变化规律与试验测试得到的变化规律基本一致。

云团爆炸分散的整个阶段，其理论计算值大于试验测试值，其主要原因在于：爆炸分散的整个阶段，爆轰产物的能量，有一部分被壳体的塑性变形及其他因素所消耗，有效驱动的能量低于理论值，因此，基于一定假设的理论模型计算结果要大于试验结果；试验测试时，高速摄影机图像捕获精度为1×10-3s，第一阶段的整个发生时间小于1×10-3s，捕获精度存在误差，因此，试验测试结果低于理论值。为提高工程计算的适用性，在建立的爆炸分散模型基础上，需增加修正项，对理论模型作进一步修正，降低初始时刻模型的计算误差。

图7 云团边界形态图像

图8 云团膨胀半径随时间变化曲线

图9 相对误差随时间变化曲线

4 结论

通过理论建模、数值计算与实验分析相结合的方法，对圆柱体发烟装置烟幕初始云团运动规律进行了分析。结果表明考虑第一阶段气动阻力的改进模型相比未改进模型，计算值误差较小且更接近试验值，能够描述烟幕初始云团膨胀半径的变化规律。