%注意求刚度的时候用k=F/delta，而不是用k=F/(B*delta)。
%数据来源：多篇参考文献
%采用方法：Webster方法（材料力学方法）
%两个都是外齿轮啮合
%原理参考《机械原理》P184和P185
clear all;
clc;
format long;
F1=170000; %齿轮作用力（摆设，后面求刚度的时候会约掉）
F2=170000;
pi=3.14;
%%%%%%%%%%%%%%齿轮副基本参数及尺寸%%%%%%%%%%%%%%%%%
z1=25;z2=39; %《直齿轮轮齿变形计算的数值积分法》颜海燕
B1=20/1000;B2=20/1000;
mm=3.0/1000;
% z1=37;z2=96; %《齿轮时变啮合刚度改进算法及刚度激励研究》李亚鹏
% mm=11/1000;
% z1=45;z2=90;
% mm=3/1000;
% B1=mm*8;B2=mm*8;
ha=1;c=0.25;
alpha=20*pi/180;
invalpha=0.01490438;
d1=mm*z1;r1=d1/2;%分度圆半径
d2=mm*z2;r2=d2/2;
% B1=d1*0.4;
% B2=d2*0.4;
da1=(z1+2*ha)*mm;ra1=da1/2;%齿顶圆半径
da2=(z2+2*ha)*mm;ra2=da2/2;
df1=(z1-2*ha-2*c)*mm;rf1=df1/2;%齿根圆半径
df2=(z2-2*ha-2*c)*mm;rf2=df2/2;
db1=d1*cos(alpha);rb1=db1/2;%基圆半径
db2=d2*cos(alpha);rb2=db2/2;
s1=pi*mm/2;s2=pi*mm/2; %齿厚
alphaa1=acos(rb1/ra1); %齿顶圆压力角
alphaa2=acos(rb2/ra2);
a=mm*(z1+z2)/2; %标准中心距
ap=a;%实际中心距
alphap=acos(a*cos(alpha)/ap);
invalphap=tan(alphap)-alphap;
dp1=d1*cos(alpha)/cos(alphap);%节圆直径1
rp1=dp1/2;
sp1=s1*rp1/r1-2*rp1*(invalphap-invalpha);%节圆齿厚1
Hp1=sp1;%节圆齿厚1
dp2=d2*cos(alpha)/cos(alphap);%节圆直径2
rp2=dp2/2;
sp2=s2*rp2/r2-2*rp2*(invalphap-invalpha);%节圆齿厚2
Hp2=sp2;%节圆齿厚2
p=pi*mm;
Pb=p*cos(alpha);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%齿轮的材料参数%%%%%%%%%%%%%%
E1=2.09e11; %单位为Pa
E2=2.09e11;
miu1=0.26;
miu2=0.26;
if B1/Hp1>5
Ee1=E1/(1-miu1^2);
else
Ee1=E1;%轮齿材料的等效弹性模量
end
if B2/Hp2>5
Ee2=E2/(1-miu2^2);
else
Ee2=E2;%轮齿材料的等效弹性模量
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%传动基本参数区%%%%%%%%%%%%%%%%%
PB1=rb1*(tan(alphaa1)-tan(alphap));
PB2=rb2*(tan(alphaa2)-tan(alphap));
B1B2=PB1+PB2;
B2C=B1B2-Pb;
N1B1=sqrt(ra1^2-rb1^2);
N1B2=N1B1-B1B2;
N2B2=sqrt(ra2^2-rb2^2);
N2B1=N2B2-B1B2;
CD=Pb-B2C;
N1C=N1B2+B2C;
N2C=N2B1+Pb;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=9;n=40;
N1Bi=zeros(1,m);
O1Bi=zeros(1,m);
rp1i=zeros(1,m);
alphap1i=zeros(1,m);
invalphap1i=zeros(1,m);
faip1i=zeros(1,m);
xp1i=zeros(1,m);
yp1i=zeros(1,m);
beta1i=zeros(1,m);
Lf1i=zeros(1,m);
rou1i=zeros(1,m);
rK1i=zeros(m,n+1);
x1i=zeros(m,n+1);
y1i=zeros(m,n+1);
Tk1i=zeros(m,n+1);
alphaK1i=zeros(m,n+1);
invalphaK1i=zeros(m,n+1);
A1i=zeros(m,n+1);
I1i=zeros(m,n+1);
equivalentA1i=zeros(m,n);
equivalentI1i=zeros(m,n);
N2Bi=zeros(1,m);
O2Bi=zeros(1,m);
rp2i=zeros(1,m);
alphap2i=zeros(1,m);
invalphap2i=zeros(1,m);
faip2i=zeros(1,m);
xp2i=zeros(1,m);
yp2i=zeros(1,m);
beta2i=zeros(1,m);
Lf2i=zeros(1,m);
rou2i=zeros(1,m);
rK2i=zeros(m,n+1);
x2i=zeros(m,n+1);
y2i=zeros(m,n+1);
Tk2i=zeros(m,n+1);
alphaK2i=zeros(m,n+1);
invalphaK2i=zeros(m,n+1);
A2i=zeros(m,n+1);
I2i=zeros(m,n+1);
equivalentA2i=zeros(m,n);
equivalentI2i=zeros(m,n);
N1Bj=zeros(1,m);
O1Bj=zeros(1,m);
rp1j=zeros(1,m);
alphap1j=zeros(1,m);
invalphap1j=zeros(1,m);
faip1j=zeros(1,m);
xp1j=zeros(1,m);
yp1j=zeros(1,m);
beta1j=zeros(1,m);
Lf1j=zeros(1,m);
rou1j=zeros(1,m);
rK1j=zeros(m,n+1);
x1j=zeros(m,n+1);
y1j=zeros(m,n+1);
Tk1j=zeros(m,n+1);
alphaK1j=zeros(m,n+1);
invalphaK1j=zeros(m,n+1);
A1j=zeros(m,n+1);
I1j=zeros(m,n+1);
equivalentA1j=zeros(m,n);
equivalentI1j=zeros(m,n);
N2Bj=zeros(1,m);
O2Bj=zeros(1,m);
rp2j=zeros(1,m);
alphap2j=zeros(1,m);
invalphap2j=zeros(1,m);
faip2j=zeros(1,m);
xp2j=zeros(1,m);
yp2j=zeros(1,m);
beta2j=zeros(1,m);
Lf2j=zeros(1,m);
rou2j=zeros(1,m);
rK2j=zeros(m,n+1);
x2j=zeros(m,n+1);
y2j=zeros(m,n+1);
Tk2j=zeros(m,n+1);
alphaK2j=zeros(m,n+1);
invalphaK2j=zeros(m,n+1);
A2j=zeros(m,n+1);
I2j=zeros(m,n+1);
equivalentA2j=zeros(m,n);
equivalentI2j=zeros(m,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xM1=rf1;
xM2=rf2;
yM1=rf1*tan(360/(2*z1)*pi/180);
yM2=rf2*tan(360/(2*z2)*pi/180);
Hf1=2*yM1;
Hf2=2*yM2;
%%%%%%%%%%%%%%%%%%%第一个双齿啮合区B2C%%%%%%%%%%%%%%%%%%%%%%%
for p=1:m
%i啮合点
%第1个外齿轮啮合点
N1Bi(1,p)=N1B2+(p-1)*(B2C/m);
O1Bi(1,p)=sqrt(N1Bi(1,p)^2+rb1^2);
rp1i(1,p)=O1Bi(1,p);
alphap1i(1,p)=acos(rb1/rp1i(1,p));%弧度
invalphap1i(1,p)=tan(alphap1i(1,p))-alphap1i(1,p);
faip1i(1,p)=s1/r1-2*(invalphap1i(1,p)-invalpha);
xp1i(1,p)=rp1i(1,p)*cos(faip1i(1,p)/2); %草图中的xp
yp1i(1,p)=0.5*(s1*rp1i(1,p)/r1-2*rp1i(1,p)*(invalphap1i(1,p)-invalpha)); %草图中的yp
beta1i(1,p)=alphap1i(1,p)-faip1i(1,p)/2; %草图中的betap
Lf1i(1,p)=xp1i(1,p)-xM1-yp1i(1,p)*tan(beta1i(1,p));
rou1i(1,p)=rb1*tan(alphap1i(1,p));
%每对应一个啮合点齿轮1分成n段
for k=1:n+1
rK1i(p,k)=rf1+(rp1i(1,p)-rf1)*(k-1)/n;
%分基圆小于齿根圆和基圆大于齿根圆两种情况讨论
if rb1<rf1
%基圆小于齿根圆时整个齿都视为渐开线齿廓
alphaK1i(p,k)=acos(rb1/rK1i(p,k)); %渐开线上任意点压力角%弧度
invalphaK1i(p,k)=tan(alphaK1i(p,k))-alphaK1i(p,k);
x1i(p,k)=rK1i(p,k); %对应草图中的xk
y1i(p,k)=0.5*(s1*rK1i(p,k)/r1-2*rK1i(p,k)*(invalphaK1i(p,k)-invalpha)); %对应草图中的yk
Tk1i(p,k)=(rp1i(1,p)-rf1)/n;
else
%基圆大于齿根圆时基圆到齿根圆部分视为矩形，矩形宽为sb
alphab1i=acos(rb1/rb1);
invalphab1i=tan(alphab1i)-alphab1i;
sb1=s1*rb1/r1-2*rb1*(invalphab1i-invalpha);
Tk1i(p,k)=(rp1i(1,p)-rf1)/n;
x1i(p,k)=rK1i(p,k);
if (rK1i(p,k)<rb1)
y1i(p,k)=0.5*sb1;
else
alphaK1i(p,k)=acos(rb1/rK1i(p,k)); %渐开线上任意点压力角%弧度
invalphaK1i(p,k)=tan(alphaK1i(p,k))-alphaK1i(p,k);
y1i(p,k)=0.5*(s1*rK1i(p,k)/r1-2*rK1i(p,k).*(invalphaK1i(p,k)-invalpha));
end
end
A1i(p,k)=2*y1i(p,k)*B1;
I1i(p,k)=B1*(2*y1i(p,k)).^3/12;
end
for k=1:n
equivalentA1i(p,k)=2*(A1i(p,k)*A1i(p,k+1)/(A1i(p,k)+A1i(p,k+1)));%小段k的当量截面积
equivalentI1i(p,k)=2*(I1i(p,k)*I1i(p,k+1)/(I1i(p,k)+I1i(p,k+1)));%小段k的当量惯性矩
end
%第2个外齿轮啮合点
N2Bi(1,p)=N2B2-(p-1)*(B2C/m);
O2Bi(1,p)=sqrt(N2Bi(1,p)^2+rb2^2);
rp2i(1,p)=O2Bi(1,p);
alphap2i(1,p)=acos(rb2/rp2i(1,p));%弧度
invalphap2i(1,p)=tan(alphap2i(1,p))-alphap2i(1,p);
faip2i(1,p)=s2/r2-2*(invalphap2i(1,p)-invalpha);
xp2i(1,p)=rp2i(1,p)*cos(faip2i(1,p)/2); %草图中的xp
yp2i(1,p)=0.5*(s2*rp2i(1,p)/r2-2*rp2i(1,p)*(invalphap2i(1,p)-invalpha)); %草图中的yp
beta2i(1,p)=alphap2i(1,p)-faip2i(1,p)/2; %草图中的betap
Lf2i(1,p)=xp2i(1,p)-xM2-yp2i(1,p)*tan(beta2i(1,p));
rou2i(1,p)=rb2*tan(alphap2i(1,p));
%每对应一个啮合点齿轮1分成n段
for k=1:n+1
rK2i(p,k)=rf2+(rp2i(1,p)-rf2)*(k-1)/n;
%分基圆小于齿根圆和基圆大于齿根圆两种情况讨论
if rb2<rf2
%基圆小于齿根圆时整个齿都视为渐开线齿廓
alphaK2i(p,k)=acos(rb2/rK2i(p,k)); %渐开线上任意点压力角%弧度
invalphaK2i(p,k)=tan(alphaK2i(p,k))-alphaK2i(p,k);
x2i(p,k)=rK2i(p,k); %对应草图中的xk
y2i(p,k)=0.5*(s2*rK2i(p,k)/r2-2*rK2i(p,k)*(invalphaK2i(p,k)-invalpha)); %对应草图中的yk
Tk2i(p,k)=(rp2i(1,p)-rf2)/n;
else
%基圆大于齿根圆时基圆到齿根圆部分视为矩形，矩形宽为sb
alphab2i=acos(rb2/rb2);
invalphab2i=tan(alphab2i)-alphab2i;
sb2=s2*rb2/r2-2*rb2*(invalphab2i-invalpha);
Tk2i(p,k)=(rp2i(1,p)-rf2)/n;
x2i(p,k)=rK2i(p,k);
if (rK2i(p,k)<rb2)
y2i(p,k)=0.5*sb2;
else
alphaK2i(p,k)=acos(rb2/rK2i(p,k)); %渐开线上任意点压力角%弧度
invalphaK2i(p,k)=tan(alphaK2i(p,k))-alphaK2i(p,k);
y2i(p,k)=0.5*(s2*rK2i(p,k)/r2-2*rK2i(p,k).*(invalphaK2i(p,k)-invalpha));
e
