马曙光，于 强，王培栋，王天舒

(清华大学 航天航空学院·北京·100084)

0 引 言

在轨服务[1]是近年来航天领域的研究热点之一。通过在轨服务可以对航天器执行加注燃料、更换零件、拆解组装和轨道转移等操作，以帮助在轨航天器保持良好工作状态和延长使用寿命，从而降低航天工程的经济成本。空间机械臂由于其技术相对成熟和通用性较强等优点，是最常用的在轨服务执行机构之一。关于空间机械臂的工程验证项目包括日本的“工程实验卫星-7”(ETS-VII)[2]和美国的“轨道快车”[3]等。与依靠关节电机驱动的机械臂相比，绳驱机械臂[4-5]的尺寸可以设计得更小。本文展示一种绳驱空间柔性机械臂，该机械臂由多根臂杆串联而成，可灵活改变形状，适合在狭小空间内进行精细在轨服务操作。

常见的机械臂动力学建模方法包括虚拟臂法[6]、等效臂法[7]和多体动力学方法[8-10]等。其中，虚拟臂法和等效臂法属于近似方法，其精确程度不如多体动力学方法。对于本文要讨论的绳驱空间机械臂，由于其结构较为复杂且自由度较高，近似方法无法精确分析其动力学特性，多体动力学方法例如Lagrange 方程、R-W方法、Kane方程和递推算法[9-10]可能会更加合适。在上述动力学方法中，递推算法由于计算效率较高(时间复杂度可达O(N)或O(logN))，比较适合自由度较高的多体系统。本文采用递推算法中的一种即空间算子代数方法[10]来对绳驱空间机械臂进行动力学建模。

数值仿真效率往往由单个时间步的计算量和时间步长大小来决定，对于本文所讨论的动力学模型而言，其数值仿真效率主要取决于时间步长大小。在多体动力学模型中，通常将动力学方程写成微分方程的形式，然后再利用微分方程的数值算法[11]进行求解。对于常见的微分方程数值算法(例如Runge-Kutta法)，由于存在稳定区域[11]对时间步长的限制，时间步长过大时仿真结果会很快发散，所以时间步长的选取存在临界值，而该临界值将限制动力学模型的数值仿真效率。本文通过对“稳定区域”的计算来分析动力学模型中影响时间步长临界值的因素，并对动力学模型做出相应改进以增大时间步长临界值，从而提高动力学模型的数值仿真效率。

1 动力学模型

1.1 机械结构

如图1所示，绳驱空间机械臂的主体机械结构分为根部、中段和端部三部分。其中，根部用于连接机械臂和漂浮基座；中段主要包含五组臂杆，每组臂杆由六根相同的臂段串联而成，全部臂段通过万向节进行相连；端部用于放置机械臂的载荷。除30根臂杆之外，机械臂还包括15根驱动绳索和60根联动绳索。如图2所示，驱动绳索穿过臂杆的过孔并带动臂杆进行运动。如图3所示，对于同一组臂段而言，其组内的万向节由联动绳索来保证转角的同步性，该同步性使得机械臂在运动时具有呈分段等曲率的几何特征，从而增强机械臂的操控性能。

图1 空间机械臂结构示意图Fig.1 Structure of the space manipulator

图2 驱动绳索示意图Fig.2 Driving ropes of the space manipulator

图3 联动绳索示意图Fig.3 Connecting ropes of the space manipulator

1.2 动力学方程

如图4所示，机械臂可视为链式多刚体系统，其动力学方程可表示为

(1)

式(1)中的广义坐标q包含漂浮基座的6个自由度和机械臂的60个转动自由度，因此q的维数为66。

式(1)中的广义力Q由三部分构成：

Q=Qc+Qd+Qe

(2)

式(2)中的Qe为漂浮基座和空间机械臂受到的外力，可通过给定的边界条件计算得出；Qc为联动绳索的等效作用力，Qd为驱动绳索的等效作用力，二者需要建立相应的绳索模型来计算得出。

图4 空间机械臂动力学模型示意图Fig.4 The model of space manipulator as a rigid multibody system

对于联动绳索，可将其中的一对联动绳索等效为扭簧：

(3)

式(3)中的ΔT为扭簧的等效扭矩，kc为扭簧的等效弹性系数，cc为扭簧的等效阻尼系数，Δθ为转角之差，式(3)的推导可通过将联动绳索类比皮带传动得到，具体的推导过程在此略去。

将式(3)同运动学方程结合即可计算出联动绳索的等效作用力Qc。

对于驱动绳索，可将其中某一根绳索视为由臂杆过孔分割的多段弹性绳连接而成，同一根绳索相邻两段的张力大小关系为

Tj+1=Tjesjθj,sj=-μjsgn(vj)

(4)

式(4)中的μj为臂杆过孔处的摩擦系数，θj为两根绳索之间的角度差，vj为绳索相对过孔的运动速度大小。

将式(4)同运动学方程结合即可计算得到所有驱动绳索的张力分布规律，进而得到驱动绳索的等效作用力Qd。

当Qc、Qd和Qe都得到之后，可将式(2)代入式(1)可以得出空间机械臂的动力学方程

(5)

式(5)是空间机械臂动力学分析的基础。需要说明的是，在求解式(5)时，无需计算出质量矩阵M后求逆，而是利用空间算子代数计算中间变量后利用递推方法求解。关于空间算子代数的算法细节可参考相关文献[10]。

式(5)在形式上可以改写为

(6)

式(6)代表一组常微分方程，可用相关数值算法例如Runge-Kutta法进行求解。本文采用四阶Runge-Kutta法来求解式中的微分方程，图5给出了相应的数值仿真流程。

图5 数值仿真流程示意图Fig.5 The scheme of the numerical simulation

2 数值仿真效率分析

动力学模型的数值仿真效率由单个时间步的计算量和时间步长大小来共同决定。在图5所示的仿真流程中，每一时间步需要重新计算各驱动绳索的张力分布规律，该部分计算量占据单个时间步总计算量的2/3左右，而且这些计算量很难通过模型简化来避免。对于本文所讨论的动力学模型而言，其数值仿真效率主要取决于时间步长大小。

图6 各阶Runge-Kutta法稳定区域Fig.6 Stability regions for the explicit Runge-Kutta methods

2.1 限制时间步长的因素

基于式(5)和式(6)，可利用摄动法来推导机械臂在平衡状态下做微振动的动力学方程

(7)

式(7)中，Mv是质量矩阵，Kv是刚度矩阵，T是系统总动能，E是绳索总的弹性势能。

式(7)中各系数的值可由动力学模型给出，其具体推导过程在这里略去。可由式(7)计算得出微振动的最高阶频率为1217Hz。式(7)可以改写为

(8)

式(8)是一阶线性微分方程组，其系数矩阵对应的特征根满足|λI-A|=0。

由图6可知，对于四阶Runge-Kutta法，系数λ和时间步长h需满足λh<2.8i，可得时间步长h与最高阶频率fmax的关系为

(9)

将fmax=1217代入到式(9)中，可计算得出时间步长临界值的理论值约为3.8×10-4s，这与通过数值仿真实验得到的经验值4×10-4s很接近。

从上述分析中，可以看出时间步长临界值主要由绳驱空间机械臂微振动的最高阶频率所决定，而该最高阶频率与振动方程中的质量矩阵和刚度矩阵有关。对于绳驱空间机械臂而言，质量矩阵主要由质量分布所决定，而刚度矩阵则由绳索弹性所决定。由于质量分布不太受动力学模型的影响，因此时间步长临界值主要取决于刚度矩阵，即影响仿真效率的主要因素为动力学模型中的绳索模型。

2.2 动力学模型的改进

通过2.1的分析可以得知时间步长的临界值主要由绳索弹性所决定。在1.2中给出的动力学模型中，联动绳索被等效成了扭簧，由于扭簧的等效刚度较高，导致最高阶振动频率较高，时间步长临界值较小。为增大时间步长的临界值，可对联动绳索的模型进行改进。

将联动绳索视为几何约束，此时系统的动力学方程可改写为

(10)

式(10)中的R(q)=0为联动绳索对应的几何约束。

基于式(10)可利用摄动法来推导空间机械臂在平衡状态下做微振动的动力学方程。计算结果给出式(10)的最高阶频率值为223Hz，对应的时间步长临界值的理论值为2.0×10-3s，这与对应的经验值2×10-3s相符。对于式(10)所代表的动力学模型，时间步长临界值增大了5倍，即仿真效率最多可以提高5倍。

为了进一步提高仿真效率，还可以对驱动绳索的模型进行改进。将驱动绳索视为刚性绳，此时动力学模型中不存在弹性项，理论上将不存在时间步长临界值，时间步长的选取将由计算精度等因素来决定。

如果将联动绳索视为几何约束，同时将驱动绳索视为刚性绳，数值仿真实验显示时间步长可以取至0.1s且仿真结果未出现发散现象。与1.2中的动力学模型相比，改进后的动力学模型仿真效率提高了200倍以上。特别的，当时间步长取0.1s时，动力学模型可以实现实时仿真。

3 结 论

本文建立了一种绳驱空间机械臂的动力学模型，并在此基础上研究了影响数值仿真效率的因素。通过分析微分方程数值算法的稳定区域，发现仿真步长的临界值主要取决于动力学模型中的绳索弹性。基于上述结论对动力学模型进行改进，仿真算例表明改进后的动力学模型相比最初的动力学模型仿真效率可提高200倍以上，并可以实现实时仿真，从而增加了该动力学模型的工程实际应用价值。