李 冉 劉書田

(大連理工大學(xué)工程力學(xué)系,工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室,遼寧大連 116024)

引言

拓?fù)鋬?yōu)化[1]是一種基于給定載荷和約束,在特定的域內(nèi)尋找結(jié)構(gòu)最優(yōu)傳力路徑及材料分布的設(shè)計(jì)方法,具有不依賴于初始構(gòu)型、設(shè)計(jì)自由度高及性能驅(qū)動(dòng)的特點(diǎn),已成為結(jié)構(gòu)與材料創(chuàng)新設(shè)計(jì)的重要工具[2-5].但由于幾何構(gòu)型較復(fù)雜,采用傳統(tǒng)的制造工藝往往加工困難.增材制造技術(shù)[6-9]是一種通過逐點(diǎn)打印、逐層累積材料來制備結(jié)構(gòu)的新興技術(shù),可實(shí)現(xiàn)高度復(fù)雜結(jié)構(gòu)的制備.將拓?fù)鋬?yōu)化方法與增材制造技術(shù)相融合,在實(shí)現(xiàn)拓?fù)鋬?yōu)化設(shè)計(jì)價(jià)值最大化的同時(shí),還可以充分發(fā)揮增材制造的制造潛力,是目前研究的重要趨勢[10-11].

但由于制造方式、加工過程及環(huán)境、設(shè)備和操作等因素的影響,通過增材制造成型的結(jié)構(gòu)件易出現(xiàn)制造缺陷[12-14].制造缺陷使得結(jié)構(gòu)件成型精度降低、表面粗糙度大,進(jìn)而導(dǎo)致結(jié)構(gòu)件表面材料分布不均勻,引起表面層異質(zhì)[15].表面層異質(zhì)會(huì)帶來結(jié)構(gòu)表面材料屬性的不確定性和表面層厚度的不確定性,致使結(jié)構(gòu)性能偏離設(shè)計(jì)、達(dá)不到預(yù)期目標(biāo).然而,當(dāng)前大多數(shù)拓?fù)鋬?yōu)化方法的研究是基于確定性的參數(shù)開展,忽略了不確定因素對結(jié)構(gòu)性能的影響,導(dǎo)致優(yōu)化后的構(gòu)型不一定是最優(yōu)的或者不符合實(shí)際,嚴(yán)重影響結(jié)構(gòu)的安全性.因此在拓?fù)鋬?yōu)化設(shè)計(jì)中,考慮不確定因素對結(jié)構(gòu)構(gòu)型的影響是十分重要的.

一般來說,考慮不確定性的拓?fù)鋬?yōu)化模型包括基于可靠度的拓?fù)鋬?yōu)化[16-17]和穩(wěn)健性拓?fù)鋬?yōu)化.其中,基于可靠度的拓?fù)鋬?yōu)化關(guān)心的是結(jié)構(gòu)在極端條件下的可靠性,而穩(wěn)健性拓?fù)鋬?yōu)化則關(guān)心的是結(jié)構(gòu)性能在其均值處的變異程度,旨在降低結(jié)構(gòu)性能對不確定因素的敏感程度.近年來,一批學(xué)者開展了針對穩(wěn)健性拓?fù)鋬?yōu)化問題的研究.Asadpoure 等[18]針對桁架結(jié)構(gòu)的隨機(jī)問題,采用攝動(dòng)技術(shù)對桁架剛度的不確定性進(jìn)行了建模和隨機(jī)響應(yīng)分析.Sigmund[19]基于圖像形態(tài)學(xué)中的腐蝕、擴(kuò)散操作[20]對制備過程中均勻變化的過刻蝕、欠刻蝕進(jìn)行了模擬,并在光子晶體波導(dǎo)的穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)中得以驗(yàn)證[21].隨后,Schevenels 等[22]通過把離散映射的控制閾值作為隨機(jī)變量,將該研究擴(kuò)展到考慮非均勻制造誤差的情況.Kang 等[23]基于水平集方法跟蹤多材料結(jié)構(gòu)界面的演化,利用級數(shù)最優(yōu)線性估計(jì)(EOLE)結(jié)合多項(xiàng)式混沌(PCE)技術(shù)對材料界面寬度的不確定性進(jìn)行了建模、離散及隨機(jī)響應(yīng)分析,建立了考慮多材料結(jié)構(gòu)界面不確定性的穩(wěn)健性形狀和拓?fù)鋬?yōu)化方法.針對具有隨機(jī)和區(qū)間混合不確定的均勻周期性微結(jié)構(gòu)動(dòng)力及結(jié)構(gòu)并行設(shè)計(jì)問題,Zheng 等[24]提出了一種基于水平集的穩(wěn)健拓?fù)鋬?yōu)化方法,采用基于二元降維方案的混合降維方法來估計(jì)不確定目標(biāo)函數(shù)的區(qū)間均值和區(qū)間方差.

研究結(jié)構(gòu)表面層異質(zhì)引起的表面層厚度不確定性,需要解決兩個(gè)關(guān)鍵問題:(1) 準(zhǔn)確識別結(jié)構(gòu)表面層;(2)對不確定性進(jìn)行傳播分析和隨機(jī)響應(yīng)估計(jì).本文建立了考慮表面層厚度不確定性的拓?fù)鋬?yōu)化模型,采用基于腐蝕操作的表面層識別技術(shù),實(shí)現(xiàn)了表面層異質(zhì)等效模型的建立.基于攝動(dòng)有限元法開展了不確定性傳播的分析和系統(tǒng)隨機(jī)響應(yīng)的預(yù)測,并推導(dǎo)了目標(biāo)函數(shù)關(guān)于設(shè)計(jì)變量的敏度信息.通過數(shù)值算例證實(shí)了考慮表面層厚度不確定性的必要性及本方法的有效性.

1 基于腐蝕操作的結(jié)構(gòu)表面層識別技術(shù)

如何準(zhǔn)確地對結(jié)構(gòu)表面層進(jìn)行建模,通常是基于密度的拓?fù)鋬?yōu)化方法中較為困難的問題.Clausen等[25]采用基于兩步過濾的操作,以較大的空間密度梯度識別了結(jié)構(gòu)的表面層,并成功應(yīng)用于基于密度的拓?fù)鋬?yōu)化問題[26-27].本文采用一種基于腐蝕操作的表面層識別技術(shù)[28]建立表面層異質(zhì)的等效模型.其主要思想是:首先對圖1(a)所示的初始結(jié)構(gòu)進(jìn)行腐蝕操作,得到圖1(b) 所示規(guī)?？s減后的被腐蝕結(jié)構(gòu);然后將初始結(jié)構(gòu)與被腐蝕結(jié)構(gòu)的不同部分,即圖1(c)所示的被腐蝕掉的部分,定義為結(jié)構(gòu)表面層.

圖1 基于腐蝕操作的表面層識別方法的基本思想Fig.1 Basic idea of the erosion-based interface identification method

腐蝕操作通過光滑過濾和離散映射兩個(gè)過程來實(shí)現(xiàn).其中,光滑過濾使得初始結(jié)構(gòu)場光滑連續(xù);離散映射使用一個(gè)較大的閾值,將光滑化的結(jié)構(gòu)場進(jìn)行截?cái)?得到被腐蝕后的結(jié)構(gòu),進(jìn)而識別出結(jié)構(gòu)表面層.圖2 給出了一維情況下腐蝕操作的一般過程(文中腐蝕操作的所有參數(shù)都用下標(biāo)e 標(biāo)記).

圖2 一維情況下腐蝕操作的一般過程Fig.2 1D example of the general process of the erosion operation

本文中,首先使用基于Helmholtz 偏微分方程的密度過濾[29-30]——PDE 技術(shù)[31],對初始結(jié)構(gòu)場μ進(jìn)行光滑化處理

式中βe為映射函數(shù)的銳度,閾值ηe=0.5+Δηe通過引入Δηe∈(0,0.5)來定義.

基于上述光滑過濾和離散映射的腐蝕操作,可以將厚度均勻可控的結(jié)構(gòu)表面層γ 識別出來.表面層厚度w與過濾半徑Re、閾值Δηe之間滿足如下關(guān)系

式中,通過合理選取Re和Δηe的值即可實(shí)現(xiàn)對表面層厚度w的精確控制.本文采取的策略是固定Δηe的值,通過改變Re的值來控制w.為了得到w和Re更規(guī)整的近似關(guān)系,文中算例依據(jù)經(jīng)驗(yàn)Δηe均取值為0.45,此時(shí)Re≈1.50w.(3)式的詳細(xì)理論推導(dǎo)參見文獻(xiàn)[28].在此基礎(chǔ)上,物理密度和剛度的插值函數(shù)可以定義為

式中ρb,Eb,ρs,Es分別表示基體和表面層的材料性質(zhì),p為懲罰因子,一般取值為3.當(dāng)μ=1 且φ=1時(shí),描述結(jié)構(gòu)的基體部分;當(dāng)μ=1 且φ=0 時(shí),描述結(jié)構(gòu)表面層;當(dāng)μ=0 時(shí),表示結(jié)構(gòu)為空.

2 基于攝動(dòng)有限元的穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)

2.1 攝動(dòng)有限元法及系統(tǒng)隨機(jī)響應(yīng)估計(jì)

隨機(jī)有限元法是將隨機(jī)分析理論與有限元方法相結(jié)合的數(shù)值分析方法.其中,基于攝動(dòng)開展的隨機(jī)有限元法得到了廣泛的應(yīng)用[18,34-35].其基本思想是:引入一個(gè)零均值的微小攝動(dòng)量Δξ,來描述隨機(jī)變量ξ 在其均值ξ0處產(chǎn)生微小攝動(dòng)所引起的不確定性,那么隨機(jī)變量ξ 可以表示為

若假設(shè)系統(tǒng)發(fā)生準(zhǔn)靜態(tài)線彈性行為,且不確定性僅存在于某些結(jié)構(gòu)參數(shù)(如幾何形狀、材料性質(zhì))中,而外荷載是確定的,那么有限元平衡方程可以寫作

式中K(ξ)和u(ξ)分別表示與隨機(jī)變量ξ 相關(guān)的剛度矩陣和位移響應(yīng)向量.在ξ0處對K(ξ)和u(ξ)進(jìn)行二階Taylor 展開

式中K0=K(ξ0),u0=u(ξ0)且Ki,Ki j,ui,uij為K和u在ξ0處對隨機(jī)變量ξ 的一階、二階偏導(dǎo)數(shù).將式(8)、式(9)代入式(7),忽略Δξ 的高階小量同時(shí)合并同階項(xiàng),可以得到

通過按自上而下的順序求解方程組(10),可以得到由式(9)表示的位移響應(yīng)u的二階近似表達(dá).若Δξ服從正態(tài)分布,那么位移響應(yīng)u的均值及協(xié)方差的二階近似估計(jì)表示為

式中σij=E[ΔξiΔξj]表示Δξ 的二階統(tǒng)計(jì)矩.

2.2 考慮表面層厚度不確定的穩(wěn)健性拓?fù)鋬?yōu)化方法

本文考慮靜力學(xué)拓?fù)鋬?yōu)化中的最小柔順性設(shè)計(jì)問題,即在確定的外載荷下,給定有限體積的材料,設(shè)計(jì)得到儲(chǔ)存應(yīng)變能最小的結(jié)構(gòu)拓?fù)?優(yōu)化列式可以描述為

式中,設(shè)計(jì)變量ρe是表示每個(gè)單元的密度,c表示結(jié)構(gòu)的柔順性,K,u,f分別表示結(jié)構(gòu)的總體剛度矩陣、位移向量和外載荷向量,V和ˉV分別表示結(jié)構(gòu)所用材料的體積和給定的體積上限.

本文第1 節(jié)中,通過基于腐蝕操作的表面層識別技術(shù),對厚度為w的結(jié)構(gòu)表面層進(jìn)行了建模.本節(jié)中以表面層厚度w為隨機(jī)變量,根據(jù)攝動(dòng)有限元法的思想,引入一個(gè)均值為、標(biāo)準(zhǔn)差為σ 且服從高斯分布的微小攝動(dòng)量δ,對表面層厚度的不確定性進(jìn)行建模.穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)旨在優(yōu)化結(jié)構(gòu)的性能,并且最小化其對不確定因素的敏感程度.那么,考慮表面層厚度不確定性的穩(wěn)健性拓?fù)鋬?yōu)化列式可以描述為

那么,結(jié)構(gòu)柔順性均值E[c]、方差Var[c]的近似估計(jì)可以表示為

3 數(shù)值算例

為了避免棋盤效應(yīng)和網(wǎng)格依賴性對優(yōu)化結(jié)果的影響,并保證得到0-1 分布的結(jié)構(gòu)場,在識別結(jié)構(gòu)表面層之前,先后對初始結(jié)構(gòu)場應(yīng)用了PDE 過濾、基于Heaviside 過濾和tanh 函數(shù)的離散映射.這一過程與基于腐蝕操作的表面層識別幾乎相同,但卻實(shí)現(xiàn)了不同的功能.文中將這一過程定義為過濾過程(所有參數(shù)都用下標(biāo)f標(biāo)記),將表面層識別的過程定義為腐蝕過程,并將過濾過程和腐蝕過程統(tǒng)稱為改進(jìn)的兩步過濾法,其實(shí)現(xiàn)流程及控制參數(shù)見圖3.此外,所有算例都使用魯棒公式[19,21-22]來實(shí)現(xiàn)結(jié)構(gòu)最小尺寸的控制,控制參數(shù)[28]為Δηf;同時(shí)通過域擴(kuò)展方法[37]保證在結(jié)構(gòu)邊界處獲得相同的腐蝕厚度.

圖3 改進(jìn)的兩步過濾法識別結(jié)構(gòu)表面層Fig.3 Improved two-step filter process for structure surface layer identification

本文算例中數(shù)值為常量的控制參數(shù)如下:懲罰系數(shù)p=3;當(dāng)計(jì)算收斂時(shí)或每50 個(gè)迭代步后,投影銳度βe=βf∈(1,20)變?yōu)楫?dāng)前值的1.5 倍.

3.1 簡支梁算例

圖4 所示為簡支梁算例的設(shè)計(jì)域及邊界條件,結(jié)構(gòu)的幾何尺寸為300 × 150,集中力F=1.表面層結(jié)構(gòu)和基體結(jié)構(gòu)的材料屬性分別為ρs=1,Es=1,ρb=0.8,Eb=0.5,即表面層為硬材料,基體為軟材料.表面層厚度為w=2,其標(biāo)準(zhǔn)差σ=0.3w.過濾過程和腐蝕過程的過濾半徑分別取值為Rf=8,Re≈1.5w,最小尺寸的控制參數(shù)Δηf=0.1.體積分?jǐn)?shù)為40%且每10 個(gè)迭代步進(jìn)行一次更新[28].目標(biāo)函數(shù)中的權(quán)系數(shù)α 選取值為4.

圖4 簡支梁算例的設(shè)計(jì)域及邊界條件Fig.4 Design domain and boundary conditions of MBB example

如圖5 所示,簡支梁結(jié)構(gòu)的確定性和穩(wěn)健性拓?fù)鋬?yōu)化結(jié)果有著不同的拓?fù)錁?gòu)型,相比于確定性設(shè)計(jì),穩(wěn)健性設(shè)計(jì)去除了細(xì)小孔洞的存在.不同結(jié)果對應(yīng)的結(jié)構(gòu)柔順性均值和標(biāo)準(zhǔn)差數(shù)值,列入表1.與確定性設(shè)計(jì)相比(σ=3.952 8),穩(wěn)健性設(shè)計(jì)(σ=3.383 3)的標(biāo)準(zhǔn)差較小,這表明穩(wěn)健性拓?fù)鋬?yōu)化結(jié)果對表面層厚度變化的敏感性更低.圖6 給出了確定性設(shè)計(jì)和穩(wěn)健性設(shè)計(jì)在優(yōu)化過程中目標(biāo)函數(shù)和體積分?jǐn)?shù)的迭代歷史,最后都達(dá)到了穩(wěn)定收斂.

圖5 簡支梁算例的拓?fù)鋬?yōu)化結(jié)果Fig.5 Topology optimization results of MBB example

表1 確定性設(shè)計(jì)和穩(wěn)健性設(shè)計(jì)結(jié)果的對比Table 1 Comparision of deterministic and rubost solutions

圖6 簡支梁算例目標(biāo)函數(shù)和體積分?jǐn)?shù)的迭代歷史Fig.6 Convergence history of objective function and volume fraction of MBB example

為了研究目標(biāo)函數(shù)中權(quán)系數(shù)α 的取值對拓?fù)鋬?yōu)化結(jié)果的影響,本文比較了當(dāng)α 分別取值為2、4、6時(shí),表1 中列出了穩(wěn)健性設(shè)計(jì)下結(jié)構(gòu)柔順性的均值和標(biāo)準(zhǔn)差的數(shù)值,圖7 所示為相對應(yīng)α 取不同值時(shí)結(jié)構(gòu)柔順性的均值與標(biāo)準(zhǔn)差的變化趨勢.可以看出結(jié)構(gòu)柔順性的均值隨著權(quán)系數(shù)α 的增大而增大,但其標(biāo)準(zhǔn)差卻呈現(xiàn)出相反的變化趨勢,這與雙準(zhǔn)則優(yōu)化問題帕累托解的基本特征是吻合的.同時(shí),標(biāo)準(zhǔn)差越小反映出結(jié)構(gòu)對不確定因素的敏感程度越低.

圖7 權(quán)系數(shù)α 取不同值時(shí)結(jié)構(gòu)柔順性的均值與標(biāo)準(zhǔn)差Fig.7 Mean and standard deviations of compliance for different weighing factor α

3.2 懸臂梁算例

圖8 所示為懸臂梁算例的設(shè)計(jì)域及邊界條件,結(jié)構(gòu)的幾何尺寸為300×100,集中力F=1.表面層結(jié)構(gòu)和基體結(jié)構(gòu)的材料屬性分別為ρs=0.7,Es=0.4,ρb=1,Eb=1,即表面層為軟材料,基體為硬材料.表面層厚度為w=2,其標(biāo)準(zhǔn)差σ=0.05w.過濾過程和腐蝕過程的過濾半徑分別取值為Rf=8,Re≈1.5w,最小尺寸的控制參數(shù)Δηf=0.1.體積分?jǐn)?shù)為40%且每10 個(gè)迭代步進(jìn)行一次更新.目標(biāo)函數(shù)中的權(quán)系數(shù)α 選取為1.

圖8 懸臂梁算例的設(shè)計(jì)域及邊界條件Fig.8 Design domain and boundary conditions of cantilever beam example

結(jié)構(gòu)的確定性和穩(wěn)健性拓?fù)鋬?yōu)化結(jié)果如圖9 所示,相比于確定性結(jié)果,穩(wěn)健性結(jié)果去除了細(xì)小桿件的存在,且結(jié)構(gòu)表面更為平滑.表2 所示為σ 取不同值時(shí)確定性設(shè)計(jì)和穩(wěn)健性設(shè)計(jì)結(jié)果對比.從表2 可以看出,穩(wěn)健性設(shè)計(jì)的標(biāo)準(zhǔn)差較小,其均值和標(biāo)準(zhǔn)差的加權(quán)和值也更低,再次驗(yàn)證了穩(wěn)健性結(jié)果在結(jié)構(gòu)整體剛度較高的同時(shí)具有更好的抵抗不確定性的能力.優(yōu)化過程中目標(biāo)函數(shù)和體積分?jǐn)?shù)的迭代歷史如圖10 所示,都達(dá)到了穩(wěn)定收斂.

圖9 懸臂梁算例的拓?fù)鋬?yōu)化結(jié)果Fig.9 Topology optimization results of cantilever beam example

圖10 懸臂梁算例目標(biāo)函數(shù)和體積分?jǐn)?shù)的迭代歷史Fig.10 Convergence history of objective function and volume fraction of cantilever beam example

本文研究了當(dāng)表面層厚度標(biāo)準(zhǔn)差取不同值時(shí),對穩(wěn)健性拓?fù)鋬?yōu)化結(jié)果的影響.在其他參數(shù)保持不變的情況下,σ=0.05w,σ=0.2w,σ=0.3w所對應(yīng)的穩(wěn)健性拓?fù)鋬?yōu)化結(jié)果分別如圖9(b)、圖11(a)、圖11(b)所示.可以看出,當(dāng)不確定性相對較小(σ=0.05w)時(shí),穩(wěn)健性設(shè)計(jì)與確定性設(shè)計(jì)有著相似的結(jié)構(gòu);而在σ=0.2w,σ=0.3w的情況下,穩(wěn)健性結(jié)果與確定性結(jié)果在桿件數(shù)量、布局等方面具有一定的差異.3 種情況下結(jié)構(gòu)柔順性的均值和標(biāo)準(zhǔn)差的數(shù)據(jù)比較見表2.數(shù)據(jù)表明,任一情況下穩(wěn)健性設(shè)計(jì)結(jié)果的標(biāo)準(zhǔn)差值、均值和標(biāo)準(zhǔn)差的和值都小于確定性結(jié)果,進(jìn)一步驗(yàn)證了本文所提出方法的有效性.

表2 σ 取不同值時(shí)確定性設(shè)計(jì)和穩(wěn)健性設(shè)計(jì)結(jié)果的對比Table 2 Comparision of deterministic and rubost solutions for different σ

圖11 σ 取不同值時(shí)懸臂梁算例的穩(wěn)健性拓?fù)鋬?yōu)化結(jié)果Fig.11 Robust topology optimization results of cantilever beam example for different σ

4 結(jié)論

針對增材制造工藝在制備過程中出現(xiàn)的結(jié)構(gòu)表面層異質(zhì)現(xiàn)象,本文以結(jié)構(gòu)柔順性均值和標(biāo)準(zhǔn)差的加權(quán)和作為優(yōu)化目標(biāo),建立了考慮結(jié)構(gòu)表面層厚度不確定的穩(wěn)健性拓?fù)鋬?yōu)化模型.采用一種基于腐蝕操作的表面層識別技術(shù)對表面層異質(zhì)進(jìn)行了等效建模,并將表面層厚度視為具有給定概率分布特征的隨機(jī)變量,通過隨機(jī)攝動(dòng)理論開展了不確定性傳播分析和系統(tǒng)隨機(jī)響應(yīng)預(yù)測.數(shù)值算例表明,所提出的方法能夠得到對結(jié)構(gòu)表面層厚度變化不敏感的合理結(jié)構(gòu),具有更強(qiáng)的抵抗不確定性的能力.因此,對增材制造結(jié)構(gòu)而言,在拓?fù)鋬?yōu)化設(shè)計(jì)中考慮表面層厚度不確定性對結(jié)構(gòu)性能的影響具有重要意義.