劉曉丹， 虞育松*， 李 永， 王淑一， 梁軍強， 宋 濤

1. 北京交通大學，機電學院氫能與航天推進技術(shù)實驗室， 100044 北京 2. 北京控制工程研究所， 100094 北京

0 引 言

在微重力環(huán)境下，充液航天器內(nèi)液體晃動及其對控制系統(tǒng)的影響是當前國內(nèi)外航天技術(shù)研究的重要課題.貯箱液體晃動對航天器的結(jié)構(gòu)強度和運動狀態(tài)都具有重要影響[1].特別是在變軌、空間對接等情況下，液體晃動明顯偏離線性區(qū)[2]，流動的非線性對于星體控制和晃動抑制技術(shù)都提出新的挑戰(zhàn).此外，隨著深空探測和載人任務(wù)的發(fā)展，推進劑占火箭和衛(wèi)星等航天器總重量的比重不斷加大.如2004年NASA發(fā)射的信使號探測器的推進劑液體質(zhì)量占航天器總質(zhì)量的56%.于是，推進劑晃動對航天器姿態(tài)運動的影響增強，通過慣性力相互作用發(fā)生的耦合現(xiàn)象更加顯著.在航天器充液量不斷增大和姿態(tài)控制精度指標不斷提高的前提下，深入開展空間環(huán)境下運動航天器內(nèi)貯箱液體晃動動力學的研究將具有重要意義[3].另外，貯箱液體晃動也是航天器總體設(shè)計、撓性航天器姿態(tài)穩(wěn)定控制問題中需要考慮的基礎(chǔ)問題.

由于受到實驗成本、數(shù)據(jù)測量等客觀條件限制，目前很難采用實驗方法進行在軌航天器貯箱液體晃動動力學行為研究.當前普遍采用數(shù)值仿真或理論方法來研究上述問題.數(shù)值方法一般采用連續(xù)流模擬方法對特定貯箱運動狀態(tài)下的液體晃動特性進行計算和分析[4-6].貯箱運動狀態(tài)均通過運動邊界條件的方式來設(shè)定，并不考慮星體動力學模塊對貯箱受力的反饋控制作用，因此與實際在軌情況存在一定差異.

在與星體動力學狀態(tài)的耦合控制方面，目前一般都采用理論方法來對貯箱晃動特性進行建模，然后與航天器控制系統(tǒng)進行耦合.基于液體小幅晃動假設(shè)的等效晃動力學模型是早期的理論模型.隨后由于需要對大幅液體晃動條件(交會對接、大角度姿態(tài)機動等)進行分析，人們又陸續(xù)發(fā)展了大幅晃動等效力學模型，如復(fù)合擺等效力學模型[7]、MPBM模型[8]、質(zhì)心面模型[9]等.最近，針對航天器剛-液耦合問題又發(fā)展了運動脈動球模型[10]，可對航天器位置-姿態(tài)跟蹤控制進行理論分析；針對外部激勵作用又提出了大幅晃動理論模型[11]，現(xiàn)實了對橫向激勵作用的分析.

綜上所述，目前在衛(wèi)星閉環(huán)姿態(tài)控制數(shù)學仿真中，一般采用等效力學理論模型對液體晃動特性進行求解[12]，較難精確模擬大幅晃動情況，也很難考慮復(fù)雜工況和貯箱內(nèi)部葉片、隔板、導(dǎo)流板等結(jié)構(gòu)對液體晃動的影響[13].為了能夠更加準確的考慮微重力下液體晃動對姿態(tài)控制性能的影響，以及貯箱結(jié)構(gòu)和星體運動對液體晃動的影響，本文將三維CFD貯箱液體晃動模型與星體動力學控制模型進行耦合，通過通訊接口開發(fā)實現(xiàn)衛(wèi)星運動與貯箱液體晃動過程的聯(lián)合仿真，并在此基礎(chǔ)上分析典型星體運動工況下的貯箱液體晃動特性，分析貯箱葉片結(jié)構(gòu)對液體晃動的影響.

1 聯(lián)合仿真方法和數(shù)學模型

1.1 星體動力學與貯箱液體晃動聯(lián)合仿真

本次研究的貯箱液體晃動聯(lián)合仿真平臺是由CFD軟件FLUENT與星體動力學控制程序構(gòu)成.FLUENT計算運動貯箱的液體晃動力、力矩、質(zhì)心結(jié)果，星體動力學程序在獲取貯箱晃動結(jié)果后通過控制算法和策略計算星體運動的線速度和角速度.兩者之間的雙向數(shù)據(jù)交換通過C/C++語言開發(fā)的接口來實現(xiàn).星體動力學模塊對貯箱受力的反饋控制作用通過接口程序傳遞給各貯箱晃動(FLUENT計算)的線速度和角速度功能來實現(xiàn).本次研究開展的聯(lián)合仿真數(shù)據(jù)交換量較小，更注重數(shù)據(jù)交互的實時性，因此采用了UDP協(xié)議作為Server端星體動力學計算和Client端液體晃動仿真數(shù)據(jù)的手段.FLUENT采用UDF基于C編譯相應(yīng)的宏來實現(xiàn)與星體動力學程序的數(shù)據(jù)交換、貯箱液體晃動并行計算.星體動力學程序與FLUENT軟件的聯(lián)合仿真平臺組成見圖1.

圖1 星體動力學程序與FLUENT軟件的數(shù)據(jù)交換流程Fig.1 Data exchange process between spacecraft dynamics program and FLUENT software

星體動力學程序是以動態(tài)鏈接庫形式進行集成，需要通過C/C++建立Server端程序?qū)⑿求w動力學動態(tài)鏈接庫導(dǎo)入；同時為星體動力學仿真定義各種初始條件，并按照指定的時間間隔導(dǎo)入貯箱晃動仿真數(shù)據(jù).與其它通信協(xié)議相比，UDP協(xié)議還具有程序架構(gòu)簡單、理論上可支持任意數(shù)量Client掛載Server等的優(yōu)點.圖2給出了星體動力學聯(lián)合仿真框架.

圖2 星體動力學聯(lián)合仿真框架Fig.2 Joint simulation framework for spacecraft dynamics

星體動力學程序與貯箱晃動CFD算例之間基于UDP通訊框架如圖3所示.

1.2 貯箱液體晃動數(shù)學模型

(1)氣液兩相流控制方程

不可壓、不相溶、等溫的氣液兩相流Navier-Stokes方程可表示為

(1)

(2)

(3)

式(1)為連續(xù)性方程，式(2)為動量守恒方程.式(2)考慮了壓力、湍流、表面張力、重力對兩相流動的影響.u為流體速度，ρ為密度，F(xiàn)和g分別為表面張力和重力，μ為粘度.流動粘度和密度可根據(jù)氣液體積分數(shù)值進行加強平均獲得.

(2)VOF方法

不同流體分子間作用力的差異將導(dǎo)致相界面上出現(xiàn)表面張力，利用體積分數(shù)方法(volume of fluid，VOF)對氣液兩相及相界面進行描述[14].引入液相體積分數(shù)ε：

(4)

根據(jù)體積分數(shù)ε可將流場分為3個區(qū)域：氣相區(qū)(ε=0)、液相區(qū)(ε=1)、氣液界面區(qū)(0<ε<1).氣液界面的流體屬性，包括流體粘度、密度采用體積分數(shù)的線性加權(quán)平均來進行計算.速度壓力耦合采用SIMPLE算法，使用PRESTO!壓力插值方案，動量方程離散采用二階迎風格式方案.為提高相間界面的清晰度，流體體積分數(shù)求解采用second-order離散方法.

2 貯箱液體晃動聯(lián)合仿真分析

2.1 PMD貯箱建模與參數(shù)設(shè)置

衛(wèi)星推進分系統(tǒng)采用板式全管理表面張力貯箱.星體4貯箱布置結(jié)構(gòu)見圖4.MMH和MON-1推進劑為對角布置.假定1、4貯箱推進劑為MMH，2和3貯箱推進劑為MON-1.

該貯箱主要由貯箱殼體和板式推進劑管理裝置(板式PMD)組成.板式PMD是該貯箱的核心部件.為了簡化結(jié)構(gòu)，降低CFD計算耗時，本次研究將導(dǎo)流板、底部蓄液葉片和蓄液器等結(jié)構(gòu)忽略，只保留4個大葉片，以分析葉片對液體晃動的影響.為了分析貯箱葉片小孔結(jié)構(gòu)對液體晃動的影響，本次研究還對比了有無孔葉片結(jié)構(gòu)下的晃動特性.帶孔和無孔葉片的貯箱計算域及網(wǎng)格劃分如圖5所示.

圖4 星體4貯箱布置Fig.4 Arrangement of 4 tanks in a spacecraft

貯箱內(nèi)流體溫度為20℃，零重力條件.貯箱內(nèi)推進劑分別為MMH和MON-1.其中，MMH的物性參數(shù)(20℃)：密度874 kg/m3、表面張力系數(shù)0.033 95 N/m、動力粘度0.9×10-3Pa·s；MON-1的物性參數(shù)(20℃)：密度1 444 kg/m3、表面張力系數(shù)0.025 6 N/m、動力粘度0.9×10-3Pa·s.本次計算設(shè)置的壁面接觸角為1°.貯箱液體晃動仿真的初始液面為重定位過程的終了平衡狀態(tài)，貯箱充液比為40%.

流體運動速度在貯箱壁面上滿足無滑移條件，等溫壁面邊界且不考慮傳熱的影響.計算時的全局Courant數(shù)限制在0.6以下.為與星體動力學計算程序保持同步計算，F(xiàn)LUENT液體晃動計算過程采用定時間步長，時間步長設(shè)定為0.000 2 s.Server端星體動力學程序設(shè)定的計算時間步長為0.001 s.因此，在具體計算時，Client端的液體晃動仿真每間隔5步，向Server送晃動結(jié)果數(shù)據(jù)；Server根據(jù)接收的晃動數(shù)據(jù)計算貯箱線速度、角速度、線加速度和角加速度值，并回發(fā)給對應(yīng)的貯箱算例.

2.2 貯箱晃動聯(lián)合仿真模型的驗證

圖6對比了俯仰0～25°機動角速度的聯(lián)合仿真與遙感數(shù)據(jù)對比.結(jié)果顯示，偏航角速度峰值的仿真偏差為-7.31%，俯仰角速度峰值的仿真偏差為-5.58%.俯仰角速度和偏航角速度的仿真結(jié)果與在軌測試結(jié)果變化規(guī)律和趨勢均非常接近，顯示了聯(lián)合仿真計算結(jié)果具有較高的精度.

2.3 貯箱晃動聯(lián)合仿真模型的驗證

本次計算將貯箱晃動的起始點設(shè)置在非沿軌跡工況的起始點，整個控制過程的時間約為80 s.該工況以X方向的角速度為主，峰值角速度約4.5°/s.圖7和圖8給出了該工況下4個貯箱的液體晃動力、力矩(以各貯箱幾何中心點為坐標系原點)和液面變動過程.星體轉(zhuǎn)動過程中4個貯箱內(nèi)的晃動力和力矩發(fā)展過程比較相似，但存在量值方面的差異.在約5 s和18 s時的Y和Z方向的液體晃動力存在峰值，2、3號貯箱的晃動力和力矩峰值大于1、4號貯箱，4個貯箱的Y向晃動力峰值在1.05～1.96 N，力矩值在0.05～0.36 N·m.

圖6 俯仰0～25°機動角速度仿真與遙感數(shù)據(jù)對比Fig.6 Comparison of simulation and remote sensing data of the pitch 0～25° maneuvering angular velocity

圖7 非沿軌跡工況下的四個貯箱液體晃動力和力矩Fig.7 Liquid slosh force and moment for four tanks under non-tracking conditions

圖8 帶孔葉片下的4貯箱液體晃動過程Fig.8 The liquid slosh processes of four tanks with perforated blades

2.4 葉片結(jié)構(gòu)對液體晃動特性的影響

圖9、圖10和圖11分別給出了貯箱內(nèi)布置帶孔葉片、無孔葉片情況下的液體質(zhì)心、晃動力和力矩結(jié)果.非沿軌跡工況下的液體晃動質(zhì)心變動均主要體現(xiàn)在X方向，不帶孔葉片的X方向液體質(zhì)心晃動幅值要小于帶孔葉片情況，但無孔葉片會導(dǎo)致Y和Z方向液體質(zhì)心變動幅度增加.這是由于貯箱轉(zhuǎn)動時，不帶孔葉片對液體的攪動程度要大于帶孔葉片.

圖9 葉片結(jié)構(gòu)對液體質(zhì)心變動的影響Fig.9 Effects of blade structure on the variation of the liquid mass center

圖10給出了20 s時帶孔葉片對貯箱內(nèi)液體的作用，可以看到葉片運動時，液體可從小孔中流過，降低了葉片運動對液體的攪動.

圖10 葉片附近氣液界面及速度分布(20 s)Fig.10 Liquid-gas surface and velocity distributionnear the blades

根據(jù)帶孔和無孔葉片條件下的液體晃動力、力矩結(jié)果對比可知(圖11和12)，非沿軌跡工況下，無孔葉片時的液體晃動力和力矩峰值顯著大于帶孔葉片情況.帶孔葉片的晃動力峰值出現(xiàn)在76 s，峰值為-0.122 N(X方向).而無孔葉片的晃動力峰值出現(xiàn)略微提前至72 s，峰值為-0.8 N(X方向).另外，葉片不開孔時，晃動力矩在晃動后期存在較大波動，力矩幅值變化較為劇烈.

圖11 葉片結(jié)構(gòu)對液體晃動力的影響Fig.11 Effects of blade structure on the liquid slosh force

圖12 葉片結(jié)構(gòu)對液體晃動力矩的影響Fig.12 Effects of blade structure on the liquid slosh moment

圖13和14給出了晃動過程中的液面形態(tài)和液面速度分布情況.晃動開始后，由于受到X方向的轉(zhuǎn)動作用，葉片轉(zhuǎn)動后帶動氣液體運動.隨著轉(zhuǎn)動速度的增加，氣泡逐漸在X方向被拉長.當星體運動至40 s時，X方向角速度開始減小，葉片轉(zhuǎn)速降低.隨后貯箱內(nèi)的液體在慣性力作用下繼續(xù)運動，與葉片和貯箱壁面發(fā)生較為激烈的作用，導(dǎo)致60 s之后的力和力矩變動較大.帶孔葉片對液體的攪動作用明顯要弱于無孔葉片.無孔葉片攪動的氣泡出現(xiàn)了明顯拉伸與離心運動(圖13)，并導(dǎo)致更大的流體速度(圖14).

圖13 葉片結(jié)構(gòu)對貯箱內(nèi)液面分布的影響Fig.13 Effects of blade structure on the liquid surface distribution in the tank

圖14 葉片結(jié)構(gòu)對貯箱內(nèi)液面晃動速度的影響Fig.14 Effects of blade structure on the liquid slosh velocity in the tank

3 結(jié) 論

本文建立了基于UDP通訊協(xié)議的貯箱液體晃動聯(lián)合仿真平臺，并開展了星體控制下的貯箱液體晃動特性研究，主要結(jié)論如下：

(1)通過對比俯仰0～25°機動角速度的聯(lián)合仿真與遙感數(shù)據(jù)，聯(lián)合仿真計算獲得的偏航角速度峰值偏差為-7.31%，俯仰角速度峰值偏差為-5.58%.俯仰角速度和偏航角速度的仿真結(jié)果與在軌測試結(jié)果變化規(guī)律和趨勢均非常接近，說明聯(lián)合仿真計算結(jié)果具有較高的精度.

(2)非沿軌跡工況下，星體轉(zhuǎn)動過程中四個貯箱內(nèi)的晃動力和力矩發(fā)展過程相似，但存在量值方面的差異.在約5 s和18 s時的Y和Z方向液體晃動力存在峰值，2、3號貯箱的晃動力和力矩峰值大于1、4號貯箱.4個貯箱的Y向晃動力峰值在1.05～1.96 N，力矩值在0.05～0.36 N·m.

(3)非沿軌跡工況下的液體晃動質(zhì)心變動均主要體現(xiàn)在X方向，不帶孔葉片的X方向液體質(zhì)心晃動幅值要小于帶孔葉片情況.葉片不開孔時，晃動力矩在晃動后期存在較大波動.帶孔葉片對液體的攪動作用明顯要弱于無孔葉片.無孔葉片對液體的攪動作用最強，導(dǎo)致了顯著的氣泡拉伸和流體流動.