登月助力航天服下肢关节迟滞模型建立与动力学仿真

李照阳,戴跃洪,唐 鹏,石三川,黄 波

(电子科技大学航空航天学院,成都611731)

1 引言

航天服最主要的功能是为航天员提供空间环境的安全防护并满足空间操作的活动性能要求[1]。与太空作业不同,月表作业要求航天服上下肢关节系统具备良好的运动性和灵活性,以便完成月面行走、标本采集、设备维修、货物运输、科学实验等舱外作业任务[2]。由于航天服真空防护所造成的余压(服内外的压差)作用,当航天员活动时,航天服不可避免地产生相当大的阻力矩[3]。该阻力矩会大大降低航天员肢体运动的灵活度、准确度和平衡性,增加体能消耗,导致航天员月表作业能力和运动能力下降。为适应未来月表低重力环境下的行走需求,航天服下肢关节系统对航天员月表作业任务的影响将作为设计与评价其助力性能的依据。

近些年,为提高航天员着服后的活动性能,常在对应的人体关节如肩、肘、髋、膝处加入软式关节结构,如平褶式[4]、波纹式[5]、橘瓣式[6]和等张立体式[7]等,但这些关节结构复杂,容易受设计和加工工艺水平影响。Sim[8]和Judnick等[9]先后提出了机械反压(Mechanical Counter Pressure)方案,但机械反压航天服对材料选取与加工工艺要求较高,且仅对人体局部起到压力防护作用。NASA和佛罗里达州人与机器认知研究机构(IHMC)共同研发X1航天员专用外骨骼服,采用联合设计模式将外骨骼与航天服材料结合,增加航天服机动与负重能力[10]。郭庆等[11]提出一种主动航天服方案,在现有航天服的基础上运用外骨骼机器人技术,并融入传感、控制、信息等技术,使其成为高度人机耦合的复杂力量随动系统。以上2种方案都未考虑航天服关节迟滞特性对控制系统的影响。

本文参考美国舱外活动单元(Extravehicular Mobile Unit,EMU)公开数据,提出一种登月助力航天服结构方案,分析航天员月表行走时下肢关节的迟滞模型和运动范围,建立航天服下肢关节阻力矩-角度映射关系,并完成二连杆动力学仿真分析,为后续随动控制提供理论基础。

2 助力航天服外骨骼结构

登月助力航天服是将力量增强型外骨骼助力技术运用到现有航天服,提高航天员的承载能力、抓举能力和运动能力,帮助航天员安全、快捷、高效、省力地实现月面行走、样本采集、装备维修等月表作业任务。其结构主要包括：①机械结构。采用仿生技术和机器学习技术,满足航天员个性化的运动要求以及舱外航天服结构尺寸,实现机械承载和助力；②感知系统。通过多传感器和数据融合技术等感知航天员的运动姿态和生理状态等,同时较为准确地预测航天员的运动意图；③控制系统。针对航天员的动作需求,通过多功能传感与柔顺耦合控制技术实现外骨骼的多自由度、助力随动的控制目标；④驱动机构。根据控制系统指令,综合应用电机、液压等驱动方式,实现轻质、大功率、小型化的驱动助力；⑤能源系统。提供各系统所需的能源。

图1为设计的登月助力航天服的下肢关节外骨骼结构图。登月助力航天服下肢关节外骨骼具体尺寸参考美国EMU舱外航天服的公开数据,仿人设计实现结构的助力和穿戴舒适,碳纤维材料的设计加工实现轻质、高强度的外骨骼；采用肌电、惯性、压力传感器等多传感器感知和数据融合技术实现步态识别在线辨识和运动意图准确预测；安装4个Maxon电机于髋、膝关节处,通过控制电机驱动实现助力。

图1 登月助力航天服下肢关节外骨骼结构Fig.1 Exoskeleton structure of lower limb joints in lunar assist space suit

3 航天服下肢关节迟滞模型

航天员在执行月表作业任务时,行走、弯腰拾物、起蹲等动作的实现都需要人体下肢各关节的参与配合。本文主要研究航天员在矢状面内的行走,其运动姿态可划分为支撑相和摆动相,主要表现为下肢关节角度的规律变化。而航天服因为充压设计(EMU,29.67 kPa)导致关节阻力矩具有极强的迟滞效应,具体表现为关节阻力矩不仅与当前关节角度有关,而且与关节历史运动角度有关,具有高度的非线性和迟滞性,使得普通函数无法实现对阻力矩的准确描述[12]。同时,登月助力航天服要实现关节助力效果,需要分析航天服下肢关节的迟滞特性,建立下肢关节阻力矩-角度映射关系,为后续的动力学方程建立提供关节阻力矩函数。

国内外对舱外航天服阻力矩的测量装置主要有NASA基于fish-scale方法的结构[13]、MIT的内置式机器人RSST(Robotic Space Suit Tester)[14]以及国内的外置式机器人[3]。本文根据fish-scale方法搭建下肢关节阻力矩测试装置,系统包括内压为29.67 kPa的充气下肢、0～50 Nm扭矩传感器、-180°～180°角度传感器、0～30 kg拉压力传感器和相应固定机械机构,其结构示意图如图2所示。

图2 下肢关节阻力矩测试装置Fig.2 Measuring device of joint resistance torque in lower limbs

选择应用最为广泛的Preisach迟滞模型[15]描述已得到的下肢关节阻力矩-角度映射关系。该模型是一种完全数学模型,利用数学公式来模拟迟滞曲线,能较精确地描述迟滞特性,其模型的数学表达式为式(1)：

式中,u(t)表示当前时刻t的输入,f(t)表示当前时刻t的输出,γα,β[u(t)]为基本迟滞变换函数,如图3所示。α和β(α≥β)分别表示向上和向下变换的变换点,μ(α,β) 是γα,β权重函数。

基本迟滞变换函数γα,β只有两个输出值+1和-1。确定方法为：当输入u(t)＞α时,γα,β[u(t)]=+1；当输入u(t)＜β时,γα,β[u(t)]=-1。 当输入介于α和β之间时,若输入u(t)沿着 ABCDE的方向逐渐上升,γα,β[u(t)]=-1；若输入沿着 EDFBA 的方向逐渐下降,γα,β[u(t)]=+1[16]。 通过γα,β的加权组合引入方向性和记忆性,从而描述迟滞现象。

图3 基本迟滞变换函数Fig.3 Elementary hysteresis transform function

其中,输出f(t)的积分区域为S,其大小为S={(α,β)α＞β,α＜αmax,β＞α0},其中α0＜αmax为积分区域α轴上的某一特定值,如图4所示。

图4 Preisach(α-β)积分区域的几何表示Fig.4 Geometrical representation of Preisach(α-β)integral area

积分区域S根据关节运动轨迹划分为S+区域(基本迟滞变换函数γα,β输出为+1)、S-区域(基本迟滞变换函数γα,β输出为-1)和边界L(t)(由坐标点(Mn,mn)构建的阶梯形),这样输出f(t)可以表示为式(2)：

定义灰色三角形区域的积分值为F(Mn,mn),那么矩形区域Qk的积分值可以表示为式(3)：

将输出f(t)的积分式(2)离散化,可得式(4)：

通过式(4),仅需要得到边界L(t)的坐标点(Mn,mn)对应的F(Mn,mn)就可以得到关节阻力矩输出值f(t)。

4 助力航天服下肢外骨骼动力学分析

4.1 登月助力航天服下肢外骨骼模型

航天员的行走平面可以分为额状面、水平面、和矢状面。在行走时,额状面和水平面的相对运动较少,为简化计算,在进行人体下肢建模时仅关注航天员行走时矢状面内髋关节、膝关节和踝关节的角度变化[17]。而人体下肢的运动是人体肌肉及韧带带动骨架绕关节的旋转运动,因此人体左右大腿、左右小腿、左右足分别可以简化为单连杆结构,各杆之间通过关节连接。但由于踝关节运动范围较小且结构复杂,故将左右足和左右小腿构成一个整体,而左右两肢受力一样,因此将下肢简化成二连杆机构,具体如图5所示。

图5 下肢二连杆简化模型Fig.5 Two-link simplified model of lower limbs

图中,i=1,2,其中qi表示关节转角,li表示连杆i的长度,lci表示上一个关节与连杆i质心的长度。

4.2 牛顿-欧拉动力学方程

相对于拉格朗日动力学建模方法,牛顿-欧拉法适用于多杆机构的动力学分析与计算机编程计算。由于登月航天服下肢模型较为简单,容易得到二连杆之间的耦合关系以及单个连杆的受力情况,可以通过使用广义坐标下qi及其导数求解连杆i的速度q·i与加速度q¨i,然后反向迭代计算连杆i的作用力fi和作用力矩τi[18]。连杆i质心加速度ac,i可以表示为式(5)：

式中,Rii-1为坐标系Oi-1到坐标系Oi的旋转矩阵,ωi为连杆i角速度,ω·i为连杆i角加速度,ri,ci为从点oi到连杆i质心oci的向量,ae,i-1为连杆i-1的末端加速度,则连杆i的末端加速度ae,i可以表示为式(6)：

式中,ri,i+1为从点oi到oi+1的向量,其代替式(5)的ri,ci。连杆i力平衡方程可以表示为式(7)：

式中,fi+1为连杆i+1施加到连杆i上的力,mi为连杆i的质量,gi为连杆i的重力加速度,可以表示为式(8)：

式中,g为地球重力加速度。连杆i力矩平衡方程可以表示为式(9)：

式中,τi为连杆i-1施加到连杆i上的力矩,Ii为连杆i的转动惯量,τH,θi为航天服 Preisach模型得到的关节阻力矩。

对于登月航天服下肢模型,已知以下向量,见式(10)、(11)：

其中,x、y、z分别表示基座坐标系的3个单位坐标向量。

通过公式(5)可得到连杆1质心加速度ac,1和连杆 2 质心加速度ac,2,见式(12),(13)：

通过公式(6)可得到连杆1末端加速度ae,1,见式(14)：

通过公式(8)可得到连杆1质心加速度g1和连杆2质心加速度g2,见式(15),(16)：

由于登月航天服下肢关节处加驱动机构,协助航天员完成月面行走任务,其驱动力大小可以通过公式(9)得到式(17)、(18)：

5 仿真分析

5.1 下肢运动参数

根据中国航天员的选拔条件,即年龄23到35岁,身高 1.6～1.72 m,体重 55～70 kg,对比国内4名男性航天员的身和体重参数[19],结合《GB/T10000-1988中国成年人人体尺寸》[20]拟定男性航天员身高为1683 mm,体重为59 kg,大腿长l1为466 mm,小腿长l2为370 mm。再参考《GB/T17245-2004成年人人体惯性参数》[21]作为人体相对质量分布、质心相对位置以及人体体段转动惯量数据来源。其中,m1表示为人体大腿质量(kg),m2表示为人体小腿质量(kg),lc1表示为人体大腿质心距离(mm),lc2表示为人体小腿质心距离(mm),Iz1表示为人体小腿长度(mm),Iz1表示为人体大腿转动惯量(kg·mm2),Iz2表示为人体小腿转动惯量(kg·mm2),具体数值如表1所示。

表1 下肢关节主要设计参数Table 1 Design parameters of the lower limb joints

利用Opensim肌肉骨骼系统软件模拟人体正常行走运动,使用系统自带的Leg6Dof9Musc模型,修改其中参数,得到航天员着服后下肢关节角度变化。 其中,髋关节为-29.75°～16.92°,膝关节为-68.69°～3.07°。根据美国EMU航天服在29.67 kPa内压情况下,关节运动范围平均减少约15%[22]。因此,拟定在1个步态周期内,髋关节运动范围为-25°～15°,膝关节运动范围为-60°～0°,如图 6 所示。

图6 下肢关节角度变化曲线Fig.6 Changes of joint angle in lower limbs

5.2 动力学仿真

根据5.1节下肢运动参数,在Creo软件中设计并装配出航天员下肢模型。由于人体骨骼十分坚硬,因此在不考虑变形的情况下,可以视人体的各部分为刚体,视关节为铰链。在自由度配置方面,设计髋关节、膝关节和踝关节各只需1个自由度,下肢模型通过保存为Parasolid格式文件导入到ADAMS中,赋予各构件的材料属性,修改质量信息,根据下肢运动参数设置模型构件的质量和转动惯性数值,在髋关节和膝关节处定义转动副,而在踝关节处定义固定副,将足部固定在小腿处。仿真模型如图7所示。

考虑到登月助力航天服中外骨骼的作用之一是帮助航天员克服航天服内余压的阻力,为了模拟出航天服对航天员运动的阻碍作用,在髋关节和膝关节处设置了2个单向力矩函数来模拟航天服的关节阻力,其中模型中的权重函数μ(α,β)需要通过试验测得,为避免构件重力的影响,在ADAMS中设置3个方向的重力加速度为0。

图7 下肢二连杆结构仿真模型Fig.7 Two-link structure simulation model of lowerlimbs

通过仿真可以得到髋关节的阻力矩随关节角度变化曲线,如图8所示,其中Flexion表示航天员髋关节从15°弯曲至-25°,Extension表示航天员髋关节从-25°伸展至15°。可以看出,髋关节前伸时阻力矩要小于后曲时的阻力矩,主要原因是因为仿真时髋关节正处于向后弯曲的动作。其中,登月助力航天服充压状态下髋关节Flexion方向力矩范围为-84.68～52.51 Nm,Extension方向力矩范围为-86.95～45.89 Nm。

膝关节的阻力矩随关节角度变化曲线如图9所示,其中Flexion表示航天员髋关节从0°弯曲至-60°,Extension表示航天员髋关节从-60°伸展至0°。可以看出,膝关节伸展时阻力矩要小于弯曲时的阻力矩,因为仿真时膝关节正处于向后弯曲的动作。其中,登月助力航天服充压状态下膝关节Flexion方向力矩范围为-14.55～13.58 Nm,Extension方向力矩范围为-13.52～13.52 Nm。由图8、9可以看出登月助力航天服下肢柔性关节在前屈/后伸时关节阻力矩相对于关节角度的变化关系呈明显的迟滞特性。同时也证明航天服柔性关节阻力矩不仅与当前关节角度有关,而且与关节运动方向有关,与文献[23]、[12]一致。

已知航天员在1个步态周期内下肢关节角度变化曲线和充压(EMU,29.67 kPa)后登月航天服关节阻力矩-角度映射关系曲线,为了实现对1个步态周期内的电机控制,需要结合上述曲线得到关节阻力矩随步态周期变化曲线,如图10所示。由于在1个步态周期内,下肢关节运动角度和运动方向实时发生变化,特别是在髋关节,运动方向变化频繁,容易导致阻力矩突变。为了实现电机的柔顺控制,对关节阻力矩-步态周期映射关系进行数据拟合处理,采用MATLAB的高阶正弦函数逼近(Sum of Sin Functions)关节阻力矩随步态周期变化曲线。其中,髋关节在1个步态周期内力矩范围为-85.54～54.15 Nm,膝关节在1个步态周期内力矩范围为-13.58～18.19 Nm。

图8 登月航天服髋关节阻力矩-角度曲线Fig.8 Damping torque-angle curves of hip joint inlunar spacesuit

图9 登月航天服膝关节阻力矩-角度曲线Fig.9 Damping torque-angle curves of knee joint in lunar spacesuit

得到登月助力航天服关节阻力矩-步态周期的拟合函数,可以在下肢髋关节、膝关节处添加2个阻力矩函数,根据登月助力航天服下肢外骨骼动力学方程求解下肢关节的驱动力矩。同时为了模拟月面低重力环境,需要修改ADAMS系统的重力加速度,其值为1/6 g。登月助力航天服下肢关节驱动力矩随步态周期变化曲线如图11所示。其中,下肢关节驱动力矩在步态周期50%左右变化较大,说明当支撑腿在将要离开地面时会造成关节驱动力矩的突变。同时髋关节在1个步态周期内驱动力矩范围为-99.65～124.43 Nm,膝关节在1个步态周期内驱动力矩范围为-27.78～32.78 Nm,对比图10,说明一个步态周期内航天服关节阻力矩在驱动力矩所占比例可达50%左右,使得航天员需要消耗更多的能量来完成相应的动作。

图10 登月航天服关节阻力矩-步态周期曲线Fig.10 Damping torque-gait cycle curves of joint in lunar spacesuit

6 结论

在未来月面正常行走过程中,航天员下肢关节的运动角度和运动方向呈现周期性变化,导致航天服关节阻力矩突变。通过迟滞模型建立与分析,髋关节在0°时,不同的运动方向导致阻力矩相差54.34 Nm,而膝关节在-30°时,不同的运动方向导致阻力矩相差8.87 Nm。通过对登月助力航天服下肢关节动力学仿真,模拟月面低重力环境,得到1个步态周期内下肢关节阻力矩和驱动力矩变化曲线,且关节阻力矩在驱动力矩所占比例可达50%左右。仿真结果为后续驱动机构的力矩控制提供理论参考。未来工作主要集中在低重力环境对航天员步态的影响以及驱动机构运动控制系统研究。

图11 下肢关节驱动力矩-步态周期曲线Fig.11 Driving torque-gait cycle curves of lower limb joint