采用自动钻铆技术实现飞机壁板的制孔与铆接是提高飞机装配质量和效率的有效途径^[1]. 自20世纪50年代起源，经过几十年的发展，自动钻铆设备已能自动完成定位、制孔、锪窝、送钉、插钉、铆接及检测等一系列工序，已是先进航空制造中不可或缺的重要装备^[2].

国外对自动钻铆设备的研究起步较早，西方航空工业发达的国家已经拥有成熟的自动钻铆设备，例如德国BROTJE公司的MPAC自动钻铆机、美国EI公司的E4000系列自动钻铆机等^[3-4]. 中国早在20世纪70年代初就尝试研制自动钻铆设备，但由于当时工业技术落后，自动钻铆设备未能研制成功. 近年来，国内各大航空制造厂在引进国外自动钻铆设备的同时，加强了与国内科研院所的合作，成功研发出一些较为先进的自动钻铆设备，例如西飞研发的基于E7000型自动钻铆设备的柔性装配托架，成飞研发的塔式五轴数控钻铆系统^[5-7].

自动钻铆系统在实际铆接时受到大载荷的循环施加，钻铆系统的结构会出现一定的弹性形变，导致加工精度下降. 自动钻铆系统需要具有较强的刚度以抵抗这种弹性形变，这就需要定量地分析钻铆系统末端的变形，研究其末端刚度，为其后的控制反馈提供一定的数据支持. Kim等^[8]研究了多体系统中工作空间内的刚度矩阵、Salgado等^[9]将机床刀具、夹具、压紧机构等机床部件的刚度矩阵合成为系统整体刚度矩阵加以分析、闫蓉等^[10]提出了一种多轴加工系统刚度场半解析方法，用以分析工作空间中的刚度特性分布规律、Gao等^[11]通过解析法建立了机床工作空间内的刚度模型，并利用该模型对不同位姿下的加工误差进行了分析、吴石等^[12]利用多体小变形理论建立了机床的综合刚度模型，并通过力椭球分析了加工各要素对加工系统综合刚度性能的影响.

以往对机床刚度的研究大多采用多体运动学的方法进行分析，很少全面地考虑机床变形的所有来源、对于本研究卧式双机联合钻铆系统，需要综合考虑机床变形的各个来源建立其末端刚度模型. 针对这一情况，本文运用机器人学相关理论建立其运动学模型，并使用雅克比矩阵法、点传递法、有限元分析法逐一对其柔性运动轴、运动关节、末端执行器等关键部件建立刚度矩阵，在此基础上基于弹性小变形原理建立多轴机床的末端综合刚度模型. 这种末端综合刚度模型可分析不同位姿下系统所产生的变形，以便优化结构设计，必要时通过控制系统给出补偿量，提高铆接精度.

卧式双机联合钻铆系统如图1所示，系统由2台带有外铆接头（用于制孔插钉）和内铆接头（用于压铆）的多轴数控机床组成. 这2台多轴机床分别位于待加工壁板的两侧，之间留有足够的加工空间，待加工的飞机壁板通过AGV车进入工位. 这种布局的自动钻铆系统具有通过性好、加工准备时间短的优点，便于实现流水化生产，极大地提高了自动钻铆的效率^[13-14].

卧式双机联合钻铆系统的刚度主要取决于两侧多轴机床中刚度较差的一个. 由于带内铆接头的多轴机床受载时变形更为明显，本文主要以带内铆接头的多轴机床作为分析对象，并视其末端刚度为卧式双机联合钻铆系统的刚度. 带内铆接头的多轴机床运动轴分布图如图2所示，该多轴机床包含X、Y、Z三个直角坐标移动轴和A、B两个正交转动轴，理论上能够以任意姿态达到加工空间内的任意位置，可满足在飞机壁板复杂曲面上的加工需求.

根据设计要求，多轴机床在工作时需承受高达5 t的压铆力，其运动关节、末端执行器（铆接头）以及柔性运动轴均会产生不同程度的变形，以往对末端坐标系的机器人刚度研究通常仅考虑运动关节变形对末端刚度的影响，运动轴及末端执行器的变形相对于传动机构的变形往往因为数值过小而不加考虑. 对于多轴机床来说，其Z轴和A轴悬伸较长，在加工时会形成悬臂结构导致其变形量较大，因此Z、A两运动轴的变形是末端变形的重要来源之一，不可忽略. 为了提高计算精度，需综合考虑运动轴、传动机构以及末端执行器的变形来建立该多轴机床的末端刚度矩阵.

多轴机床可以看作由3个棱柱联轴节以及2个转动关节所组成的五轴机械臂，因此，可以按照Denavit-Hartenberg方法建立其连杆坐标系^[15-16]. 所建立运动学模型如图3所示，图中的坐标系从参考坐标系到末端坐标系分别为Or、OX、OY、OZ、OA、OB、Ot. 多轴机床的连杆参数如表1所示。表中， ${\theta _i}$ 为两连杆夹角， ${a_{i - 1}}$ 为连杆长度，α_i−1 为连杆扭角， ${d_i}$ 为连杆距离. 在各个连杆坐标系中， ${d_X}$ 为坐标系 ${O_X}$ 的关节变量， ${d_Y}$ 为坐标系 ${O_Y}$ 的关节变量， ${d_Z}$ 为坐标系 ${O_Z}$ 的关节变量， ${\theta _A}$ 为坐标系 ${O_A}$ 的关节变量， ${\theta _B}$ 为坐标系 ${O_B}$ 的关节变量， ${L_{\rm{a}}}$ 为坐标系 ${O_Z}$ 与 ${O_A}$ 原点在 ${Z_Z}$ 轴方向上的距离.

根据机器人学相关理论，相邻坐标系之间的坐标变换可根据下式进行计算：

(1) $\mathit{T}_\mathit{i}^{\mathit{i}-1}{\bf{ = }}\left[ {\begin{array}{*{20}{c}} {{\rm c}{\theta _i}}&{ - {\rm s}{\theta _i}}&0&{{a_{i - 1}}}\\ {{\rm s}{\theta _i}{\rm c}{\alpha _{i - 1}}}&{{\rm c}{\theta _i}{\rm c}{\alpha _{i - 1}}}&{ - {\rm s}{\alpha _{i - 1}}}&{ - {d_i}{\rm s}{\alpha _{i - 1}}}\\ {{\rm s}{\theta _i}{\rm s}{\alpha _{i - 1}}}&{{\rm c}{\theta _i}{\rm s}{\alpha _{i - 1}}}&{{\rm c}{\alpha _{i - 1}}}&{{d_i}{\rm c}{\alpha _{i - 1}}}\\ 0&0&0&1 \end{array}} \right].$

式中： ${\rm{c}}$ 代表三角函数 $\cos $， ${\rm{s}}$ 代表三角函数 $\sin $.

根据如表1所示的连杆参数以及式（1），可得多轴机床从参考坐标系到坐标系的坐标变换为

(2) $\begin{aligned} &{{T}}_{{B}}^{{r}} = {{T}}_{{X}}^{{r}}{{T}}_{{Y}}^{{X}}{{T}}_{{Z}}^{{Y}}{{T}}_{{A}}^{{Z}}{{T}}_{{B}}^{{A}} =\\ &\left[ {\begin{aligned} {{\rm s}{\theta _A}{\rm c}{\theta _B}}&\quad\quad{ - {\rm s}{\theta _A}{\rm s}{\theta _B}}&\quad{{\rm c}{\theta _A}}&\quad\quad{ - {d_Z}}\\ { - {\rm c}{\theta _A}{\rm c}{\theta _B}}\!\!&\quad\quad{{\rm c}{\theta _A}{\rm s}{\theta _B}{\rm s}{\theta _B}}&\quad{{\rm s}{\theta _A}}&\quad\quad\;\;{{d_Y}}\\ { - {\rm s}{\theta _B}}\;\;&\quad\quad\quad{ - {\rm c}{\theta _B}}&\quad 0 \;\;&\quad\quad\;\;{{d_X}}\\ 0 \quad &\quad\quad\quad\quad 0 \;\;&\quad 0\;\; &\quad\quad\;\;\; 1 \end{aligned}} \right]. \end{aligned}$

式（2）描述了从参考坐标系到末端执行器坐标系之间的变换关系. 对于任意两坐标系之间的变换关系，仍可采用类似式（2）的连乘方式得到.

多轴机床在实际工作时需要承受大载荷的循环施加，其运动传递的各个环节都会产生相应的变形. 根据理论推导以及相关的实验数据，多轴机床末端变形主要由结构变形和传动变形所产生，实际表现为加工时末端执行器、柔性连杆和运动关节刚性的不足. 针对这一情况，需要首先确定末端执行器、柔性连杆和运动关节这3个环节的刚度矩阵，为计算整机的末端刚度矩阵提供必要的基础数据支持.

多轴机床中的末端执行器即为布置在该机床末端的铆接头，如图4所示为该铆接头的结构简化图. 为了减少和壁板的干涉，铆接头被设计成细长的杆状结构，在重载的情况下极易产生较大的结构变形. 铆接头的变形主要来源于铆杆、套筒和衬套的组合体，可利用有限元分析的方法对该组合体进行静力分析，得到该铆接头的刚度矩阵. 选择铆接头末端的参考点作为末端坐标系的原点，并以该参考点变形作为铆接头的变形. 运用有限元分析法，在给铆接头各个方向施加30 000 N的载荷及1 000 N·m的转矩后，即可得到铆接头各个方向上受载后的变形.

铆接头的柔度矩阵为其刚度矩阵的逆矩阵，为6×6阶矩阵. 柔度矩阵的每一列分别代表在对应方向上单位力作用下所产生的变形. 基于这个特性，可利用有限元方法对各个方向施加力得到末端参考点各个方向上的变形，依据式（3）可得到柔度矩阵上的对应列：

(3) ${{{S}}_{{i}}}{{{F}}_{{i}}}={{{D}}_{{i}}}.$

式中： ${{{S}}_{{i}}} = {\left[ {\begin{array}{*{20}{c}} {{S_x}},\;\;{{S_y}},\;\;{{S_z}},\;\;{{S_{\theta x}}},\;\;{{S_{\theta y}}},\;\;{{S_{\theta z}}}, \end{array}} \right]^{\rm T}}$ 为柔度矩阵中的列向量； ${{{F}}_{{i}}}$ 为对应方向上的单位力；D_i= ${\left[ {\begin{array}{*{20}{c}} {{d_x}},\;\;{{d_y}},\;\;{{d_z}},\;\;{{\theta _x}},\;\;{{\theta _y}},\;\;{{\theta _z}} \end{array}} \right]^{\rm T}}$ 为各个方向上的位移.

依照式（3）求取柔度矩阵各列后，考虑到铆接头在工作时不会受到转矩，同时忽略高阶微量即可得到铆接头的柔度矩阵，求取其逆后，得到刚度对角矩阵：

(4) ${{K}}{\rm h} \!=\!{\rm diag} \left[ {1.69 \!\times\! {{10}^8},\;4.43 \!\times\! {{10}^5},\;4.41 \!\times\! {{10}^5},\;0,\;0,\;0} \right].$

式中：主对角线上的前3个数值单位为N/m，后3个数值单位为N·m/rad，以下不作特殊说明均采用该单位.

在柔性运动轴中，除Z轴和A轴之外的运动轴变形较小，可以忽略，计算时只需考虑Z轴和A轴的刚度矩阵. 采用类似的有限元法计算柔性轴的柔度矩阵，在30 000 N的载荷下，得到柔性运动轴Z轴及A的应力变形云图如图5、6所示所示. 采用有限元法计算得到的Z轴、A轴刚度矩阵分别如下所示：

(5) ${{{K}}_{{Z}}} = \left[ {\begin{array}{*{20}{c}} {5.48 \times 1{0^7}}&{2.14 \times 1{0^{10}}}&0&0&0&0 \\ {2.14 \times 1{0^{10}}}&{1.48 \times 1{0^{9}}}&0&0&0&{ - 6.18 \times 1{0^{9}}} \\ 0&0&{5.20 \times 1{0^{7}}}&0&{ - 8.96 \times 1{0^{9}}}&0 \\ 0&0&0&{3.86 \times 1{0^{8}}}&0&0 \\ 0&0&{ - 8.96 \times 1{0^{9}}}&0&{1.13 \times 1{0^{8}}}&0 \\ 0&{ - 6.18 \times 1{0^{9}}}&0&0&0&{1.18 \times 1{0^{9}}} \end{array}} \right],$

(6) ${{{K}}_{{A}}} = \left[ {\begin{array}{*{20}{c}} {2.69 \times 1{0^{7}}}&0&{2.22 \times 1{0^{10}}}&0&{3.06 \times 1{0^9}}&0 \\ 0&{ - 7.96 \times 1{0^{11}}}&0&0&0&{ - 1.92 \times 1{0^{10}}} \\ {2.22 \times 1{0^{10}}}&0&{7.21 \times 1{0^8}}&0&{9.93 \times 1{0^{7}}}&0 \\ 0&0&0&{3.06 \times 1{0^9}}&0&0 \\ {3.06 \times 1{0^9}}&0&{9.93 \times 1{0^7}}&0&{ - 1.88 \times 1{0^{10}}}&0 \\ 0&{ - 1.92 \times 1{0^{10}}}&0&0&0&{ - 3.64 \times 1{0^9}} \end{array}} \right].$

多轴机床包含3个移动关节和2个转动关节，串联布置在整个机器上. 在多轴机床受载时，各个关节均会沿其运动方向产生不同程度的偏移. 所谓关节刚度，即为各个关节抵抗这种位置偏移的能力.

一般来说，各个关节的位置偏移是多轴机床末端变形的最主要来源. 为保证各个关节位置偏移求取的准确性，本文采用实验法来确定关节位移，从而确定各个关节刚度，如图7所示即为试验现场. 对于移动关节，千分表的偏移值即为加载时的关节位移. 对于转动关节，测量原理如图8所示，在转动角度较小时，所测量的直线位移可等效于对应的弧长，可由如下公式计算出偏移角度：

(7) $L = R\theta .$

式中：L为转动关节受载后转动的弧长；R为测量点到转动中心的半径； $\theta $ 为转动关节受载后转动的角度.

各个关节的关节广义力可由经典机器人理论得到，即

(8) ${{\tau}}= {{{J}}^{\rm{T}}}{{{F}}}.$

式中： ${{\tau }} = {\left[ {\begin{array}{*{20}{c}} {{f_X}},\;\;{{f_Y}},\;\;{{f_Z}},\;\;{{m_A}},\;\;{{m_B}} \end{array}} \right]^{\rm T}}$ 为关节内力广义力矩阵，对于转动关节，关节内力表现为力矩，移动关节则表现为力； ${{{J}}^{\rm{T}}}$ 为雅克比矩阵的转置； ${{{F}}}$ 为作用在末端坐标系内的广义力.

多轴机床的关节刚度矩阵为一个五阶对角矩阵，如下式所示：

(9) ${{{K}}_{\rm{j}}} = {\rm diag}\left[ {k_{\text{j}}^x,\;\;k_{\text{j}}^y,\;\;k_{\text{j}}^z,\;\;k_{\text{j}}^a,\;\;k_{\text{j}}^b} \right].$

式中： $k_{\rm j}^x$、 $k_{\rm j}^y$、 $k_{\rm j}^z$、 $k_{\rm j}^a$、 $k_{\rm j}^b$ 分别为X轴、Y轴、Z轴、A轴和B轴的关节刚度.

各个关节的关节刚度可由下式求取：

(10) $k_{\text{j}}^i = \tau _{\text{j}}^i/D_{\text{j}}^i.$

式中： $k_{\rm j}^i$ 为关节i的关节刚度； $\tau _{\rm j}^i$ 为作用于关节i上的广义力； $D_{\rm j}^i$ 为关节i受载后的关节位移.

经由试验测量及计算，在多轴机床末端加载10 000 N的情况下，得到各个关节的位移及广义力如表2所示. 根据之前求取的位移与广义力，可得到多轴机床的关节刚度矩阵为

(11) $\begin{split}{{{K}}_{\rm{j}}} =&{\rm diag} \left[ {9.679 \times {{10}^8},\;\;6.467 \times {{10}^8},\;\;5.988 \times {{10}^7},}\right.\\ &\left.{3.818 \times {{10}^4},\;\;4.291 \times {{10}^4}} \right].\end{split}$

在实际工作时，多轴机床运动传递的各个环节均会产生不同程度上的弹性变形，同时这些变形均符合弹性小变形原理，因此多轴机床末端的变形可视为各个环节变形的叠加. 这些产生弹性变形的主要环节包括运动关节的变形、柔性运动轴的变形、末端执行器的变形. 由此，总体的变形在末端坐标 ${O_t}$内可表示为

(12) ${{\delta \upsilon }}_{\rm{1}}^{{t}} = {{\delta \upsilon }}_{\rm{i}}^{{t}} + {{\delta \upsilon }}_{\rm{j}}^{{t}} + {{\delta \upsilon }}_{\rm{h}}^{{t}}.$

式中： ${{\delta \upsilon }}_{\rm{1}}^{{t}}$ 为机床末端执行器坐标系的微小位移； ${{\delta \upsilon }}_{\rm{i}}^{{t}} =$ ${\left[ \!\!\!{\begin{array}{*{20}{c}}{\delta {x_{\rm i}}},{\delta {y_{\rm i}}},{\delta {z_{\rm i}}},{\delta {\theta _{x{\rm i}}}},{\delta {\theta _{y{\rm i}}}},{\delta {\theta _{z{\rm i}}}} \end{array}} \!\!\!\right]^{\rm T}}$表示由于机床柔性运动轴变形所产生的微小位移； ${{\delta \upsilon }}_{\rm{j}}^{{t}} =$ ${\left[\!\!\!\! {\begin{array}{*{20}{c}} {\delta {x_{\rm j}}},{\delta {y_{\rm j}}},{\delta {z_{\rm j}}},{\delta {\theta _{x{\rm j}}}},{\delta {\theta _{y{\rm j}}}},{\delta {\theta _{z{\rm j}}}} \end{array}} \!\!\!\!\right]^{\rm T}}$ 表示由于各个运动关节变形所产生的微小位移； ${{\delta \upsilon }}_{\rm{h}}^{{t}} \!=\! {\left[\! \!\!{\begin{array}{*{20}{c}} {\delta {x_{\rm h}}},{\delta {y_{\rm h}}},{\delta {z_{\rm h}}},{\delta {\theta _{x{\rm h}}}},{\delta {\theta _{y{\rm h}}}},{\delta {\theta _{z{\rm h}}}} \end{array}} \!\!\!\right]^{\rm T}}$表示机床末端执行器变形所产生的微小位移.

对于任意微小位移，存在如下的映射关系：

(13) ${{\delta }}{{{\upsilon }}^{{t}}} = {{{S}}^{{t}}}{{F}}.$

式中： ${{\delta }}{{{\upsilon }}^{{t}}}$为任意微小位移在末端坐标系内的表示； ${{{S}}^{{t}}}$为引发该位移对应环节的末端柔度矩阵.

联立式（12）、（13），可得

(14) $\begin{gathered} {{\delta }}{{{\upsilon }}^{{t}}} = {{\delta \upsilon }}_{\rm{i}}^{{t}} + {{\delta \upsilon }}_{\rm{j}}^{{t}} + {{\delta \upsilon }}_{\rm{h}}^{{t}} \mathop {}\nolimits^{} = {{(S}}_{\rm{i}}^{{t}} + {{S}}_{\rm{j}}^{{t}} + {{S}}_{\rm{h}}^{{t}}{{)F}} = {{{S}}^{{t}}}{{F}} \\ \end{gathered} .$

式中： ${{S}}_{\rm{i}}^{{t}}$ 为柔性运动轴柔度矩阵在末端坐标系内的表示，为柔性运动轴刚度矩阵的逆； ${{S}}_{\rm{j}}^{{t}}$ 为运动关节柔度矩阵在末端坐标系内的表示，为运动关节刚度矩阵的逆； ${{S}}_{\rm{h}}^{{t}}$ 为末端执行器柔度矩阵在末端坐标系内的表示，为末端执行器刚度矩阵的逆.

根据式（14），分别求出柔性运动轴柔度矩阵 ${{S}}_{\rm{i}}^{{t}}$，运动关节柔度矩阵 ${{S}}_{\rm{j}}^{{t}}$，末端执行器柔度矩阵 ${{S}}_{\rm{h}}^{{t}}$ 即可得到多轴机床的综合末端柔度矩阵.

另一侧的多轴机床组成结构与带内铆接头的多轴机床类似，同为五轴结构且变形环节，也包含运动关节、末端执行器及柔性运动轴，也可由该方法得到其柔度矩阵，此处不再赘述.

对于大部分机械臂来说，运动关节在其运动方向上的变形是整个机械臂末端变形的主要来源，其刚度矩阵建立的精度对末端刚度矩阵的影响最大. 经典的关节刚度矩阵忽略了受力所产生的变形对雅克比矩阵的影响，仅适用于小载荷的情况. 本文所述的多轴机床在加工时最大会承受接近50 000 N的大载荷，使用经典的刚度矩阵会产生较大的误差. 针对这种情况，本文采用了种更为精确的算法来建立关节刚度矩阵^[17].

根据虚功原理，关节空间内的广义力做的虚功等于末端坐标系内广义力所做的虚功，即

(15) ${{\delta }}{{{\omega }}_{\rm{j}}} = {{{F}}^{\rm{T}}}{{\delta \upsilon }}_{\rm{j}}^{{t}} = {{{\tau }}^{\rm{T}}}{{Q}}.$

式中： ${{Q}} = {\left[ {\begin{array}{*{20}{c}} {{q_1}},\;{{q_2}},\;{{q_3}},\;{{q_4}},\;{{q_5}} \end{array}} \right]^{\rm T}}$ 为关节的微小位移矩阵.

根据机器人学经典理论，可得

(16) ${{\delta \upsilon }}_{\rm{j}}^{{t}} = {{J}}_{\rm{j}}^{{t}}{{Q}}.$

式中： ${{J}}_{\rm{j}}^{{t}}$ 为雅克比矩阵，表示关节空间与末端坐标系之间的微分位移关系.

联立式（15）、（16），可得广义力的关系为

(17) ${{\tau }} = {\left( {{{J}}_{\rm{j}}^{{t}}} \right)^{\rm T}}{{F}}.$

同时，各个广义力满足以下等式：

(18) ${{\tau }} = {{{K}}_{\rm{j}}}{{Q}},$

(19) ${{F}} = {{K}}_{\rm{j}}^{{t}}{{\delta \upsilon }}_{\rm{j}}^{{t}}.$

式中： ${{K}}_{{j}}^{{t}}$ 为末端坐标系内的关节刚度矩阵；对公式（17）的两侧同时求Q的偏导，可得

(20) $\frac{{\partial {{\tau }}}}{{\partial {{Q}}}} = \left( {\frac{{\partial {{\left( {{{J}}_{\rm{j}}^{{t}}} \right)}^{\rm T}}}}{{\partial {{Q}}}}} \right){{F}} + {\left( {{{J}}_{\rm{j}}^{{t}}} \right)^{\rm T}}\frac{{\partial {{F}}}}{{\partial \left( {{{\delta \upsilon }}_{\rm{j}}^{{t}}} \right)}}\frac{{\partial \left( {{{\delta \upsilon }}_{\rm{j}}^{{t}}} \right)}}{{\partial {{Q}}}}.$

式中： ${{\partial {{\tau }}}}/{{\partial {{Q}}}}$即为 ${{{K}}_{\bf{j}}}$； ${{\partial {{F}}}}/{{\partial \left( {{{\delta \upsilon }}_{\rm{j}}^{{t}}} \right)}}$即为 ${{K}}_{\rm{j}}^{{t}}$； ${{\partial \left( {{{\delta \upsilon }}_{\rm{j}}^{{t}}} \right)}}/{{\partial {{Q}}}}$即为 ${{J}}_{\rm{j}}^{{t}},$ 因此，式（20）可化简为

(21) ${{{K}}_{\bf{j}}} = \left( {\frac{{\partial {{\left( {{{J}}_{\rm{j}}^{{t}}} \right)}^{\rm T}}}}{{\partial {{Q}}}}} \right){{F}} + {\left( {{{J}}_{\rm{j}}^{{t}}} \right)^{\rm T}}{{K}}_{\rm{j}}^{{t}}{{J}}_{\rm{j}}^{{t}}.$

对于静态未加载的多轴机床，式（21）中等式右边的第一项为0，在这种情况下，可得到经典理论中的关节刚度矩阵：

(22) ${{{K}}_{\rm{j}}} = {\left( {{{J}}_{\rm{j}}^{{t}}} \right)^{\rm T}}{{K}}_{\rm{j}}^{{t}}{{J}}_{\rm{j}}^{{t}},$

(23) ${{K}}_{\rm{j}}^{{t}}{{ }} = {\left( {{{J}}_{\rm{j}}^{{t}}} \right)^{ - {\rm T}}}{{{K}}_{\rm{j}}}{\left( {{{J}}_{\rm{j}}^{{t}}} \right)^{ - 1}}.$

大多数情况下，关节变形所带来的雅克比矩阵的变化对于计算的影响较小，因此式（23）适用于绝大部分机械臂的刚度场建立. 对于本文分析的多轴机床，应当遵循式（21）来建立关节刚度场.

令式（21）中等式右边第一项为 ${{{K}}_{\rm{c}}}$，则 ${{{K}}_{\rm{c}}}$ 可表示为

(24) ${{{K}}_{\rm{c}}} \!=\! \left[ { {\displaystyle\frac{{\partial {{({{J}}_{\rm{j}}^{{t}})}^{\rm T}}}}{{\partial {q_1}}}{{F}}},\;{\displaystyle\frac{{\partial {{({{J}}_{\rm{j}}^{{t}})}^{\rm T}}}}{{\partial {q_2}}}{{F}}},\;{\displaystyle\frac{{\partial {{({{J}}_{\rm{j}}^{{t}})}^{\rm T}}}}{{\partial {q_3}}}{{F}}},\;{\displaystyle\frac{{\partial {{({{J}}_{\rm{j}}^{{t}})}^{\rm T}}}}{{\partial {q_4}}}{{F}}},\;{\displaystyle\frac{{\partial {{({{J}}_{\rm{j}}^{{t}})}^{\rm T}}}}{{\partial {q_5}}}{{F}}} } \right].$

结合式（21）、（24），可得机床关节的刚度矩阵：

(25) ${{K}}_{\rm{j}}^{{t}} = {\left( {{{J}}_{\rm{j}}^{{t}}} \right)^{ - {\rm T}}}({{{K}}_{\rm{j}}} - {{{K}}_{\rm{c}}}){\left( {{{J}}_{\rm{j}}^{{t}}} \right)^{ - 1}}.$

在多轴机床中，由于Z轴和A轴悬伸较长，产生的变形对于末端的变形来说不可忽略，将Z轴和A轴视为柔性运动轴，建立其相应的刚度矩阵.

由于各个运动轴变形属于弹性小变形，可忽略单个运动轴变形对坐标变换矩阵的影响，将运动轴变形所产生的机器人末端变形等效为各个运动轴变形的线性叠加^[18]. 基于这个假设，以柔性运动轴Z轴为例，其末端刚度矩阵推导如下.

如图9所示，考虑固连在柔性轴Z轴末端的坐标系 ${O_Z}$、多轴机床的末端坐标系 ${O_t}$ 及基坐标系 ${O_r}$，作用在末端坐标系 ${O_t}$ 上的广义力F需要变换为坐标系 ${O_Z}$ 内的广义力 ${{{F}}_{\rm{Z}}}$，其转化公式为

(26) ${{{F}}_{{Z}}} = \left[ {\begin{array}{*{20}{c}} {{{R}}_{{t}}^{{Z}}}&0 \\ 0&{{{R}}_{{t}}^{{Z}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{I}}&0 \\ {({{P}}_{{t}}^{{Z}}) \times }&{{I}} \end{array}} \right]{{F}} = {{J}}_{{\bf{F{{\rm or}}}}}^{{{{t}},{ Z}}}{{F}}.$

式中： $({{P}}_{\rm{t}}^{{Z}}) \times $ 表示两坐标系原点连线向量 ${{P}}_{{Z}}^{\rm{t}}$ 与F之间的外积； ${{R}}_{\rm{t}}^{{Z}}$ 表示坐标系 ${O_Z}$ 与 ${O_t}$ 之间的旋转变换； ${{J}}_{{{{\rm F}{\rm or}}}}^{{{t},{ Z}}}$ 为坐标系 ${O_Z}$ 与 ${O_t}$ 之间力转换关系矩阵.

在基坐标系 ${O_r}$ 中，有如下的公式成立：

(27) $ \left. {\begin{array}{*{20}{l}} {{{v}}_{{Z}}^{{r}} = {{v}}_{{t}}^{{r}} + {{w}}_{{t}}^{{r}} \times ({{R}}_{{t}}^{{r}}{{P}}_{{Z}}^{{t}}),}\\ {{{w}}_{{Z}}^{{r}} = {{w}}_{{t}}^{{r}}.} \end{array}} \right\} $

式中： ${{v}}_{{Z}}^{{r}}$ 为坐标系 ${O_Z}$ 在基坐标系内的线速度向量； ${{v}}_{{t}}^{{r}}$ 为坐标系 ${O_t}$ 在基坐标系内的线速度向量； ${{R}}_{{t}}^{{r}}$ 为坐标系 ${O_t}$ 与基坐标系之间的旋转变换矩阵； ${{w}}_{{Z}}^{{r}}$ 坐标系 ${O_Z}$ 在基坐标系内的角速度向量； ${{w}}_{{t}}^{{r}}$ 为坐标系 ${O_t}$ 在基坐标系内的角速度向量.

由式（27）可得到坐标系 ${O_Z}$ 与坐标系 ${O_t}$ 在基坐标系内的位移变换关系为

(28) ${{\delta v}}_{{Z}}^{{r}} = \left[ {\begin{array}{*{20}{c}} {{I}}& - {{{(R}}_{{t}}^{{r}}{{P}}_{{Z}}^{{t}}{{)}} \times } \\ {{0}}&{{I}} \end{array}} \right]{{\delta v}}_{{t}}^{{r}}= {{ J}}_{{\rm{Dis}}}^{{{t},{ Z},{ r}}}{{\delta v}}_{{t}}^{{r}}.$

根据机器人学相关理论，存在以下公式成立：

(29) $\left. {\begin{array}{*{20}{l}} {{{\delta }}{{{v}}_{{Z}}} = {{J}}_{{Z}}^{{r}}{{\delta v}}_{{Z}}^{{r}}} ,\\ {{{\delta }}{{{v}}_{{t}}} = {{J}}_{{t}}^{{r}}{{\delta v}}_{{t}}^{{r}}}. \end{array}} \right\}$

式中： ${{\delta }}{{{v}}_{{Z}}}$ 为坐标系 ${O_Z}$ 的微小位移在自身坐标系内表示； ${{\delta }}{{{v}}_{{t}}}$ 坐标系 ${O_t}$的微小位移在自身坐标系内表示； ${{J}}_{{Z}}^{{r}}$ 为雅克比矩阵，表示坐标系 ${O_Z}$ 与 ${O_r}$之间的微分位移关系； ${{J}}_{\rm{t}}^{\rm{r}}$ 为雅克比矩阵，表示坐标系 ${O_r}$ 与 ${O_t}$之间的微分位移关系； ${{J}}_{{\rm{Dis}}}^{{{t},{ Z}}}$ 为表示坐标系 ${O_Z}$ 与坐标系 ${O_t}$ 之间微分位移关系的矩阵.

联立式（28）、（29），可得到坐标系 ${O_Z}$ 与坐标系 ${O_t}$ 在各自坐标系内的位移变换关系为

(30) $ { {{\delta }}{{{v}}_{{Z}}} = {{J}}_{{Z}}^{{r}}{{J}}_{{\rm{Dis}}}^{{{t},{ Z},{ r}}}{\left( {{{J}}_{{t}}^{{r}}} \right)^{ - 1}}{{\delta }}{{{v}}_{{t}}} = {{J}}_{{\rm{Dis}}}^{{{t},{ Z}}}{{\delta }}{{{v}}_{{t}}} } .$

在坐标系 ${O_Z}$ 内，存在如下公式成立：

(31) ${{{F}}_{{Z}}} = {{{K}}_{{Z}}}{{\delta }}{{{v}}_{{Z}}}.$

联立式（26）、（30）（31）即可得到如下公式：

(32) ${{F}} = {\left( {{{J}}_{{\rm{For}}}^{{{t},{ Z}}}} \right)^{ - 1}}{{{K}}_{{Z}}}{{J}}_{{\rm{Dis}}}^{{{t},{ Z}}}{{\delta }}{{{v}}_{{t}}}.$

根据式（32）即可得到末端坐标系内的Z轴刚度矩阵为

(33) ${{K}}_{{Z}}^{{t}} = {\left( {{{J}}_{{\rm{For}}}^{{{t},{ Z}}}} \right)^{ - 1}}{{{K}}_{{Z}}}{{J}}_{{\rm{Dis}}}^{{{t},{ Z}}}.$

采用类似的解法可得到A轴在末端坐标系内的刚度矩阵为

(34) ${{K}}_{{A}}^{{t}} = {\left( {{{J}}_{{\rm{For}}}^{{{t},{ A}}}} \right)^{ - 1}}{{{K}}_{{A}}}{{J}}_{{\rm{Dis}}}^{{{t},{ A}}}.$

式中： ${{J}}_{{\rm{For}}}^{{{t},{ A}}}$ 为坐标系 ${O_A}$ 与坐标系 ${O_t}$ 之间力转换关系矩阵； ${{J}}_{{\rm{Dis}}}^{{{t},{ A}}}$ 为坐标系 ${O_A}$ 与坐标系 ${O_t}$ 之间微分位移关系矩阵.

由式（33）、（34），可得到柔性运动轴的在末端坐标系内的刚度矩阵为

(35) ${{K}}_{{i}}^{{t}} = {\left( {{{J}}_{{\rm{For}}}^{{{t},{ Z}}}} \right)^{ - 1}}{{{K}}_{{Z}}}{{J}}_{{\rm{Dis}}}^{{{t},{ Z}}} + {\left( {{{J}}_{{\rm{For}}}^{{{t},{ A}}}} \right)^{ - 1}}{{{K}}_{{A}}}{{J}}_{{\rm{Dis}}}^{{{t},{ A}}}.$

铆接头末端坐标系即为整个机床的末端坐标系，因此如式（4）所示的刚度矩阵无需变换，即为末端执行器在末端坐标系内的刚度矩阵.

为验证多轴机床总柔度矩阵的准确性，对不同位姿下的多轴机床的末端变形进行试验求取，如图10所示为末端变形测量试验现场. 该实验主要需要卧式双机联合钻铆系统、末端变形实验装置及激光跟踪仪这3个设备来完成. 在具体操作时，卧式双机联合钻铆系统中带内铆接头的多轴机床保持位姿不变，带外铆接头的多轴机床进行加载；末端变形试验装置安装在内铆接头的末端，其结构如图11所示. 该装置上的4个凸起内安装有磁铁，用以固定激光跟踪仪的测量球；激光跟踪仪放置在带内铆接头的多轴机床一侧，并且在合适位置完成找准工作后，需要保持该位置在整个实验过程中不变. 激光跟踪仪作为本试验的测量设备，其测量精度和距离测量球的长度有关，其精度公式如下：

(36) $\sigma = (0.015 + 0.006L_{\rm m}) \times {10^{ - 3}}.$

式中： $\sigma $为误差值，L_m为测量距离. 试验中测量距离为1~2 m，故测量精度为0.021~0.027 mm.

实际试验时，首先卧式双机联合钻铆系统的位姿参数经由其自带的控制台调整到实验所需数值，然后外铆接头一侧的多轴机床进行加载. 在保持一定压铆力的情况下，利用激光跟踪仪分别测量4个测量球在加载前及加载后的空间坐标，经计算后得到多轴机床的末端变形.

试验测得的结果与计算得到的结果如表3所示， ${F_{\rm{e}}}$ 为试验时施加于多轴机床末端的压铆力， ${D_{\rm c}}$为计算得到的多轴机床末端位移， ${D_{\rm e}}$为试验得到的多轴机床末端位移，η为相对偏差. 试验得到的末端变形结果与计算得到的末端变形很接近，其差异值不超过10%，说明所建模型具有相当高的精度，与真实情况较为符合，利用该模型计算所得的多轴机床末端变形具有相当高的可信度.

（1）在其余位置参数不变时，变换X轴位置参数对末端变形量的影响并不大；变换Y轴或者Z轴位置参数时，末端变形量变化较为明显. 因此，后续的结构优化工作可重点关注Y轴与Z轴，以提高自动钻铆设备的末端刚度.

（2）在钻铆力一定时，试验所得的末端变形量最大值与最小值的差异约为32%，表明该多轴机床具有良好的整体刚度，其末端刚度对位姿变化的影响不太敏感.

（3）利用本文所述末端综合柔度矩阵所得到的机床末端变形量与实际实验所得结果误差不超过10%，说明该末端综合柔度矩阵具有相当高的精度，并且在钻铆系统受重载的情况下仍适用. 这对于卧式双机联合钻铆系统的设计与后续的改进具有重要的参考意义.