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

孫東陽(yáng)， 陳國(guó)平

(南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國(guó)家重點(diǎn)實(shí)驗(yàn)室， 江蘇 南京 210016)

引 言

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

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

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

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

1 絕對(duì)節(jié)點(diǎn)坐標(biāo)法

1.1 基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法的梁?jiǎn)卧?/h3>

絕對(duì)節(jié)點(diǎn)坐標(biāo)法中，柔性體k的一維梁?jiǎn)卧猨，如圖1所示，該單元上任意一點(diǎn)全局位置為

rkj=Skjekj

(1)

圖1 平面梁?jiǎn)卧?/p>

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

基于上述描述，梁?jiǎn)卧膭?dòng)能可表示為

(2)

(3)

式中εkj和κkj分別為單元j中線(xiàn)上對(duì)應(yīng)點(diǎn)的應(yīng)變和曲率。

基于虛功原理，可得到該單元的動(dòng)力學(xué)方程為

(4)

(5)

(6)

(7)

1.2 模型降階

(8)

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

ek=Tkqk

(9)

(10)

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

(11)

2 自然坐標(biāo)法

自然坐標(biāo)法中，用兩個(gè)基點(diǎn)的絕對(duì)坐標(biāo)矢量描述剛體i的位置和方向，如

(12)

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

(13)

圖2 兩基點(diǎn)剛體

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

(14)

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

(15)

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

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

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

3.1 旋轉(zhuǎn)副約束

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

(16)

(17)

式中

3.2 固結(jié)約束

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

(18)

(19)

式中

3.3 消除線(xiàn)性約束方程

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

(20)

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

(21)

由式(20)可以對(duì)廣義坐標(biāo)P進(jìn)行減縮

(22)

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

(23)

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

(24)

4 數(shù)值算例

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

圖4 雙擺

圖5 C點(diǎn)Y方向絕對(duì)位移

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

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

表1 計(jì)算所耗CPU時(shí)間

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

圖6 C點(diǎn)Y方向絕對(duì)位移

圖7 柔性梁端點(diǎn)橫向變形

5 結(jié) 論

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

參考文獻(xiàn)：

[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] 田強(qiáng)，張青云，陳立平，等. 柔性多體系統(tǒng)動(dòng)力學(xué)絕對(duì)節(jié)點(diǎn)坐標(biāo)方法研究進(jìn)展[J]. 力學(xué)進(jìn)展, 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] 王文亮，杜作潤(rùn)，陳康元. 模態(tài)綜合技術(shù)短評(píng)和一種新的改進(jìn)[J]. 航空學(xué)報(bào), 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] 李彬. 作大范圍空間運(yùn)動(dòng)的柔性梁的剛-柔耦合動(dòng)力學(xué)[D]. 上海: 上海交通大學(xué), 2005.Li Bing. Rigid-flexible coupling dynamics of a flexible beam with three-dimensional large overall motion[D]. Shanghai: Shanghai Jiao Tong University, 2005.