基于增材制造各向異性的強(qiáng)度約束拓?fù)鋬?yōu)化

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

湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國家重點(diǎn)實(shí)驗(yàn)室，長沙，410082

0 引言

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

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

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

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

然而這些針對增材制造各向異性的研究大部分只考慮了結(jié)構(gòu)的整體剛度性能，而對強(qiáng)度相關(guān)問題的必要關(guān)注非常有限。因此，本文在BESO方法基礎(chǔ)上提出了一種考慮各向異性強(qiáng)度約束的拓?fù)鋬?yōu)化方法?；赥sai-Hill失效準(zhǔn)則詳細(xì)推導(dǎo)了描述各向異性結(jié)構(gòu)偏軸(考慮材料堆疊角度改變)強(qiáng)度性能的Tsai-Hill失效系數(shù)，并通過權(quán)重因子α將失效系數(shù)約束引入拓?fù)鋬?yōu)化模型。詳細(xì)推導(dǎo)了其優(yōu)化靈敏度公式，引入歸一化、敏度與密度過濾等一系列數(shù)值方法穩(wěn)定優(yōu)化歷程。最后基于兩個(gè)典型算例驗(yàn)證算法的有效性，探討各向異性的影響規(guī)律。

1 各向異性失效分析

1.1 Tsai-Hill失效準(zhǔn)則

為使優(yōu)化結(jié)果更貼近設(shè)計(jì)目標(biāo)，使用各向異性Tsai-Hill失效準(zhǔn)則以更準(zhǔn)確地評估結(jié)構(gòu)的失效風(fēng)險(xiǎn)。該準(zhǔn)則已被證明在描述增材制造結(jié)構(gòu)的各向異性力學(xué)性能方面是非常有效的[26-27]，本文只討論平面應(yīng)力問題，表達(dá)式如下：

(1)

其中，F(xiàn).I.=1表示失效臨界狀態(tài)，F(xiàn)ij為材料強(qiáng)度系數(shù)，σ1、σ2分別為材料坐標(biāo)系下1方向與2方向的主應(yīng)力分量，τ12為切應(yīng)力。材料坐標(biāo)系O′12與自然坐標(biāo)系Oxy如圖1所示。F11、F22、F12、F66可通過下式計(jì)算：

(2)

圖1 正交各向異性材料坐標(biāo)系Fig.1 Orthotropic material coordinate systems

圖2展示了平面應(yīng)力問題下Tsai-Hill準(zhǔn)則的失效空間(圖2a)與失效面(圖2b)。失效空間是由滿足準(zhǔn)則的強(qiáng)度極限所圍成的封閉區(qū)域，為保證優(yōu)化結(jié)構(gòu)不失效，需使所有單元應(yīng)力被約束其中。失效面為圖2b所示的橢圓，是失效空間在τ12=0平面(即主應(yīng)力平面)的投影，面內(nèi)包絡(luò)單元主應(yīng)力。失效面與σ1軸交點(diǎn)為強(qiáng)度極限X1，與σ2軸交點(diǎn)為X2。通過失效面與單元應(yīng)力的分布可簡單直觀地描述結(jié)構(gòu)的應(yīng)力狀態(tài)和失效風(fēng)險(xiǎn)。

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

1.2 Tsai-Hill準(zhǔn)則的坐標(biāo)變換

在增材制造工藝下，材料堆疊方向會(huì)隨著結(jié)構(gòu)的不同而做出相應(yīng)調(diào)整，使得材料坐標(biāo)系與自然坐標(biāo)系之間產(chǎn)生一定的偏轉(zhuǎn)角θ，如圖1所示。因此增材制造結(jié)構(gòu)沿xOy坐標(biāo)軸方向的強(qiáng)度極限(即偏軸強(qiáng)度極限)會(huì)隨著堆疊角θ的偏轉(zhuǎn)而發(fā)生改變，此時(shí)式(1)將不再適用，因此需進(jìn)行坐標(biāo)變換以獲取考慮堆疊角的Tsai-Hill準(zhǔn)則表達(dá)式。令

m=cosθ，n=sinθ

(3)

則σx、σy、τxy與σ1、σ2、τ12之間的關(guān)系為

(4)

式中，σx、σy、τxy為自然坐標(biāo)系Oxy下的應(yīng)力分量。

將式(4)代入式(1)求得Tsai-Hill準(zhǔn)則在Oxy坐標(biāo)系下的表達(dá)式：

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

(5)

其中，F(xiàn)xx、Fyy、Fxy、Fss、Fxs、Fys(下標(biāo)s表示剪切方向)為自然坐標(biāo)系下的強(qiáng)度系數(shù)，可由下式計(jì)算：

(6)

式(5)可重寫為

(7)

(8)

式中，ITH為Oxy坐標(biāo)系下的強(qiáng)度系數(shù)矩陣；σe為單元應(yīng)力狀態(tài)向量。

1.3 失效系數(shù)計(jì)算

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

(9)

其中，sTH定義為Tsai-Hill準(zhǔn)則的安全系數(shù)[28]，sTH≥1表示結(jié)構(gòu)不失效。為了更好地評估結(jié)構(gòu)的失效風(fēng)險(xiǎn)，定義失效系數(shù)?TH滿足

(10)

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

2 各向異性強(qiáng)度約束拓?fù)鋬?yōu)化模型

2.1 優(yōu)化模型

針對考慮各向異性的強(qiáng)度拓?fù)鋬?yōu)化問題，需保證結(jié)構(gòu)中最大Tsai-Hill失效系數(shù)滿足約束，從而得到比傳統(tǒng)柔度最小化(剛度最大化)更加安全的設(shè)計(jì)結(jié)果。其拓?fù)鋬?yōu)化數(shù)學(xué)模型可表示為

(11)

當(dāng)優(yōu)化問題涉及除體積約束以外的其他約束時(shí)，一般通過構(gòu)造一個(gè)包含原目標(biāo)函數(shù)與其他約束的新目標(biāo)函數(shù)，將有約束問題轉(zhuǎn)化為無約束問題。本文構(gòu)建的新優(yōu)化目標(biāo)函數(shù)如下：

(12)

式中，α為權(quán)重因子，α的大小控制著柔度C和失效系數(shù)約束在目標(biāo)函數(shù)η中所占的比例。

2.2 材料插值方案

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

De=xeD0

(13)

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

因此，由式(13)可知，只需考慮實(shí)體單元的應(yīng)力狀態(tài)計(jì)算，刪除的空單元應(yīng)力直接置零，表達(dá)式為

(14)

式中，Be為單元應(yīng)變矩陣；ue為單元位移向量。

2.3 全局失效系數(shù)度量

(15)

需要注意的是，P值趨近無限大時(shí)會(huì)使優(yōu)化發(fā)生強(qiáng)烈震蕩，因此P值需要適當(dāng)選取[13]。替換后目標(biāo)函數(shù)式(12)可改進(jìn)為

(16)

3 靈敏度計(jì)算方法

(17)

3.1 柔度靈敏度

(18)

式中，Re為變換矩陣，滿足關(guān)系式ue=ReU；ke為單元位移向量；k0為實(shí)體單元?jiǎng)偠染仃嚒?/p>

3.2 失效系數(shù)約束靈敏度

(19)

(20)

(21)

式(19)中的偏導(dǎo)?σe/?xe項(xiàng)滿足式(14)插值模型，因此可以按照鏈?zhǔn)椒▌t展開推導(dǎo)，得

(22)

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

(23)

(24)

(25)

式中，R為伴隨載荷。

4 優(yōu)化求解策略

4.1 權(quán)重因子的更新

(26)

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

(27)

4.2 數(shù)值處理方法

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

(28)

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

(29)

為了保證相鄰單元間靈敏度的連續(xù)性以避免棋盤格現(xiàn)象，并確保被刪除單元有機(jī)會(huì)被再次增添，需對單元靈敏度進(jìn)行過濾，實(shí)現(xiàn)敏度再分配：

(30)

為了優(yōu)化的穩(wěn)定性，當(dāng)前靈敏度的最終值也需考慮歷史平均技術(shù)：

(31)

(32)

最后對設(shè)計(jì)變量進(jìn)行密度過濾。在應(yīng)力等非線性優(yōu)化中，密度過濾可以獲得更好的結(jié)構(gòu)[13]：

(33)

4.3 優(yōu)化算法流程

所提各向異性強(qiáng)度約束的拓?fù)鋬?yōu)化算法求解流程可歸納為如下步驟：

(1)定義有限元網(wǎng)格，初始化設(shè)計(jì)域和邊界條件。

(3)通過V(l)=max(V(l-1)(1-γ)，χ)更新當(dāng)前迭代步體積分?jǐn)?shù)，使結(jié)構(gòu)體積逐漸達(dá)到目標(biāo)體積分?jǐn)?shù)χ并保持穩(wěn)定。

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

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

c.根據(jù)式(32)與式(33)更新設(shè)計(jì)變量；

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

(34)

或達(dá)到目標(biāo)體積后迭代一定步數(shù)并經(jīng)驗(yàn)性停止，輸出最優(yōu)解。

5 數(shù)值算例

本節(jié)采用兩個(gè)經(jīng)典數(shù)值算例來驗(yàn)證所提方法的有效性。數(shù)值算例中材料的彈性模量為210 MPa，泊松比為0.3；材料1方向的屈服強(qiáng)度為X1，材料2方向的屈服強(qiáng)度為X2=ψX1，其中ψ為兩正交方向的強(qiáng)度比。靈敏度過濾半徑rs=2，密度過濾半徑rd=3.5，目標(biāo)體積分?jǐn)?shù)為0.5，體積進(jìn)化率為0.02，最大體積增添率[16-17]為0.005，P范數(shù)隨優(yōu)化過程從6逐漸增加到7。

為保證探究各向異性強(qiáng)度比ψ改變對結(jié)果影響時(shí)的公平性，定義以失效空間體積不變?yōu)檎{(diào)整不同強(qiáng)度比ψ下材料方向強(qiáng)度X1、X2的依據(jù)。圖3a所示為X1=X2=100 MPa即ψ=1時(shí)所對應(yīng)失效空間，圖3b所示為X1≠X2即ψ≠1時(shí)的失效空間，已知ψ值，可通過聯(lián)立空間體積相等方程計(jì)算確定相應(yīng)的X1、X2值。

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

5.1 算例1

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

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

表1 不同失效系數(shù)約束下的設(shè)計(jì)性能對比

圖時(shí)的優(yōu)化迭代歷程Fig.6 Optimization iteration process of

(a)無約束剛度設(shè)計(jì) (b)von-Mesis應(yīng)力設(shè)計(jì)

表2 不同優(yōu)化方法設(shè)計(jì)結(jié)果性能對比

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

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

為驗(yàn)證增材制造材料堆疊角度變化對優(yōu)化結(jié)果的影響，表3右三列分別展示了材料坐標(biāo)系偏轉(zhuǎn)角θ為0°、45°、90°時(shí)的優(yōu)化結(jié)果?？梢钥闯觯询B角的改變會(huì)使結(jié)構(gòu)發(fā)生明顯變化，在滿足相同失效系數(shù)約束的條件下，不同堆疊角下的結(jié)構(gòu)剛度也會(huì)發(fā)生改變。由此可見，考慮材料堆疊角的影響對尋找增材制造結(jié)構(gòu)在特定工藝下的最優(yōu)解非常重要。

5.2 算例2

本算例采用T形梁進(jìn)行優(yōu)化驗(yàn)證，其邊界條件設(shè)置如圖9所示，結(jié)構(gòu)劃分為13 750個(gè)四節(jié)點(diǎn)平面應(yīng)力單元。左右兩端施加固定約束，載荷F=30 N均布在頂部中間的7個(gè)節(jié)點(diǎn)上，以避免施力處的應(yīng)力集中，自然坐標(biāo)系為圖9中黑色坐標(biāo)系。圖10中綠色為材料坐標(biāo)系。

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

表4 不同失效系數(shù)約束下的設(shè)計(jì)性能對比

(a)無約束

表5中右三列展示了不同材料堆疊角對T形梁優(yōu)化結(jié)果的影響。可以看出，不同角度下的優(yōu)化設(shè)計(jì)剛度相差較大，且相較于θ=0°與θ=90°時(shí)的對稱優(yōu)化結(jié)果，θ=45°時(shí)出現(xiàn)了明顯的非對稱結(jié)構(gòu)，這是由于材料低強(qiáng)度方向隨著堆疊角θ的改變而發(fā)生偏轉(zhuǎn)，為滿足設(shè)計(jì)性能需求，優(yōu)化結(jié)果將生長出更多斜向枝干以彌補(bǔ)結(jié)構(gòu)性能的不足。

6 結(jié)論

(1)本文面向增材制造工藝所引起的結(jié)構(gòu)各向異性強(qiáng)度特性，基于BESO方法提出了一種考慮各向異性強(qiáng)度約束的拓?fù)鋬?yōu)化策略。結(jié)果表明，傳統(tǒng)剛度優(yōu)化設(shè)計(jì)難以避免高失效風(fēng)險(xiǎn)區(qū)域，本文方法可以在優(yōu)化結(jié)構(gòu)剛度的同時(shí)有效抑制失效集中以保證結(jié)構(gòu)強(qiáng)度，且優(yōu)化結(jié)果邊界清晰、優(yōu)化過程迭代平穩(wěn)。通過比較發(fā)現(xiàn)，在增材制造的各向異性強(qiáng)度假設(shè)下，本文方法可以獲得比傳統(tǒng)von-Mesis應(yīng)力相關(guān)設(shè)計(jì)性能更好的結(jié)構(gòu)。

(2)材料的各向異性程度與堆疊角度參數(shù)均會(huì)對優(yōu)化產(chǎn)生很大影響。隨著各向異性程度的加劇，優(yōu)化結(jié)構(gòu)的性能也會(huì)隨之惡化，材料會(huì)趨于向強(qiáng)度較低的方向分布。與此同時(shí)，改變堆疊角度也會(huì)影響結(jié)構(gòu)性能，即優(yōu)化結(jié)果很大程度依賴于參數(shù)的調(diào)整。因此，考慮各向異性的影響有其實(shí)際必要性。