張博譞，李國岫，虞育松

（北京交通大學(xué) 機電學(xué)院，北京 100044）

0 引言

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

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

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

1 基本方程

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

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

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

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

式中：n為法向向量。

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

標(biāo)準(zhǔn)k-ε模型中湍動能k和耗散率ε可表示為

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

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

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

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

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

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

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

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

2 計算結(jié)果及分析

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

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

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

圖1 貯箱結(jié)構(gòu)Fig.1 Structure of tank

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語

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

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

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

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

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

［5］ 李俊峰，魯 異，寶音賀西，等.貯箱內(nèi)液體小幅晃動的頻率和阻尼計算［J］.工程力學(xué)，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］ 王 為，李俊峰，王天舒.航天器貯箱內(nèi)液體晃動阻尼研究（二）：數(shù)值計算［J］.宇航學(xué)報，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.