肖良，柳岳輝，陳立華，于江偉

(1.廣西大學(xué) 土木建筑工程學(xué)院，廣西 南寧 530004；2.廣西防災(zāi)減災(zāi)與工程安全重點實驗室，廣西 南寧 530004；3.中天建設(shè)集團有限公司廣西分公司，廣西 南寧 530004)

0 引言

承壓含水層抽水井流模型是地下水滲流理論的重要內(nèi)容，抽水試驗作為直接獲取研究區(qū)水文地質(zhì)參數(shù)的有效途徑，一直被學(xué)者和工程師廣泛采用。經(jīng)典泰斯模型給出了水位降深(s)與其自變量(時間t和井距r)的解析解，從而能粗略估計研究地水文地質(zhì)參數(shù)。該模型主要適用于完整井附近的達西滲流，即抽水井滲流場的雷諾系數(shù)介于1至10之間。根據(jù)室內(nèi)試驗和工程實踐發(fā)現(xiàn)，抽水井附近滲流容易出現(xiàn)非達西現(xiàn)象，即抽水井滲流場的雷諾系數(shù)小于1或者大于10，因此基于達西定律的泰斯抽水井模型在非達西滲流應(yīng)用上具有一定的局限性。此外，在發(fā)展中國家，大口井(指抽水條件下井中的儲水效應(yīng)不能忽略的抽水井)已被廣泛應(yīng)用于地下水開采中。大口井的直徑一般在2～8 m，具有進水斷面大、構(gòu)造簡單、取材容易、使用年限長等特點[1]。如何合理開發(fā)利用地下水、正確評價地下水資源、防止地下水污染等問題，已成為學(xué)者們的研究熱點[2-4]。因此在考慮井儲效應(yīng)的基礎(chǔ)上，開展完整井附近非達西滲流模型研究對于正確認識抽水試驗具有一定理論支撐作用。

當抽水速率足夠大時，在孔隙較大的多孔介質(zhì)含水層或含有較多裂隙的含水層中，地下水滲流速度與水力梯度之間不再滿足線性關(guān)系，對于此類不滿足達西線性定律的地下滲流稱之為非達西流[5]。由于達西定律不能滿足所有地下水滲流情況，地下水流非線性運動機理的研究自20世紀以來便成為水文地質(zhì)領(lǐng)域的研究熱點，廣大學(xué)者對不同介質(zhì)下的非達西流進行試驗及理論研究，并提出相應(yīng)的非達西滲流公式。其中Forchheimer于1901年提出的Forchheimer公式和Izbash在1931年提出的Izbash公式在非達西流建模中應(yīng)用最為廣泛，具體選用哪個公式一般取決現(xiàn)場水文地質(zhì)條件[6]。

目前關(guān)于大口井水力學(xué)問題的研究大多基于滲流滿足達西定律而展開[7-14]。對于考慮井儲效應(yīng)下非達西井流建模研究，土耳其學(xué)者SEN[15]運用Forchheimer公式描述非達西流，通過Boltzmann變換得到完整井附近非達西流近似解析公式，并利用特征曲線對公式中的相關(guān)參數(shù)進行特征分析。WEN等[6]采用Laplace變換以及線性化手段求得單一抽水井附近非達西流降深半解析解，通過建立數(shù)值解對該解析解進行對照驗證，發(fā)現(xiàn)線性化方法用于描述非達西流而解得的近似解在抽水前期的計算結(jié)果與數(shù)值解存在一定誤差。陳建華等[16]用滲流理論對大口井地層氣體進行模擬，將井筒內(nèi)視為“達西-非達西”兩區(qū)模型，并通過雷諾數(shù)來判定其分界點，得到井底壓力分布圖。上述研究為大口井附件非達西滲流模擬提供了一定的理論基礎(chǔ)，但是其方法計算過于復(fù)雜，依賴于計算機程序以及特征曲線擬合等，不便于水文地質(zhì)工作者直接使用。

因此，本文擬在無限承壓含水層的等速抽水試驗中，提出考慮井儲效應(yīng)的簡化非達西完整井流模型。筆者采用Izbash公式來描述非線性非達西流，采用Boltzmann變換和線性化近似相結(jié)合的方法推導(dǎo)得其簡化解析模型，并通過繪制降深-時間曲線對解析模型進行參數(shù)分析，為利用大口井開采的相應(yīng)地下工程實施提供理論依據(jù)及參考。

1 數(shù)學(xué)模型與方法

承壓含水層大口井附近水動力模型與文獻[13]中物理模型一致。如圖1所示，假設(shè)含水層中各向同性且均質(zhì)，含水層厚度不變且水平無限長，且抽水速率始終保持恒定。在基于泰斯的假設(shè)條件之外，考慮大口井自身井儲的影響，建立起數(shù)學(xué)模型如下文。

圖1 無限承壓含水層大口徑井滲流模型

流體連續(xù)性方程：

(1)

式中，r為距抽水井井中心的距離,m；t表示為抽水時間,s；q(r,t)為單位水流通量,m/s；S為儲水系數(shù)，用來表征含水層儲水能力的無量綱系數(shù);b為含水層厚度,m;s(r,t)表示為水位降深,m。

抽水前天然狀態(tài)下水力梯度為零，水平面保持一致且降深為零。模型初始條件為

s(r,0)=0。

(2)

在含水層無限遠處的邊界條件為

s(r→∞,t)=0。

(3)

考慮到抽水井自身井筒的儲水，根據(jù)質(zhì)量守恒定律，井內(nèi)部的定解關(guān)系可描述為

(4)

式中，Q為抽水速率，m3/s;rw為抽水井濾網(wǎng)有效半徑,m；rc為井套半徑,m3；sw(t)為抽水井井內(nèi)降深,m。

式(4)說明在大口井中，由于井徑較大，其儲水能力較強，抽水過程中井內(nèi)流量來自于兩方面：一是井筒內(nèi)儲水，即井儲效應(yīng)；二是含水層中對于井內(nèi)的補給流量。

運用Izbash公式來描述非達西流，其水流通量與水力梯度關(guān)系為

(5)

式中，n、k為經(jīng)驗系數(shù)，k為準水力傳導(dǎo)系數(shù),(ms-1)n，是用來描述介質(zhì)傳輸水流快慢的一個特征參數(shù);n是用來描述非達西效應(yīng)的無量綱系數(shù)，一般取值為1～2。當n等于1時，式(5)即為達西定律。

由式(5)求得q對于r的偏導(dǎo)數(shù)：

(6)

將式(6)代入式(1)可得

(7)

為除去式(7)中的非線性項，做線性化近似假設(shè)如下：

(8)

式(8)表示在抽水穩(wěn)定時，單位時間內(nèi)含水層中離井中心距離為r的任意過水斷面的流量均為Q，即2πbrq=Q。文獻[6]和文獻[17]研究已表明該近似替換引起的模型計算誤差非常小，可用于非達西流建模過程。

將式(8)代入式(7)可得非達西流線性近似為

(9)

將式(5)代入式(4)可得

(10)

Boltzmann變化在求解關(guān)于偏微分方程問題上是一種常用且便利的技術(shù)手段，因此采用Boltzmann變化將上述偏微分方程轉(zhuǎn)換為常微分方程[18]，引入?yún)⒘縱(η)和η如下所示：

(11)

(12)

式中，v(η)為s和t的相似替代項；η為r和t的相似替代項;β和α為常數(shù)。

通過應(yīng)用式(11)、式(12)，將式(7)中三個不同偏微分項可統(tǒng)一替換為關(guān)于η的方程，并如下所示：

(13)

(14)

(15)

假設(shè)β=2α以及γ=0，式(9)可表述為

(16)

由式(16)解得

(17)

將上述變化代入其初始條件及邊界條件中，式(2)、式(3)、式(10)可改寫為

v(η→∞)=0，

(18)

v(0)=0，

(19)

(20)

由式(12)、式(20)解得

(21)

從而求得

(22)

聯(lián)立式(17)、式(22)可解得

(23)

將式(23)進行化簡及積分求解可得

(24)

由式(18)、式(19)可知

(25)

因此聯(lián)立式(24)、式(25)解得

(26)

(27)

將式(27)簡化如下：

(28)

2 計算結(jié)果分析

2.1 模型準確性分析

本文將上述模型解與文獻[1]研究成果進行比較驗證提出模型的準確性。假設(shè)無井儲效應(yīng)，即取rc、rw為0，可得M=1。因此，不考慮井儲效應(yīng)前提下，由式(27)可將承壓無限含水層等速非達西完整大口井流模型可簡化為

(29)

由此可得，式(29)與文獻[17]中式(11)表達形式完全一致，由此可證明本文推導(dǎo)模型的準確性。

2.2 井儲效應(yīng)影響分析

為研究井儲效應(yīng)對抽水降深的影響，本文假設(shè)在一承壓無限含水層中，令抽水速率Q=0.108 m3/s，觀測井r取值為10 m，井筒半徑rc為1 m，濾網(wǎng)半徑取rw為0.5 m，含水層厚度b設(shè)為50 m，儲水系數(shù)S取值為0.003，將非達西系數(shù)n設(shè)為1.4，k取5×10-6(m/s)1.4，分別采用式(27)、式(29)計算降深-時間曲線，結(jié)果如圖2所示。

圖2 井儲效應(yīng)對于降深影響

由圖2可知，井儲效應(yīng)的存在導(dǎo)致其降深-時間曲線負偏離無井儲效應(yīng)的降深-時間曲線，即因井儲作用導(dǎo)致其在某一個抽水時間點的降深小于無井儲作用下的降深。抽水過程中井儲效應(yīng)對于降深的影響一般發(fā)生在抽水前期，而隨著抽水過程的進行，抽水井自身的儲水將逐漸被抽干，降深變化主要受含水層補給控制，因此井儲效應(yīng)的影響隨著抽水的進行逐漸消失。

2.3 準水力傳導(dǎo)系數(shù)k對于降深的影響

基于上述解析模型，本文將對準水力傳導(dǎo)系數(shù)k和非達西系數(shù)n對特征曲線的影響進行分析。假設(shè)抽水泵抽水流量Q取值為0.008 m3/s，觀測井r取值為10 m，井筒半徑rc為1 m，濾網(wǎng)半徑取rw為0.3 m，含水層厚度b設(shè)為50 m，儲水系數(shù)S取值為0.003，將非達西系數(shù)n設(shè)為1.4，編寫相關(guān)MATLAB程序，就k在取值為5×10-5、5×10-6、5×10-7(m/s)1.4三種情況下分別進行分析討論。

如圖3所示為不同k值情況下的觀測井水位降深-時間曲線。結(jié)果表明，隨著抽水過程的進行，承壓含水層中水位降深均隨時間增加而變大。在抽水前期，水位降深變化與準水力傳導(dǎo)系數(shù)k之間呈正相關(guān)，而在中后期，呈負相關(guān)。一般而言，k值越小含水層中介質(zhì)輸送水流速度越慢。在保持抽水泵抽水速率恒定時，由于觀測井周圍含水層對于井內(nèi)的補給速率不同會導(dǎo)致井內(nèi)降深不同。因而，在抽水中后期，k值較大時，觀測井內(nèi)水位降深相對于k值較小時降深相比較小。

圖3 不同k值下水位降深-時間曲線

2.4 非達西系數(shù)n對于降深影響

同理，取假設(shè)抽水泵抽水流量Q取值為0.008 m3/s，觀測井r取值為10 m，井筒半徑rc為1 m，濾網(wǎng)半徑取rw為0.3 m，含水層厚度b設(shè)為50 m，儲水系數(shù)S取值為0.003，準水力傳導(dǎo)系數(shù)k取值為5×10-6(m/s)n，n取值分別為1.2、1.4、1.6三種情況進行討論。

不同n取值下降深-時間曲線模擬結(jié)果如圖4所示。由此可知，非達西系數(shù)n與準水力傳導(dǎo)系數(shù)k在對水位降深影響的效應(yīng)相同；即在前期表現(xiàn)為非達西系數(shù)n與準水力傳導(dǎo)系數(shù)k之間呈現(xiàn)正相關(guān)，而在抽水中后期呈負相關(guān)。在水力學(xué)中，水頭指的是含水層中單位液體所具有的機械能，包括位置水頭、壓強水頭、流速水頭。隨著抽水的進行，井內(nèi)發(fā)生能量轉(zhuǎn)換，對于觀測井而言，主要表現(xiàn)在井內(nèi)降深(即水頭)的增大。由于流體的連續(xù)性，抽水井與含水層之間進行著能量的轉(zhuǎn)換。因此，井內(nèi)水位降深變化，側(cè)面反映著流體能量的交換。一般來說，當n值越大(n>1)，流體的紊動程度越大，流體由于慣性力占據(jù)主導(dǎo)，流體內(nèi)包含的勢能較大。紊流相對于的層流(n=1的達西流)，短期供水能力更強。因此，在同一抽水時刻，非達西系數(shù)n越大的流體，將自身勢能轉(zhuǎn)換補給井內(nèi)水頭能力越強，從而水位降深變化反而較小。此外，也可解釋另一規(guī)律：當非達西系數(shù)n的取值越大，非達西流趨于穩(wěn)定的時間越早。

圖4 不同n值下降深―時間曲線

2.5 參數(shù)敏感性分析

敏感性分析是一種量化輸入量與輸出量之間的不確定性關(guān)系的手段，參考文獻[19]提出的標準化公式，本文將作為標準化的敏感度自變量定義為

(30)

(31)

式中，ΔPj一般為10-2×Pj。運用式(31)對Izbash公式中兩個常數(shù)(n和k)進行參數(shù)敏感性分析(圖5、圖6)，選取3.2小節(jié)中同一組參數(shù)分別做出k=5×10-5(m/s)1.4、n=1.1的NSC-時間曲線圖。

圖5 n值敏感性分析

圖6 k值敏感性分析

由圖5、圖6可知，Izbash公式中的兩個經(jīng)驗系數(shù)中非達西系數(shù)n相對于準水力傳導(dǎo)系數(shù)k而言的敏感性要強得多。從Izbash公式的結(jié)構(gòu)來看，n作為反映水流非達西效應(yīng)強弱的主要系數(shù)，相對于k而言，敏感性更強。

4 結(jié)論

本文提出了一種簡化數(shù)學(xué)解析模型用于求解非達西效應(yīng)下承壓無限含水層大口井附近的水位降深，通過對所得到的的簡化解析解進行分析，得到如下結(jié)論：

① 井儲效應(yīng)僅出現(xiàn)在抽水前期，隨著井筒內(nèi)部儲水的消耗，水位降深的變化僅來自于含水層的補給。

② 本文通過對Izbash公式的兩個參數(shù)進行分析，抽水前期水位降深與準水力傳導(dǎo)系數(shù)k和非達西系數(shù)n為正相關(guān)，而在抽水后期呈負相關(guān)。且當n的取值越大，非達西流趨于穩(wěn)定的時間越早。

③ 通過對Izbash公式中的兩個參數(shù)敏感性分析，相較于準水力傳導(dǎo)系數(shù)k，非達西系數(shù)n敏感性更強。研究結(jié)果將對水文地質(zhì)參數(shù)的獲取、地下水資源評估等工作提供了一定的理論支撐。