基于邊界元法的瞬態熱傳導仿真分析

張馳，張碩

（沈陽航空航天大學民用航空學院，遼寧沈陽110136）

基于邊界元法的瞬態熱傳導仿真分析

張馳，張碩

（沈陽航空航天大學民用航空學院，遼寧沈陽110136）

首先推導了瞬態熱傳導的邊界積分方程，然后通過一系列變換得到了易求解的矩陣形式，提出用迭代法求解瞬態熱傳導問題.最后引入數值算例，計算了溫度分布及熱流密度分布，并與解析解進行比較.結果表明采用邊界元法所得的數值仿真解與解析解吻合，證明此方法的有效性.

瞬態熱傳導；邊界元法；熱流密度；溫度場

邊界元法相對有限元法能夠大大減少計算時間和存儲容量，而且邊界元不需要對區域內進行離散，因此溫度的測量點位置可以任意選擇，尤其是這種方法可以直接給出未知表面的溫度和測量困難的熱流[1-2]，因此邊界元法在傳熱學領域的工程實際應用具有重大的意義.

1 瞬態熱傳導的基本方程

瞬態熱傳導問題與穩態相比要涉及與時間相關的問題，因此相對復雜.對于無內熱源點的瞬態熱傳導問題，有

方程的基本解[3]為

其中d是空間維數.基本解（2）代入到基本方程（1）的積分方程中可得

這就是瞬態熱傳導問題的邊界積分方程.最后一項對應著t=0的初始條件.

2 邊界元法

對時間域進行劃分，假設函數T及q隨時間變化，由于二者比T*和q*變化慢得多，故近似在小的時間間隔內為常數，所以可以對（3）式分段積分，并對時間內層求積分得

其中式中D為點源i到邊界單元線的垂直距離，Ei項是指數積分函數，可以通過級數來計算[4].

通過上述計算后（3）式可以寫成

對空間域劃分，邊界Γ劃分成N個單元，域Ω劃分成M個單元，寫成矩陣形式即

式中G、H是系數矩陣，只依賴于幾何形狀、材料特性及時間步長；Pit1表示整個空間域Ω上的初始溫度分布對i點的貢獻.

在求解邊界積分方程時采用迭代解法，把時間域分成N等分，時間步長Δt=tmax/N.這樣再把空間域進行劃分，根據初始條件就可以計算出邊界上的溫度場分布和熱流密度分布.

3 數值算例

根據以上推導結果，開發C語言程序包，并采用VC++進行嚴格調試，通過數值算例對上述方法的正確性加以驗證.

方板邊長L=1，板沿x方向和y方向的熱傳導率均為k=1，比熱c=1，如圖1所示.

圖1 二維方板的邊界條件

其他邊界條件為

仿真結果如下所示.

圖2 t=0.75時刻y=0邊的溫度變化

圖3 y=1邊不同時刻沿x向的溫度變化

以上圖中邊界元解與有限元解及解析解均能較好的吻合，證明該方法的正確性.

4 結論

從上述的仿真分析中，本文所采用的邊界元法結果準確.該方法理論推導簡單，所得結果誤差小且程序編寫相對容易.

這種方法為材料的設計及應用其他相關問題如參數優化、靈敏度分析等提供算法依據，在多維坐標及彈性力學、斷裂力學、梯度材料等方面的應用也值得探討.

〔1〕吳洪潭.邊界元法在傳熱學中的應用[M].北京：國防工業出版社,2008.

〔2〕AlokSutradhar,G.H.Paulino.Thesimpleboundaryelementmethodfortransientheatconductioninfunctionally gradedmaterials[J].ComputerMethodsin AppliedMechanicsandEngineering，2004,193(2004)：4511-4539.

〔3〕A.Sutradhar,G.H.Paulion,L.J.Gray.Transientheatconductioninhomogeneousand non-homogeneousmaterialsbytheLaplace transformGalerkinboundaryelementmethod [J].EngineeringAnalysiswithBoundaryElement,2002,26(5):119-132.

〔4〕李慧.三維非穩態熱傳導邊界元方法研究及數值系統開發[D].山東科技大學,2011.

TB131；O414

A

1673-260X（2013）06-0021-02