pengmoweifen.zip_动力学仿真_单转子动力学_碰摩_简支_转子-轴承

%%%%%%%%%%%%%%%%%%%%%%% %碰撞微分方程： % %%%%%%%%%%%%%%%%%%%%%%%% clear all; clc; clf; omega=900; %想求系统随omega（转速）的变化分岔图和截面图，这里想用一个定值试试 period=2*pi; tspan = (1:period/100:500*period); y0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001]; [t,u] = ode45(@equ_elastic_without_UMP,tspan,y0,u); %%%%%微分方程 function uu=equ_elastic_without_UMP(t,flag,u) m1=4; %转子质量Kg m2=32; %轴颈质量Kg R=5*10^-2; %轴承半径m L=2*10^-2; %轴承长m miu=0.11; %平均油膜厚度 sigam=0.2; %转子静子间隙 mu=0.018; %润滑油粘度 c1=1050; %转子阻尼N.s/m c2=2100; %轴承阻尼N.s/m f=0.1; %摩擦系数 kc=6.2*10^6; %静子刚度 k=2.5*10^7; %弹性轴刚度 e=0.05; %转子质量偏心m %delta_0 = 4.5*10^-3; %均匀气隙大小 m %mu = 18*10^-3; %绝对润滑油粘度 无单位 %omega = 13; %大轴旋转角速度 %global cz %轴承间隙 global omega %%% 非线性油膜力表达式 fx fy %%%% M=m2/2; sigma = mu*omega*R*L/M*(R/miu)^2*(L/2/R)^2; %Sommerfeld修正系数%%%%%%%%改动(*2/m2) A1=u(7)+2*u(6); A2=u(5)-2*u(8); alpha=atan(A1./A2)-pi/2*sign(A1./A2)-pi/2* sign(A1); E=sqrt(u(5).^2 +u(7).^2); %轴颈轴心轨迹 E_minus=sqrt(1-E.^2); B1=u(7)*cos(alpha)-u(5)*sin(alpha); B2=u(5)*cos(alpha)+u(7)*sin(alpha); G=2./E_minus.^2*( pi/2+atan(B1./E_minus)); %此处错误，应该是E_minus.^2; V=(2 + B1 * G)./E_minus.^2; S=B2./(1-B2.^2); F_oil_same=-sqrt(A2.^2+A1.^2) ./ E_minus.^2; f_x=F_oil_same*(3*u(5)*V-sin(alpha)*G-2*cos(alpha)*S); f_y=F_oil_same*(3*u(7)*V+cos(alpha)*G-2*sin(alpha)*S); %油膜力的无量纲表达式 fx=sigma*f_x; fy=sigma*f_y; %非线性油膜力 %%%摩擦力计算%%%%%%%%%%%%%%%%%%%%%%%%% r=sqrt(u(5).^2+u(7).^2); PX=-kc* miu.*(1-sigam/r).*(u(5)-f.*u(7)); PY=-kc* miu.*(1-sigam/r).*(f.*u(5)+u(7)); g=9.8; G=g/(miu*omega^2); b=e/miu; tao=omega*t; uu = zeros(8,1); uu(1) = u(2); uu(2) = -c1/(m1*omega) * u(2) - k/(m1*omega^2) * (u(1)-u(5)) + sigma*M*fx/(m1*omega^2*miu); uu(3) = u(4); uu(4) = -c1/(m1*omega) * u(4) - k/(m1*omega^2) * (u(3)-u(7)) + sigma*M*fy/(m1*omega^2*miu)-G; uu(5) = u(6); uu(6) = -c2/(m2*omega) * u(6) - 2*k/(m2*omega^2)* (u(5)-u(1)) + PX/(m2*omega^2*miu) +b*cos(tao); uu(7) = u(8); uu(8) = -c2/(m2*omega) * u(8) - 2*k /(m2*omega^2) * (u(7)-u(3)) + PY/(m2*omega^2*miu) +b*sin(tao)-G;