基于有限差分法的馬馬崖下巖滑塌體穩定性分析

肖 蕾，駱世威，陳仁宏，路 雷

（中國電建集團貴陽勘測設計研究院有限公司，貴州 貴陽 550081）

1 工程概況

馬馬崖水電站為北盤江干流規劃梯級的第二級，屬2 等大（二）型工程，樞紐工程由碾壓混凝土重力壩、壩身開敞式溢流表孔、壩身放空底孔、左岸引水系統及左岸地下廠房等主要建筑物組成。

下巖滑塌體D2 位于馬馬崖壩址下游峽谷河段右岸，距馬馬崖壩址5.4 km。下巖滑塌體D2 分布高程490 m～1100 m，地形綜合坡角10°～16°，多處形成臺階狀地貌，分布面積3.88 km2。滑塌體下伏基巖為T2g1薄層、極薄層夾中厚層泥質白云巖、泥質灰巖夾泥巖、泥灰巖或T2g2下部的泥質灰巖、泥灰巖等，其后緣以上陡壁巖性主要為T2y 白云巖及T2g2灰巖夾泥質灰巖等，主要由孤、塊石、碎石夾粘土及滑塌層狀巖體等物質組成，鉛直深度一般60 m～110 m，總體積約3 億m3，為一個成因復雜、體積巨大的滑塌體。

根據地貌、物質組成和拉裂縫等的差異，從上往下分為Ⅰ、Ⅱ、Ⅲ三個區域。Ⅰ區為新近滑動區，Ⅱ區為古滑動區，Ⅲ區為崩塌堆積區，其中Ⅰ區下部有最新滑塌區、易滑塌區及最新開裂區。滑塌體地質平面圖及其分區見圖1。

圖1 滑塌體工程地質平面示意圖

2 研究思路及計算步驟

2.1 研究思路

對于大型邊坡而言，應力應變計算是其設計的重要內容之一。一般邊坡工程設計中，以剛體極限平衡法作為設計主要依據，以有限元、有限差分等數值計算方法為計算復核手段。本文結合該滑塌體的規模，進行有限差分數值分析，給出滑塌體的應力與變形、塑性區、拉力區、安全系數等內容。

邊坡失穩的判別主要是從以下三點進行考慮：1）塑性區貫通；2）計算不收斂判據；3）位移突變準則。對于下巖滑塌體，首先假定邊坡強度折減系數的范圍，計算分析在某個強度折減系數較大時的位移變形、塑性區及計算收斂情況，若判斷此強度折減系數下邊坡已發生塑性破壞，則減小強度折減系數進行計算，反之則增大強度折減系數，直至計算到邊坡處于臨界穩定狀態，此時的強度折減系數為邊坡穩定安全系數。

2.2 計算步驟

（1）計算天然地應力場

采用更改強度參數的彈塑性求解法求解邊坡天然自身應力場。首先將模型中材料的粘聚力和抗拉強度設為大值，計算至平衡后再將粘聚力和抗拉強度改為分析中采用的實際值計算至最終狀態，即可得到邊坡的天然應力場。

（2）利用強度折減法求解安全系數

對第（1）步計算結果的模型位移及初速度設置為0，然后將天然地應力場模型中材料更改為不同工況下相應強度折減系數對應的M-C 材料參數，進行計算達到設置的收斂參數為止。對計算結果進行位移、收斂性及塑性區分析，判斷該強度折減系數下邊坡是否穩定。分析多次強度折減系數下邊坡穩定情況，綜合確定邊坡穩定安全系數。

3 計算模型及計算參數

3.1 計算模型

選取典型剖面a-a' 為研究對象，根據計算工況分別建立計算模型。模型包含Ⅰ、Ⅱ、Ⅲ區，以X 軸為橫河向，指向岸坡為正，X 向總長3129 m；Z 軸為順河向，指向上游為正；Y 軸為鉛直向，鉛直向上為正。底部邊界Y 方向受Y 向約束，左右邊界X 方向受X 向約束，Z 向施加Z 向約束，作為平面應變問題考慮。

全部材料基本采用6 面體等參單元，局部采用三棱柱單元過渡。機組滿發工況計算模型見圖2，模型共包含20759 個單元。所取范圍內材料皆近似假定為各向同性、均勻連續彈塑性材料。

圖2 a-a' 剖面計算模型

3.2 計算參數

滑塌體、滑塌體下伏基巖、覆蓋層及基巖參數取值見表1～表3。

表1 滑塌體物理力學參數建議值表

表2 滑塌體下伏基巖物理力學參數建議值表

表3 覆蓋層及基巖計算參數

4 計算成果分析

4.1 天然地應力場

從圖3 可知：持久工況下滑塌體在自重荷載作用下大主應力基本呈層狀分布，大主應力等值線與坡面近于平行，在周圍邊界法向約束條件下，底部大主應力數值約為34 MPa；小主應力在滑塌體表面出現小范圍拉應力區。

圖3 滑塌體機組滿發水位下主應力云圖

4.2 滑塌體穩定性判別

對于機組滿發工況計算研究強度折減系數分別為0.90、0.95、0.96、0.97、0.98、0.99、1.0、1.01、1.02、1.03、1.04、1.05 時 滑塌體整體的變形、計算收斂及塑性區變化情況。

（1）位移判別

沿坡表選取21 個位移特征點，可反映整個滑坡體的變形情況。坡表特征點示意圖見圖4，特征點水平位移統計見圖5。

圖4 坡表特征點示意圖

圖5 滑塌體強度折減位移特征曲線

選取折減系數為0.99～1.01 時的位移特征曲線和位移大小綜合進行分析，安全系數為0.99～1.0 時，滑塌體下緣1～8 號點的水平向位移發生突變，初步判定滑塌體的安全系數在0.99～1.0之間。為增加安全系數精度，在0.98～1.0 之間選取8 組折減系數進行計算，進一步分析位移變化趨勢，計算結果見圖6。

圖6 強度折減位移特征曲線（0.98～1.0）

從圖6 可以看出，當強度折減系數在0.996～0.998 之間時，關鍵點水平向位移發生一定突變，因此滑塌體的安全系數在0.996～0.998 之間。

（2）不同安全系數下計算收斂情況

表4 不同強度折減系數下計算收斂情況

（3）塑性區貫通情況

由計算可知，Ⅱ、Ⅲ區穩定情況較好，在進行破壞模式分析時，僅以Ⅰ區為主要研究區域。選取強度折減系數為0.90、0.95、1.0、1.05 時的計算結果作為代表性研究對象，對其位移、塑性區、剪應變增量的變化趨勢進行分析，計算云圖見圖7～10。

圖7 F s=0.90 時計算云圖

圖8 F s=0.95 時計算云圖

圖9 F s=1.0 時計算云圖

圖10 F s=1.05 時計算云圖

由圖7～10 可知，機組滿發工況下滑塌體Ⅰ區的塑性剪切區主要分布在前緣局部、后緣拉裂縫及接觸帶區域。隨著強度折減系數增大，參數降低，塑性區沿著接觸帶從高高程向低高程逐步發育，后緣拉破壞區域逐漸增大，前緣剪切破壞范圍也增加，在河床底部塑性區達到貫通。結合塑性區、位移發展趨勢來看，滑塌體Ⅰ區接觸帶為潛在滑動面，存在沿接觸帶發生推移式破壞的可能，在推移式破壞發展過程中，滑塌體主要表現為后緣裂縫張開，前緣地形陡峭區域產生破壞。機組滿發時滑塌體的安全系數為0.996～0.998。

根據強度折減法計算成果，a-a'剖面安全系數統計見表5。通過對比可知，采用強度折減法的安全系數及滑動趨勢同剛體極限平衡法計算結果較為一致。

表5 a-a' 剖面安全系數匯總表

5 結論

（1）在各工況下，馬馬崖下巖滑塌體Ⅰ區穩定性較差，Ⅱ、Ⅲ區穩定性較好，不會發生整體破壞。機組滿發工況時，滑塌體Ⅰ區存在沿接觸帶發生推移式破壞的可能，在推移式破壞發展過程中，滑塌體主要表現為后緣裂縫張開，前緣地形陡峭區域產生破壞。其它各個工況下，塑性區范圍略有變化，但總體分布區域和發展趨勢與工況1 相近，說明水位變動不會影響滑塌體的破壞模式。

（2）強度折減法成果揭示的滑塌體安全系數及滑動趨勢與剛體極限平衡法較為一致，為實際工程設計提供了一定依據。