求助：传递矩阵法求解固有频率和振型
%求解转子系统前三个临界转速和主振型的传递矩阵法a=1;%截面剪切系数
u=0.3;%泊松比
rou=7800;
E=2.0e11;
L1=0.187;%轴的总长
G=E/(2*(1+u));
K1=2.0e7;%弹性支撑刚度
%
K=; %弹性支撑的刚度
L=*1e-3;%轴的分段
d1=*1e-3; %%%光轴尺寸
I=pi*d1.^4/64;
A=pi*d1.^2/4;
v=6*E*I./(a*G*A.*L.^2);
M=zeros(1,length(L));
Jp=zeros(1,length(L));
Jd=zeros(1,length(L));
for jk=1:length(L)
m(jk)=rou*pi*d1(jk)^2*L(jk)/8 ;%轴的质量分布
jp(jk)=rou*A(jk)*L(jk)^2*d1(jk)^2/16; %轴的极转动惯量Jp=mi*r^2/2;
end
M(:,2:length(L)-1)=m(:,2:length(L)-1)+m(:,3:length(L));
M(1)=m(1);
M(length(L))=m(length(L));
Jp(:,2:length(L)-1)=m(:,2:length(L)-1)+m(:,3:length(L));
Jp(1)=jp(1);
Jp(length(L))=jp(length(L));
Jd=Jp./2;
J=Jp-Jd;%直径转动惯量
k=0;
for w=0:10:20000;
for i=1:length(L)
T(:,:,i)=[1+(L(i)^3)*(1-v(i))*(M(i)*w^2-K(i))/(6*E*I(i)) L(i)+L(i)^2*J(i)*w^2/(2*E*I(i)) L(i)^2/(2*E*I(i)) L(i)^3*(1-v(i))/(6*E*I(i));
(L(i)^2)*(M(i)*w^2-K(i))/(2*E*I(i)) 1+L(i)*J(i)*w^2/(E*I(i)) L(i)/(E*I(i)) L(i)^2/(2*E*I(i));
L(i)*(M(i)*w^2-K(i)) J(i)*w^2 1 L(i);
M(i)*w^2-K(i) 0 0 1]; %传递矩阵
end
H=T(:,:,1);
for i2=2:length(L);
H=T(:,:,i2)*H;
end
F=H(3,1)*H(4,2)-H(3,2)*H(4,1) ; %剩余量
if F*(-1)^k < 0 %求解临界转速
k=k+1;
wi(k)=w ; %固有圆频率
w=wi(k);
ni(k)=wi(k)*30/pi ; %临界转速
end
end
ni=ni'
wi=wi'
f=ni/60
%绘制前三阶模态振型
for i1=1:3;
w=wi(i1);
for j=1:length(L)-1;
T(:,:,j)=[1+(L(j)^3)*(1-v(j))*(M(j)*w^2-K(j))/(6*E*I(j)) L(j)+L(j)^2*J(j)*w^2/(2*E*I(j)) L(j)^2/(2*E*I(j)) L(j)^3*(1-v(j))/(6*E*I(j));
(L(j)^2)*(M(j)*w^2-K(j))/(2*E*I(j)) 1+L(j)*J(j)*w^2/(E*I(j)) L(j)/(E*I(j)) L(j)^2/(2*E*I(j));
L(j)*(M(j)*w^2-K(j)) J(j)*w^2 1 L(j);
M(j)*w^2-K(j) 0 0 1];
end
H=T(:,:,1);
for j=2:length(L);
H=T(:,:,j)*H;
end
b=-H(4,1)/H(4,2);
X(:,1)=(');
for n=2:length(L)+1;
X(:,n)=T(:,:,n-1)*X(:,n-1) ; %相邻两质点右边的传递关系
end
x=zeros(1,length(L)+1);
for j1=1:length(L);
y(j1)=X(1,j1); %挠度
z(j1)=X(3,j1); %弯矩
x(j1+1)=L(j1)+x(j1);
end
y(length(L)+1)=X(1,length(L)+1);
z(length(L)+1)=X(3,length(L)+1);
y=y/max(abs(y)) ;%归一化
z=z/max(abs(z));%归一化
subplot(3,1,i1)
plot(x,y,'b-',x,z,'r:');
Tit=['第一阶频率的振型和弯矩图';'第二阶频率的振型和弯矩图';'第三阶频率的振型和弯矩图'];
title(Tit(i1,:))
xlabel('轴长'),ylabel('不平衡值')
%axis()
grid on
z;
end
legend('振型','弯矩')
振型图怎么是这个德行啊！
程序是按照论坛里面的一个例子做的，但是我没有加圆盘，就是一根光轴，计算的固有频率和ansys计算的、有限元法计算的结果相差都很大
补充内容 (2013-7-10 21:42):
转帖！ 振型图怎么能是“不平衡度”的纵坐标呢。
我没坐过转子方面的，但是如果是弯曲，纵坐标肯定是节点挠度，或者位移加速度什么的。
如果是扭转，那么肯定是扭转角度，速度之类……
仅供参考 dw04116 发表于 2013-3-25 11:14 static/image/common/back.gif
振型图怎么能是“不平衡度”的纵坐标呢。
我没坐过转子方面的，但是如果是弯曲，纵坐标肯定是节点挠度，或 ...
状态向量的第一行就是挠度啊！ 我定义的计算，第一行就是第一个节点的广义位，第二行是第二个点的……
第一列是一阶频率…… b=-H(4,1)/H(4,2);
X(:,1)=(');
for n=2:length(L)+1;
X(:,n)=T(:,:,n-1)*X(:,n-1) ; %相邻两质点右边的传递关系
end
为什么 b=-H(4,1)/H(4,2); 你这M和JP的计算公式错了。M的个数比L的多一个。
补充内容 (2013-9-7 14:43):
Jp(:,2:length(L)-1)=m(:,2:length(L)-1)+m(:,3:length(L));计算极转动惯量的公式也错了。
页:
[1]