李霞，周克民

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

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

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

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

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

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

1.2 集中桿件的選擇

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

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

1.3 分布桿件的選擇

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

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

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

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

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

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

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

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

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

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

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

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

將式（7）代入式（6）的積分表達(dá)式，得

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

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

式（10）中：s為?。弧為徑向桿件的橫截面平均面積.

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

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

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

2 算法的實(shí)現(xiàn)

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

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

1）在奇異點(diǎn)附近，桿件向奇異點(diǎn)匯交，因此該區(qū)域的桿件變化較大.因?yàn)榧辛ψ饔命c(diǎn)附近會(huì)有應(yīng)力集中，計(jì)算誤差也會(huì)較大.根據(jù)奇異點(diǎn)的性質(zhì)，在奇異點(diǎn)附近的桿件方向做特別的處理.當(dāng)桿件進(jìn)入奇異點(diǎn)附近時(shí)，桿件方向一律直接與奇異點(diǎn)相連.進(jìn)入奇異點(diǎn)附近的標(biāo)準(zhǔn)是桿件與奇異點(diǎn)的距離小于單元邊長(zhǎng)的一半.

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

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

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

圖2 奇異點(diǎn)判斷Fig.2 Judgment of singular node

2.3 確定潛在桿件的過(guò)程

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

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

3）以該線段終點(diǎn)為下一線段的起點(diǎn)畫(huà)下一段線段，直到域邊界或密度過(guò)小，得到若干直線段連接而成的折線.

4）重新沿該折線逐段計(jì)算式，累計(jì)達(dá)到指定值標(biāo)記一個(gè)“結(jié)點(diǎn)”；再計(jì)算下一段，得到若干結(jié)點(diǎn).

5）分別從這些結(jié)點(diǎn)開(kāi)始，沿橫向畫(huà)另一組折線；

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

3 算例

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

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

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

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

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