高 永，宇文雷，丁鳳林，李玉峰，黎 明，方忠堅(jiān)

（1. 北京控制工程研究所； 2. 北京市高效能及綠色宇航推進(jìn)工程技術(shù)研究中心：北京 100190；3. 蘭州空間技術(shù)物理研究所，蘭州 730010； 4. 航天東方紅衛(wèi)星有限公司，北京 100094）

0 引言

冷氣推進(jìn)系統(tǒng)是指采用高壓氣體作為工質(zhì)，經(jīng)過(guò)系統(tǒng)減壓后，由噴管加速噴出產(chǎn)生推力的一類(lèi)推進(jìn)系統(tǒng)。其優(yōu)點(diǎn)是結(jié)構(gòu)簡(jiǎn)單、性能穩(wěn)定、安全可靠、成本低、功耗小、無(wú)毒無(wú)污染，缺點(diǎn)是比沖較低；隨著化學(xué)推進(jìn)和電推進(jìn)等高比沖推進(jìn)技術(shù)的發(fā)展，冷氣推進(jìn)系統(tǒng)逐漸被淘汰[1]。然而，隨著重力類(lèi)科學(xué)探測(cè)衛(wèi)星和小型化微納衛(wèi)星的快速發(fā)展，冷氣推進(jìn)系統(tǒng)重新展現(xiàn)出其獨(dú)有的優(yōu)勢(shì)，延伸出高精度冷氣推進(jìn)和微小型模塊化冷氣推進(jìn)兩個(gè)發(fā)展方向[2-3]。

面向重力類(lèi)科學(xué)探測(cè)衛(wèi)星應(yīng)用的高精度冷氣推進(jìn)系統(tǒng)具有推力精度高和系統(tǒng)質(zhì)心穩(wěn)兩大特點(diǎn)，其推力范圍涵蓋mN 至μN(yùn) 的量級(jí)，推力精度相比傳統(tǒng)冷氣推進(jìn)系統(tǒng)提高了3～6 個(gè)數(shù)量級(jí)[3-4]。但此類(lèi)衛(wèi)星的星上有效載荷對(duì)整星質(zhì)心變化敏感，因此其推進(jìn)系統(tǒng)的運(yùn)動(dòng)部件、氣體質(zhì)量均需要進(jìn)行精確的控制和計(jì)算。其中冷氣推進(jìn)系統(tǒng)工質(zhì)剩余量是推進(jìn)系統(tǒng)的關(guān)鍵參數(shù)之一，其計(jì)算精度直接關(guān)系到衛(wèi)星質(zhì)量的計(jì)算精度，間接影響衛(wèi)星科學(xué)數(shù)據(jù)的處理精度[5]。

冷氣推進(jìn)系統(tǒng)的工質(zhì)多選用氮?dú)?、氙氣等惰性氣體，相對(duì)于其他惰性氣體，氮?dú)獾挠行П葲_較高，故應(yīng)用較多；系統(tǒng)的初始?jí)毫σ话銥?5～30 MPa。工質(zhì)剩余量的計(jì)算一般基于氣體狀態(tài)方程展開(kāi)：鄢青青等[6]在姿控推進(jìn)劑加注量的精確計(jì)算中，采用理想氣體狀態(tài)方程對(duì)增壓氣體質(zhì)量進(jìn)行評(píng)估，應(yīng)用壓力不超過(guò)0.2 MPa；李強(qiáng)等[7]對(duì)在軌冷氣推進(jìn)系統(tǒng)的泄漏估計(jì)中，同樣基于理想氣體狀態(tài)方程對(duì)氣體質(zhì)量進(jìn)行計(jì)算，應(yīng)用壓力不超過(guò)13 MPa；Zou 等[8]比較不同狀態(tài)方程計(jì)算得到最高70 MPa 壓力下的氫氣密度，發(fā)現(xiàn)在高壓下不同狀態(tài)方程的計(jì)算結(jié)果出現(xiàn)明顯差別；張濤等[9]對(duì)超高壓冷氣推進(jìn)系統(tǒng)壓力模塊的性能開(kāi)展數(shù)值研究，采用維里方程處理實(shí)際氣體行為，應(yīng)用壓力最高120 MPa。

本文針對(duì)高壓氮?dú)饫錃馔七M(jìn)系統(tǒng)的工質(zhì)剩余量計(jì)算問(wèn)題，首先比較基于實(shí)際氣體狀態(tài)方程與基于數(shù)據(jù)庫(kù)查詢的氣體密度計(jì)算誤差，從中選取適合高壓氮?dú)饷芏扔?jì)算的方法；然后提出一種顯式計(jì)算方法，用于在判讀衛(wèi)星實(shí)時(shí)遙測(cè)數(shù)據(jù)期間快速獲取氣體工質(zhì)剩余量；此外，對(duì)冷氣推進(jìn)系統(tǒng)中高壓氣瓶在不同壓力下的容積變化進(jìn)行量化評(píng)估，從而得到一種經(jīng)過(guò)容積修正的高壓冷氣推進(jìn)系統(tǒng)剩余量精確計(jì)算方法；最后通過(guò)兩件氣瓶的地面加注量試驗(yàn)，對(duì)該計(jì)算方法的適用性進(jìn)行試驗(yàn)驗(yàn)證。

1 工質(zhì)剩余量計(jì)算方法

1.1 冷氣推進(jìn)系統(tǒng)簡(jiǎn)介

冷氣推進(jìn)系統(tǒng)主要由高壓氣瓶、氣加注閥、高壓壓力傳感器、自鎖閥、減壓裝置、低壓壓力傳感器、緩沖氣容、過(guò)濾器、推力器組件以及管路和連接件組成（參見(jiàn)圖1）。

圖1 典型冷氣推進(jìn)系統(tǒng)組成Fig. 1 Sketch of satellite cold gas propulsion system

冷氣推進(jìn)系統(tǒng)的加注一般在發(fā)射場(chǎng)進(jìn)行。加注前，通過(guò)地面管路將地面氣源與加注操作臺(tái)連接，再通過(guò)氣加注閥將氣源接入推進(jìn)系統(tǒng)的高壓氣瓶，組成一套密封系統(tǒng)；加注時(shí)，由加注操作臺(tái)控制流入高壓氣瓶的氣體流量，并實(shí)時(shí)監(jiān)視氣瓶的壓力和溫度；加注完成后，需將衛(wèi)星靜置10 h 以上，待氣瓶處于平衡狀態(tài)時(shí)測(cè)量得到氣瓶真實(shí)的溫度和壓力值。氣體加注量m的計(jì)算式為

其中：氣體密度ρ根據(jù)監(jiān)測(cè)所得氣瓶溫度和壓力值利用氣體狀態(tài)方程計(jì)算或由數(shù)據(jù)庫(kù)查詢得到；氣瓶容積V通過(guò)地面試驗(yàn)測(cè)得。

1.2 氣體密度計(jì)算方法

當(dāng)氣體壓力較低時(shí)，可近似認(rèn)為其滿足理想氣體狀態(tài)方程；當(dāng)氣體壓力較高或者溫度較低時(shí)，氣體分子體積和氣體分子間作用力不可忽略，真實(shí)氣體效應(yīng)明顯。目前，描述氣體壓力、溫度和密度的真實(shí)氣體狀態(tài)方程包括Van der Waals 方程、Redliche-Kwong（RK）方程和Peng-Robinson（PR）方程等。

理想氣體狀態(tài)方程為

式(2)～式(14)中：氣體密度ρ(kg/m3)可由氣體摩爾體積Vm(m3/mol)換算得到，ρ=M/Vm，M為氣體分子質(zhì)量，kg/mol；p為系統(tǒng)壓力，Pa；T為系統(tǒng)溫度，K；R為理想氣體常數(shù)，8.314 J/(mol?K)；Tc為氣體臨界溫度，K；Pc為氣體臨界壓力，Pa；ω為偏心系數(shù)。

除了采用氣體狀態(tài)方程求解外，工程上多采用曲線圖、表格或者數(shù)據(jù)庫(kù)來(lái)表達(dá)和查詢氣體密度與溫度和壓力的關(guān)系。如美國(guó)國(guó)家標(biāo)準(zhǔn)與技術(shù)研究院（NIST）數(shù)據(jù)庫(kù)即提供各類(lèi)氣體的化學(xué)和物理性質(zhì)參數(shù)。該數(shù)據(jù)庫(kù)來(lái)自文獻(xiàn)數(shù)據(jù)的長(zhǎng)期收集和積累，是一個(gè)非常實(shí)用的數(shù)據(jù)庫(kù)。圖2 給出NIST 數(shù)據(jù)庫(kù)中的氮?dú)饷芏?壓力-溫度關(guān)系曲面（253～333 K,0～30 MPa），密度的計(jì)算精度達(dá)到0.02%[13]。

圖2 NIST 數(shù)據(jù)庫(kù)中的氮?dú)饷芏?壓力-溫度關(guān)系曲面（253～333 K, 0～30 MPa）Fig. 2 Density-pressure-temperature profile of nitrogen(data from NIST database: 253～333 K, 0～30 MPa)

然而當(dāng)需要快速求解氣體密度，例如在對(duì)衛(wèi)星在軌遙測(cè)數(shù)據(jù)進(jìn)行實(shí)時(shí)判讀期間，往往需要有氣體密度與溫度和壓力的顯式函數(shù)，以快速準(zhǔn)確確定氣體剩余量。鑒于此，本文依據(jù)NIST 物質(zhì)物性平臺(tái)，查詢氮?dú)饷芏入S壓力和溫度變化的數(shù)據(jù)，采用雙變量多項(xiàng)式、利用數(shù)據(jù)處理軟件進(jìn)行擬合，得到氮?dú)饷芏扰c溫度和壓力的顯式關(guān)系式。該擬合公式為5 次多項(xiàng)式，

其對(duì)應(yīng)的多項(xiàng)式系數(shù)見(jiàn)表1，多項(xiàng)式擬合的均方根誤差（RSME）為0.037 29。

表1 氮?dú)饷芏?溫度-壓力關(guān)系式系數(shù)Table 1 Coefficients of density-temperature-pressure fitting expression for nitrogen

圖3 和圖4 給出理想氣體狀態(tài)方程和3 種真實(shí)氣體狀態(tài)方程在常溫（293 K）、0～70 MPa 壓力下和高壓（30 MPa）、200～500 K 溫度下的氮?dú)饷芏扔?jì)算結(jié)果，同時(shí)給出基于NIST 的顯式算法得到的結(jié)果（圖中以NIST 標(biāo)注）作為對(duì)比?？梢钥吹剑? 種真實(shí)氣體狀態(tài)方程中，RK 方程對(duì)氮?dú)饷芏鹊挠?jì)算結(jié)果與NIST 結(jié)果最為接近。

圖3 不同計(jì)算模型的氮?dú)饷芏扔?jì)算結(jié)果比較（293 K、0～70 MPa）Fig. 3 Comparison among calculated results of nitrogen density by using different gas equations of state(293 K, 0 to 70 MPa)

圖4 不同計(jì)算模型的氮?dú)饷芏扔?jì)算結(jié)果比較（30 MPa、200～500 K）Fig. 4 Comparison among calculated results of nitrogen density by using different gas equations of state(30 MPa, 200 to 500 K)

1.3 氣瓶容積修正方法

氣瓶在充入高壓氣體后，其容積會(huì)有一定的膨脹。在實(shí)際工程應(yīng)用中，會(huì)對(duì)氣瓶在最大工作壓力點(diǎn)的容積進(jìn)行測(cè)量。假設(shè)容積變形量與壓力呈線性關(guān)系，則有

式中：Pmax為最高工作壓力，Pa；Vmax為最高壓力下的氣瓶容積，m3；P0為大氣壓力，Pa；V0為大氣壓力下的氣瓶容積，m3。

綜上，冷氣推進(jìn)系統(tǒng)工質(zhì)剩余量計(jì)算公式可總結(jié)為

式中的氣體密度需根據(jù)氣體狀態(tài)方程迭代求解或者根據(jù)本文顯式計(jì)算方法得到。

2 地面加注測(cè)量結(jié)果與計(jì)算結(jié)果的比較分析

選取某批次48 L 復(fù)合材料氣瓶產(chǎn)品的質(zhì)心測(cè)量試驗(yàn)結(jié)果對(duì)氮?dú)赓|(zhì)量計(jì)算方法進(jìn)行驗(yàn)證。該氣瓶為薄壁金屬鈦內(nèi)襯碳纖維纏繞的復(fù)合材料球形高壓氣瓶，額定初始工作壓力30 MPa，在軌應(yīng)用期間工質(zhì)剩余量的計(jì)算精度要求優(yōu)于0.3 kg，本批次共2 件產(chǎn)品。氣瓶的質(zhì)心測(cè)量試驗(yàn)在全自動(dòng)高精度質(zhì)心測(cè)量臺(tái)上進(jìn)行[14]，如圖5 所示，氣瓶通過(guò)兩端固定安裝在質(zhì)心測(cè)量臺(tái)上——一端為固定出口端，一端為活動(dòng)圓柱端。氣瓶的出口端固定不動(dòng)、而圓柱端可沿軸向伸縮，因此充氣加壓后，氣瓶主要向圓柱端膨脹，徑向膨脹量均勻且為相對(duì)小量，軸向膨脹伸長(zhǎng)量可由預(yù)先安裝的激光位移計(jì)測(cè)得。質(zhì)心測(cè)量臺(tái)具備實(shí)時(shí)稱重功能，可直接獲得工質(zhì)的質(zhì)量。

圖5 氣瓶質(zhì)量與質(zhì)心測(cè)量裝置Fig. 5 Experiment setup of the tank mass and mass center measurement device

2.1 氣瓶容積變化實(shí)測(cè)結(jié)果

利用質(zhì)心測(cè)量臺(tái)的位移測(cè)試功能，分別測(cè)量了不同壓力下氣瓶的軸向膨脹伸長(zhǎng)量，結(jié)果如圖6 所示?？梢钥吹?，軸向膨脹伸長(zhǎng)量與氣瓶壓力呈線性關(guān)系，擬合決定系數(shù)（R2）分別為0.972 36 和0.996 51，線性度良好。這一結(jié)果也從側(cè)面證明1.3 節(jié)中氣瓶容積變形量與壓力呈線性關(guān)系的假設(shè)成立。

圖6 氣瓶軸向膨脹伸長(zhǎng)量的測(cè)量結(jié)果Fig. 6 Measured axial deformation length of test tanks

氣瓶的容積變形測(cè)量結(jié)果如表2 所示，其中V0和V30分別為大氣壓力和30 MPa 壓力下的氣瓶容積。

表2 氣瓶容積變形測(cè)量結(jié)果Table 2 Measured results of tank volume under differentpressures

2.2 工質(zhì)剩余量的實(shí)測(cè)結(jié)果

利用質(zhì)心測(cè)量臺(tái)的稱重功能，分別測(cè)量了不同壓力下氣瓶?jī)?nèi)的氣體質(zhì)量，并與工質(zhì)剩余量計(jì)算結(jié)果進(jìn)行對(duì)比。表3 為不同計(jì)算方程的氮?dú)馐Ｓ嗔坑?jì)算結(jié)果相對(duì)于實(shí)測(cè)值的平均誤差?？梢钥吹?，基于NIST 的顯式算法平均誤差最小，不超過(guò)0.05 kg；其次為RK 方程，而Van der Waals 方程的計(jì)算結(jié)果偏差最大。

表3 不同計(jì)算方程的氮?dú)馐Ｓ嗔坑?jì)算結(jié)果的平均誤差（含容積修正）Table 3 Mean error of nitrogen residual mass values calculated by using different models

圖7 給出2 件氣瓶加注氮?dú)赓|(zhì)量試驗(yàn)的實(shí)測(cè)結(jié)果與不同計(jì)算方程的計(jì)算結(jié)果對(duì)比。從圖中可以看出，基于NIST 的顯式算法所得結(jié)果在數(shù)據(jù)絕對(duì)值和變化趨勢(shì)上都與地面實(shí)測(cè)結(jié)果最為接近。

圖7 氣瓶加注氮?dú)赓|(zhì)量的實(shí)測(cè)結(jié)果與不同計(jì)算方程的計(jì)算結(jié)果（含容積修正）對(duì)比Fig. 7 Comparison between experimental fuel mass values and calculated results with volume correction

表4 對(duì)比分析了容積修正對(duì)基于NIST 的顯式算法計(jì)算結(jié)果的相對(duì)誤差影響。可以看到，容積修正能夠有效減小氣瓶?jī)?nèi)氮?dú)馐Ｓ嗔坑?jì)算結(jié)果的誤差，并且壓力越高，容積修正的效果越明顯。

表4 容積修正對(duì)基于NIST 的顯式算法計(jì)算結(jié)果的影響Table 4 Comparison between calculated nitrogen mass errors for NIST models with and without volume correction

2.3 工質(zhì)剩余量計(jì)算誤差分析

如前所述，衛(wèi)星冷氣推進(jìn)系統(tǒng)工質(zhì)剩余量計(jì)算需依據(jù)3 個(gè)輸入條件——?dú)馄繅毫Α馄繙囟群蜌馄咳莘e。設(shè)氣瓶壓力的不確定度為u(p)、氣瓶溫度的不確定度為u(T)、氣瓶容積的不確定度為u(V)，則有工質(zhì)剩余量計(jì)算的相對(duì)不確定度u(m)[15]為

某衛(wèi)星冷氣推進(jìn)系統(tǒng)配置了2 件48 L 氣瓶，加注目標(biāo)壓力為30 MPa（20 ℃）。星上壓力傳感器采集壓力值的不確定度u(p)為0.175 MPa，氣瓶溫度測(cè)量的不確定度u(T)為0.5 ℃，氣瓶容積測(cè)量的不確定度u(V)為0.1 L，根據(jù)公式(18)，該衛(wèi)星工質(zhì)剩余量計(jì)算的不確定度為0.18 kg，相對(duì)不確定度為0.62%。

3 結(jié)束語(yǔ)

本文提出一種衛(wèi)星推進(jìn)系統(tǒng)氣體工質(zhì)剩余量的顯式計(jì)算方法，基于NIST 物質(zhì)物性平臺(tái)，給出氮?dú)饷芏?溫度-壓力的顯式擬合公式；并采用容積線性修正模型減小氣體工質(zhì)剩余量計(jì)算的平均誤差。用地面加注量試驗(yàn)數(shù)據(jù)對(duì)計(jì)算模型進(jìn)行驗(yàn)證，計(jì)算結(jié)果與測(cè)量結(jié)果間的誤差總體上不超過(guò)0.05 kg，計(jì)算方法的不確定度為0.18 kg，滿足工質(zhì)剩余量的計(jì)算精度優(yōu)于0.3 kg 的要求。