李 猛 王时龙 马 驰
夏长久重庆大学机械传动国家重点实验室，重庆,400044
大规格齿轮是航母等大型舰船、大功率风力发电机和重载工程机械等高端设备上的关键核心基础件[1]，其性能优劣直接决定装备传动系统的振动、噪声、寿命等服役性能和核心竞争力。滚齿是一种高效、低成本制齿工艺，对应的加工机床在齿轮生产制造中占据重要地位[2]。相较于中小规格机床，大规格滚齿机床的工件齿坯质量大、惯性大，X、Z轴和工作台分别采用大平面矩形镶钢导轨和多腔同步控制的静压导轨副以提高刚性。大规格滚齿机床的加工精度受机床多源误差影响，主要包括几何误差、热变形误差、切削力致误差、运动轴伺服误差等，其中几何误差对机床加工精度的影响较大且易测量、重复性好、系统性高、长时间稳定[3]，易于实施补偿，所以几何误差是机床误差补偿的重要对象。然而，大规格滚齿机床床身、立柱、工作台尺寸大，几何误差对工件精度影响更为显著且误差数量众多、耦合复杂，因此，建立大规格滚齿机床几何误差敏感性分析模型并识别出关键误差对提高补偿效率、改进机床设计具有指导意义。
几何误差模型是关键误差识别的基础，滚齿机床误差建模主要是借鉴三轴或五轴机床的建模方法，主要有多体系统理论[4-5]、齐次坐标变换[6]、旋量理论[7-8]等方法，其中多体系统理论结合齐次坐标变换(HTM)方法应用最为广泛[9]。CONG等[10]采用HTM方法建立了五轴机床的空间误差模型并预测了机床加工误差。FU等[11]基于旋量理论和运动关系建立了三轴立式加工中心的几何误差模型并通过实验验证其精度。
现有的机床关键误差识别方法主要是在机床几何误差模型的基础上，结合相关敏感性分析方法，对该模型进行全局敏感性分析，量化几何误差对加工误差的影响权重，从而识别出关键几何误差。程强等[12]基于四轴精密加工中心几何误差模型，利用矩阵微分法计算了37项几何误差敏感度系数，并识别出关键性几何误差因素。范晋伟等[13]采用基于方差的全局敏感性分析方法分析了龙门加工中心37项几何误差敏感度系数，得到了影响加工误差的关键因素。胡腾等[14]基于卧式加工中心空间误差完备模型和实际参与度，采用矩阵微分法甄别出关键误差。付国强等[15]基于误差敏感矩阵建立了五轴机床几何误差贡献值模型。杨赟等[16]基于齐次坐标建立了立式加工中心的误差模型，对该模型的21项几何误差求偏导，得到了与位置相关的误差敏感度系数矩阵，忽略次要误差，简化了误差补偿模型。CHENG等[17]基于指数积旋量理论建立了三轴立式加工中心几何误差模型，采用Morris法识别出关键几何误差。ZOU等[18]采用Sobol法对超精密金刚石车床各误差项敏感度指标进行量化，研究了各误差项对机床加工精度的影响关系。GUO等[19]考虑机床工作空间中几何误差矢量的多维输出并进行全局定量敏感性分析，确定了几何误差影响因子。FAN等[20]基于多体系统理论建立了五轴机床误差模型，将误差参数输出值表示为截断傅里叶级数的幅值，将全局敏感度指数定义为幅值方差与函数方差比值，确定了关键几何误差。GUO等[21]采用乘法降维法建立了机床几何误差一维积分模型，实现了误差敏感性分析。WANG等[22]基于微分运动关系建立了龙门导轨磨床的误差模型，采用改进的一次二阶矩法对机床加工精度进行了敏感性分析。以上关于机床几何误差敏感性分析的研究所采用的空间误差模型由于铣刀、车刀等刀具结构简单而将其简化成一点，忽略了刀具与工件间运动干涉的影响，即用刀具位姿误差表征了机床空间加工误差。同时机床运动轴的数量少，误差间的耦合作用不明显，可以忽略不计。
对于滚齿机，滚刀形状结构复杂，齿轮由滚刀多个切削刃断续切削包络而成，若用某个切削点误差来表征滚齿机床加工误差，并基于此误差模型来计算误差敏感度系数，其识别结果往往不够准确。另外，滚齿机床运动链长、误差传递关系复杂，几何误差间的耦合作用明显。因此，本文在机床刀具位姿-几何误差模型的基础上，考虑滚刀与齿轮工件的啮合原理，基于双参数包络理论，分别求解了理想状态下和误差作用下加工齿面接触迹点的坐标值，从而建立了机床齿面位姿-几何误差模型；
然后，采用Sobol法对机床齿面位姿误差模型进行全局敏感性分析，量化和分析了机床几何误差敏感度系数，并识别出关键几何误差项；
最后，通过关键误差的虚拟仿真修正并与矩阵微分法、Morris法以及基于刀具位姿误差模型识别方法的对比分析，验证了所提方法的准确性及优势。
1.1 几何误差元素分析
大规格数控滚齿机包含3个直线轴(X轴、Y轴、Z轴)和3个旋转轴(A轴、B轴、C轴)，其示意图和拓扑结构见图1。机床通过控制运动轴间联动使滚刀与工件相对运动从而实现齿轮包络成形。由于机床零部件制造和装配时存在误差，机床各轴运动偏离理想运动轨迹而产生几何误差。几何误差通常分为位置相关几何误差(PDGEs)和位置无关几何误差(PIGEs)。位置相关几何误差主要是由机床运动部件自身的制造误差引起的，它在运动过程中是因实际运动位置变化而改变的六自由度误差。以机床X轴为例，当机床沿该轴移动时会产生6项几何误差，分别是定位误差、Y向直线度误差、Z向直线度误差、滚转误差、俯仰误差和偏摆误差，如图2所示，其中δ、ε分别为位置误差和姿态误差。机床其他轴均类似，共有36项位置相关几何误差。位置无关几何误差主要是由机床装配过程中的缺陷引起的，不受运动轴位置的影响，数控滚齿机床共有15项位置无关几何误差，其定义如图3所示，其中S为垂直度误差。大规格数控滚齿机共有51项几何误差，如表1所示。
0.床身 1.C轴 2.齿轮工件 3.X轴 4.Z轴 5.A轴 6.Y轴 7.B轴 8.滚刀图1 大规格滚齿机及其拓扑结构Fig.1 Large size gear hobbing machine and its topology
图2 X轴位置相关几何误差Fig.2 The position-dependent geometric errors of X-axis
(a)垂直度误差 (b)C轴
表1 大规格滚齿机床几何误差项及编号
1.2 刀具位姿-几何误差建模
iT2,8=(T1,2T0,1)-1T0,3T3,4T4,5T5,6T6,7T7,8
(1)
当存在几何误差时，两相邻体间的实际位姿变换矩阵可表示为
(2)
实际情况下，滚刀坐标系到工件坐标系的变换矩阵为
eT2,8=(eT1,2eT0,1)-1eT0,3eT3,4eT4,5eT5,6eT6,7eT7,8
(3)
设(a,b,c,1)T、(i,j,k,0)T分别为滚刀坐标系下滚刀中心位置和姿态的齐次坐标，若δt、εt分别为刀具的位置误差和姿态误差，则机床刀具位姿-几何误差模型可表示为
(4)
1.3 齿面位姿-几何误差建模
齿轮滚齿加工是滚刀切削刃基于展成原理形成的一系列空间轨迹曲面断续切除圆柱齿坯材料，最终离散包络出齿轮齿面。为了便于分析计算，本文滚刀采用渐开线滚刀，如图4所示，其端面齿廓的参数表达式为
(5)
表2 相邻体间的位姿变换
图4 渐开线滚刀端面齿廓Fig.4 End-face profile of the involute hob
式中，φh为渐开线展角；
μh为滚刀齿廓夹角半角；
rbh为滚刀基圆半径。
令滚刀端面齿廓绕滚刀轴线做螺旋运动，它在空间扫掠形成的轨迹曲面就是滚刀的渐开螺旋齿面，其参数表达式为
rh(φh,θh)=Mhprp
(6)
(7)
式中，p为滚刀螺旋参数；
θh为回转参数。
滚刀齿面单位法向量为
(8)
基于滚齿机床运动链并考虑机床几何误差位姿变换矩阵，可以推导出滚刀处于工件坐标系中分别在理想情况和实际情况下的空间成形曲面为
(9)
(10)
根据双参数包络理论[23]，滚刀齿面在切削运动中形成包络面需要满足以下方程组：
(11)
滚齿加工时，X轴移动距离x确定，根据滚齿运动关系[24]，可以得到滚刀转角与工件齿轮转角关系，将滚刀的廓形沿其轴向离散便可得到一组参数(φhi,θhi)，然后将已知量代入式中求得滚齿啮合点坐标，获得理想情况下齿轮齿面iRg和齿面法矢iNg。当引入机床几何误差时，获得实际情况下的齿轮齿面eRg和齿面法矢eNg。因此，大规格滚齿机床齿面位姿-几何误差模型可表示为
(12)
齿面位姿误差定义为E=(δg;εg)=(δX,δY,δZ,εX,εY,εZ)T。其中，δX、δY、δZ分别为齿面点集在工件坐标系沿X、Y、Z方向的位置误差分量，描述齿面位置偏移大小；
εX、εY、εZ分别为齿面点集的单位法矢在工件坐标系沿X、Y、Z方向的误差分量，描述齿面滚转、俯仰和偏摆误差大小。
2.1 Sobol法
Sobol法主要思想是：将机床误差模型分解成递增阶数项之和，各输入量函数值方差与误差模型方差的比值作为各输入量的敏感度系数，从而得到各误差项的一阶和总体敏感度系数。机床51项几何误差为输入量，6项齿面位姿误差为输出量，机床误差模型可表示为
E=f(G)
(13)
G=(x1,x2,…,xn)T
其中，xn为几何误差项，取n=51。
由于模型有多个输出量且为标量，故可以将多输出模型分解为多个独立的单输出模型进行分析。以齿面X方向的位置误差为例，δX与机床几何误差的映射模型为
δX=f(G)
(14)
假设输入变量的空间域为一个n维的单元体I(n)=(xi|0≤xi≤1；
i=1,2,…,n)，根据Sobol法的基本原理，式(14)的分解模型为
(15)
其中，f0为基于输入参数得到f的期望值，fi(xi)表示关于xi的函数，fij(xi,xj)表示关于xi和xj的函数，其他高阶子项以此类推。
当各变量相互独立且正交时，模型f(x)的总方差和各阶子项的偏方差可分别表示为
(16)
(17)
对式(15)两边求方差后得到：
(18)
Si=Vi/VSij=Vij/VS12…n=V12…n/V
其中，Si为输入量xi的一阶敏感度系数，表示误差项xi对齿面位姿误差的影响权重，Si值越大说明影响程度越大；
Sij为输入量xi和xj的二阶敏感度系数，表示误差项xi和xj耦合作用对齿面位姿误差的影响权重，Sij值越大说明耦合作用越强；
其他高阶项以此类推。
由于误差项之间互相耦合，为衡量误差项xi自身及与其他输入误差项的耦合作用对输出方差的影响权重，引入总体敏感度系数STi，其表达式为
(19)
式中，x～i为除去xi后的所有误差项；
Exi、Vxi分别为计算xi得到的期望和方差；
Ex～i、Vx～i分别为计算x～i得到的期望和方差。
2.2 Quasi-Monte Carlo估算
为了获得上述敏感度系数，在方差计算中需要对输入参数进行采样来估计多维积分。常规蒙特卡罗采样是采用计算机内部随机函数生成的伪随机序列，不可避免地出现周期性重复。Sobol法采用偏差较小、均匀分布的拟随机序列来代替伪随机序列，从而获得更好的收敛性，提高了计算速度。采用两种方法在0.02×0.02范围内随机生成5000个点，其分布情况如图5所示。具体计算过程如下：
(1)根据Sobol法对误差项进行准蒙特卡罗采样，生成一个m×2n的随机样本矩阵D，其中m为误差项的采样个数，n为误差项的个数。
(a)伪随机序列
(20)
以Y31600CNC6型大规格滚齿机床为例进行几何误差敏感性分析。由于几何误差的取值范围对敏感性分析结果的影响较大，因此必须提前确定其概率分布。位置相关几何误差与机床运动轴位置密切相关，可通过测量工作空间中离散点的误差值来估计。使用Renishaw XL-80激光干涉仪和Renishaw QC10球杆仪对机床运动轴误差进行多次测量，测得位置误差和角度误差分别大致分布在[0,40]μm和[0,0.000 53]rad范围内。为了更准确地量化机床几何误差对齿面位姿误差的影响，需要修正原Sobol序列，即在原Sobol序列与实测几何误差之间作代数运算。采样数m影响分析计算效率和精度，经过多次重复试验，取m=100。当加工基本参数如表3所示时，机床X、Y、A轴的位置分别为616.111 mm、125 mm、9.97°，在Z轴轴向进给方向810～890 mm范围内等距离取5个位置，分别建立机床齿面位姿-几何误差模型并采用Sobol法对该模型进行敏感性分析。通过求取5个位置齿面接触迹点集的平均误差敏感度系数来代替齿轮整个齿面的误差敏感度系数。用横坐标代表几何误差序号、纵坐标代表敏感度系数，机床几何误差的敏感性分析结果如图6～图11所示。
表3 滚刀和齿轮参数
图6 齿面位姿误差分量δX敏感度系数Fig.6 The sensitivity coefficient of tooth surfaceposture errors in the δX-direction
图7 齿面位姿误差分量δY敏感度系数Fig.7 The sensitivity coefficient of tooth surfaceposture errors in the δY-direction
图8 齿面位姿误差分量δZ敏感度系数Fig.8 The sensitivity coefficient of tooth surfaceposture errors in the δZ-direction
图9 齿面位姿误差分量εX敏感度系数Fig.9 The sensitivity coefficient of tooth surfaceposture errors in the εX-direction
图10 齿面位姿误差分量εY敏感度系数Fig.10 The sensitivity coefficient of tooth surfaceposture errors in the εY-direction
图11 齿面位姿误差分量εZ敏感度系数Fig.11 The sensitivity coefficient of tooth surfaceposture errors in the εZ-direction
根据图6～图11，比较各误差项的一阶敏感度系数Si和总体敏感度系数STi，可识别出对该齿面分量影响较大的关键误差。当几何误差项的一阶敏感度系数和总体敏感度系数相差不大时，说明该误差项与其他误差的耦合作用较弱。几何误差项的一阶敏感度系数和总体敏感度系数的大小关系并不是确定的，当误差项的总体敏感度系数大于一阶敏感度系数时，误差项与其他误差的耦合作用增强了该误差项对齿面位姿误差的影响程度，当总体敏感度系数小于一阶敏感度系数时，误差项与其他误差的耦合作用减弱了该误差项对齿面位姿误差的影响程度，但误差项间的耦合作用机制并不明确，说明误差项间的耦合作用复杂。同时可通过计算运动轴多项误差的敏感度系数之和来评估敏感部件。机床51项几何误差的敏感度系数之和为1，若各项误差敏感度系数相等，则各系数约为0.02，本文以其1.5倍即0.03为阈值来评估关键误差并且将一阶敏感度系数和总体敏感度系数大于0.03的误差项均视为关键误差。以齿面位姿误差分量δX为例，对比分析各误差项敏感度系数，识别出的关键误差项为δX(X)、εY(X)、δX(Y)、δX(Z)、SZX、δX(A)、δX(C)、εY(C)、δCX、εCY、δX(B)、δBX，强耦合误差为εY(X)、εCY，敏感部件为C轴、X轴。由此可知，若减小该齿面误差，上述关键误差应重点补偿与修正。同时，也可合理提高敏感部件C轴、X轴的制造精度。同理，其他齿面位姿误差分量识别结果如表4所示。
表4 各齿面位姿误差分量的关键误差项
由敏感性分析结果可知，C轴的几何误差对六个齿面位姿误差分量都有显著影响，而Y轴和A轴的几何误差对齿面姿态误差影响较大。角度误差对齿面位姿误差的影响远大于直线度误差和定位误差的影响。由表4可知，机床的敏感部件为C、X、Y、A轴，这些运动轴是精度设计和误差补偿的重点对象。
4.1 误差识别结果仿真验证
为验证机床误差敏感性分析识别结果的准确性，对关键几何误差项进行虚拟仿真修正，即针对某一齿面位姿误差分量，将影响它的关键几何误差项修正为0，其余误差项仍为Sobol序列随机采样的数值，代入齿面位姿-几何误差模型中得到修正后的齿面位姿误差，计算关键误差修正前后各齿面位姿误差分量减小百分比并将其定义为误差消减率。其中，主误差消减率为被修正齿面位姿误差分量的误差消减率，其数值越大，表明关键误差影响权重越大，识别效果越好。经过上述修正后，所得到的误差消减率计算结果如表5所示。
表5 误差消减率
由表5可知，各齿面位姿误差分量的主误差消减率均在60%以上，最大达到98.39%,说明关键几何误差对齿面精度的影响显著，远超其他几何误差的影响，证明识别结果准确可靠。但对齿面位姿误差分量δZ的关键误差项εY(C)、εCY修正后，发现它们同样是齿面位姿误差分量δX的关键误差，而δX却反向增大，说明关键误差修正后可能与其他误差耦合导致其他齿面位姿误差分量增大，这也表现出机床几何误差耦合作用的复杂性。因此，在后续的误差补偿中，应综合考虑几何误差之间的耦合作用，从而提高机床加工精度。
4.2 关键几何误差项识别方法对比
为证明Sobol方法识别效果的优越性，以齿面位姿误差分量δX的一阶敏感度系数为参照，与基于矩阵微分法和Morris法的分析结果对比，结果如图12所示。矩阵微分法仅识别出其中4项关键误差，效果较差且无法评估误差之间的耦合强度。Morris法与Sobol法的关键误差识别结果基本一致，但是Morris法只能用作定性分析而Sobol法可以定量分析关键误差的敏感效应，因此Sobol法更加适用于机床几何误差的敏感性分析。
图12 不同方法齿面位姿误差分量δX敏感性分析结果Fig.12 Comparison of sensitivity analysis results oftooth surface posture errors in the δX-dirctionby different methods
现有机床误差敏感性分析基本是针对三轴机床以刀具位姿-几何误差模型为基础识别出关键误差。滚齿机床刀具位姿-几何误差空间模型通过机床几何误差影响滚刀轴线的位置和方向,进而改变切削点空间位姿，但忽略了滚刀复杂几何形状及其与工件的包络成形过程，识别结果往往不够准确。为进一步研究基于齿面位姿-几何误差模型敏感性分析结果的优势，以一阶敏感度系数Si为对比参数分析比较基于齿面位姿-几何误差模型识别结果与基于刀具位姿-几何误差模型识别结果，结果如图13所示。
(a)齿面位姿误差分量δX分析结果比较
由图13可知，针对误差分量δX，采用两种机床误差模型识别出的关键几何误差项基本相同，仅幅值有些差别，证明了齿面位姿-几何误差模型识别结果的可靠性。针对误差分量δY和δZ，基于齿面位姿-几何误差模型的识别方法比基于刀具位姿-几何误差模型的方法能识别出更多关键几何误差项，对误差分量δY的两种模型识别结果进行关键误差修正仿真，前者主误差消减率为68.63%，后者主误差消减率为36.71%，说明基于齿面位姿-几何误差模型的识别结果更加精准可靠。
(1)基于多体动力学和滚齿双参数包络原理，建立了大规格滚齿机床刀具位姿-几何误差模型和齿面位姿-几何误差模型以表明机床几何误差与齿面位姿误差之间的映射关系。
(2)提出了基于机床齿面位姿-几何误差模型和Sobol全局敏感性分析的机床关键误差识别方法，计算出机床各几何误差项的敏感度系数，从而识别出关键几何误差项和机床敏感部件，为后续误差补偿与精度设计提供了理论支持。
(3)通过关键误差修正仿真实验和对比分析实验，验证了识别结果的准确性和优越性。与Morris法和矩阵微分法相比，Sobol法能更准确量化几何误差的敏感系数；
与刀具位姿-几何误差模型相比，基于齿面位姿-几何误差模型的识别方法对齿面位姿误差分量δY和δZ的识别效果更好。
