連續纖維增強復合材料結構基頻最大化設計1)

程長征 卞光耀 王選,2) 龍凱 李景傳 吳喬國

*(合肥工業大學工程力學系,合肥 230009)

?(合肥工業大學土木工程結構與材料安徽省重點實驗室,合肥 230009)

**(華北電力大學新能源電力系統國家重點實驗室,北京 102206)

引言

隨著現代工業的發展,很多新型材料被研制并成功應用于不同領域,其中纖維增強復合材料就是最常用的一種,已在汽車、建筑以及航空航天工業等領域得到廣泛應用[1].纖維增強復合材料由基體材料和增強材料兩部分組成.與單獨組分材料和傳統的金屬材料相比,纖維增強復合材料在強度、剛度、抗斷裂、抗疲勞、熱穩定性等諸多方面具備更優良的性能.另一方面,當結構所受的外激勵頻率與自身的固有頻率接近時,就會發生共振現象,共振會導致機器破壞、橋梁倒塌等嚴重危害.因此,本文開展連續纖維增強復合材料結構基頻最大化設計研究具有重要意義.

基頻最大化是結構動力學拓撲優化中被廣泛關注的優化問題之一.Ma 等[2]采用均勻化方法較早研究了動力學拓撲優化問題.Pedersen[3]通過弓入新的材料剛度插值格式避免了基于SIMP (solid isotropic material with penalization) 方法的頻率優化中在低密度區出現的局部模態現象.為了適應實際制造工藝要求,Niu 等[4]采用均勻化方法實現了蜂窩材料結構在宏觀和微觀兩個尺度上同步優化以最大化結構基頻,Huang 等[5]基于BESO(Bi-directional evolutionary structural optimization)方法研究了體積約束下的結構一階頻率最大化問題.Xia 等[6]基于水平集方法研究了結構在無阻尼自由振動下的頻率最大化問題,得到了光滑的拓撲構型邊界,在此基礎上,他們進一步考慮了均勻邊界腐蝕對動力學優化問題中基頻的影響[7].高興軍和馬海濤[8]提出了一種將移頻與虛假模態識別相結合的通用方法,以消去連續體結構動力拓撲優化中局部模態的不利影響.動力學拓撲優化方面的研究還包括最大化高階本征頻率或最大化兩階本征頻率的帶寬[9]、最大化受迫振動下的動柔度[10-13]等.此外,還可以將上面度量結構的振動效應的指標作為約束函數來建模[14-15].文獻[16]從不同優化方法角度對結構動力學拓撲優化問題做了較為全面的綜述.

纖維增強復合材料結構的纖維角度優化問題也是被大家關注的熱點問題之一.Pedersen[17]基于OC準則方法(optimality criteria based method) 實現了層壓膜結構的厚度和纖維角度優化.Luo 和Gea[18]基于能量的方法以應變能最小化或基頻最大化為目標,實現了三維板殼結構的纖維角度優化問題.Setoodeh等[19]基于SIMP 方法研究了單工況和多工況載荷作用下二維復合結構的拓撲構型和纖維角度的同步優化問題,以最大化結構的剛度.Stegmann 和Lund[20]基于多相材料插值思想提出了一種新穎的DMO(discrete material optimization method)實現了一般的層壓復合殼的纖維優化.Niu 等[21]將DMO 法拓展到考慮最小聲輻射的層壓復合板結構的離散纖維角度優化問題中.Gao 等[22]將復合材料結構的離散纖維優化問題轉化為可用連續體拓撲優化方法處理的材料選擇問題,提出了一種有效的二值參數化列式.段尊義等[23]將連續化懲罰策略與Heaviside 懲罰函數弓入傳統DMO 模型中,提高了復合材料單層板結構的纖維角度優化的收斂率.他們還將所提方法拓展到纖維角度和結構拓撲構型一體化優化問題中,以最大化結構的剛度[24].Xia 和Shi[25]基于Shepard 插值方法[26]構造了一個連續分布的纖維場,研究纖維增強復合材料結構的剛度最大化問題.Yan 等[27]以最大結構剛度為目標,基于雙向漸進結構優化法提出了一種實現正交材料纖維方向和結構拓撲構型的同步優化.Xu 等[28]實現了考慮載荷不確定性的纖維增強復合結構的纖維方向的魯棒性優化.Papapetrou等[29]將SIMP 方法和水平集方法結合,實現了纖維角度和結構拓撲的同步優化,以最大化結構的剛度.Luo 等[30]針對纖維增強復合材料的剛度最大化問題,提出了離散-連續的參數化列式,可有效減小優化問題陷入局部優化解的風險.

可以看出,現有文獻大部分關注單一均質材料結構的動力學拓撲構型的優化、或單純考慮纖維角度優化問題,另外,關于纖維增強復合材料結構的拓撲構型和纖維方向的同步優化方面的研究大部分關注剛度最大化優化問題,而在動力學背景下考慮結構拓撲構型和纖維方向同步優化的研究較少.因此,本文旨在提出一種針對連續纖維增強復合材料結構無阻尼自由振動下的基頻最大化問題的拓撲優化方法,實現結構拓撲構型與纖維方向的同步優化,以最大化結構的第一階本征頻率.

1 連續纖維增強復合材料結構優化模型

1.1 復合材料結構自由振動有限元分析

考慮二維正交各向異性材料線彈性結構無阻尼自由振動分析的離散型有限元方程表示如下

其中,ui為第i階特征頻率ωi所對應的特征向量.K和M分別為結構的總體剛度陣和質量陣,由對應的單元剛度陣和質量陣組裝形成.和分別由式(2)和式(3)確定

其中,h為結構的厚度,ρ0為材料的質量密度,N為雙線性等參元形函數組成的矩陣.B為應變位移矩陣.和分別為單元e所對應的單元質量陣和剛度陣.θe為單元e所對應的表征纖維方向的角度設計變量,即建立在單元e上的纖維主方向和垂直于纖維主方向上的局部坐標X′-Y′與全局坐標X-Y之間的旋轉角度,如圖1 所示.

圖1 坐標轉換示意圖Fig.1 Schematic illustration of coordinate transformation

式(3)中的C為全局坐標下的正交本構材料矩陣,其表達式為

其中,C0為局部坐標下的正交本構材料矩陣,由沿纖維軸向的彈性剛度Ex,纖維橫向彈性剛度Ey,剪切模量Gxy,泊松比Vxy,Vyx五個材料常數組成,其表達式為

式中,Ex,Ey,Vxy,Vyx滿足VyxEx=VxyEy.

式(4)中的T(θe)為C與C0之間的轉換矩陣,其表達式由下式定義

1.2 材料插值格式

為了實現結構拓撲構型與纖維角度的同步優化,這里采用SIMP 方法來實現結構拓撲構型的優化,即給每一個單元弓入一個密度設計變量xe,然后采用SIMP 材料插值模型建立單元剛度陣和質量陣與單元密度設計變量xe之間的關系

其中,Me和Ke分別為弓入表征拓撲構型的密度變量之后的單元質量陣和剛度陣,p和q分別為剛度陣和質量陣的懲罰因子,和分別由式(2) 和式(3)確定.

眾所周知,在用SIMP 求解動力學優化問題時,材料低密度區域通常會出現局部模態[3,8-9,31],這是由于在低密度區域(如xe≤0.1),基于q=1 懲罰之后的質量陣與基于p=3 懲罰之后的剛度陣比值過高造成的[3,9].避免局部模態的常用方法有:(1)修改剛度矩陣的懲罰方式[3,5]; (2) 修改質量矩陣的懲罰方式[9-10,12].在這里,為了避免低密度區域中的局部模態,采用文獻[9]中的方法,修改質量陣的懲罰方式

這里設置剛度矩陣懲罰因子p=3,對于質量矩陣,當單元密度xe≥0.1 時,設置懲罰因子為1; 當單元密度xe＜0.1 時,設置懲罰因子q=6.通過式(8) 和式(9)可知,通過修改質量陣的懲罰方式可以很好地避免底密度區域質量陣與剛度陣的比值過高的情況,從而有效避免局部模態.

1.3 優化模型

結構拓撲優化旨在滿足某些約束條件下尋找最優的材料分布,以期獲得最佳的結構性能.與文獻[32-34] 關注的優化問題不同,本文研究連續纖維增強復合材料結構在無阻尼自由振動下的基頻最大化優化問題,以實現纖維方向角度與結構拓撲構型的同步優化,其對應的數學列式可表示如下

在此優化模型中,有兩類設計變量,其中x為表征結構拓撲構型的密度設計變量集合,θ為表征纖維方向的角度設計變量集合,Ne為用來離散設計域的單元個數.ωi為第i階特征頻率,ui為其對應的特征向量.K和M分別為弓入兩類設計變量之后結構的總體剛度陣和質量陣.Ve和V0分別為單元e和整個設計域的體積,fv為設計領域中可用實體材料的體積分數比.xmin為很小的正數,作為密度設計變量的下限以防止剛度矩陣奇異,本文取xmin=10-3.本文設置每個纖維角度變量θe的范圍為[-2π,2π],寬的纖維角度范圍允許角度旋轉到最佳位置,可以避免較窄的角度變量范圍導致的潛在的局部解問題[35].

2 靈敏度分析與靈敏度過濾

動力學優化問題可能存在重特征根問題,特別是那些設計自由度多、依賴的設計參數多的結構,在優化的過程中可能會出現兩階特征根相等的情況,重特征根會導致優化問題不連續,進而導致多重特征根關于設計變量的敏度難于計算的困難.目前也有較多優秀工作討論了這一問題,并給出了有效的解決方案[9,36].重特征根問題不是這篇文章的關注點,因此,這里假設所關注的特征根均為單根,即特征根均互不相等,這樣可以推導單特征值關于兩類設計變量的靈敏度.

2.1 特征值對密度設計變量的靈敏度

式(10)中無阻尼自由振動分析的有限元方程兩邊同時對密度設計變量xe求導,并化簡可得

對式(11)兩邊同時左乘,再經過簡單化簡可得特征值λi=關于設計變量xe的靈敏度

其中,?K/?xe和?M/?xe與材料插值格式有關,具體表達式可根據式(8)和式(9)計算得到

2.2 特征值對角度設計變量的靈敏度

類似于求特征值關于密度設計變量的靈敏度,可以推導特征值關于角度設計變量θe的靈敏度,其表達式如下

值得注意的是,這里假設質量矩陣與纖維方向無關.式(15)中剛度陣K對角度設計變量θe的導數?K/?θe可表示為

關于體積約束關于設計變量的敏度可直接求導獲得,這里省略其細節推導.

2.3 靈敏度過濾

由于有限元網格離散的影響,基于單元密度的拓撲優化方法經常遭受棋盤格和網格依賴性等數值不穩定性問題.為了獲得清晰的黑白設計,本文采用式(17)對式(12)中定義的第一階特征值(最小特征值)λ1關于密度設計變量xe的敏度進行過濾

其中權重函數Hei定義為

式中,dist(e,i) 為單元e的中心和單元i的中心之間的歐式距離,Nei是指單元中心落在以單元e為中心、半徑為rmin的圓形鄰域內的單元集合.以過濾后的敏度代替式(12)參入優化迭代.

上述優化問題可以用不同的優化算法來求解,如序列線性規劃法(SLP)[37],移動漸近線算法(MMA)[38]等,本文采用MMA 算法求解上述優化問題.數值測試表明角度設計變量靈敏度的絕對值相對于密度設計變量靈敏度的絕對值要小很多,為了避免陷入局部解,提高收斂效率,本文采用式(19)定義的變步長策略,即使設計變量移動步長隨著迭代步數的增加逐漸減小

其中,movex和moveθ 分別為單元密度變量和角度設計變量的移動步長.它們的初始值分別設置為move-x=0.5,move θ=0.1.

當優化過程滿足式(20) 定義的收斂條件時,停止迭代

式中,k為當前迭代步數,ε 為容許的收斂誤差,并取ε=1.0×10-3.

3 數值算例與討論

本節通過一個靜力學優化算例和兩個算例動力學優化算例來展示本文方法的有效性.若無特殊說明,平板厚度設置為1 mm,質量密度設置為1.6 g/cm3;本構材料陣中4 個獨立的物理參數取值為Ex=385 GPa,Ey=118 GPa,Vxy=0.3 和Gxy=84 GPa.假設初始結構為完全實體材料,即所有密度設計變量xe的初始值均為1.本文所有的角度變量均以弧度制表示.本文所有算例中MMA 算法參數取值列于表1 中.

表1 MMA 子程序中參數取值Table 1 Parameters needed for MMA subroutine

3.1 算例一

算例一考慮短懸臂梁的剛度最大化優化問題.幾何尺寸為80 mm×40 mm 的設計域和邊界條件如圖2 所示,梁左端完全固支,在右端中心點處受大小為30 N 的橫向拉力.設計域由80×40 個平面應力單元來離散.材料準許的體積分數設定為0.3.

圖2 懸臂梁的設計域和邊界條件Fig.2 Design domain and boundary condition of cantilever beam

為了驗證本文方法對靜力學優化問題的有效性,考慮以0,π/6,π/2,及[-2π,2π] 區間內產生的一個隨機數為初始纖維方向.圖3 顯示了4 種不同初始纖維方向所對應的優化結果.從結果可以看出,盡管初始纖維角度不同,但優化后的纖維基本沿軸向拉伸方向分布,增強了結構的剛度,符合工程實際情況,這驗證了本文方法對靜力學優化問題的有效性.

圖3 不同初始纖維方向對應的優化結果Fig.3 Optimization results for different initial fiber orientations

3.2 算例二

算例二考慮懸臂梁動力學優化問題.設計域與算例一相同,梁左端完全固支,在右邊中心點處附加一個大小為28.8 g 的集中質量點.設置纖維初始方向為水平方向,即設置所有角度設計變量θe的初始值為0.材料準許的體積分數設定為0.4.

圖4 顯示了連續纖維增強懸臂梁在無阻尼自由振動下結構拓撲構型的迭代過程.優化算法在迭代87 步之后收斂,此時結構對應的最小特征值為188.76 Hz.圖5 顯示了懸臂梁的纖維分布情況.可以看出,本文方法可以獲得清晰的拓撲構型,優化后纖維方向基本沿懸臂梁橫向中軸線呈對稱分布.

圖4 懸臂梁拓撲構型的迭代過程Fig.4 Iterative process of topological configuration of cantilever beam

圖5 懸臂梁纖維方向的優化結果Fig.5 Optimization result of fiber orientation for cantilever beam

為了進一步驗證所提方法的有效性,在圖5 顯示的最終纖維角度的優化結果基礎上,讓所有單元對應的角度變量逆時針偏轉角度β.在幾個典型的偏轉角度π/6,π/4,π/3,5π/12,π/2 下,結構對應的最小特征值λ1分別為128.96,96.06,75.41,64.63,61.32.可以看出,纖維方向偏轉一定角度后結構對應的最小特征值均比通過本文優化方法獲得的結果(188.76 Hz)要小,這說明本文方法能夠有效找到更優的纖維方向分布,提高結構的頻率.

圖6 顯示了各向同性材料懸臂梁在無阻尼自由振動下結構拓撲構型的優化結果,其中材料本構陣中使用的獨立物理參數為:楊氏模量E=385 GPa,泊松比V=0.3.可以看出,本文獲得的拓撲構型與文獻[15] 中的動力學優化結果類似,優化算法在迭代61 步之后收斂,此時對應的最小特征值為188.87 Hz.對于當前的懸臂梁優化算例,各向同性材料拓撲優化結果與纖維增強的正交各向異性材料結構拓撲構型結果相似.

圖6 各向同性材料懸臂梁的拓撲構型Fig.6 Topological configuration of cantilever beam with isotropic material

3.3 算例三

算例三考慮固支細長梁結構的動力學優化問題.幾何尺寸為160 mm×40 mm 的設計域和邊界條件如圖7 示,梁兩端全固定,平板中心處附加一個質量為51.2 g 的質量點.設計域由160×40 個平面應力單元來離散.設置纖維角度設計變量的初始值為θe=0.材料準許的體積分數為0.4.

圖7 固支梁的設計域和邊界條件Fig.7 Design domain and boundary condition of clamped beam

圖8 顯示了連續纖維增強固支細長梁在無阻尼自由振動下結構拓撲構型的迭代過程.圖9 顯示了固支細長梁纖維方向的優化結果.從優化結果可知,本文方法可以獲得清晰的拓撲構型,纖維方向主要沿著桿件的方向分布.圖10 顯示了前兩階特征值的迭代歷史.從結果可以看出,優化算法在迭代77 步之后到達收斂,此時結構對應的最小特征值為373.86 Hz,第2 階特征值為22 869.56 Hz,前兩階特征值在優化的過程中始終分離的,第2 階特征值始終大于第1 階特征值,未產生重特征根問題,這驗證了本文關注的特征值為單根的假定.

同樣,為了進一步驗證本文方法的有效性,在圖9 顯示的最終纖維角度的優化結果基礎上,讓所有單元對應的角度變量逆時針偏轉角度β.在幾個典型的偏轉角度π/6,π/4,π/3,5π/12,π/2 下，結構對應的最小特征值λ1的值分別為262.18,198.32,157.31,135.58,128.84.由表3 可知,結果與懸臂梁相優化問題類似,纖維方向偏轉一定角度后所得結構的最小特征值均變小,這進一步說明本文方法的有效性.

圖8 固支梁拓撲構型的迭代過程Fig.8 Iterative process of topological configuration of clamped beam

圖9 固支梁纖維方向的優化結果Fig.9 Optimization result of fiber orientation for clamped beam

圖10 固支梁前2 階特征值的迭代歷史Fig.10 Iteration history of the first two eigenvalue of clamped beam

4 結論

針對連續纖維增強復合材料結構無阻尼自由振動下的基頻最大化問題,本文提出了一種有效實現拓撲構型和纖維方向同步優化的拓撲優化方法.建立了以準許的材料用量體積分數為約束、以結構的一階特征值為目標函數的動力學拓撲優化模型,詳細推導了目標函數關于密度設計變量和角度設計變量的解析靈敏度列式.為了避免陷入局部最優解,采用了變移動步長的策略.最后,通過一個靜力學優化算例和兩個動力學優化算例驗證了本文方法的有效性.結果表明,本文方法在實現結構拓撲構型與纖維角度的一體化優化的同時,可以有效改善結構的性能.