CN105138814A - 一种高速电主轴定位预紧下角接触球轴承极限预紧力分析方法 - Google Patents

Abstract

一种高速电主轴定位预紧下角接触球轴承极限预紧力分析方法，属于轴承预紧技术领域。本方法首先确定轴承的型号、工况条件以及材料最大许用应力等参数。然后基于Hertz接触理论，考虑了轴承的预紧力、离心力和陀螺力矩的影响，建立了定位预紧下角接触球轴承动力学模型，包括角接触球轴承接触区域模型、主轴静止状态下角接触球轴承预紧模型、主轴旋转状态下角接触球轴承定位预紧模型。而后，通过牛顿法求解出工况转速条件下，轴承滚动体与轴承内外圈最大接触应力与预紧力的关系。最终，分析出材料许用应力下轴承的预紧力，既极限预紧力。

Description

一种高速电主轴定位预紧下角接触球轴承极限预紧力分析方法

技术领域

本发明是一种高速电主轴定位预紧下角接触球轴承极限预紧力分析方法，属于轴承预紧技术领域。

背景技术

高速电主轴系统的轴承一般采用角接触球轴承，角接触球轴承的动态特性直接影响着高速电主轴系统的性能。为了增加轴承的刚度、提高旋转精度、降低振动噪声、延长轴承使用寿命，需对角接触球轴承进行预紧。定位预紧是角接触球轴承的主要预紧方式之一，是通过预先选定的内外圈隔套或垫圈使组配轴承内圈和外圈之间处于某一固定位置，从而使轴承获得合适的预紧，这种预紧方式的特点是轴系刚度较强，结构简单。在实际生产过程中，为了提高高速电主轴系统的刚度，工人一般会给轴承较大的预紧力。但是，如果预紧力过大，使得轴承滚动体与轴承内外圈之间的接触应力超过材料的许用应力，则会降低轴承的寿命，影响高速电主轴系统的性能。因此，发明一种高速电主轴定位预紧下角接触球轴承极限预紧力分析方法有着重要的意义。

目前确定极限定位预紧力的方法主要有两种，一是通过工人师傅的经验确定极限预紧力，二是对轴承进行理论分析，得到极限预紧力。显然，通过工人的经验确定预紧力的方法对工人有较高的要求，而且不能够准确地得到极限预紧力。相比之下，对轴承进行理论分析的方法可以较为科学地得到极限预紧力。国内外许多专家和学者一直在轴承预紧技术领域进行不懈地探索与研究，开展了多方面的工作。例如基于Hertz接触理论建立了5自由度分析模型，利用力平衡方程推导出轴承的刚度矩阵，用数值方法得到可用预紧力的范围；基于Algor仿真分析软件，建立轴承模型，得到不同预紧力下轴承的动力学参数，确定极限预紧力；基于Ansys仿真分析软件，分析不同工况条件下轴承的动力学特性，得到极限预紧力。

上述研究提出了很多预紧力分析与确定极限预紧力的方法，基本分析方法是：确定轴承参数、建立轴承模型、确定极限预紧力。但上述研究忽略了离心力和陀螺力矩的影响，建立的模型不够准确，致使得到的极限预紧力与实际情况存在误差。因此，本分析方法基于Hertz接触理论，考虑了轴承的预紧力、离心力和陀螺力矩的影响，建立了一种新的角接触球轴承动力学模型，通过分析不同工况参数下预紧力与最大接触应力的关系，得到轴承的极限预紧力。

本发明是一种高速电主轴定位预紧下角接触球轴承极限预紧力分析方法。

发明内容

本发明的目的是建立一种复杂工况条件下，高速电主轴定位预紧下角接触球轴承极限预紧力分析方法，以便更精确地得到定位预紧下角接触球轴承的极限预紧力，为电主轴系统设计与安装调试提供理论依据。现有的研究方法主要基于Hertz接触理论建立轴承模型，或通过仿真软件进行分析，最终确定极限预紧力。但现有方法在建立的模型时忽略了离心力和陀螺力矩的影响，影响了分析结果的准确性。本分析方法首先确定轴承的型号、工况条件以及材料最大许用应力等参数。然后基于Hertz接触理论，考虑了轴承的预紧力、离心力和陀螺力矩的影响，建立了定位预紧下角接触球轴承动力学模型，包括：角接触球轴承接触区域模型、主轴静止状态下角接触球轴承预紧模型、主轴旋转状态下角接触球轴承定位预紧模型。而后，通过牛顿法求解出工况转速条件下，轴承滚动体与轴承内外圈最大接触应力与预紧力的关系。最终，分析出材料许用应力下轴承的预紧力，既极限预紧力。

如图1所示，本发明提供的一种高速电主轴定位预紧下角接触球轴承极限预紧力分析方法包括以下步骤。

S1.确定轴承的型号、工况条件以及材料最大许用应力等参数。

确定轴承的型号，得到轴承内径d、外径D、接触角α、滚动体直径D_w、滚动体数目Z，以及内外圈沟道曲率半径系数f_i和f_e等参数。

确定轴承工况条件，得到轴承工况转速n等参数。

确定轴承的材料，得到材料最大许用应力[σ]等参数。

这些数据为后续轴承建模以及确定极限预紧力提供数据。

S2.建立定位预紧下角接触球轴承动力学模型。

S2.1建立角接触球轴承接触区域模型。

法向集中力P作用在弹性半空间的(x′,y′)点，而在另一点(x,y)产生的法向位移ω(x,y)由弹性理论Boussinesq解给出：

$\mathrm{\ω}(x,y)=\frac{1-v^2}{\mathrm{\πE}}\frac{P}{\sqrt{{(x-x^{\′})}^2+{(y-y^{\′})}^2}}---\left(1\right)$

式中，E,v分别是弹性体的弹性模量和泊松比。

当P是半空间表面局部区域S′_c上的分布压力P(x′,y′)时，(x,y)点的法向位移ω(x,y)可以表示为：

$\mathrm{\ω}(x,y)=\frac{1-v^2}{\mathrm{\πE}}{\∫\∫}_{S_c^{\′}}\frac{P}{\sqrt{{(x-x^{\′})}^2+{(y-y^{\′})}^2}}d{x}^{\′}d{y}^{\′}---\left(2\right)$

压力P按半椭球函数分布，表示为：

$P(x^{\′},y^{\′})=P_0\sqrt{1-{(x^{\′}/a)}^2-{(y^{\′}/b)}^2}---\left(3\right)$

式中，a，b分别为椭圆区域S′_c上的半长轴和半短轴，P₀为椭圆中心处的最大压应力。由式(1)(2)(3)得到分布压力P(x′,y′)在(x,y)点产生的位移为：

$\mathrm{\ω}(x,y)=\frac{1-v^2}{\mathrm{\πE}}{\∫\∫}_{S_c^{\′}}\frac{1-{\mathrm{\υ}}^2}{\mathrm{\πE}}\frac{P}{\sqrt{{(x-x^{\′})}^2+{(y-y^{\′})}^2}}d{x}^{\′}d{y}^{\′}---\left(4\right)$

通过坐标变化得到：

$\mathrm{\ω}=\frac{(1-v^2)b{P}_0}{2a^{2}E}{\∫}_0^{\mathrm{\π}}\frac{a^2-x^{2}s{\mathrm{in}}^2\mathrm{\θ}-2\mathrm{xy}\sin \mathrm{\θ}\cos \mathrm{\θ}-{(a/b)}^2{y}^2{\cos }^2\mathrm{\θ}}{\sqrt{{\cos }^2\mathrm{\θ}+{(b/a)}^2{\sin }^2\mathrm{\θ}}}\mathrm{d\θ}---\left(5\right)$

式中，P₀为接触椭圆区域中心处的最大压应力。a,b分别为接触椭圆区域S′_c的半长轴和半短轴。

令：

e²＝1-(b/1)²,b<a(6)

${\&Integral;}_0^{\mathrm{\&pi;}}\frac{\mathrm{d\&theta;}}{\sqrt{{\cos }^2\mathrm{\&theta;}+{(b/a)}^2{\sin }^2\mathrm{\&theta;}}}=2K\left(e\right)---\left(7\right)$

${\&Integral;}_0^{\mathrm{\&pi;}}\frac{{\sin }^2\mathrm{\&theta;d\&theta;}}{\sqrt{{\cos }^2\mathrm{\&theta;}+{(b-a)}^2{\sin }^2\mathrm{\&theta;}}}=\frac{2}{e^2}(K\left(e\right)-E\left(e\right))---\left(8\right)$

${\&Integral;}_0^{\mathrm{\&pi;}}\frac{\sin \mathrm{\&theta;}\cos \mathrm{\&theta;d\&theta;}}{\sqrt{{\cos }^2\mathrm{\&theta;}+{(b/a)}^2{\sin }^2\mathrm{\&theta;}}}=0---\left(9\right)$

${\&Integral;}_0^{\mathrm{\&pi;}}\frac{{\cos }^2\mathrm{\&theta;d\&theta;}}{\sqrt{{\cos }^2\mathrm{\&theta;}+{(b/a)}^2{\sin }^2\mathrm{\&theta;}}}=\frac{2}{e^2}(E\left(e\right)-(1-e^2)K\left(e\right))---\left(10\right)$

式中，K(e)和E(e)分别为第一类和第二类完全椭圆积分函数，e为椭圆参数。

根据式(6)～(10)，式(5)表达为：

$\mathrm{\&omega;}=\frac{1-v^2}{E}(L-M{x}^2-N{y}^2)--\left(11\right)$

式中：

L＝P₀bK(e)(12)

$M=\frac{b{P}_0}{a^2{e}^2}(K\left(e\right)-E\left(e\right))---\left(13\right)$

$N=\frac{b{P}_0}{a^2{e}^2}(\frac{a^2}{b^2}E\left(e\right)-K\left(e\right))---\left(14\right)$

通过接触应力公式，得到外部载荷Q：

${\&Integral;\&Integral;}_{S_c^{\&prime;}}{P}_0\sqrt{1-{(x/a)}^2-{(x/b)}^2}\mathrm{dxdy}=2/3\mathrm{\&pi;ab}{P}_0=Q---\left(15\right)$

设两个弹性体V₁和V₂在未施加载荷前仅在O点出相互接触。V₁和V₂在O点处的主曲率半径分别为R₁₁、R₁₂和R₂₁、R₂₂，如图2所示。对曲率1/R_ij的正负做出如下规定：凸出的表面曲率为正，凹进的表面曲率为负。

在载荷作用下，两物体的接触区域很小，因此在O点附近用二次函数来近似描述物体的表面方程，如图3所示。设DD′是两表面之间与公切面垂直的线段，D、D′的坐标分别是(x₁,y₁)和(x₂y₂)，将DD′的距离用z表示，则这两点距离为：

$z=z_1+z_2=\frac{1}{2}(\frac{x_1^2}{R_{11}}+\frac{y_1^2}{R_{12}}+\frac{x_2^2}{R_{21}}+\frac{y_2^2}{R_{22}})---\left(16\right)$

式中，z₁、z₂分别是两个未变形物体表面对应点到到初始接触点之间的垂直距离。

因为角接触的滚动体为钢球，所以滚动体和它的接触体主曲率重合。利用坐标变化，式(16)可以变换成标准椭圆方程，因此：

$z=\frac{1}{2}(\frac{1}{R_{11}}+\frac{1}{R_{21}})x^2+\frac{1}{2}(\frac{1}{R_{12}}+\frac{1}{R_{22}})y^2={\mathrm{Ax}}^2+{\mathrm{By}}^2---\left(17\right)$

式中：

$B-A=\frac{1}{2}{\{{[(\frac{1}{R_{12}}-\frac{1}{R_{21}})+(\frac{1}{R_{22}}-\frac{1}{R_{11}})]}^2-4(\frac{1}{R_{12}}-\frac{1}{R_{11}})(\frac{1}{R_{22}}-\frac{1}{R_{21}}){\sin }^2\mathrm{\&alpha;}\}}^{1/2}---\left(18\right)$

$B+A=\frac{1}{2}(\frac{1}{R_{11}}+\frac{1}{R_{12}}+\frac{1}{R_{21}}+\frac{1}{R_{22}})---\left(19\right)$

如图4所示，接触体表面的位移δ由两部分组成，分别是接触体的刚性位移δ₁与δ₂。ω₁(x,y)、ω₂(x,y)、z₁(x,y)、z₂(x,y)为原点以外的点相对于接触平面的位移。在接触区内，满足变形协调条件：

δ＝δ₁+δ₂＝ω₁+ω₂+z₁+z₂(20)

由式(18)～(20)得椭圆率参数e的方程：

$\frac{B-A}{B+A}=F\left(\mathrm{\&rho;}\right)=\frac{2(1-e^2)}{e^2}\frac{E\left(e\right)-K\left(e\right)}{E\left(e\right)}+1---\left(21\right)$

式中：

$a={\left(\frac{2E\left(e\right)}{\mathrm{\&pi;}(1-e^2)}\right)}^{1/3}{\left(\frac{3Q}{2\mathrm{\&Sigma;\&rho;}{E}^{\&prime;}}\right)}^{1/3}---\left(22\right)$

$b={\left[\frac{2\sqrt{1-e^2}E\left(e\right)}{\mathrm{\&pi;}}\right]}^{1/3}{\left(\frac{3Q}{2\mathrm{\&Sigma;\&rho;}{E}^{\&prime;}}\right)}^{1/3}---\left(23\right)$

$\mathrm{\&delta;}=\frac{2K\left(e\right)}{\mathrm{\&pi;}}{\left(\frac{(1-e^2)\mathrm{\&pi;}}{2E\left(e\right)}\right)}^{1/3}{\left(\frac{3Q}{2\mathrm{\&Sigma;\&rho;}{E}^{\&prime;}}\right)}^{2/3}\frac{\mathrm{\&Sigma;\&rho;}}{2}---\left(24\right)$

$\mathrm{\&Sigma;\&rho;}=\frac{1}{R_{11}}+\frac{1}{R_{12}}+\frac{1}{R_{21}}+\frac{1}{R_{22}}---\left(25\right)$

$\frac{1}{E^{\&prime;}}=\frac{1-{\mathrm{\&upsi;}}_1^2}{E_1}+\frac{1-{\mathrm{\&upsi;}}_2^2}{E_2}---\left(26\right)$

式中，δ为接触体表面的位移；Σρ为曲率和函数；v₁、E₁和v₂、E₂分别是物体V₁和V₂的泊松比和弹性模量。

S2.2建立主轴静止状态下角接触球轴承预紧模型。

主轴静止状态下，角接触球轴承受到轴向预紧力F_a，每一个滚动体将承受相同的载荷并产生相同的变形。设O_i，O_e为初始状态下内外圈沟道曲率中心，α和α′为初始接触角和实际接触角，Q_n为法向接触力。若轴承外圈固定，内圈受到轴向预紧力，内圈将产生轴向位移，内圈沟道曲率中心从O_i移动到O′_i，如图5所示。

由几何关系可以得到O_iO_e和O′_iO_e的线段长度：

O_iO_e＝(f_i+f_e-1)D_w(27)

$O_i^{\&prime;}{O}_e=\frac{d}{{\cos \mathrm{\&alpha;}}^{\&prime;}}=O_i{O}_e\frac{\cos \mathrm{\&alpha;}}{{\cos \mathrm{\&alpha;}}^{\&prime;}}---\left(28\right)$

式中，f_i与f_e分别为轴承内外圈沟道曲率系数。

线段O′_iO_e与O_iO_e的距离之差是滚动体与内外滚道的法向接触变形δ_n。

δ_n＝O′_iO_e-O_iO_e(29)

法向接触载荷为：

Q_n＝K_nδ_n(30)

式中，Z为钢球数目,K_n为钢度系数。

$K_n=2.1343*{10}^5({\mathrm{\&delta;}}_i^{*}{\mathrm{\&Sigma;\&rho;}}_i^{\frac{1}{3}}+{\mathrm{\&delta;}}_o^{*}{\mathrm{\&Sigma;\&rho;}}_o^{\frac{1}{3}})---\left(31\right)$

轴承的力平衡方程：

ZQ_nsinα′＝F_a(32)

将式(27)～(31)带入式(32)：

$\frac{F_a}{{\mathrm{ZK}}_n{O}_i{O}_e}={\sin \mathrm{\&alpha;}}^{\&prime;}{(\frac{\cos \mathrm{\&alpha;}}{{\cos \mathrm{\&alpha;}}^{\&prime;}}-1)}^{1.5}---\left(33\right)$

由式(33)可以求得实际接触角α′。

钢球的曲率为：

$\frac{1}{R_{11}}=\frac{1}{R_{12}}=\frac{2}{D_w}---\left(34\right)$

对于内滚道，沟道曲率为：

$\frac{1}{R_{i1}}=-\frac{1}{f_i{D}_w}---\left(35\right)$

$\frac{1}{R_{i2}}=\frac{2\cos \mathrm{\&alpha;}}{d_m-D_w\cos \mathrm{\&alpha;}}---\left(36\right)$

对于外滚道，沟道曲率为：

$\frac{1}{R_{o1}}=-\frac{1}{f_e{D}_w}---\left(37\right)$

$\frac{1}{R_{o2}}=\frac{2\cos \mathrm{\&alpha;}}{d_m+D_w\cos \mathrm{\&alpha;}}---\left(38\right)$

将式(34)～(38)带入式(21)可以建立方程f₁,f₂：

$f_1=\frac{2(1-e_i^2)}{e_i^2}\frac{E\left(e_i\right)-K\left(e_i\right)}{E\left(e_i\right)}+1-\frac{B_i-A_i}{B_i+A_i}=0---\left(39\right)$

$f_2=\frac{2(1-e_o^2)}{e_o^2}\frac{E\left(e_o\right)-K\left(e_o\right)}{E\left(e_o\right)}+1-\frac{B_o-A_o}{B_o+A_o}=0---\left(40\right)$

S2.3建立主轴旋转状态下角接触球轴承定位预紧模型。

定位预紧的轴承在使用过程中，其内外圈相对位置是不会改变的。主轴高速旋转时，轴承受到的离心力F_ck和陀螺力矩M_gk的影响，滚动体会产生位移，如图6所示。图7中，O′为滚动体初始形心位置，O″为主轴旋转时滚动体实际形心位置，D为内圈沟道曲率中心，B为外圈沟道曲率中心。

主轴高速旋转时，轴承受到的离心力F_ck和陀螺力矩M_gk。

$F_{\mathrm{ck}}=\frac{1}{2}m{D}_w{\mathrm{\&Omega;}}^2{\left(\frac{{\mathrm{\&Omega;}}_E}{\mathrm{\&Omega;}}\right)}^2---\left(41\right)$

$M_{\mathrm{gk}}=J_b{\mathrm{\&Omega;}}^2{\left(\frac{{\mathrm{\&Omega;}}_E}{\mathrm{\&Omega;}}\right)}_k\sin {\mathrm{\&alpha;}}_k---\left(42\right)$

式中，Ω为主轴旋转角速度；Ω_E为滚动体绕着主轴公转角速度；Ω_B为滚动体自转角速度；J_b为滚动体转动惯量。

将离心力离心力F_ck和陀螺力矩M_gk按照轴向和径向进行分解得到平衡方程f₃,f₄：

$f_3=F_{\mathrm{ck}}+\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}+Q_{\mathrm{ik}}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}=0---\left(43\right)$

$f_4=Q_{\mathrm{ik}}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}+\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}=0---\left(44\right)$

由图中几何关系得到方程f₅,f₆：

f₅＝BO″*cosθ_ok+DO″*cosθ_ik-BD*cosα′＝0(45)

f₆＝BO″*sinθ_ok+DO″*sinθ_ik-BD*sinα′＝0(46)

式中，θ_ik与θ_ok分别为内外圈接触角，Q_ik与Q_ok分别为内外圈法向接触力。

DO″＝(f_i-0.5)D_w+δ_ik(47)

BO″＝(f_e-0.5)D_w+δ_ok(48)

S3.通过牛顿法求解出工况条件下最大接触应力与预紧力的关系。

将轴承内外圈椭圆率e_i，e_o，内外圈动态接触角θ_ik，θ_ok，内外圈法向接触力Q_ik，Q_ok这6个参数设为未知量，联立方程f₁到f₆得到非线性方程组。

$\left\{\begin{array}{c}{f}_1=\frac{2(1-e_i^2)}{e_i^2}\frac{E\left(e_i\right)-K\left(e_i\right)}{E\left(e_i\right)}+1-\frac{B_i-A_i}{B_i+A_i}\\ f_2=\frac{2(1-e_o^2)}{e_o^2}\frac{E\left(e_o\right)-K\left(e_o\right)}{E\left(e_o\right)}+1-\frac{B_o-A_o}{B_o+A_o}\\ f_3=F_{\mathrm{ck}}+\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}+Q_{\mathrm{ik}}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}\\ f_4=Q_{\mathrm{ik}}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}+\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}\\ f_5={\mathrm{BO}}^{\&prime;\&prime;}*\cos {\mathrm{\&theta;}}_{\mathrm{ik}}+{\mathrm{DO}}^{\&prime;\&prime;}*\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-\mathrm{BD}*\cos {\mathrm{\&alpha;}}^{\&prime;}\\ f_6={\mathrm{BO}}^{\&prime;\&prime;}*\sin {\mathrm{\&theta;}}_{\mathrm{ok}}+{\mathrm{DO}}^{\&prime;\&prime;}*\sin {\mathrm{\&theta;}}_{\mathrm{ik}}-\mathrm{BD}*\sin {\mathrm{\&alpha;}}^{\&prime;}\end{array}\right.---\left(49\right)$

应用MATLAB数值分析软件，通过牛顿法进行数值迭代，求解该方程组。

方程组的Jacobi矩阵：

$J\left(x\right)=\left[\begin{array}{cccccc}\frac{{\&PartialD;f}_1}{\&PartialD;e_i}& \frac{\&PartialD;f_1}{\&PartialD;e_o}& \frac{\&PartialD;f_1}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{\&PartialD;f_1}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{\&PartialD;f_1}{\&PartialD;Q_{\mathrm{ik}}}& \frac{\&PartialD;f_1}{\&PartialD;Q_{\mathrm{ok}}}\\ \frac{\&PartialD;f_2}{\&PartialD;e_i}& \frac{\&PartialD;f_2}{\&PartialD;e_o}& \frac{\&PartialD;f_2}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{\&PartialD;f_2}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{\&PartialD;f_2}{\&PartialD;Q_{\mathrm{ik}}}& \frac{\&PartialD;f_2}{\&PartialD;Q_{\mathrm{ok}}}\\ \frac{\&PartialD;f_3}{\&PartialD;e_i}& \frac{\&PartialD;f_3}{\&PartialD;e_o}& \frac{\&PartialD;f_3}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{\&PartialD;f_3}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{\&PartialD;f_3}{\&PartialD;Q_{\mathrm{ik}}}& \frac{\&PartialD;f_3}{\&PartialD;Q_{\mathrm{ok}}}\\ \frac{\&PartialD;f_4}{\&PartialD;e_i}& \frac{\&PartialD;f_4}{\&PartialD;e_o}& \frac{\&PartialD;f_4}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{\&PartialD;f_4}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{\&PartialD;f_4}{\&PartialD;Q_{\mathrm{ik}}}& \frac{\&PartialD;f_4}{\&PartialD;Q_{\mathrm{ok}}}\\ \frac{\&PartialD;f_5}{\&PartialD;e_i}& \frac{\&PartialD;f_5}{\&PartialD;e_o}& \frac{\&PartialD;f_5}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{\&PartialD;f_5}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{\&PartialD;f_5}{\&PartialD;Q_{\mathrm{ik}}}& \frac{\&PartialD;f_5}{\&PartialD;Q_{\mathrm{ok}}}\\ \frac{\&PartialD;f_6}{\&PartialD;e_i}& \frac{\&PartialD;f_6}{\&PartialD;e_o}& \frac{\&PartialD;f_6}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{\&PartialD;f_6}{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{\&PartialD;f_6}{\&PartialD;Q_{\mathrm{ik}}}& \frac{\&PartialD;f_6}{\&PartialD;Q_{\mathrm{ok}}}\end{array}\right]---\left(50\right)$

迭代矩阵为：

J(x^(k))d^(k)＝-F(x^(k))(51)

x^(k+1)＝x^(k)+d^(k)(52)

终止条件：

|x^(k+1)-x^(k)|<ε(53)

取初始矩阵x⁽⁰⁾，通过式(50)(51)(52)计算得到x⁽¹⁾，再将x⁽¹⁾带入式(50)(51)(52)计算得到x⁽²⁾，循环此计算，当第k+1次计算结果x^(k+1)与第k次计算结果x^(k)的差值矩阵小于给定的误差ε时，计算停止，认为矩阵x^(k+1)为所求结果。

再根据式(15)中最大接触应力P₀与内外圈法向接触力Q的关系，可以得到工况转速n下，轴承内外圈与滚动体最大接触应力P₀与预紧力F_a的关系。

S4.根据材料最大许用应力确定极限预紧力。

根据工况转速n下，轴承内外圈与滚动体最大接触应力P₀与预紧力F_a的关系，得到轴承内外圈与滚动体最大接触应力P₀为材料最大许用应力[σ]时对应的预紧力F_a，既极限预紧力F_amax。

与现有技术相比，本发明具有以下优点：

1.建立动力学模型时，本发明考虑了轴承的预紧力、离心力和陀螺力矩的影响，使所建模型更加接近实际情况，提高了本分析方法的准确性。

2.本发明适用于不同型号、不同工况条件下角接触球轴承定位预紧下极限预紧力的分析，其中建模方法与求解方法可重复性较强，进行其他型号角接触球轴承定位预紧下极限预紧力的分析时，只需重新输入参数，可以提高分析工作效率。

附图说明

图1是一种高速电主轴定位预紧下角接触球轴承极限预紧力分析方法流程图。

图2是点接触示意图。

图3是表面间距示意图。

图4是变形协调关系示意图。

图5是静止预紧下轴承位移关系图。

图6是定位预紧下轴承受力图。

图7是定位预紧下轴承位移关系图。

图8是定位预紧下轴承最大接触应力与预紧力关系图。

图9是辅助线做图法得到极限预紧力方法示意图。

具体实施方式

首先，确定轴承确定轴承的型号、工况条件以及材料最大许用应力等参数。

本文选用Si₃N₄陶瓷球轴承作为算例，Si₃N₄陶瓷球轴承的参数如下表所示，主轴工况转速n范围0r/min到12000r/min，Si₃N₄陶瓷材料最大许用应力[σ]为2000MPa。

其次，根据上述理论，建立定位预紧下角接触球轴承动力学模型，其中包括：角接触球轴承接触区域模型、主轴静止状态下角接触球轴承预紧模型、主轴旋转状态下角接触球轴承定位预紧模型。

得到椭圆率参数方程f₁,f₂：

$f_1=\frac{2(1-e_i^2)}{e_i^2}\frac{E\left(e_i\right)-K\left(e_i\right)}{E\left(e_i\right)}+1-\frac{B_i-A_i}{B_i+A_i}=0$

$f_2=\frac{2(1-e_o^2)}{e_o^2}\frac{E\left(e_o\right)-K\left(e_o\right)}{E\left(e_o\right)}+1-\frac{B_o-A_o}{B_o+A_o}=0$

得到力平衡方程f₃,f₄：

$f_3=F_{\mathrm{ck}}+\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}+Q_{\mathrm{ik}}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}=0$

$f_4=Q_{\mathrm{ik}}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}+\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}=0$

得到几何关系到方程f₅,f₆：

f₅＝BO″*cosθ_ok+DO″*cosθ_ik-BD*cosα′＝0

f₆＝BO″*sinθ_ok+DO″*sinθ_ik-BD*sinα′＝0

然后，通过牛顿法求解出工况条件下，最大接触应力与预紧力的关系

将轴承内外圈椭圆率e_i，e_o，内外圈动态接触角θ_ik，θ_ok，内外圈法向接触力Q_ik，Q_ok这6个参数设为未知量，联立方程f₁到f₆可以得到非线性方程组。

$\left\{\begin{array}{c}{f}_1=\frac{2(1-e_i^2)}{e_i^2}\frac{E\left(e_i\right)-K\left(e_i\right)}{E\left(e_i\right)}+1-\frac{B_i-A_i}{B_i+A_i}\\ f_2=\frac{2(1-e_o^2)}{e_o^2}\frac{E\left(e_o\right)-K\left(e_o\right)}{E\left(e_o\right)}+1-\frac{B_o-A_o}{B_o+A_o}\\ f_3=F_{\mathrm{ck}}+\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}+Q_{\mathrm{ik}}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}\\ f_4=Q_{\mathrm{ik}}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}+\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}\\ f_5={\mathrm{BO}}^{\&prime;\&prime;}*\cos {\mathrm{\&theta;}}_{\mathrm{ik}}+{\mathrm{DO}}^{\&prime;\&prime;}*\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-\mathrm{BD}*\cos {\mathrm{\&alpha;}}^{\&prime;}\\ f_6={\mathrm{BO}}^{\&prime;\&prime;}*\sin {\mathrm{\&theta;}}_{\mathrm{ok}}+{\mathrm{DO}}^{\&prime;\&prime;}*\sin {\mathrm{\&theta;}}_{\mathrm{ik}}-\mathrm{BD}*\sin {\mathrm{\&alpha;}}^{\&prime;}\end{array}\right.$

应用MATLAB数值分析软件，通过牛顿法进行数值迭代，求解该方程组。

再根据式(6)中最大接触应力P₀与内外圈法向接触力Q的关系，得到工况转速n下，轴承内外圈与滚动体最大接触应力P₀与预紧力F_a的关系，如图8所示。

最后，根据材料最大许用应力确定极限预紧力。

Si3N4陶瓷材料最大许用应力[σ]为2000MPa，可以在图8中做出最大接触应力P₀在2000MPa时的辅助线，主轴工况转速n范围0r/min到12000r/min，从图中可以看出预紧力F_a相同的情况下，转速n越高，最大接触应力P₀越大，所以，主轴转速n在12000r/min，最大接触应力P₀在2000MPa时对应的预紧力F_a为极限预紧力F_amax。通过做辅助线的办法可以快速得到最大预紧力F_amax为1800N，如图9所示。

Claims (1)

1.一种高速电主轴定位预紧下角接触球轴承极限预紧力分析方法，其特征在于：该方法包括以下步骤；

S1.确定轴承的型号、工况条件以及材料最大许用应力参数；

确定轴承的型号，得到轴承内径d、外径D、接触角α、滚动体直径D_w、滚动体数目Z，以及内外圈沟道曲率半径系数f_i和f_e参数；

确定轴承工况条件，得到轴承工况转速n参数；

确定轴承的材料，得到材料最大许用应力[σ]参数；

这些数据为后续轴承建模以及确定极限预紧力提供数据；

S2.建立定位预紧下角接触球轴承动力学模型；

S2.1建立角接触球轴承接触区域模型；

法向集中力P作用在弹性半空间的(x′,y′)点，而在另一点(x,y)产生的法向位移ω(x,y)由弹性理论Boussinesq解给出：

$\mathrm{\&omega;}(x,y)=\frac{1-v^2}{\mathrm{\&pi;E}}\frac{P}{\sqrt{{(x-x^{\&prime;})}^2+{(y-y^{\&prime;})}^2}}---\left(1\right)$

式中，E,v分别是弹性体的弹性模量和泊松比；

当P是半空间表面局部区域S′_c上的分布压力P(x′,y′)时，(x,y)点的法向位移ω(x,y)表示为：

$\mathrm{\&omega;}(x,y)=\frac{1-v^2}{\mathrm{\&pi;E}}{\&Integral;\&Integral;}_{S_c^{\&prime;}}\frac{P}{\sqrt{{(x-x^{\&prime;})}^2+{(y-y^{\&prime;})}^2}}{\mathrm{dx}}^{\&prime;}{\mathrm{dy}}^{\&prime;}---\left(2\right)$

压力P按半椭球函数分布，表示为：

$P(x^{\&prime;},y^{\&prime;})=P_0\sqrt{1-{(x^{\&prime;}/a)}^2-{(y^{\&prime;}/b)}^2}---\left(3\right)$

式中，a，b分别为椭圆区域S′_c上的半长轴和半短轴，P₀为椭圆中心处的最大压应力；由式(1)(2)(3)得到分布压力P(x′,y′)在(x,y)点产生的位移为：

$\mathrm{\&omega;}(x,y)=\frac{1-v^2}{\mathrm{\&pi;E}}{\&Integral;\&Integral;}_{S_c^{\&prime;}}\frac{1-{\mathrm{\&upsi;}}^2}{\mathrm{\&pi;E}}\frac{P}{\sqrt{{(x-x^{\&prime;})}^2+{(y-y^{\&prime;})}^2}}{\mathrm{dx}}^{\&prime;}{\mathrm{dy}}^{\&prime;}---\left(4\right)$

通过坐标变化得到：

$\mathrm{\&omega;}=\frac{(1-v^2){\mathrm{bP}}_0}{{2a}^{2}E}{\&Integral;}_0^{\mathrm{\&pi;}}\frac{a^2-x^2{\sin }^2\mathrm{\&theta;}-2\mathrm{xy}\sin \mathrm{\&theta;}\cos \mathrm{\&theta;}-{(a/b)}^2{y}^2{\cos }^2\mathrm{\&theta;}}{\sqrt{{\cos }^2\mathrm{\&theta;}+{(b/a)}^2{\sin }^2\mathrm{\&theta;}}}\mathrm{d\&theta;}---\left(5\right)$

式中，P₀为接触椭圆区域中心处的最大压应力；a,b分别为接触椭圆区域S′_c的半长轴和半短轴；

令：

e²＝1-(b/1)²,b<a(6)

${\&Integral;}_0^{\mathrm{\&pi;}}\frac{\mathrm{d\&theta;}}{\sqrt{{\cos }^2\mathrm{\&theta;}+{(b/a)}^2{\sin }^2\mathrm{\&theta;}}}=2K\left(e\right)---\left(7\right)$

${\&Integral;}_0^{\mathrm{\&pi;}}\frac{{\sin }^2\mathrm{\&theta;d\&theta;}}{\sqrt{{\cos }^2\mathrm{\&theta;}+{(b-a)}^2{\sin }^2\mathrm{\&theta;}}}=\frac{2}{e^2}(K\left(e\right)-E\left(e\right))---\left(8\right)$

${\&Integral;}_0^{\mathrm{\&pi;}}\frac{\sin \mathrm{\&theta;}\cos \mathrm{\&theta;d\&theta;}}{\sqrt{{\cos }^2\mathrm{\&theta;}+{(b/a)}^2{\sin }^2\mathrm{\&theta;}}}=0---\left(9\right)$

${\&Integral;}_0^{\mathrm{\&pi;}}\frac{{\cos }^2\mathrm{\&theta;d\&theta;}}{\sqrt{{\cos }^2\mathrm{\&theta;}+{(b/a)}^2{\sin }^2\mathrm{\&theta;}}}=\frac{2}{e^2}(E\left(e\right)-(1-e^2)K\left(e\right))---\left(10\right)$

式中，K(e)和E(e)分别为第一类和第二类完全椭圆积分函数，e为椭圆参数；

根据式(6)～(10)，式(5)表达为：

$\mathrm{\&omega;}=\frac{1-v^2}{E}(L-{\mathrm{Mx}}^2-{\mathrm{Ny}}^2)---\left(11\right)$

式中：

L＝P₀bK(e)(12)

$M=\frac{{\mathrm{bP}}_0}{a^2{e}^2}(K\left(e\right)-E\left(e\right))---\left(13\right)$

$N=\frac{{\mathrm{bP}}_0}{a^2{e}^2}(\frac{a^2}{b^2}E\left(e\right)-K\left(e\right))---\left(14\right)$

通过接触应力公式，得到外部载荷Q：

${\&Integral;\&Integral;}_{S_c^{\&prime;}}{P}_0\sqrt{1-{(x/a)}^2-{(x/b)}^2}\mathrm{dxdy}=2/3\mathrm{\&pi;ab}{P}_0=Q---\left(15\right)$

设两个弹性体V₁和V₂在未施加载荷前仅在O点出相互接触；V₁和V₂在O点处的主曲率半径分别为R₁₁、R₁₂和R₂₁、R₂₂；对曲率1/R_ij的正负做出如下规定：凸出的表面曲率为正，凹进的表面曲率为负；

在载荷作用下，两物体的接触区域很小，因此在O点附近用二次函数来近似描述物体的表面方程；设DD′是两表面之间与公切面垂直的线段，D、D′的坐标分别是(x₁,y₁)和(x₂y₂)，将DD′的距离用z表示，则这两点距离为：

$z=z_1+z_2=\frac{1}{2}(\frac{x_1^2}{R_{11}}+\frac{y_1^2}{R_{12}}+\frac{x_2^2}{R_{21}}+\frac{y_2^2}{R_{22}})---\left(16\right)$

式中，z₁、z₂分别是两个未变形物体表面对应点到到初始接触点之间的垂直距离；

因为角接触的滚动体为钢球，所以滚动体和它的接触体主曲率重合；利用坐标变化，式(16)变换成标准椭圆方程，因此：

$z=\frac{1}{2}(\frac{1}{R_{11}}+\frac{1}{R_{21}})x^2+\frac{1}{2}(\frac{1}{R_{12}}+\frac{1}{R_{22}})y^2={\mathrm{Ax}}^2+{\mathrm{By}}^2---\left(17\right)$

式中：

$B-A=\frac{1}{2}{\{{[(\frac{1}{R_{12}}-\frac{1}{R_{21}})+(\frac{1}{R_{22}}-\frac{1}{R_{11}})]}^2-4(\frac{1}{R_{12}}-\frac{1}{R_{11}})(\frac{1}{R_{22}}-\frac{1}{R_{21}}){\sin }^2\mathrm{\&alpha;}\}}^{1/2}---\left(18\right)$

$B+A=\frac{1}{2}(\frac{1}{R_{11}}+\frac{1}{R_{12}}+\frac{1}{R_{21}}+\frac{1}{R_{22}})---\left(19\right)$

接触体表面的位移δ由两部分组成，分别是接触体的刚性位移δ₁与δ₂；ω₁(x,y)、ω₂(x,y)、z₁(x,y)、z₂(x,y)为原点以外的点相对于接触平面的位移；在接触区内，满足变形协调条件：

δ＝δ₁+δ₂＝ω₁+ω₂+z₁+z₂(20)

由式(18)～(20)得椭圆率参数e的方程：

$\frac{B-A}{B+A}=F\left(\mathrm{\&rho;}\right)=\frac{2(1-e^2)}{e^2}\frac{E\left(e\right)-K\left(e\right)}{E\left(e\right)}+1---\left(21\right)$

式中：

$a={\left(\frac{2E\left(e\right)}{\mathrm{\&pi;}(1-e^2)}\right)}^{1/3}{\left(\frac{3Q}{2\mathrm{\&Sigma;\&rho;}{E}^{\&prime;}}\right)}^{1/2}---\left(22\right)$

$b={\left[\frac{2\sqrt{1-e^2}E\left(e\right)}{\mathrm{\&pi;}}\right]}^{1/3}{\left(\frac{3Q}{2\mathrm{\&Sigma;\&rho;}{E}^{\&prime;}}\right)}^{1/3}---\left(23\right)$

$\mathrm{\&delta;}=\frac{2K\left(e\right)}{\mathrm{\&pi;}}{\left(\frac{(1-e^2)\mathrm{\&pi;}}{2E\left(e\right)}\right)}^{1/3}{\left(\frac{3Q}{2\mathrm{\&Sigma;\&rho;}{E}^{\&prime;}}\right)}^{2/3}\frac{\mathrm{\&Sigma;\&rho;}}{2}---\left(24\right)$

$\mathrm{\&Sigma;\&rho;}=\frac{1}{R_{11}}+\frac{1}{R_{12}}+\frac{1}{R_{21}}+\frac{1}{R_{22}}---\left(25\right)$

$\frac{1}{E^{\&prime;}}=\frac{1-v_1^2}{E_1}+\frac{1-v_2^2}{E_2}---\left(26\right)$

式中，δ为接触体表面的位移；Σρ为曲率和函数；v₁、E₁和v₂、E₂分别是物体V₁和V₂的泊松比和弹性模量；

S2.2建立主轴静止状态下角接触球轴承预紧模型；

主轴静止状态下，角接触球轴承受到轴向预紧力F_a，每一个滚动体将承受相同的载荷并产生相同的变形；设O_i，O_e为初始状态下内外圈沟道曲率中心，α和α′为初始接触角和实际接触角，Q_n为法向接触力；若轴承外圈固定，内圈受到轴向预紧力，内圈将产生轴向位移，内圈沟道曲率中心从O_i移动到O′_i；

由几何关系得到O_iO_e和O′_iO_e的线段长度：

O_iO_e＝(f_i+f_e-1)D_w(27)

$O_i^{\&prime;}{O}_e=\frac{d}{\cos {\mathrm{\&alpha;}}^{\&prime;}}=O_i{O}_e\frac{\cos \mathrm{\&alpha;}}{\cos {\mathrm{\&alpha;}}^{\&prime;}}---\left(28\right)$

式中，f_i与f_e分别为轴承内外圈沟道曲率系数；

线段O′_iO_e与O_iO_e的距离之差是滚动体与内外滚道的法向接触变形δ_n；

δ_n＝O′_iO_e-O_iO_e(29)

法向接触载荷为：

Q_n＝K_nδ_n(30)

式中，Z为钢球数目,K_n为钢度系数；

$K_n=2.1343*{10}^5({\mathrm{\&delta;}}_i^{*}{\mathrm{\&Sigma;\&rho;}}_i^{\frac{1}{3}}+{\mathrm{\&delta;}}_o^{*}{\mathrm{\&Sigma;\&rho;}}_o^{\frac{1}{3}})---\left(31\right)$

轴承的力平衡方程：

ZQ_nsinα′＝F_a(32)

将式(27)～(31)带入式(32)：

$\frac{F_a}{{\mathrm{ZK}}_n{O}_i{O_e}^{1.5}}=\sin {\mathrm{\&alpha;}}^{\&prime;}{(\frac{\cos \mathrm{\&alpha;}}{\cos {\mathrm{\&alpha;}}^{\&prime;}}-1)}^{1.5}---\left(33\right)$

由式(33)求得实际接触角α′；

钢球的曲率为：

$\frac{1}{R_{11}}=\frac{1}{R_{12}}=\frac{2}{D_w}---\left(34\right)$

对于内滚道，沟道曲率为：

$\frac{1}{R_{i1}}=-\frac{1}{f_i{D}_w}---\left(35\right)$

$\frac{1}{R_{i2}}=\frac{2\cos \mathrm{\&alpha;}}{d_m-D_w\cos \mathrm{\&alpha;}}---\left(36\right)$

对于外滚道，沟道曲率为：

$\frac{1}{R_{o1}}=-\frac{1}{f_e{D}_w}---\left(37\right)$

$\frac{1}{R_{o2}}=\frac{2\cos \mathrm{\&alpha;}}{d_m+D_w\cos \mathrm{\&alpha;}}---\left(38\right)$

将式(34)～(38)带入式(21)建立方程f₁,f₂：

$f_1=\frac{2(1-e_i^2)}{e_i^2}\frac{E\left(e_i\right)-K\left(e_i\right)}{E\left(e_i\right)}+1-\frac{B_i-A_i}{B_i+A_i}=0---\left(39\right)$

$f_2=\frac{2(1-e_o^2)}{e_o^2}\frac{E\left(e_o\right)-K\left(e_o\right)}{E\left(e_o\right)}+1-\frac{B_o-A_o}{B_o+A_o}=0---\left(40\right)$

S2.3建立主轴旋转状态下角接触球轴承定位预紧模型；

定位预紧的轴承在使用过程中，其内外圈相对位置是不会改变的；主轴高速旋转时，轴承受到的离心力F_ck和陀螺力矩M_gk的影响，滚动体会产生位移，O′为滚动体初始形心位置，O″为主轴旋转时滚动体实际形心位置，D为内圈沟道曲率中心，B为外圈沟道曲率中心；

主轴高速旋转时，轴承受到的离心力F_ck和陀螺力矩M_gk；

$F_{\mathrm{ck}}=\frac{1}{2}m{D}_w{\mathrm{\&Omega;}}^2{\left(\frac{{\mathrm{\&Omega;}}_E}{\mathrm{\&Omega;}}\right)}^2---\left(41\right)$

$M_{\mathrm{gk}}=J_b{\mathrm{\&Omega;}}^2{\left(\frac{{\mathrm{\&Omega;}}_E}{\mathrm{\&Omega;}}\right)}_k\sin {\mathrm{\&alpha;}}_k---\left(42\right)$

式中，Ω为主轴旋转角速度；Ω_E为滚动体绕着主轴公转角速度；Ω_B为滚动体自转角速度；J_b为滚动体转动惯量；

将离心力离心力F_ck和陀螺力矩M_gk按照轴向和径向进行分解得到平衡方程f₃,f₄：

$f_3=F_{\mathrm{ck}}+\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}+Q_{\mathrm{ik}}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}=0---\left(43\right)$

$f_4=Q_{\mathrm{ik}}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}+\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}=0---\left(44\right)$

由几何关系得到方程f₅,f₆：

f₅＝BO″*cosθ_ok+DO″*cosθ_ik-BD*cosα′＝0(45)

f₆＝BO″*sinθ_ok+DO″*sinθ_ik-BD*sinα′＝0(46)

式中，θ_ik与θ_ok分别为内外圈接触角，Q_ik与Q_ok分别为内外圈法向接触力；

DO″＝(f_i-0.5)D_w+δ_ik(47)

BO″＝(f_e-0.5)D_w+δ_ok(48)

S3.通过牛顿法求解出工况条件下最大接触应力与预紧力的关系；

将轴承内外圈椭圆率e_i，e_o，内外圈动态接触角θ_ik，θ_ok，内外圈法向接触力Q_ik，Q_ok这6个参数设为未知量，联立方程f₁到f₆得到非线性方程组；

$\left\{\begin{array}{c}{f}_1=\frac{2(1-e_i^2)}{e_i^2}\frac{E\left(e_i\right)-K\left(e_i\right)}{E\left(e_i\right)}+1-\frac{B_i-A_i}{B_i+A_i}\\ f_2=\frac{2(1-e_o^2)}{e_o^2}\frac{E\left(e_o\right)-K\left(e_o\right)}{E\left(e_o\right)}+1-\frac{B_o-A_o}{B_o+A_o}\\ f_3=F_{\mathrm{ck}}+\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}+Q_{\mathrm{ik}}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}\\ f_4=Q_{\mathrm{ik}}\sin {\mathrm{\&theta;}}_{\mathrm{ik}}+\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-Q_{\mathrm{ok}}\sin {\mathrm{\&theta;}}_{\mathrm{ok}}-\frac{M_{\mathrm{gk}}}{D}\cos {\mathrm{\&theta;}}_{\mathrm{ok}}\\ f_5={\mathrm{BO}}^{\&prime;\&prime;}*\cos {\mathrm{\&theta;}}_{\mathrm{ik}}+{\mathrm{DO}}^{\&prime;\&prime;}*\cos {\mathrm{\&theta;}}_{\mathrm{ik}}-\mathrm{BD}*{\cos \mathrm{\&alpha;}}^{\&prime;}\\ f_6={\mathrm{BO}}^{\&prime;\&prime;}*\sin {\mathrm{\&theta;}}_{\mathrm{ok}}+{\mathrm{DO}}^{\&prime;\&prime;}*\sin {\mathrm{\&theta;}}_{\mathrm{ik}}-\mathrm{BD}*\sin {\mathrm{\&alpha;}}^{\&prime;}\end{array}\right.---\left(49\right)$

应用MATLAB数值分析软件，通过牛顿法进行数值迭代，求解该方程组；

方程组的Jacobi矩阵：

$J\left(x\right)=\left[\begin{array}{cccccc}\frac{{\&PartialD;f}_1}{{\&PartialD;e}_i}& \frac{{\&PartialD;f}_1}{{\&PartialD;e}_o}& \frac{{\&PartialD;f}_1}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_1}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{{\&PartialD;f}_1}{{\&PartialD;Q}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_1}{{\&PartialD;Q}_{\mathrm{ok}}}\\ \frac{{\&PartialD;f}_2}{{\&PartialD;e}_i}& \frac{{\&PartialD;f}_2}{{\&PartialD;e}_o}& \frac{{\&PartialD;f}_2}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_2}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{{\&PartialD;f}_2}{{\&PartialD;Q}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_2}{{\&PartialD;Q}_{\mathrm{ok}}}\\ \frac{{\&PartialD;f}_3}{{\&PartialD;e}_i}& \frac{{\&PartialD;f}_3}{{\&PartialD;e}_o}& \frac{{\&PartialD;f}_3}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_3}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{{\&PartialD;f}_3}{{\&PartialD;Q}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_3}{{\&PartialD;Q}_{\mathrm{ok}}}\\ \frac{{\&PartialD;f}_4}{{\&PartialD;e}_i}& \frac{{\&PartialD;f}_4}{{\&PartialD;e}_o}& \frac{{\&PartialD;f}_4}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_4}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{{\&PartialD;f}_4}{{\&PartialD;Q}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_4}{{\&PartialD;Q}_{\mathrm{ok}}}\\ \frac{{\&PartialD;f}_5}{{\&PartialD;e}_i}& \frac{{\&PartialD;f}_5}{{\&PartialD;e}_o}& \frac{{\&PartialD;f}_5}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_5}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{{\&PartialD;f}_5}{{\&PartialD;Q}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_5}{{\&PartialD;Q}_{\mathrm{ok}}}\\ \frac{{\&PartialD;f}_6}{{\&PartialD;e}_i}& \frac{{\&PartialD;f}_6}{{\&PartialD;e}_o}& \frac{{\&PartialD;f}_6}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_6}{{\&PartialD;\mathrm{\&theta;}}_{\mathrm{ok}}}& \frac{{\&PartialD;f}_6}{{\&PartialD;Q}_{\mathrm{ik}}}& \frac{{\&PartialD;f}_6}{{\&PartialD;Q}_{\mathrm{ok}}}\end{array}\right]---\left(50\right)$

迭代矩阵为：

J(x^(k))d^(k)＝-F(x^(k))(51)

x^(k+1)＝x^(k)+d^(k)(52)

终止条件：

|x^(k+1)-x^(k)|<ε(53)

取初始矩阵x⁽⁰⁾，通过式(50)(51)(52)计算得到x⁽¹⁾，再将x⁽¹⁾带入式(50)(51)(52)计算得到x⁽²⁾，循环此计算，当第k+1次计算结果x^(k+1)与第k次计算结果x^(k)的差值矩阵小于给定的误差ε时，计算停止，认为矩阵x^(k+1)为所求结果；

再根据式(15)中最大接触应力P₀与内外圈法向接触力Q的关系，得到工况转速n下，轴承内外圈与滚动体最大接触应力P₀与预紧力F_a的关系；

S4.根据材料最大许用应力确定极限预紧力；

根据工况转速n下，轴承内外圈与滚动体最大接触应力P₀与预紧力F_a的关系，得到轴承内外圈与滚动体最大接触应力P₀为材料最大许用应力[σ]时对应的预紧力F_a，既极限预紧力F_amax。