應用類桁架模型的連續(xù)體拓撲優(yōu)化方法

鄭偉偉，周克民

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

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

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

1.1 彈性矩陣

在類桁架連續(xù)體材料模型中，材料模型是由無限多的非均勻連續(xù)分布的桿件構成，假設任意點都是由兩組正交的桿件構成，它們在任意點的密度用t1和t2表示.假設應力σi和應變εi的線性關系為

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

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

在這種假設下，當t1＝t2時，式（2）的各項為同性材料.若假設結構坐標軸到材料坐標軸的轉角為α，那么結構坐標系中的彈性矩陣可由坐標轉換矩陣獲得，即

式（3）中：T（α）為坐標轉換矩陣，即

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

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

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

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

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

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

1.2 剛度矩陣

將單元剛度矩陣（8）帶入剛度矩陣，有

2 優(yōu)化方法

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

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

迭代過程有如下5個步驟.

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

上標0代表迭代次數(shù)，初始密度剔除標準設為0；

步驟2 通過有限元方法，計算節(jié)點處的主應力方向和主應力方向上的應變的大??；

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

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

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

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

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

3 數(shù)值算例

運用文中的算法，計算兩個經典的算例，兩個算例中：彈性模量E＝210 GPa；允許應力σp＝160 MPa.

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

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

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

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

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

第2個算例是跨中受集中荷載的簡支梁，尺寸如圖5所示.圖5中：厚度（h）為0.01 m；設計域被矩形單元劃分為60×40.

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

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

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

圖5 簡支梁尺寸Fig.5 Dimension of simply supported beam

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

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

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

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

4 結束語

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

［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］ 隋允康，彭細榮.結構拓撲優(yōu)化ICM 方法的改善［J］.力學學報，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.