含內熱源的多孔方腔自然對流非正交MRT-LBM 模擬

張瑩，包進，過海龍，連小龍，黃逸宸，李培生

（南昌大學 機電工程學院，南昌330031）

多孔介質中的流體流動與傳熱在科學和工程領域有許多應用，如土木工程、地質力學、環境科學、生物工程、材料科學及再生熱交換器、食品加工、電子設備冷卻等其他相關領域［1-2］。因此，國內外許多學者對多孔腔體內流體流動及傳熱特性進行了深入研究，并取得了重要成果。邱偉國等［3］采用有限元方法對左側部分多孔介質壁面方形封閉腔體內的自然對流傳熱進行了數值模擬。Yaacob和Hasan［4］采用有限差分的方法研究了傾斜多孔矩形方腔內的自然對流現象。

格子Boltzmann方法（LBM）作為近幾年發展起來的一種介觀數值模擬方法，與傳統的宏觀數值方法相比具有物理圖像清晰、邊界條件容易處理及并行性能好等優點［5-6］。因此，國內外許多學者采用LBM方法對多孔介質中的傳熱現象進行了數值模擬。Shokouhmand等［7］使用LBM 方法模擬了2個平行板之間的層流流動與對流傳熱問題。陸 威 等［8］基 于CLBGK （Coupled Lattice Bhatnager-Gross-Krook）模型，利用LBM 方法研究了頂蓋驅動復合方腔內的雙擴散混合對流現象，研究結果表明，多孔介質層對頂蓋驅動復合方腔內的熱質雙擴散的影響十分顯著。Huelsz和Rechtman［9］利用LBM 方法對傾斜方腔內的自然對流傳熱問題進行了模擬研究，分析了方腔傾斜角度及Rayleigh數對其流動與傳熱的影響。Zhao等［10］采用LBM方法研究了孔密度（孔尺寸）和孔隙率對自然對流的影響，結果表明，通過降低孔隙率和孔尺寸可以提高方腔內的整體傳熱效率。Bouarnouna等［11］采用多松弛格子Boltzmann方法（MRT-LBM）模擬多孔介質通道中的對流傳熱現象，表明Darcy數、Rayleigh數、多孔基材的厚度、導流板的厚度及位置對通道內的傳熱有著重要的影響。李培生等［12］采用非正交MRT方法對傾斜多孔方腔內自然對流進行數值研究，結果表明，方腔傾角對腔內傳熱影響很大，熱壁面上平均Nusselt數隨傾角的變化表現出M 型分布規律。

針對多孔方腔內流體流動與傳熱的研究，大多采用連續性的冷熱壁面條件［13-15］，然而在電子器件的冷卻和建筑物的冷卻等領域間斷式的冷卻方式有著廣泛的應用。基于此，本文主要研究了含有內熱源的多孔方腔的自然對流傳熱問題，對比了6種不同的冷源布置方案、內熱源結構形式及內熱源的位置對方腔內自然對流傳熱的影響，討論了Darcy數、Rayleigh數等參數對方腔內流體流動與傳熱特性的影響。

1 物理模型的建立

圖1 物理模型Fig.1 Physicalmodel

圖1為二維多孔方腔內含不同結構的內熱源示意圖，高溫熱塊表面溫度為Th，方腔左右壁面加粗藍線表示低溫壁面，溫度為Tc，其余壁面均為絕熱壁面。方腔內填充多孔介質為各向同性且均勻的材料。二維多孔方腔邊長為L，對應于Case 1方形高溫熱塊的邊長為H，Case 2、Case 3高溫熱塊高度為H，寬度為d，Case 1對應的H/L為0.4，Case 2、Case 3對應的H/L分別為0.3、0.5，d/L則分別為0.5、0.3，即使內置高溫塊的熱邊界周長相等。圖2給出了Case 1對應內熱源布置方式下的6種不同等溫冷源邊界的布置方案。

假設流體為不可壓縮流體，忽略黏性散熱，基于廣義的非Darcy模型，多孔介質中的流體流動和對流傳熱在局部熱平衡條件下的宏觀控制方程為

式中：ρ0為平均流體密度；u、T和p分別為流體平均速度、溫度和壓力；ε為孔隙率；υe為有效運動黏度；σ為孔隙間流體與多孔介質固體骨架之間熱容比值；αe為有效熱擴散系數；F=（Fx，Fy）為由多孔介質和其他外力引起的合外力項，可以表示為

其中：等號右端的第1項、第2項、第3項分別為線性Darcy介質阻力、非線性Forchheimer介質阻力、浮升力；υ為流體的運動黏度；K為多孔介質的滲透率；Fε為幾何函數；G為浮升力。

圖2 六種不同等溫冷源邊界布置方案Fig.2 Six different isothermal cold source boundary arrangements

式中：β為熱膨脹系數；g為重力加速度矢量；T0為參考溫度，本文取為系統的平均溫度；θ為方腔的傾斜角度，為0°；i為x軸的單位矢量，j為y軸的單位矢量。

方腔內置高溫塊熱壁面上的平均Nusselt數Nuave的計算公式為

對于流體宏觀流動涉及到的無量綱參數有：Da=K/L2，Pr=υ/αe，Ra=gβΔTL3/（υαe），J=υe/υ。其中，J為黏度比，溫差ΔT=Th-Tc。

2 耦合雙分布非正交MRT-LB模型

2.1 求解流場的M RT-LB模型

格子Boltzmann方程描述了氣體分子分布函數的時空演變，通過宏觀流動變量與分布函數之間的關系來得到宏觀流動參數。非正交轉換矩陣MRT-LB模型速度分布的演化方程為

式中：e為離散速度的方向矢量；δt為時間步長；Fx，t為合外力；fi（x，t）為密度分布函數；M 為速度空間的非正交轉換矩陣；S為對角松弛矩陣；I為單位矩陣；m和meq分別為f和feq對應的矩空間分布函數，其表達式為

對于二維多孔方腔內的流動，速度場采用D2Q9離散速度模型，各離散速度的方向矢量為

對于非正交轉換矩陣M，可以表示如下：

對于矩空間平衡態分布函數meq，定義如下：

式中：ux和uy分別為u的水平分速度和豎直分速度；ρ為流體的密度。

速度場的平衡分布函數feq

i定義如下：

力項F定義如下：

本文所采用的MTR-LBM 方法分布函數的演化與常規的LBM 方法一樣，首先是矩空間進行碰撞：

式中：m*為發生過碰撞后的矩空間分布函數。當碰撞完成后，通過公式f*=M-1·m*將矩空間轉化為速度空間，遷移仍然在速度空間進行：

流體的宏觀密度ρ和速度u定義如下：

流體的宏觀速度u可以通過引入一個臨時速度v獲得：

式中：

2.2 求解溫度場的M RT-LB模型

在對溫度場的求解時，采用計算量相對較小的D2Q5模型進行模擬，其模擬結果已經足夠準確，二維溫度場的非正交MRT-LBM方程定義如下：

式中：g為溫度分布函數；geq為溫度平衡分布函數；Q為松弛系數矩陣；N為非正交轉換矩陣，其表達式為

對于D2Q5模型，其溫度場各離散速度的方向矢量為

模型中的溫度T通過式（24）計算：

矩空間平衡態分布函數neq為

松弛系數矩陣Q定義如下：

速度空間的溫度平衡態分布函數geqi 定義為

式中：ω為常值。

松弛時間τυ和τT公式如下：

3 數值模擬結果及討論

3.1 模型驗證

圖3 不同方法流線及等溫線圖比較（Ra=106，Da=10-4，ε=0.4，Pr=1.0）Fig.3 Comparison of stream lines and isotherms with different methods（Ra=106，Da=10-4，ε=0.4，Pr=1.0）

表1 熱壁面上平均Nusselt數Table 1 Average Nusselt numbers at hot wall surface

3.2 冷源布置方案的影響

為了研究何種形式的冷卻邊界更有利于方腔內熱量的傳遞，圖4給出了ε=0.4、Ra=106、Da=10-2時所研究的6種不同冷源布置方案下方腔內溫度場及流場變化圖。圖中：Ψmax，l、Ψmax，r分別表示左右漩渦的最大值。可以看出，在給定參數下，6種不同冷源布置方案等溫線圖差異很大，方腔內等溫線只在冷壁面及熱壁面附近的薄邊界層內保持垂直，方腔中的等溫線基本保持水平，表明此時腔內對流傳熱占主導地位。通過計算方腔內高溫塊熱壁面上的平均Nusselt數及等溫線圖，部分冷源置于方腔上側腔內的對流傳熱強烈，Scheme A、Scheme B、Scheme D的布置方案更有利于熱量的傳遞。從流線圖可以看出，Scheme C由于冷源空間結構的限制，導致方腔2個漩渦均集中于方腔的下部，流體的流動在方腔上部很弱，導致方腔上部的熱量無法有效傳遞，方腔內整體傳熱效果很差。通過分析Scheme C、Scheme F兩個布置方案，Scheme F布置方案熱壁面Nuave及內熱源右側的漩渦強度明顯高于Scheme C，表明提高冷源的布置位置能夠明顯提高腔內的對流傳熱強度。由流線圖看出，方腔內產生2個流動方向相反的漩渦，分別位于內熱源的左右兩側。由于6種不同的冷源布置方案，漩渦形狀及渦心位置發生改變。冷源對稱布置時，方腔內形成的2個漩渦也對稱分布，方腔左右壁面冷源非對稱布置時，2個漩渦不再對稱分布，左右漩渦中心位置隨冷源位置移動。

為了更直觀地了解方腔內流體的流動特性，圖5給出了ε=0.4、Ra=3×106、Da=10-2時6種不同冷源布置方案下方腔中間高度上的速度分布V。以高溫塊中心點為坐標原點，向右為速度的正方向。由圖5知，在高溫熱塊左側速度首先達到負的最大值然后達到正的最大值，從而論證了方腔左側形成的是逆時針運動的漩渦。Scheme A的布置方案方腔中間高度速度變化范圍最大，而且其峰值最大，峰形更為尖銳，說明冷源雙上部的布置方式流體具有很好的流動狀態，對流傳熱作用強烈，Scheme C的布置方案其峰值最小，說明了雙下部的冷源布置方案限制了方腔內流體的流動，腔內的傳熱效果很差。圖6給出了ε=0.4、Ra=6×105時6種不同冷源布置方案下內熱源熱壁面上平均Nusselt數隨Darcy數的變化。可以看出，Da＜10-4時，熱壁面平均Nusselt數隨Darcy數的增大增加緩慢，表明方腔內的傳熱以熱傳導為主，此時Scheme B中部冷源的布置方案表現出最大的熱壁面平均Nusselt數。Darcy數從10-4增加到10-3時，熱壁面平均Nusselt數增幅劇烈，以Scheme B為例，平均Nusselt數從4.75增大到7.55，增大到原來的1.59倍。Darcy數增大到10-2后，隨著Darcy數的增大，熱壁面平均Nusselt數增加緩慢，在Da＞10-2時，Scheme A布置方案熱壁面具有最大的平均Nusselt數。這是由于Darcy數較大時，腔內流體流動受到的介質阻力較小，流動可以看做是純浮升力的誘導流，方腔內的對流傳熱占據主導地位。

圖4 六種不同等溫邊界布置下的等溫線圖和流線圖Fig.4 Isotherms and stream lines for six different isothermal boundary arrangements

圖5 六種不同等溫邊界布置下的速度分布Fig.5 Velocity profiles with six different isothermal boundary arrangements

圖6 六種不同等溫邊界布置下熱壁面平均Nusselt數隨Darcy數變化Fig.6 Variation of average Nusselt number at hot wall surface with Darcy number under six different isothermal boundary arrangements

3.3 Rayleigh數的影響

Rayleigh數作為影響方腔內自然對流傳熱的重要參數具有重要意義。當Ra＜104時，方腔內傳熱以熱傳導為主，Ra=106時，對流傳熱已經占據主導地位。為了研究Rayleigh數對6種不同冷源布置方案下方腔內對流傳熱的影響，圖7給出了ε=0.4、Da=10-2時，6種不同冷源布置方案下熱壁面平均Nusselt數隨Rayleigh數的變化情況。可以看出，Rayleigh數在1×106～6×106變化過程中，Scheme A的布置方案熱壁面平均Nusselt數均為最大值，Scheme C冷源布置方案熱壁面平均Nusselt數始終為最小值。Scheme A、Scheme B、Scheme D三種布置方案熱壁面均有較高的平均Nusselt數，Scheme A熱壁面平均Nusselt數從9.71增大到了14.49。圖8給出了6種冷源布置方案下內熱源熱壁面上的局部Nusselt數變化情況，起始于內熱源右壁面，隨之左壁面下壁面，然后終止于內熱源上部熱壁面。可以看出，方腔內部等溫塊邊角處的等溫線聚集，表明此處有較大的溫度梯度，局部Nusselt數在等溫熱塊的邊角出現尖峰。這可以從流體流動角度解釋，由于熱塊邊角不是流體流過的前緣就是后緣導致此處有較大的溫度梯度。Scheme A的布置方案在內熱源的左壁面、右壁面及上壁面均表現出最大的局部Nusselt數，解釋了其有最大的熱壁面平均Nusselt數。

圖7 六種不同等溫邊界布置下熱壁面平均Nusselt數隨Rayleigh數變化Fig.7 Variation of average Nusselt number at hot wall surface with Rayleigh number under six different isothermal boundary arrangements

3.4 內熱源結構形式及位置的影響

為了研究方腔內放置何種形式的高溫熱塊更有利于熱量的傳遞，本節比較了3種不同內熱源結構形式對腔內流體流動與傳熱的影響。3種模型內熱源的熱壁面周長相等，冷源布置方案為Scheme A采用的雙上部冷卻方式，圖9給出了ε=0.4、Ra=3×105時不同Darcy數下3種不同形狀的內熱源的等溫線圖和流線圖，流線圖給出了漩渦的最大流函數值，等溫線圖給出了熱壁面上平均Nusselt數值。由圖9可以看出，Da=10-3時，3種方案的等溫線分布均較為均勻，等溫線大都垂直分布，表明此時腔內的傳熱是導熱占主導地位，當Da=10-2時，腔中的等溫線開始發生明顯彎曲，只有冷邊界與熱塊薄邊界層內等溫線保持垂直狀態，這是由于Darcy數較大，腔內多孔介質的滲透率增大，流體受到方腔內部介質的阻力減小，腔內的對流傳熱作用增加，更有利于腔內的傳熱。由圖9看出，在高Darcy數下，腔內垂直布置的內熱源熱表面具有最大的平均Nusselt數值并且最大流函數值亦為最大，表明在此條件下Case 3內熱源的布置方式更有利于方腔內的散熱。

圖9 三種不同內熱源結構形式下方腔內的流線圖和等溫線圖Fig.9 Stream lines and isotherms in square carity under three different internal heat source structures

圖10 三種不同內熱源結構形式下熱壁面平均Nusselt數隨Rayleigh數變化Fig.10 Variation of average Nusselt number at hot wall surface with Rayleigh number in cavity under three different heat source structures

圖10給出了3種不同內熱源結構形式下熱壁面上平均Nusselt數隨Rayleigh數的變化情況。當Ra＜104時，熱壁面上平均Nusselt數增加緩慢，此時腔內以導熱為主，且水平放置的高溫熱塊（Case 2）方腔熱壁面具有更大的Nusselt數。當Ra＞105時，熱壁面上平均Nusselt數隨Rayleigh數增加劇烈，表明方腔內流體的對流傳熱作用在加強，在高Rayleigh數下，豎直放置的高溫熱塊（Case 3）腔內具有更好的對流傳熱特性。為了更好地說明在高Rayleigh數下3種不同內熱源結構形式下方腔內的對流傳熱效果，圖11給出了ε=0.4、Ra=106、Da=10-2時，3種不同內熱源結構形式下方腔中間高度上的速度分布v。以熱塊中心為坐標原點，向右為速度正方向，圖中y=0.5為方腔中間高度的水平線，由于邊界條件及內熱源的對稱布置方腔中間高度上的速度分布v亦表現出對稱特性。由圖11可以看出，Case 3內熱源結構形式具有速度的最大值，Case 2結構形式峰值最小，Case 3結構形式流動覆蓋范圍更廣，腔內具有更好的傳熱特性。

圖11 三種不同內熱源結構形式下的速度分布Fig.11 Velocity profiles in cavity under three different internal heat source structures

為了研究內熱源位置對多孔方腔內流體流動及傳熱的影響，選取了H/L=0.4的高溫熱塊為研究對象（即Case 1），高溫熱塊上壁距離方腔上壁面的無量綱長度a，其他無量綱參數設置如下：ε=0.4，Da=10-3，Ra=106。圖12給出了4種內熱源處于不同位置時（高溫熱塊中心始終處于方腔中垂線上），方腔內的等溫線圖和流線圖。圖中：左圖為等溫線圖，右圖為流線圖。由圖12可以看出，a=0時漩渦集中于方腔上部，漩渦基本對下部不產生影響，腔內下部的傳熱效果極差，隨著a值的增大，漩渦開始覆蓋方腔下部，下部的傳熱得到明顯提高。為了進一步理解內熱源位置對腔內對流傳熱的影響，圖13給出了a=0，0.1，0.2，…，0.6時，7種不同內熱源位置熱壁面局部Nusselt數的變化。熱壁面起始于左壁面，向右依次為右壁面、上壁面，終止于下壁面，橫坐標表示熱壁面的無量綱長度。由于內熱源及冷源的對稱布置，由圖13看出內熱源的左右壁面局部Nusselt數均相同，a=0時左、右熱壁面的局部Nusselt數始終最小，下部熱壁面則具有較大的局部Nusselt數。a=0.6時，內置熱塊左壁面、右壁面、上壁面均表現出最大的局部Nusselt數值。圖14則給出了內熱源熱壁面上平均Nusselt數隨無量綱長度a的變化。圖中：散點表示模擬結果，曲線是運用最小二乘法得到的關于Nusselt數的擬合關系式的曲線，擬合度為91.1%，擬合的曲線關系式為：Nuave=7.73e0.99a-1.67a2。可以看出，隨著無量綱長度a的增大，熱壁面上的平均Nusselt數表現出先增大后減小的趨勢，并且在a=0.25處，熱壁面上的平均Nusselt數達到最大值，此時方腔內的對流傳熱水平最為強烈。

圖12 不同a值下方腔內的等溫線圖和流線圖Fig.12 Isotherms and stream lines in square cavity at different values of a

圖13 不同a值下熱壁面局部Nusselt數變化Fig.13 Variation of local Nusselt number at hot wall surface under different a values

為了研究方腔位于不同水平位置時（高溫塊中心點始終位于方腔水平線上），對腔內對流傳熱的影響，選取了H/L=0.4的高溫熱塊為研究對象，高溫熱塊左壁面距離方腔左壁面的無量綱長度b=0，0.025，0.05，…，0.6，其他無量綱參數設置如下：ε=0.4，Da=10-3，Ra=106。圖15給出了不同b值下方腔內的等溫線圖和流線圖，同時給出了熱壁面上的平均Nusselt數及最大流函數值。圖中：左圖為等溫線圖，右圖為流線圖。可以看出，由于此時高溫塊的布置不再左右對稱，腔內流線圖亦表現出非對稱性，b=0時腔內形成2個大小不一的漩渦，在高溫塊的上部形成一個逆時針運動的小漩渦，右側形成一個順時針運動的大漩渦。隨著b值增加，由于方腔空間結構的改變上部漩渦開始向高溫熱塊左側移動，為了更好地理解b值對腔內傳熱的影響，圖16給出了不同b值下方腔中間高度的垂直速度分布v。以熱塊中心為坐標原點，向右為速度正方向，圖中y=0.5表示方腔中間高度的水平線，可以看出，在b由0增大至0.3的過程中，腔內左側速度峰值亦依次增大，在高溫塊的右側不同b值下的速度分布大致相同且峰值相差不大。

圖15 不同b值下方腔內的等溫線圖和流線圖Fig.15 Isotherms and stream lines in square cavity at different values of b

圖16 不同b值下方腔水平中平面的速度分布Fig.16 Velocity profiles in the horizontalm id-p lane of square cavity under different b values

圖17 熱壁面平均Nusselt數隨b值變化Fig.17 Average Nusselt number at hot wall surface versus b values

為了進一步研究內熱源位于不同水平位置對腔內自然對流傳熱的影響，圖17給出了方腔熱壁面上平均Nusselt數隨無量綱距離b的變化，同時圖18亦給出不同b值下方腔熱壁面上的局部Nusselt數變化。可以看出，熱壁面平均Nusselt數隨b的變化表現出對稱分布，對稱軸為b=0.3。因此只討論0≤b≤0.3，在b從0增大的過程中，熱壁面平均Nusselt數首先突增，在b=0.025處達到最大值，為16.0，隨后平均Nusselt數開始減小，在b=0.15處達到最小值8.68，隨后在b增大到0.3的過程中，熱壁面上的平均Nusselt數又表現出上升的趨勢。此現象可以從圖18熱壁面上的局部Nusselt數變化得到解釋，在b=0.025時由于冷熱壁面距離較近，高溫熱塊左壁面上部附近存在極大的溫度梯度，導致熱塊左壁面上側具有很大的局部Nusselt值，且高溫熱塊上部也有很大的局部Nusselt數，使得平均Nusselt數很大。b=0.15時，高溫熱塊熱壁面的右壁面和下壁面的局部Nusselt均表現出較小值，且另外2個熱壁面的局部Nusselt數處于中間水平，導致平均Nusselt值很小。

圖18 不同b值下熱壁面局部Nusselt數變化Fig.18 Variation of local Nusselt number at hot wall surface under different b values

4 結 論

本文采用非正交MRT-LBM對含內熱源的二維多孔方腔內的自然對流傳熱進行了數值模擬，研究了不同冷源布置方案、Darcy數、Rayleigh數、內熱源結構形式及其位置對方腔內對流傳熱的影響。研究結果表明：

1）當方腔內冷卻邊界及內熱源左右對陣布置時，腔內溫度場、流場亦左右對稱分布。

2）對于討論的6種冷源布置方案，雙頂部布置方式在高Rayleigh數下有更好的冷卻效果。

3）在高Rayleigh數下，垂直放置的內熱源腔內具有更好的傳熱效果。

4）當高溫方形熱源尺寸H/L=0.4為定值時，高溫熱源位置a對腔內的對流傳熱影響顯著，隨著a值的增大，熱壁面平均Nusselt數表現出先增大后減小的趨勢，且存在最佳的位置（a=0.25），使方腔內的對流傳熱效果最強。

5）當高溫方形熱源尺寸H/L=0.4為定值時，高溫熱源位置b亦對腔內的對流傳熱有著重要影響，在b從0～0.3變化過程中，熱壁面平均Nusselt數呈現出先增后減再增的趨勢。