廖 洋，劉立勝，2，劉齊文，賴 欣，池晴佳，邵愛軍

(1.武漢理工大學力學系，湖北 武漢 430070；2.武漢理工大學材料復合新技術(shù)國家重點實驗室，湖北 武漢 430070；3.湖北省江漢運河航道管理處，湖北 武漢 430032)

邊坡的穩(wěn)定性分析在邊坡防護以及施工安全等方面具有重要的意義[1-2]。傳統(tǒng)的邊坡穩(wěn)定性分析一般采用極限平衡法[3-5]和有限元強度折減法[6-8]。極限平衡法在工程中有著廣泛的應(yīng)用，但是由于其引入了條間力關(guān)系假設(shè)來消除邊坡穩(wěn)定問題的靜不定性，故只能得到特定條件下單一的邊坡安全系數(shù)。一般情況下，極限平衡法得到的既不是上限解也不是下限解。不同于極限平衡法僅對邊坡的宏觀安全特性進行描述，有限元方法通過借助強度折減技術(shù)，不僅能得到邊坡土體的應(yīng)力及變形分布規(guī)律，還可以得到邊坡的安全系數(shù)，以及邊坡失穩(wěn)時的塑性屈服面，即邊坡滑移面。有限元強度折減法在計算邊坡穩(wěn)定性問題時考慮了邊坡土體的非線性本構(gòu)關(guān)系、復雜的載荷，甚至整個施工過程，能夠較為真實地反映邊坡土體的應(yīng)力、變形隨時間的變化情況。

近場動力學方法是Silling等[9-10]提出的一種基于非局域平均思想的一類微分-積分方法，其在處理裂紋、滑移等不連續(xù)問題方面具有較好的適應(yīng)性。起初的近場動力學方法是一種基于鍵的非局域方法，比較簡單，只能模擬簡單的本構(gòu)關(guān)系；其后提出的基于非常規(guī)態(tài)的近場動力學方法通過引入矩量法[11]，可以得到質(zhì)點的任意時刻的變形梯度以及應(yīng)變。基于非常規(guī)態(tài)的近場動力學方法[12-15]作為一種非局域的無網(wǎng)格粒子類方法，不僅具備有限元方法的優(yōu)點，而且在計算非線性問題時還具有更好的收斂性以及對不連續(xù)問題的適用性。當邊坡處于臨界失穩(wěn)狀態(tài)發(fā)生滑移時，有限元方法的計算容易出現(xiàn)不收斂現(xiàn)象，而近場動力學方法不僅能計算得到邊坡達到臨界狀態(tài)的收斂解，還能模擬邊坡的滑移過程。此外，基于非常規(guī)態(tài)的近場動力學方法可以用來實現(xiàn)土體的非線性本構(gòu)模型，如Drucker-Prager屈服模型、摩爾-庫侖屈服模型等。本課題組在土體的爆炸[16]以及邊坡穩(wěn)定性分析[17]方面已做了一些探索性的工作。在此基礎(chǔ)上，本文基于摩爾-庫侖屈服準則以及非關(guān)聯(lián)流動法則，利用近場動力學強度折減法對引江濟漢工程中出口船閘的基坑邊坡的穩(wěn)定性進行了分析。

1 近場動力學強度折減法

1. 1 基于非常規(guī)態(tài)的近場動力學方法

在基于非常規(guī)態(tài)的近場動力學方法中，某一質(zhì)點的變形狀態(tài)與其自身和鄰域內(nèi)其他質(zhì)點的位移差有關(guān)。任意質(zhì)點i與其鄰域內(nèi)質(zhì)點j之間的相互作用力被稱作力態(tài)，分別為T和T′。

圖1 基于非常規(guī)態(tài)的近場動力學方法中質(zhì)點及其鄰域Fig.1 Particle in non-ordinary state-based peridynamics and its neighbourhood

任意質(zhì)點i的變形梯度張量F，可由下式計算:

(1)

(2)

式中：Y為質(zhì)點i與j之間的變形后的相對位置矢量;X為質(zhì)點i與j之間的初始相對位置矢量;V′為質(zhì)點j的體積;K為質(zhì)點i的形狀張量;w(|X|)為質(zhì)點j的影響函數(shù);H為質(zhì)點i的鄰域范圍。

根據(jù)連續(xù)介質(zhì)力學理論，物質(zhì)體在時間步n+1的形態(tài)Yn+1可由前一時間步n的形態(tài)Yn計算得到:

Yn+1=Yn+Δu

(3)

式中：Δu為當前時間步n的位移增量。

因此，從時間步n到n+1的變形梯度增量G的計算公式為

G=F··F-1

(4)

式中：F·為變形梯度F的時間導數(shù)，即速度梯度。

物質(zhì)體的運動可以分為應(yīng)變增量和剛體轉(zhuǎn)動增量兩部分：

γ=(G+GT)/2

ω=(G-GT)/2

(5)

式中：γ為應(yīng)變增量;ω為剛體轉(zhuǎn)動增量。

因此，有效應(yīng)力增量Δσ可通過彈性矩陣De與應(yīng)變增量γ進行縮并計算得到：

Δσ=De∶γ

(6)

根據(jù)Hughes-Winget理論[18]，質(zhì)點當前時刻n+1的柯西應(yīng)力張量，可由下式計算:

σn+1=σ∧n+Δσ

(7)

σ∧n=RT·σn·R

(8)

其中，σn為時刻n的應(yīng)力狀態(tài)；σ∧n為σn從時刻n到時刻n+1發(fā)生剛體轉(zhuǎn)動的結(jié)果;R為轉(zhuǎn)動張量，可以由下式計算：

R=I+(I-0.5ω)-1·ω

(9)

式中：I為二階單位張量。

根據(jù)柯西應(yīng)力張量σn+1，可計算求得第一類皮奧拉-基爾霍夫非對稱應(yīng)力張量Pn+1：

Pn+1=Jσn+1F-T

(10)

式中：J為變形梯度F的行列式。

力態(tài)T和T′的計算形式如下：

(11)

(12)

基于非常規(guī)態(tài)的近場動力學運動平衡方程為

(13)

式中：ü為質(zhì)點i在t時刻的位移二階導數(shù)，即加速度；x為質(zhì)點i的初始坐標；ρ為質(zhì)點的材料密度；b為質(zhì)點i在t時刻受到的外力矢量作用。

根據(jù)基于非常規(guī)態(tài)的近場動力學運動平衡方程，可計算得到質(zhì)點i當前時刻n的加速度，最后使用歐拉法、蛙跳法或Velocity-Verlet積分法，可求得下一時刻的質(zhì)點位移，從而完成變形體的狀態(tài)更新。

1. 2 強度折減法

強度折減法的基本原理是通過不斷減小材料的強度參數(shù)直到邊坡開始出現(xiàn)失穩(wěn)。邊坡失穩(wěn)臨界狀態(tài)所對應(yīng)的折減系數(shù)即是該邊坡的安全系數(shù)，則有

ct=cini/FOS

(14)

tanφt=tanφini/FOS

(15)

式中：FOS為折減系數(shù)；ct和φt分別為邊坡臨界狀態(tài)的內(nèi)聚力和內(nèi)摩擦角；cini和φini分別為邊坡初始狀態(tài)的內(nèi)聚力和內(nèi)摩擦角。

在有限元強度折減法中，當材料強度進行折減后，若有限元計算收斂，則認為邊坡是穩(wěn)定的；當材料強度折減到一定范圍時，有限元計算開始出現(xiàn)不收斂的現(xiàn)象，此時的折減系數(shù)一般視作邊坡的安全系數(shù)。

在近場動力學中引入強度折減法，其基本的計算思想與有限元法相同；不同之處是，近場動力學方法是以邊坡系統(tǒng)的最大位移作為判斷依據(jù)，當邊坡處于穩(wěn)定狀態(tài)時，不同折減系數(shù)下邊坡的最大位移值基本不變，一旦邊坡失穩(wěn)，邊坡的最大位移值會發(fā)生較大的變化。因此，通過觀察邊坡最大位移值的變化，可以在近場動力學強度折減法中捕捉到邊坡的安全系數(shù)。

2 摩爾-庫侖屈服準則與非關(guān)聯(lián)流動法則

2. 1 摩爾-庫侖屈服準則

摩爾-庫侖屈服準則常被用作混凝土以及巖土類材料的塑性判據(jù)。對于平面應(yīng)變問題[19]，摩爾-庫侖準則可表示為

F=-[2ccosφ-(σx+σy)sinφ]2+(σx-σy)2+(2τxy)2=0

(16)

式中：c和φ分別為土體的內(nèi)聚力和內(nèi)摩擦角，σx、σy和τxy為應(yīng)力分量。

對摩爾-庫侖屈服準則進行線性化處理[20],令X=σx-σy,Y=2τxy,R=2ccosφ-(σx+σy)cosφ，在XOY坐標系中，屈服函數(shù)可表示成一個圓，可以使用外接多邊形進行逼近。當多邊形邊數(shù)為p時，線性化后的屈服函數(shù)為

F=Mkσx+Nkσy+Pkτxy-2ccosφ

(17)

其中:

Mk=cosαk+sinφ

Nk=sinφ-cosαk

Pk=2sinαk

αk=2πk/p(k=1,2,…,p)

(18)

2. 2 基于非關(guān)聯(lián)流動法則的彈塑性矩陣

(19)

式中：dλ為塑性應(yīng)變增量的大小，它是一個非負的標量；Q為塑性勢函數(shù)；σij為應(yīng)力張量。

在非關(guān)聯(lián)流動法則中，塑性勢函數(shù)Q與屈服函數(shù)F不相等，即Q≠F。對于理想的彈塑性材料，硬化參數(shù)A為零，則應(yīng)力-應(yīng)變增量表達式為

(20)

式中：dεij為總的應(yīng)變增量;Dep為彈塑性矩陣;De、Dp為塑性矩陣，其形式[21]如下:

Dp=?F?σTDe·De?Q?σ?F?σTDe?Q?σ

(21)

2. 3 基于摩爾-庫侖屈服準則的應(yīng)力更新

若屈服函數(shù)F<0，則材料處于彈性階段，只需根據(jù)近場動力學的一般計算流程更新應(yīng)力即可；若屈服函數(shù)F≥0，則材料發(fā)生了塑性變形，此時首先需要根據(jù)公式(20)和(21)計算出當前的彈塑性矩陣Dep，用以代替公式(6)中的彈性矩陣De，其應(yīng)力增量為

Δσ=Dep∶γ

(22)

3 工程實例應(yīng)用與分析

本文以引江濟漢工程邊坡為例，基于摩爾-庫侖屈服準則以及非關(guān)聯(lián)流動法則，利用近場動力學強度折減法對引江濟漢工程中出口船閘的基坑邊坡穩(wěn)定性進行了分析。引江濟漢工程出口船閘位于潛江市高石碑鎮(zhèn)，閘址區(qū)地形平坦，地面高程為32.5 m。

3. 1 地層結(jié)構(gòu)分布

該地區(qū)的地層結(jié)構(gòu)分布如下：①亞砂土，結(jié)構(gòu)松散—較密實，不均質(zhì)，砂感明顯，黏性差，厚度為0～2.1 m，底板高程為29.8～31.0 m；②亞黏土，可塑—軟塑，黏性較強，厚度為0～2.7 m，底板高程為26.06～27.63 m；③黏土，可塑為主，局部硬塑，土質(zhì)均一細膩，黏性強，厚度為6.2～16.2 m，底板高程為8.6～11.5 m；④粉細砂，飽和，稍密—中密，厚度為2.6～11.5 m。

3. 2 邊坡土質(zhì)及計算模型

基坑邊坡地層上部為亞黏土、亞砂土、黏土，局部為淤泥質(zhì)土，中下部為細砂，邊坡土體較軟弱。閘址地面的高程為32.5 m，最大開挖深度約17 m，邊坡比取值范圍為1∶2～1∶2.5。邊坡土質(zhì)分層情況以及邊坡尺寸見圖2，邊坡土體近場動力學計算模型見圖3，邊坡土體的材料參數(shù)見表1。

圖2 邊坡土質(zhì)分層情況以及邊坡比為1∶2時的尺寸圖Fig.2 Soil layers and size of the slope when the slope ratio is 1∶2

圖3 邊坡土體的近場動力學計算模型圖Fig.3 Peridynamic model of the slope

材料土體密度/(kg·m-3)壓縮模量/(MPa)泊松比內(nèi)摩擦角/(°)內(nèi)聚力/(kPa)亞砂土19407．40．32206亞黏土18904．40．401515黏土18404．00．421316粉細砂16908．00．29261

3. 3 計算結(jié)果與分析

保持基坑的開挖深度不變，為17.0 m，將邊坡比分為1∶2和1∶2.5兩種工況，使用強度折減法針對不同的折減系數(shù)FOS(1.00、1.15、1.25、1.35、1.50)分別計算得到了邊坡的最大位移以及邊坡的塑性應(yīng)變分布情況,見圖4、圖5和圖6。

圖4 不同折減系數(shù)對應(yīng)的邊坡最大位移曲線Fig.4 Maximum displacements of the slope corresponding to different reduction factors

由圖4(a)可見，對于邊坡比為1∶2的邊坡土體，當折減系數(shù)(FOS)為1.15時，邊坡最大位移曲線就已經(jīng)出現(xiàn)了不收斂，因此可以認為此時邊坡的安全系數(shù)小于1.15。由圖4(b)可見，對于邊坡比為1∶2.5的邊坡土體，當折減系數(shù)為1.15時，邊坡最大位移曲線仍然是收斂的，但當折減系數(shù)為1.25時，邊坡最大位移曲線出現(xiàn)了不收斂，因此可以認為此時邊坡的安全系數(shù)小于1.25。

從以上計算結(jié)果來看，邊坡的安全系數(shù)較低，邊坡土體極容易發(fā)生失穩(wěn)，因此可以進一步分析邊坡的塑性屈服區(qū)域以及最大水平位移出現(xiàn)的區(qū)域。

圖5 邊坡比為1∶2、t=4.5 s時邊坡的塑性屈服區(qū)域 分布圖Fig.5 Plastic yield region of the slope when the slope ratio is 1∶2.5 and t=4.5 s

圖6 邊坡比為1∶2.5、t=4.5 s時邊坡的塑性屈服區(qū)域 分布圖Fig.6 Plastic yield region of the slope when the slope ratio is 1∶2.5 and t=4.5 s

由圖5可見，當折減系數(shù)為1.00時，邊坡坡腳就開始出現(xiàn)了塑性屈服區(qū)域，且隨著折減系數(shù)的增加，坡腳的塑性屈服區(qū)域急劇增大，同時邊坡上部也開始出現(xiàn)塑性屈服區(qū)域。對比圖5和圖6可見，隨著邊坡比的增大，邊坡坡腳的塑性屈服區(qū)域明顯減小，但是邊坡的穩(wěn)定性仍然較差，此時應(yīng)加強坡腳的防護，這一結(jié)論與實際施工的情況是符合的。

當邊坡比為1∶2或1∶2.5時，對于不同的折減系數(shù)，邊坡的水平位移分布趨勢基本相同，圖7給出了邊坡比為1∶2、折減系數(shù)為1.0時邊坡的水平位移分布圖，圖8給出了不同邊坡比情況下邊坡最大水平位移隨折減系數(shù)的變化曲線。

圖7 邊坡比為1∶2、FOS=1.0時邊坡的水平位移 分布圖Fig.7 Horizontal displacement when the slope ratio is 1∶2 and FOS=1.0

圖8 不同邊坡比情況下邊坡最大水平位移隨折減 系數(shù)的變化曲線Fig.8 Maximum horizontal displacement of the slope with the reduction factor under different slope ratios

由圖7和圖8可見，邊坡的最大水平位移出現(xiàn)在坡腳的上端，且隨著邊坡比的增大，邊坡的最大水平位移明顯減小，邊坡趨于穩(wěn)定。

4 結(jié) 論

通過上述工程實例計算與應(yīng)用，可得到如下結(jié)論：

(1) 本文所提出的基于摩爾-庫侖屈服準則的近場動力學方法可以用于描述邊坡土體材料的力學行為。

(2) 基于非常規(guī)態(tài)的近場動力學強度折減法可以用于邊坡的穩(wěn)定性分析，所得到的計算結(jié)果與實際施工情況是一致的。

[1] 楊茜,陳勝,朱桂春.基于拉格朗日乘數(shù)法的邊坡穩(wěn)定性變分法[J].安全與環(huán)境工程,2010,17(4):111-113.

[2] 郭利娜,胡斌,胡啟晨.基于Geo Slope軟件對青蓮寺邊坡的穩(wěn)定性分析[J].安全與環(huán)境工程,2011,18(6):20-24.

[3] 朱大勇,盧坤林,臺佳佳,等.基于數(shù)值應(yīng)力場的極限平衡法及其工程應(yīng)用[J].巖石力學與工程學報,2009,28(10):1969-1975.

[4] 孫聰,李春光,鄭宏.基于整體穩(wěn)定性分析法的邊坡臨界滑動面搜索[J].巖土力學,2013,34(9):2583-2588.

[5] 欒茂田,金崇磐,林皋.土體穩(wěn)定分析極限平衡法改進及其應(yīng)用[J].巖土工程學報,1992,14(S1):20-29.

[6] 張坤勇,李廣山,李旺林,等.南水北調(diào)南干渠邊坡有限元穩(wěn)定性分析[J].河北工程大學學報(自然科學版),2016,33(4):27-32.

[7] 趙尚毅,鄭穎人,時衛(wèi)民,等.用有限元強度折減法求邊坡穩(wěn)定安全系數(shù)[J].巖土工程學報,2002,24(3):343-346.

[8] 李垠,蘇凱,李杰.Mohr-Coulomb等面積圓屈服準則在邊坡穩(wěn)定分析中的應(yīng)用[J].大地測量與地球動力學,2009,29(1):135-139.

[9] Silling S A.Reformulation of elasticity theory for discontinuities and long-range forces[J].JournaloftheMechanics&PhysicsofSolids,2000,48(1):175-209.

[10]Silling S A,Epton M,Weckner O,et al.Peridynamic states and constitutive modeling[J].JournalofElasticity,2007,88(2):151-184.

[11]Li S,Qian D,Liu W K,et al.A meshfree contact-detection algorithm[J].ComputerMethodsinAppliedMechanics&Engineering,2001,190(24/25):3271-3292.

[12]Zhou X,Wang Y,Xu X.Numerical simulation of initiation,propagation and coalescence of cracks using the non-ordinary state-based peridynamics[J].InternationalJournalofFracture,2016,201(2):213-234.

[13]Kilic B,Agwai A,Madenci E.Peridynamic theory for progressive damage prediction in center-cracked composite laminates[J].CompositeStructures,2009,90(2):141-151.

[14]Amani J,Oterkus E,Areias P,et al.A non-ordinary state-based peridynamics formulation for thermoplastic fracture[J].InternationalJournalofImpactEngineering,2016,87:83-94.

[15]黃丹,章青,喬丕忠,等.近場動力學方法及其應(yīng)用[J].力學進展,2010,40(4):448-459.

[16]Lai X,Ren B,Fan H,et al.Peridynamics simulations of geomaterial fragmentation by impulse loads[J].InternationalJournalforNumerical&AnalyticalMethodsinGeomechanics,2015,39(12):1304-1330.

[17]Lai X,Liu L S,Liu Q W,et al.Slope stability analysis by peridynamic theory[J].AppliedMechanics&Materials,2015,744/745/746:584-588.

[18]Hughes T J R,Winget J.Finite rotation effects in numerical integration of rate constitutive equations arising in large-deformation analysis[J].InternationalJournalforNumericalMethodsinEngineering,1980,15(12):1862-1867.

[19]孫聰,李春光,鄭宏,等.基于四邊形網(wǎng)格的邊坡上限有限元法[J].巖石力學與工程學報,2015,34(1):114-120.

[20]Lyamin A V,Sloan S W.Upper bound limit analysis using linear finite elements and non-linear programming[J].InternationalJournalforNumerical&AnalyticalMethodsinGeomechanics,2002,26(2):61-77.

[21]張魯渝,劉東升,時衛(wèi)民.擴展廣義Drucker-Prager屈服準則在邊坡穩(wěn)定分析中的應(yīng)用[J].巖土工程學報,2003,25(2):216-219.