陆胤臣，陶 刚，张礼敬

（南京工业大学安全科学与工程学院，江苏南京210009）

球形容器内甲烷－空气爆炸特性分析与理论计算＊

陆胤臣，陶 刚，张礼敬

（南京工业大学安全科学与工程学院，江苏南京210009）

为研究甲烷－空气混合物在密闭球形容器内的爆炸特性，首先利用化学平衡计算软件确定合适的燃烧产物与化学平衡温度，估算甲烷－空气混合物的最大爆炸压力。然后基于火焰增长模型，用MATLAB编辑循环语句程序，计算了甲烷－空气混合气体爆炸的压力时程曲线，通过与实验数据对比，验证了化学平衡软件计算方法与火焰增长模型的可行性，并分析了误差的产生原因。进而利用火焰增长模型推出的经验公式计算爆燃指数，发现在当量比附近与实验结果拟合程度较好。

最大爆炸压力；爆燃指数；化学平衡计算；火焰增长模型

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

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

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

1 化学平衡计算

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

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

图1 化学平衡计算所得爆炸压力值与实验的对比Fig．1 Comparison of chemical equilibrium calculation results with experimental data

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

2．1 火焰增长模型

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

2．2 层流燃烧速率的确定

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

为使拟合公式具有工程应用价值，本文使用三次多项式进行拟合，所得表达式为y＝A＋B1x＋B2x2＋B3x3，参数的取值如下：A＝－183．12，B1＝448．52，B2＝－256．91，B3＝27．58。决定系数R2可以达到0．96954，即拟合效果较好。

2．3 绝热指数的确定

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

图2 层流燃烧速率与空气燃料比拟合曲线Fig．2 Fitted curve of laminar flame speed vs．air fuel ratio

图3 未燃部分绝热指数拟合曲线Fig．3 Fitted curve of specific heat ratio for unburned gas

图4 已燃部分绝热指数拟合曲线Fig．4 Fitted curve of specific heat ratio for unburned gas

2．4 理论计算与比较

取空气燃料比为1．055的情况进行理论计算，根据拟合公式计算可得：已燃部分的绝热指数约为1．062，未燃部分的绝热指数为1．374，层流火焰速度根据拟合公式计算为36．51cm／s；利用化学平衡计算估计可得：绝热状态下最大压力为912．87kPa，对应的量纲一压力为9．01；实验时每隔0．333ms记录一个数据点，因此理论计算时间间隔取为0．333ms，对应的量纲一时间为0．007 4；实验过程中点火方式为中心电火花点火，可以认为点火区域集中于一个半径很小的球形中，所以理论计算时点火半径取为1 cm，对应的量纲一半径为0．06；甲烷的指数参数根据文献［5－6］中数据取为0．1。用MATLAB进行编程，所有的初始值和需要提前计算的数据均已给出，将其代入式（1）～式（3）进行循环运算，直到爆炸程度达到单位1。

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

图5 火焰增长模型计算结果与实验值的比较Fig．5 Comparison of theoretical results by flame growth model with experimental data

3 爆燃指数的计算

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

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

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

表1 球形容器内甲烷－空气混合气体的爆燃指数Table 1 Deflagration index of methane－air mixture in spherical vessel

4 结 论

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

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

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

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

［1］ Bradley D，Mitcheson A．Mathematical solutions for explosions in spherical vessels［J］．Combustion and Flame，1976，26（2）：201－217．

［2］ Dahoe A E，Zevenbergen J F，Lemkowitz S M，et al．Dust explosion in spherical vessels：The role of flame thickness in the validity of the cube－root law［J］．Journal of Loss Prevention in the Process Industries，1996，18（9）：33－44．

［3］ Mashuga C V，Crowl D A．Flammability zone prediction using calculated adiabatic flame temperatures［J］．Process Safety Progress，1999，18（3）：127－134．

［4］ Bulck E V D．Closed algebraic expressions for the adiabatic limit value of the explosion constant in closed volume combustion［J］．Journal of Loss Prevention in the Process Industries，2005，18（1）：35－42．

［5］ Jo Y D，Crowl D A．Flame growth model for confined gas explosion［J］．Process Safety Progress，2009，28（2）：141－146．

［6］ Jo Y D，Crowl D A．Explosion characteristics of hydrogen－air mixtures in a spherical vessel［J］．Process Safety Progress，2010，29（3）：216－223．

［7］ Tao G，Crowl D A．Comparison of the maximum gas combustion pressure of hydrogen／oxygen／nitrogen between chemical equilibrium calculations and experimental data［J］．Procedia Engineering，2013，62：786－790．

［8］ Dandy D S．Bioanalytical microfluidics program［EB／OL］．［2015－12－28］．http：∥navier．engr．colostate．edu／tools／equil．html．

［9］ Du J G，Ma H H，Qu Z W，et al．Prediction of methanes flammability using chemical equilibrium［J］．Process Safety Progress，2015，34（1）：31－35．

［10］ Cashdollar K L，Zlochower I A，Green G M，et al．Flammability of methane，propane，and hydrogen gases［J］．Journal of Loss Prevention in the Process Industries，2000，13（3）：327－340．

［11］ Van Maaren A，Thung D S，De Goey L P H．Measurement of flame temperature and adiabatic burning velocity of methane／air mixtures［J］．Combustion Science and Technology，1994，96（4／5／6）：327－344．

［12］ Dahoe A E，De Goey L P H．On the determination of the laminar burning velocity from closed vessel gas explosions［J］．Journal of Loss Prevention in the Process Industries，2003，16（6）：457－478．

［13］ Varea E，Modica V，Vandel A，et al．Measurement of laminar burning velocity and Markstein length relative to fresh gases using a new postprocessing procedure：Application to laminar spherical flames for methane，ethanol and isooctane／air mixtures［J］．Combustion and Flame，2012，159（2）：577－590．

［14］ Chen Z．On the accuracy of laminar flame speeds measured from outwardly propagating spherical flames：Methane／air at normal temperature and pressure［J］．Combustion and Flame，2015，162（6）：2442－2453．

［15］ Bulck E V D．Closed algebraic expressions for the adiabatic limit value of the explosion constant in closed volume combustion［J］．Journal of Loss Prevention in the Process Industries，2005，18（1）：35－42．

［16］ Benedetto A D，Cammarota F，Sarli V D，et al．Anomalous behavior during explosions of CH4in oxygen－enriched air［J］．Combustion and Flame，2011，158（11）：2214－2219．

Analysis and theoretical calculation of explosion characteristics of methane－air mixture in a spherical vessel

Lu Yinchen，Tao Gang，Zhang Lijing

（College of Safety Science and Engineering，Nanjing Tech University，Nanjing210009，Jiangsu，China）

To study the characteristics of the methane－air mixture exploding in a closed spherical container，we determined the appropriate combustion products and chemical equilibrium temperature using the chemical equilibrium calculation software，thereby predicting the maximum explosion pressure of the mixture．The MATLAB program based on the flame growth model was adopted to calculate the curve showing the relationship between the explosion pressure and time．The calculation processes were verified by the comparison of the obtained results with the experimental data，and the origin of the error was also identified．Further，it is found that the empirical formula of the deflagration index KGderived from the flame growth model is well fitted with the experimental date near the chemical equivalent line．

maximum pressure；deflagration index；chemical equilibrium calculation；flame growth model

O354．1国标学科代码：1303510

A

10．11883／1001－1455（2017）04－0773－06

（责任编辑 王玉锋）

2015－12－28；

2016－03－30

江苏省2015年度普通高校研究生科研创新计划项目（KYLX15＿0774）

陆胤臣（1991－ ），男，硕士研究生；通信作者：陶 刚，taogang＠njtech．edu．cn。