基于矢量和法的三維邊坡滑面搜索研究

鄒 燁，嚴(yán)東方

(中交第二公路勘察設(shè)計(jì)研究院有限公司， 湖北 武漢 430056)

臨界滑動(dòng)面的確定是邊坡穩(wěn)定性分析過(guò)程中的基礎(chǔ)問(wèn)題，目前大部分邊坡穩(wěn)定性分析方法都是以某一滑體或滑面為對(duì)象進(jìn)行的。定義邊坡臨界滑面的方法有多種，邵國(guó)建等[1]研究了利用干擾能量等值面圖確定臨界滑動(dòng)面。孫冠華等[2]提出臨界滑面上的點(diǎn)為沿深度方向上等效塑性應(yīng)變的極大值點(diǎn)。而目前在工程中應(yīng)用最廣泛的是基于最小安全系數(shù)的臨界滑面定義方式。

自然界中的滑坡多呈現(xiàn)三維狀態(tài)，三維極限平衡分析方法直接從二維情況下發(fā)展而來(lái)，在這一過(guò)程中為求解方程組引入的假設(shè)削弱了方法的理論與適用性[3]。陳祖煜等[4-5]建立了三維極限平衡法的上、下限分析體系。鄭宏[6]提出了能適應(yīng)任意形狀滑動(dòng)面并能滿(mǎn)足所有平衡條件的嚴(yán)格三維極限平衡法。矢量和法最早由葛修潤(rùn)院士于1983年提出，該方法抓住力是矢量這一基本屬性，在極限平衡分析與數(shù)值分析中均有應(yīng)用[7-9]，且能很方便的擴(kuò)展到三維情形。

目前利用矢量和法進(jìn)行邊坡的穩(wěn)定性分析時(shí)，或結(jié)合地質(zhì)調(diào)查情況與研究者經(jīng)驗(yàn)確定潛滑面，或參照極限平衡分析法的臨界滑面，而并沒(méi)有通過(guò)滑面搜索確定矢量和法的臨界滑面。對(duì)此，本文基于邊坡應(yīng)力場(chǎng)，研究利用矢量和法進(jìn)行三維邊坡臨界滑面的搜索問(wèn)題。

1 三維矢量和法理論

矢量和法用滑面上滑動(dòng)力與抗滑力的矢量和來(lái)求解安全系數(shù)，求解圖示見(jiàn)圖1。

其安全系數(shù)定義為：

(1)

圖1三維矢量和法計(jì)算圖示

其中：R為滑面上的抗滑力矢量和在滑動(dòng)趨勢(shì)方向上的投影；T為滑面上的下滑力矢量和在滑動(dòng)趨勢(shì)方向上的投影。R與T的計(jì)算見(jiàn)式(2)、式(3)：

(2)

(3)

式中：dS為滑面上的面元。

σr與σt為滑面上任一點(diǎn)的抗滑力與下滑力，分別按式(4)、式(5)進(jìn)行計(jì)算：

σr=σs1+σn1

(4)

σt=σ·n

(5)

σt同時(shí)可以表示為：

σt=σs2+σn2

(6)

式中：σn1、σs1為該點(diǎn)的抗滑正應(yīng)力與極限抗滑剪應(yīng)力；σn2、σs2為該點(diǎn)的下滑正應(yīng)力與下滑切應(yīng)力；σ為該點(diǎn)的應(yīng)力張量；n為滑面上該點(diǎn)的外法向量(以指向滑體內(nèi)部為正)。σn2可由σt表示為：

σn2=(σt·n)n

(7)

σn1與σn2為一對(duì)作用力與反作用力，兩者大小相等，方向相反，即：

σn1=-σn2

(8)

極限抗滑剪應(yīng)力σs1的計(jì)算公式為：

σs1=(c-σn1·tanφ)·t

(9)

式中：σn1為該點(diǎn)正應(yīng)力的大小，并采用拉正壓負(fù)的假定；t為該點(diǎn)極限抗滑剪應(yīng)力的方向向量。

郭明偉等[8]對(duì)三維情形下的矢量和法極限抗滑剪應(yīng)力方向t進(jìn)行了研究。t的計(jì)算方式表述為：初始滑動(dòng)趨勢(shì)方向D0在該點(diǎn)切平面上投影方向的反方向。其中，初始滑動(dòng)趨勢(shì)方向D0定義為各dS面上實(shí)際剪應(yīng)力的合力矢方向，即：

(10)

D為滑動(dòng)趨勢(shì)方向，其定義為：滑面上極限抗滑剪應(yīng)力方向的反方向，即：

(11)

2 搜索算法

臨界滑面的確定本質(zhì)上是一個(gè)區(qū)域極值的求解問(wèn)題，目前很多研究工作集中在搜索算法的選擇上，常用搜索算法有遺傳算法[10-11]、螞蟻算法[12]、神經(jīng)網(wǎng)絡(luò)算法[13]、模擬退火法[14]等。本文采用的單純型調(diào)優(yōu)法是一種局部化的搜索方法，具有簡(jiǎn)單高效的特點(diǎn)[15]。

以安全系數(shù)Fs為目標(biāo)函數(shù)，對(duì)于確定的邊坡應(yīng)力場(chǎng)，影響安全系數(shù)大小的因素就是滑面的位置。

初始滑動(dòng)面一般通過(guò)一系列綜合分析即可得到。事實(shí)上，對(duì)實(shí)際工程來(lái)說(shuō)，利用勘察的方法往往能夠確定潛在滑動(dòng)面的大致位置。

對(duì)初始滑面X0：

X0={x1，x2，…xn}T

(12)

式中：x1，x2，xn可以為滑面上點(diǎn)的坐標(biāo)；滑面形狀已知的情況下，也可以是滑面方程中的參數(shù)。

利用初始單純型X0構(gòu)造初始單純型S0：

S0={X1，X2，…X1…Xn}

(13)

其中，

(14)

(15)

(16)

式中：t為單純型的步距，可以根據(jù)實(shí)際情況進(jìn)行調(diào)整。

在初始單純型S0的基礎(chǔ)上，經(jīng)過(guò)反射、擴(kuò)張、壓縮、收縮等手段，讓單純型進(jìn)行翻滾、變形并逐步地向目標(biāo)點(diǎn)靠攏，直到目標(biāo)函數(shù)值(即安全系數(shù)Fs)滿(mǎn)足停機(jī)準(zhǔn)則方可停止。

本文的停機(jī)準(zhǔn)則為：除去令安全系數(shù)取最大值的滑面Xh(即單純型中令目標(biāo)函數(shù)取最大值的點(diǎn))，單純型中其余各頂點(diǎn)的目標(biāo)函數(shù)值(安全系數(shù))的均方差小于等于1×10-4(該數(shù)可以根據(jù)需要進(jìn)行選取，也可選其他合適的數(shù))時(shí)，停止搜索，見(jiàn)式(17)。同時(shí)，取此時(shí)單純型中目標(biāo)函數(shù)的最小值為最小安全系數(shù)，對(duì)應(yīng)的滑面即為臨界滑面。

(17)

(18)

本文的基本搜索流程見(jiàn)圖2。

圖2搜索流程圖

本文利用單純型法，從二維邊坡入手，利用兩個(gè)算例分別進(jìn)行二維與三維邊坡的滑面搜索研究。

3 算 例

3.1 算例1-二維算例

許多邊坡穩(wěn)定性分析方法都用澳大利亞計(jì)算機(jī)應(yīng)用協(xié)會(huì)(ACADS)的邊坡考題做過(guò)驗(yàn)證。該例取自ACADS第一個(gè)考題EX1(a)。EX1(a)是一個(gè)均質(zhì)邊坡，數(shù)值計(jì)算所用的材料參數(shù)見(jiàn)表1。

表1 考題EX1(a)材料參數(shù)

首先建立模型：在有限元軟件ABAQUS中建立邊坡計(jì)算模型，采用ABAQUS中的平面應(yīng)變單元(CPE4)，單元數(shù)目為5 459個(gè)，邊坡材料服從摩爾-庫(kù)侖強(qiáng)度準(zhǔn)則?？碱}EX1(a)的有限元計(jì)算模型見(jiàn)圖3。

圖3算例1有限元模型

利用SLIDE軟件，可得利用極限平衡法求解得到的臨界滑動(dòng)面位置，并以此滑動(dòng)面為初始滑動(dòng)面進(jìn)行滑面的搜索。

本文使用多項(xiàng)式擬合二維滑面曲線，為保證精度，可采用了盡可能多的點(diǎn)。本文取初始滑面上的12個(gè)點(diǎn)來(lái)構(gòu)造初始單純型。

初始單純型即為：

X0={x1，y1，x2，y2，…x12，y12}T

(19)

t取0.5，根據(jù)式(13)—式(16)即可求得初始單純型S0。為節(jié)省篇幅，本文不再給出S0的具體計(jì)算過(guò)程。

滑面搜索過(guò)程見(jiàn)圖4與圖5。

圖4 滑面搜索過(guò)程

圖5二維滑面搜索過(guò)程收斂曲線

矢量和法與其他三種常用極限平衡條分法計(jì)算所得的臨界滑面位置對(duì)比見(jiàn)圖6。

圖6臨界滑面位置對(duì)比

安全系數(shù)的計(jì)算結(jié)果可與澳大利亞巖土工程協(xié)會(huì)公布的裁判答案進(jìn)行對(duì)比，見(jiàn)表2。

表2 二維邊坡搜索結(jié)果對(duì)比分析

需指出的是，表2中滑動(dòng)趨勢(shì)方向θ為與水平方向的夾角，而極限平衡法是邊坡穩(wěn)定性分析軟件SLIDE進(jìn)行的；筆者認(rèn)為，該軟件在計(jì)算時(shí)沒(méi)有考慮材料的彈性，利用極限平衡法求安全系數(shù)時(shí)，可能會(huì)使結(jié)果偏小。

3.2 算例2-三維算例

在二維分析的基礎(chǔ)上，采用相同的思路進(jìn)行三維邊坡滑動(dòng)面的搜索。采用Zhang[16]的三維邊坡算例，對(duì)自編的三維矢量和法程序進(jìn)行驗(yàn)證。有許多研究人員曾利用該算例來(lái)驗(yàn)證各自的邊坡計(jì)算方案，因而這里使用該算例也具有更高的可信度。

該邊坡的坡比為1∶5，坡高為12.2 m，滑動(dòng)面的形狀為一個(gè)橢球面，其方程為：

(20)

邊坡的尺寸以及滑動(dòng)面在邊坡中的位置見(jiàn)圖7；邊坡模型在Y軸方向的尺寸為200 m，滑面位于邊坡的中間位置。

圖7三維邊坡外形輪廓(單位：m)

在ABAQUS中建立有限元模型，見(jiàn)圖8。三維有限元模型采用的是八節(jié)點(diǎn)六面體單元(C3D8)，共劃分了71 190個(gè)單元。

圖8三維邊坡有限元模型

該三維邊坡為均質(zhì)邊坡，數(shù)值計(jì)算所用材料參數(shù)見(jiàn)表3。利用ABAQUS軟件進(jìn)行有限元分析得到邊坡的應(yīng)力場(chǎng)。

表3 三維邊坡算例材料參數(shù)

采用單純型優(yōu)化法搜索滑面并進(jìn)行安全系數(shù)的計(jì)算，滑面搜索仍然以橢球面進(jìn)行。

單純型法的優(yōu)化變量為橢球面的球心坐標(biāo)以及各軸的軸長(zhǎng)。考慮到對(duì)稱(chēng)性，球心Y坐標(biāo)不變，初始滑面的確定共使用5個(gè)變量，即：

X0={36.6，27.4，24.4，24.4，66.9}T

(21)

t取1.0，根據(jù)式(13)—式(16)即可求得初始單純型S0。

滑面的搜索過(guò)程見(jiàn)圖9與圖10。

圖9 三維滑面的搜索過(guò)程

圖10三維滑面搜索過(guò)程收斂曲線

將矢量和法的臨界滑面與初始滑面以及文獻(xiàn)[17]中搜索得到的滑面進(jìn)行對(duì)比，見(jiàn)圖11。

圖11三種方法臨界滑面的對(duì)比

其中，林永生等[17]利用遺傳算法進(jìn)行滑面搜索，在搜索過(guò)程中，只改變橢球體球心的X坐標(biāo)與Z坐標(biāo)，而軸長(zhǎng)均不發(fā)生改變，也就是說(shuō)，橢球體的大小及形狀均不發(fā)生改變，而只改變橢球體的位置，最終所得到的滑面方程為：

(22)

本文的搜索只對(duì)滑面的形狀做了限定，即要求滑面為一個(gè)橢球體，對(duì)軸長(zhǎng)與橢球的球心位置沒(méi)有限制。通常限定條件越少，其結(jié)果的與實(shí)際也會(huì)越接近。通過(guò)搜索最終得到的滑面方程為：

(23)

比較式(22)、式(23)，可以看出，在限定滑面為橢球面的條件下，三維矢量和法的臨界滑動(dòng)面與Zhang[16]的計(jì)算結(jié)果有一定的區(qū)別，其中矢量和法與Zhang的滑面在坡頂部分比較接近，而坡底區(qū)域，矢量和法的滑面更加的靠近坡腳。

將矢量和法計(jì)算得到的安全系數(shù)與各文獻(xiàn)所得結(jié)果進(jìn)行對(duì)比，所得結(jié)果見(jiàn)表4。從該表可以看出：當(dāng)矢量和法采用與式(20)的滑面進(jìn)行計(jì)算時(shí)，其安全系數(shù)偏差為3.86%。而經(jīng)過(guò)滑面搜索，其安全系數(shù)略有降低。值得指出的是，表4中關(guān)于誤差的計(jì)算是以Zhang的計(jì)算結(jié)果為基準(zhǔn)進(jìn)行的，但并不代表其為標(biāo)準(zhǔn)答案。

表4 三維邊坡搜索結(jié)果對(duì)比分析

事實(shí)上，鄭宏[6]指出，Zhang的方法忽略了三維條塊上四棱柱沿滑動(dòng)方向上兩個(gè)側(cè)面上的剪應(yīng)力可能導(dǎo)致其計(jì)算結(jié)果偏小?？梢哉J(rèn)為，利用矢量和法通過(guò)滑面搜索得到的計(jì)算結(jié)果與實(shí)際情況更為接近。

4 結(jié) 論

本文研究了三維矢量和法臨界滑動(dòng)面的搜索問(wèn)題，給出了滑動(dòng)面搜索的全過(guò)程，得到以下幾點(diǎn)結(jié)論：

(1) 矢量和法能很方便的進(jìn)行二維與三維邊坡的穩(wěn)定性分析，在邊坡應(yīng)力場(chǎng)求解準(zhǔn)確的情況下，該方法所得結(jié)果具有較好的可靠性。

(2) 二維情形下，利用矢量和法求解極限平衡分析臨界滑面的結(jié)果與矢量和法臨界滑面的結(jié)果相近；三維情形下，進(jìn)行滑面搜索后能更加準(zhǔn)確的進(jìn)行邊坡的穩(wěn)定性分析。