k=1.25; % 行程速比系数 theta=pi*(k-1)/(k+1); % 极位夹角 yg=250; % 摇杆长度 pusai=pi/6; % 摇杆摆角 gamin=2*pi/9; % 最小传动角 x=fsolve(@qbyg,x0);
disp ' ******** 已知条件 ********' fprintf (1,' 行程速比系数 k = %3.4f \\n',k)
fprintf (1,' 极位夹角 theta = %3.4f 度 \\n',theta*180/pi) fprintf (1,' 摇杆长度 yg = %3.4f mm \\n',yg)
fprintf (1,' 摇杆摆角 pusai = %3.4f 度 \\n',pusai*180/pi) fprintf (1,' 最小传动角 gamin = %3.4f 度 \\n',gamin*180/pi) disp ' ******** 计算结果 ********'
fprintf (1,' 曲柄长度 a = %3.4f mm \\n',x(1)) fprintf (1,' 连杆长度 b = %3.4f mm \\n',x(2)) fprintf (1,' 机架长度 d = %3.4f mm \\n',x(3))
fprintf (1,' 摇杆位置角 gamin = %3.4f 度 \\n',x(4)*180/pi)
% 铰链四杆机构非线性参数方程组 function f=qbyg(x)
k=1.25; % 行程速比系数 theta=pi*(k-1)/(k+1); % 极位夹角 yg=250; % 摇杆长度 pusai=pi/6; % 摇杆摆角 gamin=2*pi/9; % 最小传动角
% x(1)是曲柄长度;x(2)是连杆长度;x(3)是机架长度;x(4)是摇杆初始位置角
f1=(x(2)+x(1))^2+(x(2)-x(1))^2-2*(x(2)+x(1))*(x(2)-x(1))*cos(theta)-(2*yg*sin(pusai/2))^2; f2=yg^2+x(3)^2-2*yg*x(3)*cos(x(4))-(x(2)-x(1))^2;
f3=yg^2+x(3)^2-2*yg*x(3)*cos(x(4)+pusai)-(x(2)+x(1))^2; f4=yg^2+x(2)^2-2*yg*x(2)*cos(gamin)-(x(3)-x(1))^2; f=[f1;f2;f3;f4];
% 实现连架杆角位移(3组)的连杆机构运动设计 % 已知条件
f0=0;p0=0; % 连架杆初始位置角 [f]=[45 90 135]*pi/180; % 曲柄输入角 [p]=[52 82 112]*pi/180; % 摇杆输出角
% 杆件相对长度参数R1、R2和R3的系数矩阵 a1=[1 -cos(f(1)+f0) cos(p(1)+p0)]; a2=[1 -cos(f(2)+f0) cos(p(2)+p0)]; a3=[1 -cos(f(3)+f0) cos(p(3)+p0)]; a=[a1;a2;a3]
% 线性方程组右边的常数矩阵
b1=[cos(f(1)-p(1))+(f0+p0)]; b2=[cos(f(2)-p(2))+(f0+p0)]; b3=[cos(f(3)-p(3))+(f0+p0)]; b=[b1 b2 b3]'
% 线性方程组aR=b直接解法(采用求逆函数inv)
R=inv(a)*b % 或采用矩阵除法 R=a\\b % 杆件长度
d=50; % 机架长度(已知数据) x(1)=d/R(3); x(3)=d/R(2);
x(2)=sqrt(x(1)^2+x(3)^2+d^2-2*x(1)*x(3)*R(1));
% 检验解的精度(采用求解矩阵或向量范数的函数norm) en=norm(a*R-b);
disp ' ******** 计算结果 ********'
fprintf (1,' 曲柄长度 a = %3.4f mm \\n',x(1)) fprintf (1,' 连杆长度 b = %3.4f mm \\n',x(2)) fprintf (1,' 摇杆长度 c = %3.4f mm \\n',x(3)) fprintf (1,' 机架长度 d = %3.4f mm \\n',d) disp ' '
fprintf (1,' 数值解的精度 en = %3.4e \\n',en) % 连杆机构实现函数的优化设计
% (调用目标函数jfg_f和非线性约束函数cdj_g) % 1----优化设计主程序
% 设计变量的初始值和已知杆件长度(曲柄和机架) x0=[4.5;4]; qb=1;jj=5;
% 设计变量的下界与上界 lb=[1;1]; ub=[10;10];
% 线性不等式约束(g1(x)和g2(x))中设计变量的系数矩阵 a=zeros(2,2);
a(1,1)=-1;a(2,1)= 1; a(3,2)=-1;a(4,2)= 1;
% 线性不等式约束(g1(x)和g2(x))中的常数项列阵 b=[-1;10;-1;10];
% 使用约束优化命令fmincon(调用目标函数jfg_f和非线性约束函数cdj_g) % 没有等式约束,参数Aeq和beq定义为空矩阵符号“[ ]” [x,fn]=fmincon(@jfg_f,x0,a,b,[],[],lb,ub,@cdj_g);
disp ' ******** 连杆机构实现函数优化设计最优解 ********' fprintf (1,' 连杆相对长度 a = %3.4f \\n',x(1)) fprintf (1,' 摇杆相对长度 b = %3.4f \\n',x(2)) fprintf (1,' 输出角平方误差之和 f* = %3.4f \\n',fn)
% 调用约束优化非线性约束函数(cdj_g)计算最优点x*的性能约束函数值 g=cdj_g(x);
disp ' ======== 最优点的性能约束函数值 ========' fprintf (1,' 最小传动角约束函数值 g1* = %3.4f \\n',g(1)) fprintf (1,' 最大传动角约束函数值 g2* = %3.4f \\n',g(2))
% 2----连杆机构实现函数优化的目标函数(jfg_f) function f=jfg_f(x); s=50;qb=1;jj=5;fx=0;
fa0=acos(((qb+x(1))^2-x(2)^2+jj^2)/(2*(qb+x(1))*jj)); % 曲柄初始角 pu0=acos(((qb+x(1))^2-x(2)^2-jj^2)/(2*x(2)*jj)); % 摇杆初始角 for i=1:s
fai=fa0+0.5*pi*i/s;
pui=pu0+2*(fai-fa0)^2/(3*pi);
ri=sqrt(qb^2+jj^2-2*qb*jj*cos(fai));
alfi=acos((ri^2+x(2)^2-x(1)^2)/(2*ri*x(2))); bati=acos((ri^2+jj^2-qb^2)/(2*ri*jj)); if fai>0 & fai<=pi psi=pi-alfi-bati;
elseif fai>pi & fai<=2*pi psi=pi-alfi+bati; end
fx=fx+(pui-psi)^2; % 输出角平方误差之和 end f=fx;
% 3----连杆机构实现函数优化的非线性不等式约束函数(cdj_g) function [g,ceq]=cdj_g(x); qb=1;jj=5;
gamn=45*pi/180;gamm=135*pi/180;
g(1)=x(1)^2+x(2)^2-(jj-qb)^2-2*x(1)*x(2)*cos(gamn); % 最小传动角约束 g(2)=-x(1)^2-x(2)^2+(jj+qb)^2-2*x(1)*x(2)*cos(gamn); % 最大传动角约束 ceq=[];
% 曲柄摇杆机构运动分析
% (1)-----计算连杆的输出角th3和摇杆的输出角th4 % 设定各杆的长度(单位:毫米)
rs(1)=304.8; % 设定机架1长度 rs(2)=101.6; % 设定曲柄2长度 rs(3)=2.0; % 设定连杆3长度 rs(4)=177.8; % 设定摇杆4长度 dr=pi/180.0;% 角度与弧度的转换系数 % 设定初始推测的输入 % 机构的初始位置
th(1)=0.0; % 设定曲柄2初始位置角是0度(与机架1共线) th(2)=45*dr; % 连杆3的初始位置角是 45度 th(3)=135*dr; % 摇杆4的初始位置角是135度
% 摇杆4的初始位置角可以用三角形的正弦定理确定 th(3)=pi-asin(sin(th(2))*rs(3)/rs(4))
dth=5*dr; % 循环增量 % 曲柄输入角从0度变化到360度,步长为5度,计算th34 for i=1:72
[th3,th4]=ntrps(th,rs); % 调用牛顿—辛普森方程求解机构位置解非线性方程函数文件 % Store results in a matrix-th34,in degrees
% 在矩阵th34中储存结果,以度为单位;(i,:)表示第i行所有列的元素;(:,i)表示第i列所有行的元素
th34(i,:)=[th(1)/dr th3/dr th4/dr]; % 矩阵[曲柄转角 连杆转角 摇杆转角] th(1)=th(1)+dth; % 曲柄转角递增
th(2)=th3; % 连杆转角中间计算值 th(3)=th4; % 摇杆转角中间计算值 end
% 绘制输出角th(2)与th(3)—输入角th(1)的关系曲线
subplot(2,2,1) % 选择第1个子窗口 plot(th34(:,1),th34(:,2),th34(:,1),th34(:,3)) axis([0 360 0 170])
grid % 网格线 ylabel('从动件角位移/deg') title('角位移线图')
text(110,110,'摇杆4角位移') text(50,35,'连杆3角位移')
% (2)-----计算连杆的角速度om3和摇杆的角速度om4 % Setting initial conditions % 设置初始条件
om2=250; % 曲柄角速度(等速输入)
T=2*pi/om2; % 机构周期-曲柄旋转1周的时间(秒) % 曲柄输入角从0度变化到360度,步长为5度,计算om34 for i=1:72
ct(2)=i*dth;
A=[-rs(3)*sin(th34(i,2)*dr) rs(4)*sin(th34(i,3)*dr); rs(3)*cos(th34(i,2)*dr) -rs(4)*cos(th34(i,3)*dr)];
B=[om2*rs(2)*sin(ct(2));-om2*rs(2)*cos(ct(2))];
om=inv(A)*B; % 输出角速度矩阵 om3=om(1); om4=om(2);
om34(i,:)=[i om3 om4]; % 矩阵[序号 连杆角速度 摇杆角速度] t(i)=i*T/72; end
% 绘制连杆的角速度om3和摇杆的角速度om4—时间Times的关系曲线 subplot(2,2,2) % 选择第2个子窗口 plot(t,om34(:,2),t,om34(:,3)) axis([0 0.026 -190 210])
grid % 网格线 title('角速度线图')
ylabel('从动件角速度/rad/s') text(0.001,170,'摇杆4角速度') text(0.013,130,'连杆3角速度')
% (3)-----计算连杆的角加速度a3和摇杆的角加速度a4
a2=0; % 曲柄角速度是等速,角加速度a2=dom2/dt=0 % 曲柄输入角从0度变化到360度,步长为5度,计算a34 for i=1:72
c(2)=i*dth;
C=[-rs(3)*sin(th34(i,2)*dr) rs(4)*sin(th34(i,3)*dr); rs(3)*cos(th34(i,2)*dr) -rs(4)*cos(th34(i,3)*dr)]; D(1)=
a2*rs(2)*sin(c(2))+om2^2*rs(2)*cos(c(2))+om34(i,2)^2*rs(3)*cos(th34(i,2)*dr)-om34(i,3)^2*rs(4)*cos(th34(i,3)*dr);
D(2)=-a2*rs(2)*cos(c(2))+om2^2*rs(2)*sin(c(2))+om34(i,2)^2*rs(3)*sin(th34(i,2)*dr)-om34(i,3)^2*rs(4)*sin(th34(i,3)*dr);
a=inv(C)*D'; % 输出角加速度矩阵 a3=a(1); a4=a(2);
a34(i,:)=[i a3 a4]; % 矩阵[序号 连杆角加速度 摇杆加角速度] t(i)=i*T/72; end
% 绘制连杆的角加速度a3和摇杆的角加速度a4—时间Times的关系曲线 subplot(2,2,3) % 选择第3个子窗口 plot(t,a34(:,2),t,a34(:,3)) axis([0 0.026 -6*1e4 8*1e4])
grid % 网格线 title('角加速度线图') xlabel('时间/s')
ylabel('从动件加速度/rad/s^{2}') text(0.003,6.2*1e4,'摇杆4角加速度') text(0.010,3.3*1e4,'连杆3角加速度') %
% 输出1：四杆机构运动周期(0:5:360),时间,角位移,角速度,角加速度数据
disp ' 曲柄转角 连杆转角-摇杆转角-连杆角速度-摇杆角速度-连杆加速度-摇杆加速度' ydcs=[th34(:,1),th34(:,2),th34(:,3),om34(:,2),om34(:,3),a34(:,2),a34(:,3)]; disp (ydcs)
% 输出参数的数量级必须一致 %
% (4)-----运动误差分析
% 闭环矢量方程:r2+r3-r4-r1=0
% 误差矢量E=r2+r3-r4-r1的模是表示仿真有效程度的标量(ex和ey是误差分量) ex=rs(2)*cos(th34(:,1)*dth)+rs(3)*cos(th34(:,2)*dth)-rs(4)*cos(th34(:,3)*dth)-rs(1);
ey=rs(2)*sin(th34(:,2)*dth)+rs(3)*sin(th34(:,2)*dth)-rs(4)*sin(th34(:,3)*dth); ee=norm([ex ey]); % 计算误差矢量矩阵的范数(模) %
% 输出2：四杆机构运动周期(0:5:360),时间,X向误差分量,Y向误差分量 disp ' 曲柄转角 时间(秒) X向误差 Y向误差' wc=[th34(:,1),t(:),ex(:,1),ey(:,1)]; disp (wc)
fprintf (1,' 误差矢量矩阵的模 ee = %3.4f \\n',ee) %
% 绘制均方根相容性误差曲线
subplot(2,2,4) % 选择第4个子窗口 plot(t,ex(:,1),t,ey(:,1)) axis([0 0.026 -800 600])
grid % 网格线 title('均方根误差曲线') xlabel('时间/s')
ylabel('均方根误差')
text(0.012,350,'X向误差分量') text(0.003,-600,'Y向误差分量')
disp ' ****** 曲柄滑块机构的运动学分析 *******' SC=input(' 输入滑块行程的均值(mm) SC = '); P=input(' 输入曲柄轴心至滑销最远距离(mm) P = '); E=input(' 输入机构偏心距的均值(mm) E = '); RL=input(' 输入曲柄与连杆长度比的均值 RL = '); DR=input(' 输入曲柄长度偏差(mm) DR = '); DL=input(' 输入连杆长度偏差(mm) DL = '); DE=input(' 输入机构偏心距偏差(mm) DE = '); N=input(' 输入曲柄转速(r/min) N = '); L=sqrt((P-SC)^2-E^2)/(1-RL);
fprintf(1,' 连杆长度的均值(mm) L = %3.6f \\n',L) R=RL*L;
fprintf(1,' 曲柄长度的均值(mm) R = %3.6f \\n',R) CR=DR/3;CL=DL/3;CE=DE/3; EL=E/L;
fprintf(1,' 偏心距与连杆长度比的均值(mm) EL = %3.6f \\n',EL) fprintf(1,' 曲柄长度的标准离差(mm) CR = %3.6f \\n',CR) fprintf(1,' 连杆长度的标准离差(mm) CR = %3.6f \\n',CL) fprintf(1,' 偏心距的标准离差(mm) CE = %3.6f \\n',CE) W=pi*N/30;
fprintf(1,' 曲柄的角速度(mm) W = %3.6f \\n',W) CRL=sqrt((R*CL)^2+(L*CR)^2)/L^2;
fprintf(1,' 曲柄与连杆长度比的标准离差 CRL = %3.6f \\n',CRL) CEL=sqrt((E*CL)^2+(L*CE)^2)/L^2;
fprintf(1,' 偏心距与连杆长度比的标准离差 CEL = %3.6f \\n',CEL)
theta=0:10:360; hd=theta.*pi/180;
% 计算滑块位移、速度、加速度的均值
S=R.*(1-cos(hd)-EL.*sin(hd)+0.5.*RL.*sin(hd).^2); V=R.*W.*(sin(hd)-EL.*cos(hd)+0.5.*RL.*sin(2.*hd)); A=R.*W^2.*(cos(hd)+EL.*sin(hd)+RL.*cos(2.*hd)); figure(1);
subplot(1,3,1); plot(theta,S,'r')
title('\\bf \\mus 线图') subplot(1,3,2); plot(theta,V,'k')
title('\\bf \\muv 线图')
xlabel('\\bf 曲柄转角\heta(度)') subplot(1,3,3); plot(theta,A,'b')
title('\\bf \\mua 线图')
% 计算滑块位移、速度、加速度的标准离差
CS=sqrt((1-cos(hd)+(0.5.*RL.*sin(hd)-EL).*sin(hd)).^2.*CR^2+(0.5.*(CRL.*sin(hd)).^2-CEL^2).^2.*(R.*sin(hd)).^2);
CV=W.*sqrt((sin(hd)-EL.*cos(hd)+0.5.*RL.*sin(2.*hd)).^2.*CR^2+(0.5.*(CRL.*sin(hd)).^2-(CEL.*cos(hd)).^2)*R^2);
CA=W^2.*sqrt((cos(hd)+EL.*sin(hd)+RL.*cos(2.*hd)).^2.*CR^2+(R.*CRL.*sin(hd)).^2+(RL.*CEL.*cos(2.*hd)).^2); figure(2);
subplot(1,3,1); plot(theta,CS,'r')
title('\\bf \\sigmas 线图') subplot(1,3,2); plot(theta,CV,'k')
title('\\bf \\sigmav 线图')
xlabel('\\bf 曲柄转角\heta(度)') subplot(1,3,3); plot(theta,CA,'b')
title('\\bf \\sigmaa 线图')
% 计算滑块位移、速度、加速度的偏差 DS=3.*CS;DV=3.*CV;DA=3.*CA;
% 计算滑块位移、速度、加速度的最大值和最小值 SM=S+DS;SN=S-DS; VM=V+DV;VN=V-DV; AM=A+DA;AN=A-DA;
% 计算滑块位移、速度、加速度的差值 SD=2.*DS;VD=2.*DV;AD=2.*DA; figure(3);
subplot(1,3,1); plot(theta,SD,'r')
title('\\bf \\Deltas 线图') subplot(1,3,2); plot(theta,VD,'k')
title('\\bf \\Deltav 线图')
xlabel('\\bf 曲柄转角\heta(度)') subplot(1,3,3); plot(theta,AD,'b')
title('\\bf \\Deltaa 线图')
disp ' ****** 曲柄滑块机构的等影响法精度综合 *******' N=input(' 输入机构运动精度影响尺度数目 N = '); H=input(' 输入滑块行程的均值(mm) H = '); P=input(' 输入曲柄轴心至滑销最远距离(mm) P = '); DH=input(' 输入滑块位置允许误差(mm) DH = '); R=H/2;
fprintf(1,' 曲柄长度的均值(mm) R = %3.3f \\n',R) L=P-R;
fprintf(1,' 连杆长度的均值(mm) L = %3.3f \\n',L) theta=0:10:360; hd=theta.*pi/180;
% 计算曲柄长度和滑块长度的影响系数(偏导数的最大绝对值) CR=1-cos(hd);
CL=0.5.*sin(hd).^2;
CRM=max(abs(1-cos(hd)));
CLM=max(abs(0.5.*sin(hd).^2));
fprintf(1,' 曲柄长度影响系数的最大绝对值 CRM = %3.6f \\n',CRM) fprintf(1,' 连杆长度影响系数的最大绝对值 CLM = %3.6f \\n',CLM) % 计算曲柄长度和滑块长度的最大允许偏差 DRM=DH/sqrt(N)/CRM; DLM=DH/sqrt(N)/CLM;
fprintf(1,' 曲柄长度允许的最大偏差(mm) DRM = %3.6f \\n',DRM) fprintf(1,' 连杆长度允许的最大偏差(mm) DLM = %3.6f \\n',DLM) plot(theta,CR,'r')
gtext('曲柄长度影响系数曲线') title('\\bf 机构尺度影响系数线图') xlabel('\\bf 曲柄转角\heta(度)') ylabel('\\bf 尺度影响系数') hold;
plot(theta,CL,'k')
gtext('连杆长度影响系数曲线')
% 外槽轮机构运动分析
dr=pi/180.0; % 角度与弧度的转换系数
% 销轮2转角范围:-f20 lmd=sin(pi/z); % 计算曲柄2与机架1的长度比 bc=10; % 循环步长 cz=-f20/dr; % 循环初值 zz=f20/dr; % 循环终值 i=1; % 根据步长变化的运动参数矩阵cs行数计数器 for f2=cz:bc:zz % 计算槽轮角位移、类角速度、类角加速度 wy=atan(lmd*sin(f2*dr)/(1-lmd*cos(f2*dr))); sd=lmd*(cos(f2*dr)-lmd)/(1-2*lmd*cos(f2*dr)+lmd^2); jsd=-lmd*sin(f2*dr)*(1-lmd^2)/(1-2*lmd*cos(f2*dr)+lmd^2)^2; switch z % 矩阵c(i,:)表示第i行的各列元素 case 4,c4(i,:)=[f2 wy/dr sd jsd]; case 6,c6(i,:)=[f2 wy/dr sd jsd]; case 8,c8(i,:)=[f2 wy/dr sd jsd]; case 10,c10(i,:)=[f2 wy/dr sd jsd]; end i=i+1; end end % 输出外槽轮机构运动参数 ['轮槽数 z=4'] [' 销轮转角',' 槽轮角位移',' 角速度',' 角加速度'] % 矩阵c(:,j)表示第j列的各行元素 [c4(:,1),c4(:,2),c4(:,3),c4(:,4)] ['轮槽数 z=6'] [' 销轮转角',' 槽轮角位移',' 角速度',' 角加速度'] % 矩阵c(:,j)表示第j列的各行元素 [c6(:,1),c6(:,2),c6(:,3),c6(:,4)] ['轮槽数 z=8'] [' 销轮转角',' 槽轮角位移',' 角速度',' 角加速度'] % 矩阵c(:,j)表示第j列的各行元素 [c8(:,1),c8(:,2),c8(:,3),c8(:,4)] ['轮槽数 z=10'] [' 销轮转角',' 槽轮角位移',' 角速度',' 角加速度'] % 矩阵c(:,j)表示第j列的各行元素 [c10(:,1),c10(:,2),c10(:,3),c10(:,4)] % % 绘制槽轮机构运动参数曲线 figure(1); % 生成槽轮运动线图窗口 subplot(2,2,1); % 选择第1个子窗口 plot(c4(:,1),c4(:,3),c4(:,1),c4(:,4)) % 绘制z= 4的线图 title('外槽轮槽数 z=4') % 标注子窗口名称 axis([-pi/4/dr pi/4/dr -6 6]) % 定义坐标轴范围 grid % 栅格线 text(-2,4.2,'\\epsilon/\\omega^{2}') % 标注类角加速度线图 text(20,1.6,'\\omega/\\omega') % 标注类角速度线图 ylabel('槽轮运动线图') % 定义纵坐标轴名称 % subplot(2,2,2); % 选择第2个子窗口 plot(c6(:,1),c6(:,3),c6(:,1),c6(:,4)) % 绘制z= 6的线图 title('外槽轮槽数 z=6') axis([-pi/3/dr pi/3/dr -1.5 1.5]) grid text(10,-0.7,'\\epsilon/\\omega^{2}') text(30,0.7,'\\omega/\\omega') ylabel('槽轮运动线图') % subplot(2,2,3); % 选择第3个子窗口 plot(c8(:,1),c8(:,3),c8(:,1),c8(:,4)) % 绘制z= 8的线图 title('外槽轮槽数 z=8') axis([-3*pi/8/dr 3*pi/8/dr -0.8 0.8]) grid text(12,-0.3,'\\epsilon/\\omega^{2}') text(40,0.4,'\\omega/\\omega') ylabel('槽轮运动线图') % subplot(2,2,4); % 选择第4个子窗口 plot(c10(:,1),c10(:,3),c10(:,1),c10(:,4)) % 绘制z=10的线图 title('外槽轮槽数 z=10') axis([-2*pi/5/dr 2*pi/5/dr -0.5 0.5]) grid text(15,-0.2,'\\epsilon/\\omega^{2}') text(40,0.3,'\\omega/\\omega') ylabel('槽轮运动线图') % figure(2); % 生成类线图窗口 subplot(1,2,1); % 选择第1个子窗口 plot(c4(:,1),c4(:,3),c6(:,1),c6(:,3),c8(:,1),c8(:,3),c10(:,1),c10(:,3)) title('\\omega/\\omega') axis([-f20/dr f20/dr -0.1 2.5]) grid text(-10,0.35,'z=10') text(-8,0.7,'z=8') text(-8,1.1,'z=6') text(10,2.1,'z=4') ylabel('槽轮类角速度线图') % subplot(1,2,2); % 选择第2个子窗口 plot(c4(:,1),c4(:,4),c6(:,1),c6(:,4),c8(:,1),c8(:,4),c10(:,1),c10(:,4)) title('\\epsilon/\\omega^{2}') axis([-f20/dr f20/dr -5.5 5.5]) grid text(-50,0.2,'z=10') text(-30,0.9,'z=8') text(-25,1.6,'z=6') text(0,3.5,'z=4') ylabel('槽轮类角加速度线图') % 内槽轮机构运动分析 dr=pi/180.0; % 角度与弧度的转换系数 % 销轮2转角范围:-f20 lmd=sin(pi/z); % 计算曲柄2与机架1的长度比 bc=10; % 循环步长 cz=-f20/dr; % 循环初值 zz=f20/dr; % 循环终值 i=1; % 根据步长变化的运动参数矩阵cs行数计数器 for f2=cz:bc:zz % 计算槽轮角位移、类角速度、类角加速度 wy=atan(lmd*sin(f2*dr)/(1+lmd*cos(f2*dr))); sd=lmd*(cos(f2*dr)+lmd)/(1+2*lmd*cos(f2*dr)+lmd^2); jsd=lmd*sin(f2*dr)*(1-lmd^2)/(1+2*lmd*cos(f2*dr)+lmd^2)^2; switch z % 矩阵c(i,:)表示第i行的各列元素 case 4,c4(i,:)=[f2 wy/dr sd jsd]; case 6,c6(i,:)=[f2 wy/dr sd jsd]; case 8,c8(i,:)=[f2 wy/dr sd jsd]; case 10,c10(i,:)=[f2 wy/dr sd jsd]; end i=i+1; end end % 输出内槽轮机构运动参数 ['轮槽数 z=4'] [' 销轮转角',' 槽轮角位移',' 角速度',' 角加速度'] % 矩阵c(:,j)表示第j列的各行元素 [c4(:,1),c4(:,2),c4(:,3),c4(:,4)] ['轮槽数 z=6'] [' 销轮转角',' 槽轮角位移',' 角速度',' 角加速度'] % 矩阵c(:,j)表示第j列的各行元素 [c6(:,1),c6(:,2),c6(:,3),c6(:,4)] ['轮槽数 z=8'] [' 销轮转角',' 槽轮角位移',' 角速度',' 角加速度'] % 矩阵c(:,j)表示第j列的各行元素 [c8(:,1),c8(:,2),c8(:,3),c8(:,4)] ['轮槽数 z=10'] [' 销轮转角',' 槽轮角位移',' 角速度',' 角加速度'] % 矩阵c(:,j)表示第j列的各行元素 [c10(:,1),c10(:,2),c10(:,3),c10(:,4)] % % 绘制槽轮机构运动参数曲线 figure(1); subplot(2,2,1); plot(c4(:,1),c4(:,3),c4(:,1),c4(:,4)) title('内槽轮槽数 z=4') axis([-3*pi/4/dr 3*pi/4/dr -1 1]) grid text(-85,-0.2,'\\epsilon/\\omega^{2}') text(10,0.55,'\\omega/\\omega') ylabel('槽轮运动线图') % subplot(2,2,2); plot(c6(:,1),c6(:,3),c6(:,1),c6(:,4)) title('内槽轮槽数 z=6') axis([-3*pi/4/dr 3*pi/4/dr -0.6 0.6]) grid text(-85,-0.3,'\\epsilon/\\omega^{2}') text(10,0.4,'\\omega/\\omega') ylabel('槽轮运动线图') % subplot(2,2,3); plot(c8(:,1),c8(:,3),c8(:,1),c8(:,4)) title('内槽轮槽数 z=8') axis([-3*pi/4/dr 3*pi/4/dr -0.40 0.40]) grid text(-85,-0.3,'\\epsilon/\\omega^{2}') text(10,0.32,'\\omega/\\omega') ylabel('槽轮运动线图') % subplot(2,2,4); plot(c10(:,1),c10(:,3),c10(:,1),c10(:,4)) title('内槽轮槽数 z=10') % 生成槽轮运动线图窗口 % 选择第1个子窗口 % 绘制z= 4的线图 % 标注子窗口名称 % 定义坐标轴范围 % 栅格线 % 标注类角加速度线图 % 标注类角速度线图 % 定义纵坐标轴名称 % 选择第2个子窗口 % 绘制z= 6的线图 % 选择第3个子窗口 % 绘制z= 8的线图 % 选择第4个子窗口 % 绘制z=10的线图 axis([-3*pi/4/dr 3*pi/4/dr -0.35 0.35]) grid text(-80,-0.22,'\\epsilon/\\omega^{2}') text(10,0.28,'\\omega/\\omega') ylabel('槽轮运动线图') % figure(2); % 生成类线图窗口 subplot(1,2,1); % 选择第1个子窗口 plot(c4(:,1),c4(:,3),c6(:,1),c6(:,3),c8(:,1),c8(:,3),c10(:,1),c10(:,3)) title('\\omega/\\omega') axis([-3*pi/4/dr 3*pi/4/dr -0.05 0.45]) grid text(-12,0.43,'z=4') text(-12,0.35,'z=6') text(-12,0.29,'z=8') text(-15,0.21,'z=10') ylabel('槽轮类角速度线图') % subplot(1,2,2); % 选择第2个子窗口 plot(c4(:,1),c4(:,4),c6(:,1),c6(:,4),c8(:,1),c8(:,4),c10(:,1),c10(:,4)) title('\\epsilon/\\omega^{2}') axis([-3*pi/4/dr 3*pi/4/dr -0.6 0.6]) grid text(40,0.22,'z=10') text(-60,-0.15,'z=8') text(-110,-0.45,'z=6') text(60,0.05,'z=4') ylabel('槽轮类角加速度线图') % 对心直动凸轮机构压力角的计算(调用TLYLJ.M) disp ' ' disp ' ******** 对心直动凸轮机构压力角的计算 ********' disp ' ' disp ' ======== 已 知 条 件 ========' disp ' ' rb = input(' 基圆半径(mm) rb = '); h = input(' 推程升程(mm) h = '); k=h/rb;hd=pi/180; fai = input(' 推程运动角(度) fai = '); fprintf (1,' 运动结构系数 k = %3.4f \\n',k) YDGL = input(' 运动规律类型：等速-\"ZX\";等加减速-\"PW\";余弦加速-\"JX\";正弦加速-\"BX\" == '); disp ' ' if YDGL=='ZX' disp ' ======== 等速运动（直线）规律 ========' fm=0; alfm=atan(k/(fai*hd)); elseif YDGL=='PW' disp ' ======== 等加减速运动（抛物线）规律 ========' if k<=2 fm=fai*hd/2; alfm=atan(4*k/(fai*hd*(2+k))); elseif k>2 fm=fai*hd/sqrt(2*k); alfm=atan(sqrt(2*k)/(fai*hd)); end elseif YDGL=='JX' disp ' ======== 余弦加速度运动（简谐曲线）规律 ========' fm=fai*hd*acos(k/(2+k))/pi; alfm=atan(k*pi/(2*fai*hd*sqrt(1+k))); elseif YDGL=='BX' disp ' ======== 正弦加速度运动（摆线）规律 ========' x=fsolve(@TLYLJ,fai*hd/2); % 使用fsolve求解渐开线函数方程 fm=x/pi*(fai*hd); alfm=atan(k*(1-cos(2*pi*fm/(fai*hd)))/(fai*hd+k*fm-k*fai*hd*sin(2*pi*fm/(fai*hd))/(2*pi))); end fprintf (1,' 最大压力角 alfm = %3.4f 度 \\n',alfm/hd) fprintf (1,' 最大压力角的位置角 fm = %3.4f 度 \\n',fm/hd) % 压力角渐开线函数 function f=TLYLJ(x) global k % 定义全局变量 f=tan(x)-x-pi/k; disp ' ******** 偏置移动从动件盘形凸轮设计 ********' disp '已知条件:' disp ' 凸轮作逆时针方向转动,从动件偏置在凸轮轴心的右边' disp ' 从动件在推程作等加速/等减速运动,在回程作余弦加速度运动' rb = 40;rt = 10;e = 15;h = 50;ft = 100;fs = 60;fh = 90;alp = 35; fprintf (1,' 基圆半径 rb = %3.4f mm \\n',rb) fprintf (1,' 滚子半径 rt = %3.4f mm \\n',rt) fprintf (1,' 推杆偏距 e = %3.4f mm \\n',e) fprintf (1,' 推程升程 h = %3.4f mm \\n',h) fprintf (1,' 推程运动角 ft = %3.4f 度 \\n',ft) fprintf (1,' 远休止角 fs = %3.4f 度 \\n',fs) fprintf (1,' 回程运动角 fh = %3.4f 度 \\n',fh) fprintf (1,' 推程许用压力角 alp = %3.4f 度 \\n',alp) hd= pi / 180;du = 180 / pi; se=sqrt( rb^2 - e^2 ); d1 = ft + fs;d2 = ft + fs + fh; disp ' ' disp '计算过程和输出结果:' disp ' 1- 计算凸轮理论轮廓的压力角和曲率半径' disp ' 1-1 推程(等加速/等减速运动)' s = zeros(ft);ds = zeros(ft);d2s = zeros(ft); at = zeros(ft);atd = zeros(ft);pt = zeros(ft); for f = 1 : ft if f <= ft / 2 s(f)=2 * h * f ^ 2 / ft ^ 2;s = s(f); ds(f)=4 * h * f * hd / (ft * hd) ^ 2;ds = ds(f); d2s(f)=4 * h / (ft * hd) ^ 2;d2s = d2s(f); else s(f)=h - 2 * h * (ft - f) ^ 2 / ft ^ 2;s = s(f); ds(f)=4 * h * (ft - f) * hd / (ft * hd) ^ 2;ds = ds(f); d2s(f)=-4 * h / (ft * hd) ^ 2;d2s = d2s(f); end at(f)= atan(abs(ds - e) / (se + s));atd(f) = at(f) * du; p1=((se + s ) ^ 2 + (ds - e) ^ 2) ^ 1.5; p2= abs((se + s) * (d2s - se - s) - (ds - e) * (2 * ds - e)); pt(f)= p1 /p2;p = pt(f); end atm = 0; for f = 1 : ft if atd(f) > atm atm = atd(f); end end fprintf (1,' 最大压力角 atm = %3.4f 度\\n',atm) for f = 1 : ft if abs(atd(f) - atm) < 0.1 ftm = f;break end end fprintf (1,' 对应的位置角 ftm = %3.4f 度\\n',ftm) if atm > alp fprintf (1,' * 凸轮推程压力角超过许用值,需要增大基圆!\\n') end ptn = rb + h; for f = 1 : ft if pt(f) < ptn ptn = pt(f); end end fprintf (1,' 轮廓最小曲率半径 ptn = %3.4f mm\\n',ptn) for f = 1 : ft if abs(pt(f) - ptn) < 0.1 ftn = f;break end end fprintf (1,' 对应的位置角 ftn = %3.4f 度\\n',ftn) if ptn < rt + 5 fprintf (1,' * 凸轮推程轮廓曲率半径小于许用值,需要增大基圆或减小滚子!\\n') end disp ' 1-2 回程(余弦加速度运动)' s = zeros(fh);ds = zeros(fh);d2s = zeros(fh); ah = zeros(fh);ahd = zeros(fh);ph = zeros(fh); for f = d1 : d2 k = f - d1; s(f) = .5 * h * (1 + cos(pi * k / fh)); s = s(f); ds(f)=-.5 * pi * h * sin(pi * k / fh) / (fh * hd);ds = ds(f); d2s(f)= -.5 * pi ^2 * h * cos(pi * k / fh)/(fh * hd) ^2;d2s = d2s(f); ah(f)=atan(abs(ds + e) / (se + s));ahd(f) = ah(f) * du; p1=((se + s ) ^ 2 + (ds - e) ^ 2) ^ 1.5; p2= abs((se + s) * (d2s - se - s) - (ds - e) * (2 * ds - e)); ph(f)= p1 /p2;p = ph(f); end ahm = 0; for f = d1 : d2 if ahd(f) > ahm; ahm = ahd(f); end end fprintf (1,' 最大压力角 ahm = %3.4f 度\\n',ahm) for f = d1 : d2 if abs(ahd(f)- ahm) < 0.1 fhm = f;break end end fprintf (1,' 对应的位置角 fhm = %3.4f 度\\n',fhm) phn = rb + h; for f = d1 : d2 if ph(f) < phn phn = ph(f); end end fprintf (1,' 轮廓最小曲率半径 phn = %3.4f mm\\n',phn) for f = d1 : d2 if abs(ph(f) - phn) < 0.1 fhn = f;break end end fprintf (1,' 对应的位置角 fhn = %3.4f 度\\n',fhn) if phn < rt + 5 fprintf (1,' * 凸轮回程轮廓曲率半径小于许用值,需要增大基圆或减小滚子!\\n') end disp ' 2- 计算凸轮理论廓线与实际廓线的直角坐标' n = 360; s = zeros(n);ds = zeros(n);r = zeros(n);rp = zeros(n); x = zeros(n);y = zeros(n);dx = zeros(n);dy = zeros(n); xx = zeros(n);yy = zeros(n);xp = zeros(n);yp = zeros(n); xxp = zeros(n);yyp = zeros(n); for f = 1 : n if f <= ft/2 s(f) = 2 * h * f ^ 2 / ft ^ 2; s = s(f); ds(f) = 4 * h * f * hd / (ft * hd) ^ 2; ds = ds(f); elseif f > ft/2 & f <= ft s(f) = h - 2 * h * (ft - f) ^ 2 / ft ^ 2; s = s(f); ds(f) = 4 * h * (ft - f) * hd / (ft * hd) ^ 2; ds = ds(f); elseif f > ft & f <= d1 s = h;ds = 0; elseif f > d1 & f <= d2 k = f - d1; s(f) = .5 * h * (1 + cos(pi * k / fh)); s = s(f); ds(f)= -.5 * pi * h * sin(pi * k / fh) / (fh * hd); ds = ds(f); elseif f > d2 & f <= n s = 0;ds = 0; end xx(f) = (se + s) * sin(f * hd) + e * cos(f * hd); x = xx(f); yy(f) = (se + s) * cos(f * hd) - e * sin(f * hd); y = yy(f); dx(f) = (ds - e) * sin(f * hd) + (se + s) * cos(f * hd); dx = dx(f); dy(f) = (ds - e) * cos(f * hd) - (se + s) * sin(f * hd); dy = dy(f); xp(f) = x + rt * dy / sqrt(dx ^ 2 + dy ^ 2);xxp = xp(f); yp(f) = y - rt * dx / sqrt(dx ^ 2 + dy ^ 2);yyp = yp(f); r(f) = sqrt (x ^2 + y ^2 ); rp(f) = sqrt (xxp ^2 + yyp ^2 ); end disp ' 2-1 推程(等加速/等减速运动)' disp ' 凸轮转角 理论x 理论y 实际x 实际y' for f = 10 : 10 :ft nu = [f xx(f) yy(f) xp(f) yp(f)]; disp(nu) end disp ' 2-2 回程(余弦加速度运动)' disp ' 凸轮转角 理论x 理论y 实际x 实际y' for f = d1 : 10 : d2 nu = [f xx(f) yy(f) xp(f) yp(f)]; disp(nu) end disp ' 2-3 凸轮轮廓向径' disp ' 凸轮转角 理论r 实际r' for f = 10 : 10 : n nu = [f r(f) rp(f)]; disp(nu) end disp '绘制凸轮的理论轮廓和实际轮廓:' plot(xx,yy,'r-.') % 理论轮廓(红色,点划线) axis ([-(rb+h-10) (rb+h+10) -(rb+h+10) (rb+rt+10)]) % 横轴和纵轴的下限和上限 axis equal % 横轴和纵轴的尺度比例相同 text(rb+h+3,0,'X') % 标注横轴 text(0,rb+rt+3,'Y') % 标注纵轴 text(-5,5,'O') % 标注直角坐标系原点 title('偏置移动从动件盘形凸轮设计') % 标注图形标题 hold on; % 保持图形 plot([-(rb+h) (rb+h)],[0 0],'k') % 横轴(黑色) plot([0 0],[-(rb+h) (rb+rt)],'k') % 纵轴(黑色) plot([e e],[0 (rb+rt)],'k--') % 初始偏置位置(黑色,虚线) ct = linspace(0,2*pi); % 画圆的极角变化范围 plot(rb*cos(ct),rb*sin(ct),'g') % 基圆(绿色) plot(e*cos(ct),e*sin(ct),'c--') % 偏距圆(蓝绿色,虚线) plot(e + rt*cos(ct),se + rt*sin(ct),'y') % 滚子圆(黄色) plot(xp,yp,'b') % 实际轮廓(蓝色) % 旋轮线轨迹模拟 % 圆锥齿轮传动参数 m=2;z7=11;z8=36; dt7=10.5;dt8=35; bt=65.25; % 行星轮轴线与XOY平面夹角 % 搅拌杆外点旋转半径 l=65; hd=pi/180; % 圆锥齿轮几何尺寸 r7=0.5*m*z7; % r7v=r7/cos(dt7*hd); r8=0.5*m*z8; % r8v=r8/cos(dt8*hd); k1=r7+r8; k2=(r7+r8)/r7; % 1----二维旋轮线参数计算 for i=1:1:360 x=k1*cos(i*hd)+l*cos(k2*i*hd); y=k1*sin(i*hd)+l*sin(k2*i*hd); xl1(i,:)=[i x y]; end for i=1:1:720 x=k1*cos(i*hd)+l*cos(k2*i*hd); y=k1*sin(i*hd)+l*sin(k2*i*hd); xl2(i,:)=[i x y]; end for i=1:1:1080 x=k1*cos(i*hd)+l*cos(k2*i*hd); y=k1*sin(i*hd)+l*sin(k2*i*hd); xl3(i,:)=[i x y]; end for i=1:1:1440 x=k1*cos(i*hd)+l*cos(k2*i*hd); y=k1*sin(i*hd)+l*sin(k2*i*hd); xl4(i,:)=[i x y]; end % 绘制二维旋轮线 figure(1); % 生成第1个图形窗口 subplot(2,2,1); % 选择第1个子窗口 plot(xl1(:,2),xl1(:,3)) grid % 绘制网格线 title('二维旋轮线(\heta =360度)') subplot(2,2,2); % 选择第2个子窗口 plot(xl2(:,2),xl2(:,3)) grid % 绘制网格线 title('二维旋轮线(\heta =720度)') subplot(2,2,3); % 选择第3个子窗口 plot(xl3(:,2),xl3(:,3)) grid % 绘制网格线 title('二维旋轮线(\heta =1080度)') subplot(2,2,4); % 选择第4个子窗口 plot(xl4(:,2),xl4(:,3)) grid % 绘制网格线 title('二维旋轮线(\heta =1440度)') % 2----三维旋轮线参数计算 for i=1:1:360 x=(k1*cos(i*hd)+l*cos(k2*i*hd))*cos(bt*hd); y=(k1*sin(i*hd)+l*sin(k2*i*hd))*cos(bt*hd); z=sqrt(x^2+y^2)*tan(bt*hd); xlx1(i,:)=[i x y z]; end for i=1:1:720 x=(k1*cos(i*hd)+l*cos(k2*i*hd))*cos(bt*hd); y=(k1*sin(i*hd)+l*sin(k2*i*hd))*cos(bt*hd); z=sqrt(x^2+y^2)*tan(bt*hd); xlx2(i,:)=[i x y z]; end for i=1:1:1080 x=(k1*cos(i*hd)+l*cos(k2*i*hd))*cos(bt*hd); y=(k1*sin(i*hd)+l*sin(k2*i*hd))*cos(bt*hd); z=sqrt(x^2+y^2)*tan(bt*hd); xlx3(i,:)=[i x y z]; end for i=1:1:1440 x=(k1*cos(i*hd)+l*cos(k2*i*hd))*cos(bt*hd); y=(k1*sin(i*hd)+l*sin(k2*i*hd))*cos(bt*hd); z=sqrt(x^2+y^2)*tan(bt*hd); xlx4(i,:)=[i x y z]; end % 矩阵cs(:,j)表示第j列的各行元素 % 绘制三维旋轮线 figure(2); % 生成第2个图形窗口 subplot(2,2,1); % 选择第1个子窗口 plot3(xlx1(:,2),xlx1(:,3),xlx1(:,4)) grid title('三维旋轮线(\heta =360度)') subplot(2,2,2); % 选择第2个子窗口 plot3(xlx2(:,2),xlx2(:,3),xlx2(:,4)) grid title('三维旋轮线(\heta =720度)') subplot(2,2,3); % 选择第3个子窗口 plot3(xlx3(:,2),xlx3(:,3),xlx3(:,4)) grid title('三维旋轮线(\heta =1080度)') subplot(2,2,4); % 选择第4个子窗口 plot3(xlx4(:,2),xlx4(:,3),xlx4(:,4)) grid title('三维旋轮线(\heta =1440度)') % 输出数据 disp ' ' disp ' ========== 圆锥齿轮传动参数=========='; fprintf(1,' 行星轮齿数 z7 = %3.0f \\n',z7); fprintf(1,' 分度圆半径 r7 = %3.3f 度 \\n',r7); fprintf(1,' 分度圆锥角 dt7 = %3.3f 度 \\n',dt7); fprintf(1,' 中心轮齿数 z8 = %3.0f \\n',z8); fprintf(1,' 分度圆半径 r8 = %3.3f 度 \\n',r8); fprintf(1,' 分度圆锥角 dt8 = %3.3f 度 \\n',dt8); fprintf(1,' 中心轮与行星轮当量半径之比 rb = %3.0f \\n',r8/r7); disp ' ' [' 行星转',' 动点x坐标',' 动点y坐标',' 动点z坐标'] [xlx1(:,1),xlx1(:,2),xlx1(:,3),xlx1(:,4)] % 锥齿轮行星机构运动分析 i=37/12*36/10*50/14; % 前置圆柱轮系传动比 hd=pi/180;d6=14.5*hd;d7=10.25*hd; % 锥齿轮节锥角 n1=[7755 11221 16021]; % 电动机转速(第1、5、10档) for j=1:3 n6=n1(1:j)/i; % 中心锥齿轮6的绝对转速 n7=n6*sin(d6)/sin(2*d7); % 行星锥齿轮7的绝对转速 nh=n7*sin(d7)/sin(d6+d7); % 行星锥齿轮7的公转速度 n7h=n7*sin(d6+2*d7)/sin(d6+d7); % 行星锥齿轮7的自转速度 end disp ' ======== 已知条件 ========' fprintf(1,' 前置圆柱轮系传动比 i = %3.3f \\n',i); fprintf(1,' 锥齿轮节锥角 d7 = %3.3f 度 \\n',d7/hd); disp ' ======== 输出计算结果 ========' for j=1:3 fprintf(1,' 电动机转速 n1 = %3.0f r/min \\n',n1(j)); fprintf(1,' 中心锥齿轮转速 n6 = %3.3f r/min \\n',n6(j)); fprintf(1,' 行星锥齿轮转速 n7 = %3.3f r/min \\n',n7(j)); fprintf(1,' 行星锥齿转转速 nh = %3.3f r/min \\n',nh(j)); fprintf(1,' 行星锥齿轮自转转速 n7h = %3.3f r/min \\n',n7h(j)); disp ' ' end disp ' ******** 圆柱齿轮传动设计 ********' P1 = 22 ; N1 = 970 ; I = 3.5 ; K = 1.4; hd = pi/180 ; Ha = 1.0 ; Ca = 0.25; m = [2 2.5 3 4 5 6 8 10 12 16 20 25 32 40 50]; % 优先选用第一系列 disp '(注意：以下输入的齿轮材料和齿面硬度类别的标识字符要用大写,并用单引号括起)' CL = input(' 选择齿轮材料：碳钢-TG;合金钢-HG == '); CM = input(' 选择齿面硬度类别：硬齿面-YC;软齿面-RC == '); disp ' 齿宽系数的选择参考：' disp ' 对称布置--软齿面 0.8-1.4;硬齿面 0.4-0.9' disp ' 非对称布置--软齿面 0.6-1.2;硬齿面 0.3-0.6' disp ' 悬臂布置--软齿面 0.3-0.4;硬齿面 0.2-0.25' Fd = input(' 选择齿宽系数： Fd = '); disp ' ====== 已知条件 ======' fprintf (1,' 主动轮传递功率 P1 = %3.3f Kw \\n',P1) fprintf (1,' 主动轮转速 n1 = %3.3f r/min \\n',N1) fprintf (1,' 传动比 i = %3.3f \\n',I) fprintf (1,' 载荷系数 K = %3.3f \\n',K) fprintf (1,' 齿高系数 ha* = %3.3f \\n',Ha) fprintf (1,' 顶隙系数 c* = %3.3f \\n',Ca) fprintf (1,' 齿宽系数 Fd = %3.3f \\n',Fd) disp '(注意：以下输入齿面热处理硬度是数值)' if CL == 'TG' disp ' 齿轮材料-碳钢' if CM == 'YC' disp ' 齿面硬度类别-硬齿面' HRC1 = input(' 输入小齿轮感应淬火硬度HRC1 = '); HRC2 = input(' 大齿轮感应淬火硬度HRC2 = '); CHM1 = 11*HRC1+610 ; CHM2 = 11*HRC2+610; CFM1 = 380 ; CFM2 = 380; elseif CM == 'RC' disp ' 齿面硬度类别-软齿面' HBS1 = input(' 输入小齿轮调质/正火硬度HBS1 = '); HBS2 = input(' 大齿轮调质/正火硬度HBS2 = '); CHM1 = (9*HBS1+3750)/10 ; CHM2 = (9*HBS2+3750)/10; CFM1 = (2*HBS1+675)/5 ; CFM2 = (2*HBS2+675)/5; end elseif CL == 'HG' disp ' 齿轮材料-合金钢' if CM == 'YC' disp ' 齿面硬度类别-硬齿面' CHM1 = 1500 ; CHM2 = 1500; CFM1 = 460 ; CFM2 = 460; elseif CM == 'RC' disp ' 齿面硬度类别-软齿面' HBS1 = input(' 输入小齿轮调质/正火硬度HBS1 = '); HBS2 = input(' 大齿轮调质/正火硬度HBS2 = '); CHM1 = (10*HBS1+1700)/6 ; CHM2 = (10*HBS2+1700)/6; CFM1 = (19*HBS1+9700)/50; CFM2 = (19*HBS2+9700)/50; end end disp '(注意：以下输入的齿轮传动方向类别的标识字符要用大写,并用单引号括起)' FX = input(' 选择齿轮传动方向：单向传动-DX;双向传动-SX == '); if FX == 'DX' disp ' 齿轮单向传动' CFP1 = 1.4*CFM1 ; CFP2 = 1.4*CFM2; elseif FX == 'SX' disp ' 齿轮双向传动' CFP1 = CFM1 ; CFP2 = CFM2; end fprintf (1,' 小齿轮齿根弯曲许用应力 CFP1 = %3.3f MPa \\n',CFP1) fprintf (1,' 大齿轮齿根弯曲许用应力 CFP2 = %3.3f MPa \\n',CFP2) disp '(注意：以下输入齿轮齿数是数值)' if CM == 'YC' disp ' 闭式硬齿面齿轮传动小齿轮齿数范围是：17-20' Z1 = input(' 输入小齿轮齿数 z1 = '); elseif CM == 'RC' disp ' 闭式软齿面齿轮传动小齿轮齿数范围是：20-30' Z1 = input(' 输入小齿轮齿数 z1 = '); end Z2 = round(I*Z1); fprintf (1,' 大齿轮齿数 z2 = %3.0f \\n',Z2) U= Z2/Z1; disp ' 初选斜齿轮螺旋角的范围是9-10度,直齿轮是0度' disp ' 经过试算,如果螺旋角超过15度,重新初选的范围是16-17度' disp ' 经过试算,如果螺旋角超过25度,重新初选的范围是26-27度' BAT = input(' 输入齿轮螺旋角(度) Bat = '); BATR = BAT*hd; ZV1 = Z1/(cos(BATR))^3 ; ZV2 = U*ZV1; YFS1 = ZV1/(0.269118*ZV1-0.840687); YFS2 = ZV2/(0.269118*ZV2-0.840687); fprintf (1,' 小齿轮齿形系数 Yfs1 = %3.3f \\n',YFS1) fprintf (1,' 大齿轮齿形系数 Yfs2 = %3.3f \\n',YFS2) if YFS1/CFP1 > YFS2/CFP2 YFCP = YFS1/CFP1; else YFCP = YFS2/CFP2; end CHP1= 0.9*CHM1 ; CHP2 = 0.9*CHM2; fprintf (1,' 小齿轮齿面接触许用应力 CHP1 = %3.3f MPa \\n',CHP1) fprintf (1,' 大齿轮齿面接触许用应力 CHP2 = %3.3f MPa \\n',CHP2) if CHP1 > CHP2 CHP = CHP2; else CHP = CHP1; end [Aa,Ad,Am] = ADM(BAT); % 根据螺旋角查询齿轮强度计算系数 T1= 9550*P1/N1; fprintf (1,' 主动轮传递的转矩 T1 = %3.3f NM \\n',T1) if CM == 'YC' disp ' 硬齿面齿轮传动-按照齿根弯曲强度确定齿轮模数' Mj = Am*(K*T1*YFCP/(Fd*Z1^2))^(1/3); for i = 1 : 15 if Mj <= m(i) Mn = m(i) ; break end end Aj = Mn*(Z1+Z2)/(2*cos(BATR)); A = round(Aj/2+0.5)*2; elseif CM == 'RC' disp ' 软齿面齿轮传动-按照齿面接触强度确定齿轮直径' Dj = Ad*(K*T1*(U+1)/(Fd*U*CHP^2))^(1/3); Aj = Dj*(1+U)/(2*cos(BATR)); A = round(Aj/2+0.5)*2; Mj = 2*A*cos(BATR)/(Z1+Z2); for i = 1 : 15 if Mj <= m(i) Mn = m(i) ; break end end end fprintf (1,' ( 中心距计算值 Aj = %3.3f mm )\\n',Aj) fprintf (1,' ( 中心距圆整值 A = %3.3f mm )\\n',A) if BAT < 8 BATA = 0; else BA = Mn*(Z1+Z2)/(2*A); fprintf (1,' ( 螺旋角余弦值 Ba = %3.3f )\\n ',BA) if BA >= 1 disp ' 无法计算螺旋角,需要增大中心距!' DT = round((Mn*(Z1+Z2)-2*A)/2)+2; A= A+DT; BAA = Mn*(Z1+Z2)/(2*A); fprintf (1,' * 增大齿轮传动中心距后的螺旋角余弦值 = %3.3f \\n',BAA) end BATj = acos(Mn*(Z1+Z2)/(2*A)); BATA = BATj/hd; end D1= Mn*Z1/(cos(BATj)); [Aa,Ad,Am] = ADM(BATA); % 根据螺旋角查询齿轮强度计算系数 if CM == 'YC' disp ' * 按照齿面接触强度校核 !' Dx = Ad*(K*T1*(U+1)/(Fd*U*CHP^2))^(1/3); if Dx <= D1 disp ' 满足齿面接触强度!' else fprintf (1,' ( 按照齿面接触强度需要的小齿轮分度圆直径 Dx = %3.3f mm )\\n',Dx) disp ' 不满足齿面接触强度,需要重新计算!' Aj = Dx*(1+U)/(2*cos(BATR)); A = round(Aj/2+0.5)*2+2; Mj = 2*A*cos(BATR)/(Z1+Z2); for i = 1 : 15 if Mj <= m(i) Mn = m(i) ; break end end BA = Mn*(Z1+Z2)/(2*A); fprintf (1,' ( 螺旋角余弦值 Ba = %3.3f )\\n ',BA) if BA >= 1 disp ' 无法计算螺旋角,需要增大中心距!' DT = round((Mn*(Z1+Z2)-2*A)/2)+2; A = A+DT; BAA = Mn*(Z1+Z2)/(2*A); fprintf (1,' * 增大齿轮传动中心距后的螺旋角余弦值 = %3.3f \\n',BAA) end BATj = acos(Mn*(Z1+Z2)/(2*A)); BATA = BATj/hd; end elseif CM == 'RC' disp ' * 按照齿根弯曲强度校核 !' ZV1 = Z1/(cos(BATj))^3 ; ZV2 = U*ZV1; YFS1 = ZV1/(0.269118*ZV1-0.840687); YFS2 = ZV2/(0.269118*ZV2-0.840687); if YFS1/CFP1 > YFS2/CFP2 YFCP = YFS1/CFP1; else YFCP = YFS2/CFP2; end Mx = Am*(K*T1*YFCP/(Fd*Z1^2))^(1/3); if Mx <= Mj disp ' 满足齿根弯曲强度!' else fprintf (1,' ( 按照齿根弯曲强度需要的齿轮模数 Mx = %3.3f mm \\)n',Mx) disp ' 不满足齿根弯曲强度,需要重新计算!' for i = 1 : 15 if Mx <= m(i) Mn = m(i) ; break end end Aj = Mn*(Z1+Z2)/(2*cos(BATR)); A = round(Aj/2+0.5)*2; BATj = acos(Mn*(Z1+Z2)/(2*A)); BATA = BATj/hd; end end D1 = Mn*Z1/(cos(BATj)); D2 = Mn*Z2/(cos(BATj)); B = Fd*D1 ; B2 = 2*round(B/2+0.5) ; B1 = B2+6; Da1 = D1+2*Ha*Mn; Da2 = D2+2*Ha*Mn; disp ' ====== 输出齿轮传动参数的设计结果 ======' fprintf (1,' 小齿轮齿数 Z1 = %3.0f \\n',Z1) fprintf (1,' 大齿轮齿数 Z2 = %3.0f \\n',Z2) fprintf (1,' 齿轮副模数 Mn = %3.2f mm \\n',Mn) fprintf (1,' 齿轮副螺旋角 Bata = %3.3f 度 \\n',BATA) fprintf (1,' 齿轮副中心距 a = %3.3f mm \\n',A) fprintf (1,' 小齿轮分度圆直径 d1 = %3.3f mm \\n',D1) fprintf (1,' 大齿轮分度圆直径 d2 = %3.3f mm \\n',D2) fprintf (1,' 小齿轮齿顶圆直径 da1 = %3.3f mm \\n',Da1) fprintf (1,' 大齿轮齿顶圆直径 da2 = %3.3f mm \\n',Da2) fprintf (1,' 小齿轮宽度 b1 = %3.0f mm \\n',B1) fprintf (1,' 大齿轮宽度 b2 = %3.0f mm \\n',B2) % 齿轮传动计算系数函数文件ADM.m function [Aa,Ad,Am]=ADM(BAT) if BAT < 8 Aa = 483 ; Ad = 766 ; Am = 12.6; elseif BAT >= 8 & BAT < 15 Aa = 476 ; Ad = 756 ; Am = 12.4; elseif BAT >= 15 & BAT < 25 Aa = 462 ; Ad = 733 ; Am = 12.0; elseif BAT >= 25 & BAT < 35 Aa = 447 ; Ad = 709 ; Am = 11.5; End % 变位直齿圆柱齿轮参数测定 z=8; % 齿数 df0=33.43; % 齿根圆直径的测量值 % 变位齿法线长度的测量值 Wk=24.73;Wk1=39.43; % 跨齿数 k=round(z/9+0.5); if k<2 k=2; end Pb=Wk1-Wk; % 基圆齿距 alf=20;hd=pi/180; % 压力角 m=round(Pb/(pi*cos(alf*hd))); % 模数 Wkb=m*cos(alf*hd)*((k-0.5)*pi+z*0.0149044); % 标准齿法线长度 x1=(Wk-Wkb)/(2*m*sin(alf*hd)); % 变位系数 hf=(m*z-df0)/2; % 齿根高 % 齿顶高系数与顶隙系数 hc=hf/m+x1; disp ' ' fprintf(1,' 齿顶高系数与顶隙系数之和 hc = %3.2f \\n',hc); hx=1.00;cx=0.25; % 按照hc计算值确定齿制-正常齿或短齿 % 输出齿轮参数 disp ' ' disp ' ========== 变位齿轮齿轮参数 =========='; fprintf(1,' 齿数 z = %3.0f \\n',z); fprintf(1,' 压力角 alf = %3.0f 度 \\n',alf); fprintf(1,' 模数 m = %3.3f mm \\n',m); fprintf(1,' 齿顶高系数 hx = %3.2f \\n',hx); fprintf(1,' 顶隙系数 cx = %3.2f \\n',cx); fprintf(1,' 变位系数 x = %3.3f \\n',x1); disp ' ' disp ' ========== 变位齿轮测量和计算数据 ==========' fprintf(1,' 跨齿数 k = %3.0f \\n',k); fprintf(1,' 测量齿根圆直径 df0 = %3.3f mm \\n',df0); fprintf(1,' 齿根高 hf = %3.3f mm \\n',hf); fprintf(1,' 基圆齿距 Pb = %3.3f mm \\n',Pb); fprintf(1,' 测量齿法线长度 Wk = %3.3f mm \\n',Wk); fprintf(1,' 标准齿法线长度 Wkb = %3.3f mm \\n',Wkb); % 计算啮合角 Qp=2*(x1+x1)*tan(alf*hd)/(z+z)+0.0149044; % 节圆处展角弧度值 [x,f]=fsolve('tan(x)-x-0.0688793',0.0149044); % 使用fsolve求解渐开线函数方程 alfp=x/hd; % 啮合角 disp ' ' disp ' ========== 齿轮副啮合角和渐开线函数值 =========='; fprintf(1,' 啮合角 alfp = %3.3f 度 \\n',alfp); fprintf(1,' 啮合角渐开线函数值 Qp = %3.7f \\n',Qp); % 计算中心距、分离系数、齿顶变动系数与几何尺寸 a=0.5*m*(z+z); % 标准中心距 ap=a*cos(alf*hd)/cos(alfp*hd); % 实际中心距 y=(ap-a)/m; % 分离系数 sgm=x1+x1-y; % 齿顶变动系数 d=m*z; % 分度圆直径 db=d*cos(alf*hd); % 基圆直径 da=d+2*(hx+x1-sgm)*m; % 齿顶圆直径 df=d-2*(hx+cx-x1)*m; % 齿根圆直径 Wkp=Wkb+2*x1*m*sin(alf*hd); % 公法线长度 % 计算变位齿轮齿厚 alfa=acos(db/da); % 齿顶压力角 s=pi*m/2+2*x1*m*tan(alf*hd); % 分度圆齿厚 sa=s*da/d-da*(tan(alfa)-alfa-0.0149044); % 齿顶圆齿厚 sb=cos(alf*hd)*(s+d*0.0149044); % 基圆齿厚 disp ' ' disp ' ========== 变位齿轮齿厚和啮合角 =========='; fprintf(1,' 分度圆齿厚 s = %3.3f mm \\n',s); fprintf(1,' 齿顶圆齿厚 sa = %3.3f mm \\n',sa); fprintf(1,' 基圆齿厚 sb = %3.3f mm \\n',sb); fprintf(1,' 齿顶压力角 alfa = %3.3f 度 \\n',alfa/hd); fprintf(1,' 啮合角 alfp = %3.3f 度 \\n',alfp); disp ' ' disp ' ========== 变位齿轮参数和几何尺寸 =========='; fprintf(1,' 中心距分离系数 y = %3.3f \\n',y); fprintf(1,' 齿顶变动系数 sgm = %3.3f \\n',sgm); fprintf(1,' 标准中心距 a = %3.3f mm \\n',a); fprintf(1,' 实际中心距 ap = %3.3f mm \\n',ap); fprintf(1,' 齿顶圆直径 da = %3.3f mm \\n',da); fprintf(1,' 分度圆直径 d = %3.3f mm \\n',d); fprintf(1,' 基圆直径 db = %3.3f mm \\n',db); fprintf(1,' 齿根圆直径 df = %3.3f mm \\n',df); fprintf(1,' 公法线长度 Wkp = %3.3f mm \\n',Wkp); % 根据基圆齿厚、模数和压力角计算变位系数 x2=(sb/(m*cos(alf*hd))-0.5*pi-0.0149044*z)/(2*tan(alf*hd)); fprintf(1,' 变位系数 x = %3.3f \\n',x2); % 计算斜齿法线长度及其偏差 % 输入齿轮基本参数 Mn=4;z=60;an=20;bat=12.92; disp ' ' disp ' ========== 斜齿圆柱齿轮基本参数 ==========' fprintf (1,' 模数 Mn = %3.2f mm \\n',Mn) fprintf (1,' 齿数 z = %3.0f \\n',z) fprintf (1,' 压力角 an = %3.3f 度 \\n',an) fprintf (1,' 螺旋角 bat = %3.3f 度 \\n',bat) % 角度转换为弧度 hd=pi/180; anh=an*hd; bath=bat*hd; invan=tan(anh)-anh; % 计算跨齿数和公法线长度 ath=atan(tan(anh)/cos(bath)); invan=tan(anh)-anh; invat=tan(ath)-ath; zp=z*invat/invan; k=round(zp/9+0.5); Wkn=Mn*cos(anh)*(pi*(k-0.5)+zp*invan); disp ' ' disp ' ========== 计算斜齿圆柱齿法线长度及其偏差 ==========' fprintf (1,' 端面压力角 at = %3.3f 度 \\n',ath/hd) fprintf (1,' 相当齿数 zp = %3.3f \\n',zp) fprintf (1,' 跨齿数 k = %2.0f \\n',k) fprintf (1,' 公法线长度 Wkn = %3.3f mm \\n',Wkn) % 输入齿距极限偏差和齿圈径向跳动公差 fpt=0.028;Fr=0.071; % 输入齿厚极限偏差代号 H=-8;L=-16; % 计算齿厚极限偏差 Es=H*fpt;Ei=L*fpt; fprintf (1,' 齿厚上偏差 Es = %3.3f mm \\n',Es) fprintf (1,' 齿厚下偏差 Ei = %3.3f mm \\n',Ei) % 计算公法线长度极限偏差 Ews=Es*cos(anh)-0.72*Fr*sin(anh); Ewi=Ei*cos(anh)+0.72*Fr*sin(anh); fprintf (1,' 公法线长度上偏差 Ews = %3.3f mm \\n',Ews) fprintf (1,' 公法线长度下偏差 Ewi = %3.3f mm \\n',Ewi) % 轴系设计计算(斜齿圆柱齿轮传动设计计算、轴弯扭组合强度计算、滚动轴承寿命计算) % ---- 已知条件 ---- f=2000; % 输送带工作拉力(N) v=1.5; % 输送带工作速度(m/s) d=250; % 滚筒直径(mm) nu=0.97; % 斜齿圆柱齿轮传动效率 i=3.5; % 齿轮副传动比 hd=pi/180; % 角度换算成弧度的系数 disp ' ' disp ' ========== 已知条件 =========='; fprintf(1,' 输送带工作拉力 f = %3.0f N \\n',f); fprintf(1,' 输送带工作速度 v = %3.2f m/s \\n',v); fprintf(1,' 滚筒直径 d = %3.0f mm \\n',d); fprintf(1,' 齿轮传动效率 nu = %3.2f \\n',nu); fprintf(1,' 齿轮副传动比 i = %3.2f \\n',i); % ---- 1-齿轮传动设计计算(采用硬齿面齿轮传动) ---- p2=f*v/1000; % 大齿轮传递功率(kW) n2=60*v*1e3/pi/d; % 大齿轮转速(r/min) p1=p2/nu; % 小齿轮传递功率(kW) n1=i*n2; % 小齿轮转速(r/min) chm=1500; % 试验齿轮接触疲劳极限(MPa) cfm=460; % 试验齿轮弯曲疲劳极限(MPa) chp=0.9*chm; % 试验齿轮许用接触应力(MPa) cfp=1.4*cfm; % 试验齿轮许用弯曲应力(单向传动) z1=18; % 小齿轮齿数(选取) z2=round(i*z1); % 大齿轮齿数 u=z2/z1; % 齿数比 pd=0.675; % 齿宽系数 bat0=10; % 螺旋角初值 t1=9550*p1/n1; % 小齿轮传递转矩(Nm) zv1=z1/(cos(bat0*hd))^3; % 小齿轮当量齿数 zv2=u*zv1; % 大齿轮当量齿数 ysf1=4.43; % 小齿轮齿形系数 ysf2=3.88; % 大齿轮齿形系数 if ysf1>=ysf2 ysf=ysf1; % 确定计算齿形系数 else ysf=ysf2; end k=1.6; % 载荷系数 am=12.0; % 齿根弯曲强度计算系数(螺旋角范围15-25度) mnj=am*(k*t1*ysf/pd/z1^2/cfp)^(1/3); % 按照齿根弯曲强度计算模数(mm) if mnj<=2 mn=2; % 确定标准模数(mm) else mn=round(mnj+0.5) end aj=mn*z1*(1+u)/2/cos(bat0*hd); a=round(aj/5)*5+5; % 确定中心距(mm) bat=acos(0.5*mn*z1*(1+u)/a)/hd; % 确定螺旋角 disp ' ' disp ' ========== 齿轮传动设计计算 =========='; fprintf(1,' 大齿轮传递功率 p2 = %3.3f kW \\n',p2); fprintf(1,' 大齿轮转速 n2 = %3.3f r/min \\n',n2); fprintf(1,' 小齿轮传递功率 p1 = %3.3f kW \\n',p1); fprintf(1,' 小齿轮转速 n1 = %3.3f r/min \\n',n1); fprintf(1,' 试验齿轮许用接触应力 chp = %3.3f MPa \\n',chp); fprintf(1,' 试验齿轮许用弯曲应力 cfp = %3.3f MPa \\n',cfp); fprintf(1,' 小齿轮齿数 z1 = %3.0f \\n',z1); fprintf(1,' 大齿轮齿数 z2 = %3.0f \\n',z2); fprintf(1,' 齿宽系数 pd = %3.3f \\n',pd); fprintf(1,' 齿数比 u = %3.3f \\n',u); fprintf(1,' 小齿轮传递转距 t1 = %3.3f Nm \\n',t1); fprintf(1,' 小齿轮当量齿数 zv1 = %3.3f \\n',zv1); fprintf(1,' 大齿轮当量齿数 zv2 = %3.3f \\n',zv2); fprintf(1,' 小齿轮齿形系数 ysf1 = %3.3f \\n',ysf1); fprintf(1,' 大齿轮齿形系数 ysf2 = %3.3f \\n',ysf2); fprintf(1,' 齿轮模数 mn = %3.2f mm \\n',mn); fprintf(1,' 中心距 a = %3.2f mm \\n',a); fprintf(1,' 螺旋角 bat = %3.3f 度 \\n',bat); if bat>15 & bat<=25 '螺旋角在15-25度范围内,计算系数选择合适' else '螺旋角超出15-25度范围,重新选择计算系数' end d1=mn*z1/cos(bat*hd); % 计算分度圆直径(mm) d2=u*d1; han=1.0; % 正常齿制 cn=0.25; da1=d1+2*han*mn; % 计算齿顶圆直径(mm) da2=d2+2*han*mn; df1=d1-2*han*mn-2*cn; % 计算齿根圆直径(mm) df2=d2-2*han*mn-2*cn; b=pd*d1; b2=round(b/2)*2; % 确定齿宽(mm) b1=b2+6; ad=733; % 齿面接触强度计算系数(螺旋角范围15-25度) d1j=ad*(k*t1*(u+1)/pd/chp^2/u)^(1/3); % 按照齿面接触强度计算分度圆直径(mm) if d1j<=d1 '满足齿面接触强度要求' else '不满足齿面接触强度要求,需要修改设计参数' end v=pi*d1*n1/6e4; % 齿轮圆周速度(m/s) fprintf(1,' 分度圆直径 d1 = %3.3f mm \\n',d1); fprintf(1,' 小齿轮齿顶圆直径 da1 = %3.3f mm \\n',da1); fprintf(1,' 齿根圆直径 df1 = %3.3f mm \\n',df1); fprintf(1,' 齿宽 b1 = %3.0f mm \\n',b1); fprintf(1,' 大齿轮齿顶圆直径 da2 = %3.3f mm \\n',da2); fprintf(1,' 齿根圆直径 df2 = %3.3f mm \\n',df2); fprintf(1,' 齿宽 b2 = %3.0f mm \\n',b2); % ---- 2-轴的设计计算(弯扭组合) ---- c=112; % 45钢的材料系数 d0=c*(p2/n2)^(1/3)*1.05; % 按扭转估算轴径,并考虑键槽影响 d=round(d0/5)*5; alf=20; % 齿轮分度圆压力角 ft=round(2000*t1/d1); % 齿轮传递的圆周力(N) fr=round(ft*tan(alf*hd)/cos(bat*hd)); % 齿轮传递的径向力(N) fa=round(ft*tan(bat*hd)); % 齿轮传递的轴向力(N) l1=36; % 齿宽中心线到A轴承支座反力作用点的距离(mm) l2=36; % 齿宽中心线到B轴承支座反力作用点的距离(mm) fah=round(ft*l2/(l1+l2)); % A支座H面反力(N) fbh=ft-fah; % B支座H面反力(N) mch=fah*l1; % H面弯矩(Nmm) fav=round((fr*l2+fa*d2/2)/(l1+l2)); % A支座V面反力(N) fbv=fr-fav; % B支座V面反力(N) mcv1=fav*l1; % V面弯矩1(Nmm) mcv2=fbv*l2; % V面弯矩2(Nmm) mc12=mcv1-mcv2; % V面弯矩突变值(Nmm) mcm=round(fa*d2/2); % 在截面C的集中力偶矩(Nmm) mc1=round(sqrt(mch^2+mcv1^2)); % 合成弯矩1(Nmm) mc2=round(sqrt(mch^2+mcv2^2)); % 合成弯矩1(Nmm) if mc1>=mc2 % 确定最大弯矩(Nmm) mc=mc1; else mc=mc2; end t2=round(9.55*1e6*p2/n2); % 大齿轮传递转矩(Nmm) me=round(sqrt(mc^2+(0.6*t2)^2)); % 当量弯矩(Nmm) cp=60; % 对称循环许用弯曲应力(MPa) de=(me/0.1/cp)^(1/3)*1.05; % 按弯扭组合需要轴径,并考虑键槽影响 dc=48; % 危险截面C的实际直径(mm) if de<=dc '满足轴的弯扭组合强度要求' else '不满足轴的弯扭组合强度要求,需要加大轴的直径' end disp ' ========== 轴弯扭组合强度计算 =========='; fprintf(1,' 轴的最小直径 d = %3.3f mm \\n',d); fprintf(1,' 齿轮传递的圆周力 ft = %3.3f N \\n',ft); fprintf(1,' 径向力 fr = %3.3f N \\n',fr); fprintf(1,' 轴向力 fa = %3.3f N \\n',fa); fprintf(1,' H面-A支座反力 fah = %3.3f N \\n',fah); fprintf(1,' B支座反力 fbh = %3.3f N \\n',fbh); fprintf(1,' 弯矩 mch = %3.3f Nmm \\n',mch); fprintf(1,' V面-A支座反力 fav = %3.3f N \\n',fav); fprintf(1,' B支座反力 fbv = %3.3f N \\n',fbv); fprintf(1,' 弯矩1 mcv1 = %3.3f Nmm \\n',mcv1); fprintf(1,' 弯矩2 mcv2 = %3.3f Nmm \\n',mcv2); fprintf(1,' 弯矩突变值 mc12 = %3.3f Nmm \\n',mc12); fprintf(1,' 集中力偶值 mcm = %3.3f Nmm \\n',mcm); fprintf(1,' 合成弯矩1 mc1 = %3.3f Nmm \\n',mc1); fprintf(1,' 合成弯矩2 mc2 = %3.3f Nmm \\n',mc2); fprintf(1,' 大齿轮传递转距 t2 = %3.3f Nm \\n',t2); fprintf(1,' 轴危险截面的当量弯矩 me = %3.2f Nm \\n',me); fprintf(1,' 弯扭组合强度需要的轴径 de = %3.2f mm \\n',de); fprintf(1,' 轴危险截面的实际直径 dc = %3.2f mm \\n',dc); % ---- 3-滚动轴承寿命计算(圆锥滚子轴承30209,正装结构) ---- cr=67800; % 额定动载荷(N) c0r=83500; % 额定静载荷(N) e=0.40; % 判断参数 y=1.5; % 轴向系数 fra=round(sqrt(fah^2+fav^2)); % A轴承径向载荷(N) frb=round(sqrt(fbh^2+fbv^2)); % B轴承径向载荷(N) sa=round(fra/2/y); % A轴承内部轴向力(N) sb=round(frb/2/y); % B轴承内部轴向力(N) if fa+sb>=sa faa=fa+sb; % 确定紧轴承A轴向载荷(N) fab=sb; % 确定松轴承B轴向载荷(N) else fab=abs(fa-sa); % 确定紧轴承B轴向载荷(N) faa=sa; % 确定松轴承A轴向载荷(N) end if faa/fra>=e % 轴承A轴向载荷与径向载荷之比 xa=0.40; % 确定A轴承载荷折算系数X ya=y; % 确定A轴承载荷折算系数Y else xa=1; ya=0; end pa=round(xa*fra+ya*faa); % 轴承A当量动载荷(N) if fab/frb>=e % 轴承B轴向载荷与径向载荷之比 xb=0.40; % 确定B轴承载荷折算系数X yb=y; % 确定B轴承载荷折算系数Y else xb=1; yb=0; end pb=round(xb*frb+yb*fab); % 轴承B当量动载荷(N) fp=1.5; % 轴承载荷系数(减速器中等冲击) lha=round(1e6/60/n2*(cr/fp/pa)^(10/3)); % 计算轴承A寿命(h) lhb=round(1e6/60/n2*(cr/fp/pb)^(10/3)); % 计算轴承B寿命(h) disp ' ' disp ' ========== 圆锥滚子轴承寿命计算 =========='; fprintf(1,' 30209轴承额定动载荷 cr = %3.0f N \\n',cr); fprintf(1,' 额定静载荷 c0r = %3.0f N \\n',c0r); fprintf(1,' 判断参数 e = %3.3f \\n',e); fprintf(1,' 轴向系数 y = %3.3f \\n',y); fprintf(1,' A轴承径向载荷 fra = %3.0f N \\n',fra); fprintf(1,' B轴承径向载荷 frb = %3.0f N \\n',frb); fprintf(1,' A轴承内部轴向力 sa = %3.0f N \\n',sa); fprintf(1,' B轴承内部轴向力 sb = %3.0f N \\n',sb); fprintf(1,' A轴承轴向载荷 faa = %3.0f N \\n',faa); fprintf(1,' B轴承轴向载荷 fab = %3.0f N \\n',fab); fprintf(1,' [A轴承]-轴向与径向载荷之比 fae = %3.3f \\n',faa/fra); fprintf(1,' 径向载荷系数 xa = %3.2f \\n',xa); fprintf(1,' 轴向载荷系数 ya = %3.2f \\n',ya); fprintf(1,' 当量动载荷 pa = %3.0f N \\n',pa); fprintf(1,' [B轴承]-轴向与径向载荷之比 fbe = %3.3f \\n',fab/frb); fprintf(1,' 径向载荷系数 xa = %3.2f \\n',xb); fprintf(1,' 轴向载荷系数 ya = %3.2f \\n',yb); fprintf(1,' 当量动载荷 pa = %3.0f N \\n',pb); fprintf(1,' 轴承A寿命 lha = %3.0f h \\n',lha); fprintf(1,' 轴承B寿命 lha = %3.0f h \\n',lhb); '4-滚动轴承寿命计算(角接触球轴承7009C,正装结构)' c7r=25800; % 额定动载荷(N) c70r=20500; % 额定静载荷(N) s7a=round(0.4*fra); % A轴承内部轴向力(N) s7b=round(0.4*frb); % B轴承内部轴向力(N) if fa+s7b>=s7a fa7a=fa+s7b; % 确定紧轴承A轴向载荷(N) fa7b=s7b; % 确定松轴承B轴向载荷(N) else fa7b=abs(fa-s7a); % 确定紧轴承B轴向载荷(N) fa7a=s7a; % 确定松轴承A轴向载荷(N) end xda=fa7a/c70r; % 轴承A相对轴向载荷 ea=0.465; % A判断参数(插值查表) xdb=fa7b/c70r; % 轴承B相对轴向载荷 eb=0.41; % B判断参数(插值查表) if fa7a/fra>ea % 轴承A轴向载荷与径向载荷之比 x7a=0.44; % 确定A轴承载荷折算系数X y7a=1.21; % 确定A轴承载荷折算系数Y else x7a=1; y7a=0; end p7a=round(x7a*fra+y7a*fa7a); % 轴承A当量动载荷(N) if fa7b/frb>eb % 轴承B轴向载荷与径向载荷之比 x7b=0.44; % 确定B轴承载荷折算系数X y7b=1.37; % 确定B轴承载荷折算系数Y else x7b=1; y7b=0; end p7b=round(x7b*frb+y7b*fa7b); % 轴承B当量动载荷(N) fp=1.5; % 轴承载荷系数(减速器中等冲击) lh7a=round(1e6/60/n2*(c7r/fp/p7a)^3); % 计算轴承A寿命(h) lh7b=round(1e6/60/n2*(c7r/fp/p7b)^3); % 计算轴承B寿命(h) disp ' ' disp ' ========== 角接触球轴承寿命计算 =========='; fprintf(1,' 7009C轴承额定动载荷 c7r = %3.0f N \\n',c7r); fprintf(1,' 额定静载荷 c70r = %3.0f N \\n',c70r); fprintf(1,' A轴承内部轴向力 s7a = %3.0f N \\n',s7a); fprintf(1,' B轴承内部轴向力 s7b = %3.0f N \\n',s7b); fprintf(1,' [A轴承]-轴向载荷 fa7a = %3.0f N \\n',fa7a); fprintf(1,' 相对轴向载荷 xda = %3.3f \\n',xda); fprintf(1,' 判断参数 ea = %3.3f \\n',ea); fprintf(1,' 轴向与径向载荷之比 fae = %3.3f \\n',fa7a/fra); fprintf(1,' 径向载荷系数 x7a = %3.2f \\n',x7a); fprintf(1,' 轴向载荷系数 y7a = %3.2f \\n',y7a); fprintf(1,' 当量动载荷 p7a = %3.0f N \\n',p7a); fprintf(1,' 轴承寿命 lh7a = %3.0f h \\n',lh7a); fprintf(1,' [B轴承]-轴向载荷 fa7b = %3.0f N \\n',fa7b); fprintf(1,' 相对轴向载荷 xdb = %3.3f \\n',xdb); fprintf(1,' 判断参数 eb = %3.3f \\n',eb); fprintf(1,' 轴向与径向载荷之比 fbe = %3.3f \\n',fa7b/frb); fprintf(1,' 径向载荷系数 x7b = %3.2f \\n',x7b); fprintf(1,' 轴向载荷系数 y7b = %3.2f \\n',y7b); fprintf(1,' 当量动载荷 p7b = %3.0f N \\n',p7b); fprintf(1,' 轴承寿命 lh7b = %3.0f h \\n',lh7b); '5-滚动轴承寿命计算(深沟球轴承6009)' c6r=21000; % 额定动载荷(N) c60r=14800; % 额定静载荷(N) if fra>=frb '以轴承A作为计算依据' fr6=fra; else '以轴承B作为计算依据' fr6=frb; end xd6=fa/c60r; % 轴承的相对轴向载荷 ea6=0.28; % 6类轴承的判断参数(查表) if fa/fr6>ea6; % 轴承6类轴向载荷与径向载荷之比 x6a=0.56; % 确定6类轴承载荷折算系数X y6a=1.55; % 确定6类轴承载荷折算系数Y else x6a=1; y6a=0; end p6a=round(x6a*fr6+y6a*fa); % 6类轴承的当量动载荷(N) fp=1.5; % 轴承载荷系数(减速器中等冲击) lh6a=round(1e6/60/n2*(c6r/fp/p6a)^3); % 计算6类轴承的寿命(h) disp ' ' disp ' ========== 深沟球轴承寿命计算 =========='; fprintf(1,' 6009轴承额定动载荷 c6r = %3.0f N \\n',c6r); fprintf(1,' 额定静载荷 c60r = %3.0f N \\n',c60r); fprintf(1,' 相对轴向载荷 xd6 = %3.3f \\n',xd6); fprintf(1,' 判断参数 ea6 = %3.3f \\n',ea6); fprintf(1,' 轴向与径向载荷之比 fae = %3.3f \\n',fa/fr6); fprintf(1,' 径向载荷系数 x6a = %3.2f \\n',x6a); fprintf(1,' 轴向载荷系数 y6a = %3.2f \\n',y6a); fprintf(1,' 当量动载荷 p6a = %3.0f N \\n',p6a); fprintf(1,' 轴承寿命 lh6a = %3.0f h \\n',lh6a); % 主轴支承静不定结构的计算(调用zzjbd.m) disp ' ' disp ' ******** 主轴支承静不定结构的计算 ********' disp ' ' ZCLX = input(' 滚动轴承类型：角接触球轴承-\"Q\";圆锥滚子轴承-\"Z\" == '); disp ' ' if ZCLX=='Q' disp ' ======== 角接触球轴承 ========' lmd=0.3636; Qe = input(' 载荷转换参数 e = '); ctgalf=1.25/Qe; m=0.3608; n=0.6622; h=0.32; elseif ZCLX=='Z' disp ' ======== 圆锥滚子轴承 ========' lmd=0.3148; Zy = input(' 轴向载荷系数 Y = '); ctgalf=2.5*Zy; m=0.3234; n=0.6872; h=0.34; end fprintf (1,' 相对位置系数 lmd = %3.4f \\n',lmd) fprintf (1,' 接触角余切 ctgalf = %3.4f \\n',ctgalf) fprintf (1,' 曲线拟合常数 m = %3.4f \\n',m) fprintf (1,' 曲线拟合常数 n = %3.4f \\n',n) disp ' ' % 主轴结构尺寸 l=525;a=203;b=130;c=194;d=271.31;D=180; disp ' ======== 主轴的结构尺寸 ========' fprintf (1,' 支承跨度 l = %3.2f mm \\n',l) fprintf (1,' 悬臂长度 a = %3.2f mm \\n',a) fprintf (1,' 前支承宽度 b = %3.2f mm \\n',b) fprintf (1,' 齿轮位置 c = %3.2f mm \\n',c) fprintf (1,' 齿轮分度圆直径 d = %3.2f mm \\n',d) fprintf (1,' 机床切削直径 D = %3.2f mm \\n',D) disp ' ' % 齿轮传动力 Qt=4839;Qr=1820;Qa=1294;fai=128; disp ' ======== 齿轮传动力 ========' fprintf (1,' 圆周力 Qt = %3.2f N \\n',Qt) fprintf (1,' 径向力 Qr = %3.2f N \\n',Qr) fprintf (1,' 轴向力 Qa = %3.2f N \\n',Qa) fprintf (1,' 载荷方位角 fai = %3.2f 度 \\n',fai) disp ' ' % 机床切削力 Pz=6978;Py=2791;Px=4400; disp ' ======== 机床切削力 ========' fprintf (1,' 主切削力 Pz = %3.2f N \\n',Pz) fprintf (1,' 径向力 Py = %3.2f N \\n',Py) fprintf (1,' 轴向力 Px = %3.2f N \\n',Px) disp ' ' % 计算力矩参数 hd=pi/180; u=-Py*(l+a)+0.5*Px*D-(l-c)*(Qr*cos(fai*hd)-Qt*sin(fai*hd))-0.5*Qa*d*cos(fai*hd); v=-Pz*(l+a)+(l-c)*(Qr*sin(fai*hd)+Qt*cos(fai*hd))+0.5*Qa*d*sin(fai*hd); w=sqrt(u^2+v^2); disp ' ======== 力矩参数 ========' fprintf (1,' XOY面参数 u = %3.0f Nmm \\n',u) fprintf (1,' ZOY面参数 v = %3.0f Nmm \\n',v) fprintf (1,' 综合力矩参数 w = %3.0f Nmm \\n',w) Fa=Px-Qa; b1b=lmd*Fa*l*ctgalf/(w-lmd*Fa*ctgalf); if b1b<=h Fr=lmd*Fa*ctgalf/b1b; elseif b1b>h mu=w/(n*b); lo=l/(n*b)+m/n+log(Fa*ctgalf); x0=w/(l+n*b); % 前支承径向载荷Fr初值 Fr=fsolve(@zzjbd,x0); % 利用函数fsolve解非线性方程 end b1=b1b*b; % 前支承径向载荷作用位置 % 计算后支承径向载荷 Fry=u/(l+b1); % 前支承XOY面-径向载荷 Frz=v/(l+b1); % XOZ面-径向载荷 Fray=-Qr*cos(fai*hd)+Qt*sin(fai*hd)-Fry-Py; % 后支承XOY面-径向载荷 Fraz=Qr*sin(fai*hd)+Qt*cos(fai*hd)-Frz-Pz; % XOZ面-径向载荷 Fra=sqrt(Fray^2+Fraz^2); % 后支承径向载荷 disp ' ' disp ' ======== 主轴载荷和作用位置 ========' fprintf (1,' 主轴轴向载荷 Fa = %3.0f N \\n',Fa) fprintf (1,' 相对位置参数 b1b = %3.4f \\n',b1b) fprintf (1,' 前支承B-径向载荷 FrB = %3.0f N \\n',Fr) fprintf (1,' 前支承径向载荷作用位置 b1 = %3.2f mm \\n',b1) fprintf (1,' 后支承A-径向载荷 FrA = %3.0f N \\n',Fra) % 主轴支承静不定结构超越非线性方程 function f=zzjbd(x) global mu lo % 定义全局变量 f=log(x)-mu/x-lo; disp ' ****** 转轴的可靠性设计 *******' M=input(' 输入转轴危险截面上的弯矩(Nmm) M = '); T=input(' 输入转轴危险截面上的扭矩(Nmm) T = '); Kmsa=32*M/pi; fprintf (1,' 对称循环弯曲应力幅系数 Kmsa = %3.3f \\n',Kmsa) Kcsa=0.08*Kmsa; fprintf (1,' 弯曲应力幅标准离差系数 Kcsa = %3.3f \\n',Kcsa) Kmsm=16*sqrt(3)*T/pi; fprintf (1,' 稳定扭转平均应力系数 Kmsm = %3.3f \\n',Kmsm) Kcsm=0.08*Kmsm; fprintf (1,' 扭转平均应力标准离差系数 Kcsm = %3.3f \\n',Kcsm) rb=Kmsa/Kmsm; fprintf (1,' 应力幅与平均应力的比值 rb = %3.3f \\n',rb) Kmrb=sqrt(1+1/rb^2); fprintf (1,' 应力比均值系数 Kmrb = %3.3f \\n',Kmrb) Kmsf=Kmsa*sqrt(1+1/rb^2); fprintf (1,' 复合疲劳平均应力系数 Kmsf = %3.3f \\n',Kmsf) Kcsf=Kcsa*sqrt(1+1/rb^2); fprintf (1,' 复合疲劳平均应力标准离差系数 Kcsf = %3.3f \\n',Kcsf) Cb=input(' 输入转轴材料的弯曲强度极限(MPa) Cb = '); Csjdc=0.43*Cb; % 袖珍机械设计师手册(第2版),P17,表1-18,结构钢 fprintf (1,' 试件的对称循环弯曲疲劳极限(MPa) Csjbc = %3.3f \\n',Csjdc) B=input(' 输入转轴的表面质量系数 B = '); Ec=input(' 输入转轴的弯曲绝对尺寸系数 Ec = '); Et=input(' 输入转轴的扭转绝对尺寸系数 Et = '); E=(Ec+Et)/2; fprintf (1,' 转轴的弯曲绝对尺寸系数 E = %3.3f \\n',E) Kc=input(' 输入转轴的弯曲疲劳应力集中系数 Kc = '); Kt=input(' 输入转轴的扭转疲劳应力集中系数 Kt = '); Q=input(' 输入转轴的敏感系数 Q = '); Kf=1+Q*(Kc*Kt-1); fprintf (1,' 转轴的复合疲劳应力集中系数 Kf = %3.3f \\n',Kf) Cdc=Csjdc*B*E/Kf; fprintf (1,' 转轴的对称循环弯曲疲劳极限(MPa) Cdc = %3.3f \\n',Cdc) % 复合疲劳应力下强度的均值Sj,按照应力线与最佳拟合均值线的交点求出 Sj=sqrt(Cdc^2*Cb^2*(1+rb^2)/(Cb^2*rb^2+Cdc^2)); fprintf (1,' 转轴强度的均值(MPa) Sj = %3.3f \\n',Sj) Cs=0.08*Sj; fprintf (1,' 转轴强度的标准离差(MPa) Cs = %3.6f \\n',Cs) R=input(' 输入可靠度 R = '); % 根据失效概率F求联结系数z时,用累积分布反函数 z=norminv(F,mu,sigma) % 根据联结系数z求失效概率F时,用累积分布函数 F=normcdf(z,mu,sigma) % 根据联结系数z求失效频数f时,用概率密度函数 f=normpdf(z,mu,sigma) % 正态分布N(mu,sigma),x=z时,mu=0,sigma=1,为标准正态分布,mu和sigma可以省略 z=norminv(1-R); F=normcdf(z); fprintf (1,' 与可靠度R对应的失效概率 F = %3.6f \\n',F) fprintf (1,' 联结系数(可靠性系数,安全指数) z = %3.6f \\n',z) f=normpdf(z); fprintf (1,' 与联结系数z对应的失效频数 f = %3.6f \\n',f) a6=z^2*Cs^2-Sj^2; fprintf (1,' 联结方程多项式中6次方项的系数 a6 = %3.4f \\n',a6) a3=2*Sj*Kmsf; fprintf (1,' 联结方程多项式中3次方项的系数 a3 = %3.4f \\n',a3) a0=z^2*Kcsf^2-Kmsf^2; fprintf (1,' 联结方程多项式中的常数项 a0 = %3.4f \\n',a0) p=[a6 0 0 a3 0 0 a0]; % 求解多项式的根 d=roots(p); disp ' 转轴危险截面直径(联结方程多项式的根)' fprintf (1,' d1 = %3.3f mm \\n',d(1)) fprintf (1,' d2 = %3.3f mm \\n',d(2)) fprintf (1,' d3 = %3.3f mm \\n',d(3)) fprintf (1,' d4 = %3.3f mm \\n',d(4)) fprintf (1,' d5 = %3.3f mm \\n',d(5)) fprintf (1,' d6 = %3.3f mm \\n',d(6)) disp ' ****** 转轴可靠性设计的有关参数分析 *******' RR=[0.5 0.6 0.7 0.8 0.9 0.99 0.999 0.9999]; zz=norminv(1-RR); figure(1); plot(zz,RR) title('\\bf 联结系数z与可靠度R的关系曲线 \\rm R=\\int e^{-z^{2}/2}dz ') xlabel('联结系数 z') ylabel('可靠度 R') grid; aa6=zz.^2*Cs^2-Sj^2; aa3=2*Sj*Kmsf; aa0=zz.^2*Kcsf^2-Kmsf^2; p1=[aa6(1) 0 0 aa3 0 0 aa0(1)]; dr1=roots(p1); fprintf (1,' 可靠度 R=0.5000 时的联结系数 z1 = %3.6f \\n',zz(1)) fprintf (1,' 转轴直径 dr1 = %3.3f mm \\n',dr1(6)) p2=[aa6(2) 0 0 aa3 0 0 aa0(2)]; dr2=roots(p2); fprintf (1,' 可靠度 R=0.6000 时的联结系数 z2 = %3.6f \\n',zz(2)) fprintf (1,' 转轴直径 dr2 = %3.3f mm \\n',dr2(6)) p3=[aa6(3) 0 0 aa3 0 0 aa0(3)]; dr3=roots(p3); fprintf (1,' 可靠度 R=0.7000 时的联结系数 z3 = %3.6f \\n',zz(3)) fprintf (1,' 转轴直径 dr3 = %3.3f mm \\n',dr3(6)) p4=[aa6(4) 0 0 aa3 0 0 aa0(4)]; dr4=roots(p4); fprintf (1,' 可靠度 R=0.8000 时的联结系数 z4 = %3.6f \\n',zz(4)) fprintf (1,' 转轴直径 dr4 = %3.3f mm \\n',dr4(6)) p5=[aa6(5) 0 0 aa3 0 0 aa0(5)]; dr5=roots(p5); fprintf (1,' 可靠度 R=0.9000 时的联结系数 z5 = %3.6f \\n',zz(5)) fprintf (1,' 转轴直径 dr5 = %3.3f mm \\n',dr5(6)) p6=[aa6(6) 0 0 aa3 0 0 aa0(6)]; dr6=roots(p6); fprintf (1,' 可靠度 R=0.9900 时的联结系数 z6 = %3.6f \\n',zz(6)) fprintf (1,' 转轴直径 dr6 = %3.3f mm \\n',dr6(6)) p7=[aa6(7) 0 0 aa3 0 0 aa0(7)]; dr7=roots(p7); fprintf (1,' 可靠度 R=0.9990 时的联结系数 z7 = %3.6f \\n',zz(7)) fprintf (1,' 转轴直径 dr7 = %3.3f mm \\n',dr7(6)) p8=[aa6(8) 0 0 aa3 0 0 aa0(8)]; dr8=roots(p8); fprintf (1,' 可靠度 R=0.9999 时的联结系数 z8 = %3.6f \\n',zz(8)) fprintf (1,' 转轴直径 dr8 = %3.3f mm \\n',dr8(6)) dz=[dr1(6) dr2(6) dr3(6) dr4(6) dr5(6) dr6(6) dr7(6) dr8(6)]; figure(2); plot(dz,RR) title('\\bf 轴的直径d与可靠度R的关系曲线 \\it') xlabel('轴的直径 d/mm') ylabel('可靠度 R') grid; % 二维约束优化问题函数子窗口几何描述（例1分图） % 按(初值,终值,间隔数)等间隔矢量产生二维网格矩阵 xx1=linspace(-25,25,10); xx2=linspace(-20,20,10); [x1,x2]=meshgrid(xx1,xx2); % 目标函数 f=x1.^2+x2.^2-4*x1+4; subplot(2,2,1); % 选择2行2列的第1个子窗口 surfc(x1,x2,f); title('\\bf 目标函数 f=(X1)^{2}+(X2)^{2}-4(X1)+4 曲面') % 约束函数g1 g1=x1-x2+2; subplot(2,2,2); % 选择2行2列的第2个子窗口 surfc(x1,x2,g1); title('\\bf 约束函数 g1=(X1)-(X2)+2 曲面') % 约束函数g2 g2=-x1.^2+x2-1; subplot(2,2,3); % 选择2行2列的第3个子窗口 surfc(x1,x2,g2); title('\\bf 约束函数 g2=-(X1)^{2}+(X2)-1 曲面') % 约束函数g3 g3=x1; subplot(2,2,4); % 选择2行2列的第4个子窗口 surfc(x1,x2,g3); title('\\bf 约束函数 g3=(X1) 曲面') % 二维约束优化几何描述(例1组合) % 按等间隔矢量产生二维网格矩阵 xx1=linspace(-1,3,40); xx2=linspace(0,4,40); [x1,x2]=meshgrid(xx1,xx2); % 数学模型 f=x1.^2+x2.^2-4*x1+4; % 目标函数f g1=x1-x2+2; % 约束函数g1 g2=-x1.^2+x2-1; % 约束函数g2 g3=x1; % 约束函数g3 % 设计空间 figure(1); surfc(x1,x2,f); title('\\bf 目标函数和约束函数曲面'); xlabel('设计变量 \\bf X1'); ylabel('设计变量 \\bf X2'); zlabel('目标函数值和约束函数值'); hold on; % 保持图形 surfc(x1,x2,g1); surfc(x1,x2,g2); surfc(x1,x2,g3); % 设计平面 figure(2); h=contour(x1,x2,f); clabel(h); axis equal % 两坐标轴的定标因子相等 title('\\bf 设计平面'); xlabel('设计变量 \\bf X1'); ylabel('设计变量 \\bf X2'); hold on; h=contour(x1,x2,g1); clabel(h); h=contour(x1,x2,g2); clabel(h); h=contour(x1,x2,g3); clabel(h); % 二维优化几何描述(例2、例3、例4) % 按等间隔矢量产生二维网格矩阵 sx1=linspace(-6,4,30); sx2=linspace(-4,4,30); [x1,x2]=meshgrid(sx1,sx2); % 1-约束优化问题数学模型（例2） f=x1+x2.^2; % 目标函数f g1=-x1.^2-x2.^2+9; % 约束函数g1 g2=-x1-x2+1; % 约束函数g2 % 设计空间 figure(1); surfc(x1,x2,f); title('\\bf 目标函数和约束函数曲面'); xlabel('设计变量 \\bf X1'); ylabel('设计变量 \\bf X2'); zlabel('目标函数值和约束函数值'); hold on; % 保持图形 surfc(x1,x2,g1); surfc(x1,x2,g2); % 设计平面 figure(2); h=contour(x1,x2,f); clabel(h); axis equal; % 两坐标轴的定标因子相等 title('\\bf 设计平面'); xlabel('设计变量 \\bf X1'); ylabel('设计变量 \\bf X2'); hold on; h=contour(x1,x2,g1); clabel(h); h=contour(x1,x2,g2); clabel(h); % % 按等间隔矢量产生二维网格矩阵 sy1=linspace(-2,3,30); sy2=linspace(-2,4,30); [y1,y2]=meshgrid(sy1,sy2); % 2-无约束优化问题目标函数（例3） f01=y1.^4-2*y2.*y1.^2+y1.^2+y2.^2-2*y1+5; figure(3); surfc(y1,y2,f01); title('\\bf f=(X1)^{4}-2(X1)^{2}(X2)+X1^{2}+(X2)^{2}-2(X1)+5'); xlabel('设计变量 \\bf X1'); ylabel('设计变量 \\bf X2'); zlabel('目标函数值\\bf f'); figure(4); h=contour(y1,y2,f01,50); axis equal; title('\\bf f=(X1)^{4}-2(X1)^{2}(X2)+X1^{2}+(X2)^{2}-2(X1)+5 等值线'); xlabel('设计变量 \\bf X1'); ylabel('设计变量 \\bf X2'); figure(5); h=contour3(y1,y2,f01,50); title('\\bf f=(X1)^{4}-2(X1)^{2}(X2)+X1^{2}+(X2)^{2}-2(X1)+5 三维等值线'); xlabel('设计变量 \\bf X1'); ylabel('设计变量 \\bf X2'); zlabel('目标函数\\bf f'); % % 按等间隔矢量产生二维网格矩阵 sz1=linspace(-3,4,30); sz2=linspace(-2,7,30); [z1,z2]=meshgrid(sz1,sz2); % 3-无约束优化问题目标函数（例4） f02=z1.^4-2*z2.*z1.^2+z1.^2+2*z2.^2-2*z1.*z2+4.5*z1-4*z2+4; figure(6); surfc(z1,z2,f02); title('\\bf f=(X1)^{4}-2(X1)^{2}(X2)+X1^{2}+2(X2)^{2}-2(X1)(X2)-4.5(X1)-4(X2)+4'); xlabel('设计变量 \\bf X1'); ylabel('设计变量 \\bf X2'); zlabel('目标函数值\\bf f'); figure(7); h=contour(z1,z2,f02,50); axis equal; title('\\bf f=(X1)^{4}-2(X1)^{2}(X2)+X1^{2}+2(X2)^{2}-2(X1)(X2)-4.5(X1)-4(X2)+4 等值线'); xlabel('设计变量 \\bf X1'); ylabel('设计变量 \\bf X2'); figure(8); h=contour3(z1,z2,f02,50); title('\\bf f=(X1)^{4}-2(X1)^{2}(X2)+X1^{2}+2(X2)^{2}-2(X1)(X2)-4.5(X1)-4(X2)+4 三维等值线'); xlabel('设计变量 \\bf X1'); ylabel('设计变量 \\bf X2'); zlabel('目标函数\\bf f'); % 人字架体积优化设计 % 1----主程序 % 人字架优化调用目标函数文件与非线性约束文件 % 设计变量的初始值 x0=[100;800]; % 设计变量的下界与上界 lb=[20;200]; ub=[140;1200]; % 线性不等式约束(g3、g4、g5、g6)中设计变量的系数矩阵 a=zeros(4,2); a(1,1)=-1;a(2,1)= 1; a(3,2)=-1;a(4,2)= 1; % 线性不等式约束中常数项列阵 b=[-20;140;-200;1200]; % 调用约束优化函数 % 等式约束参数Aeq,beq定义为空矩阵符号“[ ]” [x,fn]=fmincon(@rzjyh_f,x0,a,b,[],[],lb,ub,@rzjyh_g); disp ' ******** 人字架体积优化设计最优解 ********' fprintf (1,' 钢管平均直径 D = %3.4f mm \\n',x(1)) fprintf (1,' 人字架高度 H = %3.4f mm \\n',x(2)) fprintf (1,' 人字架体积 V = %3.4f mm^3 \\n',fn) % 调用约束优化非线性约束函数(jsqyh_g)计算最优点x*的性能约束函数值 g=rzjyh_g(x); disp ' ======== 最优点的性能约束函数值 ========' fprintf (1,' 人字架钢管压缩强度 g1 = %3.4f \\n',g(1)) fprintf (1,' 人字架钢管稳定性 g2 = %3.4f \\n',g(2)) % 2----目标函数(rzjyh_f) function f=jsqyh_f(x); % 人字架跨距B;钢管厚度T; B=1520;T=2.5; f=2*pi*x(1)*T*sqrt((B/2)^2+x(2)^2); % 3----非线性不等式约束函数(rzjyh_g) function [g,ceq] = rzjyh_g(x); % 人字架跨距B;钢管厚度T;载荷P;弹性模量E;许用压应力Cy; B=1520;T=2.5;P=294300;E=2.119e5;Cy=690; % 钢管压缩强度条件 Q=0.5*P*sqrt((B/2)^2+x(2)^2)/x(2); % 钢管轴向压力 sgm=Q/(pi*T*x(1)); % 钢管压应力 g(1)=sgm-Cy; % 钢管稳定性条件 Cc=0.125*pi^2*E*(x(1)^2+T^2)/((B/2)^2+x(2)^2); % 稳定临界应力 g(2)=sgm-Cc; ceq=[]; % 二维优化几何描述(人字架体积优化设计) % 按等间隔矢量产生二维网格矩阵 xx1=linspace(20,140,20); % D取值范围 20- 140 xx2=linspace(200,1200,200); % H取值范围200-1200 [x1,x2]=meshgrid(xx1,xx2); % 数学模型 % 人字架跨距B;钢管厚度T;载荷P;弹性模量E;许用压应力Cy; B=1520;T=2.5;P=294300;E=2.119e5;Cy=690; f=2*pi.*x1*T.*sqrt((B/2)^2+x2.^2); % 目标函数f % 目标函数值几何描述 fh=contour(x1,x2,f); % 目标函数等高线 clabel(fh); % 标注目标函数值 title('\\bf 人字架体积优化设计平面'); xlabel('钢管平均直径 \\bf D (mm) '); ylabel('人字架高度 \\bf H (mm) '); % 钢管压缩强度条件 Q=0.5*P*sqrt((B/2)^2+x2.^2)./x2; % 钢管轴向力 sgm=Q./(pi.*x1*T); % 钢管压应力 g1=sgm-Cy; % 约束函数g1 % 钢管稳定性条件 Cc=0.125*pi^2*E.*(x1.^2+T^2)./((B/2)^2+x2.^2); % 稳定临界应力 g2=sgm-Cc; % 约束函数g2 % 约束函数几何描述 hold on; g1h=contour(x1,x2,g1); g2h=contour(x1,x2,g2); % 例1-水槽定截面时周长最小的二维无约束优化 % 1----无约束优化函数命令程序 % 初始点 x0=[25;45]; % 调用梯度法搜索 [x,Fmin]=fminunc('sc_wysyh',x0); disp ' ******** 输出最优解 ********' fprintf (1,' 截面高度h x(1)* = %3.4f mm \\n',x(1)) fprintf (1,' 斜边夹角theta x(2)* = %3.4f 度 \\n',x(2)) fprintf (1,' 截面周长s f* = %3.4f mm \\n',Fmin) % 2----二维无约束优化目标函数文件(sc_wysyh.m) function f=sc_wysyh(x) a=516;hd=pi/180; f=a/x(1)-x(1)/tan(x(2)*hd)+2*x(1)/sin(x(2)*hd); % 3----绘制水槽截面周长等高线和曲面图的程序 % 按(初值,终值,等分数)产生等间隔向量xx1,xx2 xx1=linspace(100,300,25); xx2=linspace(30,120,25); % 产生两个[5x10]的网格矩阵x1,x2 [x1,x2]=meshgrid(xx1,xx2); % 定义目标函数 a=516;hd=pi/180; f=a./x1-x1./tan(x2*hd)+2*x1./sin(x2*hd); % 将整个图形窗口分隔成2个子窗口,取左边窗口 subplot(1,2,1); % 绘制等值线并标注函数值 h=contour(x1,x2,f); clabel(h); % 定义左边窗口坐标轴刻度范围 axis([100 300 30 120]) % 标注左边窗口和坐标轴 xlabel('高度 h (mm)') ylabel('倾斜角 theta (度)') title('目标函数(截面周长)等值线') % 将整个图形窗口分隔成2个子窗口,取右边窗口 subplot(1,2,2); % 绘制曲面图 surfc(x1,x2,f); % 定义右边窗口坐标轴刻度范围 axis([100 300 30 120 600 1200]) % 标注右边窗口 title('目标函数(截面周长)曲面图') % 例2-两级斜齿轮传动中心距优化设计 % 1----减速器中心距优化设计主程序 % 设计变量的初始值 x0=[2;4;18;20;6.4;10]; % 设计变量的下界与上界 lb=[2;3.5;14;16;5.8;8]; ub=[5;6;22;22;7;15]; % 线性不等式约束(g6(x)-g17(x))中设计变量的系数矩阵 a=zeros(12,6); a(1,1)=-1;a(2,1)= 1; a(3,2)=-1;a(4,2)= 1; a(5,3)=-1;a(6,3)= 1; a(7,4)=-1;a(8,4)= 1; a(9,5)=-1;a(10,5)= 1; a(11,6)=-1;a(12,6)= 1; % 线性不等式约束(g6(x)-g17(x))中的常数项列阵 b=[-2;5;-3.5;6;-14;22;-16;22;-5.8;7;-8;15]; % 使用约束优化命令fmincon(调用目标函数jsqyh_f和非线性约束函数jsqyh_g) % 没有等式约束,参数Aeq和beq定义为空矩阵符号“[ ]” [x,fn]=fmincon(@jsqyh_f,x0,a,b,[],[],lb,ub,@jsqyh_g); disp ' ******** 两级斜齿轮传动中心距优化设计最优解 ********' fprintf (1,' 高速级齿轮副模数 Mn1 = %3.4f mm \\n',x(1)) fprintf (1,' 低速级齿轮副模数 Mn2 = %3.4f mm \\n',x(2)) fprintf (1,' 高速级小齿轮齿数 z1 = %3.4f \\n',x(3)) fprintf (1,' 低速级小齿轮齿数 z3 = %3.4f \\n',x(4)) fprintf (1,' 高速级齿轮副传动比 i1 = %3.4f \\n',x(5)) fprintf (1,' 齿轮副螺旋角 beta = %3.4f 度 \\n',x(6)) fprintf (1,' 减速器总中心距 a12 = %3.4f mm \\n',fn) % 调用约束优化非线性约束函数(jsqyh_g)计算最优点x*的性能约束函数值 g=jsqyh_g(x); disp ' ======== 最优点的性能约束函数值 ========' fprintf (1,' 高速级齿轮副接触疲劳强度约束函数值 g1 = %3.4f \\n',g(1)) fprintf (1,' 低速级齿轮副接触疲劳强度约束函数值 g2 = %3.4f \\n',g(2)) fprintf (1,' 高速级大齿轮齿根弯曲强度约束函数值 g3 = %3.4f \\n',g(3)) fprintf (1,' 低速级大齿轮齿根弯曲强度约束函数值 g4 = %3.4f \\n',g(4)) fprintf (1,' 大齿轮齿顶与轴不干涉几何约束函数值 g5 = %3.4f \\n',g(5)) % 2----两级斜齿轮减速器总中心距的目标函数(jsqyh_f) function f=jsqyh_f(x); hd=pi/180; a1=x(1)*x(3)*(1+x(5)); a2=x(2)*x(4)*(1+31.5/x(5)); cb=2*cos(x(6)*hd); f=(a1+a2)/cb; % 3----两级斜齿轮减速器优化设计的非线性不等式约束函数(jsqyh_g) function [g,ceq]=jsqyh_g(x); hd=pi/180; g(1)=cos(x(6)*hd)^3-3.079e-6*x(1)^3*x(3)^3*x(5); g(2)=x(5)^2*cos(x(6)*hd)^3-1.701e-4*x(2)^3*x(4)^3; g(3)=cos(x(6)*hd)^2-9.939e-5*(1+x(5))*x(1)^3*x(3)^2; g(4)=x(5)^2.*cos(x(6)*hd)^2-1.076e-4*(31.5+x(5))*x(2)^3*x(4)^2; g(5)=x(5)*(2*(x(1)+50)*cos(x(6)*hd)+x(1)*x(2)*x(3))-x(2)*x(4)*(31.5+x(5)); ceq=[]; % 圆柱螺旋弹簧的多目标优化设计 % 设计变量:簧丝直径x(1)=d,中径x(2)=D2,圈数x(3)=n % 1----主程序(极小极大值函数fminimax) % 设计变量的初始值 x0=[6.2;39;5]; % 设计变量的下限与上限 lb=[2.5;27.5;3]; ub=[9;51;6]; % 线性不等式约束(g6(x)-g11(x))中设计变量的系数矩阵(非零元素) a=zeros(6,3); a(1,1)=-1; a(2,1)= 1; a(3,1)=-1;a(3,2)=-1; a(4,1)= 1;a(4,2)= 1; a(5,3)=-1; a(6,3)= 1; % 线性不等式约束(g6(x)-g11(x))中的常数项列阵 b=[-2.5;9;-30;60;-3;6]'; % 没有等式约束,参数Aeq和beq定义为空矩阵符号“[ ]” L=[0.9434 42.5514 1.709e-3]; % 多目标优化函数fminimax(调目标函数TH_dmbyh_fTS和非线性约束函数TH_dmbyh_gTS) for i=1:4 H=[11*L(1)-i*L(1) 11.5*L(2)-i*L(2) 13.5*L(3)-i*L(3)]; [x,fn]=fminimax('TH_dmbyh_fTS',x0,a,b,[],[],lb,ub,'TH_dmbyh_gTS',[],L,H); for j=1:3 ff(j)=fn(j)*(H(j)-L(j))+L(j); % 将各分目标函数还原成实际值 end f1(i)=ff(1);f2(i)=ff(2);f3(i)=1./ff(3); % 将各分目标最优值赋给各数组 end disp ' ******** 圆柱螺旋弹簧多目标优化设计最优解 ********' fprintf (1,' 簧丝直径 x(1)* d = %3.4f mm \\n',x(1)) fprintf (1,' 弹簧中径 x(2)* D2 = %3.4f mm \\n',x(2)) fprintf (1,' 弹簧圈数 x(3)* n = %3.4f 圈 \\n',x(3)) fprintf (1,' 弹簧结构重量 f1* W = %3.4f N \\n',ff(1)) fprintf (1,' 弹簧自由高度 f2* H0 = %3.4f mm \\n',ff(2)) fprintf (1,' 弹簧自振频率 f3* fr = %3.4f Hz \\n',1/ff(3)) % 调用约束优化非线性约束函数(TH_dmbyh_gTS)计算最优点x*的性能约束函数值 g=TH_dmbyh_gTS(x,L,H); disp ' ######## 最优点的性能约束函数值 ########' fprintf (1,' 弹簧强度约束函数值 g1* = %3.4f \\n',g(1)) fprintf (1,' 弹簧指数约束函数值 g2* = %3.4f \\n',g(2)) fprintf (1,' 弹簧稳定性约束函数值 g3* = %3.4f \\n',g(3)) fprintf (1,' 弹簧无共振约束函数值 g4* = %3.4f \\n',g(4)) fprintf (1,' 弹簧刚度约束函数值 g5* = %3.4f \\n',g(5)) disp ' ======== 最优点的边界约束函数值 ========' fprintf (1,' 簧丝直径下限约束函数值 g6* = %3.4f \\n',-b(1)-x(1)) fprintf (1,' 上限约束函数值 g7* = %3.4f \\n',x(1)-b(2)) fprintf (1,' 弹簧外径下限约束函数值 g8* = %3.4f \\n',-b(3)-x(2)) fprintf (1,' 上限约束函数值 g9* = %3.4f \\n',x(2)-b(4)) fprintf (1,' 工作圈数下限约束函数值 g10* = %3.4f \\n',-b(5)-x(3)) fprintf (1,' 上限约束函数值 g11* = %3.4f \\n',x(3)-b(6)) % 描述多目标优化决策的极小极大值和理想平面 [fsort]=sort([f1;f2;f3]); % 对数组f1,f2,f3元素按升序进行排序 disp ' &&&&&&&& 多目标的极小极大值矩阵 &&&&&&&&' disp (fsort) subplot(1,2,1); plot3([fsort(1,1)],[fsort(2,1)],[fsort(3,1)],'go'); hold on; plot3([fsort(1,2)],[fsort(2,2)],[fsort(3,2)],'ro'); plot3([fsort(1,3)],[fsort(2,3)],[fsort(3,3)],'co'); plot3([fsort(1,4)],[fsort(2,4)],[fsort(3,4)],'b*'); grid; title('多目标函数的极小极大值'); xlabel('结构重量 f1 (N)'); ylabel('自由高度 f2 (mm)'); zlabel('自振频率 f3 (Hz)'); rotate3d; subplot(1,2,2); plot3([ff(1),0,0],[0,ff(2),0],[0,0,1/ff(3)],'r*-'); grid; title('多目标函数的理想平面') xlabel('结构重量 f1 (N)'); ylabel('自由高度 f2 (mm)'); zlabel('自振频率 f3 (Hz)'); rotate3d; % 鼠标控制三维视角变化的函数 % 2----圆柱螺旋弹簧的多目标优化设计的目标函数(TH_dmbyh_fTS) function [f]=TH_dmbyh_fTS(x,L,H); % 钢的密度7.5e-6kg/mm^3;重力加速度9.80665m/s^2 p=7.5e-6;gl=9.80665; % 对各分目标函数进行相同数量级的变换 ff=(f-L)/(H-L) f1=p*gl*pi^2*x(1)^2*x(2)*(x(3)+1.8)/4; % 弹簧重量 f(1)=(f1-L(1))/(H(1)-L(1)); f2=x(1)*(x(3)+1.3)+18.25; % 弹簧自由高度 f(2)=(f2-L(2))/(H(2)-L(2)); f3=2.809e-6*x(2)^2*x(3)/x(1); % 弹簧自振频率的倒数 f(3)=(f3-L(3))/(H(3)-L(3)); % 3----圆柱螺旋弹簧的多目标优化设计的非线性不等式约束函数(TH_dmbyh_gTS) function [g,ceq]=TH_dmbyh_gTS(x,L,H); % 载荷F=680N;许用应力t=405MPa;工作频率fr=25;工作行程h=16.59 % 弹簧曲度系数K=1.6/C^0.14(当工作循环次数N>1e3时);稳定临界高径比HD=5.3 F=680;t=405;K=1.6;HD=5.3;fr=25;h=16.59;G=8e4; g(1)=8*K*F/pi*x(2)^0.86/x(1)^2.86-t; g(2)=6-x(2)/x(1); g(3)=(x(3)+1.3)*x(1)+18.25-HD*x(2); g(4)=10*fr-3.56e5*x(1)/x(2)^2/x(3); g(5)=F-h*G/8*x(1)^4/x(2)^3/x(3); ceq=[]; % 无心磨削工艺参数优化(调用wxmx_f和wxmx_g) % 1----主程序 % 设计变量 % x(1)----工件线速度(m/min); % x(2)----工件轴向进给量(mm/r); % x(3)----磨削深度(mm). % 设计变量的初始值 x0=[60;5.3;0.016]; % 设计变量的下界与上界 vlb=[12.241;3.047;0.002]; vub=[88.579;7.621;0.030]; % 如果没有设计变量边界，则边界参数lb,ub定义为空矩阵符号“[ ]” % 六个线性不等式约束(g6、g7、g8、g9、g10、g11)中设计变量的系数矩阵 a=zeros(6,2); a(1,1)=-1;a(2,1)= 1; a(3,2)=-1;a(4,2)= 1; a(5,3)=-1;a(6,3)= 1; % 六个线性不等式约束中常数项列阵 b=[-12.241;88.579;-3.047;7.621;-0.002;0.030]; % 调用约束优化函数 % 等式约束参数Aeq,beq定义为空矩阵符号“[ ]” [x,fn]=fmincon(@wxmx_f,x0,a,b,[],[],vlb,vub,@wxmx_g); fmax=1e3*x(1)*x(2)*x(3); disp ' ' disp ' ******** 计 算 结 果 ********' disp ' ' disp ' ******** 无心磨削工艺参数的最优解 ********' fprintf (1,' 工件线速度 vw = %3.4f m/min \\n',x(1)) fprintf (1,' 工件轴向进给量 fa = %3.4f mm/r \\n',x(2)) fprintf (1,' 磨削深度 t = %3.4f mm \\n',x(3)) fprintf (1,' 金属切除率 fmax = %3.4f mm^3/min \\n',fmax) disp ' ' % 调用三维约束优化非线性约束函数(wxmx_g)计算最优点x*的性能约束函数值 g=wxmx_g(x); disp ' ======== 最优点的性能约束函数值 ========' fprintf (1,' 表面粗糙度条件 g1 = %3.4f \\n',g(1)) fprintf (1,' 防止磨削烧伤条件 g2 = %3.4f \\n',g(2)) fprintf (1,' 磨轮耐用度条件 g3 = %3.4f \\n',g(3)) fprintf (1,' 磨床主电机功率条件 g4 = %3.4f \\n',g(4)) fprintf (1,' 工件轴向速度条件 g5 = %3.4f \\n',g(5)) disp ' ' disp ' ======== 最优点的边界约束函数值 ========' fprintf (1,' 工件线速度最小值 g6 = %3.4f \\n',12.241-x(1)) fprintf (1,' 最大值 g7 = %3.4f \\n',x(1)-88.579) fprintf (1,' 工件轴向进给量最小值 g8 = %3.4f \\n',3.047-x(2)) fprintf (1,' 最大值 g9 = %3.4f \\n',x(2)-7.621) fprintf (1,' 磨削深度最小值 g10 = %3.4f \\n',0.002-x(3)) fprintf (1,' 最大值 g11 = %3.4f \\n',x(3)-0.03) % 2----目标函数(wxmx_f) function f=wxmx_f(x); % 金属切除率(mm^3/min)最大化 f=1e3/(x(1)*x(2)*x(3)); % 3----非线性不等式约束函数(wxmx_g) function [g,ceq] = wxmx_g(x); hd=pi/180; % 约束条件1-工件表面最大高度不超过表面粗糙度要求条件 Bw=28; % 工件宽度(mm); dw=55.56; % 工件直径(mm); dr=500; % 导轮直径(mm); Rz=2; % 表面粗糙度(微米); Kh=80; % 表面粗糙度系数; u=0.22; % 磨轮切削刃平均间隔(mm); Lmd=53; % 磨轮磨刃半顶角(度); n=1300; % 磨轮转速(r/min); v=1e-3*pi*dr*n; % 磨轮线速度(m/min); g(1)=1.36*Kh*u^1.2*(x(1)*x(2)/(tan(Lmd*hd)*v*Bw))^0.4*(1/dr+1/dw)^0.2-Rz; % 约束条件2-防止磨削烧伤条件 Cb=1920; % 磨削烧伤临界系数(m.mm/min); g(2)=(v-x(1))*x(3)^0.5*(dr*dw/(dr+dw))^0.5-Cb; % 约束条件3-磨轮耐用度条件 Ct=2550; % 工件材料系数; Tb=30; % 磨轮耐用度适用值(min); g(3)=Tb-Ct*dw^0.6/(x(1)*x(2))^1.82/x(3)^1.1; % 约束条件4-磨床主电机功率条件 Pc=13; % 磨床主电机功率(kW); nu=0.95; % 磨床主电机到主轴之间的传动效率; Tb=30; % 磨轮耐用度适用值(min); g(4)=0.0358*(1e3*x(1)*x(2)*x(3)/pi)^0.7-nu*Pc; % 约束条件5-工件轴向速度条件 vam=2000; % 工件轴向速度最大值(mm/min); g(5)=1e3*(x(1)*x(2))/(pi*dw)-vam; % 没有非线性等式约束条件 ceq=[]; % 4----计算结果分析 disp ' ******** 磨削工艺参数的凑整解 ********' dw=55.56; % 工件直径(mm); dr=500; % 导轮直径(mm); hd=pi/180; alf=atan(x(2)/(pi*dw)); % 导轮偏角(rad) nr=1e3*x(1)/(pi*dr*cos(alf)); % 导轮转速(r/min) nrj=[13 17 23 30 40 53 71 94]; % 磨床导轮的八级转速 for i=1:8 if nr<=nrj(i) nrd=nrj(i-1);break % 确定导轮转速 end end vwz=1e-3*pi*dr*nrd*cos(alf); % 工件线速度凑整解(m/min) fnz=(pi*dr*nrd*cos(alf))*(pi*dw*tan(alf))*x(3); % 金属切除率凑整解 alfy=2.50;nry=30;ty=0.01; % 原精磨工艺参数 fay=pi*dw*tan(alfy*hd); % 原工件轴向进给量(mm/r) vwy=1e-3*pi*dr*nry*cos(alfy*hd); % 原工件线速度(m/min) fny=1e3*vwy*fay*ty; % 原金属切除率(mm^3/min) fprintf (1,' 磨床导轮偏角 alf = %3.4f 度 \\n',alf/hd) fprintf (1,' 导轮转速 nr = %3.0f r/min \\n',nrd) fprintf (1,' 工件线速度 vwz = %3.4f m/min \\n',vwz) fprintf (1,' 金属切除率凑整解 fnz = %3.4f mm^3/min \\n',fnz) disp ' ' disp ' ******** 原磨削工艺参数 ********' fprintf (1,' 原磨床导轮偏角 alfy = %3.4f 度 \\n',alfy) fprintf (1,' 原导轮转速 nry = %3.0f r/min \\n',nry) fprintf (1,' 原工件线速度 vwy = %3.4f m/min \\n',vwy) fprintf (1,' 原工件轴向进给量 fay = %3.4f mm/r \\n',fay) fprintf (1,' 原工件磨削深度 ty = %3.4f mm \\n',ty) fprintf (1,' 原金属切除率 fny = %3.4f mm^3/min \\n',fny) disp ' ' fprintf (1,' 新旧工艺金属切除率之比 fnb = %3.4f \\n',fnz/fny) fprintf (1,' 金属切除率提高率 fnv = %3.4f \\n',(fnz-fny)/fny) % 绘制6种类型曲线 a = 5; lx=input('选择曲线类型：1-双曲线;2-幂函数;3-负指数;4-S型;5-指数;6-对数;7-直线 = '); b=input('输入正系数或负系数：b = '); fprintf (1,' 常数 a = %3.6f \\n',a) fprintf (1,' 系数 b = %3.6f \\n',b) xx = 1 : 0.1 : 10; switch lx case 1 yy = xx ./ (a .* xx + b ); case 2 yy = a * xx .^ b; case 3 yy = a * exp(b ./ xx); case 4 yy = 1 ./ (a + b .* exp(-xx)); case 5 yy = a * exp(b .* xx); case 6 yy = a + b .* log(xx) / log(10); case 7 yy = a + b .* xx; end subplot(3,2,lx) plot(xx, yy,'b-') grid % 图形标题-曲线类型和公式(字符黑体,公式斜体) switch lx case 1 title('\\bf 双曲线 \\it y=x/(ax+b)') case 2 title('\\bf 幂函数曲线 \\it y=ax^b') case 3 title('\\bf 负指数曲线 \\it y=ae^{b/x}') case 4 title('\\bf S型曲线 \\it y=1/(a+be^{-x})') case 5 title('\\bf 指数曲线 \\it y=ae^{bx}') case 6 title('\\bf 对数曲线 \\it y=a+bln(x)') case 7 title('\\bf 直线 \\it y=a+bx') end % 曲线拟合 disp '***** 圆柱齿轮复合齿形系数Yfs与当量齿数Zv的关系 ******' lx=input('选择曲线类型：1-双曲线;2-幂函数;3-负指数;4-S型;5-指数;6-对数;7-直线 = '); x0=input(' 输入插值点：x0 = '); n = 18; x=zeros(n); y=zeros(n); ys=zeros(n); ya=zeros(n); s=zeros(n); x=[12 13 14 15 16 17 18 19 20 25 30 35 40 45 50 60 70 80]; y=[5.05 4.91 4.79 4.70 4.61 4.55 4.48 4.43 4.38 4.22 4.13 4.08 4.05 4.02 4.01 3.88 3.88 3.88]; x1 = x(1); xn = x(n); for i = 1 : n ya(i) = y(i); end disp ' Zv Yfs' for i = 1 : n zy = [x(i) y(i)]; disp(zy) end % 绘制实验数据列表离散点图(红色,圆点,栅格,保持图形) plot(x,y,'r.') grid hold on; % 根据曲线类型进行变量代换 for i = 1 : n switch lx case 1 x(i) = 1 / x(i); y(i) = 1 / y(i); case 2 x(i) = log(x(i)); y(i) = log(y(i)); case 3 x(i) = 1 / x(i); y(i) = log(y(i)); case 4 x(i) = exp(-x(i)); y(i) = log(y(i)); case 5 y(i) = log(y(i)); case 6 x(i) = log(x(i)) / log(10); case 7 end end disp ' ***** 计算代换变量的累加与平方和 *****' c = 0; d = 0; xc = 0; yc = 0; xy = 0; for i = 1 : n c = c + x(i); d = d + y(i); xc = xc + x(i) ^2; yc = yc + y(i) ^2; xy = xy + x(i) * y(i); end cxy = [c d xc yc xy]; disp ' C-X C-Y C-X^2 C-Y^2 disp (cxy) dx = xc - c ^2 / n; dy = yc - d ^2 / n; dc = xy - c * d / n; lxy = [dx dy dc]; disp ' Lxx Lyy Lxy' disp (lxy) disp ' ***** 计算拟合系数和常数 *****' b = dc / dx; fprintf (1,' 拟合系数 b = %3.6f \\n',b) a = (d - b * c) / n; fprintf (1,' 代换常数 a = %3.6f \\n',a) switch lx case 1 fprintf (1,' 拟合常数 a = %3.6f \\n',a) C-(XY)^2' case 2 g = exp(a); fprintf (1,' 拟合常数 a = %3.6f \\n',g) case 3 g = exp(a); fprintf (1,' 拟合常数 a = %3.6f \\n',g) case 4 fprintf (1,' 拟合常数 a = %3.6f \\n',a) case 5 g = exp(a); fprintf (1,' 拟合常数 a = %3.6f \\n',g) case 6 fprintf (1,' 拟合常数 a = %3.6f \\n',a) case 7 fprintf (1,' 拟合常数 a = %3.6f \\n',a) end disp ' ***** 拟合效果和精度检验 *****' r = dc / sqrt(dx * dy); fprintf (1,' 相关系数 r = %3.6f \\n',r) d = 0; q = 0; yc = 0; for i = 1 : n ys(i) = a + b * x(i); switch lx case 1 ys(i) = 1 / ys(i); case 2 ys(i) = exp(ys(i)); case 3 ys(i) = exp(ys(i)); case 4 ys(i) = 1 / ys(i); case 5 ys(i) = exp(ys(i)); otherwise end s(i) = (ya(i) - ys(i)) ^ 2; d = d + ya(i); q = q + s(i); yc = yc + ya(i) ^2; end dy = yc - d * d / n; sgm = sqrt(q / (n - 2)); rr = 1 - q / dy; fprintf (1,' 剩余平方和 Q = %3.6f \\n',q) fprintf (1,' 标准误差 s = %3.6f \\n',sgm) fprintf (1,' 相关指数 R2 = %3.6f \\n',rr) switch lx case 7 f = (dy - q) * (n - 2) / q; fprintf (1,' 统计量 F = %3.6f \\n',f) otherwise end % 计算插值点函数值,根据选择的曲线类型绘制拟合曲线(蓝色,实线) xx = x1 : 1 : xn; switch lx case 1 y0 = x0 / (a * x0 + b ); yy = xx ./ (a .* xx + b ); case 2 y0 = g * x0 ^ b; yy = g * xx .^ b; case 3 y0 = g * exp(b / x0); yy = g * exp(b ./ xx); case 4 y0 = 1 / (a + b * exp(-x0)); yy = 1 ./ (a + b .* exp(-xx)); case 5 y0 = g * exp(b * x0); yy = g * exp(b .* xx); case 6 y0 = a + b * log(x0) / log(10); yy = a + b .* log(xx) / log(10); case 7 y0 = a + b * x0; yy = a + b .* xx; end fprintf (1,' 插值点函数值 y0 = %3.4f \\n',y0) plot(xx, yy,'b-') % 图形标题-拟合曲线类型和公式(字符黑体,公式斜体) switch lx case 1 title('\\bf 实验数据离散点图 / 双曲线 \\it y=x/(ax+b)') case 2 title('\\bf 实验数据离散点图 / 幂函数曲线 \\it y=ax^b') case 3 title('\\bf 实验数据离散点图 / 负指数曲线 \\it y=ae^{b/x}') case 4 title('\\bf 实验数据离散点图 / S型曲线 \\it y=1/(a+be^{-x})') case 5 title('\\bf 实验数据离散点图 / 指数曲线 \\it y=ae^{bx}') case 6 title('\\bf 实验数据离散点图 / 对数曲线 \\it y=a+bln(x)') case 7 title('\\bf 实验数据离散点图 / 直线 \\it y=a+bx') end % 最小二乘法多项式拟合 % ******** 螺旋副效率y与螺纹导程角x的关系 ******** % (1)-----选择拟合多项式拟合的阶数 x=[10 20 30 40 50 60 70]; y=[0.63 0.76 0.80 0.82 0.82 0.80 0.70]; n = 7; x1 = x(1); xn = x(n); % n个数据可以拟合(n-1)阶多项式,高阶多项式多次求导,数值特性变差 % polytool()是交互式函数,在图形上方[Degree]框中输入阶数,右击左下角的[Export]输出图形 polytool(x,y,1) % 观察多项式拟合的图形,选择置信区间最小的多项式阶数 % (2)-----计算多项式的各项系数和拟合值 m=input(' 输入多项式拟合的阶数 m = '); [p,S]=polyfit(x,y,m); disp ' 输出多项式的各项系数' fprintf (1,' a = %3.16f \\n',p) disp ' 输出多项式的有关信息 S' disp (S) [yh,delta]=polyconf(p,x,S); disp ' 观测数据 拟合数据' disp ' x y yh' for i = 1 : n xy = [x(i) y(i) yh(i)]; disp (xy) end % (3)-----绘制观测数据离散点图和多项式曲线 plot(x,y,'r.') title('\\bf 实验数据离散点图 / 多项式曲线 \\it y = a0+a1x+a2x^2+a3x^3+...') grid hold on; xi=[x1:0.1:xn]; yi=polyval(p,xi); plot(xi,yi,'k-') % (4)-----拟合效果和精度检验 Q=sum((y-yh).^2); SGM = sqrt(Q / (n - 2)); RR = sum((yh-mean(y)).^2)/sum((y-mean(y)).^2); fprintf (1,' 剩余平方和 Q = %3.6f \\n',Q) fprintf (1,' 标准误差 Sigma = %3.6f \\n',SGM) fprintf (1,' 相关指数 RR = %3.6f \\n',RR) x0=input(' 输入插值点 x0 = '); y0=polyval(p,x0); fprintf (1,' 输出插值点拟合函数值 y0 = %3.4f \\n',y0) 因篇幅问题不能全部显示，请点此查看更多更全内容
Copyright © 2019- fupindai.com 版权所有 赣ICP备2024042792号-2
违法及侵权请联系：TEL:199 1889 7713 E-MAIL:2724546146@qq.com
本站由北京市万商天勤律师事务所王兴未律师提供法律服务