基于动网格模型的液体动静压轴承刚度阻尼计算方法

本资料包含pdf文件1个，下载需要1积分

M ethod for Calculating Stiffness and Damping Coefi cients of HybridBearings Based on Dynamic M esh M odelXIONG Wanli HOU Zhiquan LO Lang YANG Xuebing YUAN Julong。

(1.National Research Center for High Eficiency Grinding，Hunan University，Changsha 410082；2.W ind Power Subsidiary Co.，Xiangtan Electric Manufacturing Group Co.Ltd.，Xiangtan 41 1 102；3.The Key Laboratory of M echanical Manufacture an d Automation of M inistry of Education，Zhejiang University ofTechnology，Hangzhou 310014)Abstract：For the reason也at Reynolds equation iS not suitable for reflecting nonlinear flow mechanics coupled by inertial effect ofoil circumferential flow，squeeze efect of the journal，and difusion efect of oil entry flow，the calculating accuracy of stifness anddamping coeficients can not sarisfy optimal design of high speed bearings.A new approach of calculating stifness and dampingcoeficients of joumal bearings is developed based on Navier-Stokes equation.The changing oil film force，the displacementperturbation，the velocity perturbation and the rotating speed are achieved by usiug user definition function programs.The movingboundary of the journal caused by rotational movement is changed into static boundary for avoiding the grid distortion.Thecoordination of moving grid caused by displacement perturbation an d velocity perturbation is updated by the spring smoothingtheoretical mode1.Comparing oil film force under steady condition with those of transient processes.the value scope of displacementperturbation and velocity perturbation are defined.The stiffness and damping coeficients of a typical bearing are calculated byapplying the approach combined with the perturbation theory，and the dyn amic development process of oil film thickness andthree-dimensiona1 pressure distribution are achieved.Finally，an experiment iS used to test the stifness coefficients of the bearingbased on the influence coeficient method，which are substantial in aeement with theoretical an alysis. The results show that thenumerical method is valid。

Key words：Hybridjournal bearing Stifness coeficients Damping coeficients Dynamic mesh国家自然科学基金(5O975O82，51275163)、教育部新世纪优秀人才支持计划fNCET-10-0368)1湖南省自然科学基金资助项目。

20120228收到初稿，20120828收到修改稿2012年 l2月 熊万里等：基于动网格模型的液体动静压轴承刚度阻尼计算方法 l190 前言液体动静压轴承以其高精度、高刚度、高阻尼吸振性和长寿命等优势，在超高速精密磨削领域获得广泛应用。目前，国家高效磨削工程技术研究中心研制出用于磨削外圆面工件的超高速液体动静压轴承电主轴，额定功率达35 kW，主轴最高转速达9 l1in；用于磨削凸凹面工件的液体动静压轴承电主轴，额定功率达 12kW，最高转速达30kr/minOJ o随着电主轴转速的进-步提高，电主轴系统因高速失稳出现刮瓦、抱轴等现象，导致电主轴系统失效的可能性逐步增加。高速稳定性已成为当前超高速精密电主轴系统动力学分析的重要课题，其难点就在于支承轴承动态刚度和阻尼的精确计算。

动静压轴承的动态刚度阻尼反映了轴颈受外载荷作用下振动位移和振动速度与油膜力之间的内在联系。1965年，LUNDt21采用4个刚度系数和4个阻尼系数建立了挠性 轴承.转子系统”模型，从理论上解释了系统涡动失稳的形成原因，揭示了 轴承.转子系统”涡动失稳的物理本质。8个系数被广泛用于表征动静压轴承的动态特性和分析轴承的高速稳定性。

计算轴承的刚度阻尼，常用的方法有差分法和偏导数法，其关键科学问题在于求解动态流体润滑方程得到油膜的压力分布。通常情况下，动态雷诺方程作为润滑方程用于计算动静压轴承的刚度阻尼p J。由于雷诺方程难以精确反映转速引起的油流周向惯性效应、轴颈环周向动态挤压效应和静压效应之间非线性耦合关系及其对三维物理场(速度尝温度惩压力场)的影响，油膜刚度阻尼的计算精度受到制约，因此有必要直接以Navier-Stokes方程为基础，精确研究轴承结构、主轴转速、位移扰动和速度扰动等对油膜刚度和阻尼的影响规律。文献[8-10采用 基 于 N-S 方程 的计 算 流 体动 力学(Computational fluid dynamics，CFD)开展了静压轴承承载特性的研究工作，证实了N-S方程在表征复杂求解域流体流动形态方面可弥补雷诺方程的不足。文献[11]进-步引入了 可动边界”的概念，采用动网格方法计算了动压轴承在外载荷作用下轴承的压力分布，但尚未解决轴颈旋转带来的网格扭曲问题。

鉴于 CFD在表征流体流动形态方面的优势及其计算轴承动态性能方面的不足，本文提出了基于CFD 动网格模型的液体动静压轴承油膜刚度阻尼系数的计算方法，该方法在求解瞬态N-s方程基础上，采用自定义程序结合动网格模型，动态反映了轴颈转速、位移扰动和速度扰动等因素对油膜压力场的影响，从而使油膜力计算结果与高速时实际工况更吻合；通过对比分析稳态计算和瞬态计算油膜力的差异，界定了采用动网格方法计算刚度阻尼时扰动位移和扰动速度的合理取值范围。通过对比典型算例和试验测试结果，验证了动网格计算方法的有效性和可行性。

1 刚度阻尼动网格计算方法基于动网格模型的轴承刚度阻尼计算方法是在CFD软件基础上，采用自定义程序实现轴颈扰动功能，结合差分计算模型，计算轴承油膜刚度阻尼的- 种新方法。该方法利用 CFD软件求解 N.S方程得到轴颈位置变动前后流体的瞬态压力场，通过自定义程序遍历轴颈表面压力得到不同方向的瞬态油膜力； 轴颈位置变动过程中出现的油膜网格变形由动网格模型更新；轴颈位置变动的三个运动参量轴颈转速、扰动位移和扰动速度由自定义程序指定，同时有效解决因轴颈旋转引起的油膜网格扭曲问题；通过对比扰动位移和扰动速度对油膜力的影响，界定计算刚度和阻尼的合理取值范围。

1.1 刚度阻尼自勺差分计算模型轴颈在平衡位置受外载荷作用，对轴承油膜产生位移扰动和速度扰动，油膜力将发生变化。在小扰动条件下，轴承油膜刚度阻尼可采用差分法求解。

根据摄动理论，将油膜力线性化处理，油膜力的变化与扰动位移和扰动速度之间的关系L4 如式(1)所示㈤(乏k A]fA y；] )㈣式中， 和 为油膜力的变化； 和 七 是直接刚度； 和 ‰ 为交叉刚度； 和c 为直接阻尼： 和c 为交叉阻尼；Ax和 为扰动位移；X和Y 为扰动速度。

采用差分模型计算轴承的刚度和阻尼系数，需要考虑位移和速度对不同方向油膜力的影响。刚度和阻尼的计算如式(2)所示-AFdxx-AFdyx： Cyx：Z"vyxJ Ax Ax譬 等 c (2)式中， ，为位移扰动引起油膜力的变化； 为速度扰动引起油膜力的变化；f、 表示X、Y中的某l20 机 械 工 程 学 报 第 48卷第 23期- 个。

由式(2)可知，采用动网格方法计算油膜刚度阻尼系数，需要解决三个问题：① 计算刚度阻尼需要计算轴颈在位移和速度变动前后瞬态油膜力，该油膜力可通过求解瞬态压力-速度耦合方程获得；② 轴颈受外载荷和自转影响，轴颈在轴承油膜中的相对位置将产生变化，油膜网络挤压变形并由动网格模型更新；③ 计算刚度阻尼需要在轴颈动边界条件中引入轴颈旋转、速度扰动和位移扰动，轴颈运动控制由白定义程序实现。

1.2 瞬态油膜力计算在 CFD软件中，油膜流动的基本动力润滑理论是N-S方程。数值求解润滑方程的基本思想是有限体积法，将流体区域离散为流体微元，计算出微元的物理参量。其核心思想是子域法加离散法。

在稳定工况下，动静压轴承间隙内流体为黏性层流，不考虑流体的粘温效应，动力润滑方程包括连续性方程和动量方程，其瞬态通用表达形式J J如式(3)所示( )div(puq)div(Fgradcp) (3)式中，P表示流体密度； 为流场速度矢量； 为通用变量；厂表示广义扩散系数； 为广义源项。

式(3)各项分别表征流体的瞬态项、对流项、扩散项和源项。对比广义雷诺方程的假设条件及其所包含的压力项、楔入项和挤压项，各项之间呈线性叠加关系。显然，CFD控制方程更适于表征间隙内流体的三维环周向流动形态。

动静压轴承间隙内流体的流场形态主要由供油压力产生的压力流和轴颈旋转引起的压力流和速度流非线性耦合形成，表现为静压与动压的强耦合效应。流体的压力-速度耦合方程如式(4)所示a n( )div(pqm)div(Fgrad)- (4)az d，7式中，r/为容积表面外法线方向。

基于有限体积法，将式(4)改写，如式(5)所示J 善( )d J diV( )d ：J diV厂 ad )d J d (5)式中，V为润滑油的体积。

采用高斯散度公式，将体积分转换为面积分，得到瞬态问题的表达式如式(6)所示(Lm，dV)dtLL，7.( ) d J△fJ ，7·(俘ad ) dH Jf sov (6)式中， 为润滑油的面积，△f为时问变化量。式(6)是 CFD软件 Fluent进行瞬态计算的理论基矗Fluent软件求解压力-速度耦合方程的基本算法是SIMPLE(Semi-implicit method for pressure linkedequations)算法。SIMPLE算法需要预定义流场压力。

初始压力和参考压力的设置对计算结果影响比较大。尤其是动网格瞬态计算过程中出现回流现象时，预先建立流场压力分布对计算的收敛性至关重要。

对比工程实践，将进口边界条件设置为供油压力，参考压力为 1 MPa，出油压力为0 MPa，即泄油端与大气连通。

当采用动网格模型处理非定常压力-速度耦合流动 问题时 ，-般采用 PISO (Pressure-implicitsplicing of operators)算法，其计算基础是 SIMPLE算法。两者的区别在于预设压力求解动量方程和连续方程后，PISO算法需要连续两次修正压力和速度方程，使压力和速度分别达到该时刻的两阶和三阶的精度。

鉴于轴颈高速对流场的影响，压力方程的插值离散格式选择 PRESTO，动量方程的插值离散格式为 QUICK。QUICK为流动的二次迎风插值格式，特别适于与流动方向对齐的结构网格，比二阶迎风格式具有更精确的计算结果。其广义的表达式如(7)所示钆 纯) 钆 纯j(s 2Sc sr 、 l面 - J (7)式中，下标P为控制体积，下标E和 分别表示东西方向与P相邻的控制体积； ， ， 分别为W，P，E相对应的控制体积的边长。伽 ，仇 ，分别为P，E，W控制体积的通量， 为P控制体积西面的通量。

求解压力-速度耦合方程，得到油膜的三维物理常通过 自定义程序循环指令遍历轴颈面域各节点的压力和黏性力并积分，得到不同方向的油膜力。

当轴颈平衡位置变动过大且转速过高时，导致油膜破裂出现负压力，即油膜空化现象。计算油膜压力时，常用的处理方法是将负值区域压力置零。该处理方法和油膜力的计算所涉及的用户定义函数主要包括宏指令 DEFINE ON DEMAND、循环函数begin f loop(， 和面积函数FAREA(A， 制。

1.3 弹簧光顺动网格模型动网格模型主要应用于解决区域几何形状随时间变化的流体动力学问题。计算动静压轴承刚度阻尼系数时，网格变形主要考虑轴承间隙内油膜随轴颈作刚体运动而产生的变化。在速度扰动和位移扰2012年 l2月 熊万里等：基于动网格模型的液体动静压轴承刚度阻尼计算方法对前面建立的计算模型，调用自定义程序使轴颈朝 方向运动，计算轴颈在不同扰动速度和扰动位移下 和 Y方向的油膜力，并与稳态计算结果比较(扰动速度和扰动位移均为零)。 鑫每个迭代时间步长，轴颈扰动位移为 1 gm，轴 岸颈在不同扰动速度 v作用下，X和 Y两个方向分力与稳态计算结果对比如图4所示。

露1御鹫寝偏心率 偏心率(a) 方向油膜力的变化 (b) 方向油膜力的变化图4 扰动速度对X和Y方向油膜承载力的影响由图4a可知，扰动速度对 X方向分力的影响不大。与稳态计算结果较靠近的速度为 2.0 gm/s 和20.0rn/s 。扰动速度为 20rn/s 时，最大差距12.5%，平均差距 7.7%；扰动速度为2.0rn/s 时，平均差距 7.5%。速度相差 10倍的情况下，平均差距相差不大。

扰动速度对Y方向分力的影响比较明显。由图4b可知，扰动速度为 0.2 mm/s 时，最大差距为24.5%，平均差距 16%。扰动速度为20.0 rn/s 时，最大差距2.5%，平均差距 1.8%；速度为2.0 gm/s 时，最大差距 0.2%。

由上可知，在同样的扰动速度和扰动位移条件下，采用动网格方法计算的X方向和Y方向分力与稳态计算结果的差距不相同，挤压运动方向的差距大于与其垂直方向的差距，这与压力-速度耦合方程的求解模型有关。稳态计算和瞬态计算的基本润滑方程不同(与时间关联)，求解模型不同，计算结果必然存在偏差。但从总体上看，在单个迭代步长内，扰动位移为 1.0 m 时，扰动速度不大于 20rn/s的计算结果与稳态计算基本吻合。这个速度可视为刚度和阻尼计算的分界点。

为减低扰动速度对油膜力的影响，单个迭代步长扰动速度取0.2 gm/s ，在不同扰动位移S作用下，X和Y方向分力与稳态计算对比如图5所示。

对于X方向分力的影响，扰动位移为2.0 m时，最大差距为 25.5%，平均差距为 19.5%；扰动位移为 0.5 um时，最大差距 13.1%，平均差距 7.8%，而位移为 1.0 gm时平均差距为5.7%。

偏心率 s(a) 方向油膜力的变化偏心率(b)y方向油膜力的变化图 5 位移扰动对X和Y方向承载力的影响对于Y方向分力的影响，位移扰动为2.0 gm时，最大差距为 1 1%，平均差距为 6.4%；扰动位移为0.5 m时，最大差距 1.3%，平均差距 0.6%，而位移为 1.0 gm时最高差距为0.3%。

对于扰动位移对 X和Y方向油膜力的影响与扰动速度不同，Y方向的差距小于 X方向的差距。由图 5可知，位移为 1.0 gm与 0.5 m时，油膜力较为接近。X方向的差距为 2.1%，Y方向的差距为 0.3%。

由上可知，扰动速度 v≤2O grn/s ，扰动位移≤1 m，扰动速度和扰动位移对油膜力的影响较小，平均差距不大于 8%。动网格的瞬态计算与稳态计算基本相 当，可视为刚度计算。扰动速度v>0.2 mm/s 时，速度对油膜力的影响比较明显。动网格的瞬态计算不能近似为稳态计算，可视为阻尼计算。在合适的速度和位移取值范围内，采用动网格计算轴承刚度和阻尼系数是切实可行的。

2 典型轴承的动网格计算采用动网格方法计算轴承的刚度阻尼系数包括两个部分：-是计算主轴在轴承间隙内的平衡位置；二是在平衡位置处计算轴承的刚度阻尼。

轴承结构如图 l所示，轴承设计参数和润滑油如前述。主轴的质量为 5.54 〃模时取轴承与轴颈的同心位置为初始位置。重力作用下的最大移动速度为 2.0 s ，最大位移为 1.0 m。计算主轴平衡位置时理论上应满足油膜力水平分力为 0，垂直分力等于55.4 N〖虑到计算效率，误差为0.001 N时停止迭代。主轴的偏位角和偏心距随转速的变化如图6所示。

计算刚度阻尼时需偏离平衡位置-定距离以适应轴颈的运动。采用Gambit建模时需要将轴颈往或Y方向偏移 1.0 gm。计算轴承刚度时，扰动速度为2.0rn/s ，扰动位移 1.0 m；计算阻尼时，扰动速度0.5 mm/s ，扰动位移1.0 m，分别计算出轴承2012年 12月 熊万里等：基于动网格模型的液体动静压轴承刚度阻尼计算方法个速度点空运行半小时后，在 和 方向先后施加20 N和 40 N作用力，每个作用力加载 3次，记录每次加载力时 和Y方向的位移偏移量，取平均值作为该速度点和作用力的有效位移。试验测试出48组数据，将结果代入式(13)，计算出轴承的刚度系数。

3.3 理论与试验对比分析将动网格数值计算结果(图 7)与试验测试结果进行比对，刚度随转速的变化趋势如图 11所示。

言 1Z蠢憋I删言Z塞甚转速n/(kr/min)(a)直接刚度O- .-F0N- 2ON]A F 40 N- ： rt-言- 星权 -转速川kr/min)(b)直接刚度 kyy转速n/(kmin) 转速n/(kr/min)(c)交叉刚度 (d)交叉刚度图l1 动网格计算结果和试验测试结果对比通过对比动网格数值计算 F0 N与试验测试F20 N和F40 N平均值可知，两者偏差不大。直接刚度 最大偏差为 10.4%，平均偏差 6.16%；直接刚度 ‰ 的最大偏差 20.9%，平均偏差 13.3%；交叉刚度 最大偏差 15%，平均偏差 9.8%；交叉刚度 ‰ 与 相差不大。

对于局部各数据点，动网格计算与试验测试结果并非完全-致，这与各 自的外部因素有关。动网格计算偏向理论分析，其结果与求解器模型、网格质量、边界条件和动网格模型以及参数设置有关；试验测试侧重于工程实践，其结果受转子质量动不平衡、轴承加工误差、轴系装配误差、转速的脉动和供油的平稳性以及加载力的平稳性和方向的平行度等因素影响，以至于F20 N和 F40 N两组数据也存在-定差距。对于总体变化趋势，动网格计算结果与试验测试数据基本吻合，平均误差不超过12.5%，体现了典型轴承动力特性系数的内在-致性。

4 结论(1)采用自定义程序实现轴颈位移扰动、速度扰动、旋转和油膜力计算等功能，并利用弹簧光顺模型更新网格的挤压变形，计算了轴承的刚度和阻尼系数。通过试验测试轴承的刚度系数，发现两者基本吻合，表明采用动网格模型计算动静压轴承刚度阻尼的数值方法是可行和有效的，从而为动静压轴承动力特性系数计算提供了-种新方法。

(2)将轴颈的旋转动边界条件转换为静边界条件，并将转速加载在静边界上，有效避免了因轴颈旋转引起的油膜网格扭曲畸变，为旋转机械流体的动边界条件加载提供了-种新途径。

(3)对比研究了扰动速度和扰动位移对油膜力的影响规律，界定了计算刚度阻尼系数时扰动速度和扰动位移的合理取值范围，解决了动网格技术计算轴承刚度阻尼系数的关键问题。

(4)再现了轴颈受位移扰动和速度扰动情况下油膜厚度的动态发展过程，反映了三维流场压力分布的时变特性，使计算过程更符合工程实际。