基于時域邊界元法的散熱結構瞬態熱傳導分析

周楓林，袁小涵，余江鴻，欽 宇，潘先云

（湖南工業大學 機械工程學院，湖南 株洲 412007）

1 研究背景

散熱是設備熱管理的核心環節，高效的散熱機制是高精密設備運行效能的重要保證，不適的工作環境會產生諸多不良影響，例如：1）電子芯片中，隨著其計算能力的不斷提升，功耗大幅增加，散熱性能已經成為制約其進一步提升的瓶頸；2）新能源汽車中，鋰電池包作為汽車的主要儲能元件，過高的溫度不僅會加速電池老化，甚至可能導致電池熱失控，引發電芯燃燒爆炸事故；3）航天飛行器中，太陽直射面會產生高溫，如果沒有配備平衡溫差的散熱器，溫差過大會導致硬件勞損，危及航天安全；4）醫療儀器的精密度和集成度相對較高，設備內部高溫容易加速元器件老化，出現性能不穩定或功能失效的故障現象，影響設備的精準度和可靠性。散熱器結構中包含許多散熱鰭片，其主要工作原理就是增加散熱面積，快速散熱，有效的熱分析有利于改進散熱系統，使設備在合適的溫度環境中運作。

有限單元法（finite element method，FEM）是一種運用非常廣泛的數值分析方法，其采用體積離散技術，將計算域劃分為有限個單元，單元與單元之間互不重疊，但是互相連接；然后在每個單元內選擇基函數，用每個單元中的基函數來逼近解析解，得出近似解，從而給出整個計算域的解。王瑞等[1]對汽車散熱器進行了有限元分析；龍俊華等[2]運用Hyperworks軟件對某公司散熱器支架進行了仿真分析；王裕林等[3]對某型飛機LED（light-emitting diode）航行燈散熱器進行了溫度散熱分析；王金龍等[4]通 過ANSYS軟 件 對CPU（central processing unit）水冷散熱器進行了數值模擬，并根據分析結果對結構進行了優化；劉貝等[5]利用有限元軟件對翅片式相變散熱器進行了仿真分析，并根據仿真結果對散熱器結構進行了優化設計，改善了其散熱性能。

上述研究中使用的軟件，其底層算法框架都是基于有限單元法[6-8]，雖然有限單元法已發展成熟，是一種可靠的數值分析方法，但有限單元法需要對整個模型進行離散，計算數據龐大，計算耗時長。隨著傳熱學的發展，數值計算方法有了極大的進步，越來越多的方法被用于傳熱學問題，常用的有邊界單元法（boundary element method，BEM ）[9-11]、有限差分法[12-13]、無網格法[14]，其中，邊界單元法只在研究區域的邊界上剖分單元，涉及的計算域較小，網格劃分較少，除了能處理有限單元法所適應的大部分問題外，還能處理有限單元法不易解決的無限域問題。如徐闖等[15]運用邊界元法對功能梯度材料進行了熱傳導分析；肖雄等[16]對正交各向異性功能梯度材料的瞬態熱傳導問題進行了邊界元分析，證明了邊界元法的可靠性；Yu H.P.等[17]運用等幾何邊界元法對電子封裝結構進行了熱分析。但對于邊界元法在三維瞬態熱傳導問題上的工程運用研究尚不多見，還有待進一步研究。故本研究將運用雙互易邊界元法和精細積分方法耦合來分析散熱結構。

本文是在熱力學理論基礎上，對散熱器結構進行熱傳導分析，推導了不含內部熱源的各項同性介質瞬態常系數熱傳導問題的邊界積分方程，將該積分方程運用到散熱器的計算邊界上，然后使用精細積分法求解邊界離散后得到的常微分方程組。經過與有限元計算結果的對比，驗證了邊界元方法的準確性，并具有計算量小的優點。

2 雙互易精細積分邊界單元法

散熱器的熱傳導分析中，散熱器的材料大部分為鋁，因而本文的工況條件設定如下：熱傳導系數不隨溫度變化，工作過程中不含內部熱源，材料為各向同性介質。于是，在其三維瞬態熱傳導問題中，其控制方程可表示如下：

邊界條件設定如下：

式中：q為熱通量；為邊界上的已知熱通量；Su表示狄利克雷型邊界；Sq表示諾依曼型邊界；Su+Sq=Γ，為計算域Ω的所有邊界；n為邊界的單位外法線方向向量。

初始條件如下：

對于不含內部熱源的三維瞬態常系數熱傳導問題的控制方程，運用加權余量法推導其邊界積分方程。引入權函數，則式（1）的加權余量式為

根據高斯散度定理和三維位勢問題的基本解，運用分部積分法將式（4）中的左邊域積分轉化為如下內部點積分方程：

當源點位于邊界時，需要對邊界進行拓撲，對于三維情況，可以得到如下邊界-域積分方程：

式中：

其中θ為邊界點處的切面所圍成的立體角。

選擇如下RBF：

式中：r為RBF插值點到源點的距離；s為形狀參數；i為總插值點數。

令

式中F為徑向基函數矩陣。

將式（8）代入式（6），得

令

將式（11）代入式（10）中，再次運用高斯散度定理和三維位勢問題的基本解，可以得到如下邊界積分方程：

式中：

在邊界積分方程式（12）中，不再涉及域積分，其中溫度對時間的一階導數項已被等效邊界積分所代替。與有限元法不同的是，因為未知量u和q都是在邊界上取值，故只需對邊界進行離散，可得：

式中：j為邊界單元個數；m為單元中的節點數，本文采用8節點二次面單元，m=8。

將源點遍及所有的邊界節點，可以得到如下矩陣形式：

式中：H和G為影響系數矩陣；矩陣上方的字母為矩陣維數，其中b為邊界節點個數，d為域內節點個數。

將式（9）代入式（15），可得：

改寫為如下矩陣形式：

考慮如下邊界條件：

式中，C1和C2為對角矩陣，對角元素的取值根據邊界條件的類型確定如下：

將式（18）引入式（17），可得：

則邊界節點上的量可表示為

式中：

此時考慮內部點，將式（16）改寫為如下形式：

將式（21）代入式（22），可得內部節點上的量：

式中：

精細積分法可用于求解該一階常系數微分方程組，式（23）的通解為

式中：Δt=tk-tk-1；tk=kΔt。

矩陣指數函數

可被細分為

式中，M為細分參數。

式中：P為截斷參數；Er為式（28）中I的后P項。則式（27）可改寫為

再次細分，通過循環語句重復計算和重新賦值，最終可得

3 山型散熱器結構熱傳導分析

本文分析用山型散熱器結構的簡化模型及尺寸（單位為m）如圖1所示。因為散熱器材料主要為6063鋁，其熱傳導系數不隨溫度而變化，故假定其熱導率、熱容和密度分別為1 W/(m·℃)，1 J/(kg·℃)和1 kg/m3。

圖1 山型散熱器結構及其尺寸示意圖Fig.1 Schematic diagram of structure and size of the gable radiator

3.1 有解析解工況下的BEM

對于前文所提到的雙互易精細積分法（precise integration-dual reciprocity boundary element method，PI-DRBEM），本節將給出一個算例，通過將PIDRBEM的結果與解析解進行對比，驗證其準確性和有效性。

定義如下相對誤差：

對于細分參數M和截斷參數P，本文分別取10和6。

模型的邊界劃分如圖2所示。

圖2 模型邊界劃分示意圖Fig.2 Boundary division diagram of the model

本算例中，Γ1和Γ2為諾依曼型邊界，其條件假定如下：

其余邊界為狄利克雷型邊界，其條件假定為

本算例不含熱源，初始條件給定為

此時根據控制方程式（1）可得解析解為

RBF插值點分布如圖3所示。

圖3 RBF插值點分布示意圖Fig.3 Distribution diagram of RBF interpolation points

整個計算域被劃分為180個矩形單元，總共800個邊界節點，92個RBF插值點?？紤]不同時間步長下域內溫度與解析解的相對誤差情況，所得結果如圖4所示。

圖4 溫度隨時間變化的相對誤差曲線Fig.4 Relative error curves of temperature with time

圖4中共給出了5種不同時間步長下，溫度隨時間變化的相對誤差曲線，由圖可知，當時間步長為0.1 s時，其相對誤差曲線位于最下方，表明同一時刻的相對誤差值最小，此時，即使是在接近初始時刻的第5 s，其相對誤差值也很小，具有較高的精確度。并且隨著時間的推進，相對誤差值越來越小，數值解也越趨近于解析解，說明BEM在分析山型散熱結構的熱傳導問題時是有效可靠的。表1給出了點A(0.600 0, 0.916 7, 0.625 0)在不同時間步長下的溫度情況，并給出其對應的溫度解析解。

表1 A點在不同時間步長下的溫度情況Table 1 Temperature of point A at different time steps ℃

由表1可以得出，在所給出的5種不同時間步長下，其所求溫度值與解析解的最大差值為2.07 ℃，出現在時間步長為2.5 s時的第50 s時刻；最小差值為0.959 ℃，出現在時間步長為0.1 s時的第5 s時刻。可見，所有時間步長下的數值計算結果和解析解幾乎吻合，進一步驗證了BEF在山型散熱器模型分析中的精確性。

3.2 無解析解工況下的BEM與FEM對比

通常情況下很難獲得系統的解析解，故本小節將在無解析解的工況下，將FEM分析軟件（Workbench）的結果與BEM結果進行對比分析。假定散熱器底面的溫度變化情況為與時間相關的二次函數，以近似于實際工況下的溫度變化情況。

底部邊界條件為狄利克雷型邊界，其條件如下：

無熱源項時，初始條件給定為常溫：

有限元分析模型的網格劃分如圖5所示。

圖5 山型散熱結構的有限元網格劃分示意圖Fig.5 Schematic diagram of finite element mesh division of mountain heat dissipation structure

圖5所示結構中，單元尺寸定義為0.05 m，共劃分為24 320個線性六面體網格，共有118 217個節點；時間步長取1 s。

對于邊界元法，時間步長同樣取1 s，然后在上一個算例的基礎上，加密邊界網格劃分，將整個計算域劃分為500個矩形單元，節點數為1 920個，RBF插值點數為164個。可以明顯看出，不論在單元個數還是節點個數上，邊界元法都遠遠少于有限元法，故其計算量相對低很多。

對于這兩種方法，考慮同一點B(0.75, 0.91, 0.52)，該點溫度在兩種數值方法的計算下，隨時間的變化情況如圖6所示。

圖6 B點PI-DRBEM和FEM的計算溫度對比曲線Fig.6 Calculated temperature comparison curve of PIDRBEM and FEM at point B

由圖6可知，在開始時刻附近，PI-DRBEM的溫度計算結果和FEM的溫度計算結果間有一點出入，但隨著時間推移，兩結果趨于吻合，對比表明，邊界元法和有限元法同樣具有很高的可靠性和精確性，并且具有計算量更小的優點。

4 結語

對于山型散熱器結構的瞬態熱傳導問題，本文運用了雙互易方法與精細積分方法耦合進行計算求解，并在理想工況下與解析解進行對比，在實際工況下與FEM的分析結果進行對比，兩種工況下的對比結果都驗證了PI-DRBEM的準確性和有效性。在與FEM進行分析對比時，可以很清楚地得知邊界元法具有計算量小、求解所需內存小、計算效率高的特點，同時還能保證計算結果的精確性。

本研究論證的邊界元法，對于散熱器結構分析來說是一種有效的方法，對于所有機械設備的散熱分析有很大的幫助，有利于改善散熱器結構布局，提高設備運作穩定性。