传递矩阵matlab程序.docx

1、传递矩阵matlab程序%main_critical.m%该程序使用Riccati传递距阵法计算转子系统的临界转速及振型%本函数中均采用国际单位制% 第一步：设置初始条件（调用函数shaft_parameters）%初始值设置包括：轴段数N，搜索次数M%输入轴段参数：内径d,外径D，轴段长度l，支撑刚度K，单元质量mm，极转动惯量JppN,M,d,D,l,K,mm,Jpp=shaft_parameters;% 第二步：计算单元的5个特征值（调用函数shaft_pra_cal）%单元的5个特征值：%m_k：:质量%Jp_k：极转动惯量%Jd_k：直径转动惯量%EI：弹性模量与截面对中性轴的惯性矩

2、的乘积%rr：剪切影响系数m_k,Jp_k,EI,rr=shaft_pra_cal(N,D,d,l,Jpp,mm);% 第三步：计算剩余量（调用函数surplus_calculate）,并绘制剩余量图%剩余量：D1for i=1:1:M ptx(i)=0; pty(i)=0;endfor ii=1:1:M wi=ii/1*2+50; D1,SS,Sn=surplus_calculate(N,wi,K,m_k,Jp_k,JD_k,l,EI,rr); D1;pty(ii)=D1;ptx(ii)=w1endylabel(剩余量)；plot(ptx,pty)xlabel(角速度red/s);grid

3、on% 第四步：用二分法求固有频率及振型图%固有频率：Critical_speedwi=50;for i=1:1:4order=i D1,SS,Sn=surplus_calculate(N,wi,k,m_k,Jp_k,Jd_k,l,EI,rr); Step=1;D2=D1;kkk=1;while kkk0 wi=wi+step; D2=D1; D1,SS,Sn=surplus_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr); end if D1*D20 wi=wi-step; step=step/2; wi=wi+step; D1,SS,Sn =surplus_

4、calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr); EndD1;Wi;If atepFen*60 for i=1:1:N+1 x(i,1)=yn(i*4-3,n)*le6; y(i,1)=yn(i*4-2,n)*le6; sitax(I,1)=yn(i+4-0,n)/pi*180 end u=F(loca*4-4+1,1); ut1=ut1;t u; xt1=xt1;t x; yt1=yt1;t y; endendrub_signsave xt1.dat xt1 ascii;save yt1.dat yt1 ascii;save ut1.dat ut1 ascii

5、;%initial_conditions.m%该程序主要进行仿真条件设置Function rub_sign loca loc_rob Famp wi r Famp1 fai=initial_conditions%需要设置的初始条件有：%rub_sign: 碰摩标志，若rub_sign=0,说明系统无碰摩故障；否则 rub_sign=1%loca: 不平衡质量的位置%loc_rub: 碰摩位置%Famp: 不平衡质量的大小 单位为：g%wi: 转速 单位为：rad% r： 偏心半径 单位为:mm%Famp1： 离心力的大小 单位为:kg,m%fai: 不平衡量的初始相位 radrub_sign=

6、0;loca=6;loc_rub=8;Famp=1;wi=3000/60*2*pi;r=30Famp1=Famp(1)/1000*wi2*r/1000;fai=30 30/180*pi %roto_parameters.m%该程序对Jeffcott转子系统进行参数设置 function N density Ef L R R0 miu=rotor_parameters%N: 划分的轴段数%density: 轴的密度 单位为：kg/m3%Ef: 轴的弹性模量 单位为:Pa%L 每个轴段的长度 单位为:m%R 每个轴段的外半径 单位为:m%R0: 每个轴段的内半径 单位为:m%miu: 每个轴段的单

7、元质量 kg/mN=11;Density=7850;Ef=2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1*lell;L=90.5 90.5 80.5 62.5 30 20 22.5 62.5 90.5 90.5 90.5/1000;R=20 20 20 20 20 90 20 20 20 20 20/2000;R0=0 0 0 0 0 0 0 0 0 0 0/2000;for i=1:1;N miu(i)=density*pi*(R(i)2-R0(i)2)end%Mst_Msr_Ks_Ge.m%该程序设置单元距阵Function Mst Msr Ks G

8、e= Mst_Msr_Ks_Ge(N,density,R,R0,L,Ef,miu)%Mst: 移动单元质量距阵%Msr: 转动单元质量距阵%Ks: 刚度单元距阵%Ge: 陀螺力距单元距阵NN0=N;NN1=NN0+1NN2=NN1+1for i=1:1:NN0 Mst(1,1,i)=156;Mst(2,1,i)=0; Mst(2,2,i)=156;Mst(3,1,i)=0; Mst(3,2,i)=-22*L(i); Mst(3,3,i)= 4*L(i)2; Mst(4,1,i)22*L(i); Mst(4,2,i)=0; Mst(4,3,i)=0;Mst(4,4,i)=4*L(i)2; Mst

9、(5,1,i)=54; Mst(5,2,i)=0Mst(5,3,i)=0; Mst(5,4,i)=13*L(i); Mst(6,1,i)=0;Mst(6,2,i)=54; Mst(6,3,i)=13*L(i); Mst(6,4,i)=0;Mst(7,1,i)=0; Mst(7,2,i)=13*L(i); Mst(7,3,i)=-3*L(i)2;Mst(7.4.i)=0; Mst(8,1,i)=-13*L(i); Mst(8,2,i)=0;Mst(8,3,i)=0; Mst(8,4,i)=-3*L(i)2; Mst(5,5,i)=156;Mst(6,5.i)=0; Mst(6,6,i)=156;

10、 Mst(7,5,i)=0; Mst(7,6,i)=22*L(i); Mst(7,7,i)=4*L(i)2;Mst(8,5,i)=-22*L(i); Mst(8,6,i)=0 Mst(8,7,i)=0Mst(8,8,i)=4*L(i)2;endfor i=1:1:NN0Msr(1,1,i)=36;Msr(2,1,i)=0; Msr(2,2,i)=36;Msr(3,1,i)=0 Msr(3,2,i)=-3*L(i); Msr(3,3,i)=4*L(i)2;Msr(4,1,i)=3*L(i); Msr(4,2,i)=0; Msr(4,3,i)=0;Msr(4,4,i)=4*L(i)2; Msr(5

11、,1,i)=36; Msr(5,2,i)=0;Msr(5,3,i)=0; Msr(5,4,i)=-3*L(i); Msr(6,1,i)=0;Msr(6,2,i)=-36; Msr(6,3,i)= 3*L(i); Msr(6,4,i)=0;Msr(7,1,i)=0; Msr(7,2,i)=-3*L(i); Msr(7,3,i)=-L(i)2;Msr(7,4,i)=0; Msr(8,1,i)=3*L(i); Msr(8,2,i)=0;Msr(8,3,i)=0; Msr(8,4,i)=-L(i)2; Msr(5,5,i)=36;Msr(6,5,i)=0; Msr(6,6,i)=36; Msr(7,5

12、,i)=0;Msr(7,6,i)=3*L(i); Msr(7,7,i)=4*L(i)2; Msr(8,5,i)=-3*L(i);Msr(8,6,i)=0; Msr(8,7,i)=0; Msr(8,8,i)=4*L(i)2;endfor i=1:1:NN0Ge(1,1,i)=0;Ge(2,1,i)=36; Ge(2,2,i)=0;Ge(3,1,i)=-3*l(i); Ge(3,2,i)=0; Ge(3,3,i)=0;Ge(4,1,i)=0; Ge(4,2,i)=-3*L(i); Ge(4,3,i)=4*L(i)2;Ge(4,4,i)=0; Ge(5,1,i)=0; Ge(5,2,i)=36;Ge

13、(5,3,i)=-3*L(i); Ge(5,4,i)=0; Ge(6,1,i)=-36;Ge(6,2,i)=0; Ge(6,3,i)=0; Ge(6,4,i)=-3*L(i0);Ge(7,1,i)=-3*L(i); Ge(7,2,i)=0; Ge(7,3,i)=0;Ge(7,4,i)=L(i)62; Ge(8,1,i)=0; Ge(8,2,i)=-3*L(i);Ge(8,3,i)=-L(i)2; Ge(8,4,i)=0; Ge(5,5,i)=0;Ge(6,5,i)=36; Ge(6,6,i)=0; Ge(7,5,i)=3*L(i);Ge(7,6,i)=0; Ge(7,7,i)=0; Ge(8,

14、5,i)=0;Ge(8,6,i)=3*L(i); Ge(8,7,i)=4*L(i)2; Ge(8,8,i)=0;endfor i=1:1:NN0Ks(1,1,i)=12;Ks(2,1,i)=0; Ks(2,2,i)=12;Ks(3,1,i)=0; Ks(3,2,i)=-6*L(i); Ks(3,3,i)=4*L(i)2;Ks(4,1,i)=6*L(i); Ks(4,2,i)=0; Ks(4,3,i)=0;Ks(4,4,i)=4*L(i)2; Ks(5,1,i)=-12; Ks(5,2,i)=0;Ks(5,3,i)=0; Ks(5,4,i)=-6*L(i); Ks(6,1,i)=0;Ks(6,2

15、,i)=-12; Ks(6,3,i)=6*L(i); Ks(6.4.i)=0;Ks(7,1,i)=0; Ks(7,2,i)=-6*L(i); Ks(7,3,i)=2*L(i)2;Ks(7,4,i)=0; Ks(8.1,i)=6*L(i); Ks(8,2,i)=0;Ks(8,3,i)=0; Ks(8,4,i)=2*L(i)62; Ks(5,5,i)=12;Ks(6,5,i)=0; Ks(6,6,i)=12; Ks(7,5,i)=0;Ks(7,6,i)=6*L(i); Ks(7,7,i)=4*L(i)2; Ks(8,5,i)=-6*L(i);Ks(8,6,i)=0; Ks(8,7,i)=0; Ks

16、(8,8,i)=4*L(i)2;endfor i=1:1:NN0 for j=1:1:8 for k=1:1:8 EI=Ef(i)*pi*(R(i)4-R0(i)4-)/4; Mst(j,k,i)=Mst(j,k,i)*miu(i)*L(i)/420; Msr(j,k,i)=Msr(j,k,i)*miu(i)*R(i)2/120/L(i); Ge(j,k,i)=-Ger(j,k,i)*2*miu(i)*R(i)2/120/L(i); Ks(j,k,i)=Ks(j,k,i)*EI/L(i)3; end endendfor i=1:1:NN0 for j=1:1:8 for k=j:1:8 Mst(j,k,i)=Mst(k,j,i); Msr(j,k,i)=Msr(k,j,i); Ks(j,k,i)=Ks(k,j,i); G