雙自由面溶質-熱毛細液層的不穩定性1)

趙誠卓 胡開鑫

(寧波大學機械工程與力學學院,浙江寧波 315211)

引言

Marangoni 對流是由于流體界面存在表面張力梯度而產生的流動[1-2].它主要存在于空間微重力環境或小尺度流動等表面張力占主導而重力效應較弱的情況中[3].其中表面張力梯度由表面溫度不均導致的稱為熱毛細對流,由表面濃度不均導致的稱為溶質毛細對流.在非等溫雙組分溶液中,液體表面往往同時存在溫度梯度和濃度梯度,這種對流被稱為溶質-熱毛細對流[4-6].它廣泛存在于晶體生長[7-8]、微流控[9-10]、合金澆筑凝固[11-12]、有機薄液膜的生長[13-15]等多種工業應用中,因而引起大量研究者的興趣.

對溶質-熱毛細對流的相關研究已有30 多年.Kanouff 和Greif[16]在研究焊接過程中鐵熔體的熱毛細流動時,發現氧化物和硫化物雜質會影響熔體的溫度分布.Yasuhiro 等[11]、Arafune 和Hirata[12]對In-Ga-Sb 合金系統凝固過程中進行大量實驗研究,證實了溶質毛細效應對熱毛細效應具有重要影響.游仁然和胡文瑞[17]研究了浮區中熱和溶質的毛細對流,發現濃度Marangoni 數對浮區中的流場和濃度場有明顯的影響,對溫度場的影響相對較小.Cr?ll等[6]通過對矩形液池內硅-鍺(Si-Ge)熔體凝固過程中界面偏析引起的熱-溶質毛細對流的實驗研究發現:溶質毛細效應和熱毛細效應會發生對抗作用.McTaggart[18]發現溶質毛細力和熱毛細力的相互作用決定了雙組分溶質-熱毛細對流由穩態轉向非穩態.Chen 和Su[19]發現由于熱擴散和濃度擴散,在液層中會出現兩種非穩態流動.

當表面張力梯度超過一定值時,Marangoni 對流會發生失穩.例如液層在交流電場不穩定性[20]、大尺寸液橋熱毛細對流失穩[21]、旋轉磁場中熱質耦合對流不穩定性[22]等.近年來一些學者研究了溶質-熱毛細對流的不穩定性.由于熱毛細效應和溶質毛細效應的相互作用,不穩定性變得非常復雜[23].在晶體生長中,溶質毛細效應可以在溶質-熱毛細流動不穩定性中發揮積極作用[24].文獻[23-24]研究了雙組分溶液的熱-溶質毛細對流,發現流動失穩的自由表面會出現多種流場.陳捷超[25]研究了毛細力比對環形液池內耦合熱-溶質對流的影響,發現隨著毛細力比的增加,流動失穩的臨界數會逐漸減小.

上述溶質-熱毛細對流的研究主要針對單自由面液層.近年來Marangoni 對流的研究已拓展到雙自由面液層.NASA 的Pettit[26]在國際空間站進行了一系列熱毛細對流實驗,結果表明雙自由面熱毛細液層能獲得一種新的材料結晶過程.在工業應用中,溶質毛細效應和熱毛細效應經常同時出現[4-6],因此非常有必要對雙自由面液層溶質-熱耦合的流動進行研究.在理論分析中,可以將單自由面液層的模型[27]拓展到雙自由面液層.Yamamoto 等[28]對Pettit[26]的實驗進行了數值模擬,發現液層的幾何形狀是導致振蕩流動結構的關鍵因素.Hu 等[29]用線性穩定性分析方法研究了具有兩個自由表面的熱毛細液體層的失穩機理.上述實驗、數值模擬和理論工作發現:雙自由面液層與單自由面液層的Marangoni對流的失穩機制、擾動流場有顯著的差別.

到目前為止還沒有關于雙自由面液層溶質-熱毛細對流的研究,本文將研究此問題的不穩定性,在第1 節中給出了物理模型、無量綱控制方程;第2 節中用模態分析的方法先后研究了不同參數下的最優模態、擾動流場、能量分析和物理機制;第3 節總結了雙自由面溶質-熱毛細對流的失穩機制.

1 控制方程

考慮具有雙自由面的液層,當液層厚度 2d 遠小于液層流動的特征長度 L,例如 d=1 mm,L=20 mm[20],此時液層可近似看作無限大,如圖1 所示.液面水平方向上的溫度梯度(T1＞T2)和濃度梯度(C1＞C2)會導致表面張力梯度,進而驅動溶質-熱毛細對流.其中 x,y 和 z 分別代表流向、展向和法向,d 為液層厚度的一半,τ13為表面剪切力.液面的表面張力為可近似的表達為溫度 T 和濃度 C 線性函數[6,12,30]

圖1 雙自由面溶質-熱毛細液層對流的示意圖Fig.1 Schematic of the solutal-thermocapillary liquid layer with two free surfaces

其中 γT為表面張力溫度系數,γC為表面張力濃度系數,.在流動模型中還考慮了Soret 效應,即溫度梯度導致的濃度擴散;忽略Dufour 效應[31],即濃度梯度會導致熱擴散.

設液層為不可壓縮的牛頓流體,其無量綱控制方程組包括連續性方程(2)、動量方程(3)、能量方程(4)、傳質方程(5)和本構方程(6)

其中 u,p,T,C 分別為速度、壓強、溫度和濃度,τ為應力張量,S 為應變率張量.在傳質方程(5) 中表示Soret 效應.Reynolds 數 Re,Prandtl 數 Pr,Marangoni 數 Ma,Schmidt 數 S c,Soret 數 χ 分別定義為

其中 ρ 為流體密度,D 為溶質擴散系數,ξ 為熱擴散率,ζ 為二元混合物的Soret 系數.

上自由面處(z=1)應力和速度的邊界條件為

式中前兩項方程表示溶質-熱毛細效應,可由式(1)推得[25,30,32].第3 項表示自由面無變形.上自由面處的熱平衡和Soret 效應分別表示為[25,30,32]

類似地,可以得到下自由面處(z=-1) 的邊界條件

上述式中畢渥數 Bi=hd/k,其中 k,h 分別是導熱系數、表面換熱系數,在本文中,取 Bi=0 ;分別為液層上、下邊界處氣體溫度.它們是強加給環境的熱通量.這兩項是為了能量平衡而引入的,可以由基本狀態解來確定.

假設基本流為平行剪切流,得

根據質量守恒,液層在同一垂直平面內通量為0,即

溫度和濃度在 x 方向的分布是線性的,可得

其中在具體實驗中,可以在一個矩形方框的兩對邊分別施加熱源和冷源,另外兩對邊為絕熱層;同理,濃度梯度[12,26,28,35]也可以這樣設置.設邊界處的溫度、濃度為

由上述條件可推出基本流的解

它們有如下關系式

以下采用模態方法分析液層的流動穩定性.在基本流場中疊加正則小擾動

式中帶下標0 的變量表示基本流,無下標的變量表示擾動.式中 σ=σr+iσi為復頻率,σr為增長率,ω=|σi|為角頻率.α,β分別表示在x,y軸上的波數,總波數為k=,方向角本文中θ是以熱毛細對流為基準.擾動方程和邊界條件在附錄A中給出.σ 可以通過Chebyshev 配點法求解,配置點數 Nc通常設置為70～ 100.

2 數值結果

Mac越大表示穩定性更強.其中流動的參數估計如下所述,雙組分溶(液的)Pr的范圍是10-2～102,Sc的范圍[6,30,32]是.Platten 等[36-37]對于Soret系數實驗研究的評述,發現雙組分有機溶液和水基溶液的Soret 系數在 [O(10-3)～O(10-4)]K-1.表1 為甲苯/正己烷混合溶液的物性參數[23].

表1 298.15 K 時甲苯/正己烷混合溶液(C0=26.17%)的物性參數Table 1 Physical properties of toluene/n-henxane mixture at T0=298.15 K and C0=26.17%

根據量級估計,取 0.01 ≤Pr ≤100,S c=100,χ=0.1.毛細比 η[7,16,38]通常情況下是負數,即溶質毛細力和熱毛細力的方向相反.取 η=-0.5,-2 這兩種情況進行計算,分別代表熱毛細力大于和小于溶質毛細力但兩者處于同一量級的情況,并且和純熱毛細對流(η=0)進行對比.

2.1 臨界曲線

圖2 畫出了 η=0,-0.5,-2 在 0.01 ≤ Pr ≤100 下的臨界曲線.其中UOW (upstream oblique waves),SSM (spanwise stationary mode),DSW (downstream streamwise wave),USW (upstream streamwise wave)分別表示為逆向斜波(90°＜θ ＜180°,σi＜0)、展向穩態模態(θ=90°,σi=0)、同向流向波(θ=0°,σi＜0)、逆向流向波(θ=180°,σi＜0).

對于純熱毛細流(η=0)的情況,圖2(a)～圖2(d)顯示在 0.01 ≤Pr ≤111,Mac、波數k 和角頻率ω 隨Pr增加而增加.在 0.01 ≤Pr ≤1.1 時,臨界模態是逆向斜波,θ 隨 Pr的增加而增加;在1 .1 ≤Pr ≤111 時,臨界模態是同向流向波.在 111 ≤Pr ≤200 時,臨界模態是逆向流向波,Mac、波數 k 隨 Pr 增加而略微增加,ω 隨 Pr 增加而略微減小.

對于 η=-0.5 的情況,圖2 顯示 a2曲線對應0.01 ≤Pr ≤45.5,臨界模態是展向穩態模態,隨 Pr 的增加而增加,k=0.b2曲線對應 45.5 ≤Pr ≤62.4,臨界模態是展向穩態模態,和 k 是隨 Pr 增加而增加.c2曲線對應 62.4 ≤Pr ≤95,臨界模態是同向流向波,Mac,ω,k 是隨 Pr 的增加而增加.d2曲線對應95 ≤Pr ≤200,臨界模態是逆向流向波,是隨 Pr 的增加而減少,ω 和 k 隨 Pr 的增加而增加.

圖2 臨界參數隨 Pr 的變化曲線:(a) Mac ;(b)波傳播角 θ ;(c)波數 k ;(d)角頻率 ω .曲線對應:UOW (a1,e3),SSM (a2,b2,d3),DSW (b3,b1,c2,c3),USW (a3,c3,d2,c1,f3)Fig.2 The variation of critical parameter with Pr .(a) Mac,(b) wave propagation angle θ,(c) wave number k and (d) angular frequency ω.The curves correspond to:UOW (a1,e3),SSM (a2,b2,d3),DSW (b3,b1,c2,c3),USW (a3,c3,d2,c1,f3)

對于 η=-2 情況,圖2(b)顯示了在0.01 ≤Pr ≤131區間內,Mac是隨 Pr的增加而增加.a3曲線對應0.01 ≤Pr ≤35.6,臨界模態是逆向流向波,0 .56 ≤k ≤0.73,0.87 ≤ω ≤1.28.b3曲線對應 35.6 ≤Pr ≤83.4,臨界模態是同向流向波,k 和 ω 隨 Pr 的增加而減少.c3曲線對應 83.4 ≤Pr ≤131,臨界模態是逆向流向波,k=0.72,ω 隨 Pr 的增加而減少.在1 31 ≤Pr≤210 區間內,Mac是隨 Pr 的增加而減小.d3曲線對應1 31 ≤ Pr ≤185,臨界模態是展向穩態模態,k 隨 Pr 的增加而減少.e3曲線對應 185 ≤Pr ≤202,臨界模態是逆向斜波,ω,θ隨 Pr 的增加而增加,k 隨 Pr 的增加而減少.f3曲線對應 202 ≤Pr ≤210,臨界模態是逆向流向波,k=0.85,ω 隨 Pr 的增加而減少.

2.2 擾動流場

本文畫出了不同臨界模態下的等溫線圖、等濃度線圖和流線圖,并以最大擾動溫度為基準進行標準化.此時橫軸方向為波的傳播方向,其數值表示擾動相位 φ=αx+βy .

圖3 和圖4 分別顯示了 η=-0.5,-2 的不同臨界模態的擾動流場和濃度場.η=-0.5 的濃度場和溫度場在同一模態下的分布是相似的,區別就是擾動的幅值大小不同.圖3(a)、圖3(c)的流場分別和純熱毛細對流(η=0)處于 a1,b1曲線段處模態的流場相似.圖3(b)～圖3(d)顯示擾動場在同一周期內沿z=0呈反對稱分布;圖3(a)顯示展向穩態模態模態的濃度場在 z 上無明顯變化,流線顯示在同一周期內4 個渦,擾動場在同一周期內沿 z=0 呈對稱分布.圖3(b)～圖3(d)顯示分別展向穩態模態、流向波、逆流向波模態的擾動流場在同一周期內呈反對稱分布且有2 個渦.圖3(b)的濃度幅值在表面附近和內部區域都有分布.圖3(c)和圖3(d)的擾動幅值只分布在內部區域.圖4(a)顯示 η=-2,Pr=1 的逆流向波模態的溫度場和濃度場分布不相似,但是溫度和濃度幅值均分布在表面區域.Pr=50,1 00 的擾動濃度場和擾動溫度場分布是相似的.圖4(a)～圖4(c)的擾動場在同一周期內沿 z=0 呈反對稱分布且有2 個渦,擾動濃度幅值分布在內部區域.圖3(c)和圖4(c)、圖3(d)和圖4(b)的波傳播角、擾動流場的方向相反.

圖3 η=-0.5 時臨界模態的擾動流場.(a) SSM (Pr=0.01),(b) SSM (Pr=50),(c) DSW (Pr=70),(d) USW (Pr=100)Fig.3 The perturbation flow field of the different preferred modes at η=-0.5.(a) SSM (Pr=30),(b) SSM (Pr=0.01),(c) DSW(Pr=70),(d) USW (Pr=100)

圖4 η=-2 時不同臨界模態所對應的擾動流場.(a) USW (Pr=1),(b) DSW (Pr=50),(c) USW (Pr=100)Fig.4 The perturbation flow field of the different preferred modes at η=-2.(a) DSW (Pr=1),(b) USW (Pr=50),(c) DSW (Pr=100)

2.3 能量分析

擾動動能的變化率可以寫成如下形式[39]

其中 N 為擾動應力做的功,N ＞0 代表耗散;M 為表面毛細做功,I 為擾動流與基本流之間的相互作用.N,M,I具體為

為了進一步研究熱毛細效應和溶質毛細效應對做功的影響,可以把 M 寫成式(29)的形式,其中 MT,MC分別為熱、溶質毛細力做功.式(29)中分別為熱、溶質毛細力

表2 不同參數下各擾動能量變化項的值Table 2 Values of perturbation energy variation terms at different parameters.

結合圖1(b),從表2 中可以看出,擾動能量的來源是表面張力做功,擾動流與基本流的相互作用可以忽略.

在 η=-0.5 情況下:(1)曲線 a2,b2對應展向穩態模態,MT＜0,MC＞0,溶質毛細力做功占比隨 Pr 增加而減少,熱毛細力做功占比隨 Pr 增加而增加,并且在 Pr 較小的情況,熱毛細力做功可以忽略;(2)c2是順向流向波,MT＞0,MC＜0 .(3) d2段是逆向流向波,MT＞0,MC隨著 Pr 的增加由負變正.

在 η=-2 情況下:(1) a3段為逆向流向波,MT＜0,MC＞0.(2) b3段為同向流向波,MC＞0,MT是隨Pr增加而減少,MT隨 Pr 增加由正變負.(3) c3段為逆向流向波,MT＜0,MC＞0.(4) d3,e3,f3分別為展向穩態模態、逆向斜波、逆向流向波,MT＞0,MC＜0.

2.4 討論

在小 Pr 數情況下,圖3(a),圖3(b)和圖4(a)顯示,擾動的幅值分布在自由面是由于對流導致的;擾動溫度的幅值小于擾動濃度的幅值,能量分析顯示表面溶質毛細力做功大于熱毛細力做功,熱毛細力和溶質毛細力相互對抗.

在大 Pr 數情況下,圖3(c),圖3(d)和圖4(b),圖4(c)擾動幅值分布在中間區域是由于對流(U0T,U0C)和垂直對流自由面的擾動是由于熱擴散和濃度擴散導致的.在 η=-2 且小 Pr 數情況下,擾動溫度場和擾動濃度場分布不同的也是由熱擴散和濃度擴散導致的.

3 結論

本文采用線性穩定性分析方法研究了不同毛細比下的溶質-熱毛細對流隨Pr 變化的對流不穩定性,并結合流場圖和能量分析,發現以下結論:

(1)溶質-熱毛細對流的臨界模態相比于純熱毛細對流有較大的差別.前者分為同向流向波、逆向流向波、展向穩態模態、逆向斜波;而后者分為逆向斜波和同向流向波;

(2) Pr 和 η 對流動穩定性的影響并不是單調的.在中小 Pr 數下,Pr 數的增加會使 Mac增加;在大Pr數下,Pr 數的增加使 Mac減小.在較大 Pr 情況中,η的減小會使 Mac減小;在其他參數中,η 的變化會使Mac發生多種變化;

(3) 不同臨界模態的擾動場在同一周期沿著z=0的水平面呈反對稱分布或對稱分布.在η=-2且 Pr 較小時,擾動溫度場和擾動濃度場分布不相似;在其他參數下,擾動溫度場和擾動濃度場分布都是相似的;

(4)表面毛細力做功為流動主要能量來源,擾動流與基本流的相互作用可以忽略不計.在 η=-0.5 情況下,當 Pr 較大時,熱毛細力做功和溶質熱毛細力共同做正功;在其他參數下,熱毛細力做功和溶質熱毛細力做功符號相反.在 η=-2 時,溶質毛細力做和熱毛細力做功符號相反.

附錄

將式(24)代入控制方程、本構方程和邊界條件中,可得線性擾動方程

在 z 方向用Gauss-Lobatto 配點進行離散.在流場內部設置為;邊界處() 設置2 個配點.變量用級數展開,例如i=1,2,···,Ncz=±1

最終式(A1),式(A5)～ (A12)可以寫成 W g=σZ g 的形式.其中 W,Z 都為 (11Nc+8) 階方陣,其中廣義特征值 σ 和廣義特征值向量 g 可以通過Matlab 中的QZ 算法算出.