張博譞，李國岫，虞育松

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

0 引言

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

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

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

1 基本方程

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

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

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

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

式中：n為法向向量。

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

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

式中：Gk為由平均速度梯度引起的湍動(dòng)能；Gb為用于浮力影響引起的湍動(dòng)能；YM為可壓速湍流脈動(dòng)膨脹對(duì)總的耗散率的影響；C1ε，C2ε，C3ε為經(jīng)驗(yàn)常數(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)確模擬剪切層中平均場(chǎng)流動(dòng)方向的改變對(duì)湍流場(chǎng)的影響。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語

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

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

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

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

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

［5］ 李俊峰，魯 異，寶音賀西，等.貯箱內(nèi)液體小幅晃動(dòng)的頻率和阻尼計(jì)算［J］.工程力學(xué)，2005，22（6）：87-90.

［6］ 王 毅，常小慶.微重力環(huán)境下推進(jìn)劑貯箱中三維氣液平衡界面的數(shù)值模擬［J］.火箭推進(jìn)，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)液體晃動(dòng)阻尼研究（二）：數(shù)值計(jì)算［J］.宇航學(xué)報(bào)，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.