碰摩转子故障仿真程序
该程序在不加碰摩力时，得出结果很好。加入碰摩力，积分结果发散，得不出仿真结果。
不平衡量在m2处，碰摩也发生在m2处。
仿真程序如下，matlab程序，newmark积分。
clear all
clc
%%----------------------------------
M=[3.5, 0,0, 0, 0, 0
0, 3.5,0, 0, 0, 0
0, 0,3.5, 0, 0, 0
0, 0,0, 3.5, 0, 0
0, 0,0, 0, 3.5, 0
0, 0,0, 0, 0, 3.5];
C=[1.05e3, 0,0, 0, 0, 0
0, 1.2e3,0, 0, 0, 0
0, 0,1.05e3, 0, 0, 0
0, 0,0, 1.05e3, 0, 0
0, 0,0, 0, 1.2e3, 0
0, 0,0, 0, 0, 1.05e3];
K=[3.1543e7, -3.0173e7,1.2342e7, 0, 0, 0
-3.0173e7, 4.3885e7,-3.0173e7, 0, 0, 0
1.2342e7, -3.0173e7,3.1543e7, 0, 0, 0
0, 0,0, 3.1543e7, -3.0173e7,1.2342e7
0, 0,0, -3.0173e7, 4.3885e7,-3.0173e7
0, 0,0, 1.2342e7, -3.0173e7,3.1543e7];
%%------------------------------------
u=zeros(6,1);
v=zeros(6,1);
x(:,1)=u; %位移
x1(:,1)=v; %速度
ee=0.13e-3; %偏心距
h=0.1e-3;; %动静间隙
kk=0.5e9; %碰摩时接触刚度
miur=0.2; %摩擦系数
delta0=0.0000001; %不对中
RR0=0.1;
delt=0.5;
alpha=0.25;
dt=2*pi/(1024*8); %积分步长
b0=1/(alpha*dt^2);
b1=delt/(alpha*dt);
b2=1/(alpha*dt);
b3=(1/(2*alpha))-1;
b4=(delt/alpha)-1;
b5=(dt/2)*((delt/alpha)-2);
b6=dt*(1-delt);
b7=delt*dt;
w=0;
t_max=4;
ii=1;
t(1)=0;
f=86;
q=zeros(6,1);
w=0;
while t(ii)<t_max
w=2*pi*f; %角速度
fux(ii)=5.5*ee*(w*w*cos(w*t(ii))); %不平衡力
fuy(ii)=5.5*ee*(w*w*sin(w*t(ii))); %不平衡力
%----------------------碰摩---------------------------
rrr(ii)=sqrt((x(2,ii))^2+(x(5,ii)-delta0)^2);
ggg(ii)=(rrr(ii)-h);
fai(ii)=((x(1,ii)*x1(5,ii)-(x(5,ii)-delta0)*x1(2,ii))/rrr(ii))+RR0*w;
if fai(ii)>0
thetar=1;
elseif fai(ii)<0
thetar=-1;
else
thetar=0;
end
if (ggg(ii))>0
N(ii)=(kk*ggg(ii)/rrr(ii));
frx(ii)=(-N(ii))*((x(2,ii))-miur*((x(5,ii))-delta0));
fry(ii)=(-N(ii))*((miur*x(2,ii))+(x(5,ii)-delta0));
else
fry(ii)=0;
frx(ii)=0;
end
Q2(:,ii)=zeros(6,1);
Q2(2,ii)=Q2(2,ii)+frx(ii);
Q2(5,ii)=Q2(5,ii)+fry(ii);
%------------------------------------------------------------
Q(:,ii)=zeros(6,1);
Q1(:,ii)=zeros(6,1);
Q1(2,ii)=Q1(2,ii)+fux(ii);
Q1(5,ii)=Q1(5,ii)+fuy(ii);
Q(:,ii)=Q1(:,ii);
Ccc=C;
K1=K+b0*M+b1*Ccc ;
if ii==1
x2(:,ii)=inv(M)*(Q(:,ii)-K*x(:,ii)-Ccc*x1(:,ii));
q(:,ii)=Q(:,ii)+M*(b0*x(:,ii)+b2*x1(:,ii)+b3*x2(:,ii))+Ccc*(b1*x(:,ii)+b4*x1(:,ii)+b5*x2(:,ii));
x(:,ii+1)=K1\q(:,ii);
x2(:,ii+1)=b0*(x(:,ii+1)-x(:,ii))-b2*x1(:,ii)-b3*x2(:,ii);
x1(:,ii+1)=x1(:,ii)+b6*x2(:,ii)+b7*x2(:,ii+1);
else
q(:,ii+1)=Q(:,ii)+M*(b0*x(:,ii)+b2*x1(:,ii)+b3*x2(:,ii))+Ccc*(b1*x(:,ii)+b4*x1(:,ii)+b5*x2(:,ii));
x(:,ii+1)=K1\q(:,ii+1);
x2(:,ii+1)=b0*(x(:,ii+1)-x(:,ii))-b2*x1(:,ii)-b3*x2(:,ii);
x1(:,ii+1)=x1(:,ii)+b6*x2(:,ii)+b7*x2(:,ii+1);
end
ii=ii+1;
t(ii)=t(ii-1)+dt;
end
%%----------------------------------------
T=1/dt;
TT=T*t_max;
t=3000:TT;
subplot(2,2,(1:2))
plot(t/T,x(2:2,t)*1000)
title('位移-1x')
grid on
subplot(2,2,3)
Nx=length(x(2:2,t));
Y=fft(x(2:2,t)*T,Nx);
pyy=Y.*conj(Y)/(1024*8);
f=T*(0:Nx/2-1)/Nx;
plot(f,pyy(1:(Nx/2)))
title('幅值谱-节点1x')
grid on
subplot(224)
plot(x(2:2,t)*1000,x(5:5,t)*1000)
title('轴心轨迹')
grid on 得出结论了
上面的程序是对的,只是积分步长不够小,改小就ok
回复 #2 zhaowu 的帖子
你好，我想问一下，是不是用有限元方法列出的转子系统方程后，直接把几个系数矩阵带入就可以了？可是我建立一个双圆盘转子的运动微分方程，为什么结出的结果不正确？ 难道我求质量矩阵，刚度矩阵，阻尼矩阵出现了问题？这位大侠能否留下联系方式，请教下问题 试着积分步长取小些
qq:120852933 本帖最后由 VibInfo 于 2016-5-18 16:29 编辑
原帖由 zhaowu 于 2007-9-18 20:01 发表
试着积分步长取小些
qq:120852933 .................. 本帖最后由 VibInfo 于 2016-5-18 16:29 编辑
原帖由 qijunshuai 于 2007-9-13 10:32 发表
你好，我想问一下，是不是用有限元方法列出的转子系统方程后，直接把几个系数矩阵带入就可以了？可是我建立一个双圆盘转子的运动微分方程，为什么结出的结果不正确？ 难道我求质量矩阵，刚度矩阵，阻尼矩阵出现 ...
直接把几个系数矩阵带入方程，形成一个二阶微分方程组，然后求解就可以了 请问楼主,这个程序碰摩力加在那里? Q2(2,ii)=Q2(2,ii)+frx(ii);
Q2(5,ii)=Q2(5,ii)+fry(ii);
frx，fry表示x,y方向碰摩力 谢谢楼主正学习中好东东当然要收藏了{:{23}:} 感谢楼主好东西哈
页:
[1]