基于FLUENT 动网格的煤矿采空区三带动态特性模拟技术

王芳

（西安科技大学，陕西 西安710000）

采空区自燃三带分布情况与工作面的推进有着密切联系。工作面推进过程中，采空区边界的移动不仅使采空区实体的几何尺寸发生了连续的动态变化，同时也影响了采空区内部的漏风供氧条件、各组分气体传输和三带分布情况。工作面的推进首先造成采空区几何边界的移动，其次由于回采工作面的两端压差提供漏风动力[1]，采空区漏风边界也随之变化并直接改变了采空区氧气浓度的传输边界和分布范围[2]，体现为采空区流场、氧浓度场等多场的动态变化[3]，因此对采空区自燃三带的数值模拟应该在工作面动态推进情况下进行模拟[4]。

1 采空区自燃三带数值模拟

由于采空区是一个封闭的区域，基本上无法进入，对采空区各场的变化规律进行数值模拟是研究采空区煤自燃危险的重要方式之一。结合采空区的实际情况，建立数学模型，并利用计算机对数学模型进行解算，分析研究采空区自燃“三带”的分布规律是是数值模拟研究采空区的一般步骤。其中，CFD（Computational Fluid Dynamics，计算流体动力学）的方法是先建立一套描述采空区风流流场、温度场及氧浓度场的数学模型，以现场实测数据为已知条件，选择合理的离散方法和计算程序，通过模拟采空区漏风流场、温度场和氧浓度场，量化分析采空区煤自燃在受到不同边界条件影响下的分布规律，最后对计算结果进行相应的处理并显示输出[5]。CFD 中的FLUENT 作为求解器，采用了多种求解方法和多重网格加速收敛技术，有较高的收敛速度和求解精度。使用FLUENT 进行求解前，应先建立模型、划分网格，根据模拟要求选择稳态或非稳态模式，选择需要解的方程（流动方程、能量方程、组分运输方程等），确定需要附加的模型（多孔介质模型），在对边界条件及材料物理性质进行设置后，对模型进行求解。对工作面动态回采情况下采空区进行数值模拟，是通过模拟开采进行一段时间内的采空区内各场变化情况实现的，各场的变化与时间相关，可以采用非稳态模式进行计算。工作面动态采情况下的采空区范围随着工作面的移动而不断扩大，采用FLUENT 中的动网格模型来实现采空区工作面动态回采，并且建立孔隙率、煤氧反应速率等与时间相关的函数通过用户自定义功能编入FLUENT 进行计算。

2 FLUENT 动网格技术

采空区范围随工作面推进而不断扩大的过程是计算区域不断扩大的过程，在Fluent 中解决该问题的方法通常为采用动网格模型。Fluent 中的动网格模型可用来模拟由于流体域边界刚性运动或者边界变形引起的流体域形状随时间变化的问题，每一个时间步上的体网格的更新是由解算器根据边界的新的位置确定的，根据边界的运动自动调节内部体网格的分布，通过指定运动区域的运动方式以及网格再生的方法可以使网格动态变化。

对于边界移动的任意控制体积V 上的一般标量的 φ的守恒方程可写为：

式中，V 为空间中大小和形状都随时间变化的控制体积，∂V为控制体积的运动边界，ug为运动网格的运动速度；ρ为流体密度；u 为流体速度；г 为耗散系数；S φ 是标量 φ的源项。

在FLUENT 中，每一个时间步上体网格的更新是由解算器根据边界的新的位置来自动完成的，即解算器可以根据边界的运动和变形自动地调节内部体网格节点的分布，使用起来非常地方便。在使用动态网格时，用户仅需提供初始网格并在模型中指定任意区域的运动即可。关于运动的指定，FLUENT 允许用户通过边界profile 文件或者用户自定义函数(UDF)或者六自由度( 6DOF )解算器来指定。当边界运动函数较复杂时，一般利用UDF 宏指定运动，然而当需要指定的边界为简单的速度- 时间关系运动时，利用Profile 文件指定运动则更为方便高效。

3 采空区自燃三带数值模拟

3.1 模拟过程

以某煤矿为例，对煤矿原型进行简化，建立模型，如图1 所示。配风量为36 m3/s，进回风巷宽4m，高4m，工作面斜长200 m，平均推进速度为4m/d。模拟研究时空范围是工作面从已回采了200 m 的位置开始之后的45 天回采期，到时采空区总长度为380m。

图1 采空区模型

所建模型含有两个部分，工作面巷道和采空区。因此划分网格时选择Multizone 网格划分方法进行划分，Multizone 可以对目标区域进行自动分区，将目标区域自动分解成多个可以扫掠或是自由划分的区域，只需要简单的指定源面、设置网格控制参数等，即可对目标区域进行得到高质量的网格。相较于传统的分割、扫掠生成六面体单元的网格划分方法，Multizone 网格划分方法省略了分割步骤，只需要进行适当参数设置即可生成高质量的六面体网格。

网格划分质量的高低直接影响了FLUENT 结果的精度以及收敛的速度，判断网格质量的高低主要依据两个指标：Element quality 和Skewness。Element quality 基于给定单元的体积与边长之间的比率，其值处于0 到1 之间，越接近于1 越好；Skewness是最基本的网格质量检查项，与Element quality 相反的是，越接近于0 表示网格质量越高。网格完成后的模型如图2 所示，共划分100512 个六面体网格，Element quality 值为0.99，Skewness 值为0.00000004，由此可见网格质量较好，可以满足计算精度的要求。

图2 采空区模型网格划分

工作面推进的过程即可视为简单的速度- 时间运动，工作面的移架、放顶虽然不是连续进行的，但是从一个很长的时间周期角度来说，开采过程中每日进刀数是较为固定的，因此可以将工作面相对于采空区的运动看作为匀速运动，可以通过Profile 文件来描述。

如回采速度为4m/d，即0.000046m/s，Profile 文件如下：

((v_y transient 4 0)

(time 0 1 2 4)

(v_y -0.000046 -0.000046 -0.000046 -0.000046))

网格再生方法主要用来计算内部网格节点的调节，采用铺层（Layering）与局部重构（Remeshing）相结合的方法进行计算。动态模型的模拟采用瞬态模式，根据所需要模拟的时长来设置时间步长以及时间步数。

时间步长大致由以下公式估计：

3.2 数值模拟结果

3.2.1 稳态结果

在不使用动网格方法的情况下进行稳态模拟，得到采空区氧气浓度场的分布情况如图3。按照氧气体积浓度指标将模拟所得氧气浓度场划分成为的自燃“三带”[6]，如图4。氧气浓度大于15%的范围为散热带，漏风较大，不利于热量累积，氧气浓度为5%～15%的范围为氧化升温带，氧气浓度小于5%为窒息带。可以看出进风侧的氧化升温带宽度明显比回风侧的宽。

3.2.2 瞬态结果

使用动网格方法进行瞬态模拟，模拟采空区从200m 的初始长度推进至380m 的工作面动态推进过程，得到采空区氧气浓度场的演变情况如图5。初期采空区内氧气浓度场变化较明显，但随着开采进行，言其浓度场逐渐趋于稳定，不再随时间发生明显变化。在氧气浓度场不再随时间发生明显变化后，按照氧气体积浓度指标划分成为的自燃“三带”，如图6 所示。

图3 采空区氧气浓度场

图4 采空区自燃“三带”

图5 工作面动态回采过程中的氧气浓度场演变

图6 采空区自燃“三带”分布图

3.2.3 结果分析

通过对稳态和瞬态结果进行对比可知，两种模拟方法所得到的结果有一定差距。相较于稳态计算结果，瞬态计算出的氧化升温带的范围更宽，而使用稳态方法进行计算的结果比较接近于使用瞬态方法的计算初期的结果。

由于瞬态计算即是在每个时间步内进行近似稳态计算，计算完整个计算域后再将所得的结果作为下一个时间步的初始值进行近似稳态计算，直至时间步的结束。稳态计算只需要最终迭代达到收敛即可，而瞬态计算则要求每个时间步内均达到收敛。若只考虑系统稳定后的状态，那么就选择使用稳态计算；若要考虑流场的演化情况，则需要使用瞬态。在实际研究中，选择稳态方法还是瞬态方法取决于研究目标的侧重点使用稳态计算的采空区自燃“三带”的分布结果，表达了采空区各场处于稳定状态时的情况，而使用FLUENT 动网格的瞬态计算，则更侧重于工作面回采过程中氧气浓度场的演变情况。

4 结论

4.1 通过对稳态和瞬态结果进行对比可知，两种模拟方法所得到的结果有一定差距。使用稳态方法进行计算的结果比较接近于使用瞬态方法的计算初期的结果。

4.2 与稳态计算结果相比，使用FLUENT 动网格的瞬态计算，则更侧重于工作面回采过程中氧气浓度场的演变情况。