陳濤,劉國西,宋飛,武蔥蘢

北京控制工程研究所,北京100190

電推進系統(tǒng)Xe物理特性計算方法

陳濤＊,劉國西,宋飛,武蔥蘢

北京控制工程研究所,北京100190

空間電推進系統(tǒng)的工質(zhì)為Xe,其工作溫度范圍為-30～45℃,該范圍覆蓋Xe的臨界點。在臨界點附近,Xe可能呈現(xiàn)出多種形態(tài),且對溫度和壓力變化十分敏感,采用傳統(tǒng)狀態(tài)方程在上述范圍內(nèi)計算其物理參數(shù)偏差最大可達30%。為解決這一問題,以CH4作為參考流體,建立了一種基于對比態(tài)原理的Xe物理特性計算方法。該方法能夠?qū)Π庀唷⒁合?、超臨界區(qū)域的所有狀態(tài)Xe物理性質(zhì)進行準確計算。試驗研究與國外數(shù)據(jù)對比的結(jié)果表明,在整個壓力-溫度范圍內(nèi),計算誤差小于0.5%。

電推進;Xe;超臨界;對比態(tài)原理;物理特性;狀態(tài)方程

電推進技術(shù)是一種利用電能加速推進工質(zhì),從而實現(xiàn)高比沖的先進航天器推進技術(shù),由于其高比沖的優(yōu)勢,在先進國家的大型衛(wèi)星平臺及微小衛(wèi)星平臺上得到日漸廣泛的應(yīng)用[1-8]。

Xe作為電推進系統(tǒng)的工質(zhì),具有分子量大、熔點高、沸點高及臨界點高等特點?？紤]到系統(tǒng)的穩(wěn)態(tài)工作、充氣加壓等過程,其工作溫度范圍為-30～45℃,壓力范圍為1～30 MPa。整個工作范圍覆蓋氙的臨界點。在臨界點附近,Xe狀態(tài)對溫度和壓力變化十分敏感,分子間的作用力不可忽略,表現(xiàn)出強烈的真實流體特性。除此之外,在上述范圍內(nèi),Xe可能呈現(xiàn)出液態(tài)、氣態(tài)、氣液混合態(tài)以及超臨界態(tài)等多種形態(tài)。

目前,比較常用的描述真實流體性質(zhì)的狀態(tài)方程可以分為兩大類:解析形式和非解析形式。解析形式的狀態(tài)方程主要特點是當(dāng)溫度T、壓強p給定時,體積V可以解析地求出,即方程中V的次數(shù)不超過4次,例如Redlich-Kwong方程[9]、Peng-Robinson方程[10]、Soave-Redlich-Kwong(SRK)方程[9]等;非解析形式的狀態(tài)方程由于V次數(shù)較高,只能采用數(shù)值方法進行求解,例如Benedict-Webb-Rubin(BWR)/Modified Benedict-Webb-Rubin(MBWR)方程、Wagner方程等。

通過比較發(fā)現(xiàn),以上幾種傳統(tǒng)狀態(tài)方程無法同時準確計算各種形態(tài)下氙單質(zhì)的p-V-T關(guān)系,計算偏差最大可達30%,精度無法滿足應(yīng)用需求。而精確獲得上述關(guān)系對于Xe的充氣加壓、發(fā)射場加注及在軌監(jiān)測都極為重要[11-12]。

因此,深入了解Xe物理特性,對電推進系統(tǒng)設(shè)計具有重要意義[13]。目前,國內(nèi)外沒有公開發(fā)表的計算方法。為了解決這一問題,本文主要針對電推進系統(tǒng)工質(zhì)Xe的物理特性開展理論分析及計算方法研究。

1 理論模型及計算方法

典型物質(zhì)的三維熱力學(xué)曲面如圖1所示[9]。按照流體所處的不同狀態(tài)可以分為固相、液相、氣相三種不同相存在。在工程應(yīng)用中,常常需要對不同狀態(tài)下物質(zhì)的物理特性進行計算。

在大量的實際應(yīng)用中,研究人員提出了各種各樣的狀態(tài)方程來描述p-V-T之間的關(guān)系。狀態(tài)方程中最簡單、最常見的是理想氣體狀態(tài)方程,主要適用于氣體密度較低(即壓力不太高、溫度不太低)分子間相互作用可以忽略不計的情形。當(dāng)壓力較高或溫度接近飽和氣體溫度時,分子間的作用力不可忽略,為了能描述真實工質(zhì)的熱物理屬性,必須采用在低溫和高壓下都能對真實流體性質(zhì)進行刻畫的狀態(tài)方程。

由于傳統(tǒng)真實流體模型在高壓區(qū)、液態(tài)區(qū)和臨界點附近計算結(jié)果偏差較大,本研究中采用對比態(tài)原理結(jié)合BWR方程[14]來計算Xe的p-V-T關(guān)系。

1.1 對比態(tài)原理

對比態(tài)原理(Corresponding State Principle, CSP)是指:對所有物質(zhì)來說,其狀態(tài)變量經(jīng)過恰當(dāng)?shù)臒o量綱化處理后,都遵從普遍的變化規(guī)律。例如,對p、V、T都以臨界點參數(shù)進行歸一化處理后(即pr=p/pc,Tr=T/Tc,Vr= V/Vc),pr-Vr-Tr關(guān)系對于所有同類極性物質(zhì)都一樣。

對比態(tài)原理是目前應(yīng)用最廣泛的物性估算方法。也就是說,在獲得參考流體的物性試驗數(shù)據(jù)后,就可以采用對應(yīng)的對比態(tài)理論來計算目標流體的物性。對于具有不同分子結(jié)構(gòu)而顯現(xiàn)出不同極性的流體,按照對比態(tài)理論中使用參數(shù)的個數(shù),可以分為2參數(shù)CSP、3參數(shù)CSP和4參數(shù)CSP。不同的模型適用于不同的物質(zhì)類型。

由于Xe屬于單原子球形分子,本文選取同極性的具有球形分子結(jié)構(gòu)的CH4作為參考流體,采用2參數(shù)對比態(tài)原理結(jié)合BWR方程來計算Xe的物理特性。

2參數(shù)CSP是指,僅用兩個特征性質(zhì)來使狀態(tài)條件無量綱化(不需引入新的無量綱特征參數(shù)),其特征參數(shù)為Tc和pc或者Tc和Vc。只有單原子物質(zhì)如Ar、Kr、Xe或簡單流體才精確地遵循2參數(shù)對比態(tài)原理,其他物質(zhì)都有一些偏離。

1.2 BWR狀態(tài)方程

BWR方程[14]能同時對氣相、液相區(qū)域的物性進行準確的估算。方程將溫度多項式與密度的冪級數(shù)和指數(shù)相結(jié)合,其具體表達形式為

其中,15個系數(shù)An的具體形式如表1所示,除了參數(shù)r外,與溫度有關(guān)的函數(shù)共涉及32個擬合參數(shù)。

表1 BWR方程中系數(shù)項Table 1 Coefficient items in BWR equation

BWR方程中共有32個擬合參數(shù),文獻基于CH4大范圍的試驗數(shù)據(jù)給出了BWR模型中對應(yīng)的32個參數(shù)取值,見表2。文獻[15]指出,采用該組參數(shù)取值對CH4物性進行計算,最大誤差不超過0.45%。

1.3 氣液狀態(tài)辨識方法

由于計算過程中涉及氣、液、氣液混合等多種狀態(tài),為了對流體狀態(tài)進行辨識,本文基于美國國家標準技術(shù)研究院(NIST)數(shù)據(jù)庫中CH4的飽和點數(shù)據(jù),來獲得CH4飽和狀態(tài)下壓力和溫度的對應(yīng)關(guān)系。擬合得到的氣化線為

式中:p1=1.383×10-8,p2=-3.04×10-6,p3=0.0002654,p4=-0.01427,p5=0.4669。對于給定的溫度,可以得到該點下的飽和蒸氣壓。若流體壓力高于該數(shù)值,則流體處于液態(tài),反之處于氣態(tài)。

在此基礎(chǔ)上,Xe的氣液狀態(tài)辨識通過式(1)結(jié)合對比態(tài)理論實現(xiàn)。

1.4 模型求解流程

采用對比態(tài)原理對Xe密度進行求解時,首先根據(jù)Xe的溫度、壓力參數(shù),無量綱化處理后,采用對比態(tài)理論獲得參考流體CH4對應(yīng)的參數(shù),在此基礎(chǔ)上,反求BWR方程得到CH4的密度,再根據(jù)對比態(tài)理論求得Xe的狀態(tài)參數(shù)。

2 試驗系統(tǒng)及試驗方案

為了對理論模型的正確性進行驗證,同時獲得Xe第一手試驗數(shù)據(jù),開展Xe物理特性試驗研究。

2.1 系統(tǒng)組成

Xe物理特性試驗設(shè)備組成及工作原理如圖2所示,設(shè)備主要包括高純置換氣體、氣體凈化器、回收氣瓶、Xe運輸罐、稱重裝置、高低溫箱、真空泵、高壓氣瓶、高精度測量與控制系統(tǒng)。其中高壓氣瓶是本試驗中的主要試驗裝置,通過充入不同質(zhì)量的Xe實現(xiàn)不同密度條件下的參數(shù)測量;高精度測量系統(tǒng)用于監(jiān)測Xe運輸罐質(zhì)量,高低溫箱中Xe的質(zhì)量、記錄壓力容器內(nèi)溫度、壓力變化數(shù)據(jù)采集頻率不低于50 Hz;高低溫箱用于控制試驗中高壓氣瓶所需恒定溫度,恒溫精度不低于0.5℃。

試驗介質(zhì)選用高純Xe,純度≥99.995%。

圖2 Xe物理特性設(shè)備原理Fig.2 Schematic of the xenon physical property experiment equipment

2.2 試驗方法

首先對高壓氣瓶容積進行標定。獲得其精確容積之后,采用定容法測定Xe的p-V-T關(guān)系。在依次完成抽真空、置換、高純置換之后,對高壓氣瓶吸入一定量的高純度Xe,后啟動高低溫箱對高壓氣瓶進行升溫、降溫的過程,依次進行該密度下Xe壓力隨溫度的變化規(guī)律測量;之后對高壓氣瓶再充入定量Xe,再重復(fù)上述降溫回溫過程,實現(xiàn)另一密度條件下的數(shù)據(jù)測量。以此類推,依次完成0.2 kg/L、0.4 kg/L、0.6 kg/L、0.8 kg/L、1.0 kg/L、1.2 kg/L、1.4 kg/L和1.6 kg/L密度條件下的數(shù)據(jù)測量。

3 計算結(jié)果及試驗對比分析

根據(jù)BWR模型和對比態(tài)理論開展Xe物理特性計算,并與NIST數(shù)據(jù)[16]及試驗結(jié)果進行對比以驗證計算方法的正確性。

Xe的臨界點參數(shù)如表3所示。當(dāng)溫度高于Tc時,Xe只處于氣體狀態(tài)。當(dāng)溫度壓力都高于臨界點參數(shù)時,此時稱為超臨界狀態(tài)。這一狀態(tài)下,流體的密度表現(xiàn)得像液態(tài),但輸運特性表現(xiàn)得像氣態(tài)。

參考流體CH4的臨界點參數(shù)見表3。

表3 Xe及參考流體CH4的臨界點參數(shù)Table 3 Critical point parameters of Xe and CH4

3.1 飽和液線及飽和氣線

當(dāng)流體從液態(tài)變?yōu)闅鈶B(tài)時,會經(jīng)歷如下幾個狀態(tài):液態(tài)、飽和液體態(tài)、氣液混合態(tài)、飽和蒸氣態(tài)和氣態(tài)。采用對比態(tài)原理計算得到的Xe氣化線如圖3所示,計算溫度范圍為170～289.74 K,覆蓋從三相點到臨界點所有工況?？梢钥吹絏e在飽和態(tài)或者氣液混合態(tài)壓力隨溫度的變化關(guān)系。對應(yīng)的飽和氣線和飽和液線分別如圖4和圖5所示。結(jié)果表明,采用本文提出的方法與NIST數(shù)據(jù)符合較好,飽和線上密度計算最大誤差小于0.1%,能較好地分辨Xe給定壓力和溫度下的物理狀態(tài)。

3.2 熱力學(xué)曲線

為了驗證程序在大范圍內(nèi)計算結(jié)果的有效性,對橫跨液相、氣液混合、氣相區(qū)的溫度和壓力組合進行密度計算,來獲得Xe的熱力學(xué)曲線,并與NIST數(shù)據(jù)[16]進行對比。溫度計算范圍170～600 K,壓力計算范圍1～30 MPa,計算結(jié)果如圖6所示?？梢钥吹?計算結(jié)果與NIST數(shù)據(jù)完全吻合,能準確辨識流體狀態(tài)。證明計算結(jié)果準確可信,可以用于給定溫度T、壓力p條件下電推進系統(tǒng)工質(zhì)Xe的密度計算。

圖3 Xe的氣化線Fig.3 Evaporating line of xenon

圖4 飽和液線Fig.4 Saturation liquid line of xenon

圖5 飽和氣線Fig.5 Saturation gas line of xenon

圖6 1～30 MPa范圍內(nèi)Xe的T-ρ關(guān)系Fig.6T-ρrelationship in 1～30 MPa for Xe

對于工程應(yīng)用中的典型密度值,將計算數(shù)據(jù)(ρth為理論計算結(jié)果)匯總,并與試驗結(jié)果進行對比,如圖7所示,平均誤差小于0.5%?？梢钥闯?各密度條件下,計算結(jié)果與試驗數(shù)據(jù)符合良好,證明了計算方法可以完全滿足電推進系統(tǒng)設(shè)計需求。

圖7 計算數(shù)據(jù)與試驗數(shù)據(jù)對比Fig.7 Comparison between calculation and the experiment

由圖7中可以看到,在不同密度條件下,當(dāng)溫度逐步降低時,壓力都呈下降趨勢。同時,高壓氣瓶中Xe密度在臨界密度(ρc=1.1029 kg/L)附近和遠離臨界點時,兩者的變化趨勢存在顯著不同。例如,ρ=2.0 kg/L時,圖7中p-T關(guān)系存在明顯的轉(zhuǎn)折點,真實流體的特性表現(xiàn)明顯;而ρ=1.0 kg/L,壓力隨溫度變化平緩,曲線較為光滑。

結(jié)合Xe三相圖在密度壓力平面的投影,即圖8中的ρ-P圖可以發(fā)現(xiàn),在不同密度條件下,隨著溫度(或壓力)的降低,Xe經(jīng)歷的相變并不相同。當(dāng)密度大于臨界密度1.102 9 kg/L時,Xe從40℃降低到-30℃過程中,會經(jīng)歷從超臨界態(tài)變?yōu)橐簯B(tài),再到氣液混合態(tài)的情形;而當(dāng)密度小于1.102 9 kg/L時,會經(jīng)歷從超臨界態(tài)變?yōu)闅鈶B(tài),再到氣液混合態(tài)的情形。這是圖7中p-T關(guān)系在不同密度條件下變化趨勢存在區(qū)別的原因。

圖8 Xeρ-p關(guān)系Fig.8 Relationship betweenρ-pfor Xe

4 結(jié)束語

本文對電推進系統(tǒng)工質(zhì)Xe的物理特性進行了理論分析,建立了對應(yīng)的真實流體狀態(tài)方程,通過開展試驗研究、并國外NIST數(shù)據(jù)庫對比分析,結(jié)果表明:

1)建立的計算方法原理簡單,計算方便,基于現(xiàn)有數(shù)據(jù),無需提供額外參數(shù)。

2)在整個壓力-溫度范圍內(nèi),該方法能夠準確對氣相、液相、超臨界區(qū)域的Xe物理性質(zhì)進行計算,誤差小于0.5%,能夠完全滿足電推進貯供系統(tǒng)設(shè)計的應(yīng)用需求。

3)研究過程中采用試驗方法合理有效,能用于Xe高精度物性測量。

4)由于臨界溫度較高,在電推進系統(tǒng)的使用溫度范圍內(nèi),氙的真實流體的特性表現(xiàn)明顯,在不同密度條件下,p-T關(guān)系存在典型的差異。

References)

[1] ZHANG T P,WANG L,WU X M.Electric propulsion development for DFH-4SP satellite platform[C]∥Proceedings of 65th IAC,Toronto,Canada,IAC-14-C4.4.2,2014.

[2] VALERIE C T,JOSEPH M M,GBROWN J M.The dawn spacecraft[J].Space Sci.Rev.,2011,163(1): 175-249.

[3] ADAM S,ALEC D G,PETER Y P.Performance of a helicon Hall thruster operating with xenon,argon and nitrogen[J].Journal of Propulsion and Power, 2014,30(3):664-671.

[4] HANI K,WEN S H,THOMAS H.Overview of the development of the solar electric propulsion technology demonstration mission 12.5-kW Hall thruster[C]∥50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference, Cleveland,OH,AIAA 2014-3898,2014.

[5] DAN L,RAANAN E,GAL A.CAM200 Hall thrusterdevelopment overview[C]∥Proceedings of 66th IAC, Jerusalem,Israel,IAC-15-C4.4.4,2015.

[6] CHIEN K R.L-3 communications ETI electric propulsion overview[C]∥The 29th International Electric PropulsionConference,Princeton University,October 31-November 4,IEPC2005-315,2005.

[7] SEMENKIN A V.Overview of electric propulsion activity in Russia[C]∥The 30th International Electric Propulsion Conference,Florence,Italy,September 17-20,IEPC 2007-275,2007.

[8] KOPPEL C R,MARCHANDISE F,PRIOUL M.The SMART-1 electric propulsion subsystem around the Moon:in flight experience[C]∥41st AIAA/ASME/ SAE/ASEE Joint Propulsion Conference&Exhibit,Joint Propulsion Conference,Tucson,Arizona,July 10-13, AIAA 2005-3671,2005.

[9] BRUCE E P,JOHN M P,JOHN P O.Properties of gases and liquids[M].New York:McGraw-Hill, 2001:111-163.

[10] PENG D Y,ROBINSON B D.A new two constant equation of state[J].Industry Engineering Chemical Fundamentals,1976,15(1):59-64.

[11] JOHN W D,JOSEPH C,ANTHONY D.Advanced Xenon feed system(AXFS)development and hot-fire testing[C]∥45th AIAA/ASME/SAE/ASEE Joint Propulsion Conference&Exhibit,Denver,Colorado, August 2-5,2009,AIAA 2009-4910.

[12] GANAPATHI G B,ENGELBRECHT C S.Performance of the xenon feed system on Deep Space One[J].Journal of Spacecraft and Rockets,2000,37(3):392-298.

[13] OSBORN M F,NETWALL C J.Higher performance xenon flow system(XFS)optimized for low mass, volume and cost[C]∥45th AIAA Joint Propulsion Conference,Denver,Colorado,August 2-5,2009: AIAA 2009-4909.

[14] JACOBSEN R T,STEWART R B.Thermodynamic properties of nitrogen including liquid and vapor phases from 63 K to 2 000 K with pressure to 10000 bar[J]. Journal of Physical and Chemical Reference Data,1973, 2(4):757.

[15] ZONG N.Modeling and simulation of cryogenic fluid injection and mixing dynamics under supercritical conditions[D].Pennsylvania:The Pennsylvania State University,2005.

[16] NIST.NIST standard reference database[EB/OL].http:∥webbok.nist.gov.

(編輯:車曉玲)

A theoretical method for xenon's physical property used in electrical propulsion system

CHEN Tao*,LIU Guoxi,SONG Fei,WU Conglong
Beijing Institute of Control Engineering,Beijing 100190,China

The possible operation temperature of propellant xenon in electrical propulsion system is between-30℃and 45℃,which covers the xenon's critical temperature.At the critical temperature,xenon's state is sensitive to the operating pressure and the ambient temperature.The propellant may exhibit different forms.These characteristics make the traditional equation of state ineffective to calculate thep-V-Trelationship in the above temperature range(max error may rise to 30%).To solve this problem,a new calculation method based on the corresponding state principle was put forward for physical property estimation.The theoretical outcomes were compared with experiment data and available database.The results show that among the whole temperature and pressure range,the method is capable to predict the xenon's physical property in gas,liquid,supercritical state with an error less than 0.5%.

electrical propulsion;xenon;supercritical;corresponding state principle; physical properties;equation of state

V43

:A

10.3780/j.issn.1000-758X.2016.0018

2015-11-09;

:2015-12-25;錄用日期:2016-01-18;< class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間:2016-02-24 13:42:48

http:∥www.cnki.net/kcms/detail/11.1859.V.20160224.1342.015.html

*

:陳濤(1986-),男,博士,工程師,364993478@163.com,主要研究方向為電推進系統(tǒng)設(shè)計

陳濤,劉國西,宋飛,等.電推進系統(tǒng)Xe物理特性計算方法[J].中國空間科學(xué)技術(shù),2016,36(1):113-119.

CHEN T,LIU G X,SONG F,et al.A theoretical method for xenon′s physical property used in electrical propulsion system[J].Chinese Space Science and Technology,2016,36(1):113-119(in Chinese).

http:∥zgkj.cast.cn