激波与液滴的相互作用问题，是一个与Richtmyer-Meshkov(RM)不稳定性有关的、典型的可压缩性气液两相流问题^[1]，对该问题的深入研究在超音速雨滴侵蚀、燃气轮机燃烧室的设计、激波抛洒以及爆轰发动机中燃料燃烧的稳定性等方面均有着重要应用。

W.Lane^[2]最早在激波管中对液滴在稳态和瞬态气流中的破碎进行了实验研究，并提出液滴的变形破碎是一个边界层不断脱落的过程。W.E.Krauss^[3]随后也研究了边界层分离与液滴破碎雾化之间的关系。O.Engel^[4]着重研究了激波对液滴的雾化作用，他认为高速气流在液滴表面生成不稳定波在液滴破碎过程中起重要作用。A.Ranger等^[5]在水平激波管中进行了激波马赫数1.3~3.5的实验，揭示了液滴在激波作用下出现压缩变形到破碎雾化的过程。发现入射激波马赫数较高时，液滴的压缩变形和破碎雾化几乎同时发生，液滴迅速形成薄雾状结构。C.Kauffman等^[6]研究了燃料液滴在氧气环境中与激波相互作用的过程：在入射激波马赫数3~5的情况下，液滴的燃烧先发生在从原始液滴脱落的雾状燃料上，在激波扫过一段时间后，整个燃料的燃烧过程才会发生。A.Wierzba等^[7]利用全息干涉技术研究了液滴在激波作用下的变形破碎过程，并把液滴脱落破碎分成4个阶段：压缩变形，横向直径增大至最大值，破碎雾化和横向直径减小至零。L.P.Hisang等^[8-9]给出了在液滴破碎模式与We(Weber数)的关系：随着We从零开始增加，液滴会依次出现无变形、无振荡变形、振荡变形、袋式破碎、混合破碎和剪切破碎模式，并当We继续增加时会出现突发破碎(catastrophic breakup)情况。

陆守香等^[10]研究了激波冲击后液滴变形和破碎的近似理论模型。耿继辉等^[11]实验观测了液滴在激波作用下的加速、变形及破碎过程。蔡斌等^[12]通过应用VOF方法和湍流模型对液滴在气流中的破碎过程进行了数值模拟，并分析了We、Oh(Ohnesorge数)以及气液密度比对液滴破碎过程的影响。姚雯等^[13]通过求解二维Euler方程，并将表面张力转换为体积力，考察了激波流场中液滴变形与破碎过程。楼建峰等^[14]将VOF方法和标准k-ε湍流模型组合，进行了液滴在气流中变形破碎过程的数值模拟，并分析了几个关键参数(Weber数、Ohnesorge数、气液密度比)对液滴破碎过程的影响。其结果表明，We对液滴变形破碎起促进作用，而Oh和液气密度比起阻碍抑制作用。

本文中，首先在水平激波管中进行激波与液滴作用的实验研究，然后利用数值模拟技术计算该实验条件下液滴破碎过程。再通过对比分析，详细描述液滴在激波作用下的变形破碎过程。

1.   实验方法

水平激波装置示意图如图 1所示。包括真空箱的激波管整体总长为12 m。激波发生系统由高压气瓶、圆形高压段、圆形低压段及方形低压段组成，实验段通过法兰与方形低压段及真空箱连接，两个压力传感器安装在方形低压段，用于测量激波马赫数。高速摄影仪放置在实验段的侧面，在实验段的前后侧板上装有透明的有机玻璃板，上、下板开有液滴下落入口和排液孔。实验段内截面尺寸为120 mm×120 mm，长为500 mm。在实验段上方的液滴下落入口插入不同直径的尖针，以便形成一个在实验段自由下落的直径可控的液滴。控制破膜时间，使入射激波刚好扫过下降的液滴，用高速摄影仪拍摄整个液滴的变形破碎过程。

图  1  激波管装置示意图

Figure  1.  Sketch of shock tube

下载: 全尺寸图片 幻灯片

2.   数学模型与计算方法

针对液滴在激波作用下的变形破碎，为了方便地模拟两相界面的变形情况，数值模拟中选用基于VOF模型的RANS方程，湍流模型使用标准k-ε双方程模型。其二维控制方程组为如下通用形式：

式中：ρ为两相平均密度，ρ=α₁ρ₁ + (1-α₁)ρ₂，其中α₁、ρ₁和ρ₂分别表示气相质量分数、气相密度和液相密度；Φ为通用因变量，其取值为1、u(x方向速度)、v(y方向速度)、E(总能)、k(湍流动能)和ε(湍流耗散)时方程分别对应连续方程、动量方程、能量方程、湍流动能输运方程和湍流耗散方程；Γ_Φ 为输运系数，包括黏性输运(根据气相质量分数加权计算^[15])和湍流输运；S_Φ为源项，主要为重力项、压力作用项、湍流作用项和相间相互作用项等组成。方程中各项具体表达式及参数取值可参见文献[15-16]。

根据液滴在气流中破碎的结构特点，选取了如图 2所示的计算区域。其中计算区域尺寸为2 cm×8 cm，液滴直径根据实验数据分别为1.64、2.16、2.46 mm。计算使用ANSYS Fluent软件中的Coupled求解器，数值离散使用二阶MUSCL格式，计算区域网格数目为42万。波后区域参数根据运动正激波解求得，也以此作为左边界入口边界参数，其他3个边界选自由出口边界。设定激波初始位置位于液滴前缘上游2 mm处。在与实验结果对比时，第1张图片的时刻均为激波通过后的0.125 ms。

图  2  计算区域示意图

Figure  2.  Schematic diagram of calculating range

下载: 全尺寸图片 幻灯片

3.   结果分析

实验采用We(Weber数)与Oh(Ohnesorge数)来分析液滴破碎变形，其具体计算公式分别为：We=ρ_au²d₀/σ_w和$O h=\mu_{\mathit{\boldsymbol{w}}} / \sqrt{\rho_{\mathit{\boldsymbol{w}}} d_{0} \sigma_{\mathit{\boldsymbol{w}}}}$，We表示惯性了力和表面张力之比，Oh表示黏性力和表面张力之比。其中，ρ_a为空气密度，u为波后气流速度；ρ_w、σ_w和μ_w分别为液滴的密度、表面张力和动力黏度；d₀为液滴的初始直径。实验和计算中各参数的具体数值，如表 1所示，表中Ma为入射激波马赫数。

表  1  实验参数

Table  1.  Experimental parameters

实验	液滴	Ma	d0/mm	We	Oh
1	水	1.10	1.64	94.48	0.002 9
2	水	1.10	2.16	124.44	0.002 6
3	水	1.10	2.46	141.72	0.002 4
4	水	1.25	1.68	518.75	0.002 9
5	甘油	1.10	2.46	167.23	3.450 0

下载: 导出CSV 

| 显示表格

实验1的实验与数值模拟结果如图 3(a)所示。左列为数值模拟结果，图示的是密度云图，右列为实验照片，上下两幅图时间间隔为0.25 ms。从图中可以看到，实验与计算的液滴变化趋势是相近的，从开始的压缩变形为一个薄圆盘状结构，其后液滴的横向直径增大，轴向直径降低；紧接着，薄圆盘状结构逐渐扩展，更细小的液滴从上下两端不断脱落，原始的球形液滴演化成了云团状，并继续沿横向和轴向扩张。由前4幅对比可以看出，激波作用后刚开始数值模拟中液滴的横向变形没有实验中的剧烈，仍显得比较“胖”，第4幅左图数值模拟结果也就没能出现类似右图中的细棍状结构。最后两张对比图片中，虽然大致上两者的雾化扩张类似，但是，数值模拟中液滴中间已经只剩下一个点状水粒，而实验照片显示液滴中间是被打散的极小液滴密集而成的云团结构。最后一张对比图中，可以明显看到实验中液滴已经雾化成云团，而数值模拟结果仍然显示液滴有未破碎的内核存在，这是因为VOF模型中没有加入液滴的破碎机制。

图  3  实验1~3的数值模拟结果(左)与实验照片(右)(Δt=0.25 ms)

Figure  3.  Numerical simulation (left) and experimental (right) results of experiments 1-3

下载: 全尺寸图片 幻灯片

图 3(b)给出了实验2的数值模拟结果和实验结果。由图可以明显看出，数值模拟的液滴变形演变速度比实验的慢。与图 3(a)不同，在图 3(b)的最后一幅图中，两相界面上出现了尖钉和气泡，这是RM不稳定性作用的结果。因为当初始液滴的直径增大时，由于液滴的体积和质量更大，其破碎变形肯定比小直径的液滴慢，所以在图 3(a)中液滴还没来得及出现RM不稳定的结构就破碎了。

图 3(c)给出了实验3的数值结果和实验结果。在前5幅图中，液滴的变化和实验照片一致。在第6幅图中，计算与实验结果都出现了RM不稳定性引起的尖钉和气泡结构，计算结果中向左侧突起的尖钉结构比较明显，而实验结果的右侧包裹有明显的气泡结构，可是计算与实验中液滴的横向直径增加趋势仍然是一致的。实验结果的前4幅图中，液滴中偏左位置的白色圆孔不是气泡结构，而是由于液滴反光造成的拍摄干扰。

图 4给出了液滴横向直径与初始直径之比d/d₀随时间的变化关系。在激波作用后液滴变形的初始线性增长阶段和中期非线性发展阶段，数值结果与实验结果符合很好，而在后期的破碎发展阶段，数值结果与实验结果偏离。尤其是实验1中最后一组数据，实验观测到液滴核心的横向直径已经因液滴周围突发破碎而减小，而数值结果则相反，原因是数值计算中没有加入突发破碎模型。

图  4  实验1~3液滴横向直径与初始直径之比

Figure  4.  Ratio of transverse diameter andinitial diameter of experiments 1-3

下载: 全尺寸图片 幻灯片

图 5给出了实验4的数值计算和实验结果对比图。实验4与实验1具有相同的Oh，而两者的We差别较大。从实验中可以看出，在激波扫过后，液滴几乎没有出现明显的压缩变形就发生突然破碎，这说明随着We增大到一定程度，液滴的破碎模式发生了根本性改变：We较小时液滴经历由变形到破碎的剪切破碎模式，而当We增大到一定值时，为突发破碎模式。一般而言，破碎模式的临界We在500左右，本文中实验也验证了这一点。图 5中的计算结果没有出现液滴突发破碎模式，说明在数值模拟中应该加入突发破碎的机制，而其他研究者还没有指出这一点。

图  5  实验4的数值模拟结果(左)与实验照片(右)(Δt=0.25 ms)

Figure  5.  Numerical simulation (left) andexperimental (right) results of experiment 4

下载: 全尺寸图片 幻灯片

综合上面的讨论可知，当Oh相同时，We越大，破碎速度越快。实验1的液滴经过1.5 ms后，横向直径才达到最大而开始裂解破碎。而实验4的液滴经过0.25 ms，就完成了横向扩张，开始裂解破碎。这表明We在液滴破碎过程中加强了液滴的不稳定，为不稳定因素。

为了得到相同We条件下，不同Oh对液滴变形和破碎的影响，对实验5进行了模拟计算，其We与实验3相近，而Oh则相差较大，结果如图 6所示。甘油液滴从初始变形上就不同于水滴，甘油液滴在波后气流的作用下，先是迎风面开始扁平化，背风面仍然保持原先的半圆状结构，同时两侧也在伸长，但后部的弧形结构明显很完整。相对而言，水滴受到波后气流作用下，其右侧半圆结构很快向内凹陷，形成开口状结构。这两个明显的差异，应该是由于甘油的黏性系数远大于水的黏性系数而造成。相同初始直径的水滴，在实验3的模拟结果中，刚到第4张图就可以观察到，水滴的微滴剥离非常明显，并且在后续图片中，已经出现了表示RM不稳定性的尖钉状结构。而在图 6中，随着甘油液滴横向扩张结束，两翼开始向后弯曲变形，但整体结构仍然保持完整，没有出现大规模的微滴剥离。并且在数值模拟结果中，整体液滴的变形演化速度慢于水滴。由此可见，Oh的减小使破碎速度显著加快，说明Oh对破碎起抑制作用。

图  6  实验5的数值模拟结果(Δt=0.25 ms)

Figure  6.  Numerical simulation results of experiment 5

下载: 全尺寸图片 幻灯片

图 7为实验3(水滴)与实验5(甘油液滴)的数值计算结果的涡量云图对比图。由图可以看出，甘油液滴算例中的涡量从开始就处于不对称的状态，之后出现了类似于卡门涡街的现象，不断有涡旋一上一下地从液滴后方脱落，向下游传播，一直持续到液滴大规模破碎为止。水滴的算例中则没有这个现象。这是因为甘油的黏性远大于水滴的黏性，其破碎速度比水滴的慢很多。

图  7  实验5和3的数值模拟涡量图(Δt=0.25 ms)

Figure  7.  Vorticity cloud of experiments 5 and 3

下载: 全尺寸图片 幻灯片

4.   结论

实验观测和数值计算结果表明：(1)液滴在激波的作用下，要经历从压缩变形、RM不稳定性变形、细小液雾剥离到雾化破碎等过程。不同参数下的结果显示，液滴变形破碎的发展趋势一致，但发展速度明显不同。数值模拟结果与实验结果基本吻合。(2)实验数据与数值模拟结果显示，在液滴破碎过程中We的增加促进液滴的破碎，为不稳定因素。当We增加到一定值后，液滴破碎模式由剪切破碎转变为突发破碎。Oh的增加对液滴的破碎起抑制作用，黏性同样会抑制液滴的破碎。

对比分析指出，在液滴变形的初期及中期，数值计算与实验结果吻合较好。在液滴变形的后期，需要在数值计算中考虑雾化模型；而随着激波马赫数的提高，还要考虑突发破碎机制，其中包括激波冲击气/液界面引起的扰动。根据计算的速度云图^[17]，可知液滴的上下承受着高剪应力(由速度梯度造成)，所以液滴的变形破碎不仅与RM不稳定性有关，还与Kelvin-Helmholtz(KH)不稳定性有关。图 6~7所示的对甘油液滴的计算结果，对于寻找RM不稳定性的实验方法有重要参考价值，即如果使用黏弹性流体(如果冻)，那就应该考虑增加激波强度，以保证其本构方程与Navier-Stokes方程相似。