气流速度振荡场中幂律液膜不稳定性分析

姚慕伟，贾伯琦，杨立军，富庆飞

北京航空航天大学 宇航学院，北京 102206

在航空航天推进系统中，燃油的喷注雾化过程常简化为平面液膜在气流中的稳定性问题。鉴于平面液膜的初始稳定性影响着最终破碎形成的液滴尺寸及分布，因此研究平面液膜在气流中的稳定性不仅对工程问题具有较大的应用价值，且作为经典的流体力学问题，也得到了国内外学者的密切关注。Squire[1]最早研究了无黏平面液膜在静止气体中的稳定性问题，提出液膜存在两种模式的不稳定，分别称为正弦模式和曲张模式。基于Squire的研究，学者们对平面液膜的稳定性问题做了进一步研究。Dombrowski[2]、Lin[3]、Li[4-5]考虑了液膜黏性的影响，研究发现表面张力使扰动增长率衰减，而气液相对速度促进液膜的初始不稳定。Barreras[6]、Dombrowski[7]、Ibrahim[8]等则分别考虑了气体黏性、密度和压缩性对液膜稳定性的影响。在火箭发动机燃烧室中，发生燃烧不稳定时的压力振荡可处理成气流速度振荡[9]，此时平面液膜与气流间的动量交换过程受振荡气流的影响，使液膜的不稳定性增加[10-11]。此时流动系统的不稳定区域可被分为Kelvin-Helmholtz(K-H)和参数不稳定区域，气流速度振荡频率的增加促进了K-H不稳定区域的不稳定性，增强了参数不稳定区域内的稳定性[12]。Sivadas[10,13]和Mulmule[14-16]等的实验也显示平面液膜在外加振荡声场的环境下更容易破碎，且破碎产生的液滴尺寸更小。以上的研究对象都是牛顿流体，最近，Jia等[17]研究了黏弹性平面液膜在气流速度振荡中的稳定性；结果显示在K-H不稳定区域扰动以行波的形式增长，在参数不稳定区域扰动的增长呈现为驻波的形式，且黏弹性液膜在所有不稳定区域中呈现出比牛顿液膜更稳定的特性。本文针对凝胶推进剂表现出的幂律性质流体所形成的平面液膜，考虑其在气流速度振荡情况下的稳定性问题，使用与Jia等[17]、Kumar和Tuckerman[18]类似的Floquet理论处理参数共振系统研究幂律性质、振荡幅值及振荡频率对液膜稳定性的影响，以期丰富平面液膜参数不稳定的理论研究。

1 物理模型及数学推导

采用的液膜物理模型如图1所示，其中，气相存在沿流向振荡的速度，且气体无黏不可压，重力及其他有势力均被忽略。建立静止坐标系，x方向为沿流向，y方向为垂直液膜流动的方向。为区分物理量，下标i=l, g分别表示液相和气相。设液膜的厚度为2a，ρl为液体的密度，σ为表面张力，pl为液体的扰动压力，τ为液体的偏应力张量，ul和vl分别为x和y方向液体的小扰动速度，Ul为液膜的基本流速度，Ug为气流速度，U0为气流定常速度，ΔU为振荡速度，ω为振荡频率，ρg为气体密度，pg为气体的扰动压力，η1和η2分别为受扰动的液膜表面不稳定波的振幅。液相的控制方程为

(1)

(2)

(3)

式中：τxx、τxy、τyx和τyy为偏应力张量的分量，其下标表示对应的方向；t为时间。

气相速度的表达式为

图1 平面幂律液膜的物理模型示意图Fig.1 Schematic of physical model of power-law liquid sheet

Ug=U0+ΔUcos (ωt)

(4)

幂律流体的本构方程为

(5)

考虑到液膜沿x方向流动，只保留x方向的法向偏应力张量τxx，并考虑线性情况，忽略非线性项，则动量方程化简为

(6)

(7)

引入速度影响因子g[19-20]，表示速度随位移的变化，单位为1/s：

(8)

将式(8)代入式(6)中，得到：

(9)

因气相假设为无黏，故气相的控制方程可写成势函数的形式：

(10)

(11)

式中：φ为气体速度势函数；下标g1表示液膜下方的气体；下标g2表示液膜上方的气体。

在考虑到气流速度振荡的条件下，不能直接将变量展开成正则模的形式。这里采用与Kumar和Tuckerman[18]类似的方法，将扰动量展开成Floquet的形式：

(ul,vl,pl,φg1,φg2,pg)=exp(βt+ikx)×

(12)

(13)

将式(12)和式(13)代入式(1)、式(7)和式(9)中，得到：

(14)

(15)

(16)

(17)

求解方程式(17)，得到液相压力的表达式：

(18)

式中：B1s和B2s均为积分常数。

同样地，将式(12)和式(13)代入式(10)和式(11) 中，并考虑气相扰动应在无穷远处为0，即

φg1→0y→-∞

(19)

φg2→0y→+∞

(20)

得到：

(21)

(22)

式中：A1s和A2s均为积分常数。要确定这些积分常数，还需要气液界面上的边界条件。在气液界面上的线性化运动边界条件为

(23)

(24)

(25)

(26)

同样地，(η1,η2)可写为

(27)

(28)

在气液界面上，应满足应力平衡条件：

(29)

而气相的压力可表达成势函数的形式：

(30)

将式(12)、式(13)、式(18)、式(19)、式(21)、式(28)和式(30)代入式(29)，可建立色散关系：

Dsηs+Es-1ηs-1+Gs+1ηs+1+Fηs-2+Fηs+2=0

(31)

色散关系可被表示成矩阵的形式：

(32)

(33)

式中：

(34)

2 结果与讨论

2.1 结果验证与阶数确定

图2 无气流速度条件下与文献[19]结果的比较Fig.2 Comparison with results obtained by Ref.[19] under condition of no gas velocity

图3 不同节点数下的增长率Fig.3 Growth rate vs node number increase

在正式讨论分析之前，还需确定色散矩阵的阶数N。理论上色散矩阵为无穷阶时计算的结果最为精确，但实际不可能取无穷阶，因此需要在误差允许的范围之内取较小的节点数s。图3做出了在不同节点数s情况下，两个不稳定区域的增长率及最大增长率Br max。从图3(a)中可以看出，在气流速度振荡条件下，不稳定区域有多个，其中波数范围最小的不稳定区域为K-H不稳定区域，随着波数的增大出现的区域为气流速度振荡引起的参数不稳定区域。在不同节点数下，增长率曲线几乎重合；借助图3(b)观察到当s=0时，只有一个K-H不稳定区域；随着节点数的增加，当s>3以后，各不稳定区域内的最大增长率数值基本不再改变。因此，从计算精度和经济性综合考虑，取s=4，因此色散矩阵的阶数为N=2s+1=9。

2.2 振荡幅值变化下幂律性质对液膜稳定性的影响

图4反映了表观雷诺数逐渐增大时无量纲振幅为0和0.56时增长率的变化情况。当ε=0时，只存在K-H不稳定区域，因为此时液膜并没有受到振荡的激励，只发生由速度剪切引起的K-H不稳定。随着ε的增加，不仅出现了多个参数不稳定区域，且K-H不稳定区域的最大增长率、主导波数及截止波数都有明显增加。这些规律与Jia等[17]在对平面液膜研究所发现的结论都是一致的，因此主要观察幂律模型的参数对增长率的影响。

测试的振荡幅值为[0,1]之间10个均匀分布的值。图5(a)统计了在不同振荡幅值下K-H不稳定区域的最大增长率与Ren的关系，发现随着Ren的增加，最大增长率逐渐增加。物理上，表观雷诺数的增加反映为稠度系数的减小，代表液膜的黏性减小，黏性耗散减小，液膜更不稳定，因此增长率更大。观察到随着ε的增加，最大增长率

图5 不同振荡幅值下表观雷诺数对K-H不稳定区域的影响Fig.5 Effects of generalized Reynolds number on K-H instability region with different oscillation amplitudes

逐渐增加，且当ε较小时增加的程度小于ε较大时增加的程度。图5(b)统计了K-H不稳定区域最大增长率对应的主导波数mdom的变化关系，发现K-H不稳定区域的主导波数在Ren的变化范围内保持不变。

图6描述了不同振荡幅值下表观雷诺数的变化对参数不稳定区域增长率的影响，其中图6(a)显示了最大增长率的变化。与K-H不稳定区域最大增长率随Ren的增加而单调增加所不同的是，参数不稳定区域的最大增长率随Ren的增加先减小后增加，且最大增长率发生转折点对应的Ren不随振荡幅值的变化而改变。图6(b)显示主导波数的变化趋势为随着Ren的增加，参数不稳定区域主导波数mdom逐渐减小。

Ω=0.125,G=5×10-6图6 不同振荡幅值下表观雷诺数对参数不稳定区域的影响Fig.6 Effects of generalized Reynolds number on parametric instability region with different oscillation amplitudes

图7展示了不同振荡幅值下无量纲速度因子G的变化对K-H不稳定区域增长率的影响。其中，G的取值为[10-7，10-4]内均匀分布的10个值。观察到随着振荡幅值的增加，最大增长率逐渐增加；虽然振幅的增加是均匀的，但最大增长率的变化却是逐渐加剧的，且主导波数的增加同样是不均匀的。随着无量纲速度因子G的增加，K-H不稳定区域内最大增长率逐渐增加，这与文献[19] 中的描述是一致的。

图8描述了不同振荡幅值下无量纲速度因子变化对参数不稳定区域增长率的影响。与幂律指数的影响规律相似，随着G的增加，在参数不稳定区域最大增长率先减小后增加，且增长率转折点所对应的G不随振荡幅值的变化而改变。对于主导波数，随着G的增加整体来看逐渐减小，表示参数不稳定区间随着G的增加逐渐向波数较小的区域移动。

图9描述了在不同振荡幅值下幂律指数变化对K-H区域增长率的影响特征，发现随着n的增大，最大增长率单调增加，同时主导波数不变。

n对参数不稳定区域的影响如图10所示，随着n的逐渐增加，最大增长率呈现先减小后增加的趋势，且振荡幅值的增加不影响最大增长率发生转折点所对应的n值，观察到此时对应的n=5。推测是因为此时n的取值最偏离n→0代表的无黏流体和n→1代表的牛顿流体之间，因而流体所呈现的非牛顿特性最强，所以液膜最稳定，增长率在各n取值中最小。与表观雷诺数和无量纲速度因子的影响所不同的是，当ε∈[0.56,0.78]时，参数不稳定区域内的主导波数随着幂律指数的增加先减小后增大，而当无量纲振荡幅值处于该区间之外时，主导波数逐渐减小至一特定值后保持不变。

图7 不同振荡幅值下无量纲速度因子对K-H不稳定区域的影响Fig.7 Effects of dimensionless velocity factor on K-H instability region with different oscillation amplitudes

图8 不同振荡幅值下无量纲速度因子对参数不稳定区域的影响Fig.8 Effects of dimensionless velocity factor on parametric instability region with different oscillation amplitudes

Ω=0.125,G=5×10-6图9 不同振荡幅值下幂律指数对K-H不稳定区域的影响Fig.9 Effects of power-law index on K-H instability region with different oscillation amplitudes

Ω=0.125,G=5×10-6图10 不同振荡幅值下幂律指数对参数不稳定区域的影响Fig.10 Effects of power-law index on parametric instability region with different oscillation amplitudes

2.3 振荡频率变化下幂律性质对液膜稳定性的影响

在对无量纲振荡频率变化的影响研究中，选取Ω=0.050～1.000，并采样10个点进行研究。通过图11(a)和图11(b)可以看出，振荡频率较小时，液膜的不稳定区间较多，在Ω=0.050的条件下，有4个不稳定区域；随着振荡频率的增大，参数不稳定区间的位置逐渐移向高波数区域，最终只存在一个K-H不稳定区域。随着振荡频率的增加，K-H不稳定区域的最大增长率、主导波数和截止波数都有一定程度的增加，这与Jia等[17]对平面液膜的观察结果是一致的。

图12统计了不同振荡频率下幂律指数变化对K-H区域增长率的影响特征，发现随着n的增大，最大增长率逐渐单调增加，同时主导波数不变。

图13描绘了不同振荡频率下幂律指数的变化对参数不稳定区域增长率的影响特征。当振荡频率较小时，最大增长率并没有出现先减小后增大的特征，而是单调增加；当Ω进一步增大，最大增长率先减小后增加，且随着Ω的增加，增长率转折点所对应的幂律指数增大。但当Ω>0.261时，在n=0.5～0.7的范围内最大增长率此时已经为负数，表现为稳定的特征，此时最大增长率的取值并没有物理意义,因此图13(a)中略去了纵坐标为负的区域。观察到当Ω>0.050、n>1时最大增长率会发生逆转：当液膜为n<1的剪切变稀流体时，振荡频率越小，最大增长率越大；当液膜为n>1的剪切变稠流体时，振荡频率越大，最大增长率越大。在图13(b)中可直观地看出当振荡频率增加时，只有在大幂律指数下存在参数不稳定区域，且随着频率的增加，发生参数不稳定时的n逐渐增加。在整个无量纲振荡频率变化范围内，主导波数随着幂律指数的增加逐渐减小或保持不变。

图11 不同振荡频率下幂律指数对液膜失稳特性的影响Fig.11 Effects of power-law index on instability characteristics of liquid sheet with different oscillation frequencies

图12 不同振荡频率下幂律指数对K-H不稳定区域的影响Fig.12 Effects of power-law index on K-H instability region with different oscillation frequencies

图13 不同振荡频率下幂律指数对参数不稳定区域的影响Fig.13 Effects of power-law index on parametric instability region with different oscillation frequencies

图14(a)和图14(b)分别讨论了当振荡频率变化时无量纲速度因子G对K-H不稳定区的影响。其中，图14(a)表示最大增长率随G的变化，图14(b)表示主导波数随G的变化。如图14(a)所示，流体的G逐渐增大会导致增长率逐渐增大，这与不施加声场扰动时液膜的增长率特征是一致的；虽然最大增长率逐渐增加，但最大增长率对应的主导波数却保持不变。

图15统计了不同振荡频率下无量纲速度因子的变化对参数不稳定区域增长率的影响，其中图15(a)表示最大增长率的变化规律。与幂律指数的影响规律相似，在Ω较小时，最大增长率呈现出单调增加的趋势；当Ω进一步增大，最大增长率随着G的增加先减小后增加。从图14中已经得知G的增加会使液膜更不稳定，物理上，如果假设G为速度通过单位长度的改变量，G越大则表示在单位长度内外界对液膜的作用越大，外界对液膜输入的能量越大，则液膜越不稳定。注意到在参数不稳定区域内振荡频率的增加使最大增长率逐渐减小，这与K-H不稳定区域内振荡频率作用的表现规律是截然不同的。当Ω很大时，参数不稳定区间将不再存在。从图15(b)主导波数随G的变化曲线中得知，在Ω=0.050时主导波数保持不变；Ω增大时，主导波数随着G的增加而减小。但当Ω=0.367时，介于[2×10-7，1×10-4]之间的G只存在K-H不稳定区间，参数不稳定区间因耗散而衰减为负值。

图14 不同振荡频率下无量纲速度因子对K-H不稳定区域的影响Fig.14 Effects of dimensionless velocity factor on K-H instability region with different oscillation frequencies

图15 不同振荡频率下无量纲速度因子对参数不稳定区域的影响Fig.15 Effects of dimensionless velocity factor on parametric instability region with different oscillation frequencies

图16显示了K-H不稳定区域最大增长率和主导波数随Ren的变化规律。图16(a)观察到随着振荡频率的增加，最大增长率单调增加。虽然振荡频率的变化间距是均匀的，但K-H不稳定区域的最大增长率却在小振荡频率下变化较敏感，当Ω>0.261以后最大增长率无明显增加。图16(b) 则是K-H不稳定区域Ren对主导波数的影响。与最大增长率类似的是，主导波数的变化也呈现出随着振荡频率的增加而增大的现象，但主导波数不随Ren的增加而改变。

图17(a)和图17(b)讨论了不同振荡频率下参数不稳定区域最大增长率和主导波数随Ren的变化情况。因为随着Ω的增加，不稳定区域逐渐减少，因此这里只讨论第一个参数不稳定区域的性质。从图17(a)中可发现，随着振荡频率的增加，最大增长率的变化程度也愈发剧烈；此时，最大增长率的转折点所对应的Ren也随着振荡频率的增加逐渐增大。随着Ω的增加，不稳定区域的个数逐渐减少，当Ω>0.578时，只存在K-H不稳定区域。这与Dighe和Gadgil[21]提到的当液膜的受迫频率超过截止频率时液膜表面的无量纲振幅不再增加一致。而类似地，无量纲振荡频率Ω=0.578时也可以视为“截止频率”，当激励频率高于该受迫频率，液膜K-H不稳定区域的最大增长率和截止波数都不再有显著增加。图17(b)则显示出主导波数在Ω=0.050下保持不变；随着Ω的增加，主导波数开始随着Ren的增加而逐渐

图16 不同振荡频率下表观雷诺数对K-H不稳定区域的影响Fig.16 Effects of generalized Reynolds number on K-H instability region with different oscillation frequencies

图17 不同振荡频率下表观雷诺数对参数不稳定区域的影响Fig.17 Effects of generalized Reynolds number on parametric instability region with different oscillation frequencies

减小。与K-H不稳定区域主导波数变化不同的是，各受迫频率下主导波数的分布较为均匀。

2.4 气液相对速度对液膜破碎长度的影响

前人的研究中[10-11]，破碎长度可以用于衡量液膜不稳定的程度。假设液膜表面振幅按照指数形式发展，当临界破碎时，表面波振幅ηb可表示为

ηb=η0exp(Br maxtb)

(35)

式中：η0为初始时刻表面波振幅；tb为破碎时间。

而液膜的破碎长度为

Lb=Ultb

(36)

因此，从式(35)中解出破碎时间，并代入式(36)，得到破碎长度与最大增长率间的关系：

(37)

在式(37)中，参考Jia等[11]的研究，将ln(ηb/η0)取为17，Br max可通过计算式(32)获得。选取Sivadas等[10]的实验数据与理论预测值进行对比，其中流动参数及气流振荡参数在文献[10]中已经给出。需要注意的是，在Sivadas等[10]的研究中工质为水，因此使用本文模型进行计算时，流变参数设置为K=1.005×10-3Pa·s、n=1，这样幂律流体即可退化为水进行计算。

ρg=1.225 kg/m3,ρl=998 kg/m3, a=2×10-4 m,σ=n=1,K=1.005×10-3 Pa·s, U0=10 m/s,ΔU=7.1 m/s,ω=750 Hz图18 气流有振荡和无振荡条件下破碎长度的理论值与实验值[10]对比Fig.18 Comparison of breakup length between theoretical value and experimental data[10] with oscillating and without oscillating gas

3 结 论

通过分析不同幂律性质下气流速度振荡幅值和振荡频率对液膜不同不稳定区域的影响，可得出以下结论：

1) 随着振荡幅值的增加，液膜出现了多个不稳定区域，K-H不稳定区域和参数不稳定区域的最大增长率及主导波数均随着振荡幅值的增加而增加。

2) 当表观雷诺数、幂律指数和无量纲速度因子逐渐增加时，K-H不稳定区域的最大增长率逐渐增加，而主导波数保持不变。在参数不稳定区域，最大增长率先减小后增加，且振荡幅值的变化不影响最大增长率转折点所对应的幂律性质。

3) 随着振荡频率的增加，参数不稳定区域逐渐向大波数范围移动，K-H不稳定区域的最大增长率、主导波数和截止波数增加；当振荡频率超过一定值后，只存在一个K-H不稳定区域。且振荡频率较小时，增长率及主导波数对振荡频率的变化较敏感。

4) 随着幂律指数和无量纲速度因子的逐渐增加，K-H不稳定区域的最大增长率逐渐增加，而主导波数保持不变；当无量纲振荡频率为0.050 时，参数不稳定区域的最大增长率单调增加，随着振荡频率的逐渐增大，出现最大增长率随着幂律指数或是无量纲速度因子的增加而先减小后增大的转折点，且此时转折点所对应的幂律指数或无量纲速度因子随着振荡频率的增加而逐渐增大。同时，主导波数随着幂律指数或无量纲速度因子的增加保持不变或逐渐减小。

5) 当表观雷诺数逐渐增加时，K-H不稳定区域的最大增长率单调增加，主导波数保持不变；而参数不稳定区域的最大增长率先减小后增加，振荡频率的增加使得最大增长率的转折点所对应的表观雷诺数逐渐增加，主导波数保持不变或逐渐减小。

6) 随着气液间相对速度的增加，液膜破碎长度逐渐减小，且相比于气流振荡的情况，无振荡条件下液膜更稳定，破碎长度更长。