基于有限元法的滚珠丝杠进给系统热特性分析*

黄 俊，袁军堂，汪振华，董香龙

(南京理工大学 机械工程学院，南京 210094)

0 引言

滚珠丝杠副传动系统，因具有传动精度高、传动效率高、传动的可逆性、同步性能好、使用寿命长等优点而被广泛应用于数控机床及加工中心，其热变形将严重影响机床的加工精度。研究表明由机床热变形引起的误差已占机床全部误差40%～70%［1］。为保证机床的加工制造精度，需要对滚珠丝杠的热变形进行误差补偿设计。因此对实际工况下的滚珠丝杠进给传动系统进行有限元热仿真分析，对于提高机床的整体精度具有指导意义。

国内外已对滚珠丝杠进给系统热分析进行了相关的研究:易学平［2］对机械系统中的丝杆机构在工作过程中摩擦引起的温升和热变形原因进行分析，得出了适用于实际工况的温度场经验公式。芮执元［3］建立中空丝杠传动过程中温度场和热变形的数学模型，通过有限元仿真模拟和丝杠冷却模拟运行试验验证该模型。王大伟［4］考虑丝杠导程引起表面积变化及螺母移动对温度场的影响，建立滚珠丝杠传动过程中的温度场和热变形模型。刘兴业［5］、曹巨江［6］对加工中心滚珠丝杠传动系统进行温度场和热变形分析，得出进给系统温度场分布及变形情况。Huang［7］用丝杠螺母、前后轴承等处热源的温度作为模型的变量，用多元线性回归的方法预测了滚珠丝杠在不同转速下的热变形。Kmisk［8］等通过有限元法建立机床滚珠丝杠进给系统的温度场。Ahn［9］结合传热学理论提出一维传热模型和分段模型来实时预测系统的温度场变化及热变形情况。

本文以滚珠丝杠直线进给系统为研究对象，在ANSYS 中对不同进给速度和不同冷却液流量下的温度场进行特性分析，并重点分析在进给速度30m/min条件下进行热变形分析，研究出滚珠丝杠整体热变形规律，为进给系统热补偿提供依据。

1 滚珠丝杠进给系统热分析模型

滚珠丝杠进给系统的主要热源是螺母的移动摩擦生热和轴承的旋转摩擦生热。丝杠与外界的热交换形式主要为热传导。

1.1 滚珠丝杠导热理论模型

在导热过程中，单位时间内通过给定截面的导热量，正比于垂直该截面方向上的温度变化率和截面面积，热量传递的方向则与温度升高的方向相反，这就是导热基本定律即傅里叶导热定律［10］。傅里叶导热定律用热流密度q表示时有下列形式:

q—热流密度，W/m2;λ—导热系数，W/m·℃;

t—物体的温度，℃;λx，λy，λz—材料沿x，y，z三个方向的导热系数;nx，ny，nz—等温面沿x，y，z三个方向的分量;gradt—空间某点的温度梯度。

丝杠工作工程中没有内热源，载荷为轴对称形式，满足热传导微分方程:

丝杠温度分布初始条件为

式中:T0—初始温度，丝杠与工作环境的初始温度均取293K。

1.2 发热量及边界条件计算

(1)轴承发热量主要是由轴承的摩擦力矩引起的，其发热量可由下式计算:

式中:n为丝杠的转速(r/min);M为滚动轴承的摩擦力矩(N/mm)。

其中M = M0+ M1

式中:M0—与轴承类型、转速和润滑油性质有关的力矩;M1—与轴承所受负载有关的力矩。

(2)丝杠螺母发热量与轴承类似，也是主要由摩擦力矩引起，其计算公式:

式中:n为丝杆的转速(r/min);M为滚珠丝杆螺母副的摩擦力矩(N/mm)。

其中M =2z(Me + Mg)cosβ

式中:Me—摩擦阻力矩;Mg—几何滑移摩擦力矩;z—滚动体数目;β—滚珠丝杆滚道的螺旋角。

式中mα，mβ—与球接触变形椭圆区偏心率有关的系数;Q—单个球体所受径向压力(N);f—滑动摩擦系数;μ1、μ2—分别是材料的泊松比;E1，E2—分别是材料的弹性模量;分别是滚道和滚珠的曲率半径;

ρ11，ρ12，ρ21，ρ22分别是两滚动体的主曲率。

(3)丝杠与空气的热传导

丝杠以一定速度旋转且丝杠具有螺旋槽特征，故其对流放热的效率应比等直径光杆要大。其准则方程为

式中:Re—雷诺数，—丝杠角速度;d—丝杠直径;Pr—普朗特数。

(4)丝杠热变形基本方程

根据材料热膨胀变形理论，滚珠丝杠进给系统中丝杠的热变形主要有温升值及热膨胀系数α 决定:

式中:α—丝杠的线膨胀系数，K-1;Δt—温度变化值;l—丝杠的长度。

结合公式(1)～(6)，及丝杠与冷却介质热对流温度分布表达式［3］，计算出不同进给速度和冷却流量下的边界条件，q1、q2分别为轴承和丝杠螺母的热流密度，如表1 和表2 所示。

表1 不同进给速度下的边界条件

表2 不同冷却液流量下的强制对流系数

2 滚珠丝杠进给系统有限元分析

滚珠丝杠进给系统的热态特性分析主要包括系统的稳态温度场分析、瞬态温度场分析以及热-结构耦合场分析。

2.1 建模及网格划分

先通过Solidworks 建立滚珠丝杆系统的模型，再导入ANSYS，其三维简化模型如图1 所示。建模时，用等直径的空心光杆代替滚珠丝杠，用等直径的空心圆柱体代替轴承以简化模型。

图1 滚珠丝杠进给系统有限元模型

绝大部分丝杠的材料主要为铬钢，其属性如下:弹性模量E=210GP，泊松比μ =0.3，密度ρ =7800kg/m3，热膨胀系数α =1.2e -5 1/℃，热导率λ =30W/(m·℃)，比热C=450J/(kg·℃)，各零件结合面的热接触传导系数= 3200(w/℃)。选择热实体单元SOLID70 划分网格。采用具有热结构特性的面接触对CONTA174 和TARAE170 来模拟丝杠系统中丝杠与螺母及轴承之间的热传导。由于丝杠的不规则，采用自由法划分网格。

2.2 移动载荷的施加

机床在空载进给时，螺母一直往复移动，为得到真实的温度场，运用APDL 语言在丝杠上施加往复移动的热源载荷，螺母每移动一次，在前一位置处删除热流密度，施加对流载荷，再在新的位置处删除其对流载荷并施加热流密度，进行求解。如此往复进行，直到丝杠、螺母和轴承温度不再升高，达到热平衡为止。

3 计算结果分析

3.1 温度场、热与结构耦合分析

以进给速度30m/min 为例，螺母循环5000 次时，对丝杠进给系统的温度变化进行定量的仿真计算，验证丝杠传动过程中的温度变化理论。

由上图2 可以看出进给系统各部分的最高温升分布不均匀，滚珠丝杠系统中两端轴承温度最高，约为45℃，温升为25℃，滚珠丝杠中间处温度相对较低，约为25℃，温升为5℃，且温度由中间向两端逐渐增大。由温度不均引起的热变形也不一样，容易产生过大的加工误差，可以通过采取散热装置使得热负载均匀。

在实际工况下，随着丝杠表面温度的升高，其与外界环境的温差增大，故对流放热系数也增大，当单位时间内流入的热量和流出的热量相等时，系统达到热平衡，并保持稳定状态。丝杠达到动态热平衡的时间与丝杠材料的传热系数、丝杠周围空气的流动速度、丝杠冷却方式、丝杠两端轴承的润滑冷却系统等因素有关。

由图3 可知，螺母需要大约800s 达到热平衡状态，轴承在前500s 温升很快，1000s 后温度逐渐趋于热平衡。因此，在结构设计时，采用合理的润滑和冷却系统，轴承的温升会迅速的减小，减小其对丝杠系统的影响。

图2 螺母循环5000次的温度场

图3 螺母、两端轴承温升曲线

同时在丝杠工作导程部分各点的温度分布和变形情况，如图4 和图5。可以看出，丝杠两端受轴承的影响温升较大，若考虑丝杠两端支撑轴承的约束，调入温度载荷分析，丝杠两端的节点位移远大于中间节点。

图4 丝杠不同节点的温度

图5 丝杠节点位移变化曲线

图6 中，丝杠总的伸长量约为25.4μm，由数学模型中Δl =α×Δt × l计算得到丝杠的伸长量为25.25μm，从而说明分析模型的准确性。

图6 滚珠丝杠热位移云图

3.2 进给速度对平衡时间的影响

在未通冷却液时，不同空载进给速度下，丝杠轴上温升最高点的温升曲线如图7 所示，30m/min，20m/min、10m/min 空载进给时热平衡的时间分别是1200s，1500s，2400s，进给速度越快，达到热平衡的时间越少。

3.3 冷却液对平衡时间的影响

中空丝杠一般只用于高转速、重负载、大传动力和高精度的工作状况，当以高速空载进给速度30m/min进给时，分别通入1L/min，2L/min，5L/min 流量的冷却油，其中油层对应的状态分别是层流、过渡、紊流状态。经过加载往复移动热源载荷后，得到丝杠温升最高点处的温升曲线，见图8。

图7 不同进给速度时系统的温升曲线

图8 不同冷却液流量下的温升曲线

当流量分别以2L/min 和5L/min 时，二者的温差值控制在2℃之内，故较适宜的冷却液流量为2L/min，丝杠轴表面的热平衡温度显著下降到13℃，达到热平衡时间仅为200s。

4 结论

根据对滚珠丝杠进给系统热有限元仿真，获得螺母、轴承的温升变化曲线，通过对丝杠各节点温度变化的对比，可以得出轴承的发热对丝杠温升影响较大，为设计合理的轴承系统提供参照。以移动载荷分析的热结果，进行热-结构耦合分析，得到丝杠节点的热变形曲线和总热伸长量，为丝杠的热补偿提供依据。对不同进给速度和冷却流量的定量分析，为进给系统的热平衡提供运行参数。

［1］陈兆年，陈子辰.机床热态特性基础［M］.北京:机械工业出版社，1989.

［2］易学平.丝杠结构热变形及近似估算探讨［J］.机械工程与自动化，2011，166(3):170 -171.

［3］芮执元，张传辉，郭俊锋.基于有限元法的中空冷却滚珠丝杠温升及热变形研究［J］. 机械设计，2011，28(12):17 -20.

［4］王大伟，刘永红，张龙，等.基于有限元法的滚珠丝杠传动过程中的温度场和热变形仿真［J］. 计算机辅助工程，2009，18(2):29 -33.

［5］刘兴业，殷国富，刘立新，等. 龙门加工中心滚珠丝杠传动系统的温度场和热变形分析［J］. 组合机床与自动化加工技术，2011，4(4):1 -4.

［6］曹巨江，李龙刚，刘言松，等. 基于有限元法的高速空心滚珠丝杠系统热特性仿真［J］. 组合机床与自动化加工技术，2010(12):21 -24.

［7］Shyh-Chour Huang. Analysis of a model to forecast thermal deformation of ball screw feed drive systems［J］. International Journal of Machine Tools and Manufacture，1995，35(8):1099 -1104.

［8］S. K. Kim，D. W. Cho，Real-time estimation of temperature distribution in a ball-screw system［J］. International Journal of Machine Tools ＆ Manufacture，1997，37(4):451 -464.

［9］Ahn，J.Y.，Chung，S.C. Real–time estimation of the temperature distribution and expansion of a ball screw system using an observer［C］. Proceedings of the Institution of Mechanical Engineers，Part B:Journal of Engineering Manufacture，2004，218 (12):1667 -1681.

［10］杨世铭，陶文铨. 传热学［M］. 北京:高等教育出版社，2006.