转子运动微分方程：MX ̈+KX ̇+CX=Q的求解

伤痕累累

回复 15 # tifanhuang 的帖子

贴出来，让大家看看

tifanhuang

求解的微分方程组（在附件中）  

其中未知量为 , , .

以下是我用的四阶龙格—库塔算法求解的程序：

第一部分是建立函数：

function dy = zunifun(t,y )
%此为两刚性支撑下悬臂转子模型
%以下为求解转轴刚度
E0=2.1e5;%转轴弹性模量
a=5.5e-1;
l=1.0e2;
d=1.0e1;%转轴直径
j=1/(64*pi*d^4);
x11=-(a^2)*(l+a)/(3*E0*j);
x12=a*(2*l+3*a)/(6*E0*j);
x21=x12;
x22=(l-3*a)/(3*E0*j);
A=[x11 x12;x21 x22];
K=inv(A);
k11=K(1,1);
k12=K(1,2);
k21=K(2,1);
k22=K(2,2);
Jp=2.7596e-11;%转子极转动惯量
Jd=0.5*Jp;%转子直径转动惯量
Cr=0.01;%转子阻尼系数
Ct=0.01;%转子阻尼系数
Cb=0.3;%外阻尼系数
m=8.64e-2;%转子质量
eu=5;%偏心距
f=20;%转子频率
Q=2*pi*f*t;
g=9.98;% 重力加速度
%转子运动微分方程组
dy=zeros(8,1)%行向量
dy(1)=y(2);
dy(2)=eu*(2*pi*f)^2*cos(Q)-((Ct+Cb)*y(2)+k11*y(1)+k12*y(6)+m*g)/m;
dy(3)=y(4);
dy(4)=eu*(2*pi*f)^2*sin(Q)-((Ct+Cb)*y(4)+k11*y(3)-k12*y(7))/m;
dy(5)=y(6);
dy(6)=(Jp*2*pi*f*y(8)-Cr*y(6)-k21*y(1)-k22*y(5))/Jd;
dy(7)=y(8);
dy(8)=-(Jp*2*pi*f*y(6)+Cr*y(8)-k21*y(3)+k22*y(7))/Jd;

end

第二部分是求解过程：

clc;
t0=0;
tN=2;
y0=[0.001;0.001;0.001;0.001;0.001;0.001;0.001;0.001];
h=0.001;
t=t0:h:tN;
N=length(t);
j=1;
for i=1:N
    t1=t0+h;
    K1=zunifun(t0,y0);
    K2=zunifun(t0+h/2,y0+h*K1/2);
    K3=zunifun(t0+h/2,y0+h*K2/2);
    K4=zunifun(t0+h,y0+h*K3);
    y1=y0+(h/6)*(K1+2*K2+2*K3+K4);
    yy1(j)=y1(1);
    yy2(j)=y1(2);
    yy3(j)=y1(3);
    yy4(j)=y1(4);
    yy5(j)=y1(5);
    yy6(j)=y1(6);
    yy7(j)=y1(7);
    yy8(j)=y1(8);
    t0=t1;
    y0=y1;
    j=j+1;
end

出来的结果全是0，希望各高手能帮忙指出其中的错误，谢谢。

伤痕累累

回复 17 # tifanhuang 的帖子

你为什么不用ode45来解呢、控制一下时间步长，也算是定步长啊。这样显得有点麻烦

tifanhuang

我也用了，可一直出现电脑在跳动

tifanhuang

如果ode45可以求解，那ode45跟四阶龙格—库塔算法也应该是差不多的，只是精度问题，而我的现在求的结果全为0

tifanhuang

以上的程序有哪位大神帮帮忙呗，是在纠结呀

hnczf739892267

我看了两遍，怎么没有找到你的负载和驱动力矩啊