計算多體系統(tǒng)動力學讀書報告

1、摘要：多體系統(tǒng)動力學的基本理論，包括多剛體系統(tǒng)動力學建模、多柔體系統(tǒng)動力學建模、多體系統(tǒng)動力學方程求解及多體系統(tǒng)動力學中的剛性（Stiff）問題。本文以建模和求解為基礎，論述了多體系統(tǒng)建模理論，多系統(tǒng)動力學數(shù)值求解以及計算多剛體系統(tǒng)動力學自動建模 等相關建模與求解技術，提到了應用于機械領域的數(shù)學建模方法一一笛卡爾方法，以及著名的微分一代數(shù)方程組（DAEs）。在多體系統(tǒng)動力學非線性方程組求解中的牛頓一拉夫森方法 以及基于超定微分代數(shù)方程組方法的解耦ODAE方法。關鍵詞：多體系統(tǒng)多柔體系統(tǒng)動力學多剛體系統(tǒng)動力學建模求解1概述多體系統(tǒng)動力學的核心問題是建模和求解問題，其系統(tǒng)研究開始于 20世紀60年

2、代。從60年代到80年代，側重于多剛體系統(tǒng)的研究，主要是研究多剛體系統(tǒng)的自動建模和數(shù)值求 解；到了 80年代中期，多剛體系統(tǒng)動力學的研究已經取得一系列成果，尤其是建模理論趨 于成熟，但更穩(wěn)定、更有效的數(shù)值求解方法仍然是研究的熱點；80年代之后，多體系統(tǒng)動力學的研究更偏重于多柔體系統(tǒng)動力學，這個領域也正式被稱為計算多體系統(tǒng)動力學。多體系統(tǒng)是指由多個物體通過運動副連接的復雜機械系統(tǒng)。多體系統(tǒng)動力學的根本目的是應用計算機技術進行復雜機械系統(tǒng)的動力學分析與仿真。它是在經典力學基礎上產生的新學科分 支，在經典剛體系統(tǒng)動力學上的基礎上，經歷了多剛體系統(tǒng)動力學和計算多體系統(tǒng)動力學兩個發(fā)展階段，目前已趨于成熟

3、。多剛體系統(tǒng)動力學是基于經典力學理論的，多體系統(tǒng)中最簡單的情況一一自由質點和 一般簡單的情況一一少數(shù)多個剛體，是經典力學的研究內容。多剛體系統(tǒng)動力學就是為多個剛體組成的復雜系統(tǒng)的運動學和動力學分析建立適宜于計算機程序求解的數(shù)學模型，并尋求高效、穩(wěn)定的數(shù)值求解方法。由經典力學逐步發(fā)展形成了多剛體系統(tǒng)動力學，在發(fā)展過程中形成了各具特色的多個流派。對于由多個剛體組成的復雜系統(tǒng)，理論上可以采用經典力學的方法，即以牛頓-歐拉方法為代表的矢量力學方法和以拉格朗日方程為代表的分析力學方法。這種方法對于單剛體或者少數(shù)幾個剛體組成的系統(tǒng)是可行的，但隨著剛體數(shù)目的增加， 方程復雜度成倍增長， 尋求其解析解往往是不

4、可能的。后來由于計算機數(shù)值計算方法的出現(xiàn)，使得面向具體問題的程序數(shù)值方法成為求解復雜問題的一條可行道路，即針對具體的多剛體問題列出其數(shù)學方程，再編制數(shù)值計算程序進行求解。對于每一個具體的問題都要編制相應的程序進行求解，雖然可以得到合理的結果，但是這個過程長期的重復是讓人不可忍受的，于是尋求一種適合計算機操作的程式化的建模和求解方法變得迫切需要了。在這個時候，也就是20世紀60年代初期,并且形成了不同派別的在航天領域和機械領域， 分別展開了對于多剛體系統(tǒng)動力學的研究,研究方法。最具代表性的幾種方法是羅伯森-維滕堡（Roberson-Wittenburg ）方法、凱恩（Kane）方法、旋量方法和變

5、分方法。計算多體系統(tǒng)動力學中所研究的多體系統(tǒng)，根據(jù)系統(tǒng)中物體的力學特性可分為多剛體系統(tǒng)、柔性多體系統(tǒng)和剛柔混合多體系統(tǒng)。多剛體系統(tǒng)是指可以忽略系統(tǒng)中物體的彈性變形而將其當作剛體來處理的系統(tǒng)，該類系統(tǒng)常處于低速運動狀態(tài)；柔性多體系統(tǒng)是指系統(tǒng)在運動 過程中會出現(xiàn)物體的大范圍運動與物體的彈性變形的耦合，從而必須把物體當作柔性體處理的系統(tǒng)，大型、輕質而高速運動的機械系統(tǒng)常屬此類；如果柔性多體系統(tǒng)中有部分物體可以當作剛體來處理，那么該系統(tǒng)就是剛柔混合多體系統(tǒng)，這是多體系統(tǒng)中最一般的模型。2多體系統(tǒng)動力學建模2.1多體系統(tǒng)建模理論對于多剛體系統(tǒng)，從 20世紀60年代到80年代，在機械領域形成了的數(shù)學建模方

6、法， 稱為笛卡爾方法。機械領域形成的笛卡爾方法是一種絕對坐標方法，即Chace和Haug提出的方法，以系統(tǒng)中每一個物體為單元，建立固結在剛體上的坐標系，剛體的位置相對于一個公共參考基進 行定義，其位置坐標（也可稱為廣義坐標）統(tǒng)一為剛體坐標系基點的笛卡爾坐標與坐標系的 方位坐標，方位坐標可以選用歐拉角或歐拉參數(shù)。單個物體位置坐標在二維系統(tǒng)中為3個，三維系統(tǒng)中為6個（如果采用歐拉參數(shù)為 7個）。對于由N個剛體組成的系統(tǒng)，位置坐標陣 q中的坐標個數(shù)為 3N （二維）或6N （或7N）（三維），由于鉸約束的存在，這些位置坐標 不獨立。系統(tǒng)動力學模型的一般形式可表示為Aq(q,t)(1-1)式中 為位置

7、坐標陣q的約束方程，q為約束方程的雅可比矩陣，為拉格朗日乘子。這類數(shù)學模型就是微分-代數(shù)方程組（DAEs - Differential Algebraic Equations ），也稱為歐拉-拉 格朗日方程組（Euler-Lagrange Equations ），其方程個數(shù)較多，但系數(shù)矩陣呈稀疏狀，適宜于計算機自動建立統(tǒng)一的模型進行處理。笛卡爾方法對于多剛體系統(tǒng)的處理不區(qū)分開環(huán)與閉環(huán)（即樹系統(tǒng)與非樹系統(tǒng)），統(tǒng)一處理。完全笛卡爾坐標方法，由Garcia和Bayo于1994年提出，是另一種形式的絕對坐標方法。這種方法的特點是避免使用一般笛卡爾方法中的歐拉角或歐拉參數(shù)，而是利用與剛體固結的若干參考點

8、和參考矢量的笛卡爾坐標描述剛體的空間位置與姿態(tài)。參考點選擇在鉸的中心，參考矢量沿鉸的轉軸或滑移軸，通?？捎啥鄠€剛體共享而使未知變量減少。完全笛卡爾坐標所形成的動力學方程與一般笛卡爾方法本質相同，只是其雅可比矩陣為坐標線性函數(shù)， 便于計算。至于柔性多體系統(tǒng)，從計算多體系統(tǒng)動力學角度看，柔性多體系統(tǒng)動力學的數(shù)學模型首 先應該和多剛體系統(tǒng)與結構動力學有一定的兼容性。當系統(tǒng)中的柔性體變形可以不計時，即退化為多剛體系統(tǒng)。當部件間的大范圍運動不存在時，即退化為結構動力學問題。柔性多體系統(tǒng)不存在連體基，通常選定一浮動坐標系描述物體的大范圍運動，物體的彈性變形將相對該坐標系定義。彈性體相對于浮動坐標系的離散將

9、采用有限單元法與現(xiàn)代模態(tài) 綜合分析方法。在用集中質量有限單元法或一致質量有限單元法處理彈性體時，用結點坐標來描述彈性變形。在用正則模態(tài)或動態(tài)子結構等模態(tài)分析方法處理彈性體時用模態(tài)坐標描述 彈性變形。這就是萊肯斯首先提出的描述柔性多體系統(tǒng)的混合坐標方法。即用坐標陣 p （qTaT）T描述系統(tǒng)的位形，其中 q為浮動坐標系的位形坐標，a為變形坐標?？紤]到多剛體系統(tǒng)的兩種流派， 在柔性多體系統(tǒng)動力學中也相應提出兩種混合坐標，即浮動坐標系的拉格朗日坐標加彈性坐標與浮動坐標系的笛卡爾坐標加彈性坐標。根據(jù)動力學基本原理推導的柔性多體系統(tǒng)動力學方程，形式同式（1-1），只是將q用p代替。即，柔性多體系統(tǒng)具有與

10、多剛體系統(tǒng)類同的動力學數(shù)學模型。2.2多體系統(tǒng)動力學數(shù)值求解多剛體系統(tǒng)笛卡爾方法產生的形如式（1-1）的動力學數(shù)學模型，是著名的微分 -代數(shù)方程組（DAEs ）。DAE問題是計算多體系統(tǒng)動力學領域的熱點問題。柔性多體系統(tǒng)的動力學數(shù)學模型，其形式與多剛體系統(tǒng)相同， 可以借鑒多剛體系統(tǒng)數(shù)學模型的求解方法。只是混合坐標中描述浮動坐標系運動的剛體坐標q通常是慢變大幅值的變量，而描述相對于浮動坐標系彈性變形的坐標a卻為快變微幅的變量，兩類變量出現(xiàn)在嚴重非線性與時變的耦合動力學方程中，其數(shù)值計算呈病態(tài)，將出現(xiàn)多剛體系統(tǒng)中見不到的數(shù)值計算困難。綜上所述，多體系統(tǒng)動力學問題的求解集中于微分-代數(shù)方程組的求解，

11、下面將簡要地介紹一下DAE問題的求解方法。1)微分-代數(shù)方程組的特性多剛體系統(tǒng)采用笛卡爾方法建模生成的微分-代數(shù)方程組為:M(q,t)q ：(q,t) Q(q,q,t) 0(q,t)0(1-2)(1-3)t R是時間，M Rnn為機械系統(tǒng)慣性矩陣,Rm n為約束雅可比矩陣,Q Rn為其中，q、q、q Rn分別是系統(tǒng)位置、速度、加速度向量，R“是拉格朗日乘子,外力向量，Rm為位置約束方程。將式(1.3)對時間求一階和二階導數(shù)，得到速度和加速度約束方程：(q,q,t) q(q,t)q (q,t)0(1-4)(q,q,q,t) q(q,t)q (q,q,t)0(1-5)其中，t(q，t)稱為速度右項

12、，(qq)qq 2 qtq tt稱為加速度右項。給定方程組初始條件：q(0) qo小、(1-6)q(0) q。2.3多剛體系統(tǒng)動力學建模計算多體系統(tǒng)動力學分析，首先在于提供一個友好方便的界面以利于建立多體系統(tǒng)的力學模型，并在系統(tǒng)內部由多體系統(tǒng)力學模型得到動力學數(shù)學模型；再者需要有一個優(yōu)良的求解器對數(shù)學模型進行求解，求解器要求效率高、穩(wěn)定性好，并具有廣泛的適應性；最后還需計算多體系統(tǒng)動力學分析，首先在于提供一個友好方便的界面以利于建立多體系統(tǒng)的力學模 型，并在系統(tǒng)內部由多體系統(tǒng)力學模型得到動力學數(shù)學模型；再者需要有一個優(yōu)良的求解器對數(shù)學模型進行求解，求解器要求效率高、穩(wěn)定性好，并具有廣泛的適應性

13、；最后還需要對求解結果提供豐富的顯示查詢手段。這其中的關鍵技術就是自動建模技術和求解器設計，所謂自動建模就是由多體系統(tǒng)力學模型自動生成其動力學數(shù)學模型，求解器的設計則必須結合系統(tǒng)的建模，以特定的動力學算法對模型進行求解。2.4計算多體系統(tǒng)動力學建模與求解一般過程一個機械系統(tǒng)，從初始的幾何模型， 到動力學模型的建立，經過對模型的數(shù)值求解，最后得到分析結果，其流程如圖1.1所示。計算多體系統(tǒng)動力學分析的整個流程，主要包括建模和求解兩個階段。建模分為物理建模和數(shù)學建模，物理建模是指由幾何模型建立物理模型，數(shù)學建模是指從物理模型生成數(shù)學模型。幾何模型可以由動力學分析系統(tǒng)幾何造型模塊所建造，或者從通用幾

14、何造型軟件導入。對幾何模型施加運動學約束、驅動約束、力元和外力或外力矩等物理模型要素，形成表達系統(tǒng)力學特性的物理模型。物理建模過程中，有時候需要根據(jù)運動學約束和初始位置條件對幾 何模型進行裝配。由物理模型，采用笛卡爾坐標或拉格朗日坐標建模方法，應用自動建模技術，組裝系統(tǒng)運動方程中的各系數(shù)矩陣，得到系統(tǒng)數(shù)學模型。 對系統(tǒng)數(shù)學模型，根據(jù)情況應用求解器中的運動學、動力學、靜平衡或逆向動力學分析算法，迭代求解，得到所需的分析結果。聯(lián)系設計目標，對求解結果再進行分析，從而反饋到物理建模過程，或者幾何模型的選擇，如此反復，直到得到最優(yōu)的設計結果。在建模和求解過程中，涉及到幾種類型的運算和求解。首先是物理建

15、模過程中的幾何模 型裝配，圖2.1中稱為“初始條件計算”，這是根據(jù)運動學約束和初始位置條件進行的，是 非線性方程的求解問題；再就是數(shù)學建模，是系統(tǒng)運動方程中的各系數(shù)矩陣自動組裝過程， 涉及大型矩陣的填充和組裝問題；最后是數(shù)值求解，包括多種類型的分析計算， 如運動學分析、動力學分析、靜平衡分析、逆向動力學分析等。運動學分析是非線性的位置方程和線性的速度、加速度方程的求解，動力學分析是二階微分方程或二階微分方程和代數(shù)方程混合問 題的求解，靜平衡分析從理論上講是一個線性方程組的求解問題，但實際上往往采用能量的方法，逆向動力學分析是一個線性代數(shù)方程組求解問題，這里面，最復雜的是動力學微分代數(shù)方程的求解

16、問題，它是多體系統(tǒng)動力學的核心問題。在多體系統(tǒng)建模與求解過程，求解器是核心，這其中涉及的所有運算和求解，如初始條件計算、方程自動組裝、各種類型的數(shù)值求解等等都由求解器所支持，它提供了所需的全部算法。圖1.1計算多體系統(tǒng)動力學建模與求解一般過程224多剛體系統(tǒng)動力學對于受約束的多體系統(tǒng)， 其動力學方程是先根據(jù)牛頓定理，給出自由物體的變分運動方程，再運用拉格朗日乘子定理，導出基于約束的多體系統(tǒng)動力學方程。與運動學分析類似， 先考慮二維多體系統(tǒng)，再討論三維多體系統(tǒng)，并對動力學三種類型的分析：正向動力學、逆向動力學和靜平衡分析分別予以討論。1.二維多剛體系統(tǒng)動力學先給出自由剛體的運動方程， 再根據(jù)拉格

17、朗日乘子定理給出約束多體系統(tǒng)帶乘子的運動方程，并討論系統(tǒng)動力學分析的三種情況和約束反力問題。(1) 自由物體的變分運動方程任意一個剛體構件i ,質量為mi，對質心的極轉動慣量為 J,，設作用于剛體的所有外力向質心簡化后得到外力矢量Fi和力矩ni，若定義剛體連體坐標系 xo y的原點o位于剛體質心，則可根據(jù)牛頓定理導出該剛體帶質心坐標的變分運動方程：廠向幾 FiiJi i nJ 0(1-7)其中，ri為固定于剛體質心的連體坐標系原點o的代數(shù)矢量，i為連體坐標系相對于全局坐標系的轉角，仃與i分別為ri與i的變分。為構件i定義的廣義坐標:r T iTqiri , i(1-8)定義廣義力：QiFiT

18、,nJT(1-9)及質量矩陣：Mi diag (mi, mi, Ji)(1-10)則可將式(2.2-61)寫作虛功原理的形式：qT(Mjqi Qi)0(1-11)式(1-11)是連體坐標系原點固定于剛體質心時用廣義力表示的剛體變分運動方程。其中廣義坐標選取、廣義力及質量矩陣計算分別按式(2-69)、(2-70)及(2-71)進行。(2) 束多體系統(tǒng)的運動方程考慮由nb個構件組成的機械系統(tǒng)，對每個構件運用式(2-72)，組合后可得到系統(tǒng)的變分運動方程為：nbqTM iqi Qi0(1-12)i 1若組合所有構件的廣義坐標矢量、質量矩陣及廣義力矢量，構造系統(tǒng)的廣義坐標矢量、質量矩陣及廣義力矢量為：

19、qr TTq1 ,q2 ,-T T , qnb (1-13)Mdiag (M1, M 2 , M nb)(1-14)QQ1T,Q：,Qnb (1-15)系統(tǒng)的變分運動方程則可緊湊地寫為：qTMq Q 0(1-16)對于單個構件，運動方程中的廣義力同時包含作用力和約束力，但在一個系統(tǒng)中，若只考慮理想運動副約束，根據(jù)牛頓第三定律，可知作用在系統(tǒng)所有構件上的約束力總虛功為零， 若將作用于系統(tǒng)的廣義外力表示為：QA Q1AT,Qf,.,QnbTT(1-17)其中:Mq :QA(1-26)QiA 下人丁，nAT , i 1,2,.,nb(1-18)則理想約束情況下的系統(tǒng)變分運動方程為：qTMq QA0(

20、1-19)式中虛位移q與作用在系統(tǒng)上的約束是一致的。系統(tǒng)運動學約束和驅動約束的組合式為：(q,t)0(1-20)注意，在動力學分析中系統(tǒng)約束方程的維數(shù)不需要與系統(tǒng)廣義坐標維數(shù)相等。如果令n 3 nb,則 q Rn，Rm，且 m n。對式(1-20)微分得到其變分形式為：q q 0(1-21)式(2-19)和(2-21)組成受約束的機械系統(tǒng)的變分運動方程，式(2-19)對所有滿足式(2.2-75)的虛位移 q均成立。為導出約束機械系統(tǒng)變分運動方程易于應用的形式，運用拉格朗日乘子定理對式(2-80)和(2-82)進行處理。bTx 0拉格朗日乘子定理：設矢量b Rn，矢量x Rn，矩陣A Rm n為

21、常數(shù)矩陣，如果有:(1-22)Ax 0對于所有滿足式(1-23 )的X條件都成立。(1-23)則存在滿足式(1-24)的拉格朗日乘子矢量Rm。(1-24)Rm n，運用拉格朗日(1-25)bT x TAx 0其中x為任意的。在式(1-19)和(1-21)中，q Rn，MRn n，QA Rn，Rm，對于任意的 q應滿足:乘子定理于式(1-19)和(1-21)，則存在拉格朗日乘子矢量Mq QAt q T q q Mq T QAT q 0由此得到運動方程的拉格朗日乘子形式:Mq :QA(1-26)式(1-26)還必須滿足位置約束方程、速度約束方程及加速度約束方程，如下：(q,t)0(1-27)(q,

22、q,t) q(q,t)q 0,t(q,t)(1-28)(q,q,q,t) q(q,t)q (q,q,t)0,( qq)qq 2 山 tt (1-29)以上三式其維數(shù)同式(1-20)。式(1-26)、(1-27)、(1-28)和(1-29)組成約束機械系統(tǒng)的完整的運動方程。將式(1-26)與(1-29)聯(lián)立表示為矩陣形式：(1-30)式(1-30)即為多體系統(tǒng)動力學中最重要的動力學運動方程，被稱為歐拉-拉格朗日方程式(1-30)還必須滿足式(1-28)和(1-29)。它是一個微分-代數(shù)方程組，不同于單純的常微分方程組 問題，其求解關鍵在于避免積分過程中的違約現(xiàn)象，此外，還要注意DAE問題的剛性問

23、題。顯然，式(1-23)有且僅有唯一解的充要條件是其系數(shù)矩陣非奇異，但這一條件不利于實 際中的判斷，可以給出更為實用的判斷。如果式(1-23)滿足如下條件：a. Rank q(q,t) m，m n ；b. 對任意 a Ker q(q,t)且 a 0， aT M (q,t)a 0。則式(1-23)中的系數(shù)矩陣是非奇異的，且q和 是唯一確定的。這是 多體系統(tǒng)運動方程 解的存在定理?？梢該?jù)此判斷，如果系統(tǒng)質量矩陣是正定的，并且約束獨立，那么運動方程就有唯一解。實際中的系統(tǒng)質量矩陣通常是正定的，只要保證約束是獨立的，運動方程就會有解。(3) 正向動力學分析、逆向動力學分析與靜平衡分析對于一個確定的約束

24、多體系統(tǒng)，其動力學分析不同于運動學分析，并不需要系統(tǒng)約束方程的維數(shù)m等于系統(tǒng)廣義坐標的維數(shù)n , m n。在給定外力的作用下，從初始的位置和速度，求解滿足位置約束式(1-27)及速度約束式(1-28)的運動方程式(1-30)，就可得到系統(tǒng)的加速度和相應的速度、位置響應，以及代表約束反力的拉格朗日乘子，這種已知外力求運動及約束反力的動力學分析，稱為正向動力學分析。如果約束多體系統(tǒng)約束方程的維數(shù)m與系統(tǒng)廣義坐標的維數(shù) n相等,m n，也就是對系統(tǒng)施加與系統(tǒng)自由度相等的驅動約束，那么該系統(tǒng)在運動學上就被完全確定,223節(jié)的約束方程、速度方程和加速度方程可求解系統(tǒng)運動。在此情況下，式(1-27)的雅可

25、比矩陣是非奇異方陣，即:q(q，t)0(1-31)展開式(1-30)的運動方程，為:Mq T QA(1-32)qq(1-33)由式(1-33)可解得q ,再由式(1-32)可求得 ，拉格朗日乘子就唯一地確定了作用在系統(tǒng)上的約束力和力矩(主要存在于運動副中)。這種由確定的運動求系統(tǒng)約束反力的動力學Mq :QA(1-26)Mq :QA(1-26)分析就是逆向動力學分析。如果一個系統(tǒng)在外力作用下保持靜止狀態(tài),也就是說，如果:那么，就說該系統(tǒng)處于平衡狀態(tài)。將式T QA(1-34)(1-34)代入運動方程式(1-26)，得到平衡方程:(1-35)由平衡方程式(1-35)及約束方程式(1-19)可求出狀態(tài)

26、q和拉格朗日乘子。這種求系統(tǒng)的平衡狀態(tài)及在平衡狀態(tài)下的約束反力的動力學分析稱為(靜)平衡分析。(4) 約束反力對于約束機械系統(tǒng)中的構件i，設其與系統(tǒng)中某構件j存在運動學約束或驅動約束，約Mq :QA(1-26)Mq :QA(1-26)束編號為k。除連體坐標系xoy夕卜，再在構件i上以某點P為原點建立一個新的固定于構 件上的坐標系x Py ，稱為運動副坐標系， 設從坐標系x Py到坐標系xoy的變換矩陣為 G，從坐標系xoy到坐標系xoy的變換矩陣為 Ai，則可導出由約束 k產生的反作用力和 力矩分別為：kT A T kT kFiCi Airi(1-36)kPT T kTkTkT (Bir i

27、)(1-37)以上兩式中，k為約束k對應的拉格朗日乘子，反作用力Fi k和力矩T k均為運動副坐標系x Py中的量。2.三維多剛體系統(tǒng)動力學三維系統(tǒng)的廣義坐標比二維系統(tǒng)復雜得多，使得問題規(guī)模更大。這里討論的是與二維多體系統(tǒng)動力學分析相應的內容，包括微分-代數(shù)混合方程組的建立，三種類型的動力學分析等。只是根據(jù)三維系統(tǒng)情況的不同，給出了角速度表示和歐拉參數(shù)表示的兩種不同形式的運動方程。(1)空間自由剛體的變分運動方程對于空間任意剛體構件i ,令其連體坐標系OiXi yi Zj原點Oi固定于剛體質心，此時連體 坐標系也稱為質心坐標系，設剛體質量為mi，其相對于質心坐標系 OiXiZj的慣性張量為Ji

28、，再設作用在剛體上的總外力Fi，外力相對于質心坐標系 Oi XiyiZi原點的力矩為ni，則相對于剛體質心坐標系的剛體牛頓 -歐拉變分運動方程為：(1-38)廠口山 Fi iTJi iiJi inJ 0其中，ri R3為剛體質心的虛位移，i R3為剛體的虛轉動，ri R3為剛體質心位移，i R3為剛體的在坐標系oi xi yi zi中表示的角速度。(2)空間約束機械系統(tǒng)的運動方程一角加速度形式考慮由nb個剛體組成的空間約束機械系統(tǒng)，系統(tǒng)的廣義坐標選為:(1-39)(1-40)(1-41)(1-42)(1-43)(1-44)r ,/.,訃丁r T TT .TP Pl , P2 , Pnb定義：r

29、 T TT ,Tr A , a , rnbM diagI3mhl 3，minb)FFlFFTTTTTL 1 )2 )nbJ diag (J1, J 2,，J nb)(1-45)nIT(1-46)1T2i.)r T TT J(1-47)(1-48)n m m，nnbdiag ( 1, 2,，nb)運用上述符號將系統(tǒng)中每個構件的牛頓-歐拉方程式(1-38)總和為：TTrTMr F TJJ n 0(1-49)同樣地，在理想約束情況下，對于一個系統(tǒng)只要考慮外力或外力矩，如此，由式(2-112)可得到約束機械系統(tǒng)的變分運動方程：rTMr FATJ J nA 0(1-50)(1-57)nIT(1-46)1

30、(1-57)nIT(1-46)1此式必須適用于所有運動學上所容許的虛位移和虛轉動。系統(tǒng)中的運動學約束和驅動約束聯(lián)立表示為：(r, P,t)K(r,P)D(r, P,t)(1-51)(1-57)nIT(1-46)1(1-57)nIT(1-46)1與運動學約束和驅動約束式子不同的是，動力學約束中允許總體上大于零的自由度出現(xiàn)，設其維數(shù)為nh(nh 6nb)。此外，采用歐拉參數(shù)還須滿足歐拉參數(shù)歸一化約束(其維數(shù)為n b)：P PT P1 1, PT P2 1,., Ps P3 1T 0(1-52)對式(1-51)微分，得到用虛位移和虛轉動表示的變分形式：r r0(1-53)對式(1-50)和(1-53

31、)應用拉格朗日乘子定理，則存在一個拉格朗日乘子矢量，滿足:rTMr FA ； TJJ nA T 0(1-54)式中r和是任意的，由此導出空間約束機械系統(tǒng)的牛頓-歐拉運動方程為：n “TlaMrrF(1-55)JT n A J(1-56)對式(1-51)分別求一次及兩次導數(shù)，得到系統(tǒng)速度方程及加速度方程：(1-57)J (1-58)3(1-63)2J (1-58)式中速度右項和加速度右項 。式(1-57)和(1-58)維數(shù)與式(1-51)相同，為nh。由式(1-55)、(1-56)及(1-58)可以得到方陣形式的系統(tǒng)運動方程為:FAA n(1-59)式(1-59)與約束方程式(1-51)及速度方

32、程式(1-57)一起組成描述 nb個剛體系統(tǒng)運動的微分-代數(shù)方程組。式(2-122)的維數(shù)關系滿足:(6nb nh) (6nb nh) (6nb nh) 1(6nb nh) 1。由于角速度是不可積的，因此應把式(1-59)看作是速度變量及代數(shù)變量的一階微分-代數(shù)混合方程。與二維約束機械系統(tǒng)類似，如果滿足:a.矩陣r,行滿秩;b.對任意 a Ker r,且 a 0,存在 aTM,Ja 0。則式(1-52)中系數(shù)矩陣非奇異，系統(tǒng)運動方程有唯一解。在實際迭代計算中，需要給定初始位置條件r(t0)和p(t0)，以及初始速度條件r(t0)和(t。)，初始條件須分別滿足位置約束方程和速度約束方程:(r,P

33、,t)0(1-60)rr(t)(t)(r,p,t)(1-61)(3)空間約束機械系統(tǒng)的運動方程一歐拉參數(shù)形式運用歐拉參數(shù)導數(shù)與角速度之間的關系，可以得出用歐拉參數(shù)導數(shù)表示的空間機械系統(tǒng)運動方程。在歐拉參數(shù)導數(shù)與角速度關系基礎上，可進一步導出歐拉參數(shù)二階導數(shù)與角加速度關 系，對于單個剛體的這些運動量間關系，在這里一并給出：2Gp(1-62)3(1-63)22Gp(1-64)4(1-65)(1-66)2Gp(1-64)42Gp(1-64)4則用歐拉參數(shù)表示式(1-57)所示的系統(tǒng)速度方程為:J pP(1-67)2Gp(1-64)42Gp(1-64)4用歐拉參數(shù)表示式(1-58)所示的系統(tǒng)加速度方程

34、為:(1-68)J pP式(1-67)和(1-68)維數(shù)同樣為nh。對歐拉參數(shù)式(1-52)微分，得到歐拉參數(shù)速度方程:2Gp(1-64)42Gp(1-64)4(1-69)2Gp(1-64)42Gp(1-64)4其維數(shù)關系為nb 4nb 4nb 1 nb 1，且其中:p2diag(p，p；,., p：)(1-70)對式(2-132)微分，得到歐拉參數(shù)加速度方程：pTTT Tpp p2g p, P2 P，.P3 p(1-71)其維數(shù)關系同式(2-132)。約束方程變分式(1-53)用歐拉參數(shù)表示為：r r p p 0(1-72)歐拉參數(shù)約束方程(2-115)的變分式為：P P 0(1-73)對于

35、用角速度表示的約束機械系統(tǒng)的變分運動方程式(2-113)，根據(jù)式(2-150)替換根據(jù)式(2-125)替換，根據(jù)式(2-127)替換，并利用2.2.3節(jié)系列公式，可得到用歐拉參數(shù)表示的系統(tǒng)變分運動方程：(1-74)(1-75)rTMr FA pT4GTJ Gp 8GTJGp 2GTn A 0其中:G diag (G1,G2,.,Gnb)式(1-74)必須對滿足式(1-72)和(1-73)的所有r和p成立。應用拉格朗日乘子定理于式(1-74)及(1-72)和(1-73),則有拉格朗日乘子 和p使下式對 任意r和p成立。rTMr FA T pT4GTJGp 8GTJGp 2GTn A ； pT p

36、0(2-139)從而得到歐拉參數(shù)形式的約束機械系統(tǒng)變分運動方程：Mr T FA(1-76)4GtJ Gpp pT p 2Gt nA 8GTJ Gp(1-77)由式(1-76)和(1-77)，再加上運動學約束和驅動約束的加速度方程式(1-68)及歐拉參數(shù)約束的加速度方程式(1-71)，可得到方陣形式的空間約束機械系統(tǒng)運動方程為:M0Tr0r04GtJ GTppTpprp000pp00p2Gt nA 8GTJGp(1-78)式(1-78)與式(1-51)和式(1-52)的運動學約束和歐拉參數(shù)歸一化約束以及式(1-67)和式-代數(shù)運動方程的混合方程(1-69)的相應速度方程一起，構成歐拉參數(shù)形式表達

37、的系統(tǒng)微分 組。其維數(shù)關系為：(8nb nh) (8nb nh) (8nb nh) 1 (8nb nh) 1式(1-78)的系數(shù)矩陣非奇異的條件同本節(jié)二維多體系統(tǒng)的運動方程解存在條件。實際求解中需要給出初始位置r(to)和p(to)及初始速度r(to)和p(to)，且初始位置條.(1-79)(1-80)件必須滿足位置約束方程，初始速度條件必須滿足速度約束方程，如下:(r, p,to)0p(P,t。)0(t。)pP(to)(r,p,t。)pp(tc) 0與角速度形式的運動方程不同，式(1-78)、(1-51)、(1-52)、(1-57)和(1-59)組成是二階微分-代數(shù)混合方程組。在角速度形式的

38、運動方程式(1-59)中，系數(shù)矩陣左上方是常數(shù)矩陣M和J，而在歐拉參數(shù)形式的運動方程式中，系數(shù)矩陣中出現(xiàn)了取決于歐拉坐標的項gtjg，這會使兩種形式的運動方程的數(shù)值性態(tài)不盡相同。(4)逆向動力學分析、平衡分析和運動副約束反力空間約束機械系統(tǒng)的逆向動力學分析、平衡分析與與平面約束機械系統(tǒng)的相應方法完全相同。逆向動力學分析，對于運動學上確定的系統(tǒng)，式(1-59)或(1-78)的系數(shù)矩陣非奇異，可以 直接求出加速度和拉格朗日乘子，再進一步得到系統(tǒng)運動狀態(tài)和約束反力。平衡分析，根據(jù)平衡狀態(tài)的定義由式(1-55)和(1-56)或由式(1-76)和(1-77)可以導出系統(tǒng)平 衡方程，再由平衡方程及約束方程

39、可以求出系統(tǒng)的平衡位置及平衡時的拉格朗日乘子，再進一步得到約束反力。運動副約束反力的計算也與平面約束機械系統(tǒng)類似。考慮一典型運動副k，運動副定義kk點為P，約束方程為0 ,相應的拉格朗日乘子為，運動副反作用力和反作用力矩在運動副定義坐標系 Pxyz中表示為Fi k和Ti k，則由拉格朗日乘子k計算運動副反作用力和反作用力矩的公式為:FiCiTAiTkTri(1-78)ct( kTP aT k T kSi Ari)(1-79)式中Ci為從運動副定義坐標系 Px y z到連體坐標系o x y z的方向余弦變換矩陣。2.2.5計算多剛體系統(tǒng)動力學自動建模系統(tǒng)的力學模型是對實際問題的力學抽象，要進行動

40、力學的求解，需要由系統(tǒng)的力學模型得到系統(tǒng)的數(shù)學模型，也就是要得到系統(tǒng)運動方程(2.2-80)(2.2-83)或(2.2-115)、(2.2-107)和(2.2-113)的具體形式，這其中的關鍵在于組裝系統(tǒng)運動方程中所有的系數(shù)矩陣。計算多體系統(tǒng)動力學是基于約束的運動學和動力學，不僅指運動的速度方程和加速度方程是在約束方程的基礎上建立，動力學的運動方程在約束方程的約束下形成微分-代數(shù)方程，也指在多體系統(tǒng)動力學分析過程中，系統(tǒng)運動方程的各種導數(shù)不是實時地采用求導算法進行計算，而是采用基于約束的計算方法。所謂基于約束的計算方法，是指對于有限的約束類型， 包括運動學約束和驅動約束，針對每一種約束計算出在

41、系統(tǒng)運動方程中所需的各種導數(shù)的相 應的代數(shù)形式，然后在建立數(shù)學模型時組裝成系統(tǒng)運動方程中各種導數(shù)的組合式，如此，在計算導數(shù)時只需代入廣義坐標、時間及其它相關參數(shù)即可，避免了導數(shù)實時計算所花的大量費用。根據(jù)系統(tǒng)運動方程(2-87)(2-90)的要求，對于所有的各種運動學約束和驅動約束，需要(q,t)、約束方程雅可比矩陣q(q,t)、速度右項給出如下幾個量的代數(shù)形式：約束方程(q,t)及加速度右項 細列出。(q,v,t)，所有基本約束相應參量的表示可以直接給出，這里不再詳主要是各種系數(shù)現(xiàn)在討論從系統(tǒng)力學模型到系統(tǒng)數(shù)學模型建立過程中的自動建模技術, 矩陣的組裝技術。如式(2-87)(2-90)所述，

42、系統(tǒng)動力學運動方程為：F(q,q,q, ,t) M (q,t)q T(q,t) f(q,q,t) 0(2-147)(q,t)0(2-148)(q,q,t) q(q,t)q (q,t)0,t(q,t)(2-149)(q,q,q,t) q(q,t)q (q,q,t)0,( qq)qq 2 qtq tt (2-150)則需要組裝的矩陣有：M(q,t)、 q(q,t)、f(q,v,t)、(q,t)、(q,t)和(q,v,t)。如果考慮到求解過程中所涉及到的一些運算，還需給出所有約束的其它更復雜參量的代數(shù)形式，如：(q,t)T q、(q,t)aq、(q,t)vq、 q(q,v,t)及 v(q,v,t)等。這些形式都可以采用本節(jié)的統(tǒng)一方法一致獲得。以二維系統(tǒng)為例，三維系統(tǒng)同理。設二維系統(tǒng)由nb個構件和nj個約束(包括運動學約束和驅動約束)組成，并對系統(tǒng)中構件施加力(矩)和力元的作用，再設系統(tǒng)廣義坐標維數(shù)為n ,約化后的約束方程維數(shù)為 的運動參量如下：m ,定義構件i的物理參量和約束j (設其聯(lián)接構件r與s)廣義坐標PP曰.質量極轉動慣量外力矢 外力矩 約束方程 速度右項 加速度右