通用形式的连续性方程为:

∂ρ∂t+∇·ρν→=0(1)

式中: ν→为速度; ρ为润滑液密度; ∇为梯度算子。

动量守恒方程为:

∂∂tρν→+∇·ρν→ν→=-∇p+∇·τ-¯+ρg→+F→(2)

式中: p为压力; g→为重力加速度; F→为外部体积力; τ-¯为应力张量。

能量方程为:

∂∂tρE+∇·ν→ρE+p=∇·k∇T+τ-¯·ν→+Sh(3)

式中: k为导热系数; Sh为体积热源项; E的表达式如下:

E=h-pρ+ν22(4)

假设液相中存在足够的气核供气泡生成, 且此时液相中随流体运动的气泡与流体间无速度差异, 则可以用Rayleigh-Plesset方程描述:

RBd2RBdt2+32dRBdt2=PB-Pρl-4νlRBRB•-2σρlRB(5)

式中: RB为气泡半径; σ为液体表面张力; ρl为液相密度; PB为气泡表面压力; P为远场压力。

采用Singhal模型描述气穴生成和溃灭的过程, 考虑液相和汽相的两相连续性方程为:

∂∂tfνρ+∇·fνρVν→=∇·Γ∇fν+Re-Rc(6)

式中: fν为汽相质量分数; Γ为扩散系数; Re、Rc分别为汽化和液化相变率。

根据Singhal模型:

当 P≤Pv时, 有:

Re=Fvapmax1.0, κ1-fν-fgσρlρv23Pv-Pρl(7)

当 P> Pv时, 有:

Rc=Fcondmax1.0, κfvσρlρv23Pv-Pρl(8)

式中: Fvap取0.02; Fcond取0.01; κ为湍动能; fg为不可溶气体质量分数; ρv为汽相密度; Psat为饱和蒸汽压力; Pv为修正后的饱和蒸汽压力, 其计算公式为:

Pv=Psat+120.39ρκ(9)

仿真计算中的边界条件设置如下:将图1中的供油槽设置为压力入口, 润滑油入口压力和温度分别如表1所示; 将轴承间隙两侧端面设为压力出口, 表面压力为0; 将轴表面设置为旋转面(Moving wall), 由于轴高速旋转, 假设轴表面温度分度均匀, 温度值与文献[29, 30]中取值相同, 如表1所示; 轴瓦固定不动, 外表面热边界条件设置为第3类边界条件, 取值见表1。

采用ANSYS Fluent 13.0软件作为计算平台, 选用3D-Double Precision-Steady-Pressure Based求解器; 选用Mixture多相流模型; 由于本文算例的 Re数为47.12~94.24, 因此黏性流模型设置为层流; 考虑流体黏性生热; 润滑油黏度采用变物性设置, 黏度-温度关系采用与文献[29, 30]中相同的表达式; 迭代求解方法采用SIMPLE算法, 各物理量的空间离散均采用Second Order Upwind设置。求解过程中, 计算开始时需适当调小能量和汽相质量分数的松弛因子, 防止计算初期的发散, 待计算残差稳定后将松弛因子调回默认值。

由于润滑油膜各个方向上都包含流动-传热信息。而润滑油膜在油膜厚度方向上尺度很小, 一般为10~100 μ m, 但油膜周向和轴向的尺度都为数十毫米。因此, 仿真计算结果对网格质量非常敏感。在不遗漏各个方向上计算结果的基础上, 同时又能提高计算收敛速度、节省计算资源, 合理的网格划分变得非常重要。本文采用如表2所示的4种网格划分方案, 选取油膜内最大压力和最高温度值进行网格无关性分析。如图2所示, 油膜内最大压力和最高温度计算结果在网格划分方案3以后已趋于稳定。因此, 本文采用表2中网格划分方案3进行后续研究, 即油膜厚度方向布置18层网格; 轴向400层网格; 周向布置360层网格。

表2

Table 2

表2(Table 2)

表2 网格无关性分析 Table 2 Grid independence analysis 网格划分方案 网格数(径× 轴× 周) 1 6× 200× 180 2 12× 200× 180 3 18× 400× 360 4 24× 400× 360

表2 网格无关性分析 Table 2 Grid independence analysis

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 网格无关性分析Fig.2 Grid independence analysis

图3给出了采用不同数值仿真方法, 轴瓦中间截面周向压力分布(逆时针方向)计算结果与试验值的对比情况。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 轴瓦中间截面周向压力分布Fig.3 Circumferential pressure distribution on middle plane of bearing bush

图3中, Reynolds方法计算结果和试验数据均来自文献[29, 30]。从图3可以看出:传统Reynolds方法的计算结果, 压力最高点与试验值的误差分别为16.1%、3.6%和8.7%; CFD方法由于采用了黏-温变物性设置, 压力最高点计算结果与试验值误差分别为4.0%、4.1%和1.2%。同时, CFD方法的压力分布计算结果在总体趋势上也更接近试验值。

图4为采用不同数值仿真方法, 轴瓦中间截面周向温度分布(逆时针方向)计算结果与试验值的对比情况。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 轴瓦中间截面周向温度分布Fig.4 Circumferential temperature distribution on middle plane of bearing bush

图4中Reynolds方法计算结果和试验数据也均来自文献[29, 30]。从图4可以看出, CFD方法获得的结果无论是在最高温度点还是整体趋势上, 都与试验值更相符, 而Reynolds方法获得的最高温度点计算结果普遍低于试验值。这为上文压力分布的计算结果提供了解释。从表1中的黏-温特性值可以发现, 润滑油黏度对温度变化非常敏感, 而润滑油黏度的空间分布会极大地影响压力分布的计算结果。由于求解时采用CFD方法, 因此获得了更为精确的流动、传热、黏性产热计算结果, 进而获得了更精确的润滑油三维温度分布。与此同时, 通过Fluent求解器中的变物性设置, 将温度分布与黏度分布实时耦合, 进而获得了更精确的压力分布结果。上述分析也在表3中得以体现, CFD方法的压力、温度计算结果均与试验值符合良好。

表3

Table 3

表3(Table 3)

表3 仿真结果与试验值标准差对比 Table 3 Standard deviation of simulation results 运行工况 p试验值 T试验值 转速/(r· min-1) 载荷/N Reynolds CFD Reynolds CFD 2000 4000 0.19 0.058 1.171 1.134 4000 6000 0.186 0.139 2.282 1.913 4000 9000 0.142 0.106 3.874 1.483

表3 仿真结果与试验值标准差对比 Table 3 Standard deviation of simulation results

从图4中还可以看出, 采用CFD方法计算温度分布, 在压力建立段(即图1中几何结构收敛段)效果较好。而在几何结构发散段, CFD方法能够模拟出温度分布的骤降, 但过于剧烈, 后段出现较明显的偏差; 传统Reynolds方法的温度分布结果则显得平缓, 两者与试验值都有一定差异。

分析这一现象的原因可能是:后段润滑油膜内的气穴破裂, 压力骤降, 因此下游低温润滑油在背压作用下, 回流进入气穴溃灭形成的真空区域, 与气穴破裂液化后形成的少量液态润滑油混合。混合阶段中, 由于倒流的润滑油在质量上占主要权重且温度偏低, 因此出现了温度的骤降。在CFD仿真中, Singhal空化模型判断油膜内压力一旦低于润滑油的气化压力, 则气穴全部立即开始溃灭, 反映在图上就形成了一个明显的温度拐点。在实际过程中, 气穴的溃灭是相对渐进的过程, 因此温度分布试验值曲线呈缓慢光滑下降。传统方法完全不能模拟气穴生成和溃灭过程, 只是采用Reynolds边界条件或者Half-Sommerfeld边界条件简单地将负压置零, 发散段传热采用简化处理, 因此温度分布过于平缓, 与温度试验结果曲线的趋势存在较大差异。

总体上看, CFD方法获得的温度分布在趋势上与试验值更为接近, 在气穴处理以及发散段温度计算中也更接近物理实际, 与传统方法相比更具优势。但现有的空化模型在判断、模拟气穴溃灭过程这一点上仍有进一步改进的空间。

图5为基于CFD方法引入全空化模型后, 润滑油膜内汽相百分比的计算结果云图, 以及产生气穴的区域与油膜压力的对应关系。图5和图6中轴颈旋转方向均为逆时针方向。从图中可看出:在润滑油膜的发散段, 出现了明显的空化区域。传统Reynolds方法, 采用Reynolds边界条件或者Half-Sommerfeld边界条件简单地将负压置零, 并不直接模拟气穴的产生和溃灭过程, 因此无法得到类似图5的计算结果。同样地, 采用CFD方法但不引入全空化模型也无法获得丰富的空化区域三维计算云图。引入全空化模型的CFD方法可为进行径向滑动轴承各项设计参数、结构参数对气穴形成区域的影响规律等后续研究提供数值仿真手段, 并可以获得丰富的两相流计算结果, 相比于Reynolds方法和以往CFD方法的局限, 体现出明显的优势。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 计算结果云图Fig.5 Contours of simulation results

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 空化模型对温度分布的影响 3.4 空化模型对温度分布计算结果影响Fig.6 Influence of cavitation models on temperature distribution

在以往文献中^{[20, 22]}, 虽然已采用了CFD方法求解径向滑动轴承流动问题, 但仍采用简化方法处理空化现象, 具体做法是利用Fluent UDF中的DEFINE_ON_DEMAND命令对可能产生气穴的区域施加Half-Sommerfeld边界条件。但在文献中只展示了采用这种方法后, 压力分布计算结果与试验值的对比, 并未讨论其对温度分布的影响。

图6为不同空化模型对温度分布的影响, 其中, 轴承转速为2000 r/min, 载荷为4 kN。从图6可以看出, 不同空化现象处理方法对径向滑动轴承几何发散区域的温度分布计算结果有明显影响:①采用Half-Sommerfeld边界, 只对压力计算结果做置零处理, 因此发散区域是理应产生气穴的区域, 在数值计算处理中, 完全由出口和间隙侧面压力边界中的回流填充, 而回流的润滑油温度取决于人为设定值, 图中上方的蓝色区域即反映了这种情况。②采用Singhal模型处理空化, 由于考虑了实际生成气穴的物理过程, 抑制了上述不合理的润滑油回流问题。

因此, 若想要获取精确的径向滑动轴承温度分布计算结果, 需要考虑气穴的影响, 引入完整的空化模型, 仅采用简化的空化现象处理方法所获得的结果与实际情况存在较大差异。