爆炸驱动燃料云团界面不稳定性发展及颗粒分散过程

申 洋， 丁 珏， 翁培奋

(上海大学上海市应用数学和力学研究所， 上海200072)

可燃气体、液固燃料空气混合物的爆炸性质对工业生产安全造成严重危害， 如金属颗粒在工厂发生的粉尘爆炸， 但同时因其爆炸毁伤效应而在国防和军事有巨大用途. 利用中心装药爆炸驱动燃料运动形成云团， 这样的燃料空气混合物， 称作燃料空气炸药(fuel air explosive，FAE)[1]. FAE 是以挥发性液体碳氢化合物或固体粉质可燃物为燃料， 以空气的氧气为氧化剂组成的非均相爆炸性混合物， 具有能量高、分布爆炸性、原料易得等特点.

燃料的爆炸驱动是提高燃料空气炸药装置威力的关键技术， 国内外学者在爆炸作用下FAE 装置壳体破裂[2-4]、燃料界面发展[5-6]及射流形成[7-9]、液体燃料运动及破碎[10-13]、FAE云团形成的规律性[14]等方面开展了很多有意义的研究.

李磊等[15]设计了一种有上下平面约束的液体爆炸分散装置， 应用平面激光诱导荧光技术得到了径向膨胀液体环内诱导的荧光照片. Kandan等[16]通过激波管实验和数值计算， 研究了冲击波加载下燃料不稳定性的发生和发展. Shi等[17]在垂直矩形激波管中， 利用半导体激光器、信号放大电路、数字示波器、高速摄影， 以及自行设计的测量系统测量了激波驱动的界面速度;并对Richtmyer-Meshkov 不稳定性(Richtmyer-Meshkov instability， RMI)的机理进行了实验研究. 结果表明: 界面速度随马赫数的增加而增大， 界面运动与混合区宽度的距离与时间呈线性关系.

Frost等[18]基于实验和数值方法研究了爆炸驱动下颗粒到达不同位置时的速率， 分析颗粒和周围流场动量和能量之间的交换. 廖斌等[19]采用实验与数值计算方法对冲击作用下液滴在环境液体中的演变过程进行了研究， 重点讨论液膜破碎之前液滴的4 种演变模式和流场.

发生两相爆轰的条件如下: 空间悬浮一定浓度的燃料， 形成燃料云团; 而爆炸驱动液体燃料扩散是液固燃料空气混合物形成的重要过程. 以往对于爆炸驱动燃料分散的研究主要集中于云团形状和状态， 而对于爆炸驱动最初阶段的力学性质， 由于受到实验或试验条件的限制，研究成果较少， 特别是涉及波、涡与燃料颗粒的相互作用机制、气液界面的不稳定演化规律尚未清晰. 因此， 本工作开展液体燃料在爆炸冲击波作用下分散机制的研究， 为工业生产的安全防治以及国防领域燃料空气炸药威力的提高提供理论指导.

1 爆炸驱动燃料云团运动的物理数学模型和数值计算方法

在工业生产、国防军事领域， 冲击波驱动燃料使得颗粒处于高温、高压和高速气流中， 速度和温度差导致两相间动量和能量的传递; 同时， 颗粒发生蒸发效应， 使两相间出现质量传递，逐渐形成由可燃气体、燃料与空气混合的云团. 本工作将重点研究爆炸驱动下流场和燃料云团的运动性质.

爆炸驱动装置通常为柱对称结构， 通过中心轴装填的高能炸药爆炸释放的能量来驱动燃料(见图1). 图2 显示燃料质量为12.5 kg 爆炸驱动装置在58 ms 时刻云团的状态.

图1 爆炸驱动装置Fig.1 Explosive driven device

图2 爆炸驱动燃料云团的形成Fig.2 Fuel cloud formation by explosion-driven

1.1 燃料云团湍流运动的物理数学模型

针对爆炸驱动燃料云团的形成和运动过程， 本工作建立相间耦合的两相流模型来描述， 即把气体当作连续介质(流体相)， 把燃料作为离散体系(颗粒相)， 用欧拉-拉格朗日方法来描述，并考虑液体燃料的蒸发效应.

1.1.1 流体相控制方程

本工作讨论的是非定常的流场， 二维多相流的控制方程为

式中：Φ 为通用变量， 分别表示 u， v， k， ε， H 和 Y 等求解变量; ρ 为密度; r 为柱坐标中的径向; Uj为j 方向的速度; Γ 为扩散系数; 在直角坐标和柱坐标下， 上标m = 0 时用于平面流动，m = 1 时用于轴对称流动; SΦ和SPΦ分别表示气相本身的项和液体燃料颗粒蒸发引起的额外源项.

为了分析气液界面不稳定的发展， 刻画冲击波诱导旋涡产生的过程， 本工作采用成熟的、适用于高雷诺数流场的标准k-ε 湍流模型， 数值研究气体的湍流性质.

气相遵守理想气体状态方程，

1.1.2 颗粒相控制方程

假设燃料颗粒的初始粒径相同， 且颗粒间的相互作用可以忽略. 控制方程为

1.1.3 液体燃料蒸发的动力学方程

本工作讨论的燃料为庚烷. 颗粒在波后加速， 由于对流传热升温发生蒸发， 形成燃料蒸气.根据颗粒蒸发模型， 其质量变化关系式为

式中: Sc 是施密特数; Sh 是Sherwood 数， 表示对流作用对传质过程的影响，

B 为 Spalding 数， 也称质传递数，

其中Ys表示颗粒表面的蒸汽质量分数， Y 为燃料蒸汽在气相中的质量分数.

颗粒表面温度变化率为

式中: Q 代表气相传递给液相燃料的热量; Nu 为努塞尔数; Ts为颗粒表面温度; T∞为环境温度; Lv为颗粒气化潜热; cp，l为液体等压比热容; D 为颗粒直径; λ 为导热系数.

1.2 燃料云团湍流运动的数值计算方法

本工作采用欧拉-拉格朗日方法， 利用压力隐式算子分裂(pressure implicit split operator，PISO)算法和压力耦合方程组半隐式方法(semi-implicit method for pressure linked equations，SIMPLE)相结合的PIMPLE 算法， 其中对流项用二阶精度的守恒型单调迎风格式(monotonic upwind scheme for conservation laws， MUSCL)， 数值研究爆炸冲击波的驱动以及燃料云团的形成过程.

为了验证数值格式的精度， 针对Sod 问题开展数值验证[20]. Riemann 类型的初始条件如下: (ρ，v，p)L=(1，0，1); (ρ，v，p)R=(0.125，0，0.100).

图3 和4 给出了t=2 时刻压力和密度数值解与精确解的对比. 由图可以看出， 数值解与精确解吻合较好， 数值格式对波和间断具有较高的分辨率.

为了刻画云团的形成性质， 开展对爆炸驱动12.5 kg 燃料云团运动过程的数值模拟. 研究结果显示， 液体燃料分散是云团在空间逐渐形成和展开的动力学过程， 即中心装药发生爆轰， 而后爆轰产物气体作用于周围介质形成冲击波， 波和爆轰产物驱动燃料运动. 在驱动初期， 爆炸作用力占主导; 随后的阶段， 气动阻力占主导， 逐渐形成燃料和空气混合的可爆性云团. 图5 和6 给出了云团尺寸的计算结果与实验数据的比较. 可以看到: 在60 ms 时间内， 计算所得的云团演化尺寸与实验数据非常接近， 即在爆炸作用力下， 云团尺寸发展较快; 当云团内的压力接近外界环境压力时， 云团膨胀运动变慢， 云团尺寸的变化趋缓. 计算所得的云团尺寸随时间演化满足指数函数规律: 半径R = 5.026 46-5.137 63e-0.079t; 高度H = 1.774 59-1.853 92e-0.084t(t : 0 ～60 ms). 按照云雾膨胀的相似率[21]， 函数中随时间t的增长系数(云团半径增长系数为0.079， 云团高度增长系数为0.084)反映了中心装药释放的总能量， 与爆炸驱动力密切相关.

图3 Sod 问题的压力曲线Fig.3 Pressure curve for Sod problem

图4 Sod 问题的密度曲线Fig.4 Density curve for Sod problem

图5 燃料云团半径随时间的变化Fig.5 Fuel cloud radius over time

图6 燃料云团高度随时间的变化Fig.6 Fuel cloud height over time

2 爆炸驱动下冲击波、湍流涡与燃料之间的相互作用

在爆炸驱动最初期， 冲击波作用于燃料， 气液界面发生演化， 由初始较为细微的扰动逐渐发展成湍流. 同时， 波与液体燃料分界面作用发生折射， 气液界面获得加速， 促进不稳定性发展， 这些对燃料爆炸抛撒、燃料混合起着重要作用.

下面将着重研究和分析最初期两相耦合的动力学性质. 图7 给出了计算区域示意图， 其中区域尺寸为 170 mm×170 mm(长×高). 为了简化， 冲击波的效应采取10 mm×170 mm(长×高)高温高压(T =1 000 K，P =1 519.875 kPa)区域来代替.Γ1为对称平面，Γ2为区域的对称轴，Γ3和Γ4为无反射边界条件.燃料为液态庚烷，质量为4.242 mg. 在数值模拟中，为了比较和分析数值研究的网格，跟踪流场中的点(0.06 m，0.04 m)处的参数. 当计算的网格尺寸分别为0.5 和0.8 mm 时，流场跟踪点处压力峰值分别为578.404 和568.707 kPa，与网格尺度1 mm(压力峰值567.106 kPa)分别相差1.70%和0.77%. 综合考虑计算量和误差， 采取尺寸为1 mm 的网格系统， 数值模拟中网格数目为28 900. 图8 给出了燃料云团的示意图， 其中蓝色箭头代表诱导涡旋转方向为顺时针， 而红色箭头代表逆时针.

图7 计算区域Fig.7 Computation field

图8 燃料云团示意图Fig.8 Schematic diagram of fuel cloud

2.1 两相介质中冲击波的性质

2.1.1 冲击波的动力学性质

波的传播是随时间演化的非定常运动过程， 波后介质会发生动力学和热力学状态的变化.当高温高压气体产生的冲击波作用于气相和燃料介质中， 在介质内存在能量传递过程. 压缩作用也由近及远地传递到介质其他部分. 在传递过程中， 扰动部分与未扰动部分的分界面是波的波阵面. 由于空气和燃料介质中波传播时的阻抗不同， 波传播速度不同， 因此波在两相中阵面发生弯曲(见图9).

图9 冲击波在两相介质中的传播Fig.9 Transmission of shock waves in two-phase media

图10 给出了在燃料介质和气相介质中冲击波压力的分布情况. 可以看到: 燃料介质中压力分布存在2 个峰值(peak1 和peak2)， 而气相介质中只存在1 个压力峰值. 分析原因主要是波在燃料中传播时发生反射， 反射波朝爆源运动时， 由于波的汇聚导致压力升高， 在爆源附近又形成一个新的波. 随后， 这个新生成的波开始在燃料介质中向前传播. 此外， 波扫过燃料， 燃料介质的动量和能量发生变化. 在40 μs 时刻前， 燃料云团中波的压力峰值较空气中更大. 由于稀疏波的影响， 压力峰值随时间逐渐降低. 在105 μs 时刻， 燃料云团中压力峰值的超压值降至0.21 MPa 左右.

图10 两相介质中冲击波压力峰值在空间的分布Fig.10 Spatial distribution of shock wave pressure peaks in two-phase media

气相和燃料介质中波的超压峰值随时间的演化如图11 所示. 选取波与燃料云团开始的21 μs 为初始作用时刻， 波即将穿过云团的105 μs 为结束时刻. 计算结果显示: 波在气相场中传播时， 由于受到燃料介质的影响， 在 21～105 μs 时间段内， 压力峰值相对稳定， 衰减较为缓慢; 然而， 在 21～105 μs 时间段内， 燃料云团介质内超压发生快速衰减. 第一道波在 66 μs 时刻后穿过云团， 第二道波于105 μs 时刻后离开云团.

图11 气相和燃料介质波超压峰值的演化Fig.11 Evolution of overpressure peaks in gas and fuel medium

冲击波在气相、燃料介质中超压峰值的演化均满足幂数律的关系，

式中: a， b， c 为系数; Z(m·kg-1/3)为比例距离， 由冲击波传播的距离和中心装药量的立方根来确定[22].

由于气相和燃料介质中波的传播特性不同， 波的衰减性质略有区别.

图12 给出了在31 μs 时刻， 波阵面位置及中轴线上流场温度、燃料颗粒温度的分布. 可以看到， 31 μs 时刻冲击波在云团内部传播(冲击波在平面中的位置用蓝色线标记)， 流场温度、颗粒温度增加， 因此颗粒的蒸发动力学效应增强.

图12 t=31 μs 时刻在y =0 m 轴线上气相场温度以及燃料蒸气浓度的分布Fig.12 Distribution of gas field temperature and fuel vapor concentration along y =0 m axis at 31 μs

2.1.2 液体燃料的运动

燃料的运动取决于流场的气体动力学与燃料云团的相互作用. 在321 μs 时间内， 燃料云团在x， y 方向上获得的平均速度和加速度演化如图13 所示. 可以看到， 在爆炸驱动最初期，冲击波作用使得燃料获得较大的径向、垂向加速度; 随后， 燃料在反射波和气液界面湍流涡的影响下， 径向加速度急剧下降， 运动速度减慢; 在200 μs 时刻后， 由于云团内侧附近流场压力的升高， 推动燃料介质继续沿径向运动， 加速度又逐渐上升. 图 13 中， 在 112 μs 时刻， x 方向上的速度几乎达到最大值22.8 m/s; 而y 方向上的速度在295 μs 时刻达到最大值8.97 m/s.

图13 燃料云团的平均速度和加速度Fig.13 Average velocity and acceleration of fuel cloud

分析计算所得的平均速度和加速度曲线可知， 在气动力作用下， 燃料介质加速度的变化存在急剧上升和下降2 个阶段.

式中: ax1为常数; axmax为加速度的最大值; b 为衰减系数. 在气动力作用下， y 方向燃料的加速度也满足上述函数关系.

分析和整理本工作的数值计算结果， 给出加速度的指数函数关系式如下.

可见， 燃料获得 x， y 方向加速度的增长系数分别为 49.92 和 8.26; 随后， 燃料在 x， y 方向加速度的衰减系数降至1.07 和1.97. 增长和衰减系数反映出气动力和大、小尺度涡对燃料运动的作用效应， 这也使得燃料介质呈现复杂的运动性质.

2.2 冲击波诱导涡与燃料界面的相互作用

在爆炸驱动的最初期阶段， 向前传播的冲击波作用于燃料云团， 界面附近产生Richtmyer-Meshkov 不稳定性(RMI)， 流场诱发出尺寸不一的旋涡， 加速了燃料的混合和扩散.

2.2.1 气液界面旋涡与液体燃料的相互作用

图14 给出了颗粒初始粒径为60 μm 的燃料云团气液界面流场涡量随时间的演化. 由图可见， 由于剪切效应占主导， 燃料界面的内侧外缘附近流场产生了大尺寸的旋涡结构. 在40 μs 时刻， 最大涡心强度达到 1.42×105s-1， 极大地促进了颗粒的分散.

此外， 本工作研究结果表明， 当冲击波驱动液体燃料时， 由于密度和压力梯度导致斜压效应， 气液界面诱导出涡量， 界面产生Richtmyer-Meshkov 不稳定性(RMI)， 并随时间扰动逐渐发展. 图15 给出了270 μs 时刻， 粒径为60 μm 的云团界面附近流场出现的尺寸不一的涡串结构. 在随后湍流涡的作用下， 扰动进一步发展， 燃料云团与周围流场发生混合.

2.2.2 诱导涡的发展和演化

下面将针对颗粒初始粒径为60 μm 的燃料云团， 波诱导出涡对随时间的演化过程进行分析.

(1) 扰动沿y 方向的发展和演化.

图16 给出了4 个时刻涡量在界面流场中的分布. 可以看到， 随时间的发展， 涡的强度逐渐降低， 即 10 和 17 μs 的涡量大约为 2 000 和 1 600 s-1， 大于 82 和 91 μs 时的涡量 1 200 和800 s-1.

图 14 60 μm 云团流场涡量随时间的演化Fig.14 Vorticity evolution of flow field over time for 60 μm radius cloud

图15 270 μs 时刻流场涡量和颗粒的分散状态Fig.15 Vorticity of flow field and dispersion of particles at 270 μs

图16 不同时刻云团内侧界面流场涡量Fig.16 Flow vorticity in the inner interface of cloud for different time

图 17 给出了 173 μs 时刻， 分别在 x=0.021 m(side1)和 x=0.024 m(side2)分界线上涡对的分布(分界线位置见图8)， 可以清晰地观察到界面附近涡对的结构. 此外， 在涡对的卷吸作用下， 液体燃料的蒸发动力学效应增强. 在y =0 m 的直线上， x=0.020 8 ～0.021 6 m 之间的区域燃料蒸汽体积的体积分数最大， 约为0.526.

图17 173 μs 时刻云团内侧界面流场涡量Fig.17 Flow vorticity in the inner interface of cloud at 173 μs

(2)扰动沿主流方向的发展和演化.

扰动在燃料介质中传播， 涡的强度逐渐降低. 在云团的外侧界面， 由于存在压力梯度和密度梯度， 界面同样出现类似的不稳定现象， 涡对结构不断演化. 图18 给出了y =0 m 轴线上流场涡沿x 轴的演化. 顺时针和逆时针旋转的涡分布在气体和液体燃料内侧、外侧界面.

2.2.3 诱导涡对液体燃料的影响

图19 给出了颗粒初始粒径分别为60 和80 μm 的云团内侧界面流场涡量的演化情况. 可以看到， 对于60 μm 颗粒组成的云团界面附近流场诱导的旋涡， 涡运动的方向为顺时针， 涡量为负值; 而对于80 μm 颗粒组成的云团界面附近流场诱导的旋涡， 涡运动的方向为逆时针， 流场性质略有不同. 颗粒分布如图20 所示.

旋涡运动沿着液体燃料界面伸展， 随着不稳定性的发展， 气液混合增长速率随小尺度涡结构的发展得到增强. 因此， 扰动产生的大尺度旋涡和小尺度旋涡使界面不断变形， 加剧了颗粒的运动和分散.

图18 y =0 m 轴线上涡量的分布Fig.18 Vorticity distribution along the axis of y =0 m

图19 y =0 m 轴线上流场涡量随时间的演化Fig.19 Flow vorticity over time along the axis of y =0 m

图 20 400 μs 时刻颗粒初始粒径分别为60， 80 μm 的燃料云团状态Fig.20 Fuel cloud with initial particle radius 60， 80 μm at 400 μs

3 结 论

爆炸驱动燃料分散过程存在冲击波、湍流、颗粒动力学等复杂物理场的相互作用， 本工作开展了爆炸作用冲击波与燃料界面相互作用过程的数值模拟， 重点分析最初期波、诱导湍流涡对燃料颗粒分散的作用机制.

(1) 爆炸驱动液体燃料分散是云团在空间展开的、逐渐形成燃料和空气混合的可爆性云团过程. 在爆炸作用力下， 云团尺寸发展较快; 当云团内的压力接近外界环境压力时， 云团运动变慢， 云团尺寸变化趋缓. 计算所得的云团尺寸随时间按照指数律发展， 满足云雾膨胀的相似率: 半径的时间演化， R = 5.026 46-5.137 63e-0.079t; 高度的时间演化，H =1.774 59-1.853 92e-0.084t(t:0 ～60 ms). 云团半径增长系数(0.079)和云团高度增长系数(0.084)反映了中心装药释放的总能量与爆炸驱动力密切相关.

(2) 在爆炸驱动最初期， 波在两相介质中传播时发生衰减， 气相和燃料介质中超压峰值的变化满足幂数律的关系. 此外， 由于受到空气和燃料介质的阻抗作用， 云团中波传播速度较慢，因此相对于空气介质， 燃料介质中波阵面发生弯曲. 同时， 波后云团状态发生变化， 燃料的运动速度提高. 如在112 μs 时间内， 初始粒径为60 μm 的燃料云团在x 方向上的运动速度增加至 22.8 m/s， 最大加速度为 5.0×105m/s2.

(3) 波绕过燃料云团的外缘产生大尺度的旋涡. 对于颗粒初始粒径为60 μm 的燃料云团，大尺度涡的强度最大达到1.42×105s-1， 驱动着颗粒向外分散. 同时， 气液界面附近流场密度和压力梯度诱导出涡量， 界面产生不稳定性现象， 并随时间逐渐发展成沿界面排列的多个涡对结构， 流场湍流得到发展.

因此， 在爆炸驱动最初阶段， 流场湍流度的增加叠加绕射波的双重作用， 加剧了气液的混合， 也为燃料颗粒的湍流分散提供了较大的初始动能.