机械原理大作业——10A班级：姓名：学号：位置方程利用两个封闭图形ABDEA和EDCGE，建立两个封闭矢量方程，由此可得：⎭⎬⎫+=++=+' s l l s l l l l56431643 (1) 把(1)式分别向x 轴、y 轴投影得：⎪⎪⎭⎪⎪⎬⎫=+=++=++=+ h l l s l l l h s l l h s l 33445334411133441123344sin sin cos cos sin sin sin cos cos cos θθθθθθθθθθ(2) 在（2）式中包含3s 、5s 、3θ、4θ四个未知数，消去其中三个可得到只含4θ一个未知数 方程：[][]{}[][]sin sin sin 2sin cos cos sin sin 244111234424224244112244111=-+--+-++-+θθθθθθθθl l h l hl h l l l h l l h(3)当1θ取不同值时，用牛顿迭代法解(3)式，可以求出每个4θ的值，再根据方程组(2)可以求出其他杆件的位置参数3s 、5s 、3θ的值：⎪⎪⎭⎪⎪⎬⎫-+=+=-= 3441113334453443sin sin sin cos cos )sin arcsin(θθθθθθθl l h s l l s l l h (4)速度方程对(2)式对时间求一次导数并把结果写成矩阵的形式得：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-----00cos sin 0cos cos 01sin sin 00cos cos sin 0sin sin cos 1111143443344334433344333θθωωωθθθθθθθθθθl l v v l l l l l s l s C e B (5)其中C v 为刨刀的水平速度，v eB 为滑块2相对于杆3的速度。

由于每个1θ对应的3s 、3θ、4θ已求出，方程组式（5）的系数矩阵均为常数，采用按列选主元的高斯消去法可求解（式5）可解得角速度ω3、ω4、eB v 、C v加速度方程把（5）对时间求导得矩阵式：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--+⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡----------=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-----00sin cos 0sin sin 00cos cos 00sin sin cos cos 0cos cos sin sin 0cos cos 01sin sin 00cos cos sin 0sin sin cos 11111114344433344433344433333344433333343443344334433344333θωθωωωωθωθωθωθωθωθωθθωθωθωθθωααl l v v l l l l l s v l s v a a θl θl θl θl θl θs θθl θs θC eBe B eBC e B（6）同样采用按列选主元的高斯消去法可求解（6）可得角加速度3α、4α、eB a 、C a %主程序开始clear;clc;l1=180; %L1=lab l3=960; %l3=lCD l4=160; %l4=lED h=900; h1=460; h2=110; du=180/pi; omega1=1; alpha1=0;theta1=linspace(0,35*pi/18,36);%定义常量和已知参数，l1代表杆1的长度，l3代表杆3的长度，l4代表杆4的长度,h 表示EG 的长度，h1表示AE 的竖直距离,h2表示AE 的水平距离，theta1表示角θ1的不同值。

theta3=zeros(1,36);theta4=zeros(1,36);s3=zeros(1,36);s5=zeros(1,36);test=zeros(1,36);vBe=zeros(1,36);vc=zeros(1,36);omega1=ones(1,36);omega3=zeros(1,36);omega4=zeros(1,36);aBe=zeros(1,36);ac=zeros(1,36);alpha1=zeros(1,36);alpha3=zeros(1,36);alpha4=zeros(1,36);A=zeros(4,4); dA=zeros(4,1);%定义最终的结果数据，当θ1取不同值时,theta3表示θ3的值，theta4表示θ4的值，s3表示BD 的长度，s5表示GC 的长度，vBe 表示B 点在杆3上运动的速度，vc 表示杆5的运动速度，即牛头刨刀的速度，omega3表示杆3的转动角速度，omega4表示杆4的转动角速度，aBe 表示B 点在杆3上运动的角加速度，ac 表示杆5的加速度，即牛头刨刀的加速度，alpha3表示杆3的角加速度，alpha4表示杆4的角加速度，矩阵A,dA 表示线性方程组的系数矩阵i=0; %i 为循环变量，在循环结构中使用 syms THETA1 THETA4 %定义符号变量，为以下计算做准备 fun1=((h1+l1*sin(THETA1)-l4*sin(THETA4))^2+(h2+l1*cos(THETA1)-l4*cos(THETA4))^2)*(l4^2*sin(THETA4)^2+h^2-2*h*l4*sin(THETA4))-l3^2*(h1+l1*sin(THETA1)-l4*sin(THETA4))^2; %定义迭代法中要求解的关于THETA4的方程。

x0=0; %定义在牛顿迭代法中的变量THEA4的初值。

for i=1:36 %用循环结构求当theta1取不同值时，theta3值。

fun2=subs(fun1,THETA1,theta1(i));%把不同的THETA1的值代入要求解的方程[theta4(i),EA,it]=NEWTON(fun2,'THETA4',x0,0.0001,1000);%用牛顿迭代法求得THEATA4，并赋值到theta4的数组中x0=theta4(i); %把这次计算的解作为下一次计算的初值。

endfor i=1:36%用循环结构求当theta1的值取不同值时，theta3、s3、s5的取值。

因为theta3的值可能的取值范围为[0,π]，对theta3求解时应分以下两种情况讨论if sign(h2+l1*cos(theta1(i))-l4*cos(theta4(i)))>0 %theta3<π/2theta3(i)=asin((h-l4*sin(theta4(i)))/l3);else theta3(i)=pi-asin((h-l4*sin(theta4(i)))/l3); %theta3>π/2endtest(i)=h1+l1*sin(theta1(i))-l4*sin(theta4(i));s5(i)=l4*cos(theta4(i))+l3*cos(theta3(i));s3(i)=(h1+l1*sin(theta1(i))-l4*sin(theta4(i)))/sin(theta3(i));endfor i=1:36%用循环结构求当theta1的值取不同值时vBe、omega3、omega4、vc的值。

A(1,1)=cos(theta3(i));A(1,2)=-s3(i)*sin(theta3(i));A(1,3)=-l4*sin(theta4(i));A(2,1)=sin(theta3(i));A(2,2)=s3(i)*cos(theta3(i));A(2,3)=l4*cos(theta4(i));A(3,2)=-l3*sin(theta3(i));A(3,3)=-l4*sin(theta4(i));A(3,4)=-1;A(4,2)=l3*cos(theta3(i));A(4,3)=l4*cos(theta4(i));dA(1,1)=-omega1(1,1)*l1*sin(theta1(i));dA(2,1)=omega1(1,2)*l1*cos(theta1(i));x=gauss(A,dA); %用按列选主元的高斯消去法求解vBe(i)=x(1);omega3(i)=x(2);omega4(i)=x(3);vc(i)=x(4);%把求得的结构赋值给各物理量endfor i=1:36 %用循环结构求当theta1的值取不同值时aBe、alpha3、alpha4、vc的值。

A(1,1)=cos(theta3(i));A(1,2)=-s3(i)*sin(theta3(i));A(1,3)=-l4*sin(theta4(i));A(2,1)=sin(theta3(i));A(2,2)=s3(i)*cos(theta3(i));A(2,3)=l4*cos(theta4(i));A(3,2)=-l3*sin(theta3(i));A(3,3)=-l4*sin(theta4(i));A(3,4)=-1;A(4,2)=l3*cos(theta3(i));A(4,3)=l4*cos(theta4(i));dA(1,1)=-omega3(i)*sin(theta3(i))*vBe(i)*2-s3(i)*omega3(i)^2*cos(theta3(i))-l4*omega4(i)^ 2*cos(theta4(i))-l1*cos(theta1(i));dA(2,1)=omega3(i)*cos(theta3(i))*vBe(i)*2-s3(i)*omega3(i)^2*sin(theta3(i))-l4*omega4(i)^2*sin( theta4(i))-l1*sin(theta1(i));dA(3,1)=-l3*omega3(i)^2*cos(theta3(i))-l4*omega4(i)^2*cos(theta4(i));dA(4,1)=-l3*omega3(i)^2*sin(theta3(i))-l4*omega4(i)^2*sin(theta4(i));%构造速度方程的系数矩阵x=gauss(A,dA); %用按列选主元的高斯消去法求解aBe(i)=x(1);alpha3(i)=x(2);alpha4(i)=x(3);ac(i)=x(4);%把求得的结构赋值给各物理量end%主程序结束%出图程序figure(1);i=1:10:360;%l3角位移图subplot(2,2,1);plot(i,theta3*du,'r');title('角位移图');xlabel('曲柄转角\theta_1/\circ'); ylabel('角位移/\circ');grid on;hold on;text(200,110,'\theta3');%l4角位移图subplot(2,2,2);plot(i,theta4*du,'r');title('角位移图');xlabel('曲柄转角\theta_1/\circ'); ylabel('角位移/\circ');grid on;hold on;text(150,10,'\theta4');%滑块2位移subplot(2,2,3);plot(i,s3,'r');title('位置');xlabel('滑块位置\s3/\circ');ylabel('毫米/\circ');grid on;hold on;text(150,500,'s3');%c点位移subplot(2,2,4);plot(i,s5,'r');title('位置');xlabel('滑块位置\s3/\circ');ylabel('毫米/\circ');grid on;hold on;text(150,0,'s5');figure(2);%l3角速度subplot(2,2,1);plot(i,omega3,'r');title('角速度图');xlabel('曲柄转角\theta_1/\circ') ylabel('角速度/rad\cdots^{-1}') grid on;hold on;text(150,0,'\omega_3');%l4角速度subplot(2,2,2);plot(i,omega4,'r');title('角速度图');xlabel('曲柄转角\theta_1/\circ') ylabel('角速度/rad\cdots^{-1}') grid on;hold on;text(150,0.2,'\omega_4');%c点速度subplot(2,2,3);plot(i,vc,'r');title('速度图');xlabel('曲柄转角\theta_1/\circ'); ylabel('速度mm/s');grid on;hold on;text(200,0,'Vc');%l4的速度subplot(2,2,4);plot(i,vBe,'r');title('速度图');xlabel('曲柄转角\theta_1/\circ'); ylabel('速度mm/s');grid on;hold on;text(200,-100,'Vbe');figure(3);%l3 角加速度图subplot(2,2,1);plot(i,alpha3,'r');title('角加速度图');xlabel('曲柄转角\theta_1/\circ') ylabel('角加速度/rad\cdots^{-2}')grid on;hold on;text(200,-0.1,'\alpha_3');%l4 角加速度图subplot(2,2,2);plot(i,alpha4,'r');title('角加速度图');xlabel('曲柄转角\theta_1/\circ')ylabel('角加速度/rad\cdots^{-2}')grid on;hold on;text(200,-0.5,'\alpha_4');%c点加速度图subplot(2,2,3);plot(i,ac,'r');title('加速度图');xlabel('曲柄转角\theta_1/\circ')ylabel('加速度/m\cdots^{-2}')grid on;hold on;text(200,150,'ac');%c点加速度图subplot(2,2,4);plot(i,aBe,'r');title('加速度图');xlabel('曲柄转角\theta_1/\circ')ylabel('加速度/m\cdots^{-2}')grid on;hold on;text(200,100,'Abe');%牛顿迭代法的函数定义function [r,ea,iter]=NEWTON(fun,x,x0,es,maxit)%定义函数名和输入输出的参数。