考慮自重影響的拓撲優化方法

宋英杰 吳超 尹書翰 黃俊斌 宋勇宏

摘 要：針對在考慮自重影響的連續體結構拓撲優化中常見的問題，本研究提出一種新的拓撲優化方法。首先，為避免寄生效應，提出一種帶懲罰的固體各向同性材料（SIMP）改進模型，建立剛度模量、質量密度和體積光滑懲罰函數三者間的合理匹配關系。同時，引入密度過濾、Heaviside映射的三場過濾方案，并提出一種低物理密度單元的體積變化率約束措施，從而極大地提高考慮自重效應的結構拓撲問題的優化計算效率，獲得清晰的0-1分布。基于上述方法，建立以結構柔順度最小為目標，同時考慮變體限約束、位移約束以及低物理密度單元的體積變化率約束等措施的優化模型。最后，采用MMA算法對優化模型進行求解，并通過算例驗證本研究所提出的方法是正確有效的。

關鍵詞：拓撲優化;連續體結構;結構自重;低密度約束

中圖分類號：TH122 ? ? 文獻標志碼：A ? ? 文章編號：1003-5168（2022）12-0030-06

DOI：10.19968/j.cnki.hnkj.1003-5168.2022.12.006

A Topology Optimization Method Considering the Effect of Self-Weight

SONG Yingjie? ? WU Chao? ? YIN Shuhan? ? HUANG Junbin? ? SONG Yonghong

（School of Automotive and Mechanical Engineering，Changsha University of Science and Technology，Changsha 410014，China）

Abstract：Aiming at the common problems in topology optimization of continuum structures considering the influence of self-weight，this paper proposes a new topology optimization design method.First，to avoid parasitic effects， an improved Solid Isotropic Material with Penalty （SIMP） model is proposed，which establishes a reasonable matching relationship between stiffness modulus and bulk smoothness penalty function.At the same time，a three-field filtering scheme of density filtering and Heaviside mapping is introduced，and a volume change rate constraint measure for low physical density cells is proposed to greatly improve the optimization calculation efficiency of structural topology problems considering self-weight effects，and obtain clear 0-1 distributionBased on the above method，an optimization model is established，which takes the minimum structural compliance as the goal，and considers the constraints of the variation limit and the volume change rate of the low physical density element.Finally，the MMA algorithm is used to solve the optimization model.And the method proposed in this paper is verified to be correct and effective by an example.

Keywords：topology optimization;continuum structure;self-weight;low density constraint

0 引言

當今新一輪工業革命方興未艾，而其核心驅動力便是智能制造[1]。智能制造技術的發展不僅要使設計出的結構具有高性能和輕量化的特點，還要考慮實際工程應用中結構自重的影響。雖然拓撲優化法在工程機械、航空航天、車輛船舶及建筑土木等領域獲得成功應用，但在早期的拓撲優化過程中，往往忽略了包括自重載荷、熱載荷、旋轉部件的離心載荷和慣性加速度載荷在內的體積力。

Rozvany等率先開展了基于自重的拓撲優化問題研究。Bruyneel等發現，對體積力在內的載荷結構進行拓撲優化時存在3個問題：①在優化過程中，柔順度目標函數不是單調函數;②優化結果的體積不受約束;③低密度區域有寄生效應[2]。近年來，雖然有學者針對上述問題進行了研究，但卻始終無法避免優化過程中的低物理密度單元對結構的影響，且整個優化過程收斂較為緩慢，直接使用固體各向同性材料懲罰模型（SIMP）并不能解決這些問題。

基于此，本研究結合各種拓撲優化方法及相關領域的最新研究進展，提出一種改進的SIMP插值方法，可有效解決考慮結構自重時所帶來的問題。本研究還提出了低密度單元體積變化率的約束方法，在進一步避免寄生效應的同時，既降低了優化過程中低密度單元對結構的影響，獲得清晰的0-1分布，又極大地提高了計算效率。

1 過濾方法和材料懲罰模型

1.1 過濾方法

為建立連續體結構的有限元模型和拓撲優化模型，本研究先將優化結構設計域離散為[N]個四節點的四邊形單元，并為設計域中的每個單元分配一個設計變量[xi]。通過密度過濾技術[3]和Heaviside映射函數[4]來形成由設計變量[xi]、中間變量[xi]、物理變量[xi]組成的三場過濾方案，既可避免拓撲優化過程中出現棋盤格現象，又能解決灰度單元問題。Heaviside映射函數見式（1）。

式中：[β]是Heaviside映射函數的曲率參數，本研究所有算例中[β]的初始值均為10-5，并按照[β（m）=max（1，2 β（m-1））]每50迭代步更新一次（即每50個外循環迭代步m增加1）;本研究中[η]取值為0.04。

1.2 改進的材料懲罰模型

SIMP法是現有拓撲優化方法中最常用的一種方法。針對考慮自重影響的拓撲優化問題，采用經典的SIMP模型來解決時，在低密度區域，結構剛度與其自重存在較大差異，使結構剛度難以承受自重，最終導致結構剛度不足和無界位移。本研究提出一種改進的固體各向同性懲罰材料模型，見式（2），用來解決拓撲優化中自重載荷引起的低密度區域的寄生效應，同時保證體積約束的有效性。

式中：[E0]是初始單元材料的彈性模量;[ρ0]為單元初始質量密度;[p]和[q]分別為單元的剛度模量懲罰因子和體積懲罰因子;[α]為剛度模量的多項式函數的懲罰參數。通過數值仿真研究發現，當[p=5]、[q=2]、[α=25]時，對寄生效應的抑制效果較好，且能獲得較好的拓撲優化結果。

2 優化模型建立及靈敏度分析

2.1 結構自重處理和目標函數

當考慮與質量密度相關的自重時，載荷向量表達式見式（3）。

[F=g+P]? ? ? （3）

式中：[F]為載荷向量;[g]為結構自重向量;[P]為外載荷向量。

考慮平面應力問題的四節點四邊形有限元和沿垂直方向施加的重力載荷，使得每個有限元四分之一的自重沿垂直方向平均分配給其四個節點，見式（4）。

式中：[Nj]為第[j]個節點相連單元編號的集合;[ag]為重力加速度;[v0i]為第[i]號單元的初始體積。

由結構剛度矩陣[K]和位移矢量[U]可以得到結構的柔順度C，見式（5）。

2.2 低密度單元體積變化率約束措施

為抑制優化過程中低密度單元的密度變量波動，并通過提升優化過程中結構清晰度來提高計算效率，本研究提出一個低密度單元體積變化率約束，見式（6）。

[h1（x）=（LV-LV（k-1））/LV（k-1）+γ1≤0? ?if? ?LV≥VLTc]

（6）

式中：[LV]為物理密度變量在0.001～0.100的所有單元體積之和;[VLTc]為是否引入低密度單元體積變化率約束的閾值;[LV（k-1）]為第k-1迭代步[LV]值。式（7）（8）（9）給出了[LV]和[VLTc]的表達式，[γ1]為一個經驗參數，本研究取值為0.05。當[LV≥VLTc]時，在后續優化模型中引入低密度單元體積變化率約束。

式中：[N]為設計域中單元的數量;[Vfull]為結構拓撲中充滿材料時設計域的總體積;[γ2]為一個經驗參數，本研究γ2取值為5。

為了使式（8）更易求導，引入Sigmoid函數[1/（1+e-ζx）]近似單位階躍函數。[LVi]和[LV]的近似平滑函數見式（10）和式（11）。

2.3 位移約束

在實際工程應用中，要在保證結構剛度的前提下，還要保證結構中指定點的位移足夠小。因此，本研究依據需要限制位移的自由度數量，引入一個或多個結構位移約束，見式（12）。

式中：[UUj]為優化結構中[j]點的位移約束限;[nj]為[j]點位移約束所指定的優化結構自由度編號;[Unj]為優化結構第[nj]號自由度的位移;[U（0）nj]為[Unj]在初始迭代步中的值;[J]為需要受到位移約束的自由度的個數，即位移約束的個數。

2.4 優化模型

優化結果體積不受約束是考慮設計變量相關載荷的結構拓撲優化中存在的三個問題之一[2]。為了進一步保持體積約束的有效性，本研究引入了體積上下限約束[V1（x）]和[V2（x）][5]。由于Heaviside曲率參數發生變化時會導致結構位移和柔順度發生較大變化，故采用Rong等[5]提出的變約束限方案，來保證結構性能在拓撲優化過程中的穩定變化。綜上所述，以結構柔順度最小為目標函數，考慮結構自重影響，包含體積相關約束和位移約束的結構拓撲優化模型見式（13）。

式中：[xmin]為一個小量;[V（0）]為初始迭代時優化結構的總體積;[VSd]為優化目標體積;[θ]為體積下限松弛參數，在本研究中θ的取值為0.004;[V（k）Sd]和[UU，（k）j]分別表示在第k個迭代步中結構體積的變化約束限和優化結構的j點位移變化約束限，計算公式見式（14）至式（17）。

式中：[ζ1]和[ζ2]為經驗參數，在本研究的算例中分別設定為0.01和0.04[6];[V（0）]為Heaviside曲率參數[β]變化步以及最接近k步的結構體積;[U（0）nj]為第[nj]號自由度處的位移。

2.5 靈敏度分析

本研究采用MMA算法[7]來求解，因此，需要求解優化模型式（13）中各個性能函數的靈敏度。由鏈式求導法則可知，結構的任意性能函數對設計變量的導數見式（18）。

由懲罰模型式（2）可得整體剛度矩陣關于物理密度變量[xi]的導數，見式（19）。

式中：[Ki]和[K0i]分別為整體坐標系下第i個單元的剛度矩陣和第i個單元的初始剛度矩陣。

與結構自重相關的等效節點力向量對物理密度變量[xi]的導數見式（20）和式（21）。

假設[Vj]為單位虛荷載作用在優化結構的第[nj]號自由度上所產生的結構位移矢量，則結構第[nj]號自由度的位移絕對值對物理密度變量[xi]的導數見式（22）。

基于式（5），優化結構柔順度對物理密度變量[xi]的導數見式（23）。

基于優化模型（13）中的兩個體積不等式約束，可求得其對物理密度變量[xi]的導數，見式（24）和式（25）。

基于式（11），可以求得其對物理密度變量[xi]的導數，見式（26）。

式（22）至式（26）是根據鏈式求導法則得到的位移、柔順度和體積相關函數對設計變量的導數。

3 算例分析

為驗證本研究所提出的拓撲優化方法的準確性，針對圖1所示的考慮自重載荷的橋梁結構進行優化求解。

圖1所示的橋梁結構的荷載和邊界條件是左右對稱的。橋梁結構設計域為一個頂部加載、底部兩端固定支撐的矩形區域，長為20 m、高為6.1 m、厚為0.006 7 m，即橋面為長20 m、高0.1 m、厚度0.006 7 m的非設計域。橋面上豎直向下施加大小為-7 500 N/m均布力。初始材料的彈性模量為[E=2×1011Pa]，其泊松比ν為0.3，其質量密度ρ0為7 850 kg/m3。整個設計域和非設計域全部離散為600×183=109 800個四節點矩形平面應力單元。設置密度過濾半徑為[rmin=2Δ]，其中[Δ]為最大單元邊長。在此示例中，目標體積與設計域初始體積的比值設定為0.5。本研究的算例在考慮自重的前提下，以柔順度最小為目標，包含體積相關約束和位移約束等多種約束條件。其中，頂部中點垂直位移規定為小于0.003 5 m。

圖2給出了使用本研究提出的方法得到的橋梁結構拓撲變量分布的優化歷程圖。其最優拓撲構型如圖2（f）和圖3所示，最優拓撲結構柔順度為113.84 N·m，橋頂中點垂向位移為0.002 78 m。圖4給出了在沒有考慮低密度單元的體積變化率約束時得到的拓撲變量分布的優化歷程圖，其最優拓撲結構如圖4（f）和圖5所示。最優拓撲結構的柔順度和橋頂中點垂向位移分別為114.07 N·m和0.002 79 m。與圖2一樣，圖4中最黑的區域為物理密度變量大于或等于0.1的單元所占區域，從而突出顯示物理密度變量在[0，0.1]范圍內的元素分布，可以有效地對比出低密度單元體積變化率約束的作用。這些數據以及最優拓撲構型圖可以充分表明本研究所提出的方法能夠獲得滿足所有約束的最佳拓撲。

盡管圖3的最優拓撲結構與圖5相似，但在相同的收斂條件下，通過本研究所提出的低密度單元體積變化率約束方法，經過350步優化迭代可獲得結構的最優拓撲構型圖及相關參數，而不受低密度單元體積變化率約束的方法需要經過550步優化迭代才能獲得最終拓撲結構。顯然，使用本研究提出的低密度單元體積變化率約束可以極大地提高考慮自重影響的結構拓撲問題的優化計算效率，在本算例中對計算效率的提升高達36.3%。

圖6給出了使用該方法得到的結構拓撲灰度值md的優化歷程，其中兩條曲線分別代表考慮與不考慮低密度單元體積變化率約束時灰度值md的變化。由于Heaviside曲率參數的變化，圖6的曲線中有一些明顯的跳躍。圖7為在使用本研究所提出的方法后得到的橋梁結構柔順度優化歷程。圖8給出通過使用本研究提出的方法得到的橋梁頂部中點垂向位移的優化歷程，其最終位移值達到了相應約束的要求。對圖2至圖8分析可以發現，雖然圖3和圖5給出了幾乎相同的最優拓撲構型圖，但本研究所提出的考慮低密度單元的體積變化率約束方法可以顯著地抑制物理密度的波動，降低優化過程中中間密度單元對結構的影響，極大地提高拓撲優化計算效率。

4 結論

本研究基于SIMP方法，提出了一種考慮自重影響的連續體結構拓撲優化方法，并通過一些具體算例證明了所提出方法的可行性、有效性和特點，并得到以下3個結論。

①本研究所提出的改進SIMP模型可有效解決基于自重影響的拓撲優化過程中的寄生效應。

②本研究通過引入的體積上下限約束和指定點的位移約束在考慮自重的拓撲優化過程中發揮作用。

③引入的Heaviside三場過濾方案和本研究提出的一種低物理密度單元的體積變化率約束措施，可極大地提高考慮自重影響的結構拓撲問題的優化計算效率，并能得到更為清晰的0-1分布。

參考文獻：

[1] 姚錫凡，景軒，張劍銘，等.走向新工業革命的智能制造[J].計算機集成制造系統，2020（9）：2299-2320.

[2] BRUYNEEL M，DUYSINX P.Note on topology optimization of continuum structures including self-weight[J].Structural and Multidisciplinary Optimization，2005（4）：245-256.

[3] BOURDIN B.Filters in topology optimization[J].International Journal for Numerical Methods in Engineering，2001（9）：2143-2158.

[4] WANG F ，LAZAROV B S，Sigmund O.On projection methods，convergence and robust formulations in topology optimization[J].Structural and Multidisciplinary Optimization，2011（6）：767-784.

[5] Rong J H，Rong X P，Peng L，et al.A new method for optimizing the topology of hinge-free and fully decoupled compliant mechanisms with multiple inputs and multiple outputs[J].International Journal for Numerical Methods in Engineering，2021（1）：2863–2890.

[6] 何一凡，趙磊.基于應力約束的框架結構拓撲優化研究[J].河南科技，2021（28）：87-93.

[7] 陳成，趙圣佞.基于Heaviside過濾和可行域調整的SIMP方法拓撲優化設計[J].河南科技，2018（34）：26-28.

[8] 李方義，盧志宏，易繼軍.考慮自重影響和屈曲約束的結構拓撲優化方法[J].長沙理工大學學報（自然科學版），2021（3）：97-104.