嚴 天 ，王 磊 ，2*，王嬌嬌 ，黃曉寧 ，馬 原 ，2，厲彥忠 ，2

（1.西安交通大學制冷與低溫工程系，西安 710049；2.航天低溫推進劑技術國家重點實驗室，北京 100028）

0 引言

在我國未來航天戰(zhàn)略布局中，月球探測及其他深空探測將成為主要攻關方向[1]。深空探測任務的順利開展要求火箭具備長時間在軌運行及變軌能力。低溫液氫、液氧推進劑具有高比沖、無毒無污染等優(yōu)點，將在探測中扮演重要角色。已有研究指出[2]，低溫推進劑長期在軌存儲與管理技術是未來航天任務發(fā)展的重要支撐技術之一。與常重力不同，在軌重力影響變弱時涉及多種弱力的綜合作用，如氣液黏性力、液體表面相對增強的張力等。因此在推進劑管理中，必須考慮低溫推進劑在軌時的運行特性及規(guī)律。在軌期間，飛行器須完成的多次變軌、姿態(tài)調整和氣液分離過程，表征為多種激勵形式，包括軸向變過載、側向激勵、慢旋等。如不對推進劑行為加以控制，將導致貯箱底部管路口暴露，氣體進入發(fā)動機，對飛行器的穩(wěn)定產生不利影響。同時箱內液體晃動和氣液相分布紊亂也不利于姿軌控制。因此，有必要對在軌狀態(tài)下推進劑受激勵后的動力學特性進行分析研究，掌握其規(guī)律，從而為推進劑管理設計及安全防護提供理論支撐。

對于推進劑晃動問題，國內外已有學者采用理論分析、數(shù)值模擬和實驗測量等不同方法開展了相關研究。理論研究領域，Abramson[3]和Dodge[4]提出了用單擺和彈簧振子兩種等效力學模型來分析容器內的流體晃動特征，并通過數(shù)學方法分析預測了流體動力學行為。在動力學數(shù)值分析方面，Himeno等[5]、Aoki等[6]、Yang等[7]和Beckman等[8]通過仿真方法證實了常重力下?lián)醢鍖α黧w晃動具有抑制作用。Godderidge等[9]比較了實驗和模擬結果的壓力峰值，建議采用非均勻模型模擬劇烈晃動。Cheng等[10]提供了一種簡化方法用于計算含擋板貯箱的動壓特性。Shekari等[11]通過有限元分析發(fā)現(xiàn)，擋板高度、隔振基礎的剛度等參數(shù)對最大流體動壓和液體振動受力等有顯著影響。Kargbo等[12]通過數(shù)值建模研究了分層晃動波特性，并按時間序列分別對界面下層和上層自由表面的液體晃動進行了快速傅里葉變換光譜分析，獲得了流體層的能量和頻率。Weidman[13]通過二維晃動特性模擬，研究了晃動頻率隨液體深度的變化規(guī)律。對推進劑晃動的熱力學行為，也有學者開展了相關研究，Zuo等[14]關注了晃動過程中擋板對低溫流體自增壓及熱分層的影響。Liu等[15-16]建立了一種非絕熱氧箱模型，使用動網(wǎng)格捕捉氣液界面運動，研究了液體初始溫度、振蕩幅度對貯箱壓力分布、振蕩力矩和界面動態(tài)響應的影響。實驗研究方面，Konopka等[17]通過實驗對比驗證了一種與擋板等效的主動控制算法的有效性。Vairamani等[18]實驗研究了磁場作用下半剛性結構膜對晃動的抑制效果。Pérez 等[19]、Himeno 等[5]、Aoki等[6]、Stephens 等[20]和Pridgen等[21]也開展了涉及推進劑晃動的實驗研究，揭示了常重力下激勵、液位高度及擋板參數(shù)對晃動過程的影響。目前對流體在外加激勵下的研究主要針對常重力或大過載環(huán)境開展，對空間微重力下流體在外加激勵下行為的實驗研究相對不足。在實驗研究方面，難以獲得同時具備較長時間和較低加速度的微重力環(huán)境[22]，對激勵發(fā)生時液面失穩(wěn)乃至傾覆的事故，也鮮有實驗數(shù)據(jù)作為支撐。

本文采用計算流體力學（CFD）軟件FLOW-3D，以某型全尺寸Cassini等[23]低溫貯箱為對象，構建可預測低溫貯箱空間多種激勵條件下流體行為的仿真模型。數(shù)值仿真空間微重力下初始傾斜界面、側向正弦波、慢旋激勵下箱內兩相流體的運動特性，系統(tǒng)揭示微重力條件對流體行為演化規(guī)律的影響，并考核箱內設置防晃隔板對流體達到新平衡的作用效果。本文研究有助于理解空間低重力下低溫貯箱內的流體運動規(guī)律，為低溫上面級控制方案的設置、時序的優(yōu)化提供理論指導。

1 模型建立

1.1 物理模型

圖1展示的貯箱由一個直徑為D，半徑為R0，高度為h2的圓柱以及分別位于圓柱頂部和底部，高度為h1的橢球形封頭組成。具體參數(shù)如表1所列。圖1坐標系原點取為貯箱底部最低一點O，xyz軸向如圖所示。液體充灌率為0.95，對應液位高度為H。箱內裝有6組環(huán)形擋板，間隙為d。擋板用于抑制上升階段液體的過量防晃。取擋板寬度為0.1D[5-6]。貯箱封頭底部設有十字交叉擋板，用于防止底部液體塌陷。此外，圖1展示了自由界面傾斜所產生的晃動，徑向激勵（文中將其簡化為正弦函數(shù)），沿軸向的慢旋激勵（文中將其簡化為恒定角速度）的作用方向。圖中標星號的位置坐標為（1.9 m，1.9 m，3.4 m），用于監(jiān)測壓力的變化。

表1 貯箱參數(shù)Tab.1 The parameters of tank

圖1 貯箱結構示意圖Fig.1 The structure of tank

1.2 CFD模型

貯箱模型及網(wǎng)格如圖2所示。在FLOW-3D中，先建立方形固體框架，再在中間挖出與貯箱等體積的空槽，通過對空槽的不同區(qū)域設置氣液相，并導入擋板等固體結構實現(xiàn)建模。貯箱內所有的流體區(qū)域均被劃分為尺寸一致的結構化網(wǎng)格。圖中，深藍色區(qū)域代表貯箱，內含流體區(qū)和擋板，由于兩部分固體結構相互內切，圖2（a）中模型只顯示外側箱壁，具體結構如圖2（b）、（c）所示。

圖2 網(wǎng)格模型及內部結構圖Fig.2 The model of mesh and the internal structure

模型設置如表2所列。其中，采用RNG k-ε湍流模型，以準確描述近壁區(qū)的強剪切效應。在軌期間，為了實現(xiàn)氣液相分離，滿足供液與排氣的需求，通常采用正推沉底發(fā)動機提供小過載正推，所產生的加速度約為10-3g0?；鸺掌陂g，不可避免地存在振動而引起箱內氣液相晃動，當火箭入軌由大過載瞬變?yōu)槲⒅亓?，升空階段的界面晃動可能對后期的兩相演化產生影響。為此，本文將其簡化為界面初始傾角機制，設置相界面與水平面夾角為θ，當不考慮該擾動作用時，θ=0°?？臻g低溫貯箱常采用多層絕熱（MLI）作為絕熱手段，對應的壁面漏熱在1 W/m2量級，在該漏熱條件下，本文研究對象所產生的流體溫升僅有3.3×10-6K，故可忽略壁面漏熱影響。

2 模型校驗

2.1 網(wǎng)格無關性校驗

依照表2中參數(shù)設置，以液面傾斜4°為初始條件，對模型進行網(wǎng)格無關性驗證，結果如圖3所示。將網(wǎng)格數(shù)由5萬逐漸增加至45萬，對應網(wǎng)格尺寸從0.1 m減小至0.045 m，并與達平衡后氣枕中監(jiān)測點的壓力進行比較。由圖可見，網(wǎng)格數(shù)到30萬后計算結果達到穩(wěn)定。為了確保計算精度，以41萬網(wǎng)格數(shù)模型開展后續(xù)的仿真計算。

表2 模型及初始參數(shù)設置Tab.2 The set of model and initial parameters

圖3 網(wǎng)格數(shù)對仿真結果的影響曲線Tab.3 The effect of mesh on the simulation results

2.2 CFD模型校驗

采用文獻[5]報道的晃動實驗校核本文所構建的CFD模型。該實驗以水為工質，在常重力下展開，施加激勵如圖4所示。貯箱沿軸線（Z方向）受常重力加速度作用，沿徑向（X方向）受外部正弦波激勵作用，初始相位角為-180°，最大振幅為12 m；頻率為1.47 Hz，振幅為5.32 m。過載加速度a與振幅A、頻率ω的關系如式（1）所示：

圖4 施加激勵的參數(shù)曲線Fig.4 The parameters of excitation

式中：t為時間；φ為相位角。

仿真結果與實驗結果的對比如圖5所示?？梢钥闯觯芗詈?，實驗貯箱和模擬貯箱內流體呈現(xiàn)相似行為，晃動的幅度和流體翻轉、破碎等行為均相近，證實本文所構建的模型能夠較精確地捕捉氣液相界面的分布及動態(tài)規(guī)律，證明了該模型對激勵工況下流體運動仿真的準確性。

圖5 仿真結果與實驗結果對比圖Fig.5 Comparison of simulation results and experimental results

3 仿真結果與分析

3.1 傾斜液位下的晃動特性

為衡量晃動的劇烈程度，本文取Z軸質心位移開展特性研究，質心由式（2）得到。定義液氧全部位于貯箱頂部與液氧全部位于底部對應的質心差為所允許的最大晃動幅度，為15.0 cm。本文以實際的質心變化范圍相較于最大幅值的比值進行分析。

式中：z為質心垂直方向坐標，m；Ω為體積，m3；ρ為流體密度，kg·m-3。

圖6（a）、（b）、（c）分別展示了10-3g0下液氧初始傾斜角度θ分別為1°、3°、5°時的質心位移演化曲線。其中，橫坐標為計算時間100 s，縱坐標為有無擋板時質心晃動幅度占最大值的百分比。結果表明，在液位恢復穩(wěn)定的過程中，質心做振幅衰減的周期性波動，并逐漸接近平衡。這是因為在微重力下，晃動時產生的湍流黏性耗散發(fā)揮作用，在分子黏性作用下動能轉化為熱能，液面逐漸靜止，趨于水平。

圖6 10-3g0下不同液面傾角下防晃效果對比曲線Fig.6 Comparison of anti-slosh effects under different liquid inclinations under 10-3g0

實際過程中，氣液界面晃動自初始狀態(tài)到完全穩(wěn)定耗時較長。為便于分析，本文以質心晃動幅值5×10-3cm（約占最大幅值的0.033%）作為達到穩(wěn)定的評判標準。取達到該標準的第一個波峰對應時間作為穩(wěn)定的時間。對比發(fā)現(xiàn)，當初始液位傾角較小時（θ=1°與θ=3°），箱內擋板對相界面晃動的抑制作用可忽略，而當θ=5°時，擋板在微重力下會加速液面從傾斜到水平的過渡速度，達到平衡的時間從77 s降到66 s，縮短了14.3%。此外，θ=1°工況下氣液界面的晃動處于較低水平，質心近似呈現(xiàn)等幅振蕩；θ=3°工況下晃動幅值可見明顯的衰減；而在θ=5°工況下晃動幅度大幅降低的同時，防晃隔板的作用凸顯。對比結果表明，箱內氣液相分布可在流體黏性耗散與防晃隔板兩種作用下達到新平衡。當初始氣液界面傾角較大時，晃動引起的慣性沖擊也相對較大，防晃隔板可有效抑制晃動幅值，加速其平復速率；而當初始傾角較小時，慣性沖擊力作用微弱，流體主要在黏性力作用下達到新的平衡，故防晃隔板作用微弱。

取液面初始傾角 5°，分別計算 10-3g0、10-4g0、10-5g0微重力工況下有擋板貯箱內的晃動特性，結果如圖7所示?？梢钥闯觯S著重力水平的降低，箱內流體晃動周期顯著延長，且所能達到的最大質心波動幅值也有所降低。重力越小，晃動周期越長，黏性耗散作用對晃動的抑制效果相應減小，故達到新的平衡位置需要更長時間。

圖7 不同重力環(huán)境下液位傾斜后質心演化曲線Fig.7 Evolution curve of centroid with liquid inclinations under different gravity environments

3.2 小過載側向激勵下的防晃特性

上面級運行時，燃料系統(tǒng)可能受到來自貯箱側向的激勵，對箱內的流體行為產生影響。為研究該激勵的特性，設激勵為作用于Y軸方向，頻率為0.144 Hz的徑向簡諧波。參考以往研究[24-26]，簡諧激勵的振幅一般取貯箱直徑的0.05～0.4倍。分別研究振幅為0.1 m、0.2 m、0.3 m的對應情況。10-3g0環(huán)境下，2.6 m、3.2 m液氧液位工況對應的質心演化規(guī)律如圖8（a）、（b）所示。

圖8 微重力下不同液位不同徑向正弦激勵防晃效果對比曲線Fig.8 Comparison of anti-slosh effects of different liquid levels and different lateral sinusoidal excitations under microgravity

由圖8可知：（1）激勵振幅越大，質心向上偏移的量越大，晃動越劇烈。圖中展示的時長約為一個周期，按設想質心曲線應呈往復趨勢，但由于液體質量很大，慣性也很大，受激勵后將立即向初始方向運動，不易向反方向運動，因此質心持續(xù)升高。而且激勵波的能量與振幅的平方成正比。振幅越大，波的能量越高，傳遞給推進劑的能量越高，轉化的位勢能越大，質心位置也越高。（2）比較圖8（a）、（b），發(fā)現(xiàn)擋板效果隨液位不同呈明顯差異。這是因為，當液位低于擋板時，液體的速度可分為兩部分，一部分平行于擋板表面，使其沿激勵作用方向運動，但不引起Z軸質心變化；另一部分垂直于擋板表面，使其做軸向運動。結合無擋板時液體質心運動趨勢，垂直方向的速度分量所占比例較大。在發(fā)生向上的位移過程中，流體運動被防晃板阻擋，故圖（a）中擋板效果顯著。當振幅為0.3 m時，可抑制質心偏移量達到2.6 cm；當液位高于擋板時，上部的液體不受擋板約束，運動隨意性較大，擋板對該部分流體幾乎不起作用，故圖（b）中抑制效果微弱，振幅為0.3 m時，質心偏移最大僅為0.3 cm。

為衡量擋板防晃效率，定義擋板效率為η，該值越大，表明擋板防晃效果越好。

擋板效率隨液位變化如圖9所示，該圖表明，當液位低于擋板最高位置時，擋板能發(fā)揮良好效果，效率較高。防晃效率與高度呈近似線性關系，重力越大，下降速率越快。10-3g0微重力下?lián)醢逍是€的降低速率比常重力更平緩，這意味著在半滿箱或更高液位時，擋板在軌階段仍能發(fā)揮一定的防晃作用。

圖9 不同重力環(huán)境下?lián)醢宸阑巫饔眯ЧS液位變化對比曲線Fig.9 Comparison curve of the baffle anti-slosh effect of different liquid levels under variablet gravity

圖10展示了圖8（b）計算周期內貯箱壓力及氣液相分布。貯箱內空白空間為氣枕區(qū)。液體顏色代表壓力水平，液體壓力由藍至紅遞增。由圖可知，在小過載下，貯箱內先出現(xiàn)壓力分層，壓力與液位深度成正比，底部液壓最大；之后，激勵傳來，氣枕沿著壁被傳送。不難推測，若以此趨勢繼續(xù)晃動，或振幅加大，氣枕將被傳送至底部輸液管口，造成輸液管含氣，再被輸至發(fā)動機將產生不利影響。因此在實際微重力工況下，應避免徑向激勵所產生的過量氣液位置偏移。

圖10 振幅為0.2 m時不同時刻貯箱壓力的分布對比圖Fig.10 Comparison of the distribution of tank pressure at different times when the amplitude is 0.2 m

3.3 空間慢旋下的防晃特性

為研究擋板對液體慢旋狀態(tài)的影響，將自旋角速度分別設置為0.1 rad/s、0.3 rad/s、0.5 rad/s，獲得10-3g0時液氧達到穩(wěn)態(tài)的氣液相界面分布，如圖11所示。隨著角速度的增大，液面內凹程度加深，壓力分布趨于均勻，中央氣體區(qū)呈近錐形，錐的高度逐漸增大，底圓直徑減小。這是由于氣液兩相間存在較大的密度差，旋轉時液相受慣性離心力更大，會被甩離軸心更遠，但受貯箱壁阻擋，表現(xiàn)為液位上升。當角速度大到一定程度，不難想象，氣錐將與底部相接。

圖11 貯箱軸向慢旋時的氣液相分布對比圖Fig.11 Comparison of gas and liquid phase distribution of tank during slow axial rotation

質心演變曲線如圖12所示。圖（a）中注明了貯箱內質心平衡位置。0.1 rad/s、0.3 rad/s工況的平衡位置相差0.65 cm，0.3 rad/s、0.5 rad/s工況的平衡位置相差1 cm。從達到平衡的時間上看，角速度越大，質心波動越大，達到平衡的時間越長。0.5 rad/s時擋板的效果對比如圖12（b）所示，從質心波動幅度判斷，有無擋板差異很小，即擋板使貯箱液體被抬高，質心位置更高。但擋板存在與否幾乎不影響流場達到穩(wěn)定建立平衡的時間。這是因為擋板主要對流體沿Z方向的位移起抑制作用，而該情況下流體受徑向的離心力為主，速度方向也與擋板作用平面平行，故防晃作用微弱。

圖12 流體質心演變曲線及防晃效果對比曲線Fig.12 Evolution curve of fluid centroid and comparison of anti-slosh effect

綜上，沿軸線的角速度有利于貯箱內實現(xiàn)快速氣液分離，而旋轉引起的尖細氣錐有可能接觸箱底，對氣體進入輸液管造成不利影響。因此有必要根據(jù)貯箱含氣率、分離完成要求等指標設計相應角速度。

3.4 防晃擋板對不同推進劑貯箱流體晃動的影響

上述結果基于液氧工質開展。而實際應用中，液氫/液氧往往成對出現(xiàn)。氫的密度、黏性力均比氧更小，因此，在地面激勵下，氫比氧需要更長的穩(wěn)定時間。在微重力下，重力的影響被削弱，為研究氫氧行為的差異性，用液氫替代液氧開展仿真研究，并比較兩種工質在不同激勵下的規(guī)律差異。

當液位傾斜時，3.1節(jié)結果已表明，工質自由振蕩具有明顯的周期性。本文Cassini貯箱內液體自由振蕩的周期由式（4）[19]計算：

式中各符號意義見1.1節(jié)，R為貯箱半徑。設θ=4°，計算氫氧工質的振蕩周期，為直觀觀察規(guī)律，取周期對數(shù)值為縱坐標，重力對數(shù)值為橫坐標，作圖如圖13所示，紅色虛線為氫氧振蕩趨勢線。

圖13 液位傾斜時氫、氧晃動周期對比曲線Fig.13 Comparison of hydrogen and oxygen slosh periods with liquid inclination

由圖可知，在本文考慮的重力范圍內，氫氧周期十分相近，但與理論計算值有一定差距。這是由于理論計算公式中未考慮擋板的存在，故當把理論公式應用于文中貯箱時，應加以修正。由圖中紅色虛線可得到重力范圍[10-5g0，10-1g0]內的簡化關系式如下：

經校驗，該式誤差不超過3.38%。

徑向激勵下流體類型對有擋板貯箱的影響如圖14所示，貯箱初始液位為3.2 m。圖中氫氧行為呈現(xiàn)較大差異。氫的振蕩幅度為6.5 cm，占最大允許值的43.3%，氧為4.0 cm，占最大允許值的16.0%。前面結果已表明，液位高于擋板時，液體晃動抑制效果較弱，結果表明，由于物性差異，氫的運動范圍遠大于氧。因此，對于徑向激勵工況，防晃效果與推進劑類型相關，應優(yōu)先關注液氫的防晃，如加高氫箱內的擋板高度等。

圖14 徑向正弦激勵下氫、氧質心演化對比曲線Fig.14 Comparison of the evolution of hydrogen and oxygen centroids under lateral sinusoidal excitation

在有擋板貯箱慢旋情況下，氫、氧兩流體行為如圖15所示。從質心演化曲線看，氫氧平衡位置幾乎等高，相較于貯箱尺寸，差異可忽略不計。因此對于慢旋工況，可不考慮氫氧箱內流體行為的差異。

圖15 慢旋下氫、氧質心演化對比曲線Fig.15 Comparison of the evolution of hydrogen and oxygen centroids under slow rotation

4 結論

本文采用CFD方法，就貯箱在空間中受多種激勵后的低溫推進劑流體運動特性開展了仿真研究，包括初始液位傾斜、徑正弦波、慢旋激勵等工況，并討論了防晃擋板的作用，結論如下：

（1）防晃擋板的作用效果受初始液位傾角大小影響，當傾角較小時，由于流體晃動的慣性沖擊較低，防晃擋板未見明顯的作用；而當初始傾角較高時，擋板可使相界面加速恢復穩(wěn)定。當過載從常重力瞬變?yōu)槲⒅亓?，傾斜液面達到動力平衡耗時約80 s。

（2）本文用有無擋板時流體質心變化量的比值來描述流體的晃動演化規(guī)律。貯箱經歷徑向激勵后，若液位處于擋板附近高度，則常、微重力下防晃擋板均能抑制晃動，且效果與液位近似成線性關系。液位越低，防晃效果越佳。同一激勵作用下，2.6 m液位工況時抑制的幅度約為2.6 cm，而3.2 m液位僅對應0.3 cm。

（3）箱內流體在徑向激勵作用后的晃動特性與其類型相關。同一激勵下，氫的質心晃動范圍可達6.5 cm，占最大允許值的43.3%，而氧的范圍為4.0 cm，僅占16.0%。對初始液位傾斜與慢旋激勵，流體類型對晃動演化的影響可忽略。

（4）微重力下，液體晃動幅度較大、慢旋速度較高均有可能造成貯箱底部暴露于氣枕，因此必須有效控制。應根據(jù)貯箱含氣率、分離完成要求等指標設計相應角速度。