请教各位大侠转子-滑动轴承的碰摩问题

yangxiaokang

建立的动力学模型在附件1:
程序在附件2、3中，
怎么计算得不出来图形，参数是参照别人的选取的，不知道判断碰摩时是否正确？
请各位大侠请教，谢谢！

yangxiaokang

方程计算得不到图形，也没有报错信息，模型和方程在附件中，请各位不吝赐教
%%%%%%%%%%%%%%%%%%%%%%%%%%
%以下是求解方程
function dy=fangcheng(t,x)
%参数
global w kc
m1=4;%轴承处集中质量
m2=32.1;%圆盘处集中质量
R=0.025;%轴承半径
L=0.012;%轴承长度
c=0.00011;%油膜厚度
mu=0.018;%油膜黏度
f=0.1;%摩擦系数
c1=1050;%轴承阻尼
c2=2100;%圆盘阻尼
k=2.5e7;%弹性轴刚度
r=0.00005;%圆盘质量偏心距
delta=0.000016;%转子与定子间隙
g=9.8;
d=r/c;
G1=g/(c*w^2);
G2=g/(c*w^2);
e=sqrt(x(5)^2+x(7)^2); %圆盘中心的径向位移
%判断是否发生碰摩
H=1/2*(sign(abs(e-delta))+sign(e-delta)); %判断是否发生碰摩，sign正数返回1，负数返回-1
%碰摩力
Px=-c*(e-delta)*kc/e*(x(5)-f*x(7))*H;
Py=-c*(e-delta)*kc/e*(f*x(5)+x(7))*H;
%非线性油膜力
s=mu*w*R*L*(R/c)^2*(L/(2*R))^2;  %短轴承油膜力的sommerfeld数
a=atan((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign(x(3)+2*x(2));
G=2/sqrt(1-x(1)^2-x(3)^2)*(pi/2+atan((x(3)*cos(a)-x(1)*sin(a))/sqrt(1-x(1)^2-x(3)^2)));
S=(x(1)*cos(a)+x(3)*sin(a))/(1-(x(1)*cos(a)+x(3)*sin(a))^2);
V=(2+(x(3)*cos(a)-x(1)*sin(a))*G)/(1-x(1)^2-x(3)^2);
fx=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(1)*V-sin(a)*G-2*cos(a)*S);
fy=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(3)*V+cos(a)*G-2*sin(a)*S);
%方程组
dy=[x(2);
    -c1/(w*m1)*x(2)-k/(w^2*m1)*(x(1)-x(5))+s/(c*w^2*m1)*fx;
    x(4);
    -c1/(w*m1)*x(4)-k/(w^2*m1)*(x(3)-x(7))+s/(c*w^2*m1)*fy-G1;
    x(6);
    -c2/(w*m2)*x(6)-2*k/(w^2*m2)*(x(5)-x(1))+Px/(c*w^2*m2)+d*cos(t);
    x(8);
    -c2/(w*m2)*x(8)-2*k/(w^2*m2)*(x(7)-x(3))+Py/(c*w^2*m2)+d*sin(t)-G2];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
下面是主程序
%主程序
clear all
clc
global w kc
w=900;   %转子角速度
kc=3.5e6;  %定子径向刚度

x0=[0;0.1;0;0.1;0;0.1;0;0.1];  %初值
options=odeset;
options.reltol=1e-6;
T=2*pi/w;
[t,y]=ode45(@fangcheng, [0:T/100:1000*T],x0);
plot(t,y(:,1))
figure
plot(y(:,5),y(:,7)) %圆盘中心轨迹图

无水1324

你把计算时间变短一点100T算一下，看是不是有输出。另外你最好将方程量纲一化，这样求解要好点，因为你计算的时候步长要求很小，已经到了2pi/90000。

yangxiaokang

谢谢，把时间变为100T后计算，还是没有输出，Y的输出为不定值NaN。方程已经无量纲化过了，另外有谁做过转子碰摩的，这样
%判断是否碰摩
H=1/2*(sign(abs(e-delta))+sign(e-delta));
%碰摩力
Px=-c*(e-delta)*kc/e*(x(5)-f*x(7))*H;
Py=-c*(e-delta)*kc/e*(f*x(5)+x(7))*H;
判断碰摩力对吗:@)

无水1324

你的方程无量纲化了？:@L 。我还是看不懂那里处理了。

yangxiaokang

首先，谢谢无水的帮助，感谢在这里有个向各位学习交流的平台，我刚刚学习转子动力学不久，只是看过一些资料，具体没做过，还只多多向各位学习。我把我的无量纲化过程说下，也许是错的，恳请各位多多指导，主要是为了解决问题嘛，谢谢！
1，碰摩力无量纲化在发件1中
2，油膜力选取的是无量纲化过的短轴承油膜力模型在附件2中
3，方程的无量纲化在附件3，4中，其中只做了轴承处x1,y1的。

无水1324

恩，现在我明白你的无量纲化了。


clear all
clc
global w kc
w=900;   %转子角速度
kc=3.5e6;  %定子径向刚度

x0=[0;0.1;0;0.1;0;0.1;0;0.1];  %初值
options=odeset;
options.reltol=1e-6;
T=2*pi/w;
[t,y]=ode45(@fangcheng, [0:T/100:1000*T],x0);
plot(t,y(:,1))
figure
plot(y(:,5),y(:,7)) %圆盘中心轨迹图


此时方程的激励周期T的计算错了，因为后面方程中都是sin t了，所以周期是2pi，你修改一下

yangxiaokang

将周期改为2pi，步长取T/1024后还是没有结果，不清楚步长该如何取合理，是不是前面参数有问题，其他人也可以发表意见啊，集思广益嘛！:@(
程序改为下面的：
%主程序
clear all
clc
global w kc
w=900;
kc=2.5e7;

x0=[0;0.1;0;0.1;0;0.1;0;0.1];  %初值
options=odeset;
options.reltol=1e-6;
T=2*pi;
[t,y]=ode45(@fangcheng, [0:T/1024:100*T],x0);
plot(t,y(:,1))
figure
plot(y(:,5),y(:,7))

chunshui2003

发生碰磨的条件是不是应该写在主程序里啊，另外
e=sqrt(x(5)^2+x(7)^2); %圆盘中心的径向位移
%判断是否发生碰摩
H=1/2*(sign(abs(e-delta))+sign(e-delta)); %判断是否发生碰摩，sign正数返回1，负数返回-1
这部分我不明白为什么这样写，如果是判断碰磨那么只要e>delta就可以了啊！也不知道理解是否正确，最近分岔程序搞得有点糊涂了。

yangxiaokang

谢谢，chunshui2003 的建议。我想判断碰摩条件是在方程中，因为碰摩力当e>delta时有，当e<或者=delta时没有，方程是随着e=sqrt(x(5)^2+x(7)^2); %圆盘中心的径向位移变化而变化的。可能是里面的取值有问题。
下面是我试着计算下
clear all
clc
global w kc
w=400;
kc=2.5e7;

x0=[0.02;0.1;0.01;0.1;0.02;0.1;0.01;0.1];  %初值
options=odeset; %确定求解的条件和要求
options.reltol=1e-6;  %限定相对误差
T=2*pi;
[t,y]=ode45(@fangcheng, [0:T/1000:100*T],x0);  %试过很多解法都不行
plot(t,y(:,1));  title('轴颈中心位移x时程图');  %轴颈中心位移x随时间的变化
figure
plot(y(:,5),y(:,7)); title('圆盘中心轨迹图');  %圆盘中心轨迹图

得到的图像在附件中，没能计算完，报警说积分最小误差不能满足，请各位帮忙看看，急啊

yangxiaokang

下面是分岔图程序，请各位高手帮忙看看，好像有问题
clear all
clc
w=400:5:600;
for i=1:length(w)
disp(w(i));
T=2*pi/w(i);
tspan=[0:T/100:100*T];
x0=[0.02;0.1;0.01;0.1;0.02;0.1;0.01;0.1];  %初值
[t,y]=ode45(@fangcheng,tspan,x0,[ ],w(i));
%y0=x(end,:);  %这里不明白
plot(w(i),y(8000:50:end,2),'markersize',20);
hold on;
xlabel('频率');
ylabel('位移');
end

chunshui2003

肯定是有问题的

tspan=[0:T/100:100*T];

plot(w(i),y(8000:50:end,2),'markersize',20);

时间间隔是100，但你取值是50，最起码应该是
plot(w(i),y(8000:100:end,2),'markersize',20);

另外，一共才10000个点，从8000取是不是有点晚了，过滤瞬态也不用这么久啊。markersize好像也太大了。
呵呵，我理解就是这样。

yangxiaokang

我按chunshui2003说的计算了下，还是不行啊，计算得到的点很小，即使取到20还是很小。分岔图应该是频率或者频率比随着位移的变化吧，请问这里面频率的变化范围怎么取合理啊，还有在做转子碰摩非线性分析时，做轨迹图、庞加莱映射图、时程图、频谱图、还有分岔图，应该有先后顺序吧，请高手赐教
w=400:5:600;
for i=1:length(w)
disp(w(i));
T=2*pi/w(i);
tspan=[0:T/100:100*T];
x0=[0.02;0.1;0.01;0.1;0.02;0.1;0.01;0.1];  %初值
[t,y]=ode45(@fangcheng,tspan,x0,[ ],w(i));
%y0=x(end,:);  %这里不明白
plot(w(i),y(4000:100:end,1),'markersize',2);%这里应该是频率随位移变化的图形

hold on;
xlabel('频率');
ylabel('位移');
end

yangxiaokang

我自己认为是先做参数随频率的分岔图，然后在分岔出现的地方，取不同频率，再做轨迹图和映射图分析，不知道对否，个人拙见

无水1324

另外不知道你的分岔图怎么样了？现在对了吗？
“我自己认为是先做参数随频率的分岔图，然后在分岔出现的地方，取不同频率，再做轨迹图和映射图分析，不知道对否”
这个想法是对的，但是有时候不好选择参数的范围的时候，可以先尝试着计算一些点的轨迹图和映射图

[ 本帖最后由 无水1324 于 2010-3-10 16:21 编辑 ]