矩形邊界越流承壓含水層中非完整井穩(wěn)定流解析解

孫前林, 譚衛(wèi)佳, 徐蓓藝, 王旭東

(南京工業(yè)大學(xué) 巖土工程研究所, 南京 211816)

0 引 言

受含水層各向異性、邊界條件和抽水井設(shè)置方式等因素的影響,抽水引起的含水層中地下水流動特征差異明顯．1863年,Dupuit[1]依據(jù)定水頭圓島模型,提出了均質(zhì)、各向同性含水層中完整井穩(wěn)定流公式．由于理想的圓島模型缺乏工程應(yīng)用的適應(yīng)性,1870年,Thiem[2]引入影響半徑概念,建立了圓形邊界承壓含水層完整井穩(wěn)定流計算Thiem模型,提出了承壓含水層完整井地下水二維穩(wěn)定流解析解．在此基礎(chǔ)上,Hantush和Jacob[3-5]進一步提出了考慮越流的承壓含水層完整井降深解析解．針對非完整井抽水的地下水流動,Hantush[6-7]研究了側(cè)向無限越流承壓含水層非完整井地下水模型,求得了井壁流量均勻分布的非完整井流解析解．為了描述圓形邊界含水層中的地下水三維流動,Javandel[8]、Chen等[9]利用有限Fourier余弦變換,在極坐標(biāo)系下分別推導(dǎo)了各向同性、正交各向異性承壓含水層非完整井抽水的三維穩(wěn)定流解析解,豐富了承壓含水層井流理論．

上述研究成果均基于極坐標(biāo)系下地下水流動微分方程、圓形邊界、中心布井的地下水計算模型,較好地反映了單井抽水的軸對稱井流特征,但極坐標(biāo)系下地下水流動解析解難以適應(yīng)矩形邊界含水層、多井、布井方式多樣性等地下水流動問題的計算分析,在工程實踐應(yīng)用中存在局限性．為此,Daly和Morel-Seytoux[10]建立了直角坐標(biāo)系下非均質(zhì)各向同性承壓含水層地下水穩(wěn)定流計算模型,利用二維有限Fourier變換,推導(dǎo)了承壓含水層混合邊界條件下無井流的地下水穩(wěn)定流水頭解析解．在此基礎(chǔ)上,Yeo等[11]提出了矩形邊界承壓含水層多個完整井抽(注)水的穩(wěn)定流降深解析解．針對均質(zhì)各向同性承壓含水層矩形邊界的多樣性,王旭東等[12]利用有限Fourier變換分別推導(dǎo)了直角坐標(biāo)系下定水頭邊界和隔水邊界條件下完整井抽水的非穩(wěn)定流降深解析解．Chan等[13]提出了矩形承壓含水層在六類邊界條件下的完整井非穩(wěn)定流降深解析解,并通過算例與改進反映法進行了對比驗證[14],表明了有限Fourier變換法求解矩形邊界承壓含水層非穩(wěn)定井流模型的可行性．

對于矩形邊界承壓含水層中的地下水井流問題,已有的研究成果在綜合考慮含水層正交各向異性、非完整井、越流等因素對地下水流動的影響上仍然存在局限性．為此,本文通過構(gòu)建直角坐標(biāo)系下考慮含水層正交各向異性、非完整井、抽水井位置、越流等因素的地下水流動數(shù)學(xué)模型,利用有限Fourier余弦變換和有限復(fù)合Fourier變換,推導(dǎo)了不同類型邊界條件下越流承壓含水層中非完整井三維穩(wěn)定流解析解,為合理開展矩形邊界含水層的減壓降水設(shè)計提供計算方法．

1 計算模型建立

直角坐標(biāo)系下,承壓含水層均質(zhì)、正交各向異性、等厚且水平,初始地下水的水頭水平．地下水流動服從Darcy定律,任意位置非完整井定流量抽水,水流沿井壁均勻進水,考慮含水層之間的越流補給,且越流強度與降深線性相關(guān)．圖1為定水頭邊界條件下越流承壓含水層非完整井的計算模型．在定水頭邊界條件下,以降深表示的越流承壓含水層非完整井穩(wěn)定流數(shù)學(xué)模型為

(1)

s|x=0=s|x=L=0,

(2)

s|y=0=s|y=B=0,

(3)

(4)

(5)

式中,s為降深,m;Kx,Ky,Kz分別為承壓含水層橫向、縱向和豎直向的滲透系數(shù),m/d;K′z為弱透水層豎直向的滲透系數(shù),m/d;L,B,M分別為承壓含水層的長度、寬度和厚度,m;M′為弱透水層的厚度,m;Q為單井抽水量,m3/d;d為抽水井濾管到含水層頂板的距離,m;lw為抽水井濾管長度,m;(xw,yw)為抽水井的位置坐標(biāo),m;δ(x-xw),δ(y-yw)為Dirac函數(shù)[15-16],用以描述抽水井的位置．

圖1 矩形邊界越流承壓含水層非完整井計算模型Fig. 1 The calculation model for the partially penetrating well in a rectangular leaky-confined aquifer

2 計算模型求解

矩形邊界越流承壓含水層中非完整井抽水引起的地下水三維流動,地下水流動控制方程(1)為關(guān)于空間變量x,y,z的二階非齊次線性偏微分方程,源匯項中流量式(5)為分段非連續(xù)函數(shù)．首先利用有限Fourier變換法[17]將二階非齊次線性偏微分方程變換為常系數(shù)方程,求得有限Fourier變換域中降深函數(shù),然后對降深函數(shù)進行逆變換得到實際空間域中的降深解析解[18]．

2.1 有限Fourier變換域降深函數(shù)

如式(5)所示,非完整井的實管段封閉,濾管段定流量Q均勻流入．根據(jù)有限復(fù)合Fourier變換的適用條件,首先需要對數(shù)學(xué)模型做關(guān)于變量z的有限Fourier余弦變換[19-20],將數(shù)學(xué)模型轉(zhuǎn)化為關(guān)于空間變量x,y的一維有限Fourier變換域中的數(shù)學(xué)模型．

關(guān)于空間變量z的有限Fourier余弦變換的表達式如下[21-22]:

(6)

由于承壓含水層的頂?shù)装寰鶠楦羲吔?如式(4)所示,z方向的特征函數(shù)χk(λkz)取為cos(kπz/M),相應(yīng)的特征值為λk=kπ/M．

對地下水流動控制方程(1)和流動邊界條件(2)、(3)進行有限Fourier余弦變換(變換過程見附錄A),經(jīng)有限Fourier變換后的二維地下水?dāng)?shù)學(xué)模型為

(7)

(8)

(9)

(10)

有限Fourier變換域中二維地下水?dāng)?shù)學(xué)模型中的地下水流動方程(7)和邊界條件(8)、(9)連續(xù)且齊次,可對地下水流動方程(7)進行有限復(fù)合Fourier變換[23-24]．關(guān)于空間變量x,y的有限復(fù)合Fourier變換的表達式如下[25]:

(11)

由于數(shù)學(xué)模型的外邊界均為定水頭邊界,x,y方向特征函數(shù)φm(αmx),ψn(βny)分別取為sin(mπx/L),sin(nπy/B),相應(yīng)的特征值分別為αm=mπ/L,βn=nπ/B．根據(jù)邊界條件(8)、(9)對控制方程(7)進行有限復(fù)合Fourier變換,有限復(fù)合Fourier變換后的降深方程為

(12)

(13)

2.2 實際空間域降深解析解

(14)

式中,K(λk),M(αm),N(βn)分別為空間變量z,x,y方向特征函數(shù)的范數(shù),范數(shù)與邊界條件有關(guān),根據(jù)特征函數(shù)確定,參見附錄B．

(15)

式中

式(15)中右側(cè)第一項表示完整井抽水的降深分布,第二項表示非完整井抽水對降深分布的貢獻．采用疊加原理,式(15)可推廣至多抽注水井問題的求解．

越流承壓含水層矩形邊界條件具有多樣性,不同類型邊界條件下越流承壓含水層非完整井穩(wěn)定流數(shù)學(xué)模型的解析解可參考本小節(jié)內(nèi)容和附錄B的推導(dǎo),相應(yīng)數(shù)學(xué)模型的降深解析解見附錄C．

3 驗證與討論

3.1 解析解退化驗證

為了驗證矩形邊界越流承壓含水層中非完整井穩(wěn)定流降深解析解的正確性, 我們在完整井抽水、 無越流條件下,ξk=0,K′z=0, 式(15)可轉(zhuǎn)化為定水頭邊界條件下水平正交各向異性承壓含水層完整井穩(wěn)定流降深解:

(16)

式(16)與文獻[11]給出的水平正交各向異性承壓含水層完整井穩(wěn)定流降深解析解完全一致．式(16)的各向同性(Kx=Ky)解析解也與文獻[12-13]提出的各向同性承壓含水層完整井非穩(wěn)定流降深長期(t→∞)解析解相一致．

上述分析表明,本文提出的矩形邊界越流承壓含水層非完整井穩(wěn)定流降深解析解,不僅可退化為考慮正交各向異性、完整井等因素的諸多已有解答,而且還能全面考慮承壓含水層的越流、非完整井、任意位置布井等復(fù)雜因素．

3.2 解析解對比驗證

為了驗證不同坐標(biāo)系下降深解析解的差異,求解越流和非完整井問題的正確性,我們利用有限元數(shù)值解和極坐標(biāo)系下相關(guān)解析解[2,5,9]等進行對比驗證．

直角坐標(biāo)系下,定水頭邊界越流承壓含水層非完整井抽水計算模型為L=B=500 m,M=20 m,M′=10 m,Kx=Ky=Kz=5 m/d,K′z=0.02 m/d,d=5 m,lw=10 m,Q=1 000 m3/d, 抽水井位置(xw,yw)=(250 m,250 m)．極坐標(biāo)系下計算區(qū)域半徑R=250 m,抽水井位于圓心．

圖2給出了計算模型的降深．Thiem解[2]和Hantush解[5]分別是極坐標(biāo)系下各向同性含水層完整井和各向同性越流含水層完整井的降深解．對比分析表明,兩種坐標(biāo)系下降深計算結(jié)果具有良好的一致性,當(dāng)含水層有越流補給時,其降深較無越流時?。畼O坐標(biāo)系下降深較直角坐標(biāo)系下降深略小,其原因是圓形邊界條件下地下水容易得到補給所致,而本文解析解與相同計算模型降深數(shù)值解的一致性較好,說明了兩種坐標(biāo)系下計算區(qū)域形狀均會對地下水流動產(chǎn)生影響．

由圖3給出的含水層頂面處的降深分布可知,降深解析解式(15)能較好地反映非完整井抽水時的降深分布,體現(xiàn)了抽水井完整性對含水層中地下水流動的影響．在距離抽水井一倍含水層厚度范圍內(nèi),完整井和非完整井抽水時降深存在明顯差異．對比分析發(fā)現(xiàn),降深分布不僅與非完整井濾管長度有關(guān),還與濾管位置d有關(guān)．在相同抽水量條件下,濾管越長,降深愈小;濾管越靠近含水層頂板,降深愈大．因此在矩形邊界承壓含水層地下水流動計算中,有必要考慮抽水井的完整性和濾管設(shè)置位置對地下水流動的影響．

圖2 降深分布曲線(越流與無越流)圖3 降深分布曲線(完整井與非完整井)Fig. 2 Drawdown distribution curves(leakage and non-leakage)Fig. 3 Drawdown distribution curves(fully penetrating well and partially penetrating well)

3.3 降深解析解計算項數(shù)取值

由于式(15)降深解析解采用級數(shù)表示,計算項數(shù)直接影響降深的計算精度．為了分析m,n取值對降深計算精度的影響,定義降深的相對誤差Δ為

(17)

式中,sSN和s300分別是計算項數(shù)m=n=TSN和計算項數(shù)m=n=300時的降深值,m．

圖4、圖5分別給出了完整井、非完整井的降深和相對誤差與計算項數(shù)TSN的關(guān)系．隨著計算項數(shù)的增加,不同位置處的降深均漸趨穩(wěn)定,相對誤差呈現(xiàn)周期性減小趨勢．計算結(jié)果顯示,相同精度條件下,完整井降深的收斂速度比非完整井降深的收斂速度要快,離抽水井距離越近,降深的收斂速度越慢．從地下水的流動特征來看,由于抽水井附近降深大,地下水三維流動明顯．由此可見,在計算這些區(qū)域降深時,需要更多的計算項數(shù)才能滿足精度要求．若以相對誤差±2%的精度確定計算項數(shù),對于完整井,計算點位于一倍含水層厚度以外取100項,一倍含水層厚度以內(nèi)取200項．對于非完整井,需通過增加計算項數(shù)以滿足計算精度要求,計算項數(shù)可分別取150項和250項．

圖4 降深及相對誤差隨計算項數(shù)的變化(完整井)圖5 降深及相對誤差隨計算項數(shù)的變化(非完整井)Fig. 4 Variations of drawdowns and relative errors with the Fig. 5 Variations of drawdowns and relative errors with the number of calculation items (fully penetrating well) number of calculation items (partially penetrating well)

3.4 邊界條件對降深的影響

圖6給出了四周定水頭邊界(附錄C中case 1)、定水頭邊界與隔水邊界對稱組合條件(附錄C中case 3)下的降深分布曲線．兩種情況下x方向均為定水頭邊界,降深分布基本一致,如圖6中右側(cè)降深曲線所示．受y方向邊界條件影響(圖6中左側(cè)降深曲線顯示),在隔水邊界上產(chǎn)生降深,隔水邊界條件下的降深較定水頭邊界條件下的降深要大,其正確性在文獻[10]中得到驗證．地下水降深的分布特征反映了含水層邊界條件對地下水流動的影響,計算結(jié)果與有限元數(shù)值解的一致性良好,也進一步說明了降深解析解的正確性．

圖6 降深分布曲線Fig. 6 Drawdown distribution curves

圖7 承壓含水層非完整井水力梯度(剖面)Fig. 7 Gradients of the partially penetrating well in a confined aquifer(profile)

3.5 非完整井及含水層正交各向異性對地下水滲流場的影響

3.6 抽水井位置對地下水滲流場的影響

圖8為矩形邊界越流承壓含水層完整井抽水的降深等值線及水力梯度．與極坐標(biāo)系下抽水井只能布置在模型中心、地下水呈現(xiàn)軸對稱流動相比,直角坐標(biāo)系下抽水井的布置更加靈活,能夠充分反映抽水井位置、含水層的定水頭邊界和隔水邊界對地下水流動的影響,降深解析解式(15)提高了工程實踐的適用性．

圖8 降深等值線及水力梯度(平面)Fig. 8 Drawdown contours and gradients(plane)

4 工程案例分析

4.1 工程概況

天津某地鐵車站場地進行了抽水試驗[26-27], 抽水目標(biāo)承壓含水層厚度M為11 m, 滲透系數(shù)K為 0.59 m/d;上部弱透水層厚度M′Ⅰ為7.2 m,下部弱透水層厚度M′Ⅱ為3.9 m,上、下弱透水層的滲透系數(shù)K′zⅠ=K′zⅡ=0.001 m/d,非完整井濾管長度lw為10 m,抽水量Q為24 m3/d．根據(jù)抽水試驗確定的影響半徑,地下水流動計算模型平面尺寸取為200 m×200 m,模型外邊界為定水頭邊界,非完整抽水井位于模型中心．抽水井和觀測井的平面布置如圖9所示,概化后的越流承壓含水層非完整井穩(wěn)定流計算模型如圖10所示．

4.2 降深對比分析

根據(jù)線性偏微分方程疊加原理,當(dāng)承壓含水層上、下均為弱透水層時,定水頭邊界條件下越流承壓含水層非完整井穩(wěn)定流降深解析解為

(18)

式中

θmn=sin(mπxw/L)sin(nπyw/B),

利用降深解析解式(18)對抽水試驗進行了降深計算,計算結(jié)果如圖11所示,降深計算值與實測值[26-27]在數(shù)值大小和分布規(guī)律上都吻合良好,計算精度能夠滿足工程設(shè)計要求,也表明了矩形邊界越流承壓含水層非完整井計算模型在工程實踐中的適用性．

圖9 抽水井和觀測井布置(單位: m)Fig. 9 The pumping well and the observation well layout(unit: m)

圖10 地下水計算模型(單位: m)圖11 降深計算值與實測值對比Fig. 10 The groundwater calculation model(unit: m)Fig. 11 Comparison of calculated and measured values of the drawdown

5 結(jié) 論

1) 本文考慮了矩形邊界越流承壓含水層的正交各向異性、越流、非完整井、邊界條件等因素對地下水流動的影響,利用有限Fourier變換和有限復(fù)合Fourier變換,提出了矩形邊界越流承壓含水層中非完整井穩(wěn)定流降深求解方法,給出了不同類型邊界條件下越流承壓含水層中非完整井穩(wěn)定流降深解析解．

2) 通過降深解析解的退化驗證,與極坐標(biāo)系下降深解析解的對比驗證,闡明了本文降深解析解的正確性．為保證降深解析解的計算精度,根據(jù)抽水井的完整性、地下水流動特征以及計算點的位置,給出了降深解析解計算項數(shù)的合理取值．

3) 矩形邊界越流承壓含水層中非完整井穩(wěn)定流降深解析解能夠考慮正交各向異性、越流、非完整井、任意井位等因素,展現(xiàn)了良好的工程適用性,可以為工程實踐中合理開展承壓含水層非完整井減壓降水設(shè)計提供計算方法,并能推廣至抽注水多井問題的地下水流動問題分析中．

致謝本文作者衷心感謝上海隧道工程有限公司科研項目(2021-SK-21)對本文的資助．

附錄A 地下水流動控制方程的有限Fourier變換

根據(jù)有限Fourier余弦變換公式(6),對地下水流動控制方程(1)的每一項進行變換,由于z方向的兩端為隔水邊界,特征函數(shù)為cos(kπz/M),有限Fourier余弦變換的結(jié)果為

(A1)

(A2)

(A3)

(A4)

Qδ(x-xw)δ(y-yw)ξk,

(A5)

式中

將每一項的變換式(A1)—(A5)代入原地下水流動控制方程(1),即可得到有限Fourier余弦變換后的地下水流動控制方程(7)．

附錄B 特征函數(shù)、特征值、范數(shù)以及變換參數(shù)

特征函數(shù)、 特征值、 范數(shù)與數(shù)學(xué)模型邊界條件有關(guān), 不同類型邊界條件下特征函數(shù)、 特征值、 范數(shù)、 變換參數(shù)的取值見表B1．

表B1 不同類型邊界條件下的特征函數(shù)、特征值、范數(shù)以及變換參數(shù)

附錄C 不同類型邊界條件下降深解析解

表C1 降深解析解(case 1～3)

表C2 降深解析解(case 4)

表C3 降深解析解(case 5,6)