toes 发表于 2006-9-3 18:00

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

球轴承的变刚度(或变柔度)振动(varying compliance vibration)在90年代曾经是一个热点,2000年后也有人继续研究。
因为前一阵正在做非线性方面的东西,所以以别人的论文为基础,将求解球轴承的变刚度振动问题作为训练手段,最后与他人的结果做了简单对比,这样可以验证相关程序并熟悉求解过程。
这几天闲来无事,^_^,将理论部分、计算程序、简单结果逐步上传,希望大家参与讨论给予意见!

一。理论基础

更正:附件图01.jpg中第一句 “如图1.2所示” 应为 “如图1.1所示”。





(其它部分将以回帖的形式陆续更新,待续)

[ 本帖最后由 yejet 于 2006-10-3 17:01 编辑 ]

toes 发表于 2006-9-5 13:25

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

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

下面以西北工业大学赵凌燕论文《滚动轴承-转子系统的非线性动力学研究》上的一组数据为例编写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*;% 质量矩阵
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!')


子程序


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);

(其余部分整理中,待续。)

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

toes 发表于 2006-9-6 18:54

三。结果分析

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

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

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

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

分岔图并不是很清楚,很多转速段可以细致再做(主要是我懒,呵呵),但从图中可以看到几处明显是周期运动或非周期运动的转速区域。比如,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

回上一楼问

他是拿别人估计是他师姐的东西重做了,我可以告诉你他没考虑轴承接触角变化
页: [1] 2
查看完整版本: 【原创】球轴承的变刚度振动