基于影響區(qū)域節(jié)點密度插值的結構拓撲優(yōu)化方法

吉先磊，張美艷，王　皓

(復旦大學 航空航天系，上海 200433)

?

基于影響區(qū)域節(jié)點密度插值的結構拓撲優(yōu)化方法

吉先磊，張美艷，王皓

(復旦大學 航空航天系，上海 200433)

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

關鍵詞：拓撲優(yōu)化； 節(jié)點密度； 區(qū)域插值

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

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

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

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

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

(1)

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

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

ρe=ωe，iρi，

(2)

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

(3)

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

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

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

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

(4)

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

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

(5)

通過鏈式求導法則，進而可以表示出目標函數對節(jié)點密度變量的靈敏度為

(6)

(7)

3數值算例

下面采用2維平面拓撲優(yōu)化算例對本文方法進行討論.

3.1算例1

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

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

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

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

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

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

表1　結果中密度值分布

3.2算例2

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

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

從圖7可以看出，本文提出的方法同樣適用于三角形網格單元，可以有效地處理非規(guī)則設計區(qū)域，并通過等密度值線獲得清晰的拓撲結構.

4結語

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

參考文獻：

[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ù)體結構拓撲優(yōu)化設計 [J].農業(yè)機械學報，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]龍凱，陳廣華.基于物質點描述的雙向漸進式拓撲優(yōu)化方法 [J].工程力學，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

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

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

中圖分類號：O 342

文獻標志碼：A

Topology Optimization Based on Regional Node Density Interpolation Method

JI Xianlei， ZHANG Meiyan， WANG Hao

(DepartmentofAeronauticsandAstronautics，FudanUniversity，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