張凱寧 中國航空工業(yè)集團公司西安航空計算技術(shù)研究所

引言

多孔介質(zhì)內(nèi)的流動與換熱是自然界中常見的一種現(xiàn)象，如地?zé)豳Y源的開發(fā)與利用、谷物 的儲存、食品工程中的熱對流以及生物技術(shù)等。隨著計算機技術(shù)的不斷發(fā)展，數(shù)值模擬已成經(jīng)為該領(lǐng)域的主要研究方法之一。國內(nèi)外許多學(xué)者運用傳統(tǒng)數(shù)值方法如有限差分法、有限容積法、有限元法、譜方法等對其進行了大量的研究。作為近二十年來興起的介觀數(shù)值模擬方法，格子Boltzmann 方法（LBM）可以非常方便地描述多孔介質(zhì)流體內(nèi)部的相互作用，因此格子Boltzmann 方法已成為研究多孔介質(zhì)內(nèi)流動與換熱現(xiàn)象的很有力的數(shù)值工具。多孔介質(zhì)內(nèi)的混合對流在工程實際中大量存在，對其開展研究具有重要意義。對于多孔介質(zhì)內(nèi)的混合對流，由于要同時考慮自然對流、強迫對流以及多孔介質(zhì)固體骨架對流動和換 熱的影響 ，因此對其實現(xiàn)準(zhǔn)確模擬較為困難。Khanafer 和 Chamkha 采用修正的 Brinkman-Darcy 模型對多孔介質(zhì)腔體內(nèi)的頂蓋驅(qū)動混合對流進行了數(shù)值模擬，研究了 Ri 數(shù) 以及 Da 數(shù)對流動換熱的影響。Kumar 等采用修正的 Brinkman-Forchheimer-Darcy 模型對多孔介質(zhì)方腔內(nèi)的混合對流進行了數(shù)值研究，研究了不同無量綱參數(shù)以及非線性阻力對流動與換熱的影響。本文將采用 LBM 對多孔介質(zhì)方腔內(nèi)的混合對流進行數(shù)值模擬，研究不同特征參數(shù)以及非線性阻力對流動與換 熱的影響。

1 不可壓單松弛格子 Boltzmann 模型

對于不可壓縮流體在多孔介質(zhì)內(nèi)的自然對流傳熱，假設(shè)流動處于層流狀態(tài)，忽略流動中的粘性熱耗散，且流體和多孔介質(zhì)固體骨架之間滿足局部熱平衡條件。在 Boussinesq 假設(shè)成 立的條件下，基于 Brinkman-Forchheimer-Darcy 模型， 結(jié)合能量方程，則多孔介質(zhì)內(nèi)熱流動 的宏觀控制方程可以寫成如下形式：

在 Boussinesq 假設(shè)成立的條件下，基于文獻[6]修正的Brinkman-Forchheimer-Darcy 模型， 多孔介質(zhì)內(nèi)混合對流換熱的無量綱宏觀控制方程可以寫成：

其中，u ， p 和 T 分別為流體的無量綱體積平均速度、壓力和溫度；ε為多孔介質(zhì)的孔隙率，Re 為有效雷諾數(shù)，Pr 為普朗特數(shù)。F為流體的總力，其表達式為：其中Gr 為格拉曉夫數(shù)，Da 為達西數(shù)，k 為 y-方向的單位向量，，u 和v 分別 為速度u 的 x 和 y 方向分量。Fε多孔介質(zhì)的幾何形狀因子， 由 Ergun 經(jīng)驗公式給出[7]：

用于不可壓縮流涕流動與換熱的熱格子Boltzmann 模型通?？梢苑譃槿悾憾嗨俣饶Ｐ汀㈦p分布函數(shù)模型和混合模型。本文將采用雙分布函數(shù)格子 Boltzmann 模型研究多孔介質(zhì)內(nèi)的混合對流換熱問題。該模型用一個密度分布函數(shù)來模擬速度場，而用一個溫度分布函數(shù) 來求解溫度場。采用 D2Q9 模型，模擬速度場和溫度場的演化方程分別如下[8]：

式中，r 為空間位置矢量，δt為時間步長，fi為密度分布函數(shù)，Ti為溫度分布函數(shù)，fieq為 對平衡態(tài)密度分布函數(shù)，eqig 為平衡態(tài)溫度分布函數(shù)，τ和τT分別為速度場和溫度場的無量 綱松弛時間，Si為外力項。平衡態(tài)密度分布函數(shù)fieq和平衡態(tài)溫度分布函數(shù)gieq分別選取如 下：

其中ωi為權(quán)系數(shù)，cs為格子聲速。D2Q9 模型的離散速度為e0=(0，0)，e1=-e3=(1，0)c，e2=-e4=(0，1)c，e5=-e74=(1，1)c，e6=-e8=(0，1)c，其中c=δx/δt，δx為網(wǎng)格步長。D2Q9 模型的權(quán)系數(shù)為ω0=4/9，ω1-4=1/9，ω5-8=1/36 ，格子聲速為。

這樣，式（6）-（9）構(gòu)成了用于求解多孔介質(zhì)內(nèi)混合對流換熱的格子Boltzmann 模型。利用Chapman-Enskog 多尺度展開方法，在不可壓限制條件下，可以準(zhǔn)確地由該模型還原出宏觀控制方程（1）-（3）。邊界處理在格子Boltzmann 方法中起著重要作用，它會對數(shù)值計算的精度、穩(wěn)定性以及計算效率產(chǎn)生很大的影響。對于速度邊界條件和溫度邊界條件，本文采用非平衡態(tài)外推方法處理，詳見文獻[8]

2 .多孔介質(zhì)方腔內(nèi)的混合對流模擬

2.1 物理模型和邊界條件

多孔介質(zhì)方腔內(nèi)的混合對流換熱如圖1 所示。方腔高為H，內(nèi)部填充飽和多孔介質(zhì)固 體骨架材料。方腔的上下壁面絕熱，右壁面溫度為Th，左壁面溫度為Tc(Th＞Tc)。方腔的 左壁面以恒定的速度沿 y 軸正方向運動，右壁面以相同大小的速度沿y 軸的負方向運動。除左右邊界外，方腔內(nèi)其余位置的初始速度為 0.

圖1 多孔介質(zhì)方腔內(nèi)混合對流換熱物理模型

該問題的速度邊界條件和溫度邊界條件可寫為：

左邊界：u=0，v=1.0，T=0 ；右邊界：u=0 ，v=-1.0，T=1.0 ；

上邊界：u=v=0， ?T/?y=0；下邊界：u=v=0， ?T/?y=0。

2.2 模擬結(jié)果與分析

為了驗證本文數(shù)值方法的可靠性，先用 LBM 模擬了二維方腔內(nèi)的頂蓋驅(qū)動混合對流換 熱問題。相關(guān)模擬參數(shù)為Pr=0.71 ，Re取100 、400 和1000 ，對應(yīng)的Ri 數(shù)分別為1×10-2，6.25×10-4和1×10-4，網(wǎng)格密度取為 128×128 。為了進行定量比較，計算了沿方腔上壁面的平均Nu 數(shù)，并與已有文獻結(jié)果進行了對比。從表1 可以看出，本文 LBM 計算結(jié)果與文獻中其它數(shù)值方法的結(jié)果吻合得很好，充分證明了 LBM 方法的可靠。

表1 LBM 模擬結(jié)果與文獻結(jié)果的比較ε=0.999，Da=106

為了研究Da 數(shù)和Gr 數(shù)對混合對流換熱的影響，本文用LBM模 擬 不 同Da 數(shù)：Da=0.001，0.01，0.1 在Gr=102和Gr=104兩 種情形下的混合對流換熱過程。網(wǎng)格密度取為256×256，Ri=0.01，=0.9。圖2 給出了Gr=102時不同Da 數(shù)下的流線和等溫線。當(dāng)Da 數(shù)較低時（Da=0.001），方腔的冷壁面和熱壁面各出現(xiàn)一個小渦，傳熱主要是由熱壁面和冷壁面的熱傳 導(dǎo)引起的，等溫線幾乎是垂直的。隨著Da 數(shù)的增加，冷壁面的渦逐漸上移并對上壁面進行冷卻，而熱壁面的渦逐漸下移并對下壁面進行加熱，方腔內(nèi)的對流換熱作用逐漸加強。當(dāng) Da=0.1 時，方腔冷壁面和熱壁面的兩個小渦合為一個大渦，傳熱機制逐漸由熱傳導(dǎo)占主導(dǎo)地位變?yōu)閷α髡贾鲗?dǎo)地位，等溫線在方腔中央逐漸變得水平，但此時Gr 數(shù)較低，方腔內(nèi)的對流換熱作用主要是由左右壁面運動產(chǎn)生的強制對流引起的。

圖3 給出了4Gr=104時不同Da數(shù)下的流線和等溫線。當(dāng)Da=0.001時，在方腔的冷壁面和熱壁面各出現(xiàn)一個小渦，等溫線在上壁面的右邊和下壁面的左邊聚集。隨著Da 數(shù)的增加（Da=0.01，0.1），方腔中央形成一個橢圓形的大渦，整個流場出現(xiàn)比較明顯的邊界層流動的特征，從而使方腔內(nèi)的換熱作用得到增強，此時對流換熱占主導(dǎo)地位。d

圖3 不同Da 數(shù)下的流線（左）和等溫線（右）（Gr=104，Ri=0.01，=0.9）

從圖4a 可以看出，沿冷壁面的Nu 數(shù)從上到下逐漸增加，由此可知隨著壁面高度的增加，對流換熱的作用不斷減少。當(dāng)Gr=104時，從 圖4b 可以看出，沿冷壁面的局部 Nu 數(shù)變化趨勢和圖4a 類似，但是在 Gr=104的情形下，由自然對流產(chǎn)生的浮升力對方腔內(nèi)的換熱起到的作用開始變得明顯。圖4 還給出了在不考慮非線性介質(zhì)阻力時（Fε=0），Gr=102和 Gr=104時不同 Da 數(shù)下沿冷壁面的局部 Nu 數(shù)。從圖4a 和圖4b 可以看出，在相同的 Da 數(shù)和 Gr 數(shù)下， 不考慮非線性介質(zhì)阻力影響時的結(jié)果比考慮非線性介質(zhì)阻力影響時的結(jié)果偏大，而在 Da 數(shù) 和 Gr 數(shù)較高的情況下，偏差會更為顯著。表2 給出了 Gr=102和 Gr=104時不同 Da 數(shù)下冷壁面的平均 Nu 數(shù)（ Nu ）計算結(jié)果。從表2 可以看出，本文 LBM 的模擬結(jié)果與文獻[5]的模擬結(jié)果吻合得很好。

圖4 不同Da 數(shù)下沿左壁面的局部Nu 數(shù)（Gr=104，Ri=0.01，=0.9）

表2 冷壁面平均 Nu 數(shù)本文模擬結(jié)果與文獻[5]模擬結(jié)果的比較ε=0.9，Ri=0.01

3 .結(jié)論

本文采用了D2Q9 離散速度模型，并以格子Boltzmann 模型進行了多孔介質(zhì)方腔內(nèi)的混合對流換熱的數(shù)值仿真計算，研究了Da 數(shù)和Gr 數(shù)對混合對流換熱的影響，Da 數(shù)值越大，多孔介質(zhì)阻力對流動的影響越小，Ri 對發(fā)熱邊界影響越顯著，數(shù)值仿真的計算結(jié)果與傳統(tǒng)計算結(jié)果對比，結(jié)果一致性較好，研究結(jié)果在對流換熱的數(shù)值仿真計算方面具有一定的指導(dǎo)意義。