雙層系統(tǒng)內(nèi)Jeffreys流體的Rayleigh-Marangoni對(duì)流不穩(wěn)定性分析

欒 致 漫

(南京航空航天大學(xué) 數(shù)學(xué)系, 南京 210016； 南京航空航天大學(xué) 飛行器數(shù)學(xué)建模與高性能計(jì)算工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室, 南京 210016)

雙層系統(tǒng)中的Rayleigh對(duì)流穩(wěn)定性問(wèn)題在工程領(lǐng)域應(yīng)用廣泛, 如地?zé)崴畮?kù)、 地下水流動(dòng)、 生物材料和核電廠的散熱等[1-7]. Nield[1]研究了底部加熱的多孔層-流體層的線(xiàn)性穩(wěn)定性問(wèn)題, 并考慮在可變形上表面的表面張力效應(yīng), 得到了恒定熱通量邊界條件的解; Chen等[2]研究了流體-多孔介質(zhì)系統(tǒng)中的指狀對(duì)流, 結(jié)果表明, Rayleigh數(shù)與波數(shù)的曲線(xiàn)為雙模態(tài), 即存在2個(gè)局部最小值, 當(dāng)厚度比較小時(shí), Rayleigh對(duì)流不穩(wěn)定性主要由多孔層主導(dǎo), 當(dāng)厚度比較大時(shí), Rayleigh對(duì)流不穩(wěn)定性主要由流體層控制； Straughan[3]利用Jones條件和Beavers-Joesph條件組成系統(tǒng)中交界面處的邊界方程, 通過(guò)改變多組物理參數(shù), 分析了雙層系統(tǒng)中的熱穩(wěn)定性； Chang[4-5]分析了雙層系統(tǒng)中Couette流和Poiseuille流的Rayleigh對(duì)流的線(xiàn)性穩(wěn)定性, 利用Chebyshev tau-QZ方法獲得相關(guān)數(shù)值結(jié)果, 并討論了水平方向與垂直方向穩(wěn)定性的差異； Samanta[6]用Darcy-Brinkman方程描述了多孔層的流動(dòng), 進(jìn)而考察了雙層系統(tǒng)的Couette流-Poiseuille流的穩(wěn)定性問(wèn)題, 并分析了厚度比對(duì)流體層和多孔層穩(wěn)定性的影響； Kolchanova[7]研究了失重狀態(tài)下雙層系統(tǒng)發(fā)生高頻縱向振動(dòng)時(shí)Rayleigh對(duì)流的產(chǎn)生, 結(jié)果表明, 對(duì)流速度與厚度比有關(guān); Block[8]認(rèn)為薄液層中的表面張力導(dǎo)致對(duì)流細(xì)胞改變; Pearson[9]認(rèn)為薄液層中的表面張力梯度產(chǎn)生了對(duì)流, 這種由表面張力驅(qū)動(dòng)的對(duì)流稱(chēng)為Marangoni對(duì)流; 劉秋生[10]比較了在流體-流體雙層系統(tǒng)中Marangoni對(duì)流與熱毛細(xì)對(duì)流的不穩(wěn)定性, 通過(guò)數(shù)值模擬可知, 穩(wěn)定狀態(tài)與界面張力溫度系數(shù)有關(guān); 劉榮等[11-12]對(duì)帶有蒸發(fā)界面的雙層系統(tǒng)進(jìn)行了Rayleigh-Marangoni對(duì)流的穩(wěn)定性研究, 結(jié)果表明, 增大蒸發(fā)速率, 不穩(wěn)定性集中于液體層的底部界面, 蒸發(fā)系數(shù)對(duì)汽-液雙層系統(tǒng)的影響較大; Liu等[13]通過(guò)考慮不同的厚度比, 研究了流體-多孔介質(zhì)系統(tǒng)中Rayleigh對(duì)流和Marangoni對(duì)流耦合模式的不穩(wěn)定性； Zhao等[14]分析了流體-多孔介質(zhì)系統(tǒng)中自由表面?zhèn)鳠釛l件對(duì)Rayleigh-Marangoni不穩(wěn)定性的熱效應(yīng), 結(jié)果表明, 不穩(wěn)定性模態(tài)轉(zhuǎn)變僅出現(xiàn)在表面張力較強(qiáng)的條件下； Yin等[15]研究了雙層系統(tǒng)中Oldroyd-B流體的熱不穩(wěn)定性, 并通過(guò)不同的參數(shù)對(duì)比分析穩(wěn)態(tài)對(duì)流與振蕩對(duì)流; 康建宏等[16]分析了多孔介質(zhì)內(nèi)黏彈性流體的Rayleigh對(duì)流穩(wěn)定性; Yin等[17]研究了雙層系統(tǒng)中黏彈性流體的熱不穩(wěn)定性; Patne等[18]研究了Jeffreys流體中的Marangoni不穩(wěn)定性; Sarma等[19-20]討論了在可變形和不可變形上表面Maxwell流體的Marangoni不穩(wěn)定性; 章紹能等[21]對(duì)溫度不均勻界面上黏彈性液滴的單層系統(tǒng)進(jìn)行了流動(dòng)穩(wěn)定性分析, 其中考慮了流體彈性因素的影響, 結(jié)果表明, 當(dāng)Prandtl數(shù)較小時(shí), 液體層趨于穩(wěn)定模態(tài); Madhukesh等[22]研究了在化學(xué)作用下Casson納米流體的Marangoni對(duì)流, 通過(guò)改變溫度、 微生物和濃度等因素分析了流動(dòng)的穩(wěn)定性; Li等[23]研究了混合納米流體中的Marangoni對(duì)流問(wèn)題, 將Darcy-Forchheimer流添加到動(dòng)量方程中, 并在熱方程中分析了具有非線(xiàn)性熱源匯與熱輻射的影響. 但上述研究?jī)H討論了Rayleigh對(duì)流和Marangoni對(duì)流的單一對(duì)流現(xiàn)象, 并未同時(shí)考慮兩類(lèi)對(duì)流.

本文在上述研究工作的基礎(chǔ)上對(duì)比了Rayleigh對(duì)流和Marangoni對(duì)流, 并將兩類(lèi)對(duì)流耦合, 研究流體-多孔介質(zhì)雙層系統(tǒng)中非牛頓流體的Rayleigh-Marangoni對(duì)流不穩(wěn)定性問(wèn)題. 通過(guò)線(xiàn)性穩(wěn)定性將3維偏微分方程組簡(jiǎn)化為常微分方程, 在振蕩對(duì)流下, 分析應(yīng)力松弛時(shí)間、 應(yīng)變弛豫時(shí)間、 上表面?zhèn)鲗?dǎo)對(duì)流特性和厚度比對(duì)Rayleigh對(duì)流和Marangoni對(duì)流的影響.

1 數(shù)學(xué)模型

圖1為流體-多孔介質(zhì)雙層系統(tǒng)的示意圖. 考慮覆蓋在多孔層z∈(-dm,0)的非牛頓流體層z∈(0,d)組成的流體-多孔介質(zhì)雙層系統(tǒng), 流體層和多孔層的交界面設(shè)為z=0.多孔層底部溫度設(shè)為T(mén)L, 流體層的上邊界是一個(gè)沒(méi)有任何變形的自由表面, 溫度為T(mén)U, 交界面處溫度為T(mén)I.溫度關(guān)系滿(mǎn)足TL>TI=Tex>TU. 流體層和多孔層充滿(mǎn)了黏彈性流體(Jeffreys流體), 其本構(gòu)方程[24-26]為

圖1 流體-多孔介質(zhì)雙層系統(tǒng)的示意圖

本文引入了Boussinesq近似, 當(dāng)流體運(yùn)動(dòng)中溫度變化很小時(shí), 其密度變化也很小, 由于運(yùn)動(dòng)由浮力誘導(dǎo)產(chǎn)生, 因此不考慮浮力時(shí)可略去密度變化.將密度函數(shù)ρ(T)進(jìn)行Taylor展開(kāi), 且保留至一階項(xiàng), 即滿(mǎn)足

ρ=ρ0[1-β(T-T0)],

(2)

其中ρ和ρ0分別是溫度為T(mén)和T0時(shí)的密度,β為常體積膨脹系數(shù).

1.1 控制方程

其中p為壓力,g為重力加速度,T為溫度,κf為流體層中的熱擴(kuò)散率.

其中um為多孔層的速度,φ為孔隙度,pm為多孔層的壓力,K為滲透率,Tm為多孔層的溫度,κm為多孔層的熱擴(kuò)散率,Gm=[φ(ρ0cp)f+(1-φ)(ρ0cp)m]/(ρ0cp)f為多孔層與流體層的比熱容比.

基本狀態(tài)方程為

(9)

為研究系統(tǒng)的穩(wěn)定性, 對(duì)系統(tǒng)施加一個(gè)小的擾動(dòng), 于是在基本狀態(tài)方程(9)中加入一個(gè)擾動(dòng)項(xiàng)

(10)

(11)

(12)

將方程(10)～(12)代入方程(3)～(8)中, 可得流體層和多孔層的線(xiàn)性化控制方程.為便于觀察, 省略擾動(dòng)的上標(biāo), 對(duì)動(dòng)量方程進(jìn)行兩次旋度運(yùn)算以消去壓力項(xiàng)π和πm.為得到無(wú)量綱方程, 分別在兩層中引入如下尺度： 在流體層中, 長(zhǎng)度的尺度為d, 速度的尺度為ν/d, 時(shí)間的尺度為d2/κf, 溫度的尺度為νΔTf/κf； 在多孔層中, 其尺度分別對(duì)應(yīng)其中ν=μ/ρ為運(yùn)動(dòng)黏性系數(shù).

因此, 對(duì)于流體層, 無(wú)量綱控制方程為

對(duì)于多孔層, 無(wú)量綱控制方程為

1.2 邊界條件

在多孔層的底部z=-1處, 邊界條件為等溫邊界條件

wm=0,Tm=0.

(17)

在流體層的頂部z=1處, 具有表面張力的邊界假設(shè)為不可變形且隔熱, 總切向力由剪切力補(bǔ)償, 于是可得

w=0,DT+BiT=0,

(18)

在交界面z=0處, 速度、 溫度和熱通量分別為

(19)

切向應(yīng)力邊界條件[27]為

(21)

無(wú)量綱參數(shù)關(guān)系為

(22)

2 簡(jiǎn)正模態(tài)

利用簡(jiǎn)正模態(tài)方法, 可將擾動(dòng)量表示為

(23)

(24)

將方程(23)代入方程(13)～(16)中, 控制方程變?yōu)?/p>

簡(jiǎn)正模態(tài)的邊界方程為

(31)

流體層的頂部有3個(gè)邊界條件, 多孔層的底部有2個(gè)邊界條件, 流體層和多孔層的交界面有5個(gè)邊界條件.這是一個(gè)由方程(25)～(28)和邊界條件(29)～(31)組成的10階微分方程, 下面求解其中的Rayleigh數(shù)和Marangoni數(shù).

3 數(shù)值方法

(32)

因此, 問(wèn)題轉(zhuǎn)化成求解廣義特征值問(wèn)題AX=RamBX或AX=MaBX, 其中A和B均為5(N+3)×5(N+3)矩陣,X=(W,Q,Θ,Wm,Θm)T,Q=(D2-a2)W.如求解特征值Ram, 矩陣A和B的表達(dá)式為

(33)

(34)

其中

4 結(jié)果與討論

不同N值的數(shù)值結(jié)果對(duì)比列于表1, 分別選取基函數(shù)個(gè)數(shù)N=5,10,16,25計(jì)算其臨界值以驗(yàn)證該算法的收斂性.結(jié)果表明,N=10后的數(shù)值結(jié)果已基本穩(wěn)定, 且隨N的增大更精確, 同時(shí)考慮處理時(shí)間和效率, 選取N=16進(jìn)行計(jì)算.

表1 不同N值的數(shù)值結(jié)果對(duì)比

4.1 應(yīng)力松弛時(shí)間對(duì)Rayleigh-Marangoni不穩(wěn)定性的影響

不同應(yīng)力松弛時(shí)間的Rayleigh對(duì)流如圖2所示. 由圖2可見(jiàn): 在振蕩對(duì)流下, 較小的應(yīng)力松弛時(shí)間可獲得較高的波數(shù)-Rayleigh數(shù)曲線(xiàn); 長(zhǎng)波和短波區(qū)間分別存在一個(gè)極小值, 即Rayleigh數(shù)曲線(xiàn)呈雙模態(tài); 由于左分支的最小值小于右分支的最小值, 因此系統(tǒng)的不穩(wěn)定性由長(zhǎng)波分支控制, 即Rayleigh對(duì)流在多孔層中開(kāi)始進(jìn)行并由其主導(dǎo)； 振蕩對(duì)流的臨界Rayleigh數(shù)隨應(yīng)力松弛時(shí)間的增大而減小, 表明彈性越大的非牛頓流體會(huì)加劇Rayleigh對(duì)流, 使系統(tǒng)變得更不穩(wěn)定.

圖2 不同應(yīng)力松弛時(shí)間的Rayleigh對(duì)流

不同應(yīng)力松弛時(shí)間的Marangoni對(duì)流如圖3所示. 由圖3可見(jiàn), 振蕩對(duì)流下的臨界Marangoni數(shù)曲線(xiàn)隨應(yīng)力松弛時(shí)間的增加而減小, 每條曲線(xiàn)在全局波數(shù)范圍內(nèi)只有一個(gè)最小值, 即Marangoni對(duì)流下的波數(shù)-Marangoni數(shù)曲線(xiàn)呈單模態(tài), 且當(dāng)Marangoni數(shù)取最小值時(shí)均在短波區(qū)間, 表明Marangoni對(duì)流只發(fā)生在流體層中, 而多孔層相對(duì)較穩(wěn)定.

4.2 應(yīng)變弛豫時(shí)間對(duì)Rayleigh-Marangoni不穩(wěn)定性的影響

不同應(yīng)變弛豫時(shí)間的Rayleigh對(duì)流如圖4所示. 由圖4可見(jiàn): 波數(shù)-Rayleigh數(shù)曲線(xiàn)隨應(yīng)變弛豫時(shí)間的增加而上升, 應(yīng)變弛豫時(shí)間提高了Rayleigh對(duì)流的穩(wěn)定性, 表明黏性越小的非牛頓流體越容易導(dǎo)致系統(tǒng)內(nèi)對(duì)流的不穩(wěn)定性; Rayleigh數(shù)曲線(xiàn)呈雙模態(tài)； Rayleigh對(duì)流不穩(wěn)定性由長(zhǎng)波分支控制, 仍在多孔層中發(fā)生并由其主導(dǎo), 流體層較穩(wěn)定.

圖4 不同應(yīng)變弛豫時(shí)間的Rayleigh對(duì)流

不同應(yīng)變弛豫時(shí)間的Marangoni對(duì)流如圖5所示. 由圖5可見(jiàn), 波數(shù)-Marangoni數(shù)曲線(xiàn)隨應(yīng)變弛豫時(shí)間的增加而上升, Marangoni曲線(xiàn)呈單模態(tài), 每條曲線(xiàn)的Marangoni數(shù)最小值均在短波分支區(qū)間, 即Marangoni對(duì)流不穩(wěn)定性在流體層中出現(xiàn)并由其控制.

圖5 不同應(yīng)變弛豫時(shí)間的Marangoni對(duì)流

4.3 導(dǎo)熱性對(duì)Rayleigh-Marangoni不穩(wěn)定性的影響

不同Biot數(shù)的Marangoni對(duì)流如圖6所示. 由圖6可見(jiàn): 與Biot數(shù)大小無(wú)關(guān), Marangoni曲線(xiàn)呈單模態(tài), 且不穩(wěn)定性均受流體層控制； 臨界Marangoni數(shù)隨Biot數(shù)的增加而增加, 當(dāng)Bi→+∞時(shí), 根據(jù)圖像可推測(cè)Ma→+∞, 即在上表面導(dǎo)熱性良好時(shí), 表面張力效應(yīng)可視為消失.

圖6 不同Biot數(shù)的Marangoni對(duì)流

圖7 不同Biot數(shù)的Rayleigh對(duì)流

4.4 厚度比對(duì)Rayleigh不穩(wěn)定性的影響

圖8 不同Biot數(shù)的Rayleigh對(duì)流

4.5 耦合模式下Rayleigh-Marangoni對(duì)流的不穩(wěn)定性

不同Marangoni數(shù)的Rayleigh-Marangoni對(duì)流耦合模式如圖9所示. 由圖9可見(jiàn): 當(dāng)Marangoni數(shù)減小時(shí), Rayleigh數(shù)曲線(xiàn)呈下降趨勢(shì), 表明Marangoni效應(yīng)減弱會(huì)增強(qiáng)耦合模式下系統(tǒng)Rayleigh對(duì)流的不穩(wěn)定性； 中性曲線(xiàn)呈雙模態(tài), 左分支均低于右分支的最小值； 當(dāng)增大Marangoni數(shù), 短波分支的Rayleigh數(shù)的最小值逐漸變大, 即流體層的Rayleigh對(duì)流會(huì)趨于穩(wěn)定, 因此Rayleigh對(duì)流不穩(wěn)定性均以長(zhǎng)波為主, 即對(duì)流發(fā)生在多孔層中.

圖9 不同Marangoni數(shù)的Rayleigh-Marangoni對(duì)流耦合模式

綜上, 本文研究了一類(lèi)物理參數(shù)對(duì)Rayleigh-Marangoni對(duì)流穩(wěn)定性的影響. 采用MATLAB對(duì)Rayleigh數(shù)和Marangoni數(shù)的求解進(jìn)行數(shù)值模擬, 討論各參數(shù)對(duì)兩類(lèi)對(duì)流穩(wěn)定性的增強(qiáng)或破壞作用. 結(jié)果表明: 黏性越大的流體和系統(tǒng)上表面具有越良好導(dǎo)熱性能均有利于加強(qiáng)對(duì)流的穩(wěn)定性, 而彈性越大的流體會(huì)降低其穩(wěn)定性; 適當(dāng)增加厚度比, Rayleigh對(duì)流不穩(wěn)定性的主導(dǎo)分支將發(fā)生轉(zhuǎn)變; 耦合模式中越大的Marangoni效應(yīng)將促進(jìn)Rayleigh對(duì)流的穩(wěn)定性.