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

即：
再根据轴系两边的边界条件，在该系统中 ，可以得到关于 的方程，对该关系式作图求出与x轴的交点就可以得到系统的固有频率。
2.
clc
clear
n=8;
d1=0.4; %圆盘直径 d1=0.4
d2=0.04; %轴直径 d2=0.04
a=0.1;%轴几何尺寸
den=7800;%密度 （轴和圆盘）
G=7.69e+10; %剪切模量
机械振动大作业
（盘轴扭振系统固有频率和主振型的计算）
学院：航空航天工程学部
班级：04040203班
姓名：李根
学号：2010040402093
2013年 5 月 12号
盘轴扭振系统
固有频率和主振型的计算
一：
该系统为非约束性盘轴扭振系统，并简化分析：
1．忽略轴的质量;
2．轴的刚度对盘的影响不做考虑;
3．将圆盘的质量集中于圆盘中心，不考虑圆盘厚度对系统的影响;
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）：
=1000*[ 0 ，0.2106 ，0.4139 ，0.7204，0.9929，1.0297，1.2204，1.3573]
引进清型矩阵： 。由于实际计算中舍入误差的存在，每次迭代后，所得的主振型中还包含前面几阶的主振型分量，因此每次计算前都要进行清型才能保证最后收敛的主振型。
2.
clc
clear
n=8;
d1=0.4; %圆盘直径 d1=0.4
d2=0.04; %轴直径 d2=0.04
a=0.1; %轴几何尺寸
den=7800; %密度 （轴和圆盘）
各阶主振型：
图2 传递矩阵法主振型
4．系统为线弹性系统，盘为刚体。
对于非约束系统，其只存在刚度矩阵，不存在柔度矩阵，即不能对刚度矩阵求逆。
二：
圆盘：
1.几何尺寸：直径 ，厚度 ；
2. 材料：杨氏模量 ，剪切模量
密度
轴：
1. 几何尺寸：直径 ，
2. 材料：杨氏模量 ，剪切模量
密度
三（1）：矩阵迭代法
1.概述
（1）：系统主振型方程为 ，引入动力矩阵 。任取一个经过归一化的假设阵型 ,用动力矩阵 前乘它，并对通过乘法运算新得到的阵型矢量进行归一化，则得： ，
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;
z(i)=0;
if mn>n
break
end
end
omega
for i=1:n
A=[1 0]';
om=omega(i);
for j=1:n
B=[1 1/k(j);-J(j)*om^2 1-J(j)*om^2/k(j)]*A;
MM(j,i)=B(1);
A=B;
end
end
MM
MMF=[0 0 0 0 0 0 0 0;MM];
break
end
A=B;
end
om1=sqrt(a)
i;
NM(:,j)=B;
om(:,j)=om1;
x=[0 0.1 0.2 0.3 0.4 0.5 0.6 1.1 1.2 1.3]';
y=[0;B;0];
%%绘制图形
plot(x,y,'k-','linewidth',1.3,'Markersize',15)
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;
x=[0 0.1 0.2 0.3 0.4 0.5 0.6 1.1 1.2 ];
z=[0 0 0 0 0 0 0 0 0 ];
%绘制图形
%作图求出与x轴的交点就可以得到系统的固有频率
plot(x,MMF,'k','linewidth',2,'Markersize',10)
hold on
plot(x,z)
0,0,0,0,-k6,k6+k7,-k7,0;
0,0,0,0,0,-k7,k7+k8,-k8;
0,Hale Waihona Puke ,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
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);
该方法的精确度不依赖于假设阵型，假设阵型的好坏只影响迭代的次数。即使假设的固有频率域一阶主振型相差很远，经过充分的迭代运算，仍可求得足够精确的基频值。
求得第一阶主振型以后，利用主振型的正交性来清除掉假设阵型中的分量，然后再进行迭代求解可以是结果收敛于第二阶主振型。同理，如果我们在假设阵型中清除掉所有前s阶这阵型分量，那么迭代的结果将得到第s+1阶固有频率及主振型。
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=[0 1 1 0 0 1 0 1]';
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
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;
取单元中的第i个单元来来进行分析：
（1）：对刚性质量元件盘进行分析可得： ，写成矩阵的形式为：
,该矩阵称为点传递矩阵。
（2）：对柔性无质量元件轴进行分析可得： ，写成矩阵的形式为：
，该矩阵称为场传递矩阵。
于是得到第i单元两端状态的传递矩阵（从第i-1个盘右侧到第i个盘右侧）
，该矩阵为单元传递矩阵。
对于有n个盘的轴系，将各个单元的单元传递矩阵反顺序连乘就可得到整个系统两端状态的矩阵形式： ,
mn=1;
t21=0;
step=0.1;
J=[J1,J2,J3,J4,J5,J6,J7,J8];
k=[k1 k2 k3 k4 k5 k6 k7 k8];
for i=1:100000
om=(i-1)*step;
Tp=eye(2);
for j=1:n
T=[1 1/k(j);-J(j)*om^2 1-J(j)*om^2/k(j)]*Tp;
式中 为新振型矢量归一化后的系数。
（2）若 ，从 开始，重复上述步骤得： ，
式中 为新振型矢量归一化后的系数。
（3）：若 ，继续重复上述步骤，进过K次矩阵乘法运算后，得到 ，在规定的有效位数内， 时停止运算，此时的 即为系统第一阶主振型 的近似值，即： ，而这时的系数 即是系统第一阶固有频率平方倒数的近似值，即： 。