RotDyn

Baffalo求解

发布时间 : 2025-03-04 18:46

DotDyn

DotDyn

介绍

RotDyn类主要用于处理单个轴的转子动力学问题，通常需要借助Commonshaft类一块使用。

原理

转动动力学方程$^{[1]}$

考虑转动效应的动力方程如下：

$$
[M]{\ddot{u}}+([G]+[C]){\dot{u}}+([K]-[K_{\mathfrak{c}}]){u}={F}
$$

式中：$[M],[C],[K]$一一分别为结构质量、阻尼和刚度矩阵；

$\left[K_{\mathrm{c}}\right]$一结构转动引起的旋转软化矩阵，它使旋转坐标系下的结构刚度降低；$[G]$一结构转动对“阻尼”矩阵的贡献，在旋转坐标系下通常称为科氏矩阵或科氏阻尼矩阵，而在固定坐标系下称为回转矩阵或回转阻尼矩阵；$[F]$一—固定坐标系下的外荷载向量。在旋转坐标系下，为外荷载向量与科氏力向量之和，科氏力向量$\langle F_{\mathrm{c}}\rangle=[G]\langle\dot{u}\rangle_{0},,\langle\dot{u}\rangle_{0}$为命令IC定义的节点速度向量。

从上式中，不难看出影响一个结构的响应主要影响在于质量、阻尼、刚度，如何获得准确的结构的相关信息，是得到准确计算的关键。

在旋转坐标系下，转动效应称为科氏效应，即生成科氏矩阵和科氏力。科氏矩阵通过转速和单元形函数求得各单元的科氏矩阵，然后累加形成结构的科氏矩阵。科氏力则为科氏矩阵与节点速度之积。

在固定坐标系下，转动效应称为回转效应，即生成回转矩阵，它通过单元回转动能和单元形函数求得。在转子动力学中，将回转效应分为两部分，即与转速相关的回转矩阵和与位移相关的回转矩阵，只是动力学方程表达形式不同而已。

在ANSYS中，无论考虑回转效应或是考虑科氏效应均采用命令CORIOLIS定义。

质量参数$^{[2]}$

首先结构的质量参数，应该在网格构建完成后便知晓，下文介绍了如何通过网格来计算结构的质量，转动惯量等参数。

设有一个旋转对称体，Ox 轴为其旋转对称轴。图１为通过旋转轴 Ox 的半平面截取此物体得到的截面轮廓图。我们在轮廓线上取 N 个点, 使得连接这些点形成的封闭折线，尽可能地逼近轮廓线。从这些点中任选一点为起点，按顺时针方向绕轮廓线顺次把这些点编号为 1、2，直到 N。

折线中的每一线段，以Ox为轴旋转一周，构成一个截头圆锥台。我们先分别计算每一圆锥台的惯性参量；然后，再按一定的规则，把各个圆锥台的惯性参数综合在一起，就得到整个旋转体的惯性参量。

一个圆锥台是计算惯性参数中的基本单元，设其母线两端点的坐标为x1,y1和x2,y2。由理论力学中知识，不难得到，该圆锥台的惯性参量为：

质量:
$$
M=\frac{1}{3}\pi\rho Ly_1^2\frac{\alpha^3-1}{\alpha-1}
$$
质量中心:
$$
x_c=x_2-\frac{L}{4}\frac{\alpha^2+2\alpha+3}{\alpha^2+\alpha+1}
$$
绕x轴的转动惯量:
$$
J_x=\frac{1}{10}\pi\rho Ly_1^4\frac{\alpha^5-1}{\alpha-1}
$$
绕中心直径的转动惯量:
$$
J_{xc}=\pi\rho(L+H)y_2^2[\frac{y_2^2}{20}+\frac{(L+H)^2}{80}+\frac{(\frac{L+H}{4}-h)^2}{3}]-\pi\rho Hy_1^2[\frac{y_1^2}{20}+\frac{H^2}{80}+\frac{(L+\frac{H}{4}-h)^2}{3}]
$$
式中,

$\rho$ — 密度,

$\alpha$ — $y_2/y_1$,

$H$ — $Ly_1/(y_1+y_2)$,

$h$ — $\frac{L}{4}\frac{\alpha^2+2\alpha+3}{\alpha^2+\alpha+1}$,

圆柱是一个特殊的圆锥台，有y1=y2，它相应的惯性参量为

质量:
$$
M=\pi\rho Ly_1^2
$$
质量中心:
$$
x_c=0.5(x_1+x_2)
$$
绕x轴的转动惯量:
$$
J_x=\frac{1}{2}\pi\rho Ly_1^4
$$
绕中心直径的转动惯量:
$$
J_{xc}=0.25\pi\rho Ly_1^2(y_1^2+\frac{L^2}{3})
$$
对于编程，首先，我们可以得到任意形状的截面网格，从而通过上述方法得到结构的转动惯量参数。

整体的质量和极转动惯量取各个圆锥台的相应量的代数和；
将各圆锥台质心，平行力系合成规则得到整体的质心；
用移轴定理求出各个圆锥台对过公共质心的直径转动惯量，然后代数相加得整体的直径转动惯量。

以一个半径50mm的钢球举例($\rho =7850kg/m^3$) , 其理论计算的公式如下 :

质量:
$$
M=\frac{4}{3}\pi\rho R^3=4.11kg
$$
绕x轴的转动惯量:
$$
J_x=\frac{2}{5}MR^2=4.11kg·m^2
$$
绕y轴的转动惯量:
$$
J_y=\frac{2}{5}MR^2=4.11kg·m^2
$$

理论计算和程序计算对比结果如下：

Items	Theory	Program	Error
Mass	4.11	4.06	1.02%
$J_x$	4.11	4.06	1.02%
$J_y$	4.11	4.06	1.02%

有一个半径50 mm ，长度200mm的轴体 , 其理论计算公式如下：

质量:
$$
M=\pi\rho R^2L=12.3kg
$$
绕x轴的转动惯量:
$$
J_x=\frac{1}{2}MR^2=15.374t·mm^2
$$
绕y轴的转动惯量:
$$
J_y=\frac{1}{4}MR^2+\frac{1}{12}ML^2=48.685kg·m^2
$$

理论计算和程序计算如下，两者基本一致。

Items	Theory	Program	Error
Mass	12.3	12.3	0%
J_x	15.374	15.374	0%
J_y	48.685	48.685	0%

阻尼参数

在ANSYS中可使用7种阻尼，即Rayleigh阻尼、材料阻尼、常阻尼比、模态阻尼、单元阻尼、材料常阻尼系数、材料结构阻尼系数。在ANSYS结构动力分析中，仅在VT法谐响应分析中可采用滞变阻尼模型，库仑阻尼模型可用于瞬态动力分析，其他均采用黏性阻尼模型。在ANSYS中可同时指定多种阻尼，程序自动形成阻尼矩阵，但不是所有阻尼均可用于各种动力分析。各种分析类型可用阻尼见下表。从表中可知，动力分析类型不同，其所能够定义的阻尼类型也不同，而不同的阻尼类型在ANSYS中的输人命令亦不相同，因此在定义阻尼时应引起足够的重视。

注：1.表中“/”表示不需要，“”不考虑，“$\sqrt{}$”可考虑。2.$①$仅可用$\beta$阻尼，$,\alpha$阻尼不可用；$②$包括超单元阻尼矩阵；$③$“MP，DAMP”定义等效材料阻尼比，用于后续的谱分析或模态叠加法分析；$④$在模态叠加法中，必须采用QRDAMP法提取模态才可考虑$\alpha$阻尼$\mathcal{S}$阻尼和材料阻尼（MP，MADP）；$⑤$在谐响应分析中，命令DMPRAT和“MP，DMPR”定义的是结构阻尼比，而非模态阻尼比；$⑥$必须用命令“TB，ELAS”定义材料性质。

在 Baffalo RotDyn中设置params.Damping为常阻尼比。

常阻尼比（ConstantDampingRatio）是定义结构阻尼的最简单方法，也就是通常简称的阻尼比$\xi$一一实际阻尼与临界阻尼之比，采用命令DMPRAT定义。常阻尼比仅适用于谐响应分析、谱分析和模态叠加法瞬态动力分析。
对于材料阻尼DAMP、材料常阻尼系数DMPR和材料结构阻尼系数SDAMP可在对应的材料里设置。

材料常阻尼系数（ConstantMaterialDampingCoefficient）仅适用于谐响应分析的完全法和QRDAMP模态叠加法，采用命令“MP，DMPR”定义。

材料结构阻尼系数是频变阻尼，采用命令“TB，SDAMP”和TBFIELD定义，仅用于变分法谐响应分析，并可选择黏性阻尼模型或滞变阻尼模型。
针对单元阻尼，可以在Support中设置，RotDyn会将其转换为combin14单元和combin214单元来计算。

Support’% Spring connected to ground

% Node number,kx,K11,K22,K12,K21,Cx,C11,C22,C12,C21

刚度参数

部件的刚度比如齿轮、轴承、联轴器均有相关的规范和计算方法，可以参考其他部件Baffalo说明。

频响函数

频响函数 ( Frequency Response Function, FRF) , 它是结构的输出响应和输入激励力之比，频响函数 FRF 具有以下性质:

频响函数定义为输入位置单位激励力引起的输出位置的响应。
频响函数是系统的固有特性, 与系统本身有关, 与激励、响应等外界因素没有关系。
频响函数具有互易性, 即 Hij = Hji, 也就是说, j 点单位激励力在 i 点引起的响应等于 i 点单位激励力在 j 点引起的响应, 这也表明频响函数矩阵是对称的。
频响函数是复值函数, 因而其可以用幅值与相位或者实部与虚部表示, 故而频响函数具有幅频、相频和实频、虚频等多种表现形式。”

在RotDyn中设置params.Type=1可设置频响函数分析。

单自由度频响函数

以单自由度的频响函数为例，

“在低频段, FRF 的幅值是 1 / k( k >> ω2 m + jωc) , 相位为 0, 表明共振频率以下的频率段主要用占主导地位的刚度项来描述。
在高频段, FRF 的幅值为 - 1 / ω2 m( ω2 m >> jωc + k) , 相位为 - 180°, 表明共振频率以上的频率段主要用占主导地位的质量项来描述。
理论上, 无阻尼固有频率处的 FRF 幅值应是无穷大, 但是由于阻尼的存在, 导致共振频率处的幅值不会无穷大, 其幅值为 1 / ωc, 相位突变 180°, 表明在共振频率处主要受阻尼控制。

AMRotor

RotDyn的时序分析和PID控制，主要基于开源转子动力学分析软件AMrotor

类结构

输入 input：

BalanceQuality : 动平衡质量
UnBalanceForce : 不平衡量me
KeyNode : 关键点
BCNode : 边界点
PointMass : 附加的节点质量
Discs : 附加转盘
Speed : 计算转速
MaterialNum : 材料编号
Shaft : 轴
PIDController : PID控制器
TimeSeries : 时序
LUTBearing : 刚度随转速变化轴承
Table : 轴承刚度和阻尼map
TorBearing: 转动方向的轴承
Bearing : 轴承包含刚度和阻尼参数
SpeedRange : 升速分析的转速范围

参数 params:

ez ：z向偏心
ey : y向偏心
Type : 求解类型 Type=1 FRF分析；Type=2 模态计算；Type=3为谐响应分析；Type=4 定转速时序分析；Type=5 s升速或降速时序分析
Rayleigh : 瑞利阻尼
PrintMode ：输出模态图
PrintCampbell : 输出坎贝尔图
Position : 轴位置
PStress : 考虑预应力
Name : 名称
NStep ：计算步数
NMode : 模态数
Modopt : 模态计算方法
HRopt: 谐响应计算方法
Material : 材料
Freq : 频率计算范围
Damping : 阻尼比
Coriols : 科式力

输出 output :

Xc : 轴重心位置
TimeSeriesResult : 时序分析结果
Time : 时间序列
Ewb : 负进动
Ewf : 正进动
eigenValue : 模态分析幅值
eigenVector ：模态分析方向
ModeResult : 模态分析结果
FRFResult ：FRF分析结果
Mass : 轴质量
SpeedUp : 加速曲线
TotalElement : 总的单元数
TotalNode ：总的节点数
CriticalSpeed : 临界转速
Shape : 轴形状
Campbell : 坎贝尔图数据
Assembly : 轴装配

案例

下图为一个等截面转子，长 l 为 1m，圆形截面的直径 d=0.04m，两端为简支支承。转子材料的质量密度为 7800kg/m 3 ，弹性模量 E 为 2.0×10 11 N/m 2 。求此转子在静止状态下的固有频率。

此转子的各阶固有圆频率 $\omega i$ 有理论解如下：

$$
\omega_i=i^2\pi^2\sqrt{\frac{EJ}{\gamma Al^4}}
$$
以截面惯性矩 $J＝πd^4 /64$，截面积 $A＝πd^2 /4$，固有频率 $f＝ \omega/2π$ 代入，得到转子的各阶固有频率 $f_i$为
$$
f_i=\frac{i^2\pi}{8}\frac{d}{l^2}\sqrt{\frac{E}{\gamma}}
$$
代入具体数据，得到前 3 阶固有频率为 79.5403、318.1612、715.8627Hz，相应的振型为半个周期、一个周期和一个半周期的正弦曲线。

% Shaft 1
inputshaft1.Length = 1000;
inputshaft1.ID = [0,0];
inputshaft1.OD = [40,40];
paramsshaft1 = struct();
obj1 = shaft.Commonshaft(paramsshaft1, inputshaft1);
obj1 = obj1.solve();

mat{1,1}=obj1.params.Material;
inputRotDyn.Shaft=obj1.output.BeamMesh;
inputRotDyn.Speed=[0,15000,30000];
inputRotDyn.MaterialNum=1;
inputRotDyn.BCNode=[1,1,1,1,0,0,0;size(obj1.output.Node,1),0,1,1,0,0,0];
paramsRotDyn.Material=mat;
paramsRotDyn.PrintCampbell=1;
paramsRotDyn.PrintMode=1;
paramsRotDyn.Freq=[0,3000];
obj2 = solve.RotDyn(paramsRotDyn,inputRotDyn);
Plot(obj2);
obj2 = obj2.solve();
ANSYSSolve(obj2.output.Assembly)
    
obj2=PlotCampbell(obj2,'NMode',8);
for i=1:10
  PlotMode(obj2,2,i,'scale',10)
end

Add Discs (Flag=2)

如图所示的刚性支承单圆盘转子，圆盘的质量$m=20kg$，半径$R=120mm$,转轴的跨度$l=750mm$，直径$d=30mm$。圆盘到左支点的距离$a=l/3=250mm$。

仅考虑轴的弯曲但不计轴的质量，考虑回转效应时的频率方程为：

$$
\omega^{4}-2\Omega\omega^{3}-2~~1340661 \times 10^{6}\omega^{2}+1~~7674781\times10^{6}\Omega\omega+1~2052387 \times 10^{11}=0
$$
其中，$\varOmega$为转速（即激励），$\omega$为待求涡动频率，两者单位均为rad/s。

ID=30;
inputshaft1.Length =750;
inputshaft1.ID = [0,0];
inputshaft1.OD = [ID,ID];
paramsshaft1.Beam_N = 16;
obj1 = shaft.Commonshaft(paramsshaft1, inputshaft1);
obj1 = obj1.solve();

mat{1,1}=obj1.params.Material;
inputRotDyn.Shaft=obj1.output.BeamMesh;
inputRotDyn.Speed=[0,400,800,1200,1600,2000];% Unit: RPM
inputRotDyn.MaterialNum=1;
inputRotDyn.BCNode=[1,1,1,1,0,0,0;size(obj1.output.Node,1),1,1,1,0,0,0];
paramsRotDyn.Material=mat;
paramsRotDyn.PrintCampbell=1;
paramsRotDyn.Position=[0,0,0,0,0,0];
Dyn = solve.RotDyn(paramsRotDyn,inputRotDyn);
[Dyn,Num1]= AddCnode(Dyn,250);
DiscMass=20*10^-3;
rou=obj1.params.Material.Dens;
OD=240;
L=DiscMass/((OD^2-ID^2)*pi/4*rou);
Dyn= AddDisc(Dyn,Num1,OD,ID,L,1);
Plot(Dyn);
Dyn = Dyn.solve();
ANSYSSolve(Dyn.output.Assembly)
PlotCampbell(Dyn,'NMode',4);

Add Userdefined Disc, output moment of inertia (Flag=3)

OD=30;
inputshaft1.Length =750;
inputshaft1.ID = [0,0];
inputshaft1.OD = [OD,OD];
paramsshaft1.Beam_N = 30;
obj1 = shaft.Commonshaft(paramsshaft1, inputshaft1);
obj1 = obj1.solve();
mat{1,1}=obj1.params.Material;
inputRotDyn.Shaft=obj1.output.BeamMesh;
inputRotDyn.Speed=[0,1500,3000];% Unit: RPM
inputRotDyn.MaterialNum=1;
inputRotDyn.BCNode=[1,1,1,1,0,0,0;size(obj1.output.Node,1),1,1,1,0,0,0];
paramsRotDyn.Material=mat;
paramsRotDyn.PrintCampbell=1;
paramsRotDyn.Position=[0,0,0,0,0,0];
Dyn = solve.RotDyn(paramsRotDyn,inputRotDyn);
[Dyn,Num1]= AddCnode(Dyn,250);
[Dyn,Num2]= AddCnode(Dyn,500);
a1=Point2D('Circle center');
a1=AddPoint(a1,0,0);
a1=AddPoint(a1,[-50;50],[0;0]);
b1=Line2D('Semi circle');
b1=AddCircle(b1,50,a1,1,'ang',180);
b1=AddLine(b1,a1,2);
S1=Surface2D(b1);
a2=Point2D('Square corner');
a2=AddPoint(a2,[-100;-100;100;100],[0;50;50;0]);
b2=Line2D('Square');
b2=AddCurve(b2,a2,1);
S2=Surface2D(b2);
Dyn= AddUserdefinedDisc(Dyn,Num1,S1,1);
Dyn= AddUserdefinedDisc(Dyn,Num2,S2,1);
Plot(Dyn);
disp(Dyn.input.PointMass)

加入自定义的截面并计算其质量、转动惯量和弯曲惯量。

Successfully revolve to solid mesh .
35.0000 0.0041 4.0578 4.0580
69.0000 0.0123 15.3742 48.6848

Add Blade (Flag=4)

OD=425;
inputshaft1.Length =1000;
inputshaft1.ID = [0,0];
inputshaft1.OD = [OD,OD];
paramsshaft1.Beam_N = 30;
obj1 = shaft.Commonshaft(paramsshaft1, inputshaft1);
obj1 = obj1.solve();
mat{1,1}=obj1.params.Material;
inputRotDyn.Shaft=obj1.output.BeamMesh;
inputRotDyn.Speed=[0,1500,3000];% Unit: RPM
inputRotDyn.MaterialNum=1;
inputRotDyn.BCNode=[1,1,1,1,0,0,0;size(obj1.output.Node,1),1,1,1,0,0,0];
paramsRotDyn.Material=mat;
paramsRotDyn.NMode=15;
paramsRotDyn.Freq=[0,3000];
paramsRotDyn.PrintCampbell=1;
paramsRotDyn.PrintMode=1;
paramsRotDyn.Position=[0,0,0,0,0,0];
Dyn = solve.RotDyn(paramsRotDyn,inputRotDyn);
[Dyn,Num1]= AddCnode(Dyn,500);
a=Point2D('Circle center');
a=AddPoint(a,[-160;-160;-155;-155;-95;-55;55;90;155;155;160;160;-160],...
[212.5;350;350;380;390;595;595;385;355;320;320;212.5;212.5]);
a=AddPoint(a,[-45;-20;18;45],...
[595;1050;1050;595]);
b1=Line2D('Section1');
b1=AddCurve(b1,a,1);
S1=Surface2D(b1);
b2=Line2D('Section2');
b2=AddCurve(b2,a,2);
S2=Surface2D(b2);
Dyn= AddUserdefinedDisc(Dyn,Num1,S1,1);
Dyn= AddBlade(Dyn,Num1,S2,1,30,3);
Plot(Dyn);
Dyn = Dyn.solve();
ANSYSSolve(Dyn.output.Assembly);
PlotCampbell(Dyn,'NMode',8);
for i=1:8
PlotORB(Dyn,2,i,'scale',100)
end

增加叶片，RotDyn会换算叶片的质量的转动惯量到轴中心。

Multi Discs & Calculate critical speed (Flag=5)

m1=0.102;% unit:ton
JT1=6377*2;% unit: ton/mm2
JD1=6377;% unit: ton/mm2
% Shaft1
OD=50;% Unit: mm
inputshaft1.Length =1200;
inputshaft1.ID = [0,0];
inputshaft1.OD = [OD,OD];
paramsshaft1.Beam_N = 30;
obj1 = shaft.Commonshaft(paramsshaft1, inputshaft1);
obj1 = obj1.solve();
mat{1,1}=obj1.params.Material;
inputRotDyn.Shaft=obj1.output.BeamMesh;
inputRotDyn.Speed=[0,1000,2000,3000];% Unit: RPM
inputRotDyn.MaterialNum=1;
inputRotDyn.BCNode=[1,1,1,1,0,0,0];
paramsRotDyn.Material=mat;
paramsRotDyn.PrintCampbell=1;
Dyn = solve.RotDyn(paramsRotDyn,inputRotDyn);
[Dyn,Num1]= AddCnode(Dyn,400);
Dyn=AddPointMass(Dyn,Num1,m1,JT1,JD1);
[Dyn,Num2]= AddCnode(Dyn,800);
Bound=[1,1,1,0,0,0];
Dyn= AddBCNode(Dyn,Num2,Bound);
[Dyn,Num3]= AddCnode(Dyn,1200);
Dyn=AddPointMass(Dyn,Num3,m1,JT1,JD1);
Plot(Dyn);
Dyn = Dyn.solve();
ANSYSSolve(Dyn.output.Assembly);
PlotCampbell(Dyn,'NMode',12);
Dyn=CalculateCriticalSpeed(Dyn);
disp(Dyn.output.CriticalSpeed);

计算临界转速

1.0e+03 *

Add Bearing (Flag=6)

如下图所示钢质圆盘转子结构，材料的弹性模量为$200\mathrm{GPa}$，密度$7830kg/m^{3}$，泊松系数为0.3，长度分别为$L_1=0.2m$，$L_2=0 ~ 3m$，$L_3=0.5m$，$L_4=0.3m$，轴径$d=0.05m$，圆盘半径分别为 $R_{1}=0.12m$ ， $R_{2}=0~2m$ ， $R_{3}=0.2m$ ,厚度分别为 $t_{1}=0.05m$ ， $t_{2}=0.05m$ ， $t_{3}=0.06m$ 。轴两端为简支边界条件，图中支承仅为示意。

L1=200;
L2=300;
L3=500;
L4=300;
t1=50;t2=50;t3=60;
R1=120;
R2=200;
R3=200;
d=50;
% Shaft1
inputshaft1.Length = [L1-t1/2;L1+t1/2;L1+L2-t2/2;L1+L2+t2/2;L1+L2+L3-t3/2;L1+L2+L3+t3/2;L1+L2+L3+L4];
inputshaft1.ID = [[0,0];[0,0];[0,0];[0,0];[0,0];[0,0];[0,0]];
inputshaft1.OD = [[d,d];[R1*2,R1*2];[d,d];[R2*2,R2*2];[d,d];[R3*2,R3*2];[d,d]];
paramsshaft1.Beam_N = 16;
paramsshaft1.N_Slice=201;
obj1 = shaft.Commonshaft(paramsshaft1, inputshaft1);
obj1 = obj1.solve();
Plot3D(obj1)
% Rigid Boundary
mat{1,1}=obj1.params.Material;
inputRotDyn.Shaft=obj1.output.BeamMesh;
inputRotDyn.Speed=[0,1000,2000,3000,4000,5000,6000];
inputRotDyn.MaterialNum=1;
inputRotDyn.BCNode=[1,1,1,1,1,0,0;size(obj1.output.Node,1),0,1,1,1,0,0];
paramsRotDyn.Material=mat;
paramsRotDyn.PrintCampbell=1;
paramsRotDyn.ShaftTorsion=1;
paramsRotDyn.PrintMode=1;
Dyn1 = solve.RotDyn(paramsRotDyn,inputRotDyn);
Dyn1 = Dyn1.solve();
Plot(Dyn1);
ANSYSSolve(Dyn1.output.Assembly);
Dyn1=PlotCampbell(Dyn1,'NMode',12);
% Bearing
mat{1,1}=obj1.params.Material;
inputRotDyn.Shaft=obj1.output.BeamMesh;
inputRotDyn.Speed=[0,1000,2000,3000,4000,5000,6000];
inputRotDyn.MaterialNum=1;
inputRotDyn.Bearing=[1,1e8,8e3,3e3,0,0,0,1e-3*8e3,2e-4*3e3,0,0;size(obj1.output.Node,1),0,8e3,3e3,0,0,0,1e-3*8e3,2e-4*3e3,0,0];
paramsRotDyn.Material=mat;
paramsRotDyn.PrintCampbell=1;
paramsRotDyn.ShaftTorsion=1;
paramsRotDyn.PrintMode=1;
inputRotDyn.BCNode=[];
Dyn2 = solve.RotDyn(paramsRotDyn,inputRotDyn);
Dyn2 = Dyn2.solve();
Plot(Dyn2);
ANSYSSolve(Dyn2.output.Assembly);
Dyn2=PlotCampbell(Dyn2,'NMode',12);
disp(Dyn1.output.Campbell);
disp(Dyn2.output.Campbell);