堰塞體形成全過程的連續離散耦合數值模擬

王 葉，周 偉，馬 剛，陳 遠，常曉林

(1. 武漢大學水資源與水電工程科學國家重點實驗室，武漢 430072；2. 武漢大學水工巖石力學教育部重點實驗室，武漢 430072)

自然界中堰塞體的形成因素很多，根據統計，有90%以上是地震誘發滑坡形成的[1]，滑坡體堆積河谷形成堰塞體，一旦失穩將引發下游洪水災害，因此深入研究地震作用下的堰塞體形成機制，對于評估堰塞體的穩定性、堰塞湖潛在次生災害可能性及其損失具有重要意義。

目前對于堰塞體形成機制的研究方法主要有有限元法(FEM)[2,3]、離散元法(DEM)[4-7]、非連續變形分析方法(DDA)[8,9]，連續-離散耦合分析方法(FDEM)[11-17]等。徐文杰等[2]在對肖家橋滑坡區詳細的地質調查基礎上，結合動力有限元分析進行，再現滑坡三維空間失穩過程；劉傳正等[3]通過有限差分方法研究了紅石巖山體在地震的動力響應特征和崩塌的形成機理。在強震作用下，巖土體發生局部大變形甚至形成強間斷面(巖石拉斷、沿結構面滑移)，此時基于連續介質力學的傳統有限元方法和有限差分法受小變形假設的限制而很難適用。曹琰波等[4]采用離散單元法，對唐家山滑坡由變形累積到破壞滑動的全過程進行了模擬，研究了地震作用下順層巖質滑坡全過程；Wu等[5]采用離散元方法分析一個地震誘發的大規?；?，研究滑坡體滑動過程和滑坡后的形態等；焦玉勇等[6]采用離散單元法模擬了邊坡失穩過程；鄔愛清等[8]采用非連續變形分析方法，以唐家山滑坡形成的堰塞壩形態和位置作為目標函數，對唐家山滑坡過程進行復演，探討地震荷載下高速滑坡的形成機制；Zhang等[9]提出附加接觸模型的改進的DDA方法，研究有黏聚力材料的邊坡模型，評估邊坡節理面處的抗剪強度；Huang等[10,11]采用考慮節理動摩擦退化的改善非連續變形分析方法(IDDA)研究地震誘發的滑坡體的速度、位移等動力學行為。

DEM和DDA能模擬邊坡失穩及滑動全過程，但是它們對裂紋、應力、應變的表征不夠直觀，因此有學者開始采用連續離散耦合分析方法來研究邊坡失穩及滑動全過程。Barla等[12]采用FDEM對位于意大利Alpetto礦山的2處邊坡進行數值模擬，與早前的研究進行比較，驗證了FDEM模擬邊坡失穩的可行性；Antolini等[13]采用FDEM方法對Torgiovannetto di Assisi滑坡體進行數值模擬，研究滑坡體的啟動機制和演化過程，為滑坡體的風險評估提供參考。Piovano等[14]研究了意大利Aosta山谷的Beauregard滑坡，表明運用FDEM方法能夠預測深層滑坡的典型不穩定面和破壞機理；常曉林等[15]采用FDEM模擬了在暴雨等情況下邊坡失穩滑移的過程，較好地模擬了抗剪斷參數降低時邊坡從小變形-大變形-局部滑動-整體滑移的全過程失效；Zhou等[16]基于FDEM方法研究了節理傾角對邊坡破壞模式、運動過程、破壞后的堆積形態等的影響。O. K. Mahabadi等[17]應用基于FDEM的Y-Geo代碼較好地模擬了山崖侵蝕的過程。但是FDEM進行滑坡全過程模擬大多限于二維，對三維滑坡體研究甚少。

在連續-離散耦合分析方法中引入界面單元，采用內聚力模型能有效地捕捉加載過程中裂紋的擴展路徑，模擬巖體漸進破壞過程，充分發揮連續介質力學方法和非連續介質力學方法各自的優勢[15-18]，實現巖體從連續狀態到非連續狀態的轉化過程。本文采用FDEM，結合紅石巖堰塞體，對三維滑坡由變形破壞到整體滑動的全過程進行模擬，以研究地震作用下巖質滑坡的失穩破壞過程，直觀地捕捉滑坡體從開裂破壞到最終的堆積狀態、滑坡體失穩過程中的速度分布情況等。

1 基于內聚力模型的巖石開裂擴展模擬

在巖體連續-離散耦合分析方法中[12-18]，巖石被離散為由實體單元和無厚度界面單元組成的系統，每個實體單元均為單獨的離散塊體。界面單元失效之前，離散塊體之間通過界面單元“黏結”成一個整體，以模擬巖石材料的連續變形。實體單元只發生彈性變形，損傷和斷裂僅發生在界面單元上。界面單元的破壞準則為帶拉斷的Mohr-Coulomb準則，界面單元的應力狀態滿足破壞準則后，采用基于斷裂能的線性損傷演化模型，損傷因子達到1.0后完全失效。界面單元失效后從模型中刪除，兩側實體單元發生接觸關系，采用罰函數法進行接觸分析。

采用內聚力模型模擬巖石材料的界面開裂行為，通過界面單元的失效模擬裂紋的萌生、擴展、交匯和貫通。在張拉應力狀態下，簡化的內聚力區應力分布見圖1。真實裂縫的尖端應力為零，內聚力區尖端應力為材料的抗拉強度，從真實裂縫尖端到內聚力區尖端應力逐漸提高。二維界面單元具有2個積分點，均位于單元的中平面上。采用共節點的方式實現界面單元與相鄰實體單元之間力、位移的傳遞。圖2給出局部模型的有限元網格拓撲結構圖，裂縫可從布置界面單位的位置擴展。

圖1 準脆性材料的內聚力區應力分布Fig.1 Illustration of cohesive zone for quasi brittle material

圖2 網格拓撲圖(實體單元面積收縮便于顯示界面單元，實際界面單元厚度為0)Fig.2 Mesh topology (the elastic elements are shrunk for illustration purpose, the actually element thickness is zero)

在界面單元出現損傷之前，假定其本構關系是線彈性：

(1)

式中：tn、ts分別為界面單元的法向應力、切向應力；Kij為剛度矩陣，本文計算時未考慮應力、應變的法向分量與切向分量之間的耦合關系，因此Kns=0 ；εn、εs分別為界面單元的法向應變、切向應變。

法向應變、切向應變定義為：

(2)

式中：δn、δt分別為界面單元沿著法向和切向產生的位移；T0為界面單元的本構厚度。

本文采用無厚度的界面單元模擬材料開裂，倘若取T0為幾何厚度，將會導致應變奇異。為了消除應變奇異，在進行界面單元的應變和應力計算時采用的是單元的本構厚度。通常定義T0=1，則界面單元的應變與相應方向上的位移大小相等。界面單元的線彈性本構方程可寫成：

(3)

式中：kn、ks分別為界面單元的法向剛度、切向剛度。

為了模擬裂紋的萌生、擴展過程，需要定義合理的裂縫起裂準則，考慮巖石等準脆性材料的破壞往往是由拉應力和剪應力共同作用導致的，本文采用二次應力準則判斷界面單元是否開裂，采用帶拉斷的Mohr-Coulomb準則計算巖石材料的抗剪強度：

(4)

式中：c為巖石材料的黏聚力；φ為巖石材料的內摩擦角；t0n為抗拉強度；t0s為抗剪強度。

二次應力準則：當各個方向的應力與其相應的臨界應力的比值的平方和達到1時，界面單元即開始出現損傷：

(5)

在損傷演化階段，界面單元的本構方程為：

(6)

ts=(1-D)ksδs

式中：D為無量綱的損傷因子，當D=0時，表明界面單元未出現損傷；當D=1時，表明界面單元完全失效，失去承載能力。

損傷因子的表達式如下：

(7)

式中：t0eff為達到破壞準則時的等效應力；δ0m為達到破壞準則時的等效位移；δmaxm為加載歷史中的最大等效位移。

在大多數情況下，界面單元處于在混合加載狀態下，同時發生法向和切向的變形，因此需要定義界面單元的等效應力teff和等效位移δm：

(8)

通過定義界面單元在損傷過程中耗散的能量，即可控制界面單元的損傷演化過程。在本文中，耗散能量也就是通常所說的斷裂能，它等于應力-分離曲線下面所包含的面積大小。本文模擬時采用基于線性軟化的Benzeggagh-Kenane準則，界面單元的應力-分離曲線見圖3。

圖3 界面單元的應力-分離曲線(Ⅰ、Ⅱ)Fig.3 Constitutive relations of the cohesive elements

Benzeggagh-Kenane準則的表達式如下：

GT=Gs+GⅠ

(9)

式中：Gshear為界面單元在剪切荷載作用下的斷裂能；GT為界面單元在混合荷載作用下的斷裂能；η為材料常數，通常由彎曲試驗測得，本文取2。

2 紅石巖滑坡體FDEM三維模擬

2.1 工程概況

紅石巖堰塞壩形成于2014年發生的云南魯甸“8·03”地震。“8·03”地震發生后，左岸滑坡堆積物表層松動并向河床滑動，右岸山體產生大規?；拢诤哟残纬裳呷８鶕F場調查，堰塞體頂部呈馬鞍形，頂部左岸高，右岸低，右岸邊緣為滑坡巖石堆積體。堰塞體組成松散，最大厚度約103 m。頂部橫河向最低高程點1 222 m，堰塞體左岸最高點為1 270 m，上游迎水面平均坡比約1∶2.5，下游面平均坡比約1∶5.5。順河向底寬約910 m，沿1 222 m高程壩軸線長度約307 m。右岸滑坡后，地形發生了較大變化，上部滑床后緣為陡崖，高度估計約150～200 m，中部形成一個橫河向近水平且傾向下游的斜面地形，在此斜面以下為陡崖。左岸坡雖有表面松散并向下滾落，但體量不大，總的地形未有大的改變。堰塞體物質組成較為復雜，主要為右岸邊坡崩滑后形成的崩塌堆積物。

2.2 FDEM模擬所需參數率定

采用連續-離散耦合分析方法進行數值模擬時，需要對所用參數進行率定。采用二維平面應變模型進行雙軸數值壓縮試驗來率定模擬參數，通過不同圍壓下的偏應力應變曲線得到應力莫爾圓從而計算出巖土體的宏觀抗剪強度參數。用于參數率定的數值試樣尺寸為80 mm×160 mm(見圖4)。剛性加載板與試樣之間定義接觸關系，采用位移控制方式加載施加于上、下剛性加載板，圍壓施加在左右加載板上。剛性加載板與數值試樣之間的摩擦系數取為0.1。

圖4 雙軸壓縮試驗試樣尺寸及加載示意圖Fig.4 Size of numerical experiments of biaxial compression and illustration of loading

紅石巖滑坡體各土層的物理力學參數見表1。在FDEM數值模擬中，界面單元的參數包括剛度、抗拉強度、黏聚力c、內摩擦角φ、斷裂能。通過反復試算，使得數值試驗得到的力學參數與表1所示的室內試驗值接近。表2對比了數值試驗得到的各土層宏觀黏聚力、內摩擦角與室內試驗各土層抗剪強度，相對誤差都在10%以內。最終確定的FDEM模擬參數見表3。

表1 紅石巖巖體室內試驗物理力學參數Tab.1 The physical and mechanical parameters of laboratory tests for Hongshiyan rocks

表2 紅石巖巖體力學參數與數值試樣計算力學參數結果對比Tab.2 The comparison of mechanical parameters ofrocks for Hongshiyan rocks between laboratorytests and numerical experiments

表3 紅石巖巖體數值模擬界面單元細觀參數Tab.3 The mesoscopic parameters of cohesive interface elements for numerical simulation of Hongshiyan rocks

2.3 計算模型

根據震前邊坡原地形圖進行滑坡前的還原。由于地震影響范圍較廣，為減小邊界反射效應，將模型邊界范圍取得足夠大。三維計算模型中，x軸方向為橫河向，由左岸指向右岸為正；y軸方向為順河向，向下游為正；z軸方向為豎直向，向上為正。橫河向長1 427.1 m，縱河向1 620.5 m，豎直向970 m。

地震荷載作用下，只考慮邊坡底部的彈性作用，消除邊坡底部對地震的放大作用，采用有剛度無質量方案進行分析，故底部約束為固定約束，側向邊界為法向約束。在盡量考慮實際地形地貌和地質結構的前提下，同時考慮計算效率，三維整體模型共剖分為341 651 個單元和376 002 個節點，其中實體單元為215 271個，界面單元為126 380 個。本文重點關注滑坡形成堰塞體的整個過程，所以只在滑坡體部分插入界面單元。整個FDEM模擬分2步進行，首先進行地應力平衡，得到初始應力場，再施加地震加速度荷載，進行地震分析。

在計算模型底部輸入3個方向的地震加速度，模擬滑坡體在地震作用下的啟動、滑動到最終堆積直至形成堰塞體的全過程。加速度時程曲線見圖5(a)、5(b)、5(c)，地震作用過程歷時20 s(10～30 s)，總的作用時間為100 s。

圖5 地震加速度動力時程曲線Fig.5 Time-history curves of earthquake acceleration

3 堰塞體形成全過程分析

圖6為堰塞體最大剖面對比圖(對比范圍相同)，實際堰塞體最大堆積厚度約為103 m，模擬的對應縱剖面最大堆積厚度為101.97 m，由于模型的簡化和模擬范圍限制，2者的堆積形態稍有不同，但大致相似。

圖6 堰塞體最大縱剖面對比Fig.6 Comparison for longitudinal section of landside dam

定義堰塞體堆積最大寬度為Lmax，堆積最大高度為Hmax，寬高比Lmax/Hmax，用這3個參數表征堆積形態。表4統計了4個典型剖面堆積形態實測值和FDEM計算值及2者的相對誤差(百分制)。由于模型簡化和網格劃分關系，表4中個別計算值與實測值有較大差別，故其相對誤差較大。圖7對比了4個典型橫剖面實測與FDEM模擬的堆積形態。從表4和圖7可以看出模擬得到的橫剖面堆積形態與實測堆積形態基本一致，進一步驗證了本文FDEM模擬邊坡失穩、滑動形成堰塞體的合理性。

表4 典型橫剖面堆積形態實測值與計算值對比Tab.4 The comparison of accumulation state of typical sections between measured value and calculated value

滑坡體不同時刻的滑動狀態見圖8。選取4個時刻滑坡典型剖面g-g′的形態[見圖9(a)～圖9(d)] 研究地震誘發滑坡失穩過程運動特征。從圖8和圖9可以看出，在地震初期由于受地震作用，表層巖體首先出現裂縫然后貫通形成整體滑裂面；隨著地震作用，滑坡體內部損傷不斷累積，t=20 s左右時，整個滑坡體與下覆基巖出現明顯的運動分離[見圖9(a)]；然后在地震和重力共同作用下滑坡體高速下滑，且在下滑過程中滑坡體不斷破碎、解體，受對岸山體的阻擋，與之發生高速撞擊，進一步破碎、解體、堆積、堵塞河道，在60 s后滑坡體逐漸穩定，逐漸堆積密實，100 s時基本靜止，形成堰塞體。

圖7 典型橫剖面實測與計算堆積形態對比Fig.7 The comparison of accumulation state of typical sections between actual measurement and calculation

圖8 滑坡體不同時刻滑動形態Fig.8 Sliding patterns of landslide at different time

圖9 滑坡體典型剖面g-g′不同時刻滑坡形態Fig.9 Sliding patterns for a typical section named “g-g′” at different time

圖10為滑坡體典型剖面g-g′不同時刻速度云圖。t=20 s左右滑坡體速度開始增大，滑塊底部首先達到8 m/s左右，在地震和重力作用下，滑坡體的速度持續增大，t=30 s左右高速下滑，達到20 m/s左右，速度的分布受邊坡坡度的影響，由于上部坡度較陡而中部坡度緩，下部與中部坡度相近，但下部是臨空面，而中部有下部巖體阻擋，故巖體失穩時滑動速度中部較小，上部和下部速度較大；t=40 s滑坡體已經解體，滑坡體在下滑過程中積聚很大的動能，故此時分布有最大滑坡速度，撞擊對岸巖體后速度迅速改變方向，不斷堆積河谷，直至逐漸穩定，60 s左右速度逐漸減小為零。

圖10 滑坡體典型剖面g-g′不同時刻速度云圖Fig.10 Cloud pictures of velocity of a typical section named “g-g′” at different time

圖11顯示了系統的能量演化曲線(E為能量值，Emax為相應能量的最大值)，包括系統總動能歷時曲線、系統摩擦耗散能歷時曲線和破碎耗散能歷時曲線。在歷時20 s左右，系統的總動能和摩擦耗散能出現突變，邊坡開始失穩，滑坡體高速下滑，在40 s左右由于撞擊對岸巖體，動能開始迅速減小。在t=60 s，各能量基本保持不變，滑坡堆積體達到穩定狀態。從摩擦耗散能和破碎耗散能歷時曲線可知，隨著滑坡體的失穩和滑動，失效界面單元不斷增加，導致滑坡體由連續體向非連續體轉換，并最終形成松散的堆積體。

圖11 系統能量演化曲線Fig.11 Evolution curves of system energy

4 結 論

(1)將FDEM應用于紅石巖的三維邊坡失穩模擬，重現了堰塞體形成的全過程。在地震作用10 s左右，系統動能出現突變，邊坡開始失穩，地震總的作用時間20 s后，滑坡體高速下滑，進而受到對岸撞擊，不斷解體并堆積在河谷。通過FDEM分析能夠直觀地捕捉滑坡體各個時刻的滑動狀態、速度分布和滑坡體的最終堆積形態。

(2)將最大縱剖面和4個典型橫剖面的數值模擬的堆積形態與現場實際堆積形態進行對比，2者吻合程度較高。對比了表征堆積形態的參數，基本一致，進一步驗證了FDEM模擬邊坡失穩和滑動堆積形成堰塞體全過程的合理性。

(3)FDEM三維模擬結果表明，紅石巖堰塞體形成機制為地震誘發巖土體內部損傷和變形累積、結構面貫通、滑坡體高速下滑并伴隨破壞、解體、撞擊對面山體，堆積河道、不斷密實，形成堰塞體。

□

[1] Costa J E, Schuster R L. The formation and failure of natural dams [J]. Geological Society of America Bulletin, 1988,100(7):1 054-1 068.

[2] 徐文杰，陳祖煜，何秉順，等. 肖家橋滑坡堵江機制及災害鏈效應研究[J]. 巖石力學與工程學報，2010，29 (5)：933-942.

[3] 劉傳正，葛永剛，江興元，等. 魯甸地震紅石巖崩塌觸發機理分析[J]. 防災減災工程學報，2016，36(4)：601-607.

[4] 曹琰波，戴福初，許 沖，等. 唐家山滑坡變形運動機制的離散元模擬[J]. 巖石力學與工程學報，2011，30(增1)：2 879-2 887.

[5] Wu Jianhong, Lin Jeenshang, Chen Chaoshi. Dynamic discrete analysis of an earthquake-induced large-scale landslide[J]. International Journal of Rock Mechanics & Mining Sciences, 2009,(46):397-407.

[6] 焦玉勇，葛修潤，劉泉聲,等. 三維離散單元法及其在滑坡分析中的應用[J]. 巖土工程學報，2000，20(1)：101-104.

[7] Zhou Wei, Lai Zhiqiang, Ma Gang. Effect of base roughness on size segregation in dry granular flows[J].Granular Matter, 2016,(18):83.

[8] 鄔愛清，林紹忠，馬貴生，等. 唐家山堰塞壩形成機制DDA模擬研究[J]. 人民長江，2008，39(22)：91-95.

[9] Zhang Yingbin, Xu Qiang, Chen Guangqi, et al. Extension of discontinuous deformation analysis and application in cohesive-frictional slope analysis[J]. International Journal of Rock Mechanics & Mining Sciences, 2014,(70):533-545.

[10] Huang Da, Song Yixiang, Cen Duofeng, et al. Numerical modeling of earthquake-induced landslide using an improved discontinuous deformation analysis considering dynamic friction degradation of joints[J]. Rock Mech Rock Eng, 2016,49:4 767-4 786.

[11] Barrett P J. The shape of rock particles, a critical review [J]. Sedimentology, 1980,27(3):291-303.

[12] Marco Barla, Giovanna Piovano, Giovanni Grasselli. Rock slide simulation with the combined finite-discrete element method [J]. Geomechanics, 2012,12(6):711-721.

[13] Francesco Antolini, Marco Barla, Andrea Giorgetti, et al. Combined finite-discrete numerical modeling of runout of the torgio-vannetto di assisi rockslide in central italy [J]. International Journal of Geo-mechanics, 2016,16(6):1-16.

[14] Giovanna Piovano, Francesco Antolini, Marco Barla, et al. Continuum-discontinuum modelling of failure and evolution mechanisms of deep seated land-slides[C]∥6th Int. Conf. on Discrete Element Method, Colorado School of Mines, Golden, CO, USA, 2013:295-300.

[15] 常曉林，胡 超，馬 剛，等. 模擬巖體失效全過程的連續-非連續變形體離散元方法及應用[J]. 巖石力學與工程學報，2011，30(10)：2 004-2 011.

[16] Zhou Wei, Yuan Wei, Ma Gang, et al. Combined finite-discrete element method modeling of rockslides [J]. Engineering Computations, 2016,33(5):1 530-1 559.

[17] O K Mahabadi, A Lisjak, A Munjiza, et al. Y-Geo: new combined finite-discrete element numerical code for geome-chanical applications[J]. Geomechanics, 2012,12(6):676-688.

[18] 馬 剛，周創兵，常曉林，等. 巖石破壞全過程的連續-離散耦合分析方法[J]. 巖石力學與工程學報，2011，30 (12)：2 444-2 454.