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

魏兆棟， 高仁璟,2， 王長生*,2

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

1 引 言

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

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

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

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

(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)所承受的機械外載，F(xiàn)ε為由于熱應(yīng)變產(chǎn)生的熱載荷向量。

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

ε0=αΔTΦT

(2)

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

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

(3)

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

(4)

以2D問題為例，對式(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è)計模型。將OC準(zhǔn)則用于設(shè)計變量的迭代更新，即分別在每個尺度中進行設(shè)計變量的優(yōu)化更新，利用能量均勻化計算每一次迭代后微結(jié)構(gòu)拓?fù)涞牡刃傩裕蛊渥鳛閮蓚€尺度之間的連接橋梁，其對應(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)文獻，當(dāng)p=3時能得到好的優(yōu)化結(jié)果。DH為利用能量均勻化方法求得的彈性張量，

(9)

4 靈敏度分析

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

C=FTU=(Fm+Fε)TU

(10)

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

(11)

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

[-1-11-111-11]T

(12)

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

(13)

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

(14)

(15)

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

(16)

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

(17)

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

(18)

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

(19)

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

(20)

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

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

(21)

式中引入一個正參數(shù)β來平滑Heaviside函數(shù)，同時

(22)

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

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

(23)

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

6 算例分析

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

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

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

圖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)化過程計算得到微結(jié)構(gòu)的等效材料屬性的精度問題導(dǎo)致并行化優(yōu)化過程中出現(xiàn)幾次波動，但不影響最終的穩(wěn)定性。為了驗證所提出的優(yōu)化方法能夠有效地提高結(jié)構(gòu)的剛度性能和散熱性能，給出并行化拓?fù)鋬?yōu)化和非并行化拓?fù)鋬?yōu)化對于宏觀結(jié)構(gòu)在x和y方向的位移以及結(jié)構(gòu)溫度分布的對比云圖，如圖5所示，左側(cè)為并行化拓?fù)鋬?yōu)化的結(jié)果，右側(cè)為非并行化的情況，可知非并行化情況下結(jié)構(gòu)的最高溫度大約在4 ℃～5 ℃，沿x方向的最大位移在0.1 cm～0.15 cm，沿y方向的最大位移大于 1 cm；然而經(jīng)過并行化處理之后的熱固耦合結(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è)計，三維懸臂梁的載荷、約束條件以及微結(jié)構(gòu)的初始設(shè)計結(jié)構(gòu)如圖6所示，結(jié)構(gòu)場的右邊界固定約束，左邊界底邊施加向下的力，溫度場頂部施加均勻分布熱源，底部溫度設(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單元。為了能夠進行微觀尺度的設(shè)計變量更新，將微結(jié)構(gòu)的內(nèi)部中心刪除2×2×2個單元。算例過濾半徑設(shè)置為1.5，其余優(yōu)化設(shè)計參數(shù)列入表2。

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

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

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

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