第４８卷第２３期２０１２年１２月
机械工程学报
ＪＯＵＲＮＡＬＯＦ
Ｖｂｌ．４８Ｄｅｃ．
ＮＯ．２３２０１２
ＭＥＣＨＡＮＩＣＡＬ
ＥＮＧＩＮＥＥＲＩＮＧ
ＤｏＩ：１０．３９０１／ＪＭ［Ｅ．２０１２．２３．１１８
基于动网格模型的液体动静压轴承
刚度阻尼计算方法木
熊万里１
侯志泉１
吕
浪１
阳雪兵２袁巨龙３
（１．湖南大学国家高效磨削工程技术研究中心长沙４１００８２；
２．湘潭电机股份有限公司风能分公司湘潭４１１１０２；
３．浙江工业大学机械制造及自动化教育部重点实验室杭州３１００１４）
摘要：针对雷诺方程不能反映高速油流周向惯性效应、轴颈周向动压效应和静压扩散效应之间的非线性耦合关系及其对三维流场特性的影响，进而导致动静压轴承刚度阻尼的计算精度难以满足高转速精密轴承设计需要之现状，提出基于纳维一斯托克斯方程的动网格计算轴承刚度阻尼新方法。该方法采用自定义程序实现轴颈的旋转速度、位移扰动和速度扰动以及扰动过程中油膜力的计算等功能，成功将轴颈的旋转动边界转换为静边界，有效避免因轴颈旋转引起的网格畸变，并利用弹簧光顺模型更新轴颈受速度扰动和位移扰动引起的油膜网格挤压变形。通过对比研究瞬态和稳态条件下油膜力的变化，界定计算刚度阻尼系数时速度扰动和位移扰动的合理取值范围。利用所提出的方法计算典型轴承的刚度阻尼系数，再现动静压轴承油膜厚度和三维流场压力分布的动态发展过程，并采用影响系数法对该轴承的刚度系数进行试验测试。通过对比试验数据和动网格计算结果发现两者基本吻合，表明所提出的计算方法是有效可行的。关键词：动静压轴承刚度系数阻尼系数动网格中图分类号：ＴＨｌ３３
Ｍｅｔｈｏｄ
ｆｏｒＣａｌｃｕｌａｔｉｎｇＳｔｉｆｆｎｅｓｓａｎｄＤａｍｐｉｎｇ
ＢｅａｒｉｎｇｓＢａｓｅｄ
ｏｎ
Ｃｏｅｆｆｉｃｉｅｎｔｓ
ＤｙｎａｍｉｃＭｅｓｈＭｏｄｅｌ
Ｌａｎ９１
ＹＡＮＧＸｕｅｂｉｎ９２
ｏｆＨｙｂｒｉｄ
ＸＩＯＮＧ
Ｗａｎｌｉｌ
ＨＯＵＺｈｉｑｕａｎｌ
ＬＯ
ＹＵＡＮＪｕｌｏｎ９３
（１．ＮａｔｉｏｎａｌＲｅｓｅａｒｃｈＣｅｎｔｅｒｆｏｒＨｉｇｈＥｆｆｉｃｉｅｎｃｙＧｒｉｎｄｉｎｇ，ＨｕｎａｎＵｎｉｖｅｒｓｉｔｙ，Ｃｈａｎｇｓｈａ４１００８２；２．ＷｉｎｄＰｏｗｅｒＳｕｂｓｉｄｉａｒｙＣｏ．，ＸｉａｎｇｔａｎＥｌｅｃｔｒｉｃＭａｎｕｆａｃｔｕｒｉｎｇＧｒｏｕｐＣｏ．Ｌｔｄ．，Ｘｉａｎｇｔａｎ４１１１０２；
３．ＴｈｅＫｅｙＬａｂｏｒａｔｏｒｙ
ｏｆＭｅｃｈａｎｉｃａｌＭａｎｕｆａｃｔｕｒｅａｎｄＡｕｔｏｍａｔｉｏｎｏｆＭｉｎｉｓｔｒｙ
ｏｆＥｄｕｃａｔｉｏｎ，
ＺｈｅｊｉａｎｇＵｎｉｖｅｒｓｉｔｙｏｆＴｅｃｈｎｏｌｏｇｙ，Ｈａｎｇｚｈｏｕ３１００１４）
Ａｂｓｔｒａｃｔ：Ｆｏｒｔｈｅ
ｒｅａｓｏｎ
ｍａｔＲｅｙｎｏｌｄｓｅｑｕａｔｉｏｎｉｓｎｏｔｓｕｉｔａｂｌｅｆｏｒｒｅｆｌｅｃｔｉｎｇｎｏｎｌｉｎｅａｒｆｌｏｗｍｅｃｈａｎｉｃｓｃｏｕｐｌｅｄｂｙｉｎｅｒｔｉａｌｅｆｆｅｃｔｏｆ
ｏｉｌｃｉｒｃｕｍｆｅｒｅｎｔｉａｌｆｌｏｗ，ｓｑｕｅｅｚｅｅｆｆｅｃｔｏｆｔｈｅｄａｍｐｉｎｇｃｏｅｆｆｉｃｉｅｎｔｓ
ｃａｌｌ
ｊｏｕｒｎａｌ．ａｎｄ
ｏｎ
ｄｉｆｆｕｓｉｏｎｅｆｆｅｃｔｏｆｏｉｌｅｎｔｒｙｆｌｏｗ，ｔｈｅｃａｌｃｕｌａｔｉｎｇａｃｃｕｒａｃｙｏｆｓｔｉｆｆｎｅｓｓａｎｄ
ｎｏｔ
ｓａｔｉｓｆｙｏｐｔｉｍａｌｄｅｓｉｇｎｏｆｈｉｇｈｓｐｅｅｄｂｅａｒｉｎｇｓ．Ａｎｅｗａｐｐｒｏａｃｈｏｆｃａｌｃｕｌａｔｉｎｇ
ｓｔｉｆｆｎｅｓｓ
ａｎｄｄａｍｐｉｎｇ
ｃｏｅｆｆｉｃｉｅｎｔｓｏｆｊｏｕｒｎａｌｂｅａｒｉｎｇｓｉｓｄｅｖｅｌｏｐｅｄｂａｓｅｄ
ｐｅｒｔｕｒｂａｔｉｏｎ，ｔｈｅｖｅｌｏｃｉｔｙｐｅｒｔｕｒｂａｔｉｏｎｂｏｕｎｄａｒｙ
ｏｆｔｈｅ
Ｎａｖｉｅｒ－Ｓｔｏｋｅｓｅｑｕａｔｉｏｎ．Ｔｈｅｃｈａｎｇｉｎｇｏｉｌｆｉｌｍｆｏｒｃｅ，ｔｈｅｄｉｓｐｌａｃｅｍｅｎｔ
ａｒｅ
ａｎｄ
ｔｈｅｒｏｔａｔｉｎｇｓｐｅｅｄａｃｈｉｅｖｅｄｂｙｕｓｉｎｇ
ｕｓｅｒ
ｄｅｆｉｎｉｔｉｏｎｆｕｎｃｔｉｏｎ
ｐｒｏｇｒａｍｓ．Ｔｈｅ
ｍｏｖｉｎｇ
ｊｏｕｒｎａｌｃａｕｓｅｄｂｙｒｏｔａｔｉｏｎａｌｍｏｖｅｍｅｎｔｉｓ
ｃｏｏｒｄｉｎａｔｉｏｎｏｆｍｏｖｉｎｇｇｒｉｄｃａｕｓｅｄｂｙｄｉｓｐｌａｃｅｍｅｎｔｔｈｅｏｒｅｔｉｃａｌｍｏｄｅｌ．Ｃｏｍｐａｒｉｎｇｏｉｌｆｉｌｍｆｏｒｃｅｐｅｒｔｕｒｂａｔｉｏｎ
ｃｈａｎｇｅｄｉｎｔｏｓｔａｔｉｃｂｏｕｎｄａｒｙｆｏｒ
ｐｅｒｔｕｒｂａｔｉｏｎａｎｄｖｅｌｏｃｉｔｙｐｅｒｔｕｒｂａｔｉｏｎｉｓ
ａｖｏｉｄｉｎｇｔｈｅｇｒｉｄｄｉｓｔｏｒｔｉｏｎ．Ｔｈｅｕｐｄａｔｅｄｂｙｔｈｅｓｐｒｉｎｇｓｍｏｏｔｈｉｎｇ
ｕｎｄｅｒｓｔｅａｄｙ
ｃｏｎｄｉｔｉｏｎｗｉｔｈｔｈｏｓｅｏｆｔｒａｎｓｉｅｎｔｐｒｏｃｅｓｓｅｓ，ｔｈｅｖａｌｕｅｓｃｏｐｅｏｆｄｉｓｐｌａｃｅｍｅｎｔ
ａｎｄｖｅｌｏｃｉｔｙ
ｐｅｒｔｕｒｂａｔｉｏｎａｒｅｄｅｆｉｎｅｄ．Ｔｈｅ
ｓｔｉｆｆｎｅｓｓ
ａｎｄｄａｍｐｉｎｇｃｏｅｆｆｉｃｉｅｎｔｓｏｆ
ａ
ｔｙｐｉｃａｌｂｅａｒｉｎｇａｒｅｃａｌｃｕｌａｔｅｄｂｙｐｒｏｃｅｓｓｏｆｏｉｌｆｉｌｍｔｈｉｃｋｎｅｓｓａｎｄ
ａｐｐｌｙｉｎｇｔｈｅａｐｐｒｏａｃｈｃｏｍｂｉｎｅｄｗｉｔｈｔｈｅｐｅｒｔｕｒｂａｔｉｏｎｔｈｅｏｒｙ，ａｎｄｔｈｅｄｙｎａｍｉｃｄｅｖｅｌｏｐｍｅｎｔ
ｔｈｒｅｅ－ｄｉｍｅｎｓｉｏｎａｌｐｒｅｓｓｕｒｅｄｉｓｔｒｉｂｕｔｉｏｎａｒｅａｃｈｉｅｖｅｄ．Ｆｉｎａｌｌｙ，ａｌｌｅｘｐｅｒｉｍｅｎｔｉｓｕｓｅｄｔｏｔｅｓｔｔｈｅｓｔｉｆｆｎｅｓｓｃｏｅｆｆｉｃｉｅｎｔｓｏｆｔｈｅｂｅａｒｉｎｇｂａｓｅｄ
ｏｎ
ｔｈｅｉｎｆｌｕｅｎｃｅｃｏｅｆｆｉｃｉｅｎｔｍｅｔｈｏｄ，ｗｈｉｃｈａｒｅｓｕｂｓｔａｎｔｉａｌｉｎａｇｒｅｅｍｅｎｔｗｉｔｈｔｈｅｏｒｅｔｉｃａｌａｎａｌｙｓｉｓ．ＴｈｅｒｅｓｕＲｓ
ｓｈｏｗｔｈａｔｔｈｅ
ｎｕｍｅｒｉｃａｌｍｅｔｈｏｄｉｓｖａｌｉｄ．
Ｋｅｙｗｏｒｄｓ：Ｈｙｂｒｉｄｊｏｕｒｎａｌｂｅａｒｉｎｇ
Ｓｔｉｆｆｎｅｓｓ
ｃｏｅｆｆｉｃｉｅｎｔｓＤａｍｐｉｎｇｃｏｅｆｆｉｃｉｅｎｔｓＤｙｎａｍｉｃｍｅｓｈ
・国家自然科学基金（５０９７５０８２，５１２７５１６３）、教育部新世纪优秀人才支持计划（ＮＣＥＴ－１０—０３６８）和湖南省自然科学基金资助项目。
２０１２０２２８收到初稿，２０１２０８２８收到修改稿
万方数据
２０１２年１２月熊万里等：基于动网格模型的液体动静压轴承刚度阻尼计算方法
１１９
０前言
液体动静压轴承以其高精度、高刚度、高阻尼吸振性和长寿命等优势，在超高速精密磨削领域获得广泛应用。目前，国家高效磨削工程技术研究中
心研制出用于磨削外圆面工件的超高速液体动静压
轴承电主轴，额定功率达３５ｋＷ，主轴最高转速达
９
ｋｒ／ｍｉｎ；用于磨削凸凹面工件的液体动静压轴承电
主轴，额定功率达１２ｋＷ，最高转速达３０
ｋｒ／ｍｉｎ【ｌＪ。
随着电主轴转速的进一步提高，电主轴系统因高速
失稳出现刮瓦、抱轴等现象，导致电主轴系统失效的可能性逐步增加。高速稳定性已成为当前超高速
精密电主轴系统动力学分析的重要课题，其难点就
在于支承轴承动态刚度和阻尼的精确计算。
动静压轴承的动态刚度阻尼反映了轴颈受外载
荷作用下振动位移和振动速度与油膜力之间的内在
联系。１９６５年，ＬＵＮＤＬ２ｊ采用４个刚度系数和４个
阻尼系数建立了挠性“轴承．转子系统”模型，从理论上解释了系统涡动失稳的形成原因，揭示了“轴承．转子系统”涡动失稳的物理本质。８个系数被广泛用于表征动静压轴承的动态特性和分析轴承的高
速稳定性。
计算轴承的刚度阻尼，常用的方法有差分法和偏导数法，其关键科学问题在于求解动态流体润滑方程得到油膜的压力分布。通常情况下，动态雷诺
方程作为润滑方程用于计算动静压轴承的刚度阻尼【ｊ√Ｊ。由于雷诺方程难以精确反映转速引起的油流周向惯性效应、轴颈环周向动态挤压效应和静压效
应之间非线性耦合关系及其对三维物理场（速度场、
温度场和压力场）的影响，油膜刚度阻尼的计算精度
受到制约，因此有必要直接以Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程为
基础，精确研究轴承结构、主轴转速、位移扰动和速度扰动等对油膜刚度和阻尼的影响规律。文献
【８—１０］采用基于Ｎ—Ｓ方程的计算流体动力学
（Ｃｏｍｐｕｔａｔｉｏｎａｌ
ｆｌｕｉｄ
ｄｙｎａｍｉｃｓ，ＣＦＤ）开展了静压轴承承载特性的研究工作，证实了Ｎ—ｓ方程在表征复
杂求解域流体流动形态方面可弥补雷诺方程的不
足。文献【１１］进一步引入了“可动边界”的概念，
采用动网格方法计算了动压轴承在外载荷作用下轴
承的压力分布，但尚未解决轴颈旋转带来的网格扭
曲问题。
鉴于ＣＦＤ在表征流体流动形态方面的优势及其计算轴承动态性能方面的不足，本文提出了基于
ＣＦＤ动网格模型的液体动静压轴承油膜刚度阻尼万方数据
系数的计算方法，该方法在求解瞬态Ｎ—Ｓ方程基础
上，采用自定义程序结合动网格模型，动态反映了轴颈转速、位移扰动和速度扰动等因素对油膜压力
场的影响，从而使油膜力计算结果与高速时实际工
况更吻合；通过对比分析稳态计算和瞬态计算油膜力的差异，界定了采用动网格方法计算刚度阻尼时扰动位移和扰动速度的合理取值范围。通过对比典型算例和试验测试结果，验证了动网格计算方法的
有效性和可行性。１
刚度阻尼动网格计算方法
基于动网格模型的轴承刚度阻尼计算方法是在
ＣＦＤ软件基础上，采用自定义程序实现轴颈扰动功能，结合差分计算模型，计算轴承油膜刚度阻尼的一种新方法。该方法利用ＣＦＤ软件求解Ｎ．Ｓ方程得到轴颈位置变动前后流体的瞬态压力场，通过自定义程序遍历轴颈表面压力得到不同方向的瞬态油膜
力；轴颈位置变动过程中出现的油膜网格变形由动
网格模型更新；轴颈位置变动的三个运动参量轴颈
转速、扰动位移和扰动速度由自定义程序指定，同时有效解决因轴颈旋转引起的油膜网格扭曲问题；
通过对比扰动位移和扰动速度对油膜力的影响，界
定计算刚度和阻尼的合理取值范围。
１．１
刚度阻尼的差分计算模型
轴颈在平衡位置受外载荷作用，对轴承油膜产
生位移扰动和速度扰动，油膜力将发生变化。在小
扰动条件下，轴承油膜刚度阻尼可采用差分法求解。根据摄动理论，将油膜力线性化处理，油膜力的变化与扰动位移和扰动速度之间的关系１４１如式（１）所示
式中，峨和峨为油膜力的变化；ｋ和Ｊｉ｝珂是直
∞］＝（乏ｋ％．Ａ㈣ａｙｊ＋（墓北）㈣接刚度；ｋ和ｋ为交叉刚度；Ｃｘｘ和．ｃ抑为直接阻尼；ｃ。和ｃ。为交叉阻尼；缸和每为扰动位移；ｘ７
和Ｙ７为扰动速度。
采用差分模型计算轴承的刚度和阻尼系数，需要考虑位移和速度对不同方向油膜力的影响。刚度
和阻尼的计算如式（２）所示
胁ｌｌ
％２言，蜴肛
印＝
％２亍刖‰
胁
＝
警堕妙％２亏
，甜’劬
细
＝
等堕少
呦２——２
心‰
ｙ。
（２）
式中，刖％为位移扰动引起油膜力的变化；刖％为速度扰动引起油膜力的变化；ｆ、，表示ｘ、Ｙ中的某
１２０
机械工程学报第４８卷第２３期
一个。
由式（２）可知，采用动网格方法计算油膜刚度阻
尼系数，需要解决三个问题：①计算刚度阻尼需要
计算轴颈在位移和速度变动前后瞬态油膜力，该油膜力可通过求解瞬态压力一速度耦合方程获得；
②轴颈受外载荷和自转影响，轴颈在轴承油膜中的
相对位置将产生变化，油膜网络挤压变形并由动网
格模型更新；③计算刚度阻尼需要在轴颈动边界条
件中引入轴颈旋转、速度扰动和位移扰动，轴颈运动控制由自定义程序实现。１．２瞬态油膜力计算
在ＣＦＤ软件中，油膜流动的基本动力润滑理论是Ｎ—Ｓ方程。数值求解润滑方程的基本思想是有限体积法，将流体区域离散为流体微元，计算出微元
的物理参量。其核心思想是子域法加离散法。
在稳定工况下，动静压轴承间隙内流体为黏性
层流，不考虑流体的粘温效应，动力润滑方程包括
连续性方程和动量方程，其瞬态通用表达形式【ｌ压¨Ｊ
如式（３）所示
ａ
鼍（ｐ缈）＋ｄｉＶ（ｐ“伊）＝ｄｉｖ（Ｆｇｒａｄ口ｏ）＋－品
（３）
式中，ｐ表示流体密度；Ｈ为流场速度矢量；缈为
通用变量；厂表示广义扩散系数；＆为广义源项。
式（３）各项分别表征流体的瞬态项、对流项、扩
散项和源项。对比广义雷诺方程的假设条件及其所包含的压力项、楔入项和挤压项，各项之间呈线性叠加关系。显然，ＣＦＤ控制方程更适于表征间隙内流体的三维环周向流动形态。
动静压轴承间隙内流体的流场形态主要由供油
压力产生的压力流和轴颈旋转引起的压力流和速度
流非线性耦合形成，表现为静压与动压的强耦合效
应。流体的压力．速度耦合方程如式（４）所示
÷（胛）＋ｄｉｖ（ｐｆｏｕ）＝ｄｉｖ（ｒｇｒａｄ伊）－署＋昂（４）
ａａｎ
ａｆ
ｏｒ／
式中，ｒ／为容积表面外法线方向。
基于有限体积法，将式（４）改写，如式（５）所示
ａ
．
ｆｒＤ（ＰｑｏｄＶ＋Ｊ矿ｄｉＶ（ｐ缈甜）ｄ矿＝
ｊｙｄｉｖ（Ｆｇｍｄ９）ｄＶ＋Ｊｒ＂Ｓ伊ｄＶ
（５）
式中，ｙ为润滑油的体积。
采用高斯散度公式，将体积分转换为面积分，
得到瞬态问题的表达式如式（６）所示Ｌ景（Ｌｐ两ｙ）出＋ＬＬ，７．（胛砧）以出＝
Ｊ＆Ｊ４，７・（厂ｇｒａｄ缈）以ｄｆ＋Ｊ出Ｊｒ品ｄｙ
（６）
万方数据
式中，Ａ为润滑油的面积，出为时间变化量。式（６）
是ＣＦＤ软件Ｆｌｕｅｎｔ进行瞬态计算的理论基础。
Ｆｌｕｅｎｔ软件求解压力．速度耦合方程的基本算法
是ＳＩＭＰＬＥ（Ｓｅｍｉ—ｉｍｐｌｉｃｉｔｍｅｔｈｏｄ
ｆｏｒｐｍｓｓｕｒｅｌｉｎｋｅｄ
ｅｑｕａｔｉｏｎｓ）算法。ＳＩＭＰＬＥ算法需要预定义流场压力。初始压力和参考压力的设置对计算结果影响比较大。尤其是动网格瞬态计算过程中出现回流现象时，预先建立流场压力分布对计算的收敛性至关重要。
对比工程实践，将进口边界条件设置为供油压力，参考压力为ｌＭＰａ，出油压力为０ＭＰａ，即泄油端与大气连通。
当采用动网格模型处理非定常压力．速度耦合
流动问题时，一般采用ＰＩＳＯ（Ｐｒｅｓｓｕｒｅ—ｉｍｐｌｉｃｉｔ
ｓｐｌｉｔｔｉｎｇｏｆ
ｏｐｅｒａｔｏｒｓ）算法，其计算基础是ＳＩＭＰＬＥ
算法。两者的区别在于预设压力求解动量方程和连
续方程后，ＰＩＳＯ算法需要连续两次修正压力和速度
方程，使压力和速度分别达到该时刻的两阶和三阶的精度。
鉴于轴颈高速对流场的影响，压力方程的插值
离散格式选择ＰＲＥＳＴＯ，动量方程的插值离散格式为ＱＵＩＣＫ。ＱＵＩＣＫ为流动的二次迎风插值格式，特别适于与流动方向对齐的结构网格，比二阶迎风
格式具有更精确的计算结果。其广义的表达式如（７）
所示
纯＝吉（击伽＋击啦）＋
ｉｌ丽％一百西细Ｊ
７
ｒ＆＋２＆
＆
、
（７）
…
式中，下标Ｐ为控制体积，下标Ｅ和形分别表示东
西方向与Ｐ相邻的控制体积；＆，＆，＆分别为
形，Ｐ，Ｅ相对应的控制体积的边长。伽，仇，伽
分别为Ｐ，Ｅ，形控制体积的通量，纯为Ｐ控制
体积西面的通量。
求解压力一速度耦合方程，得到油膜的三维物理
场。通过自定义程序循环指令遍历轴颈面域各节点的压力和黏性力并积分，得到不同方向的油膜力。当轴颈平衡位置变动过大且转速过高时，导致油膜破裂出现负压力，即油膜空化现象。计算油膜压力时，常用的处理方法是将负值区域压力置零。该处理方法和油膜力的计算所涉及的用户定义函数主要包括宏指令ＤＥＦＩＮＥＯＮＤＥＭＡＮＤ、循环函数
ｂｅｇｉｎ＿ｆ＿ｌｏｏｐ优ｆ）和面积函数Ｆ＿ＡＲＥＡ（Ａ，Ｚ矿１４１。
１．３弹簧光顺动网格模型
动网格模型主要应用于解决区域几何形状随时间变化的流体动力学问题。计算动静压轴承刚度阻尼系数时，网格变形主要考虑轴承间隙内油膜随轴
颈作刚体运动而产生的变化。在速度扰动和位移扰
２０１２年１２月熊万里等：基于动网格模型的液体动静压轴承刚度阻尼计算方法
１２１
动作用下，油膜与轴颈和油膜与轴瓦之间的边界层
度。至于大变形情况，则需要采用自定义程序解决
网格扭曲问题。
网格和油膜内部网格均产生形变。
由刚度阻尼计算模型可知，位移扰动和速度扰动产生的油膜力可线性化处理，位移扰动和速度扰
１．４油膜网格扭曲抑制
油膜网格扭曲形成的根本原因在于轴颈的旋转
运动和所采用的弹簧光顺更新模型。在确定油膜网格更新模型前提下，解决网格扭曲问题的落脚点在
动导致油膜网格产生的变形量亦可同样处理。油膜
网格运动的计算模型采用的是基于线性假设的弹性光顺Ｓｍｏｏｔｈｉｎｇ模型，其核心思想是胡克定律，网格移动与油膜力的关系【ｌ驯如式（８）所示
于轴颈的旋转运动。其关键就在于将轴颈旋转运动
参量代替实际转动嵌入至压力．速度耦合方程，计算瞬态油膜力。
Ｅ
２军如（Ｍ—Ａｄｙ）知２丽（８）
卫
．
１
解决网格扭曲问题时综合考虑了轴颈的运动。轴颈运动控制程序的基本功能包括：轴颈旋转运动和平动（位移扰动和速度扰动）。轴颈运动控制程序与动网格弹簧光顺模型相结合，实现轴颈边界网格和油膜网格的运动，进而计算得到轴颈受扰动条件
下不同位置的瞬态油膜力。
式中，Ｅ为节点ｆ的油膜力；反和ｄ，分别为节点ｉ
及相邻节点，的位置；Ａｄｉ和Ａｄｊ为节点位置变化
量；知为节点ｆ与节点．，之间的弹性系数：堵为与
节点ｉ相邻的节点数。
轴颈受外载荷作用作刚体平动。轴颈边界上的节点按胡克定律移动，通过流体内部节点传递到轴瓦边界上。在传递过程中，轴颈和轴瓦边界上网格
轴颈运动控制程序建立在Ｆｌｕｅｎｔ软件提供的用
户自定义函数基础上。轴颈的速度扰动和位移扰动由宏ＤＥＦＩＮＥ
ＣＧ
ＭＯＴＩＯＮ实现；为解决轴颈旋转
的形变对流体域网格的拉伸和挤压变形产生重要作用。为减小网格的畸变程度，需要对轴颈和轴瓦边界上网格节点的位置更新。节点控制如式（９）所示
引起的油膜网格畸变问题，将轴颈旋转的动边界转
变静边界，即将轴颈实际运动转化为物理量嵌入迭代计算，该功能通过宏ＤＥＦＩＮＥＰＲＯＦＩＬＥ实现ｕ
４ｌ。
钟州＝刃＋ＰＭ７
（９）
为验证轴颈运动控制程序成功加载动边界条件
实现轴颈对油膜网格的挤压运动，同时将轴颈旋转
式中，行和ｎ＋ｌ用于标识相邻迭代时间步的位置；∥为移动节点的松弛因子。该参数用于表征相邻内部节点的移动距离对变形边界节点移动幅度的影响程度。轴颈的运动将改变整个流体变形区域，影响流体物理场的变化。
在流体内部，网格运动时各个节点之间产生挤压或拉伸作用。节点之间的影响程度需结合物理现
的动边界条件成功转化为静边界条件，有效避免轴
颈旋转引起的油膜网格扭曲。对典型轴承利用Ｇａｍｂｉｔ建模，采用Ｆｌｕｅｎｔ进行动网格计算。
动网格计算对象为２排ｌＯ孔动静压无腔径向轴
承，如图ｌ所示。轴承基本结构参数：轴承孔径３０ｍｌＴｌ，长径比０．４，轴瓦厚１０ｍｉｌｌ，节流孔径０．６
ｍｍ，
象在ＰＩＳＯ算法中选择，默认设置对一般物理现象
通常都是有效的。流体内部节点各方向受力需取得平衡，相应的节点将产生移动，网格节点位移更新的计算模型如式（１０）所示
胁
节流孔沿轴向和周向均匀分布，轴承半径间隙２０ｐｍ。润滑油为２号主轴油，密度为８１０ｋｇ／ｍ３，动力黏度为０．００３
８５Ｐａ・Ｓ。
Ａｄ．ｍ＋１＝』Ｆ∑％
』
厶“日“Ｊ∑岛从，
（１０）
式中，ｍ为迭代次数。在轴颈旋转和平动作用下，油膜网格被动地产生形变。由于轴颈边界网格运动的特殊性（旋转和平动），弹簧光顺模型要求油膜网格类型具有非规则特征。对于规则网格，弹簧光顺模型需要特殊的指令才可以调用，但仅限于油膜网
ａ｝伦向轴承示意图
（ｂ）ｆ寺向轴嫒
图１径向轴承示意和零件图
格的挤压平动，轴颈旋转带动边界节点转动，致使
油膜网格形状扭曲，从而使计算发散。在小变形（挤压和扭曲）情况下，采用光顺模型具有较好的计算精万方数据
采用前处理软件Ｇａｍｂｉｔ建模。考虑到计算结果的精确性和网格数量等要求，流体区域采用六面体结构网格划分，使流体的流动迹线与网格分布的脉
机械工程学报第４８卷第２３期
络相适应，如图２所示。油膜厚度为１０层，油膜表面采用四边形划分，面网格采用Ｐａｖｅ模式，体网格采用Ｃｏｏｐｅｒ模式，流体域网格节点数量约为６０万。同时，设置相应的边界条件。节流孔定义为压力进
分布，偏布现象明显。这表明，在…定的偏心条件下，油压沿轴颈转动方向加强，且受轴颈转速影响明显，油膜压力场偏布表明轴颈转速已经成功加载。
口，轴承两端间隙为压力出口，轴颈和轴瓦为壁面。
置力ｐ／Ｐａ
１３０ＸｌＯ‘
ａ１流体域整体网络
｛ｈＩ局部放人示意图
（ａ）迭代１次
Ｊｔｉ力ｐ—Ｐ
蘸｛瓣
（ｂ）迭代３次
２’
ｈ－
图２滑动轴承网格模型
●３幽
１
采用Ｆｌｕｅｎｔ软件求解模型。求解器模型选择为非定常模式，压力．速度耦合方程的离散和计算模型
如前所述。边界条件设置：供油压力为３ＭＰａ，出油压力为０ＭＰａ。轴颈修改为运动壁面，其旋转速
度、扰动位移和扰动速度由轴颈运动控制程序实现，
轴瓦为静止壁面。轴承两侧出油面和间隙内流体均定义为变形区域。
轴颈的初始位置与轴承同心。调用自定义程序，对轴颈施加一个ｙ方向的速度扰动２ｌａｒｎ／ｓ，使轴颈向下运动。轴颈转速６ｋｒ／ｍｉｎ，单个迭代时间移动位移为２ｌａｍ。计算模型设置如前所述，初始化流场后，迭代６次，三维油膜压力分布和网格挤压变形
（ｃ）迭代６次
图３油膜网格受挤压变形及油膜压力分布
自定义的轴颈运动控制程序可实现轴颈的位移
扰动、速度扰动和旋转等功能，并有效避免油膜网
格的扭曲问题，从而为计算轴承油膜刚度和阻尼奠定了基础。然而，采用动网格方法计算面临一个重要问题：位移扰动和速度扰动均使轴颈产生位移，使油膜力产生变化。因此，界定位移扰动和速度扰动的取值范围就成为动网格计算油膜刚度和阻尼的又一关键问题。
１．５刚度阻尼的扰动范围界定
根据摄动理论，油膜刚度是指轴颈在平衡位置受位移扰动时承载力的变化，涉及轴颈转速和扰动位移两个物理量；油膜阻尼是指轴颈在平衡位置受速度扰动时承载力的变化，涉及轴颈转速和扰动速度两个物理量。利用雷诺方程计算时采用了虚位移
情况如图３所示（每个子图的上部表示油膜受轴颈
挤压后的压力分布，下部为受挤压后的油膜网格变形）。
轴颈平移和旋转运动均对油膜网格形变产生影响。由于采用了宏ＤＥＦＩＮＥＰＲＯＦＩＬＥ将轴颈旋转动边界转变为静边界，将运动方程嵌入动网格计算，避免了油膜网格的扭曲变形。油膜网格只在轴颈平
动方向产生挤压变形。对照不同迭代次数的网格变形情况，发现油膜网格只在轴颈平动方向产生挤压
变形，轴颈的旋转并未引起网格畸变。这表明将轴颈旋转动边界转换为静边界能避免轴颈旋转引起的网格畸变，同时说明宏ＤＥＦＩＮＥＣＧ—ＭＯＴＩＯＮ可控制轴颈的速度扰动和位移扰动。
原理，而动网格方法则反映出轴颈的真实运动。采用动网格方法计算刚度和阻尼的关键问题在于辨析
扰动位移与扰动速度对油膜力影响的内在联系。
在平衡位置，速度扰动和位移扰动对油膜承载
轴颈位移平动和旋转运动对油膜流场的影响可通过压力分布证实。轴颈向下移动和轴颈旋转均影
响油膜压力场的分布，不同之处在于轴颈旋转引起的动压效应使压力场偏布，位移平动使压力场对称分布。对照不同迭代时间的压力场可知，迭代１次时最高压力为３ＭＰａ：迭代６次时最高压力为
３．２
力的影响决定了计算刚度和阻尼时速度扰动和位移
扰动的取值范围。在瞬态计算和稳态计算的基础上，研究扰动位移与扰动速度对油膜力影响的差异。通
过对比研究动网格瞬态计算和稳态计算油膜力在平衡点的变化，分析扰动速度和扰动位移对油膜力的
影响，确定计算刚度阻尼时扰动位移和扰动速度的取值范围。
ＭＰａ，对比弦平面两侧最小间隙处油膜压力场万方数据
２０１２年１２月
熊万里等：基于动网格模型的液体动静压轴承刚度阻尼计算方法
对前面建立的计算模型，调用自定义程序使轴颈朝７方向运动，计算轴颈在不同扰动速度和扰动
位移下Ｘ和Ｙ方向的油膜力，并与稳态计算结果比较（扰动速度和扰动位移均为零）。
每个迭代时间步长，轴颈扰动位移为１ｇｍ，轴
颈在不同扰动速度ｖ作用下，ｘ和Ｙ两个方向分力与稳态计算结果对比如图４所示。
荟ｅＲ惫
涮懈簦
惹
偏心率ｓ
偏心率ｓ’
（ａ）ｘ方向油膜力的变化
（ｂ）ｙ方向油膜力的变化
图４扰动速度对ｘ和Ｙ方向油膜承载力的影响
由图４ａ可知，扰动速度对Ｘ方向分力的影响不
大。与稳态计算结果较靠近的速度为２．０ｇｒｎ／ｓ１和
２０．０
ｐＩＩｌ／ｓ１。扰动速度为２０
ｐ州ｓ１时，最大差距
１２．５％，平均差距７．７％；扰动速度为２．０ｇｒｒｄｓｌ时，平均差距７．５％。速度相差１０倍的情况下，平均差
距相差不大。
扰动速度对Ｙ方向分力的影响比较明显。由图
４ｂ可知，扰动速度为０．２ｍｍ／ｓ１时，最大差距为
２４．５％，平均差距１６％。扰动速度为２０．０ｎＩｌｌ／ｓ１时，
最大差距２．５％，平均差距１．８％；速度为２．０ｇｍ／ｓ１时，最大差距０．２％。
由上可知，在同样的扰动速度和扰动位移条件
下，采用动网格方法计算的Ｘ方向和Ｙ方向分力与稳态计算结果的差距不相同，挤压运动方向的差距大于与其垂直方向的差距，这与压力．速度耦合方程的求解模型有关。稳态计算和瞬态计算的基本润滑方程不同（与时间关联），求解模型不同，计算结果必然存在偏差。但从总体上看，在单个迭代步长内，扰动位移为１．０ｕｍ时，扰动速度不大于２０肛ｒｎ／ｓ１
的计算结果与稳态计算基本吻合。这个速度可视为
刚度和阻尼计算的分界点。
为减低扰动速度对油膜力的影响，单个迭代步长扰动速度取Ｏ．２“ＩＩｌ／ｓ１，在不同扰动位移Ｓ作用下，
Ｘ和Ｙ方向分力与稳态计算对比如图５所示。
对于ｘ方向分力的影响，扰动位移为２．０ｕｍ时，最大差距为２５．５％，平均差距为１９．５％；扰动位移为０．５ｇｍ时，最大差距１３．１％，平均差距７．８％，而位移为１．０“ｍ时平均差距为５．７％。
万方数据
冬
曩
茫
基
８
偏心率￡偏心率占
（ａ）ｘ方向油膜力的变化
（ｂ）ｙ方向油膜力的变化
图５位移扰动对Ｘ和Ｙ方向承载力的影响
对于Ｙ方向分力的影响，位移扰动为２．０“ｍ时，最大差距为１１％，平均差距为６．４％；扰动位移为０．５肛ｍ时，最大差距１．３％，平均差距０．６％，而位移为１．０ｇｍ时最高差距为０．３％。
对于扰动位移对Ｘ和Ｙ方向油膜力的影响与扰
动速度不同，Ｙ方向的差距小于ｘ方向的差距。由图５可知，位移为１．０ｇｍ与０．５９ｍ时，油膜力较为接近。Ｘ方向的差距为２．１％，Ｙ方向的差距为０．３％。
这种情况可视为刚度计算。
由上可知，扰动速度诞２０ｇｍ／ｓ１，扰动位移ｓ＜ｌ
ｇｍ，扰动速度和扰动位移对油膜力的影响较小，
平均差距不大于８％。动网格的瞬态计算与稳态
计算基本相当，可视为刚度计算。扰动速度
ｖ＞０．２
ｍｍ／ｓ１时，速度对油膜力的影响比较明显。动
网格的瞬态计算不能近似为稳态计算，可视为阻尼计算。在合适的速度和位移取值范围内，采用动网格计算轴承刚度和阻尼系数是切实可行的。
２典型轴承的动网格计算
采用动网格方法计算轴承的刚度阻尼系数包括
两个部分：一是计算主轴在轴承间隙内的平衡位置；
二是在平衡位置处计算轴承的刚度阻尼。
轴承结构如图１所示，轴承设计参数和润滑油
如前述。主轴的质量为５．５４蚝。建模时取轴承与轴
颈的同心位置为初始位置。重力作用下的最大移动速度为２．０ｇｍ／ｓ１，最大位移为１．０ｐｍ。计算主轴平衡位置时理论上应满足油膜力水平分力为０，垂直分力等于５５．４Ｎ。考虑到计算效率，误差为０．００１
Ｎ
时停止迭代。主轴的偏位角和偏心距随转速的变化
如图６所示。
计算刚度阻尼时需偏离平衡位置一定距离以适
应轴颈的运动。采用Ｇａｍｂｉｔ建模时需要将轴颈往ｘ或Ｙ方向偏移１．０ｇｍ。计算轴承刚度时，扰动速度为２．０ｇｍ／ｓ１，扰动位移１．０“ｍ；计算阻尼时，扰动速度０．５ｍｍ／ｓ１，扰动位移１．０ｇｍ，分别计算出轴承
机械工程学报第４８卷第２３期
ｆ
￥
吣
３
提翅ｇ
２
￡Ｔ≈魁尊弯
Ｏ
０
转速ｎ／（ｋｒ／ｍｉｎ）转速ｎ／（ｋｒ／ｍｉｎ）（ａ）转速对偏位角的影响
（ｂ）转速对偏心距的影响
图６转速对主轴平衡位置的影响
在速度扰动和位移扰动条件下的承载力，并利用式
（２）计算出轴承的刚度和阻尼系数，刚度阻尼系数随
转速的变化如图７、８所示。
堇１
要
＋＾ｗ—●—＾、．ｘ
｛
童
嗯。．’—日；．＋女Ｊｙ
弓
≥１０
墅
７＼√７’
蔷５
螽
曼
ｎ＇
７
Ⅻ
（ａ）转速对直接刚度的影响【ｂ）转速对交叉刚度的影响
图７转速对油膜刚度的影响
§０
（ａ）转速对直接阻尼的影响
馨
（ｂ）转速对交叉阻尼的影响
图８转速对油膜阻尼的影响
转速对油膜刚度和阻尼的影响各不相同。随着转速由１．５增至１２．０ｋｒ／ｍｉｎ，直接刚度略有下降，但交叉刚度迅速增加，由１．６增至１２．０Ｎｌｌｍ＿ｌ；直接阻尼呈下降趋势，交叉阻尼在０Ｎ／１．ｔｍ附近波动。采用动网格方法计算得到刚度阻尼随转速的变化趋势，与采用雷诺方程计算的变化趋势基本相同ｐＪ。
３试验分析
３．１试验方案
为验证采用动网格方法计算刚度和阻尼系数的可靠性，对典型轴承的刚度系数进行试验测试（刚度和阻尼由同一种方法计算得到，仅对刚度进行试验测试）。测试平台包括：双排ｌＯ孔液体动静压轴承电主轴，ＪＣ．３０００型油箱冷却器，ＶＦＤ．Ｂ型变频器，北京派莱博科技有限公司开发的振动及动态信号采集分析系统。
试验原理如图９所示。在主轴的正交平面０和
Ｙ方向，坐标与轴承一致）布置两个电容传感器，将万方数据
主轴对传感器的位移变化直接转换为电压信号，通
过电容测微仪和数据采集卡转换得到位移的变化。
主轴的转速由变频器控制。图１０为测量装置。
变频器／电器柜
三相电源几：；左丌＼９量ｌ，卜——强地７鼍兰兰｛；当
仁、ｆ≤；刁
手
液压站／冷却泵
嚣跫盟＾岸峙垤测头。人粤ｎ厂—丝叶仙
—————＜，
瘦体悬浮电生轴ｌ
、、
、匠！ｒ
——／／数据采集卡试验原理图
｛铨麓蝥０
Ｉｋ如＋七坩缸＝
Ⅲ，”…
｛嚣麓嚣
【‰勾＋‰勘２６
㈣，
”…
仁粪薹套㈣卜羔ｋ２蕊¨叫
２０１２年１２月
熊万里等：基于动网格模型的液体动静压轴承刚度阻尼计算方法
１２５
个速度点空运行半小时后，在１和＿ｙ方向先后施加
２０
Ｎ和４０Ｎ作用力，每个作用力加载３次，记录
每次加载力时ｘ和Ｙ方向的位移偏移量，取平均值作为该速度点和作用力的有效位移。试验测试出４８
组数据，将结果代入式（１３），计算出轴承的刚度系数。
３．３理论与试验对比分析
将动网格数值计算结果（图７）与试验测试结果进行比对，刚度随转速的变化趋势如图１１所示。
畲１
毒
言１
Ｚ
÷
￥
邑
曹
≥
划
刨一：麓彳—．＿Ｆ＝ＯＮ
二，ｒｔ—，，芑
—卜Ｐ＿２０Ｎ
毫
萎辎
辎柳
０
枷
转速ｎ／（ｋｒ／ｍｉｎ）转速ｎ／（ｋｒｌｍｉｎ）
（ａ）直接刚度￡“
（ｂ）直接刚度锄
言
毒Ｚ
￥
毒
魁一
蔑
以
权一
转速ｎ／（ｋｒ／ｍｉｎ）转速ｎ／（ｋｒｌｍｉｎ）（ｃ）交叉刚度≈矽（ｄ）交叉刚度盘ｐ
图１１
动网格计算结果和试验测试结果对比
通过对比动网格数值计算Ｆ＝ＯＮ与试验测试
Ｆ＝２０
Ｎ和Ｆ＝４０Ｎ平均值可知，两者偏差不大。直
接刚度ｋ最大偏差为１０．４％，平均偏差６．１６％；直接刚度‰的最大偏差２０．９％，平均偏差１３．３％；交
叉刚度‰最大偏差１５％，平均偏差９．８％；交叉刚
度ｋ与‰相差不大。
对于局部各数据点，动网格计算与试验测试结果并非完全一致，这与各自的外部因素有关。动网格计算偏向理论分析，其结果与求解器模型、网格
质量、边界条件和动网格模型以及参数设置有关；
试验测试侧重于工程实践，其结果受转子质量动不平衡、轴承加工误差、轴系装配误差、转速的脉动
和供油的平稳性以及加载力的平稳性和方向的平行
度等因素影响，以至于Ｆ＝２０Ｎ和Ｆ＝４０Ｎ两组数据
也存在一定差距。对于总体变化趋势，动网格计算结果与试验测试数据基本吻合，平均误差不超过
１２．５％，体现了典型轴承动力特性系数的内在一致性。
４结论
（１）采用自定义程序实现轴颈位移扰动、速度扰动、旋转和油膜力计算等功能，并利用弹簧光顺万方数据
模型更新网格的挤压变形，计算了轴承的刚度和阻尼系数。通过试验测试轴承的刚度系数，发现两者
基本吻合，表明采用动网格模型计算动静压轴承刚
度阻尼的数值方法是可行和有效的，从而为动静压
轴承动力特性系数计算提供了一种新方法。
（２）将轴颈的旋转动边界条件转换为静边界条件，并将转速加载在静边界上，有效避免了因轴颈旋转引起的油膜网格扭曲畸变，为旋转机械流体的动边界条件加载提供了一种新途径。
（３）对比研究了扰动速度和扰动位移对油膜力
的影响规律，界定了计算刚度阻尼系数时扰动速度
和扰动位移的合理取值范围，解决了动网格技术计
算轴承刚度阻尼系数的关键问题。
（４）再现了轴颈受位移扰动和速度扰动情况下油膜厚度的动态发展过程，反映了三维流场压力分布的时变特性，使计算过程更符合工程实际。
参考文献
［１】熊万里，林厚波，吕浪，等．超高速磨削电机内装式液
体悬浮电主轴关键技术及研发［ＥＢ／ＯＬ］．【２０１Ｉ－１２—２９］．
ｈｔｔｐ：／／ｗｗｗ．ｈｙｄｒｏ－ｓｐｉｎｄｌｅ．ｃｏｍ．
ＸＩＯＮＧ
Ｗａｎｌｉ，ＬＩＮＨｏｕｂｏ，ＬＯＬａｎｇ，ｅｔ
ａ１．Ｒｅｓｅａｒｃｈａｎｄ
ｄｅｖｅｌｏｐｍｅｎｔ
ｏｎ
ｋｅｙｔｅｃｈｎｏｌｏｇｙｏｆｍｏｔｏｒｂｕｉｌｔ・ｉｎｓｐｉｎｄｌｅ
ｓｕｐｐｏｒｔｅｄｂｙｈｙｄｒｏｄ’７ｎａｍｉｃ
ａｎｄ
ｈｙｄｒｏｓｔａｔｉｃｂｅａｒｉｎｇｓｆｏｒ
ｕｌｔｒａ—ｈｉｇｈｓｐｅｅｄｇｆｉｎｄｉｎｇ［ＥＢ／ＯＬ］．［２０１１－１２—２９］．ｈｔｔｐ：／／ｗｗｗ．ｈｙｄｒｏ－ｓｐｉｎｄｌｅ．ｃｏｍ．
［２】ＬＵＮＤ
Ｊ
Ｗ．Ｔｈｅ
ｓｔａｂｉｌｉｔｙｏｆ
ａｌｌ
ｅｌａｓｔｉｃｒｏｔｏｒｉｎ
ｊｏｕｒｎａｌ
ｂｅａｒｉｎｇｗｉｔｈｆｌｅｘｉｂｌｅｄａｍｐｅｄｓｕｐｐｏｒｔｓ［Ｊ］．ＡＳＭＥ
Ｊｏｕｒｎａｌ
ｏｆＡｐｐｌｉｅｄＭｅｃｈａｎｉｃｓ，１９６５，３２（４）：９１１－９２０．
［３】ＧＨＯＳＨＭＫ，ＭＡＪＵＭＤＡＲＢＣ，ＲＡＯＪＳ．Ｓｔｅａｄｙ・ｓｔａｔｅ
ａｎｄｄｙｎａｍｉｃｂｅｈａｖｉｏｕｒｏｆｍｕｌｔｉ－ｒｅｃｅｓｓｈｙｂｒｉｄｏｉｌｊｏｕｒｎａｌ
ｂｅａｒｉｎｇ［Ｊ］．Ｊｏｕｒｎａｌ
ｏｆＭｅｃｈａｎｉｃａｌＥｎｇｉｎｅｅｒｉｎｇＳｃｉｅｎｃｅ，
１９７９，２１（５）：３４５－３５１．
［４】ＲＯＷＥＷ
Ｂ，ＣＨＯＮＧＦＳ．Ｃｏｍｐｕｔａｔｉｏｎｏｆｄｙｎａｍｉｃ
ｆｏｒｃｅｃｏｅｆｆｉｃｉｅｎｔｓｆｏｒ
ｈｙｂｒｉｄ（ｈｙｄｒｏｓｔａｔｉｃ／ｈｙｄｒｏｄｙｎａｍｉｃ）
ｊｏｕｒｎａｌｂｅａｒｉｎｇｓｂｙｔｈｅｆｉｎｉｔｅ
ｄｉｓｔｕｒｂａｎｃｅ
ａｎｄｐｅｒｔｕｒｂａｔｉｏｎ
ｔｅｃｈｎｉｑｕｅｓ［Ｊ］．ＴｒｉｂｏｌｏｇｙＩｎｔｅｒｎａｔｉｏｎａｌ，１９８６，１９（５１：
２６０．２７１．
［５】许尚贤，陆永颐，万鹏．小孔式动静压轴承的特性分析
与实验研究［Ｊ】．机械工程学报，１９９３，２９（５）：８８—９５．
ＸＵ
Ｓｈａｎｇｘｉａｎ，ＬＵＹｏｎｇｙｉ，ＷＡＮ
Ｐｅｎｇ．Ａｎａｌｙｓｉｓａｎｄ
ｅｘｐｅｒｉｍｅｎｔａｌｉｎｖｅｓｔｉｇａｔｉｏｎ
ｏｆ
ａ
ｈｏｌｅ－ｅｎｔｒｙ
ｈｙｂｒｉｄ
ｂｅａｒｉｎｇ［Ｊ］．Ｃｈｉｎｅｓｅ
ＪｏｕｒｎａｌｏｆＭｅｃｈａｎｉｃａｌＥｎｇｉｎｅｅｒｉｎｇ，
１９９３，２９（５）：８８－９５．
［６］６郑铁生，许庆余．滑动轴承动特性系数新算法【Ｊ】．机械
１２６
机械工程学报
第４８卷第２３期
工程学报，１９９６，３２（３）：２０．２７．
ＺＨＥＮＧＴｉｅｓｈｅｎｇ，ＸＵＱｉｎｇｙｕ，Ｎｅｗ
ｎｕｍｅｒｉｃａｌｍｅｔｈｏｄｆｏｒ
ｃａｌｃｕｌａｔｉｎｇ
ｄｙｎａｍｉｃ
ｃｏｅｆｆｉｃｉｅｎｔｓ
ｏｆｏｉｌ－ｌｕｂｒｉｃａｔｅｄ
ｊｏｕｒｎａｌｂｅａｒｉｎｇｓ［Ｊ】．Ｃｈｉｎｅｓｅ
Ｊｏｕｒｎａｌｏｆ
Ｍｅｃｈａｎｉｃａｌ
Ｅｎｇｉｎｅｅｒｉｎｇ，１９９６，３２（３）：２０－２７．
［７】ＳＡＷＩＣＫＩＪＴ，ＲＡＯＴＶＶＬＮ．Ａｎｏｎｌｉｎｅａｒｍｏｄｅｌｆｏｒ
ｐｒｅｄｉｃｔｉｏｎ
ｏｆｄｙｎａｍｉｃ
ｃｏｅｆｆｉｃｉｅｎｔｓ
ｉｎ
ａ
ｈｙｄｒｏｄ）‘ｎａｍｉｃｊｏｕｒｎａｌｂｅａｒｉｎｇ［Ｊ］．Ｉｎｔｅｒｎａｔｉｏｎａｌ
Ｊｏｕｒｎａｌ
ｏｆ
Ｒｏｔａｔｉｎｇ
Ｍａｃｈｉｎｅｒｙ，２００４，１０（６）：５０７－５１３．［８］ＭＯＲＩＡ，ＭＡＫＩＮＯＴ，ＭＯＲＩＨ．Ｅｎｔｒｙ
ｆｌｏｗａｎｄｐｒｅｓｓｕｒｅ
ｊｕｍｐ
ｉｎ
ｓｕｂｍｅｒｇｅｄｍｕｌｔｉ－ｐａｄ
ｂｅａｒｉｎｇｓ
ａｎｄ
ｇｒｏｏｖｅｄ
ｂｅａｒｉｎｇｓ［Ｊ】．ＪｏｕｒｎａｌｏｆＴｒｉｂｏｌｏｇｙ，１９９２，１１４（２）：３７０－３７７．［９］ＴＵＣＫＥＲ
Ｐ
Ｇ，ＫＥＯＧＨ
ＰＳ．Ｏｎｔｈｅｄｙｎａｍｉｃｔｈｅｒｍａｌ
ｓｔａｔｅｉｎ
ａ
ｈｙｄｒｏｄｙｎａｍｉｃｂｅａｒｉｎｇｗｉｔｈ
ａ
ｗｈｉｒｌｉｎｇ
ｊｏｕｍａｌ
ｕｓｉｎｇＣＦＤｔｅｃｔｍｉｑｕｅｓ［Ｊ】．ＪｏｕｒｎａｌｏｆＴｒｉｂｏｌｏｇｙ，１９９６，
１１８（２）：３５６－３６３．［１０］ＳＡＨＬＩＮ
Ｆ，ＧＬＡＶＡＴＳＫＩＨＳ
Ｂ，ＡＬＭＱＵＩＳ
Ｔ，ｅｔａ１．Ｔｗｏ
ｄｉｍｅｎｓｉｏｎａｌＣＦＤａｎａｌｙｓｉｓｏｆｍｉｃｒｏ．ｐａＲｅｍｅｄｓｕｒｆａｃｅｓｉｎ
ｈｙｄｒｏｄｙｎａｍｉｃｌｕｂｒｉｃａｔｉｏｎ［Ｊ］．ＪｏｕｒｎａｌｏｆＴｒｉｂｏｌｏｇｙ，２００５，
１２７（１）：９６—１２７．
万方数据
ＧＥＲＴＺＯＳ
Ｋ
Ｐ，
ＮＩＫＯＬＡＫＯＰＯＵＬＯＳ
Ｐ
Ｇ，ＰＡＰＡＤＯＰＯＵＬＯＳＣＡ．ＣＦＤａｎａｌｙｓｉｓｏｆｊｏｕｒｎａｌ
ｂｅａｒｉｎｇ
ｈｙｄｒｏｄｙｎａｍｉｃ
ｌｕｂｒｉｃａｔｉｏｎ
ｂｙＢｉｎｇｈａｍｌｕｂｒｉｃａｎｔ［Ｊ］．
Ｔｒｉｂｏｌｏｇｙ
Ｉｎｔｅｒｎａｔｉｏｎａｌ，２００８，４１（１２）：１
１９０－１２０４．
［１２】王福军．计算流体动力分析－ＣＦＤ软件原理与应用［Ｍ】．
北京：清华大学出版社，２００４．ＷＡＮＧＦｕｊｕｎ．Ｃｏｍｐｕｔａｔｉｏｎａｌ
ｆｌｕｉｄ
ｄｙｎａｍｉｃｓ
ａｎａｌｙｓｉｓ－ｔｈｅ
ｐｒｉｎｃｉｐｌｅｓａｎｄａｐｐｌｉｃａｔｉｏｎｓｏｆＣＦＤｓｏｔｔｗａｒｅ［Ｍ］．Ｂｅｉｊｉｎｇ：
ＴｓｉｎｇｈｕａＵｎｉｖｅｒｓｉｔｙ
Ｐｒｅｓｓ，２００４．
［１３］
ＦＬＵＥＮＴ．ＦＬＵＥＮＴ６．３
ｕｓｅｒ’ｓ
ｇｕｉｄｅ［Ｍ］．Ｌｅｂａｎｏｎ：Ｎｅｗ
Ｈａｍｐｓｈｉｒｅ（ＵＳＡ）Ｆｌｕｅｎｔ
Ｉｎｃ．，２００６．
［１４］
ＦＬＵＥＮＴ．ＦＬＵＥＮＴ６．３
ＵＤＦ
ｍａｎｕａｌ［Ｍ］．Ｌｅｂａｎｏｎ：Ｎｅｗ
Ｈａｍｐｓｈｉｒｅ（ＵＳＡ）ＦｌｕｅｎｔＩｎｃ．，２００６．
［１５］
张直明．滑动轴承的流体动力润滑理论［Ｍ］．北京：高等教育出版社，１９８６．
ＺＨＡＮＧ
Ｚｈｉｍｉｎｇ．Ｌｕｂｒｉｃａｔｉｏｎｔｈｅｏｒｙ
ｆｏｒ
ｓｌｉｄｅ
ｂｅａｒｉｎｇ［Ｍ］．Ｂｅｉｊｉｎｇ：Ｈｉｇｈｅｒ
ＥｄｕｃａｔｉｏｎＰｒｅｓｓ，１９８６．
作者简介：熊万里（通信作者），男，１９７１年出生，博士，教授，博士研
究生导师。主要研究方向为高速电主轴系统动力学及超高速磨削装备。
Ｅ－ｍａｉｌ：ｗａｎ３６９＠ｖｉｐ．ｓｉｎａ．００ｍ
基于动网格模型的液体动静压轴承刚度阻尼计算方法
作者：作者单位：
熊万里， 侯志泉， 吕浪， 阳雪兵， 袁巨龙， XIONG Wanli， HOU Zhiquan， L(U) Lang， YANG Xuebing ， YUAN Julong
熊万里,侯志泉,吕浪,XIONG Wanli,HOU Zhiquan,L(U) Lang(湖南大学国家高效磨削工程技术研究中心 长沙410082)， 阳雪兵,YANG Xuebing(湘潭电机股份有限公司风能分公司 湘潭411102)， 袁巨龙,YUAN Julong(浙江工业大学机械制造及自动化教育部重点实验室 杭州310014)机械工程学报
JOURNAL OF MECHANICAL ENGINEERING2012,48(23)
刊名：英文刊名：年，卷(期)：
1. 熊万里;林厚波;吕浪 超高速磨削电机内装式液体悬浮电主轴关键技术及研发 2011
2. LUND J W The stability of an elastic rotor in journal bearing with flexible damped supports 1965(04)3. GHOSH M K;MAJUMDAR B C;RAO J S Steady-state and dynamic behaviour of multi-recess hybrid oil journal beating 1979(05)
4. ROWE W B;CHONG F S Computation of dynamic force coefficients for hybrid (hydrostatic/hydrodynamic)journalbearings by the finite disturbance and perturbation techniques 1986(05)5. 许尚贤;陆永颐;万鹏 小孔式动静压轴承的特性分析与实验研究 1993(05)6. 郑铁生;许庆余 滑动轴承动特性系数新算法 1996(03)
7. SAWICKIJT;RAOTVVLN A nonlinear model for prediction of dynamic coefficients in a hydrodynamic journalbearing [外文期刊] 2004(06)
8. MORIA;MAKINO T;MORI H Entry flow and pressure jump in submerged multi-pad bearings and grooved bearings1992(02)
9. TUCKER P G;KEOGH P S On the dynamic thermal state in a hydrodynamic bearing with a whirling journal usingCFD techniques[外文期刊] 1996(02)
10. SAHLIN F;GLAVATSKIH S B;ALMQUIS T Two dimensional CFD analysis of micro-patterned surfaces in hydrodynamiclubrication 2005(01)
11. GERTZOS K P;NIKOLAKOPOULOS P G;PAPADOPOULOS C A CFD analysis of journal beating hydrodynamic lubrication byBingham lubricant[外文期刊] 2008(12)
12. 王福军 计算流体动力分析-CFD软件原理与应用 200413. FLUENT FLUENT 6.3 user's guide 200614. FLUENT FLUENT 6.3 UDF manual 200615. 张直明 滑动轴承的流体动力润滑理论 1986
本文链接：http://d.g.wanfangdata.com.cn/Periodical_jxgcxb201223018.aspx
第４８卷第２３期２０１２年１２月
机械工程学报
ＪＯＵＲＮＡＬＯＦ
Ｖｂｌ．４８Ｄｅｃ．
ＮＯ．２３２０１２
ＭＥＣＨＡＮＩＣＡＬ
ＥＮＧＩＮＥＥＲＩＮＧ
ＤｏＩ：１０．３９０１／ＪＭ［Ｅ．２０１２．２３．１１８
基于动网格模型的液体动静压轴承
刚度阻尼计算方法木
熊万里１
侯志泉１
吕
浪１
阳雪兵２袁巨龙３
（１．湖南大学国家高效磨削工程技术研究中心长沙４１００８２；
２．湘潭电机股份有限公司风能分公司湘潭４１１１０２；
３．浙江工业大学机械制造及自动化教育部重点实验室杭州３１００１４）
摘要：针对雷诺方程不能反映高速油流周向惯性效应、轴颈周向动压效应和静压扩散效应之间的非线性耦合关系及其对三维流场特性的影响，进而导致动静压轴承刚度阻尼的计算精度难以满足高转速精密轴承设计需要之现状，提出基于纳维一斯托克斯方程的动网格计算轴承刚度阻尼新方法。该方法采用自定义程序实现轴颈的旋转速度、位移扰动和速度扰动以及扰动过程中油膜力的计算等功能，成功将轴颈的旋转动边界转换为静边界，有效避免因轴颈旋转引起的网格畸变，并利用弹簧光顺模型更新轴颈受速度扰动和位移扰动引起的油膜网格挤压变形。通过对比研究瞬态和稳态条件下油膜力的变化，界定计算刚度阻尼系数时速度扰动和位移扰动的合理取值范围。利用所提出的方法计算典型轴承的刚度阻尼系数，再现动静压轴承油膜厚度和三维流场压力分布的动态发展过程，并采用影响系数法对该轴承的刚度系数进行试验测试。通过对比试验数据和动网格计算结果发现两者基本吻合，表明所提出的计算方法是有效可行的。关键词：动静压轴承刚度系数阻尼系数动网格中图分类号：ＴＨｌ３３
Ｍｅｔｈｏｄ
ｆｏｒＣａｌｃｕｌａｔｉｎｇＳｔｉｆｆｎｅｓｓａｎｄＤａｍｐｉｎｇ
ＢｅａｒｉｎｇｓＢａｓｅｄ
ｏｎ
Ｃｏｅｆｆｉｃｉｅｎｔｓ
ＤｙｎａｍｉｃＭｅｓｈＭｏｄｅｌ
Ｌａｎ９１
ＹＡＮＧＸｕｅｂｉｎ９２
ｏｆＨｙｂｒｉｄ
ＸＩＯＮＧ
Ｗａｎｌｉｌ
ＨＯＵＺｈｉｑｕａｎｌ
ＬＯ
ＹＵＡＮＪｕｌｏｎ９３
（１．ＮａｔｉｏｎａｌＲｅｓｅａｒｃｈＣｅｎｔｅｒｆｏｒＨｉｇｈＥｆｆｉｃｉｅｎｃｙＧｒｉｎｄｉｎｇ，ＨｕｎａｎＵｎｉｖｅｒｓｉｔｙ，Ｃｈａｎｇｓｈａ４１００８２；２．ＷｉｎｄＰｏｗｅｒＳｕｂｓｉｄｉａｒｙＣｏ．，ＸｉａｎｇｔａｎＥｌｅｃｔｒｉｃＭａｎｕｆａｃｔｕｒｉｎｇＧｒｏｕｐＣｏ．Ｌｔｄ．，Ｘｉａｎｇｔａｎ４１１１０２；
３．ＴｈｅＫｅｙＬａｂｏｒａｔｏｒｙ
ｏｆＭｅｃｈａｎｉｃａｌＭａｎｕｆａｃｔｕｒｅａｎｄＡｕｔｏｍａｔｉｏｎｏｆＭｉｎｉｓｔｒｙ
ｏｆＥｄｕｃａｔｉｏｎ，
ＺｈｅｊｉａｎｇＵｎｉｖｅｒｓｉｔｙｏｆＴｅｃｈｎｏｌｏｇｙ，Ｈａｎｇｚｈｏｕ３１００１４）
Ａｂｓｔｒａｃｔ：Ｆｏｒｔｈｅ
ｒｅａｓｏｎ
ｍａｔＲｅｙｎｏｌｄｓｅｑｕａｔｉｏｎｉｓｎｏｔｓｕｉｔａｂｌｅｆｏｒｒｅｆｌｅｃｔｉｎｇｎｏｎｌｉｎｅａｒｆｌｏｗｍｅｃｈａｎｉｃｓｃｏｕｐｌｅｄｂｙｉｎｅｒｔｉａｌｅｆｆｅｃｔｏｆ
ｏｉｌｃｉｒｃｕｍｆｅｒｅｎｔｉａｌｆｌｏｗ，ｓｑｕｅｅｚｅｅｆｆｅｃｔｏｆｔｈｅｄａｍｐｉｎｇｃｏｅｆｆｉｃｉｅｎｔｓ
ｃａｌｌ
ｊｏｕｒｎａｌ．ａｎｄ
ｏｎ
ｄｉｆｆｕｓｉｏｎｅｆｆｅｃｔｏｆｏｉｌｅｎｔｒｙｆｌｏｗ，ｔｈｅｃａｌｃｕｌａｔｉｎｇａｃｃｕｒａｃｙｏｆｓｔｉｆｆｎｅｓｓａｎｄ
ｎｏｔ
ｓａｔｉｓｆｙｏｐｔｉｍａｌｄｅｓｉｇｎｏｆｈｉｇｈｓｐｅｅｄｂｅａｒｉｎｇｓ．Ａｎｅｗａｐｐｒｏａｃｈｏｆｃａｌｃｕｌａｔｉｎｇ
ｓｔｉｆｆｎｅｓｓ
ａｎｄｄａｍｐｉｎｇ
ｃｏｅｆｆｉｃｉｅｎｔｓｏｆｊｏｕｒｎａｌｂｅａｒｉｎｇｓｉｓｄｅｖｅｌｏｐｅｄｂａｓｅｄ
ｐｅｒｔｕｒｂａｔｉｏｎ，ｔｈｅｖｅｌｏｃｉｔｙｐｅｒｔｕｒｂａｔｉｏｎｂｏｕｎｄａｒｙ
ｏｆｔｈｅ
Ｎａｖｉｅｒ－Ｓｔｏｋｅｓｅｑｕａｔｉｏｎ．Ｔｈｅｃｈａｎｇｉｎｇｏｉｌｆｉｌｍｆｏｒｃｅ，ｔｈｅｄｉｓｐｌａｃｅｍｅｎｔ
ａｒｅ
ａｎｄ
ｔｈｅｒｏｔａｔｉｎｇｓｐｅｅｄａｃｈｉｅｖｅｄｂｙｕｓｉｎｇ
ｕｓｅｒ
ｄｅｆｉｎｉｔｉｏｎｆｕｎｃｔｉｏｎ
ｐｒｏｇｒａｍｓ．Ｔｈｅ
ｍｏｖｉｎｇ
ｊｏｕｒｎａｌｃａｕｓｅｄｂｙｒｏｔａｔｉｏｎａｌｍｏｖｅｍｅｎｔｉｓ
ｃｏｏｒｄｉｎａｔｉｏｎｏｆｍｏｖｉｎｇｇｒｉｄｃａｕｓｅｄｂｙｄｉｓｐｌａｃｅｍｅｎｔｔｈｅｏｒｅｔｉｃａｌｍｏｄｅｌ．Ｃｏｍｐａｒｉｎｇｏｉｌｆｉｌｍｆｏｒｃｅｐｅｒｔｕｒｂａｔｉｏｎ
ｃｈａｎｇｅｄｉｎｔｏｓｔａｔｉｃｂｏｕｎｄａｒｙｆｏｒ
ｐｅｒｔｕｒｂａｔｉｏｎａｎｄｖｅｌｏｃｉｔｙｐｅｒｔｕｒｂａｔｉｏｎｉｓ
ａｖｏｉｄｉｎｇｔｈｅｇｒｉｄｄｉｓｔｏｒｔｉｏｎ．Ｔｈｅｕｐｄａｔｅｄｂｙｔｈｅｓｐｒｉｎｇｓｍｏｏｔｈｉｎｇ
ｕｎｄｅｒｓｔｅａｄｙ
ｃｏｎｄｉｔｉｏｎｗｉｔｈｔｈｏｓｅｏｆｔｒａｎｓｉｅｎｔｐｒｏｃｅｓｓｅｓ，ｔｈｅｖａｌｕｅｓｃｏｐｅｏｆｄｉｓｐｌａｃｅｍｅｎｔ
ａｎｄｖｅｌｏｃｉｔｙ
ｐｅｒｔｕｒｂａｔｉｏｎａｒｅｄｅｆｉｎｅｄ．Ｔｈｅ
ｓｔｉｆｆｎｅｓｓ
ａｎｄｄａｍｐｉｎｇｃｏｅｆｆｉｃｉｅｎｔｓｏｆ
ａ
ｔｙｐｉｃａｌｂｅａｒｉｎｇａｒｅｃａｌｃｕｌａｔｅｄｂｙｐｒｏｃｅｓｓｏｆｏｉｌｆｉｌｍｔｈｉｃｋｎｅｓｓａｎｄ
ａｐｐｌｙｉｎｇｔｈｅａｐｐｒｏａｃｈｃｏｍｂｉｎｅｄｗｉｔｈｔｈｅｐｅｒｔｕｒｂａｔｉｏｎｔｈｅｏｒｙ，ａｎｄｔｈｅｄｙｎａｍｉｃｄｅｖｅｌｏｐｍｅｎｔ
ｔｈｒｅｅ－ｄｉｍｅｎｓｉｏｎａｌｐｒｅｓｓｕｒｅｄｉｓｔｒｉｂｕｔｉｏｎａｒｅａｃｈｉｅｖｅｄ．Ｆｉｎａｌｌｙ，ａｌｌｅｘｐｅｒｉｍｅｎｔｉｓｕｓｅｄｔｏｔｅｓｔｔｈｅｓｔｉｆｆｎｅｓｓｃｏｅｆｆｉｃｉｅｎｔｓｏｆｔｈｅｂｅａｒｉｎｇｂａｓｅｄ
ｏｎ
ｔｈｅｉｎｆｌｕｅｎｃｅｃｏｅｆｆｉｃｉｅｎｔｍｅｔｈｏｄ，ｗｈｉｃｈａｒｅｓｕｂｓｔａｎｔｉａｌｉｎａｇｒｅｅｍｅｎｔｗｉｔｈｔｈｅｏｒｅｔｉｃａｌａｎａｌｙｓｉｓ．ＴｈｅｒｅｓｕＲｓ
ｓｈｏｗｔｈａｔｔｈｅ
ｎｕｍｅｒｉｃａｌｍｅｔｈｏｄｉｓｖａｌｉｄ．
Ｋｅｙｗｏｒｄｓ：Ｈｙｂｒｉｄｊｏｕｒｎａｌｂｅａｒｉｎｇ
Ｓｔｉｆｆｎｅｓｓ
ｃｏｅｆｆｉｃｉｅｎｔｓＤａｍｐｉｎｇｃｏｅｆｆｉｃｉｅｎｔｓＤｙｎａｍｉｃｍｅｓｈ
・国家自然科学基金（５０９７５０８２，５１２７５１６３）、教育部新世纪优秀人才支持计划（ＮＣＥＴ－１０—０３６８）和湖南省自然科学基金资助项目。
２０１２０２２８收到初稿，２０１２０８２８收到修改稿
万方数据
２０１２年１２月熊万里等：基于动网格模型的液体动静压轴承刚度阻尼计算方法
１１９
０前言
液体动静压轴承以其高精度、高刚度、高阻尼吸振性和长寿命等优势，在超高速精密磨削领域获得广泛应用。目前，国家高效磨削工程技术研究中
心研制出用于磨削外圆面工件的超高速液体动静压
轴承电主轴，额定功率达３５ｋＷ，主轴最高转速达
９
ｋｒ／ｍｉｎ；用于磨削凸凹面工件的液体动静压轴承电
主轴，额定功率达１２ｋＷ，最高转速达３０
ｋｒ／ｍｉｎ【ｌＪ。
随着电主轴转速的进一步提高，电主轴系统因高速
失稳出现刮瓦、抱轴等现象，导致电主轴系统失效的可能性逐步增加。高速稳定性已成为当前超高速
精密电主轴系统动力学分析的重要课题，其难点就
在于支承轴承动态刚度和阻尼的精确计算。
动静压轴承的动态刚度阻尼反映了轴颈受外载
荷作用下振动位移和振动速度与油膜力之间的内在
联系。１９６５年，ＬＵＮＤＬ２ｊ采用４个刚度系数和４个
阻尼系数建立了挠性“轴承．转子系统”模型，从理论上解释了系统涡动失稳的形成原因，揭示了“轴承．转子系统”涡动失稳的物理本质。８个系数被广泛用于表征动静压轴承的动态特性和分析轴承的高
速稳定性。
计算轴承的刚度阻尼，常用的方法有差分法和偏导数法，其关键科学问题在于求解动态流体润滑方程得到油膜的压力分布。通常情况下，动态雷诺
方程作为润滑方程用于计算动静压轴承的刚度阻尼【ｊ√Ｊ。由于雷诺方程难以精确反映转速引起的油流周向惯性效应、轴颈环周向动态挤压效应和静压效
应之间非线性耦合关系及其对三维物理场（速度场、
温度场和压力场）的影响，油膜刚度阻尼的计算精度
受到制约，因此有必要直接以Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程为
基础，精确研究轴承结构、主轴转速、位移扰动和速度扰动等对油膜刚度和阻尼的影响规律。文献
【８—１０］采用基于Ｎ—Ｓ方程的计算流体动力学
（Ｃｏｍｐｕｔａｔｉｏｎａｌ
ｆｌｕｉｄ
ｄｙｎａｍｉｃｓ，ＣＦＤ）开展了静压轴承承载特性的研究工作，证实了Ｎ—ｓ方程在表征复
杂求解域流体流动形态方面可弥补雷诺方程的不
足。文献【１１］进一步引入了“可动边界”的概念，
采用动网格方法计算了动压轴承在外载荷作用下轴
承的压力分布，但尚未解决轴颈旋转带来的网格扭
曲问题。
鉴于ＣＦＤ在表征流体流动形态方面的优势及其计算轴承动态性能方面的不足，本文提出了基于
ＣＦＤ动网格模型的液体动静压轴承油膜刚度阻尼万方数据
系数的计算方法，该方法在求解瞬态Ｎ—Ｓ方程基础
上，采用自定义程序结合动网格模型，动态反映了轴颈转速、位移扰动和速度扰动等因素对油膜压力
场的影响，从而使油膜力计算结果与高速时实际工
况更吻合；通过对比分析稳态计算和瞬态计算油膜力的差异，界定了采用动网格方法计算刚度阻尼时扰动位移和扰动速度的合理取值范围。通过对比典型算例和试验测试结果，验证了动网格计算方法的
有效性和可行性。１
刚度阻尼动网格计算方法
基于动网格模型的轴承刚度阻尼计算方法是在
ＣＦＤ软件基础上，采用自定义程序实现轴颈扰动功能，结合差分计算模型，计算轴承油膜刚度阻尼的一种新方法。该方法利用ＣＦＤ软件求解Ｎ．Ｓ方程得到轴颈位置变动前后流体的瞬态压力场，通过自定义程序遍历轴颈表面压力得到不同方向的瞬态油膜
力；轴颈位置变动过程中出现的油膜网格变形由动
网格模型更新；轴颈位置变动的三个运动参量轴颈
转速、扰动位移和扰动速度由自定义程序指定，同时有效解决因轴颈旋转引起的油膜网格扭曲问题；
通过对比扰动位移和扰动速度对油膜力的影响，界
定计算刚度和阻尼的合理取值范围。
１．１
刚度阻尼的差分计算模型
轴颈在平衡位置受外载荷作用，对轴承油膜产
生位移扰动和速度扰动，油膜力将发生变化。在小
扰动条件下，轴承油膜刚度阻尼可采用差分法求解。根据摄动理论，将油膜力线性化处理，油膜力的变化与扰动位移和扰动速度之间的关系１４１如式（１）所示
式中，峨和峨为油膜力的变化；ｋ和Ｊｉ｝珂是直
∞］＝（乏ｋ％．Ａ㈣ａｙｊ＋（墓北）㈣接刚度；ｋ和ｋ为交叉刚度；Ｃｘｘ和．ｃ抑为直接阻尼；ｃ。和ｃ。为交叉阻尼；缸和每为扰动位移；ｘ７
和Ｙ７为扰动速度。
采用差分模型计算轴承的刚度和阻尼系数，需要考虑位移和速度对不同方向油膜力的影响。刚度
和阻尼的计算如式（２）所示
胁ｌｌ
％２言，蜴肛
印＝
％２亍刖‰
胁
＝
警堕妙％２亏
，甜’劬
细
＝
等堕少
呦２——２
心‰
ｙ。
（２）
式中，刖％为位移扰动引起油膜力的变化；刖％为速度扰动引起油膜力的变化；ｆ、，表示ｘ、Ｙ中的某
１２０
机械工程学报第４８卷第２３期
一个。
由式（２）可知，采用动网格方法计算油膜刚度阻
尼系数，需要解决三个问题：①计算刚度阻尼需要
计算轴颈在位移和速度变动前后瞬态油膜力，该油膜力可通过求解瞬态压力一速度耦合方程获得；
②轴颈受外载荷和自转影响，轴颈在轴承油膜中的
相对位置将产生变化，油膜网络挤压变形并由动网
格模型更新；③计算刚度阻尼需要在轴颈动边界条
件中引入轴颈旋转、速度扰动和位移扰动，轴颈运动控制由自定义程序实现。１．２瞬态油膜力计算
在ＣＦＤ软件中，油膜流动的基本动力润滑理论是Ｎ—Ｓ方程。数值求解润滑方程的基本思想是有限体积法，将流体区域离散为流体微元，计算出微元
的物理参量。其核心思想是子域法加离散法。
在稳定工况下，动静压轴承间隙内流体为黏性
层流，不考虑流体的粘温效应，动力润滑方程包括
连续性方程和动量方程，其瞬态通用表达形式【ｌ压¨Ｊ
如式（３）所示
ａ
鼍（ｐ缈）＋ｄｉＶ（ｐ“伊）＝ｄｉｖ（Ｆｇｒａｄ口ｏ）＋－品
（３）
式中，ｐ表示流体密度；Ｈ为流场速度矢量；缈为
通用变量；厂表示广义扩散系数；＆为广义源项。
式（３）各项分别表征流体的瞬态项、对流项、扩
散项和源项。对比广义雷诺方程的假设条件及其所包含的压力项、楔入项和挤压项，各项之间呈线性叠加关系。显然，ＣＦＤ控制方程更适于表征间隙内流体的三维环周向流动形态。
动静压轴承间隙内流体的流场形态主要由供油
压力产生的压力流和轴颈旋转引起的压力流和速度
流非线性耦合形成，表现为静压与动压的强耦合效
应。流体的压力．速度耦合方程如式（４）所示
÷（胛）＋ｄｉｖ（ｐｆｏｕ）＝ｄｉｖ（ｒｇｒａｄ伊）－署＋昂（４）
ａａｎ
ａｆ
ｏｒ／
式中，ｒ／为容积表面外法线方向。
基于有限体积法，将式（４）改写，如式（５）所示
ａ
．
ｆｒＤ（ＰｑｏｄＶ＋Ｊ矿ｄｉＶ（ｐ缈甜）ｄ矿＝
ｊｙｄｉｖ（Ｆｇｍｄ９）ｄＶ＋Ｊｒ＂Ｓ伊ｄＶ
（５）
式中，ｙ为润滑油的体积。
采用高斯散度公式，将体积分转换为面积分，
得到瞬态问题的表达式如式（６）所示Ｌ景（Ｌｐ两ｙ）出＋ＬＬ，７．（胛砧）以出＝
Ｊ＆Ｊ４，７・（厂ｇｒａｄ缈）以ｄｆ＋Ｊ出Ｊｒ品ｄｙ
（６）
万方数据
式中，Ａ为润滑油的面积，出为时间变化量。式（６）
是ＣＦＤ软件Ｆｌｕｅｎｔ进行瞬态计算的理论基础。
Ｆｌｕｅｎｔ软件求解压力．速度耦合方程的基本算法
是ＳＩＭＰＬＥ（Ｓｅｍｉ—ｉｍｐｌｉｃｉｔｍｅｔｈｏｄ
ｆｏｒｐｍｓｓｕｒｅｌｉｎｋｅｄ
ｅｑｕａｔｉｏｎｓ）算法。ＳＩＭＰＬＥ算法需要预定义流场压力。初始压力和参考压力的设置对计算结果影响比较大。尤其是动网格瞬态计算过程中出现回流现象时，预先建立流场压力分布对计算的收敛性至关重要。
对比工程实践，将进口边界条件设置为供油压力，参考压力为ｌＭＰａ，出油压力为０ＭＰａ，即泄油端与大气连通。
当采用动网格模型处理非定常压力．速度耦合
流动问题时，一般采用ＰＩＳＯ（Ｐｒｅｓｓｕｒｅ—ｉｍｐｌｉｃｉｔ
ｓｐｌｉｔｔｉｎｇｏｆ
ｏｐｅｒａｔｏｒｓ）算法，其计算基础是ＳＩＭＰＬＥ
算法。两者的区别在于预设压力求解动量方程和连
续方程后，ＰＩＳＯ算法需要连续两次修正压力和速度
方程，使压力和速度分别达到该时刻的两阶和三阶的精度。
鉴于轴颈高速对流场的影响，压力方程的插值
离散格式选择ＰＲＥＳＴＯ，动量方程的插值离散格式为ＱＵＩＣＫ。ＱＵＩＣＫ为流动的二次迎风插值格式，特别适于与流动方向对齐的结构网格，比二阶迎风
格式具有更精确的计算结果。其广义的表达式如（７）
所示
纯＝吉（击伽＋击啦）＋
ｉｌ丽％一百西细Ｊ
７
ｒ＆＋２＆
＆
、
（７）
…
式中，下标Ｐ为控制体积，下标Ｅ和形分别表示东
西方向与Ｐ相邻的控制体积；＆，＆，＆分别为
形，Ｐ，Ｅ相对应的控制体积的边长。伽，仇，伽
分别为Ｐ，Ｅ，形控制体积的通量，纯为Ｐ控制
体积西面的通量。
求解压力一速度耦合方程，得到油膜的三维物理
场。通过自定义程序循环指令遍历轴颈面域各节点的压力和黏性力并积分，得到不同方向的油膜力。当轴颈平衡位置变动过大且转速过高时，导致油膜破裂出现负压力，即油膜空化现象。计算油膜压力时，常用的处理方法是将负值区域压力置零。该处理方法和油膜力的计算所涉及的用户定义函数主要包括宏指令ＤＥＦＩＮＥＯＮＤＥＭＡＮＤ、循环函数
ｂｅｇｉｎ＿ｆ＿ｌｏｏｐ优ｆ）和面积函数Ｆ＿ＡＲＥＡ（Ａ，Ｚ矿１４１。
１．３弹簧光顺动网格模型
动网格模型主要应用于解决区域几何形状随时间变化的流体动力学问题。计算动静压轴承刚度阻尼系数时，网格变形主要考虑轴承间隙内油膜随轴
颈作刚体运动而产生的变化。在速度扰动和位移扰
２０１２年１２月熊万里等：基于动网格模型的液体动静压轴承刚度阻尼计算方法
１２１
动作用下，油膜与轴颈和油膜与轴瓦之间的边界层
度。至于大变形情况，则需要采用自定义程序解决
网格扭曲问题。
网格和油膜内部网格均产生形变。
由刚度阻尼计算模型可知，位移扰动和速度扰动产生的油膜力可线性化处理，位移扰动和速度扰
１．４油膜网格扭曲抑制
油膜网格扭曲形成的根本原因在于轴颈的旋转
运动和所采用的弹簧光顺更新模型。在确定油膜网格更新模型前提下，解决网格扭曲问题的落脚点在
动导致油膜网格产生的变形量亦可同样处理。油膜
网格运动的计算模型采用的是基于线性假设的弹性光顺Ｓｍｏｏｔｈｉｎｇ模型，其核心思想是胡克定律，网格移动与油膜力的关系【ｌ驯如式（８）所示
于轴颈的旋转运动。其关键就在于将轴颈旋转运动
参量代替实际转动嵌入至压力．速度耦合方程，计算瞬态油膜力。
Ｅ
２军如（Ｍ—Ａｄｙ）知２丽（８）
卫
．
１
解决网格扭曲问题时综合考虑了轴颈的运动。轴颈运动控制程序的基本功能包括：轴颈旋转运动和平动（位移扰动和速度扰动）。轴颈运动控制程序与动网格弹簧光顺模型相结合，实现轴颈边界网格和油膜网格的运动，进而计算得到轴颈受扰动条件
下不同位置的瞬态油膜力。
式中，Ｅ为节点ｆ的油膜力；反和ｄ，分别为节点ｉ
及相邻节点，的位置；Ａｄｉ和Ａｄｊ为节点位置变化
量；知为节点ｆ与节点．，之间的弹性系数：堵为与
节点ｉ相邻的节点数。
轴颈受外载荷作用作刚体平动。轴颈边界上的节点按胡克定律移动，通过流体内部节点传递到轴瓦边界上。在传递过程中，轴颈和轴瓦边界上网格
轴颈运动控制程序建立在Ｆｌｕｅｎｔ软件提供的用
户自定义函数基础上。轴颈的速度扰动和位移扰动由宏ＤＥＦＩＮＥ
ＣＧ
ＭＯＴＩＯＮ实现；为解决轴颈旋转
的形变对流体域网格的拉伸和挤压变形产生重要作用。为减小网格的畸变程度，需要对轴颈和轴瓦边界上网格节点的位置更新。节点控制如式（９）所示
引起的油膜网格畸变问题，将轴颈旋转的动边界转
变静边界，即将轴颈实际运动转化为物理量嵌入迭代计算，该功能通过宏ＤＥＦＩＮＥＰＲＯＦＩＬＥ实现ｕ
４ｌ。
钟州＝刃＋ＰＭ７
（９）
为验证轴颈运动控制程序成功加载动边界条件
实现轴颈对油膜网格的挤压运动，同时将轴颈旋转
式中，行和ｎ＋ｌ用于标识相邻迭代时间步的位置；∥为移动节点的松弛因子。该参数用于表征相邻内部节点的移动距离对变形边界节点移动幅度的影响程度。轴颈的运动将改变整个流体变形区域，影响流体物理场的变化。
在流体内部，网格运动时各个节点之间产生挤压或拉伸作用。节点之间的影响程度需结合物理现
的动边界条件成功转化为静边界条件，有效避免轴
颈旋转引起的油膜网格扭曲。对典型轴承利用Ｇａｍｂｉｔ建模，采用Ｆｌｕｅｎｔ进行动网格计算。
动网格计算对象为２排ｌＯ孔动静压无腔径向轴
承，如图ｌ所示。轴承基本结构参数：轴承孔径３０ｍｌＴｌ，长径比０．４，轴瓦厚１０ｍｉｌｌ，节流孔径０．６
ｍｍ，
象在ＰＩＳＯ算法中选择，默认设置对一般物理现象
通常都是有效的。流体内部节点各方向受力需取得平衡，相应的节点将产生移动，网格节点位移更新的计算模型如式（１０）所示
胁
节流孔沿轴向和周向均匀分布，轴承半径间隙２０ｐｍ。润滑油为２号主轴油，密度为８１０ｋｇ／ｍ３，动力黏度为０．００３
８５Ｐａ・Ｓ。
Ａｄ．ｍ＋１＝』Ｆ∑％
』
厶“日“Ｊ∑岛从，
（１０）
式中，ｍ为迭代次数。在轴颈旋转和平动作用下，油膜网格被动地产生形变。由于轴颈边界网格运动的特殊性（旋转和平动），弹簧光顺模型要求油膜网格类型具有非规则特征。对于规则网格，弹簧光顺模型需要特殊的指令才可以调用，但仅限于油膜网
ａ｝伦向轴承示意图
（ｂ）ｆ寺向轴嫒
图１径向轴承示意和零件图
格的挤压平动，轴颈旋转带动边界节点转动，致使
油膜网格形状扭曲，从而使计算发散。在小变形（挤压和扭曲）情况下，采用光顺模型具有较好的计算精万方数据
采用前处理软件Ｇａｍｂｉｔ建模。考虑到计算结果的精确性和网格数量等要求，流体区域采用六面体结构网格划分，使流体的流动迹线与网格分布的脉
机械工程学报第４８卷第２３期
络相适应，如图２所示。油膜厚度为１０层，油膜表面采用四边形划分，面网格采用Ｐａｖｅ模式，体网格采用Ｃｏｏｐｅｒ模式，流体域网格节点数量约为６０万。同时，设置相应的边界条件。节流孔定义为压力进
分布，偏布现象明显。这表明，在…定的偏心条件下，油压沿轴颈转动方向加强，且受轴颈转速影响明显，油膜压力场偏布表明轴颈转速已经成功加载。
口，轴承两端间隙为压力出口，轴颈和轴瓦为壁面。
置力ｐ／Ｐａ
１３０ＸｌＯ‘
ａ１流体域整体网络
｛ｈＩ局部放人示意图
（ａ）迭代１次
Ｊｔｉ力ｐ—Ｐ
蘸｛瓣
（ｂ）迭代３次
２’
ｈ－
图２滑动轴承网格模型
●３幽
１
采用Ｆｌｕｅｎｔ软件求解模型。求解器模型选择为非定常模式，压力．速度耦合方程的离散和计算模型
如前所述。边界条件设置：供油压力为３ＭＰａ，出油压力为０ＭＰａ。轴颈修改为运动壁面，其旋转速
度、扰动位移和扰动速度由轴颈运动控制程序实现，
轴瓦为静止壁面。轴承两侧出油面和间隙内流体均定义为变形区域。
轴颈的初始位置与轴承同心。调用自定义程序，对轴颈施加一个ｙ方向的速度扰动２ｌａｒｎ／ｓ，使轴颈向下运动。轴颈转速６ｋｒ／ｍｉｎ，单个迭代时间移动位移为２ｌａｍ。计算模型设置如前所述，初始化流场后，迭代６次，三维油膜压力分布和网格挤压变形
（ｃ）迭代６次
图３油膜网格受挤压变形及油膜压力分布
自定义的轴颈运动控制程序可实现轴颈的位移
扰动、速度扰动和旋转等功能，并有效避免油膜网
格的扭曲问题，从而为计算轴承油膜刚度和阻尼奠定了基础。然而，采用动网格方法计算面临一个重要问题：位移扰动和速度扰动均使轴颈产生位移，使油膜力产生变化。因此，界定位移扰动和速度扰动的取值范围就成为动网格计算油膜刚度和阻尼的又一关键问题。
１．５刚度阻尼的扰动范围界定
根据摄动理论，油膜刚度是指轴颈在平衡位置受位移扰动时承载力的变化，涉及轴颈转速和扰动位移两个物理量；油膜阻尼是指轴颈在平衡位置受速度扰动时承载力的变化，涉及轴颈转速和扰动速度两个物理量。利用雷诺方程计算时采用了虚位移
情况如图３所示（每个子图的上部表示油膜受轴颈
挤压后的压力分布，下部为受挤压后的油膜网格变形）。
轴颈平移和旋转运动均对油膜网格形变产生影响。由于采用了宏ＤＥＦＩＮＥＰＲＯＦＩＬＥ将轴颈旋转动边界转变为静边界，将运动方程嵌入动网格计算，避免了油膜网格的扭曲变形。油膜网格只在轴颈平
动方向产生挤压变形。对照不同迭代次数的网格变形情况，发现油膜网格只在轴颈平动方向产生挤压
变形，轴颈的旋转并未引起网格畸变。这表明将轴颈旋转动边界转换为静边界能避免轴颈旋转引起的网格畸变，同时说明宏ＤＥＦＩＮＥＣＧ—ＭＯＴＩＯＮ可控制轴颈的速度扰动和位移扰动。
原理，而动网格方法则反映出轴颈的真实运动。采用动网格方法计算刚度和阻尼的关键问题在于辨析
扰动位移与扰动速度对油膜力影响的内在联系。
在平衡位置，速度扰动和位移扰动对油膜承载
轴颈位移平动和旋转运动对油膜流场的影响可通过压力分布证实。轴颈向下移动和轴颈旋转均影
响油膜压力场的分布，不同之处在于轴颈旋转引起的动压效应使压力场偏布，位移平动使压力场对称分布。对照不同迭代时间的压力场可知，迭代１次时最高压力为３ＭＰａ：迭代６次时最高压力为
３．２
力的影响决定了计算刚度和阻尼时速度扰动和位移
扰动的取值范围。在瞬态计算和稳态计算的基础上，研究扰动位移与扰动速度对油膜力影响的差异。通
过对比研究动网格瞬态计算和稳态计算油膜力在平衡点的变化，分析扰动速度和扰动位移对油膜力的
影响，确定计算刚度阻尼时扰动位移和扰动速度的取值范围。
ＭＰａ，对比弦平面两侧最小间隙处油膜压力场万方数据
２０１２年１２月
熊万里等：基于动网格模型的液体动静压轴承刚度阻尼计算方法
对前面建立的计算模型，调用自定义程序使轴颈朝７方向运动，计算轴颈在不同扰动速度和扰动
位移下Ｘ和Ｙ方向的油膜力，并与稳态计算结果比较（扰动速度和扰动位移均为零）。
每个迭代时间步长，轴颈扰动位移为１ｇｍ，轴
颈在不同扰动速度ｖ作用下，ｘ和Ｙ两个方向分力与稳态计算结果对比如图４所示。
荟ｅＲ惫
涮懈簦
惹
偏心率ｓ
偏心率ｓ’
（ａ）ｘ方向油膜力的变化
（ｂ）ｙ方向油膜力的变化
图４扰动速度对ｘ和Ｙ方向油膜承载力的影响
由图４ａ可知，扰动速度对Ｘ方向分力的影响不
大。与稳态计算结果较靠近的速度为２．０ｇｒｎ／ｓ１和
２０．０
ｐＩＩｌ／ｓ１。扰动速度为２０
ｐ州ｓ１时，最大差距
１２．５％，平均差距７．７％；扰动速度为２．０ｇｒｒｄｓｌ时，平均差距７．５％。速度相差１０倍的情况下，平均差
距相差不大。
扰动速度对Ｙ方向分力的影响比较明显。由图
４ｂ可知，扰动速度为０．２ｍｍ／ｓ１时，最大差距为
２４．５％，平均差距１６％。扰动速度为２０．０ｎＩｌｌ／ｓ１时，
最大差距２．５％，平均差距１．８％；速度为２．０ｇｍ／ｓ１时，最大差距０．２％。
由上可知，在同样的扰动速度和扰动位移条件
下，采用动网格方法计算的Ｘ方向和Ｙ方向分力与稳态计算结果的差距不相同，挤压运动方向的差距大于与其垂直方向的差距，这与压力．速度耦合方程的求解模型有关。稳态计算和瞬态计算的基本润滑方程不同（与时间关联），求解模型不同，计算结果必然存在偏差。但从总体上看，在单个迭代步长内，扰动位移为１．０ｕｍ时，扰动速度不大于２０肛ｒｎ／ｓ１
的计算结果与稳态计算基本吻合。这个速度可视为
刚度和阻尼计算的分界点。
为减低扰动速度对油膜力的影响，单个迭代步长扰动速度取Ｏ．２“ＩＩｌ／ｓ１，在不同扰动位移Ｓ作用下，
Ｘ和Ｙ方向分力与稳态计算对比如图５所示。
对于ｘ方向分力的影响，扰动位移为２．０ｕｍ时，最大差距为２５．５％，平均差距为１９．５％；扰动位移为０．５ｇｍ时，最大差距１３．１％，平均差距７．８％，而位移为１．０“ｍ时平均差距为５．７％。
万方数据
冬
曩
茫
基
８
偏心率￡偏心率占
（ａ）ｘ方向油膜力的变化
（ｂ）ｙ方向油膜力的变化
图５位移扰动对Ｘ和Ｙ方向承载力的影响
对于Ｙ方向分力的影响，位移扰动为２．０“ｍ时，最大差距为１１％，平均差距为６．４％；扰动位移为０．５肛ｍ时，最大差距１．３％，平均差距０．６％，而位移为１．０ｇｍ时最高差距为０．３％。
对于扰动位移对Ｘ和Ｙ方向油膜力的影响与扰
动速度不同，Ｙ方向的差距小于ｘ方向的差距。由图５可知，位移为１．０ｇｍ与０．５９ｍ时，油膜力较为接近。Ｘ方向的差距为２．１％，Ｙ方向的差距为０．３％。
这种情况可视为刚度计算。
由上可知，扰动速度诞２０ｇｍ／ｓ１，扰动位移ｓ＜ｌ
ｇｍ，扰动速度和扰动位移对油膜力的影响较小，
平均差距不大于８％。动网格的瞬态计算与稳态
计算基本相当，可视为刚度计算。扰动速度
ｖ＞０．２
ｍｍ／ｓ１时，速度对油膜力的影响比较明显。动
网格的瞬态计算不能近似为稳态计算，可视为阻尼计算。在合适的速度和位移取值范围内，采用动网格计算轴承刚度和阻尼系数是切实可行的。
２典型轴承的动网格计算
采用动网格方法计算轴承的刚度阻尼系数包括
两个部分：一是计算主轴在轴承间隙内的平衡位置；
二是在平衡位置处计算轴承的刚度阻尼。
轴承结构如图１所示，轴承设计参数和润滑油
如前述。主轴的质量为５．５４蚝。建模时取轴承与轴
颈的同心位置为初始位置。重力作用下的最大移动速度为２．０ｇｍ／ｓ１，最大位移为１．０ｐｍ。计算主轴平衡位置时理论上应满足油膜力水平分力为０，垂直分力等于５５．４Ｎ。考虑到计算效率，误差为０．００１
Ｎ
时停止迭代。主轴的偏位角和偏心距随转速的变化
如图６所示。
计算刚度阻尼时需偏离平衡位置一定距离以适
应轴颈的运动。采用Ｇａｍｂｉｔ建模时需要将轴颈往ｘ或Ｙ方向偏移１．０ｇｍ。计算轴承刚度时，扰动速度为２．０ｇｍ／ｓ１，扰动位移１．０“ｍ；计算阻尼时，扰动速度０．５ｍｍ／ｓ１，扰动位移１．０ｇｍ，分别计算出轴承
机械工程学报第４８卷第２３期
ｆ
￥
吣
３
提翅ｇ
２
￡Ｔ≈魁尊弯
Ｏ
０
转速ｎ／（ｋｒ／ｍｉｎ）转速ｎ／（ｋｒ／ｍｉｎ）（ａ）转速对偏位角的影响
（ｂ）转速对偏心距的影响
图６转速对主轴平衡位置的影响
在速度扰动和位移扰动条件下的承载力，并利用式
（２）计算出轴承的刚度和阻尼系数，刚度阻尼系数随
转速的变化如图７、８所示。
堇１
要
＋＾ｗ—●—＾、．ｘ
｛
童
嗯。．’—日；．＋女Ｊｙ
弓
≥１０
墅
７＼√７’
蔷５
螽
曼
ｎ＇
７
Ⅻ
（ａ）转速对直接刚度的影响【ｂ）转速对交叉刚度的影响
图７转速对油膜刚度的影响
§０
（ａ）转速对直接阻尼的影响
馨
（ｂ）转速对交叉阻尼的影响
图８转速对油膜阻尼的影响
转速对油膜刚度和阻尼的影响各不相同。随着转速由１．５增至１２．０ｋｒ／ｍｉｎ，直接刚度略有下降，但交叉刚度迅速增加，由１．６增至１２．０Ｎｌｌｍ＿ｌ；直接阻尼呈下降趋势，交叉阻尼在０Ｎ／１．ｔｍ附近波动。采用动网格方法计算得到刚度阻尼随转速的变化趋势，与采用雷诺方程计算的变化趋势基本相同ｐＪ。
３试验分析
３．１试验方案
为验证采用动网格方法计算刚度和阻尼系数的可靠性，对典型轴承的刚度系数进行试验测试（刚度和阻尼由同一种方法计算得到，仅对刚度进行试验测试）。测试平台包括：双排ｌＯ孔液体动静压轴承电主轴，ＪＣ．３０００型油箱冷却器，ＶＦＤ．Ｂ型变频器，北京派莱博科技有限公司开发的振动及动态信号采集分析系统。
试验原理如图９所示。在主轴的正交平面０和
Ｙ方向，坐标与轴承一致）布置两个电容传感器，将万方数据
主轴对传感器的位移变化直接转换为电压信号，通
过电容测微仪和数据采集卡转换得到位移的变化。
主轴的转速由变频器控制。图１０为测量装置。
变频器／电器柜
三相电源几：；左丌＼９量ｌ，卜——强地７鼍兰兰｛；当
仁、ｆ≤；刁
手
液压站／冷却泵
嚣跫盟＾岸峙垤测头。人粤ｎ厂—丝叶仙
—————＜，
瘦体悬浮电生轴ｌ
、、
、匠！ｒ
——／／数据采集卡试验原理图
｛铨麓蝥０
Ｉｋ如＋七坩缸＝
Ⅲ，”…
｛嚣麓嚣
【‰勾＋‰勘２６
㈣，
”…
仁粪薹套㈣卜羔ｋ２蕊¨叫
２０１２年１２月
熊万里等：基于动网格模型的液体动静压轴承刚度阻尼计算方法
１２５
个速度点空运行半小时后，在１和＿ｙ方向先后施加
２０
Ｎ和４０Ｎ作用力，每个作用力加载３次，记录
每次加载力时ｘ和Ｙ方向的位移偏移量，取平均值作为该速度点和作用力的有效位移。试验测试出４８
组数据，将结果代入式（１３），计算出轴承的刚度系数。
３．３理论与试验对比分析
将动网格数值计算结果（图７）与试验测试结果进行比对，刚度随转速的变化趋势如图１１所示。
畲１
毒
言１
Ｚ
÷
￥
邑
曹
≥
划
刨一：麓彳—．＿Ｆ＝ＯＮ
二，ｒｔ—，，芑
—卜Ｐ＿２０Ｎ
毫
萎辎
辎柳
０
枷
转速ｎ／（ｋｒ／ｍｉｎ）转速ｎ／（ｋｒｌｍｉｎ）
（ａ）直接刚度￡“
（ｂ）直接刚度锄
言
毒Ｚ
￥
毒
魁一
蔑
以
权一
转速ｎ／（ｋｒ／ｍｉｎ）转速ｎ／（ｋｒｌｍｉｎ）（ｃ）交叉刚度≈矽（ｄ）交叉刚度盘ｐ
图１１
动网格计算结果和试验测试结果对比
通过对比动网格数值计算Ｆ＝ＯＮ与试验测试
Ｆ＝２０
Ｎ和Ｆ＝４０Ｎ平均值可知，两者偏差不大。直
接刚度ｋ最大偏差为１０．４％，平均偏差６．１６％；直接刚度‰的最大偏差２０．９％，平均偏差１３．３％；交
叉刚度‰最大偏差１５％，平均偏差９．８％；交叉刚
度ｋ与‰相差不大。
对于局部各数据点，动网格计算与试验测试结果并非完全一致，这与各自的外部因素有关。动网格计算偏向理论分析，其结果与求解器模型、网格
质量、边界条件和动网格模型以及参数设置有关；
试验测试侧重于工程实践，其结果受转子质量动不平衡、轴承加工误差、轴系装配误差、转速的脉动
和供油的平稳性以及加载力的平稳性和方向的平行
度等因素影响，以至于Ｆ＝２０Ｎ和Ｆ＝４０Ｎ两组数据
也存在一定差距。对于总体变化趋势，动网格计算结果与试验测试数据基本吻合，平均误差不超过
１２．５％，体现了典型轴承动力特性系数的内在一致性。
４结论
（１）采用自定义程序实现轴颈位移扰动、速度扰动、旋转和油膜力计算等功能，并利用弹簧光顺万方数据
模型更新网格的挤压变形，计算了轴承的刚度和阻尼系数。通过试验测试轴承的刚度系数，发现两者
基本吻合，表明采用动网格模型计算动静压轴承刚
度阻尼的数值方法是可行和有效的，从而为动静压
轴承动力特性系数计算提供了一种新方法。
（２）将轴颈的旋转动边界条件转换为静边界条件，并将转速加载在静边界上，有效避免了因轴颈旋转引起的油膜网格扭曲畸变，为旋转机械流体的动边界条件加载提供了一种新途径。
（３）对比研究了扰动速度和扰动位移对油膜力
的影响规律，界定了计算刚度阻尼系数时扰动速度
和扰动位移的合理取值范围，解决了动网格技术计
算轴承刚度阻尼系数的关键问题。
（４）再现了轴颈受位移扰动和速度扰动情况下油膜厚度的动态发展过程，反映了三维流场压力分布的时变特性，使计算过程更符合工程实际。
参考文献
［１】熊万里，林厚波，吕浪，等．超高速磨削电机内装式液
体悬浮电主轴关键技术及研发［ＥＢ／ＯＬ］．【２０１Ｉ－１２—２９］．
ｈｔｔｐ：／／ｗｗｗ．ｈｙｄｒｏ－ｓｐｉｎｄｌｅ．ｃｏｍ．
ＸＩＯＮＧ
Ｗａｎｌｉ，ＬＩＮＨｏｕｂｏ，ＬＯＬａｎｇ，ｅｔ
ａ１．Ｒｅｓｅａｒｃｈａｎｄ
ｄｅｖｅｌｏｐｍｅｎｔ
ｏｎ
ｋｅｙｔｅｃｈｎｏｌｏｇｙｏｆｍｏｔｏｒｂｕｉｌｔ・ｉｎｓｐｉｎｄｌｅ
ｓｕｐｐｏｒｔｅｄｂｙｈｙｄｒｏｄ’７ｎａｍｉｃ
ａｎｄ
ｈｙｄｒｏｓｔａｔｉｃｂｅａｒｉｎｇｓｆｏｒ
ｕｌｔｒａ—ｈｉｇｈｓｐｅｅｄｇｆｉｎｄｉｎｇ［ＥＢ／ＯＬ］．［２０１１－１２—２９］．ｈｔｔｐ：／／ｗｗｗ．ｈｙｄｒｏ－ｓｐｉｎｄｌｅ．ｃｏｍ．
［２】ＬＵＮＤ
Ｊ
Ｗ．Ｔｈｅ
ｓｔａｂｉｌｉｔｙｏｆ
ａｌｌ
ｅｌａｓｔｉｃｒｏｔｏｒｉｎ
ｊｏｕｒｎａｌ
ｂｅａｒｉｎｇｗｉｔｈｆｌｅｘｉｂｌｅｄａｍｐｅｄｓｕｐｐｏｒｔｓ［Ｊ］．ＡＳＭＥ
Ｊｏｕｒｎａｌ
ｏｆＡｐｐｌｉｅｄＭｅｃｈａｎｉｃｓ，１９６５，３２（４）：９１１－９２０．
［３】ＧＨＯＳＨＭＫ，ＭＡＪＵＭＤＡＲＢＣ，ＲＡＯＪＳ．Ｓｔｅａｄｙ・ｓｔａｔｅ
ａｎｄｄｙｎａｍｉｃｂｅｈａｖｉｏｕｒｏｆｍｕｌｔｉ－ｒｅｃｅｓｓｈｙｂｒｉｄｏｉｌｊｏｕｒｎａｌ
ｂｅａｒｉｎｇ［Ｊ］．Ｊｏｕｒｎａｌ
ｏｆＭｅｃｈａｎｉｃａｌＥｎｇｉｎｅｅｒｉｎｇＳｃｉｅｎｃｅ，
１９７９，２１（５）：３４５－３５１．
［４】ＲＯＷＥＷ
Ｂ，ＣＨＯＮＧＦＳ．Ｃｏｍｐｕｔａｔｉｏｎｏｆｄｙｎａｍｉｃ
ｆｏｒｃｅｃｏｅｆｆｉｃｉｅｎｔｓｆｏｒ
ｈｙｂｒｉｄ（ｈｙｄｒｏｓｔａｔｉｃ／ｈｙｄｒｏｄｙｎａｍｉｃ）
ｊｏｕｒｎａｌｂｅａｒｉｎｇｓｂｙｔｈｅｆｉｎｉｔｅ
ｄｉｓｔｕｒｂａｎｃｅ
ａｎｄｐｅｒｔｕｒｂａｔｉｏｎ
ｔｅｃｈｎｉｑｕｅｓ［Ｊ］．ＴｒｉｂｏｌｏｇｙＩｎｔｅｒｎａｔｉｏｎａｌ，１９８６，１９（５１：
２６０．２７１．
［５】许尚贤，陆永颐，万鹏．小孔式动静压轴承的特性分析
与实验研究［Ｊ】．机械工程学报，１９９３，２９（５）：８８—９５．
ＸＵ
Ｓｈａｎｇｘｉａｎ，ＬＵＹｏｎｇｙｉ，ＷＡＮ
Ｐｅｎｇ．Ａｎａｌｙｓｉｓａｎｄ
ｅｘｐｅｒｉｍｅｎｔａｌｉｎｖｅｓｔｉｇａｔｉｏｎ
ｏｆ
ａ
ｈｏｌｅ－ｅｎｔｒｙ
ｈｙｂｒｉｄ
ｂｅａｒｉｎｇ［Ｊ］．Ｃｈｉｎｅｓｅ
ＪｏｕｒｎａｌｏｆＭｅｃｈａｎｉｃａｌＥｎｇｉｎｅｅｒｉｎｇ，
１９９３，２９（５）：８８－９５．
［６］６郑铁生，许庆余．滑动轴承动特性系数新算法【Ｊ】．机械
１２６
机械工程学报
第４８卷第２３期
工程学报，１９９６，３２（３）：２０．２７．
ＺＨＥＮＧＴｉｅｓｈｅｎｇ，ＸＵＱｉｎｇｙｕ，Ｎｅｗ
ｎｕｍｅｒｉｃａｌｍｅｔｈｏｄｆｏｒ
ｃａｌｃｕｌａｔｉｎｇ
ｄｙｎａｍｉｃ
ｃｏｅｆｆｉｃｉｅｎｔｓ
ｏｆｏｉｌ－ｌｕｂｒｉｃａｔｅｄ
ｊｏｕｒｎａｌｂｅａｒｉｎｇｓ［Ｊ】．Ｃｈｉｎｅｓｅ
Ｊｏｕｒｎａｌｏｆ
Ｍｅｃｈａｎｉｃａｌ
Ｅｎｇｉｎｅｅｒｉｎｇ，１９９６，３２（３）：２０－２７．
［７】ＳＡＷＩＣＫＩＪＴ，ＲＡＯＴＶＶＬＮ．Ａｎｏｎｌｉｎｅａｒｍｏｄｅｌｆｏｒ
ｐｒｅｄｉｃｔｉｏｎ
ｏｆｄｙｎａｍｉｃ
ｃｏｅｆｆｉｃｉｅｎｔｓ
ｉｎ
ａ
ｈｙｄｒｏｄ）‘ｎａｍｉｃｊｏｕｒｎａｌｂｅａｒｉｎｇ［Ｊ］．Ｉｎｔｅｒｎａｔｉｏｎａｌ
Ｊｏｕｒｎａｌ
ｏｆ
Ｒｏｔａｔｉｎｇ
Ｍａｃｈｉｎｅｒｙ，２００４，１０（６）：５０７－５１３．［８］ＭＯＲＩＡ，ＭＡＫＩＮＯＴ，ＭＯＲＩＨ．Ｅｎｔｒｙ
ｆｌｏｗａｎｄｐｒｅｓｓｕｒｅ
ｊｕｍｐ
ｉｎ
ｓｕｂｍｅｒｇｅｄｍｕｌｔｉ－ｐａｄ
ｂｅａｒｉｎｇｓ
ａｎｄ
ｇｒｏｏｖｅｄ
ｂｅａｒｉｎｇｓ［Ｊ】．ＪｏｕｒｎａｌｏｆＴｒｉｂｏｌｏｇｙ，１９９２，１１４（２）：３７０－３７７．［９］ＴＵＣＫＥＲ
Ｐ
Ｇ，ＫＥＯＧＨ
ＰＳ．Ｏｎｔｈｅｄｙｎａｍｉｃｔｈｅｒｍａｌ
ｓｔａｔｅｉｎ
ａ
ｈｙｄｒｏｄｙｎａｍｉｃｂｅａｒｉｎｇｗｉｔｈ
ａ
ｗｈｉｒｌｉｎｇ
ｊｏｕｍａｌ
ｕｓｉｎｇＣＦＤｔｅｃｔｍｉｑｕｅｓ［Ｊ】．ＪｏｕｒｎａｌｏｆＴｒｉｂｏｌｏｇｙ，１９９６，
１１８（２）：３５６－３６３．［１０］ＳＡＨＬＩＮ
Ｆ，ＧＬＡＶＡＴＳＫＩＨＳ
Ｂ，ＡＬＭＱＵＩＳ
Ｔ，ｅｔａ１．Ｔｗｏ
ｄｉｍｅｎｓｉｏｎａｌＣＦＤａｎａｌｙｓｉｓｏｆｍｉｃｒｏ．ｐａＲｅｍｅｄｓｕｒｆａｃｅｓｉｎ
ｈｙｄｒｏｄｙｎａｍｉｃｌｕｂｒｉｃａｔｉｏｎ［Ｊ］．ＪｏｕｒｎａｌｏｆＴｒｉｂｏｌｏｇｙ，２００５，
１２７（１）：９６—１２７．
万方数据
ＧＥＲＴＺＯＳ
Ｋ
Ｐ，
ＮＩＫＯＬＡＫＯＰＯＵＬＯＳ
Ｐ
Ｇ，ＰＡＰＡＤＯＰＯＵＬＯＳＣＡ．ＣＦＤａｎａｌｙｓｉｓｏｆｊｏｕｒｎａｌ
ｂｅａｒｉｎｇ
ｈｙｄｒｏｄｙｎａｍｉｃ
ｌｕｂｒｉｃａｔｉｏｎ
ｂｙＢｉｎｇｈａｍｌｕｂｒｉｃａｎｔ［Ｊ］．
Ｔｒｉｂｏｌｏｇｙ
Ｉｎｔｅｒｎａｔｉｏｎａｌ，２００８，４１（１２）：１
１９０－１２０４．
［１２】王福军．计算流体动力分析－ＣＦＤ软件原理与应用［Ｍ】．
北京：清华大学出版社，２００４．ＷＡＮＧＦｕｊｕｎ．Ｃｏｍｐｕｔａｔｉｏｎａｌ
ｆｌｕｉｄ
ｄｙｎａｍｉｃｓ
ａｎａｌｙｓｉｓ－ｔｈｅ
ｐｒｉｎｃｉｐｌｅｓａｎｄａｐｐｌｉｃａｔｉｏｎｓｏｆＣＦＤｓｏｔｔｗａｒｅ［Ｍ］．Ｂｅｉｊｉｎｇ：
ＴｓｉｎｇｈｕａＵｎｉｖｅｒｓｉｔｙ
Ｐｒｅｓｓ，２００４．
［１３］
ＦＬＵＥＮＴ．ＦＬＵＥＮＴ６．３
ｕｓｅｒ’ｓ
ｇｕｉｄｅ［Ｍ］．Ｌｅｂａｎｏｎ：Ｎｅｗ
Ｈａｍｐｓｈｉｒｅ（ＵＳＡ）Ｆｌｕｅｎｔ
Ｉｎｃ．，２００６．
［１４］
ＦＬＵＥＮＴ．ＦＬＵＥＮＴ６．３
ＵＤＦ
ｍａｎｕａｌ［Ｍ］．Ｌｅｂａｎｏｎ：Ｎｅｗ
Ｈａｍｐｓｈｉｒｅ（ＵＳＡ）ＦｌｕｅｎｔＩｎｃ．，２００６．
［１５］
张直明．滑动轴承的流体动力润滑理论［Ｍ］．北京：高等教育出版社，１９８６．
ＺＨＡＮＧ
Ｚｈｉｍｉｎｇ．Ｌｕｂｒｉｃａｔｉｏｎｔｈｅｏｒｙ
ｆｏｒ
ｓｌｉｄｅ
ｂｅａｒｉｎｇ［Ｍ］．Ｂｅｉｊｉｎｇ：Ｈｉｇｈｅｒ
ＥｄｕｃａｔｉｏｎＰｒｅｓｓ，１９８６．
作者简介：熊万里（通信作者），男，１９７１年出生，博士，教授，博士研
究生导师。主要研究方向为高速电主轴系统动力学及超高速磨削装备。
Ｅ－ｍａｉｌ：ｗａｎ３６９＠ｖｉｐ．ｓｉｎａ．００ｍ
基于动网格模型的液体动静压轴承刚度阻尼计算方法
作者：作者单位：
熊万里， 侯志泉， 吕浪， 阳雪兵， 袁巨龙， XIONG Wanli， HOU Zhiquan， L(U) Lang， YANG Xuebing ， YUAN Julong
熊万里,侯志泉,吕浪,XIONG Wanli,HOU Zhiquan,L(U) Lang(湖南大学国家高效磨削工程技术研究中心 长沙410082)， 阳雪兵,YANG Xuebing(湘潭电机股份有限公司风能分公司 湘潭411102)， 袁巨龙,YUAN Julong(浙江工业大学机械制造及自动化教育部重点实验室 杭州310014)机械工程学报
JOURNAL OF MECHANICAL ENGINEERING2012,48(23)
刊名：英文刊名：年，卷(期)：
1. 熊万里;林厚波;吕浪 超高速磨削电机内装式液体悬浮电主轴关键技术及研发 2011
2. LUND J W The stability of an elastic rotor in journal bearing with flexible damped supports 1965(04)3. GHOSH M K;MAJUMDAR B C;RAO J S Steady-state and dynamic behaviour of multi-recess hybrid oil journal beating 1979(05)
4. ROWE W B;CHONG F S Computation of dynamic force coefficients for hybrid (hydrostatic/hydrodynamic)journalbearings by the finite disturbance and perturbation techniques 1986(05)5. 许尚贤;陆永颐;万鹏 小孔式动静压轴承的特性分析与实验研究 1993(05)6. 郑铁生;许庆余 滑动轴承动特性系数新算法 1996(03)
7. SAWICKIJT;RAOTVVLN A nonlinear model for prediction of dynamic coefficients in a hydrodynamic journalbearing [外文期刊] 2004(06)
8. MORIA;MAKINO T;MORI H Entry flow and pressure jump in submerged multi-pad bearings and grooved bearings1992(02)
9. TUCKER P G;KEOGH P S On the dynamic thermal state in a hydrodynamic bearing with a whirling journal usingCFD techniques[外文期刊] 1996(02)
10. SAHLIN F;GLAVATSKIH S B;ALMQUIS T Two dimensional CFD analysis of micro-patterned surfaces in hydrodynamiclubrication 2005(01)
11. GERTZOS K P;NIKOLAKOPOULOS P G;PAPADOPOULOS C A CFD analysis of journal beating hydrodynamic lubrication byBingham lubricant[外文期刊] 2008(12)
12. 王福军 计算流体动力分析-CFD软件原理与应用 200413. FLUENT FLUENT 6.3 user's guide 200614. FLUENT FLUENT 6.3 UDF manual 200615. 张直明 滑动轴承的流体动力润滑理论 1986
本文链接：http://d.g.wanfangdata.com.cn/Periodical_jxgcxb201223018.aspx