|
|
楼主 |
发表于 2010-9-27 19:40
|
显示全部楼层
回复 jgwang 的帖子
function out=liewenfun(t,z,flag)
omega=1000; % 轴的转速
a=0.0135; % 裂纹深度
r=0.015;L1=0.57; % r转轴半径/，L1转轴长度/
R=0.025;L=0.012; % R轴承半径/，L轴承长度/
b=0.05*10^(-3); % b圆盘偏心率/
m1=4.0;m2=32.1; % m1轴承处的等效集中质量/，m2转轴中央圆盘质量/
c=0.11*10^(-3); % c轴承间隙/
beta=pi;mu=0.018; % beta裂纹方向与偏心之间的夹角，mu润滑油粘度/
c1=1050.0;c2=2100.0; % c1转子在轴承处的结构阻尼/，c2转子在圆盘处的结构阻尼/
k=2.5*10^7;g=9.8; % k转轴刚度/
epsilon=3.3064e-007;delta=3.0244e+006; % epsilon与delta为仅与裂纹深度a有关的相关刚度参数
P=0.5*m2; % P转子圆盘重量的一半
theta=omega*t+beta; % theta裂纹截面转角
F=(1+cos(theta))/2; % F裂纹开闭规律，描述的是裂纹转角对实际裂纹轴刚度的影响规律(裂纹开闭规律模型有好几种，我用的是比较简单的余弦波模型)
s=(mu*omega*R*L*((R/c)^2)*(L/(2*R))^2)/P; % s为Sommerfeld修正系数，c平均油膜厚度
xi1=c1/(m1*omega);eta1=k/(m1*omega^2); % 方程无量纲化后的系数
xi2=c2/(m2*omega);eta2=2*k/(m2*omega^2); % 方程无量纲化后的系数
G=g/(c*omega^2);M1=(c*m1*omega^2)/(s*P); % 方程无量纲化后的系数
b1=b/c;tau=omega*t;
w=1-epsilon*delta*0.5*F;u=epsilon*0.5*F.*cos(2*tau);l=epsilon*0.5*F.*sin(2*tau);
out=[z(5);...
z(6);...
z(7);...
z(8);...
-(eta1*w+eta1*u)*z(1)-eta1*l*z(2)+(eta1*w+eta1*u)*z(3)+eta1*l*z(4)-xi1*z(5)+fx(z(1),z(2),z(5),z(6))/M1;...
-eta1*l*z(1)-(eta1*w-eta1*u)*z(2)+eta1*l*z(3)+(eta1*w-eta1*u)*z(4)-xi1*z(6)+fy(z(1),z(2),z(5),z(6))/M1-G;...
(eta2*w+eta2*u)*z(1)+eta2*l*z(2)-(eta2*w+eta2*u)*z(3)-eta2*l*z(4)-xi2*z(7)+b1*cos(tau);...
eta2*l*z(1)+(eta2*w-eta2*u)*z(2)-eta2*l*z(3)-(eta2*w-eta2*u)*z(4)-xi2*z(8)+b1*sin(tau)-G];
其中fx、fy为无量纲非线性油膜力
function fxout=fx(x,y,Dx,Dy)
if x-2*Dy~=0
alpha=atan((y+2*Dx)./(x-2*Dy))-0.5*pi*sign((y+2*Dx)./(x-2*Dy))-0.5*pi*sign(y+2*Dx);
else
alpha=atan2((y+2*Dx),(x-2*Dy))-0.5*pi*sign(y+2*Dx)-0.5*pi*sign(y+2*Dx);
end
G=(2/sqrt(1-x.^2-y.^2)).*(pi/2+atan((y.*cos(alpha)-x.*sin(alpha))./sqrt(1-x.^2-y.^2)));
V=(2+(y.*cos(alpha)-x.*sin(alpha)).*G)./(1-x.^2-y.^2);
S=(x.*cos(alpha)+y.*sin(alpha))./(1-(x.*cos(alpha)+y.*sin(alpha)).^2);
fxout=(sqrt((x-2*Dy).^2+(y+2*Dx).^2)./(1-x.^2-y.^2)).*(3.*x.*V-sin(alpha).*G-2.*cos(alpha).*S);
function fyout=fy(x,y,Dx,Dy)
if x-2*Dy~=0
alpha=atan((y+2*Dx)./(x-2*Dy))-0.5*pi*sign((y+2*Dx)./(x-2*Dy))-0.5*pi*sign(y+2*Dx);
else
alpha=atan2((y+2*Dx),(x-2*Dy))-0.5*pi*sign(y+2*Dx)-0.5*pi*sign(y+2*Dx);
end
G=(2/sqrt(1-x.^2-y.^2)).*(pi/2+atan((y.*cos(alpha)-x.*sin(alpha))./sqrt(1-x.^2-y.^2)));
V=(2+(y.*cos(alpha)-x.*sin(alpha)).*G)./(1-x.^2-y.^2);
S=(x.*cos(alpha)+y.*sin(alpha))./(1-(x.*cos(alpha)+y.*sin(alpha)).^2);
fyout=(sqrt((x-2*Dy).^2+(y+2*Dx).^2)./(1-x.^2-y.^2)).*(3.*y.*V+cos(alpha).*G-2.*sin(alpha).*S);
z0=[0 0 0 0 0 0 0 0]'; % 初始值（我自己随便给的）
[t,z]=ode45('liewenfun',[0:0.00025:0.2557],z0,[])
subplot(2,2,1)
plot(t*1000,z(:,3)) % 裂纹转子响应的时序波形
xlabel('t(milliseconds)');
ylabel('x2');
subplot(2,2,2)
plot(z(:,3),z(:,7)) % 裂纹转子响应的相轨迹
xlabel('x2');
ylabel('Dx2');
subplot(2,2,3)
plot(z(:,3),z(:,4)) % 裂纹转子响应的轴心轨迹
xlabel('x2');
ylabel('y2');
fs=4000;N=1024; % 采样频率和采样点数（我自己随便给的）
n=0:1023;
o=fft(z(:,3)); % 对信号进行快速傅里叶变换
mag=abs(o); % 求得fourier变换后的振幅
fre=n*fs/N; % 频率序列
subplot(2,2,4)
plot(fre(1:N/2,mag(1:N/2)) % 绘出Nyquist频率之前随频率变化的振幅
xlabel('频率Hz');
ylabel('振幅');
写的比较乱，麻烦大家帮忙看看啦
|
|
发表于 2010-9-26 14:44
提升卡
变色卡
显身卡
发表于 2010-9-27 22:07
