燃料消耗下充液航天器等效动力学建模与分析1)

刘峰 岳宝增 马伯乐 申云峰

*(北京理工大学宇航学院,北京 100081)

†(郑州工业应用技术学院信息工程学院,郑州 451199)

引言

在轨运行的充液航天器,通常处于较低重力(低重力或微重力)环境,当它执行轨道或姿态机动、交会对接及软着陆等任务时,贮腔内的液体可能表现出以下几种不同的主要运动模式的叠加运动:(1)横向晃动,即液体自由液面交替抬升、下沉而弓起的不流动的波动行为[1],这种晃动模式具有依赖于贮腔形状、尺寸和加速度环境等因素的固有频率和固有模态,且第一阶横向晃动主模态呈现为反对称形式;(2)旋转晃动,当贮腔受到的激励频率接近于最低阶液体晃动固有频率时,将弓起横向晃动失稳,此时横向晃动的节线发生旋转而出现旋转式流动的波[2];(3)液体自旋,即液体发生的绕其对称轴的类似于刚体的自旋运动.此外,当航天器受到的有效重力的方向不断变化或者航天器进行大角度姿态机动时,液体相对于贮腔晃动的平衡位置将时刻变化,即贮腔内的液体会发生相对于贮腔的整体性刚体运动.

虽然已有不少关于充液航天器动力学与控制方面的研究[3-26],但是还有很多新的问题需要进一步研究和解决,比如:(1)大范围机动下充液航天器发生轨道-姿态-晃动与控制系统之间的动力学耦合问题,已有的研究大多致力于研究充液航天器小角度姿态机动下的动力学建模与控制,而忽略轨道运动的影响,也未考虑大角度姿态运动下航天器系统惯量特性和液体晃动特性变化等因素; (2) 在轨微重力(或低重力)环境下充液航天器机动时,贮腔内液体可能发生复杂的非线性晃动行为,有相当研究考虑了液体大幅晃动问题,但是对于液体旋转晃动以及液体起旋时的航天器刚-液或刚-液-控耦合动力学问题鲜有研究; (3) 多贮腔布局下,不同贮腔内液体的物理特性不同,以及不同贮腔尺寸和充液比等差异导致晃动固有特性差异时,多贮腔液体组合晃动对航天器系统运动的影响尤为复杂; (4)燃料消耗弓起液体晃动特性改变时液体晃动对航天器系统动力学的复杂影响有待进一步研究.

本工作在已完成工作[27]的基础上,考虑液体燃料消耗弓起液体晃动固有特性变化,进一步提出变参数的刚体摆复合模型,模拟航天器球形贮腔内的液体晃动动力学行为,要研究燃料消耗情形下液体发生非线性晃动时充液航天器的耦合动力学问题.具体地,本工作将基于变参数的刚体摆复合模型,采用混合坐标意义下的拉格朗日方程对一类变质量充液航天器展开等效动力学建模;随后,考虑该充液航天器执行大范围机动,以大角度三轴稳定姿态机动和零冲量大范围轨道机动为例,展开变质量充液航天器大范围机动仿真和耦合动力学特性分析,为航天工程实践提供一定的理论指导和参考.

1 考虑燃料消耗的刚体摆复合模型

作者前期的工作[27]针对常/低重力环境下部分充液球形贮腔内液体非线性晃动(包括大幅横向晃动、旋转晃动以及液体自旋运动等)而提出了一种三自由度刚体摆复合模型.如图1 所示,该复合模型包含如下两部分:一是悬挂点H 位于贮腔形心的刚体摆P,它等效的是参与晃动的那部分液体,其集中质量参数mp等于参与晃动的液体质量,摆长参数lp确定的是参与晃动的液体的质心位置,它也是表征液体晃动固有频率的参数,而刚体摆关于其轴向的转动惯量Jp等效的是参与液体旋转(液体发生的旋转晃动和自旋运动统称为液体旋转) 的那部分液体所具有的转动惯量;二是集中质量点Q,它等效的是未晃动的那部分液体,其集中质量参数m0等于未晃动液体的质量,其相对于贮腔形心的距离参数l0确定的是未晃动液体的质心位置.

图1 液体晃动刚体摆复合模型及其等效参数Fig.1 The rigid-pendulum composite model for liquid slosh and its equivalent parameters

传统的弹簧-质量模型和单摆模型只考虑了液体的低阶横向晃动,球摆模型能够在研究最低阶横向晃动的同时对液体的旋转晃动加以描述; 而作者前期工作[27]中提出三自由度刚体摆复合模型不仅能够在保证力系等效性的前提下具有尽可能少的模型参数,而且能较为全面地描述航天器内液体可能表现出的相对于贮腔的多种不同运动模式的叠加运动.具体地,三自由度刚体摆复合模型的运动与贮腔内实际液体的运动之间具有以下几点相似性:一是模型中的集中质量Q(代表未晃动质量)和刚体摆P的集中质量(代表晃动质量)随着静力平衡位置变化(称为动平衡位置,可能是由于航天器发生大范围姿态运动或者等效重力加速度环境改变等因素弓起的)而发生的相对于贮腔的任意转动,以此可描述贮腔内的液体发生相对于贮腔的整体性刚体运动; 二是刚体摆P 在动平衡位置附近的有限幅度(可以是小幅的,也可以是大幅的)的往复摆动则类似于贮腔内液体自由液面的横向晃动; 三是当刚体摆P 作锥形运动时,所描述的即是液体的旋转晃动;最后是刚体摆P 以Jp为转动惯量绕其摆轴方向的自旋运动描述了液体的刚体自旋.后续对该复合模型的发展作进一步研究,使其适用于研究液体燃料消耗过程中的液体晃动动力学问题.

首先,液体燃料消耗将直接导致贮腔充液比β 减小(此时β 为时变量),进而导致参与晃动和未晃动两部分液体的质量、转动惯量及其质心位置发生变化,而这将体现在刚体摆复合模型中的等效参数具有时变性上.文献[2]中总结了常规重力或低重力环境下球形贮腔内液体晃动等效摆的质量参数及摆长参数随充液比β 的变化规律.本文通过提取文献[2]中的实验结果并将其进行拟合,如图2 和图3 所示,得到刚体摆复合模型中刚体摆的集中质量参数mp和摆长参数lp随充液比β 变化的近似关系,具体的表达式分别为

式中,R是球形贮腔的半径,ρ 是贮腔内液体的密度.

图2 刚体摆的集中质量参数随充液比的经验关系Fig.2 Empirical relationship between the lumped mass of the rigid-pendulum and the liquid-fill ratio

图3 刚体摆的摆长参数随充液比的经验关系Fig.3 Empirical relationship between the length of the rigid-pendulum and the liquid-fill ratio

然后,依据文献[27] 中指出的针对该刚体摆复合模型的参数配置原则,可得到集中质量点Q 的质量参数m0的具体表达式

以及集中质量点Q 相对于贮腔形心的距离参数l0满足如下关系式

式中,mliq是贮腔内所有液体的质量,代表液体不晃动时所有液体的质心相对于贮腔形心的距离.

最后,需要确定关于刚体摆轴向的转动惯量等效参数Jp.考虑到实际情况下参与液体旋转的液体所具有的转动惯量是一个极难确定的物理量,尤其在液体起旋阶段,参与液体旋转的那部分液体所具有的转动惯量有着比较复杂的随时间变化规律[27].这里为了简化问题,假设将液体固化后(可将液体视为一个刚体) 液体所具有的关于其对称轴的转动惯量Jsolid等于刚体摆的转动惯量等效参数Jp,即Jp=Jsolid.

2 充液航天器刚-液耦合动力学模型

采用混合坐标意义下的拉格朗日方程[28-29]来建立如图4 所示充液航天器的耦合动力模型,其贮腔内的液体燃料采用刚体摆复合模型进行等效,同时考虑液体燃料消耗的影响.图4 中,N:{n1,n2,n3}是牛顿惯性参考系,g是航天器系统受到的有效重力加速度,方向始终与-n3同向;B*是航天器主刚体(指除去液体的部分) 的质心;B:{b1,b2,b3} 是原点位于贮腔形心且与航天器主刚体固连的本体坐标系,并假定B:{b1,b2,b3}与航天器主刚体的惯性主轴坐标系(其原点位于主刚体质心) 对齐(即两个直角坐标系的对应坐标轴同向).

图4 带部分充液球形贮腔航天器及其等效物理模型示意图Fig.4 Illustration of a spacecraft with a partially liquid-filled spherical tank

需要特别指出的是由于液体燃料消耗的影响,刚体摆复合模型的所有等效参数是时变的,即整个充液航天器本质上是一类极为复杂的变参数系统.为简化后续推导,可作以下简化假设:认为所有模型参数的变化率(即时间导数)相比于系统的运动为慢变量,假设燃料消耗过程中,模型参数都是准静态变化的.因此,本文中建立的变质量充液航天器耦合动力学模型适用于燃料消耗速率比较低的情形.

2.1 系统运动描述

2.1.1 航天器主体的姿态运动

航天器在惯性参考系下的姿态运动通过本体坐标系B相对于惯性参考系N的转动来描述,转动角速度矢量可以表示为

采用欧拉四元数ε=[ε0,ε1,ε2,ε3]T来描述航天器的当前姿态[7-9],则姿态运动的发展方程为

此外,航天器的姿态矩阵(即惯性参考系N到本体坐标系B的坐标转换矩阵) 在欧拉四元数下的表达式为

2.1.2 航天器主体的轨道运动

航天器在惯性参考系下的平动运动(即轨道运动)通过本体坐标系原点B相对于惯性参考系N的位置和速度来描述,相应的位置矢量和速度矢量可分别表示为

由矢量叠加原理,航天器主刚体质心B*相对于惯性参考系的位置矢量为

式中,pBB*是由航天器整体布局决定的常矢量,设pBB*=Xb1+Yb2+Zb3.于是,航天器主刚体质心B*相对于惯性参考系的速度为

2.1.3 等效液体晃动刚体摆复合模型的运动

刚体摆复合模型中等效未晃动液体的集中质量点Q,它相对于惯性参考系N的位置矢量为

式中,pBQ始终指向重力加速度方向,于是有pBQ=-l0n3.

对式(12)两端同时求惯性参考系下对时间的微分,可得到集中质量点Q 相对于惯性参考系N的速度为

另外,刚体摆复合模型中等效晃动液体的刚体摆P 的集中质量,它相对于惯性参考系N的位置矢量为

沿用文献[27]中采用3 个欧拉角φ,θ 和ψ 来描述刚体摆相对于贮腔运动的方法,于是有

式中,a1=cos ψ sin θ+sin ψ sin φ cos θ,a2=sin ψ sin θcos ψ sin φ cos θ,a3=cos φ cos θ.

对式(14)两端同时求惯性参考系下对时间的微分,可得到刚体摆P 的集中质量相对于惯性参考系N的速度,为

2.2 系统动能

充液航天器系统的总动能包括以下3 类运动所具有的动能:航天器主刚体平动(即航天器轨道运动)和转动(即航天器姿态运动)具有的动能、代表未晃动液体的集中质量点Q 所具有的动能以及代表晃动液体的刚体摆P 所具有的动能.设航天器主刚体的质量为mhub,主刚体关于其惯量主轴坐标系的惯性张量为Ihub,则该航天器系统总动能的定义式为

将式(5)、式(9)、式(11)、式(13)和式(16)代入式(17) 中,即可推导出航天器系统总动能的具体表达式.为方便后续推导,将得到的系统总动能表达式改写成矩阵表示的形式; 并规定所有参与运算的矢量和二阶张量的坐标分量均统一在本体坐标系下给出.于是,可将系统总动能最终表示成如下的矩阵计算式

2.3 系统势能

充液航天器系统的势能为整个航天器在重力场中所具有的重力势能,其表达式如下

将式(10)、式(12)和式(14)代入式(19),可以得到如下矩阵形式的系统势能表达式

式中,c3是姿态矩阵NCB的第三列元素组成的列向量,即

2.4 耦合系统的拉格朗日方程

充液航天器系统的拉格朗日函数有如下表述

式中,zB,ε和ρ均属于广义坐标,是广义速度;而v和ω是在本体坐标系(非惯性系)下给出的速度和角速度的坐标分量矩阵,都不是真实坐标的导数,因此称它们为拟速度[30].

求解航天器轨道运动的拉格朗日方程为

式中,Fnc是航天器系统受到的非保守力(如推力、外界干扰力等)在本体坐标系下表示的坐标矩阵,其作用点位于本体坐标系原点B.当不计入航天器主刚体质量时,由牛顿第二定律,Fnc即表示液体对航天器主刚体的等效作用反力.

求解航天器姿态运动的拉格朗日方程为

式中,Tnc是航天器系统受到的关于本体坐标系原点B的非保守力矩(如控制力矩、外界干扰力矩等)在本体坐标系下表示的坐标矩阵;Tg是航天器系统受到的重力梯度力矩(一类保守力矩,当航天器处于空间微重力环境下,航天器机动过程中重力梯度力矩相比于控制力矩通常很小,此时可以忽略重力梯度力矩的影响)在本体坐标系下表示的坐标矩阵.当不计入航天器主刚体的惯量时,由力矩平衡关系,Tnc即表示液体对航天器主刚体的等效作用反力矩.

求解贮腔内液体晃动等效刚体摆运动的拉格朗日方程为

式中,Td是为了等效液体晃动时产生的黏性耗散而作用于刚体摆悬挂点的阻尼力矩,这里采用如下线性的阻尼力矩模型

2.5 充液航天器动力学方程组

将式(18)和式(20)代入式(21)后,再分别代入式(22)～式(24),可得到充液航天器轨道-姿态-晃动全耦合系统的动力学方程

3 仿真计算与分析

对充液航天器大角度姿态机动和姿态受控稳定下的轨道驱动分别进行相应的MATLAB 仿真.假定:航天器主刚体的质量为mhub=20 kg,主刚体关于其惯性主轴坐标系的惯性张量为Ihub=diag{4,6,5} kg·m2; 贮腔半径为R=0.25 m,贮腔布放位置参数为BpBB*=[0.2,-0.3,-0.5]Tm; 贮腔内液体为燃烧剂甲基肼,质量密度为ρ=874.4 kg·m-3;刚体摆的等效阻尼参数取为βφ=βθ=0.05,βψ=0;当地重力加速度大小为g=1.0 m·s-2.

燃料消耗过程中,假设贮腔充液比β 随时间的变化规律满足如下连续可导的分段函数式

式中,β1和β2分别是燃料消耗前后的充液比,这里取β1=0.6,β2=0.4;T0是燃料消耗时长,用来表示燃料消耗速度快慢的物理量.

3.1 大角度三轴稳定姿态机动

采用常用的比例微分控制(PD 控制) 对该航天器进行三轴稳定姿态机动.具体地,假定初始时刻航天器处于静止状态,即ω(t=0)=0;航天器初始姿态的欧拉四元数为

采用如下的PD 控制策略使航天器机动到目标状态

式中,和ωtar是目标状态下航天器姿态的欧拉四元数的矢量部分和航天器角速度的坐标分量矩阵,这里取=[-0.221,0.074,0.442]T,ωtar=0;Kp和Kd是控制反馈增益矩阵,设计为

此处kp=0.05,kd=0.3.

同时,航天器姿态控制系统输出的控制力矩还需要对航天器受到的重力梯度力矩进行补偿.于是,姿态控制系统输出的控制力矩为

另一方面,为克服重力,航天器还需要提供持续的推力,于是有

此外,假设初始时刻液体存在微小的残余横向晃动,对应地,MATLAB 仿真计算中取刚体摆运动的初始条件为φ(t=0)=2π/180 rad (其他状态变量均为零初始条件).本算例中,假设燃料消耗时长为T0=20 s.

图5～图9 给出了航天器进行大角度姿态机动时,贮腔内液体燃料在初始20 s 内从充液比β1=0.6消耗到充液比β2=0.4 的情形下,主刚体和等效液体晃动刚体摆的运动响应结果.

图5 航天器姿态运动的角速度响应Fig.5 Angular velocity response of the spacecraft

图6 航天器轨道运动的速度响应Fig.6 Orbital velocity response of the spacecraft

图7 航天器主刚体在惯性参考系中的位置Fig.7 Position of the spacecraft hub in the inertial reference frame

图8 等效液体晃动刚体摆的运动响应Fig.8 Motion response of the rigid-pendulum for liquid slosh

图9 刚体摆的集中质量的运动轨迹在B-b1b2 平面内的投影Fig.9 Trajectory of the lumped mass of the rigid-pendulum on B-b1b2 plane

图5 给出了航天器进行大角度三轴稳定姿态机动时的角速度响应,结果显示贮腔内液体晃动对航天器3 个方向的姿态运动均产生了局部振荡形式的扰动影响.另一方面,图6 显示了液体相对于贮腔的运动(包括液体整体性的刚体运动和液体晃动)也会弓起航天器主刚体产生微小的平动速度,但是由于燃料消耗的影响,导致机动过程中航天器轨道运动的最终速度不会收敛到零(这并不违背航天器系统动量守恒),即航天器主刚体在惯性参考系中的位置偏移会不断地增大(图7).这特别提示了燃料消耗下充液航天器在执行交会对接等空间定点任务时,应当考虑液体晃动的影响.

此外,图8 给出了本工况下描述等效液体晃动刚体摆运动的三个欧拉角的时程响应.将图8 的计算结果的数据进行转化,可以得到刚体摆的集中质量(它代表的是晃动部分液体的质心) 在贮腔坐标系中的运动轨迹,图9 是该轨迹在B-b1b2平面内的投影(xp和yp表示刚体摆的集中质量在B-b1b2平面内的投影点的坐标),结果表明此工况下晃动部分液体的运动形式表现为相对于贮腔的整体性的大范围刚体运动和横向晃动的叠加.

3.2 零冲量轨道机动

进一步研究变质量充液航天器姿态受控稳定时在零冲量轨道驱动力作用下的动力学响应.

假设航天器在受到方向始终沿惯性参考系n1轴向,大小(正值表示沿n1轴正方向,负值表示沿n1轴负方向)如图10 所示的零冲量阶跃式驱动力的作用进行轨道机动.

图10 零冲量阶跃式驱动力Fig.10 A step type propulsive force with zero impulse

此外,为了保证航天器在进行轨道机动时的姿态能够保持基本稳定,采用3.1 节中的PD 控制(式(26))来实现航天器的姿态稳定,并取kp=0.5,kd=3.这里,航天器的初始姿态和目标姿态所对应的四元数相同,假设为εv(t=0)==0,且航天器的初始角速度和目标角速度均为零.本算例中,假设燃料消耗时长为T0=70 s.

图11～图16 给出了贮腔偏心布放(贮腔形心偏离航天器主刚体质心)时,航天器在进行零冲量轨道机动过程中并保持燃料持续消耗70 s 的情况下,主刚体和等效液体晃动刚体摆的运动响应结果; 图17给出机动过程中为维持航天器姿态稳定所施加的PD控制力矩.

图11 航天器姿态运动的角速度响应Fig.11 Angular velocity response of the spacecraft

图12 航天器姿态的欧拉四元数响应Fig.12 Euler quaternions response of the spacecraft attitude

图13 航天器轨道运动的速度响应Fig.13 Orbital velocity response of the spacecraft

图14 航天器主刚体在惯性参考系中的位置Fig.14 Position of the spacecraft hub in the inertial reference frame

图15 等效液体晃动刚体摆的运动响应Fig.15 Motion response of the rigid-pendulum for liquid slosh

图16 刚体摆的集中质量的运动轨迹在B-b1b2 平面内的投影Fig.16 Trajectory of the lumped mass of the rigid-pendulum on B-b1b2 plane

图17 姿态稳定PD 控制力矩Fig.17 PD control torque for attitude stabilization of the spacecraft

仿真结果显示:航天器在轨道机动过程中,一方面,航天器在如图17 所示的控制力矩的作用下,航天器的姿态运动虽然受到了液体晃动的扰动影响(图11),但是航天器的整体姿态基本上能够维持稳定(图12); 另一方面,航天器在图10 给出的零冲量阶跃式驱动力的作用下,液体晃动不仅会对航天器轨道平动造成明显扰动影响(图13 和图14),而且由于燃料消耗的影响,会造成该充液航天器在进行零冲量轨道机动时,最终的轨道平动速度无法收敛到零,尤其是在驱动方向(惯性参考系n1轴方向) 上航天器将最终获得一个较大的速度(图13),因此,在零冲量轨道驱动力作用下,由于燃料消耗的影响,虽然航天器沿惯性参考系n1轴方向作了大范围平移,但是航天器此时仍无法保持静止状态,此后航天器将不断地偏离预期的机动位置(图14),导致机动任务失败.

特别值得注意的是,在以上工况下,贮腔内液体将发生极为复杂的非平面运动(图15 和图16),贮腔内液体会经历逐步起旋并发生明显的旋转晃动(图16),此过程中也伴随着液体横向晃动,但由于液体阻尼耗散等因素,液体的横向晃动最终会衰减并消失,贮腔内液体运动最终表现为恒定角速率的刚体自旋运动(图15).

4 结论

本文采用变参数的刚体摆复合等效力学模型模拟燃料消耗下球形贮腔内的液体非线性晃动,并基于混合坐标下的拉格朗日方程,建立了变质量充液航天器轨道-姿态-晃动全耦合的动力学模型.通过编制相应的MATLAB 仿真计算程序,对本文所建立的动力学模型进行求解,研究了该变质量充液航天器在执行大角度三轴稳定姿态机动以及姿态受控稳定下的零冲量轨道机动过程中的动力学响应.得到了对航天工程实践有一定指导意义的两点主要结论.

(1)液体相对于贮腔的运动会造成航天器主刚体位置发生偏移,当航天器在执行零冲量机动时,燃料消耗会造成航天器的轨道平动速度无法收敛到零,因此,在充液航天器执行交会对接等空间定点任务时,应当充分考虑液体晃动的影响.

(2)贮腔偏心布放时,航天器在执行轨道机动过程中贮腔内液体易发生剧烈而且形式复杂的晃动行为,进而可能造成航天器刚体运动的不稳定.