康曉宇 吳惠 劉應(yīng)開(kāi)

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

1 引言

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

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

2 一維熱傳導(dǎo)方程

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

(1)

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

圖1 一維細(xì)桿熱傳導(dǎo)示意圖

3 兩類邊界條件的一維熱傳導(dǎo)問(wèn)題

3.1 第一類邊界條件的一維熱傳導(dǎo)問(wèn)題

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

3.2 第二類邊界條件的一維熱傳導(dǎo)問(wèn)題

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

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

3.3 求解的基本方法

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

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

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

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

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

利用PDE工具箱求解導(dǎo)熱問(wèn)題的一般步驟如圖2所示.

圖2 MATLAB的計(jì)算流程圖

4 結(jié)果與討論

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

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

將坐標(biāo)軸的范圍設(shè)為0≤x≤5，-1≤y≤1，求解區(qū)域?yàn)?≤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)過(guò)5小時(shí)，20小時(shí)，250小時(shí)后第一類邊界條件的金屬桿在空間的溫度分布云圖、熱流密度分布云圖，這些圖還顯示了有限元求解的網(wǎng)格.

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

參考文獻(xiàn)

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

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

3 李燦,高彥棟,黃素逸.熱傳導(dǎo)問(wèn)題的MATLAB計(jì)算.華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版),2002,30(9):91-93

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