謝旭奎, 張 杰

(湖南致力工程科技有限公司，湖南 長沙 410208)

有限元強度折減法由Zienkiewice[1]于1975年提出，國內(nèi)外許多學(xué)者在這方面做了大量的工作[2-12]。目前對邊坡穩(wěn)定性的有限元分析大多都集中于二維模型，而二維模型分析的結(jié)果推廣到實際邊坡應(yīng)用時，往往存在較大的誤差，因此當(dāng)前多采用三維模型分析邊坡穩(wěn)定性。目前國內(nèi)外采用的邊坡有限元分析軟件主要為FlAC 3D、hyperworks、algor、nastran、abaqus等，雖然ANSYS軟件操作復(fù)雜，巖土方面運用較少，但ANSYS在模型建立和網(wǎng)格劃分優(yōu)勢明顯，因此本文采用ANSYS進(jìn)行邊坡穩(wěn)定性分析?；趶姸日蹨p法的邊坡臨界破壞狀態(tài)主要有以下傳統(tǒng)判據(jù)[7-10]:①邊坡破壞的一個重要依據(jù)是有限元計算是否收斂，當(dāng)計算收斂時，邊坡穩(wěn)定;計算不收斂，則邊坡破壞。②塑性區(qū)是否貫通是邊坡失穩(wěn)的另外一個標(biāo)志，在有限元計算中可以得到不同折減強度下的邊坡塑性應(yīng)變圖，當(dāng)塑性區(qū)貫通時，認(rèn)為邊坡破壞。③位移的急劇增加。在計算中當(dāng)發(fā)現(xiàn)邊坡的滑動面位移急劇增加時，認(rèn)為邊坡破壞。許多學(xué)者對坡體塑性區(qū)貫通狀態(tài)、有限元計算收斂性、位移值突變進(jìn)行研究，而對坡體的塑性應(yīng)力應(yīng)變值變化進(jìn)行較少研究，且對位移值突變研究僅限于定性的考慮位移突變?yōu)闊o窮大的可能性，但實際模擬中，隨著折減系數(shù)的增大，邊坡位移突變到某定值后會存在屈服下降，針對此種情況，本文提出一種新式突變判據(jù)，即在位移和塑性應(yīng)力應(yīng)變值兩類指標(biāo)上，采用新的曲線擬合方式，從理論上分析得出邊坡極限平衡狀態(tài)下的安全系數(shù)。另外對邊坡的等效塑性區(qū)變化進(jìn)行全過程模擬，在以往邊坡有限元分析中也比較少見。

1 有限元強度折減法

1.1 基本原理

有限元強度折減法的基本原理是通過不斷的調(diào)整巖土體的抗剪強度指標(biāo)，即將粘聚力和內(nèi)摩角、泊松比、彈性模量等用一個折減系數(shù)Ki進(jìn)行折減，得到折減后的新抗剪強度指標(biāo)，然后通過彈塑形有限元數(shù)值計算確定邊坡內(nèi)的應(yīng)力場、應(yīng)變場或位移場，并且基于對應(yīng)變或者位移的某些分布特征以及有限元中的某些數(shù)學(xué)特征進(jìn)行分析判斷，不斷增加折減系數(shù)Ki的值，直至邊坡失穩(wěn)，將此極限狀態(tài)下所對應(yīng)的折減系數(shù)定義為邊坡的穩(wěn)定安全系數(shù)。

1.2 參數(shù)調(diào)整

根據(jù)整體強度折減，內(nèi)摩擦角和黏聚力調(diào)整如下：

(1)

彈性參數(shù)(彈性模量和泊松比)調(diào)整強度折減時，

對于泊松比調(diào)整如下：

(2)

對于彈性模量調(diào)整如下：

(3)

1.3 屈服準(zhǔn)則

采用理想彈塑性本構(gòu)模型的Drucker-Prager屈服準(zhǔn)則。對于理想的彈塑性模型，土體進(jìn)入屈服階段就意味著破壞，將Mohr-Coulomb準(zhǔn)則作為屈服準(zhǔn)則變換得到Mohr-Coulomb屈服條件，在三維主應(yīng)力空間中的屈服函數(shù)的表達(dá)式如(4)所示，Drucker-Prager屈服準(zhǔn)則是對Mohr-Coulomb準(zhǔn)則不收斂的改進(jìn)，以此來修正Von Mises準(zhǔn)則，ANSYS自帶DP1屈服模塊(Mohr-Coulomb準(zhǔn)則外接點外接圓)，表達(dá)式(4)中的b和sy為DP1準(zhǔn)則修正后的材料常數(shù)。

(4)

1.4 邊坡失穩(wěn)判據(jù)

在通過對X方向塑性應(yīng)變、XY方向的剪切塑性應(yīng)變、塑性應(yīng)變強度，等效塑性應(yīng)變，Von Mises塑性應(yīng)變5個塑性指標(biāo)和空間位移、x方向位移、y方向位移3個位移指標(biāo)的分析的基礎(chǔ)上(將坡體整體位移值與塑性應(yīng)力應(yīng)變值中各自的最大變化值作為研究值)，采用Gauss、Lorentz、Vigot曲線進(jìn)行非線性擬合，定量獲得突變的極大值，當(dāng)最大變化值達(dá)到極大值時，邊坡發(fā)生破壞，為此當(dāng)位移或塑性指標(biāo)無限靠近極大值時，邊坡極限平衡狀態(tài)所對應(yīng)的折減系數(shù)，即為安全系數(shù)。

2 邊坡地質(zhì)條件

實例邊坡位于湖南省懷化市麻陽苗族自治縣某集鎮(zhèn)。集鎮(zhèn)分為侵蝕堆積河谷階地地貌和剝蝕構(gòu)造丘陵地貌兩種地貌區(qū)。集鎮(zhèn)沿錦江兩岸的河谷階地地貌區(qū)地段地勢較為平坦，高程約為190～200 m，河谷地段寬度約為850～1 000 m。沿錦江兩岸河谷地段之外為低山丘陵區(qū)，高程一般在230～260 m，山坡坡形呈弧形狀，自然坡度45°～55°，山間溝谷長且寬緩。

圖1 邊坡地理位置圖

該邊坡屬于人工開挖非均質(zhì)邊坡，邊坡位于錦江的轉(zhuǎn)彎處的右岸(圖1)，邊坡體縱長約135～136 m，高約41 m，坡腳高程約215 m，自然坡度約45°～50°，在坡腳有人工修建道路和開挖山體所形成的寬約51 m，高約15 m高陡邊坡，坡度為85°。邊坡地層分為4層，從上往下，第1層為新生界第4系全新統(tǒng)Qh粉質(zhì)粘土，為紫紅色，硬塑狀，以黏粒及粉粒為主，含少量高嶺土，稍有光澤，層理特征不明顯，無搖震反應(yīng)，干強度低，韌性差。該層曾發(fā)生土質(zhì)滑坡，滑坡后緣有約15～20 mm拉張裂隙，滑體兩側(cè)有剪切裂隙。第2～第4層全為中生界白堊系上統(tǒng)K1L的泥質(zhì)粉砂巖，軟弱～極軟薄厚層狀，為紫紅色，巖石硬度低，巖石抗風(fēng)化能力差，節(jié)理極發(fā)育，傾角近乎水平。第5層為基巖。該邊坡粉質(zhì)黏土層曾發(fā)生土質(zhì)滑坡，現(xiàn)基本穩(wěn)定。具體建模所需地質(zhì)參數(shù)見表1。

表1 邊坡地層參數(shù)表Table1 Formationparametertableofslop地層巖性平均層厚/m彈性模量/Pa泊松比密度/(g·m-3)粘聚力/Pa內(nèi)摩擦角/(°)1a3.01.00E+070.317002640019.92b10.01.00E+070.318502860023.93c13.01.00E+070.319333710029.54d15.02.00E+070.2820007500036.15e20.01.00E+080.25224300注:a為粉質(zhì)黏土,b為全風(fēng)化泥質(zhì)粉砂巖,c為中-強風(fēng)化泥質(zhì)粉砂巖,d為強風(fēng)化泥質(zhì)粉砂巖,e為中風(fēng)化泥質(zhì)粉砂巖;地層參數(shù)是均處于天然狀態(tài)下的。

3 ANSYS三維邊坡模型建立

取人工開挖地段的寬度為模型邊坡寬度，將模型簡化成整體高61 m，寬41 m的幾何實體，具體參數(shù)見圖2，各地層厚度按表1設(shè)置，定義單元為solid45，采用ANSYS的Drucker-Prager模塊，網(wǎng)格大小為1.5 m，網(wǎng)格掃掠(sweep)成六面體單元，生成40 018個單元(Element)，共有44 030個節(jié)點(NODE)，模型底面自由度全約束，側(cè)面法向約束，表面自由，模型法向施加重力荷載。定義材料屬性：所計算的模型為5層結(jié)構(gòu)，1～4層假定為理想彈塑性材料，第5層(中風(fēng)化泥質(zhì)粉砂巖層)假定為理想彈性材料，其中對1～4層進(jìn)行強度折減，折減系數(shù)分別為0.7、0.8、0.9、0.95、1，1.01、1.02、1.03、1.04、1.05、1.06、1.07、1.08、1.09、1.1、1.11、1.12、1.13、1.14、1.15。(對折減系數(shù)在1～1.15之間重點分析，得出邊坡每一步的內(nèi)在變形機制)根據(jù)式(1)～式(3)，對材料的初始抗剪強度參數(shù)和彈性參數(shù)進(jìn)行折減，將折減的數(shù)值代入計算。求解器設(shè)置：求解打開自動步長，時間設(shè)為1，采用PCG求解器進(jìn)行求解，求解子步數(shù)為20，打開Nonline框進(jìn)行非線性搜索，在Analysis Options打開Large Displacement Static大變形。

圖2 三維邊坡模型圖

4 計算結(jié)果分析

4.1 傳統(tǒng)判據(jù)

4.1.1計算收斂性判據(jù)

從表2可以看出，當(dāng)折減系數(shù)達(dá)到1.1時，計算仍收斂，邊坡處于平衡狀態(tài)，當(dāng)折減系數(shù)達(dá)到1.11時，計算不收斂，邊坡已經(jīng)失穩(wěn)破壞，此時的安全系數(shù)介于1.1～1.11。

表2 ANSYS計算收斂性結(jié)果表Table2 ANSYSConvergenceResultsTable折減系數(shù)0.70.80.911.011.021.031.041.051.06收斂情況收斂收斂收斂收斂收斂收斂收斂收斂收斂收斂折減系數(shù)1.071.081.091.11.111.121.131.141.15收斂情況收斂收斂收斂收斂不收斂不收斂不收斂不收斂不收斂

4.1.2塑性區(qū)貫通性判據(jù)

在巖土工程中，常用等效塑性應(yīng)變來描述巖土體的剪應(yīng)變。根據(jù)等效塑性應(yīng)變云圖(見圖3)，由于邊坡曾經(jīng)失穩(wěn)，在第一層粉質(zhì)黏土中發(fā)生過土質(zhì)滑坡。當(dāng)折減系數(shù)為0.7～0.8時，巖土體材料參數(shù)被加強，一定程度上“還原”土質(zhì)滑坡滑移狀態(tài)，最大塑性應(yīng)變出現(xiàn)在坡面上，坡腳出現(xiàn)塑性應(yīng)變區(qū)，滑坡后緣位置與現(xiàn)場情況符合；折減系數(shù)為0.9～1.0時，土質(zhì)滑坡已經(jīng)滑移，現(xiàn)基本穩(wěn)定，與現(xiàn)場情況符合，最大塑性應(yīng)變區(qū)轉(zhuǎn)移至坡腳上端1.5 m處；折減系數(shù)為1.02時等效塑性應(yīng)變區(qū)向坡體內(nèi)部延伸，因坡體內(nèi)部為泥質(zhì)粉砂巖，延伸區(qū)成線性條帶狀，折減系數(shù)為1.04～1.06時，坡體內(nèi)部出現(xiàn)塑形應(yīng)變區(qū)，與坡腳向坡體內(nèi)部延伸的塑形區(qū)連接，此時計算收斂，邊坡處于蠕變階段；折減系數(shù)為1.08～1.10時，塑形區(qū)向坡體延伸，形成貫通，貫通面尾部未到自由表面，此時計算收斂，坡體處于蠕變階段，折減系數(shù)為1.11～1.13時，塑形區(qū)形成貫通，貫通面到達(dá)自由表面，計算不收斂，坡體處于破壞變形階段，安全系數(shù)介于1.1～1.11，為巖質(zhì)滑坡，滑裂面為如圖3所示。結(jié)合等效塑性應(yīng)變矢量圖(見圖4)，邊坡先是發(fā)生土質(zhì)滑坡，隨著折減系數(shù)的增大然后發(fā)生巖質(zhì)滑坡，巖質(zhì)滑坡體方量約15 000～16 000 m3，傳統(tǒng)判據(jù)所得安全系數(shù)介于1.1～1.11。

(a)Ki=0.7

(i)Ki=1.1

4.2 新式突變判據(jù)

4.2.1位移指標(biāo)判據(jù)

從圖5中可以看出，4種位移都存在突變與屈服階段，3種曲線均有明顯的極值點，其中Z方向位移非常小，變化幅度小，可忽略不計?？臻g位移、X方向位移、Y方向位移曲線趨勢相近，其中空間位移與Y方向位移曲線數(shù)值基本相等，當(dāng)折減系數(shù)K增大到極值點之前均發(fā)生位移突變，繼續(xù)增大折減系數(shù)，位移曲線均發(fā)生了屈服，因此可以定性地以這3種位移方式下的位移曲線作為研究對象，采用最小二乘法擬合來定量確定極值點處的安全系數(shù)。假定位移與折減系數(shù)關(guān)系滿足Gauss、Lorentz、Vigot曲線方程，采用最小二乘法對該曲線進(jìn)行擬合，Gauss曲線方程為：

(4)

Lorentz曲線方程為：

(5)

Vigot曲線方程為：

Vigot曲線為Gauss曲線型卷積Lorentz曲線形型，與Gauss、Lorentz具有相同的單調(diào)性，所以當(dāng)Ki=Kic時，此時的U→極大值，即認(rèn)為邊坡處于極限平衡狀態(tài)。采用最小二乘法擬合，擬合結(jié)果的可靠度用擬合相關(guān)系數(shù)表示，相關(guān)系數(shù)越接近1，殘差平方和越接近0，擬合效果越好。

圖6～圖8給出了不同位移方式下曲線擬合情況，并根據(jù)式Ki=Kic求解對應(yīng)的安全系數(shù)，擬合結(jié)果和安全系數(shù)列于表3中。分析表中的擬合結(jié)果發(fā)現(xiàn)：Gauss、Lorentz曲線對X方向位移指標(biāo)擬合不收斂，擬合效果差，Vigot曲線對3種位移方式擬合效果都較好，特別是空間總位移，3種位移方式下得到安全系數(shù)介于1.123 6～1.133 9，相差較小。根據(jù)最佳擬合結(jié)果得出的安全系數(shù)為1.123 8，將安全系數(shù)1.123 8代入ANSYS計算所得結(jié)果不收斂，邊坡發(fā)生破壞。從Ki=0.7到Ki=1.123 8的X方向位移云圖可以看出，大位移變形區(qū)由坡體內(nèi)部過渡到坡體表面，在坡體內(nèi)部形成，明顯滑移面，位置與等效塑性應(yīng)變云圖(如圖3)貫通區(qū)吻合。

圖5 位移與折減系數(shù)的關(guān)系曲線

圖6 空間總位移曲線擬合

圖7 X方向位移曲線擬合

圖8 Y方向位移曲線擬合

(a)Ki=0.7

表3 位移曲線擬合結(jié)果Table3 Curvefitresultofdisplacement擬合對象擬合方式擬合參數(shù)U0KicAw相關(guān)系數(shù)殘差平方和安全系數(shù)Gauss擬合3.91321.12450.09060.06970.84900.03691.1245空間位移Lorentz擬合3.79831.12360.15300.07950.90710.02271.1236Vigot擬合3.79071.12380.15806.3184E-9[1]0.0826[2]0.90700.02431.1238Gauss擬合不收斂X方向位移Lorentz擬合不收斂Vigot擬合0.51941.13390.21391.5855E-06[1]0.0970[2]0.92410.02771.1339Gauss擬合3.91491.12440.08980.06900.84890.03691.1244Y方向位移Lorentz擬合3.79871.12360.15270.07930.90700.02271.1236Vigot擬合3.67981.12870.23461.7480E-06[1]0.1258[2]0.88960.02891.1287注:[1]代表系數(shù)WG;[2]代表系數(shù)WL

4.2.2塑性區(qū)指標(biāo)判據(jù)

從圖11中可以看出，9種塑性區(qū)指標(biāo)都存在突變階段與屈服階段，由于9種曲線均有明顯的極值點，其中Y、Z方向塑性應(yīng)變、XZ、YZ剪切應(yīng)變非常小，變化幅度小，可以忽略不計。對其余5種塑性指標(biāo)與折減系數(shù)關(guān)系的分析依然采用Gauss、Lorentz、Vigot曲線擬合，圖12～圖16給出了不同塑性區(qū)指標(biāo)下曲線擬合情況，并根據(jù)式Ki=Kic求解對應(yīng)安全系數(shù)，擬合結(jié)果和安全系數(shù)分析見表4，分析表中的擬合結(jié)果發(fā)現(xiàn)：Vigot曲線對5種塑性指標(biāo)都擬合，且擬合效果最優(yōu)，Gauss、Lorentz曲線只對塑性應(yīng)變強度不擬合，其余4種擬合效果較好，5種塑性指標(biāo)擬合得到的安全系數(shù)介于1.127 8～1.129 6，相差非常小，可以忽略不計，根據(jù)擬合效果，將X方向塑性應(yīng)變擬合所得結(jié)果1.129 6作為安全系數(shù)。由X方向塑性應(yīng)變云圖可以看出，隨著折減系數(shù)增大，塑性區(qū)發(fā)生貫通。

圖11 9種塑性區(qū)指標(biāo)與折減系數(shù)關(guān)系曲線

4.4 失穩(wěn)判據(jù)的對比分析

通過對比傳統(tǒng)判據(jù)和新式判據(jù)的安全系數(shù)，可以發(fā)現(xiàn)：

a.傳統(tǒng)的計算收斂性與塑性區(qū)貫通性判據(jù)：2種判據(jù)原理簡單，而且在ANSYA有限元計算過程中基本同步，折減系數(shù)都介于1～1.1。其值明顯比新式判據(jù)所得安全系數(shù)小，這是因為計算結(jié)果與網(wǎng)格劃分、收斂準(zhǔn)則的參數(shù)設(shè)置有關(guān)，具有一定的人為自由性，計算收斂性與模型參數(shù)設(shè)置也有關(guān)，如果模型參數(shù)取值不當(dāng)，也可能導(dǎo)致結(jié)果提前不收斂，則所得安全系數(shù)偏小。對于理想彈塑性單元，當(dāng)其應(yīng)力達(dá)到屈服狀態(tài)后，該單元周圍依然存在彈塑性約束，那么這些約束會限制這個單元的塑性應(yīng)變發(fā)展，坡體依然是穩(wěn)定狀態(tài)，只有繼續(xù)增大折減系數(shù)，坡體內(nèi)單元的應(yīng)變逐漸增大，塑性單元逐漸增多，才形成整體滑坡，若僅以貫通作為判斷標(biāo)準(zhǔn)，則所得安全系數(shù)偏小。

圖12 塑性應(yīng)變強度曲線擬合

圖13 Von Mises塑性應(yīng)變曲線擬合

圖14 等效塑性應(yīng)變曲線擬合

圖15 X方向塑性應(yīng)變

圖16 XY剪切塑性應(yīng)變

b.位移指標(biāo)突變與塑性指標(biāo)突變判據(jù)得到的安全系數(shù)十分接近，兩者都可以通過最小二乘法對指標(biāo)與折減系數(shù)關(guān)系進(jìn)行非線性擬合，且假定它們之間的關(guān)系滿足Gauss、Lorentz、Vigot曲線方程，從而確定邊坡位移指標(biāo)與塑性指標(biāo)與折減系數(shù)的定量關(guān)系，能夠表征邊坡在折減系數(shù)達(dá)到一定程度后產(chǎn)生多種位移與塑性應(yīng)以應(yīng)變的極大值，然后出現(xiàn)屈服階段，用來確定安全系數(shù)具有明確的物理意義，建議將兩者結(jié)合分析。

(a)Ki=1

表4 塑性指標(biāo)曲線擬合結(jié)果表Table4 Tableofplasticitycurvefittingresults擬合對象擬合方式擬合參數(shù)U0KicAw相關(guān)系數(shù)殘差平方和安全系數(shù)Gauss擬合不收斂塑性應(yīng)變強度Lorentz擬合不收斂Vigot擬合0.04361.12930.02510.0430[1]0.0085[2]0.91260.00331.1293Gauss擬合0.02921.12930.01320.04040.91190.00101.1293vonMises塑性應(yīng)變Lorentz擬合0.01171.12920.02150.04840.89820.00121.1292Vigot擬合0.02581.12930.01470.0423[1]0.0100[1]0.91260.00111.1293Gauss擬合0.02921.12930.01330.04030.91180.00111.1293等效塑性應(yīng)變Lorentz擬合0.01161.12920.02170.04830.89800.00121.1292Vigot擬合0.02591.12930.01470.0423[1]0.0097[2]0.91250.00111.1293Gauss擬合0.01801.12960.01120.04070.91570.00071.1296X方向塑性應(yīng)變Lorentz擬合0.00341.12960.01820.04880.89810.00081.1296Vigot擬合0.01661.12960.01180.0452[1]0.0052[2]0.91590.00081.1296Gauss擬合0.02441.12790.00670.03910.89770.00031.1279XY剪切塑性應(yīng)變Lorentz擬合0.01541.12780.01090.04630.89690.00031.1278Vigot擬合0.01991.12780.00860.0322[1]0.0237[2]0.90200.00031.1278注:[1]代表系數(shù)WG;[2]代表系數(shù)WL

c.采用傳統(tǒng)極限平衡法進(jìn)行計算驗證，其中摩根斯頓-普萊斯法所得安全系數(shù)為1.132，Janbu法所得安全系數(shù)為1.119，安全系數(shù)十分接近，滑面位置如圖18所示，與新式突變判據(jù)所得滑面位置吻合。

5 結(jié)論

a.采用有限元強度折減法，以等效塑性應(yīng)變?yōu)榛A(chǔ)，在一定程度上模擬了邊坡從土質(zhì)滑坡到巖質(zhì)滑坡等效塑性應(yīng)變?nèi)^程；采用傳統(tǒng)判據(jù)分析三維邊坡穩(wěn)定性，由計算收斂性所得安全系數(shù)為1.1～1.11，由塑性區(qū)貫通性判據(jù)所得安全系數(shù)為1.1～1.11，兩者結(jié)果偏于保守。

(a)摩根斯頓-普萊斯法

b.采用新式突變判據(jù)，3種位移值和5種塑性應(yīng)力應(yīng)變值都存在突變與屈服階段，采用最小二乘法對邊坡位移值和塑性應(yīng)力應(yīng)變值與折減系數(shù)關(guān)系進(jìn)行非線性擬合，假定它們之間的關(guān)系滿足Gauss、Lorentz、Vigot曲線方程，效果達(dá)到預(yù)期，其中Vigot曲線擬合效果最好，且應(yīng)用到此邊坡不存在不收斂情況，因此Vigot曲線方程擬合為該新式突變判據(jù)理想的擬合方式。

c.新式突變判據(jù)所得安全系數(shù)分別為：位移指標(biāo)突變判據(jù)的安全系數(shù)取自空間總位移曲線Vigot擬合結(jié)果，為1.123 8；塑性指標(biāo)突變判據(jù)的安全系數(shù)取自X方向塑性應(yīng)變曲線Vigot擬合結(jié)果為1.129 6。所得安全系數(shù)十分接近，且均達(dá)到極限平衡狀態(tài)采用極限平衡法對比驗證，結(jié)果可行，因此建議將兩者結(jié)合分析。