基于子域解析元素法的煤矿疏降水量预测研究

2021-07-27任晓波刘守强吴瑞芳

煤炭工程 2021年7期

任晓波，刘守强，吴瑞芳

(1.中国矿业大学(北京) 地球科学与测绘工程学院，北京 100083； 2.国家煤矿水害防治工程技术研究中心，北京 100083；3.河北省煤田地质勘查院，河北邢台 054000)

目前煤矿地下水流场的数值模拟主要采用有限差分、有限元等方法，本文旨在探索一种适用于煤矿地下水流场模拟的、效率更高且可为下一步可能需要进行的精细模拟提供指导(比如指明需要重点研究的区域以及收集哪些含水层的现场试验数据等)的模拟方法。本文首先构建基于子域解析元素法的煤矿地下水流场模型，然后通过煤矿疏降水量预测实例分析其模拟煤矿地下水流场的可行性及精度。子域解析元素法[1-3]是近年来在解析元素法[4-7]的基础上发展起来的一种比较适合处理局部地下水流问题的模拟方法。解析元素法的理论基础是通过叠加满足地下水流场控制方程和边界条件的线性解来实现对地下水流场的模拟。地下水动力学中用于模拟抽水井附近的定水头边界或隔水边界的镜像法便是比较简单的通过线性解的叠加来获得地下水流场的解析解，即通过在流场中添加各种水力要素(在镜像法中为虚拟井)以使得流场的边界条件得以满足，只不过在解析元素法中这种线性解的叠加方式更加复杂，叠加的线性解的数目更加庞大以至于难以通过人工计算来获得模拟结果。

1 子域解析元素法

1.1 线汇

经典的地下水流场分析中一般以水位h关于位置及时间的函数及偏导数来构成控制方程和边界条件表达式，在点(x，y)处，对于均质各向同性平面二维承压稳定流有：

式中，Qx、Qy为单宽流量分量，m2/d；qx、qy为达西渗流速度分量，m/d；H为含水层的厚度，m；K为含水层的渗透系数，m/d；h为含水层的水位，m。

式(1)中最后一项中的KHh为解析元素法中承压水流的流量势Φ，因此有：

在解析元素法中，流量势为描述地下水流场的基本单位之一，同时为了减少模拟计算中的变量数并充分利用复变函数中解析函数的特有性质，解析元素法对流场的模拟求解主要在复平面中进行。在复平面z=(x，y)上，式(2)中的流量势Φ(z)与其共轭调和函数Ψ(z)可组成解析函数Ω(z)。

Ω(z)=Φ(z)+iΨ(z)

(3)

由解析函数的性质可知，Φ(z)与Ψ(z)在流场内处处正交，因此可以分别用Φ(z)与Ψ(z)表示稳定流中的等水位线与流线，Ψ(z)即为流函数，Ω(z)为复势，二维承压稳定井流的复势可表示为：

式中，Q为井流量，m3/h；zw为井心坐标；C为常数。

用于模拟流场边界的一阶线汇复势表达式为[4，5]：

式中，z1、z2分别为线汇的首尾端点；Z=(X，Y)为基于线汇中心进行坐标变换后的局部坐标系；σ为线汇的强度，m2/d。

式(5)中线汇的强度σ为常数，当σ为关于位置的非线性三次函数时，线汇的复势表达式为[3]：

σ表达式的次数加1为其对应线汇的阶数，因此式(6)表示的是四阶非线性强度线汇。式(6)实部所代表的线汇流量势在跨越线汇时流量势保持连续，同时，式(6)虚部所代表的线汇流函数在跨越线汇时流函数值发生改变，这与式(5)所代表的常数强度线汇具有相同的流量势与流函数变化特征，因此高阶线汇也可用于流场边界的模拟，且较常数强度线汇在模拟边界时具有更高的灵活性和精度。

由式(2)、式(3)可知，由线汇控制的地下水流场中的水位可表示为Re(Ω)/(KH)，其中，Re(Ω)为线汇复势的实部。通过流场中两点z1、z2间(直线段z1z2不与线汇及线汇所在的负半轴相交)的流量可表示为|Im(Ω(z1))-Im(Ω(z2))|，其中，Im(Ω(z1))、Im(Ω(z2))分别为点z1、z2处线汇复势的虚部。

1.2 子域与区汇

类似于镜像法，解析元素法模型没有外边界，子域解析元素法采用解析元素法的理论基础并设置流场的外边界，在流场内部依据含水层参数(如渗透系数等)不同将其划分为若干个子域，流场的外边界和子域间的内边界均用线汇进行模拟。对于非稳定流，可将∂h/∂t有限差分为Δh/Δt[1]，二维承压平面非稳定流控制方程可表示为：

式中，γ为单位面积的承压含水层单位时间内弹性储(释)水量，m/d；S为承压含水层的储水系数，无量纲。

子域解析元素法采用区汇(area-sink)模拟由于承压含水层弹性储(释)水造成的水位变化，区汇的流量势表示如下：

区汇用径向基插值函数表示，因此有：

式中，ri为待求流量势的某点与区汇中第i(i=1，2，…n)个控制点的距离，m；r0为待求流量势的某点与区汇几何中心的距离，m；B0，bi(i=1，2，…n)分别为与区汇几何中心及第i个控制点相关的n+1个未知参数。

2 模型构建及求解

2.1 模型方程组

现以一假想研究区为例来说明如何生成模型方程组，如图1所示，示例中子域内外边界均采用二阶线汇模拟，也即每个线汇的复势表达式中含有两个关于线汇强度σ的未知参数。

图1 假想流场示意图

研究区依据含水层渗透系数不同划分为两个子域(子域1与子域2)，子域1外边界为隔水边界AB、定水头边界BC，子域1的内边界为CA，子域2外边界为隔水边界CD、法向定流量边界DA，子域2的内边界为AC，子域1的内边界CA与子域2的内边界AC为两条重合但相互独立的边界，这两条内边界分别用不同的线汇来模拟，为简便起见，子域的每条边界只用一个二阶线汇来表示，每个线汇上各设置两个控制点，然后每个线汇依照其两端端点和控制点可以划分两条法向流量控制分段。子域1与子域2中由抽水井W1产生的非稳定井流用区汇(式(7)—式(9))来模拟，每个子域对应一个不同的区汇，每个区汇设置3个控制点，同时，区汇的时间步长(式(7)中的Δt)设置为Δt1=t1-t0，Δt2=t2-t1，(t0，t1)与(t1，t2)为模拟过程经历的两个时间段。

每个子域内部及其内外边界上任意一点的水位与渗流速度由该子域的边界线汇、区汇共同决定而与子域外的任何线汇或区汇均无关，同时，为了保证在由各子域组成的研究区范围内满足达西渗流场能量与质量守恒，在子域间的内边界两侧需同时满足水位与法向流量处处相等，因此对于子域1与子域2的公共内边界AC上的控制点1、2有：

[Re(Ωt，1，AB，k+Ωt，1，BC，k+Ωt，1，CA，k)+Φt，1，k]/

(K1H1)+ht，1，k=

[Re(Ωt，2，CD，k+Ωt，2，DA，k+Ωt，2，AC，k)+Φt，2，k]/

(K2H2)+ht，2，k

t=Δt1，Δt2，k=1，2

(10)

式中，Ω为各子域的各个线汇在内边界控制点上产生的复势，Ω的四个下标分别代表线汇所处的时间段、线汇所属子域、线汇代表的边界、控制点；Φ为各子域的区汇在内边界控制点上产生的流量势，Φ的三个下标分别代表区汇所处的时间段、区汇所属子域、控制点；h为内边界控制点各时段初始水位，h的三个下标分别代表时间段、所属子域、控制点；K1、K2分别为子域1与子域2的渗透系数，m/d；H1、H2分别为子域1与子域2含水层厚度，m。

对于子域1与子域2的公共内边界AC上的法向流量控制分段Cq、qA有：

Im(Ωt，1，AB，C-Ωt，1，AB，q)+Im(Ωt，1，BC，C-

Ωt，1，BC，q)+Im(Ωt，1，CA，C-Ωt，1，CA，q)+

Ψt，1，Cq=Im(Ωt，2，CD，C-Ωt，2，CD，q)+

Im(Ωt，2，DA，C-Ωt，2，DA，q)+Im(Ωt，2，AC，C-

Ωt，2，AC，q)+Ψt，2，Cq

Im(Ωt，1，AB，q-Ωt，1，AB，A)+

Im(Ωt，1，BC，q-Ωt，1，BC，A)+

Im(Ωt，1，CA，q-Ωt，1，CA，A)+Ψt，1，qA=

Im(Ωt，2，CD，q-Ωt，2，CD，A)+Im(Ωt，2，DA，q-

Ωt，2，DA，A)+Im(Ωt，2，AC，q-Ωt，2，AC，A)+Ψt，2，qA

t=Δt1，Δt2

(11)

式中，Ψ为各子域的区汇在各个法向流量控制分段上产生的流量，Ψ的三个下标分别代表区汇所处的时间段、区汇所属子域、法向流量控制分段。

对于子域1的隔水外边界AB、子域2的隔水外边界CD，应满足通过法向流量控制分段Am、mB、Cw、wD的法向流量为零，因此有：

Im(Ωt，1，AB，A-Ωt，1，AB，m)+Im(Ωt，1，BC，A-Ωt，1，BC，m)+

Im(Ωt，1，CA，A-Ωt，1，CA，m)+Ψt，1，Am=0

Im(Ωt，1，AB，m-Ωt，1，AB，B)+Im(Ωt，1，BC，m-Ωt，1，BC，B)+

Im(Ωt，1，CA，m-Ωt，1，CA，B)+Ψt，1，mB=0

Im(Ωt，2，CD，C-Ωt，2，CD，w)+Im(Ωt，2，DA，C-Ωt，2，DA，w)+

Im(Ωt，2，AC，C-Ωt，2，AC，w)+Ψt，2，Cw=0

Im(Ωt，2，CD，w-Ωt，2，CD，D)+Im(Ωt，2，DA，w-Ωt，2，DA，D)+

Im(Ωt，2，AC，w-Ωt，2，AC，D)+Ψt，2，wD=0

t=Δt1，Δt2

(12)

对于子域2的法向定流量外边界DA应满足通过法向流量控制分段Ds、sA的法向流量为f0，因此有：

Im(Ωt，2，CD，D-Ωt，2，CD，s)+Im(Ωt，2，DA，D-Ωt，2，DA，s)+

Im(Ωt，2，AC，D-Ωt，2，AC，s)+Ψt，2，Ds=f0

Im(Ωt，2，CD，s-Ωt，2，CD，A)+Im(Ωt，2，DA，s-Ωt，2，DA，A)+

Im(Ωt，2，AC，s-Ωt，2，AC，A)+Ψt，2，sA=f0

t=Δt1，Δt2

(13)

对于子域1的定水头外边界BC上的控制点3、4应满足其水位为h0，因此有：

[Re(Ωt，1，AB，k+Ωt，1，BC，k+Ωt，1，CA，k)+Φt，1，k]/

(K1H1)+ht，1，k=h0

t=Δt1，Δt2，k=3，4

(14)

对于子域1内的区汇控制点a、b、c由式(8) —式(9)有：

bt，1，brb，a+bt，1，crc，a+Bt，1=S1Δht，1，a/

Δtt，bt，1，ara，b+bt，1，crc，b+Bt，1=S1Δht，1，b/Δtt

bt，1，ara，c+bt，1，brb，c+Bt，1=S1Δht，1，c/

Δtt，bt，1，a+bt，1，b+bt，1，c=0

t=Δt1，Δt2

(15)

式中，b为关于区汇控制点的参数，b的三个下标分别代表区汇所处的时间段、区汇所在的子域、区汇控制点；B为关于区汇函数的常数，B的两个下标分别代表区汇所处的时间段、区汇所在的子域；r为其两个下标所代表的区汇控制点间的直线距离，m；Δh为在区汇所处的时间段内某个区汇控制点水位的改变值，Δh的三个下标分别代表区汇所处的时间段、区汇所在的子域、区汇控制点，m；Δt为其下标所代表的时间段的始末时间差值，d；S1为子域1所在的含水层的储水系数。

类似地，对于子域2内的区汇控制点d、e、f由式(8)—式(9)有：

bt，2，ere，d+bt，2，frf，d+Bt，2=S2Δht，2，d/

Δtt，bt，2，drd，e+bt，2，frf，e+Bt，2=S2Δht，2，e/Δtt

bt，2，drd，f+bt，2，ere，f+Bt，2=S2Δht，2，f/

Δtt，bt，2，d+bt，2，e+bt，2，f=0

t=Δt1，Δt2

(16)

由式(10)—式(16)可知，在每个时间段，依据每个隔水外边界或法向定流量外边界的两个法向流量控制分段可以生成2个方程，依据每个定水头外边界上的两个控制点可以生成2个方程，依据内边界上的两个法向流量控制分段和两个控制点可以生成4个方程，依据每个区汇的3个控制点可以生成4个方程，每个线汇包含2个未知参数，每个区汇包含4个未知参数，因此，每个时间段的模型方程组包含20个未知参数和20个方程。

2.2 模型求解

上述生成的模型方程组中包含的未知参数分为关于线汇强度的参数和关于区汇控制点的参数两类，为线性方程组，可以用求解线性方程组的各种程序求解。求解以后，研究区内任意一点的水位可以通过该研究区内的线汇流量势及区汇流量势叠加求得。

将子域的每条边界剖分为更短的条段或在每个代表子域边界的线汇上设置更多的控制点，并设置更多的区汇控制点和更小的时间步长可以使边界条件得到更精确的满足，但模型方程组未知参数和方程数量会增加进而需要更长的求解时间。

3 现场应用

蔚县矿区单候煤矿5号煤的底板充水含水层为奥陶系灰岩岩溶裂隙含水层，该含水层的富水性强弱是决定5号煤能否安全开采的重要因素之一[8-10]。现需依照《煤矿防治水细则》中突水系数T=0.06MPa/m情形下计算奥灰含水层水位低于临界安全水位所需要疏降的水量[11-13]。研究区模型如图2所示，将模拟范围外扩至煤矿所在的水文地质单元以利用断层等作为研究区的外边界，进而提高模拟结果的可靠性[14，15]。

图2 研究区模型

依据奥陶系灰岩的岩溶裂隙发育程度和渗透系数将研究区划分为三个子域。本次模拟在AnAqSim[16]中进行，AnAqSim是专门的子域解析元素法模拟软件，可以高效地进行模型构建及求解，并进行结果分析及成图。子域的外边界采用三阶线汇，内边界采用四阶线汇，子域区汇控制点间距为1000m，模拟历时(0，10)、(10，30)、(30，70)三个时间段(单位：天)，三个时间段的时间步长数依次为20、15、10，每个时间步长历时是上个时间步长的1.5倍，每个时间步长的方程组包含1845个方程，1845个未知参数，为验证子域解析元素法的有效性及评估模拟精度，6个疏降水控制孔的位置和各孔疏降水量与采用Feflow软件模拟该研究区时一致，子域解析元素法与Feflow模拟结果对比见表1，疏降前后奥灰水位观测孔水位变化如图3所示。

表1 模拟结果对比

图3 疏降水前后奥灰水位观测孔水位变化

从与Feflow的模拟结果对比可以发现，在6个疏降水控制孔的总疏降水量为1300m3/h的条件下，其中，zk2孔疏水量550m3/h，zk3孔疏水量750m3/h，其余四孔不带压，无需疏降水，zk2孔与zk3孔的水位在57d(Feflow模拟结果为60d)时均降至奥灰临界安全水位以下，前30d子域解析元素法模拟水位较Feflow模拟水位高，但随着时间增加两者模拟水位差变小，在60d附近获得了比较接近的模拟结果，通过图3可以看出，经过疏降水四个奥灰水位观测孔水位均处于奥灰临界安全水位以下。