(完整版)动力学建模方法与解法总结

目录
1 刚体系统 (1)
2 弹性系统动力学 (6)
3 高速旋转体动力学 (10)
1 刚体系统
一般力学研究的对象，是由两个或两个以上刚体通过铰链等约束联系在一起的力学系统，为一般力学研究对象。

自行车、万向支架陀螺仪通常可看成多刚体系统。

人体在某种意义上也可简化为一个多刚体系统。

现代航天器、机器人、人体和仿生学中关于动物运动规律的研究都提出了多刚体系统的一系列理论模型作为研究对象。

多刚体系统按其内部联系的拓扑结构，分为树型和非树型（包含有闭链）；按其同外界的联系情况，则有有根和无根之别。

利用图论的工具可以一般地分析多刚体系统的构造，建立系统的数学模型和动力学方程组。

也可从分析力学中的高斯原理出发，用求极值的优化算法直接求解系统的运动和铰链反力。

依照多刚体系统动力学的理论和方法，广泛采用电子计算机对这些模型进行研究，对于精确地掌握这些对象的运动规律是很有价值的。

1.1 自由物体的变分运动方程
任意一个刚体构件i ，质量为i m ，对质心的极转动惯量为i J '，设作用于刚体的所有外力向质心简化后得到外力矢量i F 和力矩i n ，若定义刚体连体坐标系y o x '''的原点o '位于刚体质心，则可根据牛顿定理导出该刚体带质心坐标的变分运动方程：
0][][=-'+-i
i i i i i i T i n J F r m r φδφδ&&&& (1-1) 其中，i r 为固定于刚体质心的连体坐标系原点o '的代数矢量，i φ为连体坐标系相对于全局坐标系的转角，i r δ与i δφ分别为i r 与i φ的变分。

定义广义坐标：
T i T i i r q ],[φ= (1-2)
广义：
T i T i i n F Q ],[= (1-3)
及质量矩阵：
),,(i i i i J m m diag M '= (1-4)
体坐标系原点固定于刚体质心时用广义力表示的刚体变分运动方程：
0)(=-i i i T i Q q M q &&δ (1-5)
1.2 束多体系统的运动方程
考虑由nb 个构件组成的机械系统，对每个构件运用式(1-5)，组合后可得到系统的变分运动方程为：
0][1=-∑=i i i nb i T i Q q M q
&&δ (1-6)
若组合所有构件的广义坐标矢量、质量矩阵及广义力矢量，构造系统的广义坐标矢量、质量矩阵及广义力矢量为：
T T nb T T q q q q ],...,,[21= (1-7)
),...,,(21nb M M M diag M = (1-8)
T T nb T T Q Q Q Q ],...,,[21= (1-9)
系统的变分运动方程则可紧凑地写为：
0][=-Q q M q T &&δ (1-10)
对于单个构件，运动方程中的广义力同时包含作用力和约束力，但在一个系统中，若只考虑理想运动副约束，根据牛顿第三定律，可知作用在
系统所有构件上的约束力总虚功为零，若将作用于系统的广义外力表示为：
T T
A nb T A T A A Q Q Q Q ],...,,[21= (1-11) 其中：
T A T
A i A i n F Q ],[=，nb i ,...,2,1= (1-12) 则理想约束情况下的系统变分运动方程为：
0][=-A T Q q M q &&δ (1-13)
式中虚位移q δ与作用在系统上的约束是一致的。

系统运动学约束和驱动约束的组合如式(1-10)，为：
0),(=Φt q (1-14)
对其微分得到其变分形式为：
0=Φq q δ (1-15)
式(1-13)和(1-15)组成受约束的机械系统的变分运动方程。

为导出约束机械系统变分运动方程易于应用的形式，运用拉格朗日乘子定理对式(1-13)和(1-15)进行处理。

拉格朗日乘子定理：设矢量n R b ∈，矢量n R x ∈，矩阵n m R A ⨯∈为常数矩阵，如果有：
0=x b T (1-16)
对于所有满足式（1-84）的x 条件都成立。

0=Ax (1-17)
则存在满足式（1-85）的拉格朗日乘子矢量m R ∈λ。

0=+Ax x b T T λ (1-18)
其中x 为任意的。

在式(1-13)和(1-15)中，n R q ∈，n n R M ⨯∈，n A R Q ∈，n m q R ⨯∈Φ，运用拉格朗日乘子定理于式(1-13)和(1-15)，则存在拉格朗日乘子矢量m R ∈λ，对于任意的q δ应满足：
0][][=-Φ+=Φ+-q Q q M q q Q q M T A T
q q T T A δλδλδ&&&& (1-19)
由此得到运动方程的拉格朗日乘子形式：
A T q
Q q M =Φ+λ&& (1-20) 式(1-20)还必须满足式(1-10)、(1-12)和(1-14)表示的位置约束方程、速度约束方程及加速度约束方程，如下：
0),(=Φt q (1-21)
0),(),,(=-Φ=Φυq t q t q q q &&&，),(t q t Φ-=υ (1-22)
0),,(),(),,,(=-Φ=Φt q q q t q t q q q q &&&&&&&&η，tt qt q q q q q
Φ-Φ-Φ-=&&&2)(η (1-23) 以上三式其维数同式(1-14)。

式(1-20)、
(1-21)、(1-22)和(1-23)组成约束机械系统的完整的运动方程。

将式(1-20)与(1-23)联立表示为矩阵形式：
⎥⎦
⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎥⎦⎤⎢⎢⎣⎡ΦΦηλA q T q Q q M &&0 (1-24) 式(1-24)即为多体系统动力学中最重要的动力学运动方程，式(1-24)还必须满足式(1-22)和(1-23)。

它是一个微分——代数方程组，不同于单纯的常微分方程组问题，其求解关键在于避免积分过程中的违约现象，此外，还要注意DAE 问题的刚性问题。

如果系统质量矩阵是正定的，并且约束独立，那么运动方程就有唯一解。

实际中的系统质量矩阵通常是正定的，只要保证约束是独立的，运动方程就会有解。

在实际数值迭代求解过程中，需要给定初始条件，包括位置初始条件
)(0t q 和速度初始条件)(0t q &。

此时，如果要使运动方程有解，还需要满足初
值相容条件，也就是要使位置初始条件满足位置约束方程，速度初始条件
满足速度约束方程。

对于由式(1-24)及(1-21)、
(1-22)确定的系统动力学方程，初值相容条件为：
0)),((00=Φt t q (1-25)
0)),(()()),(()),(),((00000000=-Φ=Φt t q t q t t q t t q t q q υ&&& (1-26)
1.3 正向动力学分析、逆向动力学分析与静平衡分析
对于一个确定的约束多体系统，其动力学分析不同于运动学分析，并不需要系统约束方程的维数m 等于系统广义坐标的维数n ，n m <。

在给定外力的作用下，从初始的位置和速度，求解满足位置约束式(1-22)及速度约束式(1-23)的运动方程式(1-24)，就可得到系统的加速度和相应的速度、位置响应，以及代表约束反力的拉格朗日乘子，这种已知外力求运动及约束反力的动力学分析，称为正向动力学分析。

如果约束多体系统约束方程的维数m 与系统广义坐标的维数n 相等，n m =，也就是对系统施加与系统自由度相等的驱动约束，那么该系统在运动学上就被完全确定，由2.2.3节的约束方程、速度方程和加速度方程可求解系统运动。

在此情况下，雅可比矩阵是非奇异方阵，即：
0),(≠Φt q q (1-27)
展开式(1-24)的运动方程，为：
A T q
Q q M =Φ+λ&& (1-28) η=Φq q && (1-29)
由式(1-29)可解得q &&，再由式(1-28)可求得λ，拉格朗日乘子λ就唯一地确定
了作用在系统上的约束力和力矩（主要存在于运动副中）。

这种由确定的运动求系统约束反力的动力学分析就是逆向动力学分析。

如果一个系统在外力作用下保持静止状态，也就是说，如果：
0==q
q &&& (1-30) 那么，就说该系统处于平衡状态。

将式(1-30)代入运动方程式(1-20)，得到平衡方程：
A T
q Q =Φλ (1-31)
由平衡方程式(1-21)及约束方程式(1-13)可求出状态q 和拉格朗日乘子λ。

这种求系统的平衡状态及在平衡状态下的约束反力的动力学分析称为（静）平衡分析。

1.4 约束反力
对于约束机械系统中的构件i ，设其与系统中某构件j 存在运动学约束或驱动约束，约束编号为k 。

除连体坐标系y o x '''外，再在构件i 上以某点P 为原点建立一个新的固定于构件上的坐标系y P x ''''，称为运动副坐标系，设从坐标系y P x ''''到坐标系y o x '''的变换矩阵为i C ，从坐标系y o x '''到坐标系xoy 的变换矩阵为i A ，则可导出由约束k 产生的反作用力和力矩分别为：
k T
k
r T i T i k i i A C F λ
Φ-='' (1-32) k T k T k r T i T P i k i i i B s T λφ)(Φ-Φ'='' (1-33) 以上两式中，k λ为约束k 对应的拉格朗日乘子，反作用力k i F ''和力矩k
i T ''均为运动副坐标系y P x ''''中的量。

2 弹性系统动力学
由于工业机器人、机械手、弹性联动装置、带柔性附件人造卫星、直升飞机的旋翼等工程结构发展的需求， 使运动中的弹性结构的动力学分析得到了很大的进展。

运动弹性体的动力学分析属于多体系统动力学的范畴。

而导出其有限元格式的动力学方程并研究其数值解法则是计算多体系统动力学的任务。

由于弹性变形与刚体运动的耦合导致了运动弹性体的动力学方程为时变的或非线性的，因此运动中的弹性体会出现诸多非线性效应。

运动中弹性体的动力分析问题可分为两类, 其一是具有给定刚体运动的弹性体的动力分析，这类问题仅讨论弹性体的刚体运动对其弹性变形的影响，比如机械手的弹性终端杆的振动分析一般可归于此类。

第二类问题是多体系统中之刚体运动与其中的弹性体的弹性变形的相互耦合的动力分析, 在这类问题中, 弹性体的变形会受到系统刚体运动的影响， 反之弹性体的变形也会影响系统的刚体运动。

下面采用运动参考系方法并用Jourdain 动力学普遍方程导出了具有空间一般运动的弹性体之通用的有限元动力学方程，其最大的优点在于推导简单并适用于各类结构及各种单元形式。

对系统的动力学方程的数值求解, 一般可以采用直接积分法。

下面给出了对时变的运动弹性的动力学方程的Neumann 级数2直接积分解法, 该方法可以在保证计算精度的前提下很大程度地节省机时。

图2-1
图2-1 所示为一运动的弹性体B ，选用两个坐标系来定义弹性体B 的刚体运
动与弹性变形:静系—321o o o
x x ox ， 简记o 系; 原点在B 上的1o 点, 固连于B 上的动
系—32111
11o o o x x x o ，简记为1o 系。

B 的刚体移动由1o 点对于o 点的矢量1o r ，定义B 的空间转动则用1o 系对o 系的转动来定义, 而B 内任意点P 的弹性变形则用在1o 系内的弹性变形位移矢量u 来表示。

由图可见B 发生弹性变形后, 其上任意一点P 对o 系的位置矢量可以表示为:
u o p r r r +=1 （2-1）
而
u r r u += （2-2）
其中是B 未产生弹性变形时P 点在1o 系中的位置矢量，则表示P 点的弹性变形位移矢量。

把(2-2) 式代入(2-1) 式并向o 系投影, 且采用矩阵形式表示为: {}{}[]{}{}()1
111o o oo o o o
p u r A r r ++= （2-3） 其中{}o p r 和{}o o r 1分别表示p r 和1o r 向o 系的投影列阵；[]
1oo A 表示1o 系向o 系转移的方向余弦矩阵。

把(3-3) 式中{}
1o u 的用有限元的格式，表达为: {}[]{}ρN u o =1
（2-4） 其中[]N 为单元形函数矩阵，{}ρ为P 点所在单元的有限元结点位移列阵。

把(2-4) 式代入(2-3)式, 并利用公式:
[][]111oo oo oo A A Ω=⎥⎦⎤⎢⎣⎡• （2-5）
其中[]1oo Ω 是1
o 系相对于o 系转动角速度在o 系上投影的斜对称阵。

由(2-3) 式对时间分别求一次导数和二次导数可得P 点的速度{}
o
p v 和加速度{}o
p a ，进而可得到P 点的虚速度{}o p
v δ，于是P 点邻域之微元体的Jourdain 动力学普遍方程可以写作:
{}[]{}{}()0d 1T =-o
p p oo o
p a v m f A v δ （2-6）
其中: p m 为弹性体在P 点的质量密度；
{}f 是作用于P 点微元体上的全部力在1O
系上的投影。

对于{}[]{}f A v oo o
p 1
T δ可利用常规有限元的格式将它写作: {}[]{}[]{}[]{}[]⎪⎭⎫ ⎝⎛⎭⎬⎫⎩⎨⎧--⎭⎬⎫⎩⎨⎧=••ρρρδδC K F N f A v oo o p T T
T
1 （2-7） 其中: []K 和[]C 分别为单元刚度阵和单元阻力阵在P 点的值; {}F 为作用在P 点微元体上的外力在1O 系的列阵, 把求得的P 点的虚速度和加速度以及(2-7) 式代
入(2-6) 式, 并考虑到⎭
⎬⎫⎩⎨⎧•ρδ中诸元素之独立性, 可得P 点微元体的动力学方程为： []{}[]{}[]{}{}0d T T =-⎭⎬⎫⎩⎨⎧--•o
p o p p a V v m C K F N δρρ （2-8）
将(2-8) 式对单元积分便可得运动的弹性体的单元动力学方程:
[][][]{}{}
e e e e F K C M =+⎭⎬⎫⎩⎨⎧+⎭⎬⎫⎩⎨⎧•••ρρρ （2-9） 式中:
[][][]⎰=v N N m M p
e d T
[][]{}[][][][][][][][]⎰⎰+=Ω+=d
s oo oo oo p e
C C v N A A N m v N N C d 2d 111T T T μ [][][][][][]
[][][][][][]d s oo oo oo oo oo p e K K v N A A N m v B D B K +=⎪⎪⎭
⎫ ⎝⎛ΩΩ+⎥⎦⎤⎢⎣⎡Ω+=⎰⎰•d d 11111T T T {}[]{}[][][][][][][]{}{}{}d s o oo oo oo oo oo p o o oo p e
F F v r A A N m v r A N m v F N F +=⎪⎪⎭
⎫ ⎝
⎛⎪⎪⎭⎫ ⎝⎛ΩΩ+⎥⎦⎤⎢⎣⎡Ω+⎭⎬⎫⎩⎨⎧-=•••⎰⎰⎰d d d 11111111T T T T T 其中[]s C ，[]s K ，{}s F 分别是常规有限元法中的单元阻力阵、刚度阵和外力向量, 而[]d C ，[]d K ，{}d F 则分别是由于刚体运动与弹性变形的耦合而产生的附加单元动力阻尼阵、动力刚度阵和动力力向量。

而且由于它们的表达式中含有表示弹性
体空间运动量⎭
⎬⎫⎩⎨⎧••o o r 1和{}ω， 因此,通常这些动力附加项是时变的。

当弹性体的刚体运动速度特别是转动速度较大时, 弹性体受到较大的惯性力作用, 会产生变形的耦合效应。

例如转动的梁, 由于离心惯性力产生的轴向拉力会增大梁的抗弯刚度, 即所谓的“刚化效应”。

这时在(2-10) 式中的常规刚度阵[]s K 中需计入结构的几何刚度阵, 关于各类单元的几何刚度阵可参阅有关非线性有限元的书籍。

而结构的几何刚度阵往往是未知内力的函数, 这时方程(2-9)式就是一个非线性的动力方程。

但对于简单的弹性体, 如梁, 由于刚体运动的惯性力产生的轴力容易求得, 这时的几何刚度阵就变为时变阵。

本文只讨论几何刚度阵为时变阵的情况, 即方程(2-9)式为时变动力学方程时的数值解法。

显然, 若弹性体没有刚体运动, 则方程(2-9)式退化为常规的有限单元动力学方程。

把(2-9)式按常规有限元的组集方法进行组集, 便可得到对于运动弹性体的具有时变特性的、通用的有限元动力学方程:
[][][]{}{}F K C M =+⎭⎬⎫⎩⎨⎧+⎭⎬⎫⎩⎨⎧•••ρρρ
（2-10）
3 高速旋转体动力学
高速旋转体通常是由是由三个刚体──外环、内环、转子互相约束在一起而成，可使陀螺仪转子具有空间转动的三个自由度。

过去曾长期认为，高速自转的平衡对称卡登陀螺仪和单刚体陀螺仪的理论模型没有本质区别，具有所谓“定轴性。

但实际上,理论研究和精密的实验研究都已证明这个想法是错误的。

平衡对称卡登陀螺仪的空间定向大都具有里雅普诺夫意义下的不稳定性（见运动稳定性）。

卡登陀螺仪和单刚体陀螺仪模型有本质区别，只有通过多刚体系统模型的研究才能正确解释卡登陀螺仪的动力学特征。

图3-1
如图3-1所示，对于外径D 与长度l 的比值5/<l D 的转子，如多缸内燃机的曲轴、机床主轴等，这些转子的不平衡质点不是集中在同一平面内，而是分布在垂直于轴线的各个平面内。

对于这种转子动平衡问题, 一般都采用矢量法来求校
正质量b
m '、b m ''的重径积b b r m ''ρ和b b r m ''''ρ。

但是这种方法所带来问题是力多边形不易求解以及图解法不够精确。

假如采用平面解法，不仅简单正确，而且对于没有动平衡机的工厂无疑有一定的实用价值。

上述转子质量分布简图如图3-2所示，不平衡质量1m 、2m 、3m 分别分布在与回转轴线垂直的三个平面1、2、3内, 各质点距回转轴线的矢径分别为1r 、2r 、
3r 。

当转子以等角速度。

回转时, 各质点所产生的离心惯性力分别为
2111ωr m P = （3-1） 2122ωr m P = （3-2）
2133ωr m P = （3-3）
图3-2
方向如图所示。

若选择转子左、右二端面T '（过点A 与轴线垂直的平面）、T ''（过
点B 与轴线垂直的平面）作为校正平面， 在T '、T ''平面内分别加上校正质量b m '、b
m ''，矢径为b r '、b r ''，则校正质量所产生的离心惯性力为2ωb b b r m P ''='和2ωb b b r m P ''''=''，1P 、2P 、3P 、b P '和b P ''组成了空间力系。

选取三坐标轴x 、y 、z 轴如图所示，并将作用在转子上的所有力向YAZ 平面和XAY 平面投影，如图3-3所示。

图3-3
在图3-3中, 所有的力组成了平面平行力系, 列平衡方程：
0=⎪⎭
⎫ ⎝⎛∑→F m A ，02211=''+'-'l P l P l P bz z z （3-4） 0=∑z
F
，021=''+-+'bz z z bz
P P P P （3-5） 解得：
l
l P l P P z z bz
112
2'-''='' （3-6） z z z bz
P P P P 212+-''-=' （3-7） 式中：
11P P z -在Z 轴上的投影，111cos αP P z =，N ； 22P P z -在Z 轴上的投影，22cos αP P z =，N ；
b bz
P P '-'在Z 轴上的投影，β''='cos b bz P P ，N ； b bz
P P ''-''在Z 轴上的投影，β''''=''cos b bz P P ，N ； 同理，在XAY 平面内
0=⎪⎭
⎫ ⎝⎛∑→F m A ，0332211=''+'-'-'l P l P l P l P bx x x x （3-8） 0=∑z
F
，0321=''++--bx
x x x bx P P P P P （3-9） 解得：
l
l P l P l P P x x x bx
332
211'-'+'='' （3-10）
bx x x x bx
P P P P P ''-++='321 （3-11） 式中：
11P P x -在X 轴上的投影，111sin αP P x =，N ； 22P P x -在X 轴上的投影，222sin αP P x =，N ； 33P P x -在X 轴上的投影，33P P x =，N ； b bx
P P '-'在X 轴上的投影，β''=sin b bx P P ，N ；
b bx
P P ''-''在X 轴上的投影，β''''=''sin b bx P P ，N ； 在校正平面T '上，校正质量b
m '所产生的离心惯性力b P '： 22bz bx
b P P P '+'=' （3-12） bz
bx
P P ''='-1
tan β （3-13） 在校正平面T ''上，校正质量b
m ''所产生的离心惯性力b P ''： 22bz bx
b P P P ''+''='' （3-14） bz
bx
P P ''''=''-1
tan β （3-15） 当给定校正质量
b
m '和
b
m ''的失径b r 'ρ和b r ''ρ，很容易求得校正质量b m '和b m '
'的数值：
2ωb b b
r P m ''
=' （3-16）
2
ωb b b
r P m ''''='' （3-17）
b
m '和
b
m ''的方位由角度β'和β''确定，如图3-3所示。

因为方程每一项均有2ω，
故2
ω可以约去而不必计算。