一種新的考慮柔順度要求的結(jié)構(gòu)拓?fù)鋬?yōu)化方法

易繼軍 ，榮見(jiàn)華，曾韜

(1. 中南大學(xué) 機(jī)電工程學(xué)院，湖南 長(zhǎng)沙，410083；

2. 長(zhǎng)沙理工大學(xué) 汽車與機(jī)械工程學(xué)院，湖南 長(zhǎng)沙，410004)

連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化方法主要有均勻化方法[1-2](Homogenization)、實(shí)體各向同性材料懲罰法[3-4](Solid isotropic material with penalization, SIMP)、漸進(jìn)結(jié)構(gòu)優(yōu)化方法[5-6](Evolutionary structural optimization, ESO)、水平集方法[7-8](Level set method)、獨(dú)立連續(xù)拓?fù)渥兞考坝成渥儞Q方法(Independent continuous mapping,ICM)[9-10]和拓?fù)涿枋龊瘮?shù)法(Topology description function, TDF)[11]等。SIMP法及ICM方法將拓?fù)鋬?yōu)化問(wèn)題轉(zhuǎn)化為材料分布問(wèn)題來(lái)處理，因此，需要將最終求得的材料布局解釋為待求的結(jié)構(gòu)拓?fù)?。另外，采用ICM法和SIMP法，在迭代優(yōu)化過(guò)程中，被刪除的大量材料單元仍以較小的剛度和質(zhì)量存在于結(jié)構(gòu)中，一方面，易引起結(jié)構(gòu)總剛度矩陣病態(tài)以及結(jié)構(gòu)出現(xiàn)虛假的局部模態(tài)；另一方面，結(jié)構(gòu)分析的規(guī)模始終保持為最大設(shè)計(jì)域的規(guī)模，其分析計(jì)算量大。為了解決上述問(wèn)題，本文作者針對(duì)體積約束和柔順度最小的結(jié)構(gòu)拓?fù)鋬?yōu)化問(wèn)題，首先采用有理分式材料模型，建立結(jié)構(gòu)剛度矩陣與拓?fù)渥兞康年P(guān)系，然后，根據(jù)優(yōu)化準(zhǔn)則法(OC)[12]，提出一種考慮柔順度要求和基于設(shè)計(jì)空間調(diào)整的結(jié)構(gòu)拓?fù)鋬?yōu)化方法。

1 單元拓?fù)渥兞颗c過(guò)濾函數(shù)

1.1 單元拓?fù)渥兞?/h3>

類似于ICM方法，設(shè)置第i號(hào)單元的拓?fù)渥兞繛閠i。當(dāng)拓?fù)渥兞縯i=0時(shí)，表示該單元不存在；當(dāng)拓?fù)渥兞縯i=1時(shí)，表示該單元存在；當(dāng)拓?fù)渥兞?＜ti＜1時(shí)，表示該單元為從無(wú)到有的過(guò)渡狀態(tài)。

1.2 過(guò)濾函數(shù)

用過(guò)濾函數(shù) fv(ti)和 fk(ti)識(shí)別單元體積和單元?jiǎng)偠?，單元性質(zhì)參數(shù)識(shí)別采用如下公式：

其中：Vi和[Ki]分別表示拓?fù)渥兞繛閠i對(duì)應(yīng)的單元體積和單元?jiǎng)偠染仃?；分別表示單元固有體積和單元固有的剛度矩陣。這里取 fv( ti)=ti。

為了有利于設(shè)計(jì)空間減縮和擴(kuò)展的實(shí)現(xiàn)，過(guò)濾函數(shù)fk(ti)必須合理地設(shè)計(jì)。設(shè)計(jì)空間靈敏度反映新的設(shè)計(jì)變量增加對(duì)目標(biāo)函數(shù)或約束條件的影響[13]。因?yàn)樵O(shè)計(jì)空間的變化數(shù)學(xué)上是一個(gè)不連續(xù)的過(guò)程，Kim等采用增選主單元狀態(tài)和方向?qū)?shù)計(jì)算了設(shè)計(jì)空間靈敏度[14]。圖1所示為設(shè)計(jì)空間的一個(gè)增選主單元狀態(tài)，在結(jié)構(gòu)周邊引入一層接近于0的低密度材料單元，而這些單元的增添對(duì)目標(biāo)函數(shù)或約束條件影響很小。本文在設(shè)計(jì)空間擴(kuò)展時(shí)采用類似的措施(引入一層人工材料單元)，而在設(shè)計(jì)空間減縮時(shí)，將一些具有很小拓?fù)渥兞康牟牧蠁卧獜慕Y(jié)構(gòu)上去掉，這些單元的刪除對(duì)目標(biāo)函數(shù)或約束條件影響也是很小的。

圖1 增選主單元狀態(tài)概念Fig.1 Concept of pivot phase

設(shè)計(jì)空間靈敏度的定義，要求fk(ti)具有下列特性：

2 設(shè)計(jì)空間的縮減和擴(kuò)展

將結(jié)構(gòu)分為可設(shè)計(jì)和不可設(shè)計(jì)區(qū)域，設(shè)可設(shè)計(jì)區(qū)域的單元數(shù)為P，其單元編號(hào)可設(shè)為ip(p=1, 2, …, P)，可設(shè)計(jì)單元的拓?fù)渥兞?tip在迭代計(jì)算中在0到1之間變化；不可設(shè)計(jì)區(qū)域的單元數(shù)為Q，其單元編號(hào)可設(shè)為nq(q=1, 2, …, Q)，不可設(shè)計(jì)單元的拓?fù)渥兞?tnq=1在迭代計(jì)算中不變。

因此，將具有體積約束的最小柔順度的拓?fù)鋬?yōu)化問(wèn)題寫為：

其中：V0為整個(gè)設(shè)計(jì)域的初始體積； V為第 n個(gè)單元固有體積；θ(0＜θ＜1)為優(yōu)化目標(biāo)體積比；V為優(yōu)化后的結(jié)構(gòu)體積。類似于SIMP方法，這里設(shè)置非零的下限tip，以確保優(yōu)化過(guò)程中減縮有限單元模型規(guī)模后結(jié)構(gòu)仍保持非奇異。

另外，為了使迭代中的優(yōu)化拓?fù)溆休^好的0～1分布特征，同時(shí)兩相鄰迭代步的優(yōu)化拓?fù)渥兓^小，這里采用以體積約束限漸進(jìn)的方式，將模型(3)轉(zhuǎn)化為變體積約束限的系列模型(4)進(jìn)行求解。這種變體積約束限的求解模式使柔順度小量變化，又能使拓?fù)湓O(shè)計(jì)變量在某一小鄰域移動(dòng)，且具有較好的0～1分布特征。

為了易于求解具有較大規(guī)模有限元網(wǎng)格模型的優(yōu)化問(wèn)題，這里采用一些優(yōu)化迭代步將拓?fù)渥兞亢苄〉膯卧獜慕Y(jié)構(gòu)上去掉，而其他單元拓?fù)渥兞坎蛔兊牟呗?。即設(shè)置+ε(p=1, 2, …, P)，ε為一小量。把每個(gè)可設(shè)計(jì)的保留單元的拓?fù)渥兞?tip與閾值相比較來(lái)確定單元是否存在的狀態(tài)，基于方程(6)，挑選用于單元?jiǎng)h除的最潛在的候選單元，即將這些單元看成將要消失的單元。

為了確保優(yōu)化迭代中的結(jié)構(gòu)非奇異，使優(yōu)化迭代正常進(jìn)行，同時(shí)，結(jié)構(gòu)優(yōu)化時(shí)具有設(shè)計(jì)空間的擴(kuò)展(增添單元)功能。這里采用在結(jié)構(gòu)邊界和孔洞周圍附加人工材料單元的辦法來(lái)確保結(jié)構(gòu)剛度矩陣非奇異。人工材料單元附加的具體方法見(jiàn)文獻(xiàn)[5, 15]。這里將人工材料單元的彈性模量設(shè)為占據(jù)該空間的原始材料單元的彈性模量，僅將人工材料單元的拓?fù)渥兞咳閠ip，人工材料單元參與結(jié)構(gòu)參數(shù)計(jì)算和優(yōu)化迭代。根據(jù)人工材料單元拓?fù)渥兞康淖兓闆r，對(duì)設(shè)計(jì)空間進(jìn)行擴(kuò)展(即增添單元)。設(shè)增添單元的閾值為，基于方程(7)，挑選用于單元增添的最潛在的候選單元，即將這些單元看成將要增添的單元。

設(shè)計(jì)空間縮減和擴(kuò)展的準(zhǔn)則表為：

式中：nd和na分別為集合Nd和Na的單元個(gè)數(shù)；為小的經(jīng)驗(yàn)參數(shù)值。

若式(8)成立，則執(zhí)行式(10)和(11)的操作：

當(dāng)優(yōu)化迭代求解接近最優(yōu)結(jié)構(gòu)時(shí)，滿足式(6)的單元已經(jīng)很少，此時(shí)，適當(dāng)提高，對(duì)于滿足條件的單元執(zhí)行式(12)和(13)的操作：

3 靈敏度分析

在有限元分析中結(jié)構(gòu)靜平衡方程可表示成：

結(jié)構(gòu)柔順度可寫為：

對(duì)方程(15)求導(dǎo)得到：

因?yàn)閇K]-1是對(duì)稱矩陣，由方程(16)可得到：

因?yàn)樵O(shè)計(jì)變量 ti僅與第i號(hào)單元相關(guān)，因此，方程(17)可寫為：

考慮ti作為變量變化，第i號(hào)單元?jiǎng)偠染仃嚨膶?dǎo)數(shù)可表示為：

將方程(19)代入方程(18)得到：

體積對(duì)設(shè)計(jì)變量ti的導(dǎo)數(shù)為：

4 優(yōu)化模型的求解

拓?fù)鋬?yōu)化問(wèn)題式(4)的Lagrangian函數(shù)為：

當(dāng)t取極值t*時(shí)，上述Lagrangian函數(shù)滿足K-T條件：

式(23)又可寫為：

將方程(20)和(21)代入方程(24)中的第1式得：

即

將方程(26)作為優(yōu)化設(shè)計(jì)準(zhǔn)則，得到基于優(yōu)化準(zhǔn)則法的迭代公式如下：

式中：α為阻尼系數(shù)，引入α的目的是為了確保數(shù)值計(jì)算的穩(wěn)定性和收斂性。

5 拉格朗日乘子的計(jì)算

假設(shè)當(dāng)?shù)?jì)算由第k步更新到第k＋1步時(shí)，結(jié)構(gòu)體積由 Vk變化到 Vk+1。

由方程(26)，(27)和(28)整理得方程(29)，即可求解λ。

6 算例

6.1 一端面固支的短懸臂梁結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)

如圖2所示是一根短懸臂梁的經(jīng)典例子。左端固定，梁的尺寸是Lx=0.16 m，Ly=0.10 m，厚度t=0.02 m，在自由端面中部 20 mm長(zhǎng)區(qū)域鉛垂方向施加 F=50 kN/m的均布線載荷(為了避免應(yīng)力集中和計(jì)算結(jié)果與劃分網(wǎng)格的關(guān)聯(lián)度大，所有算例均采用均布載荷)；彈性模量 E＝210 GPa，泊松比ν＝0.3，密度為 7.8 t/m3。初始設(shè)計(jì)區(qū)域?yàn)?.16 m×0.10 m×0.02 m的立體區(qū)域，將其劃分為80×50×2即8 000個(gè)八節(jié)點(diǎn)六面體單元網(wǎng)格。為了評(píng)價(jià)結(jié)構(gòu)拓?fù)涞男阅?，Liang等[16]提出了性能指標(biāo)的概念。文獻(xiàn)[17]中拓?fù)涞男阅苤笜?biāo)PI=1.253 3[18]，而采用本文的方法得到拓?fù)?PI=1.275 5。這表明本文的方法得到拓?fù)浣Y(jié)構(gòu)性能比文獻(xiàn)[17]中的優(yōu)，同時(shí)它的計(jì)算效率更高。

圖2 短懸臂梁初始設(shè)計(jì)區(qū)域Fig.2 Initial design domain for a short cantilever

參數(shù)設(shè)置如下：優(yōu)化目標(biāo)體積比3.0=θ；體積約束限系數(shù)1.0=β；阻尼系數(shù)5.0=α；過(guò)濾函數(shù)系數(shù)ν = 4.5；拓?fù)渥兞肯孪拗祎i= 0 .001(i=1, 2, …, P)；拓?fù)渥兞孔兓y值η*=0.1。設(shè)置單元?jiǎng)h除閥值初始值為0.01，當(dāng)優(yōu)化結(jié)構(gòu)接近約束體積后提高到0.05，單元增加閥值

采用本文方法得到的懸臂梁拓?fù)浣Y(jié)構(gòu)演化歷程如圖3所示，共進(jìn)行77輪迭代。圖3(a)所示為第10輪迭代后得到的拓?fù)浣Y(jié)構(gòu)，結(jié)構(gòu)的體積為2.468 8×10-4mm3，柔順度為0.072 034 2；圖3(b)所示為第30輪迭代后得到的拓?fù)浣Y(jié)構(gòu)，結(jié)構(gòu)的體積為 1.350 4×10-4mm3，柔順度為0.136 74；圖3(c)所示為第60輪迭代后得到的拓?fù)浣Y(jié)構(gòu)，結(jié)構(gòu)的體積為1.0×10-4mm3，柔順度為0.133 01；圖3(d)所示為第77輪迭代后得到的拓?fù)浣Y(jié)構(gòu)，結(jié)構(gòu)的體積為 9.6×10-5mm3，柔順度為0.130 37。

圖3 本文方法得到的短懸臂梁結(jié)構(gòu)拓?fù)溲莼瘹v程Fig.3 Evolutionary history of a short cantilever structure obtained by proposed method

圖 4所示為結(jié)構(gòu)達(dá)到最優(yōu)時(shí)拓?fù)渥兞康姆植紶顟B(tài)。從圖4可以看出：結(jié)構(gòu)只存在少量的低密度單元，最優(yōu)結(jié)構(gòu)接近完全的0～1分布。

圖4 拓?fù)渥兞恐捣植紙DFig.4 Topology variable value distribution diagrams of topology structures

6.2 簡(jiǎn)支梁結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)

如圖5所示簡(jiǎn)支梁，一端固定鉸支座，另一端為可以橫向滾動(dòng)的鉸支座，梁的尺寸是Lx=0.16 m，Ly=0.03 m，厚度t=0.02 m，在梁的上邊緣中部20 mm長(zhǎng)的區(qū)域鉛垂方向施加F=50 kN/m的均布線載荷。彈性模量E=210 GPa，泊松比ν＝0.3，密度為7.8 t/m3。初始設(shè)計(jì)區(qū)域?yàn)?.16 m×0.03 m×0.02 m的立體區(qū)域，將其劃分為128×24×2即6 144個(gè)八節(jié)點(diǎn)六面體單元網(wǎng)格。

圖5 簡(jiǎn)支梁初始設(shè)計(jì)區(qū)域Fig.5 Initial design domain for a simply supported beam

采用本文方法得到的懸臂梁拓?fù)浣Y(jié)構(gòu)演化歷程如圖6所示。經(jīng)過(guò)154輪迭代步達(dá)到最優(yōu)拓?fù)?見(jiàn)圖6(d))。

圖6 本文方法得到的簡(jiǎn)支梁結(jié)構(gòu)拓?fù)溲莼瘹v程Fig.6 Evolutionary history of a simply supported beam structure obtained by proposed method

7 結(jié)論

(1) 借鑒有理分式材料模型, 建立了結(jié)構(gòu)剛度矩陣與拓?fù)渥兞康年P(guān)系，給出了不影響光滑求解算法收斂特性的設(shè)計(jì)空間調(diào)整策略和手段。

(2) 根據(jù)優(yōu)化準(zhǔn)則法，結(jié)合變體積約束限措施，提出了一種考慮柔順度要求和基于設(shè)計(jì)空間調(diào)整的結(jié)構(gòu)拓?fù)鋬?yōu)化方法。

(3) 算例表明該方法能獲得較好 0～1分布特征的優(yōu)化拓?fù)?，能解決材料分布和分析計(jì)算量大的問(wèn)題。

[1] Cheng G D, Olhoff N. An investigation eoncerning optimal design of solid elastic plates[J]. International Journal of Solids and Structures, 1981, 17(3): 305-323.

[2] 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.

[3] Belblidia F, Lee J E B, Rechak S, et al. Topology optimization of plate structures using a single- or three-layered artificial material model[J]. Advances in Engineering Software, 2001, 32(2):159-168.

[4] Giuseppe C A. Hierarchical solution of large-scale threedimensional topology optimization problems[C]//The 1996 ASME Design Engineering Technical and Computer in Engineering Conference Miechigan. Miechigan State University,1996.

[5] 榮見(jiàn)華, 唐國(guó)金, 羅銀燕, 等. 考慮位移要求的大型三維連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化數(shù)值技術(shù)研究[J]. 工程力學(xué), 2007, 24(3):20-27.RONG Jian-hua, TANG Guo-jin, LUO Yin-yan, et al. A research on the numerical topology optimization technology of large three-dimensional continuum structural considering displacement requirements[J]. Engineering Mechanics, 2007,24(3): 20-27.

[6] 榮見(jiàn)華, 姜節(jié)勝, 胡德文, 等. 基于應(yīng)力及其靈敏度的結(jié)構(gòu)拓?fù)錆u進(jìn)優(yōu)化方法[J]. 力學(xué)學(xué)報(bào), 2003, 35(5): 584-591.RONG Jian-hua, JIANG Jie-sheng, HU De-wen. A structural topology evolutionary optimization method based on stresses and their sensitivity[J]. Acta Mechanica Sinica, 2003, 35(5):584-591.

[7] Wang X, Wang M Y, Guo D. Structural shape and topology optimization in a level-set-based framework of region representation[J]. Structural and Multidisciplinary Optimization,2004, 27(1/2): 1-19.

[8] 歐陽(yáng)高飛, 張憲民. 結(jié)構(gòu)拓?fù)鋬?yōu)化中水平集方法的快速初始化研究[J]. 應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào), 2008, 16(1): 136-143.OUYANG Gao-fei, ZHANG Xian-min. Research on fast initialization for the level set method in structural topology optimization[J]. Journal of Basic Science and Engineering, 2008,16(1): 136-143.

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

[10] 隋允康, 葉紅玲, 劉建信, 等. 追究根基的結(jié)構(gòu)拓?fù)鋬?yōu)化方法[J]. 工程力學(xué), 2008, 25(A02): 7-19.SUI Yun-kang, YE Hong-ling, LIU Jian-xin, et al. A structural topological optimization method based on exploring conceptual root[J]. Engineering Mechanics, 2008, 25(A02): 7-19.

[11] 郭旭, 趙康. 基于拓?fù)涿枋龊瘮?shù)的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化方法[J]. 力學(xué)學(xué)報(bào), 2004, 36(5): 522-526.GUO Xu, ZHAO Kang. A new topology description function based approach for structural topology optimization[J]. Acta Mechanica Sinica, 2004, 36(5): 522-526.

[12] 程耿東. 工程結(jié)構(gòu)優(yōu)化設(shè)計(jì)基礎(chǔ)[M]. 北京: 水利電力出版社,1984: 161-170.CHENG Geng-dong. The base of optimal design of engineering structure[M]. Beijing: The Publishing Company of Water and Electricity, 1984: 161-170.

[13] Jang I G, Kwak B M. Design space optimization using design space adjustment and refinement[J]. Structural and Multidisciplinary Optimization. 2008, 35(1): 41-54.

[14] Kim I Y, Kwak B M. Design space optimization using a numerical design continuation method[J]. Int J Numer Methods Eng, 2002, 53: 1979-2002.

[15] Stolpe M, Svanberg K. An alternative interpolation scheme for minimum compliance topology optimization[J]. Structural and Multidisciplinary Optimization, 2001, 22(2): 116-124.

[16] Liang Q Q, Steven G P. A performance-based optimization method for topology design of continuum structures with mean compliance constraints[J]. Computer Methods in Applied Mechanics and Engineering, 2002, 191(13/14): 1471-1489.

[17] Suzuki K, Kikuchi N. A homogenization method for shape and topology optimization[J]. Computer Methods in Applied Mechanics and Engineering, 1991, 93(3): 291-318.

[18] RONG Jian-hua, LIANG Qing-quan. A level set method for topology optimization of continuum structures with bounded design domains[J]. Computer Methods in Applied Mechanics and Engineering, 2008, 197(17/18): 1447-1465.