2.txt

齿轮动力学分析

clc;
N=5;
ha=1;ns=50;
z1=40;z2=33;z3=106;m=0.003;a=22*pi/180; pb=pi*m*cos(a);
r1=m*z1/2;r2=m*z1/2;r3=m*z3/2;
ra1=m*z1/2+ha*m;ra2=m*z2/2+ha*m;ra3=m*z3/2-ha*m;
rb1=m*z1*cos(a)/2; rb2=m*z2*cos(a)/2; rb3=m*z3*cos(a)/2;
%质量
ms=4.75;m1=2.88;m2=2.88;m3=2.88;mr=16.92;
%j1=m1*rb1^2;j2=m2*rb2^2;j3=m3*rb3^2;
M=[ms,0,0,0,0;
   0,m1,0,0,0;
   0,0,m2,0,0;
   0,0,0,m3,0;
   0,0,0,0,mr];
%%%刚度
l1=sqrt(ra1^2-rb1^2)+sqrt(ra2^2-rb2^2)-(m*z1/2+m*z2/2)*sin(a);
l2=sqrt(ra2^2-rb2^2)-(sqrt(ra3^2-rb3^2)-(m*z3/2-m*z2/2)*sin(a));
e1=l1/pb;e2=l2/pb;
w=ns*z1;
T=1/w;
Ks1=1.0*(1e+9);Ks2=6.0*(1e+8);Asp1=(Ks1+Ks2)/2;Asp2=(Ks1-Ks2)/2;
Kr1=1.2*(1e+9);Kr2=8.0*(1e+8);Arp1=(Kr1+Kr2)/2;Arp2=(Kr1-Kr2)/2;
AKs=Ks1*(2-e1)+Ks2*(e1-1);
AKr=Kr1*(2-e2)+Kr2*(e2-1);
T1T2=sqrt(ra1^2-rb1^2)+(pi-a-a)*rb1+sqrt(ra3^2-rb3^2)-(r3-r2)*sin(a);
sb=pb/2;
c=mod((T1T2-sb),pb); 
%?????1/d??
d=64/T;
totaltime=100*T;
tt=totaltime*d+2;
%totaltime?????????tt????????????????d???
u=zeros(N,tt);
du=zeros(N,tt);
ddu=zeros(N,tt);
fs=zeros(3,tt);
fr=zeros(3,tt);
Gs=zeros(1,tt);
Gr=zeros(1,tt);
ks=zeros(3,tt);
kr=zeros(3,tt);
% u(:,1)=[0;0;0;0;0];
u(:,1)=[1e-7*0.8229;1e-7*(-0.4268);1e-7*(-0.4268);1e-7*(-0.4268);1e-7*(-0.0900)];
du(:,1)=[0;0;0;0;0];
ddu(:,1)=[0;0;0;0;0];
%nuwmark参数
b=0.25;   r=0.5;  dt=1/d;
a0=1/(b*(dt)^2); a1=r/(b*dt);
a2=1/(b*dt);     a3=1/(2*b)-1;
a4=r/b-1;        a5=0.5*dt*((r/b)-2);
a6=dt*(1-r);     a7=r*dt;
f=zeros(N,tt);
Ef=zeros(N,tt);
t1=1;
%initial force ???
FD=90;
FL=0;
for t=0:1/d:totaltime
    f(:,t1)=[FD;0;0;0;FL];
    t1=t1+1;
end
t1=1;
for t=0:1/d:totaltime 
ks(1,t1)=Asp1+Asp2*square(2*pi*w*t,(e1-1)*100);
ks(2,t1)=Asp1+Asp2*square(2*pi*w*(t-z1/3*T),(e1-1)*100);
ks(3,t1)=Asp1+Asp2*square(2*pi*w*(t-z1*2/3*T),(e1-1)*100);
kr(1,t1)=Arp1+Arp2*square(2*pi*w*(t-c*T),(e2-1)*100);
kr(2,t1)=Arp1+Arp2*square(2*pi*w*(t+z3/3*T-c*T),(e2-1)*100);
kr(3,t1)=Arp1+Arp2*square(2*pi*w*(t+z3*2/3*T-c*T),(e2-1)*100);
ksu=0;  kru=10000000000;
K=[ks(1,t1)+ks(2,t1)+ks(3,t1),ks(1,t1),ks(2,t1),ks(3,t1),0;
   ks(1,t1),(ks(1,t1)+kr(1,t1)),0,0,-kr(1,t1);
   ks(2,t1),0,(ks(1,t1)+kr(2,t1)),0,-kr(2,t1);
   ks(3,t1),0,0,(ks(1,t1)+kr(3,t1)),-kr(3,t1);
   0,-kr(1,t1),-kr(2,t1),-kr(3,t1),kr(1,t1)+kr(2,t1)+kr(3,t1)+kru];

%EK?????????
%  EK=K+a0*M;  
%  Ef(:,t1+1)=f(:,t1+1)+M*(a0*u(:,t1)+a2*du(:,t1)+a3*ddu(:,t1));
    C=1000*M;
    EK=K+a0*M+a1*C;
    Ef(:,t1+1)=f(:,t1+1)+M*(a0*u(:,t1)+a2*du(:,t1)+a3*ddu(:,t1))+C*(a1*u(:,t1)+a4*du(:,t1)+a5*ddu(:,t1));
    u(:,t1+1)= inv(EK)*Ef(:,t1+1);
    ddu(:,t1+1)=a0*(u(:,t1+1)-u(:,t1))-a2*du(:,t1)-a3*ddu(:,t1);
    du(:,t1+1)=du(:,t1)+a6*ddu(:,t1)+a7*ddu(:,t1+1);
    fs(1,t1)=ks(1,t1)*(u(1,t1)+u(2,t1)); fs(2,t1)=ks(2,t1)*(u(1,t1)+u(3,t1));
    fs(3,t1)=ks(3,t1)*(u(1,t1)+u(4,t1));
    fr(1,t1)=kr(1,t1)*(u(5,t1)-u(2,t1));fr(2:t1)=kr(2,t1)*(u(5,t1)-u(3,t1));
    fr(3,t1)=kr(3,t1)*(u(5,t1)-u(4,t1));
    Gs(1,t1)=fs(1,t1)/FD; Gs(2,t1)=fs(2,t1)/FD;  Gs(3,t1)=fs(3,t1)/FD;
  % Gr(1,t1)=fr(1,t1)/FD;
%  Gs(1,t1)=3*rb1*max(max(fs(1,t1),fs(2,t1)),fs(3,t1))/FD;
%  Gr(1,t1)=3*rb3*max(max(fr(1,t1),fr(2,t1)),fr(3,t1))/FL;
   t1=t1+1;
end
t=0:1/d:totaltime+1/d;

% subplot(2,2,1),plot(t/T,u(1,:)),xlabel('t/T'),ylabel('u');
% subplot(2,2,2),plot(t/T,du(1,:)),xlabel('t/T'),ylabel('du');
% subplot(2,2,3),plot(t/T,ks(1,:)),xlabel('t/T'),ylabel('ks');
% subplot(2,2,4),plot(t/T,Gs(1,:),t/T,0.3),xlabel('t/T'),ylabel('R1');
% 
% figure(1)
% plot(t/T,Gs(1,:),t/T,0.3),xlabel('t/T'),ylabel('u');
% figure(2)
% plot(t/T,Gs(2,:),t/T,0.3),xlabel('t/T'),ylabel('R');
% figure(3)
% plot(t/T,Gs(3,:),t/T,0.3),xlabel('t/T'),ylabel('ks');

%plot(t/T,Gs),xlabel('t/T'),ylabel('R1'); hold on
%plot(t/T,fs1); hold on
% plot(t/T,u(1,:)); hold on
%plot(t/T,ks(1,:),t/T,ks(2,:));

%进行fft变换
% n=totaltime*d+1;fs=d;
% y=fft(Gs(1,:),n);%进行fft变换
% mag=abs(y)*2/n;%求幅值 乘上后面的2/N得到正确幅值
% f=([0:length(y)]-1)*fs/length(y);%进行对应的频率转换
% plot(f(1:n/2),mag(1:n/2));%做频谱图
% xlabel('frequency(Hz)');
% ylabel('Amplitude');
% title('spectrum map');
grid on
clear