www.pudn.com > shichuan.rar > shichuan.m, change:2013-10-28,size:8222b


clear; 
z1=19;z2=45;%两齿轮的齿数 
m=1.75; 
b=m*z1*0.6; 
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=754;%外力 
E=2e+008;%弹性模量 
v=0.26;%泊松比 
n=100; 
step=B2C/n;%B2C为双齿啮合区 
nz1=3600;%齿轮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); 
%plot(2,2,1);plot(x,k);hold on;plot(xxx,ki2);xlabel('啮合线位移/(mm)');ylabel('线性啮合刚度k'); 
 
%plot(tt,k); 
%hold on; 
%plot(ttt,ki2); 
%hold on; 
T=0.0008772; 
 
for i=0:1; 
    A=polyfit(tt,k,2) 
    for TT=0.000048:0.0000005:0.000512; 
        K=-187835724775138*TT.*TT+80690343132.3420*TT+502073826.920299; 
        S1=6.15*0.000001*sin(2*pi*60*19*(TT+i*T)); 
        if K.*S1<0 
           P1=0; 
        else 
           P1=K.*S1;  
        end 
        plot(TT+i*T,P1,'r');xlabel('啮合时间/(s)');ylabel('激励/(N)'); 
        hold on; 
    end 
     
    for A1=0.000512:0.0000005:0.0006004; 
        A2=1898981472963.80*A1+1466425788.53; 
        S2=6.15*0.000001*sin(2*pi*60*19*(A1+i*T)); 
        if A2.*S2<0 
           P2=0; 
        else 
           P2=A2.*S2;  
        end 
        plot(A1+i*T,P2,'r'); 
        hold on; 
    end 
     
    for TTT=0.0006004:0.0000005:0.0008344; 
        B=polyfit(ttt,ki2,2); 
        KI2=-472401734368030*TTT.*TTT+659557708465.797*TTT+100570316.788086; 
        S3=6.15*0.000001*sin(2*pi*60*19*(TTT+i*T)); 
        if KI2.*S3<0 
           P3=0; 
        else 
           P3=KI2.*S3; 
        end 
        plot(TTT+i*T,P3,'r'); 
        hold on; 
    end 
     
    for A3=0.0008344:0.0000005:0.0009252; 
        A4=2020991673127.75*A3-1364307306.1; 
        S4=6.15*0.000001*sin(2*pi*60*19*(A3+i*T)); 
        if A4.*S4<0 
           P4=0; 
        else 
           P4=A4.*S4; 
        end 
        plot(A3+i*T,P4,'r'); 
    end 
end 
 
 
 
 
%subplot(2,2,3);plot(x,ki1);hold on;plot(x,kj); 
%subplot(2,2,4);plot(x,k);hold on;plot(xxx,ki2);xlabel('啮合线位移/(mm)');ylabel('线性啮合刚度k');