于贺春，岳金珂，王 进，张国庆，王文博，王仁宗

（中原工学院 机电学院，郑州 450007）

1 研究进展

高精度的加工制造对转子的回转精度有着极高的要求，空气静压轴承作为高速电主轴的主要支撑部件，具有回转精度高、低摩擦、发热小等优点［1］。转子受空气轴承的支撑作用，转子通过偏心为系统提供承载力［2］。转子在工作过程中产生的微小振动会改变轴承的承载特性，影响零件的表面加工质量［3］。因此，精确的轴承特性分析需结合转子的运动状态。

随着气体润滑技术的进步，为适应超精密技术的发展，对空气静压轴承特性的研究不断深入。轴承动态特性一般利用动态刚度与阻尼系数进行衡量，通过对气体流动状态与转子运动轨迹进行不同假设，利用不同的计算方式，得到不同状态下轴承动刚度和阻尼系数。戚社苗等［4］通过采用偏导数法求解雷诺方程，推导出动压气体轴承动态刚度和动态阻尼系数普遍适用的计算方法，从理论上解释气体轴承的动态刚度和动态阻尼系数与轴颈扰动频率的相关性。于贺春［5］等利用旋转坐标系假定转子涡动轨迹，与静态网格相结合，求出轴承在特定转子涡动轨迹下的刚度阻尼系数。陈冬菊等［6］引入稀薄效应下的流量因子，利用差分法求解流体方程，将得出的稳态结果导入扰动方程，得到的轴承动态刚度和动态阻尼系数。BOUZIDANE等［7］通过对径向轴承内流体流动进行不同形式的假设，采用数值方法对系统进行线性建模，研究了油膜厚度、凹槽压力和几何构型对静压轴颈轴承等效刚度和阻尼的影响。WANG［8］利用线性摄动理论和迭代过程解决非稳态雷诺方程，通过数值计算系统各个参数下的动刚度和动阻尼，研究了轴颈旋转和轴承表面波纹度对轴颈轴承空气静力学动态性能的影响。SHI等［9］通过联立稳态雷诺方程和多个耦合摄动雷诺方程获得轴承的动态刚度与阻尼系数，证明了气膜厚度与压力之间的非线性关系。CUI等［10］通过研究空气静压轴承角刚度对加工表面的影响，建立偏转数值模型，通过坐标系下的流场转换，结合有限元法和比例划分法，提出了轴承角刚度的计算方法，结果证明：可通过增加非平面度误差和推力偏心率的幅度来提高推力轴承的角刚度。因雷诺方程不能反映高速流场周向惯性效应、轴向动压效应与扩散效应等对三维流场的影响。为提高计算精度，GERTZOS等［11］引入“可动边界”，采用动网格技术计算了滑动轴承在动载作用下轴承的压力分布。熊万里等［12］提出基于N－S方程的动网格计算轴承的刚度阻尼新方法，采用程序编程实现轴颈的旋转速度、位移扰动和速度扰动以及扰动过程中油膜力的计算等，对比瞬态和稳态油膜力的变化，得到了在计算刚度阻尼系数时的速度扰动与位移扰动取值范围。YU等［13］利用Fluent软件进行数值模拟，对轴承增加轴向正弦摄动力，结合动网格方法，通过仿真与实验结果的对比，结果证明：轴向刚度系数随扰动频率的增加而增加，而轴向阻尼系数则随扰动频率的增加而减小。林禄生［14］通过建立轴心轨迹方程，利用动网格技术对滑动轴承进行非线性瞬态计算，得到动载荷下的系统轴心轨迹，通过动力学分析证明了流固耦合方法可以准确分析复杂转子系统的非线性轴心轨迹。崔海龙等［15］建立气体轴承的双向流固耦合数值模拟模型，通过静态数值模拟的方法的得到不同设计参数对承载力刚度的影响规律，随着动态气膜间隙逐渐增大，气膜承载力逐渐减小。毛文贵等［16］对单圆盘弹性转子系统进行分析，分别利用Hahn法和传递矩阵法计算得到了轴承轴颈处和转子部分的轴心轨迹，并采用高效微型遗传算法对系统的轴心轨迹进行优化。李威建［17］立倾斜状态下的轴承－转子系统模型，对其进行流固耦合分析，分析了阶跃载荷和斜坡载荷运动轨迹变化，并通过非接触式加载方式进行实验，对比分析系统在流固耦合计算下的准确性。唐璐阳等［18］利用动网格技术模拟流场变化，对轴承施加偏心率，求出径向轴承和止推轴承刚度的变化情况，对不同供气压力和不同偏心下的动态特性进行分析。胡辰等［19］通过Fluent动网格技术对径推—体式空气轴承进行数值仿真，使用程序控制壁面旋转与网格节点的运动，对不同状态下的轴承刚度变化进行了分析讨论。

综上所述，当前研究中多采用线性方法或摄动法，没有充分考虑转子运动及气膜力的非线性变化。在利用流固耦合方法分析时，对转子的非线性分析不足。因此在充分考虑转子运动与气膜力变化的相互影响的情况下，通过建立径向小孔轴承－转子系统动态耦合模型，将流体力学与固体力学进行耦合求解，利用时域响应曲线对转子的非线性现象进行描述，通过大量的数值计算，利用瞬态响应参数与稳态响应参数表征系统在受到阶跃负载状态后的非线性运动状态的变化，分析一定转速范围内负载突变后，系统整个过程变化趋势及特征。

2 耦合模型的建立

2.1 流体－固体耦合模型

在同一坐标系下，利用Gambit软件和SOLIDWORKS软件分别建立完整的流场模型和转子模型。初始状态下，转子轴心与轴承轴心重合，气体通过进气孔进入轴承间隙，转子与气体交界处形成动态耦合面，系统基本模型参数见表1。

表1 系统基本模型参数

如图1为小孔节流径向轴承结构，24个节流孔沿圆周方向均匀分布为两排。轴承结构参数：轴承的内径D1＝50 mm，轴承的长度L1＝60 mm，节流孔距最近端面L2＝15 mm，转子模型轴颈D2＝49.96 mm，转子总长度L3＝100 mm，气膜的平均厚度h＝20μm。

图1 系统耦合模型轴承结构示意图

2.2 网格划分

图2为流体与固体网格划分：在满足网格最大纵横比的条件下，利用Gambit软件对气膜进行网格划分，为保证网格的质量将流场区域分为3个部分：节流孔区、气腔区、气膜间隙区，如图2（a）、（b）所示。网格随流场结构的改变由密到疏逐渐变化，对结构尺寸变化较大的过度区域进行网格加密，在气膜间隙的径向方向对网格进行边界层的划分，在保证计算精度与计算效率的条件下对网格层数进行调整。

图2 网格划分示意图

系统的转子图2（c）采用尺寸相近的六面体网格进行划分，保证网格在变形过程中质量的同时，满足每一步迭代过程之间数据传递的精确度。

2.3 边界条件的设置

1）流经系统的气体为理想气体，过程为绝热过程，温度恒定；

2）气体的流动过程为层流流动；

3）轴承表面以及转子旋转表面绝对光滑；

轴承的初始进气压力P1＝0.4 MPa，出口压力P2＝0 MPa。转子在旋转过程中不发生轴向窜动。在转子的轴颈中心处设置监测点，通过记录监测点的运动轨迹得到轴颈处的运动轨迹。

3 仿真技术

3.1 流固耦合技术

流固耦合技术将流体力学与固体力学结合在一起的分析方法。依据数据传递方式分为单向传递和双向传递。单向流固耦合只考虑流体区对固体状态的改变；双向流固耦合同时考虑2个物理区域之间的相互影响，固体区在流体作用下状态发生改变，固体状态的改变反过来引起流场的变化，2者互相影响最终达到1个动态平衡，是1个动态变化过程。图3为耦合求解流程框图。

图3 耦合求解流程框图

系统进行耦合分析时，满足质量守恒、动量守恒、能量守恒，联立流体基本方程与转子动力学方程，得到轴承－转子系统交界面处的耦合控制方程，即气膜与转子在流固交界面处的应力、位移以及温度等都相等。

式中：τ为耦合交界面处的应力；r为耦合交界面处的位移；q为耦合交界面处的热流量；T为耦合交界面处的温度；f为流体；s为固体。

根据边界条件的假设，气体通过时未与周围环境进行热交换，系统温度保持恒定。

3.2 动网格技术

在系统的耦合计算过程中，转子转动过程中的位置移动会改变气膜形状和流动状态。利用CFD动网格技术对内部流场进行模拟，网格节点随着流场改变及转子运动发生变化，通过动网格方法不断对畸形网格进行光顺、重构。动网格方法保证了网格质量，以及数据在传递过程中的准确性。

4 不同状态下系统特性的耦合计算

4.1 阶跃负载下系统特性的耦合计算

当系统的进气压力P1＝0.4 MPa，转子转速n＝0 r／min时，转子轴颈处突然受到大小为400 N阶跃载荷（负载方向为Y的负方向），对该系统的稳定过程进行的耦合计算分析。

图4为突变载荷下转子运动轴心轨迹。图像表明：转子在负载和气膜力的作用下发生涡动，且涡动的轨迹由大到小，最终稳定在一定的范围内。

图4 转子运动轴心轨迹

在系统时间响应分析过程中，转子在Y方向的位移波动规律符合2阶欠阻尼系统的时间响应曲线，通过利用时间响应参数对Y方向的非线性运动进行时间响应分析。部分时间响应分析参数见表2。

表2 时间响应分析部分参数符号

在瞬态响应过程中，系统受到突变载荷发生震荡，阻尼比ξ决定了系统的震荡特性：ξ越大，系统的相对稳定性越好；转子位移由零时刻起，首次达到稳态位移值所需要的时间（上升时间tr），转子从零时刻移动至第1次幅值所需要的时间（峰值时间tp），tr、tp的大小表示系统受到突变载荷后的反应时间；最大峰值与稳态值之差与稳态值的比值（最大超调量Mp）满足：

阻尼比的大小决定最大超调量的值，ξ值越大，最大超调量越小，系统的平稳性越好。

稳态响应过程中，转子轨迹稳定在一定范围内进行波动，利用相同时间段内的位移差（波动平均差值Δx）表征系统稳态响应期间的准确性（回转精度）。

图5为转子在X、Y方向位移曲线。在瞬态响应阶段：当t＞0时，转子在负载和自身重力的作用下，向Y的负方向移动，位移量逐渐变大。随着时间的增加，当t＝tp1时，转子位移量达到最大，随后转子向正方向移动，t＝tp2时，再次向负方向移动，且波动峰值越来越小；在稳态响应阶段，转子位移稳定在－8.5μm附近。在整个过程中，转子在X方向上不受负载直接作用，最终稳定在0.3μm位置。

图5 转子在X、Y方向位移曲线

图6为不同时刻下的流场变化云图。随着时间的增加，当t＜0.5 ms时，转子在合力作用下移动，使气膜间隙变小，Y的负方向气膜被压缩，气体的高压区面积变大，t＝0.5 ms时气膜高压区面积达到波动最大值，当0.5 ms＜t＜1.6 ms时，转子所受合力方向为负载的负方向，向上运动，气膜高压区面积变小，t＝1.6 ms气膜高压区面积达到波动最小值；当t＞2 ms时，转子位移稳定在一定范围内，气膜的高压区面积保持稳定。

图6 不同时刻下的流场变化云图

图7为Y方向气膜承载力与位移随时间变化曲线。随着转子在Y方向位移的不断变化，位移越大，气体高压区面积增大，气膜承载力越大；位移变小，气膜高压区面积减小，气膜承载力越小。当t＝tp1时，转子位移最大为－12μm，对应最大气膜承载力为556 N，气膜承载力大于负载，转子位移减小，当转子位移在－8μm时，对应的气膜力为389 N，负载大于气膜承载力，转子位移增大；最终转子趋于稳定，承载力在一定范围内保持稳定。

图7 系统瞬态响应时气膜力与位移关系

图8表示系统在稳态响应时气膜力与位移的关系。系统在稳态响应期间，转子在不同时间段内的位移波动差值是不断变化的，且转子的运动总是伴随着承载力的变化，Δx增大，承载力的变化率增大。Δx减小，气膜承载力变化率减小。整个系统处于不断变化的动态平衡状态。

图8 系统稳态响应时气膜力与位移关系

综上所述：系统受阶跃负载作用后，转子位移的非线性变化、气膜的流场变化与气膜的承载力变化3者之间相互影响，最终达到1个稳定的平衡状态。

4.2 不同转速下系统特性的耦合计算

当系统进气压力P1＝0.4 MPa，平均气膜间隙0.02 mm，阶跃负载分别为200、300、400 N时，对转子转速n在0～4 000 r／min范围的系统特性进行耦合计算。

图9（a）～（c）表示瞬态响应期间气膜承载力随转子非线性运动的变化规律。系统在平稳状态下受到负载作用，随着时间的增加，不同转速下气膜承载力的波动趋势呈非线性变化，波动幅值随着时间的增加逐渐趋于减小，最终趋于稳定。

图9 瞬态响应期间承载力曲线

图10表示一定转速范围内，系统的瞬态响应参数曲线。图10（a）曲线表明：随着转速的增加，系统的上升时间逐渐增大。图10（b）曲线表明：随着转速的增加，转子到达峰值的时间逐渐增加。图10（c）曲线表明：转子峰值时间所对应的最大位移量随着转速的升高，增加量较小。图10（d）曲线表明：随着转速的升高，系统的最大超调量逐渐减小，变化率逐渐升高。

图10 耦合分析下系统瞬态响应参数曲线

图11表示一定转速范围内系统稳态响应下的波动平均差值。图中曲线表明：系统在不同负载下的波动差值随着转速增加逐渐减小。

图11 耦合分析下系统稳态响应参数曲线

综上所述：转子在瞬态响应期间，随着转子转速的升高，气体随转子周向运动加快，与此同时，转子运动轨迹变化使气膜间隙发生改变，气体流动加剧，系统的快速性降低。系统的Mp值随着转速的增加，逐渐减小，对应的系统的阻尼比ξ值增大，系统的平稳性越高。在稳态响应期间，转子偏心率相对稳定，转速的增加所形成的楔形效应使气膜的承载力增加，系统的自定心作用逐渐增加。

4.3 不同负载下系统特性的耦合计算

当系统进气压力P1＝0.4 MPa，平均气膜间隙为0.02 mm，转子转速分别为0、2 000、4 000 r／min，对阶跃负载在200～400 N范围的系统特性进行耦合计算。

图12（a）～（c）表示瞬态响应期间气膜承载力在不同负载作用下的变化规律：气膜承载力的波动幅值随着负载的增加不断增加，气膜承载力的变化率随负载的增加而逐渐增大，承载力稳定值随着负载的增加而不断提高。

图12 瞬态响应期间承载力曲线

图13表示一定负载范围内，系统在耦合作用下的瞬态响应参数变化。图13（a）曲线表明：随着负载的增加，系统在不同转速下的上升时间逐渐增加。图13（b）曲线表明：随着负载的增加，系统在不同转速下到达峰值的时间逐渐增大。图13（c）曲线表明：随着负载的增加，转子在不同转速下到达峰值的位移量逐渐增大。图13（d）曲线表明：随着负载的增加，系统在不同转速下的最大超调量逐渐减小。

图13 耦合分析下系统瞬态响应参数曲线

图14为系统在稳态响应下的参数变化。曲线表明，随着转速的增加，系统波动的平均差值逐渐减小，且变化率逐渐减小。

图14 耦合分析下系统稳态响应参数曲线

综上所述：转子在瞬态响应期间，受不同的负载作用发生震荡，转子位移的非线性变化，使气膜厚度发生变化，气膜承载力发生非线性变化，负载越大，位移峰值越大。相同转速下，负载增加使气膜楔形效应会得到加强，系统的Mp值随着负载的增加，逐渐减小，对应的系统的阻尼比ξ值增大，系统的稳定性提高。在系统在稳态响应期间：系统在不同外部载荷作用下的转子偏心率相对稳定，偏心作用下形成楔形效应随之增加，负载越大，偏心率越大，楔形效应越强，气膜的高压区面积越大，气膜承载力越强，系统的稳定性越高。

5 结论

1）通过数值模拟方法，利用流固耦合将非线性转子动力学参数与流体计算结果进行交换，采用六面体结构化网格对轴承－转子系统进行网格划分，基于Fluent动网格计算方法，对转子非线性运动下的流场状态进行模拟，实现系统非线性变化状态下的耦合过程计算。

2）基于时域分析方法，结合转子的非线性运动轨迹与气膜力非线性变化，利用瞬态响应参数分析系统的快速性与稳定性；一定范围内，增加转速对系统的峰值位移影响较小；增加转速、增加负载会使系统上升时间、峰值时间增加，影响系统的快速性；在一定转速范围内通过提高转速和负载可以降低系统的超调量，提高阻尼比，有利于提高系统的平稳性。

3）转子波动平均差值反应系统在稳态响应期间的位移波动变化情况，在一定条件下增加负载，提高转速可以使转子波动平均差值减低，提高系统回转精度。