1.7 ADAMS求解器算法介紹

1.7.1 ADAMS數值算法簡介

運動學、靜力學分析需求解一系列的非線性和線性代數方程。ADAMS采用了修正的Newton-Raphson迭代算法求解非線性代數方程，以及基于LU分解的CALAHAN方法和HARWELL方法求解線性代數方程。

對動力學微分方程，根據機械系統特性，選擇不同的積分算法；對剛性系統，采用變系數的BDF（Backwards Differentiation Formulation）剛性積分程序，它是自動變階、變步長的預估校正法（Predict-Evaluate-Correct-Evaluate, PECE），并分別為I3、SI2、SI1積分格式，在積分的每一步采用了修正的Newton-Raphson迭代算法。

對高頻系統（High-Frequencies），采用坐標分塊法（Coordinate-Partitioned Equation）將微分-代數（DAE）方程簡化為常微分（ODE）方程，并分別利用ABAM（ADAMS-Bashforth-ADAMS-Moulton）方法和龍格-庫塔（RKF45）方法求解。

在ADAMS中，具體方法如下：

● 線性求解器（求解線性方程），采用稀疏矩陣技術，以提高效率。

? CALAHAN求解器。

? HARWELL求解器。

● 非線性求解器(求解代數方程)，采用Newton-Raphson迭代算法。

● DAE求解器（求解微分-代數方程），采用BDF剛性積分法。

? SI2:GSTIFF、WSTIFF與CONSTANT_BDF。

? SI1: GSTIFF、WSTIFF與CONSTANT_BDF。

? I3:GSTIFF、WSTIFF、 Dstiff與CONSTANT_BDF。

● ODE求解器（求解非剛性常微分方程）。

? ABAM求解器。

? RKF45求解器。

1.7.2 動力學求解算法介紹

微分-代數（DAE）方程的求解算法過程ADAMS中DAE方程的求解采用了BDF剛性積分法，有如下幾個操作步驟。

（1）預估階段

用Gear預估-校正算法有效地求解微分-代數方程。首先，根據當前時刻的系統狀態矢量值用泰勒級數預估下一時刻系統的狀態矢量值：

其中，時間步長h=tn 1+-tn。這種預估算法得到的新時刻的系統狀態矢量值通常不準確，由Gear k+1階積分求解程序（或其他向后差分積分程序）來校正。

其中，yn 1+ 為y(t)在t=tn+1時的近似值，β0和αi為Gear積分程序的系數值。

式（1-60）經過整理，可表示為：

（2）校正階段

求解系統方程G，若G( y,,t)=0，則方程成立，此時的y為方程的解，否則繼續。

求解Newton-Raphson線性方程，得到Δy，以更新y，使系統方程G更接近于成立。

JΔy=G( y,,tn+1)

其中，J為系統的雅可比矩陣。

利用Newton-Raphson迭代更新y：

yk 1+ =yk+Δyk

重復以上步驟，直到Δy足夠小。

（3）誤差控制階段

預估積分誤差并與誤差精度進行比較，若積分誤差過大則舍棄此步。

計算優化的步長h和階數n。

若達到仿真結束時間則停止，否則t=t+Δt，重新進入第一步。

1.7.3 坐標縮減的微分方程求解過程算法

ADAMS程序提供ABAM（ADAMS-Bashforth and ADAMS-Moulton）和RKF45積分程序，采用坐標分離算法，將微分-代數方程減縮成用獨立廣義坐標表示的純微分方程，然后用ABAM或RKF45程序進行數值積分。下面以ABAM為例介紹其求解過程。

坐標減縮微分方程的確定及其數值積分過程按以下步驟進行。

Step01 坐標分離。將系統的約束方程進行矩陣的滿秩分解，可將系統的廣義坐標列陣{q}分解成獨立坐標列陣{qi}和非獨立坐標列陣{qd}，即{q}=。

Step02 預估。用ADAMS-Bashforth顯式公式根據獨立坐標前幾個時間步長的值預估tn 1+時刻的P獨立坐標值，P表示預估值。

Step03 校正。用ADAMS-Moulton隱式公式對上面的預估值根據給定的收斂誤差限進行校正，以得到獨立坐標的校正值，C表示校正值。

Step04 確定相關坐標。確定獨立坐標的校正值之后可由相應公式計算出非獨立坐標和其他系統狀態變量值。

Step05 積分誤差控制。與上面預估-校正算法積分誤差控制過程相同，如果預估值與校正值的差值小于給定的積分誤差限就接受該解，進行下一時刻的求解；否則減小積分步長，重復第二步開始的預估步驟。

1.7.4 動力學求解算法特性比較

（1）I3積分格式僅監控位移和其他微分方程的狀態變量的誤差。當積分步長變小時Jacobian矩陣不能保持穩定，會出現奇異，積分易發散。積分過程不能監控速度和約束反力，因而速度、加速度、約束反力計算精度差一些。

（2）SI2積分格式中考慮了速度約束方程，控制拉氏乘子的誤差、速度誤差，仿真結果更精確，給出速度、加速度較為精確的解。Jacobian矩陣在步長很小時仍能保持穩定。Jacobian矩陣小步長不會奇異、病態，增加了校正器在小步長時的穩定性和魯棒性。校正階段不會像I3積分格式那樣容易失敗。

SI2積分格式精確處理高頻問題，但比I3積分格式慢，驅動約束為速度時，輸入必須可微、光滑。非光滑驅動約束運動輸入會產生無限加速度，從而導致SI2積分失敗。位移驅動約束輸入不能是變量的函數，速度、加速度輸入是變量的函數，而I3驅動約束輸入是變量的函數，會給仿真帶來不便。

（3）SI1積分格式中考慮了速度約束方程，但并沒有引入加速度約束方程，相對應引入了拉氏乘子的導數而使方程降階，控制拉氏乘子的誤差、速度誤差，仿真結果很精確，Jacobian矩陣在步長很小時仍能保持穩定，增加了校正器在小步長時的穩定性和魯棒性。

SI1積分格式給出速度、加速度較為精確的解，監控所有狀態變量，比如位移、速度和拉氏乘子，比SI2精度高，但對具有摩擦、接觸的模型很敏感。

上述3種積分方式的比較如表1-1所示。

1.7.5 求解器的特點比較

（1）Gstiff求解器為剛性穩定算法，采用多步、變階（最高階為6）、變步長、固定系數算法，可直接求解DAE方程，有I3、SI2、SI1三種積分格式。在預估中采用泰勒級數，而且其系數是假設步長不變而得到的固定系數，因而當步長改變時會產生誤差。其特點是計算速度快，位移精度高，I3格式時速度（尤其是加速度）會產生誤差，通過最大步長來控制求解中步長的變化，從而提高精度，使仿真運行在定步長狀態。當步長小時，Jacobian矩陣是步長倒數的函數會變成病態，SI2及SI1積分格式時Jacobian矩陣步長很小時仍能保持穩定。該算法適用于很多仿真分析問題。

（2）Wstiff求解器為剛性穩定算法，采用多步、變階（最高階為6）、變步長、變系數算法，可直接求解DAE方程，有I3、SI2、SI1三種積分格式。在預估中采用NDF（Newton Divided Difference）公式，根據步長信息修改相應階的系數，而且步長改變并不影響精度，因而更具健壯性、更穩定，但仿真時間比Gstiff長。

（3）Dstiff求解器為剛性穩定算法，采用多步、變階（最高階為6）、變步長、變系數（固定第一個系數）算法，可直接求解DAE方程，ADAMS中僅有I3一種積分格式。在預估中采用NDF（Newton Divided Difference）公式，固定第一個系數，從而使第一個系數與步長無關，其他可變系數隨步長變化而變化，根據步長信息修改相應階的系數，比較穩定，但仿真時間比Gstiff長。Dstiff求解器基于DASSL積分器，是由Petzold開發的。

（4）Constant_BDF求解器為剛性穩定算法，采用多步、變階（最高階為6）、固定步長算法，可直接求解DAE方程，有I3、SI2、SI1三種積分格式。在預估中采用NDF（Newton Divided Difference）公式，在SI2積分格式時小步長非常穩定健壯，可解Gstiff失敗的問題，位移、速度求解精度高，而且對加速度和力的不連續性沒有Gstiff求解器敏感，有些問題沒有Gstiff、Wstiff快，Hmax太大會導致結果不準、太小則速度太慢。

（5）ABAM求解器為非剛性穩定算法，采用多步、變階算法（最高階為12）、變步長算法，適合求解低阻尼、瞬態系統，尤其適合求解非剛性系統但存在突變或高頻的系統。ABAM利用坐標分塊技術將DAE方程變為ODE方程，僅獨立坐標被積分求解，其他非獨立坐標利用約束方程（代數方程）求解。ABAM求解器是由L.F.Shampine和M.K.Gordon開發的。

（6）RKF45非剛性穩定算法采用單步算法，是以上多步算法的補充，但在積分計算時計算導數費時，而且與其他算法相比不能給出高精度結果，且速度比ABAM積分器慢。L.FShampine H.A.Watts開發了DDERKF積分器。

1.7.6 剛性問題求解算法選擇

數值剛性問題是指系統的特征值分布廣泛，存在低頻、高頻，而且對應的高頻部分具有較高阻尼，因而當系統有可能高頻振動時高頻阻尼會使其很快散掉。剛度比為系統隱藏的最高頻率（對應較高阻尼）與系統表現出的最低頻率（對應較低阻尼）的比值。

一般剛度比為200時為稱為剛性系統，剛度比為20以下時稱為非剛性系統。非剛性系統的最高頻率一定對應較小阻尼而被激發出。例如，具有柔性體的系統，柔體的高頻都具有高阻尼，一般不會被激發，都是低頻被激發，系統的高頻被激發時系統則變為非剛性系統。剛性積分器對數值剛性系統的微分方程進行有效的積分，剛性積分器中積分器步長被限制為最高主動頻率（系統表現出的最高頻率）的倒數，而非剛性積分器中積分器步長被限制為最高頻率（系統所有頻率中的最高頻率，包含隱藏頻率）的倒數，這樣非剛性積分器對數值剛性系統的微分方程積分的效率非常低。

在ADAMS中，如果一個系統是非數值剛性系統，就采用ABAM或RK45積分器，也采用Gstiff、Wstiff、Dstiff、Constant_BDF積分器；如果系統是數值剛性系統，采用了ABAM或RK45，系統就不會收斂或計算速度奇慢。

數值剛性系統除在剛度方面存在較大差異外，還有一個很重要的特征是對應高頻的阻尼較大，使較高頻率被基本阻尼掉，而低頻則處于未阻尼狀態。當數值剛性系統采用ADAMS非剛性數值算法（如ABAM或RKF45）時會出現數值困難，很難收斂，而用剛性數值算法（如Gstiff、Wstiff、Dstiff或Constant_BDF）時則很快收斂。去掉阻尼后的物理剛性系統若高頻沒有被阻尼掉則為高頻系統，采用非剛性數值算法（如ABAM或RKF45）以及剛性數值算法（如Gstiff、Wstiff、Dstiff或Constant_BDF）都較快收斂。這個例子說明數值剛性系統必須采用專用于求解剛性問題的數值方法。