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

卿 菁， 朱 蕾， 陳 昊， 謝 鵬

(湖北工業大學土木建筑與環境學院， 湖北 武漢 430068)

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

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

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

1 滑坡地質特征

1.1 基本概況

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

1.2 物質組成及結構特征

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

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

滑床大部分表現為中風化狀態，主要由三疊系巴東組泥巖、砂巖、泥灰巖組成。產狀均為140°∠30°，通過勘查知道：巖石內部有正在發育的裂隙且有鈣質填充，巖石強度較高，內見泥巖條帶，巖石層面清晰。詳細剖面情況如圖2所示。

圖 2 盧家沱滑坡地質剖面圖

2 滑坡二維數值模擬分析

2.1 滑坡二維計算方法

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

根據盧家沱滑坡的地質結構，選取I-I′主剖面為計算剖面并建立出網格劃分模型，確定了單元數為1493，節點數為1544，如圖3所示。

圖 3 盧家沱滑坡網格圖

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

表1 盧家沱滑坡計算參數

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

表2 計算工況

2.2 滑坡工況對比分析

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

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

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

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

3 滑坡滲流場與應力場流固耦合分析

3.1 滑坡物理參數確定

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

表3 盧家沱滑坡有限元計算物理力學參數取值范圍

3.2 巖土體本構模型與屈服準則

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

3.3 邊界條件

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

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

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

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

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

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

圖 5 模型計算網格

3.5 滑坡工況計算

為更清楚的分析滑坡的應力應變等物理變化情況及規律模式，可以采用Abaqus軟件來進行不同工況組合下的流固耦合數值模擬分析。現取兩種工況進行詳細的分析，主要從這兩種相對特殊的工況來對比研究滑坡穩定性的狀態。

1)應力計算結果及分析

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

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

圖 6 工況1條件下主應力場分布 kPa

圖 7 工況2條件下主應力場分布 kPa

2)位移計算結果及分析

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

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

(a)工況1

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

3)塑性區計算結果及分析

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

圖9 塑性區分布圖

4)孔壓計算結果及分析

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

圖10 孔壓值顯示圖

5)滑坡整體穩定性

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

4 結論

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

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

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

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