多孔介質(zhì)內(nèi)嵌熱源引起的對流換熱格子Boltzmann模擬

劉子陽，陳姿穎，王存海

(北京科技大學能源與環(huán)境工程學院，北京 100083)

多孔介質(zhì)是由固體骨架和由骨架分隔成大量密集的微小空隙所構(gòu)成的物質(zhì)，是地球生物圈的基石之一.土壤、動物的微細血管網(wǎng)絡(luò)、植物的根莖枝葉等都是天然多孔體，陶瓷、磚瓦、玻璃纖維等人造多孔介質(zhì)材料也參與構(gòu)成了人類社會的方方面面.多孔介質(zhì)內(nèi)的對流換熱過程普遍存在于先進材料制造[1]、電子器件冷卻[2]，以及生物醫(yī)學[3]等相關(guān)技術(shù)與工程應(yīng)用中，因此多孔介質(zhì)內(nèi)對流換熱的研究具有重要意義.然后，由于多孔介質(zhì)結(jié)構(gòu)的復(fù)雜性，多孔介質(zhì)內(nèi)對流換熱的實驗研究難度較大且費時費力，因此數(shù)值模擬方法在多孔介質(zhì)內(nèi)對流換熱問題的研究取得了廣泛關(guān)注：Bayta等[4]采用一種隱式的有限差分法(Finite Difference Method，F(xiàn)DM)求解了達西方程和能量方程，數(shù)值模擬了梯形多孔介質(zhì)內(nèi)的自然對流過程并研究了瑞利數(shù)對溫度場和流場的影響；Chu等[5]采用直接數(shù)值模擬(Direct Numerical Simulation，DNS)方法研究了多孔介質(zhì)中的對流換熱問題并分析了雷諾數(shù)(Reynolds number，Re)的影響；Habbachi等[6]應(yīng)用有限體積法(Finite Volume Method，F(xiàn)VM)分析了填充多孔介質(zhì)的三維方腔中的自然對流，發(fā)現(xiàn)固體與流體的熱導率是影響對流換熱的重要因素.

近年來，從格子氣自動機[7]發(fā)展而來的格子-Boltzmann方法(Lattice Boltzmann Method，LBM)[8-14]是一種利用介觀動力學模型來模擬復(fù)雜宏觀輸運現(xiàn)象的數(shù)值方法.該方法具有清晰的物理意義并且易于并行運算，因此廣受關(guān)注并被成功地應(yīng)用于多孔介質(zhì)對流換熱問題[15-20]的研究.目前，針對多孔介質(zhì)中對流換熱的研究多為壁面加熱的情況[21-23]，但是在電子元件冷卻器[24]、化學反應(yīng)器[25]和熱交換器[26]等諸多工程設(shè)施往往同時涉及多孔介質(zhì)和內(nèi)部熱源.因此，研究多孔介質(zhì)內(nèi)嵌熱源引起的對流換熱問題有助于對相關(guān)工程設(shè)備的優(yōu)化設(shè)計.然而，由于底部加熱帶來的熱流分岔等非線性問題[27]以及內(nèi)部熱壁面之間的相互影響，多孔介質(zhì)內(nèi)嵌熱源引起的對流換熱問題[28 ，29]的求解和分析仍需進一步研究.

本文針對底部內(nèi)嵌熱源的多孔介質(zhì)模型，基于Guo和Zhao[30]發(fā)展的求解流場和溫度場的分布函數(shù)，充分考慮內(nèi)嵌熱源和方腔壁面的邊界條件，建立了多孔介質(zhì)內(nèi)嵌熱源引起的對流換熱的格子-Boltzmann模型.通過實驗結(jié)果對比，驗證了本文數(shù)值模型和代碼的正確性，進一步分析了瑞利數(shù)(Rayleigh number，Ra)、達西數(shù)(Darcy number，Da)和介質(zhì)孔隙率(Medium Porosity，ε)等特征參數(shù)對多孔介質(zhì)內(nèi)流場和溫度場的影響.

1 物理模型與數(shù)值方法

1.1 物理模型

考慮如圖1所示的底部存在內(nèi)嵌熱源的多孔介質(zhì).方腔的尺寸為L×H，左右壁面為T=Tc的定溫壁面，上下壁面為絕熱壁面.內(nèi)嵌熱源位于底邊中心位置，尺寸為l×h，溫度恒定為T=Th(Th>Tc).考慮到方腔內(nèi)的流體流速較小，該問題可看作是不可壓縮的層流問題.

圖1 物理模型

1.2 控制方程

在表征單元體系(Representative Elementary Volume，REV)尺度下，多孔介質(zhì)對流換熱問題的宏觀控制方程表示為[31-33]

?·u=0，

(1)

(2)

(3)

公式中：u為流體速度；T為溫度；p為壓力；ε為多孔介質(zhì)的孔隙率；ρ為流體密度；υe為有效運動粘度；αe為有效熱擴散系數(shù)；σ為熱容比，此處為1；F為多孔介質(zhì)和其他外力產(chǎn)生的總作用力，可表示為[32]

(4)

G=gβ(T-T0)j，

(5)

公式中：G為浮力，υ為運動粘度，g為重力加速度，β表示熱膨脹系數(shù)，T0為基準溫度，j為垂直方向的單位矢量；Fε為幾何形狀因子，K為滲透率，根據(jù)Ergun的經(jīng)驗關(guān)系[34]可以表示為[35]

(6)

(7)

對于上述控制方程所描述的物理問題，可以由無量綱參數(shù)表征，即達西數(shù)Da、瑞利數(shù)Ra、普朗特數(shù)Pr、努賽爾數(shù)Nu與粘度比J為

(8)

(9)

同時，無量綱速度和溫度可以表示為

(10)

1.3 格子-Boltzmann模型

本文采用多分布函數(shù)的熱格子-Boltzmann模型，邊界條件采用非平衡外推格式.流場部分使用含外力項的D2Q9格子的單松弛模型處理，溫度場部分使用D2Q5格子的單松弛模型處理，其演化方程與對應(yīng)的平衡分布函數(shù)為[30]

(11)

(12)

(13)

(14)

公式中：Fi為外部體力項，可表示為

(15)

宏觀變量可表示為

(16)

(17)

公式中：引入的中間速度v與修正參數(shù)c0和c1為

(18)

(19)

1.4 邊界條件

對于流場，采用非平衡外推邊界，即將邊界節(jié)點處的分解為平衡部分與非平衡部分，邊界節(jié)點的離散分布函數(shù)fi(xs，t)表示為[36]

(20)

公式中：xf=xs+eiδt表示和固體邊界節(jié)點xs相鄰的流體節(jié)點xf.

溫度場的邊界條表示為[37]

左邊界

(21)

右邊界

(22)

公式中：對于外部方腔的左右邊界，無量綱溫度取θ=θc；對于內(nèi)嵌熱源邊界取θ=θh.內(nèi)嵌熱源頂壁的邊界條件為

(23)

封閉方腔的上下邊界是絕熱的，故可以根據(jù)下式求取未知的分布函數(shù).

上邊界

g4(x，ym)=g4(x，ym-1)，

(24)

下邊界

g2(x，y0)=g2(x,y1)，

(25)

公式中：下標m和0分別表示方腔頂壁和底壁在y方向的節(jié)點索引.

2 模型驗證

2.1 網(wǎng)格無關(guān)性驗證

首先設(shè)置孔隙率和達西數(shù)分別為ε=0.999 9和Da=108，將本文發(fā)展的多孔介質(zhì)模型應(yīng)用于底邊局部加熱自然對流過程[38]的求解：加熱尺寸為d=0.8 L，模擬參數(shù)設(shè)置為Pr=0.71、Ra=1.836×105.計算所得的熱壁面平均努賽爾數(shù)Nuavg如圖2所示，模擬選用的網(wǎng)格類型分別有Nx×Ny=40×40，60×60，80×80，100×100，120×120，140×140和160×160.從圖中可以看出，當網(wǎng)格數(shù)量達到120×120時，努賽爾數(shù)的變化逐漸趨于平穩(wěn).為確保收斂，Nx×Ny=300×300時的工況也被計算，其數(shù)值結(jié)果為Nuavg=7.501，以該結(jié)果為基準，采用120×120網(wǎng)格所得結(jié)果的相對誤差小于1.0%.綜合考慮計算效率和精度，本文接下來的數(shù)值模擬所選取的網(wǎng)格劃分均為120×120.

圖2 網(wǎng)格無關(guān)性驗證

2.2 正確性驗證

首先將采用本文模型得到的模擬結(jié)果與文獻[38]的實驗圖像進行比對.由圖3可知，本文采用的LB模型精確地復(fù)現(xiàn)了該方腔內(nèi)的溫度分布特征，驗證了本文模型的準確性.

圖3 底部局部加熱方腔內(nèi)溫度分布

為定量分析模型的數(shù)值可靠性，本文模擬了不同瑞利數(shù)時底邊加熱的自然對流，分析高溫壁面上的平均努賽爾數(shù)隨瑞利數(shù)的變化.表1中列出了在加熱尺寸分別為d=0.2L和d=0.8L時Ra=103、104、105和106的數(shù)值結(jié)果，與文獻[38]中的數(shù)據(jù)符合程度良好.以上對比結(jié)果表明本文采用的LB模型能夠準確地模擬自然對流換熱現(xiàn)象.

表1 本文LB得到的加熱壁面上的Nuavg和實驗結(jié)果[38]對比

3 結(jié)果與討論

3.1 瑞利數(shù)變化對換熱的影響

瑞利數(shù)是與浮升力直接相關(guān)的無量綱準則數(shù).為了研究瑞利數(shù)對具有底部矩形熱源的多孔介質(zhì)方腔內(nèi)部自然對流的影響，本節(jié)將模擬瑞利數(shù)為Ra=101～107時的溫度場至穩(wěn)態(tài).相關(guān)模擬參數(shù)設(shè)置為Da=10-3，ε=0.4，Pr=0.71，l=0.2L，h=0.2H，部分計算結(jié)果如圖4所示.

圖4 不同瑞利數(shù)條件下多孔介質(zhì)內(nèi)的等溫線(左)和流線(右)

由圖4(左)所示的溫度分布結(jié)果可知，熱量從熱源壁面逐漸擴散到多孔介質(zhì)的內(nèi)部空間，隨著Ra的增加，方腔頂部溫度梯度更大.觀察圖4(a)可以發(fā)現(xiàn)，當Ra較低時，等溫線在熱源壁面附近密集，此時的換熱性質(zhì)主要為熱傳導；隨著Ra的增加，熱傳導的作用減弱，熱量逐漸從熱源壁面附近向多孔方腔的頂部擴散：在Ra=106時，對流換熱效果十分顯著.圖4(右)為對應(yīng)圖4(左)工況下的流線圖.根據(jù)流線分布可以發(fā)現(xiàn)，多孔方腔中形成了沿x=0.5L對稱的反向旋渦：當Ra=104時，旋渦的中心位置較低，同時流線主要在方腔底部密集；當Ra增大時，旋渦中心位置隨之升高，在方腔頂部能夠觀察到更為密集的流線分布：多孔介質(zhì)內(nèi)的整體流動情況得到明顯改善.

不同Ra下的垂直速度分布曲線與Nuavg變化曲線如圖5所示.根據(jù)圖5(a)展示的多孔介質(zhì)腔體在x=0～0.5L、y=0.5H處的垂直速度分布可知：在低Ra時，速度幾乎為零，流動處于停滯狀態(tài)；而Ra=105～107時，由于浮升驅(qū)動力的增強，速度曲線產(chǎn)生明顯的變化：最大速度出現(xiàn)在熱源上方位置，最小速度則靠近低溫側(cè)壁，流體的流動行為較為顯著.由此可見，增大Ra能夠有效促進多孔介質(zhì)中流體的運動，但是需要達到臨界值：在本文的研究工況和數(shù)值范圍內(nèi)，該臨界值為Ra=105.從圖5(b)中可以發(fā)現(xiàn)，熱源各壁面的Nuavg與Ra呈現(xiàn)正相關(guān)性，即Ra的增大能夠提升多孔介質(zhì)的整體對流換熱強度.但考慮到高溫熱源存在多個熱壁面，且壁面間的熱交換可能相互影響，故還需要對熱源的各個壁面進行對比分析.

圖5 不同瑞利數(shù)條件下：(a)介質(zhì)中心線y/H=0.5上流體垂直速度分布；(b) 熱源壁面上的Nuavg

當Ra較低，例如Ra=101～104時，沿各壁面處的Nuavg均低于1.0，說明多孔介質(zhì)內(nèi)的換熱形式以熱傳導為主.同時可以發(fā)現(xiàn)，沿熱源側(cè)壁(含熱源左側(cè)壁面與右側(cè)壁面)處的Nuavg始終低于沿頂壁處的Nuavg.隨著Ra的提升，沿熱源側(cè)壁處的Nuavg呈現(xiàn)增長趨勢，沿頂壁處的Nuavg則逐漸降低，直到Ra=104時兩者幾乎相等.這意味著，在換熱性質(zhì)主要為熱傳導時，頂壁的換熱情況并未隨著浮升力的增強得到改善，反而是隨側(cè)壁換熱能力的提升而出現(xiàn)減弱.由此可見，熱源頂壁處的換熱實際上受到了側(cè)壁換熱的影響：隨著浮升驅(qū)動力的增大，熱源側(cè)壁向上方傳遞的熱量增多.但由于熱傳導的影響，熱源頂壁附近熱邊界層較厚，熱流體處于流動滯止區(qū)，頂壁處的溫度梯度低，從而使Nuavg的數(shù)值下降.

當Ra繼續(xù)增加，Nuavg出現(xiàn)明顯提升，同時沿熱源側(cè)壁處的Nuavg開始高于頂壁，這在Ra=107時尤為顯著.由于Nuavg遠高于1.0，熱對流取代熱傳導成為了多孔介質(zhì)內(nèi)換熱的主要形式.觀察圖像還可以發(fā)現(xiàn)，在臨界值Ra=105處，沿頂壁處的Nuavg繼續(xù)降低，而此時沿側(cè)壁處的換熱性質(zhì)開始向?qū)α鬓D(zhuǎn)變.同時在Ra>105時，沿側(cè)壁處的Nuavg隨Ra增加而大幅提升；熱源頂壁處的換熱情況雖然得到改善，但其Nuavg仍遠低于側(cè)壁.即在高Ra下，熱源側(cè)壁處的換熱是影響多孔介質(zhì)內(nèi)對流換熱強度的重要因素，頂壁處的對流換熱仍受到側(cè)壁換熱的影響，導致Nuavg仍處于較低水平.綜上，對于本文討論的多孔介質(zhì)方腔對流換熱問題，熱壁面間的相互作用是需要考慮的關(guān)鍵因素.該研究能夠為針對性地改善多孔介質(zhì)內(nèi)部換熱提供相關(guān)的理論指導.

3.2 達西數(shù)變化對換熱的影響

達西數(shù)是表征多孔介質(zhì)滲透能力的無量綱準則數(shù).在本文的工作中，討論了Da=10-7～10-1時的多孔介質(zhì)對流換熱，其它模擬參數(shù)設(shè)置為Ra=106，ε=0.4，Pr=0.71，l=0.2L，h=0.2H，模擬溫度場至穩(wěn)態(tài)的部分結(jié)果如圖6所示.

圖6 不同達西數(shù)條件下多孔介質(zhì)內(nèi)的等溫線(左)和流線(右)：(a) Da=10-4；(b) Da=10-3；(c) Da=10-2

根據(jù)模擬結(jié)果可知，Da較低時，流體的流動強度受限于多孔介質(zhì)的滲透能力，導致對流換熱強度處于較低水平；隨著Da的增加，多孔介質(zhì)的滲透能力逐漸增強，流體的流動強度逐漸增大.在Da=10-2時，由于多孔介質(zhì)對流體部分的影響較小，溫度與流線分布中存在顯著的對流效應(yīng)，整體流動換熱情況與純流體類似.故可以發(fā)現(xiàn)Da對多孔介質(zhì)對流換熱的影響與Ra類似，是決定多孔介質(zhì)內(nèi)部換熱形式的主要因素.

Da的變化對垂直速度分布和Nuavg的影響如圖7所示.根據(jù)速度曲線的變化可以發(fā)現(xiàn)，Da與Ra類似，它的增加能夠強化流體的運動.但在Da>10-4后，隨著Da的繼續(xù)增加，方腔中的流體運動得到了整體增強，而Ra的增加則僅會使方腔的x/L=0與x/L=0.5附近出現(xiàn)顯著的速度梯度.由圖7(b)可知，熱源各壁面的Nuavg隨Da的增加而出現(xiàn)提升：在低Da的情況下，熱源頂壁的換熱情況優(yōu)于側(cè)壁，此時的換熱性質(zhì)主要為熱傳導.隨著Da的增加，熱源頂壁處的Nuavg逐漸降低，同時在Da>10-5時側(cè)壁處的Nuavg開始高于頂壁處的Nuavg；當換熱形式以熱對流為主時，頂壁的換熱能力開始與Da的提升成正比.結(jié)合上節(jié)的分析可以發(fā)現(xiàn)，在換熱形式為熱對流時，熱源頂壁處的換熱主要取決于Da和Ra.而在熱傳導時，熱源側(cè)壁的熱交換成為了影響頂壁換熱情況的主要因素.

圖7 不同達西數(shù)條件下：(a)介質(zhì)中心線y/H=0.5上流體垂直速度分布；(b)熱源壁面上的Nuavg

3.3 孔隙率變化對換熱的影響

孔隙率ε是由于孔隙的存在而引申出的密實固體所不曾具備的、多孔介質(zhì)所特有的屬性.為了研究不同的孔隙率對多孔介質(zhì)內(nèi)對流換熱的影響，選取孔隙率為ε=0.1 ～ 0.9，同時設(shè)置對流換熱參數(shù)為Da=10-3，Ra=106.其余模擬參數(shù)設(shè)置為Pr=0.71，l=0.2L，h=0.2H，模擬溫度場至穩(wěn)態(tài)的部分結(jié)果如圖8所示.

圖8 不同孔隙率條件下多孔介質(zhì)內(nèi)的等溫線(左)和流線(右)：(a)ε=0.2；(b)ε=0.4；(c)ε=0.6

從圖8中可以發(fā)現(xiàn)，孔隙率ε的變化對溫度場的影響相對較小，但在孔隙率較高時，能夠觀察到溫度更為有力的擴散現(xiàn)象，方腔頂部具有面積更大的高溫區(qū)域.從圖8的流線分布來看，多孔介質(zhì)中的流體運動得到增強，但方腔內(nèi)對稱反向旋渦的形狀與旋渦的中心位置幾乎一致，整體的流動性質(zhì)未能產(chǎn)生明顯變化.由此可見，孔隙率的變化對于流場不具有顯著影響，這與Ra和Da的情況不同.

垂直速度分布和Nuavg曲線隨孔隙率的變化如圖9所示.從圖9中可以發(fā)現(xiàn)，隨著孔隙率的增加，多孔介質(zhì)中的流體運動得到顯著強化，同時沿熱源各個壁面的Nuavg出現(xiàn)明顯提升，即孔隙率的增加能夠有效改善多孔介質(zhì)內(nèi)對流換熱效果.特別的是，由于在孔隙率較低時仍能夠觀察到明顯的對流效應(yīng)，這說明孔隙率的變化未能使多孔介質(zhì)內(nèi)的換熱性質(zhì)發(fā)生轉(zhuǎn)變.故孔隙率并不是影響多孔介質(zhì)內(nèi)部的流動換熱情況的主導因素.

圖9 不同孔隙率條件下：(a)介質(zhì)中心線y/H=0.5上流體垂直速度分布；(b)熱源壁面上的Nuavg

4 結(jié) 論

本文采用LBM模擬底邊中心處存在矩形熱源的多孔介質(zhì)方腔內(nèi)的自然對流現(xiàn)象，討論了不同的Ra、Da以及孔隙率ε對多孔介質(zhì)對流換熱特性的影響.研究結(jié)果表明：

1.當Ra>105時，多孔介質(zhì)內(nèi)對流換熱水平隨Ra增加而得到顯著提升；當Ra≤105時，Ra對換熱的影響可以忽略.

2.當Da>10-5時，Da的增加能有效改善多孔介質(zhì)內(nèi)的流動換熱情況；當Da≤10-5時，Nuavg不隨Da而產(chǎn)生明顯變化.

3.孔隙率的增加能夠強化流動換熱，但不會改變多孔介質(zhì)內(nèi)換熱性質(zhì).

4.當換熱性質(zhì)以熱傳導為主時，隨Ra或Da的增加，沿熱源側(cè)壁處的Nuavg逐漸提升，并使熱源頂壁處的熱邊界層變厚，導致沿頂壁處的Nuavg小幅降低.