【论坛出品】转子系统临界转速Matlab计算程序（1.1.0版）

论坛出品 · 发表于 2008-10-13 14:31

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

%=====================================================
%转子系统临界转速计算程序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 输出的临界转速，单位为转/分

论坛出品 · 发表于 2008-10-13 14:32

算例一：如图所示两端简支的光轴系统，轴长为2m，轴直径为0.1m，转子弹性模量为2.095e11Pa，转子材料密度为7.85e3kg/m3。
绘图1.jpg

很显然，该转子可以求得临界转速的理论解，理论解为：
一阶临界转速：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 = [0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25]; %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1; %各轴段外径
ShaftDO = ones(1,Nshaft)*0.0; %各轴段内径；
SPtype= 1; %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
LocationF=[1,9]; %支承所在节点编号；
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 = [0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25]; %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1; %各轴段外径
ShaftDO = ones(1,Nshaft)*0.0; %各轴段内径；
SPtype= 2; %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
LocationF=[1,9]; %支承所在节点编号；
SPstiffness=[1.0e9 1.0e9]; %支承刚度，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 = [0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25]; %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1; %各轴段外径
ShaftDO = ones(1,Nshaft)*0.0; %各轴段内径；
SPtype= 2; %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
LocationF=[1,9]; %支承所在节点编号；
SPstiffness=[2.0e7 2.0e7]; %支承刚度，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 编辑 ]

论坛出品 · 发表于 2008-10-13 14:33

算例二：如图所示转子，转子具体尺寸见程序中的输入参数，转子弹性模量为2.095e11Pa，转子材料密度为7.85e3kg/m3，转子支承刚度 5.0e7 N/m。

绘图2.jpg

计算程序如下：

clc
clear all;
Nshaft=5; %轴段数量；
%%
RotorE = 2.095e11; %转子弹性模量；
RotorM = 7.85e3; %转子材料密度
ShaftL = [0.55,0.45,0.45,0.55,0.5]; %各轴段长度
ShaftDI = ones(1,Nshaft)*0.1;ShaftDI(2)=0.15;ShaftDI(3)=0.15; %各轴段外径
ShaftDO = ones(1,Nshaft)*0.0; %各轴段内径；
LocationF=[1,5]; %支承所在节点编号；
SPtype= 2; %支承类型选择， 1 -- 刚性支承， 2 -- 弹性支承；
LocationF=[1,5]; %支承所在节点编号；
SPstiffness=[5.0e7 5.0e7]; %支承刚度，SPtype = 1时，无需进行任何设置；
%%
addtionN = [2,3,4,6]; %附加轮盘编号
addtionM = [144.1157,288.2313,144.1157,144.1157]; %附加轮盘质量
addtionJ = [4.7933,18.4648,4.7933,4.7933]; %附加轮盘转动惯量
%%
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 = [0.5,0.1,0.4,0.1,0.4,0.1,0.5,0.4,0.1]; %各轴段长度
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=[1,8]; %支承所在节点编号；
SPstiffness=[5.0e7 5.0e7]; %支承刚度，SPtype = 1时，无需进行任何设置；
%%
addtionN = [2,3,4,5,6,7,9,10]; %附加轮盘编号
addtionM = [70.1311,70.1311,144.1157,144.1157,70.1311,70.1311,70.1311,70.1311]; %附加轮盘质量
addtionJ = [2.3888,2.3888,9.2324,9.2324,2.3888,2.3888,2.3888,2.3888]; %附加轮盘转动惯量
%%
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，因此低阶临界转速的计算结果要略高一点
实际计算时，如果要获得比较准确的结果，需要考虑轮盘对转子临界转速的影响，其考虑办法可以参考西安交通大学编写的《透平零件结构和强度计算》一书中的相关内容进行修正。

对第二种分段模式节点进行加密，如图
绘图4.jpg

计算程序如下：

clc
clear all;
Nshaft=14; %轴段数量；
%%
RotorE = 2.095e11; %转子弹性模量；
RotorM = 7.85e3; %转子材料密度
ShaftL = [0.25,0.25,0.1,0.2,0.2,0.1,0.2,0.2,0.1,0.25,0.25,0.2,0.2,0.1]; %各轴段长度
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=[1,12]; %支承所在节点编号；
SPstiffness=[5.0e7 5.0e7]; %支承刚度，SPtype = 1时，无需进行任何设置；
%%
addtionN = [2,3,6,7,9,10,14,15]; %附加轮盘编号
addtionM = [70.1311,70.1311,144.1157,144.1157,70.1311,70.1311,70.1311,70.1311]; %附加轮盘质量
addtionJ = [2.3888,2.3888,9.2324,9.2324,2.3888,2.3888,2.3888,2.3888]; %附加轮盘转动惯量
%%
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 编辑 ]

jgwang · 发表于 2008-10-25 13:23

顶一下,再慢慢品位

oiloil · 发表于 2008-10-27 11:05

请问程序中的转动惯量是极转动惯量Jp=1/2mR^2，还是过轴心的直径转动惯量Jd=1/2Jp？

oiloil · 发表于 2008-10-27 20:56

感觉不错，可是还是不会用啊，是不是要建立一个代码那样m文件阿？然后如何调用p文件呢？操作了一下，没有得到结果

tcy003 · 发表于 2008-11-12 14:51

楼主，你的这个.P文件用什么工具打开啊？

jinfengshi_2005 · 发表于 2008-12-1 13:12

不知道有什么用！

wbl617 · 发表于 2008-12-5 14:59

我正在研究阶梯型转轴的临界转速，现在进展很缓慢，到处搜集资料，今天有幸看见该贴，麻烦楼主把源程序给我看一下，邮箱83886895@qq.com，谢谢了！！

[ 本帖最后由 wbl617 于 2008-12-5 15:00 编辑 ]

余方 · 发表于 2008-12-7 13:26

:lol :@) :handshake

gujianhai1985 · 发表于 2009-1-14 10:59

请问程序中的转动惯量是极转动惯量Jp=1/2mR^2，还是过轴心的直径转动惯量Jd=1/2Jp？

shanepeng · 发表于 2009-1-14 13:19

請教一下,這裡的支撐剛性只考慮單一垂直之承方向嗎??
是否需要考慮水準剛性呢??
因為油膜軸承的水準和垂直剛性會有所不同

zz80 · 发表于 2009-2-17 21:36

是matlab程序吗？

xxgy2000 · 发表于 2009-2-19 23:23

要仔细研究下，不过先谢谢哈！

haobaobaos · 发表于 2009-2-21 16:12

:victory: :victory: :victory: :victory: 超好的资料

账号		自动登录	找回密码
密码			我要加入

[转子动力学] 【论坛出品】转子系统临界转速Matlab计算程序（1.1.0版）

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

本帖被以下淘专辑推荐:

点评

看得稀里糊涂