基于變密度法的結構強度拓撲優化策略

丁 卯，耿 達，周明東, 來新民

(上海交通大學 機械與動力工程學院， 上海 200240)

復雜工況下的結構設計必須考慮結構在外力作用下應力的分布情況，以避免局部應力過高導致結構失效等問題[1].現代工業的不斷發展對具有高比剛度及高比強度結構的需求愈發迫切.結構拓撲優化方法基于最優化理論，結合計算機輔助工程(CAE)技術，自動地設計符合性能要求的工程結構，近年來已作為新興的結構設計方法在航空航天及汽車等領域的高比剛度、高比強度結構設計中得到了廣泛的研究與應用[2-3].目前結構比剛度問題的拓撲優化已得到充分發展，而考慮結構強度的拓撲優化方法仍存在許多問題[4-6].

基于變密度法的結構拓撲優化方法優化過程中結構拓撲變化靈敏，被廣泛應用于實際工程結構的優化設計.但傳統基于變密度法的結構強度拓撲優化方法受優化問題非線性影響，結構近似最大單元應力難以準確計算，應力變化趨勢難以精準評估[7-8],且優化過程對優化參數的選擇較敏感，不易收斂.同時最終優化結果往往存在較多無實際物理意義的灰度單元，后處理前后結構強度性能偏差較大，常需反復試錯才能得到符合強度設計需求的結構.Zhou等[9]采用基于單純復形的結構拓撲優化方法，使用二階單元劃分結構網格，并選取較高應力指數系數以構建結構最大單元應力近似計算公式[10]，使優化結果的應力水平得以降低.Lian等[11]采用“結構拓撲開洞-結構形狀優化-可變復型網格跟新”的迭代循環方法進行考慮結構應力的拓撲優化，優化結果不存在灰度單元，且優化過程中結構應力水平近似計算較準確.

在此基礎上，本文針對基于變密度法考慮結構強度拓撲優化問題，采用過濾-投影結構參數化方法，通過研究主要優化參數對結構比強度優化問題優化過程的影響規律，提出一種新的優化策略，實現了在拓優化過程中對結構應力水平及其變化趨勢的準確控制，同時保證優化過程的穩定性.最后通過典型結構的算例實驗對該方法的合理性及實用性進行驗證.

1 考慮結構強度的變密度法拓撲優化 模型

基于變密度法的結構強度拓撲優化研究一直以來是工程應用關注的對象和學界的研究熱點.本文基于固體各向同性材料懲罰(SIMP)法，對以結構體積為約束，最小化結構最大應力的拓撲優化問題進行研究[12].其優化列式可表述為

s.t.KU=F

(1)

0≤ρe≤1,e=1,2,…,ne

式中：ρ={ρ1,ρ2,…,ρe}表征各單元的相對密度，其中ρe在0與1之間連續變化，ρe=0和ρe=1分別表示該位置為空單元及實體單元，0<ρe<1則表示此處單元具有中間密度材料；σmax為結構單元最大應力，表征結構整體應力水平，本文采用Von Mises屈服準則對單元應力進行表達(后文亦稱σmax為最大單元Von Mises應力)；K、U及F分別為結構整體剛度矩陣、位移向量及外力向量；ve為對應單元的體積；V*為所允許的結構體積上限；ne為結構內單元總數量設計變量場.

(2)

式中：Ωe為單元e附近所有形心落在過濾區域內的單元集合；加權函數ωj為相應單元中心離距的線性函數；ρj為該對應單元的相對密度；單元e的過濾區域為以該單元形心為圓心，過濾半徑為r0的圓形區域；rj為Ωe內單元j與單元e的中心離距.

本文采用Von Mises屈服準則進行單元應力計算[13-14]，并構建如下所示的類SIMP法單元應力插值函數，以解決單元應力奇異性問題：

(3)

(4)

e=1,2,…,ne

式中：σ11、σ22及σ12為單元e應力矢量的3個分量.

單元應力是局部物理量，優化過程中需對各單元應力進行控制以實現對結構整體應力水平的預測.而將每個單元應力均作為優化約束時，存在約束過多、難以求解的問題.因此，本文采用如下所示的p-norm法近似計算：

(5)

2 結合投影法的參數化方法

基于傳統變密度拓撲優化方法所設計的結構中，中間密度單元使結構應力在對應位置處無法準確計算，而為獲得單元密度非0即1的離散結構，一般對優化結果采用如下所示的投影后處理，以去除結構中間密度單元：

(6)

(7)

圖1 投影函數曲線Fig.1 Curves of projection function

針對上述具有較高非線性的結構應力優化問題，本文采用移動漸近線(MMA)法[16]作為優化器對其進行尋優.優化器中設計變量變化閾值(ML)表征每次優化迭代中各設計變量的變化范圍，亦即優化算法中的搜索步長[12]， 其同p-norm法的p和投影函數的β耦合，綜合影響優化過程的穩定性和收斂性.本文基于式(1)的優化模型，通過控制變量法對上述3個主要優化參數的影響規律進行研究.

以圖2所示的L型梁作為優化對象，圖中R為半徑.采用 6 400 個邊長為1 mm的一階四節點四邊形單元對其進行網格劃分，L型梁左上端固定，為避免產生應力集中，將豎直向下大小為5 N的力分解加載在結構右端扇形區域內的所有節點上，實體材料彈性模量E=1.0 MPa，泊松比ν=0.3，優化過程中結構體積約束值為0.3.

圖2 二維L型梁有限元模型(mm)Fig.2 Finite element model of a L-shaped beam (mm)

圖3～6所示為ML、p及β取不同值時的結構拓撲、應力分布云圖及結構最大應力.結構材料分布圖下方的應力云圖中，顏色從藍到紅表征單元應力由0逐漸增大至結構最大單元應力值，如圖7所示.圖3～6中，以相同p和ML作為一組算例，通過對比各組算例的優化結果可知，β越大，優化結果的中間密度單元比例越低，越能避免如圖3(q)、3(r)、3(t) 等所示的結構斷裂現象.而由各圖中的(a)、(e)、(i)、(m)及(q)可知，隨著p不斷增大，優化結果的最大單元Von Mises應力有逐漸下降的趨勢.當p取12及30等較高值時，優化過程受優化問題非線性影響較大，結構最大單元應力值有上升趨勢.由此可知，在實際應用中，當p和β取較小值時，優化問題非線性較弱，結構內中間密度單元比例較大，此時設計變量變化閾值ML宜取較大值，以使結構拓撲變化更靈敏，從而更快地獲得具有較低應力水平的優化結果.當為保證優化過程中對結構應力水平和最大應力變化趨勢的控制精度而采用較高p值，同時為降低優化結果中間密度單元比例采用而較高β值時，應使用較低的設計變量變化閾值ML.

圖3 β=2時優化結果Fig.3 Optimized results at β=2

圖4 β=4時優化結果Fig.4 Optimized results at β=4

圖5 β=8時優化結果Fig.5 Optimized results at β=8

圖6 β=16時優化結果Fig.6 Optimized results at β=16

圖7 優化結果應力分布圖Fig.7 Stress distribution map of optimized results

3 優化策略

基于上一節所研究的主要優化參數對優化過程及優化結果的影響規律，本文提出了先結構拓撲優化，后近似形狀優化的優化策略，以同時提升優化過程中結構應力水平的控制精度和優化過程的穩定性，并使最終設計結果的強度性能滿足要求.首先在優化過程的拓撲優化階段，通過采用較小β及p作為初始值，降低優化問題非線性影響，以實現結構拓撲的快速變化.待符合收斂條件后，增大p至30以提高優化過程中結構應力水平控制精度.滿足收斂條件后增大β至64以大幅降低結構中間密度單元比例，使結構變化主要發生在邊界，以構建近似形狀優化的結構優化過程，最終采用后處理投影方法，獲得非黑即白的設計結果.

為探明拓撲優化階段p及整個優化過程中ML的選取方法，同時驗證本文優化策略的有效性，以式(1)為優化模型，并以圖2所示的L型梁為研究對象進行優化實驗.為提高單元應力計算精度，降低非線性對優化問題的影響，采用二階矩形單元對L型梁結構進行網格劃分，網格數量、尺寸、荷載、邊界條件及體積約束等設置與前文相同.在不同的初始p值(圖中p′)及ML值下，最終設計結果的σmax變化曲線如圖8～10所示.

圖8 ML=0.05時優化結果的σmaxFig.8 σmax of optimized results at ML=0.05

圖9 ML=0.1時優化結果的σmaxFig.9 σmaxof optimized results at ML=0.1

圖10 ML=0.2時優化結果的σmaxFig.10 σmax of optimized results at ML=0.2

根據圖8～10可知，當p′≤9時，優化結果的最大Von Mises應力較小.當p′>9時，優化問題的非線性程度較高，優化結果的最大Von Mises值較采用低p′值時更高，整體呈增大趨勢，且存在較大振蕩現象.當ML取0.1時，易于獲得具有更低應力水平的最終設計結果，并保證結構拓撲的靈敏變化.優化策略如圖11所示.

圖11 優化策略流程Fig.11 Flowchart of optimization strategy

為了驗證上述方法的合理性與一般性，以圖12所示的Portal結構[11]作為驗證實驗的優化對象.采用邊長為1 mm的二階正方形單元對該結構進行網格劃分，結構左下角5個節點的所有自由度固定，右下角5個節點僅固定豎直方向的自由度，大小為5 N的力平均作用在結構頂部中間5 mm范圍的節點上.采取與L型梁算例相同優化參數設置構建對比實驗，設計結果的最大單元Von Mises應力變化曲線如圖13～15所示.

圖12 Portal結構尺寸及邊界條件(mm)Fig.12 Size and boundary condition of the portal structure (mm)

圖13 ML=0.05時優化結果的σmaxFig.13 σmax of optimized results at ML=0.05

根據圖中驗證實驗的優化結果可知，設計結果最大單元應力的變化趨勢以及優化參數對設計結果強度的影響規律與L型梁實驗一致，該系列優化結果以圖16為例，表明上文所總結的優化參數影響規律及提出的參數選取策略具有一定的普適性及實用性，可用于指導工程實際應用.

圖14 ML=0.1時優化結果的σmaxFig.14 σmax of optimized results at ML=0.1

圖15 ML=0.2時優化結果的σmaxFig.15 σmaxof optimized results at ML=0.2

圖16 Portal結構優化結果Fig.16 Optimized results of portal structure

5 結語

本文采用基于過濾-投影的參數化方法，實現在迭代優化過程中結構中間密度單元比例的不斷下降，進而降低后處理前后優化結果應力水平的偏差，使其滿足工程實際的設計需求.通過對經典L型梁結構的比強度優化實驗結果的分析，研究了關鍵優化參數對優化過程及優化結果應力水平的影響規律，并在此基礎上，提出了先結構拓撲優化，后近似形狀優化的新優化策略，同時在優化過程中保證了對結構應力水平的控制精度以及優化結果的穩定性.在上述優化策略下對典型Portal結構進行優化設計，并對優化結果的應力水平進行了分析，闡明了該方法的合理性及實用性，可為實際工程中結構比強度拓撲優化設計提供指導.