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

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

關(guān)鍵詞：拓?fù)鋬?yōu)化;連續(xù)體結(jié)構(gòu);結(jié)構(gòu)自重;低密度約束

中圖分類號：TH122 ? ? 文獻(xiàn)標(biāo)志碼：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 引言

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

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

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

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

1.1 過濾方法

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

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

1.2 改進(jìn)的材料懲罰模型

SIMP法是現(xiàn)有拓?fù)鋬?yōu)化方法中最常用的一種方法。針對考慮自重影響的拓?fù)鋬?yōu)化問題，采用經(jīng)典的SIMP模型來解決時，在低密度區(qū)域，結(jié)構(gòu)剛度與其自重存在較大差異，使結(jié)構(gòu)剛度難以承受自重，最終導(dǎo)致結(jié)構(gòu)剛度不足和無界位移。本研究提出一種改進(jìn)的固體各向同性懲罰材料模型，見式（2），用來解決拓?fù)鋬?yōu)化中自重載荷引起的低密度區(qū)域的寄生效應(yīng)，同時保證體積約束的有效性。

式中：[E0]是初始單元材料的彈性模量;[ρ0]為單元初始質(zhì)量密度;[p]和[q]分別為單元的剛度模量懲罰因子和體積懲罰因子;[α]為剛度模量的多項式函數(shù)的懲罰參數(shù)。通過數(shù)值仿真研究發(fā)現(xiàn)，當(dāng)[p=5]、[q=2]、[α=25]時，對寄生效應(yīng)的抑制效果較好，且能獲得較好的拓?fù)鋬?yōu)化結(jié)果。

2 優(yōu)化模型建立及靈敏度分析

2.1 結(jié)構(gòu)自重處理和目標(biāo)函數(shù)

當(dāng)考慮與質(zhì)量密度相關(guān)的自重時，載荷向量表達(dá)式見式（3）。

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

式中：[F]為載荷向量;[g]為結(jié)構(gòu)自重向量;[P]為外載荷向量。

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

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

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

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

為抑制優(yōu)化過程中低密度單元的密度變量波動，并通過提升優(yōu)化過程中結(jié)構(gòu)清晰度來提高計算效率，本研究提出一個低密度單元體積變化率約束，見式（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]的表達(dá)式，[γ1]為一個經(jīng)驗參數(shù)，本研究取值為0.05。當(dāng)[LV≥VLTc]時，在后續(xù)優(yōu)化模型中引入低密度單元體積變化率約束。

式中：[N]為設(shè)計域中單元的數(shù)量;[Vfull]為結(jié)構(gòu)拓?fù)渲谐錆M材料時設(shè)計域的總體積;[γ2]為一個經(jīng)驗參數(shù)，本研究γ2取值為5。

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

2.3 位移約束

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

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

2.4 優(yōu)化模型

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

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

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

2.5 靈敏度分析

本研究采用MMA算法[7]來求解，因此，需要求解優(yōu)化模型式（13）中各個性能函數(shù)的靈敏度。由鏈?zhǔn)角髮?dǎo)法則可知，結(jié)構(gòu)的任意性能函數(shù)對設(shè)計變量的導(dǎo)數(shù)見式（18）。

由懲罰模型式（2）可得整體剛度矩陣關(guān)于物理密度變量[xi]的導(dǎo)數(shù)，見式（19）。

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

與結(jié)構(gòu)自重相關(guān)的等效節(jié)點力向量對物理密度變量[xi]的導(dǎo)數(shù)見式（20）和式（21）。

假設(shè)[Vj]為單位虛荷載作用在優(yōu)化結(jié)構(gòu)的第[nj]號自由度上所產(chǎn)生的結(jié)構(gòu)位移矢量，則結(jié)構(gòu)第[nj]號自由度的位移絕對值對物理密度變量[xi]的導(dǎo)數(shù)見式（22）。

基于式（5），優(yōu)化結(jié)構(gòu)柔順度對物理密度變量[xi]的導(dǎo)數(shù)見式（23）。

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

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

式（22）至式（26）是根據(jù)鏈?zhǔn)角髮?dǎo)法則得到的位移、柔順度和體積相關(guān)函數(shù)對設(shè)計變量的導(dǎo)數(shù)。

3 算例分析

為驗證本研究所提出的拓?fù)鋬?yōu)化方法的準(zhǔn)確性，針對圖1所示的考慮自重載荷的橋梁結(jié)構(gòu)進(jìn)行優(yōu)化求解。

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

圖2給出了使用本研究提出的方法得到的橋梁結(jié)構(gòu)拓?fù)渥兞糠植嫉膬?yōu)化歷程圖。其最優(yōu)拓?fù)錁?gòu)型如圖2（f）和圖3所示，最優(yōu)拓?fù)浣Y(jié)構(gòu)柔順度為113.84 N·m，橋頂中點垂向位移為0.002 78 m。圖4給出了在沒有考慮低密度單元的體積變化率約束時得到的拓?fù)渥兞糠植嫉膬?yōu)化歷程圖，其最優(yōu)拓?fù)浣Y(jié)構(gòu)如圖4（f）和圖5所示。最優(yōu)拓?fù)浣Y(jié)構(gòu)的柔順度和橋頂中點垂向位移分別為114.07 N·m和0.002 79 m。與圖2一樣，圖4中最黑的區(qū)域為物理密度變量大于或等于0.1的單元所占區(qū)域，從而突出顯示物理密度變量在[0，0.1]范圍內(nèi)的元素分布，可以有效地對比出低密度單元體積變化率約束的作用。這些數(shù)據(jù)以及最優(yōu)拓?fù)錁?gòu)型圖可以充分表明本研究所提出的方法能夠獲得滿足所有約束的最佳拓?fù)洹?/p>

盡管圖3的最優(yōu)拓?fù)浣Y(jié)構(gòu)與圖5相似，但在相同的收斂條件下，通過本研究所提出的低密度單元體積變化率約束方法，經(jīng)過350步優(yōu)化迭代可獲得結(jié)構(gòu)的最優(yōu)拓?fù)錁?gòu)型圖及相關(guān)參數(shù)，而不受低密度單元體積變化率約束的方法需要經(jīng)過550步優(yōu)化迭代才能獲得最終拓?fù)浣Y(jié)構(gòu)。顯然，使用本研究提出的低密度單元體積變化率約束可以極大地提高考慮自重影響的結(jié)構(gòu)拓?fù)鋯栴}的優(yōu)化計算效率，在本算例中對計算效率的提升高達(dá)36.3%。

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

4 結(jié)論

本研究基于SIMP方法，提出了一種考慮自重影響的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化方法，并通過一些具體算例證明了所提出方法的可行性、有效性和特點，并得到以下3個結(jié)論。

①本研究所提出的改進(jìn)SIMP模型可有效解決基于自重影響的拓?fù)鋬?yōu)化過程中的寄生效應(yīng)。

②本研究通過引入的體積上下限約束和指定點的位移約束在考慮自重的拓?fù)鋬?yōu)化過程中發(fā)揮作用。

③引入的Heaviside三場過濾方案和本研究提出的一種低物理密度單元的體積變化率約束措施，可極大地提高考慮自重影響的結(jié)構(gòu)拓?fù)鋯栴}的優(yōu)化計算效率，并能得到更為清晰的0-1分布。

參考文獻(xiàn)：

[1] 姚錫凡，景軒，張劍銘，等.走向新工業(yè)革命的智能制造[J].計算機集成制造系統(tǒng)，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] 何一凡，趙磊.基于應(yīng)力約束的框架結(jié)構(gòu)拓?fù)鋬?yōu)化研究[J].河南科技，2021（28）：87-93.

[7] 陳成，趙圣佞.基于Heaviside過濾和可行域調(diào)整的SIMP方法拓?fù)鋬?yōu)化設(shè)計[J].河南科技，2018（34）：26-28.

[8] 李方義，盧志宏，易繼軍.考慮自重影響和屈曲約束的結(jié)構(gòu)拓?fù)鋬?yōu)化方法[J].長沙理工大學(xué)學(xué)報（自然科學(xué)版），2021（3）：97-104.