[已答复] 滚珠球轴承转子系统计算出现:求教。代码在附件 警告: 在 t=2.521361e-03 处失败。在时间 t 处,步长大小必须降至所允许的最小值(6.938894e-18)以下,才...

[复制链接]
okzhu 发表于 2022-4-25 17:56:29
滚珠球轴承转子系统计算出现:求教。代码在附件
警告: 在 t=2.521361e-03 处失败。在时间 t 处,步长大小必须降至所允许的最小值(6.938894e-18)以下,才能达到积分容差要求。
> In ode45 (line 360)
  In main (line 13)
clear;
clc;
global n
n=300;      %转速
y0=[0;1e-4;0;1e-4;0;1e-4;0;1e-4;0;1e-4;0;1e-4];
w=2*pi*(n/60);      %转速和角速度换算
T=2*pi/w;
dt=T/512;
tspan=0:dt:1.2;    %时间范围
%tspan=[0:0.001:1.2];
%options = odeset('MaxStep',1.75e-15);
options = odeset('RelTol', 1e-6,'AbsTol',1e-6);%相对误差和绝对误差
[t,y]=ode45('fangcheng',tspan,y0 ,options);
figure
plot(t,y(:,1));

function [fx, fy]=Force_bearing(t,x,y,r0,wcage)
global  Nb cb
theta_j=zeros(1,Nb);
for i =1:Nb
theta_j(i)=wcage*t+2*pi*(i-1)/Nb;               %各滚珠的各时转角
end
delta_j=x.*cos(theta_j)+y.*sin(theta_j)-r0;   %各滚珠的法向接触变形量
Heaviside_1 =  heaviside(delta_j);          %用于筛选 delta_j>0的变形量
useful_delta_j=delta_j.* Heaviside_1;
fjx=sqrt(useful_delta_j.^3).*cos(theta_j);%轴承支撑力的x方向分量
fjy=sqrt(useful_delta_j.^3).*sin(theta_j); %轴承支撑力的y方向分量
fy=cb*sum(fjy);
fx=cb*sum(fjx);

function dy=fangcheng (t,y)
dy=zeros(12,1);
global n Nb cb;
Nb=8;               %滚珠个数
cb=13.34e+09;       %轴承滚珠和滚道的接触刚度系数  单位 N/m^3/2
r=0.0401;           %轴承内圈半径  单位m
R=0.0639;           %轴承外圈半径  单位m
w=2*pi*(n/60);      %转速和角速度换算
wcage=w*r/(r+R);     %轴承保持架的角速度
r0=4e-6;          %轴承游隙
%右侧轴承的非线性作用力;轴承几何中心坐标为(y(5),y(7))
[fxr fyr]=Force_bearing(t,y(5),y(7),r0,wcage);
%左侧轴承的非线性作用力;轴承几何中心坐标为(y(9),y(11))
[fxl fyl]=Force_bearing(t,y(9),y(11),r0,wcage);

mrp=8;             %圆盘处集中质量
mrl=2;mrr=2;       %左轴承处、右轴承处集中质量
g=9.8;
crp=2940;           %圆盘处阻尼
crb=200;            %轴承处阻尼
k=2.5e8;            %轴系刚度系数
e=0.00001;          %质量偏心距离


dy(1)=y(2);
dy(2)=e*(w^2)*cos(w*t)-crp*y(2)/mrp-k*(y(1)-y(5))/mrp-k*(y(1)-y(9))/mrp;
dy(3)=y(4);
dy(4)=e*(w^2)*sin(w*t)-g-crp*y(4)/mrp-k*(y(3)-y(7))/mrp-k*(y(3)-y(11))/mrp;
dy(5)= y(6);
dy(6)=fxr/mrr-crb*y(6)/mrr-k*(y(5)-y(1))/mrr;
dy(7)= y(8);
dy(8)=fyr/mrr-g-crb*y(8)/mrr-k*(y(7)-y(3))/mrr;
dy(9)= y(10);
dy(10)=fxl/mrl-crb*y(10)/mrl-k*(y(9)-y(1))/mrl;
dy(11)=y(12);
dy(12)= fyl/mrl-g-crb*y(12)/mrl-k*(y(11)-y(3))/mrl;


fangcheng.m

1.22 KB, 下载次数: 3

main.m

385 Bytes, 下载次数: 3

fangcheng.m

1.22 KB, 下载次数: 3

2 条回复


lvcggg 发表于 2022-4-28 00:06:28
程序正负号有误

mzj2 发表于 4 天前
楼主,解决了吗,同问,谢谢了!!
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表