基于應力約束的框架結構拓撲優化研究

何一凡 趙磊

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

關鍵詞：拓撲優化;框架結構;應力約束;應力松弛;凝聚函數法

中圖分類號：TB21 文獻標識碼：A 文章編號：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

拓撲優化作為結構創新設計的重要手段，已在機械、土木及航空航天等工程領域得到了廣泛應用。目前，大多數研究仍關注連續體結構剛度相關的優化設計[1-4]，而強度問題是工程結構設計需要考慮的首要條件之一。因此，考慮強度的結構拓撲優化設計獲得了越來越多的關注[5-8]。

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

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

1 框架結構的懲罰模型和等效應力

1.1 框架結構的等效應力

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

對框架結構而言，其構件截面上應力與截面中點的位置坐標有關。將如圖1（b）所示的梁單元端面黑色標注點作為該截面的臨危點[15]，采用臨危點應力校核結構強度設計的安全性。若b、h和t分別表示梁單元橫截面的寬度、高度及厚度參數，則單元橫截面積可表示為：

在如圖1（a）所示的局部坐標系（x，y，z）下，梁單元的內力可表示為：

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

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

梁單元截面上任意一點的應力狀態為：

其中：

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

第i號梁單元第j截面上臨危點的Von Mises等效應力可表示為：

1.2 框架結構的懲罰模型

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

式中：Ki0為第i號單元局部坐標系中固有單元剛度矩陣;p為懲罰參數。為了解決結構單元灰度問題，采用Heaviside映射技術，形成設計變量與物理變量ρi映射的表達式：

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

式中：vi0和表示第i號單元的固有體積和j（j=1，2）截面臨危點固有應力;αv、p和q為懲罰參數;為結構固有的最大許用應力。

2 優化模型

考慮到工程框架結構的強度需求，可構建應力及體積約束的框架結構柔順度最小化的拓撲優化模型，如式（17）所示：

式中：是拓撲設計變量向量，m為結構中梁單元的數目;為第l組載荷工況下的自適應權重系數;Ctol為組合的目標柔順度;Cl為第l組載荷工況下的柔順度，由式（19）給出;K為整體剛度矩陣;Ul為第l組載荷工況下整體坐標系中結構的位移矢量;Fl為第l個載荷工況;n為載荷工況數;ρi和vi0分別為第i個物理變量和第i號單元的初始體積;V（0）為初始迭代時結構的總體積;為在第l組載荷工況下第i號單元的j號截面的Von Mises等效應力值;和分別為拓撲設計變量上下限，本文取;V*為目標體積。自適應權重系數的表達式為：

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

為了解決大量局部應力約束導致的計算量大的困難，基于qp應力松弛技術和p范數凝聚函數法，將優化模型式（17）中的大量局部應力約束近似等效為下列凝聚應力約束[15]。

式中：可以通過凝聚參數pn（pn>1）的值使得第l組載荷工況下的結構凝聚應力約束近似等效于第l組載荷工況下結構大量局部應力約束;表示第l組載荷工況下的結構第i號單元的j（j=1，2）截面臨危點固有應力;為第w迭代步的值。但是，pn取值過大將使凝聚應力約束的非線性程度增大，會導致優化求解過程反復震蕩;當pn取值較小時，凝聚函數無法取代包絡對象的最大值。相關參考文獻[16]引入系數cp，使得式（20）為緊約束，其中：

3 靈敏度分析

3.1 結構柔順度及體積的靈敏度

由優化模型式（17）的結構柔順度表達式可推出目標函數Ctol關于物理變量ρk的靈敏度：

式中：Ul（w）為第l組載荷工況下整體坐標系中結構第w步的位移矢量。

3.2 應力約束的靈敏度

由鏈式求導法則可得凝聚應力約束對物理變量的靈敏度為：

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

式中：m的初始值為1，并且其值每50個外循環迭代步增加1。

式中：和分別為第w步和第w-1步設計變量的值;md為灰度指標;ε1、ε2為兩個較小的經驗參數值，本文取ε1=0.01，ε2=0.01。

4 算例分析

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

結構的初始基結構設計域如圖3所示，由1 656個梁單元構成。采取4種不同的應力約束方案。方案一=3.2 MPa，方案二=3.5 MPa，方案三=5 MPa，方案四未考慮應力約束。

如圖2所示，在工況1和工況2的共同作用下，結構在不同應力約束限下的優化拓撲及其應力云圖分別如圖4至圖11所示。可知，框架結構的優化拓撲隨應力約束限的不同產生細微變化，而相應的應力云圖隨應力約束限的不同產生較大變化。由圖4至圖11的對比分析還可以看出：有無考慮應力約束的優化結果存在明顯不同。

結構在不同約束限下的柔順度、最大Von Mises應力及結構體積的迭代歷程曲線分別如圖12至圖15所示，可知：優化過程穩定收斂，且其優化結構的柔順度和最大應力都隨桿件減少而增大。對應圖4至圖9中優化構型的最大應力分別為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中未考慮應力約束的優化構型的最大應力為3.014 2 MPa，其最終柔順度為11.088 7 N·m。由此可知，優化結構的最大應力也隨不同應力約束限而不同，且優化結構的最大應力隨應力約束限的增大而增大，而其柔順度未發生明顯變化。同時，考慮應力約束的優化結構的最大應力明顯小于未考慮應力約束的優化結果。可見，提出的方法有效控制了優化結構的最大應力，得到的拓撲優化結果是滿足應力約束且合理的。

5 結論

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

①本方法能獲得穩定收斂的解，且能得到滿足應力約束的清晰優化拓撲。

②結構的優化構型隨應力約束限的不同而發生變化，且其柔順度和最大應力都隨桿件的消除而增大。

③考慮應力約束的優化結構的最大應力明顯低于未考慮應力約束的最大應力。提出的方法有效控制了優化結構的最大應力，得到的拓撲優化結果滿足應力約束且合理。

參考文獻：

[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]冷國俊， 張卓， 保宏， 等.考慮重疊過濾及穩定性約束的桁架拓撲優化方法[J].工程力學， 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].力學學報， 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]龍凱， 王選， 吉亮.面向應力約束的獨立連續映射方法[J].力學學報， 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