陳大偉, 斯小琴

(安徽建筑大學(xué)城市建設(shè)學(xué)院 a. 教育處; b. 基礎(chǔ)部, 安徽 合肥 238076)

如果物體內(nèi)各點(diǎn)的溫度不同,就有熱量從溫度較高地方流向溫度較低的地方.在不少生產(chǎn)實(shí)際問(wèn)題中,經(jīng)常需要考察物體上各點(diǎn)的溫度隨時(shí)間分布狀態(tài),這就要求解各種定解條件下的熱傳導(dǎo)方程.一維熱傳導(dǎo)方程是熱傳導(dǎo)方程的最簡(jiǎn)單形式,但是求解和對(duì)結(jié)果的理解說(shuō)明并不簡(jiǎn)單,這也使得此問(wèn)題成為該領(lǐng)域的一個(gè)熱點(diǎn)問(wèn)題.當(dāng)前有通過(guò)分離變量法與其他方法相結(jié)合獲得解析解,該類方法避免不了繁瑣的積分微分運(yùn)算,并且對(duì)復(fù)雜的結(jié)果很難做易理解的說(shuō)明[1-3];也有很多通過(guò)差分法和其他一些延伸方法可給出較高精度的數(shù)值解,但都涉及到高級(jí)計(jì)算機(jī)程序或較復(fù)雜的數(shù)學(xué)計(jì)算,這就使得在工程實(shí)際中直接運(yùn)用比較困難,并且對(duì)所獲得結(jié)果未能作出較全面和直觀的解釋[4-12].本文通過(guò)辦公軟件Microsoft Excel,可不用進(jìn)行語(yǔ)言編程即可快速獲得一維熱傳導(dǎo)方程高精度數(shù)值解,并通過(guò)OriginPro 8.0將巨大量的數(shù)值解數(shù)據(jù)模擬成溫度的時(shí)空分布曲面圖,方法簡(jiǎn)單實(shí)用,結(jié)果形象直觀[13].

1 模型與方法

考察一個(gè)簡(jiǎn)單模型,如均勻細(xì)桿內(nèi)熱量傳播過(guò)程.設(shè)細(xì)桿側(cè)面絕熱,熱量只能沿軸向傳導(dǎo).又設(shè)細(xì)桿橫截面面積為常數(shù)A且很小,以致在任何時(shí)刻,都可以把橫截面上的溫度視作相同,因此,是一個(gè)一維情形.圖1是密度為ρ、比熱為c、熱傳導(dǎo)系數(shù)為k的均勻細(xì)桿.

圖1 傳熱細(xì)桿模型Fig.1 The model of heat transfer rod

取x軸與細(xì)桿重合,以u(píng)(x,t)表示x點(diǎn)t時(shí)刻溫度,則有

ut=a2uxx+f(x,t).

(1)

現(xiàn)將溫度u(x,t)在(x0,t0)點(diǎn)沿x向前、向后h展開(kāi)為泰勒級(jí)數(shù),有

由式(2)和式(3),可得

(4)

將溫度u(x,t)在(x0,t0)點(diǎn)沿t向前τ展開(kāi)為泰勒級(jí)數(shù),有

將式(5)整理,得

(6)

聯(lián)立式(1)、式(4)和式(6),有

2 數(shù)值模擬結(jié)果及分析

考慮一實(shí)際問(wèn)題,一根長(zhǎng)為0.5 m的勻質(zhì)傳熱樞軸,熱擴(kuò)散系數(shù)a2=4.5×10-5,它的初溫為常數(shù)100 ℃,其兩端的溫度保持為0 ℃,考察該樞軸上溫度的分布情況.

首先,該問(wèn)題為一混合問(wèn)題,即為:

(9)

取τ=1 s、h=0.01 m,其熱傳導(dǎo)方程的差分方程為

0.5 m的空間長(zhǎng)度和3 000 s的時(shí)間長(zhǎng)度組成的時(shí)空平面被分割成15 000個(gè)時(shí)空節(jié)點(diǎn).在Excel中編寫(xiě)計(jì)算公式,并利用其下拉迭代計(jì)算功能迅速得到各節(jié)點(diǎn)的溫度值.從細(xì)桿端點(diǎn)0 m處開(kāi)始每隔0.1 m取一個(gè)特征位置點(diǎn),一直取到終點(diǎn),共取6個(gè)特征位置點(diǎn).將所有特征位置點(diǎn)各時(shí)刻的溫度值通過(guò)OriginPro 8.0模擬成溫度-時(shí)間曲線(如圖2所示);從0 s開(kāi)始每隔100 s取一個(gè)特征時(shí)間點(diǎn),一直取到400 s,共取5個(gè)特征時(shí)間點(diǎn).將各特征時(shí)間點(diǎn)各位置點(diǎn)的溫度值通過(guò)OriginPro 8.0模擬成溫度-位置曲線(如圖3所示).

圖2 各特征位置對(duì)應(yīng)的溫度-時(shí)間曲線

圖3 各特征時(shí)間對(duì)應(yīng)的溫度-位置曲線

圖2中0 m和0.5 m分別為端點(diǎn)位置.從圖中可以明顯看出,除兩個(gè)端點(diǎn)位置溫度恒定外,其他各考察點(diǎn)的溫度都隨著時(shí)間推移而降低,靠近端點(diǎn)處溫度降低更快,而從時(shí)間角度考察,起初一段較短時(shí)間溫度下降較“劇烈”,后期逐漸平緩;另外0和0.5、0.1和0.4、0.3和0.2 m位置處的時(shí)間-溫度曲線重合,反映出兩端點(diǎn)邊界條件相同、傳熱樞軸勻質(zhì),各點(diǎn)溫度變化也對(duì)稱的性質(zhì).而從圖3中可以發(fā)現(xiàn),0 s之后各時(shí)刻傳熱樞軸溫度呈現(xiàn)出兩端低中間高,而且隨著時(shí)間推移,各點(diǎn)(兩端點(diǎn)除外)溫度都在越來(lái)越“平緩”的降低;另外,各特征時(shí)刻位置-溫度曲線左右對(duì)稱形狀與圖2兩條對(duì)應(yīng)曲線重合相互驗(yàn)證.

圖2、圖3所反映出的傳熱樞軸各點(diǎn)溫度隨時(shí)間推移而降低且左右對(duì)稱,與熱量從高溫的中部向低溫的兩端對(duì)稱流失的實(shí)際完全相符.為將溫度隨時(shí)間、空間變化反映在一個(gè)圖上,通過(guò)OriginPro 8.0將數(shù)據(jù)轉(zhuǎn)化成矩陣,并通過(guò)該軟件的Plot 3D Surface功能模擬三維圖像.

由于溫度在前期變化較快,為突出溫度變化的影響因素,取前400 s所有數(shù)據(jù)繪制三維圖像,如圖4所示.為考察熱傳導(dǎo)后期特點(diǎn)每隔200 s取數(shù)據(jù)直至3 000 s,并將其模擬三維圖像,如圖5所示.

圖4 前400 s溫度時(shí)空分布曲面

圖5 前3 000 s溫度時(shí)空分布曲面

在圖4中,沿時(shí)間軸向前看,傳熱樞軸溫度整體降低,且表現(xiàn)“前急后緩”的整體特征;沿位置軸向前看,傳熱樞軸溫度兩端低中間高,且表現(xiàn)“左右對(duì)稱”的整體特征.這些特征與圖2、圖3所反映出的特征是一致的,但圖4對(duì)特征信息的反映更整體、全面、集成,且通過(guò)不同灰度對(duì)比,更直觀形象.

圖5是圖4在時(shí)間上的繼續(xù)延伸,從圖中可以發(fā)現(xiàn)溫度“越來(lái)越平緩”地下降,以至于到1 800 s后整個(gè)樞軸上各點(diǎn)溫度相差無(wú)幾,各位置點(diǎn)溫度隨時(shí)間變化也幾乎為零,溫度的時(shí)空分布近似為平面,此時(shí)接近于“穩(wěn)態(tài)”,這一特征與工程實(shí)際相吻合.

3 結(jié) 論

(1) 只要滿足穩(wěn)定條件,且時(shí)間、空間步距足夠小,利用熱傳導(dǎo)方程的差分格式就可以獲得一定精度要求的熱傳導(dǎo)方程數(shù)值解.

(2) 步距取值小,將時(shí)空平面分成的網(wǎng)格節(jié)點(diǎn)數(shù)量大,計(jì)算量就很大(本文計(jì)算15 000個(gè)時(shí)空節(jié)點(diǎn)).但是利用基本辦公軟件EXCEL的下拉迭代計(jì)算功能可以快速、簡(jiǎn)便地完成大量節(jié)點(diǎn)的溫度計(jì)算.

(3) 利用Origin將產(chǎn)生的大量數(shù)據(jù)繪制成三維曲面圖,能直觀地得到熱傳導(dǎo)媒質(zhì)上的溫度時(shí)空分布,該方法所得到的結(jié)論與工程實(shí)際相符.

(4) 對(duì)于傳熱媒質(zhì)中有熱源即熱傳導(dǎo)方程中有非齊次項(xiàng)f(x,t)時(shí),此方法同樣適用.

(5) 本文利用最簡(jiǎn)單辦公軟件獲得一維熱傳導(dǎo)數(shù)值解,并通過(guò)Origin將數(shù)值解數(shù)值模擬成相應(yīng) 圖像.該方法獲得與其他方法同樣的結(jié)果[15-16],但過(guò)程簡(jiǎn)單、快捷、方便,在解決此類工程實(shí)際問(wèn)題中有一定的參考和使用價(jià)值.