滑动轴承做时间历程、轴心轨迹、相图、庞加莱映射程序

dj15175771160

我以滑动轴承为例做时间历程、轴心轨迹、相图、庞加莱映射图，程序如下，但是变换转速和偏心率，做出来的图总感觉不对，希望大家运行下看看你，滑动轴承采用的是capone油膜力，对建立的转子动力学方程进行了无量纲化。
function  jeffcottzhuanzixitong_solution
clear all
clc
n=10000;
w=pi*n/30;
%w=[240,300,400,550,620,700,900];
%x0=[0.42;0;0;0];     
x0=[0;0.1;0;0.1];
options=odeset('RelTol',1.e-7,'AbsTol',[1.e-8 1.e-8 1.e-8 1.e-8] );
% %for j=1:length(w)
%     T=2*pi;     %一个周期
%     [t,X]=ode45('force11',[0:T/100:140*T],x0,options);
%     figure    % 1时间历程x-tau
%     %subplot(221)
%     plot(t,X(:,3))
%     title('时间历程y-t')
%     xlabel('t'),  ylabel('Y')    %  grid on

    % 2相图x-x'
%     figure
%     %subplot(222)
%     T1=2*pi/w;
%     [t,X]=ode45('force11',[0:T1/100:14000*T1],x0,options)
%     plot(X(1150000:end,3),X(1150000:end,4),'b')
%     title('相图')
%     xlabel('Y'),  ylabel('dY')         %  grid on
   % axis([-1 1 -1 1])
% % 3轴心轨迹
    figure
    %subplot(223)
    T1=2*pi/w;
    [t,X]=ode45('force11',[0:T1/100:14000*T1],x0,options)
    plot(X(1000:end,1),X(1000:end,3),'b')
    title('轴心轨迹')
    xlabel('X'),ylabel('Y')          %  grid on
    axis([-1 1 -1 1])

% % 4 ***绘制poincae截面图
%     figure
%     %subplot(224)
%     T1=2*pi;  %周期
%     [t1,X1]=ode45('force11',[0:T1/100:14000*T1],x0,options);
%     i=1000000:100:length(X1);
%     plot(X1(i,1),X1(i,2),'*')
%     axis([-1 1 -1 1])
%     title('poincare截面')
%   
%     % 5***频谱分析fft对数值解作傅立叶变换，功率谱分析
%     disp('对数值解作傅立叶变换，傅立叶功率分析')
%     T1=2*pi/w;
%     [t,X]=ode45('force11',[0:T1/100:14000*T1],x0,options)
%     tt=100/T1;
%     x=X(:,3);
%     Y=fft(x);   %%对数值解作傅立叶变换%傅立叶功率分析
%     Y(1)=[];    %%去掉零频分量
%     
%     n=length(Y)/2;   %计算不同频率的个数
%     power=abs(Y(1:n)).^2/(length(Y).^2);  %计算功率谱
%     freq=tt*(1:n)/length(Y);         %%计算频率，因为步长为0.01，而不是1，故乘以100
%     power1=power ;               %%寻找两个最大的频率
%     [id1,daa1]=max(power1);
%     power1(daa1)=0;
%     [id2,daa2]=max(power1);
%     WF=2*pi*[freq(daa1),freq(daa2)];%求圆频率
%
%    figure %画功率谱图
%     plot(2*pi*freq(1:100),power(1:100))   %%画功率谱图
%     %hold on
%     title('功率谱')
%     xlabel('\omega/s^{-1}');  ylabel('p/m^2');
%    % plot(WF,[power(daa1),power(daa2)],'r.','markersize',30)
%


function dxdt=force11(t,x,flag, w)
r=0.0118;L=0.022; c=0.0001;
m=0.6689/2;   g=9.8;   
p=0.25;      %p=0.068，0.1， 0.2， 0.3， 0.4 ，0.5; p为偏心率
u=0.018;     %u=η为润滑油动力粘度
a1=c/r;        % w=ω为转速，
n=5000;
w=pi*n/30;
M=m*c*w*a1^2/(L*r*u);    %令 x(1)=x, x(2)=x', x(3)=y, x(4)=y'
G1=g/(c*w^2);  %a1为间隙比，L为轴瓦宽度，c为轴承半径间隙，r为轴承半径
a=atan((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign(x(3)+2*x(2));
G=2/sqrt(1-x(1)^2-x(3)^2)*(pi/2+atan((x(3)*cos(a)-x(1)*sin(a))/sqrt(1-x(1)^2-x(3)^2)));
S=(x(1)*cos(a)+x(3)*sin(a))/(1-(x(1)*cos(a)+x(3)*sin(a))^2);   
V=(2+(x(3)*cos(a)-x(1)*sin(a))*G)/(1-x(1)^2-x(3)^2);
fx=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(1)*V-sin(a)*G-2*cos(a)*S);
fy=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(3)*V+cos(a)*G-2*sin(a)*S);
dxdt=zeros(4,1);
dxdt=[ x(2);
    -fx/M+p*sin(t);
     x(4);
     -fy/M+p*cos(t)+G1];

由于我总也是不会排版，每次上传程序就会变了模样，请大家多担待，顺便指导下怎么排版我的问题如下：
1.capone油膜力本身就是无量纲化的公式，后面对其转子动力学方程有进行了无量纲化，我编写的force11程序段是否正确？
2.后面对其画时间历程图、轴心轨迹、功率谱、庞加莱映射、相图，利用龙格库塔法解方程，我采用的是tspan=[0:T/100:14000*T],T=2*pi/w,是否合理？
3.对于非线性油膜力，是不是给定的初值不同，最后出来的结果也不相同？

非线性动力学

楼主，问题你解决了吗？我也在研究这个问题，可以交流下~

turkeyhazel

楼主你的问题解决了么 这个我也在做

zhangzy

如图所示，楼主点击那里，然后把代码贴进去  格式就好了

wjt1713573225

初值用rand 函数赋值 一般不会影响太大的 只是非线性油膜力注意无量纲化