熱傳導方程的積分變換求解與MATLAB實現(xiàn)

丁 波 歐志華 奉瑞萍

熱傳導方程的積分變換求解與MATLAB實現(xiàn)

丁 波 歐志華 奉瑞萍

（湖南工業(yè)大學 土木工程學院，湖南 株洲 412007）

熱傳導方程是物理學中經(jīng)典方程之一，反映熱的作用規(guī)律。有關熱傳導方程的解法有分離變量法、延拓法、特殊函數(shù)法、積分變換法等。本文首先介紹Fourier變換和Laplace變換在求解熱傳導方程中的應用，然后以拋物型方程Heat Equation為例，給出了通過MATLAB中的偏微分方程工具箱PDE Toolbox進行建模求解及模擬直觀。通過空間中的溫度場分布和梯度方向，發(fā)現(xiàn)溫度是由溫度高的一面向低的一面變化，而且最靠近熱源的分子首先被加熱，經(jīng)過一段時間，溫度場中每一點都非常接近與它所能到達的最高溫度，逐漸趨于穩(wěn)態(tài)。

熱傳導方程；Fourier變換；Laplace變換；MATLAB實現(xiàn)

熱的作用服從于一些不變的規(guī)律，如果不借助數(shù)學工具分析就不能精確地發(fā)現(xiàn)和了解這些規(guī)律。法國數(shù)學家和物理學家傅里葉(Jean-Baptiste-Joseph Fourier，1768-1830)，最早建立熱傳導方程來研究熱的作用規(guī)律，1807年，他向科學院提交的論文《關于熱傳導的研究報告》中建立了關于不連續(xù)的物質(zhì)和特殊形狀的連續(xù)體的熱擴散（即熱傳導）方程。和重力、磁場一樣，熱貫穿于宇間的一切物質(zhì)之中，熱傳導方程是物理學經(jīng)典方程之一，反映了溫度在空間中如何變化[1]。熱傳導方程是二階線性偏微分方程，解的呈現(xiàn)形式復雜且多樣，孫桂榮[2]等運用Nevanlinna值分布理論討論了二階線性微分方程解的徑向振蕩問題，楊琰琰[3]等運用Nevanlinna值分布理論和方法考察了一類針對非線性差分方程解，并得到解的增長性和表示，對于求解熱傳導方程的數(shù)學方法有很多，其中針對齊次方程的有分離變量法、無窮級數(shù)法等[4]，針對非齊次方程的有函數(shù)變易法、齊次化原理[5-6]、特殊函數(shù)法[7]等等。積分變換法[8-10]通過將偏微分方程轉(zhuǎn)化為常微分方程或代數(shù)方程，而且不論是求解齊次方程還是非齊次方程，以及對于無界問題和半無界問題的求解，都避開了復雜的函數(shù)處理技巧與求解運算，十分方便得到方程的解。本文主要介紹Fourier變換和Laplace變換在求解熱傳導方程中的應用，另外利用編程和可視化極強的軟件——MATLAB求解偏微分方程[11]，可以直觀地反映空間中的溫度場是如何變化的，并可以得到溫度梯度的方向，這對于指導實際工程有很大的指導意義。

1　熱傳導方程及其定解問題的導出

1.1 熱傳導方程

其中(,,)稱為物體在點(,,)處的熱傳導系數(shù)(＞0)，等式右邊的負號表示熱流的方向與溫度梯度相反，由于熱量總是從溫度高的一側(cè)流向低的一側(cè)。

在物體G內(nèi)任取一封閉曲面Γ，它所包圍的區(qū)域記為Ω，從時刻1到2流進此閉曲面的全部熱量為

流入的熱量使物體內(nèi)部溫度發(fā)生變化，則在同一時間段(1,2)內(nèi)，它所吸收的熱量為

其中為比熱，為密度。若物體內(nèi)沒有熱源，則根據(jù)熱量守恒定律，得1=2，即

假設函數(shù)關于空間變量,,具有二階連續(xù)偏導數(shù)，關于時間具有一階連續(xù)偏導數(shù)，利用高斯公式[12]將第二型曲面積分轉(zhuǎn)化為三重積分，上式整理為

此為非均勻各向同性體的熱傳導方程。若物體均勻，、及均為常數(shù)，記/=2，有

上式稱為三維的齊次熱傳導方程，若物體內(nèi)部有熱源，設熱源密度即單位時間內(nèi)、單位體積內(nèi)發(fā)出的熱量為(,,)，記(,,,)=(,,)/，得非齊次方程[9]

1.2 定解問題

Fourier所做的工作就是將諸多不同情況的熱擴散（或熱傳導）問題，轉(zhuǎn)化為在不同邊界條件和不同初始條件下的求解。從物理學的角度來看，如果知道了物體在邊界上的溫度分布（或熱交換狀況）和物體在初始時刻的溫度，就可以完全確定物體在接下來任意時刻的溫度分布，因此熱傳導方程的定解問題就是在已給定的初始條件與邊界條件下求問題的解。對于熱傳導方程的初始條件的提法顯然為

其中(,,)為已知函數(shù)，表示物體在=0時的溫度分布。三類邊界條件如下：

1）物體的表面的溫度是已知的，為

其中Γ表示物體的邊界曲面，(,,,)是定義在(,,)∈Γ，0≤≤T上的已知函數(shù)，此稱為熱傳導方程的第一類邊界條件(又稱狄利克雷(Dirichlet)邊界條件)；

2）溫度函數(shù)在表面上的外法向方向?qū)?shù)是已知的，為

此稱為熱傳導方程的第二類邊界條件(又稱諾伊曼(Neumann)邊界條件)；

3）溫度函數(shù)及其外法向方向?qū)?shù)的現(xiàn)行組合在界面上是已知的，為

其中為已知正數(shù)，此稱為熱傳導方程的第三類邊界條件。

上述的三類邊界條件加上方程構(gòu)成的定解問題分別稱為：第一邊值問題(Dirichlet問題)，第二類邊值問題(Neumann問題)和第三類邊值問題(Robin問題)。

2 Fourier變換求解熱傳導方程

2.1 Fourier變換定義及其性質(zhì)

定義：設()是定義在(-∞,+∞)上的絕對可積函數(shù)，且在任意有限區(qū)間上滿足Dirichlet收斂條件，則稱函數(shù)

為()的Fourier變換，記為[]；稱函數(shù)

相比較于分離變量法中用Fourier級數(shù)來求解熱傳導方程，F(xiàn)ourier變換可以看作是Fourier級數(shù)的極限形式，是將一種代表某種關于時間或空間的信號（signal）從時域轉(zhuǎn)換到頻域內(nèi)表達的手段。

性質(zhì)1：線性性質(zhì)：對于任意復數(shù)α，β以及函數(shù)1，2,成立

2.2 一維熱傳導方程的初值問題(無界問題)

定解問題(1)中的方程和初始條件兩邊關于作Fourier變換，記

利用2.1節(jié)Fourier變換的性質(zhì)，(2)代入(1)中得

方程(3)是一階線性非齊次常微分方程，可由常數(shù)變易法得

然后對方程(4)兩邊求反演，記

則

因為

很容易得

同理可得

利用2.1節(jié)Fourier變換性質(zhì)2和疊加原理得

此即為用Fourier變換法求得一維熱傳導方程問題的解。

同理，也可得到二維熱傳導方程柯西問題

解的表達式為[13]

3 Laplace變換求解熱傳導方程

3.1 Laplace變換定義及其性質(zhì)

因為Fourier變換要求函數(shù)在(-∞,+∞)上有定義，如果我們研究的半無界問題（混合問題），則不僅不能對進行Fourier變換，而且也不能對進行變換(≥0)，為此，對于定義在≥0中的函數(shù)()，考慮Laplace變換，通過Laplace變換將微分方程轉(zhuǎn)化為常微分方程或代數(shù)方程[14]。

定義：對于函數(shù)()，函數(shù)

為()的Laplace變換，記為[]；稱函數(shù)

為()到()的Laplace逆變換，記為-1[].

性質(zhì)1：線性性質(zhì)：設?1、2，由Laplace變換的定義，得

性質(zhì)2：微分性質(zhì)：設函數(shù)()的導函數(shù)的Laplace變換存在，則有

一般地，有

性質(zhì)3：卷積性質(zhì)：

3.2 Laplace變換求解熱傳導方程的混合問題

方程(10)與邊值條件(11)兩邊關于作Laplace變換，記

利用3.1節(jié)Laplace變換的性質(zhì)1和2，得

可見，通過Laplace變換使原本偏微分方程(9)轉(zhuǎn)化成常微分方程(11)，易得其特解為

記

因為

代入(12)式中，并由卷積性質(zhì)3，得

其中

由(,)的表達式(13)及(15)，得

因此，利用微分性質(zhì)2，得[15]

4 MATLAB中求解熱傳導方程

根據(jù)積分變換法得出熱傳導方程的解含有空間坐標函數(shù)和時間函數(shù)，從解的形式上來看，很難直觀看出熱的傳導規(guī)律，溫度場的變化和溫度梯度的方向。MATLAB是具有極高編程效率和強大作圖功能的科技應用軟件之一，通過MATLAB中的偏微分方程工具箱(PDE Toolbox)對拋物型方程(Heat Equation為例)進行求解。為研究熱在介質(zhì)中的擴散規(guī)律，現(xiàn)考慮一個帶有矩形孔(孔長0.8，寬0.1)的金屬板(板長1.6，寬1)的熱傳導問題，板的左邊界保持在100℃，板的右邊界熱量從板向環(huán)境空氣定常流動，其余各邊界保持絕緣。初始條件=0時板的溫度為0℃，由上述條件可得到該定解問題為：

因為在MATLAB中對于拋物型方程的一般形式為：

其中=(,)，、、和是標量復函數(shù)形式的系數(shù)，所以在PDE Specification中設置方程為拋物型（Parabolic），=1，=1，=0，=0；左邊界為狄利克雷（Dirichlet）條件：=1，=100；右邊界為諾伊曼（Neumann）條件：=－1，=0；其余邊界為諾伊曼（Neumann）條件：=0，=0。

方程建立、設置邊界條件和網(wǎng)格劃分，經(jīng)過計算，分別給出Time=0，1，3，5，10，100時刻的溫度場和梯度圖：

圖1 建模、邊界條件建立

圖2 網(wǎng)格劃分

圖3 Time=0時刻的溫度場和梯度方向

圖4 Time=1時刻的溫度場和梯度方向

圖5 Time=3時刻的溫度場和梯度方向

圖6 Time=5時刻的溫度場和梯度方向

圖7 Time=10時刻的溫度場和梯度方向

圖8 Time=100時刻的溫度場和梯度方向

可以發(fā)現(xiàn)板內(nèi)溫度是由溫度高的一面向低的一面變化，而且最靠近熱源的分子將首先被加熱，經(jīng)過一段時間后，金屬板的每一點都將獲得非常接近于它所能達到的最高溫度。在Time=10之后至Time=100板內(nèi)溫度變化不大，也就是說隨著時間越來越長，板內(nèi)的溫度也逐漸趨于穩(wěn)態(tài)。每個點的最高溫度取決于距離熱源的距離，且溫度隨離熱源越遠而越低。

5 結(jié)束語

熱傳導方程是反映了溫度變化的規(guī)律，通過積分變換法應用于熱傳導方程的求解十分簡便，積分變換的核心是求象函數(shù)的逆變換，通過卷積性質(zhì)可以方便的求出函數(shù)(,)，且積分變換法對于其他偏微分方程例如波動方程，也是十分適用。通過MATLAB軟件的模擬直觀地反映了溫度在介質(zhì)中的變化規(guī)律，對解決實際生活中的工程問題具有指導意義。

[1]傅立葉.熱的解析理論[M].桂質(zhì)亮,譯.北京:北京大學出版社,2008.

[2]孫桂榮,楊琰琰.二階復域微分方程解的徑向振蕩[J].蘇州科技大學學報(自然科學版),2018,35(02):11-14.

[3]楊琰琰,黃志剛,胡夢薇.一類非線性差分方程解的若干性質(zhì)[J].蘇州科技大學學報(自然科學版),2018,35(03):23-26.

[4]王暢,王翹楚,劉偉,史慶藩.數(shù)學物理方程求解中的創(chuàng)新思維探源[J].大學物理,2018,37(09):56-59.

[5]劉子建.一類二階常系數(shù)非齊次線性微分方程初值問題的若干種解法[J].黑龍江科技信息,2016(34):54.

[6]樊龍,李高.利用齊次化原理求解常系數(shù)非齊次線性方程初值問題[J].大學數(shù)學,2017,33(02):111-113.

[7]朱一超.數(shù)學物理方程與進階分析工具[M].北京:科學出版社,2020.

[8]代瑩,肖冰.熱傳導方程中傅氏積分與傅氏變換的應用[J].新疆師范大學學報(自然科學版),2019,38(02):5-15.

[9]數(shù)學物理方程/谷超豪等編著.—3版.—北京：高等教育出版社,2012.7

[10]張景軍,郭艷鳳.拉普拉斯變換法求解超短脈沖方程[J].數(shù)學的實踐與認識,2018,48(21):220-225.

[11]趙迎春,布仁滿都拉.偏微分方程的積分變換法及其MATLAB解算[J].現(xiàn)代計算機(專業(yè) 版),2016(29):53-55.

[12]周軒,崔虹云.電磁場中的高斯定理證明及其巧妙解題運用[J].廣東化工,2019,46(21):168-169.

[13]何維明.熱傳導方程Cauchy問題某些特殊情形的簡便算法[J].長沙交通學院學報,1994(04):14-19.

[14]郭維斌.拉普拉斯(Laplace)變換法解常微分方程的初值問題[J].數(shù)學學習與研究,2014(15):124.

[15]崔海波.Laplace變換在偏微分方程中的應用[J].教育教學論壇,2017(04):219-220.

O175.24

A

1673-2219（2021）03-0003-05

2021-03-11

湖南省建設廳科技計劃基金資助項目（項目編號C10108）；湖南工業(yè)大學研究生科研創(chuàng)新項目（項目編號CX1924）。

丁波（1996－），男，江蘇鹽城人，碩士研究生，研究方向為建筑材料。

歐志華（1975－），男，湖南常寧人，博士，高級工程師，研究方向為高性能混凝土水泥基材料。

（責任編校：文春生）