ode45求滚动轴承支撑下碰摩转子动力学方程，程序报错。
%%%%动力学方程程序function dx=gd_pm1(t,x)
dx=zeros(12,1);
global w
%%求滚动轴承非线性作用力
global Nb Kb
Nb=8; %滚珠个数
Kb=13.34e+09; %轴承滚珠和滚道的接触刚度系数单位 N/m^3/2 8.9176
rb=0.0401; %轴承内圈半径单位m
Rb=0.0639; %轴承外圈半径单位m
w0=2*pi*(w/60); %转速和角速度换算公式
wb=w0*rb*(rb+Rb); %轴承保持架的角速度
c=0.000008; %轴承游隙
%左侧轴承1的非线性作用力；轴承1几何中心坐标（x2,y2）写为（x（5）,x（7））
= Force_bearing(t,x(5),x(7),c,wb);
%右侧轴承2的非线性作用力；轴承2几何中心坐标（x3,y3）写为（x（9）,x（11））
= Force_bearing(t,x(9),x(11),c,wb);
%%系统相关系数
m1=32.1; %圆盘处集中质量
m2=4; %轴承1处、轴承2处集中质量
g=9.8;
c1=2100; %圆盘处阻尼
c2=1050; %轴承处阻尼
k0=0.85e7; %轴系刚度系数
e=0.00005; %质量偏心距离
%%无量纲化
% wn=(k0/m1)^0.5;
kesai1=c1/(m1*w); %无量纲阻尼
kesai2=c2/(m2*w);
numda1=k0/(m1*w^2);
numda2=k0/(m2*w^2);
gg=g/(c*w^2);
U=e/c;
% omega=w/wn;
beta=0;
%%求碰摩力
ku=0.1; %摩擦系数
kc=3.5e7; %碰摩刚度系数
r=sqrt(x(1)^2+x(3)^2);%圆盘几何中心位移
deta=0.00002; %碰摩间隙
H=0.5*(sign(abs(r-deta))+sign(r-deta));%判断是否发生碰摩，sigh正数为1，负数为-1
px=H*(-c*(r-deta)*kc*(x(1)-ku*x(3))/r);
py=H*(-c*(r-deta)*kc*(ku*x(1)+x(3))/r);
%%方程
dx(1)=x(2);
dx(2)=-kesai1*x(2)-numda1*(x(1)-x(5))-numda1*(x(1)-x(9))+U*cos(t-beta)+px/(w^2*m1*c);
dx(3)=x(4);
dx(4)=-kesai1*x(4)-numda1*(x(3)-x(7))-numda1*(x(3)-x(11))+U*sin(t-beta)+py/(w^2*m1*c)-gg;
dx(5)=x(6);
dx(6)=-kesai2*x(6)-numda2*(x(5)-x(1))+Fx_2/(w^2*m2*c);
dx(7)=x(8);
dx(8)=-kesai2*x(8)-numda2*(x(7)-x(3))+Fy_2/(w^2*m2*c)-gg;
dx(9)=x(10);
dx(10)=-kesai1*x(10)-numda1*(x(9)-x(1))+Fx_3/(w^2*m2*c);
dx(11)=x(12);
dx(12)=-kesai1*x(12)-numda1*(x(11)-x(3))+Fy_3/(w^2*m2*c)-gg;
%%%%求解程序
tic
clear all
clc
global w
w=400;
T=2*pi;
delt_T=T/100;
x0=;
tspan=;
options = odeset('RelTol', 1e-6,'AbsTol',1e-6);
=ode45('gd_pm1',tspan,x0,options);
x0=x(end,:);
tspan=;
=ode45('gd_pm1',tspan,x0,options);
figure
%%time domain plot
subplot(4,1,1)
u=[];
u=x(1:end,1)-mean(x(1:end,1));
plot(t(1:end),u(1:end,1));%,set(gca,'xlim',)
xlabel('t/s');
ylabel('x');
title('(a)')
subplot(4,1,2)
%%axis orbit
plot(x(1:end,1),x(1:end,3))
xlabel('x');
ylabel('y');
title('(b)')
%% FFT
fs=100/(2*pi);
N = length(x(1:end,1));
df = 0:fs/N:(N-1)*fs/N;
a= abs(fft(x(1:end,1)-mean(x(1:end,1))))*2/N;
DF=2*pi*df;
subplot(4,1,3)
plot(DF,a),set(gca,'xlim',)
xlabel('f/fr');
ylabel('\rmAmplitude');
title('(c)')
%%poincare map
m=[];
n=[];
m=x(1:100:end,1);
n=x(1:100:end,2);
subplot(4,1,4)
plot(m,n,'k.')
xlabel('x');
ylabel('x’');
title('(d)')
%%%%轴承力程序
function = Force_bearing(t,x,y,c,wb)%该函数用于求滚动轴承的非线性支撑力
globalNb Kb
theta_ball=zeros(1,Nb);
for i=1:Nb
theta_ball(i) = 2*pi*(i-1)/Nb; %各转动体初始角
end
theta_change=wb*t; %转动体角度变动量
theta_ball=theta_ball+theta_change;%各滚珠的各时变转角
sigma_ball= x*cos(theta_ball)+y*sin(theta_ball)-c;%各滚珠的法向接触变形量
Heaviside_1 =heaviside(sigma_ball); % heaviside 函数用于筛选 sigma_ball>0的变形量
useful_sigma = sigma_ball.*Heaviside_1;
sigma_1=useful_sigma.^3;
sigma = sqrt(sigma_1);
fix_bearing = (sigma.*cos(theta_ball));%轴承支撑力的x方向分量
fiy_bearing = (sigma.*sin(theta_ball));%轴承支撑力的y方向分量
Fx_bearing=Kb * sum(fix_bearing);
Fy_bearing=Kb * sum(fiy_bearing);
end
将三段程序分别保存为3个m文件，运用求解程序，出现报错。求助坛内大神，小弟应该做如何修改。
先行谢过了！！！
页:
[1]