随着各类旋转机械性能要求的提高，水润滑轴承以其润滑介质无污染、来源广泛、节省能源等优点，自上世纪70年代以来越来越受到人们的广泛关注^[1-6]。自上世纪80年代末期开始，前苏联对水润滑流体动静力轴承性能进行了深入的研究；同时，英国、德国和日本等国家，也在水润滑轴承的研究方面做了大量的工作，发明并使用了一些新的水润滑轴承。从上世纪80年代中期开始，我国相继开展了水润滑轴承的理论探索和试验研究工作^[5-7]，并相继在船用离心泵和轴流泵中采用了水润滑轴承。从上世纪90年代中期开始，重庆奔腾科技发展有限公司与重庆大学机械传动国家重点实验室合作，已开发了多种规格的水润滑复合橡胶轴承，相关研究成果逐步得到发表^[8-9]。

然而截至目前，水润滑轴承的研究还有许多问题有待解决。例如，现有水润滑轴承的研究还没有进行轴承―水介质―轴(轴承系统)的整体热弹流性能分析，以更真实地反映轴承系统的性能。在较大的载荷下，显著的轴承和轴变形(包括弹性变形和热变形)可能会发生。而现有的水润滑轴承研究，大都集中在轴承或轴单独元件的性能研究上。这种方法最大的缺点是，轴承、水介质和轴界面处的热边界条件(对轴承性能具有显著影响)的确定具有人为性，以至于仿真结果难以真实地反映该类轴承的性能，从而限定了其应用。

针对上述问题，作者通过把轴―水介质―轴承作为一个整体，来研究水润滑轴承系统的热弹流性能。这种方法的优点是，把轴、水介质和轴承界面处的通常难以确定的热边界条件转化为系统内部边界条件，从而避免了过多人为假设条件的引入，数值仿真结果从而更接近于轴承系统的真实工况下的性能。由此，将进行轴―水介质―轴承整体的热弹流三维有限元分析。在有限元分析时，使用了多目标优化法来加快轴承水膜提供的承载力和轴承外力的平衡，而弹性变形与热变形将使用影响系数法来求解，以缩短计算时间。

1 控制方程

由于水润滑轴承结构众多，被研究的水润滑轴承系统结构可简化为图 1所示。图 1中，r为径向坐标；ϕ为周向坐标；z为的轴向坐标。外力F施加在轴上，轴以角速度Ω周向旋转，从而可使轴与轴承之间水膜产生动压力以抵消施加在轴上的外载荷。

由于水润滑介质的低密度和低粘度，同油润滑轴承比较更容易发生变形和接触，因此需使用合理的方程进行它的性能研究。

1.1 雷诺方程

图 1所示的水润滑轴承的润滑机理，可用雷诺润滑方程来描述。假设水是牛顿液体、连续且各项同性，在忽略水介质的密度、体积力和惯性力作用下，描述水介质润滑的方程可写为

$ \begin{array}{*{20}{c}} {\frac{1}{{{R^2}}}\frac{\partial }{{\partial \varphi }}\left( {{\varphi _x}\frac{{{h^3}}}{{12\mu }}\frac{{\partial p}}{{\partial \varphi }}} \right) + \frac{\partial }{{\partial z}}\left( {{\varphi _z}\frac{{{h^3}}}{{12\mu }}\frac{{\partial p}}{{\partial z}}} \right) = }\\ {{\varphi _c}\frac{U}{2}\frac{{\partial h}}{{R\partial \varphi }} + \sigma \frac{U}{2}\frac{{\partial {\varphi _s}}}{{{R_i}\partial \varphi }}。} \end{array} $

(1)

式(1)中，R为轴承内径，φ为周向位置，z为轴向的宽度，U为轴周向速度，μ和ρ分别为水的粘度和密度，流量因子φ_x、φ_y以及接触因子φ_c表达式分别见文献[10-12]。式(1)中，水膜厚方程为

$ h = c + e\cos \varphi + \alpha \left( {z - 0.5} \right)\cos \left( {\varphi - \theta } \right) + \delta 。$

(2)

式(2)中，α为轴与轴承中心线间的夹角(错位角)，c为轴与轴承间的间隙，e和θ分别为轴的偏心距及偏位角，δ为轴或轴承的变形。式(1)中的粘度使用Roelands粘度公式来计算

$ \mu = {\mu _0}\exp \left\{ {{A_1}\left[ {{{\left( {1 + {A_2}p} \right)}^z}{{\left( {{A_3}T - {A_4}} \right)}^{ - {s_0}}} - 1} \right]} \right\}。$

(3)

式(3)中，μ₀和T₀分别为水的环境粘度和环境温度。常数A₁=lnμ₀+9.67，A₂=5.1×1.0e-9，A₃=1/(T₀-138)，A₄=138/(T₀-138)，z=α/(A₁A₂)，s₀=γ/(A₁A₃)。粘压系数α取为2.2×10^-8，粘温系数γ取为4.76×10^-2。

为求解方程(1)，使用下列压力边界条件

$ \left. \begin{array}{l} p = 0,z = 0,{L_z},\\ \frac{{\partial p}}{{\partial \varphi }} = 0,\varphi = 0,\vartheta 。\end{array} \right\} $

(4)

其中，L_Z为轴承的宽度，ϑ为油膜破裂处的周向位置。

1.2 能量方程

水润滑轴承的承载能力常因较低的水介质粘度而不足，进一步轴承和轴承间的接触常会发生，因而接触产生的热对轴承系统性能的影响应加以考虑。实践表明，水介质与轴和轴承交界处的热边界条件对轴承系统性能有明显影响。本研究中，把轴、水介质和轴承作为一个有机的整体来进行研究。这样，水介质与轴和轴承界面处的热边界条件便转化为系统内部边界条件，避免了通常把它们分开处理时，交界处的热边界难以确定的现象。当然，这种做法要求有限元分析时，需要大量的网格数，从而增大了计算量，但对于目前的计算机的计算能力，做到这点并不难。研究水润滑轴承热作用所使用的能量方程如下

$ \begin{array}{*{20}{c}} {\rho Cp\left[ {{V_r}\frac{{\partial T}}{{\partial r}} + {V_\theta }\frac{{\partial T}}{{R\partial \varphi }} + {V_z}\frac{{\partial T}}{{\partial z}}} \right] = }\\ {\frac{\partial }{{R\partial r}}\left( {r{k_r}\frac{{\partial T}}{{\partial r}}} \right) + \frac{\partial }{{R\partial \varphi }}\left( {{k_\theta }\frac{{\partial T}}{{R\partial \varphi }}} \right) + \frac{\partial }{{\partial z}}\left( {{k_z}\frac{{\partial T}}{{\partial z}}} \right) + \mathit{\Phi }。} \end{array} $

(5)

式(4)中，T为轴、轴承或水介质温度，ρ和C_p分别为它们的材料密度和比热，r为径向坐标，V_r、V_θ及V_z分别为它们半径方向、周向以及宽度方向的速度。k_r、k_θ和k_z分别为它们半径方向、周向和宽度方向的热传导系数。需要指出的是，对于轴及轴承而言，式(5)中的热源项Φ置为零；对于润滑剂而言，Φ可由水的粘度、厚度和压力等求得。

轴承的外表面采用了对流热边界条件

$ k\frac{{\partial T}}{{\partial n}} = - {h_h}\left( {T - {T_\infty }} \right)。$

(6)

式(6)中，k为热传导系数，即式(5)中的k_r、k_θ或k_z；h_h为换热系数，T_∞为环境温度。

1.3 变形方程

轴与轴承的变形分别由它们各自的径向弹性变形δ_E与热变形δ_T组成。它可统一表达为

$ \delta \left( {\theta ,z,\Delta T,p} \right) = {\delta _E}\left( {\theta ,z,p} \right) + {\delta _T}\left( {\theta ,z,\Delta T} \right)。$

(7)

旋转轴的温度可近似看成均匀，其热弹性变形通过式(8)求得

$ {\delta _T}\left( {{\theta _B},\Delta T} \right) = {\alpha _J} \cdot \Delta T \cdot R\left( {1 + \varepsilon \cos \left( {\varphi - \theta } \right)} \right)。$

(8)

式(8)中，α_J为轴热膨胀系数，ΔT是平均温升，ε(ε=e/c)为轴的偏心率。

这样，水润滑轴承热弹流性能求解的关键便集中在轴与轴承的弹性变形，以及轴承的热变形的计算。由于变形计算通常耗费大量计算时间，因此寻求快速的变形计算方法显得尤为重要。对于轴和轴承，它们的变形可用影响系数方得到。这种方法是变形计算(包括弹性变形和热变形分析)前，用其他的数值方法先求出计算变形的影响系数，然后在变形计算时，把求得的影响系数分别乘以压力或温升，便得到弹性变形或热变形。轴或轴承的表面径向弹性变形弹性影响系数，是指在固体内(文中为轴或轴承)指定点(θ_ξ, z_η)处的单位力，引起其附近点(θ_j, z_k)的表面法向弹性变形，假设它用τ_E(θ_j, z_k, θ_ξ, z_η)，在本研究中用有限元方法获得。那么，轴或轴承表面指定点(θ_ξ, z_η)处，由油膜动压作用提供的承载力w_h(θ_ξ, z_η)以及它们表面上微突体提供的承载力w_asp(θ_ξ, z_η)，引起其附近点(θ_j, z_k)的径向弹性变形便可通过式(9)获得

$ \begin{array}{l} {\delta _E}\left( {{\theta _j},{z_k}} \right) = \\ \;\;\;\sum\limits_\xi {\sum\limits_\eta {{\tau _E}\left( {{\theta _j},{z_k},{\theta _\xi },{z_\eta }} \right)\left[ {{w_h}\left( {{\theta _\xi },{z_\eta }} \right) + {w_{asp}}\left( {{\theta _\xi },{z_\eta }} \right)} \right]} } 。\end{array} $

(9)

同理，轴承热变形，可通过热影响系数方法求得。热影响系数是指在固体内(本文为轴承)指定点(r_ζ, θ_ξ, z_η)处的单位温升，引起其附近点(θ_j, z_k)的表面径向变形，假设它用τ_T(θ_j, z_k, r_ζ, θ_ξ, z_η)表示，在本研究中它也用有限元方法提前获得。那么，轴承表面指定点(r_ζ, θ_ξ, z_η)处的温升ΔT(r_ζ, θ_ξ, z_η)引起其附近点(θ_j, z_k)的热变形便可通过式(10)评估

$ {\delta _T}\left( {{\theta _j},{z_k}} \right) = \sum\limits_\zeta {\sum\limits_\xi {\sum\limits_\eta {{\tau _T}\left( {{\theta _j},{z_k},{r_\zeta },{\theta _\xi },{z_\eta }} \right)\Delta T\left( {{r_\zeta },{\theta _\xi },{z_\eta }} \right)} } } 。$

(10)

上述影响系数法的使用也出现在参考文献[13]中。

1.4 平衡方程

水润滑轴承系统的平衡，涉及到施加在轴上的外力、水膜提供的承载力的平衡，轴承稳态运行条件为

$ \left| {F - W} \right|/F \le \beta 。$

(11)

式(11)中：W为水膜提供的承载力，F为施加在轴上的外力，β为收敛精度，本研究取1.0e×10^-3。

由于常规的方法难以满足大载荷下水润滑轴承稳态运行条件，为加快收敛速度，采用了非线性优化方法中的多目标函数法, 而目标函数用有限单元法求得。现选取下列函数作为目标函数

$ f = \left| {1 - \frac{W}{F}} \right| - \beta 。$

(12)

这样，当f→0时，收敛条件方程(12)便得到满足，从而便能得到轴承稳定的性能解。

在求解方程(11)时, 通常选取对压力有显著影响的因素作为设计变量, 从以上的控制方程看, 偏心率ε对轴承性能有较大的影响, 故可选取它们作为设计变量，而约束条件的上下限可人为指定。值得说明的是, 用有限元法求得这种目标函数的值是近似值, 所以最后求得的设计变量并不一定是最优值, 但可借助这种方法来改善计算收敛速度。

2 数值方案

轴—水介质—轴承热弹流性能分析时，首先把方程(1)离散成有限元表达式，进一步在式(2)~(4)的基础上，使用八节点等参单元来求解方程(1)，从而获得压力。通过把得出的压力作为输入参数，利用热边界条件式(6)来求解能量方程(5)以获得温度。在先前求解得到的压力和温度，以及预先利用影响系数法求得的压力影响系数和温度影响系数的基础上，通过式(8)和式(9)便可获得轴或轴承的弹性变形和热变形，以及它们的温升。反复进行压力和温度的迭代计算，直到它们分别满足各自的收敛精度(分别为1‰和2.5‰)，然后进行平衡外力的计算。若计算得到的油膜动压承载力与微凸体承载力之总和W，与预先指定的外载荷F的差值(见式(11))，小于一定的收敛精度(本研究中为2‰)，则进行后处理操着，并结束计算。否则，重复以上计算，并通过式(12)来优化偏心率及油膜间隙，直到上述收敛标准均得到满足。上述求解流程见图 2。

3 计算结果与比较 3.1 程序代码的有效性

在使用流固耦合有限元代码研究水润滑轴承系统热弹流性能前，有必要验证作者编制的程序的有效性。验证时，在相同的几何和工况参数下，把由该程序计算和由文献[13]得到的温度结果进行了比较，这种比较是针对油膜轴承得到的。图 3给出了轴的转速Ω为1 000 r/min、施加在轴上的外载荷F为6 672 N时，程序与文献[13]计算得到的轴承前端面(z=0)的温度比较。可见，两种方法得到的温度结果相差很小。这表明，本程序能有效地进行轴承系统性能的分析。如今，相关程序已经交付给工程机械公司使用，并通过了该公司的实验验证。这样，通过本程序代码及前述公式，便可进行轴承—水介质—轴的热弹流性能分析。

3.2 热弹流性能分析

使用上述开发的程序可进行不同工况下的水润滑轴承系统的热弹流性能仿真研究。水润滑轴承系统性能仿真时，采用的轴承系统尺寸和工况输入参数见表 1。

图 4给出了轴的转速度Ω为166 r/min时不同F下的水膜的无量纲压力(P=pc²/(μ₀ΩR²))。图中，X(X=x/R)和Y(Y=y/L_z)为无量纲周向和轴向坐标，x、y为有量纲周向和轴向坐标(以下各图相同)。由图 4(a)和4(b)可以看出，由水膜产生的动压力区集中在轴承的底端(周向位置π处)附近。实际上，只有这种分布下的压力构成的水膜合力才能抵消外力。比较图 4(a)和4(b)可以看出，外载荷的增加会带来压力峰值增大。实际上由式(11)可以可知，要平衡增大的外载荷，水膜压力必须增大。

图 5给出了外载荷F为200 kN、轴的转速度Ω为166 r/min时的轴承内表面和轴的外表面无量纲弹性变形(δ_E/c)。从图 5可以看出，周向位置π处附近，轴承和轴的弹性变形比较严重，这是因为此处的压力较其他位置的压力大。这种弹性变形会导致轴和轴承间隙的增大，致使轴承承载力下降。因此，传统的水轴承性能研究中，没有考虑弹性变形显然是不可取的。

图 6给出了外载荷F为200 kN时，轴的转速Ω为166 r/min时的轴承的无量纲的热变形。从图 6可以看出，周向位置π处附近，轴承热变形((δ_T/c))比较严重。由图 5分析可知，此处的水膜压力较其他周向位置的压力大，进一步式(5)中的热源项Φ较大，因而轴承的温升较其他位置高，其热变形较大。热变形会导致轴和轴承间隙的减小，从而削弱弹性变形的作用。一般说来热变形小于弹性变形，(这点可由图 6和图 5(a)比较看出)，热变形的削弱作用是有限的。因此，在水润滑轴承研究中，也应考虑热变形的作用。

图 7给出了外载荷F为200 kN时的轴承与轴间接触压力(P=p_ac²/(μ₀ΩR²))。其中，p_a为接触压力，其具体求解方法可参见文献[14]。从图 7可以看出，在200 kN这样大的外载荷下，轴承与轴间的接触会发生。同时可以看出，接触压力主要集中在周向位置π处附近，这是因为此处的压力较其他位置的压力大的缘故。这从理论上解释了传动的水润滑轴承容易发生接触疲劳失效的原因。而且可以看出，接触压力主要集中在轴端(y=0)处。这是因为计算时此处给定了一定的错位角(见表 1)。进一步比较7(a)和7(b)可以看出，轴转速的增加会使接触压力增大。这种高的接触压力的产生可能带来轴和轴承的疲劳失效，因此在水润滑轴承设计时应当加以考虑。

水润滑系统的温升是另一个关心的重点。图 8给出了不同无量纲油膜厚度下的轴承内表面的温度。由图 8可以看出，在膜厚比h/σ小于3下，轴承温升均很显著，尤其以周向位置π附近(油膜破裂处)轴承内表面温升最显著。这是因为在这样小的膜厚比下，大量摩擦产生的热还不能由轴与水介质运动充分带走。较高的温升会导致轴承系统的热失效，因此，水润滑轴承系统性能研究时应考虑温升的变化，特别是小膜厚比下的温升。需要指出的是，若轴的转速更高或外载荷更大，弹性变形、热变形以及温升将更加明显。

4 结论

1) 使用自编的流固耦合轴承有限元程序，通过影响系数法、有限元法和非线性优化方法，对某水润滑轴承—介质—轴的热弹流性能进行了三维有限元分析。同时，程序的有效性被证实。

2) 一定工况下的水润滑轴承会存在明显的弹性变形、热变形以及一定程度接触压力。因此，这些变形在水润滑轴承性能研究中应当考虑。

3) 小膜厚比下水润滑轴承会出现较大的温升(尤其在周向油膜破裂处)，从而有可能引起水润滑轴承润滑失效。因此，进一步研究如何通过轴承设计参数的改变来控制其温升，是水润滑轴承研究必要的工作。