顯示多體動(dòng)力學(xué)方程自動(dòng)編制技術(shù)

李瑾岳，胡權(quán)，張景瑞

北京理工大學(xué) 宇航學(xué)院，北京 100081

顯示多體動(dòng)力學(xué)方程自動(dòng)編制技術(shù)

李瑾岳，胡權(quán)*，張景瑞

北京理工大學(xué) 宇航學(xué)院，北京 100081

為高效得到復(fù)雜柔性多體航天器的動(dòng)力學(xué)方程解析表達(dá)式，研究了系統(tǒng)動(dòng)力學(xué)方程自動(dòng)編制技術(shù)。首先，基于Kane方法推導(dǎo)了構(gòu)型簡(jiǎn)單的鏈狀多體系統(tǒng)動(dòng)力學(xué)顯式方程；然后，根據(jù)方程中各項(xiàng)物理意義和分布規(guī)律，開發(fā)了適用于樹狀構(gòu)型多體航天器的解析動(dòng)力學(xué)方程輸出程序，并運(yùn)用Latex軟件對(duì)輸出文本進(jìn)行二次編譯，增強(qiáng)輸出復(fù)雜數(shù)學(xué)公式的可讀性。所得方法極大地簡(jiǎn)化了對(duì)此類航天器的建模工作。

多柔體系統(tǒng)；機(jī)械臂；顯式動(dòng)力學(xué)方程；自動(dòng)編寫；Kane方程

隨著航天任務(wù)的復(fù)雜化，航天器多為由多個(gè)部件組成的多體系統(tǒng)。建立此類系統(tǒng)的動(dòng)力學(xué)方程成為系統(tǒng)特性分析和控制設(shè)計(jì)的前提。Hooker和Margulies首先將多體動(dòng)力學(xué)的概念運(yùn)用于航天器系統(tǒng)建模1]，基于牛頓/歐拉法推導(dǎo)了任意數(shù)量剛體組成的樹形多體系統(tǒng)動(dòng)力學(xué)方程，但涉及的約束力消元成為該方法的難點(diǎn)。之后，Hooker使用拉格朗日方法避免方程中出現(xiàn)未知約束力2]，但拉格朗日法分析系統(tǒng)計(jì)算量大，得到的方程形式復(fù)雜，不易于分析。同一時(shí)期，達(dá)郎貝爾原理、漢密爾頓方法、波茲曼/海默方法等其他方法也都曾被用于系統(tǒng)分析3]。1960年起，Kane發(fā)表了系列論文4-6]，提出基于新的運(yùn)動(dòng)學(xué)變量——偏速度，簡(jiǎn)化了方程的推導(dǎo)，僅需要將偏速度與已知向量(慣性力、慣性力矩、接觸力、接觸力矩)點(diǎn)乘，即可得到系統(tǒng)的運(yùn)動(dòng)方程。

以上方法可用于數(shù)值推導(dǎo)或符號(hào)推導(dǎo)系統(tǒng)動(dòng)力學(xué)方程。對(duì)系統(tǒng)運(yùn)動(dòng)的解算，依賴方程的數(shù)值計(jì)算與求解；但是，若能獲得系統(tǒng)顯式的動(dòng)力學(xué)方程，則能夠?qū)ο到y(tǒng)特性分析提供模型參考。早期對(duì)于計(jì)算機(jī)符號(hào)分析法應(yīng)用的嘗試未能取得成功，這是由于在對(duì)方程進(jìn)行計(jì)算時(shí)，出現(xiàn)了大量的中間表達(dá)數(shù)，致使計(jì)算機(jī)內(nèi)存不足7]。隨著Kane方法的運(yùn)用，這一問題得到了解決。Levinson首先展示了Kane算法與符號(hào)分析法結(jié)合的可行性8]；幾乎同時(shí)，Schiehlen和Kreuzer宣布了基于符號(hào)分析法及牛頓/歐拉方程研發(fā)的多體動(dòng)力學(xué)代碼，NEWEUL9]。不久之后Rosenthal及Sherman基于凱爾方法開發(fā)了商業(yè)軟件SD/EXACT10]。近些年，另一些符號(hào)動(dòng)力學(xué)程序也被開發(fā)了出來(lái)，如AUTOLEV11]和SYMKIM12]。這些基于符號(hào)分析法的代碼非常適合運(yùn)用于具有少量剛體的系統(tǒng)，但并不能很好地處理具有多個(gè)柔性體的模型。而大多數(shù)可處理具有大量柔性體系統(tǒng)的商業(yè)軟件仍然使用數(shù)值分析法13]。

本研究基于已有的顯式多體系統(tǒng)動(dòng)力學(xué)方程，發(fā)現(xiàn)鏈狀及樹狀多柔體系統(tǒng)的顯式動(dòng)力學(xué)方程形式上的規(guī)律，以及動(dòng)力學(xué)方程形式與系統(tǒng)構(gòu)型間的聯(lián)系，提出一種顯式系統(tǒng)動(dòng)力學(xué)方程高效編寫方法。此方法適用于具有樹狀結(jié)構(gòu)的多柔體系統(tǒng)，并可將柔體退化為剛體。同時(shí)為了提高輸出方程的可讀性，基于Latex排版程序?qū)敵鑫谋具M(jìn)行了二次編輯，便于后續(xù)解析的分析各項(xiàng)的物理含義。最后給出了應(yīng)用實(shí)例。

1 多柔體系統(tǒng)動(dòng)力學(xué)方程的推導(dǎo)

以帶有一部機(jī)械臂系統(tǒng)的空間鏈狀多體系統(tǒng)為例，說(shuō)明系統(tǒng)方程內(nèi)部的固有規(guī)律。如圖1所示，中心體與機(jī)械臂間、機(jī)械臂的3節(jié)臂桿之間采用鉸鏈鉸接。假定每個(gè)鉸鏈具有1～3個(gè)相對(duì)轉(zhuǎn)動(dòng)自由度。

圖1 目標(biāo)系統(tǒng)的構(gòu)型示意Fig.1 Structure of target system

1.1 系統(tǒng)的矢量定義

為方便描述系統(tǒng)及各臂桿運(yùn)動(dòng)，引入中心體體坐標(biāo)系，其原點(diǎn)位于復(fù)雜連接體中心Ob,Xb軸、Yb軸與Zb軸為復(fù)雜連接體中心體的慣性主軸，且構(gòu)成右手坐標(biāo)系；機(jī)械臂臂桿坐標(biāo)系，原點(diǎn)位于機(jī)械臂系統(tǒng)中第j節(jié)機(jī)械臂臂桿與其內(nèi)接體的鉸鏈Oj處，坐標(biāo)系各軸與第j節(jié)臂桿固連；在本系統(tǒng)中，j=1,2,3。

根據(jù)以上定義的位置矢量，可以得到各體質(zhì)量微元在慣性系中的位置矢量14]：

(1)

(2)

式中：Aj,k表示從單位基向量bk到bj的坐標(biāo)轉(zhuǎn)換矩陣。

位置矢量對(duì)時(shí)間求一階導(dǎo)數(shù)，便可求得系統(tǒng)中各體參考質(zhì)量點(diǎn)的速度矢量。

(3) 2)第1節(jié)臂桿質(zhì)量微元：

(4)

3)第2節(jié)臂桿質(zhì)量微元：

(5)

4)第3節(jié)臂桿質(zhì)量微元：

(6)

為簡(jiǎn)化后續(xù)推導(dǎo)，定義：

(7)

(8)

(9)

速度矢量中，uj,j+1為因彈性變形而產(chǎn)生的位移。

速度矢量對(duì)時(shí)間再次求導(dǎo)即可得加速度矢量，這里不再給出加速度矢量的具體形式。

1.2 動(dòng)力學(xué)方程的推導(dǎo)

Kane方程15-16]的一般形式：

(10)

在上述廣義力公式中，Gh為偏速度向量，系統(tǒng)廣義速度表示為v，則有：

(11)

式中：Wh為選定的系統(tǒng)廣義偏速度15]。

由式(11)可得，偏速度向量Gh可表示為速度向量v對(duì)各階偏速度Wh求導(dǎo)所得向量，即：

(12)

將各體質(zhì)量微元的偏速度、加速度代入Kane方程，可得各體的動(dòng)力學(xué)方程。

1)中心體的平動(dòng)方程：

(13)

2)中心體的轉(zhuǎn)動(dòng)方程：

(14)

3)第k節(jié)機(jī)械臂的轉(zhuǎn)動(dòng)方程：

(14)

由以上各項(xiàng)的物理含義，可以發(fā)現(xiàn)顯式系統(tǒng)動(dòng)力學(xué)方程的形式與系統(tǒng)的構(gòu)型有關(guān)，且在形式上規(guī)整，具有很強(qiáng)的規(guī)律性。

2 樹狀構(gòu)型多體系統(tǒng)動(dòng)力學(xué)方程

當(dāng)系統(tǒng)為樹狀構(gòu)型時(shí)，動(dòng)力學(xué)方程形式會(huì)發(fā)生變化，原有顯式方程形式上的規(guī)律也將隨系統(tǒng)構(gòu)型發(fā)生變化。

以圖 2所示系統(tǒng)為例，說(shuō)明方程各項(xiàng)的修正方法。為減小方程推導(dǎo)計(jì)算量，將系統(tǒng)退化為剛性。同上節(jié)推導(dǎo)思路，直接給出系統(tǒng)各體的顯式動(dòng)力學(xué)方程，并找出規(guī)律。

圖2 存在樹狀構(gòu)型的示例系統(tǒng)Fig.2 Target system with tree structure

1)中心體的平動(dòng)方程：

(15)

(16)

3)第1系統(tǒng)中第1節(jié)機(jī)械臂的轉(zhuǎn)動(dòng)方程：

(17)

4)第1系統(tǒng)中第2節(jié)機(jī)械臂的轉(zhuǎn)動(dòng)方程：

(18)

5)第2系統(tǒng)中第1節(jié)機(jī)械臂的轉(zhuǎn)動(dòng)方程：

(19)

通過觀察推導(dǎo)方程形式，得出結(jié)論：

樹狀構(gòu)型系統(tǒng)顯式動(dòng)力學(xué)方程形式與鏈狀構(gòu)型系統(tǒng)顯式動(dòng)力學(xué)方程形式類似，方程中描述中心體運(yùn)動(dòng)的相關(guān)項(xiàng)不變，僅描述機(jī)械臂的運(yùn)動(dòng)的部分項(xiàng)發(fā)生變化，由原來(lái)描述單一機(jī)械臂運(yùn)動(dòng)的單一項(xiàng)變?yōu)槊枋龆鄠€(gè)機(jī)械臂運(yùn)動(dòng)的多個(gè)項(xiàng)的疊加。

3 程序框架介紹

由上節(jié)推導(dǎo)多體航天器系統(tǒng)動(dòng)力學(xué)方程，發(fā)現(xiàn)顯式動(dòng)力學(xué)方程在形式上與系統(tǒng)構(gòu)型關(guān)聯(lián)，形式規(guī)整具有規(guī)律性。通過運(yùn)用顯式方程形式上的規(guī)律，編寫了顯式動(dòng)力學(xué)方程自動(dòng)編制程序。圖3給出了程序流程。

圖3 程序工作流程Fig.3 Flowchart of the program

程序分為兩大部分，一部分存儲(chǔ)動(dòng)力學(xué)方程的所有子項(xiàng)。另一部分負(fù)責(zé)處理輸入構(gòu)型。使用程序時(shí)，根據(jù)提示輸入系統(tǒng)構(gòu)型，包括機(jī)械臂數(shù)量及臂桿數(shù)量，并聲明中心體及每個(gè)臂桿的剛、柔性等其他系統(tǒng)信息。完成信息輸入后，程序會(huì)根據(jù)輸入信息，調(diào)用相應(yīng)的子項(xiàng)，并根據(jù)輸入信息對(duì)子項(xiàng)編號(hào)(即程序自動(dòng)標(biāo)注子項(xiàng)中的上下角標(biāo)信息)。完成子項(xiàng)處理后，根據(jù)輸入系統(tǒng)構(gòu)型信息，結(jié)合上節(jié)得出的方程構(gòu)成規(guī)律，自動(dòng)組合子項(xiàng)，得出目標(biāo)系統(tǒng)的動(dòng)力學(xué)方程。

程序采用C語(yǔ)言編制，因此具有眾多上下角標(biāo)及特殊符號(hào)的數(shù)學(xué)方程難以辨識(shí)。為解決此問題，使程序輸出包含符合Latex程序語(yǔ)言規(guī)則的“.tex”文本文檔，輸出符合數(shù)學(xué)公式閱讀習(xí)慣的文本。這里使用了Windows操作系統(tǒng)上的基于Latex的Ctex編譯器進(jìn)行二次編譯。最終文檔以PDF文檔形式輸出。Latex編輯輸出的文檔相比DOS直接的輸出文檔更易于閱讀。

使用Kane方程推導(dǎo)該系統(tǒng)動(dòng)力學(xué)方程，所得到的方程與程序輸出方程一致，證明了該程序的正確性。

4 軟件應(yīng)用實(shí)例

航天器帶有一部6自由度機(jī)械臂及一部可轉(zhuǎn)動(dòng)天線，航天器中中心體視為柔性體，機(jī)械臂中第1、2、3節(jié)臂桿視為柔性體，4、5、6節(jié)臂桿視為剛體。通過自動(dòng)編制程序?qū)懗銎滹@式動(dòng)力學(xué)方程。將柔性天線視為一部具有1節(jié)臂桿的機(jī)械臂系統(tǒng)。輸入系統(tǒng)構(gòu)型之后，程序開始自動(dòng)編寫，編譯完成后生成equation.pdf文檔。給出部分文檔內(nèi)容,如圖4所示。

5 結(jié)束語(yǔ)

本文研究了樹狀構(gòu)型多體系統(tǒng)動(dòng)力學(xué)顯式方程的自動(dòng)編制方法。首先基于Kane方程推導(dǎo)了鏈狀構(gòu)型多柔體系統(tǒng)的動(dòng)力學(xué)方程，通過觀察、分析方程中各項(xiàng)的物理意義和規(guī)律，編寫了顯式多體動(dòng)力學(xué)方程自動(dòng)編制程序。程序避免了人工推導(dǎo)復(fù)雜多體航天器系統(tǒng)顯式方程，為后續(xù)動(dòng)力學(xué)建模和分析提供了解析的模型基礎(chǔ)。

References)

[1] HOOKER W W,MARGULIES G. The dynamical attitude equation for an n-body satelliteJ]. Journal of the Astronautical Sciences, 1965, 12(4):123-128.2] HOOKER W W. A set of r dynamical attitude equations for an arbitrary n-body satellite having r rotational degrees of freedomJ].AIAA Journal, 1970, 8(7):1205-1207.

[3] KANE T R, LEVINSON D A. Formulation of equations of motion for complex spacecraftJ]. Journal of Guidance and Control, 1980, 3(2):99-112

[4] KANE T R. Partial rates of change of position and of orientationJ]. Bulletin of Mechanical Engineering Education, 1960, 2(17):17-24.

圖4 經(jīng)Latex編譯后的部分文檔示例Fig.4 An example of text compiled by Latex

[5] KANE T R.Dynamics of nonholonomic systemsJ]. Journal of Applied Mechanics, 1961, 28(4):574-578.

[6] KANE T R, WANG C F. On the derication of equations of motionJ]. Journal of the Society for Industrial and Applied Mathematics, 1965, 13(2):487-492.

[7] ARUN K B. Contributions of multibody dynamics to space flight: a brief reviewJ]. Journal of Guidance, Control, and Dynamics, 2003, 26(3):385-394.

[8] LEVINSON D A. Equations of motion for multiple-rigid-body systems via symbolic manipulationJ]. Journal of Spacecraft and Rockets, 1977, 14(8), 479-487.

[9] SCHIEHLEN W, KREUZER E. Symbolic computerized derivation of equations of motionM]∥MAGNUSK. Dynamics of Multibody Systems.Berlin: Springer-Verlag,1978:290-305.

[10] ROSENTHAL D E, SHERMAN M A. High performance multibody simulations via symbolic equation mainpulation and Kane′s methodJ]. Journal of the Astronautical Sciences, 1986, 34(3):223-239.

[11] LEVINSON D A, KANE T R. AUTOLEV—a new approach to multibody dynamicsM]∥SCHIEHLENW. Multibody Systems Handbook. Berlin: Springer-Veralg, 1990:81-102.

[12] KECSKEMETHY A, KRUPP T, HILLER M. Symbolic processing of multiloop mechanism dynamics using closed-form kinematics solutionsJ]. Multibody System Dynamics, 1997, 1(1):23-45.

[13] BODLEY C S, DEVER A D, PARK A C,et al.A digital computer program for thedynamic interaction simulation of controls and structureJ].Oranization Studies,1978(8):149-168.

[14] 黃文虎，邵成勛. 多柔體動(dòng)力學(xué)M]. 北京：科學(xué)出版社，1996:15-18.

[15] PETER C H. Spacecraft attitude dynamicsM]. New York: Dover Publications, 2004.

[16] 胡權(quán)，賈英杰，徐世杰. 多體系統(tǒng)動(dòng)力學(xué)Kane方法的改進(jìn)J]. 力學(xué)學(xué)報(bào)，2011(5):968-972.

HU Q, JIA Y J,XU S J. An improved Kanes method for multibody dynamicsJ]. Chinese Journal of Theoretical and Applied Mechanics, 2011(5):968-972(in Chinese).

(編輯：高珍)

Automatically generating explicit equations of motion for multi-body systems

LI Jinyue,HU Quan*, ZHANG Jingrui

SchoolofAerospaceEngineering,BeijingInstituteofTechnology,Beijing100081,China

The development of the complex multibody spacecraft raises higher requirement to multi-body dynamics modeling and control technology. The explicit dynamics equations of the system are always needed to analyze system′s characteristics and design controllers, but the complexity of system brings difficulty to deriving the equations. In order to solve this problem,the explicit form of dynamics equations′ analytical expression was inferred by using Kane′s method, the chain structure system′s explicit dynamics equations;then, by analyzing this equation, physical meaning and regularities of distribution of each equation item were established;using those regularities, a program that can output the explicit dynamics equations of spacecraft system with tree structure was developed. This program reduces the workload of modeling this type of spacecraft significantly.

multi-body system; robotic arm; explicit formation of dynamics;auto writing; Kane equation

10.16708/j.cnki.1000-758X.2017.0005

2016-07-12；

2016-10-10；錄用日期：2016-11-24；

時(shí)間：2017-02-16 16:48:33

http:∥www.cnki.net/kcms/detail/11.1859.V.20170216.1648.004.html

國(guó)家自然科學(xué)基金(11502018，11572037)

李瑾岳(1993-)，男，碩士研究生，ljy19930530@126.com，研究方向?yàn)榭臻g機(jī)器人姿態(tài)動(dòng)力學(xué)與控制

*通訊作者：胡權(quán)(1987-)，男，助理教授，huquan2690@163.com，研究方向?yàn)楹教炱髯藨B(tài)動(dòng)力學(xué)與控制

李瑾岳,胡權(quán),張景瑞.顯示多體動(dòng)力學(xué)方程自動(dòng)編制技術(shù)J].中國(guó)空間科學(xué)技術(shù), 2017,37(1)：49-56.LIJY,HUQ,ZHANGJR.Automaticallygeneratingexplicitequationsofmotionformulti-bodysystemsJ].ChineseSpaceScienceandTechnology, 2017,37(1)：49-56(inChinese).

V43

A

http:∥zgkj.cast.cn