一种基于变形几何的采空区漏风流场动态数值模拟方法
【专利摘要】本发明提供一种基于变形几何的采空区漏风流场动态数值模拟方法，包括：根据以每日平均推进速度开采一天后所形成的采空区，建立采煤工作面开采的初始采空区几何模型；设定采空区自由变形区域、固定边界、移动边界；设置采空区漏风流场；将初始采空区几何模型划分网格并开展采空区漏风场数值模拟计算；以每日平均推进速度逐步推进采空区几何模型并重新剖分网格并开展采空区漏风场数值模拟计算，直至推进至停采线或未至停采线但达到设定停止计算时间。本方法采用变形几何控制采煤工作面和采空区边界移动和区域变化的动态演化过程，采用自由和多孔介质流体流动控制工作面和采空区的流场，二者耦合计算描述采空区漏风流场动态变化和分布规律。
【专利说明】
一种基于变形几何的采空区漏风流场动态数值模拟方法
技术领域
[0001] 本发明涉及采空区漏风模拟技术领域，具体涉及一种基于变形几何的采空区漏风 流场动态数值模拟方法。
【背景技术】
[0002] 采空区漏风是影响工作面安全开采的重要因素，漏风导致采空区遗煤自燃并诱导 采空区有毒有害气体大量向工作面涌出，因此，必须明确采空区漏风规律为煤矿安全生产 提供理论依据和技术支持。常用的采场风流场数值解算与仿真方法较早的主要有20世纪80 ~90年代的基于FORTRAN语言的自编程序采场流场计算方法，21世纪以来的G3软件仿真方 法、双坐标系方法、移动坐标系方法，还有常用的商业软件PHOENICS、fluent、comsol等方 法，以上各种方法对采空区流场的仿真计算一定程度上提高了采场流场理论计算的水平为 工作面安全生产提供了指导，但这些方法部分不能实现工作面推进过程中采空区的动态变 化，部分基于自编程序开发而得技术和理论难度较大在本行业的推广和应用上具有一定局 限。
【发明内容】
[0003] 本发明的目的在于提供一种基于变形几何的采空区漏风流场动态数值模拟方法。
[0004] 本发明的技术方案是：
[0005] -种基于变形几何的采空区漏风流场动态数值模拟方法，包括：
[0006] 根据以每日平均推进速度开采一天后所形成的采空区，建立采煤工作面开采的初 始采空区几何模型：以每日平均推进速度为采空区长度增加速率，基于变形几何控制采空 区几何模型随采煤工作面的回采后推而延长，设置采煤工作面的进风巷和回风巷；
[0007] 设定采空区自由变形区域、采空区固定边界、采空区移动边界；
[0008] 设置采空区漏风流场:将采空区设置为多孔介质流动区域，将采煤工作面设置为 自由流体流动区域，将采空区岩石碎胀系数为随采煤工作面推进速度、推进时间、y坐标控 制的变量，将采空区冒落煤岩的孔隙率、渗透率、遗煤厚度为随采空区岩石碎胀系数变化的 变量;采煤工作面的进风巷为风速控制边界，采煤工作面的回风巷为风压控制边界，采用风 压差逐点沿程控制采空区外边界的分压分布；
[0009] 将初始采空区几何模型划分网格并开展采空区漏风场数值模拟计算，包括采空区 岩石碎胀系数计算、采空区冒落煤岩的孔隙率、渗透率、遗煤厚度、风压、漏风速度计算；
[0010] 以每日平均推进速度逐步推进采空区几何模型并重新剖分网格并开展采空区漏 风场数值模拟计算，直至推进至停采线或未至停采线但达到设定停止计算时间时完成采空 区漏风流场数值计算。
[0011] 将初始采空区几何模型设定为所述采空区自由变形区域;将开切眼设定为所述固 定边界，作为全局坐标系。
[0012] 所述采空区移动边界包括X方向移动、y方向固定的单向移动边界和y方向固定、X 方向由所述单向移动边界自由向前推进的自由移动边界。
[0013] 所述单向移动边界的X方向移动速度为采煤工作面回采速度，所述单向移动边界 的X方向位移采用网格位移量依据推进速度和推进时间乘积指定。
[0014] 所述单向移动边界包括采煤面以及采空区的上下帮。
[0015] 有益效果：
[0016] 本方法基于采煤工作面连续推进的采空区边界动态移动和采空区区域动态增长 的现场开采实际，采用comsol软件中的内置的变形几何方法控制采煤工作面和采空区边界 移动和区域变化的动态演化过程，采用自由和多孔介质流体流动控制工作面和采空区的流 场，二者耦合计算和仿真描述采场动态演化过程中的采空区漏风流场动态变化和分布规 律。通过本方法可以实现工作面推进过程中采空区漏风流场的动态描述，既符合工程现场 实际也能提高理论计算水平，能为采空区动态演化过程中漏风流场动态变化过程的预测和 描述提供更便利和准确的方法，还能为采空区有毒有害气体的涌出、自然发火等问题的预 测与深入研究和治理方法的提出提供基础支撑。
【附图说明】
[0017] 图1是本发明【具体实施方式】的初始采空区几何模型；
[0018] 图2是本发明【具体实施方式】的采煤工作面推进不同时间时采空区冒落煤岩碎胀系 数，（a)为采煤工作面推进20天时采空区冒落煤岩碎胀系数，（b)为采煤工作面推进40天时 采空区冒落煤岩碎胀系数，（c)为采煤工作面推进60天时采空区冒落煤岩碎胀系数，（d)为 采煤工作面推进90天时采空区冒落煤岩碎胀系数；
[0019] 图3是本发明【具体实施方式】的距开切眼30m沿倾向碎涨系数随开采时间延续的变 化规律；
[0020] 图4是本发明【具体实施方式】的开切眼y = 85m、45m沿程碎胀系数随开采时间延续的 变化规律；
[0021] 图5是本发明【具体实施方式】的采空区遗煤厚度随开采时间延续的变化规律，（a)采 煤工作面推进20天时采空区遗煤厚度，（b)采煤工作面推进40天时采空区遗煤厚度，（c)采 煤工作面推进60天时采空区遗煤厚度，（d)采煤工作面推进90天时采空区遗煤厚度；
[0022] 图6是本发明【具体实施方式】的距开切眼30m沿倾斜方向遗煤厚度随开采时间延续 的变化规律，（a)为沿程全局分布规律，（b)为两道附近的遗煤厚度分布规律，（c)为采空区 中部的遗煤厚度分布规律；
[0023]图7是本发明【具体实施方式】的开切眼处y = 85m、45m沿程遗煤厚度随开采时间延续 的变化规律；
[0024] 图8是本发明【具体实施方式】的开采推进速度3m/d生产条件下不同时期采空区冒落 煤岩的孔隙率分布规律，（a)为采煤工作面推进20天时采空区冒落煤岩的孔隙率分布规律， (b)为采煤工作面推进40天时采空区冒落煤岩的孔隙率分布规律，（c)为采煤工作面推进60 天时采空区冒落煤岩的孔隙率分布规律，（d)为采煤工作面推进90天时采空区冒落煤岩的 孔隙率分布规律；
[0025] 图9是本发明【具体实施方式】的距离开切眼30m处冒落煤岩的孔隙率随时间变化规 律；
[0026] 图10是本发明【具体实施方式】的距离开切眼10m，y = 45、85m处冒落煤岩的孔隙率变 化规律；
[0027] 图11是本发明【具体实施方式】的冒落煤岩的孔隙率衰减规律；
[0028] 图12是本发明【具体实施方式】的渗透率衰减规律；
[0029]图13是本发明【具体实施方式】的开切线7 = 5111，7 = 45111，7 = 85111处3个点随采煤工作 面推进和时间延长漏风速度变化规律与常规方法下采空区区域和边界不变化时随采煤工 作面推进和时间延长漏风速度变化规律;其中（a)、（c)、（e)分别描述点位置为开切线y = 5111，7 = 45111，7 = 85111处3个点随采煤工作面推进和时间延长漏风速度变化规律，（13)、（(1)、（〇 分别开切线7 = 5111，7 = 451]1，7 = 851]1处3个点在常规方法下采空区区域和边界不变化；
[0030] 图14是本发明【具体实施方式】的不同推进速度下采空区漏风规律，（a)为风速0.8m/ s下采空区漏风规律，（b)为风速lm/s下采空区漏风规律，（c)为风速1.4m/s下采空区漏风规 律，（d)为风速1.8m/s下采空区漏风规律。
【具体实施方式】
[0031] 下面结合附图对本发明的【具体实施方式】做详细说明。
[0032] -种基于变形几何的采空区漏风流场动态数值模拟方法，包括：
[0033] 步骤1、根据以每日平均推进速度开采一天后所形成的采空区，建立采煤工作面开 采的初始采空区几何模型：以每日平均推进速度为采空区长度增加速率，基于变形几何控 制采空区几何模型随采煤工作面的回采后推而延长，设置采煤工作面的进风巷和回风巷；
[0034] 本实施方式中初始采空区几何模型如图1所示，采煤工作面长90m，进回风巷道宽 3.5m、高4m，回采掌子宽约7m，即煤壁与采空区最大距离7m为例，如图1所示，LI = 3.5m，L2 = 7m。采用变形几何进行采空区控制，即以1天为时间步长，以每日平均推进速度为采空区长 度增加速率随采煤工作面的后退延长采空区，初始设置时第一天如图1中的域A所示的采空 区已形成，L3 = vt · l[d](vt为每日平均推进速度，单位是m/d)，后续开采过程中采空区长 度不断增加，其长度为L3 = vt · t[d]，式中d表示天。图中线3-3为始采的开切眼位置，模拟 中由于采煤工作面采用回退回采方法，所以设置该线处于采煤工作面移动几何坐标系横坐 标的零点，即X1 = 0m，依次类推。x2 = -L3 = -vt · 1 [d]，x3 = x2-7m，图中x3、x4之间的距离为 模拟过程中的进回风巷道长度，为了简化模型和节省计算资源本实施方式中将该段长度L4 设置为3m，则x4 = x3-3m。域C1、域C2分别为回风巷道和进风巷道部分，域B为采煤工作面空 间，其中省略采煤机、电缆槽、两排液压支架等设备，将煤壁、顶底板以及各种设备对通风的 影响考虑入风阻中即可相对降低由于省略以上设备所带来的误差。采煤工作面移动几何坐 标系中y轴的y〇 = 〇m，yl = 3.5m，y2 = 86.5m，y3 = 90m。
[0035] 几何变形的域选择与控制:变形几何总控制区域为域Cl、域C2、域B、域A;设定几何 变形节次为2即几何变形控制函数为二阶导函数，自由变形设定网格平滑控制方法采用 Wins low网格平滑方法。
[0036] 步骤2、设定采空区自由变形区域、采空区固定边界、采空区移动边界；
[0037] 将初始采空区几何模型设定为所述采空区自由变形区域;将开切眼设定为固定边 界，作为全局坐标系。
[0038] 采空区移动边界包括X方向移动、y方向固定的单向移动边界和y方向固定、X方向 由所述单向移动边界自由向前推进的自由移动边界。
[0039] 单向移动边界的X方向移动速度为采煤工作面回采速度，所述单向移动边界的X方 向位移采用网格位移量依据推进速度和推进时间乘积指定。单向移动边界包括采煤面以及 采空区的上下帮。
[0040] 本实施方式具体设定如下：
[0041]设定采空区自由变形区域：由于在采煤工作面回采后推过程中这个采煤工作面和 采空区都是向前移动的，因此将域C1、C2,域B，域A均设置为自由变形区域。自由变形区域设 定网格初始位移dx〇 = 0m，dyo = 0m。
[0042] 设定采空区固定边界：图1中的开切眼设置3-3设定为固定边界，该边界不因采煤 工作面和采空区的移动而变化，该固定边界的坐标系为全局坐标系，即该边界的材料和几 何坐标系都是固定不移动和变形的，指定其x、y方向网格位移一直为〇m。
[0043] 设走米空区移动边界一移动边界1和移动边界2:
[0044] (1)设定移动边界1:随采煤工作面后退和采空区长度逐步增加，取时间步长为1 天，采煤工作面后推过程中边界2-1、3_1为X方向移动、y方向固定的单向移动边界;X方向移 动速度为采煤工作面回采速度。
[0045] (2)设定移动边界2:采煤工作面推进和采空区延长的过程中还存在另一类y方向 固定不移动、X方向由上述移动边界1自由向前推进的移动边界，该类边界主要集中在采煤 工作面以及采空区的上下帮，如图1中的1-2、1-3、1-4、1-5、1-6、2-2、2-3、3-2、3-4,这类边 界的X方向随工作面移动不设置速度，y方向不移动。
[0046] 步骤3、设置采空区漏风流场:将采空区设置为多孔介质流动区域，将采煤工作面 设置为自由流体流动区域，将采空区岩石碎胀系数为随采煤工作面推进速度、推进时间、y 坐标控制的变量，将采空区冒落煤岩的孔隙率、渗透率、遗煤厚度为随采空区岩石碎胀系数 变化的变量;采煤工作面的进风巷为风速控制边界，采煤工作面的回风巷为风压控制边界， 采用风压差逐点沿程控制采空区外边界的分压分布；
[0047] 设置自由和多孔介质流体流动区域：自由流体流动区域为域C1、域C2、域B，多孔介 质流体流动区域为域A;
[0048] (1)自由流体流动区域
[0049]自由流体流动的控制域为域C1、域C2、域B，数值计算采用方程（1)~（2)控制该区 域的流体自由流动：
[0052]式中：F为自由流体流动空间中由采煤工作面的风流温度、湿度、密度、采煤工作面 上下道高差、重力等引起的上浮力，单位是N/m3;P为热气流密度，单位是kg/m3;p2为采煤工 作面绝对风压，即自由湍流的气压，Pa;u为气流的矢量流速，m/s;u为气体的动力粘度，Pa/ s2; ▽为流速梯度。
[0053] (2)多孔介质流体流动区域
[0054]多孔介质流体流动的控制域为域A，其内部流体流动的方程采用（3)~(4)进行控 制：
[0057] 式中：ερ为冒落煤岩的孔隙率;κ为多孔介质的渗透率;Qbr为多孔介质内气流质量 流量，单位是kg/m 3/s。
[0058] 采空区冒落煤岩的孔隙率εΡ计算方法如式(1)，渗透率κ计算方法如式(2)。
[0059] ερ= l-1/kp (1)
[0060] κ = ερ3?ρ2/150(1-ερ)2 (2)
[0061] 式中:kp为碎胀系数，其计算方法如式(3) ;dp为采空区遗煤的平均粒径，单位是m。
[0062] kp = kp〇+(kpmax_kp〇)exp(ai(vt · t+L3+x)(l_exp(a2 · a3 · min(y，L_y)))) (3)
[0063] 式中，kpo为初始碎胀系数;kpmax为最大碎胀系数;ai、a2、a3分别为依据现场情况确 定的常数项;vt为采煤工作面每日平均推进速度，单位是m/d。
[0064] 采空区中遗煤/浮煤厚度hi采用(3)控制。
[0067]式中：hi为遗煤/浮煤厚度，m;ml为可采厚度，本实施方式根据瑞安公司014ΝΓ1综 放工作面实际取其为10.4m;切为碎胀系数，计算方法如式(3)，Π 为回采率，本实施方式根 据瑞安公司014ΝΓ1综放工作面实际按照式(4)取数。
[0068] 无渗漏壁（即无漏风边界）为边界1-3、2-2、3-2、3-3、3 -4、2-3、1-4、1-6、2-1、1 -5。
[0069] 设置风速控制边界与风压控制边界:将边界1-2设置为工作面进风巷入口，该入口 采用法向流入速度控制，即采煤工作面的进风巷为风速控制边界;将边界1-1设置为工作面 回风出口，采用出口压力控制，即采煤工作面的回风巷为风压控制边界;将采煤工作面和采 空区的共同边界3-1采用风压差逐点控制边界风压。
[0070] 步骤4、将初始采空区几何模型划分网格并开展采空区漏风场数值模拟计算，包括 采空区岩石碎胀系数计算、采空区冒落煤岩的孔隙率、渗透率、遗煤厚度、风压、漏风速度计 算；
[0071] 划分网格分采用普通物理场进行控制，网格尺寸采用特别细化处理，最大单元尺 寸0.5m，最小单元尺寸0.02m，最大单元生成率1.13,曲折因子0.3。
[0072] 步骤5、以每日平均推进速度逐步推进采空区几何模型并重新剖分网格并开展采 空区漏风场数值模拟计算，直至推进至停采线或未至停采线但达到设定停止计算时间时完 成采空区漏风流场数值计算。
[0073] 通过本实施方式的基于变形几何的采空区漏风流场数值模拟方法，可以分析出采 空区材料参数及流场变化规律：
[0074] (1)采空区碎胀系数及遗煤厚度变化规律
[0075] (1.1)碎胀系数
[0076]在不考虑开切眼附近只采不放煤区域时，图2(a)~（d)分别为采煤工作面推进20 天、40天、60天时、90天时采空区冒落煤岩碎胀系数，采煤工作面附近及上下两道附近冒落 煤岩的碎胀系数较大，开切眼附近与工作面附近类似，碎胀系数大冒落煤岩的孔隙率必然 大，因此综采放顶煤开采工艺条件下"两道"、"两线"为易自然发火区域。
[0077]图2中采空区深部和中部的碎胀系数较低，且随着开采时间的持续，采空区长度逐 渐加大，某一区域的碎胀系数也将发生变化，如图2~图3。
[0078]图3为距开切眼30m沿倾向碎胀系数随开采时间延续的变化规律，在采煤工作面倾 斜方向上碎胀系数两道大中部低关于中线对称(近水平采煤工作面可视为对称）；该区域与 开切眼相对位置保持不变，采煤工作面开采时间的持续，该区域逐步远离采煤工作面，碎胀 系数逐渐降低，且在开采初期碎胀系数较低较快，开采持续时间较长时如70~90天时，采空 区中部区域的碎胀系数基本一致，此段时间该区域已经埋入采空区190~250m，根据矿压现 象规律该区域也已基本稳定压实。综上，采煤工作面向后推进，采空区长度不断延长的过程 中，采空区的每一个区域都经历着碎胀系数由大到小逐渐衰减最终在该区域被充分压实之 后碎胀系数维持在最小过程，所以进行采煤工作面推进条件下的采空区数值模拟不能忽略 这一原则性关键规律。
[0079]图4为描述开切眼处，y = 45m和y = 85m两点碎胀系数随时时间的变化规律，图中D 表示采用本发明方法，C表示不采用变形几何方法的常规方法，后续图中出现具有同样意 义，图中显示采用本发明方法的条件下随时间的持续碎胀系数由大到小逐步衰减，采空区 中部区域衰减快于两道，开采时间较长之后此点的碎胀系数稳定到充分压实状态下的最小 值，y = 85m、45m处最小碎胀系数分别为1.081U.05;而示不采用变形几何方法的常规方法 的条件下，采空区的尺寸是预先设置固定的，所以任一点的碎胀系数只于其坐标有关，而常 规方法中材料坐标系不能移动，从而导致碎胀系数不随时间发生变化，与时间无关，碎胀系 数随开采时间延续分别一直保持为1.081、1.05，因此在采煤工作面推进过程中采空区相关 问题时如果不考虑动态过程，便不能实现采空区材料框架和几何框架的共同移动将不能如 实反应各物理场真实的变化规律。
[0080] (1.2)遗煤厚度变化规律
[0081] 图5(a)~(d)为采煤工作面推进不同时间后采空区遗煤厚度分布规律，图中显示 两道遗煤较厚，且靠近两道壁帮的区域遗煤厚度由于碎胀的影响将高于可采煤层厚度；同 时采空区遗煤分布也是非连续的，两道影响区域以外的遗煤厚度和影响范围内的遗煤厚度 差距较大。采空区遗煤厚度时间变化规律如图6~图7。
[0082]图6为采空区距离开切眼30m沿倾斜方向的遗煤分布规律，（a)为沿程全局分布规 律，（b)、（c)分别为两道附近和采空区中部的遗煤厚度分布规律，图中显示两道附近遗煤较 厚最大为10.92m大于10.4m的可采厚度，这是由于两道附近只采不放回采率底且碎胀系数 较大的缘故。如图6(a)所示，综放工作面采空区遗煤沿倾斜方向分布是非均匀的，尤其在两 道附件开采进行90d时期遗煤厚度在8-10.92m范围内，而中部区域在1.1~1.15m范围内，且 由于采空区中部压实程度较好，采空区中部遗煤厚度最小。由于顶板压力的原因，随开采时 间延续采空区各处遗煤厚度也逐渐降低，如图6(b)、（c)所示。
[0083]图7描述了在开切眼处取y = 45、85m处两点的浮煤厚度变化规律，图中显示切实考 虑采空区的动态形成过程时，该两点的遗煤厚度随时间的延续逐渐降低，其主要原因为随 时间延续在顶板压力作用下采空区冒落煤岩反复压实，碎胀系数变小导致遗煤厚度较低， 90天时y = 85、45m处的最小遗煤厚度分别为8.664、1.092m;两道的遗煤厚度随时间降低的 速率低于采空区中部;结合图6(c)，开采持续时间较长以后顶板逐步稳定压力也趋于稳定， 遗煤厚度也逐渐稳定且随时间延长而降低的趋势基本消失。图7中不考虑采空区的动态变 化而在数值模拟建模过程中直接将采空区的形状固定条件下的采空区遗煤厚度是不随时 间发生变化，y = 85、45m处的遗煤厚度分别保持为8.664、1.092m，原因与碎胀系数不发生变 化的原因类似；即不考虑采空区动态增长过程而在模拟时直接将一定推进速度下一定时间 内形成的采空区区域直接固定边界和区域大小，得到的采空区遗煤厚度只是上述设定开采 时间的最后时刻的遗煤厚度，这种数值计算方法忽略遗煤厚度在工作面推进和采空区增长 过程中的变化，属于在采空区瞬间形成的假设条件下开展的，是不科学的。
[0084] (2)采空区冒落煤岩的孔隙率、渗透率变化规律
[0085]图8(a)~(d)分别为开采推进速度3m/d生产条件下，采煤工作面推进20、40、60、90 天时采空区冒落煤岩的孔隙率分布，图中显示，采空区冒落煤岩的孔隙率在两道较高采空 区中部较低;对比不同开采时期下各云图显示任意区域在工作面推进一定时间后埋入采空 区深部，该区域的冒落煤岩的孔隙率随时间延长而变小。采空区边界移动、区域增长过程的 冒落煤岩的孔隙率变化详细规律如图9~图12所示。
[0086]图9为工作面推进10、30、50、70、90天时距离开切眼30m沿采煤工作面倾斜方向冒 落煤岩的孔隙率沿程分布及变化情况，图中显示采空区冒落煤岩的孔隙率深部较小两道较 大，且随开采时间延续该线所在区域冒落煤岩的孔隙率降低;采煤工作面推进30天以后冒 落煤岩的孔隙率随时间降低的趋势变弱，说明采空区中任意区域的冒落煤岩的孔隙率随采 煤工作面的推进及采空区区域的增长而降低，埋入采空区一定深度后趋于稳定，如图10所 不。
[0087] 图10更明确地对采空区冒落煤岩的孔隙率的变化规律进行描述，为距离开切眼 10m，y = 45、85m处冒落煤岩的孔隙率变化规律。图中显示，本发明的数值模拟下采空区任意 点（以距离开切眼l〇m，y = 45、85m处为例）的冒落煤岩的孔隙率随时间是衰减的，采空区中 间冒落煤岩的孔隙率衰减速率快于两道，在本图例选用的推进速度3m/d条件下，采空区中 心区域在开采60d后冒落煤岩的孔隙率开始稳定在0.05左右，上下道附近在工作面推进90d 时仍具有较大孔隙率并仍在持续衰减。而图中常规方法(采空区边界不移动、区域不增长） 所得冒落煤岩的孔隙率只是本方法的最后时刻稳定状态下的冒落煤岩的孔隙率，是不随工 作面推进而变化的，分别为〇.〇48、0.16，因为采用常规传统方法时首先建立3111/(1*90(1 = 270m长的采空区，任意点的坐标都是固定的，距离采煤工作面越远其冒落煤岩的孔隙率越 小，但与采煤工作面的距离不能随时间发生变化，而本方法中，采空区从无到有，由短到长 逐步按照推进速度增长，且采空区的工作面边界也依据推进速度和时间逐步移动。
[0088] 常规方法下某一点的冒落煤岩的孔隙率是一个固定常数，仅与其所处采空区的位 置有关，与时间无关，动态形成过程中该孔隙率随时间发生变化，因随工作面的推进采空区 形成并延长的过程中深部区域逐渐压实因此其孔隙率逐渐降低。按照"〇"形圈理论，冒落煤 岩的孔隙率与采空区某点相对位置有关(x，y)，而对于采空区中已经形成和存在任意点与 开切眼的距离是固定的，但与工作面的距离是逐渐增加的，必然要经历逐步和反复压实的 过程，待其进入采空区较深区域时顶板压力和压实程度才能基本稳定，因此，在不同开采速 度下，相同的时间内与采煤工作面的距离不同，推进速度越慢与工作面越近且增长率相对 较低，所以埋入采空区较深区域需要的时间相对较长，而推进速度大时，任意点或区域与采 煤工作面之间的距离快速增加，埋入采空区的速度较快所需时间较短，因此推进速度较大 时各项材料参数衰减较快，如图11~图12分别为距切眼倾向中点冒落煤岩的孔隙率和渗透 率衰减规律。
[0089] 根据图11所示的冒落煤岩的孔隙率衰减规律，采煤工作面推进90天时推进速度 2m/d条件下开切眼倾向中点冒落煤岩的孔隙率降低为0.04843,推进速度3~12m/d情况下 90天时冒落煤岩的孔隙率已经降低为0.04766，渗透率降低为3.5077e_ 8m2，且在90天之前某 一时间已经稳定在该值。中位粒径一定的条件下采空区冒落煤岩的孔隙率受碎胀系数控 制，渗透率受冒落煤岩的孔隙率控制，因此在碎胀系数稳定的前提下采空区冒落煤岩的孔 隙率和渗透率也应当是稳定的，推进速度3、5、6、8、10、12m/d时，孔隙率和渗透率达到稳定 所需时间为 89、53、45、34、27、23d。
[0090] 常规方法将采空区边界和区域范围进行固定，虽然采用迎风算法、移动坐标系法 等方法，但只能解决气氛浓度场、化学反应场、渗流场、气氛传递与扩散场、热场等场在时间 上的先后问题，不能解决或忽视了如碎胀系数、冒落煤岩的孔隙率、渗透率等关键材料参数 的瞬态变化问题，由于这些基本的关键材料参数不能实现瞬态变化必定导致数值模拟结果 的错误。
[0091] (3)采空区漏风规律
[0092] 常规方法数值模拟过程中采空区的风压风分布类似稳态分布，其不与时间构成显 著的相关关系，动态方法每一个计算时间步长采空区的长度都增加，所以每个时间步长中 采空区的冒落煤岩的孔隙率、渗透率等物理参数重新分布，其风压等也重新分布，由于风压 的重新分布导致漏风的重新分布，所以动态方法与常规方法的漏风分布是不同的，这也是 更符合采空区实际的。
[0093] 以工作面推进速度3m/d为例，工作面开采90d时采空区长度为270m，这270m是90d 内逐步形成的，因此采空区内漏风速度也将随这个过程发生变化。图13(a)、（c)、（e)分别描 述点位置为开切线7 = 5111，7 = 45111，7 = 85111处3个点随采煤工作面推进和时间延长漏风速度 变化规律，（b)、（d)、（f)分别开切线7 = 5111，7 = 451]1，7 = 851]1处3个点在常规方法下采空区区 域和边界不变化，采空区长270m不变。图13(a)、（c)、（e)显示依据采空区形成规律进行设置 的数值模拟结果反映出符合理论和实践的结果，即随采煤工作面推进采空区的逐步增长， 采空区任意点的漏风速度逐渐降低，当该点埋入采空区更深区域时漏风速度稳定在一个较 小值。
[0094]图13(b)、（d)、（f)为常规方法所得结果，图中显示，在不考虑采空区动态变化的客 观实际时，采空区任意点的漏风风速不随时间变化，与该方法下的各材料参数类似都置于 该点在采空区中的坐标有关，同时图13(b)、（d)、（f)的风速值与左列漏风速度稳定时期的 风速值相同，充分说明不考虑采空区动态过程的数值模拟等还不是标准的瞬态问题，只能 描述描述稳态结果。
[0095]图14(a)~（d)为风速0.8、1、1.4、1.8!11/8时不同推进速度下开切眼中点漏风速度 变化趋势，图中显示随开采时间延续漏风速度由大衰减并稳定在一个较低值，且推进速度 越大漏风速度衰减越快，即工作面推进速度越大采空区漏风达到稳态所需的时间越短。开 采时间相同和风速一定时，采空区漏风速度随推进速度增加而降低，这也是推进速度加快 有利于控制漏风，有利于控制采空区自然发火和有毒有害气体向工作面涌出的重要原因。
[0096] 采煤工作面在开采期间不断向后推进，采空区从无到有并逐步延长，即采空区的 动态变化是边界逐步移动、区域范围渐增的过程，因此对采空区相关问题进行分析不能将 采空区的边界和区域范围现行固定，必须以实际开采速度和时间为依据客观描述采空区的 逐步变化过程。
[0097] 本实施方式在建立的采空区动态演化数值模拟基础上，对采空区动态演化过程中 的材料参数如冒落煤岩碎胀系数、遗煤厚度、冒落煤岩的孔隙率、渗透率等以及采空区风压 和漏风分布及变化规律进行分析，得出以下主要结论：
[0098] (1)采空区动态演化过程中材料参数变化规律
[0099] 采空区动态演化和扩大过程中，碎胀系数、遗煤厚度、冒落煤岩的孔隙率、渗透率 等关键材料参数由大到小随推进时间衰减，埋入采空区较深区域后以上材料参数才能稳定 到较低值，达到稳态状态;采煤工作面推进速度越大，材料参数衰减越快;动态演化模型能 够全面呈现动态变化过程直至达到稳态。不考虑采空区动态演化过程的数值模拟描述的仅 是稳态结果，不能如实描述达到稳态之前的动态变化过程。
[0100] (2)采空区动态演化过程中漏风变化规律
[0101]漏风速度随采煤工作面推进动态变化，随时间延续逐渐衰减并至稳态，采煤工作 面开采时间相同时，采煤工作面推进速度越大漏风速度越小且衰减更快，有利于采空区防 火和控制瓦斯等有毒有害气体向采煤工作面涌出。
[0102] (3)基于采空区瞬间假设开展的数值模拟由于采空区边界和区域被先行固定不能 实现真正的采空区动态数值模拟，从而导致基于以上假设所开展的模拟结果在小推进速度 下误差大，在大推进速度下误差相对较低；同时指出基于采空区瞬态形成假设的方法仅适 用于采煤工作面长期停采(停采时间必须远大于模拟区域形成所需时间）条件下由压实区 外端至工作面这一相对较小范围内的相关问题，如停采撤架期间的自然发火、气体涌出、漏 风规律等。
[0103] 以上结果表明，采场动态演化数值模拟能如实描述采空区材料参数(碎胀系数、遗 煤厚度、冒落煤岩的孔隙率、渗透率)及风压和漏风速度的动态变化过程，而一般现行将采 空区边界及区域进行固定的数值模拟方法都是基于采空区瞬间形成，从而不能有效真实的 描述上述问题，只能呈现稳态结果，不是真实意义的采空区动态数值模拟方法。
[0104] 采场动态演化数值模拟方法及模型能真实呈现采空区边界移动、区域增长过程中 各变量及参数的变化规律。采空区材料如碎胀系数、遗煤厚度、冒落煤岩的孔隙率、渗透率 等在采空区动态演化过程中随时间逐步衰减后达稳定状态，采煤工作面推进速度越大衰减 越快;材料参数的衰减与变化导致采空区漏风速度随时间动态衰减，推进速度越大采空区 漏风衰减越快，漏风量越小，越利于采空区防火和控制有毒有害气体涌出。以上采空区材料 参数、漏风速度的动态变化规律表明对采空区进行数值模拟不能基于采空区瞬间形成的假 设，不能将采空区边界和区域先行固定，必须实现边界的动态移动和区域动态增长。
【主权项】
1. 一种基于变形几何的采空区漏风流场动态数值模拟方法，其特征在于，包括： 根据以每日平均推进速度开采一天后所形成的采空区，建立采煤工作面开采的初始采 空区几何模型：以每日平均推进速度为采空区长度增加速率，基于变形几何控制采空区几 何模型随采煤工作面的回采后推而延长，设置采煤工作面的进风巷和回风巷； 设定采空区自由变形区域、采空区固定边界、采空区移动边界； 设置采空区漏风流场:将采空区设置为多孔介质流动区域，将采煤工作面设置为自由 流体流动区域，将采空区岩石碎胀系数为随采煤工作面推进速度、推进时间、y坐标控制的 变量，将采空区冒落煤岩的孔隙率、渗透率、遗煤厚度为随采空区岩石碎胀系数变化的变 量;采煤工作面的进风巷为风速控制边界，采煤工作面的回风巷为风压控制边界，采用风压 差逐点沿程控制采空区外边界的分压分布； 将初始采空区几何模型划分网格并开展采空区漏风场数值模拟计算，包括采空区岩石 碎胀系数计算、采空区冒落煤岩的孔隙率、渗透率、遗煤厚度、风压、漏风速度计算； 以每日平均推进速度逐步推进采空区几何模型并重新剖分网格并开展采空区漏风场 数值模拟计算，直至推进至停采线或未至停采线但达到设定停止计算时间时完成采空区 漏风流场数值计算。2. 根据权利要求1所述的基于变形几何的采空区漏风流场动态数值模拟方法，其特征 在于，将初始采空区几何模型设定为所述采空区自由变形区域;将开切眼设定为所述固定 边界，作为全局坐标系。3. 根据权利要求1所述的基于变形几何的采空区漏风流场动态数值模拟方法，其特征 在于，所述采空区移动边界包括X方向移动、y方向固定的单向移动边界和y方向固定、X方向 由所述单向移动边界自由向前推进的自由移动边界。4. 根据权利要求3所述的基于变形几何的采空区漏风流场动态数值模拟方法，其特征 在于，所述单向移动边界的X方向移动速度为采煤工作面回采速度，所述单向移动边界的X 方向位移采用网格位移量依据推进速度和推进时间乘积指定。5. 根据权利要求3或4所述的基于变形几何的采空区漏风流场动态数值模拟方法，其特 征在于，所述单向移动边界包括采煤面以及采空区的上下帮。
【文档编号】G06F17/50GK105868472SQ201610187845
【公开日】2016年8月17日
【申请日】2016年3月29日
【发明人】齐庆杰, 董子文, 王欢, 于文惠, 李兴华, 秦凯
【申请人】辽宁工程技术大学