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

利用传递矩阵法和Riccati 传递矩阵法分析转子临界转速
一、
所需求解转子参数
将转子简化为如下所示：
三个盘的参数为：1232
2212322
2
1
230.0160.050.0160.0120.0250.012P P P d d d I kg m I kg m I kg m I kg m I kg m I kg m ⎪
=⋅=⋅=⋅⎨⎪=⋅=⋅=⋅⎩ 另，阶梯轴的三段轴的截面惯性矩分别为： 414243
1.73.20.9J cm J cm J cm ⎧=⎪
=⎨⎪=⎩
三段轴的单位长度轴段的质量分别为：123
2.45/
3.063/1.587/m kg m m kg m m kg m =⎧⎪
=⎨⎪=⎩
二、 试算转轴的传递矩阵
取试算转速1200/p rad s ω== ; 则，各轴段的传递矩阵分别为： 第1段
840.061.7102.45/l m J m m kg m -=⎧⎪=⨯⎨⎪=⎩
1 1.0006e+000 6.0007e-00
2 5.2943e-007 1.0588e-008 3.7356e-002 1.0006e+000 1.7649e-005 5.2943e-007 6.3506e+00
3 1.2701e+002 1.0006e+000 6.0007e-002 2.1170e+005 6.3506e+003 3.7356e-002 H = 1.0006e+000
⎧⎪⎪⎨⎪⎪⎩
第2段
840.153.2103.063/l m J m m kg m -=⎧⎪=⨯⎨⎪=⎩
2 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+00
3 1.0145e+000 1.5044e-001 6.6353e+005 4.9669e+00
4 3.8782e-001 H = 1.0145e+000
⎧⎪⎪⎨⎪⎪⎩
第3段
840.053.2103.063/l m J m m kg m -=⎧⎪=⨯⎨⎪=⎩
3 1.0002e+000 5.0002e-002 1.9531e-007 3.2552e-009 1.4358e-002 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 H = 1.0002e+000
⎧⎪⎪⎨⎪⎪⎩
第4段
840.033.2103.063/l m J m m kg m -=⎧⎪=⨯⎨⎪=⎩
4 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.0000e-002 1.3232e+00
5 1.9848e+003 3.1013e-003 H = 1.0000e+000
⎧⎪⎪⎨⎪⎪⎩
第5段
840.10.9101.587/l m J m m kg m -=⎧⎪=⨯⎨⎪=⎩
5 1.0053e+000 1.0011e-001 2.7788e-00
6 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 1.1430e+004 2.1163e-001 H = 1.0053e+000
⎧⎪⎪⎨⎪⎪⎩
第6段
840.060.9101.587/l m J m m kg m -=⎧⎪=⨯⎨⎪=⎩
6 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 H = 1.0007e+000
⎧⎪⎪⎨⎪⎪⎩
此6段传递矩阵均采用MATLAB 编程求解，MATLAB 的源文件为H.m 三、
采用传递矩阵法进行各段轴的状态参数的传递
初始参数列阵为:0101010120101012011P
d d X X I I p M p I Q mp x θθωθ⎛⎫
⎛⎫ ⎪
⎪ ⎪ ⎪ ⎪=⎛⎫ ⎪-- ⎪
⎪ ⎪ ⎪
⎝⎭⎝⎭ ⎪⎝⎭ 令011X =，则初始矩阵可化为：010*******.046e θθ⎛⎫
⎪
⎪ ⎪ ⎪⎝⎭
以初始矩阵乘第一轴段的传递矩阵，则可得第一段轴的终端状态参数：
1011011011010.06306+ 1.0541.102 +5890 3.0885 26566.0 5..7062556k k k k e e X M Q θθθθθ⎛⎫⎛⎫
⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ++⎪⎝⎭⎝⎭
由于考虑支座的支撑刚度系数变化从5101*101*10，先取51*10，那么
100001000010001KK k ⎡⎤⎢⎥⎢
⎥=⎢⎥⎢⎥
-⎣⎦
，此处510k =，则可得支座A 后第2段的起始端参数阵为：
0201020102015
60201 0.06306 + 1.0541.102 + 5890.0 3.088*10260.2 5.149*2.01076X M Q θθθθθ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ++⎪⎝⎭⎝⎭
用第2段的传递矩阵乘此矩阵，可得第2段终端参数：
2012012012016
6 0.2402+ 2.4721.282 + 19.11900.0 1.147*1099133.0 6.170477*1k k k k X M Q θθθθθ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭
++
用中间圆盘的传递矩阵乘第2段终端参数阵，即可得第3段起始端参数：
03010306
103016307010.2402 + 2.472 1.282 + 1958022.0 1.848*102.52*10 3.11*10.47X M Q θθθθθ+⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎭⎝⎭+⎪⎝
用第3段传递矩阵乘其始端参数矩阵：
'013'013'013'6356017 0.3238 + 3.9092.231+ 401.855*10 3.419*102.581*10 3.178*10.02k k k k X M Q θθθθθ⎛⎫⎛⎫ ⎪ ⎪
⎪ ⎪= ⎪ ⎪
⎪ ⎪
⎪⎝
⎭++⎝⎭
用上式乘以支座刚度矩阵，得其终端参数：
013056
67130130130.3238 + 3.9092.231 + 40.01.855*10 3.419*102.549*10 3.1392*10k k k k X M Q θθθθθ⎛⎫⎛⎫
⎪ ⎪
⎪ ⎪= ⎪ ⎪ ⎪ ⎪
⎝
⎭⎭+⎪+⎝
则，根据可得： ,则可得支座B 后第4段的起始端参数阵为：
56
0104010401040104670.3238 + 3.9092.231 + 40.01.855*10 3.419*102.549*10 3.139*102X M Q θθθθθ⎛⎫⎛⎫
⎪ ⎪
⎪ ⎪= ⎪ ⎪ ⎪ ⎪
⎪⎝
⎭⎝++⎭
同上，用此段轴的传递函数乘其起始端的状态参数，可得：
40140140140156
67 0.4056 + 5.372 3.281 + 582.626*10 4.369*102.597*10 3..172*20k k k k X M Q θθθθθ+⎛⎫⎛⎫
⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭+
则，根据40k M =可得：01-16.64θ= 则，可得第5段的起始参数矩阵：
050750505-1.3753.69601.12*10X M Q θ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎭⎝-⎪⎝
⎭ 其中，5θ为铰链处的转角。

用第5段的传递矩阵乘此参数矩阵，即得第5段的终端参数：
555556
5575 0.1001* - 2.05 1.005* - 2380.9* 1.135*1011433.07.* 1.11753*0k k k k X M Q θθθθθ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭
--
用上式乘以支座刚度矩阵可得第6段的初始状态参数阵：
0650650650766
50.1001* - 2.051.005* - 2380.9* 1.135*101420.0* 1.17.733*10X M Q θθθθθ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭
--
则，用第6段的传递矩阵乘此状态参数即可得其终端的参数阵：
656565656
70.1609 - 5.075 1.025 - 7960.9 1.827*1019300.0 1.178*610.99k k k k X M Q θθθθθ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭
--
根据最右边盘得传递矩阵，可得转子终端的状态参数：
6766
5660.1609*b - 5.075 1.025*b - 76863.0* 2.27*107.144*10* 3.371*106.99k k k k b b X M Q θ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝
⎭-⎝-⎭
则根据终端的自由状态，则应该60k M =;60k Q =
通过令60k Q =解出c R ，并将其带入到6k M 的表达式中，可得：
-1.9463e+006∆=；
此处使用的MATLAB 源程序为calculate.m
在MATLAB 中使用线性插值法寻找最佳p 值使得逼近于0。

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

四、
采用Riccati 传递矩阵法进行各段轴的状态参数的传递
根据Riccati 传递矩阵法的原理，只需在传递矩阵法的基础之上求得各截面的Riccati 传递矩阵。

将转子截面的状态参数分组：{}i i
M f Q ⎧⎫
=⎨⎬⎩⎭,{}i i X e θ⎧⎫=⎨⎬⎩⎭
因为左端和右端均为自由端，故{}00f =，{}00e ≠；{}0n f =，{}0n e ≠; 所以，我们可得到左端截面的Riccati 传递矩阵[]00S =；
根据第i+1截面{}f 、{}e 之间的Riccati 变换公式：
[][][]1
111221221i i i
S T S T T S T -+=++
可得（同样，试算转速选为1200/p rad s =）：
左盘右边截面的Riccati 传递矩阵：[]6
1 5.04*1005760.00S ⎡
⎤
=⎢⎥⎣⎦； 第1轴段末的Riccati 传递矩阵：[]5
2
653.273*105. 826*10
-13403.0.100273*S -⎡⎤
=⎢⎥⎣⎦
； 刚性支承5101*101*10在此处的处理因为涉及到刚度，取51*10的情况，所以在获取其支承左边的Riccati 传递矩阵[]2S 后，需将{}[]{}222f S e =转换为
{}[]{}1
222e S f -=(即第二类Riccati 变换)，然后再代入到普通传递矩阵的式子：
3210000100
0100
01f k
f e e ⎡⎤⎢
⎥-⎧⎫⎧⎫⎢⎥=⎨⎬⎨⎬⎢⎥⎩⎭⎩⎭⎢⎥⎣⎦
{}{}{}{}{}322
32f f K e e e =+⎧⎪⎨=⎪⎩
;000K k ⎡⎤
=⎢⎥-⎣⎦
可得：{}[][]{
}
{}1
1
13322e S I K S f ---=+。

最终可得：
{}[]{}(5)
(7)(35)1
32
0.0001871.072*10 4.3865*101.072*10e S f ------⎡⎤
=⎢
⎥⎣⎦
,即
[](5)
(7)1
(5)3
0.00018761.072*10 4.385*101.072*10S ------⎡⎤=⎢
⎥⎣⎦;
此处的[]1
3
S -即为刚性支承右端截面的第二类Riccati 传递矩阵。

则，第2段轴段末的Riccati 传递矩阵：
[][]()[](
)
1
1
1
21221112
43
3S H S H H S H ---=++
可得：[]54
5658.213*10 1.632*103.97*10
8.213*10S ⎡⎤=-⎢⎣-⎥⎦ 通过中间圆盘后，可得圆盘右边的Riccati 传递矩阵：
[]55
5658.213*10 1.992*106.11*10
8.213*10S ⎡⎤
=⎢⎣-⎥⎦; 则，第3轴段右边的Riccati 传递矩阵为：
[]6
67.642*1-47877.090100.047877.00S ⎡⎤
=⎢⎥⎣⎦
;
到达第二个刚性支承处，同样采用第二类Riccati 变换，并带入Riccati 传递矩阵公式：
[][][]{}
1
11
11i i
i
S S I K S ----+=+
即可得：[](8)(7)(5)81
)7
(7.019*10 1.322*101.106*107.019*10S -----⎡⎤
=⎢⎥⎣⎦
- 继续进行传递，第4轴段末的Riccati 传递矩阵为：
[][]()[](
)
1
1
1
21221112
87
7S H S H H S H ---=++
可得：[]5
86
51.615*107.709*10 1.615*160100.00S ⎡⎤
=⎢⎥⎣⎦
-
由于，第4轴段末是球头联轴器，故，在此也要进行另外的推导，由于球头联轴器的力矩刚性系数趋近于0，则在此利用弹性铰链的传递矩阵：
11100
010100
0100
1i i e e k f f +⎡⎤⎢⎥⎢⎥⎧⎫
⎧⎫=⎢⎥⎨⎬⎨⎬⎢⎥⎩⎭⎩⎭⎢⎥⎢⎥⎣⎦
{}{}{}{}{}11i i i
i i e e K f f f ++=+⎧⎪⎨=⎪⎩
；000K k ⎡⎤=⎢⎥⎣⎦；11k k = 把{}[]{}i i i f S e =带入上式中，可得： [][][]()1
1i i i S S I K S -+=+
则，可得： []6
98.143*00010S ⎡
⎤=⎢⎥⎣⎦
再根据第5轴段的传递矩阵，可得第5轴段末的Riccati 传递矩阵：
[]510
6761.589*10 1.586*101.584*10 1.589*10S ⎡⎤
=⎢⎥⎣-⎦
- 第5轴段末为一刚性支承，则同样采用第二种Riccati 变换，可得刚性支承右端的Riccati 传递矩阵：
[](5)1
11
0.0004185 0.004199
-0.04.177*10004185S --⎡⎤=⎢⎣-⎥⎦
同样，根据第6轴段的传递矩阵，可得第6轴段末的Riccati 传递矩阵：
[]512
652.501*101.473*140211.0 2.*100501S --⎡⎤
=⎢⎥⎣⎦
最后，通过最后一个圆盘的传递，可以得到转子最右端的Riccati 传递矩阵：
[]5
13
652.501*102.847*10 45999.2.51001*0S -⎡⎤
=⎢⎥⎣⎦
即，最右边截面上应该满足：{}[]{}131313f S e =；又，由于{}1300f ⎧⎫
=⎨⎬⎩⎭；所
以，只有当[]13det 0S =时，才能取到{}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/rad s
取61*10时，一阶临界转速值为3.2390e+002/rad s 取71*10时，一阶临界转速值为9.3859e+002/rad s 取81*10时，一阶临界转速值为1.8631e+003/rad s 取91*10时，一阶临界转速值为2.1843e+003/rad s 取101*10时，一阶临界转速值为2.2241e+003/rad s
因此随着刚度的增加，一阶临界转速的值越来越大，而当不考虑支座的刚度变化，假设为完全刚性的话，一阶临界转速值为2228.7/rad s ，因此当取101*10时一阶临界转速值已相当接近完全刚性的情况。

附录：
H.m
function output = H(in1,in2,in3,in4)
% in1=l in2=J in3=m in4=p
syms S T U V k a E J m p l;
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=[S T/a U/(a^2*E*J) V/(a^3*E*J);a*V S T/(a*E*J)
U/(a^2*E*J);a^2*E*J*U a*E*J*V S T/a;a^3*E*J*T a^2*E*J*U a*V S];
Calculate.m
function output=calculate(in)
syms p a b
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);
o1=[1 a -(1-0.016/0.012)*0.012*p^2*a 3.5*p^2]';
o1=H1*o1;
KK=[1 0 0 0;0 1 0 0;0 0 1 0;-10^8 0 0 1];
o2=KK*o1;
o2=H2*o2;
K1=[1 0 0 0;0 1 0 0;0 0.025*p^2 1 0;7*p^2 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;
o5=H5*o5;
o6=KK*o5;
o6=H6*o6;
K2=[1 0 0 0;0 1 0 0;0 -(1-0.016/0.012)*0.012*p^2 1 0;3*p^2 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 p k;
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=[0 0;0 0];
K1=[1 0 0 0;0 1 0 0;0 -(1-0.016/0.012)*0.012*p^2 1 0;3.5*p^2 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);H1(4,3:4)];
S2=(T11*S1+T12)*inv(T21*S1+T22);
iS2=inv(S2);
KK=[0 0;-10^5 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*p^2 1 0;7*p^2 0 0 1];
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=[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)];
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=[1 0 0 0;0 1 0 0;0 -(1-0.016/0.012)*0.012*p^2 1 0;3*p^2 0 0 1];
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.m
syms a b;
a=1200;
b=2300;
a1=a+0.5*(b-a);
f1=calculate(a);
f2=calculate(b);
fa=calculate(a1);
for i=1:10000
if abs(fa)<1e-2
disp('the critical speed is'); disp(a1);
break;
else if f1*fa<0
f2=fa;
b=a1;
a1=a+0.5*(b-a);
fa=calculate(a1);
else if f2*fa<0
f1=fa;
a=a1;
a1=a+0.5*(b-a);
fa=calculate(a1);
end
end
end
end。