液体静压主轴热态特性多物理场耦合仿真与实验研究*

(中国工程物理研究院机械制造工艺研究所 四川绵阳 621900)

液体静压主轴具有高刚度、高回转精度、高阻尼抗振、长寿命等特性，广泛应用于各类精密超精密加工机床[1-2]。液体静压主轴在工作过程中，由于内置电机的发热和轴承的发热，会导致系统整体温度升高，且会因温度分布不均匀而产生应力集中和变形，从而导致加工误差。在超精密加工过程中，各种误差源都会影响最终的加工精度，其中热误差是其中的主要因素，占总误差的40%～70%[3]。随着制造技术的提升，过去影响较大的几何误差已经逐渐被控制，目前越是精密的机床，热误差影响就越大。因此，如何减少热误差对机床加工精度的影响，改善机床的热态特性已经成为机械制造领域急需解决的问题。

目前，针对液体静压主轴的热态特性，国内外学者已经进行了卓有成效的研究。SATISH等[4-5]利用有限元法数值迭代方法求解雷诺方程，研究了狭缝节流式和小孔节流式液体静压主轴的静态性能和动态性能随温度变化的规律，重点探讨了润滑油黏度随温度变化对主轴性能的影响。LU和MA[6]分析了主轴的热源和散热方式，利用有限元法建立了主轴的热态数学模型，计算了主轴温度场的分布云图。颜超英等[7]以主轴热变形量最小为优化目标，对液体静压主轴的结构参数进行了优化，得到了最优的轴承长径比和封油边的长度。黄智等人[8]基于有限元方法分析了多种参数对主轴系统发热与变形的影响。孙久伟等[9]研究了止推油膜厚度和偏心率对主轴温升和承载的影响。赵春明等[10]基于单向流-固耦合方法，对液体静压主轴的压力场和温度场进行了建模计算。郭传杜和崔怡[11]建立了主轴轴承的单向热-流-固耦合模型，计算了油膜压力场和温度场以及轴承的温度场和变形场。马建刚等[12]针对主轴系统进行了整体的仿真与实验分析。以上针对各类液体静压主轴的各种建模方法并未考虑油膜区域的散热边界条件或对边界条件进行了简化，模型计算准确性欠佳。另一方面，大部分研究由于缺乏实验数据支撑，模型的正确性也无法得到验证。

在前人研究的基础上，针对目前研究中存在的不足，本文作者以一种超精密单点金刚石卧式辊筒加工机床的液体静压主轴为研究对象，考虑油膜区域散热条件建立了“流-热-固”耦合的热态性能仿真模型。首先建立了主轴系统的三维几何模型，分析了其热源以及传热方式，进而推导了主轴的热态性能模型；然后利用一种间接耦合手段，基于有限体积法和有限单元法，以Fluent和Ansys Workbench为仿真平台，分别对油膜和主轴系统进行了仿真计算，分析了工作到稳态下油膜的压力场和温度场，主轴系统的温度场和产生的热变形；同时通过实验平台对仿真结果进行了验证，验证了模型的正确性和准确性；最后应用该模型分析主轴转速、供油压力、液压油黏度、油膜间隙以及轴向封油边长度对主轴温升的影响。

1 主轴热态特性模型

1.1 主轴的几何模型

文中研究的对象是一个电机后置式液体静压主轴，如图1所示。其主要由转子、径向止推节流器、电机、止推板、密封板、壳体等部件组成。液体静压主轴的工作原理是利用一套专用液压站，将润滑油通过管道泵入主轴的径向止推节流器中，液压油经过节流器的节流小孔后，在均压腔中形成了具有一定压力的润滑油膜，利用不同均压腔之间的压力差使主轴浮起。

图1 液体静压主轴系统结构图

主轴工作时主要发热来源于电机损耗和轴承摩擦。文中主轴采用一种径向止推一体化的滑动轴承，其结构如图2所示。其径向与止推均为等面积六腔静压轴承，节流方式采用小孔节流。主轴旋转时，黏性液压油在间隙中摩擦生热。轴承生热模型的相关计算参数如表1所示。下面分别建立不同热源的生热模型以及主轴系统的散热模型。

图2 径向止推一体化滑动轴承结构示意图

参数名称参数值参数名称参数值主轴的半径R/mm60止推轴承内径R1/mm70节流孔直径d/mm0.4止推轴承油腔内径R2/mm85供油压力ps/MPa1止推轴承油腔外径R3/mm115径向轴承油腔的长度l/mm100止推轴承外径R4/mm130径向轴承轴向长度L/mm130止推轴承油腔的夹角θ2/(°)40油腔的深度z/mm2止推油膜厚度h2/mm0.035径向轴承油腔的夹角θ1/(°)40液压油密度ρ0/(kg·m-3)855周向封油边宽度b1/mm10.47液压油动力黏度η/(Pa·s)0.022径向油膜厚度h1/mm0.03液压油比热容C/(J·kg-1·℃-1)2 150油腔数目N6液压油导热系数λ/(W·m-1·K-1)0.127

1.2 主轴生热与散热数学模型

1.2.1 电机生热模型

内置电机的损耗功率可以由下式计算：

(1)

式中：U为电压；I为电流；cosφ为功率因素；η为电机效率。

假设电机损耗全部转换为发热，其中定子产生2/3的热量，动子产生1/3的热量[13]。则电机的生热率为

(2)

式中：P为电机损耗功率；V为热源体积。

1.2.2 油膜生热模型

油膜生热功率包括由主轴旋转引起的摩擦生热功率和液压油供油泵的输出功率，需要分别进行计算。

对于六腔径向轴承，总的内摩擦力可以由下式计算：

(3)

A2=2Rlθ1

式中：A1为封油面的面积；A2为油腔的面积；n为主轴转速；Ff为径向轴承总的内摩擦力。其余参数见表1。

则径向轴承的发热功率为

(4)

对于六腔止推轴承，其发热模型相对于径向轴承更为复杂一些。由于其结构具有周期性，取其一个周期即1/6模型进行分析，如图3所示，可以将1/6模型分为5个区域，4号区域为均压腔，其余区域为止推间隙，由于主轴旋转时不同区域的线速度不同，因此需要分别积分计算每个区域的发热功率。止推油膜的发热功率可以由下式计算：

(5)

式中:Nfi(i=1,2,3,4,5)分别为每个区域的发热功率;Nf为止推油膜1/6模型整体发热功率。其余参数见表1。

供油泵的输出功率等于供油压力与输出流量之积，可以由下式计算：

(6)

式中:ps为供油压力;Q为输出流量;K0为流量系数;pb0为均压腔内压力;β为节流比。其余参数见表1。

图3 止推油膜1/6模型示意图

1.2.3 散热边界条件

主轴系统产生的热量通过其表面和空气、液压油的对流换热散出，不同区域的对流换热系数是重要的边界条件。对流换热在液体静压主轴的实际情况下可以分为自然对流换热、强制对流换热和旋转轴对流换热3种不同的形式，不同区域需要分别进行计算。主轴换热区域可以分为壳体表面与空气的自然对流换热，内部主轴和轴瓦与润滑油的强制对流换热，主轴空心孔与空气的旋转对流换热，以及润滑油在管道内的强制对流换热。

换热系数可以通过努谢尔特准则进行计算，对于不同换热形式努谢尔特数的计算方法不同。不同换热状态的换热系数和相应的努谢尔特数可以由下式计算：

(7)

其中，Nu=0.59(GrPr)0.25(自然对流)

Nu=0.133Re0.67Pr0.33(旋转轴对流)

式中:h为对流换热系数;Pr为普朗克数;Prw为壁面温度时的普朗特数；Gr为格拉晓夫准数;g为重力加速度;Nu为努谢尔特数;λ为流体导热系数;μ为流体的黏度;L为换热壁面定型几何尺寸;β为流体膨胀系数;Re为雷诺数。

由式(1)—(7)，即建立了主轴的生热模型和传热模型，可以依此计算有限元仿真所需的边界条件。

2 主轴系统仿真分析

主轴系统可以分为油膜流体区域和主轴固体区域。油膜区域的温度场十分复杂，并非均匀分布，且缺乏测量手段，过去在仿真过程中往往利用经验公式[14-16]对主轴油膜区域的温升进行计算，并不准确。因此，需要先对油膜区域的温度场进行计算。考虑到计算收敛性与速度的问题，采用间接耦合的方式将油膜的温度场条件赋予到结构场表面，计算主轴结构场的温度场以及产生的热变形。

2.1 油膜流体区域压力场与温度场仿真分析

主轴的止推油膜和径向油膜都是重要的发热源，利用Fluent软件分别对主轴的径向油膜轴承和止推油膜轴承进行仿真。在对流场进行仿真的时候，需建立流场的几何模型并划分网格，如图4所示。在计算时由于止推轴承具有周期性结构，因此取其一个周期即1/6模型进行计算，以减小其计算量；径向轴承因为存在偏心，又要考虑旋转，因此建立其完整模型进行计算。流体计算网格质量要求较高，需要尽量划分结构化网格。网格划分完毕止推轴承1/6模型节点数为3 322 817，单元数为3 172 972，径向轴承完整模型节点数为9 003 386，单元数为8 369 590，网格质量较好，满足计算要求。

图4 油膜区域有限元模型

Fig 4 Finite element model of oil film region(a)regional geometric model of complete oil film;(b) 1/6 finite element model of thrust bearing oil film;(c) finite element model of radial bearing oil film

在Fluent中选择稳态计算，考虑能量方程，考虑油膜黏性发热，液压油相关性质参数按照表1设置，压力入口边界设置为0.7 MPa，初始温度为293.15 K。压力出口设置为0。采用PISO算法，设置压力松弛因子0.5，密度松弛因子0.5，体积力松弛因子1，动量松弛因子1，能量的松弛因子0.5。边界条件采用旋转壁面，转速为250 r/min，热边界条件设为对流换热，换热系数由式(7)计算为42 W/(m2·K)。迭代计算1 000步，最终残差小于10-4，计算收敛。在250 r/min下计算得到油膜的压力场和温度场如图5所示。

图5 油膜区域压力场与温度场仿真结果

由图5(a)、(c)所示的压力场仿真结果可以看出小孔和狭缝明显的节流效果，在小孔和封油边均有明显的压降，均压腔内压力十分稳定，确实起到了均压的效果，压力从入口到出口逐渐降低至大气压。径向油膜设置有5 μm的偏心量，可以明显看出下部油腔压力高于上部油腔，从而起到支承作用，在该条件下计算得到的承载力为693.59 N。

由图5(b)、(d)所示的温度场仿真结果可以发现，止推轴承主轴温升集中在内径和外径边缘，且由于内外线速度的不同，导致外径温升高于内径的温升。均压腔内温度均匀且温升较小，由于均压腔的深度远远大于油膜的厚度，因此导致均压腔内发热量远远小于油膜区域，因此温升较小，这与式(5)的计算结果吻合，整个止推油膜的平均温升为4.42 ℃。径向油膜由于线速度一致，因此温度场更加均匀，其均压腔内的温升依旧低于油膜区域，入口处温度最低，温升最大的区域集中在封油边附近，整个径向油膜的平均温升为2.41 ℃。由于止推油膜的平均直径大于径向轴承，径向油膜的平均温升低于止推油膜的平均温升。将Fluent仿真结果作为边界条件附加在结构场中，可以对整机的温度场和变形进行仿真计算。

2.2 主轴固体区域温度场与结构场仿真分析

对于主轴的整机模型，其主要热源来自于电机损耗和油膜的黏性耗散。依据前文的分析，建立主轴系统的有限元模型，如图6所示。经统计网格单元数为3 170 673，节点数为5 177 534，网格质量良好，满足计算要求。

图6 主轴系统有限元模型

由式(1)、(2)计算得到电机的发热功率为17 296.2 W/m3，在转速250 r/min条件下计算得到径向油膜的发热功率为11.54 W，止推油膜的发热功率为34.76 W。油膜区域对流换热系数为42 W/(m2·K)，主轴壳体表面与空气的自然对流换热系数为7.5 W/(m2·K)，转子内通孔与空气的旋转对流换热系数为15 W/(m2·K)，液压油在管道内强制对流换热系数为20 W/(m2·K)。

主轴系统初始温度为20 ℃，环境温度为20 ℃，管道内液压油温度设为液压站出口温度19.4 ℃，油膜区域温度设置为Fluent仿真计算所得温度。分别设置好边界条件，采用稳态求解模块计算主轴的温度场，结果如图7所示。主轴的整体最高温升出现在电机的定子处，温升2.75 ℃，最低温升出现在前密封板处，温升0.579 ℃。观察主轴壳体表面温度分布情况，电机与轴承对应位置为温升最大区域，其中后轴承靠近电机，因此温升比前轴承更高；进油口处由于有低温的液压油通过，温升相对较低。

图7 主轴系统温度场仿真结果

将主轴系统温度场作为热载荷导入到结构场中，考虑重力影响，计算主轴系统在温升下的变形，结果如图8所示。主轴系统的总体最大变形量为7.731 μm，其中轴向变形为主轴系统的主要变形，最大轴向变形量为5.286 μm。

图8 主轴结构场仿真结果

3 试验验证

为了验证前文建立的主轴系统完整的热态性能分析模型以及仿真结果的正确性，设计了稳态测量试验，对主轴到达稳态时关键点的温升以及轴向的变形量进行测量。设计的稳态测量实验，主轴前端夹持标准棒，液压油通过外接油冷机保持入口处油温恒定为19.4 ℃，环境温度为20 ℃，主轴在稳定转速下持续运转，直至达到稳定状态。在此过程中，同时测量主轴的温度变化和轴向发生的热变形量。

采用PT100电阻温度传感器在主轴壳体表面选点进行温度测量。PT100传感器测量范围为-50～150 ℃，分辨率为0.1 ℃，满足测量要求。在温度测点的选择上，选择靠近径向轴承和止推轴承位置，温度传感器的布点位置如图9所示，其详细说明如表2所示。

图9 主轴温度测点示意图

表2 主轴温度测点

对于主轴轴向变形量的测量，由于主轴旋转，选择了非接触式电容位移传感器进行测量，其量程为800 μm，分辨率为0.01 μm，精度为0.4 μm，满足测量要求。温度传感器与位移传感器的信号通过信号采集器进行采集和记录。温度传感器采样频率为1 Hz，位移传感器采样频率为50 Hz。实验装置的具体设置如图10所示。

图10 实验测量设备示意图

实验室环境为20 ℃恒温，主轴转速设为250 r/min。试验3 h左右，温度传感器以及位移传感器的数据不再变化，主轴达到热平衡状态。稳态温度和位移的测量结果如图11所示。可以看出壳体表面确实出现了明显的温升，初期温升较快，后逐渐平缓直至稳定。比较不同测点的温升结果，可以看出由于电机的影响，靠近电机一侧的壳体表面温度会略高于远离电机的一侧，证明了不同热源之间的耦合作用。在对称位置上，其温升基本相同。最高温度在靠近电机的壳体顶部，此处距离止推油膜最高温升区十分接近，又靠近电机，距热源最近，因此温升最高。

图11 温度与位移试验测量结果

从图11所示的位移测量结果来看，主轴到达稳态时，轴向伸长量为5.41 μm，其结果与仿真结果误差为2.5%。从图12可以看出，达到稳态情况下，温度测量结果与温度场仿真结果相比，各点的温升的最大误差为17.5%，最小误差为0.3%，平均误差10.5%。试验结果验证了仿真结果的正确性，可以利用此模型对不同工况下的主轴温升与变形进行预测。

图12 温度场仿真结果与试验结果比较

4 主轴热态性能分析

液体静压主轴的温升对其加工精度的影响很大，在满足主轴性能的前提下，主轴设计时需要选择设计参数使得主轴的发热量最小。通过理论模型分析可知，影响电机发热的主要因素是主轴的转速，影响主轴油膜发热量的主要因素包括主轴转速、供油压力、液压油黏度、油膜间隙以及轴向封油边的长度。

根据建立的理论模型，采取控制变量法研究上述参数对主轴热态性能的的影响，以径向轴承为对象进行分析，结果如图13所示。

由图13(b)、(c)可以看出，提高供油压力、选择较小黏度的液压油，可以有效地降低油膜区域的温升。其原因是在供油压力增大、液压油黏度减小的情况下，液压油的流量会增大很多，而液压油的温度较低，实际上对主轴系统起到了冷却作用，因此其流量越大，从主轴区域带走的热量也就越多，主轴温升也就越低。由图13(a)、(d)可以看出，在设计轴承时选择较大的油膜间隙、较短的封油边，也可以有效地降低主轴的发热量。这是因为较大的油膜间隙使得液压油流量增大，从主轴区域带走的热量越多，主轴温升也就越低；较短的封油边可以使得油膜剪切区域的面积减小，从而发热量减小。由图13(e)可以看出，随主轴转速的增大，主轴温升先缓慢增大后快速增大，因此在满足主轴性能的前提下，应采用较低的主轴转速。

图13 不同设计参数对主轴平均温升的影响

综上，在液体静压主轴的设计阶段，应当综合考虑各项因素对主轴性能的影响，选取最优参数，而后通过仿真的手段进一步验证，再进行主轴的加工，从而最大限度上保证主轴的性能。

5 结论

(1)分析主轴发热散热机制，指出电机和油膜是主轴的主要热源，两热源相互耦合影响，在主轴工作时会导致主轴的温升和形变，影响主轴的加工精度。

(2)建立了较为准确的主轴系统的温度场与结构场仿真模型，通过该模型可对不同工况下的主轴温度场以及变形量进行预测。该模型的温度场仿真结果与试验结果的平均误差为10.5%，结构场仿真结果与试验结果误差为2.5%，试验结果和仿真结果具有较好的一致性，验证了仿真结果的正确性。

(3)采用该模型分析主轴转速、供油压力、液压油黏度、油膜间隙以及轴向封油边长度对主轴温升的影响，结果表明，主轴温升随主轴转速、液压油黏度、封油边长度的增大而增大，随供油压力、油膜间隙的增大而减小。因此，为了控制主轴的温升，应当在满足主轴性能的前提下，使用较低的主轴转速、较大的供油压力，选择小黏度的液压油、较大的油膜间隙以及较短的封油边长度。