天然气作为目前广泛使用的燃料，其主要成分甲烷与空气的混合气体在密闭球形容器内的爆炸特性，如最大爆炸压力、爆燃指数等是非常重要的物质参数，可以用于风险分析与泄放设计。采用实验方法获得这些数据需要昂贵的实验器材、训练有素的操作员和大量的时间，因此，进行相关的理论和计算研究工作具有重要意义。

国内外研究者已经在爆炸特性方面进行了较多的研究工作。D.Bradley等^[1]建立了同时包含已燃区域和未燃区域的两区域模型，为后来的理论和实验工作奠定了基础；A.E.Dahoe等^[2]研究了密闭容器中的粉尘爆炸，推导了粉尘爆燃指数K_st的公式，对可燃性气体爆燃指数K_G的推导有很大参考意义；C.V.Mashuga等^[3]利用绝热火焰温度对爆炸极限进行预测，为爆炸极限和最大爆炸压力的研究提供了新思路；E.V.D.Bulck^[4]在两区域模型基础上推导出了爆燃指数的计算公式，该公式形式简单，并明确给出爆燃指数与火焰燃烧速度联系；Y.D.Jo等^[5-6]建立了火焰增长模型，并将其用于H₂和空气混合气体压力时程曲线的计算，发现计算结果与实验符合较好，进而推导出爆燃指数的表达式；G.Tao等^[7]利用化学平衡计算方法对H₂/O₂/N₂混合气体的爆炸压力进行了估算，通过与实验结果的对比发现，该结果具有较好的预测效果。

本文中，首先使用化学平衡方法计算得到最大爆炸压力，然后基于火焰增长模型，根据收集的参数数据，用Tablecurve2D拟合经验公式，并利用MATLAB编写程序，通过理论计算得到甲烷-空气混合气体的压力时程曲线和爆燃指数，以期为最大爆炸压力及压力时程曲线的快速预测提供简便方法。

1.   化学平衡计算

化学平衡计算方法是近年来出现的预测爆炸极限和最大爆炸压力的方法，在氢气研究中已经取得较好的效果。爆炸平衡参数可以利用STANJAN、GASEQ等化学平衡计算软件得到，假设燃烧满足绝热和固定体积条件，通过计算可得平衡态时的压力、温度和组分。本文中使用D.S.Dandy教授制作的在线软件进行化学平衡计算^[8]，设初始压力为101325 Pa，初始温度为298 K。在该条件下，估算得到的平衡温度为1450 K，平衡压力为689.48 kPa，接近实际的最大爆炸压力。J.G.Du等^[9]指出1450 K是可燃性混合物的合适限值，即绝热火焰温度高于该值时可以认为是可燃的，低于该值则认为是不可燃的，与我们的估算相符。在一定范围内，平衡压力估算值的改变不会对计算结果造成明显的影响。

计算最大的难点在于产物的确定，通过反复实验消除振荡影响后发现，反应后12种产物包括CH₄、N₂、O₂、C、H、O、N、OH、CO、CO₂、H₂O、H₂，其中C和H₂对计算结果有较大的影响，必须列入产物中，而CH₂、CH₂OH、CH₃O、H₂O₂、NH等物质的含量极低，只有10^-9量级，将其加入到产物中或去除对结果影响较小，甚至没有影响。空气线上化学平衡计算结果与实验值的比较如图 1所示，其中“p”为压力，“α”为甲烷和空气物质的量之比(以下称为空气燃料比)。从图 1中可以看出：实验结果在当量比附近(相当于在α=1.0附近取值)与计算结果符合较好；而在上、下限附近，两者差别较大。这主要是因为当量比附近化学反应速度快，火焰传播迅速，反应完全，与估计的产物相符程度较高；而在上、下限部分因为燃烧不充分，燃料没有燃完，反应速率慢，化学产物成分复杂，不易估计。另一重要的原因是重力对火焰传播和压力的影响，在不考虑重力或微重力情况下，对所有可燃预混气体，火焰都能基本保持球形对称向外传播，在上、下限附近火焰微弱，火焰形状容易受到影响。在实际重力条件下，该情况会更加明显，当浮力引起的气流上升速度大于燃烧速度时，火焰只能向上传播，造成容器下半部分的燃料无法燃烧，造成压力比理想状态低很多，成分复杂。K.L.Cashdollar等^[10]也得到了相似的结果。

图  1  化学平衡计算所得爆炸压力值与实验的对比

Figure  1.  Comparison of chemical equilibrium calculation results with experimental data

下载: 全尺寸图片 幻灯片

2.   参数确定及压力时程曲线的理论计算

2.1   火焰增长模型

火焰增长模型是一种描述球形容器内的爆炸过程的理论。该模型假设在中心点火后，火焰前沿以球形向外扩展，并且以一个很小的微元步长一层一层向外反应，直至火焰前沿达到容器壁。与之前模型不同的是，该模型区分已燃部分和未燃部分的绝热指数，火焰增长情况更加接近实际燃烧过程，预测的压力值也更加准确。火焰增长模型的假设如下：(1)火焰前沿为球形；(2)气体压力均一；(3)气体浮力和火焰厚度可以忽略。基于以上假设，火焰增长模型中压力、燃烧爆炸程度、任意时刻半径可由以下方程得到，即：

式中：t为时间；t为量纲一时间，且t=tS_u/r_a，其中S_u为层流燃烧速率，r_a为罐的半径；ξ为爆炸程度；γ_u和γ_b分别为未燃部分和已燃部分的绝热指数；r_b=r_b/r_a为已燃部分的量纲一半径，r_b为已燃部分的半径；p=p/p₀为量纲一压力，其中p为某一时刻压力，p₀为初始压力；p_E为固定体积情况下的平衡爆炸压力；n为燃烧速率的指数参数。初始条件取为：t=0，p=1，ξ=0，r_b=ε₀。其中ε₀为非常小的正值，可以根据点燃器的尺寸确定。

对于一个给定浓度的甲烷空气混合气体，计算p和t曲线的过程如下：首先用化学平衡计算方法和拟合公式计算得到量纲一的平衡爆炸压力和绝热指数，然后将初始状态代入式(1)~式(3)进行迭代计算，过程如下：将量纲一半径代入式(1)，计算得到爆炸程度；将所得爆炸程度代入式(2)，计算得到量纲一爆炸压力；将量纲一爆炸压力代入式(3)，得到量纲一半径；再次将量纲一半径代入式(1)，如此循环计算。计算过程根据步长逐步推进，当爆炸程度达到单位1时，计算完成并输出结果。

2.2   层流燃烧速率的确定

图  2  层流燃烧速率与空气燃料比拟合曲线

Figure  2.  Fitted curve of laminar flame speed vs. air fuel ratio

下载: 全尺寸图片 幻灯片

火焰增长模型中，层流火焰速度S_u的定义为平坦的、无拉伸的、绝热的、预混的火焰前锋传播到未燃气体时的火焰速度。层流火焰速度是燃烧混合物的重要参数，决定着燃烧程度和实际燃烧过程中火焰的稳定性。层流火焰速度是化学机理研究的重要对象，许多科学家对甲烷/空气混合物的层流燃烧速率进行过测量。本文中对有代表性的数据进行分析^[11-14]，并通过数值拟合得出层流火焰速度与空气燃料比α之间的关系。

为使拟合公式具有工程应用价值，本文使用三次多项式进行拟合，所得表达式为y=A+B₁x+B₂x²+B₃x³，参数的取值如下：A=-183.12，B₁=448.52，B₂=-256.91，B₃=27.58。决定系数R²可以达到0.96954，即拟合效果较好。

2.3   绝热指数的确定

绝热指数是求解气体膨胀程度的重要参数，在火焰增长模型中，该参数在爆炸程度和压力时程曲线计算中有很大作用。其他模型将已燃部分和未燃部分的绝热指数取为相等的值，以简化计算，与实际情况有较大的不同，会造成计算上的误差。E.V.D.Bulck^[15]通过实验方法得到了甲烷已燃和未燃部分绝热指数数据。基于文献[15]数据，本文利用Tablecurve 2D软件进行数值拟合，并从Tablecurve2D上千个拟合多项式中挑选出形式简单且拟合效果较好的结果，如图 3和图 4所示。从图 3和图 4可以看出；γ_u和γ_b拟合曲线的决定系数均非常接近1，说明拟合效果很好。

图  3  未燃部分绝热指数拟合曲线

Figure  3.  Fitted curve of specific heat ratio for unburned gas

下载: 全尺寸图片 幻灯片

图  4  已燃部分绝热指数拟合曲线

Figure  4.  Fitted curve of specific heat ratio for unburned gas

下载: 全尺寸图片 幻灯片

2.4   理论计算与比较

取空气燃料比为1.055的情况进行理论计算，根据拟合公式计算可得：已燃部分的绝热指数约为1.062，未燃部分的绝热指数为1.374，层流火焰速度根据拟合公式计算为36.51 cm/s；利用化学平衡计算估计可得：绝热状态下最大压力为912.87 kPa，对应的量纲一压力为9.01；实验时每隔0.333 ms记录一个数据点，因此理论计算时间间隔取为0.333 ms，对应的量纲一时间为0.0074；实验过程中点火方式为中心电火花点火，可以认为点火区域集中于一个半径很小的球形中，所以理论计算时点火半径取为1 cm，对应的量纲一半径为0.06；甲烷的指数参数根据文献[5-6]中数据取为0.1。用MATLAB进行编程，所有的初始值和需要提前计算的数据均已给出，将其代入式(1)~式(3)进行循环运算，直到爆炸程度达到单位1。

图 5给出了基于火焰增长模型计算结果与实验数据的比较。由图 5可知，当爆炸程度达到1时，最大压力为902.73 kPa，与用化学平衡计算得到的爆炸压力912.87 kPa基本一致，也出现略高于实验值的情况，但更为精确。从图 3还可以看出，当爆炸程度达到1时，量纲一半径已经略微超过1。这主要是因为在用MATLAB进行迭代计算时必须给定一个很小的量纲一半径初始值，而爆炸程度初始值为零，虽然初始值较小，但经过接近300步的迭代和求根计算累积的误差，造成量纲一半径和爆炸程度不可能同时达到单位1；另一方面，火焰前沿接触到容器壁时会出现拉伸、反射、压缩的现象，会造成压力进一步的升高。总体来说，甲烷-空气混合物在空气燃料比为1.055时，理论计算得到的压力时程曲线与实验曲线的拟合效果较好，但略高于实验值。这主要是因为理论计算假设是理想状态，不考虑热量散失、浮力等造成的。与Y.D.Jo等^[5-6]所得结果相比，本文采用化学平衡计算预测最大压力，假设平衡温度为1450 K，参与反应过程元素假设与其并不完全相同，部分元素对结果几乎没有影响。对于层流燃烧速率、已燃部分绝热指数、未燃部分绝热指数，本文总结了较多数据并进行了拟合，有更好的工程应用价值。此外，通过循环计算也能得到爆炸程度与时间的关系、爆炸压力上升速度与时间的关系。

图  5  火焰增长模型计算结果与实验值的比较

Figure  5.  Comparison of theoretical results by flame growth model with experimental data

下载: 全尺寸图片 幻灯片

3.   爆燃指数的计算

应用火焰增长模型还可计算爆燃指数。爆燃指数是估计泄放面积、预防事故发生方面的重要参数，火焰增长模型推导出的爆燃指数计算公式如下：

式(4)假设火焰前沿接触爆炸容器球形表面时的状态已知，未考虑爆炸结束段火焰表面积的变化。由于三次方定律，爆燃指数与体积无关，但事实上爆炸釜体积的变化会对爆燃指数产生影响。当火焰前沿厚度不可忽略时，假设就会失效。

表 1给出了不同爆炸程度条件下球形容器的爆燃指数，其中K_G, exp为爆燃指数的实验值。由表 1可知：当ξ=1时，利用式(4)计算得到的爆燃指数K_G=8.24091 MPa·m·s^-1，比实验值高出较多；当爆炸程度为0.95和0.90时，计算结果与实验值更加接近。上述结果表明，一部分的火焰前沿在爆炸完成(即爆炸程度ξ=1)前就已经接触到了球形爆炸容器的表面，造成火焰表面积有很大程度的变化，对爆炸压力上升速度有较大影响，从而影响爆燃指数值。因此，爆炸程度ξ=0.90时，爆燃指数的计算值与实验值符合更好。

表  1  球形容器内甲烷-空气混合气体的爆燃指数

Table  1.  Deflagration index of methane-air mixture in spherical vessel

空气燃料比	KG, exp/(MPa·m·s-1)	KG/(MPa·m·s-1)
		ξ<1	ξ<0.95	ξ<0.90
0.828	2.136	4.44118	4.21257	3.88981
0.945	6.305	7.50427	7.00937	6.40062
1.055	7.417	8.24091	7.87229	7.17593
1.058	7.365	8.39706	7.84194	7.14657
1.184	5.910	7.66787	7.17321	6.53207

下载: 导出CSV 

| 显示表格

4.   结论

(1) 利用化学平衡计算方法计算甲烷-空气最大爆炸压力时，在当量比附近有较好的预测结果，在上、下限附近，由于燃烧机理的复杂性、燃烧速度较低，造成预测结果不理想。

(2) 在当量比附近，应用火焰增长模型得到的压力时程曲线与实验数据符合性较好。

(3) 应用火焰增长模型计算得到的爆燃指数K_G与实验值符合性较好，并且使用爆炸程度为0.9时的压力计算得到的爆燃指数更加准确，说明在釜中燃料燃烧结束前，爆炸压力上升速度已达最大值。

感谢江苏研究生科研创新计划项目及Michigan Technological University的资助，同时感谢Daniel A.Crowl教授的帮助。