ADAMS软件算法基本理论陈立平机械系统动力学分概要Word下载.docx

1、软件用刚体Bi的质心笛卡尔坐标和反映刚体方位的欧拉角作为广义坐标，即 qi二x,y,zi,Y f，q =q1，q2,，q【r。由于采用了不独立的广义坐标，系统动力学方程虽然是最大数量，但却是高度稀疏耦合的微分代数方程，适用于稀疏矩阵的方法高 效求解。4. 2 ADAMS运动学分析4.2.1ADAMS运动学方程利用ADAMS建立机械系统仿真模型时，系统中构件与地面或构件与构件之间存在运 动副的联接，这些运动副可以用系统广义坐标表示为代数方程，这里仅考虑完整约束。设 表示运动副的约束方程数为 nh，则用系统广义坐标矢量表示的运动学约束方程组为：K（q）珂皆（q）,站（q）,，吒（q）T =0 （4

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

3、q,t） =0 （4-4）若令 t （q,t），则速度方程为：G（q,q,t）二：i（q,t）q-吟-0 （4-5）对式（4-4）求导，可得加速度方程：$（q,q,q；t） =6q（q,t）q +（ q（q,t）q）qq + 2 s（q,t）q + tt（q,t） =0 （4-6）若令 =C：qq）qq -2Gqtqtt，则加速度方程为：J（q,q,q,t） q（q,t）q （q, q, t） 0 （4-7）矩阵Gq,为雅可比矩阵，如果 门的维数为m, q维数为n,那么G q维数为m n矩 阵，其定义为（q）（i,j）= /闵）。在这里q为ncFc （nh个运动学约束，nc nh个 驱动约束，

4、nc个广义坐标）的方阵。4.2.2ADAMS运动学方程的求解算法在ADAMS仿真软件中，运动学分析研究零自由度系统的位置、速度、加速度和约束 反力，因此只需求解系统的约束方程：G（q,tn） =0 （4-8）运动过程中任一时刻tn位置的确定，可由约束方程的 Newton-Raphson迭代法求得：叫j.q ：（qj,tj =0 （4-9）其中，q =qj 4 -qj，表示第j次迭代。tn时刻速度、加速度可以利用线性代数方程的数值方法求解， ADAMS中提供了两种线性代数方程求解方法： CALAHAN 方法（由Michigan大学 Donald Calahan教授提出） 与HARWELL方法（由

5、HARWELL 的Ian Duff教授提出 ），CALAHAN 方法不能处理冗 余约束问题，HARWELL方法可以处理冗余约束问题， CALAHAN方法速度较快。q =-叮 t （4-10）dW-Gq，_C：-qd）qq W tt （4-11）4. 3 ADAMS动力学分析4.3.1ADAMS动力学方程ADAMS 中用刚体 B的质心笛卡尔坐标和反映刚体方位的欧拉角作为广义坐标，即q 二x,y,乙门，T，令 R Jx,y,z， - i ，q =RT, TT。构件质心参考坐标系与地面坐标系间的坐标变换矩阵为：cos- cos sin cost sin cos- sin sincost cos si

6、n- sin |Agi = sin即 cos + cos即 cos日 sin sin即 sin 0+cos即 cos日 cos cos即 sin日 sin 日 sin sin 日 cos cos 日（4-12）定义一个欧拉转轴坐标系，该坐标系的三个单位矢量分别为上面三个欧拉转动的轴， 因而三个轴并不相互垂直。该坐标系到构件质心坐标系的坐标变换矩阵为：sin）sin 0 cosB= sin 日cos 0 -sin 日 （4-13）cos e 1 o _构件的角速度可以表达为：B （ 4-14）ADAMS中引入变量-e为角速度在欧拉转轴坐标系分量:（4-15）考虑约束方程，ADAMS利用带拉格朗日

7、乘子的拉格朗日第一类方程的能量形式得到 如下方程：（4-16）dt （中y二dt ；q i 1 :qT为系统广义坐标表达的动能， qj为广义坐标，Qj为在广义坐标qj方向的广义力,qj方向的约束反力。（4-17）最后一项涉及约束方程和拉格朗日乘子表达了在在广义坐标ADAMS中近一步引入广义动量：简化表达约束反力为:这样方程（4-16）可以简化为:动能可以近一步表达为：T JRTmR 1 |TBTJB J2 2其中M为构件的质量阵，J为构件在质心坐标系下的惯量阵。 将（4-19）分别表达为移动方向与转动方向有：=Qr _Cr其中 pR = d t d mR = mV ， = 0。R dtl/MR

8、 丿 dt 喩（4-21）式可以简化为：mV = Qr -cr待 TP = =BTJB ,，由于B中包含欧拉角，为了简化推导， ADAMS回Y丿步推导P，而是将其作一个变量求解。（4-18）（4-19）（4-20）（4-21）（4-22）（4-23）中并没有进这样ADAMS中每个构件具有如下 15个变量（而非12个）和15个方程（而非12个）。 变量：厂 一 JV=M，Vy，Vz（4-24）R = lx, y, zf 宀尸冋P&P点 COe =灼屮召灼Y=即月冲T方程:其中，P为系统的广义动量；H为外力的坐标转换矩阵。为了更好地说明ADAMS的建模过程下面以一个单摆为例进行建模推导。图4-1单

9、摆示意图如图所示，单摆的质量为 M、惯量为I,杆长为2L,并在0点以转动副与大地相连接约 束在大地的 OXY平面内。在单摆质心处建立单摆的跟随坐标系局部构件参考坐标系Op Xp Yp,其坐标在地面坐标系 OXY中为（x, y）,单摆的姿态角为 0。 系统的动能表达式：T MX2 My2 （4-27）广义动量表达式：% MX外力表达式:cX（4-28 ）（4-29 ）（4-30 ）0 1 HTF = Mg0 j约束方程：x- Leos： - 0y - Lsin）- 0约束方程的雅克比矩阵:（4-31）Lsin vL cosB约束对应的拉格朗日乘子:（4-32）力、力矩平衡方程：一 mvL 十打=

10、o 1p） + qT?“+H TF =0n MVy +入2 -Mg =0 （4-33）L祐十人L si n日一九2 L cos日=0动量矩表达式：P-汀妒P = I .运动学关系方程：Vx -殳=0 1U=d= Vy_y=0Qg -d = 0其方程集成表达为：一 mVx+入=0 IMVy + 打 _Mg =0P + 人L sin 日一毎 L cos日=0I PT日=。Vx - X=0I Vy-y=0co 0 tl =0x - L cos日=0y Lsi n 日=0其中系统需求解变量为：x y 二 Vx Vy r Pj 1 2（4-34）（4-35）（4-36）（4-37）432初始条件分析在进

11、行动力学、静力学分析之前， ADAMS会自动进行初始条件分析，以便在初始系统模型中各物体的坐标与各种运动学约束之间达成协调，这样可以保证系统满足所有的约 束条件。初始条件分析通过求解相应的位置、速度、加速度的目标函数的最小值得到。（1 ）对初始位置分析，需满足约束最小化问题1 TMinimize : C q - q0 W（q -q0）Subject to : ：j q =q为构件广义坐标，W为权重矩阵，qo为用户输入的值，如果用户输入的值为精确值, 则相应权重较大，并在迭代中变化较小。可以利用拉格朗日乘子将上述约束最小化问题变 为如下极值问题：1 . T ” TLSq-q 570+ 9 （4-

12、38）闩 FIL取最小值，则由 0, . 0得:cq cK其中，q为用户设定的准确的或近似的初始速度值，或者为程序设定的缺省速度值;w为对应q0的权重系数矩阵。同样可以利用拉格朗日乘子将上述约束最小化问题变为如下极值问题:q为已知，该方程为线性方程组可求解如下方程:（4-43）（3）对初始加速度、初始拉氏乘子的分析，可直接由系统动力学方程和系统约束方程 的两阶导数确定。4.3.3ADAMS动力学方程的求解对于式（4-26）微分代数方程的求解，ADAMS采用两种方式求解，第一种为对 DAE方程的直接求解，第二种为DAE方程利用约束方程将广义坐标分解为独立坐标和非独立坐 标然后化简为ODE方程求解

13、。关于具体求解器将在 4.5节介绍。DAE方程的直接求解将二 阶微分方程降阶为一阶微分方程来求解，通过引入 u = q，将所有拉格朗日方程均写成一阶微分形式，该方程为 Index 3微分代数方程。I3积分格式：运用一阶向后差分公式，上述方程组对 （u q ）求导，可得其Jacobian矩阵,然后利用 Newt on-Raps on求解。可以看出，当积分步长 h减小并趋近于0时，上述Jacobian矩阵呈现病态。为了有效地监测速度积分的误差，可采用降阶积分方法（ In dexreduction methods ）。通常来说，微分方程的阶数越少，其数值求解稳定性就越好。ADAMS 还采用两种方法来

14、降阶求解，即 SI2（Stabilized-lndex Two）和 SI1（Stabilized-lndex One）方法。SI2积分格式：（4-45）q,t =0i q,u,t =0F 二 f （u,q,t）上式能同时满足 和去求解不违约，且当步长 h趋近于0时，Jacobian矩阵不会呈现病态现象。SI1积分格式:上式中，为了对方程组降阶，引入 .和i来替代拉格朗日乘子，即 二，i =I。这种变化有效地将上述方程组的阶数降为 1。因为只需要微分速度约束方程一次来显示地计算表达式；和i。运用SI1积分器，能够方便地监测 q , u, n和Z的积分误差，系统的 加速度也趋向于更加精确。 但在处

15、理有明显的摩擦接触问题时， SI1积分器十分敏感并具有挑剔性。4. 4 ADAMS静力学及线性化分析441静力学分析在进行静力学、准静力学分析时，对动力学方程的速度、加速度设置为零，则得到静（4-47）力学方程如下:u =0G q,t =0F = f （u,q,t）4.4.2线性化分析在系统的某点处，q=q”,u=u”可对系统的动力学方程进行线性化,（4-48）MU Cu Kq = 0 u =qM，C，K为常数阵可对（4 .4 - 1）式求解得到系统的频率和振动模态。4. 5 ADAMS求解器算法介绍4.5.1ADAMS数值算法简介运动学、静力学分析需求解一系列的非线性代数方程、线性代数方程，

16、 ADAMS采用了修正的Newton-Raphson迭代算法求解非线性代数方程， 以及基于LU分解的CALAHAN 方法和HARWELL方法求解线性代数方程。对动力学微分方程，根据机械系统特性，选择不同的积分算法；对刚性系统，采用变系数的 BDF（Backwards Differentiation Formulation）刚性积分程序，它是自动变阶、变步长的预估校正法（PECE , Predict-Evaluate-Correct-Evaluate），并分别为 Index3、SI2、SI1 积分格式，在积分的每一步 采用了修正的Newton-Raphson迭代算法；对高频系统（High-Fre

17、quencies），采用坐标分块法 （Coordinate-Partitioned Equation）将微分代数（DAE ）方程简化为常微分（ ODE）方程分 别利用 ABAM （Adams-Bashforth-Adams-Moulton ）方法和龙格库塔 （RKF45）方法求解。在ADAMS中具体如下：线性求解器（求解线性方程），采用稀疏矩阵技术以提高效率。CALAHAN 求解器与HARWELL求解器。非线性求解器（求解代数方程），采用了 Newton-Raphson迭代算法。DAE求解器（求解微分代数方程），采用BDF刚性积分法。SI2: GSTIFF、WSTIFF 与 CONSTANT_

18、BDF。SI1 :I3 : GSTIFF、WSTIFF、 Dstiff 与 CONSTANT_BDF。ODE求解器（求解非刚性常微分方程）ABAM求解器与RKF45求解器。4.5.2动力学求解算法介绍1微分代数（DAE ）方程的求解算法过程ADAMS中DAE方程的求解采用了 BDF刚性积分法，以下为其步骤：（1 ）预估阶段用Gear预估-校正算法可以有效地求解微分 -代数方程。首先，根据当前时刻的系统状态矢量值，用泰勒级数预估下一时刻系统的状态矢量值：yn 1 = yn n h 川 （4-49）Ct 2! Ct其中，时间步长h二tn 1 -tn。这种预估算法得到的新时刻的系统状态矢量值通常不准

19、确，可以由Geark 1阶积分求解程序（或其他向后差分积分程序）来校正。kyn 1 二-h Tn7 ryn4 1 （4-50 ）i ：其中，yn 1为y（t）在t =tn1时的近似值； r和J为Gear积分程序的系数值。_1 k上式经过整理，可表示为：yn 1 yn1 二 （4-51）i 二（2 ）校正阶段求解系统方程 G，如G（y,y）t） & ，则方程成立，此时的y为方程的解，否则继续;求解Newton-Raphson线性方程，得到二y，以更新y ,使系统方程G更接近于成立。J y二G（y, y,tn J，其中J为系统的雅可比矩阵。k -J1 k k利用 Newt on-R aphson

20、迭代，更新 y: y = y =y重复以上步骤直到 Ay足够小。（3）误差控制阶段预估计积分误差并与误差精度比较，如积分误差过大则舍弃此步。计算优化的步长h和阶数n。如达到仿真结束时间，则停止，否则 t ：t，重新进入第一步。2坐标缩减的微分方程求解过程算法ADAMS 程序提供 ABAM（Adams Bashforth and Adams-Moulton）和 RKF45 积分程序， 采用坐标分离算法，将微分 -代数方程减缩成用独立广义坐标表示的纯微分方程，然后用ABAM或RKF45程序进行数值积分。以下以 ABAM为例介绍其求解过程。坐标减缩微分方程的确定及其数值积分过程按以下步骤进行：（1

21、）坐标分离将系统的约束方程进行矩阵的满秩分解， 可将系统的广义坐标列阵 分解成独立坐标列阵tq和非独立坐标列阵 gd，即 g j = g qd j。（2）预估用Adams-Bashforth显式公式，根据独立坐标前几个时间步长的值，预估 tn 1时刻的独立坐标值iq f , P表示预估值。（3）校正用Adams-Moulton隐式公式对上面的预估值，根据给定的收敛误差限进行校正，以得 到独立坐标的校正值 T :C ,C表示校正值。（4） 确定相关坐标确定独立坐标的校正值之后，可由相应公式计算出非独立坐标和其他系统状态变量值。（5） 积分误差控制与上面预估-校正算法积分误差控制过程相同，如果预估

22、值与校正值的差值小于给定 的积分误差限，接受该解，进行下一时刻的求解。否则减小积分步长，重复第二步开始的 预估步骤。4.5.3动力学求解算法特性比较1微分代数（DAE ）方程的求解三种积分格式比较I3积分格式仅监控位移和其它微分方程的状态变量的误差。 当积分步长变小时Jacobian矩阵不能保持稳定，会出现奇异，积分易发散。积分过程不能监控速度和约束反力。因而 速度、加速度、约束反力计算精度差一些。SI2积分格式中考虑了速度约束方程， 可以控制拉氏乘子的误差、 速度误差，仿真结果更精确，可以给出速度、加速度较为精确解。 Jacobian矩阵在步长很小时仍能保持稳定，Jacobian矩阵小步长不

23、会奇异、病态，增加了校正器在小步长时的稳定性和鲁棒性。校正 阶段不会象I3积分格式那样容易失败。可以精确处理高频问题。但比 I3积分格式慢，驱动约束为速度时，输入必须可微、光滑。非光滑驱动约束运动输入会产生无限加速度，而 导致SI2积分失败。位移驱动约束输入不能是变量的函数，速度、加速度输入可以是变量 的函数，而I3驱动约束输入可以是变量的函数，这给仿真带来不便。SI1积分格式中考虑了速度约束方程， 但并没有引入加速度约束方程，相对应引入了拉氏乘子的导数而使方程降阶，可以控制拉氏乘子的误差、速度误差，仿真结果很精确， Jacobian矩阵在步长很小时仍能保持稳定，增加了校正器在小步长时的稳定性

24、和鲁棒性。 可以给出速度、加速度较为精确解，可以监控所有状态变量如位移、速度、拉氏乘子，比 SI2精度高，但对具有摩擦、接触的模型很敏感。三种积分方式比较如下表：表4-1三种积分方式的比较In dex 3SI2SI1求解精度位移精度咼位移，速度，加速度 精度咼位移，速度，加速度， 拉氏乘子精度咼求解稳定性一般好求解速度快般处理高频问题中低频问题适合高频适合2.求解器的特点比较（1） GstiffGstiff求解器为刚性稳定算法，采用多步、变阶（最高阶为 6）、变步长、固定系数算法。可直接求解 DAE方程，有13、SI2、SI1三种积分格式。在预估中采用泰勒级数，而 且其系数是假设步长不变而得到

25、的固定系数，因而当步长改变时会产生误差。其奇特点是 计算速度快，位移精度高，I3格式时速度、尤其加速度会产生误差，可以通过控制最大步 长来控制求解中步长的变化，从而提高精度使仿真运行在定步长状态。 当步长小时，Jacobian矩阵是步长倒数的函数会变成病态， SI2及SI1积分格式时Jacobian矩阵可以步长很小时仍能保持稳定。该算法可以适应很多仿真分析问题。（2） WstiffWstiff求解器为刚性稳定算法，采用多步、变阶（最高阶为 6）、变步长、变系数算法。可直接求解 DAE方程，有I3、SI2、SI1三种积分格式。在预估中采用 NDF （Newton Divided Differe nee）公式用于预估，可以根据步长信息修改相应阶的系数，而且步长改变并不影响 精度，因而更具健壮