多体动力学仿真算例

2 1 2
4
BUAA
一般情况下系统 动力学基本关系
小圆盘质心切向力：
多体动力学数值仿真大作业2
二、用牛顿欧拉方法建模： 由动量矩定理：
J11 Fx L k1 c11 M(t )
J 22 M(t ) c22
y
M
m2 L1 Fx m2 g sin 1
即铰链约束力切向分量：
M
拉格朗日函数：
L 1 1 1 1 2 2 m1 R1212 m2 L212 m2 R2 2 k12 m2 g cos 1 4 2 4 2
mg
O
x
非有势力的广义力： Q 1 c11 M ( t )
Q 2 M (t ) c2 2
图2 双圆盘系统主动力图
14
BUAA
未修正整理得：
多体动力学数值仿真大作业2
违约修正
1 ( A B ) 1 (c11 k1 c2 2
m2 gL sin 1 B w sin wt )
1 2
解得：
( A B) 1 ( 2 max sin t
2 1 w2max sin wt
2
max 1.8rad,
初始条件：
2.0rad/s
M(N*m)
1
0
1 0,1 1,1 0,1 1
粉、绿、蓝线分别对应1、2、3情况
-1
-2
0
2
4
6
8
10
t(s)
2015-6-12 16
BUAA
大圆盘转角与阻尼系数
系统参数：
多体动力学数值仿真大作业2
2
c1 0.2Nms/rad, c2 0.5Nms/rad; c1 0.5Nms/rad, c2 0.5Nms/rad; c1 0.5Nms/rad, c2 1Nms/rad.
给静稳定点
变化如图中红线所示。蓝线表示无
初角速度的情形。
1.6 1.4 1.2
11 3，圆盘1转角变化如图中
2.2 2 1.8 1.6
11
1.2757一个大扰
红线所示。蓝线是无初角速度情形。
1源自文库
1(rad)
0.8 0.6 0.4
1(rad)
0 5 10 15 20
1.4 1.2 1 0.8 0.6
BUAA
多体动力学数值仿真大作业2
双圆盘系统运动 的动力学特性分析
姓名：XXX 专业：XXXXX
2015-6-12
1
BUAA
多体动力学数值仿真大作业2
质量为m1的均质圆盘1与刚度系数为k的扭簧连接绕轴O转动， 质量为m2的均质圆盘2由圆盘1上的电机驱动绕C2轴转动，两圆 盘转动时受到与其角速度成正比的阻尼力偶作用，阻尼系数为 c1和c2。 两轴间的距离为L。
0.2 0
0.4 0.2
0
5
10
15
20
t(s)
t(s)
2015-6-12
8
BUAA
五、给定驱动力偶
多体动力学数值仿真大作业2
水平状况下解析解与数值解对照
给定周期驱动力偶下，大圆盘的运动方程如下： 若令g=0,即圆盘水平，运动方程是可以求出解析解的。取c1=0.2，有： 当 1 0,1 0 ，解得：
（5）
B12 c22
其中：
( 2 1 max cos t ) ( 2 1 max sin t )
c2 2 B (k1 c11 ) A)
A 2(m R 2m2 L )
2 1 1
2 1 B 2(m2 R2 )
代入数值：
（4）
把平衡解代入：
2V 2V 1.5 2 cos1.2757 0 2 2 11 12 2V 1.5 2 0 2 13
2sin 1 1.51 0
用fzero函数解非线性方程，得3个平衡点：
11 1.2757,12 1.2757,13 0
3
2015-6-12
BUAA
整理得：
多体动力学数值仿真大作业2
一般情况下运动微分方程
1 A(c11 k1 M(t ) m2 gL sin1 )
2 B(M(t ) c22 )
其中：
（2）
A 2(m1R12 2m2 L2 )1
B 2(m2 R )
2015-6-12
多体动力学数值仿真大作业2
五、给定运动约束条件
约束条件： f 2 1 max sin t 0

拉格朗日方程与约束条件联立：：
1 m1 R121 m2 L21 k1 c11 m2 gl sin 1 2
广义力： Q 11 k11 c111 m2 gl sin 11
初始条件：
2.0rad/s
0 -0.5 -1 -1.5 -2
1 0,1 0,2 0,2 0
1 0,1 0.001,2 0, 2 0
蓝线表示无初角速度的情况，红线
0
10
20
30
40
50
60
t(s)
表示角速度为0.001时的情况。 微小扰动后20秒后出现了大差异，可判断出现了混沌现象
（3）
Fx m2 g sin 1 m2 L1
又圆盘转动惯量：
2 R2
1 J1 m1 R12 2 1 2 J 2 m2 R2 2
将约束力切向分量和转动惯量均代入式 子（3），易看出与式子（1）完全相符。 即两种建模方法的结果是相吻合的。
2015-6-12
R1
1
O
x
图3 双圆盘系统受力分析图
一、用第二类拉格朗日方程建模
1 m1 R121 m2 L21 m2 gL sin 1 k1 c11 M ( t ) 2

1 2 m2 R2 2 M ( t ) c2 2 2
y
势能：
V 1 2 k1 m2 gL cos 1 2
R1
2 R2 1
( J1 m2 L2 )1 c11 k1 mgL sin 1 M maxsint
0.8 0.6 0.4 0.2
0.161 0.21 1.51 0.5sin2t
1 e0.5t (0.3077cos 2.6926t
1(rad)
0.4571sin 2.6926t ) -0.5385sin2t 0.3077cos2t
Q 1 k1 c11 m2 gl sin 1
Q 2 2 c222 Q 2 c2 2
2015-6-12
1 2 m2 R2 2 c2 2 2
（4）
f 2 1 max sin t 0
解出的 即圆盘的驱动力偶矩
1 2 V k1 m2 gL cos 1 势能： 2
2015-6-12
故静稳定点为：
静不稳点为：
11 ,12
13
7
BUAA
c1,c2分别取0.2和0.5
多体动力学数值仿真大作业2
0 一个微
动
给静不稳定平衡点
小的扰动
13
0.01，圆盘1转角
13
数值法验证平 衡点的稳定性
经试算，在c1取大致0.236-0.265范围内会出现此现象
2015-6-12 13
BUAA
用第一类拉格朗日 方程建模求约束力
动能：
T
2 1 1 1 J112 m2 L1 J1 22 2 2 2 1 1 1 m1 R1212 m2 L212 m2 R22 22 4 2 4
-1.4 -1.6 -1.8
0
5
10
15
20
t(s)
10
BUAA
大圆盘周期运动情况
系统参数：
多体动力学数值仿真大作业2
c1 0.35Nms/rad, c 2 0.5Nms/rad, M max 0.5rad,
初始条件：
0.2 0 -0.2
2.0rad/s
1(rad)
-0.4 -0.6 -0.8 -1 -1.2
i (k 1) i (k ) h i (k )
i (k 1) i (k ) (i (k 1) i (k )) h / 2
（4）令t(k)=t(k+1),判断t(k)是否达到终止条件tf， 若达到则终止计算，未达到则回到步骤（2）。
2015-6-12 6
BUAA
1 0,1 0,2 0,2 0
1 0,1 0.001 , 2 0, 2 0
蓝线表示无初角速度的情况，红线 表示角速度为0.001时的情况。
2015-6-12
-1.4 -1.6 -1.8
0
5
10
15
20
25
30
t(s)
11
BUAA
大圆盘周期运动情况
系统参数：
蓝线表示无初角速度的情况，红线 表示角速度为0.001时的情况。
2015-6-12
0
10
20
30
40
50
60
t(s)
12
BUAA
大圆盘混沌运动情况
系统参数：
多体动力学数值仿真大作业2
2 1.5 1 0.5
c1 0.25Nms/rad, c 2 0.5Nms/rad,
1(rad)
M max 0.5rad,
2 1
违约修正方程：
d2 f df f 0 2 dt dt
2015-6-12
1 A( c11 k1 m2 gL sin 1 )
（6）
2 B( c22 )
15
BUAA
驱动力偶与阻尼系数
系统参数：
多体动力学数值仿真大作业2
4
c1 0.2Nms/rad, c2 0.5Nms/rad; 3 c1 0.5Nms/rad, c2 0.5Nms/rad; c1 0.5Nms/rad, c2 1Nms/rad.
多体动力学数值仿真大作业2
解析求系统静平 衡点及其稳定性
对位移求二阶导有：
2V k m2 gL cos 1 2 1 1.5 2cos 1
四、不加驱动对系统进行研究
当系统静止，则各角速度和角加速 度均为0，驱动力M也为0，则由式 子（1）有：
m2 gL sin 1 k1 0
2015-6-12
1
2 R2
C2
O
x
图1 双圆盘系统位形坐标图
2
BUAA
一般情况下系统 动力学基本关系
动能：
T
2 1 1 1 J112 m2 L1 J1 22 2 2 2 1 1 1 m1 R1212 m2 L212 m2 R22 22 4 2 4
多体动力学数值仿真大作业2
多体动力学数值仿真大作业2
2 1.5 1 0.5
c1 0.2Nms/rad, c 2 0.5Nms/rad, M max 0.5rad,
初始条件：
2.0rad/s
1(rad)
0 -0.5 -1 -1.5 -2
1 0,1 0,2 0,2 0
1 0,1 0.001,2 0, 2 0
数值解和解析解的对比如图， 红线为数值解，蓝线为解析解
0 -0.2 -0.4 -0.6 -0.8
0
1
2
3
4
5
6
7
8
9
10
t(s)
2015-6-12
9
BUAA
大圆盘周期运动情况
系统参数：
多体动力学数值仿真大作业2
c1 0.5Nms/rad, c 2 0.5Nms/rad, M max 0.5rad,
设系统参数： m
1
2.0kg,
m2 1.0kg,
R1
y
R1 0.4m, R2 0.1m, OC2 L 0.2m, k 1.5Nm/rad, c1 0.2 0.6Nms/rad, c 2 0.5 1.0Nms/rad,
研究目标：用解析法和数值法分析其动 力学特性
5
BUAA
三、数值算法
多体动力学数值仿真大作业1
x1 1, x2 1, x3 2 , x4 2
（1）设定初值x（0）、时间步长h、时间上限tf （2）根据式子（2）把x（k）代入求 1 , 2 （3）增加步长：
t (k 1) t (k ) h
相应有角速度： 角度：
初始条件：
0.2 0 -0.2
2.0rad/s
1(rad)
-0.4 -0.6 -0.8 -1 -1.2
1 0,1 0,2 0,2 0
1 0,1 0.001,2 0, 2 0
蓝线表示无初角速度的情况，红线 表示角速度为0.001时的情况。
2015-6-12
(rad)
1.5
1
max 1.8rad,