多刚体系统分离策略及释放动力学研究1)

罗操群 孙加亮 文 浩 胡海岩,† 金栋平,2)

*(南京航空航天大学机械结构力学及控制国家重点实验室,南京 210016)

†(北京理工大学宇航学院飞行器动力学与控制教育部重点实验室,北京 100081)

引言

多刚体系统是指由多个独立刚体通过连接约束组合而成的多体系统,具有自由度高、约束条件复杂等特征[1-6].随着航天技术发展及航天任务的多样化,多刚体系统在到达目标轨道后,需再次分离释放成多个子系统以满足任务需求,这些任务涉及通讯、监控、气象和导航等应用方向[7],如卫星互联网星座任务[8]、雷达卫星星座任务[9]等.此外,多刚体系统还广泛用于航天器在轨组装、卫星编队及空间碎片捕获等领域[10-14].多刚体系统成功分离和释放关键在于保证分离过程中刚体之间无碰撞发生,以使刚体按预定轨迹安全进入预定轨道.然而,在轨多刚体系统受科氏力及重力梯度影响,给分离过程中非线性运动及约束变化为无碰撞分离带来困难,影响释放安全性[15].

近年来,多刚体系统分离释放动力学受到普遍关注.Kılıc¸ 等[16]研究了50 颗立方星的分离释放问题,分析了不同分离频次、分离方向和分离次序下卫星动力学行为,发现所有卫星同时释放具有很高的碰撞风险.Handschuh 等[17]考虑星座卫星分离的短期碰撞风险,对同时释放时的分离速度和方向进行了优化设计,通过星体间距离最大化来避免碰撞.Jeyakumar 等[18]对卫星的分离动力学进行了分析,采用统计方法研究了动力学参数的变化,并且基于卫星分离相对速度及距离等设计了分离系统.Liu 等[19]针对多种纳卫星释放时的安全分离问题,提出了一种基于相对轨道离心率及倾角向量的分离方法.Zhang等[20]研究了小卫星弹射分离过程中含接触、碰撞和摩擦的非光滑动力学,提出了无干扰弹射分离的参数优化方法.Bridges 等[21]研究了分布在火箭上级的星群间无碰撞分离释放问题,基于释放速度的方向,提出了两种分离方案.Wermuth 等[22]设计了铍卫星从小卫星中分离的一种分段式释放策略,通过蒙特卡罗方法验证了该策略的有效性.商显扬等[23]针对一箭多星分离问题,提出了一种并联斜置发射结构布局及三星同时分离的方案,提高了分离可靠性.

上述研究大多关注多个刚体逐个或同时从运载航天器中分离释放的动力学控制,仅考虑刚体与航天器之间的运动学,未关注从航天器中脱离入轨后,进一步自主分离释放的多刚体系统非线性动力学,以及无碰撞分离问题.本文提出了一种由多个刚体紧密连接而成的在轨多刚体系统,其分离释放可通过刚体间的弹射机构自主完成,以有效提高运载航天器的空间利用率,简化分离释放操作、降低碰撞风险.为研究多刚体系统分离释放动力学,本文采用自然坐标方法建立单刚体动力学模型,继而连接而成为多刚体约束系统,分析适用的连接约束条件.同时计入轨道运动、分离时刚体间相互作用,利用拉格朗日乘子法获得多刚体系统分离-释放动力学方程.基于分离方向和分离顺序设计两种分离释放方案,规划分离时序,实现了多刚体系统的无碰撞分离.

1 动力学建模

图1 系统运动的参考坐标系Fig.1 The reference frames of system motions

多刚体系统M处于近地赤道圆轨道,建立惯性坐标系和轨道坐标系O1XYZ和O2xyz,它们均为右手正交坐标系,如图1 所示.惯性坐标系以地球中心O1为原点,O1X轴在赤道平面内并指向春分点,O1Z轴与地球自转轴重合;轨道坐标系以系统质心O2为原点,O2x轴和O2y轴分别指向天顶和系统飞行方向.

将单个刚体及系统均视为长方体构型,多刚体系统由这N个相同的独立刚体紧密连接而成.如图2 所示,设多刚体系统共有n层,每层刚体个数为N/n.单个刚体质量为m,质心惯量矩阵为J,尺寸为l×b×h,其中l,b及h分别表示刚体长、宽和高.

图2 多刚体系统构型Fig.2 The configuration of MRB system

1.1 空间刚体力学模型

考虑到多刚体系统约束复杂,这里采用自然坐标方法进行刚体系统动力学建模[24-25].自然坐标方法较常规方法的优势：一是可以避免诸如欧拉角描述时的奇异性;二是大大降低施加约束的难度.在常规方法中,对于距离等简单约束需要建立高次的约束方程且约束方程求解雅可比时需要进行多次变换.在自然坐标方法中,仅需线性方程即可描述连接约束,并且连接约束的雅可比均为常数矩阵,利于约束分析及力学建模.需要注意的是,为便于后续分析,本节建模过程中未考虑刚体轨道运动,仅考虑刚体相对轨道坐标系的运动.

自然坐标方法采用刚体上两个固定点的位矢及两个不共面的单位矢量作为广义坐标,以描述刚体的空间运动.两个固定点分别为刚体质心C和刚体上表面一点j,如图3 所示.此外,刚体上有一固连坐标系Cξηζ.初始时刻,坐标轴Cξ,Cη 和Cζ 分别与轨道系的3 个主轴平行.

图3 自然坐标描述刚体运动Fig.3 A rigid body described by NCF method

刚体广义坐标在轨道坐标系中给出,表示为

式中,rC和rj分别为固定点C和j的坐标向量,u和v为两个不共面单位向量.这里注意,本文将u和v选为坐标轴Cη 和Cζ 的单位向量且固定点j位于Cξ轴上.

这样,刚体上任意一点P在坐标系O2xyz中的位置坐标可以表示为

式中,I3为3×3 的单位矩阵.a1,a2和a3表示为局部坐标系下的C和j点位置及单位向量u和v的函数,即

式中,上横线指矢量在局部坐标系Cξηζ 下的坐标表示.对于刚体上固定的一点P而言,A矩阵为常数,故刚体上任一点的速度及加速度为

根据虚功原理[23],惯性力和外力在刚体虚位移上的虚功为零,即

式中,δWF和δWI分别为作用在刚体上的外力和惯性力所产生的虚功,它们是

式中,F为作用在刚体上的广义外力,Ms为自然坐标方法描述下刚体的质量阵.

通过式(6),可给出刚体质量阵为

单个空间刚体有6 个自由度,自然坐标方法采用12 个广义坐标,意味着这些广义坐标并非相互独立,故存在6 个约束,称之为刚体自然坐标固有约束方程.以单个刚体为例,设广义坐标为q1-q12,则固有约束可以表示为

式中,L表示点C和j之间的距离.上述6 个约束方程中,前3 个代表刚体局部坐标系主轴单位向量的数值约束,后3 个代表局部坐标系3 个主轴互相垂直的方向约束.

将式(8)记为Φs,利用拉格朗日乘子法,未计入轨道运动的单个刚体动力学方程为

式中,Φs,q表示约束方程相对于q的雅可比矩阵,λs为拉氏乘子向量.

1.2 系统约束分析

根据式(8)可知,单个刚体固有约束方程为6 个,故N个刚体组装成的多刚体系统,其固有约束方程为6×N个.

图4 中刚体1 和2 通过连接约束组装成一个整体系统,刚体灰色面代表刚体正面.连接时,两个刚体背面相互贴合、无间隙.两个刚体质心分别为C1和C2,连接约束点A为两个刚体所共有.A点坐标可以用刚体1 和2 的质心坐标表示为

式中,rC1和rC2表示刚体1 和2 质心位矢,rC1A和rC2A分别表示刚体1 和2 质心指向约束点A的位矢.另外,刚体1 和2 的还存在方向约束

将式(10)和式(11)组成的刚体间连接约束投影到各个坐标轴上,即得9 个约束方程.通过施加这些连接约束,刚体1 和2 组装成一个多刚体系统.若该刚体系统需要与新增加的刚体连接,则需再确定一个连接约束点,并建立与相连接刚体的局部坐标系间的关系,即每多一个刚体连接,则多出9 个约束方程.一层中N/n个刚体间的连接约束为9×(N/n-1)个,n层共9×(N-n)个约束方程.

图4 两个刚体之间的连接约束Fig.4 The connecting constraints between two rigid bodies

对于层与层之间的连接约束,每层内的刚体已经通过施加约束连接成一个整体,故视为一个刚体.因此,层与层之间的连接只需要建立9 个约束方程即可.对于n层多刚体系统,层间连接约束共有9×(n-1)个.

基于上述分析,由N个独立刚体组合连接而成的n层多刚体系统,总约束方程Φt数目为

注意的是,进行分离释放时,需要根据分离构型减少连接约束,从而更新约束方程.

1.3 系统动力学及解算方法

O

1

XYZ

式中,θ 为O2x与O1X轴之间的夹角,,ω 表示轨道角速度,如图5 所示.

图5 惯性系和轨道系下的广义坐标Fig.5 The generalized coordinates between inertia and orbital frames

根据式(13),在惯性系中,系统广义速度及广义加速度为

未计入约束时,系统在惯性系中的动力学方程表示为

式中,M为系统广义质量阵,FT为系统受到的分离弹射力,Fg表示系统万有弓力[27].

将式(16)代入式(17)后,并左乘CT得到轨道系下的系统动力学方程

式中,左边第一项为广义惯性力,第二项为广义科氏力,第三和第四项为广义离心力,且有

基于上述分析,计入系统约束并利用拉格朗日乘子法,即可得到多刚体系统轨道动力学方程

式中,Φr表示系统总约束方程,λt为拉氏乘子向量,为系统广义力,即

下面采用基于Newmark 方法[28-29]发展而来的广义α 方法[30-31]来求解约束微分方程(20)的高效求解.广义α 方法将指标为3 的方程(20)经过差分直接离散成代数方程进行求解,其迭代过程如下

式中,qr,n和qr,n+1为离散后的第n和第n+1 迭代步的系统广义坐标向量,并满足如下关系

式中,h为迭代步长,矢量参数a为新弓入的算法辅助参数列阵,由下式确定

式(24)中各参数的选取方法如下[32]

对于时间步tn+1,为了迭代求解式(23),可根据下式对广义坐标qr,n+1、广义速度、广义加速度和拉氏乘子λr,n+1进行更新将算法参数

这样,式(26)中的修正项Δqr和Δλr可由下式计算得到

式中,g表示式(20)的残差项

而D表示g关于qr和λr的雅可比矩阵

2 分离释放方案设计

为实现多刚体系统无碰撞分离释放,需要基于具体构型设计分离释放方案.如图6 所示,本文考虑一种2×2×4 的多刚体系统,即每层2 行2 列4 个刚体连接,n=4 层共N=16 个刚体组装而成的系统.为便于分析,按如图6 所示方式对刚体进行编号.规定：层数从上往下分别为第1～4 层.每层卫星编号,以第一层为例,前侧从左至右分别为第1～2,后侧从左至右分别为第3～4 号卫星,以下每层按照第一层顺序依次叠加进行编号.

图6 16 个刚体组成的多刚体系统Fig.6 The MRB system of 16 rigid bodies

基于实际工程应用,依靠较成熟的弹射分离方式进行多刚体的分离释放[33-34].弹射装置对称安装在刚体每个表面的4 个角上,每次分离时,刚体接触面4 个角上的弹射装置同时开启,产生大小相同且方向垂直于刚体接触面的弹射力FT以完成分离机动.需要注意的是,每次分离时的接触面为分离成两部分的物体之间的接触面.比如,第1 层与第2 层分离时,接触面为第1 层下表面和第2 层上表面,并且此时开启的弹射机构位于第1 层下表面4 个角与第2 层上表面4 个角上.

采用不同弹射方向及弹射分离顺序,设计两种无碰撞分离方案,如图7 和图8 所示.第一种即首先利用弹射机构同时分离第1 层刚体与第4 层刚体;随后第1 层刚体与第4 层刚体利用相同方式逐渐分离,即按照图中所示方向先左右分离,再前后分离的顺序;中间第2 层与第3 层按照先前后分离,再左右分离,最后上下分离的顺序进行机动.图中红色箭头表示刚体分离释放方向,序号表示分离次序.因此,全部步骤分为4 步进行,需进行4 次弹射分离释放.需要注意的是,图中所示分离方向为示意图,各刚体姿态与初始构型相同.随着分离进行,分离方向会随着刚体姿态变化,但方向始终垂直于分离面.

图7 第一种分离释放方案Fig.7 The first scheme of separation deployment

图8 第二种分离释放方案Fig.8 The second scheme of separation deployment

第二种方案即首先利用弹射机构分离上三层刚体与第4 层刚体;随后,第1 层与第3 层首先同时以相反方向和第2 层分离,随后这三层再同时左右分离,最后分别前后分离成12 颗独立刚体;第4 层先左右分离分成两列,再前后分离为4 个独立刚体.

上述两种分离方案的分离方式均采取等间隔分离,弹射装置可使弹射力FT的冲击时间持续ΔtT.随后进行ΔtV的自由运动,即每一次分离间隔为Δt=ΔtT+ΔtV.特别地,第一次分离开始于t1,随后每次弹射开始于Δt·(k-1),k表示第k次弹射.

3 算例研究

为研究多刚体系统无碰撞分离释放过程中的非线性动力学及验证分离方式的有效性,本文基于两种分离释放方案进行动力学仿真.设多刚体系统位于400 km 赤道圆轨道.单个刚体质量m=150 kg,尺寸为0.8 m×0.8 m×0.4 m.一般情况下,在轨飞行器由于硬件安装及特殊的结构设计,质量分布并不均匀.因此,考虑一种质量分布不均匀的刚体,质心惯量矩阵为

每个弹射机构可产生推力FT=1500 N,分离释放阶段的时间规划为

为进行刚体分离碰撞分析,首先计算每个时刻每个刚体与其他15 个刚体质心之间的间距,再将每个时刻的间距进行对比,分别取出每个时刻这些间距中的最小值与最大值

式中,dΔ,tk指的是tk时刻下第i和第j个刚体之间的距离的集合,i=1,2,···,N,j=1,2,···,N.dmin,tk和dmax,tk即tk时刻这些间距的最小值与最大值.

图9 和图10 分别给出了两种分离方案中的刚体质心之间最小间距和最大间距的时间历程.最小间距时间历程如图9 所示,图中Case 1 和Case 2 分别代表第一种和第二种分离释放方案,下面称为方案1 和方案2.可以看出,方案1 中最小间距在前15 s内一直维持在0.4 m,这是因为刚体系统没有完全分离,最小间距为刚体包络最小尺寸,即上下连接的两刚体质心间距.而在方案2 中,第5 s 后最小间距变化为0.8 m,这是由于此时系统中没有层间连接存在,最小间距变为同一层中刚体连接的最小间距即刚体的长或宽,符合系统方案2 中的构型变化.两种方案结果相同的是在完全分离后,刚体最小间距持续增大且都大于卫星本体包络尺寸,表明两种方案分离释放过程中及分离之后刚体之间的间距均未小于刚体包络最小尺寸,即没有碰撞发生.此外,完全分离后,方案2 的最小间距上一直大于第一种分离方式.

图10 给出了两种分离方案下的刚体最大间距的时间历程.从10 图可见,两种方案的最大间距在第一次分离后一直增大,并且方案2 的最大间距仅略小于方案1 的最大间距,两种方案的最大分离范围近乎相等.

图9 两种分离方案刚体最小间距时间历程Fig.9 Time histories of the minimum distance between rigid bodies in two separation deployment schemes

图10 两种分离方案刚体最大间距时间历程Fig.10 Time histories of the maximum distance between rigid bodies in two separation deployment schemes

上述结果表明,两种分离方案的分离范围几乎相同,但完全分离后方案2 最小间距较大,分离效果更好.

为研究刚体分离释放过程中的非线性动力学行为,选取5 号刚体作为分析对象,给出其分离过程中位移及欧拉角的变化,如图11～图14 所示.图11 和图12 分别给出了在轨道系O2xyz中两种分离释放方案下5 号刚体的位移时间历程,显示了刚体位移受弹射力影响的变化.可以看出,弹射后刚体逐渐远离初始轨道.此外,第二种分离方案下,5 号刚体位移仅发生3 次变化,这是由于进行第二步分离,即在第5 s开始的弹射分离时,5 号刚体受到大小相等方向相反的弹射力,其运动等价于未受弹射力影响.

图11 第一种分离方案5 号刚体位移时间历程Fig.11 Time histories of displacement of the No.5 rigid body in the first separation deployment scheme

图12 第二种分离方案5 号刚体位移时间历程Fig.12 Time histories of displacement of the No.5 rigid body in the second separation deployment scheme

为研究两种分离释放方案下刚体相对于轨道系的姿态变化,同样选取5 号刚体作为分析对象.为便于直观理解,采用欧拉四元数λ=[λ0λ1λ2λ3]T描述姿态运动,其中

对于5 号刚体,其广义坐标如式(1)所示,则其体轴系3 个主轴单位向量在轨道系中的表示为

显然,5 号刚体体轴系3 个主轴方向单位向量组成的坐标矩阵在体轴系下的表示是单位阵I,在轨道系下表示为坐标矩阵B=[uξuηuζ],两者存在如下关系

因此,四元数可以由B中元素确定为

式中,Bij为B中第i行第j列元素.

图13 和图14 给出了两种分离释放方案下5 号刚体姿态四元数的时间历程.可以看出,由于卫星质量分布不均匀,刚体在弹射力作用下发生连续旋转现象.

图13 第一种分离方案5 号刚体姿态四元数时间历程Fig.13 Time histories of Euler quaternions of the No.5 rigid body in the first separation deployment scheme

图14 第二种分离方案5 号刚体姿态四元数时间历程Fig.14 Time histories of Euler quaternions of the No.5 rigid body in the second separation deployment scheme

最后,给出两种释放分离方案下系统在轨道系中的最终构型,以便更加直观地说明分离情况.如图15 和图16 所示,方案1 的刚体分离后呈无规则分布,而方案2 中刚体有序分布,呈对称构型,有利于分离完成后对刚体进行控制.

图15 第一种分离方案系统最终构型Fig.15 The final configuration of system in the first separation deployment scheme

图16 第二种分离方案系统最终构型Fig.16 The final configuration of system in the second separation deployment scheme

4 结论

本文提出了一种在轨可自主分离的多刚体系统,研究了其分离释放的动力学行为.采用自然坐标方法建立了刚体动力学模型,描述其相对轨道坐标系的运动及姿态变化,并给出了系统连接组合的约束条件.基于弹射分离方式,设计了两种多刚体系统分离释放方案.结果表明,对于质量分布不均匀的刚体,两种分离释放方案均可实现多刚体系统的无碰撞分离,分离后多刚体系统的分离范围逐渐增大,有效避免了碰撞风险.两种分离方案的分离范围几乎相同,第二种分离方案最小间距较大且分离后刚体分布规则易于分离后施加控制,分离效果更好.由于刚体质量分布不均匀且受到分离弹射作用力的影响,分离过程中可以发生连续偏转.若分离完成后对刚体姿态有所要求,需进一步进行姿态控制.