【保姆级教学】手把手教你完成机械设计基础大作业，运用解析法设计滚子直动凸轮，并校验最大压力角

核心代码，记得看注释，哪些地方要改注释已经标明

clc;clear;close all;
%%
pic_num=1;
%滚子直动从动件盘形凸轮
phi=120;
%回程运动角
phis=60;
%远休止角
phi_=120;
%推程运动角
phis_=60;
%近休止角
h=14;
%升程
r0=33;
%基圆半径
omega=0.001;
%角速度
r1=12;
%滚子半径
e=-12;
%偏置，正值为右偏置

%从角度转化为弧度
phi=phi*pi/180;
phis=phis*pi/180;
phi_=phi_*pi/180;
phis_=phis_*pi/180;

s0=(r0^2-e^2)^(1/2);
beta0=atan(e/s0);
beta=0;
i=1;thetaR=0;s=0;
while(thetaR<2*pi)
    %推程运动
if thetaR-beta+beta0<=phi
    s=spatial(thetaR-beta+beta0,h,phi,s);
    v=velocity(omega,h,phi,thetaR-beta+beta0,s);
    beta=atan(e/(s0+s));
    p=sqrt((s+s0)^2+e^2);
    alpha(i)=atan(abs(v/omega+e)/(s+s0));
    theta(i)=thetaR;
    randa=alpha(i)+beta;
    pT=sqrt(p^2+r1^2-2*p*r1*cos(randa));
    delta=atan(r1*sin(randa)/(p-r1*cos(randa)));
    thetaT=thetaR+delta;
    x(i)=pT*cos(thetaT+beta0);
    y(i)=pT*sin(thetaT+beta0);
    xd(i)=(s+s0)*cos(thetaR-beta+beta0)-e*sin(thetaR-beta+beta0);
    yd(i)=(s+s0)*sin(thetaR-beta+beta0)+e*cos(thetaR-beta+beta0);
    dot(i,:)=[x(i),y(i)];
    
    i=i+1;
end
    %远休止
if thetaR-beta+beta0>=phi&&thetaR-beta+beta0<=phis+phi
    s=14;
    v=0;
    beta=atan(e/(s0+s));
    p=sqrt((s+s0)^2+e^2);
    alpha(i)=atan(abs(v/omega+e)/(s+s0));
    theta(i)=thetaR;
    randa=alpha(i)+beta;
    pT=sqrt(p^2+r1^2-2*p*r1*cos(randa));
    delta=atan(r1*sin(randa)/(p-r1*cos(randa)));
    thetaT=thetaR+delta;
    x(i)=pT*cos(thetaT+beta0);
    y(i)=pT*sin(thetaT+beta0);
    xd(i)=(s+s0)*cos(thetaR-beta+beta0)-e*sin(thetaR-beta+beta0);
    yd(i)=(s+s0)*sin(thetaR-beta+beta0)+e*cos(thetaR-beta+beta0);
    dot(i,:)=[x(i),y(i)];
    i=i+1;
end
    %回程运动
    if thetaR-beta+beta0>=phis+phi&&thetaR-beta+beta0<=phis+phi+phi_
           s=spatial_(h,thetaR-beta+beta0,phis,phi_,phi);
    v=velocity_(h,omega,phi_,thetaR-beta+beta0,phi,phis);
    beta=atan(e/(s0+s));
    p=sqrt((s+s0)^2+e^2);
    alpha(i)=atan(abs(v/omega+e)/(s+s0));
    theta(i)=thetaR;
    randa=alpha(i)+beta;
    pT=sqrt(p^2+r1^2-2*p*r1*cos(randa));
    delta=atan(r1*sin(randa)/(p-r1*cos(randa)));
    thetaT=thetaR+delta;
    x(i)=pT*cos(thetaT+beta0);
    y(i)=pT*sin(thetaT+beta0);
    xd(i)=(s+s0)*cos(thetaR-beta+beta0)-e*sin(thetaR-beta+beta0);
    yd(i)=(s+s0)*sin(thetaR-beta+beta0)+e*cos(thetaR-beta+beta0);
    dot(i,:)=[x(i),y(i)];
    i=i+1;
    end
    %近休止
    if thetaR-beta+beta0>=phis+phi+phi_&&thetaR-beta+beta0<=phis+phi+phi_+phis_
    s=0;
    v=0;
    beta=atan(e/(s0+s));
    p=sqrt((s+s0)^2+e^2);
    alpha(i)=atan(abs(v/omega+e)/(s+s0));
    theta(i)=thetaR;
    randa=alpha(i)+beta;
    pT=sqrt(p^2+r1^2-2*p*r1*cos(randa));
    delta=atan(r1*sin(randa)/(p-r1*cos(randa)));
    thetaT=thetaR+delta;
    x(i)=pT*cos(thetaT+beta0);
    y(i)=pT*sin(thetaT+beta0);
    xd(i)=(s+s0)*cos(thetaR-beta+beta0)-e*sin(thetaR-beta+beta0);
    yd(i)=(s+s0)*sin(thetaR-beta+beta0)+e*cos(thetaR-beta+beta0);
    dot(i,:)=[x(i),y(i)];
    i=i+1;
    end
    thetaR=thetaR+omega;
end

[max,position]=max(alpha);   %输出最大压力角及其所在位置
disp(max*180/pi);
x(i)=x(2);
y(i)=y(2);
dot(i,:)=[x(i),y(i)];
disp(theta(position)*180/pi);     %输出最大压力角对应的转角
min=10000;
for j=2:i-2
    d(j)=get_r(dot(j,:),dot(j+1,:),dot(j+2,:));
    if d(j)>15  &&d(j)<min
        min=d(j);
    end
end
disp(min);      %输出最小曲率半径
plot(x,y)%绘图

辅助的四个函数，以下四个函数为简谐运动，其他题目要求的等加速或者正弦加速公式请根据高等教育出版社的《机械设计基础》第七版第40,41页来修改这四个函数

一、推程位移函数，distance为化为弧度制的推程转角，height为升程，phine为phine转角，distance_为化为弧度制的回程转角，omega为角速度，distance_s为化为弧度制的远休止角转角

function output = spatial(phine,height,distance)
%UNTITLED 此处显示有关此函数的摘要
%   此处显示详细说明
output=height*(1-cos(phine*pi/distance))/2;
end

二、推程速度函数，

function output = velocity(omega,height,distance,phine)
output=height*pi*omega*sin(pi*phine/distance)/(2*distance);
end

三、回程位移函数

function output = spatial_(height,phine,distance_s,distance_,distance)
%UNTITLED 此处显示有关此函数的摘要
%   此处显示详细说明
output=height*(1+cos(pi/distance_*(phine-distance-distance_s)))/2;
end

四、回程速度函数

function output = velocity_(height,omega,distance_,phine,distance,distance_s)
%VELOCITY_ 此处显示有关此函数的摘要
%   此处显示详细说明
output=-height*pi*omega/(2*distance_)*sin(pi/distance_*(phine-distance-distance_s));
end

获得曲率半径

function r = get_r(p1,p2,p3)
%GET_R 此处显示有关此函数的摘要
%   此处显示详细说明
a=[p1(1),p1(2),1;p2(1),p2(2),1;p3(1),p3(2),1];
b=[-p1(1)^2-p1(2)^2,-p2(1)^2-p2(2)^2,-p3(1)^2-p3(2)^2]';
x=a\b;
D=x(1);E=x(2);F=x(3);
r=sqrt(D^2+E^2-4*F)/2;
end