三維“體點”散熱問題建模仿真分析

徐仲強，王遠弟

(上海大學數學系，上海 200444)

1 引言

隨著電子元器件尺寸越來越小，功率也不斷攀升，使得元器件散熱變得愈發重要，元器件的傳熱優化問題受到廣泛關注[1-3]。

現在這個問題的有效解決辦法之一是在電子元器件內部構造高傳熱材料傳導路徑，通過插入高傳熱性材料(如金剛石或碳纖維)，可以更有效地排出內部的熱量。如何利用有限數量的高傳熱性材料來構造熱傳導路徑，將平均溫度最小化轉化為一個優化問題，在熱傳導優化中，稱為“體點”問題(VP)。“體點”問題是熱傳導優化的基礎問題，最初由Bejan定義。

Bejan等人于1996年[4]提出了基于構型理論的樹狀分叉網絡，討論了其在微電子器件冷卻散熱優化方面的應用，對平面體點問題進行了優化。他的基本思路是從最小的二維面積單元開始(通常與工藝上能達到的最小制造單元相當)，對該面積單元利用優化方法使其熱阻最小，達到其最優外形；然后利用優化方法對這個最小的面積單元進行第一次組合，得到第一次組合體結果，從第二次組合體開始，新的組合體一定包含兩個較小的兩個優化后的組合體。按照這樣的優化方法進行下去，直到經過若干次的組合優化后得到的最優組合體能覆蓋住所給的面積。Bejan的構形理論建立在主次干道相互垂直的基礎上，并對每一干道的寬度和長度進行優化，得到了在內熱源均勻、高傳熱材料和基體材料的傳熱系數比值較大的情況下高導材料最優布置方案。

過增元等人以提高傳熱效率為優化目標，基于生命演化的自然原理提出并發展了仿生優化方法[7]，定義了火積(描述熱傳導能力的物理量)[5]、熱量傳遞勢容和熱量傳遞勢容耗散函數[6]，提出了最小熱量傳遞勢容耗散原理，得出了當傳熱系數最佳分布時，全場的溫度梯度應處處相等。仿生優化方法所遵循的梯度均勻化原則的理論推導利用了傳熱系數的連續性[6]，對內熱源均勻且高傳熱材料和基體材料的傳熱系數比值較小、非均勻內熱源[8]、具有相變的非穩態狀態[9]等問題都有較好的優化結果。

2018年，王遠弟和張俊頂等人[10]在討論二維傳熱問題時，以全場溫度均值最低為優化目標，提出并證明了高傳熱系數材料最優分布時的曲面面積極小化原則，最低全場溫度均值時高傳熱材料布置對應最小溫度場曲面面積，利用數值模擬說明了極小曲面法的合理有效性，并對比了極小曲面法和仿生優化方法在不同條件下的優化效果。

受二維仿生優化以及極小曲面啟發，本文的優化方向分別是提升傳熱效率和最小化全場溫度均值，對于三維“體點”問題，得出了三維高導材料填充標準。在此基礎上對兩種方法填充，對比剖析兩種方法模擬出的構造，可知三維模型的填充結構與二維類似，受到高傳熱材料傳熱比以及填充量的影響，得到的三維模型高傳熱材料結構對于實際問題具有借鑒意義。

2 數學模型

底部中心模型如圖1所示。

圖1 底部中心模型

圖1左圖模型a和右圖模型b描述了一個空間元器件材料“體點”傳熱問題。立方體邊長為L，內均勻填充熱源q。高導材料和初始材料的傳熱系數分別為Kp和K0，高導材料填充量為Vp，填充率為η(Vp/V)。模型a底部有長度為δ(δ?L)的正方形開口，邊界絕熱，開口溫度恒定U0。模型b中心鏤空，鏤空處保持恒定溫度U0，邊界絕熱。 求區域內高傳熱系數材料最優填充，使得全場溫度均值最低。一直處在運行中的器件，可以視為穩態，所以考慮平衡態下的傳熱問題。設u=u(x，y，z)為溫度函數，為了最小化全場溫度均值。記

(1)

為全場溫度均值。其中|Ω|為空間區域Ω的度量即體積，易知|Ω|恒定，省略它對結果無影響。于是三維“體點”傳熱問題的數學模型a、b可化簡為

(2)

在邊界面各個方向上，對于底部開口的模型a：

(3)

對于中間鏤空的模型b

(4)

其中Γ11、Γ21、Γ22、Γ23分別為模型a的開口、前后絕緣壁、左右絕緣壁和上下絕緣壁，Γ12、Γ24、Γ25、Γ26分別為模型b的開口、前后絕緣壁、左右絕緣壁和上下絕緣壁，式 (2)為優化目標，式(3)(4)分別為模型a、b的邊界條件。取傳熱口處不變溫度為u=0，否則做平移變換即可得。

3 理論分析

這里討論兩種三維模型填充方法，即仿生優化和極小曲面方法。注意到模型a和b僅僅只是在中間溫度的區別，外部邊界條件類似，所以這里以模型a為例作理論推導。

3.1 極小曲面方法

數學上，把R3中平均曲率為零的曲面稱為極小曲面，該函數u=u(x，y)滿足偏微分方程

(5)

相應地，三維空間上的函數u=u(x，y)滿足極小曲面方程

(6)

(7)

(8)

計算導數

(9)

其中v={cosα，cosβ，cosγ}是?Ω的外法向量，由w∈W的任意性，可得

(10)

由式(10)不難計算出當全場溫度均值最低時，溫度函數滿足極小曲面方程。

(11)

3.2 仿生優化方法

對于三維VP問題，同樣引入了火積[6]

Zdis=?k|?u|2dA

(12)

其中，Zdis表示熱量傳遞過程中的火積耗散，k為傳熱系數，u=u(x，y，z)是傳熱區域Ω上平衡態下的溫度函數， ?u為溫度梯度。火積耗散越小，則溫度梯度場越均勻，給定區域的溫度也越低，因而傳熱優化程度就越高。因此，在對元器件模型作網塊劃分后，高傳熱材料應該首先放在火積最大的微元位置，這些點相當于火積最“突出”的點，當把高導材料填充在這些位置點時，必然最大程度地減少整體火積。

3.2.1 三維球狀元件的散熱問題

最小熱量傳遞勢容耗散原理指出傳熱系數或密度為最佳徑向分布時，熱量傳遞勢容耗散最小，即滿足

(14)

3.2.2 三維體點散熱仿生優化

設u=u(x，y，z)是傳熱區域Ω上的溫度函數，以火積耗散最小為優化目標，只需分析min?Ωk|?u|2dV，這里的|Ω|指空間區域Ω 的度量即體積。滿足第2節中模型a的邊界條件。這里dV是體積微元，u=u(x，y，z)，u=u(x，y，z)∈

kuxx+kuyy+kuzz+kxux+kyuy+kzuz=0

(15)

當k為常數時，溫度函數u=u(x，y，z)滿足拉普拉斯方程。

4 高傳熱系數材料優化準則

4.1 極小曲面

根據3.1的理論推導可知，最低全場溫度均值時溫度函數圖像對應的超曲面是“極小曲面”。因而，三維超曲面極小化原則指明高導材料最優布置應該使得溫度函數圖像的“面積”趨于極小化。由于每一塊高導材料的放置都會影響先前的溫度場，溫度函數圖像更新。故極小曲面原則即要求不斷地極小化溫度函數圖像的“面積”。

然而溫度函數圖像“面積”極小化難以量化，高導材料的最優布置只通過極小曲面原則是無法明確的，所以需要進一步簡化超曲面“面積”極小化原則。

數值模擬過程中利用含有源項的穩態熱傳導方程

?(K(x，y，z)?u)+q=0

(16)

設內熱源為q=10000W/m3，圖2左圖是初始溫度場，其中鄰接網塊點間溫差小于1K。圖2右圖是填充高導材料后的溫度場，存在溫度為338.5959K的網塊點，其鄰接上、下、左、右、前、后六個網格點上的溫度分別是338.9062K，338.3786K，338.6786K，338.4355K，338.6786K和338.4355K。即ux<1，uy<1，uz<1所以有uxuyuxy?1，uxuzuxz?1，uyuzuyz?1。這是鄰接網塊點間溫差小于的一個例子，填充后溫度場各網塊均符合此現象。易知填充高導材料后溫差變小，式(11)表達的極小曲面方程可近似看作調和方程

uxx+uyy+uzz=0

(17)

圖2 填充前模型a溫度場(左圖)

模型a極小曲面法傳熱比為 300，填充率為8%溫度場(右圖)

通過將最低均值溫度時的溫度函數近似為調和函數，便可利用調和函數的平均值定理[11]，即

(18)

上式為調和函數平均值定理滿足的必要條件，記

(19)

4.2 仿生優化

仿生優化方法[7]將高導材料的布置分為進化與退化。在溫度梯度最大位置進化，梯度最小處退化淘汰，溫度梯度趨向均勻化。對(12)式作離散化處理，對立方體元器件作步長為n的分割后，目標函數轉化為

(20)

圖3 (i，j，z)微元方塊與其右側方塊界面

對于微元方塊(i，j，z)處，其梯度模的二次方

(21)

對于內部微元方塊，比如圖3中微元方塊(i，j，z)，考察其右界面，根據傅里葉定律(Q=-kAdT/dx)有

(22)

即

(23)

同理其它非邊緣界面上相應溫度

對于邊緣界面，比如右邊緣微元方塊，其右邊界溫度Tw=T(i，j，z)。

5 數值計算

5.1 數值計算細化

因為模型a和b區別僅僅在于開口構造不同，數值計算本質上原理趨同，不失一般性，這里以模型a為例闡明數值模擬細則。取基體材料傳熱系數處處為1，開口邊長占元器件邊長的1/5，開口溫度為300K。

極小曲面方法根據式(19)，采用中心點周圍相鄰一層的球體鄰域上的積分平均值來進行數值模擬。數值模擬 (參考圖19) 可知，按球體鄰域積分均值所得全場溫度均值效果優于正方體鄰域，這與球體鄰域相鄰方塊更能代表中心體溫度信息有關，正方體鄰域中，四角相對較遠的方塊會沖淡相鄰方塊含有的中心體溫度變化信息，從而使得正方體鄰域效果相對不如球體鄰域，所以后面以球體鄰域進行計算。

填充過程中，利用對稱性，空間上每個高導材料填充點有四個對稱點。于是，三維“體點”傳熱問題的數值模擬過程可細化為：

1)預定義參數，包括傳熱系數，網塊步長；

2)根據(16)式求溫度場散布；

3)利用(18)式求得中心微元體與其球體鄰域溫度均值的差值(絕對值)場；

4)在差值最大處放置高導材料，更新傳熱系數值；

5)重新利用步驟2)更新溫度場，當達到填充比例η時終止填充。

仿生優化方法中，整個模擬計算過程具體操作步驟不同之處在于上述步驟3) 4)，仿生優化方法中，特殊的操作步驟為：

3′)根據式(20)計算每個微元體的火積，得到火積場；

4′)根據放置原則，由火積場選擇火積最大的微元塊作為高傳熱系數材料最優填充位置，并改變選取位置的傳熱系數值。

5.2 優化結果

本文是基于模型a和b采用極小曲面方法和仿生優化方法的數值模擬。三維模型相比二維模型，計算量激增，具體到式(11)求解溫度場的線性方程組系數矩陣從二維的n2×n2躍升為三維的n3×n3維度(其中n表示每一邊長剖分后的網格/塊數)。

下文分析中，導熱系數比值較大以高導熱系數材料與基體材料導熱系數比值等于300為例，導熱系數比值較小以導熱系數比值等于4為例，填充率取8%。

5.2.1 模型a

模型a極小曲面

傳熱系數比值比較大時，圖2左圖表示初始溫度場，溫度最高值為344.6795K，在上方四角，全場溫度均值為T=340.4682K。圖2右圖表示極小曲面原則優化后溫度場，全場溫度均值為T=301.0257K，同比下降39.4425K；最高溫度為301.5658K，同比下降43.1137K。可見高導材料放置后，最高溫度和全場溫度均值都有明顯降低。

圖4 模型a極小曲面法傳熱比為 300，填充率為 8%三視圖：主視圖(左上)、左視圖(右上)、俯視圖(左下)、三維圖(右下)

傳熱系數比值比較小時：優化結果如圖6所示，圖5為相應的溫度場散布圖像。優化后全場溫度均值T=314.0793K，傳熱區域內最高溫度為317.654K。

模型a仿生優化

傳熱系數比值比較大時：梯度均勻化原則優化結果如圖7所示，優化后全場溫度均值僅為T=301.0739K，傳熱區域內最高溫度為301.7303K。

圖5 模型a極小曲面法傳熱比為 4，填充率為 8%溫度場

圖6 模型a極小曲面法傳熱比為 4，填充率為 8%三視圖：主視圖(左上)、左視圖(右上)、俯視圖(左下)、三維圖(右下)

圖7 模型a仿生優化法傳熱比為 300，填充率為 8%三視圖：主視圖(左上)、左視圖(右上)、俯視圖(左下)、三維圖(右下)

傳熱系數比值比較小時：優化結果如圖8所示，優化后全場溫度均值T=313.9304K，最高溫度為317.5472K。

圖8 模型a仿生優化法傳熱比為 4，填充率為 8%三視圖：主視圖(左上)、左視圖(右上)、俯視圖(左下)、三維圖(右下)

5.2.2 模型b

模型b極小曲面

傳熱系數比值比較大時 圖9左圖是模型b填充前的溫度場，由于模型b中間鏤空開口溫度恒定，故未填充之前Z 軸方向上每一層XY平面溫度相等，且從內向外溫度越來越高，模型b最高溫度出現在到散熱豎直開口距離最遠的四角處，大小為305.6639K，全場溫度均值為T=304.3392K，注意到模型a優化前溫度最高處大小為344.6795K，全場溫度均值為T=340.4682K，因為模型b鏤空開口比較大，所以比模型a的散熱效果要好很多，在不填充傳熱材料時也能達到較低的體平均溫度；另外，由于模型b未填充前每一層溫度完全相同，這時極小曲面方法會有短暫的失效，這時根據其它填充經驗，第一個點填充在開口附近便能得到比較理想的填充效果，于是便采用第一點(準確地說是根據對稱性預填充四個點)預填充人為得改變初始溫度場，這時候極小曲面方法便可在此基礎上達到最終填充完整效果。填充效果如圖10所示，圖9右圖是模型b填充后的溫度場。根據極小曲面原則優化后全場溫度均值僅為T=301.5272K，傳熱區域內最高溫度為301.9917K。

傳熱系數比值比較小時：優化結果如圖11所示，圖12是模型b填充后的溫度場。根據極小曲面優原則優化后全場溫度均值為T=303.3659K，最高溫度為304.6382K。通過填充效果圖發現，這時候的填充效果并不是十分理想，在后面的分析中，會對比仿生優化方法說明模型b極小曲面方法局限性。

圖9 模型b填充前溫度場(左圖)、模型b極小曲面法傳熱比為 300，填充率為 8%溫度場(右圖)

圖10 模型b極小曲面法傳熱比為 300，填充率為 8%三視圖：主視圖(左上)、左視圖(右上)、俯視圖(左下)、三維圖(右下)

圖11 模型b極小曲面法傳熱比為 4，填充率為 8%三視圖：主視圖(左上)、左視圖(右上)、俯視圖(左下)、三維圖(右下)

圖12 模型b極小曲面法傳熱比為 4，填充率為 8%溫度場

模型b仿生優化

傳熱系數比值比較大時 優化結果如圖13所示，優化后全場溫度均值T=301.1784K，最高溫度為301.6163K。

圖13 模型b仿生優化法傳熱比為 300，填充率為 8%三視圖：主視圖(左上)、左視圖(右上)、俯視圖(左下)、三維圖(右下)

圖14 模型b仿生優化法傳熱比為 4，填充率為 8%三視圖：主視圖(左上)、左視圖(右上)、俯視圖(左下)、三維圖(右下)

傳熱系數比值比較小時 優化結果如圖14所示，優化后全場溫度均值T=302.8593K，最高溫度為304.0459K。

5.3 結果分析

以上僅僅對部分結果進行了展示，本節重點對結果進行分析：

5.3.1 兩種方法對比分析

在5.2節中，展示了兩個模型兩種方法的填充率為8%，傳熱比分別為300 和4 時候的情況，對于其它傳熱系數比值以及填充比下兩個模型兩種方法的模擬對比見下文分析。

優化效果(體平均溫度變化曲線)

模型a：圖15～圖18分別描述了傳熱比在2～400，1～126，127～269，270～400，填充率在3%～8%時兩種方法的溫度均值差，這里橫軸為傳熱比，縱軸為溫度均值差。虛線以下部分表示極小曲面方法所得溫度均值更低。

圖15 傳熱比1～400，填充率3%～8%下兩種方法計算的溫度均值差值曲線

圖16 傳熱比1～126，填充率3%～8%下兩種方法計算的溫度均值差值曲線

圖17 傳熱比127～269，填充率3%～8%下兩種方法計算的溫度均值差值曲線

圖18 傳熱比270～400，填充率3%～8%下兩種方法計算的溫度均值差值曲線

從圖15～圖18四個圖中的計算結果可以看出，不同填充條件下，極小曲面和仿生優化方法互有所長：①圖15表明優化曲線穩定后，傳熱比值在170～400，填充率為5%～8%時，極小曲面原則效果較好。②如圖16所示，當傳熱比值偏小，大約低于126K時仿生優化效果更佳。③圖17和圖18表明，除了當傳熱系數比值較大(127K以上)時填充率為3%和4%時整體上仿生優化效果好；填充率為5%～8% 時，傳熱比127～269時，兩種方法優化效果基本相同，傳熱比270～400時，極小曲面方法優化效果優于仿生優化方法，兩者的最大全場溫度均值差值在0.1K以內。

模型b：圖19表示傳熱系數比值在2～400之間，填充率為3%～8%時極小曲面方法采用球形鄰域與方形鄰域的溫度均值差值曲線，顯然，傳熱系數比值在50以上時，球形鄰域優化效果要比方形鄰域好，最大差值達到了0.5K，所以在取平均值時選擇的鄰域為球形鄰域；圖20表示傳熱系數比值在2～400之間，填充率為3%～8%時極小曲面方法采用球形鄰域極小曲面方法與仿生優化方法的溫度均值差值曲線，由圖可見，模型b中，仿生優化效果優于極小曲面方法，且隨著高導材料比值越大，兩者差異越小。

圖19 球形鄰域方形鄰域溫度均值差值曲線

圖20 球形鄰域仿生優化溫度均值差值曲線

公共填充率

在本節中，把兩種方法優化結果共同填充點在總填充點中的比例作為公共填充率，公共填充率描述了兩種方法共同填充點，在某種程度上，可以作為填充優化的置信度。

模型a：圖21表示傳熱系數比值在1～400之間，填充率為3%～8%時極小曲面方法和仿生優化方法優化結果中填充共同點的比例，由圖可見，當傳熱系數大于一定程度時，共同填充率趨于平穩；圖22表示傳熱系數比值在300～400之間，填充率為3%～8%時極小曲面方法和仿生優化方法優化結果中填充共同點的比例，共同填充率在30%～65%之間，且填充率為5%時共同填充率相對最高。

圖21 模型a傳熱比1～400，填充率3%～8%兩種方法填充共同點比例

圖22 模型a傳熱比300～400，填充率3%～8%兩種方法填充共同點比例

模型b 圖23表示傳熱系數比值在1～400，填充率為3%～8%時極小曲面方法和仿生優化方法優化結果中填充共同點比例，由圖可見，整體上共同填充率比較平穩，填充率為3%時，兩種方法會間斷出現無共同填充點現象，這與填充比較低時填充點少有一定關系；圖24表示傳熱系數比值在300～400，填充率為3%～8%時極小曲面方法和仿生優化方法優化結果中填充共同點比例在3%～40%之間，且填充率為7%時共同填充率相對最高。

圖23 模型b傳熱比1～400，填充率3%～8%兩種方法填充共同點比例

圖24 模型b傳熱比300～400，填充率3%～8%兩種方法填充共同點比例

由公共填充率分析可知兩個模型兩種方法優化結果的共同填充率比較平穩，整體上在20%以上，結合兩種方法的優化效果，可以相信兩種方法都能達到比較理想的優化效果。

5.3.2 20-50維度對比分析

圖25左右圖分別表示模型a傳熱系數比為18，填充率為3%時極小曲面方法和仿生優化方法方法的填充效果，此時，兩種填充微元均集中在開口附近，這與20 維度三維效果一致。

圖25 模型a傳熱比為 18，填充率為3%極小曲面法填充(左圖)、仿生優化法填充(右圖)

圖26左圖表示模型a傳熱系數比值為300，填充率為5%時極小曲面方法的填充效果，圖26右圖表示模型b傳熱系數比值為150，填充率為5%時極小曲面方法的填充效果，此時，兩種填充微元均比較分散，呈細長分布狀，這與20 維度三維效果一致。

圖26 模型a傳熱比為300，填充率為5%極小曲面法填充(左圖)、模型b傳熱比為150，填充率為5%極小曲面法填充(右圖)

5.3.3 2D-3D對比分析

圖27 二維傳熱系數比為300，填充率為8%極小曲面法填充(左圖)、仿生優化法填充(右圖)

圖27左圖和右圖分別表示二維傳熱系數比值為300，填充率為8%時極小曲面和仿生優化方法填充，對比圖4以及圖7發現，傳熱比值偏大情況下，三維空間上高導材料有類似二維優化結果的長條，這與高傳熱系數材料能顯著改變區域內的傳熱效果有關，長條形態使熱量沿著高傳熱長條更快得流向散熱口處。

圖28正上圖、左下圖、右下圖分別表示三維模型a傳熱比為300，填充率為8% 極小曲面法50 維度填充結果的Z=1半橫切面、X=23縱切面以及Y=3縱切面，可見，在三維空間元器件的二維切面上，也有類似二維優化結果的細長分布長條，這在數值上說明三維極小曲面和仿生優化理論與二維結論有一定相似性。

圖28 模型a極小曲面法傳熱比為300，填充率為8%，Z=1半橫切面(正上圖)、X=23縱切面(左下圖)、Y=3縱切面(右下圖)

圖29 二維傳熱系數比值為4，填充率為8%極小曲面法填充(左圖)、仿生優化法填充(右圖)

圖29左圖和右圖分別表示二維傳熱系數比值為4，填充率為8%時極小曲面和仿生優化方法填充，對比圖6以及圖8發現，當傳熱系數比值較小時，三維空間上高傳熱系數材料有類似二維優化效果，高傳熱材料集中在開口附近。

圖30 模型a極小曲面法傳熱比為4，填充率為8%，Z=2半橫切面(正上圖)、X=18縱切面(左下圖)、Y=25縱切面(右下圖)

圖30正上圖、左下圖、右下圖分別表示三維模型a傳熱比為4，填充率為8%極小曲面法50 維度填充結果的Z=2半橫切面、X=18縱切面以及Y=25縱切面，可見，在三維空間元器件的二維切面上，也有類似二維優化結果的開口集中形狀，直觀得說明了數值上三維極小曲面和仿生優化理論與二維結果有一定相似性。

6 結論

本文借鑒二維“體點”模型研究方法，對三維“體點”問題中高導材料的最佳填充問題進行研究，根據三維極小曲面方法和三維仿生優化理論找到了三維填充優化準則。數值計算結果表明，三維“體點”傳熱問題中高導材料填充問題可以通過極小曲面和仿生優化方法加以解決，有如下結論：

1)三維“體點”最低全場溫度均值時高導材料分布對應最小溫度場曲面面積，即極小曲面的面積。

2)三維“體點”傳熱問題中高導材料的填充準則與二維最優分布填充準則一致，可根據超曲面面積極小化原則(極小曲面原則)以及火積耗散最小原則(仿生優化原則)來確定。

3)三維“體點”傳熱問題中高傳熱系數材料最優區域分布與二維“體點”傳熱問題中高傳熱系數材料最優區域分布具有很大的相似性。從三維立體結構規律上看，極小曲面和仿生優化方法得到的高導材料最優填充相像。傳熱比值偏大時高導材料趨向長條狀；傳熱系數比值偏小時高導材料開口處聚集。從三維“體點”的二維切面上看，高傳熱系數材料最優區域分布中存在一些切面與二維“體點”傳熱問題中高傳熱系數材料最優區域分布具有很大的相似性，滿足二維分布規律。

4)通過對比可知，三維“體點”傳熱問題中超曲面面積極小化方法和仿生優化方法對“體點”傳熱問題的最優分布有一定比例的共同點，且兩種方法的優化效果各有優勢。當開口豎直鏤空時，仿生優化方法優化效果好；當開口僅存在底部時，填充率為5%～8%，高傳熱系數比值在127～400時極小曲面方法優化效果較好。