Bingham 流體雙自由面熱毛細液層的穩(wěn)定性分析1)

王勝 胡開鑫

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

引言

當液體層表面存在水平溫度梯度時,其誘導的表面張力梯度將驅(qū)動液體運動.這種流動被稱為熱毛細對流.由于其在晶體生長[1]中的重要作用,熱毛細對流得到了廣泛的研究.Smith等[2]采用線性穩(wěn)定性方法研究了液層的熱毛細不穩(wěn)定.Patne等[3]對傾斜溫度梯度液層進行穩(wěn)定性分析.Chen等[4]研究了在零重力條件下由矩形池的等溫側(cè)壁引起的瞬態(tài)熱毛細對流.

近年來熱毛細效應的研究從單一自由面發(fā)展到雙自由面.美國宇航局宇航員Pettit[5]在國際空間站上進行一系列熱毛細流的微重力實驗.其中環(huán)形水的液膜暴露在不均勻的溫度分布中,這些實驗表明,具有兩個自由表面的液層有可能獲得一種新的材料結(jié)晶過程[6].Ueno等[7]通過數(shù)值和實驗方法研究在矩形孔液體薄膜中的熱毛細流動.Messmer等[8]研究在熱毛細作用力下的的雙自由表面薄膜,表明在低Marangoni 數(shù)下有兩種基本的流動結(jié)構(gòu).Toshiki等[9]研究具有兩個氣液界面的自由液體膜中的熱毛細流動,發(fā)現(xiàn)該流動表現(xiàn)出從二維穩(wěn)流態(tài)向三維振蕩態(tài)的轉(zhuǎn)變.Hu等[10]用線性穩(wěn)定性分析方法研究雙自由面熱毛細液層的失穩(wěn)機理.趙誠卓等[11]研究雙自由面溶質(zhì)-熱毛細液層的不穩(wěn)定性,發(fā)現(xiàn)在大Pr時,溶質(zhì)毛細力的出現(xiàn)使流動更加穩(wěn)定;在其他參數(shù)下,溶質(zhì)毛細力會減弱穩(wěn)定性.

非牛頓流體的熱毛細對流廣泛存在于鍍膜[12]、薄膜干燥[13]、脫濕[14]、光刻[15]、噴墨打印[16]和微重力聚合物加工[17]等應用中.非牛頓流體與牛頓流體在流動特性上有很大的差異.黏彈性熱毛細液層已被許多作者[18-20]研究過,結(jié)果表明:彈性導致的正應力會顯著影響流動穩(wěn)定性.Hu等[21]對剪切稀化流體的熱毛細液層進行了線性穩(wěn)定性分析.

黏塑性流體出現(xiàn)在許多工業(yè)應用和自然環(huán)境中,如聚合物加工[22]、熔漿[23]等.黏塑性流體的主要特征是其屈服應力,在應力較高時表現(xiàn)出類液體行為,在應力較低時表現(xiàn)出類固體行為.Bingham 流體是黏塑性流體的一種理想化模型.當剪切力未達到屈服值之前,Bingham 流體表現(xiàn)出彈性固體的性質(zhì)而沒有流動性,當超過屈服值之后,才發(fā)生流動且流體的黏度保持不變,類似牛頓流體[24].

近幾十年來,對Bingham 流體進行了許多的研究.Nouar等[25]從特征值敏感性、最小缺陷和轉(zhuǎn)捩標度律方面研究Bingham 流體在通道中流動的穩(wěn)定性.Hu等[26]研究無限大液層中Bingham 流體熱毛細對流的不穩(wěn)定性.王世芳等[27]研究Bingham 流體在低滲透多孔介質(zhì)中球向滲流的分形模型.韓建國等[28]使用同軸雙圓柱流變儀測試Bingham 流體流變參數(shù).Baioumy等[29]研究圓形管道入口區(qū)的Bingham 流體的流動.Fusi[30]提出一種基于有限差分的數(shù)值格式來研究Bingham 流體在通道中的非定常運動(平面泊肅葉流).Maurya等[31]研究Bingham 流體在垂直T 通道的子通道中的流動反轉(zhuǎn)情況.Esmaeili等[32]分析Bingham 流體在非平行板間的擠壓流動.Zhang等[33]利用射頻加熱法對稠油油藏開發(fā)中的Bingham 流體流動進行數(shù)值研究.

近年來,中國科學家開始對聚合物空間3D 打印技術(shù)進行試驗.Dou等[34]提出在微重力條件下設計和制造高精度陶瓷組件的方法,其中制備的陶瓷漿料的流變特性類似Bingham 流體.本文以陶瓷漿料的空間3D 打印為應用背景,研究Bingham 流體在雙自由面液體層中的熱毛細對流及其穩(wěn)定性.本文的組織結(jié)構(gòu)如下.在第1 節(jié)中,給出了該問題的物理模型和數(shù)值描述.推導出了基本流的解析解和控制方程.然后在第2 節(jié)中,得到了不同Bingham 數(shù)(B)和垂直溫差(Q)時的臨界參數(shù);顯示了臨界模態(tài)的流場圖,研究了能量機制.最后在第3 節(jié)進行總結(jié).

1 控制方程

本文研究的流動模型如圖1 所示.從狹縫中射出的聚合物流體液層,其表面的水平溫度梯度會導致內(nèi)部的熱毛細對流.d是液體層深度的一半,τ13是表面的剪應力,U0是速度場.在該參考系中,流體流向的平均流速為0.具有兩個自由表面的液層在其上下表面上受到恒定的溫度梯度b=-dT/dx>0的影響.其中x,y,z為流向、展向和法向.其中表面張力σ?隨溫度T的變化滿足 σ ?=σ?0-γ(T-T0) 的關(guān)系.由于層內(nèi)部的剪切速率小于表面附近的值,因此液層中間會有一個栓塞區(qū)域.該流動由3 個區(qū)域組成,III 是栓塞區(qū),I和II 是屈服區(qū)域,其范圍分別是(-z0,z0),[-1,-z0],[z0,1].

圖1 Bingham 流體雙自由面熱毛細液層示意圖Fig.1 The thermocapillary liquid layer with two free surfaces for a Bingham fluid

流動的無量綱控制方程組由下式給出,分別為連續(xù)性方程,動量方程和能量方程

其中u,p,T分別為速度、壓強和溫度,τ為應力張量.Reynolds 數(shù)Re,Marangoni 數(shù)Ma和Prandtl 數(shù)Pr分別定義為

采用Bingham 流體模型來反映聚合物的本構(gòu)關(guān)系.其本構(gòu)方程為

z=±1的界面為液層自由面,邊界條件為

其中速度場只有x方向的分量,Tb為垂直方向上的溫度分布.在液層整體參考系中,任意垂直截面上的質(zhì)量流量為0,可推出基本流,詳見附錄A.圖2是基本流的速度分布和垂直溫度分布,從圖2(a)可知,B越大,栓塞區(qū)域越大;從圖2(b)可知,隨著B的增加,垂直溫度梯度越來越小.

圖2 基本流的(a)速度分布和(b)垂直溫度分布(Ma=100,Q=0)Fig.2 (a) The velocity distribution and(b) vertical temperature distribution at Ma=100for the basic flow

式(13)中Q是上、下自由面垂直溫度差

設當Q=0時,Tb(±1)=0,根據(jù)邊界條件可得到Q1,Q2和Q的關(guān)系

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

定義臨界Marangoni數(shù)Mac為所有波數(shù)下中性Marangoni數(shù)MaN(σr=0)的最小值

本章討論了Q=0和Q＞0兩種情況,得到了不同參數(shù)下的臨界曲線,畫出了對應的擾動流場圖,并分析了其中的能量機制.

2.1 臨界曲線

2.1.1Q=0

本節(jié)畫出了Q=0,Bi=0時不同B下的臨界曲線.其中Mac隨Pr的變化如圖3(a)所示.可以發(fā)現(xiàn),Mac隨著B的增加而增大;當B=0時,流體為牛頓流體,圖中顯示有斜波和流向波兩種臨界模態(tài),在較小Pr,最不穩(wěn)定的模態(tài)是斜波,在大Pr下,流向波比斜波更不穩(wěn)定,因此臨界模態(tài)發(fā)生變化.臨界曲線的交叉處是兩種波中性MaN(σr=0)相等的情況.當B=0.01 時,作為Bingham 流體的臨界參數(shù)發(fā)生明顯突變,且在大Pr下,模態(tài)由同向流向波(θ=0°)變成了逆向斜波 θ ∈(90°,180°) ;在中小Pr下的臨界模態(tài)是逆向斜波,而在大Pr下則為逆向流向波(θ=180°);當B=0.6 時,臨界模態(tài)存在逆向斜波、同向斜波θ ∈(0°,90°)以及同向流向波3 種.結(jié)合圖3(b)綜合來看,k隨著B的增加而增大.在中小Pr下,牛頓流體的波數(shù)變化很小,而Bingham流體的波數(shù)隨著Pr的增加而明顯增大.圖3(c)為波傳播角 θ隨Pr變化的曲線,可以看出牛頓流體在Pr=1.2 的時候,θ 降為0°,之后隨著B的增加而變大.在大Pr的情況下,B=0.2 時,θ 變成180°,而對于B=0.6,θ 則變?yōu)?°.

圖3 Q=0,Bi=0時,不同B 下,(a) Mac隨Pr 的變化曲線以及(a)所對應的(b)波數(shù)和(c)波傳播角.曲線對應:a,c,d,e,g,h為斜波,b,f,i,j為流向波Fig.3 (a) The variation of Macwith Pr under different B at when Q=0and Bi=0,(b) wave number and(c) wave propagation angle.The curves correspond to:a,c,d,e,g,h oblique wave andb,f,i,j streamwise wave

此外可以發(fā)現(xiàn),Bi的增大并不會影響臨界模態(tài)的種類,而只會增強流動的穩(wěn)定性.

2.1.2Q＞0

本小節(jié)考慮上、下自由面有溫度差時的臨界模態(tài).圖4為B=0.2,Ma=100時垂直方向上的溫度分布,可以看出,當Q=0時,Tb是關(guān)于z=0對稱,Q的增加導致了Tb的不對稱.

圖4 B=0.2,Ma=100時垂直方向上的溫度分布Fig.4 Vertical temperature distribution at B=0.2 and Ma=100

圖5為Q=0.05,Bi=0時不同B下的臨界曲線.可以看出,牛頓流體存在逆向斜波、同向流向波、展向穩(wěn)態(tài)模態(tài)3 種臨界模態(tài),而Bingham 流體雖然在B=0.01 時存在與牛頓流體相同的臨界模態(tài),但在B=0.2時,臨界模態(tài)只有逆向斜波和展向穩(wěn)態(tài)模態(tài)兩種.對于斜波和流向波模態(tài),Mac基本隨著Pr的增加而增加,而展向穩(wěn)態(tài)模態(tài)則與之相反.此外對于同種模態(tài),B的增加對 θ 的影響較小.

圖5 Q=0.05 時,(a) Mac隨Pr 的變化曲線以及(a)所對應的(b)波數(shù)和(c)波傳播角.曲線對應:a,d,g為斜波,b,e為流向波,c,h,f為展向穩(wěn)態(tài)模態(tài)Fig.5 (a) The variation of Macwith Pr under Q=0.05,(b) wave number and(c) wave propagation angle.The curves correspond to:a,d,g oblique wave,b,e streamwise wave,c,h,f spanwise stationary mode

2.2 擾動流場

本節(jié)畫出了不同臨界模態(tài)下的等溫線和流線圖,并將最大擾動溫度進行歸一化.為了比較不同參數(shù)下的流場,將擾動的波長固定為 2π,此時橫軸數(shù)值表示擾動的相位,橫軸方向為波的傳播方向.

圖6 顯示了Q=0,Bi=0時臨界模態(tài)所對應的擾動流場.從圖6(b) 圖可以發(fā)現(xiàn),在B=0.2,Pr=0.01 時,擾動溫度分布在屈服區(qū)和栓塞區(qū),且呈對稱分布;在圖6(d)中可知,B=0.4,Pr=20時的擾動溫度呈反對稱分布,且擾動溫度只存在于屈服區(qū).

圖6 Q=0,Bi=0時不同臨界模態(tài)所對應的擾動流場:逆向斜波(B=0.2,Pr=0.01)的(a)上屈服面擾動,(b)溫度和流線;逆向斜波(B=0.4,Pr=20)的(c)上屈服面擾動,(d)溫度和流線Fig.6 The perturbation flow field of the different preferred modes at Q=0and Bi=0:(a) upper yield surface perturbation,(b) temperature and streamlines of the upstream oblique wave(B=0.2,Pr=0.01);(c) upper yield surface perturbation,(d) temperature and streamlines of the upstream oblique wave(B=0.4,Pr=20)

圖7 顯示了Q=0,Bi=5 時的擾動流場.從圖中可以看出,在小Pr下,擾動溫度呈對稱分布,隨著B的增加,擾動溫度在栓塞區(qū)呈截斷的狀態(tài).

圖7 Q=0,Bi=5 時不同臨界模態(tài)所對應的擾動流場:逆向斜波(B=0.4,Pr=0.01)的(a)上屈服面擾動,(b)溫度和流線;逆向斜波(B=0.6,Pr=0.01)的(c)上屈服面擾動,(d) 溫度和流線Fig.7 The perturbation flow field of the different preferred modes at Q=0and Bi=5:(a) upper yield surface perturbation,(b) temperature and streamlines of the upstream oblique wave(B=0.4,Pr=0.01);(c) upper yield surface perturbation,(d) temperature and streamlines of the upstream oblique wave(B=0.6,Pr=0.01)

圖8 顯示了Q=0.05,Bi=0時的擾動流場.從圖8(a)中可以看出牛頓流體的擾動溫度在整個流場中分布不規(guī)則,且主要集中在下半部分流場區(qū)域;而從圖8(c)中可知,Bingham 流體的擾動溫度分布在下半部分流場區(qū)域內(nèi).此外由于栓塞區(qū)的存在,Bingham 流體中無法找到如圖8(a)所示的貫穿液層的渦胞結(jié)構(gòu).

圖8 Q=0.05,Bi=0時不同臨界模態(tài)所對應的擾動流場:同向流向波(B=0,Pr=10)的(a)溫度和流線;同向斜波(B=0.2,Pr=0.1)的(b)上屈服面擾動,(c)溫度和流線Fig.8 The perturbation flow field of the different preferred modes at Q=0.05 and Bi=0:(a) temperature and streamlines of the downstream streamwise wave(B=0,Pr=10);(b) upper yield surface perturbation;(c) temperature and streamlines of the downstream oblique wave(B=0.2,Pr=0.1)

2.3 能量分析

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

式中N為擾動應力做的功,M為Marangoni 力在表面做的功,I為擾動流與基本流之間的相互作用.這里將擾動動能進行歸一化處理,即表1中給出了不同參數(shù)下各擾動能量變化項的值.可以發(fā)現(xiàn)在Q=0,Bi=0時,Bingham 流體的I在小Pr數(shù)下比牛頓流體的值大很多,因此不能忽略.然而Pr=100時,I則可以忽略不計,N＞0代表耗散,說明擾動能量的主要來源是Marangoni 力在表面做的功.

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

3 結(jié)論

本文采用線性穩(wěn)定性理論研究了Bingham 流體雙自由面熱毛細液層的穩(wěn)定性,分析了Prandtl 數(shù)、Bingham 數(shù)、垂直溫差(Q)和Biot 數(shù)(Bi)對流動穩(wěn)定性的影響,并結(jié)合流場圖和能量分析發(fā)現(xiàn)以下結(jié)論:

(1)從牛頓流體(B=0)變?yōu)锽ingham 流體(B＞0)時,因為流場產(chǎn)生了栓塞區(qū),臨界參數(shù)有明顯的突變,且在大Pr下,模態(tài)由同向流向波變?yōu)槟嫦蛐辈?牛頓流體臨界模態(tài)中貫穿液層的渦胞結(jié)構(gòu)也消失了;

(2)Mac隨著B和Bi的增加而增加,垂直溫差Q＞0時,Pr的增加會減弱流動的穩(wěn)定性;

(3)在小Pr情況下,擾動溫度存在于整個流場區(qū)域,在大Pr情況下,擾動溫度分布在屈服區(qū)域,栓塞區(qū)擾動溫度為0;

(4)擾動動能的主要能量來源是表面張力做功,但小Pr下基本流也有一定貢獻.

附錄A 基本流推導

將基本流形式代入動量方程式(2),在I和II 區(qū)域中,在x方向和z方向的動量方程如下

對于本構(gòu)方程,擾動應力由兩部分組成,一部分由擾動應變率引起,另一部分則由擾動黏性引起,即

當Q=0時,計算臨界參數(shù)將雙自由面分成對稱性和反對稱性兩種情況.其中對稱性邊界條件為

而反對稱性邊界條件則為