针对双转台五轴数控机床旋转轴几何误差辨识精度问题，提出一种在机床坐标系下的旋转轴位置无关几何误差(PIGEs)辨识模型。建立了刀具球与工件球在机床坐标系下的实际初始坐标，通过逆矩阵得到工件球在所测轴的实际初始位置，基于4种测量模式建立了包含球杆仪安装误差与旋转轴PIGEs的杆长变化量数学模型。仿真分析了PIGEs对测量模式的影响，结果表明刀具球在两旋转轴轴线交点位置时，平行度误差不影响球杆仪杆长变化量。最后，通过实验辨识出旋转轴8项PIGEs，并对旋转轴PIGEs的4项位置误差进行补偿。实验结果表明，补偿后的位置误差最大绝对值由203.5μm减小至5.1μm，所提出的辨识模型可以有效提高五轴机床精度。
关键词：五轴数控机床;旋转轴;位置无关几何误差;机床坐标系;球杆仪
论文《双转台五轴数控机床旋转轴位置无关几何误差辨识》发表在《中国机械工程》，版权归《中国机械工程》所有。本文来自网络平台，仅供参考。
0 引言
五轴机床在加工形状复杂的零部件时有独特的优势，提高五轴机床的精度是提高加工精度的关键[1]。在加工过程中，影响加工精度的因素包括几何误差[2]、热误差[3]、伺服误差[4]和切削力[5]引起的误差。根据标准ISO230-1[6]、ISO230-7[7]和ISO10791-1[8]中的定义，机床几何误差是加工过程中影响加工精度的主要来源。目前已知五轴机床的几何误差多达41项，其中旋转轴的几何误差成为加工过程中的主要误差源，包括12项位置相关几何误差(position-dependent geometric errors, PDGEs)以及8项位置无关几何误差(position-independent geometric errors, PIGEs)[9]。PDGEs与旋转轴本身的制造缺陷有关，而PIGEs主要是旋转轴组装时引起的装配误差[10]。在零件加工过程中，PDGEs与运动轴的空间位置有关，因此其大小随着被控轴的控制命令改变而改变。而PIGEs具有恒定的特性，即不会随被控轴的控制命令改变而发生变化。所以对于五轴机床，PIGEs是决定机床精度的关键因素[11]，因此针对旋转轴的PIGEs开展测量与辨识研究尤为重要。
为了测量和辨识PIGEs，学者们做了很多工作，提出了各种识别PIGEs的数学模型和测量方法，如基于激光跟踪仪的测量方法[12-14]、球杆仪测量方法[10,15-17]、R-test装置测量方法[18-19]、触控触发探头[20-21]和试件加工测试[22-23]等测量方法。与其他方法相比，球杆仪的测量方法更加容易建立机床PIGEs与球杆仪杆长变化量之间的数学模型，因此球杆仪被广泛应用于五轴机床旋转轴几何误差的测量。
目前，针对旋转轴PIGEs的定义有“绝对表示”和“相对表示”这两种定义方法。IBARAKI等[24]总结了“相对表示”的定义方法在简化运动学的描述方面具有优势，但不利于找到最小的定位误差集。用“相对表示”误差定义法建立每个旋转轴的数学模型时，两旋转轴之间的误差耦合影响测量的难度，且误差测量上的精确度取决于刀具球与工件球能否放入理想位置。项四通[25]通过控制平动轴静止和两旋转轴同时运动的方法建立了球杆仪与机床旋转轴PIGEs的数学模型，最终得到8项PIGEs的解析解，但其精确度取决于刀具球球心能否精确放置在两旋转轴轴心位置。TSUTSUMI等[16]利用球杆仪测量时由多轴同时控制，使它在每个旋转轴的轴向、切向和径向敏感方向上保持不变，最后成功识别出旋转轴8项PIGEs，但多轴控制会增加平动轴与旋转轴几何误差之间的耦合影响，对准确辨识旋转轴PIGEs会产生一定影响。
相对于“相对表示”的定义方法，“绝对表示”法可以轻松找到最小的定位误差集，通过控制单一旋转轴运动可以在模型的建立上分离出机床拓扑结构中末端运动体对前一级运动体的影响。LEE等[10]通过控制单旋转轴运动的方法建立了PIGEs数学模型，测量结果的标准不确定度分析验证了模型的有效性，该方法在提高辨识精度的同时增加了测量的难度，须使球杆仪安装误差最小化且该方法在C轴辨识模型上没有考虑A轴几何误差的影响，认为C轴的几何误差不受A轴影响。OSEI等[26]通过控制单一旋转轴和采取三种不同的测量方法建立了球杆仪杆长变化量的数学模型，通过实验验证了该模型对PIGEs识别精度高，能有效地提高五轴机床的精度。在模型的建立上该方法在A轴辨识模型中考虑了C轴几何误差的影响，以及球杆仪杆长计算没有保证工件球与刀具球的坐标一致性(一个是真实坐标，另一个是理想坐标)，这会对PIGEs的识别精度有影响。LEE等[15]考虑了A轴与C轴建模时两者间的耦合关系，在A轴模型的建立上分离了C轴误差的影响，同时在杆长变化量上保证了工件球与刀具球的坐标计算一致性，提高了辨识精度，但忽略了工件球端在所测轴的实际位置上的测量难度。综上所述，无论是“相对表示”还是“绝对表示”误差定义的方法建模，在计算球杆仪杆长变化量时会存在刀具球与工件球的坐标参考系不一致问题，以及引入平动轴几何误差的问题，这些问题会影响机床旋转轴PIGEs的辨识精度。除此之外，在实验测量上还要考虑操作的简便性，以便显著缩减测量的时间。
针对以上问题，本文提出了一种双转台五轴机床在机床坐标系下的旋转轴位置无关误差(PIGEs)辨识模型。在机床坐标系下建立刀具球与工件球的实际初始坐标，通过逆矩阵反求工件球在所测轴的实际初始位置，保证刀具球与工件球的坐标计算在同一个坐标系中。然后控制单一旋转轴运动，采用齐次坐标变换理论基于4种测量模式建立了包含球杆仪安装误差与旋转轴PIGEs的杆长变化量数学模型。在A轴测量模式的数学模型中分离出C轴几何误差的影响，在C轴测量模式的数学模型中考虑A轴几何误差的影响，并对4种数学模型进行仿真，分析了PIGEs对各个测量模式的影响。最后通过实验数据对比文献[14]中建模方法的拟合效果，得出旋转轴PIGEs的大小后对PIGEs的4个位置偏差进行补偿实验验证，并与不同文献的误差补偿效果进行了对比。
1 机床结构及误差定义
本文实验采用球杆仪测量五轴机床旋转轴PIGEs。在建立旋转轴PIGEs辨识模型之前需要分析机床的结构，设定坐标系及定义误差类型。
1.1 机床结构和旋转轴PIGEs定义
本文研究的双转台五轴数控机床结构如图1a所示，由三个平动轴X、Y、Z和两个旋转轴A、C组成，运动链结构为X-Y-Z-A-C。在理想情况下，将机床坐标系(MCS)S_M(O_M X_M Y_M Z_M)定义为原点位于A轴和C轴交点的参考系统，A轴坐标系(ACS)S_A(O_A X_A Y_A Z_A)原点对应于S_M的坐标原点，C轴坐标系(CCS)S_C(O_C X_C Y_C Z_C)在转盘中心上，令工件坐标系(WCS)S_W(O_W X_W Y_W Z_W)与C轴坐标系S_C重合，其中Z_CM为机床到C轴之间的距离(即A轴到工作台之间的距离)。当机床处于机床零点时，所有线性轴坐标系定义为与机床坐标系S_M共轴且原点相同。
图1 坐标系的建立 Fig. 1 Establishment of coordinate systems
(a)五轴机床模型;(b)坐标系示意图
旋转轴PIGEs是由机床在组装时旋转轴发生微小偏差产生的。每个旋转轴有4个PIGEs，包括2个位置误差和2个平行度误差，如图2所示。其中，δ_YA、δ_ZA分别为A轴轴线在机床坐标系S_M的Y_M和Z_M方向上的位置误差;S_YA、S_ZA分别为A轴轴线依次绕着Y_M轴和Z_M轴的平行度误差;δ_XC、δ_YC分别为C轴轴线在A轴坐标系S_A的X_A和Y_A方向上的位置误差;S_XC、S_YC分别为C轴轴线依次绕着X_A轴和Y_A轴的平行度误差。
图2 旋转轴误差 Fig. 2 Rotary axis error
(a)A轴位置无关误差;(b)C轴位置无关误差
1.2 球杆仪安装误差定义
球杆仪的安装误差主要集中在主轴端刀具球[27]。主轴一端的刀具杯经过刀柄夹紧时偏离主轴轴线，使得球杆仪的刀具球中心分别在绝对坐标系S_M的X和Y方向上产生安装误差e_X和e_Y，工件球端的中心座经过校准球校准锁紧后，所产生的误差可以忽略不计[28]。如图3所示，安装球杆仪时使主轴侧刀具球球心在Z方向与坐标系S_M原点在同一高度。
图3 球杆仪安装误差 Fig. 3 Ballbar installation error
2 旋转轴PIGEs辨识及仿真分析
2.1 测量模式
本文在旋转轴PIGEs测量的过程中采用的是控制1个旋转轴旋转且其他4个轴保持不动的方法，减少了机床其他运动轴几何误差对此旋转轴测量结果的影响，提高了测量结果的精确度。
测量模式如图4所示，4种测量模式与文献[15]中方法相同，每个旋转轴有2种测量模式，分别测量了1个旋转轴的2个位置误差和2个平行度误差。测量每个旋转轴的2个位置误差时，安装在主轴侧的刀具球球心与2个旋转轴轴心理想相交点重合且固定不动，此外安装在工作台侧的工件球随着旋转轴转动，如图4a、图4c所示。测量每个旋转轴的2个平行度误差时，安装在主轴侧的刀具球球心与所测旋转轴理想轴心重合且固定不动，刀具球球心B₁到C轴轴线和A轴轴线的距离分别为Z_B₁C=100mm和Z_B₁A=100mm，此外安装在工作台侧的工件球随着旋转轴转动，如图4b、图4d所示。4种路径可以测量出2个旋转轴的8项PIGEs。
图4 测量模式 Fig. 4 Measurement mode
(a)模式一测量A轴位置误差;(b)模式二测量A轴平行度偏差;(c)模式三测量C轴位置误差;(d)模式四测量C轴平行度偏差
2.2 误差建模
本文研究的双转台五轴机床的结构如图1a所示，基于齐次变换矩阵(HTM)理论建立了包含球杆仪安装误差与旋转轴PIGEs的球杆仪杆长变化量(ΔL)的数学模型。首先在建立模型时以理想情况下A、C轴轴线相交点O_M为机床坐标系S_M原点，如图1所示，在球杆仪安装误差与旋转轴的影响下建立刀具球球心B₁和工件球球心B₂在S_M坐标系下的实际初始坐标，通过逆矩阵反求工件球在所测轴的实际初始位置。然后控制单一旋转轴运动，通过齐次坐标变换求出工件球旋转后在机床坐标系下的实际坐标。最后根据实验数据拟合可得出PIGEs的大小。用这种方式建模，使刀具球与工件球的坐标计算保证在同一个坐标系中，同时在测量时不用反复寻找C轴的准确位置，减少了测量时间。
在机床坐标系下，由于C旋转轴和工作台合为一体并安装在A旋转轴上，C轴的姿态坐标受A轴误差的影响，A轴姿态不受C轴误差的影响。因此在A轴建模时不受C轴几何误差的影响，相反，在C轴建模时受A轴几何误差的影响。在这种建模方法下，可以明确2个旋转轴PIGEs之间的影响。
球杆仪的测量原理是利用传感器将刀具球与工件球之间的长度变化实时传送出来。在球杆仪测量时，工具杯固定在主轴一端，磁力中心座固定在工作台上。刀具球安装在工具杯上不动，且安装在中心座上随着旋转轴旋转。由1.2节可知，球杆仪在工具杯端存在安装误差，使工具杯分别在S_M坐标系的X、Y方向偏移e_X、e_Y，和工件端的安装误差忽略不计。在S_M坐标系下，B₁的坐标不变，B₂的坐标随着旋转轴的转动而改变。
模式一测量时使B₁位置与S_M坐标系原点位置在Z方向位于同一高度上。理想情况下B₂绕着B₁做圆周运动，实际测量过程中B₂是绕着A轴做圆周运动的，杆长B₁B₂的变化量随着B₂绕A轴旋转而变化，球杆仪的实际测量路径如图5所示。
图5 A轴实际测量路径 Fig. 5 Actual measurement path for rotary axis A
模式一下B₂在S_M坐标系下的初始坐标为：
[^{M} P_{B_{2\_origin}}^{(1)}=left[egin{array}{c}X_{1\_origin} \ Y_{1\_origin} \ Z_{1\_origin} \ 1 end{array} ight]=left[egin{array}{c}e_{X} \ e_{Y}-L \ 0 \ 1 end{array} ight]]
式中，L为球杆仪的公称长度。
S_A与S_M的坐标系变换矩阵为：
[
egin{aligned}
^M_A T & =left[egin{array}{cccc}
1 & -S_{Z A} & S_{Y A} & 0 \
S_{Z A} & 1 & 0 & delta_{Y A} \
-S_{Y A} & 0 & 1 & delta_{Z A} \
0 & 0 & 0 & 1
end{array} ight] cdot \
& left[egin{array}{cccc}
1 & 0 & 0 & 0 \
0 & cos alpha_{A} & -sin alpha_{A} & 0 \
0 & sin alpha_{A} & cos alpha_{A} & 0 \
0 & 0 & 0 & 1
end{array} ight]
end{aligned}
]
式中，α_A为A轴旋转的角度。
模式一下B₂在S_A坐标系下的初始坐标为：
[^{A} P_{B_{2\_origin}}^{(1)}=({^M_A T})^{-1} {^M} P_{B_{2\_origin}}^{(1)} quad (3)]
由于S_YA、S_ZA、δ_ZA、δ_YA、e_X、e_Y等的值都是微小量，所以去掉二阶微小量可得B₂在S_M参考坐标系下的坐标约为：
[
egin{gathered}
^M P_{B_{2}}^{(1)}=left[egin{array}{c}
X_{1} \
Y_{1} \
Z_{1} \
1
end{array} ight]={^M_A T} {^A} P_{B_{2\_origin}}^{(1)} approx \
left[egin{array}{c}
e_{X}-L S_{Z A}+Lleft(S_{Y A} sin alpha_{A}-S_{Z A} cos alpha_{A} ight) \
left(e_{Y}-L-delta_{Y A} ight) cos alpha_{A}+delta_{Z A} sin alpha_{A}+delta_{Y A} \
left(e_{Y}-L-delta_{Y A} ight) sin alpha_{A}-delta_{Z A} cos alpha_{A}+delta_{Z A} \
1
end{array} ight] quad (4)
end{gathered}
]
刀具球是固定不变的，因此模式一下B₁在S_M参考坐标系下的坐标为：
[^M P_{B_{1}}^{(1)}=left[egin{array}{c}e_{X} \ e_{Y} \ 0 \ 1end{array} ight] quad(5)]
球杆仪工具球与工件球的坐标都在S_M坐标系下，这样可以提高模型的辨识精度。
球杆仪的读数对应于B₁与B₂之间的实际距离，则球杆仪的长度可以通过下式确定：
[
egin{gathered}
L_{1}=L+Delta L_{1}=left|P_{B_{1}} P_{B_{2}} ight|=left|^M P_{B_{2}}^{(1)}-^M P_{B_{1}}^{(1)} ight|= \
left{L^{2}left(S_{Z A}+S_{Y A} sin alpha_{A}-S_{Z A} cos alpha_{A} ight)^{2}+ ight. \
{left[left(e_{Y}-L-delta_{Y A} ight) cos alpha_{A}+delta_{Z A} sin alpha_{A}+delta_{Y A}-e_{Y} ight]^{2}+} \
left.left[left(e_{Y}-L-delta_{Y A} ight) sin alpha_{A}-delta_{Z A} cos alpha_{A}+delta_{Z A} ight]^{2} ight}^{frac{1}{2}} quad(6)
end{gathered}
]
式中，L₁为模式一下的球杆仪杆长，ΔL₁为模式一下的球杆仪杆长变化量。
同理，在模式二下只改变B₁在S_A坐标系下X方向上的距离(即B₁到C轴的距离)Z_B₁C=100mm，其他条件皆同模式一，因此B₂在S_M坐标系下的初始位置坐标为：
[^M P_{B_{2\_origin}}^{(2)}=left[egin{array}{c}X_{2\_origin} \ Y_{2\_origin} \ Z_{2\_origin} \ 1end{array} ight]=left[egin{array}{c}Z_{B_{1} C}+e_{X} \ e_{Y}-L \ 0 \ 1end{array} ight] quad(7)]
B₁在S_M坐标系下的坐标为：
[^M P_{B_{1}}^{(2)}=left[egin{array}{c}Z_{B_{1} C}+e_{X} \ e_{Y} \ 0 \ 1end{array} ight] quad(8)]
模式二下球杆仪的长度可以由下式确定：
[
egin{gathered}
L_{2}=L+Delta L_{2}=left|P_{B_{1}} P_{B_{2}} ight|=left|^M P_{B_{2}}^{(2)}-^M P_{B_{1}}^{(2)} ight| approx \
left{L^{2}left(S_{Z A}+S_{Y A} sin alpha_{A}-S_{Z A} cos alpha_{A} ight)^{2}+ ight. \
{left[left(e_{Y}-L-delta_{Y A}-Z_{B_{1} C} S_{Z A} ight) cos alpha_{A}+Z_{B_{1} C} S_{Z A}- ight.} \
left.left(Z_{B_{1} C} S_{Y A}-delta_{Z A} ight) sin alpha_{A}+delta_{Y A}-e_{Y} ight]^{2}+ \
{left[left(e_{Y}-L-delta_{Y A}-Z_{B_{1} C} S_{Z A} ight) sin alpha_{A}+delta_{Y A}-e_{Y} ight]^{2}+} \
left.left(Z_{B_{1} C} S_{Y A}-delta_{Z A} ight) cos alpha_{A}-Z_{B_{1} C} S_{Y A}+delta_{Z A} ight]^{2} ight}^{frac{1}{2}} quad(9)
end{gathered}
]
式中，ΔL₂为模式二下的球杆仪杆长变化量。
在模式三下测量C轴时，刀具球球心位置与模式一中刀具球球心的位置相同。由于存在机床旋转轴装配误差及球杆仪安装误差，因此实际测量过程中B₂绕着C轴中心做圆周运动，杆长B₁B₂的变化量随着球B₂绕C轴旋转而变化。C轴实际测量路径如图6所示。
图6 C轴实际测量路径 Fig. 6 Actual measurement path for rotary axis C
由于C轴受到A轴的影响，因此在模型的建立上要考虑A轴几何误差，模式三下B₂在S_M坐标系下的初始坐标为：
[^M P_{B_{2\_origin}}^{(3)}=left[egin{array}{c}X_{3\_origin} \ Y_{3\_origin} \ Z_{3\_origin} \ 1end{array} ight]=left[egin{array}{c}L+e_{X} \ e_{Y} \ 0 \ 1end{array} ight] quad(10)]
S_C与S_A的坐标系变换矩阵为：
[
egin{aligned}
^A_C T & =left[egin{array}{cccc}
1 & 0 & S_{Y C} & delta_{X C} \
0 & 1 & -S_{X C} & delta_{Y C} \
-S_{Y C} & S_{X C} & 1 & 0 \
0 & 0 & 0 & 1
end{array} ight] cdot \
& left[egin{array}{cccc}
cos alpha_{C} & -sin alpha_{C} & 0 & 0 \
sin alpha_{C} & cos alpha_{C} & 0 & 0 \
0 & 0 & 1 & -Z_{C M} \
0 & 0 & 0 & 1
end{array} ight]
end{aligned}
]
式中，α_C为C轴旋转的角度。
模式三下B₂在S_C坐标系下的初始坐标为：
[^C P_{B_{2\_origin}}^{(3)}=({^M_A T} {^A_C T})^{-1} {^M} P_{B_{2\_origin}}^{(3)} quad (12)]
模式三下B₁在S_M坐标系下的坐标为：
[^M P_{B_{1}}^{(3)}=left[egin{array}{c}e_{X} \ e_{Y} \ 0 \ 1end{array} ight] quad(13)]
C轴旋转后，模式三下B₂在S_M坐标系下的坐标为：
[^M P_{B_{2}}^{(3)}=left[egin{array}{c}X_{3} \ Y_{3} \ Z_{3} \ 1end{array} ight]={^M_A T} {^A_C T} {^C} P_{B_{2\_origin}}^{(3)} quad (14)]
所以在模式三中，球杆仪的长度可以由下式求出：
[
egin{gathered}
L_{3}=L+Delta L_{3}=left|P_{B_{1}} P_{B_{2}} ight|=left|^M P_{B_{2}}^{(3)}-^M P_{B_{1}}^{(3)} ight| approx \
left{left[left(L+e_{X}-delta_{X C} ight) cos alpha_{C}+left(delta_{Y A}+delta_{Y C}+e_{Y} ight) sin alpha_{C}+ ight. ight. \
left.left.delta_{X C}-e_{X} ight]^{2}+left[left(L+e_{X}-delta_{X C} ight) sin alpha_{C}- ight. ight. \
left.left(delta_{Y A}+delta_{Y C}-e_{Y} ight) cos alpha_{C}+delta_{Y A}+delta_{Y C}-e_{Y} ight]^{2}+ \
L^{2}left[left(S_{Y A}+S_{Y C} ight)left(1-cos alpha_{C} ight)+S_{X C} sin alpha_{C} ight]^{2} ight}^{frac{1}{2}} quad(15)
end{gathered}
]
式中，ΔL₃为模式三下的球杆仪杆长变化量。
同理，在模式四中只改变B₁在S_C坐标系下Z轴方向上的距离(即B₁到A轴的距离)Z_B₁A=100mm，其他条件皆同模式三，则B₂在S_M坐标系下的初始坐标为：
[^M P_{B_{2\_origin}}^{(4)}=left[egin{array}{c}X_{4\_origin} \ Y_{4\_origin} \ Z_{4\_origin} \ 1end{array} ight]=left[egin{array}{c}L+e_{X} \ e_{Y} \ Z_{B_{1} A} \ 1end{array} ight] quad(16)]
模式四下B₁在S_M坐标系下的坐标为：
[^M P_{B_{1}}^{(4)}=left[egin{array}{c}e_{X} \ e_{Y} \ Z_{B_{1} A} \ 1end{array} ight] quad(17)]
因此，C轴旋转后，模式四下球杆仪的长度可由下式求出：
[
egin{gathered}
L_{4}=L+Delta L_{4}=left|P_{B_{1}} P_{B_{2}} ight|=left|^M P_{B_{2}}^{(4)}-^M P_{B_{1}}^{(4)} ight| approx \
left(left{left[L+e_{X}-delta_{X C}-Z_{B_{1} A}left(S_{Y A}+S_{Y C} ight) ight] cos alpha_{C}+ ight. ight. \
left(delta_{Y A}+delta_{Y C}-e_{Y}-Z_{B_{1} A} S_{X C} ight) sin alpha_{C}+ \
Z_{B_{1} A}left(S_{Y A}+S_{Y C} ight)+delta_{X C}-e_{X} ight}^{2}+left(left[L+e_{X}- ight. ight. \
left.delta_{X C}-Z_{B_{1} A}left(S_{Y A}+S_{Y C} ight) ight] sin alpha_{C}- \
left(delta_{Y A}+delta_{Y C}-e_{Y}-Z_{B_{1} A} S_{X C} ight) cos alpha_{C}- \
Z_{B_{1} A} S_{X C}+delta_{Y A}+delta_{Y C}-e_{Y} ight}^{2}+L^{2}left[S_{Y A}+S_{Y C}+ ight. \
left.left.S_{X C} sin alpha_{C}-left(S_{Y A}+S_{Y C} ight) cos alpha_{C} ight]^{2} ight)^{frac{1}{2}} quad(18)
end{gathered}
]
式中，ΔL₄为模式四下的球杆仪杆长变化量。
综上，式(6)、式(9)、式(15)、式(18)给出了球杆仪的杆长计算数学模型，在此基础上辨识旋转轴PIGEs。
2.3 各误差对不同测量模式球杆仪变化量的影响
为了研究8项PIGEs误差对4种测量模式中ΔL的影响，本文进行了仿真分析。表1列出了用于仿真的机床基本结构参数和误差数值，所有的PIGEs位置误差设为±10μm，而所有的平行度误差设为±0.005°。
表1 仿真参数 Tab. 1 Parameters for simulation
| 参数 | 数值 |
| 球杆仪的标称长度L (mm) | 100 |
| 模式2中刀具球球心到C轴的距离Z_B₁C (mm) | 100 |
| A轴的旋转范围 | [-70°, 15°] |
| C轴的旋转范围 | [0°, 360°] |
| A轴到工作台的距离Z_CM (mm) | 50 |
| 模式4中刀具球球心到A轴的距离Z_B₁A (mm) | 100 |
| 位置误差δ_YA, δ_ZA, δ_XC, δ_YC (μm) | ±10 |
| 平行度误差S_YA, S_ZA, S_XC, S_YC (°) | ±0.005 |
由于C旋转轴和工作台合为一体并安装在A旋转轴上，C轴的姿态坐标受A轴几何误差的影响，A轴姿态不受C轴几何误差的影响，因此表2中模式一和模式二的仿真没有C轴的几何误差。
表2 8项PIGEs对四种测量模式的影响 Tab. 2 Influences of the 8 PIGEs on four measurement models
| 测量模式 | S_YA=±0.005° | S_ZA=±0.005° | δ_YA=±10μm | δ_ZA=±10μm | S_XC=±0.005° | S_YC=±0.005° | δ_XC=±10μm | δ_YC=±10μm |
| 模式一 | - | - | ● | ● | - | - | - | - |
| 模式二 | ● | ● | ● | ● | - | - | - | - |
| 模式三 | ● | - | ● | - | - | - | ● | ● |
| 模式四 | ● | - | ● | - | ● | ● | ● | ● |
注：-无误差;●正误差;▲负误差
模式一中，C轴及S_YA、S_ZA对ΔL无影响，δ_YA、δ_ZA对ΔL有影响，因此可以用模式一辨识出A轴的两个位置误差δ_YA和δ_ZA。模式二中，C轴对ΔL无影响，S_YA、S_ZA和δ_YA、δ_ZA对ΔL都有影响，结合模式一、模式二可以辨识出A轴的两个平行度误差S_YA和S_ZA。模式三中，空白处表示该误差对ΔL无影响，S_YA、S_XC、S_YC对ΔL无影响，δ_YA、δ_XC、δ_YC对ΔL有影响，且δ_YA、δ_YC两者的影响是一致的，因此可以利用模式一辨识出δ_YA以及利用模式三可以辨识出C轴的两个位置误差δ_XC和δ_YC。模式四中，空白处表示该误差对ΔL无影响，S_YA、δ_YA、S_XC、S_YC和δ_XC、δ_YC对ΔL都有影响，结合模式一、模式二、模式三、模式四可以辨识出C轴PIGEs的两个平行度误差S_XC和S_YC。
仿真结果表明，刀具球在两旋转轴轴线交点位置时，平行度误差不会影响球杆仪的杆长变化量，以此可辨析出旋转轴PIGEs中的位置误差。
2.4 几何误差的识别
分别对4个测量路径进行实验，每次实验时旋转轴转动5°记录一次球杆仪的杆长L_{A/C}。由于不可能完全补偿平动轴误差，而在主轴运动100mm时难免会出现一点小误差，因此球杆仪的公称长度L为0°初始值示数。
对式(6)等号两边平方，去掉二阶以上的误差值可得：
[
egin{gathered}
Delta L_{1}=left(e_{Y}-delta_{Y A} ight) cos alpha_{A}-delta_{Z A} sin alpha_{A}+delta_{Y A}-e_{Y}= \
left[egin{array}{lll}
cos alpha_{A} & sin alpha_{A} & 1
end{array} ight]left[egin{array}{c}
e_{Y}-delta_{Y A} \
-delta_{Z A} \
delta_{Y A}-e_{Y}
end{array} ight] quad(19)
end{gathered}
]
由表2可知，模式一中只有δ_YA、δ_ZA对ΔL有影响，根据路径一实验数据结合式(6)、式(19)，运用最小二乘算法通过MATLAB自定义拟合函数Fittype求出δ_YA、δ_ZA。
由式(9)可得：
[
egin{gathered}
Delta L_{2}=left(delta_{Y A}-e_{Y}+Z_{B_{1} C} S_{Z A} ight)+left(e_{Y}-delta_{Y A}- ight. \
left.Z_{B_{1} C} S_{Z A} ight) cos alpha_{A}+left(Z_{B_{1} C} S_{Y A}-delta_{Z A} ight) sin alpha_{A}= \
{left[egin{array}{lll}
cos alpha_{A} & sin alpha_{A} & 1
end{array} ight]left[egin{array}{c}
e_{Y}-delta_{Y A}-Z_{B_{1} C} S_{Z A} \
Z_{B_{1} C} S_{Y A}-delta_{Z A} \
delta_{Y A}-e_{Y}+Z_{B_{1} C} S_{Z A}
end{array} ight] quad(20)}
end{gathered}
]
模式二中S_YA、S_ZA、δ_YA、δ_ZA对ΔL都有影响，根据路径二实验数据结合式(9)、式(20)，将模式一求出的δ_YA、δ_ZA代入式，运用最小二乘算法通过MATLAB自定义拟合函数Fittype求出S_YA、S_ZA。
由式(15)可得：
[
egin{gathered}
Delta L_{3}=left(e_{X}-delta_{X C} ight)+left(delta_{X C}-e_{X} ight) cos alpha_{C}+ \
left(delta_{Y A}+delta_{Y C}-e_{Y} ight) sin alpha_{C}= \
{left[egin{array}{lll}
cos alpha_{C} & sin alpha_{C} & 1
end{array} ight]left[egin{array}{c}
delta_{X C}-e_{X} \
delta_{Y A}+delta_{Y C}-e_{Y} \
e_{X}-delta_{X C}
end{array} ight] quad(21)}
end{gathered}
]
模式三中只有δ_XC、δ_YA、δ_YC对ΔL都有影响，根据路径三实验数据结合式(15)、式(21)，将模式一求出的δ_YA代入式，运用最小二乘算法通过MATLAB自定义拟合函数Fittype求出δ_XC、δ_YC。
由式(18)可得：
[
egin{gathered}
Delta L_{4}=left[e_{X}-delta_{X C}-Z_{B_{1} A}left(S_{Y A}+S_{Y C} ight) ight]+ \
{left[delta_{X C}-e_{X}+Z_{B_{1} A}left(S_{Y A}+S_{Y C} ight) ight] cos alpha_{C}+} \
left(delta_{Y A}+delta_{Y C}-e_{Y}-Z_{B_{1} A} S_{X C} ight) sin alpha_{C}= \
{left[egin{array}{lll}
cos alpha_{C} & sin alpha_{C} & 1
end{array} ight]left[egin{array}{c}
delta_{X C}-e_{X}+Z_{B_{1} A}left(S_{Y A}+S_{Y C} ight) \
delta_{Y A}+delta_{Y C}-e_{Y}-Z_{B_{1} A} S_{X C} \
e_{X}-delta_{X C}-Z_{B_{1} A}left(S_{Y A}+S_{Y C} ight)
end{array} ight]} \
(22)
end{gathered}
]
模式四中S_XC、S_YA、S_YC、δ_XC、δ_YA、δ_YC对ΔL都有影响，根据路径四实验数据结合式(18)、式(22)，将模式一求出的δ_YA、S_YA及模式三求出的δ_XC、δ_YC代入式，运用最小二乘算法通过MATLAB自定义拟合函数Fittype求出S_XC、S_YC。
使用Fittype函数自定义拟合函数，可以满足线性拟合和非线性拟合。Fittype函数具有很灵活的配置，基本满足各种复杂场景，在复杂函数中可以简单快速求得问题答案。
3 实验验证与补偿
3.1 球杆仪的安装误差测量
本实验在捷甬达V545Ⅲ五轴数控机床上进行了旋转轴PIGEs测量，机床结构如图1a所示，实验验证模式如图4中4种测量模式所示。在测量旋转轴之前，通过丝杆螺距补偿将五轴机床3个平动轴方向的定位误差补偿进去。表3给出了测试所用的球杆仪规格，实验参数与表1中仿真参数一致。
表3 球杆仪规格 Tab. 3 Specifications of double ballbar
| 球杆仪型号 | Renishaw QC20GW |
| 标称长度 (mm) | 100 |
| 校准长度 (mm) | 99.9957 |
| 传感器分辨力 (μm) | 0.1 |
| 测量精度 (μm) | ±0.5 (20℃);±(0.8 + 0.4%读数)(0~40℃) |
| 测量范围 | ±5mm |
本次实验测试前，机床要按照标准预热程序达到热稳定性，实验室温度应调至20℃。安装完球杆仪后，利用千分表和球杆仪来测量球杆仪安装误差，如图7所示。
图7 安装误差测量 Fig. 7 Measurement for installation errors
(a)Z_t测量;(b)e_X和e_Y测量
刀具球B₁置于主轴端，工件球B₂距刀具球B₁在X方向上的距离为100mm，通过千分表测量工具杯的高度Z_t，使B₁球心置于绝对坐标系Z轴原点上，转动主轴使B₁绕主轴轴线旋转，记录球杆仪B₁B₂在主轴旋转角度为0°时的长度示数L和主轴转动过程中最大值L_max、最小值L_min及其对应主轴旋转的角度θ₁、θ₂(θ₁和θ₂的大小相差180°)，重复测量两次，数值取平均值。安装误差的计算原理如图8所示，根据下式可计算得出球杆仪的安装误差e_X、e_Y：
[R=frac{L_{max}-L_{min}}{2} quad(23)]
[
cos left(angle O B_{1} B_{2} ight)=frac{R^{2}+L^{2}-left(L_{min}+R ight)^{2}}{2 R L} quad(24)
]
[e_{X}=-R cos left(angle O B_{1} B_{2} ight) quad(25)]
[
e_{Y}= egin{cases}R sin left(angle O B_{1} B_{2} ight) & 0^{circ}< heta_{2} leq 180^{circ} \ -R sin left(angle O B_{1} B_{2} ight) & 180^{circ}< heta_{2}<360^{circ}end{cases}
quad(26)
]
图8 安装误差测量原理 Fig. 8 Principle of installation error measurement
通过球杆仪和千分表测量球杆仪安装误差，操作简单方便，缩短了实验测量时间。根据测量所得，主轴侧工具杯的高度为Z_t=205.257mm，工具杯的安装误差e_X、e_Y分别为-0.0403mm、-0.0765mm。
3.2 PIGEs实验测量与补偿验证
以C轴为例，不补偿主轴侧工具杯安装误差而直接测量旋转轴，控制A轴静止，C轴从0°旋转到360°，重复两次，记录C轴每旋转5°时球杆仪长度的读数L_C，最后数据处理时取平均值，测量过程如图9所示。
图9 旋转轴测量实验 Fig. 9 Rotary axis measurement experiment
(a)测量C轴位置偏差;(b)测量C轴平行度误差
由2.4节可知，将所测得的C轴旋转轴转动角度和L_C数据，根据式(15)、式(21)编程出来，通过使用MATLAB中Fittype函数拟合可得未知量。
对比文献[15]，当工件球端安装误差忽略不计时，此时实验数据拟合如图10所示，可以看出，当工件球端安装误差忽略不计时，文献[15]的数据拟合效果不佳，说明文献[15]在建模的过程中不能忽略工件球端的安装误差，会影响机床旋转轴PIGEs的辨识精度。而本文在建模时，忽略工件球端的安装误差，结果表明模型函数拟合效果与实验数据相吻合，模型函数拟合效果优于文献[15]。
图10 模型函数拟合效果对比 Fig. 10 Comparison of model function fitting effects
通过对机床旋转轴进行实验，补偿前和补偿后测得的数据如表4所示。由于机床PIGEs平行度误差比较难补偿，因此本文实验只补偿4个PIGEs位置偏差量。
表4 位置误差在补偿前和补偿后的数值大小 Tab. 4 Values of the position error before and after compensation
| 误差类型 | 补偿前 | 补偿后 |
| δ_XC (μm) | 203.5 | -5.16 |
| δ_YC (μm) | 2.61 | 0.905 |
| S_XC (°) | -0.0017 | - |
| S_YC (°) | -0.0547 | - |
| δ_YA (μm) | -24.89 | 3.94 |
| δ_ZA (μm) | -93.82 | -2.12 |
| S_YA (°) | -0.0025 | - |
| S_ZA (°) | -0.0032 | - |
对比补偿前和补偿后的球杆仪杆长变化量ΔL变化结果如图11所示。
图11 球杆仪伸缩量变化 Fig. 11 Double ballbar expansion changes
(a)模式一ΔL变化;(b)模式二ΔL变化;(c)模式三ΔL变化;(d)模式四ΔL变化
→补偿前;■补偿后;▲理想状态
对比表4与图11，平行度误差S_YA、S_ZA对模式一无影响，平行度误差S_XC、S_YC对模式三无影响，补偿后的ΔL变化和理想状态下的ΔL变化一致，验证了误差模型构建的正确性及误差补偿的可行性。由图11d和表4可知，由于补偿前S_YC较大，本实验只补偿了旋转轴的位置误差，因此补偿后的曲线与理想状态的曲线不一致。由图11a、b、c可知，由于补偿前平行度误差S_YA、S_ZA较小，当补偿完A轴位置误差后，补偿后的曲线与理想状态的曲线一致。实验数据表明本文误差建模方法可行，有利于提高机床旋转轴几何误差的辨识精度。
与其他建模方法对比误差补偿前后的改善情况(只对比“绝对表示”类建模方法中PIGEs的位置误差)，如表5所示。
表5 不同建模方法误差补偿前后的对比 Tab. 5 Comparison of different modeling methods before and after error compensation
| 误差 | 文献[15] | | | 文献[26] | | | 本文方法 | | |
| | 补偿前(μm) | 补偿后(μm) | 提升比(%) | 补偿前(μm) | 补偿后(μm) | 提升比(%) | 补偿前(μm) | 补偿后(μm) | 提升比(%) |
| δ_YA | 25.4 | 1.0 | 96.1 | 27.13 | 8.01 | 70.5 | 24.9 | 3.9 | 84.3 |
| δ_ZA | 9.1 | 0.2 | 97.8 | 1.90 | 0.68 | 64.2 | 93.82 | -2.12 | 97.7 |
| δ_XC | 22.3 | 4.0 | 82.1 | 23.36 | 8.09 | 65.4 | 203.5 | 5.1 | 97.5 |
| δ_YC | 3.3 | 2.0 | 39.4 | 214.20 | 68.12 | 68.2 | 2.6 | 0.9 | 65.4 |
| 整体提升比(%) | - | - | 78.85 | - | - | 67.08 | - | - | 86.5 |
通过与现有“绝对表示”方法的辨识模型进行对比，本文模型在补偿后的误差辨识结果上整体提升比明显高于其他两种方法的整体提升比，位置误差补偿后精度得到更好的改善，补偿后的位置误差最大绝对值由203.5μm减小至5.1μm，所提出的误差辨识模型可以有效提高五轴机床几何误差辨识精度。
4 结论
(1) 本文在机床坐标系下建立刀具球与工件球的实际初始坐标，利用逆矩阵反求工件球在所测轴的实际初始坐标，基于齐次坐标变换理论建立了球杆仪杆长变化量的数学模型，解决了坐标计算不一致性问题，同时在实验上避免了反复寻找旋转轴轴线准确位置的问题。
(2) 通过仿真分析位置无关几何误差(PIGEs)对测量模式中球杆仪杆长变化量的影响，研究结果表明刀具球在两旋转轴轴线交点位置时平行度误差不影响杆长变化量，以此辨识出PIGEs的位置误差。
(3) 通过对比实验数据拟合及不同建模方法的误差补偿效果，补偿后的位置误差最大绝对值由203.5μm减小至5.1μm，本文的建模方法可以有效提高五轴机床旋转轴几何误差辨识精度。
参考文献
[1] BI Qingzhen, HUANG Nuodi, SUN Chao, et al. Identification and Compensation of Geometric Errors of Rotary Axes on Five-axis Machine by On-machine Measurement[J]. International Journal of Machine Tools and Manufacture, 2015, 89: 182-191.
[2] SCHWENKE H, KNAPP W, HAITJEMA H, et al. Geometric Error Measurement and Compensation of Machines: an update[J]. CIRP Annals, 2008, 57(2): 660-675.
[3] LIU Kuo, SUN Mingjia, ZHU Tiejun, et al. Modeling and Compensation for Spindle’s Radial Thermal Drift Error on a Vertical Machining Center[J]. International Journal of Machine Tools and Manufacture, 2016, 105: 58-67.
[4] YANG Jixiang, ALTINTAS Y. A Generalized On-line Estimation and Control of Five-axis Contouring Errors of CNC Machine Tools[J]. International Journal of Machine Tools and Manufacture, 2015, 88: 9-23.
[5] WAN Min, DANG Xuebin, ZHANG Weihong, et al. Optimization and Improvement of Stable Processing Condition by Attaching Additional Masses for Milling of Thin-walled Workpiece[J]. Mechanical Systems and Signal Processing, 2018, 103: 196-215.
[6] ISO. ISO230-1 Test Code for Machine Tools - Part1: Geometric Accuracy of Machines Operating Under No-load or Quasi-static Conditions[S]. Geneva: ISO, 2012.
[7] ISO. ISO230-7 Test Code for Machine Tools - Part7: Geometric Accuracy of Axes of Rotation[S]. Geneva: ISO, 2015.
[8] ISO. ISO10791-1 Test Conditions for Machining Centres - Part1: Geometric Tests for Machines with Horizontal Spindle and with Accessory Heads (Horizontal Z-axis)[S]. Geneva: ISO, 1998.
[9] LIU Yang, WAN Min, XING Wanjing, et al. Identification of Position Independent Geometric Errors of Rotary Axes for Five-axis Machine Tools with Structural Restrictions[J]. Robotics and Computer-Integrated Manufacturing, 2018, 53: 45-57.
[10] LEE K I, YANG S H. Robust Measurement Method and Uncertainty Analysis for Position-independent Geometric Errors of a Rotary axis Using a Double Ball-bar[J]. International Journal of Precision Engineering and Manufacturing, 2013, 14(2): 231-239.
[11] 郭世杰, 张东升. 五轴机床旋转轴几何误差分析与补偿[J]. 工程科学与技术, 2020, 52(2): 130-139. GUO Shijie, ZHANG Dongsheng. Geometric Error Analysis and Compensation of Rotary Axes of Five-axis Machine Tools[J]. Advanced Engineering Sciences, 2020, 52(2): 130-139.
[12] SUN Tao, ZHAI Yapu, SONG Yimin, et al. Kinematic Calibration of a 3-DoF Rotational Parallel Manipulator Using Laser Tracker[J]. Robotics and Computer-Integrated Manufacturing, 2016, 41: 78-91.
[13] WANG Jindong, GUO Junjie. The Identification Method of the Relative Position Relationship between the Rotary and Linear Axis of Multi-axis Numerical Control Machine Tool by Laser Tracker[J]. Measurement, 2019, 132: 369-376.
[14] DENG Ming, LI Huimin, XIANG Sitong, et al. Geometric Errors Identification Considering Rigid-body Motion Constraint for Rotary Axis of Multi-axis Machine Tool Using a Tracking Interferometer[J]. International Journal of Machine Tools and Manufacture, 2020, 158: 103625.
[15] LEE K I, YANG S H. Measurement and Verification of Position-independent Geometric Errors of a Five-axis Machine Tool Using a Double Ball-bar[J]. International Journal of Machine Tools and Manufacture, 2013, 70: 45-52.
[16] TSUTSUMI M, SAITO A. Identification and Compensation of Systematic Deviations Particular to 5-Axis Machining Centers[J]. International Journal of Machine Tools and Manufacture, 2003, 43(8): 771-780.
[17] LEE D G M, ZHU Z, LEE K I, et al. Identification and Measurement of Geometric Errors for a Five-axis Machine Tool with a Tilting Head Using a Double Ball-bar[J]. International Journal of Precision Engineering and Manufacturing, 2011, 12(2): 337-343.
[18] IBARAKI S, OYAMA C, OTSUBO H. Construction of an Error Map of Rotary Axes on a Five-axis Machining Center by Static R-test[J]. International Journal of Machine Tools and Manufacture, 2011, 51(3): 190-200.
[19] BRINGMANN B, KNAPP W. Model-based “Chase-the-ball” Calibration of a 5-Axes Machining Center[J]. CIRP Annals, 2006, 55(1): 531-534.
[20] CHEN Y, MORE P, LIU C. Identification and Verification of Location Errors of Rotary Axes on Five-axis Machine Tools by Using a Touch-trigger Probe and a Sphere[J]. The International Journal of Advanced Manufacturing Technology, 2018, 100(9/12): 2653-2667.
[21] HOLUB M, JANKOVYCH R, ANDRS O, et al. Capability Assessment of CNC Machining Centres as Measuring Devices[J]. Measurement, 2018, 118: 52-60.
[22] IBARAKI S, SAWADA M, MATSUBARA A, et al. Machining Tests to Identify Kinematic Errors on Five-axis Machine Tools[J]. Precision Engineering, 2010, 34(3): 387-398.
[23] PEZESHKI M, AREZOO B. Kinematic Errors Identification of Three-axis Machine Tools Based on Machined WorkPieces[J]. Precision Engineering, 2016, 43: 493-504.
[24] IBARAKI S, KNAPP W. Indirect Measurement of Volumetric Accuracy for Three-axis and Five-axis Machine Tools: a Review[J]. International Journal of Automation Technology, 2012, 6: 110-124.
[25] 项四通. 五轴数控机床空间误差测量、建模与补偿技术研究[D]. 上海: 上海交通大学, 2016. XIANG Sitong. Volumetric Error Measuring, Modeling and Compensation Technique for Five-axis Machine Tools[D]. Shanghai: Shanghai Jiao Tong University, 2016.
[26] OSEI S, WANG W R, DING Q C. A New Method to Identify the Position-independent Geometric Errors in the Rotary Axes of Five-axis Machine tools[J]. Journal of Manufacturing Processes, 2023, 87: 46-53.
[27] 何坤, 何晓虎, 徐凯, 等. 基于参数化建模的球杆仪安装误差与直线轴几何误差分离方法[J]. 计算机集成制造系统, 2021, 27(12): 3475-3483. HE Kun, HE Xiaohu, XU Kai, et al. Separation Method for Installation Errors of Double Ball Bar and Geometric Errors of Linear Axes Based on Parametric Modeling[J]. Computer Integrated Manufacturing System, 2021, 27(12): 3475-3483.
[28] 夏长久, 王时龙, 徐凯, 等. 基于球杆仪单轴运动测量的旋转轴几何误差辨识[J]. 仪器仪表学报, 2020, 41(7): 19-28. XIA Changjiu, WANG Shilong, XU Kai, et al. Geometric Error Identification of Rotary Axes Based on Uniaxial Motion Measurement of Double Ball Bar[J]. Chinese Journal of Scientific Instrument, 2020, 41(7): 19-28.
