袁 科
(中鐵第四勘察設計院集團有限公司,湖北武漢 430063)
水庫水位變化和降雨會導致坡體穩定性的下降。對于在這些條件下坡體穩定性計算方法,首先要確定浸潤線的位置,然后根據浸潤線來確定滲透壓力,再進行穩定性分析。然而目前浸潤線的確定大多根據經驗,這樣可能造成治理工程的不安全。因此,浸潤線的計算有必要進一步研究。本文對在庫水位等速下降時岸坡地下水浸潤線的公式進行推導計算。
①含水層各向同性、均質,側向無限延伸,具有水平不透水層。
②潛水流為一維流,水位下降前原始潛水面水平。
③水位以v0等速下降。
依據以上假設條件,可得到各向同性土體的一維非穩定滲流運動的基本微分方程[1]
(1)
上式為二階非線性偏微分方程,求其解析解,通常是采用簡化方法,對其進行線性化,然后再求解。簡化方法一般是把含水層厚度H看作為一常量,即用始、末時段潛水流的平均厚度hm來代替,可得到簡化的滲流運動方程為
(2)

(3)
式中k——滲透系數/(m/d);
hm——含水層的平均厚度/m;
μ——給水度或貯水率;
t——水位下降持續時間/d。
浸潤線計算如圖1所示。初始時刻時,即t=0,由假設條件②可知,區內各點初始水位均為h0,0。設距岸坡x處t時刻的地下水位變幅為u(x,t)=h0,0-hx,t,t=0時刻時的水位變幅為:u(x,0)=0。由假設條件③,當水位以v0勻速下降時,在x=0處,t時刻地下水變幅為:u(0,t)=v0t;在x=∞處,可以認為水位無變化,u(∞,t)=0。

圖1 浸潤線計算
根據上述的邊界和初始條件,把上述水位下降的半無限含水層中地下水非穩定滲流歸結為以下模型方程
(4)

(5)
由定積分的分部積分法可知式(5)左邊
(u(x,0)=0)
式(5)右邊
因此
(6)
將方程組(4)中第二個方程兩端同乘e-st并積分得
(7)
因此由以上變換可將模型方程組(4)變換為如下方程組
(8)

(9)
對式(9)進行拉普拉斯(Laplace)逆變換可得
(10)

最終可得水位下降時浸潤線隨時間變化曲線方程
h(x,t)=h0,0-u(x,t)
(11)

從R(λ)的計算式可以看出,R(λ)的計算十分復雜,需要積分才能求得,不利于工程應用。為了得到便于應用的表達式,對圖2所示的曲線進行多項式擬和,得到的擬和公式如下

圖2 λ與R(λ)關系曲線
R(λ)=
(12)
于是可得到庫水位下降時浸潤線計算的簡化公式,其表達式為
hx,t=
(13)
為了驗證所得簡化公式的適用性,運用3D-Modflow對前述公式進行驗證。Visual-Modflow是由加拿大Waterloo水文地質公司在Modflow基礎上,應用現代可視化技術研制而成,是目前世界上最通行且被各國同行一致認可的可視化專業軟件系統。該系統自1994年8月首次在國際上公開發行以來,在工程建設、環境保護、城鄉發展規劃、水資源利用等許多領域的科研和生產中得到了廣泛的應用[3~5],后期經過不斷完善,目前已經成為最為普及的地下水運動數值模擬程序。
為了方便驗證,采用Visual-Modflow分析軟件,對不同滲透系數、水位下降速度及給水度等因素的影響效應進行對比驗證分析。
岸坡巖土體為均質直立岸坡,各向滲透系數相同,其值為0.25 m/d,庫水位下降高度為30 m,下降速度為1 m/d,給水度為0.03,孔隙率為0.25,含水層厚度取升降前后的平均值為40 m,模型的范圍為400 m×60 m×60 m,長度400 m。
利用程序與簡化公式計算對模型進行地下水變動預測,預測時間為30 d,滲透系數按0.05 m/d,0.1 m/d,0.5 m/d,1.0 m/d,2.0 m/d計算,水位變化計算結果對比曲線見圖3。

圖3 V-M程序與簡化公式計算結果
經對比分析可知,對于庫水位等速下降時,滲透系數越小,二者預測浸潤線相差越小;相反,滲透系數越大,此時在坡體前部相差不大,在坡體后部二者有較明顯差異,簡化公式計算結果偏高,但總體而言二者擬合度較好。
參數設置同前,庫水位下降速度按0.25 m/d,0.5 m/d,1.0 m/d,2.0 m/d,5.0 m/d計算,計算結果見圖4。

圖4 V-M程序與簡化公式計算結果
由計算結果可知,二者擬合度較高,其偏差只是在坡體后部偏差要大些。
參數設置同前,給水度分別按0.005,0.008,0.03,0.08,0.12進行計算,計算結果見圖5。

圖5 V-M程序與簡化公式計算結果
由計算結果可知,二者計算結果基本吻合。隨給水度增大,預測結果與簡化公式計算結果偏差降低。
綜合可知,文中得到的簡化公式可以滿足計算精度的要求,且應用計算方便。
由上述分析,令水位的下降高度Δh=v0t,則時間t=Δh/v0,此時有
(14)
從圖2可以看出,隨λ的增大,R(λ)值越小,R(λ)為減函數。也就是說,λ越小,坡體自由水面下降越快;反之,λ越大,坡體自由水面下降越慢。當λ=0,R(λ)=1,坡體中的自由水面與庫水位同步升降;當λ=∞,R(λ)=0,即坡體中自由水面不受庫水位升降的影響。當λ≥2時,R(λ)已非常接近0,此時庫水位的變化對地下水浸潤線的影響可以忽略。因此可以把該處作為庫水位升降對地下水浸潤線產生影響的水平范圍,即
由前面的計算分析可知,浸潤線的確定與坡體滲透系數、給水度、含水層平均厚度、水位升降速度及水位升降高度有密切關系,因此在計算坡體浸潤線時,這些因素的確定就顯得十分重要。
給水度是一個十分重要常用的水文地質參數,它的大小應通過實際測試加以確定。我們通常把它定義為地下水位下降一個單位深度,從地下水位延伸到地表面的單位水平面積巖土柱體,在重力作用下釋出的水的體積。在當實際工作中當無實驗資料時,往往借助經驗公式。
別申斯基(1960)根據砂礫石土樣(最小粒經d=0.25 cm,最小滲透系數k=1.87×10-2cm/s)的試驗結果分析得出給水度μ的簡單經驗公式為
式中:k為滲透系數/(m/d),此式只適合粗顆粒土,不適用于顆粒小于細砂的土質。
給水度代表巖土體的排泄水量的能力,它的大小除了跟滲透系數密切相關外,還與土的密實度有關,即與孔隙率密切相關。文獻[2]根據國內外砂礫石和黏性土的試驗資料,分析求得給水度的經驗公式為
μ=1.137n(0.000 117 5)0.607(6+lgk)
式中n——孔隙率;
k——滲透系數/(cm/s)。
在無試驗資料時,可用此式來確定給水度的大小。
由前面分析可知,為了求解非穩定滲流的基本微分方程,通常采用始、末段潛水流厚度的平均值hm代替式中H,這樣的假定對于升降高度比較小的情況適用,但對于升降高度比較大時誤差就比較大,為此可采用下面的方法來確定含水層平均厚度。
(1)不透水層為水平面
根據滲流分析研究成果,坡體中的浸潤線可簡化為一條拋物線(拋物線方程y2=ax+b),這個假定在堤壩的分析研究中有應用。為了得到浸潤線的平均厚度,可用拋物線模型來確定潛水流的含水層平均厚度,如圖6為庫水位下降和上升兩種情況下含水層厚度計算簡圖。以水位下降為例,假設初始含水層厚h0,水位下降高度為h,影響范圍為s,把點0(0,0),P(s,h)代入拋物線方程,則得到下降后的浸潤線方程為
(15)
圖中浸潤線osp的面積為
將含水層轉化為寬度為s的矩形,可得到庫水位下降時含水層平均厚度hm,即
(16)

圖6 含水層厚度計算
同理可以推導出庫水位上升時含水層的平均厚度為
(17)
式中:h為水位升降高度。
從得到的公式可以看出,當水位升降時含水層的平均厚度僅與水位升降高度和初始含水層厚度有關。在這種情況下很容易確定含水層的厚度,在實際工程應用中十分方便。
(2)不透水層為不規則面的情況
在實際工程中經常遇到不透水層形狀不規則(圖7),此時含水層平均厚度一般可采用如下公式來進行近似計算
(18)
由于一般地質剖面都是用AutoCAD繪制,在AutoCAD里確定OABC區的面積是非常容易的,因此計算平均含水層厚度也就變的十分簡單。

圖7 含水層厚度計算
(1)根據非穩定滲流運動基本微分方程和邊界條件,推導庫水位等速下降條件下坡體內浸潤線的計算公式,并對其進行多項式擬合得到簡化公式,便于實際工程應用。
(2)通過Visual-Modflow預測結果與簡化公式的預測結果對比分析可知,雖然簡化公式在考慮諸多因素情況下,與有限元模擬結果有一定偏差,但多數情況下偏差很小,說明該簡化公式的應用是可行的。但考慮到各因素的綜合影響效應,在應用中應需要對其進行適當的修正。
[1]顧慰慈.滲流計算原理及其應用[M].北京:中國建材工業出版社,2000:12-22
[2]毛昶熙.滲流計算分析與控制(第二版)[M].北京:中國水利水電出版社,2003:352-354
[3]魏玉虎,許 模,等.長江三峽工程奉節庫岸斜坡地下水滲流場模擬分析[J].四川地質學報,2003,23(1):18-22
[4]李 平,盧文喜,等.Visual MODFLOW在地下水數值模擬中的應用[J].工程勘察,2006(3):24-27
[5]周念清,等.MODFLOW在宿遷市地下水資源評價中的應用[J].水文地質與工程地質,2000(6):9-13
[6]王立久.尾礦壩滲流浸潤線計算的實用方法[J].金屬礦山,1987(5)
[7]吳漢杰,林建南,侯曙光.降雨入滲條件下公路邊坡滲流數值分析[J].勘察科學技術,2005(3)
[8]鄭穎人,時衛民,孔位學.庫水位下降時滲透力及地下水浸潤線的計算[J].巖石力學與工程學報,2004(18)
[9]張友誼,胡卸文.庫水位等速上升作用下岸坡地下水浸潤線的計算[J].水文地質工程地質,2007(5)