基于MATLAB 的傳熱學課程虛擬仿真實驗平臺設計

張程賓，韓 群，陳永平,2

（1. 東南大學 能源與環境學院，江蘇 南京 210096；2. 蘇州科技大學 環境科學與工程學院，江蘇 蘇州 215009）

傳熱學是一門研究由溫差引起的熱能傳遞規律的科學。傳熱學的研究方法主要有實驗測定、理論分析和數值模擬[1]，其中測定實驗教學因受經費、教學學時和實驗儀器的限制，傳熱實驗在可控和實驗現象可視化方面往往不能滿足教學要求。由于實際傳熱問題的復雜性，用理論分析方法也很難得出能量方程這一類偏微分方程的解析解[2]。數值模擬方法因具有教學成本低、可視化效果好、便于學生理解和運用等優點，在傳熱學課程教學上得到了重視。

在對傳熱問題的數值模擬軟件中，MATLAB 具有顯著的數值分析和圖形處理能力，適用于求解傳熱問題的偏微分方程[3]。為此，本文采用MATLAB 軟件對一維和二維熱傳導問題的有限差分法進行求解，通過MATLAB 圖形用戶界面（graphics user interface，GUI）開發了傳熱問題數值求解的虛擬實驗軟件，對熱傳導問題實現數據輸入、回調運算、數據輸出和結果顯示功能。通過對不同條件下一維和二維熱傳導問題進行數值模擬，強化了學生在傳熱學課程上對該知識點的學習，提高了傳熱學課程的可視化教學效果。

1 傳熱方程的數值求解

傳熱問題數值求解的基本思路，是將導熱物體在時間和空間上連續的溫度場用有限個離散點的值代替，并采用有限差分和有限體積等數值方法建立離散值集合的代數方程，通過求解代數方程來獲得離散點上溫度的值[4-5]。如圖1 所示，用與坐標軸平行的網格線劃分求解區域，網格線間的交點稱為節點，相鄰節點之間的距離稱為步長。

圖1 網格劃分

1.1 建立數學模型和區域離散

二維非穩態導熱微分方程的一般形式為[6]：

式中，ρ為微元體的密度、c為微元體的比熱容、λ為微元體的導熱系數、τ為時間、t為溫度、φ為源項，x和y表示水平和豎直方向的節點坐標。當左端時間項設置為0 時，為穩態導熱微分方程；只考慮x方向的熱傳導時，為一維非穩態的導熱微分方程。對于內部節點離散方程，可以采用泰勒級數展開法進行建立：

相應的代數方程求解格式分別為顯格式和隱格式。

取水平和豎直的網格步長Δx=Δy，對于顯格式求解方法，聯立方程（1）、（2）、（3）、（4），可得內節點的離散方程

同理，對于隱格式求解方法，聯立方程（1）、（2）、（3）、（5）可得內節點離散方程：

1.2 邊界條件處理

由于含有未知的邊界溫度，內節點建立的代數方程組是不封閉的，需要給定邊界節點方程。常見的導熱問題有3 類邊界條件：規定邊界的溫度（第一類邊界條件）、規定邊界的熱流密度（第二類邊界條件）和規定表面傳熱系數h及周圍流體的溫度tf（第三類邊界條件）[7]。

邊界上的節點（m,n）如圖2 所示。由能量守恒定律可得離散方程組

圖2 邊界上的節點

當邊界條件為規定溫度值時，可以對方程組直接求解并得到溫度場的分布。當給定的邊界條件為第二類或第三類邊界條件時，需要分別進行討論：

（1）第二類邊界條件：規定傳入計算區域的熱量為正，此時將給定的熱流密度值q代入式（8）進行求解；當給定的是絕熱邊界時，可令式中

2 編程依據的算法流程

如圖 3 所示，基于 MATLAB 的傳熱學課程虛擬仿真實驗平臺的軟件算法[8-9]求解基本流程如圖3 所示。

圖3 求解基本流程圖

（1）根據實際物理問題建立數學控制方程，例如在一維非穩態導熱問題中，左端為恒定熱流，右端為與外界環境對流換熱求解，數值求解不同時間下的溫度分布。

（2）確定邊界條件和初始條件，在GUI 界面上輸入初始溫度、熱擴散率、環境溫度、長度和時間等參數數據[10-12]。

（3）劃分網格并對區域進行離散化，建立代數方程。

（4）設立初值并進行迭代求解，當滿足收斂要求時將顯示計算結果，如圖4 所示。

圖4 仿真實驗計算結果

3 功能模塊介紹

傳熱學課程的虛擬仿真實驗平臺分為2 個模塊：一維區域傳熱模塊和二維區域傳熱模塊，每個模塊又根據導熱類型分為穩態傳熱和非穩態傳熱，根據邊界條件分為3 類熱邊界條件（恒定溫度、恒定熱流和對流換熱）。在各個模塊下設立彈出式選單，通過模塊的選項在面板中設置參數，從而進入選定的模型中進行虛擬仿真模型計算。

3.1 一維區域傳熱模塊

一維區域傳熱模塊包括一維穩態傳熱模塊和一維非穩態傳熱模塊，在選定區域和導熱類型后，設置熱邊界條件、物性參數和空間區域(非穩態需要劃分時間)，再通過點擊“計算”按鈕以對一維傳熱問題進行求解，計算的溫度分布圖顯示在坐標軸上[13-14]。

3.1.1 數學模型

以一維非穩態問題為例，取邊界條件為恒定溫度的一維非穩態導熱，相應的數學表達式如式（9）所示：

在數值計算中需要采用的物理參數如表1 所示。

表1 一維非穩態導熱問題的物理參數表

對于本問題的解析解為：

式中a為熱擴散系數：

取時間t=0.1 s 時的溫度分布進行數值計算結果與精確解析解的對比。如圖5 所示，二者的結果非常接近，表明了數值計算結果的正確性。

圖5 t=0.1 s 時數值解與分析解對比圖

3.1.2 GUI 運行結果與分析

在數值計算過程中，分別設定時間為 0.005 s、0.025 s、0.05 s、0.1 s 和 100 s。在 GUI 的坐標軸上溫度分布圖如圖6 所示。

圖6 一維區域非穩態傳熱的溫度分布

選取每個設定時間最后的運行結果進行對比分析，其結果如圖7 所示。仿真結果表明：在一維非穩態導熱的初期，受邊界低溫的影響，物體靠近邊界的溫度較低，中間段的溫度較高。隨著時間的推移，非穩態導熱過程繼續進行，物體的溫度持續下降，最終物體內部各處的溫度無限接近于給定的邊界溫度。

圖7 不同時間的數值計算結果

通過本模塊的演示，學生可以學習到一維穩態和非穩態傳熱的基本規律，認識到一維區域內溫度的演化過程，掌握邊界條件和物性參數的因素對傳熱效果的影響。

3.2 二維區域傳熱模塊

二維區域傳熱模塊包括二維穩態傳熱模塊和二維非穩態傳熱模塊，相關的GUI 數值模擬計算步驟與一維區域傳熱模塊的計算步驟類似：選定區域和導熱類型，設置熱邊界條件、物性參數和空間區域（非穩態需要劃分時間），再點擊“計算”按鈕以對二維傳熱問題進行求解[15-16]。

3.2.1 數學模型

以二維非穩態傳熱問題為例，二維區域傳熱模塊的功能的物理量參數如表2 所示。

表2 數值計算中的物理量參數表

3.2.2 GUI 運行結果與分析

在 GUI 的時間劃分面板上，分別設定時間t為2.5×10-4s、5×10-4s、0.001 s、0.005 s、0.025 s 和 0.05 s進行仿真。二維區域非穩態傳熱的溫度分布模擬結果如圖8 所示。圖中矩形區域內的溫度分布隨不同時間設定變化，說明溫度的傳遞是從左端逐漸向內部傳遞，并通過一定時間的熱傳導，溫度場的分布達到穩定狀態。這是因為在矩形二維非穩態導熱問題中，右邊界設置為第三類邊界條件（固體表面與周圍流體間換熱），下表面設置為絕熱邊界，上邊界和左邊界為第二類邊界條件（恒定熱流密度）且左邊界熱流遠大于上邊界熱流，所以二維區域的熱源輸入主要由左端提供。在整個傳熱過程中，越靠近左邊界的溫度越高。

圖8 非穩態傳熱條件下矩形區域的溫度分布

隨著時間的推移，在矩形區域內通過熱傳導方式傳遞熱量，溫度由左邊界向內部傳遞，最后達到穩定的溫度分布。

通過本模塊的學習，學生可以仔細觀察二維條件下溫度場演化過程，掌握邊界條件對導熱的影響規律。

4 結語

基于 MATLAB 的傳熱學課程虛擬仿真實驗平臺實現了一維和二維區域傳熱的穩態/非穩態傳熱數值模擬和結果可視化動態演示，通過與解析解對比，驗證了該虛擬仿真實驗平臺運算的可靠性。

虛擬仿真傳熱學實驗具有可視化效果好、操作簡單、價格低廉、運算準確、高效的優勢。將虛擬仿真實驗軟件應用于傳熱學的課程中，有助于學生對傳熱問題的數值求解思路和方法的理解，同時也學習了MATLAB 在數值求解上的編程應用和 GUI 的設計，具有較高的教學實用性。