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

(1.西安交通大學(xué) 能源與動(dòng)力工程學(xué)院 動(dòng)力工程多相流國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049； 2.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院 陜西省先進(jìn)核能技術(shù)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049)

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

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

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

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

1 熱物性計(jì)算模型

1.1 分子氫熱物性

1.1.1狀態(tài)方程 分子氫氣作為實(shí)際氣體考慮，狀態(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；下標(biāo)cr表示臨界狀態(tài)。

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

h=h0(T)+pV-RgT-

(7)

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

(8)

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

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

(9)

(10)

式中：ui、vi為系數(shù)；下標(biāo)p表示定壓。

比定容熱容cV為：

(11)

式中，下標(biāo)T、V分別表示定溫、定容。

1.1.4聲速 聲速Csound計(jì)算如下：

(12)

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

1.1.5黏度 黏度μ計(jì)算采用Muzny模型[4]：

(13)

(14)

(15)

1.1.6導(dǎo)熱系數(shù) 導(dǎo)熱系數(shù)λ計(jì)算采用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擴(kuò)散系數(shù) 擴(kuò)散系數(shù)D計(jì)算模型[20]如下：

(20)

1.2 原子氫熱物性

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

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

pV=RgT

(21)

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

1.2.2比焓 比焓計(jì)算模型[21]如下：

(22)

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

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

(23)

1.2.4聲速 聲速計(jì)算如下：

(24)

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

1.2.5黏度 黏度計(jì)算模型[20]如下：

(25)

1.2.6導(dǎo)熱系數(shù) 導(dǎo)熱系數(shù)計(jì)算模型[20]如下：

(26)

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

(27)

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

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

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

(28)

(29)

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

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

(30)

值得注意的是，T≤1 000 K時(shí)Kp?10-10，因此本研究認(rèn)為T(mén)≤1 000 K時(shí)的Kp=0。

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

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

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

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

(31)

式中：wH為原子氫質(zhì)量分?jǐn)?shù)；下標(biāo)H-H2表示混合氫工質(zhì)。

1.3.3比焓、比熱容、聲速 比焓、比熱容、聲速加權(quán)計(jì)算方式[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為標(biāo)準(zhǔn)摩爾生成焓，即1 mol H2在標(biāo)準(zhǔn)狀態(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黏度 黏度計(jì)算采用文獻(xiàn)[7-8]中的模型：

(37)

(38)

(39)

(40)

1.3.5導(dǎo)熱系數(shù) 導(dǎo)熱系數(shù)計(jì)算采用文獻(xiàn)[7-8]中的模型：

λH-H2=λf+λr

(41)

(42)

(43)

(44)

(45)

(46)

2 程序開(kāi)發(fā)

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

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

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

2.1 程序組成與功能

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

2.2 程序計(jì)算流程

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

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

3 計(jì)算結(jié)果與驗(yàn)證

3.1 H2熱物性驗(yàn)證

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

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

圖4 H2熱物性對(duì)比Fig.4 Comparison of H2 thermophysical property

(47)

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

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

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

3.2 H-H2熱物性驗(yàn)證

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

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

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

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

4 總結(jié)

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

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