關(guān) 宏 徐世杰

北京航空航天大學宇航學院，北京 100191

?

基于Fluent計算的火箭離軌姿態(tài)運動仿真與分析

關(guān) 宏 徐世杰

北京航空航天大學宇航學院，北京 100191

在航天器姿態(tài)運動分析中常遇到液體的流動、晃動及流/固耦合等問題，針對這類問題建立精確的數(shù)學模型是很困難的，通常需要利用數(shù)值計算和仿真手段進行分析。利用用戶互動功能(UDF)，引入自定義的變量源函數(shù)，應(yīng)用Fluent對運載火箭離軌姿態(tài)控制進行剩余燃料的流場數(shù)據(jù)分析。通過在流體運動方程中加入相關(guān)的牽連運動，得到在軌運行狀態(tài)下貯箱內(nèi)剩余液體的運動參數(shù)和對貯箱的干擾力矩，為運載火箭系統(tǒng)的姿態(tài)運動分析和仿真提供運算參數(shù)，使得流場的變化與運載火箭的姿態(tài)運動相關(guān)聯(lián)，分析各個時刻流場運動狀態(tài)對箭體姿態(tài)的影響。

Fluent;離軌;剩余液體;用戶互動功能(UDF)

隨著運載火箭技術(shù)的進步，使用液體推進劑的火箭已經(jīng)成為航天器的主要運載工具。一些一次性使用的運載火箭和重復使用的航天飛行器，在主發(fā)動機關(guān)機后，不可用推進劑、安全儲備、以及飛行混合比偏差所引起的剩余推進劑都將存留在貯箱里。這些殘存的液體燃料在貯箱里的運動將會影響到運載火箭的姿態(tài)運動，甚至引起火箭姿態(tài)失穩(wěn)。因此對貯箱里殘存液體的運動分析和計算是十分必要的。實際工程中，運載火箭在分離后的回收過程中常遇到剩余液體晃動產(chǎn)生的問題。一般是將剩余液體的運動作為對火箭的干擾來處理，根據(jù)剩余液體不同的分布規(guī)律，采用偏微分方程描述其運動變化規(guī)律[1-2]，并作為姿態(tài)運動的流體分析干擾模型。但實際上，當殘留液體體積遠小于貯箱容積時，此方法難以反映貯箱-液體系統(tǒng)參數(shù)的瞬態(tài)分布特征，如果對貯箱和液體建立完整的動態(tài)數(shù)學模型，則存在建模困難、計算量大的問題，而且往往不能反映符合實際的參數(shù)分布特征。

目前，隨著計算流體力學、分布式仿真技術(shù)的發(fā)展[3]，利用成熟的CFD仿真軟件Fluent，對工程中常遇到的具有非定常分布模型的流場能夠有效地進行計算，并且其建模簡單、計算精度和可靠性都很高，在工程預研、仿真驗證中得到了廣泛的應(yīng)用。Fluent軟件具有很多優(yōu)點，但還存在一定的局限性：一方面，僅靠流場仿真無法揭示航天器姿態(tài)運動和貯箱內(nèi)流體運動的相互作用；另一方面，需要將Fluent加載到其他系統(tǒng)仿真軟件中，若接口設(shè)計不當或網(wǎng)絡(luò)連接不穩(wěn)定將會影響仿真效率[4]。

本文根據(jù)Fluent軟件本身的特點，利用其用戶互動功能(UDF)[4]，通過系統(tǒng)參數(shù)的實時同步傳輸，使姿態(tài)動力學與運動學模塊能夠?qū)崟r地影響Fluent計算參數(shù)，從而建立了研究分布參數(shù)的閉環(huán)控制仿真平臺，為含有流體運動的姿態(tài)運動仿真提供新的途徑。通過利用這種仿真手段，本文對某運載火箭離軌前的姿態(tài)穩(wěn)定控制進行了數(shù)值計算與仿真，驗證了該方法的可用性和有效性，分析了含有剩余燃料的情況下箭體離軌姿態(tài)的變化規(guī)律。

1 仿真系統(tǒng)的建立

本文選擇某運載火箭末級推進器為研究對象，在完成軌道任務(wù)后氧化劑與燃料貯箱內(nèi)均剩余少量推進劑，體積遠小于貯箱容積。離軌過程中由于姿態(tài)調(diào)整以及箭體自旋運動，液體的分布難以確知，運動存在不可控因素，故采用Fluent流場模型對貯箱內(nèi)液體進行仿真，同時通過UDF功能實現(xiàn)箭體姿態(tài)與液體運動的閉環(huán)仿真。

此處需要說明的是當箭體與液體作為整體系統(tǒng)看待時，不受外界干擾力矩作用，但是為了分析流體運動對箭體姿態(tài)的影響，需要將箭體與液體作為兩個子系統(tǒng)，以分析其相互作用。

1.1 箭體系統(tǒng)動力學分析

首先以箭體系統(tǒng)作為研究對象，其轉(zhuǎn)動慣量沿自旋軸對稱，箭體本體坐標系Sb定義如圖1所示，坐標原點Ob在箭體質(zhì)心，Xb軸沿本體縱向?qū)ΨQ軸方向，Yb軸沿箭體縱切面向上，Zb軸按右手法則定義。

圖1 箭體本體坐標系定義

定義系統(tǒng)慣量主軸沿坐標軸方向：繞橫軸的轉(zhuǎn)動慣量為Iy，Iz，且Iy=Iz，繞縱軸(自旋軸)的轉(zhuǎn)動慣量為Ix，且Ix

(1)

其中ωx,ωy,ωz分別為繞Xb,Yb,Zb軸的角速度。

方程左側(cè)的Tx,Ty,Tz為箭體系統(tǒng)受到的外干擾力矩，本文即為剩余液體對箭體的作用力矩。

由式(1)以及軸對稱剛體運動的特性可知：

1)剛體自旋軸Xb繞角動量H的圓錐運動稱為空間章動，其空間章動速率Ω=H/It；

2)章動角γ為動量矩H與自旋軸之間的夾角，滿足關(guān)系式cosγ=Ixωx/H；

圖2 本體錐和空間錐

定義軌道參考系S系，與t=0時刻的本體系重合，x軸軌道切向，y軸沿t=0時刻箭體縱切面向上，z軸按右手法則定義。定義空間章動參考系S′系，坐標原點在航天器質(zhì)心，x′軸沿本體縱向?qū)ΨQ軸方向，y′軸沿章動角速度Ω在本體系的投影方向，t=0時刻箭體縱切面向上，隨空間錐運動，z′軸按右手法則定義。在t=0時刻，此3坐標系重合。運動Δt時刻后，y軸在空間指向不變，y′軸運動α角，Yb軸運動到β角位置，α=∫ΔtΩdt為箭體空間章動運動角度，β=∫Δtωxdt為箭體自旋運動角度。即

1.2 液體子系統(tǒng)運動

得到箭體空間運動規(guī)律，再考慮液體系統(tǒng)的相應(yīng)晃動影響。

液體子系統(tǒng)在微重力環(huán)境下運動，受擾動后不能斷定其確切位置，液氣界面穩(wěn)定性差。一旦航天器進行姿態(tài)控制(起旋至穩(wěn)定旋轉(zhuǎn))，剩余液體產(chǎn)生的擾動就會影響到姿態(tài)控制的精度，破壞航天器的運動穩(wěn)定性。又因為液體運動規(guī)律不易簡化，如不能準確仿真其運動狀態(tài)，便無法找到適當?shù)目刂撇呗?，會對回收造成惡劣的影響[6]。

本文流場模型選擇圖2中S′系為參考系，該參考系隨箭體在空間不斷運動，始終保證液體受空間錐離心加速度方向在y′軸上，網(wǎng)格坐標系原點為箭體質(zhì)心位置，便于對Fluent定義相對加速度以及對質(zhì)心求作用力矩。

基于質(zhì)量守恒定律、動量定理、能量守恒定律、熱力學定律以及流體的本身物性，在流體力學中存在一組制約流體運動的基本方程組，對于黏性不可壓縮流體，滿足：

(2)

其中，ρ為流體密度，v為流體相對運動速度，ρ=const，div(P)表示單位體積上應(yīng)力張量的散度，f1為單位質(zhì)量上的質(zhì)量力分布函數(shù)。

從動量定理出發(fā)，任取一體積τ的流體，它的邊界為A，其中流體動量的變化率等于作用于該體積上的質(zhì)量力和面力之和。以f1表示作用在單位質(zhì)量上的質(zhì)量力分布函數(shù)，pn為作用在單位面積上的面力分布函數(shù)，則作用在τ上和A上的總質(zhì)量力和面力為∫τρfδτ及∫SpnδA，而體積τ內(nèi)的動量是∫τρvδτ，于是動量定理寫成

(3)

由于∫ApnδA=∫τdiv(P)δτ，故得到式(2)中右側(cè)第2項，又ρf1表示單位體積上的質(zhì)量力。根據(jù)相對運動學原理，流體運動方程將寫成：

(4)

在原有的流體運動方程中，流場網(wǎng)格模型可以計算出式(4)中的ar，即流體相對于S′系的相對加速度。但是要考慮與箭體運動的耦合運動，就需要補充由于旋轉(zhuǎn)相對運動產(chǎn)生的慣性牽連加速度ae，以及柯氏加速度2(ω×vr)。本文采用S′系作為流場仿真的參考系，其本身就具有相對本體系旋轉(zhuǎn)的相對運動，因此只需要考慮牽連運動，即此處ae。

如前節(jié)中指出箭體上各點相對于角動量H方向的向心加速度a=Ω2r，此處的ae=a。通過對箭體的姿態(tài)運動學求解，得到各個時刻的相對姿態(tài)角、姿態(tài)角速度，根據(jù)式(1)以及軸對稱剛體運動的特性，即可通過瞬時姿態(tài)信息，解得空間章動速率Ω、章動角γ，通過三角函數(shù)關(guān)系便可輕松的得到各點到角動量H的距離r。由此得到需要嵌入到網(wǎng)格計算模型中的自定義加速度a。

1.3 Fluent-UDF接口的設(shè)計

為了解決流場模型與箭體姿態(tài)運動模型數(shù)據(jù)交換問題，在保證流場模型和姿態(tài)運動模型實時同步運行的前提下，數(shù)據(jù)準確交換，提出使用Fluent-UDF數(shù)據(jù)交換接口實現(xiàn)交互仿真。UDF與Fluent數(shù)據(jù)交互流程如圖3所示。

圖3 使用UFD的流體模型與姿態(tài)運動模型交互示意圖

首先，根據(jù)流體網(wǎng)格中的初始化邊界條件(其中包括了對流場加速度a的初始化定義)，進行1個步長的網(wǎng)格計算。得到這一步中網(wǎng)格內(nèi)剩余液體對箭體的干擾力矩Tx，Ty和Tz。下面就進入到UDF程序計算中：

第1步：由Fluent流場模型解決流體運動的分布計算問題，根據(jù)流體對箭體(貯箱)的作用情況，定義UDF程序的入口：液體子系統(tǒng)輸出對模型壁面的作用力、作用力矩；通過流場模型得到剩余液體對于箭體系統(tǒng)的干擾Tx，Ty和Tz，即為UDF程序的調(diào)用入口參數(shù)；

第2步：當UDF從入口處讀出Tx，Ty和Tz后，再根據(jù)UDF程序中編寫的姿態(tài)控制方程，計算包括了仿真對象的姿態(tài)動力學與運動學方程的姿態(tài)運動狀況；

第3步：仍然是在UDF程序包中，調(diào)用加速度計算和轉(zhuǎn)換程序。根據(jù)相對作用原理，得到流場各點的下一個時刻的牽連加速度a1，通過UDF輸出接口賦值到Fluent網(wǎng)格計算環(huán)境中。

此時，完成一個步長的箭體與液體系統(tǒng)的交互仿真計算，根據(jù)此時刻的網(wǎng)格場內(nèi)加速度a1，進行下一個步長的網(wǎng)格參數(shù)計算，回到第1步。其中主要涉及的UDF源文件定義如下：

對應(yīng)第1步：源文件1：

DEFINE_INIT()

{定義箭體姿態(tài)運動初始條件}

對應(yīng)第2步：源文件2：

DEFINE_ADJUST()

{加載所需參數(shù)，包括作用力矩，姿態(tài)運動信息，姿態(tài)運動變化率，中間參數(shù)；

判斷循環(huán)條件；

姿態(tài)動力學與運動學運動方程計算、疊加；

返回姿態(tài)運動變化參數(shù)，姿態(tài)運動信息；}

對應(yīng)第3步，源文件3：

DEFINE_SOURCE()

{Real 加載所需參數(shù)；

調(diào)用加速度定義函數(shù)；

Return 加速度；}

仿真開始后，當網(wǎng)格計算一步長時，源文件1運行，將姿態(tài)運動方程的初始化條件Tx，Ty和Tz載入；之后，源文件2運行，計算姿態(tài)動力學和運動學方程，得出這一步箭體的姿態(tài)運動狀態(tài)；最后，源文件3調(diào)用上一步得到的姿態(tài)運動狀態(tài)參數(shù)，計算流場模型需要加載慣性加速度。每一步UDF計算開始時，都需要從網(wǎng)格載入上次循環(huán)時在流體系統(tǒng)計算中得到的作用力矩Tx，Ty和Tz，通過上述3個源文件再得到下一步的網(wǎng)格場加速度。將下一個時刻的網(wǎng)格場加速度a1輸出到網(wǎng)格邊界條件中，進行下一步網(wǎng)格計算。循環(huán)這一過程，就實現(xiàn)了Fluent計算和姿態(tài)動力學與運動學計算之間的實時同步數(shù)據(jù)傳輸過程。

2 仿真結(jié)果與分析

圖4 液體晃動對航天器的作用力矩

由于液體的振動以及章動運動，S′系中y′軸方向上的作用力矩振蕩上升，始終在其正方向上。即，在Fluent仿真的網(wǎng)格計算模型中，通過自定義源函數(shù)加載的慣性加速度始終指向離心方向，符合客觀規(guī)律。通過Fluent流體系統(tǒng)仿真模型還可以看到運動過程中各個時間點的流場分布情況，如圖5。

圖5 不同時刻流場分布圖

隨著貯箱內(nèi)的液體不斷地向一側(cè)壁面碰撞，產(chǎn)生對該側(cè)的作用力。圖7中左側(cè)的貯箱離質(zhì)心較遠，這種碰撞作用力產(chǎn)生的力矩效果更加明顯。右側(cè)的貯箱內(nèi)包含有坐標原點，液體因表面張力的作用環(huán)繞壁面流動，產(chǎn)生對壁面的剪應(yīng)力，以及對質(zhì)心的作用力矩。

圖5中是根據(jù)時間順序選擇的3幅典型的流場分布圖，其中左側(cè)貯箱內(nèi)，如果對照圖5中各個狀態(tài)下的流場分布情況，便可輕易地分析圖4(a)中的作用力矩變化的原因。左側(cè)貯箱的液體越來越多地與貯箱的一側(cè)碰撞，使碰撞力的作用增大；右側(cè)貯箱的液體分布逐漸均勻，與貯箱壁面的接觸面積增大，摩擦力矩也隨之增大。因此，圖4(a)中的力矩變化剛好與圖5中的流場分布變化相吻合。

原來羅漠與任何人都是一樣的，他的獨特，他的深情，都只為一人而已。他的愛那么忠誠，經(jīng)得起各種各樣的考驗，她其實是滿意的，可是她忍不住流了一路的眼淚。

在這種作用下，箭體的姿態(tài)運動受到的干擾應(yīng)當是越來越大的，下面再對箭體系統(tǒng)進行分析。圖4(a)中將S′系中的作用力矩投影到本體系Sb中，其方向與本體系當時的相對位置有關(guān)。箭體系統(tǒng)受到的液體干擾力矩就是圖4(b)中的作用力矩。根據(jù)本文第1.1節(jié)中箭體動力學方程，得到在剩余液體干擾下箭體的姿態(tài)運動狀態(tài)，如圖6所示。

圖6 液體干擾作用下的箭體自旋軸空間指向

從圖6分析整個系統(tǒng)的運動狀態(tài)：在仿真時間內(nèi)，箭體系統(tǒng)章動逐漸增加，因此旋轉(zhuǎn)離心作用增大，必然導致液體相對于箭體運動加強。根據(jù)液體運動對質(zhì)心的力矩分析可知，液體干擾作用力矩也應(yīng)當是逐漸增強的，與Fluent網(wǎng)格模型的仿真結(jié)果一致。并且隨著箭體運動穩(wěn)定性的惡化，繞最小慣量軸(自旋軸Xb軸)的自旋運動最終轉(zhuǎn)化為繞最大慣量軸(Yb軸或Zb軸)的不規(guī)則旋轉(zhuǎn)運動，因而箭體運動發(fā)散后在S′系中y′軸方向上的作用力矩出現(xiàn)下降并且幅值也產(chǎn)生不規(guī)則的變化。

本文選擇的仿真方法不僅可以測得流場與箭體的相互影響參數(shù)，通過編寫合理的UDF程序，還可以選擇模型中任意位置作為觀測點，測量其在流場中所受的作用力矩，或者對流場內(nèi)的壓強、溫度、過載進行中間控制，通過適當?shù)目刂瞥绦蛟谙鄳?yīng)時間修改這些參數(shù)。并且，也可以將包含控制判斷算法的程序以UDF程序的形式添加到閉環(huán)系統(tǒng)中，形成含有實時控制系統(tǒng)的姿態(tài)耦合運動仿真模型。

這種方法的優(yōu)點在于可以直觀形象地觀察到仿真各個時刻系統(tǒng)的變化和狀態(tài)，便于分析影響姿態(tài)運動的各種因素。同時，這種選擇離心運動坐標系為仿真參考系的方法可以簡化箭體姿態(tài)運動對液體影響的描述，通過簡單的投影關(guān)系就可以得到清晰、解耦的作用關(guān)系。

由此可知，F(xiàn)luent-UDF技術(shù)適用于更復雜的分布參數(shù)系統(tǒng)，可以推廣到更復雜的包含流體計算的系統(tǒng)控制模型中。

3 結(jié)論

[1] S.Dutta, M.K.Laha.Analysis of the Small Amplitude Sloshing of a Liquid in a Rigid Container of Arbitrary Shape Using a Low-order Boundary Element Method[J].International Journal for Numberical Methods in Engineering, 2000, 1(47): 1633-1648.

[2] 王照林, 李磊.航天器內(nèi)部液體晃動對交會對接動力學與控制的影響[J].航天控制, 1991,(2): 24-32.(WANG Zhaolin, LI Lei.Effect on Dynamics and Control of Rendezvous and Docking due to Liquid Sloshing Internal to Spacecraft[J].Aerospace Control, 1991, (2): 24-32.)

[3] 耿浩, 蔣志文.基于網(wǎng)格的分布式仿真技術(shù)[J].航天控制, 2004, 22 (6): 46-48.(GENG Hao, JIANG Zhiwen.Distributed Simulation Technology Based on the Grid[J].Aerospace Control, 2004,22 (6): 46-48.)

[4] 鮑文, 李偉鵬, 等.基于Fluent/Matlab接口的分布參數(shù)系統(tǒng)閉環(huán)控制仿真[J].北京:系統(tǒng)仿真學報, 2008, 20 (11):2851-2854.(BAO Wen, LEI Weipeng, et al.Close-loop Control Simulation Technology of Distributed Parameter System Based on FLUENT/MATLAB Interface[J].Journalof System Simulation, 2008,20 (11): 2851-2854.)

[5] 肖業(yè)倫.航天器飛行動力學原理[M].北京: 宇航出版社, 1995: 212-217.(XIAO YeLun.Theory of Spacecraft Dynamics[M].Chinese Astronautics Press, 1995: 212-217.)

[6] 廖少英.低重力或失重環(huán)境下剩余推進劑的管理和排放技術(shù)研究[J].上海航天, 1993, 5: 13-18,40.(LIAO Shaoying.Research of the Residual Propellant Dumping System under Micro-gravity[J].Aerospace Shanghai, 1993, 5: 13-18,40.)

《航天控制》雜志

□歡迎訂閱 □歡迎刊登廣告 □歡迎投稿

開本：大16開 雙月刊 頁碼：96頁 定價：20.00元(全年120元)

國內(nèi)郵發(fā)代號：80-338 國外發(fā)行代號：BM4668

國內(nèi)統(tǒng)一連續(xù)出版物號：CN11-1989/V 國際標準連續(xù)出版物號： ISSN 1006-3242

編輯部地址：北京142信箱402分箱 郵政編碼：100854

電話：(010) 68762264、68388585 E-mail: ht12@httx.com.cn

網(wǎng)址：http://htkz.chinajournal.net.cn

The Analysis and Simulation of Attitude Motion for a De-orbit Rocket Based on Fluent Flow Field Calculation

GUAN Hong XU Shijie

School of Astronautics, Beihang University, Beijing 100191, China

Duetotheattitudemotionanalysisofaspacecraft,thefluxionandsloshingofliquidarefrequentlyencountered.Sinceitisdifficulttoestablishaprecisemathematicalmodelforsuchissues,researchersprefertousenumericalcalculationandsimulationforanalysis.Inthispaper,theUDFisusedtointroducetheuserdefinedvariablefunctionandtherocketde-orbitcontrolwithflowfieldsimulationofresidualfuelbasedonfluentistokenforexampletoanalyze.Byaddingtransportmotiontothefluidmotionequations,theparametersanddisturbingtorqueofflueundertheon-orbit-statecanbeobtained.Theoperationparametersofattitudeanalysisandsimulationforlaunchvehiclesystemareprovidedsothatthechangesinflowandtheattitudemotionofthelaunchvehiclearecorrelated.

Fluent;De-orbit;Residualfuel;User-definefunction(UDF)

2011-08-16

關(guān) 宏(1985-)，女，北京人，博士研究生，主要研究方向為航天飛行器動力學與控制；徐世杰(1951-)，男，吉林人，教授，博士生導師，主要研究方向為航天飛行器動力學與控制。

V525

A

1006-3242(2012)03-0093-06