基于Fluent計算的火箭離軌姿態運動仿真與分析

關 宏 徐世杰

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

?

基于Fluent計算的火箭離軌姿態運動仿真與分析

關 宏 徐世杰

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

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

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

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

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

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

1 仿真系統的建立

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

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

1.1 箭體系統動力學分析

首先以箭體系統作為研究對象，其轉動慣量沿自旋軸對稱，箭體本體坐標系Sb定義如圖1所示，坐標原點Ob在箭體質心，Xb軸沿本體縱向對稱軸方向，Yb軸沿箭體縱切面向上，Zb軸按右手法則定義。

圖1 箭體本體坐標系定義

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

(1)

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

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

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

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

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

圖2 本體錐和空間錐

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

1.2 液體子系統運動

得到箭體空間運動規律，再考慮液體系統的相應晃動影響。

液體子系統在微重力環境下運動，受擾動后不能斷定其確切位置，液氣界面穩定性差。一旦航天器進行姿態控制(起旋至穩定旋轉)，剩余液體產生的擾動就會影響到姿態控制的精度，破壞航天器的運動穩定性。又因為液體運動規律不易簡化，如不能準確仿真其運動狀態，便無法找到適當的控制策略，會對回收造成惡劣的影響[6]。

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

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

(2)

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

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

(3)

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

(4)

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

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

1.3 Fluent-UDF接口的設計

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

圖3 使用UFD的流體模型與姿態運動模型交互示意圖

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

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

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

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

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

對應第1步：源文件1：

DEFINE_INIT()

{定義箭體姿態運動初始條件}

對應第2步：源文件2：

DEFINE_ADJUST()

{加載所需參數，包括作用力矩，姿態運動信息，姿態運動變化率，中間參數；

判斷循環條件；

姿態動力學與運動學運動方程計算、疊加；

返回姿態運動變化參數，姿態運動信息；}

對應第3步，源文件3：

DEFINE_SOURCE()

{Real 加載所需參數；

調用加速度定義函數；

Return 加速度；}

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

2 仿真結果與分析

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

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

圖5 不同時刻流場分布圖

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

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

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

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

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

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

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

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

由此可知，Fluent-UDF技術適用于更復雜的分布參數系統，可以推廣到更復雜的包含流體計算的系統控制模型中。

3 結論

[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] 王照林, 李磊.航天器內部液體晃動對交會對接動力學與控制的影響[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] 耿浩, 蔣志文.基于網格的分布式仿真技術[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接口的分布參數系統閉環控制仿真[J].北京:系統仿真學報, 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] 肖業倫.航天器飛行動力學原理[M].北京: 宇航出版社, 1995: 212-217.(XIAO YeLun.Theory of Spacecraft Dynamics[M].Chinese Astronautics Press, 1995: 212-217.)

[6] 廖少英.低重力或失重環境下剩余推進劑的管理和排放技術研究[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元)

國內郵發代號：80-338 國外發行代號：BM4668

國內統一連續出版物號：CN11-1989/V 國際標準連續出版物號： ISSN 1006-3242

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

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

網址：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

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

V525

A

1006-3242(2012)03-0093-06