吉先磊，張美艷，王　皓

(復(fù)旦大學(xué) 航空航天系，上海 200433)

?

基于影響區(qū)域節(jié)點密度插值的結(jié)構(gòu)拓?fù)鋬?yōu)化方法

吉先磊，張美艷，王皓

(復(fù)旦大學(xué) 航空航天系，上海 200433)

摘要：變密度法是結(jié)構(gòu)拓?fù)鋬?yōu)化中一種常用的處理方法，但在優(yōu)化過程中會出現(xiàn)棋盤格、網(wǎng)格依賴性、灰度單元等數(shù)值不穩(wěn)定現(xiàn)象.為解決單元變密度法中出現(xiàn)的棋盤格問題，在優(yōu)化過程中引入節(jié)點密度作為設(shè)計變量.在提出基于影響區(qū)域節(jié)點密度插值優(yōu)化方法的同時，給出了一種基于等密度線的后處理方法.計算結(jié)果表明該方法在消除棋盤格現(xiàn)象的同時，處理了網(wǎng)格依賴性問題，并消除了結(jié)果中的灰度單元現(xiàn)象.

關(guān)鍵詞：拓?fù)鋬?yōu)化； 節(jié)點密度； 區(qū)域插值

連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化是在一個特定的設(shè)計區(qū)域內(nèi)，在邊界、載荷等約束條件已知的條件下，為達(dá)到結(jié)構(gòu)剛度最大、基頻最大等目標(biāo)而尋求材料最優(yōu)分布的一種結(jié)構(gòu)設(shè)計方法.拓?fù)鋬?yōu)化設(shè)計方法是目前連續(xù)體結(jié)構(gòu)優(yōu)化設(shè)計研究的熱點，取得了迅速發(fā)展，成為結(jié)構(gòu)設(shè)計的一種有效方法.

在拓?fù)鋬?yōu)化眾多的求解方法中，變密度法是一種常見的處理方法.在變密度法中通常設(shè)定單元密度為設(shè)計變量，通過單元密度的迭代獲得最優(yōu)的拓?fù)浣Y(jié)構(gòu).然而由于有限元單元密度相互獨立，普遍存在棋盤格和網(wǎng)格依賴性等問題[1]，在單元法中需要借助額外的方法如靈敏度過濾[2]、周長約束[3]、局部梯度[4]等進(jìn)行處理.為了從根本上避免棋盤格問題的出現(xiàn)，可以將設(shè)計變量定義在節(jié)點位置，使得單元密度場本身具有連續(xù)性.基于節(jié)點變量處理棋盤格問題的思想可以追溯到Y(jié)oun等[5]提出的基于節(jié)點的密度再分配后處理方法，該方法在單元密度迭代之后，定義節(jié)點密度為環(huán)繞單元密度的平均值，再由節(jié)點密度插值出新的單元密度，通過保證相鄰單元密度的連續(xù)性避免棋盤格的出現(xiàn).之后Rahmatalla等[6]，Guest等[7]在此基礎(chǔ)上提出了將有限元網(wǎng)格節(jié)點密度作為設(shè)計變量的優(yōu)化方法.Kang等[8]，占金青等[9]，Yi等[10]等推廣了該方法，將設(shè)計變量定義在非網(wǎng)格節(jié)點位置進(jìn)行處理.龍凱等[11]進(jìn)一步提出了基于物質(zhì)點描述的優(yōu)化方法.節(jié)點密度法可以避免棋盤格的出現(xiàn)，然而目前對于節(jié)點密度下處理網(wǎng)格依賴性和灰度單元問題的研究還較為有限.

本文將設(shè)計變量定義在網(wǎng)格節(jié)點位置，通過影響區(qū)域節(jié)點密度插值的方法在消除棋盤格的同時，避免了網(wǎng)格依賴性的出現(xiàn).

1影響區(qū)域節(jié)點密度插值方法

變密度法拓?fù)鋬?yōu)化以有限元單元作為優(yōu)化的基本單位，通過單元偽密度的引入將0，1離散優(yōu)化問題轉(zhuǎn)換為[0，1]上的連續(xù)優(yōu)化問題.變密度法假定單元的彈性模量與單元偽密度之間的關(guān)系滿足

(1)

其中：ρe為單元的偽密度，ρe在[0，1]上連續(xù)變化，ρe為0時表示該單元物質(zhì)不存在，ρe為1時表示該單元為實體密度；E0為材料固體部分彈性模量；Ee為單元的等效的彈性模量；p為懲罰參數(shù)，文中取為3.

如圖1(a)所示(見第312頁)，在目前常用的節(jié)點密度法中，通常將單元的密度定義為該單元包含節(jié)點密度的插值函數(shù).然而隨著設(shè)計區(qū)域網(wǎng)格劃分的不同，往往會產(chǎn)生嚴(yán)重的網(wǎng)格依賴性問題，主要表現(xiàn)在隨著網(wǎng)格的細(xì)分，優(yōu)化結(jié)果中結(jié)構(gòu)最小尺寸會隨之變小.在優(yōu)化過程中為獲得更精確的結(jié)果往往需要細(xì)化網(wǎng)格，而實際應(yīng)用中又需要對結(jié)構(gòu)中的最小尺寸做出限制.為了達(dá)到這樣的目標(biāo)，這里引入影響區(qū)域內(nèi)節(jié)點密度插值的方法.如圖1(b)所示(見第312頁)，在節(jié)點密度作為設(shè)計變量的基礎(chǔ)上，定義單元e的偽密度值為該單元插值區(qū)域上節(jié)點偽密度的加權(quán)平均，即：

ρe=ωe，iρi，

(2)

其中：ρe為加權(quán)平均的單元偽密度；ρi為節(jié)點偽密度值，是優(yōu)化過程中的設(shè)計變量；ωe，i為對應(yīng)節(jié)點變量的權(quán)函數(shù)，且

(3)

其中：rmin為限制的結(jié)構(gòu)最小尺寸；dist(e，i)表示單元e中心到節(jié)點i的距離；Ωr為距e單元中心距離小于rmin的子區(qū)域.

同時為便于獲得最終0，1分布的結(jié)構(gòu)，在優(yōu)化的結(jié)果中引入“密度等值線”作為參考.如圖2(見第312頁)所示的區(qū)域，其4個端點位置密度分別為1.0，1.0，1.0，0.001，定義其0.5密度等值線為圖中虛線所示線段，定義大于等密度線的部分為實體結(jié)構(gòu)，小于等密度線的部分為孔洞結(jié)構(gòu)，幫助獲得最終的優(yōu)化結(jié)構(gòu).

2優(yōu)化模型及靈敏度

基于節(jié)點變密度法，以結(jié)構(gòu)剛度最大化為目標(biāo)，體積為約束，連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的模型可以表示為

(4)

其中： C為目標(biāo)結(jié)構(gòu)的總體應(yīng)變能；ue為單元位移向量；U為整體位移向量；F為外在載荷；K為整體剛度矩陣；n為有限元單元數(shù)目；f為目標(biāo)體積分?jǐn)?shù)；V0為設(shè)計區(qū)域總體積；ρe為單元偽密度；ρi為節(jié)點偽密度變量.

基于變密度法，由(1)式，目標(biāo)函數(shù)對有限單元e偽密度的靈敏度表示為

(5)

通過鏈?zhǔn)角髮?dǎo)法則，進(jìn)而可以表示出目標(biāo)函數(shù)對節(jié)點密度變量的靈敏度為

(6)

(7)

3數(shù)值算例

下面采用2維平面拓?fù)鋬?yōu)化算例對本文方法進(jìn)行討論.

3.1算例1

考慮如圖3所示的梁結(jié)構(gòu)剛度最大化問題，結(jié)構(gòu)設(shè)計區(qū)域的尺寸為60cm×10cm×0.1cm，載荷在上邊界中點位置，大小為F=3500N，材料彈性模量為E=210GPa，泊松比為v=0.3，目標(biāo)體積為整個設(shè)計區(qū)域體積的50%.

采用平面四邊形單元對設(shè)計區(qū)域進(jìn)行離散，網(wǎng)格數(shù)量為120×20，網(wǎng)格尺寸為0.5cm×0.5cm.分別采用單元節(jié)點密度插值、本文的影響區(qū)域節(jié)點密度插值方法進(jìn)行優(yōu)化，兩種節(jié)點密度方法下的拓?fù)鋬?yōu)化結(jié)果對比如圖4所示.

如圖4(a)所示，在采用單元節(jié)點密度插值時，得到結(jié)構(gòu)的最小尺寸與網(wǎng)格的粗細(xì)相關(guān)聯(lián)，不可做進(jìn)一步處理.圖4(b-1)、圖4(b-2)、圖4(b-3)為采用本文影響區(qū)域節(jié)點密度插值方法獲得的結(jié)果，通過不同插值區(qū)域的選取，在網(wǎng)格固定的條件下，可以明顯改變結(jié)果中的最小結(jié)構(gòu)尺寸，從而避免網(wǎng)格依賴性對優(yōu)化結(jié)果的影響.

由于節(jié)點密度定義的變量值為節(jié)點密度值，這里將圖4中采用節(jié)點密度法得到的結(jié)構(gòu)中的節(jié)點密度進(jìn)行分類，討論最終結(jié)構(gòu)中節(jié)點密度的0，1分布問題.

從表1可以看出，采用影響區(qū)域節(jié)點密度插值法獲得的結(jié)果中，隨著插值半徑的擴大，“1”密度節(jié)點個數(shù)逐漸增多，即在高密度區(qū)域材料分布更為集中.而在低密度區(qū)域，“0”密度節(jié)點變少，中間密度節(jié)點變多，即在這些區(qū)域材料分布更為分散化.為此，可以在較大影響半徑下采用較高等密度線的方法獲得最終的優(yōu)化結(jié)構(gòu).如對圖4中影響區(qū)域半徑為1.5cm的優(yōu)化結(jié)果，采用等密度線獲得的最終結(jié)構(gòu)如圖5所示.

從圖5中可以看出，在采用影響區(qū)域插值方法時，較大插值半徑下，材料更加集中在一些核心區(qū)域，在此意義上，等密度線可以在結(jié)構(gòu)設(shè)計后期對最終結(jié)構(gòu)提供參考.

表1　結(jié)果中密度值分布

3.2算例2

為考察本文方法在三角形單元中的應(yīng)用，考慮如圖6圓盤結(jié)構(gòu)剛度最大化優(yōu)化問題，圓盤中間固定，半徑為50cm，載荷在圓周角30°、90°、150°、210°、270°、330°位置，方向為切向方向，大小為F=100N，材料彈性模量為E=210GPa，泊松比為v=0.3.

以目標(biāo)體積不超過設(shè)計區(qū)域體積的60%作為優(yōu)化目標(biāo)，采用三角形單元進(jìn)行離散，采用上文的方法獲得的優(yōu)化結(jié)果如圖7所示.

從圖7可以看出，本文提出的方法同樣適用于三角形網(wǎng)格單元，可以有效地處理非規(guī)則設(shè)計區(qū)域，并通過等密度值線獲得清晰的拓?fù)浣Y(jié)構(gòu).

4結(jié)語

針對變密度法拓?fù)鋬?yōu)化中出現(xiàn)的網(wǎng)格依賴性的現(xiàn)象，本文討論了基于影響區(qū)域節(jié)點密度插值的優(yōu)化方法.將設(shè)計變量定義在有限元網(wǎng)格節(jié)點位置，采用影響區(qū)域插值的方法獲得有限元單元的偽密度值，通過插值半徑的選取控制結(jié)構(gòu)中的最小尺寸，從而處理網(wǎng)格依賴性問題.并提出了基于密度等值線的清晰化方法，得到了0，1分布的優(yōu)化結(jié)果.經(jīng)過具體的數(shù)值算例，證明了本文方法可以在四邊形單元和三角形單元中有效地處理網(wǎng)格依賴性問題，獲得清晰的拓?fù)浣Y(jié)構(gòu).

參考文獻(xiàn)：

[1]GEORGE R， TOMASZ L. Topology optimization in structural and continuum mechanics [M]. Germany： Springer， 2014.

[2]SIGMUND O. Design of material structures using topology optimization [D]. Denmark： Department of Solid Mechanics， Technical University of Denmark，1994.

[3]HABER R B， BENDSE M P， JOG C. A new approach to variable-topology shape design using a constraint on the perimeter [J].StructuralOptimization， 1996，11(1)： 1-12.

[4]PETERSSON J， SIGMUND O. Slope constrained topology optimization [J].InternationalJournalforNumericalMethodsinEngineering， 1998，41(8)： 1417-1434.

[5]SUNG K Y， SANG H P. A study on the shape extraction process in the structural topology optimization using homogenized material [J].ComputersandStructures， 1997，62(3)： 527-538.

[6]RAHMATALLA S F， SWAN C C. A Q4/Q4 continuum structural topology optimization implementation [J].StructMultidiscOptim， 2004，27(1)： 130-135.

[7]GUEST J K， PREVOST J H， BELYTSCHKO T. Achieving minimum length scale in topology optimization using nodal design variables and projection functions [J].InternationalJournalforNumericalMethodsinEngineering， 2004，61(2)： 238-254.

[8]KANG Z， WANG Y Q. Structural topology optimization based on non-local shepard interpolation of density field [J].ComputerMethodsinAppliedMechanicsandEngineering， 2011，200(49)： 3515-3525.

[9]占金青，楊康，黃志超.基于節(jié)點變量法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計 [J].農(nóng)業(yè)機械學(xué)報，2014，45(9)： 329-333.

[10]YI J J， ZENG T， RONG J H，etal. A topology optimization method based on element independent nodal density [J].JCentSouthUniv， 2014，21(2)： 558-566.

[11]龍凱，陳廣華.基于物質(zhì)點描述的雙向漸進(jìn)式拓?fù)鋬?yōu)化方法 [J].工程力學(xué)，2012，29(8)： 308-312.

[12]SIGMUND O. A 99 line topology optimization code written in matlab [J].StructMultidiscOptim， 2001，21(2)： 120-127.

[13]BENDSOE M P. Optimization of structural topology， shape and material [M]. Germany： Springer， 1995.

文章編號：0427-7104(2016)03-0310-05

收稿日期：2016-01-05

基金項目：國家自然科學(xué)基金(11572089)

作者簡介：吉先磊(1990—)，男，碩士研究生；王皓(1963—)，男，副教授，通訊聯(lián)系人，E-mail： wanghao@fudan.edu.cn.

中圖分類號：O 342

文獻(xiàn)標(biāo)志碼：A

Topology Optimization Based on Regional Node Density Interpolation Method

JI Xianlei， ZHANG Meiyan， WANG Hao

(DepartmentofAeronauticsandAstronautics，F(xiàn)udanUniversity，Shanghai200433，China)

Abstract：Variable density method is a common approach in topology optimization， while there will be problems like checkerboard， grid dependence， gray cells and other instability occurs in the optimization procedure. In order to solve the checkerboard problem in topology optimization using variable density method， the node density is introduced in the optimization process as a design variable. A new method based on regional node density interpolation is proposed， as well as a node density based on clarity approach is introduced. The results show that the presented node density method can effectively eliminate the grid dependence problem as well as the checkerboard phenomenon in topology optimization.

Keywords：topology optimization； nodal density； regional interpolation