液膜冷却是液体火箭发动机常用的热防护方式，冷却剂的一种重要导入方式是通过多孔结构进入燃烧室，在热壁面形成冷却液膜^[1]。如果液膜波动很大，使高温气体直接接触室壁，会对发动机造成严重损害。因此，对液膜沿多孔壁面流动的稳定性研究具有重要意义。Benjamin^[2]和Yih^[3]首先对斜面上的液膜流动进行了线性稳定性分析，通过求解长波扰动近似下的Orr-Sommerfield方程获得不稳定的临界雷诺数Re_cr，即中性曲线上雷诺数的最小值。对于壁面为多孔结构的情况，普遍认为多孔介质中的流动满足达西定律。Beavers和Joseph^[4]通过实验研究发现多孔壁面内存在边界层，并提出用滑移边界条件替代壁面内边界层效果，这一条件被称为Beavers-Joseph滑移边界条件。采用孔隙特征长度远小于液膜厚度的假设，Pascal^[5]将牛顿流体层和多孔层的动力学进行解耦，使问题得到极大简化。采用Yih^[3]提出的求解方法进行稳定性分析，得到斜面渗透率的增大对液膜流动有不稳定化作用的结论。Pascal^[6]将这一工作扩展到非牛顿流体沿多孔斜面流动的稳定性研究，分析了幂律流体层表面滚动波的产生及结构演化。Sadiq和Usha^[7-8]利用长波理论获得了牛顿流体及幂律流体液膜厚度的非线性演化Benney方程。其线性稳定性分析结果与Pascal^[5-6]一致，将多孔渗透性的不稳定作用解释为流体层底部的速度梯度减小，从而使流体与壁面的摩擦减小。文献[9]采用切比雪夫配点法求解液膜沿多孔层流动的耦合问题，揭示了一种多孔模式的不稳定性，研究了厚度比、达西数和Beavers-Joseph系数对不稳定性的影响。使用非模态稳定性理论研究了三维扰动的线性稳定性^[10]。Ogden等^[11]针对基底为波纹状、加热环境和基底有渗透性等复杂情况的液膜贴壁流动，建立了权重残差模型，并通过数值模拟验证了其模型的正确性。Praveen Kumar等^[12]探究了倾斜多孔层上流动的双层液膜的稳定性。Deepu等^[13]考虑了多孔壁面渗透性的各向异性及不均匀性对液膜流动稳定性的影响。Barletta和Celli^[14]分析了倾斜多孔槽道中的对流不稳定性。综上所述，大部分研究集中在小雷诺数(Re＜Re_cr)。而推力室中的冷却液膜存在中等雷诺数(Re_cr＜Re＜100)状态。因此，本文将对中等雷诺数时多孔壁面上的液膜流动进行线性稳定性分析。
1 数学模型 考虑如图 1所示的二维流动，倾角为θ的无限长多孔斜面上有一层液膜，多孔壁面内充满同种牛顿流体。气液界面有小幅波动，虚线表示未受扰动的界面。液体有黏，气体无黏，气、液不可压。建立如图 1所示的直角坐标系。

图 1 液膜沿多孔壁面流动示意图 Fig. 1 Schematic diagram of liquid films flowing down a porous wall

图选项

将壁面上的液膜流动分解为定常基本流和扰动量，如下所示：

其中:u、v、p和h分别为速度的x向分量、y向分量、压强和液膜厚度；U、V、P和H为基本流的物理量；u′、v′、p′和h′为扰动量。
定常基本流速度沿x方向，控制方程(下标为x，y，t时表示求偏导)为

(1)

式中：ρ为流体密度；g为重力加速度；μ为流体动力黏度。
气液界面的边界条件为

(2)

壁面处采用Beavers-Joseph滑移边界条件^[4]：

式中：u_p和v_p分别为多孔介质中x方向和y方向的达西平均渗透速度；α为与多孔介质的结构相关的无量纲参数；κ为多孔介质的渗透率。
多孔介质的孔隙特征长度远小于液膜厚度时，应用Pascal^[5]的尺度分析将壁面条件简化为

则基本流的壁面边界条件为

(3)

由式(1)~式(3)解得基本流的速度和压强分布为

(4)

考虑施加小扰动(u′, v′, p′, h′)后的非定常流动，忽略表面张力，线性化的控制方程为

(5)

气液界面处，线性化的运动边界条件为

(6)

线性化后动力边界条件的切向及法向方程为

(7)

壁面无滑移边界条件为

(8)

利用下列变换式对式(4)~式(8)进行无量纲化：

其中：为壁面不可渗透时液膜基本流的平均速度；上标“*”表示无量纲量。得到无量纲化的基本流解及扰动方程组：

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

式中：为液膜雷诺数；为表征多孔介质渗透性的无量纲参数。式(9)~式(17)即多孔壁面上液膜流动的小扰动数学模型。
2 线性稳定性分析 研究者广泛采用的长波近似方法限制在小雷诺数的范围。本文采用动量积分方法分析中等雷诺数时液膜的稳定性。考虑时间模式，假设正则模形式的扰动解为

(18)

式中：k为实数，表示波数；ω为复数，其实部ω_r为频率，虚部ω_i为时间增长率；c.c.表示复共轭。
已求得基本流轴向速度型为半抛物线型，认为扰动解的轴向速度型与之类似，为二次函数型：

(19)

由连续方程式(10)得

(20)

将式(18)~式(20)代入边界条件式(13)~式(16)中，可得常系数c₁、c₂和c₃的表达式分别为

(21)

对动量方程式(12)求y的定积分，从1积到y，并利用法向动力边界条件式(17)求得p′为

(22)

将式(22)代入动量方程式(11)中，求y的定积分，从0积到1，得液膜厚度内的平均动量方程。将式(21)代入，得色散方程为

(23)

式中：

3 结果分析 3.1 色散关系 利用已求得的色散方程式(23)，使用MATLAB中的“solve”函数可求得每一k值下的ω，取增长率最大的解，进而获得如图 2所示的小雷诺数下色散曲线。固定斜面倾角θ和多孔介质参数β的值。横坐标为波数k，纵坐标为时间增长率ω_i。

图 2 小雷诺数的色散曲线(θ=45°, β=0.1) Fig. 2 Dispersion curves with small Reynolds number(θ=45°, β=0.1)

图选项

时间增长率表示由小扰动引起的液膜表面波动幅值随时间增长的快慢。由图 2看出，随波数k增大，时间增长率ω_i呈现出2种趋势：在雷诺数Re很小的范围内，ω_i随k单调减小，始终保持负值，液膜对任意波长的扰动都是稳定的；当Re大于某个值时，ω_i先增大为正值，后减小到负值，表明液膜对某一范围内的波数是不稳定的。Re=0.69作为这2种趋势的一个临界情形，被称作临界雷诺数Re_cr，是判定流动稳定性的重要参数。这一数值与Pascal^[6]利用长波理论得到的同种工况的Re_cr≈0.67很接近。
对于中等雷诺数的色散曲线(见图 3)，增长率的最大值不再单调地随雷诺数增大而增大。当Re≥6时，雷诺数越大，最大增长率越小。截止波数，即增长率从正值变为负值时的波数，随雷诺数的变化是单调递增的。雷诺数越大，不稳定波数区域越大。

图 3 中等雷诺数的色散曲线(θ=45°, β=0.1) Fig. 3 Dispersion curves with moderate Reynolds number(θ=45°, β=0.1)

图选项

图 4为最大增长率ω_im随雷诺数Re变化的曲线。ω_im在小于临界雷诺数的小雷诺数范围单调增大，在中等雷诺数范围有极值。β越大，极值越大，且达到极值的雷诺数越小。

图 4 最大增长率随雷诺数变化(θ=45°) Fig. 4 Change of maximum growth rate with Reynolds number(θ=45°)

图选项

判别稳定性最直观的曲线是中性稳定曲线，曲线上每一点的增长率都为零。如图 5为利用色散关系做出的中性稳定曲线。

图 5 中性稳定曲线(θ=45°) Fig. 5 Neutral stability curves(θ=45°)

图选项

中性稳定曲线包围的右下方代表不稳定区域，其余为稳定区域。曲线上雷诺数的最小值是临界雷诺数。如图 5所示，每条中性稳定曲线的临界雷诺数都在曲线与横坐标的交点处。3条曲线分别对应不同的β值。β越大，临界雷诺数越小。因此，多孔壁面的渗透性对液膜流动起着不稳定化的作用。在雷诺数较大(Re＞9)时，β越大，对应的不稳定波数范围越小。
ω的实部ω_r除以波数k是扰动行波的相速度c。图 6为波速c随波数k变化的曲线。从整体趋势来看，c随k增大而减小，表明波长越长的波移动越快。考虑多孔介质渗透性的影响，对比不同β值的曲线，β越大，c越大。因此多孔壁面的渗透性加快了行波的移动。

图 6 波速随波数变化(θ=45°, Re=10) Fig. 6 Change of wave speed with wave number(θ=45°, Re=10)

图选项

3.2 能量分析 为探究各因素对流动不稳定性的作用机理，在获得速度场的基础上进行能量分析^[15]。用扰动速度乘以动量方程，并在一个波长和一个液膜厚度的区域内做平均，得到能量方程为

(24)

式中：KIN、REY、HYD、SHE和DIS分别为动能变化率、雷诺应力做功功率、重力做功功率、剪切应力做功功率和黏性耗散，表达式分别为

其中：λ为波长。
图 7为能量分析图，为能量变化率。由图 7(a)可得，剪切应力做功功率与黏性耗散相较于动能变化率和其他2种力的做功功率大得多。剪切应力做功来自于气液界面的波动，始终为正值，促进了扰动动能的增大。黏性耗散始终为负值，抑制失稳。图 7(b)所示为放大后其余3种能量的变化情况。扰动动能变化率KIN表征了流动稳定情况。k小于截止波数时KIN大于零，流动不稳定，k大于截止波数时KIN小于零，流动稳定。雷诺应力做功表示基本流和扰动之间的能量输运。雷诺应力和重力均在动能变化率大于零时做负功，在动能变化率小于零时做正功，这两者始终抑制扰动动能的改变。

图 7 能量分析图(θ=45°, Re=5, β=0.1) Fig. 7 Energy analysis chart(θ=45°, Re=5, β=0.1)

图选项

为探究多孔介质渗透性对液膜流动稳定性的影响机制，对比不同β值的动能变化率和各力做功功率，如图 8所示。图中箭头所指表示随着β增大，该种能量变化率的变化趋势。

图 8 不同β值时动能变化率及各力做功功率随波数的变化(θ=45°, Re=5) Fig. 8 Rates of change of kinetic energy and change of power of various forces with wave number under different β values(θ=45°, Re=5)

图选项

动能变化率随β增大而增大，对应于3.1节得到的结论，多孔壁面的渗透性促进流动失稳。造成动能变化率改变的原因是式(24)等号右边部分的改变。雷诺应力和重力做负功的功率是增大趋势，剪切应力做正功的功率是减小趋势，它们都有使动能变化率减小的作用；负的黏性耗散减小，起到增大动能变化率的效应。也就是说，多孔壁面渗透性的增强使得粘性耗散减小得如此剧烈以至于比其他3种力做功功率的变化之和还要大，从而导致动能变化率增大。这可能是由于多孔壁面处流体速度存在滑移，使得速度梯度下降，但平均流速增大，所以只与平均流速相关的雷诺应力和重力做功功率增大，而与速度梯度相关的剪切应力做功功率和黏性耗散减小。黏性耗散的变化最大，掩盖了其他效应，最终导致动能变化率的增大。
4 结论 1) 多孔壁面渗透性促进了液膜流动的不稳定，加快了液膜表面波的移动。
2) 中等雷诺数范围的最大增长率随雷诺数增大具有先增大后减小的趋势，不同于小雷诺数时单调增长的情况。
3) 渗透性越强的壁面，最大增长率的极值越大，达到极值的雷诺数越小。
4) 不稳定波数区域随雷诺数增大单调增大，随壁面渗透性增强具有先增大后减小趋势。
5) 多孔壁面渗透性使黏性耗散减小剧烈，超过其他能量的变化，使得综合效果为动能变化率增大，流动不稳定性增强。

参考文献

[1]	周红玲, 杨成虎, 刘犇. 液体火箭发动机液膜冷却研究综述[J].载人航天, 2012, 18(4): 8–13. ZHOU H L, YANG C H, LIU B. Review of research on liquid film cooling for liquid-propellant rocket engine[J].Manned Spaceflight, 2012, 18(4): 8–13.(in Chinese)
[2]	BENJAMIN T B. Wave formation in laminar flow down an inclined plane[J].Journal of Fluid Mechanics, 1957, 2(6): 554–574.DOI:10.1017/S0022112057000373
[3]	YIH C S. Stability of liquid flow down an inclined plane[J].Physics of Fluids, 1963, 6(3): 321–334.DOI:10.1063/1.1706737
[4]	BEAVERS G S, JOSEPH D D. Boundary conditions at a naturally permeable wall[J].Journal of Fluid Mechanics, 1967, 30(1): 197–207.DOI:10.1017/S0022112067001375
[5]	PASCAL J P. Linear stability of fluid flow down a porous inclined plane[J].Journal of Physics D:Applied Physics, 1999, 32(4): 417–422.DOI:10.1088/0022-3727/32/4/011
[6]	PASCAL J P. Instability of power-law fluid flow down a porous incline[J].Journal of Non-Newtonian Fluid Mechanics, 2006, 133(2): 109–120.
[7]	SADIQ I M R, USHA R. Thin Newtonian film flow down a porous inclined plane:Stability analysis[J].Physics of Fluids, 2008, 20(2): 022105.DOI:10.1063/1.2841363
[8]	SADIQ I M R, USHA R. Effect of permeability on the instability of a non-Newtonian film down a porous inclined plane[J].Journal of Non-Newtonian Fluid Mechanics, 2010, 165(19): 1171–1188.
[9]	LIU R, LIU Q S. Instabilities of a liquid film flowing down an inclined porous plane[J].Physical Review E Statistical Nonlinear & Soft Matter Physics, 2009, 80(3 Pt 2): 036316.
[10]	LIU R, LIU Q S. Instabilities and transient behaviors of a liquid film flowing down a porous inclined plane[J].Physics of Fluids, 2010, 22(7): 074101.DOI:10.1063/1.3455503
[11]	OGDEN K A, D'ALESSIO S J D, PASCAL J P. Gravity-driven flow over heated, porous, wavy surfaces[J].Physics of Fluids, 2011, 23(12): 122102.DOI:10.1063/1.3667267
[12]	PRAVEEN KUMAR A A, USHA R, BANERJEE T, et al. Instabilities of a free bilayer flowing on an inclined porous medium[J].Physical Review E Statistical Nonlinear & Soft Matter Physics, 2013, 88(6): 063012.
[13]	DEEPU P, SRINIVAS K, PRATEEK A, et al. Stability of a liquid film flowing down an inclined anisotropic and inhomogeneous porous layer:An analytical description[J].Journal of Fluid Mechanics, 2016, 807: 135–154.DOI:10.1017/jfm.2016.613
[14]	BARLETTA A, CELLI M. Instability of combined forced and free flow in an inclined porous channel[J].International Journal of Computational Methods, 2016, 13(2): 1640001.DOI:10.1142/S0219876216400016
[15]	KALLIADASIS S, RUYER-QUIL C, SCHEID B, et al. Falling liquid films[M].Berlin: Springer, 2012: 57-60.