一類非線性熱傳導方程的數值解法

高忠社

（天水師范學院 數學與統計學院，甘肅 天水741001）

1 非線性熱傳導方程

熱量的傳播是發生在存在溫度差的同一物體或者物體之間，熱量從高溫物體傳到低溫物體的過程，就叫做熱傳導。而在實際應用中有時需要考慮隔熱問題，有時需要考慮散熱問題，從而盡量有效控制熱量的傳播。這種熱傳導現象利用數學方法描述，就得到了熱傳導方程。

在實際應用中，對于單個設備的熱傳導問題相對比較容易進行模擬與仿真，但是關于結構比較復雜的熱傳導問題、包含大量設備的復雜的多尺度系統的導熱問題等，其模擬仿真仍是十分困難的課題，比如復雜電路板的熱傳導問題。

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

該文將研究具有非線性導熱系數的一維熱傳導問題，假設熱源在最左端，熱量沿著細棒傳播，不考慮其他形式熱量傳播，如圖1所示。

圖1 熱傳導問題分析示意圖

在科學工程領域中，為了控制熱量的傳播而建立了一類具有非線性的熱傳導方程[1－3]

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

由方程（1）可得

如果k（u）和u無關，有，方程（1）是經典的線性拋物型偏微分方程，方程（2）為非線性的拋物型偏微分方程。

在半導體材料的導熱等很多熱傳導問題中，由于熱傳導系數隨著溫度的改變而改變，所以熱傳導系數通常是一個與溫度有關的函數。 如半導體材料的熱容量遠遠小于導熱系數，在實驗測定中發現，熱傳導系數近似于指數的變化，設

由于k0是非零的正常數，因此，只需考慮下列形式的一類非線性熱傳導方程

對方程（3）化簡，可得

設方程（1）在有限閉區間[a，b]上，滿足初邊值條件

其中邊值條件（5）給出了在區間兩端關于時間的溫度分布函數，初始條件（6）給出了在初始時刻的溫度分布情況。

2 時間方向的隱式Euler方法離散

對于?k（u）/?u＝0線性的熱傳導方程已經有很多數值方法[5-7]，如有限差分法、有限元法等。在文獻[3]的啟發下，將對方程（4）及初邊值條件（5）、（6），在時間方向使用隱式的Euler方法進行離散化。

設時間方向的步長為τ，時間t＞0，使用等距的劃分

對方程（4）利用隱式Euler方法進行離散化，則有

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

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

則方程（8）變形為

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

這樣式（11）、（12）、（13）組成了一個關于未知量為un的非線性兩點邊值問題。

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

設在區間[a，b]上等距劃分[8-10]，即

對式（10）中的二階導數，用二階中心差分算子

方程（10）可表示為

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

而

其中

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

4 差分格式的誤差分析

根據文獻[11-12]，對于上述的差分格式進行誤差估計，由式（8）、（15）、（18），有

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

對于上式中各項在（xi，tn）處Taylor級數展開

則有

將由式（21）、（22）、（23）代入式（20）整理，將式（20）與式（8）相減可得，截斷誤差

滿足

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

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

其中

則有

其中

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

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

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

記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關于un，i，vn，i的偏導數，即

其中

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

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

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

6 算例

考慮均質細桿的熱傳導問題，以細桿建立x軸，假設均質細桿介于x=0與x=2之間，熱源在坐標原點處，不考慮其他形式熱量傳播，設物體密度ρ和熱容比Cp都是常數，熱傳導系數和物體的溫度有關，即

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

初始溫度分布滿足

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

圖2 導熱參數k=-1熱傳導三維圖形和u對x平面圖形

圖3 導熱參數k=0熱傳導三維圖形和u對x平面圖形

圖4 導熱參數k=1熱傳導三維圖形和u對x平面圖形

7 結語

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