1.2 ADAMS多體系統動力學的建模、分析和計算方法

ADAMS采用世界上廣泛流行的多剛體系統動力學理論中的拉格朗日方程方法建立系統的動力學方程。它選取系統內每個剛體質心在慣性參考系的3個直角坐標和確定剛體方位的3個歐拉角作為笛卡兒廣義坐標，用帶乘子的拉格朗日方程處理具有多余坐標的完整約束系統或非完整約束系統，導出以笛卡兒廣義坐標為變量的運動學方程。

ADAMS的計算程序應用了吉爾（Gear）的剛性積分算法以及稀疏矩陣技術，大大提高了計算效率。

1.2.1 廣義坐標的選擇

動力學方程的求解速度在很大程度上取決于廣義坐標的選擇。研究剛體在慣性空間中的一般運動時，用它的連體基的原點（一般與質心重合）確定位置，用連體基相對慣性基的方向余弦矩陣確定方位。

為了解析和描述方位，必須規定一組轉動廣義坐標表示方向余弦矩陣。

● 第一種方法是用方向余弦矩陣本身的元素作為轉動廣義坐標，但是變量太多，同時還要附加6個約束方程。

● 第二種方法是用歐拉角或卡爾登角作為轉動坐標，它的算法規范，缺點是在逆問題中存在奇點，在奇點位置附近數值計算容易出現困難。

● 第三種方法是用歐拉參數作為轉動廣義坐標，它的變量不太多，由方向余弦計算歐拉角時不存在奇點。

ADAMS軟件用剛體的質心笛卡兒坐標和反映剛體方位歐拉角作為廣義坐標，由于采用了不獨立的廣義坐標，因此系統動力學方程雖然是最大數量，但卻是高度稀疏耦合的微分代數方程，適用于稀疏矩陣的高效求解。

1.2.2 多體系統動力學研究狀況

多體系統動力學的核心問題是建模和求解問題，其系統研究開始于20世紀60年代。從20世紀60年代到80年代，多體系統動力學側重于多剛體系統的研究，主要是研究多剛體系統的自動建模和數值求解。

到了20世紀80年代中期，多剛體系統動力學的研究已經取得一系列成果，尤其是建模理論趨于成熟，不過更穩定、更有效的數值求解方法仍然是研究的熱點。

20世紀80年代之后，多體系統動力學的研究更偏重于多柔體系統動力學。這個領域也正式被稱為計算多體系統動力學，至今仍然是力學研究中最有活力的分支之一，但已經遠遠超過一般力學的涵義。

1．多體系統動力學研究的發展

機械系統動力學分析與仿真是隨著計算機技術的發展而不斷成熟的，多體系統動力學是其理論基礎。計算機技術自誕生以來，幾乎滲透到科學計算和工程應用的每一個領域。

數值分析技術與傳統力學的結合曾在結構力學領域取得了輝煌的成就，出現了以ANSYS、NASTRAN等為代表的應用極為廣泛的結構有限元分析軟件。

計算機技術在機構的靜力學分析、運動學分析、動力學分析以及控制系統分析上的應用在20世紀80年代形成了計算多體系統動力學，并產生了以ADAMS和DADS為代表的動力學分析軟件。兩者共同構成計算機輔助工程（CAE）技術的重要內容。

多體系統是指由多個物體通過運動副連接的復雜機械系統。多體系統動力學的根本目的是應用計算機技術進行復雜機械系統的動力學分析與仿真。它是在經典力學基礎上產生的新學科分支，在經典剛體系統動力學上的基礎上經歷了多剛體系統動力學和計算多體系統動力學兩個發展階段，目前已趨于成熟。

多剛體系統動力學是基于經典力學理論的，多體系統中最簡單的情況（自由質點）和一般簡單的情況（少數多個剛體）是經典力學的研究內容。

多剛體系統動力學就是為多個剛體組成的復雜系統的運動學和動力學分析建立適宜于計算機程序求解的數學模型，并尋求高效、穩定的數值求解方法。經典力學逐步發展成多剛體系統動力學，在發展過程中形成了各具特色的多個流派。

早在1687年，牛頓就建立起牛頓方程，解決了質點的運動學和動力學問題。剛體的概念最早由歐拉于1775年提出，他采用反作用力的概念隔離剛體，以描述鉸鏈等約束，并建立了經典力學中的牛頓-歐拉方程。

1743年，達朗貝爾研究了約束剛體系統，區分了作用力和反作用力，達朗貝爾將約束反力稱為“丟失力”，并形成了虛功原理的初步概念。

1788年，拉格朗日發表了《分析力學》，系統地研究了約束機械系統。他系統地考慮了約束，并提出了廣義坐標的概念，利用變分原理考慮系統的動能和勢能，得出第二類拉格朗日方程——最少數量坐標的二階常微分方程（ODE）；并利用約束方程與牛頓定律得出帶拉格朗日乘子的第一類拉格朗日方程——最大數量坐標的微分代數方程（DAE）。

虛功形式的動力學普遍方程尚不能解決具有非完整約束的機械系統問題，1908年若丹給出了若丹原理——虛功率形式的動力學普遍方程，利用若丹原理可方便地討論碰撞問題和非完整系統的動力學問題。

對于由多個剛體組成的復雜系統，理論上采用經典力學的方法，即以牛頓-歐拉方法為代表的矢量力學方法和以拉格朗日方程為代表的分析力學方法。

這種方法對于單剛體或者少數幾個剛體組成的系統是可行的，但隨著剛體數目的增加，方程復雜度成倍增長，尋求其解析解往往是不可能的。后來計算機數值計算方法的出現使得面向具體問題的程序數值方法成為求解復雜問題的一條可行道路，即針對具體的多剛體問題列出其數學方程，再編制數值計算程序進行求解。

對于每一個具體的問題都要編制相應的程序進行求解，雖然可得到合理的結果，但是這個過程長期的重復是讓人不可忍受的，于是尋求一種適合計算機操作的程序化的建模和求解方法就變得非常迫切了。在這個時候，也就是20世紀60年代初期，航天領域和機械領域分別展開了對于多剛體系統動力學的研究，并且形成了不同派別的研究方法。

最具代表性的幾種方法是羅伯森-維滕堡（Roberson-Wittenburg）方法、凱恩（Kane）方法、旋量方法和變分方法。

（1）羅伯森與維滕堡于1966年提出一種分析多剛體系統的普遍性方法，簡稱為R/W方法。這種方法的主要特點是利用圖論的概念及數學工具描述多剛體系統的結構，以鄰接剛體之間的相對位移作為廣義坐標，導出適合于任意多剛體系統的普遍形式動力學方程，并利用增廣體概念對方程的系數矩陣做出物理解釋。

R/W方法以十分優美的風格處理了樹結構多剛體系統；對于非樹系統，通過鉸切割或剛體分割方法將非樹系統轉變成樹系統進行處理。

（2）凱恩方法是在1965年左右形成的分析復雜系統的一種方法，利用廣義速率代替廣義坐標描述系統的運動，直接利用達朗伯原理建立動力學方程，并將矢量形式的力與達朗伯慣性力直接向特定的基矢量方向投影以消除理想約束力，兼有矢量力學和分析力學的特點，既適用于完整系統，也適用于非完整系統。

（3）旋量方法是一種特殊的矢量力學方法（或牛頓-歐拉方法，簡稱為N/E方法）。其特點是將矢量與矢量矩合為一體，采用旋量的概念，利用對偶數作為數學工具，使N/E方程具有極其簡明的表達形式，在開鏈和閉鏈空間機構的運動學和動力學分析中得到廣泛運用。

（4）變分方法是不同于矢量力學或分析力學的另一類分析方法。高斯最小拘束原理是變分方法的基本原理。保保夫和里洛夫從這一原理出發，發展了兩種不同風格的計算方法。該方法有利于結合控制系統的優化進行綜合分析，而且其不受鉸的約束數目的影響，適用于帶多個閉環的復雜系統。

這幾種方法構成了早期多剛體系統動力學的主要內容，借助計算機數值分析技術，解決由多個物體組成的復雜機械系統動力學分析問題，但是多體系統動力學在建模與求解方面的自動化程度相對于結構有限元分析的成熟來說相差甚遠。

正是為了解決多體系統動力學建模與求解的自動化問題，美國Chace和Haug于20世紀80年代提出了適宜于計算機自動建模與求解的多剛體系統笛卡兒建模方法。這種方法不同于以羅伯森-維滕堡方法為代表的拉格朗日方法，以系統中的每個物體為單元，建立固結在剛體上的坐標系。剛體的位置相對于一個公共參考基進行定義，其位置坐標統一為剛體坐標系基點的笛卡兒坐標與坐標系的方位坐標，再根據鉸約束和動力學原理建立系統的數學模型進行求解。

20世紀80年代，Haug等人確立了“計算多體系統動力學”這門新的學科，使多體系統動力學的研究重點由多剛體系統走向側重多柔體系統、柔性多體系統動力學也成為計算多體系統動力學的重要內容。

柔性多體系統動力學在20世紀70年代逐漸引起人們的注意。一些系統（如高速車輛、機器人、航天器、高速機構、精密機械等）中柔性體的變形對系統的動力學行為產生了很大影響。

二十多年來，柔性多體系統動力學一直是研究熱點，這期間產生了許多新的概念和方法，有浮動標架法、運動-彈性動力學方法、有限段方法以及絕對節點坐標法等，其中浮動標架法最早是在航天領域研究中提出來的。

計算多體系統動力學是指用計算機數值手段來研究復雜機械系統的靜力學分析、運動學分析、動力學分析以及控制系統分析的理論和方法。相比于多剛體系統，對于柔性體和多體與控制混合問題的考慮是其重要特征。其具體任務為：

● 建立復雜機械系統運動學和動力學程序化的數學模型。要開發實現這個數學模型的軟件系統，用戶只需輸入描述系統的最基本數據，借助計算機就能自動進行程序化處理。

● 開發和實現有效的處理數學模型的計算方法與數值積分方法，自動得到運動學規律和動力學響應。

● 實現有效的數據后處理，采用動畫顯示、圖表或其他方式提供數據處理結果。

計算多體系統動力學的產生極大地改變了傳統機構動力學分析的面貌，使工程師從傳統的手工計算中解放了出來，只需根據實際情況建立合適的模型，就可由計算機自動求解，并可提供豐富的結果分析和利用手段；對于原來不可能求解或求解極為困難的大型復雜問題，也可利用計算機的強大計算功能順利求解。而且現在的動力學分析軟件提供了與其他工程輔助設計或分析軟件的強大接口功能，與其他工程輔助設計和分析軟件一起提供了完整的計算機輔助工程（CAE）技術。

2．多體系統動力學研究活動

自20世紀60年代以來，國內外在多體系統動力學方面多次召開了深具意義的會議。國際范圍內的會議有1977年由國際理論和應用力學學會（International Union of Theoretical and Applied Mechanics - IUTAM）發起在德國慕尼黑由Magnus主持召開的第一次多剛體系統動力學討論會。自此，國際多體系統動力學研究活動如雨后春筍般涌現。

國內有影響的教材和專著主要有：

● 劉延柱，洪嘉振，楊海興．多剛體系統動力學[M]．北京：高等教育出版社，1989

● 黃文虎，邵成勛．多柔體系統動力學[M]．北京：科學出版社，1996

● 陸佑方．柔性多體系統動力學[M]．北京：高等教育出版社，1996

● 洪嘉振．計算多體系統動力學[M]．北京：高等教育出版社，1999

技巧提示以上教材僅提供參考，有關多體力學的書籍很多，讀者可自行查閱學習。

3．多體系統動力學研究現狀

計算多體系統動力學中所研究的多體系統根據系統中物體的力學特性可分為多剛體系統、柔性多體系統和剛柔混合多體系統。

多剛體系統是指忽略系統中物體的彈性變形而將其當作剛體來處理的系統，該類系統常處于低速運動狀態。柔性多體系統是指系統在運動過程中會出現物體的大范圍運動與物體的彈性變形的耦合，從而必須把物體當作柔性體處理的系統，大型、輕質而高速運動的機械系統常屬此類。如果柔性多體系統中有部分物體當作剛體來處理，那么該系統就是剛柔混合多體系統，這是多體系統中最一般的模型。

1.2.3 多體系統建模理論

對于多剛體系統，從20世紀60年代到80年代，在航天和機械兩個領域形成了兩類不同的數學建模方法，分別稱為拉格朗日方法和笛卡兒方法；20世紀90年代，在笛卡兒方法的基礎上又形成了完全笛卡兒方法。這幾種建模方法的主要區別在于對剛體位形描述的不同。

航天領域形成的拉格朗日方法是一種相對坐標方法，以Roberson-Wittenburg方法為代表，它以系統每個鉸的一對鄰接剛體為單元，以一個剛體為參考物，另一個剛體相對該剛體的位置由鉸的廣義坐標（又稱拉格朗日坐標）來描述，廣義坐標通常為鄰接剛體之間的相對轉角或位移。

這樣開環系統的位置完全可由所有鉸的拉格朗日坐標陣q所確定。其動力學方程的形式為拉格朗日坐標陣的二階微分方程組，即

這種形式首先在解決拓撲為樹的航天器問題時推出。其優點是方程個數最少，樹系統的坐標數等于系統自由度，而且動力學方程易轉化為常微分方程組（Ordinary Differential Equations, ODES）。但方程呈嚴重非線性，為使方程具有程序化與通用性，在矩陣A與B中常常包含描述系統拓撲的信息，其形式相當復雜，而且在選擇廣義坐標時需人為干預，不利于計算機自動建模。不過目前對于多體系統動力學的研究比較深入，采用拉格朗日方法的幾種應用軟件也取得了較好的效果。

對于非樹系統，拉格朗日方法要采用切割鉸的方法以消除閉環，這引入了額外的約束，使得產生的動力學方程為微分代數方程，不能直接采用常微分方程算法去求解，需要專門的求解技術。

機械領域形成的笛卡兒方法是一種絕對坐標方法，即Chace和Haug提出的方法，以系統中每一個物體為單元，建立固結在剛體上的坐標系，剛體的位置相對于一個公共參考基進行定義，其位置坐標（也可稱為廣義坐標）統一為剛體坐標系基點的笛卡兒坐標與坐標系的方位坐標，方位坐標選用歐拉角或歐拉參數。

單個物體位置坐標在二維系統中為3個、三維系統中為6個（如果采用歐拉參數為7個）。對于由N個剛體組成的系統，位置坐標陣q中的坐標個數為3N（二維）或6N（或7N）（三維）。由于鉸約束的存在，因此這些位置坐標不獨立。系統動力學模型的一般形式可表示為

式中，Φ為位置坐標陣q的約束方程，Φq為約束方程的雅可比矩陣，λ為拉格朗日乘子。這類數學模型就是微分-代數方程組（Differential Algebraic Equations, DAES），也稱為歐拉-拉格朗日方程組（Euler-Lagrange Equations），其方程個數較多，但系數矩陣呈稀疏狀，適宜于計算機自動建立統一的模型進行處理。笛卡兒方法對于多剛體系統的處理不區分開環與閉環（樹系統與非樹系統），統一處理。目前國際上最著名的兩個動力學分析商業軟件ADAMS和DADS都采用這種建模方法。

完全笛卡兒坐標方法由Garcia和Bayo于1994年提出，是另一種形式的絕對坐標方法。這種方法的特點是避免使用一般笛卡兒方法中的歐拉角或歐拉參數，而是利用與剛體固結的若干參考點和參考矢量的笛卡兒坐標描述剛體的空間位置與姿態。

參考點選擇在鉸的中心，參考矢量沿鉸的轉軸或滑移軸，通常可由多個剛體共享而使未知變量減少。完全笛卡兒坐標所形成的動力學方程與一般笛卡兒方法本質相同，只是其雅可比矩陣為坐標線性函數，便于計算。

至于柔性多體系統，從計算多體系統動力學角度看，柔性多體系統動力學的數學模型首先應該和多剛體系統與結構動力學有一定的兼容性。當系統中的柔性體變形不計時就退化為多剛體系統，當部件間的大范圍運動不存在時就退化為結構動力學問題。

柔性多體系統不存在連體基，通常選定一個浮動坐標系描述物體的大范圍運動，物體的彈性變形將相對該坐標系定義。彈性體相對于浮動坐標系的離散將采用有限單元法與現代模態綜合分析方法。

在用集中質量有限單元法或一致質量有限單元法處理彈性體時用結點坐標來描述彈性變形。

在用正則模態或動態子結構等模態分析方法處理彈性體時用模態坐標描述彈性變形。這就是萊肯斯首先提出的描述柔性多體系統的混合坐標方法，即用坐標陣p=(qTaT)T描述系統的位形，其中q為浮動坐標系的位形坐標、a為變形坐標。考慮到多剛體系統的兩種流派，在柔性多體系統動力學中也相應提出兩種混合坐標，即浮動坐標系的拉格朗日坐標加彈性坐標與浮動坐標系的笛卡兒坐標加彈性坐標。

根據動力學基本原理推導的柔性多體系統動力學方程形式同式（1-1）和式（1-2），只是將q用p代替，即柔性多體系統具有與多剛體系統類似的動力學數學模型。

1.2.4 多體系統動力學數值求解

多剛體系統拉格朗日方法產生的形如式（1-1）的動力學數學模型是形式復雜的二階常微分方程組（ODES），系數矩陣包含描述系統拓撲的信息。對于該類問題的求解，通常采用符號-數值相結合的方法或者全數值的方法。

符號-數值方法是先采用基于計算代數的符號計算方法進行符號推導，得到多剛體系統拉格朗日模型系數矩陣簡化的數學模型，再用數值方法求解ODE問題。鑒于計算機技術的發展，目前全數值方法也較為流行，就是將多剛體系統拉格朗日數學模型當作一般ODE問題進行求解，這方面的技術已經較為成熟。

多剛體系統笛卡兒方法產生的形如式（1-2）的動力學數學模型是著名的微分-代數方程組（DAES）。DAE問題是計算多體系統動力學領域的熱點問題。

柔性多體系統的動力學數學模型形式與多剛體系統相同，借鑒多剛體系統數學模型的求解方法。只是混合坐標中描述浮動坐標系運動的剛體坐標 q 通常是慢變大幅值的變量，而描述相對于浮動坐標系彈性變形的坐標 a 卻為快變微幅的變量，兩類變量出現在嚴重非線性與時變的耦合動力學方程中，其數值計算呈病態，將出現多剛體系統中見不到的數值計算困難。

綜上所述，多體系統動力學問題的求解集中于微分-代數方程組的求解。下面將簡要地介紹一下DAE問題的求解方法。

1．微分-代數方程組的特性

多剛體系統采用笛卡兒方法建模生成的微分-代數方程組為：

其中，q、∈Rn分別是系統位置、速度、加速度向量，λ∈Rm是拉格朗日乘子，t∈R是時間，M∈Rn×n為機械系統慣性矩陣，Φq∈Rm×n為約束雅可比矩陣，Q∈Rn為外力向量，Φq∈Rm為位置約束方程。

將式（1-4）對時間求一階和二階導數，得到速度和加速度約束方程：

其中，υ=- Φtq t( , )稱為速度右項，η=-(Φq-2 Φqt-Φtt稱為加速度右項。

給定方程組初始條件：

微分-代數方程組的特性和需要注意的問題有：微分-代數方程問題不是常微分方程（ODE）問題；由式（1-3）和（1-4）組成的微分-代數方程組是指標3問題，通過對約束方程求導，化為由式（1-3）～（1-6）組成的微分-代數方程組后，其指標降為1；微分-代數方程數值求解的關鍵在于避免積分過程中代數方程的違約現象；初值式（1-7）與位置約束式（1-4）及速度約束式（1-5）的相容性；微分-代數方程組的剛性問題。

2．微分-代數方程組積分技術

自20世紀70年代以來，國際上對微分-代數方程問題做了大量的研究，新的算法不斷涌現。根據對位置坐標陣和拉格朗日乘子處理技術的不同，將微分-代數方程組問題的處理方法分為增廣法和縮并法。

（1）增廣法

傳統的增廣法是把廣義坐標加速度和拉格朗日乘子λ作為未知量同時求解，再對加速度進行積分求出廣義坐標速度及廣義坐標位置q，包括直接積分法和約束穩定法。近十年來，在傳統增廣法的基礎上又發展形成了超定微分-代數方程組（ODAEs）方法等新的算法。

● 直接積分法：將式（1-3)和（1-6）聯立在一起，同時求出與λ，然后對積分得和q。該方法未考慮式（1-4）和（1-5）的坐標和速度違約問題，積分過程中誤差積累嚴重，極易發散。在實際的數值計算過程中，并不直接采用直接積分法，但在直接積分法的基礎上發展了一系列控制違約現象的數值方法。

● 約束穩定法：將控制反饋理論引入微分-代數方程組的數值積分過程以控制違約現象。通過把式（1-6）右邊量替換為含位置約束和速度約束的參數式，保證位置約束和速度約束在式（1-3）和（1-6）聯立求解時恒滿足。該方法穩定性好，響應快，但如何選擇參數式中速度項和位置項適當的系數是一個問題。

● 超定微分-代數方程組（ODAEs）法：將系統速度作為變量引入微分-代數方程組，從而將原來的二階DAE化為超定的一階DAE，再為所得方程組引入未知參數，根據模型的相容性消除系統的超定性，如此可使數值計算的穩定性明顯改變；或者將系統位置、速度、加速度向量和拉格朗日乘子向量聯立作為系統廣義坐標，再將由式（1-3）、（1-4）、（1-5）和（1-6）組成的微分-代數方程組及速度與位置、加速度與速度的微分關系式作為約束，化二階DAE為超定的一階DAE，再根據系統相容性引入兩個未知參數，消除超定性，這樣所得的最終約化模型更為簡單，但方程組要多n個。在ODAE方法的基礎上產生了一系列新的更為有效的算法。

● 解耦ODAE法：在ODAE方法的基礎上，發展形成了一類解耦思想，就是在ODAEs基礎上對常用的隱式ODE方法采用預估式，再按加速度、速度和位置的順序進行求解。后來進一步發展形成了無須對隱式ODE方法利用預估式的解耦思想，進一步提高了效率。

（2）縮并法

縮并法就是通過各種矩陣分解方法將描述系統的n個廣義坐標用p個獨立坐標表達，從而將微分-代數方程組從數值上化為與式（1-1）類似的數學模型，以便于用ODE方法進行求解。傳統的縮并法包括LU分解法、QR分解法、SVD分解法以及可微零空間法等，后來在傳統縮并法的基礎上產生了局部參數化縮并方法等新的算法。縮并法中的這些具體方法分別對應著約束雅可比矩陣的不同分解。

● LU分解法：又稱為廣義坐標分塊法。把廣義位置坐標q用相關坐標u和獨立坐標v分塊表示，再將約束雅可比矩陣Φq用LU分解法分塊，得到廣義坐標速度、加速度用獨立坐標速度、加速度表達的式子。將這兩個表達式代入式（1-3），就可得到形如式（1-1）的關于獨立坐標加速度的二階微分方程。該算法可靠、精確，并可控制誤差，但效率稍低。

● QR分解法：通過對約束雅可比矩陣Φq正交分解的結果做微分流型分析，得到可選作受約束系統獨立速度的，并將微分-代數方程組化作二階微分方程，如此可保證在小時間間隔內由積分引起的廣義坐標的變化不會導致大的約束違約。

● SVD分解法：把約束雅可比矩陣Φq做奇異值分解所得結果分別用于式（1-3）和（1-6），得到縮并后的系統動力學方程。在該方法推導過程中沒有用到式（1-4）和（1-5），所以也存在位置和速度違約問題，可用約束穩定法改善其數值性態。

● 可微零空間法：通過Gram-Schmidt正交化過程自動產生約束雅可比矩陣Φq的可微、唯一的零空間基來對系統方程降階。

具體做法是對由Φq∈R m×n和任意矩陣B∈R(n-m ×n)構造的矩陣P∈R n×n采用Gram-Schmidt正交化過程，將P化為正交非奇異矩陣V。

再引入新的速度矢量∈Rn，使滿足=VT，將新速度矢量和加速度矢量按正交化結果分塊，得到新的獨立速度矢量和加速度矢量。如此可將微分-代數方程組化為關于新的獨立加速度矢量的動力學方程。

● 局部參數化縮并方法：先將式（1-3）～（1-6）改寫為等價的一階形式，再用微分流形理論的切空間局部參數化方法將等價的歐拉-拉格朗日方程降為參數空間上的常微分方程。

總的說來，微分-代數方程組數值求解的方法都可歸為增廣法或縮并法。除了上面所介紹的這些增廣法和縮并法所運用的增廣和縮并技術外，近幾年來還出現了不少獨具特色的處理算法，或者是在數值求解算法中獨具匠心，或者是針對某些具體情況做了專門研究。

3．相容性問題和剛性問題

● 初值相容性問題：在微分-代數方程組的數值求解過程中，給定的位置和速度初始條件與微分-代數方程組中的位置和速度約束的相容性是值得注意的一個問題。相容性是微分-代數方程組有解的必要條件。

● 剛性問題：現代機械系統的復雜性會由于系統的耦合而使所得到的微分-代數方程組呈現剛性特性。對于剛性問題的求解，目前最常用的方法是隱式方法。隱式方法不但用于求解剛性問題，而且相比于顯式方法具有更好的穩定性和計算精度。近幾年來，無論是在LU分解法基礎上發展起來的新縮并法，還是基于ODAE方法的增廣法或是基于多體系統正則方程的解法，應用的無不是隱式方法。

1.2.5 計算多剛體系統動力學自動建模

系統的力學模型是對實際問題的力學抽象。要進行動力學的求解，需要由系統的力學模型得到系統的數學模型，這其中的關鍵就在于組裝系統運動方程中所有的系數矩陣。

計算多體系統動力學是基于約束的運動學和動力學，不僅指運動的速度方程和加速度方程是在約束方程的基礎上建立，動力學的運動方程在約束方程的約束下形成微分-代數方程，也指在多體系統動力學分析過程中系統運動方程的各種導數不是實時采用求導算法進行計算，而是采用基于約束的計算方法。

所謂基于約束的計算方法，是指對于有限的約束類型，包括運動學約束和驅動約束，針對每一種約束計算出在系統運動方程中所需的各種導數的相應代數形式，然后在建立數學模型時組裝成系統運動方程中各種導數的組合式。如此，在計算導數時只需代入廣義坐標、時間及其他相關參數即可，避免了導數實時計算所花的大量費用。

1.2.6 多體系統動力學中的剛性問題

剛性（Stiff）問題存在于多剛體系統動力學的某些情形中，更普遍地存在于多柔體系統動力學中，是多體系統動力學的一個重要問題。剛性首先在常微分方程求解理論中提出，并形成了完整的定義和求解理論。常微分方程剛性理論是多體系統動力學中剛性問題的理論基礎，這里先介紹常微分方程剛性問題，再討論多體系統動力學剛性問題。

微分方程的剛性問題是微分方程的一個重要問題，微分-代數方程（DAE）中同樣存在剛性問題。微分-代數方程早期的數值求解中并沒有考慮到這個問題，采用的大多是顯式方法，到了20世紀80年代才發現一些隱式方法不但具有更好的適應性，而且可用于求解剛性問題。

1．剛性方程與剛性穩定性

為描述剛性方程的性質，先考慮線性系統：

其中，y(t)=[y1(t), ..., ym(t)]T∈Rm為解向量函數，?(t)=[?1(t), ...,?m(t)]T∈Rm為已知向量函數，A∈R m×m為常系數矩陣，設其特征值為λj=αj+iβj, j=1,..., m，相應的特征向量為ξj。

定義如果在線性系統中，A的特征值λj(j=1,..., m)應滿足：

，j=1,..., m

則稱式（1-8）為剛性方程，比值s為剛性比，通常剛性比s達到O(10 p)(p≥≥1) 就認為是剛性的。

對于非線性系統：

令(t)為式（1-9）滿足初始條件y(a)=y0的精確解，在解(t)的鄰域內對方程做線性逼近：

或

其中，J(t)是在點(t,(t))處f(t, y)的雅可比（Jacobi）矩陣?f(t, y)/?y的值，如果矩陣J(t)的特征值λj=λj(t)( j=1,..., m)滿足：

則稱方程（1-9）為剛性方程，s(t)稱為在t處的局部剛性比。

在剛性方程中，剛性比s＞＞ 1，矩陣A或J(t)是病態的，故剛性方程也稱為病態方程或壞條件方程。

剛性方程數值積分過程中存在快變分量max|Re(λj)|和慢變分量min|Re(λj)|的差別，快變分量要求積分步長很小，而慢變分量則使得在該步長條件下計算步數很多，舍入誤差較大，這就是剛性方程數值解的困難所在。

考慮到實際問題中可能會出現單個方程情形，或者矩陣A的特征值有實部為零或實部為很小正數的情形，給出與實際問題更為接近的剛性方程的定義。

定義 若線性系統滿足條件：

（1）矩陣A的所有特征值實部小于不大的正數。

（2）A至少有一個特征值的實部是很大的負數。

（3）對應于具有最大負實部的特征值的解分量變化是緩慢的。

對于剛性方程數值穩定性的討論，一般是針對試驗方程進行的。

定義 一個數值方法以定步長h解試驗方程(1-12)，得到線性差分方程的解yn，當n→∞時，若yn→0，則稱該方法對步長h是絕對穩定的。

定義 一個數值方法稱為 A 穩定的，如果它的絕對穩定域包含整個左半平面Re(hλ)＜0。

對于A穩定，存在如下結論：

（1）任何顯式線性多步法（包括顯式RK方法）不可能是A穩定的。

（2）A穩定的隱式線性多步法的階不超過2。

（3）具有最小誤差常數的2階A穩定隱式線性多步法是梯形法。

A 穩定的條件過于苛刻，滿足條件的數值方法太少。為了突破這個限制，放寬穩定性條件，給出A(α)穩定的定義。再從剛性問題特點出發，給出剛性穩定性定義。

定義 一個數值方法稱為A(α)穩定的，如果它的絕對穩定域包含無限的楔形區域Wα={hλ|-α＜π-arg(hλ)＜α}, α∈(0, π/2)。

定義 一個數值方法稱為剛性穩定的，如果它是收斂的，并存在正常數 D、α、θ使在區域R 1={hλ|Re(hλ)≤-D}是絕對穩定的，而在區域R2={hλ|-D<Re(hλ)<α, |Im(hλ)|<θ}上具有高精度且是絕對穩定或相對穩定的。

2．剛性微分方程的數值方法

常微分方程組初值問題：

其中，y=[y1,..., ym]T∈Rm為解向量，f =[f1,..., fm]T∈Rm為已知向量函數，y0∈Rm為已知初始向量。

對于常微分方程組初值問題，常用的方法有3種。

（1）線性多步法（LMM）

當β0=0時為顯式公式，當β0≠0時為隱式公式。當k=1時稱為單步法，當k＞1時稱為多步法。

（2）預估校正法（PECE）

式（1-15）為顯式預估公式，式（1-16）為隱式校正公式。

（3）龍格-庫塔法（RKM）

如果當j≥i時aij=0，式（1-17）和（1-18）就是顯式RK方法。如果當j＞i時aij=0，而對角元素aii不全為0，式（1-17）和（1-18）就稱為半隱式RK方法；若此時aii均相等，則稱為對角隱式RK方法，簡稱為DIRK方法。

在這些求解常微分方程初值問題數值方法的基礎上產生了求解剛性常微分方程的幾類方法，分別是以BDF方法為代表的線性多步法和隱式龍格-庫塔方法。BDF方法是一類特殊的隱式線性多步法。一般的隱式龍格-庫塔方法計算量較大，這里只給出一類特殊的隱式龍格-庫塔方法：對角隱式RK方法。

（1）向后差分公式（BDF）

BDF方法是隱式線性多步法，為k步k階方法。當k=1,2時，BDF方法是A穩定的，當k=3～6時，BDF方法是A(α)穩定和剛性穩定的。實用的BDF方法只能取k=1～6。

（2）對角隱式龍格-庫塔方法（DIRK）

對角隱式龍格-庫塔方法常用的有2級3階和3級4階兩個公式，都是A穩定的。此外，還有A穩定的4級4階公式，但給不出A穩定的更高階DIRK公式。DIRK方法解高頻振蕩的問題比Gear方法（BDF方法）好，但對高精度問題比不上Gear方法好，且計算量比Gear方法大。