張君茹， 程耿東

(大連理工大學(xué) 工程力學(xué)系，大連 116024)

1 引 言

旋轉(zhuǎn)結(jié)構(gòu)在工業(yè)中有著廣泛的應(yīng)用，如渦輪發(fā)動機、直升機、風(fēng)力發(fā)動機葉片和旋轉(zhuǎn)太陽帆。當(dāng)旋轉(zhuǎn)結(jié)構(gòu)工作時，可能會受到外界環(huán)境的干擾或者沖擊，如異物撞擊或者砂粒等[1,2]，從而產(chǎn)生偏離正常運行的振動，影響旋轉(zhuǎn)結(jié)構(gòu)的正常功能，研究沖擊荷載下的旋轉(zhuǎn)結(jié)構(gòu)減振問題是十分有必要的。

旋轉(zhuǎn)結(jié)構(gòu)動力學(xué)理論研究工作從研究旋轉(zhuǎn)梁的動力學(xué)特性和響應(yīng)起步，單桿機械臂或者葉片的結(jié)構(gòu)分析和優(yōu)化是以梁結(jié)構(gòu)的研究為基礎(chǔ)[3-5]，但是梁結(jié)構(gòu)不能很好地模擬具有低展弦比的結(jié)構(gòu)，在很多應(yīng)用中，研究旋轉(zhuǎn)板的動力學(xué)特性和優(yōu)化方法更合理必要。旋轉(zhuǎn)板建模方法的主流方法有浮動坐標(biāo)法和絕對節(jié)點法兩類?；诟幼鴺?biāo)法，文獻[6,7]研究了旋轉(zhuǎn)薄板的動力學(xué)響應(yīng)和模態(tài)特性，Hashemi等[8]研究了旋轉(zhuǎn)厚板的動力學(xué)響應(yīng)和模態(tài)特性。趙將等[9]基于絕對節(jié)點法研究了旋轉(zhuǎn)薄板的頻率和模態(tài)特性，但是在處理斜率不連續(xù)問題時仍然需要進一步研究[10]，優(yōu)化設(shè)計過程中出現(xiàn)相鄰單元密度不相等造成斜率不連續(xù)情況。本文采用浮動坐標(biāo)法建立平面旋轉(zhuǎn)板的多體動力學(xué)數(shù)學(xué)模型。

旋轉(zhuǎn)板結(jié)構(gòu)的減振優(yōu)化設(shè)計和研究工作相對少見。劉利軍等[11]研究了轉(zhuǎn)動ACLD懸臂板的動力學(xué)建模及主動約束層阻尼層的振動控制。孫家亮等[12]基于移動可變形組件法(MMC)和絕對節(jié)點法(ANCF)研究了旋轉(zhuǎn)薄板頻率拓撲優(yōu)化方法。

動力學(xué)減振優(yōu)化的研究按照目標(biāo)函數(shù)或者約束函數(shù)分為兩類。一類是基于結(jié)構(gòu)動力學(xué)特性如頻率或者模態(tài)，這方面的工作很多；另一類是基于動力學(xué)響應(yīng)。Kang等[13]綜述了基于響應(yīng)的優(yōu)化工作，指出基于結(jié)構(gòu)響應(yīng)的優(yōu)化問題需要求解每個時間積分步上的動力學(xué)響應(yīng)，求解每個時間步響應(yīng)的敏度十分耗費時間。基于動力學(xué)響應(yīng)但不考慮每一個時刻的響應(yīng)，趙君鵬等[14]提出基于響應(yīng)的積分形式目標(biāo)函數(shù)，閻琨等[15,16]使用積分形式的二次指標(biāo)作為目標(biāo)函數(shù)，使用李亞普洛夫方法將該指標(biāo)簡化，并使用伴隨法計算該指標(biāo)的敏度,求解了受沖擊荷載板結(jié)構(gòu)阻尼拓撲優(yōu)化。

本文結(jié)合多體動力學(xué)理論和有限元方法推導(dǎo)了常速旋轉(zhuǎn)板的面內(nèi)振動動力學(xué)方程，并采用二次型性能指標(biāo)作為目標(biāo)函數(shù)，結(jié)合李亞普洛夫方法，建立以單元人工密度作為設(shè)計變量的SIMP(Solid Isotropic Material with Penalty)結(jié)構(gòu)拓撲優(yōu)化列式，研究在給定材料體積的條件下，受到?jīng)_擊荷載作用的旋轉(zhuǎn)板的拓撲優(yōu)化減振設(shè)計。本文還討論了以二次型性能指標(biāo)和以基頻為目標(biāo)的優(yōu)化設(shè)計的異同，發(fā)現(xiàn)在遠端有環(huán)向支撐時，板基頻最大的設(shè)計是重頻率設(shè)計，但是，以二次型性能指標(biāo)作為目標(biāo)函數(shù)的優(yōu)化設(shè)計沒有重頻。

2 垂直于旋轉(zhuǎn)軸的平板動力學(xué)方程

2.1 浮動坐標(biāo)法下廣義位置坐標(biāo)

考慮板以Z軸為旋轉(zhuǎn)軸在平面內(nèi)大范圍轉(zhuǎn)動，并假設(shè)發(fā)生的面內(nèi)變形是小變形，使用平面應(yīng)力板單元來離散板結(jié)構(gòu)和模擬板的變形[8]。任取板中一點P，點P到板結(jié)構(gòu)坐標(biāo)系ObXbYb原點Ob的位置矢量為uP，點P到全局坐標(biāo)系OXY原點O的位置矢量為r，點Ob在全局坐標(biāo)系OXY的位置矢量為R，旋轉(zhuǎn)板繞Zb軸旋轉(zhuǎn)θ角，點P發(fā)生的變形為u。則P點全局位置矢量r表示為

r=R+T(uP+u)

(1)

(2)

點P變形表示為

u=Nqe

(3)

式中平面應(yīng)力板單元自由度qe和形函數(shù)N的表達式參見文獻[17]。旋轉(zhuǎn)板結(jié)構(gòu)上任意點P的全局位置矢量r由廣義位置坐標(biāo)qP唯一確定表示

(4)

式中結(jié)構(gòu)自由度q由單元自由度qe組集而成。

圖1 板中任取點P的位置

2.2 旋轉(zhuǎn)板的質(zhì)量和剛度矩陣

點P速度由該點的位置矢量對時間求導(dǎo)得

(5)

旋轉(zhuǎn)板的動能和質(zhì)量矩陣表示為

(6)

式中M(qP)為結(jié)構(gòu)有限元模型的單元質(zhì)量矩陣，其分量表達式參見文獻[17]。

板的彈性勢能和剛度矩陣可以表示為

(7)

式中應(yīng)變矩陣B和彈性矩陣D參見文獻[18]。

2.3 穩(wěn)定旋轉(zhuǎn)板受沖擊擾動動力學(xué)方程

動力學(xué)問題的拉格朗日方程為

(8)

(9)

(10)

(i=1,2,…)

(11)

式中結(jié)構(gòu)的阻尼比ξ是與材料的性質(zhì)或者應(yīng)力的狀態(tài)有關(guān)的，一般在0.5%～15%之間[19]。本文采用α和β為 10和1e -5。使得對于本文算例，與基頻相應(yīng)的阻尼比落在此區(qū)間。

假定旋轉(zhuǎn)板從初始啟動后經(jīng)過一段時間進入穩(wěn)定旋轉(zhuǎn)，轉(zhuǎn)速ω為常數(shù)，在轉(zhuǎn)動中設(shè)想受到外物碰撞、打擊或者爆炸等作用時間極短的沖擊，假設(shè)沖擊后結(jié)構(gòu)還處于未變形狀態(tài)，并將沖擊荷載簡化為在板沖擊點處施加沖擊速度[20]。受到?jīng)_擊穩(wěn)定旋轉(zhuǎn)板的擾動動力學(xué)方程可表示為

(12)

將式(12)轉(zhuǎn)化為狀態(tài)方程，可表示為

(13)

式中A為狀態(tài)矩陣，x為由擾動速度和擾動位移組成的狀態(tài)向量。

3 旋轉(zhuǎn)板的減振拓撲優(yōu)化

研究均勻厚度為h的旋轉(zhuǎn)板的拓撲優(yōu)化，待設(shè)計的板結(jié)構(gòu)所在的設(shè)計域為Ω。假定板的材料性質(zhì)和邊界條件、允許使用的板的材料體積、沖擊荷載的位置和大小以及觀察點的位置與方向都已經(jīng)給定。穩(wěn)定旋轉(zhuǎn)板在受沖擊后發(fā)生振動，本文目的優(yōu)化結(jié)構(gòu)的拓撲以降低觀察點的振動。

3.1 以二次型性能指標(biāo)為目標(biāo)的優(yōu)化列式

(14a)

(14b)

(14c)

(e=1,2,…,n) (14d)

(15)

式中uD為人為選取受沖擊板上觀察點D在體坐標(biāo)系下Yb方向的擾動位移，沖擊點和觀察點可以相同或者不同，算例1中觀察點D的位置如圖2所示。W為權(quán)系數(shù)矩陣，除了對角線上W(uD,uD)的值等于1，其余元素均等于0。在拓撲優(yōu)化過程中，材料的密度分布在不斷發(fā)生變化，為了保證結(jié)構(gòu)受到外界沖量恒定，在沖擊點處設(shè)置不可設(shè)計域。

3.2 拓撲優(yōu)化數(shù)值方法

由于優(yōu)化列式(14)允許單元人工密度取0～1的中間值，所以可以使用基于梯度的高效優(yōu)化方法求解，但是最優(yōu)結(jié)構(gòu)的密度應(yīng)該是0～1分布，為此，SIMP方法通常假定單元材料彈性模量和人工密度的關(guān)系為

(16)

為了盡量減少在優(yōu)化結(jié)果中出現(xiàn)大量的中間密度單元，通常將懲罰系數(shù)p取為3，但是，對于靜力柔度優(yōu)化問題，已經(jīng)證明p=3時拓撲優(yōu)化問題非凸，優(yōu)化迭代可能落入局部最優(yōu)解。此外，對于動力學(xué)問題，p=3也導(dǎo)致局部振型，需要特殊的處理方法。為了避免優(yōu)化解出現(xiàn)棋盤格式等問題，需要采用線性密度過濾和非線性密度過濾等一系列變量變換，這樣的方法稱為三場SIMP(three field SIMP)法[21]。

本文采用p=3及上述方法求解了一些例題，但是，由于本文研究的旋轉(zhuǎn)結(jié)構(gòu)動力學(xué)問題不同于結(jié)構(gòu)靜力分析，優(yōu)化往往導(dǎo)致局部最優(yōu)解。本文采用p=1，除了減少了落入局部最優(yōu)解的可能，還避免了局部振動模式，得到的中間密度在本問題中也可以解釋為板的厚度取了中間值。為了解決優(yōu)化解有灰色單元的情況，引入文獻[14]包含灰度信息的加權(quán)目標(biāo)函數(shù)為

(17)

式中w為權(quán)系數(shù)，該值在優(yōu)化過程中固定，為了統(tǒng)一量級，該值的計算方式取為初始設(shè)計的二次型指標(biāo)除以灰度信息。

4 基于李亞普洛夫方法的拓撲優(yōu)化列式和方法

為了求解式(14)的優(yōu)化問題，使用基于敏度的優(yōu)化方法求解優(yōu)化問題，需要高效的計算目標(biāo)函數(shù)及其靈敏度的方法。式(15)目標(biāo)函數(shù)是一個在無窮時間區(qū)間上的積分，需要時程分析和數(shù)值積分，計算目標(biāo)及其靈敏度十分耗費時間。

4.1 李亞普洛夫方法第二方法

根據(jù)李亞普洛夫第二方法，對于漸近穩(wěn)定的動力學(xué)系統(tǒng)，結(jié)構(gòu)在初始時刻受到x(0)的沖擊響應(yīng)后，由于阻尼而使系統(tǒng)的振動逐漸減小，無窮時刻結(jié)構(gòu)位移為0，因此x(∞)取為0。經(jīng)簡單推導(dǎo)可以將式(14a)的目標(biāo)函數(shù)轉(zhuǎn)化為

J=x(0)TPx(0)

(18)

式中P矩陣是對稱正定且滿足式(19)的代數(shù)黎卡提(Ricatti)方程，

ATP+PA=-Q

(19)

式中A矩陣為狀態(tài)方程(13)的狀態(tài)矩陣。利用式(18，19)可以把式(15)改寫為目標(biāo)函數(shù)容易計算的優(yōu)化列式，不再需要時程積分。

4.2 振型疊加法

當(dāng)結(jié)構(gòu)劃分為很多單元時，結(jié)構(gòu)矩陣維度高，采用振型疊加法將高維矩陣縮減為低維矩陣，可縮短求解P矩陣的時間。振型疊加法先求解模態(tài)方程，

(i=1,2,…)(20)

式中ωi為第i階頻率，φi為第i階模態(tài)。約定模態(tài)矩陣Φ=[φ1，φ2，…]相對于質(zhì)量矩陣M正交歸一化，即ΦTMΦ為單位陣。從模態(tài)矩陣Φ取出前s階模態(tài)，稱為縮減模態(tài)，形成矩陣Φs并引入模態(tài)坐標(biāo)φ，使用縮減模態(tài)可將板的響應(yīng)近似為

(21)

將式(20)代入式(13)得到縮減的狀態(tài)方程為

(22)

(23)

結(jié)合李亞普洛夫方法和振型疊加法，旋轉(zhuǎn)板拓撲優(yōu)化列式(14)可轉(zhuǎn)化為

(24a)

(24b)

(14c～14d)

(24c)

5 算 例

在動力學(xué)減振優(yōu)化中常采用以最大化基頻為目標(biāo)函數(shù)減小結(jié)構(gòu)的振動，本文則采用最小化二次型性能指標(biāo)為目標(biāo)函數(shù)?；l反映的是整個板結(jié)構(gòu)的振動情況，二次型性能指標(biāo)直接反映的是觀察點處的振動情況。下面的第一個算例將比較和討論不同轉(zhuǎn)速下基于二次型性能指標(biāo)的優(yōu)化設(shè)計和基于最大化基頻的優(yōu)化設(shè)計。

5.1 算例參數(shù)、阻尼參數(shù)和縮減模態(tài)階數(shù)選擇

圖2給出了旋轉(zhuǎn)板的幾何尺寸和工況。板長為 0.4 m，寬為 0.1 m，厚為 0.05 m，密度為 2766.67 kg/m3，彈性模量為6.89×1010Pa，泊松比為 0.3，體分比為 0.5。薄板結(jié)構(gòu)繞A點以轉(zhuǎn)速ω旋轉(zhuǎn)，板的左端固定，右端點B和點C為滑動支座，允許板沿相對體坐標(biāo)系的水平方向移動，即允許徑向移動，但不允許其沿環(huán)向移位。在外徑處有圓環(huán)提供支承的輪盤葉片就可以近似模擬成這樣的結(jié)構(gòu)。假定在板穩(wěn)定旋轉(zhuǎn)后，點D受到沿環(huán)向的初始沖擊速度v0，v0設(shè)置為1 m/s，初始擾動位移設(shè)置為0。為了減小撞擊點的位移，選取撞擊點D為目標(biāo)函數(shù)的減振觀察點，如式(15)所示。有限元分析時矩形板劃分為120×30，在撞擊點處設(shè)置長寬各為2個單元、密度為1的不可設(shè)計域。

圖2 旋轉(zhuǎn)板的幾何尺寸和工況

本文采用振型疊加法縮減黎卡提方程，縮減模態(tài)階數(shù)選取越多，目標(biāo)函數(shù)越精確，但選擇過少的模態(tài)數(shù)會導(dǎo)致基于模態(tài)縮減的多體系統(tǒng)動力學(xué)拓撲優(yōu)化中優(yōu)化算法不收斂和灰色區(qū)域的出現(xiàn)[22]。另一方面，選擇的模態(tài)過多，縮減的結(jié)構(gòu)矩陣和狀態(tài)矩陣的規(guī)模越大，黎卡提方程求解時間越長。由前面討論可知，高階模態(tài)的頻率如果使得相應(yīng)的阻尼比ξ大于臨界阻尼比1，相應(yīng)模態(tài)的振動不再振蕩，并不影響結(jié)構(gòu)的功能。算例中選取的模態(tài)階數(shù)使得選取模態(tài)的阻尼比都小于1。經(jīng)過計算，本算例前30階模態(tài)的阻尼比都小于1，故優(yōu)化時選取30階縮減模態(tài)。

5.2 優(yōu)化結(jié)果

初始設(shè)計取滿足約束條件的可行設(shè)計，采用國際單位制，目標(biāo)函數(shù)在數(shù)值計算時乘以1020。表2的第2～4行給出了不同轉(zhuǎn)速下的最優(yōu)指標(biāo)值、優(yōu)化結(jié)構(gòu)的頻率及其構(gòu)型。第5行max Frquency給出了基頻最大化的優(yōu)化結(jié)構(gòu)，這是一個基頻為重頻率的優(yōu)化設(shè)計。最大化基頻問題優(yōu)化列式可表示為

(25a)

β≤ωi(i=1,2,…,nr)

(25b)

(14c～14d)

(25c)

式中nr為計算時控制的特征值數(shù)，本文選為3，為了避免模態(tài)切換，采用了重特征值的靈敏度分析方法，具體細節(jié)參見文獻[23]。

由表2可知，隨著轉(zhuǎn)速的增加，指標(biāo)優(yōu)化得到具有不同基頻特性的指標(biāo)設(shè)計。由表2第4列可知，轉(zhuǎn)速0 rad/s和2000 rad/s的優(yōu)化設(shè)計及其指標(biāo)值或者基頻都非常接近。當(dāng)轉(zhuǎn)速增加到4000 rad/s，優(yōu)化設(shè)計端部的材料加強明顯。最大基頻設(shè)計在轉(zhuǎn)速0 rad/s，2000 rad/s和4000 rad/s下的指標(biāo)值分別為 1854.85，1442.69 和 1564.16，比指標(biāo)最小化設(shè)計的指標(biāo)值大約50～60倍，說明對于本問題，基頻優(yōu)化設(shè)計并不是在沖擊荷載下點D振動最小的設(shè)計。圖3給出了初始設(shè)計和轉(zhuǎn)速為2000 rad/s的指標(biāo)設(shè)計在點D的殘余振動的時程曲線，可以看出指標(biāo)設(shè)計比初始設(shè)計的點D殘余振動衰減要快，說明采用的最小化指標(biāo)能夠減小旋轉(zhuǎn)結(jié)構(gòu)的振動。由表2可知，基頻最大設(shè)計和指標(biāo)最小化設(shè)計的基頻和構(gòu)型不一樣，這說明以二次型性能指標(biāo)作為目標(biāo)函數(shù)的優(yōu)化設(shè)計沒有重頻。

表2 長寬比為4∶1時指標(biāo)設(shè)計和基頻設(shè)計的比較

圖3 轉(zhuǎn)速為2000 rad/s時初始設(shè)計和指標(biāo)設(shè)計點D的Y方向位移響應(yīng)對比

5.3 邊長比、沖擊點、觀察點及支撐條件對優(yōu)化結(jié)果的影響

為了考慮邊長比、沖擊點、觀察點及支撐條件對優(yōu)化結(jié)果的影響，研究了圖4列出的不同問題。圖4(a)的邊界條件不同于圖2結(jié)構(gòu)，該板遠端自由，且在該端中點施加集中質(zhì)量0.5 kg，其左端考慮了只有1/3邊界固支和全部固支兩種情況。表3 給出了相應(yīng)于這兩種情況的優(yōu)化設(shè)計，可以看出當(dāng)左端部分固支時，相較于全部固支，指標(biāo)設(shè)計的材料會向固支點移動。圖4(b，c)的問題采用 圖2 的算例參數(shù)，但是沖擊點和觀察點的位置改變。圖4(b)的問題，沖擊點和觀察點位置不同。圖4(c)的問題，觀察點和沖擊點都在板的中心，表4 給出了相應(yīng)的優(yōu)化設(shè)計。

圖4 轉(zhuǎn)速為2000 rad/s旋轉(zhuǎn)板不同工況或者邊界條件

由表4可知，當(dāng)沖擊點和指標(biāo)點不重合時，指標(biāo)設(shè)計很難得到清晰的優(yōu)化結(jié)果，這樣的問題本質(zhì)上和柔性機構(gòu)優(yōu)化設(shè)計相似，后者的荷載作用點和位移控制點不一致，獲得黑白結(jié)構(gòu)比較困難；當(dāng)沖擊點與指標(biāo)點重合移動到其他點時，即圖4(c)的情況也很難得到清晰的優(yōu)化結(jié)果，還需要進一步研究其改進方法。此外，本文還計算了長寬比為8∶1的旋轉(zhuǎn)板指標(biāo)優(yōu)化，該算例采用表1的荷載參數(shù)和材料參數(shù)，板的寬度變?yōu)?.05 m，長度仍為0.4 m，網(wǎng)格劃分為160×20，優(yōu)化結(jié)果列入表5。與長寬比為4∶1不同的是，長寬比為8∶1的指標(biāo)設(shè)計基頻和指標(biāo)不是隨轉(zhuǎn)速增加而增加，值得進一步研究。

表3 加集中質(zhì)量的板左端全部固支和部分固支的優(yōu)化結(jié)果比較

表4 選取指標(biāo)中不同觀察點或者沖擊點位置時的優(yōu)化結(jié)果比較

表5 長寬比為8∶1時指標(biāo)設(shè)計和基頻設(shè)計的比較

6 結(jié) 論

本文研究了繞定軸旋轉(zhuǎn)平板結(jié)構(gòu)拓撲優(yōu)化設(shè)計，采用二次型性能指標(biāo)為目標(biāo)函數(shù)以減小受到?jīng)_擊荷載時平板結(jié)構(gòu)的振動，得到了比基頻最大的拓撲優(yōu)化設(shè)計好很多的設(shè)計，且指標(biāo)最小化設(shè)計的拓撲構(gòu)型與基頻最大的設(shè)計相差較大。除此，指標(biāo)設(shè)計隨著轉(zhuǎn)速的增加，端部材料得到加強且基頻增加，間接地反映了旋轉(zhuǎn)帶來的慣性效應(yīng)，還受到邊界條件、沖擊荷載的位置、觀察點的位置和板的長寬比等的影響得到不同的構(gòu)型，其中沖擊點或者指標(biāo)點改變有可能很難得到清晰的優(yōu)化構(gòu)型，需要進一步研究。