采空区多参数流-固耦合温度场模拟研究＊

张九零 陈庆亚 王月红

（1.河北联合大学矿业工程学院，河北省唐山市，063009；2.矿山热动力灾害与防治教育部重点实验室，辽宁省阜新市，123000）

煤的自燃是环境散热和煤体放热相互作用的物理化学过程，煤自燃过程会产生大量的有毒有害气体，严重威胁着工作面附近工作人员的安全。煤在自燃的过程中，煤体内的各种气体浓度和含量、温度都在不断变化，且温度的分布不均匀，进而造成耗氧速率和放热强度的变化，进一步造成采空区高温区域的温度、位置随时间不断变化。近年来，随着科学技术的进步和煤矿管理水平的提高，煤矿自燃火灾百万吨发火率虽然每年有所下降，但是我国煤矿自燃火灾仍较为严重，对煤矿的安全生产构成严重威胁，所以研究影响采空区自然发火的因素是很有意义的。

1 基于有限体积法的采空区内动态温度场数学模型

根据采煤工作面随着开采推进不断移动的特点，依据有限体积法，利用传热学理论、能量守恒定律及多孔介质理论对采空区的热力分布进行分析，同时，由于采空区内遗煤和气体温度之间存在差异，故将采空区内的遗煤温度和气体温度进行耦合关联研究，建立移动坐标条件下的采空区温度场多参数流-固耦合数学模型，如式 （1）所示。移动坐标的引进将采空区温度场复杂的动态变化问题转变为静态问题进行处理，避免了非稳态方程组复杂、计算量较大的问题。

式中：vx——气流沿x轴方向的分量，m／s；

vy——气流沿y轴方向的分量，m／s；

λy——采空区内遗煤导热系数，W／m·℃；

λg——采空区内气体导热系数，W／m·℃；

tg——采空区内气体温度，K；

ts——采空区内遗煤温度，K；

ρg——采空区内气体密度，kg／m3；

ρs——采空区内遗煤密度，kg／m3；

cs——采空区内遗煤比热容，kJ／ （kg·K）；

cg——采空区内气体比热容，kJ／ （kg·K）；

Ke——对流换热系数，J／ （m2·s·K）；

Sn——单元体内遗煤与气体对流换热表面积，

m2；

q （t）——单位时间单位体积内采空区遗煤的放热量，kJ／ （m3·s）；

a——单元体的比表面积，1／m；

Δx——单元体的长，m；

Δy——单元体的宽，m；

Δz——单元体的高，m；

v0——工作面推进速度，m／s；

F——单元体的面积，m2；

Γ——单元体的边界。

2 采空区动态多参数温度场流-固耦合模拟

2.1 自燃模型数值计算范围

为了确定边界条件，需测定采空区内上下两巷位置的温度，利用传热范围边界来替代采空区实际边界，假定一个采空区的传热范围，则这个假定的传热范围边界属于绝热边界来处理，因为采空区的热量传递不只限于在采空区的实际边界内进行，也与采空区相邻的煤层顶底板、采空区四周的煤壁等进行热量交换，所以在采空区的边界温度或热流通量不能明确的情况下，采空区温度场的计算需要将采空区的事实边界进行扩大，将采空区温度的边界外推到几乎没有热流通量的地方，这样就可以定义温度场解算的边界条件 （热流通量为0）。根据经验将采空区自然发火模型的温度场边界条件分别向上和向下各推20m，以达到在此边界属于绝热边界，如图1中虚线所确定的范围为采空区温度场模型的解算范围。

图1 采空区坐标系及数值求解范围

2.2 网格划分

对采空区进行三角形网格剖分，为在实际情况中减少计算量，在保证程序计算精度前提下，对采空区采用疏密不一的网格划分方式，对采空区容易发生自燃的 “两道两线” （进风道、回风道、切眼线和停采线）处，因变化率较大 （本文未考虑停采的情况），进行网格加密处理，变化率较小的地方网格划分较疏散，在不影响计算精度的情况下减小了计算量。采空区划分后，共产生1000 个节点，划分为4851个小单元。

2.3 参数选取

工作面长度100m，工作面倾角0°，工作面通风阻力为15Pa，工作面掘进速度为2.4m／d，采空区遗煤平均厚度为0.7m。采空区遗煤原始温度30℃，工作面的进风温度25℃，煤的密度1400kg／m3，临界温度下煤的平均比热1200J／（kg·℃）。

2.4 主要采掘参数对采空区内温度场的影响

影响采空区温度的因素有很多，除了本身煤体的性质外，还和采煤方法及采煤工艺等因素有关。在采空区自然发火解算模型中提到了工作面长度、走向长度、推进速度、遗煤厚度等影响因素，由于篇幅有限，只对影响采空区温度场变化的工作面长度、推进速度、通风阻力3个采掘参数进行分析。

2.4.1 工作面长度对采空区温度的影响

为了研究采空区工作面长度对采空区温度场的影响，其他参数不改变，对工作面长度50 m、80 m、120m、150m、200m 下采空区的温度场进行解算，并用软件Tecplot10.0 对解算的数据处理，以图像的形式对采空区内气体温度进行显示，如图2所示。

图2 不同工作面长度时采空区温度场分布示意图

由图2 （a）可知，采空区气体温度沿着x、y轴的方向逐渐升高，且沿着x 轴方向温度变化剧烈。沿着x 轴方向在距离工作面100 m 处温度由开始设定的25℃快速升高到41℃，在距离采空区130 m 处温度逐渐降低到39℃，之后变化平稳。沿着y 轴方向温度逐渐升高，在距离回风口28m处温度达到最大值41℃。图2 （b）可知，采空区的温度规律和图2 （a）相似，都是沿着x、y 轴的方向逐渐升高，沿着x 轴方向温度变化剧烈，工作面长度为200m 时采空区的温度达到了45℃。

图3 不同工作面长度下采空区内温度分布示意图

图3为不同工作面长度下采空区温度分布示意图，可以看出，随着工作面长度的增加，温度在不断升高。由图3 （a）可以看出，在工作面长度为50m 时采空区出现的最高温度为41℃，随着距离工作面由20 m 到200 m 时，在距离工作面20 m的条件下工作面的长度为0 m 时，采空区的温度为25℃，到距离工作面200 m 的条件下工作面的长度为0m 时，采空区的温度升高为37℃。由图3（b）可以看出，在工作面长度为200 m 时采空区出现最高温度46℃，随着距离工作面由20 m 到200m 时，采空区温度由30℃升高到46℃，出现最高温度的位置向采空区深部移动。采空区最高温度随着工作面长度的增加而升高，虽然采空区出现最高温度的地点距离工作面由100 m 延后到了150m，但是采空区煤自燃发生的可能性还是很大，所以选择适当的工作面的长度可以在一定程度上减小煤自燃的几率。

工作面长度为50m 时，在采空区浅部130m处对温度影响最大，温度也达到最高，之后温度下降，在采空区深部超过150m 时，温度变化缓慢，温度的变化不再受工作面长度的变化。工作面长度为200 m 时，在采空区浅部150 m 处影响最大，在采空区深部超过150 m 时，温度变化缓慢，温度的变化不再受工作面长度的变化。由此可以得出，工作面长度越长，采空区内遗煤氧化的热量就越不容易排出，为采空区遗煤升温氧化自燃的发生提供了充足的时间。由此得出，工作面的长度越长，采空区遗煤的自燃速度就越快。

2.4.2 推进速度与采空区温度场相互关系

为了研究推进速度对采空区温度场的影响，在其他参数不变的情况下，增加采空区采煤的推进速度，速度由开始的1m／d分别增加到2m／d、3m／d、4m／d、5m／d。对采空区温度场进行解算，并用软件Tecplot 10.0对解算数据进行处理，对采空区内气体温度以图像的形式显示，如图4所示。

图4 不同推进速度下采空区温度场分布示意图

由图4 （a）可知，采空区气体温度沿着x、y轴的方向逐渐升高，且沿着x 轴方向温度变化剧烈。沿着x 轴方向在距离工作面104 m 处温度由开始设定的25℃快速升高到66℃，在距离采空区140 m 处温度逐渐降低到50℃，之后变化平稳。沿着y 轴方向温度逐渐升高，在距离回风口60m处温度达到最大值66℃。图4 （b）中，采空区气体温度沿着x、y 轴的方向逐渐升高，沿着x 轴方向在距离工作面150 m 处温度由开始设定的26℃升高到37℃，之后变化平稳。

图5 不同推进速度下采空区温度分布

图5 为不同推进速度下采空区温度分布示意图，在图5 （a）中，随着距离工作面越远，采空区的温度逐渐升高，温度由开始的28℃升高到66℃。在图5 （b）中，推进速度由1m／d增加到5 m／d时，采空区内气体出现的最高温度明显的降低，最高温度由66.7℃降低到37.8℃，出现最高温度的地点距离工作面的距离分别为104.5 m 和153m，速度越快则出现最高温度的地点离工作面越来越远。

推进速度为1 m／d 时，在采空区浅部104 m处对采空区内温度影响最大，温度也达到最高，之后温度下降，在采空区深部超过120 m 时，温度变化缓慢，温度的变化不再受工作面长度变化的影响。推进速度为5 m／d 时，在采空区浅部150 m处影响最大，在采空区深部超过160 m 时，温度变化缓慢，温度的变化不再受工作面长度变化的影响。

随着工作面推进速度的加快，采空区内最高温度急剧降低，推进速度对采空区最高温度影响较大，这是由于随着工作面移动，流进采空区低温气流和采空区遗煤增多，促进与采空区之间的热传导，对采空区起到冷却的作用，致使采空区温度降低，有效地降低了采空区煤自然发火的发生几率。

2.4.3 通风阻力对采空区温度场的影响

为了研究通风阻力对采空区温度场的影响，在其他参数不变的情况下，增加采空区的通风阻力，由开始的10Pa分别增加到20Pa、30Pa、40Pa、50Pa。利用采空区自然发火模型对采空区的温度场进行解算，用软件Tecplot10.0对解算数据进行处理，并对采空区内气体温度以图像的形式进行显示，如图6所示。

图6 不同通风阻力时采空区温度场分布示意图

由图6 （a）可见，随着通风阻力的增大，采空区的最高温度有所升高，采空区气体温度沿着x、y 轴的方向逐渐升高，且沿着x 轴方向温度变化剧烈。沿着x轴方向在距离工作面122m 处温度由开始设定的25℃快速升高到45℃，沿着y轴方向温度逐渐升高，在距离回风口48m 处温度达到最大值45℃。图6 （b）中采空区的温度沿着x、y轴的方向逐渐升高，在距离工作面160 m 处温度由开始设定的28℃快速升高到49℃，之后变化平稳。

图7 为不同通风阻力下采空区温度分布示意图，在图7 （a）显示随着距离工作面越远，采空区的温度逐渐升高，温度由开始的28℃升高到45℃，图7 （b）显示温度由开始的28℃升高到49℃。通风阻力由10Pa增加到50Pa时，采空区内气体出现的最高温度有所升高，由45℃升高到49℃，出现最高温度的地点距离工作面的距离由122m 延后到了158m。

图7 不同通风阻力下采空区温度分布示意图

通风阻力为10Pa时，在采空区浅部120m 处对温度影响最大，温度也达到最高，之后温度下降，在采空区深部超过130m 时，温度变化缓慢，温度的变化不再受工作面长度变化的影响。通风阻力为50Pa时，在采空区浅部160m 处影响最大，在采空区深部超过170 m 时，温度变化缓慢，温度的变化不再受工作面长度变化的影响。

虽然增大通风阻力后采空区内的温度增加的幅度较小，但通风阻力较大时，采空区的漏风量较大，增加了煤的氧化放热几率，使得采空区内的温度逐渐升高，且向离工作面远的方向移动。尽管通风阻力会使采空区出现最高温度的地点向远离工作面的方向移动，但是温度会逐渐升高，所以为了避免煤自燃的发生，可以选择合适的通风阻力来保障安全生产。

3 结论

（1）通过有限体积法 （FVM）来建立采空区内动态温度场多参数耦合数学模型，将动态问题转化为静态问题，避免了计算量过大的问题，并且确定了采空区自然发火模型计算范围，将采空区自然发火模型的温度场边界条件分别向上和向下各推20m来达到绝热边界。对采空区进行三角形划分后选取参数，利用Tecplot10.0软件对解算数据进行处理，并以图像显示采空区温度场的变化情况，使得结果显示简单明了。

（2）当工作面长度增加时，采空区内的最高温度升高，虽然采空区出现最高温度的地点逐渐远离工作面，但是出现自燃的可能性还是很大，所以选择适当的工作面长度可以有效降低采空区内温度，防止自燃事故的发生。

（3）当采空区采煤的推进速度增加时，采空区出现的最高温度逐渐变低，由1 m／d到5 m／d温度变化较为显著，当推进速度超过5 m／d时，温度变化较小，而且出现最高温度的地点一直远离工作面，所以要保证生产，降低采空区煤自燃的发生必须选择合适的推进速度。

（4）随着工作面通风阻力的增加，采空区出现的最高温度有所升高，并且采空区出现高温的地点向远离工作面的方向移动，所以选择适当的通风阻力也可以有效保证煤矿的安全生产。

［1］ 纪玉龙，时国庆等.综放面采空区瓦斯抽放条件下氧浓度场的CFD 模拟 ［J］.中国煤炭，2012 （4）

［2］ 徐精彩，张辛亥等.煤层自燃胶体防灭火理论与技术 ［M］.北京：煤炭工业出版社，2003

［3］ 罗海珠，梁运涛.煤自然发火预测预报技术的现状与展望 ［J］.中国安全科学学报，2003 （3）

［4］ 邓军，徐精彩等.国内外煤炭自然发火预测预报技术综述 ［J］.西安矿业学院学报，1999 （4）

［5］ 李宗翔，韦涌清等.非均质采空区气-固耦合温度场风有限元求解 ［J］ .昆明理工大学学报，2004（2）

［6］ 邬剑明，迟克勇等.煤自然发火过程的模拟实验研究 ［J］.中国煤炭，2008 （4）

［7］ 张国枢.通风安全学 ［M］.徐州：中国矿业大学出版社，2000

［8］ 张国枢，戴广龙.煤炭自燃理论与防治实践 ［M］.北京：煤炭工业出版社，2002

［9］ 牛立东，王海平 .水地庄煤矿井下自燃案例探讨［J］.中国煤炭，2011 （7）

［10］ 王同科.一维二阶椭圆和抛物线微分方程的高精度有限体积元方法 ［J］ .数值计算与计算机应用，2002 （4）

［11］ 尹哲 .多孔介质渗流问题的对称有限体积方法［D］.山东：山东大学，2007