煤堆非均质分布及其自热研究

李 坤 林 琦 杨红阳

（安徽工业大学，安徽 马鞍山 243002）

自燃是自然界中较为常见的灾害现象，给人类的生产、生活、健康，乃至大气环境等方面均带来了极大的危害[1-2]。仅以煤炭行业为例，我国近十年发生的195 起煤矿火灾事故中九成以上是由自燃引起的，造成了共计200 余人死亡，损毁的煤炭资源高达千万吨[3-4]。

自燃具有潜伏期长、隐蔽性强、熄灭困难等特点，在潜伏期内进行人工干预的效果明显优于自燃发生后的灭火治理，故对潜伏期内自热演化过程的准确掌握尤为重要[5]，特别是谷仓、煤堆等大型堆积态系统。现阶段，大型系统自燃规律预测的相关研究工作仍停留在宏观层面，多对系统内部结构采用均质假设。自燃系统往往是典型的非均质多孔结构，其内部的空隙变化会直接影响气流组织分布，进而改变了系统内部各区域的生热和蓄热能力。因此，探究堆积态系统内部的空隙率分布规律是自燃现象准确预测的重要前提条件。

当前，工程尺度下的堆体内部非均质空隙率分布测试仍是一个难题。虽然CT 扫描等可以在不损坏堆体的前提下获得三维堆体的真实结构[6-8]，但其样本范围多限于厘米量级。而对于大型的堆体，颗粒的松散滑动使得空隙率数据较难完好无误地提取。堆积态系统是通过颗粒从高空倾倒、下落形成，若能通过模拟仿真方法来重现物料颗粒的真实下落堆积过程，就可获得堆体内部结构，实现堆体几何结构的数值重构[9]，进而为堆体内部空隙率分布特性分析提供数据支撑。离散元法恰好可以模拟大量颗粒碰撞、迁移，追踪颗粒运动，确定下落颗粒实时位置，实现堆体的数值重构。

该文基于离散元方法模拟真实煤粒下落堆积过程，研究不同当量粒径和下落高度工况下的堆体非均质空隙分布，并将堆体空隙分布变化规律和与空隙率相关的变量编成自定义函数进行非均质煤自热模拟，使得堆体自燃研究从“黑箱”均质假设扩展到“灰箱”非均质。

1 堆积计算模型

1.1 物理模型

工程中的物料堆积多通过传送带送至特定高度，从落料口自由下落并堆积。本研究为简化过程，采用颗粒从漏斗落下过程来近似物料送出落料口的过程。在倾倒初期，总质量为10 kg 颗粒从漏斗颈部小口流出，进行自由落体运动，直至接触到底部圆盘，发生碰撞、反弹，进而缓慢堆积，最后积叠成一个小堆体。为保证模拟结果更贴近现实，通过查询Rocky 帮助手册确定煤粒的形状参数，其垂直纵横比为1.25，水平纵横比为0.7，平滑度为1，转角数为11，超二次曲面角为6.2。煤体与平台的物性参数[10]中漏斗与地面所用的钢材剪切模量为8.00×1010Pa，煤的剪切模拟量为4.70×108Pa，钢与煤的泊松比均为0.3，煤与煤的恢复系数为0.64，煤与钢的恢复系数为0.65。

1.2 数学模型

1.2.1 假设

为使问题简化以便数学建模，作如下假设：

1）忽略颗粒之间的黏性作用；

2）忽略颗粒的破碎作用；

3）煤粒非刚体，接触过程中允许出现微小变形，以弹簧模型来近似反映颗粒的挤压。

1.2.2 运动控制方程

式中：mi为颗粒i的质量，kg；Vi为颗粒i的速度，m/s；t为时间，s；ki为与颗粒i接触的颗粒数；Fij为颗粒i与j的接触力，N；g为重力加速度，m/s2。Ii为颗粒i的转动惯量，kg·m2；wi为颗粒i的角速度，rad/s；Tij为颗粒j作用于颗粒i的力矩，N·m。

1.2.3 碰撞模型

颗粒相互碰撞，其间的接触力可分解为法向与切向接触力。

式中：Fn为颗粒i与j的接触力的法向分量，即法向接触力；Fτ为颗粒i与j的接触力的切向分量，即切向接触力。颗粒间的法向接触力可用滞回线性弹簧模型描述，而颗粒间的切向接触力可用线性弹簧库仑极限模型描述。

1.3 模型验证

为了验证模型的准确性，本文模拟分析了实验台尺度下的煤粒下落堆积过程，并与实验结果[10]进行对比分析，发现煤体颗粒积叠形态的模拟和实验堆积角分别为28.27°和28.51°，仅为0.84%的相对误差，模拟与实验结果的堆积轮廓线基本吻合，呈现出较好的一致性。

2 空隙分布

收集分析1.5 cm 当量粒径颗粒从0.55、0.65、0.75、0.85、0.95 m 的高度下落以及1.1、1.3、1.5、1.7、1.9 cm 当量粒径颗粒从0.65 m 高度下落所形成的堆体空隙率数据。结果发现，堆体的垂直方向上的空隙率随垂直高度的增大而先减小后增大，堆体水平方向上的空隙率基本保持不变，而下落高度与颗粒粒径对堆体的空隙分布影响较小。

3 煤自热模拟

煤堆水平方向空隙基本恒定，随垂直高度的增大而先减小后增大，则煤堆整个区域可用垂直方向先减小后增大的变化趋势规律描述，应用此规律研究非均质煤堆自热过程。

3.1 几何模型及简化假设

考虑实际煤场中煤体堆积要防止滑坡以及再运输方便，多将其堆积为横截面为等腰梯形的大坝体。如图1（a）所示，该模型为二维，由两个域组成：

图1 气流作用下的煤堆（m）

1）多孔域煤堆为等腰梯形，底边长30 m，高10 m，底角为45°；

2）剩余区域代表煤堆周围开放空间的流场，长为90 m，宽为30 m。

将该模型进行网格划分，如图1（b）所示，采用全局尺寸为1 m 的四边形网格填充外部流场，局部尺寸为0.5 m 的三角形网格填充多孔域。煤堆与流场交界区域质量、热量交换最为频繁，是自热首先发生的地方，故在此处建立5 层厚度均为0.1 m的附面层薄层。

煤的自燃过程涉及多种反应，因此很难建立精确的数值模型，故必须进行一些简化[11]：

1）煤堆物性参数不随自热过程发生改变；

2）采用一步全局氧化反应来研究煤自燃过程；

3）忽略煤堆中影响较小热辐射效应和水分的作用；

4）流场与煤堆是一个整体，质量、热量交换均在其内部。

3.2 控制方程

煤一步全局氧化反应可简单表示为[12]

其 中：vO2、vCO2、vCO和（2vO2-2vCO2-vCO）代 表O2、CO2、CO 和煤氧复杂固体产物的化学计量系数。在30～170 ℃范围内，根据不同的反应机理将煤氧化反应分为三个阶段：30～60 ℃为第一阶段，60～100 ℃为第二阶段，100～170 ℃为第三阶段。具体参数设置见表1。

表1 数值模拟中煤氧化的设置

3.2.1 能量方程

煤在低温下的氧化速率极为缓慢，将其产生的热量定为源项，写成用户定义函数（UDF）。固体煤颗粒和气流之间的温差对煤堆内部的热相互作用影响颇大，需要将各自储存的能量以及二者的热能交换量分别表示出来。

1）固体能量守恒

式中：n为煤堆内部的空隙率；ρc为煤粒的密度，kg/m3；Cpc为煤粒的定压比热容，J/（kg·K）；Tc为煤体的温度，K；t为时间，s；λc为煤的导热系数，W/（m·K）；d为煤粒的直径，m；h为煤堆的对流传热系数，W/（m2·K）；Tg为煤堆内部气体的温度，K；r为耗氧速率，kmol/（m3·s），ΔH为消耗单位摩尔氧气产生的氧化热，kJ/mol。方程左边表示煤颗粒的内能增长，右边第一项为固体煤中的热扩散，第二项为内部气体相互作用的热对流，第三项为煤氧化产生的热量源项。

2）气体能量守恒

式中：ρg为煤堆内部气体的密度，kg/m3；Cpg为内部气体的定压比热容，J/（kg·K）；v为气流速度矢量；P为压力，Pa；λg为内部气体的导热系数，W/（m·K）。方程左边第一项表示气流的瞬态能量上升，第二项表示气流的热对流；方程右边第一项表示气流中的热扩散，第二项表示与固体煤相互作用的热对流。

3.2.2 组分守恒方程

在整个自热过程中，氮气既不消耗也不产生，因此主要研究反应物氧气以及生成物二氧化碳、一氧化碳的组分守恒。

氧气：

式中：Co、Ci、Cj分别表示氧气、二氧化碳、一氧化碳的气体浓度，kmol/m3；Do、Di、Dj分别表示气体混合物中氧气、二氧化碳、一氧化碳的质量扩散系数，kg·m2/s；DT,o、DT,i、DT,j分别表示气体混合物氧气、二氧化碳、一氧化碳的热扩散系数，kg·m2/s。

3.2.3 多孔介质中的动量平衡

加入粘性项对N-S 方程进行修正，求解多孔介质煤堆的动量平衡，得到：

其中多孔介质中两个较为重要的力为粘性力与惯性力，其表达式：

式中：τ为黏滞应力张量，N/m2；Fb为体积力矢量，N/m3；μ为气体粘性，kg·m/s2；k为煤基质的渗透率，m2；C1为粘性阻力，m-2；C2为惯性阻力，m-1。

综上可知，空隙率是众多自热因素的因变量，空隙率变化将会使得这些参数相应地发生改变。因此，必须将空隙率变化规律以及与空隙率相关的多孔介质煤堆内的粘性阻力、惯性阻力、氧气的消耗速率、二氧化碳和一氧化碳的生成速率等参量编译成用户自定义函数并进行代码调试，直至确定无误后，方可进行非均质煤堆的自热模拟计算。

3.3 煤自热模拟结果

将参数和求解的边界条件设置完毕，采用耦合算法、二阶迎风格式进行瞬态计算。在保持计算稳定性与收敛性的前提下尽可能地减少计算时间，将计算时间步长定为1 h，监测煤堆在热失控前的温度场、压力场以及组分场的变化。

如图2 所示，煤体10 d 升高3.5℃， 30 d 升高11℃， 60 d 升高26.5℃，90 d 升高100℃，随时间的推移，升温的速度极速加快，逐渐接近热失控。温度场的整体温度随时间的推移逐渐升高，温度较为明显地分成三个区域：煤堆的边界处由于与氧气接触最为频繁、反应产生热量最多而温度最高；再者是煤堆的中部区域；核心半圆处由于氧气最少、反应产热不足而温度最低。随着时长的增加，高温区域逐渐从与风流接触的外边界缩小至迎风侧、背风侧，再慢慢收缩成靠近两侧与底部的两个高温点而成为煤堆安全的最大威胁。

如图3 所示，煤堆内部的压力场主要可分为两个区域，靠近边界的压力加大，内部的压力较小。随着时间的推移，低压区域逐渐减小，由近似梯形变为半椭圆形，半椭圆的高度也是逐渐减小，到最后几乎贴着底部边界。

图3 煤堆流场随时间的压力变化图

关于煤堆内的各类别气体浓度变化，其氧气场与温度场的变化相近，中心有低氧区域较低温区偏大，二氧化碳与一氧化碳的浓度变化则与氧气浓度场变化相反。这是因为氧气要想通过重重煤颗粒的阻挡进入到中心区域是颇为困难的，进入到中心区域的氧气本就很少并快速与周边的煤粒反应生成二氧化碳、一氧化碳以及煤氧复合固体；相对而言，热量传导就容易得多，中心外部发生反应会将产生的热量往中心传导，致使低氧区中的低温区面积减小。

4 结论

为了消除现有堆体自燃预测中均质假设对结果准确性的影响，该文基于离散元方法模拟研究了大量颗粒从高空下落、堆积形成堆体的过程，获取了空隙分布的基本信息，并借此研究非均质煤的自热过程。分析概括系统内非均质性分布及其煤自热情况，得如下结论：

1）堆体内的空隙随垂直高度增加而先减小后增大，在水平方向上近似均匀。

2）研究范围内的倾倒高度大小变化以及粒径粗细变化对堆体内部空隙率影响较小。

3）非均质煤自燃中高温区域逐渐从与风流接触的整个外边界缩小至煤堆的迎风侧、背风侧，再慢慢收缩成两个靠近迎风侧、背风侧与底部接触的较其他区域温度高出很多的点，这两个“高温点”将成为煤堆安全的最大威胁。