【原创】球轴承的变刚度振动

引用 toes 2006-9-5 13:25

二。程序实现
（本来与第一部分写在一起，后来发现太长了，可能很少有人耐心看完，所以将这一部分分开发布。）

理论基础部分只是阐述了非线性轴承力的产生和推导，实际计算的时候是考虑了一个转子-轴承系统中轴承非线性力对整个系统的影响。
假设转子为刚性，两端由球轴承支撑，盘在轴段的中间。转子轴承系统的动力学方程如下：

05-2

下面以西北工业大学赵凌燕论文《滚动轴承－转子系统的非线性动力学研究》上的一组数据为例编写Matlab程序，其具体数据见程序，不再一一列出。
程序主要是求解一个2元2阶的常微分方程组，数值积分得到系统的响应。其中，为了绘制分岔图，选择响应取点的周期为激励力的周期（这里即变刚度激励力的周期1/fvc）。

主程序

1. 
   
2. function BallBrg_NonL_Forum
   
3. % 求解外圈固定球轴承的变柔度(VC-Varying Compliance)振动（基于赵凌燕的论文）
   
4. % 程序有一些不合理、甚至错误的地方，可以用更好的代码代替，由于时间关系没有修改，
   
5. % 如有人感兴趣可以把修改的程序发布出来。
   
6. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
7. % 作者：toes
   
8. % 版本：论坛发布版
   
9. % 相关程序：BallBrg_NonL_Sub_Forum
   
10. % 调试环境：Matlab7.0   WinXP SP2
    
11. % 参考文献：
    
12. % 1.赵凌燕.滚动轴承－转子系统的非线性动力学研究.西北工业大学硕士论文.2003.3.
    
13. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
14. clear
    
15. clc
    
16. 
    
17. %% 参数设置
    
18. % 用了全局变量来传递一些变量，不推荐，但是懒得改了，好心人优化一下。
    
19. global w d D Nb gama kn M C F
    
20. 
    
21. % 为了方便绘制分岔图而设置的参数
    
22. n_One_T = 100;% 每个周期的采样点数
    
23. n_T = 100;% 采样时间占几个周期
    
24. 
    
25. % 61903/P5(17*30*7) 球轴承参数
    
26. d=0.0173;% 内滚道直径
    
27. D=0.0265;% 外滚道直径
    
28. Nb=9;     % 滚子数
    
29. 
    
30. n_n = 0;
    
31. w_limit1=100;% 最低转速（rpm）
    
32. w_limit2=20000;% 最高转速（rpm）
    
33. w_step = 100;% 转速变化步长(rpm)
    
34. 
    
35. q_initial(1:4,1) = 1e-11;% 初始值
    
36. gama = 0.00002;% 间隙(m)
    
37. F = 6;% 径向力(N)
    
38. kn = 7.055e9*0.001^1.5;
    
39. % 滚子与滚道之间接触力与变形量的关系(N/mm^1.5)。赵的论文给出。
    
40. M=0.6*[1 0;0 1];% 质量矩阵
    
41. C=200*[1 0; 0 1];% 阻尼矩阵
    
42. 
    
43. %% 响应计算循环
    
44. for w_rpm=w_limit1:w_step:w_limit2
    
45. 
    
46.     n_n = n_n+1 % 计数变量
    
47.     disp(w_rpm)
    
48.     w = w_rpm*pi/30;% 转化为rad/s单位
    
49.    
    
50.     wi = w;% 内圈角速度
    
51.     wo = 0;% 外圈角速度
    
52.    
    
53.     w_cage = ( wi*d/2+wo*D/2 )/2/((D+d)/4);% 保持架
    
54.     w_vc = w_cage*Nb/2/pi; % 变刚度频率(vc频率)。单位Hz
    
55.     T_vc = 1/w_vc;% vc周期
    
56. 
    
57.     dt=T_vc/n_One_T;% 取点时间步长，dt随转速变化。
    
58.     time=n_T*T_vc;% 总的时间
    
59.    
    
60.     n = round(time/dt);% 离散点数
    
61.     t_span(1:n) = linspace(0,time,n);% 时间数组
    
62.    
    
63.     [t,q]= ode23tb('BallBrg_NonL_Sub_Forum', t_span, q_initial);
    
64.     % 至于用什么ode函数求解合适需要比较验证
    
65. 
    
66.     Q{n_n}=q;
    
67.     save Q.mat Q; % 存储数据
    
68. end
    
69. disp('Calculation is done!')
    

复制代码

子程序

1. 
   
2. function dq = BallBrg_NonL_Sub_Forum(t,q)
   
3. % BallBrg_NonL调用的微分方程子程序
   
4. % 求解外圈固定球轴承的变柔度(VC)振动
   
5. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
6. % 作者：toes
   
7. % 版本：论坛发布版
   
8. % 相关程序：BallBrg_NonL_Forum
   
9. % 参考文献：
   
10. % 1.赵凌燕.滚动轴承－转子系统的非线性动力学研究.西北工业大学硕士论文.2003.3.
    
11. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
12. global w d D Nb gama kn M C F
    
13. 
    
14. wi = w;
    
15. wo = 0;
    
16. w_cage=( wi*d+wo*D )/4/((D+d)/4);% 保持架转速（rad/s）
    
17. 
    
18. fq=zeros(2,1);% 轴承力初值
    
19. 
    
20. diff_1_3 = q(1,1);% 水平方向位移
    
21. diff_2_4 = q(2,1);% 垂直方向位移
    
22. 
    
23. % 求轴承的非线性反力
    
24. for No_ball=1:Nb
    
25.     sita(No_ball) = 2*pi/Nb*(No_ball-1) + w_cage*t;% 第No_ball个滚珠的位置角
    
26.     Clearance(No_ball,1) = diff_1_3*sin( sita(No_ball) ) ...
    
27.         + diff_2_4*cos( sita(No_ball) ) - gama;% 滚珠与内滚道的间隙变化。
    
28.     % 判断哪几个滚动体受到接触力
    
29.     if Clearance(No_ball)<=0;
    
30.        Clearance(No_ball) = 0;
    
31.     end
    
32.     fs = abs( (1000*Clearance(No_ball))^1.5 );
    
33. 
    
34.     fq(1,1) = fq(1,1)+kn*fs*sin(sita(No_ball));
    
35.     fq(2,1) = fq(2,1)+kn*fs*cos(sita(No_ball));
    
36. end
    
37. 
    
38. F_m1d1_cos = 0;% 不平衡力在水平方向的投影。本例不考虑。
    
39. F_m1d1_sin = 0;% 不平衡力在垂直方向的投影。本例不考虑。
    
40. 
    
41. Fq(1,1)= - fq(1,1) + F_m1d1_cos;% 水平方向外力
    
42. Fq(2,1)= - fq(2,1) + F_m1d1_sin - F;% 垂直方向外力
    
43. 
    
44. K = [0 0; 0 0];% 刚性转子，轴段为刚性。
    
45. 
    
46. % 动力学微分方程
    
47. dq(3:4,1)=inv(M)*(Fq-K*q(1:2,1)-C*q(3:4,1));% x和y方向加速度
    
48. dq(1:2,1)=q(3:4,1);
    

复制代码

（其余部分整理中，待续。）

[ 本帖最后由 toes 于 2006-9-5 13:32 编辑 ]

引用 toes 2006-9-6 18:54

三。结果分析

应用第二部分的MATLAB代码求解，可以得到转子－轴承系统的非线性响应。
我们可以对非线性响应进行分析，来研究响应的特点，也能由此判断系统的运行状况。一般，分析会从时域波形、频谱图、庞加莱图或分岔图等方法开始着手。这里，就频谱图和分岔图各举一例。

1。下图为转子－轴承系统10000rpm时的频谱图：

FFT01

可以看到，图中主要以变刚度频率fvc和它的次谐波频率为主，包括fvc、1/5fvc等等；而且图中杂乱的频率很多，幅值也较大，就像有“噪声”干扰一样，这说明系统已不是在做周期运动了，而是一种非周期的状态。

2。再从分岔图来分析。下图为系统100～20000rpm的分岔图：

bifurcation02

分岔图并不是很清楚，很多转速段可以细致再做（主要是我懒，呵呵），但从图中可以看到几处明显是周期运动或非周期运动的转速区域。比如，12000～15000rpm之间系统处于周期1的运动状态，而在频谱图中分析的10000rpm左右为非周期运动（具体是混沌还是准周期还可从庞加莱图判断，但从其附近的状态来讲，应该是混沌，实际上其频谱图也证明了这一点）。

3。需要补充说明的是：
这个算例只是为了考察轴承变刚度对转子系统动力特性的影响，所以其响应中的主要频率成分是变刚度频率fvc。对于工程实际中的转子系统，转子的不平衡量是无法去除的，在转速较高的情况下，响应中的频率成分会以不平衡力频率为主，所以变刚度的影响相对会很小，fvc不容易被试验测到。另外在存在不平衡量的情况下，应该用多频频闪法做庞加莱图和分岔图，虽然我认为多频法在轴承的变刚度分析中并很不实用，呵呵。

（全文完，盼望有人扔西红柿，谢绝砖头。）

[ 本帖最后由 toes 于 2006-9-6 19:10 编辑 ]

引用 hongking111 2006-9-6 21:34

动力学方程中的系数Kn是定常的吗？变刚度体现在哪？
还有这个系统有解析解吗？

引用 yejet 2006-9-6 21:47

原帖由 hongking111 于 2006-9-6 21:34 发表
动力学方程中的系数Kn是定常的吗？变刚度体现在哪？
还有这个系统有解析解吗？

这里的Kn应该是个常数，变刚度应该体现在公式中的几个三角函数上

这个系统解析解估计是不太可能得到，不过近似解倒是有可能

上面仅是个人看法，希望楼主给出确切答案

引用 toes 2006-9-6 23:40

原帖由 hongking111 于 2006-9-6 21:34 发表
动力学方程中的系数Kn是定常的吗？变刚度体现在哪？
还有这个系统有解析解吗？

kn为常数。
变刚度体现在公式1-10。由这个公式得到的轴承当量刚度是随时间发生变化的。
好像没有看到有人求得它的解析解和近似解。

引用 NASA 2006-9-7 21:15

原帖由 toes 于 2006-9-6 23:40 发表


kn为常数。
变刚度体现在公式1-10。由这个公式得到的轴承当量刚度是随时间发生变化的。
好像没有看到有人求得它的解析解和近似解。

如果要考虑近似解的话大家不知道有什么建议？比如用什么方法比较合适？

引用 toes 2006-9-8 15:41

估计有难度，因为这个模型中的非线性因素至少有两个：
1。单个滚珠与单个滚道之间力与变形的关系是非线性的，一般是3/2次方。
2。就算滚珠与滚道之间的刚度是常数，随着轴承的运转，整个轴承的当量刚度也是非线性的。

所以两个非线性因素叠加在一起的话，求近似解是有相当困难的。

引用 glise 2006-9-9 19:58

原帖由 toes 于 2006-9-8 15:41 发表
估计有难度，因为这个模型中的非线性因素至少有两个：
1。单个滚珠与单个滚道之间力与变形的关系是非线性的，一般是3/2次方。
2。就算滚珠与滚道之间的刚度是常数，随着轴承的运转，整个轴承的当量刚度也是非线 ...

当量刚度指的是哪一部分？

引用 toes 2006-9-9 20:18

指的是轴承径向载荷与径向位移之间的关系，由二者的比值可以得到一个当量刚度。

引用 laosy 2006-9-11 16:15

精细地研究这个问题应该相当复杂！
如果转速较高，并且载荷较大的时候。
楼主用的Hertz定接触并不适用，此时接触刚度是动态变化的 ，需要考虑接触碰撞过程的能量耗散！

引用 yeziliu 2006-9-11 17:30

知道质量，刚度，阻尼比，怎么计算阻尼？

引用 vicky_boy 2006-12-8 10:50

你建立动力学方程时，有没有考虑接触角的变化？你的接触角是否为常量
我也在做球轴承的非线性振动的研究，欢迎大家经常交流，加qq：23126773，验证码：轴承

引用 sssxxx921 2006-12-11 09:26

什么叫“多频频闪法” ，怎么没听说过，查百度也只你这有，你是西工大袁如的学生吧？

引用 sssxxx921 2006-12-11 09:30

他是拿别人估计是他师姐的东西重做了，我可以告诉你他没考虑轴承接触角变化

引用 sssssxxxxx921 2007-6-6 21:06

你这个最后怎么画图呀  麻烦把画图那部分的程序也发上来吧 :handshake

引用 tianya7071 2007-9-17 15:42

楼主请你的上面两个图的程序代码也发上来吧
邮箱tianyarrrbbb@com

引用 toes 2010-1-11 17:22

本帖最后由 VibInfo 于 2016-5-6 14:58 编辑

原帖由 sssxxx921 于 2006-12-11 09:30 发表
他是拿别人估计是他师姐的东西重做了，我可以告诉你他没考虑轴承接触角变化

很久没来了。
是拿别人的东西重做了，主要是赵凌燕的（其实还有其他人的论文），我在前面正文和程序部分都说明了，看了她的论文理论部分编写的程序，就是自己练练手。
没有考虑接触角变化。
她不是我师姐，呵呵。

[ 本帖最后由 toes 于 2010-1-11 17:50 编辑 ]

引用 bingkafei7155 2011-11-12 00:45

群主你好，我也是研究滚动轴承的，能否将你最后两图的代码发给我啊？不知你还保存着没啊？万分感谢！837043771@qq.com

引用 冬夜 2012-4-10 17:35

还有没有研究这一方面的研究，大家一起交流：2582968900

查看全部评论(29)