劉 峰，岳寶增，唐 勇

(北京理工大學(xué)宇航學(xué)院，北京 100081)

0 引 言

隨著現(xiàn)代航天器內(nèi)的液體推進(jìn)劑的質(zhì)量占比不斷提高，液體晃動(dòng)動(dòng)力學(xué)問題和有關(guān)的控制研究已經(jīng)成為航天器總體設(shè)計(jì)中的核心問題之一[1-2]。在軌充液航天器內(nèi)的液體晃動(dòng)特別是大幅液體晃動(dòng)具有比較低的固有頻率，很有可能與航天器上的大型柔性附件及控制系統(tǒng)發(fā)生動(dòng)力學(xué)耦合共振，這可能導(dǎo)致航天器結(jié)構(gòu)的破壞或者姿控系統(tǒng)失效而使航天器姿態(tài)翻滾[3-4]。因此，現(xiàn)代航天工程中迫切需要建立準(zhǔn)確、高效且適用性強(qiáng)的液體大幅晃動(dòng)類航天器耦合系統(tǒng)動(dòng)力學(xué)模型。

考慮到推進(jìn)劑包含燃燒劑和助燃劑兩種以上不同液體，航天器上往往需要攜帶多個(gè)卡西尼腔(或者球腔)。與單腔不同，各個(gè)貯腔尺寸及安裝位置差異、推進(jìn)劑消耗不均、不同液體特性等諸多因素，使得多個(gè)貯腔內(nèi)的液體晃動(dòng)與航天器主體之間的動(dòng)力學(xué)耦合特性更為復(fù)雜。而關(guān)于這方面現(xiàn)有的研究極少，更多的研究工作有待進(jìn)一步的完善。本文所建立的四充液貯腔航天器耦合動(dòng)力學(xué)模型具備較強(qiáng)的適用性，模型中不僅考慮了航天器主體軌道平動(dòng)與姿態(tài)轉(zhuǎn)動(dòng)，而且還包含有諸如貯腔的幾何尺寸、布放位置、充液比及液體物理特性等反映耦合系統(tǒng)動(dòng)力學(xué)特征的重要的工程參數(shù)。

Berry和Tegart[5]首次提出了等效液體晃動(dòng)質(zhì)心面模型，并基于此模型編制了相應(yīng)的大幅晃動(dòng)(Large amplitude slosh, LAMPS)Fortran程序，通過大量二維情形的低重力晃動(dòng)試驗(yàn)驗(yàn)證了LAMPS程序在預(yù)測(cè)液體晃動(dòng)力上的準(zhǔn)確性。另一方面，文獻(xiàn)[6-7]將二維質(zhì)心面模型發(fā)展到三維情形，考慮單個(gè)貯腔作水平橫向簡(jiǎn)諧運(yùn)動(dòng)時(shí)，通過對(duì)比由Fluent軟件模擬得到的液體晃動(dòng)力驗(yàn)證了質(zhì)心面模型用于航天器貯腔大幅液體晃動(dòng)分析的有效性；隨后，采用牛頓-歐拉法將質(zhì)心面模型整合到充液航天器耦合姿態(tài)動(dòng)力學(xué)模型中，編制相應(yīng)的Fortran程序?qū)钨A腔航天器姿態(tài)機(jī)動(dòng)問題進(jìn)行了仿真計(jì)算，并與歐洲航天局2005年發(fā)射的液體晃動(dòng)實(shí)驗(yàn)衛(wèi)星(單貯腔衛(wèi)星)的在軌測(cè)量結(jié)果進(jìn)行了對(duì)比，驗(yàn)證了質(zhì)心面模型應(yīng)用于大幅晃動(dòng)類充液航天器耦合動(dòng)力學(xué)建模的可靠性。此外，文獻(xiàn)[8] 進(jìn)一步改進(jìn)和發(fā)展了質(zhì)心面模型。盡管如此，等效液體晃動(dòng)質(zhì)心面模型在等效液體黏性耗散力以及液體運(yùn)動(dòng)約束等方面尚有很大的研究空間[5]。

本文結(jié)合中國(guó)某型號(hào)航天器預(yù)研要求，基于已有的質(zhì)心面等效模型，采用拉格朗日方法，系統(tǒng)地建立了任意復(fù)雜激勵(lì)下四貯腔航天器剛-液耦合動(dòng)力學(xué)精確模型，盡管在準(zhǔn)確度上可能不及CFD-剛體耦合動(dòng)力學(xué)模型[9]，但是由于該模型所需求解的自由度遠(yuǎn)遠(yuǎn)低于CFD-剛體耦合動(dòng)力學(xué)模型，能夠?qū)崿F(xiàn)與航天器控制系統(tǒng)進(jìn)行快速高效的實(shí)時(shí)數(shù)據(jù)交互。然后，編制了相應(yīng)的MATLAB仿真計(jì)算程序，考慮失重環(huán)境下航天器受復(fù)雜激勵(lì)情形，通過對(duì)比Flow3D求解所得的液體推進(jìn)劑產(chǎn)生的晃動(dòng)力及晃動(dòng)力矩結(jié)果，驗(yàn)證了所建模型的可靠性；同時(shí)，得到了特定工況下質(zhì)心面模型的等效摩擦特性系數(shù)隨貯腔充液比的經(jīng)驗(yàn)關(guān)系。最后，考慮此類航天器在進(jìn)行三軸穩(wěn)定姿態(tài)機(jī)動(dòng)時(shí)，液體晃動(dòng)將可能對(duì)航天器的姿態(tài)運(yùn)動(dòng)帶來強(qiáng)烈的擾動(dòng)。為此，基于本文建立的耦合動(dòng)力學(xué)模型，初步設(shè)計(jì)了相應(yīng)的姿態(tài)機(jī)動(dòng)補(bǔ)償控制器，結(jié)果顯示該補(bǔ)償控制器能很好地抑制航天器姿態(tài)角速度的振蕩。

1 四充液儲(chǔ)箱航天器剛-液耦合動(dòng)力學(xué)建模

1.1 質(zhì)心面模型與儲(chǔ)箱布局介紹

1975年，Berry和Tegart[5]提出了一種適用于求解微重力環(huán)境下液體晃動(dòng)力的大范圍運(yùn)動(dòng)及大幅晃動(dòng)模型，該模型將貯腔內(nèi)的所有液體等效為點(diǎn)質(zhì)量并沿特定的約束表面運(yùn)動(dòng)，因此也稱為質(zhì)心面模型。如圖1所示，針對(duì)航天器工程中主要采用的兩類液體貯腔：卡西尼腔和球腔，質(zhì)心約束表面(質(zhì)心面)為一旋轉(zhuǎn)橢球面，在腔體坐標(biāo)系下質(zhì)心面的曲面方程有如下表示：

(1)

特別地，對(duì)于球腔，有a=b。

圖1 質(zhì)心面模型示意圖Fig.1 Illustration of the constraint surface model forliquid sloshing

圖2 航天器四貯腔布局與坐標(biāo)系Fig.2 Four tanks layout in spacecraft and frames

1.2 系統(tǒng)運(yùn)動(dòng)描述

1)航天器相對(duì)于慣性參考系的轉(zhuǎn)動(dòng)

采用航天器整體坐標(biāo)系相對(duì)于慣性參考系發(fā)生的轉(zhuǎn)動(dòng)來描述航天器的姿態(tài)運(yùn)動(dòng)，并使用無(wú)奇異的歐拉四元數(shù)描述法來定義航天器姿態(tài)矩陣

(2)

設(shè)航天器的轉(zhuǎn)動(dòng)角速度矢量在整體系中的坐標(biāo)矩陣為

并由此規(guī)定任意矢量在整體系中的坐標(biāo)表示方法。

由四元數(shù)表示的航天器姿態(tài)運(yùn)動(dòng)學(xué)方程則由如下關(guān)系式定義：

(3)

此外，重力加速度矢量在整體系下的坐標(biāo)表示為：

(4)

2)航天器相對(duì)于慣性參考系的平動(dòng)

航天器主剛體質(zhì)心在慣性參考系下的位置和速度所對(duì)應(yīng)的坐標(biāo)矩陣分別記為

(5)

(6)

于是主剛體質(zhì)心的絕對(duì)速度矢量在整體系下的坐標(biāo)表示為

(7)

3)液體等效質(zhì)心點(diǎn)的運(yùn)動(dòng)

記第i個(gè)液體質(zhì)心點(diǎn)相對(duì)于慣性系的位置矢量在整體系下的坐標(biāo)表示為

1ri=1rOi/N+1ρi

(8)

式(8)對(duì)時(shí)間求絕對(duì)導(dǎo)數(shù)，可得質(zhì)心點(diǎn)的絕對(duì)速度在整體系中的坐標(biāo)表示：

記Ri=1rO/C+1rOi/O+1ρi，質(zhì)心點(diǎn)的絕對(duì)速度在整體系中的坐標(biāo)矩陣最后可表示為

(9)

1.3 剛-液耦合系統(tǒng)動(dòng)力學(xué)方程

本文采用混合坐標(biāo)意義下的Lagrange方程[10-11]建立四充液貯腔航天器耦合動(dòng)力學(xué)模型。建模過程需要計(jì)算航天器主剛體、儲(chǔ)腔內(nèi)液體及動(dòng)量輪組成的耦合系統(tǒng)的動(dòng)能和勢(shì)能。設(shè)航天器主剛體的質(zhì)量為mhub，關(guān)于其質(zhì)心坐標(biāo)系的慣性張量為Ihub；動(dòng)量輪的慣性張量為Iw。

系統(tǒng)的動(dòng)能由如下定義式給出

(10)

式中：mi是第i個(gè)質(zhì)心點(diǎn)的質(zhì)量，數(shù)值上等于對(duì)應(yīng)貯腔內(nèi)液體總質(zhì)量；Ω是動(dòng)量輪相對(duì)整體系的轉(zhuǎn)動(dòng)速度。

將式(9)代入式(10)并整理成矩陣形式：

(11)

只考慮系統(tǒng)的重力勢(shì)能

(12)

將式(4)、(5)、(8)代入式(12)得到系統(tǒng)勢(shì)能的矩陣形式：

(13)

于是，該耦合系統(tǒng)的Lagrange函數(shù)有如下表述：

(14)

1)航天器主剛體平動(dòng)的動(dòng)力學(xué)方程

以航天器主剛體質(zhì)心的位置坐標(biāo)和速度坐標(biāo)表示的Lagrange方程為

(15)

式中：F是航天器系統(tǒng)受到的外力的合力在慣性參考系中表示的坐標(biāo)矢量，其作用點(diǎn)位于主剛體質(zhì)心。當(dāng)不計(jì)入航天器主剛體及貯腔的質(zhì)量時(shí)，由牛頓第二定律，F(xiàn)即表示四個(gè)貯腔內(nèi)液體對(duì)航天器主剛體總的等效作用反力。

把式(11)、(13)同時(shí)代入式(14)后再代入式(15)，最終得到描述航天器沿軌道平動(dòng)的動(dòng)力學(xué)方程：

(16)

2)航天器主剛體轉(zhuǎn)動(dòng)的動(dòng)力學(xué)方程

偽速度(ω)下的廣義Lagrange方程為

(17)

式中：M是航天器系統(tǒng)受到關(guān)于主剛體質(zhì)心的外部作用力矩的合力矩在整體系中表示的坐標(biāo)矢量。當(dāng)不計(jì)入航天器主剛體及貯腔的慣量并不考慮動(dòng)量輪作用時(shí)，由力矩平衡關(guān)系，M即表示四個(gè)貯腔內(nèi)液體對(duì)航天器主剛體總的等效作用反力矩。

最終得到描述航天器姿態(tài)轉(zhuǎn)動(dòng)的動(dòng)力學(xué)方程：

(1ω)×[Imb(1ω)+Iw(1ω+1Ω)+

(18)

3)質(zhì)心點(diǎn)運(yùn)動(dòng)的動(dòng)力學(xué)方程

以第i個(gè)質(zhì)心點(diǎn)的位置坐標(biāo)和速度坐標(biāo)表示的Lagrange方程為

(19)

最終得到的描述質(zhì)心點(diǎn)運(yùn)動(dòng)的控制方程為：

(20)

此外，動(dòng)量輪的輸出方程有如下表述：

(21)

1.4 求解質(zhì)心面約束力

對(duì)于第i個(gè)貯腔，質(zhì)心面在腔體坐標(biāo)系下的曲面方程如式(1)表示，當(dāng)質(zhì)心點(diǎn)被約束在質(zhì)心面上運(yùn)動(dòng)(聯(lián)系運(yùn)動(dòng))時(shí)，其沿曲面外法向ni的速度投影應(yīng)始終為0，即有

(22)

并且有

(23)

式(22)對(duì)時(shí)間求微分可得聯(lián)系運(yùn)動(dòng)的約束條件為

(24)

將式(23)代入式(24)，得

(25)

再將式(20)代入式(25)，可得質(zhì)心點(diǎn)受到的法向支持力大小為

(26)

2 復(fù)雜激勵(lì)下等效晃動(dòng)力與等效晃動(dòng)力矩驗(yàn)證

為了驗(yàn)證基于等效液體晃動(dòng)質(zhì)心面模型，根據(jù)我國(guó)某型號(hào)航天器預(yù)研的需要建立四充液貯腔航天器耦合動(dòng)力學(xué)模型，這里，在航天器主剛體質(zhì)心處施加全自由度的加速度激勵(lì)，將液體所產(chǎn)生的作用于主剛體質(zhì)心的晃動(dòng)力及晃動(dòng)力矩的模型計(jì)算結(jié)果與Flow3D計(jì)算結(jié)果進(jìn)行對(duì)比。其中Flow3D計(jì)算方案主要包括：推進(jìn)劑處理為不可壓縮、黏性流體；使用單相流模型，氣體部分僅提供環(huán)境壓力，其密度相對(duì)液體很小，不進(jìn)行動(dòng)量計(jì)算；考慮表面張力；液固界面采用Navier滑移邊界條件；先計(jì)算得到表面張力作用下穩(wěn)態(tài)時(shí)的液面構(gòu)型，再以此為初始條件導(dǎo)入給定激勵(lì)進(jìn)行計(jì)算。

給定主剛體質(zhì)心同時(shí)受到1000 s的三維平動(dòng)加速度激勵(lì)與三維角加速度激勵(lì)，它們的時(shí)程圖分別如圖3和圖4所示。

圖3 三維平動(dòng)加速度激勵(lì)時(shí)程Fig.3 Acceleration excitation in three directions

圖4 三維角加速度激勵(lì)時(shí)程Fig.4 Angular acceleration excitation in three directions

表1 推進(jìn)劑特性參數(shù)Table 1 Characteristic parameters of propellant

表2 充液儲(chǔ)腔球心位置坐標(biāo)Table 2 Coordinates of the center of the tanks

本文給出了以下兩種工況下的晃動(dòng)力和晃動(dòng)力矩的對(duì)比結(jié)果：工況1為各推進(jìn)劑貯腔充液比均為20%(圖5、圖6)；工況2為各推進(jìn)劑貯腔充液比均為60%(圖7、圖8)。計(jì)算結(jié)果表明，充液比較高時(shí)，模型的計(jì)算結(jié)果更為準(zhǔn)確，主要原因是由于當(dāng)貯腔內(nèi)液體越少時(shí)，同等水平激勵(lì)下晃動(dòng)的幅度越大，并且Flow3D計(jì)算顯示液體出現(xiàn)更為明顯的破碎。總的來說，晃動(dòng)力和晃動(dòng)力矩的對(duì)比結(jié)果體現(xiàn)出較高的吻合度，說明本文建立的四充液貯腔剛-液耦合動(dòng)力學(xué)模型具有應(yīng)用于航天工程實(shí)際的重大價(jià)值。

圖5 晃動(dòng)力:質(zhì)心面模型vs.Flow3D(工況1)Fig.5 Slosh force: model results vs. Flow3Dresults(Case 1)

圖6 晃動(dòng)力矩:質(zhì)心面模型vs.Flow3D(工況1)Fig.6 Slosh torque: model results vs. Flow3D results(Case 1)

圖7 晃動(dòng)力:質(zhì)心面模型vs.Flow3D(工況2)Fig.7 Slosh force: model results vs. Flow3Dresults(Case 2)

圖8 晃動(dòng)力矩:質(zhì)心面模型vs.Flow3D(工況2)Fig.8 Slosh torque: model results vs. Flow3D results(Case 2)

圖9 等效摩擦特性系數(shù)與充液比經(jīng)驗(yàn)關(guān)系Fig.9 Empirical relationship between μi andliquid-filled ratio

通過更多充液比工況下的對(duì)比研究，對(duì)于這里所給的球腔尺寸、液體黏性特性以及加速度環(huán)境，本文得到了質(zhì)心面模型中的等效摩擦特性系數(shù)μi與充液比的近似經(jīng)驗(yàn)關(guān)系如圖9所示。

此外，本文基于質(zhì)心面等效模型建立的多貯腔航天器剛-液耦合動(dòng)力學(xué)模型在評(píng)估液體晃動(dòng)力和晃動(dòng)力矩時(shí)的計(jì)算效率遠(yuǎn)高于Flow3D計(jì)算流體動(dòng)力學(xué)仿真軟件。使用3.4 GHz的PC機(jī)，計(jì)算時(shí)長(zhǎng)為1000 s，編制MATLAB四階-五階Runge-Kutta(ode 45)算法解算本文所建立的動(dòng)力學(xué)模型只需2 h左右，而使用Flow3D仿真軟件則需3天甚至更多時(shí)間(主要取決于計(jì)算網(wǎng)格尺寸相關(guān))。

3 三軸穩(wěn)定姿態(tài)機(jī)動(dòng)設(shè)計(jì)

本節(jié)中，不考慮航天器軌道運(yùn)動(dòng)，啟動(dòng)動(dòng)量輪對(duì)航天器進(jìn)行文獻(xiàn)[4]中的“rest-to-rest”形式(初始角速度與末態(tài)角速度均為零)的姿態(tài)機(jī)動(dòng)，并研究機(jī)動(dòng)過程中的剛-液耦合動(dòng)力學(xué)特性。假定各個(gè)貯腔的充液比均為60%；航天器主剛體不受到任何外部作用力矩，即M=0。

研究發(fā)現(xiàn)，當(dāng)航天器處于完全失重環(huán)境，且不受到任何強(qiáng)干擾(或激勵(lì))作用時(shí)，對(duì)航天器進(jìn)行三軸穩(wěn)定姿態(tài)機(jī)動(dòng)過程中，主體姿態(tài)運(yùn)動(dòng)與液體晃動(dòng)之間的動(dòng)力學(xué)耦合效應(yīng)不顯著。進(jìn)一步考慮在軌航天器的如下工況，即航天器主發(fā)動(dòng)機(jī)點(diǎn)火(認(rèn)為是一種強(qiáng)干擾)產(chǎn)生一個(gè)持續(xù)的推力使得航天器獲得一個(gè)恒定的加速度，即航天器將處于一個(gè)等效的重力加速度環(huán)境，假定g=1 m/s2。航天器主剛體關(guān)于其質(zhì)心坐標(biāo)系的慣性張量為Ihub=diag(2700,1300,2300) kg·m2。

圖10為采用文獻(xiàn)[4]中的線性PD控制器：

Tw=kpεv+kd(1ω)

取kp=250,kd=2000，進(jìn)行三軸穩(wěn)定姿態(tài)機(jī)動(dòng)時(shí)的航天器主體角速度響應(yīng)。結(jié)果顯示姿態(tài)機(jī)動(dòng)過程中，貯腔內(nèi)液體晃動(dòng)將對(duì)航天器的姿態(tài)運(yùn)動(dòng)帶來強(qiáng)烈的擾動(dòng)，這種影響的主要原因是推進(jìn)劑大幅晃動(dòng)將造成航天器整體的慣性特性顯著變化，同時(shí)產(chǎn)生較大的晃動(dòng)擾動(dòng)力和擾動(dòng)力矩。這有可能引起航天器的結(jié)構(gòu)破壞，甚至造成航天器姿態(tài)翻滾等嚴(yán)重后果，在工程實(shí)際中應(yīng)設(shè)法解決或避免。

為此，初步設(shè)計(jì)了相應(yīng)的姿態(tài)機(jī)動(dòng)補(bǔ)償控制器：

同樣取kp=250,kd=2000，圖11為機(jī)動(dòng)過程中相應(yīng)的角速度響應(yīng)。結(jié)果表明，該類控制器能很好地抑制航天器姿態(tài)角速度的振蕩，從而能有效避免實(shí)際工程中的潛在風(fēng)險(xiǎn)。

圖10 PD控制下三軸穩(wěn)定姿態(tài)機(jī)動(dòng)角速度響應(yīng)Fig.10 Angular velocity response under PD controller

圖11 補(bǔ)償控制下三軸穩(wěn)定姿態(tài)機(jī)動(dòng)角速度響應(yīng)Fig.11 Angular velocity under compensation controller

4 結(jié) 論

本文基于已有的質(zhì)心面液體晃動(dòng)等效力學(xué)模型，采用拉格朗日方法，系統(tǒng)地建立了計(jì)及航天器質(zhì)心平動(dòng)、姿態(tài)轉(zhuǎn)動(dòng)與多貯腔內(nèi)液體推進(jìn)劑大幅晃動(dòng)的剛-液耦合動(dòng)力學(xué)精確模型。通過對(duì)比Flow3D求解所得的完全失重、復(fù)雜激勵(lì)下液體產(chǎn)生的晃動(dòng)力及晃動(dòng)力矩結(jié)果，驗(yàn)證了所建模型的可靠性；同時(shí)，得到了特定工況下質(zhì)心面模型等效摩擦特性系數(shù)與充液比之間的經(jīng)驗(yàn)關(guān)系，具體呈現(xiàn)為二次函數(shù)的近似關(guān)系。最后，基于本文中建立的耦合動(dòng)力學(xué)模型，研究了該類航天器在進(jìn)行三軸穩(wěn)定姿態(tài)機(jī)動(dòng)時(shí)的推進(jìn)劑晃動(dòng)與航天器主體姿態(tài)運(yùn)動(dòng)之間的動(dòng)力學(xué)耦合特性，主要表現(xiàn)為推進(jìn)劑大幅晃動(dòng)造成航天器整體的慣性特性顯著變化，同時(shí)產(chǎn)生較大的晃動(dòng)擾動(dòng)力和擾動(dòng)力矩，從而給航天器姿態(tài)運(yùn)動(dòng)帶來明顯的振蕩影響；為此，本文設(shè)計(jì)了一類抑制姿態(tài)振蕩的補(bǔ)償控制器，并且計(jì)算結(jié)果驗(yàn)證了該控制器的有效性。