袁 科

(中鐵第四勘察設(shè)計院集團(tuán)有限公司,湖北武漢 430063)

水庫水位變化和降雨會導(dǎo)致坡體穩(wěn)定性的下降。對于在這些條件下坡體穩(wěn)定性計算方法,首先要確定浸潤線的位置,然后根據(jù)浸潤線來確定滲透壓力,再進(jìn)行穩(wěn)定性分析。然而目前浸潤線的確定大多根據(jù)經(jīng)驗,這樣可能造成治理工程的不安全。因此,浸潤線的計算有必要進(jìn)一步研究。本文對在庫水位等速下降時岸坡地下水浸潤線的公式進(jìn)行推導(dǎo)計算。

1 基本假設(shè)及計算模型

1.1 基本假設(shè)

①含水層各向同性、均質(zhì),側(cè)向無限延伸,具有水平不透水層。

②潛水流為一維流,水位下降前原始潛水面水平。

③水位以v0等速下降。

依據(jù)以上假設(shè)條件,可得到各向同性土體的一維非穩(wěn)定滲流運(yùn)動的基本微分方程[1]

(1)

上式為二階非線性偏微分方程,求其解析解,通常是采用簡化方法,對其進(jìn)行線性化,然后再求解。簡化方法一般是把含水層厚度H看作為一常量,即用始、末時段潛水流的平均厚度hm來代替,可得到簡化的滲流運(yùn)動方程為

(2)

(3)

式中k——滲透系數(shù)/(m/d)；

hm——含水層的平均厚度/m；

μ——給水度或貯水率；

t——水位下降持續(xù)時間/d。

1.2 計算模型建立

浸潤線計算如圖1所示。初始時刻時,即t=0,由假設(shè)條件②可知,區(qū)內(nèi)各點初始水位均為h0,0。設(shè)距岸坡x處t時刻的地下水位變幅為u(x,t)=h0,0-hx,t,t=0時刻時的水位變幅為：u(x,0)=0。由假設(shè)條件③,當(dāng)水位以v0勻速下降時,在x=0處,t時刻地下水變幅為：u(0,t)=v0t；在x=∞處,可以認(rèn)為水位無變化,u(∞,t)=0。

圖1 浸潤線計算

根據(jù)上述的邊界和初始條件,把上述水位下降的半無限含水層中地下水非穩(wěn)定滲流歸結(jié)為以下模型方程

(4)

1.3 微分方程求解

(5)

由定積分的分部積分法可知式(5)左邊

(u(x,0)=0)

式(5)右邊

因此

(6)

將方程組(4)中第二個方程兩端同乘e-st并積分得

(7)

因此由以上變換可將模型方程組(4)變換為如下方程組

(8)

(9)

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

(10)

最終可得水位下降時浸潤線隨時間變化曲線方程

h(x,t)=h0,0-u(x,t)

(11)

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

圖2 λ與R(λ)關(guān)系曲線

R(λ)=

(12)

于是可得到庫水位下降時浸潤線計算的簡化公式,其表達(dá)式為

hx,t=

(13)

2 簡化公式驗證分析

為了驗證所得簡化公式的適用性,運(yùn)用3D-Modflow對前述公式進(jìn)行驗證。Visual-Modflow是由加拿大Waterloo水文地質(zhì)公司在Modflow基礎(chǔ)上,應(yīng)用現(xiàn)代可視化技術(shù)研制而成,是目前世界上最通行且被各國同行一致認(rèn)可的可視化專業(yè)軟件系統(tǒng)。該系統(tǒng)自1994年8月首次在國際上公開發(fā)行以來,在工程建設(shè)、環(huán)境保護(hù)、城鄉(xiāng)發(fā)展規(guī)劃、水資源利用等許多領(lǐng)域的科研和生產(chǎn)中得到了廣泛的應(yīng)用[3～5],后期經(jīng)過不斷完善,目前已經(jīng)成為最為普及的地下水運(yùn)動數(shù)值模擬程序。

為了方便驗證,采用Visual-Modflow分析軟件,對不同滲透系數(shù)、水位下降速度及給水度等因素的影響效應(yīng)進(jìn)行對比驗證分析。

2.1 滲透系數(shù)的影響效應(yīng)

岸坡巖土體為均質(zhì)直立岸坡,各向滲透系數(shù)相同,其值為0.25 m/d,庫水位下降高度為30 m,下降速度為1 m/d,給水度為0.03,孔隙率為0.25,含水層厚度取升降前后的平均值為40 m,模型的范圍為400 m×60 m×60 m,長度400 m。

利用程序與簡化公式計算對模型進(jìn)行地下水變動預(yù)測,預(yù)測時間為30 d,滲透系數(shù)按0.05 m/d,0.1 m/d,0.5 m/d,1.0 m/d,2.0 m/d計算,水位變化計算結(jié)果對比曲線見圖3。

圖3 V-M程序與簡化公式計算結(jié)果

經(jīng)對比分析可知,對于庫水位等速下降時,滲透系數(shù)越小,二者預(yù)測浸潤線相差越小；相反,滲透系數(shù)越大,此時在坡體前部相差不大,在坡體后部二者有較明顯差異,簡化公式計算結(jié)果偏高,但總體而言二者擬合度較好。

2.2 水位下降速度的影響效應(yīng)

參數(shù)設(shè)置同前,庫水位下降速度按0.25 m/d,0.5 m/d,1.0 m/d,2.0 m/d,5.0 m/d計算,計算結(jié)果見圖4。

圖4 V-M程序與簡化公式計算結(jié)果

由計算結(jié)果可知,二者擬合度較高,其偏差只是在坡體后部偏差要大些。

2.3 給水率的影響效應(yīng)

參數(shù)設(shè)置同前,給水度分別按0.005,0.008,0.03,0.08,0.12進(jìn)行計算,計算結(jié)果見圖5。

圖5 V-M程序與簡化公式計算結(jié)果

由計算結(jié)果可知,二者計算結(jié)果基本吻合。隨給水度增大，預(yù)測結(jié)果與簡化公式計算結(jié)果偏差降低。

綜合可知,文中得到的簡化公式可以滿足計算精度的要求,且應(yīng)用計算方便。

3 浸潤線影響因素

由上述分析,令水位的下降高度Δh=v0t,則時間t=Δh/v0,此時有

(14)

從圖2可以看出,隨λ的增大,R(λ)值越小,R(λ)為減函數(shù)。也就是說,λ越小,坡體自由水面下降越快；反之,λ越大,坡體自由水面下降越慢。當(dāng)λ=0,R(λ)=1,坡體中的自由水面與庫水位同步升降；當(dāng)λ=∞,R(λ)=0,即坡體中自由水面不受庫水位升降的影響。當(dāng)λ≥2時,R(λ)已非常接近0,此時庫水位的變化對地下水浸潤線的影響可以忽略。因此可以把該處作為庫水位升降對地下水浸潤線產(chǎn)生影響的水平范圍,即

4 給水度與含水層厚度確定

由前面的計算分析可知,浸潤線的確定與坡體滲透系數(shù)、給水度、含水層平均厚度、水位升降速度及水位升降高度有密切關(guān)系,因此在計算坡體浸潤線時,這些因素的確定就顯得十分重要。

4.1 給水度的確定

給水度是一個十分重要常用的水文地質(zhì)參數(shù),它的大小應(yīng)通過實際測試加以確定。我們通常把它定義為地下水位下降一個單位深度,從地下水位延伸到地表面的單位水平面積巖土柱體,在重力作用下釋出的水的體積。在當(dāng)實際工作中當(dāng)無實驗資料時,往往借助經(jīng)驗公式。

別申斯基(1960)根據(jù)砂礫石土樣(最小粒經(jīng)d=0.25 cm,最小滲透系數(shù)k=1.87×10-2cm/s)的試驗結(jié)果分析得出給水度μ的簡單經(jīng)驗公式為

式中：k為滲透系數(shù)/(m/d),此式只適合粗顆粒土,不適用于顆粒小于細(xì)砂的土質(zhì)。

給水度代表巖土體的排泄水量的能力,它的大小除了跟滲透系數(shù)密切相關(guān)外,還與土的密實度有關(guān),即與孔隙率密切相關(guān)。文獻(xiàn)[2]根據(jù)國內(nèi)外砂礫石和黏性土的試驗資料,分析求得給水度的經(jīng)驗公式為

μ=1.137n(0.000 117 5)0.607(6+lgk)

式中n——孔隙率；

k——滲透系數(shù)/(cm/s)。

在無試驗資料時,可用此式來確定給水度的大小。

4.2 含水層厚度的確定

由前面分析可知,為了求解非穩(wěn)定滲流的基本微分方程,通常采用始、末段潛水流厚度的平均值hm代替式中H,這樣的假定對于升降高度比較小的情況適用,但對于升降高度比較大時誤差就比較大,為此可采用下面的方法來確定含水層平均厚度。

(1)不透水層為水平面

根據(jù)滲流分析研究成果,坡體中的浸潤線可簡化為一條拋物線(拋物線方程y2=ax+b),這個假定在堤壩的分析研究中有應(yīng)用。為了得到浸潤線的平均厚度,可用拋物線模型來確定潛水流的含水層平均厚度,如圖6為庫水位下降和上升兩種情況下含水層厚度計算簡圖。以水位下降為例,假設(shè)初始含水層厚h0,水位下降高度為h,影響范圍為s,把點0(0,0),P(s,h)代入拋物線方程,則得到下降后的浸潤線方程為

(15)

圖中浸潤線osp的面積為

將含水層轉(zhuǎn)化為寬度為s的矩形,可得到庫水位下降時含水層平均厚度hm,即

(16)

圖6 含水層厚度計算

同理可以推導(dǎo)出庫水位上升時含水層的平均厚度為

(17)

式中：h為水位升降高度。

從得到的公式可以看出,當(dāng)水位升降時含水層的平均厚度僅與水位升降高度和初始含水層厚度有關(guān)。在這種情況下很容易確定含水層的厚度,在實際工程應(yīng)用中十分方便。

(2)不透水層為不規(guī)則面的情況

在實際工程中經(jīng)常遇到不透水層形狀不規(guī)則(圖7),此時含水層平均厚度一般可采用如下公式來進(jìn)行近似計算

(18)

由于一般地質(zhì)剖面都是用AutoCAD繪制,在AutoCAD里確定OABC區(qū)的面積是非常容易的,因此計算平均含水層厚度也就變的十分簡單。

圖7 含水層厚度計算

5 結(jié)論

(1)根據(jù)非穩(wěn)定滲流運(yùn)動基本微分方程和邊界條件,推導(dǎo)庫水位等速下降條件下坡體內(nèi)浸潤線的計算公式,并對其進(jìn)行多項式擬合得到簡化公式,便于實際工程應(yīng)用。

(2)通過Visual-Modflow預(yù)測結(jié)果與簡化公式的預(yù)測結(jié)果對比分析可知,雖然簡化公式在考慮諸多因素情況下,與有限元模擬結(jié)果有一定偏差,但多數(shù)情況下偏差很小,說明該簡化公式的應(yīng)用是可行的。但考慮到各因素的綜合影響效應(yīng),在應(yīng)用中應(yīng)需要對其進(jìn)行適當(dāng)?shù)男拚?/p>

[1]顧慰慈.滲流計算原理及其應(yīng)用[M].北京：中國建材工業(yè)出版社,2000：12-22

[2]毛昶熙.滲流計算分析與控制(第二版)[M].北京：中國水利水電出版社，2003：352-354

[3]魏玉虎,許 模，等.長江三峽工程奉節(jié)庫岸斜坡地下水滲流場模擬分析[J].四川地質(zhì)學(xué)報,2003,23(1)：18-22

[4]李 平,盧文喜，等.Visual MODFLOW在地下水?dāng)?shù)值模擬中的應(yīng)用[J].工程勘察,2006(3)：24-27

[5]周念清，等.MODFLOW在宿遷市地下水資源評價中的應(yīng)用[J].水文地質(zhì)與工程地質(zhì),2000(6)：9-13

[6]王立久.尾礦壩滲流浸潤線計算的實用方法[J].金屬礦山,1987(5)

[7]吳漢杰,林建南,侯曙光.降雨入滲條件下公路邊坡滲流數(shù)值分析[J].勘察科學(xué)技術(shù),2005(3)

[8]鄭穎人,時衛(wèi)民,孔位學(xué).庫水位下降時滲透力及地下水浸潤線的計算[J].巖石力學(xué)與工程學(xué)報,2004(18)

[9]張友誼,胡卸文.庫水位等速上升作用下岸坡地下水浸潤線的計算[J].水文地質(zhì)工程地質(zhì),2007(5)