內部裂隙對節理巖質邊坡失穩影響的數值分析

馬鵬飛，郭德龍，許文年，陳 新*，夏 棟

(1.河海大學巖土力學與堤壩工程教育部重點實驗室，江蘇 南京 210024；2.中鐵十四局集團有限公司，山東 濟南 250014；3.三峽大學水泥基生態修復技術湖北省工程研究中心，湖北宜昌 443002；4.三峽大學水利與環境學院，湖北 宜昌 443002)

在形成現實巖體的過程中,自然界巖石綜合體經受風化作用、冷熱交替、水的破壞及改造,且在構造變動、工程建設、卸荷作用等條件狀態下,促使了巖石內部各種宏細觀裂紋裂隙的形成和發育[1]。這些損傷降低了巖層的強度并且不易被觀察,是地下硐室群開挖、地基工程、巖質邊坡等巖石工程的潛在威脅。地基工程和地下硐室工程基本上在有限變形理論領域里實施,其巖體或圍巖變形量是受限的；而人工開挖邊坡為一個有側向臨空面的半開放系統[2],在環境應力改變(如降雨、坡頂荷載)后既有瞬時變形又有與時間有關的流變屬性,邊坡一旦發生類似無限變形的顯著位移即可造成塌方、滑坡和堵江等嚴重的地質災害。因此,賦存于一定地質物理環境且作為工程對象的巖質邊坡的穩定安全性問題,不僅自身屬于自然科學問題,更是整個巖石力學界及工程地質學科最為關心的應用領域之一。

發育著不連續裂隙缺陷體的巖質邊坡,在內外界因素誘發下,內部裂隙源將沿一定的路徑發生拓展、匯合、互相銜接貫通直至最后宏觀裂縫(或破碎帶)的產生。工程實踐及大量反復的試驗已證實,巖體工程的失事大都與巖石節理的發育、原有結構面的衍生及其不穩定擴展演化過程相關聯[1,3-4]。也正是因為實際巖石中含有具備不同特性且以不同組合方式組成的結構面,巖質邊坡在失穩時,受裂隙的宏觀結構特征及其力學性質的影響嚴重,造成其破壞模式繁雜多樣,而節理面常被看作邊坡關鍵性破壞面的組成部分。所以，對含內部裂隙巖質邊坡滑移破壞機制與穩定性的進一步闡明具有重要而鮮明的指示意義。

由于巖體工程各種相互作用、材料性質及初始條件等問題的復雜性,以固體力學和經典數學為理論基礎的理論分析無法解決現有的許多工程難題[5],故而要想給出服務于工程需求的結果，采用數值方法在不少情況下已成為一種不可缺少或者不得已的手段。巖土數值計算通過對工程實體的計算機模擬,能從整體上獲得一些原本可能不了解的邊坡變形失效細節,甚至在計算中也可以了解如何去分析邊坡產生過度變形或破壞的原因,從而節約科研時間和經費。近些年新興的基于斷裂力學的動力擴展有限元法(XFEM)在處理物體內部不連續介質力學問題時特別有效[6-7],該方法不僅較為精準,并且還秉承了常規有限元法的優點,尤其在模擬剪切帶不連續場、裂隙開裂以及裂紋擴展方面具有明顯的優勢。

回顧先前的研究,主要存在著以下不足：①針對邊坡穩定性問題,巖土力學工作者采用強度折減有限元法基本對各種類型的邊坡的穩定都進行了廣泛的計算,但評價內部裂隙對邊坡滑動面的影響及其演化規律的研究仍偏少,目前的邊坡穩定性分析多是結合勻質邊坡模型展開的,亦或將邊坡模型建立為均質材料層的“堆砌體”,忽略了巖土體的裂隙性；②有些研究雖考慮了裂隙的影響,卻直接將帶裂隙服役的巖體通過設置很弱的強度參數來等效裂隙對邊坡整體穩定性的不利作用,導致數值仿真結果失真[8-9],即無法模擬出預期理想的邊坡介質材料開裂/破斷效果。實質上,上述做法視邊坡內部基巖或基土整體為連續體,未把裂隙等不連續界面夾雜其中,這也恰恰暴露了傳統連續方法如有限單元法(FEM)、有限差分法(FDM)等的局限性。為此,本文以Abaqus程序軟件環境下的XFEM二維模式為基礎,結合溫控參數的強度折減方法原理,通過幾個算例闡介了將XFEM運用于裂隙巖質邊坡相關問題研究的切實可行性，并簡要探討了不同的裂隙面排布狀況及含量對巖質邊坡的穩定系數、滑裂面演化路徑以及滑坡最終破壞形態的控制,揭示工程巖體失穩的漸進破壞現象,以期為指導水利水電、公路、采礦工程中邊坡的設計、穩定建設與科學修繕提供理論支撐。

1 基于擴展有限單元的強度折減法求解

1.1 動力擴展有限元法(XFEM)基本理論

有限單元法(FEM)在模擬微裂隙時,要預先確定裂隙擴展方向,并且劃分的單元邊界須遵從結構內部幾何或物理界面,特別是在裂隙端部,計算網格的布置要求苛刻[6],形成了裂尖附近區域極其密集其余區域稀疏的不均勻網格分布,從計算機的運行成本方面考慮是不合適的。XFEM由力學家Belytschko和Black[10]于1999年首次提出,該方法采用水平集法描述裂隙位置,是在標準有限元框架內的重要延伸。隨后,Mo?s等[11]進一步將階躍函數作為結點富集項引入,彌補了傳統有限元每個裂隙擴展步中都需要對裂尖區域重構網格以適應位移強、弱間斷等不連續邊界演化問題的缺點。

基于單元分解思想,引入反映裂隙局部特性的附加函數對有限元的逼近位移模式進行改進,將含裂隙物體的位移u(x)拆分為連續和非連續兩項之和,得到XFEM的近似位移場如下[7,12-13]：

(1)

在此位移模式基礎上,仿照常規有限元的虛功原理就可寫出XFEM的支配方程(有時又叫作求解方程)[14-15],然后求解總剛矩陣、整體等效節點荷載矩陣,繼而求解節點位移列向量。將XFEM應用于實際問題的分析中,裂隙的張開、閉合、生長和偏折行為完全與計算網格相互獨立。

1. 2 溫度相關參數強度折減理論及其方法實現

自Dawson等[16]提出強度折減法(SRM)以來,該方法就被普遍運用在工程穩定性系數(FoS，簡稱穩定系數)的自動搜尋中,其優點之一是不用事先知道破壞面。遺憾的是,Abaqus/CAE早年版本僅支持通過手動修改inp文件的方式來實現強度折減功能,不甚方便[17]。這里將等同于強度折減系數FV概念的溫度Temp置于邊坡體中(見圖1),通過相對簡潔的溫度相關參數強度折減計算得出FoS。應當說明的是：①溫度相關參數強度折減法是指由溫度變量關聯操縱整個強度折減過程的方法,其在使用時比較便捷,節省了在關鍵字中插入關于FV變化的描述語句這部分工作量,因此作者團隊將該方法命名為“convenient strength reduction method”，簡稱C-SRM；②Abaqus新版本是否內置了強度折減法,筆者不敢斷言；③下文各處出現的NT11就為圖1中的Temp。

圖1 抗剪斷強度折減定義窗口(示例)

按量變到質變理論,由下面公式(2)對邊坡巖體或土體所能提供的宏觀抗剪強度黏聚力c和內摩擦角φ實施連續的線性折減[17]：

ck=c/NT11,tanφk=tanφ/NT11

(2)

式中:ck、φk分別為維持邊坡平衡所需要或實際發揮的材料黏聚力(kPa)和內摩擦角(°);NT11表示坡體溫度(量綱為1)。

為了防止邊坡可能發生初始破壞,令NT11以固定步長從某一比1小的數值開始增加(溫度初始值一般習慣設定為0.25、0.5或0.7,折減步長的選擇則是依據對計算結果準確性的期望、對邊坡安全性的預估及一些書籍上的通常做法等),將不同的ck、φk代入計算模型進行計算,判斷邊坡系統是否處于平衡狀態，直至折減到邊坡到達破壞或臨界靜力失穩狀態為止,此時坡體的溫度NT11即為平常人們所說的邊坡安全系數(嚴格來講應該叫作邊坡穩定系數)[18]，且此時SRM自動識別出的最危險滑面所在的位置也即實際滑動面[19]。

C-SRM在ABAQUS有限元軟件中的實施流程要點如下：在Property模塊中定義與溫度相關的c、φ；進入Load模塊,先選擇分析步initial,對集合rock創建溫度分布場(宜事先對rock和crack-i部件各建立一個“Set”),溫度值填入溫度初始值,再選擇分析步reduce,將溫度進一步修正為某一足夠大的終值。

強度折減數值計算中邊坡失穩狀態的判定標準主要有：①邊坡域內最先形成連續的貫通區或塑性點呈大面積分布；②邊坡特殊關鍵部位的位移量出現突增現象；③有限元應力-應變方程在有限迭代次數內不收斂。

在進行各類裂隙發育邊坡的穩定性求解時,裂隙的存在誘致了所研究問題的極端非線性。數值計算結果不收斂情況有兩種：第一種是在邊坡達到極限平衡狀態之前單元就畸變,并引發模型過程過早停止,事實上,因模擬方法不夠完善這種現象較易出現,此時上述邊坡失穩的判定依據①～③均無用武之地,所以務必讓強度折減的整個過程持續到正常不收斂才可停止[措施包括：修改模型計算進程的默認分析參數和控制參數(分析步增量步長、允許的不連續增量數目I0、IR以及迭代不收斂次數IA等)；嘗試重置網格單元邊長]；另一種是在塑性應變區貫通后有限元分析不收斂,但邊坡內特征點位置處此刻已產生了無窮大的形變(或位移-折減倍數曲線末端出現幾近豎直的線段),說明邊坡已或者早已破壞,故據此得到的邊坡穩定系數偏大。因此，本文認為將邊坡失穩狀態的判定標準③作為裂隙巖質邊坡失穩破壞標準是不可取、不明智的。

1. 3 屈服準則

巖土介質的非線性變形規律比鋼材等其他材料要復雜得多,但邊坡巖土體的屈服和破壞符合Mohr應力空間中的經典M-C準則,同時M-C準則因其形式和使用簡單，在諸巖土失效準則中至今仍居于統領地位[20-22]，故文本在分析巖質邊坡巖塊體的強度問題時選用M-C條件,并認為材料屈服只取決于大、小主應力σ1、σ3,不考慮中間主應力σ2,其表達式為

(3)

理想彈塑性模型的庫侖-莫爾屈服條件可以表示為[20]

(4)

式中：I1和J2分別為應力張量的第一不變量和第二偏應力不變量，這里I1的單位即應力的單位，J2的單位則是應力單位的平方;θσ為應力羅德角(-30°≤θσ≤30°)；F為屈服面函數。

1. 4 數值模擬的基本假定條件

本次數值模擬的基本假定條件如下：

(1) 邊坡縱向視為足夠長,在確保計算效率的前提下進行平面應變狀態下的分析,對一些外輪廓形狀沿縱向變化緩慢的邊坡,多數情況下二維處理能近似滿足工程要求[23-24]。

(2) 對節理巖體而言,其包含著節理裂隙面,在模型結構上類似于復合材料結構,節理巖體裂隙界面的力學性質弱于基巖部分的力學性質。為了降低研究難度,將裂隙簡化為不占體積的空材料,并忽略其厚度,基質體按介質分布連續、均勻的理想彈塑性體考慮。

(3) 實際的工程巖體中存在規模尺度不一的大小裂隙,在現有技術條件下,很難把小微裂隙網絡也納入模型坡體,故暫只考慮內部長大裂隙,并將各種形態、特點的非連續結構面概化為坡頂裂縫和坡內斜直裂隙。

(4) 巖石材料的損傷起始判據應用最大主應力(Maxps)準則,其損傷演化規律則基于能量、線性軟化、剛度衰減方面的規則。Maxps損傷準則認為裂紋/裂隙沿著最大周向拉應力達到最大值的方向發展,待滿足界面裂紋起裂擴展準則時裂紋擴展。

(5) 滑體采用小位移計算模式,Abaqus-CAE前處理的Step模塊中幾何大變形開關仍撥置off,巖質邊坡某處斷裂的圖形表示效果可通過變形放大倍數調控。

(6) 斜坡面為單一直線。

2 算例分析

相比同規模土質邊坡,巖質邊坡滑坡能量更大,毀壞性更強,所以巖質邊坡穩定性評價成為廣大同行共同關注的熱點。本次數值模擬方案是建立多種人工設計節理巖質邊坡算例背景并對算例背景中的邊坡巖體在相互作用裂隙端部塑性區擴展影響下的穩定性響應(包括等效塑性應變、變形位移、折斷過程等)逐一進行初探,以此研究XFEM在含巖體結構面邊坡工程及在類似工程(特別如膨脹土邊坡[25]、崩崗治理區崩壁土坡[26-27]等)中的應用。

2. 1 對比算例分析

2.1.1 含4條非共線裂隙(背景一)的巖質邊坡算例

圖2簡示了某一假想的巖質邊坡地質斷面幾何模型,設該邊坡內部存在原始裂隙1～4(表面張裂隙1和埋藏型裂隙2、3、4)，數值分析模型中裂隙空間展布信息，見表1。

圖2 含4條非共線裂隙(背景一)巖質邊坡算例穩 定性計算的概化模型及裂隙布置(單位：m)

表1 含4條非共線裂隙巖質邊坡算例裂隙模型幾何參數

模型的單元庫類型使用CPE4R(平面四節點四邊形元,縮減積分),與完全積分單元相比,有意識地少取積分點是對高斯數值積分的一種處理,反而會使精度更高。計算網格大小以在模型分區邊線上布結點數的方式設置,由于研究側重于滑面附近及滑面左側巖塊的變形,故圍繞含裂隙重點計算域hbcg的邊緣上結點較稠密，以提高計算精度,而在擴展控制域外的gdeh區域分網較粗。模型有限元網格的結點數和單元數分別為1 395、1 317,但值得注意的是,不要對裂隙部件進行網格劃分。地基的底邊在X、Y方向上皆固支,在地基的左右兩側和巖質邊坡的右側安裝輥軸支撐,坡面自由約束。

在建模分析時不需要輸入任何關于內部結構面的力學參數,而輸入的巖體(結構塊)參數有：泊松比ν為0.29,彈性模量E為7.22 GPa,密度d為3.08 g/cm3,黏聚力c為100 kPa,內摩擦角φ為20°,塑性勢角ψ為0°(忽視體積變形),能量釋放率GC為500 J/m2,巖石所能提供的斷裂能KIC為151.1 N/m,巖石損傷穩定黏度系數為1×10-5，重力加速度g按10 N/kg計。為了使含裂隙巖質邊坡在自重一定的條件下產生變形破壞,借助材料降強的方式將邊坡帶入極限狀態[21],得到如圖3所示的該巖質邊坡失穩演化中塑性區的變化情況。

圖3 Abaqus塑性區云圖顯示的含4條非共線裂隙巖 質邊坡算例滑裂面延展過程

由圖3分析可知：①強度折減初期坡體的溫度NT11(亦稱降強系數)較低,也就是前期的邊坡巖石抗剪強度指標值遠高于實際黏聚力和內摩擦角,單元應力沒達到屈服面,因此就未形成塑性區；②塑性區最先在裂隙3的上端部位被隱約發現[見圖3(a)],緊接著是在裂隙3、4和2下端局部薄弱地帶出現[見圖3(b)]；③隨強度折減程度的加強,裂隙擴展前端塑性區繼續向前擴展,當強度折減到第15次時,裂隙尖端塑性區加速向對方延伸,即裂隙3、4之間塑性區明顯接觸,同時由于裂隙3靠近易遭破壞的坡腳,其下尖端與坡趾之間形成連通的滑弧通道，當強度折減到第17次時,裂隙4上端與裂隙2之間塑性區也互相連通；④伴隨坡體溫度NT11繼續遞增,裂隙2上部尖端塑性區逐漸向上發展,最終在后緣被張拉裂隙撕裂后,與裂縫1下部尖端相連接,坡內發育形成了從坡腳逐個連通裂縫,直至貫通坡頂的塑性區，此時的溫度值約為1.916、強度折減次數為21,按照邊坡工程失穩判據第①條,確定本算例邊坡的穩定系數FoS為1.916。

圖4 終極破壞時的含4條非共線裂隙巖質邊坡 算例狀態(兩級破壞形態)

由圖4可見,巖質邊坡內含有多條裂隙的條件下,初始破壞面(指最先形成的貫穿滑面)形態并不呈圓弧狀,而是相當不規則。分析其原因認為：巖體中的層面、節理、片理面、裂隙、裂縫、斷層等地質界面的力學強度一般比巖石材料弱得多,這些優勢結構面具控穩機制,體現在其對邊坡的滑坡型式、變形規模和速率、位移場展布以及附近巖石中應力傳遞都具有控制作用[28]。待巖質邊坡達到真正意義上的破壞狀態時,更大規模的次級滑動趨勢面已經醞釀完成。

邊坡特征點處的位移值和坡體整體運動具有相同的趨勢,為了結合位移突變理論得出邊坡符合失穩判據②時的穩定系數FoS,圖5中記錄了該巖質邊坡臨空面端點c水平位移(U1)隨降強系數(NT11)變化的發展趨勢。

圖5 非共線4條裂隙巖質邊坡算例坡頂監測位點c的 水平位移(U1)隨降強系數(NT11)的變化曲線

由圖5可見：該巖質邊坡臨空面端點c的水平位移(U1)在坡體溫度NT11由初始值增加到1.837的過程中只有微微的變化(當前階段特點是沒有滑坡跡象[29])；隨巖石材料塑性參數c、φ值的持續下降,臨空面端點c的水平位移U1產生了數量級的變化,位移切線角愈來愈大直到變形曲線上位移突變點顯現后,-U1出現大幅度激增,表示滑坡過程由漸變階段轉化到了紅色預警的急劇加速階段[29],實際上此時邊坡退化成了一個機構。

圖5中位移拐點難以辨識(邊坡穩定系數不易于合理地推求),但依舊可借助“雙切線法”來輔助確定：過U1-NT11關系曲線末端及曲線上降強系數NT11為1.837的點分別做2條切線,交于點R,再過點R做鉛直線,其與位移曲線的交點就是位移量拐點,對應的橫坐標在2.083附近。所以在邊坡失穩判據②下求得的邊坡穩定系數FoS≈2.083,與采用邊坡失穩判據①求得的FoS之間相差8.72%。考慮到邊坡失穩判據①可操性強、直觀性強、獲得的FoS值偏于保守,且模型中裂隙之間應變突變區動態連通過程即為邊坡漸進失穩的過程,故后文中統一依據塑性區分布情況來判別邊坡的失穩狀況。

2.1.2 含3條非共線裂隙(背景二)的巖質邊坡算例

為了論證內部裂隙數量對巖質邊坡穩定性與安全性的影響,對背景一巖質邊坡算例模型稍作改變，即剔除裂隙4后重新計算。圖6清晰地展示了6個典型強度折減分析時步時的模型等效塑性區發展云圖。

圖6 含3條非共線裂隙巖質邊坡算例的失穩過程

由圖6可以看出：含3條非共線裂隙巖質邊坡算例的塑性區首先在裂隙3的兩尖端處形成,然后裂隙2端部也出現了不明顯的塑性區；隨著強度折減過程的推進,邊坡坡腳處出現了局部剪壞面,并與裂隙3下端連通,且按弧線狀逐漸向上開展，但是由于受裂隙的影響迫使失穩滑面沿兩個方向擴展,一條是與裂隙2尖端相連接的主滑面(一級滑裂面),另一條為按原路線延伸的二級滑裂面[圖6(d)、6(e)]；在第20次強度折減完成時,邊坡域內塑性應變已貫通,但該巖質邊坡還未進入整體破壞,且在此后邊坡塑性應變繼續增大、二級滑裂面繼續拓展；當塑性區發展到一定程度(第450次強度折減),有限元計算不收斂,該巖質邊坡也徹底破壞,次級破壞面發展為與裂隙2交匯的第2條貫通塑性面[見圖6(f)和圖7]。

圖7 等比例放大1 717倍的含3條非共線裂隙巖質 邊坡算例網格變形圖

總體言之,若裂隙之間距離較遠,則兩裂隙中間部分可視作無裂隙擴展，換言之，隨著坡體溫度的升高(在現實中表現為作用于邊坡的自然力及其作用時間的增加,或人類應力擾動的日益加強),域內已形成的廣義塑性應變區向外傳播,直至塑性面延展到與邊坡內某個與之最靠近的預制裂隙(“nearest-fracture,簡記N-F”)之間距離足夠小時,塑性面之后的擴展才會顯然受到N-F的影響(對應于本文,當強度折減到圖6(c)的狀態時從坡腳引出的滑面才開始分叉)；相反,若內部裂隙之間距離較近,兩裂隙端部塑性應變則易發生連通(如裂縫1、2相隔較近,兩者之間塑性區先于裂隙2與3間的塑性區連通)。

2.1.3 含2條非共線裂隙(背景三)的巖質邊坡算例

設上述巖質邊坡內僅有裂隙1和2,其余計算條件均不變,重新計算。為了研究分析，圖8展示了施加加速度g引起模型在第13、17、20、48次強度折減時的積分等效塑性應變云圖。

由圖8可以看出：該巖質邊坡破壞區在裂隙兩側端面最早出現,這一點在前述諸背景算例中也得到體現，這是由于事實上裂紋/裂隙是應力釋放和應力產生奇異的關鍵部位,其局部高應力必會致使材料屈服,直接影響著巖體的變形和穩定，此外,眾所周知脆性巖石抗壓不抗拉,且邊坡滑動時巖石的抗拉強度未能充分發揮,故塑性區首先出現在圍繞裂隙尖端一個小范圍拉應力集中區域；而后,坡趾附近受力較大也漸漸被征服,當強度折減到第20次時求得了最弱的滑動面[見圖8(c)],它是剛剛貫通的臨界破壞區,此時邊坡轉動破壞機制為坡面破壞；主滑移面貫通后如強度折減過程繼續下去,第二條塑性面也將貫通,直至邊坡發生最終的坡腳破壞[見圖8(d)]。

圖9給出的是該巖質邊坡滑動位移矢量圖。

圖9 含2條非共線裂隙巖質邊坡算例的最終破壞 形態(變形放大595倍)

由圖9可大概判斷出該巖質邊坡穩定失效時的滑體范圍，但因裂隙2端部距離坡腳的位置太遠,或者說裂隙2邊緣與坡腳之間的部分相當于無裂隙區,坡內預制裂隙幾乎影響不到二次貫通塑性區的發展路徑,所形成的過坡趾的整體性潛在失穩滑面(與之對應的降強系數為2.331)與裂隙1、2都不相交。

2.1.4 不含裂隙(背景四)的巖質邊坡算例

在背景一巖質邊坡算例的基礎上,刪除計算模型中的所有裂隙。因均質巖質邊坡不含裂隙,滑動面剪入口自然就不會承受后緣裂隙的影響,為了降低邊界效應,根據圣維南原理(僅僅對周圍一定區域有影響,對更遠的地方影響可忽略不計),現將圖2中的cd邊、of邊往右加長了437.5 cm。不難看出,邊界范圍已給圖8(d)示出的計算結果帶來了一定影響。

采用傳統有限元法對此模型的穩定性進行分析,利用最后2個強度折減步(Frame)得到該均質巖質邊坡增量位移分布云圖[見圖10(a)],可判定其潛在滑帶大致呈弧形,符合完整邊坡滑動面的形態特點,這與考慮裂隙發育巖質邊坡的破壞機制存在較大的差異。溫控強度折減過程中該均質巖質邊坡最早明顯貫通的破壞面，見圖10(b)。

圖10 均質巖質邊坡內部大致的滑動面形狀及位置

2.1.5 含2條共線非貫通裂隙(背景五)的巖質邊坡算例

圖11 含2條共線非貫通裂隙的巖質邊坡示意圖

設定的巖質邊坡巖石物理力學性質同第2.1.1節,圖12大致展示了模型的運行結果。

圖12 利用XFEM模擬的含2條共線非貫通裂隙(背景 五)巖質邊坡連續塑性破壞面和裂縫的演化

由圖12(b)可以看到,計算域內出現一條通過非連續弱面且貫穿的潛在滑面,在裂隙1′起裂以后邊坡將啟動平面滑動型失穩模式。

結合軟件中輸出參數PHILSM(裂紋面水平集函數)，可查看計算不收斂時裂隙擴展的結果,結合響應量STATUSXFEM(裂隙擴展狀態參量)可得兩裂隙之間的損傷值(1代表該單元完全開裂；0代表該單元無裂隙),見圖12(c)和12(d)。即,在邊坡系統到達最后失穩破壞狀態時,裂隙5、6在各自兩尖端處沿著與巖橋呈某一鈍角方向有所擴展(長度約6～12 cm),但萌生出的次生或翼裂隙較為短小。

綜上,表2列出了不同背景算例下計算得到的巖質邊坡最小穩定系數FoS。

表2 不同背景算例下計算得到的巖質邊坡最小穩定系數FoS對比

由表2可知：一方面，只要巖質邊坡內存在裂隙,含裂隙巖質邊坡的FoS就明顯低于均質巖坡的FoS；另一方面,當裂隙數目一樣時,即使裂隙6比裂隙3長了1 m,但背景二下巖質邊坡的FoS比背景五下巖質邊坡的FoS還小,這可能的原因是裂隙尖端位置和兩條平行埋藏裂隙的傾角不同。

假設有這樣一種背景算例(命名為**)：將背景二中裂隙2和裂隙3繞各自的下部尖端逆時針旋至與bc面平行,則旋轉后的裂隙2、3也接近共線,其他計算條件均不改變。可計算出背景**下巖質邊坡的FoS為2.290(>2.086)。可見,節理邊坡穩定系數對滑動層面(結構面)傾角的變化比較靈敏,因此在滑坡治理工程中,針對內部含不同傾角、產狀裂隙的巖質邊坡應給予區別對待。

該巖質邊坡(指圖11模型)的最終破壞形式及相應的總位移分布特征，見圖13。

圖13 含2條共線非貫通裂隙(背景五)巖質邊坡的 最終破壞形式及相應的總位移分布特征

由圖13可見，沿著被該巖質邊坡裂隙切割出的復合滑動面上方單元巖塊位移明顯大于其下方滑床位移,滑裂面的上部和中部處均有脫離母巖現象,該模擬研究結果與預期變形情況(即巖質邊坡已滑出失穩并向坡底滑落)的符合度較高。

2.2 獨立算例(滑移路徑呈階梯狀的巖質邊坡)分析

在中國西部山區,由于巖體結構經歷了多期表生演化作用,階梯式貫通變形破壞常見于受淺表部陡、緩兩組結構面控制的卸荷巖質高邊坡的破壞模式中,其破裂面形態整體上為陡-緩相間的階梯狀[30]。為了進一步考察XFEM的可用性與推廣性,本文以文獻[31]中的平行斷續節理巖質邊坡模型為實例(與第2.1節無關),就更為復雜難解的裂隙性邊坡變形穩定問題嘗試開展額外的仿真和分析。該巖質邊坡剖面地形結構如圖14所示,內含有一組與邊坡傾向一致的平行等間距主控節理裂隙,巖石主要物理力學參數的取值見表3。

圖14 平行斷續節理巖質邊坡模型(90°巖橋 )[31]

表3 巖石主要物理力學特性

對圖14模型進行擴展有限元計算,得到該巖質邊坡在程序正常不收斂時的裂隙擴展情況、運動矢量及相關等值線圖(涉及塑性應變、剪應力S12、水平位移U1、最大主應力),分別列于圖15。經數次反復的溫控強度折減試算可知：當塑性區剛剛從坡腳到坡頂貫通時對應的巖質邊坡模型溫度(即FoS)約為5.255；如果清除該巖質邊坡內部的節理裂隙,由FEM求出該巖質邊坡的FoS約為8.17。

圖15 獨立算例的動力擴展有限元模擬結果匯總

由圖15可得出幾點有參考價值的信息：①觀察圖15(a),在原始裂隙尖端處,裂隙發生了擴展,同時伴有翼裂隙,巖橋失效,斷續結構面間在宏觀層面形成連通；②坡體整體失穩后的塑性應變帶由4個相互連通的近橢圓狀的塑性區域組成,見圖15(b),并且橢圓長軸的兩端點基本位于原始主控裂隙面的中心,橢圓短軸的兩端與原始裂隙尖端重合；③播放圖15(d)的動畫發現,坡體先沿臨近坡腳的節理產生局部滑移,而后再沿著由斷續節理系以及擴展萌生的裂隙所組成的滑動面發生整體朝下的階梯狀貫通滑裂；④以貫通面為界線,邊坡位移場分布在模型內產生明顯分異(與圖13相一致)，且對于圖15(e),X軸方向位移值在可滑動塊體的腳部位置處達到最大,坡體中部位移較大,坡肩部位位移次之,這說明階梯狀變形破壞機制具有“牽引漸進式”破壞性質；⑤如圖15(f)、(c)所示,節理裂隙的存在將引起巖石應力重分布，裂隙尖端部位應力較為集中,說明開裂擴展還可能會繼續下去，而相對于邊坡其他地帶,節理周圍處剪應力大小也較大,這與周子涵等[32]利用FLAC3D得到的數值計算結果很類似。

3 動力擴展有限元法(XFEM)的優劣分析

以上建立的各背景算例巖質邊坡數值模型大體上達到了預期效果,有限元計算所得的一系列成果圖與現有文獻中的研究結果可良好吻合(如邊坡主滑動區集中出現在巖體結構面上部坡體區域)。基于XFEM的溫度相關參數強度折減法(C-SRM)意義明確,不僅可以一次性地讓Abaqus算出斷續裂隙巖質邊坡的滑裂面位置及最小穩定系數值,而且有限元計算中無需裂隙力學參數的幫助就可正確地反映巖質邊坡的變形及滑動面孕育過程(若在以往,要通過假定、室內測試、經驗公式等方式事先獲得節理裂隙的強度參數),同時本文提出的方法能針對折線滑面以及同時包含折線和曲線的滑面進行分析,并能較好地描述裂隙對邊坡穩定狀態的削弱作用,間接驗證了基于XFEM的C-SRM具有更加高效、合理與簡潔等優越性。

但是，本文方法也存在一些局限性：無法預測內部含完全貫通節理(包括節理從坡面延伸到坡頂和從坡趾延伸到坡頂兩種情況)巖質邊坡的穩定系數；無法對邊坡坡內的彎折裂隙及交叉型裂隙進行模擬；此外,XFEM計算的收斂性尚需提高。

4 結 論

不是使用過去習慣的傳統有限元強度折減法,而是在提出基于溫度場變量的便捷降強計算方法的基礎上,本文改為采用XFEM(從而使內部含非連續弱面這一類型邊坡的穩定問題能夠被有效地處理和評價)研討和開拓了滑坡全過程中多裂隙巖質邊坡內部裂尖塑性區的萌生、發展規律,以及幾種簡單的節理巖質邊坡模型在靜力荷載下可能的失穩模式。通過研究獲得了如下結論：

(1) Abaqus中的XFEM提供了直觀的裂隙形式,有別于條分法等經典的極限平衡法,采用XFEM得到的巖質邊坡主滑面更符合滑坡形成特點,由于日后可考慮某具體邊坡工程后緣張裂隙的存在,故不會高估巖質或土質邊坡的穩定系數；而且將XFEM與強度折減法結合還能尋找出兩個潛在的失穩滑面。因此,本文所做工作為XFEM在裂隙邊坡穩定性計算領域的應用提供了有益參考，建議將XFEM在相近工程實踐中加以推廣應用。例如,可考慮利用XFEM模擬南方7省區崩崗區花崗巖殘積土裂隙的動態行為,這無疑對崩壁崩塌機制的深層次和多角度發掘具有重要的意義。

(2) 裂隙的存在縮短塑性區貫通的路徑長度(就好比裂隙擴展必須要沿著指定的薄弱面),促使邊坡的抗力下降,因此一般情況下,裂隙數目越多,邊坡的穩定系數就越小。

(3) 完整巖石坡體在外部荷載驅動下的失穩破壞面是穿過坡頂的典型圓弧形或對數螺線形滑裂面；含裂隙節理巖質邊坡的破壞始于裂隙尖端形成的塑性區,其失穩破壞面是通過斷續裂隙的不規則滑移面,且滑面形態與裂隙的發育程度及排布息息相關；階梯狀滑移方式形成的連通塑性區是由若干離心率較小(豐滿)的橢圓塑性帶與原有節理面共同組合而成的。

(4) 受內部裂隙的影響,巖體邊坡處于極限狀態時的位移場特征形成了以貫通滑面為界線的明顯分區；在順著裂隙長度方向上,邊坡巖石應力分布規律也呈現出很強的統一性。

限于篇幅,下一步的研究工作擬開展含連續巖體節理(完全貫通)、內部裂隙傾角變化范圍更廣(如優勢結構面逆向坡、反傾層狀邊坡)、坡頂存在多條豎直裂隙等工況作用下巖質邊坡穩定性及其滑移規律的數值模擬研究,并與已有解析解(含有解析解的直線滑面的邊坡算例)做對照分析。

致謝:謹向參加本論文審稿并提出建設性修改建議的專家、老師致以最真誠的謝意！