通風及內熱源參數的方腔內混合對流模擬

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

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

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

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

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

1 方腔內混合對流模型的建立

1.1 問題描述

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

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

1.2 計算模型

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

(1)

(2)

(3)

(4)

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

(5)

(6)

式中：n為沿熱表面外法線方向的矢量；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處的速度分布函數矢量；m(r,t)、m(eq)(r,t) 為矩和矩對應的平衡態；M為正交轉換矩陣，并滿足m=Mf；S為馳豫系數矩陣；δt、δx為單位長度的離散時間步長和格子間距；c=δx/δt為格子速率。

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

(8)

(9)

(10)

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

(11)

(12)

圖2 D2Q9模型Fig.2 Model of D2Q9

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

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

(13)

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

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

(14)

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

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

(15)

(16)

(17)

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

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

1.3 模型驗證

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

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

表1 模擬結果對比Table 1 Comparison of simulation results

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

2 對流換熱的計算結果與分析

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

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

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

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

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

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

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

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

3 結論

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

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

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