24小时热门版块排行榜

返回列表

欧阳球飞

铜虫 (小有名气)

应助: 0 (幼儿园)
金币: 12.6
散金: 25
帖子: 62
在线: 42.6小时
虫号: 3799827
注册: 2015-04-10
专业: 机械工程

[求助] 齿轮－转子－滚动轴承传动系统的弯扭耦合振动已有3人参与

使用matlab编程，利用龙格库塔求解齿轮－转子－滚动轴承传动系统的弯扭耦合振动模型，程序出现问题，没找出原因求指点。
图片6是程序出现的问题：  试图访问 dtx(8)；由于 numel(dtx)=3，索引超出范围。
主程序：
function dtx = zhichilun(t,x)

%齿轮参数
m1=667.859; m2=141.038;
J1=44.355;  J2=0.207;
Td=100; Tl=100;
p1=2e-6; p2=2e-5;%偏心距
z1=100; z2=25;
R1 = (z1*8*cos(pi/9))*(1/2)/1000 ;%基圆半径
R2 = (z2*8*cos(pi/9))*(1/2)/1000;
me=J1*J2/(J1*R2*R2+J2*R1*R1);%等效质量
km=6e8;
cm=2*0.1*sqrt(km*me);
e0=2e-5; er=3e-5;
we=2000;
w1=2*pi*500*100/60;
w2=2*pi*2000*25/60;
et=e0+er*sin(we*t);%齿轮啮合误差
det=er*we*cos(we*t);
F=cm*(x(6)-x(12)+R1*x(8)-R2*x(14)+p1*w1*cos(w1*t+x(7))-p2*w2*cos(w2*t+x(13))-det)+...
  km*(x(5)-x(11)+R1*x(7)-R2*x(13)+p1*sin(w1*t+x(7))-p2*cos(w2*t+x(13))-et);
%轴承
mb1=115.2; mb2=115.2;
mb3=73.2; mb4=73.2;
kbx1=6e10; kby1=6e10;
kbx2=6e10; kby2=6e10;
cbx1=2*0.01*sqrt(kbx1*mb1*mb2/(mb1+mb2)); cby1=2*0.01*sqrt(kby1*mb1*mb2/(mb1+mb2));
cbx2=2*0.01*sqrt(kbx2*mb3*mb4/(mb3+mb4)); cby2=2*0.01*sqrt(kby2*mb3*mb4/(mb3+mb4));

%轴
ks1=2.513e8; ks2=0.157e8;%弯曲刚度及阻尼
cs1=2*0.07*sqrt(ks1*m1*mb1/(m1+mb1));
cs2=2*0.07*sqrt(ks2*m2*mb3/(m1+mb3));
Jd=17.883; Jl=5.271;
kt1=4.027e8; kt2=0.254e8;%扭转刚度及阻尼
ct1=2*0.07*sqrt(kt1*Jd*J1/(Jd+J1));
ct2=2*0.07*sqrt(kt2*Jl*J2/(Jl+J2));
s1=0.5; s2=0.5; s3=0.5; s4=0.5;%齿轮处于传动轴的位置

Fx1=0;Fy1=0; Fx2=0;Fy2=0; Fx3=0;Fy3=0; Fx4=0;Fy4=0;

Nb1=14;

R1=0.25; r1=0.2;
r01=0.25e-6;
kb1=13.34e9;
%左1轴承
for i=1:Nb1
theta(i)=(w1*r1)/(R1+r1)*t+2*pi/Nb1*(i-1);
distance(i)=x(15)*cos(theta(i))+x(17)*sin(theta(i))-r01;%%
if distance(i)<=0;
   distance(i)=0;
end
Fx1=Fx1+kb1*distance(i)^1.5*cos(theta(i));
Fy1=Fy1+kb1*distance(i)^1.5*sin(theta(i));
end

%右1轴承
for j=1:Nb1
theta(j)=(w1*r1)/(R1+r1)*t+2*pi/Nb1*(j-1);
distance(j)=x(19)*cos(theta(j))+x(21)*sin(theta(j))-r01;%%
if distance(j)<=0;
   distance(j)=0;
end
Fx2=Fx2+kb1*distance(j)^1.5*cos(theta(j));
Fy2=Fy2+kb1*distance(j)^1.5*sin(theta(j));
end

Nb2=10;

R2=0.125; r2=0.1;
r02=0.5e-6;
kb2=10.34e9;
%左2轴承
for u=1:Nb2
theta(u)=(w2*r2)/(R2+r2)*t+2*pi/Nb2*(u-1);
distance(u)=x(23)*cos(theta(u))+x(25)*sin(theta(u))-r02;%%
if distance(u)<=0;
   distance(u)=0;
end
Fx3=Fx3+kb2*distance(u)^1.5*cos(theta(u));
Fy3=Fy3+kb2*distance(u)^1.5*sin(theta(u));
end

%右2轴承
for v=1:Nb2
theta(v)=(w2*r2)/(R2+r2)*t+2*pi/Nb2*(v-1);
distance(v)=x(23)*cos(theta(v))+x(25)*sin(theta(v))-r02;%%
if distance(v)<=0;
   distance(v)=0;
end
Fx4=Fx4+kb2*distance(v)^1.5*cos(theta(v));
Fy4=Fy4+kb2*distance(v)^1.5*sin(theta(v));
end

%%齿轮系统模型
%输入端扭转角
dtx(1)=x(2);
dtx(2)=(Td-ct1*(x(2)-x(8))-kt1*(x(1)-x(7)))/Jd;
%主动齿轮沿 x,y的振动位移，扭转角
dtx(3)=x(4);
dtx(4)= (m1*p1*dtx(8)*sin(w1*t+x(7))+m1*p1*(w1+x(8))^2*cos(w1*t+x(7))-cs1*(x(4)-s2*x(16)-s1*x(20))-ks1*(x(3)-s2*x(15)-s1*x(19)))/m1;
dtx(5)=x(6);
dtx(6)= (m1*p1*dtx(8)*cos(w1*t+x(7))-m1*p1*(w1+x(8))^2*sin(w1*t+x(7))-F-m1*g-cs1*(x(6)-s2*x(18)-s1*x(22))...
   -ks1*(x(5)-s2*x(17)-s1*x(21)))/m1;
dtx(7)=x(8);
dtx(8)= (m1*p1*sin(w1*t+x(7))*dtx(4)+m1*p1*cos(w1*t+x(7))*dtx(6)-F*rb1-ct1*(x(8)-x(2))-kt1*(x(7)-x(1)))/(J1+m1*p1*p1);
%从动齿轮沿 x,y的振动位移，扭转角
dtx(9)=x(10);
dtx(10)= (m2*p2*dtx(14)*sin(w2*t+x(13))+m2*p2*(w2+x(14))^2*cos(w2*t+x(13))-cs2*(x(10)-s4*x(24)-s3*x(28))...
   -ks2*(x(9)-s4*x(23)-s3*x(27)))/m2;
dtx(11)=x(12);
dtx(12)= (-m2*p2*dtx(14)*cos(w2*t+x(13))+m2*p2*(w2+x(14))^2*sin(w2*t+x(13))+F-m2*g-cs2*(x(12)-s4*x(26)-s3*x(30))...
   -ks2*(x(11)-s4*x(25)-s3*x(29)))/m2;
dtx(13)=x(14);
dtx(14)= (m2*p2*sin(w2*t+x(13))*dtx(10)-m2*p2*cos(w2*t+x(13))*dtx(12)+F*rb2-ct2*(x(14)-x(32))-kt2*(x(13)-x(31)))/(J2+m2*p2*p2);
%轴承1振动平移位移
dtx(15)=x(16);
dtx(16)= (Fx1-cs1*s2*(-x(4)+s2*x(16)+s1*x(20))-cbx1*x(16)-ks1*s2*(-x(3)+s2*x(15)+s1*x(19)))/mb1;
dtx(17)=x(18);
dtx(18)= (Fy1-mb1*g-cs1*s2*(-x(6)+s2*x(18)+s1*x(22))-cby1*x(18)-ks1*s2*(-x(5)+s2*x(17)+s1*x(21)))/mb1;
%轴承2振动平移位移
dtx(19)=x(20);
dtx(20)= (Fx2-cs1*s1*(-x(4)+s2*x(16)+s1*x(20))-cbx2*x(20)-ks1*s1*(-x(3)+s2*x(15)+s1*x(19)))/mb2;
dtx(21)=x(22);
dtx(22)= (Fy2-mb2*g-cs1*s1*(-x(6)+s2*x(18)+s1*x(22))-cby2*x(22)-ks1*s1*(-x(5)+s2*x(17)+s1*x(21)))/mb2;
%轴承3振动平移位移
dtx(23)=x(24);
dtx(24)= (Fx3-cs2*s4*(-x(10)+s4*x(24)+s3*x(28))-cbx3*x(24)-ks2*s4*(-x(9)+s4*x(23)+s3*x(27)))/mb3;
dtx(25)=x(26);
dtx(26)= (Fy3-mb3*g-cs2*s4*(-x(12)+s4*x(26)+s3*x(30))-cby3*x(26)-ks2*s4*(-x(11)+s4*x(25)+s3*x(29)))/mb3;
%轴承4振动平移位移
dtx(27)=x(28);
dtx(28)= (Fx4-cs2*s3*(-x(10)+s4*x(24)+s3*x(28))-cbx4*x(28)-ks2*s3*(-x(9)+s4*x(23)+s3*x(27)))/mb4;
dtx(29)=x(30);
dtx(30)= (Fy4-mb4*g-cs2*s3*(-x(12)+s4*x(26)+s3*x(30))-cby4*x(30)-ks2*s3*(-x(11)+s4*x(25)+s3*x(29)))/mb4;
%输出端扭转
dtx(31)=x(32);
dtx(32)= (-Tl-ct2*(x(32)-x(14))-kt2*(x(31)-x(13)))/Jl;
dtx=dtx';

求解程序：
clear all;
close all;
clc;
x0 = rand(1,32)/1e30;tic,
[t,y]=ode45(@zhichilun,[0,2],x0);

figure
plot(t,y(:,1))
xlabel('时间/s'),ylabel('扭转角/rad'),title('输入端扭转角')
figure
plot(t,y(:,3))
xlabel('时间/s'),ylabel('幅值/m'),title('主动轮x向响应')

1.jpg