雙自由面溶質(zhì)-熱毛細液層的不穩(wěn)定性1)

趙誠卓 胡開鑫

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

引言

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

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

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

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

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

1 控制方程

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

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

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

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

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

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

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

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

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

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

假設(shè)基本流為平行剪切流,得

根據(jù)質(zhì)量守恒,液層在同一垂直平面內(nèi)通量為0,即

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

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

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

它們有如下關(guān)系式

以下采用模態(tài)方法分析液層的流動穩(wěn)定性.在基本流場中疊加正則小擾動

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

2 數(shù)值結(jié)果

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

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

根據(jù)量級估計,取 0.01 ≤Pr ≤100,S c=100,χ=0.1.毛細比 η[7,16,38]通常情況下是負數(shù),即溶質(zhì)毛細力和熱毛細力的方向相反.取 η=-0.5,-2 這兩種情況進行計算,分別代表熱毛細力大于和小于溶質(zhì)毛細力但兩者處于同一量級的情況,并且和純熱毛細對流(η=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)、展向穩(wěn)態(tài)模態(tài)(θ=90°,σi=0)、同向流向波(θ=0°,σi＜0)、逆向流向波(θ=180°,σi＜0).

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

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

圖2 臨界參數(shù)隨 Pr 的變化曲線:(a) Mac ;(b)波傳播角 θ ;(c)波數(shù) 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區(qū)間內(nèi),Mac是隨 Pr的增加而增加.a3曲線對應0.01 ≤Pr ≤35.6,臨界模態(tài)是逆向流向波,0 .56 ≤k ≤0.73,0.87 ≤ω ≤1.28.b3曲線對應 35.6 ≤Pr ≤83.4,臨界模態(tài)是同向流向波,k 和 ω 隨 Pr 的增加而減少.c3曲線對應 83.4 ≤Pr ≤131,臨界模態(tài)是逆向流向波,k=0.72,ω 隨 Pr 的增加而減少.在1 31 ≤Pr≤210 區(qū)間內(nèi),Mac是隨 Pr 的增加而減小.d3曲線對應1 31 ≤ Pr ≤185,臨界模態(tài)是展向穩(wěn)態(tài)模態(tài),k 隨 Pr 的增加而減少.e3曲線對應 185 ≤Pr ≤202,臨界模態(tài)是逆向斜波,ω,θ隨 Pr 的增加而增加,k 隨 Pr 的增加而減少.f3曲線對應 202 ≤Pr ≤210,臨界模態(tài)是逆向流向波,k=0.85,ω 隨 Pr 的增加而減少.

2.2 擾動流場

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

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

圖3 η=-0.5 時臨界模態(tài)的擾動流場.(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 時不同臨界模態(tài)所對應的擾動流場.(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具體為

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

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

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

在 η=-0.5 情況下:(1)曲線 a2,b2對應展向穩(wěn)態(tài)模態(tài),MT＜0,MC＞0,溶質(zhì)毛細力做功占比隨 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分別為展向穩(wěn)態(tài)模態(tài)、逆向斜波、逆向流向波,MT＞0,MC＜0.

2.4 討論

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

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

3 結(jié)論

本文采用線性穩(wěn)定性分析方法研究了不同毛細比下的溶質(zhì)-熱毛細對流隨Pr 變化的對流不穩(wěn)定性,并結(jié)合流場圖和能量分析,發(fā)現(xiàn)以下結(jié)論:

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

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

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

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

附錄

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

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

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