基于子域解析元素法的煤礦疏降水量預測研究

任曉波，劉守強，吳瑞芳

(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]?，F需依照《煤礦防治水細則》中突水系數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可以看出，經過疏降水四個奧灰水位觀測孔水位均處于奧灰臨界安全水位以下。

4 結 論

1)通過與Feflow的模擬結果對比以及疏降水前后奧灰水位觀測孔水位變化說明子域解析元素法是一種有效的可以用于煤礦地下水流場模擬的方法，同時，增加各時間段的步長數量或采用更高階的線匯模擬子域內外邊界以及將邊界剖分為更短長度可進一步提高模擬精度。

2)對于同等復雜程度模型，子域解析元素法的模型求解耗時較短，且可靈活地進行模型調整，例如改變邊界的性質或位置時只需在模型中調整與邊界對應的線匯屬性即可。

3)子域解析元素法的模擬結果可以揭示研究區流場的總體特征，進而指明需要重點研究的區域及需要收集的水文地質勘探資料，可為下一步精細復雜模擬做好準備，同時應用子域解析元素法時需要在模擬精度和計算效率上做好平衡，因此，子域解析元素法是煤礦地下水流場模擬方法的有益補充。