工程科学学报，第38卷，第10期：1350-1358,2016年10月 Chinese Journal of Engineering,Vol.38,No.10:1350-1358,October 2016 D0l:10.13374/j.issn2095-9389.2016.10.002:http://journals..ustb.edu.cn 非均质孔隙率采空区氧化升温规律四维动态模拟 周佩玲”，张英华”，黄志安区，袁飞2》，高玉坤)，王辉”，孙倩” 1)北京科技大学金属矿山高效开采与安全教有部重点实验室，北京1000832)北京科技大学治金与生态工程学院，北京100083 3)北京科技大学钢铁治金新技术国家重点实验室，北京100083 ☒通信作者，E-mail:huang-za@qg.com 摘要运用Fuet动网格模型实现采空区的四维动态变化，并用用户自定义函数将煤低温氧化动力学机理及非均质孔隙 率函数编入Fut中，结合时间和空间，对U+L型通风系统采空区升温规律进行四维动态模拟研究.研究表明：非均质孔隙 率四维动态模型能更真实地反应孔隙率的空间与时间变化，空间某一位置的孔隙率随时间呈负指数递减：工作面推进速度越 大，采空区升温速率越小，推进速度为3.6m"d1时平均升温速率仅为推进速度为1.2m·d时的1/5：然而，推进速度越大，高 温点的深度越大，不利于自燃的预防：尾巷的存在使得温度场范围扩大，温度升高，C0主要从尾巷流出，尾巷释放的C0量是 回风巷C0释放量的10倍.最后利用现场实测的数据对结果进行验证，表明模拟结果是正确可信的. 关键词采空区：非均质孔隙率：煤燃烧：氧化升温：四维：动态模拟 分类号TD75·2 4D dynamic simulation of coal oxidation heating law in gobs with heterogeneous porosity ZHOU Pei--Hing》,ZHANG Ying-hua》,HUANG Zhi--ane,YUAN Fei2,3》,GA0Yu-hun',WANG Hui,SUN Qian” 1)Key Laboratory of High-Efficient Mining and Safety of Metal Mines (the Ministry of Education),University of Science and Technology Beijing, Beijing 100083,China 2)School of Metallurgical and Ecological Engineering,University of Science and Technology Beijing,Beijing 100083,China 3)State Key Laboratory of Advanced Metallurgy,University of Science and Technology Beijing,Beijing 100083,China Corresponding author,E-mail:huang_za@qq.com ABSTRACT Using the Fluent dynamic mesh model to realize the 4D dynamic movement of a gob and inputting the kinetic mecha- nism of coal low-temperature oxidation and the dynamic change of heterogeneous porosity into Fluent through the user defined function, 4D dynamic simulation is performed on the spontaneous heating law in a gob for U+L ventilation of a certain mine.The results show that unsteady heterogeneous porous media can be more realistic in response to the change of porosity,which exponentially decreases with time.The greater the advancing speed of the working face,the smaller the heating rate,and the average heating rate at an advan- cing speed of 3.6md is 1/5 of that at an advancing speed of 1.2m'd.However,the higher the advancing speed,the deeper the depth of the high temperature region,which is unfavorable for the prevention from spontaneous combustion of coal.Because of the existence of a tail roadway in the gob,the temperature field expands,the temperature rises,the main way to release CO is the tail roadway,and the amount of CO from the tail roadway is 10 times that from the outlet.Finally,the results are verified by the field test data,indicating that the simulation results are correct and reliable. KEY WORDS gobs;heterogeneous porosity:coal combustion:oxidation heating:four-dimensional:dynamic simulation 收稿日期：201603-29 基金项目：国家自然科学基金资助项目(51474017)：新疆维吾尔族自治区基金资助项目(2014211013)
工程科学学报，第 38 卷，第 10 期: 1350--1358，2016 年 10 月 Chinese Journal of Engineering，Vol． 38，No． 10: 1350--1358，October 2016 DOI: 10． 13374 /j． issn2095--9389． 2016． 10． 002; http: / /journals． ustb． edu． cn 非均质孔隙率采空区氧化升温规律四维动态模拟 周佩玲1) ，张英华1) ，黄志安1) ，袁 飞2，3) ，高玉坤1) ，王 辉1) ，孙 倩1) 1) 北京科技大学金属矿山高效开采与安全教育部重点实验室，北京 100083 2) 北京科技大学冶金与生态工程学院，北京 100083 3) 北京科技大学钢铁冶金新技术国家重点实验室，北京 100083 通信作者，E-mail: huang_za@ qq． com 摘 要 运用 Fluent 动网格模型实现采空区的四维动态变化，并用用户自定义函数将煤低温氧化动力学机理及非均质孔隙 率函数编入 Fluent 中，结合时间和空间，对 U + L 型通风系统采空区升温规律进行四维动态模拟研究． 研究表明: 非均质孔隙 率四维动态模型能更真实地反应孔隙率的空间与时间变化，空间某一位置的孔隙率随时间呈负指数递减; 工作面推进速度越 大，采空区升温速率越小，推进速度为 3. 6 m·d － 1时平均升温速率仅为推进速度为 1. 2 m·d － 1时的 1 /5; 然而，推进速度越大，高 温点的深度越大，不利于自燃的预防; 尾巷的存在使得温度场范围扩大，温度升高，CO 主要从尾巷流出，尾巷释放的 CO 量是 回风巷 CO 释放量的 10 倍． 最后利用现场实测的数据对结果进行验证，表明模拟结果是正确可信的． 关键词 采空区; 非均质孔隙率; 煤燃烧; 氧化升温; 四维; 动态模拟 分类号 TD75 + 2 收稿日期: 2016--03--29 基金项目: 国家自然科学基金资助项目( 51474017) ; 新疆维吾尔族自治区基金资助项目( 2014211B013) 4D dynamic simulation of coal oxidation heating law in gobs with heterogeneous porosity ZHOU Pei-ling1) ，ZHANG Ying-hua1) ，HUANG Zhi-an1) ，YUAN Fei2，3) ，GAO Yu-kun1) ，WANG Hui1) ，SUN Qian1) 1) Key Laboratory of High-Efficient Mining and Safety of Metal Mines ( the Ministry of Education) ，University of Science and Technology Beijing， Beijing 100083，China 2) School of Metallurgical and Ecological Engineering，University of Science and Technology Beijing，Beijing 100083，China 3) State Key Laboratory of Advanced Metallurgy，University of Science and Technology Beijing，Beijing 100083，China Corresponding author，E-mail: huang_za@ qq． com ABSTＲACT Using the Fluent dynamic mesh model to realize the 4D dynamic movement of a gob and inputting the kinetic mechanism of coal low-temperature oxidation and the dynamic change of heterogeneous porosity into Fluent through the user defined function， 4D dynamic simulation is performed on the spontaneous heating law in a gob for U + L ventilation of a certain mine． The results show that unsteady heterogeneous porous media can be more realistic in response to the change of porosity，which exponentially decreases with time． The greater the advancing speed of the working face，the smaller the heating rate，and the average heating rate at an advancing speed of 3. 6 m·d － 1 is 1 /5 of that at an advancing speed of 1. 2 m·d － 1 ． However，the higher the advancing speed，the deeper the depth of the high temperature region，which is unfavorable for the prevention from spontaneous combustion of coal． Because of the existence of a tail roadway in the gob，the temperature field expands，the temperature rises，the main way to release CO is the tail roadway，and the amount of CO from the tail roadway is 10 times that from the outlet． Finally，the results are verified by the field test data，indicating that the simulation results are correct and reliable． KEY WOＲDS gobs; heterogeneous porosity; coal combustion; oxidation heating; four-dimensional; dynamic simulation
周佩玲等：非均质孔隙率采空区氧化升温规律四维动态模拟 ·1351· 采空区遗煤自燃是煤矿灾害的主要来源之一，严 尼乌斯活化能，由化学动力学实验可以得到：T是温 重威胁煤矿安全生产-.采空区遗煤的自燃是受漏 度，K;n是反应级数，煤和其他含碳材料低温氧化，n 风、氧化放热、瓦斯涌出、工作面推进等多因素作用的 的取值范围为0.5~1，取一般经验值为1，这也符合我 结果.其中，工作面的动态推进是影响采空区遗煤自 国的煤体参数.因此，煤氧反应的速率方程为式 燃的重要因素.目前，国内外学者在煤氧复合反应 (2).通过用户自定义函数开源编程将反应速率加载 的基础上，建立了三维采空区自燃数学模型，运用多 至Fluent多场耦合数学模型中. 场耦合的方法，研究了U型通风三维采空区自然发火 r=ACo.exp (-E/RT). (2) 规律.李宗翔等B,刀开发了G3程序对采空区进行动态 煤氧反应放出热量，采空区散热量主要是由于漏 模拟，描述了U型采空区工作面推进过程中采空区内 风流的存在，受进回风巷压差及风量的影响，煤氧低温 漏风流态等分布规律；时国庆等圆利用Fluent中的用 氧化初期，则能量输运方程四为 户自定义函数(UDF)导入随时间变化的动属性，在定 网格的条件下对采空区的氧浓度场进行数值模拟；对 epn()p. dt 于有尾巷存在的U+L型通风采空区的研究，杨胜强 和邵昊等9-0运用Fluent软件就尾巷对采空区流场和 +"哥+u）= dy 瓦斯涌出的影响方面进行了研究 上述对采空区流场或温度场的研究都建立在静 ()+0c,n. (3) 态的基础上，有关于工作面推进的研究也是假设采 式中：是孔隙率：P和p.分别是气体和固体的密度， 空区实体是静态的，与时间相关的属性只是在定区 gcm3:c和ce分别是气体和固体的比热容，J·g· 域内的变化，即可研究采空区内固相煤体并没有变 K1;1为时间，s;入是煤基有效导热系数，入，和入.分别 化，是实际情况的一种假设：同时，研究对象都是· 是空气和煤的导热系数，J小cmsK:Q(C,T)是 型通风采空区，对于+L通风采空区的研究仅做到 煤氧反应的反应热，与温度和氧气浓度及比表面积有 对流场的稳态研究，没有对升温过程的进一步研究 关，所以有 成果.然而，实际的采空区是受空间和时间动态变化 入m=e入.+(1-e）·A.,Q(C,T)=H..C.r-S 影响的非稳态结果，工作面后煤岩不断垮落，采空区 其中H,为标准生成焓，kJ·mol:对于单个球形小颗粒 不断扩大，孔隙率不仅在空间上具有非均质性，而且 随时间不断变化，从而影响采空区内的传热和传质 的比表面积5。==合D是粒子的直径，V和S分 过程.工作面的动态推进对温度场及氧浓度场的影 别是理想球体的表面积和体积，而对于单元体而言，比 响最大，而对其他场的影响可以忽略不计2-切.因 此，本文基于山西某矿的现场要求，利用Fluent中的 表面积=1-8)二， ,所以S=61-ea D 动网格模型，实现采空区工作面的动态推进，该模型 煤矿现场广泛采用C0作为煤自燃早期识别与预 具备动网格动属性及动边界条件：同时，运用用户自 测的标志性气体.根据传质学理论，采空区内气体组 定义函数将煤氧反应动力学机理及孔隙率与空间、 分方程为国： 时间的函数导入luent中，建立具有尾巷的U+L型 a(pC,) div (puC,)div (D,grad (pC,))+P.(4) 采空区，对该矿孔隙率的非均质性和动态推进采空 at 区氧化升温过程进行研究：最后，通过现场测温数据 式中，C为s组分气体的质量分数：u是流经单元体 对模拟结果进行验证 的气体流速，m·s;D,为该组分的扩散系数，cm2· 1采空区非稳态四维动态模型 s:P,为该组分在体系内单位时间单位体积的质量 变化率（即0，或C0的变化量）.通过用户自定义函 1.1采空区煤低温氧化动力学数学模型的建立 数将P,作为源项编入数学模型中进行解算，本文将 煤氧发生低温氧化，消耗氧气，产生热量，产生 煤与氧气在低温条件下的氧化视为主要是化学反 C0,C0,等，是一个复杂的氧化动力学过程4.根 应，所以O2的消耗和C0的产生符合以下公式： 据煤氧化学反应动力学中的反应速率方程： Souce_O2 =-0.032r:Source_CO =0.056r,r(1) r=ACo.exp (E/RT). (1) 的反应速率. 式中：化学反应速率r是反应物的函数，kmol·m3· 为了研究采空区遗煤的低温氧化规律，假定采空 s:A是反应指前因子，C。,是氧气质量分数；E为阿累 区是由破碎的岩体和遗煤组成，数学模型为：
周佩玲等: 非均质孔隙率采空区氧化升温规律四维动态模拟 采空区遗煤自燃是煤矿灾害的主要来源之一，严 重威胁煤矿安全生产［1--2］． 采空区遗煤的自燃是受漏 风、氧化放热、瓦斯涌出、工作面推进等多因素作用的 结果． 其中，工作面的动态推进是影响采空区遗煤自 燃的重要因素［3--5］． 目前，国内外学者在煤氧复合反应 的基础上，建立了三维采空区自燃数学模型［6］，运用多 场耦合的方法，研究了 U 型通风三维采空区自然发火 规律． 李宗翔等［3，7］开发了 G3 程序对采空区进行动态 模拟，描述了 U 型采空区工作面推进过程中采空区内 漏风流态等分布规律; 时国庆等［8］利用 Fluent 中的用 户自定义函数( UDF) 导入随时间变化的动属性，在定 网格的条件下对采空区的氧浓度场进行数值模拟; 对 于有尾巷存在的 U + L 型通风采空区的研究，杨胜强 和邵昊等［9--11］运用 Fluent 软件就尾巷对采空区流场和 瓦斯涌出的影响方面进行了研究． 上述对采空区流场或温度场的研究都建立在静 态的基础上，有关于工作面推进的研究也是假设采 空区实体是静态的，与时间相关的属性只是在定区 域内的变化，即可研究采空区内固相煤体并没有变 化，是实际情况的一种假设; 同时，研究对象都 是 U 型通风采空区，对于 U + L 通风采空区的研究仅做到 对流场的稳态研究，没有对升温过程的进一步研究 成果． 然而，实际的采空区是受空间和时间动态变化 影响的非稳态结果，工作面后煤岩不断垮落，采空区 不断扩大，孔隙率不仅在空间上具有非均质性，而且 随时间不断变化，从而影响采空区内的传热和传质 过程． 工作面的动态推进对温度场及氧浓度场的影 响最大，而对其他场的影响可以忽略不计［12--13］． 因 此，本文基于山西某矿的现场要求，利用 Fluent 中的 动网格模型，实现采空区工作面的动态推进，该模型 具备动网格动属性及动边界条件; 同时，运用用户自 定义函数将煤氧反应动力学机理及孔隙率与空间、 时间的函数导入 Fluent 中，建立具有尾巷的 U + L 型 采空区，对该矿孔隙率的非均质性和动态推进采空 区氧化升温过程进行研究; 最后，通过现场测温数据 对模拟结果进行验证． 1 采空区非稳态四维动态模型 1. 1 采空区煤低温氧化动力学数学模型的建立 煤氧发生低温氧化，消耗氧气，产生热 量，产 生 CO、CO2 等，是一个复杂的氧化动力学过程［14--15］． 根 据煤氧化学反应动力学中的反应速率方程: r = ACn O2 exp( － E /ＲT) ． ( 1) 式中: 化学 反 应 速 率 r 是 反 应 物 的 函 数，kmol·m － 3· s － 1 ; A 是反应指前因子，CO2 是氧气质量分数; E 为阿累 尼乌斯活化能，由化学动力学实验可以得到; T 是温 度，K; n 是反应级数，煤和其他含碳材料低温氧化，n 的取值范围为 0. 5 ～ 1，取一般经验值为 1，这也符合我 国的煤体参数［16］． 因此，煤氧反应的速率方程为式 ( 2) ． 通过用户自定义函数开源编程将反应速率加载 至 Fluent 多场耦合数学模型中． r = ACO2 exp( － E /ＲT) ． ( 2) 煤氧反应放出热量，采空区散热量主要是由于漏 风流的存在，受进回风巷压差及风量的影响，煤氧低温 氧化初期，则能量输运方程［12］为 ［ερg cpg + ( 1 － ε) ρc cpc ］T t + ρg cpg ( μ T x + ν T y + ω T ) z = λeff ( 2 T x 2 + 2 T y 2 + 2 T z 2 ) + Q( C，T) ． ( 3) 式中: ε 是孔隙率; ρg和 ρc分别是气体和固体的密度， g·cm － 3 ; cpg和 cpc分别是气体和固体的比热容，J·g － 1· K － 1 ; t 为时间，s; λeff是煤基有效导热系数，λg和 λc分别 是空气和煤的导热系数，J·cm － 1·s － 1·K － 1 ; Q( C，T) 是 煤氧反应的反应热，与温度和氧气浓度及比表面积有 关，所以有 λeff = ελg + ( 1 － ε)·λc，Q( C，T) = Hs ·C·r·Sex ． 其中 Hs为标准生成焓，kJ·mol － 1 ; 对于单个球形小颗粒 的比表面积 Sex = S' V' = 6 D ，D 是粒子的直径，V'和 S'分 别是理想球体的表面积和体积，而对于单元体而言，比 表面积 S V = ( 1 － ε)·S' V'，所以 Sex = 6( 1 － ε) D ［14］ ． 煤矿现场广泛采用 CO 作为煤自燃早期识别与预 测的标志性气体． 根据传质学理论，采空区内气体组 分方程为［13］: ( ρCs ) t + div( ρuCs ) = div( Dsgrad( ρCs ) ) + Ps ． ( 4) 式中，Cs为 s 组分气体的质量分数; u 是流经单元体 的气体流 速，m·s － 1 ; Ds 为该 组 分 的 扩 散 系 数，cm2 · s － 1 ; Ps为该组分在体系内单位时间单位体积的质量 变化率( 即 O2 或 CO 的变化量) ． 通过用户自定义函 数将 Ps作为源项编入数学模型中进行解算，本文将 煤与氧气在低温条件下的氧化视为主要是化学反 应，所 以 O2 的 消 耗 和 CO 的 产 生 符 合 以 下 公 式: Souce_O2 = － 0. 032r; Source_CO = 0. 056r，r 为( 1) 中 的反应速率． 为了研究采空区遗煤的低温氧化规律，假定采空 区是由破碎的岩体和遗煤组成，数学模型为: · 1531 ·
·1352· 工程科学学报，第38卷，第10期 [+0+0=0, r0.2e-a23w+m+0.1）(e-a15s+1）x1.05, x≤L2: (0.2e-a23r+m+0.1）(e-asu-0+)x1.05,x>l/2. aCo+Q. aCo+Q ay +0C= (6) 式中：是工作面走向长度，m:初始孔隙率空间分布函 )-m 数如图1所示，其中e,e,和e分别是沿工作面走向(x 轴)，采空区走向(y轴)和采空区高度方向(z轴)上的 -cn+(g++)- 孔隙率，孔隙率分布不仅是空间坐标的函数，同时是时 间的函数，空间某位置的孔隙率随着工作面的推进不 (@盟+0+0) 断减小.由Carman-Kozeny公式得出，渗透率k与孔隙 是(k兴)+(K兴)+是(K）=0 率ε的关系为 do 62 (5) k=0.24(1-8 (7) 该多元方程组中，x、y、z为三维空间坐标，Q为漏 式中，d。是粒径的平方，量级可取l0°m2.将孔隙率代 风流量，K为绝对渗透率，假设多孔介质为各向同性， 入其中可以得到渗透率的空间分布及随时间的变化规 根据实验数据K=K,=K=1×10-6m2,Co,为02的 律，使用用户自定义函数中DEFINE_PROFILE将孔隙 质量分数，取24%，D.为氧气在煤体中的扩散系数， 率及渗透率编入Fluent进行解算. V(T)为实际耗氧速率，mol·sl·cm3;H为热力分压， Pa;Q(C,T)为煤氧化反应放热源项- 0.6 0.5 综上所述，影响煤氧反应速率的因素主要有温 0.4 度、氧含量及孔隙结构，化学反应中组分的变化及放 不03 0.2 热量又是反应速率的函数.建立U+L型通风方式 0. 的采空区几何模型，将煤氧反应的机理编写成用户 工作面 050100150200250050100150200250 自定义函数程序，导入模拟软件Fluent中，为数值计 x/m 算增加化学反应的组分源项和能量源项，研究采空 图1初始孔隙率空间分布函数曲面 区的氧化升温规律，力求减小计算误差，提高仿真模 Fig.I Initial porosity function space distribution 拟的真实性. 本文以山西某矿U+L通风采空区为原型建立几 1.2U+L通风采空区四维动态模型的建立 何模型.工作面长度为200m,采空区走向取250m, 采空区工作面推进是一个动态过程，使用计算流 进、回风巷及联络巷宽4m,高2.5m,长10m.其中4 体力学(CD)动网格的方法，进行网格动态重构，建立 煤层赋存于山西组下部，煤厚0.84~6.05m,平均 采空区四维动态模型，使采空区范围不断扩大.使用 2.98m,含夹矸1~4层.4和3煤层之间的距离为 用户自定义函数中宏RP_Get_Real“flow-time”获取工 9.32m,3和2煤层之间的距离为6.94m,属于近距离 作时间，并以此建立工作面推进速度、孔隙率和渗透率 煤层.根据该矿的实际情况，结合模拟要求作以下 随时间变化的函数，编入Fluent进行非稳态解算 假设： 实际中孔隙率是随着工作面的推进、周期来压和 (1)采空区02的消耗主要以煤氧化学反应为主， 煤岩特征不断变化的.采空区应力、时间及孔隙率的 C0为煤氧反应的主要产物，遗煤分布在底板上方0~ 关系近似符合负指数关系0，E(e,l,y）=aeo,, 2.5m的空间，上部为岩体，认为不发生反应，只进行 o1为顶板压力，a为最大孔隙率，v是推进速度，b是待 风流的输运和热量的传递，模型高度为10m 定系数.孔隙率主要受两个因素影响，一个是冒落矸 (2)每隔40m设置一个联络巷，由于动网格计算 石压力，由于此压力与基本顶压力相差较大，所以可以 量过大，所以针对联络巷相对工作面距离的周期变化， 忽略；另一个是岩体粒径，所以认为孔隙率在采空区纵 只进行一个周期的计算，对煤氧低温氧化的初级阶段 深方向上主要受工作面推进速度的影响.在对文 进行研究 献6]的孔隙率分布函数的基础上，本文在孔隙率函 (3)假设采空区内同一坐标位置气体温度和煤温 数的y轴方向添加源项，该源项是关于(y,,t)的函 相等，工作面推进时新垮落的煤岩温度为初温. 数，从而得到四维空间孔隙率分布函数： 计算模型如图2所示，煤与空气的物理参数见 ￡=8.，：= 表1
工程科学学报，第 38 卷，第 10 期 Qx x + Qy y + Qz z = 0， CO2 t + Qx ·CO2 x + Qy ·CO2 y + Qz ·CO2 z = De ( 2 CO2 x 2 + 2 CO2 y 2 + 2 CO2 z 2 ) － V( T) ， ρc cpc T t = Q( C，T) + λc ( 2 T x 2 + 2 T y 2 + 2 T z 2 ) － ρg cpg ( Qx ·T x + Qy ·T y + Qz ·T ) z ， ( x Kx H ) x + ( y Ky H ) y + ( z Kz H ) z = 0 ． ( 5) 该多元方程组中，x、y、z 为三维空间坐标，Q 为漏 风流量，K 为绝对渗透率，假设多孔介质为各向同性， 根据实验数据 Kx = Ky = Kz = 1 × 10 － 6 m2 ，CO2 为 O2 的 质量分数，取 24% ，De 为氧气在煤体中的扩散系数， V( T) 为实际耗氧速率，mol·s － 1·cm － 3 ; H 为热力分压， Pa; Q( C，T) 为煤氧化反应放热源项［17--20］． 综上所述，影响煤氧反应速率的因 素 主 要 有 温 度、氧含量及孔隙结构，化学反应中组分的变化及放 热量又是反应速率的函数． 建立 U + L 型通风方式 的采空区几何模型，将煤氧反应的机理编写成用户 自定义函数程序，导入模拟软件 Fluent 中，为数值计 算增加化学反应的组分源项和能量源项，研究采空 区的氧化升温规律，力求减小计算误差，提高仿真模 拟的真实性． 1. 2 U + L 通风采空区四维动态模型的建立 采空区工作面推进是一个动态过程，使用计算流 体力学( CFD) 动网格的方法，进行网格动态重构，建立 采空区四维动态模型，使采空区范围不断扩大． 使用 用户自定义函数中宏 ＲP_Get_Ｒeal“flow-time”获取工 作时间，并以此建立工作面推进速度、孔隙率和渗透率 随时间变化的函数，编入 Fluent 进行非稳态解算． 实际中孔隙率是随着工作面的推进、周期来压和 煤岩特征不断变化的． 采空区应力、时间及孔隙率的 关系近似符合负指数关系［21］，ε( v，t，y) = ae － bσ1( xi，yi ，t) ， σ1 为顶板压力，a 为最大孔隙率，v 是推进速度，b 是待 定系数． 孔隙率主要受两个因素影响，一个是冒落矸 石压力，由于此压力与基本顶压力相差较大，所以可以 忽略; 另一个是岩体粒径，所以认为孔隙率在采空区纵 深方向 上 主 要 受 工 作 面 推 进 速 度 的 影 响． 在 对 文 献［6］的孔隙率分布函数的基础上，本文在孔隙率函 数的 y 轴方向添加源项，该源项是关于( y，v，t) 的函 数，从而得到四维空间孔隙率分布函数: ε = εxεyεz = ( 0. 2e － 0. 022 3( y + vt) + 0. 1) ( e － 0. 15x + 1) × 1. 05z ， x≤L /2; ( 0. 2e － 0. 022 3( y + vt) + 0. 1) ( e － 0. 15( L － x) + 1) × 1. 05z { ， x ＞ L /2． ( 6) 式中: L 是工作面走向长度，m; 初始孔隙率空间分布函 数如图 1 所示，其中 εx、εy和 εz分别是沿工作面走向( x 轴) ，采空区走向( y 轴) 和采空区高度方向( z 轴) 上的 孔隙率，孔隙率分布不仅是空间坐标的函数，同时是时 间的函数，空间某位置的孔隙率随着工作面的推进不 断减小． 由 Carman--Kozeny 公式得出，渗透率 k 与孔隙 率 ε 的关系为 k = d0 0. 241 ε2 ( 1 － ε) 2 ． ( 7) 式中，d0是粒径的平方，量级可取10 － 9 m2 ． 将孔隙率代 入其中可以得到渗透率的空间分布及随时间的变化规 律，使用用户自定义函数中 DEFINE_PＲOFILE 将孔隙 率及渗透率编入 Fluent 进行解算． 图 1 初始孔隙率空间分布函数曲面 Fig． 1 Initial porosity function space distribution 本文以山西某矿 U + L 通风采空区为原型建立几 何模型． 工作面长度为 200 m，采空区走向取 250 m， 进、回风巷及联络巷宽 4 m，高 2. 5 m，长 10 m． 其中 4# 煤层 赋 存 于 山 西 组 下 部，煤 厚 0. 84 ～ 6. 05 m，平 均 2. 98 m，含夹矸 1 ～ 4 层． 4# 和 3# 煤层之间的距离为 9. 32 m，3# 和 2# 煤层之间的距离为 6. 94 m，属于近距离 煤层． 根据该矿的实际情况，结合模拟要求作 以 下 假设: ( 1) 采空区 O2 的消耗主要以煤氧化学反应为主， CO 为煤氧反应的主要产物，遗煤分布在底板上方 0 ～ 2. 5 m 的空间，上部为岩体，认为不发生反应，只进行 风流的输运和热量的传递，模型高度为 10 m． ( 2) 每隔 40 m 设置一个联络巷，由于动网格计算 量过大，所以针对联络巷相对工作面距离的周期变化， 只进行一个周期的计算，对煤氧低温氧化的初级阶段 进行研究． ( 3) 假设采空区内同一坐标位置气体温度和煤温 相等，工作面推进时新垮落的煤岩温度为初温． 计算模型如 图 2 所 示，煤与空气的物理参数见 表 1． · 2531 ·
周佩玲等：非均质孔隙率采空区氧化升温规律四维动态模拟 ·1353· 0s020160300 200 150 2计算结果及分析 100 50 0 10 2.1非均质孔隙率分布及四维变化 进风巷 在建立了采空区四维动态变化模型的基础上，对 岩层 煤层 回风巷 采空区模型进行解算，得到孔隙率时空演化过程如 联络巷 y 图3所示.图3是推进速度为2.4m·d时，1=2、9 图2计算模型示意图（单位：m) 和16d时孔隙率的时空演化过程.可见孔隙率在空间 Fig.2 Schematic diagram of the calculation model (unit:m) 以非均质形式分布，工作面附近及两巷位置由于支护 表1煤与空气的物理参数 及巷道的作用，孔隙率较大，中间压实区域孔隙率较 Table 1 Physical parameters of coal and air 小，符合“0”形圈理论 量名称 数值 单位 图4为采空区距底板1m高截面孔隙率动态变 煤的密度 1300 kg'm-3 化.随着开采的不断进行，采空区纵深增长，孔隙率的 煤的比热容 1003.2 J小kglK 同一等值线随着工作面向前移动，深部孔隙率不断减 煤的热传导率 0.1998 W.m-1.K-1 小，根据相对运动原理，相当于近工作面前端100m范 反应热(1mol02) 300 ◇ 围内孔隙率始终不变.如图5所示，采空区某一空间 空气的比热容 1.0 J.kg-1.K- 点(10,10,1)的孔隙率随时间以负指数形式不断减 空气的热传导率 2.6×10-2 Wml…Kl 小，这与实际采空区孔隙率不断被压实变化更接近 空气动力黏度 1.8×10-5 kg'm-1.s-1 2.2动态推进采空区氧化升温规律 2.2.1工作面推进速度对采空区温度场的影响 气体扩散系数 1.5×10-5 m2s-1 由于工作面的不断推进，采空区是一个动态地不 初始温度 300 K 断扩大地空间，所以影响煤氧反应的主要因素有遗煤 进口风速 1.4 厚度、孔隙率、工作面风量、工作面推进速度等.其中， 孔隙率 (b) ■0.76 0.63 0.52 0.38 d 0.25 0.13 0 图3孔隙率时空演化图.(a)t=2d,z=7m截面：(b)t=2d,空间图：(ct=9d,z=7m截面：(d)t=9d,空间图：(e)t=16d,z=7m 截面：(0t=16d,空间图 Fig.3 Porosity temporal and spatial evolution:(a)t=2d,:=7m,section:(b)t=2d,space:(c)t=9d,z=7m,section:(d)t=9d,space: (e)t=16d,z=7m,section:(f)t=16d,space 遗煤厚度、工作面风量等对采空区升温的影响已有前 风侧.工作面推进速度越大采空区升温速率越小，温 人做出部分研究，本文主要针对工作面推进速度对采 度越低，v=1.2m·d-时最高温度为313.39K,v= 空区升温规律的影响进行研究 3.6md-时最高温度为305.46K,主要是因为推进相 工作面推进速度是影响煤体氧化升温过程的主要 同的距离，如果推进速度大则采空区内物质输运速度 因素，决定了高温点的温度及位置.下面分别对推进 增加，工作面附近低温煤体温度增加速度加快，同时温 速度为1.2、2.4、3.6和6md-时煤氧低温氧化升温 度较高的煤体将迅速进入室息带，从而抑制了温度的 规律进行模拟计算。 升高，对预防遗煤的自燃具有显著作用.然而，工作面 图6是工作面推进18m时，距底板1m高度处采 推进速度越大，高温区域深度越大，如表2所示，高温 空区温度场等值线分布.可知高温区域主要集中在进 区域深度的增加又不利于对遗煤自燃的防治，所以对
周佩玲等: 非均质孔隙率采空区氧化升温规律四维动态模拟 图 2 计算模型示意图( 单位: m) Fig． 2 Schematic diagram of the calculation model ( unit: m) 表 1 煤与空气的物理参数 Table 1 Physical parameters of coal and air 量名称 数值 单位 煤的密度 1300 kg·m － 3 煤的比热容 1003. 2 J·kg － 1·K － 1 煤的热传导率 0. 1998 W·m － 1·K － 1 反应热( 1 mol O2 ) 300 kJ 空气的比热容 1. 0 J·kg － 1·K － 1 空气的热传导率 2. 6 × 10 － 2 W·m － 1·K － 1 空气动力黏度 1. 8 × 10 － 5 kg·m － 1·s － 1 气体扩散系数 1. 5 × 10 － 5 m2 ·s － 1 初始温度 300 K 进口风速 1. 4 m·s － 1 2 计算结果及分析 2. 1 非均质孔隙率分布及四维变化 在建立了采空区四维动态变化模型的基础上，对 采空区模型进行解算，得到孔隙率时空演化过程如 图 3 所示． 图 3 是推进速度 v 为 2. 4 m·d － 1 时，t = 2、9 和 16 d 时孔隙率的时空演化过程． 可见孔隙率在空间 以非均质形式分布，工作面附近及两巷位置由于支护 及巷道的作用，孔隙率较大，中间压实区域孔隙率较 小，符合“O”形圈理论． 图 4 为采空区距底板 1 m 高截面孔隙率动态变 化． 随着开采的不断进行，采空区纵深增长，孔隙率的 同一等值线随着工作面向前移动，深部孔隙率不断减 小，根据相对运动原理，相当于近工作面前端 100 m 范 围内孔隙率始终不变． 如图 5 所示，采空区某一空间 点( 10，10，1) 的孔隙率随时间以负指数形式不断减 小，这与实际采空区孔隙率不断被压实变化更接近． 2. 2 动态推进采空区氧化升温规律 2. 2. 1 工作面推进速度对采空区温度场的影响 由于工作面的不断推进，采空区是一个动态地不 断扩大地空间，所以影响煤氧反应的主要因素有遗煤 厚度、孔隙率、工作面风量、工作面推进速度等． 其中， 图 3 孔隙率时空演化图 . ( a) t = 2 d，z = 7 m 截面; ( b) t = 2 d，空间图; ( c) t = 9 d，z = 7 m 截面; ( d) t = 9 d，空间图; ( e) t = 16 d，z = 7 m 截面; ( f) t = 16 d，空间图 Fig． 3 Porosity temporal and spatial evolution: ( a) t = 2 d，z = 7 m，section; ( b) t = 2 d，space; ( c) t = 9 d，z = 7 m，section; ( d) t = 9 d，space; ( e) t = 16 d，z = 7 m，section; ( f) t = 16 d，space 遗煤厚度、工作面风量等对采空区升温的影响已有前 人做出部分研究，本文主要针对工作面推进速度对采 空区升温规律的影响进行研究． 工作面推进速度是影响煤体氧化升温过程的主要 因素，决定了高温点的温度及位置． 下面分别对推进 速度为 1. 2、2. 4、3. 6 和 6 m·d － 1时煤氧低温氧化升温 规律进行模拟计算． 图 6 是工作面推进 18 m 时，距底板 1 m 高度处采 空区温度场等值线分布． 可知高温区域主要集中在进 风侧． 工作面推进速度越大采空区升温速率越小，温 度越低，v = 1. 2 m·d － 1 时最 高 温 度 为 313. 39 K，v = 3. 6 m·d － 1时最高温度为 305. 46 K，主要是因为推进相 同的距离，如果推进速度大则采空区内物质输运速度 增加，工作面附近低温煤体温度增加速度加快，同时温 度较高的煤体将迅速进入窒息带，从而抑制了温度的 升高，对预防遗煤的自燃具有显著作用． 然而，工作面 推进速度越大，高温区域深度越大，如表 2 所示，高温 区域深度的增加又不利于对遗煤自燃的防治，所以对 · 3531 ·
·1354. 工程科学学报，第38卷，第10期 孔隙率 0.30m ■0.64 =2.4m… 0.58 0.51 △/=4.8m 0.28 ◆ 0.45 3d .3 0.26 8 6d .1 0.13 9d 0.06 进风口 0.22 12d 0.20 ■ 0.18 2 6 1012 联络巷、 时间/d 出风口 100 250278.8 图5空间某点孔隙率与时间的关系 采空区走向长度/m Fig.5 Relationship between porosity and time 图4工作面推进及孔隙率动态变化图 于推进速度的制定要综合考虑其对采空区升温速率的 Fig.4 Dynamic change of porosity when the working face is ad- 影响、高温区域的深度等因素 vanced b 最高温度为313.39K 最高温度为：306.67K 312K 306K 301K 302K 302K 303K 回 303K04K 0 100 200 巷 100 200 50 150 50 150 d 最高温度为：305.46K 最高温度为：305.83K 305.83K 305K 301K 301K 302K 30263K304K 303K304K 必 回 卷 100 200 100 巷 150 50 150 图6不同推进速度下采空区温度场分布.(a)r=1.2md:(b)=2.4md-1:(c)m=3.6md1:(d)m=2.4md1,U型通风 Fig.6 Temperature field distribution in the gob at different advancing speeds:(a)=1.2md-;(b)=2.4m'd-:(c)=3.6m'd-!;(d) =2.4m'd-1,U-ype ventilation 本文研究的是U+L型通风采空区，对比图6(b) 于瓦斯的抽放和防治，但同时会对遗煤自燃的防治产 和()，在同一推进速度条件下，尾巷的联络巷具有热 生影响 量的传输，近联络巷处温度较高，使得U+L通风采空 2.2.2动态推进采空区温度场四维动态演化规律 区温度场范围大于U型通风采空区，温度也增加，最 推进速度为2.4m·d时，在采空区建立监测点 高温度也高于U型采空区近1K,虽然设置尾巷有利 阵，对采空区温度场的时空变化过程进行研究
工程科学学报，第 38 卷，第 10 期 图 4 工作面推进及孔隙率动态变化图 Fig． 4 Dynamic change of porosity when the working face is advanced 图 5 空间某点孔隙率与时间的关系 Fig． 5 Ｒelationship between porosity and time 于推进速度的制定要综合考虑其对采空区升温速率的 影响、高温区域的深度等因素． 图 6 不同推进速度下采空区温度场分布． ( a) v = 1. 2 m·d － 1 ; ( b) v = 2. 4 m·d － 1 ; ( c) v = 3. 6 m·d － 1 ; ( d) v = 2. 4 m·d － 1，U 型通风 Fig． 6 Temperature field distribution in the gob at different advancing speeds: ( a) v = 1. 2 m·d － 1 ; ( b) v = 2. 4 m·d － 1 ; ( c) v = 3. 6 m·d － 1 ; ( d) v = 2. 4 m·d － 1，U-type ventilation 本文研究的是 U + L 型通风采空区，对比图 6( b) 和( d) ，在同一推进速度条件下，尾巷的联络巷具有热 量的传输，近联络巷处温度较高，使得 U + L 通风采空 区温度场范围大于 U 型通风采空区，温度也增加，最 高温度也高于 U 型采空区近 1 K，虽然设置尾巷有利 于瓦斯的抽放和防治，但同时会对遗煤自燃的防治产 生影响． 2. 2. 2 动态推进采空区温度场四维动态演化规律 推进速度为 2. 4 m·d － 1 时，在采空区建立监测点 阵，对采空区温度场的时空变化过程进行研究． · 4531 ·