四阶变步长龙格库塔 齿轮系统动力学

clear;
z1=10;
z2=10;
m=5;
d1=m*z1;
d2=m*z2;
R1=d1/2;
R2=d2/2;
S=pi*m/2;
ha=m;
hf=1.25*m
h=2.25*m;
B=0.4*d1;
Alpha=20*pi/180;
Rb1=R1*cos（Alpha）;
Rb2=R2*cos（Alpha）;
Ra1=R1+ha;
Ra2=R2+ha;
Rf1=R1-hf;
Rf2=R2-hf;
Pb=pi*m*cos（Alpha）;	
thi0=tan（20*pi/180）-Alpha;
n=100
for i=1:n
    
    
r1k（i）=Rb1+（i-1）*（Ra1-Rb1）/（n-1）;
Ak1（i）=acos（Rb1/r1k（i））;
thi_k1（i）=tan（Ak1（i））-Ak1（i）;
phi1（i）=S/（2*R1）-（thi_k1（i）-thi0）;
x1（i）=r1k（i）*cos（phi1（i））;
y1（i）=r1k（i）*sin（phi1（i））;
S1（i）=r1k（i）*phi1（i）*2;
r2k（i）=Rb2+（i-1）*（Ra2-Rb2）/（n-1）;
Ak2（i）=acos（Rb2/r2k（i））;
thi_k2（i）=tan（Ak2（i））-Ak2（i）;
phi2（i）=S/（2*R2）-（thi_k2（i）-thi0）;
x2（i）=r2k（i）*cos（phi2（i））;
y2（i）=r2k（i）*sin（phi2（i））;
S2（i）=r2k（i）*phi2（i）*2;
end
figure（‘color‘[1 1 1]）;
figure（1）
plot（r1kS1‘r‘）
xlabel（‘轮齿上任意圆周/mm‘）
ylabel（‘齿厚/mm‘）
figure（‘color‘[1 1 1]）;
figure（2）
plot（x2y2）;hold on
plot（r2kS2‘r‘）
A_a1=acos（Rb1/Ra1）;
A_a2=acos（Rb2/Ra2）;
epsilon=（z1*（tan（A_a1）-tan（Alpha））+z2*（tan（A_a2）-tan（Alpha）））/（2*pi）;
n1=500;
Tz=60/（z1*n1）;
t1=Tz*（epsilon-1）;
t2=Tz*（2-epsilon）;
P=5500;
T=P/（2*pi*n1/60）;
Fn=2*T/（d1*cos（Alpha））;
Ft=Fn*cos（Alpha）;
m1=2;
ratio=z2/z1;
B1B2=pi*m*epsilon*cos（Alpha）;
PN2=Rb2*tan（Alpha）;
PN1=Rb1*tan（Alpha）;
N2B2=Rb2*tan（A_a2）;
N1B1=Rb1*tan（A_a1）;
PB2=N2B2-PN2;
PB1=N1B1-PN1;
N1B2=PN1-PB2;
N2B1=PN2-PB1;
CB2=B1B2-Pb;
Rs1=sqrt（Rb1^2+N1B2^2）;
Rs2=sqrt（Rb2^2+N2B1^2）;
for i=1:n
    x（i）=（i-1）*B1B2/（n-1）;
    if x（i）        SN（i）=PB2-x（i）;
    else
        SN（i）=x（i）-PB2;
    end
    um（i）=pi*n1/（30*ratio）*（R2*sin（Alpha）+SN（i）/2*（1-ratio））;
    Ru1（i）=N1B2+x（i）;
    Ru2（i）=N2B2-x（i）;
    Ru（i）=Ru1（i）* Ru2（i）/（Ru1（i）+ Ru2（i））;
end
n2=100;
E=200e9;
I=B*S^3/12;
for i=1:n
    %  x（i）=（i-1）*B1B2/（n2-1）;
     Rk1（i）=sqrt（Rb1^2+（N1B2+x（i））^2）;
     Rk2（i）=sqrt（Rb2^2+（N2B2+x（i））^2）;
     h1（i）=Rk1（i）-Rf1+0.25*m;
     h2（i）=Rk2（i）-Rf2+0.25*m;
     kk10（i）=3*E*I/h1（i）^3;
     kk20（i）=3*E*I/h2（i）^3;
     K_t1（i）=kk10（i）*kk20（i）/（kk10（i）+kk20（i））;
     K_t（i）= K_t1（i）;
     if x（i）         kk=i
     end
     K_t2（i-kk+1）=K_t1（i）;
     F_L（i）=Fn;
     if x（i）         kk1=i;
     end
end
% for i=1:kk
%