土石混合體邊坡細觀特征對滑面形成影響研究*

陳 曉 石 崇 楊俊雄

（①河海大學，巖石力學與堤壩工程教育部重點實驗室，南京210098，中國）

（②河海大學巖土工程科學研究所，南京210098，中國）

0 引 言

土石混合體是由粒徑不同、強度較高的巖塊和軟弱的土體顆粒組成，一般存在于第四系松散堆積層中（Xu et al.，2009，2011；徐文杰等，2013）。其各種力學性質(zhì)及特性和普通勻質(zhì)土體差異較大，因此，由此類土體構成的巖土混合體邊坡的穩(wěn)定性影響因素較多，與普通勻質(zhì)邊坡的滑面發(fā)展和破壞機制有所不同（楊忠平等，2017；Wei et al.，2018；朱晨光等，2019）。這種類型邊坡是大型土木工程的潛在威脅，因此研究該類介質(zhì)滑動機理與穩(wěn)定性具有重要意義。

土石混合體作為一種廣泛分布在工程領域的重要地質(zhì)材料，其力學性質(zhì)及物理特性直接影響工程項目的安全穩(wěn)定。Shakoor et al.（1990）研究了塊石含量及尺寸對土石混合體抗壓特性的影響。Wang et al.（2018）借助CT技術對三軸壓縮試驗中土石混合試樣內(nèi)部的微結構變化進行了觀測。徐文杰等（2007a，2012）利用數(shù)值圖像技術充分研究了土石混合體細觀結構，將細觀結構定量化，為土石混合體離散元數(shù)值模型建立提供基礎；并進一步采用離散元及有限元方法模擬了土石混合體直剪試驗，發(fā)現(xiàn)塊石的存在使內(nèi)部細觀應力狀態(tài)有結構性特征，含石量及試樣尺寸均影響土石混合體的宏觀抗剪強度和剪脹性（徐文杰等，2007a，2007b，2014）。李世海等（2004）等建立三維離散元隨機計算模型，分析了土石混合體內(nèi)部應力場分布與土石級配及巖石塊度的關系。王環(huán)玲等（2019）通過顆粒流方法模擬了雙軸壓縮試驗，研究了含石率、土體孔隙率、不均勻系數(shù)對土石混合體強度特性的影響。

土石混合體邊坡廣泛分布在中國西南地區(qū)，其穩(wěn)定性對于水利水電等大型工程的建設尤為重要。由于土石混合體具有復雜的力學性質(zhì)，Vessia et al.（2017）提出采用多維度理論模型計算土石邊坡的安全系數(shù)；Medley et al.（2004）通過分析簡單的土石混合地質(zhì)模型，提出塊石含量是影響邊坡安全系數(shù)的重要因素；Napoli et al.（2018）建立了不同有限元邊坡模型，并結合極限平衡方法，研究了不同塊石含量邊坡的穩(wěn)定性。在國內(nèi)的相關研究中，徐文杰（2009），徐文杰等（2008）以數(shù)字圖像處理技術建立土石混合體邊坡細觀結構模型，并采用有限元強度折減法分析了塊石對于邊坡穩(wěn)定性的影響；李亮等（2017）運用Flac3D模型，探討了土石混合體邊坡內(nèi)部的塑性區(qū)擴展規(guī)律及土石界面強度參數(shù)對邊坡安全系數(shù)的影響。劉順青等（2019）提出了一種考慮不同含石率及塊石隨機分布的邊坡穩(wěn)定性分析方法。盡管目前在土石混合體邊坡方面已有較多研究，但相關的研究更多集中于有限元方法（邵帥等，2014；張森等，2016；龔健等，2017），采用顆粒離散元的相關研究相對欠缺；且在以上研究中，均忽略了塊石形態(tài)的影響，將塊石形態(tài)簡化為矩形或圓形。

顆粒離散元方法以顆粒（ball）和接觸（contact）作為基本單元，可以直觀地觀測到邊坡內(nèi)部裂隙擴展及邊坡失穩(wěn)過程（程升等，2018；鄭博寧等，2019）?；诖?，本文以某一典型土石混合體邊坡為例，建立土石混合體邊坡顆粒流細觀模型，依據(jù)室內(nèi)試驗數(shù)據(jù)標定土體材料參數(shù)。在合理細觀模型的基礎上，首先模擬均質(zhì)邊坡與土石混合體邊坡的失穩(wěn)過程，然后分析不同含石量邊坡滑面發(fā)展過程的差異，探討細觀結構對滑面發(fā)展及邊坡穩(wěn)定性的影響。

1 土石混合體細觀模型構造

1.1 土石混合體模型構建方法

土石混合體細觀特征的構建，重在骨架顆粒的外輪廓和微觀裂隙統(tǒng)計，在現(xiàn)場多采用統(tǒng)計窗方式。預先規(guī)劃好地質(zhì)統(tǒng)計窗，利用參照物拍照，然后借助AUTOCAD等工程軟件繪制輪廓，進一步分析其細觀特征。

數(shù)字圖像處理方法是一種被廣泛采用的土石混合體細觀特征提取方法，但是數(shù)字圖像方法的像素往往較高，如果將每一個像素都按照相應位置轉(zhuǎn)化為數(shù)值模型，則單元、節(jié)點多，計算工作量大。

實際上，圖像識別是一種有損識別方法，由于光照、陰影、拍照角度的差異，每一幅圖片中土石區(qū)分均有差異，因此可以只采用數(shù)字圖像識別的骨架顆粒輪廓線構造土石混合體的細觀特征。故本文提出如下塊石輪廓構造方法：

如圖1a，將數(shù)字圖像識別的骨架顆粒輪廓線作為邊界線，對每一條邊界進行讀?。╬olyline），處于任一多段線內(nèi)的像素屬于塊石、而不在任一多段線內(nèi)的像素屬于膠結物。這樣即可將每一像素的性質(zhì)（土或石）區(qū)分開，并借助這些多段線數(shù)據(jù)開展顆粒粒徑、形狀等信息的統(tǒng)計。最終，將這些信息導入PFC中作為模型生成的邊界（圖1b）。

圖1 模型邊界（徐文杰等，2008）Fig.1 Model boundaries（Xu et al.，2008）

1.2 土石混合體邊坡模型的建立

首先，按照圖1b所示的邊界生成墻體（wall），再在邊界內(nèi)采用基本的顆粒生成命令形成完整的顆粒（ball）體系。生成的初始模型中顆粒是松散的，顆粒間的間隙較大，且應力分布極不均勻。此時，我們需要利用伺服機制（Cundall et al.，2008；Itasca.Consulting Group Inc.，2008）使初始模型達到一種密實且應力分布均勻的狀態(tài)。

對邊界wall施加一定的速度以表示恒定約束力，墻體移動產(chǎn)生的作用于墻體上的力可以表示為：

式中：kn為平均接觸剛度；ΔUn為單位計算步內(nèi)位移增量；Nc為與wall接觸顆?？倲?shù)目；u·為墻體的運動速度；Δt為單位計算步時長。作用在墻體上的應力的變化值可以采用以下公式計算：

式中：A為與墻體接觸作用的面積，二維情況下，通常取墻體的長度；于是可得wall的速度與單位計算步內(nèi)應力增量的關系：

式中：G為伺服系數(shù)，不同計算步下，與墻體接觸的總數(shù)目不同，每一計算步更新得到不同的伺服系數(shù)。伺服的過程就是通過調(diào)整墻體的運動使墻體與模型之間的接觸力不斷趨向于目標應力值的大小，墻體的接觸力與目標應力值的差值可以表示為：

式中：σm為當前計算時步監(jiān)測得到的墻體接觸應力；σr為伺服結果需要實現(xiàn)的最終應力。由當前計算步的應力差值，根據(jù)式（5）和式（6）可以計算得到下一計算時步的墻體運動速度，不斷迭代更新，直至Δσ′小于某一限定值。圖2a為伺服前的模型圖，圖2b為伺服結束時的模型圖。可以看出伺服后，模型的接觸（contact）數(shù)目及分布均勻程度均更接近真實狀態(tài)。

圖2 接觸分布形態(tài)Fig.2 The distribution of contact

在得到穩(wěn)定的初始模型后，根據(jù)圖1b所示的多邊形邊界，獲取所有顆粒的位置信息，判斷顆粒位置與多邊形邊界的關系，當判斷位于多邊形內(nèi)時，則確定該顆粒屬于巖石材料，位于同一多邊形區(qū)域內(nèi)的顆粒封裝形成塊石簇（clump）；并將邊界內(nèi)的其他顆粒定義為基質(zhì)土材料。

通過以上模型構造方法，建立了如圖3所示的顆粒離散元邊坡數(shù)值模型。模型尺寸如圖3a所示。模型共包含222i633個顆粒（ball），顆粒（ball）半徑為0.5～1 cm?？紤]到后續(xù)模擬滑坡過程的需要，在模型的左邊界和下邊界設置了邊界條件，固定邊界顆粒的運動，如圖3中紅色區(qū)域所示。

1.3 模型細觀參數(shù)標定

圖3 邊坡細觀結構模型Fig.3 Numerical model of slope mesoscopic structure

在顆粒離散元方法中，采用圓球（ball）模擬材料顆粒，并通過接觸（contact）模擬顆粒之間的相互作用，通過為接觸（contact）賦予不同的細觀參數(shù)模擬材料的不同力學性質(zhì)。故細觀參數(shù)的標定是顆粒離散元數(shù)值模擬中的重要環(huán)節(jié)，合理細觀參數(shù)的選取是數(shù)值模擬結果可信的基礎（韓振華等，2019）。

在顆粒離散元中存在不同的接觸本構模型（Itasca Consulting Group Inc.，2008），在本文中分別采用線性平行黏結模型（Linear Parallel Bond Model）和線性接觸黏結模型（Linear Contact Bond Model）模擬巖石和土的力學性質(zhì)。由于在滑坡過程中，塊石不會發(fā)生破壞，故巖石的參數(shù)根據(jù)經(jīng)驗取值。對土的細觀參數(shù)進行了標定。

在接觸細觀參數(shù)中，對宏觀力學行為影響顯著的主要包括黏結強度和有效模量。經(jīng)相關研究發(fā)現(xiàn)（Potyondy，2004；Cho et al.，2007），接觸有效模量直接影響材料的宏觀楊氏模量；法向切向剛度比控制材料的泊松比；黏結的拉伸與剪切強度比與試樣的破壞模式相關聯(lián)?；谝陨弦?guī)律，并根據(jù)土體的室內(nèi)試驗曲線，如圖4中虛線所示，首先確定強度參數(shù)和模量參數(shù)，保證不同圍壓下，試樣的抗壓強度與試驗曲線基本一致。而后控制拉伸與剪切強度比保證試樣表現(xiàn)出較強的塑性破壞特征。在標定過程中發(fā)現(xiàn)，當接觸抗拉強度與剪切強度基本相等時，試樣表現(xiàn)為明顯的脆性破壞；隨著抗拉強度與抗剪強度比值的增加，試樣的塑性破壞特征增強；并且在低強度試樣中，有效模量參數(shù)不僅影響試樣的宏觀楊氏模量，同時會對試樣的強度產(chǎn)生影響。最終，標定得到的細觀參數(shù)如表1所示，不同圍壓下的數(shù)值試驗曲線如圖4中實線部分所示。

除以上標定的強度特性參數(shù)外，在模擬滑坡過程中，阻尼參數(shù)的選取也尤為重要。在顆粒離散元中，共包含兩種阻尼參數(shù)：局部阻尼和黏性阻尼。前者可以加快數(shù)值模擬的計算平衡，后者則反映顆粒碰撞過程中的能量消散。為了反映滑坡的真實過程，在模擬滑坡過程中不考慮局部阻尼作用，法向和切向黏性阻尼比分別取0.4和0.2。

圖4 不同圍壓下土體應力－應變曲線Fig.4 Stress-strain curve of soil under different confining pressure

表1 土石混合體主要細觀力學參數(shù)Table 1 The main mesomechanical parameters of soil-rock mixtures

1.4 邊坡穩(wěn)定性計算方法

強度折減法是由Dawson et al.（1999），Griffiths et al.（1999）提出的“包圍方法”自動尋找安全系數(shù)。將宏觀強度參數(shù)，如黏聚力和摩擦角，除以折減系數(shù)K作為新的材料參數(shù)，計算不同的折減參數(shù)并判斷系統(tǒng)是否處于平衡狀態(tài)，監(jiān)測穩(wěn)定與不穩(wěn)定的臨界狀態(tài)，從而得到模型破壞時的安全系數(shù)。折減后ck、φk為：

在顆粒離散元中，對巖土體賦予的不同接觸模型的細觀力學參數(shù)并非宏觀參數(shù)，因此需對該微觀參數(shù)進行折減。其主要實現(xiàn)方法是通過對細觀張拉強度、剪切強度以及摩擦系數(shù)進行參數(shù)折減，當模型的不平衡力比（ratio）能夠達到1e－5時，則認為模型能夠達到穩(wěn)定狀態(tài)，不平衡力比不能達到要求時，則認為模型為非穩(wěn)定狀態(tài)。具體折減參數(shù)如下：

2 模型計算結果

為了研究土石混合體細觀結構對邊坡滑面形成機制的影響，分別模擬純土體邊坡和土石混合體邊坡的失穩(wěn)過程，計算兩種邊坡的安全系數(shù)，記錄邊坡內(nèi)部應力變化與顆粒間的傳力機制，以及裂隙和滑面發(fā)展演變過程。對比兩種邊坡破壞機制的差異，分析土石混合介質(zhì)對邊坡滑面形成的影響。

2.1 初始邊坡模型

分別將細觀參數(shù)賦值于邊坡細觀模型，并在自重作用下平衡。圖5是純土體邊坡自重平衡下的接觸力鏈分布圖?？梢钥闯鲈谥亓ψ饔孟拢麦w表層的力鏈分布稀疏，且坡體表層的接觸力明顯要小于坡體內(nèi)部的接觸力，坡腳位置力鏈密集區(qū)與稀疏區(qū)有較為明顯的分界線，可以判定這是邊坡的潛在滑面位置。

圖6是土石混合體邊坡在自重作用下平衡后的接觸力鏈分布圖，從圖中可以看出邊坡內(nèi)部接觸力鏈的分布情況與圖5中純土體邊坡的力鏈分布明顯不同，圖6中力鏈的分布明顯比純土質(zhì)邊坡的力鏈分布復雜，尤其是在塊石周圍存在明顯的剪切環(huán)，由于塊石的存在，邊坡內(nèi)部接觸力鏈的分布在遇到石塊時會繞開石塊，形成沿塊體邊緣的剪切閉環(huán)傳力路徑。且淺層的力鏈分布明顯比純土質(zhì)邊坡的接觸力鏈分布密集。

2.2 邊坡模型穩(wěn)定性計算

圖5 純土體邊坡接觸力鏈分布Fig.5 The distribution of force chain of pure soil slope

圖6 土石混合體邊坡接觸力鏈分布Fig.6 The distribution of force chain of soil-rock slope

為了準確描述該土石混合體邊坡的穩(wěn)定性，在邊坡模型的淺層坡體內(nèi)選取了5個監(jiān)測點（圖7），分別監(jiān)測并記錄不同強度折減系數(shù)下5個測點的位移值，進而求得土石混合體邊坡的安全系數(shù)。監(jiān)測結果如圖8所示，其中圖8a～圖8d分別對應強度折減系數(shù)為1.0、1.05、1.10、1.15 4個工況下的監(jiān)測點位移變化情況。從位移曲線中可以看出，強度折減系數(shù)為1.0時，5個監(jiān)測點的位移曲線均在某一個水平上下浮動，且位移值均小于1 cm，故邊坡處于穩(wěn)定狀態(tài)；強度折減系數(shù)為1.05時，測點1的位移明顯較大，且其余測點位移也開始有上升的趨勢，但最終位移均趨近于穩(wěn)定，即邊坡仍處于穩(wěn)定狀態(tài)；當強度折減系數(shù)達到1.10時，位移曲線形狀發(fā)生突變，不再收斂于某一固定值。測點的位移均隨計算時間增加而增加，且測點1的位移始終是最大的，故邊坡失穩(wěn)，發(fā)生滑坡；強度折減系數(shù)為1.15時，邊坡的滑動更加明顯，表層完全錯動，即形成滑坡。從整個過程中可以得知：此邊坡的安全系數(shù)在1.05～1.10之間。由于坡腳位移最大且最早發(fā)生，因此可推斷此滑坡屬于牽引式滑坡。采用同樣的方法，經(jīng)過計算發(fā)現(xiàn)，此純土邊坡的安全系數(shù)約為1.05。故此土石混合體邊坡的安全系數(shù)要略高于相同條件下純土體邊坡的安全系數(shù)。

圖7 邊坡位移監(jiān)測點Fig.7 Monitoring point of slope displacement

圖8 不同強度折減系數(shù)下監(jiān)測點位移Fig.8 Displacement of monitoring points under different strength reduction coefficient

2.3 邊坡模型滑面特征

從以上計算結果可以看出，土石混合體邊坡模型與純土體邊坡模型存在明顯差異，為了分析土石混合體細觀特征對滑面形成機制的影響，分別模擬純土體邊坡和土石混合體邊坡的滑坡過程。在數(shù)值模型內(nèi)部，微觀裂隙的衍生過程即為邊坡模型的失穩(wěn)破壞過程。通過追蹤模型內(nèi)部的裂隙擴展路徑，即可表征邊坡內(nèi)的滑面形成機制。圖9、圖10分別記錄了純土體邊坡模型和土石混合體邊坡模型計算25萬步（圖9a）、50萬步（圖9b，圖10b）、75萬步（圖9c，圖10c）、100萬步（圖9d，圖10d）時的裂隙發(fā)育形態(tài)。

圖9 純土體邊坡失穩(wěn)過程Fig.9 The landslide process of pure soil slope

圖10 土石混合體邊坡失穩(wěn)過程Fig.10 The landslide process of soil-rock mixture slope

從圖9中可以看出，邊坡表層的顆粒由于無墻體約束，在強度折減后黏結力降低，發(fā)生淺層的顆?；瑒印Ｓ嬎愕?0萬步時邊坡體內(nèi)開始發(fā)生破壞，隨著計算步數(shù)的增加，邊坡內(nèi)部滑面從底部尖端向上發(fā)展。與圖10相比，純土體邊坡滑面的形成過程與土石混合體邊坡模型滑面的形成過程存在顯著不同。純土體邊坡滑面發(fā)展方向單一，微觀破壞裂隙的分布沿滑面的發(fā)展方向相對分散。相較于純土體邊坡內(nèi)微觀破壞的離散性，土石混合體邊坡內(nèi)的微觀破壞局部集中，且由于塊石的存在，坡體內(nèi)部的破壞形態(tài)明顯發(fā)生改變，滑面的位置更趨近于表層，在破壞速度上，也明顯低于純土體邊坡，相同計算時間內(nèi)，純土體邊坡模型內(nèi)的微裂隙數(shù)目約為土石混合體邊坡模型微觀裂隙數(shù)目的1.5倍，故塊石的存在明顯改變了邊坡內(nèi)部的細觀結構，增強了坡體抵抗破壞的能力，提高了邊坡的穩(wěn)定性。同時，從圖10b和圖10c中可以發(fā)現(xiàn)，破壞滑面明顯沿著兩個方向發(fā)展，一條主滑面，一條次滑面，兩條滑面均穿過塊石之間的間隙，且裂隙的分布在其發(fā)展方向上相對集中。圖10d中邊坡體內(nèi)發(fā)育形成一條貫通的可視滑面，最終的滑面位置存在于塊石之間的間隙。

圖11為與圖1所示邊坡相似的一處土石混合體邊坡，地質(zhì)情況相似，滑坡的形狀也接近，此圖片為邊坡滑坡后的情形，從圖中可以看出滑面的位置與滑坡體形態(tài)均與上述模擬結果基本吻合。

圖11 現(xiàn)場滑面Fig.11 The sliding surface on site

2.4 不同含石率邊坡滑面破壞特征

圖12 不同含石率下邊坡滑面擴展過程Fig.12 The spreading process of slope sliding surface under different rock content

為了進一步論證土石混合體細觀特征對邊坡破壞機理的影響，分別進行了土石混合體邊坡中含石量為10%、20%、30%、50%的失穩(wěn)過程模擬（圖12）。含石率為10%時，滑坡體表現(xiàn)為典型的圓弧狀破壞形式，滑面表現(xiàn)為自坡腳向坡頂貫穿的圓弧形態(tài)。含石率為20%時，在坡體底部，微裂隙沿塊石間的薄弱位置發(fā)育。由于在坡體頂部，塊石分布密集，滑面形成受阻，在塊石周圍散布少許微裂隙，且未完全貫通。故在邊坡體內(nèi)形成一條不完整的圓弧狀滑動帶。含石率為30%時，邊坡破壞形式發(fā)生明顯改變，裂隙在坡腳位置集中，且基本分布在邊坡淺層附近，邊坡的破壞形式由整體滑動變?yōu)槠履_局部滑動，在坡腳的表層形成一條局部的滑面。含石率為50%時，微裂隙的分布基本密布在坡腳位置，在坡腳的位置形成一條局部的圓弧滑面，其滑面半徑約是含石率10%時的圓弧狀滑面半徑的1／4。就不同含石率下的滑坡體體積而言，含石率10%與20%時，滑坡體體積基本相同，無較大差異。含石率30%時，滑坡體體積明顯減少，約為含石率10%時滑坡體的1／8。含石率繼續(xù)增加至50%時，滑坡體體積繼續(xù)減少，約為含石率30%時滑坡體體積的1／2。

故土石混合體細觀特征對于邊坡穩(wěn)定性及滑面形成過程影響顯著。不同含石率下，滑面的形態(tài)各異，尤其含石率超過30%時，邊坡的失穩(wěn)形式發(fā)生明顯變化，滑面的形態(tài)也由整體向局部變化。

3 結 論

本文基于土石混合體邊坡實例開展顆粒離散元數(shù)值模擬研究，探討了純土體邊坡和土石混合體邊坡的穩(wěn)定性、變形承力機制和滑面破壞機理的差異，得到以下主要結論：

（1）塊石的存在直接影響土體邊坡內(nèi)部接觸力鏈的分布。相較于純土體邊坡，土石混合體邊坡內(nèi)部的接觸力鏈分布更加復雜，在塊石周圍形成剪切閉環(huán)傳力路徑。

（2）在純土體邊坡內(nèi)部，形成單一的邊坡滑面，且在滑面方向上，微裂隙呈發(fā)散分布態(tài)勢。在土石混合體邊坡內(nèi)，存在明顯的繞石現(xiàn)象。由于塊石的存在，滑面同時向兩個方向擴展，且微裂隙發(fā)展相對集中。

（3）塊石的存在會改變邊坡內(nèi)部的細觀結構，增強坡體抵抗破壞的能力，降低邊坡失穩(wěn)的速度。且在不同含石率下，滑面的形態(tài)各異，尤其含石率超過30%時，邊坡的失穩(wěn)形式發(fā)生明顯變化，滑面的形態(tài)也由整體向局部變化。

可見，土石混合體邊坡與純土體邊坡存在明顯差距，塊石的存在會改變邊坡內(nèi)力的分布與傳遞及滑面的形成機制，且對于提高邊坡穩(wěn)定性有明顯作用。因此，在對土石混合體邊坡進行研究時，必須充分考慮內(nèi)部塊石的大小、分布，考慮塊石對于邊坡力學性質(zhì)的影響。