多體系統(tǒng)動力學(xué)的常用積分器算法*

任輝 周平

（哈爾濱工業(yè)大學(xué)航天學(xué)院，黑龍江 150001）

引言

多體系統(tǒng)動力學(xué)得到的動力學(xué)方程組是微分代數(shù)方程組（Differential-Algebraic Equations，簡稱DAEs），其一般形式特別復(fù)雜.為了簡單起見，本文以最典型的index-3的完整約束DAEs：

為主要對象，介紹一下常見的通用型軟件積分器的常用算法與實現(xiàn)細(xì)節(jié).通常來說，通用型多體動力學(xué)軟件的積分器主要強調(diào)以下幾個方面的特性：

（1）精確性.通過控制積分過程的局部誤差，可以大致保證計算結(jié)果的正確性和精度；

（2）高效性.在滿足給定的誤差條件下，自適應(yīng)地選擇盡可能大的積分步長，保證計算速度；

（3）魯棒性.積分器應(yīng)當(dāng)可以穩(wěn)定求解所有常見類型的多體系統(tǒng)動力學(xué)方程組.

通用型軟件中的動力學(xué)方程組在形式上要比方程組（1）繁瑣得多.除了完整約束方程組之外，很多實際系統(tǒng)中還包含非完整約束方程組、控制器方程組（一階常微分方程組或代數(shù)方程組）、離散時間方程組等等.但從數(shù)值方法的角度來看，最復(fù)雜也最難求解的還是類似方程組（1）的index-3的微分代數(shù)方程組.因此，本文主要以方程組（1）為分析對象，而只在文中簡要討論一下其他情況，主要是帶非完整約束的DAEs的求解方法.

本文只強調(diào)適合大規(guī)模多體系統(tǒng)動力學(xué)仿真軟件的通用型積分器技術(shù)，而不詳細(xì)討論針對某些專門問題的特殊積分器技巧.在通用型軟件中，積分器一般是作為獨立模塊開發(fā)的，主要用來宏觀控制求解進(jìn)程，而盡量少與其它模塊在底層交互.軟件中的方程組模型一般是提前設(shè)定好的.當(dāng)軟件的積分器模塊開發(fā)調(diào)試完成后，一般很少再作大的修改.隨著求解問題領(lǐng)域和規(guī)模的擴大，當(dāng)軟件架構(gòu)需要擴展升級時，積分器模塊往往也需要整體重新設(shè)計和開發(fā)，此時亟需整理和消化前人開發(fā)過程中發(fā)展出來的技術(shù).因此積分器模塊中求解技術(shù)的維護和傳承非常重要.還應(yīng)該指出的是，盡管積分器的理論對程序的開發(fā)有絕對的指導(dǎo)作用，但好的積分器中的絕大部分有效的技術(shù)細(xì)節(jié)往往都來自實踐中的反復(fù)測試、改進(jìn)和提供，而非理論公式推導(dǎo)本身.

本文盡量整理了作者在科研和程序開發(fā)中學(xué)習(xí)和體會到的積分器的一些技術(shù)細(xì)節(jié).文中絕大部分積分器都經(jīng)作者親自開發(fā)并測試過，希望能對我國自主研發(fā)多體系統(tǒng)動力學(xué)軟件有所幫助.

動力學(xué)方程組的特性分類

動力學(xué)過程是非線性時變的.積分器必須能夠自適應(yīng)地追蹤系統(tǒng)的時變特性，并選擇合適的時間步長來復(fù)現(xiàn)真實的物理過程.要做到既不能錯過重要的物理現(xiàn)象，也沒必要采用過密的時間節(jié)點而浪費大量計算資源.盡管待求解的方程組都是方程（1）的形式，但在仿真計算過程中，由于具體系統(tǒng)的動力學(xué)過程千差萬別，積分器的表現(xiàn)也是非常不同的.從動力學(xué)特性上來說，多體系統(tǒng)仿真中常見的DAE系統(tǒng)大致可以分為以下幾類：

（1）剛體為主導(dǎo)的系統(tǒng).這類系統(tǒng)在大型多體系統(tǒng)動力學(xué)中最為常見.由于建模方法不同，所得方程組的數(shù)學(xué)特征也稍有不同.

?通過遞歸（recursive）公式建立的動力學(xué)方程組，其主要特點是：方程的非線性很強，表達(dá)式相當(dāng)復(fù)雜；其Jacobian矩陣經(jīng)常為滿陣，但其規(guī)模相對較小.

?由絕對坐標(biāo)描述，通過拉格朗日方程建立的方程組，其主要特點是：約束方程的個數(shù)很多，可能與動力學(xué)方程的個數(shù)為同一量級；剛體之間的相互作用絕大部分由作用力和約束方程來描述.

一般來說，多剛體系統(tǒng)方程的剛性問題往往并不嚴(yán)重，而且光滑性比較好，往往可以用高階格式并結(jié)合比較大的時間步長進(jìn)行數(shù)值積分.

（2）柔性體動力學(xué)為主導(dǎo)的系統(tǒng).一般來說，在這種多體系統(tǒng)中，剛體的個數(shù)仍然遠(yuǎn)多于柔性體的個數(shù).但每個剛體只有6個自由度，而每個柔性體的模態(tài)坐標(biāo)個數(shù)可能比較多.更重要的是，模態(tài)坐標(biāo)對應(yīng)的固有頻率可能很高，會在DAE系統(tǒng)中引入比較強的剛性，導(dǎo)致Jacobian矩陣的條件數(shù)比較大.實踐中發(fā)現(xiàn)，積分器本身的耗散往往不能有效地衰減模態(tài)坐標(biāo)對應(yīng)的高頻振動，因此推薦在模態(tài)坐標(biāo)中引入少許線性阻尼或結(jié)構(gòu)阻尼以提高仿真計算效率.

（3）大變形體為主導(dǎo)的系統(tǒng).通常包含用幾何精確方法或者絕對節(jié)點坐標(biāo)方法建立的大變形有限元模型，需要引入幾何非線性和/或材料非線性的影響.這類系統(tǒng)的動力學(xué)方程個數(shù)通常遠(yuǎn)遠(yuǎn)多于約束方程個數(shù)，且方程組主要由二階動力學(xué)方程組成，具有特殊的稀疏格式.一般說來，從結(jié)構(gòu)動力學(xué)中發(fā)展出來的方法，例如廣義α方法或者HHT方法，在求解這類問題中被證明是很有效的.

（4）接觸碰撞過程為主導(dǎo)的系統(tǒng).在散體顆粒系統(tǒng)或者傳動系統(tǒng)的動力學(xué)仿真中，經(jīng)常會遇到這一類系統(tǒng).首先，由于方程結(jié)構(gòu)的變拓?fù)湫裕∈杈仃嚽蠼馄鞯姆朙U分解不再能復(fù)用，而根據(jù)具體條件進(jìn)行矩陣LU分解是計算中所必需的.其次，常用的罰函數(shù)方法得到的接觸碰撞動力學(xué)方程組的剛性和高頻振蕩問題都比較嚴(yán)重.最后，這類系統(tǒng)仿真過程對于接觸體的幾何形狀和接觸模型比較敏感，采用不同的算法求出的接觸力可能差異比較大.

（5）有特殊要求的一些多體系統(tǒng).例如在虛擬現(xiàn)實仿真中，要求做到實時仿真，需要用到實時積分；高速旋轉(zhuǎn)系統(tǒng)中，需要對轉(zhuǎn)動進(jìn)行有效精確描述；長時間運行的無耗散系統(tǒng)，對能量守恒要求比較高，需要用到守恒格式的積分器；對最優(yōu)控制設(shè)計得到的系統(tǒng)，可能需要求解的并非初值問題，而是邊值問題等等.這些問題將不在本文中詳細(xì)討論.

實際中遇到的大型多體系統(tǒng)仿真問題，往往是上述幾種類型的組合.因此實踐中需要的通用型積分器，要求能同時求解所有這些類型的方程組.此外，具體的應(yīng)用場合也非常重要，而不同場合考驗的是積分器不同方面的特性.例如，很多設(shè)計過程需要反復(fù)進(jìn)行大量的動力學(xué)仿真以優(yōu)化產(chǎn)品的動力學(xué)性能；而在另一些場合中，實時仿真的要求可能更加重要.一個優(yōu)秀的通用型積分器應(yīng)該在求解常見應(yīng)用問題的時候具有高效性，而同時在不常見的特殊應(yīng)用中也應(yīng)該具有魯棒性.

常見通用型DAE積分器簡介

通用型DAE積分器需要有能力做到：

（1）自適應(yīng)地選擇積分步長；對于變階積分器，還需要自適應(yīng)地選擇階數(shù)；

（2）有效地進(jìn)行誤差估計，既可用于保證計算精確性，又可用于變階變步長；

（3）不連續(xù)性和不穩(wěn)定性現(xiàn)象的檢測與處理，自動進(jìn)行算法調(diào)整和事件解決；

（4）有效處理由于數(shù)值原因?qū)е碌乃俣?、加速度、拉格朗日乘子曲線的不光滑性；

（5）在奇異位形、不光滑點處積分器遇到無法解決的問題時可以自動重啟.

所有的DAE積分器都源自相應(yīng)的常微分方程（Ordinary Differential Equations，簡稱 ODEs）積分器.但DAEs與ODEs在數(shù)值特性上有很大的不同［1］，而求解算法的理論和程序?qū)崿F(xiàn)細(xì)節(jié)也不同.

常用的通用型DAE積分器大部分都是隱式格式，其主流算法通常認(rèn)為主要有三大類：基于向后差分公式（Backward Difference Formula，簡稱BDF）的積分器族，基于廣義 α（Generalized-α）方法的積分器族，與基于隱式龍格庫塔（Implicit Runge-Kutta，簡稱IRK）方法的積分器族.這些方法可以直接求解DAE方程組（1）或其變形.此外，將DAEs轉(zhuǎn)化為ODEs并采用ODE積分器來進(jìn)行計算也是重要的DAE求解思路.

BDF族積分器包含了一系列變階變步長積分器.其具體的公式有好幾種變形格式［2-4］；求解的DAE方程組的index可以為1、2或者3；截斷誤差的階數(shù)可以從1階到5階.對于光滑性比較好的多體系統(tǒng)，例如沒有接觸碰撞的多剛體系統(tǒng)，可以用高階格式結(jié)合大的時間步長進(jìn)行積分，求解效率非常高.這類方法的缺點是其低階格式對低頻運動的阻尼比較大，因此如果系統(tǒng)的光滑性不太好，其動力學(xué)響應(yīng)可能會有比較大的衰減.Gear［5］最早將剛性O(shè)DEs求解中常用的BDF方法引入到index-1的微分代數(shù)方程組的求解中，并將其應(yīng)用于電路系統(tǒng)和化學(xué)過程的動力學(xué)仿真.Orlandea等［6］將其進(jìn)一步推廣到了多體系統(tǒng)動力學(xué)的應(yīng)用場合.經(jīng)過深入研究，Brenan等［7］從數(shù)學(xué)角度證明了BDF方法在一般index-1和index-2 DAE系統(tǒng)中的收斂性；同時指出［8］，對于一般的index-3的DAE系統(tǒng)，該方法的收斂性不能從理論上得到保證.但大型商業(yè)軟件的實踐證明，在一定限制條件下［9］，BDF方法完全可以用來求解index-3的多體動力學(xué)方程組（1），而且到目前為止該方法在所有方法中似乎是相當(dāng)高效的.工業(yè)界中，BDF方法仍然是一些大型商業(yè)軟件的缺省求解器，甚至在很多結(jié)構(gòu)動力學(xué)的應(yīng)用場合，其性能仍然非常卓越.直接用BDF積分器求解方程（1）遇到的最大問題是速度變量的誤差估計比較困難［10，11］.這是因為Jacobian矩陣的條件數(shù)很大，求解帶來的舍入誤差會嚴(yán)重污染速度變量的截斷誤差，往往會導(dǎo)致對速度變量無法進(jìn)行有效的誤差估計.這一問題有幾種解決方案，投影濾波方法［11］即是其中之一.Gear等［12］證明，index-3的多體動力學(xué)系統(tǒng)可以等價地轉(zhuǎn)化為數(shù)值特性更好的index-2的DAE系統(tǒng)，因此可以更穩(wěn)定地利用BDF方法求解.實踐表明，這樣求解的結(jié)果精度更高，尤其是速度和加速度的連續(xù)性和光滑性更好.但同時，得到的index-2的DAE系統(tǒng)的方程規(guī)模會更大，因此計算量會稍微大一些.BDF積分器在實踐中還有很嚴(yán)重的問題；三階以上的BDF格式不是絕對穩(wěn)定的，因此在求解寬頻動力學(xué)問題時很容易有頻率進(jìn)入不穩(wěn)定性區(qū)域，導(dǎo)致階數(shù)和步長的估計不再可靠.在某些應(yīng)用場合甚至需要手動設(shè)置最大階數(shù)為2才能順利計算下去，非常不方便.如何自適應(yīng)地選擇積分器的階數(shù)來避免不穩(wěn)定性問題，這一方面的研究也已取得了一些進(jìn)展［13，14］.作者認(rèn)為，BDF 族積分器的性能還存在著一定的提升空間.

近年來，廣義 α 族方法［15］，包括作為同類格式的HHT方法［16］，受到了更多人的研究和使用.這類方法是傳統(tǒng)Newmark方法的高精度擴展，具有統(tǒng)一的二階截斷誤差精度.與BDF方法不同，這種方法計算中用到了位置、速度與加速度之間的關(guān)系，可以直接離散二階動力學(xué)方程組；計算中修正的量本質(zhì)上是加速度和拉格朗日乘子.這類方法具有低頻不衰減而高頻衰減的優(yōu)良特性［17，18］，魯棒性很好，特別適合求解結(jié)構(gòu)動力學(xué)問題或者柔性多體系統(tǒng)動力學(xué)問題等.此外，這種方法程序設(shè)計比較簡單，且可以推廣到 index-2格式［19］.盡管經(jīng)過推廣，這種方法也被應(yīng)用到了一階動力學(xué)方程組，但其數(shù)值精度有所損失.最后，該族積分器的精度保持為一階或二階，對于光滑的動力學(xué)系統(tǒng)，計算效率偏低.

隱式龍格庫塔算法［20，21］也可以用來求解大型多體系統(tǒng).與廣義α方法類似，這種方法也是單步法，其自動選步長的算法比較簡單；同時，這種方法可以實現(xiàn)變階，對于光滑的動力學(xué)系統(tǒng)，理論上也可以采用高階格式結(jié)合大的時間步長來計算［22］.隱式龍格庫塔算法既可以直接求解index-3的DAE系統(tǒng)［23］例如方程組（1），又可以求解index-2的DAE系統(tǒng)［24］，很適合多體系統(tǒng)動力學(xué)仿真的應(yīng)用［25］.這種方法的缺點是，高階格式的方程規(guī)模比較大，對于大型系統(tǒng)計算效率不高.因此實踐中的隱式龍格庫塔算法大多采用對角格式（Diagonally Implicit Runge-Kutta，簡稱 DIRK）［26，27］.數(shù)值實踐發(fā)現(xiàn)，大部分高階隱式龍格庫塔格式在實際計算中事實上達(dá)不到理論上的階數(shù)［28］，而具有階數(shù)飽和現(xiàn)象［29］.最近這方面的研究有一些新的進(jìn)展［30］，其有效性還有待進(jìn)一步的測試.到目前為止，大型通用型軟件很少選擇隱式龍格庫塔方法作為主要的積分器算法，因此本文中也不再詳細(xì)討論其計算細(xì)節(jié).

以上積分器求解DAEs時都會遇到如下困難：

（1）Jacobian矩陣的相關(guān)問題.例如在約束系統(tǒng)的奇異位形下，Jacobian矩陣可能是奇異的或近似奇異的；在剛性問題中，可能由于Jacobian矩陣條件數(shù)太大［31］，帶來計算精度的降低、誤差估計的失敗或其它數(shù)值問題；

（2）對方程組（1）的直接計算仿真過程中，隱含速度約束往往不能得到有效滿足，而加速度和拉格朗日乘子的曲線往往會有強烈的不連續(xù)性（spikes現(xiàn)象）.此外，約束方程本身的不光滑性也可能會帶來數(shù)值計算中的嚴(yán)重困難；

（3）微分代數(shù)方程解的連續(xù)性往往弱于外激勵的連續(xù)性［1］，高頻激勵常常會導(dǎo)致計算結(jié)果出現(xiàn)一些數(shù)值不連續(xù)現(xiàn)象；而輸入變量或者作用力在時間上的不光滑性會對積分器的魯棒性帶來很大的挑戰(zhàn)；在接觸碰撞等問題中，不連續(xù)性與高頻激勵同時出現(xiàn)，甚至?xí)?dǎo)致計算無法順利開展下去；

（4）大多數(shù)積分器的誤差估計和自適應(yīng)步長選擇策略都是基于動力學(xué)響應(yīng)解的光滑性假設(shè).在各種非光滑情形，步長估計可能無效，往往引起大量的計算失敗，大大增加整體計算量.

除了上面的直接方法以外，將DAEs轉(zhuǎn)化為ODEs并用ODE積分器來進(jìn)行仿真的方法也由來已久.通常認(rèn)為這類方法計算效率不高，魯棒性也較差，因此很少應(yīng)用于商業(yè)軟件的實踐.如果約束方程組足夠光滑，Baumgarte降階方法［32，33］常常計算效果也不錯.其它的一些降階方法要么結(jié)合積分步長選擇罰因子［34］，要么將結(jié)果在約束流形上投影［35］或用罰函數(shù)拉回［36］，要么采用修正的動力學(xué)方程組［37］，更詳細(xì)的內(nèi)容可以在綜述文章［38］中找到.通常這些降階方法會引入一些高頻的自由度，需要用剛性積分器或者特殊的顯式積分器［39，40］來求解.

坐標(biāo)分解方法也是一種相當(dāng)有效的方法，而且計算精度較高.利用約束方程將一部分廣義坐標(biāo)用和自由度個數(shù)相等的其它廣義坐標(biāo)表示出來，從而將動力學(xué)方程組（1）轉(zhuǎn)化為 ODEs［41］，進(jìn)而進(jìn)行求解.一般的多體動力學(xué)系統(tǒng)很難找到在仿真過程中一直適用的廣義坐標(biāo)，必須在某些時刻進(jìn)行坐標(biāo)轉(zhuǎn)換.其廣義坐標(biāo)的選擇可以采用流形方法［42］或者投影方法［43］.Shabana教授團隊提出的兩步隱式稀疏矩 陣 積 分 法［44-46］（two-loop implicit sparse matrix numerical integration，簡稱TLISMNI方法），在每一步積分中都要進(jìn)行坐標(biāo)分解，并保證其同時滿足幾何約束、速度約束和加速度約束.這類積分器計算精度很高，但每步計算量太大.如果能證實其積分步長可以取得相當(dāng)大，足以攤平每步的計算成本，該方法也將適合用于通用型積分器.作者尚未對該類積分器的效果進(jìn)行測試，因此這部分內(nèi)容也將不在本文中多做討論.

1 微分代數(shù)方程組的初始條件

與ODEs不同，DAEs的初始條件不能任意給出 .理論分析表明［8，10，21］，DAEs的初始條件下約束方程需要非常高精度地滿足，計算才能有效地開展下去.方程組（1）中的約束方程組為：

隱含速度約束

和加速度條件C?（q，t）=0，即

實際問題提供的初始條件q0和q?0不一定滿足這些約束方程，其原因可能有以下幾點：

（1）q0和 q?0無法準(zhǔn)確給出，或者由于人為原因?qū)е洛e誤的初值設(shè)定；

（2）實踐中，初值往往是測量出來的.而一般測量誤差遠(yuǎn)大于數(shù)值計算中要求的誤差量級，導(dǎo)致初始約束方程不能滿足計算的精度要求；

（3）設(shè)定的初始數(shù)據(jù)，其中有些量可能相對比較精確，而其它的可能并不精確.因而需要調(diào)整不精確的初始變量以滿足約束方程；

（4）求解過程中可能只確保了幾何約束（2）能精確滿足，而速度約束和加速度約束不一定滿足，導(dǎo)致積分器重啟后的初始條件不合適.

初始條件分析的目的是優(yōu)化初始數(shù)據(jù)，使各種約束條件都能精確滿足，并算出初始加速度和約束力.在此之前，首先要確認(rèn)約束方程組是適定的.

冗余約束分析

約束方程可能是冗余的，甚至是互相矛盾的.在建模過程中，有時也可能故意用冗余約束的方法來保證系統(tǒng)安裝的完整性.冗余約束分析的目的是找出數(shù)值特性最好的一組約束方程組，而刪除所有其它的冗余約束.冗余約束分析的基本方法是對Cq（q0，0）進(jìn)行全選主元的LU分解，并選擇一組最大程度線性無關(guān)的約束方程組，而其余的約束方程即被認(rèn)為是冗余約束而舍棄.

初始奇異位形檢測

多體系統(tǒng)動力學(xué)仿真偶爾會陷入奇異位形中.在一般的計算過程中，由于舍入誤差的存在，多體系統(tǒng)滑入奇異位形的可能性比較小，一般不需要專門的處理.更常見的情形是，多體系統(tǒng)的初始位形本身就是奇異位形，導(dǎo)致在冗余約束分析中，有些非冗余約束可能會被誤判為冗余約束而舍去.此時，仿真計算的模型并不是需要的物理模型，因此有必要進(jìn)行初始奇異位形檢測，從而將這些非冗余約束撿回系統(tǒng).

初始奇異位形檢測的算法是對初始位形進(jìn)行攝動，然后重新進(jìn)行冗余約束分析，驗證得到的獨立約束個數(shù)是否與原來相同.如果不同，說明初始位形有奇異性.發(fā)生了這種情況，需要給出警告，要求用戶輸入更多信息來選擇初始位置.如果用戶不能提供其它信息，就隨機選擇若干個與初始位形非常接近的攝動位形分別出發(fā)進(jìn)行仿真計算，并驗證它們的最終計算結(jié)果是否相近.

初始位形分析

給出初始廣義坐標(biāo)q0，希望求得適合約束方程的初始位形.這可以簡單地通過求解優(yōu)化問題：

而得出.其中Wq為對角權(quán)重矩陣，其權(quán)重可根據(jù)q0的先驗知識而調(diào)整.一般對已知精確的分量，其對角權(quán)重取為106，而對于不確定的分量，其權(quán)重取為1.從方程（6）得到的優(yōu)化方程組為

這可以用牛頓迭代法求解，初始迭代向量取為q0.迭代收斂后，用計算出來的q取代初始變量q0.初始位形分析一般只涉及完整約束和/或狀態(tài)變量本身的約束表達(dá)式，而不需要用到非完整約束等其它形式的約束方程.

初始速度分析

給出初始廣義坐標(biāo)q?0，為了求得適合速度約束方程的初始速度，可以簡單地求解優(yōu)化問題

其中Wv為對角速度權(quán)重矩陣，其權(quán)重系數(shù)也可根據(jù)先驗知識進(jìn)行調(diào)整.得到的優(yōu)化方程組為

它完全是線性方程組，可以直接求解出滿足速度約束方程組的初始速度.然后，用計算出來的q?取代初始變量q?0.應(yīng)該強調(diào)的是，在一般情況下，方程（8）中也應(yīng)該包括動力學(xué)系統(tǒng)所有非完整約束方程，因此通常也需要通過迭代才能求解.

初始加速度和約束力計算

很多積分器需要提供初始加速度和拉格朗日乘子的值，這可以由方程（1）和（4）計算出來：

通過初始值分析，既精簡了約束方程，又得到了更合適的初始坐標(biāo) q0、速度 q?0、加速度 q?0和拉氏乘子λ0，是開展動力學(xué)仿真前的基礎(chǔ)準(zhǔn)備.

2 運動學(xué)問題求解方法

運動學(xué)問題是動力學(xué)系統(tǒng)的自由度為0的特殊情形.此時，約束方程組（2）足以決定整個系統(tǒng)的運動學(xué).另外，此時Jacobian矩陣Cq（q，0）為方陣，且除了個別奇異位形以外是可逆的.運動學(xué)問題實際上就是用Newton迭代求解方程（2）

求解的過程要注意以下幾點：

（2）理論上，迭代時間步長h可以任意選取，但考慮到初始估計（12）必須足夠精確以保證Newton迭代（11）的收斂性，h也不能太大.

（3）由于Jacobian矩陣的生成和LU分解計算都比較昂貴，實用中常用擬Newton方法迭代.

計算得到了tn時刻的qn之后，可以進(jìn)而求出tn時刻的速度，只需代入方程（3），即

求解得出q?n.更進(jìn)一步，要求解tn時刻的加速度，需要用到方程（4），即

最后，很多應(yīng)用中還需要算出約束力，只需求解

即可求出任意時刻的拉格朗日乘子λn，進(jìn)而得到約束力.方程（13）、（14）和（15）中可使用同一Cq矩陣的LU分解.當(dāng)然，運動學(xué)仿真在實用中只占很少份額，絕大部分問題還是需要動力學(xué)仿真方法，這將是下文的主要內(nèi)容.

3 廣義α族積分器

廣義α方法最早是在結(jié)構(gòu)動力學(xué)仿真中提出來的［17］，它與稍早提出的HHT方法［18］實際上是等價的，而它們都是傳統(tǒng)的Newmark方法的推廣.這些方法都是單步法，其一般提法為：

問題 1：已知 tn時刻的狀態(tài)變量 qn、q?n、q?n和 λn，及步長的估計h之后，求出tn+1=tn+h時刻對應(yīng)的變量 qn+1、q?n+1、q?n+1以及 λn+1，并作出誤差估計以判斷其是否滿足誤差條件要求.如果滿足，則估計下一步的時間步長；如果不滿足，則估計更合適的時間步長，并重新開始計算.

傳統(tǒng)Newmark方法可用來求解常微分方程

假設(shè)已知 tn時刻的 qn、q?n和 q?n，則在估計出步長 h 之后，只需將q?n+1當(dāng)作未知量，取近似

其中參數(shù)β和γ用來調(diào)節(jié)計算精度和穩(wěn)定性，且

代入方程（16），得到關(guān)于x的非線性方程組

用Newton-Raphson方法迭代求解出修正量x，進(jìn)而代入（19）中得到 qn+1、q?n+1和q?n+1.

而tn+1時刻的狀態(tài)變量假設(shè)為

其中公式（17）中的加速度替代為稍早時刻的值

而動力學(xué)方程組在較早的時間節(jié)點tn+1-α上離散

其中未知量為 an+1=q?n+1-α，而

HHT積分器直接求解的并不是方程組（1），而是其變形格式（24），而且求解時間節(jié)點也不完全在tn+1時刻，利用an+1近似地插值得到tn+1時刻的加速度，就可以直接求解tn+1時刻的動力學(xué)方程組（1），得到的即為廣義α積分器.這兩種積分器的公式稍有不同，廣義α積分器的參數(shù)穩(wěn)定區(qū)域大于HHT積分器.除此之外，二者的各項性能都極其接近，誤差估計與變步長策略也完全相同.本文將以廣義 α積分器作為主要討論對象.

在廣義α積分器中，每步需求解的未知變量仍然是公式（23）中的an+1，但此時假設(shè)

其中αm和αf兩個參數(shù)一般并不相互獨立，可以用一個自由參數(shù)ρ表示出來，取值 ρ∈[0，1].令

從方程（22）和（26）得出tn+1時刻的狀態(tài)變量

此外，假設(shè)乘子λ具有一定的連續(xù)性，即取

將表達(dá)式（29）和（31）代入方程（1），得到關(guān)于更新量x和z的非線性方程組

用Newton-Raphson方法求解，其Jacobian矩陣為

更新量x可以用來進(jìn)行誤差估計和步長選擇.借鑒Negrut等［16］的推導(dǎo)，可以得出該方法的截斷誤差主項為ch2x；其中c為一個量級為1的常數(shù).此外，由于廣義α方法是二階精度的，其廣義坐標(biāo)的截斷誤差正比于h3；而步長的選擇應(yīng)該滿足每個分量的相對誤差都小于但接近于設(shè)定的誤差限，即

本文中的范數(shù)‖?‖定義如下

如果Ξ>1，則說明誤差條件沒有滿足，需要減小步長；而如果Ξ太小，則說明計算步長太小，會降低計算效率，最好放大步長.一個合適的步長選擇策略是

其中系數(shù)0.9是安全系數(shù).由于這里的計算是用本步得出的后驗誤差來估計下一步的積分步長，其正確性依賴于系統(tǒng)參數(shù)的緩變特性，因此安全系數(shù)是必要的.如果預(yù)測失敗，即Ξ>1，則可以利用新得出的Ξ和方程（36）設(shè)置本步的計算步長重新進(jìn)行計算.由于該估計誤差為后驗誤差，得到的hnew一般能很好地滿足誤差條件.如果仍不滿足誤差要求，直接將時間步長減半重新計算；如果再不滿足，說明在tn和tn+h之間出現(xiàn)了某種間斷現(xiàn)象，需要進(jìn)行間斷檢測，發(fā)現(xiàn)間斷位置tn+h*，運行到該時刻然后重啟積分器.廣義α方法保持二階精度，而且穩(wěn)定性很好.以上的自適應(yīng)步長選擇策略失敗率低，魯棒性好，對很多間斷情形也不需要積分器重啟.由于間斷性檢測程序較復(fù)雜，在附錄的算法3中將暫時忽略.

上述方法計算得出的加速度曲線和拉格朗日乘子曲線上經(jīng)常會出現(xiàn)很多尖刺（Spikes），這種現(xiàn)象是非物理的隨機誤差，與所求解的方程組是index-3的DAEs有關(guān).方程組（1）等價于

該方程需要在整個約束流形上滿足，但數(shù)值離散格式不是在流形上開展的，而是在擴展的相空間（q，v）上進(jìn)行的，這將導(dǎo)致如下問題：

1.即使方程q?-v=0在約束流形的切叢上是點點滿足的，該方程的離散格式往往已經(jīng)稍稍脫離了約束流形，而進(jìn)入了相空間；

2.盡管約束方程的引入強制保證了計算出的q?近似位于約束流形上，但隱含的速度約束方程（3）卻沒有保證會被滿足.

將q?和v作為獨立變量；根據(jù)流形上的常微分方程理論，可以將方程q?-v=0換為約束流形上的等價格式；并將隱含的速度約束顯式地引入方程組，得到新的動力學(xué)方程組

其未知量為q、v、λ 和 μ.該方程組是index-2的微分代數(shù)方程組，其與方程組（1）的等價性早已經(jīng)被Gear等人［12］所證明.下面用廣義α方法來離散方程組（37）.此時，方程（26）和（28）中的速度項仍然被離散為

此時速度項vn+1與q?n+1不再相同，但二者差別不應(yīng)該很大，不妨假設(shè)

計算中，取an+1與an+1其預(yù)測值相同，且假設(shè)

進(jìn)而類似方程（29）可以得出，

代入方程（37）得到

其未知量為x、y、λ和μ，迭代Jacobian矩陣為

由方程（38）可以看出，x是與μ同量級的量.用截斷誤差分析可以證明［47］，該積分格式是一階的.該方法關(guān)于q的截斷誤差近似為hx，而關(guān)于v的截斷誤差近似為hy.同樣記

則Ξ≤1為誤差判定條件；而新步長估計方法仍然基于公式（36）.Jay和Negrut的分析表明［47］，如果采用等時間步長積分，求解index-2方程組（37）的格式（38）是二階精度的；但如果采用變時間步長積分，則離散格式（38）只有一階精度.另外，他們還提出一種修正方法，可以將對應(yīng)的HHT格式變?yōu)槎A精度的，而且該方法似乎也很容易推廣到廣義α方法.但作者還沒有親自測試證實，尚不知該方法的實踐效果如何.

4 BDF族積分器

BDF族積分器有好幾種公式實現(xiàn)［8］.它一般用來求解常微分方程組和index-1、index-2的微分代數(shù)方程組［48］，且結(jié)合一些特別的技巧后也能用來計算index-3的微分代數(shù)方程組（1）.與廣義 α方法不同，BDF方法是多步法，而且是一類變階數(shù)、變步長的方法.多步法的一般提法為：

問題2：已知之前 k+1個相繼時刻tn?i（i=0，1，…，k）對應(yīng)的廣義坐標(biāo) qn-i、速度 q?n-i、加速度q?n-i以及拉格朗日乘子 λn-i，在給出時間步長 h之后，求出 tn+1=tn+h 時刻對應(yīng)的變量 qn+1、q?n+1、q?n+1以及λn+1，并給出計算誤差的估計，判斷計算結(jié)果是否滿足誤差條件.如果滿足，則估計下一步的階數(shù)和時間步長；如果不滿足，則更新本步階數(shù)和時間步長，并重新進(jìn)行計算.

下面以常微分方程為例，介紹兩套常用的BDF公式.常微分方程組的形式為：

準(zhǔn)步長公式

假設(shè)給定一系列等步長為h的時間節(jié)點tn-k，tn-k+1，…，tn，tn+1，以及其對應(yīng)的數(shù)據(jù)xn-k，xn-k+1，…，xn，xn+1.根據(jù)在以 h為步長的均勻時間網(wǎng)格上，后差分算子?如下遞歸定義：

其中j=0，1，2，….用中值定理可以證明，存在某個ξ滿足tn+1-kh≤ξ≤tn+1，使得

給定k≥1，記

遞歸地使用（41），可得?kxn+1=?kxn+δ，?k-1xn+1=?k-1xn+?kxn+δ，…，其通項為

特別當(dāng)j=0時有

另外，可以直接驗證［49］多項式

為唯一滿足 ph（tn-j）=xn-j的k階多項式，其中j=0，1，2，…，k；對此方程求微分可得

將近似表達(dá)式（48）中的n替換為n+1，得到

其中的常數(shù)為

為了與后文符號的統(tǒng)一，記

以及

則方程（46）和（50）可以改為

其中矩陣Rk和Uk為k× k的矩陣，其（i，j）分量分別為

參考公式（47），分別將t=tn-iζh（其中i=0，1，…，k）代入方程（56）得到k個線性方程：

計算流程

（1）利用公式（53）得出 xn+1和 hx?n+1的初始估計，然后將（45）和（49）代入微分方程組得到離散格式，用Newton迭代求出更新量δ；

a）?k+1xn+1=δ已經(jīng)求得；

b）?k+2xn+1=δ-?k+1xn；

c）依次計算 ?jxn+1=?jxn+?j+1xn+1，其中j=k，k-1，…，1；

變步長公式

在這套公式下，均勻網(wǎng)格是不需要的，但理論推導(dǎo)要用到Newton插值公式：[xn]=xn

其中j=1，2，3，….可以用數(shù)學(xué)歸納法和中值定理證明，存在某個ξ，滿足tn-k≤ξ≤tn

為了符號方便，對j=1，2，…，k+1，記

并且記（注意它與公式（41）中的定義沒有關(guān)系）

利用這個公式遞歸可得

以及

是唯一滿足p（tn+1-m）=xn+1-m的k階多項式，其中.直接微分上式可得

將公式（60）代入方程（62）中，簡化得到

其中j=1，2，…，k+1.進(jìn)一步，記

于是公式（61）和（63）可以寫成

計算流程

變步長積分器的步進(jìn)計算為：

（2）利用公式（69）得出xn+1和hx?n+1的初始估計，然后將（68）代入微分方程組得到其離散格式，再用Newton迭代求出更新量δ；

（5）令n←n+1，進(jìn)入下一步循環(huán)的起始步1.

方程組求解細(xì)節(jié)

將表達(dá)式（45）和式（49）或者式（68）代入方程（40）得到關(guān)于δ的非線性方程組

求解方程組（1）用到q和v?q?的差分矩陣

從而這兩個差分矩陣可以算出變量的估計表達(dá)式

因此，廣義坐標(biāo)的BDF預(yù)測-校正公式為

而對應(yīng)速度變量的BDF預(yù)測-校正公式為

可以用Newton迭代法求解，其Jacobian矩陣為

求解得到的δ可以直接用來估計q的誤差，但由方程（76）得出的ε用來估計q?的誤差卻可能會被嚴(yán)重放大［10］，需特殊處理后才能應(yīng)用.這里介紹一種更新誤差估計中的ε的方法，即投影濾波法.深入分析［11］發(fā)現(xiàn)，方程（76）計算出來的速度變量誤差并不是沿所有方向都會被放大；而誤差真正被放大的方向與約束流形的切空間正交.如果將ε投影到約束流形切空間內(nèi)，即

為了改進(jìn)速度變量的誤差估計，并減小加速度和拉氏乘子的計算誤差，可以將方程組（1）轉(zhuǎn)化為與其等價的index-2的方程組［12］

然后用BDF積分器進(jìn)行求解.該方程組的規(guī)模比直接求解方程組（1）大了一倍，所需的計算量也大了很多.具體計算步驟為：將變量的離散格式

以及時間導(dǎo)數(shù)項的離散格式

代入方程組（81），得到離散代數(shù)方程組為

在求解以上Index-2微分代數(shù)方程組并得到ε和δ后，就可以同時判斷 ε和 δ相對于q和v的誤差是否滿足計算要求.計算得出的qn+1和vn+1都是準(zhǔn)確的，但拉格朗日乘子λn+1和加速度

可能還是誤差較大.為了進(jìn)一步提高計算精度，可以將 q、v和 a=v?都變成 index-1 的變量，即方程（81）修改為

該方程組仍然是index-2的DAEs，但其對應(yīng)的index-2的變量為較次要的λ、μ和τ，而所有重要的變量q、v和a都是index-1的變量，在計算中都可以精確地計算出來并進(jìn)行誤差控制.

變步長變階

BDF方法是變階變步長的，可以根據(jù)本步計算的后驗誤差進(jìn)行下一步的階數(shù)和步長的估計.早期的 BDF 族 ODE 積分器，例如 DIFSUB［50］或者LSODI［51］中，通常建議k階格式等步長運行k+2步之后再變階變步長.但在通用型DAE積分器中，為了及時捕捉系統(tǒng)的時變特性，建議每一步都要變階和變步長.為了方便討論，記

作為下一步計算步長仍為h時對應(yīng)的j階格式的誤差預(yù)測，利用某種規(guī)則選擇下一步最合適的階數(shù)和時間步長.計算所需的相關(guān)數(shù)據(jù)都含在矩陣

中，其中備選的階數(shù)為 j=1，2，···，k+1.因此，每步計算中可以隨時降多階，但至多只能升1階.這樣一來，積分器既可以在處理意外（不光滑）事件過程中做到自適應(yīng)降階；又可以在求解光滑的動力學(xué)過程中根據(jù)需要將階數(shù)逐漸增加，提高仿真效率；達(dá)到了高效性與魯棒性之間的有效平衡.

階數(shù)和步長的選擇有很多方法，最簡單的思路是選擇階數(shù)使下一步的預(yù)測步長最大，即選擇

和對應(yīng)步長hnew=0.8?h/εknew，其中0.8為安全系數(shù).這種方法實現(xiàn)容易，其各種變形或修正格式在通用型軟件的積分器中仍被廣泛采用.

一般來說，BDF低階格式計算效率低，而高階格式（特別是k≥3的格式）的穩(wěn)定性較差.此外，如果保持階數(shù)和步長，則上一步的Jacobian矩陣往往可以復(fù)用，會節(jié)省不少計算量.因此在新的階數(shù)和步長的選擇中，一般采用如下原則：盡量不變階也不變步長；升階時候盡量保守.例如：

（1）在Gear的程序DIFSUB［49］中，每步要求升/降階次都不超過1.如果本步的階數(shù)為k，在進(jìn)行下一步階數(shù)選擇時只考慮k-1、k、k+1這三種情況，從中選擇可以在下一步取最大步長的階數(shù).在該程序中，上述選階原則是通過調(diào)整安全系數(shù)來實現(xiàn)的：

還要指出的是，如果變步長［53］或者變階［54］算法不合適，也可能帶來穩(wěn)定性問題，因此通常步長和階數(shù)不能增加得太快.作者經(jīng)過測試發(fā)現(xiàn)，Gear等［53］給出的每一階的最大步長增加量偏于保守，而Skelboe［13］給出的步長增加量在實踐中更為合理，其步長增加上界fk如表1所示.

表1 BDF格式變步長的最大穩(wěn)定增加量Table 1 The maximum stable increment of the variable step-size BDF format

高階格式穩(wěn)定性檢測與自適應(yīng)處理

不同階數(shù)的BDF格式具有不同的穩(wěn)定區(qū)域，其中只有前兩階是絕對穩(wěn)定的，而更高階的算法都在虛軸上有不穩(wěn)定區(qū)域.在寬頻動力學(xué)系統(tǒng)，特別是在含有柔性體或者大變形體的多體系統(tǒng)仿真中，某些頻率可能會滑到這些不穩(wěn)定區(qū)域內(nèi)，造成誤差的嚴(yán)重放大.此時上文的誤差估計方法就不再正確，且階數(shù)和步長選擇策略也不再具有魯棒性，在計算中表現(xiàn)為一種積分階數(shù)的頻繁跳躍現(xiàn)象，并伴隨出現(xiàn)大量的積分步失敗.為了避免這一問題，很多程序（例如LSODI［51］）中建議手動設(shè)置最高階數(shù)為2.為了進(jìn)一步解決這一問題，實踐中還發(fā)展出了幾種其他解決方案：

表2 BDF高階格式的失穩(wěn)判斷參數(shù)Table 2 Instability parameters of high-order BDF format

（3）Stewart［14］進(jìn)一步考慮了步長變化的影響，提出當(dāng)?shù)趎步結(jié)束時，如果下列條件滿足：

則下一步暫時不考慮升階，其中

作者經(jīng)過測試發(fā)現(xiàn)，將Skelboe和Stewart的方案結(jié)合起來，雖然最有效地克服了失穩(wěn)問題，但在實踐中似乎偏于保守，會影響計算效率.DASSL采用的策略雖然也可以在一定程度上避免失穩(wěn)現(xiàn)象，但卻不能保證不會出現(xiàn)失穩(wěn).其工作原理大致如下：如果x（t）至少k+1階光滑可導(dǎo)，一般有

這是因為若系統(tǒng)可以用Taylor展開表達(dá)式逐階近似，高階Taylor項應(yīng)當(dāng)比低階項更小.系統(tǒng)的數(shù)值光滑度可定義為選擇最大的k滿足如下條件

如果 x（t）達(dá)不到k階光滑度，就不滿足BDF算法的基本假設(shè)，只能采用

5 隱式積分器的其他計算細(xì)節(jié)

隱式積分器每步的計算量比較大，但其時間步長較大，總體來說計算效率較高；其穩(wěn)定特性很好，可以很容易地處理方程組的剛性問題.但是，時間步長太大，可能會在計算中不能有效地處理局部的一些突變問題，而需要一些特殊的技巧.

5.1 擬Newton迭代與Jacobian復(fù)用

在求解常微分方程（20）、（70）或者微分代數(shù)方程（32）、（38）、（78）和（84）的過程中，甚至在后文中顯式積分器的約束流形投影過程中，都需要用到Newton迭代或者擬Newton迭代.動力學(xué)仿真中的非線性方程組的求解具有以下特點：

（1）多.每個時間步都要用擬Newton迭代求解；

（2）快.狀態(tài)變量的預(yù)測值通常非常準(zhǔn)確，往往只需2-3次擬Newton迭代即可收斂；

（3）大.待求解方程組的規(guī)模比較大；

（4）規(guī).Jacobian矩陣往往具有特定的稀疏性.

Jacobian矩陣的相關(guān)計算往往是整個仿真計算的瓶頸問題，這是由以下兩個原因?qū)е碌模?/p>

（1）Jacobian矩陣生成過程的計算量很大，無論是用解析表達(dá)式還是數(shù)值差分方法；

（2）Jacobian矩陣的LU分解計算復(fù)雜度比其它計算過程要高，即使采用稀疏矩陣計算也是如此.

一般情況下多體動力學(xué)計算仿真中都用的是擬Newton迭代，以避免大量的Jacobian矩陣計算，并在所有地方都盡量復(fù)用Jacobian矩陣的LU分解.實踐證明，這樣做往往可以將計算速度提高5到10倍.還應(yīng)該指出的是，幾乎在所有的仿真過程中，總會遇到少數(shù)擬Newton迭代收斂很慢甚至失敗的情形，給整個仿真帶來問題.如果這些問題不解決，會導(dǎo)致整個計算的中止.而擬Newton迭代的失敗，大多是由于Jacobian矩陣不再精確造成的.因此，為了解決這一問題，可以設(shè)一個最大迭代步數(shù)（例如4到5次），如果達(dá)到最大迭代步數(shù)還不收斂，就可以認(rèn)為擬Newton迭代失敗，需要重新生成Jacobian矩陣并繼續(xù)進(jìn)行迭代.

擬Newton迭代的一般過程為，為了求解非線性方程組F（x）=0，采用迭代格式

其中x*為真實解.該公式可以用來判斷迭代是否收斂，因為右邊項可以計算，其中σ可以由下式估計得出

為了保證擬Newton迭代快速收斂，需要監(jiān)測收斂速度σ.計算中如果發(fā)現(xiàn)σ>0.9，最好重新計算Jacobian并估計σ.另外，實際計算中，條件（87）可能是不滿足的，特別是在以下情形中：

（2）強剛性問題中，Jacobian矩陣的條件數(shù)很大，近似Jacobian與真實Jacobian矩陣之間很小的差異，也會導(dǎo)致條件（87）不滿足［31］；

（3）在很多問題中，真實解析Jacobian矩陣的生成太過煩瑣，或者其稀疏性不太好，往往用近似Jacobian計算，可能不滿足條件（87）；

（4）在某些場合，非線性問題本身的收斂域可能很小，而這種情況是完全無法提前預(yù)知.

檢測到擬Newton迭代不收斂時需要立即重新計算Jacobian矩陣，并進(jìn)行LU分解.一般說來，擬Newton迭代的收斂速度低于Newton迭代，但計算量遠(yuǎn)小于后者.實際計算中需要結(jié)合二者的優(yōu)點，盡可能地加快計算速度.

總體說來，廣義 α方法計算步長較小，擬Newton方法的收斂性基本上可以保證；除非在某些情形，Jacobian矩陣接近奇異，或者出現(xiàn)某些不連續(xù)現(xiàn)象.而BDF族方法一般步長較大，因此擬Newton迭代不收斂，或者收斂速度太慢的問題也更加嚴(yán)重.作者建議如果在5步以內(nèi)擬Newton迭代不收斂的話，最好重新計算Jacobian矩陣.

另一方面，由于一般動力學(xué)系統(tǒng)的質(zhì)量矩陣、剛度矩陣和約束矩陣等隨時間變化比較慢.因此，如果步長變化不大，上一步用到的Jacobian矩陣就可以在下一步被重復(fù)使用 .作者經(jīng)過比較［8，56］，建議采用如下判據(jù)：如果新步長與原步長相比，

則可以重復(fù)使用Jacobian矩陣，并采用松弛因子

最后，即使擬Newton迭代本身收斂，但得到的修正更新值也可能不滿足預(yù)設(shè)的誤差條件.總體說來，誤差判斷失敗的原因可能有以下幾點：

（1）動力學(xué)中遇到非光滑事件（接觸碰撞、機構(gòu)奇點、參數(shù)或模型突變等），導(dǎo)致狀態(tài)預(yù)測值不準(zhǔn)確；

（2）由于參數(shù)變化，用上一步仿真的后驗誤差估計得到的時間步長h太大，不滿足誤差限要求；

（3）系統(tǒng)剛性太強，或者積分步長太小，導(dǎo)致Jacobian矩陣條件數(shù)很高，舍入誤差太大.

其解決方案就是用得到的后驗修正更新量重新設(shè)置迭代步長，并在必要的情況下降低積分器階數(shù).這一過程和變步長變階策略類似，只是不會用到升階，因而不用進(jìn)行k+1階的誤差估計.

5.2 數(shù)值奇點檢測與處理

數(shù)值奇點是仿真中經(jīng)常遇到的問題.與物理奇異點不同，數(shù)值奇點指的是因為數(shù)值原因?qū)е碌钠娈愇恍?，例如用歐拉角等參數(shù)描述旋轉(zhuǎn)矩陣時遇到的奇點.數(shù)值奇異性的處理流程為：

（1）奇點檢測.對于每種可能遇到的奇異性都要設(shè)置接近于奇點處的報警功能；

（2）奇點消除技術(shù).選擇新的坐標(biāo)系統(tǒng)將數(shù)值奇點處的位形變換到新坐標(biāo)系統(tǒng)下的非奇異位形；

（3）積分器重啟技術(shù).在新的坐標(biāo)下重啟積分器.

數(shù)值奇點檢測與處理的缺點在于隨著部件數(shù)量的增加，不但自由度個數(shù)會增加，而且奇點出現(xiàn)的機率也會增加，而積分器重啟的可能性也會增加.自由度個數(shù)增加，本來計算量就會增加；而積分器重啟的機率也會增加，給仿真帶來的效率損失也會增加.這種損失并非物理系統(tǒng)本身的特性導(dǎo)致，而是由于所選取坐標(biāo)系統(tǒng)的數(shù)值特性不合適導(dǎo)致，只能用數(shù)值方法來解決.

數(shù)值奇點附近光滑性保持問題對于單步法比較簡單的，例如對于廣義 α方法，只要在奇點附近換一套局部非奇異坐標(biāo)，換算出新坐標(biāo)下的q、q?、q?和a就可以直接開始下一步的計算.積分步長只是局部發(fā)生了一些變化，而不會影響整體的仿真效果.而對于多步法，情況比較復(fù)雜.一種思路是總是從最低階（1階）用很小的步長重啟并開始計算，其缺點是重啟過程太慢，一般需要10-20步才能恢復(fù)正常仿真階數(shù)和步長.對于很多部件組成的系統(tǒng)，奇點出現(xiàn)機會可能比較高，反復(fù)重啟可能速度太慢，影響仿真效率.另一種思路是在奇點嚴(yán)重影響仿真效率之前，盡早進(jìn)行坐標(biāo)變換，并將前幾步得到的結(jié)果也進(jìn)行同樣變換，保持高階積分格式.

5.3 不光滑性檢測與處理

不光滑性出現(xiàn)一般有如下幾個原因：接觸碰撞、模型拓?fù)渥兓?、不連續(xù)外力作用等等.其具體表現(xiàn)為某些物理量或者其時間導(dǎo)數(shù)、高階導(dǎo)數(shù)的不連續(xù)性.不連續(xù)性問題在很多研究者看來似乎不甚重要，因為很多自適應(yīng)步長的積分器往往可以通過多次試錯突破不連續(xù)性的障礙，但其代價一般還是比較大的.好在對于絕大部分問題來說，不連續(xù)點的位置比較少，因此這些代價相當(dāng)于整個計算成本來說還是比較小.但事實上，確實有不少工業(yè)界遇到的問題，是由于不連續(xù)性太強而求解失敗.因此作為商業(yè)軟件級別的積分器，還是應(yīng)該包含不光滑性自動檢測和定位的功能，從而可以保證更好的效率和魯棒性.Gear等［57］提出一種ODE的不連續(xù)的檢測和定位方法，可以用來將積分器的時間節(jié)點置于定位得到的不連續(xù)點，而不連續(xù)性問題可以有效地通過積分器重啟得到解決.這一處理不連續(xù)性的方法也可以推廣到DAE中去［58］，并可利用以下性質(zhì)：

（1）質(zhì)量矩陣M隨時間的變化一般緩慢光滑，而外力項Q可能有一些不連續(xù)現(xiàn)象，而Q的不光滑性檢測可以給下一步的預(yù)測提供信息；

（2）除了在有限個間斷點之外，約束C（q，ts）都是至少二階連續(xù).事實上，由于在所有計算過程中，都希望做到C≈0、C?≈0以及C?≈0，因此不光滑性測試不能直接作用在它們上面，而選取實際中最關(guān)鍵的矩陣Cq中的非零元素作為不光滑性測試的量.

6 顯式積分器簡介

一般說來，顯式積分器不適合用于求解代數(shù)微分動力學(xué)系統(tǒng)，但少數(shù)情形例外.在用遞歸方法得到的動力學(xué)方程組中，約束方程的個數(shù)遠(yuǎn)少于自由度個數(shù)；同時，如果系統(tǒng)的剛性問題不太嚴(yán)重，或者高頻振動不能有效地得到衰減（例如對于大規(guī)模的接觸碰撞問題的仿真），此時顯式積分器的計算性能可能會優(yōu)于隱式積分器.另外，顯式積分器的計算量基本可控，因此比較適用于實時仿真.還有，由于顯式積分器無需迭代，在多積分器合作仿真中實現(xiàn)起來方便得多.

約束修正法

從形式上看，約束修正法有些類似直接求解index-3的微分代數(shù)方程組.在這種意義下，大部分顯式常微分方程積分器都可以用來求解方程組（1）.首先，將方程組（1）改寫成常微分方程組

其中γ和方程（5）中的一樣.直接求解方程（88）會引起約束違約現(xiàn)象，因為公式（88）的第二個方程并不等價于公式（1）中的約束方程.方程（88）是常微分方程組，可以用高階顯式的常微分方程積分器來從tn時刻到tn+1時刻進(jìn)行計算，所用積分器既可以是單步法又可以是多步法.計算得到的解不一定滿足約束，因此需要進(jìn)行修正，將每一步求得的q和q?在約束流形上投影.顯式積分器求解方程（88）的計算過程需要用到

在多體系統(tǒng)動力學(xué)中，由于M（q）和Cq（q，t）變化緩慢.因而J也是時間緩變的，可以采用各種近似和復(fù)用來提高計算效率.求解常微分方程的計算過程可以用算子在形式上表示為

和速度修正

這兩種修正計算都需要通過迭代完成，計算中可以復(fù)用方程（89）中的J，以節(jié)省計算量.

罰函數(shù)法

形式上，罰函數(shù)法有些類似求解index-2的微分代數(shù)方程組（81），只是做了一些修改

罰函數(shù)的意思是利用很大的正數(shù)χ（罰因子）將index-2的微分代數(shù)方程組中的約束方程組

中形如f=0的約束改寫成方程f?+χf=0的形式，得到

可以從中求解出q?（q，v，t）和v?（q，v，t），而且求解過程用到方程（89）中的J.得到的常微分方程組可以用顯式積分器進(jìn)行仿真，當(dāng)然并非所有常微分方程積分器都適合開展罰函數(shù)法，有一類特殊的Runge-Kutta-Chebyshev類型的積分器［39-40］，其在負(fù)半實軸上穩(wěn)定域很長（見圖1），因此適合用來顯式求解方程組（91）和（92）.本文推薦采用四階的RKC積分器，其穩(wěn)定區(qū)域是遞歸擴展的［40］，且負(fù)半軸上穩(wěn)定區(qū)域近似與遞歸階的平方成正比.

圖1 RKC積分器的穩(wěn)定區(qū)域示意圖Fig.1 Schematic of the Stable Domain for RKC Integrators

該積分器在虛軸上收斂區(qū)域有限，不適合用來求解剛性問題，但對一般非剛性的多剛體系統(tǒng)動力學(xué)遞歸格式，顯式積分格式求解（91）和（92）的效果非常好.要指出的是，在計算之前可以從系統(tǒng)最大特征值基本確定罰因子χ，因為χ必須遠(yuǎn)遠(yuǎn)大于系統(tǒng)的最大特征值，進(jìn)而可以定出遞歸階.此外，提高計算效率的關(guān)鍵在于如何有效地從微分方程組（91）和（92）中求解出q?和v?，可以根據(jù)所求問題的特點加進(jìn)這一過程.

7 非完整約束DAEs的積分技術(shù)

一般來說，非完整約束的DAEs比完整約束的DAEs簡單，因為對應(yīng)的變量一般是index-2的.帶非完整約束的DAEs的一般格式為［46］

該方程組可以分別被上面介紹的幾種方法離散：

（1）廣義α族方法.同樣假設(shè)q、q?和q?的表達(dá)式為方程（29）和（30），代入方程組（93），得到關(guān)于x、λ和μ的代數(shù)方程組

可以用Newton-Raphson方法或者其它迭代方法求解，其對應(yīng)Jacobian矩陣為

（2）BDF族方法.將BDF方法的公式（74）、（75）和（77）代入方程組（93）得到離散格式

其對應(yīng)Jacobian矩陣為

在求解出位置修正量δ后，可以計算出速度修正量ε；同樣的，ε不能直接用來估計計算誤差，而需要對幾何約束作投影濾波修正.

帶非完整約束的index-2格式的DAEs為

該方程組同樣可以用廣義α族方法和BDF族方法來進(jìn)行積分，可以避免spikes現(xiàn)象，計算出更光滑更高精度的速度、加速度和約束力曲線.

此外，在實際應(yīng)用中，多體動力學(xué)系統(tǒng)往往和控制算法結(jié)合起來應(yīng)用，因此在仿真中也要將控制律方程組同時引入進(jìn)行仿真.控制方程組往往由一階微分方程組和代數(shù)方程組組成，因此可以通過一階微分方程積分器，例如BDF方法等進(jìn)行有效地時間離散.另外，動力學(xué)仿真中往往還包含一些離散動力學(xué)方程組，其積分方法大致有兩種，一種是將離散系統(tǒng)連續(xù)化，另一種是交互仿真技術(shù).前者比較簡單，已經(jīng)普遍應(yīng)用于商業(yè)軟件可以解決自適應(yīng)步長與離散動力學(xué)時間節(jié)點的不協(xié)調(diào)問題；后者需要用到積分器與離散系統(tǒng)動力學(xué)之間的交互仿真技術(shù)，在本文中將不做贅述.

8 典型算例

下面用兩個典型多體系統(tǒng)動力學(xué)測試算例來說明文中討論的一些技術(shù)細(xì)節(jié).應(yīng)該指出的是，不同族的積分器的局部誤差估計在仿真精度控制實踐中有著不同的意義，這既與積分器的階數(shù)有關(guān)又與局部誤差到全局誤差的傳導(dǎo)機制有關(guān).仿真中的整體計算精度其實是由全局誤差來控制，但全局誤差很難估計.大量測試發(fā)現(xiàn)，廣義α族積分器的局部相對誤差限需要比BDF族積分器低兩個數(shù)量級才能達(dá)到類似的整體仿真精度.例如在商業(yè)軟件ADAMS中［59］，HHT 積分器的缺省相對誤差限為10-5，而BDF族積分器的對應(yīng)誤差限為10-3，才能保證二者的整體計算效果大致接近.在以下算例計算中，對廣義 α族積分器，取局部相對誤差限為10-6；而對BDF族積分器，取局部相對誤差限為10-4.計算中，選擇index-3的積分器作為廣義 α族積分器中的代表；選擇變步長公式的index-2和index-3格式的積分器作為BDF族積分器中的代表；而選擇四階顯式積分器作為顯式積分器族中的代表開展計算.

8.1 七剛體機構(gòu)

七剛體模型［60］是多剛體系統(tǒng)動力學(xué)的標(biāo)準(zhǔn)測試程序，其機構(gòu)簡圖為圖2.該機構(gòu)與地面有四個連接點：點O位于原點；點A的坐標(biāo)為（xA，yA）；點B的坐標(biāo)為（xB，yB）；以及點C的坐標(biāo)為（xC，yC）.每個剛體的幾何形狀如圖1所示，其質(zhì)心坐標(biāo)系用箭頭標(biāo)識出來；定常驅(qū)動力矩 τ作用于原點；C點連接的彈簧原長為 l0，其剛度系數(shù)為 k0.這些機構(gòu)的幾何參數(shù)及驅(qū)動力矩見表3，而每個剛體的慣性參數(shù)見表4.

表3 七剛體機構(gòu)的參數(shù)Table 3 Parameters of the seven body mechanism

表4 七剛體機構(gòu)的慣性參數(shù)表Table 4 Inertia parameters of the seven body mechanism

圖2 七剛體機構(gòu)系統(tǒng)Fig.2 Schematic of the seven body mechanism

上述模型分別用廣義α積分器、顯式積分器、BDF格式的index-3積分器進(jìn)行仿真，其計算結(jié)果如圖3所示.

圖3 七剛體機構(gòu)仿真中的曲線圖Fig.3 Diagram for the simulations of the seven body mechanism

仿真過程中各種積分器自適應(yīng)計算的步數(shù)和時間比較見表5.其結(jié)果顯示，對于上述多剛體系統(tǒng)仿真來說，BDF積分器一般要比廣義α積分器快；而對于非剛性問題，顯式積分器也比廣義α積分器快.這主要是由于七剛體模型滿足高階格式的連續(xù)性條件，可以采用大的時間步長進(jìn)行積分.顯式積分器基本不能用于一般的剛性問題求解，因為此時決定其積分步長的是穩(wěn)定性而非精度，而此時穩(wěn)定性條件導(dǎo)致顯式積分器的步長非常小.

表5 七剛體機構(gòu)仿真結(jié)果比較Table 5 The result comparison of the seven body mechanism

圖4給出了計算中BDF積分器的階數(shù)變化，證實了上述模型的仿真中，BDF積分器在整個計算過程中其自適應(yīng)選擇的格式階數(shù)基本都高于二階，這也是其可以采用大的時間步長來計算的主要原因.

圖4 七剛體機構(gòu)仿真中的BDF階數(shù)曲線圖Fig.4 BDF order profiles for the simulations of the Seven Body Mechanism

8.2 曲柄滑塊機構(gòu)

曲柄滑體機構(gòu)模型是測試積分器性能的柔性多體系統(tǒng)動力學(xué)模型，因為實踐表明，三階以上的BDF積分器會在該模型中遇到穩(wěn)定性的考驗.

模型描述如圖5所示，系統(tǒng)由剛性桿1、柔性桿2以及滑塊3組成.剛性桿1的質(zhì)量為0.36 kg，質(zhì)心轉(zhuǎn)動慣量為0.002727 kg?m2，長度l1=0.15 m；柔性桿的密度為7570 kg/m3，截面為8×8 mm2的矩形，楊氏模量為200Gpa，未變形時長度0.30 m；滑塊3的質(zhì)量為0.7555 kg.忽略重力作用，假設(shè)給定旋轉(zhuǎn)驅(qū)動為φ（t）=ωt，其中ω=150 rad/s.計算中取柔性桿前3階橫向振動模態(tài)和第1階縱向振動模態(tài)的模態(tài)坐標(biāo)來描述其變形.假設(shè)機構(gòu)的初始位形為φ=0和x=0.45 m，以及橫向振動模態(tài)坐標(biāo)q1=q2=0，q3=1.033×10-5和縱向振動模態(tài)坐標(biāo)q4=1.69×10-5；速度初始值由φ?=ω可以由DAE初始條件分析算出其它廣義初速度.

圖5 曲柄-滑體機構(gòu)系統(tǒng)Fig.5 Schematic of the crank slider mechanism

該模型具有很強的剛性，只能用隱式積分器來計算.用廣義 α積分器、BDF格式的Index-2和Index-3積分器的計算結(jié)果如圖6所示.從中可以看出，幾個積分器的結(jié)果基本完全吻合；另外，直接求解index-3方程組得到的加速度曲線圖中發(fā)生了誤差放大的尖刺現(xiàn)象.

圖6 曲柄滑體機構(gòu)仿真中的曲線圖Fig.6 Diagram for the simulations of the crank slider mechanism

由于BDF積分器里使用了自適應(yīng)穩(wěn)定性算法，在各階頻率組合導(dǎo)致高階格式可能會發(fā)生不穩(wěn)定性的場合，程序自適應(yīng)地找到最合適的低階格式（一般為二階格式，如圖7所示）來進(jìn)行仿真，而無需手動修改設(shè)置，這樣做為程序的使用帶來了很大方便.另外，即使在如此不利的情況下，BDF積分器的計算效率也與廣義 α方法大致接近，如表6所示，這也驗證了BDF族積分器的高效性和魯棒性.

表6 曲柄滑塊機構(gòu)用各種積分器仿真的效率比較Table 6 Comparison of simulation efficiency of various integrators for the crank slider mechanism

圖7 曲柄滑塊機構(gòu)仿真中的BDF階數(shù)曲線圖Fig.7 BDF order profiles for the simulations of the crank slider mechanism

9 幾種積分器的性能比較

積分器的性能比較應(yīng)該是在同樣的計算精度下進(jìn)行.在同樣的誤差要求下，index-2的DAE積分器的速度變量計算結(jié)果比index-3積分器的結(jié)果更準(zhǔn)確；而且，index-2的DAE積分器可以解決index-3積分器在加速度和拉格朗日乘子曲線上的尖刺狀誤差.此外，index-1的DAE積分器計算得到的加速度曲線更光滑，而且可以對它們進(jìn)行誤差控制.其次從計算規(guī)模上來說，index-1的方程組規(guī)模大于index-2的方程組，更大于index-3的方程組；而其Jacobian矩陣的規(guī)模和復(fù)雜程度也是同樣順序.因此從計算效率上來說，index-3的DAE積分器大于index-2的DAE積分器，更大于index-1的DAE積分器.最后，從計算魯棒性來說，大部分時候BDF族的index-2格式積分器優(yōu)于index-3格式積分器；而index-3的廣義 α積分器的魯棒性似乎更佳.

從計算速度上來說，一般情況下認(rèn)為在同樣的精度要求下，BDF方法快于廣義α方法.顯式積分器不能求解強剛性問題，而根據(jù)目前的測試結(jié)果，它在弱剛性問題和一些中等剛性問題（剛性<104）的動力學(xué)仿真中比廣義α方法要快.從程序設(shè)計角度來說，顯式積分器最容易編程實現(xiàn)，因為無需迭代求解非線性方程組；廣義α族積分器其次，而且穩(wěn)定性比較好；BDF族積分器的結(jié)構(gòu)最復(fù)雜，因為既需要實現(xiàn)變階變步長，還要特別處理高階格式的穩(wěn)定性問題；特別是對于index-3的BDF族積分器還需要特殊處理速度變量的誤差估計問題.從計算魯棒性來說，顯式積分器強于隱式積分器，而廣義α方法優(yōu)于BDF方法，因為大步長積分可能在連續(xù)性較差的問題中會引起計算失敗.另外，廣義α方法除了幾個簡單的參數(shù)調(diào)節(jié)以外，可以改進(jìn)的地方很少，而BDF積分器中的Dk+2n+1矩陣所含信息量很大，如果能更好地利用可以有效提高其計算性能和穩(wěn)定性.最后，顯式積分器無需迭代求解大規(guī)模非線性方程組，在連續(xù)性較差的問題和并行仿真上應(yīng)用潛力似乎更大.

就BDF方法的兩種公式來說，一般認(rèn)為，準(zhǔn)定步長公式稍快于變步長公式，但由于后者直接采用原始數(shù)據(jù)而無需進(jìn)行插值或者步長轉(zhuǎn)換，在程序設(shè)計方便程度和魯棒性方面優(yōu)于前者.另外，廣義 α方法的index-2格式盡管可以完全滿足速度約束，但在很多實際問題應(yīng)用中計算效率似乎不夠高，值得進(jìn)一步深入研究.

10 展望

微分代數(shù)方程積分器是多體系統(tǒng)動力學(xué)軟件的核心模塊.在上世紀(jì)70至90年代，伴隨多體系統(tǒng)動力學(xué)軟件的開發(fā)熱潮，積分器技術(shù)的發(fā)展進(jìn)入了黃金時期.很多著名的研究團隊，通過大量的理論分析和數(shù)值實驗，提出了很多思路來解決動力學(xué)仿真中的各種數(shù)值問題.這些研究結(jié)果也大量被商業(yè)軟件采用、驗證和改進(jìn).同時，為了使自己的軟件能更有效地求解各種工程問題，各大商業(yè)公司也投入成本，展開了大量的測試研究，從工程實踐中提煉出了很多需求，并在軟件升級中一一解決.目前可以認(rèn)為，在比較有名的通用型軟件中，傳統(tǒng)單個積分器的技術(shù)已經(jīng)比較成熟.

這幾年，李群積分器和辛積分器的研究引起了很多學(xué)者的重視.李群積分器在求解高速旋轉(zhuǎn)問題中有無與倫比的優(yōu)勢，其計算結(jié)果的精確性、算法的高效性和魯棒性得到研究者們的贊賞.辛積分器在長時間仿真的應(yīng)用中性能卓越.完善這些積分器的功能并開發(fā)出高效計算程序是目前正在開展的重要研究課題.

隨著硬件技術(shù)的進(jìn)步、計算規(guī)模的擴大和學(xué)科的交叉，多體系統(tǒng)動力學(xué)中的積分器技術(shù)也面臨著新的挑戰(zhàn)和機遇.積分器的計算復(fù)雜度一般在O（n2）與O（n4）之間，取決于系統(tǒng)的物理特性，其中n為系統(tǒng)的自由度.因此，積分器技術(shù)的一個重要研究方向是如何盡可能地降低通用型大規(guī)模多體動力學(xué)仿真中的計算復(fù)雜度.

本文討論的都是動力學(xué)仿真中的初值問題.在最優(yōu)控制的動力學(xué)系統(tǒng)仿真中，往往生成DAE的邊值問題，即對某些狀態(tài)變量，初始值給定；而對其它狀態(tài)變量，最終值給定.一般多體系統(tǒng)動力學(xué)的邊值問題要比初值問題復(fù)雜得多，其Newton迭代可能不收斂.發(fā)展適合大型多體系統(tǒng)動力學(xué)邊值問題的積分器，也是需要重點研究的方向.

多物理場仿真目前來看是多體系統(tǒng)動力學(xué)的重點研究領(lǐng)域，而工業(yè)需求為積分器的研究提出了新的研究課題.如何將多體系統(tǒng)動力學(xué)中發(fā)展出來的積分器和其他學(xué)科的積分器整合起來求解多物理場問題，并要保證計算過程的精確性、高效性與魯棒性，也是積分器技術(shù)的重要研究方向.

聲明與致謝

多體系統(tǒng)動力學(xué)積分器的開發(fā)需要理論的指導(dǎo)，更需要詳實的細(xì)節(jié)實現(xiàn).在工作和學(xué)習(xí)中，作者閱讀和研究過從Gear的DIFSUB到Petzold等人的DASPK等十幾種ODE和DAE積分器的代碼，從中受益非淺.本文主要以廣義α族和BDF族積分器為對象，并結(jié)合作者實踐中的一些體會，盡量系統(tǒng)化地介紹積分器程序開發(fā)中的一些關(guān)鍵技術(shù).這些技術(shù)也可以用于其他族積分器的實現(xiàn).作者聲明，在本文內(nèi)絕不包含任何非開源代碼的技術(shù)細(xì)節(jié)，以尊重商業(yè)軟件的專利和知識產(chǎn)權(quán).

本文初稿完成后，北京理工大學(xué)胡海巖院士、清華大學(xué)任革學(xué)教授、青島大學(xué)潘振寬教授、美國同事宋培林博士等閱讀了相關(guān)內(nèi)容，并給作者提出了寶貴的修改意見.作者深表感謝，也盡量按照幾位前輩的意見作了相應(yīng)的修改.作者在完成和修改本文過程中，曾就文中很多細(xì)節(jié)與北京理工大學(xué)劉鋮研究員和哈工大魏承副教授反復(fù)交流討論，最終才得以定稿，在此一并致謝.另外，還要感謝我的幾個學(xué)生-楊凱、楊思銘、袁騰飛等，最初本文只是團隊內(nèi)部的簡單討論筆記，各位同學(xué)向作者指出了初稿中的一些含糊不清之處，進(jìn)而幫助作者完善了本文中的一些細(xì)節(jié)內(nèi)容.