多孔介質(zhì)方腔內(nèi)混合對(duì)流換熱的仿真模擬

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

引言

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

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

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

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

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

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

式中，r 為空間位置矢量，δt為時(shí)間步長，fi為密度分布函數(shù)，Ti為溫度分布函數(shù)，fieq為 對(duì)平衡態(tài)密度分布函數(shù)，eqig 為平衡態(tài)溫度分布函數(shù)，τ和τT分別為速度場(chǎng)和溫度場(chǎng)的無量 綱松弛時(shí)間，Si為外力項(xiàng)。平衡態(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)混合對(duì)流換熱的格子Boltzmann 模型。利用Chapman-Enskog 多尺度展開方法，在不可壓限制條件下，可以準(zhǔn)確地由該模型還原出宏觀控制方程（1）-（3）。邊界處理在格子Boltzmann 方法中起著重要作用，它會(huì)對(duì)數(shù)值計(jì)算的精度、穩(wěn)定性以及計(jì)算效率產(chǎn)生很大的影響。對(duì)于速度邊界條件和溫度邊界條件，本文采用非平衡態(tài)外推方法處理，詳見文獻(xiàn)[8]

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

2.1 物理模型和邊界條件

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

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

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

左邊界：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é)果與分析

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

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

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

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

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

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

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

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

3 .結(jié)論

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