王平凡， 畢金鋒， 羅先啟

(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240)

0 引 言

堤壩結(jié)構(gòu)或邊坡工程長期受到庫水位或降雨的影響，土體結(jié)構(gòu)由于滲流影響會發(fā)生內(nèi)部失穩(wěn)，內(nèi)侵蝕作用已經(jīng)成為導(dǎo)致潰壩及滑坡的主要因素之一。如1998年，長江九江段4號閘與5號閘之間決堤，潰口發(fā)生前堤防下游坡后有較多滲透變形跡象，隨后發(fā)展為堤防潰口。

目前，在諸多研究土體內(nèi)侵蝕的試驗(yàn)中，主要研究的是影響內(nèi)侵蝕的因素和細(xì)顆粒流失規(guī)律。Benamar等[1]在不同的流動條件和土體力學(xué)參數(shù)下進(jìn)行滲流試驗(yàn)，研究內(nèi)侵蝕的發(fā)生機(jī)理和動力學(xué)特性；李喜安等[2]利用解析方法推導(dǎo)內(nèi)侵蝕發(fā)生的臨界條件，指出內(nèi)侵蝕發(fā)生的臨界條件應(yīng)由具體受力情況來分析；Marot等[3]研究內(nèi)侵蝕過程中土體密度、孔隙壓力和被侵蝕顆粒的質(zhì)量變化，表明細(xì)顆粒流失導(dǎo)致孔隙壓力增加，進(jìn)而出現(xiàn)局部失穩(wěn)；Chang等[4]對內(nèi)侵蝕試樣進(jìn)行三軸壓縮試驗(yàn)，測試土體力學(xué)特性受內(nèi)侵蝕作用的影響程度，當(dāng)細(xì)顆粒流失量達(dá)到一定量級后，土體由剪脹變?yōu)榧艨s。Bonelli等[5]對管涌中細(xì)顆粒流失規(guī)律進(jìn)行研究，建立流體、被侵蝕顆粒和土體整體的質(zhì)量守恒方程；Sterpi[6]提出一種計(jì)算細(xì)顆粒流失速率的方法，通過可被侵蝕細(xì)顆粒的質(zhì)量計(jì)算在任意時(shí)刻下被侵蝕的細(xì)顆粒的質(zhì)量；吳夢喜等[7]對細(xì)顆粒流失量不同的土體進(jìn)行三軸和側(cè)限壓縮試驗(yàn)，將土體的應(yīng)變增量分成由應(yīng)力和模量衰化引起的應(yīng)變增量。

另外，還有學(xué)者借助有限元模擬內(nèi)侵蝕過程來分析內(nèi)侵蝕特性。Sibille等[8]利用離散元和格子波爾茲曼法來模擬內(nèi)侵蝕過程中細(xì)顆粒從骨架顆粒脫離并隨水流運(yùn)移，通過流體作用在顆粒上的功來耦合兩者之間的相互作用；Zhang等[9]根據(jù)熱力學(xué)原理和多孔介質(zhì)理論建立內(nèi)侵蝕結(jié)構(gòu)模型，解釋管涌現(xiàn)象在臨界水力梯度的區(qū)域加劇原因；Sterpi[6]結(jié)合被侵蝕顆粒的質(zhì)量守恒方程和基于滲流實(shí)驗(yàn)得到的侵蝕規(guī)律，建立起細(xì)顆粒侵蝕與運(yùn)移的有限元模型；劉忠玉等[10]建立臨界水頭梯度與細(xì)顆粒流失的關(guān)系，采用有限元方法分析細(xì)顆粒流失對工程穩(wěn)定性的影響。羅玉龍等[11]結(jié)合潛蝕本構(gòu)方程和土體三相質(zhì)量守恒方程，建立有限元數(shù)值模型研究管涌過程中的內(nèi)侵蝕作用。

在上述文獻(xiàn)中，試驗(yàn)所采用的樣本是重塑土，這與天然土存在較大的差異，試驗(yàn)得到的內(nèi)侵蝕開始和結(jié)束條件與真實(shí)情況有些許差別。而且有限元模擬采用的控制方程沒有將質(zhì)量守恒方程、連續(xù)性方程、滲流控制方程耦合起來，只是單獨(dú)考慮一個(gè)方面或只是耦合了其中兩個(gè)控制方程。本文參考其他學(xué)者所做的細(xì)顆粒流失試驗(yàn)和內(nèi)侵蝕控制方程，以錦屏呷爬滑坡為例，采用有限元的方法，耦合質(zhì)量守恒方程、連續(xù)性方程和滲流方程建立改進(jìn)后的內(nèi)侵蝕數(shù)值模型，研究內(nèi)侵蝕發(fā)生時(shí)不同的水頭對土體的孔隙率、被侵蝕顆粒的質(zhì)量、細(xì)顆粒流失速率的影響。

1 內(nèi)侵蝕形成條件

水流流過土體內(nèi)部而將細(xì)顆粒帶走的現(xiàn)象稱為內(nèi)侵蝕。內(nèi)侵蝕過程中，土體內(nèi)部細(xì)顆粒在滲透力作用下發(fā)生松動并脫離主骨架，從而造成土體的局部孔隙率增大，隨著水對土體的內(nèi)侵蝕作用的加劇，會在局部形成優(yōu)勢路徑或孔洞，土體骨架遭到破壞，局部強(qiáng)度削弱，容易造成土體整體在優(yōu)勢路徑所在平面上發(fā)生剪切破壞。

沙金煊[12]通過對作用在單個(gè)土顆粒上的滲透力的計(jì)算，提出滲透侵蝕臨界水力坡降的近似公式

(1)

式中：r為可被侵蝕顆粒粒徑;k為滲透系數(shù);φ為孔隙率。

《水利水電工程地質(zhì)勘察規(guī)范》中對于間斷級配的土體，以間斷區(qū)間的中間值作為可被侵蝕的細(xì)顆粒的粒徑。根據(jù)錦屏一級水電站壩址相關(guān)勘察報(bào)告及工程經(jīng)驗(yàn)，確定呷爬滑坡可被侵蝕的細(xì)顆粒最小粒徑為1 mm?；麦w的滲透系數(shù)約為0.3 mm/s，孔隙率約為0.3，故發(fā)生內(nèi)部侵蝕的臨界水力坡降約為0.398。

圖1為庫水位線處的水力梯度變化圖，從圖中可以看出，隨著庫水位的變化，水位線處的水力梯度大致在0.52～0.65范圍內(nèi)變化，大于臨界水力坡降，所以在庫水位以下的滑坡體內(nèi)會發(fā)生內(nèi)侵蝕現(xiàn)象。為了模擬出滑坡體在內(nèi)侵蝕過程中土體的孔隙率、被侵蝕顆粒的質(zhì)量、細(xì)顆粒流失速率等，需要建立相應(yīng)的內(nèi)侵蝕控制方程。

圖1 水力梯度變化圖

2 內(nèi)侵蝕控制方程

2.1 內(nèi)侵蝕下的質(zhì)量守恒方程

在內(nèi)侵蝕過程中，細(xì)顆粒從土體骨架脫離出來并隨水運(yùn)動。假設(shè)其運(yùn)動速度與水的流速一致，Cividini等[13]利用質(zhì)量守恒定律得到內(nèi)侵蝕作用的連續(xù)性方程：

(2)

式中：ρtr為被侵蝕的細(xì)顆粒的密度；vi為被侵蝕顆粒在i方向上的移動速度；qe為單位體積細(xì)顆粒被侵蝕的速率；qdp為被侵蝕顆粒沉積到土骨架的速率。

Vardoulakis等[14]將飽和多孔介質(zhì)單元體分為未被侵蝕的土體、水和被侵蝕的細(xì)顆粒三部分,質(zhì)量分別為dMs，dMw，dMfs，體積分別為dVs，dVw，dVfs。當(dāng)有細(xì)顆粒被侵蝕時(shí)，細(xì)顆粒隨水一起流動，被侵蝕的顆粒和水組成的流動體密度為

(3)

式中：ρs為土顆粒的密度；ρw為水的密度；c為被侵蝕細(xì)顆粒的體積分?jǐn)?shù)。被侵蝕顆粒的密度為

(4)

將式(4)代入(2)，則內(nèi)侵蝕作用的連續(xù)性控制方程可以改寫為

(5)

2.2 內(nèi)侵蝕下的滲流方程

發(fā)生內(nèi)侵蝕時(shí)，土體內(nèi)滲流狀態(tài)十分復(fù)雜，達(dá)西流、非達(dá)西流、優(yōu)勢流可同時(shí)存在，且描述土體滲流狀態(tài)的方程也不盡相同。但在內(nèi)侵蝕發(fā)生的初期，土體中的滲流以達(dá)西流為主，故本文采用達(dá)西定律來表示內(nèi)侵蝕下的滲流方程，即

Vw=k▽p

(6)

Kozeny-Carman公式是描述多孔介質(zhì)的滲透率與孔隙率關(guān)系精度最高的經(jīng)驗(yàn)公式，其表達(dá)式為

(7)

式中：η為Kozeny常數(shù)，將初始滲透率k0和初始孔隙率φ0代入上式，可得

(8)

則式(5)可以寫為

(9)

隨著內(nèi)侵蝕的發(fā)生，流速增大，雷諾數(shù)隨之增大，流體的滲流速度與水力梯度之間將不再呈線性關(guān)系，達(dá)西定律不再適用，可改用福希海默方程[15]等描述非達(dá)西流的滲流方程。

2.3 內(nèi)侵蝕速率

Cividini等[16]根據(jù)對重塑土樣的室內(nèi)試驗(yàn)，認(rèn)為在一定水力梯度下一直對土體進(jìn)行內(nèi)侵蝕作用，最終也不會將全部的細(xì)顆粒帶走，而是會有細(xì)顆粒殘留在土體內(nèi)。在滲透作用下，土體內(nèi)部的細(xì)顆粒被侵蝕，被侵蝕的顆粒粒徑由小變大。所以設(shè)土體的骨架可被侵蝕的最大顆粒為rmax，可被侵蝕的顆粒最小粒徑為rmin。如果rrmin時(shí)，則細(xì)顆粒有可能被侵蝕出來。所以，本文考慮在呷爬滑坡土體級配曲線上，只有顆粒粒徑在rmin和rmax之間的顆粒才考慮內(nèi)侵蝕速率。

(10)

式中，土體的體積變化率為

(11)

(12)

式中：φ為侵蝕后土體孔隙率，

φ=dVfs/dVs+dVw

3 呷爬滑坡內(nèi)侵蝕特性分析

3.1 基本情況

錦屏呷爬滑坡位于錦屏一級水電站壩址上游11.5 km，前緣高程1.655 km，與枯期河水位持平，后緣高程2.120 km，滑坡縱長約880 m，寬約260～300 m，滑體平均厚度約60 m，體積約1.3×107m3?；w物質(zhì)主要為塊碎石土，滑帶厚1～8 m，主要由灰黑色泥夾碎石組成[17]。

前期勘察表明：呷爬滑坡除水庫蓄水至正常蓄水位(1.880 km高程)或迭加VII度地震時(shí)，滑體有滑動的可能，其余工況均能保持穩(wěn)定或基本穩(wěn)定，根據(jù)滑體物質(zhì)組成及現(xiàn)今地形特點(diǎn)分析，預(yù)計(jì)失穩(wěn)方式將以逐級牽引式滑塌為主。

水庫蓄水期，呷爬滑坡前緣見小規(guī)模次級滑塌，整體未見明顯變形破壞跡象。2015年11月水庫第2次蓄水至1.880 km后調(diào)查，發(fā)現(xiàn)呷爬滑坡產(chǎn)生整體蠕滑，在高程約1.900、2.050 km及后緣地帶出現(xiàn)貫通性拉裂縫。直至2017年6月庫水位消落至1.800 km后調(diào)查，發(fā)現(xiàn)呷爬滑坡仍在整體滑移，坡體原有拉張裂縫的下錯幅度有明顯增加，其中，最后緣貫通裂縫最大下座高度已大于200 cm，較2016年1.800 km水位調(diào)查時(shí)增加約50～70 cm，表明呷爬滑坡的變形仍在增加。該呷爬滑坡的地形圖和地質(zhì)剖面圖如圖2和圖3所示。

圖3 呷爬滑坡地質(zhì)剖面圖

3.2 內(nèi)侵蝕特性分析

本文利用COMSOL Multiphysics軟件建立有限元模型，求解滲流和侵蝕耦合控制方程組，即式(5)、(6)和(12)。

根據(jù)圖2的滑坡地形圖和圖3的地質(zhì)剖面圖,忽略滑坡的中上部裂縫與后緣裂縫，建立數(shù)值分析的模型。在數(shù)值模型中，認(rèn)為水為不可壓縮流體，固體顆粒的壓縮量極小可以忽略，并且不考慮流體中含有細(xì)顆粒時(shí)對黏滯系數(shù)的影響，根據(jù)實(shí)驗(yàn)結(jié)果和級配曲線，各參數(shù)取值如下:ρw=1 000 kg/m3,Es=80 MPa,μw=1.0 mPa·s,νs=0.35,φbr=0.2,φs=0.348 3,Kbr=0.01 mm/s,ρs=2 100 kg/m3,k0=0.4 mm/s,rmin=1 mm/s,φ0=0.3,rmax=6 mm。

整個(gè)模擬試驗(yàn)過程持續(xù)525 d，水頭邊界條件根據(jù)實(shí)測的變化情況施加在模型邊界上?；麦w表面為透水邊界，滑帶為不透水邊界，初始水頭邊界條件根據(jù)實(shí)測水位取值；孔隙率的初始條件設(shè)置為φ0；試驗(yàn)開始時(shí)的水頭與試樣的上邊界一致，初始速度場v為0；在考慮重力的情況下，模型內(nèi)部的初始壓力為重力產(chǎn)生的壓力場，即p0=ρwgz。

根據(jù)數(shù)值模擬的結(jié)果可以得到滑坡體在庫水位作用下的滲透率k、被侵蝕顆粒的體積分?jǐn)?shù)φ以及內(nèi)侵蝕速率qe的變化情況。

圖4為1.80 km和1.84 km高程處在庫水位變化情況下，滑坡體內(nèi)滲透率k的變化情況。1.80 km是測量周期內(nèi)的最低水位，此處始終受到內(nèi)侵蝕的作用。從圖中可以看出，在內(nèi)侵蝕持續(xù)作用的整個(gè)過程中，不斷有細(xì)顆粒被侵蝕，滑坡體的滲透率一直在增大且逐漸趨于平穩(wěn)。內(nèi)侵蝕發(fā)生的初期0～200 d，滲透率的斜率較大，表明有大量粒徑大于rmin的可被侵蝕顆粒被水流沖刷帶走；中期200～450 d內(nèi)，由于小粒徑的顆粒逐漸變少以及顆粒在滲流通道內(nèi)的累積，滲透率的變化趨勢變緩；450 d以后，大部分可被侵蝕顆粒都被侵蝕掉了，但仍有少量粒徑接近rmax的顆粒被侵蝕，此時(shí)內(nèi)侵蝕作用效果逐漸減弱，滲透率也趨于平穩(wěn)，最終將保持不變。在1.84 km高程處的滲透率與1.80 km處的滲透率變化規(guī)律相似，但在75～165 d及405 d以后，庫水位在1.84 km以下，所以此時(shí)高程1.84 km處不受庫水位的作用，內(nèi)侵蝕作用暫時(shí)停止，被侵蝕顆粒的量為0，在圖中表現(xiàn)為滲透率曲線保持水平狀態(tài)。當(dāng)庫水位再次超過1.84 km時(shí)，內(nèi)侵蝕又再次發(fā)生。由于越往下的部分受到庫水作用的時(shí)間就越長，受到內(nèi)侵蝕的時(shí)間越來越長，滲透率越大，故1.80 km高程處的滲透率比1.84 km高程處的滲透率大。

圖4 滑坡體內(nèi)滲透率的變化

圖5為不同高程下滑坡體內(nèi)可被侵蝕顆粒的侵蝕量占土體的體積分?jǐn)?shù)φ，與滲透率的變化規(guī)律基本一致。1.80 km高程處完全處于水下，受到的內(nèi)侵蝕時(shí)間最長，被侵蝕的細(xì)顆粒越多?？杀磺治g顆粒占土體的體積分?jǐn)?shù)的定義與孔隙率的定義類似，對于顆粒粒徑小于rmax的顆?？梢钥醋髟趦?nèi)侵蝕發(fā)生之前就已經(jīng)被侵蝕掉了，此時(shí)的可被侵蝕顆粒的體積分?jǐn)?shù)相當(dāng)于土體的初始孔隙率0.3，對于粒徑處于rmin和rmax之間的顆粒，則需要在庫水的作用下被侵蝕掉。從圖中的趨勢線可以看出，如果在525 d之后庫水位超過1.84 km，那么1.84 km高程處的可被侵蝕顆粒的φ將沿著原來的增長趨勢繼續(xù)增加，最終將和1.80 km高程處的φ值趨于同一個(gè)穩(wěn)定的φmax。當(dāng)所有可被侵蝕的顆粒全部被侵蝕完之后，內(nèi)侵蝕結(jié)束。

圖5 滑坡體內(nèi)可被侵蝕顆粒體積分?jǐn)?shù)的變化

圖6為庫水位對內(nèi)侵蝕速率qe的影響規(guī)律。從圖中數(shù)值模擬結(jié)果可以看出，由于初始孔隙率的存在，臨界粒徑rmin的顆粒很容易通過土體內(nèi)部原有的孔隙被水流沖刷帶出，所以在內(nèi)侵蝕發(fā)生的初期的0～30 d內(nèi)，內(nèi)侵蝕速率是增加的。從總的趨勢來看，隨著細(xì)顆粒不斷被侵蝕，可被侵蝕的顆粒越來越少，內(nèi)侵蝕速率是持續(xù)下降的，直到所有可被侵蝕的顆粒被侵蝕完之后，內(nèi)侵蝕速率變?yōu)?，內(nèi)侵蝕過程結(jié)束。其中，在75～165 d及405 d以后，1.84 km高程處的水位是低于1.84 km的，此時(shí)的內(nèi)侵蝕速率為0，滲透率k和可被侵蝕顆粒的φ保持不變，這與圖4和5中表現(xiàn)的規(guī)律是一致的。本文所考慮的內(nèi)侵蝕速率是由可被侵蝕但仍未被侵蝕的顆粒含量決定的，同一時(shí)間，高程越低的地方受庫水影響的時(shí)間就越長，被侵蝕顆粒的量越多，可被侵蝕但還未被侵蝕的顆粒就越少，導(dǎo)致內(nèi)侵蝕速率越來越小，故在圖6中體現(xiàn)出的是在同一時(shí)間下，絕對水位是一定的，高程高的地方，由于受庫水的影響時(shí)間少，未被侵蝕的顆粒含量仍很高，故內(nèi)侵蝕速率較高程低的地方大。

圖6 內(nèi)侵蝕速率的變化

4 結(jié) 論

本文對庫水位作用下土體內(nèi)侵蝕特性進(jìn)行了研究，總結(jié)如下：

(1) 根據(jù)質(zhì)量守恒方程，建立連續(xù)性方程、滲流方程和內(nèi)侵蝕速率。初步考慮土體的顆粒級配對內(nèi)侵蝕影響，以可被侵蝕的最小顆粒粒徑rmin和可被侵蝕的最大顆粒粒徑rmax來控制受內(nèi)侵蝕作用影響的范圍。當(dāng)rrmin時(shí)，則細(xì)顆粒有可能被侵蝕出來。

(2) 以錦屏呷爬滑坡滑體為研究對象，利用COMSOL Multiphysics軟件建立有限元模型，耦合上述土體內(nèi)侵蝕過程中被侵蝕顆粒的質(zhì)量守恒方程、內(nèi)侵蝕的滲流方程和考慮級配影響的內(nèi)侵蝕速率。模擬結(jié)果表明土體孔隙率和被侵蝕顆粒的體積分?jǐn)?shù)在內(nèi)侵蝕作用下不斷增加，當(dāng)所有可被侵蝕的顆粒都被沖刷掉之后，細(xì)顆粒的流失速率最終趨于穩(wěn)定。而且某一高程處的內(nèi)侵蝕作用會在水位低于此處高程時(shí)暫時(shí)停止，當(dāng)水位重新超過此處高程時(shí)，內(nèi)侵蝕作用會繼續(xù)沿著原來的趨勢進(jìn)行下去。但是該模型未考慮非達(dá)西流情況下，即滲透速率與水力梯度呈現(xiàn)非線性增長時(shí)的變化情況，以及未考慮水位變化速率的影響。

(3) 該呷爬滑坡的前緣長期處于水下，數(shù)值模擬結(jié)果顯示處于水下的部分始終受到內(nèi)侵蝕的影響，細(xì)顆粒在不斷的流失，孔隙率持續(xù)增大，直至滑坡失穩(wěn)。因此本文的研究可以用于大致確定滑坡體正處于內(nèi)侵蝕的哪個(gè)階段和水下部分的細(xì)顆粒流失程度，對實(shí)時(shí)監(jiān)測滑坡的穩(wěn)定性有一定的實(shí)際意義。