基于時(shí)域邊界元法的散熱結(jié)構(gòu)瞬態(tài)熱傳導(dǎo)分析

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

（湖南工業(yè)大學(xué) 機(jī)械工程學(xué)院，湖南 株洲 412007）

1 研究背景

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

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

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

本文是在熱力學(xué)理論基礎(chǔ)上，對(duì)散熱器結(jié)構(gòu)進(jìn)行熱傳導(dǎo)分析，推導(dǎo)了不含內(nèi)部熱源的各項(xiàng)同性介質(zhì)瞬態(tài)常系數(shù)熱傳導(dǎo)問(wèn)題的邊界積分方程，將該積分方程運(yùn)用到散熱器的計(jì)算邊界上，然后使用精細(xì)積分法求解邊界離散后得到的常微分方程組。經(jīng)過(guò)與有限元計(jì)算結(jié)果的對(duì)比，驗(yàn)證了邊界元方法的準(zhǔn)確性，并具有計(jì)算量小的優(yōu)點(diǎn)。

2 雙互易精細(xì)積分邊界單元法

散熱器的熱傳導(dǎo)分析中，散熱器的材料大部分為鋁，因而本文的工況條件設(shè)定如下：熱傳導(dǎo)系數(shù)不隨溫度變化，工作過(guò)程中不含內(nèi)部熱源，材料為各向同性介質(zhì)。于是，在其三維瞬態(tài)熱傳導(dǎo)問(wèn)題中，其控制方程可表示如下：

邊界條件設(shè)定如下：

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

初始條件如下：

對(duì)于不含內(nèi)部熱源的三維瞬態(tài)常系數(shù)熱傳導(dǎo)問(wèn)題的控制方程，運(yùn)用加權(quán)余量法推導(dǎo)其邊界積分方程。引入權(quán)函數(shù)，則式（1）的加權(quán)余量式為

根據(jù)高斯散度定理和三維位勢(shì)問(wèn)題的基本解，運(yùn)用分部積分法將式（4）中的左邊域積分轉(zhuǎn)化為如下內(nèi)部點(diǎn)積分方程：

當(dāng)源點(diǎn)位于邊界時(shí)，需要對(duì)邊界進(jìn)行拓?fù)?，?duì)于三維情況，可以得到如下邊界-域積分方程：

式中：

其中θ為邊界點(diǎn)處的切面所圍成的立體角。

選擇如下RBF：

式中：r為RBF插值點(diǎn)到源點(diǎn)的距離；s為形狀參數(shù)；i為總插值點(diǎn)數(shù)。

令

式中F為徑向基函數(shù)矩陣。

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

令

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

式中：

在邊界積分方程式（12）中，不再涉及域積分，其中溫度對(duì)時(shí)間的一階導(dǎo)數(shù)項(xiàng)已被等效邊界積分所代替。與有限元法不同的是，因?yàn)槲粗縰和q都是在邊界上取值，故只需對(duì)邊界進(jìn)行離散，可得：

式中：j為邊界單元個(gè)數(shù)；m為單元中的節(jié)點(diǎn)數(shù)，本文采用8節(jié)點(diǎn)二次面單元，m=8。

將源點(diǎn)遍及所有的邊界節(jié)點(diǎn)，可以得到如下矩陣形式：

式中：H和G為影響系數(shù)矩陣；矩陣上方的字母為矩陣維數(shù)，其中b為邊界節(jié)點(diǎn)個(gè)數(shù)，d為域內(nèi)節(jié)點(diǎn)個(gè)數(shù)。

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

改寫(xiě)為如下矩陣形式：

考慮如下邊界條件：

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

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

則邊界節(jié)點(diǎn)上的量可表示為

式中：

此時(shí)考慮內(nèi)部點(diǎn)，將式（16）改寫(xiě)為如下形式：

將式（21）代入式（22），可得內(nèi)部節(jié)點(diǎn)上的量：

式中：

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

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

矩陣指數(shù)函數(shù)

可被細(xì)分為

式中，M為細(xì)分參數(shù)。

式中：P為截?cái)鄥?shù)；Er為式（28）中I的后P項(xiàng)。則式（27）可改寫(xiě)為

再次細(xì)分，通過(guò)循環(huán)語(yǔ)句重復(fù)計(jì)算和重新賦值，最終可得

3 山型散熱器結(jié)構(gòu)熱傳導(dǎo)分析

本文分析用山型散熱器結(jié)構(gòu)的簡(jiǎn)化模型及尺寸（單位為m）如圖1所示。因?yàn)樯崞鞑牧现饕獮?063鋁，其熱傳導(dǎo)系數(shù)不隨溫度而變化，故假定其熱導(dǎo)率、熱容和密度分別為1 W/(m·℃)，1 J/(kg·℃)和1 kg/m3。

圖1 山型散熱器結(jié)構(gòu)及其尺寸示意圖Fig.1 Schematic diagram of structure and size of the gable radiator

3.1 有解析解工況下的BEM

對(duì)于前文所提到的雙互易精細(xì)積分法（precise integration-dual reciprocity boundary element method，PI-DRBEM），本節(jié)將給出一個(gè)算例，通過(guò)將PIDRBEM的結(jié)果與解析解進(jìn)行對(duì)比，驗(yàn)證其準(zhǔn)確性和有效性。

定義如下相對(duì)誤差：

對(duì)于細(xì)分參數(shù)M和截?cái)鄥?shù)P，本文分別取10和6。

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

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

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

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

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

此時(shí)根據(jù)控制方程式（1）可得解析解為

RBF插值點(diǎn)分布如圖3所示。

圖3 RBF插值點(diǎn)分布示意圖Fig.3 Distribution diagram of RBF interpolation points

整個(gè)計(jì)算域被劃分為180個(gè)矩形單元，總共800個(gè)邊界節(jié)點(diǎn)，92個(gè)RBF插值點(diǎn)?？紤]不同時(shí)間步長(zhǎng)下域內(nèi)溫度與解析解的相對(duì)誤差情況，所得結(jié)果如圖4所示。

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

圖4中共給出了5種不同時(shí)間步長(zhǎng)下，溫度隨時(shí)間變化的相對(duì)誤差曲線，由圖可知，當(dāng)時(shí)間步長(zhǎng)為0.1 s時(shí)，其相對(duì)誤差曲線位于最下方，表明同一時(shí)刻的相對(duì)誤差值最小，此時(shí)，即使是在接近初始時(shí)刻的第5 s，其相對(duì)誤差值也很小，具有較高的精確度。并且隨著時(shí)間的推進(jìn)，相對(duì)誤差值越來(lái)越小，數(shù)值解也越趨近于解析解，說(shuō)明BEM在分析山型散熱結(jié)構(gòu)的熱傳導(dǎo)問(wèn)題時(shí)是有效可靠的。表1給出了點(diǎn)A(0.600 0, 0.916 7, 0.625 0)在不同時(shí)間步長(zhǎng)下的溫度情況，并給出其對(duì)應(yīng)的溫度解析解。

表1 A點(diǎn)在不同時(shí)間步長(zhǎng)下的溫度情況Table 1 Temperature of point A at different time steps ℃

由表1可以得出，在所給出的5種不同時(shí)間步長(zhǎng)下，其所求溫度值與解析解的最大差值為2.07 ℃，出現(xiàn)在時(shí)間步長(zhǎng)為2.5 s時(shí)的第50 s時(shí)刻；最小差值為0.959 ℃，出現(xiàn)在時(shí)間步長(zhǎng)為0.1 s時(shí)的第5 s時(shí)刻?？梢?jiàn)，所有時(shí)間步長(zhǎng)下的數(shù)值計(jì)算結(jié)果和解析解幾乎吻合，進(jìn)一步驗(yàn)證了BEF在山型散熱器模型分析中的精確性。

3.2 無(wú)解析解工況下的BEM與FEM對(duì)比

通常情況下很難獲得系統(tǒng)的解析解，故本小節(jié)將在無(wú)解析解的工況下，將FEM分析軟件（Workbench）的結(jié)果與BEM結(jié)果進(jìn)行對(duì)比分析。假定散熱器底面的溫度變化情況為與時(shí)間相關(guān)的二次函數(shù)，以近似于實(shí)際工況下的溫度變化情況。

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

無(wú)熱源項(xiàng)時(shí)，初始條件給定為常溫：

有限元分析模型的網(wǎng)格劃分如圖5所示。

圖5 山型散熱結(jié)構(gòu)的有限元網(wǎng)格劃分示意圖Fig.5 Schematic diagram of finite element mesh division of mountain heat dissipation structure

圖5所示結(jié)構(gòu)中，單元尺寸定義為0.05 m，共劃分為24 320個(gè)線性六面體網(wǎng)格，共有118 217個(gè)節(jié)點(diǎn)；時(shí)間步長(zhǎng)取1 s。

對(duì)于邊界元法，時(shí)間步長(zhǎng)同樣取1 s，然后在上一個(gè)算例的基礎(chǔ)上，加密邊界網(wǎng)格劃分，將整個(gè)計(jì)算域劃分為500個(gè)矩形單元，節(jié)點(diǎn)數(shù)為1 920個(gè)，RBF插值點(diǎn)數(shù)為164個(gè)?？梢悦黠@看出，不論在單元個(gè)數(shù)還是節(jié)點(diǎn)個(gè)數(shù)上，邊界元法都遠(yuǎn)遠(yuǎn)少于有限元法，故其計(jì)算量相對(duì)低很多。

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

圖6 B點(diǎn)PI-DRBEM和FEM的計(jì)算溫度對(duì)比曲線Fig.6 Calculated temperature comparison curve of PIDRBEM and FEM at point B

由圖6可知，在開(kāi)始時(shí)刻附近，PI-DRBEM的溫度計(jì)算結(jié)果和FEM的溫度計(jì)算結(jié)果間有一點(diǎn)出入，但隨著時(shí)間推移，兩結(jié)果趨于吻合，對(duì)比表明，邊界元法和有限元法同樣具有很高的可靠性和精確性，并且具有計(jì)算量更小的優(yōu)點(diǎn)。

4 結(jié)語(yǔ)

對(duì)于山型散熱器結(jié)構(gòu)的瞬態(tài)熱傳導(dǎo)問(wèn)題，本文運(yùn)用了雙互易方法與精細(xì)積分方法耦合進(jìn)行計(jì)算求解，并在理想工況下與解析解進(jìn)行對(duì)比，在實(shí)際工況下與FEM的分析結(jié)果進(jìn)行對(duì)比，兩種工況下的對(duì)比結(jié)果都驗(yàn)證了PI-DRBEM的準(zhǔn)確性和有效性。在與FEM進(jìn)行分析對(duì)比時(shí)，可以很清楚地得知邊界元法具有計(jì)算量小、求解所需內(nèi)存小、計(jì)算效率高的特點(diǎn)，同時(shí)還能保證計(jì)算結(jié)果的精確性。

本研究論證的邊界元法，對(duì)于散熱器結(jié)構(gòu)分析來(lái)說(shuō)是一種有效的方法，對(duì)于所有機(jī)械設(shè)備的散熱分析有很大的幫助，有利于改善散熱器結(jié)構(gòu)布局，提高設(shè)備運(yùn)作穩(wěn)定性。