劉成治，李春曦，周靜宜，葉學(xué)民

（華北電力大學(xué)能源動力與機(jī)械工程學(xué)院， 河北 保定 071003）

引 言

降膜流動廣泛用于精密涂層[1]、燃燒室冷卻[2]、降膜換熱器[3]、生物醫(yī)療[4]等諸多領(lǐng)域，其流動穩(wěn)定性對于提高涂膜質(zhì)量、保障降膜換熱器換熱性能和促進(jìn)藥物吸收具有重要作用。因此，深入研究降膜流動的穩(wěn)定性問題具有重要的實際應(yīng)用價值[5-7]。

降膜流動是一種開式流體動力系統(tǒng)，具有顯著的對流不穩(wěn)定特征[8-10]。重力驅(qū)動下的降膜系統(tǒng)穩(wěn)定性可由臨界Reynolds 數(shù)Rec判別，當(dāng)液膜流動的Reynolds 數(shù)Re超過Rec后，擾動致使液膜流動不穩(wěn)定，最終導(dǎo)致液膜失穩(wěn)甚至混沌[11-13]。液膜受熱不均或添加表面活性劑后，在溫度梯度或活性劑濃度梯度的作用下，液膜表面將會產(chǎn)生表面張力梯度，該類現(xiàn)象分別稱為熱毛細(xì)效應(yīng)與溶質(zhì)Marangoni 效應(yīng)，此類效應(yīng)對液膜流動的穩(wěn)定性產(chǎn)生重要影響[14-18]，本文僅關(guān)注溶質(zhì)Marangoni效應(yīng)的影響。

目前，對于受溶質(zhì)Marangoni效應(yīng)影響的降膜流動穩(wěn)定性研究主要以線性穩(wěn)定性理論為主。Ji 等[19]通過求解O-S 方程，研究了表面活性劑對垂直壁面上降膜流動線性穩(wěn)定性的影響，結(jié)果表明：在中等Reynolds 數(shù)條件下，溶質(zhì)Marangoni 效應(yīng)在液膜表面引起了擴(kuò)散波，并證明了表面活性劑的解吸附是擴(kuò)散波產(chǎn)生的主要原因。Blyth 等[20]分析了傾斜平板上含表面活性劑的降膜流動線性穩(wěn)定性，結(jié)果表明：表面活性劑具有促進(jìn)液膜穩(wěn)定的作用，且添加表面活性劑可有效提高Rec。Pereira 等[21]考察了溶質(zhì)Marangoni效應(yīng)對降膜流動非線性區(qū)的影響，發(fā)現(xiàn)表面波的擾動演化為孤立脈沖，其振幅和速度隨溶質(zhì)Marangoni 效應(yīng)增強(qiáng)而減小。Georgantaki 等[22]實驗研究了SDS 活性劑對降膜流動穩(wěn)定性的影響，指出加入活性劑后，入口初始擾動顯著衰減，一次失穩(wěn)的Rec提高了一個數(shù)量級，并在一定參數(shù)下擾動會觸發(fā)二次失穩(wěn)。Hu 等[23]基于加權(quán)殘差模型研究了表面活性劑影響下的降膜流動線性穩(wěn)定性，并考慮了表面彈性和表面黏度對流動的影響，結(jié)果表明：提高表面彈性可降低擾動增長率，增大Rec，同時抑制表面波的色散效應(yīng)；提高表面黏度則降低了擾動增長率與截斷波數(shù)，但不影響Rec大小。

綜上所述，表面活性劑對降膜系統(tǒng)的穩(wěn)定性有重要影響。在以往的研究中[20-21,23]，所得臨界Reynolds數(shù)與Weber數(shù)無關(guān)，并極少考慮Peclet數(shù)對線性穩(wěn)定性的影響。而有研究表明[24]，無表面活性劑情形下，提高Weber 數(shù)可增大降膜流動穩(wěn)定區(qū)。因此，本文將在Blyth 等[20]研究基礎(chǔ)上，進(jìn)一步考慮Weber數(shù)與Peclet數(shù)對臨界Reynolds數(shù)的影響，提出包括Weber數(shù)、Marangoni數(shù)、壁面傾角和Peclet數(shù)在內(nèi)的臨界Reynolds 數(shù)表達(dá)式，進(jìn)而分析相關(guān)參數(shù)對流動線性穩(wěn)定性的影響，并通過數(shù)值模擬探究溶質(zhì)Marangoni 效應(yīng)作用下表面波的時空演化特征以及相關(guān)參數(shù)對流動非線性特征的影響。

1 理論模型

如圖1 所示，在傾角為θ的固體基底斜面上均勻附著一層厚度為h*的薄液膜，初始液膜厚度為，液膜受到擾動后表面呈現(xiàn)波狀，波長為l*。該液膜在重力作用下沿斜面向下流動，重力加速度為g，液膜表面分布有不溶性表面活性劑。假設(shè)該液膜為不可壓縮牛頓流體，其密度為ρ*，動力黏度為μ*。忽略液膜表面外界的氣體黏度與液膜蒸發(fā)的影響，不考慮分子間作用力，上角標(biāo)“*”代表有量綱量。

圖1 含不溶性活性劑薄液膜沿斜面下降流動示意圖Fig.1 Schematic diagram of a liquid film flow over an inclined plane in the presence of an insoluble surfactant

本文采用長波近似法研究中小Reynolds數(shù)范圍內(nèi)液膜的穩(wěn)定性，其優(yōu)勢是避免直接求解O-S 方程特征根所帶來的計算煩瑣的問題[8,21]。在長波近似法的假設(shè)下，本文基于Navier-Stokes 方程、連續(xù)性方程與活性劑擴(kuò)散方程[10]，并在基底處采用無滑移無滲透邊界條件，自由界面處考慮切向應(yīng)力、法向應(yīng)力平衡邊界條件以及運(yùn)動學(xué)邊界條件[21]，最終推導(dǎo)出如下形式的液膜厚度與活性劑濃度演化方程：

式中，Re為Reynolds 數(shù)，中小Reynolds 數(shù)條件下Re取值范圍為0～10；Pe為Peclet 數(shù)，表征表面活性劑的擴(kuò)散程度，其值與表面活性劑的種類有關(guān)，取值范圍為10-1～103；We為Weber 數(shù)，表征液膜表面張力與慣性力之比，本文引入S=ε2We表示經(jīng)縮小后的Weber數(shù)，取值范圍為0～10；M為Marangoni數(shù)，表征溶質(zhì)Marangoni 效應(yīng)的強(qiáng)度，取值范圍為0～10；U*為特征速度，U*=ρ*gh*2sinθ/3μ*；ε=h*/l*?1，表示液膜厚度與表面波波長的比值；Ds*為表面活性劑的擴(kuò)散系數(shù)，表征表面活性劑的擴(kuò)散強(qiáng)度。演化方程推導(dǎo)過程中假定表面張力與活性劑濃度之間符合線性關(guān)系式[25]：σ*(Γ*)=-α*(Γ*-)，其中是參考濃度，σm*為Γm*所對應(yīng)的表面張力，α*=-(?σ*/?Γ*)為表面張力系數(shù)，表征表面張力隨活性劑濃度變化的劇烈程度。

2 含表面活性劑的降膜流動線性穩(wěn)定性分析

2.1 色散方程的推導(dǎo)

采用攝動法對基態(tài)流施加小擾動可得：

將無窮小擾動按照正則模法展開為(η,γ) =(η?,γ?)eikx+ωt，其中復(fù)波速ω=ωr+iωi，ωr表示擾動增長率，ωi表示擾動頻率，實數(shù)k表示擾動波數(shù)，(η?,γ?)代表液膜厚度與活性劑濃度的擾動振幅，將展開后的無窮小擾動代入式(5)和式(6)可得色散關(guān)系式：

式(7)為以復(fù)波速ω為未知數(shù)的一元二次方程，其復(fù)數(shù)根的實數(shù)部分與虛數(shù)部分分別定義為擾動增長率ωr與擾動頻率ωi，具體如下：

2.2 液膜流動中性穩(wěn)定狀態(tài)

式中，kc表示截斷波數(shù)。液膜流動在0＜k＜kc為失穩(wěn)區(qū)，其他區(qū)為穩(wěn)定區(qū)。

2.3 線性穩(wěn)定性分析

下文主要關(guān)注Re、M、θ、Pe和S對液膜流動穩(wěn)定性的影響，若無特殊說明，Re= 1，M= 1，θ= 60°，Pe=100，S= 4.5。

圖2 和圖3 為最大擾動增長率ωr,m和最危險波數(shù)km隨Re的變化，表征流動的最不穩(wěn)定狀態(tài)。圖2表明：當(dāng)Re＜Rec時，ωr,m= 0，即擾動在任意波數(shù)下均不隨時間增長，流動保持穩(wěn)定；當(dāng)Re＞Rec時，ωr,m＞0，ωr,m隨Re的增大近似呈線性增長，衍生出失穩(wěn)區(qū)并隨Re增加而擴(kuò)大，最危險模式下液膜擾動增長越來越劇烈。增大M與減小θ降低了ωr,m的線性增長率，說明溶質(zhì)Marangoni 效應(yīng)的增強(qiáng)與慣性因素的減弱抑制了最危險模式下擾動的增長，有利于液膜的穩(wěn)定流動。

圖2 最大擾動增長率ωr,m隨M和θ 的變化（“ ”表示臨界狀態(tài)）Fig.2 Variation in the maximum perturbation growth rate ωr,m with M and θ (“ ” indicates a critical state)

圖3 表明：當(dāng)Re≤Rec時，km= 0，此時ωr,m≤ 0 在任意波數(shù)k時恒成立，液膜流動處于穩(wěn)定狀態(tài)。當(dāng)Re＞Rec時，最危險波數(shù)km隨Re增大而提高，液膜流動的最危險模式向短波區(qū)擴(kuò)張，表明增大Re不利于液膜的穩(wěn)定流動；增大θ與減小M均減弱了Marangoni應(yīng)力對慣性力的抑制作用，并使流動最危險模式向短波區(qū)移動，因此縮小了液膜的穩(wěn)定區(qū)，降低了液膜流動的穩(wěn)定性。例如無表面活性劑(M= 0)或基底豎直擺放時(θ= 90°)，Rec最小，此時km最大，液膜的失穩(wěn)區(qū)范圍最大。

圖3 最危險波數(shù)km隨M和θ 的變化（“ ”表示臨界狀態(tài)）Fig.3 Variation in the most dangerous wave number km with M and θ (“ ” indicates a critical state)

圖4 反映了Marangoni 數(shù)M對ωr+與相速度ci的影響。圖4 表明：提高M(jìn)減小了液膜在失穩(wěn)區(qū)的擾動增長率ωr+，同時使液膜截斷波數(shù)kc向長波區(qū)移動，從而減小了失穩(wěn)區(qū)內(nèi)液膜擾動增長與失穩(wěn)區(qū)面積，表明溶質(zhì)Marangoni效應(yīng)抑制了失穩(wěn)狀態(tài)下擾動的發(fā)展。在無表面活性劑(M= 0)情形下，相速度ci為一常數(shù)，表明表面波是非色散波，其相速度不隨擾動波數(shù)的變化而變化；添加表面活性劑后，隨波數(shù)k的增大相速度ci先降低后升高，表面波的演化受色散效應(yīng)影響。在長波區(qū)(k＜ 0.8)內(nèi)，ci在最危險波數(shù)km附近存在一極小值點，此時表面波的色散效應(yīng)最強(qiáng)，液膜擾動的發(fā)展最顯著。隨著M的增大，相速度曲線左移，ci極小值增大，表明溶質(zhì)Marangoni 效應(yīng)的進(jìn)一步增大將削弱表面波的色散。

圖4 擾動增長率ωr+和相速度ci隨M的變化Fig.4 Variation in perturbation growth rate ωr+ and phase speed ci with M

圖5 為不同參數(shù)下的中性穩(wěn)定曲線。由圖5 可知，增大M使中性穩(wěn)定性曲線下移，Rec增大，液膜失穩(wěn)區(qū)逐漸收縮，體現(xiàn)了溶質(zhì)Marangoni效應(yīng)對液膜流動的穩(wěn)定作用。提高θ則導(dǎo)致中性穩(wěn)定曲線上移，Rec減小，失穩(wěn)區(qū)擴(kuò)大。提高Weber 數(shù)S將減小失穩(wěn)區(qū)面積，促使流動穩(wěn)定，但不影響Rec的大小，這與Dholey 等[24]在研究無表面活性劑時降膜流動穩(wěn)定性中的結(jié)論一致。Pe表征活性劑單體擴(kuò)散速度的大小，與液膜表面活性劑分子的擴(kuò)散系數(shù)成反比。Pe越大，活性劑分子在液膜表面的擴(kuò)散速度越慢，Marangoni應(yīng)力對慣性力的阻礙作用越弱，液膜越趨于失穩(wěn)。

圖5 不同參數(shù)下的中性穩(wěn)定曲線（“ ”表示臨界狀態(tài)）Fig.5 Neutral stability under different parameters (“ ” indicates a critical state)

圖6對比了本文結(jié)果與Blyth等[20]基于直接求解O-S 方程所得結(jié)果，實線表示本文結(jié)果，虛線表示Blyth 等[20]的結(jié)果。圖6(a)表明，在Stokes 流假設(shè)下(Re= 0)，M= 0 時，兩種方法所得到的ωr+均小于0，液膜流動始終處于穩(wěn)定狀態(tài)，且兩條曲線在小波數(shù)范圍內(nèi)(k＜ 0.5)的數(shù)據(jù)基本吻合。圖6(b)表明，采用長波近似法所得出的中性穩(wěn)定曲線相較于直接求解O-S 方程得出的結(jié)果而言，曲線整體略微向x軸負(fù)半軸移動，流動穩(wěn)定區(qū)有所收縮，但變化幅度不大，M= 0與M= 1時對應(yīng)的Rec分別下降了0.41與0.50。

圖6 本文結(jié)果與文獻(xiàn)[20]結(jié)果的對比（“ ”表示臨界狀態(tài)）Fig.6 Comparison of the results of this paper with that of Ref.[20](“ ” indicates a critical state)

3 表面波演化歷程

3.1 初始條件與邊界條件

為進(jìn)一步探究溶質(zhì)Marangoni 效應(yīng)對表面波的影響，采用偏微分方程數(shù)值求解軟件“FreeFEM++v4.7”求解非線性演化方程組式(1)和式(2)，該軟件具有網(wǎng)格自適應(yīng)和移動網(wǎng)格生成的特點，目前已廣泛用于相關(guān)領(lǐng)域的研究中[26-29]。為模擬擾動的演化，根據(jù)演化尺度選取0＜x＜ 200的一維區(qū)域進(jìn)行計算，并設(shè)置初始擾動振幅為基本流膜厚的1%，如式(11)所示，相應(yīng)的邊界條件如式(12)所示。經(jīng)過網(wǎng)格無關(guān)性驗證后，在滿足精度與減少計算時間的條件下，設(shè)置計算域網(wǎng)格數(shù)為20000，時間步長為0.1。

初始液膜厚度與活性劑濃度分布[30]為：

平板兩端液膜厚度與活性劑濃度分布邊界條件為：

3.2 表面波演化特征

圖7 為溶質(zhì)Marangoni 效應(yīng)對表面波演化的影響。無表面活性劑(M=0)時，初始擾動最終演化為孤立波。添加表面活性劑(M＞0)后，表面波由于色散效應(yīng)的影響，衍生出快波與慢波。一方面，快波演化受溶質(zhì)Marangoni 效應(yīng)抑制，其振幅隨M的增大而減小，這與線性穩(wěn)定性結(jié)論一致；另一方面，在色散效應(yīng)影響下，慢波隨時間演化出了多峰波狀結(jié)構(gòu)，Georgantaki 等[22]的實驗同樣觀測到了這一現(xiàn)象。進(jìn)一步增大M則抑制了慢波的發(fā)展，使其振幅衰減，即溶質(zhì)Marangoni 效應(yīng)抑制了初始擾動的發(fā)展。

圖8 展示了溶質(zhì)Marangoni 效應(yīng)對活性劑濃度分布的影響。由圖8可知，隨著液膜的流動，活性劑分子逐漸向下游區(qū)域擴(kuò)散并呈現(xiàn)波狀分布。對比圖8(a)、(b)可知，隨著M的增加，活性劑濃度變化逐漸趨于平緩，即溶質(zhì)Marangoni效應(yīng)促進(jìn)了活性劑分子的重新分布。觀察圖7 可知，在M= 0 時，上游區(qū)液膜厚度分布均勻，下游區(qū)液膜初始擾動演化為快波；而當(dāng)M＞ 0時，隨著液膜的流動，上游區(qū)出現(xiàn)了慢波并且慢波振幅隨時間逐漸增長，而下游區(qū)的快波振幅則隨時間逐漸衰減，這實質(zhì)上是活性劑濃度分布不均所產(chǎn)生的Marangoni 應(yīng)力作用的結(jié)果。M相同時，圖7 中慢波區(qū)液膜厚度的波峰基本對應(yīng)圖8中活性劑濃度的波谷，二者相位差略小于180°，而圖7 中快波區(qū)液膜厚度的波峰對應(yīng)圖8 中活性劑濃度的波峰，Marangoni應(yīng)力驅(qū)動流體從高濃度處流向低濃度處，這一流動疊加液膜整體從左向右流動之后，導(dǎo)致慢波振幅隨時間增加，而快波振幅隨時間衰減，且慢波與快波間的距離不斷增大。

圖7 溶質(zhì)Marangoni效應(yīng)對表面波演化的影響（實線箭頭表示慢波，虛線箭頭表示快波）Fig.7 Influence of solutal Marangoni effect on surface waves evolution(solid arrow represents slow waves, and dashed arrow represents fast waves)

圖8 溶質(zhì)Marangoni效應(yīng)對活性劑濃度分布的影響Fig.8 Influence of solutal Marangoni effect on the distribution of surfactant concentration

圖9 為Peclet 數(shù)Pe對表面波演化的影響。圖9表明，提高Pe增大了慢波振幅，促使慢波后續(xù)演化出多峰波狀結(jié)構(gòu)。提高Pe即減緩了活性劑分子的擴(kuò)散速度，阻礙了活性劑分子在液膜表面的擴(kuò)散，且從活性劑濃度演化方程式(2)可知，Pe直接影響活性劑濃度Γ的二階導(dǎo)數(shù)項，因此在活性劑濃度變化劇烈的慢波區(qū)，Pe增大導(dǎo)致活性劑濃度梯度作用被放大，局部Marangoni 應(yīng)力提高，從而增大了慢波振幅，降低了液膜流動的穩(wěn)定性。Pe的變化對快波的演化基本無影響，這是因為快波區(qū)活性劑分布相對慢波區(qū)較為均勻，活性劑濃度Γ的二階導(dǎo)數(shù)項較小，因而Pe變化對快波的影響并不顯著。

圖9 Peclet數(shù)對表面波演化的影響（實線箭頭表示慢波，虛線箭頭表示快波）Fig.9 Influence of Peclet number on surface waves evolution (solid arrow represents slow waves, and dashed arrow represents fast waves)

圖10 為Weber 數(shù)S對表面波演化的影響。圖10 表明，降低S增強(qiáng)了慢波振幅，提高了慢波的波數(shù)，不利于液膜的穩(wěn)定。從液膜厚度和活性劑濃度的演化方程式(1)和式(2)可知，S直接影響液膜厚度的四階導(dǎo)數(shù)項，因此對液膜厚度變化劇烈的慢波區(qū)影響顯著，對液膜厚度變化平緩的快波區(qū)影響較小。由前文分析可知，在慢波區(qū)，液膜厚度與活性劑濃度存在略小于180°的相位差，同時，根據(jù)本文S的定義式，減小S意味著減小表面張力的作用，導(dǎo)致表面波演化更劇烈，在后續(xù)演化過程中慢波的波數(shù)將逐漸增加，因此慢波會演化為多峰結(jié)構(gòu)，慢波振幅也逐漸加大。這一點與圖7 中減小M的作用類似，但減小S的效果更為顯著。

圖10 Weber數(shù)對表面波演化的影響（實線箭頭表示慢波，虛線箭頭表示快波）Fig.10 Influence of Weber number on surface waves evolution (solid arrow represents slow waves, and dashed arrow represents fast waves)

4 結(jié) 論

本文基于長波近似法推導(dǎo)出綜合考慮Weber數(shù)S、傾角θ、Peclet 數(shù)Pe、Marangoni 數(shù)M的色散關(guān)系式與臨界Reynolds 數(shù)Rec表達(dá)式，研究了溶質(zhì)Marangoni效應(yīng)作用下的降膜流動線性穩(wěn)定性，并通過數(shù)值模擬探究了溶質(zhì)Marangoni 效應(yīng)作用下表面波的時空演化特征，主要結(jié)論如下。

（1）對于受溶質(zhì)Marangoni 效應(yīng)影響的降膜流動，提高M(jìn)即增大了溶質(zhì)Marangoni效應(yīng)對液膜流動的影響。這將導(dǎo)致Rec的增大并降低失穩(wěn)區(qū)內(nèi)的擾動增長率與截斷波數(shù)kc，由此減小液膜失穩(wěn)區(qū)范圍，抑制失穩(wěn)區(qū)內(nèi)擾動波的增長，有利于液膜流動的穩(wěn)定，而增大Re與θ則降低了液膜流動的穩(wěn)定性。

（2）S與Pe均對液膜流動的穩(wěn)定性產(chǎn)生重要影響：提高S可增大Marangoni 應(yīng)力對液膜流動的穩(wěn)定作用，縮小失穩(wěn)區(qū)范圍，但并不改變Rec；增大Pe可減緩活性劑分子在液膜表面的擴(kuò)散速度，阻礙活性劑分子的鋪展，不利于液膜的穩(wěn)定。

（3）表面波的演化結(jié)果表明：與無活性劑時液膜流動相比，溶質(zhì)Marangoni效應(yīng)引起了表面波的色散效應(yīng)，使初始擾動分離為快波與慢波，而慢波在后續(xù)演化中進(jìn)一步發(fā)展為多峰波狀結(jié)構(gòu)；進(jìn)一步增強(qiáng)溶質(zhì)Marangoni效應(yīng)則減小了活性劑濃度梯度，同時減小了快波與慢波的振幅，由此抑制了表面波的演化；提高Pe導(dǎo)致活性劑濃度梯度作用被放大，局部Marangoni應(yīng)力提高，從而降低了液膜流動的穩(wěn)定性；減小S意味著減小了表面張力的作用，導(dǎo)致表面波演化更劇烈，不利于液膜的穩(wěn)定流動。