基于類桁架材料模型優(yōu)化方法計算穩(wěn)定性研究

李霞,周克民

(華僑大學土木工程學院,福建 廈門 361021)

0 引言

在結(jié)構(gòu)優(yōu)化領(lǐng)域,Michell的理論具有里程碑的意義[1],不僅揭示了優(yōu)化結(jié)構(gòu)的類桁架性質(zhì),而且一直作為檢驗其它方法的可靠依據(jù).但是Michell理論采用解析方法求解困難,得到的拓撲優(yōu)化結(jié)構(gòu)是非均勻各向異性連續(xù)體,并非工程意義上的桁架,難以在工程實際中應(yīng)用.隨著計算機技術(shù)、有限元方法和數(shù)學優(yōu)化理論的發(fā)展,結(jié)構(gòu)拓撲優(yōu)化研究取得了許多重要進展,主要集中在基于有限元的數(shù)值方法.人們試圖通過在等厚板中尋找孔的優(yōu)化分布實現(xiàn)拓撲優(yōu)化,包括均勻化方法[2],進化結(jié)構(gòu)優(yōu)化方法[3],水平集方法[4]和ICM方法[5]等.為了得到等厚板,普遍需要抑制中間密度,抑制中間密度直接導致了單元鉸接,棋盤格現(xiàn)象、單元依賴性等數(shù)值不穩(wěn)定問題[6-7].雖然自由材料優(yōu)化方法不存在數(shù)值不穩(wěn)定問題,但優(yōu)化結(jié)果缺乏與實際結(jié)構(gòu)的明確對應(yīng)關(guān)系[8-9].

按照Michell理論,拓撲優(yōu)化結(jié)構(gòu)一般是由無限多致密桿件構(gòu)成的非均勻各向異性連續(xù)體,各種拓撲優(yōu)化數(shù)值方法試圖用均勻等厚帶孔板逼近這種非均勻各向異性連續(xù)體.在有限元理論的基本假設(shè)下,這種黑白相間的“棋盤格”狀態(tài)是很好的近似.隨著單元網(wǎng)格的增加,拓撲描述清晰度的提高,出現(xiàn)更多的細部結(jié)構(gòu),出現(xiàn)優(yōu)化結(jié)果的網(wǎng)格依賴性.因此,出現(xiàn)這些數(shù)值不穩(wěn)定問題是必然的.目前也有一些方法,如周長控制、梯度控制等可以避免或減輕這些現(xiàn)象.理論上的拓撲優(yōu)化結(jié)構(gòu)是各向異性連續(xù)體,實際需要均勻等厚帶孔板,這兩者之間的矛盾是出現(xiàn)這些問題的根本原因.

本作者曾提出基于類桁架材料模型的拓撲優(yōu)化方法[10],試圖通過優(yōu)化類桁架材料的分布場實現(xiàn)結(jié)構(gòu)的拓撲優(yōu)化,而不是直接優(yōu)化等厚板中孔的分布.采用有限元的數(shù)值分析方法構(gòu)造設(shè)計域內(nèi)桿件非均勻連續(xù)分布場,從而直接找到理論上的最優(yōu)結(jié)構(gòu)-Michell桁架.由于不抑制中間密度,所以完全沒有數(shù)值不穩(wěn)定問題.為了滿足工程實際需要,再將類桁架連續(xù)體轉(zhuǎn)化為離散結(jié)構(gòu).一些算例表明這種離散的誤差不大[11-12].本文中目的在于分析該方法的穩(wěn)定性,即對有限元網(wǎng)格的無關(guān)性以及棋盤格現(xiàn)象等,這是傳統(tǒng)方法普遍存在的問題[13].

1 類桁架材料模型

為了構(gòu)造類桁架連續(xù)體結(jié)構(gòu),首先建立類桁架連續(xù)體材料模型.假設(shè)在弱的基材料內(nèi)非均勻連續(xù)分布致密的桿件,桿件的分布密度和方向在設(shè)計域內(nèi)連續(xù)變化.基材料不受力,力完全由分布桿件承擔.優(yōu)化的目的是建立這種桿件分布場.對于本文中研究的單工況應(yīng)力約束體積最小問題,假設(shè)桿件在任意點沿兩個正交方向分布桿件.該問題是最經(jīng)典的結(jié)構(gòu)拓撲優(yōu)化問題,經(jīng)常作為驗證各種拓撲優(yōu)化方法的標準算例,其優(yōu)化結(jié)果常稱為Michell桁架.

1.1 彈性矩陣

σi=Etiεi,i=1,2

(1)

在單工況下,考慮到類桁架材料中相鄰平行桿件之間沒有相互作用,所以橫向變形系數(shù)和剪切剛度均為零.為了避免結(jié)構(gòu)剛度矩陣奇異,假設(shè)切應(yīng)力τ12和切應(yīng)變γ12之間的關(guān)系為,

τ12=E(t1+t2)γ12/4

(2)

如此假設(shè)條件下,該材料模型當t1=t2時可以表示各向同性材料.實際上在優(yōu)化結(jié)構(gòu)中,由于主應(yīng)力沿材料主軸方向,所以材料主方向無切應(yīng)變,剪切剛度不起作用.因此,材料主軸方向的應(yīng)力應(yīng)變的可以表示為

[σ1σ2τ12]=D(t1,t2,0)[ε1ε2γ12]

(3)

式中D(t1,t2,0)是材料主軸方向的彈性矩陣,D(t1,t2,0)=E·diag[t1t2(t1+t2)/4]

(4)

diag表示對角矩陣.

在多工況下,桿件一般不正交,也存在切應(yīng)變.作為一種近似,本文中仍采用了以上假設(shè).

除對申請書進行評審之外，二審專家組的另一任務(wù)是在評審結(jié)束后將各領(lǐng)域的相對資助率進行綜合，平衡各研究領(lǐng)域的資助力度。

(5)

可以建立結(jié)構(gòu)坐標系下的彈性矩陣D(t1,t2,α)=TT(α)D(t1,t2,0)T(α)

(6)

為了簡單具體寫出式(6),定義幾個常數(shù)矩陣,

(7)

和函數(shù)矩陣

g(α)=[cos2αsin2α1]

(8)

(9)

式中sbr和gr分別是式(7)和式(8)定義矩陣的分量.

1.2剛度矩陣如果桿件在結(jié)點j位置的密度和方向分別為t1j,t2j和αj,彈性矩陣記作D(t1j,t2j,αj),單元內(nèi)任意點(ξ,η)的彈性矩陣De(ξ,η)可以通過形函數(shù)Nj(ξ,η)插值計算

(10)

式中Se為屬于單元e的結(jié)點集合.定義與材料分布無關(guān)的常數(shù)矩陣

(11)

(12)

2 優(yōu)化方法

以桿件在結(jié)點位置的密度和方向t1j,t2j,αj(j=1,2,…,J)為優(yōu)化設(shè)計變量,單元內(nèi)部任意點的桿件密度由結(jié)點位置的桿件密度插值得到.將桿件密度在設(shè)計域內(nèi)積分得到結(jié)構(gòu)體積

(13)

式中Sj是圍繞結(jié)點j的單元集合,J是結(jié)點總數(shù),zj是形函數(shù)在圍繞結(jié)點j的單元內(nèi)桿件體積積分之和

(14)

式中Ve是單元面積.本文中研究應(yīng)力約束體積最小問題

(15)

式中l(wèi)表示工況,L為工況總數(shù),σp為允許應(yīng)力.

2.1單工況優(yōu)化采用滿應(yīng)力準則法優(yōu)化桿件在結(jié)點位置的密度

(16)

2.2多工況優(yōu)化多工況下仍假設(shè)優(yōu)化結(jié)構(gòu)中任意點有兩個方向的正交桿件.首先在各單工況下分別得到桿件優(yōu)化分布t1l,t2l,αl.按照坐標變換的概念可以寫出任意方向的剛度矩陣

D(φ;t1l,t2l,αl)=D(t1l,t2l,φ-αl),

(17)

取對角元第1行第1列的元素為度量基礎(chǔ),定義為方向剛度,各工況下方向剛度的包絡(luò)線方程為

(18)

多工況下任意點位置的桿件密度和方向記作x1,x2,θ,定義多工況下方向剛度與各單工況下方向剛度的包絡(luò)線的差值為

(19)

假設(shè)多工況下的優(yōu)化結(jié)構(gòu)與各單工況下的方向最大剛度差值最小,即取住值

(20)

可以解得多工況下桿件密度和方向角

x1+x2=2I0x1-x2=4(I1sin2θ+I2cos2θ), tan2θ=I1/I2

(21)

(22)

3 數(shù)值算例

下面通過簡支結(jié)構(gòu)、懸臂結(jié)構(gòu)、多工況3個經(jīng)典算例分析不同網(wǎng)絡(luò)劃分，即不同單元密度對優(yōu)化結(jié)果的影響.由于問題與量綱無關(guān)，這里的所有數(shù)據(jù)均無量綱，力和彈性模量均取1.圖1是矩形域下邊兩端簡支，中點作用豎向集中力的簡支結(jié)構(gòu)，圖2(a)、圖2(b)是簡支結(jié)構(gòu)采用2種不同的4結(jié)點矩形規(guī)則單元密度的優(yōu)化結(jié)果，圖2(c)是解析解[14].圖3是矩形域左邊固定，右邊中點作用豎向集中力的懸臂結(jié)構(gòu)，圖4(a)、圖4(b)、圖4(c)是懸臂結(jié)構(gòu)采用3種不同的規(guī)則單元密度的優(yōu)化結(jié)果，圖4(d)是解析解[15].圖5是懸臂結(jié)構(gòu)采用2種不同的不規(guī)則單元的優(yōu)化結(jié)果.圖6是2個力分別作用在右邊上下2個角點的多工況，圖7是多工況采用2種不同的規(guī)則單元密度的優(yōu)化結(jié)果.表1給出了優(yōu)化結(jié)構(gòu)所使用的材料體積及迭代次數(shù).

表1 單元數(shù)量及優(yōu)化結(jié)果

圖1 簡支結(jié)構(gòu)

圖2 簡支結(jié)構(gòu)不同密度規(guī)則單元的優(yōu)化結(jié)果(a)20×12;(b)40×24;(c)解析解.

圖3 懸臂結(jié)構(gòu)

圖4 懸臂結(jié)構(gòu)不同密度規(guī)則單元的優(yōu)化結(jié)果(a)10×16;(b)20×32;(c)40×64;(d)解析解.

圖5 懸臂結(jié)構(gòu)不同的不規(guī)則單元的優(yōu)化結(jié)果(a)10×16;(b)10×16.

圖6 多工況結(jié)構(gòu)

圖7 多工況不同密度規(guī)則單元的優(yōu)化結(jié)果(a)20×32;(b)40×64.

比較這些優(yōu)化結(jié)果可以看出,其數(shù)據(jù)均非常接近,單元網(wǎng)格的密度和形狀對優(yōu)化結(jié)果影響很小.所有算例均未出現(xiàn)數(shù)值不穩(wěn)定問題,驗證了基于類桁架連續(xù)體材料模型優(yōu)化方法數(shù)值計算結(jié)果的穩(wěn)定性.

4 結(jié)語

數(shù)值不穩(wěn)定問題是各種數(shù)值拓撲優(yōu)化方法普遍存在的問題,雖然可以采取一些方法減輕這種問題,但不可能完全解決,且增加了計算量.基于類桁架連續(xù)體材料模型的拓撲優(yōu)化方法沒有單元鉸接、棋盤格現(xiàn)象,其結(jié)果受網(wǎng)格劃分影響很小.

[1] Michell A. The limits of economy of material in frame structure[J]. Phil Mag,1904,8(6):589-597.

[2] Bendsoe M, Kikuchi N. Generating optimal topologies in structural design using a homogenization method[J]. Comput Methods Appl Mech Engrg,1988,71(2):19-224.

[3] Querin O M, Young V, Steven G P, et al. Computational efficiency and validation of bi-directional evolutionary structural optimization[J]. Comput Methods Appl Mech Engrg,2000,189(2):559-573.

[4] Sethian J, Wiegmann A. Structural boundary design via level set and immersed interface methods[J]. J Comput Phys,1999,163:489-528.

[5] Sui Y, Yang D. A new method for structural topological optimization based on the concept of independent continuous variables and smooth model[J]. Acta Mechanica Sinica,1998,14(2):179-185.

[6] Sigmund O, Petersson J. Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima[J]. Struct Optim,1998,16(1):68-75.

[7] Cheng K, Olhoff N. An investigation concerning optimal design of solid elastic plates[J]. Int J Solids Structures,1981,17(3):305-323.

[8] Matsui K, Terada K. Continuous approximation of material distribution for topology optimization[J]. Int J Numer Meth Engng,2004,59(4):1925-1944.

[9] Bendsoe M, Sigmund O. Material interpolations in topology optimization[J]. Arch Appl Mech,1999,69(9-10):635-654.

[10] Zhou K, Li X. Topology optimization of structures under multiple load cases using fiber-reinforced composite material model[J]. Comput Mech,2005,38(2):163-170.

[11] Zhou K, Li J. The exact weight of discretized michell trusses for a central point load[J]. Struct Optim,2004,28(1):69-72.

[12] Prager W. A note on discretized michell structures[J]. Comput Methods Appl Mech Eng,1974(3):349-355.

[13] Eschenauer H, Olhoff N. Topology optimization of continuum structures: a review[J]. Applied Mechanics Reviews,2001,54(4):331-389.

[14] Hemp W. Optimum structure[M]. Oxford, Clarendon Press,1973.

[15] Lewinski T, Zhou M, Rozvany G. Extended exact solutions for least-weight truss layouts — paper Ⅰ: cantilever with a horizontal axis of symmetry[J]. Int J Mech Sci,1994,36(5):375-398.