應(yīng)用類桁架模型的連續(xù)體拓?fù)鋬?yōu)化方法

鄭偉偉，周克民

（華僑大學(xué) 土木工程學(xué)院，福建 廈門361021）

結(jié)構(gòu)拓?fù)鋬?yōu)化的研究包括均勻化方法［1］、演化優(yōu)化算法［2］、水平集方法［3］、ICM 法［4］等.為了獲得清晰的結(jié)構(gòu)，抑制中間密度單元，這會(huì)引起了數(shù)值不穩(wěn)定現(xiàn)象，如棋盤格現(xiàn)象［5］，需要格外的技術(shù)解決這個(gè)問題［6-8］.Michell［9］揭示了拓?fù)鋬?yōu)化結(jié)構(gòu)是非均勻質(zhì)各向異性類桁架連續(xù)體，即由無限密的桿件構(gòu)成，這一特性啟發(fā)了從類桁架連續(xù)體中獲取最優(yōu)拓?fù)涞南敕ǎ?0-11］.工程中更多地使用帶孔等厚板，通過剔除密度小于一定值的節(jié)點(diǎn)形成孔洞.演化優(yōu)化算法［12-13］有采用過這種做法，但是由于演化優(yōu)化算法采用的是各向同性材料，且設(shè)計(jì)變量為單元的密度，因此，該算法得到的結(jié)果會(huì)出現(xiàn)棋盤格現(xiàn)象且與網(wǎng)格密度有關(guān).本文通過優(yōu)化孔洞的分布形成帶孔等厚板的優(yōu)化算法，在優(yōu)化迭代過程中逐步由類桁架連續(xù)體演化為帶孔板，使優(yōu)化結(jié)果更接近工程需要.

1 類桁架連續(xù)體材料模型

1.1 彈性矩陣

在類桁架連續(xù)體材料模型中，材料模型是由無限多的非均勻連續(xù)分布的桿件構(gòu)成，假設(shè)任意點(diǎn)都是由兩組正交的桿件構(gòu)成，它們?cè)谌我恻c(diǎn)的密度用t1和t2表示.假設(shè)應(yīng)力σi和應(yīng)變?chǔ)舏的線性關(guān)系為

式（1）中：E是彈性模量.

式（2）中：diag［·］表示對(duì)角矩陣.

在這種假設(shè)下，當(dāng)t1＝t2時(shí)，式（2）的各項(xiàng)為同性材料.若假設(shè)結(jié)構(gòu)坐標(biāo)軸到材料坐標(biāo)軸的轉(zhuǎn)角為α，那么結(jié)構(gòu)坐標(biāo)系中的彈性矩陣可由坐標(biāo)轉(zhuǎn)換矩陣獲得，即

式（3）中：T（α）為坐標(biāo)轉(zhuǎn)換矩陣，即

彈性矩陣（3）可以表示為

式（5）中：sb，r和gr（α）分別為常數(shù)矩陣和函數(shù)矩陣分量，Ar為常數(shù)矩陣，有

節(jié)點(diǎn)j（j＝1，2，…，J）處的桿件密度t1，j，t2，j和方向αj作為設(shè)計(jì)變量，J是節(jié)點(diǎn)的總數(shù)，設(shè)計(jì)變量節(jié)點(diǎn)j處的彈性矩陣可以表示為

有限元內(nèi)部任意一點(diǎn)的彈性矩陣可以由該單元節(jié)點(diǎn)位置的彈性矩陣插值得到

式（7）中：Nj（ξ，η）是型函數(shù)；ξ，η是局部坐標(biāo)；Se是屬于單元e的節(jié)點(diǎn)的集合.

將式（6）帶入式（7），得到

1.2 剛度矩陣

將單元?jiǎng)偠染仃嚕?）帶入剛度矩陣，有

2 優(yōu)化方法

將節(jié)點(diǎn)處的桿件密度和方向作為設(shè)計(jì)變量，材料的體積作為目標(biāo)函數(shù)，每一節(jié)點(diǎn)施加應(yīng)力約束，優(yōu)化問題的列式為

式（12）中：σb，j為節(jié)點(diǎn)處桿件b的應(yīng)力；σp為允許應(yīng)力.

迭代過程有如下5個(gè)步驟.

步驟1 設(shè)計(jì)域被劃分為有限個(gè)單元，單元采用4節(jié)點(diǎn)矩形單元，桿件初始密度和方向設(shè)置為

上標(biāo)0代表迭代次數(shù)，初始密度剔除標(biāo)準(zhǔn)設(shè)為0；

步驟2 通過有限元方法，計(jì)算節(jié)點(diǎn)處的主應(yīng)力方向和主應(yīng)力方向上的應(yīng)變的大??；

步驟3 利用滿應(yīng)力準(zhǔn)則，優(yōu)化桿件的密度，桿件的方向與主應(yīng)力的方向一致，即

步驟4 為了防止剛度矩陣奇異，設(shè)置低于密度剔除標(biāo)準(zhǔn)的桿件密度為t0（這里取10-6），而不是直接設(shè)置為0，高于平均密度的桿件的密度將設(shè)置為

步驟5 密度剔除標(biāo)準(zhǔn).隨迭代次數(shù)而增加，tic＋1＝tic＋Δt，Δt是密度剔除標(biāo)準(zhǔn)的增量.

通過大量的數(shù)值計(jì)算，取10-4比較合適.太小則孔洞形成的較慢，需要更多步的迭代；太大則形成的拓?fù)渑c最優(yōu)拓?fù)湎嗖钶^大.重復(fù)步驟2～5，直到?jīng)]有密度的節(jié)點(diǎn)數(shù)占總節(jié)點(diǎn)數(shù)的60%，剔除材料比例可以由用戶根據(jù)需要自行決定.

圖1 懸臂梁尺寸Fig.1 Dimension of cantilever beam

3 數(shù)值算例

運(yùn)用文中的算法，計(jì)算兩個(gè)經(jīng)典的算例，兩個(gè)算例中：彈性模量E＝210 GPa；允許應(yīng)力σp＝160 MPa.

第1個(gè)算例是一個(gè)右端中點(diǎn)受集中荷載的懸臂梁，尺寸如圖1所示.圖1中：厚度（h）為0.01 m；設(shè)計(jì)域被劃分為64×40矩形單元.所有節(jié)點(diǎn)密度（兩個(gè)方向的桿件密度之和）的等值線填充圖，如圖2，3所示.圖2，3中：等值線的間距為最大密度的5%.圖2是迭代31次后的結(jié)果，設(shè)計(jì)域開始出現(xiàn)一些孔.圖3是迭代62次后的結(jié)果，孔不斷擴(kuò)展、溶合，結(jié)構(gòu)拓?fù)溱呌诤?jiǎn)單.

圖2 迭代31次后的密度等值線圖Fig.2 Isogram of density after 31 iterations

圖3 迭代62次后的密度等值線圖Fig.3 Isogram of density after 62 iterations

大部分區(qū)域的密度分布較均勻，只有支座和荷載作用點(diǎn)附近區(qū)域密度較大，但沒有出現(xiàn)嚴(yán)重的應(yīng)力集中現(xiàn)象.體積比的定義為優(yōu)化后結(jié)構(gòu)總體積占初始體積的百分比.體積比（Rv）和保留的節(jié)點(diǎn)數(shù)百分比（Rn）隨著迭代次數(shù)（N）的變化曲線，如圖4所示.

從圖4可知：由于荷載比較小，優(yōu)化后纖維密度較初始密度小很多，優(yōu)化后的體積只占初始體積很小一部分.保留的節(jié)點(diǎn)百分比指密度大于臨界值（這里取2t0）的節(jié)點(diǎn)占總節(jié)點(diǎn)數(shù)的百分比，反映了類桁架連續(xù)體演化為帶孔板的進(jìn)程，當(dāng)其等于40%時(shí)，迭代停止.滿應(yīng)力優(yōu)化準(zhǔn)則可以在少量的迭代下使結(jié)構(gòu)體積達(dá)到最小，通過密度限值約束，去除小密度桿件，結(jié)構(gòu)離散得越來越簡(jiǎn)單，但體積只增加了一點(diǎn).

第2個(gè)算例是跨中受集中荷載的簡(jiǎn)支梁，尺寸如圖5所示.圖5中：厚度（h）為0.01 m；設(shè)計(jì)域被矩形單元?jiǎng)澐譃?0×40.

圖4 體積比及節(jié)點(diǎn)數(shù)比的變化Fig.4 Variation of volume ratio and nodes ratio

未經(jīng)密度限值約束迭代100次得到的優(yōu)化結(jié)果，如圖6所示.圖6中：為了使荷載集中點(diǎn)的密度不至于大其他區(qū)域的密度太多倍，仍施加了密度上限為平均密度的約束.從圖6可知：未經(jīng)密度限值約束的優(yōu)化結(jié)果形成連續(xù)材料場(chǎng)，而不能形成桿件.

增加密度限值約束后迭代56次的結(jié)果，如圖7所示.對(duì)比圖6，7可以看出：增加密度限值約束后，在設(shè)計(jì)域的中間部分形成了幾根桿件，而且與圓拱的密度基本一致，比較均勻.

圖5 簡(jiǎn)支梁尺寸Fig.5 Dimension of simply supported beam

圖6 沒有密度限值約束的優(yōu)化結(jié)果Fig.6 Optimal result without density restriction

圖7 經(jīng)密度限值約束的優(yōu)化結(jié)果Fig.7 Optimal result with density restriction

體積比和保留節(jié)點(diǎn)數(shù)百分比這兩項(xiàng)指標(biāo)隨迭代次數(shù)的變化，及施加密度限值約束前后的對(duì)比，如圖8所示.從圖8可知：施加密度限值約束后的體積比會(huì)比施加前稍微大一些，且體積比會(huì)隨著迭代次數(shù)逐步提高；對(duì)比施加密度限值約束前后節(jié)點(diǎn)數(shù)比，非零密度（大于2t0）節(jié)點(diǎn)比在施加密度限值約束后比施加前減少快很多.隨著迭代次數(shù)的變化線形，這一指標(biāo)決定于密度剔除標(biāo)準(zhǔn)的確定方法，因此，可以通過改變密度剔除標(biāo)準(zhǔn)的函數(shù)形式（線性或非線性）來控制類桁架連續(xù)體向帶孔板演化的過程.

圖8 施加密度限值約束前后指標(biāo)變化及對(duì)比Fig.8 Indexes variation and contrast before and after density restriction

4 結(jié)束語(yǔ)

研究了基于類桁架連續(xù)體材料模型的帶孔板拓?fù)鋬?yōu)化方法，通過在每一迭代步中剔除小于一定值的密度形成孔洞，限定大于一定值的密度達(dá)到密度均勻的目的，最終得到帶孔板.下一步工作是對(duì)連續(xù)體施加最小寬度的約束，從而可以得到孔數(shù)量不同的連續(xù)體拓?fù)鋬?yōu)化結(jié)果.

［1］ BENDS?E M P，KIKUCHI N.Generating optimal topologies in structural design using a homogenization method［J］.Comput Methods Appl Mech Eng，1988，71（2）：197-224.

［2］ XIE Yin-min，STEVEN G P.A simple evolutionary procedure for structural optimization［J］.Comput Struct，1993，49（5）：885-896.

［3］ SETHIAN J A，WIEGMANN A.Structural boundary design via level set and immersed interface methods［J］.J Comput Phys，1999，163（2）：489-528.

［4］ 隋允康，彭細(xì)榮.結(jié)構(gòu)拓?fù)鋬?yōu)化ICM 方法的改善［J］.力學(xué)學(xué)報(bào)，2005，37（2）：190-198.

［5］ SIGMUND O.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.

［6］ HABER R B，JOG C S.BENDSφE M P.A new approach to variable-topology design using a constraint on the per-imeter［J］.Structural and Multidisciplinary Optimization，1996，11（1）：1-12.

［7］ PETERSSON J，SIGMUND O.Slope constrained topology optimization［J］.International Journal for Numerical Methods in Engineering，1998，41（8）：1417-1434.

［8］ Sigmund O.Materials with prescribed constitutive parameters：An inverse homogenization problem［J］.International Journal of Solids and Structures，1994，31（17）：2313-2329.

［9］ MICHELL A G M.The limits of economy of material in frame structure［J］.Phil Mag，1904，8（6）：589-597.

［10］ ZHOU Ke-ming，LI Xia.Topology optimization of structures under multiple load cases using a fiber-reinforced composite material model［J］.Comput Mech，2005，38（2）：163-170.

［11］ ZHOU Ke-ming，LI Xia.Topology optimization for minimum compliance under multiple loads based on continuous distribution of members［J］.Struct Multidisc Optim，2008，37（1）：49-56.

［12］ XIE Yin-min，STEVEN G P.Optimal design of multiple load case structures using an evolutionary procedure［J］.Engineering Computations，1994，11（4）：295-302.

［13］ CHU D N，XIE Yin-min，HIRA A，et al.Evolutionary structural optimization for problems with stiffness constraints［J］.Finite Elements in Analysis and Design，1996，21（4）：239-251.