sjjshguiaghw_齿轮_动力学_行星齿轮_行星齿轮箱

%% newmark method;linear systems;five degree freedom，纽马克方法；线性系统；八个自由度 %% A、初始计算 %振动常微分非齐次方程: M*x''+C*x'+K*x=F clc clear all close x0=[0 0 0 0 0 0]; %初位移 v0=[0 0 0 0 0 0]; %初速度 x(:,1)=x0; v(:,1)=v0; dt=1/2000; %时间步长 time=20; beta=0.25; %参数设置 beta>=1/4*(delta+1/2)^2 delta=0.5; %delta>=0.5 %1、初始计算 mp=157;ms=538;mr=1422;jc=146/(0.326)^2;jr=1422/(0.538)^2; jp=5/(0.213)^2; js=538/(0.114)^2;zr=104; aa=pi*5/36;%压力角 M=diag([jc+3*mp jr js jp jp jp]); kct=0;kst=0;krt=1*10^9; Ts=1579;Tc=2400; %% b0=1/(beta*dt^2); %积分常数 b1=delta/(beta*dt); b2=1/(beta*dt); b3=1/(2*beta)-1; b4=delta/beta-1; b5=dt/2*(delta/beta-2); b6=dt*(1-delta); b7=dt*delta; ksm=5*10^8;krm=ksm*10; cs=2*0.07*sqrt(ksm*ms*mp/(ms+mp)); cr=2*0.07*sqrt(krm*mr*mp/(mr+mp)); %% B、每一时间增量计算 k=1;ratio=1; for tt=0:dt:time-dt w=15/60; ws=w*5.7; T=1./(w*zr); kb=diag([kct krt kst 0 0 0]); %[ks11, ks22, ks33]= meshstiffness(tt,T); ks11=(5+5/pi*sin(2*pi*(2-tt/T)))*10^8;%-5/pi*sin(2*pi*(0.5-tt/T)))*10^8; %dks1=(5+5/pi*sin(2*pi*(2-tt*ws))-5/pi*sin(2*pi*(0.5-tt*ws)))*10^8; ks1=ks11%+(1-ratio)*dks1; kr1=(5+5/pi*sin(2*pi*(1-tt/T)))*10^9;%-5/pi*sin(2*pi*(0.5-tt/T)))*10^9; ks22=(5+5/pi*sin(2*pi*(2-tt/T)))*10^8%-5/pi*sin(2*pi*(0.5-tt/T)))*10^8; %dks2=(5+5/pi*sin(2*pi*(2-tt*ws))-5/pi*sin(2*pi*(0.5-tt*ws)))*10^8; ks2=ks22%+(1-ratio)*dks2; kr2=(5+5/pi*sin(2*pi*(-1-tt/T)))*10^9%-5/pi*sin(2*pi*(-1.5-tt/T)))*10^9; ks33=(5+5/pi*sin(2*pi*(3-tt/T)))*10^8%-5/pi*sin(2*pi*(2.5-tt/T)))*10^8; %dks3=(5+5/pi*sin(2*pi*(3-tt*ws))-5/pi*sin(2*pi*(2.5-tt*ws)))*10^8; ks3=ks33%+(1-ratio)*dks3; kr3=(5+5/pi*sin(2*pi*(-3-tt/T)))*10^9%-5/pi*sin(2*pi*(-3.5-tt/T)))*10^9; % k11=(ks1+ks2+ks3)*(cos(aa))^2+(kr1+kr2+kr3)*(cos(aa))^2; k12=-(kr1+kr2+kr3)*cos(aa);k13=-(ks1+ks2+ks3)*cos(aa);k14=kr1*cos(aa)-ks1*cos(aa); k15=kr2*cos(aa)-ks2*cos(aa);k16=kr3*cos(aa)-ks3*cos(aa); k22=(kr1+kr2+kr3);k23=0;k24=-kr1;k25=-kr2;k26=-kr3;k33=(ks1+ks2+ks3);k34=ks1;k35=ks2;k36=ks3; k44=kr1+ks2;k45=0;k46=0;k55=kr2+ks2;k56=0;k66=kr2+ks3; k21=k12;k31=k13;k32=k23;k41=k14;k42=k24;k43=k34;k51=k15;k52=k25;k53=k35;k54=k45; k61=k16;k62=k26;k63=k36;k64=k46;k65=k56; Kt=[k11 k12 k13 k14 k15 k16;k21 k22 k23 k24 k25 k26;k31 k32 k33 k34 k35 k36; ... k41 k42 k43 k44 k45 k46;k51 k52 k53 k54 k55 k56;k61 k62 k63 k64 k65 k66]; % c11=3*cs*(cos(aa))^2+3*cr*(cos(aa))^2;c12=-3*cr*cos(aa);c13=-3*cs*cos(aa);c14=cr*cos(aa)-cs*cos(aa); c15=cr*cos(aa)-cs*cos(aa);c16=cr*cos(aa)-cs*cos(aa); c22=3*cr;c23=0;c24=-cr;c25=-cr;c26=-cr;c33=3*cs;c34=cs;c35=cs;c36=cs; c44=cr+cs;c45=0;c46=0;c55=cr+cs;c56=0;c66=cr+cs; c21=c12;c31=c13;c32=c23;c41=c14;c42=c24;c43=c34;c51=c15;c52=c25;c53=c35;c54=c45; c61=c16;c62=c26;c63=c36;c64=c46;c65=c56; C=[c11 c12 c13 c14 c15 c16;c21 c22 c23 c24 c25 c26;c31 c32 c33 c34 c35 c36; ... c41 c42 c43 c44 c45 c46;c51 c52 c53 c54 c55 c56;c61 c62 c63 c64 c65 c66]; % K1=b0*M+Kt+b1*C; %等效刚度矩阵 f0(:,k)=[-Tc 0 0 0 0 0]; %外扰力 f0(:,k+1)=[-Tc 0 0 0 0 0]; a(:,1)=M\(f0(:,1)-Kt*x(:,1)-C*v(:,1)); f(:,k+1)=f0(:,k+1)+M*(b0*x(:,k)+b2*v(:,k)+b3*a(:,k))+C*(b1*x(:,k)+b4*v(:,k)+b5*a(:,k)); %等效载荷 x(:,k+1)=K1\f(:,k+1); a(:,k+1)=b0*(x(:,k+1)-x(:,k))-b2*v(:,k)-b3*a(:,k); v(:,k+1)=v(:,k)+b6*a(:,k)+b7*a(:,k+1); k=k+1; end t1=time/2;t2=time; t=0:dt:time; pli=0*t; figure; plot(t,a(3,:),t,pli,'r'); xlabel('时间(s)');ylabel('加速度(m/s2)'); title('太阳轮扭振加速度'); xlim([0 8]); %%傅立叶分析 fs=1/dt; ss=floor(1/dt); n=length(a(3,:)); N = 2^nextpow2(n); n1 = 0:N-1; f1=n1*fs/N; y1=fft(a(3,:),N); mag1=abs(y1)*2/N; figure; plot(f1(1:N/2),mag1(1:N/2),'r'); xlabel('频率(Hz)');ylabel('幅值(m/s2)'); title('太阳轮扭振加速度频谱');