声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1940|回复: 14

[分形与混沌] 看看我画的分岔图那里出错了?

[复制链接]
发表于 2007-9-13 16:12 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
感觉系统一直处于混沌状态,不应该是这样的,找不出原因所在
file:///C:/DOCUME~1/lxn/LOCALS~1/Temp/msohtml1/01/clip_image002.jpg
回复
分享到:

使用道具 举报

 楼主| 发表于 2007-9-13 16:30 | 显示全部楼层
[localimg=265,205]1[/localimg]
 楼主| 发表于 2007-9-13 16:31 | 显示全部楼层
888.bmp
发表于 2007-9-13 16:55 | 显示全部楼层

回复 #3 tianya7071 的帖子

时间步是不是太大了,你缩小试试看
发表于 2007-9-13 17:42 | 显示全部楼层

回复 #3 tianya7071 的帖子

很难猜测你的程序、及图的坐标。
发表于 2007-9-14 07:58 | 显示全部楼层
肯定不是混沌的!建议楼主将有关参数,如步长、误差精度进行调整再计算!
发表于 2007-9-14 08:41 | 显示全部楼层
没有程序,没有系统,什么都看不出来
 楼主| 发表于 2007-9-14 09:02 | 显示全部楼层
横坐标是转子的转速,纵坐标是响应的位移
画的是位移随转速变化时的分岔图
是不是步长在大了啊
我取的是周期的1/200,我觉得步长不大啊
发表于 2007-9-14 09:38 | 显示全部楼层
程序能不能发上来
发表于 2007-9-14 09:54 | 显示全部楼层

回复 #8 tianya7071 的帖子

感觉你的采样周期没有选好,是一周期计算200个点,那么你画图也应该是间隔200取一次点
 楼主| 发表于 2007-9-14 15:48 | 显示全部楼层

下面是我的程序代码

function dq = BallBrg_NonL_Sub_lyl(t,q)
%子程序

% 求解球轴承的变柔度(VC)振动
global w1 w2 Ri Re Nb gama kn  w2_rpm

n_One_T = 200;% 每个周期的采样点数
n_T = 300;% 采样时间占几个周期
%  滚子 轴承参数
Ri=0.0585;% 内滚道直径
Re=0.0665;% 外滚道直径
Nb=34;     % 滚子数
F = 3400;% 径向力(N)
q_initial(1:8,1) = 1e-11;% 初始值
gama = 0.00003;% 间隙(m)
kn =6.2e8;
w_cage = ( w2*Ri+w1*Re )/(Re+Ri);% ;% 保持架转速(rad/s)
M=[7.86 0 0 0;0 7.86 0 0;0 0 11.63 0;0 0 0 11.63];% 质量矩阵
    C=1200*[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1];% 阻尼矩阵
fq=zeros(2,1);% 轴承力初值
diff1_1_3 = q(1,1);% 内圈水平方向位移
diff1_2_4 = q(2,1);% 内圈垂直方向位移
diff12_1_3 = q(3,1);% 外圈水平方向位移
diff12_2_4 = q(4,1);% 外圈垂直方向位移
% 求轴承的非线性反力
for No_ball=1:Nb
    sita(No_ball) = 2*pi/Nb*(No_ball-1) + w_cage*t;% 第No_ball个滚珠的位置角
    Clearance(No_ball,1) = (diff1_1_3-diff12_1_3)*sin( sita(No_ball) ) ...
        + (diff1_2_4-diff12_2_4)*cos( sita(No_ball) ) - gama;% 滚珠与内滚道的间隙变化。
    % 判断哪几个滚动体受到接触力
    if Clearance(No_ball)<=0;
       Clearance(No_ball) = 0;
    end
    fs = abs( (Clearance(No_ball))^1.1 );
    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
Fq(1,1)= - fq(1,1); % 内圈水平方向外力
Fq(2,1)= - fq(2,1); ;% 内圈垂直方向外力
Fq(3,1)= fq(1,1); % 外圈水平方向外力
Fq(4,1)= fq(2,1)+F ;% 外圈垂直方向外力
% 动力学微分方程
dq(5:8,1)=inv(M)*(Fq-C*q(5:8,1));% x和y方向加速度
dq(1:4,1)=q(5:8,1);
function BallBrg_NonL_Forum
%主程序
% 求解轴承的变柔度(VC)振动
clear
clc
%% 参数设置
global w1 w2 Ri Re Nb gama kn  w2_rpm
n_One_T = 200;% 每个周期的采样点数
n_T = 300;% 采样时间占几个周期
%  滚子 球轴承参数
Ri=0.0585;% 内滚道直径
Re=0.0665;% 外滚道直径
Nb=34;     % 滚子数
F = 3400;% 径向力(N)
n_n = 0;
w2_limit1=4100;% 内圈最低转速(rpm)
w2_limit2=10600;% 内圈最高转速(rpm)
w_step = 100;% 转速变化步长(rpm)
q_initial(1:8,1) = 1e-11;% 初始值
gama = 0.00003;% 间隙(m)
kn =6.2e8;
% 滚子与滚道之间接触力与变形量的关系(N/mm^1.1)。
    w1 = 14400*pi/30;% 外圈角速度
    M=[7.86 0 0 0;0 7.86 0 0;0 0 11.63 0;0 0 0 11.63];% 质量矩阵
    C=1200*[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1];% 阻尼矩阵
    %% 响应计算循环
for  w2_rpm=w2_limit1:w_step:w2_limit2  %内圈转速变化
    n_n = n_n+1 % 计数变量
    disp(w2_rpm)
   w22 = w2_rpm*pi/30;% 转化为rad/s单位
    w2 = w22;% 内圈角速度
    w_cage = ( w2*Ri+w1*Re )/(Re+Ri);% 保持架
    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);% 时间数组
    [t,q]= ode45('BallBrg_NonL_Sub_lll', t_span, q_initial);   
plot(w2_rpm,q(50000:100:60000,1),'k.','markersize',3)
  hold on

[ 本帖最后由 tianya7071 于 2007-9-14 15:50 编辑 ]
发表于 2007-9-14 16:17 | 显示全部楼层
你不是一周期取200吗,怎么这里plot(w2_rpm,q(50000:100:60000,1)是100
发表于 2007-9-14 16:43 | 显示全部楼层
很不好意思的说,我不知道这个200怎么取得,没看到周期,呵呵
 楼主| 发表于 2007-9-17 09:42 | 显示全部楼层
我已经改过来了
我再运算一下看看
画分叉图一般取多少个周期的点画出来的效果最好啊
发表于 2007-9-17 10:06 | 显示全部楼层
多少个不好说,只要稳定了之后的差不多200就可以了,或者更多
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-14 14:27 , Processed in 0.090089 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表