0 引　言

向心关节轴承是一种球面滑动轴承，通过凸面内圈和凹面外圈的相互转动进行工作^[1]，主要承受径向载荷，因其承载能力大、耐冲击、工作平稳可靠等特点，被广泛的应用于工程液压油缸、锻压机床、工程机械等方面。从某型船舶液压传动设备承受重载实际工况入手，该船舶设备暴露于海水中，不断承受海流和波浪的水动力，进而传至关节轴承，所以向心关节轴承作为液压传动的关键部件，影响着机构及设备安全性，稳定性和磨损寿命。

协调接触（Hertz接触力模型）的应用前提是接触区域宽度要远远小于接触物体的宏观尺寸^[2]，适用于间隙大于运动副宏观尺寸约1%的大间隙接触或滚动接触碰撞过程描述，其接触应力是随着嵌入深度递减的二次函数，远离接触点的应力会迅速衰减，求解点接触和线接触问题具有良好的精确度；协调接触是指间隙一般仅为运动副宏观尺寸的0.01%～1%，甚至更小，运动副的接触面积接近于构件宏观尺寸，此时基于弹性半空间理论的Hertz假设不再适用。对于向心关节轴承的内外圈面的接触而言，属于球面协调接触。

对于承受重载，球面协调接触的向心关节轴承，其分布不均匀的接触应力，直接影响轴承磨损及寿命，甚至影响船用液压传动设备的安全和稳定性。因此，接触应力分析对轴承优化设计至关重要。本文以某型船用向心关节轴承为研究对象，为解决其受径向载荷时非完整球面协调接触问题，进行公式推导并建立小间隙高共形球面接触的近似通用模型，采用理论计算和有限元分析相结合的方法，验证接触应力分布规律，明确关键参数如径向间隙$\Delta $R，径向载荷Q，外圈半宽C对接触应力及接触区域的影响，为进一步开展大型关节轴承结构改进优化及更新换代提供参考意义。

1 关节轴承接触力学解析模型

向心关节轴承，即为球面滑动轴承（Spherical plain bearing），主要承受径向载荷，结构如图1所示，Q为向心关节轴承所受径向荷载，轴承内圈半径为$ {R}_{1} $，外圈半径为$ {R}_{2} $，外圈和内圈间的接触为同曲率协调接触。为使关节轴承能完成摆动、倾斜和旋转运动，内外圈之间具有游隙（包括初始制造游隙及使用磨损游隙）。

关节轴承在承受径向载荷时，接触面产生接触应力。接触应力是轴承设计过程中的一个重要依据，尤其是分布不均匀轴承接触应力直接影响轴承磨损及寿命。因此，接触应力分析对轴承优化设计及了解轴承使用状况至关重要。

1.1 完整球面协调接触力学解析模型

向心关节轴承承受径向载荷时，其内圈和外圈的相互作用是轴对称弹性体的共形协调接触^[3]，由于接触点产生局部变形而形成接触区截面半径为$ l_{\mathrm{max}} $的球冠，其球顶K为最初的接触点。协调接触中，半径$ {R}_{1} $的内圈与半径为$ {R}_{2} $的外圈均不能简化为半无限大弹性体，接触区域不再能简化成平面，而是以初始接触点K为中心的轴对称三维有限曲面，如图2所示。

其中，Q为径向载荷；$ {F}_{n} $为法向接触力，可分解为沿公法线方向（径向）的力$ {F}_{s} $和为沿水平方向的力$ {F}_{z} $；ψ为接触点和球心的连线与公法线（径向）之间的夹角。径向载荷为0时，两球体仅在K点接触，随着载荷的增加，接触区半径随之增大,其最大接触力仍发生在K点，接触应力分布式为：

$ P\left(l\right)={P}_{0}{\left(1-\frac{{l}^{2}}{{l}_{\mathrm {\max}}^{2}}\right)}^{n}。$

(1)

式中：$ {P}_{0} $为点K处的接触应力，MPa；$ l $为接触曲面上的点与z轴的水平距离，mm；n为压力分布指数，在非协调接触中，n一般取 0.5，但当接触区域半径$ l\mathrm{_{\max}} $与外圈半径之比大于 0.6 时，n随着$ l\mathrm{_{\max}} $的增加而增加，通过大量的有限元仿真计算^[4]，可以用以下指数关系表示：

$ n=0.75+0.9\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-10\left(1-\frac{{l}_{\mathrm {\max}}}{{R}_{2}-T}\right)\right]。$

(2)

为得到式(1)中 $ {P}_{0} $和$ {l}_{\max} $等2个未知参数，取接触曲面上某微小接触单元 dA 为研究对象，分析可得：

$ \left\{\begin{aligned} & \mathrm{d}F_n=2\text{π}lP\left(l\right){\mathrm{d}}l，\\ & \mathrm{d}F_Z=\mathrm{d}F_n\times \mathrm{cos}\phi=2\text{π}lP\left(l\right)\mathrm{cos}\phi {\mathrm{d}}l{\mathrm{d}}\phi。\end{aligned}\right. $

(3)

由此，径向载荷Q按照下式进行积分：

$ \begin{split} & Q=\int_{0 }^{l_{\mathrm{max}}}\int_{ 0}^{\mathrm{cos}^{-1}\frac{\sqrt{R_1^2-l_{\mathrm{max}}^2}}{R_1}}\mathrm{d}F_Z= \\ & \int_{ 0}^{l_{\mathrm{max}}}\int_{ 0}^{\mathrm{cos}^{-1}\frac{\sqrt{R_1^2-l_{\mathrm{max}}^2}}{R_1}}2\mathrm{{\text{π}}}lP\left(l\right)\mathrm{c}\mathrm{o}\mathrm{s}\phi\mathrm{d}l\mathrm{d}\phi。\end{split} $

(4)

联立式(1)和式(4)，可得接触曲面上的最大接触应力$ {P}_{0} $为：

$ P_0=(n+1)\frac{Q}{\mathrm{{\text{π}}}l_{\mathrm{max}}^2\times \mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{cos}^{-1}\dfrac{\sqrt{R_1^2-l_{\mathrm{max}}^2}}{R_1}\right)}。$

(5)

根据弹性半空间体受集中载荷理论求解得到未知参数$ l\mathrm{_{\max}} $的解析式：

$ l_{\mathrm{max}}=\sqrt[3]{\frac{4R_2R_1\left(k_1+k_2\right)K\left(n+\displaystyle\frac{1}{2}\right)\left(n+1\right)Q}{\mathrm{{\text{π}}}\left(R_2-R_1\right)}} 。$

(6)

式中：$ {k}_{1}=({1-{{\mu }_{1}}^{2}})/{\mathrm{{\text{π}} }{E}_{1}} $；$ {k}_{2}=({1-{{\mu }_{2}}^{2}})/{\mathrm{{\text{π}}}{E}_{2}} $，$ {\mu }_{1} $、$ {\mu }_{2} $及$ {E}_{1} $、$ {E}_{2} $分别为内圈，外圈泊松比及其弹性模量；$ K= ({\sqrt{\mathrm{{\text{π}} }}\Gamma(n+1)})/\left[{2\Gamma\left(n+\dfrac{3}{2}\right)}\right] $，$\Gamma $为Gamma函数，其函数表达式为$ \Gamma\left(x\right)=\displaystyle{\int }_{0}^{\infty }{\mathrm{e}}^{-t}{t}^{x-1}\mathrm{d}t $。

推导过程中加入了$ \mathrm{cos}\mathrm{\theta } $进行修正并采用级数分离变量进行近似，因此，结果存在一定的误差。为此，引入参数t和w对式(6)中的半径差$ \left({R}_{2}{-R}_{1}\right) $进行修正，修正后方程为：

$ {l}_{\mathrm{m}\mathrm{a}\mathrm{x}}=\sqrt[3]{\frac{4{R}_{2}{R}_{1}\left({k}_{1}+{k}_{2}\right)K\left(n+\displaystyle\frac{1}{2}\right)\left(n+1\right)Q}{\mathrm{{\text{π}} }[t\left({R}_{2}{-R}_{1}\right)+\omega ]}} 。$

(7)

1）当$ {R}_{1}={R}_{2} $，内外圈半径相同，随着载荷的增加，接触区半径随之增大，但仍不超过内圈（外圈）半径，即$ l_{\mathrm{max}}=R_2 $，解得：

$ \omega =\frac{4\left({k}_{1}+{k}_{2}\right)K\left(n+\displaystyle\frac{1}{2}\right)\left(n+1\right)Q}{\mathrm{{\text{π}} }{R}_{2}}。$

(8)

2）当$ {R}_{2}\gg {R}_{1} $，$ {l}_{\max}/{R}_{2}\to 0 $，接触退化为Hertz非协调接触，其接触区域半径如下：

$ {l}_{\mathrm{m}\mathrm{a}\mathrm{x}}=\sqrt[3]{\frac{3{R}_{2}{R}_{1}\left({{k}}_{1}+{{k}}_{2}\right)Q}{4\left({R}_{2}{+R}_{1}\right)}}。$

(9)

对比Hertz模型接触半径公式解得：

$ t=\frac{16K(n+1)\left(n+\displaystyle\frac{1}{2}\right)}{3{\mathrm{{\text{π}} }}^{2}} 。$

(10)

1.2 非完整球面协调接触力学解析模型

以上公式建立的前提是内外圈是完整的球面接触，而向心关节轴承的外圈宽度通常小于内圈宽度，当接触区半径超过轴承宽度时，接触区形状并不是完整的球面。此外，由于关节轴承外环中球面接触区附近的材料在实际应用中没有任何位移约束，构成自由边界。自由边界在接触变形过程中将产生自由边界效应，自由边界效应会在靠近结构边界的部位产生较大的应力梯度和应力集中，可能会导致结构或材料的失效。所以需进一步针对完整球面的球面协调模型进行修正。假定外圈宽度为2C，接触首先发生在中心点，径向载荷增加，接触区域随之增大，接触区域的投影半径将经历2个阶段：0＜$ l\mathrm{_{\max}} $≤C和$ l\mathrm{_{\max}} $＞C。

1）当0＜$ l\mathrm{_{\max}} $≤C时，接触区在垂直于径向载荷平面内的投影仍为完整圆形，轴承外圈宽度对接触区域影响不大，接触区的应力分布可以用完整球面接触解析模型表示。

2）当$ l\mathrm{_{\max}} $＞C时，接触区的应力分布将受到外圈宽度的限制，如图3所示。将受径向载荷大小为$ Q\mathrm{_{total}} $的完整球冠去除两个完全相同的边界接触区域，剩下的部分即为轴承外圈接触区。轴承在承受径向载荷时，外圈的接触应力分布用等效力法则求解。定义关节轴承接触区$ {S}_{0} $上承受的法向载荷大小定义为Q。

根据式(5)，完整球冠接触区上的接触应力分布可以表示为：

$ {P}_{o}=\frac{(n+1){Q}_\mathrm{total}}{{{\text{π}} }{{l}_\mathrm{max}}^{2}\times \mathrm{s}\mathrm{i}\mathrm{n}\left({\mathrm{cos}}^{-1}\dfrac{\sqrt{{R1}^{2}-{{l}_\mathrm{max}}^{2}}}{R1}\right)} 。$

(11)

对轴承触区域$ {S}_{0} $上的接触应力进行积分，可得到$ {S}_{0} $上承受的径向载荷之和Q_o：

$ \begin{aligned}Q= & \int_{ }^{ }\mathrm{d}F_Z=4\int_{ 0}^C\mathrm{d}x\int_{ 0}^{\sqrt{l_{\mathrm{max}}^2-x^2}}P(x,y)\mathrm{d}y= \\ & 4P_0\int_{ 0}^C\int_{ 0}^{\sqrt{l_{\mathrm{max}}^2-x^2}}\left(1-\frac{x^2+y^2}{l_{\mathrm{max}}^2}\right)^n\mathrm{d}y\end{aligned} $

(12)

结合式(10)和式(11)，由轴承接触区的径向力得到整体球面上承受的径向载荷为：

$ (n+1)\displaystyle\frac{Q_{\mathrm{total}}}{\mathrm{{\text{π}}}l_{\mathrm{max}}^2}=\frac{Q}{4u}。$

(13)

其中，令$ u=\displaystyle\int_{ 0}^C\displaystyle\int_{ 0}^{\sqrt{l_{\mathrm{max}}^2-x^2}}\left(1-\displaystyle\frac{x^2+y^2}{l_{\mathrm{max}}^2}\right)^n\mathrm{d}y $。

方鑫^[5]的原模型中，为了修正球面协调接触解析模型的引入的参数$ t $与$ \omega $分别为$ t={2}/{\mathrm{{\text{π}}}}+\left({l_{\mathrm{max}}}/{R_2}\right)^2 $与$ \omega = {3.8304\left({k}_{1}+{k}_{2}\right)KQ}/{\mathrm{{\text{π}} }{R}_{2}} $，将本文所改进的完整球面与非完整球面的协调接触模型以及方鑫提出的解析模型与有限元仿真结果的误差对比如表1所示。本文的改进模型与有限元仿真的误差在受径向载荷较小时较为波动，集中在10%之内，当径向载荷超过500 kN时，相对误差集中在1%左右，证明了改进的模型在受重载，大轴承的接触应力问题具有较强的适用性。

2 接触应力有限元仿真

借助Abaqus有限元仿真分析软件，验证上述解析模型建立和理论公式推导的合理性。以某船用向心关节轴承为例，在Abaqus中建立轴承物理模型，针对其内圈和外圈之间的面-面接触，分析研究接触应力的分布情况，主要包括最大接触应力及接触区域面积。轴承内圈材料为9Cr18，弹性模量203 GPa，泊松比0.28，密度7700 kg·m⁻³；外圈材料为17-4 PH，弹性模量210 GPa，泊松比0.24，密度7780 kg·m⁻³；加载板/芯轴材料为45#钢，弹性模量196 GPa，泊松比0.24，密度7780 kg·m⁻³。轴承几何结构如图1所示，其中内圈半径$ {R}_{1} $为120 mm，宽度$ {D}_{1} $为142 mm，外圈半径$ {R}_{2} $为150 mm，宽度$ {D}_{2} $为40 mm，芯轴直径185 mm。

2.1 仿真模型的建立

为分析向心关节轴承在恒定径向载荷作用下的接触压力分布以及验证前述的内，外圈协调接触解析模型的正确性，针对上述向心关节轴承建立有限元分析模型。为避免单元过度扭曲和减小计算误差，采用C3D8I单元类型，六面体（Hex）单元形状进行网格划分；为获得形状规则的网格，采用进阶算法（advancing front）和扫略网格（Sweep）的划分方法；在接触区域进行局部细化布种，保证计算精度的同时节约了计算成本。

在轴承加载板上表面施加沿y轴正方向的径向载荷Q，固定芯轴圆柱形外表面，分别在内圈内表面和芯轴圆柱形外表面，外圈外表面和加载板内表面建立2对绑定约束（Tie），在外圈内表面（设为主面）和内圈外表面（设为从面）建立面-面接触，采用摩擦系数表示接触面之间的摩擦特性，为利于计算收敛，设置自动过盈接触限度（automatic overclosure tolerance）^[6]。具体如图4所示。

2.2 接触应力分布对比及分析

通过Abaqus分析，径向载荷为3000 kN时的内圈外表面（从面）接触面的接触应力云图如图5所示。共建立了10种不同径向载荷下的有限元模型，得到不同载荷下轴承的接触状态，如表1所示。提取接触区域沿周向中线（x轴）的接触应力及理论计算的应力分布曲线如图6所示，提取轴向中线 （z轴）节点的接触应力分布曲线如图7所示。

接触压力分布是协调接触分析的关键。图5表明轴承在受径向载荷时接触应力的扩散状况，受外圈宽度制约，其扩散并非圆形而是椭圆形。由图6可知，在不同径向载荷下建立的解析计算及有限元仿真得到的轴承接触应力分布是高度吻合的，以及表1中最大接触应力$ {P}_{0} $的误差分析，证明了修正的解析模型在能够准确地计算球体与球壳小变形协调接触问题上具有通用性。且接触应力分布规律与式(5)呈同一规律，即随着与初始接触点K的距离增加，接触应力减小速率缓慢，当距离超过(2/3)$ {R}_{1} $后，其减小速率加快；随着径向载荷增加，接触应力增大，接触区域半径也增加，符合实际情况。图8描述了受外圈宽度的限制，宽度方向上接触应力的分布产生自由边界效应，导致外圈端面承受的接触应力较大，发生倾斜时，易达到极限而破坏，外圈宽度应合理选取。

可知，向心关节轴承在受不同径向载荷时，本文的改进模型及方鑫提出原模型的最大接触应力解析解与有限元仿真结果的相对误差，改进的模型误差在受径向载荷较小时较波动，集中在10%之内，当径向载荷超过500 kN时，相对误差集中在1%左右，证明了改进的模型在受重载，大轴承的接触应力问题具有较强的适用性。相较方鑫提出的原模型，其承受重载时，相对误差集中在30%～40%，说明其对于大轴承受重载的工况并不适用。

2.3 径向载荷对接触应力及接触区域半径的影响

由式(5)和式(6)分析，轴承最大接触应力P₀随着径向载荷Q增大而线性增大；而接触区域半径$ l\mathrm{_{\max}} $则与径向载荷Q是1/3幂函数，随着径向载荷Q增加$ l\mathrm{_{\max}} $也随之增加,但增长速度会逐渐变慢。图9为有限元仿真和解析模型的最大接触应力P₀与径向载荷Q之间的关系，呈线性变化；图10所示在有限元仿真结果中，$ l\mathrm{_{\max}} $与径向载荷Q之间的非线性关系，当载荷较小时，$ l\mathrm{_{\max}} $随着径向载荷Q增加的速度较快，即接触面积快速增加，而当Q超过某一值后，$ l\mathrm{_{\max}} $缓慢变化最终趋于外圈半径，同样很好地吻合了解析模型中$ l\mathrm{_{\max}} $与径向载荷Q的非线性特征。从轻载到重载较大范围内，$ l\mathrm{_{\max}} $都逼近外圈半径的值，也说明本文的第一个极限假设正确，此假设已在理论推导的过程中证明（当$\Delta $R = 0时，$ l\mathrm{_{\max}} $ ≡ $ {R}_{2} $）。

3 影响因素对接触应力和接触区域的影响

基于前文的向心关节轴承协调接触通用模型，能够分析物理参数$ {P}_{0} $、Q、 $ {l}_\mathrm{max} $，游隙ΔR(即为R₂−R₁)，外圈半宽C之间的相互关系，为下一步产品设计和接触磨损分析提供参考。下面分别描述了游隙ΔR和外圈半宽C与轴承的最大接触应力P₀及接触区域半径$ {l}_\mathrm{max} $的关系。

3.1 游隙ΔR

图11(a)为游隙ΔR与接触区域半径l_max的关系。不论重载或轻载，当ΔR = 0时，l_max都逼近R₂，再次说明协调接触解析模型的第一个极限假设是正确的（当ΔR = 0时，a ≡ R₂）。ΔR较小时，l_max随ΔR的增加而减小的速率较快^[7]，而ΔR超过0.2 mm后，l_max变化较缓；图11(b)描述了P₀在游隙为0.01～0.05 mm范围内较小，而超过0.1 mm后，P₀显著增大，为游隙为0.01 mm的1.5倍，所以控制向心关节轴承内外圈游隙≤0.05 mm为最佳。

3.2 外圈半宽C

图12(a)描述了外圈半宽C与$ {l}_\mathrm{max} $的关系，当载荷较大时，$ {l}_\mathrm{max} $随C的增加而减小，说明随着C增大，外圈宽度限制越小，轴向（z向）也能够承受接触应力，接触面逐渐由椭圆趋近于圆，接触面积不变的情况下，接触区域半径自然减小。图12(b)为C与$ {P}_{0} $的非线性关系，$ {P}_{0} $随着C的增加而急剧减小，当超过0.08 m后，$ {P}_{0} $缓慢变化并趋于某极限值。对于承受重载的轴承，保证C在2/3$ {R}_{1} $能够消除宽度的限制且保证倾斜时外圈的强度。图13所示同样证明了控制在C=2/3$ {R}_{1} $时，接触应力减小的效果显著。

4 结　语

采用修正的球面协调接触模型及原模型分别对向心关节轴承进行接触应力分析，通过有限元分析验证，对比改进模型与原模型与仿真模型的相对误差，表明新修正解析模型在某船用向心关节轴承承受重载球面协调接触问题中具有适用性。结论如下：

1）向心关节轴承通过有限元仿真得到的最大接触应力与径向载荷呈线性正相关，与接触区域半径呈非线性相关的关系，此与非完整球面协调接触力学模型的解析公式相吻合。

2）最大接触应力在游隙为0.01～0.05 mm范围内较小，超过0.1 mm后，$ {P}_{0} $与ΔR呈线性正相关，且显著增大, 游隙为0.1 mm时的接触应力是游隙为0.01 mm时的1.5倍，游隙为0.2 mm的接触应力则是其的3倍，所以控制内外圈游隙≤0.05 mm为最佳。

3）由于外圈半宽C限制，内外圈接触面积为椭圆，随着C增大而逐渐趋于正圆。最大接触应力$ {P}_{0} $与C呈非线性关系，当超过0.08 m（即(2/3)$ {R}_{1} $）后，$ {P}_{0} $缓慢变化并趋于某极限值，尤其受重载的轴承，保证外圈半宽在2/3$ {R}_{1} $能够消除其限制且保证倾斜时外圈的强度。