月面探测器圆形薄膜太阳翼展开动力学建模与分析

辛鹏飞，吴跃民，荣吉利，危清清，刘 宾，刘 鑫

（1. 空间智能机器人系统技术与应用北京市重点实验室，北京 100094；2. 北京空间飞行器总体设计部，北京 1000941；3. 北京理工大学 宇航学院，北京 100081）

引 言

载人登月是人类深空探测活动中的重要一步，选择高效、可靠、大收纳比的太阳翼能够为月面着陆器及未来的月球基地提供充足的能源保证[1]，这一直以来是载人登月研究的重点之一。尤其是在月球极区环境下，光照条件差，对高效太阳翼的需求更加显著。而圆形薄膜太阳翼具有结构紧凑、质量轻、功率质量比高、低转动惯量、展开刚度高、可重复展收和扩展性好等特点[2]，非常适应用于深空探测。美国国家航空航天局（National Aeronautics and Space Administration，NASA）在新的载人登月路线图中，公布了搭载圆形薄膜太阳翼的月面探测器构型，效果如图1所示。

圆形薄膜太阳翼最早由美国AEC-able公司（现为诺格公司）研发，由移动箱板、固定箱板、展开机构、薄膜分片等组成，如图2所示，收拢体积小、展开可靠性高。到目前为止，世界上仅美国开展了该领域的深入研究，并成功发射了多组圆形薄膜太阳翼。太阳翼工程样机由三角形柔性薄膜分片组成[3]，在2007年发射的“凤凰号”（Phoenix）火星着陆器以及2018年发射的“洞察号”（Insight）火星探测器上均配套了2个直径为2.1 m的圆形薄膜太阳翼（展开效果如图3所示）并成功实现了展开应用。从2015年开始，NASA一直在为“天鹅座”-OG5飞船设计配套的圆形薄膜太阳翼。

在薄膜结构的动力学建模分析方面，许多学者分别基于流固耦合方法[4]、多刚体动力学方法[5]以及能量–动量模型方法[6]对薄膜结构开展研究，展现出在特殊情况下的充气薄膜动力学特性。Liu等[7]采用绝对节点坐标方法对薄膜结构进行建模，研究了薄膜结构的褶皱、屈曲等动力学特性，仿真结果与其他理论研究方法及相关实验结果较好吻合。2007年，NASA在圆形薄膜太阳翼实验中，利用单点激振以及激光测振仪在真空环境下对半径为1.6 m的圆形薄膜太阳翼缩比模型进行了模态测试[4]；之后又利用单点激振以及加速度传感器对采集的实验数据进行了进一步分析；2013年，NASA利用有限元软件对圆形薄膜太阳翼进行模态分析[9]，充分考虑到重力悬吊系统及空气作用的影响，得到了太阳翼前三阶整体模态。尽管美国对该型薄膜太阳翼进行了大量的理论分析和试验分析，但测试数据、研究成果发表极少。

图1 NASA展示的载人登月效果图Fig. 1 The artist’s rendering NASA released

图2 圆形薄膜太阳翼构型Fig. 2 Configuration of circular

图3 “洞察号”上的圆形薄膜太阳翼Fig. 3 Insight’s circular membrane solar array

薄膜结构在月面低重力环境下展开复杂，地面试验难以准确模拟，若采用实际任务测试实验的方法在成本上难以承受[10]；同时圆形薄膜太阳翼展开动力学试验较难操作，实验数据获取困难。因此只有精确地完成月面环境下圆形薄膜太阳翼结构的有序展开动力学特性分析，才能有效提高圆形薄膜太阳翼的展开可靠性和展开精度[11]，建模仿真分析工作具有重要的理论意义和实际工程价值。

1 展开动力学建模

1.1 数值建模基础

本文采用绝对节点坐标法搭建圆形薄膜太阳翼动力学模型，该方法具有简洁的动力学方程形式[13]，可以完成对复杂结构精确的刚柔耦合动力学分析。

一个绝对节点坐标全参数梁单元如图4所示，其单元上任一点的全局位置矢量定义为

其中：rP为任一点P的全局位置矢量；x、y和z是点P的局部坐标；S为单元形函数。

单元节点坐标e可以表示为

其中：ri，x表示r对x的偏导数，依次类推。

由式（2）可知，该单元任一节点坐标均包含一个位置矢量与3个斜率矢量，共计采用24个广义坐标来描述。

图4 基于绝对节点坐标法的全参数梁单元Fig. 4 Full parameterization beam element based on ANCF

与此类似，本文采用考虑Kirchhoff假设的绝对节点坐标矩形缩减薄板单元搭建大面积薄膜结构数值模型[14]。如图5所示，该单元类型略去了板沿厚度方向的变形，即节点广义坐标不包含沿厚度方向的梯度向量，以A点为例，该点的节点广义坐标为

图5 缩减绝对节点坐标矩形薄板单元Fig. 5 Reduced rectangular thin plate element based on ANCF

每个缩减薄板单元包含4个节点，合计36个节点自由度，具体形函数形式等可参考文献[14]中的相关内容。基于以上绝对节点坐标单元类型，本文中圆形薄膜太阳翼的支撑肋条结构采用绝对节点坐标梁单元进行建模，薄膜太阳翼结构采用缩减绝对节点坐标薄板单元进行建模，搭建完整结构的动力学数值仿真模型，如图6所示。

图6 完整结构动力学仿真模型示意图Fig. 6 Demonstration of full structure dynamics model

通过第一类拉格朗日方程，可以推导系统方程为

其中：M为系统质量矩阵；q为系统广义坐标；C为系统约束方程；Cq为系统约束方程对广义坐标的雅克比矩阵；λ为拉格朗日乘子；Q（q）为系统广义外力矩阵；F（q）为系统弹性力矩阵。本文中，采用广义alpha方法[15]求解该方程组。

1.2 接触碰撞检测算法

在太阳翼展开过程中，薄膜分片间以及薄膜与肋条之间存在复杂的接触碰撞现象。为了高效检测薄膜间的接触碰撞，将整个接触检测过程分为两步执行，即全局检测阶段和局部检测阶段[16]。在全局检测中，基于层次包围盒思想，构造包围盒树状层次结构，如图7所示，快速匹配彼此靠近的四边形接触块，判断接触可能性，形成潜在接触对；在局部检测中，通过检测点与四边形接触块之间的接触，得到点与曲面间接触的粗略碰撞信息，以此作为初值进行迭代，得到薄膜结构间的精确碰撞点。

图7 层次包围树结构示意图Fig. 7 Demonstration of hierarchical bounding volumes

图8为两个可能发生接触的绝对节点坐标薄板单元，通过全局检测，得到加粗的曲线，表示的2个四边形接触块（潜在测试对）在彼此靠近。

此时需要在绝对节点坐标单元2上找到一点Q，使得接触块a上的接触点P到Q的距离为点P到单元2的最短距离，一次判断具体接触情况。设Q在单元2上的局部坐标为（ξQ，ηQ），若点Q为单元2内一点，则有

其中：S为薄膜单元的形函数；rP为P点的全局位置矢量；e2为单元2的广义坐标。

图8 局部检测示意图Fig. 8 Demonstration of local detection

略去高阶小量后，式（5）可以快速求解单元间的最短距离，进而获得接触点及接触力信息[17]。

2 关节展开轨迹规划

2.1 关节空间轨迹规划定义

圆形薄膜太阳翼结构的展开机构仅有一个自由度，定义电机驱动轨迹函数，实际上就是定义了圆形薄膜太阳翼的展开轨迹。

通常称电机函数设计空间为关节空间。一般地，在起点到终点的运动过程中，由静止状态开始，在规定的运动时间内运动到目标位置后保持静止。

为了更好、更灵活地控制轨迹曲线，本文采用贝塞尔（Bézier）曲线作为电机驱动函数，贝塞尔曲线是能够描述复杂形状的一种曲线，现在已成为计算机图形学相当重要的参数曲线，在计算机辅助设计（Computer Aided Design，CAD）领域广泛应用[18]。选择贝塞尔曲线进行轨迹规划有以下优点：

1）轨迹能够满足运动学约束，轨迹、速度和加速度连续；

2）轨迹方便可调，且不影响上述准则。

贝塞尔曲线由数据点和控制点两部分组成，数据点控制曲线的起点和终点，控制点控制曲线轨迹的弯曲程度。一般地，数据点和控制点也统称为控制点。假设曲线为n次曲线，则控制点数量为n+1，控制点编号为Pi（i= 0，1，···，n），其中：P0为曲线起点；直线P0P1为起点切线方向；Pn为曲线终点；直线Pn-1Pn为终点切线方向。除首尾控制点外，其他控制点通常不在曲线上。

因此，定义空间n+1点的位置Pi，则n+1阶（即n次）贝塞尔曲线的描述为

其中：Bi,n（u）是n次Bernstein基函数，且满足

其中i= 0，1，2，···，n。

根据式（6）定义，期望初始及终止时刻关节速度均为零，需要贝塞尔曲线的前两个控制点重合，后两个控制点重合，以使得两个位置上的斜率均为零，即速度为零。另外还需要至少两个控制点控制曲线的弯曲程度，本文中采用的轨迹规划函数为五次贝塞尔函数（共6个控制点），即电机的驱动函数。根据上述仿真条件，得到的电机驱动转角随时间变化曲线如图9所示。

图9 贝塞尔曲线转角轨迹Fig. 9 Rotation trajectory based on Bézier curve

2.2 关节驱动展开数值仿真

电机驱动展开圆形薄膜太阳翼的数值仿真条件设置为：初始状态为时间零点、转角值为零度；在本算例中，为凸显结构展开的振动效果，运动规划时间和仿真总时间均设置较小，分别设为10 s和15 s，后5 s的时间用于考察残余振动效果；终止状态为时间结束、转角值为2π（完全展开）；转角中π出现的时刻为5 s时刻。

圆形薄膜仿真模型主要参考美国“凤凰号”（Phoenix）太阳翼测试样机数据。系统内刚体包括箱板结构、短粗连接件机构等，数量为30，柔性肋条数目为9，弹性模量为70 GPA，泊松比为0.33，每根肋条长度为1.05 m，均划分为45个绝对节点坐标梁单元；薄膜分片数量为10，每个分片薄膜共包含126个缩减绝对节点坐标薄板单元。系统总单元个数为1 695，广义坐标数量为55 440。

NASA选用了聚酰亚胺Kapton作为特种薄膜材料。该材料具有优良的化学稳定性、耐高温性、坚韧性、耐磨性、阻燃性、电绝缘性等，目前广泛应用于航空航天器领域，其基本材料参数如表1所示。

表1 薄膜材料参数Table 1 Parameters of membrane material

为了评估所采用的贝塞尔函数方法在轨迹规划效果以及残余振动抑制方面的作用，采用工程中常用的Sine轨迹函数作为电机驱动轨迹，进行展开控制对比分析。

在数值仿真分析中，分析结果考察柔性中间肋条在大范围运动中的转动轨迹误差，轨迹误差φ定义如图10所示，为肋条实际构型末端点切线与预定构型的夹角。

图10 误差转角定义Fig. 10 Definition of rotation angle error

图11为柔性中间肋条末端偏离预定轨迹的误差曲线，图12为单独提取的10 ～ 15 s残余振动误差曲线。从图11～12中观察得到以下结论：所提出的五次贝塞尔函数转动规划轨迹在初始阶段和终止阶段均较为平缓，对电机驱动的要求低，适合工程应用；在轨迹运动过程中，贝塞尔曲线控制下的肋条振动较大；但在残余振动过程中，Sine轨迹函数控制的残余振动关节角误差幅值最大为0.026 rad，而贝塞尔轨迹误差幅值最大为0.006 rad，误差比例为1∶0.23，即在残余振动过程中，贝塞尔曲线规划函数能有效地降低残余振动，控制轨迹误差。

图11 中间肋条轨迹误差曲线Fig. 11 Curves of trajectory error of middle ribs

图12 中间肋条残余振动误差对比Fig. 12 Comparison of residual vibration error of middle ribs

2.3 贝塞尔曲线优化设计

贝塞尔曲线还可以通过局部调整控制点的位置来调整曲线局部的特征。经不断调整，在仅调节中间两个控制点的情况下，图13为控制点调整后的贝塞尔曲线图。可以看出在起始和终止时刻，曲线更加平缓，这意味着这两段时刻，电机转动角速度更小，对工程应用更加友好。

图14为关节角误差结果显示曲线，可以得到以下结论：调整后的贝塞尔曲线在运动过程中的轨迹跟踪误差方面表现不佳，最大关节角误差由0.098 rad增加到0.107 rad，变化9.18%。图15为残余振动阶段关节角误差结果显示曲线，调整后的贝塞尔曲线在残余振动阶段抑制残余振动方面有所提升，最大关节角振动幅值由0.006 rad进一步降低到0.005 3 rad，变化11.67%。结果从幅值来看降低幅度不大，效果并不明显；调整贝塞尔曲线的控制点对残余振动抑制方面影响不大，通常还伴随着运动过程中关节角误差增大。

图13 贝塞尔曲线调整Fig. 13 Adjustment of Bezier curve

图14 中间肋条轨迹误差曲线Fig. 14 Curves of trajectory error of middle ribs

图15 中间肋条残余振动误差对比Fig. 15 Comparison of residual vibration error of middle ribs

单独采用五次贝塞尔曲线规划函数驱动圆形薄膜太阳翼结构实际展开，数值仿真结果表明该轨迹规划函数能够提高柔性肋条末端轨迹的运动精度，降低柔性振动幅值，使得结构展开过程更为可控和精准。但由于函数本身的特点以及结构柔性的存在，运动过程中实际轨迹与规划轨迹之间的误差较大。单独优化贝塞尔曲线本身对运动过程跟踪和残余振动抑制作用不大。

3 前馈–反馈联合展开控制

3.1 前馈–反馈联合展开控制算法

轨迹跟踪控制研究一般需要轨迹规划方法与控制算法相结合进行处理[19]。而在现有的轨迹跟踪控制理论中，应用最广泛的两类控制方法分别是前馈控制和反馈控制，这两种控制方法都有一定的优点，同时也伴随着一定的缺点。

在此基础上，本文采用前馈–反馈联合控制方法对圆形薄膜太阳翼结构进行展开控制研究。针对基于绝对节点坐标方法建立的数值仿真模型，首先求解得到驱动约束对应的约束反力/力矩，即为所需要的前馈力/力矩。而反馈控制采用工程中常用的PID（Proportion Integral Derivative）控制方法，易于实现，且对系统模型精确度要求不高。

文献[20]中推导的关节转动副前馈力/力矩公式形式简单，计算方便，不仅适用于转动关节，也适用于其他驱动关节的求解。关节转动轴ξ方向即为力矩方向，力矩大小为

式（8）的推导方法也能推广至其他约束形式对应的前馈力/力矩的求解。

为了得到光滑的前馈力/力矩从而避免引起系统过大的振动，Liu等[20]进行了深入研究，结果表明，高弹性模量的模型本质上更接近于多刚体模型，以高频振动为主，因此建议采用多刚体模型的前馈力/力矩。本文所研究的圆形薄膜太阳翼结构柔性肋条及箱板均属于高弹性模量材料，采用上述结论，利用刚体模型进行逆动力学分析，从而得到前馈力/力矩。反馈控制方法采用工程上容易实现、效果良好的PID控制方法。前馈–反馈联合控制方法的方框图如图16所示。

图16 前馈–反馈联合控制方法流程Fig. 16 Flow chart of forward-feedback joint control

基于以上所述的前馈力/力矩和PID控制策略，前馈–反馈联合控制方法在驱动关节处施加的力/力矩大小为

其中：Tf为计算的前馈力/力矩；θp和θ分别是规划的转角轨迹和实际的转角轨迹；Kp、Ki、Kd分别为PID反馈增益。

3.2 前馈–反馈联合展开控制算例

圆形薄膜太阳翼结构几何参数和材料参数同样如前文所示，关节转动副的结构参数与几何参数如表2所示，所述的关节均为单自由度、绕转轴转动的关节，齿轮副中的固定齿轮和关节的轴瓦固结，轴瓦、固定齿轮与固定箱板均固结于初始位置。

表2 关节转动副几何及材料参数Table 2 Geometric and material parameters of joint

关节转动副的转角规划函数如公式相同。Kp、Ki、Kd分别取值为0.465、16.3、4.03。为了对比分析前馈–反馈联合控制方法的控制精度问题，另外分别采用仅前馈控制、仅PID反馈控制对圆形薄膜太阳翼结构展开过程进行控制对比研究。

基于搭建的动力学仿真数值模型，进行月面重力环境下的圆形薄膜展开过程动力学分析。求解得到圆形薄膜展开过程构型如所示，结果表明，采用的展开控制策略能够驱动结构有序展开。展开过程构型如图17所示。

图17 月球重力工况展开过程示意图Fig. 17 Deployment under micro-gravity atmosphere

图18～19结果曲线显示：①3种控制策略都能有效跟踪关节转动的预定轨迹，但控制精度有所差别，联合控制策略控制精度在展开开始时误差较大，但之后迅速降低，精度高、衰减快；②仅采用前馈控制的情况下，实际轨迹存在振荡情况，并没有稳定，误差与不施加任何控制方式的结果基本相同；③采用PID反馈控制需要多次试验才能得到合适的反馈增益，展开控制精度相对较高；④而前馈–反馈联合控制方法对增益系数并不敏感，少量试验后即可得到较好的控制效果；⑤在残余振动阶段，采用前馈–反馈联合控制方法得到的控制误差几乎为零，显示了该方法在残余振动阶段具有显著的抑振作用。

图18 轨迹跟踪关节角误差Fig. 18 Curves of trajectory tracking errors

图19 残余振动关节角误差对比Fig. 19 Comparison of residual vibration error

太阳翼驱动展开过程中，面向展开过程的控制精度和残余振动的抑制能力等任务需求，可以采用前馈–反馈联合控制方法配合电机轨迹规划函数，能够经过少量仿真试验得到较好的控制参数，展开过程和残余振动阶段都获得较好的控制精度。

获得的电机驱动力矩如图20所示，随着结构展开，电机的驱动力矩也逐渐增大，最大值达到287 N·m。完全展开太阳翼结构后，驱动力矩迅速降低，直至为零。

图20 电机驱动力矩时间历程曲线Fig. 20 Times history of motor drive torque

4 结 论

对于圆形薄膜太阳翼结力学仿真分析困难、薄膜展开过程复杂的问题，本文利用绝对节点坐标方法搭建了NASA样机尺寸的圆形薄膜太阳翼动力学数值分析模型，基于贝塞尔曲线对关节空间轨迹进行规划，并结合前馈–反馈联合控制策略进行了结构展开动力学数值仿真分析。仿真结果表明：该展开驱动策略能够高精度、有序、稳定地展开圆形薄膜太阳翼结构，并对残余振动具有显著地抑制作用；展开驱动力矩最大值为287 N·m；薄膜运动复杂，张紧–回弹现象明显，最终随结构展开而趋于稳定。