一類(lèi)非線性熱傳導(dǎo)方程的數(shù)值解法

高忠社

（天水師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，甘肅 天水741001）

1 非線性熱傳導(dǎo)方程

熱量的傳播是發(fā)生在存在溫度差的同一物體或者物體之間，熱量從高溫物體傳到低溫物體的過(guò)程，就叫做熱傳導(dǎo)。而在實(shí)際應(yīng)用中有時(shí)需要考慮隔熱問(wèn)題，有時(shí)需要考慮散熱問(wèn)題，從而盡量有效控制熱量的傳播。這種熱傳導(dǎo)現(xiàn)象利用數(shù)學(xué)方法描述，就得到了熱傳導(dǎo)方程。

在實(shí)際應(yīng)用中，對(duì)于單個(gè)設(shè)備的熱傳導(dǎo)問(wèn)題相對(duì)比較容易進(jìn)行模擬與仿真，但是關(guān)于結(jié)構(gòu)比較復(fù)雜的熱傳導(dǎo)問(wèn)題、包含大量設(shè)備的復(fù)雜的多尺度系統(tǒng)的導(dǎo)熱問(wèn)題等，其模擬仿真仍是十分困難的課題，比如復(fù)雜電路板的熱傳導(dǎo)問(wèn)題。

文獻(xiàn)[1-4]中H.S.Carslaw，J.Lienemann等人給出了熱傳導(dǎo)問(wèn)題的非線性模型，S.M.Filipov等人研究了一般形式的非線性熱傳導(dǎo)方程的計(jì)算問(wèn)題。 由于不同材料導(dǎo)熱系數(shù)隨著溫度的變化而不同，筆者將研究一類(lèi)具有指數(shù)型熱傳導(dǎo)系數(shù)的熱傳導(dǎo)方程的計(jì)算問(wèn)題，對(duì)于時(shí)間方向使用隱式Euler差分格式，空間方向使用二階中心差分格式。 利用Taylor級(jí)數(shù)方法得到該差分格式的誤差階為O（τ+h2），對(duì)于離散化后的代數(shù)方程組利用Newton迭代法進(jìn)行求解，最后通過(guò)數(shù)值算例分析討論了該方法的可行性。

該文將研究具有非線性導(dǎo)熱系數(shù)的一維熱傳導(dǎo)問(wèn)題，假設(shè)熱源在最左端，熱量沿著細(xì)棒傳播，不考慮其他形式熱量傳播，如圖1所示。

圖1 熱傳導(dǎo)問(wèn)題分析示意圖

在科學(xué)工程領(lǐng)域中，為了控制熱量的傳播而建立了一類(lèi)具有非線性的熱傳導(dǎo)方程[1－3]

其中u（x，t）表示在（x，t）處的溫度，ρ為材料密度，Cp是材料的熱容量，k（u）（u表示溫度）為熱傳導(dǎo)系數(shù)，即它與溫度有關(guān)，所得方程是典型的非線性熱傳導(dǎo)方程。

由方程（1）可得

如果k（u）和u無(wú)關(guān)，有，方程（1）是經(jīng)典的線性?huà)佄镄推⒎址匠?，方程?）為非線性的拋物型偏微分方程。

在半導(dǎo)體材料的導(dǎo)熱等很多熱傳導(dǎo)問(wèn)題中，由于熱傳導(dǎo)系數(shù)隨著溫度的改變而改變，所以熱傳導(dǎo)系數(shù)通常是一個(gè)與溫度有關(guān)的函數(shù)。 如半導(dǎo)體材料的熱容量遠(yuǎn)遠(yuǎn)小于導(dǎo)熱系數(shù)，在實(shí)驗(yàn)測(cè)定中發(fā)現(xiàn)，熱傳導(dǎo)系數(shù)近似于指數(shù)的變化，設(shè)

由于k0是非零的正常數(shù)，因此，只需考慮下列形式的一類(lèi)非線性熱傳導(dǎo)方程

對(duì)方程（3）化簡(jiǎn)，可得

設(shè)方程（1）在有限閉區(qū)間[a，b]上，滿(mǎn)足初邊值條件

其中邊值條件（5）給出了在區(qū)間兩端關(guān)于時(shí)間的溫度分布函數(shù)，初始條件（6）給出了在初始時(shí)刻的溫度分布情況。

2 時(shí)間方向的隱式Euler方法離散

對(duì)于?k（u）/?u＝0線性的熱傳導(dǎo)方程已經(jīng)有很多數(shù)值方法[5-7]，如有限差分法、有限元法等。在文獻(xiàn)[3]的啟發(fā)下，將對(duì)方程（4）及初邊值條件（5）、（6），在時(shí)間方向使用隱式的Euler方法進(jìn)行離散化。

設(shè)時(shí)間方向的步長(zhǎng)為τ，時(shí)間t＞0，使用等距的劃分

對(duì)方程（4）利用隱式Euler方法進(jìn)行離散化，則有

其中un=un（x），un-1＝un-1（x）分別表示u（x，tn）和u（x，tn-1）的近似值。

為了求解方程（8），設(shè)vn=dun/dx

則方程（8）變形為

根據(jù)式（5），邊值條件可以離散為

這樣式（11）、（12）、（13）組成了一個(gè)關(guān)于未知量為un的非線性?xún)牲c(diǎn)邊值問(wèn)題。

3 空間方向的二階離散方法

設(shè)在區(qū)間[a，b]上等距劃分[8-10]，即

對(duì)式（10）中的二階導(dǎo)數(shù)，用二階中心差分算子

方程（10）可表示為

其中un，i+1，un，i，un，i－1分別是un（xi+1），un（xi），un（xi-1）的近似值。即溫度函數(shù)u（x，t）在點(diǎn)（xi，tn）的近似值。

而

其中

根據(jù)式（13）及邊值條件（16），得到下列方程組

4 差分格式的誤差分析

根據(jù)文獻(xiàn)[11-12]，對(duì)于上述的差分格式進(jìn)行誤差估計(jì)，由式（8）、（15）、（18），有

其中un，i表示u（xi，tn）的近似值。

對(duì)于上式中各項(xiàng)在（xi，tn）處Taylor級(jí)數(shù)展開(kāi)

則有

將由式（21）、（22）、（23）代入式（20）整理，將式（20）與式（8）相減可得，截?cái)嗾`差

滿(mǎn)足

5 非線性方程組的牛頓迭代法求解

對(duì)式（19）用非線性牛頓迭代法進(jìn)行求解，建立N行的列向量

其中

則有

其中

給定初值un（0），對(duì)非線性方程組（26）利用牛頓迭代法求解

其中Jn（k）是Gn的Jacobian矩陣，即

方程組（25）中第二式中f（un，i，vn，i；un，i－1），i＝2，3，…，N－1，與un，i，vn，i有關(guān)，為了分析Jacobian矩陣元素，先求f（un，i，vn，i；un，i－1），i＝2，3，…，N－1關(guān)于un，i，vn，i的偏導(dǎo)數(shù)。

記fn＝f（un，i，vn，i；un，i－1），gn＝g（un，i，vn，i；un，i－1），而qn＝q（un，i，vn，i；un，i－1）和pn＝p（un，i，vn，i；un，i－1）分別表示fn關(guān)于un，i，vn，i的偏導(dǎo)數(shù)，即

其中

將式（32）、（33）代入式（30）、（31），則有

綜合分析，可得Jacobian矩陣的元素

式（28）是一步兩層迭代格式，給定初值un（0），可得近似迭代解{un（k+1）}。如果此序列收斂，可得該向量序列的極限是非線性代數(shù)方程組的解。在實(shí)際的計(jì)算中通常使用||un（k+1）-un（k）||＜ε結(jié)束迭代過(guò)程。

6 算例

考慮均質(zhì)細(xì)桿的熱傳導(dǎo)問(wèn)題，以細(xì)桿建立x軸，假設(shè)均質(zhì)細(xì)桿介于x=0與x=2之間，熱源在坐標(biāo)原點(diǎn)處，不考慮其他形式熱量傳播，設(shè)物體密度ρ和熱容比Cp都是常數(shù)，熱傳導(dǎo)系數(shù)和物體的溫度有關(guān)，即

假定物體密度ρ=1、熱容比Cp=1和k0=0.1，物體在兩端的溫度為常數(shù)，即

初始溫度分布滿(mǎn)足

則式（4）、（38）、（39）構(gòu)成非線性熱傳導(dǎo)方程的初邊值問(wèn)題。根據(jù)文中上述方法，可得該問(wèn)題的數(shù)值解，其中給定空間方向x∈[0，2]，步長(zhǎng)h＝0.05，時(shí)間步長(zhǎng)τ=0.5，時(shí)間t∈[0，15]，指數(shù)參數(shù)k=-1，0，1時(shí)給出了三維圖形和u關(guān)于x的圖形，如圖2、圖3、圖4所示。

圖2 導(dǎo)熱參數(shù)k=-1熱傳導(dǎo)三維圖形和u對(duì)x平面圖形

圖3 導(dǎo)熱參數(shù)k=0熱傳導(dǎo)三維圖形和u對(duì)x平面圖形

圖4 導(dǎo)熱參數(shù)k=1熱傳導(dǎo)三維圖形和u對(duì)x平面圖形

7 結(jié)語(yǔ)

熱傳導(dǎo)方程是一個(gè)經(jīng)典的拋物方程，同時(shí)該方程也廣泛應(yīng)用于很多科學(xué)工程領(lǐng)域。 文中根據(jù)文獻(xiàn)[3]討論的一般類(lèi)型非線性熱傳導(dǎo)問(wèn)題，討論了一類(lèi)具有指數(shù)型熱傳導(dǎo)系數(shù)的熱傳導(dǎo)方程的計(jì)算問(wèn)題。由于很多材料的導(dǎo)熱系數(shù)隨著溫度的不同而不同，文中主要研究了熱傳導(dǎo)系數(shù)為k（u）=k0eku的非線性熱傳導(dǎo)方程的數(shù)值解法。 其中時(shí)間方向利用隱式Euler差分格式、空間方向利用二階中心差分格式進(jìn)行離散化，同時(shí)得到該差分格式的誤差階為O（τ+h2），離散化后的代數(shù)方程組用Newton迭代法進(jìn)行求解。 最后通過(guò)數(shù)值算例，說(shuō)明該方法具有一定的可行性和實(shí)用性。