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

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

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

氫分子H2是一種分子量很小的雙原子分子，摩爾質(zhì)量為2.015 88 g/mol。氫分子存在兩種自旋異構(gòu)體，即兩個核自旋平行的正氫(Ortho-hydrogen，O-H2)和兩個核自旋反平行的仲氫(Para-hydrogen，P-H2)，兩者物理性質(zhì)有所差異，化學性質(zhì)相同，可相互轉(zhuǎn)化。熱平衡狀態(tài)下的氫氣稱為標準氫(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。

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

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

1 熱物性計算模型

1.1 分子氫熱物性

1.1.1狀態(tài)方程 分子氫氣作為實際氣體考慮，狀態(tài)方程采用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為氫氣氣體常數(shù)，Rg=4 124.5 J·kg-1·K-1；b、c為常數(shù)，m3·kg-1；α0、α為系數(shù)；ω為偏心因子，-0.219；下標cr表示臨界狀態(tài)。

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

h=h0(T)+pV-RgT-

(7)

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

(8)

式中，ai為系數(shù)。

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

(9)

(10)

式中：ui、vi為系數(shù)；下標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導熱系數(shù) 導熱系數(shù)λ計算采用Assael模型[3]：

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

(16)

(17)

(18)

(19)

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

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

(20)

1.2 原子氫熱物性

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

1.2.1狀態(tài)方程 狀態(tài)方程采用理想氣體模型：

pV=RgT

(21)

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

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

(22)

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

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

(23)

1.2.4聲速 聲速計算如下：

(24)

式中，κ為絕熱指數(shù)，κ=5/3。

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

(25)

1.2.6導熱系數(shù) 導熱系數(shù)計算模型[20]如下：

(26)

1.2.7擴散系數(shù) 擴散系數(shù)計算模型[20]如下：

(27)

1.3 氫氣熱解及混合氫工質(zhì)熱物性

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

圖1 氫工質(zhì)相圖Fig.1 Phase diagram of hydrogen

(28)

(29)

平衡常數(shù)只與溫度相關(guān)，與壓強、物質(zhì)濃度無關(guān)。本文采用H2熱解平衡常數(shù)模型[5]為：

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

(30)

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

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

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

1.3.2密度 密度采用質(zhì)量分數(shù)加權(quán)計算：

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

(31)

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

1.3.3比焓、比熱容、聲速 比焓、比熱容、聲速加權(quán)計算方式[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在標準狀態(tài)下(壓力為0.1 MPa)反應(yīng)生成2 mol H產(chǎn)生的反應(yīng)焓變；γ為絕熱指數(shù)，γH-H2=wHγH+(1-wH)γH2，γH=5/3，γH2=cp,H2/cV,H2；R為通用氣體常數(shù)。

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

(37)

(38)

(39)

(40)

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

λH-H2=λf+λr

(41)

(42)

(43)

(44)

(45)

(46)

2 程序開發(fā)

公開文獻調(diào)研到的相關(guān)氫工質(zhì)熱物性參數(shù)計算程序列于表1，這些程序的計算范圍從三相點到104K、101～109Pa不等。

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

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

2.1 程序組成與功能

Prop_H_H2為模塊化結(jié)構(gòu)設(shè)計，可作為其他程序的子程序。Prop_H_H2由輸入模塊、H熱物性計算模塊、H2熱物性計算模塊、H-H2熱物性計算模塊、熱解模塊、輸出模塊組成，各熱物性計算模塊內(nèi)設(shè)置相應(yīng)的熱物性參數(shù)子函數(shù)。Prop_H_H2計算速度快，可實現(xiàn)單狀態(tài)點計算，也可實現(xiàn)曲線計算結(jié)果輸出，具有良好的人機交互繪圖、數(shù)據(jù)儲存、編程接口等特點。

2.2 程序計算流程

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

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

3 計算結(jié)果與驗證

3.1 H2熱物性驗證

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

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

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

(47)

式中，X代表相應(yīng)的熱物性參數(shù)。

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

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

3.2 H-H2熱物性驗證

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

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

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

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

4 總結(jié)

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

本程序可為氫工質(zhì)相關(guān)的航天推進、物理學、能源動力等行業(yè)的科研和應(yīng)用提供支持借鑒。