刚-液耦合航天器系统的Hamilton结构及稳定性分析*

易中贵, 岳宝增, 刘 峰, 卢 涛, 邓明乐

(1.北京理工大学 数学与统计学院, 北京 100081;2.北京理工大学 宇航学院, 北京 100081;3.北京宇航系统工程研究所, 北京 100076;4.中国空间技术研究院通信与导航卫星总体部, 北京 100094)

0 引 言

现代航天器通常都需要携带大量的液体燃料推进剂.当这些复杂航天器或者探测器在执行任务时,如姿态机动、轨道转移、交会对接、悬停与避障等,复杂航天器系统的液体推进剂与部件的相互运动将会严重影响该耦合系统的稳定性.

自20世纪60年代开始,研究人员们就发现可以采用刚体来等效封闭储腔内部分填充连续介质液体的动力学行为[1-2].自此,得益于动力学描述的简化、计算复杂度的减小以及在舰实时控制的可实现性等优势,相比于计算流体动力学(computational fluid dynamics,CFD)方法,等效力学模型(包括质量弹簧与平面摆模型[1]、球摆模型[3-6]、运动脉动球模型[7]等)在航天工程领域受到了越来越多工程人员的青睐.然而,对于大幅非平面晃动的情形(即液体首先会经历一个伴随有液体起旋的大幅横向晃动,进而会出现明显的旋转晃动和液体自旋运动,来自Tang和Yue[8]工作的数值结果很好地解释了这一物理过程),Liu和Yue等[9-10]提出了一个复合3D刚体摆模型,该模型被证明是一个有效且可行的方案.通过与解析解、实验解和数值解进行对比,该模型成功预测了球形储腔中的晃动力.在此之前,据研究者统计,还没有发现研究液体起旋的相关工作.由于控制方程中出现的奇异性,即便是球摆模型也不能完成此项工作[10].

几何力学是从现代微分几何(流形)的观点对经典Lagrange和Hamilton力学的现代描述[11].这些现代的、内禀的几何技术为人们提供了一个全局和无坐标描述的方案,这种方案可以避免选择局部坐标时容易遇到的繁冗的计算以及致命的奇异性问题[12],如刚体姿态表示里面的Euler角就是一个例子(这实际上是Lie群中的特殊正交群So(3))[10].

Arnold将著名的Lyapunov稳定性理论向前推广,从而得到能量-Casimir方法(该术语是由文献[13]首次提出的).通过对能量-Casimir函数取一阶导数即可得到相对平衡点,再结合二阶导数以及凸分析即可得到相应力学系统的稳定性[14].Ozkazanc[15]使用Lagrange力学推导了带全充液储腔的航天器系统的动力学方程,并使用能量-Casimir方法分析了该系统的稳定性.该航天器系统被建模为一个刚体携带有不可压缩、无黏、均质液体的全充液系统.因此,该系统实际上没考虑液体晃动的影响.但是该文献对该航天器系统的非正则Hamilton结构、Lie-Poisson描述与Euler-Poincaré描述、相对平衡态以及控制问题都做了非常详细的研究.Ardakani等[16]提出了一种推导储腔内带有自由液面的二维不可压缩旋转流体流运动的变分原理,并采用Euclide群表示刚性储腔的运动.Gasbarri等[17]采用多体建模方法研究了刚-液-柔耦合航天器系统的动力学建模以及稳定性分析,并采用球摆模型来等效储腔内部分填充液体的晃动行为.Salman和Yue[18]使用Lyapunov以及Casimir能量函数研究了充液航天器系统的运动稳定性问题.该文献中,他们采用平面摆来等效液体的晃动问题,因此该模型只适用于小幅线性的横向晃动问题.此外,他们还研究了该耦合系统的分叉以及混沌问题.闫玉龙[19]采用能量-Casimir方法研究了刚-液-柔耦合航天器系统的姿态稳定性问题.其中液体晃动被等效为一个平面摆模型,柔性附件采用的是线性剪切梁.

然而,当人们尝试把能量-Casimir方法推广到几何精确杆(或者三维弹性板或壳)时,则被证明这是行不通的[20].这是由于约化空间里表示的Casimir函数很难表示或实际上根本就不存在[21-22].幸运的是,Simo等[20]针对这一问题做了改进,并由此提出了能量-动量方法.文献[20]研究的另外一个非常关键且重要的结果就是分块对角化技术,它把能量动量函数的二阶变分分块对角化为整体刚性运动和内部振动两块,即将对称群产生的“旋转”扰动与“内部(变形)”扰动的互补空间分离,其中前一个分块即对应著名的Arnold形式.

据此,笔者从几何力学的Lagrange角度出发,系统研究了刚-液-柔耦合航天器系统的全局和无坐标描述动力学方程、相对平衡态的寻找及其相对平衡特性的证明、稳定性准则的建立和分析[23].模型中的柔性附件采用的是可以在三维空间做任意运动的几何精确杆(包括拉伸、剪切、扭转和弯曲).因此采用能量-动量方法与分块对角化技术研究了该刚-液-柔耦合航天器系统的稳定特性.

在本文中,针对刚-液耦合航天器系统液体推进剂的非线性晃动行为,给出了适用于球腔及柱腔的3D刚体摆等效力学模型.由此,采用内禀的几何技术研究了该耦合系统的Hamilton结构,推导了其在约化空间上的约化Poisson括号.最后,采用能量-动量方法以及分块对角化技术,建立了该耦合系统的相对平衡态和相应的自旋稳定性准则.

1 刚-液耦合系统的Hamilton结构

1.1 Lagrange函数

燃料推进剂在球形储腔或者圆柱储腔中的非线性晃动行为可通过一个刚体摆来等效.如图1所示,液体燃料的晃动部分通过刚体摆来等效其力学行为,而未晃动部分(即“冻结”燃料部分)则等效为固定在主刚体上的一个集中质量点.而且当h=0时,对应于球形储腔的等效,否则对应于圆柱储腔的等效.关于等效原则(包括静态属性和动态属性)的叙述可参考文献[1]中的3.2节,较详细的可参考NASA组织研究的,关于液体晃动动力学行为的文献[24]中的第6章.对于球腔和圆柱腔中等效力学模型的各等效参数也可参考文献[1](P48)中提供的经验公式,或者文献[24](P204)中提供的经验公式,两者的区别只是坐标系选取的差异.等效方式可参考文献[10]中2.1节的详细叙述.如果不考虑未晃动部分的“冻结”燃料,那么图1中的力学模型就完全等价于经典的双刚体系统[25-30].

图1 刚-液耦合航天器系统的等效力学模型

图1中的物理参数解释如下:OI为空间惯性坐标系OIE1E2E3的原点;O0为系统质心;O1为主刚体平台在质心处的连体坐标系O1e1e2e3的原点,也是储腔的形心;O2为刚体摆在质心处的连体坐标系O2f1f2f3的原点;Of为“冻结”燃料集中质量点;Oh为等效液体晃动刚体摆的虚拟悬挂点;h为主刚体平台上连体坐标系中,从点O1到点Oh的矢量;l为等效刚体摆上连体坐标系中,从点Oh到点O2的矢量;d为主刚体平台上连体坐标系中,从点O1到点Of的矢量;φ0为空间惯性坐标系中,从点OI到点O0的矢量;φ1为空间惯性坐标系中,从点OI到点O1的矢量;φ2为空间惯性坐标系中,从点OI到点O2的矢量;φf为空间惯性坐标系中,从点OI到点Of的矢量;B1为主刚体平台从连体坐标系到空间坐标系的姿态矩阵;B2为等效刚体摆从连体坐标系到空间坐标系的姿态矩阵.

根据图1,可以得到如下的运动关系:

φ2=φ1+B1h+B2l,

(1a)

φf=φ1+B1d.

(1b)

根据上式与系统质心关系,可得

M=m1+m2+mf,

(2a)

(2b)

式中m1,m2,mf分别为主刚体平台、等效刚体摆和“冻结”燃料的质量.

该刚-液耦合航天器系统的构型流形可由如下定义的Cartesius积表示:

Q=3×So(3)×So(3)={q=(B1,B1,B2)}.

(3)

假设q1∈3是主刚体上的任意一点在连体坐标系O1e1e2e3下的位置矢量,那么其在空间坐标系OIE1E2E3下的空间惯性位置可表示为φq1=φ1+B1q1.因此主刚体动能可写为

(4)

(5)

(6)

同理,根据相同的思路可以推导得到等效晃动刚体摆的动能,并直接写为

(7)

对于“冻结”燃料部分(即未晃动部分的固定质量)的动能也可以根据相同的推导思路直接写为

(8)

其中Jf=mf(‖d‖213-d⊗d).

根据式(5)—(7),便可得到图1中表示的刚-液耦合航天器系统的总动能,并表示如下:

T=T1+T2+Tf=

(9)

由于本文中不考虑系统所处的重力环境,因此系统没有势能项.所以图1中表示的刚-液耦合航天器系统的Lagrange函数L:TQ→可表示为

L(φ1,B1,B2,v1,Ω1,Ω2)=

(10)

1.2 刚-液耦合系统的约化

接下来,考虑特殊Euclide群SE(3)里面的一个群单元g:

(11)

式中A∈So(3),a∈3×1,0∈1×3.其中特殊Euclide群SE(3)可以理解为一个Cartesius积[12],SE(3)=So(3)×3.

根据Lie群在光滑流形上的左作用的定义[11-12],便可定义特殊Euclide群SE(3)在式(3)中介绍的构型流形Q上的左作用,因此可得

Lg(q)=g·q=(Aφ1+a,AB1,AB2).

(12)

通过验算可以发现式(10)中的Lagrange函数在此作用下是不变的,所以该刚-液耦合航天器系统的Hamilton函数在此群作用下也是不变的.由此,就可通过此群对该系统做约化处理.此约化可分为两步:第一步先介绍3约化(即系统的平移不变性),这也对应于系统的总线动量不变性;第二步再介绍So(3)约化(即系统的旋转不变性),这也对应于系统的总角动量不变性.

Q=3×So(3)×So(3)={q=(φ0,B1,B2)}.

(13)

系统的总线动量P可以表示为

(14)

因此根据此式,式(10)中的Lagrange函数就可重新改写为包含系统质心矢量的线动量的函数:

(15)

式(12)中引入的左作用的动量映射就可定义为

(16)

其可通过标准函数Jξ(vq)=〈FL(vq),ξQ(q)〉计算得到[11-12].这是因为由式(15)可知FL(vq)=∂L/∂v0=Mv0=P,因此Jξ=〈P,ξ〉,从而可得J(vq)=P.Legendre变换FL在空间T3×TSo(3)×TSo(3)上诱导了一个辛结构.

在常数点P对应的约化空间为J-1(P)/3=TSo(3)×TSo(3).此时在约化空间上,式(15)中Lagrange函数的第一项就是一个常数项,因此可以去掉,从而可得约化后的Lagrange函数为

(17)

接下来再介绍So(3)约化(即旋转约化),这对应于系统的总角动量不变性.这里的理论依据为Poisson约化理论[25,32].为了获得系统角动量的显式表达式,这里将式(17)中的Lagrange函数改为如下二次形式的函数:

(18)

从而将Legendre变换FL应用在此二次形式的Lagrange函数,可得

(19)

(20)

式中Π=[Π1;Π2].由式(15)可知相对姿态B包含在J-1里的J12中.

先定义So(3)在C=So(3)×So(3)上如下的一个左Lie群作用:

Φ:So(3)×C→C,(R,(B1,B2))(RB1,RB2).

(21)

相对该作用在T*(So(3)×So(3))上的余切提升为

(22)

从而在约化空间T*(So(3)×So(3))/So(3)里的每个等价类的表示可写为

(23)

(24)

是一个Poisson映射,π1=B1Π1,π2=B2Π2称为空间角动量.

FP=F∘P

(25)

定义在空间T*(So(3)×So(3))上的一个函数FP,从而使得约化空间上的函数满足如下关系:

(26)

又由于在T*(So(3)×So(3))上的正则Poisson括号[1]的定义如下:

(27)

其中函数HP∈C∞(T*(So(3)×So(3))).从而根据约化Poisson结构的定义可得

(28)

{F,H}(Π1,Π2,B)=

(29)

1.3 刚-液耦合系统的Poisson括号

在继续推导上式中各项的具体表达式之前,先引入特殊正交群So(3)的二次切丛及其对偶空间上任意元素全局表示的概念.

根据文献[33]中5.2.3小节的介绍,二次切丛TTSo(3)及余切丛T*TSo(3)上任意元素的全局表示可分别写为

(30a)

(30b)

式中a,u,b,w,Ω均为3中的单元,而群元素R∈So(3).根据此公式,现在就可以推导TT*(So(3)×So(3))与其对偶T*T*(So(3)×So(3))上的任意单元的全局表示,具体推导步骤如下:

(31)

(32)

(33)

(34)

结合式(33)及式(34),即可得到

(35)

{F,H}(Π1,Π2,B)=

(36)

该式可进一步简化为

{F,H}(Π1,Π2,B)=

(37)

根据文献[33]中的式(2.13)和上式,可得图1描述的刚-液耦合航天器系统的如下约化Poisson括号形式的动力学方程:

(38)

此外需要说明的是,式(38)中推导的约化Poisson括号可以推广到带柔性附件的情况.针对线性剪切梁的情形[33],只需将式(20)中描述的系统的构形流形改为C=So(3)×So(3)×M[34],其中M表示从区间[0,L]到的函数,L表示梁的长度.根据文献[33]的推导过程,刚-液-柔耦合航天器系统的约化Poisson括号可写为其中式(5.50)的形式.而针对更为一般的非线性弹性体(包括杆和板等)的Hamilton结构,感兴趣的读者可以阅读经典文献[35].对于这里提到的带有线性剪切梁的柔性附件的情形,文献[14,19,36]对能量-Casimir方法的构造与其在相对平衡态的稳定性分析上均做了非常详细的介绍和研究.

2 刚-液耦合系统的稳定性分析

本节将采用能量-动量方法研究刚-液耦合航天器系统的稳定性特征.

2.1 系统的相对平衡态

文献[20]通过构造一个能量-动量函数,为描述相对平衡点提供了一种变分方案,这点可从其中2.4小节“相对平衡理论”的式(2.22)中得以体现.重构能量-动量函数便可得到修正势能的定义,而此修正势能的极值点恰恰就是相对平衡点.这一事实被称为对称临界原理.本小节将采用此对称临界原理来推导刚-液耦合航天器系统的相对平衡点.为了探讨系统自旋角速度对系统稳定性的影响,这里在系统主刚性平台上安装一个匀速旋转转子,此时的系统称为陀螺力学系统,更多的介绍可参考文献[33]中的式(2.81)和式(5.83).

根据式(13)和(17)可知,系统的构形流形可重新表示为Q=So(3)×So(3)={q=(B1,B2)},并带有如下定义的Riemann度量:

《vq,vq》=〈vu,Jvu〉,

(39)

系统中剩余的陀螺场及对称群分别为[23]

(40)

式中矢量y1,y2的表达式由下式确定:

(41)

式中Jd,ωφ分别表示转子的常数惯性并矢及转动角速度.

从而Lie群G=So(3)左作用于该构形流形Q时有

Ψ:G×(So(3)×So(3))→(So(3)×So(3)),(R,(B1,B2))(RB1,RB2).

(42)

(43)

根据能量-动量方法的分块对角化技术[33],定义如下的锁定惯性张量Ilock(q):g→g*及在g*中诱导出的陀螺-动量IY(q):

(44)

再根据文献[33]中式(2.98)的定义,便可写出上述刚-液耦合航天器系统的如下修正势能项:

(45)

(46)

根据前面关于对称临界原理的叙述,可知式(45)的极值点即为系统的相对平衡点,因此这里需要对修正势能Vξ(q)进行微分.根据切空间TqQ上的矢量在Q上生成的曲线可诱导出如下微分:

(47)

将式(46)代入上式即可得到

(48)

式中s=εcd-εah,并且这里假设J1是对称的,后文2.2小节中式(53)给出的相对平衡点可以验证该假设是成立的.由此式即可得到系统平衡态的构形流形(B1,e,B2,e)满足如下的平衡态条件:

(49)

系统处于该相对平衡态时,可以证明:

1)图1中虚拟悬挂点矢量的空间表示B1,ehe、摆杆矢量的空间表示B2,ele和系统空间角速度ξ三者之间还满足如下的共面条件:

B1,ehe×B2,ele·ξe=0.

(50)

该式的详细证明过程可参考文献[23]中定理3.3条件(i)的证明, 并且可以发现该共面条件不受陀螺场的影响.

2)如果不考虑动量轮,则系统空间角速度ξe就是锁定惯性并矢I12f,e的特征矢量,满足条件

ξe×I12f,eξe=0,

(51)

即I12f,eξe=λξe.如果考虑动量轮,那么上式应改为

ξe×(I12f,eξe+B1,eκe)=ξe×αe=0,

(52)

式中αe表示系统在相对平衡点处的总角动量.该条件的具体证明过程可参考文献[23]中定理3.3条件的证明,其中锁定惯性并矢I12f,e表达式的定义可参考文献[23]中的式(3.33),不考虑柔性附件部分.

2.2 系统的稳定性分析

结合能量-动量方法以及2.1小节中的相对平衡点的介绍和讨论,接下来将应用这些方法来研究图1中描述的刚-液耦合航天器系统的稳定性.

本文将研究T型构形均匀旋转的特殊情形.假设刚-液耦合航天器系统中其他等效部件的位置构形是共线的,即 “冻结”质量点的方向与等效刚体摆的摆杆方向是共线的.系统空间角速度的旋转轴与等效刚体摆的摆杆方向垂直,且通过系统质心.即对应于图1中摆杆方向位于e1轴正方向时构成的状态,旋转轴沿e3轴方向.相关参数设置如下:

(53)

式中{e1,e2,e3}, {f1,f2,f3}分别对应于图1中主刚体平台、等效刚体摆的连体坐标系.由上式可看出,摆杆与悬挂高度处于一种类似于“折叠”的状态.另外可以验证这些参数是满足式(49)中的相对平衡方程的.

[(J1,e,3+J2,e,3+2εcdl-2εahl)ξ+κ]B1,ee3,

(54)

这里假设式中各惯性量J1,J2,J12均只在各自的主对角线上存在元素,如J1=diag[J1,1,J1,2,J1,3].

(55)

V={δq∈TqeQ|《δq,ηQ(qe)》=0,∀η∈gμe}=

{(J1,e,3+εcdl-εahl)u1,e,3+(J2,e,3+εcdl-εahl)u2,e,3=0}.

(56)

再根据文献[33]关于分解空间VRIG的定义即可得

(57)

根据文献[21]中式(2.31)的叙述可知,在相对平衡点处,限制在空间V上的修正势能的二次微分的正定性则隐藏着形式稳定性,即D2Vξ(qe)|V×V>0.因此,由式(45)和式(46)可得

(εcdl-εahl)ξ2(u1,e,3-u2,e,3)2≥

(58)

式中的系数σi(i=1,2,…,5)分别对应于前一不等式中各变量前的系数.具体推导过程与式(47)相似,这里将不再赘述.

同理,为了后续计算简便,这里也对式(56)中的等式做如下简化:

γ1u1,e,3+γ2u2,e,3=0,

(59)

其中系数γ1,γ2分别对应式(56)中等式里各变量前的系数.

从而将上式代入式(58),展开得到的表达式,并收集同类项后可得限制在空间V上的修正势能的二次微分为

(60)

综上所述,将前面各式中各相应的系数代入上式中替换对应的参数即可得到刚-液耦合航天器系统的如下自旋稳定性条件:

(61)

式中σ5=εcdl-εahl,而各惯性参数以及各约化质量可参考式(15)里的具体定义.

这里需要强调的是,假设上述等效刚-液航天器系统中仅由主刚体平台构成,那么上式中的稳定性条件即变为如下经典的刚体稳定性准则:

(62)

此时对应于l=0.对于仅有等效刚体摆时同理.该刚体稳定性准则可从文献[37]中采用的保结构算法计算得到的动量球上的余伴随轨道(图8)得以清晰体现.其中式(16)解释部分提到的辛结构[38]是力学系统中经常讲的保结构[39]的一种内蕴的几何结构,其可用于构造相应力学系统优异的数值迭代格式[40-42].

(63)

(64)

结合文献[33]中的式(2.142)以及本文中的式(54),即可得到

(65)

从而结合上面两式即可得到式(63)中的矩阵表达形式:

(66)

由于耦合系统的空间旋转角速度ξ旋转轴与e3轴的方向相同,从而应将上式中的主对角线上的第三项去掉.因此该耦合系统的稳定性就等价于要求上式由剩下的项构成的矩阵是正定的,即要求主对角线上的第一项以及第二项要同时大于零.再将式(64)和式(65)代入上式,即可得到刚-液耦合系统的如下Arnold形式的稳定性条件:

(67)

通过与式(61)对比容易发现,该Arnold形式的稳定性条件包含在式(61)中表述的稳定性条件中.该稳定性条件不能保证系统的稳定性,即不能保证式(61)中的稳定性条件成立.因此,对于刚-液耦合航天器系统的稳定性分析,不能将耦合系统的运动视为一个整体刚性运动,必须考虑耦合效应对稳定性的影响.

2.3 稳定性的数值仿真

根据2.2小节中式(61)给出的刚-液耦合航天器系统的自旋稳定性条件,本小节将采用图形的方式来体现该系统在具体参数下的稳定域.这里将液体推进剂的非线性晃动行为通过一个高度与半径相等的等体积圆柱体来等效,如图1所示.

模型中各参数的具体数值设置如下,其中部分参数引自文献[23],部分数据通过文献[1]中(P48)给出的经验公式而得

(68)

式中,ρL表示液体推进剂的质量密度,pr与ph分别表示圆柱储腔的半径与高度.式中的坐标系如式(53)所示.此外,液体燃料的晃动部分与“冻结”部分的质量满足:它们的质心与刚体平台的质心重合.

由式(61)可知,耦合系统的稳定性主要受液体推进剂的充液比Ξ、系统的空间旋转角速度ωs以及动量轮转速ωφ的影响,仿真后的稳定域如图2中的阴影区域所示.

图2 刚-液耦合航天器系统的稳定域

图2(a)是系统充液比与系统自旋角速度的关系,它由式(67)Arnold形式的稳定性条件绘制而得.由此可看出,随着充液比的逐渐增加,系统临界自旋角速度是逐渐增大的,即可允许的自旋角速度空间是逐渐增大的.同理,将式(68)中的系统参数代入式(61)绘制而得的稳定性边界为一带状区域ωs≤1 rad/s,通过与图2(a)中的临界值进行比较可知,这也再次验证了Arnold形式的稳定性条件是包含在式(61)中的稳定性条件中的.而图2(b)是系统充液比与转子角速度的关系,由图可知,只要转子角速度在图中曲线之上系统均是稳定的,或者说在不同充液比条件下转子需要满足的最小转速.

3 结 论

本文从几何力学角度出发,针对刚-液耦合航天器系统的3D刚体摆等效力学模型的动力学问题,系统推导了该等效模型的Hamilton结构,详细研究了系统的自旋稳定性特征.论文首先介绍了系统的平移不变性约化和旋转不变性约化,根据正则Poisson括号推导了该刚-液耦合航天器系统在约化空间上的约化Poisson括号形式的动力学方程.接着根据对称临界原理推导了该系统的相对平衡态,并且发现在平衡态条件下:悬挂点矢量、摆杆矢量和系统自旋角速度三者在空间上是共面的特性,且该特性与陀螺项无关.在平衡点条件下,根据能量-动量方法和分块对角化技术推导了系统的自旋稳定性条件与Arnold形式的稳定性边界,并由此可清晰地看到耦合效应对稳定性的影响.最后以图形的方式给出了具体模型参数下的自旋稳定域.