单兵防暴喷射器是一类将刺激剂喷射到有生目标，使其暂时失去抵抗能力的单兵自卫非致命武器，与动能打击类非致命武器相比，不存在致命、致伤风险，实战应用更加灵活，在警卫、守卫、巡逻等勤务中能有效增强执法人员的个体防护和现场处置能力，具有广阔的应用前景。而现有单兵防暴喷射器多采用二元内腔结构，以高压气体为动力，平时高压贮存，容易出现喷头堵塞、液体泄漏等问题，给临战使用和贮存保管带来不便，限制了其作战效能的发挥。
对此，本文采用非贮压喷射技术，设计了一种新型防暴喷射管，其结构如图1所示。该喷射管以火药为动力源，战剂平时装填在发射管内，能有效解决贮存保管问题，勤务性能大为提高。然而，要使其具有良好的喷射性能，就需要对其喷射过程进行深入研究。
图1 某型防暴喷射管结构示意图
采用CFD的方法实现流场的可视化研究，具有成本低、周期短、能模拟出射流发展细节等优点，是目前常用的科学研究方法[1]。陆朝晖[2]就借助该方法深入研究了冲击挤压式、截断式脉冲射流形成过程的流场变化特点；刘来国[3]针对自激喷嘴建立数学模型，揭示了自激喷嘴对射流产生的增压特性，并对其结构进行了优化设计；张德荣等[4]借助Fluent软件对水力喷砂射孔器喷嘴结构进行优选，来指导其结构优化。贾凯等[5]为解决现代火炮身管膛线加工排屑不畅的问题，通过Fluent仿真试验，对所设计的新型双曲线射流通道进行优化，得出了最优解。上述研究都是单独针对射流的管内流场或管外流场进行分析研究，缺乏对喷射过程的整体考量。课题组前期在某脉冲防暴水炮近流场[1]的研究中，就采用了以模拟管内流动为管外射流提供初始条件的方法，揭示了其发展过程。然而，本文的设计涉及到了更多物理场，以往的数值模型不适用于本文的研究。
因此，本文拟采用多物理场耦合[6]的方法建立数值计算模型来深入分析该喷射管喷射过程，以期为下一步火药选型、喷射管和喷嘴结构优化等提供基础和借鉴。
1 工作原理分析
该喷射管由底火、上下隔板、活塞、战剂、喷嘴及管体组成，其工作原理是：当底火被刺燃后，底火内火药瞬间产生高温高压气体，达到一定压力后冲破下隔板，推动活塞使管内的战剂受压，通过喷嘴的整型作用从喷嘴喷出形成射流，射流在空气阻力作用下逐渐雾化，喷射到目标。
可以看出，其喷射过程涉及到火药燃烧、高温高压气体膨胀、隔板破裂、活塞运动、射流流动、雾化多个复杂的物理过程，而其中火药能量、活塞速度、射流变化是关注的重点问题，因此在数值模拟过程中需要在合理假设的基础上进行简化，以便于深入分析和研究。
2 数值计算模型
2.1 基本假设
1) 火药充分燃烧完毕后冲破下隔板；
2) 火药气体膨胀过程满足理想气体状态方程，且不考虑传热；
3) 活塞为线弹性材料，与管壁无摩擦，且在运动过程中无变形；
4) 不考虑上、下隔板的破裂过程；
5) 战剂为不可压缩牛顿流体。
基于以上假设，该喷射管喷射过程简化为高压气体膨胀、活塞推动以及射流喷射三个过程，其数值计算模型可由高压气体的压力控制模型、基于流固耦合的管内流动模型和基于水平集方法的射流喷射模型构成。
2.2 高压气体的压力控制模型
由基本假设1)和2)可知，火药气体膨胀满足理想气体状态方程，即：
PV=nRT
(1)
其中: P为理想气体的压强；V为理想气体的体积；n为气体物质的量； T为热力学温度； R为气体常数。
由于不考虑传热，则在气体膨胀过程中P0V0=P1V1，故
P1=P0V0/V1=P0V0/(V0+S×h)
(2)
其中： P0为火药燃烧后初始压强；V0为气室体积；P1为火药气体膨胀过程中的压强；V1为火药气体膨胀过程中的体积；S为喷射管截面积； h为活塞的位移。
2.3 基于流固耦合的管内流动模型
活塞推动战剂的管内流动过程为流固耦合过程，由牛顿第二定律，活塞运动的动量平衡方程可描述为[6]：
(3)
式中： ρ1为材料密度；usolid为活塞运动位移场；S为第二类Piola-Kirchhoff应力张量; F为变形梯度张量； Fv表示体积力； I是单位张量。
由基本假设(5)可知，战剂的管内流动可考虑为不可压缩层流，其N-S方程和连续性方程如式(4)所示。
(4)
式中：u2为战剂速度场；p为战剂压力； ρ2为战剂密度； μ是战剂动力粘度； Fext是作用在战剂上的外力。
流固耦合发生在活塞与战剂接触的边界上，使用任意拉格朗日-欧拉(ALE)法，将使用欧拉描述和空间坐标系的流体流动[式(4)]和使用拉格朗日描述和材料坐标系表示的固体力学[式(3)]组合进行求解。
2.4 基于水平集方法的射流喷射模型
战剂的管外喷射为不可压缩两相湍流过程。在多种湍流模型中，由于Realizable k-ε模型在圆口射流模拟中能获得更好的湍流比率精度[7]，适合于本文的研究，因此，采用该模型来描述外部流场，即
(5)
式中，
在上述方程中，湍流粘性系数
是平均应变率与旋度的函数，Pk表示由于平均速度梯度引起的湍动能产生，u为流体速度场，p为流体压力，ρ为流体密度，μ是流体动力粘度，Fext是作用在流体上的外力。作为默认常数，c2=1.9，湍动能k与耗散率ε的湍流普朗特数分别为σk=1.0，σε=1.2。
为追踪射流界面，本文采用水平集方法[8]，即：
(6)
其中，φ是平滑阶跃函数，在气液交界面，其值从0到1平滑过渡。γ为重新初始化参数，取1 m/s，εls为界面厚度控制参数。
流体的密度ρ和动力粘度μ由下式给出：
(7)
其中，ρ气、ρ2分别是空气和战剂的密度; μ气、μ2分别是空气和战剂的动力粘度。
式(2)～式(7)即为本文所建立的喷射管喷射过程数值计算模型。为保证数值模型的鲁棒性，其实现过程采用先计算管内流动过程，得出射流瞬时出口平均速度，而后以此为初始条件进行管外两相流动计算的方式进行。
3 数值计算过程
如前所述，该喷射管的喷射过程涉及多个物理场，Comsol Multiphysics平台不仅提供了多种物理场模块，而且具有定义和耦合任意数量偏微分方程的能力，是当今国际上主流的多物理场耦合数值模拟软件[9]，因此，本文基于该平台，采用有限元法进行数值计算。
3.1 物理模型的建立和网格划分
为减少对计算能力的需求，该喷射管的喷射流场可简化为二维轴对称模型，其几何结构如图2所示。
图2 喷射管物理模型示意图
在Comsol中采用“用户控制网格”进行网格划分，其中管内部分采用四边形网格，最小单元大小设定为0.007 5 mm，最大单元大小设定为0.65 mm，管外部分则采用自由三角形网格，网格总量约为15万个单元，网格平均质量为0.968 7。其中喷嘴内外局部网格划分如图3所示。
图3 网格划分局部示意图
3.2 管内流动过程计算
管内流动过程的数值计算由动网格、固体力学、层流、流固耦合四大模块组成。
1) 动网格模块
动网格控制流体区域，其变形域指定为活塞运动的区域，采用Yeoh平滑类型；其法向网格位移边界指定为活塞运动经过的管壁。
2) 固体力学模块
固体力学模块中，刚性域设置控制活塞，并指定活塞只能横向运动，其初始位移、速度均设置为0，其载荷边界选定为与气室接触的边，载荷类型定义为公式(2)中的P1。
3) 层流模块
层流域指定为活塞与喷嘴喷口之间的流体域，初始速度场和压力均设置为0，指定喷射管管壁为无滑移壁面，喷口处设置为开放边界，其法向应力设置为0。
使用Comsol的流固耦合模块，指定流体和结构耦合结构，采用全耦合类型。将气室内初始压强P0，即火药燃烧后初始压强，设定为10 MPa，战剂属性设置为Water，活塞材料设置为Steel。求解时间步设置为0.0000 1 s，速度场、压力和空间网格位移计算采用PARDISO求解器，位移场计算采用MUMPS求解器，对管内流动过程进行瞬态计算。
3.3 管外流动过程计算
管外流动过程的数值计算由湍流、Level-Set模块、两相流模块组成。
1) 湍流模块
指定整个外部流场为湍流计算域，初始速度场和压力均设置为0，指定外流场边界为开放边界，将3.3计算得到的喷口处流体速度指定为入口边界条件。
2) Level-Set模块
Level-Set模块的计算域为整个外部流场域，并将其初始相指定为空气，设置喷口处为入口边界条件，初始相指定为战剂，外流场的其余边界指定为出口边界条件。
使用Comsol两相流模块，设定相应的流体流动和移动截面接口。战剂和气体材料分别设置为Water和Air，求解时间步设定为0.000 1 s，水平集变量、速度场、压力、湍流变量均采用PARDISO求解器，对管外流动过程进行瞬态计算。
4 结果分析
4.1 管内流动过程分析
1) 活塞运动分析
图4为活塞运动速度及位移随时间变化曲线。可以看出，活塞速度在气体膨胀推动作用下急剧增大，在0.31 ms内达到最大值8.81 m/s，而后受管流阻力的影响迅速下降并趋于平缓，活塞在18.47 ms到达喷嘴处，此时速度约为2.02 m/s。
活塞的动力来自于火药气体膨胀做功，气室内初始压强越大，活塞运动速度越快，运动时程也就越短，赋予射流的能量也就越大。当然这也会增加对活塞及喷射管管壁的强度要求，这就需要综合考虑射流速度、射程等指标要求确定火药初始压强，进而对火药进行选型设计。
图4 活塞速度及位移随时间变化曲线
2) 管内压力变化分析
图5为5 ms、10 ms、15 ms和18 ms时刻管内战剂压力分布云图。可以看出：在活塞冲击作用下，喷射管内战剂产生压力，以弹性波的形式向喷嘴传播[10]。由于喷嘴结构的突变，战剂在由喷嘴短圆柱段进入长圆柱段时形成压力梯度，由于射流的卷吸作用，在喷嘴长圆柱段末端出现负压。
图5 管内战剂压强分布云图
图6为喷射管内和喷嘴圆柱段入口处战剂平均压强随时间变化曲线。由于火药压力瞬间释放推动活塞运动，喷射管内战剂压强达到最大值约9.12 MPa后迅速降低并趋于平缓，喷嘴入口处压强也呈现类似的变化趋势，由于喷嘴直径的减少，战剂流速增加，使得喷嘴入口压强比喷射管内压强要低，其最大值约4.61 MPa。
图6 战剂压强-时间变化曲线
3) 管内战剂流速分析
图7为5 ms、10 ms、15 ms和18 ms时刻管内战剂流速分布云图。可以看出：活塞的推动赋予了战剂流动的动力，在每一时刻，喷射管内战剂流动速度比较均匀，喷嘴直径要小于喷射管直径，使得战剂在喷嘴处加速，在喷嘴的短圆柱段形成速度梯度，再经过长圆柱段的整形加速从喷口喷出。
图7 管内战剂速度云图
图8为喷嘴长圆柱段出入口的速度随时间变化曲线。在喷射过程中，战剂速度呈现与活塞运动一致的趋势：先急剧增加而后迅速降低并趋于平缓。通过喷嘴的整形加速，战剂出口速度平均增加了约8.7 m/s。显然，在一定长度范围内，随着长圆柱段的加长，出口速度增加的也就越多。超出了这个范围，随着管流阻力的增加，喷嘴的增速作用会减弱，这就需要对喷嘴的结构进行优化，使得在同样条件下战剂速度达到最优。
图8 战剂速度-时间变化曲线
4.2 管外喷射过程分析
图9和图10给出了不同时刻射流气液体积分数和速度分布云图。射流的发展过程可以分为3个阶段：
1) 伞状结构阶段。在射流初期，低速射流喷出喷口，在后续高速射流的推动及前端空气摩擦阻力[10]的作用下侧向发展偏离中心射流，在失去加速动力后相对中心区域向后滑移，在0.5 ms左右形成伞状结构。
2) 水团形成及破裂阶段。随着射流出口速度的逐渐降低，前端的高速射流向前运动，由于空气的摩擦作用，表面出现波纹，其波幅逐渐增大，约在1 ms时形成大块水团；其后端核心部分呈紧密状态。而后水团吸入空气，核心部分断面积越来越小，水团约在2 ms时破碎，并随着吸入空气量的增多进一步变成水滴。
3) 射流继续发展阶段。管内战剂以较低的速度继续从喷口喷出直至喷完，在管外形成细长射流。可以预见，射流在继续运动过程中，由于空气摩擦阻力的作用，逐渐雾化形成液滴。
图9 不同时刻射流气液体积分数分布云图
图10 不同时刻气液速度分布云图
5 结论
1) 针对某型防暴喷射管高压气体膨胀、活塞推动和射流喷射3个过程分别建立了数值模型，采用多物理场耦合的方法实现了模拟计算，对脉冲射流的形成和发展研究提供了可视化的计算分析方法。
2) 管内流动分析表明：在初始压力在10 MPa的条件下，活塞在0.31 ms内速度达到最大值8.81 m/s，在18.47 ms时刻运动到喷嘴处；管内战剂的压力以弹性波的形式向喷嘴传播，经过长圆柱段的整形加速从喷口喷出，最大出口速度为70 m/s。
3) 管外喷射过程可分为伞状结构、水团形成及破裂和射流继续发展3个阶段，空气摩擦阻力是射流破碎的主要因素。
[1] 赵法栋,庄弘炜,战仁军,等.脉冲防暴水炮近流场的数值模拟[J].计算机工程,2011,37(15):274-275，282.
[2] 陆朝晖.高压脉冲水射流流场结构的数值模拟及破硬岩机理研究[D].重庆：重庆大学,2012.
[3] 刘来国.基于CFD的振荡脉冲水射流压力特性及自激喷嘴的结构设计研究[D].北京：北京工业大学,2014.
[4] 张德荣,吴思梦.基于FLUENT的水力喷砂射孔器喷嘴优化设计[J].石油机械,2017,45(12):71-77.
[5] 贾凯,关波,刘修配,等.负压装置喷射流通道结构优化设计仿真[J/OL].兵器装备工程学报:1-4[2019-09-29].http://kns.cnki.net/kcms/detail/ 50.1213.tj.20190611.1418.005.html.
[6] 孙培德,杨东全,陈奕柏.多物理场耦合模型及数值模拟导论[M].北京：中国科学技术出版社,2007.
[7] ATANOV G A,KOKIN V V.Control of the jet of a pulsed flywheel water cannon[J].Journal of Mathematical Sciences,1995,76(3):2424-2427.
[8] OLSSON E,KREISS G.A Conservative Level Set Method for Two Phase Flow[J].J.Comput.Phys.,210:225-246,2005.
[9] 中仿科技公司.COMSOL Multiphysics有限元法多物理场建模与分析[M].北京：人民交通出版社，2007.
[10] REHBINDER G.Investigation of water jet pulses generated by an impact piston device[J].Applied Scientific Research,1983,40(1):7-37.
[11] 李晓红,卢义玉,向文英.水射流理论及在矿业工程中的应用[M].重庆：重庆大学出版社,2007.