多孔介質方腔內混合對流換熱的仿真模擬

張凱寧 中國航空工業集團公司西安航空計算技術研究所

引言

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

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

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

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

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

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

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

其中ωi為權系數，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為網格步長。D2Q9 模型的權系數為ω0=4/9，ω1-4=1/9，ω5-8=1/36 ，格子聲速為。

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

2 .多孔介質方腔內的混合對流模擬

2.1 物理模型和邊界條件

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

圖1 多孔介質方腔內混合對流換熱物理模型

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

左邊界：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 模擬結果與分析

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

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

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

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

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

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

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

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

3 .結論

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