劉俊誼， 謝興坤， 洪 娟， 黃 恒， 計(jì) 淞

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

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

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

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

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

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

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

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

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

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

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

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

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

(4)

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

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

(5)

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

式中：α(0k)、β(0k)、γ(0k)分別為剛體Bk相對(duì)慣性系運(yùn)動(dòng)的卡爾丹角。

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

(7)

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

i,r,m=1,2,3

(9)

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

(10)

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

點(diǎn)的速度和體的角速度可表示為廣義坐標(biāo)導(dǎo)數(shù)的線性函數(shù)組合，這些導(dǎo)數(shù)的系數(shù)稱為“偏速度”和“偏角速度”。剛體Bk在慣性參考系中的絕對(duì)角速度ω(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)系的原點(diǎn)，Qk為剛體Bk與其內(nèi)接剛體Bj的連接點(diǎn)并固定在Bj上。矢量qk給定連接點(diǎn)Qk與原點(diǎn)oj(j=L(k))的相對(duì)位置，用剛體Bj的連體坐標(biāo)系表示。矢量sk表示原點(diǎn)ok與連接點(diǎn)Qk的相對(duì)位置，描述剛體Bk相對(duì)其內(nèi)接剛體Bj的移動(dòng)，用剛體Bj的連體坐標(biāo)系表示。矢量rk給定質(zhì)心Ck與原點(diǎn)ok的相對(duì)位置，用剛體Bk的連體坐標(biāo)系表示。

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

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

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

(15)

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

式中：qvn和svn分別為通路上連接點(diǎn)矢量qv和相對(duì)位移矢量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對(duì)于原點(diǎn)位于其質(zhì)心處連體坐標(biāo)系的慣量矩陣為[ICk]，則Bk相對(duì)于慣性參考系的慣量矩陣為

[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)部光滑鉸和光滑面接觸等理想約束成對(duì)出現(xiàn)，其約束力對(duì)廣義主動(dòng)力的貢獻(xiàn)為零，因而可不予考慮。

(24)

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

(25)

式中的系數(shù)為

2 接點(diǎn)約束力的求解方法

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

2.1 基于Kane方法求解

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

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

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

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é)的方法求解

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

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

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

(32)

(33)

3 算例分析

3.1 參數(shù)選擇

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

(34)

3.2 結(jié)果分析

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

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

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

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

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

圖6 平面三級(jí)擺各鉸接點(diǎn)x方向約束力

圖7 平面三級(jí)擺各鉸接點(diǎn)z方向約束力

4 結(jié)論

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

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