鄒 燁,嚴東方
(中交第二公路勘察設計研究院有限公司, 湖北 武漢 430056)
臨界滑動面的確定是邊坡穩定性分析過程中的基礎問題,目前大部分邊坡穩定性分析方法都是以某一滑體或滑面為對象進行的。定義邊坡臨界滑面的方法有多種,邵國建等[1]研究了利用干擾能量等值面圖確定臨界滑動面。孫冠華等[2]提出臨界滑面上的點為沿深度方向上等效塑性應變的極大值點。而目前在工程中應用最廣泛的是基于最小安全系數的臨界滑面定義方式。
自然界中的滑坡多呈現三維狀態,三維極限平衡分析方法直接從二維情況下發展而來,在這一過程中為求解方程組引入的假設削弱了方法的理論與適用性[3]。陳祖煜等[4-5]建立了三維極限平衡法的上、下限分析體系。鄭宏[6]提出了能適應任意形狀滑動面并能滿足所有平衡條件的嚴格三維極限平衡法。矢量和法最早由葛修潤院士于1983年提出,該方法抓住力是矢量這一基本屬性,在極限平衡分析與數值分析中均有應用[7-9],且能很方便的擴展到三維情形。
目前利用矢量和法進行邊坡的穩定性分析時,或結合地質調查情況與研究者經驗確定潛滑面,或參照極限平衡分析法的臨界滑面,而并沒有通過滑面搜索確定矢量和法的臨界滑面。對此,本文基于邊坡應力場,研究利用矢量和法進行三維邊坡臨界滑面的搜索問題。
矢量和法用滑面上滑動力與抗滑力的矢量和來求解安全系數,求解圖示見圖1。
其安全系數定義為:
(1)
圖1三維矢量和法計算圖示
其中:R為滑面上的抗滑力矢量和在滑動趨勢方向上的投影;T為滑面上的下滑力矢量和在滑動趨勢方向上的投影。R與T的計算見式(2)、式(3):
(2)
(3)
式中:dS為滑面上的面元。
σr與σt為滑面上任一點的抗滑力與下滑力,分別按式(4)、式(5)進行計算:
σr=σs1+σn1
(4)
σt=σ·n
(5)
σt同時可以表示為:
σt=σs2+σn2
(6)
式中:σn1、σs1為該點的抗滑正應力與極限抗滑剪應力;σn2、σs2為該點的下滑正應力與下滑切應力;σ為該點的應力張量;n為滑面上該點的外法向量(以指向滑體內部為正)。σn2可由σt表示為:
σn2=(σt·n)n
(7)
σn1與σn2為一對作用力與反作用力,兩者大小相等,方向相反,即:
σn1=-σn2
(8)
極限抗滑剪應力σs1的計算公式為:
σs1=(c-σn1·tanφ)·t
(9)
式中:σn1為該點正應力的大小,并采用拉正壓負的假定;t為該點極限抗滑剪應力的方向向量。
郭明偉等[8]對三維情形下的矢量和法極限抗滑剪應力方向t進行了研究。t的計算方式表述為:初始滑動趨勢方向D0在該點切平面上投影方向的反方向。其中,初始滑動趨勢方向D0定義為各dS面上實際剪應力的合力矢方向,即:
(10)
D為滑動趨勢方向,其定義為:滑面上極限抗滑剪應力方向的反方向,即:
(11)
臨界滑面的確定本質上是一個區域極值的求解問題,目前很多研究工作集中在搜索算法的選擇上,常用搜索算法有遺傳算法[10-11]、螞蟻算法[12]、神經網絡算法[13]、模擬退火法[14]等。本文采用的單純型調優法是一種局部化的搜索方法,具有簡單高效的特點[15]。
以安全系數Fs為目標函數,對于確定的邊坡應力場,影響安全系數大小的因素就是滑面的位置。
初始滑動面一般通過一系列綜合分析即可得到。事實上,對實際工程來說,利用勘察的方法往往能夠確定潛在滑動面的大致位置。
對初始滑面X0:
X0={x1,x2,…xn}T
(12)
式中:x1,x2,xn可以為滑面上點的坐標;滑面形狀已知的情況下,也可以是滑面方程中的參數。
利用初始單純型X0構造初始單純型S0:
S0={X1,X2,…X1…Xn}
(13)
其中,
(14)
(15)
(16)
式中:t為單純型的步距,可以根據實際情況進行調整。
在初始單純型S0的基礎上,經過反射、擴張、壓縮、收縮等手段,讓單純型進行翻滾、變形并逐步地向目標點靠攏,直到目標函數值(即安全系數Fs)滿足停機準則方可停止。
本文的停機準則為:除去令安全系數取最大值的滑面Xh(即單純型中令目標函數取最大值的點),單純型中其余各頂點的目標函數值(安全系數)的均方差小于等于1×10-4(該數可以根據需要進行選取,也可選其他合適的數)時,停止搜索,見式(17)。同時,取此時單純型中目標函數的最小值為最小安全系數,對應的滑面即為臨界滑面。
(17)
(18)
本文的基本搜索流程見圖2。
圖2搜索流程圖
本文利用單純型法,從二維邊坡入手,利用兩個算例分別進行二維與三維邊坡的滑面搜索研究。
許多邊坡穩定性分析方法都用澳大利亞計算機應用協會(ACADS)的邊坡考題做過驗證。該例取自ACADS第一個考題EX1(a)。EX1(a)是一個均質邊坡,數值計算所用的材料參數見表1。
表1 考題EX1(a)材料參數
首先建立模型:在有限元軟件ABAQUS中建立邊坡計算模型,采用ABAQUS中的平面應變單元(CPE4),單元數目為5 459個,邊坡材料服從摩爾-庫侖強度準則??碱}EX1(a)的有限元計算模型見圖3。
圖3算例1有限元模型
利用SLIDE軟件,可得利用極限平衡法求解得到的臨界滑動面位置,并以此滑動面為初始滑動面進行滑面的搜索。
本文使用多項式擬合二維滑面曲線,為保證精度,可采用了盡可能多的點。本文取初始滑面上的12個點來構造初始單純型。
初始單純型即為:
X0={x1,y1,x2,y2,…x12,y12}T
(19)
t取0.5,根據式(13)—式(16)即可求得初始單純型S0。為節省篇幅,本文不再給出S0的具體計算過程。
滑面搜索過程見圖4與圖5。
圖4 滑面搜索過程
圖5二維滑面搜索過程收斂曲線
矢量和法與其他三種常用極限平衡條分法計算所得的臨界滑面位置對比見圖6。
圖6臨界滑面位置對比
安全系數的計算結果可與澳大利亞巖土工程協會公布的裁判答案進行對比,見表2。
表2 二維邊坡搜索結果對比分析
需指出的是,表2中滑動趨勢方向θ為與水平方向的夾角,而極限平衡法是邊坡穩定性分析軟件SLIDE進行的;筆者認為,該軟件在計算時沒有考慮材料的彈性,利用極限平衡法求安全系數時,可能會使結果偏小。
在二維分析的基礎上,采用相同的思路進行三維邊坡滑動面的搜索。采用Zhang[16]的三維邊坡算例,對自編的三維矢量和法程序進行驗證。有許多研究人員曾利用該算例來驗證各自的邊坡計算方案,因而這里使用該算例也具有更高的可信度。
該邊坡的坡比為1∶5,坡高為12.2 m,滑動面的形狀為一個橢球面,其方程為:
(20)
邊坡的尺寸以及滑動面在邊坡中的位置見圖7;邊坡模型在Y軸方向的尺寸為200 m,滑面位于邊坡的中間位置。
圖7三維邊坡外形輪廓(單位:m)
在ABAQUS中建立有限元模型,見圖8。三維有限元模型采用的是八節點六面體單元(C3D8),共劃分了71 190個單元。
圖8三維邊坡有限元模型
該三維邊坡為均質邊坡,數值計算所用材料參數見表3。利用ABAQUS軟件進行有限元分析得到邊坡的應力場。
表3 三維邊坡算例材料參數
采用單純型優化法搜索滑面并進行安全系數的計算,滑面搜索仍然以橢球面進行。
單純型法的優化變量為橢球面的球心坐標以及各軸的軸長??紤]到對稱性,球心Y坐標不變,初始滑面的確定共使用5個變量,即:
X0={36.6,27.4,24.4,24.4,66.9}T
(21)
t取1.0,根據式(13)—式(16)即可求得初始單純型S0。
滑面的搜索過程見圖9與圖10。
圖9 三維滑面的搜索過程
圖10三維滑面搜索過程收斂曲線
將矢量和法的臨界滑面與初始滑面以及文獻[17]中搜索得到的滑面進行對比,見圖11。
圖11三種方法臨界滑面的對比
其中,林永生等[17]利用遺傳算法進行滑面搜索,在搜索過程中,只改變橢球體球心的X坐標與Z坐標,而軸長均不發生改變,也就是說,橢球體的大小及形狀均不發生改變,而只改變橢球體的位置,最終所得到的滑面方程為:
(22)
本文的搜索只對滑面的形狀做了限定,即要求滑面為一個橢球體,對軸長與橢球的球心位置沒有限制。通常限定條件越少,其結果的與實際也會越接近。通過搜索最終得到的滑面方程為:
(23)
比較式(22)、式(23),可以看出,在限定滑面為橢球面的條件下,三維矢量和法的臨界滑動面與Zhang[16]的計算結果有一定的區別,其中矢量和法與Zhang的滑面在坡頂部分比較接近,而坡底區域,矢量和法的滑面更加的靠近坡腳。
將矢量和法計算得到的安全系數與各文獻所得結果進行對比,所得結果見表4。從該表可以看出:當矢量和法采用與式(20)的滑面進行計算時,其安全系數偏差為3.86%。而經過滑面搜索,其安全系數略有降低。值得指出的是,表4中關于誤差的計算是以Zhang的計算結果為基準進行的,但并不代表其為標準答案。
表4 三維邊坡搜索結果對比分析
事實上,鄭宏[6]指出,Zhang的方法忽略了三維條塊上四棱柱沿滑動方向上兩個側面上的剪應力可能導致其計算結果偏小??梢哉J為,利用矢量和法通過滑面搜索得到的計算結果與實際情況更為接近。
本文研究了三維矢量和法臨界滑動面的搜索問題,給出了滑動面搜索的全過程,得到以下幾點結論:
(1) 矢量和法能很方便的進行二維與三維邊坡的穩定性分析,在邊坡應力場求解準確的情況下,該方法所得結果具有較好的可靠性。
(2) 二維情形下,利用矢量和法求解極限平衡分析臨界滑面的結果與矢量和法臨界滑面的結果相近;三維情形下,進行滑面搜索后能更加準確的進行邊坡的穩定性分析。