基于fft的深沟球轴承和角接触球轴承弹性变形快速计算方法
【专利摘要】本发明涉及一种基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算方法，属于滚动轴承优化设计技术领域。该方法基于深沟球轴承和角接触球轴承的点接触应力分布和影响系数矩阵，将离散循环卷积定理和FFT技术应用于轴承弹性变形计算中，计算各滚动体与内圈或外圈之间的弹性变形量及各滚动体与内外圈总的弹性变形量。本发明提供的方法能够有效解决传统方法计算效率不高的问题，因为使用传统方法计算深沟球轴承和角接触球轴承弹性变形时需要进行复杂的四重循环计算，严重耗费计算机内存且计算量很大，而本方法只要正确运用离散循环卷积定理和FFT技术，便能准确、快速地获得计算区域内的轴承弹性变形。
【专利说明】
基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算 方法
技术领域
[0001]本发明属于滚动轴承优化设计技术领域，涉及一种基于FFT(Fast Fourier Transform:快速傅里叶变换)的深沟球轴承和角接触球轴承弹性变形快速计算方法。
【背景技术】
[0002] 作为重要的动力传递部件，深沟球轴承和角接触球轴承在工程机械、轨道交通、航 空航天及新能源等领域有着广泛的应用。深沟球轴承和角接触球轴承各滚动体与内圈或外 圈之间弹性变形的快速求解有助于轴承动态特性（刚度、阻尼）分析和油膜厚度计算。在经 典的轴承分析中普遍采用刚性套圈假设，即假定轴承套圈是不发生弹性变形的刚体，弹性 变形只发生在滚动体上，滚动体与套圈的接触只会产生局部的接触变形而不会改变套圈的 整体形状和尺寸。这种处理方法虽然为套圈的位移和轴承弹性变形分析带来了很大的方 便，但是它无法准确反映出轴承内外套圈的变形量，从而无法准确预测轴承的油膜刚度、接 触刚度和阻尼，也无法用于分析考虑轴承表面加工精度时的轴承弹性变形。此外，使用传统 的Boussinesq方法计算深沟球轴承和角接触球轴承弹性变形时需要进行复杂的四重循环 运算，严重耗费计算机内存且计算量很大，尤其是当网格密度较大时，计算量将呈指数级增 长。
[0003] 所以，传统的滚动轴承弹性变形计算方法更多的是具有统计与经验意义，其精度 上不够准确、效率上不够高效。
【发明内容】
[0004] 有鉴于此，本发明的目的在于提供一种基于FFT的深沟球轴承和角接触球轴承弹 性变形快速计算方法，该方法能够解决传统方法中存在的精度差、效率低的问题。
[0005] 为达到上述目的，本发明提供如下技术方案：
[0006] -种基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算方法，其特征在于： 包括以下步骤：
[0007] S1:给定/修正轴承载荷分布范围参数e ;
[0008] S2:查询载荷分布积分表并通过线性插值得到载荷分布积分Jr( e);
[0009] S3:由式Sr = Pd/[2(l-2e)]计算得到内、外圈在径向载荷方向上的相对位移Sr，其 中Pd表示径向游隙；
[0010] S4:由式?:。= 21(11(心-?(1/2)1勺:（￡)计算得到计算外载荷?:。，其中2表示滚动体数 目，1表示滚动体与内、外圈之间总的负荷-变形常数；
[0011] S5:判断计算外载荷Fr。与给定的径向载荷Fr之间的关系是否满足收敛条件，若 (Frc_Fr)/Fr | > |，则按照式enew = eold(Frc/Fr) -e修正e并重复步骤S1 -步骤S4，其中，|一般取 0 ? 00 1，在0 ? 00 1与0 ? 0 1之间取任意值；若| (Frc-Fr ) / Fr |彡|，则依据式
计算各滚球的载荷分布么^…，々，其中七表示各滚动体的方位 角，i = l，???，Z，执行下一步；
[0012] S6:依据各滚球的载荷分布，运用Hertz接触理论，计算接触椭圆长、短半轴长a、b;
[0013] S7:计算各滚球在内、外圈接触区上的接触应力分布p，将接触应力分布函数p(x, y)离散为(M+l) X (N+1)的矩阵p，然后补零使其成为2MX 2N的应力矩阵p' ；
[0014] S8:以坐标原点为中心，在4倍的计算区域内求得影响系数矩阵C，维数为2MX2N;
[0015] S9:对矩阵p'和C分别进行FFT变换，得到其FFT变换矩阵r和^
[0016] S10: F和e的对应项相乘得到频域内的弹性变形矩阵6,对?作IFFT变换（Inverse Fast Fourier Transform:快速傅里叶逆变换)得到时域上的弹性变形矩阵D;
[0017] 311:令￥^ = 0[^,_(1=0，1，...凡」=0，1，...，?，获得待求计算域内的弹性变 形V^，分别在各滚球与内、外圈接触区上进行步骤S6-步骤S11，得到各滚球与内、外圈间的 弹性变形分布。
[0018]进一步，在步骤S8中，影响系数矩阵C的求解采用Green函数法
[0019] 本发明的有益效果在于:本方法能用于同时计算滚动体的弹性变形和内、外圈的 弹性变形，也能用于考察滚动体和内、外滚道加工精度对轴承弹性变形的影响，能够使计算 结果更加符合实际情况，也使得轴承刚度、阻尼和油膜厚度的计算结果更加准确，避免了传 统Boussinesq方法中四重循环的运算，尤其是网格密度较大时，能极大程度地提高计算效 率，本发明提出的方法对轴承的优化设计和轴承系统的受力分析是有益的。
【附图说明】
[0020] 为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行 说明：
[0021 ]图1为本发明的迭代计算方法流程图；
[0022] 图2为深沟球轴承和角接触球轴承的滚球与内圈接触椭圆示意图；
[0023] 图3为影响系数矩阵求解说明示意图；
[0024] 图4为影响系数矩阵的计算区域示意图。
【具体实施方式】
[0025]下面将结合附图，对本发明的优选实施例进行详细的描述。
[0026] 图1为本发明的迭代计算方法流程图，对于给定了径向游隙Pd和径向载荷Fr的深沟 球轴承和角接触球轴承，能够用图1流程图右半部分所示的控制误差的试解法求解各滚球 的载荷分布：
[0027] 首先假定轴承载荷分布范围参数e的值，查询载荷分布积分表并通过线性插值得 到载荷分布积分Jr(e)。然后由式S r = Pd/[2(l_2e)]计算得到内、外圈在径向载荷方向上的 相对位移Sr，由式F rc = ZKn( Sr-Pd/2)h5Jr(e )计算得到计算外载荷Frc，其中Z表示滚动体数 目，1表示滚动体与内、外圈之间总的负荷-变形常数。最后，判断计算外载荷Fr。是否满足 (F rc-Fr)/Fr|彡若不满足，则按照式enew=ecild(F rc/Fr)-e修正e并重复上述过程，其中彳一 般取0 . 0 0 1，0在0 . 0 0 1与0 . 0 1之间取任意值；若满足收敛条件，则依据式
计算各滚球的载荷分布^^其中，也i表示各滚动体的方位角（i = 1，…，Z) 〇
[0028]由于考虑的是中低速轴承，故没有考虑离心力和陀螺力矩的影响，各滚球与内、外 圈的接触角相同，以上计算出来的滚球载荷既表示滚球与内圈间的接触载荷也表示滚球与 外圈间的接触载荷。此载荷为接触点上的单值，后续接触应力分布P的计算需用到此值。 [00 29]根据Hertz接触理论，接触应力在接触区内按椭球体规律分布。当接触椭圆的长轴 方向与x轴相重合时，接触应力分布p为：
⑴
[0031]图2为深沟球轴承和角接触球轴承的滚球与内圈接触椭圆示意图，其中，如图2，a、 b分别表示接触椭圆的长、短半轴长，采用下列公式能够确定：
(2)
[0034]其中，ka、kb为接触椭圆尺寸系数，E'表示两接触体的综合弹性模量，常数A、B与两 接触体的几何形状有关，其值由下列两方程确定：
[0037]如图2，以滚球和内滚道接触为例，上两方程中的变量Rlx、Rly、R2x、R 2y下标中的第1 个数字1和2分别表示滚球和内滚道，下标中的第2个字母x和y分别表示轴向和周向，也即： Rlx表示滚球的轴向曲率半径，Rly表示滚球的周向曲率半径，R2x表示内滚道的轴向曲率半 径，R 2y表示内滚道的周向曲率半径。
[0038] 求出常数A、B之后，根据式C〇S0 = (B-A)/(A+B)计算0值，然后根据0值查询"接触椭 圆尺寸系数表"能得到接触椭圆尺寸系数ka、k b的值，对于表中未列的0值，插值可得到。
[0039] 至此，每个滚球与内、外圈间的接触应力分布p均能依据上述方法得到。
[0040] 深沟球轴承和角接触球轴承的弹性变形求解问题可等价为求解半无限体表面在 上述求出的接触应力分布P下的弹性变形求解问题。其弹性变形方程可表达为：
(6)
[0042]使用传统的Boussinesq方法求解上述弹性变形方程时，需要将式(6)离散化，这种 计算存在两个困难:首先，式(6)右端是奇异积分，奇点为(x，y) = U，n);其次，为求某一点 的弹性变形，就需要对整个求解域计算积分(需两重循环），而为了求解整个计算域的弹性 变形，又需用两重循环，故每一次弹性变形计算都需要进行四重循环运算。这严重耗费计算 机内存且计算量很大，尤其是当网格密度较大时，计算量将呈指数级增长。例如，当网格密 度M X N为32 X 32时，循环需进行1048576次，当网格密度M X N仅增大1倍为64 X 64时，循环次 数便迅速增加到16777216次。可见，使用一般的数值积分方法，计算工作量太大，而且误差 的频繁累积导致难以达到要求的精度。
[0043]在本发明方法中，基于线弹性叠加原理，为方便采用数值方法求解，式(6)可离散 为：
(7)
[0045] 其中，C(xi-xk，yj-yi)为影响系数，简记为CJj，定义为在节点(k，1)上作用的单位力 在节点（i，j)上产生的弹性变形，如图3所示。从式(7)能发现弹性变形的计算分为两步:影 响系数的计算和卷积和的计算：
[0046] 第一步:弹性变形影响系数的计算
[0047] 根据定义，点接触弹性变形影响系数可按下式计算：
[0048] C；j -<5, v,(8)
[0049] 其中，A Q为(Xk，yi)点处的单元网格，如图3所示阴影部分。wU，n)为插值函数，用 来近似单元网格上的接触应力分布，g ( x，y )称为响应函数（亦称G r e e n函数）， 当采用等间距网格时，頌只和两点之间的距离有关且存在对称性，因此影 响系数的计算可只针对坐标原点进行，可进一步简记为Qu，故式(8)可简化为：
[0050] Ck, = J5；lt(9)
[0051] 这里XhhS力作用点和待求变形点间距离的坐标分量，如图3所示。
[0052] 对于式(9)，可用Green函数法求解：
[0053]插值函数w(|，n)在单元网格内取1，同时假设Green函数在单元格内不变化，取节 点值h (xi，y j)作为其近似值，则影响系数可简写为：
[0054] Ci, j = gi, j ? A x A y (10)
[0055] Green函数在原点存在一个奇异点，可通过以原点为中心的单元网格上的计算均 值来代替原点处的离散值，即当I Xi |彡2 A x且|yj |彡2 A y时，以式（11)计算值来代替点(xi, yj)处的离散值：
{〇)
[0057]当给定目标区域(计算区域)后，影响系数应在4倍的目标区域内计算，即在如图4 所示的扩展区域内计算。当采用等间距网格时，只需计算出扩展区域内各节点相对于坐标 原点的影响系数。根据对称性，只需计算一个象限内的影响系数便能确定其他象限内的影 响系数。
[0058]第二步:弹性变形卷积和的计算
[0059]直接计算卷积和存在大量的重复运算，而运用离散循环卷积定理和FFT是一种高 效的较为理想的计算方法。在计算出各滚球在内、外圈接触区上的接触应力分布P和影响系 数矩阵C之后，基于离散循环卷积定理和FFT技术的深沟球轴承和角接触球轴承弹性变形计 算如下：
[0060] (1)将接触应力分布函数p(x，y)离散为(M+1)X(N+1)的矩阵P，然后补零使其成为 2MX2N的应力矩阵p' ；
[0061] (2)以坐标原点为中心，在4倍的计算区域内求得影响系数矩阵C，维数为2MX2N;
[0062] (3)对矩阵p'和C分别进行FFT变换，得到其FFT变换矩阵F和亡；
[0063] (4)#和的对应项相乘得到频域内的变形矩阵D，对i作IFFT变换得到时域上的 变形矩阵D;
[0064] (5)令￥1,」=〇[^,_(1 = 0，1，...，」=0，1，...，《，获得待求计算域内的弹性变 形Vi, j 0
[0065]综上，如图1所示，本发明基于FFT的深沟球轴承和角接触球轴承弹性变形快速计 算方法的具体计算过程为：
[0066] (1)给定/修正轴承载荷分布范围参数e ;
[0067] (2)通过线性插值得到载荷分布积分Jr(e);
[0068] (3)由式Sr = Pd/[2(l-2e)]计算得到内、外圈在径向载荷方向上的相对位移Sr;
[0069] (4)由式Frc = ZKn( Sr-Pd/2)1.5 Jr (e)计算得到计算外载荷Frc;
[0070] (5)判断计算外载荷Fr。是否收敛，若I (Frc-Fr)/Fr| >|，则按照 Frre修正e并重复步骤（1 )_步骤（4);若I (Fn-Fr )/Fr I彡|，则依据式
计算各滚球的载荷分布仏，执行下一步；
[0071] (6)依据各滚球的载荷分布，运用Hertz接触理论，计算接触椭圆长、短半轴长a、b;
[0072] (7)计算各滚球在内、外圈接触区上的接触应力分布p，将接触应力分布函数p(x, y)离散为(M+l) X (N+1)的矩阵p，然后补零使其成为2MX 2N的应力矩阵p' ；
[0073] (8)以坐标原点为中心，在4倍的计算区域内求得影响系数矩阵C，维数为2MX2N;
[0074] (9)对矩阵p'和C分别进行FFT变换，得到其FFT变换矩阵，和幻
[0075] (10)F和6的对应项相乘得到频域内的弹性变形矩阵i，对iH乍IFFT变换得到时域 上的弹性变形矩阵D;
[0076] (11)令Vi,j = DM+i,N+j(i=0,1，. . .，M, j = 0,1，. . .，N)，获得待求计算域内的弹性变 形V^，分别在各滚球与内、外圈接触区上进行步骤(6)_步骤(11)，得到各滚球与内、外圈间 的弹性变形分布。
[0077]最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通 过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在 形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。
【主权项】
1. 一种基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算方法，其特征在于:包 括以下步骤： Sl:给定/修正轴承载荷分布范围参数ε ; S2:查询载荷分布积分表并通过线性插值得到载荷分布积分Jr (ε); S3:由式Sr = Pd/ [ 2 (1 -2 ε)]计算得到内、外圈在径向载荷方向上的相对位移Sr，其中Pd表 示径向游隙； S4:由式Fr。= ZKn(Sr-Pd/2)15Jr(ε)计算得到计算外载荷F r。，其中Z表示滚动体数目，Kn 表示滚动体与内、外圈之间总的负荷-变形常数； S5:判断计算外载荷Frc与给定的径向载荷Fr之间的关系是否满足收敛条件，若I (Frc-Fr)/Fr| >ξ，则按照式8_=￡。1(1^。/^)_{!修正6并重复步骤31-步骤34，其中，|一般取 O · OO 1，β在O · OO 1与O · O 1之间取任意值；若I (Frc-Fr ) / Fr I彡ξ，则依据式 私C0Sft)r计胃求白勺^荷分#仏，=1，…⑷，队表示云加本白勺方{立 角，i = l，···，Ζ，执行下一步； S6:依据各滚球的载荷分布，运用Hertz接触理论，计算接触椭圆长、短半轴长a、b; S7 :计算各滚球在内、外圈接触区上的接触应力分布p，将接触应力分布函数p(X，y)离 散为(M+l) X (N+1)的矩阵p，然后补零使其成为2MX 2N的应力矩阵p' ； S8:以坐标原点为中心，在4倍的计算区域内求得影响系数矩阵C，维数为2MX2N; S9:对矩阵p '和C分别进行FFT变换，得到其FFT变换矩阵少和6 ; SlO: 和6的对应项相乘得到频域内的弹性变形矩阵对A作IFFT变换得到时域上的 弹性变形矩阵D; S11:令Vi,j = DM+i,N+j(i = 0，l，. . .，M, j = 0,l, . . . ,N)，获得待求计算域内的弹性变形 V1,J，分别在各滚球与内、外圈接触区上进行步骤S6-步骤S11，得到各滚球与内、外圈间的弹 性变形分布。2. 根据权利要求1所述的基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算方 法，其特征在于:在步骤S8中，影响系数矩阵C的求解采用Green函数法。
【文档编号】G06F17/14GK105893327SQ201610203061
【公开日】2016年8月24日
【申请日】2016年3月31日
【发明人】孟凡明, 刘越, 程志涛, 浦超
【申请人】重庆大学, 中国北方车辆研究所