一种滑动轴承油膜性能计算的动网格方法

2020-01-03尹雪梅李蒙蒙李奕君

中国机械工程 2019年24期

吴超尹雪梅李蒙蒙李奕君王文

1.郑州轻工业大学机电工程学院，郑州，4500022.郑州轻工业大学能源与动力工程学院，郑州，4500023.上海大学机电工程与自动化学院，上海，200072

0 引言

计算油膜轴承性能的关键问题就是求解动态流体润滑方程，得到油膜的压力分布[1-4]。通常情况下，采用动态雷诺方程计算轴承的性能时，难以精确反映转速所引起的油流周向惯性效应、动态挤压效应和静压效应之间的线性耦合关系及其对油膜三维压力场、温度场和速度场的影响，因此，有必要直接求解Navier-Stokes方程，精确研究轴承参数对油膜性能的影响[5]。文献[6-8]采用计算流体动力学(computational fluid dynamics，CFD)开展了静压轴承承载特性的研究，证实了在表征复杂求解域流体流动形态方面Navier-Stokes方程可弥补雷诺方程的不足。文献[9-10]利用Fluent软件，采用静网格方法，研究了气穴现象对滑动轴承性能的影响，但静网格方法不能计算变载荷下轴心的轨迹。文献[11-12]将轴颈的旋转动边界转换为静边界，并利用Smoothing动网格模型，成功计算了油膜轴承的动特性系数。文献[13-15]自定义动网格更新程序，提出了基于弹性变形的动网格调整法，能用来求解瞬态轴心轨迹，并利用弱耦合算法研究了轴承刚度随方向的变化以及涡动中心与载荷、不平衡量的关系。文献[16]针对求解复杂转子-轴承系统非线性动力学特性的问题，基于Smoothing动网格技术，提出了一种计算流体力学和计算转子动力学的流固耦合新方法，计算结果表明，该方法能够得到精确的轴心轨迹，并能准确分析复杂转子-轴承系统非线性动力学特性。文献[17]利用非定常动网格技术建立了考虑轴颈涡动频率与涡动轨迹的滑动轴承动力特性求解模型，研究了不同的轨迹下轴颈涡动频率和偏心率对滑动轴承动力特性的影响。以上在油膜轴承性能计算中运用的动网格方法是基于Smoothing模型或其变形模型基础上进行的，能够方便地求解出油膜轴承的特性参数，分析轴承-转子系统的动力学性能。

Smoothing动网格模型[18]的特点是不改变网格节点间的拓扑关系，只进行网格形状的改变，能够计算变载荷下轴承的轴心轨迹，但Smoothing模型网格位移的更新方法容易造成网格畸变，当偏心率过大时，垂直于轴颈的网格线就会出现严重倾斜，甚至出现负网格的情况，导致计算发散。鉴于Fluent软件在表征流体流动形态方面的优势及其计算油膜轴承动态性能方面的不足，本文提出了一种基于Fluent方法的滑动轴承油膜性能计算的动网格更新方法。在它的每一时间步内，所有网格节点更新方式是依据轴颈中心坐标确定的，使径向网格线始终垂直于轴颈表面，保证了移动前后膜厚方向的网格线沿圆周方向上等均分布，网格不发生扭曲变形，避免了网格节点更新产生的累计误差。通过与典型算例和实验结果对比分析来验证该方法的有效性，并考察了进油压力和载荷对转子静平衡位置的影响。

1 动网格更新算法原理

1.1 已有动网格方法

动网格模型可以用来模拟由于边界运动而引起的流域随时间的变化。边界运动可以是主动的规定动作(例如，用户可以定义刚体中心的速度和角速度)，也可以是被动的待求解的动作(例如，六自由度模型中，已知运动状态的刚体边界，受合力后确定下一时刻的速度和加速度)。根据每一个时间步的运动边界的位置，Fluent软件自动对动网格更新处理。对于动网格模型，用户需要确定初始网格和运动区域。3种常用的动网格更新算法[18]分别是Smoothing、Layering和Remeshing。Smoothing方法是不改变网格节点间的拓扑关系，只进行网格形状的改变，但当网格出现畸变时计算容易发散。Layering方法是随着动边界的移动，在边界处发生网格的增加或合并。Remeshing方法是将控制区内的所有网格重新划分，一般适用于非结构化网格。理论上，考虑到油膜几何尺寸的特点，应当选择Smoothing动网格更新方式，但在轴承性能计算过程中也出现了网格畸变，显然Smoothing方法也不能很好地解决轴承油膜的瞬态计算问题。

文献[13]提出的基于Smoothing的网格更新模型示意图见图1。网格节点位移按公式Δxi=niΔx/n和 Δyi=niΔy/n进行计算(n为膜厚方向网格的层数，ni为网格节点所在的层数，Δx、Δy为轴颈中心在该时间步的位移)。利用该方法，垂直于轴颈表面的网格线就会出现倾斜，导致计算结果误差较大，特别是多步计算后的位置累计误差。

图1 文献[13]网格更新方式Fig.1 Mesh updating method in references

1.2 自定义动网格更新算法原理

为解决上述难点，本文提出了一种新的网格更新算法，在网格更新前后，膜厚方向的网格线始终指向轴颈的中心，从而避免了出现网格倾斜。本文网格结构更新算法的示意图见图2，网格节点位置计算原理见图3。

图2 本文提出的新的网格更新方式Fig.2 New mesh updating method in this paper

图3 网格节点更新原理图Fig.3 Schematic diagram of mesh node updating

为使网格不发生扭曲变形，就需要指向轴颈中心的网格线方向始终不改变，即始终垂直于轴颈表面，同时要保证在移动前后膜厚方向的网格线沿圆周方向上等均分布(膜厚网格线的夹角相等)，即更新前后网格线平行，也就是图3中的向量O1N1与O2N2共线。具体方法如下。

记轴承的中心为O，移动前后的轴颈的中心分别为O1和O2，油膜网格在膜厚方向上分成k层网格，假设N1为t1时刻网格上任意一个待更新节点，节点N1以及与其相对应的更新后的节点N2均在第ki层网格上，根据此刻轴颈所受的油膜力、轴颈坐标、轴颈中心位移速度、轴颈质量可以得到下一时刻t2的轴颈中心坐标，根据时刻t2的轴颈中心坐标和该节点相对轴心的单位向量(与向量O1N1平行或者共线)，作射线分别交轴颈圆和轴承圆于点A2、B2。又知待求点N2所在的层数为ki，即已知LA2N2=(ki/k)LA2B2，LOB2为轴承半径，LO2A2为轴颈半径。N2坐标求解算法步骤如下。

(1)求出LOO2；

(2)已知向量O1N1的单位法向量e，由于向量O1N1与向量O2N2共线，故可得到向量O2N2的单位法向量e；

(3)求解向量O2O与向量O2B2的夹角余弦cosα；

(4)利用三角形关系式求解LO2B2：

(5)求出LA2B2=LO2B2-LO2A2；

已知轴颈中心移动后的坐标O2，根据上述算法可求出节点位置N1移动到新位置的坐标N2。

2 轴心静平衡位置的求解方法

2.1 数学模型

根据牛顿第二定律，运动方程[13]为

(1)

式中，m为转子质量的一半；S为转子中心距原点的位移，G=mg；F为转子受的油膜合力。

已知时间步Δt，上一时间步对应的速度v0，根据牛顿第二定律求解每一个时间点轴颈中心的相对位移和速度计算公式如下：

(2)

v=v0+(G+F)Δt/m

(3)

求解油膜力F的分量Fx和Fy的公式如下：

(4)

(5)

式中，Fx、Fy分别为x轴和y轴方向的油膜力；A为某一网格单元的面积；p为根据N-S方程计算出的轴承油膜各点的压力。

2.2 求解轴心静平衡位置流程及宏函数的选用

基于Fluent软件，利用该方法计算油膜轴承所支撑轴颈静平衡位置的流程如图4所示，首先设定轴颈初始速度和初始位移，通过udf宏DEFINE_GRID_MOTION计算出油膜力Fx和Fy，再根据式(2)和式(3)计算出轴心坐标和速度。同时按照新提出的算法原理进行网格更新，在该时间步内迭代完后进入下一个时间步。下一个时间步用到上一时间步的轴心坐标和速度(读取储存于文件2和文件3最后一次更新的数值)。一直循环下去，直到轴颈中心趋近于一定点，即可停止计算，该定点即稳态下径向油膜轴承所支撑轴颈中心的静平衡位置。

图4 轴心静平衡位置计算流程图Fig.4 Computational flow chart about static balance position of the rotor

一般来说，动网格节点的位移由坐标确定，但滑动轴承油膜是一种长径比很大的网格，用坐标判断节点位置容易出现错误，可采用节点全局编号判定节点的位置，强制精确控制每一个节点的位移变化。利用DEFINE_GRID_MOTION宏对网格节点位置进行强制性精确定义，可有效降低网格畸变的可能性，由于轴颈表面各点的速度不一样，可通过编写DEFINE_PROFILE宏程序控制轴颈表面各点的速度。按照以上步骤，可以利用动网格方法实现对滑动轴承油膜状态的仿真。

3 算法的正确性和有效性验证

本文提出的网格更新方法与文献[13]所使用的网格更新方法的主要区别就是更新后油膜厚度方向的网格线是否垂直于轴颈表面。要验证所提动网格算法的正确性和有效性，需要进行两方面的验证：用于静态油膜力计算的精度验证；用于轴颈中心静平衡位置求解的累积误差验证。

3.1 算例1：用于静态油膜力计算的精度验证(单次油膜力计算)

利用所提出的网格结构和文献[13-15]的网格结构，分别计算轴颈转动角速度为500 rad/s、偏心率为0.5和0.9的圆柱轴承的油膜力，来对比这两种网格在静态油膜力计算(用于静网格计算时，单次网格计算)的精度。算例选取的轴承几何参数和润滑油物性参数见表1。轴承供油方式采用两侧双向进油并且设置轴向油槽。

由于油膜厚度尺寸很小，故油膜厚度方向的网格划分对计算结果影响最大，需要先进行网格独立性验证。网格独立性计算时选择多项流混合模型，气穴模型选用Singhal全空化模型，连续性方程、动量方程和气穴方程均选用一阶迎风格式进行求解，速度耦合格式选择SIMPILE格式。轴

表1 轴承和润滑油物性参数

承油膜网格初始划分的网格数(膜厚方向×周向×轴向)为5×600×80，通过改变膜厚方向的网格层数，并以净流量小于进口和出口流量两者之间的最小值的1%作为判断收敛的标准。膜厚方向上设置不同的网格层数，以承载力作为观测量，网格独立性验证结果见表2。结果表明，本文所用的网格密度计算结果的最大偏差均小于4%。

表2 网格独立性验证结果

综合考虑计算效率和计算精度，把油膜膜厚方向网格划分为6层(6×600×80)，见图5。算例中压力-速度耦合格式分别选SIMPILE、SIMPILEC、PISO方式，其他设置同上。采用瞬态计算方法，利用文献中的网格方法和本文的网格方法计算得到的承载力结果和圆周方向的压力分布结果(取轴承宽度中间L/2处，偏心率为0.9)分别见表3和图6。

图5 油膜网格Fig.5 Oil film mesh

表3 两种静网格结构计算的承载力对比

(b) 本文方法的计算结果

由表3和图6可以看出，基于两种网格方法得到的轴承承载力和瓦块圆周方向的压力计算结果相同，均与文献[19]结果相接近，表明两种网格(油膜厚度方向的网格线垂直于或者不垂直于轴颈表面)用于轴承油膜稳态油膜力计算时，均能满足精度要求(<9%)。由于在计算过程中网格不更新(静网格法)，无论油膜厚度方向的网格线是否倾斜，累计计算误差都不存在，故两种网格的计算结果基本无差别。这就证明了Smoothing网格方法和本文提出的网格方法来计算轴承稳态油膜力(静网格法)的精度是满足要求的。

3.2 算例2：用于静平衡位置求解的累计误差验证(多次油膜力计算)

根据文献[13]的算例，选取轴颈角速度为500 rad/s，转子质量为17.527 kg(171.76 N)，轴承、润滑油的参数和油膜网格划分同算例1，对轴颈中心轨迹的收敛过程及静平衡位置进行了仿真，结果见图7，计算得到其偏心率为0.221。根据文献[19](未给出进油压力，其他条件完全一样)，在相同轴承参数下，据此载荷计算出轴承偏心率为0.267，与本文计算结果相比，相对误差是17.2%。文献[13]和本文在计算时考虑了气穴的影响，而文献[19]没有考虑气穴的影响。

(a) 二维图

(b) 三维图

文献[13]的静平衡位置所对应的偏心率为0.046。将文献[13]结果和本文的计算结果分别与文献[18]相比较，文献[13]的计算结果差别更大，偏心率明显不在同一数量级，表明文献[13]的网格经多次更新后，出现了较大的累计误差。这与文献[13]在网格更新时沿膜厚方向的网格线与轴颈表面不垂直有关，也与其网格节点坐标更新的计算方式与轴颈表面的位置有关，从而导致多次迭代后求解轴颈中心静平衡位置时出现了较大的误差。

为验证所编制程序的稳定性，从3个不同起始位置开始分别计算相同条件下的轴心静平衡位置，轴心轨迹收敛曲线如图8所示。可以看出，3条曲线最终收敛于一点，表明计算的起始位置对轴颈的静平衡位置并无影响，验证了所编制的网格更新程序的稳定性。

图8 不同起点位置的轴心收敛轨迹Fig.8 Axial convergence trajectory of different starting positions

3.3 算例3：进油压力对轴心静平衡位置的影响

为了分析进油压力对轴承所支撑转子轴心静平衡位置的影响，在给定气化压力29 185 Pa、转子质量17.527 kg、角速度500 rad/s和 000 rad/s的条件下，计算了供油压力在0.1～0.5 MPa之间时转子轴心的静平衡位置。轴承的偏心率和偏位角随进油压力的变化关系如图9所示。算例的其他计算参数见表1。

(a) 偏心率

(b) 偏位角

由图9可以看出，随着进油压力的增大，偏心率减小(<7%)，偏位角增大(<9%)。这是因为在载荷和转速不变的前提下，随着进油压力的增大，轴承油膜的最大压力增大，承载能力增强，偏心率减小。同时由于进油口在水平方向上，故随着进油压力增大，偏位角增大。计算结果表明进油压力对轴承的承载性能有影响，但影响率小于10%。结合算例2的结果，进一步证明了所提出的动网格更新方法用于轴心轨迹计算时累计误差相对较小。

通过以上算例验证了所提出的动网格方法在油膜性能计算的有效性、可行性和稳定性。