ADAMS软件算法基本理论陈立平机械系统动力学分概要.docx
《ADAMS软件算法基本理论陈立平机械系统动力学分概要.docx》由会员分享,可在线阅读,更多相关《ADAMS软件算法基本理论陈立平机械系统动力学分概要.docx(19页珍藏版)》请在冰豆网上搜索。
ADAMS软件算法基本理论陈立平机械系统动力学分概要
第4章ADAM软件基本算法
本章主要介绍ADAMS软件的基本算法,包括ADAMS建模中的一些基本概念、运动学分析算法、动力学分析算法、静力学分析及线性化分析算法以及ADAMS软件积分器介绍。
通过本章的学习可以对ADAMS软件的基本算法有较深入的了解,为今后合理选择积分器进行仿真分析提供理论基础,为更好地使用ADAMS打下良好的理论基础。
4.1ADAMS建模基础
ADAMS利用带拉格朗日乘子的第一类拉格朗日方程导出――最大数量坐标的微分-代数方程(DAE)。
它选取系统内每个刚体质心在惯性参考系中的三个直角坐标和确定刚体方位的三个欧拉角作为笛卡尔广义坐标,用带乘子的拉格朗日第一类方程处理具有多余坐标的完整约束系统或非完整约束系统,导出以笛卡尔广义坐标为变量的动力学方程。
4.1.1参考标架
在计算系统中构件的速度和加速度时,需要指定参考标架,作为该构件速度和加速度的参考坐标系。
在机械系统的运动分析过程中,有两种类型的参考标架——地面参考标架和构件参考标架。
地面参考标架是一个惯性参考系,它固定在一个“绝对静止”的空间中。
通过地面参考标架建立机械系统的“绝对静止”参考体系,属于地面标架上的任何一点的速度和加速度均为零。
对于大多数问题,可以将地球近似为惯性参考标架,虽然地球是绕着太阳旋转而且地球还有自转。
对于每一个刚性体都有一个与之固定的参考标架,称为构件参考标架,刚性体上的各点相对于该构件参考标架是静止的。
4.1.2坐标系的选择
机械系统的坐标系广泛采用直角坐标系,常用的笛卡尔坐标系就是一个采用右手规则的直角坐标系。
运动学和动力学的所有矢量均可以用沿3个单位坐标矢量的分量来表示。
坐标系可以固定在一个参考标架上,也可以相对于参考框架而运动。
合理地设置坐标系可以简化机械系统的运动分析。
在机械系统运动分析过程中,经常使用3种坐标系:
(1)地面坐标系(GroundCoordinateSystem)。
地面坐标系又称为静坐标系,是固定在地面标架上的坐标系。
ADAMS中,所有构件的位置、方向和速度都用地面坐标系表示。
(2)局部构件参考坐标系(LocalPartReferenceFrame,LPRF)。
这个坐标系固定在构件上并随构件运动。
每个构件都有一个局部构件参考坐标系,可以通过确定局部构件参考坐标系在地面坐标系的位置和方向,来确定一个构件的位置和方向。
在ADAMS中,局部构件参考坐标系缺省与地面坐标系重合。
(3)标架坐标系(MarkerSystem)。
标架坐标系又称为标架,是为了简化建模和分析
在构件上设立的辅助坐标系,有两种类型的标架坐标系:
固定标架和浮动标架。
固定标架固定在构件上,并随构件运动。
可以通过固定标架在局部构件参考坐标系中的位置和方向,确定固定标架坐标系的位置和方向。
固定标架可以用来定义构件的形状、质心位置、作用力和反作用力的作用点、构件之间的连接位置等。
浮动标记相对于构件运动,在机械系统的运动分析过程中,有些力和约束需要使用浮动标架来定位。
动力学方程的求解速度很大程度上取决于广义坐标的选择。
研究刚体在惯性空间中的一般运动时,可以用它的质心标架坐标系确定位置,用质心标架坐标相对地面坐标系的方向余弦矩阵确定方位。
为了解析地描述方位,必须规定一组转动广义坐标表示方向余弦矩阵。
第一种方法是用方向余弦矩阵本身的元素作为转动广义坐标,但是变量太多,同时还要附加六个约束方程;第二种方法是用欧拉角或卡尔登角作为转动坐标,它的算法规范,缺点是在逆问题中存在奇点,在奇点位置附近数值计算容易出现困难;第三种方法是用欧拉参数作为转动广义坐标,它的变量不太多,由方向余弦计算欧拉角时不存在奇点。
ADAMS
软件用刚体Bi的质心笛卡尔坐标和反映刚体方位的欧拉角作为广义坐标,即qi二[x,y,zi,Yf,q=[q1,q2,…,q【r。
由于采用了不独立的广义坐标,系统动力
学方程虽然是最大数量,但却是高度稀疏耦合的微分代数方程,适用于稀疏矩阵的方法高效求解。
4.2ADAMS运动学分析
4.2.1ADAMS运动学方程
利用ADAMS建立机械系统仿真模型时,系统中构件与地面或构件与构件之间存在运动副的联接,这些运动副可以用系统广义坐标表示为代数方程,这里仅考虑完整约束。
设表示运动副的约束方程数为nh,则用系统广义坐标矢量表示的运动学约束方程组为:
「K(q)珂皆(q),站(q),…,吒(q)]T=0(4-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(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-2Gqtq…tt,则加速度方程为:
:
:
J(q,q,q,t)q(q,t)q—(q,q,t)0(4-7)
矩阵Gq,为雅可比矩阵,如果门的维数为m,q维数为n,那么Gq维数为mn矩阵,其定义为(①q)(i,j)=/闵)。
在这里①q为ncFc(nh个运动学约束,nc—nh个驱动约束,nc个广义坐标)的方阵。
4.2.2ADAMS运动学方程的求解算法
在ADAMS仿真软件中,运动学分析研究零自由度系统的位置、速度、加速度和约束反力,因此只需求解系统的约束方程:
G(q,tn)=0(4-8)
运动过程中任一时刻tn位置的确定,可由约束方程的Newton-Raphson迭代法求得:
叫j.q:
:
」(qj,tj=0(4-9)
其中,「q=qj4-qj,表示第j次迭代。
tn时刻速度、加速度可以利用线性代数方程的数值方法求解,ADAMS中提供了两种
线性代数方程求解方法:
CALAHAN方法(由Michigan大学DonaldCalahan教授提出)与HARWELL方法(由HARWELL的IanDuff教授提出),CALAHAN方法不能处理冗余约束问题,HARWELL方法可以处理冗余约束问题,CALAHAN方法速度较快。
q=-叮①t(4-10)
dW-Gq,_C:
-qd)qqW"tt(4-11)
4.3ADAMS动力学分析
4.3.1ADAMS动力学方程
ADAMS中用刚体B的质心笛卡尔坐标和反映刚体方位的欧拉角作为广义坐标,即
q二[x,y,乙门,]T,令RJx,y,z「,-i,q=[RT,T]T。
构件质心参考
坐标系与地面坐标系间的坐标变换矩阵为:
cos'-cos—sincostsin€cos'-sin—sincostcossin'-sin|
Agi=sin即cos©+cos即cos日sin©—sin即sin0+cos即cos日cos©—cos即sin日sin日sin©sin日cos©cos日
(4-12)
定义一个欧拉转轴坐标系,该坐标系的三个单位矢量分别为上面三个欧拉转动的轴,因而三个轴并不相互垂直。
该坐标系到构件质心坐标系的坐标变换矩阵为:
sin)sin0cos
B=sin日cos©0-sin日(4-13)
'cose1o_
构件的角速度可以表达为:
~B\(4-14)
ADAMS中引入变量-'e为角速度在欧拉转轴坐标系分量:
(4-15)
考虑约束方程,ADAMS利用带拉格朗日乘子的拉格朗日第一类方程的能量形式得到如下方程:
(4-16)
dt(中』y「'二
dt;qi1:
q
T为系统广义坐标表达的动能,qj为广义坐标,Qj为在广义坐标qj方向的广义力,
qj方向的约束反力。
(4-17)
最后一项涉及约束方程和拉格朗日乘子表达了在在广义坐标
ADAMS中近一步引入广义动量:
简化表达约束反力为:
这样方程(4-16)可以简化为:
动能可以近一步表达为:
TJRTmR1|TBTJBJ
22
其中M为构件的质量阵,J为构件在质心坐标系下的惯量阵。
将(4-19)分别表达为移动方向与转动方向有:
=Qr_Cr
其中pR=dtdmR=mV,——=0。
Rdtl/MR丿dt喩
(4-21)式可以简化为:
mV=Qr-cr
待'T
P=—=BTJB,,由于B中包含欧拉角,为了简化推导,ADAMS
回Y丿
步推导P,而是将其作一个变量求解。
(4-18)
(4-19)
(4-20)
(4-21)
(4-22)
(4-23)
中并没有进
这样ADAMS中每个构件具有如下15个变量(而非12个)和15个方程(而非12个)。
变量:
厂一J
V=M,Vy,Vz]
(4-24)
R=lx,y,zf宀尸冋P&P点COe=[灼屮⑷召灼
Y=即月冲]T
方程:
其中,P为系统的广义动量;H为外力的坐标转换矩阵。
为了更好地说明ADAMS的建模过程下面以一个单摆为例进行建模推导。
图4-1单摆示意图
如图所示,单摆的质量为M、惯量为I,杆长为2L,并在0点以转动副与大地相连接约束在大地的OXY平面内。
在单摆质心处建立单摆的跟随坐标系——局部构件参考坐标系
Op—Xp—Yp,其坐标在地面坐标系OXY中为(x,y),单摆的姿态角为0。
系统的动能表达式:
TMX2My2(4-27)
广义动量表达式:
%MX
外力表达式:
cX
(4-28)
(4-29)
(4-30)
■01HTF=Mg
■0j
约束方程:
x-Leos:
-0
y-Lsin)-0
约束方程的雅克比矩阵:
(4-31)
Lsinv
—LcosB
约束对应的拉格朗日乘子:
(4-32)
力、力矩平衡方程:
一mvL十打=o1
p)—+①qT?
“+HTF=0nMVy+入2-Mg=0(4-33)
L祐十人Lsin日一九2Lcos日=0
动量矩表达式:
P-汀妒P=I.
运动学关系方程:
「Vx-殳=01
U=d=Vy_y=0
[Qg-d=0
其方程集成表达为:
一mVx+入=0I
MVy+打_Mg=0
P+人Lsin日一毎Lcos日=0
IP^T⑷日=。
Vx-X=0
IVy-y=0
co0—tl=0
x-Lcos日=0
y—Lsin日=0
其中系统需求解变量为:
xy二VxVyrPj1'2
(4-34)
(4-35)
(4-36)
(4-37)
432初始条件分析
在进行动力学、静力学分析之前,ADAMS会自动进行初始条件分析,以便在初始系
统模型中各物体的坐标与各种运动学约束之间达成协调,这样可以保证系统满足所有的约束条件。
初始条件分析通过求解相应的位置、速度、加速度的目标函数的最小值得到。
(1)对初始位置分析,需满足约束最小化问题
1T
Minimize:
Cq-q0W(q-q0)
Subjectto:
:
jq=°
q为构件广义坐标,W为权重矩阵,qo为用户输入的值,如果用户输入的值为精确值,则相应权重较大,并在迭代中变化较小。
可以利用拉格朗日乘子将上述约束最小化问题变为如下极值问题:
1.T”T
LSq-q°570+①9'(4-38)
闩FI
L取最小值,则由0,.0得:
cqcK
其中,q°为用户设定的准确的或近似的初始速度值,或者为程序设定的缺省速度值;
w为对应q0的权重系数矩阵。
同样可以利用拉格朗日乘子将上述约束最小化问题变为如下极值问题:
q为已知,该方程为线性方程组可求解如下方程:
(4-43)
(3)对初始加速度、初始拉氏乘子的分析,可直接由系统动力学方程和系统约束方程的两阶导数确定。
4.3.3ADAMS动力学方程的求解
对于式(4-26)微分—代数方程的求解,ADAMS采用两种方式求解,第一种为对DAE
方程的直接求解,第二种为DAE方程利用约束方程将广义坐标分解为独立坐标和非独立坐标然后化简为ODE方程求解。
关于具体求解器将在4.5节介绍。
DAE方程的直接求解将二阶微分方程降阶为一阶微分方程来求解,通过引入u=q,将所有拉格朗日方程均写成一
阶微分形式,该方程为Index3微分代数方程。
I3积分格式:
运用一阶向后差分公式,上述方程组对(uq)求导,可得其Jacobian矩阵,
然后利用Newton-Rapson求解。
可以看出,当积分步长h减小并趋近于0时,上述
Jacobian矩阵呈现病态。
为了有效地监测速度积分的误差,可采用降阶积分方法(Index
reductionmethods)。
通常来说,微分方程的阶数越少,其数值求解稳定性就越好。
ADAMS还采用两种方法来降阶求解,即SI2(Stabilized-lndexTwo)和SI1(Stabilized-lndexOne)方法。
SI2积分格式:
(4-45)
q,t=0
iq,u,t=0
F二f(u,q,t)
上式能同时满足①和去求解不违约,且当步长h趋近于0时,Jacobian矩阵不会呈现
病态现象。
SI1积分格式:
上式中,为了对方程组降阶,引入.和i来替代拉格朗日乘子,即「二,,i=・I。
这
种变化有效地将上述方程组的阶数降为1。
因为只需要微分速度约束方程一次来显示地计
算表达式;和i。
运用SI1积分器,能够方便地监测q,u,n和Z的积分误差,系统的加速度也趋向于更加精确。
但在处理有明显的摩擦接触问题时,SI1积分器十分敏感并具有
挑剔性。
4.4ADAMS静力学及线性化分析
441静力学分析
在进行静力学、准静力学分析时,对动力学方程的速度、加速度设置为零,则得到静
(4-47)
力学方程如下:
u=0
Gq,t=0
F=f(u,q,t)
4.4.2线性化分析
在系统的某点处,q=q”,u=u”可对系统的动力学方程进行线性化,
(4-48)
MUCuKq=0u=q
M,C,K为常数阵
可对(4.4-1)式求解得到系统的频率和振动模态。
4.5ADAMS求解器算法介绍
4.5.1ADAMS数值算法简介
运动学、静力学分析需求解一系列的非线性代数方程、线性代数方程,ADAMS采用
了修正的Newton-Raphson迭代算法求解非线性代数方程,以及基于LU分解的CALAHAN方法和HARWELL方法求解线性代数方程。
对动力学微分方程,根据机械系统特性,选择
不同的积分算法;对刚性系统,采用变系数的BDF(BackwardsDifferentiationFormulation)
刚性积分程序,它是自动变阶、变步长的预估校正法(PECE,Predict-Evaluate-Correct-Evaluate),并分别为Index3、SI2、SI1积分格式,在积分的每一步采用了修正的Newton-Raphson迭代算法;对高频系统(High-Frequencies),采用坐标分块法(Coordinate-PartitionedEquation)将微分—代数(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求解器。
4.5.2动力学求解算法介绍
1•微分—代数(DAE)方程的求解算法过程
ADAMS中DAE方程的求解采用了BDF刚性积分法,以下为其步骤:
(1)预估阶段
用Gear预估-校正算法可以有效地求解微分-代数方程。
首先,根据当前时刻的系统状
态矢量值,用泰勒级数预估下一时刻系统的状态矢量值:
2
yn1=ynnh「》『川(4-49)
Ct2!
Ct
其中,时间步长h二tn1-tn。
这种预估算法得到的新时刻的系统状态矢量值通常不准
确,可以由Geark1阶积分求解程序(或其他向后差分积分程序)来校正。
k
yn1二-hTn「7ryn41(4-50)
i:
±
其中,yn1为y(t)在t=tn1时的近似值;r和J为Gear积分程序的系数值。
_1k
上式经过整理,可表示为:
yn1「[yn1二%—](4-51)
i二
(2)校正阶段
•求解系统方程G,如G(y,y)t)&,则方程成立,此时的y为方程的解,否则继续;
•求解Newton-Raphson线性方程,得到二y,以更新y,使系统方程G更接近于成立。
Jy二G(y,y,tnJ,其中J为系统的雅可比矩阵。
k-J1kk
•利用Newton-Raphson迭代,更新y:
y=y=y
•重复以上步骤直到Ay足够小。
(3)误差控制阶段
•预估计积分误差并与误差精度比较,如积分误差过大则舍弃此步。
•计算优化的步长h和阶数n。
•如达到仿真结束时间,则停止,否则t:
t,重新进入第一步。
2•坐标缩减的微分方程求解过程算法
ADAMS程序提供ABAM(Adams—BashforthandAdams-Moulton)和RKF45积分程序,采用坐标分离算法,将微分-代数方程减缩成用独立广义坐标表示的纯微分方程,然后用
ABAM或RKF45程序进行数值积分。
以下以ABAM为例介绍其求解过程。
坐标减缩微分方程的确定及其数值积分过程按以下步骤进行:
(1)坐标分离
将系统的约束方程进行矩阵的满秩分解,可将系统的广义坐标列阵分解成独立坐
标列阵tq"』和非独立坐标列阵gd』,即gj=g"qdj。
(2)预估
用Adams-Bashforth显式公式,根据独立坐标前几个时间步长的值,预估tn1时刻的独
立坐标值iq"f,P表示预估值。
(3)校正
用Adams-Moulton隐式公式对上面的预估值,根据给定的收敛误差限进行校正,以得到独立坐标的校正值T:
C,C表示校正值。
(4)确定相关坐标
确定独立坐标的校正值之后,可由相应公式计算出非独立坐标和其他系统状态变量值。
(5)积分误差控制
与上面预估-校正算法积分误差控制过程相同,如果预估值与校正值的差值小于给定的积分误差限,接受该解,进行下一时刻的求解。
否则减小积分步长,重复第二步开始的预估步骤。
4.5.3动力学求解算法特性比较
1微分—代数(DAE)方程的求解三种积分格式比较
I3积分格式仅监控位移和其它微分方程的状态变量的误差。
当积分步长变小时Jacobian
矩阵不能保持稳定,会出现奇异,积分易发散。
积分过程不能监控速度和约束反力。
因而速度、加速度、约束反力计算精度差一些。
SI2积分格式中考虑了速度约束方程,可以控制拉氏乘子的误差、速度误差,仿真结果
更精确,可以给出速度、加速度较为精确解。
Jacobian矩阵在步长很小时仍能保持稳定,
Jacobian矩阵小步长不会奇异、病态,增加了校正器在小步长时的稳定性和鲁棒性。
校正阶段不会象I3积分格式那样容易失败。
可以精确处理高频问题。
但比I3积分格式慢,驱
动约束为速度时,输入必须可微、光滑。
非光滑驱动约束运动输入会产生无限加速度,而导致SI2积分失败。
位移驱动约束输入不能是变量的函数,速度、加速度输入可以是变量的函数,而I3驱动约束输入可以是变量的函数,这给仿真带来不便。
SI1积分格式中考虑了速度约束方程,但并没有引入加速度约束方程,相对应引入了拉
氏乘子的导数而使方程降阶,可以控制拉氏乘子的误差、速度误差,仿真结果很精确,Jacobian矩阵在步长很小时仍能保持稳定,增加了校正器在小步长时的稳定性和鲁棒性。
可以给出速度、加速度较为精确解,可以监控所有状态变量如位移、速度、拉氏乘子,比SI2精度高,但对具有摩擦、接触的模型很敏感。
三种积分方式比较如下表:
表4-1三种积分方式的比较
Index3
SI2
SI1
求解精度
位移精度咼
位移,速度,加速度精度咼
位移,速度,加速度,拉氏乘子精度咼
求解稳定性
一般
好
好
求解速度
快
般
般
处理高频问题
:
中低频问题适合
高频适合
高频适合
2.求解器的特点比较
(1)Gstiff
Gstiff求解器为刚性稳定算法,采用多步、变阶(最高阶为6)、变步长、固定系数算
法。
可直接求解DAE方程,有13、SI2、SI1三种积分格式。
在预估中采用泰勒级数,而且其系数是假设步长不变而得到的固定系数,因而当步长改变时会产生误差。
其奇特点是计算速度快,位移精度高,I3格式时速度、尤其加速度会产生误差,可以通过控制最大步长来控制求解中步长的变化,从而提高精度使仿真运行在定步长状态。
当步长小时,Jacobian
矩阵是步长倒数的函数会变成病态,SI2及SI1积分格式时Jacobian矩阵可以步长很小时仍
能保持稳定。
该算法可以适应很多仿真分析问题。
(2)Wstiff
Wstiff求解器为刚性稳定算法,采用多步、变阶(最高阶为6)、变步长、变系数算法。
可直接求解DAE方程,有I3、SI2、SI1三种积分格式。
在预估中采用NDF(NewtonDividedDifferenee)公式用于预估,可以根据步长信息修改相应阶的系数,而且步长改变并不影响精度,因而更具健壮