张海涛-2013424064-第四次作业+Baumgarte违约修正

弹簧阻尼四连杆机构动力学分析—Baumgarte违约修正

Qe??0?m1gM0?m2g0FxFy0?

T7、Baumgarte违约修正法

在上述系统的动力学方程中加入Baumgarte违约修正法，如下式。为系统的约束方程，为对

求雅克比矩阵，再乘以。为系统约束方程直接对时间求导。

系统的约束方程为：

5

弹簧阻尼四连杆机构动力学分析—Baumgarte违约修正

即为：

x4-l2/2*cos(x6);

x5-l2/2*sin(x6);

x4+l2/2*cos(x6)-x7+l3/2*cos(x9); x5+l2/2*sin(x6)-x8+l3/2*sin(x9); x7+l3/2*cos(x9)-x10+l4/2*cos(x12); x8+l3/2*sin(x9)-x11+l4/2*sin(x12); x10+l4/2*cos(x12)-l5; x11+l4/2*sin(x12);

x4 x5 x6 x7 x8 x9 x10 x11 x12分别代表二杆、三杆、四杆的位移量。

为：

y4 y5 y6 y7 y8 y9 y10 y11 y12分别代表二杆、三杆、四杆的速度量。

将x、y转化成编程所用的数组q(i)、q_v(i)得到： 为：

q(1,i)-l2/2*cos(q(3,i));

q(2,i)-l2/2*sin(q(3,i));

q(1,i)+l2/2*cos(q(3,i))-q(4,i)+l3/2*cos(q(6,i)); q(2,i)+l2/2*sin(q(3,i))-q(5,i)+l3/2*sin(q(6,i)); q(4,i)+l3/2*cos(q(6,i))-q(7,i)+l4/2*cos(q(9,i)); q(5,i)+l3/2*sin(q(6,i))-q(8,i)+l4/2*sin(q(9,i)); q(7,i)+l4/2*cos(q(9,i))-l5; q(8,i)+l4/2*sin(q(9,i));

6

弹簧阻尼四连杆机构动力学分析—Baumgarte违约修正

为：

q_v(1,i)+q_v(3,i)*l2/2*sin(q(3,i)); q_v(2,i)-q_v(3,i)*l2/2*cos(q(3,i));

q_v(1,i)-q_v(3,i)*l2/2*sin(q(3,i))-q_v(4,i)- q_v(6,i)*l3/2*sin(q(6,i));

q_v(2,i)+ q_v(3,i)*l2/2*cos(q(3,i))-q_v(5,i)+ q_v(6,i)*l3/2*cos(q(6,i));

q_v(4,i)- q_v(6,i)*l3/2*sin(q(6,i))-q_v(7,i)- q_v(9,i)*l4/2*sin(q(9,i));

q_v(5,i)+ q_v(6,i)*l3/2*cos(q(6,i))-q_v(8,i)+ q_v(9,i)*l4/2*cos(q(9,i));

q_v(7,i)-q_v(9,i)*l4/2*sin(q(9,i)); q_v(8,i)+q_v(9,i)*l4/2*cos(q(9,i))

分别乘以和加到

-(l2*q_v(3,i)^2*cos(q(3,i)))/2-2*a*(q_v(1,i)+q_v(3,i)*l2/2*sin(q(3,i)))-b^2*(q(1,i)-l2/2*cos(q(3,i)));

-(l2*q_v(3,i)^2*sin(q(3,i)))/2-2*a*(q_v(2,i)-q_v(3,i)*l2/2*cos(q(3,i)))-b^2*(q(2,i)-l2/2*sin(q(3,i)));

(l2*cos(q(3,i))*q_v(3,i)^2)/2+(l3*cos(q(6,i))*q_v(6,i)^2)/2-2*a*(q_v(1,i)-q_v(3,i)*l2/2*sin(q(3,i))-q_v(4,i)-q_v(6,i)*l3/2*sin(q(6,i)))-b^2*(q(1,i)+l2/2*cos(q(3,i))-q(4,i)+l3/2*cos(q(6,i)));

(l2*sin(q(3,i))*q_v(3,i)^2)/2+(l3*sin(q(6,i))*q_v(6,i)^2)/2-2*a*(q_v(2,i)+q_v(3,i)*l2/2*cos(q(3,i))-q_v(5,i)+q_v(6,i)*l3/2*cos(q(6,i)))-b^2*(q(2,i)+l2/2*sin(q(3,i))-q(5,i)+l3/2*sin(q(6,i)));

(l3*cos(q(6,i))*q_v(6,i)^2)/2+(l4*cos(q(9,i))*q_v(9,i)^2)/2-2*a*(q_v(4,i)-q_v(6,i)*l3/2*sin(q(6,i))-q_v(7,i)-q_v(9,i)*l4/2*sin(q(9,i)))-b^2*(q(4,i)+l3/2*cos(q(6,i))-q(7,i)+l4/2*cos(q(9,i)));

(l3*sin(q(6,i))*q_v(6,i)^2)/2+(l4*sin(q(9,i))*q_v(9,i)^2)/2-2*a*(q_v(5,i)+q_v(6,i)*l3/2*cos(q(6,i))-q_v(8,i)+q_v(9,i)*l4/2*cos(q(9,i)))-b^2*(q(5,i)+l3/2*sin(q(6,i))-q(8,i)+l4/2*sin(q(9,i)));

(l4*q_v(9,i)^2*cos(q(9,i)))/2-2*a*(q_v(7,i)-q_v(9,i)*l4/2*sin(q(9,i)))-b^2*(q(7,i)+l4/2*cos(q(9,i))-l5);

(l4*q_v(9,i)^2*sin(q(9,i)))/2-2*a*(q_v(8,i)+q_v(9,i)*l4/2*cos(q(9,i))

7

弹簧阻尼四连杆机构动力学分析—Baumgarte违约修正

)-b^2*(q(8,i)+l4/2*sin(q(9,i)));];

使用欧拉迭代法加Baumgarte违约修正法就可以计算出各个曲线，并改变不同的

，具体

Matlab

程序见文件夹中的

sigan_donglixue_spring_Baumgarte_duibi.m和sigan_donglixue_spring_Baumgarte.m程序。

二、说明

1、由于本程序采用欧拉迭代法求解，所以计算过程有误差。当仿真时间越长时，与ADAMS曲线相差越大，但曲线形状一致，只有相位滞后。通过Baumgarte违约修正法可以明显看到各仿真曲线更接近于ADAMS的仿真曲线，改变不同的正结果。通过观察发现

，可以得到不同的修

仿真结果基本与ADAMS完全重合。

2、由于在ADAMS中用常规方法无法测量三杆的角度。所以再建立两个辅助杆，并设这两个辅助杆的质量和转动惯量都为零，用于测量三杆的角度，如图2中辅助杆1和辅助杆2。再就可以测量由E、F、G点组成的角度作为三杆的角度。

点G 点F 点E 辅助杆2 辅助杆1

图2 辅助测量杆

8