石川法求齿轮啮合刚度（仅供参考）
clear;z1=45;z2=90;%两齿轮的齿数
m=3;
b=m*8;
ha=m;
c=0.25*m;
d1=m*z1;d2=m*z2;
r1=m*z1/2;r2=m*z2/2;%分度圆半径
hf1=1.25*m;
hf2=1.25*m;
alpha=20*pi/180;%分度圆压力角
invalpha=tan(alpha)-alpha;
db1=d1*cos(alpha);
db2=d2*cos(alpha);
rb1=db1/2;
rb2=db2/2;%基圆半径
rf1=r1-hf1;
rf2=r2-hf2;%齿根圆半径
da1=d1+2*ha;
da2=d2+2*ha;
ra1=da1/2;
ra2=da2/2;%齿顶圆半径
alpha_a1=acos(rb1/ra1);
alpha_a2=acos(rb2/ra2);%齿顶圆压力角
alpha_f1=acos(rb1/rf1);
alpha_f2=acos(rb2/rf2);%齿根圆压力角
s=pi*m/2;%分度圆弧齿厚
e=s;%分度圆齿槽宽
sk1=ra1*(s/r1+2*((tan(alpha)-alpha)-(tan(alpha_a1)-alpha_a1)));
sk2=ra2*(s/r2+2*((tan(alpha)-alpha)-(tan(alpha_a2)-alpha_a2)));%齿顶圆齿厚
PB1=r1*cos(alpha)*(tan(alpha_a1)-tan(alpha));
PB2=r2*cos(alpha)*(tan(alpha_a2)-tan(alpha));
B1B2=PB1+PB2;
Pb=pi*m*cos(alpha);%基圆齿距
Epsilona=B1B2/Pb;
N1B1=sqrt(ra1^2-rb1^2);
N1B2=N1B1-B1B2;
N2B2=sqrt(ra2^2-rb2^2);
N2B1=N2B2-B1B2;
rF2=sqrt(rb2^2+N2B1^2);
rF1=sqrt(rb1^2+N1B2^2);%有效齿根圆半径
alpha_F1=acos(rb1/rF1);
alpha_F2=acos(rb2/rF2);
sf1=2*rF1*sin(pi/2/z1+invalpha-tan(alpha_F1)+alpha_F1);
sf2=2*rF2*sin(pi/2/z2+invalpha-tan(alpha_F2)+alpha_F2);
h1=sqrt(ra1^2-(sk1/2)^2)-sqrt(rf1^2-(sf1/2)^2);
h2=sqrt(ra2^2-(sk2/2)^2)-sqrt(rf2^2-(sf2/2)^2);
hr1=sqrt(rF1^2-(sf1/2)^2)-sqrt(rf1^2-(sf1/2)^2);
hr2=sqrt(rF2^2-(sf2/2)^2)-sqrt(rf2^2-(sf2/2)^2);
hi1=(h1*sf1-hr1*sk1)/(sf1-sk1);
hi2=(h2*sf2-hr2*sk2)/(sf2-sk2);
N2C=N2B1+Pb;
B2C=B1B2-Pb;%双齿啮合区
CD=Pb-B2C;%单齿啮合区
N1C=N1B1-Pb;
Fn=1000;%外力
E=2e+008;%弹性模量
v=0.26;%泊松比
n=100;
step=B2C/n;%B2C为双齿啮合区
nz1=3000;%齿轮1转速r/min
Tz=60/z1/nz1;%循环的周期
t1=B2C/Pb*Tz;
t2=CD/Pb*Tz;
step2=t1/n;
step4=t2/n;
%双齿啮合区设一啮合为i点，一啮合点为j点。
for i=1:n
x(i)=i*step;
tt(i)=i*step2;
xx(i)=Pb+i*step;
N1Bi(i)=N1B2+i*step;%双齿啮合区i啮合点公式中具体参数的计算
O1Bi(i)=sqrt(N1Bi(i)*N1Bi(i)+rb1^2);
ai1(i)=acos(rb1/O1Bi(i));
gamai1(i)=pi/2/z1+invalpha-tan(ai1(i))+ai1(i);
miui1(i)=ai1(i)-gamai1(i);
rxi1(i)=O1Bi(i);
hxi1(i)=rxi1(i)*cos(gamai1(i))-sqrt(rf1^2-(sf1/2)^2);
N2Bi(i)=N2B2-i*step;
O2Bi(i)=sqrt(N2Bi(i)^2+rb2^2);
ai2(i)=acos(rb2/O2Bi(i));
gamai2(i)=pi/2/z2+invalpha-tan(ai2(i))+ai2(i);
miui2(i)=ai2(i)-gamai2(i);
rxi2(i)=O2Bi(i);
hxi2(i)=rxi2(i)*cos(gamai2(i))-sqrt(rf2^2-(sf2/2)^2);
N1Bj(i)=N1Bi(i)+Pb;%双齿啮合区j啮合点公式中具体参数的计算
O1Bj(i)=sqrt(N1Bj(i)^2+rb1^2);
aj1(i)=acos(rb1/O1Bj(i));
gamaj1(i)=pi/2/z1+tan(alpha)-alpha-tan(aj1(i))+aj1(i);
miuj1(i)=aj1(i)-gamaj1(i);
rxj1(i)=O1Bj(i);
hxj1(i)=rxj1(i)*cos(gamaj1(i))-sqrt(rf1^2-(sf1/2)^2);
N2Bj(i)=N2Bi(i)-Pb;
O2Bj(i)=sqrt(N2Bj(i)^2+rb2^2);
aj2(i)=acos(rb2/O2Bj(i));
gamaj2(i)=pi/2/z2+tan(alpha)-alpha-tan(aj2(i))+aj2(i);
miuj2(i)=aj2(i)-gamaj2(i);
rxj2(i)=O2Bj(i);
hxj2(i)=rxj2(i)*cos(gamaj2(i))-sqrt(rf2^2-(sf2/2)^2);
sigmabri1(i)=12*Fn*cos(miui1(i))^2*(hxi1(i)*hr1*(hxi1(i)-hr1)+hxi1(i)^3/3)/b/E/sf1^3;
sigmabri2(i)=12*Fn*cos(miui2(i))^2*(hxi2(i)*hr2*(hxi2(i)-hr2)+hxi2(i)^3/3)/b/E/sf2^3;
sigmabti1(i)=6*Fn*cos(miui1(i))^2*((hi1-hxi1(i))/(hi1-hr1)*(4-(hi1-hxi1(i))/(hi1-hr1))-2*log((hi1-hxi1(i))/(hi1-hr1))-3)*(hi1-hr1)^3/b/E/sf1^3;
sigmabti2(i)=6*Fn*cos(miui2(i))^2*((hi2-hxi2(i))/(hi2-hr2)*(4-(hi2-hxi2(i))/(hi2-hr2))-2*log((hi2-hxi2(i))/(hi2-hr2))-3)*(hi2-hr2)^3/b/E/sf2^3;
sigmasi1(i)=2*(1+v)*Fn*cos(miui1(i))^2*(hr1+(hi1-hr1)*log((hi1-hr1)/(hi1-hxi1(i))))/b/E/sf1;
sigmasi2(i)=2*(1+v)*Fn*cos(miui2(i))^2*(hr2+(hi2-hr2)*log((hi2-hr2)/(hi2-hxi2(i))))/b/E/sf2;
sigmagi1(i)=24*Fn*hxi1(i)*cos(miui1(i))^2/pi/b/E/sf1^2;
sigmagi2(i)=24*Fn*hxi2(i)*cos(miui2(i))^2/pi/b/E/sf2^2;
sigmap=4*Fn*(1-v^2)/pi/b/E;
sigmabrj1(i)=12*Fn*cos(miuj1(i))^2*(hxj1(i)*hr1*(hxj1(i)-hr1)+hxj1(i)^3/3)/b/E/sf1^3;
sigmabrj2(i)=12*Fn*cos(miuj2(i))^2*(hxj2(i)*hr2*(hxj2(i)-hr2)+hxj2(i)^3/3)/b/E/sf2^3;
sigmabtj1(i)=6*Fn*cos(miuj1(i))^2*((hi1-hxj1(i))/(hi1-hr1)*(4-(hi1-hxj1(i))/(hi1-hr1))-2*log((hi1-hxj1(i))/(hi1-hr1))-3)*(hi1-hr1)^3/b/E/sf1^3;
sigmabtj2(i)=6*Fn*cos(miuj2(i))^2*((hi2-hxj2(i))/(hi2-hr2)*(4-(hi2-hxj2(i))/(hi2-hr2))-2*log((hi2-hxj2(i))/(hi2-hr2))-3)*(hi2-hr2)^3/b/E/sf2^3;
sigmasj1(i)=2*(1+v)*Fn*cos(miuj1(i))^2*(hr1+(hi1-hr1)*log((hi1-hr1)/(hi1-hxj1(i))))/b/E/sf1;
sigmasj2(i)=2*(1+v)*Fn*cos(miuj2(i))^2*(hr2+(hi2-hr2)*log((hi2-hr2)/(hi2-hxj2(i))))/b/E/sf2;
sigmagj1(i)=24*Fn*hxj1(i)*cos(miuj1(i))^2/pi/b/E/sf1^2;
sigmagj2(i)=24*Fn*hxj2(i)*cos(miuj2(i))^2/pi/b/E/sf2^2;
ki1(i)=Fn/(sigmabri1(i)+sigmabri2(i)+sigmabti1(i)+sigmabti2(i)+sigmasi1(i)+sigmasi2(i)+sigmagi1(i)+sigmagi2(i)+sigmap);%双齿啮合时i点刚度
kj(i)=Fn/(sigmabrj1(i)+sigmabrj2(i)+sigmabtj1(i)+sigmabtj2(i)+sigmasj1(i)+sigmasj2(i)+sigmagj1(i)+sigmagj2(i)+sigmap);%双齿啮合时j点刚度
k(i)=ki1(i)+kj(i);
end
step3=CD/n;%CD为单齿啮合区
for i=1:n
xxx(i)=B2C+i*step3;
ttt(i)=t1+i*step4;
N1Bi(i)=N1C+i*step3;
O1Bi(i)=sqrt(N1Bi(i)^2+rb1^2);
ai1(i)=acos(rb1/O1Bi(i));
gamai1(i)=pi/2/z1+invalpha-tan(ai1(i))+ai1(i);
miui1(i)=ai1(i)-gamai1(i);
rxi1(i)=O1Bi(i);
hxi1(i)=rxi1(i)*cos(gamai1(i))-sqrt(rf1^2-(sf1/2)^2);
N2Bi(i)=N2C-i*step3;
O2Bi(i)=sqrt(N2Bi(i)^2+rb2^2);
ai2(i)=acos(rb2/O2Bi(i));
gamai2(i)=pi/2/z2+invalpha-tan(ai2(i))+ai2(i);
miui2(i)=ai2(i)-gamai2(i);
rxi2(i)=O2Bi(i);
hxi2(i)=rxi2(i)*cos(gamai2(i))-sqrt(rf2^2-(sf2/2)^2);
sigmabri1(i)=12*Fn*cos(miui1(i))^2*(hxi1(i)*hr1*(hxi1(i)-hr1)+hxi1(i)^3/3)/b/E/sf1^3;
sigmabri2(i)=12*Fn*cos(miui2(i))^2*(hxi2(i)*hr2*(hxi2(i)-hr2)+hxi2(i)^3/3)/b/E/sf2^3;
sigmabti1(i)=6*Fn*cos(miui1(i))^2*((hi1-hxi1(i))/(hi1-hr1)*(4-(hi1-hxi1(i))/(hi1-hr1))-2*log((hi1-hxi1(i))/(hi1-hr1))-3)*(hi1-hr1)^3/b/E/sf1^3;
sigmabti2(i)=6*Fn*cos(miui2(i))^2*((hi2-hxi2(i))/(hi2-hr2)*(4-(hi2-hxi2(i))/(hi2-hr2))-2*log((hi2-hxi2(i))/(hi2-hr2))-3)*(hi2-hr2)^3/b/E/sf2^3;
sigmasi1(i)=2*(1+v)*Fn*cos(miui1(i))^2*(hr1+(hi1-hr1)*log((hi1-hr1)/(hi1-hxi1(i))))/b/E/sf1;
sigmasi2(i)=2*(1+v)*Fn*cos(miui2(i))^2*(hr2+(hi2-hr2)*log((hi2-hr2)/(hi2-hxi2(i))))/b/E/sf2;
sigmagi1(i)=24*Fn*hxi1(i)*cos(miui1(i))^2/pi/b/E/sf1^2;
sigmagi2(i)=24*Fn*hxi2(i)*cos(miui2(i))^2/pi/b/E/sf2^2;
sigmap=4*Fn*(1-v^2)/pi/b/E;
ki2(i)=Fn/(sigmabri1(i)+sigmabri2(i)+sigmabti1(i)+sigmabti2(i)+sigmasi1(i)+sigmasi2(i)+sigmagi1(i)+sigmagi2(i)+sigmap);
end
o1=polyfit(tt,k,3);
o2=polyfit(ttt,ki2,3);
o3=polyfit(tt,ki1,3);
o4=polyfit(tt,kj,3);
subplot(2,2,1),plot(x,k),hold on,plot(xxx,ki2),xlabel('啮合线位移/(mm)'),ylabel('线性啮合刚度k')
subplot(2,2,2),plot(tt,k),hold on,plot(ttt,ki2),xlabel('啮合时间/(s)'),ylabel('线性啮合刚度k')
subplot(2,2,3),plot(x,ki1),hold on,plot(x,kj)，figure
plot(x,k),hold on,plot(xxx,ki2),xlabel('啮合线位移/(mm)'),ylabel('线性啮合刚度k') 谢谢楼主分享！！！！！！ 忘了说一件事，弹性模量E本来该取E=2e+011，由于齿轮的变形是以mm为单位，所以数量级做了调整，最后求得的刚度单位是N/m，不是N/mm 谢谢，楼主分享，正研究这个呢？学习了 不错，赞一下搂主。 LZ可否给一份石川法的文字资料.
xiaoke24@yahoo.com.cn
Q：20026347
thx~ sf1=2*rF1*sin(pi/2/z1+invalpha-tan(alpha_F1)+alpha_F1);
sf2=2*rF2*sin(pi/2/z2+invalpha-tan(alpha_F2)+alpha_F2); 什么意思啊？？？
sf1=2*rF1*sin(pi/2/z1+invalpha-tan(alpha_F1)+alpha_F1);
sf2=2*rF2*sin(pi/2/z2+invalpha-tan(alpha_F2)+alpha_F2); 回复 8 # 569340337 的帖子
朱秋玲的《齿轮系统动力学分析及计算机仿真》里有公式说明。
顺便说一下公式里的b是齿宽，我取错了，齿宽应该为b=d1(轮1分度圆直径）*齿宽系数（取0.4~0.9）。
程序编的匆忙，大家仅作参考 你这个是什么程序？ 斜齿轮怎么考虑
谢谢楼主分享 请问楼主，是否在计算双齿啮合区时有误？我认为你少计算了一半 太感动了 现在正在学这个呢不知道这位高手有没有解齿轮系统的方程的程序啊？？ 回复 13 # nrg 的帖子
双齿啮合区有两个啮合点，我分别设为i点，j点，两点都计算了