基于應(yīng)力約束的框架結(jié)構(gòu)拓?fù)鋬?yōu)化研究

何一凡 趙磊

摘 要：針對(duì)工程框架結(jié)構(gòu)的強(qiáng)度需求，建立了應(yīng)力和體積約束下框架結(jié)構(gòu)柔順度最小化的拓?fù)鋬?yōu)化模型。首先，為解決應(yīng)力優(yōu)化過程中的應(yīng)力奇異和大量局部約束問題，利用qp應(yīng)力松弛技術(shù)和p范數(shù)凝聚函數(shù)法構(gòu)建了應(yīng)力約束的歸一化等效應(yīng)力約束方案。其次，提出了基于凝聚應(yīng)力和變體積約束限措施的修正方案，以克服優(yōu)化過程中應(yīng)力約束嚴(yán)重非線性和最大局部應(yīng)力波動(dòng)等問題。再次，導(dǎo)出了目標(biāo)函數(shù)和應(yīng)力約束的靈敏度公式，并采用移動(dòng)漸近線方法算法進(jìn)行優(yōu)化求解。最后，給出優(yōu)化算例，驗(yàn)證了本文方法的正確性與可行性。

關(guān)鍵詞：拓?fù)鋬?yōu)化;框架結(jié)構(gòu);應(yīng)力約束;應(yīng)力松弛;凝聚函數(shù)法

中圖分類號(hào)：TB21 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1003-5168（2021）28-00-07

Abstract： For the strength requirements of the engineering frame structures， this article establishes a topology optimization model that minimizes the flexibility of the frame structure under stress and volume constraints. In order to solve the stress singularity and deal with a large number of local constraints during the stress optimization process， a normalized equivalent stress constraint scheme for stress constraints is constructed by using the qp stress relaxation technique and the p-norm aggregation function method. Subsequently， a modified scheme based on the aggregation stress and the variable volume constraint limit measures is proposed to overcome the nonlinear problem of aggregation stress constraints and the maximum local stress fluctuation problem during an optimization process. The sensitivity formulas of the objective function and stress constraints are derived， and the MMA algorithm is adopted to optimize the model. Finally， the optimization example given verifies the correctness and feasibility of the proposed method.

Keywords： topology optimization;frame structure;stress constraint;stress relaxation;aggregate function method

拓?fù)鋬?yōu)化作為結(jié)構(gòu)創(chuàng)新設(shè)計(jì)的重要手段，已在機(jī)械、土木及航空航天等工程領(lǐng)域得到了廣泛應(yīng)用。目前，大多數(shù)研究仍關(guān)注連續(xù)體結(jié)構(gòu)剛度相關(guān)的優(yōu)化設(shè)計(jì)[1-4]，而強(qiáng)度問題是工程結(jié)構(gòu)設(shè)計(jì)需要考慮的首要條件之一。因此，考慮強(qiáng)度的結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)獲得了越來越多的關(guān)注[5-8]。

近年來，涉及強(qiáng)度需求的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)方法取得了較大進(jìn)展，但由于應(yīng)力約束的局部特性，優(yōu)化模型中存在大量的局部約束導(dǎo)致優(yōu)化求解困難，限制了相關(guān)研究的發(fā)展。為此，許多學(xué)者采用p范數(shù)或K-S（Kreisselmeier-Steinhauser）凝聚函數(shù)法對(duì)應(yīng)力約束進(jìn)行等效化處理[9]，以解決應(yīng)力優(yōu)化相關(guān)的大規(guī)模局部約束問題。然而，由于應(yīng)力約束本身具有非線性特征，加之結(jié)構(gòu)中尖端、凹槽等區(qū)域出現(xiàn)的應(yīng)力集中現(xiàn)象和應(yīng)力約束凝聚處理導(dǎo)致了應(yīng)力約束的強(qiáng)非線性問題，使得優(yōu)化結(jié)果對(duì)優(yōu)化模型及優(yōu)化求解器的參數(shù)高度敏感。為了規(guī)避優(yōu)化結(jié)構(gòu)中的應(yīng)力集中現(xiàn)象和解決具有強(qiáng)非線性優(yōu)化模型的求解問題，需要對(duì)凝聚參數(shù)合理取值，并修正優(yōu)化模型，以自適應(yīng)方式抑制設(shè)計(jì)變量大的變化。同時(shí)，應(yīng)力水平的精確評(píng)估與基于梯度的優(yōu)化算法的選用對(duì)克服該強(qiáng)非線性問題至關(guān)重要[10-11]。為此，張維聲等采用兩種全局應(yīng)力評(píng)估方法（包括應(yīng)力梯度和局部曲率信息）來控制結(jié)構(gòu)的局部應(yīng)力[12]。XIA等采用一種自適應(yīng)的有限元技術(shù)來避免優(yōu)化過程中出現(xiàn)的人工弱材料，以精確計(jì)算局部應(yīng)力水平[13]。PICELLI等提出一種應(yīng)力平均技術(shù)以提高應(yīng)力求解精確性[14]，并采用自適應(yīng)歸一化方法以準(zhǔn)確約束局部應(yīng)力。此外，優(yōu)化過程中低密度區(qū)域會(huì)引起應(yīng)力奇異現(xiàn)象，進(jìn)一步導(dǎo)致應(yīng)力相關(guān)優(yōu)化求解存在較大困難。針對(duì)該問題，許多學(xué)者提出了ε松弛和qp應(yīng)力松弛技術(shù)，并成功應(yīng)用到了應(yīng)力相關(guān)的拓?fù)鋬?yōu)化中[11]。但上述研究仍局限于連續(xù)體結(jié)構(gòu)的優(yōu)化設(shè)計(jì)，而工程中由梁、桿等組成的框架結(jié)構(gòu)應(yīng)用普遍。為推進(jìn)工程框架結(jié)構(gòu)的實(shí)際應(yīng)用，需要進(jìn)一步發(fā)展考慮應(yīng)力約束的框架結(jié)構(gòu)的拓?fù)鋬?yōu)化方法。

考慮到工程中框架結(jié)構(gòu)的強(qiáng)度需求，筆者提出了考慮應(yīng)力及變體積限約束的框架結(jié)構(gòu)柔順度最小化的拓?fù)鋬?yōu)化設(shè)計(jì)方法。為規(guī)避優(yōu)化結(jié)構(gòu)中的細(xì)長桿件，首先基于Heaviside映射函數(shù)構(gòu)建了設(shè)計(jì)變量與物理變量的合理映射關(guān)系。其次，采用qp應(yīng)力松弛技術(shù)和p范數(shù)凝聚函數(shù)對(duì)應(yīng)力約束進(jìn)行等效化處理，以克服應(yīng)力奇異問題和大量局部約束問題。再次，基于變體積限約束方法和凝聚函數(shù)的修正方案，建立了考慮應(yīng)力約束的框架結(jié)構(gòu)拓?fù)鋬?yōu)化的等效優(yōu)化模型。最后，通過移動(dòng)漸近線方法（Method of Moving Asymptotes，MMA）算法進(jìn)行了優(yōu)化求解，并通過兩個(gè)數(shù)值優(yōu)化算例說明了優(yōu)化設(shè)計(jì)方法的正確性與可行性。

1 框架結(jié)構(gòu)的懲罰模型和等效應(yīng)力

1.1 框架結(jié)構(gòu)的等效應(yīng)力

選取研究對(duì)象為由空心矩形截面的Euler-Bernoulli梁單元[如圖1（a）所示]組成的框架結(jié)構(gòu)。

對(duì)框架結(jié)構(gòu)而言，其構(gòu)件截面上應(yīng)力與截面中點(diǎn)的位置坐標(biāo)有關(guān)。將如圖1（b）所示的梁單元端面黑色標(biāo)注點(diǎn)作為該截面的臨危點(diǎn)[15]，采用臨危點(diǎn)應(yīng)力校核結(jié)構(gòu)強(qiáng)度設(shè)計(jì)的安全性。若b、h和t分別表示梁單元橫截面的寬度、高度及厚度參數(shù)，則單元橫截面積可表示為：

在如圖1（a）所示的局部坐標(biāo)系（x，y，z）下，梁單元的內(nèi)力可表示為：

式中：Ui、Ki和Fi分別為單元整體坐標(biāo)系中第i號(hào)單元節(jié)點(diǎn)位移矢量、在局部坐標(biāo)系下第i號(hào)單元的固有剛度矩陣及其單元內(nèi)力;Ti為局部坐標(biāo)系與整體坐標(biāo)系之間的變換矩陣，由單元局部坐標(biāo)軸和結(jié)構(gòu)整體坐標(biāo)軸的方向余弦構(gòu)成。在局部坐標(biāo)系下，第i號(hào)梁單元兩端截面內(nèi)力可以表示為：

式中：qi1和qi2分別為第i號(hào)梁單元兩端截面的內(nèi)力;r11=[I6×6 06×6]，r2=[06×6 I6×6]，I6×6和06×6分別為維度為6×6的單位矩陣及零矩陣;Fx，ij表示第i號(hào)梁單元j（j=1，2）截面的軸力;Fy，ij和Fz，ij分別為相應(yīng)截面平行于y和z軸的剪力;Mx，ij、My，ij和Mz，ij分別為相應(yīng)截面受到的繞x軸、y軸和z軸的彎矩。

梁單元截面上任意一點(diǎn)的應(yīng)力狀態(tài)為：

其中：

式中：σF是軸向力Fx引起的正應(yīng)力;σMy、σMz是彎矩My和Mz引起的彎曲應(yīng)力;τMy是由扭矩Mx引起的扭轉(zhuǎn)切應(yīng)力;τFy、τFz是由剪力Fy和Fz引起的彎曲切應(yīng)力。

第i號(hào)梁單元第j截面上臨危點(diǎn)的Von Mises等效應(yīng)力可表示為：

1.2 框架結(jié)構(gòu)的懲罰模型

采用固體各向同性材料懲罰模型（Solid Isotropic Material with Penalization，SIMP）插值模型，第i號(hào)單元的剛度矩陣Ki表示為：

式中：Ki0為第i號(hào)單元局部坐標(biāo)系中固有單元?jiǎng)偠染仃?p為懲罰參數(shù)。為了解決結(jié)構(gòu)單元灰度問題，采用Heaviside映射技術(shù)，形成設(shè)計(jì)變量與物理變量ρi映射的表達(dá)式：

式中：β為Heaviside映射曲率參數(shù)。借鑒相關(guān)文獻(xiàn)[11]的應(yīng)力奇異處理措施，可將單元體積vi、截面應(yīng)力σ（jvm，t）和應(yīng)力限分別表示為：

式中：vi0和表示第i號(hào)單元的固有體積和j（j=1，2）截面臨危點(diǎn)固有應(yīng)力;αv、p和q為懲罰參數(shù);為結(jié)構(gòu)固有的最大許用應(yīng)力。

2 優(yōu)化模型

考慮到工程框架結(jié)構(gòu)的強(qiáng)度需求，可構(gòu)建應(yīng)力及體積約束的框架結(jié)構(gòu)柔順度最小化的拓?fù)鋬?yōu)化模型，如式（17）所示：

式中：是拓?fù)湓O(shè)計(jì)變量向量，m為結(jié)構(gòu)中梁單元的數(shù)目;為第l組載荷工況下的自適應(yīng)權(quán)重系數(shù);Ctol為組合的目標(biāo)柔順度;Cl為第l組載荷工況下的柔順度，由式（19）給出;K為整體剛度矩陣;Ul為第l組載荷工況下整體坐標(biāo)系中結(jié)構(gòu)的位移矢量;Fl為第l個(gè)載荷工況;n為載荷工況數(shù);ρi和vi0分別為第i個(gè)物理變量和第i號(hào)單元的初始體積;V（0）為初始迭代時(shí)結(jié)構(gòu)的總體積;為在第l組載荷工況下第i號(hào)單元的j號(hào)截面的Von Mises等效應(yīng)力值;和分別為拓?fù)湓O(shè)計(jì)變量上下限，本文取;V*為目標(biāo)體積。自適應(yīng)權(quán)重系數(shù)的表達(dá)式為：

式中：Cl（w）表示第w迭代步第l（l=1，2，…，n）組載荷工況下的結(jié)構(gòu)柔順度;Cl為第l組載荷工況下的結(jié)構(gòu)柔順度，其可表示為：

為了解決大量局部應(yīng)力約束導(dǎo)致的計(jì)算量大的困難，基于qp應(yīng)力松弛技術(shù)和p范數(shù)凝聚函數(shù)法，將優(yōu)化模型式（17）中的大量局部應(yīng)力約束近似等效為下列凝聚應(yīng)力約束[15]。

式中：可以通過凝聚參數(shù)pn（pn>1）的值使得第l組載荷工況下的結(jié)構(gòu)凝聚應(yīng)力約束近似等效于第l組載荷工況下結(jié)構(gòu)大量局部應(yīng)力約束;表示第l組載荷工況下的結(jié)構(gòu)第i號(hào)單元的j（j=1，2）截面臨危點(diǎn)固有應(yīng)力;為第w迭代步的值。但是，pn取值過大將使凝聚應(yīng)力約束的非線性程度增大，會(huì)導(dǎo)致優(yōu)化求解過程反復(fù)震蕩;當(dāng)pn取值較小時(shí)，凝聚函數(shù)無法取代包絡(luò)對(duì)象的最大值。相關(guān)參考文獻(xiàn)[16]引入系數(shù)cp，使得式（20）為緊約束，其中：

3 靈敏度分析

3.1 結(jié)構(gòu)柔順度及體積的靈敏度

由優(yōu)化模型式（17）的結(jié)構(gòu)柔順度表達(dá)式可推出目標(biāo)函數(shù)Ctol關(guān)于物理變量ρk的靈敏度：

式中：Ul（w）為第l組載荷工況下整體坐標(biāo)系中結(jié)構(gòu)第w步的位移矢量。

3.2 應(yīng)力約束的靈敏度

由鏈?zhǔn)角髮?dǎo)法則可得凝聚應(yīng)力約束對(duì)物理變量的靈敏度為：

基于上述靈敏度分析，結(jié)合MMA算法，對(duì)本文所提出的優(yōu)化模型進(jìn)行優(yōu)化求解，當(dāng)滿足如式（37）和式（38）所示的收斂條件時(shí)則終止迭代，得到優(yōu)化解。為了進(jìn)一步解決嚴(yán)重的非線性求解問題，在MMA迭代求解過程中采用xml對(duì)設(shè)計(jì)變量進(jìn)行控制。

式中：m的初始值為1，并且其值每50個(gè)外循環(huán)迭代步增加1。

式中：和分別為第w步和第w-1步設(shè)計(jì)變量的值;md為灰度指標(biāo);ε1、ε2為兩個(gè)較小的經(jīng)驗(yàn)參數(shù)值，本文取ε1=0.01，ε2=0.01。

4 算例分析

一高聳空心結(jié)構(gòu)高80 m，最底部外尺寸為4 m×4 m的正方形區(qū)域，最頂部外尺寸為直徑4 m的圓形區(qū)域，采用TOMáS等[17]提出的方式生成該模型。中間空心區(qū)域是底面半徑為3 m、高80 m的圓柱形區(qū)域。初始框架基結(jié)構(gòu)中，梁構(gòu)件的截面均為0.1 m×0.1 m且壁厚0.01 m的空心矩形截面。結(jié)構(gòu)底面固定且受到如圖2所示的兩組方向相反的1 000 kN外載荷作用，其中F1沿y軸的反方向（工況1），F(xiàn)2沿y軸的正方向（工況2）。結(jié)構(gòu)的彈性模量E=210 GPa，泊松比v=0.3，結(jié)構(gòu)的初始最大應(yīng)力為3.974 1 MPa。本算例中取SIMP懲罰參數(shù)為p=3，應(yīng)力約束限的懲罰參數(shù)q=2.5，凝聚參數(shù)pn=8，目標(biāo)體積為0.06，變體積限的移動(dòng)步長為0.02。

結(jié)構(gòu)的初始基結(jié)構(gòu)設(shè)計(jì)域如圖3所示，由1 656個(gè)梁單元構(gòu)成。采取4種不同的應(yīng)力約束方案。方案一=3.2 MPa，方案二=3.5 MPa，方案三=5 MPa，方案四未考慮應(yīng)力約束。

如圖2所示，在工況1和工況2的共同作用下，結(jié)構(gòu)在不同應(yīng)力約束限下的優(yōu)化拓?fù)浼捌鋺?yīng)力云圖分別如圖4至圖11所示?？芍蚣芙Y(jié)構(gòu)的優(yōu)化拓?fù)潆S應(yīng)力約束限的不同產(chǎn)生細(xì)微變化，而相應(yīng)的應(yīng)力云圖隨應(yīng)力約束限的不同產(chǎn)生較大變化。由圖4至圖11的對(duì)比分析還可以看出：有無考慮應(yīng)力約束的優(yōu)化結(jié)果存在明顯不同。

結(jié)構(gòu)在不同約束限下的柔順度、最大Von Mises應(yīng)力及結(jié)構(gòu)體積的迭代歷程曲線分別如圖12至圖15所示，可知：優(yōu)化過程穩(wěn)定收斂，且其優(yōu)化結(jié)構(gòu)的柔順度和最大應(yīng)力都隨桿件減少而增大。對(duì)應(yīng)圖4至圖9中優(yōu)化構(gòu)型的最大應(yīng)力分別為2.880 9 MPa、2.6345 MPa及2.529 7 MPa，其最終柔順度分別為10.284 6 N·m、11.231 6 N·m和10.435 1 N·m。圖10和圖11中未考慮應(yīng)力約束的優(yōu)化構(gòu)型的最大應(yīng)力為3.014 2 MPa，其最終柔順度為11.088 7 N·m。由此可知，優(yōu)化結(jié)構(gòu)的最大應(yīng)力也隨不同應(yīng)力約束限而不同，且優(yōu)化結(jié)構(gòu)的最大應(yīng)力隨應(yīng)力約束限的增大而增大，而其柔順度未發(fā)生明顯變化。同時(shí)，考慮應(yīng)力約束的優(yōu)化結(jié)構(gòu)的最大應(yīng)力明顯小于未考慮應(yīng)力約束的優(yōu)化結(jié)果?？梢?，提出的方法有效控制了優(yōu)化結(jié)構(gòu)的最大應(yīng)力，得到的拓?fù)鋬?yōu)化結(jié)果是滿足應(yīng)力約束且合理的。

5 結(jié)論

針對(duì)框架結(jié)構(gòu)的強(qiáng)度需求，基于SIMP插值模型和Heaviside映射函數(shù)建立了應(yīng)力及體積約束下框架結(jié)構(gòu)輕量化拓?fù)鋬?yōu)化模型。采用qp應(yīng)力松弛技術(shù)和p范數(shù)凝聚函數(shù)法構(gòu)建了等效應(yīng)力約束函數(shù)，通過凝聚應(yīng)力約束修正方案和變體積約束限方法形成了相應(yīng)的等效優(yōu)化模型。此外，進(jìn)一步解決了嚴(yán)重的非線性模型求解問題，對(duì)每一步設(shè)計(jì)變量的變化范圍進(jìn)行了控制。最后，通過數(shù)值算例來說明所提出優(yōu)化方法的正確性與可行性，并得到如下結(jié)論。

①本方法能獲得穩(wěn)定收斂的解，且能得到滿足應(yīng)力約束的清晰優(yōu)化拓?fù)洹?/p>

②結(jié)構(gòu)的優(yōu)化構(gòu)型隨應(yīng)力約束限的不同而發(fā)生變化，且其柔順度和最大應(yīng)力都隨桿件的消除而增大。

③考慮應(yīng)力約束的優(yōu)化結(jié)構(gòu)的最大應(yīng)力明顯低于未考慮應(yīng)力約束的最大應(yīng)力。提出的方法有效控制了優(yōu)化結(jié)構(gòu)的最大應(yīng)力，得到的拓?fù)鋬?yōu)化結(jié)果滿足應(yīng)力約束且合理。

參考文獻(xiàn)：

[1]DUYSINX P， BENDSOE M P. Topology optimization of continuum structures with stress constraints[J]. International Journal for Numerical Methods in Engineering， 1998（8）：1453-1478.

[2]ZHU J H， GUO W J ， ZHANG W H， et al. Integrated layout and topology optimization design of multi-frame and multi-component fuselage structure systems[J]. Structural and Multidisciplinary Optimization， 2017（1）：21-45.

[3]冷國俊， 張卓， 保宏， 等.考慮重疊過濾及穩(wěn)定性約束的桁架拓?fù)鋬?yōu)化方法[J].工程力學(xué)， 2013（2）：8-13.

[4]CUI H， AN H， HUANG H . Truss topology optimization considering local buckling constraints and restrictions on intersection and overlap of bar members[J]. Structural and Multidisciplinary Optimization， 2018（2）：575-594.

[5]CHANGIZI N， KABOO DA NIAN H， JALALPOUR M. Stress-based topology optimization of frame structures under geometric uncertainty[J]. Computer Methods in Applied Mechanics and Engineering， 2017， 315：121-140.

[6]JEONG S H， CHOI D H， YOON G H . Separable stress interpolation scheme for stress-based topology optimization with multiple homogenous materials[J]. Finite Elements in Analysis & Design， 2014（5）：16-31.

[7]GAO X J， LI YX， MA H T， et al. Improving the overall performance of continuum structures： a topology optimization model considering stiffness， strength and stability[J]. Computer Methods in Applied Mechanics and Engineering， 2020，359：112660.

[8]LE C， NORATO J， BRUNS T， et al. Stress-based topology optimization for continua[J]. Structural and Multidisciplinary Optimization， 2010（4）：605-620.

[9]王選， 劉宏亮， 龍凱， 等.基于改進(jìn)的雙向漸進(jìn)結(jié)構(gòu)優(yōu)化法的應(yīng)力約束拓?fù)鋬?yōu)化[J].力學(xué)學(xué)報(bào)， 2018（2）：385-394.

[10]PARíS J， NAVARRINA F， COLOMINAS I， et al. Topology optimization of continuum structures with local and global stress constraints[J]. Chinese Journal of Theoretical & Applied Mechanics， 2009（4）：98-104.

[11]MOON S， YOON G. A newly developed qp-relaxation method for element connectivity parameterization to achieve stress-based topology optimization for geometrically nonlinear structures [J]. Computer Methods in Applied Mechanics and Engineering， 2013（3）：226-241.

[12]ZHANG W， GUO X， WANG M， et al. Optimal topology design of continuum structures with stress concentration alleviation via level set method [J]. International Journal for Numerical Methods in Engineering， 2013（9）：942-959.

[13]XIA Q， SHI T， LIU S， et al. A level set solution to the stress-based structural shape and topology optimization [J]. Computers & Structures， 2012（1）：55-64.

[14]PICELLI R， TOWNSEND S， BRAMPTON C， et al. Stress-based shape and topology optimization with the level set method [J]. Computer Methods in Applied Mechanics and Engineering， 2018（1）：1-23.

[15]ZUO W， YU J， SAITOU K. Stress sensitivity analysis and optimization of automobile body frame consisting of rectangular tubes[J]. International Journal of Automotive Technology， 2016（5）：843-851.

[16]龍凱， 王選， 吉亮.面向應(yīng)力約束的獨(dú)立連續(xù)映射方法[J].力學(xué)學(xué)報(bào)， 2019（2）：320-329.

[17]TOMáS ZEGARD， PAULINO G H . GRAND3— Ground structure based topology optimization for arbitrary 3D domains using MATLAB[J]. Structural & Multidisciplinary Optimization， 2015（6）：1161-1184.

3424500338229