干煤矿通风安全或者做采空区数值模拟的人大概率都遇到过这种纠结现场测点数据摆在那但是采空区内部的情况看不到只能靠经验猜。我也一样最开始用二维剖面算氧气和瓦斯的浓度分布算出来的结果总觉得差了点什么——尤其是上下隅角附近和高顶区域二维假设把垂向流动给抹平了很多关键现象根本显示不出来。后来我改用COMSOL搭三维采空区模型在通风条件下同时求解氧气浓度和瓦斯浓度分布整个物理图像一下子清楚了很多。这篇文章就从建模思路、方程设置、边界条件、求解技巧和结果判读这几个方面把这个三维模拟的完整过程写出来给正在做类似课题的人一个可以直接参考的路线。1. 先理清物理图景采空区里氧气和瓦斯到底怎么走的1.1 采空区本质上是高温差、强漏风的非均匀多孔介质采空区不是空的它是煤层开采后顶板垮落堆积形成的碎胀体由大小不一的岩块、碎煤、裂隙网络组成。从流动角度看它是一个典型的多孔介质区域只不过这个多孔介质极其不均匀靠近工作面的地方垮落岩块松散、孔隙大往里推进一段距离后上方岩层逐渐压实孔隙度和渗透率都明显下降。这种不均匀性直接决定了漏风风流在采空区内部的路径进而控制氧气和瓦斯的输运。做数值模拟第一步必须把这个物理图景立住。我自己的习惯是先画一张采空区的分区示意松散区、压实区、充分压实区分别对应不同的渗透率区间。这和我们常说的自燃三带散热带、自燃带、窒息带是直接相关的。1.2 两种气体方向相反却会汇合在同一条漏风通道里氧气和瓦斯在采空区内的行为可以简单概括为来源不同、命运纠缠氧气来自漏风。工作面进风巷的风流通过破碎煤体缝隙渗入采空区沿阻力最小路径向回风侧流动。沿途氧气不断被煤体低温氧化消耗浓度从21%逐渐下降。瓦斯来自采空区深部和垮落煤岩。残留煤体不断解吸释放甲烷在压力梯度和浓度梯度驱动下向漏风通道扩散、对流最终汇入漏风流被带回工作面回风隅角。一个往里走、一个往外走它们在采空区中间某个区域交汇。这个交汇区恰恰是最危险的甲烷浓度如果落入5%~16%的爆炸界限同时氧气浓度又高于维持燃烧所需的临界值一般认为低于8%就基本不具备氧化自燃条件甲烷爆炸则需要氧气浓度大于12%左右那就同时具备了爆炸和自燃的潜在条件。三维模型要解决的问题就是把这个交汇区的位置和范围在空间里精确地画出来。2. 三维几何建模尺寸、分区和参数一次给定2.1 走向、倾向、高度三个方向的尺寸怎么取COMSOL建模前尺寸不能拍脑袋。以我做过的一个走向长度200m、工作面倾向长度150m、采高3m的典型长壁工作面为例采空区三维几何大致是这样取的走向X工作面推进方向取100m到200m。采空区深部压实后浓度变化梯度非常小取太长意义不大反而增加网格量。倾向Y平行工作面方向直接取工作面长度150m两端分别延伸到上下顺槽。高度Z垂向这是二维模型最容易忽略的。垮落带高度一般按采高的2.5~5倍估算所以采高3m时垮落带上限可以取8~10m。再加上上方裂隙带的一部分模型总高取15~30m更有代表性。一个需要注意的地方是很多文献建模把采空区想成规则的矩形盒但这种模型对垂向流场低估得很明显。我更建议把垮落带做成靠近工作面高、深部逐渐压实变低的楔形或者至少不同高度分层给不同渗透率。这样上下隅角附近的流场形态才贴近实际。COMSOL里几何搭建本身不复杂用块Block或者拉伸Extrude就能完成。关键是要把下面的分区参数留好后面赋材料时才有区分度。2.2 材料参数分区渗透率和孔隙率的取值参考采空区多孔介质的两个核心参数是渗透率k单位m²和孔隙率εp无量纲。我给的参考值如下表具体数值要根据岩性、采高、开采深度和经验公式修正分区位置特征渗透率参考范围m²孔隙率参考范围松散垮落区散热带距工作面0~20m1×10⁻⁸ ~ 1×10⁻⁷0.25~0.40压实过渡区自燃带/氧化带距工作面20~60m1×10⁻¹⁰ ~ 1×10⁻⁹0.15~0.25深部压实区窒息带距工作面60m以远1×10⁻¹² ~ 1×10⁻¹⁰0.05~0.12COMSOL里不需要把这个表做成连续函数最简单的办法是建立三个域分别赋不同的材料参数。如果想更精细也可以用解析表达式写成分段函数甚至用插值函数从实测渗透率数据生成空间分布。建议至少细分三层以上否则漏风路径过于理想化。3. 风流场与浓度场耦合控制方程与边界条件3.1 用达西定律模拟漏风流场而不是N-S方程采空区漏风速度极低一般每秒只有几毫米到几厘米对流项惯性效应弱本质上属于蠕动流。这种情形用完整的Navier-Stokes方程纯属浪费计算资源COMSOL里的达西定律接口Darcys Law才是正确选择。达西定律控制方程为[ \mathbf{u} -\frac{k}{\mu} \nabla p ]其中u是达西速度矢量m/sk是渗透率m²μ是气体动力粘度Pa·sp是压力Pa。再结合连续性方程[ \frac{\partial \rho \epsilon_p}{\partial t} \nabla \cdot (\rho \mathbf{u}) Q_m ]Q_m是质量源项瓦斯涌出就体现在这里。稳态问题时第一项可以去掉。这里要重点理解一点达西速度u不是孔隙内部的实际流速而是宏观面通量。组分输运方程里用的对流速度应该用达西速度除以孔隙率后得到的真实流速即v u / εp。COMSOL的多孔介质稀物质传递接口已经内置了这个换算使用时要确认勾选了多孔介质属性里的孔隙率。3.2 氧气和瓦斯的输运方程对流-扩散-源项浓度场我用COMSOL的多孔介质中的稀物质传递Transport of Diluted Species in Porous Media接口。对组分i氧气O₂甲烷CH₄控制方程是[ \frac{\partial (\epsilon_p c_i)}{\partial t} \nabla \cdot (\mathbf{u} c_i) \nabla \cdot (\epsilon_p D_{e,i} \nabla c_i) S_i ]这里c_i是组分摩尔浓度mol/m³D_{e,i}是有效扩散系数m²/sS_i是源项。甲烷的源项S_ch4来自采空区残留煤体解吸释放单位是mol/(m³·s)。一个常用的估算方式是根据工作面的绝对瓦斯涌出量比如5~20 m³/min除以采空区体积折算一个体积平均涌出强度。更精细的做法是只在深部压实区赋源项因为浅部垮落区内散煤释放快深部煤体持续释放是长期瓦斯的补给来源。氧气的源项相对复杂。如果只做物理输运可以先不设氧消耗项得到的是保守结果——实际氧化耗氧会让氧浓度进一步降低。如果想同时间接模拟氧化带范围可以加一个简单的线性耗氧项[ S_{O2} - k_r c_{O2} ]k_r是耗氧速率系数1/s取值需要靠煤样自燃倾向性实验来标定不能瞎给。工程上先不加反应项、把输运结果作为上限来看往往是更稳的做法。3.3 边界条件怎么设进口、出口、围岩一个都不能漏边界条件是这类模拟最影响结果的部分。我的设置习惯是边界位置类型设置内容进风侧工作面进风隅角压力边界给定相对压力p_in氧浓度0.21 mol/mol甲烷浓度接近0回风侧工作面回风隅角压力边界给定相对压力p_out低于p_in一般压差设50~200Pa采空区顶部、底部、走向深部无通量边界法向流速为0组分法向通量为0瓦斯释放区域源项域甲烷体积源项按涌出强度折算这里有个关键进、回风边界不能都设成固定浓度。出口侧应该让组分按对流流出也就是设置成流出条件否则甲烷会在出口边界人为堆积浓度虚高。COMSOL稀物质传递接口里默认的出口条件对流通量就是这个逻辑不需要改成固定浓度。压力边界压差的设定决定了漏风量建议先用单相流场单独跑一遍观察采空区的漏风速度分布是否在合理量级工作面附近毫米到厘米每秒。如果速度过大或过小调渗透率比调压差更符合实际。4. 求解策略与常见坑稳态、网格和收敛4.1 稳态还是瞬态各自解决什么问题这个选择取决于你想回答什么稳态模拟假设采空区风流稳定、瓦斯涌出稳定用稳态求解器同时解流场和两个浓度场。优点是计算快内存占用小几十万网格几分钟到几十分钟就能算完。适用于日常的三带划分、危险区域静态识别。瞬态模拟需要关注瓦斯涌出波动、停风/恢复通风、采空区气体累积过程时使用。例如模拟掘进面停风后采空区瓦斯向外扩散积聚的过程时间步长要取得够密通常0.5~1h一个输出点跑上几天墙钟时间都正常。我的路线是先稳态再用稳态结果当瞬态初值这样比瞬态从零开始收敛快得多。COMSOL里直接用研究1的稳态结果作为研究2的瞬态初值这个操作在扩展研究时非常好用。4.2 网格剖分与求解器配置的实例网格策略直接影响收敛。采空区几何规则我一般这样剖走向和倾向方向用自由四面体最大单元边长控制在5~8m在进风侧和回风侧隅角附近加密因为浓度梯度最大垂向方向至少剖4~6层尤其垮落带上部要单独设定边界层否则瓦斯垂向积聚现象剖不出来总网格量控制在30万~80万之间再多就要考虑服务器内存了。求解器配置上稳态流场用固定阻尼或PARDISO直接求解器都很稳定。如果流场和浓度场双向耦合不强实际也确实不强可以分步先单独算达西流场冻结流场后再算两个组分输运。这种冻结速度场的策略能把很多莫名其妙的收敛问题直接绕开。COMSOL里用研究步骤分割或者先禁用组分接口跑流场再启用组分接口并继承流场结果就行。4.3 不收敛和浓度振荡的排查思路做这类多孔介质输运模拟最容易遇到两类问题第一类流场不收敛。原因排在前三的都是渗透率差异跨度过大松散区到压实区差了5~6个数量级、压差太大导致局部达西速度畸变、网格质量差。排查方法很简单把渗透率分布用等值面和切片显示出来再看流线是否出现不合理的锯齿。如果流线在某阶梯处骤变那大概率是渗透率分区赋错了域或者边界层网格没剖好。第二类浓度场负值或振荡。一般出现在网格太粗的区域特别是甲烷源项强而网格没能分辨浓度边界层时。COMSOL的PARDISO可以快速算但流的对流占优问题时可能需要打开流线扩散streamline diffusion。默认设置下我习惯把它设为0.1~0.2倍的特征对流尺度能明显压低振荡。这是个经验值具体还要看模型和网格。5. 三维结果怎么读自燃三带与瓦斯积聚区5.1 三维氧浓度切片怎么切才不遗漏算完第一件事是看氧气浓度的三维切片而不是只看一个平面。常用的切片包括水平切片取采空区底面、垮落带中部和顶部分别切看氧浓度在垂向上的差异垂直切片沿工作面走向切一条、沿倾向切几条用来分析上下隅角的局部流场三维等值面把氧浓度18%和8%两个临界值做成等值面直接生成散热带/氧化带/窒息带的界面。我印象最深的一次是同一个走向断面上底部氧浓8%的位置和顶部氧浓8%的位置差了将近15m。如果只用二维剖面去划自燃带边界这个空间差异就完全丢失了。三维等值面一出来上部氧化带明显更深入采空区深处这就是三维模拟的价值所在。5.2 甲烷浓度在垂向上的积聚规律瓦斯密度小于空气在有垮落带碎胀空间里天然有向上运动的倾向。三维模型算出来后甲烷浓度在垂向上的分层清晰可见顶板附近甲烷浓度往往比底板高出不少尤其是在采空区深部那里漏风微弱、对流不强甲烷只能靠扩散传播积聚效应更加明显。这直接关系到工作面回风隅角瓦斯的治理很多瓦斯超限不是来自工作面直接涌出而是采空区高浓度瓦斯从顶部裂隙带流向回风隅角。三维模型里用流线甲烷浓度彩色显示可以看到非常直观的瓦斯向回风隅角汇聚的过程。这时再去考虑Y型通风U型通风高位抽采等方案时位置选点偏上还是偏下、抽采钻孔打在哪个高度就对得上了。5.3 把爆炸危险区域看出来数值模拟最终要服务于安全判断。我的做法是在COMSOL里新增一个自定义表达式[ \text{Risk_CH4} \text{if}(c_{CH4}0.05 \text{ and } c_{CH4}0.16 \text{ and } c_{O2}0.12, 1, 0) ]把甲烷爆炸界限体积分数5%~16%和氧气临界值大于12%同时满足的区域标成1其余为0。渲染出来就是一组爆炸危险体元。这个三维危险区往往并不是一团连续的云而是断断续续的孤岛区分布在回风隅角高风险通道附近。真实工作面的抽采钻孔、监测束管布置应该重点覆盖这些位置。我还会把它和氧浓度18%等值面、8%等值面叠加在一起同时判断煤自燃氧化带和瓦斯爆炸区的空间关系。两套危险区域如果重叠那这个采空区的灾害耦合风险就很高必须优先预警。6. 模型参数的敏感性、现场校验与时间成本6.1 渗透率和瓦斯涌出强度是两个最大的变量在所有输入参数里渗透率k和瓦斯涌出源项S_ch4对结果的影响最显著。渗透率决定了漏风的通道和流量差一个数量级氧化带整体就移动好几个位置瓦斯涌出强度则直接决定爆炸界限区域的范围。我建议做模拟时不要只跑一组参数。至少算三组渗透率取参考值、放宽、收紧瓦斯涌出量取低值、均值、高值。然后把危险区域的范围列个表对比一下看看结论对哪个参数最敏感。通常你会发现在采空区深部无论怎么调参数氧气都被消耗到窒息带水平结论稳定而在中部氧化带结论可能完全翻转。这时候现场监测方案就要往最敏感区域加测点。6.2 用现场实测数据做粗校核模拟再漂亮也要和现场对上。采空区气体监测最常用的是束管监测把取样头埋进采空区定期抽气分析氧和甲烷浓度。COMSOL模型算完后可以和束管监测的浓度剖面做逐点对比。误差在10%~20%以内就算合理不必追求严丝合缝。一个有效的方法是沿着束管布置的方向在模型里设一条截线导出氧浓度和甲烷浓度的曲线再叠加上实测点。如果趋势一致、绝对数值偏差大优先检讨瓦斯涌出强度的标定如果趋势都不一致多半是渗透率分区方向搞错了或者模型边界条件设置与真实通风系统不符。这一步是模型可信度的试金石。6.3 一些降低试错成本的实际心得最后分享几个我用COMSOL做采空区三维模拟攒下来的细节经验都是常规软件教程里不太会写的先跑二维预判参数。在几乎所有参数没标定之前先用二维切面模型把达西流场跑通确认压差、渗透率、源项的量级合理再建三维。三维模型的调试成本高多花这一步能省大量时间。几何尺寸别贪大。采空区走向取太深深部完全窒息区对目标问题毫无贡献还要平白增加网格。按照你关心的危险区域一般就是工作面上方和上下隅角附近反推模型范围效率最高。把进出口压差当成主动变量来标定。先固定渗透率调整压差使模型中的漏风量和工作面实测的漏风率一致之后再动渗透率。参数之间的纠缠是这种多参数模型最难的部分每次只调一个变量的原则必须坚持。后处理和结果导出尽早自动化。COMSOL的结果导出可以保存为数据文件切片图、等值面图都先输出来再统一画图。否则每改一次参数都要重新处理一遍可视化繁琐而且容易出错。三维采空区通风模拟里氧气和瓦斯浓度分布这两个量一旦解耦分析很多安全问题都能在模型里提前暴露。我自己的体会是这类模型永远做不到100%精准但它的价值在于把现场看不见的危险区域从猜变成算给通风设计和灾害防治一个可以辩论的基准。这篇文章里写的参数、步骤和踩坑记录都是我在实际课题里反复试过的希望能让后来的人少走一段弯路。