邊坡負(fù)壓排水非穩(wěn)定流研究

雷 怡, 帥飛翔, 孫紅月, 張澤坤, 熊 超

浙江大學(xué)海洋學(xué)院，浙江 舟山 316021

0 引言

我國東部地區(qū)的滑坡多為降雨誘發(fā)型，有效的邊坡排水措施是此類滑坡防治的關(guān)鍵[1-5]。目前，邊坡虹吸方法已經(jīng)被廣泛應(yīng)用在工程實(shí)踐當(dāng)中[6-14]。自啟動(dòng)負(fù)壓排水技術(shù)是基于邊坡虹吸排水的一種改進(jìn)技術(shù)，相比虹吸排水而言，其具有自啟動(dòng)、流速快、不易淤堵的優(yōu)勢[15]。自啟動(dòng)負(fù)壓排水方法的工作原理是：當(dāng)坡體內(nèi)地下水位上升引起排水管進(jìn)水口(圖1中點(diǎn)B)的總水頭大于鉆孔的孔口(圖1中點(diǎn)A)高程時(shí)，透水鉆孔段內(nèi)的地下水就會(huì)在水頭差的作用下從虹吸排水管自然流出，利用虹吸作用在空腔內(nèi)形成負(fù)壓，坡體內(nèi)的地下水會(huì)加速流向空腔內(nèi);當(dāng)透水鉆孔段空腔內(nèi)的地下水排干時(shí)，排水管的進(jìn)水口會(huì)進(jìn)入空氣，排水管的虹吸作用消失，一次排水過程結(jié)束; 隨著坡體地下水入滲充滿鉆孔空腔，引起排水管進(jìn)水口的水頭再次大于鉆孔的孔口高程時(shí)，排水過程再次進(jìn)行。

A.鉆孔的孔口；B.排水管進(jìn)水口。

Dupuit于1863年運(yùn)用常微分方程描述水流連續(xù)性條件，分析地下水的穩(wěn)定流動(dòng)，提出了潛水和承壓水井的流量公式，被稱為Dupuit模型[6]。地下水動(dòng)力學(xué)中關(guān)于潛水運(yùn)動(dòng)的研究[6]始興于此，包括：Forchheimer運(yùn)用Dupuit模型，于1901年建立了潛水面穩(wěn)定分布的偏微分方程，被稱為Dupuit-Forchheimer方程；Bossinesq于1904年進(jìn)一步考慮潛水波動(dòng)過程的質(zhì)量守恒特征，提出了著名的Boussinesq方程，其也是第一個(gè)描述地下水非穩(wěn)定運(yùn)動(dòng)的偏微分方程；Thies于1935年借用熱傳導(dǎo)理論，把滲流場類比為溫度場，流量類比為熱流量、水頭類比為溫度，推導(dǎo)得到的井流公式稱為Thies公式。隨著Theis公式的提出，國內(nèi)外的學(xué)者開始對抽水引起的地下水非穩(wěn)定流動(dòng)進(jìn)行深入的研究工作[6],如：Hantush 等考慮越流條件下建立了非穩(wěn)定流微分方程，并得到了相應(yīng)的解析解；Boulton考慮井流中的滯后給水效應(yīng)，得到了相應(yīng)條件下的非穩(wěn)定井流解析解。

本文在前人非穩(wěn)定井流理論的基礎(chǔ)上，針對自啟動(dòng)負(fù)壓排水技術(shù)的特點(diǎn)，提出了相應(yīng)的計(jì)算模型，模型中假定潛水層均質(zhì)、各向同性，考慮自啟動(dòng)負(fù)壓排水源匯項(xiàng)為點(diǎn)源的情況，引入狄拉克函數(shù)，采用二維傅里葉正逆變換求解負(fù)壓排水模型的定解問題，得到負(fù)壓排水滲流場解析解；再使用Geo-Studio 2012有限元分析軟件建立自啟動(dòng)負(fù)壓排水方法模型，在SEEP/W滲流模塊對其滲流場進(jìn)行數(shù)值模擬，得到二維邊坡剖面模型滲流場分布；最后通過對比數(shù)值模擬，得到的邊坡浸潤線與Dupuit假設(shè)簡化條件下解析解下的浸潤線，用以分析模型的有效性。

1 數(shù)學(xué)模型的建立

1.1 問題描述與基本假設(shè)

自啟動(dòng)負(fù)壓排水理論模型簡圖見圖2。在實(shí)際應(yīng)用中，邊坡負(fù)壓排水技術(shù)的實(shí)施是在滑坡體內(nèi)布設(shè)沿邊坡軸線排狀分布的俯傾孔。在理論模型中，該排俯傾孔可簡化為貫穿y軸的排水溝。由于集水腔體尺度與邊坡尺度相比可忽略不計(jì)，單個(gè)的負(fù)壓排水孔相當(dāng)于坡體里的一個(gè)點(diǎn)匯。無排水干擾時(shí)，可引入Dupuit假定，忽略垂向滲流速度；有負(fù)壓排水干擾時(shí)，由于點(diǎn)匯的存在，垂向滲流速度不可忽略，即不能利用Dupuit假定，故考慮剖面二維流的點(diǎn)源情況。

h1，h2分別為左水頭邊界、右水頭邊界。

1.2 控制方程

式(1)為二維潛水運(yùn)動(dòng)的控制方程[6]：

(1)

式中：x，y分別為水平、豎直方向距離變量，m；K為滲透系數(shù)，m/s；h為水頭，m；W為源匯項(xiàng)，m/s；μ為彈性給水度，無量綱；t為時(shí)間，s。由于其為一非線性的拋物型方程，尚無直接求解方法；因此，目前在地下水滲流的數(shù)值模擬研究中，其常見的做法是引入含水層平均厚度hm，使其線性化。本模型考慮的剖面二維潛水非穩(wěn)定流方程，在有源匯時(shí)，式(1)可轉(zhuǎn)換為式(2)：

(2)

式中：z為垂直方向距離變量，即位置水頭，m;μd為重力給水度。

上文提到，在負(fù)壓排水模型中，集水腔體相對于整個(gè)坡體尺度很小，故簡化為一個(gè)點(diǎn)，此處用狄拉克函數(shù)表征點(diǎn)匯，即上式中源匯項(xiàng)函數(shù)W定義為

W=-Qδ(x-x0)δ(z-z0)。

(3)

式中：Q為源匯量,m3/s；δ為狄拉克函數(shù)符號(hào)；(x0,z0)為匯點(diǎn)坐標(biāo)。

1.3 方程的定解條件

1.3.1 初始條件

描述所研究問題初始時(shí)刻(t=0)研究區(qū)域D內(nèi)各點(diǎn)處水頭分布的情況，稱為初始條件[16]。以水平、豎直、垂直三維流動(dòng)為例，通常表述為

h(x,y,z,t)|t=0=h0(x,y,z)，(x,y,z)∈D。

(4)

式中，h0為t=0條件下的水頭函數(shù)h(x,y,z,t)。

1.3.2 邊界條件

邊坡地下水的流動(dòng)問題主要有3類邊界條件。

1)第一類邊界條件，給定水頭的邊界條件[16]。

已知邊界上水頭分布的邊界條件即屬于此類，以三維流為例,可以表示為

h|B1=hb(x,y,z,t)，(x,y,z)∈B1。

(5)

式中:B1為研究區(qū)D上的第一類邊界；hb為B1上的已知水頭函數(shù)。若邊界上的水頭或水頭函數(shù)不隨時(shí)間改變，則可稱為定水頭邊界，即可表示為

h|B1=hb(x,z)，(x,z)∈B1。

(6)

2)第二類邊界條件，給定流量的邊界條件[16]。

邊界上單寬流量q或滲流速度v已知，或水力坡度已知的邊界條件稱為第二類邊界條件。對于本研究所關(guān)注的剖面二維流則可以表示為

(7)

3)第三類邊界條件，潛水面邊界條件。

地下水的補(bǔ)給來源主要是大氣降水入滲，平原淺埋地區(qū)地下水的蒸發(fā)是其主要排泄形式。因此，潛水面作為地下水入滲補(bǔ)給和蒸發(fā)排泄的邊界，是任何一個(gè)完整地下水排泄系統(tǒng)必須刻畫的極其重要的邊界條件。

如果研究的區(qū)域條件不允許引入Dupuit假設(shè)，則潛水面要作為上邊界條件來刻畫，由于潛水面是移動(dòng)邊界，其形狀和位置未知，因此建立潛水面邊界條件比上文所敘述的其他邊界條件要復(fù)雜得多。

本文采用類似潛水運(yùn)動(dòng)控制方程推導(dǎo)的質(zhì)量守恒法推導(dǎo)(三維流情況)，在潛水層取一尺度為Δx、Δy、Δz的微元(圖3)。

vx|(x,y,z,t)、vx|(x+Δx,y,z,t)分別為水平方向流入和流出微元的地下水流速。

(8)

式中：z為位置水頭，m；p為壓強(qiáng)，Pa；ρ為密度，kg/m3；g為重力加速度，m/s2。

潛水面高程為zwt，考慮三維流，其為x、y、t的函數(shù)，即

zwt=zwt(x,y,t)。

則邊界條件可以描述為，當(dāng)點(diǎn)處于潛水面上，即點(diǎn)的高程z等于zwt時(shí)，水頭函數(shù)值h也等于zwt(當(dāng)基準(zhǔn)取在隔水地板)，即

H=h(x,y,z,t)=zwt(x,y,t)，

(x,y,z)∈Bwt或z=zwt。

(9)

式中:H為潛水面各點(diǎn)水頭；Bwt為潛水面邊界。

由水均衡原理，控制體各側(cè)面流入及流出的結(jié)果將導(dǎo)致zwt的變化，則

[vx|(x,y,z,t)zwt-vx|(x+Δx,y,z,t)zwt]ΔyΔt+

[vy|(x,y,z,t)zwt-vy|(x,y+Δy,z,t)zwt]ΔxΔt=

μd[zwt|(x,y,t+Δt)-zwt|(x,y,t)]ΔyΔx，

z=zwt。

(10)

方程兩端同時(shí)除以ΔyΔxΔt，則得到

(11)

由Δx→0、Δy→0及Δt→0可得

(12)

又有：

(13)

則可將式(12)改寫為式(14)：

(14)

對于有源匯項(xiàng)的情況，式(14)可改寫為式(15)：

(15)

引入Darcy定律，則式(15)可轉(zhuǎn)化為

(16)

而對于本文所討論的土體各向同性、含源匯項(xiàng)的剖面二維流問題，邊界條件方程為

(17)

這里需要特別注意的是，式(10)—(17)是建立與浸潤線函數(shù)zwt=zwt(x,y,t)有關(guān)的方程，而非水頭函數(shù)h=h(x,y,z,t)。垂向上的滲流會(huì)使水頭沿z分布不等，即等水頭線不是垂向上的直線，而是曲線；而垂向上的滲流顯然并不會(huì)引起浸潤線的高程變化。應(yīng)注意浸潤線函數(shù)和水頭函數(shù)的意義與區(qū)別。

2 問題的求解

綜合上文推導(dǎo)的控制方程、初始條件和邊界條件方程，本研究中負(fù)壓排水模型的定解問題可由下列方程組描述：

h(x,z,t)|t=0=h0(x)=

(18)

式中，L為左、右水頭邊界間的水平距離。

h(0,z,t)=h1，h(+∞,z,t)=h2；

(19)

以下通過疊加原理來解決本研究的問題。

如果地下水流控制方程采用線性偏微分方程，則這種方程滿足二階線性偏微分方程的疊加原理。設(shè)Hi是方程

的解，而Hj是方程

的解。令

H(x,z,t)=aiHi(x,z,t)+ajHj(x,z,t)，

(20)

其中ai、aj為常數(shù)，則H(x,z,t)必是方程

(21)

的解。

特別地，如果H0是以下齊次方程

的解，則H(x,z,t)=a0H0(x,z,t)+aiHi(x,z,t)+ajHj(x,z,t)也是方程(21)的解。其中，a0為常數(shù)。

通過將水頭函數(shù)分為初始水頭函數(shù)和變水頭函數(shù)的方法，本文將本模型定解問題分解為兩個(gè)子問題。

子問題一：

(22)

式中：h0為水頭的初始分布函數(shù)；s為變水頭函數(shù)。則總水頭函數(shù)h=h0+s。這樣分離的目的是使子問題二的邊界條件齊次化，大大有利于方程的求解。

子問題二：

s(x,z,t)|t=0=0；

s(∞,z,t)=0；

(23)

對控制方程(22)做x、z的二維傅里葉變換(可依次變換)。對x做傅里葉變換：

對z做傅里葉變換：

接下來求解二維像函數(shù)F2與t的偏微分方程，得

(24)

帶入同樣變換后的初始條件，得C0=0，即

(25)

接下來對F2作逆變換。觀察函數(shù)形式，我們發(fā)現(xiàn)對像函數(shù)積分使用極坐標(biāo)系更為方便。則

y=rsinθ；

ξ=ρsinφ。

則

(26)

將原函數(shù)也用極坐標(biāo)表達(dá)，即

x=rcosθ，y=rsinθ。

則

(27)

若條件允許簡化s，使其關(guān)于匯點(diǎn)(x0,z0)具有圓對稱性，即s(x,z,t)=s(r,θ,t)=s(r,t)，利用貝塞爾函數(shù)關(guān)系，即

(28)

式中，J0(a1)是一個(gè)0階貝塞爾函數(shù)。則式(28)可轉(zhuǎn)換為

(29)

與潛水/承壓水二元結(jié)構(gòu)完整井流的降深公式(34)模型及解析解公式結(jié)構(gòu)類似。

(30)

其中：

2λ1=αtη(1+x2)；

由上文所述疊加原理，總水頭函數(shù)即為水頭初始分布與變水頭函數(shù)的疊加。

3 浸潤線的數(shù)值模擬分析

數(shù)值模擬模型：水平方向長度為65 m，邊坡左側(cè)高35 m，右側(cè)高15 m。土層為均質(zhì)、各向同性的介質(zhì)。對模型進(jìn)行網(wǎng)格劃分，網(wǎng)格密度為1 m×1 m，圖4所示。完成了基礎(chǔ)設(shè)置及網(wǎng)格的劃分步驟之后接著對各處材料進(jìn)行定義，并對其參數(shù)進(jìn)行取定。其中，土體的飽和滲透系數(shù)為8×10-6cm/s。

圖4 邊坡模型及網(wǎng)格劃分示意圖

使用SEEP/W模塊做滲流分析時(shí)，首先要確定邊坡體各層材料的非飽和滲透函數(shù)和水-土特征曲線。通常而言，我們可以通過物理實(shí)驗(yàn)來測量飽和土體的體積含水量，也可以利用統(tǒng)計(jì)方法，通過粒度分布曲線來預(yù)測。而對于非飽和土體滲透系數(shù)的測定方法則較困難，特別是當(dāng)基質(zhì)吸力較大時(shí)[17-23]。本研究中使用的體積含水量函數(shù)及滲透函數(shù)是在利用Geo-Studio軟件中SEEP/W模塊提供的幾個(gè)典型的體積含水率曲線及滲透系數(shù)函數(shù)的基礎(chǔ)上，根據(jù)實(shí)際土體假設(shè)進(jìn)行了適當(dāng)修正，為數(shù)值模擬程序提供與模型更加符合的水-土特征曲線[24-26]。

完成材料及參數(shù)的定義后，開始對邊界條件進(jìn)行定義和完善。左右邊界取第一類邊界條件中的定水頭邊界，左邊界初始水頭為28 m，右邊界初始水頭為8 m，坡面為自由排水邊界，底面為第二類邊界條件，是不透水邊界。值得注意的是，考慮設(shè)置邊界條件以模擬排水孔對邊坡滲流場的影響，在排水結(jié)束后的終態(tài)穩(wěn)定情況下，排水孔壓力水頭為一定值。

當(dāng)有自啟動(dòng)負(fù)壓排水干擾后，邊坡總水頭等值線圖如圖5所示。

箭頭方向代表地下水流速方向，箭頭長度代表該點(diǎn)流速大小。

設(shè)置兩種工況分別計(jì)算了自啟動(dòng)負(fù)壓排水方法下數(shù)值模擬結(jié)果和理論解得到的邊坡浸潤線曲線，并進(jìn)行對比分析，結(jié)果如圖6、7所示。

1)其他邊界條件不變，改變排水孔位置。在左、右水頭邊界為h1=28 m,h2=8 m條件下，對排水孔位置在(34，15)(24，15)(24，10)3種條件分別進(jìn)行數(shù)值模擬。將浸潤線數(shù)據(jù)導(dǎo)出，得到數(shù)值模擬與浸潤線理論解對比圖(圖6)。由圖6可知，3種條件下，數(shù)值模擬結(jié)果與解析求得邊坡浸潤線均吻合良好。值得注意的是，排水孔位置越深，理論解誤差越大。圖6中排水孔位置為(24，10)時(shí)，靠近排水孔區(qū)域誤差最大?？紤]是由于排水孔位置越深，忽略流速垂分量造成的影響越大，故理論解誤差越大。且排水孔越深，實(shí)際排水效果越難達(dá)到理論降深。

圖6 排水孔位置改變浸潤線對比圖

2)排水孔位置不變，改變其他邊界條件。對左、右水頭邊界條件分別為h1=28 m、h2=8 m，h1=25 m、h2=7 m，h1=22 m、h2=6 m三種條件分別進(jìn)行數(shù)值模擬，將浸潤線數(shù)據(jù)導(dǎo)出，得到其與浸潤線解析解對比圖(圖7)。由圖7可見：在排水孔位置固定的情況下，左、右水頭同時(shí)越高，兩曲線之間差距越大，如h1=28 m、h2=8 m條件下；反之，左、右水頭值同時(shí)越低，兩曲線之間差距越小，數(shù)值模擬與理論值吻合越好，如h1=22 m、h2=6 m條件下。這是由于當(dāng)排水孔位置固定，左、右水頭越高，排水孔相對自由水面越深，則相對排水前降深越大，忽略流速垂分量的影響越大，故理論解誤差越大。

圖7 左右水頭邊界改變浸潤線對比圖

4 結(jié)論與建議

1)通過解析法，理論推導(dǎo)并創(chuàng)新地得到了潛水面邊界條件。將剖面二維潛水運(yùn)動(dòng)方程作為控制方程，引入狄拉克函數(shù)表述負(fù)壓排水模型中點(diǎn)匯的數(shù)理特性，結(jié)合排水前坡體滲流場作為初始條件，推得考慮垂向滲流速度，即不利用Dupuit假設(shè)的剖面二維流模型下負(fù)壓排水模型的定解方程組。

2)通過傅里葉正逆變換等數(shù)學(xué)手段推求分析負(fù)壓排水方法干擾后邊坡滲流場，經(jīng)過適當(dāng)簡化后可驗(yàn)證其降深公式與陳崇希解函數(shù)形式一致。為今后的物理實(shí)驗(yàn)、數(shù)值模擬到工程應(yīng)用起到非常積極的作用。

3)進(jìn)行數(shù)值模擬，得到邊坡負(fù)壓排水不同工況下邊坡浸潤線，與簡化的理論浸潤線擬合較好。揭示了邊坡負(fù)壓排水模型的適用性，以及在降深較大的情況下解析解的局限性。

邊坡負(fù)壓排水技術(shù)作為一種新興的高效、經(jīng)濟(jì)的自啟動(dòng)排水技術(shù)，在其內(nèi)在機(jī)理、工程實(shí)踐、優(yōu)化設(shè)計(jì)上還有許多細(xì)節(jié)值得深入研究。