高溫散體對流換熱數(shù)值模擬Matlab編程實現(xiàn)

龐潤芳，鄭坤燦

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

0 引言

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

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

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

1 高溫散體對流換熱模型的建立

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

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

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

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

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

1)連續(xù)性方程

2)能量方程

①氣體能量方程

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

②固體能量方程

③動量方程

高溫散體多采用Ergun公式計算，即：

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

H—對應壓降的高溫散體高度，m；

de—當量直徑，m。

④狀態(tài)方程

式中：N—氣體摩爾數(shù)，mol；

R—氣體常數(shù)，8.314J/(mol·K)；

M—氣體摩爾質量，kg/mol。

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

3)定解條件

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

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

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

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

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

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

式中：n—時間節(jié)點編號；

k—空間節(jié)點編號。

2)氣體能量方程離散

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

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

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

3)固體能量方程離散

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

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

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

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

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

3.1 數(shù)值求解思路

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

圖2 程序設計總框圖

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

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

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

2)求解各時刻各空間節(jié)點溫度和速度

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

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

圖3 初始時刻溫度計算流程圖

②其余時刻各空間節(jié)點溫度的計算

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

部分代碼如下：

4 結果分析處理

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

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

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

圖6 環(huán)冷機內(nèi)部燒結礦的溫度場

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

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

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

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

5 結論

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