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

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

图17 二杆质心加速度变化曲线

图18 三杆质心加速度变化曲线

17

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

图19 四杆质心加速度变化曲线

18

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

附录一：求雅克比矩阵的Matlab程序

clc clear

syms l2 l3 l4 l5; syms m2 m3 m4;

syms x4 x5 x6 x7 x8 x9 x10 x11 x12 g=[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); ]

cq=jacobian(g,[x4 x5 x6 x7 x8 x9 x10 x11 x12]) m=[m2 0 0 0 0 0 0 0 0; 0 m2 0 0 0 0 0 0 0;

0 0 m2/12*(l2)^2 0 0 0 0 0 0; 0 0 0 m3 0 0 0 0 0; 0 0 0 0 m3 0 0 0 0;

0 0 0 0 0 m3/12*(l3)^2 0 0 0; 0 0 0 0 0 0 m4 0 0; 0 0 0 0 0 0 0 m4 0;

0 0 0 0 0 0 0 0 m4/12*(l4)^2] o=zeros(8,8); cq1=cq.'; A=[m cq1;cq o]

附录二：求

clc clear

的Matlab程序

syms l2 l3 l4 l5; syms m2 m3 m4;

19

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

syms x4 x5 x6 x7 x8 x9 x10 x11 x12 syms y4 y5 y6 y7 y8 y9 y10 y11 y12 h=[y4+y6*l2/2*sin(x6); y5-y6*l2/2*cos(x6);

y4-y6*l2/2*sin(x6)-y7-y9*l3/2*sin(x9); y5+y6*l2/2*cos(x6)-y8+y9*l3/2*cos(x9); y7-y9*l3/2*sin(x9)-y10-y12*l4/2*sin(x12); y8+y9*l3/2*cos(x9)-y11+y12*l4/2*cos(x12); y10-y12*l4/2*sin(x12); y11+y12*l4/2*cos(x12);; ]

qd1=jacobian(h,[x4 x5 x6 x7 x8 x9 x10 x11 x12]); q=[y4 y5 y6 y7 y8 y9 y10 y11 y12].'; qd=-qd1*q

附录三：用欧拉迭代法+Baumgarte违约修正进行动力学分析的Matlab程序

%对加弹簧阻尼四杆机构进行动力学分析 clear clc

m2=1.1289713258; %三个杆的质量kg m3=3.3132513258;

m4=3.0012113258;

g=9806.65; %重力加速度9806.65mm/s2 l2=150; %二杆长度mm l3=500; %三杆长度mm l4=450; %四杆长度mm

l5=400; %两个底座之间的距离mm M=10000000; %二杆扭矩为100000N.m k=4000; %弹簧刚度为4N/mm c=400; %阻尼为0.5 L0=416.08; %弹簧原长mm dt=0.00001; %步长

t=0:dt:5; %设置时间为5s n=length(t) %仿真步数

q=zeros(9,n); %产生9行n列的零矩阵 q_v=zeros(9,n); q_ac=zeros(9,n);

q(1:9,1)=[0 75 pi/2 200 300 0.643501109 400 225 1.5*pi]; %设置迭代初值，也是三个杆的初始位置

20