微納尺度體點導熱的拓撲優(yōu)化*

李含靈 曹炳陽

(清華大學工程力學系,熱科學與動力工程教育部重點實驗室,北京 100084)

體點導熱問題是電子器件散熱優(yōu)化方面的基礎(chǔ)問題之一,已有研究大多建立在傅里葉導熱定律的基礎(chǔ)上,但隨著電子器件的特征尺度降低到微納米量級,導熱優(yōu)化需要考慮非傅里葉效應(yīng).本文結(jié)合聲子玻爾茲曼方程的數(shù)值解和對聲子平均自由程進行插值的固體各向同性材料懲罰方法,發(fā)展了微納米尺度下彈道-擴散導熱的拓撲優(yōu)化方法.在彈道-擴散導熱機制下,體點導熱的拓撲優(yōu)化得到的材料分布明顯不同于擴散導熱下的樹狀分布,且會隨努森數(shù)的變化而變化,與拓撲優(yōu)化的插值方式和聲子彈道輸運有關(guān).隨著彈道效應(yīng)的增強,尺寸效應(yīng)使得材料分布中微小結(jié)構(gòu)的等效熱導率低于粗壯結(jié)構(gòu),因而拓撲優(yōu)化結(jié)果朝著增多粗壯結(jié)構(gòu)、減少微小結(jié)構(gòu)的方向發(fā)展.彈道效應(yīng)足夠強時,填充材料聚集在低溫邊界附近,主干和枝合并,呈團狀分布.

1 引 言

以芯片為代表的電子器件,在當代社會中扮演著不可替代的重要角色.隨著電子器件向小型化、集成化的方向發(fā)展,器件的特征尺寸已經(jīng)降低至微納米量級[1],器件內(nèi)的功率密度大幅增加.2012年,芯片內(nèi)的平均熱流密度便已達到250 W/cm2左右[2].高功率密度會讓電子器件溫度升高,產(chǎn)生局部熱點,導致其可靠性和壽命顯著下降[3].因此,電子器件的熱管理已經(jīng)成為相關(guān)領(lǐng)域發(fā)展面臨的關(guān)鍵挑戰(zhàn)之一[4-7].根據(jù)電子器件內(nèi)熱點分布和傳熱結(jié)構(gòu)布局等因素,進行散熱仿真和優(yōu)化改進,對實際產(chǎn)品的設(shè)計和制造有著重要意義.

在芯片內(nèi)部填充高導熱材料以加強導熱是提高芯片散熱能力的主要途徑,考慮到芯片內(nèi)空間寶貴,高導熱材料的添加數(shù)量必須受到限制,這樣的芯片散熱問題可以被抽象為體點(volume-to-point,VP)導熱問題[8],如圖1所示.給定體積的區(qū)域代表芯片,區(qū)域內(nèi)的均勻內(nèi)熱源代表焦耳加熱,產(chǎn)生的熱量通過邊界處一小段溫度為 T0的低溫熱沉流出,其余邊界絕熱.現(xiàn)提供數(shù)量有限的高導熱填充材料,希望尋找合適的填充材料分布方式,構(gòu)造高導熱通道,讓整個系統(tǒng)的溫度(平均溫度或最高溫度)降到最低.關(guān)于體點導熱問題,人們采用多種優(yōu)化方法進行了研究,包括構(gòu)型理論[8]、仿生優(yōu)化[9]、拓撲優(yōu)化[10-13]、模擬退火和遺傳算法[14]等.盡管這些方法都能有效降低系統(tǒng)溫度、提高散熱能力,但所得的高導熱填充材料分布方式不盡相同,這與優(yōu)化問題自身帶有的局部最優(yōu)性和多解性有關(guān).當優(yōu)化問題改變時,不同方法的性能優(yōu)劣也可能發(fā)生變化.在各類優(yōu)化方法中,拓撲優(yōu)化方法因其能夠改變結(jié)構(gòu)的拓撲構(gòu)型、具備理論上的最高設(shè)計自由度而受到人們的青睞.利用拓撲優(yōu)化,設(shè)計者不僅可以實現(xiàn)設(shè)計對象某方面性能的顯著提升,而且所花費的設(shè)計時間更少,能得到傳統(tǒng)經(jīng)驗之外的具有啟發(fā)性的新結(jié)果[15].體點導熱問題的拓撲優(yōu)化會得到填充材料充分伸入內(nèi)部區(qū)域、整體呈樹狀的材料分布,其散熱性能明顯優(yōu)于其他方法所預(yù)測的結(jié)構(gòu)[11].

圖1 體點導熱問題示意圖Fig.1.The schematic diagram of the VP problem.

已有的導熱拓撲優(yōu)化方法,都建立在經(jīng)典的傅里葉導熱定律的基礎(chǔ)上[13],這在宏觀尺度下是合理的,但對微納米尺度下的導熱過程,傅里葉導熱定律不再成立[16-22],熱仿真及熱優(yōu)化必須考慮非傅里葉效應(yīng).聲子是電子器件半導體中的主要載熱子[23],聲子彈道輸運和聲子相干是微納米尺度下發(fā)生非傅里葉導熱的兩個重要原因[5].聲子彈道輸運的強度可由努森數(shù) Kn 來表征,其定義為材料聲子平均自由程l和系統(tǒng)特征尺寸L的比值,即Kn=l/L.宏觀尺度下,系統(tǒng)尺寸遠大于聲子平均自由程(Kn→0),聲子間散射充分,聲子按擴散方式傳遞,導熱過程遵循傅里葉導熱定律.但當系統(tǒng)尺寸與聲子平均自由程相當(Kn～1)甚至更小(Kn＜1)時,聲子不經(jīng)歷散射而直接抵達邊界,這樣的過程稱為彈道輸運,會導致傅里葉導熱定律失效[24].在微納米結(jié)構(gòu)中,部分聲子發(fā)生彈道輸運,其余聲子仍按擴散方式運動,對應(yīng)的導熱機制稱為彈道-擴散導熱[25,26].聲子彈道輸運會引起熱導率的尺寸效應(yīng),表現(xiàn)為系統(tǒng)的等效熱導率與系統(tǒng)尺寸正相關(guān),尺寸減小會造成等效熱導率降低[27,28],系統(tǒng)熱點溫度升高[29].此外,在給定溫度邊界處,由于彈道輸運的非平衡本質(zhì),邊界處的溫度并不連續(xù),會出現(xiàn)溫度跳躍現(xiàn)象[30,31].聲子相干則主要出現(xiàn)在特征尺寸和聲子波長相當?shù)南到y(tǒng)中[32],例如石墨烯這樣的二維材料和超晶格這樣的特殊結(jié)構(gòu)[33-35].聲子相干通過改變聲子群速度、弛豫時間和態(tài)密度等因素來影響導熱過程[36].對常用的半導體材料硅,室溫下聲子平均自由程約為300 nm[37],聲子波長小于10 nm[17],而目前大多數(shù)電子器件的特征尺寸在100 nm量級[1],故本文主要關(guān)注聲子彈道輸運引起的非傅里葉導熱.關(guān)于微納米尺度下彈道-擴散導熱過程的優(yōu)化,前人曾結(jié)合動力學理論和拓撲優(yōu)化來處理[38],但其優(yōu)化目標只是特定點之間的溫差,無法考慮系統(tǒng)整體的溫度.目前,尚缺少適用于微納米尺度下彈道-擴散導熱過程的一般性優(yōu)化方法,關(guān)于彈道效應(yīng)對優(yōu)化結(jié)果的影響也缺乏系統(tǒng)的認識.

本文結(jié)合聲子玻爾茲曼輸運方程(Boltzmann transport equation,BTE)和固體各向同性材料懲罰(solid isotropic material with penalization,SIMP)方法,發(fā)展了微納尺度下彈道-擴散導熱的拓撲優(yōu)化方法.用聲子BTE代替傅里葉導熱定律作為導熱本構(gòu)方程,以反映聲子彈道輸運; 用SIMP方法對聲子平均自由程的倒數(shù)進行插值以引入拓撲優(yōu)化,并通過對相對密度變量全局梯度施加顯示約束的方法來確保解的適定性和網(wǎng)格無關(guān)性.將此方法用于體點導熱優(yōu)化,發(fā)現(xiàn)彈道-擴散導熱機制下,拓撲優(yōu)化預(yù)測的填充材料分布明顯不同于傅里葉定律下經(jīng)典的樹狀分布,且會隨 Kn 的變化而變化,這與SIMP方法的插值方式和聲子彈道輸運有關(guān).

2 優(yōu)化方法

2.1 聲子玻爾茲曼輸運方程

聲子的運動規(guī)律符合聲子BTE,穩(wěn)態(tài)、弛豫時間近似下的聲子BTE寫作

其中f是聲子分布函數(shù); vg是群速度矢量; f0是平衡態(tài)聲子分布函數(shù),即玻色-愛因斯坦分布; τ 是弛豫時間.聯(lián)合求解(1)式和能量守恒方程,就可以獲得聲子的微觀運動規(guī)律和分布函數(shù),進一步可計算得到溫度、熱導率等宏觀熱性質(zhì).

本文選擇離散坐標法(discrete ordinate method,DOM)[39-41]來求解聲子BTE,此方法最早用于輻射輸運方程的求解,相關(guān)計算方法的發(fā)展較為成熟.根據(jù)聲子和光子的類比,可定義聲子輻射強度其中ω表示聲子頻率,D(ω)表示態(tài)密度,“”表示對聲子極化的求和.利用此定義,可以將聲子BTE轉(zhuǎn)化為聲子的輻射輸運方程(equation of phonon radiative transfer,EPRT)[42]

(3)式中包含N+1個等式,與未知的聲子輻射強度數(shù)量相對應(yīng),方程組封閉.

得到聲子輻射強度后便可計算溫度.聲子發(fā)生彈道輸運時,當?shù)責崃W條件遠離平衡態(tài),傳統(tǒng)的在熱平衡態(tài)下定義的溫度失去意義[43],而如何在這樣的非平衡態(tài)下定義溫度一直存在爭議[44].本文采用等效平衡溫度的概念[45,46],認為空間內(nèi)某個離散單元在發(fā)射能量時,產(chǎn)生的聲子服從玻色-愛因斯坦分布,根據(jù)發(fā)射能量反解出的分布函數(shù)中的溫度,便是等效平衡溫度.使用DOM離散輻射強度后,等效平衡溫度的計算式為其中 σphonon是聲子斯特藩-玻爾茲曼參數(shù)[47].使用這種等效平衡溫度的好處是,在接近擴散導熱時,解聲子BTE所得結(jié)果能夠與傅里葉導熱定律的結(jié)果相吻合[39].

2.2 拓撲優(yōu)化

拓撲優(yōu)化本質(zhì)上屬于大規(guī)模0-1整數(shù)規(guī)劃問題,區(qū)域內(nèi)的每個點要么是空洞,要么是實體.這樣的離散優(yōu)化問題在數(shù)學上是病態(tài)的,解的存在性、收斂性以及求解算法的實現(xiàn)都存在著巨大的困難[48].實現(xiàn)拓撲優(yōu)化的關(guān)鍵是對優(yōu)化問題進行適當?shù)恼齽t化處理,將設(shè)計變量由離散的變?yōu)檫B續(xù)的,從而能夠直接使用針對連續(xù)設(shè)計變量優(yōu)化問題的求解方法.在這方面,SIMP方法[49]是經(jīng)典而常用的技巧,其核心思想是將材料性質(zhì)設(shè)定為關(guān)于相對密度的冪函數(shù),從而將設(shè)計變量轉(zhuǎn)變?yōu)檫B續(xù)的相對密度.考慮到(2)式作為描述非傅里葉導熱的本構(gòu)方程,僅有聲子消光系數(shù) κ 是與材料類型有關(guān)的性質(zhì)[38],于是本文提出對 κ 進行插值的SIMP方法來實現(xiàn)拓撲優(yōu)化.空間中任意位置的消光系數(shù)寫作

其中 ρ 是值在[0,1]區(qū)間上連續(xù)分布的相對密度變量; κs和 κf分別是基底材料和填充材料的消光系數(shù),亦即各自聲子平均自由程的倒數(shù); p是值大于1的懲罰因子.懲罰因子的作用是在優(yōu)化過程中對介于0和1之間的中間密度值進行懲罰,使得 ρ 的值逐漸趨向于0和1兩個端點值,從而在優(yōu)化結(jié)果中得到清晰的空洞(ρ=0)和實體(ρ=1)分布.關(guān)于SIMP方法插值方式的選擇,將在下一節(jié)結(jié)合優(yōu)化結(jié)果做進一步討論.

即使采用SIMP方法,拓撲優(yōu)化仍會面臨解的不存在性、不收斂性和不唯一性等問題,導致優(yōu)化結(jié)果會出現(xiàn)棋盤格、網(wǎng)格依賴性等問題[50].為了得到可靠的優(yōu)化結(jié)果,需要引入更多的正則化處理以抑制數(shù)值不穩(wěn)定性.本文采用對相對密度變量 ρ 的全局梯度施加顯式約束的方法[15],即在目標函數(shù)中增加關(guān)于 ρ 的全局梯度的權(quán)重懲罰項,以約束ρ的變化程度.至此,可寫出通過求解EPRT來獲得溫度分布的二維體點導熱拓撲優(yōu)化的數(shù)學模型為：

其中 Tave是系統(tǒng)平均溫度,Tave,ref是人為選取的參考溫度,α 是相對密度變量全局梯度的權(quán)重系數(shù),是無量綱的相對密度變量梯度,A=a2是二維區(qū)域的面積,? 是預(yù)先給定的高導熱填充材料占整個區(qū)域的比例.對平均溫度和設(shè)計變量梯度進行無量綱化的目的是使二者通過線性求和構(gòu)成實際的目標函數(shù)g.現(xiàn)在的模型以降低系統(tǒng)平均溫度為優(yōu)化目標,如果目標是減小系統(tǒng)的最高溫度,則可以將優(yōu)化對象變更為文獻[51]中提到的幾何平均溫度,其他設(shè)置不變.需要說明的是,因為本文的關(guān)注點是如何在導熱優(yōu)化中考慮聲子彈道輸運并分析彈道輸運對優(yōu)化結(jié)果的影響,所以優(yōu)化模型中并未考慮不同材料之間的界面熱阻,這樣可以更直觀地反映彈道效應(yīng)的影響.界面熱阻對現(xiàn)有優(yōu)化結(jié)果的影響將在第3節(jié)中進行討論.

求解(5)式的拓撲優(yōu)化的流程圖如圖2所示,包括以下步驟：1)初始化,對區(qū)域進行離散并任意給定一組 ρ 的初始值; 2)系統(tǒng)重分析,根據(jù)SIMP方法插值得到當前設(shè)計點對應(yīng)的消光系數(shù),再求解EPRT和穩(wěn)態(tài)能量守恒方程獲得聲子輻射強度,進而計算平均溫度和目標函數(shù); 3)靈敏度分析,計算目標函數(shù)和約束函數(shù)對設(shè)計變量的導數(shù); 4)優(yōu)化求解,根據(jù)導數(shù)信息,確定滿足約束且減小目標函數(shù)的設(shè)計變量演化方向(可行下降方向),并計算最佳步長,更新當前設(shè)計變量值; 5)收斂判斷,滿足收斂條件時結(jié)束迭代,輸出結(jié)果; 否則重復步驟2)-4).

圖2 體點導熱拓撲優(yōu)化的流程圖Fig.2.The flow chart of topology optimization for the VP problem.

3 結(jié)果與分析

在體點導熱問題中,設(shè)置低溫邊界與區(qū)域邊長的比值為 c/a=0.1 ,低溫邊界的溫度為 T0=300K ,填充材料體積占比為 ?=0.1 ,設(shè)計變量初始值取為 ρ=? ,即材料均勻分布.定義無量綱空間坐標為 X=x/a ,Y=y/a ,參考文獻[12]定義無量綱溫度其中 ks是基底材料熱導率,這樣無量綱化的好處是消除了內(nèi)熱源功率、熱源面積、熱導率和邊界溫度具體值對優(yōu)化結(jié)果的影響,使填充材料和基底材料間導熱性能差異成為優(yōu)化結(jié)果的決定因素.在數(shù)值計算中,用數(shù)量為n×n的方形網(wǎng)格離散計算域,網(wǎng)格尺寸為 h=a/n.DOM的離散方向數(shù)設(shè)定為 N=24 ,此值既能避免求解溫度場時出現(xiàn)過強的“射線”效應(yīng)[41],又不至于讓計算量過大.參考文獻[52],相對密度變量梯度的無量綱化方式為方法的懲罰系數(shù)取經(jīng)典值 p=3.為了提高求解效率,對聲子BTE使用灰體近似,迭代過程中的靈敏度分析使用伴隨方法[53],優(yōu)化求解則采用構(gòu)造原問題局部近似模型的移動漸進方法(method of mo ving asymptotes,MMA)[54].收斂準則設(shè)定為其中上標j表示迭代次數(shù),認為滿足此條件時設(shè)計變量、目標函數(shù)值都達到穩(wěn)定.考慮到大部分區(qū)域為基底材料,用基底材料的平均自由程來計算努森數(shù),即 Kn=ls/a.

3.1 擴散導熱的體點優(yōu)化

在進行求解聲子BTE的拓撲優(yōu)化之前,先在傅里葉定律適用的條件下實現(xiàn)體點導熱的拓撲優(yōu)化,以檢驗所發(fā)展的拓撲優(yōu)化方法.盡管凡是在實體和空材料間構(gòu)造連續(xù)函數(shù)并對中間密度值進行懲罰的方法都可稱為SIMP方法,但通常都是對本構(gòu)方程的系數(shù)進行插值[15].對基于傅里葉定律的拓撲優(yōu)化,本構(gòu)方程中與材料性質(zhì)有關(guān)的系數(shù)是熱導率,對應(yīng)SIMP方法的插值公式為k(ρ)=ks+(kf-ks)ρp,其中 kf是填充材料熱導率.此時,填充材料和基底材料導熱能力的差異通過熱導率比kf/ks來體現(xiàn),kf/ks越大,傳導同樣熱量所需要的高導熱材料橫截面積越小,因而在相同的填充量下填充材料能夠更遠地延伸到區(qū)域內(nèi)部,形成的樹狀結(jié)構(gòu)的“樹枝”更為細長,優(yōu)化效果也更好[12].對本文提出的優(yōu)化模型,本構(gòu)方程中與材料性質(zhì)有關(guān)的系數(shù)是聲子消光系數(shù) κ ,對應(yīng)的插值公式為(4)式.根據(jù)動力學理論,熱導率和聲子平均自由程間的關(guān)系為其中C是材料比熱,v是聲子群g速度的大小.材料的比熱、群速度通常是確定的,于是有 kf/ks∝lf/ls,意味著不同材料導熱性能的差別主要源自聲子平均自由程的差別[55].此時(4)式中關(guān)于聲子消光系數(shù) κ (=l-1)的插值,對應(yīng)于傅里葉導熱定律下對熱導率倒數(shù) k-1的插值,而非傳統(tǒng)的對熱導率k的插值.在 kf/ks=500 ,α=0.05 ,n=80的條件下,分別對k和 k-1進行插值,基于傅里葉導熱定律的拓撲優(yōu)化得到的結(jié)果如圖3所示.設(shè)定 kf/ks=500 是為了讓填充材料和基底材料的導熱能力存在較大差異,更直觀地體現(xiàn)優(yōu)化的效果,且實際電子器件中絕緣材料和導電材料的熱導率通常也相差2個數(shù)量級[1].文獻中關(guān)于傅里葉導熱定律下體點導熱拓撲優(yōu)化的研究也大多采用此熱導率比值.圖3左側(cè)表示 ρ 在區(qū)域內(nèi)的分布情況,白色代表基體材料(ρ=0),黑色表示填充材料(ρ=1),顏色越深意味著該位置的 ρ 值越接近于1;右側(cè)則是無量綱溫度 T?的分布.插值k的拓撲優(yōu)化得到了充分伸入內(nèi)部區(qū)域、含多個分岔和枝的樹狀填充材料分布,如圖3(a)所示,符合文獻[10-13]中的結(jié)果,對應(yīng)的無量綱溫度平均值是插值 k-1的優(yōu)化結(jié)果則由圖3(b)給出,填充材料集中在區(qū)域下半平面,與低溫邊界相接觸,向區(qū)域內(nèi)伸展的過程中只形成了2處分岔,此時與圖3(a)相比,圖3(b)中填充材料的主干更短、更寬,枝數(shù)量更少,對應(yīng)的平均溫度和最大溫度更高.在擴散導熱條件下,不同形狀、尺寸的填充材料熱導率相等,具有相同的構(gòu)建高導熱通道的能力.為了降低溫度,填充材料要設(shè)法伸入內(nèi)部區(qū)域,縮短高導熱通道和區(qū)域內(nèi)熱源間的距離,所以圖3(a)的結(jié)果要比圖3(b)更好.受制于SIMP方法的插值方式,圖3(b)得到的是體點導熱優(yōu)化的某個局部最優(yōu)解,但填充材料仍然有向區(qū)域內(nèi)伸展的趨勢.圖3的結(jié)果表明,SIMP方法的插值方式會影響拓撲優(yōu)化的效果,對本構(gòu)方程的系數(shù)進行插值是更合理的選擇.以傅里葉導熱定律為例,盡管插值 k-1的拓撲優(yōu)化也能收斂,但對應(yīng)的是優(yōu)化問題的某個局部最優(yōu)解,優(yōu)化效果不如插值k的結(jié)果好.類似地,當導熱的本構(gòu)方程變成(2)式時,SIMP方法應(yīng)當對方程的系數(shù) κ 進行插值.

圖3 擴散導熱下不同插值方式對材料分布和溫度分布的影響 (a)插值k; (b)插值k-1Fig.3.The effect of interpolation ways on the material distributions and temperature distributions in diffusive heat conduction：(a) Interpolating k; (b) interpolating k-1.

圖4 擴散導熱下拓撲優(yōu)化得到的材料分布隨 α 和n的變化Fig.4.The material distributions obtained by topology optimization in diffusive heat conduction varying with α and n.

除了SIMP方法的插值方式,α 的取值也很重要,因為 α 顯著響著拓撲優(yōu)化解的數(shù)值穩(wěn)定性[15].保持 kf/ks=500 不變,并使用對k的插值,擴散導熱條件下拓撲優(yōu)化得到的填充材料分布隨 α 和n的變化如圖4所示.整體來看,填充材料的形狀都呈樹狀,與圖3(a)的結(jié)果大致相同,但細節(jié)上存在不少差異.當 α=0 時,隨著網(wǎng)格密度的增大,高導熱通道的主干周圍會出現(xiàn)越來越多狹長的末梢,這是因為拓撲優(yōu)化的數(shù)值解天然地存在著網(wǎng)格依賴性.當 α=0.05 時,由于引入了對設(shè)計變量全局梯度的約束,拓撲優(yōu)化結(jié)果的網(wǎng)格依賴性得到有效抑制,網(wǎng)格細化并未導致太多狹長末梢.進一步增大到 α=0.1 ,不同網(wǎng)格密度下拓撲優(yōu)化結(jié)果間的差異更小,網(wǎng)格無關(guān)性更好.但是,比較相同網(wǎng)格密度的拓撲優(yōu)化結(jié)果發(fā)現(xiàn),增大 α 的值后,中間密度值對應(yīng)的灰色區(qū)域增多,填充材料和基底材料之間的邊界變得模糊,這是控制設(shè)計變量全局梯度的代價.綜合考慮后,本文的拓撲優(yōu)化選擇 α=0.05 ,n=80,此時所得結(jié)果具有較好的網(wǎng)格無關(guān)性和較清晰的材料邊界,且計算量相對較少.

3.2 彈道-擴散導熱的體點優(yōu)化

根據(jù)上文確定的參數(shù)設(shè)置,對(5)式進行數(shù)值迭代,獲得了體點導熱在彈道-擴散導熱機制下的拓撲優(yōu)化解.設(shè)定 lf/ls=500 以便和傅里葉導熱定律下的結(jié)果對比,不同 Kn 下優(yōu)化得到的材料分布和溫度分布如圖5所示.總體來看,填充材料集中在區(qū)域下半平面,與低溫邊界充分接觸,其形狀完全不同于圖3(a)所示的伸入內(nèi)部區(qū)域且含多個枝的樹狀結(jié)構(gòu),而與圖3(b)的材料分布更為接近.如3.1節(jié)所述,本文采用的插值 κ 的方法,對應(yīng)于傅里葉定律下對 k-1的插值,因而優(yōu)化結(jié)果與插值k-1的結(jié)果相似.可以預(yù)見的是,在 Kn→0 時,導熱過程接近純擴散導熱,優(yōu)化會趨向于傅里葉定律下插值 k-1的結(jié)果.但受聲子彈道輸運的影響,圖5與圖3(b)仍然存在明顯的差異.當 Kn=0.002 時,低溫邊界處幾乎無溫度跳躍,區(qū)域內(nèi)無量綱溫度的最大值略高于圖3(b)中的0.23,表明大部分聲子以擴散方式導熱,但仍有部分聲子發(fā)生了彈道輸運,這是因為填充材料的聲子平均自由程已和區(qū)域尺寸相當(lf～a).此時填充材料在低溫邊界處連通,主干呈“U”型,兩側(cè)共有6處外凸的細枝.增大 Kn 至0.01,靠近低溫邊界處 T?=0.06 ,發(fā)生了一定的溫度跳躍,則是0.44,約為圖3(b)結(jié)果的2倍,表明彈道效應(yīng)增強.此時填充材料分布的“U”形主干開口變大、寬度變厚,枝數(shù)量減少為2個.當 Kn=0.1 時,填充材料的體聲子平均自由程已比區(qū)域尺寸大一個數(shù)量級,邊界處無量綱溫度的跳躍值接近0.6,區(qū)域內(nèi)峰值溫度幾乎是圖3(b)結(jié)果的6倍,彈道效應(yīng)已十分顯著.此時,填充材料聚集在低溫邊界附近,除了頂部有凹陷外,其它方向的伸展長度基本一致,主干和枝合并形成團狀分布.

為了檢驗所提出的拓撲優(yōu)化方法的效果,分別在不同努森數(shù)下,對圖3和圖5中出現(xiàn)的5種材料分布進行熱仿真,得到的值列于表1,括號內(nèi)的數(shù)據(jù)表示的值,即相對于材料均勻分布時的變化程度.表1顯示,不同努森數(shù)下,本文發(fā)展的拓撲優(yōu)化方法所對應(yīng)的值最低,如表中加粗數(shù)據(jù)所示,表明所預(yù)測的材料分布方式確實是最優(yōu)的,驗證了聲子BTE和拓撲優(yōu)化相結(jié)合的優(yōu)化方法的正確性.以 Kn=0.1 時為例,解傅里葉定律的拓撲優(yōu)化所得的材料分布,對應(yīng)的值大概為80%,而解聲子BTE的拓撲優(yōu)化得到的材料分布,對應(yīng)的值在70%左右,進一步下降了10個百分點,其中又以Kn=0.1條件下的拓撲優(yōu)化所得的構(gòu)型效果最好,此外,對相同的材料分布方式,隨著 Kn 的增大,的值增大,表明添加高導熱材料所能帶來的優(yōu)化效果變差.例如對圖3(a)所示的樹狀結(jié)構(gòu)(擴散導熱條件下插值k的拓撲優(yōu)化結(jié)果),Kn=0.002 時,是33.7%,但Kn=0.01時此值上升到53.2%,Kn=0.1 時只有86.0%.圖5和表1的結(jié)果共同表明,彈道-擴散導熱機制下,體點導熱問題的最佳材料分布會隨著彈道效應(yīng)強度的變化而變化,基于傅里葉定律的拓撲優(yōu)化所得的材料分布不再是最優(yōu)的.因此,對微納米尺度下的散熱優(yōu)化問題,在設(shè)計階段必須要考慮聲子彈道輸運的影響.值得一提的是,如果將現(xiàn)有優(yōu)化模型中對 κ 的插值改為對其倒數(shù)l的插值,優(yōu)化模型的收斂性會嚴重惡化,材料分布會出現(xiàn)嚴重的棋盤格問題和大量的中間密度值.

圖5 彈道-擴散導熱下拓撲優(yōu)化的材料分布和溫度云圖隨 Kn 的變化 (a) Kn=0.002; (b) Kn=0.01; (c) Kn=0.1Fig.5.The topology optimization obtained material distributions and temperature distributions varying with Kn in ballistic-diffusive heat conduction：(a) Kn=0.002; (b) Kn=0.01; (c) Kn=0.1.

進一步分析圖5中拓撲優(yōu)化得到的最優(yōu)材料分布,發(fā)現(xiàn)隨著 Kn 的增大,彈道效應(yīng)增強,填充材料從低溫邊界向內(nèi)部區(qū)域發(fā)展的趨勢減弱,且分布形狀中微小的枝所占的比例減少,粗壯的主干所占的比例增多,這可以由彈道輸運產(chǎn)生熱導率的尺寸效應(yīng)來解釋.不同于擴散導熱,彈道-擴散導熱機制下,結(jié)構(gòu)的等效熱導率會隨尺寸的減小而降低.因為微小枝的尺寸比粗壯主干的尺寸更小,所以前者的等效熱導率比后者更低,前者在構(gòu)建高導熱通道方面的性能不如后者; 且隨著彈道效應(yīng)的增強,這種能力上的差別會增大.在 Kn 較小時,不同尺寸的填充材料間導熱性能差異較小,影響優(yōu)化效果的主要還是熱源和高導熱通道的距離,所以圖5(a)的結(jié)果像圖3(b)那樣有一些分岔,且伸入內(nèi)部區(qū)域較多.隨著 Kn 的增大,彈道效應(yīng)增強,結(jié)構(gòu)等效熱導率的降低對優(yōu)化效果的影響加強,拓撲優(yōu)化傾向于產(chǎn)生等效熱導率相對較大的粗壯結(jié)構(gòu),而非伸入內(nèi)部區(qū)域的微小分散結(jié)構(gòu).Kn=0.1 時,盡管基底材料聲子平均自由程僅為系統(tǒng)特征尺寸的0.1倍,但填充材料的聲子平均自由程已經(jīng)是系統(tǒng)特征尺寸的50倍,整個系統(tǒng)的彈道效應(yīng)已十分顯著,團狀的填充材料分布在構(gòu)建高導熱通道方面的性能最優(yōu),但由于填充材料比例限制,通道只能覆蓋低溫邊界附近的區(qū)域.彈道導熱極限下,拓撲優(yōu)化的結(jié)果中填充材料聚集在低溫邊界附近,沿各個方向的伸展長度基本相同,呈半圓形團狀分布.圖5(c)的結(jié)果已經(jīng)較為接近彈道極限下的結(jié)果.溫度分布中低溫區(qū)域的輪廓反映了填充材料所形成的高導熱通道的效果.在圖3和圖5中,低溫區(qū)域的輪廓都與填充材料的形狀大體一致,表明所有填充材料都有效地形成了高導熱通道.根據(jù)上述分析,在彈道效應(yīng)較強時,如果填充材料的形狀中含有較多的微小結(jié)構(gòu),由于這些小尺寸結(jié)構(gòu)在構(gòu)建高導熱通道方面是低效的,對應(yīng)區(qū)域的溫度將無法得到有效降低,溫度輪廓將勢必不同于填充材料的形狀.作為檢驗,給出擴散導熱條件下插值 k-1的拓撲優(yōu)化所預(yù)測的材料分布(圖3(b))在 Kn=0.1 時的溫度場如圖6所示.盡管材料分布與圖3(b)相同,但圖6的溫度輪廓明顯不同于圖3(b),反而與圖5(c)的溫度輪廓更為接近,說明填充材料中粗壯的主干形成了高效的導熱通道,但細小的枝因為等效熱導率更低,所以在導熱方面是相對低效的,未能有效降低所在區(qū)域的溫度.

表1 不同材料分布在不同努森數(shù)下對應(yīng)的無量綱溫度平均值Table 1.The average dimensionless temperature of different material distributions at different Knudsen numbers.

對體材料熱導率比值和界面熱阻對拓撲優(yōu)化結(jié)果的影響做簡單說明.盡管上述分析建立在kf/ks=500的條件下,但因為彈道輸運主要通過熱導率的尺寸效應(yīng)來影響拓撲優(yōu)化結(jié)果,所以改變kf/ks的值不影響體點導熱拓撲優(yōu)化結(jié)果隨彈道效應(yīng)強度的定性變化.定量上,體材料熱導率比的值越大,發(fā)生彈道輸運時不同尺寸結(jié)構(gòu)等效熱導率的差異越大,相同努森數(shù)下彈道輸運對拓撲優(yōu)化結(jié)果的影響越顯著.界面熱阻會使得材料間界面處出現(xiàn)溫度不連續(xù),導致系統(tǒng)總熱阻增大,平均溫度升高[4].雙層材料間界面熱阻的影響會隨著材料尺寸的減小而增大[56],如果在拓撲優(yōu)化中考慮界面熱阻,為了降低界面熱阻對導熱的阻礙作用,優(yōu)化結(jié)果會朝著增大結(jié)構(gòu)尺寸的方向發(fā)展.圖5中隨著彈道效應(yīng)的增強,填充材料構(gòu)型也明顯表現(xiàn)出粗壯結(jié)構(gòu)增多、微小結(jié)構(gòu)減少的趨勢,這表明雖然現(xiàn)在的優(yōu)化模型未考慮界面熱阻,但所得的結(jié)果有利于抑制界面熱阻的作用.

圖6 擴散導熱下插值 k-1 的拓撲優(yōu)化所得的材料分布在 Kn=0.1 時的溫度分布Fig.6.The temperature distribution for the material distribution obtained by topology optimization which interpolates k-1in diffusive heat conduction at Kn=0.1.

4 結(jié) 論

針對微納米電子器件的散熱優(yōu)化設(shè)計,結(jié)合聲子玻爾茲曼方程和固體各向同性材料懲罰方法,發(fā)展了適用于微納尺度下彈道-擴散導熱的拓撲優(yōu)化方法.聲子玻爾茲曼方程被轉(zhuǎn)化為聲子輻射輸運方程,然后用離散坐標法求解以獲得某個設(shè)計方案下的溫度場; 拓撲優(yōu)化則通過對聲子平均自由程的倒數(shù)進行插值的固體各向同性材料懲罰方法來實現(xiàn),并采用對相對密度變量全局梯度施加顯示約束的方法來確保解的適定性和網(wǎng)格無關(guān)性.

成功實現(xiàn)了體點導熱問題在彈道-擴散導熱機制下的拓撲優(yōu)化解,發(fā)現(xiàn)彈道-擴散導熱機制下,拓撲優(yōu)化預(yù)測的最佳材料分布明顯不同于傅里葉定律下經(jīng)典的樹狀結(jié)構(gòu),且會隨努森數(shù)Kn的變化而變化,表明微納米尺度下的導熱優(yōu)化問題必須要考慮聲子彈道輸運的影響.在Kn較小時,解聲子玻爾茲曼方程的拓撲優(yōu)化的結(jié)果會趨向傅里葉定律下于對熱導率倒數(shù)進行插值的拓撲優(yōu)化結(jié)果,而非傳統(tǒng)的插值熱導率所對應(yīng)的樹狀結(jié)構(gòu).隨著Kn的增大,拓撲優(yōu)化結(jié)果中填充材料逐漸向低溫邊界附近聚集,其構(gòu)型中粗壯的主干結(jié)構(gòu)所占比例增大,微小的枝結(jié)構(gòu)所占的比例降低,且優(yōu)化后系統(tǒng)平均溫度和優(yōu)化前系統(tǒng)平均溫度的比值增大.

拓撲優(yōu)化結(jié)果和Kn的相關(guān)性與聲子彈道輸運產(chǎn)生的熱導率尺寸效應(yīng)有關(guān).擴散導熱機制下,不同尺寸的填充材料熱導率相同,均能有效構(gòu)建高導熱通道; 但隨著彈道效應(yīng)的增強,尺寸效應(yīng)使得材料分布中微小結(jié)構(gòu)的等熱導率比粗壯結(jié)構(gòu)低,微小結(jié)構(gòu)在提升系統(tǒng)導熱性能方面的作用不如粗壯結(jié)構(gòu),因而拓撲優(yōu)化會向增多粗壯結(jié)構(gòu)、減少微小結(jié)構(gòu)的方向發(fā)展.當彈道效應(yīng)足夠強時,填充材料會聚集在低溫邊界附近,主干和枝合并,呈團狀分布.