牛头刨床运动分析 机械原理

——10A

班级：机械113

姓名：姚小龙

学号：201106263

位置方程

利用两个封闭图形ABDEA和EDCGE，建立两个封闭矢量方程，由此可得：

 (1)

把(1)式分别向轴、轴投影得：

 (2)

在（2）式中包含、、、四个未知数，消去其中三个可得到只含一个未知数

方程：

  (3)

当取不同值时，用牛顿迭代法解(3)式，可以求出每个的值，再根据方程组(2)可以求出其他杆件的位置参数、、的值：

                        (4)

速度方程

对(2)式对时间求一次导数并把结果写成矩阵的形式得：

 (5)

其中为刨刀的水平速度，为滑块2相对于杆3的速度。由于每个对应的、、已求出，方程组式（5）的系数矩阵均为常数，采用按列选主元的高斯消去法可求解（式5）可解得角速度ω3、ω4、、

加速度方程

把（5）对时间求导得矩阵式：

（6）

同样采用按列选主元的高斯消去法可求解（6）可得角加速度、、、

%主程序开始

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);                         %把这次计算的解作为下一次计算的初值。

end

for i=1:36 

%用循环结构求当theta1的值取不同值时，theta3、s3、s5的取值。因为theta3的值可能的取值范围为[0,π]，对theta3求解时应分以下两种情况讨论

if sign(h2+l1*cos(theta1(i))-l4*cos(theta4(i)))>0 %theta3<π/2

    theta3(i)=asin((h-l4*sin(theta4(i)))/l3);

else theta3(i)=pi-asin((h-l4*sin(theta4(i)))/l3); %theta3>π/2

    end

test(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));   

end

for 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); 

    %把求得的结构赋值给各物理量

end

for 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('曲柄转角 \heta_1/\\circ');

ylabel('角位移/\\circ');

gridon;hold on;

text(200,110,'\heta3');

%l4角位移图

subplot(2,2,2);

plot(i,theta4*du,'r');

title('角位移图');

xlabel('曲柄转角 \heta_1/\\circ');

ylabel('角位移/\\circ');

gridon;hold on;

text(150,10,'\heta4');

%滑块2位移

subplot(2,2,3);

plot(i,s3,'r');

title('位置');

xlabel('滑块位置 \\s3/\\circ');

ylabel('毫米/\\circ');

gridon;hold on;

text(150,500,'s3');

%c点位移

subplot(2,2,4);

plot(i,s5,'r');

title('位置');

xlabel('滑块位置 \\s3/\\circ');

ylabel('毫米/\\circ');

gridon;hold on;

text(150,0,'s5');

figure(2);

%l3角速度

subplot(2,2,1);

plot(i,omega3,'r');

title('角速度图');

xlabel('曲柄转角 \heta_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('曲柄转角 \heta_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('曲柄转角 \heta_1/\\circ');

ylabel('速度 mm/s');

grid on;

hold on;

text(200,0,'Vc');

%l4的速度

subplot(2,2,4);

plot(i,vBe,'r');

title('速度图');

xlabel('曲柄转角 \heta_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('曲柄转角 \heta_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('曲柄转角 \heta_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('曲柄转角 \heta_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('曲柄转角 \heta_1/\\circ')

ylabel('加速度/m\\cdots^{-2}')

grid on;

hold on;

text(200,100,'Abe');

%牛顿迭代法的函数定义

function [r,ea,iter]=NEWTON(fun,x,x0,es,maxit)

%定义函数名和输入输出的参数。输出参数为r,ea,iter。其中，r代表方程的解，ea代表最终解r代入方程的误差值，iter代表在运算过程中迭代的次数。fun,x,x0,es,maxit为输入参数。其中，fun代表要求解的方程，x代表要求解的未知数名称，x0代表求解过程取的初值，es表示求解要求的精度，maxit表示最大迭代步数。

if nargin<3,error('请输入包括函数名，函数变量，变量初始值在内的至少三个参数'),end

ifnargin<4||isempty(es),es=0.005;end

ifnargin<5||isempty(maxit),maxit=300;end

%对调用函数时输入的参数进行检查，如果调用时输入的参数不足3个，则报错，如果用户没输入es和maxit的值，则设置默认的求解精度为0.005，默认的迭代步数为300。

iter=0;                                                     %初始迭代次数为0。

t=sym('t','real');                      %定义符号变量，来取代输入函数的变量。

y=sym('y','real');                            %定义符号变量，来取代输入函数。

y=subs(fun,x,t);

y1=inline(y);                  %定义内联函数，表达式与y的表达式一致，方便计算。

y2=inline(t-y/diff(y,t));                     %定义内联函数，用于牛顿迭代。

r=x0;

for i=0:maxit

%定义一个循环结构来产生迭代过程，当循环次数大于最大迭代次数时，循环结束。

    r=y2(r);                                                        %牛顿迭代

iter=iter+1;                              %每循环一次，迭代次数iter加一。

if abs(y1(r))<=es||iter>=maxit,break;end

    %如果误差小于允许误差或循环次数大于最大迭代次数，迭代停止。

end

ea=y1(r);                                           %把最终迭代误差赋值给ea。

end

%NEWTON函数定义完成。

%按列选主元的高斯消去法的函数定义

function x = gauss( A,b )

%定义函数名和输入输出的参数。输出参数为列向量x，即线性方程的解向量。A,b为输入参数。其中，A为与解向量x维数相同的方阵，b为与解向量维数相同的列向量，这个函数的作用是解线性方程组‘AX=b’。

[m,n]=size(A); 

%获得A的行数和列数，其中m代表矩阵的行数，n代表矩阵的列数。

if m~=n, error('A必须是方阵'); end                        %检查A是否是方阵。

B=[A,b];                                       %把方阵A和向量b组成增广矩阵B。

for k=1:n-1                                        %用嵌套的循环结构进行消元。

[big,i]=max(abs(B(k:n,k)));

%找出B(k，k)、B（k+1，k）、B（k+2，k）…B(n,k)中的最大值。

    u=i+k-1;                                  %u为第k列中的最大元素所在的列。

if u~=k

       a=B(k,:);B(k,:)=B(u,:);B(u,:)=a;

       %把第u行的元素与第k行的元素位置互换。

end

for i=k+1:n

factor=B(i,k)/B(k,k);

        B(i,k:n+1)=B(i,k:n+1)-factor*B(k,k:n+1);

        %用初等行变换对第k+1行至第n行进行消元。

end

end

x=zeros(n,1);

x(n)=B(n,n+1);

for i=n-1:-1:1

x(i)=(B(i,n+1)-B(i,i+1:n)*x(i+1:n))/B(i,i);

%的回代过程

end

end

%gauss函数定义完成。

位置-时间曲线

速度-时间曲线

加速度-时间曲线