剛?cè)狁詈隙囿w系統(tǒng)動力學模型降階

孫東陽， 陳國平

(南京航空航天大學機械結(jié)構(gòu)力學及控制國家重點實驗室， 江蘇 南京 210016)

引 言

隨著空間技術(shù)的發(fā)展，航天器柔性構(gòu)件的尺寸越來越大，柔性構(gòu)件的變形對航天器的動力學行為產(chǎn)生很大的影響，傳統(tǒng)多剛體系統(tǒng)動力學理論已經(jīng)不能滿足人們對設(shè)備精度的要求。最近幾十年，考慮構(gòu)件柔性的柔性多體系統(tǒng)動力學，已經(jīng)成為國內(nèi)外眾多學者研究的熱點并取得了大量的研究成果。考慮構(gòu)件彈性變形與大范圍剛體運動的耦合，Likins提出了浮動坐標方法[1]，該方法將構(gòu)件的位形認為是浮動坐標系大范圍剛體運動與相對于該局部坐標系的變形的疊加。1987年，Kane在研究梁的高速旋轉(zhuǎn)運動時第一次發(fā)現(xiàn)了動力剛化現(xiàn)象[2]。為解決動力剛化問題，Haering等在多體系統(tǒng)動力學建模過程中考慮了高階項[3,4]。

然而，浮動節(jié)點坐標法是基于小變形假設(shè)，對于存在大變形的多體系統(tǒng)已經(jīng)不再適用。Shabana提出了目前廣泛應(yīng)用于分析大變形柔性多體系統(tǒng)動力學的絕對節(jié)點坐標方法(ANCF)[5]。該方法中單元節(jié)點的坐標定義在全局坐標系下, 采用斜率矢量代替?zhèn)鹘y(tǒng)有限單元中的節(jié)點轉(zhuǎn)角坐標，推導(dǎo)建立的多體系統(tǒng)微分-代數(shù)方程的質(zhì)量矩陣為常數(shù)矩陣，且具有不存在科氏力和離心力項的優(yōu)點。Berzeri等提出了幾種基于不同假設(shè)的一維梁的彈性力簡化模型[6]，并做了比較研究。Escalona等首先將該單元應(yīng)用于柔性大變形多體系統(tǒng)動力學的研究[7]。Omar等放寬梁中線的切線矢量與梁截面的法線方向重合的的假設(shè)[8]，將梁的剪切變形考慮到梁單元中,首次提出了一種平面應(yīng)變剪變梁單元。該單元由于彎曲應(yīng)變與軸向應(yīng)變不一致而帶來剪變閉鎖問題。為了解決剪變閉鎖問題，Kerkk?nen等通過改變單元的運動學描述[9,10],提出了一些可有效減輕剪變閉鎖問題的線性剪變梁單元。考慮到絕對節(jié)點坐標體系下剛度矩陣的強非線性，導(dǎo)致采用絕對坐標方法建立的微分-代數(shù)方程的計算效率比較低, García-Vallejo提出不變矩陣法[11], 該方法將非線性剛度矩陣分解為常數(shù)剛度矩陣與廣義坐標相關(guān)的剛度矩陣之和。

絕對節(jié)點坐標法雖然能夠精確描述多體系統(tǒng)的剛性運動，但是即使是剛體也要劃分單元[12]，這必然導(dǎo)致該方法較難處理剛?cè)狁詈隙囿w系統(tǒng)動力學問題。Sugiyama等結(jié)合浮動節(jié)點坐標法能有效處理剛體運動和絕對節(jié)點坐標法能描述柔性體大變形的特點[13～15]，使存在柔性體大變形的剛?cè)狁詈隙囿w系統(tǒng)得到了解決。然而，浮動節(jié)點坐標法建立的多體系統(tǒng)動力學方程的質(zhì)量矩陣與廣義坐標相關(guān)，因而得到的剛?cè)狁詈隙囿w系統(tǒng)的微分-代數(shù)方程的質(zhì)量矩陣也與廣義坐標相關(guān)，每次計算都需要對質(zhì)量陣進行計算，會大大降低計算效率。自然坐標法以其具有建立的多體系統(tǒng)微分-代數(shù)方程的質(zhì)量矩陣為常數(shù)矩陣的優(yōu)點而成為剛?cè)狁詈隙囿w系統(tǒng)中替代浮動節(jié)點坐標法用于描述剛體構(gòu)件的方法。García-Vallejo用自然坐標法與絕對節(jié)點坐標法的混合方法對平面剛?cè)狁詈隙囿w系統(tǒng)動力學進行了研究[16]，并且對構(gòu)件的各種連接情況進行了討論。在此基礎(chǔ)上，García-Vallejo進一步對空間剛?cè)狁詈隙囿w系統(tǒng)動力學問題進行了研究[17]。

上述研究隨著柔性體單元數(shù)量的增加，剛?cè)狁詈隙囿w系統(tǒng)動力學方程的計算效率將會比較低。為了提高計算效率，需要對絕對節(jié)點坐標法建立的柔性多體系統(tǒng)進行模型降階。傳統(tǒng)的模型降階方法主要有子結(jié)構(gòu)方法[18～21]，Krylov子空間方法等[22,23]。Hurty首先提出了模態(tài)綜合法的概念[18]，并且應(yīng)用于對大規(guī)模線性系統(tǒng)進行模型降階。R R Craig和C C Bampton對此方法作了部分修正[19]，形成了現(xiàn)在的固定界面模態(tài)綜合法。隨后，Kobayashi成功將Craig-Bampton固定界面法應(yīng)用于基于絕對節(jié)點坐標法的柔性多體系統(tǒng)的模型降階[24]。本文采用Craig-Bampton固定界面法對基于絕對節(jié)點坐標法和自然坐標法建立的剛?cè)狁詈隙囿w系統(tǒng)進行模型降價。

1 絕對節(jié)點坐標法

1.1 基于絕對節(jié)點坐標法的梁單元

絕對節(jié)點坐標法中，柔性體k的一維梁單元j，如圖1所示，該單元上任意一點全局位置為

rkj=Skjekj

(1)

圖1 平面梁單元

式中Skj為單元形函數(shù)，ekj為單元節(jié)點的廣義坐標矢量，可表示為

基于上述描述，梁單元的動能可表示為

(2)

(3)

式中εkj和κkj分別為單元j中線上對應(yīng)點的應(yīng)變和曲率。

基于虛功原理，可得到該單元的動力學方程為

(4)

(5)

(6)

(7)

1.2 模型降階

(8)

用Craig-Bampton方法進行減縮時，約束模態(tài)僅取前nc階，則廣義坐標ek可減縮為

ek=Tkqk

(9)

(10)

將式(9)代入式(8)，等式兩邊左乘TkT，則有

(11)

2 自然坐標法

自然坐標法中，用兩個基點的絕對坐標矢量描述剛體i的位置和方向，如

(12)

式中di為包含基點C和D的剛體i的坐標矢量(如圖2所示)，該坐標矢量是非獨立的，C和D之間存在距離約束，有

(13)

圖2 兩基點剛體

基于上述剛體描述，自然坐標法建立的剛體系統(tǒng)動力學方程為

(14)

固結(jié)在剛體上的動坐標系的旋轉(zhuǎn)矩陣可表示為

(15)

3 剛?cè)狁詈隙囿w系統(tǒng)

柔性體之間、剛體之間、柔性體與剛體之間存在各種約束。本文僅描述剛體與柔性體之間的旋轉(zhuǎn)副約束和固結(jié)約束(如圖3所示)，其他相關(guān)約束可參考文獻[16]。

圖3 柔性體與剛體之間的約束

3.1 旋轉(zhuǎn)副約束

假設(shè)剛體i的P點與柔性體k的節(jié)點n存在旋轉(zhuǎn)副約束(如圖3(a)所示)，約束方程可表示為

(16)

(17)

式中

3.2 固結(jié)約束

假設(shè)剛體i的P點與柔性體k的節(jié)點n存在固結(jié)約束(如圖3(b)所示)，則約束點除了存在位置約束外還存在方向約束，即

(18)

(19)

式中

3.3 消除線性約束方程

將約束節(jié)點作為柔性體的邊界，則柔性體的邊界節(jié)點廣義坐標可以用連接剛體的自然坐標和柔性體邊界的減縮坐標表示。假設(shè)剛體i的P點與柔性體k的節(jié)點n存在旋轉(zhuǎn)副約束，則柔性體k經(jīng)Craig-Bampton方法減縮后的廣義坐標可表示為

(20)

由式(11)和(14)得剛?cè)狁詈隙囿w系統(tǒng)動力學方程為

(21)

由式(20)可以對廣義坐標P進行減縮

(22)

將式(22)代入式(21)，等式兩邊左乘TT，則

(23)

假設(shè)剛體i的P點與柔性體k的節(jié)點n存在固結(jié)約束，同理可得消除邊界約束后的剛?cè)狁詈隙囿w系統(tǒng)動力學方程為

(24)

4 數(shù)值算例

圖4為受重力的平面雙擺。相關(guān)參數(shù)參考文獻[25],A點與B點均為旋轉(zhuǎn)鉸，梁AB與梁BC的長度均為1.8 m，截面積為2.5×10-4m2，密度為2.766 67×103kg/m3。梁AB與梁BC都可以作為剛體或柔性體。本文首先基于3種情況對該雙擺系統(tǒng)進行動力學分析：(1)梁AB與梁BC均為柔性體，用絕對節(jié)點坐標法對其進行研究；(2)梁AB與梁BC均為剛體，用自然坐標法對其進行研究(NCF)；(3)梁AB為柔性體，梁BC為剛體，用文獻[16]的平面剛?cè)狁詈隙囿w系統(tǒng)動力學方法對其進行研究(NCF-ANCF)。如果梁為柔性體，將梁等分為10個單元，楊氏模量為6.895×109Pa，截面慣量矩為1.302×10-9m4，取重力加速度為9.81 m/s2，仿真時間為2.5 s。在一臺具有Intel Pentium 3.2 GHz處理器及3GB RAM的PC機上運行。3種情況下，C點Y方向的絕對位移如圖5所示

圖4 雙擺

圖5 C點Y方向絕對位移

由圖5可知，3種情況下C點Y方向的位移存在差異，這說明彈性變形對雙擺端點C的運動會產(chǎn)生影響，但3種情況下C點位移相差不大，因此，在精度要求不高的情況下，可以把梁作為剛體考慮。

多體系統(tǒng)按照構(gòu)件在運行過程中的變形，可以將構(gòu)件分為剛體構(gòu)件和柔性體構(gòu)件，其中，隨著柔性體構(gòu)件劃分單元數(shù)量的增加，剛?cè)狁詈隙囿w系統(tǒng)動力學方程會存在計算時間長，計算效率低的問題。因此，本文將模型降階的方法引入剛?cè)狁詈隙囿w系統(tǒng)動力學，提出基于模態(tài)綜合法的剛?cè)狁詈隙囿w系統(tǒng)減縮方法，解決了傳統(tǒng)存在大變形的剛?cè)狁詈隙囿w系統(tǒng)動力學方程計算效率低的問題。為了驗證該方法的正確性和有效性，本文以圖4的雙擺為例，將梁AB作為柔性體，取其楊氏模量為6.895×108Pa，梁BC作為剛體，用該方法對此系統(tǒng)進行研究。在選取主模態(tài)集時，分別取前5階模態(tài)(C-N-A(5))，前7階模態(tài)(C-N-A(7))和前9階模態(tài)(C-N-A(9))，將其與原模型(N-A)進行對比。減縮模型和原模型C點Y方向位移和柔性梁端點橫向變形分別如圖6，7所示,計算所耗CPU時間如表1所示。

表1 計算所耗CPU時間

由圖6可知，取較少的主模態(tài)就能使C點Y方向絕對位移誤差很小，這是因為柔性梁的振動主要為低頻振動，因此，只用選取較少的模態(tài)就能粗略描述柔性梁的彈性變形。圖7中，原模型的柔性梁端點橫向變形曲線與取前9階模態(tài)時的柔性梁端點橫向變形曲線幾乎重合，而與取前5階、前7階模態(tài)時的柔性梁端點橫向變形曲線存在明顯差異。由此可見，只有適當選取更多模態(tài)才能更好地描述柔性體的彈性變形。由表1可知，減縮模型的計算時間都比原模型的計算時間要少，而且選擇的模態(tài)數(shù)量越少越節(jié)省計算時間。結(jié)合圖7和表1可以得到：減縮模型僅選取前9階模態(tài)時，能夠在滿足計算精度的情況下使計算時間僅為原模型計算時間的34.9%。由上述分析可知，適當選取模態(tài)就可以在保證計算精度的情況下，減少剛?cè)狁詈隙囿w系統(tǒng)動力學方程的計算時間，從而提高計算效率。

圖6 C點Y方向絕對位移

圖7 柔性梁端點橫向變形

5 結(jié) 論

本文考慮到存在大變形的剛?cè)狁詈隙囿w系統(tǒng)動力學方程計算效率比較低，提出了基于模態(tài)綜合法的剛?cè)狁詈隙囿w系統(tǒng)模型降階方法。通過雙擺系統(tǒng)對該方法的正確性和有效性進行了研究。由數(shù)值仿真可知，減縮模型隨著選擇模態(tài)數(shù)量的減少，同原模型相比越節(jié)省計算時間，而且只要適當選擇減縮模態(tài)就可以保證計算精度。這說明該方法能夠提高剛?cè)狁詈隙囿w系統(tǒng)動力學方程的計算效率。

參考文獻：

[1] Likins P W. Finite element appendage equations for hybrid coordinate dynamic analysis[J]. Journal of Solids and Structures, 1972, 8: 709—731.

[2] Kane T R, Ryan R R, Banerjee A K. Dynamics of a cantilever beam attached to a moving base[J]. Journal of Guidance, Control, and Dynamics, 1987, 10(2): 139—151.

[3] Haering W J, Ryan R R, Scott R A. A new flexible body dynamic formulation for beam structures undergoing large overall motion[A]. 33rd Structural Dynamics and Materials Conference[C]. Dallas: AIAA, 19921415-1423.

[4] Mayo J M, García-Vallejo D, Domínguez J. Study of the geometric stiffening effect: comparison of different formulations[J]. Multibody System Dynamics, 2004, 11: 321—341.

[5] Shabana A A. Definition of the slopes and the finite element absolute nodal coordinate formulation[J]. Multibody System Dynamics, 1997, 1: 339—348.

[6] Berzeri M, Shabana A A. Development of simple models for the elastic forces in the absolute nodal coordinate formulation[J]. Journal of Sound and Vibration, 2000, 235(4): 539—565.

[7] Escalona J L, Hussien H A, Shabana A A. Application of the absolute nodal coordinate formulation to multibody system dynamics[J]. Journal of Sound and Vibration, 1998, 214(5): 833—851.

[8] Omar M A, Shabana A A. A two-dimensional shear deformable beam for large rotation and deformation problems[J]. Journal of Sound and Vibration, 2001, 243(3): 565—576.

[9] Kerkk?nen K S, Sopanen J T, Mikkola A M. A linear beam finite element based on the absolute nodal coordinate formulation[J]. ASME Journal of Mechanical Design, 2005, 127: 621—630.

[10] García-Vallejo D, Mikkola A M, Escalona J L. A new locking-free shear deformable finite element based on absolute nodal coordinates[J]. Nonlinear Dynamics, 2007, 50: 249—264.

[11] García-Vallejo D, Mayo J, Escalona J L, et al. Efficient evaluation of the elastic forces and the Jacobian in the absolute nodal coordinate formulation[J]. Nonlinear Dynamics, 2004, 35: 313—329.

[12] 田強，張青云，陳立平，等. 柔性多體系統(tǒng)動力學絕對節(jié)點坐標方法研究進展[J]. 力學進展, 2010, 40(2): 189—202.TIAN Qiang, ZHANG Yunqing, CHEN Liping, et al. Advances in the absolute nodal coordinate method for the flexible multibody dynamics[J]. Advances in Mechanics, 2010, 40(2): 189—202.

[13] Lee S, Park T, Seo J, et al. The development of a sliding joint for very flexible multibody dynamics using absolute nodal coordinate formulation[J]. Multibody System Dynamics, 2008, 20: 223—237.

[14] Hussein B A, Weed D, Shabana A A. Clamped end conditions and cross section deformation in the finite element absolute nodal coordinate formulation[J]. Multibody System Dynamics, 2009, 21: 375—393.

[15] Sugiyama H, Yamashita H. Spatial joint constraints for the absolute nodal coordinate formulation using the non-generalized intermediate coordinates[J]. Multibody System Dynamics, 2011, 26: 15—36.

[16] García-Vallejo D, Escalona J L, Mayo J, et al. Describing rigid-flexible multibody systems using absolute coordinates[J]. Nonlinear Dynamics, 2003, 34: 75—94.

[17] García-Vallejo D, Mayo J, Escalona J L, et al. Three-dimensional formulation of rigid-flexible multibody systems with flexible beam elements [J]. Multibody System Dynamics, 2008, 20: 1—28.

[18] Hurty W C. Dynamic analysis of structural systems using component modes[J]. AIAA Journal, 1965, 3(4): 678—685.

[19] Craig R R, Bampton M C C. Coupling of substructures for dynamic analysis[J]. AIAA Journal, 1968, 6(7): 1 313—1 319.

[20] Hou S N. Review of modal synthesis techniques and a new approach[J]. Shock and Vibration Bulletin, 1969, 40(4): 25—30.

[21] 王文亮，杜作潤，陳康元. 模態(tài)綜合技術(shù)短評和一種新的改進[J]. 航空學報, 1979, 3(2): 32—51.

[22] Antoulas A C. Approximation of Large-scale Dynamical Systems[M]. Philadelphia: SIAM, 2004.

[23] Bai Z. Krylov subspace techniques for reduced-order modeling of large-scale dynamical systems[J]. Applied Numerical Mathematics, 2002, 43: 9—44.

[24] Kobayashi N, Wago T, Sugawara Y. Reduction of system matrices of planar beam in ANCF by component mode synthesis method[J]. Multibody System Dynamics, 2011, 26: 265—281.

[25] 李彬. 作大范圍空間運動的柔性梁的剛-柔耦合動力學[D]. 上海: 上海交通大學, 2005.Li Bing. Rigid-flexible coupling dynamics of a flexible beam with three-dimensional large overall motion[D]. Shanghai: Shanghai Jiao Tong University, 2005.