滑动轴承做时间历程、轴心轨迹、相图、庞加莱映射程序
本帖最后由 dj15175771160 于 2014-3-17 20:00 编辑我以滑动轴承为例做时间历程、轴心轨迹、相图、庞加莱映射图，程序如下，但是变换转速和偏心率，做出来的图总感觉不对，希望大家运行下看看你，滑动轴承采用的是capone油膜力，对建立的转子动力学方程进行了无量纲化。
functionjeffcottzhuanzixitong_solution
clear all
clc
n=10000;
w=pi*n/30;
%w=;
%x0=;
x0=;
options=odeset('RelTol',1.e-7,'AbsTol', );
% %for j=1:length(w)
% T=2*pi; %一个周期
% =ode45('force11',,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;
% =ode45('force11',,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;
=ode45('force11',,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;%周期
% =ode45('force11',,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;
% =ode45('force11',,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 ; %%寻找两个最大的频率
% =max(power1);
% power1(daa1)=0;
% =max(power1);
% WF=2*pi*;%求圆频率
%
% figure %画功率谱图
% plot(2*pi*freq(1:100),power(1:100)) %%画功率谱图
% %hold on
% title('功率谱')
% xlabel('\omega/s^{-1}');ylabel('p/m^2');
% % plot(WF,,'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=,T=2*pi/w,是否合理？
3.对于非线性油膜力，是不是给定的初值不同，最后出来的结果也不相同？
楼主，问题你解决了吗？我也在研究这个问题，可以交流下~ 楼主你的问题解决了么 这个我也在做 如图所示，楼主点击那里，然后把代码贴进去格式就好了
初值用rand 函数赋值 一般不会影响太大的 只是非线性油膜力注意无量纲化
页:
[1]