近場動(dòng)力學(xué)方法在邊坡穩(wěn)定性分析中的應(yīng)用

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

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

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

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

1 近場動(dòng)力學(xué)強(qiáng)度折減法

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

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

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

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

(1)

(2)

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

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

Yn+1=Yn+Δu

(3)

式中：Δu為當(dāng)前時(shí)間步n的位移增量。

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

G=F··F-1

(4)

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

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

γ=(G+GT)/2

ω=(G-GT)/2

(5)

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

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

Δσ=De∶γ

(6)

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

σn+1=σ∧n+Δσ

(7)

σ∧n=RT·σn·R

(8)

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

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

(9)

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

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

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

(10)

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

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

(11)

(12)

基于非常規(guī)態(tài)的近場動(dòng)力學(xué)運(yùn)動(dòng)平衡方程為

(13)

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

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

1. 2 強(qiáng)度折減法

強(qiáng)度折減法的基本原理是通過不斷減小材料的強(qiáng)度參數(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)摩擦角。

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

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

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

2. 1 摩爾-庫侖屈服準(zhǔn)則

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

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

(16)

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

對摩爾-庫侖屈服準(zhǔn)則進(jìn)行線性化處理[20],令X=σx-σy,Y=2τxy,R=2ccosφ-(σx+σy)cosφ，在XOY坐標(biāo)系中，屈服函數(shù)可表示成一個(gè)圓，可以使用外接多邊形進(jìn)行逼近。當(dāng)多邊形邊數(shù)為p時(shí)，線性化后的屈服函數(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)流動(dòng)法則的彈塑性矩陣

(19)

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

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

(20)

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

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

(21)

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

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

Δσ=Dep∶γ

(22)

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

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

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

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

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

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

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

圖3 邊坡土體的近場動(dòng)力學(xué)計(jì)算模型圖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粉細(xì)砂16908．00．29261

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

保持基坑的開挖深度不變，為17.0 m，將邊坡比分為1∶2和1∶2.5兩種工況，使用強(qiáng)度折減法針對不同的折減系數(shù)FOS(1.00、1.15、1.25、1.35、1.50)分別計(jì)算得到了邊坡的最大位移以及邊坡的塑性應(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的邊坡土體，當(dāng)折減系數(shù)(FOS)為1.15時(shí)，邊坡最大位移曲線就已經(jīng)出現(xiàn)了不收斂，因此可以認(rèn)為此時(shí)邊坡的安全系數(shù)小于1.15。由圖4(b)可見，對于邊坡比為1∶2.5的邊坡土體，當(dāng)折減系數(shù)為1.15時(shí)，邊坡最大位移曲線仍然是收斂的，但當(dāng)折減系數(shù)為1.25時(shí)，邊坡最大位移曲線出現(xiàn)了不收斂，因此可以認(rèn)為此時(shí)邊坡的安全系數(shù)小于1.25。

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

圖5 邊坡比為1∶2、t=4.5 s時(shí)邊坡的塑性屈服區(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時(shí)邊坡的塑性屈服區(qū)域 分布圖Fig.6 Plastic yield region of the slope when the slope ratio is 1∶2.5 and t=4.5 s

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

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

圖7 邊坡比為1∶2、FOS=1.0時(shí)邊坡的水平位移 分布圖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é) 論

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

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

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

[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] 朱大勇,盧坤林,臺(tái)佳佳,等.基于數(shù)值應(yīng)力場的極限平衡法及其工程應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào),2009,28(10):1969-1975.

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

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

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

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

[8] 李垠,蘇凱,李杰.Mohr-Coulomb等面積圓屈服準(zhǔn)則在邊坡穩(wěn)定分析中的應(yīng)用[J].大地測量與地球動(dòng)力學(xué),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]黃丹,章青,喬丕忠,等.近場動(dòng)力學(xué)方法及其應(yīng)用[J].力學(xué)進(jìn)展,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].巖石力學(xué)與工程學(xué)報(bào),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]張魯渝,劉東升,時(shí)衛(wèi)民.擴(kuò)展廣義Drucker-Prager屈服準(zhǔn)則在邊坡穩(wěn)定分析中的應(yīng)用[J].巖土工程學(xué)報(bào),2003,25(2):216-219.