李 斌，李本文，2，陳元元，許學(xué)成

(1.武漢科技大學(xué)材料與冶金學(xué)院，湖北 武漢，430081；2.大連理工大學(xué)能源與動力學(xué)院，遼寧 大連，116024)

多孔介質(zhì)內(nèi)流體的自然對流現(xiàn)象廣泛存在于自然界及工程技術(shù)領(lǐng)域，近幾十年來，國內(nèi)外眾多研究者通過實驗、理論分析和數(shù)值模擬等圍繞該現(xiàn)象展開了大量研究。多孔介質(zhì)結(jié)構(gòu)具有很強的隨機性、復(fù)雜性，得益于計算機技術(shù)的不斷發(fā)展，研究者常采用數(shù)值方法來研究其中的流動問題。如Kaviany[1]采用有限差分法模擬了等溫平板間多孔通道中的層流流動，通過引入多孔介質(zhì)形狀參數(shù)γ，分析了γ在極限范圍(0，∞)內(nèi)對多孔介質(zhì)通道流體的速度分布及壓力損失的影響。Nithiarasu等[2]采用有限元法模擬了變孔隙率條件下多孔介質(zhì)中流體的等溫流動，利用半隱式時間離散格式并結(jié)合速度修正方法推導(dǎo)出瞬態(tài)多孔介質(zhì)中流場的動量方程。楊勃等[3]采用有限容積法模擬多孔介質(zhì)通道內(nèi)的流體流動過程，研究了不同孔隙率下流場的變化規(guī)律。常煥靜[4]采用有限差分法分別模擬了水和冪律流體在多孔介質(zhì)中的流動，并分析了多孔介質(zhì)不同的排列方式及顆粒直徑對流體的壓降和阻力系數(shù)的影響。需要指出的是，有限差分法、有限元法和有限容積法等經(jīng)典數(shù)值方法雖然在針對多孔介質(zhì)內(nèi)流體流動特性的研究中運用比較廣泛，但它們均為低階收斂方法，相比之下，譜方法具有指數(shù)級收斂精度且易實施、穩(wěn)定性好[5]，因而應(yīng)用前景廣闊。本課題組基于Chebyshev-Fourier配置點譜方法已穩(wěn)定地求解了圓柱坐標(biāo)系下三維瞬態(tài)Navier-Stokes方程[6]，并采用Chebyshev配置點譜方法成功模擬了方形多孔腔內(nèi)的自然對流問題[7-9]。在數(shù)值研究中，往往易忽略模型入口段長度和邊界層厚度，導(dǎo)致速度分布和傳熱計算產(chǎn)生嚴(yán)重誤差，因此，在本課題組前期研究基礎(chǔ)上，本文基于Chebyshev配置點譜方法借助MATLAB模擬多孔介質(zhì)平板通道內(nèi)的流體流動特性，并對程序可靠性進行驗證，重點探討達西數(shù)(Da)、雷諾數(shù)(Re)及孔隙率(ε)對流體流動速度分布、邊界層厚度及入口長度的影響。

1 物理與數(shù)學(xué)模型

以典型的多孔介質(zhì)平板通道流動模型為研究對象，圖1所示為2維軸對稱多孔介質(zhì)平板通道流動物理模型，其高度為H/2，長度L為2H，流體以速度u(y)進入通道，平均速度為u0。

·V=0

(1)

(2)

τ=0,0≤X≤2,0≤Y≤0.5∶U=V=0

(3)

(4)

其它邊界條件如壓力等見計算方法部分。

2 計算方法

2.1 時間離散

首先采用準(zhǔn)隱式格式[12]對動量方程(式(2))進行時間離散，其中時間導(dǎo)數(shù)項采用二階向后差分格式，壓力項、擴散項和源項均采用隱式格式，對流項采用顯式Adams-Bashforth格式，則式(2)可離散為

(5)

式中：τs為無量綱時間步長；n-1、n和n+1表示時間離散級數(shù)。

2.2 改進的投影算法

采用由Hugues等[13]提出的IPS(improved projection scheme)處理速度場和壓力場的耦合求解問題，執(zhí)行步驟如下：

(1)壓力預(yù)估步，求解預(yù)估壓力P*。對式(5)兩邊取散度并用連續(xù)性條件限制，得到預(yù)估壓力泊松方程

ΔP*=

(6)

(7)

式中：n表示邊界的法向量。

為了提高穩(wěn)定性，將拉普拉斯算子分解為一個螺旋矢量場和一個無旋矢量場[14]，并用連續(xù)性條件限制后得到

2ΔVn-ΔVn-1=-2××Vn+××Vn-1

(8)

(2)速度預(yù)估步，求解預(yù)估速度V*。根據(jù)P*和式(5)，獲得亥姆霍茲方程

(9)

·Vn+1=0

(10)

(11)

(12)

Vn+1=V*-

(13)

(14)

2.3 空間離散

(15)

(16)

(17)

(18)

式中：T表示矩陣轉(zhuǎn)置。

2.4 系數(shù)矩陣的確定

AΦ+ΦB=F

(19)

式中：A、B分別為(Ny+1)、(Nx+1)階系數(shù)矩陣；Φ、F為(Ny+1)×(Nx+1)系數(shù)矩陣，Φ表示φ(Xi,Yj)的值。事實上，未知邊界向量可以表示為已知邊界值和函數(shù)內(nèi)部節(jié)點值φ(Xi,Yj)(i=1,2,…,Nx-1，j=1,2,…,Ny-1)的函數(shù)關(guān)系式，將矩陣方程(19)改寫為

(20)

(21)

式中：hx+,j、gx-,j、hy+,i、hy-,i均為已知邊界函數(shù)。然后根據(jù)離散邊界條件和已知邊界值可獲得未知邊界值表達式

(22)

1≤j,k≤Ny-1

(23)

(24)

1≤i≤Nx-1，1≤j≤Ny-1

(25)

2.5 計算步驟

多孔介質(zhì)平板通道流動問題求解的計算步驟如下：

(2)給無量綱變量及各參數(shù)賦值，初始化速度場、壓力場和中間變量場，并令φ-1=φ0；

(4)設(shè)置收斂標(biāo)準(zhǔn)，開始時間步進；

(8)模擬結(jié)果如果滿足收斂標(biāo)準(zhǔn)就終止循環(huán)，否則，推進一個時間步，令φn-1=φn、φn=φn+1，返回步驟(5)；

(9)輸出結(jié)果。

計算中，取無量綱時間步長Δτ=1×10-6，無量綱速度和無量綱壓力的收斂標(biāo)準(zhǔn)分別為‖Un+1-Un‖max<1×10-6、‖Vn+1-Vn‖max<1×10-6、‖Pn+1-Pn‖max<1×10-6。

2.6 網(wǎng)格無關(guān)性檢驗

通過網(wǎng)格加密，考察距離對稱軸最近的一條水平網(wǎng)格線上U(Xi,Y1)的分布曲線變化趨勢以驗證網(wǎng)格無關(guān)性。x、y方向的網(wǎng)格節(jié)點數(shù)(Nx+1)×(Ny+1)從13×4增加到73×19。計算中取ε=0.9、Da=1×10-4、Re=100，檢驗結(jié)果如圖2所示。從圖2中可以看出，當(dāng)網(wǎng)格節(jié)點數(shù)從13×4增加到53×14時，無量綱速度U(Xi,Y1)變化較為明顯，繼續(xù)增加網(wǎng)格節(jié)點數(shù)，無量綱速度U(Xi,Y1)幾乎不再變化，表明網(wǎng)格節(jié)點數(shù)為53×14時數(shù)值解滿足網(wǎng)格無關(guān)性要求，故而，在本文后續(xù)計算中x、y方向的網(wǎng)格節(jié)點數(shù)取53×14。

圖2 網(wǎng)格無關(guān)性檢驗

2.7 程序驗證

不同γ條件下，流體在多孔介質(zhì)平板通道出口速度分布以及在充分發(fā)展區(qū)域最大速度的精確解表達式[1]分別為

U|X=2=[1-e-2γ-(1-e-γ)(eγ(Yj-1)+e-γYj)]×

[1-e-2γ-2(1-e-γ)2γ-1]-1

(26)

(27)

(a) U

(b) Umax

3 結(jié)果與分析

3.1 流體速度分布

當(dāng)ε為0.9、Da為1×10-4、Re為100時，流體在多孔介質(zhì)平板通道內(nèi)沿來流方向不同截面處的無量綱速度分布如圖4所示。由圖4可見，沿著來流方向，由于多孔介質(zhì)的彌散作用，近壁面處的速度梯度迅速增大，當(dāng)X趨近于0.590時，除貼近壁面的極短區(qū)域外，其它區(qū)域的流體速度分布接近于常數(shù)值。

圖4 流體速度分布

3.2 Da對多孔介質(zhì)平板通道內(nèi)流體流動的影響

當(dāng)ε為0.9、Re為100時，在不同Da條件下，多孔介質(zhì)平板通道中的流體在出口(X=2)和對稱軸(Y=0)上的無量綱速度分布如圖5所示。根據(jù)圖5模擬結(jié)果，可求得不同Da條件下，多孔介質(zhì)平板通道流體流速在充分發(fā)展區(qū)域出口處的無量綱邊界層厚度δ和入口長度xe如表1所示。結(jié)合圖5和表1可知，隨著Da的減小，近壁面處流體速度梯度逐漸增大，無量綱邊界層厚度δ及入口長度xe不斷減小，但當(dāng)Da<1×10-4時，繼續(xù)減小Da，δ及xe均不再有明顯變化。

(a) 出口

(b) 對稱軸

表1 不同Da時多孔介質(zhì)平板通道流邊界層厚度和入口長度

3.3 Re對多孔介質(zhì)平板通道內(nèi)流體流動的影響

當(dāng)ε為0.9、Da為0.01時，在不同Re條件下，多孔介質(zhì)平板通道中的流體在出口和對稱軸上的無量綱速度分布如圖6所示，并根據(jù)圖6模擬結(jié)果，求得相應(yīng)條件下多孔介質(zhì)平板通道流在充分發(fā)展區(qū)域出口處的速度邊界層厚度δ和入口段長度xe如表2所示。結(jié)合圖6和表2可以看出，隨著Re的增大，多孔介質(zhì)平板通道流體在充分發(fā)展區(qū)域出口處的速度邊界層不斷變薄，而入口段長度相應(yīng)增大。

(a) 出口

(b) 對稱軸

表2 不同Re時多孔介質(zhì)平板通道流的邊界層厚度和入口長度

3.4 ε對多孔介質(zhì)平板通道內(nèi)流體流動的影響

當(dāng)Da為0.01、Re為10時，在不同ε條件下，多孔介質(zhì)平板通道中的流體在出口和對稱軸上的無量綱速度分布如圖7所示，相應(yīng)條件下多孔介質(zhì)平板通道流體在充分發(fā)展區(qū)域出口處的速度邊界層厚度δ和入口段長度xe的計算結(jié)果列于表3。結(jié)合圖7和表3可知，隨著ε的減小，多孔介質(zhì)平板通道流體在充分發(fā)展區(qū)域出口處的速度邊界層厚度及入口段長度均呈現(xiàn)出增加趨勢。

(a) 出口

(b) 對稱軸

表3 不同ε時多孔介質(zhì)平板通道流邊界層厚度和入口長度

4 結(jié)語

本文基于Chebyshev配置點譜方法，借助MATLAB利用數(shù)值模擬分析了多孔介質(zhì)平板通道內(nèi)流體流動的問題。針對動量方程的離散，空間上和時間上分別采用Chebyshev配置點譜方法和準(zhǔn)隱式格式離散，同時采用IPS處理速度場和壓力場的耦合求解問題。將數(shù)值模擬結(jié)果與精確解進行比較，發(fā)現(xiàn)兩者具有很好的一致性，這證實了使用Chebyshev配置點譜方法求解該類問題有效可行。同時，數(shù)值模擬還發(fā)現(xiàn)，達西數(shù)(Da)、雷諾數(shù)(Re)及孔隙率(ε)對多孔介質(zhì)平板通道流體在充分發(fā)展區(qū)域的速度邊界層厚度和入口長度具有明顯影響作用：當(dāng)其它條件確定時，速度邊界層厚度和入口長度隨著Da的減小而減小，但當(dāng)Da減小到1×10-4以下時，這些變化不再明顯；隨著Re的增大，速度邊界層厚度不斷變薄，入口長度逐漸增大；速度邊界層厚度和入口長度隨著ε的減小均呈現(xiàn)出增大趨勢。本文研究結(jié)果為今后進一步研究多孔介質(zhì)內(nèi)流體的流動、換熱與燃燒問題提供了參考。