【论坛出品】转子系统临界转速Matlab计算程序（1.1.0版）
%=====================================================%转子系统临界转速计算程序1.1.0版；
%本程序暂时不考虑陀螺效应的影响；
%程序的后续版本将逐步补充完善相关功能。
%=====================================================
%本程序由中国振动联盟网站(http://forum.vibunion.com)开发；
%本程序版权属中国振动联盟所有，未经中国振动联盟网站同意，禁止将本程序用于任何商业行为；
%本程序仅供联盟会员学习和验算使用，中国振动联盟不对程序的正确性作任何承诺；
%感谢yejet在程序编写过程中做的大量的工作，感谢心灯提供的参考程序；
%转载本程序请注明程序来源：中国振动联盟网http://forum.vibunion.com。
%=====================================================
本程序正在逐步开发中，本站将根据开发进度，阶段性的发布最新版本的程序
由于多方面的原因，源代码不能公布，如果有特殊需要，大家可以回帖说明，我们将在能力范围内为大家扩展相应的功能
1.1.0 版在1.0.0版的基础上对支承系统的处理进行了改进
1.0.0 版见帖子http://forum.vibunion.com/forum/thread-70692-1-1.html
改进后，程序可以根据需要选择刚性支承或者弹性支承。
调用格式：
CriticalSpeeds=Chinavib_CriticalSpeeds(Nshaft,RotorE,RotorM,ShaftL,ShaftDI,ShaftDO,SPtype,LocationF,SPstiffness,addtionN,addtionM,addtionJ,CSN)
参数说明：（所有输入单位均采用国际制单位）
Nshaft%轴段数量；
%%
RotorE %转子弹性模量；
RotorM %转子材料密度
ShaftL %各轴段长度
ShaftDI %各轴段外径
ShaftDO%各轴段内径；
%%
SPtype %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
LocationF %支承所在节点编号；
SPstiffness %支承刚度，SPtype = 1时，无需进行任何设置；
%%
addtionN%附加轮盘编号
addtionM%附加轮盘质量
addtionJ %附加轮盘转动惯量
CSN %输出的临界转速数量
CriticalSpeeds 输出的临界转速，单位为 转/分 算例一：如图所示两端简支的光轴系统，轴长为2m，轴直径为0.1m，转子弹性模量为2.095e11Pa，转子材料密度为7.85e3kg/m3。
很显然，该转子可以求得临界转速的理论解，理论解为：
一阶临界转速：3043转/分
二阶临界转速：12172 转/分
三阶临界转速：27387 转/分
四阶临界转速：48689 转/分
五阶临界转速：76076 转/分
刚性支承下本程序的计算结果为：
一阶临界转速：3041转/分
二阶临界转速：12138 转/分
三阶临界转速：27234 转/分
四阶临界转速：48289 转/分
五阶临界转速：75341 转/分
计算结果和理论解基本完全吻合，尤其是低阶临界转速
从整体上看，结果本程序计算结果略低于理论解
这是由于刚性支承的处理方式造成的，今后将改进相关的刚度处理方式，以进一步提高解的准确性
该算例具体计算方法如下：
将该光轴等分成8段，显然个轴段长度为0.25m，各轴段外径为0.1m，各轴段内经为0m，支承所在节点为1，9号节点
然后将各轴段的相关参数输入程序，并调用本程序，具体如下：clc
clear all;
Nshaft=8; %轴段数量；
%%
RotorE = 2.095e11; %转子弹性模量；
RotorM = 7.85e3; %转子材料密度
ShaftL = ; %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1; %各轴段外径
ShaftDO = ones(1,Nshaft)*0.0; %各轴段内径；
SPtype= 1; %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
LocationF=; %支承所在节点编号；
SPstiffness=[]; %支承刚度，SPtype = 1时，无需进行任何设置；
%%
addtionN = []; %附加轮盘编号
addtionM = []; %附加轮盘质量
addtionJ = []; %附加轮盘转动惯量
%%
CSN = 5; %输出的临界转速数量
CriticalSpeeds=Chinavib_CriticalSpeeds(Nshaft,RotorE,RotorM,ShaftL,ShaftDI,ShaftDO,SPtype,LocationF,SPstiffness,addtionN,addtionM,addtionJ,CSN)支承刚度为 1e9时，计算程序和计算结果如下：clc
clear all;
Nshaft=8; %轴段数量；
%%
RotorE = 2.095e11; %转子弹性模量；
RotorM = 7.85e3; %转子材料密度
ShaftL = ; %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1; %各轴段外径
ShaftDO = ones(1,Nshaft)*0.0; %各轴段内径；
SPtype= 2; %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
LocationF=; %支承所在节点编号；
SPstiffness=; %支承刚度，SPtype = 1时，无需进行任何设置；
%%
addtionN = []; %附加轮盘编号
addtionM = []; %附加轮盘质量
addtionJ = []; %附加轮盘转动惯量
%%
CSN = 5; %输出的临界转速数量
CriticalSpeeds=Chinavib_CriticalSpeeds(Nshaft,RotorE,RotorM,ShaftL,ShaftDI,ShaftDO,SPtype,LocationF,SPstiffness,addtionN,addtionM,addtionJ,CSN)CriticalSpeeds =
1.0e+004 *
0.3033
1.2016
2.6620
4.6314
7.0321
支承刚度为 2e7时，计算程序和计算结果如下：clc
clear all;
Nshaft=8; %轴段数量；
%%
RotorE = 2.095e11; %转子弹性模量；
RotorM = 7.85e3; %转子材料密度
ShaftL = ; %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1; %各轴段外径
ShaftDO = ones(1,Nshaft)*0.0; %各轴段内径；
SPtype= 2; %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
LocationF=; %支承所在节点编号；
SPstiffness=; %支承刚度，SPtype = 1时，无需进行任何设置；
%%
addtionN = []; %附加轮盘编号
addtionM = []; %附加轮盘质量
addtionJ = []; %附加轮盘转动惯量
%%
CSN = 5; %输出的临界转速数量
CriticalSpeeds=Chinavib_CriticalSpeeds(Nshaft,RotorE,RotorM,ShaftL,ShaftDI,ShaftDO,SPtype,LocationF,SPstiffness,addtionN,addtionM,addtionJ,CSN)CriticalSpeeds =
1.0e+004 *
0.2703
0.7849
1.3099
2.2146
3.8483
[ 本帖最后由 论坛出品 于 2008-10-13 14:58 编辑 ] 算例二：如图所示转子，转子具体尺寸见程序中的输入参数，转子弹性模量为2.095e11Pa，转子材料密度为7.85e3kg/m3，转子支承刚度 5.0e7 N/m。
计算程序如下：
clc
clear all;
Nshaft=5; %轴段数量；
%%
RotorE = 2.095e11; %转子弹性模量；
RotorM = 7.85e3; %转子材料密度
ShaftL = ; %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1;ShaftDI(2)=0.15;ShaftDI(3)=0.15; %各轴段外径
ShaftDO = ones(1,Nshaft)*0.0; %各轴段内径；
LocationF=; %支承所在节点编号；
SPtype= 2; %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
LocationF=; %支承所在节点编号；
SPstiffness=; %支承刚度，SPtype = 1时，无需进行任何设置；
%%
addtionN = ; %附加轮盘编号
addtionM = ; %附加轮盘质量
addtionJ = ; %附加轮盘转动惯量
%%
CSN = 5;
CriticalSpeeds=Chinavib_CriticalSpeeds(Nshaft,RotorE,RotorM,ShaftL,ShaftDI,ShaftDO,SPtype,LocationF,SPstiffness,addtionN,addtionM,addtionJ,CSN)
计算结果如下：
CriticalSpeeds =
1.0e+004 *
0.1314
0.2050
0.4457
0.8980
1.2281
改变分段方式如下：
计算程序如下：clc
clear all;
Nshaft=9; %轴段数量；
%%
RotorE = 2.095e11; %转子弹性模量；
RotorM = 7.85e3; %转子材料密度
ShaftL = ; %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1;ShaftDI(2)=0.15;ShaftDI(3)=0.15; %各轴段外径
ShaftDI(4)=0.15;ShaftDI(5)=0.15;ShaftDI(6)=0.15;
ShaftDO = ones(1,Nshaft)*0.0; %各轴段内径；
SPtype= 2; %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
LocationF=; %支承所在节点编号；
SPstiffness=; %支承刚度，SPtype = 1时，无需进行任何设置；
%%
addtionN = ; %附加轮盘编号
addtionM = ; %附加轮盘质量
addtionJ = ; %附加轮盘转动惯量
%%
CSN = 5;
CriticalSpeeds=Chinavib_CriticalSpeeds(Nshaft,RotorE,RotorM,ShaftL,ShaftDI,ShaftDO,SPtype,LocationF,SPstiffness,addtionN,addtionM,addtionJ,CSN)计算结果如下：
CriticalSpeeds =
1.0e+004 *
0.1378
0.2166
0.4518
0.9169
1.1840
两种不同分段形式所得到的结果不同，主要是因为分段时考虑轴段刚度不一致造成的
在第一种分段模式中，在2节点处轴段左侧轴段直径为0.1m，第4节点处右侧的轴段直径为0.1m
而第二种分段模式中，上述两个位置轴段的直径为0.15m，因此低阶临界转速的计算结果要略高一点
实际计算时，如果要获得比较准确的结果，需要考虑轮盘对转子临界转速的影响，其考虑办法可以参考西安交通大学编写的《透平零件结构和强度计算》一书中的相关内容进行修正。
对第二种分段模式节点进行加密，如图
计算程序如下：clc
clear all;
Nshaft=14; %轴段数量；
%%
RotorE = 2.095e11; %转子弹性模量；
RotorM = 7.85e3; %转子材料密度
ShaftL = ; %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1;ShaftDI(3)=0.15;ShaftDI(4)=0.15; %各轴段外径
ShaftDI(5)=0.15;ShaftDI(6)=0.15;ShaftDI(7)=0.15;ShaftDI(8)=0.15;ShaftDI(9)=0.15;
ShaftDO = ones(1,Nshaft)*0.0; %各轴段内径；
SPtype= 2; %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
LocationF=; %支承所在节点编号；
SPstiffness=; %支承刚度，SPtype = 1时，无需进行任何设置；
%%
addtionN = ; %附加轮盘编号
addtionM = ; %附加轮盘质量
addtionJ = ; %附加轮盘转动惯量
%%
CSN = 5;
CriticalSpeeds=Chinavib_CriticalSpeeds(Nshaft,RotorE,RotorM,ShaftL,ShaftDI,ShaftDO,SPtype,LocationF,SPstiffness,addtionN,addtionM,addtionJ,CSN)计算结果如下：
CriticalSpeeds =
1.0e+004 *
0.1400
0.2180
0.4441
0.8866
1.0966
从结果上看，低阶临界转速计算结果略有增加，高阶临街转速略有降低，但是影响不大
[ 本帖最后由 论坛出品 于 2008-10-13 14:50 编辑 ]
看得稀里糊涂
顶一下,再慢慢品位 请问程序中的转动惯量是极转动惯量Jp=1/2mR^2，还是过轴心的直径转动惯量Jd=1/2Jp？ 感觉不错，可是还是不会用啊，是不是要建立一个代码那样m文件阿？然后如何调用p文件呢？操作了一下，没有得到结果 楼主，你的这个.P文件用什么工具打开啊？ 不知道有什么用！ 我正在研究阶梯型转轴的临界转速，现在进展很缓慢，到处搜集资料，今天有幸看见该贴，麻烦楼主把源程序给我看一下，邮箱83886895@qq.com，谢谢了！！[ 本帖最后由 wbl617 于 2008-12-5 15:00 编辑 ] :lol :@) :handshake 请问程序中的转动惯量是极转动惯量Jp=1/2mR^2，还是过轴心的直径转动惯量Jd=1/2Jp？ 請教一下,這裡的支撐剛性只考慮單一垂直之承方向嗎??
是否需要考慮水準剛性呢??
因為油膜軸承的水準和垂直剛性會有所不同 是matlab程序吗？ 要仔细研究下，不过先谢谢哈！ :victory: :victory: :victory: :victory: 超好的资料