基于雙曲正弦函數(shù)的拓?fù)鋬?yōu)化在熱傳導(dǎo)結(jié)構(gòu)設(shè)計(jì)中的應(yīng)用

薛紅軍，陶才勇，黨思娜

(西北工業(yè)大學(xué) 航空學(xué)院， 西安 710072)

0 引 言

飛機(jī)電子設(shè)備不斷增多，其散熱結(jié)構(gòu)是航空工程結(jié)構(gòu)設(shè)計(jì)過(guò)程中必須考慮的重要問(wèn)題之一[1]，飛機(jī)中的傳熱結(jié)構(gòu)設(shè)計(jì)既要滿足一定的強(qiáng)度和剛度要求以達(dá)到特定的功能，又要滿足散熱要求，目的是保證關(guān)鍵設(shè)備不會(huì)因?yàn)闇囟冗^(guò)高而出現(xiàn)故障。通過(guò)熱力學(xué)計(jì)算設(shè)計(jì)的傳統(tǒng)散熱裝置，已無(wú)法滿足對(duì)熱防護(hù)結(jié)構(gòu)具有高性能要求的實(shí)際工程需求[2]，因此，尋求一種新的散熱裝置設(shè)計(jì)方法成為解決熱傳導(dǎo)問(wèn)題的關(guān)鍵。

針對(duì)熱傳導(dǎo)問(wèn)題，H.Rodrigues等[3-4]和C.Jog[5]將拓?fù)鋬?yōu)化應(yīng)用在熱彈性結(jié)構(gòu)中，初次實(shí)現(xiàn)了拓?fù)鋬?yōu)化在熱傳導(dǎo)中的應(yīng)用；A.Bejan[6]在熱傳導(dǎo)的背景下介紹了熱傳導(dǎo)通道構(gòu)造理論，建立了所謂的“體-點(diǎn)(VP)問(wèn)題”或“通道問(wèn)題”，并討論了在發(fā)熱域(例如CPU)中布置固定數(shù)量的材料的需要；Li Q等[7]明確地將ESO用于熱傳導(dǎo)問(wèn)題。

現(xiàn)階段應(yīng)用于熱傳導(dǎo)結(jié)構(gòu)的拓?fù)鋬?yōu)化方法主要有：漸進(jìn)結(jié)構(gòu)優(yōu)化方法(ESO)、雙向漸進(jìn)結(jié)構(gòu)優(yōu)化方法(BESO)[8]、變密度方法[9]、水平集方法[10]、拓?fù)鋵?dǎo)數(shù)方法[11]。就目前的研究成果而言，應(yīng)用比較廣泛的拓?fù)鋬?yōu)化方法是變密度法。

基于變密度方法的散熱結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì)主要使用：帶懲罰的各向同性固體材料模型[12](Solid Isotropic Material with Penalization，簡(jiǎn)稱SIMP)和材料屬性的理性近似模型[13](Rational Approximation of Material Properties，簡(jiǎn)稱RAMP)。F.H.Burger等[14]利用SIMP和移動(dòng)漸近線(MMAs)的方法，探討了體-點(diǎn)問(wèn)題的3D解決方案；T.Van Oevelen等[15]采用SIMP模型演示了用雙層簡(jiǎn)化模型表示全三維解的共軛傳熱問(wèn)題的解；G.Marck等[16]使用SIMP方法進(jìn)行了多目標(biāo)優(yōu)化；Li Jiachun等[17]使用了材料特性的理性近似模型(RAMP)，進(jìn)行了熱傳導(dǎo)結(jié)構(gòu)優(yōu)化設(shè)計(jì)的探索。兩種模型對(duì)于導(dǎo)熱結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì)都能得到比較好的布局，但是對(duì)于“體-點(diǎn)”問(wèn)題[18]，SIMP方法會(huì)產(chǎn)生很多樹(shù)狀末梢結(jié)構(gòu)，這樣的結(jié)果有一定的缺點(diǎn)，如可制造性差，并且結(jié)果并不是最靠近最優(yōu)值的結(jié)果[19]。因此有學(xué)者通過(guò)RAMP方法得到了一定的改進(jìn)，能過(guò)獲取相對(duì)更加良好的拓?fù)浣Y(jié)果[20]。但是，這兩種方法有其局限性，在解決二維、三維均勻熱源的散熱問(wèn)題時(shí)，SIMP方法不夠精確，RAMP方法迭代步長(zhǎng)過(guò)多。

本文在變密度法的基礎(chǔ)上，建立基于sinh函數(shù)的插值模型，與SIMP和RAMP模型的計(jì)算結(jié)果進(jìn)行對(duì)比，并將此模型應(yīng)用于機(jī)載LRM模塊導(dǎo)熱拓?fù)鋬?yōu)化設(shè)計(jì)中進(jìn)行驗(yàn)證。

1 基于sinh函數(shù)插值模型的拓?fù)鋬?yōu)化在熱傳導(dǎo)結(jié)構(gòu)設(shè)計(jì)中的模型建立

基于變密度法的拓?fù)鋬?yōu)化的關(guān)鍵是選擇合適的插值函數(shù)和懲罰技術(shù)，將問(wèn)題的物理量表示為連續(xù)設(shè)計(jì)變量的函數(shù)。本文采用基于sinh函數(shù)的插值模型[19]進(jìn)行求解計(jì)算。

1.1 基于sinh函數(shù)的插值模型

基于sinh函數(shù)的插值模型，其實(shí)質(zhì)是在插值模型中應(yīng)用了雙曲正弦函數(shù)(sinh函數(shù))：

(1)

可以簡(jiǎn)化為

(2)

sinh函數(shù)插值模型采用雙曲正弦函數(shù)來(lái)對(duì)中間密度進(jìn)行懲罰，其函數(shù)關(guān)系如圖1所示。

圖1 sinh插值模型

sinh函數(shù)插值模型的微分形式為

(3)

當(dāng)ρ=0時(shí)，

當(dāng)ρ=1時(shí)，

1.2 sinh插值模型靈敏度分析

計(jì)算導(dǎo)數(shù)最好的辦法是使用伴隨方法，在此也采用伴隨方法推導(dǎo)sinh函數(shù)插值模型的目標(biāo)函數(shù)的靈敏度。

本文以熱柔度為目標(biāo)函數(shù)，基于sinh函數(shù)插值模型的目標(biāo)函數(shù)為

(4)

式中：K0為單元導(dǎo)熱系數(shù)矩陣。

對(duì)于上面給出的熱柔度目標(biāo)函數(shù)，通過(guò)添加零函數(shù)，將其重寫為

C=QTT-λT(KT-Q)

(5)

式中：λ為拉格朗日乘數(shù)，是一個(gè)任意的、但是固定的實(shí)數(shù)向量。

然后將其對(duì)設(shè)計(jì)變量求導(dǎo)得到：

(6)

當(dāng)λ滿足伴隨方程

QT-λTK=0

(7)

目標(biāo)函數(shù)靈敏度可寫為

(8)

根據(jù)伴隨方程，可以直接得到λ=T，因此目標(biāo)函數(shù)靈敏度為

(9)

其中，對(duì)于傳熱系數(shù)矩陣對(duì)設(shè)計(jì)變量的導(dǎo)數(shù)為

(10)

因此，當(dāng)傳熱系數(shù)矩陣采用sinh函數(shù)插值模型對(duì)其進(jìn)行插值，則目標(biāo)函數(shù)的靈敏度為

(11)

式中：K0為具有單元熱傳導(dǎo)矩陣，當(dāng)材料的導(dǎo)熱系數(shù)為k0時(shí)，對(duì)于二維的四節(jié)點(diǎn)矩形單元，單元的尺寸為單位長(zhǎng)度，K0可寫為

(12)

對(duì)于三維八節(jié)點(diǎn)正六面體單元，單元的尺寸為單位長(zhǎng)度，單元導(dǎo)熱系數(shù)矩陣K0可寫為

(13)

當(dāng)懲罰因子p=3時(shí)，靈敏度隨著設(shè)計(jì)變量的變化曲線如圖2所示，可以看出：所有的靈敏度值

都是負(fù)值，也就意味著隨著導(dǎo)熱材料的增加，目標(biāo)函數(shù)在不斷地減小，導(dǎo)熱性能更好。

圖2 靈敏度變化曲線

為了避免棋盤格現(xiàn)象和網(wǎng)格依賴現(xiàn)象，在計(jì)算中采用靈敏度過(guò)濾技術(shù)，將靈敏度改寫為

(14)

式中：Hei為單元i的一個(gè)集合；γ為一個(gè)比較小的常數(shù)，為避免分母為0，可取γ=10-3；單元i中心到單元e的距離Δ(e,i)小于最小濾波半徑rmin。

1.3 基于sinh函數(shù)插值模型的優(yōu)化準(zhǔn)則

本文以結(jié)構(gòu)的散熱弱度最小為目標(biāo)函數(shù)，體積分?jǐn)?shù)、結(jié)構(gòu)平衡方程和設(shè)計(jì)變量約束為約束條件，建立基于sinh函數(shù)插值模型的熱傳導(dǎo)結(jié)構(gòu)的拓?fù)鋬?yōu)化數(shù)學(xué)模型為

min:

(15)

式中：C為散熱弱度；K為結(jié)構(gòu)的總體傳熱系數(shù)矩陣；F為溫度載荷向量；T為節(jié)點(diǎn)溫度向量；Ve為單元體積；V0為結(jié)構(gòu)初始狀態(tài)下總的體積；f為體積分?jǐn)?shù)，是優(yōu)化后的結(jié)構(gòu)總體積與初始狀態(tài)下總體積之比；ximin為設(shè)計(jì)變量取值的下限，通常取一個(gè)較小的數(shù)，如0.001；ximax為設(shè)計(jì)變量取值的上限，取為ximax=1。

利用拉格朗日乘數(shù)法將多為高次多約束問(wèn)題轉(zhuǎn)換為無(wú)約束最優(yōu)化問(wèn)題，對(duì)應(yīng)于拓?fù)鋬?yōu)化問(wèn)題的拉格朗日函數(shù)為

(16)

式中：λ1為拉格朗日乘子；λ2，λ3，λ4為拉格朗日乘子向量；bi和ci為松弛因子；x為xi組成的列向量。

(17)

(18)

考慮上式等于0的情況，并且代入C=TTK(X)T有：

(19)

利用導(dǎo)熱系數(shù)矩陣的對(duì)稱性得出：

(20)

式中：λ2為列矢量，其取值無(wú)限制。取λ2=-2T代入式(20)得出：

(21)

進(jìn)而得出：

(22)

因此可以得到基于優(yōu)化準(zhǔn)則法的迭代公式為

(23)

式中：m為移動(dòng)極限，通常取一個(gè)較小的正數(shù)；η為阻尼系數(shù)，引入η的目的是為了保證數(shù)值計(jì)算的穩(wěn)定性和收斂性。

(24)

(25)

將式(21)代入式(25)得到：

(26)

通過(guò)求解上述方程便可以求解出λ1，通常采用二分法進(jìn)行求解。

設(shè)計(jì)變量通過(guò)優(yōu)化準(zhǔn)則更新，其收斂性的檢查方法如式(27)所示：

|max (xk+1-xk)|<ε

(27)

式中:xk+1為更新的設(shè)計(jì)變量；xk為舊的設(shè)計(jì)變量；ε為可行性誤差，即每個(gè)單元更新后的密度值與前一次的密度值之差的絕對(duì)值不超過(guò)規(guī)定的誤差值。

2 模型算法實(shí)現(xiàn)以及算例分析

2.1 算法實(shí)現(xiàn)

本文在MATLAB上進(jìn)行程序編譯，實(shí)現(xiàn)模型功能，整體的實(shí)現(xiàn)步驟如圖3所示。

圖3 基于sinh的熱傳導(dǎo)結(jié)構(gòu)拓?fù)鋬?yōu)化方法算法流程

2.2 算例分析

2.2.1 二維算例

(1) 體-點(diǎn)問(wèn)題

算例1：整塊板均勻加熱，左邊界中心處給定溫度T=0 ℃，其余邊界絕熱，其結(jié)構(gòu)示意圖如圖4所示。

圖4 導(dǎo)熱板均勻加熱，左側(cè)中心處冷卻模型示意圖

高熱導(dǎo)率材料的導(dǎo)熱系數(shù)為k0=1 W/m·℃的假想材料，低熱導(dǎo)率材料的導(dǎo)熱系數(shù)為kmin=0.001 W/m·℃的假想材料，模型離散網(wǎng)格為180×180。優(yōu)化參數(shù)設(shè)置為：體積分?jǐn)?shù)0.4；懲罰因子p=3；采用敏度濾波，濾波半徑r=2。分別采用SIMP模型，RAMP模型和基于sinh函數(shù)模型對(duì)其進(jìn)行優(yōu)化計(jì)算，其優(yōu)化結(jié)果如圖5所示，對(duì)比結(jié)果如表1所示。

(a) SIMP (b) RAMP (c) sinh模型

圖5 導(dǎo)熱板均勻加熱，左側(cè)中心處冷卻模型優(yōu)化結(jié)果

Fig.5 Optimization results of heat-conducting plate evenly heating and cooling at the left center

表1 導(dǎo)熱板均勻加熱，左側(cè)中心處冷卻模型優(yōu)化數(shù)據(jù)對(duì)比

(2) 均勻熱源問(wèn)題其他邊界條件

算例2：整塊板均勻加熱，左邊界給定溫度T=0 ℃，其余邊界絕熱，其結(jié)構(gòu)示意圖如圖6所示。高熱導(dǎo)率材料的導(dǎo)熱系數(shù)為k0=1 W/m·℃的假想材料，低熱導(dǎo)率材料的導(dǎo)熱系數(shù)為kmin=0.001 W/m·℃的假想材料，優(yōu)化參數(shù)設(shè)置為：體積分?jǐn)?shù)0.5；懲罰因子p=3；采用敏度濾波，濾波半徑r=3。分別采用SIMP模型，RAMP模型和基于sinh函數(shù)模型對(duì)其進(jìn)行優(yōu)化計(jì)算，其優(yōu)化結(jié)果如圖7所示，對(duì)比結(jié)果如表2所示。

圖6 導(dǎo)熱板均勻加熱，左側(cè)邊界冷卻模型示意圖

(a) SIMP (b) RAMP (c) sinh模型

圖7 導(dǎo)熱板均勻加熱，左側(cè)邊界冷卻模型優(yōu)化結(jié)果

Fig.7 Optimization results of heat-conducting plate evenly heating, and the left edge cooling

表2 導(dǎo)熱板均勻加熱，左側(cè)邊界冷卻模型優(yōu)化數(shù)據(jù)對(duì)比

(3) 集中熱源問(wèn)題優(yōu)化分析

算例3：平板中心加熱，四個(gè)角點(diǎn)給定溫度T=0 ℃，其結(jié)構(gòu)示意圖如圖8所示,其優(yōu)化參數(shù)與算例2一致。分別采用SIMP模型，RAMP模型和基于sinh函數(shù)模型進(jìn)行優(yōu)化計(jì)算，優(yōu)化結(jié)果如圖9所示，對(duì)比結(jié)果如表3所示。

圖8 導(dǎo)熱板的中心加熱，四個(gè)頂點(diǎn)冷卻模型示意圖

(a) SIMP (b) RAMP (c) sinh模型

圖9 導(dǎo)熱板的中心加熱，四個(gè)頂點(diǎn)冷卻模型優(yōu)化結(jié)果

Fig.9 Optimization Results of heat-conducting plate center heated and four apexes cooled

表3 導(dǎo)熱板的中心加熱，四個(gè)頂點(diǎn)冷卻模型優(yōu)化結(jié)果數(shù)據(jù)對(duì)比

算例4：板上有5個(gè)集中熱源，四周邊界給定溫度T=0 ℃，其結(jié)構(gòu)示意圖10如圖所示,其優(yōu)化參數(shù)與算例2一致。高熱導(dǎo)率材料的導(dǎo)熱系數(shù)為k0=1 W/m·℃的假想材料，低熱導(dǎo)率材料的導(dǎo)熱系數(shù)為kmin=0.001 W/m·℃的假想材料。優(yōu)化參數(shù)設(shè)置為：體積分?jǐn)?shù)0.5；懲罰因子p=3；采用敏度濾波，濾波半徑r=2.5。分別采用SIMP模型，RAMP模型和基于sinh函數(shù)模型進(jìn)行優(yōu)化計(jì)算，優(yōu)化結(jié)果如圖11所示，對(duì)比結(jié)果如表4所示。

圖10 導(dǎo)熱板上5個(gè)集中熱源加熱，四周邊界散熱模型示意圖

(a) SIMP (b) RAMP (c) sinh模型

圖11 導(dǎo)熱板上5個(gè)集中熱源加熱，四周邊界散熱模型優(yōu)化結(jié)果

Fig.11 Optimization results of the heating of 5 centralized heat sources on the thermal plate and the heat dissipation at the surrounding boundary

表4 導(dǎo)熱板上5個(gè)集中熱源加熱，四周邊界散熱模型優(yōu)化結(jié)果數(shù)據(jù)對(duì)比

2.2.2 三維算例

(1) 均勻熱源問(wèn)題

算例5：設(shè)計(jì)域下面均勻加熱，上平面四個(gè)頂點(diǎn)溫度固定T=0 ℃，模型示意圖如圖12所示。高熱導(dǎo)率材料的導(dǎo)熱系數(shù)為k0=1 W/m·℃的假想材料，低熱導(dǎo)率材料的導(dǎo)熱系數(shù)為kmin=0.001 W/m·℃的假想材料。優(yōu)化參數(shù)設(shè)置為：體積分?jǐn)?shù)0.3；懲罰因子p=3；采用密度濾波，濾波半徑r=1.4。分別采用SIMP模型，RAMP模型和基于sinh函數(shù)模型進(jìn)行優(yōu)化計(jì)算，優(yōu)化結(jié)果如圖13所示，對(duì)比結(jié)果如表5所示。

圖12 下表面均勻加熱，上表面四個(gè)頂點(diǎn)冷卻模型示意圖

(a) SIMP (b) RAMP (c) sinh模型

圖13 下表面均勻加熱，上表面四個(gè)頂點(diǎn)冷卻模型優(yōu)化結(jié)果

Fig.13 Optimization results of uniform heating of the lower surface and cooling of the four vertices of the upper surface

表5 下表面均勻加熱，上表面四個(gè)頂點(diǎn)冷卻模型優(yōu)化結(jié)果數(shù)據(jù)對(duì)比

(2) 集中熱源問(wèn)題

算例6:正方體中心有集中熱源，八個(gè)頂點(diǎn)溫度固定為T=0 ℃，模型示意圖如圖 14所示。高熱導(dǎo)率材料的導(dǎo)熱系數(shù)為k0=1 W/m·℃的假想材料，低熱導(dǎo)率材料的導(dǎo)熱系數(shù)為kmin=0.001 W/m·℃的假想材料。優(yōu)化參數(shù)設(shè)置為：體積分?jǐn)?shù)0.3；懲罰因子p=3；采用密度濾波，濾波半徑r=1.4。分別采用SIMP模型，RAMP模型和基于sinh函數(shù)模型進(jìn)行優(yōu)化計(jì)算，優(yōu)化結(jié)果如圖15所示，對(duì)比結(jié)果如表6所示。

圖14 正方體中心加熱，八個(gè)頂點(diǎn)冷卻模型示意圖

(a) SIMP (b) RAMP (c) sinh模型

圖15 正方體中心加熱，八個(gè)頂點(diǎn)冷卻模型優(yōu)化結(jié)果

Fig.15 Optimization results of cube center heating and eight vertices cooling

表6 正方體中心加熱，八個(gè)頂點(diǎn)冷卻模型優(yōu)化結(jié)果數(shù)據(jù)對(duì)比

通過(guò)算例將SIMP模型和RAMP模型與基于sinh函數(shù)的插值模型進(jìn)行對(duì)比，得出基于sinh函數(shù)的插值模型較SIMP插值模型精確，較RAMP插值模型的迭代步長(zhǎng)少，計(jì)算效率高。

3 機(jī)載LRM模塊導(dǎo)熱拓?fù)鋬?yōu)化設(shè)計(jì)

機(jī)載LRM(現(xiàn)場(chǎng)可更換模塊)模塊已成為飛機(jī)電子系統(tǒng)設(shè)備的硬件標(biāo)準(zhǔn),其優(yōu)點(diǎn)是標(biāo)準(zhǔn)程度高、模塊功能化強(qiáng)、維修性好、成本低。由于機(jī)載LRM模塊是由超大規(guī)模集成電路構(gòu)成的，包含大量的VLSI和VHSIC芯片，而這些設(shè)備在工作時(shí)會(huì)產(chǎn)生大量熱量，使得LRM模塊溫度急劇上升，因此，機(jī)載LRM模塊輕量化的冷卻成為設(shè)計(jì)的首要問(wèn)題。

由于LRM模塊空間的限制，傳導(dǎo)冷卻方法被優(yōu)先考慮。此方法是通過(guò)將熱量從LRM模塊芯片傳至導(dǎo)熱體，再經(jīng)由導(dǎo)熱體傳至機(jī)架冷板，最后通過(guò)通風(fēng)冷卻或通液冷卻等散熱形式將熱量帶走。簡(jiǎn)化后的LRM模塊導(dǎo)熱模型如圖16所示，導(dǎo)熱體材料為銅，導(dǎo)熱系數(shù)為377 W/m·℃。

(a) 側(cè)視圖

(b) 頂視圖

通過(guò)拓?fù)鋬?yōu)化進(jìn)行導(dǎo)熱路徑優(yōu)化設(shè)計(jì)，體積分?jǐn)?shù)0.5；懲罰因子p=3；采用敏度濾波，濾波半徑r=1.2，分別采用SIMP模型、RAMP模型和基于sinh函數(shù)的模型的優(yōu)化結(jié)果如圖17所示，對(duì)比結(jié)果如表7所示。

(a) SIMP

(b) RAMP

(c) sinh模型

表7 LRM模塊導(dǎo)熱優(yōu)化數(shù)據(jù)對(duì)比

Table 7 LRM module thermal conductivity optimization data comparison

通過(guò)對(duì)LRM模塊的導(dǎo)熱體進(jìn)行拓?fù)鋬?yōu)化設(shè)計(jì)，得出本文提出的基于sinh函數(shù)的插值模型較SIMP模型和RAMP模型能夠更加高效快捷地獲取良好的導(dǎo)熱路徑。

4 結(jié) 論

(1) 通過(guò)二維、三維均勻熱源和集中熱源問(wèn)題的算例分析，與SIMP模型和RAMP模型進(jìn)行對(duì)比，得出基于sinh函數(shù)插值模型較SIMP插值模型精確，較RAMP插值模型的步長(zhǎng)少。

(2) 本文采用基于sinh函數(shù)插值模型能夠更加精確快捷地解決散熱問(wèn)題，同時(shí)為傳熱問(wèn)題提供一種新思路。