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

1.青島大學機電工程學院，青島，2660712.青島大學動力集成及儲能系統(tǒng)工程技術(shù)中心，青島，266071

0 引言

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

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

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

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

1 熱傳導周期性功能梯度結(jié)構(gòu)描述

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

圖1 熱傳導周期性功能梯度結(jié)構(gòu)Fig.1 Heat conduction periodic functional gradient structure

2 周期性功能梯度結(jié)構(gòu)的建模方法

2.1 基于SIMP的傳統(tǒng)熱傳導結(jié)構(gòu)建模

針對熱傳導結(jié)構(gòu)設(shè)計，基于傳統(tǒng)SIMP法建立的相對密度插值模型為

(1)

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

建立穩(wěn)態(tài)熱傳導的結(jié)構(gòu)拓撲優(yōu)化數(shù)學模型：

(2)

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

2.2 功能梯度結(jié)構(gòu)建模和求解

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

圖2 分層梯度設(shè)置示意圖Fig 2 Layered gradient setting diagram

針對功能梯度結(jié)構(gòu)設(shè)計，建立SIMP插值模型，可表示為

(3)

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

其拓撲優(yōu)化數(shù)學模型可表示為

(4)

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

熱傳導結(jié)構(gòu)拓撲優(yōu)化設(shè)計靈敏度求解過程中，利用伴隨變量法推導出目標函數(shù)C相對于單元密度xi，j的靈敏度，即

(5)

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

(6)

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

(7)

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

(8)

2.3 設(shè)計變量更新和靈敏度過濾

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

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

(9)

(10)

(11)

2.4 周期性結(jié)構(gòu)設(shè)置

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

(12)

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

(13)

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

3 算例

3.1 2D算例

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

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

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

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

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

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

(a)My=3；Mx=1

(a)My=3；Mx=1

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

表2 2D拓撲構(gòu)型中梯度層內(nèi)子結(jié)構(gòu)Tab.2 Gradient layer sub-structure in 2D topological configurations

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

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

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

(a)(45，45)mm (b)(1，90)mm (c)(90，90)mm圖10 結(jié)構(gòu)設(shè)計域Fig.10 Structural design domain

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

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

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

3.2 3D算例

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

表3 功能梯度結(jié)構(gòu)散熱性能對比Tab.3 Comparison of heat dissipation performance of functional gradient structures

圖13 三維宏觀設(shè)計域和結(jié)構(gòu)拓撲構(gòu)型Fig.13 Three-dimensional macro design domain and macro-structural topology configurations

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

圖15 3D拓撲構(gòu)型中梯度層內(nèi)子結(jié)構(gòu)Fig.15 Gradient layer sub-structure in 3D topological configurations

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

以圖14a為例，三維熱傳導結(jié)構(gòu)的散熱曲線隨迭代次數(shù)變化情況見圖16，可以看出，當?shù)?0步后，變化曲線趨于平穩(wěn)。不同周期性約束下功能梯度熱傳導結(jié)構(gòu)的散熱弱度變化不大，與二維設(shè)計相對應(yīng)，進一步驗證了功能梯度結(jié)構(gòu)具有良好的散熱性能。

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

4 結(jié)論

(1)提出了一種考慮周期性約束的功能梯度結(jié)構(gòu)穩(wěn)態(tài)熱傳導拓撲優(yōu)化設(shè)計方法，對比探討了全局周期以及周期性功能梯度拓撲優(yōu)化，通過2D與3D數(shù)值算例驗證了所提方法的有效性。

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

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

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