康曉宇 吳惠 劉應開

(云南師范大學物理與電子信息學院 云南 昆明 650092)(昭通師專物理系 云南 昭通 657000)

1 引言

隨著計算機技術(shù)的發(fā)展和數(shù)值計算方法的改進，求解偏微分方程組等早期不能完成的任務現(xiàn)在逐步變?yōu)楝F(xiàn)實.在數(shù)值計算的各種軟件和軟件平臺中，MATLAB簡單易學且具有強大的可視化功能，為科學研究提供了極大的方便.特別是MATLAB中的PDE工具箱，無需編程則可直接求解具有特殊邊界條件的導熱問題，并將結(jié)果可視化，從而對物理過程和有關特征量的演化有一個直觀的認識.

在熱學中，熱傳導問題是一個基礎性的課題，對其進行研究有助于理解熱傳導的基本規(guī)律，澄清基本物理概念.在數(shù)學物理方法中，熱傳導問題是由熱傳導方程來描述.熱傳導方程的求解因邊界條件不同而難度各異.有的問題沒有嚴格的解析解，只能求得數(shù)值解；有的則存在嚴格的解析解，其解的形式可能是積分形式.也可能是級數(shù)形式，從這些解中很難準確知道溫度、熱流密度等物理量在空間的分布和隨時間的演化規(guī)律.本文利用MATLAB中的PDE工具箱求解兩類一維熱傳導問題，它們分別屬于第一類和第二類邊界條件.通過求解兩類問題中的溫度場、熱流密度分布與演化規(guī)律，揭示熱傳導的本質(zhì)，從而澄清穩(wěn)定態(tài)與平衡態(tài)的區(qū)別與聯(lián)系.

2 一維熱傳導方程

均勻金屬桿可看作一維傳熱系統(tǒng).如圖1所示，桿的左端作為坐標原點，則右端的坐標為L，桿上各點的溫度T(x,t)滿足偏微分方程

(1)

其中a=κ/cρ，κ，c和ρ分別為金屬桿的熱導率、比熱容和密度.

圖1 一維細桿熱傳導示意圖

3 兩類邊界條件的一維熱傳導問題

3.1 第一類邊界條件的一維熱傳導問題

長為L均勻鐵制金屬桿，桿的初始溫度為20℃，其左端與0℃的恒溫熱源接觸，右端與100℃恒溫熱源接觸，桿的側(cè)面與外界環(huán)境絕熱.現(xiàn)討論該金屬桿的溫度場分布、熱流密度分布及其演化情況，其方程如下

3.2 第二類邊界條件的一維熱傳導問題

長為L均勻鐵制金屬桿，初始時刻桿的左端溫度為0℃，右端溫度為100℃，桿上溫度梯度均勻分布，整個金屬桿與外界絕熱.現(xiàn)討論該金屬桿的溫度場分布、熱流密度分布及其演化情況，其方程如下

其中a=κ/cρ，計算中鐵制桿的密度，比熱容和熱導分別取ρ=7 800 kg/m3,c=500 J/(kg·℃),κ=48.6 W/(m·℃)[1].

3.3 求解的基本方法

問題的求解是基于有限元法[2]，其基本原理是把計算區(qū)域劃分成一系列的三角形單元，每個單元上取一個節(jié)點，選定一個形狀函數(shù)(拋物線型或雙曲線型)，并通過單元中節(jié)點上的被求解變量值表示該函數(shù).有限元法的最大優(yōu)點是對不規(guī)則區(qū)域的適應性好，故用MATLAB方法求解的結(jié)果在邊界上也比較精確[3].

MATLAB中的PDE工具箱在傳熱學中主要用于解如下兩類偏微分方程

式中，u為求解變量，t為時間變量，d，a，f，c為常數(shù)或變量.

MATLAB中的PDE工具箱定義了兩類邊界條件，即

式中n為垂直于邊界的單位矢量，h,r,q和t為常量或與u有關的變量[1，4].

利用PDE工具箱求解導熱問題的一般步驟如圖2所示.

圖2 MATLAB的計算流程圖

4 結(jié)果與討論

下面我們用MATLAB中的PDE工具箱來求它們的數(shù)值解.

4.1 第一類邊界條件的一維熱傳導問題中的溫度場、熱流密度分布與演化規(guī)律

將坐標軸的范圍設為0≤x≤5，-1≤y≤1，求解區(qū)域為0≤x≤4,-0.2≤y≤0.2.

矩形的上下邊界為Neuman邊界條件，取q=0，g=0， 左右邊界均為Dirichlet邊界條件， 分別取

h=1，r=0和h=1，r=100.

方程的類型為拋物線類型，其參數(shù)分別取c=48.6，a=0，d=7 800×500，f=0，初始條件為u(t0)=20[2].

圖3是經(jīng)過5小時，20小時，250小時后第一類邊界條件的金屬桿在空間的溫度分布云圖、熱流密度分布云圖，這些圖還顯示了有限元求解的網(wǎng)格.

圖3 經(jīng)過不同時間后第一類邊界條件的金屬桿溫度空間分布云圖及熱流密度分布云圖

由圖3(a1)可知，經(jīng)過5小時后，金屬桿兩端溫度梯度較大，中間溫度梯度較小.圖3(a2)為，隨著時間的增加，金屬桿不同部位的溫度梯度差異縮小，但溫度分布面仍為曲面；圖3(a3)為，當時間進一步延長，溫度分布面變得平直，說明溫度梯度趨于相等，金屬桿的溫度在空間的分布形成穩(wěn)定梯度，不同空間位置的溫度分布不再隨時間而發(fā)生明顯的變化.

由圖3(b1)可知，經(jīng)過5小時后，金屬桿不同部位熱流密度的分布是不同的，低溫端熱流密度比較平穩(wěn)，熱流密度小，而2 m～4 m區(qū)域，由于溫度從20℃～100℃變化，熱流密度變化較大，大量的熱由高溫熱源輸運到金屬桿中；圖3(b2)所示，金屬桿不同部位的熱流密度分布存在一定差異，即隨著時間的增加，空間各處熱流密度差異減??；圖3(b3)所示，表明經(jīng)過足夠長的時間金屬桿各部位的熱流密度趨于相等，對于某一給定區(qū)域來說，熱流流入多少，就幾乎流出多少，也就是說，經(jīng)過足夠長的時間后熱流穩(wěn)定地從高溫熱源輸運到低溫熱源，這與溫度在空間的穩(wěn)定分布一致.

4.2 第二類邊界條件的一維熱傳導問題的溫度場、熱流密度分布與演化規(guī)律

在使用PDE工具箱求解時，只需將上面求解過程中的左右邊界條件改為Neuman邊界條件，即取g=0，q=0和初始條件設為u(t0)=25×x即可，其余條件不變[2].

圖4 經(jīng)過不同時間后第二類邊界條件的金屬桿溫度空間分布云圖及熱流密度分布云圖

圖4是經(jīng)過5小時，20小時，250小時后第二類邊界條件的金屬桿在空間的溫度分布云圖、熱流密度分布云圖，這些圖還顯示了有限元求解的網(wǎng)格.

由圖4(a1)可知，經(jīng)過5小時后，金屬桿的溫度分布已發(fā)生變化，溫度梯度不再是常量，空間的溫度分布平面變彎曲；圖4(a2)為，隨著時間的增加，斜面變得較平坦，溫度梯度進一步減?。粓D4(a3)為，隨著時間的進一步增加，即經(jīng)過250小時后，溫度分布平面變得幾乎是平行xOy平面，溫度約為50℃，說明金屬桿上各點溫度近似相等，不存在明顯的溫度梯度.即經(jīng)過足夠長的時間后，金屬桿上的溫度為一個恒定的值50℃.

由圖4(b1)可知，經(jīng)過5小時后，金屬桿的熱流密度在桿的兩端近似為零，而桿的中間部位的熱流密度相當大并且近似為1 200 W/m2；由圖4(b2)可知，經(jīng)過20小時后，整個金屬桿中的熱流密度已經(jīng)大大減小，即隨著時間的增加，金屬桿各部位的熱流密度在減小；由圖4(b3)可知，經(jīng)過250小時后，金屬桿的熱流密度逐漸趨于零，也就是說，經(jīng)過足夠長的時間后金屬桿內(nèi)不存在能量的輸入與輸出，這與溫度在空間的均勻平衡分布是一致的.

5 結(jié)論

綜上所述，本文分別對第一類邊界條件下和第二類邊界條件下的金屬桿的一維熱傳導問題進行了可視化分析，從它們的溫度分布云圖和熱流密度云圖可以得出：第一類邊界條件下的金屬桿經(jīng)過足夠長的時間后，將處在一個宏觀性質(zhì)不隨時間變化的穩(wěn)定狀態(tài)，即空間各點的溫度是定值，但仍有穩(wěn)定的熱流從高溫熱源流向低溫熱源，能量在空間各點沒有堆積現(xiàn)象，是一個穩(wěn)定態(tài).而第二類邊界條件下的金屬桿經(jīng)過足夠長的時間后，金屬桿上各點的溫度均為50℃，不會隨時間的而變化，而且空間各點的熱流密度也趨于零，是一個平衡態(tài).從二者的比較可知：

(1)有外界影響的情況下，經(jīng)過足夠長的時間后，金屬桿內(nèi)空間各點的溫度有穩(wěn)定分布，但熱流不為零，屬于穩(wěn)定態(tài).

(2)無外界影響的情況下，經(jīng)過足夠長的時間后，金屬桿由非平衡態(tài)轉(zhuǎn)化為平衡態(tài)，空間各點溫度相等，也無熱流通過.

二者比較更能清楚地區(qū)分穩(wěn)定態(tài)與平衡態(tài)的區(qū)別與聯(lián)系，有助于對穩(wěn)定態(tài)和平衡態(tài)概念認識的深化.

參考文獻

1 彭芳麟.數(shù)學物理方程的MATLAB解法與可視化.北京:北京師范大學出版社,2004

2 陶文銓.數(shù)值傳熱學(第2版) .西安:西安交通大學出版社,2001

3 李燦,高彥棟,黃素逸.熱傳導問題的MATLAB計算.華中科技大學學報(自然科學版),2002,30(9):91-93

4 徐梓斌,閔劍青.基于PDE tool的熱傳導數(shù)值計算.佳木斯大學學報(自然科學版),2006,24(4):270-272