基于Abaqus 的航天服波纹式肘关节有限元仿真平台开发

董长林， 尚 坤， 周仕明， 李道奎*

（1.国防科技大学空天科学学院， 长沙 410073； 2.空天任务智能规划与仿真湖南省重点实验室， 长沙 410073；3.中国航天员科研训练中心人因工程重点实验室， 北京 100094）

1 引言

舱外航天服主要用于航天员舱外作业时维持其所必须的生存环境，但航天服充气后，由于内外压差、运动过程中体积变化等因素的影响，关节处会产生较大的阻力矩，从而降低航天员的操作能力。 因此关节阻力矩大小是评价航天服性能的主要指标，确定和预测关节阻力矩在航天服设计与应用中具有重要意义。

目前，研究人员主要通过实验、理论、数值模拟等方法确定关节阻力矩。 Vykukal 等介绍了一种航天服关节阻力矩的测量装置，多应用于外置法测量。 Matty提出了一种水平面内的阻力矩测量方法，能够消除重力的影响，并在工程中得到广泛应用。 基于该方法，Meyen 等针对航天服软关节进行了内置法实验，得到了航天服关节阻力矩。 Schmidt将航天服分别简化为膜模型和梁模型，并进行了实验验证，结果表明航天服力学特性更符合膜模型的特点，即关节阻力矩特性主要由压缩气体导致。 刘文樵等提出了一种利用龙门式模台滑组的外置法实验方法，保证了测量过程中的关节匀速转动及转角测量的准确性。 王晓东等开展了航天员舱外作业上肢关节实验，发现关节阻力矩不仅具有迟滞特性，而且具有与运动历程相关的特性。 上述内置法及外置法实验方法，NASA及中国航天员科研训练中心均有采用。

由于实验方法较仿真成本高，随着计算机科学和仿真技术的发展，数值模拟方法被应用于求解航天服关节阻力矩。 Furuya 等分别针对充气状态的圆筒关节模型和波纹管关节模型进行了关节弯曲仿真，并分析了航天服关节弯曲过程中褶皱出现的原因以及褶皱的发展。 尚坤等针对平褶式航天服关节进行了有限元仿真分析，并通过能量法分析发现关节阻力矩是由于气体压缩、材料变形及材料摩擦共同引起。 王鲁豫等针对波纹式髋关节单历程关节阻力矩特性进行了仿真分析和优化，降低了关节阻力矩。 在理论方面，Schmidt建立了航天服关节变形的物理模型，对航天服阻力矩进行了计算，但由于体积计算存在误差且只考虑体积变化引起的阻力矩，导致阻力矩计算仅适用于关节在小范围运动的情况。李广利等基于阿布拉莫夫对软关节内外平衡关系的假设，利用能量法推导了线弹性范围内软关节阻力矩方程，但并未考虑材料厚度变化、材料非线性以及限制带等实际情况，因此只适用于小转动。

由于航天服关节的实际运动复杂，理论方法适用范围较小，难以直接应用；实验方法成本较高，难以用于优化设计；有限元仿真计算耗时较长。 因此大量学者转向对航天服关节复杂运动历程的阻力矩预测方法进行研究。 王昊等、赵京东等通过外置法测量得到了航天服关节阻力矩，并通过样条曲线对实验数据进行了拟合，但外插精度有限。 张新军等采用了一种物理对称的Jiles⁃Atherton 迟滞模型对阻力矩数据进行了拟合，由于模型本身的对称性，因此对阻力矩的拟合结果局部存在较大误差，在航天服阻力矩领 域 应 用 存 在 局 限 性。 Schmidt首 次 将Preisach 模型用于预测关节阻力矩特性，之后的大部分研究均是基于该模型展开。 念龙生等从数学的角度对Preisach 模型进行了改进，在此基础上，刘文樵等引入了新的内插方法并将其运用到了航天服关节阻力矩的预测上。 文剑将BP 神经网络与Preisach 模型相结合，进行了阻力矩预测；王晓东等采用NNOPM（RBF 神经网络优化的Preisach 模型）对关节阻力矩进行了预测。 但采用Preisach 模型时，需要较多实验数据才能保证预测的精度，数据较少则会导致预测出现局部误差；若直接进行参数辨识则需要长时间的计算才能获得较高精度的预测模型。

航天服人-服耦合特性下模型复杂，参数化难度大，复杂运动仿真分析耗时长。 因此有必要开发一套专用于航天服关节的有限元仿真平台，并将关节阻力矩预测方法集成到软件内部以实现快速预测。 人-服耦合运动的求解是一个强非线性问题，Abaqus 较其他有限元软件具有明显优势，同时为提高开发效率并保障可靠性，现有研究多基于Abaqus 进行二次开发。 目前针对Abaqus进行二次开发的方法主要有3 种：①通过Python进行参数化建模与分析，但该方法可视化程度低，要求使用者具备一定的专业基础和编程能力；②基于Abaqus 所提供的RSG（Really Simple GUI）进行简单的插件开发，但各插件之间的关联性和系统的封装性均较差，且该方法继承性也较差；③基于C++、VS、MFC 等进行图形界面开发，再通过命令流对Abaqus 进行控制，由于利用外部软件进行开发，交互性差，功能有限，输出结果形式固定。 因此有必要开发可视化程度高、可操作性强、封装型好、集成性强的航天服关节有限元仿真平台。

航天服肘关节是完成上肢作业的主要关节，且一般采用波纹式结构来降低肘关节运动时产生的阻力矩。 本文以航天服波纹式肘关节为例，针对上述问题，考虑体积与压力变化、材料非线性和人-服耦合特性，建立参数化的仿真模型；并对Preisach 模型进行改进，提出复杂运动历程中的关节阻力矩预测方法；基于参数化仿真模型与阻力矩快速预测方法，开发一套航天服波纹式肘关节有限元仿真平台。

2 仿真平台设计

2.1 平台整体框架

如图1 所示，航天服波纹式肘关节仿真平台由仿真分析与阻力矩快速预测2 个功能模块和1个GUI 图形交互界面组成。 仿真分析模块的主要功能为：对航天服肘关节进行建模和计算，得到肘关节的阻力矩-转角迟滞特性曲线。 阻力矩快速预测模型主要功能为：根据仿真分析得到迟滞特性曲线，对复杂运动下关节阻力矩的快速预测。GUI 图形交互界面的主要功能为：将上述2 个功能模块进行集成，并实现各模块之间的交互和调用。

图1 航天服波纹式肘关节仿真平台模块功能执行流程图Fig.1 Module function and execution flow figure of spacesuit corrugated elbow joint simulation platform

2.2 仿真分析模块

人-服耦合仿真分析流程主要包括以下步骤：①在CAD 软件中建立航天服波纹式肘关节的几何模型；②将几何模型（含人体）导入到CAE软件中，并根据各部分的结构特点进行网格划分；③建立航天服肘关节两端的端盖以形成封闭腔体，进而定义边界条件和载荷（关节转角）；④通过修改Inp 文件的方式输入实验数据来定义材料属性；⑤进行显式动力学分析；⑥结果后处理。

2.2.1 几何建模

舱外航天服分为3 层结构：外部的防护层、中间的限制层以及内部的气密层。 其中，限制层为织物材料，是主要的承载结构，因此以限制层作为主要研究对象。

本文通过SolidWorks 软件建立了如图2 所示的波纹式肘关节模型，包括波纹管结构和限制带结构。 同时，根据GJB 36A-2008《飞行员人体模板设计和使用要求》的中号模板，建立了170 cm人体的大臂和小臂关节模型。

图2 波纹式肘关节和人体手臂模型Fig.2 Corrugated elbow joint and human arm model

2.2.2 有限元建模

根据Schmidt 的研究，航天服波纹管厚度较薄，性能更接近于薄膜结构；而限制带较厚，更接近于壳体结构。 因此，对限制带和波纹管分别进行网格划分，如图3 所示。 因人体刚度远大于织物材料刚度，故可将人体简化为刚体模型。 两端端盖与波纹式航天服形成封闭腔体，由于不考虑端盖所引起的体积变化，因此将两端端盖也设置为刚体。

图3 航天服波纹式肘关节的限制带和波纹管网格模型Fig.3 Corrugated elbow joint restraint band and corrugated pipe for spacesuit mesh model

在之前的实验研究中，获得了织物材料的力学性能。 由于织物材料的屈曲、非线性特性以及经线与纬线之间的摩擦，织物材料拉伸和剪切力学性能呈现出明显非线性，加载卸载曲线不一致，即加载卸载过程存在能量损耗，如图4 所示。因此，在定义材料属性时难以采用Abaqus 软件中自带的材料模型。 为解决该问题，一般有2 种思路：一是建立本构模型，通过VUMAT 来定义材料特性；二是修改Inp 文件，通过输入实验数据来定义材料特性。 为降低对工程人员在使用过程中的技术要求，本文采用后者。

图4 织物材料拉伸、剪切力学性能曲线图Fig.4 Graph of tensile and shear mechanical proper⁃ties of fabric materials

采用内置法的人体驱动方式，以模拟人体驱动下的充压航天服波纹式肘关节弯曲情况，即人体着服后，大臂及上端盖固定约束，小臂绕肘关节旋转一定的角度以驱动航天服肘关节弯曲。 肘关节的转动中心如图5 所示。 对小臂进行加载时，先加载至目标角度再卸载到初始位置。 航天服内部通过流体腔法充压以模拟航天服内气体环境，航天服内部压力曲线和肘关节转角曲线如图6 所示。 同时，定义人服接触以及航天服自接触。

图5 航天服装配关系及关节旋转中心图Fig.5 Diagram of spacesuit assembly relation and joint rotation center

图6 航天服内压与关节驱动加载图Fig.6 Diagram of spacesuit internal pressure and joint actuation

2.2.3 分析与求解

航天服肘关节运动仿真分为2 步：第一步为充气，对航天服充气以模拟航天服所需的气体环境；第二步为运动，小臂绕肘关节中心转动。 由于本文采用流体腔法来模拟航天服的内部气体环境，因此进行显式求解，以获得航天服关节阻力矩-转角迟滞特性曲线。

2.3 阻力矩预测模块

以仿真获得的阻力矩-转角迟滞特性曲线作为输入，采用Preisach 模型对复杂运动的阻力矩进行预测。 Preisach 模型积分定义如式（1）所示：

式中，（） 为关节转角，（，） 为迟滞函数，为符号函数，在正、负积分区域的取值分别为+1和-1。 根据的正负性可将式（1）整理为式（2）：

定义整个积分区域为，对式（2）进行整理，并将阻力矩转为正积分或负积分区域的积分形式，有式（3）：

当关节运动时，关节转角的变化将导致积分平面发生变化，如图7 所示。 积分平面正负积分区域边界拐点记录加载卸载历程中的极值点，即通过该方法对阻力矩进行预测时能够考虑运动历程的影响，如图8 所示为复杂运动历程下正积分区域分解图，能够用于计算积分值。

图7 复杂运动对应积分平面变化图Fig.7 Diagram of complex motion corresponds to the integral plane change

考虑到阻力矩可以通过积分或微分方式进行描述。 微分描述迟滞函数的方法受实验数据噪声影响较大，因此利用积分几何意义，采用权重系数法对预测点位置的积分值进行预测。

当输入数据较少时，预测点的临近数据点的权重会发生变化，进而使得预测值产生较大的局部误差。 利用图8 中积分平面边界＝上的点在正积分区域积分为0 的特性，对发生了权重偏移的点进行修正，以减少数据离散所引起的局部误差。

图8 复杂运动积分平面分解图Fig.8 Diagram of complex motion integral plane de⁃composition

2.4 GUI 图形交互界面

GUI 图形交互界面如图9 所示。 仿真平台主界面主要包括主菜单、可拖拽工具栏、功能菜单、建模与分析工具栏、辅助工具箱、子界面、显示框、信息提示框和命令行。 主菜单提供模型文件保存、视图设定、帮助文档查询等操作；可拖拽工具栏提供与建模相关的辅助功能；功能菜单可实现各功能之间转换；建模与分析工具栏包括航天服波纹式肘关节的快速建模、分析及阻力矩预测等按钮；辅助工具箱与可拖拽工具栏的功能类似，主要提供模型分割、载荷设置等辅助功能；子界面提供参数定义、模型可视化等功能，或进行错误提示；显示框进行建模过程及结果的显示；信息提示框将执行结果及辅助功能操作结果进行显示；命令行提供命令流的输入功能。

图9 航天服仿真平台及各功能模块图Fig.9 Diagram of spacesuit simulation⁃platform and function modules

图形界面允许工程人员随时介入分析流程，增强了可操作性，调用时设置了报错和中断机制，避免模型参数不匹配或模型干涉等问题导致有限元建模过程出现错误。

2.5 内部执行流程

仿真平台内部执行流程如图10 所示。 首先通过bat 文件调用航天服波纹式肘关节阻力矩仿真平台，执行App.py 脚本，创建应用程序并打开主窗口；其次，运行MainWindow.py 进行主窗口上层布局初始化，界面如图9 所示；然后，执行Mod⁃uleGui.py 和ToolboxButtonsGui.py 设置菜单栏以及功能键并关联执行事件的py 文件；接着，通过Form.py 脚本和DB.py 脚本启动核文件Kernel.py所需参数的输入界面并进行参数传递，同时设置参数报错机制，参数验证无误后将参数传递至VB文件并调用SolidWorks 建立参数化模型，导入软件并进行前处理和仿真分析，计算完成后输出阻力矩-转角迟滞特性曲线并进行显示；最后，通过GUI 图形交互界面输入关节转角变化曲线，读取阻力矩-转角迟滞特性曲线，调用核文件Kernel.py 内置的阻力矩预测方法，实现复杂运动下航天服关节阻力矩的预测。

图10 仿真平台内部执行流程图Fig.10 Internal overall frame diagram of the simula⁃tion platform

3 仿真平台使用流程

根据上述设计思路，完成了仿真平台的搭建。结合各功能模块，对仿真平台使用流程进行说明。

双击执行文件打开仿真平台，仿真平台的主要功能流程主要分为以下6 个步骤：①点击“Step1”按钮，根据模型结构提示设置参数，完成参数化建模；②点击“Step2”按钮，进行材料设置，默认的航天服材料为织物材料，也可下拉菜单选择“新材料”的选项，实现新材料的输入；③点击“Step3”按钮，根据工况类型设置载荷及边界条件，接触无需手动定义，仿真平台将自动定义；④点击“Step4”按钮，设置工作路径并保存模型；⑤点击“Step5”按钮，进行显式计算，计算完成后将自动进行后处理，并显示结果；⑥点击“Step6”按钮，输入关节运动历程曲线，计算该历程下关节阻力矩变化情况，计算完成后将自动显示并保存结果。 各步骤的交互界面如图11 所示。

图11 建模与分析参数输入界面Fig.11 Input interfaces of modeling and analyzing parameters

4 仿真平台分析算例

4.1 模型建立

利用本文开发的航天服波纹式肘关节有限元仿真平台，按照使用流程对肘关节的运动过程进行计算。 航天服肘关节的结构参数为：中间上部织物材料高度＝13 mm，中间层分隔单元高度＝6 mm，中间下部织物材料高度＝13 mm，上部织物材料与水平面夹角＝5.5°，下部织物材料与水平面夹角＝5.5°，上外部织物材料截面圆角半径＝12 mm，下外部织物材料截面圆角半径＝12 mm，内部织物材料截面圆角半径＝12 mm 以及肘关节关节半径＝82 mm。 生成的单个波纹管模型如图12 所示。

图12 单个波纹管模型Fig.12 Corrugated elbow unit of spacesuit

4.2 结果与讨论

读取仿真结果（odb 文件），得到航天服波纹式肘关节的应力云图，如图13 所示。 航天服肘关节在运动过程中阻力矩具有迟滞特性，其关节阻力矩-转角迟滞特性曲线（含位移云图）如图14 所示。

图13 航天服波纹式肘关节应力云图Fig.13 Stress cloud of spacesuit corrugated elbow

图14 关节阻力矩-转角迟滞特性曲线Fig.14 Hysteresis characteristic curve of joint resist⁃ance torque rotation

基于仿真平台对文献［7］中关节运动阻力矩进行预测，得到图15 所示航天服关节阻力矩预测结果，与文献中实验数据对比，方差不超过0.04（Nm），验证了阻力矩预测方法的有效性和正确性。

图15 关节阻力矩预测模型验证结果Fig.15 Validation results of a prediction model for joint resistance torque

同时对任意给定的关节转角在时间历程下的变化曲线进行预测，如图16 所示。 图中（a）段人服刚开始接触，阻力矩较小；（c）段关节转角未发生变化，阻力矩基本保持不变。 因此，主要针对（b）段展开研究，该段曲线有加载后卸载段、卸载后加载段以及卸载后过加载段，因此图中运动历程能够充分体现关节运动过程中的阻力矩迟滞特性以及运动历程的相关性。 采用有限元方法计算结果对阻力矩预测模型进行验证。 阻力矩预测曲线与仿真结果对比如图17 所示，方差为0.4220（Nm），从图中还可看到，仿真结果在6.3 s 处存在局部波动，此时误差最大，为6.13%；其余各处最大误差为4.83%，误差较小，验证了关节阻力矩预测模型的有效性和准确性。

图16 任意关节转角时间历程曲线Fig.16 Time history curve of arbitrary joint rotation

图17 关节阻力矩时间历程曲线Fig.17 Time history curve of joint resistance torque

5 结论

本文开发了航天服波纹式肘关节仿真平台，针对有人状态下的航天服肘关节进行了运动分析，并对复杂运动阻力矩进行了预测，得到以下结论：

1） 仿真分析模块能够实现波纹式肘关节快速建模与分析，并自动输出仿真结果。 显著减少了设计优化人员的建模及分析工作量，降低了对设计人员所需掌握专业知识的要求。

2） 阻力矩预测模块能够降低数据点减少引起的局部误差，通过实验结果和有限元计算结果验证了改进后Preisach 模型的精度。 对复杂运动历程关节阻力矩预测方法的研究有一定意义。

3） 仿真平台实现了对外部软件的调用和先进算法的内部集成，具有可视化程度高、操作性强、封装性好、集成性强的特点，可用于提高航天服肘关节结构设计的高效性和便捷性。