动力学建模方法与解法总结
目 录
1 刚体系统 ....................................................................................... 1 2 弹性系统动力学 ............................................................................ 6 3 高速旋转体动力学 ...................................................................... 10
1 刚体系统
一般力学研究的对象，是由两个或两个以上刚体通过铰链等约束联系在一起的力学系统，为一般力学研究对象。自行车、万向支架陀螺仪通常可看成多刚体系统。人体在某种意义上也可简化为一个多刚体系统。现代航天器、机器人、人体和仿生学中关于动物运动规律的研究都提出了多刚体系统的一系列理论模型作为研究对象。多刚体系统按其内部联系的拓扑结构，分为树型和非树型（包含有闭链）；按其同外界的联系情况，则有有根和无根之别。利用图论的工具可以一般地分析多刚体系统的构造，建立系统的数学模型和动力学方程组。也可从分析力学中的高斯原理出发，用求极值的优化算法直接求解系统的运动和铰链反力。依照多刚体系统动力学的理论和方法，广泛采用电子计算机对这些模型进行研究，对于精确地掌握这些对象的运动规律是很有价值的。
1.1 自由物体的变分运动方程
任意一个刚体构件i，质量为mi，对质心的极转动惯量为Ji'，设作用于刚体的所有外力向质心简化后得到外力矢量Fi和力矩ni，若定义刚体连体坐标系x'o'y'的原点o'位于刚体质心，则可根据牛顿定理导出该刚体带质心坐标的变分运动方程：
' δriT[mi ri-Fi]+δφi[Jiφi-ni]=0 (1-1)
其中，ri为固定于刚体质心的连体坐标系原点o'的代数矢量，φi为连体坐标系相对于全局坐标系的转角，δri与δφi分别为ri与φi的变分。
定义广义坐标：
qi=[riT,φi]T (1-2)
广义：
Qi=[FiT,ni]T (1-3)
及质量矩阵：
Mi=diag(mi,mi,Ji') (1-4)
体坐标系原点固定于刚体质心时用广义力表示的刚体变分运动方程：
i-Qi)=0 (1-5) δqiT(Miq
1.2 束多体系统的运动方程
考虑由nb个构件组成的机械系统，对每个构件运用式(1-5)，组合后可得到系统的变分运动方程为：
∑δq
i=1
nb
Ti i-Qi]=0 (1-6) [Miq
若组合所有构件的广义坐标矢量、质量矩阵及广义力矢量，构造系统的广义坐标矢量、质量矩阵及广义力矢量为：
TTTT
q=[q1,q2,...,qnb] (1-7)
M=diag(M1,M2,...,Mnb) (1-8)
TTTTQ=[Q1,Q2,...,Qnb] (1-9)
系统的变分运动方程则可紧凑地写为：
-Q]=0 (1-10) δqT[Mq
对于单个构件，运动方程中的广义力同时包含作用力和约束力，但在一个系统中，若只考虑理想运动副约束，根据牛顿第三定律，可知作用在 系统所有构件上的约束力总虚功为零，若将作用于系统的广义外力表示为：
AQA=[Q1A,Q2A,...,Qnb]T (1-11)
T
T
T
其中：
QiA=[FiA,nA]T，i=1,2,...,nb (1-12)
T
则理想约束情况下的系统变分运动方程为：
-QA]=0 (1-13) δqT[Mq
式中虚位移δq与作用在系统上的约束是一致的。
系统运动学约束和驱动约束的组合如式(1-10)，为：
Φ(q,t)=0 (1-14)
对其微分得到其变分形式为：
Φqδq=0 (1-15)
式(1-13)和(1-15)组成受约束的机械系统的变分运动方程。
为导出约束机械系统变分运动方程易于应用的形式，运用拉格朗日乘子定理对式(1-13)和(1-15)进行处理。
拉格朗日乘子定理：设矢量b∈Rn，矢量x∈Rn，矩阵A∈Rm⨯n为常数矩阵，如果有：
bTx=0 (1-16)
对于所有满足式（1-84）的x条件都成立。
Ax=0 (1-17)
则存在满足式（1-85）的拉格朗日乘子矢量λ∈Rm。
bTx+λTAx=0 (1-18)
其中x为任意的。
在式(1-13)和(1-15)中，q∈Rn，M∈Rn⨯n，QA∈Rn，Φq∈R
m⨯n
，运用
拉格朗日乘子定理于式(1-13)和(1-15)，则存在拉格朗日乘子矢量λ∈Rm，对于任意的δq应满足：
AT
-QA]Tδq+λTΦqδq=[Mq +ΦT[Mqqλ-Q]δq=0 (1-19)
由此得到运动方程的拉格朗日乘子形式：
A
+ΦTMqλ=Q (1-20) q
式(1-20)还必须满足式(1-10)、(1-12)和(1-14)表示的位置约束方程、速度约束方程及加速度约束方程，如下：
Φ(q,t)=0 (1-21)
(q,q ,t)=Φq(q,t)q -υ=0，υ=-Φt(q,t) (1-22) Φ
(q,q )qq -2Φqtq -Φtt (1-23) ,q ,t)=Φq(q,t)q -η(q,q ,t)=0，η=-(ΦqqΦ
以上三式其维数同式(1-14)。
(1-21)、(1-22)和(1-23)组成约束机械系统的完整的运动方程。 式(1-20)、
将式(1-20)与(1-23)联立表示为矩阵形式：
⎡M⎢⎢Φq⎣⎤ ΦT⎡QA⎤q⎡q⎤
⎥⎢⎥=⎢⎥ (1-24) 0⎦⎥⎣λ⎦⎣η⎦
式(1-24)即为多体系统动力学中最重要的动力学运动方程，式(1-24)还必须满足式(1-22)和(1-23)。它是一个微分——代数方程组，不同于单纯的常微分方程组问题，其求解关键在于避免积分过程中的违约现象，此外，还要注意DAE问题的刚性问题。
如果系统质量矩阵是正定的，并且约束独立，那么运动方程就有唯一解。实际中的系统质量矩阵通常是正定的，只要保证约束是独立的，运动方程就会有解。
在实际数值迭代求解过程中，需要给定初始条件，包括位置初始条件
(t0)。此时，如果要使运动方程有解，还需要满足初q(t0)和速度初始条件q
值相容条件，也就是要使位置初始条件满足位置约束方程，速度初始条件(1-22)确定的系统动力学方程，满足速度约束方程。对于由式(1-24)及(1-21)、初值相容条件为：
Φ(q(t0),t0)=0 (1-25)
(q(t0),q (t0),t0)=Φq(q(t0),t0)q (t0)-υ(q(t0),t0)=0 (1-26) Φ
1.3 正向动力学分析、逆向动力学分析与静平衡分析
对于一个确定的约束多体系统，其动力学分析不同于运动学分析，并不需要系统约束方程的维数m等于系统广义坐标的维数n，m
如果约束多体系统约束方程的维数m与系统广义坐标的维数n相等，
m=n，也就是对系统施加与系统自由度相等的驱动约束，那么该系统在运
动学上就被完全确定，由2.2.3节的约束方程、速度方程和加速度方程可求解系统运动。在此情况下，雅可比矩阵是非奇异方阵，即：
Φq(q,t)≠0 (1-27)
展开式(1-24)的运动方程，为：
A
+ΦTMqqλ=Q (1-28)
=η (1-29) Φqq
，再由式(1-28)可求得λ，拉格朗日乘子λ就唯一地确定由式(1-29)可解得q
了作用在系统上的约束力和力矩（主要存在于运动副中）。这种由确定的运动求系统约束反力的动力学分析就是逆向动力学分析。
如果一个系统在外力作用下保持静止状态，也就是说，如果：
=q =0 (1-30) q
那么，就说该系统处于平衡状态。将式(1-30)代入运动方程式(1-20)，得到平衡方程：
AΦTqλ=Q (1-31)
由平衡方程式(1-21)及约束方程式(1-13)可求出状态q和拉格朗日乘子
λ。这种求系统的平衡状态及在平衡状态下的约束反力的动力学分析称为（静）平衡分析。
1.4 约束反力
对于约束机械系统中的构件i，设其与系统中某构件j存在运动学约束或驱动约束，约束编号为k。除连体坐标系x'o'y'外，再在构件i上以某点P为原点建立一个新的固定于构件上的坐标系x''Py''，称为运动副坐标系，设从坐标系x''Py''到坐标系x'o'y'的变换矩阵为Ci，从坐标系x'o'y'到坐标系
xoy的变换矩阵为Ai，则可导出由约束k产生的反作用力和力矩分别为：
k
Fi''k=-CiTAiTΦkriλ (1-32) kk
(1-33) Ti''k=(si'PBiTΦk-Φ)λriφi
T
T
T
T
以上两式中，λk为约束k对应的拉格朗日乘子，反作用力Fi''和力矩Ti''均为运动副坐标系x''Py''中的量。
kk
2 弹性系统动力学
由于工业机器人、机械手、弹性联动装置、带柔性附件人造卫星、直升飞机的旋翼等工程结构发展的需求， 使运动中的弹性结构的动力学分析得到了很大的进展。运动弹性体的动力学分析属于多体系统动力学的范畴。而导出其有限元格式的动力学方程并研究其数值解法则是计算多体系统动力学的任务。由于弹性变形与刚体运动的耦合导致了运动弹性体的动力学方程为时变的或非线性的，因此运动中的弹性体会出现诸多非线性效应。
运动中弹性体的动力分析问题可分为两类, 其一是具有给定刚体运动的弹性体的动力分析，这类问题仅讨论弹性体的刚体运动对其弹性变形的影响，比如机械手的弹性终端杆的振动分析一般可归于此类。第二类问题是多体系统中之刚体运动与其中的弹性体的弹性变形的相互耦合的动力分析, 在这类问题中, 弹性体的变形会受到系统刚体运动的影响， 反之弹性体的变形也会影响系统的刚体运动。
下面采用运动参考系方法并用Jourdain 动力学普遍方程导出了具有空间一般运动的弹性体之通用的有限元动力学方程，其最大的优点在于推导简单并适用于各类结构及各种单元形式。对系统的动力学方程的数值求解, 一般可以采用直接积分法。下面给出了对时变的运动弹性的动力学方程的Neumann 级数2直接积分解法, 该方法可以在保证计算精度的前提下很大程度地节省机时。
图2-1
图2-1 所示为一运动的弹性体B，选用两个坐标系来定义弹性体B的刚体运
123动与弹性变形:静系—oxoxoxo， 简记o系; 原点在B上的o1点, 固连于B
上的动
123
系—o1xo简记为o1系。B的刚体移动由o1点对于o点的矢量ro1，定义B的xx，1o1o1
空间转动则用o1系对o系的转动来定义, 而B内任意点P的弹性变形则用在o1系内的弹性变形位移矢量u来表示。
由图可见B发生弹性变形后, 其上任意一点P对o系的位置矢量可以表示为:
rp=ro1+ru
（2-1）
而
ru=r+u （2-2）
其中是B未产生弹性变形时P点在o1系中的位置矢量，则表示P点的弹性变形位移矢量。把(2-2) 式代入(2-1) 式并向o系投影, 且采用矩阵形式表示为:
{r}={r}+[A]({r}+{u}) （2-3）
o
p
oo1
oo1
o1
o1
ooo1
rr其中rp和roo分别表示和向 系的投影列阵；表示o1系向o系转移Aoop11
{}{}[]
的方向余弦矩阵。把(3-3) 式中uo1的用有限元的格式，表达为:
{}
{u}=[N]{ρ} （2-4）
o1
其中[N]为单元形函数矩阵，{ρ}为P点所在单元的有限元结点位移列阵。把(2-4) 式代入(2-3)式, 并利用公式:
（2-5）⎡∙oo1⎤oo1oo1
A=ΩA ⎢⎥⎣⎦
[][]
其中Ωoo1 是o系相对于o系转动角速度在o系上投影的斜对称阵。
1
由(2-3) 式对时间分别求一次导数和二次导数可得P点的速度vop和加速度
[]
{}
{a}，进而可得到P点的虚速度δ{v}，于是P点邻域之微元体的Jourdain 动力
op
op
学普遍方程可以写作:
oo
δ{vop}([A
T
1
]{f}-mdv{a})=0 （2-6）
p
o
p
{f}是作用于P点微元体上的全部力在O1其中: mp为弹性体在P点的质量密度；
系上的投影。
oo1
{f}可利用常规有限元的格式将它写作: 对于δ{vop}A
T
[]
δ{v
oT
p
}[A]
oo1
⎧∙⎫⎛T⎧∙⎫⎫{f}=δ⎨ρ⎬ [N]{F}-[K]{ρ}-[C]⎨ρ⎬⎪ （2-7）⎩⎭⎝⎩⎭⎭
T
其中: [K]和[C]分别为单元刚度阵和单元阻力阵在P点的值; {F}为作用在P点微元体上的外力在O1系的列阵, 把求得的P点的虚速度和加速度以及(2-7) 式代
∙⎧入(2-6) 式, 并考虑到δ⎨ρ⎫ ⎬中诸元素之独立性, 可得P点微元体的动力学方程为：⎩⎭
⎧∙⎫
[N]{F}-[K]{ρ}-[C]⎨ρ⎬-mpdvδVpo
⎩⎭
T
{}{a}=0
T
op
（2-8）
将(2-8) 式对单元积分便可得运动的弹性体的单元动力学方程:
[]
e
∙
⎧∙∙⎫e⎧⎫ M⎨ρ⎬+C⎨ρ⎬+Ke{ρ}=Fe
（2-9）⎩⎭⎩⎭
[][]{}
式中:
[M]=⎰m[N][N]dv
[C]
e
T
p
e
[K]
e
=⎰[N]{μ}[N]dv+2⎰mp[N]A
T
T
[][Ω][A][N]dv=[C]+[C]
oo1T
oo1
oo1
s
d
oo1T
=⎰[B][D][B]dv+⎰mp[N]A
T
T
[]
=[Ks]+[Kd]
⎛⎡∙oo1⎤
⎢Ω⎥+Ωoo1Ωoo1
⎦⎝⎣
⎪[A][N]dv[][]⎫⎪
oo1
⎭
{F}
e
∙∙⎛T⎧o⎫TToo1
=⎰[N]{F}dv- ⎰mp[N]Aro1⎬dv+⎰mp[N]Aoo1⎨ ⎩⎭⎝={Fs}+{Fd}
T
[][]
T
⎛⎡∙oo1⎤
⎢Ω⎥+Ωoo1Ωoo1
⎦⎝⎣
[][][]{}
⎫oo1o1⎫
⎪Ardv⎪⎪⎪⎭⎭
其中[Cs]，[Ks]，{Fs}分别是常规有限元法中的单元阻力阵、刚度阵和外力向量, 而[Cd]，[Kd]，{Fd}则分别是由于刚体运动与弹性变形的耦合而产生的附加单元动力阻尼阵、动力刚度阵和动力力向量。而且由于它们的表达式中含有表示弹性
⎧∙∙⎫
体空间运动量⎨roo⎬和{ω}， 因此,通常这些动力附加项是时变的。当弹性体的刚1
⎩⎭体运动速度特别是转动速度较大时, 弹性体受到较大的惯性力作用, 会产生变形的耦合效应。例如转动的梁, 由于离心惯性力产生的轴向拉力会增大梁的抗弯刚度, 即所谓的“刚化效应”。这时在(2-10) 式中的常规刚度阵[Ks] 中需计入结构的几何刚度阵, 关于各类单元的几何刚度阵可参阅有关非线性有限元的书籍。而结构的几何刚度阵往往是未知内力的函数, 这时方程(2-9)式就是一个非线性的动力方程。但对于简单的弹性体, 如梁, 由于刚体运动的惯性力产生的轴力容易求得, 这时的几何刚度阵就变为时变阵。本文只讨论几何刚度阵为时变阵的情况, 即方程(2-9)式为时变动力学方程时的数值解法。
显然, 若弹性体没有刚体运动, 则方程(2-9)式退化为常规的有限单元动力学方程。把(2-9)式按常规有限元的组集方法进行组集, 便可得到对于运动弹性体的具有时变特性的、通用的有限元动力学方程:
⎧∙∙⎫⎧∙⎫
[M]⎨ρ⎬+[C]⎨ρ⎬+[K]{ρ}={F}（2-10） ⎩⎭⎩⎭
3 高速旋转体动力学
高速旋转体通常是由是由三个刚体──外环、内环、转子互相约束在一起而成，可使陀螺仪转子具有空间转动的三个自由度。过去曾长期认为，高速自转的平衡对称卡登陀螺仪和单刚体陀螺仪的理论模型没有本质区别，具有所谓“定轴性。但实际上,理论研究和精密的实验研究都已证明这个想法是错误的。平衡对称卡登陀螺仪的空间定向大都具有里雅普诺夫意义下的不稳定性（见运动稳定性）。卡登陀螺仪和单刚体陀螺仪模型有本质区别，只有通过多刚体系统模型的研究才能正确解释卡登陀螺仪的动力学特征。
图3-1
如图3-1所示，对于外径D与长度l的比值D/l
'rb'和mb'、mb''的重径积mb''rb''。但是这种方法所带来问题是力多边形不正质量mb
易求解以及图解法不够精确。假如采用平面解法，不仅简单正确，而且对于没有动平衡机的工厂无疑有一定的实用价值。
上述转子质量分布简图如图3-2所示，不平衡质量m1、m2、m3分别分布在与回转轴线垂直的三个平面1、2、3内, 各质点距回转轴线的矢径分别为r1、r2、
r3。当转子以等角速度。回转时, 各质点所产生的离心惯性力分别为
2
P1=m1r1ω （3-1）
P2=m2r1ω2 （3-2）
P3=m3r1ω2 （3-3）
图3-2
T''（过方向如图所示。若选择转子左、右二端面T'（过点A与轴线垂直的平面）、
'、T''平面内分别加上校正质量mb点B与轴线垂直的平面）作为校正平面， 在T'、''，矢径为rb'、rb''，则校正质量所产生的离心惯性力为Pb'=mb'rb'ω2和mb
''rb''Pb''=mbω2，P1、P2、P3、Pb'和Pb''组成了空间力系。
选取三坐标轴x、y、z轴如图所示，并将作用在转子上的所有力向YAZ平面和XAY平面投影，如图3-3所示。
图3-3
在图3-3中, 所有的力组成了平面平行力系, 列平衡方程：
⎛→⎫'2zl2'+Pbz''l=0 （3-4） F⎪=0，P∑mA 1zl1-P⎝⎭
∑F
解得：
z
'+P''=0，Pbz1z-P2z+Pbz=0 （3-5）
''=Pbz
''-P1zl1'P2zl2
（3-6） l
Pbz'=-P2''z-P1z+P2z 式中：
P1z-P1在Z轴上的投影，P1z=P1cosα1，N； P2z-P2在Z轴上的投影，Pz=P2cosα2，N； Pbz
'-Pb'在Z轴上的投影，Pbz'=Pb'cosβ'，N； Pbz
''-Pb''在Z轴上的投影，Pbz''=Pb''cosβ''，N； 同理，在XAY平面内
∑m⎛→⎫A ⎝F⎪⎭
=0，P1xl1'-P2xl2'-P3xl3'+Pbx''l=0 ∑F
z
=0，Pbx-P1x-P2x+P3x+Pbx''=0 解得：
P''=P1xl1'+P2xl'2
-P3xl3'bx
l
Pbx'=P1x+P2x+P3x-Pbx'' 式中：
P1x-P1在X轴上的投影，P1x=P1sinα1，N； P2x-P2在X轴上的投影，P2x=P2sinα2，N； P3x-P3在X轴上的投影，P3x=P3，N； Pbx
'-Pb'在X轴上的投影，Pbx=Pb'sinβ'，N； （3-7）
（3-8） （3-9）
3-10）
3-11）
（ （
''=Pb''sinβ''，N； ''-Pb''在X轴上的投影，PbxPbx
'所产生的离心惯性力Pb'： 在校正平面T'上，校正质量mb
'2+Pbz'2 （3-12） Pb'=Pbx
β'=tan-1
'Pbx
（3-13） 'Pbz
''所产生的离心惯性力Pb''： 在校正平面T''上，校正质量mb
''2+Pbz''2 （3-14） Pb''=Pbx
β''=tan-1
''Pbx
（3-15） ''Pbz
''''''
当给定校正质量mb和mb的失径rb'和rb''，很容易求得校正质量mb和mb的数值：
'=mb
Pb'
rb'ω2 （3-16） Pb''rb''ω2 （3-17）
''=mb
2'mb''mb
和的方位由角度β'和β''确定，如图3-3所示。因为方程每一项均有ω，
故ω可以约去而不必计算。
2