Weber公式计算齿轮啮合刚度 结果出错请各位指教

xiaowu

以下是小弟用Weber公式计算齿轮啮合刚度 结果为何与教材上的图形曲线不一样？
请各位高手指教，谢谢

计算公式及参数意义

[localimg=508,563]1[/localimg]

[localimg=465,403]3[/localimg]


计算结果

[localimg=562,421]2[/localimg]




计算程序：
z1=45;z2=90;%两齿轮的齿数
m=3;
b=m*8;
ha=1;
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;
rf1=r1-hf1;
rf2=r2-hf2;%齿根圆半径
alpha=20*pi/180;%分度圆压力角
invalpha=tan(alpha)-alpha;
db1=d1*cos(alpha);
db2=d2*cos(alpha);
rb1=db1/2;
rb2=db2/2;%基圆半径
da1=d1+2*ha*m;
da2=d2+2*ha*m;
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;%分度圆齿槽宽
sL1=rf1*(s/r1+2*((tan(alpha)-alpha)-(tan(alpha_f1)-alpha_f1)));
sL2=rf2*(s/r2+2*((tan(alpha)-alpha)-(tan(alpha_f2)-alpha_f2)));  %齿根圆齿厚
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;
N2C=N2B1+Pb;
B2C=B1B2-Pb;%双齿啮合区
CD=Pb-B2C;%单齿啮合区
N1C=N1B1-Pb;
FL=1000;%外力
E=20000;%弹性模量
t=linspace(0,0.46,100);
x1=rb1./cos(t).*sin(pi/2/z1-tan(t)+t+invalpha);
y1=rb1./cos(t).*cos(pi/2/z1-tan(t)+t+invalpha)-rb1*cos(pi/2/z1+invalpha);%齿轮1轮廓曲线
x2=rb2./cos(t).*sin(pi/2/z2-tan(t)+t+invalpha);
y2=rb2./cos(t).*cos(pi/2/z2-tan(t)+t+invalpha)-rb2*cos(pi/2/z2+invalpha);%齿轮2轮廓曲线
n=100;
n1=100;
step=B2C/n;%B2C为双齿啮合区
temp1=0;
temp2=0;
temp3=0;
temp4=0;
temp5=0;
temp6=0;
temp7=0;
temp8=0;
%双齿啮合区设一啮合为i点，一啮合点为j点。
for i=1:n
    x(i)=i*step;
    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);
    sitai1(i)=ai1(i)-gamai1(i);
    roui1(i)=N1Bi(i);
    hi1(i)=roui1(i)-tan(sitai1(i))*rb1;
    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);
    sitai2(i)=ai2(i)-gamai2(i);
    roui2(i)=N2Bi(i);
    hi2(i)=roui2(i)-tan(sitai2(i))*rb2;
    yLi1(i)=rb1/cos(sitai1(i))-rb1*cos(pi/2/z1+invalpha);
    yLi2(i)=rb2/cos(sitai2(i))-rb2*cos(pi/2/z2+invalpha);
    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);
    sitaj1(i)=aj1(i)-gamaj1(i);
    rouj1(i)=N1Bj(i);
    hj1(i)=rouj1(i)-tan(sitaj1(i))*rb1;
    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);
    sitaj2(i)=aj2(i)-gamaj2(i);
    rouj2(i)=N2Bj(i);
    hj2(i)=N2Bj(i)-tan(sitaj2(i))*rb2;
    yLj1(i)=rb1/cos(sitaj1(i))-rb1*cos(pi/2/z1+invalpha);
    yLj2(i)=rb2/cos(sitaj2(i))-rb2*cos(pi/2/z2+invalpha);
    step1=yLi1(i)/n1;
    step2=yLi2(i)/n1;
    step3=yLj1(i)/n1;
    step4=yLj2(i)/n1;
for j=1:n1
    Y1(j)=j*step1;
    X1(j)=interp1(y1,x1,Y1(j));
    temp1=temp1+1/(2*X1(j))*step1;%积分运算
    temp2=temp2+(yLi1(i)-Y1(j))^2/(8*X1(j)^3)*step1;%积分运算
    Y2(j)=j*step2;
    X2(j)=interp1(y2,x2,Y2(j));
    temp3=temp3+1/(2*X2(j))*step2;%积分运算
    temp4=temp4+(yLi2(i)-Y2(j))^2/(8*X2(j)^3)*step2;%积分运算
    Y1(j)=j*step3;
    X1(j)=interp1(y1,x1,Y1(j));
    temp5=temp5+1/(2*X1(j))*step3;%积分运算
    temp6=temp6+(yLj1(i)-Y1(j))^2/(8*X1(j)^3)*step3;%积分运算
    Y2(j)=j*step4;
    X2(j)=interp1(y2,x2,Y2(j));
    temp7=temp7+1/(2*X2(j))*step4;%积分运算
    temp8=temp8+(yLj2(i)-Y2(j))^2/(8*X2(j)^3)*step4;%积分运算
end
    sigmaHi(i)=0.5793*FL*(log(E*b*hi1(i)*hi2(i)*(roui1(i)+roui2(i))/0.5793*FL/roui1(i)/roui2(i))-0.4286)/E/b;
    sigmaTi1(i)=FL*cos(sitai1(i))^2*(12*temp2+(3+tan(sitai1(i))^2)*temp1)/E/b;
    sigmaTi2(i)=FL*cos(sitai2(i))^2*(12*temp4+(3+tan(sitai2(i))^2)*temp3)/E/b;
    sigmaAi1(i)=FL*cos(sitai1(i))^2*(5.214*(yLi1(i)^2/sL1)^2+1.04*yLi1(i)/sL1+1.39*(1+tan(sitai1(i))^2)+0.4485*tan(sitai1(i))^4)/E/b;
    sigmaAi2(i)=FL*cos(sitai2(i))^2*(5.214*(yLi2(i)^2/sL2)^2+1.04*yLi2(i)/sL2+1.39*(1+tan(sitai2(i))^2)+0.4485*tan(sitai2(i))^4)/E/b;
    kci(i)=FL/(sigmaHi(i)+sigmaTi1(i)+sigmaTi2(i)+sigmaAi1(i)+sigmaAi2(i));%双齿啮合时啮合点i的刚度
    sigmaHj(i)=0.5793*FL*(log(E*b*hj1(i)*hj2(i)*(rouj1(i)+rouj2(i))/0.5793*FL/rouj1(i)/rouj2(i))-0.4286)/E/b;
    sigmaTj1(i)=FL*cos(sitaj1(i))^2*(12*temp6+(3+tan(sitaj1(i))^2)*temp5)/E/b;
    sigmaTj2(i)=FL*cos(sitaj2(i))^2*(12*temp8+(3+tan(sitaj2(i))^2)*temp7)/E/b;
    sigmaAj1(i)=FL*cos(sitaj1(i))^2*(5.214*(yLj1(i)^2/sL1)^2+1.04*yLj1(i)/sL1+1.39*(1+tan(sitaj1(i))^2)+0.4485*tan(sitaj1(i))^4)/E/b;
    sigmaAj2(i)=FL*cos(sitaj2(i))^2*(5.214*(yLj2(i)^2/sL2)^2+1.04*yLj2(i)/sL2+1.39*(1+tan(sitaj2(i))^2)+0.4485*tan(sitaj2(i))^4)/E/b;
    kcj(i)=FL/(sigmaHj(i)+sigmaTj1(i)+sigmaTj2(i)+sigmaAj1(i)+sigmaAj2(i));%双齿啮合时啮合点j的刚度
    kc(i)=kci(i)+kcj(i);
end
step5=CD/n;%CD为单齿啮合区
temp1=0;
temp2=0;
temp3=0;
temp4=0;
for i=1:n
    xxx(i)=B2C+i*step5;
    N1Bi(i)=N1C+i*step5;
    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);
    sitai1(i)=ai1(i)-gamai1(i);
    roui1(i)=N1Bi(i);
    hi1(i)=roui1(i)-tan(sitai1(i))*rb1;
    N2Bi(i)=N2C-i*step5;
    O2Bi(i)=sqrt(N2Bi(i)*N2Bi(i)+rb2^2);
    ai2(i)=acos(rb2/O2Bi(i));
    gamai2(i)=pi/2/z2+invalpha-tan(ai2(i))+ai2(i);
    sitai2(i)=ai2(i)-gamai2(i);
    roui2(i)=N2Bi(i);
    hi2(i)=roui2(i)-tan(sitai2(i))*rb2;
    yLi1(i)=rb1/cos(sitai1(i))-rb1*cos(pi/2/z1+invalpha);
    yLi2(i)=rb2/cos(sitai2(i))-rb2*cos(pi/2/z2+invalpha);
    step1=yLi1(i)/n1;
    step2=yLi2(i)/n1;
for j=1:n1
    Y1(j)=j*step1;
    X1(j)=interp1(y1,x1,Y1(j));
    temp1=temp1+1/(2*X1(j))*step1;%积分运算
    temp2=temp2+(yLi1(i)-Y1(j))^2/(8*X1(j)^3)*step1;%积分运算
    Y2(j)=j*step2;
    X2(j)=interp1(y2,x2,Y2(j));
    temp3=temp3+1/(2*X2(j))*step2;%积分运算
    temp4=temp4+(yLi2(i)-Y2(j))^2/(8*X2(j)^3)*step2;%积分运算
end
     sigmaH(i)=0.5793*FL*(log(E*b*hi1(i)*hi2(i)*(roui1(i)+roui2(i))/0.5793*FL/roui1(i)/roui2(i))-0.4286)/E/b;
     sigmaT1(i)=FL*cos(sitai1(i))^2*(12*temp2+(3+tan(sitai1(i))^2)*temp1)/E/b;
     sigmaT2(i)=FL*cos(sitai2(i))^2*(12*temp4+(3+tan(sitai2(i))^2)*temp3)/E/b;
     sigmaA1(i)=FL*cos(sitai1(i))^2*(5.214*(yLi1(i)^2/sL1)^2+1.04*yLi1(i)/sL1+1.39*(1+tan(sitai1(i))^2)+0.4485*tan(sitai1(i))^4)/E/b;
     sigmaA2(i)=FL*cos(sitai2(i))^2*(5.214*(yLi2(i)^2/sL2)^2+1.04*yLi2(i)/sL2+1.39*(1+tan(sitai2(i))^2)+0.4485*tan(sitai2(i))^4)/E/b;
     kc1(i)=FL./(sigmaH(i)+sigmaT1(i)+sigmaT2(i)+sigmaA1(i)+sigmaA2(i));
end
plot(x,kc);
hold on;
plot(xxx,kc1);
hold on;
plot(xx,kc);
xlabel('啮合线位移/(mm)'),ylabel('线性啮合刚度k/(N/m)');

xiaowu

计算公式与计算结果

糖窝窝

你很有毅力，我实践了很久的，最后还是放弃weber，用石川法计算了。加油，希望你能解决问题。

xiaowu

weber其实也不是太复杂，这个结果不符合实际，还是编程有问题。

tao_seu

正在关注这方面的理论...

569340337

课本上的正确答案是什么样的？你贴出来，我给你检查一下

569340337

我的569340337，我现在也做这方面的研究，可以加我一起探讨一下

powershift

齿轮啮合内部动态激励数值模拟 一文中提到建立了齿面任意啮合位置接触有限元模型的算法，并开发了相应的程序，不知道哪位高人能提供一下这样的程序吗？让我们这些后来人也学习学习吧，

fawcgzmg

没考虑螺旋角，斜齿轮，只是直齿轮啊

wxyygn

"齿轮啮合内部动态激励数值模拟 一文中提到建立了齿面任意啮合位置接触有限元模型的算法，并开发了相应的程序，不知道哪位高人能提供一下这样的程序吗？让我们这些后来人也学习学习吧"，我也需要，哪位能共享下，非常感谢

donghaofane

还不错 我是用LTCA 计算的

meiqiuchixiaoji

我正在研究齿轮这方面的知识，我刚入门，希望一块探讨下 2417433971

yangbosen

糖窝窝 发表于 2010-8-6 15:51
**** 作者被禁止或删除 内容自动屏蔽 ****

跪求程序呀，我是新手，能不能加您qq或者微信