CN105893327A - 基于fft的深沟球轴承和角接触球轴承弹性变形快速计算方法 - Google Patents

Abstract

本发明涉及一种基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算方法，属于滚动轴承优化设计技术领域。该方法基于深沟球轴承和角接触球轴承的点接触应力分布和影响系数矩阵，将离散循环卷积定理和FFT技术应用于轴承弹性变形计算中，计算各滚动体与内圈或外圈之间的弹性变形量及各滚动体与内外圈总的弹性变形量。本发明提供的方法能够有效解决传统方法计算效率不高的问题，因为使用传统方法计算深沟球轴承和角接触球轴承弹性变形时需要进行复杂的四重循环计算，严重耗费计算机内存且计算量很大，而本方法只要正确运用离散循环卷积定理和FFT技术，便能准确、快速地获得计算区域内的轴承弹性变形。

Description

基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算方法

技术领域

本发明属于滚动轴承优化设计技术领域，涉及一种基于FFT(Fast Fourier Transform：快速傅里叶变换)的深沟球轴承和角接触球轴承弹性变形快速计算方法。

背景技术

作为重要的动力传递部件，深沟球轴承和角接触球轴承在工程机械、轨道交通、航空航天及新能源等领域有着广泛的应用。深沟球轴承和角接触球轴承各滚动体与内圈或外圈之间弹性变形的快速求解有助于轴承动态特性(刚度、阻尼)分析和油膜厚度计算。在经典的轴承分析中普遍采用刚性套圈假设，即假定轴承套圈是不发生弹性变形的刚体，弹性变形只发生在滚动体上，滚动体与套圈的接触只会产生局部的接触变形而不会改变套圈的整体形状和尺寸。这种处理方法虽然为套圈的位移和轴承弹性变形分析带来了很大的方便，但是它无法准确反映出轴承内外套圈的变形量，从而无法准确预测轴承的油膜刚度、接触刚度和阻尼，也无法用于分析考虑轴承表面加工精度时的轴承弹性变形。此外，使用传统的Boussinesq方法计算深沟球轴承和角接触球轴承弹性变形时需要进行复杂的四重循环运算，严重耗费计算机内存且计算量很大，尤其是当网格密度较大时，计算量将呈指数级增长。

所以，传统的滚动轴承弹性变形计算方法更多的是具有统计与经验意义，其精度上不够准确、效率上不够高效。

发明内容

有鉴于此，本发明的目的在于提供一种基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算方法，该方法能够解决传统方法中存在的精度差、效率低的问题。

为达到上述目的，本发明提供如下技术方案：

一种基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算方法，其特征在于：包括以下步骤：

S1：给定/修正轴承载荷分布范围参数ε；

S2：查询载荷分布积分表并通过线性插值得到载荷分布积分J_r(ε)；

S3：由式δ_r＝P_d/[2(1-2ε)]计算得到内、外圈在径向载荷方向上的相对位移δ_r，其中P_d表示径向游隙；

S4：由式F_rc＝ZK_n(δ_r-P_d/2)^1.5J_r(ε)计算得到计算外载荷F_rc，其中Z表示滚动体数目，K_n表示滚动体与内、外圈之间总的负荷-变形常数；

S5：判断计算外载荷F_rc与给定的径向载荷F_r之间的关系是否满足收敛条件，若|(F_rc-F_r)/F_r|＞ξ，则按照式ε_new＝ε_old(F_rc/F_r)^-β修正ε并重复步骤S1-步骤S4，其中，ξ一般取0.001，β在0.001与0.01之间取任意值；若|(F_rc-F_r)/F_r|≤ξ，则依据式计算各滚球的载荷分布其中，ψ_i表示各滚动体的方位角，i＝1,…,Z，执行下一步；

S6：依据各滚球的载荷分布，运用Hertz接触理论，计算接触椭圆长、短半轴长a、b；

S7：计算各滚球在内、外圈接触区上的接触应力分布p，将接触应力分布函数p(x,y)离散为(M+1)×(N+1)的矩阵p，然后补零使其成为2M×2N的应力矩阵p’；

S8：以坐标原点为中心，在4倍的计算区域内求得影响系数矩阵C，维数为2M×2N；

S9：对矩阵p’和C分别进行FFT变换，得到其FFT变换矩阵和

S10：和的对应项相乘得到频域内的弹性变形矩阵对作IFFT变换(Inverse FastFourier Transform：快速傅里叶逆变换)得到时域上的弹性变形矩阵D；

S11：令V_i,j＝D_M+i,N+j(i＝0,1,...,M,j＝0,1,...,N)，获得待求计算域内的弹性变形V_i,j，分别在各滚球与内、外圈接触区上进行步骤S6-步骤S11，得到各滚球与内、外圈间的弹性变形分布。

进一步，在步骤S8中，影响系数矩阵C的求解采用Green函数法

本发明的有益效果在于：本方法能用于同时计算滚动体的弹性变形和内、外圈的弹性变形，也能用于考察滚动体和内、外滚道加工精度对轴承弹性变形的影响，能够使计算结果更加符合实际情况，也使得轴承刚度、阻尼和油膜厚度的计算结果更加准确，避免了传统Boussinesq方法中四重循环的运算，尤其是网格密度较大时，能极大程度地提高计算效率，本发明提出的方法对轴承的优化设计和轴承系统的受力分析是有益的。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为本发明的迭代计算方法流程图；

图2为深沟球轴承和角接触球轴承的滚球与内圈接触椭圆示意图；

图3为影响系数矩阵求解说明示意图；

图4为影响系数矩阵的计算区域示意图。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

图1为本发明的迭代计算方法流程图，对于给定了径向游隙P_d和径向载荷F_r的深沟球轴承和角接触球轴承，能够用图1流程图右半部分所示的控制误差的试解法求解各滚球的载荷分布：

首先假定轴承载荷分布范围参数ε的值，查询载荷分布积分表并通过线性插值得到载荷分布积分J_r(ε)。然后由式δ_r＝P_d/[2(1-2ε)]计算得到内、外圈在径向载荷方向上的相对位移δ_r，由式F_rc＝ZK_n(δ_r-P_d/2)^1.5J_r(ε)计算得到计算外载荷F_rc，其中Z表示滚动体数目，K_n表示滚动体与内、外圈之间总的负荷-变形常数。最后，判断计算外载荷F_rc是否满足|(F_rc-F_r)/F_r|≤ξ，若不满足，则按照式ε_new＝ε_old(F_rc/F_r)^-β修正ε并重复上述过程，其中，ξ一般取0.001，β在0.001与0.01之间取任意值；若满足收敛条件，则依据式计算各滚球的载荷分布其中，ψ_i表示各滚动体的方位角(i＝1,…,Z)。

由于考虑的是中低速轴承，故没有考虑离心力和陀螺力矩的影响，各滚球与内、外圈的接触角相同，以上计算出来的滚球载荷既表示滚球与内圈间的接触载荷也表示滚球与外圈间的接触载荷。此载荷为接触点上的单值，后续接触应力分布p的计算需用到此值。

根据Hertz接触理论，接触应力在接触区内按椭球体规律分布。当接触椭圆的长轴方向与x轴相重合时，接触应力分布p为：

$p(x,y)=\frac{3Q_{{\mathrm{\ψ}}_i}}{2\mathrm{\π}ab}(1-\frac{x^2}{a^2}-\frac{y^2}{b^2})---\left(1\right)$

图2为深沟球轴承和角接触球轴承的滚球与内圈接触椭圆示意图，其中，如图2，a、b分别表示接触椭圆的长、短半轴长，采用下列公式能够确定：

$a=k_a{\[\frac{3Q_{{\mathrm{\ψ}}_i}}{2E^{\′}(A+B)}\]}^{1/3}---\left(2\right)$

$b=k_b{\[\frac{3Q_{{\mathrm{\ψ}}_i}}{2E^{\′}(A+B)}\]}^{1/3}---\left(3\right)$

其中，k_a、k_b为接触椭圆尺寸系数，E'表示两接触体的综合弹性模量，常数A、B与两接触体的几何形状有关，其值由下列两方程确定：

$A+B=\frac{1}{2}(\frac{1}{R_{1x}}+\frac{1}{R_{1y}}+\frac{1}{R_{2x}}+\frac{1}{R_{2y}})---\left(4\right)$

$B-A=\frac{1}{2}{\[{(\frac{1}{R_{1x}}-\frac{1}{R_{1y}})}^2+{(\frac{1}{R_{2x}}-\frac{1}{R_{2y}})}^2+2(\frac{1}{R_{1x}}-\frac{1}{R_{1y}})(\frac{1}{R_{2x}}-\frac{1}{R_{2y}})\]}^{1/2}---\left(5\right)$

如图2，以滚球和内滚道接触为例，上两方程中的变量R_1x、R_1y、R_2x、R_2y下标中的第1个数字1和2分别表示滚球和内滚道，下标中的第2个字母x和y分别表示轴向和周向，也即：R_1x表示滚球的轴向曲率半径，R_1y表示滚球的周向曲率半径，R_2x表示内滚道的轴向曲率半径，R_2y表示内滚道的周向曲率半径。

求出常数A、B之后，根据式cosθ＝(B-A)/(A+B)计算θ值，然后根据θ值查询“接触椭圆尺寸系数表”能得到接触椭圆尺寸系数k_a、k_b的值，对于表中未列的θ值，插值可得到。

至此，每个滚球与内、外圈间的接触应力分布p均能依据上述方法得到。

深沟球轴承和角接触球轴承的弹性变形求解问题可等价为求解半无限体表面在上述求出的接触应力分布p下的弹性变形求解问题。其弹性变形方程可表达为：

$V(x,y)=\frac{2}{{\mathrm{\πE}}^{\′}}\∫{\∫}_{\mathrm{\Ω}}\frac{p(\mathrm{\ξ},\mathrm{\η})}{\sqrt{{(x-\mathrm{\ξ})}^2-{(y-\mathrm{\η})}^2}}d\mathrm{\ξ}d\mathrm{\η}---\left(6\right)$

使用传统的Boussinesq方法求解上述弹性变形方程时，需要将式(6)离散化，这种计算存在两个困难：首先，式(6)右端是奇异积分，奇点为(x,y)＝(ξ,η)；其次，为求某一点的弹性变形，就需要对整个求解域计算积分(需两重循环)，而为了求解整个计算域的弹性变形，又需用两重循环，故每一次弹性变形计算都需要进行四重循环运算。这严重耗费计算机内存且计算量很大，尤其是当网格密度较大时，计算量将呈指数级增长。例如，当网格密度M×N为32×32时，循环需进行1048576次，当网格密度M×N仅增大1倍为64×64时，循环次数便迅速增加到16777216次。可见，使用一般的数值积分方法，计算工作量太大，而且误差的频繁累积导致难以达到要求的精度。

在本发明方法中，基于线弹性叠加原理，为方便采用数值方法求解，式(6)可离散为：

$V(x_i,y_j)=\frac{2}{{\mathrm{\πE}}^{\′}}\underset{k=0}{\overset{M-1}{\Σ}}\underset{l=0}{\overset{N-1}{\Σ}}C(x_i-x_k,y_j-y_l)p(x_k,y_l)---\left(7\right)$

其中，C(x_i-x_k,y_j-y_l)为影响系数，简记为定义为在节点(k,l)上作用的单位力在节点(i,j)上产生的弹性变形，如图3所示。从式(7)能发现弹性变形的计算分为两步：影响系数的计算和卷积和的计算：

第一步：弹性变形影响系数的计算

根据定义，点接触弹性变形影响系数可按下式计算：

$C_{k,l}^{i,j}=\∫{\∫}_{\mathrm{\Δ}\mathrm{\Ω}}g(x_i-\mathrm{\ξ},y_j-\mathrm{\η})w(\mathrm{\ξ},\mathrm{\η})d\mathrm{\ξ}d\mathrm{\η}---\left(8\right)$

其中，ΔΩ为(x_k,y_l)点处的单元网格，如图3所示阴影部分。w(ξ,η)为插值函数，用来近似单元网格上的接触应力分布，g(x,y)称为响应函数(亦称Green函数)，当采用等间距网格时，只和两点之间的距离有关且存在对称性，因此影响系数的计算可只针对坐标原点进行，可进一步简记为C_k,l，故式(8)可简化为：

$C_{k,l}={\∫}_{X_k-\mathrm{\Δ}x/2}^{X_k+\mathrm{\Δ}x/2}{\∫}_{Y_l-\mathrm{\Δ}y/2}^{Y_l+\mathrm{\Δ}y/2}g(\mathrm{\ξ},\mathrm{\η})w(\mathrm{\ξ},\mathrm{\η})d\mathrm{\ξ}d\mathrm{\η}---\left(9\right)$

这里X_k、Y_l为力作用点和待求变形点间距离的坐标分量，如图3所示。

对于式(9)，可用Green函数法求解：

插值函数w(ξ,η)在单元网格内取1，同时假设Green函数在单元格内不变化，取节点值h(x_i,y_j)作为其近似值，则影响系数可简写为：

C_i,j＝g_i,j·ΔxΔy (10)

Green函数在原点存在一个奇异点，可通过以原点为中心的单元网格上的计算均值来代替原点处的离散值，即当|x_i|≤2Δx且|y_j|≤2Δy时，以式(11)计算值来代替点(x_i,yj)处的离散值：

$g(x_i,y_j)=\frac{1}{\mathrm{\Δ}x\mathrm{\Δ}y}{\∫}_{x_i-\mathrm{\Δ}x}^{x_i+\mathrm{\Δ}x}{\∫}_{y_j-\mathrm{\Δ}y}^{y_j+\mathrm{\Δ}y}g(x,y)dxdy---\left(11\right)$

当给定目标区域(计算区域)后，影响系数应在4倍的目标区域内计算，即在如图4所示的扩展区域内计算。当采用等间距网格时，只需计算出扩展区域内各节点相对于坐标原点的影响系数。根据对称性，只需计算一个象限内的影响系数便能确定其他象限内的影响系数。

第二步：弹性变形卷积和的计算

直接计算卷积和存在大量的重复运算，而运用离散循环卷积定理和FFT是一种高效的较为理想的计算方法。在计算出各滚球在内、外圈接触区上的接触应力分布p和影响系数矩阵C之后，基于离散循环卷积定理和FFT技术的深沟球轴承和角接触球轴承弹性变形计算如下：

(1)将接触应力分布函数p(x,y)离散为(M+1)×(N+1)的矩阵p，然后补零使其成为2M×2N的应力矩阵p’；

(2)以坐标原点为中心，在4倍的计算区域内求得影响系数矩阵C，维数为2M×2N；

(3)对矩阵p’和C分别进行FFT变换，得到其FFT变换矩阵和

(4)和的对应项相乘得到频域内的变形矩阵对作IFFT变换得到时域上的变形矩阵D；

(5)令V_i,j＝D_M+i,N+j(i＝0,1,...,M,j＝0,1,...,N)，获得待求计算域内的弹性变形V_i，j。

综上，如图1所示，本发明基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算方法的具体计算过程为：

(1)给定/修正轴承载荷分布范围参数ε；

(2)通过线性插值得到载荷分布积分J_r(ε)；

(3)由式δ_r＝P_d/[2(1-2ε)]计算得到内、外圈在径向载荷方向上的相对位移δ_r；

(4)由式F_rc＝ZK_n(δ_r-P_d/2)^1.5J_r(ε)计算得到计算外载荷F_rc；

(5)判断计算外载荷F_rc是否收敛，若|(F_rc-F_r)/F_r|＞ξ，则按照式ε_new＝ε_old(F_rc/F_r)^-β修正ε并重复步骤(1)-步骤(4)；若|(F_rc-F_r)/F_r|≤ξ，则依据式计算各滚球的载荷分布执行下一步；

(6)依据各滚球的载荷分布，运用Hertz接触理论，计算接触椭圆长、短半轴长a、b；

(7)计算各滚球在内、外圈接触区上的接触应力分布p，将接触应力分布函数p(x,y)离散为(M+1)×(N+1)的矩阵p，然后补零使其成为2M×2N的应力矩阵p’；

(8)以坐标原点为中心，在4倍的计算区域内求得影响系数矩阵C，维数为2M×2N；

(9)对矩阵p’和C分别进行FFT变换，得到其FFT变换矩阵和

(10)和的对应项相乘得到频域内的弹性变形矩阵对作IFFT变换得到时域上的弹性变形矩阵D；

(11)令V_i,j＝D_M+i,N+j(i＝0,1,...,M,j＝0,1,...,N)，获得待求计算域内的弹性变形V_i,j，分别在各滚球与内、外圈接触区上进行步骤(6)-步骤(11)，得到各滚球与内、外圈间的弹性变形分布。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。

Claims (2)

1.一种基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算方法，其特征在于：包括以下步骤：

S1：给定/修正轴承载荷分布范围参数ε；

S2：查询载荷分布积分表并通过线性插值得到载荷分布积分J_r(ε)；

S3：由式δ_r＝P_d/[2(1-2ε)]计算得到内、外圈在径向载荷方向上的相对位移δ_r，其中P_d表示径向游隙；

S4：由式F_rc＝ZK_n(δ_r-P_d/2)^1.5J_r(ε)计算得到计算外载荷F_rc，其中Z表示滚动体数目，K_n表示滚动体与内、外圈之间总的负荷-变形常数；

S5：判断计算外载荷F_rc与给定的径向载荷F_r之间的关系是否满足收敛条件，若|(F_rc-F_r)/F_r|＞ξ，则按照式ε_new＝ε_old(F_rc/F_r)^-β修正ε并重复步骤S1-步骤S4，其中，ξ一般取0.001，β在0.001与0.01之间取任意值；若|(F_rc-F_r)/F_r|≤ξ，则依据式计算各滚球的载荷分布其中，ψ_i表示各滚动体的方位角，i＝1,…,Z，执行下一步；

S6：依据各滚球的载荷分布，运用Hertz接触理论，计算接触椭圆长、短半轴长a、b；

S7：计算各滚球在内、外圈接触区上的接触应力分布p，将接触应力分布函数p(x,y)离散为(M+1)×(N+1)的矩阵p，然后补零使其成为2M×2N的应力矩阵p’；

S8：以坐标原点为中心，在4倍的计算区域内求得影响系数矩阵C，维数为2M×2N；

S9：对矩阵p’和C分别进行FFT变换，得到其FFT变换矩阵和

S10：和的对应项相乘得到频域内的弹性变形矩阵对作IFFT变换得到时域上的弹性变形矩阵D；

S11：令V_i,j＝D_M+i,N+j(i＝0,1,...,M,j＝0,1,...,N)，获得待求计算域内的弹性变形V_i,j，分别在各滚球与内、外圈接触区上进行步骤S6-步骤S11，得到各滚球与内、外圈间的弹性变形分布。

2.根据权利要求1所述的基于FFT的深沟球轴承和角接触球轴承弹性变形快速计算方法，其特征在于：在步骤S8中，影响系数矩阵C的求解采用Green函数法。