高压氢气压缩因子计算模型研究

氢气燃烧发热量高且绿色环保；氢能是当今世界最热门的新型能源之一。随着对氢能的开发和利用，对高压氢气的精确计量要求也越来越高。

氢气的计量仪表包括科里奥利流量计、气体涡轮流量计、气体腰轮流量计、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数据与各模型计算结果进行了比较。

确定气体热工状态的主要参数有比容v、压力p和温度T，它们之间有确定的关系，称为气体状态方程，即f(p, v, T)=0。对于理想气体(或称“完全气体”)，气体状态方程遵循Clapeyron方程：

式中：p为气体压力，Pa；v为比容，m3/kg；T为温度，K；R为气体常数，J/(mol·K)。

在分子运动论中，所谓理想气体是指：气体分子几何尺寸无限小，分子与分子之间没有吸引力与排斥力。在压力足够低，密度足够小的情况下，各种气体都接近理想气体性质，式(1)表示的Clapeyron理想气体状态方程具有一定的精确性。

在气体压力较高、密度较大(温度较低)时，气体中分子运动规律复杂。本研究所讨论的高压氢气实际条件(2~100 MPa、50~298 K)即属于这种工况。此时，理想气体状态方程不能用来描述氢气实际状态下的p-v-T关系，而需要用较复杂的实际气体状态方程。一般而言，实际气体状态方程是在理想气体状态方程的基础上，通过理论和实验的方法，对理想气体状态方程进行修正取得的。引入压缩因子Z=f(p, T)，对理想气体状态方程用压缩因子进行修正是实际气体状态方程常用的表达方式。

气体压缩因子Z表示实际气体偏离理想气体特性的程度，是对氢气、天然气等实际气体计量中的一个重要物性参数。压缩因子Z的定义为[3]：

从式(2)知，压缩因子的实质是实际气体或气体在实际状态下对理想气体状态方程(Z=1)的修正系数。

美国燃气协会(AGA)提出了AGA8-92DC方程，采用气体摩尔组成计算方法，给出了工作状态下的天然气压缩因子方程[4-5]，但AGA8-92DC不支持计算氢气的压缩因子。

本研究归纳了几种高压氢气压缩因子的计算模型，并进行了分析讨论。

等容线性方程是指在等容条件下，氢气的p和T满足线性关系：

式中：A为常数项；K为系数。

式(3)中部分A、K值列于表 1中。式(3)计算简单，但误差较大[6]。

对于密度不大的气体，可以采用Virial方程，即：

式中：RM为通用气体常数，J/(mol·K)；VM为摩尔体积，L/mol；B、C为第二、第三Virial系数，可以由式(5)、式(6)确定[6]：

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

式中：c0=1.810 5×10-8 m6/mol2(文献[6]中该数据为1.310 5 ×10-8 m6/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%。

考虑分子自身占有的容积和分子之间的相互作用力，Van der Waals对理想气体状态方程作了修正，提出了V-d-W方程，即：

式中：下标c表示临界状态，对于氢气，其临界参数pc=1.296 6 MPa，Tc=33.24 K，vc=0.032 2 m3/kg，Zc=0.304。

将式(7)写成对比态形式：

式中：$p_{\mathrm{r}}=\frac{p}{p_{\mathrm{c}}}, T_{\mathrm{r}}=\frac{T}{T_{\mathrm{c}}}$，二者称为对比态参数。

式(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方程计算高压氢气压缩因子。

R-K方程改进了V-d-W方程，因此，R-K方程可以获得更精确的压缩因子。R-K方程如式(11)所示：

式中：a0、b0、c0、b为与气体临界参数有关的常数；n为与物质有关的常数。

将氢气的临界参数代入上式，计算得到a0=609 9.48 J/mol，b0=9×10-3 m3/kg，c0=3.5×10-6 m3/kg。可见c0＜＜b0，因此，取b≈b0。

将式(10)写成对比态形式：

对于天然气，取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所示。

将NIST给出的高压氢气的密度ρ根据式(2)换算成其压缩因子Z。采用“平滑散点图”方式绘出Z-p-T曲线图，如图 7所示。

若以p/T作为横坐标，则当p/T≥1 MPa/K时，Z-p/T可以拟合成一直线方程，如图 8所示，拟合的直线方程为：

式中：p为氢气压力，MPa；T为氢气温度，K。

式(17)计算结果显示：当p/T=1~2 MPa/K时，线性误差 < 0.55%；当p/T < 1 MPa/K时，各温度下的压缩因子曲线严重发散(见图 9)，因此，该区间内不能用式(17)或其他形式的线性方程来计算压缩因子。

(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时，由于线性误差较大，因此不能拟合成任何形式的线性方程。