基于ANSYS经典界面的偏心圆盘启动的转子动力学分析－２

这是＜基于ANSYS经典界面的偏心圆盘启动的转子动力学分析－１＞的续篇，为方便阅读，重抄题目如下．
【问题描述】
一个转子如下图所示。

X
该转子两端简支，在距离下端三分之一处有一个轴承支撑；而在距离上端三分之一处有一个刚性圆盘，该圆盘存在一个偏心质量。
现在该转轴从静止开始转动，经过４秒转速达到５０００ＲＰＭ。要求对该转子的启动过程进行仿真，考察偏心圆盘所在处的轴位移和弯曲应力随时间变化的过程。
其它已知条件如下：
（１）几何参数
轴的长度: 0.4 m
轴的半径: 0.01 m
（２）圆盘的惯性参数:
质量= 16.47 kg
惯性矩 (XX,YY) = 9.47e-2 kg.m²
惯性矩 (ZZ) = 0.1861 kg.m²
偏心质量 (0.1g)位于圆盘上，距离轴心的距离为0.15 m .
（３）材料参数：
弹性模量(E) = 2.0e+11 N/m²
泊松比 (υ) = 0.3
密度 = 7800 kg/m³
《注》该算例来自于ANSYS APDL转子动力学部分的帮助实例。
【问题分析】
１.这是一个瞬态动力学问题。
２.该瞬态动力学问题要考虑科里奥利效应，因此在分析中要打开coriolis开关。
３.使用ＢＥＡＭ １８８建模转子，使用ＭＡＳＳ２１建模圆盘，使用两个一维的弹簧单元ＣＯＭＢＩＮ１４建模轴承，对这两个弹簧设置不同的刚度，以模拟非轴对称的轴承。
４.轴上所有节点限制轴向位移以及围绕轴转动的位移。轴的两个端点进一步限制住所有平移，而代表轴承与机架连接的ＣＯＭＢＩＮ１４的机架端节点位移全固定。
５.使用数组存储不平衡力的两个分量。这里要使用一系列中间变量。具体算法请参见＜基于ANSYS经典界面的偏心圆盘启动的转子动力学分析－１＞。
６.做瞬态动力学分析，在加载时，对节点施加集中力，并把不平衡力数组设置为其参数，从而施加变化的力。
７.使用时间历程后处理器，首先提取所需节点的位移和弯曲应力，然后绘制其图形。
８.本文使用命令流的方式进行讲解。
【求解步骤】
一 建模
１.定义参数
在命令窗口中输入
length = 0.4
ro_shaft = 0.01
ro_disk = 0.15
md = 16.47
id = 9.427e-2
ip = 0.1861
kxx = 2.0e+5
kyy = 5.0e+5
beta = 2.e-4
这里定义了轴的长度等变量，这些变量在下面会用到。至于其具体含义，读者在下面的命令中发现有不知道含义的变量出现时，可以返回来查找，就可以确定其确切含义。
命令运行后，在参数对话框中显示如下

２.定义材料模型
继续输入
/prep7
mp,ex,1,2.0e+11
mp,nuxy,1,.3
mp,dens,1,7800
其中，
第１行：进入到前处理器
第２－４行：定义材料的弹性模量，泊松比和密度。
３.定义单元类型
（１）定义建模轴的单元
继续输入
et,1,188
sect,1,beam,csolid
secdata,ro_shaft,20
其中，
第１行：定义单元类型是ＢＡＥＭ１８８
第２行：设置该单元类型是实心圆截面
第３行：定义该圆截面的半径是ro_shaft。查阅前面的参数，知道ro_shaft = 0.01 ，即其半径是１０ｍｍ
（２）定义建模圆盘的质量单元
继续输入
et,2,21
r,2,md,md,md,id,id,ip
其中，
第１行：定义单元类型是ＭＡＳＳ２１
第２行：设置该质量单元的惯性量，其中参数md,id,ip 的具体数据请参见前面的参数定义。
（３）定义建模轴的两中ＣＯＭＢＩＮ１４单元
继续输入
et,3,14,,1
r,3,kxx,beta*kxx
其中，
第１行：定义单元类型是ＣＯＭＢＩＮ１４.
第２行：定义该单元的刚度和阻尼。刚度是kxx；阻尼是beta*kxx 。其大小请参见前面的参数定义。
继续输入
et,4,14,,2
r,4,kyy,beta*kyy
其中，
第１行：定义单元类型是ＣＯＭＢＩＮ１４.
第２行：定义该单元的刚度和阻尼。刚度是kyy；阻尼是beta*kyy。其大小请参见前面的参数定义。
之所以要定义两个ＣＯＭＢＩＮ１４单元，是因为该轴承两个方向的刚度和阻尼不相同。
此时如果打开实常数对话框，可以看到

可见，已经输入了３个实常数。
４.创建轴
继续输入
type,1
secn,1
mat,1
k,1
k,2,,,length
l,1,2
lesize,1,,,9
lmesh,all
其中，
第１－３行：设定下面要用到的单元类型，截面以及材料。
第４－５行：创建轴的两个端点
第６行：连接这两个端点，得到一条直线
第７行：指定将该直线分为９个单元
第８行：划分网格。
为了能够看清楚结果
继续输入
/PNUM,KP,1
/PNUM,LINE,1
/PNUM,NODE,1
/PNUM,ELEM,1
/ANG,1
/REP,FAST
其中，
第１－４行：说明打开关键点，直线，节点和单元的编号开关
第５－６行：说明查看仰视图，并重新绘制。
上述命令执行后，主窗口显示如下

可见，该直线已经被划分为９个单元，它有１０个节点。
５.创建圆盘
继续输入
type,2
real,2
e,5
其含义是：选择２号质量单元，及其对应的实常数，在５号节点上创建一个质量单元。
命令执行后，主窗口显示如下

可见，５号节点处又多了一个１０，它意味着１０号单元创建在此处。
６.创建轴承
继续输入
n,21,-0.05,,2*length/3
type,3
real,3
e,8,21
其含义是；定义一个节点２１号，然后在８－２１号节点之间创建一个弹簧单元ＣＯＭＢＩＮ１４，并使用３号实常数。
命令执行后，主窗口显示如下

可见，一个新的弹簧单元已经创建。
继续输入
type,4
real,4
e,8,21
其含义是；在８－２１号节点之间再创建一个弹簧单元ＣＯＭＢＩＮ１４，并使用４号实常数。
命令执行后，主窗口显示如下

可见，在１１上面又覆盖了１２这个数据，意味着这两个单元位置是重合的。
虽然位置重合，但是其方向是不同的。
７.施加位移边界条件
继续输入
dk,1,ux,,,,uy
dk,2,ux,,,,uy
其含义是：对于轴的两个端点固定其ＵＸ，ＵＹ两个自由度。
命令执行后，主窗口显示如下

继续输入
d,all,uz
d,all,rotz
其含义是：固定所有节点的ＵＺ，ＲＯＴＺ自由度。
命令执行后，主窗口显示如下

继续输入
d,21,all
finish
其含义是：把轴承单元机架端的节点全固定。
然后退出前处理器。
命令执行后，主窗口显示如下

８.创建不平衡力数组
继续输入
pi = acos(-1)
spin = 5000*pi/30
tinc = 0.5e-3
tend = 4
spindot = spin/tend
nbp = nint(tend/tinc) + 1
unb = 1.e-4
f0 = unb*ro_disk
上述变量是为了下面得到惯性力的分量而设置的，其含义在下面编程的时候会用到。
继续输入
*dim,spinTab,table,nbp,,,TIME
*dim,rotTab, table,nbp,,,TIME
*dim,fxTab, table,nbp,,,TIME
*dim,fyTab, table,nbp,,,TIME
*vfill,spinTab(1,0),ramp,0,tinc
*vfill,rotTab(1,0), ramp,0,tinc
*vfill,fxTab(1,0), ramp,0,tinc
*vfill,fyTab(1,0), ramp,0,tinc
这里首先定义了４个数组，然后为这四个数组填充了数据。
这四个数组：
spinTab－启动过程中各个时刻的转速
rotTab－　启动过程中各个时刻的转角
fxTab－启动过程中各个时刻的不平衡惯性力的Ｘ分量
fyTab－启动过程中各个时刻的不平衡惯性力的Ｙ分量
继续输入
tt = 0
*do,iloop,1,nbp
spinVal = spindot*tt
spinTab(iloop,1) = spinVal
spin2 = spinVal**2
rotVal = spindot*tt**2/2
rotTab(iloop,1) = rotVal
sinr = sin(rotVal)
cosr = cos(rotVal)
fxTab(iloop,1)= f0*(-spin2*sinr + spindot*cosr)
fyTab(iloop,1)= f0*( spin2*cosr + spindot*sinr)
tt = tt + tinc
*enddo
fini
一共有１４行。其算法思想请参见博文＜基于ANSYS经典界面的偏心圆盘启动的转子动力学分析－１＞，下面稍作解释。
第１行：为下面的循环设置时间的初始值。
第２行，倒数第２行：循环的开始和结束。
第３行，计算此时的角速度。
第４行，将该时刻的角速度存储到数组spinTab中。
第５行，计算角速度的平方。
第６行，计算此时的转角。
第７行，将该时刻的转角度存储到数组rotTab中。
第８－９行，计算转角的正弦和余弦。
第１０－１１行，这是该段程序的终极目的，计算不平衡惯性力的两个分量，并存储到数组fxTab，fyTab中。
第１２行，时间变量增加一个固定值，时间往前递进。相当于是转子往前又转了一个角度。
二 仿真
１.定义分析类型和参数
继续输入
/solu
antype,transient
time,tend
deltim,tinc,tinc/10,tinc*10
kbc,0
outres,all,all
其中，
第１行：进入求解器
第２行：设置要做瞬态动力学分析
第３行：设定仿真的结束时间
第４行：设定仿真的载荷子步数目
第５行：设定是斜坡激励。就是说，在前后两个载荷子步之间不平衡力是线性增加的。
第６行：请ＡＮＹＳ输出所有结果。
２.打开coriolis开关
继续输入
coriolis,on,,,on
omega,,,spin
其中，
第１行：指定考虑科里奥利效应，这是转子动力学的核心命令。
第２行：指定转子的转速，这里指定的是最终的转速。
３.施加不平衡力
继续输入
f,5,fx,%fxTab%
f,5,fy,%fyTab%
这是本分析最重要的地方。
它说明，在５号节点上施加两个方向的力，这两个力取自于前面所填充的数组fxTab，fyTab。
在计算过程中，转子每往前转一个载荷步，ＡＮＳＹＳ就会调用该数组中对应的载荷进行加载。
４.求解
继续输入
solve
fini
其含义是，开始计算，算完后退出求解器。
由于涉及到８０００个子步的计算，时间会较长，请读者耐心等待。
三 后处理
１.得到５号节点的位移
计算完毕后，继续输入
/post26
nsol,2,5,U,X,UXdisk
prod,3,2,2
nsol,4,5,U,Y,UYdisk
prod,5,4,4
add,6,3,5
sqrt,7,6,,,Ampl_At_Disk
其中，
第１行：进入到时间历程后处理器
第２行：取出５号节点的Ｘ位移，定义为２号变量，名称是UXdisk
第１行：将２号变量平方，存为３号变量。
第３行：取出５号节点的Ｙ位移，定义为４号变量，名称是Uｙdisk
第４行：将４号变量平方，存为５号变量。
第５行：把３号，５号变量相加，存为６号变量。
第６行：６号变量开方得到７号变量。
实际上，上述程序做的事情相当简单。它使用下列公式

计算节点位移。
２．绘制５号节点的位移
继续输入
/axlab,y,Displacement (m)
plvar,7
第一行设置了下面绘图的Ｙ坐标的名字，并绘制了７号变量随时间变化的曲线。

从该曲线可以看出，大致在２.７秒左右，该点的位移出现极值。此处所对应的转速应该是轴的临界转速。
３.得到５号节点的弯曲应力
继续输入
esol,8,4,5,smisc,32,Sy_At_Disk
esol,9,4,5,smisc,34,Sz_At_Disk
其含义是，取出５号节点的Ｙ，Ｚ两个方向的弯曲应力，并分别存在变量８,９中。
４.绘制５号节点的弯曲应力
继续输入
/axlab,y,Bending Stresses (N/m2)
plvar,8,9
其含义是，设置新出现的图形的Ｙ坐标的标签是Bending Stresses (N/m2)，然后绘制５号节点Ｙ，Ｚ两个方向的弯曲应力随时间变化的过程。
命令执行后，结果如下

有两条曲线，这两条曲线的变化规律与该点的位移保持一致。在２.７秒左右，应力出现最大值，其中Ｙ方向应力的最大值达到５ＭＰＡ左右。这远小于一般钢材的屈服极限，基本上可以断言，该轴不存在强度方面的问题。