滑动轴承动力特性的数值计算方法

李 强,许伟伟,王振波,金有海,郑水英

(1.中国石油大学化学工程学院,山东青岛 266580;2.中国石油大学储运与建筑工程学院,山东青岛 266580; 3.浙江大学化工机械研究所,浙江杭州 310027)

滑动轴承动力特性的数值计算方法

李 强1,许伟伟2,王振波1,金有海1,郑水英3

(1.中国石油大学化学工程学院,山东青岛 266580;2.中国石油大学储运与建筑工程学院,山东青岛 266580; 3.浙江大学化工机械研究所,浙江杭州 310027)

将计算流体动力学与简谐激励法应用于滑动轴承动力特性系数的求解,通过采用全新的变流域动网格技术提出一种基于瞬态流场计算的滑动轴承动特性的计算方法。利用所提出的方法计算典型滑动轴承的刚度、阻尼系数,并与已有的经典计算结果进行比较。结果表明:两种方法的计算结果基本吻合,提出的新方法不仅考虑了计算初值的影响,而且计入了实际存在的油膜破裂现象,满足精度要求,有效可行。

滑动轴承;计算流体力学(CFD);动特性系数;动网格;数值计算

滑动轴承广泛地应用于大型、高速旋转机械,其油膜力既是转子－轴承系统阻尼的主要来源,也是导致机组稳定性下降的重要原因,滑动轴承动态特性对转子振动有重要影响。能精确描述滑动轴承动力性能的数学模型是非线性的,在非线性模型还不完善的前提下,在线性假定的基础上采用4个刚度系数和4个阻尼系数表示小扰动下油膜的动态特性比较准确,可以满足通常的工程实际需要。为此国内外学者从理论和实验两方面做了大量工作[1-5]。通过计算流体力学方法(CFD)直接求解N-S方程来分析滑动轴承的3D润滑流场,该方法便于求解复杂区域上的问题,对于事先未知的自由边界或求解区域内部不同介质的交界面比较容易处理,因此基于CFD技术的3D流场计算方法广泛地应用于滑动轴承性能的研究[6-9]。考虑到现有动网格方法在计算滑动轴承瞬态流场时网格会出现较大的畸变,笔者利用CFD-FLUENT软件对液体动压滑动轴承进行建模,采用全新的变流域动网格技术对滑动轴承内瞬态流场进行计算分析,提出一种滑动轴承动特性系数的计算方法,并将计算结果和已有的数值计算结果进行对比,验证提出方法的有效性和可行性。

1 滑动轴承CFD计算模型

实际工作时滑动轴承在轴承开扩区有油和气混合存在,因此为了更准确地描述滑动轴承动力特性,本文建立的CFD计算模型可以考虑两相流对滑动轴承性能的影响。

1.1 控制方程

滑动轴承中润滑流场可以用质量、动量守恒方程和空穴边界方程来描述。

在边界运动的任意控制体V内积分形式的非定常不可压黏性流的连续性方程和动量方程的通式为

式中,ρm为气液两相的混合密度,kg/m3;φ为通量; dV为控制体V的边界;vm为气液两相的速度矢量;vs为动网格的网格变形速度矢量;Γ为扩散系数;Sφ为通量φ的源项。

为了充分地考虑实际运行中滑动轴承内润滑油的分离和重新合并的情况,本文中选用的空穴模型[10]是基于质量守恒边界条件的,

式中,vv为气相速度矢量;f为气相质量分数;γ为有效交换系数;Re和Rc分别为空穴的生成和凝聚率。

1.2 建模及网格划分

以圆柱滑动轴承为研究对象,图1给出了滑动轴承的3D网格结构。轴承水平中分面两侧各有1个进油槽,润滑油由两侧进油口进入轴承间隙内,由轴向两端流出。参照文献[11],轴承直径D=50 mm,轴承宽度B=25 mm,轴颈半径间隙c=0.05 mm,油槽包角α=30°,油槽轴向长度b=3.75 mm,轴颈偏心率ε=0.5,转速Ω=1 000 rad/s,润滑油的动力黏度为μ=1.25×10－2N·s/m2,气态润滑油参数取空气参数,进油压力p=103 kPa,流体流动状态为层流。

利用前处理软件Gambit对滑动轴承进行网格划分,由于滑动轴承润滑流场在空间3个方向的尺寸大小相差较大,须采用结构化网格来提高数值计算精度和节省计算时间。对不同网格密度下CFD计算模型的计算结果进行比较分析后发现,油膜间隙处不同的网格数导致计算结果差别较大;而周向和轴向网格密度对计算结果影响不大,因此这两个方向的网格尺寸主要从长宽比角度考虑(最佳长宽比为1∶1)。选用轴承间隙径向5层网格,轴向和周向网格密度为0.2。

图1 滑动轴承网格模型Fig.1 Grid structure of journal bearing model

1.3 FLUENT计算参数设置

对于滑动轴承的润滑流场,选择3D耦合求解器,隐式格式求解;流场计算的边界条件要用到进口、出口和壁面3种边界条件。滑动轴承的进口为压力边界条件,给定进口总压力;滑动轴承的出口同样为压力边界条件,压力根据实际情况设定为环境大气压;轴瓦表面为固定无滑移边界,近壁面应用标准壁面函数,轴颈表面绕偏心位置旋转。

采用有限体积法离散控制方程、连续方程,动量方程采用一阶迎风格式,压力差值格式采用PRESTO格式,压力速度耦合采用SIMPLEC算法。计算时所有方程的残差都小于1×10－4,并在计算过程中对进、出口流体质量流量进行监控,数值基本相等时计算收敛。

2 滑动轴承动力特性的计算

2.1 动网格方法

由于轴承径向间隙与其他方向的结构尺寸相差比较大,利用FLUENT软件自带的动网格方法进行网格更新将产生较大的网格畸变,影响瞬态计算的准确性,因此在进行伴随轴颈移动的网格更新时采用自行开发的变流域动网格技术。首先将轴承间隙流场用结构化网格进行划分,当轴颈移动时,轴承间隙流场最内圈上的所有节点都随轴颈一起移动,最外圈和油槽上的节点保持不动,中间的节点按照结构化网格划分原理进行移动,因此可以通过编程计算出每个节点坐标的变化值,然后通过FLUENT的UDF接口将各个节点移至新的位置。图2给出了轴颈移动到一定位置后轴承网格变化情况(图中油膜间隙被大大地夸大)。从图2可以看出,通过该方法移动后的网格即使在很大的轴颈偏心下依然能保证良好的光滑性和规整性,而且基本不会出现负体积。

轴颈中心的位置随着轴颈的移动不断变化,因此轴颈表面的旋转速度不能按常规的绕固定点旋转的方法设定,本文中通过FLUENT－UDF程序捕捉到每一步的轴颈中心,然后将旋转速度分解为x、y两个方向的速度施加到轴颈表面。

图2 轴颈移动后网格节点截面图Fig.2 Front view of grids in moved journal bearing

2.2 刚度阻尼系数计算方法

轴颈在静平衡位置受到外载荷作用时,对轴颈油膜产生位移和速度扰动,油膜力将发生变化,力的变化和扰动之间的关系一般是非线性的,为简化分析,当扰动较小时,常将这种关系线性化,如图3所示。油膜力可近似由下列线性关系表达:

式中,Fx、Fy为油膜力在x、y方向的分量;Fx0、Fy0为静态平衡位置时,油膜力在x、y方向的分量;kij、cij分别为轴颈中心在静平衡位置(x0,y0)时油膜的刚度(N/m)和阻尼系数,N·s/m。

由式(3)可以得到油膜力的增量,即动态油膜力,写成矩阵的形式为

式中,[K]为刚度系数矩阵;[C]为阻尼系数矩阵。

图3 动力特性系数示意图Fig.3 Schematic diagram of dynamic coefficients

动力特性系数决定了小扰动下转子－轴承系统的动特性。提出了一种基于CFD动网格技术的动特性系数计算新方法,首先利用轴颈在x、y两个方向的平移求出刚度系数,然后在此基础上利用简谐激振原理计算求得阻尼系数。

图4 动力特性系数求解流程图Fig.4 Schematic diagram of solution of dynamic coefficients

图4给出了利用新方法求解刚度和阻尼系数的过程。首先利用位移增量原理求出滑动轴承的刚度系数,使轴颈从平衡位置水平位移Δx(这时Δy= 0),计算时须在轴颈上加力ΔFx和ΔFy。这时可根据定义求得刚度系数kxx=ΔFx/Δx,kyx=ΔFy/Δx。同理通过使轴颈从平衡位置竖直位移Δy(这时Δx =0),可得到kxy=ΔFx/Δy,kyy=ΔFy/Δy。

同样以滑动轴承的静平衡位置为基准,分别在x、y方向施加频率为ω,相位差90°的已知简谐位移:

通过自编FLUENT-UDF程序积分求解得到这时作用在轴颈上的油膜力:

假定转子的质量忽略不计,则将式(5)、(6)代入牛顿第二定律可以很容易得到阻尼系数:

除了利用轴颈平移先求得刚度系数外,还可以对轴颈分别施加两次线性独立的简谐位移激励,从而直接求解出所需的4个刚度系数和阻尼系数。

3 数值计算结果

3.1 计算初值的影响

鉴于CFD求解方法能更好地考虑轴承内部3D瞬态流场的影响和表征流体流动形态,目前利用CFD技术求解滑动轴承动特性的方法主要有两种:一是利用相对坐标系的方法将瞬态流场近似为稳态流场求解,该方法仅局限于涡动轨迹为圆的圆柱轴承,若涡动轨迹或轴承结构稍作变化则无法使用;第二种较常用的方法是利用小扰动法直接瞬态求解,但该方法没有考虑计算初始值的影响,结果的准确性有待商榷。

非稳态计算属于微分方程的初值计算问题,不同的初值在计算初期得到的计算结果并不一样。图5给出了轴颈做指定的直线运动时油膜力的变化情况。从图5可以看出,数值计算初期,结果振荡比较大,随着计算的进行,数值结果逐渐趋于稳定。因此,传统的不考虑初值影响而只根据轴颈移动时的计算初始值求解得到的动特性系数是不准确的。

图5 位移扰动下油膜力的变化过程Fig.5 Oil film versus displacement step under displacement perturbation

3.2 计算结果

为了验证上述计算滑动轴承动力特性系数方法的适用性和准确性,对文献[11]中提供的典型滑动轴承模型进行计算,并与该文献提供的计算结果进行对比。其中VT-FAST,DyRoBes-BePerf和VT-EXPRESS都是目前应用比较广泛的基于经典Reynolds方程求解的转子动力学计算软件。

3.2.1 刚度系数

为了保证刚度系数求解的准确性,轴颈的平移距离要足够小,分别取位移步长(0.1、0.2 μm)使轴颈从静平衡位置沿x正方向和y负方向移动,得到的油膜力和拟合曲线如图6所示。从图6可知,随着轴颈在x和y方向的移动,所需施加在轴颈上的力与位移呈明显的线性比例关系,该比值即为刚度系数(N/m):

图6 刚度系数的拟合求解Fig.6 Fitting of stiffness coefficients

3.2.2 阻尼系数

不同于刚度系数的求解,阻尼系数是与轴颈扰动速度有关的量,因此本文中利用简谐激励原理进行求解。首先令轴颈以下列形式绕静平衡位置做圆轨迹涡动:

涡动轨迹如图7所示。

图7 轴颈圆形涡动轨迹Fig.7 Circle whirling orbit of journal

然后根据计算值拟合转子表面的油膜力,可得拟合公式为

将上述拟合值代入式(7)即可得到4个阻尼系数为

油膜力及其拟合值如图8所示。为了消除初始值的影响,通常需要将初期的计算结果舍弃。

3.2.3 结果对比

表1、2中的第4组数据分别为根据图6、8拟合计算出的刚度系数和阻尼系数。从表中可以看出,针对简单滑动轴承结构,本文计算的刚度和阻尼值与传统计算方法得到的值基本一致,虽然交叉系数的符号不一样,但经对比发现,符号与文献[12]中的经典数据完全一致,因此本文提出的基于CFD技术的滑动轴承动特性系数求解方法是可行的,且求解精度满足要求。该方法的优点是针对多油楔、带沟槽等复杂结构滑动轴承的动特性系数求解不需要重新推导复杂公式,并且可以很容易地考虑油膜本身质量的影响,尤其针对高速轻型转子轴承系统,惯性项的影响是至关重要的。

图8 x、y方向的油膜力变化及其拟合值Fig.8 Oil film force in x-axis and y-axis direction during circle whirling

表1 圆柱滑动轴承刚度系数Table 1 Stiffness coefficients of journal bearing

表2 圆柱滑动轴承阻尼系数Table 2 Damping coefficients of column journal bearing

4 结 论

(1)提出的新的变流域动网格方法理论上适合大部分结构的滑动轴承,实时地反映轴承润滑流场的动态变化。

(2)新方法的计算结果与已有的经典数据基本吻合,验证了该计算方法的可行性和有效性。

(3)CFD瞬态计算位移扰动下的润滑流场时,计算初值扰动比较大,会影响瞬态计算结果的准确性。

[1] TIEU A K,QIU Z L.Identification of sixteen dynamic coefficients of two journal bearings from experimental unbalance response[J].Wear,1994,177(1):63-69.

[2] ZHENG T,HASEBE N.Calculation of equilibrium position and dynamic coefficients of a journal bearing using free boundary theory[J].ASME Journal of Tribology, 2000,122:616-621.

[3] PERRINATO B,FLACK R D.Test results for a highly preloaded three-lobe journal bearing-effect of load orientation on static and dynamic characteristics[J].Journal of the Society of Tribology and Lubrication Engineers,2001, 57(9):23-30.

[4] WANG J K,KHONSARI M M.A new derivation for journal bearing stiffness and damping coefficients in polar coordinates[J].Journal of Sound and Vibration,2006, 290:500-507.

[5] BRAD A M,SAMUEL A H.Identifying bearing rotor-dynamic coefficients using an extended Kalman filter[J].Tribology Transactions,2009,52:671-679.

[6] GANDJALIKHAN N S A.Inertia effect on the thermohydrodynamic characteristics of journal bearings[J].Journal of Engineering Manufacture,2005,219(6):459-467.

[7] SHEIKH N Y,GANDJALIKHAN N S A.Three-dimensional numerical analysis of hydrodynamic characteristics of axial groove journal bearings running with ferrofluids under magnetic field[J].Journal of Engineering Tribology,2010,224(7):609-619.

[8] PANDAY K M,CHOUDHURY P L,KUMAR N P.Numerical unsteady analysis of thin film lubricated journal bearing[J].International Journal of Engineering and Technology,2012,4(2):185-191.

[9] DELIGANT M,PODEVIN P,DESCOMBES G.CFD model for turbocharger journal bearing performances[J].Applied Thermal Engineering,2011,31:811-819.

[10] SINGHAL A K,ATHAVALE M M,LI H Y,et al.Mathematical basis and validation of the full cavitation model[J].Journal of Fluids Engineering,2002,124 (3):617-624.

[11] GUO Z,HIRANO T,KIRK R G.Application of CFD analysis for rotating machinery:part 1-hydrodynamic,hydrostatic bearings and squeeze film damper[J].ASME Journal of Engineering for Gas Turbines and Tower, 2005,127:445-451.

[12] LUND J W.Rotor-bearing dynamics[R].Lecture Notes:Technical University of Denmark,1979.

(编辑 沈玉英)

Numerical calculation method of dynamic characteristics of journal bearing

LI Qiang1,XU Weiwei2,WANG Zhenbo1,JIN Youhai1,ZHENG Shuiying3
(1.College of Chemical Engineering in China University of Petroleum,Qingdao 266580,China; 2.College of Transport&Storage and Civil Engineering in China University of Petroleum,Qingdao 266580,China; 3.Institute of Chemical Machinery,Zhejiang University,Hangzhou 310027,China)

The computational fluid dynamics(CFD)and harmonic excitation method were applied to the numerical calculation of dynamic characteristics of journal bearing.By employing a new mesh movement approach based on structured grid,a new approach for calculating the dynamic characteristics of journal bearing was proposed based on the transient flow calculation.The stiffness and damping coefficients of a typical bearing were calculated by applying the new approach.The results obtained from the method were compared with previous classic computation results.The results show that the computation results of two methods are consistent.The effects of the computational initial value and the oil film fracture phenomenon are considered in this method,which is suitable for most of the journal bearing structures.The numerical method has good accuracy,and the method is valid.

journal bearing;computational fluid dynamic(CFD);dynamic coefficient;dynamic mesh;numerical calculation

TH 133

A

1673-5005(2014)05-0165-07

10.3969/j.issn.1673-5005.2014.05.024

2013－12－22

国家“863”高技术研究发展计划项目(2009AA04Z413);国家自然科学基金项目(51275452)

李强(1984－),男,讲师,博士,从事转子动力学、振动测试和故障诊断研究。E－mail:liqiangsydx@163.com。

李强,许伟伟,王振波,等.滑动轴承动力特性的数值计算方法[J].中国石油大学学报:自然科学版, 2014,38(5):165-171.

LI Qiang,XU Weiwei,WANG Zhenbo,et al.Numerical calculation method of dynamic characteristics of journal bearing [J].Journal of China University of Petroleum(Edition of Natural Science),2014,38(5):165-171.