高溫散體對流換熱數值模擬Matlab編程實現

龐潤芳，鄭坤燦

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

0 引言

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

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

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

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

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

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

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

1.2 高溫散體對流換熱數學模型

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

1)連續性方程

2)能量方程

①氣體能量方程

結合連續性方程可以簡化為：

②固體能量方程

③動量方程

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

式中：ΔP—高溫散體內部壓降，Pa；

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

de—當量直徑，m。

④狀態方程

式中：N—氣體摩爾數，mol；

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

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

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

3)定解條件

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

2 一維非穩態的數值求解

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

2.1 離散代數方程組

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

1)連續方程離散

式中：n—時間節點編號；

k—空間節點編號。

2)氣體能量方程離散

a)內部節點

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

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

3)固體能量方程離散

①內部節點

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

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

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

3 數值求解思路及Matlab編程實現

3.1 數值求解思路

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

圖2 程序設計總框圖

3.2 數值求解編程實現

1)網格和變量初始化

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

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

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

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

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

②其余時刻各空間節點溫度的計算

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

部分代碼如下：

4 結果分析處理

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

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

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

圖6 環冷機內部燒結礦的溫度場

圖7 環冷機內部冷卻空氣的溫度場

圖8 計算數據與工業實測數據比較（湍流關聯式）

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

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

5 結論

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