考慮位移約束的結構布局優化方法

易繼軍 李艷梅 榮見華 曾 韜

1.中南大學，長沙，410083 2.長沙理工大學，長沙，410004 3.湖南水利水電職業技術學院，長沙，410131

0 引言

結構布局優化設計包括結構拓撲優化、結構形狀優化和結構尺寸優化［1］。多數的優化方法是在概念設計階段進行拓撲優化，在詳細設計階段進行尺寸和形狀優化。此類優化算法由于設計空間的隔離，可能無法獲得最優解［2］。連續體結構尺寸和拓撲組合優化是在優化運算中同時考慮結構尺寸和拓撲變量，存在著變量類型和尺寸度量上的差異，直接對它們進行同時處理有一定困難。現有的組合優化主要針對桁架結構。Steven等［3］應用漸進優化方法（ESO）對離散結構進行了拓撲和尺寸的組合優化。劉濤等［4］提出了一種桁架拓撲、形狀和尺寸組合優化的漸進優化方法，該方法采用拓撲優化和形狀尺寸優化兩個子問題分層處理。Rajan等［5］利用遺傳算法（GA）分別進行了6節點桁架和14節點桁架的優化。Chan等［6］采用準則法和遺傳算法的混合算法對鋼框架結構進行了拓撲和尺寸優化。在桁架形狀與尺寸組合優化方面，常用的方法可以分為數學規劃法［7］和漸進優化方法［4］兩類。在連續體組合優化方面，Afonso等［8］提出了基于準則優化法的盤類結構拓撲和尺寸的分級優化策略。王偉等［9］為解決機翼結構布局優化設計問題，提出了一種兩級三層優化設計方法。榮見華等［10］基于ESO方法，提出了板梁組合的連續體結構優化方法。Zhou等［2］提出了一種拓撲、尺寸和形狀集成的優化算法，由于沒有考慮變量尺寸度量上的差異，提出的算法存在不足。

本文提出一種連續體結構組合優化方法，同時考慮了設計變量類型和尺度上的差異，實現了拓撲和尺寸變量的協同優化。

1 組合結構拓撲變量及優化數學模型

在初始設計域中，整個結構被分成兩大類區域。在一類子區域中，單元的尺寸（如厚度）作為設計變量，即尺寸設計變量；在另外一類子區域中，進行拓撲設計，拓撲設計單元的拓撲變量作為設計變量。在尺寸設計區域，不能進行拓撲設計的單元用p（p＝1，2，…，P）編號，尺寸設計變量用z（z＝1，2，…，H）編號，Pz和Sz分別表示第z個子區域的單元數量和面積，尺寸設計變量用hz表示。在拓撲設計區域，Q表示拓撲設計變量的個數，拓撲設計單元用q（q＝1，2，…，Q）編號，Sq表示其單元面積，它們的拓撲變量用ρq（q＝1，2，…，Q）表示。

1.1 拓撲變量和尺寸變量

類似于SIMP（solid isotropic microstructure with penalization）方法，采用材料屬性的有理近似模型（rational approximation of material properties，RAMP）的過濾方法，將離散的拓撲變量轉換成［0，1］區間的連續變量。通過懲罰因子對設計變量在（0，1）之間的中間密度值進行懲罰，使連續變量的拓撲優化模型能很好地逼近傳統的0/1離散變量的拓撲優化模型。使用RAMP過濾函數fk（ρiq）識別單元剛度，用過濾函數fv（ρiq）識別單元體積，單元特性參數識別采用如下公式：

對于兩類變量在尺寸度量上的差異，采用歸一化處理方法將尺寸變量轉化成［0，1］區間的連續變量優化問題。即，用hz/代替厚度變量hz，這樣厚度變量的約束區間變為［hz/，1］。

1.2 優化數學模型

位移約束的最小體積的布局優化問題可以表示為

2 近似優化模型和求解方法

2.1 變位移約束限的序列近似優化模型

為了使式（2）中的位移近似函數在每次迭代循環中保持近似成立，建立了變位移約束限的近似系列模型。這些變化的位移約束限起到了設計變量的置信區間的作用。通過這個方法，能夠保持設計變量的小量變化，從而確保迭代過程求解的有效性［9-10］。模型如下：

2.2 近似優化模型的求解方法

定義

這里ju（i）表示在第j組載荷作用下，第i號單元的位移矢量。

對于式（5），省略常數項，則式（3）轉換成：

這里

將目標函數用二階泰勒展開近似代替，那么式（6）的求解能夠轉換成對下式的求解：

其中

根據對偶理論，類似于文獻［11-12］，式（7）轉換成下面的對偶規劃問題：

這里

這里采用交互式的循環迭代策略來近似求解式（7）和式（8）的最優解（i＝1，2，…，Q＋H）和拉格朗日乘子λ。為了求解拉格朗日乘子，將目標函數φ（λ）≈L（x＊，λ）進行近似二階泰勒展開。對L（x＊，λ）求關于λ的一階和二階導數：

對式（8）應用K－T條件，有

并且，由式（11）的第一式和第三式，有

根據式（9）～式（15），容易得到下面的方程：

如果φ（λ）用其二階近似表達，同時省略常數項，我們能夠得到下面的二次規劃方程：

這里，H為L×J階向量；D為（L×J）×（L×J）矩陣。

通過二次規劃法求解式（18），得到拉格朗日乘子。然后通過式（12）得到變量（i＝1，2，…，Q＋H）的新的近似值。用新的（i＝1，2，…，Q＋H）更新系數矩陣D和H。再次求解式（18），得到新的（i＝1，2，…，Q＋H）。這個循環過程一直反復進行，直到滿足收斂條件 ‖x（s＋1）－x（s）‖/‖x（s）‖ ≤ε3（s是求解二次規劃的循環迭代次數）。

3 優化策略和設計空間調整

3.1 優化策略

在本文研究的方法中，整個優化過程被分成兩個優化調整階段和一個相變轉換階段。在第一個階段獲得包括灰色區域和尺寸改變的拓撲結構，當設計結構到達優化目標的鄰域，如果滿足下式優化進入相變轉換階段：

在此階段，進行設計空間調整，移除部分灰色區域，然后優化進程進入第二個優化調整階段，在此階段，進一步移除灰色區域，在滿足位移約束的前提下，如果滿足收斂準則，則迭代終止。收斂準則為

這里k表示當前循環迭代步，τVmax是最大允許誤差。N是一個整數，在本文中N＝3。這意味著在6個收斂迭代步中達到穩定［13］。

優化流程如圖1所示。

圖1 優化流程圖

3.2 設計空間調整

為了移除灰色區域，同時為了易于求解具有較大規模有限元網格模型的優化問題，這里采用在一些優化迭代步將拓撲變量很小的單元從結構上去掉，而其他單元拓撲變量不變的策略。即設置，使得＋ε（i＝1，2，…，Q），ε為一小量。把每個可設計的保留單元的拓撲變量與閾值相比較來確定單元是否存在的狀態，基于式（22），挑選用于單元刪除的候選單元，即將這些單元看成將要消失的單元。

為了確保優化迭代中的結構非奇異，同時結構優化時具有設計空間的擴展（增添單元）功能，這里采用在結構邊界和孔洞周圍附加人工材料單元的辦法［11，14］來確保結構剛度矩陣非奇異。將人工材料單元的彈性模量設為占據該空間的原始材料單元的彈性模量，僅將人工材料單元的拓撲變量取為ρq。人工材料單元參與結構參數計算和優化迭代。根據人工材料單元拓撲變量的變化情況，進行設計空間的擴展。設增添單元的閾值為，基于式（23），挑選出用于單元增添的候選單元，即將這些單元看成將要增添的單元。

設計空間縮減和擴展的準則表為

式中，nd、na分別為集合Nd和Na的單元個數，上述和都是小的經驗參數值。

如果式（24）成立，則執行以下操作：

當執行式（25）和式（26）的操作后，非零的所對應的單元材料特性參數編號變為保留單元的特性參數編號，為零的所對應的單元材料特性參數編號變為無材料單元所對應的特性參數編號。這一輪的拓撲優化就完成了設計空間的減縮和擴展，否則，該輪拓撲優化不進行設計空間的減縮和擴展。對于相變轉換階段的設計空間調整參照式（25）和式（26）進行。

當拓撲設計區域發生改變時，如果不調整厚度，將會使兩類變量變化不一致，因此有必要擬定厚度調整方案。根據迭代過程中拓撲變化區域和厚度變化區域體積小量變化以及相關變量在前后迭代步的連續性原則，可以建立以下近似式：

4 算例

4.1 算例1—— 立柱

初始的結構模型和最大的設計域如圖2所示。整個結構被分成了8個子區域，其中子區域6～8為拓撲設計區域，其單元采用四節點平面應力單元，子區域1～5為尺寸變量設計區域，四節點平面應力單元厚度作為設計變量，且在一個子區域中單元厚度變量相同。兩組載荷被施加到底部固定的結構上。一組載荷P1＝500N，方向沿水平向右，施加到結構左上角；另一組載荷P2＝500N，方向水平向左，施加到結構右上角。最大設計區域1.0m×3.0m×0.002m （Lx＝1.0m，Ly＝3.0m）離散成50×150共7500個四節點平面單元，它們的初始厚度是0.002m，允許的最小厚度是0.0002m。位移約束點在結構頂部的中間位置，水平方向的初始位移是0.153mm，允許值是0.400mm。彈性模量E＝210GPa，泊松比ν＝0.3。在第一階段單元刪除閾值和增加閾值都為0.001，在第二階段單元刪除閾值和增加閾值都為0.01，在階段轉換迭代步刪除閾值為0.1。

圖2 模型初始結構和最大的設計域（算例1）

結構的拓撲優化歷程如圖3所示。在第70步滿足收斂條件式（20），優化迭代進入相變轉換階段，移除部分拓撲優化的灰色區域并調整厚度，圖3b為執行相變轉換后的結構拓撲。在第90步（圖3c）滿足收斂條件式（21），優化迭代結束，其約束點的位移是0.400mm。

圖3 結構的拓撲優化歷程（算例1）

結構的體積和約束點優化歷程如圖4、圖5所示。

圖4 體積優化歷程（算例1）

圖5 位移優化歷程（算例1）

當第一個優化調整階段結束的時候，得到的拓撲構型接近最優結構，通過階段轉換步的操作，結構設計區域中的黑/白分布清晰，最優結構布局在第二階段獲得。最優解的結構體積是0.002 538 7m3，結構位移是0.400mm，結構的厚度分布如表1所示。結構的厚度優化歷程見圖6。

表1 優化后的厚度值 mm

圖6 厚度迭代歷程（算例1）

4.2 算例2——簡支梁

初始的結構模型和最大的設計域如圖7所示。整個結構被分成了9個子區域，其中子區域7～9為拓撲設計區域，其單元采用四節點平面應力單元，子區域1～6為尺寸變量設計區域，四節點平面應力單元厚度作為設計變量，且在一個子區域中單元厚度變量相同。一組載荷P＝－2000N被施加到結構中部（如圖7a所示），方向向下。結構底部左端固定，右端簡支。最大設計區域1.0m×3.0m×0.002m （Lx＝1.0m，Ly＝3.0m）離散成50×150×1共7500個四節點平面單元，它們的初始厚度是2mm，允許的最小厚度是0.2mm。位移約束點在結構頂部的中間，垂直向下初始位移是－0.0647mm，允許值是0.300mm。在第一階段單元刪除閾值和增加閾值都為0.001，在第二階段單元刪除閾值和增加閾值都為0.01，在階段轉換迭代步刪除閾值為0.05。

圖7 模型初始結構和最大設計域（算例2）

結構的拓撲優化歷程如圖8所示。在第135步滿足收斂條件（式（20）），優化迭代進入相變轉換階段，移除部分拓撲優化的灰色區域并調整厚度，圖8b為執行相變轉換后的結構拓撲。在第296步（圖8c所示）滿足收斂條件（式（21）），優化迭代結束，其約束點的位移是0.300mm。

圖8 結構的拓撲優化歷程（算例2）

結構的體積優化和約束點優化歷程如圖9、圖10所示。

圖9 體積優化歷程（算例2）

圖10 位移優化歷程（算例2）

當第一個優化調整階段結束的時候，得到的拓撲構型接近最優結構，通過階段轉換步的操作，結構設計區域中的黑/白分布清晰，最優結構布局在第二階段獲得。最優解的結構體積是0.002 538 7m3，結構位移是－0.300mm，結構的厚度分布如表2所示。結構的厚度優化歷程見圖11。

表2 優化后的厚度值 mm

圖11 厚度迭代歷程（算例2）

5 結論

本文研究了拓撲變量和尺寸變量組合的結構優化設計方法。采用材料屬性的有理近似模型的過濾方法，將離散的拓撲變量轉換成［0，1］區間的連續變量；通過歸一化的處理方法使兩類變量變成同一尺度。將組合優化問題轉化為單一類型變量優化問題。采用變位移約束限和二次規劃法求解優化問題，通過設計空間調整的方法移除灰色區域，同時調整厚度，得到了清晰的拓撲構型和優化的厚度值。給出的算例說明了本文方法的正確性和有效性。

［1］Rozvany G I N，Zhou M.Advances in Overcoming Computational Pitfalls in Topology Optimization［C］//6th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization.Reston，1996：1122－1132.

［2］Zhou M，Pagaldipti N，Thomas H L，et al.An Integrated Approach to Topology，Sizing and Shape Optimization［J］.Structural and Multidisciplinary Optimization，2004，26（5）：308－317.

［3］Steven G，Querin O，Xie M.Evolutionary Structural Optimisation （ESO）for Combined Topology and Size Optimisation of Discrete Structures［J］.Computer Methods in Applied Mechanics and Engineering，2000，188（4）：743－754.

［4］劉濤，鄧子辰.桁架結構尺寸和形狀、拓撲的漸進優化方法［J］.西北工業大學學報，2004，22（6）：739－743.

［5］Rajan S D.Sizing，Shape，and Topology Design Optimization of Trusses Using Genetic Algorithm［J］.Journal of Structural Engineering，1995，121（10）：1480－1487.

［6］Chan C－M，Wong K－M.Structural Topology and Element Sizing Design Optimisation of Tall Steel Frameworks Using a Hybrid OC– GA Method［J］.Structural and Multidisciplinary Optimization，2008，35（5）：473－488.

［7］劉軍偉，姜節勝.桁架動力學形狀優化的統一設計變量方法［J］.振動工程學報，2000，13（1）：84－88.

［8］Afonso S B，Belblidia F，Hinton E，et al.A Combined Structural Topology and Sizing Optimization Procedure for Optimum Plates Design［C］//European Congress on Computational Methods in Applied Sciences and Engineering.Barcelona，2000：125－131.

［9］王偉，楊偉，趙美英.大展弦比飛翼結構拓撲、形狀與尺寸綜合優化設計［J］.機械強度，2008，30（4）：596－600.

［10］Rong J H，Jiang J S，Xie Y M.Evolutionary Structural Topology Optimization for Continuum Structures with Structural Size and Topology Variables［J］.Advances in Structural Engineering，2007，10（6）：681－695.

［11］榮見華，邢曉娟，鄧果.一種變位移約束限的結構拓撲優化方法［J］.力學學報，2009，41（3）：431－439.

［12］榮見華，張強，葛森，等.基于設計空間調整的結構拓撲優化方法［J］.力學學報，2010，42（2）：256－267.

［13］Huang X，Xie Y M.Optimal Design of Periodic Structures Using Evolutionary Topology Optimization［J］.Structural and Multidisciplinary Optimization，2008，36（6）：597－606.

［14］Stolpe M，Svanberg K.An Alternative Interpolation Scheme for Minimum Compliance Topology Optimization［J］.Structural and Multidisciplinary Optimization，2001，22（2）：116－124.