sssssxxxxx921 发表于 2007-6-8 22:42

怎么在Matlab求解后画分岔图

这个是一个Matlab求解非线性方程的程序,不明白求解完后怎么画分岔图,还有庞加莱图,求助高手指点,小弟不胜感激

主程序

[ - ]CODE:

function BallBrg_NonL_Forum
% 求解外圈固定球轴承的变柔度(VC-Varying Compliance)振动(基于赵凌燕的论文)
% 程序有一些不合理、甚至错误的地方,可以用更好的代码代替,由于时间关系没有修改,
% 如有人感兴趣可以把修改的程序发布出来。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 作者:toes
% 版本:论坛发布版
% 相关程序:BallBrg_NonL_Sub_Forum
% 调试环境:Matlab7.0 WinXP SP2
% 参考文献:
% 1.赵凌燕.滚动轴承-转子系统的非线性动力学研究.西北工业大学硕士论文.2003.3.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc

%% 参数设置
% 用了全局变量来传递一些变量,不推荐,但是懒得改了,好心人优化一下。
global w d D Nb gama kn M C F

% 为了方便绘制分岔图而设置的参数
n_One_T = 100;% 每个周期的采样点数
n_T = 100;% 采样时间占几个周期

% 61903/P5(17*30*7) 球轴承参数
d=0.0173;% 内滚道直径
D=0.0265;% 外滚道直径
Nb=9; % 滚子数

n_n = 0;
w_limit1=100;% 最低转速(rpm)
w_limit2=20000;% 最高转速(rpm)
w_step = 100;% 转速变化步长(rpm)

q_initial(1:4,1) = 1e-11;% 初始值
gama = 0.00002;% 间隙(m)
F = 6;% 径向力(N)
kn = 7.055e9*0.001^1.5;
% 滚子与滚道之间接触力与变形量的关系(N/mm^1.5)。赵的论文给出。
M=0.6*;% 质量矩阵
C=200*;% 阻尼矩阵

%% 响应计算循环
for w_rpm=w_limit1:w_step:w_limit2

n_n = n_n+1 % 计数变量
disp(w_rpm)
w = w_rpm*pi/30;% 转化为rad/s单位

wi = w;% 内圈角速度
wo = 0;% 外圈角速度

w_cage = ( wi*d/2+wo*D/2 )/2/((D+d)/4);% 保持架
w_vc = w_cage*Nb/2/pi; % 变刚度频率(vc频率)。单位Hz
T_vc = 1/w_vc;% vc周期

dt=T_vc/n_One_T;% 取点时间步长,dt随转速变化。
time=n_T*T_vc;% 总的时间

n = round(time/dt);% 离散点数
t_span(1:n) = linspace(0,time,n);% 时间数组

= ode23tb('BallBrg_NonL_Sub_Forum', t_span, q_initial);
% 至于用什么ode函数求解合适需要比较验证

Q{n_n}=q;
save Q.mat Q; % 存储数据
end
disp('Calculation is done!')

子程序

[ - ]CODE:

function dq = BallBrg_NonL_Sub_Forum(t,q)
% BallBrg_NonL调用的微分方程子程序
% 求解外圈固定球轴承的变柔度(VC)振动
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 作者:toes
% 版本:论坛发布版
% 相关程序:BallBrg_NonL_Forum
% 参考文献:
% 1.赵凌燕.滚动轴承-转子系统的非线性动力学研究.西北工业大学硕士论文.2003.3.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global w d D Nb gama kn M C F

wi = w;
wo = 0;
w_cage=( wi*d+wo*D )/4/((D+d)/4);% 保持架转速(rad/s)

fq=zeros(2,1);% 轴承力初值

diff_1_3 = q(1,1);% 水平方向位移
diff_2_4 = q(2,1);% 垂直方向位移

% 求轴承的非线性反力
for No_ball=1:Nb
sita(No_ball) = 2*pi/Nb*(No_ball-1) + w_cage*t;% 第No_ball个滚珠的位置角
Clearance(No_ball,1) = diff_1_3*sin( sita(No_ball) ) ...
+ diff_2_4*cos( sita(No_ball) ) - gama;% 滚珠与内滚道的间隙变化。
% 判断哪几个滚动体受到接触力
if Clearance(No_ball)<=0;
Clearance(No_ball) = 0;
end
fs = abs( (1000*Clearance(No_ball))^1.5 );

fq(1,1) = fq(1,1)+kn*fs*sin(sita(No_ball));
fq(2,1) = fq(2,1)+kn*fs*cos(sita(No_ball));
end

F_m1d1_cos = 0;% 不平衡力在水平方向的投影。本例不考虑。
F_m1d1_sin = 0;% 不平衡力在垂直方向的投影。本例不考虑。

Fq(1,1)= - fq(1,1) + F_m1d1_cos;% 水平方向外力
Fq(2,1)= - fq(2,1) + F_m1d1_sin - F;% 垂直方向外力

K = ;% 刚性转子,轴段为刚性。

% 动力学微分方程
dq(3:4,1)=inv(M)*(Fq-K*q(1:2,1)-C*q(3:4,1));% x和y方向加速度
dq(1:2,1)=q(3:4,1);

[ 本帖最后由 eight 于 2007-6-8 23:18 编辑 ]

xjzuo 发表于 2007-6-9 07:59

确定控制参数,作循环即可.

sssssxxxxx921 发表于 2007-6-9 15:00

回复 #2 xjzuo 的帖子

我在主程序循环后用wc=[];wc=作为控制参数;
最后plot(wc,q(0:T_vc:time,1))
也总是出错不知道怎么回事请指点

happy 发表于 2007-6-9 16:44

针对不同的wc分别求解微分方程,得到相应的时间序列

在特定的wc下,每周期取一个点画图,即可得到庞加莱图

以wc为横坐标,对应的将各wc下庞加莱图时用的点画上就是分岔图了

详细的例子我记得一般力学版有好多个例子

sssssxxxxx921 发表于 2007-6-10 15:05

回复 #4 happy 的帖子

你的一意思时不时说   先将点映射到庞加来截面上 ,然后在画分岔图啊?

直接用wc做横坐标,q做纵坐标画出来的图不是分岔图吗?

庞加来图到底怎么画啊,对于我这个问题,它的截面选取是怎么一回事啊

你能帮我画一个做个示范吗

谢谢了

happy 发表于 2007-6-11 02:11

poincare在帖子http://forum.vibunion.com/forum/thread-12312-1-1.html一贴中已经讨论的挺清楚地了
建议好好看看,里边也有例子可以参考

happy 发表于 2007-6-11 02:19

原帖由 sssssxxxxx921 于 2007-6-10 15:05 发表 http://www.chinavib.com/forum/images/common/back.gif
直接用wc做横坐标,q做纵坐标画出来的图不是分岔图吗?

这个显然不是分岔图

sssssxxxxx921 发表于 2007-6-11 09:39

回复 #7 happy 的帖子

我又仔细看了你说的,你的意思和我用直接用wc做横坐标,q做纵坐标画分岔图,应该一样啊

你的意思是:先用转速得到频率,也就是时间序列,
               
               然后在某一转速(在特定的wc下),在求出的q中每周期取一个点画图,即可得到庞加莱图
                  
                  以wc为横坐标,对应的将各wc下庞加莱图时用的点画上就是分岔图

是不是这意思啊

sssssxxxxx921 发表于 2007-6-11 09:42

回复 #7 happy 的帖子

如果我说的不对,你能不能帮我看看这个程序,给我画个庞加来图和分岔图,

我想对着图和程序自己琢磨琢磨比较好,也省的这么麻烦你

谢谢了

appleseed05 发表于 2007-6-13 14:28

简单的说:
所谓分岔图就是隔周期选取一个响应的数据点,然后用想要研究的参数做为x轴,一般用转速比较多,在同一转速下隔周期取得的所有位移结果做为y轴画图。
所谓庞加莱图,一般就是将位移做为x轴,同一时刻速度作为y轴画图。
页: [1]
查看完整版本: 怎么在Matlab求解后画分岔图