一種鏈式多體系統接點約束力的求解方法

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

(陸軍工程大學 野戰工程學院，江蘇 南京 210007)

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

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

目前，基于Kane方法所開展的多體系統動力學分析中，主要側重于系統中各物體本身的運動和動力響應的求解，而對于物體間的接點約束力的研究相對較少。對此，本文以鏈式多體系統為研究對象，基于Kane方法建立其動力特性分析的數值計算模型，給出其接點約束力的求解方法，并通過與多體系統離散時間傳遞矩陣法和經典力學方法相關結果的對比分析，驗證所建模型的合理性和有效性。

1 鏈式多體系統動力特性分析的數學模型

Kane方法兼有矢量力學與分析力學的優點，Huston方法是Kane方法在多剛體系統中的具體應用，它采用了獨創的低序體陣列描述多體系統結構，并引入了偏速度、偏角速度、廣義速率和廣義力等概念。運用Huston方法進行多體系統動力學分析的核心是計算由偏速度、偏角速度及廣義速率等表達的廣義慣性力和廣義主動力，進而建立Kane方程并進行求解。下面按照Kane-Huston方法[7]的基本思想，建立鏈式多體系統動力特性分析的數學模型。

1.1 坐標系的建立及坐標轉換矩陣

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

圖1 鏈式多體系統示意圖

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

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

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

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

(4)

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

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

(5)

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

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

變換矩陣的導數可通過矩陣的乘積表示為

(7)

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

i,r,m=1,2,3

(9)

1.2 廣義速率和廣義坐標

(10)

1.3 偏速度、偏角速度陣列及其導數

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

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

(11)

式中：m、l為啞標；n0m為慣性坐標系的矢量基。根據廣義速率的定義，偏角速度ωklm的表達式為

式中δm，l表示Kronecker’s delta函數

(13)

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

圖2 通路上的物體和相關矢量

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

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

(15)

偏速度分量vklm的表達式為

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

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 鏈式多體系統的凱恩方程

通過前面的分析可以發現：任意剛體Bk質心的速度v(0k)和加速度ak以及Bk的角速度ω(0k)和角加速度αk可用偏速度和偏角速度陣列及其導數表示為

若剛體Bk對于原點位于其質心處連體坐標系的慣量矩陣為[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)

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

(24)

將廣義力的表達式代入式(24)，同時考慮到鏈式多體系統物體間無相對移動，化簡后可得到鏈式多體系統的基本運動力學方程

(25)

式中的系數為

2 接點約束力的求解方法

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

2.1 基于Kane方法求解

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

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

鏈式多體系統中，暫時去掉接點約束可得方程為

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

(27)

(28)

由于假定鏈式多體系統中的鉸鏈均為光滑鉸，所以鉸鏈處只有約束力而不存在約束力矩(阻尼力矩)的作用。因此，系統的廣義約束力為

(29)

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

(30)

(31)

2.2 基于經典力學的方法求解

對鏈式多體系統中任意剛體Bk進行力學分析如圖3所示。

圖3 N級擺系統內桿件受力示意圖

根據牛頓第二定律和第三定律可得

(32)

(33)

3 算例分析

3.1 參數選擇

取N級擺節數N=3；重力加速度g=9.81 m/s2；積分時間步長Δt=0.001 s；積分時長t=10 s；各桿件長度l=1 m；半徑R=0.1 m；質量m=1 kg；由于各桿件為均質桿件，根據上述對連接矢量的定義，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平面內運動，則有αk=γk=0(k=1,2,3)、β1=π/3,β2=β3=0；各桿件在以其質心為原點的連體坐標系中的慣量矩陣為

(34)

3.2 結果分析

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

圖4 能量守恒原理對平面三級擺系統數值仿真精確性的驗證

圖5 平面三級擺的角位移計算結果

對比表明，兩種方法的計算結果一致性很好，說明了Huston方法解決多剛體動力學問題的有效性及數值計算的可信性。對兩種確定多體系統內約束力的方法進行比較，結果分別如圖6和圖7所示。

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

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

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

4 結論

本文首先運用Kane-Huston方法對鏈式多體系統的運動響應進行了建模，并提出兩種確定系統內(鉸)接點約束力方法。然后，以三級擺為算例對多體系統的運動響應及(鉸)接點約束力進行了求解。通過能量守恒原理及與文獻[18]中的計算結果對比驗證了Huston方法數值計算的有效性，進而將基于Kane方法和經典力學方法求解(鉸)接點約束力的結果對比，驗證了兩種方法的有效性和準確性。

本文所建模型及相關結果對利用Kane方法開展多體系統動力學分析具有一定的參考意義。但同時，本文的研究也存在一些不足。例如，算例中的鉸鏈假定為光滑鉸、系統內各物體假定無初速度且不受除重力以外的外力作用。進一步針對更為復雜的工況條件開展鏈式多體系統動力學分析具有重要的工程實踐意義。