鄒元杰
(北京空間飛行器總體設(shè)計(jì)部,北京 100094)
航天器結(jié)構(gòu)或機(jī)構(gòu)中存在大量的桿系結(jié)構(gòu)和連接機(jī)構(gòu)[1-2],考慮這些結(jié)構(gòu)的特點(diǎn),在力學(xué)分析中有時(shí)采用彈簧元模擬比較方便。航空航天領(lǐng)域常用的有限元分析軟件M SC.N AS TRAN,既提供了標(biāo)量彈簧單元(Scalar Spring),也提供了6 個(gè)自由度廣義彈簧單元(Bush)。其中,標(biāo)量彈簧單元剛度系數(shù)物理意義明確,應(yīng)用起來比較容易。而廣義彈簧單元Bush 的單元原理和參數(shù)含義尚不明確。
作者在使用該單元的過程中遇到諸多問題,難以合理解釋。首先,按照M SC.NAS TRAN 的用戶手冊,Bush 元的剛度陣K0(對應(yīng)單個(gè)節(jié)點(diǎn)的6 個(gè)自由度,同下文(5)式中的Κ11或K22)形式為[3]:
其中, ku、kv、kw分別代表三個(gè)平動自由度的剛度系數(shù), kθx、kθy、kθz分別代表三個(gè)轉(zhuǎn)動自由度的剛度系數(shù)。按照這個(gè)解釋,Bush 元的剛度陣(對應(yīng)單個(gè)節(jié)點(diǎn))其實(shí)是對角陣,其6 個(gè)自由度顯然是不耦合的。而實(shí)際計(jì)算發(fā)現(xiàn),在Bush 元節(jié)點(diǎn)上施加垂直于單元軸線方向的力,會產(chǎn)生彎曲方向的轉(zhuǎn)角;同樣,施加彎矩,也會在垂向產(chǎn)生平動位移。這說明,6個(gè)自由度的剛度并非解耦,即Bush 單元剛度陣不可能如(1)式所示,必然存在耦合項(xiàng)。其次,工程上對于桿系和連接機(jī)構(gòu)建模,也常采用梁單元(主要包括歐拉梁元和剪切梁元),然而梁單元與Bush 元的對應(yīng)關(guān)系尚不清楚。再次,在計(jì)算中發(fā)現(xiàn),標(biāo)量彈簧單元的剛度僅與輸入的剛度系數(shù)有關(guān),其長度可以為0,即定義單元的兩個(gè)節(jié)點(diǎn)可以重合,但Bush 元的長度取為0 則出錯(cuò)。
欲解開上述困惑,唯一的途徑是對Bush 元的原理進(jìn)行研究。為此,我們開展了以下研究:1)設(shè)法找到Bush 元剛度陣與輸入的6 個(gè)剛度系數(shù)的關(guān)系;2)搞清楚6 個(gè)剛度系數(shù)的真實(shí)物理意義,進(jìn)而分析Bush 元與標(biāo)量彈簧單元、梁單元的關(guān)系;3)結(jié)合工程算例,對上述Bush 有限元原理進(jìn)行驗(yàn)證。
由于大多數(shù)商用軟件對用戶來說都是“黑箱子”, 很難確切知道內(nèi)核程序和算法, 因此, 得到NAS TRAN 中Bush 元剛度陣的具體表達(dá)式是非常困難的。我們只有通過大量的計(jì)算、比對,并結(jié)合相關(guān)的有限元理論知識,由數(shù)值結(jié)果反推單元剛度陣的表達(dá)式。利用DMAP 語言,我們可以輸出Bush元的剛度陣,然后,通過改變輸入?yún)?shù),觀察剛度陣的變化,進(jìn)而獲取剛度陣的解析表達(dá)式。
圖1 Bush 單元示意圖Fig.1 Sketch of Bush element
對于圖1 所示的Bush 單元,其兩端節(jié)點(diǎn)的位移向量q 和對應(yīng)的力向量R 為
單元的靜力平衡方程為
其中,單元剛度陣K 可以按兩個(gè)節(jié)點(diǎn)分解為4個(gè)分塊矩陣
設(shè)Bush 元的6 個(gè)輸入?yún)?shù)為ki(i =1, ...,6), 單元長度為L 。通過數(shù)值計(jì)算最終反推得到的各分塊矩陣解析表達(dá)式如下:
上述表達(dá)式已經(jīng)通過一些算例驗(yàn)證。從公式(6)~(8)以及相關(guān)算例的計(jì)算結(jié)果可以看出:1)對應(yīng)單個(gè)節(jié)點(diǎn)的剛度陣(如K11)并非對角陣,存在垂向位移和轉(zhuǎn)角的耦合剛度項(xiàng),而且6 個(gè)對角線元素也并非對應(yīng)輸入?yún)?shù)ki(i =1, ...,6)。2)對于純拉壓或純扭轉(zhuǎn)問題,剛度陣中只有k1 或k4 起作用,此時(shí),采用標(biāo)量彈簧元和Bush 元建模是等效的。而對于更復(fù)雜的力學(xué)問題,由于存在耦合項(xiàng),Bush元的剛度不可能由6 個(gè)相互獨(dú)立的標(biāo)量彈簧元(分別對應(yīng)節(jié)點(diǎn)6 個(gè)自由度的剛度)來代替,二者是無法等效的。3)Bush 元剛度陣元素不僅和輸入的剛度系數(shù)ki(i =1, ...,6)有關(guān),而且是Bush 單元長度L 的函數(shù)。由于Nastran 程序默認(rèn)Bush 元需要考慮耦合項(xiàng)影響,因此,耦合項(xiàng)不能為0(將k3或k2取為0 的情況除外),即式(6)~(8)中的L 不能取為0。而標(biāo)量彈簧元的剛度陣與單元長度無關(guān),因此,單元長度L 可以取為0,即兩端節(jié)點(diǎn)可以重合。
為了找到Bush 單元與梁單元的關(guān)系,對比了兩類單元的剛度陣表達(dá)式,試圖用梁單元的輸入?yún)?shù)來表示Bush 元的輸入?yún)?shù)。
梁單元的相關(guān)研究多年來持續(xù)不斷[4-9]。一般單元剛度陣的推導(dǎo)都應(yīng)用虛位移原理,但是對于梁,應(yīng)用結(jié)構(gòu)力學(xué)中的解析方法(直接法)推導(dǎo)更為方便,只要列出梁單元兩端的位移要素與外力之間的關(guān)系即可求出。剪切梁(Timoshenko 梁)單元的剛度陣如下(節(jié)點(diǎn)位移和外力的定義同Bush 元,見圖1)[10]:
在式(10)~(12)中, E 為材料的彈性模量, G為材料的剪切模量, A 為截面面積, L 為梁單元的長度, Iy、Iz分別為繞y 、z 軸的彎曲慣性矩, J 為扭轉(zhuǎn)慣性矩,φy、φz分別為y 向和z 向的剪切修正系數(shù)。
比較式(10)~(12)和式(6)~(8)可以看出,Bush 元剛度陣的非零元素和剪切梁單元是一一對應(yīng)的,說明二者所表征的物理含義有可能相同。將式(9)和式(5)逐項(xiàng)對應(yīng),可以將Bush 元的剛度參數(shù)用剪切梁的參數(shù)表達(dá)。
Bush 元各剛度系數(shù)與剪切梁參數(shù)的對應(yīng)關(guān)系為
若將Bush 元的剛度參數(shù)表達(dá)為(13)~(18)式,則式(9)和式(5)中的剛度陣完全一樣。進(jìn)一步分析可以看出,兩個(gè)剛度陣各有7 個(gè)獨(dú)立參數(shù)表征:在Bush 元的剛度陣中,7 個(gè)參數(shù)分別為ki(i =1, ...,6)和L ;在剪切梁元中,7 個(gè)參數(shù)為EA 、EI y 、EIz 、GJ 、φy 、φz 和L 。雖然用于表達(dá)的參數(shù)不同,但Bush 元和剪切梁元的剛度陣卻可以完全等價(jià)。
Bush 元各參數(shù)的物理意義也不難理解。k1 、k4分別是拉壓剛度(軸向)和扭轉(zhuǎn)剛度,這兩個(gè)參數(shù)的測量也比較容易。k5、k6分別是z 向和y 向的彎曲剛度,只要使梁處于純彎曲狀態(tài),這兩個(gè)參數(shù)容易測出。k2、k3的物理意義解釋稍為復(fù)雜,這里以k2為例加以說明。如圖1 所示,對于Bush 單元i~j ,若i 端固支, j 端約束繞z 軸轉(zhuǎn)角(θzj=0),同時(shí)施加y 向外力Fyj,此時(shí), j 端y 向位移為uyj。由式(4)可以看出:也就是說, k2 對應(yīng)懸臂梁自由端轉(zhuǎn)角被約束時(shí)自由端的垂向(y 向)剛度。如果要直接測得k2或k3,則需要約束懸臂梁自由端的轉(zhuǎn)角,測量垂向載荷和垂向位移,進(jìn)而換算得到k2 或k3 。
對于歐拉梁,由于忽略其剪切效應(yīng)影響,剪切修正系數(shù)φy、φz應(yīng)取為0。此時(shí),其對應(yīng)的Bush 元6個(gè)剛度參數(shù)中, k2 與k6 、k3 與k5 是相互關(guān)聯(lián)的。所以,與歐拉梁對應(yīng)的Bush 元只有四個(gè)剛度參數(shù)是獨(dú)立的。如果采用試驗(yàn)方法測得6 個(gè)剛度系數(shù),則k2與k6(或k3與k5)可以相互校驗(yàn)測量結(jié)果的準(zhǔn)確性。
本節(jié)針對典型的航天器結(jié)構(gòu), 應(yīng)用Bush 單元進(jìn)行建模和分析。通過對比采用Bush 元的分析結(jié)果與采用梁單元的計(jì)算結(jié)果,可以進(jìn)一步驗(yàn)證Bush元參數(shù)與剪切梁參數(shù)對應(yīng)關(guān)系的準(zhǔn)確性。
對太陽翼展開狀態(tài)的動力學(xué)分析而言,鉸鏈的建模十分關(guān)鍵。通常工程上采用梁單元進(jìn)行建?!,F(xiàn)改用Bush 元建立鉸鏈模型,替換原模型中的梁單元,而其它結(jié)構(gòu)單元屬性均不改變。利用式(13)~(18),由梁單元參數(shù)確定Bush 元參數(shù),研究計(jì)算結(jié)果的變化。由于Bush 單元沒有質(zhì)量屬性,采用Bush 元時(shí), 需要在節(jié)點(diǎn)增加集中質(zhì)量(Lumped mass)單元,集中質(zhì)量大小取為梁單元質(zhì)量的一半。動力分析時(shí),兩種模型的質(zhì)量陣均采用集中質(zhì)量處理。圖2 為某太陽翼的有限元模型。
圖2 太陽翼有限元模型Fig.2 Finite element model for solar array
首先計(jì)算靜力響應(yīng)。將太陽翼根鉸與太陽翼驅(qū)動機(jī)構(gòu)(SADA)連接處固定支撐,在整個(gè)結(jié)構(gòu)上施加1 gn的慣性加速度載荷(垂直于面板的方向,即Z向),計(jì)算太陽翼最外端,節(jié)點(diǎn)898 處的靜力響應(yīng)。對鉸鏈分別采用Bush 元和剪切梁元建模時(shí)節(jié)點(diǎn)898 處的靜力響應(yīng)對比情況見表1。從表中數(shù)據(jù)看,兩者的計(jì)算結(jié)果是完全一致的。T1, T2, T3分別對應(yīng)于沿X、Y 、Z 軸方向的線位移,R1, R2,R3分別對應(yīng)于繞X 、Y 、Z 軸的轉(zhuǎn)角。
表1 太陽翼的靜力分析結(jié)果Table 1 Static analysis results for solar array
然后,進(jìn)行模態(tài)分析。邊界條件同靜力分析。分析結(jié)果表明:兩種方法計(jì)算的各階模態(tài)頻率非常接近(見表2),振型也完全一致。
下面計(jì)算一個(gè)更為復(fù)雜的桁架結(jié)構(gòu)。針對某大型桁架結(jié)構(gòu)(見圖3),對所有的連接桿結(jié)構(gòu)分別采用梁單元和Bush 元建模,對比兩種方法的靜力和模態(tài)分析結(jié)果。該桁架結(jié)構(gòu)需要采用2 816 個(gè)Bush 元或梁單元建模。
表2 太陽翼的模態(tài)分析結(jié)果Table 2 Modal analysis results for solar array
首先在桁架左端固支,在右端某節(jié)點(diǎn)加靜載荷(Z 向)1N,采用兩種單元處理方法計(jì)算加載點(diǎn)的靜力響應(yīng)。加載點(diǎn)的靜力響應(yīng)見表3。T1,T2,T3分別對應(yīng)于沿X、Y 、Z 軸方向的線位移,R1,R2,R3分別對應(yīng)于繞X 、Y 、Z 軸的轉(zhuǎn)角。從數(shù)據(jù)看,Bush 元的計(jì)算結(jié)果與剪切梁單元非常接近。另外,進(jìn)行了桁架結(jié)構(gòu)的模態(tài)分析,分析結(jié)果表明:兩種方法計(jì)算的各階模態(tài)頻率(見表4)和振型均吻合較好。
圖3 桁架結(jié)構(gòu)圖Fig.3 Structure of truss
表3 桁架靜力分析結(jié)果Table 3 Static analysis results for truss
表4 桁架模態(tài)分析結(jié)果Table 4 Modal analysis results for truss
上述兩個(gè)工程應(yīng)用算例對Bush 元的靜力響應(yīng)和動力特性進(jìn)行了分析,計(jì)算結(jié)果表明本文對Bush元剛度陣的解析表達(dá)式推導(dǎo)是準(zhǔn)確的,同時(shí)所獲得的Bush 單元與梁單元的對應(yīng)關(guān)系也是正確的。
本文利用數(shù)值計(jì)算結(jié)果確定了Bush 元剛度陣的解析表達(dá)式,進(jìn)而,通過對比Bush 元與剪切梁元的剛度陣,推導(dǎo)了Bush 元參數(shù)與剪切梁剛度參數(shù)的對應(yīng)關(guān)系,并利用展開狀態(tài)太陽翼結(jié)構(gòu)和大型桁架結(jié)構(gòu)的靜力和動力分析驗(yàn)證了相關(guān)推導(dǎo)。主要結(jié)論如下:
1)Bush 元和剪切梁元各有7 個(gè)獨(dú)立的參數(shù)。雖然兩組參數(shù)不同,但Bush 元和剪切梁元的剛度陣可以完全等價(jià)。
2)Bush 元對應(yīng)單個(gè)節(jié)點(diǎn)的剛度陣并非對角陣,存在垂向位移和轉(zhuǎn)角的耦合剛度項(xiàng),而且對角線元素也并非輸入?yún)?shù)k i(i =1, ...,6);由于存在耦合項(xiàng),一般來說,Bush 元的剛度不可能由6 個(gè)標(biāo)量彈簧元(分別對應(yīng)節(jié)點(diǎn)6 個(gè)自由度的剛度)來代替,二者無法等效。
3)對于歐拉梁來說,剪切效應(yīng)可忽略,此時(shí),與其剛度等價(jià)的Bush 元只有4 個(gè)獨(dú)立的剛度系數(shù)。如果采用試驗(yàn)方法測得6 個(gè)剛度系數(shù),則6 個(gè)剛度系數(shù)中有兩對系數(shù)(k2與k6、k3與k5)可以相互校驗(yàn)測量結(jié)果的準(zhǔn)確性。
)
[1]陳烈民.航天器結(jié)構(gòu)與機(jī)構(gòu)[M].北京:中國科學(xué)技術(shù)出版社,2005
[2]Thomas P S, Wiley J L.Spacecraft structures and mechanisms-from concept to launch [M].Torrance,California:Microcosm, Inc., 1995
[3]MSC.NAS TRAN 2001 Reference Manual[Z].MSC.Software Corporation, 2002
[4]COOK R D.Concepts and applications of finite element analysis(Second Edition)[M].H oboken:John Wiley &Sons,1981
[5]張阿舟,諸德超,姚起杭, 等.實(shí)用振動工程(1)——振動理論與分析[M].北京:航空工業(yè)出版社, 1996
[6]邱吉寶, 向樹紅, 張正平.計(jì)算結(jié)構(gòu)動力學(xué)[M].合肥:中國科學(xué)技術(shù)大學(xué)出版社, 2009
[7]李國強(qiáng), 李進(jìn)軍.Timoshenko-Euler 楔形梁有限元[J].力學(xué)季刊,2002, 23(1):64-69
[8]朱炳麒,陳學(xué)宏.理性Timoshenko 梁單元及其應(yīng)用[J].力學(xué)與實(shí)踐, 2008, 30 (1):31-34
[9]王勖成, 邵敏.有限單元法基本原理和數(shù)值方法(第2版)[M].北京:清華大學(xué)出版社, 1997
[10]金咸定,趙德有.船體振動學(xué)[M].上海:上海交通大學(xué)出版社,1997