ADMAS的理论基础.doc

1、ADMAS 的理论基础ADAMS 利用带拉格朗日乘子的第一类拉格朗日方程导出最大数量坐标的微分代数方程（DAE）。它选取系统内每个刚体质心在惯性参考系中的三个直角坐标和确定刚体方位的三个欧拉角作为笛卡尔广义坐标，用带乘子的拉格朗日第一类方程处理具有多余坐标的完整约束系统或非完整约束系统，导出以笛卡尔广义坐标为变量的动力学方程。全部动力学问题归结为求解一个动力学普遍方程。是最高度的概括。动力学普遍方程存在着致命的弊病：在求多自由度时会遇到极大困难，几乎难以下手。原因：方程中 个笛卡尔坐标不是独立的。对于完整系统只有n3个是独立的。snN3拉格朗日开创的分析力学，就是为克服动力学普遍方程的弱点，解

2、决多自由度、非自由系统的动力学问题而发展完善的。第一类拉格朗日方程：我们引入符号 kzfjyixfriikk对约束方程两边变分 ),.321(0),.(21 stfnk ,.1krfiki 实际这也是虚位移应当满足的约束。引入拉格朗日乘子 将上式两端乘 并对 求和),.2(skk011 iiksniiknisk rfrf 将上式与动力学普遍方程相减，可得 011 iikskiini rfrmF独立坐标有 个，对于不独立坐标，我们可选取适当的 使上式等于s3 k零，从而有01. ikskii rfrmF),.21(ni这就是拉格朗日方程乘子的动力学方程，即第一拉格朗日方程。共有 个未知量，可与个

3、 约束方程联立求解。sn3s固定标架可以用来定义构件的形状、质心位置、作用力和反作用力的作用点、构件之间的连接位置等。浮动标记相对于构件运动，在机械系统的运动分析过程中，有些力和约束需要使用浮动标架来定位。动力学方程的求解速度很大程度上取决于广义坐标的选择。研究刚体在惯性空间中的一般运动时，可以用它的质心标架坐标系确定位置，用质心标架坐标相对地面坐标系的方向余弦矩阵确定方位。为了解析地描述方位，必须规定一组转动广义坐标表示方向余弦矩阵。第一种方法是用方向余弦矩阵本身的元素作为转动广义坐标，但是变量太多，同时还要附加六个约束方程；第二种方法是用欧拉角或卡尔登角作为转动坐标，它的算法规范，缺点是在

4、逆问题中存在奇点，在奇点位置附近数值计算容易出现困难；第三种方法是用欧拉参数作为转动广义坐标，它的变量不太多，由方向余弦计算欧拉角时不存在奇点。ADAMS 软件用刚体 的质心笛卡尔坐标和反映刚体方位的欧拉角作为广义坐iB标，即 。由于采用了不独立的广义坐标，TnTi qqzyxq,.,21系统动力学方程虽然是最大数量，但却是高度稀疏耦合的微分代数方程，适用于稀疏矩阵的方法高效求解。利用 ADAMS 建立机械系统仿真模型时，系统中构件与地面或构件与构件之间存在运动副的联接，这些运动副可以用系统广义坐标表示为代数方程，这里仅考虑完整约束。设表示运动副的约束方程数为 ，则用系统广义坐标矢量nh表示的

5、运动学约束方程组为：(1) 0)(),.(,)(21 TKnhKqq考虑运动学分析，为使系统具有确定运动，要使系统实际自由度为零，为系统施加等于自由度 的驱动约束：)(nhc(2)0),(tqD在一般情况下，驱动约束是系统广义坐标和时间的函数。驱动约束在其集合内部及其与运动学约束合集中必须是独立和相容的，在这种条件下，驱动系统运动学上是确定的，将作确定运动。由式(1)表示的系统运动学约束和式(2) 表示的驱动约束组合成系统所受的全部约束：(3)0),(),(tqtDK式（3）为 个广义坐标的 个非线性方程组，其构成了系统位置方程。ncnc对式(3) 求导，得到速度约束方程：(4)0),(),(

6、),( tqttqq若令 ，则速度方程为：),(t(5),(),( ttq对式(4)求导，可得加速度方程：(6) 0),(),(2),(),(),( tqtttqtq qqq若令 ，则加速度方程为：tt2(7)0),(),(),( ttq矩阵 ，为雅可比矩阵，如果 的维数为 ， 维数为 ，那么 维数q mqnq为 矩阵，其定义为 。在这里 为 ( 个运动学约nmjijiq/)(, ch束， 个驱动约束， 个广义坐标）的方阵。hcncADAMS 运动学方程的求解算法 在 ADAMS 仿真软件中，运动学分析研究零自由度系统的位置、速度、加速度和约束反力，因此只需求解系统的约束方程:(8)0),(n

7、tq运动过程中任意时刻 位置的确定，可有约束方程的 Newton-Raphson 迭代nt法求得：(9)0),(njjqt其中， ，表示第 j 次迭代。jjjq1时刻速度、加速度可以利用线性代数方程的数值方法求解，ADAMS 中nt提供了两种线性代数方程求解方法：CALAHAN 方法（由 Michigan 大学Donald Calahan 教授提出）与 HARWELL 方法（由 HARWELL 的 Ian Duff 教授提出），CALAHAN 方法不能处理冗余约束问题，HARWELL 方法可以处理冗余约束问题，CALAHAN 方法速度较快。(10) tq1(11)2)(1ttqADAMS 动力

8、学方程ADAMS 中用刚体 B 的质心笛卡尔坐标和反应刚体方位的欧拉角作为广义坐标。即 ，令 。Tzyxq, , TTTRqzyxR构件质心参考坐标系与地面坐标系间的坐标变换矩阵为：(12) conAgi sinsin siscoisicocsncoio定义一个欧拉转轴坐标系，该坐标系的三个单位矢量分别为上面三个欧拉转动的轴，因而三个轴并不相互垂直。该坐标系到构件质心坐标系的坐标变换矩阵为:(13)01cossinincoB构件的角速度可以表达为：(14)ADAMS 中引入变量 为角速度在欧拉转轴坐标系分量:e(15)e考虑约束方程，ADAMS 利用带拉格朗日乘子的拉格朗日第一类方程的能量形式

9、得到如下方程:(16)jnijjj qQqTdt 1)(T 为系统广义坐标表达的动能， 为广义坐标， 为在广义坐标 方向j j jq的广义力，最后一项涉及约束方程和拉格朗日乘子表达了在在广义坐标 方j向的约束反力。ADAMS 中进一步引入广义动量：(17)jjqTp简化表达约束反力为:(18)jnijqC1这样方程(16) 可以简化为：(19)jjjjQqTP动能可以进一步表达为：(20)JBMRTT21其中 M 为构件的质量阵，J 为构件在质心坐标系下的惯量阵。将(19)分别表达为移动方向和转动方向有：(21)RRCQqTP(22)其中 , 。MVRdtqTdtPRR)()( 0RqT式可简

10、化为：(23)RCQ，由于 B 中包含欧拉角，为了简化推导，ADAMS 中并JqTP)(没有进一步推导 ,而是将其作为一个变量求解。这样 ADAMS 中每个构件具有如下 15 个变量（而非 12 个）和 15 个方程（而非 12 个）。变量：(24)TTTzyxPRV,方程：(25) eeTRJBPCQqVM集成约束方程 ADAMS 可自动建立系统的动力学方程微分代数方程:(26) ),(00tqufFTPFHq其中，P 为系统的广义动量；H 为外力的坐标转换矩阵。初始条件分析 在进行动力学、静力学分析之前，ADAMS 会自动进行初始条件分析，以便在初始系统模型中各物体的坐标与各种运动学约束之

11、间达成协调，这样可以保证系统满足所有的约束条件。初始条件分析通过求解相应的位置、速度、加速度的目标函数的最小值得到。(1)对初始位置分析，需满足约束最小化问题:Minimize： )()(2100qWqCTSubject to： 为构件广义坐标， 为权重矩阵， 为用户输入的值，如果用户输入的q 0值为精确值，则相应权重较大，并在迭代中变化较小。可以利用拉格朗日乘子上述约束最小化问题变为如下极值问题：(27)TTqWqL)()(2100取最小值，则由 得：L,(28)0)(qT因约束函数中存在广义坐标，该方程为非线性方程须用 Newton-Raphson 迭代求解，迭代方程如下：(29)(00q

12、WqqWTT 对初始速度分析，需满足约束最小化问题Minimize： )()(2100CTSubject to： tq其中， 为用户设定的准确的或近似的初始速度值，或者为程序设定的缺0省速度值； 为 对应的权重系数矩阵。Wq同样可以利用拉格朗日乘子将上述约束最小化问题变为如下极值问题：(30)TT tqqWqL)()(2100取最小值，得：(31)0)(0tqT为已知，该方程为线性方程组可求解如下方程：q(32)tWqqWT00对初始加速度、初始拉氏乘子的分析，可直接由系统动力学方程和系统约束方程的两阶导数确定。ADAMS 动力学方程的求解 对于式(26)微分代数方程的求解，ADAMS 采用两

13、种方式求解，第一种为对 DAE 方程的直接求解，第二种为 DAE 方程利用约束方程将广义坐标分解为独立坐标和非独立坐标然后化简为 ODE 方程求解。 DAE 方程的直接求解将二阶微分方程降阶为一阶微分方程来求解，通过引入，将所有拉格朗日方程均写成一阶微分形式，该方程为 Index 3 微分代数方程。I3 积分格式：(33),(00tqufFTPFHq运用一阶向后差分公式，上述方程组对 求导，可得其 Jacobi 矩阵，),(qu然后利用 Newton-Rapson 求解。可以看出，当积分步长 减小并趋近于 0 时，h上述 Jacobi 矩阵呈现病态。为了有效地监测速度积分的误差，可采用降阶积分

14、方法（Index reduction methods）。通常来说，微分方程的阶数越少，其数值求解稳定性就越好。ADAMS 还采用两种方法来降阶求解，即 SI2(Stabilized-Index Two)和 SI1(Stabilized-Index One)方法。SI2 积分格式：(34),(0),();0tqufFtqTPFH上式能同时满足 和 求解不违约，且当步长 趋近于 0 时，Jacobi 矩阵 h不会呈现病态现象。SI1 积分格式：(35),(0),(0tqufFtqTPFH上式中，为了对方程组降阶，引入 和 来替代拉格朗日乘子，即 ,。这种变化有效地将上述方程组的阶数降为 1。因为只

15、需要微分速度约束方程一次来显示地计算表达式 和 。运用 SI1 积分器，能够方便地监测 , , qu和 的积分误差，系统的加速度也趋向于更加精确。但在处理有明显的摩擦接触问题时，SI1 积分器十分敏感并具有挑剔性。静力学分析 在 进 行 静 力 学 、 准 静 力 学 分 析 时 ， 对 动 力 学 方 程 的 速 度 、 加 速 度 设 置为 零 ， 则 得 到 静 力 学 方 程 如 下 :(36),(00tqufFTPFHq该方程为非线性代数方程利用 Newton-Rapson 迭 代 求 解 求 解 。线性化分析 在 系 统 的 某 点 处 ， 可 对 系 统 的 动 力 学 方 程

16、进 行 线 性 化 ，*,uq(37)KCM0为 常 数 阵 可 对 （ 36） 式 求 解 得 到 系 统 的 频 率 和 振 动 模 态 。KCM,ADAMS 求 解 器 算 法 介 绍 ADAMS 数值算法简介 运 动 学 、 静 力 学 分 析 需 求 解 一 系 列 的 非 线 性 代 数 方 程 、 线 性 代 数 方 程 ，ADAMS 采 用 了 修 正 的 Newton-Raphson 迭 代 算 法 求 解 非 线 性 代 数 方 程 ， 以 及基 于 LU 分 解 的 CALAHAN 方 法 和 HARWELL 方 法 求 解 线 性 代 数 方 程 。 对 动 力 学微

17、分 方 程 ， 根 据 机 械 系 统 特 性 ， 选 择 不 同 的 积 分 算 法 ； 对 刚 性 系 统 ， 采 用 变系 数 的 BDF(Backwards Differentiation Formulation)刚 性 积 分 程 序 ， 它是 自 动 变 阶 、 变 步 长 的 预 估 校 正 法 (PECE,Predict-Evaluate-Correct-Evaluate)， 并 分 别 为 Index3、 SI2、 SI1 积 分 格 式 ， 在 积 分 的 每 一 步 采 用了 修 正 的 Newton-Raphson 迭 代 算 法 ； 对 高 频 系 统 (High-F

18、requencies),采用 坐 标 分 块 法 (Coordinate-Partitioned Equation)将 微 分 代 数 （ DAE）方 程 简 化 为 常 微 分 (ODE)方 程 分 别 利 用 ABAM（ Adams-Bashforth-Adams-Moulton） 方 法 和 龙 格 库 塔 （ RKF45） 方 法 求 解 。在 ADAMS 中 具 体 如 下 ： 线 性 求 解 器 (求 解 线 性 方 程 )， 采 用 稀 疏 矩 阵 技 术 以 提 高 效 率 。 CALAHAN 求 解 器 与 HARWELL 求 解 器 。 非 线 性 求 解 器 (求 解 代

19、 数 方 程 )， 采 用 了 Newton-Raphson 迭 代 算 法 。 DAE 求 解 器 (求 解 微 分 代 数 方 程 )， 采 用 BDF 刚 性 积 分 法 。SI2： GSTIFF、 WSTIFF 与 CONSTANT_BDF。SI1： GSTIFF、 WSTIFF 与 CONSTANT_BDF。I3： GSTIFF、 WSTIFF、 DSTIFF 与 CONSTANT_BDF。 ODE 求 解 器 (求 解 非 刚 性 常 微 分 方 程 ) ABAM 求 解 器 与 RKF45 求 解 器 。动力学求解算法介绍 微分代数（DAE） 方 程 的 求 解 算 法 过 程A

20、DAMS 中 DAE 方程的求解采用了 BDF 刚性积分法，以下为其步骤：(1)预 估 阶 段 用 Gear 预估-校正算法可以有效地求解微分-代数方程。首先，根据当前时刻的系统状态矢量值，用泰勒级数预估下一时刻系统的状态矢量值：(38).!2121 htytynnn其中，时间步长 。这种预估算法得到的新时刻的系统状态矢量值通常nth1不准确，可以由 Gear 阶积分求解程序(或其他向后差分积分程序)来校正。k(39)kiinnnyhy1101其中， 为 在 时的近似值， 和 为 Gear 积分程序的系数值。1n)(t1nt0i上式经过整理，可表示为：(40)1110kiinnnyhy（2）

21、校 正 阶 段求 解 系 统 方 程 ，如 ，则方程成立，此时 的为方程的解，G),(tyy否则继续；求解 Newton-Raphson 线性方程，得到 ，以更新 ，使系统方程 更接yG近于成立。，其中为系统的雅可比矩阵。),(1ntyGJ利用 Newton-Raphson 迭代，更新 :ykky1重 复 以 上 步 骤 直 到 足够小。y（3） 误 差 控 制 阶 段预估计积分误差并与误差精度比较，如积分误差过大则舍弃此步。计 算 优 化 的 步 长 和阶数 。hn如达到仿真结束时间，则停止，否则 ，重新进入第一步。tt坐 标 缩 减 的 微 分 方 程 求 解 过 程 算 法ADAMS 程

22、 序 提 供 ABAM(AdamsB ashforth and Adams-Moulton)和 RKF45积分程序，采 用 坐 标 分 离 算 法 ， 将 微 分 -代 数 方 程 减 缩 成 用 独立 广 义 坐 标 表 示 的 纯 微 分 方 程 ， 然 后 用 ABAM 或 RKF45 程序进行数值积分。以下以 ABAM 为例介绍其求解过程。坐标减缩微分方程的确定及其数值积分过程按以下步骤进行：坐 标 分 离 将 系 统 的 约 束 方 程 进 行 矩 阵 的 满 秩 分 解 ， 可 将 系 统的 广 义 坐 标 列 阵 分 解 成 独 立 坐 标 列 阵 和 非 独 立 坐 标 列qi

23、q阵 ， 即 。dqTdi（1）预 估 用 Adams-Bashforth 显 式 公 式 ， 根 据 独 立 坐 标 前 几 个 时 间 步 长 的值 ， 预 估 时 刻 的 独 立 坐 标 值 ， 表 示 预 估 值 。1nt Piq（2）校正 用 Adams-Moulton 隐 式 公 式 对 上 面 的 预 估 值 ， 根 据 给 定 的 收 敛 误 差 限 进行 校 正 ， 以 得 到 独 立 坐 标 的 校 正 值 ， 表 示 校 正 值 。Ciq（4） 确 定 相 关 坐 标确 定 独 立 坐 标 的 校 正 值 之 后 ， 可 由 相 应 公 式 计 算 出 非 独 立 坐 标 和 其 他 系统 状 态 变 量 值 。（5） 积 分 误 差 控 制与 上 面 预 估 校 正 算 法 积 分 误 差 控 制 过 程 相 同 ， 如 果 预 估 值 与 校 正 值 的差 值 小 于 给 定 的 积 分 误 差 限 ， 接 受 该 解 ， 进 行 下 一 时 刻 的 求 解 。 否 则 减 小 积分 步 长 ， 重 复 第 二 步 开 始 的 预 估 步 骤 。