重庆大学矩阵论课程研究报告

“矩阵论”课程研究报告科目：矩阵理论及其应用教师：曾理姓名：ss 学号：20150xxxxxxx专业：机械工程类别：学术上课时间：2015 年9 月至2015 年12 月考生成绩：阅卷评语：阅卷教师(签名)平面四杆机构运动分析摘要： 采用计算机研究机械结构已成为现代机械设计分析的新方法。

本文采用矩阵法研究了平面四杆机构的运动情况，通过MATLAB 计算并输出计算结果及图表，对结果进行分析，得出机构的运行情况。

关键词：矩阵法，平面四杆机构，MATLAB一、问题描述机构运动分析的任务是在已知结构尺寸及原动件运动规律的情况下，确定机构中其他构件上某些点的轨迹、位移、速度及加速度和构件的角位移、角速度及角加速度。

上述这些内容，无论是设计新的机械，还是为了了解新得机械的运动性能，都是十分必要的，而且它还是研究机械动力性能的必要条件。

本文研究如下图2.1所示的四杆机构，其各构件的长度和原动件杆AB 的运动规律已知，即1θ为已知，而杆DA 固定，故4=0θ，故可求得两个未知方位角2θ及3θ。

同时给定1ω可求得在任一1θ下的杆BC 和杆CD 的角速度和角加速度。

二、方法简述机构运动分析的方法很多，主要有图解法和解析法。

采用解析法并借助计算机，不仅可获得很高的计算精度及一系列位置的分析结果，并能绘出机构相应的运动线图，同时还可把机构分析和机构综合问题联系起来，以便机构优化设计。

本文采取的矩阵法分析机构运动是一种利于计算机处理的方法，已成为现代机械机构分析的主要方法。

首先需列出机构的封闭矢量方程，如图2.1所示，先建立直角坐标系。

图2.1 机构封闭矢量图设杆AB 的长度为1l ，其方位角为1θ，1l 为杆AB 的矢量，即1l =AB。

机构中其余构件均可表示为相应的矢量，这样就形成由各矢量组成的的一个封闭矢量多边形，即ABCDA 。

在这个封闭矢量多边形中，有12430+--=l l l l（2.1）（1）位置分析将机构的封闭矢量方程向两坐标轴投影，得如下形式：2233411223311cos cos cos sin sin sin l l l l l l l θθθθθθ-=-⎫⎬-=-⎭（2.2）这个方程组为非线性二元方程组，代入已知量1θ，解此方程可得方位角3θ、2θ。

（2）速度分析将式2.2对时间取一次导数，可得222333111222333111sin sin sin cos cos cos l l l l l l ωθωθωθωθωθωθ-+=⎫⎬-=-⎭（2.3）写成矩阵形式为2233211122333111sin sin sin cos cos cos l l l l l l θθωωθθθωωθ-⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦通过上式可得12223311132233111sin sin sin cos cos cos l l l l l l ωθθωθωθθωθ--⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦ （2.4）（3）加速度分析将式2.4对时间t 求导，按照函数矩阵的求导法[1]则得22223332223321112222333322333111cos cos sin sin cos sin sin cos cos sin l l l l l l l l l l ωθωθωθθαωθωθωθωθθαωθ--⎡⎤⎡⎤⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦（2.5）整理得122223322233321112322332223333111sin sin cos cos cos cos cos sin sin sin l l l l l l l l l l αθθωθωθωωθαθθωθωθωωθ-⎛⎫--⎡⎤⎡⎤⎡⎤⎡⎤⎡⎤=- ⎪⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥ ⎪--⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦⎝⎭（2.6）三、实验数据和结果如下表所示参数为原始参数：表1 原始参数现需要求得杆BC ，杆CD 的角速度与角加速度变化图 采用MATAB 进行计算，取1θ=N ,=0,12,3360N ⋅⋅⋅,将10=θ，10.06AB l l ==m ，20.12BC l l ==m ，30.09CD l l ==m ，40.12CD l l ==m 代入式2.2，得23230.12cos 0.09cos 0.060.12sin 0.09sin 0-=⎫⎬-=⎭θθθθ解得246.5674=︒θ，3=75.5224︒θ将246.5674=︒θ，3=75.5224︒θ，110=ωrad/s 代入式2.4，得122233111322331111sin sin sin cos cos cos 0.12sin(46.5674)0.09sin(75.5224)0.0610sin(0)100.12cos(46.5674)-0.09cos(75.5224)0.0610cos(0)10l l l l l l ---⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦-︒︒⨯⨯-⎡⎤⎡⎤⎡⎤==⎢⎥⎢⎥⎢︒︒-⨯⨯-⎣⎦⎣⎦⎣⎦ωθθωθωθθωθ(/)rad s ⎥将246.5674=︒θ，3=75.5224︒θ，231010-⎡⎤⎡⎤=⎢⎥⎢⎥-⎣⎦⎣⎦ωω代入式2.6得122223322233321112322332223333111sin sin cos cos cos cos cos sin sin sin 0.12sin46.56740.09sin 75.5224=0.12cos 46.56740.09l l l l l l l l l l -⎛⎫--⎡⎤⎡⎤⎡⎤⎡⎤⎡⎤=- ⎪⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥ ⎪--⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦⎝⎭-⨯︒⨯︒⨯︒-⨯（）（）（）αθθωθωθωωθαθθωθωθωωθ1220.0610cos 0cos 75.52240.0610sin 00.12-10cos 46.56740.09-10cos 75.5224-1051.6397 -0.12-10sin 46.56740.09-10sin 75.5224-10189.3458-⎛⎡⎤⨯⨯⎡⎤ ⎢⎥⎢⎥ ︒⨯⨯⎣⎦⎣⎦⎝-⨯⨯︒⨯⨯︒⎫⎡⎤⎡⎤⎡=⎪⎢⎥⎢⎥⎢-⨯⨯︒⨯⨯︒⎣⎦⎣⎦⎣⎭（）（）（）（）（）（）（）2(/)rad s ⎤⎥⎦同理可得其他N 值对应的角速度，角加速度。

将所得数据绘制成曲线得如下图3.1所示杆BC 与CD 的角速度随时间变化图和图3.2杆BC 与CD 的角加速度随时间变化图。

图3.1 两杆角速度变化图图3.2 角加速度图四、结果分析与说明图3.1显示了杆BC 和杆CD 在原动件杆AB 的带动下其各自角速度随杆AB 的角位移的变化情况。

二者的均有先朝正方向增大，再向负方向增大的变化趋势，最终随着杆AB 完成一周的转动，二者的角速度返回初始值，即杆AB 位于初始位置时，如图所示，均为-10 rad/s 。

通过理论力学分析可知，在初始时刻，杆BC 和杆CD 的角速度大小等于杆AB ，方向为负，考虑到此结构速度不会突变，故在一段时间内，2θ与3θ将减小，比较图2.1与图4.1可验证这一点。

注意到在角位移350︒附近，两个杆的速度均有较快的变化，可见在此时机构有较大的冲击载荷作用于铰链上。

由于原动件的角速度不变为10 rad/s ，故图的X 横坐标除以10，就得到了角速度与时间的关系，考虑到速度的正负，可以看出2条曲线与横轴所围面积小于360，这意味着杆BC 和CD 均作摆动运动。

有平面四杆机构的基本知识[2]可知，该机构最短杆参与构成周转副，其余均为摆动副。

图4.1 1=0θ时，各杆的位置示意图图3.2显示了杆BC 和杆CD 在原动件杆AB 的带动下其各自角加速度随杆AB 的角位移的变化情况。

二者的变化趋势大致相同，在开始和结束存在较大的加速度变化，存在“凸点”。

中间的加速度变化较为平稳，可见机构运动较为平稳。

最终随着杆AB 完成一周的转动，二者的角加速度返回初始值。

综上，采用矩阵法分析计算得到的结果符合其他理论分析的结果和现实。

相较于其他方法，矩阵法便于计算机运算，得到的结果可以处理成图表，体现了矩阵法的优越性。

参考资料[1]李新，何传江. 矩阵理论及其应用. 第1版. 重庆：重庆大学出版社，2005. 96页[2]孙桓，陈作模，葛文杰. 机械原理. 第8版. 北京：高等教育出版社，2013. 131页附录MATLAB代码clcclearsyms x1 x y l1 l2 l3 l4 w1 cos2 cos3 theta m w2 w3;w=([-l2*sin(x),l3*sin(y);l2*cos(x),-l3*cos(y)])\(w1*[l1*sin(x1),-1*l1*cos(x1)]');alpha=([-1*l2*sin(x),l3*sin(y);l2*cos(x),-1*l3*cos(y)])\([w1^2*l1*cos(x1);w1^2*l1*sin(x1)]-[-1*w2*l2*cos(x),w3*l3*cos(y);-1*w2*l2*sin(x),w3*l3*sin(y)]*[w2;w3]);theta1=[0:1:360];n=1;while n<=361[c2,c3]=solve(0.12*cos2-0.09*cos3-0.12+0.06*cos(theta1(n)/180*pi),0.12*sqrt(1-cos2^2)-0.09* sqrt(1-cos3^2)+1*0.06*sin(theta1(n)/180*pi),cos2,cos3);coss2(n)=double(c2);coss3(n)=double(c3);n=n+1;endcoss2;coss3;a2=acos(coss2)/pi*180;a3=acos(coss3)/pi*180;omega=double(subs(w,{x,y,w1,x1,l1,l2,l3,l4},{acos(coss2),acos(coss3),10,theta1*pi/180,0.06,0. 12,0.09,0.12}))alph=double(subs(alpha,{x,y,w1,x1,l1,l2,l3,l4,w2,w3},{acos(coss2),acos(coss3),10,theta1*pi/18 0,0.06,0.12,0.09,0.12,omega(1,:),omega(2,:)}))figure(1)plot(theta1,omega(1,:),theta1,omega(2,:))xlabel('角位移(°)','FontSize',40)ylabel('角速度(rad/s)','FontSize',40)text(theta1(25),omega(1,25),' \leftarrow 杆BC','FontSize',40)text(theta1(25),omega(2,25),' \leftarrow 杆CD','FontSize',40)grid onhold onfigure(2)plot(theta1,alph(1,:),theta1,alph(2,:))xlabel('角位移(°)','FontSize',40)ylabel('角加速度(rad/s^2)','FontSize',40)text(theta1(25),alph(1,25),' \leftarrow 杆BC','FontSize',40) text(theta1(25),alph(2,25),' \leftarrow 杆CD','FontSize',40) grid onhold on。