基于改進BESO算法的帶隔板薄壁方管耐撞性拓撲優(yōu)化

王章駿，許平，邢杰，趙紫亮，蔡雋堃

(1.中南大學 交通運輸工程學院，湖南 長沙410075；2.中南大學 軌道交通安全教育部重點實驗室，湖南 長沙410075)

交通安全一直以來都是人們關注的熱點話題。吸能裝置作為被動防護的關鍵組件之一，在發(fā)生碰撞事故時，具有保護乘員的重要作用。薄壁管有著優(yōu)越的吸能特性，已經(jīng)廣泛地作為吸能裝置被應用到各種車輛中[1]。在薄壁管中加入橫隔板是改善吸能管吸能特性的有效方式之一。唐智亮[2]受到竹子節(jié)隔膜的啟發(fā)，結(jié)合仿生學的思想，將橫隔板加入到了非凸薄壁管中，并研究了隔板數(shù)量對非凸薄壁管能量吸收性能的影響。GAO等[3]利用實驗與仿真方法研究了帶隔板方管的耐撞性，研究表明，隔板能改變薄壁管的變形模式，提高其吸能量。近來，由于帶隔板吸能管優(yōu)秀的吸能特性，將其作為吸能裝置的軌道車輛出現(xiàn)了持續(xù)增長，因此，對帶隔板方管的優(yōu)化設計產(chǎn)生了迫切的需求。XU等[4]研究一種帶橫隔板的嵌套管并利用遺傳算法優(yōu)化了關鍵的結(jié)構(gòu)參數(shù)，提高了吸能管的耐撞性。在隨后的研究中，XU等[5]將橫隔板插入到軸對稱的方管中，結(jié)論表明，隔板可以讓這種新型結(jié)構(gòu)產(chǎn)生穩(wěn)定的變形，對特定的幾何參數(shù)進行優(yōu)化設計，可以提高吸能管的耐撞性。回顧過去優(yōu)化薄壁管幾何結(jié)構(gòu)的相關研究，大多數(shù)都是作者根據(jù)經(jīng)驗提出優(yōu)化參數(shù)，或者根據(jù)經(jīng)驗預設結(jié)構(gòu)形狀。拓撲優(yōu)化是結(jié)構(gòu)優(yōu)化的一種，能夠找到材料在指定條件下的最優(yōu)分布，比之傳統(tǒng)結(jié)構(gòu)優(yōu)化方法基于設計師的創(chuàng)造力和經(jīng)驗，以及反復試驗的過程，拓撲優(yōu)化允許創(chuàng)造性和根本性的設計修改[6]，且不需要設計師在原結(jié)構(gòu)上進行反復設計。具體到薄壁管的優(yōu)化設計，已經(jīng)有大量研究表明，在薄壁管的適當位置引入合適形狀的通孔可以改善薄壁管的耐撞性能[7?9]。傳統(tǒng)結(jié)構(gòu)優(yōu)化方法只能根據(jù)經(jīng)驗預定義通孔的位置和幾何形狀再進行優(yōu)化，這大大局限了薄壁管的優(yōu)化設計。而使用拓撲優(yōu)化方法，可以有效地實現(xiàn)材料分布的自動化設計，從而實現(xiàn)合理的通孔分布，提高薄壁管的耐撞性能。盡管拓撲優(yōu)化在諸多工程領域已經(jīng)有了長足的發(fā)展，但其應用于結(jié)構(gòu)耐撞性的研究起步較晚。HUANG等[10]提出將單位體積的能量吸收最大化作為優(yōu)化目標，并以此進行了實心簡支梁的吸能優(yōu)化，取得了良好效果。但這項工作沒有考慮單元平面外的屈曲，這與薄壁管的壓潰模式相去甚遠。DUDDECK等[11]提出了一種改進的HCA算法，將其用于薄壁管結(jié)構(gòu)，得到了內(nèi)墻的最佳幾何分布與厚度分布。但這項研究沒有將管壁考慮在內(nèi)，不均勻的厚度分布也給制造帶來困難。BAHRAMIAN等[12]用改進的BESO算法對簡單方管進行了拓撲優(yōu)化，得到了最佳的沿管壁的材料分布，本文的研究在其基礎上進行。本文關注于帶橫隔板的薄壁方管(the thin-walled square tube with diaphragms，以下簡稱TSTD)，為了在滿足其總吸能量要求的基礎上實現(xiàn)輕量化的設計目標，首次嘗試采用拓撲優(yōu)化來解決這一問題，得到其管壁和隔板的最佳材料分布。由于TSTD管在壓潰過程中隔板和管壁的不同作用，現(xiàn)有算法不能滿足TSTD拓撲優(yōu)化設計的需要，本文對BESO算法主要做出了以下改進：1)提出了多設計域優(yōu)化目標分離的BESO算法，實現(xiàn)隔板與管壁的同步優(yōu)化，滿足兩者的不同優(yōu)化目標，提高計算效率；2)提出了克隆體法，讓所有迭代步的模型信息完美保存，克服了傳統(tǒng)拓撲優(yōu)化算法在同一模型上修改導致的迭代步信息覆蓋缺陷。

1 初始結(jié)構(gòu)和有限元模型

1.1 初始幾何結(jié)構(gòu)

本文所研究的TSTD如圖1所示。主體結(jié)構(gòu)由2塊U形板焊接而成，形成管壁；兩端分別焊接端板，內(nèi)部均勻分布著4塊橫隔板。TSTD管段總長為504 mm，方管寬為150 mm，端板也為正方形，寬180 mm。隔板為缺了4個角的正方形，總寬為150 mm，與管壁焊接的長度為120 mm，沿軸向均勻分布，兩兩間隔100 mm。所有的板材厚度均為2 mm。

圖1 原始幾何模型Fig.1 Original geometric model

1.2 材料性質(zhì)

本文研究的實驗試件材料全為Q345低碳鋼。為了得到更加準確的仿真結(jié)果，進行了準靜態(tài)拉伸試驗來獲得材料的性質(zhì)。本次試驗所用的試驗機為MTS646液壓機，試驗方案依據(jù)國標GB/T 228.1—2010擬定。試驗共進行3次，計算3次結(jié)果的平均值，得到工程應力應變曲線。將工程應力應變曲線轉(zhuǎn)化為真實應力應變曲線。試驗過程及結(jié)果如圖2所示。所用到的Q345低碳鋼的其他材料參數(shù)參見表1。

表1 材料參數(shù)Table 1 Mechanical parameters

圖2 材料試驗與特性Fig.2 Material tests and properties

1.3 有限元模型及實驗驗證

本文的有限元模型用商業(yè)軟件ABAQUS建立。對于TSTD的管壁和隔板，采用減縮積分的殼單元(SR4)來模擬，這種單元不僅可以節(jié)省計算資源，而且已經(jīng)被證實可以有效地模擬薄壁管的壓潰模式[2]。根據(jù)之前章節(jié)描述的幾何模型，建立的有限元模型如圖3所示。除TSTD之外，加入了一塊固定的剛性墻，使TSTD遠離剛性墻一側(cè)的端板以0.5 m/s的速度撞擊剛性墻，來模擬TSTD的準靜態(tài)壓潰過程。特別的，在最初的0.08 s的時間內(nèi)，端板的從0 m/s開始平滑過渡到0.5 m/s，這樣可以減少動態(tài)效應影響仿真的準確性。有限元模型采用通用接觸，靜態(tài)摩擦因數(shù)為0.3，動態(tài)摩擦因數(shù)為0.2。

圖3 有限元模型Fig.3 Finite element model

為了讓有限元模型的網(wǎng)格劃分合理，進行了網(wǎng)格的獨立性驗證。當縮減網(wǎng)格尺寸至5 mm后，再縮減尺寸，有限元結(jié)果無較大變化。選定網(wǎng)格大小為5 mm。

為保證有限元仿真的有效性與準確性，對有限元仿真進行了實驗驗證。實驗采用INSTRON萬能試驗機在壓縮速度為10 mm/min的條件下進行。從圖4(a)的變形模式對比上來看，仿真與試驗都發(fā)生了相同的變形模式，產(chǎn)生的褶皺數(shù)量也相同。從圖4(b)力位移曲線的對比上來看，仿真得到的曲線很好地貼合了試驗結(jié)果。綜上所述，仿真與試驗的結(jié)果相吻合，驗證了仿真的有效性。

由于文本所研究的結(jié)構(gòu)具有對稱性，為節(jié)約計算資源，減少計算量，采用1/4有限元模型進行計算。模型中A面與C面相對，B面與D面相對。保留模型的相鄰A面與B面的一半和1～4號隔板的1/4進行拓撲優(yōu)化，如圖3中(b)所示。

2 算法及其實現(xiàn)

2.1 問題描述

BESO是最常用的拓撲優(yōu)化算法之一，算法依據(jù)單元的靈敏度信息(Sensitivity Number)，在每一個迭代步中，根據(jù)進化率(Evolutionary Ratio，簡稱ER)所得到的目標體積，對單元進行刪除或添加。優(yōu)化過程中的目標體積定義如下：

式中：vk+1表示下一個迭代步中的目標體積；vk表示當前迭代步中的結(jié)構(gòu)體積。

本文研究的問題陳述如下，在TSTD軸向壓潰至原長度的75%的條件下，最優(yōu)化TSTD的吸能效率。與普通薄壁管不同的是，TSTD中有起控制薄壁管變形的橫隔板存在。橫隔板的作用并不是能量吸收，在薄壁管的壓潰過程中不會發(fā)生大的塑性變形[2]，因此，最大化單位體積的能量吸收不宜作為橫隔板的優(yōu)化目標。

為此，在BAHRAMIAN研究的基礎上，本文提出了多設計域優(yōu)化目標分離的BESO算法來解決這一問題。即STD的管壁與橫隔板視為2個設計域，為兩者設置不同的優(yōu)化目標和不同的進化率，分開計算不同區(qū)域的靈敏度信息。但考慮大變形有限元分析帶來的巨大計算量，將兩者的拓撲優(yōu)化同時進行。數(shù)學模型陳述如下。

對于管壁：

式中：E為吸能量；V為體積；Umax為最大壓潰長度，相應地，U*為最大壓潰長度約束，本文中設置為TSTD原長的75%，即375 mm。Etarget為吸能約束，設置為原結(jié)構(gòu)吸能量的90%。

加入能量約束是為了維持吸能結(jié)構(gòu)的吸能量，以克服單元的刪除/添加可能會誘發(fā)薄壁管產(chǎn)生新的壓潰模式，從而導致總吸能量發(fā)生較大變化的難題。將E/V作為優(yōu)化目標同時考慮了能量吸收和結(jié)構(gòu)體積，可以使優(yōu)化具有更好的魯棒性[12]。

對于橫隔板：

式中：C為結(jié)構(gòu)柔度；t積分終止時刻，對應于最終狀態(tài)；F和u?分別為載荷向量與速度向量；V*為結(jié)構(gòu)體積約束；vj為第j個單元的體積。

將最小化結(jié)構(gòu)柔度，即最大化結(jié)構(gòu)剛度作為橫隔板的優(yōu)化目標是由橫隔板的作用決定的。為了讓橫隔板能夠影響薄壁管的壓潰模式，增強薄壁管漸進壓潰的穩(wěn)定性，就必須保證其有足夠的剛度來抵抗變形。因此對于橫隔板的體積分數(shù)而言，設置過小，將會使橫隔板的強度不足；太大又會使得拓撲優(yōu)化的輕量化設計效果不明顯。保守考慮后，設置橫隔板的體積優(yōu)化目標V*為原結(jié)構(gòu)體積的0.5。

2.2 更新準則

BESO算法是根據(jù)單元的靈敏度信息來更新材料的分布的。本文對管壁和橫隔板采取了不同的優(yōu)化目標，因此，也存在不同的靈敏度表達。

對于管壁：

式中：αi，Vi，Ei分別代表第i個單元的靈敏度，體積和應變能。V與E代表結(jié)構(gòu)的總體積和總應變能，詳細推導過程參見文獻[13]。

對于橫隔板：

為了改善棋盤格現(xiàn)象，減小目標函數(shù)振蕩的行為，算法也加入了靈敏度過濾技術(shù)和靈敏度進行歷史平均[14]。

2.3 計算流程

如圖3(b)所示，本文將計算模型分為3個區(qū)域，端板區(qū)(圖中白色區(qū)域)，管壁區(qū)(圖中藍色區(qū)域)和隔板區(qū)(圖中黃色區(qū)域)。其中端板區(qū)參與有限元計算但不參與優(yōu)化計算，管壁區(qū)采用針對管壁的優(yōu)化目標，隔板區(qū)采用針對橫隔板的優(yōu)化目標。

本文使用的BESO算法與傳統(tǒng)BESO算法另一個不同之處在于迭代步驟的終止條件不是達到收斂條件，而是達到預定義的最大循環(huán)次數(shù)。這是因為采用了殼單元雖然很好地反應了薄壁管壓潰的過程，但也給拓撲優(yōu)化帶來了挑戰(zhàn)，即在拓撲優(yōu)化過程中，對管壁引入新的材料空洞所導致的壓潰形式不可預測，因此BESO的優(yōu)化過程中的目標曲線不太可能會平滑過渡，拓撲優(yōu)化的終止步也不太可能是最優(yōu)步。相關問題在文獻[12]中也有詳盡的闡述。

由于符合約束條件的相對最優(yōu)結(jié)構(gòu)可能出現(xiàn)在拓撲優(yōu)化的過程中，這就需要保存所有迭代步的模型信息。為了適應這一要求，本文提出了克隆體法，具體策略闡述如下：在完成單元靈敏度的更新后，與傳統(tǒng)方法不同，不在原有模型上進行單元的刪除或添加，而是生成上一個迭代步的克隆體，克隆體繼承上一個迭代步的所有信息。然后將上一個迭代步的模型信息保存，在新生成的克隆體上進行單元的刪除或添加工作。在完成單元變更后，將克隆體模型作為本次迭代步的最終模型，并傳遞到下一個迭代步中。改進BESO算法的具體流程如圖5所示。

圖5 程序流程圖Fig.5 Program flow chart

程序計算所要設定的參數(shù)統(tǒng)計如表2。

表2 相關參數(shù)Table 2 Relevant parameters

3 計算結(jié)果與討論

使用ABAQUS運行Python代碼后得到結(jié)果，迭代歷史記錄如圖6所示。

必須說明的是，由于兩端板并未參與拓撲優(yōu)化，所以迭代歷史的統(tǒng)計計算中未加入兩端板的材料體積。由圖6(a)可知，隨著迭代的進行，結(jié)構(gòu)的體積分數(shù)在ER的控制下穩(wěn)步下降，并在第9個迭代步時出現(xiàn)第1個拐點，這是因為隔板的體積分數(shù)已經(jīng)達到目標值0.5，隔板的體積分數(shù)停止下降；第2個拐點出現(xiàn)在第29個迭代步，此時，管壁的體積分數(shù)也下降到了給定的界限值0.75，在往后的迭代步中，其值不出現(xiàn)下降。

圖6 迭代歷史記錄Fig.6 Iteration history

目標函數(shù)的迭代歷程如圖6(b)所示。為了更加直觀地與原始結(jié)構(gòu)作對比，將原始結(jié)構(gòu)的單位體積吸能量設為單位1。從圖6中可以看到，在整個迭代過程中，目標函數(shù)的值不斷震蕩，這種現(xiàn)象符合前文所闡述的變形模式的不確定性。但從整體趨勢來看，隨著迭代的進行，其值呈現(xiàn)上升趨勢，這表明了優(yōu)化程序的有效性。在整個迭代的過程中，相對體積吸能量最大的迭代步出現(xiàn)在第27步，是初始結(jié)構(gòu)的1.39倍。圖6(c)顯示了所有迭代步的能量吸收量。從圖6(c)可以觀察到，引入材料空洞，TSTD的吸能量不一定下降，如第9個迭代步，其吸能量是初始結(jié)構(gòu)的1.03倍，這也佐證了合理的材料空洞分布可能引導更有效的壓潰模式。目標函數(shù)表現(xiàn)最好的第27個迭代步的吸能量也在紅色虛線之上，滿足吸能量的約束。綜合判斷，將第27個迭代步作為最終的拓撲優(yōu)化結(jié)果。

將第27個迭代步的模型取出后進行對稱操作，將模型還原，最終結(jié)構(gòu)如圖7所示。但是優(yōu)化直接得到的拓撲結(jié)構(gòu)將會給工程制造帶來困難，必須對空洞特征做平滑處理，并刪除一些微小特征，修繕后的工程結(jié)構(gòu)如圖8所示。將工程結(jié)構(gòu)進行CAD建模后也導入ABAQUS進行壓潰仿真分析。

圖7 最終拓撲結(jié)構(gòu)Fig.7 Final topological structure

圖8 工程結(jié)構(gòu)Fig.8 Engineering structure

為了更好地評價TSTD的耐撞性能，將評價指標比吸能[15](specific energy absorption,SEA)引入，表達式為：

式中：Et是吸能管吸收的能量；m是吸能管的整體質(zhì)量。

統(tǒng)計原始結(jié)構(gòu)、拓撲結(jié)構(gòu)與工程結(jié)構(gòu)的性能指標在表3中。

表3 不同結(jié)構(gòu)的性能比較Table 3 Performance comparison of different structures

為了了解耐撞性產(chǎn)生變化的原因，將原始結(jié)構(gòu)、拓撲結(jié)構(gòu)與工程結(jié)構(gòu)的壓潰過程做對比，變形模式對比如圖9所示，力位移曲線對比如圖10所示。

從圖10可以觀察到拓撲構(gòu)型和工程結(jié)構(gòu)的初始峰值力比之原始結(jié)構(gòu)出現(xiàn)明顯下降。各褶皺產(chǎn)生時的峰值力比之原始結(jié)構(gòu)也有不同程度的下降，整個壓潰過程的力位移曲線更加平穩(wěn)，這是沿管壁分布的通孔導致的。觀察各峰值力的出現(xiàn)位置可以發(fā)現(xiàn)，拓撲結(jié)構(gòu)和工程結(jié)構(gòu)的各峰值力出現(xiàn)位置較原始結(jié)構(gòu)更靠前，間隔也更??；從力位移曲線的峰值力數(shù)量上看，原始結(jié)構(gòu)有5個峰值，而拓撲結(jié)構(gòu)和工程結(jié)構(gòu)擁有6個峰值，結(jié)合圖9的變形模式，觀察到原始結(jié)構(gòu)產(chǎn)生了5個褶皺，拓撲結(jié)構(gòu)和工程結(jié)構(gòu)產(chǎn)生了6個褶皺，這表明拓撲優(yōu)化后的結(jié)構(gòu)在壓潰過程中表現(xiàn)得更加緊湊。綜上所述，雖然拓撲結(jié)構(gòu)的各峰值力出現(xiàn)明顯下降，但由于拓撲優(yōu)化后的模型在整個壓潰過程中保持著一個較高的平均壓潰力，擁有更緊湊的壓潰模式，其能量吸收量并未出現(xiàn)大幅下降。對比拓撲結(jié)構(gòu)與工程結(jié)構(gòu)，力位移曲線變化微小，但由于改變了通孔的形狀，使得壓潰的觸發(fā)順序有所改變，此外，由于改變了拓撲優(yōu)化后的材料分布，使得總吸能量出現(xiàn)小幅下降，下降幅度為4.6%，比吸能下降了1.6%。總而言之，回歸工程化后的結(jié)構(gòu)的耐撞性能比之拓撲結(jié)構(gòu)出現(xiàn)了微幅下降，但優(yōu)于原始結(jié)構(gòu)。

圖9 壓潰模式對比Fig.9 Comparison of crushing mode

圖10 力位移曲線對比Fig.10 Comparison of force displacement curve

4 結(jié)論

1)提出了多設計域優(yōu)化目標分離的BESO算法，并將其用于TSTD在大變形壓潰下的拓撲優(yōu)化設計，優(yōu)化結(jié)果與原始結(jié)構(gòu)的對比表明了所提算法的有效性。

2)沿管壁引入空洞可以有效地降低薄壁管在軸向壓潰過程中的峰值力。

3)合適的材料空洞分布可以誘導吸能效率更高的變形模式，本文為使用拓撲優(yōu)化實現(xiàn)TSTD管的輕量化設計提供了思路。