CN114329837A - 一种数控机床主轴热误差建模及加工精度可靠性分析方法 - Google Patents

Abstract

一种数控机床主轴热误差建模及加工精度可靠性分析方法，属于机床精度可靠性分析领域，其特征在于：基于数控机床热误差产生机理和免疫算法，建立了数控机床的热误差模型；基于多体系统理论，建立了同时包含几何误差和主轴热误差的加工精度模型；基于Latin超立方抽样蒙特卡洛方法，对数控机床加工精度可靠进行了分析。本发明建立的热误差模型具有更高的预测精度，可为数控机床热误差的预测提供方法支持；本发明建立的加工精度模型考虑了多源误差的影响，使得加工精度的预测更接近实际情况；本发明提出的加工精度可靠性分析方法，将机床加工精度可靠性评估从传统的静态范畴扩展到动态范畴，为机床加工精度可靠性分析提供了有效的理论指导。

Description

一种数控机床主轴热误差建模及加工精度可靠性分析方法

技术领域

本发明涉及到一种数控机床主轴热误差建模及加工精度可靠性分析方法，属于机床精度可靠性分析领域。

背景技术

机床作为机械设备生产的技术装备，是国家制造产业的核心设备，其发展水平是国家工业化程度的主要标志之一，随着现代科学技术的迅猛发展，精密超精密加工技术已经成为现代机械制造业发展的主要趋势。数控机床是一种高精度、高效率、高技术的现代机电设备，作为先进制造技术的基础与核心设备，在机械生产中的应用越来越广泛，并制约着制造领域和各高新科技的发展。当前数控机床高速发展的进程中，存在可靠性低、精度保持性差和性能达不到要求的问题，尤其是机床可靠性问题，已成为长期制约数控机床产业发展的关键问题之一。

数控机床的加工精度主要受到几何误差和热误差的影响。机床的几何误差是指由于机床设计、制造、装配等中的缺陷，使得机床中各组成环节或部件的实际几何参数和位置相对于理想几何参数和位置发生偏离。该误差一般与机床各组成环节或部件的几何要素有关，是机床本身固有的误差。机床的热误差是由于机床的电主轴单元在高速旋转的状态下，由于内部旋转部件的摩擦产热，从而导致主轴单元产生的热变形。目前学者所做的大量研究中，主要基于经验公式的热误差模型、温度-热误差模型和多变量-热误差模型。上述热误差建模方法均没有深入到热误差产生机理层面，导致热误差的预测精度较低。因此，建立一个预测精度高的热误差模型是亟待解决的问题。

结构机构可靠性反映了在规定的时间和规定的条件下，结构机构完成规定功能的能力。机床的加工精度可靠性作为评价机床能力的一个重要指标，反映了机床在规定的时间和规定的条件下，完成规定加工精度的能力。大多数学者的研究主要集中在静态范畴上对加工精度可靠性进行分析，其主要考虑的是单一误差源对加工精度可靠性的影响。因此，如何分析机床在几何误差和主轴热误差共同影响下的加工精度可靠性是评估机床能力的关键因素，同时也为机床的设计提供有效的理论支持。

发明内容

本发明的目的是提供一种数控机床主轴热误差建模及加工精度可靠性分析方法。通过对数控机床热误差产生机理进行分析，并基于免疫算法建立机床的热误差模型；基于多体系统理论建立同时包含几何误差和热误差的动态加工精度模型；基于Latin超立方抽样蒙特卡洛方法分析机床的加工精度可靠性。

为了实现上述目的，本发明采用如下技术方案：

一种数控机床主轴热误差建模及加工精度可靠性分析方法，其特征在于：包括如下步骤：

步骤1：建立所述数控机床的热误差模型；

步骤1.1：分析所述数控机床热误差的产生机理；

数控机床的电主轴单元作为机床的核心部件，其包含有许多发热部件。例如主轴单元的前轴承、后轴承和电机。这些发热部件会在主轴单元高速运转时，由于摩擦和电损耗而产生大量的热量，引起主轴单元温度上升从而产生变形。通常情况下，发热部件的周围会设置螺旋式冷却水套，水套中的冷却液用以加速热量的耗散。当电主轴处于工作状态时，其内部的温度梯度使得热量从高温区域向低温区域传递，而与空气接触的部件以热辐射的形式向空气传递热量。

步骤1.1：计算所述主轴单元的轴承发热功率，如公式(1)所示：

W_b＝1.047×10^-4ω·(M_v+M_e) (1)

公式(1)中，ω为轴承的角速度，M_v为与粘性润滑相关的摩擦力矩，M_e为与施加的载荷相关的摩擦力矩。

步骤1.2：计算所述主轴单元的电机发热功率，如公式(2)所示：

公式(2)中，

为电机发热功率，Q₁为电机转子的发热量，Q₂电机定子的发热量，V₁为电机转子的体积，V₂为电机定子的体积。

步骤1.3：计算所述冷却液的换热系数，如公式(3)所示：

公式(3)中，Re为雷诺数，Pr为普朗特数，l_a为主轴单元结构的横截面周长。

步骤1.4：基于免疫算法对所述的发热功率和换热系数进行优化；

建立发热功率和换热系数的优化方程，如公式(4)所示：

公式(4)中，U_{F_B_M}分别为前轴承、后轴承和电机的发热功率，h_{f_n}分别为强制对流和自然对流换热系数，

分别为发热功率和换热系数的比例系数，

分别为发热功率和换热系数的偏差修正系数。

建立优化设计变量的矢量形式，如公式(5)所示：

步骤2：基于多体系统理论建立所述数控机床的加工精度模型；

步骤2.2：建立所述数控机床的特征矩阵；

将所述数控机床抽象为多体系统，该系统主要由8个典型体组成，根据多体系统理论，对应机床各组成部件建立相应的“体”，由D_j(j＝1,2…8)表示。那么可以获得两条分支：床身—X轴运动部件—Y轴运动部件—主轴箱—刀具分支；床身—Z轴运动部件—B轴旋转部件—工件分支。

所述数控机床共有X,Y,Z轴和一个旋转轴B轴，共产生30项几何误差，包括定位误差、直线度误差、角度误差、平行度误差和垂直度误差；共产生3项热误差，包括X，Y，Z方向的热误差。其相邻体间的理想静止、运动特征矩阵和相邻体间的静止、运动误差特征矩阵如表1所示。

表1相邻体间的理想静止、运动特征矩阵和体间静止、运动误差特征矩阵

表1中，p为静止下标，s为运动下标，Δ为相对误差符号，

为体D_i和D_j之间的理想静止特征矩阵，

为体D_i和D_j之间的理想运动特征矩阵，

为体D_i和D_j之间的静止误差特征矩阵；

为D_i和D_j之间的运动误差特征矩阵。

特征矩阵中误差参数分别表示了机床X，Y，Z轴的各部件之间的几何误差和热误差，如表2-表6所示：

表2四轴数控机床X轴6项几何误差和热误差(单位：mm)

表3四轴数控机床Y轴6项几何误差和热误差(单位：mm)

表4四轴数控机床Z轴6项几何误差和热误差(单位：mm)

表5四轴数控机床B轴6项几何误差(单位：mm)

表6四轴数控机床单元间姿态几何误差(单位：mm)

步骤2.3：建立所述数控机床的加工精度模型；

设刀具成形点在刀具坐标系内的坐标为，如公式(6)所示：

则工件成形点在工件坐标系内的坐标为，如公式(7)所示：

H_w＝[H_wx H_wy H_wz 1]^T (7)

公式(6)和公式(7)中，H_tx，H_ty，H_tz分别为刀具成形点在刀具坐标系X轴，Y轴，Z轴上的坐标值，H_wx，H_wy，H_wz分别为工件上成形点在工件坐标系X轴，Y轴，Z轴上的坐标值。

当数控机床在理想情况下空转时，刀具形成点的坐标与工件形成点的坐标将重叠一起。因此，可以得到在理想情况下的加工精度方程，如公式(8)所示：

在得到理想情况下的加工精度方程后，可以将机床的理想刀具形成点由理想工件形成点来表示，如公式(9)所示：

在实际情况中，每个相邻体之间的相对运动产生的几何误差，都会影响加工精度模型的建立。则实际情况下的加工精度方程为，如公式(10)所示：

步骤2.4：所述数控机床的几何误差的测量和主轴热误差的建模；

利用双频激光干涉仪测量数控机床的30项几何误差。

主轴热误差的建模过程为：建立所述的主轴单元的三维有限元模型；将主轴单元的三维模型导入有限元分析软件中进行基本条件的设置；将优化后的发热部件的发热功率和换热系数作为热载荷和边界条件施加给主轴单元；获得主轴单元热稳态温度场；基于主轴单元的温度场得到变形场。

为了验证本发明提出的热误差模型，采用粒子群算法对发热部件的发热功率和冷却液的换热系数进行优化，将优化后的发热功率和换热系数用于主轴单元的热稳态仿真分析中，将实验获得的温度与基于免疫算法和粒子群算法仿真获得的温度进行对比。

步骤2.5：所述数控机床的加工精度的表征；

基于测量获得的30项几何误差和建模获得的热误差，表征同时包含几何误差和主轴热误差的加工精度，将30项几何误差和热误差分别带入公式(10)，均匀选取可以整个工件的30个测量点，表征包含热误差的30个测量点的加工精度和不包含热误差的30个测量点的加工精度，并进行实验获得30个测量点的实测加工精度。

步骤3：所述数控机床的加工精度可靠性分析；

步骤3.1：Latin超立方抽样蒙特卡洛方法分析；

Latin超立方抽样先确定所需的抽样的次数，并对变量的概率进行分层，在累积概率尺度[0,1]上把累积分布曲线分成相等的区间。然后从每个区间或分层中抽取样本，将抽取的样本用作区间的代表值，这些代表值将用来重建变量的概率分布。

用Latin超立方抽样对随机向量V＝(V₁,V₂,...,V_K)^T抽取K个样本v＝(v_i1,v_i2,...,v_in)^T(i＝1,2,...,K)的方法如下：

(1)将每个随机变量V_j(j＝1,2,...,n)的范围分成K个相等的等概率区间，即将变量V_j的累积分布函数

的值域[0,1]等分成K个相互不重叠的子区间[0,1/K],[1/N,2/K],...,[1-1/K,1]。

(2)对每个变量V_j，在它的所有K个子区间内各抽取一个样本。每个区间仅产生一个随机数，作为该区间的代表值。对第i区间，当随机选取其区间代表值u_i时，可先产生U(0,1)的随机数u，则u_i＝(i-1+u)/N；当取作区间中心时，u_i＝(i-1/2)/N。

(3)对变量V_j的K个样本值按其所属区间的序号进行随机排序，然后按变量顺序放一起。这相当于构造一个K行n列抽样序数矩阵R＝[r_ij]_K×n，变量的顺序就是列序，每一列都是序数1,2,…K的随机排列，且没有相同的列，每个变量的K个样本按照各列中的数字排列。

Latin超立方抽样蒙特卡洛方法只需要采用较少的样本就可以反映影响机床加工精度可靠性的随机变量的分布特征，这种特点可以有效地模拟机床的失效可能性。

步骤3.2：加工精度可靠性分析；

当机器结构的某一部分或者整体在超过某一特定状态，不能满足设计规定的某一功能要求和条件时，此特定状态称为结构的极限状态，极限状态是结构工作可靠和不可靠运行的临界状态。加工精度可靠性体现的是机床在规定的功能和条件下，达到机床所规定的加工精度的能力，用加工精度的失效概率来表示。根据结构极限状态的标志，建立功能函数，如公式(11)所示：

η＝g(ζ)＝(ζ₁,ζ₂,...,ζ_n) (11)

在公式(11)中，ζ＝(ζ₁,ζ₂,...,ζ_n)^T为影响加工精度可靠性的随机变量，η＜0表示结构处于失效状态，η＝0表示结构处于临界状态。

函数η是连续随机变量，可以得到结构的失效概率，如公式(12)所示：

在公式(10)中，f_η(z)为随机变量η的概率密度函数，F_η(z)为随机变量η的累积分布函数。

结构的失效概率p_f取决于功能函数η的分布形式，在功能函数η服从正态分布的情况下，其均值为

标准差为

这时可以得到η的概率密度函数f_η(z)，如公式(13)所示：

通过变换

将η转化为标准正态分布，此时可获得其概率密度函数和累积分布函数，如公式(14)所示：

根据公式(14)，正态分布变量η的概率密度函数和累积分布函数可以进行转化，如公式(15)所示：

最终可获得结构的失效概率，如公式(16)所示：

步骤3.3：基于Latin超立方抽样蒙特卡洛方法的加工精度可靠性分析；

基于Latin超立方抽样蒙特卡洛方法的原理和加工精度可靠性的分析过程，计算所述的30测量点的失效概率，并分别使用改进的一次二阶矩方法和蒙特卡洛方法进行加工精度可靠性的分析，用以验证本发明所提出的加工精度可靠性分析方法的有效性。基于Latin超立方抽样蒙特卡洛方法、改进的一次二阶矩方法和蒙特卡洛方法对所述的数控机床的加工精度可靠性的分析过程均在仿真软件中进行。

为了验证本发明提出的加工精度可靠性分析方法，对工件进行实际加工，并选取坐标点，以坐标点的加工精度是否满足要求为依据，计算实际的失效概率。

本发明取得的有益效果如下：

1.本发明提出的一种数控机床主轴热误差建模及加工精度可靠性分析方法，通过分析热误差产生机理，计算发热部件的发热功率和冷却液的换热系数；基于遗传算法对发热功率和换热系数进行优化，优化后的发热功率和换热系数作为热载荷和边界条件应用在主轴单元的热稳态分析来预测其热误差。该模型比现有的经验方法和温度-热误差模型更能揭示主轴热误差的内在机理，且具有更高的热误差预测精度，可以为数控机床热误差的预测提供方法支持。

2.本发明提出的一种数控机床主轴热误差建模及加工精度可靠性分析方法，基于多体系统理论建立机床了的加工精度模型，该模型考虑了几何误差和热误差的综合影响，将加工精度的预测从单一的误差源扩展到多误差源，使得加工精度的预测更接近实际情况。

3.本发明提出的一种数控机床主轴热误差建模及加工精度可靠性分析方法，基于Latin超立方抽样蒙特卡洛方法对机床的加工精度可靠性进行分析，该方法考虑了多个误差源和误差的随机性，将机床加工精度可靠性评估从传统的静态范畴扩展到动态范畴，更接近实际加工过程，为机床加工精度可靠性分析提供了有效的理论指导。

附图说明

图1四轴数控机床主轴单元内部传热过程图；

图2免疫算法的优化过程图；

图3四轴数控机床的整机结构图；

图4四轴数控机床的拓扑结构图；

图5四轴数控机床主轴单元有限元网格结构图；

图6四轴数控机床主轴单元在热稳态下的变形云图：(Ⅰ)主轴单元总体变形，(Ⅱ)主轴单元X方向变形，(Ⅲ)主轴单元Y方向变形，(Ⅳ)主轴单元Z方向变形；

图7四轴数控机床主轴单元的芯轴轴向变形；

图8四轴数控机床30个测量点的加工精度：(Ⅰ)X方向加工精度，(Ⅱ)Y方向加工精度，(Ⅲ)Z方向加工精度；

图9本发明的整体结构框架图；

图10四轴数控机床30个测量点的加工精度可靠性：(Ⅰ)基于蒙特卡洛方法，(Ⅱ)基于Latin超立方抽样蒙特卡洛方法，(Ⅲ)基于改进的一次二阶矩方法。

具体实施方式

通过上面的描述并结合附图说明，本发明会更加清晰，附图说明用于解释本发明方法及实施过程。

本发明以某四轴数控机床为例，对本发明提出的方法进行验证。

具体实施步骤如下：

步骤1：建立所述数控机床的热误差模型；

步骤1.1：分析所述数控机床热误差的产生机理；

数控机床的电主轴单元作为机床的核心部件，其包含有许多发热部件。例如主轴单元的前轴承、后轴承和电机。这些发热部件会在主轴单元高速运转时，由于摩擦和电损耗而产生大量的热量，引起主轴单元温度上升从而产生变形。通常情况下，发热部件的周围会设置螺旋式冷却水套，水套中的冷却液用以加速热量的耗散。当电主轴处于工作状态时，其内部的温度梯度使得热量从高温区域向低温区域传递，而与空气接触的部件以热辐射的形式向空气传递热量。

步骤1.1：计算所述主轴单元的轴承的发热功率，如公式(17)所示。

W_b＝1.047×10^-4ω·(M_v+M_e) (17)

公式(17)中，ω为轴承的角速度，M_v为与粘性润滑相关的摩擦力矩，M_e为与施加的载荷相关的摩擦力矩。

步骤1.2：计算所述主轴单元的电机的发热功率，如公式(18)所示。

公式(18)中，

为电机的发热功率，Q₁为电机转子的发热量，Q₂电机定子的发热量，V₁为电机转子的体积，V₂为电机定子的体积。

步骤1.3：计算所述冷却液的换热系数，如公式(19)所示：

公式(19)中，Re为雷诺数，Pr为普朗特数，l_a为主轴单元结构的横截面周长。

步骤1.4：基于免疫算法对所述的发热功率和换热系数进行优化，如图2所示。

建立发热功率和换热系数的优化方程，如公式(20)所示：

公式(20)中，U_{F_B_M}分别为前轴承、后轴承和电机的发热功率，h_{f_n}分别为强制对流和自然对流换热系数，

分别为发热功率和换热系数的比例系数，

分别为发热功率和换热系数的偏差修正系数。

建立优化设计变量的矢量形式，如公式(21)所示：

步骤2：基于多体系统理论，建立所述数控机床的加工精度模型；

步骤2.2：建立所述数控机床的特征矩阵；

基于多体系统理论建立所述数控机床的加工精度模型，如图3所示。将所述数控机床抽象为多体系统，该系统主要由8个典型体组成，根据多体系统理论，对应机床各组成部件建立相应的“体”，由D_j(j＝1,2…8)表示。那么可以获得两条分支：床身(D₁)—X轴运动部件(D₂)—Y轴运动部件(D₃)—主轴箱(D₄)—刀具(D₅)分支；床身(D₁)—Z轴运动部件(D₆)—B轴旋转部件(D₇)—工件(D₈)分支，如图4所示。

所述数控机床共有X，Y，Z和一个旋转轴B轴，共产生30项几何误差，包括定位误差、直线度误差、角度误差、平行度误差和垂直度误差；共产生3项热误差，包括X，Y，Z方向的热误差。其相邻体间的理想静止、运动特征矩阵和相邻体间的静止、运动误差特征矩阵如表7所示。

表7相邻体间的理想静止、运动特征矩阵和相邻体间的静止、运动误差特征矩阵

表7中，p为静止下标，s为运动下标，Δ为相对误差符号，

为体D_i和D_j之间的理想静止特征矩阵，

为体D_i和D_j之间的理想运动特征矩阵，

为体D_i和D_j之间的静止误差特征矩阵；

为D_i和D_j之间的运动误差特征矩阵。

特征矩阵中误差参数分别表示了机床X，Y，Z轴的各部件之间的几何误差和热误差，如表8-表12：

表8四轴数控机床X轴6项几何误差和热误差(单位：mm)

表9四轴数控机床Y轴6项几何误差和热误差(单位：mm)

表10四轴数控机床Z轴6项几何误差和热误差(单位：mm)

表11四轴数控机床B轴6项几何误差(单位：mm)

表12四轴数控机床单元间姿态几何误差(单位：mm)

步骤2.3：所述四轴数控机床的加工精度建模；

设刀具成形点在刀具坐标系内的坐标为，如公式(22)所示：

则工件成形点在工件坐标系内的坐标为，如公式(23)所示：

H_w＝[H_wx H_wy H_wz 1]^T (23)

公式(4)和公式(5)中，H_tx，H_ty，H_tz分别为刀具成形点在刀具坐标系X轴，Y轴，Z轴上的坐标值，H_wx，H_wy，H_wz分别为工件上成形点在工件坐标系X轴，Y轴，Z轴上的坐标值。

当数控机床以理想情况空转时，刀具形成点的坐标与工件形成点的坐标会重叠一起。因此，可以得到在理想情况下的加工精度方程，如公式(24)所示：

在得到理想情况下的加工精度方程后，可以将机床的理想刀具形成点由理想工件形成点来表示，如公式(25)所示：

在实际情况中，每个相邻体之间的相对运动产生的几何误差，都会影响加工精度模型的建立。则实际情况下的加工精度方程为，如公式(26)所示：

步骤2.4：所述数控机床的几何误差的测量和热误差的建模；

利用双频激光干涉仪测量数控机床的30项几何误差。

热误差的建模过程为：建立所述的主轴单元的三维有限元模型并进行网格划分，如图5所示；将网格模型导入有限元分析软件中进行基本条件的设置；将优化后的发热部件的发热功率和换热系数作为热载荷和边界条件施加给主轴单元；获得主轴单元热稳态温度场；基于主轴单元的温度场得到变形场。如图6和图7所示。

为了验证本发明提出的热误差模型，采用粒子群算法对发热部件的发热功率和冷却液的换热系数进行优化，将优化后的发热功率和换热系数用于主轴单元的热稳态仿真分析中，将实验获得的温度与基于免疫算法和粒子群算法仿真获得的温度进行对比。如表13所示。

表13预测温度与实验温度的比较

由表13中可知：基于免疫算法获得的仿真温度值比基于粒子群算法获得的仿真温度值更接近于实验温度。由此验证了本发明提出的热误差建模方法的优越性。

步骤2.5：所述数控机床的加工精度表征；

基于测量获得的30项几何误差和建模获得的热误差，表征同时包含几何误差和主轴热误差的加工精度，将30项几何误差和热误差分别带入公式(26)，均匀选取可以整个工件的30个测量点，表征包含热误差的30个测量点的加工精度和不包含热误差的30个测量点的加工精度，并进行实验获得30个测量点的实测加工精度，如图8所示。

图8中，包含热误差的加工精度更接近实测加工精度，说明热误差对加工精度具有显著的影响。

步骤3：所述的数控机床的加工精度可靠性分析；

步骤3.1：Latin超立方抽样蒙特卡洛方法分析；

Latin超立方抽样先确定所需的抽样的次数，并对变量的概率进行分层，在累积概率尺度[0,1]上把累积分布曲线分成相等的区间。然后从每个区间或分层中抽取样本，将抽取的样本用作区间的代表值，这些区间的代表值将用来重建变量的概率分布。

用Latin超立方抽样对随机向量V＝(V₁,V₂,...,V_K)^T抽取K个样本v＝(v_i1,v_i2,...,v_in)^T(i＝1,2,...,K)的方法如下：

(1)将每个随机变量V_j(j＝1,2,...,n)的范围分成K个相等的等概率区间，即将变量V_j的累积分布函数

的值域[0,1]等分成K个相互不重叠的子区间[0,1/K],[1/N,2/K],...,[1-1/K,1]。

(2)对每个变量V_j，在它的所有K个子区间内各抽取一个样本。每个区间仅产生一个随机数，作为该区间的代表值。对第i区间，当随机选取其区间代表值u_i时，可先产生U(0,1)的随机数u，则u_i＝(i-1+u)/N；当取作区间中心时，u_i＝(i-1/2)/N。

(3)对变量V_j的K个样本值按其所属区间的序号进行随机排序，然后按变量顺序放一起。这相当于构造一个K行n列抽样序数矩阵R＝[r_ij]_K×n，变量的顺序就是列序，每一列都是序数1,2,…K的随机排列，且没有相同的列，每个变量的K个样本按照各列中的数字排列。

Latin超立方抽样蒙特卡洛方法只需要采用较少的抽样就可以反映了影响机床加工精度可靠性的随机变量的分布特征，这种特点可以有效地模拟机床的失效可能性。

步骤3.2：加工精度可靠性分析；

当机器结构的某一部分或者整体在超过某一特定状态，不能满足设计规定的某一功能要求和条件时，此特定状态称为结构的极限状态，极限状态是结构工作可靠和不可靠运行的临界状态。加工精度可靠性体现的是机床在规定的功能和条件下，达到机床所规定的加工精度的能力，用加工精度的失效概率来表示。根据结构极限状态的标志，建立功能函数，如公式(27)所示：

η＝g(ζ)＝(ζ₁,ζ₂,...,ζ_n) (27)

在公式(27)中，ζ＝(ζ₁,ζ₂,...,ζ_n)^T为影响加工精度可靠性的随机变量，η＜0表示结构处于失效状态，η＝0表示结构处于临界状态。

函数η是连续随机变量，可以得到结构的失效概率，如公式(28)所示：

在公式(28)中，f_η(z)为随机变量η的概率密度函数，F_η(z)为随机变量η的累积分布函数。

结构的失效概率p_f取决于功能函数η的分布形式，在功能函数η服从正态分布的情况下，其均值为

标准差为

这时可以得到η的概率密度函数f_η(z)，如公式(29)所示：

通过变换

将η转化为标准正态分布，此时可获得其概率密度函数和累积分布函数，如公式(30)所示，

根据公式(30)，正态分布变量η的概率密度函数和累积分布函数可以进行转化，如公式(31)所示：

最终可获得结构的失效概率，如公式(32)所示：

本发明内容的总体框架如图9所示。

步骤3.3：基于Latin超立方抽样蒙特卡洛方法的加工精度可靠性分析；

基于Latin超立方抽样蒙特卡洛方法的原理和加工精度可靠性的分析过程，计算所述的数控机床30测量点的失效概率，并分别使用改进的一次二阶矩方法和蒙塔卡洛方法进行加工精度可靠性的分析，用以验证本发明所提出的加工精度可靠性分析方法的有效性，基于Latin超立方抽样蒙特卡洛方法、改进的一次二阶矩方法和蒙塔卡洛方法对所述的数控机床的加工精度可靠性的分析过程在仿真软件中进行，分析结果如图10所示。

图10中，基于Latin超立方抽样蒙特卡洛方法的加工精度可靠性预测值偏离蒙特卡洛方法的程度比改进的一次二阶矩方法更小，说明本发明提出的方法的加工精度可靠性预测精度更高。

为了验证本发明提出的加工精度可靠性分析方法，对工件进行实际加工，并选取坐标点，以坐标点的加工精度是否满足要求为依据，计算实际的失效概率，如表14所示。

表14不同方法下的测量点的失效概率

表14中，在不同的坐标下分别加工300个工件，可以看出，基于Latin超立方抽样蒙特卡洛法预测的不同坐标下的加工精度可靠性与一次二阶矩法预测的不同坐标下的加工精度可靠性相比具有更高的精度，验证了本发明提出的加工精度可靠性分析方法的有效性。

以上所述仅是本发明优选的实施方式，但本发明的保护范围并不局限于此，任何基于本发明所提供的技术方案和发明构思进行的改造和替换都应包含在本发明的保护范围内。

Claims (1)

1.一种数控机床主轴热误差建模及加工精度可靠性分析方法，其特征在于：包括如下步骤：

步骤1：建立所述数控机床的热误差模型；

步骤1.1：分析所述数控机床热误差的产生机理；

数控机床的电主轴单元作为机床的核心部件，其包含有许多发热部件。例如主轴单元的前轴承、后轴承和电机。这些发热部件会在主轴单元高速运转时，由于摩擦和电损耗而产生大量的热量，引起主轴单元温度上升从而产生变形。通常情况下，发热部件的周围会设置螺旋式冷却水套，水套中的冷却液用以加速热量的耗散。当电主轴处于工作状态时，其内部的温度梯度使得热量从高温区域向低温区域传递，而与空气接触的部件以热辐射的形式向空气传递热量。

步骤1.1：计算所述主轴单元的轴承发热功率，如公式(1)所示：

W_b＝1.047×10^-4ω·(M_v+M_e) (1)

公式(1)中，ω为轴承的角速度，M_v为与粘性润滑相关的摩擦力矩，M_e为与施加的载荷相关的摩擦力矩。

步骤1.2：计算所述主轴单元的电机发热功率，如公式(2)所示：

公式(2)中，

为电机发热功率，Q₁为电机转子的发热量，Q₂电机定子的发热量，V₁为电机转子的体积，V₂为电机定子的体积。

步骤1.3：计算所述冷却液的换热系数，如公式(3)所示：

公式(3)中，Re为雷诺数，Pr为普朗特数，l_a为主轴单元结构的横截面周长。

步骤1.4：基于免疫算法对所述的发热功率和换热系数进行优化；

建立发热功率和换热系数的优化方程，如公式(4)所示：

公式(4)中，U_{F_B_M}分别为前轴承、后轴承和电机的发热功率，h_{f_n}分别为强制对流和自然对流换热系数，

分别为发热功率和换热系数的比例系数，

分别为发热功率和换热系数的偏差修正系数。

建立优化设计变量的矢量形式，如公式(5)所示：

步骤2：基于多体系统理论建立所述数控机床的加工精度模型；

步骤2.2：建立所述数控机床的特征矩阵；

将所述数控机床抽象为多体系统，该系统主要由8个典型体组成，根据多体系统理论，对应机床各组成部件建立相应的“体”，由D_j(j＝1,2…8)表示。那么可以获得两条分支：床身—X轴运动部件—Y轴运动部件—主轴箱—刀具分支；床身—Z轴运动部件—B轴旋转部件—工件分支。

所述数控机床共有X,Y,Z轴和一个旋转轴B轴，共产生30项几何误差，包括定位误差、直线度误差、角度误差、平行度误差和垂直度误差；共产生3项热误差，包括X，Y，Z方向的热误差。其相邻体间的理想静止、运动特征矩阵和相邻体间的静止、运动误差特征矩阵如表1所示。

表1 相邻体间的理想静止、运动特征矩阵和体间静止、运动误差特征矩阵

表1中，p为静止下标，s为运动下标，Δ为相对误差符号，

为体D_i和D_j之间的理想静止特征矩阵，

为体D_i和D_j之间的理想运动特征矩阵，

为体D_i和D_j之间的静止误差特征矩阵；

为D_i和D_j之间的运动误差特征矩阵。

特征矩阵中误差参数分别表示了机床X，Y，Z轴的各部件之间的几何误差和热误差，如表2-表6所示：

表2 四轴数控机床X轴6项几何误差和热误差(单位：mm)

表3 四轴数控机床Y轴6项几何误差和热误差(单位：mm)

表4 四轴数控机床Z轴6项几何误差和热误差(单位：mm)

表5 四轴数控机床B轴6项几何误差(单位：mm)

表6 四轴数控机床单元间姿态几何误差(单位：mm)

步骤2.3：建立所述数控机床的加工精度模型；

设刀具成形点在刀具坐标系内的坐标为，如公式(6)所示：

则工件成形点在工件坐标系内的坐标为，如公式(7)所示：

H_w＝[H_wx H_wy H_wz 1]^T (7)

公式(6)和公式(7)中，H_tx，H_ty，H_tz分别为刀具成形点在刀具坐标系X轴，Y轴，Z轴上的坐标值，H_wx，H_wy，H_wz分别为工件上成形点在工件坐标系X轴，Y轴，Z轴上的坐标值。

当数控机床在理想情况下空转时，刀具形成点的坐标与工件形成点的坐标将重叠一起。因此，可以得到在理想情况下的加工精度方程，如公式(8)所示：

在得到理想情况下的加工精度方程后，可以将机床的理想刀具形成点由理想工件形成点来表示，如公式(9)所示：

在实际情况中，每个相邻体之间的相对运动产生的几何误差，都会影响加工精度模型的建立。则实际情况下的加工精度方程为，如公式(10)所示：

步骤2.4：所述数控机床的几何误差的测量和主轴热误差的建模；

利用双频激光干涉仪测量数控机床的30项几何误差。

主轴热误差的建模过程为：建立所述的主轴单元的三维有限元模型；将主轴单元的三维模型导入有限元分析软件中进行基本条件的设置；将优化后的发热部件的发热功率和换热系数作为热载荷和边界条件施加给主轴单元；获得主轴单元热稳态温度场；基于主轴单元的温度场得到变形场。

为了验证本发明提出的热误差模型，采用粒子群算法对发热部件的发热功率和冷却液的换热系数进行优化，将优化后的发热功率和换热系数用于主轴单元的热稳态仿真分析中，将实验获得的温度与基于免疫算法和粒子群算法仿真获得的温度进行对比。

步骤2.5：所述数控机床的加工精度的表征；

基于测量获得的30项几何误差和建模获得的热误差，表征同时包含几何误差和主轴热误差的加工精度，将30项几何误差和热误差分别带入公式(10)，均匀选取可以整个工件的30个测量点，表征包含热误差的30个测量点的加工精度和不包含热误差的30个测量点的加工精度，并进行实验获得30个测量点的实测加工精度。

步骤3：所述数控机床的加工精度可靠性分析；

步骤3.1：Latin超立方抽样蒙特卡洛方法分析；

Latin超立方抽样先确定所需的抽样的次数，并对变量的概率进行分层，在累积概率尺度[0,1]上把累积分布曲线分成相等的区间。然后从每个区间或分层中抽取样本，将抽取的样本用作区间的代表值，这些代表值将用来重建变量的概率分布。

用Latin超立方抽样对随机向量V＝(V₁,V₂,...,V_K)^T抽取K个样本v＝(v_i1,v_i2,...,v_in)^T(i＝1,2,...,K)的方法如下：

(1)将每个随机变量V_j(j＝1,2,...,n)的范围分成K个相等的等概率区间，即将变量V_j的累积分布函数

的值域[0,1]等分成K个相互不重叠的子区间[0，1/K]，[1/N，2/K]，...，[1-1/K，1]。

(2)对每个变量V_j，在它的所有K个子区间内各抽取一个样本。每个区间仅产生一个随机数，作为该区间的代表值。对第i区间，当随机选取其区间代表值u_i时，可先产生U(0,1)的随机数u，则u_i＝(i-1+u)/N；当取作区间中心时，u_i＝(i-1/2)/N。

(3)对变量V_j的K个样本值按其所属区间的序号进行随机排序，然后按变量顺序放一起。这相当于构造一个K行n列抽样序数矩阵R＝[r_ij]_K×n，变量的顺序就是列序，每一列都是序数1,2,…K的随机排列，且没有相同的列，每个变量的K个样本按照各列中的数字排列。

Latin超立方抽样蒙特卡洛方法只需要采用较少的样本就可以反映影响机床加工精度可靠性的随机变量的分布特征，这种特点可以有效地模拟机床的失效可能性。

步骤3.2：加工精度可靠性分析；

当机器结构的某一部分或者整体在超过某一特定状态，不能满足设计规定的某一功能要求和条件时，此特定状态称为结构的极限状态，极限状态是结构工作可靠和不可靠运行的临界状态。加工精度可靠性体现的是机床在规定的功能和条件下，达到机床所规定的加工精度的能力，用加工精度的失效概率来表示。根据结构极限状态的标志，建立功能函数，如公式(11)所示：

η＝g(ζ)＝(ζ₁,ζ₂,...,ζ_n) (11)

在公式(11)中，ζ＝(ζ₁,ζ₂,...,ζ_n)^T为影响加工精度可靠性的随机变量，η＜0表示结构处于失效状态，η＝0表示结构处于临界状态。

函数η是连续随机变量，可以得到结构的失效概率，如公式(12)所示：

在公式(10)中，f_η(z)为随机变量η的概率密度函数，F_η(z)为随机变量η的累积分布函数。

结构的失效概率p_f取决于功能函数η的分布形式，在功能函数η服从正态分布的情况下，其均值为

标准差为

这时可以得到η的概率密度函数f_η(z)，如公式(13)所示：

通过变换

将η转化为标准正态分布，此时可获得其概率密度函数和累积分布函数，如公式(14)所示：

根据公式(14)，正态分布变量η的概率密度函数和累积分布函数可以进行转化，如公式(15)所示：

最终可获得结构的失效概率，如公式(16)所示：

步骤3.3：基于Latin超立方抽样蒙特卡洛方法的加工精度可靠性分析；

基于Latin超立方抽样蒙特卡洛方法的原理和加工精度可靠性的分析过程，计算所述的30测量点的失效概率，并分别使用改进的一次二阶矩方法和蒙特卡洛方法进行加工精度可靠性的分析，用以验证本发明所提出的加工精度可靠性分析方法的有效性。基于Latin超立方抽样蒙特卡洛方法、改进的一次二阶矩方法和蒙特卡洛方法对所述的数控机床的加工精度可靠性的分析过程均在仿真软件中进行。

为了验证本发明提出的加工精度可靠性分析方法，对工件进行实际加工，并选取坐标点，以坐标点的加工精度是否满足要求为依据，计算实际的失效概率。

Images