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）都较快收敛。这个例子说明数值刚性系统必须采用专用于求解刚性问题的数值方法。