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

弹簧阻尼四连杆机构动力学分析—Baumgarte违约修正弹簧阻尼四连杆机构动力学分析——Baumgarte违约修正目录一、动力学分析 21、外力矢量 22、外力所做的虚功 33、系统的约束方程 34、雅克比矩阵Cq 45、Qd 46、系统的动力学方程 47、Baumgarte违约修正法 5二、说明 7三、Matlab与ADAMS仿真曲线对比 8仿真分析 12附录一：求雅克比矩阵的Matlab程序 19附录二：求Qd的Matlab程序 19附录三：用欧拉迭代法+Baumgarte违约修正进行动力学分析的Matlab程序 20

使用ADAMS建立如图1所示的四连杆机构，二杆长150mm，三杆长500mm，四杆长450mm，二杆所受的转矩为100N∙m；三个杆都受重力，二杆初始角度为90度；在四杆的质心和地面之间加上弹簧和阻尼，弹簧原长为416.08mm，弹簧刚度为4N/mm，阻尼为0.4N/(mm/s)。运用欧拉法，并加入Baumgarte违约修正法，通过Baumgarte违约修正法可以明显看到各仿真曲线更接近于ADAMS的仿真曲线，matlab程序见文件夹中的sigan_donglixue_spring_Baumgarte.m文件。改变不同的∝和β值，可以得到不同的修正结果，见sigan_donglixue_spring_Baumgarte_duibi.m文件，通过改变不同的∝转矩弹簧阻尼四杆三杆二杆转矩弹簧阻尼四杆三杆二杆图1加弹簧阻尼四杆机构一、动力学分析1、外力矢量（1）弹簧力为：fs=k*(sqrt((q(7,i)-750)^2+(q(8,i)-0)^2)-L0);（2）阻尼力为：fc=-c*q_v(7,i)*(750-q(7,i))/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2)+c*q_v(8,i)*(q(8,i)-0)/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2)（3）四杆质心x轴方向所受的外力：Fx=fs*(750-q(7,i))/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2)+fc*(750-q(7,i))/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2);（4）四杆质心y轴方向所受的外力：Fy=-m4*g-fs*(q(8,i)-0)/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2)-fc*(q(8,i)-0)/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2);2、外力所做的虚功3、系统的约束方程R4、雅克比矩阵Cq由系统的约束方程，求雅克比矩阵得Cq。具体程序见附录一，可得雅克比矩阵C5、Qd具体程序见附录二，可得Qd6、系统的动力学方程其中：M=7、Baumgarte违约修正法在上述系统的动力学方程中加入Baumgarte违约修正法，如下式。C为系统的约束方程，Cqq为对系统约束方程C求雅克比矩阵，再乘以q。C系统的约束方程为：R即C为：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);x4x5x6x7x8x9x10x11x12分别代表二杆、三杆、四杆的位移量。Cqy4+(l2*y6*sin(x6))/2y5-(l2*y6*cos(x6))/2y4-y7-(l2*y6*sin(x6))/2-(l3*y9*sin(x9))/2y5-y8+(l2*y6*cos(x6))/2+(l3*y9*cos(x9))/2y7-y10-(l3*y9*sin(x9))/2-(l4*y12*sin(x12))/2y8-y11+(l3*y9*cos(x9))/2+(l4*y12*cos(x12))/2y10-(l4*y12*sin(x12))/2y11+(l4*y12*cos(x12))/2y4y5y6y7y8y9y10y11y12分别代表二杆、三杆、四杆的速度量。将x、y转化成编程所用的数组q(i)、q_v(i)得到：C为：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));利用yakebi.mq_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))Cqq、C分别乘以-2∝-(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)))-b^2*(q(8,i)+l4/2*sin(q(9,i)));];使用欧拉迭代法加Baumgarte违约修正法就可以计算出各个曲线，并改变不同的∝和β值观察对仿真结果的影响，具体Matlab程序见文件夹中的sigan_donglixue_spring_Baumgarte_duibi.m和二、说明1、由于本程序采用欧拉迭代法求解，所以计算过程有误差。当仿真时间越长时，与ADAMS曲线相差越大，但曲线形状一致，只有相位滞后。通过Baumgarte违约修正法可以明显看到各仿真曲线更接近于ADAMS的仿真曲线，改变不同的∝和β值，可以得到不同的修正结果。通过2、由于在ADAMS中用常规方法无法测量三杆的角度。所以再建立两个辅助杆，并设这两个辅助杆的质量和转动惯量都为零，用于测量三杆的角度，如图2中辅助杆1和辅助杆2。再就可以测量由E、F、G点组成的角度作为三杆的角度。辅助杆1辅助杆2点E点F点G辅助杆1辅助杆2点E点F点G图2辅助测量杆三、Matlab与ADAMS仿真曲线对比局部放大局部放大图3弹簧阻尼器所受的力图4上图局部放大的曲线对比局部放大局部放大图5二杆位移变化曲线图6上图局部放大的曲线对比局部放大局部放大图7铰链2处作用力变化曲线图8上图局部放大的曲线对比仿真分析由图4、6、8可以看出：当∝=5，无修正的绿色曲线与ADAMS相差最大；∝=1，β=10时比无修正时好，但比∝=5可以得出结论Baumgarte违约修正法，如果选取∝、图9铰链2处作用力变化曲线图10铰链3处作用力变化曲线图11铰链4处作用力变化曲线图12三杆质心位置变化曲线图13四杆质心位置变化曲线图14二杆质心速度变化曲线图15三杆质心速度变化曲线图16四杆质心速度变化曲线图17二杆质心加速度变化曲线图18三杆质心加速度变化曲线图19四杆质心加速度变化曲线

附录一：求雅克比矩阵的Matlab程序clcclearsymsl2l3l4l5;symsm2m3m4;symsx4x5x6x7x8x9x10x11x12g=[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,[x4x5x6x7x8x9x10x11x12])m=[m200000000;0m20000000;00m2/12*(l2)^2000000;000m300000;0000m30000;00000m3/12*(l3)^2000;000000m400;0000000m40;00000000m4/12*(l4)^2]o=zeros(8,8);cq1=cq.';A=[mcq1;cqo]附录二：求Qd的Matlab程序clcclearsymsl2l3l4l5;symsm2m3m4;symsx4x5x6x7x8x9x10x11x12symsy4y5y6y7y8y9y10y11y12h=[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,[x4x5x6x7x8x9x10x11x12]);q=[y4y5y6y7y8y9y10y11y12].';qd=-qd1*q附录三：用欧拉迭代法+Baumgarte违约修正进行动力学分析的Matlab程序%对加弹簧阻尼四杆机构进行动力学分析clearclcm2=1.1289713258;%三个杆的质量kgm3=3.3132513258;m4=3.0012113258;g=9806.65;%重力加速度9806.65mm/s2l2=150;%二杆长度mml3=500;%三杆长度mml4=450;%四杆长度mml5=400;%两个底座之间的距离mmM=10000000;%二杆扭矩为100000N.mk=4000;%弹簧刚度为4N/mmc=400;%阻尼为0.5L0=416.08;%弹簧原长mmdt=0.00001;%步长t=0:dt:5;%设置时间为5sn=length(t)%仿真步数q=zeros(9,n);%产生9行n列的零矩阵q_v=zeros(9,n);q_ac=zeros(9,n);q(1:9,1)=[075pi/22003000.6435011094002251.5*pi];%设置迭代初值，也是三个杆的初始位置force=zeros(8,n);F=zeros(1,n);%%用欧拉法对四连杆机构进行位置分析，速度分析，加速度分析。动力学方程的左边矩阵，AX=B，求解XA=[MCqT;Cq0]fori=1:nA=[m2,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,00,m2,0,0,0,0,0,0,0,0,1,0,1,0,0,0,00,0,(l2^2*m2)/12,0,0,0,0,0,0,(l2*sin(q(3,i)))/2,-(l2*cos(q(3,i)))/2,-(l2*sin(q(3,i)))/2,(l2*cos(q(3,i)))/2,0,0,0,00,0,0,m3,0,0,0,0,0,0,0,-1,0,1,0,0,00,0,0,0,m3,0,0,0,0,0,0,0,-1,0,1,0,00,0,0,0,0,(l3^2*m3)/12,0,0,0,0,0,-(l3*sin(q(6,i)))/2,(l3*cos(q(6,i)))/2,-(l3*sin(q(6,i)))/2,(l3*cos(q(6,i)))/2,0,00,0,0,0,0,0,m4,0,0,0,0,0,0,-1,0,1,00,0,0,0,0,0,0,m4,0,0,0,0,0,0,-1,0,10,0,0,0,0,0,0,0,(l4^2*m4)/12,0,0,0,0,-(l4*sin(q(9,i)))/2,(l4*cos(q(9,i)))/2,-(l4*sin(q(9,i)))/2,(l4*cos(q(9,i)))/21,0,(l2*sin(q(3,i)))/2,0,0,0,0,0,0,0,0,0,0,0,0,0,00,1,-(l2*cos(q(3,i)))/2,0,0,0,0,0,0,0,0,0,0,0,0,0,01,0,-(l2*sin(q(3,i)))/2,-1,0,-(l3*sin(q(6,i)))/2,0,0,0,0,0,0,0,0,0,0,00,1,(l2*cos(q(3,i)))/2,0,-1,(l3*cos(q(6,i)))/2,0,0,0,0,0,0,0,0,0,0,00,0,0,1,0,-(l3*sin(q(6,i)))/2,-1,0,-(l4*sin(q(9,i)))/2,0,0,0,0,0,0,0,00,0,0,0,1,(l3*cos(q(6,i)))/2,0,-1,(l4*cos(q(9,i)))/2,0,0,0,0,0,0,0,00,0,0,0,0,0,1,0,-(l4*sin(q(9,i)))/2,0,0,0,0,0,0,0,00,0,0,0,0,0,0,1,(l4*cos(q(9,i)))/2,0,0,0,0,0,0,0,0];fs=k*(sqrt((q(7,i)-750)^2+(q(8,i)-0)^2)-L0);%弹簧力fc=-c*q_v(7,i)*(750-q(7,i))/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2)+c*q_v(8,i)*(q(8,i)-0)/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2);%阻尼力F(1,i)=fs+fc;%弹簧力加阻尼力B=[0;-m2*g;M;0;-m3*g0;fs*(750-q(7,i))/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2)+fc*(750-q(7,i))/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2);%四杆质心所受的x轴方向的力-m4*g-fs*(q(8,i)-0)/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2)-fc*(q(8,i)-0)/sqrt((q(7,i)-750)^2+(q(8,i)-0)^2);%四杆质心所受的y轴方向的力0;-(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)))-b^2*(q(8,i)+l4/2*sin(q(9,i)));];temp=A\B;q_ac(:,i)=temp(1:9);%解的前9项为杆的加速度和角加速度force(:,i)=temp(10:17);%解的后8项为铰链处的力ifi==n;breakendq_v(:,i+1)=q_v(:,i)+dt*q_ac(:,i);q(:,i+1)=q(:,i)+dt*q_v(:,i);endh0=figure(14);set(h0,'name','弹簧阻尼器所受的力');plot(t,-F(1,:)/1000)holdondata0=load('force.tab');plot(data0(:,1),data0(:,2),'r-');xlabel('时间t/s')ylabel('力/N')title('弹簧阻尼器所受的力')legend('matlab','adams')%%二杆质心位移%二杆质心X位移变化h1=figure(1);set(h1,'name','二杆质心位移');subplot(311)plot(t,q(1,:),'b-')holdon;data1=load('gan2_xy.tab');plot(data1(:,1),data1(:,2),'r-');xlabel('时间t/s')ylabel('位移/m')title('二杆质心X位移变化')legend('matlab','adams')%二杆质心Y位移变化subplot(312)plot(t,q(2,:),'b-')holdon;plot(data1(:,1),data1(:,3),'r-');xlabel('时间t/s')ylabel('位移/m')title('二杆质心Y位移变化')legend('matlab','adams')%二杆质心角度变化subplot(313)plot(t,q(3,:),'b-')holdon;data0=load('gan234_a.tab');plot(data0(:,1),data0(:,2),'r-');xlabel('时间t/s')ylabel('角度rad')title('二杆质心角度变化')legend('matlab','adams')%%三杆质心位移%三杆质心X位移变化h2=figure(2);set(h2,'name','三杆质心位移');subplot(311)plot(t,q(4,:),'b-')holdon;data2=load('gan3_xy.tab');plot(data2(:,1),data2(:,2),'r-');xlabel('时间t/s')ylabel('位移/m')title('三杆质心X位移变化')legend('matlab','adams')%三杆质心Y位移变化subplot(312)plot(t,q(5,:),'b-')holdon;plot(data2(:,1),data2(:,3),'r-');xlabel('时间t/s')ylabel('位移/m')title('三杆质心Y位移变化')legend('matlab','adams')%三杆质心角度变化subplot(313)plot(t,q(6,:),'b-')holdon;plot(data0(:,1),data0(:,3),'r-');xlabel('时间t/s')ylabel('角度rad')title('三杆质心角度变化')legend('matlab','adams')%%四杆质心位移%四杆质心X位移变化h3=figure(3);set(h3,'name','四杆质心位移');subplot(311)plot(t,q(7,:),'b-')holdon;data3=load('gan4_xy.tab');plot(data3(:,1),data3(:,2),'r-');xlabel('时间t/s')ylabel('位移/m')title('四杆质心X位移变化')legend('matlab','adams')%四杆质心Y位移变化subplot(312)plot(t,q(8,:),'b-')holdon;plot(data3(:,1),data3(:,3),'r-');xlabel('时间t/s')ylabel('位移/m')title('四杆质心Y位移变化')legend('matlab','adams')%四杆质心角度变化subplot(313)plot(t,q(9,:),'b-')holdon;plot(data0(:,1),data0(:,4)+pi,'r-');xlabel('时间t/s')ylabel('角度rad')title('四杆质心角度变化')legend('matlab','adams')%%二杆质心速度变化%二杆质心X速度变化h4=figure(4);set(h4,'name','二杆质心速度');subplot(311)plot(t,q_v(1,:),'b-')holdon;data4=load('gan2v_xya.tab');plot(data4(:,1),data4(:,2),'r-');xlabel('时间t/s')ylabel('速度m/s')title('二杆质心X速度变化')legend('matlab','adams')%二杆质心Y速度变化subplot(312)plot(t,q_v(2,:),'b-')holdon;plot(data4(:,1),data4(:,3),'r-');xlabel('时间t/s')ylabel('速度m/s')title('二杆质心Y速度变化')legend('matlab','adams')%二杆质心角速度变化subplot(313)plot(t,q_v(3,:),'b-')holdon;plot(data4(:,1),data4(:,4),'r-');xlabel('时间t/s')ylabel('角速度rad/s')title('二杆质心角速度变化')legend('matlab','adams')%%三杆质心速度变化%三杆质心X速度变化h5=figure(5);set(h5,'name','三杆质心速度');subplot(311)plot(t,q_v(4,:),'b-')holdon;data5=load('gan3v_xya.tab');plot(data5(:,1),data5(:,2),'r-');xlabel('时间t/s')ylabel('速度m/s')title('三杆质心X速度变化')legend('matlab','adams')%三杆质心Y速度变化subplot(312)plot(t,q_v(5,:),'b-')holdon;plot(data5(:,1),data5(:,3),'r-');xlabel('时间t/s')ylabel('速度m/s')title('三杆质心Y速度变化')legend('matlab','adams')%三杆质心角速度变化subplot(313)plot(t,q_v(6,:),'b-')holdon;plot(data5(:,1),data5(:,4),'r-');xlabel('时间t/s')ylabel('角速度rad/s')title('三杆质心角速度变化')legend('matlab','adams')%%四杆质心速度变化%四杆质心X速度变化h6=figure(6);set(h6,'name','四杆质心速度');subplot(311)plot(t,q_v(7,:),'b-')holdon;data6=load('gan4v_xya.tab');plot(data6(:,1),data6(:,2),'r-');xlabel('时间t/s')ylabel('速度m/s')title('四杆质心X速度变化')legend('matlab','adams')%四杆质心Y速度变化subplot(312)plot(t,q_v(8,:),'b-')holdon;plot(data6(:,1),data6(:,3),'r-');xlabel('时间t/s')ylabel('速度m/s')title('四杆质心Y速度变化')legend('matlab','adams')%四杆质心角速度变化subplot(313)plot(t,q_v(9,:),'b-')holdon;plot(data6(:,1),data6(:,4),'r-');xlabel('时间t/s')ylabel('角速度rad/s')title('四杆质心角速度变化')legend('matlab','adams')%%二杆质心加速度%二杆质心X加速度变化h7=figure(7);set(h7,'name','二杆质心加速度');subplot(311)plot(t,q_ac(1,:),'b-')holdon;data7=load('gan2ac_xya.tab');plot(data7(:,1),data7(:,2),'r-');xlabel('时间t/s')ylabel('加速度m/s^2')title('二杆质心X加速度变化')legend('matlab','adams')%二杆质心Y加速度变化subplot(312)plot(t,q_ac(2,:),'b-')holdon;plot(data7(:,1),data7(:,3),'r-');xlabel('时间t/s')ylabel('加速度m/s^2')title('二杆质心Y加速度变化')legend('matlab','adams')%二杆质心角加速度变化subplot(313)plot(t,q_ac(3,:),'b-')holdon;plot(data7(:,1),data7(:,4),'r-');xlabel('时间t/s')ylabel('角加速度rad/s^2')title('二杆质心角加速度变化')legend('matlab','adams')%%三杆质心加速度%三杆质心X加速度变化h8=figure(8);set(h8,'name','三杆质心加速度');subplot(311)plot(t,q_ac(4,:),'b-')holdon;data8=load('gan3ac_xya.tab');plot(data8(:,1),data8(:,2),'r-');xlabel('时间t/s')ylabel('加速度m/s^2')title('三杆质心X加速度变化')legend('matlab','adams')%三杆质心Y加速度变化subplot(312)plot(t,q_ac(5,:),'b-')holdon;plot(data8(:,1),data8(:,3),'r-');xlabel('时间t/s')ylabel('加速度m/s^2')title('三杆质心Y加速度变化')legend('matlab','adams')%三杆质心角加速度变化subplot(313)plot(t,q_ac(6,:),'b-')holdon;plot(data8(:,1),data8(:,4),'r-');xlabel('时间t/s')ylabel('角加速度rad/s^2')title('三杆质心角加速度变化')legend('matlab','adams')%%四杆质心加速度%四杆质心X加速度变化h9=figure(9);set(h9,'name','四杆质心加速度');subplot(311)plot(t,q_ac(4,:),'b-')holdon;data9=load('gan3ac_xya.tab');plot(data9(:,1),data9(:,2),'r-');xlabel('时间t/s')ylabel('加速度m/s^2')title('四杆质心X加速度变化')legend('matlab','adams')%四