基于隨機(jī)裂隙的長(zhǎng)距離調(diào)水工程渠道邊坡穩(wěn)定性模擬

王棟 方攀博 韓延成 黃兆虎 鄧曉川 趙嘉誠(chéng)

關(guān)鍵詞：裂隙巖體；隨機(jī)裂隙網(wǎng)絡(luò)；流固耦合；強(qiáng)度折減法；邊坡穩(wěn)定性

很多長(zhǎng)距離調(diào)水工程如南水北調(diào)中線和東線、引黃濟(jì)青工程、膠東調(diào)水工程等采用明渠形式輸水。明渠經(jīng)過(guò)裂隙發(fā)育地區(qū)時(shí)，裂隙和孔隙會(huì)對(duì)明渠邊坡內(nèi)部的滲流場(chǎng)、位移場(chǎng)產(chǎn)生重要的影響，導(dǎo)致渠道產(chǎn)生滲透破壞，使邊坡穩(wěn)定性降低，嚴(yán)重影響明渠的輸水能力和安全。

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

近年來(lái)因巖體裂隙而造成渠道邊坡破壞的事故較多，裂隙滲透對(duì)渠道邊坡穩(wěn)定性的影響受到廣泛關(guān)注。筆者以山東省膠東調(diào)水工程萊州段輸水明渠為例，基于蒙特卡洛原理，生成符合實(shí)際巖體裂隙分布的二維隨機(jī)裂隙網(wǎng)絡(luò)，將達(dá)西定律模塊和固體力學(xué)模塊進(jìn)行耦合，建立了渠道邊坡二維隨機(jī)裂隙模型，采用強(qiáng)度折減法對(duì)比研究了兩種降雨條件（歷時(shí)長(zhǎng)而強(qiáng)度弱、歷時(shí)短而強(qiáng)度強(qiáng)）和改變?cè)嫉叵滤粭l件下兩種介質(zhì)模型（裂隙模型和孔隙模型）渠道的邊坡穩(wěn)定性。

１研究區(qū)概況

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

２隨機(jī)裂隙網(wǎng)絡(luò)生成

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

３模型構(gòu)建及原理

３．１滲流－應(yīng)力耦合模型

３．１．１滲流微分方程

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

３．１．２滲流－應(yīng)力耦合原理

１）滲流場(chǎng)對(duì)應(yīng)力場(chǎng)的影響。巖體中滲透壓力對(duì)巖體的變形有重要影響，動(dòng)水壓力導(dǎo)致裂隙開度增大，巖體變形影響巖體的穩(wěn)定性［１５］。Ｌｏｕｉｓ［１６］通過(guò)分析測(cè)試認(rèn)為滲透壓力對(duì)裂隙巖體的作用由３種體積力構(gòu)成，即在一組平行節(jié)理中由水的黏性產(chǎn)生的切向力、靜水壓力或阿基米德壓力、動(dòng)水壓力或滲透力。

２）應(yīng)力場(chǎng)對(duì)滲流場(chǎng)的影響。當(dāng)巖體受滲透體積力和其他外部荷載作用時(shí)，巖體中應(yīng)力場(chǎng)發(fā)生改變，從而巖體裂隙寬度改變，進(jìn)而影響滲透系數(shù)和巖體滲流場(chǎng)。根據(jù)達(dá)西定律，流體的滲透系數(shù)可表示為［１５］

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

３．２邊界條件及相關(guān)參數(shù)

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

在力學(xué)模型中，坡面上層邊界、渠道坡面和渠底邊界均為自由邊界；模型左側(cè)邊界和右側(cè)邊界均設(shè)置為輥支撐邊界（橫向位移為０）；坡體底部邊界為固定邊界。對(duì)于滲流邊界，當(dāng)?shù)叵滤簧仙蛘咄馏w達(dá)到飽和狀態(tài)時(shí)，水從邊坡內(nèi)不斷滲出。邊界設(shè)置如圖４所示（黃色為滲流邊界條件、綠色為力學(xué)邊界條件）。

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

３．３強(qiáng)度折減法失穩(wěn)判據(jù)

目前，用強(qiáng)度折減法來(lái)判斷巖土邊坡的破壞失穩(wěn)判據(jù)主要有３種：１）有限元計(jì)算結(jié)果不收斂，說(shuō)明此時(shí)邊坡已經(jīng)發(fā)生破壞，相應(yīng)的應(yīng)力應(yīng)變無(wú)法確定，則前一個(gè)輸入的折減系數(shù)對(duì)應(yīng)的狀態(tài)為此邊坡的極限平衡狀態(tài)［１９］；２）如果塑性變形區(qū)域貫通了整個(gè)土層，則表示該區(qū)域即將失穩(wěn)，在塑性變形貫穿前的強(qiáng)度折減系數(shù)即為邊坡的安全系數(shù)［２０］；３）邊坡位移改變率發(fā)生突變時(shí)折減系數(shù)就是邊坡的安全系數(shù)［２１］。

４結(jié)果分析

４．１模擬方案

共設(shè)置１２種模擬工況，見表４。

４．２滲流分析

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

滲流模擬結(jié)果表明，裂隙模型的入滲深度大于孔隙模型的，裂隙模型的體積含水率小于孔隙模型的。這是因?yàn)橛晁M(jìn)入裂隙后迅速下滲，入滲深度迅速增大，裂隙內(nèi)的入滲水量下滲至地下水，所以入滲區(qū)內(nèi)體積含水率較小，孔隙水壓力減小。裂隙模型的入滲率大于孔隙模型的，隨著降雨時(shí)長(zhǎng)的增長(zhǎng)，兩種介質(zhì)模型的入滲率分布趨近于穩(wěn)定值。

４．３不同降雨條件邊坡穩(wěn)定性分析

根據(jù)上述滲流模擬結(jié)果，進(jìn)行不同介質(zhì)模型邊坡塑性變形和位移數(shù)值模擬結(jié)果對(duì)比。

４．３．１塑性變形

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

４．３．２位移

模擬得到在極限平衡狀態(tài)不同降雨強(qiáng)度下裂隙和孔隙模型的位移分布云圖，如圖７、圖８所示。由圖７、圖８可知，裂隙模型和孔隙模型得到的位移場(chǎng)大致相同，邊坡整體的滑動(dòng)趨勢(shì)都是從坡頂逐漸向坡腳移動(dòng)且坡腳處位移最大。２０ｍｍ／ｈ降雨強(qiáng)度下，裂隙模型坡腳處最大位移為７．４３ｃｍ，安全系數(shù)為１．４８；孔隙模型坡腳處最大位移為３．２３ｃｍ，安全系數(shù)為１．５０。５ｍｍ／ｈ降雨強(qiáng)度下，裂隙模型坡腳處最大位移為３ｃｍ，安全系數(shù)為１．５０；孔隙模型坡腳處最大位移為２．２９ｃｍ，安全系數(shù)為１．５２?？v向?qū)Ρ瓤芍?，雨?qiáng)較小時(shí)邊坡模型的安全系數(shù)較大，邊坡更不易發(fā)生失穩(wěn)；橫向?qū)Ρ瓤傻茫嗤涤陾l件下，裂隙模型的邊坡安全系數(shù)始終小于孔隙模型的，孔隙模型邊坡更穩(wěn)定。

４．３．３最大位移與安全系數(shù)的關(guān)系

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

４．４不同地下水位邊坡穩(wěn)定性分析

４．４．１塑性變形

通過(guò)數(shù)值模擬得到不同地下水位下邊坡塑性變形情況，數(shù)據(jù)見表５。由表５可知，兩種介質(zhì)模型邊坡均形成了明顯的貫通坡腳和坡頂?shù)乃苄宰冃螀^(qū)域，且裂隙模型塑性變形區(qū)域大于孔隙模型的。最大塑性變形發(fā)生在高地下水位條件下，裂隙模型和孔隙模型對(duì)應(yīng)的最大值分別為５．１４、４．４１ｃｍ。隨著地下水位的降低，折減系數(shù)逐漸增大，坡腳處的最大塑性變形區(qū)域減小。裂隙模型的邊坡塑性變形均大于孔隙模型的。

４．４．２位移

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

４．４．３最大位移與安全系數(shù)的關(guān)系

基于孔隙模型和裂隙模型模擬得到改變初始地下水位條件下邊坡最大位移與安全系數(shù)的關(guān)系，如圖１０、圖１１所示?？芍冢玻埃恚恚杞涤陱?qiáng)度下，位移增長(zhǎng)率最大的是工況１，邊坡穩(wěn)定性最好的是工況６，工況６安全系數(shù)最大且最大位移最??；在５ｍｍ／ｈ的降雨強(qiáng)度下，位移增長(zhǎng)率最大的是工況１１，邊坡穩(wěn)定性最好的是工況４。發(fā)生最大位移的分別是工況１１和工況１２（原因已在４．４．２中解釋）。邊坡最大位移隨著安全系數(shù)的增大而增大，相同地下水位條件下裂隙模型邊坡最大位移均大于孔隙模型的。

５結(jié)論

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

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

２）兩種模型模擬結(jié)果均表明，在相同降雨強(qiáng)度條件下，高地下水位時(shí)邊坡穩(wěn)定性較差；在相同地下水位條件下，降雨強(qiáng)度較大時(shí)邊坡穩(wěn)定性較差。

３）裂隙模型在高雨強(qiáng)條件下位移較大的部分主要集中在邊坡中下部，失穩(wěn)時(shí)可能發(fā)生牽引式滑動(dòng)。

因此，對(duì)裂隙發(fā)育地區(qū)的渠道邊坡穩(wěn)定性分析，按常規(guī)忽略裂隙的孔隙模型計(jì)算，結(jié)果可能與實(shí)際不符，結(jié)果偏于安全，會(huì)給工程運(yùn)行帶來(lái)潛在的安全風(fēng)險(xiǎn)。