石油与天然气化工  2023, Vol. 52 Issue (1): 103-109
引用本文
屈宏强, 邢静芳, 郑雪晶, 胡方舒. 基于燃烧法的天然气发热量直接测量不确定度评估[J]. 石油与天然气化工, 2023, 52(1): 103-109.   DOI: 10.3969/j.issn.1007-3426.2023.01.017
基于燃烧法的天然气发热量直接测量不确定度评估
屈宏强1,2 , 邢静芳1,2 , 郑雪晶3 , 胡方舒3     
1. 河北省计量监督检测研究院;
2. 河北省气体计量与大数据分析重点实验室;
3. 天津大学环境科学与工程学院
收稿日期：2022-03-01；修回日期：2022-06-22
基金项目：国家市场监督管理总局科技计划项目“天然气发热量测量方法及其溯源研究”(2020MK011)
作者简介：屈宏强，1965年生，正高级工程师，主要从事流量计量、天然气计量研究工作。E-mail：qhqganen@sina.com.
通信作者：邢静芳，1985年生，高级工程师，主要从事流量计量、天然气计量研究工作。E-mail：xjfhbjl@126.com.
摘要：目的 准确测量天然气发热量，实现天然气贸易交接公平准确。方法 以定质量式开放火焰燃烧法为原理设计了一种扩展不确定度优于0.3%(k=2)的天然气发热量直接测量实验台；在发热量理论公式的基础上，建立了天然气发热量测量数学模型，综合考虑经济成本，选定了具有适当技术指标的测量设备，以甲烷为气体样品，对天然气发热量直接测量实验台进行了不确定度评估，确定了该实验台研制的可行性，并分析了各不确定度分量的贡献率。结果 该实验台测量甲烷发热量的不确定度为0.27%(k=2)。结论 满足常用的1.0级及其以下的热值仪的发热量的量值溯源，为实验台搭建奠定了基础。
关键词：不确定度    发热量    甲烷    燃烧法    直接测量    
Evaluation of direct measurement uncertainty of natural gas calorific value based on combustion method
Qu Hongqiang1,2 , Xing Jingfang1,2 , Zheng Xuejing3 , Hu Fangshu3     
1. Institute of Metrology of Hebei Province, Shijiazhuang, Hebei, China;
2. Hebei Key Laboratory of Gas Metrology and Big Data Analysis, Shijiazhuang, Hebei, China;
3. School of Environmental Science and Engineering, Tianjin University, Tianjin, China
Abstract: Objective The aims are to measure the calorific value of natural gas accurately and achieve a fair and accurate natural gas trade. Methods Based on the principle of the fixed-mass open flame combustion method, a direct measurement experimental platform of natural gas calorific value is designed, which has an extended uncertainty better than 0.3% (k=2). Based on the theoretical formula of calorific value, a mathematical model for measuring the calorific value of natural gas is established. Considering the economic cost comprehensively, the measurement equipment with appropriate technical indicators is selected. The uncertainty of the direct measurement experiment platform of natural gas calorific value is evaluated by taking methane as the gas sample, the feasibility of the development of the experimental platform is confirmed, and the contribution rate of each uncertainty component is analyzed. Results The uncertainty of measuring methane calorific value is 0.27%(k=2). Conclusions The uncertainty of measuring methane calorific value meets the value traceability of the calorific value of the commonly used 1.0 level and below calorific value meter, which lays a foundation for the establishment of the experimental platform.
Key words: uncertainty    calorific value    methane    combustion method    direct measurement    

天然气作为一种清洁能源，在全球范围内的需求量持续增加，相应地，天然气贸易量亦持续增长，合理的计量方式是维持天然气贸易准确公平的重要保证。目前，天然气贸易中的计量方式主要有3种，分别为体积计量、质量计量和能量计量[1]。天然气的使用价值在于提供热量，相同体积或质量的天然气可能由于气体组分不同，从而发热量不同，若按体积或质量进行计量与计价则会造成天然气的使用价值与价格不符。2019年5月，国家发展改革委、国家能源局、住房城乡建设部和市场监管总局联合印发的《油气管网设施公平开放监管办法》要求国家推行天然气能量计量计价[2]。

为确保政策的实施，须准确测量天然气发热量，而甲烷作为天然气的主要组分，准确测得甲烷发热量是天然气发热量准确测定的基础。研究表明，基于定质量式开放火焰燃烧法的等周热量计在燃气发热量测量方面具有最高的准确度[3]。该类型热量计最早是在20世纪30年代由美国国家标准局的Rossini提出的，研究表明[4]，基于Rossini测量原理的装置的不确定度可以达到0.05%(k=2)，由这套装置测得的纯甲烷气体的摩尔发热量为890.578 kJ/mol, 与ISO 6976-1995《Natural gas-Calculation of calorific values, density, relative density and Wobbe index from composition》中给出的890.63 kJ/mol非常接近[5-8]。Cutler-Hammer热量计是ISO 15971-2008《Natural gas-Measurement of properties-Calorific value and Wobbe index》推荐的1级基准装置，最高可以达到0.17%(k=2)的不确定度水平[9]。

GB/T 18603-2014《天然气计量系统技术要求》中A级规定的发热量测定不确定度≤1%[10]。目前，我国未能设计并搭建出完整的符合准确度等级要求的天然气发热量直接测量实验台。为此，河北省计量监督检测研究院设计了一种高准确度燃气发热量直接测量实验台，为确定该实验台的准确度，本研究以甲烷为气体样品，对该实验台进行了详细的不确定度分析与计算。

1 天然气发热量直接测量实验原理

图 1所示为天然气发热量直接测量实验台的原理图，其原理为将一定量的甲烷、氧气与氩气通入到燃烧室内进行缓慢定压燃烧，燃烧过程中产生的热量被量热容器内的吸热介质(水)所吸收，通过测量甲烷燃烧过程中吸热介质的温升、热量计的当量热容及参与燃烧的甲烷质量等相关参数来计算得到甲烷的发热量。

图 1     天然气发热量直接测量实验台原理图 1-甲烷气瓶; 2-氩气气瓶; 3-氧气气瓶; 4~6-质量流量控制器; 7-电子天平; 8-储气小球; 9-参考球; 10-热敏电阻; 11-搅拌器; 12-导流筒; 13-燃烧室; 14-量热容器内壳体; 15-量热容器外壳体; 16-恒温槽; 17-干燥管; 18-气泵; 19-红外线气体分析仪

天然气发热量直接测量实验台的供气系统如图 1所示，由甲烷气瓶、氩气气瓶、氧气气瓶、质量流量控制器及相应阀门、管道、管件等组成。各供气管道上的质量流量控制器可用于控制天然气、氩气、氧气以一定的流量通入到燃烧室中进行燃烧反应。氩气虽然不参与反应，但可以起到稳定火焰的作用[11]。天然气质量测量系统由电子天平、储气小球与参考球构成。天然气发热量测量主体由热量计、计算机、数字万用表组成。其中，热量计包括燃烧室、导流筒、搅拌器、热敏电阻、量热容器、恒温槽。

2 数学模型

甲烷发热量的计算公式如式(1)所示。

$ H_{{\mathrm{CH}}_4}=\frac{C_{{\text {cal }}} \Delta T_{{\mathrm{ad}}, {\mathrm{comb}}}+K}{m_{{\mathrm{CH}}_4}} $ (1)

式中：HCH4为甲烷发热量，J/g；Ccal为热量计热容，J/℃；ΔTad, comb为甲烷燃烧过程中吸热介质的绝热温升，℃；mCH4为参与燃烧的甲烷质量，g；K为能量修正，J。

2.1 吸热介质绝热温升

甲烷燃烧过程中吸热介质的温度变化见图 2。吸热介质的绝热温升计算公式如式(2)~式(4)所示。

$ \Delta T_{{\mathrm{ad}}}=T_{{\mathrm{e}}}-T_{{\mathrm{b}}}-T_{\mathtt{δ}} $ (2)
$ T_{\mathtt{δ}}=k_{{\mathrm{c}}}\left(T_{\infty}-T_{{\mathrm{m}}}\right)\left(t_{{\mathrm{e}}}-t_{{\mathrm{b}}}\right) $ (3)
$ T_{{\mathrm{m}}}=\frac{1}{t_{{\mathrm{e}}}-t_{{\mathrm{b}}}} \int_{t_{{\mathrm{b}}}}^{t_{{\mathrm{e}}}} T_2(t) {\mathrm{d}} t $ (4)
图 2     吸热介质温度变化曲线

式中：ΔTad为吸热介质绝热温升，℃；Te为主升阶段结束时吸热介质的温度，℃；Tb为主升阶段开始时吸热介质的温度，℃；Tδ为温度修正，℃；kc为冷却系数，s-1；T∞为时间足够长后吸热介质达到稳态时的温度，℃；Tm为主升阶段吸热介质的平均温度，℃；te为主升阶段结束时刻，s；tb为主升阶段开始时刻，s；t为燃烧时间，s。

2.2 热量计热容

为保证发热量测量准确，在燃烧实验前后分别进行一次热容测定试验，取两次的算术平均值作为热量计热容的最终测量结果。热容的计算公式如式(5)~式(6)所示。

$ C_{{\mathrm{cal}}}=\frac{E_{{\mathrm{cal}}, 1}}{2 \Delta T_{{\mathrm{ad}}, {\mathrm{cal}}, 1}}+\frac{E_{{\mathrm{cal}}, 2}}{2 \Delta T_{{\mathrm{ad}}, {\mathrm{cal}}, 2}} $ (5)
$ E_{\mathrm {cal }}=\frac{\sum\limits_{i=1}^n U_{\mathrm {heat }, i} U_{\mathrm {ref }, i} {\mathrm{~d}} t_i}{R_{\mathrm {ref }}} $ (6)

式中：Ecal, 1为燃烧前热量计热容测定实验时的电加热量，J；Ecal, 2为燃烧后热量计热容测定实验时的电加热量，J；ΔTad, cal, 1为燃烧前热量计热容测定实验时吸热介质的绝热温升，℃；ΔTad, cal, 2为燃烧后热量计热容测定实验时吸热介质的绝热温升，℃；Uheat, i为第i个采样点的电加热带两端的电压，V；Uref, i为第i个采样点的标准电阻两端的电压，V；dti为第i个采样点时间，s；Rref为标准电阻的阻值，Ω；n为采样数量, 个。

2.3 参与燃烧的甲烷质量

参与燃烧的甲烷质量计算公式如式(7)~式(10)所示。

$ m_{{\mathrm{CH}}_4}=\Delta m_{{\text {before }}}-\Delta m_{{\text {after }}}-m_{\mathtt{δ}} $ (7)
$ \Delta m_{{\text {before }}}=\frac{m_{{\mathrm{M}} 1 {\text {, before }}}+m_{{\mathrm{M}} 2, {\text { before }}}}{2}-\frac{m_{{\mathrm{C}} 1 {\text {, before }}}+m_{{\mathrm{C}} 2, {\text { before }}}}{2} $ (8)
$ \Delta m_{{\text {after }}}=\frac{m_{{\mathrm{M}} 1, {\text { after }}}+m_{{\mathrm{M}} 2 {\text {, after }}}}{2}-\frac{m_{{\mathrm{C}} 1, {\text { after }}}+m_{{\mathrm{C}} 2 {\text { after }}}}{2} $ (9)
$ m_{\mathtt{δ}}=\varphi_{{\mathrm{CH}}_4} \rho_{{\mathrm{CH}}_4} V_{{\text {total }}} \times 10^{-6} $ (10)

式中：Δmbefore为燃烧前充好甲烷的储气小球与参考球的质量之差，g；Δmafter为燃烧后储气小球与参考球的质量之差，g；mM1, before为燃烧前第1次测得的充好甲烷的储气小球质量，g；mM1, after为燃烧后第1次测得的储气小球质量，g；mM2, before为燃烧前第2次测得的充好甲烷的储气小球质量，g；mM2, after为燃烧后第2次测得的储气小球质量，g；mC1, before为燃烧前第1次测得的参考球质量，g；mC1, after为燃烧后第1次测得的参考球质量，g；mC2, before为燃烧前第2次测得的参考球质量，g；mC2, after为燃烧后第2次测得的参考球质量，g；mδ为未参与燃烧的甲烷质量，g；φCH4为储气钢瓶中甲烷的体积分数，10-6；ρCH4为标况下的甲烷密度，g/L；Vtotal为标况下储气钢瓶内混合气的体积，L。

2.4 能量修正

能量修正计算公式如式(11)所示。

$ K=E_{{\mathrm{w}}}+E_{{\mathrm{g}}}-E_{{\mathrm{l}}}-E_{{\mathrm{f}}} $ (11)

式中：Ew为燃烧产物水蒸气未完全冷凝引入的能量修正，J；Eg为进出燃烧室气体温度不同引入的能量修正，J；El为点火瞬间点火能量引入的能量修正，J；Ef为搅拌器搅拌吸热介质过程中引入的能量修正，J。

3 不确定度分析与计算
3.1 甲烷发热量A类不确定度

本节参考国外测得的甲烷发热量实验数据进行甲烷发热量的A类不确定度计算。甲烷发热量的实验数据如表 1所列[12]，其平均值为55 507.996 J/g。

表 1    甲烷发热量实验数据 

甲烷发热量A类不确定度计算公式如式(12)所示：

$ u_{{\mathrm{A}}}\left(H_{{\mathrm{CH}}_4}\right)=\sqrt{\frac{1}{n_1-1} \sum\limits_{i=1}^{n_1}\left(H_{{\mathrm{CH}}_4, i}-\overline{H_{{\mathrm{CH}}_4}}\right)^2} $ (12)

式中：uA(HCH4)为甲烷发热量A类不确定度, J/g; HCH4, i为第i次实验的甲烷发热量，J/g；$\overline{H_{{\mathrm{CH}}_4}}$为甲烷发热量的平均值，J/g；n1为实验次数。

根据表 1中的数据与式(12)，计算得到甲烷发热量A类不确定度为16.441 J/g，考虑到实验台为项目组首次搭建，存在不确定度因素。因此，在上述结果的基础上，将A类不确定度进行了扩大处理，取值为原来的2倍，即32.882 J/g。

3.2 甲烷发热量B类不确定度

根据式(1)、式(5)与不确定度传播律[12]，可得甲烷发热量B类不确定度，其计算公式如式(13)~式(15)所示。

$ \begin{array}{r} u_{{\mathrm{B}}}^2\left(H_{{\mathrm{CH}}_4}\right)=\left(\frac{E_{{\mathrm{cal}}, 1}}{2 m_{{\mathrm{CH}}_4}}\right)^2 u_{{\mathrm{B}}}^2\left(Z_1\right)+\left(\frac{E_{{\mathrm{cal}}, 2}}{2 m_{{\mathrm{CH}}_4}}\right)^2 u_{{\mathrm{B}}}^2\left(Z_2\right)+ \\ \left(\frac{Z_1}{2 m_{{\mathrm{CH}}_4}}\right)^2 u_{{\mathrm{B}}}^2\left(E_{{\mathrm{cal}}, 1}\right)+\left(\frac{Z_2}{2 m_{{\mathrm{CH}}_4}}\right)^2 u_{{\mathrm{B}}}^2\left(E_{{\mathrm{cal}}, 2}\right)+ \\ \left(\frac{H_{{\mathrm{CH}}_4}}{m_{{\mathrm{CH}}_4}}\right)^2 u_{{\mathrm{B}}}^2\left(m_{{\mathrm{CH}}_4}\right)+\left(\frac{1}{m_{{\mathrm{CH}}_4}}\right)^2 u_{{\mathrm{B}}}^2(K) \end{array} $ (13)
$ Z_1=\frac{\Delta T_{{\mathrm{ad}}, {\mathrm{comb}}}}{\Delta T_{{\mathrm{ad}}, {\mathrm{cal}}, 1}} $ (14)
$ Z_2=\frac{\Delta T_{{\mathrm{ad}}, {\mathrm{comb}}}}{\Delta T_{{\mathrm{ad}}, {\mathrm{cal}}, 2}} $ (15)

式中：uB(HCH4)为甲烷发热量B类不确定度，J/g；Z1为甲烷燃烧实验中吸热介质绝热温升与燃烧前热量计热容测定时吸热介质绝热温升之间的比例；Z2为甲烷燃烧实验中吸热介质绝热温升与燃烧后热量计热容测定时吸热介质绝热温升之间的比例。

由于燃烧实验前后的热量计热容测定实验中采用的实验设备与流程均一致，因此, Z1≈Z2、uB(Z1)≈uB(Z2)、Ecal, 1≈Ecal, 2、u(Ecal, 1)≈u(Ecal, 2)。

令：

$ u\left(Z_1\right) \approx u\left(Z_2\right)=u(Z) $ (16)
$ u\left(E_{{\mathrm{cal}}, 1}\right) \approx u\left(E_{{\mathrm{cal}}, 2}\right)=u\left(E_{{\mathrm{cal}}}\right) $ (17)

则甲烷发热量B类不确定度计算公式见式(18)：

$ \begin{aligned} & u_{{\mathrm{B}}}\left(H_{{\mathrm{CH}}_4}\right)= \\ & \sqrt{\frac{\frac{E_{{\mathrm{cal}}}{ }^2 u_{{\mathrm{B}}}{ }^2(Z)}{2}+\frac{Z^2 u_{{\mathrm{B}}}^2\left(E_{{\mathrm{cal}}}\right)}{2}+H_{{\mathrm{CH}}_4}{ }^2 u_{{\mathrm{B}}}{ }^2\left(m_{{\mathrm{CH}}_{{\mathrm{4}}}}\right)+{u_{{\mathrm{B}}}}^2(K)}{m_{{\mathrm{CH}}_4}^2}} \\ & \end{aligned} $ (18)

甲烷发热量B类不确定度的影响因素分析图见图 3。

图 3     B类不确定度的影响因素分析图

下面对式(18)中的各B类不确定度分项进行分析与计算。

3.2.1 绝热温升比例引入的B类不确定度

根据不确定度传播律[12]，由绝热温升比例引入的B类不确定度uB2(Z)计算公式如式(19)所示：

$ \begin{aligned} u_{{\mathrm{B}}}^2(Z)= & \frac{1}{\Delta T_{{\mathrm{ad}}, {\mathrm{cal}}}{ }^2} u_{{\mathrm{B}}}^2\left(\Delta T_{{\mathrm{ad}}, {\text { comb }}}\right)+ \\ & \frac{\Delta T_{{\mathrm{ad}}, {\mathrm{comb}}}{ }^2}{\Delta T_{{\mathrm{ad}}, {\mathrm{cal}}}{ }^4} u_{{\mathrm{B}}}^2\left(\Delta T_{{\mathrm{ad}}, {\mathrm{cal}}}\right)- \\ & 2 \frac{\Delta T_{{\mathrm{ad}}, {\text { comb }}}}{\Delta T_{{\mathrm{ad}}, {\text { cal }}}{ }^3} u_{{\mathrm{B}}}\left(\Delta T_{{\mathrm{ad}}, {\text { comb }}}, \Delta T_{{\mathrm{ad}}, {\text { cal }}}\right) \end{aligned} $ (19)

式中：ΔTad, comb为燃烧过程中的绝热温升，℃；Tb, comb为燃烧过程吸热介质的起始温度，℃；Te, comb为燃烧过程吸热介质的最终温度，℃；ΔTad, cal为电加热过程中的绝热温升，℃；Tb, cal为电加热过程吸热介质的起始温度，℃；Te, cal为电加热过程吸热介质的最终温度，℃。

根据F Haloua等的实验数据，Z=1.016[12]，据此假设ΔTad, comb=3.00 ℃、Tb, comb=Tb, cal=23.50 ℃，则ΔTad, cal=2.95 ℃、Te, comb=26.50 ℃、Te, cal=26.45 ℃，则最终Z=1.017。

3.2.1.1 uB(ΔTad, comb)与uB(ΔTad, cal)分析与计算

根据式(2)与不确定度传播律[12]，由绝热温升引起的B类不确定度u2(ΔTad)如式(20)所示：

$ \begin{aligned} u_{{\mathrm{B}}}{ }^2\left(\Delta T_{{\mathrm{ad}}}\right)= & u_{{\mathrm{B}}}{ }^2\left(T_{{\mathrm{e}}}\right)+u_{{\mathrm{B}}}{ }^2\left(T_{{\mathrm{b}}}\right)+u_{{\mathrm{B}}}{ }^2\left(T_{\mathtt{δ}}\right)- \\ & 2 u_{{\mathrm{B}}}\left(T_{{\mathrm{e}}}, T_{{\mathrm{b}}}\right)-2 u_{{\mathrm{B}}}\left(T_{{\mathrm{e}}}, T_{\mathtt{δ}}\right)+ \\ & 2 u_{{\mathrm{B}}}\left(T_{{\mathrm{b}}}, T_{\mathtt{δ}}\right) \end{aligned} $ (20)

式中: Te为加热过程吸热介质的最终温度，℃；Tb为加热过程吸热介质的起始温度，℃；Tδ为温度修正，℃。

式(20)中，由温度引起的B类不确定度见式(21)。

$ u_{{\mathrm{B}}}^2(T)=u_{{\text {ther }}}^2(T)+u_{{\text {multi }}}^2(T)+u_{{\text {hom }}}^2(T) $ (21)

式中：uther(T)为热敏电阻引入的不确定度，℃；umulti(T)为数字多用表引入的不确定度，℃；uhom(T)为吸热介质温度均匀性引入的不确定度，℃。

对于式(20)中的其他项，由于还未进行实验，缺少数据，参考相关文献中的数据[12]，即：

$ \begin{aligned} & u_{{\mathrm{B}}}^2\left(T_{\mathtt{δ}}\right)-2 u_{{\mathrm{B}}}\left(T_{{\mathrm{e}}}, T_{{\mathrm{b}}}\right)-2 u_{{\mathrm{B}}}\left(T_{{\mathrm{e}}}, T_{\mathtt{δ}}\right)+2 u_{{\mathrm{B}}}\left(T_{{\mathrm{b}}}, T_{\mathtt{δ}}\right) \\ & =-4.921 \times 10^{-5} \end{aligned} $

综上，根据绝热温升不确定度计算公式，可得uB(ΔTad, comb)与uB(ΔTad, cal) 均为1.887×10-3 ℃。

3.2.1.2 uB(ΔTad, comb, ΔTad, cal)分析与计算

为了最大化测量不确定度，将uB(ΔTad, comb, ΔTad, cal)设为0。

综上，根据上述计算结果与式(19)计算得绝热温升比例引入的B类不确定度为9.062×10-4 ℃。

根据式(19)，计算得uB(Z)为1.337×10-3。

3.2.2 电加热量引入的B类不确定度

为简化符号，令Uheat, iUref, idt=Ai，则式(6)可转化为下式：

$ E_{{\mathrm{cal}}}=\frac{\sum\limits_{i=1}^n A_i}{R_{{\mathrm{ref}}}} $ (22)

根据式(22)与不确定度传播律[12]，由电加热量引入的B类不确定度如式(23)所示：

$ \begin{aligned} u_B^2\left(E_{{\mathrm{cal}}}\right)= & \frac{n}{R_{{\mathrm{ref}}}{ }^2} u_B^2\left(A_i\right)+\frac{\left(\sum\limits_{i=1}^n A_i\right)^2}{R_{{\mathrm{ref}}}{ }^4} u_B^2\left(R_{{\mathrm{ref}}}\right)+ \\ & \frac{n(n-1)}{R_{{\mathrm{ref}}}{ }^2} u_B\left(A_i, A_j\right) \end{aligned} $ (23)

式中：j=i+1。

3.2.2.1 uB(Ai)分析与计算

根据不确定度传播律，由于多次测量引入的B类不确定度uB2(Ai)如式(24)所示：

$ \begin{aligned} u_{{\mathrm{B}}}^2\left(A_i\right)= & \left(\overline{U_{{\text {ref }}}} \cdot \overline{{\mathrm{d}} t}\right)^2 u_{{\mathrm{B}}}^2\left(U_{{\text {heat }}, i}\right)+ \\ & \left(\overline{U_{{\text {heat }}}} \cdot \overline{{\mathrm{d}} t}\right)^2 u_{{\mathrm{B}}}^2\left(U_{{\text {ref }}, i}\right)+ \\ & \left(\overline{U_{{\text {heat }}}} \cdot \overline{U_{{\text {ref }}}}\right)^2 u_{{\mathrm{B}}}^2\left({\mathrm{~d}} t_i\right) \end{aligned} $ (24)

考虑各个设备仪器的精度，按照均匀分布计算得到uB(Uheat, i)、uB(Uref, i)和uB(dti):

$ \begin{array}{l} u_{{\mathrm{B}}}\left(U_{{\text {heat }}, i}\right) =\frac{41.355 \times 0.004\;5 \%+100 \times 0.000\;6 \%}{\sqrt{3}} \\ =2.734 \times 10^{-4} {\mathrm{~V}}\\ u_{{\mathrm{B}}}\left(U_{{\text {ref }}, i}\right) =\frac{0.115 \times 0.004 \%+1 \times 0.000\;7 \%}{\sqrt{3}} \\ =1.289 \times 10^{-6} {\mathrm{~V}} \\u_{{\mathrm{B}}}\left({\mathrm{d}} t_i\right)=\frac{0.01 \%}{\sqrt{3}}=5.774 \times 10^{-5} {\mathrm{~s}} \end{array} $

根据相关文献中的实验数据[12]，可知：

$ \overline{U_{{\text {heat }}}}=41.355 {\mathrm{~V}} ; \overline{U_{{\text {ref }}}}=0.115 {\mathrm{~V}} ; \overline{{\mathrm{d}} t}=1 {\mathrm{~s}} $

根据式(24)，可得uB(Ai)为3.176× 10-4 J·Ω。

3.2.2.2 uB(Rref)分析与计算

标准电阻选用BZ3-0.1Ω，其准确度等级为0.01%，则标准电阻引入的不确定度为：

$ u_{{\mathrm{B}}}\left(R_{{\mathrm{ref}}}\right)=\frac{0.1 \times 0.000\;1}{\sqrt{3}}=5.774 \times 10^{-6} \;\Omega $
3.2.2.3 uB(Ai, Aj)分析与计算

uB(Ai, Aj)计算公式为：

$ u_B\left(A_i, A_j\right)={\overline{U_{{\text {ref }}}}}^2 \overline{{\mathrm{d}} t}^2 u_B\left(U_{{\text {heat }}, i}, U_{{\text {heat }}, j}\right)+\\{\overline{U_{{\text {heat }}}}}^2{\overline{{\mathrm{d}} t}}^2 u_B\left(U_{{\text {ref }}, i}, U_{{\text {ref }}, j}\right)+\\{\overline{U_{{\text {heat }}}}}^2{\overline{U_{{\text {heat }}}}}^2{\overline{U_{{\text {ref }}}}}^2 u_B\left({\mathrm{~d}} t_i, {\mathrm{~d}} t_j\right) $ (25)

uB(Ai, Aj)分析与计算：

$ u_{{\mathrm{B}}}\left(U_{{\text {heat }}, i}, U_{{\text {heat }}, j}\right)=u_{{\mathrm{B}}}^2\left(U_{{\text {heat }}}\right) $ (26)
$ u_{{\mathrm{B}}}\left(U_{{\text {ref }}, i}, U_{{\text {ref }}, j}\right)=u_{{\mathrm{B}}}^2\left(U_{{\text {ref }}}\right) $ (27)
$ u_{{\mathrm{B}}}\left({\mathrm{d}} t_i, {\mathrm{~d}} t_j\right)=u_{{\mathrm{B}}}^2({\mathrm{~d}} t) $ (28)

根据式(26)~式(28)，计算得uB(Uheat, i, Uheat, j)为7.477×10-8 V2，uB(Uref, i, Uref, j)为1.661×10-12 V2，uB(dti, dtj)为2.500×10-9 V2。

综上，根据uB(Ai, Aj)计算公式，可得uB(Ai, Aj)为6.037×10-8 J2·Ω2。

根据实验数据[13]：

$ \sum\limits_{i=1}^n A_i=5\;696.486 {\mathrm{~J}} \cdot \Omega ; n=1\;200 $

根据式(23)，计算得uB(Ecal)为4.417 J。

3.2.3 参与燃烧的甲烷质量引入的B类不确定度

根据式(7)与不确定度传播律[12]，由参与燃烧的甲烷质量引入B类不确定度uB2(mCH4)的计算公式, 如式(29)所示：

$ u_{{\mathrm{B}}}^2\left(m_{{\mathrm{CH}}_4}\right)=u_{{\mathrm{B}}}^2\left(\Delta m_{{\text {before }}}\right)+u_{{\mathrm{B}}}^2\left(\Delta m_{{\text {after }}}\right)+u_{{\mathrm{B}}}^2\left(m_{\mathtt{δ}}\right) $ (29)

式中：Δmbefore为燃烧前充好甲烷的储气小球与参考球之间的质量之差，g；Δmafter为燃烧后充好甲烷的储气小球与参考球之间的质量之差，g；mδ为未参与燃烧的甲烷质量，g。

3.2.3.1 uB(Δm)分析与计算

为进行符号简化，将式(8)、式(9)简化为下式：

$ \Delta m=\frac{m_{{\mathrm{M}} 1}+m_{{\mathrm{M}} 2}}{2}-\frac{m_{{\mathrm{C}} 1}+m_{{\mathrm{C}} 2}}{2} $ (30)

式中：mM1为充好甲烷的储气小球质量的第1次测量值，g；mM2为充好甲烷的储气小球质量的第2次测量值，g；mC1为参考球质量的第1次测量值，g；mC2为参考球质量的第2次测量值，g。

根据式(30)与不确定度传播律[12]，由参考球与储气小球质量差引入B类不确定度uB2(Δm)计算公式, 如式(31)所示：

$ u_{{\mathrm{B}}}^2(\Delta m)=\frac{u_{{\mathrm{B}}}^2\left(m_{{\mathrm{M}} 1}\right)+u_{{\mathrm{B}}}^2\left(m_{{\mathrm{M}} 2}\right)+u_{{\mathrm{B}}}^2\left(m_{{\mathrm{C}} 1}\right)+u_{{\mathrm{B}}}^2\left(m_{{\mathrm{C}} 2}\right)}{4} $ (31)

采用MS204TS分析天平测量甲烷储气小球和参考球的质量，其可读性为0.1 mg，最大允许误差为1 mg，按均匀分布计算，可得：

$ \begin{gathered} u_{{\mathrm{B}}}\left(m_{{\mathrm{M}} 1}\right)=u_{{\mathrm{B}}}\left(m_{{\mathrm{M}} 2}\right)=u_{{\mathrm{B}}}\left(m_{{\mathrm{C}} 1}\right)=u_{{\mathrm{B}}}\left(m_{{\mathrm{C}} 2}\right)= \\ \frac{1 \times 10^{-3}}{\sqrt{3}}=5.774 \times 10^{-4} {\mathrm{~g}} \end{gathered} $

因此，根据式(31)，uB(Δmbefore)与uB(Δmafter)为5.774×10-4 g。

3.2.3.2 uB(mδ)分析与计算

根据式(10)与不确定度传播律[12]，由未参与燃烧的甲烷质量引入的B类不确定度uB2(mδ)计算公式, 如式(32)所示：

$ u_{{\mathrm{B}}}^2\left(m_{\mathtt{δ}}\right)=\left(\varphi_{{\mathrm{CH}}_4} \rho_{{\mathrm{CH}}_4} \times 10^{-6}\right)^2 u_{{\mathrm{B}}}^2\left(V_{{\text {total }}}\right) $ (32)

由于缺乏甲烷含量等参数，因此采用相关文献中的数据[12]，并将文献中流量计的准确度等级为0.5%改为本实验台所采用流量计的准确度1.0%，经计算，得uB(mδ)为1.468×10-4 g。

根据式(29)，计算得到uB(mCH4)为8.296×10-4 g。

3.2.4 能量修正引入的B类不确定度

根据不确定度传播律[13]，由能量修正引入的B类不确定度如式(33)所示：

$ u_{{\mathrm{B}}}^2(K)=u_{{\mathrm{B}}}^2\left(E_{{\mathrm{w}}}\right)+u_{{\mathrm{B}}}^2\left(E_{{\mathrm{g}}}\right)+u_{{\mathrm{B}}}^2\left(E_{{\mathrm{l}}}\right)+u_{{\mathrm{B}}}^2\left(E_{{\mathrm{f}}}\right) $ (33)

式中：Ew为水蒸气未完全冷凝引入的能量修正，J；Eg为进出燃烧室气体温度不同引入的能量修正，J；El为点火瞬间点火能量引入的能量修正，J；Ef为搅拌器引入的能量修正，J。

由于许多参数未知，因此, 参考相关文献中的结论，即能量修正项引入的不确定度占总合成不确定度的0.028%，因此, 推导得出计算公式如式(34)所示：

$ \sqrt{\frac{14}{99\;972}\left[E_{{\mathrm{cal}}}{ }^2 u_{{\mathrm{B}}}{ }^2({\mathrm{Z}})+Z^2{u_{{\mathrm{B}}}}^2\left(E_{{\mathrm{cal}}}\right)+2 H_{{\mathrm{CH}}_{{\mathrm{4}}}}{ }^2 u_{{\mathrm{B}}}{ }^2\left(m_{{\mathrm{CH}}_{{\mathrm{4}}}}\right)\right]} $ (34)

根据式(34)，计算得uB(K)为1.415 J。

参考相关数据[12]，可得：Ecal=56 956.888 J; mCH4=1.057 g; HCH4=55 507.996 J/g。

根据式(18)与上述各不确定度分项计算结果，得出uB(HCH4)为67.102 J/g。

3.3 甲烷发热量合成不确定度

发热量HCH4各分量不确定度及贡献率如表 2所列。

表 2    发热量HCH4各分量不确定度及贡献率

甲烷发热量合成标准不确定度计算公式为：

$ u\left(H_{{\mathrm{CH}}_4}\right)=\sqrt{u_{{\mathrm{A}}}^2\left(H_{{\mathrm{CH}}_4}\right)+u_{{\mathrm{B}}}^2\left(H_{{\mathrm{CH}}_4}\right)} $ (35)

甲烷发热量相对标准不确定度计算公式为：

$ u_{{\mathrm{r}}}\left(H_{{\mathrm{CH}}_4}\right)=\frac{u\left(H_{{\mathrm{CH}}_4}\right)}{\overline{H_{{\mathrm{CH}}_4}}} $ (36)

甲烷发热量扩展相对不确定度计算公式为：

$ U_{{\mathrm{r}}}\left(H_{{\mathrm{CH}}_4}\right)=k u_{{\mathrm{r}}}\left(H_{{\mathrm{CH}}_4}\right) $ (37)

式中：k为包含因子，取2。

根据式(35)~式(37)，计算得甲烷发热量合成标准不确定度u(HCH4)、相对标准不确定度ur(HCH4)与扩展相对不确定度Ur(HCH4)分别为74.726 J/g、0.135%与0.270%。

各项影响因素对于甲烷B类不确定度的贡献率计算公式如式(38)所示[12]：

$ {\mathrm{B}} 类不确定度贡献率 = \frac{{灵敏系数 \times 相应标准不确定度}}{{{\mathrm{B}} 类合成标准不确定度}} $ (38)
4 结论与建议

以甲烷为气体样品，对河北省计量监督检测研究院设计的高准确度天然气发热量直接测量实验台进行了详细的不确定度分析与计算。

甲烷发热量由A类不确定度和B类不确定度构成，经计算与评估，A类不确定度为32.882 J/g，B类不确定度为67.102 J/g，将其合成之后，得到甲烷发热量合成标准不确定度与相对标准不确定度分别为74.726 J/g、0.135%，在包含因子为2的条件下，得到甲烷发热量扩展相对不确定度为0.270%，优于0.3%，满足常用的1.0级及其以下的热值仪的发热量的量值溯源[10]。

通过计算甲烷发热量合成标准不确定度，根据表 2数据可知，由绝热温升比例和参与燃烧的甲烷的质量引入的B类不确定度贡献率最大。由于温度的测量难度较大，涉及到提供恒温环境的恒温槽、测量温度的热敏电阻和测温仪3个设备，环节和影响因素多, 并且质量测量能力相对电参数测量能力较低，总体来说，电加热量的贡献率较低，能量修正本身就是一个很小的量，因此, 其贡献率最小。

根据以上分析，为了进一步优化测量不确定度，在不考虑成本的前提下，可以从温度测量、电加热量测量以及天然气质量测量等方面，采用测量准确度等级更高的计量设备，以提高B类不确定度。同时，还要考虑到环境温度控制和国内机械加工的误差控制，以降低重复性引入的A类不确定度。

参考文献
[1]
黄维和, 罗勤, 黄黎明, 等. 天然气能量计量体系在中国的建设和发展[J]. 石油与天然气化工, 2011, 40(2): 103-108.
[2]
油气管网设施公平开放监管办法[J]. 中国计量, 2019(7): 13-17.
[3]
HALOUA F, HAY B, FILTZ J R. New French reference calorimeter for gas calorific value measurements[J]. Journal of Thermal Analysis and Calorimetry, 2009, 97(2): 673-678. DOI:10.1007/s10973-008-9701-z
[4]
SCHLEY P, BECK M, UHRIG M, et al. Measurements of the calorific value of methane with the new GERG reference calorimeter[J]. International Journal of Thermophysics, 2010, 31(4): 665-679.
[5]
ROSSINI F D. The heats of combustion of methane and carbon monoxide[J]. Bureau of Standards Journal of Research, 1931, 6: 37-49. DOI:10.6028/jres.006.002
[6]
ROSSINI F D. The heat of formation of water and the heats of combustion of methane and carbon monoxide. A correction[J]. Bureau of Standards Journal of Research, 1931, 7: 329-330. DOI:10.6028/jres.007.017
[7]
ARMSTRONG G T, JOBE JR T L. Heating values of natural gas and its components[R]. Washington: National Bureau of Standards, 1982.
[8]
ISO. Natural gas—Calculation of calorific values, density, relative density and Wobbe index from composition: ISO 6976: 1995[S]. ISO, 1995.
[9]
ISO. Natural gas — Measurement of properties—Calorific value and Wobbe index: ISO 15971: 2008[S]. ISO, 2008.
[10]
中华人民共和国国家质量监督检验检疫总局, 中国国家标准化管理委员会. 天然气计量系统技术要求: GB/T 18603-2014[S]. 北京: 中国标准出版社, 2015.
[11]
DALE A, LYTHALL C, AUCOTT J, et al. High precision calorimetry to determine the enthalpy of combustion of methane[J]. Thermochimica Acta, 2002, 382(1/2): 47-54.
[12]
耿维明. 测量误差与不确定度评定[M]. 北京: 中国质检出版社, 2011.