高 永，宇文雷，丁鳳林，李玉峰，黎 明，方忠堅

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

0 引言

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

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

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

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

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

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

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

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

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

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

1.2 氣體密度計算方法

當(dāng)氣體壓力較低時，可近似認(rèn)為其滿足理想氣體狀態(tài)方程；當(dāng)氣體壓力較高或者溫度較低時，氣體分子體積和氣體分子間作用力不可忽略，真實(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ù)庫來表達(dá)和查詢氣體密度與溫度和壓力的關(guān)系。如美國國家標(biāo)準(zhǔn)與技術(shù)研究院（NIST）數(shù)據(jù)庫即提供各類氣體的化學(xué)和物理性質(zhì)參數(shù)。該數(shù)據(jù)庫來自文獻(xiàn)數(shù)據(jù)的長期收集和積累，是一個非常實(shí)用的數(shù)據(jù)庫。圖2 給出NIST 數(shù)據(jù)庫中的氮?dú)饷芏?壓力-溫度關(guān)系曲面（253～333 K,0～30 MPa），密度的計算精度達(dá)到0.02%[13]。

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

其對應(yīng)的多項(xiàng)式系數(shù)見表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ú)饷芏扔嬎憬Y(jié)果，同時給出基于NIST 的顯式算法得到的結(jié)果（圖中以NIST 標(biāo)注）作為對比?？梢钥吹?，3 種真實(shí)氣體狀態(tài)方程中，RK 方程對氮?dú)饷芏鹊挠嬎憬Y(jié)果與NIST 結(jié)果最為接近。

圖3 不同計算模型的氮?dú)饷芏扔嬎憬Y(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 不同計算模型的氮?dú)饷芏扔嬎憬Y(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 氣瓶容積修正方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表4 對比分析了容積修正對基于NIST 的顯式算法計算結(jié)果的相對誤差影響?？梢钥吹?，容積修正能夠有效減小氣瓶內(nèi)氮?dú)馐Ｓ嗔坑嬎憬Y(jié)果的誤差，并且壓力越高，容積修正的效果越明顯。

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

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

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

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

3 結(jié)束語

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