邊坡穩(wěn)定性分析的離散元極限平衡法研究

陳 頡，曾亞武

(1.武漢大學(xué) 土木建筑工程學(xué)院， 湖北 武漢 430072； 2.湖北工業(yè)大學(xué) 計(jì)算機(jī)學(xué)院， 湖北 武漢 430065)

極限平衡法是最古老的邊坡穩(wěn)定性分析方法，其基本理論是假定邊坡為剛體，可能沿著剛體中的某一滑動面滑動，計(jì)算時(shí)將可能滑動部分切分成條塊，并對條塊之間的作用力進(jìn)行合理的假設(shè)，通過迭代試算獲得邊坡的安全系數(shù)。根據(jù)條間力的假定不同，可將極限平衡法分為瑞典圓弧法、畢肖普法、簡布法、斯賓塞法、摩根斯坦-普萊斯法、莎爾瑪法、不平衡推力法，以及通用條分法等[1]。這些方法都是對滑動面和條間力的三要素(力的大小、方向、作用點(diǎn))進(jìn)行了合理假設(shè)，降低了問題復(fù)雜度并簡化了計(jì)算。但是，幾乎每種方法都有遇到失效或者迭代計(jì)算不收斂的情況[2]，而且隨著計(jì)算模型精密程度提高，計(jì)算量也大為提高。

有限元法[3]是最早引入邊坡穩(wěn)定分析的一種數(shù)值方法，其中以基于強(qiáng)度折減的有限元法[4-5]使用最為廣泛，其基本思想也與極限平衡法中的強(qiáng)度折減相類似，不同的是該方法考慮了邊坡巖土體的受力變形特征。實(shí)際的使用過程中，如何判定邊坡的極限狀態(tài)、確定準(zhǔn)確的安全系數(shù)仍然是值得研究的問題。此外，在計(jì)算精度要求較高的情況下，有限元法的計(jì)算量較大，尤其在進(jìn)行可靠度分析時(shí)，計(jì)算量巨大。為此李典慶等[6]將隨機(jī)模型引入到有限元法進(jìn)行邊坡穩(wěn)定可靠度分析時(shí)，結(jié)合了極限平衡法思想，有效的減少了計(jì)算時(shí)間，是一種有益的嘗試。

離散元法引入邊坡穩(wěn)定分析中以后[7]，很快顯示出對傳統(tǒng)的極限平衡法以及有限元法的互補(bǔ)性優(yōu)勢。由于離散元法的物理模型簡單，數(shù)學(xué)理論較為完備，在將邊坡的宏觀土力學(xué)物理特征參數(shù)合理的轉(zhuǎn)化為細(xì)觀顆粒物理力學(xué)參數(shù)后，能比較準(zhǔn)確的預(yù)測邊坡在各種工況條件下的運(yùn)動趨勢。但是，離散元法計(jì)算邊坡穩(wěn)定性一般都需要消耗大量的計(jì)算時(shí)間，而且計(jì)算時(shí)間隨著顆粒數(shù)量呈指數(shù)級增長。即使在計(jì)算機(jī)硬件性能不斷提高的情況下，離散元軟件的計(jì)算時(shí)間并沒有顯著降低，其本質(zhì)原因是離散元模型是一種在時(shí)間和空間上串行的計(jì)算模型。針對離散元模型計(jì)算效率不高的缺陷，結(jié)合計(jì)算機(jī)硬件和軟件的優(yōu)化方案也是一個(gè)當(dāng)前研究的熱點(diǎn)。常新正[8]將GPU并行計(jì)算引入離散元計(jì)算中，取得了一定的效果，但是該方案主要是針對顆粒在空間的碰撞檢測過程進(jìn)行優(yōu)化，計(jì)算機(jī)效率提高程度受到了制約。田健[9]將Open-CL并行計(jì)算引入離散元計(jì)算中，使得加速過程與計(jì)算機(jī)硬件的相關(guān)性大大降低，適用性大為提高，但優(yōu)化的主要方向依然是加速顆粒在空間碰撞的檢測過程。

考慮到離散元方法和極限平衡方法各自的特點(diǎn)，本文提出了一種新的邊坡穩(wěn)定可靠度分析方法——離散元極限分析法。該方法結(jié)合了離散元法和極限平衡法的優(yōu)點(diǎn)，模型離散化方面采用了離散元的思想，而在計(jì)算過程以及邊坡穩(wěn)定性判定方案的設(shè)計(jì)上運(yùn)用了極限平衡法的主要思想，因而大大降低了計(jì)算過程的復(fù)雜程度，節(jié)約計(jì)算時(shí)間，且保證計(jì)算精度基本不變，可以為類似的工程案例提供一種參考。

1 邊坡穩(wěn)定性分析的離散元極限平衡法

離散元極限平衡法主要分為四步實(shí)現(xiàn)：第一步是將實(shí)體邊坡截面離散化為許多個(gè)直徑較大的圓盤，同時(shí)保證邊坡的土力學(xué)性質(zhì)，比如土密度，力學(xué)邊界的剛度，空隙率等，圓盤與周邊圓盤以及圓盤與剛性邊界保持接觸關(guān)系；第二步是利用極限平衡理論，計(jì)算出各個(gè)圓盤接觸點(diǎn)上的法向壓力與切向靜摩擦力的大小，同時(shí)計(jì)算出切向摩擦力與法向壓力的比值——將該比值定義為相對摩擦系數(shù)，并按照相對摩擦系數(shù)的大小對圓盤接觸點(diǎn)進(jìn)行排序，進(jìn)而勾勒出潛在的可能滑動帶；第三步是將可能滑動帶以外的區(qū)域設(shè)置成剛性塊體，僅將潛在滑動帶區(qū)域用半徑更小的圓盤進(jìn)行覆蓋，以減少圓盤的數(shù)量，利用剛性塊體和小圓盤的靜力平衡，再次計(jì)算出相對摩擦系數(shù)較大的新接觸點(diǎn)集合，勾勒出最危險(xiǎn)滑動線(面)；第四步是按照兩次計(jì)算出的滑動帶(面)中各個(gè)圓盤接觸點(diǎn)的相對摩擦系數(shù)的分布律，綜合考慮這些接觸點(diǎn)摩擦系數(shù)在整體安全系數(shù)計(jì)算中的權(quán)重值，按照正態(tài)分布模型迅速的計(jì)算出在極限平衡滑動帶區(qū)域內(nèi)在統(tǒng)計(jì)意義上的最小邊坡穩(wěn)定系數(shù)。針對土質(zhì)邊坡的計(jì)算結(jié)果表明，該方法計(jì)算的精度是可以保證的，與傳統(tǒng)的極限平衡法相比誤差很小，比單獨(dú)采用離散元法或者有限元法進(jìn)行計(jì)算時(shí)間更短，效率更高。

1.1 實(shí)體邊坡的離散化模型

邊坡離散化的主要目的是將邊坡按照有利于迅速找到潛在滑動區(qū)域帶，有利于加速后續(xù)的穩(wěn)定性計(jì)算過程的原則，劃分為規(guī)則的若干區(qū)域，在保證土力學(xué)屬性以及力學(xué)模型的前提下降低力學(xué)計(jì)算的復(fù)雜程度。在綜合比對了諸多極限平衡法和離散元法的計(jì)算實(shí)例后[10]，提出一種綜合兩者優(yōu)勢的離散化模型，分兩步計(jì)算出潛在滑動線(面)。第一步用直徑較大的相互接觸良好的圓盤覆蓋整個(gè)計(jì)算剖面區(qū)域，如圖1所示，并按照靜力計(jì)算結(jié)果初步確定潛在滑移帶區(qū)域；第二步在第一步計(jì)算結(jié)果的基礎(chǔ)上采用較小直徑的相互接觸良好的圓盤覆蓋潛在滑動帶區(qū)域，而帶外則采用剛性塊體替代，如圖2所示，再次計(jì)算并細(xì)化潛在滑移區(qū)域帶，確定最危險(xiǎn)滑動線(面)。

計(jì)算時(shí)，首先采用如圖1所示的離散模型，計(jì)算出各圓盤接觸點(diǎn)的法向壓力和切向摩擦力，以及相對摩擦系數(shù)(切向摩擦力和法向壓力之間的比值)，按照相對摩擦系數(shù)大小排序，初步定位出潛在滑動帶區(qū)域；然后將邊坡計(jì)算剖面劃分為潛在滑動帶和潛在滑動帶以外的區(qū)域，滑動帶以外的區(qū)域按照極限平衡的思想整體化為一個(gè)剛性塊體，而潛在滑動帶區(qū)域則用直徑較小的圓盤覆蓋，整個(gè)過程中要保持圓盤與圓盤、圓盤與剛體之間的相切接觸。根據(jù)需要，可以設(shè)置自由邊界(沒有約束的邊界)，如圖1所示。

圖1 邊坡計(jì)算剖面的第一次離散化(整體覆蓋)

圖2邊坡計(jì)算剖面的第二次離散化(分區(qū)覆蓋)

這種計(jì)算方法可以減少圓盤的數(shù)量，簡化計(jì)算模型進(jìn)而加速計(jì)算。需要指出的是，在采用圓盤進(jìn)行覆蓋時(shí)，應(yīng)該采用穩(wěn)定的三角形排列方式，而不是正方形排列方式，因?yàn)榍罢咴诹W(xué)上是穩(wěn)定的，后者不穩(wěn)定。

在圓盤覆蓋的過程中，對于不同的邊界，圓盤的放置有一些規(guī)則。對于潛在滑動帶與剛體之間的分界，圓盤要盡量保持與剛體的相切接觸，對于垂直邊界和水平邊界也應(yīng)如此，這些一般都可以通過調(diào)整圓盤直徑和剛體尺寸來實(shí)現(xiàn)。針對個(gè)別不能嚴(yán)格相切的區(qū)域，采取一個(gè)折中方法使得相互接觸的條件能夠滿足：用折線邊界與鄰近圓盤實(shí)質(zhì)性接觸，此法相對于圓盤填充方法，可以在不損失計(jì)算精度的前提下大量節(jié)省計(jì)算時(shí)間[11-12]。

1.2 離散化后各單元接觸點(diǎn)作用力計(jì)算

離散后的模型各單元(圓盤或剛體)之間沒有粘結(jié)，只有接觸。對于圓盤單元來說，有兩種接觸，一種是圓盤與圓盤之間的接觸，另一種是圓盤與剛體邊界或界面之間的接觸。根據(jù)本文約定的圓盤排列規(guī)則，一個(gè)圓盤最多與其他六個(gè)圓盤接觸，且此圓盤屬于“中間”圓盤，即只與圓盤接觸，如圖3(a)所示；除此以外，與剛體邊界或界面接觸或相鄰的圓盤，稱作“邊界”單元，其接觸數(shù)應(yīng)小于六個(gè)，如圖3(b)所示。

圖3圓盤單元的接觸

對于剛體來說，可以與若干個(gè)圓盤單元接觸。由于要保證邊界單元盡可能與剛體接觸，剛體邊界可以是直線，也可以是折線，理論上也可以是曲線，本文為簡化計(jì)算，不考慮曲線邊界情況。如圖4(a)所示為剛體直線邊界與多個(gè)圓盤的接觸，圖4(b)所示為剛體折線邊界與多個(gè)圓盤的接觸。

圖4剛體與圓盤的接觸

所有單元除了重力外，還受到相鄰單元的作用力，該作用力作用在相應(yīng)的接觸點(diǎn)上，可分解為法向作用力和切向摩擦力，且法向作用力只能為壓力，不能為拉力。單元受力分析見圖5。

圖5單元受力示意圖

假定所有單元(圓盤和剛體)在重力、接觸力和邊界約束力的作用下保持靜力平衡狀態(tài)，即假定單元的強(qiáng)度及單元之間的摩擦系數(shù)無窮大，單元不會破壞，也不會產(chǎn)生相對運(yùn)動。按照靜力平衡要求，寫出每個(gè)單元水平方向、垂直方向的受力平衡方程和相對于單元形心的力矩平衡方程，有

(1)

(2)

(3)

式中：i為單元接觸編號，i=1，2，……，6；j為單元編號，j=1，2，……，n，其中n為單元總數(shù)。

對于圓盤單元而言，其形心為圓盤的圓心，重力作用于圓心，法向力都指向圓心，所以都不產(chǎn)生力矩，所有力矩均由接觸點(diǎn)的摩擦力產(chǎn)生，且力臂都等于圓盤半徑，因此圓盤單元的力矩平衡方程式(3)可以寫成摩擦力的平衡方程，見式(4)。

(4)

將所有單元的靜力平衡方程聯(lián)立，得到整個(gè)計(jì)算模型的代數(shù)方程組，其中方程個(gè)數(shù)為3n，未知數(shù)個(gè)數(shù)為2m，m為單元與單元之間、單元與邊界之間的接觸點(diǎn)總數(shù)。由于每個(gè)單元至少與臨近單元或邊界有2個(gè)接觸點(diǎn)(否則不穩(wěn)定)，即m>2n，因此未知數(shù)個(gè)數(shù)2m應(yīng)大于方程總數(shù)3n，即離散模型是超靜定系統(tǒng)，需要補(bǔ)充其他條件或假定才能求解。

因此，為了使得整個(gè)離散模型的代數(shù)方程組可解，參照極限平衡法思想，做出如下三個(gè)假定：

(1) 圓盤與圓盤之間在水平方向的接觸點(diǎn)具有相離的趨勢，因此假定其法向作用力和切向摩擦力為零。

(2) 假定圓盤與垂直剛性邊界之間的法向壓力等于土力學(xué)中的靜止土壓力，即

N=(1-sinφ)×∑G

(5)

式中：N為圓盤與垂直剛性邊界之間的法向壓力(即水平作用力)；φ為內(nèi)摩擦角；G為接觸圓盤及其正上方所有圓盤的重力之和；并且假定圓盤與垂直剛性邊界接觸點(diǎn)的摩擦力為零(光滑接觸)。

(3) 假定圓盤與水平剛性邊界之間的壓力等于該圓盤及其上部相應(yīng)范圍內(nèi)圓盤的重力之和。

通過上述三個(gè)假定，可以使整個(gè)離散模型的代數(shù)聯(lián)立方程組中的未知量大大減少，變成靜定問題，從而可解，求出單元與單元之間、單元與邊界之間接觸點(diǎn)的作用力。在計(jì)算單個(gè)圓盤的受力時(shí)，先從邊界的單元開始計(jì)算，進(jìn)而通過迭代求解得到周圍的圓盤的受力情況。

1.3 邊坡潛在滑動帶的搜索

在獲得了邊坡離散體系中每個(gè)接觸點(diǎn)的法向壓力和切向摩擦力后，就可以利用接觸點(diǎn)的作用力來確定邊坡潛在的滑動帶區(qū)域。

首先，將邊坡離散體系中接觸點(diǎn)的最大允許切向力與法向力的比值定義為極限摩擦系數(shù)μmax，根據(jù)文獻(xiàn)[13]中關(guān)于離散元細(xì)觀參數(shù)與土力學(xué)宏觀參數(shù)的實(shí)驗(yàn)結(jié)果，該極限摩擦系數(shù)與土體的內(nèi)摩擦角φ之間存在如下的函數(shù)關(guān)系

μmax=φ/37.292-0.1044

(6)

文獻(xiàn)[12]的實(shí)驗(yàn)顯示極限摩擦系數(shù)μmax與土體黏聚力c關(guān)系不大，而且在數(shù)值上極限摩擦系數(shù)大于內(nèi)摩擦角的正切值，說明實(shí)驗(yàn)結(jié)果中黏聚力的作用已經(jīng)體現(xiàn)在了極限摩擦系數(shù)的計(jì)算結(jié)果中。而黏聚力為零的沙質(zhì)土情況下的極限摩擦系數(shù)與其自身的級配、孔隙比以及含水率等有著密切關(guān)系，不在本文討論范圍內(nèi)。

其次，將計(jì)算結(jié)果中每個(gè)接觸點(diǎn)的切向力和法向力的比值定義為相對摩擦系數(shù)μre，即

μre=f/N

(7)

式中：f為接觸點(diǎn)的切向摩擦力絕對值；N為接觸點(diǎn)的法向壓力值。

根據(jù)整個(gè)離散體系靜力平衡方程，可求解出所有接觸點(diǎn)的作用力，進(jìn)而獲得相對摩擦系數(shù)在邊坡二維平面的一個(gè)分布。計(jì)算出相對摩擦系數(shù)可能大于、等于或小于極限摩擦系數(shù)，顯然，相對摩擦系數(shù)大于極限摩擦系數(shù)的接觸點(diǎn)，在受力過程中會發(fā)生破壞；相對摩擦系數(shù)等于極限摩擦系數(shù)的接觸點(diǎn)，處于極限平衡狀態(tài)；相對摩擦系數(shù)小于極限摩擦系數(shù)的接觸點(diǎn)，則處于安全狀態(tài)。也就是說，接觸點(diǎn)相對摩擦系數(shù)的大小以及它們與極限摩擦系數(shù)的接近程度，反映了接觸點(diǎn)的安全狀態(tài)，而整個(gè)離散體系接觸點(diǎn)相對摩擦系數(shù)的分布，則反映了整個(gè)邊坡的安全程度。為此，進(jìn)一步定義每個(gè)接觸點(diǎn)的相對摩擦系數(shù)與極限摩擦系數(shù)之比為極限滑移比K：

K=μre/μmax

(8)

顯然極限滑移比K值越大，則該接觸點(diǎn)發(fā)生滑動破壞的可能性越大。將所有接觸點(diǎn)的極限滑移比按照數(shù)值大小排序，并在整個(gè)邊坡離散體系中搜索一系列極限滑移比較大，相互毗鄰的圓盤接觸點(diǎn)。這些圓盤在空間分布上一般是從坡面的某一位置一直連續(xù)的延伸到坡頂。如圖6所示，給出了一個(gè)由極限滑移比較大的圓盤組成的潛在滑動帶。

圖6通過極限滑移比大小排序搜索獲得的邊坡潛在滑動帶

只要圓盤單元直徑相對于邊坡尺度足夠小，連接潛在滑動帶中相鄰圓盤中心的折線即可以認(rèn)為是極限平衡法中的潛在滑動線(面)。但由于圓盤單元直徑過小，會導(dǎo)致整個(gè)離散體系單元數(shù)量過多，計(jì)算量過大，難以快速獲得精確結(jié)果。因此，采用本文提出的計(jì)算方案，先以直徑較大、數(shù)量較少的圓盤單元來建立邊坡的離散體系，給出邊坡潛在滑移帶的初步范圍；然后再在初步計(jì)算的基礎(chǔ)上，將潛在滑移帶區(qū)域用更小的圓盤單元來替代，而該區(qū)域外則采用剛性塊體來替代，以減少離散單元體系的單元數(shù)量，再進(jìn)行精確計(jì)算，確定邊坡最危險(xiǎn)滑動線(面)。在第一次初定潛在滑移帶范圍后，第二次用小圓盤加剛性塊體模型計(jì)算獲得的更精細(xì)的邊坡潛在滑移線(面)。

1.4 潛在滑動帶內(nèi)安全系數(shù)的計(jì)算

獲得邊坡的潛在滑動線(面)后，最簡單的方法就是利用傳統(tǒng)的極限平衡法求解獲得邊坡的穩(wěn)定安全系數(shù)。

實(shí)際上，從概率統(tǒng)計(jì)的角度來看，可以認(rèn)為，邊坡離散體系中各接觸點(diǎn)的極限滑移比與邊坡的整體安全性是相關(guān)的。首先定義各接觸點(diǎn)的極限滑移比的倒數(shù)為該接觸點(diǎn)的安全系數(shù)Ks

Ks=1/K

(9)

顯然Ks越大則該接觸點(diǎn)產(chǎn)生滑移的概率越低，反之亦然。按照離散體系中約定的圓盤排列方式，一個(gè)圓盤最多可能和周邊圓盤都接觸產(chǎn)生六個(gè)接觸點(diǎn)，如圖3(a)所示。這些接觸點(diǎn)中并不是每個(gè)接觸點(diǎn)都被納入計(jì)算滑動帶內(nèi)安全系數(shù)的統(tǒng)計(jì)樣本中，因?yàn)橛行┙佑|點(diǎn)的安全系數(shù)較大，將會導(dǎo)致邊坡的安全系數(shù)計(jì)算值偏大，這是偏于不安全的。為此，以第二次計(jì)算得到的最危險(xiǎn)滑動線(面)相關(guān)的圓盤為中心，再加上與之相接觸的兩層圓盤為整體考慮的區(qū)域，將這個(gè)區(qū)域內(nèi)圓盤接觸點(diǎn)作為一個(gè)初步統(tǒng)計(jì)樣本，但是要去除區(qū)域內(nèi)圓盤與區(qū)域外圓盤之間的接觸點(diǎn)。這些接觸點(diǎn)的安全系數(shù)可以構(gòu)成一個(gè)集合，計(jì)算其均值、方差、中位數(shù)值、最大值、最小值等。采用那些大于中值與中位數(shù)值，且與最大值相比大于三分之二的圓盤接觸點(diǎn)作為最終需要在計(jì)算邊坡安全系數(shù)中考慮的集合，該集合應(yīng)該是滑動帶內(nèi)圓盤接觸點(diǎn)的一個(gè)子集合，將這個(gè)集合定義為P。

對于集合P中的所有接觸點(diǎn)的安全系數(shù)，可以計(jì)算出它們的平均值K′以及方差值σ。這些值符合正態(tài)分布，可以將其規(guī)則化為一個(gè)標(biāo)準(zhǔn)的正態(tài)分布，這些安全系數(shù)點(diǎn)的分布可以計(jì)算出安全系數(shù)。

2 關(guān)鍵問題的討論

2.1 建模中剛性塊體的設(shè)定

按照傳統(tǒng)的極限平衡理論分析，在邊坡達(dá)到極限平衡狀態(tài)即將產(chǎn)生滑移失穩(wěn)的瞬間，邊坡在滑動面區(qū)域以外的應(yīng)力分布并沒有同時(shí)到達(dá)極限狀態(tài)，其形狀的整體性以及其它土力學(xué)性質(zhì)基本是保持不變的。應(yīng)力達(dá)到極限狀態(tài)的區(qū)域大致在滑動面附近一個(gè)不大的區(qū)域內(nèi)，因此將其滑動帶以外大部分區(qū)域按照極限平衡法的思想設(shè)定為若干個(gè)剛性塊體是基本符合邊坡的土力學(xué)特征的，這樣可以縮小圓盤所需覆蓋的范圍，簡化進(jìn)一步的計(jì)算。

從土力學(xué)原理分析，在潛在滑移帶以外的區(qū)域有一部分是位于滑動面曲線的上部的(滑移帶的凹向)，該部分在發(fā)生邊坡失穩(wěn)時(shí)有一定的整體性，因此可以整體上將這一個(gè)部分設(shè)置為剛性塊體，這樣在第二次建模計(jì)算中可大大減少圓盤單元數(shù)量，加快運(yùn)算速度。但是如果采用三角形的剛性塊體，在靜力學(xué)分析中，其尖角受力分析比較困難，在文獻(xiàn)[14]中特別設(shè)置了性質(zhì)復(fù)雜的力學(xué)模型，這樣無疑會加大計(jì)算量。為此，本文建立的是一個(gè)梯形剛體，計(jì)算分析表明是一個(gè)比較合理的方案。

2.2 圓盤直徑選擇及其與邊界接觸的處理

按照PFC2D的試驗(yàn)分析，圓盤的直徑越小，得到的應(yīng)力分析結(jié)果越準(zhǔn)確，但計(jì)算量也會呈指數(shù)級增加。而且，圓盤直徑減小到一定尺度后，進(jìn)一步的減小不會實(shí)質(zhì)性的提高計(jì)算準(zhǔn)確度[15]。顯然，模型中的圓盤直徑越大，覆蓋區(qū)域所需的圓盤越少，計(jì)算速度越快，但是計(jì)算精度不高；反之亦然。在分析實(shí)驗(yàn)結(jié)果的基礎(chǔ)上，第一次計(jì)算時(shí)，模型中采用邊坡高度值的2%至5%作為圓盤直徑可以保證后續(xù)計(jì)算結(jié)果的精確度。

由于邊坡離散化方案中采用的是三角形穩(wěn)定覆蓋，因此在圓盤覆蓋的過程中，底部的圓盤可以與底部的水平邊界完全接觸；但是對于垂直邊界及傾斜剛體邊界，若采用直線邊界，則難以做到與邊界上的所有圓盤相接觸，此時(shí)圓盤受力分析可能會產(chǎn)生較大的誤差。傳統(tǒng)的處理方法是用更小的圓盤進(jìn)行填充，使邊界盡可能與邊界圓盤接觸，這樣做會增加計(jì)算量，且不符合本文的離散方案。為此本文考慮采用折線邊界來代替直線邊界，或者用矩形小剛塊來契合邊界和鄰近的圓盤(可以簡單的將這些小矩形剛性塊體理解為剛性邊界的突出部分)。這樣處理，可以在不增加計(jì)算量的情況下滿足計(jì)算精度要求。

3 計(jì)算程序與流程

上述離散元極限平衡法計(jì)算邊坡穩(wěn)定性的方法已基于MATLAB平臺編制了相應(yīng)的邊坡潛在滑動帶安全系數(shù)計(jì)算程序Slope-zone，其基本計(jì)算流程如圖7所示。

圖7離散元極限平衡法計(jì)算邊坡穩(wěn)定性的基本流程

4 算例分析

4.1 均質(zhì)土坡穩(wěn)定性計(jì)算

本算例為澳大利亞邊坡委員會1986年給出的一個(gè)經(jīng)典算例。某均質(zhì)土層邊坡，邊坡高度10 m，上坡長度20 m，下坡長度40 m，土重度20 kN/m3，內(nèi)摩擦角19.6°，黏聚力為3 kPa，建議的安全系數(shù)標(biāo)準(zhǔn)答案為1.00，也可以按照式(10)對安全系數(shù)進(jìn)行預(yù)估。

(10)

式中：φ為內(nèi)摩擦角；c為黏聚力；h為邊坡高度；γ為土的重度；β為邊坡的傾斜角度。

采用本文方法進(jìn)行分析(第一次采用了3 256個(gè)圓盤，圓盤直徑為邊坡高度的4%，第二次采用4 230個(gè)圓盤，圓盤直徑為邊坡高度的2.5%，兩塊剛體區(qū)域)，并與畢肖普法、式(10)、離散元法結(jié)果進(jìn)行比較，結(jié)果如表1所示。

表1 均質(zhì)土坡算例計(jì)算結(jié)果

由表1的計(jì)算結(jié)果可以看出，本文方法的計(jì)算精度優(yōu)于離散元法，且計(jì)算時(shí)間遠(yuǎn)小于離散元法。

4.2 水平層狀邊坡穩(wěn)定性計(jì)算

某水平層狀土質(zhì)邊坡[16-18]，由上下兩層土體組成，坡高24 m，坡角37°，兩層土體的厚度分別為18 m和6 m，內(nèi)摩擦角分別為30°和20°，黏聚力分別為7 kPa和2 kPa，如圖8所示。

圖8分層邊坡示意圖

利用極限平衡法中的畢肖普法按照0.618法搜索潛在滑動圓弧得到的最小安全系數(shù)為1.424，利用PFC2D的計(jì)算結(jié)果為1.472，用本文本文方法計(jì)算得到的結(jié)果為1.441。其中，第一次采用了3 565個(gè)圓盤，圓盤直徑為邊坡高度的4%；第二次采用4 623個(gè)圓盤，圓盤直徑為邊坡高度的2.5%。計(jì)算結(jié)果如表2所示。

表2 非均質(zhì)土坡算例計(jì)算結(jié)果

5 結(jié) 語

本文提出了一種結(jié)合離散元法和極限平衡分析法各自優(yōu)點(diǎn)的離散元極限平衡法，給出了分析邊坡穩(wěn)定性的步驟，并基于MATLAB平臺編制了相應(yīng)的計(jì)算程序。采用本文方法分別針對均質(zhì)土坡和水平分層邊坡進(jìn)行了穩(wěn)定性分析，并與常用的極限平衡法和離散元法計(jì)算結(jié)果進(jìn)行了比較，結(jié)果表明本文方法不僅能保證較高的精度，而且比離散元法效率更高。

(1) 結(jié)合離散元法和極限平衡法的核心思想，提出了一種用于邊坡穩(wěn)定性分析的離散元極限平衡方法，將搜索邊坡最不利滑動面的復(fù)雜工作轉(zhuǎn)化為離散單元間相對摩擦系數(shù)的排序，計(jì)算復(fù)雜度大大降低，且能保證較高的精度。

(2) 本文方法的核心思想之一是將邊坡最危險(xiǎn)滑移面(線)的搜索分兩步計(jì)算完成。第一步是初步計(jì)算，得到潛在的滑移帶區(qū)域；第二步是精細(xì)計(jì)算，將第一次計(jì)算確定的潛在滑移帶區(qū)域以外的部分以剛性塊體替代，而在潛在滑移帶區(qū)域以更小的圓盤單元進(jìn)行覆蓋，從而減少圓盤單元數(shù)目，提高計(jì)算效率。

(3) 針對本文提出的離散元極限平衡法，基于MATLAB平臺編制了相應(yīng)的計(jì)算程序，可用于簡單邊坡和復(fù)雜邊坡的穩(wěn)定性分析。針對均質(zhì)土坡和水平層狀邊坡的算例結(jié)果表明，相較于離散元法，本文方法不僅具有更高的計(jì)算效率，而且具有更高的計(jì)算精度。

參考文獻(xiàn)：

[1] 陳祖煜.土質(zhì)邊坡穩(wěn)定分析——原理·方法·程序[M].北京:中國水利水電出版社,2003.

[2] 陳祖煜.巖質(zhì)邊坡穩(wěn)定分析——原理·方法·程序[M].北京:中國水利水電出版社,2003.

[3] 鄭穎人,趙尚毅,李安紅,等.有限元極限分析法及其在邊坡中的應(yīng)用[M].北京:人民交通出版社,2011.

[4] 馬曉雨.基于強(qiáng)度折減有限元法的土坡穩(wěn)定性分析研究[D].大連：大連理工大學(xué),2008.

[5] 劉漢東，賈聿頡.基于FLAC3D的邊坡穩(wěn)定性分析自編強(qiáng)度折減程序的修正[J].華北水利水電大學(xué)學(xué)報(bào)(自然科學(xué)版)，2015,36(6):59-62.

[6] 李典慶，肖 特，曹子君，等.基于極限平衡法和有限元法的邊坡協(xié)同式可靠度分析[J].巖土工程學(xué)報(bào)，2016,38(6):1004-1012.

[7] 王泳嘉，邢紀(jì)波.離散單元法及其在巖土力學(xué)中的應(yīng)用[M].沈陽:東北工學(xué)院出版社,1991.

[8] 常新正.基于GPU的顆粒離散元計(jì)算方法研究[D].大連：大連理工大學(xué)，2013.

[9] 田 健.基于OpenCL的離散元方法仿真優(yōu)化軟件改進(jìn)研究[D].吉林：吉林大學(xué)，2015.

[10] Itasca Consulting Group, Inc. PFC2D(Particle code in 2 Dimensions), Version 2.0[M]. Minneapolis:ICG, 2002.

[11] 王秀菊，石 崇，李德杰，等.基于距離和局部Delaunay三角化控制的顆粒離散元模型填充方法研究[J].巖土力學(xué)，2015，36(5)：2081-2087.

[12] 孫小三.邊坡穩(wěn)定性的條分法和有限元法耦合分析[D].杭州：浙江大學(xué),2005.

[13] 羅 勇.土工問題的顆粒流模擬及應(yīng)用研究[D].杭州：浙江大學(xué),2007.

[14] 秦忠國，張向陽.一種邊坡穩(wěn)定分析的新方法[J].巖土工程學(xué)報(bào),2016,38(5)946-951.

[15] 邵 帥.土石混合體邊坡穩(wěn)定性的FEM-DEM分析[D].大連：大連理工大學(xué),2012.

[16] 趙尚毅.用有限元強(qiáng)度折減法求邊坡穩(wěn)定安全系數(shù)[J].巖土工程學(xué)報(bào),2002,24(3):343-346.

[17] 趙尚毅，鄭穎人，張玉芳.極限分析有限元法講座——Ⅱ有限元強(qiáng)度折減法中邊坡失穩(wěn)的判據(jù)探討[J].巖土力學(xué),2005,26(2):332-336.

[18] 欒茂田，武亞軍，年延凱.強(qiáng)度折減有限元法中邊坡失穩(wěn)的塑性區(qū)判據(jù)及其應(yīng)用[J].防災(zāi)減災(zāi)工程學(xué)報(bào),2003,23(3):1-8.