查看: 875|回复: 6|关注: 0

[已答复] MATLAB ode45求解轴承动力学微分方程组,显示积分公差无法满足

[复制链接]

新手

5 麦片

财富积分


050


3

主题

14

帖子

0

最佳答案
发表于 2019-1-11 19:07:49 | 显示全部楼层 |阅读模式
function fun=ballbearings6(t,y)
global Mri  delta_r theta_2 alpha_2 Cbl alpha_1 Cbr delta_l theta_1 Chol  Chor g  Fx Fy Fz B1 B2 Mhol  Mhor  Khxl  Khyl  Khzl  Khxr  Khyr  Khzr  Ix  Iy  Iz  l1  l2  l  D  f01  v01  Ns   Nb  wc  d    Ke  alpha01  A   deltaep Fcr Fcl  Flx  Fly  Flz  FILx  FILy1  FILy2 Frx Fry Frz FIRx FIRy1 FIRy2
Mri = 3.5;
Cbl = 2.9911*10^3;   %轴承的阻尼
Cbr = 2.9911*10^3;
Chol = 2.9911*10^3;
Chor = 2.9911*10^3;
g = 10;
Fx= 100*sin(wc*t) ;  %径向力
Fy= 0;
Fz= 0;    %轴向力
Mhol = 4.2;   %外圈和轴承座的质量总和
Mhor = 4.2;
Khxl = 2.9911*10^7;   %轴承座刚度
Khyl = 2.9911*10^7;
Khzl = 2.9911*10^7;
Khxr = 2.9911*10^7;
Khyr = 2.9911*10^7;
Khzr = 2.9911*10^7;
Ix = 0.5177;    % 转动惯量
Iy = 0.5177;
Iz = 0.5177;
l1 = 0.0875;   %左边轴承到转子重心的距离
l2 = 0.1275;   %右边轴承到转子质心的距离
l= 0.174;  %轴向力径向力到左边轴承的距离
D = 0.065; %轴承节圆直径
f01 = 1.7; %计算摩擦力时的润滑参数
v01= 0.37;   %油或润滑脂基础油的运动粘度
Ns = 2000;  %轴承转速,以每分钟转数为单位rmp
Nb = 8;     %轴承球的个数
d = 0.015081;  %球的直径
Ke =5.26*10^5; %接触刚度,需要换成油膜总力计算出来的刚度
alpha01= 15;  %初始接触角度
A = 0.588159;   %两个座圈的曲率半径和减去球的直径,需要重新计算
deltaep = 4.48*10e-6;   %需要重新计算,在预载荷作用下的接触变形
B1=0.092718;   %sin(alpha_p),需要重新计算
B2=0.9956924;   %cos(alpha_p),需要重新计算
Zl = deltaep*B1-y(5);
Xl = y(1) + l1*sin(y(21));
Yl = y(3) - l1*sin(y(19));
Zr = deltaep*B1+y(5);
Xr = y(1) - l2*sin(y(21));
Yr= y(3) + l2*sin(y(19));
wc = 776;
Flx = 0;   %初始力都设为零
Fly = 0;
Flz = 0;
FILx = 0;
FILy1 = 0;
FILy2 = 0;
Frx = 0;
Fry = 0;
Frz = 0;
FIRx = 0;
FIRy1 = 0;
FIRy2 = 0;
Mf=2.21*10^(-4); %摩擦力,经过了中间计算,需要重新计算
for i = 1:Nb
   theta_1(i) = wc*t+2*pi*(i-1)/Nb;  %球的位置角
   alpha_1(i)=atan((A*sind(alpha01)+Zl+0.5*D*(y(21)*cos(theta_1(i))+y(19)*sin(theta_1(i))))/(A*B2+deltaep*B2+Xl*cos(theta_1(i))+Yl*sin(theta_1(i))));  %球和座圈的接触角,可以直接用矩阵表示
   delta_l(i)=((A*sind(alpha01)+Zl+0.5*D*(y(21)*cos(theta_1(i))+y(19)*sin(theta_1(i))))^2+(A*B2+deltaep*B2 +Xl*cos(theta_1(i))+Yl*sin(theta_1(i)))^2)^0.5-A; %球和座圈之间的接触变形,需要加入时变位移激励
   Fcl=Ke*(abs(delta_l(i)))^1.5; %接触力,需要用刚度阻尼计算
   Flx= Flx+Fcl*cos(alpha_1(i))*cos(theta_1(i));
   Fly= Fly+Fcl*cos(alpha_1(i))*sin(theta_1(i));
   Flz= Flz+Fcl*sin(alpha_1(i));
   FILx= FILx+Fcl*sin(alpha_1(i))*sin(theta_1(i));
   FILy1= FILy1+Fcl*cos(alpha_1(i))*cos(theta_1(i));
   FILy2= FILy2+Fcl*sin(alpha_1(i))*cos(theta_1(i));
end
for j = 1:1:Nb
  theta_2(j) = wc*t+2*pi*(j-1)/Nb;
  alpha_2(j)=atan((A*sind(alpha01)+Zr-0.5*D*(y(21)*cos(theta_2(j))+y(19)*sin(theta_2(j))))/(A*B2+deltaep*B2+Xr*cos(theta_2(j))+Yr*sin(theta_2(j))));
  delta_r(j)=((A*sind(alpha01)+Zr-0.5*D*(y(21)*cos(theta_2(j))+y(19)*sin(theta_2(j))))^2+(A*B2+deltaep*B2+Xr*cos(theta_2(j))+Yr*sin(theta_2(j)))^2)^0.5-A;
  Fcr=Ke*(abs(delta_r(j)))^1.5;
  Frx= Frx+Fcr*cos(alpha_2(j))*cos(theta_2(j));
  Fry= Fry+Fcr*cos(alpha_2(j))*sin(theta_2(j));
  Frz= Frz+Fcr*sin(alpha_2(j));
  FIRx= FIRx+Fcr*sin(alpha_2(j))*sin(theta_2(j));
  FIRy1= FIRy1+Fcr*cos(alpha_2(j))*cos(theta_2(j));
  FIRy2= FIRy2+Fcr*sin(alpha_2(j))*cos(theta_2(j));

end

      fun(1,1)= y(2);
      fun(2,1)=(Cbl*y(2)-Flx-Frx+Mri*g-Fx)/Mri;
      fun(3,1)=y(4);
      fun(4,1)=(Cbl*y(4)-Fly-Fry-Fy)/Mri;
      fun(5,1)=y(6);
      fun(6,1)=(Cbl*y(6)+Flz-Frz-Fz)/Mri;
      fun(20,1)=((Iy-Iz)*y(22)*y(24)+l1*Fly-l2*Fry+(l+l1)*Fy+0.5*D*FILx-0.5*D*FIRx+Mf)/Ix;
      fun(21,1)=y(22);
      fun(22,1)=((Iz-Ix)*y(24)*y(20)-l1*FILy1+l2*FIRy1-(l+l1)*Fx-0.5*D*FILy2+0.5*D*FIRy2+Mf)/Iy;
      fun(23,1)=y(24);
      fun(24,1)=((Ix-Iy)*y(20)*y(22)+Mf)/Iz;
end
主程序:
clear
clc
globalMri  Cbl Cbr  Chol  Chor g Fx Fy Fz B1 B2 Mhol  Mhor  Khxl Khyl  Khzl  Khxr Khyr  Khzr  Ix Iy  Iz  l1 l2  l  D f01  v01  Ns  Nb  wc  d   Ke  alpha01  A  deltaep  Fcl  Flx Fly  Flz  FILx FILy1  FILy2 Frx Fry Frz FIRxFIRy1 FIRy2  delta_l theta_1 alpha_1 alpha_2theta_2 delta_r t
Mri = 3.5;
Cbl =2.9911*10^3;
Cbr =2.9911*10^3;
Chol =2.9911*10^3;
Chor =2.9911*10^3;
g = 10;
Fx= 100;
Fy= 0;
Fz= 0;
Mhol =4.2;
Mhor =4.2;
Khxl =2.9911*10^7;
Khyl =2.9911*10^7;
Khzl =2.9911*10^7;
Khxr =2.9911*10^7;
Khyr =2.9911*10^7;
Khzr =2.9911*10^7;
Ix =0.05177;
Iy =0.05177;
Iz =0.05177;
l1 =0.0875;
l2 =0.1275;
l= 0.174;
D = 0.065;
f01 = 1.7;
v01= 0.37;
Ns = 2000;
Nb = 8;
d = 0.015081;
Ke =5.26*10^5;
alpha01= 15;
A =0.588159;
deltaep =4.48*10e-6;
B1=0.092718;
B2=0.9956924;
y0=[10^-6,0,10^-6,0,10^-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2000];
[t,y]=ode45('ballbearings6',[0:1],y0);

plot(t,y(:,1));

警告: 在 t=3.208764e-02 处失败。在时间 t 处,若不将步长降至允许的最小值(1.110223e-16)以下,积分公差要求无法满足。
> In ode45 at 308
方程较多,微分方程组哪里,拿掉了几个方程,不影响结果,如果问题出自程序本身,帮忙指教一下,要是数据模型的问题,我再去修改修改。

新手

5 麦片

财富积分


050


3

主题

14

帖子

0

最佳答案
 楼主| 发表于 2019-1-11 20:39:40 | 显示全部楼层
自己顶一个

新手

5 麦片

财富积分


050


3

主题

14

帖子

0

最佳答案
 楼主| 发表于 2019-1-11 20:52:46 | 显示全部楼层
再顶一个

新手

5 麦片

财富积分


050


3

主题

14

帖子

0

最佳答案
 楼主| 发表于 2019-1-12 19:06:00 | 显示全部楼层
顶   顶

入门

51 麦片

财富积分


50500


0

主题

7

帖子

0

最佳答案
发表于 2019-1-13 09:52:01 来自手机 | 显示全部楼层
我之前也出现过类似情况,建议检查微分方程或者先从一个微分方程开始调试,调好了再往上加

新手

5 麦片

财富积分


050


3

主题

14

帖子

0

最佳答案
 楼主| 发表于 2019-1-13 10:36:05 | 显示全部楼层
减少了几个方程,然后发现初值矩阵表达式不对,,应该是个列矩阵

论坛优秀回答者

权威

3516 麦片

财富积分



2

主题

3724

帖子

788

最佳答案
  • 关注者: 165
发表于 2019-1-13 21:36:55 | 显示全部楼层
不会写程序的猪 发表于 2019-1-13 10:36
减少了几个方程,然后发现初值矩阵表达式不对,,应该是个列矩阵

关于这个问题,好像已经有很多的贴。
其实一般的可能是:
1、微分方程本身的问题,比如分母为零、很快增长。
2、初值的问题,看看初值的设置是否正确。
    这个与1有点关联,对于某些初值,积分可能会发散。
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 返回顶部 返回列表