溶質Marangoni效應對降膜流動穩定性的影響

劉成治，李春曦，周靜宜，葉學民

（華北電力大學能源動力與機械工程學院， 河北 保定 071003）

引 言

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

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

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

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

1 理論模型

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

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

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

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

2 含表面活性劑的降膜流動線性穩定性分析

2.1 色散方程的推導

采用攝動法對基態流施加小擾動可得：

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

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

2.2 液膜流動中性穩定狀態

式中，kc表示截斷波數。液膜流動在0＜k＜kc為失穩區，其他區為穩定區。

2.3 線性穩定性分析

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

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

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

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

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

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

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

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

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

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

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

3 表面波演化歷程

3.1 初始條件與邊界條件

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

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

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

3.2 表面波演化特征

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

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

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

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

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

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

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

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

4 結 論

本文基于長波近似法推導出綜合考慮Weber數S、傾角θ、Peclet 數Pe、Marangoni 數M的色散關系式與臨界Reynolds 數Rec表達式，研究了溶質Marangoni效應作用下的降膜流動線性穩定性，并通過數值模擬探究了溶質Marangoni 效應作用下表面波的時空演化特征，主要結論如下。

（1）對于受溶質Marangoni 效應影響的降膜流動，提高M即增大了溶質Marangoni效應對液膜流動的影響。這將導致Rec的增大并降低失穩區內的擾動增長率與截斷波數kc，由此減小液膜失穩區范圍，抑制失穩區內擾動波的增長，有利于液膜流動的穩定，而增大Re與θ則降低了液膜流動的穩定性。

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

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