直动式固定凸轮与连杆机构的设计-毕业设计.doc

湖南工业大学
直动式固定凸轮与连杆机构的设计

设计者
所在院(系)：湖南工业大学
专业：机械设计制造及其自动化
班级
学号：
指导老师：
时间：2015年12月27日

目录
一、课程设计的目的 1
二、设计内容与步骤 1
1、设计内容 1
2.设计步骤 2
三、设计要求 2
四、设计指导 3
1、概述 3
2、基本参数 6
3、设计步聚 7
1）确定驱动方案 7
2）确定e 8
3）确定h 8
4）确定 8
5）确定 8
6） 求算b1、b2 9
7）设计凸轮廊线 10
8）检验压力角 11
五、参数优化 13
六、结论 14
七、参考文献 14
八、附图 14

摘要

包装设计课程设计是在完成机械设计课程学习后，一次重要的实践性教学环节。是高等工科院校大多数专业学生第一次较全面的设计能力训练，也是对机械设计课程的全面复习和实践。其目的是培养理论联系实际的设计思想，训练综合运用机械设计和有关选修课程的理论，结合生产实际分析和解决工程实际问题的能力，巩固、加深和扩展有关机械设计方面的知识。
本次设计的题目是直动式固定凸轮与连杆机构的设计。根据题目要求和机械设计的特点作者做了以下几个方面的工作：①根据有关参数进行计算或编写有关设计计算程序；②利用程序设计的方法输出结果并自动生成图形；③画出装配图及其主要零件图；④完成设计计算说明书。

一、 课程设计的目的

《包装机械设计》课程设计是本课程各教学环节中重要的一环,它让学习者联系实际进一步深入理解、掌握所学的理论知识。其基本目的是：培养理论联系实际的设计思想，训练综合运用包装机械和有关先修课程的理论，结合生产实际分析和解决工程实际问题的能力，巩固、加深和扩展有关包装机械设计方面的知识。通过制订设计方案，合理选择裹包机中块状物品推送机构和零件类型，正确计算零件工作能力、确定尺寸和选择材料，以及较全面地考虑制造工艺、使用和维护等要求，之后进行结构设计，达到了解和掌握机械零件、包装机械经常采用的机构的设计过程和方法。进行设计基本技能的训练。例如计算、绘图、熟悉和运用设计资料（手册、图册、标准和规范等）以及使用经验数据、进行经验估算和处理数据的能力。

二、设计内容与步骤
1、设计内容
　　以裹包机中块状物品推送机构的典型机构——固定凸轮与连杆组合机构为题。课程设计通常包括如下内容：读懂块状物品推送机构典型机构——固定凸轮与连杆组合机构，了解设计题目要求；分析该块状物品推送机构设计的可能方案；具体计算和设计该方案中机构的基本参数；进行机体结构及其附件的设计；绘制装配图及零件工作图；编写计算说明书以及进行设计答辩。
2.设计步骤
（1）设计准备
认真研究设计任务书，明确设计要求、条件、内容和步骤；通过阅读有关资料、图纸、参观实物或模型、观看电视教学片、挂图以及推送机构进行拆装实验等，了解设计对象；复习有关课程内容，熟悉零部件的设计方法和步骤；准备好设计需要的图书、资料和用具；拟定设计计划等。
（2）推送机构装置的总体设计
决定推送机构装置的方案；选择机构的类型，计算机构装置的运动参数。
（3）装配图设计
计算和选择机构的参数；确定机体结构和有关尺寸；绘制装配图草图；选择计算轴承和进行支承结构设计；进行机体结构及其附件的设计；完成装配图的其他要求；审核图纸。
（4）零件工作图设计
（5）整理和编写计算说明书
（6）设计总结和答辩

三、设计要求
在课程设计之前，准备好必要的设计手册或参考资料，以便在设计过程中逐步去学习查阅资料。确定设计题目后，至少应复习在课程中学过的相关内容。完成本课程设计的具体要求如下：
1、设计说明书要全面反映设计思想、设计过程和结论性认识。其工艺设计要有文字、计算、公式来源、参数选取的资料名称或代号、图表（草图）。说明书用A4纸打印，约20页左右，并装订成册。
2、设计图样按“机械制图”、“公差与配合”等国家标准完成。
3、零件图按生产图样要求完成，零件的有关精度和技术要求要有合理的标注或说明。
　　设计过程中，提倡独立思考、深入钻研，主动地、创造性地进行设计，反对不求甚解、照抄照搬或依赖老师。要求设计态度严肃认真、有错必改，反对敷衍塞责，容忍错误的存在。只有这样，才能保证课程设计达到教学基本要求，在设计思想、设计方法和设计技能等方面得到良好的训练。

四、设计指导
裹包机所包装的产品，绝大多数是单件或多件集合而成的块状物品。包装作业线中前后机之间物品的输送、换向、排列组合，及单机内部的物品移动等，需要用各种各样的机构或装置完成。以下是几种典型的推送块状物品的组合机构一一固定凸轮与连杆组合机构。
1、概述
图1所示，是该机构的结构简图，用于香皂、糖果等裹包机中，将物品向上推送较大距离。原动杆件AB按逆时针方向转动，驱动铰销C上的滚动轴承6在固定槽凸轮4的槽内运动，再通过连杆CD使推送杆（即滑块）2按预定规律作上下往复移动。这种直动从杆类型的固定凸轮与连杆组合机构相当于连杆长度可变的曲柄滑块机构，曲柄为AB，滑块为推送杆，连杆为BD，在运动过程中连杆BD的长度是变化的。

图1　直动从动杆类型的固定凸轮和连杆组合的推送机构结构简图
1－推送板　2－推送杆　3－导轨　4－固定槽凸轮　5－支座
6－滚动轴承　7－导轨

图2　固定凸轮与连杆组合机构示意图
1-推料板 2-推料杆 3-固定凸轮 4-滚子
以所示推送机构，还有2种运动方式:一是是曲柄AB推着杆件BC运动，杆件BC承受压力；二是曲柄AB拉着BC杆运动，杆件BC承受拉力。这是两种不同的驱动方案。当然，无论是前者还是后者，都可以在两种驱动方案中任意选择，本题我们选择直动从动杆类型的固定凸轮和连杆组合机构。
2、基本参数
为研究方便，特规定：以曲柄回转中心A为坐标的原点，并作x、y轴。对于直动从动类型（见图1所示），y轴与从动杆的运动方向平行；对于摆动从动杆类型（见图3所示）y轴与铰销D的两个运动极限位置之连线D0D1平行。考虑到曲柄有两种转向，又规定y轴的正轴逆着曲柄转向旋转900后所得轴为轴的正轴，于是，前者 轴的正轴向右，而手者则向左。基本参数有：
e—y轴与D0D1线的间距，简称偏心距；
h—铰销D至轴的最小距离；
a—曲柄AB长；
b1、b2—杆件BC、CD长；
—从动杆升程运动起始时刻的曲柄位置AB0和y轴负轴的夹角，=1800　－∠B0AY。
铰销B和D的距离用b表示，b=BD，它的最大值和最小值分别用bmax、bmin表示。
已知参数：
从动杆的行程

升程和降程的运动规律
100

升程和降程按五次多项式运动
表 1 初始参数表

固定凸轮与连杆组合机构的特点是，从动杆的运动可以象凸轮机构的从动杆那样实现停留和按照定规律（如五次多项式）运动。从动杆的行程、动停时间、运动速度由工艺要求预先给定。这样，当参数e、h、 、，确定后，每一运动时刻的b值及、值也随之确定。显然，b1、b2应满足下式
(1)
因此，应根据从动杆的运动规律和确定的c、h、、、l值，先计算出bmax、bmin，然后用下式求算b1、b2值：
(2)
3、设计步聚
1）确定驱动方案
它对凸轮的压力角机构的传动效率影响较大。应根据运动要求确定之。用下列符号表示运动要求：
Sm—分别为直动总行程；
1—升程运动对应的曲柄转角；
—最高位置停留对应的曲柄转角；
—降程运动对应的曲柄转角；
—最低位置停留对应的曲柄转角；
=3600，当时，先用曲柄AB拉着杆件BC运动的方案；当时，应选用曲柄AB推着BC杆运动的方案。时，可任选其中一种方案。
2）确定e
直动从动杆，取e=0
3）确定h
从结构紧凑和减小凸轮压力角考虑，应将h值取小些。但h值愈小，对从动杆驱动力的压力角也愈大。通常取h≥Sm,h=3sm。
4）确定
若a值过小，会使凸轮压力角明显增大，甚至不能实现预期动动。可取a=0.8Sm。
5）确定
其值对凸轮的压力角影响极大，过小，尤其是过大，会使压力角急剧增加。在前述参数确定后，最好将优化，目标函数为

式中a1m为凸轮的最大压力角。
δ=-4
6） 求算b1、b2
须先求算bmax、bmin。
（1）直动从动杆类型
参阅图4，依据铰销B、D的坐标，可建立它们之间距离的公式。B的坐标为
（3）

D的坐标为：
（4）
式中 —曲柄转角，取升程起始时的 =0°；S 为与相对应的从动杆位移，即铰销D至其最低位置的距离。S值分为升程（=0～1）、最高位置停留（=1～1+2）、降程（=1+2～1+2+3）、最低位置停留（=1+2+3～360°）四个阶段求算。b值为:
b= (5)
将式（3）、（4）代入式（5），求算bnax、bmin，然后用式（2）算得 b1、 b2。
因此,应根据从动杆的运动规律和确定的e , h , a 值,先计算出bmax , bmin ,利用MA TL AB 强大的可视化功能,绘出b 随φ的曲线图,见图4 ,程序如附录程序一.
可以很方便地得到: bmax = 397.8748mm , bmin =256.9556mm。

图4 　B D 长度曲线图
然后再利用（2）式求得、分别为：b1 =70.4613mm , b2 = 327.4171 mm
7）设计凸轮廊线
固定凸轮的理论廊线就是滚子中心C的运动轨迹线，根据铰销B、D的位置及b1、b2值可确定C的位置。令铰销B , D 的连线B D 与D0 D1 线(或y 轴) 的夹角为, BD与CD的夹角为β,则:
(6)
（7）
显然, xB > x D 时θ为正值,反之则为负值,而β始终为正值。铰销C 的坐标为:
（8）
该式运算符号“ + ”和“ —”的确定原则是:令b = bmax时的为 , b = bmin时的为,在～区间,取“ - ”号;在～和～区间,取“ + ”。由图4可得到φm =82.8496°(1.4460 弧度) ,φ′m= 216.7442°(3.7829弧度) 。
在MATLAB中，编程计算得到凸轮理论轮廓即C点的轨迹上各点的坐标，并绘出C点的轨迹线，见图5,程序如附录程序二。:

图5　凸轮轮廓线图

8）检验压力角
凸轮机构设计中,初始参数选择不当会导致凸轮轮廓设计不合理,从而导致从动件运动不能按预定规律进行,或者运行中动力学性能欠佳。
（1）凸轮的压力角
参阅图2，为PC和Vc的夹角。PC为驱动铰销运动的力(不考虑摩擦力)，BC重合，vc为铰链C的运动方向,，与C点的凸轮廊线切线重合。用K1、K2分别表示PC、vc的斜率，则
K1 = （9）
　 （10）
（11）
计算式K2 时需求函数或,的一阶导数,这些函数很难用显性方程表示,求解一阶导数更加繁琐,因此本文利用MA TL AB 方便的向量运算功能来求解K2 。将凸轮细分为若干个点,第i 点的斜率可采用式(12) 来求解。
(12)
在MATLAB 中对式(12) 进行求解非常方便,其程序代码简短易懂,如下:
For i=2:length(fa)-1 %fa凸轮某一段的转角
=((i-1))-(i+1)/((i-1))-(i+1)
end
应保证的最大值不超过许用值，即。可取利用附录程序三，通过matlab得到该机构的 = 44.1427°< 45°,设计合理。

（2）从动杆的压力角
为CD杆对铰销D的驱动力（不考虑摩擦力）与D的运动方向的夹角。对于直动从动杆，显然
　　 （13）
对于AB推动BC的驱动方案，在区间取“—”号，其它区间取“+”号。应使的最大值。直动从动杆。利用附录程序四，通过matlab计算得到该机构的
通过对压力角的检验可知,以上参数的设计合理。

五、参数优化
δ值是从动杆升程运动起始时刻的曲柄位置A B 0 和y 轴负轴的夹角,通过改变δ值,观察凸轮压力角和从动杆压力角的变化 ,见表2。

5
4
3
2
1
0
-1
-2
-3

61.30
58.67
56.31
54.18
52.22
50.41
48.72
47.12
45.60

25.35
25.7
25.79
26.02
26.24
26.45
26.67
26.88
27.09

-4
-5
-6
-7
-8
-9
-10
-11
-12

44.14
44.22
45.02
45.82
46.62
47.61
48.43
49.44
50.45

27.30
27.51
27.75
27.92
28.12
28.23
28.42
28.61
28.83

表 2 与压力角的关系
δ值对从动杆压力角影响很小, 但对凸轮压力角影响极大。当δ为5°时,凸轮压力角超出了许用范围;随着δ的减小, 变小;当δ为- 4°时, 达到最小值44.14°;随着δ继续减小, 增大。可见当以上其它参数不变的情况下,δ为- 4°,该机构的动力性能达到最优。

六、结论
1) 运用MA TL AB 强大的编程功能实现了凸轮连杆机构的计算机辅助设计计算,得到了连杆的长度和凸轮曲线的轨迹,方便地检验了凸轮压力角和从动杆压力角,并分析了关键参数δ值与压力角的关系,最后找到最优方案,从而提高了机构的动力性能。
2) 借用MA TL AB 简单的数组计算功能,提出了凸轮压力角有效的计算方法。

七、参考文献
[1]许林成 彭国勋.包装机械[M].长沙:湖南大学出版社,1989 .
[2]许林成.包装机械设计与原理[M].上海:上海科学技术出版社, 1988 .
[3]孙智慧 晏祖根.包装机械概论[M].北京：印刷工业出版社,2012,(7).
[4]肖乾 周新建.凸轮机构的计算机辅助设计与运动仿真分析[J]华东交通大学学报,2006 ,23 (4) :103 - 105 .
[5]周开利 邓春晖.MATLAB应用教程[M].北京：北京大学出版社，2007,3
[6]张义智 郭连考.给纸机递纸吸嘴机构共轭凸轮设计[J].包装工程,2007 ,28 (1) :104 - 105 .
[7]王家文 王皓 刘海.MATLAB7.0 编程基础[ M ].北京:机械工业出版社,2005 .
[8]孙智慧 高德.包装机械[M].北京：中国轻工业出版社,2010,(1).
[9]刘筱霞.包装机械[M].北京：化学工业出版社,2007,(1).
[10]高德.包装机械设计[M].北京:化学工业出版社,2005,(8).
[11]刘月花 郭谆钦.机械设计基础[M].长春：东北师范大学出版社，2010，（3）.

八、附录
附录程序一：
clc;clear;
global b1 b2 FAMAX FAMIN sita E
fa1=130*pi/180;fa2=0*pi/180;fa3=130*pi/180;fa4=100*pi/180;sm=120;
E=0;h=1.5*sm;a=0.8*sm;
sita=-8*pi/180;
fa01=0:0.001:fa1;
faa=fa01/fa1;
s=sm.*(6.*(faa).^5-15.*(faa).^4+10.*(faa).^3);
xb=a.*sin(sita+fa01);
yb=-a.*cos(sita+fa01);
xd=E;
yd=h+s;
b01=sqrt((xb-xd).^2+(yb-yd).^2);
figure;
plot(fa01,b01);
hold on;
fa02=fa2+fa1:0.001:fa3+fa1+fa2;
s=sm.*(1-(6.*((fa02-fa1-fa2)/fa3).^5-15.*((fa02-fa1-fa2)/fa3).^4+10.*((fa02-fa1-fa2)/fa3).^3));xb=a*sin(sita+fa02);
yb=-a*cos(sita+fa02);
xd=E;
yd=h+s;
b02=sqrt((xb-xd).^2+(yb-yd).^2);
plot(fa02,b02);
hold on;
bmin=min([b01,b02])
bmax=max([b01,b02])
b1=(bmax-bmin)/2
b2=(bmax+bmin)/2
famax=b01(1);
for i=2:2095
if famax<b01(i)
famax=b01(i);
FAMAX=fa01(i);

end
end
famax,FAMAX
famin=b02(1);
for i=2:2095
if famin>b02(i)
famin=b02(i);
FAMIN=fa02(i);

end
end
famin,FAMIN

附录程序二：
famax=1.3990;famin=3.8449;
b1=70.4613;b2=327.4171;
fa1=130*pi/180;fa2=0*pi/180;fa3=130*pi/180;fa4=100*pi/180;sm=100;
e=0;h=3*sm;a=0.8*sm;
sita=-4*pi/180;
fa01=0:0.002:famax-0.002;
faa=fa01/fa1;
s=sm.*(6.*(faa).^5-15.*(faa).^4+10.*(faa).^3);
xb=a.*sin(sita+fa01);
yb=-a.*cos(sita+fa01);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta1=asin((xb-xd)./b);
beta1=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta1+beta1);
yc=yd-b2*cos(theta1+beta1);
plot(xc,yc,'r');
hold on;
fa02=famax+0.002:0.002:fa1;
faa=fa02/fa1;
s=sm.*(6.*(faa).^5-15.*(faa).^4+10.*(faa).^3);
%plot(fa02,s)
xb=a.*sin(sita+fa02);
yb=-a.*cos(sita+fa02);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta2=asin((xb-xd)./b);
beta2=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta2-beta2);
yc=yd-b2*cos(theta2-beta2);
plot(xc,yc,'r');hold on;
fa03=fa2+fa1:0.002:famin;
s=sm.*(1-(6.*((fa03-fa1-fa2)/fa3).^5-15.*((fa03-fa1-fa2)/fa3).^4+10.*((fa03-fa1-fa2)/fa3).^3));%plot(fa03,s)
xb=a*sin(sita+fa03);
yb=-a*cos(sita+fa03);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta=asin((xb-xd)./b);
beta3=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta-beta3);
yc=yd-b2*cos(theta-beta3);
plot(xc,yc,'r');hold on;
fa04=famin+0.002:0.002:fa1+fa2+fa3;
s=sm.*(1-(6.*((fa04-fa1-fa2)/fa3).^5-15.*((fa04-fa1-fa2)/fa3).^4+10.*((fa04-fa1-fa2)/fa3).^3));%plot(fa04,s)
xb=a*sin(sita+fa04);
yb=-a*cos(sita+fa04);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta=asin((xb-xd)./b);
beta4=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta+beta4);
yc=yd-b2*cos(theta+beta4);
plot(xc,yc,'r');hold on;
fa05=fa1+fa2+fa3:0.002:2*pi;
s=0;%plot(fa05,s);
xb=a*sin(sita+fa05);
yb=-a*cos(sita+fa05);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta=asin((xb-xd)./b);
beta5=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta+beta5);
yc=yd-b2*cos(theta+beta5);
plot(xc,yc,'r');hold on;

附录程序三：
clc;clear;
global sita E h a sm
famax=1.4090;famin=3.7329;
b1=67.8307;b2=182.2473;
fa1=130*pi/180;fa2=0*pi/180;fa3=130*pi/180;fa4=100*pi/180;sm=100;
e=0;h=1.5*sm;a=0.8*sm;
sita=-2*pi/180;
fa01=0:0.001:famax;
faa=fa01/fa1;
s=sm.*(6.*(faa).^5-15.*(faa).^4+10.*(faa).^3);
xb=a.*sin(sita+fa01);
yb=-a.*cos(sita+fa01);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta1=asin((xb-xd)./b);
beta1=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta1+beta1);
yc=yd-b2*cos(theta1+beta1);
hold on;
k11=(yc-yb)./(xc-xb);
for i=2:length(fa01)-1
k21(i)=(yc(i-1)-yc(i+1))./(xc(i-1)-xc(i+1));
end
k21(1)=k21(2);
for i=1:length(fa01)-1
a1(i)=abs(atan((k21(i)-k11(i))./(1+k11(i)*k21(i))))*180/pi;
end
fa01(length(fa01))=[];
a1(length(fa01))=a1(length(fa01)-1);
a1;
plot(fa01,a1);hold on;
fa02=famax:0.001:fa1;
faa=fa02/fa1;
s=sm.*(6.*(faa).^5-15.*(faa).^4+10.*(faa).^3);
%plot(fa02,s)
xb=a.*sin(sita+fa02);
yb=-a.*cos(sita+fa02);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta2=asin((xb-xd)./b);
beta2=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta2-beta2);
yc=yd-b2*cos(theta2-beta2);
hold on;
k12=(yc-yb)./(xc-xb);
for i=2:length(fa02)-1
k22(i)=(yc(i-1)-yc(i+1))./(xc(i-1)-xc(i+1));
end
k22(1)=k22(2);
for i=1:length(fa02)-1
a2(i)=abs(atan((k22(i)-k12(i))./(1+k12(i)*k22(i))))*180/pi;
end
fa02(length(fa02))=[];
a2(length(fa02))=a2(length(fa02)-1);
a2;
plot(fa02,a2);hold on;
fa03=fa2+fa1:0.001:famin;
s=sm.*(1-(6.*((fa03-fa1-fa2)/fa3).^5-15.*((fa03-fa1-fa2)/fa3).^4+10.*((fa03-fa1-fa2)/fa3).^3));%plot(fa03,s)
xb=a*sin(sita+fa03);
yb=-a*cos(sita+fa03);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta=asin((xb-xd)./b);
beta3=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta-beta3);
yc=yd-b2*cos(theta-beta3);
hold on;

k13=(yc-yb)./(xc-xb);
for i=2:length(fa03)-1
k23(i)=(yc(i-1)-yc(i+1))./(xc(i-1)-xc(i+1));
end
k23(1)=k23(2);
for i=1:length(fa03)-1
a3(i)=abs(atan((k23(i)-k13(i))./(1+k13(i)*k23(i))))*180/pi;
end
fa03(length(fa03))=[];
a3(length(fa03))=a3(length(fa03)-1);
a3;
plot(fa03,a3);hold on;
fa04=famin:0.001:fa1+fa2+fa3;
s=sm.*(1-(6.*((fa04-fa1-fa2)/fa3).^5-15.*((fa04-fa1-fa2)/fa3).^4+10.*((fa04-fa1-fa2)/fa3).^3));%plot(fa04,s)
xb=a*sin(sita+fa04);
yb=-a*cos(sita+fa04);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta=asin((xb-xd)./b);
beta4=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta+beta4);
yc=yd-b2*cos(theta+beta4);
hold on;
k14=(yc-yb)./(xc-xb);
for i=2:length(fa04)-1
k24(i)=(yc(i-1)-yc(i+1))./(xc(i-1)-xc(i+1));
end
k24(1)=k24(2);
for i=1:length(fa04)-1
a4(i)=abs(atan((k24(i)-k14(i))./(1+k14(i)*k24(i))))*180/pi;
end
fa04(length(fa04))=[];
a4(length(fa04))=a4(length(fa04)-1);
a4;
plot(fa04,a4);hold on;
fa05=fa1+fa2+fa3:0.001:2*pi;
s=0;%plot(fa05,s);
xb=a*sin(sita+fa05);
yb=-a*cos(sita+fa05);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta=asin((xb-xd)./b);
beta5=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta+beta5);
yc=yd-b2*cos(theta+beta5);
hold on;

k15=(yc-yb)./(xc-xb);
for i=2:length(fa05)-1
k25(i)=(yc(i-1)-yc(i+1))./(xc(i-1)-xc(i+1));
end
k25(1)=k25(2);
for i=1:length(fa05)-1
a5(i)=abs(atan((k25(i)-k15(i))./(1+k15(i)*k25(i))))*180/pi;
end
fa05(length(fa05))=[];
a5(length(fa05))=a5(length(fa05)-1);
a5;
plot(fa05,a5);hold on;
max([a1,a2,a3,a4,a5])
min([a1,a2,a3,a4,a5])

附录程序四：
clc;clear;
global sita E h a sm
famax=1.3990;famin=3.8449;
b1=70.4613;b2=327.4171;
fa1=130*pi/180;fa2=0*pi/180;fa3=130*pi/180;fa4=100*pi/180;sm=100;
e=0;h=3*sm;a=0.8*sm;
sita=-4*pi/180;
fa01=0:0.002:famax-0.002;
faa=fa01/fa1;
s=sm.*(6.*(faa).^5-15.*(faa).^4+10.*(faa).^3);
xb=a.*sin(sita+fa01);
yb=-a.*cos(sita+fa01);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta1=asin((xb-xd)./b);
beta1=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta1+beta1);
yc=yd-b2*cos(theta1+beta1);
hold on;
%+++++++++++++++++++++++++++++++++++从动件压力角——红色
c1=abs(theta1-beta1)*180/pi;
plot(fa01,c1,'r--*');hold on;
fa02=famax+0.002:0.002:fa1;
faa=fa02/fa1;
s=sm.*(6.*(faa).^5-15.*(faa).^4+10.*(faa).^3);
%plot(fa02,s)
xb=a.*sin(sita+fa02);
yb=-a.*cos(sita+fa02);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta2=asin((xb-xd)./b);
beta2=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta2-beta2);
yc=yd-b2*cos(theta2-beta2);
hold on;
c2=abs(theta2+beta2)*180/pi;
plot(fa02,c2,'r--*');hold on;
fa03=fa2+fa1:0.002:famin;
s=sm.*(1-(6.*((fa03-fa1-fa2)/fa3).^5-15.*((fa03-fa1-fa2)/fa3).^4+10.*((fa03-fa1-fa2)/fa3).^3));%plot(fa03,s)
xb=a*sin(sita+fa03);
yb=-a*cos(sita+fa03);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta3=asin((xb-xd)./b);
beta3=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(theta3-beta3);
yc=yd-b2*cos(theta3-beta3);
hold on;
c3=abs(theta3+beta3)*180/pi;
plot(fa03,c3,'r--*');hold on;
fa04=famin+0.002:0.002:fa1+fa2+fa3;
s=sm.*(1-(6.*((fa04-fa1-fa2)/fa3).^5-15.*((fa04-fa1-fa2)/fa3).^4+10.*((fa04-fa1-fa2)/fa3).^3));%plot(fa04,s)
xb=a*sin(sita+fa04);
yb=-a*cos(sita+fa04);
xd=e;
yd=h+s;
b=sqrt((xb-xd).^2+(yb-yd).^2);
theta4=asin((xb-xd)./b);
beta4=acos((b.^2+b2.^2-b1.^2)./(2.*b.*b2));
xc=xd+b2*sin(t