基于黏聚力裂縫模型的反傾層狀巖質邊坡傾倒破壞模擬

馬文著，徐 衍，李曉雷，陳 敏

(1.河北工業(yè)大學土木與交通學院，天津 300401；2.北京科技大學土木與資源工程學院，北京 100083; 3.中水北方勘測設計研究有限責任公司，天津 300202；4.宿遷市宿城區(qū)水利局，江蘇 宿遷 223800)

反傾層狀巖質邊坡的傾倒破壞是常見的一種地質災難，尤其是在坡腳開挖或河谷下切作用下，很容易發(fā)生傾倒失穩(wěn)。如Wali滑坡等[1]，正是由于傾倒變形誘發(fā)的大規(guī)模滑坡。很多礦山工程和水利水電工程就面臨這樣的工程問題，如三峽庫區(qū)龔家方滑坡[2]、湖北鶴峰紅蓮池鐵礦山體崩塌[3]等。因此，研究此類邊坡的傾倒破壞機制和影響因素對滑坡災害的控制和治理具有非常重要的意義。

前人對巖質邊坡的破壞問題進行了大量數值計算與模擬研究。在研究反傾層狀巖質邊坡的破壞機制時，能夠考慮大變形的離散元程序(如UDEC)已經被不少學者用來進行層狀巖質邊坡的分析。劉云鵬等[4]采用UDEC軟件，以都汶公路一處邊坡為例，研究了其在地震響應下的變形破壞機制；黃波林等[5]對三峽庫區(qū)廖家坪危巖體進行了離散元模擬分析；王章瓊等[6]借助離散元方法以紅蓮礦山為背景，分析了多因素下坡體變形破壞特征；Gu等[7]對龔家方滑坡進行離散元模擬，概化出“懸臂梁模型”，用以解釋坡體的破壞機理。Alzoubi等[8]采用UDEC損傷模型(UDEC-DM)對反傾層狀巖質邊坡進行離心破壞模擬，結果表明：巖層的抗拉強度是導致坡體發(fā)生傾倒破壞的關鍵因素。

在研究反傾層狀巖質邊坡傾倒破壞在連續(xù)小變形條件下的成因機理方面則主要采用有限元模型來進行計算。王宵等[9]以某梯級電站廠后邊坡為例，構建三維有限元模型，深入闡明了“似層狀”巖質邊坡的具體演化過程；代仲海等[10]通過對反傾層狀邊坡的有限元分析，揭示了其破壞面發(fā)展的時空規(guī)律；王宇等[11]運用節(jié)理有限元法，構建出反傾層狀巖質邊坡有限元模型，探討了此類坡體的變形破壞機制及穩(wěn)定性影響因素。

基于F-DEM方法(有限元-離散元耦合)的連續(xù)-離散數值模型既能研究巖石彈性階段的小變形機理又能夠分析大變形情況下巖石的破壞模式。針對邊坡變形破壞問題，陶志剛等[12]以戒臺寺古滑體為背景建立連續(xù)-離散模型對其開裂破壞過程進行研究，判斷出古滑體最危險區(qū)域。陳小婷等[13]以巫山縣長江左岸危巖體為例，證明了將連續(xù)-離散模型用于邊坡破壞研究的可行性。有限元黏聚力裂縫模型(Cohesive Crack Model, CCM)是連續(xù)-離散模型中的一種，該模型最早由Dugdale[14]和Barenblatt[15]提出，并在晶體、巖石等脆性材料裂紋擴展研究中得到廣泛應用。相較于其他連續(xù)-離散模型，CCM更注重模擬物質的微觀破裂形態(tài)。Jiang等[16]運用此方法進行了刀齒破巖的數值模擬研究，得出不同形狀壓頭在侵入巖石時侵入力與侵入距離的關系。Jaime等[17]采用此方法進行了巖石切割的模擬研究，結果表明：此方法可以合理地模擬巖石切割過程。Wang等[18]將此方法與擴展有限元法進行耦合，分析了水力作用下巖石的裂隙擴展過程，得出巖石裂縫幾何形狀對裂紋擴展應力的影響。CCM在模擬巖石破壞行為方面應用廣泛，但針對邊坡破壞問題的算例極少。

本文在通用有限元軟件ABAQUS平臺上，通過全局嵌入零厚度黏聚力單元，建立反傾層狀巖質邊坡CCM。經過參數標定和對比得到與鄭達等[19]離心試驗模型相似的CCM。基于以上模型及參數，研究邊坡傾倒破壞的演化過程，分析邊坡的層間應力分布特征。最后，探究層間剪切強度對邊坡巖層初始折斷位置、垮塌范圍、破裂面傾角的影響，簡要探討巖層層厚對滑坡最終形態(tài)的影響。

1 邊坡數值模型

1.1 黏聚力裂縫模型(CCM)

1.1.1CCM基本原理

如圖1所示，黏聚力單元將材料界面分離破壞的過程用實體單元兩個面之間的“相對分離位移-力”關系表達，這種關系被稱為牽引-分離準則[20-21]。該準則定義了黏聚力單元的力學行為，適用于內部存在黏結作用的材料(巖石、混凝土等)，將0厚度的黏聚力單元插入實體單元后構成CCM，此模型通過黏聚力單元的刪除，較準確地模擬巖石的斷裂行為。該準則假定沿著破裂面的法線方向存在一個牽引力σn，沿著破裂面的切線方向存在兩個牽引力σs、σt，這些牽引力隨破裂面相對位移的增加，先增大后減小。當這些牽引力減小至零時，單元刪除。

圖1 黏聚力單元示意圖Fig.1 Schematic diagram of the cohesive element

本文采用如圖2所示的混合型牽引-分離準則定義節(jié)理面以及巖層內的零厚度黏聚力單元。其中tn0、ts0、tt0、tm0分別為黏聚力單元達到初始損傷時的法向應力、剪切應力以及混合模式的應力，δn0、δs0、δt0、δm0分別為黏聚力單元達到初始損傷應力時所對應破裂面法線方向、剪切方向及混合模式方向的相對位移，δnf、δsf、δtf、δmf分別為黏聚力單元失效刪除時所對應的破裂面法線方向和剪切方向的相對位移以及混合模式方向的相對位移。定義損傷因子D(0≤D≤1)來表示黏聚力單元損傷過程的發(fā)展，在黏聚力單元達到初始損傷強度后，損傷因子D的值從0開始逐漸增加，當損傷因子D的值達到1時黏聚力單元被刪除，巖塊發(fā)生斷裂。

圖2 黏聚力單元破壞模式Fig.2 Failure mode of the cohesive element

顯然，這個過程也是單元剛度衰減過程。黏聚力單元三種破壞模式下的應力損傷關系式為：

(1)

式中：tn——黏聚力單元法向應力；

ts、tt——黏聚力單元切向應力。

黏聚力單元損傷演化過程中的法向及剪切剛度可表示為：

kn=(1-D)kn0

ks=(1-D)ks0

kt=(1-D)kt0

(2)

式中：kn——黏聚力單元法向剛度；

ks、kt——黏聚力單元切向剛度；

kn0——黏聚力單元法向初始損傷剛度；

ks0、kt0——黏聚力單元切向初始損傷剛度。

對于混合模式下黏聚力單元損傷演化，使用等效總位移來描述，可表示為:

(3)

式中：〈〉——Macaulay括號(當δn>0時〈δn〉=δn,否則〈δn〉=0)；

δn——破裂面法向相對位移；

δs、δt——破裂面切向相對位移。

損傷因子D在損傷演化階段應滿足以下關系：

(4)

式中：δmm——加載過程中允許的最大位移。

黏聚力單元在損傷演化過程中應當遵循最大主應力原則，該原則可以表示為：

(5)

式中：〈〉——Macaulay括號。

如圖2所示，根據黏聚力單元面法向牽引力tn和切向牽引力ts、tt與其對應方向的初始損傷強度tn0、ts0、tt0之間的關系以及其組合公式，大致分為以下三種破壞模式：

①若沿著單元面的法向牽引力tn先達到初始損傷強度tn0，那么單元最終的破壞模式則為張拉破壞；

②若沿著單元面剪切方向切向牽引力ts、tt先達到與其對應的初始損傷強度ts0、tt0，那么單元最終的破壞模式則為剪切破壞；

③如果符合式(6)，那么最終的破壞模式則為拉剪破壞。

(6)

層狀巖體中的節(jié)理面和巖層內的巖塊之間都考慮這三種破壞模式。

1.1.2CCM構建過程

由于通過ABAQUS軟件的界面操作在每一個實體單元間插入0厚度黏聚力單元是困難的，因而需要借助程序對有限元模型進行前處理。首先建立原始幾何模型并劃分網格，導出模型文件。之后通過MATLAB程序對模型文件中單元編號以及節(jié)點坐標進行更改，大致思路如下：

①將原始模型中相同的節(jié)點分別賦予兩個不同的序號，目的是將單元進行離散，為零厚度單元的插入做準備；

②將產生的新編號的節(jié)點數據以及開裂區(qū)單元節(jié)點的變化添加到模型文件中；

③添加黏聚力單元的數據信息，即在源文件的關鍵字中聲明單元類型的更改；

④為這些新加入的單元設置材料和截面屬性并再次導入ABAQUS中進行邊界條件等的設定。

1.2 反傾層狀巖質邊坡CCM建立

1.2.1數值模型信息及加載

將3 586個零厚度黏聚力單元嵌入到2 455個實體單元之間，并根據離心機模型(圖3)按1∶1比例構建相對于邊坡的CCM(圖4)。模型尺寸為800 mm×670 mm，巖層傾角70°，層厚15 mm；模型分階段開挖三次，開挖坡比均為1∶0.18，開挖后坡腳均為80°，模型Z方向僅劃分一層單元并約束模型所有節(jié)點在厚度方向上的自由度。鄭達等[19]在試驗過程中并未使用機械手在離心機轉動時進行開挖，而是將離心機停轉后通過人工手段實現開挖工況。為保證離心機試驗模型和相應的數值計算模型在開挖條件上一致，數值模型加載過程如下：

通過“單元殺死”模擬開挖條件，將原始邊坡穩(wěn)步加載至120g，穩(wěn)定轉動1 200 s后，加速度降至0g，進行第一級開挖；再穩(wěn)步加載至120g，如此往復直至第三級開挖。

圖3 離心機模型示意Fig.3 Centrifuge model

1.2.2黏聚力裂縫模型中的時間

使用ABAQUS中Dynamic、Explicit分析步進行計算，此分析步為基于動力學方程的動態(tài)顯式算法。此算法中時間即為運動方程中的時間，是真實時間。模型設置的最小時間步長為0.01 s，并通過定義幅值曲線的方式對加速度時程曲線進行輸入。

圖4 黏聚力裂縫模型Fig.4 Cohesive crack model

1.3 黏聚力裂縫模型參數及標定

如圖4所示，模型中的黏聚力單元按嵌入位置不同，分為在巖層之間的層間黏聚力單元與巖層里的基質體(非層間黏聚力單元)。由于CCM模型屬細觀唯象模型，其參數不具備嚴格的宏觀物理意義，因而需要與試驗進行標定和對比，以確定黏聚力單元的具體參數。如果經標定和對比后的CCM模擬結果與試驗結果基本吻合，那么就可以認為模型建立兩者相似[16]，也就驗證了CCM的正確性。具體標定過程如下：

(1)巖石材料及基質體參數標定

通過模擬抗壓強度試驗可以對基質體參數進行標定[16]。建立直徑80 mm、高120 mm的圓柱形試件模擬抗壓強度試驗，與鄭達等[19]給出的巖層試驗材料的彈性模量和單軸抗壓強度進行標定，得出基質體內的黏聚力單元參數(表1)。

表1 數值模擬中各單元的力學參數

(2)巖層面參數標定

通過將數值模型監(jiān)測點位移曲線與鄭達等[20]離心機試驗結果中的位移曲線進行對比(圖5)，對CCM中的層間黏聚力單元參數進行校核，最終得到巖層面(巖層間黏聚力單元)參數(表1)。

圖5 位移的比較Fig.5 Displacement with time

2 數值模擬結果

2.1 數值計算與模型試驗對比

圖6 傾倒變形特征的比較Fig.6 Comparison between the numerical toppling deformation and experimental results

圖6表明在離心機試驗模型及相應的數值計算模型中，各階段的坡體變形特征基本一致，驗證了CCM的正確性。原始模型階段，邊坡并未產生明顯變形，大多數的變形產生于自重下邊坡頂部的沉降，且邊坡后緣范圍內層間發(fā)生較小的相對位移。一級開挖階段，坡體后緣范圍內的層間相對位移繼續(xù)增加，巖體層間裂隙進一步增多，表現為坡體距表面約150 mm的斷續(xù)裂隙，坡表巖體有向臨空面“點頭哈腰”傾倒的趨勢，但未見傾倒。二級開挖階段，延續(xù)上一階段的變形，邊坡表層巖體繼續(xù)向臨空面變形，層間相對位移較之前增大，層間裂隙發(fā)育更為明顯。三級開挖階段，較以上兩階段有明顯變化，邊坡巖層出現傾倒現象，坡體后緣拉開，邊坡自坡腳至坡頂產生貫通破裂面，形成一級破裂面，且坡腳處出現局部垮塌。最終破壞階段，坡表巖體沿一級彎折帶發(fā)生失穩(wěn)破壞形成滑坡，并同時出現二級破裂面，這說明了此方法在模擬層狀巖體傾倒變形方面的合理性。

圖7 演化過程圖Fig.7 Picture of the evolution process

2.2 邊坡破壞演化過程

由于原始模型階段以及一級開挖、二級開挖階段，坡體變形較小，因此本節(jié)只針對模型破壞演化主要發(fā)生的三級開挖和最終破壞階段進行演化過程的分析，結果見圖7。隨著重力加速度的提高，坡體變形加劇，離心加速度達到42g左右時(圖7a)，層間已發(fā)生明顯錯動，坡腳處出現一條較緩的近乎水平的破裂面，且出現初始折斷，此時坡腳巖體并未壓碎擠出。離心加速度增至60g左右時，距離坡體后緣約40 mm 處出現明顯的張裂縫，坡腳處巖體被少量壓碎但并未被大范圍擠出，層間錯動加劇(圖7b)。離心機加速度100g左右，坡腳處巖體碎裂區(qū)域進一步擴大，層間錯動進一步加劇，且局部巖體被擠出，破裂面進一步向坡頂方向擴展將坡體貫通，形成一級破裂面(圖7c)。值得注意的是，隨著層間錯動的增加，深層坡腳破碎區(qū)萌生出距表層較遠的斷續(xù)裂隙，這說明二級破裂面正在發(fā)育。離心加速度達到120g左右時(圖7d)，坡腳處的破碎巖體被大范圍擠出，一級破裂面貫通形成滑坡，同時二級破裂面貫通。這說明，反傾邊坡最終的破壞形式是沿著某條貫通的破裂面發(fā)生整體失穩(wěn)破壞，且層間錯動對破裂面形成有促進作用。

綜上得知，反傾層狀邊坡的演化是一個動態(tài)過程：從坡腳處開始，坡腳處巖體斷裂，隨后坡體后緣出現張裂縫，使層間錯動加劇，致使坡腳處巖體被碾碎擠出；而坡腳處巖體的擠出又反過來加劇了后緣的張開以及層間的錯動，兩者相互作用，最終導致滑坡災害發(fā)生。綜合原始模型階段與一級二級開挖階段模型的變形過程(圖6)可知，坡體臨空條件的改變是滑坡災害發(fā)生的重要因素，開挖實質上是為坡體巖層發(fā)生層間錯動提供了變形空間。坡體的整個破壞演化主要分為兩個時期：第一時期主要表現為臨空條件改變導致的層間錯動；第二個時期為邊坡的傾倒變形破壞，主要表現為破裂面的發(fā)育。

2.3 邊坡層間應力分布特征

以往的反傾層狀邊坡研究經常簡化為疊合懸臂梁模型[22-23],因而巖體層間作用力特別是法向作用力的作用形式對于反傾層狀邊坡的破壞模式具有重要影響。

圖8 監(jiān)測路徑示意Fig.8 Monitoring path

圖9 層間法向應力圖Fig.9 Normal stresses between layers

設定9個局部坐標系(圖8)，基于每一個局部坐標系，設定數據提取路徑，從坡腳至坡頂共9條，定義為1 #路徑、2 #路徑、……、9 #路徑。離心機加載至100g時,破裂面接近完全貫通，提取此時各個監(jiān)測路徑上的層間法向應力繪于圖9。如圖所示，坡腳位置以1#路徑為例，從坡體表面到坡底0～150 mm范圍內，出現兩個主要的法向應力峰值，且法向應力波動明顯，隨著距離的增加法向應力并未出現明顯增加。坡體中部位置以4 #路徑為例，從坡體表面至坡底 450 mm 范圍內，出現兩個主要的法向應力峰值，距離坡體表面約200 mm范圍內應力波動較為明顯,且沿著監(jiān)測路徑方向的深度的增加，應力未出現明顯增長。坡體后緣以8#路徑為例，從坡體表面至坡底550 mm范圍內，僅出現一個主要的法向應力峰值，應力波動較小，且隨著距離的增加，法向應力表現出不斷增大的趨勢。從路徑應力峰值大小的角度進行分析，層間法向應力整體上呈現由坡頂至坡腳逐漸增大的趨勢，這說明距離坡腳越近巖體破壞越嚴重。

此外隨著坡高的增加，層間法向應力有向坡體內部增大的趨勢，尤其是在8 #和9 #監(jiān)測路徑表現明顯，產生此現象的主要原因是隨著距離坡體表面深度的增加，巖體自重應力產生的法向分量增大。對各個監(jiān)測路徑中層間法向應力峰值出現的位置進行分析可知，坡腳處監(jiān)測路徑(1—4 #路徑)有2個較為明顯的法向應力峰值，而坡體中部與坡體后緣位置的監(jiān)測路徑(5—9 #路徑)只有1個較為明顯的法向應力峰值，且峰值的位置恰好對應1—4 #路徑范圍內的兩條破裂面的位置。繼續(xù)沿著坡頂的方向，將5—9# 監(jiān)測路徑法向應力峰值所在位置在坡體對應位置上標出，恰好就是破裂面的所在位置。也就是說，坡體的破裂面基本會沿巖層間法向應力峰值所對應位置的連線發(fā)育。

3 層間剪切強度對邊坡破壞特征的影響

3.1 層間剪切強度對初始折斷位置的影響

初始折斷位置影響滑坡形態(tài)。為探究層間剪切強度對初始折斷位置的影響，按不同層間剪切強度將模型分為a(0.25 MPa)、b(0.15 MPa)、c(0.05 MPa)三組，其他條件不變。將巖層初始折斷時的位移云圖繪制于圖10。層間強度取0.25 MPa(a組)，巖層在t=650 s左右時發(fā)生初始折斷，折斷位置接近坡腳(圖10a)。層間強度取0.15 MPa(b組)，巖層在t=550 s前后發(fā)生初始折斷，折斷位置接近坡體中部(圖10b)。層間強度取0.05 MPa(c組)，巖層大約在t=450 s時發(fā)生初始折斷，折斷位置接近坡體后緣(圖10c)。

由此可知，層間剪切強度不同，坡體初始折斷時間和初始折斷位置不同。隨著層間剪切強度的增大，巖層初始折斷位置自坡腳至坡頂逐漸升高。巖層初始折斷時間隨層間剪切強度的增大而延后。巖層層間抗剪強度越大，坡體發(fā)生初始折斷時的整體位移越小。

圖10 不同層間剪切強度下位移云圖及初始折斷區(qū)位置Fig.10 Displacement and location of initial fracture zone under different shear strength between layers

3.2 層間剪切強度對滑坡形態(tài)及破裂面的影響

將不同層間剪切強度條件下模型的計算結果繪于圖11。層間剪切強度為0.25 MPa時(a組)，巖層于坡腳處垮塌且邊坡垮塌區(qū)域的范圍較小(圖11a)。層間剪切強度為0.15 MPa時(b組)，坡體下半段垮塌區(qū)域的范圍與a組相比明顯增大。層間剪切強度 0.05 MPa 時(c組)，坡體垮塌區(qū)域的范圍進一步增大，發(fā)生整體垮塌(圖11c)。分析得知，隨著巖層層間剪切強度的降低，垮塌區(qū)域的范圍由坡腳向坡頂逐漸增大，直至發(fā)生整體垮塌。

將圖10所示初始折斷區(qū)的位置在圖11中對應標出，可以發(fā)現初始折斷區(qū)的位置在垮塌區(qū)域上半段。這說明初始折斷位置對滑坡形態(tài)影響明顯，并且有一定預測作用。

在圖11中，α和β分別為一級破裂面和二級破裂面與模型X正方向的夾角。由圖11(a)可見層間剪切強度為0.25 MPa條件下破裂面的最終發(fā)育形態(tài)：一級破裂面和二級破裂面與X正方向夾角分別為116°和127°。圖11(b)是層間剪切強度為0.15 MPa條件下的破裂面最終形態(tài)，α和β的值分別為130°和135°。圖11(c)是層間剪切強度為0.05 MPa 條件下的破裂面傾角和位置，一級破裂面和二級破裂面與X正方向夾角分別為137°和146°。隨著坡體巖層層間剪切強度的降低，破裂面與模型X正方向的夾角逐漸增大。與一級破裂面相比，距坡表更遠位置的二級破裂面角度變化更為明顯。

圖11 不同層間剪切強度下滑坡形態(tài)Fig.11 Landslides under different shear strength between layers

4 討論

巖層厚度對此類邊坡的穩(wěn)定性也具有較大的影響[11,23]。為探討巖層厚度對反傾邊坡滑坡形態(tài)的影響，分別建立含不同巖層厚度(10,15,20 mm)的反傾邊坡模型并進行對比分析，巖層傾角和坡腳分別固定為70°和80°。圖12給出了不同巖層厚度的邊坡模型計算結果。結果表明：隨著巖層厚度的增大，坡體由三級傾倒(巖層厚度為10 mm)逐漸整合成一級傾倒(巖層厚度為20 mm)，坡體的破碎程度逐漸降低，且一級傾倒破裂面逐漸向坡體深部發(fā)展。由此可見：巖層厚度越大，坡體的整體性增強；一級傾倒破裂面的深度越大，坡體傾倒部分提供的彎矩越大，說明了反傾邊坡具有更大的抗彎承載力。

圖12 不同坡體層厚下滑坡形態(tài)Fig.12 Landslides under different layers thickness

5 結論

(1)經過參數標定得到的CCM，能夠與試驗結果基本吻合，證明了CCM可用于層狀巖質邊坡的離心模擬分析，并能作為處理這類層狀巖體邊坡的有效計算方法。

(2)通過坡體演化過程的分析可知，反傾層狀巖質邊坡的演化是一個動態(tài)過程，巖層內部破壞是誘發(fā)滑坡的直接因素，層間相互錯動加速巖層內部的破壞，而臨空條件的改變則促進層間相互錯動。

(3)破裂面位置可以通過監(jiān)測層間法向應力獲得，破裂面大致沿層間法向應力峰值位置的連線發(fā)育。

(4)層間剪切強度對坡體破壞形態(tài)影響顯著，層間剪切強度越大，巖層初始折斷區(qū)域越接近坡腳，巖層垮塌區(qū)域范圍越小，破裂面傾角越緩。

(5)坡體層厚越大，坡體的整體性越強，抗彎承載力越大。