二。程序实现
（本来与第一部分写在一起，后来发现太长了，可能很少有人耐心看完，所以将这一部分分开发布。）
理论基础部分只是阐述了非线性轴承力的产生和推导，实际计算的时候是考虑了一个转子-轴承系统中轴承非线性力对整个系统的影响。
假设转子为刚性，两端由球轴承支撑，盘在轴段的中间。转子轴承系统的动力学方程如下：
05-2
下面以西北工业大学赵凌燕论文《滚动轴承－转子系统的非线性动力学研究》上的一组数据为例编写Matlab程序，其具体数据见程序，不再一一列出。
程序主要是求解一个2元2阶的常微分方程组，数值积分得到系统的响应。其中，为了绘制分岔图，选择响应取点的周期为激励力的周期（这里即变刚度激励力的周期1/fvc）。
主程序
- 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*[1 0;0 1];% 质量矩阵
- C=200*[1 0; 0 1];% 阻尼矩阵
- %% 响应计算循环
- 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);% 时间数组
-
- [t,q]= ode23tb('BallBrg_NonL_Sub_Forum', t_span, q_initial);
- % 至于用什么ode函数求解合适需要比较验证
- Q{n_n}=q;
- save Q.mat Q; % 存储数据
- end
- disp('Calculation is done!')
子程序
- 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 = [0 0; 0 0];% 刚性转子，轴段为刚性。
- % 动力学微分方程
- 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);
（其余部分整理中，待续。）
[ 本帖最后由 toes 于 2006-9-5 13:32 编辑 ] |
发表于 2006-9-3 18:00
提升卡
变色卡
显身卡
楼主
发表于 2006-9-6 21:47
