Bush 有限單元原理及其在航天器結構建模中的應用

鄒元杰

(北京空間飛行器總體設計部,北京 100094)

1 引言

航天器結構或機構中存在大量的桿系結構和連接機構[1-2],考慮這些結構的特點,在力學分析中有時采用彈簧元模擬比較方便。航空航天領域常用的有限元分析軟件M SC.N AS TRAN,既提供了標量彈簧單元(Scalar Spring),也提供了6 個自由度廣義彈簧單元(Bush)。其中,標量彈簧單元剛度系數物理意義明確,應用起來比較容易。而廣義彈簧單元Bush 的單元原理和參數含義尚不明確。

作者在使用該單元的過程中遇到諸多問題,難以合理解釋。首先,按照M SC.NAS TRAN 的用戶手冊,Bush 元的剛度陣K0(對應單個節點的6 個自由度,同下文(5)式中的Κ11或K22)形式為[3]:

其中, ku、kv、kw分別代表三個平動自由度的剛度系數, kθx、kθy、kθz分別代表三個轉動自由度的剛度系數。按照這個解釋,Bush 元的剛度陣(對應單個節點)其實是對角陣,其6 個自由度顯然是不耦合的。而實際計算發現,在Bush 元節點上施加垂直于單元軸線方向的力,會產生彎曲方向的轉角;同樣,施加彎矩,也會在垂向產生平動位移。這說明,6個自由度的剛度并非解耦,即Bush 單元剛度陣不可能如(1)式所示,必然存在耦合項。其次,工程上對于桿系和連接機構建模,也常采用梁單元(主要包括歐拉梁元和剪切梁元),然而梁單元與Bush 元的對應關系尚不清楚。再次,在計算中發現,標量彈簧單元的剛度僅與輸入的剛度系數有關,其長度可以為0,即定義單元的兩個節點可以重合,但Bush 元的長度取為0 則出錯。

欲解開上述困惑,唯一的途徑是對Bush 元的原理進行研究。為此,我們開展了以下研究:1)設法找到Bush 元剛度陣與輸入的6 個剛度系數的關系;2)搞清楚6 個剛度系數的真實物理意義,進而分析Bush 元與標量彈簧單元、梁單元的關系;3)結合工程算例,對上述Bush 有限元原理進行驗證。

2 Bush 單元剛度陣表達式的獲取

由于大多數商用軟件對用戶來說都是“黑箱子”, 很難確切知道內核程序和算法, 因此, 得到NAS TRAN 中Bush 元剛度陣的具體表達式是非常困難的。我們只有通過大量的計算、比對,并結合相關的有限元理論知識,由數值結果反推單元剛度陣的表達式。利用DMAP 語言,我們可以輸出Bush元的剛度陣,然后,通過改變輸入參數,觀察剛度陣的變化,進而獲取剛度陣的解析表達式。

圖1 Bush 單元示意圖Fig.1 Sketch of Bush element

對于圖1 所示的Bush 單元,其兩端節點的位移向量q 和對應的力向量R 為

單元的靜力平衡方程為

其中,單元剛度陣K 可以按兩個節點分解為4個分塊矩陣

設Bush 元的6 個輸入參數為ki(i =1, ...,6), 單元長度為L 。通過數值計算最終反推得到的各分塊矩陣解析表達式如下:

上述表達式已經通過一些算例驗證。從公式(6)～(8)以及相關算例的計算結果可以看出:1)對應單個節點的剛度陣(如K11)并非對角陣,存在垂向位移和轉角的耦合剛度項,而且6 個對角線元素也并非對應輸入參數ki(i =1, ...,6)。2)對于純拉壓或純扭轉問題,剛度陣中只有k1 或k4 起作用,此時,采用標量彈簧元和Bush 元建模是等效的。而對于更復雜的力學問題,由于存在耦合項,Bush元的剛度不可能由6 個相互獨立的標量彈簧元(分別對應節點6 個自由度的剛度)來代替,二者是無法等效的。3)Bush 元剛度陣元素不僅和輸入的剛度系數ki(i =1, ...,6)有關,而且是Bush 單元長度L 的函數。由于Nastran 程序默認Bush 元需要考慮耦合項影響,因此,耦合項不能為0(將k3或k2取為0 的情況除外),即式(6)～(8)中的L 不能取為0。而標量彈簧元的剛度陣與單元長度無關,因此,單元長度L 可以取為0,即兩端節點可以重合。

3 Bush 單元與梁單元的對應關系

為了找到Bush 單元與梁單元的關系,對比了兩類單元的剛度陣表達式,試圖用梁單元的輸入參數來表示Bush 元的輸入參數。

梁單元的相關研究多年來持續不斷[4-9]。一般單元剛度陣的推導都應用虛位移原理,但是對于梁,應用結構力學中的解析方法(直接法)推導更為方便,只要列出梁單元兩端的位移要素與外力之間的關系即可求出。剪切梁(Timoshenko 梁)單元的剛度陣如下(節點位移和外力的定義同Bush 元,見圖1)[10]:

在式(10)～(12)中, E 為材料的彈性模量, G為材料的剪切模量, A 為截面面積, L 為梁單元的長度, Iy、Iz分別為繞y 、z 軸的彎曲慣性矩, J 為扭轉慣性矩,φy、φz分別為y 向和z 向的剪切修正系數。

比較式(10)～(12)和式(6)～(8)可以看出,Bush 元剛度陣的非零元素和剪切梁單元是一一對應的,說明二者所表征的物理含義有可能相同。將式(9)和式(5)逐項對應,可以將Bush 元的剛度參數用剪切梁的參數表達。

Bush 元各剛度系數與剪切梁參數的對應關系為

若將Bush 元的剛度參數表達為(13)～(18)式,則式(9)和式(5)中的剛度陣完全一樣。進一步分析可以看出,兩個剛度陣各有7 個獨立參數表征:在Bush 元的剛度陣中,7 個參數分別為ki(i =1, ...,6)和L ;在剪切梁元中,7 個參數為EA 、EI y 、EIz 、GJ 、φy 、φz 和L 。雖然用于表達的參數不同,但Bush 元和剪切梁元的剛度陣卻可以完全等價。

Bush 元各參數的物理意義也不難理解。k1 、k4分別是拉壓剛度(軸向)和扭轉剛度,這兩個參數的測量也比較容易。k5、k6分別是z 向和y 向的彎曲剛度,只要使梁處于純彎曲狀態,這兩個參數容易測出。k2、k3的物理意義解釋稍為復雜,這里以k2為例加以說明。如圖1 所示,對于Bush 單元i～j ,若i 端固支, j 端約束繞z 軸轉角(θzj=0),同時施加y 向外力Fyj,此時, j 端y 向位移為uyj。由式(4)可以看出:也就是說, k2 對應懸臂梁自由端轉角被約束時自由端的垂向(y 向)剛度。如果要直接測得k2或k3,則需要約束懸臂梁自由端的轉角,測量垂向載荷和垂向位移,進而換算得到k2 或k3 。

對于歐拉梁,由于忽略其剪切效應影響,剪切修正系數φy、φz應取為0。此時,其對應的Bush 元6個剛度參數中, k2 與k6 、k3 與k5 是相互關聯的。所以,與歐拉梁對應的Bush 元只有四個剛度參數是獨立的。如果采用試驗方法測得6 個剛度系數,則k2與k6(或k3與k5)可以相互校驗測量結果的準確性。

4 Bush 單元的工程應用

本節針對典型的航天器結構, 應用Bush 單元進行建模和分析。通過對比采用Bush 元的分析結果與采用梁單元的計算結果,可以進一步驗證Bush元參數與剪切梁參數對應關系的準確性。

4.1 太陽翼結構的分析

對太陽翼展開狀態的動力學分析而言,鉸鏈的建模十分關鍵。通常工程上采用梁單元進行建模。現改用Bush 元建立鉸鏈模型,替換原模型中的梁單元,而其它結構單元屬性均不改變。利用式(13)～(18),由梁單元參數確定Bush 元參數,研究計算結果的變化。由于Bush 單元沒有質量屬性,采用Bush 元時, 需要在節點增加集中質量(Lumped mass)單元,集中質量大小取為梁單元質量的一半。動力分析時,兩種模型的質量陣均采用集中質量處理。圖2 為某太陽翼的有限元模型。

圖2 太陽翼有限元模型Fig.2 Finite element model for solar array

首先計算靜力響應。將太陽翼根鉸與太陽翼驅動機構(SADA)連接處固定支撐,在整個結構上施加1 gn的慣性加速度載荷(垂直于面板的方向,即Z向),計算太陽翼最外端,節點898 處的靜力響應。對鉸鏈分別采用Bush 元和剪切梁元建模時節點898 處的靜力響應對比情況見表1。從表中數據看,兩者的計算結果是完全一致的。T1, T2, T3分別對應于沿X、Y 、Z 軸方向的線位移,R1, R2,R3分別對應于繞X 、Y 、Z 軸的轉角。

表1 太陽翼的靜力分析結果Table 1 Static analysis results for solar array

然后,進行模態分析。邊界條件同靜力分析。分析結果表明:兩種方法計算的各階模態頻率非常接近(見表2),振型也完全一致。

4.2 桁架結構的分析

下面計算一個更為復雜的桁架結構。針對某大型桁架結構(見圖3),對所有的連接桿結構分別采用梁單元和Bush 元建模,對比兩種方法的靜力和模態分析結果。該桁架結構需要采用2 816 個Bush 元或梁單元建模。

表2 太陽翼的模態分析結果Table 2 Modal analysis results for solar array

首先在桁架左端固支,在右端某節點加靜載荷(Z 向)1N,采用兩種單元處理方法計算加載點的靜力響應。加載點的靜力響應見表3。T1,T2,T3分別對應于沿X、Y 、Z 軸方向的線位移,R1,R2,R3分別對應于繞X 、Y 、Z 軸的轉角。從數據看,Bush 元的計算結果與剪切梁單元非常接近。另外,進行了桁架結構的模態分析,分析結果表明:兩種方法計算的各階模態頻率(見表4)和振型均吻合較好。

圖3 桁架結構圖Fig.3 Structure of truss

表3 桁架靜力分析結果Table 3 Static analysis results for truss

表4 桁架模態分析結果Table 4 Modal analysis results for truss

上述兩個工程應用算例對Bush 元的靜力響應和動力特性進行了分析,計算結果表明本文對Bush元剛度陣的解析表達式推導是準確的,同時所獲得的Bush 單元與梁單元的對應關系也是正確的。

5 結論

本文利用數值計算結果確定了Bush 元剛度陣的解析表達式,進而,通過對比Bush 元與剪切梁元的剛度陣,推導了Bush 元參數與剪切梁剛度參數的對應關系,并利用展開狀態太陽翼結構和大型桁架結構的靜力和動力分析驗證了相關推導。主要結論如下:

1)Bush 元和剪切梁元各有7 個獨立的參數。雖然兩組參數不同,但Bush 元和剪切梁元的剛度陣可以完全等價。

2)Bush 元對應單個節點的剛度陣并非對角陣,存在垂向位移和轉角的耦合剛度項,而且對角線元素也并非輸入參數k i(i =1, ...,6);由于存在耦合項,一般來說,Bush 元的剛度不可能由6 個標量彈簧元(分別對應節點6 個自由度的剛度)來代替,二者無法等效。

3)對于歐拉梁來說,剪切效應可忽略,此時,與其剛度等價的Bush 元只有4 個獨立的剛度系數。如果采用試驗方法測得6 個剛度系數,則6 個剛度系數中有兩對系數(k2與k6、k3與k5)可以相互校驗測量結果的準確性。

)

[1]陳烈民.航天器結構與機構[M].北京:中國科學技術出版社,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]張阿舟,諸德超,姚起杭, 等.實用振動工程(1)——振動理論與分析[M].北京:航空工業出版社, 1996

[6]邱吉寶, 向樹紅, 張正平.計算結構動力學[M].合肥:中國科學技術大學出版社, 2009

[7]李國強, 李進軍.Timoshenko-Euler 楔形梁有限元[J].力學季刊,2002, 23(1):64-69

[8]朱炳麒,陳學宏.理性Timoshenko 梁單元及其應用[J].力學與實踐, 2008, 30 (1):31-34

[9]王勖成, 邵敏.有限單元法基本原理和數值方法(第2版)[M].北京:清華大學出版社, 1997

[10]金咸定,趙德有.船體振動學[M].上海:上海交通大學出版社,1997