基于隨機裂隙的長距離調水工程渠道邊坡穩定性模擬

王棟 方攀博 韓延成 黃兆虎 鄧曉川 趙嘉誠

關鍵詞：裂隙巖體；隨機裂隙網絡；流固耦合；強度折減法；邊坡穩定性

很多長距離調水工程如南水北調中線和東線、引黃濟青工程、膠東調水工程等采用明渠形式輸水。明渠經過裂隙發育地區時，裂隙和孔隙會對明渠邊坡內部的滲流場、位移場產生重要的影響，導致渠道產生滲透破壞，使邊坡穩定性降低，嚴重影響明渠的輸水能力和安全。

引起邊坡失穩的原因有很多，主要原因有降雨入滲、地下水位變化、邊坡地質構造等。降雨因素特別是歷時長、強度大的暴雨入滲到邊坡內部，影響土體骨架的力學性質，使基質吸力下降，抗剪強度降低，從而引起邊坡失穩［１］。地下水位的變化會影響土體強度，使坡面產生不穩定變化［２］。馬慧敏等［３］研究了工程特性較特殊的膨脹土（巖）渠段的邊坡穩定性，并給出了分析邊坡穩定性時參數取值的經驗方法。很多調水工程邊坡存在裂隙且現有的邊坡穩定性分析軟件不具備模擬巖體中隨機裂隙網絡的功能，使得模擬結果并不能完全反映實際情況下渠道的邊坡穩定性情況，故模擬含裂隙邊坡的穩定性較為困難。常峻嶺［４］根據蒙特卡羅（Ｍｏｎｔｅ－Ｃａｒｌｏ）法模擬隨機裂隙網絡，研究了幾何參數、力學參數對邊坡穩定性的影響。鄧東平等［５］研究了張裂縫位置及張裂縫水位對邊坡穩定性的影響。目前邊坡穩定性分析方法主要有極限平衡法和有限元強度折減法，因極限平衡法并不能求出變形情況，且不能反映實際條件下邊坡所受內力，還需要假設滑動面位置及形狀，故較多學者采用有限元強度折減法進行邊坡穩定性分析［６－９］?？撞环驳龋郏保埃莺投鹩竦龋郏保保輰㈦x散單元法和強度折減法結合模擬邊坡穩定性，同樣取得了較好的成果。

近年來因巖體裂隙而造成渠道邊坡破壞的事故較多，裂隙滲透對渠道邊坡穩定性的影響受到廣泛關注。筆者以山東省膠東調水工程萊州段輸水明渠為例，基于蒙特卡洛原理，生成符合實際巖體裂隙分布的二維隨機裂隙網絡，將達西定律模塊和固體力學模塊進行耦合，建立了渠道邊坡二維隨機裂隙模型，采用強度折減法對比研究了兩種降雨條件（歷時長而強度弱、歷時短而強度強）和改變原始地下水位條件下兩種介質模型（裂隙模型和孔隙模型）渠道的邊坡穩定性。

１研究區概況

山東省膠東調水工程萊州段輸水明渠位于萊州市虎頭崖鎮西北部的秀東村至后趴埠村。渠道邊坡巖性主要為花崗巖，受巖體邊坡中裂隙影響，強降雨期間渠道的地下水位變幅較大。

２隨機裂隙網絡生成

通過水文地質勘測資料及渠道兩岸鉆孔調查，對裂隙各參數分布規律進行計算分析，將研究區的裂隙分為兩組，得到各組裂隙的初始幾何參數，見表１?；诿商乜寤驹恚郏保玻?，利用Ｍａｔｌａｂ編制二維離散裂隙網絡生成程序，生成的裂隙網絡模型（如圖１所示）尺寸為６０ｍ×３０ｍ（長×深）。山東省膠東調水工程萊州段（樁號６２＋８０９）輸水明渠邊坡單側斷面如圖２所示。處理后最終生成隨機裂隙分布圖（見圖３）。

３模型構建及原理

３．１滲流－應力耦合模型

３．１．１滲流微分方程

由于降雨對滲流及邊坡穩定性具有很大的影響，因此滲流計算采用非飽和滲流Ｒｉｃｈａｒｄｓ方程［１３］：

３．１．２滲流－應力耦合原理

１）滲流場對應力場的影響。巖體中滲透壓力對巖體的變形有重要影響，動水壓力導致裂隙開度增大，巖體變形影響巖體的穩定性［１５］。Ｌｏｕｉｓ［１６］通過分析測試認為滲透壓力對裂隙巖體的作用由３種體積力構成，即在一組平行節理中由水的黏性產生的切向力、靜水壓力或阿基米德壓力、動水壓力或滲透力。

２）應力場對滲流場的影響。當巖體受滲透體積力和其他外部荷載作用時，巖體中應力場發生改變，從而巖體裂隙寬度改變，進而影響滲透系數和巖體滲流場。根據達西定律，流體的滲透系數可表示為［１５］

可知，影響巖土體滲透性能的因素主要為流體性質和巖土體骨架性質［１７］。

３．２邊界條件及相關參數

在滲流模型中，上邊界和坡面上層邊界為大氣邊界（在降雨條件下為降雨邊界）；坡面下層和坡底位置即渠道襯砌部分實測滲漏量非常小，可以忽略不計；渠道內無水時渠道邊坡內地下水水位與渠道內水位之差最大，揚壓力最大，為最不利工況，因此渠道底部概化為無流動邊界；模型左側邊界為水頭邊界即初始地下水位，右側按年庚乾等［１８］研究的滲流邊界處理，滲流邊界方程如下：

在力學模型中，坡面上層邊界、渠道坡面和渠底邊界均為自由邊界；模型左側邊界和右側邊界均設置為輥支撐邊界（橫向位移為０）；坡體底部邊界為固定邊界。對于滲流邊界，當地下水位上升或者土體達到飽和狀態時，水從邊坡內不斷滲出。邊界設置如圖４所示（黃色為滲流邊界條件、綠色為力學邊界條件）。

在滲流計算模塊中飽和－非飽和滲流計算采用達西定律匹配理查茲方程作為基本方程，裂隙模型和孔隙模型除了飽和滲透系數不同外，其他采用相同的水力參數。非飽和水力參數采用ＶａｎＧｅｎｕｃｈｔｅｎ模型計算。固體力學模塊計算模型采用線彈性本構模型搭配土壤塑性模型，采用Ｄｒｕｃｋｅｒ－Ｐｒａｇｅｒ準則匹配Ｍｏｈｒ－Ｃｏｕｌｏｍｂ準則作為破壞準則。相關計算參數由地勘結果、監測井抽水試驗、室內力學試驗得到，見表２、表３。

３．３強度折減法失穩判據

目前，用強度折減法來判斷巖土邊坡的破壞失穩判據主要有３種：１）有限元計算結果不收斂，說明此時邊坡已經發生破壞，相應的應力應變無法確定，則前一個輸入的折減系數對應的狀態為此邊坡的極限平衡狀態［１９］；２）如果塑性變形區域貫通了整個土層，則表示該區域即將失穩，在塑性變形貫穿前的強度折減系數即為邊坡的安全系數［２０］；３）邊坡位移改變率發生突變時折減系數就是邊坡的安全系數［２１］。

４結果分析

４．１模擬方案

共設置１２種模擬工況，見表４。

４．２滲流分析

本文用Ｃｏｍｓｏｌ有限元數值計算軟件分別進行基于孔隙模型和裂隙模型的非飽和滲流模擬，得到任意時刻的滲流場。圖５和圖６為降雨強度為２０ｍｍ／ｈ和５ｍｍ／ｈ條件下７０ｈ和２８０ｈ的滲流場（云圖代表孔隙水壓力，面上流線代表達西速度場）。

滲流模擬結果表明，裂隙模型的入滲深度大于孔隙模型的，裂隙模型的體積含水率小于孔隙模型的。這是因為雨水進入裂隙后迅速下滲，入滲深度迅速增大，裂隙內的入滲水量下滲至地下水，所以入滲區內體積含水率較小，孔隙水壓力減小。裂隙模型的入滲率大于孔隙模型的，隨著降雨時長的增長，兩種介質模型的入滲率分布趨近于穩定值。

４．３不同降雨條件邊坡穩定性分析

根據上述滲流模擬結果，進行不同介質模型邊坡塑性變形和位移數值模擬結果對比。

４．３．１塑性變形

通過數值模擬得到不同降雨條件邊坡塑性變形的情況?？芍?，在２０ｍｍ／ｈ的降雨強度下，裂隙模型最大折減系數為１．４８，最大塑性變形為３．０１ｃｍ；孔隙模型最大折減系數為１．５０，最大塑性變形為１．７１ｃｍ。５ｍｍ／ｈ降雨強度下，裂隙模型最大折減系數為１．５０，最大塑性變形為１．６２ｃｍ；孔隙模型最大折減系數為１．５２，最大塑性變形為１．２４ｃｍ。由此可得：隨著折減系數的增大，塑性變形區域從坡腳逐漸向坡頂移動，折減系數最大時從坡腳至坡頂形成了貫通區域，最大塑性應變發生在坡腳處?？紫赌Ｐ偷陌踩禂荡笥诹严赌Ｐ偷?，且塑性變形小于裂隙模型的。

４．３．２位移

模擬得到在極限平衡狀態不同降雨強度下裂隙和孔隙模型的位移分布云圖，如圖７、圖８所示。由圖７、圖８可知，裂隙模型和孔隙模型得到的位移場大致相同，邊坡整體的滑動趨勢都是從坡頂逐漸向坡腳移動且坡腳處位移最大。２０ｍｍ／ｈ降雨強度下，裂隙模型坡腳處最大位移為７．４３ｃｍ，安全系數為１．４８；孔隙模型坡腳處最大位移為３．２３ｃｍ，安全系數為１．５０。５ｍｍ／ｈ降雨強度下，裂隙模型坡腳處最大位移為３ｃｍ，安全系數為１．５０；孔隙模型坡腳處最大位移為２．２９ｃｍ，安全系數為１．５２?？v向對比可知，雨強較小時邊坡模型的安全系數較大，邊坡更不易發生失穩；橫向對比可得，相同降雨條件下，裂隙模型的邊坡安全系數始終小于孔隙模型的，孔隙模型邊坡更穩定。

４．３．３最大位移與安全系數的關系

基于孔隙模型和裂隙模型模擬得到不同降雨條件下邊坡最大位移與安全系數的關系，如圖９所示。由圖９可知，邊坡最大位移隨安全系數的增大而增大，且增長趨勢大致相同。當安全系數小于１．４６時，降雨強度和模型介質是否存在裂隙對邊坡的最大位移影響很小，兩種介質模型的最大位移幾乎相同。安全系數為１．４６～１．４８時，最大位移增長率從大到小對應的降雨強度工況是裂隙模型２０ｍｍ／ｈ、裂隙模型５ｍｍ／ｈ、孔隙模型２０ｍｍ／ｈ、孔隙模型５ｍｍ／ｈ。安全系數為１．４８～１．５０時，５ｍｍ／ｈ降雨強度下，裂隙模型位移增長率大于孔隙模型的。綜上，隨著安全系數的增大，裂隙模型的位移增長率和最大位移均大于孔隙模型的，裂隙模型邊坡較孔隙模型更易破壞且失穩時位移更大。

４．４不同地下水位邊坡穩定性分析

４．４．１塑性變形

通過數值模擬得到不同地下水位下邊坡塑性變形情況，數據見表５。由表５可知，兩種介質模型邊坡均形成了明顯的貫通坡腳和坡頂的塑性變形區域，且裂隙模型塑性變形區域大于孔隙模型的。最大塑性變形發生在高地下水位條件下，裂隙模型和孔隙模型對應的最大值分別為５．１４、４．４１ｃｍ。隨著地下水位的降低，折減系數逐漸增大，坡腳處的最大塑性變形區域減小。裂隙模型的邊坡塑性變形均大于孔隙模型的。

４．４．２位移

模擬得到在極限平衡狀態下不同地下水位裂隙和孔隙模型的位移等值線分布云圖?？芍?，位移較大的部分從坡頂一直貫通至坡腳處。隨著地下水位的降低，兩種介質模型的安全系數逐漸增大，邊坡穩定性逐漸增強，低地下水位條件下模型的邊坡穩定性最佳。在２０ｍｍ／ｈ的降雨強度下，裂隙模型和孔隙模型的邊坡安全系數都為１．５２，最大位移分別為７．４３、３．２３ｃｍ，裂隙模型位移較大的部分主要集中在邊坡中下部，失穩時有可能發生牽引式滑動，其原因是降雨強度較大，降雨入滲時，雨水幾乎全部通過裂隙向下滲透，入滲速度較快，入滲深度較大，坡腳處先達到飽和，從而使坡腳處的穩定性大幅減弱。在５ｍｍ／ｈ的降雨強度下，裂隙模型和孔隙模型邊坡安全系數都為１．５２，最大位移分別為１４、１１ｃｍ，是１２種工況中最大的。原因是降雨歷時太長且地下水位較高時，入滲水較快到達地下水位后上升至坡腳，致使坡腳處一定范圍內巖體浸水全部達到飽和，坡腳長時間浸于水中承受較大孔隙水壓力而產生較大的位移。裂隙模型安全系數始終小于等于孔隙模型的，在相同安全系數的前提下，裂隙模型的邊坡最大位移大于孔隙模型的，孔隙模型的邊坡穩定性更佳。相同地下水位、雨強較小時，邊坡安全系數更大、更穩定。

４．４．３最大位移與安全系數的關系

基于孔隙模型和裂隙模型模擬得到改變初始地下水位條件下邊坡最大位移與安全系數的關系，如圖１０、圖１１所示?？芍冢玻埃恚恚杞涤陱姸认?，位移增長率最大的是工況１，邊坡穩定性最好的是工況６，工況６安全系數最大且最大位移最?。辉冢担恚恚璧慕涤陱姸认?，位移增長率最大的是工況１１，邊坡穩定性最好的是工況４。發生最大位移的分別是工況１１和工況１２（原因已在４．４．２中解釋）。邊坡最大位移隨著安全系數的增大而增大，相同地下水位條件下裂隙模型邊坡最大位移均大于孔隙模型的。

５結論

本文基于裂隙模型和孔隙模型，采用強度折減法對比研究了不同降雨條件、不同地下水位條件下兩種模型對膠東調水工程萊州段輸水明渠渠道邊坡穩定性的影響，得到如下結論：

１）裂隙模型的安全系數小于等于孔隙模型的，且在相同安全系數的前提下，裂隙模型邊坡的塑性變形和最大位移大于孔隙模型的。裂隙模型在位移突變時的位移增長率大于孔隙模型的，表明裂隙模型的邊坡穩定性更差。

２）兩種模型模擬結果均表明，在相同降雨強度條件下，高地下水位時邊坡穩定性較差；在相同地下水位條件下，降雨強度較大時邊坡穩定性較差。

３）裂隙模型在高雨強條件下位移較大的部分主要集中在邊坡中下部，失穩時可能發生牽引式滑動。

因此，對裂隙發育地區的渠道邊坡穩定性分析，按常規忽略裂隙的孔隙模型計算，結果可能與實際不符，結果偏于安全，會給工程運行帶來潛在的安全風險。