1 引 言

随着精确制导武器的广泛使用，发展与之对抗的干扰理论及相关技术，提高部队战场生存能力，受到世界各国军队的重视^[1]。烟幕干扰具有成本低、效费比高等优点，是干扰领域发展的重要方向。其中，分析烟幕初始云团的运动规律，计算烟幕初始云团的参数，是烟幕作战效能评估研究的热点之一^[2]。将烟幕初始云团近似看成不断扩张的球体，其最大半径用于表征烟幕遮蔽面积，是烟幕评估中的一项重要指标^[3,4,5]。

为了探索烟幕初始云团运动规律，提高爆炸分散的理论模型精度和适用性，近年来，学者们做了大量工作，如陈宁等^[6]研究了真空环境中烟幕的形成机理，并建立了膨胀模型；朱晨光等^[7]分析了赤磷烟幕初始云团膨胀过程，并指出在膨胀过程中空气阻力影响较大，建立了基于爆轰压力和空气阻力作用下的烟幕初始云团膨胀模型，利用常数变异法建立了初始云团运动速度与半径的表达式；许兴春等^[8]基于等熵膨胀基本理论和球形装药模型，对爆轰压力进行了数值计算，建立烟幕初始云团半径的求解模型，并利用四阶龙格—库塔法得到了初始云团半径的变化规律。陈嘉琛等^[9]利用Fluent仿真软件对燃烧型云雾抛撒过程进行了数值模拟，得到了径向速度和云团范围随时间的变化规律且与实验结果较为吻合。

上述研究针对烟幕初始云团形成及爆炸分散过程运动规律的不同方面进行了研究，但仍存在以下问题：（1）理论模型大多基于一定的假设前提，未考虑壳体膨胀变形和热交换损耗，数值计算精度低；（2）理论模型采用四阶龙格—库塔法求解数值解的效率低，收敛性和平滑性存在不足；（3）限于空气与固体颗粒在有限元建模时，一般采用Euler和SPH方法无法有效耦合，数值模拟方法在爆炸分散型烟幕中应用较少。针对以上问题，为消除假设前提引入的误差，对爆炸分散过程中爆轰压力的计算作进一步改进，采用Truegrid^[10]和Autodyn^[11]仿真建立壳体模型，采用SPH方法对干扰剂和空气进行建模，对中心装药起爆到壳体破裂过程进行仿真模拟，借助于MATLAB^[12]分段线性插值函数，建立爆轰压力随时间变化关系，作为爆炸分散模型求解的已知参量，利用欧拉法^[13]对烟幕初始云团爆炸分散理论模型进行数值计算，并与四阶龙格—库塔法^[14]进行对比分析。通过野外实爆实验，对所提方法进行分析与验证。

2 数值模拟

2.1 爆轰压力数值模拟

2.1.1 计算模型建立

依据特种弹药战斗部装药结构设计原理，建立了一种发烟装置的结构模型，其外形为圆柱体，如图1所示。壳体材料为合金钢，外径62.50 mm，壁厚7.00 mm，高127.50 mm。为使烟幕初始云团成形良好，在壳体四周对称分布了四个“V”型刻槽。刻槽深度4.50 mm，夹角120°。内部装填碳基干扰剂，装填密度2.30 g·cm^-3，装药半径56.50 mm。中心为黑索今（RDX）药柱，装填密度为1.70 g·cm^-3，装药半径9.00 mm，雷管位于RDX药柱上端，用于起爆RDX药柱。

a. longitudinal profile    b. horizontal profile

图1 发烟装置模型结构简图

Fig.1 Structure diagrams of smoke pot model

针对四周刻槽的复杂几何模型，利用Truegrid软件建立1/4壳体网格模型，采用辅助平面映射的方法，建立了“V”型刻槽，如图2所示。网格划分情况见表1。将该网格模型导入Autodyn中，对其余部分采用相同尺寸网格和Langrange算法进行建模计算。碳基干扰剂为粉末物质，且爆炸分散开后，依然呈现粉末形态，因此，采用SPH算法进行仿真建模。限于Euler算法与SPH算法无法耦合的问题，空气采用SPH算法建立仿真计算模型，如图3所示。

图2 Truegrid网格模型

Fig.2 The mesh model of Truegrid

表2 RDX炸药的JWL状态方程参数

Table 2 Parameters of JWL equation of state for RDX explosive

parameter	A / GPa	B / GPa	R1	R2	w	density / g·cm-3
index	669.9	12.9	4.3	1.2	0.3	1.7

图3 Autodyn仿真模型

Fig.3 The simulation model of Autodyn

2.1.2 材料模型及参数

（1） 炸药材料模型及参数采用JWL状态方程对RDX物理化学性质进行描述，JWL方程的p‑V关系^[15]为：

（2） 含碳基干扰剂药柱模型及参数RDX药柱四周为压装的碳基干扰剂，数值计算采用Autodyn自带SAND本构模型，反应材料压实状态的部分利用线性压实状态方程描述，通过测试一定质量的材料，在不同压力状态下的压实密度，建立p‑ρ关系。为有效模拟烟幕云团爆炸分散过程最大半径的膨胀规律，利用锻压系统对含碳基干扰剂药柱的压力与密度关系进行测试。

（3）壳体材料模型及参数外部壳体采用30CrMnSiA合金钢材料，数值计算采用的模型为Johnson‑Cook模型，其本构方程^[15]为：

式中,σ为材料受到的应力，MPa;A为材料的屈服应力，MPa; ε 为材料的应变; ε 0 为材料的参考应变;T_τ为参考温度，K;T_m为材料熔点温度，K;B为应变硬化常量，kPa;n为应变硬化指数;m为热软化指数;C为应变率敏感特性参数。各参数具体数值见表3^[16]。

表3 30CrMnSiA合金钢的参数^[16]

Table 3 Parameters of 30CrMnSiA alloy steel ^[16]

parameter	density / g·cm-3	G / GPa	E / GPa	A / MPa	B / MPa	n	C	M	Tm / K	Tτ / K
index	7.83	77	206	1200	683	0.26	0.014	1.03	1793	294

G is shear modulus. E is elasticity modulus.

2.2 爆炸分散模型建立

发烟装置起爆后，产生爆轰波，并迅速形成大量高温高压气体（中心RDX药柱爆炸形成）向四周膨胀做功，驱动周围的碳基干扰剂向外挤压，壳体内部压力逐渐增大，当超过壳体承受极限后，发烟装置发生爆炸，形成一个近似球体且不断膨胀的高温高压云团。爆炸分散过程实际为爆轰产物在空气介质中流动并驱动碳基干扰剂向外膨胀的过程，烟幕初始云团是爆轰产物与碳基干扰剂形成的混合物。将烟幕初始云团在内部压力衰减为零之前的过程称为爆炸分散的第一阶段，从爆轰压力衰减为零到烟幕粒子速度衰减为零为爆炸分散的第二阶段^[7]。

（1）第一阶段模型建立。取云团边缘一个小微元体，其质量为dm，受爆轰波压力作用面积为dS（微元体体积远小于云团体积），受到的作用力主要是爆轰波压力和气动阻力，则有^[6]：

式中,a为微元体膨胀过程的加速度，可表达为a=d²R/dt²，m·s^-2。p为爆轰波压强，Pa。k′为气动阻力系数。ρ'为标准条件下大气密度，kg·m^-3。v为微元体的速度，m·s^-1。

将第一阶段模型整理得到如下微分方程组：

（2）第二阶段模型建立。在第一阶段，爆轰作用力在云团向外膨胀过程中不断衰减。当压力衰减为零时，进入自由扩散阶段。这一阶段，微元体主要受气动阻力作用，则^[7]：

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

3 实验部分

3.1 实验仪器

上海良平FA‑100电子天平，合肥江淮铸造YH33‑100锻压系统，美国VEO‑400高速摄影机。

3.2 实验过程

（1）含碳基干扰剂药柱的锻压实验。称取40 g碳基干扰剂，装入直径为60 mm的装药模具，并置于压铸机工作台上。根据实验要求，在控制台上设定好下行慢速位、下限位等工艺参数，同时，设定压力分别为20，40，60，80，100，120，140，150 MPa。为确保测试结果的准确性，每个测试点分别进行5 min保压，设置完毕后开始下行，到慢下位时，滑块以设定速度缓慢下行，系统逐渐增压，直到压力传感器检测达到设定值，根据系统设定参数，位移传感器全过程检测滑块的位移。保压时间结束，滑块继续下行，到达下一个测试压力点位置，直到所有测试点测试完毕后，滑块上行结束测试。

（2）发烟装置实爆实验。为有效捕获烟幕初始云团的形成及爆炸分散过程，根据测试要求，在地面风速小于3 m·s^-1的自然条件下，开展了发烟装置野外实爆实验，采用“摄像法”测定烟幕初始云团最大半径的变化规律，图4为实验系统的示意图。实验的基本原理为：利用高速摄像机对烟幕爆炸分散过程进行记录，通过测距仪测量吊架两侧距离。利用图像数据，采用灰度变换、二值化处理、阀值分割、边缘检测方法，获取吊架两侧的像素点坐标和云团边界两侧的像素点坐标，利用二者的比例关系，求解云团径向膨胀的最大半径，其表达式为：

图4 发烟装置实爆实验系统示意图

Fig.4 Schematic diagram of smoke pot explosion testing system

式中,L为吊架两侧的距离，m; Δx₁为图像中吊架两侧像素点距离; Δx₂为图像中云团两侧像素点的最大距离。

4 实验结果分析

4.1 干扰剂本构关系建立

为有效模拟含碳基干扰剂药柱的爆轰压力衰减情况，将试验测得的药柱高度进一步计算得到药柱密度。以压力为横坐标，药柱密度为纵坐标，得到含碳基干扰剂药柱p‑ρ曲线如图5所示。

图5 含碳基干扰剂药柱p‑ρ曲线

Fig.5 The p‑ρ curve of carbon‑containing interference agent

由图5可知，在压力初期，药柱密度增加较快，随压力逐渐增大药柱密度逐步趋于定值。因此，p‑ρ曲线具有对数函数的特点，对压力取对数，得到lnp‑ρ曲线，如图6所示。基于最小二乘法拟合得到碳基干扰剂的压力与密度关系式如下：

图6 含碳基干扰剂药柱lnp‑ρ曲线

Fig.6 The lnp‑ρ curve of carbon‑containing interference agent

4.2 烟幕初始云团运动规律分析

（1）爆轰压力数值模拟。选择仿真模型中靠近刻槽附近的多个干扰剂粒子作为参考点，为便于观察刻槽位置破裂情况，将空气模型暂时隐藏，仿真中某一时刻的仿真形态如图7所示。通过仿真得到烟幕初始云团压力衰减曲线，如图8所示。

图7 发烟装置某一时刻的仿真形态

Fig.7 The simulation form of a smoke device at some time

图8 参考点初始爆轰阶段p‑t曲线

Fig.8 The p‑t curve of reference point initial detonation state

从图8可知，直到约7×10^-6 s，压力值发生突变，这一过程为中心RDX药柱发生爆轰反应并驱动含碳基干扰剂药柱向四周压实的过程。之后短时间内压力急剧增大，壳体合金钢材料达到屈服极限后，迅速破裂，发烟装置发生爆炸，压力在较短时间增大到最大爆轰压力，并迅速衰减为零。因此，与基于等熵膨胀理论的爆炸分散模型相比，基于Autodyn仿真考虑了干扰剂压实和壳体塑性膨胀过程，为了减少插值节点数和误差，借助于MATLAB提供的分段线性插值函数，以1×10^-9 s为时间步长，建立第一阶段模型的p‑t关系。

（2）数值计算与分析。式（4）和式（6）为二阶非线性微分方程组，可采用数值积分法求得数值解。常用的数值积分法包括欧拉法和四阶龙格‑库塔法等。一般针对低阶微分方程，为确保数值解的光滑性，采用欧拉法进行计算。其截断误差为O（h²），计算公式如下：

为进一步计算模型的数值解，首先需确定初值。根据中心装药的体积和装药密度等参数，根据文献[17]，计算爆压：

式中,p'为炸药爆轰压力，GPa；ρ_D为炸药装药密度，g·cm^-3；∑N为炸药氮当量。

通过式（10）计算得到瞬时爆轰压力p'=27.76 GPa。而30CrMnSiA壳体材料破裂的强度为1.20 GPa，因此，可确定求解理论模型爆炸瞬时的初值为：p₀=1.20 GPa，R₀=56.50 mm。

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

a. first stage

b. second stage

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

Fig.9 History of smoke cloud expansive acceleration

a. first stage

b. second stage

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

Fig10 History of smoke cloud expansive rate

a. first stage

b. second stage

图11 云团膨胀最大半径随时间变化曲线

Fig.11 History of smoke cloud maximum expansive radius

为进一步证明欧拉法的适用性，编写了四阶龙格‑库塔法MATLAB程序，并与欧拉法进行对比分析，如图12所示。由图12可知，采用较小时间步长时，四阶龙格‑库塔法计算效率低，得到的数值解短时间内难以收敛，选择较大步长时，误差较大，不同时间步长得到的数值解相差较大；相比图11而言，采用欧拉法时间步长为1×10^-9 s，计算精度和效率较为合理，更利于此类问题的求解。

图12 基于四阶龙格‑库塔法云团膨胀半径随时间变化曲线

Fig.12 History of smoke cloud expansive radius based on the fourth order Runge‑Kutta method

4.3 实验结果与讨论

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

图13 发烟装置实爆试验图

Fig.13 Test charts of smoke equipment explosion

通过4.1节的图像处理算法得到的烟幕初始云团边缘曲线较为清晰，如图14所示，利用式（7）分别计算各时间节点烟幕初始云团最大半径。同时，采用欧拉法，分别对传统理论计算和基于Autodyn仿真与理论相结合的方法计算得到对应时刻云团最大半径并进行了比较与分析，如图15所示。由图15可知，基于Autodyn与理论计算相结合的改进方法求解得到的膨胀半径随时间变化曲线比传统理论计算^[9]结果更接近实验结果，烟幕初始云团膨胀半径，在爆炸发生短时间内，急剧增大。之后，逐渐趋于平缓，其变化规律与试验测试得到的变化规律基本一致。为进一步证明改进方法的适用性，计算了两种方法的相对误差曲线，如图16所示。误差随时间增长逐步趋于稳定，且在0.2 s以后，均小于10%；基于Autodyn仿真与理论相结合的方法计算结果相对误差较小，相比传统理论计算结果，误差降低了3%~8%。

图14 烟幕初始云团边界形态图像

Fig.14 Image of smoke cloud boundary form

图15 烟幕初始云团膨胀最大半径随时间变化曲线

Fig.15 History of smoke cloud maximum expansive radius

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

Fig.16 History of relative error

[5] 姚禄玖, 高钧麟, 肖凯涛, 等.烟幕理论与测试技术[M]. 北京：国防工业出版社, 2004：97-104.[6] 陈宁, 潘功配, 陈厚和, 等.真空环境中烟幕云团形成阶段的膨胀模型[J]. 火工品, 2006(1): 1-5.[7] 朱晨光, 潘功配, 关华, 等.烟幕云团形成初期的流动规律研究[J]. 含能材料, 2007, 15(5): 540-543.[8] 许兴春, 高欣宝, 李天鹏, 等.烟幕初始云团半径变化规律理论模型及实验研究[J]. 爆炸与冲击, 2016, 36(2): 183-187.[9] 陈嘉琛, 张奇, 马秋菊, 等.固体与液体混合燃料抛撒过程数值模拟[J]. 兵工学报, 2014, 35(7): 973-976.[13] 庄小兰, 王琦.一阶时滞微分方程欧拉法的动力性[J]. 广东工业大学学报, 2018, 35(1):47-49.[14] 沈艳, 张丽玲, 张琦智, 等. 基于三阶和四阶龙格库塔法的GM（1,1）模型优化及应用[J]. 数学的实践与认识, 2016, 46(7):169-173.[15] 石少卿, 汪敏, 孙波, 等. Autodyn工程动力分析及应用实例[M]. 北京：中国建筑工业出版社, 2012: 38-54.[16] 夏开文, 刘文彦, 唐志平.30CrMnSiA钢高温动态力学性能的实验研究[J]. 爆炸与冲击, 1998, 18(4): 311-316.[17] 刘玉存, 王作山, 吕春玲, 等. 炸药的输出钢凹值与爆压的关系研究[J]. 华北工学院学报, 2001, 22(4): 304-307.

YAO Lu‑jiu, GAO Jun‑lin, XIAO Kai‑tao, et al. Theory and testing technigue of smoke[M]. Beijing: National Defence Industry Press, 2004： 97-104.CHEN Ning, PAN Gong‑pei, CHEN Hou‑he, et al. Expansive model of smoke cloud forming course in vacuum[J]. Initiators & Pyrotechnics, 2006(1): 1-5.ZHU Chen‑guang, PAN Gong‑pei, GUAN Hua, et al. Initial flow ability of smoke cloud forming[J]. Chinese Journal of Energetic Materials （Hanneng Cailiao）, 2007, 15(5): 540-543.XU Xing‑chun, GAO Xin‑bao, LI Tian‑peng, et al. Theoretical model and experiment of radius variation of initial smoke cloud[J]. Explosion and Shock Waves, 2016, 36(2): 183-187.CHEN Jia‑chen, ZHANG Qi, MA Qiu‑ju, et al. Numerical simulation of dispersal process of solid‑liquid mixed fuel[J]. Acta Armament‑rii, 2014, 35(7): 973-976.[10] Concannon B. Truegrid[CP]. XYZ Scientific Applications, 2001.[11] Naury B. Autodyn[CP]. Century Dynamics, 1985.[12] Cleve M. MATLAB[CP]. Mathworks, 2014.ZHUANG Xiao‑lan, WANG Qi. Dynamics of Euler method for the first order delay differential equation [J]. Journal of Guangdong University of Technology, 2018, 35(1):47-49.SHEN Yan, ZHANG Li‑ling, ZHANG Qi‑zhi, et al. Optimization and its application for GM(1,1) model based on the third and the fourth order Runge‑Kutta method[J]. Mathematics in Practice and Theory, 2016, 46(7):169-173.SHI Shao‑qing, WANG Min, SUN Bo, et al. Autodyn engineering dynamic analysis and application examples[M]. Beijing: China Building Industry Press, 2012: 38‑54.XIA Kai‑wen, LIU Wen‑yan, TANG Zhi‑ping. Experimental study of dynamic properties of 30CrMnSiA steel at high temperature[J]. Explosion and Shock Waves, 1998, 18(4): 311-316.LIU Yu‑cun, WANG Zuo‑shan, LÜ Chun‑ling, et al. The relation between steel dent data and detonation pressure of explosive sample[J]. Journal of North China Institute of Technology, 2001, 22(4): 304-307.

5 结 论

基于圆柱体发烟装置，采用理论建模、仿真模拟、数值计算与实验验证相结合的方法，对烟幕初始云团运动规律进行了分析与研究，得到如下结论：

（1） 基于含碳基干扰剂药柱的锻压实验，建立了干扰剂的本构关系，为基于Autodyn含碳基干扰剂药柱的仿真建模提供依据。采用SPH方法建立含碳基干扰剂药柱模型和空气模型，模拟爆炸分散过程中空气阻力对烟幕初始云团的作用，有效解决了采用SPH方法建立的干扰剂模型与Euler方法建立的空气模型无法耦合的问题。

（2） 采用Truegrid与Autodyn混合仿真，模拟烟幕初始云团爆炸分散过程中的爆轰压力，作为已知量代入爆炸分散模型中，并与理论计算相结合求解烟幕初始云团的最大半径，此方法可降低基于等熵膨胀理论假设的单一理论计算方法引起的误差。所提出的计算方法与传统理论计算方法相比，烟幕初始云团最大半径随时间变化曲线的计算误差降低了3%~8%，且膨胀半径变化规律与试验情况基本一致。

（3） 采用欧拉法计算理论模型的数值解，相比四阶龙格‑库塔法，收敛性好、计算效率高，更适合此类问题的求解。

参考文献