趙誠卓 胡開鑫

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

引言

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

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

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

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

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

1 控制方程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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).

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

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

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

2.2 擾動流場

本文畫出了不同臨界模態(tài)下的等溫線圖、等濃度線圖和流線圖,并以最大擾動溫度為基準(zhǔn)進(jìn)行標(biāo)準(zhǔn)化.此時(shí)橫軸方向?yàn)椴ǖ膫鞑シ较?其數(shù)值表示擾動相位 φ=αx+βy .

圖3 和圖4 分別顯示了 η=-0.5,-2 的不同臨界模態(tài)的擾動流場和濃度場.η=-0.5 的濃度場和溫度場在同一模態(tài)下的分布是相似的,區(qū)別就是擾動的幅值大小不同.圖3(a)、圖3(c)的流場分別和純熱毛細(xì)對流(η=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 個(gè)渦,擾動場在同一周期內(nèi)沿 z=0 呈對稱分布.圖3(b)～圖3(d)顯示分別展向穩(wěn)態(tài)模態(tài)、流向波、逆流向波模態(tài)的擾動流場在同一周期內(nèi)呈反對稱分布且有2 個(gè)渦.圖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 個(gè)渦,擾動濃度幅值分布在內(nèi)部區(qū)域.圖3(c)和圖4(c)、圖3(d)和圖4(b)的波傳播角、擾動流場的方向相反.

圖3 η=-0.5 時(shí)臨界模態(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 時(shí)不同臨界模態(tài)所對應(yīng)的擾動流場.(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 為擾動應(yīng)力做的功,N ＞0 代表耗散;M 為表面毛細(xì)做功,I 為擾動流與基本流之間的相互作用.N,M,I具體為

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

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

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

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

在 η=-2 情況下:(1) a3段為逆向流向波,MT＜0,MC＞0.(2) b3段為同向流向波,MC＞0,MT是隨Pr增加而減少,MT隨 Pr 增加由正變負(fù).(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)顯示,擾動的幅值分布在自由面是由于對流導(dǎo)致的;擾動溫度的幅值小于擾動濃度的幅值,能量分析顯示表面溶質(zhì)毛細(xì)力做功大于熱毛細(xì)力做功,熱毛細(xì)力和溶質(zhì)毛細(xì)力相互對抗.

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

3 結(jié)論

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

(1)溶質(zhì)-熱毛細(xì)對流的臨界模態(tài)相比于純熱毛細(xì)對流有較大的差別.前者分為同向流向波、逆向流向波、展向穩(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í),擾動溫度場和擾動濃度場分布不相似;在其他參數(shù)下,擾動溫度場和擾動濃度場分布都是相似的;

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

附錄

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

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

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