周期性功能梯度結構穩態熱傳導拓撲優化設計

李信卿 趙清海， 張洪信 張鐵柱 陳建良

1.青島大學機電工程學院，青島，2660712.青島大學動力集成及儲能系統工程技術中心，青島，266071

0 引言

功能梯度結構作為一種新型非均勻結構，具有緩沖吸振、隔熱降噪等優異的綜合性能，已得到航空航天、軌道交通、通信電子等行業的廣泛關注，近年來逐漸成為跨結構、材料領域的研究熱點。自20世紀80年代BENDS?E和KIKUCHI[1]首次提出連續體結構拓撲優化概念以來，拓撲優化方法已發展出均勻化法、漸近結構優化法[2]、變密度法[3]、獨立連續映射法[4]、水平集法[5]和移動變形組件法[6-7]等，各類方法相輔相成推動了宏/微觀結構的協同發展。PAULINO等[8]率先通過準則法搭建了功能梯度材料的固體各向同性微結構懲罰(solid isotropic microstructure with penalization，SIMP)模型，探討了功能梯度結構拓撲優化設計的有效性。此后，功能梯度宏觀結構拓撲優化不斷得到發展和完善。邱克鵬等[9]采用凸規劃求解策略，獲得了功能梯度夾層結構中夾芯的拓撲構型。楊旭靜等[10]結合節點梯度位移函數，建立了功能梯度材料結構拓撲優化近似顯式模型。CHENG等[11]提出了一種應力約束下基于結構拓撲優化設計的梯度點陣結構新方法，以生成具有可預測屈服性能的輕量化點陣結構。上述研究拓展了功能梯度結構拓撲優化方法，但未考慮微觀材料的結構、布置對宏觀結構的影響，極大限制了微觀材料性能的充分利用。

鑒于此，部分學者在宏觀結構拓撲優化中施加周期性約束，為發掘微觀材料結構提供了新思路。SIGMUND[12-13]運用逆均勻化技術實現了具有極值力學性能(如零剪切性能、負泊松比)的周期性材料微結構設計。除均勻化法外，焦洪宇等[14]提出了一種基于變密度理論SIMP模型的周期性拓撲優化方法。杜義賢等[15]以周期性單胞為研究對象，構建出了基于宏觀力學性能的細觀點陣結構優化模型。付志方等[16]通過搭建線彈性結構周期性穩健拓撲優化模型，得到了具有良好穩定性的周期性結構。此外，周期性設置在熱傳導結構上的應用已逐漸成為新的研究熱點。龍凱等[17]通過建立周期性結構拓撲優化模型，實現了傳熱微觀復合材料優化設計。賈嬌等[18]基于SIMP周期性結構模型，驗證了具有宏觀導熱性能的周期性材料的設計可行性問題。趙清海等[19]考慮周期性約束下的多材料結構穩態熱傳導拓撲優化設計模型，探討了一種周期性多材料結構拓撲設計方法。上述文獻探討了宏觀結構中最優微觀材料結構，但并未考慮宏/微觀結構之間的最佳性能匹配設計，往往難以實現結構的綜合性能最優。

在此背景下，周期性功能梯度結構拓撲優化設計應運而生，它可綜合考慮宏/微觀結構的性能匹配。ZHOU等[20]借助兩相梯度微結構的逆均勻化方法，確定了周期性基底單元沿平行于梯度的變化方向。DENG等[21]搭建并行優化模型應用到輕量化結構拓撲和材料微觀結構設計，探究了功能梯度結構的多目標設計問題。GUO等[22]研究了多尺度框架下有界載荷不確定性材料和結構的魯棒并行優化問題。此類設計方法主要針對由單一材料微結構組成的宏觀結構設計，研究具有多層微結構的功能梯度拓撲優化時，則根據材料性能需要在不同區域分布不同的微結構，這可增加結構的多樣性。張衛紅等[23]針對功能梯度材料，給出了與尺寸關聯的結構和材料拓撲優化設計結果。COELHO等[24]基于層次計算方法探討了特定局部微結構的等效材料特性和層次優化設計問題，并實現了三維功能梯度結構設計。XIA等[25]實現了功能梯度結構材料性能和拓撲布局的并行優化設計。HUANG等[26]采用雙向進化結構優化方法對材料單元進行靈敏度分析，得到了多孔材料和復合材料的各向異性微結構。LI等[27]借助蜂窩復合材料和功能梯度材料的優勢，提出了一種用于集成材料和結構設計的分層多尺度的拓撲優化方法。上述研究集中在結構場下力學性能對宏觀結構和微觀構型的影響，有關熱傳導功能梯度方面的研究報道較少。相對于力場問題，受熱載荷影響熱場結構的拓撲優化問題則更加復雜[28]。在散熱結構設計中，應綜合考慮周期性與功能梯度約束，實現裝配簡單、生產成本低、易于模塊化設計等周期性結構制造加工的同時，并通過對結構進行梯度分層來保證結構具有良好的散熱性能。

本文提出了一種周期性約束下的功能梯度結構穩態熱傳導拓撲優化方法。采用基于SIMP模型的材料插值方法，利用熱傳導結構提取最優構型中各預設梯度層的體積分數，通過重新分配散熱弱度來實現周期性約束設置，同時采用移動漸近線法推導出設計變量的迭代公式，借助2D和3D數值算例，探討周期性約束下功能梯度穩態熱傳導結構的材料分布、各區域的宏微觀結構以及優化結果。

1 熱傳導周期性功能梯度結構描述

對于熱傳導周期性功能梯度結構，宏觀設計域被劃分為多個梯度層(L1，…，Lj，…，Ln)，其中j(j=1，2，…，n)為梯度層序號，n為劃分的總梯度層數，每個梯度層由周期性材料微結構組成。首先給定熱傳導結構的Dirichlet邊界ГD、Neumann邊界ГN和內熱源等初始條件，對宏觀結構進行拓撲優化設計，獲得各梯度層材料密度最優結果；進而對梯度層宏觀結構和微觀材料協同優化設計，獲得周期性功能梯度材料微結構最佳分布，實現功能梯度拓撲優化設計。熱傳導周期性功能梯度結構如圖1所示。

圖1 熱傳導周期性功能梯度結構Fig.1 Heat conduction periodic functional gradient structure

2 周期性功能梯度結構的建模方法

2.1 基于SIMP的傳統熱傳導結構建模

針對熱傳導結構設計，基于傳統SIMP法建立的相對密度插值模型為

(1)

式中，xi為單元i的相對密度(即設計變量)，xi∈[0，1]；λ為插值后的熱導率；λ0、λv分別為材料初始熱導率和孔洞區域熱導率；p為懲罰因子。

建立穩態熱傳導的結構拓撲優化數學模型：

(2)

式中，C為結構散熱弱度；Q、T和Kλ分別為熱載荷向量、節點溫度向量和整體熱傳導剛度矩陣；Ne為設計域單元總數；vi為單元i的體積；V、f、V0分別為優化后體積、設計域的體積比和總體積；xmin為單元相對密度的下限值。

2.2 功能梯度結構建模和求解

為計算體積約束下功能梯度結構的總體布局，借助熱傳導結構拓撲優化模型，可提取最優拓撲構型中各預設梯度層的體積分數。將二維矩形設計域劃分為n個梯度層(L1，L2，…，Lj，…，Ln)，分層梯度設置如圖2所示。

圖2 分層梯度設置示意圖Fig 2 Layered gradient setting diagram

針對功能梯度結構設計，建立SIMP插值模型，可表示為

(3)

式中，xi，j為第j梯度層單元i的相對密度，xi，j∈[0，1]。

其拓撲優化數學模型可表示為

(4)

式中，N為梯度層內的單元數；ti，j、kλ分別為第j梯度層單元i的節點溫度向量和單元熱傳導剛度矩陣；vi，j為第j梯度層單元i的體積；fj、Vj分別為第j梯度層的體積比和總體積。

熱傳導結構拓撲優化設計靈敏度求解過程中，利用伴隨變量法推導出目標函數C相對于單元密度xi，j的靈敏度，即

(5)

不考慮溫度載荷的相關性，則?QT/?xi，j=0，根據平衡方程KλT=Q對設計變量xi，j求導，可得

(6)

將式(6)代入式(5)可得

(7)

體積約束V相對于單元密度xi，j的導數為

(8)

2.3 設計變量更新和靈敏度過濾

目前常用的優化求解算法包括數學規劃法和優化準則法，由于移動漸近線法(method of moving asymptotes，MMA)[29]對復雜目標和多約束拓撲優化問題具有更好的適定性，因此，本文選取MMA法進行功能梯度結構拓撲優化設計求解。

對于拓撲優化結果中存在的數值不穩定現象，目前常見的解決方法包括周長約束法、過濾法等。基于Helmholts方程的偏微分方程(partial differential equation，PDE)過濾方法具有高效的計算效率，適合用于大規模的優化求解問題。本文借助PDE過濾方法進行靈敏度過濾，其過濾方程組為

(9)

(10)

(11)

2.4 周期性結構設置

為獲得周期性功能梯度結構拓撲優化形式，可將每個梯度層劃分為Mxj×Myj個相同的子區域(其中Mxj、Myj分別為第j梯度層在x與y軸方向的子區域數)，使得各梯度層不同子結構在相同位置處的單元具有相同的材料屬性。如圖3所示，梯度層內周期性結構的數學模型可描述為

(12)

圖3 周期性分層設置示意圖Fig.3 Schematic diagram of periodic layered setup

(13)

圖4 周期性功能梯度拓撲優化流程圖Fig.4 Flowchart of periodic functional gradient structure

3 算例

3.1 2D算例

算例1 以二維熱傳導問題為例，宏觀設計域如圖5a所示，選取幾何尺寸為180 mm×180 mm 的正方形結構，將其離散為180×180個四節點四邊形單元。在整個設計域施加強度為Q=0.01 W/mm2的均布熱源，左邊界中間位置給定溫度為零，其余邊界絕熱，宏觀設計域的體積分數為0.55，懲罰因子p=3，過濾半徑rmin=1.2 mm。分層設計域如圖5b所示，沿x軸方向分別對設計域進行梯度預設，給定梯度層數n=3，5，10，由宏觀拓撲優化可獲得三種預設梯度層的體積分數如表1所示。

(a)宏觀設計域 (b)分層設計域圖5 二維設計域Fig.5 Two-dimensional design domain

表1 各梯度層的體積分數Tab.1 Volume fraction of gradient layers

全局周期性約束下的拓撲構型如圖6所示。周期性約束為My=3，5，10；Mx=1，以此進行周期性功能梯度拓撲優化研究。功能梯度結構的拓撲構型以及溫度分布分別如圖7和圖8所示，并將梯度層子結構的優化結果列于表2。

由圖6可知，對于全局周期性約束下的拓撲構型，隨子區域數量的增加，其微觀子區域的細節特征趨于簡潔，宏觀結構則呈現類似“框架”的形式。由圖7可知，拓撲構型呈現一種由絕熱邊界向零溫點邊界聚攏的“樹枝”狀分布，且在周期性約束相同的情況下，隨梯度層數的增加，結構的拓撲構型差異愈發明顯。由表2可以發現，每個梯度層的子區域均呈現出材料最優化分布。圖8中，由藍色到黃色的溫度分布表示溫度由低到高。由圖8可知，零溫處的溫度最低，離零溫處越遠其溫度呈梯度升高，具有一定的周期性特征，且溫度分布較為均勻。

(a)My=Mx =1 (b)My=Mx =3

(a)My=3；Mx=1

(a)My=3；Mx=1

圖9為全局周期(圖6c)和周期性分層梯度(圖7a～圖7c中梯度層數劃分為n=5時)設置下散熱弱度隨迭代次數的變化曲線，可以發現分層梯度會對結構的散熱效果產生影響，在周期性約束相同的情況下，結構的散熱性能隨梯度層數的增加而降低；而對比全局周期以及周期性分層梯度設置下結構的散熱性能發現，后者的散熱性能明顯優于前者的散熱性能，說明周期性功能梯度結構具有優越性。由圖9還可知，分層梯度設置下目標函數曲線經過迭代更新可收斂到最優解，且收斂速度較快，驗證了所提方法的有效性。

表2 2D拓撲構型中梯度層內子結構Tab.2 Gradient layer sub-structure in 2D topological configurations

圖9 目標性能隨迭代次數的變化曲線Fig.9 Variation curve of target performance with number of iterations

算例2 選取幾何尺寸為180 mm×180 mm的正方形結構，分別離散為180×180個四節點四邊形單元，將強度為Q=324 W/mm2的內熱源施加在設計區域(45 mm，45 mm)、(1 mm，90 mm)和(90 mm，90 mm)位置點，四個角點的溫度恒為零，四邊界絕熱宏觀設計域的體積分數為0.55，懲罰因子p=3，過濾半徑rmin=1.2 mm。

宏觀結構設計域如圖10所示，沿x軸方向分別對設計域進行梯度預設，給定梯度層數n=1，3，6，由宏觀拓撲優化可獲得三種預設梯度層的體積分數。周期性約束為My=3；Mx=1，不同熱源位置下功能梯度結構優化構型結果及其溫度分布分別如圖11和圖12所示，功能梯度結構的散熱性能對比見表3。

(a)(45，45)mm (b)(1，90)mm (c)(90，90)mm圖10 結構設計域Fig.10 Structural design domain

(a)(45 mm，45 mm)位置

(a)(45 mm，45 mm)位置

由圖11可知，不同受熱位置下熱傳導結構的拓撲構型有明顯差異。由表3可以發現，同一分層梯度設置下(即梯度層數相同時)，不同熱源位置結構的散熱性能會有較大差異，這說明熱源位置可能會影響周期性功能梯度結構的最短散熱路徑。由圖12可知，溫度最高點出現在熱源位置，整個設計區域的溫度分布較為均勻，整體溫度梯度較小，內熱源附近的溫度梯度較大。

3.2 3D算例

以三維熱傳導問題為例，如圖13所示，選取幾何尺寸為120 mm×6 mm×120 mm的長方體結構，將其離散為120×6×120個八節點六面體單元，三維結構下端面中間位置邊界溫度設為零，其余邊界表面絕熱。將強度為Q=0.01 W/mm2的均布熱源施加在整個三維區域，體積分數約束為0.55，懲罰因子p=3，過濾半徑rmin=1.4 mm，沿z軸方向分別對設計域進行梯度預設，給定梯度層數n=3，6，由宏觀拓撲優化可獲得兩種預設梯度層的體積分數。周期性約束為Mx=3，6；My=Mz=1，功能梯度結構拓撲構型如圖14所示，圖14b中梯度層數劃分為n=3時梯度層內的子結構如圖15所示。

表3 功能梯度結構散熱性能對比Tab.3 Comparison of heat dissipation performance of functional gradient structures

圖13 三維宏觀設計域和結構拓撲構型Fig.13 Three-dimensional macro design domain and macro-structural topology configurations

(a)My=Mz=1；Mx=3

圖15 3D拓撲構型中梯度層內子結構Fig.15 Gradient layer sub-structure in 3D topological configurations

由圖14可以看出，周期性功能梯度結構子區域構型呈現出了與圖13相似的“樹枝狀”拓撲構型。在周期性約束相同的情況下，隨著梯度層數的增加，結構拓撲構型越發明顯且邊界清晰合理，散熱性能逐漸降低，從而證明了所提方法可實現熱傳導結構周期性功能梯度拓撲優化設計，保證了該方法在工程上的適用性。

以圖14a為例，三維熱傳導結構的散熱曲線隨迭代次數變化情況見圖16，可以看出，當迭代到60步后，變化曲線趨于平穩。不同周期性約束下功能梯度熱傳導結構的散熱弱度變化不大，與二維設計相對應，進一步驗證了功能梯度結構具有良好的散熱性能。

圖16 目標性能隨迭代次數的變化曲線Fig.16 Variation curve of target performance with number of iterations

4 結論

(1)提出了一種考慮周期性約束的功能梯度結構穩態熱傳導拓撲優化設計方法，對比探討了全局周期以及周期性功能梯度拓撲優化，通過2D與3D數值算例驗證了所提方法的有效性。

(2)通過對宏觀設計區域進行梯度分層設置，并結合固體各向同性微結構懲罰(SIMP)插值方法和重新分配散熱弱度基值，可獲得功能梯度結構周期性熱傳導構型。

(3)數值算例結果表明，各梯度層不同設計變量設置下，結構在宏微觀構型方面具有差異性，表明了宏觀性能約束下，梯度層微觀結構的改變會對宏觀拓撲構型產生影響。

(4)采用所提方法獲得的周期性功能梯度結構具有構型材料分布清晰、方法設置簡便的優點，且有廣闊的工程適用前景。