李霞，周克民

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

Michell桁架是精確的理論解，經(jīng)常被用來檢驗各種數(shù)值優(yōu)化方法的正確性，但由于沒有一般的求解方法，求解困難［1-3］.因此，許多數(shù)值方法大多采用有限元數(shù)值分析方法［4-7］.為了克服數(shù)值不穩(wěn)定問題，陸續(xù)提出了周長控制、梯度控制等方法［8］.這些方法不僅增加了計算量，而且計算過程中的一些控制參數(shù)事先難以估計，不適當(dāng)?shù)膮?shù)可能會得不到有意義的結(jié)果.Michell理論已經(jīng)揭示了拓撲優(yōu)化結(jié)構(gòu)的類桁架性質(zhì)，拓撲優(yōu)化結(jié)構(gòu)理論上一般是非均質(zhì)各向異性連續(xù)體.因此，上述優(yōu)化方法所采用的各向同性材料無法精確描述這種拓撲優(yōu)化結(jié)構(gòu).一些學(xué)者［9-11］將問題完全放松，采用一般各向異性材料模型.但是這種材料與工程結(jié)構(gòu)沒有明確的對應(yīng)關(guān)系，后期處理困難，而且也沒有反映類桁架結(jié)構(gòu)的本質(zhì).本文提出有限元優(yōu)化類桁架連續(xù)體方法，解決了求解困難［12-15］.

1 類桁架連續(xù)體離散化

1.1 類桁架連續(xù)體中桿件的分布性質(zhì)

按照Michell理論，拓撲優(yōu)化結(jié)構(gòu)是由無限細、無限密桿件構(gòu)成的類桁架連續(xù)體.習(xí)慣上將桿件在單工況應(yīng)力約束下的優(yōu)化分布區(qū)域分為5種：單向拉伸（R+）、單向壓縮（R-）、各向均勻拉伸（S+）、各向均勻壓縮（S-）和兩向分別拉壓（T）.如果不區(qū)分拉壓，去掉上角標中的符號，可歸納為3類：單向拉壓（R，桿件沿拉壓方向），各向均勻拉壓（S，桿件沿任意方向）和兩向分別拉壓（T，桿件沿拉壓兩個方向）.對于多工況或其他約束，桿件優(yōu)化分布性質(zhì)有所不同.將以上桿件優(yōu)化分布劃分方式推廣到更一般的情況，其中：S區(qū)域仍然表示桿件任意分布；R區(qū)域表示桿件沿某個單一方向分布；而T區(qū)域表示桿件沿某幾個確定方向（不限于2個，也不一定正交）分布.這種劃分方式可以包括桿件所有分布情況.S區(qū)域由于桿件優(yōu)化方向是任意分布的，所以不需要研究其優(yōu)化方向，實際優(yōu)化問題中也不常見.但是，實際優(yōu)化問題中會經(jīng)常遇到S區(qū)域退化為1個孤立的點的特殊情況，這是一個比較特殊的情況.例如當(dāng)許多桿件匯交于一點時就會出現(xiàn)這種情況，文中將這樣的點稱為“奇異點”.

在優(yōu)化的桿件分布場中存在分布桿件和集中桿件兩種情況.分布桿有無限多，不可能都保留.但集中桿數(shù)量有限，全部是平衡必需的，應(yīng)該全部保留.

1.2 集中桿件的選擇

一個力學(xué)問題有力和位移兩種邊界條件，在有限元計算中，各種荷載一般都要等效到結(jié)點上，形成結(jié)點集中力.在有限元分析完成后，位移邊界條件也可以轉(zhuǎn)化為結(jié)點集中力.一個優(yōu)化問題實際上就是設(shè)計這些集中力的優(yōu)化傳遞路徑.因此，文中只討論結(jié)點集中力的情況.結(jié)點集中力需要集中桿件傳遞，或者需要無限多匯交于一點的分布桿件.當(dāng)無限多的分布桿件匯交于一點時，桿件在匯交點位置的方向不確定，該點就是奇異點（S區(qū)）.因此，集中力作用點（包括位移約束結(jié)點，以下不再特別說明）的位置必然有集中桿件或是奇異點.從平衡的角度看，由于集中桿件也可以理解為集中力，所以集中桿件的端部必然是集中力或奇異點，選擇所有集中力作用點和奇異點作為集中桿可能經(jīng)過的位置，在這些位置上布置桿件就可以將所有集中桿件選擇上.

確定奇異點位置的桿件方向比較困難，需要進一步分析.分布桿的奇異點位置的桿件方向不確定，由于集中桿的兩端只能是集中力或奇異點，而集中力和奇異點的數(shù)量是有限的，所以為了將兩個端點都位于奇異點位置的集中桿件也選上，在每個奇異點位置分別選擇指向其他每個奇異點方向的桿件.

1.3 分布桿件的選擇

無論是集中桿件還是分布桿件，由于僅考慮軸力而不考慮彎矩作用，曲桿僅在桿端力的作用下不能平衡，必須借助橫向分布桿件.因此，在曲桿的橫向應(yīng)該有分布桿件.

在形成桿系結(jié)構(gòu)過程中，分布桿件選擇要適當(dāng).桿件選多了可以減少誤差，更接近理論解，但過多的桿件并不實用.選擇的標準是用最少的桿件實現(xiàn)最小的誤差.

分布桿與集中桿結(jié)構(gòu)的體積誤差，如圖1所示.圖1（a）中：曲線AB表示一段集中曲桿；曲率半徑為R；圓心角為2α；圓心為O；A，B兩端的集中力F分別為沿圓弧在A，B點的切向.如果曲桿AB不受彎矩作用，在扇形ABO區(qū)域內(nèi)應(yīng)有徑向分布桿，即Michell桁架的一個扇形段，使得圓弧桿的橫向受到分布應(yīng)力σ的作用，它傳遞A，B和O3點集中力在理論上的最優(yōu)傳力路徑［13］.根據(jù)圓弧桿AB的豎向平衡條件，可以知道這些分布桿件的合力作用在O點，大小為2Fsinα.假設(shè)這個扇形區(qū)域內(nèi)有均勻徑向應(yīng)變ε，材料的彈性模量為E，那么，所需材料體積為Vm＝（4FR／Eε）α.

1.4 離散桿件系統(tǒng)體積誤差估計

采用圖1（b）所示的三角形離散桁架來代替圖1（a）所示的理論上的分布桿，由平衡關(guān)系，圖1（b）中各桿件的軸力為

圖1 分布桿與集中桿結(jié)構(gòu)的體積誤差Fig.1 Volume error between distributed and concentrated members

在同樣的應(yīng)變場下，三角形離散桁架的體積為

由式（2）減式（1），得到兩個結(jié)構(gòu)的體積誤差，即

如果離散結(jié)構(gòu)的桿件較多，其圓心角不會太大，可以利用三角函數(shù)的展開式，即

將式（4）代入式（3），得

由圖1（a）的平衡關(guān)系，得

式（6）中：t是徑向桿件在環(huán)向截面上的桿件密度.作為近似估計，假設(shè)密度是線性變化，則

將式（7）代入式（6）的積分表達式，得

再將式（8）代入式（6）中，得

最后，將式（9）代入式（5）中，得

式（10）中：s為??；ˉt為徑向桿件的橫截面平均面積.

結(jié)點i位置的徑向桿件在環(huán)向截面的密度記作為ti，角度為αi.結(jié)點i與結(jié)點i-1之間的單元邊界長為si.結(jié)點1到k之間徑向桿件的平均環(huán)向截面面積的計算式為

將式（11）代入式（10）中，得

由離散桿系代替連續(xù)分布桿件引起的體積增加量（ΔV）作為控制條件，決定離散桿件的選擇.

2 算法的實現(xiàn)

2.1 離散桿系結(jié)構(gòu)

由于集中桿件的端部一定是集中力或奇異點，所以在每個集中力作用點布置桿件，并且所有奇異點之間連接桿件就可以保證所有集中桿件都被選擇上.由于分布桿件無限多，不能全部保留，只能保留相距一定間距的部分桿件.將其余的分布桿件集中到被保留的桿件上.確定桿件間距的依據(jù)就是使式（12）為常數(shù).此外，還有2點需要說明.

1）在奇異點附近，桿件向奇異點匯交，因此該區(qū)域的桿件變化較大.因為集中力作用點附近會有應(yīng)力集中，計算誤差也會較大.根據(jù)奇異點的性質(zhì)，在奇異點附近的桿件方向做特別的處理.當(dāng)桿件進入奇異點附近時，桿件方向一律直接與奇異點相連.進入奇異點附近的標準是桿件與奇異點的距離小于單元邊長的一半.

2）由于存在數(shù)值計算誤差，沒有材料部分的桿件密度會比0大一些.特別是對于位移約束結(jié)點，設(shè)定一個閥值，當(dāng)密度低于閥值時認為沒有桿件了.

2.2 奇異結(jié)點的判斷

在形成離散結(jié)構(gòu)過程中，判斷奇異點是一個比較困難的問題.從理論上講，桿件匯交點就是奇異點.但是，由于數(shù)值計算誤差的存在，計算結(jié)果并不能保證相鄰幾個結(jié)點的桿件恰好相交于奇異點.奇異點判斷，如圖2所示.圖2中：假設(shè)結(jié)點i是奇異點；圍繞結(jié)點i的所有結(jié)點（稱為相鄰結(jié)點）為j，j∈Ji；相鄰結(jié)點與結(jié)點i同屬于一個單元，包括沒有邊聯(lián)接的結(jié)點2；密度足夠大的桿件都匯交結(jié)點i.奇異點是匯交點，所以密度應(yīng)該比周圍的密度大.結(jié)點j密度最大方向桿件所在直線到結(jié)點i的距離應(yīng)足夠小.兩個連續(xù)的相鄰點和奇異點就構(gòu)成屬于S區(qū)域的一個子域，為了避免誤判，要求至少存在2個這樣的連續(xù)子域.

圖2 奇異點判斷Fig.2 Judgment of singular node

2.3 確定潛在桿件的過程

1）選擇所有集中力作用結(jié)點包括位移約束結(jié)點，沿著桿件密度足夠大的方向畫線段，交于單元的邊界作為線段終點.由該線段終點所在單元邊界兩端的結(jié)點處的桿件密度和方向，插值得到線段終點處的桿件密度和方向.

2）為了提高計算精度，取線段起點和終點的桿件方向的平均值，從起點重新畫該直線，得到修正的線段.

3）以該線段終點為下一線段的起點畫下一段線段，直到域邊界或密度過小，得到若干直線段連接而成的折線.

4）重新沿該折線逐段計算式，累計達到指定值標記一個“結(jié)點”；再計算下一段，得到若干結(jié)點.

5）分別從這些結(jié)點開始，沿橫向畫另一組折線；

6）重復(fù)過程1）～3），得到一個曲線網(wǎng)絡(luò).曲線網(wǎng)絡(luò)相交點作為結(jié)點，用直線代替兩個結(jié)點之間的折線構(gòu)成離散桁架.

3 算例

力學(xué)模型，如圖3所示.圖3中：矩形設(shè)計域左邊固定；上邊作用兩個集中力.按照滿應(yīng)力準則，任意位置的應(yīng)變不超過允許應(yīng)變.當(dāng)離散后的結(jié)構(gòu)桿件足夠多，應(yīng)變差異應(yīng)該不大由于拓撲優(yōu)化結(jié)果與力和尺寸大小無關(guān).集中力（n）的大小分別取-2，-1，0，2，當(dāng)n從0到2之間變化時，拓撲優(yōu)化結(jié)果差別不大.因此，沒有給出n＝1的結(jié)果.采用40×20矩形單元，優(yōu)化的桿件分布場和離散桁架，如圖4所示.圖4中：左列給出桿件的優(yōu)化分布；線段長度表示桿件密度；線段方向表示桿件方向；右列給出對應(yīng)的拓撲優(yōu)化離散桿系結(jié)構(gòu).離散桿系結(jié)構(gòu)中的桿件數(shù)量可以根據(jù)需要選擇，粗線表示壓桿，細線表示拉桿.結(jié)果與理論解［2，4］比較接近.

圖3 力學(xué)模型Fig.3 Mechanics model

圖4 優(yōu)化的桿件分布場和離散桁架Fig.4 Optimum truss-like continua and their discrete truss

4 結(jié)束語

研究了基于類桁架材料模型的桿系結(jié)構(gòu)拓撲優(yōu)化方法.類桁架優(yōu)化過程中沒有抑制中間密度，完全避免了數(shù)值不穩(wěn)定問題.桿系結(jié)構(gòu)通過選擇類桁架中的部分桿件形成，桿件的數(shù)量可以直觀地控制，從而得到滿足工程需要的結(jié)構(gòu).將桁架結(jié)果按照截面相等的原則轉(zhuǎn)化為等厚帶孔板，結(jié)果同樣會比均勻化等以單元表示結(jié)構(gòu)拓撲的方法更精確，效率更高，具體實現(xiàn)方法是下一步要進行的工作.

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

［2］ LEWINSKI T，ZHOU M，ROZVANY G I N.Extended exact solutions for least-weight truss layouts（Part I）：Cantilever with a horizontal axis of symmetry［J］.International Journal of Mechanical Sciences，1994，36（5）：375-398.

［3］ ESCHENAUER H A，OLHOFF N.Topology optimization of continuum structures：A review［J］.Applied Mechanics Reviews，2001，54（4）：331-389.

［4］ ROZVANY G I N.Exact analytical solutions for some popular benchmark problems in topology optimization［J］.Structural Optimization，1998，15（1）：42-48.

［5］ ROZVANY G I N，BENDSOE M P.KIRSCH U.Layout optimization of structures［J］.Applied Mechanics Reviews，1995，48（2）：41-119.

［6］ BENDSOE M P，KIKUCHI N.Generating optimal topologies in structural design using a homogenization method［J］.Computer Methods in Applied Mechanics and Engineering，1988，71（2）：197-224.

［7］ XIE Yin-min，STEVEN G P.A simple evolutionary procedure for structures optimization［J］.Computers and Structures，1993，49（5）：885-896.

［8］ SIGMUND O，PETERSSON J.Numerical instabilities in topology optimization：A survey on procedures dealing with checkerboards，mesh-dependencies and local minima［J］.Structural Optimization，1998，16（1）：68-75.

［9］ HORNLEIN H R E M，KOCVARA M，WERNER R.Material optimization：bridging the gap between conceptual and preliminary design［J］.Aerospace Science and Technology，2001，5（8）：541-554.

［10］ HSU Y，SHO M，CHEN Chuan-tang.Interpreting results from topology optimization using density contours［J］.Computers and Structures，2001，79（10）：1049-1058.

［11］ MATSUI K，TERADA K.Continuous approximation of material distribution for topology optimization［J］.International Journal for Numerical Methods in Engineering，2004，59（14）：1925-1944.

［12］ ZHOU Ke-min，LI Xia.Topology optimization for minimum compliance using fiber-reinforced composite material model［J］.Structural Optimization，2006，37（1）：49-56.

［13］ ZHOU Ke-min，LI Jun-feng.Topology optimization of structures under multiple load cases using fiber-reinforced composite［J］.Structural Optimization，2008，38（5）：525-532.

［14］ PRAGER W.A note on discretized Michell structures［J］.Computer Methods in Applied Mechanics and Engineering，1974，3（3）：349-355.

［15］ ZHOU Ke-min，LI Jun-feng.The exact weight of discretized Michell trusses for a central point load［J］.Structural Optimization，2004，28（1）：69-72.