滑坡抗剪強度參數反演數值模擬研究

劉光華，熊 超，趙 鵬

(1.重慶交通大學巖土工程研究所，重慶400074;2.重慶市地質災害防治工程勘查設計院，重慶400700)

滑坡是丘陵山區主要的地質災害類型及災害性地貌過程，具有強致災性和毀滅性的特點。丘陵山區滑坡災害頻繁發生，往往導致交通中斷、航道阻塞，嚴重威脅居民生命財產安全。

在進行滑坡穩定性評價及防治工程設計時，滑坡抗剪強度參數的選擇十分困難，而滑坡推力對c，φ值很敏感，強度參數發生微小的變化，其滑坡推力可能成倍增加，直接影響滑坡穩定性評價及治理工程設計的合理性、安全性、經濟性。目前c，φ值的確定常根據室內試驗或現場試驗，結合有關工程經驗數值，由極限平衡原理，建立計算公式進行反演求解。目前常采用以下3種方法：①假定剩余下滑力為0，安全系數K=1，滑坡處于臨界狀態，根據滑面土體條件結合類似土質資料或者室內試驗資料，定出其中一個指標反演另一個指標，求解強度參數;②在一次滑動中，找出兩個相鄰近的瞬間滑動計算斷面，列出兩個邊界條件相同的方程式，聯立求解;③根據同一斷面位置，不同時間但條件相似的兩次滑動瞬間計算斷面建立兩個反演式聯立求解［1-2］。極限平衡原理是將滑體假定為剛性體，相鄰條塊間只傳遞推力，因此，根據上述3種方法進行c，φ的反演有一定的局限性。

筆者運用FLAC3D有限差分的數值方法，建立計算模型，采用自編強度折減法進行穩定性求解，根據求解數值與滑坡所處穩定狀態進行比較，判斷其c，φ值，并與極限平衡法反演值進行比較，研究穩定系數對c，φ值的敏感性，綜合確定c，φ值，為滑坡抗剪強度參數反演提出新的方法。

1 滑坡穩定性分析的強度折減法

1.1 強度折減法的定義

根據庫倫公式，土體的抗剪強度τ=σtanφ+c，由于軟弱面上的c，φ值低于周圍土體，易于產生應力集中，因此，當軟弱面上某點或某部分的應力達到和超過其抗剪強度時，首先發生剪切破壞，該點或該部分稱為滑坡源。如果軟弱面的抗剪強度不斷降低，則進入剪切屈服的部位也不斷擴展，直到這些部位發展成為一條連貫的屈服面，形成整個滑動面［3］。將土體的強度參數c，φ分別乘以1/K(圖1)，并用于計算抗剪強度τ值，逐漸增加K使抗剪強度逐漸降低，直至失穩破壞。

圖1 抗剪強度與正應力的關系Fig.1 The relationship of normal stress and shear strength

1.2 強度折減法在FLAC3D中的實現

FLAC3D采用“二分法”［4］求解穩定系數。采用該法需定義解所在區間，合理選定區間值，可以減少計算時間。FLAC3D內置強度折減法，穩定系數的初始值分別定義為0和64，迭代計算復雜，運算時間長，不利于工程應用，因此采用自編強度折減法［5］，減少計算時間。其迭代過程分下述4個步驟：

1)假定K1，K2的值，取折減系數為 Ks=(K1+K2)/2，相應的強度為 c/Ks、tanφ/Ks，進行迭代計算;

2)若單次迭代計算收斂，但|K1－K2|＞ε，則K1=Ks、K2=K2帶入 1)中;

3)若單次迭代計算不收斂，則K2=Ks、K1=K1帶入1)中;

4)若迭代|K1－K2|＜ε時，計算截止，求得穩定系數Ks=(K1+K2)/2。

經過不斷的改變Ks，改變強度參數，進行循環計算，采用解的不收斂作為破壞標準，最終求得穩定系數。

1.3 穩定系數與強度參數反演的關系

所求穩定系數與滑坡體所處穩定狀態相互比較，若穩定系數與滑坡體所處穩定狀態相一致，則迭代計算中土體的c、φ值即為滑坡的強度參數;若所求穩定系數比滑坡體所處穩定狀態偏大，則減小c，φ值后進行計算;若所求穩定系數比滑坡體所處穩定狀態偏小，則增大c，φ值后進行計算。

2 工程應用

某路基填方邊坡，原地形坡度較陡，在填方施工中，在路基邊坡中下部設置擋土墻，坡表采用素混凝土格構護坡，但擋土墻未嵌入基巖。在暴雨作用下，表層土體呈飽水狀態，大量雨水滲入坡體，使土體的重度增加，強度參數降低，土體沿土巖交界面推動擋墻發生順坡向位移，導致后緣形成弧形裂縫，格構在土體的不均勻變形下發生斷裂，且前緣擋墻未嵌入基巖，擋土墻順滑坡體發生斷裂變形，失去支撐作用。為保證該路基滑坡強度參數的準確性，對該滑坡采用有限差分強度折減法，反演其強度參數，用于滑坡推力的計算。筆者只考慮自重荷載，取該路基滑坡的強變形段進行分析，由于該路基邊坡已經失穩，因此其穩定性系數取 0.95［6］。

2.1 計算模型的建立

根據鄭穎人，等［7］研究認為：在均質巖土體中，坡腳到左端邊界的距離為坡高的1.5倍，坡頂到右端邊界的距離為坡高的2.5倍，且上下邊界總高不低于2倍坡高，計算精度較為理想。但根據典型工程地質剖面圖，滑坡土體為土巖組合體(圖2)，其坡腳為土巖交界面，因此取坡腳到左端邊界的距離為10 m;由于路基邊坡剖面近似為梯形，因此取坡頂到右端邊界的距離為20 m;整體長度為65 m，橫向寬度取10 m。在ANSYS中建立三維地質模型，劃分地層單元及網格后生成符合FLAC3D的節點命令［8］。在FLAC3D生成的模型如圖3。

圖2 工程地質剖面Fig.2 Geological profile of engineering

圖3 滑坡有限差分模型Fig.3 FD model of landslide

2.2 材料本構模型及計算方法

巖土體的材料本構模型采用理想彈塑性模型，屈服準則采用摩爾-庫倫等面積圓屈服準則。由于土體沿基巖面發生滑動失穩，取內摩擦角的變化范圍為［10°，17°］，步長為 1，取黏聚力的變化范圍為［9 kPa，23 kPa］，步長為 2，進行敏感性分析，其余物理力學參數見表1。

表1 滑坡材料參數Table 1 Material parameters of landslide

3 計算結果與分析

采用自編強度折減法對不同的c，φ值進行數值分析，穩定性計算結果如表2。采用極限平衡法計算的穩定性計算結果如表3。

表2 FLAC3D穩定性計算結果Table 2 Stability calculation results of FLAC3D

表3 極限平衡法穩定性計算結果Table 3 Stability calculation results of LEM

對比表2、表3，極限平衡法計算的結果大于強度折減法，兩種方法所求的穩定系數的相對差值總體上逐漸增大，在1.145% ～7.705%之間。由于Fs=0.95，對應于穩定性計算表，在FLAC3D穩定性計算表中相應的c=21 kPa，φ=15°，而極限平衡法穩定性計算表中c=19 kPa，φ=15°。以此為基準值進行敏感性分析。敏感性系數選用如下公式：

式中：S為敏感系數;ΔX為某因素變化量;ΔFs為Fs對應ΔX的變化量;Fso為Fs的基準值;Xmax－Xmin為某因素最大變化量。

3.1 滑坡穩定性對φ值的敏感性分析

3.1.1 強度折減法計算的Fs對φ值的敏感性分析

根據表2的穩定性系數可知，Fs-φ的關系如圖4，Fs-φ 近似于線性變化，取 φ =15°，c=9 ～23 kPa進行敏感性分析，其敏感性計算結果如表4。φ對滑坡穩定性敏感系數為13.7% ～54.7%，平均值為29.6%。

圖4 Fs-φ關系曲線Fig.4 Curve of correlation between Fsand φ

表4 Fs-φ的敏感性Table 4 Sensitivity of Fs-φ

3.1.2 極限平衡法計算的Fs對φ值的敏感性分析

根據表3的穩定性系數可知，Fs-φ的關系如圖5，Fs-φ 近似于線性變化，取 φ =15°，c=9 ～23 kPa進行敏感性分析，其敏感性計算結果如表5。φ對滑坡穩定性敏感系數為19.3% ～46.0%，平均值為33.4%。

圖5 Fs-φ關系曲線Fig.5 Curve of correlation between Fsand φ

表5 Fs-φ的敏感性Table 5 Sensitivity of Fs-φ

3.2 c值對滑坡穩定性的敏感性分析

3.2.1 強度折減法計算的Fs對c值的敏感性分析

根據表3的穩定性系數可知，Fs-c的關系如圖6，Fs-c近似于線性變化，取 c=21 kPa，φ =10°～17°進行敏感性分析，其敏感性計算結果如表6。φ對滑坡穩定性敏感系數為21.6% ～36.1%，平均值為28.8%。

圖6 Fs-c關系曲線Fig.6 Curve of correlation between Fsand c

表6 Fs-c的敏感性Table 6 Sensitivity of Fs-c

3.2.2 極限平衡法計算的Fs對c值的敏感性分析

根據表4的穩定性系數可知，Fs-c的關系如圖7，Fs-c近似于線性變化，取 c=19 kPa，φ =10°～17°進行敏感性分析，其敏感性計算結果如表7。φ對滑坡穩定性敏感系數為12.1% ～43.8%，平均值為26.6%。

圖7 Fs-c關系曲線Fig.7 Curve of correlation between Fsand c

表7 Fs-c的敏感性Table 7 Sensitivity of Fs-c

3.3 強度參數的選取

根據以上分析結果表明，φ值相對于c值對滑坡的穩定性更敏感。參數值按照FLAC3D∶極限平衡法 =1∶1 選取，因此 c=20.0 kPa，φ =15°。

3.4 數值驗證

將以上分析得到的c，φ值進行數值模擬。在該值下，典型的應變增量云圖如圖8，在滑坡中下部其剪應變增量最大，在后緣剪應變增量相對較小，滑坡的主要破壞形式是剪切破壞，從運動形式分為牽引式滑坡。其水平位移等值線如圖9，其最大位移為0.45 m，與現場實際情況基本一致，此參數可用于滑坡支擋結構設計。

圖8 剪應變增量云圖Fig.8 Contour of shear-strain increment

位移等值線Fig.9 Displacement equivalence value maps

3.5 對比分析

將反算后的c，φ值用傳遞系數法計算，滑坡的剩余下滑力為235.32 kN，穩定系數為0.971;根據滑坡土體室內快剪試驗得出的巖土參數(c=19.2 kPa，φ =14.5°)計算出滑坡的剩余下滑力為 258.57 kN，穩定系數為0.962。二者值所計算的結果誤差在10%以內，能達到工程治理設計的精度。因此，反演計算的結果能用于工程實踐，指導工程設計。

4 結論

1)筆者基于自編強度折減法對滑坡抗剪強度參數進行反算，結果表明強度折減法計算所得的穩定系數與極限平衡法所得結果較一致，其相對誤差主要是由于極限平衡法假定滑塊之間為剛性體，不產生變形，數值模擬方法采用的是彈塑性材料進行巖土體的模擬，允許變形。通過強度折減法與極限平衡法的比較，按照一定比例選取c，φ值進行滑坡的治理工程設計可行。

2)強度折減法計算的穩定系數對c，φ值進行敏感性分析認為φ值相對于c值更敏感，與極限平衡法分析結果相一致。

3)在滑坡中下部剪應變增量最大，在后緣剪應變增量相對較小，滑坡的主要破壞形式是剪切破壞，從運動形式分為牽引式滑坡。

4)筆者僅僅是將一個特定的填方路基滑坡作為研究對象，通過改變c，φ值得出穩定系數，通過與極限平衡法的相互比較得出c，φ值，并通過室內試驗印證了該方法能滿足工程精度的要求。該方法對于其他滑坡同樣適用，可供大家參考借鑒，解決工程實際問題。

［1］王松齡.滑坡區巖土工程勘察與整治［M］.北京：中國鐵道出版社，2001.Wang Songling.Landslide Area of Geotechnical Engineering and Renovation［M］.Beijing：China Railway Publishing Press，2001.

［2］劉天韻，呂和林.滑坡土體抗剪強度指標反算研究［J］.西部探礦工程，2004(3)：185-187.Liu Tianyun，Lv Helin.Research on the back calculation of the landslide’s anti-sbear strength index［J］.West-China Exploration Engineering，2004(3)：185-187.

［3］潘家錚.建筑物的抗滑穩定和滑坡分析［M］.北京：中國水利水電出版社，1980.Pan Jiazheng.Analysis of the Building Sliding Stability and Landslides［M］.Beijing：China Water Power Press，1980.

［4］顏慶津.數值分析［M］.北京：北京航空航天大學出版社，2006.Yan Qingjin.Numerical Analysis［M］.Beijing：Beihang University Press，2006.

［5］陳育民，徐鼎平.FLAC/FLAC3D基礎與工程實例［M］.北京：中國水利水電出版社，2009.Chen Yumin，Xu Dingping.Basis and Engineering Application of FLAC/FLAC3D［M］.Beijing：China Water Power Press，2009.

［6］DZ/T 0219—2006滑坡防治工程設計與施工技術規范地質礦產行業標準［S］.北京：中國標準出版社，2006.DZ/T 0219—2006 Specification of Design and Construction for Landslide Stabilization［S］.Beijing：China Planning Press，2006.

［7］鄭穎人，趙尚毅.有限元強度折減法在土坡與巖坡中的應用［J］.巖石力學與工程學報，2004，23(19)：3381-3388.Zheng Yingren，Zhao Shangyi.Application of strength reduction fem in soil and rock slope［J］.Chinese Journal of Rock Mechanics and Engineering，2004，23(19)：3381-3388.

［8］廖秋林，曾錢幫，劉彤，等.基于 ANSYS平臺復雜地質體FLAC3D模型的自動生成［J］.巖石力學與工程學報，2005，24(6)：1010-1013.Liao Qiulin，Zeng Qianbang，Liu Tong，et al.Automatic model generation of complex geologic body with FLAC3D based on ANSYS platform［J］.Chinese Journal of Rock Mechanics and Engineering，2005，24(6)：1010-1013.