高溫氫工質熱物理性質計算分析

房玉良，王成龍，田文喜，蘇光輝，秋穗正

(1.西安交通大學 能源與動力工程學院 動力工程多相流國家重點實驗室，陜西 西安 710049； 2.西安交通大學 核科學與技術學院 陜西省先進核能技術重點實驗室，陜西 西安 710049)

氫元素是最簡單、最輕的化學元素，也是宇宙中最多的元素，氫原子由1個質子和1個核外電子組成。作為一種綠色清潔能源，氫工質主要用作燃料、火箭推進劑、冷卻劑和化工原料等，在未來能源系統轉型、宇航推進動力中將具有廣闊的發展前景。

氫分子H2是一種分子量很小的雙原子分子，摩爾質量為2.015 88 g/mol。氫分子存在兩種自旋異構體，即兩個核自旋平行的正氫(Ortho-hydrogen，O-H2)和兩個核自旋反平行的仲氫(Para-hydrogen，P-H2)，兩者物理性質有所差異，化學性質相同，可相互轉化。熱平衡狀態下的氫氣稱為標準氫(Normal hydrogen，N-H2)，在室溫時，N-H2中的P-H2與O-H2之比為3∶1。N-H2的三相點為13.957 K，沸點為20.369 K，臨界壓力為1.296 4 MPa，臨界溫度為33.145 K，臨界密度為31.262 kg/m3，偏心因子為-0.219。

基于氫工質在航天動力推進、基礎物理學、等離子體物理學等方面的應用，研究人員開展了寬溫度、壓力范圍內的氫工質熱力學性質與輸運性質的基礎研究[1-9]。相關實驗研究和理論分析為物性計算程序開發積累了可靠數據、計算模型支持。早期程序[10-14]開發主要基于Fortran語言，采用數據表格插值或物性計算模型求解方法，通過輸入壓力、溫度(或比焓)獲取氫工質液態、氣態、熱解平衡態等條件下的密度、比焓(或溫度)、比熵、聲速、比熱容、比熱容比、黏度、導熱系數等相關熱物性參數。這些程序的計算范圍從三相點到104K、101～109Pa不等。美國國家標準局總結了20世紀80年代之前氫的熱物性研究，發布了兩本重要專著[15-16]，部分成果被當作科研與工業應用領域的標準數據庫。進入21世紀以來，氫的熱物性研究進入了精細化、精確化階段[17-18]，程序用戶界面更加友好、編程環境與接口更加多樣，使用方便快捷，如CoolProp、REFPROP等。

現有相關物性計算軟件主要存在年代久遠、模型與數據未及時更新，部分程序僅能計算分子氫物性等問題。本文基于最新的熱物性模型，開展了高溫氫工質熱物理性質計算模型研究和程序開發，計算求解密度、比焓、比熱容、聲速、黏度、導熱系數等熱物性參數，可為氫相關行業科研及應用提供借鑒支持。

1 熱物性計算模型

1.1 分子氫熱物性

1.1.1狀態方程 分子氫氣作為實際氣體考慮，狀態方程采用Aungier-Redlich-Kwong (ARK)模型[19]：

(1)

(2)

n=0.498 6+1.173 5ω+0.475 4ω2

(3)

(4)

(5)

(6)

式中：p為壓力，Pa；T為溫度，K；V為比容，m3·kg-1；Rg為氫氣氣體常數，Rg=4 124.5 J·kg-1·K-1；b、c為常數，m3·kg-1；α0、α為系數；ω為偏心因子，-0.219；下標cr表示臨界狀態。

1.1.2比焓 比焓h采用實際氣體狀態方程可推導如下：

h=h0(T)+pV-RgT-

(7)

理想氣體比焓h0(T)計算模型[2]如下：

(8)

式中，ai為系數。

1.1.3比定壓熱容 比定壓熱容cp采用實際氣體狀態方程可推導如下：

(9)

(10)

式中：ui、vi為系數；下標p表示定壓。

比定容熱容cV為：

(11)

式中，下標T、V分別表示定溫、定容。

1.1.4聲速 聲速Csound計算如下：

(12)

式中：ρ為密度，kg·m-3；下標s表示定熵。

1.1.5黏度 黏度μ計算采用Muzny模型[4]：

(13)

(14)

(15)

1.1.6導熱系數 導熱系數λ計算采用Assael模型[3]：

λ(ρ,T)=λo(T)+Δλe(ρ,T)+Δλcr(ρ,T)

(16)

(17)

(18)

(19)

式中，f1,i、f2,i、g1,i、g2,i、E1～E3為系數。

1.1.7擴散系數 擴散系數D計算模型[20]如下：

(20)

1.2 原子氫熱物性

原子氫是一種由氫原子組成的亞穩態物質，在高溫條件下視為理想氣體。氫原子的摩爾質量為1.007 947 g/mol，原子氫氣體熱力學行為簡單，只有空間3個方向上的平移，沒有旋轉和振動。

1.2.1狀態方程 狀態方程采用理想氣體模型：

pV=RgT

(21)

式中，Rg為原子氫氣體常數，Rg=8 248.9 J·kg-1·K-1。

1.2.2比焓 比焓計算模型[21]如下：

(22)

式中：k1～k5為系數；M為摩爾質量，M=1.007 947 g·mol-1；Tre為參考溫度，Tre=T/1 000，K。

1.2.3比定壓熱容 比定壓熱容計算模型[21]如下：

(23)

1.2.4聲速 聲速計算如下：

(24)

式中，κ為絕熱指數，κ=5/3。

1.2.5黏度 黏度計算模型[20]如下：

(25)

1.2.6導熱系數 導熱系數計算模型[20]如下：

(26)

1.2.7擴散系數 擴散系數計算模型[20]如下：

(27)

1.3 氫氣熱解及混合氫工質熱物性

1.3.1分子氫熱解 高溫狀態下，氫分子之間的共價鍵斷裂形成兩個氫原子，發生熱解或熱離反應。當溫度足夠高時，氣體還會發生電離，電子脫離原子核束縛產生等離子體。在高溫條件下，氫工質熱解產生單原子氫與分子氫會使得熱物理性質變化顯著[1]。常壓狀態下， H2的熱解溫度約為1 500 K，在溫度達到5 000 K以上時，H2幾乎全部熱解變成原子氫氣體，如圖1所示。分子氫熱解過程如下：

圖1 氫工質相圖Fig.1 Phase diagram of hydrogen

(28)

(29)

平衡常數只與溫度相關，與壓強、物質濃度無關。本文采用H2熱解平衡常數模型[5]為：

lgKp=-2.379 43×104/T+6.331 53

(30)

值得注意的是，T≤1 000 K時Kp?10-10，因此本研究認為T≤1 000 K時的Kp=0。

由上述研究可知，原子氫的份額是溫度和壓力的函數，該模型計算結果與文獻[5-6]數據對比如圖2所示，兩者相對偏差在±3%以內。

圖2 平衡態下原子氫的份額對比Fig.2 Comparison of atomic hydrogen fraction in equilibrium state

1.3.2密度 密度采用質量分數加權計算：

ρH-H2=wH·ρH(T,p)+(1-wH)·ρH2(T,p)

(31)

式中：wH為原子氫質量分數；下標H-H2表示混合氫工質。

1.3.3比焓、比熱容、聲速 比焓、比熱容、聲速加權計算方式[6]如下：

hH-H2=α·hH(T,p)+(1-α)·hH2(T,p)

(32)

MH-H2·cp,H-H2=α·MH·cp,H(T,p)+

(33)

(34)

ΔH(T)=2MH·hH(T)-MH·hH2(T,1 bar)

(35)

MH-H2=xHMH+(1-xH)MH2

(36)

式中：α為熱解度；ΔH為標準摩爾生成焓，即1 mol H2在標準狀態下(壓力為0.1 MPa)反應生成2 mol H產生的反應焓變；γ為絕熱指數，γH-H2=wHγH+(1-wH)γH2，γH=5/3，γH2=cp,H2/cV,H2；R為通用氣體常數。

1.3.4黏度 黏度計算采用文獻[7-8]中的模型：

(37)

(38)

(39)

(40)

1.3.5導熱系數 導熱系數計算采用文獻[7-8]中的模型：

λH-H2=λf+λr

(41)

(42)

(43)

(44)

(45)

(46)

2 程序開發

公開文獻調研到的相關氫工質熱物性參數計算程序列于表1，這些程序的計算范圍從三相點到104K、101～109Pa不等。

表1 氫工質熱物性計算程序Table 1 Thermophysical property code of hydrogen

基于上述計算模型，本研究采用MATLAB程序開發平臺編制了氫工質熱物性參數計算程序Prop_H_H2。Prop_H_H2可計算100～3 500 K、104～5×107Pa范圍內的H、H2及H-H2混合氫工質的密度、比焓、比熱容、聲速、黏度、導熱系數、擴散系數等熱物性參數。

2.1 程序組成與功能

Prop_H_H2為模塊化結構設計，可作為其他程序的子程序。Prop_H_H2由輸入模塊、H熱物性計算模塊、H2熱物性計算模塊、H-H2熱物性計算模塊、熱解模塊、輸出模塊組成，各熱物性計算模塊內設置相應的熱物性參數子函數。Prop_H_H2計算速度快，可實現單狀態點計算，也可實現曲線計算結果輸出，具有良好的人機交互繪圖、數據儲存、編程接口等特點。

2.2 程序計算流程

Prop_H_H2計算流程如圖3所示，根據物質種類(H、H2、H-H2)、所需計算的熱物性參數代碼，通過輸入溫度、壓力從而計算所選物質的相關物性參數。

圖3 Prop_H_H2程序框圖Fig.3 Block diagram of Prop_H_H2 program

3 計算結果與驗證

3.1 H2熱物性驗證

本文采用REFPROP程序計算結果進行H2物性驗證。REFPROP是由美國國家標準與技術研究院研制開發的工質熱力學和輸運性質計算軟件，該程序提供包括MATLAB在內的多種編程接口。REFPROP計算N-H2熱物性參數范圍為13.957～1 000 K、約2 000 MPa，但無法計算熱解平衡態的氫工質熱物性參數。本文通過REFPROP計算模型外推獲取1 000～3 500 K范圍內的N-H2熱物性參數。

Prop_H_H2與REFPROP熱物性參數計算結果對比如圖4所示，兩程序計算結果符合較好。熱物性參數相對偏差計算公式為：

圖4 H2熱物性對比Fig.4 Comparison of H2 thermophysical property

(47)

式中，X代表相應的熱物性參數。

圖5示出H2熱物性計算的相對偏差。結果表明：密度、比定壓熱容的最大相對偏差分別為6.80%、6.76%；黏度、導熱系數的最大相對偏差分別為13.05%、13.21%。出現最大相對偏差的狀態區域是在低溫、高壓區(100～200 K，1×107～5×107Pa)，此區域Prop_H_H2計算結果偏高。這是由于溫度越低、壓力越高，H2狀態越靠近液相，本文中采用的ARK模型在該區域適用性相對較差，計算出的密度偏高，進而導致其他熱物性在該區域處計算結果偏高。在200～3 500 K、1×104～2×107Pa范圍內，Prop_H_H2與REFRPROP熱物性參數計算結果的相對偏差在±1%以內。

圖5 H2熱物性計算偏差Fig.5 Deviation of H2 thermophysical property

3.2 H-H2熱物性驗證

有關熱解氫熱物性的研究主要集中在21世紀之前，相關研究資料與數據相對有限。本文根據搜集到的文獻數據[5-6，8-9，15，24]與Prop_H_H2計算結果進行驗證對比。這些文獻數據主要為20世紀90年代之前的數據，由美國國家標準與技術研究院、宇航局、洛斯阿拉莫斯國家實驗室等機構提供，部分熱物性參數與現有研究結果存在一定偏差。H-H2熱物性參數對比結果如圖6所示。由圖6可見，Prop_H_H2計算的熱解氫工質H-H2熱物性參數偏高，偏差較大的區域主要集中在程序計算范圍的邊界區域，如高溫、低壓區和低溫、高壓區。

圖6 H-H2熱物性對比Fig.6 Comparison of H-H2 thermophysical property

在計算范圍的大部分區域內，Prop_H_H2計算結果與文獻值相差較小，相對偏差在±5%之內。Prop_H_H2在計算比熱容、導熱系數的過程中，考慮了分子氫熱解過程中反應熱的影響。因此，在熱解反應程度較大，即原子氫摩爾份額占比較大時，比熱容、導熱系數隨溫度的變化有較大的波峰出現，圖中發現其峰值出現位置在原子氫熱解摩爾份額占比為50%附近。

綜上，Prop_H_H2計算H2、H、H-H2的密度、比定壓熱容、黏度、導熱系數等熱物性參數在100～3 500 K、104～5×107Pa范圍內是合理可靠的，計算值較其他程序、文獻的結果偏高，最大相對偏差為±15%；在200～3 000 K、104～107Pa范圍內，Prop_H_H2計算值會更加準確，相對偏差在±5%左右。

4 總結

本文開展了高溫氫工質熱力學與輸運性質研究，建立了原子態氫、分子態氫、熱解平衡態氫的熱物性計算模型，并基于MATLAB語言開發了熱物性計算程序Prop_H_H2。Prop_H_H2可適用于計算100～3 500 K、104～5×107Pa范圍內的H、H2及H-H2混合氫工質的密度、比用性相對較差，計算出的密度偏高，進而導致其他熱物性在該區域處計算結果偏高。在200～3 500 K、1×104～2×107Pa范圍內，Prop_H_H2與REFRPROP熱物性參數計算結果的相對偏差在±1%以內。

本程序可為氫工質相關的航天推進、物理學、能源動力等行業的科研和應用提供支持借鑒。