庫水位波動與降雨作用下的滑坡流固耦合分析

卿 菁， 朱 蕾， 陳 昊， 謝 鵬

(湖北工業(yè)大學(xué)土木建筑與環(huán)境學(xué)院， 湖北 武漢 430068)

隨著我國工程建設(shè)的不斷發(fā)展，諸多地質(zhì)災(zāi)害問題也日益顯著。三峽庫區(qū)沿岸滑坡與水的影響密不可分，以庫水位升降、地下水作用和庫區(qū)降雨為主要影響因素的滑坡災(zāi)害多年來一直都是巖土工程界的重點研究對象。多種因素影響下的滑坡變形破壞是一個較為復(fù)雜的演變過程，這些過程里既存在有巖土體物理參數(shù)影響且改變周期相對較長的內(nèi)因，同時也存在水位變化、降雨、人工影響等綜合作用且改變周期相對較短的外因[1-3]。根據(jù)調(diào)查發(fā)現(xiàn)，絕大部分滑坡變形破壞起因都是因為地下水位的變化[4-7]。

在信息高速發(fā)展的今天，眾多學(xué)者廣泛采用了實驗分析、數(shù)值模擬等方法對滑坡水的作用進(jìn)行細(xì)致調(diào)查和剖析。例如，Rubin[8]針對二維滲流的相關(guān)問題，以有限差分法的手段求解了Richards方程的數(shù)值理論解；王恭先[9]對三峽地區(qū)滑坡的各種類型及其對應(yīng)的分布的特征行為進(jìn)行了歸納總結(jié)，同時對于滑坡的形成機(jī)理與穩(wěn)定性狀態(tài)進(jìn)行了分析；王思敬[10]根據(jù)水巖作用及水環(huán)境改變的類型將庫區(qū)滑坡分為水巖作用滑坡與天然形成的滑坡兩大類；而張桂榮[11]針對三峽庫區(qū)對某滑坡分析了水作用下的十種工況，并采用SEEP/W軟件模擬進(jìn)行分析。除此之外，三維有限元方法也日漸成熟，學(xué)者王媛[12]提出在Biot理論上，提出了一種假設(shè)孔隙水壓力與位移量來確定出滲流場并與應(yīng)力場耦合的方法。后來，陳慶中等[13]探究了滲流場計算方法，建立出新的流固耦合數(shù)學(xué)模型；劉建軍[14]針對地下水與地層結(jié)構(gòu)的相互作用關(guān)系，建立地下水流固模型來求解相關(guān)滲流問題。

現(xiàn)階段大多數(shù)的滑坡破壞機(jī)理研究只進(jìn)行了二維數(shù)值模擬或者是三維流固耦合分析，本文在前人已完成的眾多工作的基礎(chǔ)上，以秭歸縣盧家沱滑坡為例，結(jié)合二維分析有限元分析軟件Geo-Studio數(shù)值模擬分析和基于Abaqus建模的三維流固耦合分析，進(jìn)一步開展水作用對滑坡相關(guān)的影響行為研究，對該滑坡進(jìn)行破壞機(jī)理和致災(zāi)原因等方面的剖析，旨在為今后開展相關(guān)滑坡防治的研究工作者們提供進(jìn)一步的參考。

1 滑坡地質(zhì)特征

1.1 基本概況

盧家沱滑坡前緣高程135 m，后緣高程370 m，前緣175～195 m高程帶為一平臺，臨江處為坡度45°的陡坎，公路從滑體320 m高程處穿過。盧家沱滑左右邊界以基巖山脊為界，整體坡度約23°?；w縱向長約540 m，平均寬約280 m，均厚約22 m，面積約1.34×105m2，體積約2.95×106m3。具體滑坡平面情況詳見圖1。

1.2 物質(zhì)組成及結(jié)構(gòu)特征

盧家沱滑坡下部巖層組層多為泥巖，在庫水位升降及降雨等多方面因素影響下，易產(chǎn)生滑移變形，為動水壓力型滑坡?；w物質(zhì)分為兩部分：上部物質(zhì)為碎石土，褐黃色或灰褐色，碎石成分為灰?guī)r，粒徑一般為30～90 mm，最大250 mm，塊石呈中～強(qiáng)風(fēng)化狀態(tài)。土質(zhì)為砂質(zhì)粘土，土石分布不均勻，可塑，稍密-密實，土石比為6∶4～4∶6，表層土較松散，向下碎塊石含量增大。下部物質(zhì)為塊石土，塊石成分為泥質(zhì)灰?guī)r，塊徑大小不一，通常為30～120 mm，塊徑最大約為260 mm，土質(zhì)為土石比5∶5～3∶7的可塑性粉質(zhì)粘土層。

滑帶為角礫土，灰褐色，巖芯呈土柱狀，巖芯內(nèi)含角礫，礫徑為1～4 mm，土質(zhì)為石比為8:2的可塑性粉質(zhì)粘土層，地表未見滑帶露頭。

滑床大部分表現(xiàn)為中風(fēng)化狀態(tài)，主要由三疊系巴東組泥巖、砂巖、泥灰?guī)r組成。產(chǎn)狀均為140°∠30°，通過勘查知道：巖石內(nèi)部有正在發(fā)育的裂隙且有鈣質(zhì)填充，巖石強(qiáng)度較高，內(nèi)見泥巖條帶，巖石層面清晰。詳細(xì)剖面情況如圖2所示。

圖 2 盧家沱滑坡地質(zhì)剖面圖

2 滑坡二維數(shù)值模擬分析

2.1 滑坡二維計算方法

以Geo-Studio軟件中的SEEP/W模板對盧家沱滑坡二維模型進(jìn)行不同水位工況組合分析，進(jìn)而確定出對應(yīng)的滲流場分布情況；在得到滲流場的分布之后，再由SLOPE/W模板來施加滑體自重及庫水壓力大小，并結(jié)合Morgenstern-Price對滑體進(jìn)行穩(wěn)定性計算。

根據(jù)盧家沱滑坡的地質(zhì)結(jié)構(gòu)，選取I-I′主剖面為計算剖面并建立出網(wǎng)格劃分模型，確定了單元數(shù)為1493，節(jié)點數(shù)為1544，如圖3所示。

圖 3 盧家沱滑坡網(wǎng)格圖

依照滑坡地勘報告中試驗得到物理參數(shù)及相關(guān)建議值，本文確定出盧家沱滑坡計算參數(shù)如表1所示。

表1 盧家沱滑坡計算參數(shù)

針對盧家沱滑坡的計算工況考慮水位下降和降雨組合兩類，根據(jù)三峽庫區(qū)相關(guān)技術(shù)要求通常175 m水位下降至159 m速率為0.13 m/d，而}水位降至145 m時下降速率可在0.6～1.2 m/d中進(jìn)行調(diào)控，同時結(jié)合降雨條件，本文采用的計算工況如表2所示。

表2 計算工況

2.2 滑坡工況對比分析

根據(jù)SEEP/W分析出滲流分布情況后，用摩根斯坦-普瑞斯(Morgenstern-Price)法進(jìn)行計算，得出滑坡的穩(wěn)定性系數(shù)分布，并在此基礎(chǔ)上得出基于兩種工況組合的盧家沱滑坡穩(wěn)定系數(shù)。

從工況1的組合中可以看到：當(dāng)水位由175 m跌落至145 m的階段里，滑坡的整體穩(wěn)定性系數(shù)隨庫水位的跌落而持續(xù)降低，且在145 m水位線處達(dá)到最小值。同時，滑體的穩(wěn)定性也受到不同水位下降速率的影響，降速越大線性變化幅度越大，但從總體數(shù)據(jù)上看，穩(wěn)定系數(shù)值變化的范圍影響較小，在0.005內(nèi)波動。由此可看出，滑坡的穩(wěn)定性受庫水位降幅變化的影響不大。

從工況組合2中可以看出：在庫水位從175 m跌落的過程中，對施加降雨強(qiáng)度21.7 mm/d的工況進(jìn)行分析，當(dāng)在水位處于155～151.4 m區(qū)域時，施加降雨強(qiáng)度為55.33 mm/d的50年一遇暴雨進(jìn)行工況組合，直至145 m水位線處。受到不同水位降速的影響，盧家沱滑坡穩(wěn)定系數(shù)變化程度也均不相同，與工況1一樣，曲線的降幅與水位下降速率相關(guān)，降速越大降幅越大。當(dāng)水位線處于155～151.4 m區(qū)間中，由于施加50年一遇暴雨強(qiáng)度組合條件，曲線陡然下跌，待降雨停止后，曲線減小幅度明顯降低，由此可判斷滑坡穩(wěn)定性受到降雨因素的影響，且具有一定的滯后性。由整體曲線分析，穩(wěn)定系數(shù)值變化的范圍影響較小，在0.005內(nèi)?？傻贸鼋Y(jié)論：滑坡的穩(wěn)定性受庫水位降幅變化的影響不大。

工況組合2中疊加50年一遇三日暴雨，其穩(wěn)定性系數(shù)比不疊加降雨的穩(wěn)定性系數(shù)有所降低，可以確定出在暴雨條件下，滑坡的穩(wěn)定性系數(shù)會有一定程度的降低，但影響較小，且隨著水位降速增大的同時施加降雨，會明顯降低滑坡的整體穩(wěn)定性，所以當(dāng)水位降速為1.2 m/d時疊加強(qiáng)降雨為最不利組合情況，該情況下的滑坡穩(wěn)定性系數(shù)為1.057，說明盧家沱滑坡仍可處于基本穩(wěn)定狀態(tài)。

3 滑坡滲流場與應(yīng)力場流固耦合分析

3.1 滑坡物理參數(shù)確定

由相關(guān)勘察及監(jiān)測資料可以判斷，盧家沱滑坡的形成是由于地表特征、巖層結(jié)構(gòu)、地形地貌等內(nèi)因與水位升降、降雨、人工活動等外因共同作用下的結(jié)果，因此在數(shù)值模擬的過程中應(yīng)充分考慮相關(guān)因素影響下的物理參數(shù)取值，綜合分析，確定住有限元參數(shù)取值如表3所示。

表3 盧家沱滑坡有限元計算物理力學(xué)參數(shù)取值范圍

3.2 巖土體本構(gòu)模型與屈服準(zhǔn)則

滑坡體在受到水位升降、人工活動及降雨影響下，會發(fā)生各種程度的彈塑性變形，而在變形過程中會伴隨著張拉和壓剪破壞的產(chǎn)生。本文針對滑坡體在水位變化和降雨條件下的破壞模式，在考慮巖土體的屈服破壞狀態(tài)時，選擇采用拉破壞準(zhǔn)則及Mohr-Columb準(zhǔn)則共同來判斷。

3.3 邊界條件

根據(jù)飽和非恒定滲流與應(yīng)力耦合理論來針對盧家沱滑坡進(jìn)行流固耦合分析，所建立的邊坡邊界上的應(yīng)力條件及滲流條件如圖4所示。

圖 4 盧家沱滑坡邊界條件示意圖

3.4 計算工況及荷載組合的確定

為重點體現(xiàn)出水位波動與降雨條件下變化過程，本文僅采取普通工況和極端工況(表4中工況1、2)對滑坡進(jìn)行流固耦合數(shù)值模擬分析。

表4 計算及物理模型試驗工況及荷載組合

按照滑坡對應(yīng)的地形平面及剖面(圖1、圖2)條件，確定了滑體、滑帶及基巖的分布情況，建立出盧家沱滑坡三維模型尺寸如圖5所示，劃分出了4004個單元，以及12 418個相關(guān)節(jié)點數(shù)。

圖 5 模型計算網(wǎng)格

3.5 滑坡工況計算

為更清楚的分析滑坡的應(yīng)力應(yīng)變等物理變化情況及規(guī)律模式，可以采用Abaqus軟件來進(jìn)行不同工況組合下的流固耦合數(shù)值模擬分析?，F(xiàn)取兩種工況進(jìn)行詳細(xì)的分析，主要從這兩種相對特殊的工況來對比研究滑坡穩(wěn)定性的狀態(tài)。

1)應(yīng)力計算結(jié)果及分析

從圖6中工況1的應(yīng)力計算結(jié)果中可以發(fā)現(xiàn)：當(dāng)處于175 m水位時，拉應(yīng)力與壓應(yīng)力的分布區(qū)域較為明顯。對于滑體而言，最大主應(yīng)力位于滑體前緣部分，分布區(qū)間為-0.986～0.010 MPa；最小主應(yīng)力位于滑體后緣部分，分布區(qū)間為-2.286～0.014 MPa。對于滑帶而言，最大主應(yīng)力位于滑帶中后方區(qū)域，分布區(qū)間為-1.078～0.839 MPa；最小主應(yīng)力位于滑帶后部區(qū)域，分布區(qū)間為-2.386～0.123 MPa。

從圖7中工況2的應(yīng)力計算結(jié)果中可以發(fā)現(xiàn)：水位由162 m跌落至145 m同時遭遇50年一遇強(qiáng)降雨條件下，滑坡應(yīng)力場發(fā)生了明顯的改變。其中，前緣滑體主要受水位升降作用的影響較大，應(yīng)力大小發(fā)生了改變；中后緣滑體主要受降雨條件的影響較大，拉應(yīng)力與壓應(yīng)力的分布均發(fā)了較大的變化。對于滑體而言，在降雨停止之后，最大主應(yīng)力位于滑體中前緣部分，分布區(qū)間為-1.007～0.075 MPa；最小主應(yīng)力位于滑體后緣部分，分布區(qū)間為-2.341～0.021 MPa。對于滑帶而言，最大主應(yīng)力位于滑帶中部段區(qū)域，分布區(qū)間為-1.103～0.833 MPa；最小主應(yīng)力也位于滑帶中部段區(qū)域，分布區(qū)間為-2.537～0.077 MPa。由于應(yīng)力場與滲流場存在相互耦合的關(guān)系，因此導(dǎo)致應(yīng)力發(fā)生變化的重要影響因素就是水位升降聯(lián)合降雨作用促使滑體內(nèi)部滲流場發(fā)生了變化。

圖 6 工況1條件下主應(yīng)力場分布 kPa

圖 7 工況2條件下主應(yīng)力場分布 kPa

2)位移計算結(jié)果及分析

對于圖8中工況1的位移組合情況，可以得到：在庫水位保持175 m時，水平位移最大為2.568 cm，出現(xiàn)在滑體中后部范圍內(nèi)。其中，后緣滑體的水平位移區(qū)間為0.297～2.568 cm，中段滑體的水平位移區(qū)間為0.287～2.448 cm，前緣滑體的水平位移區(qū)間為0.153～0.183 cm。

對于圖8中工況2的位移組合情況，可以得到：最大水平位移位于滑體中后緣區(qū)域，大小為3.253 cm；在降雨條件停止后，后緣滑體的水平位移區(qū)間為0.266～3.189 cm，中段滑體的水平位移區(qū)間為0.271～3.253 cm；前緣滑體的水平位移區(qū)間為0.153～1.834 cm。經(jīng)上述結(jié)構(gòu)可以得出：相對于降雨作用，庫水位升降對滑坡的位移影響結(jié)果相對較小。隨著水位下降，水平位移增量先增加后處于減小趨勢，水平位移加速度逐漸減小，最后趨于穩(wěn)定，垂直位移增量變化以及垂直位移變化速率與水平方向位移變化規(guī)律一致，而垂直位移加速度先增大后不斷減小。

(a)工況1

(a)工況2圖 8 位移場分布 m

3)塑性區(qū)計算結(jié)果及分析

從圖9a中工況1塑性區(qū)分布情況可以看到滑體不存在明顯的塑性破壞范圍，基本可以確定當(dāng)蓄水至175 m水位時，滑坡沒有發(fā)現(xiàn)顯著的破壞跡象,一直處于穩(wěn)定狀態(tài)。從圖9b中工況2塑性區(qū)分布情況可以看到，當(dāng)水位從162 m降至145 m且施加50年一遇強(qiáng)降雨工況條件下，滑體不存在明顯的塑性破壞范圍，此時滑體沒有發(fā)現(xiàn)明顯的破壞跡象，處于基本穩(wěn)定狀態(tài)。

圖9 塑性區(qū)分布圖

4)孔壓計算結(jié)果及分析

由圖10a可知：175 m水位以下范圍內(nèi)，水位變化直接影響著孔壓的分布，前緣最大孔壓值大小為0.353 MPa，分布在滑體140 m高程處；地下水位對滑坡后緣孔壓影響較大，最大孔壓值為2.286 MPa。由圖10b可知：庫水位升降及降雨對滑坡前緣地下水的分布影響較大，當(dāng)水位區(qū)間在162 m至145 m跌落時，地下水位也相對降低且具有一定的滯后性特征；庫水位的升降作用對后緣滑坡的影響程度不大，但是受到降雨條件的影響，后緣孔壓值會有一定的增加。

圖10 孔壓值顯示圖

5)滑坡整體穩(wěn)定性

當(dāng)庫水位升降和大氣降雨聯(lián)合作用條件下，滑坡前緣至后緣均出現(xiàn)了一定程度的變形跡象，不同條件下的應(yīng)力分布及位移量均存在差異，局部變形較為明顯，可能出現(xiàn)局部滑動的情況，但在滑體內(nèi)部未出現(xiàn)大面積的破壞區(qū)域，說明盧家沱滑坡目前仍保持基本穩(wěn)定狀態(tài)。

4 結(jié)論

以盧家沱滑坡為研究對象，通過Geo-Studio軟件計算在不同水位、下降速率及降雨等組合工況條件下的滑坡穩(wěn)定性系數(shù)分布；同時根據(jù)Abaqus軟件模擬了盧家沱滑坡在庫水位變動聯(lián)合降雨組合工況下的滲流場、應(yīng)力場、位移場及塑性區(qū)等流固耦合及三維穩(wěn)定性情況，得到具體結(jié)論如下。

1)從二維穩(wěn)定性系數(shù)曲線圖可以看出:滑坡穩(wěn)定隨庫水位下降速率增大而不斷降低,疊加降雨條件情況下滑坡穩(wěn)定性出現(xiàn)最小值,滑坡前緣坍岸將會加劇,但由于滑坡前緣土體為中-強(qiáng)透水性,前緣內(nèi)的地下水滲出較為容易,滑坡穩(wěn)定性受水位下降影響不大,失穩(wěn)的可能性較小。

2)從三維數(shù)值模擬獲得的位移、應(yīng)力、孔壓及塑性區(qū)在各工況下的計算結(jié)果表明:在庫水位和降雨耦合作用下，滑坡會產(chǎn)生一定的變形,但整體仍處于基本穩(wěn)定狀態(tài),僅局部區(qū)域(滑體后部及滑體前緣)在庫水作用下變形相對較大,且同時考慮降雨影響,滑坡整體位移量相對較大。

3)由二維模型計算結(jié)果從宏觀方向上分析，可以得出在庫水位升降作用下，盧家沱滑坡的穩(wěn)定性逐漸降低，但水位下降速率對滑坡的影響不大，若聯(lián)合降雨條件作用下，對滑坡整體穩(wěn)定性影響相對較大；由三維流固耦合模型結(jié)果從微觀上分析，可以看出庫水位變化聯(lián)合降雨條件將會改變滑坡內(nèi)部的滲流場，導(dǎo)致滑坡發(fā)生較大的變形，若遭遇不利暴雨條件，盧家沱滑坡存在局部失穩(wěn)的可能。