局部有限壽命疲勞約束條件下的結構拓撲優化方法

江旭東, 武子旺, 滕曉艷

(1. 哈爾濱理工大學 機械動力工程學院,哈爾濱 150080;2. 哈爾濱工程大學 機電工程學院,哈爾濱 150001)

面向金屬增材制造的結構拓撲優化技術極大地豐富了工程結構的設計空間,兼顧了復雜結構和高性能構件成形需求,在航空、航天、交通、核電等領域具有廣闊的應用前景和發展空間[1-2]。目前,多數的拓撲優化研究工作專注于體積約束下的結構剛度最大化問題,而工程結構往往服役于交變載荷作用的工作環境,剛度最優的設計一般不能完全滿足抗疲勞要求,從而嚴重影響結構在全壽命周期內運行的可靠性。為此,在結構的輕量化設計中考慮疲勞性能約束,對豐富工程結構強度設計理論、提高服役性能有著重要的科學價值和工程意義。

為了抑制應力集中和高應力值引起的結構斷裂,應力約束下的拓撲優化問題引起了高度的關注。應力約束拓撲優化設計的難點在于大量局部約束增加了敏度分析的計算代價。為了減少應力約束數量,一般采用P范數[3-5]和K-S函數[6]等凝聚方法將局部應力約束轉化為一個全局約束。為了保持應力約束的局部性本質,代替上述函數凝聚技術,Silva等[7-9]利用增廣拉格朗日方法,將多約束問題轉化為一系列無約束子問題,求解了局部應力約束的拓撲優化問題。最近,Oliver等[10]將非結構化多邊形網格技術[11]與增廣拉格朗日方法相融合,提出了適用于復雜邊界工程結構的局部應力約束拓撲優化設計方法。與應力約束拓撲優化問題不同,疲勞失效與加載歷史具有強相關性,局部特征更為顯著,導致疲勞約束拓撲優化問題呈現高度非線性,因而相關研究工作相對較少。研究學者利用Palmgren-Miner[12-14]、修正的Goodman[15-16]、Sines[17]及考慮初始缺陷的Murakami[18]等疲勞失效準則,結合準靜態分析法[19-20]、頻域[21-23]或時域等效靜載荷等動態分析法[24-25]預測結構疲勞損傷,通過P范數方法凝聚疲勞約束,初步提出了考慮疲勞性能的連續體結構拓撲優化設計方法。盡管P范數方法能夠大幅減少疲勞約束及其敏度分析的計算成本,但是其近似屬性并不能精確地控制疲勞約束。

由此,本文基于增廣拉格朗日方法提出局部疲勞約束問題的拓撲優化框架。考慮變幅比例載荷作用,通過Palmgren-Miner疲勞準則評定結構的疲勞損傷,以材料用量為目標函數,以疲勞性能為約束條件,提出基于局部疲勞約束的結構輕量化設計方法。此外,將非結構化多邊形網格技術融入到拓撲優化模型,實現復雜幾何邊界結構的抗疲勞拓撲優化設計。對比現有P范數方法的優化結果,驗證提出的非凝聚方法處理疲勞約束拓撲優化問題的有效性。

1 優化模型

在本文提出的拓撲優化框架內,以結構輕量化為目標,分別求解非凝聚的局部疲勞約束問題(Q1)和P范數凝聚的全局疲勞約束問題(Q2),對比分析兩個問題的優化解,驗證提出方法的有效性。

(1)

1.1 疲勞失效分析

為了估計變幅比例載荷作用下結構的高周疲勞損傷,本文通過雨流計數法確定多軸應力狀態下的平均應力和應力幅,采用Sines準則評估疲勞等效應力,最后基于Palmgren-Miner線性累積損傷模型評估結構的疲勞失效。工程結構高周疲勞一般表現為低應力失效,因而可采用線性準靜態分析法,通過線性疊加參考載荷下的結構響應,獲得變幅載荷作用下的結構疲勞響應。

基于修正的固體各向同性材料懲罰模型(solid isotropic material with penalization,SIMP),材料插值函數mE[y(z)]表示為[26]

mE[y(z)]=ε+(1-ε)mV[y(z)]p

(2)

(3)

式中:Ersatz數ε?1;p為懲罰因子;mV[y(z)]為Heaviside投影體積插值函數;α為門檻密度;η為投影控制參數。

為了加強優化問題的適定性,在規范化映射下將設計變量函數與線性“帽子”核函數卷積[27],使單元密度變量與設計變量具有如下關系

y=Pz

(4)

(5)

(6)

式中:P為過濾矩陣;R為過濾半徑; ‖x-xk‖2為單元與k形心位置矢量x,xk的歐式距離;q為濾波指數。

為了能夠求解具有復雜設計域結構的拓撲優化問題,利用非結構化多邊形網格技術離散空間設計域,如圖1所示,對于n邊形參考單元,其節點i處的Wachspress形函數定義為

圖1 多邊形單元Fig.1 Polygon element definition

(7)

(8)

式中:αi(ξ)為插值函數;ξ=[ξ1ξ2]為n邊形單元的內點;Ai(ξ),Ai+1(ξ)為相應陰影三角形的面積(如圖1(a)所示),可根據式(9)計算獲得

(9)

式中,pi-1=[p1,i-1p2,i-1]和pi=[p1,ip2,i]為n邊形單元的相鄰節點位置矢量。

考慮如圖2所示的變幅載荷作用,采用準靜態分析法求解結構在參考載荷下的疲勞響應,則有

圖2 雨流計數法Fig.2 Rainflow-counting method

K[y(z)]u=Fr

(10)

(11)

式中:K[y(z)]為結構插值總體剛度矩陣;Ke0為密度變量為1時的單元剛度矩陣。

為了避免優化問題的奇異解現象,采用修正的qp松弛技術處理單元應力。對于無預應力的線彈性結構,參考載荷作用下的單元參考應力表示為

(12)

為了預測工程結構在變幅載荷下的疲勞損傷,往往通過雨流計數法從應力譜中提取峰值與谷值(見圖2),確定不同循環特征的應力循環,由此獲得應力幅縮放因子與平均應力縮放因子,根據式(13)估算任一應力循環的應力狀態。

(13)

式中:σea,i,σem,i分別為單元e在第i個循環的應力幅和平均應力列陣;cai,cmi分別為第i個循環的應力幅和平均應力縮放因子。

根據Sines有限壽命疲勞準則,平面應力狀態下有限壽命的等效單軸應力可由交變的八面體剪切應力和靜水平均應力依據式(14)獲得[28]

(14)

Basquin方程描述了單軸等效應力與疲勞壽命的關系,表示為

(15)

式中:σf為疲勞強度系數;Nei為應力循環i下結構的疲勞壽命;b為疲勞強度指數。

由此,根據式(14)、式(15),采用Palmgren-Miner線性累積損傷模型可以評定整個載荷譜作用下的任一單元的疲勞損傷,工程結構須滿足如下疲勞約束,方能避免疲勞斷裂發生。

(16)

式中:De為單元e的累計損傷;cD為縮放參數(一般大于1);ni為載荷循環i的循環次數;nRF為不同應力循環特征的組合總數。

1.2 疲勞約束

對于優化問題Q2, P范數疲勞約束可表示為

(17)

式中,P為P范數因子。

高的P范數因子有利于估計函數的最大值,但是,也會增加優化問題的非線性,不利于問題求解和造成收斂困難。為了使P范數損傷靠近實際損傷的最大值,一般采用自適應約束縮放技術[29]修正P范數損傷。

(18)

如圖3所示,多項式約束函數式(18)在約束違背域的懲罰幅度遠大于線性約束函數,而在Λj≈0處類似于傳統線性約束。盡管式(18)的立方項增加了疲勞約束的非線性,但是,數值試驗表明:相對于線性約束函數,對疲勞約束違背域施加非線性懲罰能夠驅動優化解更加快速地向低損傷方向迭代。

圖3 傳統線性約束與三次方約束的比較Fig.3 Comparison between traditional linear and cubic constraints

1.3 敏度分析

依據鏈式求導法則,優化問題Q1的目標函數對于設計變量的靈敏度為

(19)

根據優化問題式(1),目標函數對于材料插值函數和體積插值函數的偏導數為

(20)

給定材料剛度參數E=mE(y)和體積分數V=mV(y),它們對于設計變量的偏導數為

(21)

多項式疲勞性能約束對于設計變量的靈敏度表示為

(22)

由式(18),多項式疲勞性能約束對于材料插值函數和體積插值函數的偏導數為

(23)

(24)

(25)

根據Palmgren-Miner線性累積損傷模型式(16)和式(15),損傷變量與疲勞壽命的偏導數分別為

(26)

將式(12)代入式(13),可獲得應力幅列陣和平均應力列陣的偏導數,表示為

(27)

(28)

(29)

2 增廣拉格朗日方法

圖4 增廣拉格朗日優化框架流程圖Fig.4 Schematic flowchart of the AL-based topology optimization framework

當單元的數量足夠大時,懲罰項P(k)(z,u)可能遠大于目標函數f(z),從而引起算法失效。由此,引入歸一化的懲罰項[31]避免優化問題的奇異性,則歸一化的拉格朗日子問題表示為

(30)

(31)

(32)

μ(k+1)=min[αμ(k),μmax]

(33)

式中:α>1為懲罰因子更新參數,上確界μmax用于抑制算法的數值失穩問題。拉格朗日乘子按照式(34)更新

(34)

基于梯度的全局收斂移動漸近線算法[32](globally convergent method of moving asymptotes, GCMMA)具有較好的穩定性和魯棒性,本文采用GCMMA求解上述疲勞約束拓撲優化問題。按照鏈式求導法則,歸一化的拉格朗日目標函數的靈敏度為

(35)

由式(31),P(k)與體積分數V不直接相關,而通過局部疲勞性能約束與剛度參數E相關,它對于上述參數的偏導數為

(36)

(37)

對于變幅載荷作用的多工況問題,伴隨方法可有效求解優化問題式(30)的靈敏度。將參考載荷(分解為應力幅載荷與平均應力載荷)作用下的平衡方程以殘差形式引入到式(37)中,形成增廣形式的懲罰項敏度分析方程,則有

(38)

Rai=Kuai-Fr,ai=0,Rmi=Kumi-Fr,mi=0

(39)

式中:Rai,Rmi分別為應力幅和平均應力平衡方程的殘差表達式;ξai,ξmi為相應的伴隨矢量。

將式(39)對材料剛度參數求偏導數,則有

(40)

觀察式(38),隨著應力循環數的增加,伴隨矢量的計算將消耗大量的時間。因而,通過式(41) 縮放參考載荷作用下的位移,表示任意循環載荷作用下的幅值位移與平均位移,可以有效縮減伴隨矢量的解算規模。

uai=caiu,umi=cmiu

(41)

由此,式(38)的伴隨項表示為

(42)

另外,為了避免直接求解?u/?E,通過合理選擇伴隨矢量滿足式(43),使其在式(42)中消失,則有

(43)

由此,確定所有縮放的伴隨矢量后,結合式(44)則可獲得懲罰項對于剛度參數的偏導數。

(44)

(45)

(46)

最后,由式(12)和式(13),?σea,i/?u和?σem,i/?u可顯式獲得,表示為

(47)

綜上,編制有限壽命疲勞分析的多邊形有限元計算程序、敏度分析程序、設計變量濾波程序與GCMMA求解器,實現局部有限壽命約束條件下的結構優化。

3 結果與討論

3.1 L型板

如圖5所示,L型支架的尺寸參數L=1 m,彈性模量E0=70 GPa,泊松比μ=0.3,疲勞強度系數σf=1 350 MPa,疲勞強度指數b=-0.086。L型支架上端施加固定約束,右側角點施加橫向變幅載荷F(見圖2),其參考載荷Fr=20 kN。過濾半徑R=0.065,懲罰因子p=3,門檻密度α=0.5,投影控制參數η=2～10(增量為0.7,更新頻率為5)。

圖5 L型支架設計域與邊界條件Fig.5 Design domain and boundary conditions for a L-bracket

設計域離散為N=50 176個單元,橫向載荷施加于垂向棱邊的8個單元上,避免在加載區域的應力集中。上述單元不包含在設計域中,它們的密度變量設置為1。施加隨機載荷Pk=rand(Fr,-Fr),其中k=1 000,載荷譜縮放參數cD=750。圖6為L型支架采用P范數凝聚方法(P范數因子P=8)和非凝聚方法獲得的優化拓撲及損傷分布。對比兩種方法的分析結果,盡管設計構型存在一定的相似性,但是,局部約束作用的優化拓撲在載荷施加區域附近具有更加合理的結構形式,結構拓撲接近于等強度設計,因而材料的抗疲勞性能得到充分利用。另外,從材料用量來看,P范數凝聚方法的目標值為0.305,非凝聚化方法的目標值為0.252,后者質量減輕17.4%,驗證了所提方法的有效性。

圖6 L型支架的優化拓撲與損傷分布Fig.6 Optimal topology and damage distribution for a L-bracket

3.2 天線支架

如圖7所示,天線支架的尺寸參數L=4 m,H=6 m,彈性模量E0=120 GPa,泊松比μ=0.3,疲勞強度系數σf=1 350 MPa,疲勞強度指數b=-0.086。天線支架左端施加固定約束,右側角點施加橫向變幅載荷F(見圖2),其參考載荷Fr=400 kN。過濾半徑R=0.06,懲罰因子p=3,門檻密度α=0.5,投影控制參數η=2～10(增量為1.0,更新頻率為5)。

圖7 天線支架設計域與邊界條件Fig.7 Design domain and boundary conditions for an antenna bracket

設計域離散為N=50 000個單元,橫向載荷施加于斜向棱邊的8個單元上,避免在加載區域的應力集中。上述單元不包含在設計域中,它們的密度變量設置為1。施加隨機載荷Pk=rand(Fr,-Fr),其中k=1 000,載荷譜縮放參數cD=30 000。圖8為天線支架采用P范數凝聚方法(P范數因子P=10)和非凝聚方法獲得的優化拓撲及損傷分布。對比兩種方法的分析結果,局部約束作用的優化拓撲改變了內部加強筋的布局,使得在P范數約束優化拓撲中高損傷區域處的疲勞損傷有所下降,導致更少的單元承受疲勞失效,因而具有更優的抗疲勞性能。另外,從材料用量來看,P范數凝聚方法的目標值為0.382,非凝聚化方法的目標值為0.256,后者質量減輕33.0%,有利于工程結構的輕量化設計。

圖8 天線支架的優化拓撲與損傷分布Fig.8 Optimal topology and damage distribution for an antenna bracket

3.3 門式鋼架

如圖9所示,門式鋼架的尺寸參數L=1.6 m,H=0.8 m,l=0.15 m,彈性模量E0=100 GPa,泊松比μ=0.3,疲勞強度系數σf=1 350 MPa,疲勞強度指數b=-0.086。門式鋼架底端施加固定約束,上表面中心處施加橫向變幅載荷F(見圖2),其參考載荷Fr=70 kN。過濾半徑R=0.035,懲罰因子p=4.5,門檻密度α=0.5,投影控制參數η=1～10(增量為3,更新頻率為5)。

圖9 門式鋼架設計域與邊界條件Fig.9 Design domain and boundary conditions for a portal frame

設計域離散為N=50 000個單元,橫向載荷施加于上表面棱線中間的12個單元上,避免在加載區域的應力集中。上述單元不包含在設計域中,它們的密度變量設置為1。施加隨機載荷Pk=rand(Fr,-Fr),其中k=1 000,載荷譜縮放參數cD=10 000。圖10為門式鋼架采用P范數凝聚方法(P范數因子P=10)和非凝聚方法獲得的優化拓撲及損傷分布。對比兩種方法的分析結果,局部約束作用的優化拓撲優化了載荷施加區域附近的傳遞路徑,祛除了對于削弱疲勞損傷作用較小的冗余結構,最優拓撲圖像結構清晰,傳力路徑明確,具有更優的抗疲勞性能。另外,從材料用量來看,P范數凝聚方法的目標值為0.286,非凝聚化方法的目標值為0.226,后者質量減輕21.0%,驗證了所提方法的有效性。

圖10 門式鋼架的優化拓撲與損傷分布Fig.10 Optimal topology and damage distribution for a portal frame

3.4 圓孔矩形板

如圖11所示,圓孔矩形板的尺寸參數L=1.2 m,H=0.4 m,D=0.24 m,彈性模量E0=100 GPa,泊松比μ=0.3,疲勞強度系數σf=1 350 MPa,疲勞強度指數b=-0.086。孔矩形板的左端施加固定約束,右端中心處施加橫向變幅載荷F(見圖2),其參考載荷Fr=15 kN。過濾半徑R=0.05,懲罰因子p=4.5,門檻密度α=0.5,投影控制參數η=2～10(增量為0.7,更新頻率為5)。

圖11 圓孔矩形板設計域與邊界條件Fig.11 Design domain and boundary conditions for a rectangular plate with circular holes

設計域離散為N=30 000個單元,橫向載荷施加于右端垂向棱邊的6個單元上,避免在加載區域的應力集中。上述單元不包含在設計域中,它們的密度變量設置為1。施加隨機載荷Pk=rand(Fr,-Fr),其中k=1 000,載荷譜縮放參數cD=10 000。圖12為圓孔矩形板采用P范數凝聚方法(P范數因子P=8)和非凝聚方法獲得的優化拓撲及損傷分布。對比兩種方法的分析結果,盡管設計構型存在一定的相似性,但是,局部約束作用的優化拓撲在傳力路徑上有效減少了冗余材料,使其與P范數凝聚的優化拓撲具有相近的疲勞損傷,因而具有更優的材料分布。另外,從材料用量來看,P范數凝聚方法的目標值為0.416,非凝聚化方法的目標值為0.386,后者質量減輕7.2%。與實體域結構相比,由于結構中增加的圓孔削弱了抗疲勞能力,導致可減少的材料用量受到限制。

圖12 圓孔矩形板的優化拓撲與損傷分布Fig.12 Optimal topology and damage distribution for a rectangular plate with circular holes

3.5 圓孔曲線板

如圖13所示,圓孔曲線輪廓板的尺寸參數L=2.8 m,D1=0.35 m,D2=0.6 m,R1=0.3 m,R2=0.5 m,彈性模量E0=100 GPa,泊松比μ=0.3,疲勞強度系數σf=1 350 MPa,疲勞強度指數b=-0.086。圓孔曲線輪廓板右端施加固定約束,左端圓孔施加均布橫向變幅載荷F(見圖2),其參考載荷Fr=10 kN。過濾半徑R=0.03,懲罰因子p=3,門檻密度α=0.5,投影控制參數η=1～10(增量為3,更新頻率為5)。

設計域離散為N=30 000個單元,橫向均布載荷施加于左側圓孔下方所有單元上。上述單元不包含在設計域中,它們的密度變量設置為1。施加隨機載荷Pk=rand(Fr,-Fr),其中k=1 000,載荷譜縮放參數cD=11 000。圖14為圓孔曲線輪廓板采用P范數凝聚方法(P范數因子P=10)和非凝聚方法獲得的優化拓撲及損傷分布。對比兩種方法的分析結果,考慮局部約束作用的優化拓撲在兩圓孔間增加了多條加強筋,有效地阻礙了疲勞損傷的傳遞,使得P范數凝聚方法的高損傷區域轉變為中度損傷區域,提高了結構的抗疲勞性能。另外,從材料用量來看,P范數凝聚方法的目標值為0.258,非凝聚化方法的目標值為0.228,后者質量減輕11.6%,驗證了所提方法的優越性。

圖14 圓孔曲線輪廓板的優化拓撲與損傷分布Fig.14 Optimal topology and damage distribution for a curved plate with circular holes

4 結 論

(1) 針對變幅載荷作用下的高周疲勞結構的輕量化問題,利用Palmgren-Miner損傷準則在有限單元上逐一構建多項式局部疲勞性能約束,避免任意材料評估點處的疲勞失效;同時,基于增廣拉格朗日方法將原問題轉化為一系列無約束優化子問題,有效地解決了拓撲優化問題中的大量局部約束所帶來的挑戰。

(2) 與疲勞約束凝聚化的優化結果相比,局部疲勞約束條件所獲得的設計構型與其具有很大的差異,優化結構的疲勞損傷和材料用量都顯著減小,因而考慮局部疲勞約束有益于材料的充分利用和結構抗疲勞性能的改善。

(3) 與現有的基于單元的疲勞問題拓撲優化方法相比,本文采用非結構化的多邊形有限元法完成連續體結構的疲勞分析與損傷評估,能夠實現具有任意曲線邊界設計域結構的拓撲優化設計。

綜上,所提方法為實現考慮疲勞性能約束的連續體結構的輕量化設計提供了一種新的思路。未來工作將以本文為基礎,考慮慣性和阻尼效應的影響,研究動態載荷作用下結構疲勞約束問題的設計優化。