變軌條件下衛(wèi)星貯箱內液體推進劑晃動特性三維數(shù)值模擬

張博譞，李國岫，虞育松

（北京交通大學 機電學院，北京 100044）

0 引言

在微重力或零重力環(huán)境中，充液航天器內液體晃動及其對控制系統(tǒng)的影響是當前國內外航天技術研究的重要課題。微重力條件下的液體流動行為研究始于20世紀60年代，研究對象是火箭推進劑罐內的液體晃動動力學特性。隨后，低重力條件下的在軌空間飛行器燃料貯箱內的液體晃動特性被重點研究。隨著航天事業(yè)的發(fā)展，推進劑占火箭、衛(wèi)星等航天器總質量的比重不斷加大，推進劑晃動與箭體姿態(tài)運動通過慣性力的相互作用力而直接耦合?；蝿硬环€(wěn)定，晃動幅值不斷增大，火箭姿態(tài)角就會因晃動不斷加大，從而導致箭體姿態(tài)運動發(fā)散［1］。充液航天器總體設計中必須考慮液體晃動，推進劑的運動對在軌航天器的穩(wěn)定、定位、對接的影響相當關鍵［2］。在星體充液量不斷增大和控制飛行器姿態(tài)的指向精度條件下，并考慮零重力或微重力環(huán)境，充液體晃動動力學的研究顯得更重要和復雜［3］。

在液體的晃動對航天器的穩(wěn)定性研究中，為便于分析，大量研究采用了等效單擺模型和質量－彈簧力學模型［4］。分析中進行了簡化，忽略了部分力的影響，故與微重力或零重力下的實際狀態(tài)相比，仍有一定差異。對微重力或零重力下的液體晃動的研究多為定體積，而航天器在變軌過程中必然使用液體推進劑，使貯箱內液體推進劑體積減小，且微重力或零重力下的液體晃動是一個非線性問題，這導致變軌過程中的液體晃動問題的研究更復雜。工程中最關心晃動頻率和阻尼兩個參數(shù)。在充液衛(wèi)星姿態(tài)控制中，計算晃動頻率以在設計中避免壁面液體晃動、姿態(tài)運動和彈性附件振動的共振，晃動過程中的阻尼與章動發(fā)散時間直接相關。在衛(wèi)星入軌后，姿態(tài)控制系統(tǒng)必須在一定時間內及時開啟，防止衛(wèi)星章動發(fā)散失穩(wěn)［5］。采用三維數(shù)值模擬對液體流動過程進行直接模擬，不作簡化假設，對液體晃動的阻尼比和頻率的預測更準確，同時能得到更直觀的流場、壓力、速度分布圖，便于對流體流動過程做更進一步的研究。

本文用Fluent三維數(shù)值模擬軟件，以某衛(wèi)星液體推進劑貯箱為研究對象，在零重力環(huán)境中，選取不同充液比，對衛(wèi)星變軌過程液體的分布及重定位進行預測。

1 基本方程

不相容氣液兩相牛頓流體的可壓縮N-S方程為

式中：ρ為流體密度；p為流體微元上的壓力；I為單位對稱矩陣；αl，αv分別為兩相的體積分數(shù)；ρl，ρv分別為兩相的密度；u為速度；F為表面張力；g為重力。

瞬態(tài)、不可壓、不相溶、等溫、定黏度且存在表面張力（不相溶）的兩相牛頓流體Navier-Stokes方程為

式中：τ為流體的黏性應力；S為應變率張量。式（6）、（7）考慮了壓力、湍流、表面張力和重力對流動的影響。由于氣液兩相不相溶，相與相間存在一個界面，故相界面上存在法向的F。F與液體性質，液面曲率相關，有

式中：n為法向向量。

為建立準確的近壁面數(shù)學模型，引入重整化群（RNG）k-ε雙方程湍流模型。

標準k-ε模型中湍動能k和耗散率ε可表示為

式中：Gk為由平均速度梯度引起的湍動能；Gb為用于浮力影響引起的湍動能；YM為可壓速湍流脈動膨脹對總的耗散率的影響；C1ε，C2ε，C3ε為經(jīng)驗常數(shù)，且C1ε＝1.41～1.45，C2ε＝1.90～1.92，C3ε＝0.8。

標準k-ε模型也有一定的局限性，主要表現(xiàn)在：采用了Boussinesq假定，即采用了雷諾應力與平均速度梯度線性關聯(lián)的梯度型和湍流黏性系數(shù)各向同性的概念，使k-ε模型難以準確模擬剪切層中平均場流動方向的改變對湍流場的影響。

重整化群k-ε模型是對瞬時的 Navier-Stokes方程用重整化群的數(shù)學方法推導而得。模型中的常數(shù)與標準k-ε模型不同，且方程中也出現(xiàn)了新的函數(shù)或項。其湍動能及耗散率方程的形式與標準k-ε模型相似。重整化群k-ε模型有助于處理低雷諾數(shù)和近壁流動問題的模擬。RNGk-ε模型的可信度和精度較標準k-ε模型更高。因此，本文的數(shù)值模擬研究采用RNGk-ε模型。

不相容的氣液兩相，相與相間存在界面。不同流體分子間作用力的差異將導致相界面上出現(xiàn)張力。利用體積分數(shù)（VOF）模型能只用一組N-S方程就可描述氣液兩相及相界面，使問題得以簡化。位于相界面處的單元內氣液兩相同時共存，非相界面單元只存在單相。因此，引入k相的體積分數(shù)εk，某個單元內εk可表示為

整個流場單元可被分為3個區(qū)域：εk（cell）＝0時，單元內不存在k相流體；εk（cell）＝1時，單元內只存在k相流體；0＜εk（cell）＜1時，單元內有相界面存在。

在微重力或零重力環(huán)境中，毛細力作用占主導地位，液體表面高度彎曲，且貯箱結構對液體定位的影響非常大。當邦德數(shù)Bo≤1時，毛細力占主導地位。Bo為慣性力與毛細力的比值，有

式中：σ為液體表面張力；a為加速度；R為貯箱特征尺寸［5］。文獻［7］給出了圓柱容器內液體小幅晃動第一階模態(tài)阻尼比的經(jīng)驗公式

式中：υ為液體運動黏性系數(shù)；R為容器半徑；h為容器內的液深。當h／R＞1時，阻尼比基本不隨液深變化。

2 計算結果及分析

推進劑甲基肼（MMH）的密度為875kg／m3（20℃），黏性系數(shù)為0.000 855Pa·s，表面張力系數(shù)為0.034N／m。某衛(wèi)星推進劑貯箱結構如圖1所示。考慮模型為對稱結構，為節(jié)省計算時間，減小計算區(qū)域，將模型縱剖，推進劑貯箱的一半作為計算區(qū)域。將剖面設置為對稱面。對衛(wèi)星某次變軌過程行數(shù)值模擬，變軌過程初始條件為：變軌過程加速度3.78×10－3m／s2；微重力方向為－Z向；充液比0.89；變軌過程中的液體流出流量34.066g／s。

用Fluent三維數(shù)值模擬軟件，對零重力環(huán)境中貯箱內液體推進劑分布及重定位過程進行了數(shù)值模擬。在零重力條件下，對貯箱有液體推進劑流出時的變軌過程進行了數(shù)值模擬計算。

變軌工況下液體質心變化的計算結果如圖2所示。由圖可知：整個流動過程的時間約4 000s；變軌過程初期（0～1 000s），貯箱內液體發(fā)生較大振幅的阻尼振蕩，液體晃動幅值逐漸減小，同時液體質心位置隨出箱內液體量的減少（液體的Z向質心位置線性下降）。變軌過程后期（1 000～4 000s），液體晃動幅度逐漸趨近于零，液體質心由于液體量減少而呈線性下降。

圖1 貯箱結構Fig.1 Structure of tank

圖2 變軌工況下液體質心（Z向）變化Fig.2 Change of gravity center of liquid（Zdirection）

變軌過程中，時刻t＝0，136，1 336，4 180s和變軌結束后液面穩(wěn)定狀態(tài)下的軸線截面以及三維液面透視圖如圖3所示。

由于氣液固界面處存在毛細力作用，在氣液固界面處液面會與固體壁面呈一定角度，稱為接觸角。接觸角是在固、液、氣三相接觸達到平衡時，三相接觸周邊的任一點上，液氣界面切線與固體表面間形成的并包含液體的夾角。該值與固體、液體和氣體屬性均有關系。不同情況下，接觸角可為0°～180°。文獻［8］用光學張量測量法得到MMH與鋁制容器間的接觸角在10°附近。本文將MMH與鋁制容器間的接觸角定在10°［8］。

因存在接觸角，微重力或零重力下貯箱內液體推進劑的液面不會為水平液面，液面會在毛細力與表面張力的共同作用下沿壁面爬升，直至固液氣三相接觸達到平衡。上述算例中的變軌零時刻液面為數(shù)值模擬形成，即將初始液面設置為水平液面，使液面在完全失重條件下沿壁面自由爬升，達到穩(wěn)定平衡狀態(tài)（質心穩(wěn)定不再變化，且液面基本停止波動），并將此時刻作為變軌過程的初始時刻（圖3（a））。

表1 兩次變軌過程的液體晃動阻尼比和頻率Tab.1 Damping ratio and frequency of liquid sloshing during two orbit changes

變軌過程中貯箱內液體不斷流出，導致液面持續(xù)下降（貯箱內液深持續(xù)變化），且加上變軌初期液面有較大幅度晃動，故數(shù)值模擬方法與式（13）求得的阻尼比均為變動值。

不同液深時刻阻尼比的數(shù)值模擬計算值與經(jīng)驗公式的計算值如圖4所示。由圖可知：本文的計算結果與經(jīng)驗公式計算結果非常接近。國內學者曾對液體晃動的阻尼比進行簡化計算，利用速度勢函數(shù)推導的計算液體晃動頻率和阻尼的特征值方程，簡化處理轉化為一般的廣義特征值。其結果為充液比為0.5時，ζ為5.67×10－4［5］。此結果與由阻尼比經(jīng)驗公式得出的值相差非常大，故可認為簡化計算對阻尼的預測結果并不可靠。本文模擬能較準確地得出的變軌過程中液體晃動的阻尼比。

圖4 不同液面高度晃動阻尼比計算結果對比Fig.4 Calculation results of damping ratios with various liquid height ratios

變軌過程初始時刻貯箱內流線如圖5所示。由圖可知：貯箱內液體推進劑甲基肼與氦氣的流動趨勢，甲基肼為Z軸負向流動。這一方面是貯箱模型底部有甲基肼質量流量出口，另一方面是變軌過程加速度為Z軸負向，導致液體推進劑會在變軌過程中有下沉的趨勢。

變軌過程時刻136s貯箱內流線和速度場分別如圖6、7所示。由圖6、7可更直觀地觀察在變軌過程中貯箱內甲基肼與氦氣的流動狀態(tài)，甲基肼在貯箱近壁面處與貯箱中部的速度方向相反，即貯箱內甲基肼已開始在Z軸方向呈現(xiàn)震蕩。

圖5 變軌初期貯箱內流線Fig.5 Path line of flow in tank at begin of orbital transfer

圖6 變軌時刻136s貯箱內流線Fig.6 Path line of flow in tank at 136sof orbital transfer

圖7 變軌時刻136s貯箱內速度場Fig.7 Velocity field in tank at 136sof orbital transfer

甲基肼對貯箱的作用力如圖8所示。由圖可知：變軌開始后貯箱受到的作用力迅速增至約1.6N，作用力為Z軸負向；隨后作用力出現(xiàn)衰減振蕩；貯箱受到的作用力在變軌初期振蕩衰減，隨著變軌過程的進行以及貯箱內甲基肼體積的減少，作用力大小呈線性降低。

圖8 變軌貯箱Z向受力變化過程Fig.8 Change of force of tank inZdirection during orbit change

貯箱內甲基肼在變軌前后質心的變化如圖9所示。一方面，變軌過程中甲基肼以恒定的質量流量流出貯箱，導致變軌后甲基肼的體積較變軌前不斷減少，另一方面，考慮表面張力的作用，液體在晃動中會有發(fā)散趨勢，即隨著時間的增加，自由液面不再回到其初始靜平衡位置［3］。

圖9Z向質心穩(wěn)定值Fig.9 Stable value of gravity center inZdirection

3 結束語

本文以某衛(wèi)星作為研究對象，對變軌條件下貯箱內液體推進劑晃動特性進行了三維數(shù)值模擬。獲得了衛(wèi)星變軌過程中，液體推進劑的質心變化、液體推進劑對貯箱作用力變化、液體推進機晃動一階頻率以及液體推進劑的一階晃動阻尼比，通過貯箱內的流場分布圖和速度場圖，直觀地給出變軌過程中貯箱內推進劑的流動過程及流動的變化趨勢。由本文數(shù)值模擬結果可知，考慮表面張力時液體晃動具有發(fā)散的趨勢，即隨著時間的增加，自由液面不再回到其初始靜平衡位置，與文獻［3］中的結論吻合。隨時間的增加，貯箱壁面處的液體平衡位置有所增加，這說明將有更多液體依附壁面，從而引起液體晃動質量的降低。將本文得出的阻尼比數(shù)值模擬計算值與由實驗值擬合出的經(jīng)驗公式得出的阻尼比計算值進行對比?？砂l(fā)現(xiàn)兩者非常接近，即在對衛(wèi)星變軌過程中，衛(wèi)星貯箱內液體推進劑的阻尼比預測準確度較高。

［1］ 尹立中，王本利，鄒經(jīng)湘.航天器液體晃動與液固耦合動力學研究概述［J］.哈爾濱工業(yè)大學學報，1999，31（2）：118-122.

［2］ 陳存蕓.運載火箭三級無動力飛行段晃動穩(wěn)定性研究［J］.上海航天，2004，21（4）：29-33.

［3］ 岳寶增，劉延柱.低重力環(huán)境下三維液體非線性晃動的數(shù)值模擬［J］.宇航學報，2000，21（4）：25-30.

［4］ 包光偉，劉延柱.三軸定向充液衛(wèi)星的姿態(tài)穩(wěn)定性［J］.空間科學學報，1993，13（1）：31-38.

［5］ 李俊峰，魯 異，寶音賀西，等.貯箱內液體小幅晃動的頻率和阻尼計算［J］.工程力學，2005，22（6）：87-90.

［6］ 王 毅，常小慶.微重力環(huán)境下推進劑貯箱中三維氣液平衡界面的數(shù)值模擬［J］.火箭推進，2007，33（3）：31-35.

［7］ ABRAMSON H N.The dynamic behavior of liquids in moving containers［R］.NASA SP 106，1966.

［8］ 王 為，李俊峰，王天舒.航天器貯箱內液體晃動阻尼研究（二）：數(shù)值計算［J］.宇航學報，2006，27（2）：177-180.

［9］ KREPPEL S.Scaling and modeling of propellant sloshing and zero gravity equilibrium for the orion service module propellant tanks［D］.Kenosha：Carthage College，2010.