Comparison of different methods for source terms in detonation simulation
-
摘要: 为解决爆震燃烧模拟中出现的刚性问题，对处理刚性源项问题常见的一步法、逼近法、拟稳态逼近法(α qusai steady state, αQSS)和点隐方法进行对比，从稳定性等方面分析源项处理方法应满足的时间步长要求，并探索各方法之间的联系以及适应化学反应特征变化的能力，进一步通过球头激波诱导燃烧算例比较每种方法的计算效率。理论分析和数值计算表明：一步法在积分刚性源项时，积分步长需小于或等于2倍最短反应特征时间，而逼近法、αQSS法和点隐方法对时间步长取值没有影响；αQSS法可根据化学反应特征的变化自动调整α值和时间步长，适用范围较广，而一步法和逼近法则是αQSS方法的特例。点隐等隐式方法在求解数学意义上的刚性问题时稳定性很好，但计算效率较低。相比而言，αQSS法在计算稳定性和适应化学反应变化方面都具有良好的性能，且针对激波诱导燃烧算例，αQSS法消耗的CPU时间仅为点隐方法的一半，是处理刚性源项较好的选择。Abstract: In this study, to solve the stiff source terms resulting from chemical reactions in detonation simulation, we examined the one step method, the asymptotic approach, the α quasi steady state method (αQSS) and the point implicit and compared their performances in coping with the stiff source term problems. We studied the limitations of each method using stability analysis, and investigated their relationships and capabilities in adapting to the changes in chemical reactions, with the shock-induced combustion simulated to compare their efficiencies. The results indicate that the one step method requires at least two times of the smallest time scale while the other three methods have no constraint on the time step. The αQSS can adjust the value of α and the time step with different reaction characteristics, and the one step method and the asymptotic method are the special cases of the αQSS with a constant α. An implicit approach has a better performance in mathematically solving the stiff equations but its low computation efficiency from the inversion of the matrix is sometimes unacceptable. The αQSS method can only consume a half of the CPU time that with the point implicit in shock-induced combustion simulation. In general, the αQSS is a good choice for dealing with stiff source term problems.
-
Key words:
- stiff source /
- stability analysis /
- α quasi steady state /
- point implicit method /
- asymptotic method /
- one step method
-
随着工业的快速发展，可燃气体、粉尘共存的工况不断出现，如煤炭开采中的瓦斯与煤尘、制药行业的溶剂蒸气与药粉、印刷行业的溶剂蒸气与色素粉、塑料行业的乙烯与聚乙烯粉等[1-2]。气/粉两相混合体系爆炸事故已在工业粉尘爆炸事故总数中占相当比例[3]。因此，气/粉两相体系爆炸特性逐渐引起研究人员的关注。Cashdollar等[4]、Dufaud等[5]、Sanchirico等[6]通过实验发现体积分数低于单相气体爆炸下限的可燃气体和可燃粉尘混合后可以被点燃，仍具有爆炸的危险性。Bartknecht等[7]开展实验探索甲烷、丙烷和丁烷气体对纤维素粉爆炸强度作用规律，结果表明3种烷烃气体均能引起粉尘最大爆炸压力的上升。Denkevits等[8]通过实验发现氢气/铝粉两相爆炸压力及压力上升速率与单相爆炸相比都显著提升。Amyotte等[9-10]实验研究了乙烯/聚乙烯等两相体系的爆炸特性，得出乙烯的添加能够引起聚乙烯粉尘爆炸指数显著提升。Liu等[11]利用3.2 L长方形容器研究甲烷/煤粉混合物的火焰传播，得出甲烷和煤粉共存显著加快了爆炸火焰传播速度，并且提高了最高火焰温度。上述研究从宏观上初步明确了可燃气体能够显著提高粉尘爆炸感度、强度等参数，即提高粉尘爆炸危险性。但现有的研究还无法解释气/粉两相体系爆炸参数的变化规律，需要继续深入开展研究。并且已开展的研究多关注低体积分数气体对粉尘爆炸特性的影响，而对高体积分数气体的影响关注较少。基于此，本文中采用改进的20 L球形粉尘爆炸装置，研究相比单相爆炸，气/粉两相体系爆炸强度特性参数的变化规律，以补充和完善气/粉两相爆炸理论及爆炸防治技术。
1. 实验装置
实验基于改进的标准20 L球形粉尘爆炸装置开展。如图 1所示，装置主要由爆炸容器、点火系统、配气系统、测量及数据采集系统等组成。本装置已在多个研究项目中使用，其可靠性已得到验证[12-13]。
爆炸容器为符合国际标准E1226的20 L球形容器，内部布置反弹式喷嘴、点火系统，底部安装气/粉两相阀。在20 L球形容器内使用10 kJ点火头容易出现过驱动效应，导致压力上升速率等测量结果明显偏大[14-15]，且高能量点火头容易掩盖可燃气体对粉尘爆炸特性的影响[16]。因此，采用能量为0.5 kJ的化学点火头进行点火。该点火头引起的压力峰值约为8 kPa，相对较低，在很大程度上避免了点火头对混合体系爆炸压力变化规律的影响，并且该点火头在实验过程中能够实现爆炸介质的稳定点火。
在粉尘仓开口并连接甲烷气瓶使之满足气/粉两相体系测试要求。两相体系的配制流程为：(1)将一定量粉尘(根据粉尘质量浓度计算)置于粉尘仓，并密封。(2)将甲烷注入粉尘仓。注入的量依据道尔顿分压定律并结合压力数值计算。(3)启动真空泵, 将20 L球形容器抽真空。(4)向粉尘仓充入高压空气至2 MPa，甲烷与空气在粉尘仓中实现初次预混。(5)启动两相阀，一定体积分数的甲烷/空气携带粉尘高速喷入容器内，实现甲烷、石松子粉尘、空气的均匀混合。(6)60 ms后，点火头点火点燃两相混合物。
本实验中采用预混气体扬尘的方式配置气/粉混合物，更有利于可燃气体、粉尘、空气的均匀混合。采用高频压力传感器测量球体主压力，其量程为0~2 MPa，采集频率为5 kHz，精度等级为0.25级。
采用分散性和流动性好的石松子粉尘和甲烷气体(纯度99.99%)作为研究对象。石松子粉尘中位直径为39 μm，其显微结构和粒径分布如图 2~3所示。选取质量浓度为50、100、150、200、250、500、750、1 000、1 250 g/m3的石松子粉和体积分数为2%、4%、6%、8%、10%、12%的甲烷进行配比，开展实验。为了保证实验数据的可靠性，每组实验至少重复3次。
图 3 石松子粉尘粒径分布(D[50]=39 μm)Figure 3. Particle size distribution of lycopodium dust with D[50]=39 μm2. 结果与分析
2.1 甲烷对石松子粉尘最大爆炸压力的影响
图 4为甲烷体积分数为0%、2%、4%、6%、8%、10%和12%时，石松子粉尘爆炸压力随粉尘质量浓度的变化。由图 4可知，当粉尘质量浓度为250 g/m3时，向粉尘中添加体积分数为2%、4%、6%和8%的甲烷后，粉尘爆炸压力显著提高。这是因为低质量浓度的石松子粉尘爆炸属于贫燃料燃烧过程，其爆炸压力主要由反应物的数量控制，甲烷的加入：一方面，提高了反应物的数量；另一方面，加速粉尘分解，诱导粉尘燃烧更完全，释放出更多热量。当粉尘质量浓度升高至750~1 250 g/m3时，加入甲烷后，粉尘爆炸压力降低。这是因为高质量浓度石松子粉尘爆炸属于富燃料燃烧，其爆炸压力主要由体系中的氧气控制，甲烷的加入不仅提高了体系中的燃料含量，还降低了体系中的氧气含量，诱导体系中可燃介质的不完全燃烧，使燃烧释放出的热量减少。
图 4 甲烷体积分数为0%、2%、4%、6%、8%、10%和12%时石松子粉尘爆炸压力Figure 4. Explosion pressure of lycopodium dust at the methane mass fractions of 0%, 2%, 4%, 6%, 8%, 10% and 12%图 5为两相体系爆炸压力(平均值)随甲烷体积分数的变化。图中空心数据点为不同体积分数纯甲烷气体的爆炸压力，半实心圆形图标代表各甲烷体积分数对应的甲烷/石松子粉两相混合体系的最大爆炸压力，灰色区域上边界线对应的压力值为单相石松子粉尘最大爆炸压力。由图 5可知，单相石松子粉最大爆炸压力为681.2 kPa，对应的最佳爆炸质量浓度为750 g/m3。当甲烷体积分数分别为2%、4%、6%时，两相体系最大爆炸压力分别为632.4、656.4、679.6 kPa，对应的最佳爆炸质量浓度分别为750、500和250 g/m3。因此，两相体系的最大爆炸压力并未随甲烷含量的增加而显著提高。这是因为，对某一粉尘，其最大爆炸压力的定义为所有质量浓度下爆炸压力的最大值。尽管甲烷的加入能够提高低质量浓度粉尘的爆炸压力，但也会降低高质量浓度粉尘的爆炸压力，在所有质量浓度范围内，石松子粉尘最大爆炸压力受甲烷影响并不明显。当甲烷体积分数超过8%时, 如图 4(b)所示，两相体系的最大爆炸压力为单相甲烷爆炸压力，对应的最佳爆炸质量浓度为0 g/m3。此时，体系中的甲烷体积分数接近或大于当量体积分数(9.5%)，该体系中的石松子粉尘更多的是作为吸热体来吸收甲烷燃烧放出的热量，因此向该体系中添加任何质量浓度的石松子都将降低体系爆炸压力。由图 5还可以发现，单相石松子及两相混合体系的最大爆炸压力均小于单相甲烷的最大爆炸压力。
图 5 甲烷/石松子混合体系爆炸压力Figure 5. Explosion pressure of hybrid mixtures of methane and lycopodium dust由前述可知，甲烷对石松子粉尘最佳爆炸质量浓度有显著影响。基于不同甲烷体积分数条件下的石松子粉尘最佳爆炸质量浓度，得出石松子粉尘最佳爆炸质量浓度ρo随甲烷体积分数的变化规律如图 6所示。由图 6可知，体系中加入体积分数为2%的甲烷时，石松子粉尘最佳爆炸质量浓度保持在750 g/m3，即低体积分数甲烷对石松子粉尘最佳爆炸质量浓度的影响并不明显。随着甲烷体积分数进一步提高，石松子粉尘最佳爆炸质量浓度明显降低。因此，甲烷能够显著提高粉尘的爆炸敏感性和危险性。
图 6 甲烷/石松子混合体系的最佳爆炸质量浓度随甲烷体积分数的变化Figure 6. Optimum mass concentration of methane/lycopodium mixtures varying with methane volume fraction2.2 甲烷对石松子粉尘最大爆炸升压速率的影响
图 7为石松子粉尘最大爆炸升压速率随甲烷体积分数的变化规律，图中数据均为3次实验的平均值。由图 7可以看出，当甲烷体积分数为0%时，单相石松子粉尘最大爆炸升压速率(dp/dt)max为18.81 MPa/s。石松子粉尘最大爆炸升压速率随甲烷体积分数的增大呈现升高趋势。当甲烷体积分数为6%且石松子粉尘质量浓度为250 g/m3时，甲烷/石松子体系的最大爆炸升压速率达到极值，高达28.21 MPa/s，与单相石松子粉尘相比提升了50%。当甲烷体积分数在8%~12%范围内时，体系最大爆炸升压速率均为单相甲烷对应的爆炸升压速率(101.08、126.94和71.7 MPa/s)，远大于两相体系最大爆炸升压速率。即在该甲烷体积分数范围内，石松子粉尘作为吸热体存在，向甲烷中添加任何质量浓度的石松子粉尘均将诱导体系最大爆炸升压速率降低。
图 7 石松子粉尘最大爆炸升压速率随甲烷体积分数的变化Figure 7. Maximum explosion pressure rise rate of lycopodium dust at different methane volume fractions石松子粉尘燃烧需要经过吸热、分解、与空气混合、引燃、火焰传播等过程，甲烷的加入：一方面，能够改变石松子粉尘燃烧过程中的限速反应，进而改变粉尘燃烧的动力学特征，加快粉尘燃烧速度；另一方面，气相甲烷与体系中的氧气首先发生氧化反应并释放出热量，该热量作用于体系中石松子粉尘，加速石松子粉尘热解，进而提高粉尘燃烧速度，燃烧速度的提高导致升压速率的升高。因此，甲烷的加入导致石松子粉尘爆炸升压速率明显升高，且该升高幅度随甲烷含量的增加而增大。
2.3 甲烷/石松子粉两相体系爆炸强度评估
爆炸指数Kst是直接评估粉尘爆炸强度的另一个重要参数。根据NFPA68[17]和ASTM E1226[18]标准，可通过立方根定律计算出甲烷、石松子粉尘及甲烷/石松子粉混合体系爆炸指数，结果见表 1，表中数据均为平均值。由表 1可知，两相体系爆炸指数随甲烷含量的增加而增大，当甲烷体积分数φ为6%且石松子粉尘质量浓度ρ为250 g/m3时，甲烷/石松子粉两相体系爆炸指数最高，为7.66 MPa·m·s-1。
表 1 甲烷、石松子粉尘及甲烷/石松子混合体系爆燃指数Table 1. Explosive deflagration indices of methane, lycopodium and methane/lycopodium mixturesφ/% ρ/(g·m-3) (dp/dt)max/(MPa·s-1) Kst(KG)/(MPa·m·s-1) 0(pure dust) 750 18.81 5.11 2 750 19.98 5.42 4 500 21.16 5.74 6 250 28.21 7.66 10(pure CH4) 0 126.94 34.46 结合前述分析可知：单相石松子粉尘和甲烷/石松子两相混合体系的最大爆炸压力、最大爆炸升压速率、爆炸指数均低于单相甲烷的；3种可爆体系爆炸强度高低顺序为：单相甲烷爆炸强度＞甲烷/石松子混合体系爆炸强度＞单相石松子粉尘爆炸强度。因此，工业生产过程中应尽量避免可燃粉体中混入可燃气体，以降低粉尘爆炸危险性。
3. 结论
基于改进的标准20 L球形爆炸装置，系统实验研究了甲烷/石松子粉尘两相体系爆炸强度的变化规律。根据实验结果，得到以下结论：(1)甲烷的添加能够显著提高低质量浓度石松子粉尘的爆炸压力而降低高质量浓度石松子粉尘的爆炸压力，但甲烷对石松子粉尘的最大爆炸压力没有显著影响。(2)石松子粉尘的最大爆炸升压速率却随甲烷含量的增加而升高，最佳爆炸质量浓度随甲烷含量的增加而显著降低，即甲烷诱导石松子粉尘爆炸强度和敏感性增大。(3)石松子、甲烷/石松子和甲烷3种体系的爆炸强度满足如下关系：单相甲烷爆炸强度＞甲烷/石松子混合体系爆炸强度＞单相石松子粉尘爆炸强度。因此，工业生产过程中应尽量避免可燃粉体中混入可燃气体，以降低粉尘爆炸危险性。
-
图 2 一步法处理组分质量分数随轴线分布曲线
Figure 2. Mass fraction distribution of species along the axis by one step method
3(a) αQSS方法处理组分质量分数随轴线分布曲线
3(a). Mass fraction distribution of species along the axis by αQSS method
3(b) 点隐法处理组分质量分数随轴线分布曲线
3(b). Mass fraction distribution of species along the axis by point implicit method
3(c) 逼近法处理组分质量分数随轴线分布曲线
3(c). Mass fraction distribution of species along the axis by asymptotic method
表 1 实验数据与理论模型结果对比
Table 1. Comparison of experimental and theoretical data
方法 f/Hz Ma=4.48 Ma=4.79 实验 425.0 712.0 αQSS 424.0 721.5 点隐法 426.1 725.7 一步法 424.1 721.5 逼近法 425.2 721.5 
下载: 导出CSV
表 2 不同源项处理方法消耗CPU时间
Table 2. The CPU time in different stiff source term methods
方法 CPU time/s αQSS 77 590 点隐法 130 007 一步法 154 774 逼近法 62 372
下载: 导出CSV
-
[1] 刘君, 周松柏, 徐春光.超声速流动中燃烧现象的数值模拟方法及应用[M].长沙:国防科技大学出版社, 2008:76-79. [2] BUSSING T R A. A finite volume method for the Navier-Stokes equations with finite rate chemistry[D]. Cambridge: Massachusetts Institute of Technology, 1985. http: //ci. nii. ac. jp/ncid/BB03917781 [3] BUSSING T R A, Murman E M. Finite-volume method for the calculation of compressible chemically reacting flows[J]. AIAA Journal, 1988, 26(9):1070-1078. doi: 10.2514/3.10013 [4] ZHONG X. Additive semi-implicit Runge-Kutta methods for computing high-speed non-equilibrium reactive flows[J]. Journal of Computational Physics, 1996, 128(1):19-31. doi: 10.1006/jcph.1996.0193 [5] ROSENBROCK H H. Some general implicit processes for the numerical solution of differential equations[J]. The Computer Journal, 1963, 5(4):329-330. doi: 10.1093/comjnl/5.4.329 [6] LEVEQUE R J, YEE H C. A study of numerical methods for hyperbolic conservation laws with stiff source terms[J]. Journal of Computational Physics, 1990, 86(1):187-210. doi: 10.1016/0021-9991(90)90097-K [7] ORAN E S, BORIS J P. Numerical simulation of reactive flow[M]. Cambridge: Cambridge University Press, 2005:114-158. [8] YOUNG T R, BORIS J P. A numerical technique for solving stiff ordinary differential equations associated with the chemical kinetics of reactive-flow problems[J]. The Journal of Physical Chemistry, 1977, 81(25):2424-2427. doi: 10.1021/j100540a018 [9] CHIANG T, HOFFMANN K. Determination of computational time step for chemically reacting flows[C]//Proceedings of AIAA 20th Fluid Dynamics, Plasma Dynamics and Laser Conference. Buffalo, New York, USA, 1989. [10] MOTT D R, ORAN E S, VAN L B. A quasi-steady-state solver for the stiff ordinary differential equations of reaction kinetics[J]. bJournal of Computational Physics, 2000, 164(2):407-428. doi: 10.1006/jcph.2000.6605 [11] MOTT D R, ORAN E S. CHEMEQ2: A solver for the stiff ordinary differential equations of chemical kinetics[R]. Naval Research Lab, Washington D C, 2001. [12] 刘瑜. 化学非平衡流的计算方法研究及其在激波诱导燃烧现象模拟中的应用[D]. 长沙: 国防科技大学, 2008. http: //cdmd. cnki. com. cn/Article/CDMD-90002-2009213254. htmLIU Yu. Investigations into numerical methods of chemical non-equilibrium flow and its application to simulation of shock-induced combustion phenomenon[D]. Changsha: National University of Defense Technology, 2008. http: //cdmd. cnki. com. cn/Article/CDMD-90002-2009213254. htm [13] 刘君, 张涵信, 高树椿.一种新型的计算化学非平衡流动的解耦方法[J].国防科技大学学报, 2000, 22(5):19-23. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=gfkjdxxb200005005LIU Jun, ZHANG Hanxin, GAO Shuchun. A new uncoupled method for numerical simulation of non-equilibrium flow[J]. Journal of National University of Defense Technology, 2000, 22(5):19-23. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=gfkjdxxb200005005 [14] 刘世杰, 林志勇, 孙明波, 等.采用不同化学反应源项处理方法的胞格爆轰数值研究[J].国防科技大学学报, 2010, 32(5):1-6. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=gfkjdxxb201005001LIU Shijie, LIN Zhiyong, SUN Mingbo, et al. Numerical simulation of cellular detonation using different chemical reaction source term methods[J]. Journal of National University of Defense Technology, 2010, 32(5):1-6. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=gfkjdxxb201005001 [15] JACHIMOWSKI C J. An analytical study of the hydrogen-air reaction mechanism with application to scramjet combustion[R]. NASA TechnicaI Paper 2791, 1988. [16] 潘沙. 高超声速气动热数值模拟方法及大规模并行计算研究[D]. 长沙: 国防科技大学, 2010. http: //cdmd. cnki. com. cn/Article/CDMD-90002-2010271147. htmPAN Sha. Hypersonic aerothermal numerical simulation method and massive parallel computation research[D]. Changsha: National University of Defense Technology, 2010. http: //cdmd. cnki. com. cn/Article/CDMD-90002-2010271147. htm [17] TORO E F. Riemann solvers and numerical methods for fluid dynamics: A practical introduction[M]. Springer Science & Business Media, 2009: 531-542. [18] LEHR H F. Experiments on shock-induced combustion[J]. Astronautica Acta, 1972, 17:589-597. https://www.researchgate.net/publication/279938305_Experiments_on_shock-induced_combustion [19] MCVEY J B, TOONG T Y. Mechanism of instabilities of exothermic hypersonic blunt-body flows[J]. Combustion Science and Technology, 1971, 3(2):63-76. doi: 10.1080/00102207108952273 [20] CHOI J Y, JEUNG I S, YOON Y. Computational fluid dynamics algorithms for unsteady shock-induced combustion: Part 1: Validation[J]. AIAA Journal, 2000, 38(7):1179-1187. doi: 10.2514/2.1112 [21] 刘世杰, 孙明波, 林志勇, 等.钝头体激波诱导振荡燃烧现象的数值模拟[J].力学学报, 2010, 42(4):597-606. doi: 10.6052/0459-1879-2010-4-lxxb2010-086LIU Shijie J, SUN Mingbo, LIN Zhiyong, et al. Numerical research on blunt body shock-induced oscillation combustion phenomenon[J]. Chinese Journal of Theoretical and Applied Mechanics, 2010, 42(4):597-606. doi: 10.6052/0459-1879-2010-4-lxxb2010-086 期刊类型引用(0)
其他类型引用(1)
-
下载:
计量
- 文章访问数: 5282
- HTML全文浏览量: 1730
- PDF下载量: 175
- 被引次数: 1