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

《2022年利用传递矩阵法和Riccati传递矩阵法分析转子临界转速 .pdf》由会员分享，可在线阅读，更多相关《2022年利用传递矩阵法和Riccati传递矩阵法分析转子临界转速 .pdf（13页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、利用传递矩阵法和Riccati传递矩阵法分析转子临界转速一、所需求解转子参数将转子简化为如下所示：ABC1m2m3m1J2J3J三个盘的参数为：1232221232221233.5730.0160.050.0160.0120.0250.012PPPdddmkgmkgmkgIkg mIkg mIkg mIkg mIkg mIkg m另，阶梯轴的三段轴的截面惯性矩分别为：4142431.73.20.9JcmJcmJcm三段轴的单位长度轴段的质量分别为：1232.45/3.063/1.587/mkgmmkgmmkgm二、试算转轴的传递矩阵取试算转速1200/prads ; 则，各轴段的传递矩阵分别为

2、：第 1 段840.061.7 102.45/lmJmmkgm精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 13 页1 1.0006e+000 6.0007e-002 5.2943e-007 1.0588e-008 3.7356e-002 1.0006e+000 1.7649e-005 5.2943e-007 6.3506e+003 1.2701e+002 1.0006e+000 6.0007e-002 2.1170e+005 6.3506e+003 3.7356e-002 H1.0006e+000第 2 段840.153.2 103

3、.063/lmJmmkgm2 1.0145e+000 1.5044e-001 1.7595e-006 8.7927e-008 3.8782e-001 1.0145e+000 2.3506e-005 1.7595e-006 4.9669e+004 2.4821e+003 1.0145e+000 1.5044e-001 6.6353e+005 4.9669e+004 3.8782e-001 H1.0145e+000第 3 段840.053.2 103.063/lmJmmkgm3 1.0002e+000 5.0002e-002 1.9531e-007 3.2552e-009 1.4358e-002

4、1.0002e+000 7.8128e-006 1.9531e-007 5.5135e+003 9.1890e+001 1.0002e+000 5.0002e-002 2.2054e+005 5.5135e+003 1.4358e-002 H1.0002e+000第 4 段840.033.2103.063/lmJmmkgm4 1.0000e+000 3.0000e-002 7.0313e-008 7.0313e-010 3.1013e-003 1.0000e+000 4.6875e-006 7.0313e-008 1.9848e+003 1.9848e+001 1.0000e+000 3.00

5、00e-002 1.3232e+005 1.9848e+003 3.1013e-003 H1.0000e+000第 5 段840.10.9 101.587/lmJmmkgm精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 13 页5 1.0053e+000 1.0011e-001 2.7788e-006 9.2607e-008 2.1163e-001 1.0053e+000 5.5614e-005 2.7788e-006 1.1430e+004 3.8094e+002 1.0053e+000 1.0011e-001 2.2877e+005

6、 1.1430e+004 2.1163e-001 H1.0053e+000第 6 段840.060.9 101.587/lmJmmkgm6 1.0007e+000 6.0008e-002 1.0000e-006 2.0000e-008 4.5706e-002 1.0007e+000 3.3338e-005 1.0000e-006 4.1137e+003 8.2272e+001 1.0007e+000 6.0008e-002 1.3714e+005 4.1137e+003 4.5706e-002 H1.0007e+000此 6 段传递矩阵均采用MATLAB 编程求解， MATLAB 的源文件为

7、H.m三、采用传递矩阵法进行各段轴的状态参数的传递初始参数列阵为 :0101010120101012011PddXXII pMp IQmp x令011X，则初始矩阵可化为：0101157605.04 6e以初始矩阵乘第一轴段的传递矩阵，则可得第一段轴的终端状态参数：1011011011010.06306+ 1.0541.102 +58903.088 5 26566.05.706255 6kkkkeeXMQ由于考虑支座的支撑刚度系数变化从5101*101*10，先取51*10，那么精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 13 页1

8、00001000010001KKk，此处510k，则可得支座 A后第 2 段的起始端参数阵为：020102010201560201 0.06306 + 1.0541.102 + 5890.03.088*10260.25.149*2. 0107 6XMQ用第 2 段的传递矩阵乘此矩阵，可得第2 段终端参数：20120120120166 0.2402+ 2.4721.282 + 19.11900.01.147*1099133.06.170477*1kkkkXMQ用中间圆盘的传递矩阵乘第2 段终端参数阵，即可得第3 段起始端参数：03010306103016307010.2402 + 2.472 1

9、.282 + 1958022.01.848*102.52*103.11*10.47XMQ用第 3 段传递矩阵乘其始端参数矩阵：0130130136356017 0.3238 + 3.9092.231+ 401.855*103.419*102.581*103.178*10.02kkkkXMQ用上式乘以支座刚度矩阵，得其终端参数：01305667130130130.3238 + 3.9092.231 + 40.01.855*103.419*102.549*103.1392*10kkkkXMQ则，根据可得： , 则可得支座 B后第 4 段的起始端参数阵为：560104010401040104670.

10、3238 + 3.9092.231 + 40.01.855*103.419*102.549*103.139*102XMQ同上，用此段轴的传递函数乘其起始端的状态参数，可得：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 13 页4014014014015667 0.4056 + 5.372 3.281 + 582.626*104.369*102.597*103.172*20kkkkXMQ则，根据40kM可得：01-16.64则，可得第 5 段的起始参数矩阵：050750505-1.3753.69601.12*10XMQ其中，5为铰链处的

11、转角。用第 5 段的传递矩阵乘此参数矩阵，即得第5 段的终端参数：5555565575 0.1001* - 2.05 1.005* - 2380.9*1.135*1011433.07.*1.11753* 0kkkkXMQ用上式乘以支座刚度矩阵可得第6段的初始状态参数阵：065065065076650.1001* - 2.051.005* - 2380.9*1.135*101420.0*1.17.733*10XMQ则，用第 6 段的传递矩阵乘此状态参数即可得其终端的参数阵：65656565670.1609 - 5.075 1.025 - 7960.91.827*1019300.01.178*61

12、0.99kkkkXMQ根据最右边盘得传递矩阵，可得转子终端的状态参数：67665660.1609*b - 5.075 1.025*b - 76863.0*2.27*107.144*10 *3.371*106.99kkkkbbXMQ则根据终端的自由状态，则应该60kM;60kQ精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 13 页通过令60kQ解出cR，并将其带入到6kM的表达式中，可得：-1.9463e+006；此处使用的 MATLAB 源程序为 calculate.m在 MATLAB 中使用线性插值法寻找最佳p值使得逼近于 0。 其

13、程序为 rotor.m经计算，考虑支撑刚度变化5101*101*10之间时，取51*10时，一阶临界转速值为 1.0344e+002/rads取61*10时，一阶临界转速值为3.2390e+002/rads取71*10时，一阶临界转速值为9.3859e+002/rads取81*10时，一阶临界转速值为1.8631e+003/rads取91*10时，一阶临界转速值为2.1843e+003/rads取101*10时，一阶临界转速值为2.2241e+003/rads因此随着刚度的增加， 一阶临界转速的值越来越大， 而当不考虑支座的刚度变化，假设为完全刚性的话， 一阶临界转速值为 2228.7/rad

14、s， 因此当取101*10时一阶临界转速值已相当接近完全刚性的情况。四、采用Riccati传递矩阵法进行各段轴的状态参数的传递根据 Riccati传递矩阵法的原理， 只需在传递矩阵法的基础之上求得各截面的 Riccati传递矩阵。将转子截面的状态参数分组：iiMfQ,iiXe因为左端和右端均为自由端，故00f，00e；0nf，0ne; 所以，我们可得到左端截面的Riccati传递矩阵00S；根据第 i+1 截面f、 e 之间的 Riccati变换公式：1111221221iiiST STT ST可得（同样，试算转速选为1200/prads） ：精选学习资料 - - - - - - - - -

15、名师归纳总结 - - - - - - -第 6 页，共 13 页左盘右边截面的 Riccati传递矩阵：615.04*1005760.00S；第 1 轴段末的 Riccati传递矩阵：52653.273*105.826*10 -13403.0.100273*S；刚性支承5101*101*10在此处的处理因为涉及到刚度，取51*10的情况，所以在获取其支承左边的Riccati传递矩阵2S后，需将222fSe转换为1222eSf( 即第二类 Riccati变换) ，然后再代入到普通传递矩阵的式子：32100001000100001fkfee32232ffKeee;000Kk可得：1113322eS

16、IK Sf。最终可得：(5)(7)(35)1320.0001871.072*104.3865*101.072*10eSf,即( 5)( 7)1( 5)30.00018761.072*104.385*101.072*10S;此处的13S即为刚性支承右端截面的第二类 Riccati传递矩阵。则，第 2 段轴段末的 Riccati传递矩阵：11121221112433SHSHHSH可得：545658.213*101.632*103.97*108.213*10S通过中间圆盘后，可得圆盘右边的Riccati传递矩阵：555658.213*101.992*106.11*108.213*10S;则，第 3轴

17、段右边的 Riccati 传递矩阵为：667.642*1-47877.090100.047877.00S;精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 13 页到达第二个刚性支承处， 同样采用第二类 Riccati 变换，并带入Riccati 传递矩阵公式：11111iiiSSIK S即可得：( 8)( 7)( 5)81)7(7.019*101.322*101.106*107.019*10S继续进行传递，第 4轴段末的 Riccati 传递矩阵为：11121221112877SHSHHSH可得：58651.615*107.709*10

18、1.615*160100.00S由于，第 4轴段末是球头联轴器，故，在此也要进行另外的推导，由于球头联轴器的力矩刚性系数趋近于0，则在此利用弹性铰链的传递矩阵：111000101000100001iieekff11iiiiieeKfff；000Kk；11kk把iiifSe 带入上式中，可得：11iiiSSIK S则，可得：698.143*00010S再根据第 5轴段的传递矩阵，可得第 5轴段末的 Riccati 传递矩阵：5106761.589*101.586*101.584*101.589*10S第5轴段末为一刚性支承， 则同样采用第二种 Riccati 变换，可得刚性支承右端的Riccat

19、i 传递矩阵：( 5)1110.0004185 0.004199-0.04.177*10004185S精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 13 页同样，根据第 6 轴段的传递矩阵，可得第6 轴段末的 Riccati传递矩阵：512652.501*101.473*140211.02.*100501S最后， 通过最后一个圆盘的传递， 可以得到转子最右端的Riccati传递矩阵：513652.501*102.847*10 45999.2.51001* 0S即，最右边截面上应该满足：131313fSe；又，由于1300f；所以，只有

20、当13det0S时，才能取到13e的值，所以只有临界转速值p 才能使13det S尽量接近于 0。此处；13det-1.9345e+011S此处的算法 MATLAB 程序为 calculate2.m 利用线性插值法来求取最终的临界转速值p，只需将第一次求解的rotor.m的程序中的算法程序calculate.m改为 calculate2.m即可。考虑支撑刚度变化5101*101*10之间时，取51*10时，一阶临界转速值为1.0344e+002/rads取61*10时，一阶临界转速值为3.2390e+002/rads取71*10时，一阶临界转速值为9.3859e+002/rads取81*10时

21、，一阶临界转速值为1.8631e+003/rads取91*10时，一阶临界转速值为2.1843e+003/rads取101*10时，一阶临界转速值为2.2241e+003/rads因此随着刚度的增加， 一阶临界转速的值越来越大， 而当不考虑支座的刚度变化，假设为完全刚性的话， 一阶临界转速值为 2228.7/rads， 因此当取101*10时一阶临界转速值已相当接近完全刚性的情况。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 13 页附录：H.m function output = H(in1,in2,in3,in4)% in1=l

22、in2=J in3=m in4=psyms S TU V kaE Jm pl ;E=2e11;l=in1;J=in2;m=in3;p=in4;k=l*(m*(p2)/(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=S T/a U/(a2*E*J) V/(a3*E*J);a*V S T/(a*E*J) U/(a2*E*J);a2*

23、E*J*U a*E*J*V S T/a;a3*E*J*T a2*E*J*U a*V S;Calculate.m function output=calculate(in)syms pabp=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=1 a -(1-0.016/0.012)*0.012*p2*a 3.5*p2;o1=H1*

24、o1;KK=1 0 0 0;0 1 0 0;0 0 1 0;-108 0 0 1;o2=KK*o1;o2=H2*o2;K1=1 0 0 0;0 1 0 0;0 0.025*p2 1 0;7*p2 0 0 1;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;精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 13 页o5=H5*o5;o6=KK*o5;o6=H6*o6;K2=1 0 0 0;0 1 0 0;0 -(1-0.016/

25、0.012)*0.012*p2 1 0;3*p2 0 0 1;oo=K2*o6;RC=solve(oo(4);ee=subs(oo(3),b,RC);output=double(ee);ee=vpa(ee);Calculate2.m function output=calculate2(in)syms pk;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.

26、9e-8,1.587,p);S=0 0;0 0;K1=1 0 0 0;0 1 0 0;0 -(1-0.016/0.012)*0.012*p2 1 0;3.5*p2 0 0 1;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);

27、H1(4,3:4);S2=(T11*S1+T12)*inv(T21*S1+T22);iS2=inv(S2);KK=0 0;-105 0;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=1 0 0 0;0 1 0 0;0 0.025*p2 1 0;7*p2 0 0 1;精选学习资料 - - - - - - - - - 名

28、师归纳总结 - - - - - - -第 11 页，共 13 页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);

29、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=0 0;1/k 0;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);S

30、10=(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=1 0 0 0;0 1 0 0;0 -(1-0.016/0.012)*0.012*p2 1 0;3*p2 0 0 1;T22=K3(1,1:2);K3(2,1:2);T21=K3(1,3

31、: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.m syms ab;精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 13 页a=1200;b=2300;a1=a+0.5*(b-a);f1=calculate(a);f2=calculate(b);fa=calculate(a1);for i=1:10000if abs(fa)1e-2 disp(the critical speed is); disp(a1);break;elseif f1*fa0 f2=fa; b=a1; a1=a+0.5*(b-a); fa=calculate(a1);elseif f2*fa0 f1=fa; a=a1; a1=a+0.5*(b-a); fa=calculate(a1);endendendend精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 13 页