羅天培，孫 德，張家仙，李 茂

（1.北京航天試驗技術(shù)研究所，北京，100074；2.北京市航天試驗技術(shù)與裝備工程技術(shù)研究中心，北京，100074）

0 引 言

在液體火箭動力系統(tǒng)、液體火箭發(fā)動機試驗系統(tǒng)的運行過程中，推進劑貯箱以及增壓氣瓶等的充、放氣是常見的操作之一。尤其對于地面試驗臺來說，在研制階段為了考核發(fā)動機在不同的泵前壓力或者燃燒壓力下的工作性能等，充放氣過程會更為頻繁。同時為了滿足發(fā)動機、動力系統(tǒng)工作時序以及安全使用等的要求，貯箱或氣瓶的充壓和泄壓需要在嚴格的規(guī)定時間內(nèi)完成。以某型號火箭動力系統(tǒng)試驗為例，箭上氣瓶的充氣過程需在點火前-9 h到-3.5 h內(nèi)完成，考慮到充氣過程溫升帶來的強度、疲勞等問題，該時間不僅不能延長，同時也不能顯著縮短；點火前推進劑貯箱的預(yù)增壓需在-4 min到-3 min之內(nèi)完成，快速增壓一方面保證了發(fā)動機泵前壓力入口壓力的需求，另一方面也使得推進劑溫度基本保持不變，溫度品質(zhì)好。因此準確計算容器的充、放氣過程對于火箭動力系統(tǒng)的時序安排以及氣瓶的強度設(shè)計等方面具有重要意義。

容器的充、放氣是多變指數(shù)和傳熱系數(shù)變化的復(fù)雜過程[1]，同時，受諸如氣管長度、通徑、內(nèi)壁粗糙度等很多因素的影響[2～4]，其中高壓容器的充、放氣過程相對來說更為復(fù)雜，以放氣為例：放氣初始階段由于放氣口前后壓比一般遠超臨界壓力比，放氣口處流動會達到壅塞狀態(tài)，流動速度為當?shù)芈曀?，放氣過程導(dǎo)致的容器壓力降低會使得流量逐漸下降；當降到臨界壓力比以下時，流動變?yōu)閬喡曀贍顟B(tài)，流量進一步降低。在放氣的同時，氣體的膨脹過程還會使得容器內(nèi)溫度下降。

許多學(xué)者曾針對這一過程開展過大量的研究：在充氣方面，Wilson[5]等研究了恒壓氣源向剛性容器充氣的一般過程；安剛[6]基于氫燃料電池汽車的儲氫需求，對高壓氫氣瓶充氣時的溫升效應(yīng)開展了分析；McMurray 等[7]對橢圓形截面氣瓶充氣過程中氣瓶的變形和破裂進行了分析。在放氣方面，Bu 等[8]建立了天然氣貯罐放氣模型，研究了放氣過程中容器內(nèi)溫度的變化以及流量的變化；Todd 等[9]研究了汽車油箱的放氣過程，尤其對于該過程里油箱內(nèi)可能發(fā)生的相變進行了探索；Wang 等[10]研究了放氣過程對容器本身結(jié)構(gòu)特性的影響，針對某次放氣過程導(dǎo)致的管道破裂開展了分析；Marek 等[11]還曾基于CFD 仿真技術(shù)對非牛頓流體從貯箱中的泄出過程開展過研究。

本文為了尋求一種可準確、高效的計算容器充、放氣過程的方法，進而滿足動力系統(tǒng)工作時序設(shè)計的要求，基于開口系能量方程、理想氣體狀態(tài)方程以及不同狀態(tài)下的質(zhì)量流量方程等，利用集總參數(shù)法編制了計算程序，并針對常溫容器及低溫貯箱的充、放氣過程開展了計算模擬，通過與試驗數(shù)據(jù)的對比，驗證了程序的合理性，并開展了相應(yīng)的分析。

1 數(shù)學(xué)模型及計算方法

1.1 假設(shè)條件

本文數(shù)學(xué)模型的建立基于如下假設(shè)條件：

a）容器內(nèi)氣體符合理想氣體狀態(tài)方程；

b）氣體的定壓比熱容Cp及定容比熱容CV為常數(shù)，不隨溫度發(fā)生變化；

c）由于充、放氣過程相對迅速，整個過程視為絕熱過程，不考慮氣體與容器的熱交換以及容器與外界環(huán)境的熱交換；

d）容器內(nèi)的壓力、溫度和密度是均勻分布的，不考慮溫度分層等現(xiàn)象；

e）將排氣口橫截面由流動邊界層等因素導(dǎo)致的壓力、溫度和密度的不均勻統(tǒng)一用一個流量系數(shù)Cd來折算；

f）忽略氣體重力的影響。

1.2 數(shù)學(xué)模型建立

將容器及排氣口整體建立開口系能量方程，充、放氣過程中容器內(nèi)部內(nèi)能的變化量等于流入、流出氣體焓的差值，即：

式中min為進入容器內(nèi)的氣體質(zhì)量；mout為流出容器的氣體質(zhì)量；hin為流入氣體的焓值；hout為流出氣體的焓值；ECV為容器內(nèi)氣體整體的內(nèi)能；dECV為容器內(nèi)部內(nèi)能的變化量，計算方法如下：

式中u為氣體的內(nèi)能；TCV為容器內(nèi)氣體的溫度。

容器內(nèi)氣體的壓力、密度、溫度符合理想氣體狀態(tài)方程：

式中Rg為氣體常數(shù)。

氣體在排氣口的流動狀態(tài)由臨界壓力比βcr確定，由下式計算：

式中γ為理想氣體比熱比，γ=1.4；求得βcr的大小為0.528。

當排氣口下游與上游的壓力比小于或等于βcr時流動達到壅塞狀態(tài)，此時流動的質(zhì)量流量qm為

式中Cd為流量系數(shù)，一般取值范圍為0.5～0.9；下標us 表示上游；q(λ)在壅塞狀態(tài)時為1；A為流出口的面積；K對于特定的氣體來說為常數(shù)，K的計算方法為

當排氣口下游與上游的壓力比大于βcr時流動為亞聲速狀態(tài)，此時流動的質(zhì)量流量qm為

式中 下標ds 表示下游。

本文將以上各式通過C++語言編程聯(lián)立求解，其中時間步長Δt經(jīng)敏感性測試，取為0.1 ms，利用式（7）或式（9）計算得出的qm與Δt相乘即可得到式（1）與式（2）中的δmin和δmout，各式在一個時間步長內(nèi)逐一更新后進入下一個步長內(nèi)運算，直至達到預(yù)設(shè)時間后跳出循環(huán)，計算結(jié)束。

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

2.1 常溫容器充、放氣過程計算

如圖1 所示，試驗系統(tǒng)由高壓室、開關(guān)閥及連接管道以及低壓室組成，其中低壓室上有一個排氣口，試驗介質(zhì)為氮氣。試驗開始后迅速打開開關(guān)閥，高壓室會向低壓室充氣，同時低壓室壓力升高后會從排氣口放氣，各設(shè)備的參數(shù)如表1 所示，試驗對低壓室的壓力變化進行了采集。

圖1 常溫容器充放氣試驗原理Fig.1 Schematic of Vessel Charge and Discharge Test at Normal Temperature

表1 常溫容器充放氣試驗參數(shù)Tab.1 Parameters of Vessel Charge and Discharge Test at Normal Temperature

圖2 給出了低壓室壓力變化的結(jié)果，圖中test 1 及test 2 為2 次試驗測得的結(jié)果，simulation 為計算所得到的結(jié)果。

圖2 低壓室壓力隨時間的變化Fig.2 Pressure Variations in Low Pressure Chamber Against Time

由圖2 可知，在0.5 s 左右低壓室內(nèi)壓力達到了其峰值1.4 MPa 左右，試驗和計算吻合度較高，由于低壓室本身同時包含了充氣及放氣的過程，兩者符合較好可充分證明計算模型的合理性和準確度。

圖3 中兩條藍線分別代表了高壓室（High Pressure Chamber）、低壓室（Low Pressure Chamber）壓力隨時間的變化。兩條紅線分別代表了高壓室向低壓室排氣的流量以及低壓室向外界大氣環(huán)境排氣的流量，為了便于顯示，流量坐標軸采用對數(shù)分布形式表示。由圖3可知，在t=0 時刻，高壓室初始時刻壓力為15 MPa，初始排放流量達到了約300 kg/s，低壓室初始壓力為0，與外界壓力平衡，因此初始排放流量也為0。隨著閥門開啟，高壓室迅速向低壓室充氣，同時低壓室向外界排氣的流量從0 開始增加，低壓室的壓力以及向外界排氣的流量在0.5 s 左右達到峰值，隨后開始減小；在0.6 s 左右高壓室向低壓室排氣的流量曲線與低壓室向外界排氣的流量曲線相交，此時低壓室內(nèi)氣體質(zhì)量達到峰值，隨后向外界環(huán)境排氣流量大于由高壓室補充進來的流量；在1 s 左右高壓室、低壓室壓力平衡，該平衡狀態(tài)一直維持到4.5 s 左右，即兩者壓力與外界環(huán)境壓力平衡后排氣過程結(jié)束。

圖3 高、低壓室壓力以及排氣流量隨時間的變化Fig.3 Variations of Pressure and Discharge Flow Rate in High and Low Pressure Chambers Against Time

2.2 低溫貯箱放氣過程計算

低溫貯箱放氣模型如圖4 所示，貯箱內(nèi)承裝液體為液氧，增壓氣體為氮氣，增壓前貯箱通過敞口靜置過程將氣枕壓力排至表壓接近0，因此增壓后相比于氮氣，氣枕內(nèi)氧氣濃度很小，計算時予以忽略，各物性參數(shù)全部選為氮氣的參數(shù)；另外，取氣枕區(qū)域為計算對象，忽略液體對排氣過程的影響，具體計算條件如表2 所示，初始氣枕溫度取180 K。

表2 低溫容器放氣計算參數(shù)Tab.2 Computation Parameters of Cryogenic Tank Discharge

圖4 低溫貯箱放氣模型Fig.4 Schematic of Cryogenic Tank Discharge

圖5 給出低溫貯箱放氣過程試驗和計算的對比結(jié)果，由圖5 可知，兩者基本趨勢一致，但在不同時間段內(nèi)，兩者放氣速度不完全一致：在約0～1 s 內(nèi)，計算放氣速度高于試驗結(jié)果，在約1～5 s 內(nèi)，計算放氣速度低于試驗結(jié)果，5 s 之后放氣速度再次高于試驗結(jié)果，并隨著時間推移，兩者差別逐漸減小。

圖5 低溫貯箱放氣過程試驗和計算結(jié)果的對比Fig.5 Comparison of Test and Computation of Cryogenic Tank Discharge Process

造成這一區(qū)別的主要可能原因分析如下：

a）初始階段（0～1 s）：由于貯箱壓力顯著高于外界壓力，放氣初始階段流動處于壅塞狀態(tài)，放氣流量計算滿足式（7），此式中Tus計算程序選用為氣枕平均溫度，由于實際貯箱里氣枕溫度會有分層現(xiàn)象，從上到下溫度分布由高到低，因此實際排氣口溫度會高于平均溫度，即計算程序里給的溫度偏低，導(dǎo)致排氣流量高于實際情況，即放氣速度高于試驗結(jié)果。

b）中間階段（1～5 s）：隨著初始放氣階段的進行，放氣口溫度逐漸接近容器內(nèi)平均溫度，放氣流量的計算誤差逐漸減小，此時放氣過程帶來的降溫效應(yīng)起了主導(dǎo)作用。對于氣枕內(nèi)主要的成分氮氣來說，貯罐內(nèi)初始壓力為2.7 MPa，此壓力對應(yīng)的飽和蒸氣溫度約為120.3 K，而氣液交界面溫度與液氧溫度基本一致（90 K），這表明在增壓過程中在氣液交界面附近勢必有一部分氮氣會發(fā)生相變形成液氮，放氣帶來的降溫效應(yīng)會促進這一過程，相當于容器氣枕內(nèi)一部分質(zhì)量減少了，導(dǎo)致了試驗實際放氣速度“高于”計算速度。

c）完成階段（5 s 之后）：此階段溶解在液氧中的氮氣逸出對容器壓力的變化起了主導(dǎo)作用，氣體在液體中的溶解受壓力、溫度的影響，在溫度一定的情況下，壓力越小溶解度越低，對于容器內(nèi)的液氧來說，放氣導(dǎo)致的降溫主要影響氣枕溫度，對液氧溫度影響很小，而隨著容器內(nèi)的壓力降低，會有氮氣不斷從液氧中逸出，而放氣過程帶來的擾動會不斷加劇這一過程，相當于氣枕內(nèi)的氣體質(zhì)量增壓了，反映到放氣曲線上即計算的放氣速度高于實際試驗放氣速度。

d）假設(shè)條件與實際的偏差：本文計算程序采用理想氣體模型，該模型常溫、常壓下較為合理，在高壓、低溫等條件下會與實際存在一定偏差；同時，計算中設(shè)定氣體定壓比熱容Cp及定容比熱容CV為常數(shù)，而實際兩者會隨溫度發(fā)生變化，這些假設(shè)本身會造成計算的誤差。

另外，需要指出的是，分析中定義的初始、中間、完成各階段發(fā)生的物理過程并不是在該時段內(nèi)獨立發(fā)生，而是對應(yīng)時段內(nèi)該過程起了主導(dǎo)作用。

總體來說，計算誤差不大，采用該模型可以計算低溫容器的放氣過程。

3 結(jié) 論

本文基于開口系能量方程、理想氣體狀態(tài)方程以及不同流動狀態(tài)下的質(zhì)量流量方程等，利用集總參數(shù)法編制了計算程序，對常溫容器的充、放氣過程以及低溫容器的放氣過程進行了計算，并與試驗數(shù)據(jù)進行了對比分析，結(jié)果表明：

a）該方法對于常溫容器的充、放氣過程可以準確模擬；

b）由于未考慮貯箱氣枕內(nèi)的熱分層現(xiàn)象，放氣初始階段計算質(zhì)量流量偏大；

c）放氣過程所帶來的降溫效應(yīng)會促進氣枕內(nèi)氮氣的液化，同時氣枕內(nèi)壓力的降低會減少氮氣在液氧中的溶解，兩者對氣枕內(nèi)總氣體質(zhì)量產(chǎn)生了影響，進而造成了計算和試驗貯箱壓力變化的偏差。