ToqurVibratory.rar_ode45 齿轮_动力学_齿轮_齿轮动力_齿轮动力学‘

function dy=ToqurVibratory(t,y) dy=zeros(24,2); t %齿轮轴承质量KG mb1=32; mb2=27; mb3=43; mb4=97; m1=21; m2=82; m3=221;g=9.8; %主动轮转速 n1=12; km=7.8e8; %齿轮基础参数 单位mm z1=20; z2=48; mn=10; n2=n1*z1/z2; e0=2.0*10^-5; er=3.0*10^-5; p1=1.0*10^-5; p2=p1; b=1.0*10^-5;p3=1.0*10^-5;%偏心距m %齿轮角度参数% betab=βb基圆螺旋角 alphal=120*pi/180; alphat=15.4*pi/180; betaa=15*pi/180; betab=13.8*pi/180; r1=0.001*mn*z1*cos(alphat); r2=0.001*mn*z2*cos(alphat); %转动惯量kg.m2 Jd=9.7; Jl=10.2; J1=0.087; J2=8.1; %刚度系数 kt1=1.2e8; kt2=5.1e9; ksx1=0.71e8; ksy1=0.71e8; ksz1=1.3e8; ksx2=7.3e8; ksy2=.3e8; ksz2=9.4e8; ksx3=2.1e9; ksy3=2.1e9;ksz3=3.5e9; ksx4=2.3e9; ksy4=2.3e9; ksz4=3.2e9; ksx5=5.4e9; ksy5=5.4e9; kbx1=0.41e8; kbx2=0.31e8; kbx3=1.1e8; kbx4=1.7e8; kby1=0.28e8; kby2=0.16e8; kby3=0.4e8; kby4=0.8e8; kbz1=0.13e8; kbz2=0; kbz3=0; kbz4=1.2e8; %阻尼系数 ct1=2*0.07*((kt1/(1/Jd+1/J1))^0.5); ct2=2*0.07*((kt2/(1/J2+1/Jl))^0.5); csx1=2*0.07*((ksx1/(1/mb1+1/m1))^0.5); csy1=2*0.07*((ksy1/(1/mb1+1/m1))^0.5); csz1=2*0.07*((ksz1/(1/mb1+1/m1))^0.5); csx2=2*0.07*((ksx2/(1/mb2+1/m1))^0.5); csy2=2*0.07*((ksy2/(1/mb2+1/m1))^0.5); csz2=2*0.07*((ksz2/(1/mb2+1/m1))^0.5); csx3=2*0.07*((ksx3/(1/mb3+1/m2))^0.5); csy3=2*0.07*((ksy3/(1/mb3+1/m2))^0.5); csz3=2*0.07*((ksz3/(1/mb3+1/m2))^0.5); csx4=2*0.07*((ksx4/(1/mb3+1/m2))^0.5); csy4=2*0.07*((ksy4/(1/mb3+1/m2))^0.5); csz4=2*0.07*((ksz4/(1/mb3+1/m2))^0.5); cbx1=2*0.01*((kbx1*mb1*m1/(mb1+m1))^0.5); cby1=2*0.01*((kby1*mb1*m1/(mb1+m1))^0.5); cbz1=2*0.01*((kbz1*mb1*m1/(mb1+m1))^0.5); cbx2=2*0.01*((kbx2*mb2*m1/(mb2+m1))^0.5); cby2=2*0.01*((kby2*mb2*m1/(mb2+m1))^0.5); cbz2=2*0.01*((kbz2*mb2*m1/(mb2+m1))^0.5); cbx3=2*0.01*((kbx3*mb3*m2/(mb3+m2))^0.5); cby3=2*0.01*((kby3*mb3*m2/(mb3+m2))^0.5); cbz3=2*0.01*((kbz3*mb3*m2/(mb3+m2))^0.5); cbx4=2*0.01*((kbx4*mb4*m2/(mb4+m2))^0.5); cby4=2*0.01*((kby4*mb4*m2/(mb4+m2))^0.5); cbz4=2*0.01*((kbz4*mb4*m2/(mb4+m2))^0.5); %角度参数ω=omm φ=phi α=alp omm1=2*pi*n1; omm2=2*pi*n2; ommm=2*pi*n1*z1/60; %输入轴扭矩N.m Tdm=320; Tgm=400; Tdr=200; Tgr=200; phid=omm1*t+y(1); phig=omm2*t+y(43); phi1=omm1*t+y(15); phi2=omm2*t+y(35); Td=Tdm+Tdr*sin(omm1*t+phid); Te=Tgm+Tgr*sin(omm2*t+phig); %齿轮力计算 mm=J1*J2/(J2*(r1^2)+J1*(r2^2)); cm=2*0.1*((mm*km)^0.5);%km为平均啮合刚度 et=e0+er*sin(ommm*t);%公式4 delta=(r1*y(17)-r2*y(35))+(y(9)+p1*cos(phi1)-(y(29)+p2*cos(phi2)))*cos(alphal-alphat)+... (y(11)-p1*sin(phi1)-(y(31)+p2*sin(phi2)))*sin(alphal-alphat)+(y(13)-y(33))*tan(betaa)-et;%公式3 deltat=(r1*y(18)-r2*y(36))+(y(10)+p1*cos(phi1)-(y(30)+p2*cos(phi2)))*cos(alphal-alphat)+... (y(12)-p1*sin(phi1)-(y(32)+p2*sin(phi2)))*sin(alphal-alphat)+(y(14)-y(34))*tan(betaa)-et; F=km*delta+cm*deltat;%公式5 Fx=F*cos(alphal-alphat); %公式13 Fy=F*sin(alphal-alphat); Fz=F*tan(betaa); %轴承力计算 %轴承参数 Rb1=0.25; rb1=0.2; rb2=0.1; Rb2=0.125; Db1=0.05; Db2=0.025;%轴承参数 Nb1=14; Nb2=10; Pd1=0.25e-6; Pd2=0.5e-6; fi=0.525; f0=0.515; dm1=0.5*(rb1+Rb1); A01=(f0+fi-1)*Db1; dm2=0.5*(rb2+Rb2); A02=(f0+fi-1)*Db2; ommi1=omm1; ommi2=omm2; ommcage1=ommi1*rb1/(rb1+Rb1); ommcage2=ommi2*rb2/(rb2+Rb2); alpha01=acos(1-Pd1/(2*A01)); alpha02=acos(1-Pd2/(2*A02)); Rj1=0.5*dm1+(fi*Db1-0.5*Db1)*cos(alpha01); Rj2=0.5*dm2+(fi*Db2-0.5*Db2)*cos(alpha02); gam1=Db1*cos(alpha01)/dm1; gam2=Db2*cos(alpha02)/dm2; Fpi1=(1/fi+2*gam1/(1-gam1))/(4-1/fi+2*gam1/(1-gam1));%公式30 Fpi2=(1/fi+2*gam2/(1-gam2))/(4-1/fi+2*gam2/(1-gam2));%公式30 deltai1=-327.6145+1883.338*Fpi1-3798.1121*Fpi1^2+3269.6154*Fpi1^3-1026.96*Fpi1^4;%公式29 deltai2=-327.6145+1883.338*Fpi2-3798.1121*Fpi2^2+3269.6154*Fpi2^3-1026.96*Fpi2^4;%公式29 Fp01=(1/f0-2*gam1/(1-gam1))/(4-1/f0-2*gam1/(1-gam1));%公式30 Fp02=(1/f0-2*gam2/(1-gam2))/(4-1/f0-2*gam2/(1-gam2));%公式30 delta01=-327.6145+1883.338*Fp01-3798.1121*Fp01^2+3269.6154*Fp01^3-1026.96*Fp01^4;%公式29 delta02=-327.6145+1883.338*Fp02-3798.1121*Fp02^2+3269.6154*Fp02^3-1026.96*Fp02^4;%公式29 Sigp01=(1/Db1)*(4-1/f0+2*gam1/(1-gam1));%公式28 Sigp02=(1/Db2)*(4-1/f0+2*gam2/(1-gam2));%公式28 Sigpi1=(1/Db1)*(4-1/fi+2*gam1/(1-gam1));%公式28 Sigpi2=(1/Db2)*(4-1/fi+2*gam2/(1-gam2));%公式28 ki1=2.15*(10^5)*(Sigpi1^(-0.5))*(deltai1^(-1.5));%公式27 ki2=2.15*(10^5)*(Sigpi2^(-0.5))*(deltai2^(-1.5));%公式27 k01=2.15*(10^5)*(Sigp01^(-0.5))*(delta01^(-1.5));%公式27 k02=2.15*(10^5)*(Sigp02^(-0.5))*(delta02^(-1.5));%公式27 Kc1=(1/((1/ki1)^1.5+(1/k01)^1.5))^1.5;%公式26 Kc2=(1/((1/ki2)^1.5+(1/k02)^1.5))^1.5;%公式26 for i=1:1:Nb1 thetai1=ommcage1*t+2*pi*(i-1)/Nb1; alpha10=atan((A01*sin(alpha01)+y(7)+Rj1*y(15)*cos(thetai1))/(A01*cos(alpha01)+y(3)*cos(thetai1)+y(5)*sin(thetai1)));%公式19 A10=((A01*sin(alpha01)+y(7)+Rj1*y(15)*cos(thetai1))^2+(A01*cos(alpha01)+y(3)*cos(thetai1)+y(5)*sin(thetai1))^2)^0.5; deltasum1=A10-A01;%公式14 if deltasum1>=0 H=1; else H=0; end Q1=Kc1*(deltasum1^1.5)*H; FBX1(1,i)=Q1*cos(thetai1)*cos(alpha10); FBY1(1,i)=Q1*cos(thetai1)*sin(alpha10); FBZ1(1,i)=Q1*sin(alpha10); end Fbx1=sum(FBX1); Fby1=sum(FBY1); Fbz1=sum(FBZ1); for i=1:1:Nb1 thetai2=ommcage1*t+2*pi*(i-1)/Nb1;%公式23 alpha20=atan((A01*sin(alpha01)+y(21)+Rj1*y(15)*cos(thetai2))/(A01*cos(alpha01)+y(17)*cos(thetai2)+y(19)*sin(thetai2)));%公式19 A20=((A01*sin(alpha01)+y(21)+Rj1*y(15)*cos(thetai2))^2+(A01*cos(alpha01)+y(17)*cos(thetai2)+y(19)*sin(thetai2))^2)^0.5; deltasum2=A20-A01;%公式14 if deltasum2>=0 H=1; else H=0; end Q2=Kc1*(deltasum2^1.5)*H;%公式24 FBX2(1,i)=Q2*cos(thetai2)*cos(alpha20); FBY2(1,i)=Q2*cos(thetai2)*sin(alpha20); FBZ2(1,i)=Q2*sin(alpha20); end Fbx2=sum(FBX2); Fby2=sum(FBY2); Fbz2=sum(FBZ2); for i=1:1:Nb2 thetai3=ommcage2*t+2*pi*(i-1)/Nb2; alpha30=atan((A01*sin(alpha01)+y(27)+Rj1*y(15)*cos(thetai2))/(A01*cos(alpha01)+y(23)*cos(thetai2)+y(25)*sin(thetai2)));%公式19 A30=((A01*sin(alpha01)+y(27)+Rj1*y(15)*cos(thetai2))^2+(A01*cos(alpha01)+y(23)*cos(thetai2)+y(25)*sin(thetai2))^2)^0.5; deltasum3=A30-A02;%公式14 if deltasum3>=0 H=1; else H=0; end Q3=Kc2*(deltasum3^1.5)*H;%公式24 FBX3(1,i)=Q3*cos(thetai3)*cos(alpha30); FBY3(1,i)=Q3*cos(thetai3)*sin(alpha30); FBZ3(1,i)=Q3*sin(alpha30); end Fbx3=sum(FBX3); Fby3=sum(FBY3); Fbz3=sum(FBZ3); for i=1:1:Nb2 thetai4=ommcage2*t+2*pi*(i-1)/Nb2;%公式23 alpha40=atan((A01*sin(alpha01)+y(41)+Rj1*y(15)*cos(thetai2))/(A01*cos(alpha01)+y(37)*cos(thetai2)+y(39)*sin(thetai2)));%公式19 A40=((A01*sin(alpha01)+y(41)+Rj1*y(15)*cos(thetai2))^2+(A01*cos(alpha01)+y(37)*cos(thetai2)+y(39)*sin(thetai2))^2)^0.5; deltasum4=A40-A02;%公式14 if deltasum4>=0 H=1; else H=0; end Q4=Kc2*(deltasum4^1.5)*H;%公式24 FBX4(1,i)=Q4*cos(thetai4)*cos(alpha40); FBY4(1,i)=Q4*cos(thetai4)*sin(alpha40); FBZ4(1,i)=Q4*sin(alpha40); end Fbx4=sum(FBX4); Fby4=sum(FBY4); Fbz4=sum(FBZ4); %ω=omm φ=phi α=alp dy=[y(2); (Td-ct1*(y(2)-y(44))-kt1*(y(1)-y(43)))/Jd; y(4); (Fbx1-csx1*(y(4)-y(10))-cbx1*y(4)-ksx1*(y(3)-y(9)))/mb1; y(6); (Fby1-mb1*g-csy1*(y(6)-y(12))-cby1*y(6)-ksy1*(y(5)-y(11)))/mb1; y(8); (Fbz1-csz1*(y(8)-y(14))-cbz1*y(8)-ksz1*(y(7)-y(13)))/mb1; y(10); (-Fx+m1*p1*dy(44)*sin(phi1)+m1*p1*cos(phi1)*((omm1+y(44))^2)-csx1*(y(10)-y(4))-csx2*(y(10)-y(18))... -ksx1*(y(9)-y(3))-ksx2*(y(9)-y(17)))/m1; y(12); (-m1*g-Fy-m1*p1*cos(phi1)*dy(44)+m1*p1*sin(phi1)*((omm1+y(44))^2)-csy1*(y(12)-y(6))... -csy2*(y(12)-y(20))-ksy1*(y(11)-y(5))-ksy2*(y(11)-y(19)))/m1; y(14); (-Fz-csz1*(y(14)-y(8))-csz2*(y(14)-y(20))-ksz1*(y(13)-y(7))-ksz2*(y(13)-y(21)))/m1; y(16); (m1*p1*dy(10)*sin(phi1)-m1*p1*dy(12)*cos(phi1)-F*rb1-ct1*(y(44)-y(2))-kt1*(y(43)-y(1)))/(J1+m1*p1^2); y(18); (Fbx2-csx2*(y(18)-y(10))-ksx2*(y(17)-y(9))-cbx2*y(18))/mb2; y(20); (Fby2-mb2*g-csy2*(y(20)-y(12))-