matlab六杆机构运动分析

以图示六杆机构为例，已知构件1的运动，确定机构中其它构件的运动（包括位移、速度、加速度）

参考书籍：《机械原理matlab辅助设计》

一.建立数学模型

对六杆机构进行运动分析时，可以将其拆分成两个四杆机构，采用封闭矢量多边形法求解。首先建立机构封闭位置矢量方程式，之后对位置方程求一次导和二次导得到其速度和加速度方程，求解即可。

二.程序设计框图

每个平面连杆机构运动分析MATLAB程序都由主程序和子函数两部分组成，其程序设计流程如下图。

三.确立数学解析式

为了对机构进行运动分析，以A为原点建立直角坐标系，并将各构件表示为杆矢。假设x轴正方向为水平向右，分别代表杆1，2，4，5与X轴正方向的距离（变量命名详见函数注释）。

1.位置分析

如图所示，由封闭图形ABCA可写出机构各杆矢所构成的封闭矢量方程

将矢量方程分别投影到x轴和y轴可得

同理写出封闭图形CDEC各杆矢所构成的封闭矢量方程，联立可得六杆机构位置方程。

2.速度分析

将位置矢量方程对时间t求一次导数，得速度关系

若用矩阵形式来表示，则上式可写为

求解可得角速度w2和滑块速度v3。

3. 加速度分析

将式（1-2）对时间t求二次导数，可得加速度关系表达式

同理，加速度矩阵方程如下

求解可得角加速度a2和滑块加速度a3。

由封闭图形CDEC可写出机构各杆矢所构成的封闭矢量方程，位置分析、速度分析与加速度分析同理；联立两个四杆机构分析方程即可解出机构中每一构件的运动。

四.matlab程序分析

六杆机构的图形如上图片，假设杆1的长度为101.6mm，杆2的长度为254mm，杆3的长度为177.8mm，杆4的长度为304.8mm，原动件1以匀角速度w1 = 250 rad/s 逆时针转动，进行分析：

% 调用fsolve函数求解位移方程
function f=pos_6(theta,theta1,length1,length2,length4,length5,L1,L2)
f(1)=length1*cos(theta1)+length2*cos(theta(1))-theta(4);
f(2)=length1*sin(theta1)+length2*sin(theta(1));
f(3)=length4*cos(theta(2))+length5*cos(theta(3))-L1+theta(4);
f(4)=length5*sin(theta(3))+length4*sin(theta(2))-L2;
end
% theta1、theta2、theta4、theta5分别代表杆1，2，4，5与X轴正方向的距离
% omega2、omega4、omega5代表杆2，4，5的角速度
% length3表示滑块3和A点的距离
% velocity表示滑块3的线速度
% alpha_line表示滑块3的加速度
% theta2=theta(1);
% theta4=theta(2);
% theta5=theta(3);
% length3=theta(4);

这里我采用的是matlab用来求解非线性方程组的fsolve函数，求解原理就是利用最小二乘法解方程组。其中fsolve返回值为函数的解，pos_6是用于定义需求解的非线性方程组的函数文件名，X0是求根过程的初值，options为最优化工具箱的选项设定。Display选项决定函数调用时中间结果的显示方式，其中off为不显示，optimset(‘Display’,‘off’)将设定Display选项为off其中初值X0的取值。

function [theta,omega,alpha] = analysis_data_6(theta1,omega1,alpha1,length1,length2,length4,length5,L1,L2)
options=optimset('display','off');
x0=[-1,0,1,100];
theta=fsolve(@pos_6,x0,options,theta1,length1,length2,length4,length5,L1,L2);
% theta=theta.';   % 矩阵转置
theta2=theta(1);
theta4=theta(2);
theta5=theta(3);
length3=theta(4);% 计算角速度 omega
A=[-length2*cos(theta2),0,0,0;0,(length4)*sin(theta4),(length5)*sin(theta5),-1;0,length4*cos(theta4),length5*cos(theta5),0;-(length2)*sin(theta2),0,0,-1];
B=[length1*cos(theta1);0;0;(length1)*sin(theta1)];
omega = A\(B*omega1);
omega2 = omega(1);
omega4 = omega(2);
omega5 = omega(3);
velocity = omega(4);   % 滑块3的线速度% 计算角加速度alpha
A1 = [length2*cos(theta2),0,0,0;0,length4*sin(theta4),length5*sin(theta5),-1;0,length4*cos(theta4),length5*cos(theta5),0;-length2*sin(theta2),0,0,-1];
At = [-omega2*length2*sin(theta2),0,0,0;0,omega4*length4*cos(theta4),omega5*length5*cos(theta5),0;0,-omega4*length4*sin(theta4),-omega5*length5*sin(theta5),0;-omega2*length2*cos(theta2),0,0,0];
Bt = [omega1*length1*sin(theta1);0;0;omega1*length1*cos(theta1)];
alpha = A1\(-At*omega+omega1*Bt);
end

主程序如下

length1=0.025;
length2=0.071;
length4=0.038;
length5=0.030;
L1=0.1;
L2=0.04;
omega1 = 10;
alpha1 = 0;
hd = pi/180;
du = 180/pi;%调用函数计算六杆机构的位移、速度、加速度
for n1 = 1:361theta1 = (n1-1)*hd;        % 将角度转换为弧度制[theta,omega,alpha] = analysis_data_6(theta1,omega1,alpha1,length1,length2,length4,length5,L1,L2);theta2(n1) = theta(1);        % 杆2的角位移theta4(n1) = theta(2);        % 杆4的角位移theta5(n1) = theta(3);        % 杆5的角位移length3(n1) = theta(4);omega2(n1) = omega(1);        % 杆2的角速度omega4(n1) = omega(2);        % 杆4的角速度omega5(n1) = omega(3);        % 杆5的角位移velocity(n1) = omega(4);alpha2(n1) = alpha(1);    % 杆2的角加速度alpha4(n1) = alpha(2);    % 杆4的角加速度alpha5(n1) = alpha(3);  % 杆5的角加速度alpha_line(n1) = alpha(4);
endn1 = 1:361; % 绘制图像subplot(2,2,1);
plot(n1,theta2,n1,theta4,n1,theta5,'k');    % 杆2、杆4、杆5随角度变化
xlabel('曲柄转角 \theta1/\circ');
ylabel('角位移/\circ');
legend('theta2','theta4','theta5');
grid on;
hold on;
% 以构件5的分析为例
subplot(2,2,2);
% plot(n1,(theta5-theta5(1))*length5,'k');    % 杆5角位移随角度变化
plot(n1,theta5*du,'k');    % 杆5角位移随角度变化
xlabel('曲柄转角 \theta1/\circ');
ylabel('角位移/\circ');
legend('theta5');
grid on;
hold on;subplot(2,2,3);     % 杆5角速度变化
plot(n1,omega5,'k');
xlabel('曲柄转角 \theta1/\circ');
ylabel('角速度/ rad\cdots^{-1}');
legend('w5');
grid on;
hold on;subplot(2,2,4);     %杆5角加速度变化
plot(n1,alpha5,'k');
xlabel('曲柄转角 \theta1/\circ');
ylabel('角加速度/ m\cdots^{-1}');
legend('a5');
grid on;
hold on;

运行结果如下图

最后做了一个机构运动仿真模拟六杆机构的运动

j=0;
for n1= 1:5:360
clf;
j=j+1;
x(1)=0;
y(1)=0;
x(2)=length1* cos((n1 - 1)*hd);
y(2)=length1* sin((n1 - 1)*hd);
x(3)=length3(n1);
y(3)=0;
x(4)=length3(n1)+length4*cos(theta4(n1));
y(4)=length4*sin(theta4(n1));
x(5)=L1;
y(5)=L2;
x(6)=0;
y(6)=0;
plot(x,y);
grid on;hold on;
plot(x(1),y(1),'o');
plot(x(2),y(2),'o');
plot(x(3),y(3),'o');
plot(x(4),y(4),'o');
plot(x(5),y(5),'o');
plot(x(6),y(6),'o');
axis([-0.05 0.14,-0.05 0.14]);
xlabel('mm');
ylabel('mm');
m(j) = getframe;        %getframe函数的作用是捕获坐标区或图窗作为影片帧。
%     f(:,j)=getframe;
end
movie(m);

仿真结果就不上传了，主要是没整明白怎么上传动图