★ FEATURED ARTICLE 矿井里最难发现的火不在明面上而在你看不见的角落。煤体深部、采空区垮落带、矸石堆的缝隙里温度正在缓慢爬升。没有火焰、没有浓烟等你通过束管气体监测或者钻孔测温察觉到异常时氧化带往往已经蔓延了相当大的范围。这种无焰的缓慢燃烧就是阴燃它最麻烦的地方在于氧气靠着漏风通道一点一点往里送反应热又被周围松散的冒落介质捂着散不出去温度就在这种“供氧不足但散热也差”的条件下慢慢积累直到越过临界点。这个阶段几乎没法靠现场实测拿到完整数据。采空区垮落带是人员不可能进入的区域测点布不进去缩尺实验又很难还原垮落体的非均质结构。所以行业里普遍的做法是数值模拟。我最近一直在用COMSOL做这一块的二维模拟思路很直接用达西定律接口算采空区气体速度场用传递类接口算氧气浓度和瓦斯浓度用固体传热接口算温度场再通过化学反应源项把四个场耦合起来。刚跑完的一个二维模型输出了速度分布、氧气浓度分布、瓦斯浓度分布和温度分布四套结果对判断阴燃氧化升温带的位置和自燃风险区域参考价值已经很高了。这篇文章把我从几何搭建到模型收敛的完整过程记录下来包括中间踩过的坑和最终的调参顺序。如果你也在做煤自燃或者阴燃相关的模拟工作照着这个思路能省掉不少试错时间。1. 为何只能用模型来盯阴燃实验手段够不到的地方1.1 阴燃和明火的本质差别我最早接触阴燃这个概念时第一反应是它和普通燃烧有什么区别。后来在模型里反复跑温度分布才真正理解阴燃是低温、无焰、以氧化反应为主的缓慢放热过程典型温度区间比明火燃烧低得多很多时候连烟都没有。它的自维持靠的是两个条件——持续但微弱的氧气补给以及周围介质把反应热捂住不让它散掉。采空区天然具备这两个条件。垮落体孔隙率大漏风可以渗透进去但对流不强煤岩体的导热系数又不高氧化产生的热量很难快速传导出去。于是阴燃可以像一个慢炖锅一样几个小时甚至几天里温度从常温慢慢爬到几百摄氏度。一旦温度积累到某个阈值反应速率指数式上升就可能转向明火。问题是等你看到明火前面漫长的阴燃阶段已经错过了扑救窗口。这个特点决定了阴燃阶段的信息在时间和空间上都是高度隐匿的你想知道它在哪里发生又想知道它发展到什么温度了但现场根本没有手段直接看到内部。所以模型在这里不是锦上添花而是唯一能提供“全空间场”信息的手段。1.2 现场监测为什么喂不饱一个模型有人可能会问现场不是有束管监测和钻孔测温吗确实有但采空区的测点在空间上是极其稀疏的。钻孔只能打在有限的几个位置束管抽气分析的是某个取样断面的混合气体这些数据描述的是“点”和“线”不是分布场。更要命的是气体在地下运移的路径很复杂。采空区内部的漏风风流会沿着大孔隙通道走捷径束管抽到的样气可能来自十几米甚至几十米外的区域测到的氧浓度并不代表你关心的那一点的真实状态。温度测点同样受位置影响极大垮落体的非均质性导致温度在很短的距离内就能出现明显差别一个钻孔测到低温完全不代表旁边没在升温。所以现场数据可以用来标定模型的边界条件和验证整体趋势但没法直接回答“氧化带在哪、温升速率多少、氧和瓦斯的空间分布长什么样”这些问题。要回答这些只能靠求解控制方程把整个区域的速度场、浓度场和温度场都算出来。1.3 二维模型到底能回答哪些工程问题这是我入坑阴燃模拟后经常被问到的。二维模型确实有简化但很多工程问题在这个尺度上是能回答清楚的氧化升温带在采空区深度方向上的大致位置以及它随漏风强度的迁移规律氧气浓度从边界向深部的衰减曲线判断氧化的极限深度瓦斯释放对氧浓度稀释作用的大小评估瓦斯积聚与阴燃之间的相互影响温度场的峰值位置和量级为防灭火钻孔设计提供平面参考当然二维模型算不出三维空间里局部冒落不均形成的涡流和绕流也在捕捉复杂漏风路径上有先天局限。但对做参数敏感性分析和趋势研判来说二维模型的性价比高得惊人。这一点后面我会专门展开。2. 二维几何搭建与参数标定把垮落体压缩成一张能算的图2.1 几何域的切法走向剖面而不是俯视平面我刚做第一个采空区二维模型时想当然地画了一个俯视平面图工作面在左边采空区往右延伸进风巷和回风巷分列上下。跑出来的结果总觉得怪怪的——温度分布对高度方向完全没有分辨能力而实际上垮落体在竖直方向上的空隙率变化恰恰是决定冒落区氧化自燃的关键因素。后来换成了走向剖面x方向从工作面沿走向深入采空区典型长度100到200米y方向从底板往上覆盖冒落带高度取30到40米。这样的几何虽然简单但能自然地把“沿深度方向的漏风衰减”和“沿高度方向的空隙率变化”同时装进模型里。左侧进风侧右侧回风侧顶底板都保留一层岩体参与热传导避免边界上直接把热量流失掉导致结果失真。二维剖面模型有个额外的好处它本质上就是一个等价的“研究单元”可以把三维采空区的空间变异问题转化为沿两个方向的变化规律物理图像清晰调试起来也直观。我自己建议如果你是从零开始先用这种剖面模型跑通流程再考虑更复杂的三维局部模型。2.2 空隙率与渗透率的分区赋值采空区垮落体的空隙率绝不是常数。靠近工作面的区域垮落充分、碎胀系数大空隙率可以到0.3以上随着深度增加上覆岩层压力逐渐压实矸石空隙率会掉到0.1甚至更低。我最初偷懒全模型给了一个均匀空隙率0.2结果氧气在模型里输运得异常顺畅氧化升温带直接不存在。因为氧被均匀地带到了深部反应热也被均匀稀释了温度根本集中不起来。后来我按照采空区“散热带—氧化升温带—窒息带”的经典三带划分把几何域切成了三块每块独立赋值。空隙率的基准分布我用的是指数衰减函数配合分区修正从工作面侧约0.30沿深度衰减到深部约0.08再在高度方向按冒落压实程度调整。渗透率则是通过空隙率换算用多孔介质常用的Kozeny-Carman类关系把空隙率和渗透率关联起来颗粒特征尺寸取垮落矸石的平均块度。这里有个实操经验渗透率的量级对速度场的影响比空隙率大得多而空隙率又直接参与浓度方程和传热方程里的体积平均计算。所以赋值时两者要一起调不能只调一个。二维模型的好处是每一次试算都很快可以多跑几组参数感受变化再把组合固定下来。2.3 边界条件的给法速度入口、压力出口和顶底板的坑边界条件直接决定模型的“说服力”。我一开始用固定压力边界模拟全采空区漏风发现速度场分布太平滑深部几乎没什么风后来换成在工作面侧给定速度入口倒是接近现场实际情况了但又出现入口附近的局部流速过大氧化带被“吹散”了。最后还是回到实用组合入口给达西速度边界数值取漏风实测或经验估算出口用开放压力边界上边界和底边界处理成不可渗透但热边界绝不能简单设为绝热——不然温度场会出现虚假高热区。我现在的做法是让几何域外侧多画一层岩体岩体边界再给对流换热条件这样温度场的边界失真会小很多热量传导路径也更贴近真实。另外提醒一个容易被忽略的细节如果二维剖面模型中进风侧和回风侧的高度不一致或者模型本身不对称速度场会出现明显的偏流这种偏流对氧气和瓦斯分布的影响很大。画几何之前先想清楚漏风路径的主方向别把模型画得太理想否则后面的浓度场和温度场全会被带偏。3. 四个场的咬合关系速度、氧气、瓦斯、温度怎么互相牵制3.1 速度场是所有输运的基础采空区里的气流速度很低深部漏风速度常在小到每秒毫米甚至更低气流在这种多孔介质的细小孔隙里流动惯性力占比很小黏性力主导。所以达西定律作为动量方程在这个场景下是完全合理的COMSOL的达西定律接口直接能给出速度和压力分布。我在模型里把入口达西速度设在一定范围内数值中心化在每秒几十毫米左右模拟结果是漏风在上部松散区速度稍高越往深处衰减越快这跟工程经验是一致的。速度场是四个场中的“骨架”。氧气的对流输运要靠它瓦斯气体从采空区深部向回风侧的运移也要靠它温度场的对流散热量更是直接取决于速度大小。我在调试时第一条原则就是先把速度场单独跑稳确认风向、风速量级合理再往下走。如果速度场不对后面所有浓度场和温度场都是空中楼阁。3.2 氧气浓度场对流和扩散谁说了算氧气在采空区的输运方程由对流、扩散和反应消耗三项组成。采空区近边界处对流占主导氧气跟着漏风往里走到了深部流速已经非常小分子扩散和机械弥散开始起主要作用。判断哪个机制主导可以用Péclet数估算——对流输运和扩散输运的比值。我在模型中用低流速和中等弥散系数时算出来的Péclet数在深部区域下降到接近1甚至更低说明深部的氧供给确实依赖于扩散和低速对流这也从机理上解释了为什么阴燃的氧化带总是出现在某个深度区间而不是无限延伸。浓度场和速度场的耦合不是单向的。如果使用物理上更严谨的浓物质传递而不是稀物质近似浓度变化会影响混合物的密度和黏度反过来又会影响压力场和速度场。不过在实际操作中氧气体积分数从边界处的约21%降到深部的接近0%变化范围很大严格说稀物质近似存在物理失真。我个人的做法是先用稀物质模型快速进行趋势判断在关键工况复核时再切换到浓物质或者手动修正对流项兼顾效率和可靠性。还有一个和氧气模型强相关的量温度场的边界漏风蓄热。漏风温度和采空区环境温度不一致时风流本身就会带来热量的传递。如果漏风温度比煤岩体温度高氧化热像叠加了一层保温如果相反它又变成散热项。我在模型中把入口温度单独设置并参与热方程耦合跑出来跟不考虑漏风温度的结果差别很明显。这一项很多人做模拟时会忽略。3.3 瓦斯浓度场它不是旁观者而是氧气稀释剂采空区瓦斯来源主要是遗煤解析和围岩涌出。在二维模型中我处理的瓦斯浓度场其实很有讲究瓦斯本身是不可燃的气氛组分在阴燃阶段它更大的作用是把氧气稀释掉压缩氧化反应的空间可它又是潜在的可燃气体如果温度足够高且氧气恢复了供应就可能引发更为凶险的瓦斯参与燃烧。我在模型里设置瓦斯释放源项时用了拾取区域的方式把氧化升温带附近的煤体标记成源区瓦斯释放强度按照煤的瓦斯含量和解析特性的经验数据给一个分布式源项。跑出来的结果显示在瓦斯释放较强的区域氧气浓度下降的速度明显加快原本可能发生氧化的区域被“挤占”掉了。这个结果在工程上是有意义的瓦斯积聚区和氧化升温带在空间上的错开与否直接关系自燃与瓦斯耦合灾害的风险。有意思的是有些模型会把瓦斯当成单纯的可燃组分参与燃烧反应这在我看来是值得商榷的。因为采空区在阴燃阶段的供氧条件极差瓦斯很难在缺氧情况下持续燃烧而它对氧气的稀释作用反而是更普遍的过程。所以瓦斯浓度场和温度场的耦合在阴燃模拟里应该以“稀释氧”为主导而不是“燃烧放热”。3.4 温度场放热和蓄热的一场拔河温度场方程的核心是热平衡氧化放热、对流散热、导热扩散再加上边界换热。氧化放热速率随温度按指数增长这是整个模型的非线性发动机。我做模型时最关注的就是这个正反馈回路温度升高加速反应、反应放热继续推高温度如果不加以控制或没有足够散热解就会一路冲向“爆炸”数值上也表现为发散。蓄热条件在阴燃模拟里是被低估的一个概念。多孔介质的热物性参数按空隙率做了体积平均等效导热系数通常会低于完整煤岩体。空隙率越大介质越“蓬松”导热性反而越差这正是阴燃能闷在里头持续升温和蓄热的原因。我测试过把等效导热系数调为常数和调为随空隙率变化两种方式后者得到的温度分布会形成更明显的局部高温核也更符合阴燃“点状升温”的特征。温度场对速度场的反作用同样不能忽略。温度升高会改变气体的黏度和密度在采空区这种低速区域浮力驱动的自然对流有时候能改变局部流场方向。这个效应在二维走向剖面上表现得比较明显高氧化带温度升高后热气会向上浮升抽吸周围的空气形成局部的循环气流反过来增强氧化带的氧供应。我第一次看到这个现象时还挺意外的因为它说明阴燃在空间上是会“自己去找氧气”的。4. 化学反应源项整个模型最敏感、最容易翻车的一环4.1 低温氧化速率别追基元反应用好表观动力学煤炭低温氧化过程是一长串自由基链式反应真正把所有基元步骤都耦合进CFD模型是不现实的行业里通用的做法是使用总体表观动力学反应速率。我用的形式是经典的Arrhenius幂律R A · C_O2^n · exp(−Ea / RT)其中A是指前因子C_O2是氧气浓度n是反应级数Ea是表观活化能R是气体常数T是绝对温度。煤炭低温氧化文献里Ea常见取40到70 kJ/mol量级n在0.5到1之间。这里的难点在于不同煤种的参数差异很大如果你直接把文献里的A和Ea抄到自己的模型里大概率会出现两个极端温度完全起不来或者几分钟内爆表。我做过一组计算验证这个敏感性Ea只差5 kJ/mol同一个入口风速下温度峰值能差出上百摄氏度。原因是Arrhenius项对温度本来就是指数级响应参数稍有偏差正反馈的强度就完全不在一个量级上。我的标定思路是先固定Ea在常见范围调整A去匹配煤样绝热氧化实验的起始升温速率然后再做一组敏感性扫描观察温度场的变化范围。这种“先定趋势再精调量级”的思路比较稳健。4.2 放热和耗氧必须来自同一套速率模型里最容易出现的错误是把氧气消耗和温度放热当成两回事独立设置。实际上氧气消耗速率和反应放热速率必须同时由同一个反应表达式给出它们之间只差一个反应焓系数。我在COMSOL里用一个全局变量或表达式统一控制两个场方程的源项这样就能从机制上杜绝“一面猛放热、一面氧浓度纹丝不动”这种物理上荒谬的结果。具体到量级上单位质量的煤氧化反应焓放热大约在数百千焦每摩尔量级对应到耗氧和放热的配比我会在数学关系里同时乘上一个放热比例常数。这里要特别留神单位换算因为浓度方程用的是mol/m³或体积分数而传热方程用的是W/m³两者差一个反应焓的换算算错一个量级整个温度场都会面目全非。这也是为什么很多阴燃模型看起来结构完整结果却经不起推敲的原因。我在实际操作中还会把放热源项做一个平滑处理避免它在数值上过于剧烈。反应速率在低温区很小但一旦温度越过拐点会急剧上升如果不加平滑很容易导致单步内源项从极小跳到巨大求解器自然就崩了。平滑函数不会改变物理趋势但能显著提高数值稳定性。4.3 反应源项与网格尺度的恩怨反应源项是高度的空间集中量它只在氧化升温带附近显著不为零。如果网格太粗这个源项会被稀散到整个单元里局部温升根本成形不了如果网格太细局部反应强度跨越网格时波动太大时间步长会被逼得很小计算效率急剧下降。我的实践做法是分两步走先在整个域用粗网格把流动场、浓度场跑通确认整体趋势然后把氧化带可能经过的区域加密网格做一次网格无关性检查。加密后如果温度峰值还在变说明网格仍然不够如果温度和反应速率随网格加密趋于稳定说明当前网格尺度是可靠的。二维模型做这种检查成本很低多跑几组网格没问题。这个步骤别省省了之后换参考数据时出了问题你根本说不清是参数的问题还是网格的问题。5. 收敛调试实录从发散到稳定运行的完整调参顺序5.1 我踩过的三次翻车第一次翻车是初始条件不对。我一开始就把全场的反应源项完整打开温度迅速从环境温度直接飙升到虚高的高温紧接着氧气浓度在局部跌成负值模型瞬间崩掉。问题出在初始条件没有渐进过程冷态无风场的初始状态直接面对强放热反应数值上等同于一个阶跃冲击。第二次翻车是时间步长失控。阴燃的物理过程跨度很大背景流场是秒级到分钟级反应动力学特征是分钟级但整体温度演化要模拟几小时到几十小时。COMSOL自动时间步长有时候会冲过头跨过了反应速率的快速变化区间导致浓度出现振荡甚至负值。后来我把最大时间步长卡死在物理过程允许的范围才基本消除了这个现象。第三次翻车是入口氧浓度的“无限供氧”。我把入口边界设成固定氧体积分数21%但由于计算域太长这个边界条件导致整个模型区域的氧气被持续“喂饱”深部的氧化带始终等不到耗尽。后来我调整了模型的计算域长度和入口到氧化带的距离并且允许出口侧回流的空气组分参与混合问题才得到改善。5.2 我的调参顺序先流场、再输运、最后开反应经过这些翻车之后我总结了一套固定的激活顺序现在跑新工况都按这个来第一步只求解速度场稳态收敛。确认漏风路径、速度量级和压力分布合理后再往下走。这一步通常几十秒就完成也是后续一切的基础。第二步打开氧气和瓦斯的输运方程但先不加反应源项。看看氧气在纯对流和扩散作用下能输运到多深瓦斯的积累区域在哪里。这一步能提前暴露边界条件和弥散参数的不合理。第三步把反应速率从最终值的10%开始逐步提高。我通常分四档10%、30%、60%、100%。每一档都跑到温度场稳定后再提高相当于给非线性正反馈一个逐步加码的适应过程。第四步在反应全开的基础上调整热边界条件和等效导热系数观察温度峰值是否合理氧化带是否和实际工况对得上。这套流程看起来啰嗦但它有一个核心价值每一段不稳定都能立刻定位原因。如果不按这个顺序一次性全部打开模型崩了你根本不知道是流场问题、扩散问题还是反应问题只能靠猜那就彻底陷入黑暗调试了。5.3 用二维模型做敏感性分析是真香模型收敛之后二维模型最大的价值就体现出来了跑得快、改参数方便、可以系统地做扫描。我最近用固定模型做了一组工况扫描分别调整入口速度、空隙率衰减系数和反应活化能。一套十组工况跑下来一个晚上就能出结果这在三维模型里是根本不敢想的量级。这组扫描带来的一个典型认识是漏风速度增大时氧化升温带会向采空区深部迁移同时温度峰值会下降因为对流散热和氧输运同步增强。这个趋势如果单独看一组结果是体会不出来的只有批量扫描才看得出来。工程上这个结论直接关系着注氮、注浆和风量调节措施的参数选择。当然二维模型不能替代三维精细化研究。我的个人习惯是先用二维模型做敏感性定位确定关键参数和重点区域然后针对性建立三维局部模型细化验证。这个工作流既保证了效率又能在关键点上给出可信的空间细节。最后分享一个和长时间模拟有关的小经验阴燃模拟的时间跨度大如果每个输出时间步都保存所有场变量结果文件会膨胀得非常快。我在输出设置里固定每600秒保存一帧既保留了温度演化的完整过程又不会把磁盘撑爆。另外模型的实用性很大程度上取决于边界条件是否反映真实工况建议一开始就选定一组已知的现场升温数据作为基准标定工况先把这组数据复现了再去拓展参数研究不然所有后续分析都没法站稳脚跟。
阅读完成 · 觉得有帮助？
