龐潤(rùn)芳，鄭坤燦

(1.內(nèi)蒙古科技大學(xué)工程訓(xùn)練中心，內(nèi)蒙古包頭 014010；2.內(nèi)蒙古科技大學(xué)能源與環(huán)境工程學(xué)院，內(nèi)蒙古包頭 014010)

0 引言

高溫散體廣泛存在于化工、環(huán)保、冶金、建材、礦業(yè)等許多工業(yè)過程，而且在航空航天、核反應(yīng)堆等領(lǐng)域也有著廣泛應(yīng)用，因此高溫散體的高效熱回收和強(qiáng)化傳熱問題近年來一直是多孔介質(zhì)對(duì)流換熱研究的熱點(diǎn)和難點(diǎn)。人們對(duì)高溫散體對(duì)流換熱已經(jīng)進(jìn)行了一定的實(shí)驗(yàn)和數(shù)值模擬研究：馬淑杰等[1]建立了燒結(jié)礦冷卻過程的多孔介質(zhì)模型，探討了燒結(jié)礦熱物性參數(shù)分別采用經(jīng)驗(yàn)值和實(shí)驗(yàn)值對(duì)數(shù)值模擬的影響；夏建芳等[2]應(yīng)用VB.NET語言,基于CFD平臺(tái),開發(fā)了環(huán)式鼓風(fēng)冷卻機(jī)燒結(jié)礦冷卻過程自動(dòng)仿真軟件；劉立鈞[3]采用有限差分法對(duì)燒結(jié)礦豎冷窯內(nèi)氣-固換熱過程進(jìn)行了數(shù)值模擬研究；倪文杰[4]將富氧和噴吹焦?fàn)t煤氣與煙氣循環(huán)聯(lián)合使用,形成綜合燒結(jié)工藝,通過數(shù)值模擬研究煙氣循環(huán)條件下富氧和焦?fàn)t煤氣噴吹參數(shù)對(duì)燒結(jié)工藝和產(chǎn)品質(zhì)量的影響;劉瑋寅[5]等通過仿真模擬軟件COMSOL Multiphysics 對(duì)燒結(jié)礦豎罐式冷卻過程進(jìn)行數(shù)值模擬研究；馮軍勝[6]等利用Fluent 軟件及其二次開發(fā)平臺(tái)對(duì)燒結(jié)礦余熱回收豎罐內(nèi)氣固傳熱過程進(jìn)行了數(shù)值分析；果晶晶[7]等采用CFD 軟件對(duì)燒結(jié)礦余熱罐內(nèi)的氣固換熱過程進(jìn)行了數(shù)值模擬研究;馮軍勝[8]等借助Fluent模擬計(jì)算平臺(tái),利用正交實(shí)驗(yàn)方法對(duì)燒結(jié)礦豎罐內(nèi)氣固傳熱過程進(jìn)行了數(shù)值模擬與優(yōu)化；吳禮忠[9]等基于Fluent 模擬軟件,通過使用UDS 和UDF 分別對(duì)函數(shù)方程和源項(xiàng)進(jìn)行定義與修改,采用多孔介質(zhì)模型對(duì)燒結(jié)礦層冷卻過程進(jìn)行模擬。本文選擇了Matlab 語言來實(shí)現(xiàn)對(duì)高溫散體對(duì)流換熱過程的數(shù)值模擬。

Matlab 最早是專門針對(duì)矩陣簡(jiǎn)化運(yùn)算開發(fā)的，所以被稱為矩陣實(shí)驗(yàn)室(Matrix Laboratory) 。該軟件具有強(qiáng)大高效的多維數(shù)據(jù)處理能力，可以進(jìn)行數(shù)值分析、矩陣計(jì)算、科學(xué)數(shù)據(jù)可視化和非線性動(dòng)態(tài)系統(tǒng)的建模與仿真等，具有代碼簡(jiǎn)潔、編程效率高等特點(diǎn),因此比其他同類軟件操作簡(jiǎn)單方便[10-11]。同時(shí)軟件還具有巨量的庫函數(shù)，而且也可以方便定義自己庫函數(shù)。

高溫散體對(duì)流換熱模型的數(shù)值模擬的基本步驟包括物理數(shù)學(xué)模型的建立、空間和數(shù)學(xué)方程的離散差分、代數(shù)方程的求解、計(jì)算機(jī)編程和結(jié)果分析處理。

1 高溫散體對(duì)流換熱模型的建立

1.1 高溫散體冷卻過程物理模型

目前工業(yè)高溫散體對(duì)流換熱常見工藝有氣固逆流和氣固錯(cuò)流兩大類，一般固體為緩慢移動(dòng)，本文以后者作為研究對(duì)象。假設(shè)某高溫散體，溫度在300～1 500℃之間，在冷空氣中緩慢移動(dòng)冷卻直到降到某一特定溫度下，冷卻空氣從高溫散體下部垂直吹入散體料層，從上部進(jìn)入煙罩或排放。整個(gè)過程可以簡(jiǎn)化為如圖1所示的物理模型。

圖1 高溫散體冷卻過程物理模型

1.2 高溫散體對(duì)流換熱數(shù)學(xué)模型

針對(duì)圖1 所示的高溫散體冷卻過程，在料層較寬的情況下忽略邊壁效應(yīng)、寬度維數(shù)、長(zhǎng)度維數(shù)，再假設(shè)高溫散體由大小不同的球形顆粒隨機(jī)堆積而成且各向同性，忽略對(duì)環(huán)境熱損失、氣體導(dǎo)熱和輻射傳熱，高溫散體三維數(shù)學(xué)模型就可以簡(jiǎn)化為一維非穩(wěn)態(tài)模型，如式(1)到式(6)所示。

1)連續(xù)性方程

2)能量方程

①氣體能量方程

結(jié)合連續(xù)性方程可以簡(jiǎn)化為：

②固體能量方程

③動(dòng)量方程

高溫散體多采用Ergun公式計(jì)算，即：

式中：ΔP—高溫散體內(nèi)部壓降，Pa；

H—對(duì)應(yīng)壓降的高溫散體高度，m；

de—當(dāng)量直徑，m。

④狀態(tài)方程

式中：N—?dú)怏w摩爾數(shù)，mol；

R—?dú)怏w常數(shù)，8.314J/(mol·K)；

M—?dú)怏w摩爾質(zhì)量，kg/mol。

式(5-24)也可以寫為另一種形式：PM=ρRT。

3)定解條件

初始時(shí)刻(t=0)，料層入口處初始料溫為Ts0，空氣入口處(z=0)，空氣入口速度ug0.，空氣溫度為環(huán)境溫度Tg0。

2 一維非穩(wěn)態(tài)的數(shù)值求解

數(shù)值求解就是將時(shí)空離散化，把數(shù)學(xué)微分方程組轉(zhuǎn)換為時(shí)空節(jié)點(diǎn)上的代數(shù)方程，最后解出所有節(jié)點(diǎn)的代數(shù)方程組，得到相應(yīng)的溫度時(shí)空分布。

2.1 離散代數(shù)方程組

數(shù)學(xué)方程離散方法有很多，而目前傳熱過程主要采取有限容積法，內(nèi)部節(jié)點(diǎn)用有限差分，邊角節(jié)點(diǎn)用熱平衡法。針對(duì)一維非穩(wěn)態(tài)模型，為了保證差分方程的穩(wěn)定，采用隱式差分格式。

1)連續(xù)方程離散

式中：n—時(shí)間節(jié)點(diǎn)編號(hào)；

k—空間節(jié)點(diǎn)編號(hào)。

2)氣體能量方程離散

a)內(nèi)部節(jié)點(diǎn)

b)入口邊界點(diǎn)：k=0，Tn,k=Tg0，Tg0為環(huán)境空氣溫度。

c) 出口邊界點(diǎn)：k=kmax，根據(jù)能量平衡，忽略熱損失，方程為：

3)固體能量方程離散

①內(nèi)部節(jié)點(diǎn)

②入口邊界點(diǎn)：k=0，根據(jù)能量平衡，忽略熱損失，方程為：

式中，k=0時(shí)，T0,k=Ts0，Ts0為入料溫度。

③出口邊界點(diǎn)：k=kmax，根據(jù)能量平衡，忽略熱損失，方程為：

3 數(shù)值求解思路及Matlab編程實(shí)現(xiàn)

3.1 數(shù)值求解思路

數(shù)值求解主要包括前處理（網(wǎng)格劃分）、代數(shù)方程求解和后處理（結(jié)果分析處理），總程序框圖如圖2所示。

圖2 程序設(shè)計(jì)總框圖

3.2 數(shù)值求解編程實(shí)現(xiàn)

1)網(wǎng)格和變量初始化

劃分時(shí)間網(wǎng)格和空間網(wǎng)格，時(shí)間用t表示，空間高度用z表示，節(jié)點(diǎn)編號(hào)為1、2、3……設(shè)置空氣、散體物性參數(shù)、運(yùn)動(dòng)參數(shù)和散體結(jié)構(gòu)參數(shù)等，如：

2)求解各時(shí)刻各空間節(jié)點(diǎn)溫度和速度

①求解初始時(shí)刻各空間節(jié)點(diǎn)溫度和速度

由于定解條件只知道t=0時(shí)的固體溫度和z=0處的氣體溫度，初始時(shí)刻z≠0處的氣體溫度未知，所以先要假設(shè)初始溫度場(chǎng)和壓力場(chǎng)，進(jìn)而算出速度場(chǎng)，采用迭代求解值反復(fù)修改初始溫度場(chǎng)和壓力場(chǎng)，直至溫度場(chǎng)和壓力場(chǎng)穩(wěn)定，此時(shí)即得到初始時(shí)刻的真實(shí)溫度，在此基礎(chǔ)上就可以計(jì)算以后各個(gè)時(shí)刻的溫度場(chǎng)了。該求解過程具體如圖3所示。

圖3 初始時(shí)刻溫度計(jì)算流程圖

②其余時(shí)刻各空間節(jié)點(diǎn)溫度的計(jì)算

初始時(shí)刻的溫度、速度和壓力真值求出后，根據(jù)狀態(tài)方程、初始?jí)毫?、初始溫度求出初始第一?jié)點(diǎn)的密度，再結(jié)合初始速度利用Ergun 方程求出下一個(gè)節(jié)點(diǎn)的壓力，該壓力再通過狀態(tài)方程可以得到該節(jié)點(diǎn)的密度和速度，于是利用TDMA追趕法求出該時(shí)刻該節(jié)點(diǎn)的溫度。就這樣反復(fù)采用該方法就可以求出下一個(gè)空間節(jié)點(diǎn)，下下個(gè)節(jié)點(diǎn)直至所有空間節(jié)點(diǎn)的溫度、壓力、速度和溫度。此處要注意的是物性均隨溫度而變化，在該節(jié)點(diǎn)溫度未求出時(shí)物性是按上一時(shí)刻的溫度來計(jì)算的，所以誤差就會(huì)產(chǎn)生，尤其是時(shí)間間隔取得較大時(shí)更為明顯。改變的方法是再加一層循環(huán)，即把現(xiàn)在計(jì)算出的溫度再代回去重新計(jì)算，直到新舊溫度間的差值小于指定精度為止。

部分代碼如下：

4 結(jié)果分析處理

對(duì)結(jié)果的分析處理一般稱為后處理，采用圖形顯示更為直觀，最后可以輸出不同位置處料溫變化和空氣溫度變化如圖4 和圖5 所示，也可以輸出燒結(jié)礦和空氣的時(shí)空溫度場(chǎng)如圖6 和圖7 所示，還可以結(jié)合實(shí)驗(yàn)數(shù)據(jù)輸出對(duì)比驗(yàn)證結(jié)果如圖8所示等。

圖4 燒結(jié)礦不同高度溫度沿環(huán)冷機(jī)周向方向上的變化

圖5 燒結(jié)礦不同高度處風(fēng)溫變化規(guī)律

圖6 環(huán)冷機(jī)內(nèi)部燒結(jié)礦的溫度場(chǎng)

圖7 環(huán)冷機(jī)內(nèi)部冷卻空氣的溫度場(chǎng)

圖8 計(jì)算數(shù)據(jù)與工業(yè)實(shí)測(cè)數(shù)據(jù)比較（湍流關(guān)聯(lián)式）

固體溫度變化代碼（氣體代碼略）：

legend('Zukauskas 計(jì)算值','實(shí)測(cè)值')%,'Wakao and Kaguei','Kreith,Frank','Zukauskas','分形self',2)%,'nu7');,'nu6'

5 結(jié)論

論文利用Matlab 編寫了高溫散體對(duì)流換熱通用計(jì)算程序，通過TDMA 追趕法求解了氣流速度、氣體溫度和散體溫度的時(shí)空離散解，實(shí)現(xiàn)了氣固換熱過程流動(dòng)與傳熱的數(shù)值模擬。最后利用Matlab 強(qiáng)大的圖形處理功能對(duì)數(shù)值模擬結(jié)果進(jìn)行了可視化處理，獲得了不同高度的風(fēng)溫和料溫變化規(guī)律，繪制了全場(chǎng)的風(fēng)溫和料溫分布圖，展示了模擬計(jì)算結(jié)果和實(shí)驗(yàn)值之間的差別。