高壓氫氣壓縮因子計算模型研究

王新跟 沈昱明

上海理工大學

氫氣燃燒發熱量高且綠色環保；氫能是當今世界最熱門的新型能源之一。隨著對氫能的開發和利用，對高壓氫氣的精確計量要求也越來越高。

氫氣的計量儀表包括科里奧利流量計、氣體渦輪流量計、氣體腰輪流量計、V錐流量計等。低溫或常溫下的高壓氫氣已不滿足理想氣體狀態方程。文獻[1]采用CFD(computational fluid dynamics)技術，對以高壓氫氣為介質的音速噴嘴的流出系數進行了系統研究。結果表明，氫氣的可壓縮因子和比熱比對音速噴嘴的流出系數影響顯著，且隨著噴嘴上游滯止壓力的增大，這種影響也越大。pVTt是氣體流量計量的原級基準，高壓氫氣的壓縮因子將直接影響pVTt的計量精確度。文獻[2]給出了包含壓縮因子的氫氣體積換算公式，其原理就是運用實際氣體狀態方程。可見，必須引入氫氣壓縮因子才能較精確計算氣罐中的氫氣質量。但目前，研究高壓氫氣壓縮因子的文獻尚不多見。

本研究對高壓氫氣壓縮因子的計算模型進行了討論，分別給出了Virial(維里)方程、Van der Waals方程(以下簡稱V-d-W方程)和經過修正的Redlich and Kwong方程(以下簡稱R-K方程)。并對美國國家標準技術研究所(The National Institute of Standards and Technology，以下簡稱NIST)給出的氫氣密度數據進行了回歸分析，且將NIST數據與各模型計算結果進行了比較。

1 氣體壓縮因子

1.1 氣體熱工狀態方程

確定氣體熱工狀態的主要參數有比容v、壓力p和溫度T，它們之間有確定的關系，稱為氣體狀態方程，即f(p,v,T)=0。對于理想氣體(或稱“完全氣體”)，氣體狀態方程遵循Clapeyron方程：

pv=RT

(1)

式中：p為氣體壓力，Pa；v為比容，m3/kg；T為溫度，K；R為氣體常數，J/(mol·K)。

在分子運動論中，所謂理想氣體是指：氣體分子幾何尺寸無限小，分子與分子之間沒有吸引力與排斥力。在壓力足夠低，密度足夠小的情況下，各種氣體都接近理想氣體性質，式(1)表示的Clapeyron理想氣體狀態方程具有一定的精確性。

1.2 實際氣體狀態方程

在氣體壓力較高、密度較大(溫度較低)時，氣體中分子運動規律復雜。本研究所討論的高壓氫氣實際條件(2～100 MPa、50～298 K)即屬于這種工況。此時，理想氣體狀態方程不能用來描述氫氣實際狀態下的p-v-T關系，而需要用較復雜的實際氣體狀態方程。一般而言，實際氣體狀態方程是在理想氣體狀態方程的基礎上，通過理論和實驗的方法，對理想氣體狀態方程進行修正取得的。引入壓縮因子Z=f(p,T)，對理想氣體狀態方程用壓縮因子進行修正是實際氣體狀態方程常用的表達方式。

1.3 壓縮因子

氣體壓縮因子Z表示實際氣體偏離理想氣體特性的程度，是對氫氣、天然氣等實際氣體計量中的一個重要物性參數。壓縮因子Z的定義為[3]：

(2)

從式(2)知，壓縮因子的實質是實際氣體或氣體在實際狀態下對理想氣體狀態方程(Z=1)的修正系數。

美國燃氣協會(AGA)提出了AGA8-92DC方程，采用氣體摩爾組成計算方法，給出了工作狀態下的天然氣壓縮因子方程[4-5]，但AGA8-92DC不支持計算氫氣的壓縮因子。

本研究歸納了幾種高壓氫氣壓縮因子的計算模型，并進行了分析討論。

2 計算模型

2.1 等容線性方程

等容線性方程是指在等容條件下，氫氣的p和T滿足線性關系：

p=A+KT

(3)

式中：A為常數項；K為系數。

式(3)中部分A、K值列于表1中。式(3)計算簡單，但誤差較大[6]。

表1 式(3)中的常數A和K值VM/(10-6 m3·mol-1)A/105 PaK/(105 Pa·K-1)26132.79.8228166.28.6130174.3 7.58①32169.26.6434157.15.8236142.95.1038126.74.45 注:①文獻[6]中該數據為7.08。

2.2 Virial方程

對于密度不大的氣體，可以采用Virial方程，即：

(4)

式中：RM為通用氣體常數，J/(mol·K)；VM為摩爾體積，L/mol；B、C為第二、第三Virial系數，可以由式(5)、式(6)確定[6]：

(5)

式中：B1=42.464，B2=-37.117 2，B3=-2.298 2，B4=-3.048 4，T0=109.781 K。

(6)

式中：c0=1.810 5×10-8m6/mol2(文獻[6]中該數據為1.310 5 ×10-8m6/mol2)，c1=2.148 6；T’0=20.615 K。

式(4)的計算結果如圖1所示。由圖1可以看到，氫氣密度隨著壓力的增大而增大，而隨著壓力的增大，Virial方程誤差也逐漸增大。將計算結果與NIST數據比較發現[7]：當T<200 K、p=10～100 MPa時，氫氣密度較大，Virial方程誤差較大，Virial方程已不再適用；當T>200 K，p=10～70 MPa時，最大誤差<0.9%。

2.3 V-d-W方程

考慮分子自身占有的容積和分子之間的相互作用力，Van der Waals對理想氣體狀態方程作了修正，提出了V-d-W方程，即：

(7)

(8)

(9)

式中：下標c表示臨界狀態，對于氫氣，其臨界參數pc=1.296 6 MPa，Tc=33.24 K，vc=0.032 2 m3/kg，Zc=0.304。

將式(7)寫成對比態形式：

(10)

式(10) 不包含任何表征物質特殊性的常數。根據熱力學相似性質，式(10)適用于一切符合V-d-W方程的物質，包括氫氣。根據V-d-W方程繪制的壓縮因子曲線圖如圖2所示。

對于pr<7、Z<1的應用場合，V-d-W方程與實驗得出的一些實際氣體，包括CH4、CmHn、CO2、N2等的Z與pr、Tr的關系式是一致的，如圖3所示。

但是對于p>10 MPa的高壓氫氣，pr>7，Z>>1，式(10)計算的壓縮因子與NIST數據比較，存在很大的偏差，如圖4所示。因此，不建議用V-d-W方程計算高壓氫氣壓縮因子。

2.4 R-K方程

R-K方程改進了V-d-W方程，因此，R-K方程可以獲得更精確的壓縮因子。R-K方程如式(11)所示：

(11)

(12)

(13)

b=b0-c0

(14)

(15)

式中：a0、b0、c0、b為與氣體臨界參數有關的常數；n為與物質有關的常數。

將氫氣的臨界參數代入上式，計算得到a0=609 9.48 J/mol，b0=9×10-3m3/kg，c0=3.5×10-6m3/kg。可見c0<

將式(10)寫成對比態形式：

(16)

對于天然氣，取n=0.47，可以獲得較高的計算精確度。例如，AGA8報告中給出的1#、4#、6#3組氣樣的計算結果與AGA8報告給出的實測值的誤差小于1.2%[8]。

對于高壓氫氣，本研究取n=0.31，計算結果如圖5所示。

計算結果顯示， R-K方程計算結果與NIST高壓氫氣壓縮因子數據有很好的一致性；在溫度為50～293 K、壓力為10～100 MPa時，二者最大誤差<5%；當溫度為200～ 293 K、壓力為10～60 MPa時，最大誤差<1%。

常溫、常壓下，R-K方程計算的氫氣壓縮因子與NIST數據也能很好吻合。計算結果顯示，當溫度為298 K、壓力為2～50 MPa時，誤差<0.4%，壓縮因子曲線如圖6所示。

3 NIST數據

將NIST給出的高壓氫氣的密度ρ根據式(2)換算成其壓縮因子Z。采用“平滑散點圖”方式繪出Z-p-T曲線圖，如圖7所示。

若以p/T作為橫坐標，則當p/T≥1 MPa/K時，Z-p/T可以擬合成一直線方程，如圖8所示，擬合的直線方程為：

(17)

式中：p為氫氣壓力，MPa；T為氫氣溫度，K。

式(17)計算結果顯示：當p/T=1～2 MPa/K時，線性誤差<0.55%；當p/T<1 MPa/K時，各溫度下的壓縮因子曲線嚴重發散(見圖9)，因此，該區間內不能用式(17)或其他形式的線性方程來計算壓縮因子。

4 結論

討論了幾種計算氣體壓縮因子的模型，特別是高壓氫氣的壓縮因子計算模型，其結論如下：

(1) 當T>200 K時，用Virial方程可以獲得較高的計算精確度。計算表明，當T>200 K，p=10～70 MPa時，最大誤差<0.9%。

(2) V-d-W方程和R-K方程都是根據熱力學對應態定律而導出的實際氣體狀態方程。在pr<7、Tr=1～2時，用V-d-W方程計算CH4、CmHn、CO2、N2、天然氣等實際氣體的壓縮因子都具有很好的一致性。

(3) R-K方程改進了V-d-W方程，可以獲得更精確的壓縮因子。取n=0.47，R-K方程計算天然氣壓縮系數的計算精度<1.2%；R-K方程計算高壓氫氣時，取n=0.31，在較寬的溫度、壓力范圍內，其計算精度<1%。

(4) 用NIST密度數據換算為壓縮因子數據。僅當p/T≥1 MPa/K時，Z-p/T可以擬合成一個直線方程(式(17))，線性誤差<0.55%；但當p/T<1 MPa/K時，由于線性誤差較大，因此不能擬合成任何形式的線性方程。