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

任曉波，劉守強(qiáng)，吳瑞芳

(1.中國礦業(yè)大學(xué)(北京) 地球科學(xué)與測繪工程學(xué)院，北京 100083； 2.國家煤礦水害防治工程技術(shù)研究中心，北京 100083；3.河北省煤田地質(zhì)勘查院，河北 邢臺 054000)

目前煤礦地下水流場的數(shù)值模擬主要采用有限差分、有限元等方法，本文旨在探索一種適用于煤礦地下水流場模擬的、效率更高且可為下一步可能需要進(jìn)行的精細(xì)模擬提供指導(dǎo)(比如指明需要重點(diǎn)研究的區(qū)域以及收集哪些含水層的現(xiàn)場試驗(yàn)數(shù)據(jù)等)的模擬方法。本文首先構(gòu)建基于子域解析元素法的煤礦地下水流場模型，然后通過煤礦疏降水量預(yù)測實(shí)例分析其模擬煤礦地下水流場的可行性及精度。子域解析元素法[1-3]是近年來在解析元素法[4-7]的基礎(chǔ)上發(fā)展起來的一種比較適合處理局部地下水流問題的模擬方法。解析元素法的理論基礎(chǔ)是通過疊加滿足地下水流場控制方程和邊界條件的線性解來實(shí)現(xiàn)對地下水流場的模擬。地下水動力學(xué)中用于模擬抽水井附近的定水頭邊界或隔水邊界的鏡像法便是比較簡單的通過線性解的疊加來獲得地下水流場的解析解，即通過在流場中添加各種水力要素(在鏡像法中為虛擬井)以使得流場的邊界條件得以滿足，只不過在解析元素法中這種線性解的疊加方式更加復(fù)雜，疊加的線性解的數(shù)目更加龐大以至于難以通過人工計(jì)算來獲得模擬結(jié)果。

1 子域解析元素法

1.1 線匯

經(jīng)典的地下水流場分析中一般以水位h關(guān)于位置及時間的函數(shù)及偏導(dǎo)數(shù)來構(gòu)成控制方程和邊界條件表達(dá)式，在點(diǎn)(x，y)處，對于均質(zhì)各向同性平面二維承壓穩(wěn)定流有：

式中，Qx、Qy為單寬流量分量，m2/d；qx、qy為達(dá)西滲流速度分量，m/d；H為含水層的厚度，m；K為含水層的滲透系數(shù)，m/d；h為含水層的水位，m。

式(1)中最后一項(xiàng)中的KHh為解析元素法中承壓水流的流量勢Φ，因此有：

在解析元素法中，流量勢為描述地下水流場的基本單位之一，同時為了減少模擬計(jì)算中的變量數(shù)并充分利用復(fù)變函數(shù)中解析函數(shù)的特有性質(zhì)，解析元素法對流場的模擬求解主要在復(fù)平面中進(jìn)行。在復(fù)平面z=(x，y)上，式(2)中的流量勢Φ(z)與其共軛調(diào)和函數(shù)Ψ(z)可組成解析函數(shù)Ω(z)。

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

(3)

由解析函數(shù)的性質(zhì)可知，Φ(z)與Ψ(z)在流場內(nèi)處處正交，因此可以分別用Φ(z)與Ψ(z)表示穩(wěn)定流中的等水位線與流線，Ψ(z)即為流函數(shù)，Ω(z)為復(fù)勢，二維承壓穩(wěn)定井流的復(fù)勢可表示為：

式中，Q為井流量，m3/h；zw為井心坐標(biāo)；C為常數(shù)。

用于模擬流場邊界的一階線匯復(fù)勢表達(dá)式為[4，5]：

式中，z1、z2分別為線匯的首尾端點(diǎn)；Z=(X，Y)為基于線匯中心進(jìn)行坐標(biāo)變換后的局部坐標(biāo)系；σ為線匯的強(qiáng)度，m2/d。

式(5)中線匯的強(qiáng)度σ為常數(shù)，當(dāng)σ為關(guān)于位置的非線性三次函數(shù)時，線匯的復(fù)勢表達(dá)式為[3]：

σ表達(dá)式的次數(shù)加1為其對應(yīng)線匯的階數(shù)，因此式(6)表示的是四階非線性強(qiáng)度線匯。式(6)實(shí)部所代表的線匯流量勢在跨越線匯時流量勢保持連續(xù)，同時，式(6)虛部所代表的線匯流函數(shù)在跨越線匯時流函數(shù)值發(fā)生改變，這與式(5)所代表的常數(shù)強(qiáng)度線匯具有相同的流量勢與流函數(shù)變化特征，因此高階線匯也可用于流場邊界的模擬，且較常數(shù)強(qiáng)度線匯在模擬邊界時具有更高的靈活性和精度。

由式(2)、式(3)可知，由線匯控制的地下水流場中的水位可表示為Re(Ω)/(KH)，其中，Re(Ω)為線匯復(fù)勢的實(shí)部。通過流場中兩點(diǎn)z1、z2間(直線段z1z2不與線匯及線匯所在的負(fù)半軸相交)的流量可表示為|Im(Ω(z1))-Im(Ω(z2))|，其中，Im(Ω(z1))、Im(Ω(z2))分別為點(diǎn)z1、z2處線匯復(fù)勢的虛部。

1.2 子域與區(qū)匯

類似于鏡像法，解析元素法模型沒有外邊界，子域解析元素法采用解析元素法的理論基礎(chǔ)并設(shè)置流場的外邊界，在流場內(nèi)部依據(jù)含水層參數(shù)(如滲透系數(shù)等)不同將其劃分為若干個子域，流場的外邊界和子域間的內(nèi)邊界均用線匯進(jìn)行模擬。對于非穩(wěn)定流，可將?h/?t有限差分為Δh/Δt[1]，二維承壓平面非穩(wěn)定流控制方程可表示為：

式中，γ為單位面積的承壓含水層單位時間內(nèi)彈性儲(釋)水量，m/d；S為承壓含水層的儲水系數(shù)，無量綱。

子域解析元素法采用區(qū)匯(area-sink)模擬由于承壓含水層彈性儲(釋)水造成的水位變化，區(qū)匯的流量勢表示如下：

區(qū)匯用徑向基插值函數(shù)表示，因此有：

式中，ri為待求流量勢的某點(diǎn)與區(qū)匯中第i(i=1，2，…n)個控制點(diǎn)的距離，m；r0為待求流量勢的某點(diǎn)與區(qū)匯幾何中心的距離，m；B0，bi(i=1，2，…n)分別為與區(qū)匯幾何中心及第i個控制點(diǎn)相關(guān)的n+1個未知參數(shù)。

2 模型構(gòu)建及求解

2.1 模型方程組

現(xiàn)以一假想研究區(qū)為例來說明如何生成模型方程組，如圖1所示，示例中子域內(nèi)外邊界均采用二階線匯模擬，也即每個線匯的復(fù)勢表達(dá)式中含有兩個關(guān)于線匯強(qiáng)度σ的未知參數(shù)。

圖1 假想流場示意圖

研究區(qū)依據(jù)含水層滲透系數(shù)不同劃分為兩個子域(子域1與子域2)，子域1外邊界為隔水邊界AB、定水頭邊界BC，子域1的內(nèi)邊界為CA，子域2外邊界為隔水邊界CD、法向定流量邊界DA，子域2的內(nèi)邊界為AC，子域1的內(nèi)邊界CA與子域2的內(nèi)邊界AC為兩條重合但相互獨(dú)立的邊界，這兩條內(nèi)邊界分別用不同的線匯來模擬，為簡便起見，子域的每條邊界只用一個二階線匯來表示，每個線匯上各設(shè)置兩個控制點(diǎn)，然后每個線匯依照其兩端端點(diǎn)和控制點(diǎn)可以劃分兩條法向流量控制分段。子域1與子域2中由抽水井W1產(chǎn)生的非穩(wěn)定井流用區(qū)匯(式(7)—式(9))來模擬，每個子域?qū)?yīng)一個不同的區(qū)匯，每個區(qū)匯設(shè)置3個控制點(diǎn)，同時，區(qū)匯的時間步長(式(7)中的Δt)設(shè)置為Δt1=t1-t0，Δt2=t2-t1，(t0，t1)與(t1，t2)為模擬過程經(jīng)歷的兩個時間段。

每個子域內(nèi)部及其內(nèi)外邊界上任意一點(diǎn)的水位與滲流速度由該子域的邊界線匯、區(qū)匯共同決定而與子域外的任何線匯或區(qū)匯均無關(guān)，同時，為了保證在由各子域組成的研究區(qū)范圍內(nèi)滿足達(dá)西滲流場能量與質(zhì)量守恒，在子域間的內(nèi)邊界兩側(cè)需同時滿足水位與法向流量處處相等，因此對于子域1與子域2的公共內(nèi)邊界AC上的控制點(diǎn)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)

式中，Ω為各子域的各個線匯在內(nèi)邊界控制點(diǎn)上產(chǎn)生的復(fù)勢，Ω的四個下標(biāo)分別代表線匯所處的時間段、線匯所屬子域、線匯代表的邊界、控制點(diǎn)；Φ為各子域的區(qū)匯在內(nèi)邊界控制點(diǎn)上產(chǎn)生的流量勢，Φ的三個下標(biāo)分別代表區(qū)匯所處的時間段、區(qū)匯所屬子域、控制點(diǎn)；h為內(nèi)邊界控制點(diǎn)各時段初始水位，h的三個下標(biāo)分別代表時間段、所屬子域、控制點(diǎn)；K1、K2分別為子域1與子域2的滲透系數(shù)，m/d；H1、H2分別為子域1與子域2含水層厚度，m。

對于子域1與子域2的公共內(nèi)邊界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)

式中，Ψ為各子域的區(qū)匯在各個法向流量控制分段上產(chǎn)生的流量，Ψ的三個下標(biāo)分別代表區(qū)匯所處的時間段、區(qū)匯所屬子域、法向流量控制分段。

對于子域1的隔水外邊界AB、子域2的隔水外邊界CD，應(yīng)滿足通過法向流量控制分段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應(yīng)滿足通過法向流量控制分段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上的控制點(diǎn)3、4應(yīng)滿足其水位為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內(nèi)的區(qū)匯控制點(diǎn)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為關(guān)于區(qū)匯控制點(diǎn)的參數(shù)，b的三個下標(biāo)分別代表區(qū)匯所處的時間段、區(qū)匯所在的子域、區(qū)匯控制點(diǎn)；B為關(guān)于區(qū)匯函數(shù)的常數(shù)，B的兩個下標(biāo)分別代表區(qū)匯所處的時間段、區(qū)匯所在的子域；r為其兩個下標(biāo)所代表的區(qū)匯控制點(diǎn)間的直線距離，m；Δh為在區(qū)匯所處的時間段內(nèi)某個區(qū)匯控制點(diǎn)水位的改變值，Δh的三個下標(biāo)分別代表區(qū)匯所處的時間段、區(qū)匯所在的子域、區(qū)匯控制點(diǎn)，m；Δt為其下標(biāo)所代表的時間段的始末時間差值，d；S1為子域1所在的含水層的儲水系數(shù)。

類似地，對于子域2內(nèi)的區(qū)匯控制點(diǎn)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)可知，在每個時間段，依據(jù)每個隔水外邊界或法向定流量外邊界的兩個法向流量控制分段可以生成2個方程，依據(jù)每個定水頭外邊界上的兩個控制點(diǎn)可以生成2個方程，依據(jù)內(nèi)邊界上的兩個法向流量控制分段和兩個控制點(diǎn)可以生成4個方程，依據(jù)每個區(qū)匯的3個控制點(diǎn)可以生成4個方程，每個線匯包含2個未知參數(shù)，每個區(qū)匯包含4個未知參數(shù)，因此，每個時間段的模型方程組包含20個未知參數(shù)和20個方程。

2.2 模型求解

上述生成的模型方程組中包含的未知參數(shù)分為關(guān)于線匯強(qiáng)度的參數(shù)和關(guān)于區(qū)匯控制點(diǎn)的參數(shù)兩類，為線性方程組，可以用求解線性方程組的各種程序求解。求解以后，研究區(qū)內(nèi)任意一點(diǎn)的水位可以通過該研究區(qū)內(nèi)的線匯流量勢及區(qū)匯流量勢疊加求得。

將子域的每條邊界剖分為更短的條段或在每個代表子域邊界的線匯上設(shè)置更多的控制點(diǎn)，并設(shè)置更多的區(qū)匯控制點(diǎn)和更小的時間步長可以使邊界條件得到更精確的滿足，但模型方程組未知參數(shù)和方程數(shù)量會增加進(jìn)而需要更長的求解時間。

3 現(xiàn)場應(yīng)用

蔚縣礦區(qū)單候煤礦5號煤的底板充水含水層為奧陶系灰?guī)r巖溶裂隙含水層，該含水層的富水性強(qiáng)弱是決定5號煤能否安全開采的重要因素之一[8-10]?，F(xiàn)需依照《煤礦防治水細(xì)則》中突水系數(shù)T=0.06MPa/m情形下計(jì)算奧灰含水層水位低于臨界安全水位所需要疏降的水量[11-13]。研究區(qū)模型如圖2所示，將模擬范圍外擴(kuò)至煤礦所在的水文地質(zhì)單元以利用斷層等作為研究區(qū)的外邊界，進(jìn)而提高模擬結(jié)果的可靠性[14，15]。

圖2 研究區(qū)模型

依據(jù)奧陶系灰?guī)r的巖溶裂隙發(fā)育程度和滲透系數(shù)將研究區(qū)劃分為三個子域。本次模擬在AnAqSim[16]中進(jìn)行，AnAqSim是專門的子域解析元素法模擬軟件，可以高效地進(jìn)行模型構(gòu)建及求解，并進(jìn)行結(jié)果分析及成圖。子域的外邊界采用三階線匯，內(nèi)邊界采用四階線匯，子域區(qū)匯控制點(diǎn)間距為1000m，模擬歷時(0，10)、(10，30)、(30，70)三個時間段(單位：天)，三個時間段的時間步長數(shù)依次為20、15、10，每個時間步長歷時是上個時間步長的1.5倍，每個時間步長的方程組包含1845個方程，1845個未知參數(shù)，為驗(yàn)證子域解析元素法的有效性及評估模擬精度，6個疏降水控制孔的位置和各孔疏降水量與采用Feflow軟件模擬該研究區(qū)時一致，子域解析元素法與Feflow模擬結(jié)果對比見表1，疏降前后奧灰水位觀測孔水位變化如圖3所示。

表1 模擬結(jié)果對比

圖3 疏降水前后奧灰水位觀測孔水位變化

從與Feflow的模擬結(jié)果對比可以發(fā)現(xiàn)，在6個疏降水控制孔的總疏降水量為1300m3/h的條件下，其中，zk2孔疏水量550m3/h，zk3孔疏水量750m3/h，其余四孔不帶壓，無需疏降水，zk2孔與zk3孔的水位在57d(Feflow模擬結(jié)果為60d)時均降至奧灰臨界安全水位以下，前30d子域解析元素法模擬水位較Feflow模擬水位高，但隨著時間增加兩者模擬水位差變小，在60d附近獲得了比較接近的模擬結(jié)果，通過圖3可以看出，經(jīng)過疏降水四個奧灰水位觀測孔水位均處于奧灰臨界安全水位以下。

4 結(jié) 論

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

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

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