随着空间技术的快速发展，大尺寸、轻量化的空间机械臂已经被广泛地用于空间站的组装、修理和检测，航天员舱外活动，空间飞行器捕获等，其可靠性对于整个空间任务至关重要。空间机械臂动力学分析涉及运动过程中大变形和间隙碰撞问题，且结构参数存在不确定性，因此，要建立精确的空间机械臂动力学模型就需要考虑各种因素。

Flores等^[1]对多种间隙模型进行了研究，通过对含间隙机构进行动力学分析发现间隙严重影响机构动力学特性。Zhao等^[2]对含间隙空间机器人进行了动力学特性研究，研究发现间隙尺寸会显著影响机器人动力学特性。Erkaya^[3]建立了含间隙焊接机械臂的动力学模型。Varedi-Koulaei等^[4]采用粒子群算法对含间隙的3-RRR并联机械手进行了优化设计，以减小间隙碰撞对该机械手动力学特性的影响。Megahed等^[5]对含多间隙机构进行了研究，结果表明间隙位置，运行速度和间隙大小都会显著影响机构动力学特性。谷勇霞等^[6]研究了间隙尺寸和间隙关节数量对漂浮空间机械臂动力学特性的影响。考虑到构件柔性，Khemili等^[7]对含间隙平面柔性机构进行了动力学分析，结果显示构件柔性能够有效控制碰撞力峰值。任武等^[8]建立了考虑间隙的多柔体臂架的刚柔混合模型。Ben Abdallah等^[9]对含多间隙的柔性机构进行了研究。考虑到机械臂在运动过程中存在大变形，Tian等^[10]基于绝对节点坐标法^[11]对含PD控制器的双连杆柔性机械臂进行了控制精度分析。然而，在建模过程中忽略了间隙对机械臂控制精度的影响。

上述研究都认为机械臂为理想情况，然而由于材料不均匀、加工误差、摩擦、磨损等，真实的机械臂中普遍存在不确定性。Li等^[12]研究了存在随机不确定性的含间隙机械臂动力学特性和运动精度问题，建模过程中忽略了机械臂存在的变形。Wang等^[13]采用Chebyshev多项式方法对存在区间不确定性参数的刚柔耦合空间机械臂进行了动力学分析，建模过程中忽略了间隙碰撞。笔者将系统考虑大变形、间隙碰撞和不确定性对机械臂动力学特性和控制精度的影响，研究含区间铰间隙的柔性机械臂动力学建模和控制精度分析方法。

1 旋转铰间隙建模 1.1 旋转铰间隙数学模型

如图 1所示，体i, j之间存在平面间隙转动副，轴承是体i的部件，轴是体j的部件。整体坐标系为xy。点P_i和P_j分别为轴承和轴的中心，对应的位置矢量分别为r_i和r_j。轴与轴承的中心距矢量e_ij可表示为

$ {\mathit{\boldsymbol{e}}_{ij}} = {\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}。$

(1)

轴与轴承中心距为

$ {e_{ij}} = \sqrt {\mathit{\boldsymbol{e}}_{ij}^{\rm{T}}{\mathit{\boldsymbol{e}}_{ij}}} 。$

(2)

垂直轴承与轴接触表面的单位法向量为

$ {\mathit{\boldsymbol{n}} = \frac{{{\mathit{\boldsymbol{e}}_{ij}}}}{{{e_{ij}}}}}。$

(3)

当轴承与轴碰撞时，嵌入深度为

$ {\delta = {e_{ij}} - c,} $

(4)

式中c为间隙大小，c=R_i-R_j，R_i和R_j分别为轴承和轴半径。

体i, j的接触点分别为Q_i和Q_j，基于绝对节点坐标法接触点的全局位置为

$ \mathit{\boldsymbol{r}}_k^Q = {\mathit{\boldsymbol{r}}_k} + {R_k}\mathit{\boldsymbol{n}}\quad (k = i,j)。$

(5)

全局坐标系下接触点Q_i和Q_j的速度可以通过将式(5)对时间求导得到

$ \mathit{\boldsymbol{\dot r}}_k^Q = {{\mathit{\boldsymbol{\dot r}}}_k} + {R_k}\mathit{\boldsymbol{n}}\quad (k = i,j)。$

(6)

接触点的法向相对速度v_N和切向相对速度v_T可表示为

$ {{v_{\rm{N}}} = {{\left( {\mathit{\boldsymbol{\dot r}}_j^Q - \mathit{\boldsymbol{\dot r}}_i^Q} \right)}^{\rm{T}}}\mathit{\boldsymbol{n}},} $

(7)

$ {{v_{\rm{T}}} = {{\left( {\mathit{\boldsymbol{\dot r}}_j^Q - \mathit{\boldsymbol{\dot r}}_i^Q} \right)}^{\rm{T}}}\mathit{\boldsymbol{t}},} $

(8)

式中切向量t通过向量n逆时针旋转90°得到。

1.2 接触碰撞力模型

Hertz接触力模型是应用最广泛的接触力模型，该模型认为接触表面为完全弹性碰撞，因此不能描述碰撞过程中的能量耗散。为了描述接触碰撞过程中的能量损耗，基于Hertz模型，Lankarani，Herbert和Lee等提出了多种含能量耗散项的接触碰撞力模型^[14]。Liu等^[15]研究发现Hertz接触力模型仅仅适用于间隙相对较大且载荷较小的情况。对于小间隙情况，计算结果将无法满足精度要求。针对上述问题，Bai等^[16]构建了一种基于Lankarani-Nikravesh碰撞力模型和WEF(Winkler elastic foundation)模型的混合碰撞力模型，该碰撞力模型可以表示为

$ {F_{\rm{n}}} = {K_{\rm{n}}}{\delta ^n} + {D_{{\rm{mod}}}}\dot \delta , $

(9)

式中：F_n为法向碰撞力， $\dot{\delta}$为撞击点的相对速度。n的值取决于碰撞面的材料，对于金属材料通常选取n=1.5。K_n和D_mod分别是非线性接触刚度系数和修改后的阻尼系数，其表达式如下

$ {{K_{\rm{n}}} = \frac{1}{8}{\rm{ \mathsf{ π} }}{E^*}\sqrt {\frac{{2\delta {{\left( {3{R_i} - 3{R_j} + 2\delta } \right)}^2}}}{{{{\left( {{R_i} - {R_j} + \delta } \right)}^3}}}} ,} $

(10)

$ {{D_{{\rm{mod}}}} = \frac{{3{K_{\rm{n}}}\left( {1 - c_{\rm{e}}^2} \right){{\rm{e}}^{2\left( {1 - {c_{\rm{e}}}} \right)}}{\delta ^n}}}{{4{{\dot \delta }^{( - )}}}},} $

(11)

式中： $\dot{\delta}$^(-)为撞击点的初始相对速度，c_e为恢复系数，E^*是复合弹性模量，可表示为如下形式

$ \frac{1}{{{E^*}}} = \frac{{1 - v_i^2}}{{{E_i}}} + \frac{{1 - v_j^2}}{{{E_j}}}, $

(12)

式中：v_k(k=i, j)和E_k(k=i, j)分别是接触材料的泊松比和弹性模量。该间隙碰撞力模型能够普遍适用于大间隙和小间隙情况，通过与有限元仿真结果比较，该碰撞力模型的精度高于Hertz接触力模型。

1.3 摩擦力模型

Coulomb摩擦模型是最简单的摩擦力模型。该模型的摩擦力方向与切向速度方向相反，并且摩擦力大小与法向接触力成正比。Coulomb摩擦模型在切向速度趋近于零时会产生不连续的摩擦力，这将有可能导致计算不收敛。为了解决该问题，Ambrósio^[17]将速度引入Coulomb摩擦模型，建立了如下形式的摩擦力模型，

$ {F_{\rm{T}}} = - {c_{\rm{f}}}{c_{\rm{d}}}{F_{\rm{N}}}\frac{{{\mathit{\boldsymbol{v}}_{\rm{T}}}}}{{{v_{\rm{T}}}}}, $

(13)

式中：c_f为摩擦系数，v_T为相对切向速度。c_d为动态修正系数，可表示为

$ {c_{\rm{d}}} = \left\{ {\begin{array}{*{20}{l}} 0&{{\rm{ if }}\quad {v_{\rm{T}}} \le {v_0},}\\ {\frac{{{v_{\rm{T}}} - {v_0}}}{{{v_1} - {v_0}}}}&{{\rm{ if }}\quad {v_0} \le {v_{\rm{T}}} \le {v_1},}\\ 1&{{\rm{ if }}\quad {v_{\rm{T}}} \ge {v_1},} \end{array}} \right. $

(14)

式中：v₀和v₁是为计算动态修正系数给定的速度极限值。Ambrósio摩擦力模型如图 2所示。

2 含间隙柔性机械臂动力学建模

含间隙柔性机械臂模型如图 3所示。连杆OA_i在O点处的运动副为理想平面转动副，连杆OA_i与连杆A_jP_j在A处的运动副为平面间隙转动副。τ₁和τ₂分别为作用于连杆OA_i和连杆A_jP_j的驱动力矩，n₁和t₁为连杆OA_i在O点处的单位法向量和切向量，n₂和t₂为连杆OA_i在A_i点处的单位法向量和切向量，n₃和t₃为连杆A_jP_j在A_j点处的单位法向量和切向量，θ₁为连杆OA_i在O点处的转角，θ₃为连杆OA_i在A_i点处的转角，θ₂为连杆A_jP_j在A_j点处的转角。考虑到运动过程中连杆OA_i和连杆A_jP_j存在大变形，采用绝对节点坐标法对连杆OA_i和连杆A_jP_j进行动力学建模。

基于绝对节点坐标法^[9]，平面梁单元上任意一点全局位置可表示

$ \mathit{\boldsymbol{r}} = \mathit{\boldsymbol{S}}{\mathit{\boldsymbol{e}}^k}, $

(15)

式中S为由参数x决定的单元形函数，参数x∈[0  l]表示未变形时特定点的轴向坐标，l为单元初始长度，单元形函数表示为：

$ \mathit{\boldsymbol{S}} = \left[ {\begin{array}{*{20}{l}} {{S_1}{\mathit{\boldsymbol{I}}_2}}&{{S_2}{\mathit{\boldsymbol{I}}_2}}&{{S_3}{\mathit{\boldsymbol{I}}_2}}&{{S_4}{\mathit{\boldsymbol{I}}_2}} \end{array}} \right], $

(16)

式中，I₂为2×2的单位矩阵； $S_{1}=1-3\left(\frac{x}{l}\right)^{2}+2\left(\frac{x}{l}\right)^{3}$, $S_{2}=l\left[\frac{x}{l}-2\left(\frac{x}{l}\right)^{2}+\left(\frac{x}{l}\right)^{3}\right]$, $S_{3}=3\left(\frac{x}{l}\right)^{2}-2\left(\frac{x}{l}\right)^{3}$, $S_{4}=l\left[\left(\frac{x}{l}\right)^{3}-\left(\frac{x}{l}\right)^{2}\right]$；e^k为单元节点的广义坐标矢量，表示为

$ {\mathit{\boldsymbol{e}}^k} = {\left[ {\mathit{\boldsymbol{r}}_m^{\rm{T}}\quad {{\left( {\frac{{\partial {\mathit{\boldsymbol{r}}_m}}}{{\partial x}}} \right)}^{\rm{T}}}\quad \mathit{\boldsymbol{r}}_n^{\rm{T}}\quad {{\left( {\frac{{\partial {\mathit{\boldsymbol{r}}_n}}}{{\partial x}}} \right)}^{\rm{T}}}} \right]^{\rm{T}}}。$

(17)

r_i(i=m, n)表示节点i(i=m, n)位置矢量，∂r_i/∂x(i=m, n)为节点i(i=m, n)的位置矢量对x方向的二维偏导数矢量。

基于上述描述，梁单元的动能可表示为

$ {T_{\rm{e}}} = \frac{1}{2}\int_V \rho {{\mathit{\boldsymbol{\dot r}}}^{\rm{T}}}\mathit{\boldsymbol{\dot r}}{\rm{d}}V = \frac{1}{2}{\left( {{{\mathit{\boldsymbol{\dot e}}}^k}} \right)^{\rm{T}}}{{\bf{M}}_{\rm{e}}}{{\mathit{\boldsymbol{\dot e}}}^k}, $

(18)

式中M_e为梁单元的常数质量矩阵，可通过式(19)计算得到，

$ {\mathit{\boldsymbol{M}}_{\rm{e}}} = \frac{1}{2}\int_V \rho {\mathit{\boldsymbol{S}}^{\rm{T}}}\mathit{\boldsymbol{S}}{\rm{d}}V, $

(19)

式中：ρ和V分别为单元密度和体积。

梁单元的弹性势能可表示为

$ {U_{\rm{e}}} = \frac{1}{2}\int_V {{\mathit{\boldsymbol{\sigma }}^{\rm{T}}}} \mathit{\boldsymbol{\varepsilon }}{\rm{d}}V, $

(20)

式中：σ和ε分别为应力矢量和应变矢量。

驱动力矩τ₁和τ₂所做的虚功可以表示为

$ \delta W = {\tau _1}\delta {\theta _1} + {\tau _2}\delta {\theta _4} = {\tau _1}\delta {\theta _1} + {\tau _2}\delta {\theta _2} - {\tau _2}\delta {\theta _3}。$

(21)

对于一维两节点柔性梁θ₁的变分可以表示为

$ \delta {\theta _1} = \frac{{\partial r_{Ox}^i/\partial x\delta \left( {\partial r_{Oy}^i/\partial x} \right) - \partial r_{Oy}^i/\partial x\delta \left( {\partial r_{Ox}^i/\partial x} \right)}}{{f_1^2}}{\rm{ , }} $

(22)

式中： $f_{1}=\sqrt{\left(\partial r_{O x}^{1} / \partial x\right)^{2}+\left(\partial r_{O y}^{1} / \partial x\right)^{2}}$， $\left[\partial r_{0 x}^{1} / \partial x \quad \partial r_{0 y}^{1} / \partial x\right]^{\mathrm{T}}$为O点处位置矢量对x方向的二维偏导数矢量，则驱动力矩τ₁作用在O点的等效广义力表达为：

$ {\left. {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{Q}}_1} = \left[ {\begin{array}{*{20}{l}} 0&0&{ - {\tau _1}\left( {\partial r_{Oy}^1/\partial x} \right)/f_1^2} \end{array}} \right.}&{{\tau _1}\left( {\partial r_{Ox}^1/\partial x} \right)/f_1^2}&0&0&0&0 \end{array}} \right]^{\rm{T}}}。$

(23)

同理可以得到τ₂的等效广义力。基于虚功原理和Baumgarte方法^[18]，该含间隙柔性机械臂的动力学方程表示为

$ \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{M}}&{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_\mathit{\boldsymbol{q}}^{\rm{T}}}\\ {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_\mathit{\boldsymbol{q}}}}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\ddot q}}}\\ \mathit{\boldsymbol{\lambda }} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{B\tau }} - \mathit{\boldsymbol{F}} + {\mathit{\boldsymbol{F}}_{\rm{c}}}}\\ {\gamma - 2\alpha \mathit{\boldsymbol{ \boldsymbol{\dot \varPhi} }} - {\beta ^2}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} \end{array}} \right], $

(24)

式中：M为机械臂的常数质量矩阵，λ为拉格朗日乘子，F和F_c分别为广义弹性力和广义间隙碰撞力，Φ和Φ_q分别为约束方程及其对广义坐标q的偏导数矩阵， $\gamma=-\left(\boldsymbol{\mathit{\Phi}}_{q} \dot{\boldsymbol{q}}\right)_{q}-2 \boldsymbol{\mathit{\Phi}}_{q t} \dot{\boldsymbol{q}}-\boldsymbol{\mathit{\Phi}}_{t t}$，α和β为稳定性系数，其值依靠经验选取。Bτ为驱动力矩的等效广义力，τ=[τ₁  τ₂]^T由PID控制方法^[19]得到。

3 含区间铰间隙柔性机械臂控制精度和动力学特性分析

由于材料不均匀、加工误差、磨损等原因，机械臂广泛存在不确定性，如间隙尺寸、部件杨氏模量等。这些不确定性因素会显著的影响机械臂动力学特性和控制精度，考虑到间隙尺寸和杨氏模量数据的不完备性，这里将间隙尺寸和杨氏模量处理为区间变量，式(24)变为

$ \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{M}}&{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_\mathit{\boldsymbol{q}}^{\rm{T}}}\\ {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_\mathit{\boldsymbol{q}}}}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot {\tilde q}}}}\\ {\mathit{\boldsymbol{\tilde \lambda }}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{B\tilde \tau }} - \mathit{\boldsymbol{\tilde F}} + {{\mathit{\boldsymbol{\tilde F}}}_{\rm{c}}}}\\ { - 2\alpha \mathit{\boldsymbol{\dot { \boldsymbol{\tilde \varPhi} }}} - {\beta ^2}\mathit{\boldsymbol{ \boldsymbol{\tilde \varPhi} }}} \end{array}} \right], $

(25)

式中上标表示该值与区间变量(c, E)有关。式(25)为含区间不确定性变量和控制器的微分代数方程，对其进行求解将非常困难。这里将提出该动力学模型的求解方法，图 4列出了该方法的主要步骤，具体为：

1) 选择区间不确定性的抽样数量N。

2) 从每一个区间参数中选择一个样本(c_i, E_i)。

3) 在给定区间参数(c_i, E_i)的情况下，通过计算动力学方程(24)得到柔性机械臂末端控制误差e_i(t)=(x_is(t)-x_id(t)，y_is(t)-y_id(t))((x_is, y_is)和(x_id, y_id)分别为实际和理想情况末端坐标)、末端加速度 $\left(\ddot{x}_{i \mathrm{s}}, \ddot{y}_{i \mathrm{s}}\right)$和驱动力矩(τ_i1, τ_i2)。

4) 判断区间不确定性的N个样本是否已经全部分析。如果否，回到步骤(2)；如果是，继续步骤(5)。

5) 找出末端控制误差、末端加速度和驱动力矩上边界和下边界，该上边界和下边界围成的区域即为考虑区间不确定性时末端控制误差、末端加速度和驱动力矩可能的区域。

4 仿真分析

该含间隙柔性机械臂末端点进行圆轨迹跟踪控制，如图 5所示。连杆OA_i和连杆A_jP_j长度分别为0.4 m和1.6 m。连杆OA_i划分为2个单元，连杆A_jP_j划分为4个单元。连杆OA_i和A_jP_j具有相同的截面形状和材料参数，截面积为2.5×10^-3 m²，截面惯性矩为5.208×10^-7 m⁴，杨氏模量为30 GPa，密度为6 000 kg/m³。该机械臂的材料属性及动力学分析相关参数如表 1所示。仿真时间为10 s，跟踪圆轨迹的角速度为

$ \omega = - \frac{{{\rm{6 \mathsf{ π} }}}}{{500}}\left( {{t^2} - 10t} \right)。$

(26)

为了研究间隙对机械臂动力学特性和控制精度的影响，分别取间隙c=0.01 mm，c=0.2 mm和c=0.5 mm，摩擦系数c_f=0.1，将上述多种间隙情况计算得到的响应结果与理想铰接副的响应结果进行比较。图 6为3种间隙情况下的末端控制误差图、末端加速度图和驱动力矩图。由图 6可知：当间隙大小为0.01 mm时，间隙对末端控制误差、末端加速度和驱动力矩几乎没有影响；当间隙大小为0.2 mm和0.5 mm时，间隙对末端控制误差的影响非常显著，而且间隙越大控制精度越差；当间隙大小为0.2 mm和0.5 mm时，间隙对末端加速度和驱动力矩的影响依然不显著。通过对图 7(d)进一步分析发现初始时刻机械臂末端y方向的加速度扩展较大，其原因在于：初始时刻间隙铰中的轴与轴承存在剧烈碰撞，导致初始时刻机械臂末端y方向的加速度较大。

为了分析部件柔性对机械臂动力学特性和控制精度的影响，取连杆杨氏模量为3，30，300 GPa，3种工况下的间隙c=0.5mm。末端控制误差、末端加速度和驱动力矩如图 7所示。可以发现：当杨氏模量为3 GPa时，末端控制误差、末端加速度和驱动力矩与杨氏模量为30 GPa的分析结果存在明显差异，而且弹性模量较低时，末端剧烈抖动持续时间将增加；当杨氏模量为30 GPa和300 GPa时，两种工况的末端控制误差、末端加速度和驱动力矩差异不显著。因此，当连杆弹性模量较小时，连杆柔性会显著的影响控制误差和机械臂的动力学特性，而连杆弹性模量较大时，连杆柔性对该机械臂的控制精度和动力学特性的影响将不显著。

由于间隙尺寸和连杆柔性会显著影响该机械臂控制精度和动力学特性，这里考虑间隙尺寸和杨氏模量为区间变量，其取值范围分别为c=[0.1   0.5mm]，E=[20   30]GPa。采用本文中提出的方法对机械臂进行动力学特性和控制精度分析。末端控制误差、末端加速度和驱动力矩如图 8所示。可以发现：考虑不确定性时，末端控制误差将扩张，从而导致控制精度降低；末端加速度和驱动力矩也会扩展，从而影响该机械臂的动力学特性。

5 结语

1) 当间隙尺寸较小时，间隙对末端控制误差、末端加速度和驱动力矩几乎没有影响。当间隙尺寸较大时，间隙对末端控制误差的影响非常显著，而且间隙越大控制精度越差。然而，较大间隙情况下，间隙对末端加速度和驱动力矩的影响依然不显著。

2) 部件柔性会显著的影响末端控制误差、末端加速度和驱动力矩。当部件弹性模量较低时，末端会剧烈抖动。

3) 参数不确定性会显著的降低该机械臂的控制精度，而且也会影响该机械臂的动力学特性。