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

（1. 東南大學(xué) 能源與環(huán)境學(xué)院，江蘇 南京 210096；2. 蘇州科技大學(xué) 環(huán)境科學(xué)與工程學(xué)院，江蘇 蘇州 215009）

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

在對(duì)傳熱問(wèn)題的數(shù)值模擬軟件中，MATLAB 具有顯著的數(shù)值分析和圖形處理能力，適用于求解傳熱問(wèn)題的偏微分方程[3]。為此，本文采用MATLAB 軟件對(duì)一維和二維熱傳導(dǎo)問(wèn)題的有限差分法進(jìn)行求解，通過(guò)MATLAB 圖形用戶(hù)界面（graphics user interface，GUI）開(kāi)發(fā)了傳熱問(wèn)題數(shù)值求解的虛擬實(shí)驗(yàn)軟件，對(duì)熱傳導(dǎo)問(wèn)題實(shí)現(xiàn)數(shù)據(jù)輸入、回調(diào)運(yùn)算、數(shù)據(jù)輸出和結(jié)果顯示功能。通過(guò)對(duì)不同條件下一維和二維熱傳導(dǎo)問(wèn)題進(jìn)行數(shù)值模擬，強(qiáng)化了學(xué)生在傳熱學(xué)課程上對(duì)該知識(shí)點(diǎn)的學(xué)習(xí)，提高了傳熱學(xué)課程的可視化教學(xué)效果。

1 傳熱方程的數(shù)值求解

傳熱問(wèn)題數(shù)值求解的基本思路，是將導(dǎo)熱物體在時(shí)間和空間上連續(xù)的溫度場(chǎng)用有限個(gè)離散點(diǎn)的值代替，并采用有限差分和有限體積等數(shù)值方法建立離散值集合的代數(shù)方程，通過(guò)求解代數(shù)方程來(lái)獲得離散點(diǎn)上溫度的值[4-5]。如圖1 所示，用與坐標(biāo)軸平行的網(wǎng)格線劃分求解區(qū)域，網(wǎng)格線間的交點(diǎn)稱(chēng)為節(jié)點(diǎn)，相鄰節(jié)點(diǎn)之間的距離稱(chēng)為步長(zhǎng)。

圖1 網(wǎng)格劃分

1.1 建立數(shù)學(xué)模型和區(qū)域離散

二維非穩(wěn)態(tài)導(dǎo)熱微分方程的一般形式為[6]：

式中，ρ為微元體的密度、c為微元體的比熱容、λ為微元體的導(dǎo)熱系數(shù)、τ為時(shí)間、t為溫度、φ為源項(xiàng)，x和y表示水平和豎直方向的節(jié)點(diǎn)坐標(biāo)。當(dāng)左端時(shí)間項(xiàng)設(shè)置為0 時(shí)，為穩(wěn)態(tài)導(dǎo)熱微分方程；只考慮x方向的熱傳導(dǎo)時(shí)，為一維非穩(wěn)態(tài)的導(dǎo)熱微分方程。對(duì)于內(nèi)部節(jié)點(diǎn)離散方程，可以采用泰勒級(jí)數(shù)展開(kāi)法進(jìn)行建立：

相應(yīng)的代數(shù)方程求解格式分別為顯格式和隱格式。

取水平和豎直的網(wǎng)格步長(zhǎng)Δx=Δy，對(duì)于顯格式求解方法，聯(lián)立方程（1）、（2）、（3）、（4），可得內(nèi)節(jié)點(diǎn)的離散方程

同理，對(duì)于隱格式求解方法，聯(lián)立方程（1）、（2）、（3）、（5）可得內(nèi)節(jié)點(diǎn)離散方程：

1.2 邊界條件處理

由于含有未知的邊界溫度，內(nèi)節(jié)點(diǎn)建立的代數(shù)方程組是不封閉的，需要給定邊界節(jié)點(diǎn)方程。常見(jiàn)的導(dǎo)熱問(wèn)題有3 類(lèi)邊界條件：規(guī)定邊界的溫度（第一類(lèi)邊界條件）、規(guī)定邊界的熱流密度（第二類(lèi)邊界條件）和規(guī)定表面?zhèn)鳠嵯禂?shù)h及周?chē)黧w的溫度tf（第三類(lèi)邊界條件）[7]。

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

圖2 邊界上的節(jié)點(diǎn)

當(dāng)邊界條件為規(guī)定溫度值時(shí)，可以對(duì)方程組直接求解并得到溫度場(chǎng)的分布。當(dāng)給定的邊界條件為第二類(lèi)或第三類(lèi)邊界條件時(shí)，需要分別進(jìn)行討論：

（1）第二類(lèi)邊界條件：規(guī)定傳入計(jì)算區(qū)域的熱量為正，此時(shí)將給定的熱流密度值q代入式（8）進(jìn)行求解；當(dāng)給定的是絕熱邊界時(shí)，可令式中

2 編程依據(jù)的算法流程

如圖 3 所示，基于 MATLAB 的傳熱學(xué)課程虛擬仿真實(shí)驗(yàn)平臺(tái)的軟件算法[8-9]求解基本流程如圖3 所示。

圖3 求解基本流程圖

（1）根據(jù)實(shí)際物理問(wèn)題建立數(shù)學(xué)控制方程，例如在一維非穩(wěn)態(tài)導(dǎo)熱問(wèn)題中，左端為恒定熱流，右端為與外界環(huán)境對(duì)流換熱求解，數(shù)值求解不同時(shí)間下的溫度分布。

（2）確定邊界條件和初始條件，在GUI 界面上輸入初始溫度、熱擴(kuò)散率、環(huán)境溫度、長(zhǎng)度和時(shí)間等參數(shù)數(shù)據(jù)[10-12]。

（3）劃分網(wǎng)格并對(duì)區(qū)域進(jìn)行離散化，建立代數(shù)方程。

（4）設(shè)立初值并進(jìn)行迭代求解，當(dāng)滿(mǎn)足收斂要求時(shí)將顯示計(jì)算結(jié)果，如圖4 所示。

圖4 仿真實(shí)驗(yàn)計(jì)算結(jié)果

3 功能模塊介紹

傳熱學(xué)課程的虛擬仿真實(shí)驗(yàn)平臺(tái)分為2 個(gè)模塊：一維區(qū)域傳熱模塊和二維區(qū)域傳熱模塊，每個(gè)模塊又根據(jù)導(dǎo)熱類(lèi)型分為穩(wěn)態(tài)傳熱和非穩(wěn)態(tài)傳熱，根據(jù)邊界條件分為3 類(lèi)熱邊界條件（恒定溫度、恒定熱流和對(duì)流換熱）。在各個(gè)模塊下設(shè)立彈出式選單，通過(guò)模塊的選項(xiàng)在面板中設(shè)置參數(shù)，從而進(jìn)入選定的模型中進(jìn)行虛擬仿真模型計(jì)算。

3.1 一維區(qū)域傳熱模塊

一維區(qū)域傳熱模塊包括一維穩(wěn)態(tài)傳熱模塊和一維非穩(wěn)態(tài)傳熱模塊，在選定區(qū)域和導(dǎo)熱類(lèi)型后，設(shè)置熱邊界條件、物性參數(shù)和空間區(qū)域(非穩(wěn)態(tài)需要?jiǎng)澐謺r(shí)間)，再通過(guò)點(diǎn)擊“計(jì)算”按鈕以對(duì)一維傳熱問(wèn)題進(jìn)行求解，計(jì)算的溫度分布圖顯示在坐標(biāo)軸上[13-14]。

3.1.1 數(shù)學(xué)模型

以一維非穩(wěn)態(tài)問(wèn)題為例，取邊界條件為恒定溫度的一維非穩(wěn)態(tài)導(dǎo)熱，相應(yīng)的數(shù)學(xué)表達(dá)式如式（9）所示：

在數(shù)值計(jì)算中需要采用的物理參數(shù)如表1 所示。

表1 一維非穩(wěn)態(tài)導(dǎo)熱問(wèn)題的物理參數(shù)表

對(duì)于本問(wèn)題的解析解為：

式中a為熱擴(kuò)散系數(shù)：

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

圖5 t=0.1 s 時(shí)數(shù)值解與分析解對(duì)比圖

3.1.2 GUI 運(yùn)行結(jié)果與分析

在數(shù)值計(jì)算過(guò)程中，分別設(shè)定時(shí)間為 0.005 s、0.025 s、0.05 s、0.1 s 和 100 s。在 GUI 的坐標(biāo)軸上溫度分布圖如圖6 所示。

圖6 一維區(qū)域非穩(wěn)態(tài)傳熱的溫度分布

選取每個(gè)設(shè)定時(shí)間最后的運(yùn)行結(jié)果進(jìn)行對(duì)比分析，其結(jié)果如圖7 所示。仿真結(jié)果表明：在一維非穩(wěn)態(tài)導(dǎo)熱的初期，受邊界低溫的影響，物體靠近邊界的溫度較低，中間段的溫度較高。隨著時(shí)間的推移，非穩(wěn)態(tài)導(dǎo)熱過(guò)程繼續(xù)進(jìn)行，物體的溫度持續(xù)下降，最終物體內(nèi)部各處的溫度無(wú)限接近于給定的邊界溫度。

圖7 不同時(shí)間的數(shù)值計(jì)算結(jié)果

通過(guò)本模塊的演示，學(xué)生可以學(xué)習(xí)到一維穩(wěn)態(tài)和非穩(wěn)態(tài)傳熱的基本規(guī)律，認(rèn)識(shí)到一維區(qū)域內(nèi)溫度的演化過(guò)程，掌握邊界條件和物性參數(shù)的因素對(duì)傳熱效果的影響。

3.2 二維區(qū)域傳熱模塊

二維區(qū)域傳熱模塊包括二維穩(wěn)態(tài)傳熱模塊和二維非穩(wěn)態(tài)傳熱模塊，相關(guān)的GUI 數(shù)值模擬計(jì)算步驟與一維區(qū)域傳熱模塊的計(jì)算步驟類(lèi)似：選定區(qū)域和導(dǎo)熱類(lèi)型，設(shè)置熱邊界條件、物性參數(shù)和空間區(qū)域（非穩(wěn)態(tài)需要?jiǎng)澐謺r(shí)間），再點(diǎn)擊“計(jì)算”按鈕以對(duì)二維傳熱問(wèn)題進(jìn)行求解[15-16]。

3.2.1 數(shù)學(xué)模型

以二維非穩(wěn)態(tài)傳熱問(wèn)題為例，二維區(qū)域傳熱模塊的功能的物理量參數(shù)如表2 所示。

表2 數(shù)值計(jì)算中的物理量參數(shù)表

3.2.2 GUI 運(yùn)行結(jié)果與分析

在 GUI 的時(shí)間劃分面板上，分別設(shè)定時(shí)間t為2.5×10-4s、5×10-4s、0.001 s、0.005 s、0.025 s 和 0.05 s進(jìn)行仿真。二維區(qū)域非穩(wěn)態(tài)傳熱的溫度分布模擬結(jié)果如圖8 所示。圖中矩形區(qū)域內(nèi)的溫度分布隨不同時(shí)間設(shè)定變化，說(shuō)明溫度的傳遞是從左端逐漸向內(nèi)部傳遞，并通過(guò)一定時(shí)間的熱傳導(dǎo)，溫度場(chǎng)的分布達(dá)到穩(wěn)定狀態(tài)。這是因?yàn)樵诰匦味S非穩(wěn)態(tài)導(dǎo)熱問(wèn)題中，右邊界設(shè)置為第三類(lèi)邊界條件（固體表面與周?chē)黧w間換熱），下表面設(shè)置為絕熱邊界，上邊界和左邊界為第二類(lèi)邊界條件（恒定熱流密度）且左邊界熱流遠(yuǎn)大于上邊界熱流，所以二維區(qū)域的熱源輸入主要由左端提供。在整個(gè)傳熱過(guò)程中，越靠近左邊界的溫度越高。

圖8 非穩(wěn)態(tài)傳熱條件下矩形區(qū)域的溫度分布

隨著時(shí)間的推移，在矩形區(qū)域內(nèi)通過(guò)熱傳導(dǎo)方式傳遞熱量，溫度由左邊界向內(nèi)部傳遞，最后達(dá)到穩(wěn)定的溫度分布。

通過(guò)本模塊的學(xué)習(xí)，學(xué)生可以仔細(xì)觀察二維條件下溫度場(chǎng)演化過(guò)程，掌握邊界條件對(duì)導(dǎo)熱的影響規(guī)律。

4 結(jié)語(yǔ)

基于 MATLAB 的傳熱學(xué)課程虛擬仿真實(shí)驗(yàn)平臺(tái)實(shí)現(xiàn)了一維和二維區(qū)域傳熱的穩(wěn)態(tài)/非穩(wěn)態(tài)傳熱數(shù)值模擬和結(jié)果可視化動(dòng)態(tài)演示，通過(guò)與解析解對(duì)比，驗(yàn)證了該虛擬仿真實(shí)驗(yàn)平臺(tái)運(yùn)算的可靠性。

虛擬仿真?zhèn)鳠釋W(xué)實(shí)驗(yàn)具有可視化效果好、操作簡(jiǎn)單、價(jià)格低廉、運(yùn)算準(zhǔn)確、高效的優(yōu)勢(shì)。將虛擬仿真實(shí)驗(yàn)軟件應(yīng)用于傳熱學(xué)的課程中，有助于學(xué)生對(duì)傳熱問(wèn)題的數(shù)值求解思路和方法的理解，同時(shí)也學(xué)習(xí)了MATLAB 在數(shù)值求解上的編程應(yīng)用和 GUI 的設(shè)計(jì)，具有較高的教學(xué)實(shí)用性。