孫夢楠， 宋桂秋， 周世華， 董祉序

(1.東北大學(xué) 機械工程與自動化學(xué)院，遼寧 沈陽 110819；2.沈陽工業(yè)大學(xué) 機械工程學(xué)院，遼寧 沈陽 110870)

在針對換熱器[1]、空調(diào)系統(tǒng)[2]、電子設(shè)備冷卻[3]等眾多工程設(shè)計和能源問題的研究中，被分析流體往往受到外力和溫度的共同影響，多以混合對流形式存在。合理設(shè)計和布置設(shè)備結(jié)構(gòu)不僅可以有效加大能源利用率，還能夠減少重要零件的損耗。利用方腔模型可以分析多種復(fù)雜場內(nèi)流體的運動特性和傳熱規(guī)律[4]。在針對混合對流的研究中，也多采用具有內(nèi)部加熱元件的開口方腔模型進行模擬。Radhakrishnan等[5]通過實驗證明了模型的可靠性，并指出采用數(shù)值模擬方法有利于優(yōu)化熱流表現(xiàn)。孫猛等[6]在不同的格拉曉夫數(shù)、瑞利數(shù)以及縱橫比的條件下，研究了外界來流對方腔內(nèi)混合對流換熱的影響。Karimi等[7-8]建立了具有多個加熱單元的復(fù)雜混合對流模型，并對熱源的大小和位置進行了討論。Chamkha等[9-10]借助于方腔模型對多種結(jié)構(gòu)因素進行分析，研究結(jié)果表明，進、出流口的位置對混合對流的影響很大，在改善換熱效率方面起著關(guān)鍵作用。

隨著流體流動和換熱機理研究的深入，新的數(shù)值計算方法也在不斷涌現(xiàn)。格子Boltzmann方法(LBM)作為介觀模擬算法的代表，在近幾十年里發(fā)展迅速。與傳統(tǒng)計算方法相比，LBM使用簡化的動力學(xué)模型來處理微觀過程，這使其在具有出色并行計算能力的同時能夠更直觀地表現(xiàn)粒子間的相互作用[11]。多參數(shù)弛豫時間(MRT)的引入也進一步提高了LBM的計算精度和穩(wěn)定性[12]。目前，LBM主要通過多速度法[13]、多分布法[14]和混合法[15]對多場問題進行模擬，并在混合對流、輻射傳熱、雙擴散對流等問題上得到了驗證。

綜上可以發(fā)現(xiàn)，在分析腔內(nèi)熱源對混合對流的影響時，研究工作主要集中在熱元件的布置方式和尺寸上，而對熱元件的形狀特征很少提及。此外，內(nèi)熱源結(jié)構(gòu)和外部方腔的開口位置共同變化下的流體和換熱表現(xiàn)也很少被討論。因此，本文采用結(jié)合MRT-LBM與有限差分法(FDM)的耦合算法[15]，通過模擬流函數(shù)線圖、等溫線圖和努賽爾數(shù)，探討了不同進、出流口位置下熱源形狀變化對方腔內(nèi)混合對流的流動、溫度及換熱特性的影響。

1 方腔內(nèi)混合對流模型的建立

1.1 問題描述

本文研究的含內(nèi)熱源的開口方腔結(jié)構(gòu)如圖1所示。計算基于笛卡爾坐標(biāo)系(x,y)，坐標(biāo)原點位于計算域的左下角。正方形腔體與等長寬中心加熱元件的尺寸分別定義為L和d=0.3L。方腔內(nèi)壁為具有溫度Tc的恒溫壁面，熱源外表面的溫度為Th(Th>Tc)。

具有速度uin和溫度Tin=Tc的流體通過左腔壁上尺寸為w=0.2L的進流口流入，最終從對側(cè)壁面上相同尺寸的出流口以速度uout和溫度Tout流出，流動過程中所受重力加速度為g。以進、出流口下邊緣高度l和l′ 來定義開口位置，并分別在底部、中部和頂部開口下進行討論。同時，引入圓角半徑R來討論加熱元件的形狀變化，R的值取為0、0.25d和0.5d。由圖1中可見，R=0和R=0.5d分別表示熱源具有正方形截面和圓形截面。

1.2 計算模型

本文采用普朗特數(shù)Pr為0.71的空氣作為工作流體，認(rèn)為它是滿足Boussinesq近似的牛頓流體，并以穩(wěn)態(tài)、層流運動。因此，可以得到混合對流的無量綱控制方程為：

(1)

(2)

(3)

(4)

引入努賽爾數(shù)來表征熱源表面的對流換熱強度。由局部溫度梯度可求得的局部努賽爾數(shù)Nu和平均努賽爾數(shù)Nuav分別為：

(5)

(6)

式中：n為沿?zé)岜砻嫱夥ň€方向的矢量；A為熱源廓線周長。

用于求解速度場的MRT-LBM采用如圖2所示的9速度模型(D2Q9)，其演化方程為：

f(r+ciδt,t+δt)-f(r,t)=

-M-1S(m(r,t)-m(eq)(r,t))+F

(7)

式中：f(r,t)=(f0,f1, …,f8)T為位置r處的速度分布函數(shù)矢量；m(r,t)、m(eq)(r,t) 為矩和矩對應(yīng)的平衡態(tài)；M為正交轉(zhuǎn)換矩陣，并滿足m=Mf；S為馳豫系數(shù)矩陣；δt、δx為單位長度的離散時間步長和格子間距；c=δx/δt為格子速率。

各個速度方向上的粒子速度集c，權(quán)重系數(shù)集ω和體積力F分別定義為：

(8)

(9)

(10)

在LBM中，宏觀物理量被視為分布函數(shù)的平均特性，因此流體密度ρ和速度u可以表示為：

(11)

(12)

圖2 D2Q9模型Fig.2 Model of D2Q9

在MRT-LBM的算式中，方腔內(nèi)壁面和進、出流口位置采用非平衡態(tài)外推邊界條件[16]。對于曲線邊界與網(wǎng)格不完全一致的熱源表面，有些粒子在被壁面反彈后會脫離節(jié)點，落到間隙位置(圖3中的xg位置)，則采用線性插值邊界條件[17]。如圖3所示，該方法首先假設(shè)分布函數(shù)f1(xf,t)可以穿過邊界節(jié)點xw到達實體節(jié)點xs，由線性插值得到：

f1(xw,t+δt)=(1-q)f1(xf,t+δt)+qf1(xs,t+δt)

(13)

式中q為流體節(jié)點到邊界間的距離與網(wǎng)格長度之比，q=|xf-xw|/|xf-xs|。再根據(jù)動量守恒和固體邊界的無滑移條件，最終得到在實體節(jié)點xs處的分布函數(shù)為：

(1-q)f1(xff,t)]

(14)

式(14)具有二階精度，并可以避免按q≤1/2和q>1/2劃分時所引起的數(shù)值振蕩。

在計算溫度場時，采用有限差分方法在與MRT-LBM相同的網(wǎng)格上對方程(4)進行求解[18]。這里采用一階向前和二階中心差分方法來分別求解時間和空間的偏微分方程，即有：

(15)

(16)

(17)

圖3 曲面插值邊界條件Fig.3 Interpolation boundary condition for curved wall

不同位置的無量綱熱邊界條件定義為：在方腔內(nèi)壁面和進流口位置，θ=0；在熱源表面，θ=1；在出流口位置，?θ/?X=0。

1.3 模型驗證

為確認(rèn)所采用算法的可靠性和準(zhǔn)確性，以文獻[9]中所描述的具有中部進流口、頂部出流口以及方形截面內(nèi)熱源的方腔混合對流模型作為算例進行驗證。通過對比圖4中給出的流場和溫度場分布(Re=200，Ri=1)可以看出，采用本文方法所得到的模擬結(jié)果與文獻中的結(jié)果具有很好的一致性。同時，表1中的對比表明，當(dāng)Ri變化時，2種計算的Nuav具有相同的趨勢，最大百分比差異為3.52%，說明本文的計算方法是具有足夠精度的。

圖4 本文和文獻[9]的流場和溫度場的算例驗證Fig.4 Comparison of streamlines and isotherms between the present work and that of ref.[9]

表1 模擬結(jié)果對比Table 1 Comparison of simulation results

考慮到算法對網(wǎng)格的敏感性，分別采用100×100、150×150、200×200和250×250的網(wǎng)格分布對本文的混合對流模型進行計算。當(dāng)網(wǎng)格超過200×200時，繼續(xù)加密網(wǎng)格對模擬結(jié)果的改進很小，因此在下述分析中均采用200×200的網(wǎng)格。

2 對流換熱的計算結(jié)果與分析

為保證混合對流換熱中始終存在強制對流和自然對流，且方腔內(nèi)的流體為層流狀態(tài)，本文在Re=200，Ri=1的工況下對參數(shù)進行討論。不同形狀的加熱元件在不同進、出流口位置下的流場、溫度場分布如圖5～7所示。

圖5給出了當(dāng)進流口位于底部(l=0)時腔內(nèi)混合對流的流函數(shù)線圖和等溫線圖。在圖5(a)中，當(dāng)出流口也位于底部(l′=0)時，流體從方腔底部直接流過，慣性力使加熱元件被逆時針環(huán)流包圍。由于強迫對流和自然對流的共同作用，在加熱元件左側(cè)和方腔右上角還會出現(xiàn)2個渦流。R的增大會使環(huán)流增強，渦流被擠壓。當(dāng)l′=0.4L時(圖5(b))，冷來流與加熱元件的右下角發(fā)生直接接觸，腔內(nèi)右部的流動發(fā)生變化。由于進流口附近的流體受到出流口位置的影響較小，所以熱源左側(cè)與腔壁間仍存在渦流。當(dāng)R=0時，有小部分來流會從左側(cè)繞經(jīng)熱源后再流至腔外，但隨著R的增大，腔內(nèi)逐漸出現(xiàn)與圖5(a)中相似的環(huán)流結(jié)構(gòu)。當(dāng)l′=0.8L時(圖5(c))，熱源被整個包含在來流中，R的增大會減少上部繞流，但也為渦流的流動提供了更大的空間。在出流口位置由l′=0升高到l′=0.8L的過程中，冷來流的上移使腔內(nèi)右部的等溫線向上收縮，熱源底部的等溫線也變得越來越密集，這說明該位置的具有較大的溫差。R的增大會使熱源附近的等溫線發(fā)生變化，但對整體溫度場的影響不大。

圖5 進流口位于底部時，R和出流口位置對流函數(shù)線和等溫線的影響Fig.5 Effects of R and outlet location on the streamlines and isotherms for cavity with bottom inlet

圖6給出的是當(dāng)進流口位于中部(l=0.4L)時腔內(nèi)混合對流的流函數(shù)線圖和等溫線圖。此時，加熱元件正對進流口，在進行熱交換的同時也阻礙了流體運動。流入的流體會在進流口兩側(cè)產(chǎn)生上大下小的兩個低溫角渦，并從熱源上下兩側(cè)流過，最終在出流口處匯合。等溫線圖表明，腔體左側(cè)溫度主要由來流決定，而當(dāng)流體流至右側(cè)時熱源的作用開始變得明顯起來。隨著出流口位置l′的升高，右側(cè)腔體內(nèi)的流函數(shù)線和等溫線受到較大影響，但腔內(nèi)流體的整體流動模式變化很小。R的增大會使渦流的運動區(qū)域變大，并使對應(yīng)位置的等溫線向右傾斜。

進流口位于頂部(l=0.8L)時腔內(nèi)混合對流的流函數(shù)線和等溫線如圖7所示。在圖7(a)和(b)中(l′=0和l′=0.4L)，腔內(nèi)的繞流雖然占主導(dǎo)地位，但它們難以覆蓋進流口下方的空間，因而產(chǎn)生了順時針渦流。等溫線聚集在熱源上表面，可以看出上部繞流對溫度場的影響要大于下部繞流。隨著R的增大，流體更多地從腔體上部流過，下部繞流減少，并逐漸脫離熱源表面。當(dāng)進流口和出流口都位于圖7(c)頂部(l=l′=0.8L)時，流體不再受加熱元件阻礙，幾乎水平地流過腔體頂部。在冷來流未經(jīng)過的區(qū)域內(nèi)，溫度也分布得較為均勻。R的增大會改變腔內(nèi)其他位置的流體運動，使部分渦流開始向環(huán)流轉(zhuǎn)變。由于重力的存在，腔體內(nèi)的流函數(shù)線和等溫線分布與圖5(a)并不對稱。

從上面的討論可以看出，進、出流口位置選擇對腔內(nèi)流場的影響要大于加熱元件的形狀變化，R的變化對溫度場的影響并不明顯。實際上，熱源形狀的特征主要反映在加熱表面的對流換熱上。

圖8給出了9種不同方腔結(jié)構(gòu)下的熱源表面局部努賽爾數(shù)Nu，圖中以熱源上表面與其中心線的交點為橫坐標(biāo)零點，沿逆時針方向繪制了Nu曲線，R為0、0.25d和0.5d時所對應(yīng)的熱表面的無量綱長度分別為1.2、1.07和0.94。當(dāng)R=0時，在熱源的4個直角位置會出現(xiàn)尖銳的峰值，這表明該位置的溫度梯度很大，在相應(yīng)的溫度場中也會出現(xiàn)了密集的等溫線。隨著R的增大，4個峰值開始下降并變得平緩，至熱源截面變?yōu)閳A(R=0.5d)時，Nu曲線中只剩下一個明顯的主峰值。更改進、出流口的位置可以改變峰值的大小和位置分布，相比而言，進流口的影響要大于出流口。

圖9中給出了不同方腔結(jié)構(gòu)下熱源表面的平均努賽爾數(shù)Nuav，圖中的大寫字母B、M和T分別表示底部、中部和頂部進流口；B′、M′和T′分別表示底部、中部和頂部出流口。與局部努賽爾數(shù)的最大值出現(xiàn)在R=0時相反，平均努賽爾數(shù)是隨著R的增大而增大的，這說明雖然一些特殊位置的局部對流換熱能力有所下降，但整體表面的換熱卻是增強的，且Nuav的值并不與換熱表面的面積大小成正比。另外可以看出，當(dāng)方腔的進、出流口在同一位置時，冷流體直接從頂部或底部通道流出腔體，Nuav較??；而當(dāng)方腔具有中部進流口和頂部出流口時，熱源與來流正面接觸，換熱充分，Nuav取得最大值。

3 結(jié)論

1)流場和溫度場分布與熱源形狀關(guān)系不大，而主要受進、出流口位置影響。當(dāng)進流口位于腔體底部或頂部時，隨著出流口位置與進流口位置之間距離的增大，熱源周圍的環(huán)流被擠壓，并逐漸被外部來流所取代；而當(dāng)進流口位于中部時，外部來流始終流經(jīng)熱源的上、下表面，出流口位置不改變腔內(nèi)流體的流動模式。在溫度場中，由于受到外部強制對流的影響，等溫線會向進流口附近的熱源表面以及出流口位置聚集。

2)當(dāng)熱源截面為圓形時(R=0.5d)，熱源表面的局部努賽爾數(shù)Nu曲線具有平緩變化的單峰值；隨著R的減小，Nu曲線上會逐漸出現(xiàn)多個尖銳峰值。Nu的最大值出現(xiàn)在熱源表面靠近進流口位置的一側(cè)，而出流口位置對Nu的影響不大。當(dāng)進、出流口分別位于頂部和底部時，可在正方形熱源截面(R=0)的左上角獲得最佳局部對流換熱位置。

3)當(dāng)進、出流口位置固定時，熱源表面的平均努賽爾數(shù)Nuav隨著R的增大而增大。當(dāng)R不改變時，Nuav的最大值在方腔具有中部進流口和頂部出流口時出現(xiàn)。 因此，在進、出流口分別位于中部和頂部時，采用圓形截面的熱源(R=0.5d)可以使熱源整體表面的對流換熱最強。