基于增材制造各向異性的強度約束拓撲優化

何智成 楊丁丁 姜 潮 伍 毅 江和昕

湖南大學汽車車身先進設計制造國家重點實驗室，長沙，410082

0 引言

目前，隨著航空航天、汽車等高端裝備領域的發展，設計中應用優化手段特別是拓撲優化技術的需求越來越強烈[1]。與此同時，增材制造技術的發展與其所具有的超高制造自由度特點，極大程度上推動了拓撲優化在復雜關鍵零部件設計中的應用[2-4]。然而，增材制造工藝因其沿路徑逐層堆疊材料的物理特性，不可避免地會使制造結構具有各方向不同的力學特性[5-7]，傳統基于各向同性假設的設計方法將不再適用。特別是當面對實際的復雜應用工況環境時，如果設計中缺乏對增材制造關鍵零部件失效強度的有效考量，其結構服役的安全性與可靠性將受到極大挑戰。因此，研究考慮增材制造工藝影響的強度優化設計方法具有重要的實際價值與指導意義。

結構的斷裂、疲勞、塑性變形等失效問題大多由局部應力集中引起，所以為了將優化方法推廣到強度設計中，已經發展了許多以結構應力為約束或目標的拓撲優化方法。YANG等[8]與DUYSINX等[9]在早期的連續體結構應力拓撲優化研究中闡明了其面臨的三個主要難點：一是應力的“奇異”現象；二是應力的局部性導致的約束數量龐大；三是計算非線性程度高。為解決這些問題，ε松弛與qp松弛方法[10-11]、全局應力聚合法(如P聚合、KS聚合函數)[8，12]、分區應力聚合[13]和改進的過濾與平均策略被相繼提出，并取得了顯著的效果。

目前，大多數應力拓撲優化研究都基于固體各向同性材料懲罰法(solid isotropic material with penalization，SIMP)[11-13]與水平集方法(level set method，LSM)[14-15]展開，并且它們的有效性已經在各種應用場景中得以證明。近幾年，雙向漸進結構優化方法(bi-directional evolutionary structural optimization，BESO)也逐漸被應用于解決應力拓撲優化問題。XIA等[16]、FAN等[17]基于BESO方法提出了一系列有效求解應力拓撲優化問題的新框架；王選等[18]改進了BESO方法以有效降低其處理應力約束問題時產生的震蕩現象。值得注意的是，與SIMP法相比，BESO法與LSM法由于沒有中間密度單元的概念，本身就具備避免應力“奇異”現象的潛力，但LSM法由于對初始設計較為依賴，容易受到一定限制[19]。因此基于BESO方法展開應力相關的拓撲優化研究具有很大的潛力。

盡管現有的應力拓撲優化工作能有效處理大部分涉及強度的優化問題，但這些研究幾乎均基于各向同性的材料假設。近年來，針對增材制造工藝下的各向異性拓撲優化問題，已經有部分學者從不同角度開展了一些研究工作。CHIU等[20]研究了剛度各向異性與角度對優化結構最小化體積的影響。YANG等[5]基于光固化(Stereolithography，SLA)增材制造工藝下的材料試驗實現了各向異性的剛度優化設計。YU等[21]在材料各向異性和沉積角度的基礎上，實現了打印沉積路徑規劃與結構優化的并行設計。MIRZENDEHDEL等[22]提出了一種基于LSM的增材制造各向異性強度優化設計方法，并通過實驗驗證了方法的有效性。LI等[23]、WU等[24]針對增材制造結構的各向異性斷裂行為，開發了一種增強結構抗彈塑性斷裂的拓撲優化方法。除此之外，諸如多材料多打印方向優化、各向異性點陣結構優化等[25]也得到了廣泛關注。

然而這些針對增材制造各向異性的研究大部分只考慮了結構的整體剛度性能，而對強度相關問題的必要關注非常有限。因此，本文在BESO方法基礎上提出了一種考慮各向異性強度約束的拓撲優化方法。基于Tsai-Hill失效準則詳細推導了描述各向異性結構偏軸(考慮材料堆疊角度改變)強度性能的Tsai-Hill失效系數，并通過權重因子α將失效系數約束引入拓撲優化模型。詳細推導了其優化靈敏度公式，引入歸一化、敏度與密度過濾等一系列數值方法穩定優化歷程。最后基于兩個典型算例驗證算法的有效性，探討各向異性的影響規律。

1 各向異性失效分析

1.1 Tsai-Hill失效準則

為使優化結果更貼近設計目標，使用各向異性Tsai-Hill失效準則以更準確地評估結構的失效風險。該準則已被證明在描述增材制造結構的各向異性力學性能方面是非常有效的[26-27]，本文只討論平面應力問題，表達式如下：

(1)

其中，F.I.=1表示失效臨界狀態，Fij為材料強度系數，σ1、σ2分別為材料坐標系下1方向與2方向的主應力分量，τ12為切應力。材料坐標系O′12與自然坐標系Oxy如圖1所示。F11、F22、F12、F66可通過下式計算：

(2)

圖1 正交各向異性材料坐標系Fig.1 Orthotropic material coordinate systems

圖2展示了平面應力問題下Tsai-Hill準則的失效空間(圖2a)與失效面(圖2b)。失效空間是由滿足準則的強度極限所圍成的封閉區域，為保證優化結構不失效，需使所有單元應力被約束其中。失效面為圖2b所示的橢圓，是失效空間在τ12=0平面(即主應力平面)的投影，面內包絡單元主應力。失效面與σ1軸交點為強度極限X1，與σ2軸交點為X2。通過失效面與單元應力的分布可簡單直觀地描述結構的應力狀態和失效風險。

(a)Tsai-Hill失效空間 (b)Tsai-Hill失效面圖2 Tsai-Hill失效空間與失效面Fig.2 Tsai-Hill failure space and failure surface

1.2 Tsai-Hill準則的坐標變換

在增材制造工藝下，材料堆疊方向會隨著結構的不同而做出相應調整，使得材料坐標系與自然坐標系之間產生一定的偏轉角θ，如圖1所示。因此增材制造結構沿xOy坐標軸方向的強度極限(即偏軸強度極限)會隨著堆疊角θ的偏轉而發生改變，此時式(1)將不再適用，因此需進行坐標變換以獲取考慮堆疊角的Tsai-Hill準則表達式。令

m=cosθ，n=sinθ

(3)

則σx、σy、τxy與σ1、σ2、τ12之間的關系為

(4)

式中，σx、σy、τxy為自然坐標系Oxy下的應力分量。

將式(4)代入式(1)求得Tsai-Hill準則在Oxy坐標系下的表達式：

2Fxsσxτxy+2Fysσyτxy=1

(5)

其中，Fxx、Fyy、Fxy、Fss、Fxs、Fys(下標s表示剪切方向)為自然坐標系下的強度系數，可由下式計算：

(6)

式(5)可重寫為

(7)

(8)

式中，ITH為Oxy坐標系下的強度系數矩陣；σe為單元應力狀態向量。

1.3 失效系數計算

將式(7)改寫為如下方程：

(9)

其中，sTH定義為Tsai-Hill準則的安全系數[28]，sTH≥1表示結構不失效。為了更好地評估結構的失效風險，定義失效系數?TH滿足

(10)

分析式(10)可知，?TH≤1代表結構不失效。

2 各向異性強度約束拓撲優化模型

2.1 優化模型

針對考慮各向異性的強度拓撲優化問題，需保證結構中最大Tsai-Hill失效系數滿足約束，從而得到比傳統柔度最小化(剛度最大化)更加安全的設計結果。其拓撲優化數學模型可表示為

(11)

當優化問題涉及除體積約束以外的其他約束時，一般通過構造一個包含原目標函數與其他約束的新目標函數，將有約束問題轉化為無約束問題。本文構建的新優化目標函數如下：

(12)

式中，α為權重因子，α的大小控制著柔度C和失效系數約束在目標函數η中所占的比例。

2.2 材料插值方案

在傳統BESO框架中，設計域被離散成N個設計變量，設計變量的存在與否由1和xmin表示，且不存在中間狀態。為了避免低密度單元可能引起的應力奇異問題，采用以下插值策略：

De=xeD0

(13)

式中，D0為材料彈性矩陣；De為單元彈性矩陣。

因此，由式(13)可知，只需考慮實體單元的應力狀態計算，刪除的空單元應力直接置零，表達式為

(14)

式中，Be為單元應變矩陣；ue為單元位移向量。

2.3 全局失效系數度量

(15)

需要注意的是，P值趨近無限大時會使優化發生強烈震蕩，因此P值需要適當選取[13]。替換后目標函數式(12)可改進為

(16)

3 靈敏度計算方法

(17)

3.1 柔度靈敏度

(18)

式中，Re為變換矩陣，滿足關系式ue=ReU；ke為單元位移向量；k0為實體單元剛度矩陣。

3.2 失效系數約束靈敏度

(19)

(20)

(21)

式(19)中的偏導?σe/?xe項滿足式(14)插值模型，因此可以按照鏈式法則展開推導，得

(22)

將式(20)～式(22)代入式(19)中得

(23)

(24)

(25)

式中，R為伴隨載荷。

4 優化求解策略

4.1 權重因子的更新

(26)

式中，q為區間控制參數，取值范圍推薦為[0.9，1]。

(27)

4.2 數值處理方法

與式(13)插值策略類似，對其單元靈敏度進行以下處理：

(28)

當xe=xmin時，低密度單元的靈敏度將趨于0，有效杜絕了應力奇異問題。另外有研究表明，BESO方法在處理非線性問題時，結構的進化方向容易受到前步中值較大的“主導靈敏度”影響而產生震蕩[29]，因此需對靈敏度進行歸一化處理，即

(29)

為了保證相鄰單元間靈敏度的連續性以避免棋盤格現象，并確保被刪除單元有機會被再次增添，需對單元靈敏度進行過濾，實現敏度再分配：

(30)

為了優化的穩定性，當前靈敏度的最終值也需考慮歷史平均技術：

(31)

(32)

最后對設計變量進行密度過濾。在應力等非線性優化中，密度過濾可以獲得更好的結構[13]：

(33)

4.3 優化算法流程

所提各向異性強度約束的拓撲優化算法求解流程可歸納為如下步驟：

(1)定義有限元網格，初始化設計域和邊界條件。

(3)通過V(l)=max(V(l-1)(1-γ)，χ)更新當前迭代步體積分數，使結構體積逐漸達到目標體積分數χ并保持穩定。

a.依據式(28)計算單元初始靈敏度；

b.依據式(29)～式(31)計算修正后的單元靈敏度值；

c.根據式(32)與式(33)更新設計變量；

(6)重復執行步驟(3)～(5)，直至滿足下式的收斂條件：

(34)

或達到目標體積后迭代一定步數并經驗性停止，輸出最優解。

5 數值算例

本節采用兩個經典數值算例來驗證所提方法的有效性。數值算例中材料的彈性模量為210 MPa，泊松比為0.3；材料1方向的屈服強度為X1，材料2方向的屈服強度為X2=ψX1，其中ψ為兩正交方向的強度比。靈敏度過濾半徑rs=2，密度過濾半徑rd=3.5，目標體積分數為0.5，體積進化率為0.02，最大體積增添率[16-17]為0.005，P范數隨優化過程從6逐漸增加到7。

為保證探究各向異性強度比ψ改變對結果影響時的公平性，定義以失效空間體積不變為調整不同強度比ψ下材料方向強度X1、X2的依據。圖3a所示為X1=X2=100 MPa即ψ=1時所對應失效空間，圖3b所示為X1≠X2即ψ≠1時的失效空間，已知ψ值，可通過聯立空間體積相等方程計算確定相應的X1、X2值。

(a)ψ=1 (b)ψ=0.65圖3 不同ψ值下的Tsai-Hill失效空間Fig.3 Tsai-Hill failure space for different values of ψ

5.1 算例1

本算例選取經典算例L形梁進行算法的驗證研究，其邊界條件設置如圖4所示，結構劃分為10 000個四節點平面應力單元，上端施加固定約束。為避免施力點的應力集中，載荷F=30 N均布在右側5個節點上，自然坐標系標注在左下角。圖5中綠色坐標系為材料坐標系，坐標軸的方向和長度分別代表材料1、2兩主方向與對應方向上材料的強度極限X1、X2，右上角繪制相應優化結果的主應力失效面與約束狀態，云圖中紅色代表高失效系數，藍色代表低失效系數。

圖4 L形梁問題的邊界條件Fig.4 Boundary conditions for the L-beam problem

表1 不同失效系數約束下的設計性能對比

圖時的優化迭代歷程Fig.6 Optimization iteration process of

(a)無約束剛度設計 (b)von-Mesis應力設計

表2 不同優化方法設計結果性能對比

為了研究不同強度比ψ對優化結果的影響，表3左三列分別展示了ψ為0.65、1、1.5的優化結果。表3所示的失效系數分布云圖中主應力的單位均為MPa。當ψ=1即材料展現各向同性強度特性時，其失效面橢圓(云圖右上角)的長軸與主應力平面45°夾角方向重合，面內可包絡更多剛度最大化設計的主應力點(圖5a右上角)，因此較易滿足約束且優化結構剛度最好。隨著各向異性程度即ψ值的減小(ψ=0.65)或增大(ψ=1.5)，失效面橢圓長軸會向主應力平面橫軸或縱軸靠攏，導致包絡區域偏離剛度最大化設計的主應力點，約束較難滿足且結構剛度均有所下降，較大的失效系數逐漸向著設計結果的縱向或橫向分布。

(a)von-Mesis應力約束剛度設計 (b)Tsai-Hill失效系數約束剛度設計圖8 von-Mesis應力約束與Tsai-Hill失效系數約束的主應力失效面對比Fig.8 Comparison of von-Mesis stress constraint and Tsai-Hill failure coefficient constraint failure surface

為驗證增材制造材料堆疊角度變化對優化結果的影響，表3右三列分別展示了材料坐標系偏轉角θ為0°、45°、90°時的優化結果。可以看出，堆疊角的改變會使結構發生明顯變化，在滿足相同失效系數約束的條件下，不同堆疊角下的結構剛度也會發生改變。由此可見，考慮材料堆疊角的影響對尋找增材制造結構在特定工藝下的最優解非常重要。

5.2 算例2

本算例采用T形梁進行優化驗證，其邊界條件設置如圖9所示，結構劃分為13 750個四節點平面應力單元。左右兩端施加固定約束，載荷F=30 N均布在頂部中間的7個節點上，以避免施力處的應力集中，自然坐標系為圖9中黑色坐標系。圖10中綠色為材料坐標系。

圖9 T形梁問題的邊界條件Fig.9 Boundary conditions for the T-beam problem

表4 不同失效系數約束下的設計性能對比

(a)無約束

表5中右三列展示了不同材料堆疊角對T形梁優化結果的影響。可以看出，不同角度下的優化設計剛度相差較大，且相較于θ=0°與θ=90°時的對稱優化結果，θ=45°時出現了明顯的非對稱結構，這是由于材料低強度方向隨著堆疊角θ的改變而發生偏轉，為滿足設計性能需求，優化結果將生長出更多斜向枝干以彌補結構性能的不足。

6 結論

(1)本文面向增材制造工藝所引起的結構各向異性強度特性，基于BESO方法提出了一種考慮各向異性強度約束的拓撲優化策略。結果表明，傳統剛度優化設計難以避免高失效風險區域，本文方法可以在優化結構剛度的同時有效抑制失效集中以保證結構強度，且優化結果邊界清晰、優化過程迭代平穩。通過比較發現，在增材制造的各向異性強度假設下，本文方法可以獲得比傳統von-Mesis應力相關設計性能更好的結構。

(2)材料的各向異性程度與堆疊角度參數均會對優化產生很大影響。隨著各向異性程度的加劇，優化結構的性能也會隨之惡化，材料會趨于向強度較低的方向分布。與此同時，改變堆疊角度也會影響結構性能，即優化結果很大程度依賴于參數的調整。因此，考慮各向異性的影響有其實際必要性。