盧文，李博勇，李 星，于 浩，曾 俊

(1. 武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072；2. 武漢大學 水工巖石力學教育部重點實驗室，武漢 430072)

0 引 言

考慮到復(fù)雜的地質(zhì)歷史(例如地層出現(xiàn)的先后順序)、組成成分和構(gòu)造作用的影響，實際的巖土工程往往面臨具有可變水力特性的非均質(zhì)地層[1,2]。以我國滇西南地區(qū)為例，該地區(qū)高邊坡中常分布透水性差且遇水易軟化的軟弱夾層[3],類似含夾層的邊坡在白鶴灘、溪洛渡等樞紐工程的壩址區(qū)均有大量分布。飽和/非飽和滲流分析方法是查明巖土體工程地質(zhì)、水文地質(zhì)等問題的基礎(chǔ)，在巖土體穩(wěn)定分析、開挖涌水預(yù)測、污染物遷移等領(lǐng)域有較廣泛的應(yīng)用[4, 5]。特別當區(qū)域邊界條件不明確、降雨/蒸發(fā)對地下水流系統(tǒng)影響顯著時，采用非飽和滲流分析方法才能反映出更真實的水文地質(zhì)條件。

隨著計算機技術(shù)的發(fā)展，利用數(shù)值方法求解飽和/非飽和滲流的研究越來越豐富。Neuman[6]最早將Galerkin型有限單元法應(yīng)用于求解二維飽和-非飽和滲流的問題，他提出了用零壓面替代活動的自由面的方法，從而實現(xiàn)飽和區(qū)與非飽和區(qū)的統(tǒng)一計算。赤井浩一[7]通過砂槽模型試驗?zāi)M吸濕與脫濕情形下的非飽和滲流，驗證了Newman提出的數(shù)值解法。吳良驥等[8]基于有限差分法提出了飽和-非飽和滲流數(shù)學模型并編寫了相應(yīng)的計算程序，提升了質(zhì)量平衡精度。高驥等[9]建立飽和區(qū)與非飽和區(qū)的整體數(shù)學模型，應(yīng)用到堤壩滲流分析中。張家發(fā)[10]建立了考慮降雨入滲邊界條件的三維飽和/非飽和滲流數(shù)學模型，并用赤井浩一的砂槽試驗進行驗證，編寫的有限元計算程序應(yīng)用到三峽船閘高邊坡滲流分析中。孫役[11]建立了有暴雨入滲下三維裂隙網(wǎng)絡(luò)飽和-非飽和滲流數(shù)學模型，并將其用于龍灘水電站左岸高邊坡的滲流分析中。胡云進[12]將編制的三維飽和/非飽和滲流有限元計算程序用于模擬溪洛渡電站水墊塘區(qū)高邊坡的滲流。然而，飽和/非飽和滲流數(shù)值求解中的一大困難在于邊界具有高度的非線性，許多學者對此展開了大量的探索工作。朱岳明等[13, 14]提出了包含降雨與蒸發(fā)的非飽和滲流逸出邊界的數(shù)值處理方法，并將其用于三峽2號永久船閘高邊坡的滲流分析中。Borsi等[15]通過數(shù)學方法證明了互補邊界的統(tǒng)一格式。胡冉等[16]建立了拋物線型變分不等式(PVI)提法，有效解決了邊界高度非線性問題，可實現(xiàn)對飽和-非飽和滲流邊界條件的精細化求解，擴大了飽和-非飽和滲流問題數(shù)值模擬的適用范圍。

土水特征曲線是描述非飽和區(qū)飽和度與基質(zhì)吸力之間本構(gòu)關(guān)系的曲線。戚國慶等[17]對已有的土水特征曲線數(shù)學模型進行研究，將其劃分為4種類型，提出冪函數(shù)多項式形式的通用模型，形式簡單且擬合精度高。徐紹輝[18]等對常用的幾個經(jīng)驗?zāi)Ｐ瓦M行數(shù)值求解，證實在使用V-G模型[19]擬合粗粒土及基質(zhì)吸力低的細粒土的土水特征曲線時，得到的結(jié)果精度較高。

本文基于水體質(zhì)量守恒方程及V-G模型，采用拋物線型變分不等式(PVI)方法對含夾層邊坡進行飽和-非飽和滲流模擬。通過與Rulon JJ等[4]在試驗中觀測到的水頭及流量結(jié)果對比，驗證了算法的可靠性，并探究了現(xiàn)場尺度下夾層滲流特性、發(fā)育程度及位置對邊坡滲流場的影響。

1 飽和/非飽和滲流理論

1.1 飽和/非飽和滲流數(shù)學模型

由水體質(zhì)量守恒方程(也稱作Richards方程)，可建立區(qū)域Ω內(nèi)巖土體的飽和/非飽和滲流控制方程：

(1)

式中：h=pw/γw為壓力水頭；pw為孔隙水壓力；γw為水的容重；C=?Se/?h為容水度，在非飽和區(qū)可表示為壓力水頭h或基質(zhì)吸力s=pa-pw的函數(shù)，在飽和區(qū)為0；Se為飽和度；pa為孔隙氣壓力；Ss為單位貯存量，對非飽和區(qū)巖土體，其值為零；β在飽和區(qū)等于1，非飽和區(qū)等于0。

巖土體中的飽和/非飽和滲流規(guī)律可用廣義Darcy定律表示為：

v=-kr(h)K▽φ

(2)

式中：φ=h+z為總水頭；z為垂直坐標分量；K為巖土體的飽和滲透張量；kr為相對滲透率，可表示為飽和度Se、基質(zhì)吸力s或壓力水頭h的函數(shù)。

為了確定區(qū)域Ω內(nèi)唯一的滲流場分布，還需要給定初始條件及邊界條件。

式(1)的初始條件如下：

h(x,y,z,t0)=h0(x,y,z) (inΩ)

(3)

式中：h0為初始水頭場。

式(1)還應(yīng)滿足如下邊界條件：

(1) 水頭邊界條件。

(4)

(2) 流量邊界條件。

(5)

(3) 第三類邊界條件。此外，飽和溢出邊界、降雨邊界和蒸發(fā)邊界的形式比較復(fù)雜，但經(jīng)過嚴格的數(shù)學方法可表示為統(tǒng)一的互補形式[15, 16]，將其稱作第三類邊界條件?；パa形式如下：

(6)

(onΓT=Γs∪Γue∪Γui)

1.2 土水特征曲線

質(zhì)量守恒方程及運動方程涉及壓力水頭h、流速v、飽和度Se3個未知量，因此求解方程還需給出巖土體飽和度Se與壓力水頭h(或基質(zhì)吸力s)之間的關(guān)系，即土水特征曲線?；贛ualem理論，van Genuchten[19]給出了一種光滑連續(xù)的土水曲線模型，稱作VG模型：

(7)

式中：α與進氣狀態(tài)相關(guān)，約等于進氣壓力值的倒數(shù)；n與孔徑分布相關(guān)，m與土水特征曲線的整體對稱性相關(guān)，本文取m=1-1/n。相應(yīng)的kr～Se表達式為：

(8)

2 飽和/非飽和滲流有限元模擬方法及實驗驗證

2.1 有限元格式

通過對式(1)進行空間域離散和時間域差分，并選取水頭為基本未知量，可建立飽和/非飽和滲流問題的有限元數(shù)值計算格式如下：

(C+ΔtθK)φi+1=[C-Δt(1-θ)K]φi+Δtq

(9)

式中：

(10)

(11)

(12)

式中：N為有限單元的形函數(shù)矩陣；B為有限元模型的幾何矩陣；i為時間步；Δt為i時步到i+1時步的時間步長；θ為積分常數(shù)(θ≥0.5)。

由于互補邊界具有高度的非線性特征，在求解時具有較大的難度。采用胡冉等[16]建立的拋物線型變分不等式(PVI)提法，可以將互補邊界條件轉(zhuǎn)化成自然邊界條件，大大減小數(shù)值求解的難度，改善數(shù)值計算的收斂性和穩(wěn)定性。

2.2 驗 證

為了對上述有限元方法的正確性進行驗證，本文模擬了Rulon JJ等[4]設(shè)計的含夾層邊坡非飽和滲流實驗，有限元模型及輸入?yún)?shù)依據(jù)該物理實驗確定。如圖1(a)所示，該物理模型長2.44 m，高1 m，厚0.1 m，由中砂和細砂兩種材料填充，其中細砂夾層的厚度為0.1 m。ED邊界接受降雨，降雨量為2.1×10-4m/s(1.26 cm/min)；AF邊界為定水頭邊界，水位高0.3 m；EF邊界為溢出邊界，其余邊界為隔水邊界。本文建立同等尺寸的有限元模型，劃分出8 000 個單元，16 362 個節(jié)點，見圖1(b)。

圖1 實驗?zāi)Ｐ图坝邢拊嬎隳Ｐ?單位：m)Fig.1 Schematic diagram of the model and finite element model

依據(jù)Rulon JJ等[4]測量的參數(shù)值，并考慮到裝砂方式及測量儀器安裝將導(dǎo)致試樣參數(shù)發(fā)生一定程度的改變，經(jīng)過一系列模擬確定中砂和細砂的參數(shù)取值如表1所示。

表1 材料參數(shù)表Tab.1 Parameters of medium sand and fine sand

數(shù)值模擬得到的地下水位線及實驗觀測到的水位線繪制在圖2中，結(jié)果顯示在細砂層上下自由面發(fā)生分層現(xiàn)象，且水位線計算值與測量值吻合較好。為了量化對比結(jié)果，圖3將Rulon JJ等[4]觀測到的55個測點的壓力水頭值與模擬的壓力水頭值進行對比，數(shù)據(jù)點均分布在直線ψ觀測值=ψ計算值附近，這表明數(shù)值模擬的結(jié)果精確度很高。

圖2 地下水位線對比圖Fig.2 Comparison of simulated and observed water table

圖3 測點壓力水頭對比圖Fig.3 Comparison of simulated and observed pressure head

此外，實驗結(jié)果顯示出明顯的分層地下水流現(xiàn)象，因而對上層及斜坡底部的流量分別進行統(tǒng)計。且考慮到實驗中流量收集器可能存在一定誤差，以流量所占比例代替流量值進行對比，如圖4所示。其中，QU代表上層溢出量所占比例，QL代表斜坡基部滲流流量所占比例。結(jié)果顯示，流量占比模擬值與監(jiān)測值較為接近，再次驗證了本文計算方法的正確性。

圖4 流量對比圖Fig.4 Comparison of simulated and observed outflow rates

3 現(xiàn)場尺度邊坡參數(shù)敏感性分析

水平或緩傾角的泥質(zhì)夾層在我國西南地區(qū)發(fā)育廣泛。見圖5(a)所示，西南某水電站左岸壩肩邊坡發(fā)育5條傾角約為9°～30°的軟弱夾層，夾層主要含凝灰?guī)r泥質(zhì)物[20]。以白鶴灘樞紐區(qū)某勘探線[如圖5(b)所示]為例，該剖面發(fā)育的水平結(jié)構(gòu)面C2、C3-1泥化、軟化現(xiàn)象明顯，其間泥化帶在垂直層面方向有不同程度的阻水性。

然而，由于巖土體的滲流特性具有空間變異性，能代表現(xiàn)場尺度的滲流參數(shù)很難確定。因此，本節(jié)建立了一個含相對弱透水夾層的理想邊坡模型，模擬邊坡在歷經(jīng)多年降雨循環(huán)之后的天然滲流場分布情況，探究夾層的非飽和水力參數(shù)(α和n)、與相鄰巖土體的滲透系數(shù)比及其厚度與位置對邊坡滲流場的影響程度，對于分析含類似結(jié)構(gòu)巖土體的區(qū)域水文地質(zhì)條件及確定其滲流參數(shù)有重要的指導(dǎo)意義。

圖5 西南地區(qū)某邊坡剖面[20]Fig.5 Some slope profiles in southwest China[20]

如圖6所示，本文對上述坡體結(jié)構(gòu)進行概化，建立了一個現(xiàn)場尺度的考慮單一水平夾層的邊坡模型，該理想邊坡兩側(cè)分別以河床中心線和分水嶺為界，故AB及CD邊界取為隔水邊界；FGA邊界為定水頭邊界，水頭值60 m；DEF邊界為互補邊界，降雨量參考Rulon JJ等[4]定為2.54×10-8m/s(800 mm/y)；底部為隔水邊界。共劃分20 120 個單元，40 854 個節(jié)點。材料參數(shù)設(shè)計如表2所示，材料1的滲透系數(shù)參照Rulon JJ等[4]的取值，對夾層和材料1的滲透系數(shù)比進行敏感性分析，K2/K1分別取1/20，1/25和1/30。同時對于夾層材料土水特征曲線的擬合參數(shù)α和n分別進行敏感性分析，其取值范圍參照van Genuchten[19]獲取的粉質(zhì)壤土及冰磧土的土水特征曲線擬合參數(shù)。以夾層α=0.06，n=3.2作為正常工況。初始步長取0.07 d，最大步長10 d，計算時長為30 a，初始地下水位線與坡頂重合。

圖6 含夾層的邊坡模型(單位：m)Fig.6 Model of layered slope

材料ks/(m·s-1)θsθrα/m-1nSs/m-1材料11×10-70.460.180.063.21.54×10-5材料23.3×10-9,4×10-9,5×10-90.400.050.02,0.04,0.06,0.08,1.02.0,2.6,3.2,3.8,4.42.38×10-4

由圖7～圖9可見，當邊坡內(nèi)發(fā)育一低滲透性的夾層時，地下水在夾層處形成分層水流。圖7顯示減小夾層與鄰近巖土體的滲透系數(shù)比，夾層的相對弱滲透性越明顯，對地下水下滲的阻隔作用越顯著，因此上層地下水位升高，而下層地下水位降低。由圖8可知，增大夾層的非飽和參數(shù)α，上層地下水位呈現(xiàn)升高的趨勢，相反下層地下水位降低；夾層內(nèi)非飽和區(qū)在靠近坡表處范圍減小，但向山體內(nèi)側(cè)延伸至更深處。這是由于隨著夾層α值增大，進氣值減小，吸濕與脫濕過程都更容易發(fā)生，臨近坡表處雨水更容易入滲，因此該處夾層的非飽和范圍減小。但在相同的負壓條件下，α越大飽和度越低，相對滲透系數(shù)也就越小，故上層水位升高，而下層水位降低，下層的非飽和區(qū)也向山體深處擴展。然而，圖9的結(jié)果表明增大夾層的非飽和參數(shù)n，該邊坡內(nèi)滲流場表現(xiàn)出與增大α相反的變化規(guī)律。這是由于在低負壓區(qū)，n值越大，飽和度隨基質(zhì)負壓的變化反而更慢，即吸濕與脫濕過程更難發(fā)生，因此呈現(xiàn)出與增大α相反的水位分布規(guī)律。

此前的一些研究表明非飽和土的抗剪強度隨著基質(zhì)吸力的增加而增加[21]。由此可見當?shù)叵滤謱蝇F(xiàn)象越明顯，邊坡上層的非飽和區(qū)越小，相同位置的基質(zhì)吸力越小，導(dǎo)致其穩(wěn)定性越差；而下層非飽和區(qū)范圍越大，故深部巖土體的抗剪強度越高，穩(wěn)定性越好。

圖7 含夾層的邊坡模型(單位：m)Fig.7 Model of layered slope

此外，圖8顯示了分別改變夾層α(或n)取值時，上層溢出流量(除去坡面降雨補給量)的變化規(guī)律，當夾層的非飽和參數(shù)α增大(或n減小)時，上層溢出流量增大，與地下水位線變化規(guī)律一致。

圖9顯示了當夾層高度變化時，該邊坡壓力水頭的分布情況，圖9中h表示夾層中心線所在高程。結(jié)果顯示，隨著夾層高度的增加，上層水位升高，而下層的非飽和區(qū)范圍更大，地下水分層現(xiàn)象更為顯著。圖10為夾層高度不變，不同厚度對應(yīng)的壓力水頭分布圖。結(jié)果表明夾層的厚度越厚，其對邊坡滲流場的影響更顯著：隨著夾層厚度的增加，夾層上的飽和范圍增大，邊坡最大滲透壓力減小，地下水的分層現(xiàn)象也更明顯。

圖8 滲流場分布對比圖(改變夾層α值)(單位：m)Fig.8 Effect of the interlayer’s α on the seepage field

圖9 滲流場分布對比圖(改變夾層n值)(單位：m)Fig.9 Effect of the interlayer’s n on the seepage field

圖10 上層溢出流量對比圖(α或n變化)Fig.10 Comparison of the outflow rates across the upper seepage faces when α or n changes

圖11 壓力水頭分布對比圖(改變夾層高度)(單位：m)Fig.11 Effect of the elevation of interlayer on the pressure head distribution

圖12 壓力水頭分布對比圖(改變夾層厚度)(單位：m)Fig.12 Effect of the thickness of the interlayer on the pressure head distribution

4 結(jié) 語

(1)本文采用有限單元法模擬含夾層邊坡飽和-非飽和滲流問題，能精確且高效地模擬飽和/非飽和溢出邊界，模擬結(jié)果與實驗數(shù)據(jù)吻合良好，證實了該算法的可靠性。

(2)歷經(jīng)多年降雨后，含低滲透性夾層的邊坡易形成多層地下水流的現(xiàn)象。夾層的滲流特性與發(fā)育程度直接影響該區(qū)域地下水的分布情況：夾層的非飽和參數(shù)α增大(或n減小)，夾層與鄰近巖土體滲透系數(shù)比減小，高度(或厚度)增大，將導(dǎo)致地下水分層現(xiàn)象更加顯著。

(3)采用飽和/非飽和滲流分析方法模擬出更真實可靠的地下水分布情況也給邊坡穩(wěn)定分析提供了參考依據(jù)。具體表現(xiàn)為地下水分層現(xiàn)象越明顯，越不利于邊坡上層的穩(wěn)定，但有利于深層巖土體的穩(wěn)定。