基于密度法的傳熱結構拓撲優化設計

崔天福， 丁曉紅， 侯麗園

（上海理工大學 機械工程學院，上海 200093）

具有散熱通道的傳熱結構的傳熱性能與散熱通道的材料、形態分布和外部熱環境密切相關.當材料和外部熱環境一定時，傳熱結構中的散熱通道就對整個結構的熱量傳遞起到主要作用，因此，如何在傳熱結構中合理布置散熱通道對提高結構的傳熱性能、延長其使用壽命至關重要.在發熱量比較集中、空間有限的集成電路中，可采用由高導熱材料構成的散熱通道直接嵌入到元器件內部，將熱量導出到外部環境中，再采用其它高效冷卻方式進行冷卻，可以有效解決空間限制問題和冷卻效率問題.傳統的傳熱結構設計一般基于傳熱學基本理論及實踐經驗，然而這種方法很難解決復雜傳熱結構的設計問題.針對這種復雜模型，可使用結構拓撲優化方法，計算出高自由度的散熱通道拓撲形態，在此形態的基礎上，通過經驗改進以及進一步的尺寸優化獲得最佳的結構尺寸.

結構拓撲優化以材料分布為優化對象，在設計空間內尋求最佳的分布方案，以得到結構某種性能最優.在拓撲優化方法中，理論相對成熟且有一定工程應用的方 法 主 要 有 密 度 法［1］、均 勻 化 方 法［2－3］及Level Set方法［4－5］等.近年來這些方法在傳熱結構拓撲優化設計中已有所應用.如Bends?e等［6］將結構拓撲優化中的密度法直接拓展到了傳熱結構優化當中，并對簡單的穩態熱傳導問題進行了研究.Iga［7］在均勻化理論的基礎上建立了以結構總勢能為目標函數的優化數學模型，確定了結構具有最佳傳熱效應時，在Dirichlet邊界、Neumann邊界、對流邊界情況下目標 函 數 的 最 大 值、最 小 值 問 題.Yamada［8－9］以Lveleset方法為基礎建立了熱擴散最大為目標的傳熱拓撲優化數學模型，并對熱力耦合情況下考慮熱擴散效果的結構剛度最大化問題進行了研究.對于體－點熱 傳 導 問 題，Bejan［10］提 出 構 形 理 論（constructal theory）進行求解，該方法首先對初始單元進行優化，獲得高導熱通道的尺寸參數，然后對由初始單元組裝而成的一級裝配體進行優化，獲得高導熱通道的尺寸參數，逐級裝配直至將高導熱材料組成的散熱通道覆蓋至所需要的散熱面積.Gersborg－Hansen［11］提出了在有限體積法基礎上的傳熱結構拓撲優化的方法.丁曉紅［12－13］通過研究植物根系的成長機理，提出自適應成長法來解決薄板結構的加強筋分布的設計問題，并將自適應成長法應用于體－點散熱問題，該方法通過對局部分枝的靈敏度分析，使靈敏度大的區域分枝生長迅速，靈敏度小的區域分枝生長緩慢或者退化，從而在迭代過程中實現分枝結構的自適應功能，使得最終分枝的拓撲結構具有最優的形態布局.

本文以變密度法為理論基礎，并以SIMP（solid isortopic material with penalization model）插值模型［14－15］為材料插值模型，建立結構散熱弱度為目標函數的優化數學模型，推導出滿足KKT（Karush－Kuhn－Tucker）條件的傳熱結構拓撲優化迭代公式.采用高階單元方法來解決優化過程中存在的數值不穩定現象，然后通過不同邊界組合下的數值算例，研究了拓撲形態分布上的一些特性.

1 理論分析與優化模型

1.1 傳熱結構的拓撲優化數學模型

如圖1所示，在設計區域Ωd內有生熱率載荷Q，存 在Dirichlet 邊 界ΓT，邊 界 溫 度 為T0.在Neumann邊界Γq，熱流密度為q.傳熱結構的拓撲優化問題為：如何在區域Ωd中布置高導熱材料，形成散熱通道，將區域內的熱量傳送到邊界.

定義 TTP 為結構的散熱弱度［16］，其中，T 和P是離散結構中的節點溫度矢量和等效節點載荷矢量，等效節點載荷矢量包括生熱率載荷和熱流邊界引起的等效載荷.結構的散熱弱度是結構對外部熱載荷的響應，由于散熱弱度與熱勢能具有一致性，當散熱弱度最小時，在相同的外部熱載荷作用下所能傳遞的熱量最大，使結構具有最佳的溫度分布［7］.最佳的溫度分布體現在結構的平均溫度梯度小、平均溫度低等方面.因此，可以通過使結構的散熱弱度最低來確定熱傳導問題中高導熱材料的最優分布.

圖1 二維熱傳導問題Fig.1 Two－dimensional heat conduction problem

在基于密度法的結構拓撲優化設計中，SIMP插值模型是廣泛采用的模型，如圖2所示（見下頁）.在該插值模型中引入了一個偽密度變量的概念，每一個離散單元都對應一個偽密度變量.對于傳熱結構拓撲優化問題，可使偽密度變量關聯的材料屬性為離散單元的熱傳導系數.圖2的SIMP 插值關系可以表示為

式中，xi為離散單元偽密度變量，即優化設計變量；λi表示偽密度變量xi的單元熱傳導系數；λmax表示傳熱結構中高導熱材料的熱傳導系數；λmin表示低導熱材料的熱傳導系數；p 為懲罰系數.

引入p 的目的是在優化過程中避免中間密度單元的形成，使優化后的傳熱結構中盡可能只存在λmax和λmin兩種導熱系數的材料.當p＜1時，對中間密度具有強化作用；當p＞1時，對中間密度具有懲罰作用；當p＝1時，不對中間密度進行懲罰或者強化.通常p 取2～5范圍內的值.

圖2 SIMP插值曲線Fig.2 SIMP interpolation curve

為了便于表達數學模型，引入一個高、低導熱材料的熱傳導系數比h＝λmax／λmin.式（1）可以表示為

式中，mi為單元導熱比例系數，表示偽密度變量xi的導熱系數λi與低導熱材料導熱系數λmin之間的比值.

由單元熱傳導系數與單元熱傳導矩陣的線性關系可推出

式中，kmin是當熱傳導系數為λmin時的單元熱傳導矩陣；ki是導熱系數為λi時的單元熱傳導矩陣.

以結構的散熱弱度C＝TTP 為目標函數，在SIMP插值模型下，其拓撲優化模型可表示為

式中，ti為第i個單元的節點溫度矢量；vi為第i個單元的體積；K 為優化模型的總熱傳導矩陣；V＊為體積約束上限值；V 為優化后的體積；V0為設計域體積；f 為體積約束百分率；xmin為偽密度下限值；n為離散單元總數.

1.2 基于最優準則法的迭代公式

由上述目標函數、體積約束、平衡方程約束、設計變量約束構造拉格朗日函數，當該構造函數滿足KKT 條件時，目標函數C 具有最小值.綜合考慮優化設計變量的邊界約束以及熱傳導矩陣的對稱性，對于離散單元應滿足KKT 條件的表達式為

式中，λl為Lagrangian體積乘子.

由（6）式可得

其中

式中，w 為當前迭代次數；si為第i個單元的熱耗散（heat dissipation）；為第w 次迭代過程中第i個單元的跌倒系數.

將式（8）和式（9）代入式（7），整理得

將式（10）作為優化設計準則，得出基于該準則下的迭代公式.

式中，η是為了保證數值計算的穩定性而引入的一個阻尼系數；z為密度變化上限值.

基于最優準則法實現的熱傳導體結構拓撲優化程序流程如圖3所示.首先建立有限元模型，包括幾何模型的建立以及材料、熱載荷、邊界約束的定義.然后初始化設計變量，將材料參數用設計變量來表示.計算優化模型的目標值以及敏度值.計算滿足體積約束條件的Lagrangian體積乘子，并通過當前迭代步中的體積乘子和設計變量來更新下一個迭代步中的設計變量.當設計變量滿足收斂條件時，優化過程結束.否則，進入下一個循環過程，直到滿足收斂條件.對于迭代過程中的收斂性判定，采用相鄰迭代的最大密度變化值為判定依據，表達式為

式中，Xw為第w 次迭代的偽密度矢量；ε是很小的密度參考值，一般取0.01.

圖3 優化程序流程圖Fig.3 Flow of topology optimization

與以目標變化值為收斂判定條件時相比，采用密度變化值判定條件時，所得到的拓撲形態密度分布更為集中，中間密度單元也相對較少，但迭代次數也相應增多.

1.3 數值不穩定處理方法

在拓撲優化過程中，由于場函數不穩定或數值奇異解的存在導致優化后易出現數值不穩定現象，如棋盤格現象、網格依賴性及局部極值等［17］.棋盤格和網格依賴這兩種現象一般同時出現在優化結果中，而能夠有效去除棋盤格的方法通常也能有效克服網格依賴現象.棋盤格的產生與分析單元的選擇有關.研究表明［18－19］，合理選擇高階單元或采用非協調元，可有效降低或消除棋盤格.如圖4 所示，8節點的高階單元比4節點的低階單元在每個邊中心點處多出了1個中間節點，使節點溫度函數的一階導數在單元節點更加圓滑，可有效降低棋盤格現象.

圖4 8節點單元Fig.4 8node element

2 數值算例

2.1 平面問題的Dirichlet邊界算例

首先，考慮如圖5所示的設計問題.設計空間是0.1m×0.1m 的正方形.在內部設計域Ωd存在均勻 生 熱 率Q1＝6×104W／m3，在 下 邊 界 存 在Dirichlet邊界ΓT，邊界溫度為T0＝0 ℃，邊界長度L1＝0.01m，低導熱系數λmin＝1 W／（m·K），熱傳導系數比h＝400，體積約束分別為10%，20%，30%，50%時的傳熱結構拓撲形態如圖6所示（見下頁）.設計域內黑色部分表示高導熱材料，白色部分表示低導熱材料.由拓撲形態可知，高導熱材料主要分布在底部溫度邊界附近，并以枝狀形態向四周分歧延伸.隨著高導熱材料體積約束值的提高，出現了枝的形態逐漸變粗、分歧細化、高導熱材料覆蓋區域擴大等現象.

圖5 均勻生熱率載荷Dirichlet邊界設計模型Fig.5 Uniform heat generation model satisfying Dirichlet boundary condition

圖6 不同體積約束下均勻生熱率載荷Dirichlet邊界的設計結果Fig.6 Topology forms of uniform heat generation models under different volume fraction

圖7 非均勻生熱率載荷Dirichlet邊界設計模型Fig.7 Non－uniform heat generation design model satisfying Dirichlet boundary condition

表1 生熱率載荷Tab.1 Uniform heat generation rate （W·m－3）

其次，考慮在Dirichlet邊界下，非均勻生熱率載荷的算例模型，如圖7所示（見下頁）.設計空間是由面積相等的4個子域Ω1，Ω2，Ω3，Ω4組成的正方形區域，每個子域中分別存在均勻生熱率載荷Q1，Q2，Q3，Q4，如表1所示（見下頁）.通過不同的生熱率載荷組合來實現不同的非均勻生熱率載荷工況.工況1和工況2生熱率載荷為上下分布，工況1的上半部分Ω1和Ω2是高熱載荷區域，下半部分Ω3和Ω4是低熱載荷區域，工況2與工況1的載荷分布相反.工況3生熱率載荷為對角分布，Ω2和Ω3是高載荷區域，Ω1和Ω4是低載荷區域.工況4生熱率載荷為對角分布，Ω1和Ω3是高載荷區域，Ω2和Ω4是低載荷區域.體積約束均為30%，低導熱系數λmin＝1 W／（m·K），熱傳導系數比均為400.在不同生熱率組合下的傳熱結構拓撲形態如圖8所示.在工況1中，高導熱材料主要集中在上半部的高載荷區域和高載荷區與溫度邊界的最短路徑上.在工況2中，高導熱材料主要集中在下半部高載荷區域，并且有一部分延伸到了上半部的低載荷區域.在工況3中，高導熱材料主要集中在對角形的高載荷區域，在載荷區域高導熱材料未被填充.在工況4 中，大部分高導熱材料集中在左側的高載荷區，只有一少部分分支結構延伸到了低載荷區.在4個工況中都明顯出現了高導熱材料往高熱載區域集中的趨勢.

圖8 相同體積約束下非均勻生熱率載荷Dirichlet邊界設計結果Fig.8 Topology forms of non－uniform heat generation models under same volume fraction

圖9（a）的算例模型是四邊Dirichlet邊界ΓT，并在內部存在均勻的生熱率載荷 Qd＝3×105W／m3.由圖9（b）拓撲形態可知，高導熱材料從四邊向中心區域擴散，并且具有四邊對稱分布.圖10（a）的算例模型是四角Dirichlet邊界ΓT，并在內部存在均勻的生熱率載荷Qd＝3×105W／m3.由圖10（b）拓撲形態可知，高導熱材料從四角向中心區域擴散，且在模型中出現了十字型的空白區域.圖11（a）的算例模型是四角Dirichlet邊界ΓT，并在中心點存在集中生熱率載荷Q1＝3×106W／m3.由圖11（b）拓撲形態可知，高導熱材料由四角和中心點形成了一個X 型的材料分布.圖12（a）的算例模型是四邊Dirichlet邊界ΓT，并在中心點存在集中生熱率載荷Q1＝3×106W／m3.由圖12（b）拓撲形態可知，高導熱材料由四角和中心點形成了一個十字型的材料分布.綜上所述，圖9（b）和圖10（b）均明顯出現了有分歧狀的拓撲結構，而圖11（b）和圖12（b）出現了清晰的X和十型的拓撲結構，這是由于集中生熱時，熱量流動不易受到分支結構的引導，所以，拓撲結構中高導熱材料都集中在生熱點與熱量出口的最短路徑上.

圖10 四角溫度邊界－均勻生熱Fig.10 Four corners temperature boundary with uniform heat generation

圖11 四角溫度邊界－中心生熱Fig.11 Four corners temperature boundary with centre point heat generation

2.2 三維問題的Dirichlet邊界算例

圖12 四邊溫度邊界－中心生熱Fig.12 Four sides temperature boundary with centre point heat generation

對于薄壁管狀體，由于管壁相對于管徑而言小很多，所以，此類問題可以采用殼單元進行計算.圖13是薄壁管狀體算例模型，在管的表面有尺寸0.04m×0.04 m 方 形Dirichlet 邊 界，邊 界 溫 度T0＝30 ℃，管體表面有生熱率載荷Q1＝10 W／m3，高導熱系數λmax＝80 W／（m·K），低導熱系數λmin＝0.2 W／（m·K），體積約束為20%.圖14（見下頁）為薄壁管狀體傳熱結構拓撲形態的不同視角顯示，圖14（a）是未旋轉時的軸測圖，圖14（b）是以管體中心為旋轉點繞X 軸順時針旋轉60°時的視角，圖14（c）是旋轉180°的視角，可以看出，高導熱材料并未出現在方形邊界內部，這是因為在邊界內部高導熱材料的有無并不能改變邊界上的溫度，所以，邊界上的目標敏度在優化過程中始終保持為零，從而高導熱材料并未填充到Dirichlet邊界內.在拓撲形態上，與平面算例的分布相似，但由于是管體結構，所以，在軸向上并沒有實現形態閉合.

圖13 三維薄壁管狀體結構模型Fig.13 Thin walled tube model

2.3 平面問題的混合邊界算例

圖14 薄壁管狀體結構拓撲結果（不同視角）Fig.14 Different views of the topology form

圖15 混合邊界模型Fig.15 Model satisfying both Dirichlet and neumann boundary conditions

熱傳導問題中存在Dirichlet邊界和Neumann邊界的混合邊界的算例模型如圖15所示（見下頁）.設計空間規格為0.1m×0.1m 的正方形，左下角有寬度為L1的流入熱量的Dirichlet邊界ΓT，右上角有寬度為L2的流出熱量的Neumann邊界Γq，4種工況參數如表2所示.體積約束均為30%，低導熱系數為λmin＝1 W／（m·K），熱傳導系數比為400.圖16是4種工況參數下的拓撲優化結果.從圖16（a）中可以看出，高導熱材料集中在兩種邊界的最短路徑上，隨著Dirichlet邊界長度L1的增加，中間材料向兩邊分開，如圖16（b）所示.當L2的邊界長度也由0.01m 增加到0.10m 時，高導熱材料越來越多地集中在熱量的入口邊界上，如圖16（c）所示.由此可知，當模型中只有Dirichlet邊界和Neumann邊界時，高導熱材料主要分布在Dirichlet邊界和Neumann邊界的最短路徑上.

表2 工況參數Table.2 Parameters of working condition

圖16 混合邊界下的優化結果Fig.16 Results of different length of boundary

3 結 論

應用密度法設計傳熱結構中散熱通道的分布，建立了在穩態熱傳導下的拓撲優化數學模型，推出了滿足KKT 條件的傳熱結構拓撲優化的最優準則法迭代公式，通過對Dirichlet邊界和Neumann邊界的不同組合下的數值算列，驗證了算法的有效性.在拓撲結果的形態上，Dirichlet邊界和非集中生熱率載荷時，高導熱材料分布形態為類似植物分枝結構的拓撲形態；Dirichlet邊界和集中生熱率載荷時，高導熱材料分布在集中載荷點與Dirichlet邊界的最短路徑上；Dirichlet邊界和Neumann邊界的混合邊界時，高導熱材料分布在兩種邊界的最短路徑上，且集中在Neumann邊界附近.

［1］Bends?e M P.Optimal shape design as a material distribution problem［J］.Structural Optimization，1989，4（1）：193－202.

［2］Bends?e M P，Kikuchi N.Generating optimal topologies instructuraldesign using a homogenizationmethod［J］.Computer Methods in Applied Mechanics and Engineering，1988，71（2）：197－224.

［3］Fujii D，Chen B C，Kikuchi N.Composite material design of two－dimensional structures using the homogenization design method［J］.International Journal for Numerical Methods in Engineering，2001，50（9）：2031－2051.

［4］Wang M，Wang X，Guo D M.A level set method for structural topology optimization ［J］.Computer Methods in Applied Mechanics and Engineering，2003，192（1／2）：227－246.

［5］Allaire G，Jouve F，Toader A.Structural optimization using sensitivity analysis and a level－set method［J］.Journal of Computational Physics，2004，194（1）：363－393.

［6］Bends?e M P，Sigmund O.Topology optimization theory，methods and applications［M］.Berlin：Springer Verlag，2003.

［7］Iga A，Nishiwaki S，Izui K，et al.Topology optimization for thermal problems based on assumed continuous approximation of material distributions［J］.Transacions of the Japan Society of Mechanical Engineers，Series C，2007，73（733）：2426－2433.

［8］Yamata T，Nishiwaki S.Level set－based topolog optimization method for thermal problems［J］.Transacions of the Japan Society of Mechanical Engineers，Series C，2009，75（759）：2868－2876.

［9］Iga A，Yamada T，Nishiwaki S，et al.Topogy optimization for coupled thermal and srtuctural problems using the level set method［J］.Transactions of the Japan Society of Mechanical Engineers，Series C，2010，76（761）：36－43.

［10］Bejan A.Constructal－theory network of conducting paths for cooling a heat generating volume［J］.International Journal of Heat and Mass Transfer，1997，40（4）：799－816.

［11］Gersborg－Hansen A，Bends?e M P，Sigmund O.Topology optimization of heat conduction problems using the finite volume method［J］.Structural and Multidisciplinary Optimization，2004，31（4）：251－259.

［12］Ding X H，Yamazaki K.Stiffener layout design for plate structures by growing and branching tree model（application to vibration－proof design）［J］.Structural and Multidisciplinary Optimization，2004，26（2）：99－110.

［13］Ding X H，Yamazaki K.Constructal design of cooling channel in heat transfer system by utilizing optimality of branch systems in nature［J］.Journal of Heat Transfer，2007，129（3）：245－255.

［14］Bends?e M，Sigmund O.Material interpolation schemes in topology optimization［J］.Archive of Applied Mechanics，1999，69（9／10）：635－654.

［15］Rozvany G I N.Aims，scope，methods，history and unified terminology of computer－aided topology optimization in structural mechanics［J］.Structural and Multidisciplinary Optimization，2001，21（2）：90－108.

［16］左孔天，陳立平，張云清，等.用拓撲優化方法進行熱傳導散熱體的結構優化設計［J］.機械工程學報，2005，41（4）：13－16.

［17］Sigmund O，Petersson J.Numerical instabilities in topology optimization：a survey on procedures dealing with checkerboards，mesh－dependencies and local minima［J］.Structural Optimization，1998，16（1）：68－75.

［18］Diaz A，Sigmund O.Checkerboard patterns in layout optimization［J］.Structural Optimization，1995，10（1）：40－45.

［19］袁振，吳長春.采用非協調元的連續體拓撲優化設計［J］.力學學報，2003，35（2）：176－180.