基于種群動力學模型的連續體結構拓撲優化仿生方法

開依沙爾·熱合曼，熱合買提江·依明江 ，買買提明·艾尼

(1.新疆大學 機械工程學院，烏魯木齊 830047; 2.新疆大學 數學與系統科學學院，烏魯木齊 830046)

連續體結構拓撲優化問題因數學模型建立困難、所設計的變量較多,計算量大而被公認為是當前結構優化領域內具有挑戰性的課題之一 。自從Bends?e等[1]奠基性論文提出基于均勻化理論的結構拓撲優化設計方法開創了連續體結構拓撲優化研究的新篇章。在過去的幾十年里，連續體結構拓撲優化被廣泛地應用于機械工程、建筑工程等工程設計領域。目前連續體結構拓撲優化方法主要有均勻化方法[1]、變密度方法[2]、實體各向同性材料懲罰法[3](Solid Isotropic Material with Penalization, SIMP)、漸進結構優化方法(Evolutionary Structural Optimization, ESO)[4]、雙向漸進結構優化方法(Bi-directional Evolutionary Structural Optimization, BESO)[4]、水平集方法(Level set method)[5]以及獨立連續拓撲變量及映射變換方法 (Independent Continuous Mapping, ICM)[6]等。有關連續體結構拓撲優化方面的綜述可見文獻[7-8]。

上面所述的拓撲優化問題均來自工程領域，一般都用最優化方法來解決。而自然界中如:樹干中的不同木質分布、動物骨骼中材料的分布等,同樣存在自優化現象,值得借鑒。因此本文著重考慮尋求骨骼重建機理為手段，種群動力學模型為基礎，連續體結構為研究目的的拓撲優化仿生方法。

骨骼在生物力學環境下通過新陳代謝成骨、吸骨、再成形的方式優化自身的拓撲結構，創造出各種骨骼體系的動態穩定的結構和拓撲形狀。骨骼體系的這種動態平衡是通過骨重建并隨著成骨細胞和破骨細胞再吸過程來決定骨形成過程中的微妙平衡[9]。Wolff[10]曾指出骨骼具有功能適應能力，它受到外界載荷作用時會通過改變自身的物質分布以最少的材料獲得最大的結構強度。2000年Huiskes等[11]在Nature上發表論文提出了考慮骨骼應力信息及其成骨細胞和破骨細胞耦合反饋的骨骼功能適應性模型，并指出在合適的外載荷條件下，骨細觀結構維持最佳拓撲結構的結論。Tezuka等[12]在基于Turing模型的基礎上建立了iBone模型并骨骼宏觀結構進行數值模擬得到與真實骨骼較相似的結果。Bagge[13]從Wolff 法則出發，以骨的力學調控性為基礎，提出了骨結構的適應性功能可作為結構拓撲優化目標的重要結論。Zhu等[14]拓撲優化方法引入骨生物力學，架起了生物學與工程學的橋梁。Mattheck[15]根據樹木生長的自優化現象通過有限元方法提出了計算機輔助形狀優化方法。Andr’es[16]根據骨重建理論提出了混合元胞自動機拓撲優化方法。Cai等[17]將骨重建過程比擬為待優化連續體對連續體結構進行了拓撲優化。Kaveh等[18]將蟻群的方法引入到結構優化設計提出了仿生拓撲優化方法- 蟻群優化方法。Lee等[19]將蜂群算法引入到結構優化設計中，提出了人工蜂群拓撲優化算法。買買提明．艾尼等[20-21]采用有限元法耦合反應-擴散骨重建模型對骨宏觀微觀結構進行數值模擬。開依沙爾.熱合曼等[22]從骨骼的功能適應性出發，采用交叉型反應擴散模型提出了骨骼功能適應性模型，并對連續體結構進行了拓撲優化計算。

本文從骨骼功能自適應性法則出發，通過種群動力學模型和有限元方法的耦合建立了骨骼重建的數學模型。該模型通過單元添加(骨材料形成)和單元刪除(骨材料被吸收)準則提出了連續體結構具有仿生骨骼功能的拓撲優化計算方法。

本文在如下四個部分展開研究。第一部分，建立了以種群動力學模型為基礎的骨骼重建假想模型；第二部分，通過有限元方法和像素單元的添加與刪除準則對連續體結構提出了仿生拓撲優化計算方法；第三部分，通過兩個典型的數值算例來驗證了本文方法的有效性；第四部分，給出了三種不同邊界條件下對長懸臂梁模型和四種不同邊界條件下拱橋模型，并它們進行拓撲優化計算。最后給出了文章的結論。

1 基于種群動力學模型的骨重建模型

20世紀初，Lotka 和 Volterra 運用種群動力學方法解決生物種群之間的相互演化關系，引起人們的關注。后來人們發現生物進化論的規律、細胞的相互作用，以及細胞的增長規律，都可以用種群動力學的方法來描述。種群動力學從物理學、力學到生命科學等方面得到很好的應用[23]。兩種群在一個共同的自然環境中生存，他們之間的相互作用的最簡單的模型可以寫成如下方程式(1)。

(1)

其中:x和y分別表示兩種群的相對密度，ci為常數(i=1,2,…,6)。當c2≤0，c4≤0時，模型(1)成為強競爭系統[23]。

圖1 基于種群動力學模型的骨骼重建假想模型

本文中將利用種群動力學模型(1)中的兩種群假設為骨骼形成活化因子和骨骼吸收抑制因子，對骨重建做如下兩種假設(見圖1)：骨的成形及吸收與活化因子(Activator)和抑制因子(Inhibitor)各自的局部濃度成比例；反之局部有效應力又影響活化因子的濃度，當活化因子的局部濃度大于抑制因子時，大量成骨細胞聚集在這局部領域形成新的骨，當抑制因子的局部濃度大于活化因子時，破骨細胞聚集骨被吸收。本文從Wolff的骨骼功能適應原則出發,啟發Huiskes等人的工作，將種群動力學模型(1)中引入應力傳感器σi建立了式(2)所示的骨重建數學模型

(2)

其中:Ai和Ii表示第i單元中活性因子和抑制因子的局部濃度ai(i=1,2,…,5)為活性因子和抑制因子反饋參數，aS為應力補償參數。當有外力作用于這種模型時，受到激勵的骨骼通過沉積和吸收的方式使這個模型的形狀很快地適應外力。方程參數ai和aS取為沒有應力項aSδi下使活性因子和抑制因子濃度均勻分布的一族參數:a1=0.04，a2=-0.08，a3=0.5，a4=0.06，a5=0.01，aS=0.02

因為a2≤0，a4≤0，所以模型(2)為構成骨骼形成活化因子和骨骼吸收抑制因子的強競爭系統。

2 連續體結構拓撲優化仿生方法

為了實現連續體結構的拓撲優化仿生方法，把骨骼的形成與被吸收功能轉化為材料的形成和材料的被吸收過程，對連續結構提出了拓撲優化仿生方法，因為骨重建過程就是骨材料密度經過適當分布達到應變能密度均勻分布狀態的過程,因此在優化方法中取應變能密度的均勻分布作為優化準則來更新材料分布，一直到總應變能趨于穩定為止，反復迭代基于種群動力學骨重建模型公式(2)。實際拓撲優化計算時作者編寫了基于種群動力學模型的仿生拓撲優化軟件，軟件由前處理，主處理及后處理三個模塊組成的。前處理和后處理部分由實驗室自開發的基于波形法的有限元前后處理軟件(Finite element Analysis Support Tools with Wave Mesh Generator: FAST_wmg)來完成[24]。包括建立有限元初始網格模型，設定相應的邊界條件，材料屬性，結果可視化等。 主處理模塊是在Linux 平臺下用C++語言編寫了8個節點有限元分析模塊程序和基于種群動力學模型的拓撲優化計算程序來完成的。

{δk}=[D]k[B]k[uk]

(3)

其中:[D]k和[B]k分別表示第k單元的彈性模量矩陣和幾何矩陣。對于靜力學問題，考慮各向同性材料，Von Mises應力 一直是最常用的準則之一。它定義為：

(4)

其中:δxx和δyy為x和y方向的正應力，τxy為剪應力。

應變能密度為

(5)

其中:E表示表面彈性模量，USED表面應變能密度。

(1) 根據給定的初始設計區域，通過FAST_wmg軟件來建立有限元初始網格模型,給定相應邊界條件和材料的屬性參數；

(2) 對有限元初始網格模型進行有限元分析，計算各個節點上的有效應力σi的近似值；

(3) 通過骨重建模型(2)來計算各個節點上的A和I，并(Ai-bIi)的正負來判斷決定添加或刪除的單元；

(4) 根據添加和刪除的單元，決定新的模型；

(5) 判斷模型形狀是否收斂；

(6) 如果不收斂用新模型繼續到第二步并進行優化計算，直到應變能密度的均勻分布為止。

圖2 拓撲優化計算過程流程圖

圖3 材料形成與吸收及像素單元的添加和刪除準則

3 拓撲優化仿生方法的驗證

本文通過兩個連續體結構拓撲優化設計中常用的經典算例進行拓撲優化計算來驗證本文拓撲優化仿生方法的有效性，并且將其拓撲優化結果與其它被廣泛應用的幾種拓撲優化方法結果進行比較。優化問題可被描述為:在滿足材料體積約束下使結構總應變能最小化。算例中材料參數取為：彈性模量70 GPa、泊松比0.3。模型創建和邊界條件設定和結果的顯式由Fast_wmg軟件[24]來完成，所有算例由基于種群動力學模型的仿生拓撲優化軟件來完成。

算例1：短懸臂梁拓撲優化設計

圖4(a)為初始設計區域為1.6L×L的短懸臂梁。結構左邊界上下角上的兩個端點固定，右邊界中點作用于集中載荷，其值為F=100 N。結構邊長度L=100 mm。初始設計區域被劃分為160×100的pixel單元，如圖4(b)所示。每個pixel單元由1×1的正方形單元組成。為了檢驗設計域上孔洞的形狀對優化結果的影響，分別設計域上設計了均布的圓形孔洞和方形孔洞(直徑均為11 mm)如圖4(c)和圖4(e)所示,體積比V設定為45%。圓形孔洞有限元初始網格模型和方形孔洞有限元初始網格模型，通過本文拓撲優化仿生方法計算,經過40次循環后的優化結果如圖4(d)和圖4(f)所示。從圖4(d)可看出圓形孔洞所得到的最終拓撲結果，分別蟻群拓撲優化[18](圖5(b)), BESO方法[4](圖5(c)), SIMP方法[3]( 圖5(d)),變密度方法[2](圖5(e))和水平及方法[5](圖5(f))的結果較相似。從圖4(f) 可看出方形孔洞得到的最終拓撲結果經典的拓撲優化方法-均勻化方法[1](圖5(a))的結果較相似。從圖4(d)和圖4(f)可看出孔洞的形狀對最終拓撲優化結果有一定的影響,但是無論是圓形孔洞還是方形孔洞所得到的最終拓撲結果具有一定的規則性和對稱性。為了檢查本文方法對網格是否有依賴性，對初始區域分別被劃分為40×25pixel單元和80×50pixel單元,并對它們進行優化計算，40×25pixel單元的網格模型經過12次循環后的優化結果如圖6(a)所示，而80×50pixel單元的網格模型經過27次循環后的優化結果如圖6(b)所示。從圖6可看出本文方法幾乎無網格依賴性現象，但是初始設計區域上的網格對計算迭代次數和最終拓撲結果的光滑性有一定的影響。為了檢查體積比對拓撲優化結果的影響，體積比V分別設定為40%和30%。圖7為體積比為40%和30%時圓形孔洞對應的最終拓撲結果圖。圖8為體積比為40%和30%時方形孔洞對應的最終拓撲結果圖。從圖7和圖8可看出圓形孔洞和方形孔洞體積比為40%和30%時的最終拓撲形狀與其它傳統的拓撲優化方法的結果也較相似。因此體積比設定的大小對拓撲優化結果幾乎無影響。圖9為體積比V=45%時圓形孔洞和方形孔洞對應的應變能密度和體積比,隨著迭代次數的演化歷程,從圖9可看出首先方孔對應的應變能密度和體積比，比圓孔對應的應變能密度和體積比略大一些，但是通過12步迭代后保持持平，之后的迭代步驟應變能密度和體積比,逐步下降, 35～40步驟之后基本保持不變,都能收斂到體積比45%。

圖4 本文所得到的短懸臂梁模型的拓撲優化結果

圖5 其它方法獲得拓撲優化結果

圖6 不同網格節點下的拓撲優化結果

圖7 不同的體積下圓形孔洞的最終拓撲結果

圖8 不同的體積下方形孔洞的最終拓撲結果

圖9 短懸臂梁的應變能密度與體積比,隨迭代次數的演化歷程

算例2:簡支的Michell型結構拓撲優化設計

圖10(a)為初始設計區域2L×L的一端鉸支、一端滾珠支承的傳統Michell型結構，結構下邊緣中部鉛垂方向施加F=100 N的集中載荷。結構邊長度L=100 mm。初始設計區域被劃分為200×100的pixel單元，如圖10(b)所示。本算例中同算例1，檢驗設計域上孔洞的形狀對優化結果的影響，分別設計域上設計了均布的圓形孔洞和方形孔洞如圖10(c)和圖10(e)所示。體積比V設定為50%。圓形孔洞網格模型和方形孔洞網格模型，拓撲優化仿生方法計算40次循環后的優化結果如圖10(d)和圖10(f) 所示。從圖10(d)和圖10(f)可看出孔洞的形狀對最終拓撲優化結果有一定的影響。分別為蟻群方法[18](圖11(a))、人工蜂群優化算法[19](圖11(b))、混合元胞自動機優化方法[16](圖11(c))、元胞自動機優化方法[26](圖11(d)) 、SIMP方法[3](圖11(e))和水平集方法[25](圖11(f))方法所得到的拓撲優化結果圖11所示。從圖10和圖11比較可看出本文方法所得的拓撲優化結果其它六種拓撲優化結果比較相似。圖12為體積比V=50%時圓形孔洞和方形孔洞對應的應變能密度和體積比,隨迭代次數的演化歷程，從圖12可看出首先方孔對應的應變能密度比圓孔對應的應變能密度大一些，而相應的體積比小一些，通過15步迭代后兩種模型對應的應變能密度保持持平，之后迭代步驟應變能密度逐步下降，而體積比經30步迭代后保持持平，37～40步驟之后體積比都能收斂到50%。

圖10 本文所得到的Michell型結構最終拓撲優化結果

圖11 其它方法獲得拓撲優化結果

圖12 Michell型結構的應變能密度與體積比,隨迭代次數的演化歷程

4 拓撲優化仿生方法的應用

本小節中分別對三種不同邊界條件下初始設計區域為3L×L長懸臂梁和四種不同邊界條件下初始設計區域為3L×L的拱橋模型進行拓撲優化計算。

算例3和算例4中結構邊長度取為：L=100 mm，初始設計區域被劃分為300×100共30 000個的pixel單元。體積比V都設為45%。優化迭代步驟設定為50次。材料參數同算例1和2一樣取為：彈性模量70 GPa,泊松比0.3。

算例3：長懸臂梁不同邊界條件下的拓撲優化設計

圖13(a)為結構上下邊中點的兩個端點固定，左右邊中點處,分別作用于F=100 N集中載荷的長懸臂梁初始模型，通過仿生拓撲優化得到的最終拓撲優化結果如圖13(b)所示。

圖13 本文所得到的長懸臂梁的最終拓撲結果

圖14(a)為一端鉸支、一端滾珠支承,結構上邊緣中部鉛垂方向施加F=100 N集中載荷的長懸臂梁初始模型，通過仿生拓撲優化得到的最終拓撲優化結果如圖14(b)所示。

圖14 本文所得到的長懸臂梁的最終拓撲結果

圖15 本文所得到的長懸臂梁的最終拓撲結果

圖15(a)為一端鉸支、一端滾珠支承,結構下邊緣相等距離的三個點上鉛垂方向施加F=100 N集中載荷的長懸臂梁初始模型，通過仿生拓撲優化計算得到的最終拓撲優化結果如圖15(b)所示。

從圖13(b)～圖15(b)可看出通過本文方法所得到的長懸臂梁結構的最終拓撲優化結果具有規則性和對稱性。

算例4：拱橋的拓撲優化設計

本算例中分別給出了初始設計區域為300 mm×100 mm的上、中、下承式拱橋橋梁模型，結構下部兩角部位設固支約束，橋面系設置在拱圈之上，橋面之上作用于F=1 N的均布荷載。

圖16(a)為橋面系位于設計域上面位置的上承式拱橋的初始設計區域，圖16(b)為通過仿生拓撲優化計算所得到的上承式拱橋最終拓撲優化結果。

圖16 本文所得到的上承式拱橋模型的最終拓撲優化結果

圖17(a)為橋面系位于設計域的正中央位置的中承式拱橋的初始設計區域，圖17(b)為通過仿生拓撲優化計算所得到的中承式拱橋最終拓撲優化結果。

圖17 本文所得到的中承式拱橋模型的最終拓撲優化結果

圖18(a)為橋面系位于設計域1/4高度位置的中承式拱橋的初始設計區域，圖18(b)為通過仿生拓撲優化計算所得到的中承式拱橋最終拓撲優化結果。

圖18 本文所得到的下中承式拱橋模型的最終拓撲優化結果

圖19(a)為橋面系位于設計域,最下面位置的下承式拱橋的初始設計區域，圖19(b)為通過仿生拓撲優化計算所得到的下承式拱橋最終拓撲優化結果。

圖19 本文所得到的下承式拱橋模型的最終拓撲優化結果

從圖16～19可看出本文方法所得的上、中、下承式拱橋橋梁模型最終拓撲優化結果具有規則性和對稱性。

5 結 論

本文根據骨骼功能適應性原理，通過種群動力學模型和有限元方法的耦合對連續體結構提出了仿生拓撲優化計算方法，并編寫了基于種群動力學模型的仿生拓撲優化軟件。

通過對兩個連續體結構拓撲優化計算中被廣泛應用的經典數值例子進行拓撲優化計算并將其優化結果與其它幾種拓撲優化方法比較，驗證了本文拓撲優化仿生方法的有效性，最后對三種不同邊界條件下的長懸臂梁模型和四種不同邊界條件下的上、中、下承式拱橋橋梁模型進行拓撲優化計算，結果表明本文提出的仿生拓撲優化方法所得到的最優拓撲較為合理，容易被工程界采納，具有一定的使用價值。

[1] Bends?e 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.

[2] Mlejnek H P, Schirrmacher R. An engineer’s approach to optimal material distribution and shape finding [J]. Computer Methods in Applied Mechanics and Engineering, 1993, 106(1): 1-26.

[3] Bends?e M P, Sigmund O. Topology optimization: theory, methods, and applications[M]. New York: Springer, 2003.

[4] Huang X, Xie Y M. Evolutionary topology optimization of continuum structures: methods and applications[M]. United Kingdom: John Wiley & Sons Ltd , 2010.

[5] Sethian J A, Wiegmann A. Structural boundary design via level set and immersed interface methods [J]. Journal of Computational Physics, 2000, 163(2): 489-528.

[6] Yunkang S, Deqing Y. A new method for structural topological optimization based on the concept of independent continuous variables and smooth model [J]. Acta Mechanica Sinica, 1998, 14(2): 179-185.

[7] Eschenauer H A, Olhoff N. Topology optimization of continuum structures: A review*[J]. Applied Mechanics Reviews, 2001, 54(4): 331-390.

[8] 夏天翔，姚衛星. 連續體結構拓撲優化方法評述[J].航空工程進展,2011,2(1):1-12.

XIA Tian-xiang，YAO Wei-xing . A survey of topology optimization of continuum structure [J].Advances in Aeronautical Science and Engineering, 2011,2(1):1-12.

[9] Mow Van C, Huiskes Rik. Basic orthopaedic biomechanics and mechano-biology [M].3rd Edition, Philadelphia: Lippincott Williams & Wilkins, 2005.

[10] Wolff J．The law of bone remodeling[M]．Berlin：Springer 1986.

[11] Huiskes R, Ruimerman R, Van Lenthe G H, et al. Effects of mechanical forces on maintenance and adaptation of form in trabecular bone[J]. Nature, 2000, 405(6787): 704-706.

[12] Tezuka K, Wada Y, Takahashi A, et al. Computer-simulated bone architecture in a simple bone-remodeling model based on a reaction-diffusion system[J]. Journal of Bone and Mineral Metabolism, 2005, 23(1): 1-7.

[13] Bagge M. A model of bone adaptation as an optimization process [J].Journal of Biomechanics, 2000,33(11):1349-1357.

[14] Zhu X H, He G, Bingzhao G. The application of topology optimization on the quantitative description of the external shape of bone structure [J]. Journal of Biomechanics, 2005, 38(8): 1612-1620.

[15] Mattheck C. Trees: the mechanical design [M]. Springer-Verlag, 1991.

[16] Andr’es T. Bone remodeling as a hybrid cellular automaton optimization process[D].Notre Dame：University of Notre Dame，2004．

[17] Cai K, Chen B S, Zhang H W. Topology optimization of continuum structures based on a new bionics method[J]. International Journal for Computational Methods in Engineering Science and Mechanics, 2007, 8(4): 233-242.

[18] Kaveh A, Hassani B, Shojaee S, et al. Structural topology optimization using ant colony methodology[J]. Engineering Structures, 2008, 30(9): 2559-2565.

[19] Lee C H, Park J Y, Park J Y, et al. Application of artificial bee colony algorithm for structural topology optimization[C]//Natural Computation (ICNC), 2012 Eighth International Conference on. IEEE, 2012: 1023-1025.

[20] 買買提明．艾尼. 用有限元法耦合反應擴散模型的骨重建和仿生拓撲優化方法研究[J].新疆大學學報(自然科學版) , 2009,26(4):402-407.

Mamtimin.Geni .Study on bone remodeling and bionic topology optimization method by using reaction-diffusion model coupled with FEM [J]. Journal of Xinjiang University (Natural science) 2009,26(4):402-407.

[21] 買買提明·艾尼，賈麗華，開依沙爾·熱合曼,等. 受多向載荷骨組織的仿生拓撲優化[C]. 鄭州:中國力學學會學術大會2009 (CCTAM2009),2009:824-826.

[22] Kaysar.Rahman, Mamtimin.Geni, et al. A new bionic topology optimization method based model of bone adaptation [J]. Applied Mechanics and Materials, 2013, 433: 2254-2259.

[23] 陸征一，王穩地.生物數學前沿[M].北京：科學出版社,2008,6:98-112.

[24] 買買提明．艾尼. 基于波形法的有限元分析輔助軟件系統[p].中國：2007-0019, 2007,12.

[25] Yamada T, Izui K, Nishiwaki S, et al. A topology optimization method based on the level set method incorporating a fictitious interface energy [J]. Computer Methods in Applied Mechanics and Engineering, 2010, 199(45): 2876-2891.

[26] Bochenek B, Tajs-Zielińska K. Novel local rules of cellular automata applied to topology and size optimization [J]. Engineering Optimization, 2012, 44(1): 23-35.