ode方程右端含左侧部分导数的求解问题

chunshui2003

模型和具体内容来自于文献“褚福磊 and 唐云, 碰摩转子系统的稳定性. 清华大学学报：自然科学版, 2000. 40(4): p. 119-123.”。系统为两段轴承，中间转子的常见结构。对时间t进行无量纲化，并作出假设如下：

则运动微分方程变为

其中，碰摩力为

因为s2’和s4’中含有s5’和s6’，所以首先计算s5’和s6’，即ode方程先让s5’和s6’出现，表达如下

function uu = equ_0724(T,u)

m = 2;

c = 580;

k = 60*10^3;

e0 = 40*10^-6;

delta = 0.1*10^-3;

f = 0.2;

kc = 9*10^6;

mu = 15*10^-3;

R = 40*10^-3;

L = 10*10^-3;

c1 = 100*10^-6;

g = 9.81;

global w

e = sqrt(u(1).^2 + u(3).^2);

H = 1/2 * (sign(abs(e - delta)) + sign(e - delta));

Fx = -(e - delta)*kc./e * (u(1) - f.* u(3)) * H;

Fy = -(e - delta)*kc./e * (f.* u(1) + u(3)) * H;

D1 = pi*mu*R*L^3;

D2 = (c1^2 - u(5).^2 - u(6).^2).^(3/2);

D3 = (c1^2 - u(5).^2 - u(6).^2).^(5/2);

D4 = 2*c*w*D2*D3 + 2*w*D1*D3 + 6*w*D1*D2*u(5).*u(5);

D5 = 6*w*D1*D2*u(5)*u(6);

D6 = 2*c*w*D2*D3*u(2) - 2*D2*D3*k*(u(5) - u(1)) - w*D1*D3*u(6);

D7 = 2*c*w*D2*D3 + 2*w*D1*D3 + 6*w*D1*D2*u(6).*u(6);

D8 = 2*c*w*D2*D3*u(4) - 2*D2*D3*k*(u(6) - u(3)) + w*D1*D3*u(5);

uu = zeros(6,1);

uu(5) = (D6.*D7 - D5.*D8) ./ (D4.*D7 - D5.*D5);

uu(6) = (D4.*D8 - D5.*D6) ./ (D4.*D7 - D5.*D5);

uu(1) = u(2);

uu(2) = -c/(m*w) * (u(2) - uu(5)) - k/(m*w^2) * (u(1) - u(5)) + Fx/(m*w^2) + e0*cos(T);

uu(3) = u(4);

uu(4) = -c/(m*w) * (u(4) - uu(6)) - k/(m*w^2) * (u(3) - u(6)) + Fy/(m*w^2) + e0*sin(T) - g/w^2;

  

求解方程

period = 2*pi;

step = period/200;

tspan = 0:step:1000*period;

y0 = [0.05*10^-3 0.05*10^-3 0.05*10^-3 0.05*10^-3 0.05*10^-3 0.05*10^-3];

global w

w = 5300*(2*pi/60);

[t,u] = ode45('equ_0724',tspan,y0);

但是据此计算得到的轨迹图和映射图同文献相比有着比较明显的差距，特别是映射图（我采用的是频闪法，文献的方法不清楚），基本对不上号。如果按正常顺序把s5’和s6’放在末尾，即uu(5)和uu(6)放在方程末尾，同样无法得到满意的结果。

我的结果：没有经过无量纲化，量纲基准为delta = 0.1mm。数量级可以同文献对应上。

文献的结果

想请教一下大家，我的方程写的有没有问题。如果没有的话为什么结果差距这么大，是哪一个环节出现了问题。谢谢。

伤痕累累

哥哥，你可是这方面的专家了，关注中......