拓撲優(yōu)化均勻化方法的改進迭代算法

孫 蓓,蘇 超

(河海大學水利水電工程學院,江蘇南京 210098)

拓撲優(yōu)化方法是目前結(jié)構(gòu)優(yōu)化領域研究的熱點[1-5].均勻化方法[6]是連續(xù)體拓撲優(yōu)化方法中最常見的一種方法,該方法的一般迭代算法對所有連續(xù)體結(jié)構(gòu)都是可用的.若材料特性和荷載情況較為簡單,結(jié)構(gòu)規(guī)模較小,該方法迭代步較少,計算時間也較短,是一個極為有效的方法,但當結(jié)構(gòu)規(guī)模較大,尤其是一些復雜的空間結(jié)構(gòu)(如拱壩),其材料特性和荷載都會出現(xiàn)較為復雜的情況,計算量變大,單步的迭代時間會明顯增加,迭代步數(shù)也會增多,總計算負荷會很大.為了減少計算時間,提高計算效率和精度,本文對結(jié)構(gòu)拓撲優(yōu)化均勻化方法的一般迭代算法進行了改進,并引入了用于控制計算過程的2個密度閾值.

1 改進迭代算法的原理

結(jié)構(gòu)拓撲優(yōu)化均勻化方法[7]在采用有限元方法對設計區(qū)域進行離散化[8-9]的基礎上,將整個設計空間假設成類似“氣孔分布”的微結(jié)構(gòu)單元(單胞),單胞在優(yōu)化開始前分布均勻,大小相同.在拓撲優(yōu)化過程中,單胞密度分布發(fā)生變化,即高應力區(qū)域單胞密度變大,低應力區(qū)域單胞密度變小.優(yōu)化過程中形成了一種承重結(jié)構(gòu),這種結(jié)構(gòu)在高應力區(qū)域氣孔“密集”,在低應力區(qū)域氣孔密度較低[10].當?shù)嬎闳客瓿珊?再定義一個合理的密度最小值,然后剔除設計空間中單胞密度低于這個最小值的區(qū)域,將產(chǎn)生一個材料效應最高的重量優(yōu)化承重結(jié)構(gòu)[11].

以彈性板殼結(jié)構(gòu)為例,根據(jù)有限單元基本理論,可得彈性板殼單元的剛度矩陣

其中

式中:ki——單元剛度矩陣;Bi——單元轉(zhuǎn)換矩陣;Di——單元常量矩陣;ξ,η,ζ——形函數(shù)變量;Di(x,y,z)——整體坐標系下單元常量矩陣;Di(ξ,η,ζ)——單元局部坐標系下單元常量矩陣.

采用一般迭代算法時材料密度的求解公式為

其中

式中:ηi——第i個微結(jié)構(gòu)單元密度;W——初始結(jié)構(gòu)材料質(zhì)量;W*——刪除單元的總質(zhì)量;m——空洞單元數(shù)量;ui——單元節(jié)點位移矩陣;vi——單元體積;S——迭代步數(shù);vj——空洞單元體積;ε——極小量;C——常量.

對任何迭代計算過程來說,當某一計算數(shù)值足夠小時,計算機將視其為0,因此引入一個極小量ε,將空洞用密度為ε的單元來替代.材料體積約束(或質(zhì)量約束)計算時不包括這部分單元的體積(或質(zhì)量).采用一般迭代算法對結(jié)構(gòu)進行拓撲優(yōu)化時,那些迭代伊始密度就較小的單元,整個迭代過程單元密度始終呈減小趨勢,將在最終采用的結(jié)構(gòu)中刪除.這些單元的密度由于沒有小到計算機視為0的程度,因此每一步迭代計算都會消耗大量的計算資源.而那些迭代伊始密度就較大的單元,整個迭代過程單元密度亦始終保持較大或不斷增大趨勢,將在最終采用的結(jié)構(gòu)中保留.

假設一個較小的密度閾值,若優(yōu)化產(chǎn)生小于這一閾值的密度單元時,即將這些單元密度以極小值ε替代,視其為空洞,這樣就可以提高迭代計算的收斂速度.實際上,這些單元并不是空洞單元,是有一定密度的,若將其人為設定為空洞單元,則材料質(zhì)量約束計算時這部分質(zhì)量將不存在,這將使得替代后的總質(zhì)量大大小于替代前的總質(zhì)量.雖然從約束上講,這是可以接受的,但與實際情況不符,且計算時極有可能出現(xiàn)“早熟”等畸形收斂,計算結(jié)果將難以接受.因此,再給出一個較大的密度閾值,將小于這一閾值的單元密度以這一閾值替代.這種替代有2個優(yōu)點:一是解決了刪除較小密度單元帶來的替代后總質(zhì)量過小的問題;二是替代后單元密度值更為集中,剩余單元間密度差異減小,避免了收斂速度過快造成的畸形收斂,這就可在加速迭代的基礎上最大程度地保證結(jié)構(gòu)拓撲優(yōu)化迭代計算的穩(wěn)定性.為了區(qū)別2個密度閾值,本文將較小的密度閾值稱為MFS,較大的密度閾值稱為MFL.

以圖1說明改進迭代算法的原理.圖1(a)顯示了微結(jié)構(gòu)單元的密度曲線L,橫坐標n為微結(jié)構(gòu)單元數(shù)量,縱坐標 ηi為第i個微結(jié)構(gòu)單元密度,單元密度在0～1之間.圖1(b)顯示了替代前微結(jié)構(gòu)單元的密度曲線L,陰影部分顯示的面積為總的單元質(zhì)量.圖1(c)顯示了用MFS和MFL2個密度閾值進行替代的微結(jié)構(gòu)單元的密度曲線L,陰影部分顯示的面積仍為總的單元質(zhì)量,但此總的單元質(zhì)量要小于或等于圖1(b)所示的總的單元質(zhì)量.

圖1 改進迭代算法的原理Fig.1 Principles of modified iterative method

2 改進迭代算法的流程和算式

改進迭代算法的流程如圖2所示.圖中:Ws0為初始做結(jié)構(gòu)材料質(zhì)量;Wsl為使用密度閾值替代后微結(jié)構(gòu)材料質(zhì)量.

對式(2)進行修正,得到改進迭代算法的算式,即

圖2 改進迭代算法的流程Fig.2 Flow chart of modified iterative method

3 算例及分析

3.1 算例1

圖3(a)所示為2m×2m的矩形平面結(jié)構(gòu),材料彈性模量為20GPa,左側(cè)邊界固定,右側(cè)下部作用一鉛直向下的集中荷載P=1kN.圖3(b)所示為有限元網(wǎng)格,結(jié)構(gòu)共劃分為40×40個矩形單元.

刪除率選擇為80%,極小值ε取0.001.首先采用結(jié)構(gòu)拓撲優(yōu)化均勻化方法的一般迭代算法進行計算,計算結(jié)果如圖4所示,保留的微結(jié)構(gòu)單元密度為0.1～1;然后采用改進迭代算法進行計算,計算時 MFS取0.1,MFL取0.4,計算結(jié)果如圖5所示.

圖3 算例1Fig.3 Plan and meshes of Example 1

圖4 算例1一般迭代算法的計算結(jié)果Fig.4 Results of Example1 of general iterative method

圖5 算例1改進迭代算法的計算結(jié)果Fig.5 Results of Example 1 of modified iterativemethod

采用一般迭代算法,共需迭代48步.為了能夠清楚地看出計算過程,圖中顯示了密度大于0.1的單元.而實際上,在一般迭代算法的計算中,除了密度小于極小值的單元外,其余所有單元都是參與計算的.對于密度為0.1～1的微結(jié)構(gòu)單元,隨著迭代的不斷進行,其有限單元的數(shù)量會不斷減少,達到優(yōu)化目標時迭代結(jié)束.

采用改進迭代算法,共需迭代19步.由于MFS取為0.1,因此選擇顯示的單元密度為0.1～1.

比較圖4(a)和圖5(a)可知,當S=6步時,改進迭代算法的結(jié)果比一般迭代算法的結(jié)果外形更為光滑.圖5(b)顯示,當S=11步時,改進迭代算法的結(jié)果已非常接近于最終優(yōu)化結(jié)果.而對于一般迭代算法來說,迭代30步的結(jié)果也與最終的優(yōu)化結(jié)果有較大差異.圖4(c)與圖5(c)相比,外形上兩者結(jié)果基本一致,僅結(jié)構(gòu)個別部位有微小差異.隨著結(jié)構(gòu)復雜性和荷載的增加以及MFS和MFL取值的不同,這些微小差異會有所增大.由于采用了不同的算法,這些微小差異是不可避免的,因此,當出現(xiàn)這些微小差異時,可認為優(yōu)化計算結(jié)果一致.

從本例來看,改進迭代算法和一般迭代算法的優(yōu)化計算結(jié)果一致,而其迭代步數(shù)僅為一般迭代算法的40%左右,計算時間較一般迭代算法減少50%以上,且計算時占用的計算資源大大減少,這使得復雜的大型結(jié)構(gòu)在拓撲優(yōu)化計算時可以采用更細小的網(wǎng)格.

3.2 算例2

采用不同的MFS和MFL進行計算,并根據(jù)計算結(jié)果給出MFS和MFL的取值范圍,以避免因密度閾值的不確定而過多地耗費時間.

圖6所示為1.6m×1m的矩形平面結(jié)構(gòu),材料彈性模量為20GPa,3個集中荷載(P1=P2=P3=1kN)作用于結(jié)構(gòu)的下邊界.結(jié)構(gòu)共劃分為64×40個矩形板單元,下方左右兩端受固定約束.刪除率取為80%,極小值ε取為0.001.圖7給出了一般迭代算法的計算結(jié)果.不同MFS和M FL取值改進迭代算法的計算結(jié)果如圖8所示.

從圖8(a)～(d)可知,在 MFL相同的情況下,隨著 MFS的增大,迭代步明顯減少,且當MFS和MFL取值都較大時,迭代收斂情況明顯優(yōu)于MFS和MFL取值較小時的收斂情況.

從圖8(b)～(e)可知,在 MFS相同的情況下,隨著 MFL的增大,迭代步也相應減少,減少的幅度較M FL相同而M FS增大情況下的幅度小很多.圖8(b)與(c)顯示的結(jié)果可認為是一致的.圖8(d)與(e)的 MFS和MFL取值均比圖8(b)和(c)大,但從優(yōu)化結(jié)果來看,并非僅有微小差別,由于算例結(jié)構(gòu)和荷載的簡單性,這里不能再認為兩者的優(yōu)化結(jié)果一致.圖8(e)的結(jié)果從單元數(shù)量來說,明顯比圖8(d)的要少.這一差別主要體現(xiàn)在結(jié)構(gòu)中部受荷載區(qū)域.該區(qū)域內(nèi)的V形結(jié)構(gòu)較為單薄,這一結(jié)果是不可用的.由于其脆弱部分很少,可以人為地進行彌補.

圖8(f)給出了一個MFS和MFL都較大的極端情況,在此情況下,拓撲優(yōu)化計算迭代僅需9步,但結(jié)果卻很不理想.其主要問題是優(yōu)化得到的結(jié)構(gòu)V形區(qū)域部分十分脆弱,有較多單元都是通過角連接的,不能為設計人員所接受.

圖6 算例2平面結(jié)構(gòu)Fig.6 Plan of Example2

圖7 算例2一般迭代算法的計算結(jié)果Fig.7 Results of Example 2 of general iterative Method

圖8 算例2改進迭代算法的計算結(jié)果Fig.8 Results of Example 2 of modified iterativemethod

4 結(jié) 論

a.當M FS和M FL取值合理時,拓撲優(yōu)化均勻化方法的改進迭代算法所需要的迭代步數(shù)遠小于一般迭代算法,收斂速度也大為提高.

b.從理論上講,一般迭代算法的計算過程是沒有人為干預的,所得優(yōu)化結(jié)果應該是最合理的;而改進迭代算法的計算過程受到人為的干預,干預程度與MFS和MFL的取值大小有關(guān),這種干預勢必會影響優(yōu)化結(jié)果的合理性.但是從算例的計算結(jié)果來看,當MFS取0.1,MFL取0.3時,改進迭代算法得到的優(yōu)化結(jié)果與一般迭代算法得到的優(yōu)化結(jié)果相似程度最高.因此,當M FS和M FL取值合理時,改進迭代算法計算過程干預對優(yōu)化結(jié)果的影響可以忽略.

c.當 MFS和 MFL取值合理時,增大 MFS和 MFL均會提高收斂速度,而增大 MFS對提高收斂速度的作用要比增大MFL對提高收斂速度的作用大.當MFL相同時,MFS的微小變化也會對迭代步數(shù)產(chǎn)生較大影響;而當M FS相同時,增大M FL盡管也可以加速收斂,但是對迭代步數(shù)不會產(chǎn)生決定性影響.

d.M FS取值范圍定為0.1～0.2,M FL的取值范圍定為0.3～0.45比較合理.對于大型計算來說,迭代過程中MFS和MFL的取值可以根據(jù)計算實際情況在一定范圍內(nèi)進行調(diào)整,這可以加速收斂.

e.本文對拓撲優(yōu)化均勻化方法的一般迭代算法進行了改進,為了加強對迭代過程的控制,引入了2個密度閾值,并給出了密度閾值的建議取值范圍.算例表明,這一改進迭代算法可行,當采用合理的密度閾值時,拓撲優(yōu)化計算收斂速度會明顯提高.

[1]TALICHI C,PAULINO G H,LE CH.Honeycomb wachspress finite elements for structural topology optimization[J].Structural and Multidisciplinary Optimization,2009,37(6):569-583.

[2]GUO X,YAMAZAKIK,CHENGG.A new approach for the solution of singular optimain truss topology optimizationwith stress and local buckling constraints[J].Structural and Multidisciplinary Optimization,2001,22(5):364-373.

[3]MAIERG,BUACA C,PATTAVANA A.Topology-information periodic updates in multi-domain ASON networks with topology aggregation[J].Fiber and Integrated Optics,2008,27(4):265-277.

[4]FITZWATERL,KHALIL R,HUNTER E,et al.Topology optimization risk reduction[C]//Annual Forum Proceedings-AHSInternational,AHS International 64th Annual Forum.Montreal,Canada:American Helicopter Society,2008:543-556.

[5]張忠中,蘇超.結(jié)構(gòu)拓撲優(yōu)化方法在拱式渡槽拱軸線設計中的應用[J].河海大學學報:自然科學版,2007,35(4):444-447.(ZHANGZhong-zhong,SUChao.Continuum topology optimizationmethod for arch axisdesign of arciform aqueducts[J].Journal of Hohai University:Natural Sciences,2007,35(4):444-447.(in Chinese)).

[6]BENDEOEM P,KIKUCHI N.Generatingoptimal topologies in structural design usinga homogenization method[J].Computer Methods in Applied Mechanics and Engineering,1988,71(2):197-224.

[7]BENDSOE M P,MOTA S.Topology design of structures[M].NATO ASI:Kluwer Academic Publishers,1993.

[8]HASSANI B,HINTON E.A review of homogenization and topology optimization I:homogenization theory for media with periodic[J].Computers&Structures,1998,69:707-717.

[9]BENDSOEM P,SIGMUNDO.Material interpolations in topology optimization[J].Archiveof Applied Mechanics,1999,69:635-654.

[10]YOUN SK,PARK S H.A study on the shape extraction process in the structural topology optimization using homogenized material[J].Computers&Structures,1997,62(3):527-538.

[11]BREMICKERM,KIKUCHI N,CHIREHDAST M,et al.Integrated topology and shape optimization in structural design[J].Journal of Structural Mechanics,1991,19(4):551-587.