采空区气体流场多尺度数值模拟及热蒸发效应分析

杜 斌，赵 磊

（1.中煤科工集团沈阳研究院有限公司，辽宁 抚顺 113122；2.煤矿安全技术国家重点实验室，辽宁 抚顺 113122；3.平顶山天安煤业股份有限公司，河南 平顶山 467099）

采空区遗煤自然发火，威胁井下安全生产，还可能引发采空区瓦斯爆炸，其条件之一是采空区漏风持续地为遗煤供氧。因此国内外许多学者通过现场实测、数值模拟等方法对采空区内气体流动的速度场、压力场、氧浓度场开展了研究。

采空区内气体流动规律在宏观尺度服从N-S方程，气体流动驱动力主要是工作面上下隅角的气压差。然而，煤与氧的物理吸附、化学吸附、化学反应过程主要在煤基质的表面（气-固界面）进行，包括煤块与采空区气体的宏观气-固界面和煤块中孔隙气体与基质的介观气-固界面。因此除了需要对采空区宏观气体流动规律进行研究，还需要对煤孔隙中的介观气体流动规律进行研究。然而由于研究尺度不同，煤孔隙中的气体流动将不再像宏观气体流动那样服从N-S 方程，而是服从达西定律。梁运涛等[1-2]将发生自然发火的遗煤按照从小到大的层次划分为5 个尺度：反应粒子、颗粒内孔、单煤颗粒、颗粒聚团、堆积媒体，提出使用关联型多尺度方法，对自然发火过程进行跨尺度研究。但是目前尚未有学者对2 种尺度下多孔介质内的气体流动进行关联型分析。另外，煤孔隙中的残余水分占据了气体流动通道，甚至存在“水锁”效应，阻碍了气体流动，同时，气体流动产生的蒸发效应，加快了残余水分的蒸发，减小了气体流动通道内阻碍物质的体积，增加了有效渗透率，利于气体流动，因此需要考虑残余水分与气体流动间的耦合影响，即气液两相流之间的耦合影响；同时水分蒸发吸热使煤体温度降低，又反作用于气体流动和水分的蒸发速率，使气体分子的热运动和蒸发速率降低。因此，对介观尺度下煤中热-流-湿多场耦合机制进行研究。

1 采空区多尺度气体流动规律

在矿井采场空间范围内共存在5 个层次的研究尺度，其中边界尺度为反应粒子、单煤颗粒、堆积煤场，介于边界尺度之间的为颗粒内孔（介尺度I）和颗粒聚团（介尺度II）。从宏观尺度入手，按照尺度逐渐减小的顺序，对采空区松散煤体及之间的空隙、遗煤基质及孔隙裂隙内的气体流动情况展开研究。

1.1 采空区宏观尺度气体流动规律

选取长度300 m、宽度180 m 的采空区堆积煤体范围作为研究对象，其中0～40 m 为自然堆积区、40～100 m 为载荷影响区、100～300 m 为压实稳定区[3]。进风巷和回风巷宽度为5 m。将采空区内岩石垮落堆积状态及其空隙视为多孔介质，相关模拟参数为：①空气摩尔质量：29 g/mol；②采空区空隙率：0.29；③动力黏度：2.01×10-5Pa·s；④气体扩散系数：2.88×10-2；⑤采空区区域：300 m×180 m；⑥自然堆积区渗透率：5.38×10-6m2；⑦载荷影响区渗透率：2.6×10-6m2；⑧压实稳定区渗透率：1.3×10-6m2；⑨CO2生成速率：1.026×10-3mol/（m3·s）；⑩CO 生成速率：3.4×10-4mol/（m3·s）。

根据实际情况，在COMSOL Multiphysics 数值模拟软件中设定采空区入口风速为0.9 m/s，分析采空区内气体流速分布情况。采空区气体流场示意图如图1。

图1 采空区气体流场示意图Fig.1 Distribution of gas flow field in goaf

从图1 可以发现：在流线密集处，采空区漏风风速较大；在流线稀疏处，采空区漏风风速较小；流线密集处主要集中在采空区的进风侧、回风侧以及采空区的浅部。为了更详细具体地分析采空区内气体流场的分布情况，进而根据采空区漏风风速划分采空区自燃“三带”的分布范围，在采空区中部沿工作面走向方向设置了1 条二维截线，借此分析采空区漏风风速在二维截线上各点的分布情况；二维截线的起点坐标为（0，80），终点坐标为（260，80）；二维截线的位置如图1 中红色线段。

根据国内外学者对采空区漏风的研究，采空区“三带”的范围可以根据采空区漏风流速划分[5]。采空区漏风风速沿二维截线的分布情况如图2。

图2 沿采空区二维截线（y=80 m）风流速度分布曲线Fig.2 Curve of gas velocity in goaf along two-dimensional section line（y＝80 m）

从图2 中可以看出：在0～40 m 之间风速较大，在0～40 m 范围内的遗煤由于散热速率大于产热速率，煤体氧化产生的热量大部分都被漏风携带走，不具备良好的蓄热条件，因此将这一区间划分为“散热带”；而在135 m 之后，漏风风速较小，煤体氧化产生的热量大部分留存在松散的煤体内，散热速率小于产热速率，具有良好的蓄热条件，但却由于漏风风速偏小，不能及时地为煤提供足量的氧气，造成煤的氧化反应进程受到抑制，因此将这一区间划分为“窒息带”；在40～135 m 之间的采空区遗煤由于漏风风速适中，既能提供较充足氧气供给煤的氧化反应进程，又具有较好的储热条件，因此将这一区间划分为“氧化升温带”。

1.2 采空区介观尺度气体流动规律

采空区某处遗煤及附近区域几何结构示意图如图3。

图3 采空区某处遗煤及附近区域几何结构示意图Fig.3 Schematic diagram of geometric structure of residual coal and nearby area in goaf

为了进一步研究介观尺度下颗粒聚团和颗粒内孔内部气体流动状态，选取采空区氧化升温带内某处（x=40 m，y=80 m）遗煤及附近区域作为研究对象。研究区域的长度为80 cm，高度为50 cm，假设其中的煤块为椭圆形，煤块长轴为24 cm，短轴为8 cm，与采空区底板的夹角为40°。实验测得孔隙率为9.5%，渗透率为1×10-16m2[6]。根据关联型多尺度研究方法，煤块作为颗粒聚团尺度，其入口风速由上一级尺度所决定，颗粒聚团处的风速为0.04 m/s。

模拟后得到颗粒聚团及附近区域气体流动的速度场和压力场分布情况如图4 和图5。

图4 采空区某颗粒聚团附近气体流动速度场分布示意图Fig.4 Distribution diagram of gas velocity field near a particle cluster in goaf

图5 采空区某颗粒聚团附近气体压力场分布示意图Fig.5 Distribution diagram of gas pressure field near a particle cluster in goaf

1.2.1 介观尺度垂直方向气体流动规律

为了进一步研究颗粒内孔尺度下，孔隙内部的气体在垂直方向的流动情况，在图3 中设置1 条垂直于底板且纵贯颗粒聚团的二维截线，函数表达式为x=40 cm（6 cm＜y＜14.5 cm），颗粒聚团内部垂直方向二维截线位置图如图6，垂直方向二维截线如图6中黄色线段。

图6 颗粒聚团内部垂直方向二维截线位置图Fig.6 Location diagram of two-dimensional vertical section line in a particle cluster

颗粒内孔气体流动速度沿垂直二维截线的分布曲线如图7。从图7 可以发现：越靠近颗粒聚团的表面，颗粒内孔中气体流动的速度越快。

图7 颗粒内孔气体流动速度沿垂直二维截线的分布曲线Fig.7 Distribution curve of gas velocity in pores along vertical two-dimensional section line

1.2.2 介观尺度水平方向气体流动规律

为了进一步研究颗粒内孔尺度下孔隙内部的气体在水平方向的流动情况，在图3 中设置1 条平行于底板且贯穿颗粒聚团的二维截线，函数表达式为y=7 cm（32 cm＜x＜42 cm），颗粒聚团内部水平方向二维截线位置示意图如图8，水平方向二维截线如图8 中黄色线段。

图8 颗粒聚团内部水平方向二维截线位置示意图Fig.8 Location diagram of two-dimensional horizontal section line in a particle cluster

颗粒内孔气体流动速度沿水平二维截线的分布曲线如图9。从图9 可以发现：越靠近颗粒聚团的迎风侧，颗粒内孔的气体流动速度越快。通过对介尺度层面颗粒聚团内部孔隙气体流动速度进行分析可以发现：颗粒内孔气体流动速度的单位大约在纳米级别，这和颗粒内孔的孔径是在一个尺度层面的，证明了模拟结果正确性。另外可以发现：越靠近颗粒聚团表面，越靠近颗粒聚团的迎风侧，颗粒孔隙中的气体流动速度越大。

图9 颗粒内孔气体流动速度沿水平二维截线的分布曲线Fig.9 Distribution curve of gas velocity in pores along horizontal two-dimensional section line

2 采空区气体流场对遗煤湿分的热蒸发作用分析

从图7 和图9 中可以看出：颗粒内孔中气体流动的速度大部分在10～20 pm/s 之间，且越靠近气-固交界面，气体流速越快。煤作为一种孔隙-裂隙双孔隙率多孔介质，孔隙-裂隙不仅是气体发生流动的主要结构，水分也主要在这些孔隙结构中进行储存和输运。滕腾[6]认为煤体裂隙中的残余水分在表面张力的作用下会镶嵌在迂曲基质块表面的沟堑和凹槽里，或者黏附在基质表面，也即裂隙表面；残余水分的存在不仅降低了基质表面的裸露面积，同时还减小了裂隙的有效开度；并且认为裂隙网络中的非连续残余水分会随着气体流动发生蒸发迁移。贺琼琼[7]通过采用NMR 方法分析褐煤中水分的赋存情况，发现胜利褐煤中的水分主要集中在0.003～0.3 μm 孔隙范围内，昭通褐煤中的水分主要集中在0.01～1 μm 孔隙范围内，褐煤中1 个有效亲水点位形成的水簇尺寸大约由6～7 个水分子组成；然后使用N2吸附差值法研究了介孔中水分的分布情况，发现大部分水分集中在26～100 nm 的孔隙中；大孔范围和裂隙并非储水结构，而是仅作为脱水过程中的传质通路。介观尺度下的气体也主要在这些传质通路中流动。介观尺度下煤中水分的储水及输运模型如图10。

图10 基质内孔与裂隙结构的“树杈”型储水及输运模型[2]Fig.10“Tree branch”type water storage and transport model of pore and fracture structure in matrix

当发生氧化蓄热升温现象时，颗粒内孔中的水分子会在热蒸发效应下发生相变成为气相水分子，之后通过扩散作用混合进入干燥气体中，并沿着输运通路进行运移。

3 热蒸发过程物理模型

采空区气体流场与遗煤之间的热湿传递过程主要涉及3 种物理过程：流体流动、传热、液体的相变以及与空气之间的输运。在自由流动域主要发生的是空气的流动、水分与空气之间的输运等物理效应；在多孔介质域主要发生的是水分的流动输运、相变吸热等物理效应。

3.1 液 相

在多孔介质域存在气相和液相2 种不同的流动相，且满足以下约束：

式中：Sg为气相的饱和度，对应于多孔介质中的湿空气；Sl为液相的饱和度，对应多孔介质中的水分。

多孔介质孔隙中气相为不饱和状态流动，所以使用Brinkman 方程计算多孔介质中湿空气的流场ug和压力分布pg。与湿空气速度相比，液相速度较小，因此根据气相压力梯度定义的Darcy 定律[8]可用于计算水流速度ul：

式中：κl为液相在煤体中的有效渗透率，m2；μl为液相的动力黏度系数，（N·s）/m2；▽pg为气相压力分布梯度，Pa/m。

多孔介质中的液相为不饱和流动，其有效渗透率与煤体绝对渗透率、相对渗透率、饱和度有关，因此式（2）可以进一步表示为：

式中：κ 为煤体绝对渗透率，m2；κrl为液相的相对渗透率；μw为水的动力黏度系数，（N·s）/m2；ε 为孔隙率。

液相的饱和度与水的浓度呈线性关系：

式中：cw为液相在煤体中的浓度，mol/m3；Mw为水的摩尔质量，g/mol；ρw为水的密度，kg/m3。

相对渗透率通常根据经验或者实验来确定，其形式在很大程度上取决于多孔材料属性和液体本身。根据参考文献［9］，将相对渗透率始终定义为正：

式中：κrg为气相的相对渗透率；eps 为浮点相对精度。

式中：Sli为残留在多孔介质中的液相的不可约饱和度。

3.2 水蒸气相

多孔介质中水蒸气的传递涉及2 个过程：第1个传递机制是根据Brinkman 方程求出的流体速度场，将其值赋给蒸汽相；第2 个传递机制是水蒸气及其气相干燥空气的二元扩散产生的速度分量。这2个组分之间的有效扩散系数Deff通常可以采用Millington & Quirk 方程进行描述：

式中：Dva为蒸汽-空气扩散率，取2.6×10-5m2/s。

将这2 种扩散效应得出的速度场应用于水蒸气传递方程，可以得到：

式中：uv为湿空气的速度场；ug为水蒸气的速度场；Ma为干燥空气的摩尔质量，g/mol；Mma为湿空气的摩尔质量，g/mol；ρma为湿空气的密度，kg/m3。

将湿空气视为理想气体，湿空气的密度和压强之间具有以下关系：

式中：ρg为气相的密度，kg/m3；pg为气相的压强，Pa；R 为理想气体常数，8.314 J/（mol·K）；Tg为气相的温度，K；Mg为气相的摩尔质量，g/mol。

3.3 相 变

为了计算蒸发到空气中的水量，并考虑到液态水的减少和湿空气中水蒸气浓度的增加，定义了以下计算式：

式中：mevap为蒸发量，mol/（m3·s）；K 为蒸发率，1/s；aw为蒸汽系数；csat为饱和情况下的蒸汽浓度，mol/m3；c 为当前蒸汽浓度，mol/m3。

将水蒸气视为理想气体，可得其饱和水蒸气浓度[10]为：

式中：csat为饱和水蒸气浓度，mol/m3；psat（T）为饱和水蒸气分压，Pa；T 为温度，K。

不同温度时所对应的饱和水蒸气分压近似为[7]：

对饱和蒸气分压的数学表达式进行数学分析可以发现其变化趋势主要取决于指数部分。将指数部分做如下数学变换：

当温度在［273.15 K，373.15 K］之间变化时，饱和蒸气分压psat（T）的变化区间为［610.7 Pa，100 kPa］。

4 采空区气体流场对遗煤的热蒸发效应数值模拟

为了研究采空区内气体流场和湿度场之间的相互关系，使用多物理场模拟软件，首先对采空区气体流动状态开展稳态研究，然后在稳态的基础上，研究热湿传递过程的瞬态变化。采空区内气体流动稳态结果如图4 和图5，分别表示采空区内气体流动的速度和压力分布情况。瞬态模拟的时间跨度为［0 h，24 h］，研究在此时间区间内受气体流动和蒸发吸热影响下，自由流动域和多孔介质域内的水蒸气浓度和湿度的变化情况。

4.1 采空区气体流场对遗煤中水分的影响

4.1.1 多孔介质域中水分含量的变化情况

多孔介质域水分含量随时间的变化情况如图11。多孔介质域中迎风侧某点水分含量的变化曲线如图12。多孔介质域内含水量随时间变化情况如图13。

图11 多孔介质域水分含量随时间的变化情况Fig.11 Variation of water content with time in porous media domain

图12 多孔介质域中迎风侧某点水分含量的变化曲线Fig.12 Variation curve of moisture content on the windward side in porous media domain

从图11 可以发现：煤中水分的蒸发首先从煤块外表面开始，之后逐步影响煤块内部；而且迎风侧比背风侧水分含量小。

从图12 中可以看到：观测点处的水分含量变化趋势为“下降-上升-下降”，且2 个下降阶段的曲线斜率大致相同，所以认为这2 个阶段中水分含量下降主要是由于热蒸发效应导致，中间阶段水分含量上升主要是受多孔介质内流场的影响，将迎风侧的水分向观测点处堆积迁移所导致的。

从图13 可以看出：通过对多孔介质域的水分含量分别在0 h 和24 h 进行积分运算，可以得到在［0 h，24 h］的时间区间内，水分含量由于蒸发效应减少了0.046 mol/m。这与刘剑等[11]之前使用热重分析得到的结果是相同的。即在煤常温氧化的早期阶段，煤重量的减少主要是由于煤中水分的蒸发所导致的。

4.1.2 自由流动域水蒸气浓度的变化情况

采空区内水蒸气浓度随时间变化情况如图14。

图14 采空区内水蒸气浓度随时间变化情况Fig.14 Variation of water vapor concentration in goaf with time

从图14 可以看出：在0.5 h 时自由流动域中水蒸气浓度最高处位于煤块背风侧风流尾部，在6、12、18 h 时自由流动域中水蒸气浓度最高处位于煤块的外表面，说明蒸发现象早期主要发生在煤块表层，导致了自由流动域气-固交界面处水蒸气浓度较低；后期水分才逐渐从煤块内部迁移至外表面处蒸发，导致了外表面的水蒸气浓度较高，说明在早期的表层蒸发和后期的水分迁移之间存在1 个启动过程。

4.2 热蒸发效应对采空区遗煤温度的影响

在多孔介质域，由于多孔介质域内的水分蒸发吸热，会使多孔介质域内的温度分布发生改变，多孔介质域中某点的蒸发吸热量Q 与蒸发量有关：

式中；Hevap为蒸发潜热，J/mol；mevap为蒸发量，mol/（m3·s）。而多孔介质域内的温度服从固体热传导规律，符合傅立叶传热方程。为了计算多孔介质域内的温度分布情况，需要将蒸发热作为源项插入到傅立叶传热方程中。需要注意的是，采空区氧化升温带内的遗煤在蓄热条件良好的情况下会将氧化蓄热积存在煤体内部，导致煤体温度缓慢升高，但这一过程不在研究范围内，此处仅考虑由于蒸发吸热导致的煤体温度变化情况。

多孔介质域蒸发吸热对煤体温度变化的影响如图15。

图15 多孔介质域蒸发吸热对煤体温度变化的影响Fig.15 Effect of evaporation and heat absorption in porous media on temperature of coal

模拟设定的环境温度为293.15 K，在2 h 时，煤块内部温度高于环境温度，在煤块的上部和下部温度较高。在6～12 h 时，煤块内部温度依然高于环境温度，但是从表层开始温度逐渐降低，煤块内部高于环境温度的区域越来越小。在18、24 h 时，煤块内部温度由于蒸发吸热，会出现温度低于环境温度的情况。

5 结 语

1）利用关联型多尺度研究方法，基于采空区气体流动在宏观尺度服从的N-S 方程和在介观尺度服从的达西定律，使用多物理场模拟软件研究了采空区流场速度和压力的跨尺度分布情况。

2）基于多孔介质模型总结了采空区遗煤在热蒸发效应下所涉及的气-液两相非饱和流动、相变、蒸发吸热等物理过程，构建了关于流体场、湿度场、温度场的多物理场耦合模型。并使用多物理场模拟软件研究了采空区多孔介质域孔隙内水分含量的变化情况以及自由流动域内水蒸气浓度的变化情况。

3）初步探讨了多孔介质域由于蒸发吸热导致的煤体温度分布变化情况。