利用传递矩阵法和Riccati传递矩阵法分析转子临界转速

利用传递矩阵法和Riccati传递矩阵法分析转子临界转速所需求解转子参数将转子简化为如下所示：三个盘的参数为：QUOTEm1=m2=另，阶梯轴的三段轴的截面惯性矩分别为：QUOTEJ1=1.7cm4J2三段轴的单位长度轴段的质量分别为：QUOTEm1=2.47kgmm试算转轴的传递矩阵取试算转速QUOTEp=ω=1000rads;则，各轴段的传递矩阵分别为：第1段QUOTEl=0.06mJ=1.7×10-8m4m=2.47第2段QUOTEl=0.15mJ=3.2×10-8m4第3段第4段第5段第6段此6段传递矩阵均采用MATLAB编程求解，MATLAB的源文件为HYPERLINKcalculate.m在MATLAB中使用线性插值法寻找最佳p值使得∆逼近于0。其程序为rotor.m经计算，考虑支撑刚度变化之间时，取时，一阶临界转速值为取时，一阶临界转速值为取时，一阶临界转速值为取时，一阶临界转速值为取时，一阶临界转速值为取时，一阶临界转速值为因此随着刚度的增加，一阶临界转速的值越来越大，而当不考虑支座的刚度变化，假设为完全刚性的话，一阶临界转速值为，因此当取时一阶临界转速值已相当接近完全刚性的情况。采用Riccati传递矩阵法进行各段轴的状态参数的传递根据Riccati传递矩阵法的原理，只需在传递矩阵法的基础之上求得各截面的Riccati传递矩阵。将转子截面的状态参数分组：,QUOTEfi={MQ}因为左端和右端均为自由端，故QUOTEf0=0，QUOTEe0≠0；QUOTEfN=0，QUOTEeN≠0;所以，我们可得到左端截面的Riccati传递矩阵QUOTES0=0；根据第i+1截面、QUOTEf和e之间的Riccati变换公式：可得（同样，试算转速选为QUOTEp=1000rads）：左盘右边截面的Riccati传递矩阵：QUOTES1=0250007e60第1轴段末的Riccati传递矩阵：QUOTES2=4.063e5-3.23e38.07e6-4.063e5刚性支承在此处的处理因为涉及到刚度，取QUOTE∞的情况，所以在获取其支承左边的Riccati传递矩阵后，需将QUOTEf2=S2e2转换为QUOTEe2=S2-1f2(即第二类Riccati变换)，然后再代入到普通传递矩阵的式子：;可得：。最终可得：,即QUOTE{e}3=00-3.0962e-40f3，即S3-1=00则，第2段轴段末的Riccati传递矩阵：可得：QUOTE[S]4=-8.309e51.270e5-5.134e6通过中间圆盘后，可得圆盘右边的Riccati传递矩阵：;则，第3轴段右边的Riccati传递矩阵为：;到达第二个刚性支承处，同样采用第二类Riccati变换，并带入Riccati传递矩阵公式：即可得：QUOTES7-1=00继续进行传递，第4轴段末的Riccati传递矩阵为：可得：QUOTE[S]8=-2.337e76.604e5-9.146e8由于，第4轴段末是球头联轴器，故，在此也要进行另外的推导，由于球头联轴器的力矩刚性系数趋近于0，则在此利用弹性铰链的传递矩阵：；；把QUOTEfi=SieiQUOTESi+1=Si(I+K则，可得：QUOTES9=00-8.765e7再根据第5轴段的传递矩阵，可得第5轴段末的Riccati传递矩阵：第5轴段末为一刚性支承，则同样采用第二种Riccati变换，可得刚性支承右端的Riccati传递矩阵：同样，根据第6轴段的传递矩阵，可得第6轴段末的Riccati传递矩阵：最后，通过最后一个圆盘的传递，可以得到转子最右端的Riccati传递矩阵：即，最右边截面上应该满足：QUOTEf13=S13e13；又，由于QUOTEf13={00}；所以，只有当QUOTEdetS13=0时，才能取到QUOTEe13的值，所以只有临界转速值p才能使QUOTEdetS13尽量接近于0。此处QUOTEdetS13=1.2e12；此处的算法MATLAB程序为calculate2.m利用线性插值法来求取最终的临界转速值p，只需将第一次求解的rotor.m的程序中的算法程序calculate.m改为calculate2.m即可。考虑支撑刚度变化之间时，取时，一阶临界转速值为取时，一阶临界转速值为取时，一阶临界转速值为取时，一阶临界转速值为取时，一阶临界转速值为取时，一阶临界转速值为因此随着刚度的增加，一阶临界转速的值越来越大，而当不考虑支座的刚度变化，假设为完全刚性的话，一阶临界转速值为，因此当取时一阶临界转速值已相当接近完全刚性的情况。附录：H.mfunctionoutput=H(in1,in2,in3,in4)%in1=lin2=Jin3=min4=psymsSTUVkaEJmpl;E=2e11;l=in1;J=in2;m=in3;p=in4;k=l*(m*(p^2)/(E*J))^0.25;S=0.5*(0.5*(exp(k)+exp(-k))+cos(k));T=0.5*(0.5*(exp(k)-exp(-k))+sin(k));U=0.5*(0.5*(exp(k)+exp(-k))-cos(k));V=0.5*(0.5*(exp(k)-exp(-k))-sin(k));a=k/l;output=[ST/aU/(a^2*E*J)V/(a^3*E*J);a*VST/(a*E*J)U/(a^2*E*J);a^2*E*J*Ua*E*J*VST/a;a^3*E*J*Ta^2*E*J*Ua*VS];Calculate.mfunctionoutput=calculate(in)symspabp=in;H1=H(0.06,1.7e-8,2.45,p);H2=H(0.15,3.2e-8,3.063,p);H3=H(0.05,3.2e-8,3.063,p);H4=H(0.03,3.2e-8,3.063,p);H5=H(0.1,0.9e-8,1.587,p);H6=H(0.06,0.9e-8,1.587,p);o1=[1a-(1-0.016/0.012)*0.012*p^2*a3.5*p^2]';o1=H1*o1;KK=[1000;0100;0010;-10^8001];o2=KK*o1;o2=H2*o2;K1=[1000;0100;00.025*p^210;7*p^2001];o3=K1*o2;o3=H3*o3;o4=KK*o3;o4=H4*o4;d=solve(o4(3));o5=subs(o4,'a',d);o5(2)=o5(2)+b;o5=H5*o5;o6=KK*o5;o6=H6*o6;K2=[1000;0100;0-(1-0.016/0.012)*0.012*p^210;3*p^2001];oo=K2*o6;RC=solve(oo(4));ee=subs(oo(3),'b',RC);output=double(ee);ee=vpa(ee);Calculate2.mfunctionoutput=calculate2(in)symspk;p=in;H1=H(0.06,1.7e-8,2.45,p);H2=H(0.15,3.2e-8,3.063,p);H3=H(0.05,3.2e-8,3.063,p);H4=H(0.03,3.2e-8,3.063,p);H5=H(0.1,0.9e-8,1.587,p);H6=H(0.06,0.9e-8,1.587,p);S=[00;00];K1=[1000;0100;0-(1-0.016/0.012)*0.012*p^210;3.5*p^2001];T22=[K1(1,1:2);K1(2,1:2)];T21=[K1(1,3:4);K1(2,3:4)];T12=[K1(3,1:2);K1(4,1:2)];T11=[K1(3,3:4);K1(4,3:4)];S1=(T11*S+T12)*inv(T21*S+T22);T22=[H1(1,1:2);H1(2,1:2)];T21=[H1(1,3:4);H1(2,3:4)];T12=[H1(3,1:2);H1(4,1:2)];T11=[H1(3,3:4);H1(4,3:4)];S2=(T11*S1+T12)*inv(T21*S1+T22);iS2=inv(S2);KK=[00;-10^50];I=eye(2);iS3=iS2*inv(I+KK*iS2);T11=[H2(1,1:2);H2(2,1:2)];T12=[H2(1,3:4);H2(2,3:4)];T21=[H2(3,1:2);H2(4,1:2)];T22=[H2(3,3:4);H2(4,3:4)];S4=(T21*iS3+T22)*inv(T11*iS3+T12);K2=[1000;0100;00.025*p^210;7*p^2001];T22=[K2(1,1:2);K2(2,1:2)];T21=[K2(1,3:4);K2(2,3:4)];T12=[K2(3,1:2);K2(4,1:2)];T11=[K2(3,3:4);K2(4,3:4)];S5=(T11*S4+T12)*inv(T21*S4+T22);T22=[H3(1,1:2);H3(2,1:2)];T21=[H3(1,3:4);H3(2,3:4)];T12=[H3(3,1:2);H3(4,1:2)];T11=[H3(3,3:4);H3(4,3:4)];S6=(T11*S5+T12)*inv(T21*S5+T22);iS6=inv(S6);iS7=iS6*inv(I+KK*iS6);T11=[H4(1,1:2);H4(2,1:2)];T12=[H4(1,3:4);H4(2,3:4)];T21=[H4(3,1:2);H4(4,1:2)];T22=[H4(3,3:4);H4(4,3:4)];S8=(T21*iS7+T22)*inv(T11*iS7+T12);KKK=[00;1/k0];S9=S8*inv(I+KKK*S8);S9=limit(S9,k,0);T22=[H5(1,1:2);H5(2,1:2)];T21=[H5(1,3:4);H5(2,3:4)];T12=[H5(3,1:2);H5(4,1:2)];T11=[H5(3,3:4);H5(4,3:4)];S10=(T11*S9+T12)*inv(T21*S9+T22);iS10=inv(S10);iS11=iS10*inv(I+KK*iS10);T11=[H6(1,1:2);H6(2,1:2)];T12=[H6(1,3:4);H6(2,3:4)];T21=[H6(3,1:2);H6(4,1:2)];T22=[H6(3,3:4);H6(4,3:4)];S12=(T21*iS11+T22)*inv(T11*iS11+T12);K3=[1000;0100;0-(1-0.016/0.012)*0.012*p^210;3*p^2001];T22=[K3(1,1:2);K3(2,1:2)];T21=[K3(1,3:4);K3(2,3:4)];T12=[K3(3,1:2);K3(4,1:2)];T11=[K3(3,3:4);K3(4,3:4)];S13=(T11*S12+T12)*inv(T21*S12+T22);ee=det(S13);output=double(ee);rotor.msymsab;a=1200;b=2300;a1=a+0.5*(b-a);f1=calculate(a);f2=calculate(b);fa=calculate(a1);fori=1:10000ifabs(fa)<1e-2disp('thecriticalspeedis');disp(a1);break;elseiff1*fa<0f2=fa;b=a1;a1=a