燃料消耗下充液航天器等效動力學建模與分析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)貯腔偏心布放時,航天器在執行軌道機動過程中貯腔內液體易發生劇烈而且形式復雜的晃動行為,進而可能造成航天器剛體運動的不穩定.