基于位移突變準則的顆粒流強度折減法

戴健健， 朱振學， 葉 楠， 蘇 超*， 張 恒

(1.河海大學水利水電學院， 南京 210098； 2.吉林省水利水電勘測設計審查總站， 長春 130033； 3.吉林水利科學研究院， 長春 130022)

邊坡穩定分析作為巖土工程中重要的研究課題，其評價合理與否直接關系到邊坡工程的成敗。邊坡失穩破壞時往往存在巖土體的滑動、平移與轉動等現象，這是一個復雜的過程且在宏觀上表現出非連續性。顆粒離散元[1]作為一種典型的非連續介質分析方法，在邊坡工程中得到廣泛應用，李龍起等[2]采用顆粒流程序PFC(particle flow code)模擬降雨作用下土質邊坡的變形特征；戎澤鵬等[3]聯合PFC3D和Rockfall對危巖體進行分析；張志飛等[4]使用PFC研究反傾層狀巖質邊坡的變形及破壞過程；Wu等[5]利用PFC對永嘉山體滑坡的動力過程進行分析。在邊坡穩定分析中，強度折減法[6]由于概念明確、容易實現，周健等[7]將其引入到顆粒流中，采用累計位移或不平衡力作為臨界破壞準則，為顆粒流用于評價邊坡的穩定性提供了新思路。采用顆粒流強度折減法評價邊坡穩定性時，能否直觀合理地判斷邊坡處于臨界破壞狀態顯得至關重要。

近年來關于顆粒流強度折減法在邊坡穩定性方面的研究較多。王培濤等[8]選取了黃金分割點作為折減比例，采用不平衡力收斂準則對鐵礦現場巖質邊坡開挖擾動前后的變形特性和穩定性進行了研究。王鑫等[9]以坡頂位置的顆粒作為特征點，通過繪制不同折減系數與特征點x方向位移的關系曲線，以曲線出現明顯拐點時的折減系數作為穩定系數，對沿海吹填土堤邊坡的穩定性進行了顆粒流模擬。田雷等[10]在邊坡頂部設置測量圓，繪制測量圓內顆粒的平均位移與折減系數關系曲線，以曲線位移突變時的折減系數作為穩定系數，分析了植被護坡對邊坡穩定性的影響。汪華安等[11]通過繪制坡面測點的最大位移與折減系數關系曲線，采用曲線轉折點處的折減系數作為邊坡穩定系數。金磊等[12]提出一種基于能量演化的方法來判別邊坡是否失穩破壞，當顆粒做功和摩擦耗能超過顆粒接觸應變能，且保持一定速率持續增長時則判定為失穩狀態，并將該方法用于土石混合邊坡穩定性分析中。

上述關于顆粒流強度折減法失穩破壞判據的研究多集中在顆粒不平衡力、位移和顆粒體系能量方面。顆粒不平衡力比值以及位移臨界值的選取均具有主觀性與經驗性，目前尚無統一的標準。采用顆粒位移與折減系數關系曲線求解穩定系數直觀明了，但關于特征顆粒位置以及位移的選擇有待進一步探究。對位移突變準則進行改進，通過繪制坡面不同位置的特征顆粒位移與折減系數關系曲線求解邊坡穩定系數，并與有限元強度折減法和極限平衡分析的結果進行對比，為進一步驗證該方法的合理性與可靠性，模擬了滑坡的破壞形態，并從變形和應力特征兩個不同角度分析其破壞機理。

1 邊坡算例及土體細觀參數標定

1.1 邊坡算例

圖1為某土質邊坡工程的斷面示意圖，土體力學參數如表1所示。

表1 土體力學參數Table 1 Soil mechanical parameters

圖1 邊坡算例斷面圖Fig.1 Cross-sectional view of slope calculation example

1.2 土體細觀參數標定

利用PFC2D程序對土體進行雙軸壓縮試驗標定其細觀參數，構建的土體雙軸壓縮試驗的模型尺寸為3.0 m×1.5 m(高×寬)如圖2所示。土體顆粒采用接觸黏結模型[13]，該模型能真實地反映黏土類黏性材料的宏觀力學特性。采用試錯法經大量計算最終確定土體的細觀參數如表2所示。選取100、200、300 kPa 3個不同的圍壓，對土體試樣進行雙軸壓縮試驗，可得試樣的峰值強度分別為206.32、379.12、543.48 kPa。通過MATLAB回歸分析可求得土體的凝聚力為15.08 kPa，內摩擦角為14.79°，土體在不同圍壓下的偏應力-軸向應變曲線如圖3所示，其莫爾圓與強度包絡線如圖4所示。

圖2 雙軸壓縮試驗模型Fig.2 Biaxial compression test model

圖3 偏應力-軸向應變曲線Fig.3 Deviator stress-axial strain curve

τ為剪應力；σ為正應力；c為凝聚力；φ為內摩擦角圖4 莫爾應力圓與強度包絡線Fig.4 Mohr’s stress circle and strength envelope

表2 土體顆粒的細觀參數Table 2 Micro-parameters of soil particles

2 邊坡顆粒流模型構建

根據該邊坡工程的斷面圖，采用PFC2D與AutoCAD相結合的方法建立邊坡模型，運用geometry import與wall import geometry命令將dxf格式文件導入生成邊坡的輪廓與墻體，運用ball distribute命令生成土體顆粒，在建立顆粒流計算模型的過程中，以顆粒的平均比率小于1.0×10-5時作為平衡條件，在模型達到平衡狀態后進行顆粒速度和位移的清零，賦予土體顆粒的細觀參數。最終建立的邊坡顆粒流模型如圖5所示，計算區域取坡頂向右延伸10 m，坡高12 m，顆粒總數為120 888個，兩側和底部采用墻體約束。此外，在邊坡內部4個不同位置布設半徑為1.0 m的測量圓，以記錄坡體在滑動過程中應力的變化情況。

a～d為測量圓圖5 邊坡顆粒流模型Fig.5 Slope particle flow model

3 顆粒流強度折減法求解邊坡穩定系數

顆粒流強度折減法[6]的基本思路為將顆粒的法向、切向黏結強度和摩擦系數同時進行折減，其計算公式為

(1)

式(1)中：Fs為滑坡的穩定系數；cbtenf為顆粒的法向黏結力；cbshef為顆粒的切向黏結力；fric為顆粒間的摩擦系數；cbtenfcr、cbshefcr和friccr分別為顆粒的臨界法向黏結力、切向黏結力和摩擦系數。

對于同一邊坡而言，不同研究者因選取臨界破壞判斷準則的不同將得到完全不同的計算結果，這將導致將該方法用于邊坡穩定性的評價具有強烈的主觀性與經驗性。為更加直觀、合理地求解邊坡穩定系數，同時為減少人為的主觀性與經驗性，本文通過繪制特征顆粒位移與折減系數關系曲線求解邊坡的穩定系數，具體思路為：①對土體顆粒的法向、切向黏結強度和摩擦系數同時進行折減；②選取一系列位于坡面位置的特征顆粒，記錄計算完成時特征顆粒位移值；③繪制選取的特征顆粒位移與折減系數關系曲線；④對于每個特征顆粒，其位移有明顯突變趨勢時有一個對應的折減系數，取所有對應的折減系數的最小值作為邊坡的穩定系數。

3.1 特征顆粒的選取

采用位移與折減系數關系曲線求解邊坡穩定系數直觀明了，在顆粒流中，位移可采用顆粒的位移或測量圓內顆粒的平均位移。由于測量圓的位置是固定的，若在坡頂布設，坡頂會出現明顯張拉縫隙，坡體下滑導致測量圓內的顆粒數量在不斷變化。測量圓位置布設在滑動面區域最為合適，但滑動面事先是未知的。在邊坡發生破壞時，坡面的顆粒隨坡體一起滑動，有十分明顯的位移變化，因此本文選用坡面位置的顆粒作為特征顆粒。

在滑坡過程中，考慮到坡腳處的顆粒會出現滑移堆積以及坡頂處滑坡體易出現二次破壞(圖6)，因此選取坡腳至坡頂中間處6個坡面位置的顆粒作為特征顆粒，選取原則為選用邊坡輪廓線處的顆粒以及顆粒的位置分布均勻，特征顆粒的位置如圖7所示。

藍線為邊坡的輪廓圖6 坡腳和坡頂處特征顆粒Fig.6 Characteristic particles at the toe and top of the slope

圖7 特征顆粒位置示意圖Fig.7 Schematic diagram of the location of characteristic particles

3.2 特征顆粒位移與折減系數關系曲線

經大量的滑坡顆粒流計算并綜合考慮計算時間成本，最終確定顆粒流強度折減法的計算時步取為5×105步。對邊坡施加重力，記錄不同折減系數計算完成時特征顆粒的位移值，繪制的坡面位置特征顆粒位移與折減系數關系曲線如圖8所示。

圖8 特征顆粒位移與折減系數關系曲線Fig.8 The relationship curve between characteristic particle displacement and reduction coefficient

3.3 邊坡穩定系數的求解

對于每個選定的特征顆粒，其位移有明顯突變趨勢時有一個對應的折減系數。該折減系數的確定方法為：①當曲線有明顯拐點時，取拐點對應的折減系數；②當特征顆粒位移與折減系數關系曲線有明顯加速趨勢段時，通過繪制關系曲線在位移拐點附近的兩條切線確定。

根據圖8所示的特征顆粒位移與折減系數關系曲線，確定的折減系數如表3所示，從安全的角度出發，確定邊坡的穩定系數為1.46。

表3 特征顆粒對應的折減系數Table 3 Reduction factor corresponding to characteristic particles

為進一步探究上述方法求解邊坡穩定系數的適用性，首先基于強度折減法采用有限元軟件ABAQUS求解邊坡的穩定系數。邊坡有限元計算模型共有12 700個單元，13 005個結點，采用平面應變單元CPE4，如圖1所示。采用位移突變準則求得的穩定系數為1.586，計算不收斂時的塑性區如圖9所示。其次對該邊坡進行極限平衡分析，采用Bishop法求得的穩定系數為1.612。圖10為相應的滑動面位置。采用顆粒流強度折減法求得的穩定系數比有限元強度折減法和Bishop法的更小，這是由于顆粒流作為一種非連續數值方法，通過力-位移關系模擬顆粒間的接觸，當法向力或切向力超過允許的強度后，顆粒間的接觸黏結可發生斷裂導致的。

圖9 計算不收斂時的塑性區Fig.9 Plastic zone when calculation doesn’t converge

圖10 Bishop法滑動面Fig.10 Bishop method sliding surface

采用改進的位移突變準則，通過繪制特征顆粒位移與折減系數關系曲線求得的邊坡穩定系數為1.46。對于巖質邊坡只需采用可描述巖石力學特征的接觸模型，亦可采用上述方法求解其穩定系數；而對于三維邊坡，可在不同剖面位置選取坡面不同位置的顆粒作為特征顆粒，采用同樣的方法進行求解穩定系數。

4 滑坡破壞過程模擬及其機理分析

4.1 變形特征分析

根據上述方法求得的穩定系數對土體顆粒的法向、切向黏結強度和摩擦系數同時進行折減1.46倍，對邊坡施加重力進行破壞形態的模擬，計算時步取1×106步。

由顆粒流數值模擬的結果，在計算初期，坡腳處的土體顆粒變形不斷累加，在重力作用下，在5×104步時坡頂向下沉陷，坡腳處顆粒向外擠出，部分顆粒的接觸黏結斷裂。1×105步時坡腳處發生擠壓破壞變形加劇，坡內顆粒接觸黏結斷裂形態向圓弧形式發展，土體發生剪切破壞。在上覆壓力作用下，斷裂的接觸黏結貫穿坡頂，貫通面形成，坡頂失去支撐的作用，2×105步時坡頂產生明顯張拉縫隙，坡體發生失穩破壞，沿貫通的圓弧式斷裂黏結滑移[圖11(b)]。隨著時間的推移，斷裂的接觸黏結繼續發展，破壞的坡體繼續下滑，3×105步時滑坡體顆粒間部分黏結破壞，出現二次破壞，由于土體具有一定凝聚力，滑坡體在一些局部位置發生張拉破壞，滑坡體局部破碎，并滑移堆積至坡腳處。5×105步時坡體滑動面與Bishop法求得的滑動面以及有限元強度折減法計算不收斂時的塑性區形態一致[圖11(d)]，這表明采用本文方法求解邊坡穩定系數的合理性與可靠性。圖11(e)為7×105步時滑坡的形態，在1×106步時滑坡體趨于穩定。由模擬的滑坡過程可知，顆粒流強度折減法很好地再現了整個滑坡破壞過程，而有限元強度折減法只能給出塑性區的分布，極限平衡分析如Bishop法只能給出圓弧滑動面的位置。

紫色線條為原狀邊坡的輪廓線圖11 邊坡滑坡破壞形態Fig.11 Slope landslide failure mode

4.2 應力特征分析

采用上述求得的穩定系數對土體顆粒的細觀參數進行折減，對邊坡施加重力進行1×106步的計算。該邊坡滑坡破壞過程中記錄的應力(以拉為正)變化曲線如圖12所示，坡頂位置土體顆粒由于接觸黏結斷裂產生張拉裂縫和滑移，測量圓a內x方向和y方向應力由起初增大的拉應力變成壓應力，當坡體下滑時壓應力逐漸平緩變化。隨著黏結斷裂的增加，坡體持續向下滑移，測量圓b內x方向和y方向壓應力變化平緩。測量圓c內y方向壓應力變化比x方向大，這是由于在重力作用下，斷裂的接觸黏結貫穿坡頂，隨著坡體的滑移，滑坡體顆粒間的部分黏結亦破壞，滑坡體發生張拉破壞和擠推作用所導致的。測量圓d位于坡腳位置，由于坡體下滑時的擠推和堆積作用，使得x方向壓應力比y方向更大。

圖12 邊坡測量圓的應力變化曲線Fig.12 The stress change curve of the slope measurement circle

5 結論

(1)考慮到在滑坡過程中坡腳處顆粒會出現滑移堆積，坡頂處滑坡體易出現二次破壞，建議選取坡腳至坡頂中間處坡面位置的顆粒作為特征顆粒。

(2)采用本文方法對邊坡穩定性進行評價時，可模擬顆粒間接觸黏結斷裂導致土體發生剪切破壞，斷裂的黏結逐漸貫穿坡頂，坡體發生失穩破壞，很好地再現了整個滑坡破壞形態，這是有限元強度折減法和極限平衡分析無法模擬的。

(3)本文方法模擬得到的邊坡破壞形態與有限元強度折減法和Bishop法的結果一致，驗證了該方法的合理性與可靠性。