邊坡負壓排水非穩定流研究

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

浙江大學海洋學院，浙江 舟山 316021

0 引言

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

A.鉆孔的孔口；B.排水管進水口。

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

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

1 數學模型的建立

1.1 問題描述與基本假設

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

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

1.2 控制方程

式(1)為二維潛水運動的控制方程[6]：

(1)

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

(2)

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

上文提到，在負壓排水模型中，集水腔體相對于整個坡體尺度很小，故簡化為一個點，此處用狄拉克函數表征點匯，即上式中源匯項函數W定義為

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

(3)

式中：Q為源匯量,m3/s；δ為狄拉克函數符號；(x0,z0)為匯點坐標。

1.3 方程的定解條件

1.3.1 初始條件

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

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

(4)

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

1.3.2 邊界條件

邊坡地下水的流動問題主要有3類邊界條件。

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

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

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

(5)

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

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

(6)

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

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

(7)

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

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

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

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

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

(8)

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

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

zwt=zwt(x,y,t)。

則邊界條件可以描述為，當點處于潛水面上，即點的高程z等于zwt時，水頭函數值h也等于zwt(當基準取在隔水地板)，即

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

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

(9)

式中:H為潛水面各點水頭；Bwt為潛水面邊界。

由水均衡原理，控制體各側面流入及流出的結果將導致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)

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

(11)

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

(12)

又有：

(13)

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

(14)

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

(15)

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

(16)

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

(17)

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

2 問題的求解

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

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

(18)

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

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

(19)

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

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

的解，而Hj是方程

的解。令

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

(20)

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

(21)

的解。

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

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

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

子問題一：

(22)

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

子問題二：

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

s(∞,z,t)=0；

(23)

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

對z做傅里葉變換：

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

(24)

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

(25)

接下來對F2作逆變換。觀察函數形式，我們發現對像函數積分使用極坐標系更為方便。則

y=rsinθ；

ξ=ρsinφ。

則

(26)

將原函數也用極坐標表達，即

x=rcosθ，y=rsinθ。

則

(27)

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

(28)

式中，J0(a1)是一個0階貝塞爾函數。則式(28)可轉換為

(29)

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

(30)

其中：

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

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

3 浸潤線的數值模擬分析

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

圖4 邊坡模型及網格劃分示意圖

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

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

當有自啟動負壓排水干擾后，邊坡總水頭等值線圖如圖5所示。

箭頭方向代表地下水流速方向，箭頭長度代表該點流速大小。

設置兩種工況分別計算了自啟動負壓排水方法下數值模擬結果和理論解得到的邊坡浸潤線曲線，并進行對比分析，結果如圖6、7所示。

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

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

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

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

4 結論與建議

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

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

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

邊坡負壓排水技術作為一種新興的高效、經濟的自啟動排水技術，在其內在機理、工程實踐、優化設計上還有許多細節值得深入研究。