史鵬鈺，宗一杰，滕開慶，劉健軍，肖 良

(1.廣西大學(xué) 土木建筑工程學(xué)院，廣西 南寧 530004；2.廣西防災(zāi)減災(zāi)與工程安全重點(diǎn)實(shí)驗(yàn)室，廣西 南寧 530004；3.廣西大學(xué) 工程防災(zāi)與結(jié)構(gòu)安全教育部重點(diǎn)實(shí)驗(yàn)室，廣西 南寧 530004)

定水頭抽水試驗(yàn)是一種特殊的抽水試驗(yàn)，其井中水頭在抽水時需保持穩(wěn)定[1]，常被應(yīng)用于滲透系數(shù)較低的地層以及抽水量難以維持穩(wěn)定的工況[2]。定水頭抽水試驗(yàn)對于確定水力參數(shù)或解決垃圾填埋場引起的環(huán)境破壞問題具有重要意義[2-3]。自1950 年以來，水文地質(zhì)學(xué)家一直在研究由定水頭抽水試驗(yàn)引起的抽水井流問題，C.E.Jacob 等[4]第一個提出了用于研究定水頭抽水試驗(yàn)期間的達(dá)西流流動特性的解析解。此后，水文學(xué)家們在定水頭抽水井流的理論和物理模型方面進(jìn)行了大量有價值的工作[2,5-8]。

一般情況下，抽水引起的地下水流動遵循達(dá)西定律，符合線性流動假設(shè)[9]。然而，許多學(xué)者發(fā)現(xiàn)，當(dāng)水力梯度較大或較小時，流量與水力梯度會呈非線性變化關(guān)系。這種非線性抽水井流被定義為非達(dá)西流[10-11]，通常由Izbash 方程或Forchheimer 方程描述，并根據(jù)現(xiàn)場條件在2 個方程之間進(jìn)行選擇[12-13]。在過去10 年中，大量學(xué)者基于上述2 種方法，通過原位試驗(yàn)研究了定水頭抽水引起的非達(dá)西流流動特性問題[5,8,12,14-15]。例如，P.M.Quinn 等[15]在加拿大安大略省圭爾夫市附近的裂縫性白云巖和砂巖中進(jìn)行了定水頭原位試驗(yàn)，并指出使用線性特征的達(dá)西流來擬合實(shí)驗(yàn)數(shù)據(jù)可能會產(chǎn)生很大的誤差。Dan Hancheng 等[6]基于定水頭試驗(yàn)的降深數(shù)據(jù)，對無黏結(jié)級配骨料的參數(shù)進(jìn)行評估，發(fā)現(xiàn)即使在低水力梯度下非達(dá)西流現(xiàn)象也很普遍。Liu Zhongyu 等[8]使用改進(jìn)的滲透率測試裝置在飽和黏土含水層中進(jìn)行定水頭測試，發(fā)現(xiàn)抽水試驗(yàn)會導(dǎo)致水流具有明顯的非線性特性。這些研究認(rèn)為，忽略非達(dá)西效應(yīng)會使抽水井流模擬和確定含水層參數(shù)產(chǎn)生較大的誤差。然而，現(xiàn)有的定水頭抽水模型大多是基于現(xiàn)場或室內(nèi)試驗(yàn)數(shù)據(jù)擬合得出，適用地質(zhì)范圍較為狹窄，且不具備理論依據(jù)。到目前為止，關(guān)于定水頭抽水引起的非達(dá)西流數(shù)值模型研究很少，提出可通用于常見含水層系統(tǒng)的相關(guān)數(shù)值模型是具備現(xiàn)實(shí)價值的。

由文獻(xiàn)評述筆者發(fā)現(xiàn)，現(xiàn)有定水頭抽水試驗(yàn)研究通?；诔樗畬邮浅袎簩拥募僭O(shè)上提出的。然而，世界上許多抽水井都處于非承壓含水層中[16-17]，并且關(guān)于非承壓含水層定水頭抽水井流的研究一直較為缺乏。到目前為止，現(xiàn)有的非承壓含水層抽水井流模型主要分為考慮和不考慮垂向流分量2 類。在考慮垂直流分量模型中，重力釋水不是瞬間釋放的，這進(jìn)一步導(dǎo)致水頭上方的非飽和帶會對下方飽和帶進(jìn)行補(bǔ)給[18]?；诤畬又亓λ窟h(yuǎn)大于彈性儲水量的事實(shí)，在抽水井內(nèi)降深較小或抽水時間足夠長的定水頭試驗(yàn)中，垂直流的影響有限，可以合理使用不考慮垂向流分量的Dupuit 假設(shè)，該假設(shè)認(rèn)為重力儲存為瞬時釋放[19-21]。根據(jù)Dupuit 的調(diào)查，當(dāng)非承壓含水層降深變化較緩時，可以忽略垂直流的影響，并假設(shè)所有補(bǔ)給均來自水平方向。因此，它可用于簡化無約束流動問題，使其在解析上易于處理。在非承壓含水層中，當(dāng)水位降深相對初始水位較小時，可以優(yōu)先考慮使用Dupuit 假設(shè)來模擬流動，以滿足實(shí)際工程的需要[2,13,20,22-24]。

從理論上說，由于抽水井附近的流速足夠大，會對相應(yīng)區(qū)域的流態(tài)造成影響。足夠大的流速可以使雷諾數(shù)大于臨界雷諾數(shù)，從而導(dǎo)致抽水井附近水流呈現(xiàn)非達(dá)西流態(tài)。而隨著徑向距離的增加，雷諾數(shù)隨流速逐漸減小。因此，它可以進(jìn)一步引起非達(dá)西流向達(dá)西流的轉(zhuǎn)變，導(dǎo)致含水層中呈現(xiàn)出兩區(qū)流特性[25-29]。迄今為止，在非承壓含水層中通過定水頭試驗(yàn)引發(fā)的這種兩區(qū)流動機(jī)制仍不清楚，相關(guān)研究不足。

針對上述問題，筆者提出一種非承壓含水層定水頭抽水試驗(yàn)導(dǎo)致的非達(dá)西-達(dá)西兩區(qū)滲流數(shù)值模型。當(dāng)含水層儲水系數(shù)較小時，通過Izbash 方程描述比流量和水力梯度之間的關(guān)系，以進(jìn)行數(shù)學(xué)建模?；贒upuit 假設(shè)，利用有限差分法推導(dǎo)出易于實(shí)際操作的簡潔數(shù)值模型。所提出模型通過與COMSOL Multiphysics 的數(shù)值解進(jìn)行比較，驗(yàn)證滲流模型的可靠性，以期為兩區(qū)流特性分析奠定基礎(chǔ)。

1 兩區(qū)流數(shù)學(xué)模型

圖1 為在非承壓含水層中定水頭抽水試驗(yàn)導(dǎo)致的兩區(qū)滲流模型示意圖，模型考慮了非達(dá)西和達(dá)西2 個區(qū)域的水流流動，其數(shù)學(xué)模型的假設(shè)包括：(1) 該含水層是非承壓、均質(zhì)、各向同性的非承壓含水層，初始水頭為一定值；(2) 抽水井完全穿透非承壓含水層，抽水井內(nèi)的水頭自抽水開始后保持不變；(3) 非達(dá)西區(qū)域(簡稱N 區(qū))為有限區(qū)域，主要發(fā)生在抽水井周圍；達(dá)西區(qū)域(簡稱D 區(qū))為半無限區(qū)域，位于N 區(qū)外圍；(4) 抽水井半徑與有效半徑相同，為一個固定正數(shù)[23]。

圖1 考慮兩區(qū)轉(zhuǎn)換的非承壓含水層定水頭抽水模型Fig.1 Schematic of the constant-head test in an unconfined aquifer considering the two-region transform

基于質(zhì)量守恒定律和Dupuit 假設(shè)，根據(jù)Boussinesq方程將N 區(qū)和D 區(qū)中的抽水流量分別用下式描述。

式中：q1和q2分別為N 區(qū)和D 區(qū)中的比流量，m/s；h1和h2分別為N 區(qū)和D 區(qū)中的水頭，m；上述參數(shù)均為時間t和水平距離r(指某一點(diǎn)與抽水井中心的水平距離)的函數(shù)；S1和S2分別為N 區(qū)和D 區(qū)的儲水系數(shù)。

假設(shè)N 區(qū)和D 區(qū)之間的轉(zhuǎn)換界面位于r=R處，根據(jù)流體連續(xù)性假設(shè)得到邊界處的邊界條件如下：

達(dá)西定律描述了水力坡度(I)和比流量(q)之間的線性關(guān)系。對于多孔介質(zhì)流，研究表明其線性關(guān)系在以下2 種情況下是不成立的。一種是在較低水力梯度下通過低滲介質(zhì)的流動，另一種則是在較高水力梯度下通過粗粒高滲介質(zhì)的流動。J.Bear[30]提到，基于平均直徑的雷諾數(shù)(Re)在Re<1～10 時，達(dá)西定律是適用的。在此范圍內(nèi)，滲流流動可以被歸納為層流，稱為達(dá)西流，否則即可被稱為非達(dá)西流。此外，Re較高時(Re>10)的流動為湍流。對于抽水井附近的流動，水力梯度通常很高，流速隨之增大[31]。

1.1 非達(dá)西區(qū)域

Izbash 方程可用于描述非達(dá)西流的I和q之間的關(guān)系。使用該非線性函數(shù)的一個假設(shè)是，至少在某些流動階段，所涉及的非達(dá)西系數(shù)n和準(zhǔn)滲透系數(shù)kc與空間和時間無關(guān)[32]。其公式表示為：I=cqn,c為常數(shù)，當(dāng)參數(shù) 1

當(dāng)Re>10 時，可認(rèn)為該區(qū)域處于N 區(qū)，其抽水井流為非達(dá)西滲流[25]，N 區(qū)的比流量可用Izbash 方程表示。

式中：K1為N 區(qū)中的準(zhǔn)滲透系數(shù)，(m/s)n1(n1條件下)；n1為N 區(qū)中的非達(dá)西湍流系數(shù)。當(dāng)n1<1 時滲流被稱為前達(dá)西流；當(dāng) 1

模型的初始條件是：

在定水頭抽水試驗(yàn)中，代表完整井的第一類邊界條件如下：

式中：Q為抽水流量，m3/s。

將式(5)代入式(1)可得：

已知N 區(qū)中的流量等于比流量乘以過流斷面面積，如下式所述：

式中：A=2πrh1為井壁處過流斷面面積，m2；在抽水井中h1為已知參數(shù)。

假設(shè)非承壓含水層的降深遠(yuǎn)遠(yuǎn)小于b，則可以基于Dupuit 假設(shè)，在式(8)和式(9)中近似引入來線性化一小部分推導(dǎo)過程。研究表明采用該近似方法可使建模過程大為簡化，同時該近似方法所引起的計(jì)算誤差相對較小，可以忽略不計(jì)[34-35]。因此，式(9)和式(10)可以近似改寫為：

結(jié)合式(11)與式(12)，N 區(qū)中的控制方程可表示如下：

1.2 達(dá)西區(qū)域

當(dāng) 1

式中：K2為D 區(qū)中的滲透系數(shù)，m/s。

初始條件為：

無窮遠(yuǎn)處的邊界條件為：

將式(14)代入式(2)可得：

同理，為了簡化建模過程，假設(shè)非承壓含水層的降深遠(yuǎn)遠(yuǎn)小于b，基于Dupuit 假設(shè)，引入近似關(guān)系來線性化式(17)，可得D 區(qū)中的控制方程可表示如下：

2 基于有限差分的數(shù)值解

有限差分法簡稱FDM，常用于求解大型非線性地下水流運(yùn)動問題，F(xiàn)DM 在不規(guī)則含水層邊界和含水層內(nèi)區(qū)域參數(shù)的緊密空間近似方面具有靈活性[36]。因此，本研究采用FDM 基本思想按時間步長和空間步長對定解區(qū)域進(jìn)行網(wǎng)格劃分[37]，把從無窮遠(yuǎn)處到抽水井處的空間區(qū)域及從抽水開始到抽水結(jié)束的時間區(qū)域用有限個網(wǎng)格節(jié)點(diǎn)代替。

采用FDM 推求上述方程的數(shù)值解，需先對抽水井流模型進(jìn)行如下兩步預(yù)處理。第一步，針對無窮遠(yuǎn)邊界，本文擬采用恒定不變且取值較大的半徑值re來代替，使抽水井流模型可被認(rèn)為是一個因變量為半徑r和時間t的一維非穩(wěn)定流模型。第二步，將時間區(qū)間[0,t] 離散為N個步長，徑向空間 [rw,re] 離散為J個步長。其中，任何一個時間坐標(biāo)上節(jié)點(diǎn)m處的時間值tm必滿足 0 ≤tm≤t,m=0,1,···,N，任何一個徑向坐標(biāo)上節(jié)點(diǎn)j處的徑向距離rj需滿足rw≤rj≤re,j=0,1,···,J。其步長大小的具體值表示為：Δt=t/N，Δr=(re-rw)/J。由于在抽水井附近水位降深變化較大，因此在離抽水井較近的地方空間步長應(yīng)盡可能小，而離抽水井較遠(yuǎn)的地方空間步長可以適當(dāng)加大。對N 區(qū)和D 區(qū)的控制方程及其所有邊界條件進(jìn)行離散化，采用向后差分的方法，依次得出N 區(qū)和D 區(qū)的水頭。

2.1 非達(dá)西區(qū)域

N 區(qū)控制方程式(13)為一個二階齊次微分方程，基于FDM 的基本理論，將控制方程離散化，得到離散化后的控制方程如下：

其中：

對于N 區(qū)的邊界條件式(7)和式(8)，離散化后的邊界條件如下：

2.2 達(dá)西區(qū)域

同理，D 區(qū)控制方程式(18)為一個二階齊次微分方程，基于FDM 的基本理論離散化后可得：

其中：

對于D 區(qū)的邊界條件式(15)和式(16)，離散化后的邊界條件如下：

利用Matlab 軟件，將式(19)-式(24)與式(25)-式(30)分別導(dǎo)入，輸入合適的步長，即可由軟件自動計(jì)算并分別得出N 區(qū)及D 區(qū)各個節(jié)點(diǎn)上的流量及水頭值。通過在rj≤R的徑向節(jié)點(diǎn)上利用式(19)以及在rj>R的徑向節(jié)點(diǎn)上利用式(25)計(jì)算得出的結(jié)果，調(diào)用Matlab 中名為Equations and System Solver 的內(nèi)置模塊，即可得到在假設(shè)的定水頭試驗(yàn)條件下任意點(diǎn)的流量及水頭情況，N 區(qū)與D 區(qū)求解矩陣分別如下：

定水頭抽水試驗(yàn)可以使抽水井附近的水頭在抽水早期迅速下降，并在短時間內(nèi)接近hw。在此之后，水頭變化的幅度會大大減小。因此，除了在抽水剛開始時的短時間內(nèi)，其余階段基于Dupuit 假設(shè)提出的解是可以接受的。

在定水頭抽水試驗(yàn)中，由于流動的湍流程度會隨著與泵井徑向距離的增加而減弱，因此可將N 區(qū)定義為Re大于臨界雷諾數(shù)Rec的區(qū)域，并且通過判斷Re=Rec|r→R來確定該時刻的R值。此外，還應(yīng)該注意的是，所提出的解也可通過假設(shè)R=rw來評估純達(dá)西流或設(shè)定R=∞ 調(diào)查純非達(dá)西流的降深。計(jì)算雷諾數(shù)的公式[25]為Re=qd/υ，其中d為介質(zhì)的特征直徑，m；υ 為水的運(yùn)動黏度，m2/s，其值通常為 1 ×106m2/s?？梢缘玫竭_(dá)西區(qū)和非達(dá)西區(qū)轉(zhuǎn)換界面處的Rec表達(dá)式如下：

根據(jù)相關(guān)研究，通常認(rèn)為Rec=10?；谟?jì)算得到的水頭h值，可以通過Matlab 進(jìn)而推求兩區(qū)流界面位置R值。擬通過計(jì)算得出的達(dá)西區(qū)水頭h2進(jìn)行迭代計(jì)算，其具體步驟如下：

(1) 選取初始值。假設(shè)在抽水前，井流為達(dá)西流，即tm=0 時，時間步長數(shù)m=0，兩區(qū)流交界面位置為R=Rm=0。

(2) 計(jì)算流速。計(jì)算第m+1 個時間步長tm+1時的水頭并計(jì)算每個節(jié)點(diǎn)處的流速及雷諾數(shù)Re，然后根據(jù)臨界雷諾數(shù)Rec確定兩區(qū)流交界面位置Rmc。

(3) 收斂性判斷。δ為一個給定的判別標(biāo)準(zhǔn)，若 |Rm-Rmc|<δ，則認(rèn)為第m個時間步長的兩區(qū)流交界面位置R=Rm。然后計(jì)算下一時間步長兩區(qū)流交界面位置：令m=m+1，若mδ，令Rm=Rmc。

(4) 在新交界面處計(jì)算流速。根據(jù)計(jì)算得出的水頭值，計(jì)算每個時間節(jié)點(diǎn)處的流速與雷諾數(shù)，由Rec確定兩區(qū)流交界面位置Rmc，返回步驟(3)進(jìn)行判斷，進(jìn)入下一步迭代計(jì)算，直至滿足δ<5% 的精度要求，迭代終止。

3 模型驗(yàn)證與對比

鑒于目前缺少非承壓含水層定水頭抽水計(jì)算模型的相關(guān)研究，難以找到相關(guān)工況完全一致的模型研究或?qū)崪y案例，為驗(yàn)證本文提出有限差分?jǐn)?shù)值模型的可靠性，本章將提出的非承壓定水頭抽水井流模型退化為承壓定水頭達(dá)西抽水井流模型，并與現(xiàn)有的經(jīng)典承壓定水頭達(dá)西抽水井流解析模型(Jacob 和Lohman 模型)進(jìn)行對比分析[4]。同時，通過假設(shè)案例，對所提出有限差分模型與基于COMSOL 有限元模型的計(jì)算結(jié)果進(jìn)行對比分析，以驗(yàn)證所提出的有限差分模型的可靠性。為方便計(jì)算，本文基于Matlab 軟件實(shí)現(xiàn)對提出模型的自動化求解。

3.1 與Jacob 和Hantush 解析解比較

在抽水過程中，通過控制含水層壓力水頭始終大于含水層厚度，將模型退化成承壓含水層定水頭抽水井流模型。然后，再令n1=1 將模型退化成達(dá)西流，采用Jacob 和Lohman 模型的參數(shù)，將比流量q及抽水時間t無量綱處理并代入本模型中，可得到對比結(jié)果如圖2 所示。從圖中可以看出，本文在定水頭抽水的承壓含水層達(dá)西流情況下的有限差分解與Jacob 和Lohman 模型解高度吻合，說明在模型推導(dǎo)過程中由數(shù)值差分帶來的計(jì)算誤差可忽略不計(jì)[38]。

圖2 有限差分解與Jacob 和Lohman 解對比Fig.2 Comparison between the finite difference solution with Jacob and Lohman solution

3.2 與COMSOL 有限元解比較

非承壓含水層非達(dá)西定水頭抽水徑流模型采用有限元模型基于COMSOL Multiphysics 軟件建立，該軟件已被證明是最可靠的多物理場仿真軟件之一[39]。COMSOL 模型是基于N 區(qū)的式(6)-式(9)及D 區(qū)的式(15)-式(17)推導(dǎo)而來。

需要注意的是，式(16)中無窮遠(yuǎn)邊界的距離在COMSOL 解中用r=1 000 m 表示，假設(shè)非承壓含水層是砂質(zhì)介質(zhì)，基于現(xiàn)有研究可知非承壓含水層的儲水系數(shù)為 0 .01～0.03，非承壓含水層中細(xì)砂的滲透系數(shù)為 1～5 m/d[20]。因此，案例A 參數(shù)選取見表1。考慮到Dupuit 假設(shè)的要求，抽水井內(nèi)部降深需要取一個相對較小的值，故取b=12 m,hw=9 m。表1 中的案例B、C 和D 用于下文討論部分的參數(shù)選取。

表1 假設(shè)案例的參數(shù)值Table 1 Parameters of each hypothetical cases

圖3 展示了兩區(qū)流(R為 變量)、純達(dá)西流(R=rw=0.2 m) 和純非達(dá)西流(R=∞)的水頭-時間曲線對比。純非達(dá)西流COMSOL 解近似為R=1 000 m。應(yīng)該注意到，K1和S1用于純非達(dá)西流動的模擬，而K2和S2用于純達(dá)西流動的模擬。這一對比清楚地表明，在所有3 種流動條件下，模擬結(jié)果與有限差分解之間的偏差都很小。在抽水開始后第一個小時內(nèi)出現(xiàn)一定偏差，可認(rèn)為是忽略了所提出的有限差分解的垂直流動而引起，這種偏差隨著時間的推移逐漸消失，模擬結(jié)果進(jìn)一步驗(yàn)證了所提出的有限差分模型用于兩區(qū)流動的水頭模擬的有效性。

圖3 不同R 值下有限差分解與數(shù)值解的水頭-時間曲線對比Fig.3 Comparison of head-time curves between finite difference solution and numerical solution for different R

此外，兩區(qū)流的水頭明顯小于純非達(dá)西流的水頭，而大于純達(dá)西流的水頭(圖3)，說明非達(dá)西流對水頭下降有負(fù)影響。這是因?yàn)镽值越大表示非達(dá)西范圍越大，抽水井附近的非承壓含水層湍流范圍更廣。相比于純達(dá)西流，純非達(dá)西流的湍流程度更高。紊流會導(dǎo)致過水?dāng)嗝媪魉僮兇螅鄳?yīng)比流量變大，含水層從外界得到水的能力增強(qiáng)，無窮遠(yuǎn)處對抽水井附近含水層補(bǔ)給變大，進(jìn)而延緩因抽水導(dǎo)致的含水層水位下降[40]。因此，在抽水后期，純非達(dá)西流的水頭能夠明顯大于純達(dá)西流的水頭。

圖4 展示了基于不同模型R-t曲線的對比。結(jié)果表明，不同時間下本文提出的有限差分模型解與COMSOL 有限元模型解擬合程度較高，驗(yàn)證了通過迭代法計(jì)算出的R值的可靠性。而R值隨時間推移不斷減小，說明非達(dá)西流對兩區(qū)流動的影響隨著抽水的繼續(xù)而減弱。在抽水開始后，非承壓含水層水頭逐漸下降，接近于hw，它進(jìn)一步導(dǎo)致水力梯度和比流量的下降。比流量越小，流速越慢，雷諾數(shù)越小，最終導(dǎo)致N 區(qū)消失。當(dāng)抽水時間足夠長時，N 區(qū)可以完全消失，兩區(qū)流可以轉(zhuǎn)化為純達(dá)西流。

圖4 有限差分解界面距離-時間曲線與數(shù)值解的比較Fig.4 Interface distance-time curve of finite difference solution vs numerical solution

3.3 與實(shí)測數(shù)據(jù)比較

L.Jones 等[41]在美國威斯康星州開展了定水頭抽水試驗(yàn)，所調(diào)查的非承壓含水層由厚度為2.5 m 的含泥量低、含砂量高的淺風(fēng)化冰磧層組成。抽水井有效半徑為0.051 m，且完全貫穿非承壓含水層，其底部由礫石充填。抽水共進(jìn)行24 h，抽水井內(nèi)定水頭為1.5 m[42]。兩組降深數(shù)據(jù)分別在徑向距離為0.87、3.62 m 的觀測井中采集。根據(jù)C.S.Chen 等[2]的研究，水平方向的給水度和滲透系數(shù)分別為 0.014～0.042 和0.241～0.362 m/d。因此，基于式(19)、式(25)經(jīng)對比分析，選取參數(shù)S1=S2=0.031，K1=0.285 m/d，K2=0.320 m/d，n1=1.3時，可得到實(shí)測數(shù)據(jù)與本文模型解的對比曲線擬合度良好，如圖5 所示。但與圖3 的結(jié)果類似，有限差分解模擬曲線在抽水期間與實(shí)測數(shù)據(jù)的偏差很小，僅在抽水中期會出現(xiàn)一定程度的誤差，且這種偏差隨著時間的增加逐漸消失。圖5 進(jìn)一步驗(yàn)證了所提出的有限差分模型在實(shí)際工程中的適用性。

圖5 有限差分解與實(shí)測數(shù)據(jù)對比Fig.5 Finite difference solution vs measured data

4 討論

在實(shí)際抽水試驗(yàn)中，一般認(rèn)為抽水附近區(qū)域流態(tài)為達(dá)西流，以方便概化模型及預(yù)測井內(nèi)降深。本研究提出的定水頭抽水兩區(qū)流數(shù)值解，將抽水井附近的有限區(qū)域視為非達(dá)西區(qū)域，遠(yuǎn)離抽水井的半無限區(qū)域視為達(dá)西區(qū)域，從而將非承壓含水層的地下水滲流根據(jù)流態(tài)區(qū)分為2 個流域。相比于純達(dá)西流模型，該兩區(qū)流模型更符合實(shí)際工況，計(jì)算結(jié)果精度更高，但相應(yīng)計(jì)算過程也較復(fù)雜。該模型適用于沙土或碎石比例較大、滲透性較強(qiáng)的含水層，其滲流流速相對較快，Re相對較高，在抽水井附近極易形成非達(dá)西流。但在黏土等滲透能力弱的含水層，滲流流速較小，非達(dá)西范圍相對小，該模型將不再適用。在實(shí)際工程應(yīng)用中，通常會著重觀測井內(nèi)水頭及抽水量的動態(tài)變化。因此，該模型有必要討論偏離達(dá)西程度對流態(tài)的影響，進(jìn)而判斷不同非達(dá)西湍流系數(shù)n1對水頭及抽水井抽水量的影響，并討論井內(nèi)水頭hw對抽水井抽水速率的影響。準(zhǔn)確評估抽水流量Q的動態(tài)變化對于設(shè)計(jì)定水頭抽水試驗(yàn)至關(guān)重要。在下文中，基于假設(shè)案例B、C 和D(參數(shù)見表1)揭示非達(dá)西流的n1、兩區(qū)流轉(zhuǎn)換界面位置R和井內(nèi)水頭hw對抽水速率動態(tài)變化的影響，以及徑向距離對斷面流量的影響。

圖6 展示了案例B 中不同時刻hw的Q-t曲線的模擬結(jié)果。結(jié)果表明，固定時刻下Q隨著hw的增加不斷降低，在抽水早期不同曲線對應(yīng)的Q之間偏差較大，隨著抽水的繼續(xù)，流量偏差逐漸減小。Q的變化與hw呈負(fù)相關(guān)關(guān)系，在抽水早期Q-t曲線下降速度較快，隨著抽水的繼續(xù)Q下降的趨勢越來越平緩，符合hw的小幅降低會使流速大幅度升高的事實(shí)。流量變化實(shí)際上是由抽水井內(nèi)部水頭下降而引起的水力梯度變化造成的，因此，較大的hw意味著較小的比流量，從而導(dǎo)致抽水流量較小。

圖6 泵井內(nèi)不同水頭下流量的動態(tài)變化Fig.6 Dynamic development of the flow rates with different hydraulic head inside the pumping well

圖7 展示了案例C 中非達(dá)西湍流系數(shù)對泵送流量的影響，其中圖7a 和圖7b 分別為固定徑向距離r=2 m 時的Q-t曲線和相同抽水時間t=1 d 時的h-r曲線。從圖7a 中可以發(fā)現(xiàn)，在抽水早期，固定時間點(diǎn)的Q隨著n1的增加而增加。而隨著抽水的繼續(xù)，不同n1對應(yīng)的流量偏差不斷減小。從理論上講，這是可以理解的。抽水井附近含水層補(bǔ)給水量以重力釋水為主，如式(10)所示，假設(shè)其等于比流量乘以過流斷面面積，比流量越大，抽水流量越大。由此可見，在抽水早期，抽水井附近區(qū)域含水層的重力儲水逐漸釋放并補(bǔ)充到抽水井內(nèi)，n1越大，湍流程度越大，導(dǎo)致非承壓含水層儲水的釋放越快。隨著抽水的繼續(xù)，N 區(qū)重力儲水的釋放逐漸趨于穩(wěn)定，抽水井中的水主要通過D 區(qū)的水流進(jìn)行補(bǔ)給，導(dǎo)致Q-t曲線逐漸趨同。

圖7 不同非達(dá)西湍流系數(shù)對泵送流量的影響Fig.7 Effects of different non-Darcy coefficients on pumping flow

圖7b 體現(xiàn)了t=1 d 時不同徑向距離的非達(dá)西湍流的發(fā)展特征。結(jié)果表明，固定r下h隨著n1的增大而增大，在r較小時不同曲線對應(yīng)的h之間的偏差較大，隨著r的增大，h偏差逐漸減小。非達(dá)西系數(shù)與水頭之間存在正相關(guān)關(guān)系，在r較小時h-r曲線上升的速度較快，隨著r的增大，曲線斜率逐漸減小。n1越大，表示水力梯度越大，對應(yīng)流量越大。在抽水井附近，水頭波動較大，但沒有出現(xiàn)不合理的突變，證明了所提有限差分解的適用性。

圖8 展示了徑向距離對斷面流量的影響，其中圖8a 和圖8b 分別用于表示抽水前期 (t=0.1 d) 和抽水后期 (t=100.0 d) 的Q1-r曲線(Q1為斷面流量)，其水力參數(shù)見表1 的案例D。對比圖8a 和圖8b 可以發(fā)現(xiàn)，t=0.1 d 時Q1隨r的增加而快速下降，其Q1較t=100.0 d 時更大。隨著抽水時間的增加，R逐漸變小，該結(jié)論符合圖4 中R-t曲線的趨勢。Q1隨著r的增大不斷減小，在N 區(qū)部分Q1曲線隨距離下降的速度較快，并在經(jīng)過轉(zhuǎn)換界面R后出現(xiàn)突變，在D 區(qū)Q1曲線隨距離下降的速率明顯變慢。從理論上來說，定水頭抽水過程中，在井壁處存在Q1變化的最大值，隨著r的增大，距離抽水井越遠(yuǎn)，對抽水井補(bǔ)給量越小，斷面流量越小。在抽水前中期，含水層中距離抽水井較遠(yuǎn)的部分無法及時對抽水井附近進(jìn)行充足的補(bǔ)給，抽水井附近區(qū)域含水層的重力儲水會逐漸釋放并補(bǔ)充到抽水井內(nèi)。隨著抽水時間的增長，在抽水后期，遠(yuǎn)處的水流逐漸補(bǔ)充到抽水井附近，抽水井周邊區(qū)域的水位和斷面流量的下降速度逐漸變小。相比于N 區(qū)，D 區(qū)水力梯度更小，其比流量隨徑向距離的下降速率更慢，其Q1-r曲線的轉(zhuǎn)折點(diǎn)出現(xiàn)在r=R處。

圖8 不同時間下徑向距離對斷面流量的影響Fig.8 Influence of radial distance on cross-section flow at different time

5 結(jié)論

a.所提出的有限差分解可用于在由定水頭抽水試驗(yàn)引起的兩區(qū)流、純達(dá)西流和純非達(dá)西流的情況下進(jìn)行降深模擬。在兩區(qū)流動的情況下，所提出的解也可用于評價抽水流量和非達(dá)西區(qū)域的動態(tài)發(fā)展。

b.忽略非達(dá)西區(qū)域的影響可能會給兩區(qū)流動的水頭模型帶來誤差。有限非達(dá)西區(qū)域的湍流會分別導(dǎo)致兩區(qū)流中水頭較純非達(dá)西流和純達(dá)西流的水頭偏大和偏小，且隨抽水時間的增加逐漸變大。

c.利用Izbash 方程，發(fā)現(xiàn)瞬態(tài)抽水速率-時間變化關(guān)系可以由抽水井內(nèi)水頭hw和非達(dá)西滯流系數(shù)n1顯著控制。在抽水過程中，減小hw或增大n1可以提高抽水速率，這種影響隨著抽水時間的增加會不斷減弱并在抽水結(jié)束時消失。

d.在抽水過程中，Q1隨著r的增大不斷減小，并在r=R處出現(xiàn)轉(zhuǎn)折點(diǎn)，具體表現(xiàn)為：在N 區(qū)中Q1曲線隨距離下降的速度較快，在經(jīng)過轉(zhuǎn)換界面R處后出現(xiàn)突變，在D 區(qū)中Q1曲線隨距離下降的速度明顯變慢。