基于熱固耦合問(wèn)題的多尺度拓?fù)鋬?yōu)化設(shè)計(jì)

魏兆棟， 高仁璟,2， 王長(zhǎng)生*,2

(1.大連理工大學(xué) 汽車工程學(xué)院，工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024;2.大連理工大學(xué) 寧波研究院，寧波 315033)

1 引 言

拓?fù)鋬?yōu)化是一種按照既定的設(shè)計(jì)要求和約束條件，在給定的設(shè)計(jì)區(qū)域內(nèi)對(duì)材料分布進(jìn)行優(yōu)化的數(shù)學(xué)方法和工具，改變了依據(jù)經(jīng)驗(yàn)和直覺(jué)的傳統(tǒng)設(shè)計(jì)方法。隨著理論研究的不斷完善，拓?fù)鋬?yōu)化的研究涉及到多學(xué)科多物理領(lǐng)域并取得了巨大的成果，如熱傳導(dǎo)[1-3]問(wèn)題和靜力學(xué)[4-6]問(wèn)題等。然而在實(shí)際生產(chǎn)過(guò)程中，機(jī)械結(jié)構(gòu)不僅承受外載荷，還經(jīng)常面臨高溫的工作環(huán)境，因此考慮如何對(duì)這類結(jié)構(gòu)進(jìn)行合理拓?fù)鋬?yōu)化是十分必要的。近些年來(lái)，利用逆均勻化方法按照優(yōu)化意圖對(duì)材料微結(jié)構(gòu)進(jìn)行優(yōu)化設(shè)計(jì)來(lái)充分挖掘材料微結(jié)構(gòu)的使用潛力受到關(guān)注，如何在材料微結(jié)構(gòu)和宏觀拓?fù)浣Y(jié)構(gòu)之間建立橋梁，充分發(fā)揮材料和結(jié)構(gòu)的互補(bǔ)優(yōu)勢(shì)并應(yīng)用到熱固耦合結(jié)構(gòu)問(wèn)題是值得研究的熱點(diǎn)。在微結(jié)構(gòu)優(yōu)化設(shè)計(jì)方面，Sigmund[7]首次采用逆均勻化方法，基于拓?fù)鋬?yōu)化手段對(duì)復(fù)合材料的微結(jié)構(gòu)極端性能進(jìn)行研究設(shè)計(jì)，實(shí)現(xiàn)了泊松比在區(qū)間[-1,1]的微結(jié)構(gòu)的最優(yōu)設(shè)計(jì)，這項(xiàng)研究開啟了設(shè)計(jì)具有特定彈性屬性微結(jié)構(gòu)的先河，此后又設(shè)計(jì)了具有極端熱膨脹系數(shù)的三相材料微結(jié)構(gòu)設(shè)計(jì)。Zhang等[8]以材料體積分?jǐn)?shù)為約束條件，基于拓?fù)鋬?yōu)化方法對(duì)復(fù)合材料極端熱傳導(dǎo)性能進(jìn)行優(yōu)化，獲得在不同體分比下微結(jié)構(gòu)的優(yōu)化結(jié)果。以上文獻(xiàn)是采用均勻化方法對(duì)材料微結(jié)構(gòu)有效屬性進(jìn)行預(yù)測(cè)，針對(duì)預(yù)測(cè)材料等效屬性問(wèn)題，Xia等[9]提出了一種基于能量均勻化方法設(shè)計(jì)微結(jié)構(gòu)的方法。Gao等[10]提出了三維的能量均勻化方法，并結(jié)合參數(shù)化水平集優(yōu)化方法對(duì)三維微結(jié)構(gòu)材料的彈性性能進(jìn)行了拓?fù)鋬?yōu)化設(shè)計(jì)。能量均勻化方法和漸進(jìn)均勻化方法在預(yù)測(cè)微結(jié)構(gòu)材料等效屬性方面的結(jié)果是相似的，兩者區(qū)別在于測(cè)試應(yīng)變場(chǎng)處理，漸進(jìn)均勻化方法的應(yīng)變場(chǎng)是由單位測(cè)試應(yīng)變場(chǎng)和由于單胞的非均質(zhì)性引起的應(yīng)變場(chǎng)疊加形成，而基于能量均勻化方法的應(yīng)變場(chǎng)是由直接施加在邊界上單位測(cè)試應(yīng)變場(chǎng)引起的，在實(shí)際編程操作和公式推導(dǎo)中，能量均勻化方法對(duì)于處理微結(jié)構(gòu)的周期性以及材料的應(yīng)變場(chǎng)問(wèn)題是比較清晰和簡(jiǎn)潔的，這也是本文選用該方法的考慮點(diǎn)。在多物理場(chǎng)耦合拓?fù)鋬?yōu)化領(lǐng)域，Sigmund[11,12]研究了涉及電、熱和機(jī)械等多物理場(chǎng)耦合時(shí)柔順機(jī)構(gòu)的拓?fù)湓O(shè)計(jì)方法。Deaton等[13]提出了基于應(yīng)力約束的熱固耦合結(jié)構(gòu)拓?fù)鋬?yōu)化模型。李冬梅等[14]提出了一種基于概率約束的可靠性熱固耦合拓?fù)鋬?yōu)化設(shè)計(jì)方法。Zuo等[15]探討熱固耦合場(chǎng)結(jié)構(gòu)的控制方程，建立了基于SIMP法的耦合場(chǎng)拓?fù)淠Ｐ汀ao等[16]基于雙材料進(jìn)行了質(zhì)量約束下的熱固耦合結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)。Zhao等[17]建立了基于Ordered-RAMP模型的多材料插值模型，對(duì)熱固耦合結(jié)構(gòu)進(jìn)行多材料拓?fù)鋬?yōu)化設(shè)計(jì)，建立了綜合散熱弱度與結(jié)構(gòu)柔度的多目標(biāo)加權(quán)設(shè)計(jì)方案。占金青等[18]提出一種非均勻溫度場(chǎng)下的熱固耦合多相材料的拓?fù)鋬?yōu)化設(shè)計(jì)方法，通過(guò)集成多相材料的方式來(lái)提升結(jié)構(gòu)的性能。孟慶軒等[19]提出了一種在溫度場(chǎng)均勻改變作用下的熱固耦合結(jié)構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化方法。上述文獻(xiàn)的優(yōu)化尺度均停留在宏觀結(jié)構(gòu)方面，沒(méi)有考慮材料微結(jié)構(gòu)的優(yōu)化層次，這對(duì)于充分發(fā)掘材料性能是一種損失。Anaya[20]基于漸進(jìn)均勻化方法計(jì)算微結(jié)構(gòu)的等效屬性，同時(shí)結(jié)合雙向進(jìn)化結(jié)構(gòu)優(yōu)化方法(BESO)提出了一種多尺度熱固耦合拓?fù)鋬?yōu)化算法。與其相比，本文采用能量均勻化方法計(jì)算材料等效屬性，并將其嵌入SIMP優(yōu)化方法中建立一種熱固耦合并行化優(yōu)化方法，對(duì)于宏觀結(jié)構(gòu)的溫度變化量的處理，本文不再將其設(shè)置為一個(gè)常數(shù)，而是依賴于離散單元節(jié)點(diǎn)的溫度；此外進(jìn)一步研究了該并行化優(yōu)化方法對(duì)宏觀結(jié)構(gòu)位移和溫度的影響，以期改善結(jié)構(gòu)的力學(xué)性能。總體來(lái)說(shuō)，針對(duì)熱固耦合結(jié)構(gòu)優(yōu)化問(wèn)題，如何利用拓?fù)鋬?yōu)化手段且不使用優(yōu)異性能材料以有效地改善結(jié)構(gòu)的力學(xué)性能仍是目前研究的難點(diǎn)和熱點(diǎn)。

本文提出一種基于穩(wěn)態(tài)熱源作用下的熱固耦合結(jié)構(gòu)多尺度拓?fù)鋬?yōu)化設(shè)計(jì)，基于能量均勻化結(jié)合 SIMP插值方法建立并行化拓?fù)鋬?yōu)化模型，以結(jié)構(gòu)剛度作為目標(biāo)函數(shù)，以材料體積分?jǐn)?shù)作為約束條件，為了便于數(shù)值計(jì)算，利用直接法推導(dǎo)目標(biāo)函數(shù)及約束的靈敏度，同時(shí)采用Heaviside非線性密度過(guò)濾技術(shù)抑制優(yōu)化過(guò)程的數(shù)值不穩(wěn)定現(xiàn)象，將OC準(zhǔn)則用于優(yōu)化問(wèn)題求解。

2 熱固耦合結(jié)構(gòu)分析

對(duì)非均勻溫度場(chǎng)下的熱固耦合結(jié)構(gòu)進(jìn)行分析，溫度場(chǎng)和結(jié)構(gòu)場(chǎng)系統(tǒng)用有限元方程(1)進(jìn)行描述：

(1)

式中Km和Kt分別為結(jié)構(gòu)的整體剛度矩陣和熱剛度矩陣，T為結(jié)構(gòu)的溫度向量，U為結(jié)構(gòu)的位移向量，F(xiàn)t為加載在結(jié)構(gòu)上的熱通量向量，F(xiàn)m為結(jié)構(gòu)所承受的機(jī)械外載，F(xiàn)ε為由于熱應(yīng)變產(chǎn)生的熱載荷向量。

由文獻(xiàn)[21]可知，結(jié)構(gòu)受熱膨脹引起的初始熱應(yīng)變?chǔ)?可表示為

ε0=αΔTΦT

(2)

式中α為材料的熱膨脹系數(shù)，ΔT為單元溫度變化量,在2D情形中，向量Φ=[1,1,0]，3D問(wèn)題中Φ=[1,1,1,0,0,0]。

由溫度場(chǎng)變化產(chǎn)生的熱載荷向量Fε表示為

(3)

式中N為離散單元的數(shù)目，B為單元的應(yīng)變矩陣，D0為彈性矩陣。假設(shè)材料為各向同性，平面應(yīng)力問(wèn)題的彈性矩陣具體形式如下，

(4)

以2D問(wèn)題為例，對(duì)式(3)展開得熱載荷向量為

[-1-11-111-11]T

(5)

3 多尺度拓?fù)鋬?yōu)化數(shù)學(xué)模型

利用拓?fù)鋬?yōu)化技術(shù)結(jié)合能量均勻化方法，建立以結(jié)構(gòu)剛度為優(yōu)化目標(biāo)的熱固耦合并行化拓?fù)鋬?yōu)化設(shè)計(jì)模型。將OC準(zhǔn)則用于設(shè)計(jì)變量的迭代更新，即分別在每個(gè)尺度中進(jìn)行設(shè)計(jì)變量的優(yōu)化更新，利用能量均勻化計(jì)算每一次迭代后微結(jié)構(gòu)拓?fù)涞牡刃傩裕蛊渥鳛閮蓚€(gè)尺度之間的連接橋梁，其對(duì)應(yīng)的數(shù)學(xué)模型如下，

(i=1,2,…,NM;j=1,2,…,Nm)

s.t.α(uM,vM,DM)=l(vM), ?vM∈H(ΩM,d)

α(um,vm,Dm)=l(vm), ?vm∈H(Ωm,d)

(6)

(7)

式中f為宏觀結(jié)構(gòu)的體力，h為宏觀結(jié)構(gòu)的Neumann邊界ΓM上的牽引力。宏觀和微觀彈性張量DM和Dm可基于能量均勻化方法求解，根據(jù)修正的SIMP插值方法[4]定義如下，

(8)

式中為避免數(shù)值奇異，Emin為孔洞材料彈性模量，p為懲罰因子，根據(jù)相關(guān)文獻(xiàn)，當(dāng)p=3時(shí)能得到好的優(yōu)化結(jié)果。DH為利用能量均勻化方法求得的彈性張量，

(9)

4 靈敏度分析

為了便于數(shù)值計(jì)算，利用有限元技術(shù)對(duì)式(6)的目標(biāo)函數(shù)顯式化，根據(jù)對(duì)熱平衡方程和等效熱載荷的分析，將剛度優(yōu)化目標(biāo)函數(shù)改寫為

C=FTU=(Fm+Fε)TU

(10)

采用直接法[22]計(jì)算目標(biāo)函數(shù)對(duì)宏觀設(shè)計(jì)變量的敏度信息，同時(shí)結(jié)構(gòu)外載荷Fm是一個(gè)確定常數(shù)，與設(shè)計(jì)變量無(wú)關(guān)，

(11)

以2D情形為例，?Ft h/?xM的求解具體如下,

[-1-11-111-11]T

(12)

對(duì)于離散單元的溫度變化量處理如下，

(13)

式中ti為單元節(jié)點(diǎn)溫度，n為離散單元的節(jié)點(diǎn)數(shù)，t0為參考溫度，通常情況下取值為0。根據(jù)文獻(xiàn)[23]描述單元節(jié)點(diǎn)溫度ti表達(dá)如下，

(14)

(15)

此時(shí)可以得出離散單元溫度變化對(duì)設(shè)計(jì)變量的導(dǎo)數(shù)為

(16)

將式(12,16)代入式(11)整理得到2D情形的目標(biāo)函數(shù)靈敏度,

(17)

式中ke為單元的剛度矩陣，ue為單元的位移向量。根據(jù)式(6)求解宏觀設(shè)計(jì)變量的體積靈敏度,

(18)

對(duì)于微結(jié)構(gòu)設(shè)計(jì)變量的目標(biāo)靈敏度和體積靈敏度求解具體如下，

(19)

對(duì)于?DH(xm)/?xm的求解如下,

(20)

5 過(guò)濾技術(shù)及優(yōu)化算法

為了獲得清晰0/1離散分布的拓?fù)鋱D并抑制宏微觀尺度在優(yōu)化過(guò)程中的數(shù)值不穩(wěn)定現(xiàn)象，根據(jù)文獻(xiàn)[24]采用Heaviside非線性密度過(guò)濾技術(shù)，修正后的單元密度為

(21)

式中引入一個(gè)正參數(shù)β來(lái)平滑Heaviside函數(shù)，同時(shí)

(22)

式中Ne為與單元e的質(zhì)心距離小于過(guò)濾半徑所有單元f的數(shù)目，權(quán)重系數(shù)He f具體如下,

He f=max.{0,rmin-dist(e,f)}

(23)

式中rmin為過(guò)濾半徑，dist(e,f)為單元e和單元f的質(zhì)心距離。

6 算例分析

通過(guò)對(duì)熱固耦合并行化拓?fù)鋬?yōu)化算例進(jìn)行分析，驗(yàn)證本文方法拓?fù)鋬?yōu)化的有效性。根據(jù)文獻(xiàn)[9]可知，通常初始設(shè)計(jì)域由均勻分布的密度場(chǎng)組成，以避免出現(xiàn)局部最小值的設(shè)計(jì)，然而由于應(yīng)用的周期邊界條件將形成均勻分布的敏度場(chǎng)，不能持續(xù)更新變量，因此需要在微結(jié)構(gòu)的初始設(shè)計(jì)域定義一些簡(jiǎn)單的空白區(qū)域，微結(jié)構(gòu)的初始設(shè)計(jì)結(jié)構(gòu)形狀如圖1所示。宏觀結(jié)構(gòu)的初始設(shè)計(jì)域長(zhǎng)度和寬度分別設(shè)為40 cm和15 cm，離散為80×30的網(wǎng)格，設(shè)計(jì)域頂部施加均勻分布的發(fā)熱熱源，底邊中點(diǎn)加載方向向下的集中載荷，左右兩邊界固定且絕熱，同時(shí)底邊溫度設(shè)置為0 ℃，具體載荷和邊界條件如圖2所示。微結(jié)構(gòu)劃分為50×50的網(wǎng)格，每一個(gè)微結(jié)構(gòu)的尺寸均為5 mm，微結(jié)構(gòu)的尺寸遠(yuǎn)小于宏觀結(jié)構(gòu)的尺寸，符合均勻化理論的要求。過(guò)濾半徑設(shè)置為1.4，其余設(shè)計(jì)參數(shù)列入表1。

該算例的停止準(zhǔn)則定義為優(yōu)化迭代步數(shù)不小于100步且相鄰迭代的宏微觀設(shè)計(jì)變量最大差值的絕對(duì)值均不超過(guò)0.05。

圖1 微結(jié)構(gòu)初始設(shè)計(jì)形狀

圖2 結(jié)構(gòu)載荷及約束條件

圖3給出并行化拓?fù)鋬?yōu)化結(jié)構(gòu)和微結(jié)構(gòu)的彈性張量矩陣和熱傳導(dǎo)系數(shù)矩陣；圖4為非并行化和并行化拓?fù)鋬?yōu)化的優(yōu)化迭代圖，由于微觀拓?fù)湓趦?yōu)化過(guò)程計(jì)算得到微結(jié)構(gòu)的等效材料屬性的精度問(wèn)題導(dǎo)致并行化優(yōu)化過(guò)程中出現(xiàn)幾次波動(dòng)，但不影響最終的穩(wěn)定性。為了驗(yàn)證所提出的優(yōu)化方法能夠有效地提高結(jié)構(gòu)的剛度性能和散熱性能，給出并行化拓?fù)鋬?yōu)化和非并行化拓?fù)鋬?yōu)化對(duì)于宏觀結(jié)構(gòu)在x和y方向的位移以及結(jié)構(gòu)溫度分布的對(duì)比云圖，如圖5所示，左側(cè)為并行化拓?fù)鋬?yōu)化的結(jié)果，右側(cè)為非并行化的情況，可知非并行化情況下結(jié)構(gòu)的最高溫度大約在4 ℃～5 ℃，沿x方向的最大位移在0.1 cm～0.15 cm，沿y方向的最大位移大于 1 cm；然而經(jīng)過(guò)并行化處理之后的熱固耦合結(jié)構(gòu)的最高溫度為0.3 ℃，在x方向位移大約2.5×10-3cm～3×10-3cm，沿y方向最大位移約為 0.025 cm，數(shù)據(jù)結(jié)果表明，并行化處理之后，結(jié)構(gòu)的溫度分布和位移形變得到大幅度的降低，由此可知所提出的優(yōu)化方法能夠顯著改善結(jié)構(gòu)剛度性能和散熱性能。

表1 優(yōu)化參數(shù)

圖3 2D并行化拓?fù)鋬?yōu)化結(jié)果

圖4 拓?fù)鋬?yōu)化迭代

在2D的基礎(chǔ)上開展3D情形的拓?fù)鋬?yōu)化設(shè)計(jì)，三維懸臂梁的載荷、約束條件以及微結(jié)構(gòu)的初始設(shè)計(jì)結(jié)構(gòu)如圖6所示，結(jié)構(gòu)場(chǎng)的右邊界固定約束，左邊界底邊施加向下的力，溫度場(chǎng)頂部施加均勻分布熱源，底部溫度設(shè)置為0 ℃。宏觀結(jié)構(gòu)的尺寸為30 cm×10 cm×1 cm，劃分為60×20×2單元；微結(jié)構(gòu)尺寸為5 mm×5 mm×5 mm，劃分為 10×10×10單元。為了能夠進(jìn)行微觀尺度的設(shè)計(jì)變量更新，將微結(jié)構(gòu)的內(nèi)部中心刪除2×2×2個(gè)單元。算例過(guò)濾半徑設(shè)置為1.5，其余優(yōu)化設(shè)計(jì)參數(shù)列入表2。

該算例的停止準(zhǔn)則定義為優(yōu)化迭代步數(shù)不小于100步，且相鄰迭代的宏微觀設(shè)計(jì)變量最大差值的絕對(duì)值均不超過(guò)0.02。

圖7展示了3D情形的并行化拓?fù)浣Y(jié)構(gòu)，包括微結(jié)構(gòu)的彈性張量矩陣以及熱傳導(dǎo)系數(shù)矩陣。圖8 是并行化優(yōu)化過(guò)程迭代圖。由于并行化的緣故，宏觀結(jié)構(gòu)的拓?fù)湫螤顚?duì)比于非并行化的拓?fù)浣Y(jié)構(gòu)有很大改變，如圖9所示。微結(jié)構(gòu)整體上是有兩塊直板構(gòu)成，中間沒(méi)有材料連接的結(jié)構(gòu)，但是如果將材料的彈性模量和熱傳導(dǎo)系數(shù)增大，如彈性模量設(shè)置為10 GPa，熱傳導(dǎo)系數(shù)設(shè)置為10 W/(m·℃)，其余優(yōu)化參數(shù)不變，具體優(yōu)化結(jié)果如圖10所示。對(duì)比圖7的拓?fù)浣Y(jié)構(gòu)，宏觀結(jié)構(gòu)和微結(jié)構(gòu)有很大的變化，由于提高了材料的力學(xué)屬性，宏觀拓?fù)浣Y(jié)構(gòu)明顯少了一些分支，同時(shí)微結(jié)構(gòu)的中間部分有材料連接，因此改變材料參數(shù)會(huì)使得結(jié)構(gòu)的拓?fù)湫螤畎l(fā)生較大的變化，側(cè)面說(shuō)明了所提出的優(yōu)化設(shè)計(jì)方法能夠有效地進(jìn)行優(yōu)化設(shè)計(jì)。

圖5 并行化拓?fù)鋬?yōu)化與非并行化的結(jié)果對(duì)比

表2 優(yōu)化參數(shù)

圖6 懸臂梁載荷和約束條件以及微結(jié)構(gòu)初始結(jié)構(gòu)

圖7 3D并行化拓?fù)鋬?yōu)化結(jié)果

圖8 3D并行化優(yōu)化迭代

圖9 非并行化宏觀拓?fù)浣Y(jié)構(gòu)

圖10 并行化拓?fù)鋬?yōu)化結(jié)果

7 總 結(jié)

本文利用能量均勻化方法結(jié)合拓?fù)鋬?yōu)化技術(shù)建立了基于穩(wěn)態(tài)熱源作用下的熱固耦合結(jié)構(gòu)多尺度拓?fù)鋬?yōu)化設(shè)計(jì)模型，以結(jié)構(gòu)剛度為優(yōu)化目標(biāo)，體分比為約束，通過(guò)開展2D和3D的典型數(shù)值算例優(yōu)化設(shè)計(jì)，驗(yàn)證了所提出的優(yōu)化方法能夠有效地進(jìn)行優(yōu)化設(shè)計(jì)，并且通過(guò)與非并行化熱固耦合拓?fù)鋬?yōu)化做對(duì)比，證明了該優(yōu)化方法能夠顯著地改善結(jié)構(gòu)的剛度性能和散熱性能，順應(yīng)了現(xiàn)代化生產(chǎn)對(duì)工業(yè)產(chǎn)品高強(qiáng)度輕量化的設(shè)計(jì)趨勢(shì)，為基于多物理場(chǎng)耦合問(wèn)題下開展結(jié)構(gòu)/材料多尺度一體化拓?fù)鋬?yōu)化設(shè)計(jì)提供了一定的參考價(jià)值和指導(dǎo)意義。