机械振动--盘轴扭振系统固有频率和主振型的计算

机械振动大作业

（盘轴扭振系统固有频率和主振型的计算）

学院：航空航天工程学部

班级：04040203班

姓名：李根

学号：2010040402093

2013年5月12号

盘轴扭振系统固有频率和主振型的计算

一：简化简化分析

分析该系统为非约束性盘轴扭振系统，并简化分析分析：

：1．忽略轴的质量;

2．轴的刚度对盘的影响不做考虑;

3．将圆盘的质量集中于圆盘中心，不考虑圆盘厚度对系统的影响;4．系统为线弹性系统，盘为刚体。

对于非约束系统，其只存在刚度矩阵，不存在柔度矩阵，即不能对刚度矩阵求逆。

二：条件

圆盘：

1.几何尺寸：直径10.4d m =，厚度0.02h m =；

2.材料：杨氏模量112210(/)E N m =×，剪切模量1027.6910(/)

G N m =×密度37800(/)

kg m ρ=轴：

1.几何尺寸：直径20.04d m =，0.1a m

=2.材料：杨氏模量112210(/)E N m =×，剪切模量1027.6910(/)

G N m =×密度37800(/)

kg m ρ=三（1）：矩阵迭代法

1.1.概

概述（1）：系统主振型方程为{}[]{}21A M K A ω???=??，引入动力矩阵[][]1

D M K ???=??。任

取一个经过归一化的假设阵型{}0A ,用动力矩阵[]D 前乘它，并对通过乘法运算新得到的阵型矢量进行归一化，则得：{}110[]{}D A a A =，式中1a 为新振型矢量归一化后的系数。

（2）若{}10{}A A ≠，从1{}A 开始，重复上述步骤得：{}121[]{}D A a A =，式中2a 为新振型矢量归一化后的系数。

（3）：若{}21{}A A ≠，继续重复上述步骤，进过K 次矩阵乘法运算后，得到

{}1[]{}k k k D A a A ?=，在规定的有效位数内，{}1{}k k A A ?=时停止运算，此时的{}1k A ?即

为系统第一阶主振型(1){}A 的近似值，即：{}(1)1{}k A A ?≈，而这时的系数k a 即是系统第一阶固有频率平方倒数的近似值，即：211/k a ω≈。

该方法的精确度不依赖于假设阵型，假设阵型的好坏只影响迭代的次数。即使假设的固有频率域一阶主振型相差很远，经过充分的迭代运算，仍可求得足够精确的基频值。

求得第一阶主振型以后，利用主振型的正交性来清除掉假设阵型中的分量，然后再进行迭代求解可以是结果收敛于第二阶主振型。同理，如果我们在假设阵型中清除掉所有前s 阶这阵型分量，那么迭代的结果将得到第s+1阶固有频率及主振型。

引进清型矩阵：[]()j 1

{}{}[][]j T s

j j

A A M Q I M ==?∑（）。由于实际计算中舍入误差的存在，每次迭

代后，所得的主振型中还包含前面几阶的主振型分量，因此每次计算前都要进行清型才能保证最后收敛的主振型。

2.2.计算程序

计算程序clc clear n=8;

d1=0.4;%圆盘直径d1=0.4d2=0.04;%轴直径d2=0.04a=0.1;%轴几何尺寸

den=7800;%密度（轴和圆盘）G=7.69e+10;%剪切模量h=0.02;%圆盘厚度

J1=0.5*pi*den*h*(d1/2)^4;%转动惯量

Ip=pi*d2^4/32;%质心

k=G*Ip/a;%刚度

%假设阵型中刚度

k1=k;

k2=k;

k3=k;

k4=k;

k5=k;

k6=k;

k7=k/5;

k8=k;

J2=J1;

J3=J1;

J4=J1;

J5=J1;

J6=J1;

J7=J1;

J8=J1;

J=diag([J1,J2,J3,J4,J5,J6,J7,J8])

%刚度矩阵

K=[k2,-k2,0,0,0,0,0,0;

-k2,k2+k3,-k3,0,0,0,0,0;

0,-k3,k3+k4,-k4,0,0,0,0;

0,0,-k4,k4+k5,-k5,0,0,0;

0,0,0,-k5,k5+k6,-k6,0,0;

0,0,0,0,-k6,k6+k7,-k7,0;

0,0,0,0,0,-k7,k7+k8,-k8;

0,0,0,0,0,0,-k8,k8;]

D=inv(J)*K;%动力矩阵D

%a为新振型矢量归一化后的系数。

%ai是系统第一阶固有频率平方倒数的近似值

%在假设阵型中清除掉所有前s阶这阵型分量，

%那么迭代的结果将得到第s+1阶固有频率及主振型for j=1:n

A=[01100101]';

if j>1

D=D*(eye(n)-(B*B'*J)/(B'*J*B));

end

for i=1:1000

B=D*A;

a=B(n);

B=B/a;

if max(abs(B-A))<0.000001

break end A=B;end om1=sqrt(a)i;NM(:,j)=B;om(:,j)=om1;

x=[00.10.20.30.40.50.61.11.21.3]';y=[0;B;0];%%绘制图形

plot(x,y,'k-','linewidth',1.3,'Markersize',15)hold on end grid on NM om

3.3.结论及图像

结论及图像固有频率（单位rad/s ）：

ω

=1000*[1.3573，1.2204，1.0297，0.9929，0.7204，0.4139，0.2106，0+0.0000i ]各阶主振型：

????????????

?

???????

?????? 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9395 0.9100 0.6524 0.0529- 1.0000- 1.1508- 2.0213- 2.7374- 0.3294 0.0504 2.2192- 5.0389- 1.0000- 0.4713 13.4060 29.7304 0.1271 0.1261- 2.0222- 0.7306- 1.0000 0.2179- 24.0115- 74.8917-0.1533- 0.2912- 1.1224- 4.3469 1.0000 0.4384- 11.1157 100.38930.4820- 0.4301- 0.1675 4.8476 1.0000- 0.2841 12.6595 99.5286-0.7870- 0.5303- 1.3992 0.2442 1.0000- 0.3956 24.0443- 72.5356 0.9737- 0.5827- 2.1446 4.6163- 1.0000 0.3437- 11.8958 26.4976-= NM

1矩阵迭代法主振型

三（2）：传递矩阵法

1.概述

对于盘轴扭振系统，可以将该整体系统分割成结构形式相同的单元，即一段轴连接一个圆盘，将盘视为刚性质量元件、轴段视为柔性无质量元件，每个单元的扭转运动取决于轴段梁福安的转角和扭矩，只有这两个主要因素，因此取各元件的转角和扭矩作为状态变量，表示为[][],T

X T θ=，对整个盘轴扭振系统可以理解为状态变量是通过各个单元从系统的左侧传递到系统的右侧，我们通过对一个单元进行分析，得出各个单元的传递矩阵，进而可以求出整个系统的传递矩阵，然后再根据边界条件就可以求得关于2ω的方程，由于方程过于复杂，我们可以将该表达式的图形作出来，找出使该表达式为0的点，该点即是我们要求的ω值，从而求出该系统的各阶固有频率。

取单元中的第i 个单元来来进行分析：（1）：对刚性质量元件盘进行分析可得：

2

R L

i i R

L

i i i i

T T J θθωθ=?=?，写成矩阵的形式为：

1

11/01L

R

P i i i i k S T T θθ???????==??????

??????,该矩阵称为点传递矩阵。（2）：对柔性无质量元件轴进行分析可得：

11

1

()L R R

i i i L

R

i i k T T T

θθ????==，写成矩阵的形式为：

21

01R

L

F i i i i

S T J T θθω??????==??????

???????，该矩阵称为场传递矩阵。于是得到第i 单元两端状态的传递矩阵（从第i-1个盘右侧到第i 个盘右侧）

22

1

11/(1)R R

P F i i i i i i

k S S S J T T J k θθωω??????

????

===???????????????

?，该矩阵为单元传递矩阵。对于有n 个盘的轴系，将各个单元的单元传递矩阵反顺序连乘就可得到整个系统两端状态的矩阵形式：1100R R R i i i X S S S X SX ?==?,

即：22011122202122

()()()()R R n R R n S S T T S S θθωωωω????

??=????????????再根据轴系两边的边界条件，在该系统中00L R n T T ==，可以得到关于2ω的方程，对该关系式作图求出与x 轴的交点就可以得到系统的固有频率。2.2.计算程序

计算程序clc clear n=8;

d1=0.4;%圆盘直径d1=0.4d2=0.04;%轴直径d2=0.04a=0.1;%轴几何尺寸

den=7800;%密度（轴和圆盘）G=7.69e+10;%剪切模量h=0.02;%圆盘厚度

J1=0.5*pi*den*h*(d1/2)^4;%转动惯量Ip=pi*d2^4/32;%质心k=G*Ip/a;%刚度

%刚度

k1=k;

k2=k;

k3=k;

k4=k;

k5=k;

k6=k;

k7=k/5;

k8=k;

J2=J1;

J3=J1;

J4=J1;

J5=J1;

J6=J1;

J7=J1;

J8=J1;

%质量矩阵

M=diag([J1,J2,J3,J4,J5,J6,J7,J8])

%刚度矩阵

K=[k2,-k2,0,0,0,0,0,0;

-k2,k2+k3,-k3,0,0,0,0,0;

0,-k3,k3+k4,-k4,0,0,0,0;

0,0,-k4,k4+k5,-k5,0,0,0;

0,0,0,-k5,k5+k6,-k6,0,0;

0,0,0,0,-k6,k6+k7,-k7,0;

0,0,0,0,0,-k7,k7+k8,-k8;

0,0,0,0,0,0,-k8,k8];

D=inv(M)*K;%动力矩阵D

[A lam]=eig(D);

for i=1:n

omt(i)=sqrt(abs(lam(i,i)));

end

omt

A

%Transfer matrix method for natural frequencies mn=1;

t21=0;

step=0.1;

J=[J1,J2,J3,J4,J5,J6,J7,J8];

k=[k1k2k3k4k5k6k7k8];

for i=1:100000

om=(i-1)*step;

Tp=eye(2);

for j=1:n

T=[11/k(j);-J(j)*om^21-J(j)*om^2/k(j)]*Tp;

Tp=T;

end

if T(2,1)==0

omega(mn)=om;

mn=mn+1;

end

if T(2,1)*t21<0

omega(mn)=om-step*abs(T(2,1))/(abs(t21)+abs(T(2,1)));

mn=mn+1;

end

t21=T(2,1);

x(i)=om;

y(i)=T(2,1);

z(i)=0;

if mn>n

break

end

end

omega

for i=1:n

A=[10]';

om=omega(i);

for j=1:n

B=[11/k(j);-J(j)*om^21-J(j)*om^2/k(j)]*A;

MM(j,i)=B(1);

A=B;

end

end

MM

MMF=[00000000;MM];

x=[00.10.20.30.40.50.61.11.2];

z=[000000000];

%绘制图形

%作图求出与x轴的交点就可以得到系统的固有频率

plot(x,MMF,'k','linewidth',2,'Markersize',10)

hold on

plot(x,z)

set(gca,'YTick',[-4:0.2:4])

set(gca,'XTick',[0:0.1:1.2])

axis([0,1.4,-4,4])

grid on

3.结果表述

固有频率（单位（单位rad/s

rad/s ）：ω=1000*[0，0.2106，0.4139，0.7204，0.9929，1.0297，1.2204，1.3573]

各阶主振型：

????????????

????????

??????0.0376- 0.0841 2.9095- 1.0000 0.2166- 0.4663 1.7161- 1.00000.1033 0.1699- 3.3484 1.0000- 0.0115 0.3042 1.5616- 1.0000 1.1220- 1.1270 1.3711- 1.0000- 1.0915 1.0348- 0.0864- 1.0000 2.8263 2.0185- 0.6340 1.0000 0.1583 0.9429- 0.2164 1.0000 3.7886- 0.9344 1.2755 1.0000 0.9416- 0.5234- 0.4997 1.0000 3.7561 1.0642 0.8264- 1.0000- 1.0501- 0.0781 0.7381 1.0000 2.7374- 2.0213- 1.1508- 1.0000- 0.0529- 0.6524 0.9100 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000= MM 00.10.20.30.40.50.60.70.80.91 1.1 1.2 1.3 1.4

-4

-3.8-3.6-3.4-3.2-3-2.8-2.6-2.4-2.2-2-1.8-1.6-1.4-1.2-1-0.8-0.6-0.4-0.200.20.40.60.811.21.41.61.822.22.42.62.833.23.43.63.8

4

2传递矩阵法主振型