舒曉暉，穆朗楓，牛 駿，丁 健

(1.中國石油化工股份有限公司，北京100083;2.密蘇里科技大學(xué)，美國羅拉密蘇里州65401;3.中國石油大學(xué)，北京102249;4.中國海油 伊拉克有限公司，北京100010)

粘性和非線性的存在使得多孔介質(zhì)中流體熱對流研究極為復(fù)雜，但是相關(guān)研究成果在石油工程、生物工程等方面均有著廣泛的應(yīng)用前景［1－4］。因此，目前國內(nèi)外有大量學(xué)者正在該領(lǐng)域開展研究工作?，F(xiàn)有成果對于多孔介質(zhì)內(nèi)牛頓流體的流動(dòng)與傳熱已經(jīng)進(jìn)行了廣泛深入分析［5－9］，例如 Kubitschek［7］等通過引入畢奧數(shù)對3類底部加熱邊界條件下多孔介質(zhì)中牛頓流體的熱對流啟動(dòng)進(jìn)行了深入研究，并總結(jié)了各參數(shù)對于熱對流啟動(dòng)點(diǎn)的影響。但是對于非牛頓流體在多孔介質(zhì)內(nèi)的流動(dòng)與傳熱研究還相對匱乏。

Kim等［10］對上、下兩邊均處于封閉且等溫加熱邊界條件下的多孔介質(zhì)內(nèi)粘彈性流體熱對流進(jìn)行了線性和非線性穩(wěn)定性分析，結(jié)果熱對流啟動(dòng)時(shí)總是振蕩對流，振蕩對流存在超臨界和次臨界分叉現(xiàn)象。Jun等［11］對封閉邊界多孔介質(zhì)內(nèi)黏彈性流體在不同底部加熱邊界條件下的熱對流啟動(dòng)進(jìn)行了線性穩(wěn)定性分析，結(jié)果等熱流加熱邊界條件下熱對流啟動(dòng)最早。但目前極少有針對非封閉多孔介質(zhì)內(nèi)黏彈性流體熱對流啟動(dòng)特征的研究。

文中通過引入修正后的Darcy定律，對上端開口、底部等溫加熱條件下多孔介質(zhì)內(nèi)粘彈性流體熱對流啟動(dòng)進(jìn)行了線性穩(wěn)定新分析，得出了熱對流啟動(dòng)時(shí)的臨界條件及黏彈性特征對于啟動(dòng)點(diǎn)的影響。

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

考慮一個(gè)高度為d，長度為a的二維矩形多孔介質(zhì)。垂直兩側(cè)面絕熱，并且流體不能透過。多孔介質(zhì)上下邊界維持恒溫T1，T2，且T2－T1=ΔT。上邊界開口，并與外界連通，壓強(qiáng)為p0。假設(shè)多孔介質(zhì)的滲透率為K，并且其中充滿了不可壓粘彈性流體，該流體動(dòng)力學(xué)粘度為μ，熱膨脹系數(shù)為β，密度為ρ，整個(gè)系統(tǒng)的熱擴(kuò)散系數(shù)為κ。通過運(yùn)用修正的Darcy定律和Boussinesq近似，可以得到粘彈性流體在多孔介質(zhì)中運(yùn)動(dòng)的控制方程組:

該模型的邊界條件為:

其中邊條件(5b)為等壓邊界條件，用于描述多孔介質(zhì)上端開口邊界。將其帶入(2)得

為了方便計(jì)算，對上述方程組中的變量進(jìn)行了無量綱化。這里選擇的無量綱化參量為:長度d，時(shí)間d2/κ，速度 κ/d，壓強(qiáng) κμ/K，溫度 ΔT。

2 線性穩(wěn)定性分析

將控制方程組和邊界條件無量綱化后，引入溫度擾動(dòng)量θ和垂直速度擾動(dòng)量w得:

其中，λ，ε為無量綱化后的松弛時(shí)間和延遲時(shí)間，Ra=Kρ1βgdΔT/μκ 為瑞利數(shù)。邊界條件為:

根據(jù)流動(dòng)模式分析，溫度、速度擾動(dòng)量在水平方向上均存在周期性:

其中，α=mπ/a，m為正整數(shù)，Θ(z)為關(guān)于z的函數(shù)。

將方程(16)代入(12)得垂直速度擾動(dòng)量

式中:D=?/?z為微分算子。

將方程(16)，(17)代入方程(12)—方程(14)得關(guān)于z的常微分方程:

分析該微分方程形式易知其通解為:

這里，Aj(j=1，2，3，4)是待定常數(shù)。將(19)代回原方程(18)得kj滿足方程:

由邊界條件易得A2=－A1，A4=－A3。

由于方程(20)中僅含有k的偶數(shù)次冪，因此k的4個(gè)解滿足k2=－k1，k4=－k3。進(jìn)一步求解(20)可得:

由邊界條件式(13)和式(14)知 A1，A3和 k1，k3滿足下列方程組:

當(dāng)多孔介質(zhì)內(nèi)開始發(fā)生對流時(shí)，內(nèi)部流體的速度分布不能處處為零。所以線性方程組(17)和(18)必須有非零解，即方程組系數(shù)矩陣的行列式必須為零:

上式中當(dāng)σ=0時(shí)代表的是穩(wěn)態(tài)對流。而對于粘彈性流體，熱對流發(fā)生時(shí)為振蕩對流，此時(shí)σ應(yīng)該是一個(gè)純虛數(shù)iω，這里ω定義為振蕩頻率。在此情況下，kj(j=1，2，3，4)與 σ 耦合在一起，需要同時(shí)進(jìn)行求解，因此無法得出解析解形式。我們使用了數(shù)值方法對振蕩對流情況進(jìn)行編程求解。

3 結(jié)果與討論

為了研究λ，ε對多孔介質(zhì)中的粘彈性流體對流產(chǎn)生和流動(dòng)模式的影響，選取下面四組粘彈性參數(shù)值作為算例:(1)λ=0.003，ε=0.001;(2)λ =0.3，ε=0.2;(3)λ =0.2，ε =0.1;(4)λ =0.3，ε =0.1。

3.1 粘彈性對熱對流啟動(dòng)的影響

算例(1)中的松弛時(shí)間和延遲時(shí)間值均很小，其熱對流啟動(dòng)時(shí)的臨界瑞利數(shù)Rac隨多孔介質(zhì)邊長變化如圖1中所示。隨著邊長的增加，Rac值振蕩衰減。其全局最小瑞利數(shù)值為29.298，該值與前人得出的開口多孔介質(zhì)中牛頓流體熱對流啟動(dòng)時(shí)的臨界瑞利數(shù)29.3非常接近。由于牛頓流體熱對流啟動(dòng)時(shí)總是穩(wěn)態(tài)對流，因此可以推斷當(dāng)粘彈性參數(shù)都很小時(shí)，微弱的粘彈性效應(yīng)無法引起振蕩對流，首先出現(xiàn)的是穩(wěn)態(tài)對流。

圖1 算例(1)中臨界瑞利數(shù)隨多孔介質(zhì)邊長變化Fig.1 Variation of the critical Rayleigh number as a function of the medium geometry for case(1)

圖2 算例(2)，(3)和(4)中臨界瑞利數(shù)隨多孔介質(zhì)邊長變化Fig.2 Variation of the critical Rayleigh number as a function of the medium geometry for cases(2)，(3)and(4)respectively

圖2a—c顯示算例(2)，(3)和(4)對應(yīng)的臨界瑞利數(shù)隨多孔介質(zhì)邊長的變化，圖中也標(biāo)出了對應(yīng)不同階段最佳模態(tài)m的取值。3個(gè)算例的全局最小臨界瑞利數(shù)為28.953，27.567和19.112，均小于穩(wěn)態(tài)對流啟動(dòng)時(shí)的臨界瑞利數(shù)29.298，計(jì)算結(jié)果也表明對流啟動(dòng)時(shí)均為振蕩對流。由此得出，當(dāng)流體的粘彈性參數(shù)足夠大時(shí)，振蕩對流優(yōu)先啟動(dòng)，并且其啟動(dòng)點(diǎn)早于穩(wěn)態(tài)對流。此外，從圖2中還可以發(fā)現(xiàn)，隨著多孔介質(zhì)長度的增加，速度場和溫度分布中的最佳模態(tài)在不斷增大。

算例(2)和(4)具有相同的松弛時(shí)間0.3s，延遲時(shí)間不同。對比兩例的計(jì)算結(jié)果，圖2a和圖2c顯示，對每一個(gè)固定邊長a，臨界瑞利數(shù)Rac隨延遲時(shí)間ε值的增大而增大，說明較大的延遲時(shí)間能夠推遲振蕩對流的發(fā)生。兩幅子圖中Rac的兩個(gè)相鄰極小值之間的距離相同，這個(gè)距離隨著ε值的減小而縮短。

算例(3)和(4)中的延遲時(shí)間相同，均為0.1，松弛時(shí)間不同，對比圖2b和圖2c，對于任意一個(gè)固定的邊長a，Rac值隨著λ值的增大而減小。因此，較大的松弛時(shí)間值能夠促進(jìn)粘彈性流體在多孔介質(zhì)內(nèi)振蕩對流的發(fā)生。上述兩個(gè)算例對應(yīng)的臨界瑞利數(shù)值不盡相同，但是兩條曲線的變化趨勢非常相似，且對應(yīng)的兩幅圖中的兩個(gè)相鄰極值間的距離一致。由此可得出結(jié)論，此時(shí)松弛時(shí)間λ的取值對于最優(yōu)模態(tài)的變化沒有影響。

3.2 粘彈性對溫度分布的影響

取邊長為3.2，時(shí)間 t=2nπ/ω，圖 3 a，b 和 c分別為算例(1)，(2)和(4)中熱對流啟動(dòng)時(shí)的溫度擾動(dòng)量等值線圖，其流動(dòng)模式分別為m=2，3，4。其中亮斑表示溫度較高區(qū)域，暗斑表示溫度較低區(qū)域?？梢钥闯?，亮斑與暗斑交替出現(xiàn)，亮、暗斑總數(shù)與流動(dòng)模式有關(guān)。在上下兩邊界處的溫度等值線均與邊界平行，這與上下兩邊的等溫邊界條件一致。算例(2)，(4)對應(yīng)的振蕩對流溫度等值線圖明顯比算例(1)對應(yīng)的穩(wěn)態(tài)對流溫度等值線圖復(fù)雜。由于松弛時(shí)間λ的取值對于最優(yōu)模態(tài)的變化沒有影響，因此算例(3)對應(yīng)的熱對流啟動(dòng)時(shí)的溫度等值線圖形式應(yīng)與算例(4)一致。

3.3 與閉口多孔介質(zhì)內(nèi)計(jì)算結(jié)果的對比

圖4中顯示的是在開口和閉口邊條件下對應(yīng)算例(2)的臨界瑞利數(shù)隨多孔介質(zhì)邊長變化曲線，其中閉口情況對于曲線一直位于開口情況對應(yīng)曲線上方，二者的全局最小臨界瑞利數(shù)分別為28.953和32.643。余下3個(gè)算例(1)，(3)和(4)在開口和閉口條件下的臨界瑞利數(shù)計(jì)算結(jié)果也呈現(xiàn)出相同的規(guī)律。上述結(jié)果說明開口邊界條件對于多孔介質(zhì)中的流動(dòng)起到擾動(dòng)作用，有利于熱對流的啟動(dòng)。

圖3 算例(1)，(2)和(4)中，取邊長為3.2、時(shí)間為2nπ/ω時(shí)熱對流啟動(dòng)時(shí)的溫度擾動(dòng)量等值線Fig.3 Contours of the disturbed temperature at the onset point for cases(1)，(2)and(4)respectively，with the length of the medium being 3.2 and the time being 2nπ/ω

圖4 算例(2)對應(yīng)開口和閉口邊界條件時(shí)的臨界瑞利數(shù)隨多孔介質(zhì)邊長變化Fig.4 Variation of the critical Rayleigh number with the length of porous media corresponding the open-top boundary conditions and closed-top boundary conditions for case(2)

4 結(jié)論

1)對于上端開口底部受等溫加熱多孔介質(zhì)內(nèi)的粘彈性流體，當(dāng)粘彈性參數(shù)都很小時(shí)微弱的粘彈性效應(yīng)無法引起振蕩對流，首先出現(xiàn)的是穩(wěn)態(tài)對流。

2)當(dāng)流體的黏彈性參數(shù)足夠大時(shí)，振蕩對流優(yōu)先啟動(dòng)，且振蕩對流啟動(dòng)時(shí)的臨界瑞利數(shù)隨多孔介質(zhì)邊長增大而呈振蕩衰減。

3)松弛時(shí)間的增大及延遲時(shí)間的減小能夠使得粘彈性流體在多孔介質(zhì)內(nèi)的振蕩對流提前啟動(dòng)。但是松弛時(shí)間的取值對于最優(yōu)模態(tài)的變化沒有影響。

4)相對于閉口邊界條件，開口邊界條件對多孔介質(zhì)內(nèi)流動(dòng)具有更強(qiáng)的擾動(dòng)作用，能夠提前熱對流啟動(dòng)。

符號注釋

d:多孔介質(zhì)高度，m;a:多孔介質(zhì)長度，m;T1:多孔介質(zhì)上邊界溫度，K;T2:多孔介質(zhì)下邊界溫度，K;p0:多孔介質(zhì)上邊界壓力，Pa;K:多孔介質(zhì)滲透率，m2;μ:粘彈性流體粘度，Pa·s;β:粘彈性流體熱膨脹系數(shù)，1/K;ρ:粘彈性流體密度，kg/m3;κ:多孔介質(zhì)內(nèi)粘彈性流體熱擴(kuò)散系數(shù)多孔介質(zhì)內(nèi)的速度矢量，m/s;u*:水平速度分量，m/s;v*:垂直速度分量，m/s;p*:粘彈性流體壓力，Pa;g:重力加速度，m/s2:延遲時(shí)間，s:松弛時(shí)間，s;:垂直方向單位矢量;ρ1:溫度在T1時(shí)的流體密度，kg/m3;θ:無量綱溫度擾動(dòng)量;w:無量綱垂直速度擾動(dòng)量;ε:無量綱延遲時(shí)間;λ:無量綱松弛時(shí)間;Ra:瑞利數(shù);m:任意整數(shù);Θ(z):關(guān)于z的函數(shù);D:微分算子;ω:振蕩頻率。

［1］王補(bǔ)宣.多孔介質(zhì)傳熱研究的意義和現(xiàn)狀［J］.中國科學(xué)基金，1991，1:32 －34.Wang Buxuan.Meaning and of research on heat transfer in porous media［J］.China Science Foundation，1991，1:32 －34.

［2］劉偉.多孔介質(zhì)傳熱傳質(zhì)理論與應(yīng)用［M］.科學(xué)出版社，2006.Liu Wei.Theory and application and heat and mass transfer in porous media［M］.Science Press，2006.

［3］Khaled A R A，Vafai K.The role of porous media in modeling flow and heat transfer in biological tissues［J］.International Journal of Heat and Mass Transfer，2003，46(26):4989 －5003.

［4］Kuznetsov A V，Avramenko A A.Stability analysis of bioconvection of gyrotactic motile microorganisms in a fluid saturated porous media［J］.Transport in Porous Media，2003，53(1):95 －104.

［5］Horton C W，Rogers F T.Convection current in porous medium heated from below［J］.Journal of Applied Physics，1945，6:367 －370.

［6］Wang C Y.Onset of convection in a fluid-saturated rectangular box bottom heated by constant flux［J］.Physics of Fluids，1999，11:1673 －1675.

［7］Kubitschek J P，Weidman P D.Stability of a fluid － saturated porous medium heated from below by forced convection［J］.International Journal of Heat and Mass Transfer，2003，46(9):3697 －3705.

［8］Hayat T，Javed T，Abbas Z.Slip flow and heat transfer of a second grade fluid past a stretching sheet through a porous space［J］.International Journal of Heat and Mass Transfer，2008，51(17):4528－4534.

［9］Fetecau C，Hayat T，F(xiàn)etecau C.Starting solutions for oscillating motions of Oldroyd-B fluids in cylindrical domains［J］.Journal of Non-Newtonian Fluid Mechanics，2006，153(2):191 －201.

［10］Kim MC，et al.Thermal instability of viscoelastic fluids in porous media［J］.International Journal of Heat and Mass Transfer，2003，46(26):5065－5072.