主程序调用子程序提示以下错误
Jeffcott刚性转子，采用Capone油膜力，分析它的时间历程子程序如下：
function dy=fuhuan140922(t,x,flag,wj,p)
uo=0.018; %外油膜粘度,pa.s
ui=0.018; %内油膜粘度,pa.s
co=0.000114; %外油膜间隙,m
ci=0.000093; %内油膜间隙,m
Rro=0.018830; %浮环外径,m
Rri=0.014093; %浮环内径,m
Rj=0.014; %轴颈半径,m
L=0.018; %轴承宽度,m
mj=0.018257/2; %轴子质量之半,kg
mr=0.000079/2; %浮环质量之半,kg
g=9.8; %重力加速度,N/kg(m/s^2)
simga=0.018; %润滑油的动力粘度
wr=wj/(1+(uo/ui)*(ci/co)*(Rro/Rri)^3); %浮环的转速
P=mj*g;
Po=mr*g;
detai=simga*wj*L*Rj*(Rj/ci)^2*(L/(2*Rj))^2; %转子Sommerfeld修正数
detao=simga*wr*L*Rro*(Rro/co)^2*(L/(2*Rro))^2; %浮环Sommerfeld修正数
Mj=mj*ci*wj.*wj./detai;
Mr=mr*co*wr.*wr./detao;
Mk=mr*co*wr.*wr./detai;
Gi=g./(ci*wj.^2);
Go=g./(co*wr.^2);
Xj=x(1);dXj=x(2);Yj=x(3);dYj=x(4);
Xr=x(5);dXr=x(6);Yr=x(7);dYr=x(8);
%以下为求解“外层”油膜力在x,y方向上的无量纲分力——————————————————
ppp1=(Yr+2.00*dXr)/(Xr-2.00*dYr);
sign1=sign(ppp1);
ppp2=Yr+2.00*dXr;
sign2=sign(ppp2);
alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
alphaa=atan((Yr*cos(alpha)-Xr*sin(alpha))/sqrt(abs(1.00-abs(Xr*Xr)-abs(Yr*Yr))));
fgo=2.00*(pi/2.00+alphaa)/sqrt(abs(1.00-abs(Xr*Xr)-abs(Yr*Yr)));
fvo=(2.00+(Yr*cos(alpha)-Xr*sin(alpha))*fgo)/(1.00-abs(Xr*Xr)-abs(Yr*Yr));
fso=(Xr*cos(alpha)+Yr*sin(alpha))/(1.00-abs((Xr*cos(alpha)+Yr*sin(alpha))*(Xr*cos(alpha)+Yr*sin(alpha))));
fao=sqrt(abs(abs((Xr-2.00*dYr)*(Xr-2.00*dYr))+abs((Yr+2.00*dXr)*(Yr+2.00*dXr))))/(1.00-abs(Xr*Xr)-abs(Yr*Yr));
fxo=1.00*fao*(3.00*Xr*fvo-sin(alpha)*fgo-2.00*cos(alpha)*fso);
fyo=1.00*fao*(3.00*Yr*fvo+cos(alpha)*fgo-2.00*sin(alpha)*fso);
%以上为求解“外层”油膜力在x,y方向上的无量纲分力——————————————————
%以下为求解“内层”油膜力在x,y方向上的无量纲分力——————————————————
ppp3=(1.34*Yj+2.00*dXj)/(1.34*Xj-2.00*dYj);
sign3=sign(ppp3);
ppp4=1.34*Yj+2.00*dXj;
sign4=sign(ppp4);
alpha1=atan(ppp3)-pi/2.00*(sign3+sign4);
alphaa1=atan((Yj*cos(alpha1)-Xj*sin(alpha1))/sqrt(abs(1.00-abs(Xj*Xj)-abs(Yj*Yj))));
fgi=2.00*(pi/2.00+alphaa1)/sqrt(abs(1.00-abs(Xj*Xj)-abs(Yj*Yj)));
fvi=(2.00+(Yj*cos(alpha1)-Xj*sin(alpha1))*fgi)/(1.00-abs(Xj*Xj)-abs(Yj*Yj));
fsi=(Xj*cos(alpha1)+Yj*sin(alpha1))/(1.00-abs((Xj*cos(alpha1)+Yj*sin(alpha1))*(Xj*cos(alpha1)+Yj*sin(alpha1))));
fai=sqrt(abs(abs((Xj-2.00*dYj)*(Xj-2.00*dYj))+abs((Yj+2.00*dXj)*(Yj+2.00*dXj))))/(1.00-abs(Xj*Xj)-abs(Yj*Yj));
fxi=-1.00*fai*(3.00*Xj*fvi-sin(alpha1)*fgi-2.00*cos(alpha1)*fsi);
fyi=-1.00*fai*(3.00*Yj*fvi+cos(alpha1)*fgi-2.00*sin(alpha1)*fsi);
%以上为求解“内层”油膜力在x,y方向上的无量纲分力——————————————————
fxi=fxi./detai;fyi=fyi./detai;
fxo=fxo./detao;fyo=fyo./detao;
% 以下为浮环轴承系统动力学模型进行无量纲化处理———————————————————
dy=zeros(8,1);
dy=[x(2);
fxi./Mj+p*sin(t);
x(4);
fyi./Mj+p*cos(t)+Gi;
x(6);
fxo./Mr-fxi./Mk;
x(8);
fyo./Mr-fyi./Mk+Go];
主程序如下：
functionfuhuanjeffcott_solution
p=0; %偏心率
n=12000; %转子的转速
wj=2*pi*n/60; %转子的角频率
%uo=0.018; %外油膜粘度,pa.s
%ui=0.018; %内油膜粘度,pa.s
%co=0.000114; %外油膜间隙,m
%ci=0.000093; %内油膜间隙,m
%Rro=0.018830; %浮环外径,m
%Rri=0.014093; %浮环内径,m
%Rj=0.014; %轴颈半径,m
%L=0.018; %轴承宽度,m
%wr=wj/(1+(uo/ui)*(ci/co)*(Rro/Rri)^3); %浮环的转速
x0=[-0.2 0 -0.2 0 -0.2 0 -0.2 0]; %赋初值,另一组初值;
options=odeset('RelTol',1.e-7,'AbsTol', );%控制精度
T=2*pi; %一个周期
T1=2*pi/wj;
%
figure % 1时间历程
%=ode45('fuhuan140922',,x0,options，wj,p);
=ode45('fuhuan140922',,x0,options,wj,p);
subplot(121)
plot(t(100000:end),X(100000:end,1))
xlabel('t')
ylabel('Yj') %grid on
title('轴颈时间历程x-t')
提示错误如下：
Warning: Failure at t=6.250756e-004.Unable to meet integration tolerances without
reducing the step size below the smallest value allowed (1.734723e-018) at time t.
> In ode45 at 371
In fuhuanjeffcott_solution at 24
可能是你option中的收敛误差设置的太小了，积分残差不满足设置的值！
页:
[1]