劉俊誼， 謝興坤， 洪 娟， 黃 恒， 計 淞

(陸軍工程大學(xué) 野戰(zhàn)工程學(xué)院，江蘇 南京 210007)

作為當(dāng)今力學(xué)領(lǐng)域的研究熱點和難點之一，多體系統(tǒng)動力學(xué)為航空航天、兵器、機(jī)械、機(jī)器人和近海工程等領(lǐng)域中相關(guān)系統(tǒng)的動態(tài)性能評估和優(yōu)化設(shè)計提供了重要的理論基礎(chǔ)和技術(shù)支撐[1-3]。20世紀(jì)70年代，Wittenburg[4]出版了第一本多體系統(tǒng)動力學(xué)專著，利用圖論方法奠定了多剛體系統(tǒng)動力學(xué)拉格朗日(Lagrange)方法的基礎(chǔ)。之后，Schichlen[5]基于有限元系統(tǒng)和連續(xù)系統(tǒng)相互等價的假定，提出了Schichlen方法。Kane等[6]建立了兼有矢量力學(xué)和分析力學(xué)優(yōu)點的凱恩(Kane)方法，并研究了該方法在航天器動力學(xué)上的應(yīng)用。在此基礎(chǔ)上，Huston[7]提出了適于數(shù)值計算的Kane-Huston方法。該方法采用低序體陣列描述多體系統(tǒng)的拓?fù)浣Y(jié)構(gòu)，利用歐拉(Euler)參數(shù)描述系統(tǒng)物體間的相對方位，避免了數(shù)值計算過程的奇異性困難。

Kane方程的形式較為簡單，非常適合創(chuàng)建復(fù)雜系統(tǒng)的動力學(xué)模型，因此很多學(xué)者將其應(yīng)用于空間機(jī)械臂[8]、機(jī)器人[9]及近海工程結(jié)構(gòu)物[10]等的動態(tài)響應(yīng)問題研究。Yang等[11-12]基于Kane動力學(xué)方程分別建立了水下蛇形機(jī)器人和水下四足行走機(jī)器人的動力學(xué)模型，并對它們的運動響應(yīng)進(jìn)行了數(shù)值分析。杜曉旭等[13]利用Kane方法建立了波浪驅(qū)動的水下航行器動力學(xué)方程，研究了纜索長度和航行深度對系統(tǒng)運動的影響。通過將纜索簡化為一系列具有不同物理特性的相互鉸接的剛性桿，李曉平等[14-15]基于Kane方法建立了水下拖曳柔索系統(tǒng)的三維有限段模型，并利用模型試驗驗證了計算模型的合理性。王磊等[16]將Kane方法與模態(tài)疊加法相結(jié)合，建立了陸上風(fēng)電機(jī)組整機(jī)剛?cè)狁詈辖Y(jié)構(gòu)動力學(xué)模型。此外，Liu等[17]利用Kane方法對走錨消能式攔阻系統(tǒng)受船舶撞擊后的動態(tài)響應(yīng)進(jìn)行了研究，并通過多體系統(tǒng)動力學(xué)軟件機(jī)械系統(tǒng)動力學(xué)自動分析(automatic dynamic analysis of mechanical systems，ADAMS)驗證了所建立的多體動力學(xué)模型的有效性及開發(fā)程序的可信性。

目前，基于Kane方法所開展的多體系統(tǒng)動力學(xué)分析中，主要側(cè)重于系統(tǒng)中各物體本身的運動和動力響應(yīng)的求解，而對于物體間的接點約束力的研究相對較少。對此，本文以鏈?zhǔn)蕉囿w系統(tǒng)為研究對象，基于Kane方法建立其動力特性分析的數(shù)值計算模型，給出其接點約束力的求解方法，并通過與多體系統(tǒng)離散時間傳遞矩陣法和經(jīng)典力學(xué)方法相關(guān)結(jié)果的對比分析，驗證所建模型的合理性和有效性。

1 鏈?zhǔn)蕉囿w系統(tǒng)動力特性分析的數(shù)學(xué)模型

Kane方法兼有矢量力學(xué)與分析力學(xué)的優(yōu)點，Huston方法是Kane方法在多剛體系統(tǒng)中的具體應(yīng)用，它采用了獨創(chuàng)的低序體陣列描述多體系統(tǒng)結(jié)構(gòu)，并引入了偏速度、偏角速度、廣義速率和廣義力等概念。運用Huston方法進(jìn)行多體系統(tǒng)動力學(xué)分析的核心是計算由偏速度、偏角速度及廣義速率等表達(dá)的廣義慣性力和廣義主動力，進(jìn)而建立Kane方程并進(jìn)行求解。下面按照Kane-Huston方法[7]的基本思想，建立鏈?zhǔn)蕉囿w系統(tǒng)動力特性分析的數(shù)學(xué)模型。

1.1 坐標(biāo)系的建立及坐標(biāo)轉(zhuǎn)換矩陣

首先，對鏈?zhǔn)蕉囿w系統(tǒng)(N級擺)中的各物體進(jìn)行編號并建立如圖1所示坐標(biāo)系，與固定端o0相連的桿件記為B1，然后沿其連接方向依次編號為B2,…,BN。慣性系的原點位于系統(tǒng)的固定端o0處，x0軸水平向右，z0軸豎直向上，y0軸與x0軸、z0軸滿足右手定則。剛體Bk與其內(nèi)接剛體(鄰接較低序號的剛體)的鉸接點即是Bk的連體坐標(biāo)系原點ok，xk軸沿桿件方向，yk軸與y0軸方向相同，zk軸與xk軸、yk軸同樣滿足右手定則。圖中θk表示剛體Bk與x0軸正向的夾角。

圖1 鏈?zhǔn)蕉囿w系統(tǒng)示意圖

物體偏速度、偏角速度的計算與坐標(biāo)系間的變換矩陣密切相關(guān)。因此，首先需要明確變換矩陣的相關(guān)計算。當(dāng)使用剛體Bk相對其內(nèi)接剛體Bj運動的卡爾丹角αk、βk、γk描述剛體姿態(tài)時，剛體Bk相對其內(nèi)接剛體Bj的坐標(biāo)變換矩陣可寫為

式中：C表示余弦運算cos，S表示正弦運算sin。為了避免使用卡爾丹角進(jìn)行數(shù)值計算時可能遇到的奇異性問題，采用歐拉參數(shù)描述剛體轉(zhuǎn)動。由變換矩陣求解歐拉參數(shù)εki(i=1,2,3,4)的公式為

S(jk)可用歐拉參數(shù)εki(i=1,2,3,4)寫為

初始時刻，各物體間的夾角已知，因而可以通過式(1)求得變換矩陣S(jk)。然后，可以利用式(2)求解此時刻相應(yīng)的歐拉參數(shù)，之后的數(shù)值計算中，每一時刻的歐拉參數(shù)均可通過求解式(4)的一階微分方程得到。

(4)

式中：ωk1、ωk2、ωk3為剛體Bk相對其內(nèi)接剛體的角速度分量。將歐拉參數(shù)的計算結(jié)果代入式(3)即可得到相應(yīng)的剛體間變換矩陣。剛體Bk相對慣性坐標(biāo)系的變換矩陣S(0k)可由Bk至零剛體(慣性系)通路上各對連接剛體間的變換矩陣連乘得到，即

v=Lt(k),w=Lt+1(k),Lu(k)=1

(5)

L(k)為鄰接較低序號算子，在此基礎(chǔ)上求解剛體Bk相對于慣性系轉(zhuǎn)動運動響應(yīng)的公式為

式中：α(0k)、β(0k)、γ(0k)分別為剛體Bk相對慣性系運動的卡爾丹角。

變換矩陣的導(dǎo)數(shù)可通過矩陣的乘積表示為

(7)

式中：W(0k)為剛體Bk角速度ω(0k)的對偶矩陣。其表達(dá)式為

i,r,m=1,2,3

(9)

1.2 廣義速率和廣義坐標(biāo)

(10)

1.3 偏速度、偏角速度陣列及其導(dǎo)數(shù)

點的速度和體的角速度可表示為廣義坐標(biāo)導(dǎo)數(shù)的線性函數(shù)組合，這些導(dǎo)數(shù)的系數(shù)稱為“偏速度”和“偏角速度”。剛體Bk在慣性參考系中的絕對角速度ω(0k)可用偏角速度表示為

k=1,2,…,N;l=1,2,…,6N;m=1,2,3

(11)

式中：m、l為啞標(biāo)；n0m為慣性坐標(biāo)系的矢量基。根據(jù)廣義速率的定義，偏角速度ωklm的表達(dá)式為

式中δm，l表示Kronecker’s delta函數(shù)

(13)

如圖2所示，令ok為剛體Bk連體坐標(biāo)系的原點，Qk為剛體Bk與其內(nèi)接剛體Bj的連接點并固定在Bj上。矢量qk給定連接點Qk與原點oj(j=L(k))的相對位置，用剛體Bj的連體坐標(biāo)系表示。矢量sk表示原點ok與連接點Qk的相對位置，描述剛體Bk相對其內(nèi)接剛體Bj的移動，用剛體Bj的連體坐標(biāo)系表示。矢量rk給定質(zhì)心Ck與原點ok的相對位置，用剛體Bk的連體坐標(biāo)系表示。

圖2 通路上的物體和相關(guān)矢量

同理，剛體Bk的絕對速度v(0k)可用偏速度表示為

k=1,2,…,N;l=1,2,…,6N;m=1,2,3

(15)

偏速度分量vklm的表達(dá)式為

式中：qvn和svn分別為通路上連接點矢量qv和相對位移矢量sv在Bv內(nèi)接剛體連體坐標(biāo)系中的投影值；rkn為rk在剛體Bk連體坐標(biāo)系中的投影值。

v=Lt(k),w=Lt+1(k),Lu(k)=1,

l=1,2,…,3k;k=1,2,…,N;t,h,m=1,2,3

l=3k+1,3k+2,…,3N;k=1,2,…,N

l=3N+1,3N+2,…,6N;k=1,2,…,N

(17)

1.4 鏈?zhǔn)蕉囿w系統(tǒng)的凱恩方程

通過前面的分析可以發(fā)現(xiàn)：任意剛體Bk質(zhì)心的速度v(0k)和加速度ak以及Bk的角速度ω(0k)和角加速度αk可用偏速度和偏角速度陣列及其導(dǎo)數(shù)表示為

若剛體Bk對于原點位于其質(zhì)心處連體坐標(biāo)系的慣量矩陣為[ICk]，則Bk相對于慣性參考系的慣量矩陣為

[Ik]=[S(0k)]T[ICk][S(0k)]

(19)

(21)

Fl=vklmFkm+ωklmTkm

l=1,2,…,6N;k=1,2,…,N;m=1,2,3

(23)

系統(tǒng)內(nèi)部光滑鉸和光滑面接觸等理想約束成對出現(xiàn)，其約束力對廣義主動力的貢獻(xiàn)為零，因而可不予考慮。

(24)

將廣義力的表達(dá)式代入式(24)，同時考慮到鏈?zhǔn)蕉囿w系統(tǒng)物體間無相對移動，化簡后可得到鏈?zhǔn)蕉囿w系統(tǒng)的基本運動力學(xué)方程

(25)

式中的系數(shù)為

2 接點約束力的求解方法

根據(jù)上述的Huston方法，可以求得每一時刻每個物體相對靜止坐標(biāo)系的運動響應(yīng)，包括角速度、角加速度、速度、加速度等物理量，然而，計算廣義主動力時并沒有考慮光滑鉸處的約束力，因此若要確定它們的大小則需進(jìn)一步地分析計算。下面分別基于Kane方法和經(jīng)典力學(xué)的方法對其進(jìn)行分析。

2.1 基于Kane方法求解

首先，根據(jù)前面分析中已經(jīng)消去的方程計算廣義約束力(作用于各剛體質(zhì)心處約束力主矢和主矩的廣義力)。其過程如下：

將接點約束暫時去掉，而保留約束力和力矩分量，且視之為外部作用的力和力矩。系統(tǒng)將具有與被去掉的約束相應(yīng)的增加的自由度，也就出現(xiàn)了與這些增加的自由度相應(yīng)的增加的動力學(xué)方程(即原來被消去的那些方程)。待求的廣義約束力將出現(xiàn)在這些方程中，而且，由于采用體間的相對角速度分量或相對位移速率作為廣義速率，廣義約束力將在這些方程中單獨(無耦合)出現(xiàn)。即在每個方程中將僅有廣義約束力的1個分量(如用方位角導(dǎo)數(shù)作為廣義速率，將不是這種情況)。將已知的廣義速率和廣義坐標(biāo)代入這些方程，則各方程即成為求解廣義約束力的解耦的代數(shù)方程。

鏈?zhǔn)蕉囿w系統(tǒng)中，暫時去掉接點約束可得方程為

l=3N+1,3N+2,…,6N;p=1,2,…,3N

(27)

(28)

由于假定鏈?zhǔn)蕉囿w系統(tǒng)中的鉸鏈均為光滑鉸，所以鉸鏈處只有約束力而不存在約束力矩(阻尼力矩)的作用。因此，系統(tǒng)的廣義約束力為

(29)

l=1,2,…,6N;k=1,2,…,N;m=1,2,3

(30)

(31)

2.2 基于經(jīng)典力學(xué)的方法求解

對鏈?zhǔn)蕉囿w系統(tǒng)中任意剛體Bk進(jìn)行力學(xué)分析如圖3所示。

圖3 N級擺系統(tǒng)內(nèi)桿件受力示意圖

根據(jù)牛頓第二定律和第三定律可得

(32)

(33)

3 算例分析

3.1 參數(shù)選擇

取N級擺節(jié)數(shù)N=3；重力加速度g=9.81 m/s2；積分時間步長Δt=0.001 s；積分時長t=10 s；各桿件長度l=1 m；半徑R=0.1 m；質(zhì)量m=1 kg；由于各桿件為均質(zhì)桿件，根據(jù)上述對連接矢量的定義，q1=[0,0,0]T、qk=[l,0,0]T(k=2,3,…,N)、rk=[l/2,0,0]T(k=1,2,…,N)。初始時刻各桿件靜止且與x0軸正向夾角均為-π/3，各擺僅在x0o0z0平面內(nèi)運動，則有αk=γk=0(k=1,2,3)、β1=π/3,β2=β3=0；各桿件在以其質(zhì)心為原點的連體坐標(biāo)系中的慣量矩陣為

(34)

3.2 結(jié)果分析

選取參數(shù)后，按照Huston方法的求解過程編程計算，即可得到系統(tǒng)內(nèi)各擺的運動響應(yīng)。下面首先確立數(shù)值計算的有效性，基于能量守恒原理的驗證及與多體系統(tǒng)離散時間傳遞矩陣法計算結(jié)果[18]的對比分別如圖4和圖5所示。

圖4 能量守恒原理對平面三級擺系統(tǒng)數(shù)值仿真精確性的驗證

圖5 平面三級擺的角位移計算結(jié)果

對比表明，兩種方法的計算結(jié)果一致性很好，說明了Huston方法解決多剛體動力學(xué)問題的有效性及數(shù)值計算的可信性。對兩種確定多體系統(tǒng)內(nèi)約束力的方法進(jìn)行比較，結(jié)果分別如圖6和圖7所示。

觀察圖6和圖7可以發(fā)現(xiàn)：兩種方法計算接點約束力的結(jié)果十分吻合，相互驗證了它們各自的有效性和準(zhǔn)確性。兩種方法的效率和精度均基本一致，但從計算便捷性的角度看，Kane方法略有優(yōu)勢，因為其不需要對每個物體進(jìn)行單獨的受力分析。

圖6 平面三級擺各鉸接點x方向約束力

圖7 平面三級擺各鉸接點z方向約束力

4 結(jié)論

本文首先運用Kane-Huston方法對鏈?zhǔn)蕉囿w系統(tǒng)的運動響應(yīng)進(jìn)行了建模，并提出兩種確定系統(tǒng)內(nèi)(鉸)接點約束力方法。然后，以三級擺為算例對多體系統(tǒng)的運動響應(yīng)及(鉸)接點約束力進(jìn)行了求解。通過能量守恒原理及與文獻(xiàn)[18]中的計算結(jié)果對比驗證了Huston方法數(shù)值計算的有效性，進(jìn)而將基于Kane方法和經(jīng)典力學(xué)方法求解(鉸)接點約束力的結(jié)果對比，驗證了兩種方法的有效性和準(zhǔn)確性。

本文所建模型及相關(guān)結(jié)果對利用Kane方法開展多體系統(tǒng)動力學(xué)分析具有一定的參考意義。但同時，本文的研究也存在一些不足。例如，算例中的鉸鏈假定為光滑鉸、系統(tǒng)內(nèi)各物體假定無初速度且不受除重力以外的外力作用。進(jìn)一步針對更為復(fù)雜的工況條件開展鏈?zhǔn)蕉囿w系統(tǒng)動力學(xué)分析具有重要的工程實踐意義。