一种数控机床滚珠丝杠进给系统热误差预测方法

摘要

本发明属于热误差预测方法技术领域，提出了一种数控机床滚珠丝杠进给系统热误差预测方法。利用自适应实时模型(ARTM)预测滚珠丝杠进给系统的热误差。本方法首先建立实验测量工作台表面温度随时间的变化，然后建立一种自适应实时模型用来预测滚珠丝杠的瞬时温度分布和热误差分布。利用有限元结合蒙特卡洛方法(MC)确定滚珠丝杠进给系统两个轴承、移动丝母、滑轨的发热率。然后依据有限元计算数据，建立进给速度与时间的指数方程，用来描述测量表面点和运动对中心之间的温度差随时间的变化，最后建立一种数值预测算法用来预测滚珠丝杠进给系统的热误差。通过监测两个轴承座和运动螺母侧面的表面温度，使用数值预测算法快速高精度地预测相应的热误差。

说明书

技术领域

本发明属于热误差预测方法技术领域，涉及数控机床滚珠丝杠进给系统热误差建模方法。

背景技术

滚珠丝杠作为数控机床进给系统的关键部件，工作时轴承、丝杠螺母等部件摩擦生热引起丝杠温度的上升，其温度的变化导致机械结构的变形,从而产生热误差降低了进给系统的定位精度。由于热误差的存在，单件的加工精度易不合格；批量加工产品的一致性较差。因此，机床热误差问题一直是近期国内外研究的重点，已经成为了提高数控机床定位精度的重要研究方向。

近年来，国内外研究学者建立了几种滚珠丝杠热误差模型的有效方法，但是都存在其局限性。相关研究主要集中两个方面：第一种，对滚珠丝杠进给系统的结构进行简化，应用传热学理论为基础近似计算系统的热边界条件，建立滚珠丝杠进给系统热误差有限元模型。第二种，实测大量数据进行处理，建立滚珠丝杠进给系统热误差经验模型。通过对目前研究的分析。

发明内容

为了克服现有技术中存在的问题，本发明提出了一种数控机床滚珠丝杠进给系统热误差预测方法。该方法利用自适应实时模型(ARTM)预测滚珠丝杠进给系统的热误差。本方法首先建立实验测量工作台表面温度随时间的变化，然后建立一种自适应实时模型用来预测滚珠丝杠的瞬时温度分布和热误差分布。利用有限元结合蒙特卡洛方法(MC)确定滚珠丝杠进给系统两个轴承、移动丝母、滑轨的发热率。然后依据有限元计算数据，建立进给速度与时间的指数方程，用来描述测量表面点和运动对中心之间的温度差随时间的变化，最后建立一种数值预测算法用来预测滚珠丝杠进给系统的热误差。

本发明的具体技术方案为：

步骤1，滚珠丝杠进给系统温度测量

确定进给系统表面温度检测点，其中，温度检测点包括两个丝杠支撑轴承的表面点、丝母法兰的表面点和丝杠轴一侧导轨表面均匀分布点；T型贴片热电偶通过磁性底座分别固定在两个丝杠支撑轴承的表面点、丝母法兰的表面点上，按周期进行连续采样；通过红外线测量仪对测量点进行定时检测；

步骤2，建立其有限元模型；

步骤2-1，根据目标数控机床进给系统的几何尺寸，建立其有限元模型，模型包括伺服电机、两个支撑轴承、丝杠、螺母和两个滑轨，各部分均选择固体单元，而丝杠与两个支撑轴承、丝母之间的连接选用接触单元；

步骤2-2，确定系统的热源；系统包括固定热源和移动热源；固定热源为伺服电机、轴承一和轴承二；依据系统结构，将伺服电机发热率合并到轴承一，记为Q_B1，轴承二的发热率记为Q_B2；移动热源为两个滑轨和一个丝母，在此只考虑各热源传入模型部分的发热率，两个滑轨发热率为Q_g1＝Q_g2＝Q_g，丝母的发热率Q_n；

步骤2-3，对流换热系数的确定；

在丝杠旋转过程，丝杠表面与空气产生对流传热，则对流换热系数为：

其中，N_u为努塞尔数；λ_fluid为空气的导热系数；d为丝杠的直径。

步骤3，用蒙特卡洛方法实现发热率的有限元计算；

利用x₁，x₂，x₃和x₄分别表示发热率Q_B1，Q_B2，Q_n和Q_g，则蒙特卡洛模拟计算为优化如下的目标函数：

其中，T_ij是在第j个采样时间步长的给定点i的温度，i＝1,2,…,15，j＝1,2,…,N，上标EM和MC分别代表实验测量值和有限元结合蒙特卡洛的模拟值；

蒙特卡洛方法具体为：

首先确定采样次数为k，并给定每个参数x_l的取值[x_Ll,x_Ul]区间；在每一次采样过程中，确定各变量的数值：

x_l＝x_Ll+random(0,δ_l),δ_l＝x_Ul-x_Ll,l＝1,2,3,4.(5)

每次采样后进行有限元的计算，提取各检测点的温度，依据式(4)求出目标函数F，并将k次采样计算的最小目标函数记为F^*＝min(F₀,F₁,…,F_k-1)，而对应参数记为和完成k次蒙特卡洛模拟后，将F^*与给定的计算精度ε进行比较，如果达到给定计算精度，则输出发热率和结束模拟计算；否则，判断取值区间大小：

min(δ₁,δ₂,δ₃,δ₄)<δ₀(6)

其中，δ₀为给定常数。

如果式(6)满足，则计算失败，计算过程结束；否则，修改各变量的随机抽样区间如式(7)，并进行重新模拟计算。

步骤4，运动副中心与表面温度关系的确定

描述测量表面点和运动副对中心温度差：

ΔT(t,v)＝av(1-e^-t/bv)(8)

其中，v是工作台的进给速度，t是工作时间，a和b是通过曲线拟合得到的常数；

在机床工作过程中，检测支撑轴承或丝母表面的温度，计算得到其中心温度：

T_c＝T_m+ΔT(t,v)(9)

其中，T_c运动对的中心温度，T_m是其表面上给定点的测量温度；

步骤5，滚珠丝杠实时热误差预测；

步骤5-1建立滚珠丝杠热传导模型

将丝杠简化为一维杆，其热传导方程：

其中，T(x,t)是时间t和位置x的函数，表示t时刻丝杠上距离热源为x的位置点的温度变化；κ是导热系数；ρ是滚珠丝杠的材料密度；c是比热容；h是对流系数；T_f(t)是工作环境空气的温度；

对滚珠丝杠T(x,t)的函数划分空间步长s和时间步长τ，x_k＝k·s，k＝1,2,…,M，t_j＝j·τ，j＝1,2,…,N；用表示一阶导数和二阶导数的有限差分方程如下：

把差分方程式(11)和式(12)代入热传导方程式(10)，得到：

其中A＝ρc/κ；B＝4h/(κd)，当(As²-Bs²+2τ)/(As²)≥0且(τ/s²)≤(A-B)/2时，等式成立；

依据式(9)，两个轴承中心的温度如下式表示：

T_cbl＝T_mbl+ΔT_bl(t,v)l＝1,2.(14)

移动螺母中心的温度表示为：

T_cn＝T_mn+ΔT_n(t,v)(15)

步骤5-2，滚珠丝杠热误差计算；

在加工过程中，由于摩擦热引起的内部温度上升导致滚珠丝杠轴的伸长；热伸长值如下式进行计算:

ΔL＝L·α·ΔT(16)

其中，L是滚珠丝杠的长度，因此滚珠丝杠在X轴的热伸长量：

其中ΔL(x)表示t时刻的热延伸量，ΔT(x,t)表示滚珠丝杠的温度变化，α是线性热膨胀系数。

本发明的有益效果：

本文提出利用蒙特卡洛模拟所得发热率有限元计算，可用于预测滚珠丝杠进给驱动系统的发热率，温度分布和热误差。提出的自适应实时模型(ARTM)可用于预测滚珠丝杠进给驱动系统的热误差。通过监测两个轴承座和运动螺母侧面的表面温度，使用数值预测算法快速高精度地预测相应的热误差。

附图说明

图1是滚珠丝杠进给系统温度测量点位置图；

图2是结合蒙特卡洛法的有限元模拟流程图；

图3是试验检测温度与利用蒙特卡洛法进行有限元计算所得温度结果变化图；EM代表实验测量数据，FEM代表有限元计算数据；

图4是轴承1在三种进给速度下温度差计算有限元结果和拟合曲线的对比。

图5是滚珠丝杠进给系统热源简化图；

图6是滚珠丝杠温度测量点位置图；

图7是进给速度为10m/min时，有限元方法；自适应实时模型预测计算；试验结果三者之间的温升对比图；

图8是进给速度为10m/min时，有限元方法，自适应实时模型预测计算，试验结果三者之间热误差对比图。

具体实施方式

以下结合附图和具体实施例对技术方案进行详细阐述。

一种数控机床滚珠丝杠进给系统热误差预测方法，包括如下步骤：

步骤1，滚珠丝杠进给系统温度测量

确定进给系统表面温度检测点，如图1所示。其中，T1和T8为两个丝杠支撑轴承表面点，T9为丝母法兰表面点；其余在丝杠轴一侧导轨表面均匀分布点，如T2-T7和T10-T15。用磁性底座分别将三个T型贴片热电偶固定在T1，T8和T9点，有一台笔记本电脑对其以20s的采样周期进行连续采样。使用高精度红外线测量仪对测量点T2-T7和T10-T15进行定时检测。

步骤2，建立进给系统的有限元模型

步骤2-1，所研究数控机床进给系统的几何尺寸，利用ANSYS软件建立其有限元模型。模型包括伺服电机、两个支撑轴承、丝杠和螺母和两个滑轨五部分，各部分均选择SOLID87固体单元，而丝杠与两个支撑轴承、丝母之间的连接选用CONTA174接触单元。

步骤2-2，确定系统的热源。系统包括三个固定热源和三个移动热源。三个固定热源为伺服电机、轴承1和轴承二。依据系统结构，将伺服电机发热率合并到轴承1，记为Q_B1，轴承二的发热率记为Q_B2。三个移动热源为两个滑轨和一个丝母，在此只考虑各热源传入模型部分的发热率，两个导轨发热率为Q_g1＝Q_g2＝Q_g，丝母的发热率Q_n。

步骤2-3，对流换热系数的确定。

在丝杠旋转过程，丝杠表面与空气产生对流传热，丝杠表面空气的雷诺数为：

其中，ω为丝杠旋转角速度，ν_fluid为空气中的运动粘度，d为丝杠的直径。

努塞尔数为：

其中，P_r是普朗特数：

其中，C_fluid为空气的比热容，μ_fluid为空气的动力黏度，λ_fluid为空气的导热系数，则对流换热系数为：

因为床鞍的进给速度相对较低，空气流速对对流换热系数的影响可以忽略，所以床鞍和轴承座和壳体的表面使用12W/m²K的对流换热系数。

步骤3，用蒙特卡洛方法实现发热率的有限元计算

有限元计算目标是依据步骤1检测得到的表面温度确定各热源的发热率。在此，采用蒙特卡洛模拟实现系统各热源发热率的计算。为方便描述，利用x₁，x₂，x₃和x₄分别表示发热率Q_B1，Q_B2，Q_n和Q_g，则蒙特卡洛模拟计算为优化如下的目标函数：

其中，T_ij是在第j个采样时间步长的给定点i的温度，i＝1,2,…,15，j＝1,2,…,N，上标EM和MC分别代表实验测量值和有限元结合蒙特卡洛的模拟值。

蒙特卡洛模拟是一个采样过程。为了保证模拟计算精度，首先确定一个足够大的采样次数为k，并给定每个参数x_l的取值[x_Ll,x_Ul]区间。在每一次采样过程中，利用ANSYS>

x_l＝x_Ll+random(0,δ_l),δ_l＝x_Ul-x_Ll,l＝1,2,3,4.(5)

每次采样后进行有限元的计算，提取各检测点的温度，依据式(4)求出目标函数F，并将k次采样计算的最小目标函数记为F^*＝min(F₀,F₁,…,F_k-1)，而对应参数记为和完成k次蒙特卡洛模拟后，将F^*与给定的计算精度ε进行比较，如果达到给定计算精度，则输出发热率和结束模拟计算；否则，判断取值区间大小：

min(δ₁,δ₂,δ₃,δ₄)<δ₀(6)

其中，δ₀为给定常数。

如果式(6)满足，则计算失败，计算过程结束；否则，修改各变量的随机抽样区间如式(7)，并进行重新模拟计算。

蒙特卡洛方法实现发热率的有限元计算程序流程图如图2所示。

步骤4，运动副中心与表面温度关系的确定。

对某数控车床X轴进给系统进行温度检测试验，其进给速度分别为5m/min、10m/min和15m/min。利用步骤2计算方法对试验检测得的温度数据进行计算，得到不同进给速度各热源的发热率。再将其代入的有限元模型进行计算，提取两个轴承和丝母中心及表面检测点的温度，如图3所示。由图3可见，表面检测温度与利用蒙特卡洛模拟所得发热率有限元计算所得结果吻合比较好，证明了步骤2算法的有效性。比较图3中各运动副表面温度与中心温度的差值，可以用式(8)描述测量表面点和运动副对中心温度差：

ΔT(t,v)＝av(1-e^-t/bv)(8)

其中，v是工作台的进给速度，t是工作时间，a和b是通过曲线拟合得到的常数。表1为利用有限元计算结果进行曲线拟合所得的轴承一和二、丝母的参数及对应的最大误差。图4为轴承1在三种进给速度下温度差计算有限元结果和拟合曲线的对比。由图4和表1可看出，式(8)的进给速度和时间的幂函数能够很好地描述运动副中心与表面检测点之间的温度差。

表1指数拟合曲线的参数和测量表面点与运动对中心之间的最大温度误差

在机床工作过程中，只要检测支撑轴承或丝母表面的温度，即可计算得到其中心温度：

T_c＝T_m+ΔT(t,v)(9)

其中，T_c运动对的中心温度，T_m是其表面上给定点的测量温度。

步骤5，滚珠丝杠实时热误差预测

步骤5-1建立滚珠丝杠热传导模型。

如图5所示，滚珠丝杠进给系统中的支承轴承和移动丝母是滚珠丝杠的热源，导致其滚珠丝杠的热变形。滚珠丝杠轴长度远大于其直径，为了便于研究，忽略径向温差，并且滚珠丝杠进给系统中，向下一端固定向上一端自由，如图1和图5所示，所以将丝杠简化为一维杆，其热传导方程：

其中，T(x,t)是时间t和位置x的函数，表示t时刻丝杠上距离热源为x的位置点的温度变化；κ是导热系数；ρ是滚珠丝杠的材料密度；c是比热容；h是对流系数。T_f(t)是工作环境空气的温度。

对滚珠丝杠T(x,t)的函数划分空间步长s和时间步长τ，x_k＝k·s，k＝1,2,…,M，t_j＝j·τ，j＝1,2,…,N。用表示，所以一阶导数和二阶导数的有限差分方程如下表示：

把差分方程式(11)和式(12)代入热传导方程式(10)，得到：

其中A＝ρc/κ；B＝4h/(κd)，当(As²-Bs²+2τ)/(As²)≥0且(τ/s²)≤(A-B)/2时，等式成立。

依据式(9)，两个轴承中心的温度可以如下式表示：

T_cbl＝T_mbl+ΔT_bl(t,v)l＝1,2.(14)

移动螺母中心的温度可以表示为：

T_cn＝T_mn+ΔT_n(t,v)(15)

步骤5-2，滚珠丝杠热误差计算；

在加工过程中，由于摩擦热引起的内部温度上升导致滚珠丝杠轴的伸长。热伸长值如下式进行计算:

ΔL＝L·α·ΔT(16)

其中，L是滚珠丝杠的长度，因此滚珠丝杠在X轴的热伸长量：

其中ΔL(x)表示t时刻的热延伸量，ΔT(x,t)表示滚珠丝杠的温度变化，α是线性热膨胀系数。

步骤6，数控机床进给系统热误差预测方法验证试验；

为了验证热误差预测方法的准确性，对一种数控机床进行了试验，滚珠丝杠的参数如表2所示。温度测量试验已经在步骤1中进行了详细阐述。然后使用激光干涉仪测量X轴定位误差。机器参考原点设置为该测量的起始点，并且整个行程范围为220mm，确定进给系统表面位移误差检测点，如图6所示，即位于点P1-P6的工作台上，每隔44mm测量定位误差。首先，初始几何误差是在机床初始接通时常温下测量的。然后，X轴工作台以给定的进给速度沿其行程范围移动，此外，滚珠丝杠和两个导轨用热敏摄相机以10分钟的间隔进行拍摄。在机床加热10分钟，20分钟，30分钟，40分钟和50分钟后，同步测量定位误差。在这项研究中，分别以5,10和15m/min的进给速度进行了三次空转试验，每次试验在环境温度条件下启动，工作台在整个行程范围内前后移动。在工作台的工作过程中，由于应用程序的运行速度，计算机的位置采样周期为48ms，所以步骤3的时间步长假定为τ＝0.048s。

表2数控机床X轴进给系统参数规格

如图7所示，是进给速度为10m/min工作过程中，有限元方法；ARTM(自适应实时模型)；试验结果三者之间的温升对比图，从图7可以看出，自适应实时模型预测计算和有限元计算的结果与试验测得的结果十分吻合。

如图8所示是ARTM和FEM计算的定位误差结果与试验测量结果，可以看出两个模型计算的定位误差结果与实际测量误差结果十分吻合，并且他们的精度都满足要求。但是有限元计算时间较长，但自适应实时模型ARTM的计算时间仅为120毫秒。因此，本文提出的ARTM适用于快速高精度地预测定位误差。