基于連續(xù)-非連續(xù)方法的高填方邊坡變形穩(wěn)定性分析

司憲志，寧 宇，石 崇*，黃青富，張一平

(1.河海大學(xué) 巖土工程研究所，江蘇 南京 210098；2.中國電建集團(tuán)昆明勘測設(shè)計(jì)研究院有限公司，云南 昆明 650051)

我國西南地區(qū)多為山地丘陵地貌，在這些區(qū)域進(jìn)行機(jī)場建造時(shí)，為制造出合適的場地條件和凈空條件，勢必會進(jìn)行深挖高填工作，產(chǎn)生大量的高填方邊坡[1]，同時(shí)山地及丘陵地帶復(fù)雜的地質(zhì)條件，也給機(jī)場建設(shè)的安全性帶來了挑戰(zhàn)，填方地基及高填方體的穩(wěn)定性問題成為山區(qū)機(jī)場建設(shè)中存在的核心問題[2]，加強(qiáng)對該高填方邊坡穩(wěn)定性的研究，尋求科學(xué)、合理、有效的評價(jià)及防治措施，具有重要的現(xiàn)實(shí)意義。數(shù)值模擬由于其成本低、效果明顯的特點(diǎn)被廣泛地應(yīng)用于邊坡穩(wěn)定性分析[3]中，在連續(xù)數(shù)值仿真方法方面：陳正東[4]等基于有限差分軟件FLAC3D軟件，建立了三維邊坡數(shù)值計(jì)算模型，研究了開挖后的邊坡體的位移變形特征，并據(jù)此對開挖后的邊坡進(jìn)行了穩(wěn)定性評價(jià)；潘網(wǎng)生[5]等基于FLAC3D拉格朗日快速差分算法模擬煤層開挖，揭示邊坡變形過程及其滑動(dòng)機(jī)理；言志信[6]等以含軟弱層的錨固巖體邊坡為研究對象，利用FLAC3D軟件模擬分析了地震作用下含軟弱層巖體邊坡兩錨固界面剪切作用及其演化。但連續(xù)數(shù)值方法中單元之間需要公用節(jié)點(diǎn)才能合理描述荷載力的傳遞，所能模擬的變形量有限，當(dāng)變形增大到一定程度時(shí)，單元會發(fā)生畸變導(dǎo)致剛度矩陣無法求解，導(dǎo)致計(jì)算終止。在非連續(xù)數(shù)值仿真方法方面：李龍起[7]等利用顆粒流軟件PFC2D對土質(zhì)滑坡在持續(xù)降雨作用下的變形發(fā)展過程進(jìn)行模擬研究，揭示了該滑坡的發(fā)展及破壞模式；楊忠平[8]等用振動(dòng)臺物理模型試驗(yàn)，結(jié)合UDEC離散元分析方法，研究了頻發(fā)微震作用下典型上覆軟弱巖體邊坡的累計(jì)損傷過程、動(dòng)力響應(yīng)特征及破壞模式，Hung[9]等通過有限元分析和離散元分析，探討了同震滑坡的破壞前機(jī)制和破壞后運(yùn)動(dòng)過程與初始時(shí)間和運(yùn)動(dòng)擺度相關(guān)的顯著特征。離散單元方法中塊體或顆?？梢苑蛛x，通過每個(gè)時(shí)間步不斷地接觸判斷更新塊體或顆粒的運(yùn)動(dòng)方程，不受變形量限制，但是大量的接觸判斷導(dǎo)致計(jì)算速度降低。

綜上所述，如果能將連續(xù)-非連續(xù)耦合數(shù)值模擬方法相互結(jié)合，潛在破壞區(qū)域采用非連續(xù)方法模擬，對于基巖部分變形量小不會出現(xiàn)大變形破壞的區(qū)域則采用連續(xù)數(shù)值模擬方法分析，則既能滿足計(jì)算效率的要求，又不受變形量限制。本文基于連續(xù)-非連續(xù)理論，通過邊界墻耦合方法建立了某高填方邊坡連續(xù)-非連續(xù)數(shù)值模型，研究了該邊坡在天然工況下攔擋壩的應(yīng)力分布情況，將強(qiáng)度折減法應(yīng)用于細(xì)觀模擬，進(jìn)行了邊坡在降雨工況下的穩(wěn)定性分析，并求解了邊坡的安全系數(shù)，分析了邊坡在極端工況下可能的變形破壞情況。

1 連續(xù)非連續(xù)數(shù)值模擬方法

1.1 連續(xù)-非連續(xù)邊界耦合

在連續(xù)-非連續(xù)相互耦合分析中，有限差分軟件FLAC用來從宏觀上模擬連續(xù)區(qū)域內(nèi)的力學(xué)行為，而顆粒流軟件PFC用來從細(xì)觀上模擬離散區(qū)域內(nèi)介質(zhì)的力學(xué)行為，兩者間的相互耦合依賴于連續(xù)區(qū)域與離散區(qū)域的接觸邊界，不同區(qū)域間的計(jì)算數(shù)據(jù)是借助于Socket O/I接口來進(jìn)行相互傳輸與交換[10-11]，連續(xù)-非連續(xù)耦合計(jì)算原理見圖1。

圖1 連續(xù)-非連續(xù)耦合計(jì)算原理Fig.1 Calculation principle of continuous-discontinuous coupling method

本文采用基于邊界控制墻體的連續(xù)-非連續(xù)耦合方法，在離散區(qū)域和連續(xù)區(qū)域之間創(chuàng)建邊界墻體，離散區(qū)域的顆粒運(yùn)動(dòng)過程中作用于墻體上的接觸力和接觸彎矩，采用等效力的方法分配到邊界墻體的頂點(diǎn)上，進(jìn)而將力的信息傳遞給連續(xù)區(qū)域的單元，這些力參與到連續(xù)區(qū)域的有限差分法分析，同樣連續(xù)區(qū)域的運(yùn)動(dòng)也會帶動(dòng)邊界墻體的運(yùn)動(dòng)，進(jìn)而將位置和速度通過邊界墻傳遞給離散區(qū)域，由此實(shí)現(xiàn)離散區(qū)域與連續(xù)區(qū)域的耦合計(jì)算。

1.2 強(qiáng)度折減法

強(qiáng)度折減法通常應(yīng)用于計(jì)算邊坡的安全系數(shù)，它是通過逐漸減小邊坡體材料的強(qiáng)度參數(shù)，從而使邊坡體進(jìn)入到臨界破壞狀態(tài)來求解安全系數(shù)。對于Mohr-Coulomb破壞準(zhǔn)則來說，安全系數(shù)F根據(jù)下面的方程來定義：

(1)

(2)

式中c′為折減后的粘聚力；φ′為折減后的內(nèi)摩擦角；F′為折減系數(shù)。

通過不斷調(diào)整邊坡的巖土體強(qiáng)度指標(biāo)，即內(nèi)聚力c和內(nèi)摩擦角φ，然后對邊坡進(jìn)行穩(wěn)定性分析，隨后通過不斷增加折減系數(shù)，進(jìn)行一系列的計(jì)算直至邊坡達(dá)到極限平衡狀態(tài)，這時(shí)候得到的折減系數(shù)即為安全系數(shù)[12]。

2 工程應(yīng)用

2.1 工程概況

現(xiàn)擬建某高填方機(jī)場，場區(qū)高差近400 m，海拔1 500～1 900 m，本期規(guī)劃機(jī)場設(shè)計(jì)標(biāo)高1 743 m，挖方量約1 655×104m3，最大挖方高度約70余米，填方量約1 457×104m3，最大填方高度70余米。由于填方高度較大，存在的主要工程地質(zhì)問題為高填方邊坡的穩(wěn)定性問題以及不均勻沉降問題。

2.2 計(jì)算模型與條件

基于以上工程地質(zhì)狀況通過FLAC3D與PFC對NPH3剖面建立連續(xù)非連續(xù)耦合模型，三維模型尺寸為280 m×8 m×162 m，如圖2所示，模型分為基巖、壩體、風(fēng)化帶和土石堆積四部分。土石堆積體用離散元顆粒表示，生成顆粒33 600個(gè)，顆粒的直徑分布在0.4～0.8 m之間，土石混合體部分顆粒的孔隙率為0.55，顆粒間的接觸采用接觸粘結(jié)模型。擋土壩和基巖由于變形較小，故采用連續(xù)單元模型，由于單元的尺寸會對計(jì)算結(jié)果產(chǎn)生影響，網(wǎng)格密度越高計(jì)算結(jié)果越精細(xì)，但密度過高的網(wǎng)格，其計(jì)算效率也越低，為保證合理的計(jì)算效率，依據(jù)文獻(xiàn)[13-15]中相同尺度的模型取單元尺寸大小，單元網(wǎng)格尺寸應(yīng)取2～5 m，本模型中取3 m，通過ANSYS劃分網(wǎng)格后再將模型導(dǎo)入FLAC3D中，共生成單元網(wǎng)格7 897個(gè)，單元節(jié)點(diǎn)12 225個(gè)，對連續(xù)區(qū)域的模型使用Mohr-Coulomb準(zhǔn)則，模型的約束條件設(shè)置為左右兩側(cè)約束水平位移，底邊固定，模型的宏細(xì)觀參數(shù)標(biāo)定如表1所示。

圖2 三維數(shù)值模型Fig.2 3D numerical model

表 1 數(shù)值模型巖土體宏細(xì)觀參數(shù)

3 結(jié)果分析

3.1 天然工況

對建立的連續(xù)-非連續(xù)耦合模型賦予表1所示的宏細(xì)觀參數(shù)，并使其在自重下平衡，最終得到的天然工況下土石堆積體力鏈分布如圖3所示，由于自重作用，堆積體底部力鏈較為粗壯，頂部力鏈相對稀疏，在堆積體和風(fēng)化巖的接觸界面，由于接觸界面并不為水平，堆積體有相對滑動(dòng)趨勢，但是在壩體的作用下并無位移，因此堆積體與壩體接觸部位力鏈相對粗壯，且堆積體右側(cè)力鏈出現(xiàn)部分張拉力鏈。天然工況下堆積體的大主應(yīng)力分布

圖3 天然工況下堆積體力鏈圖Fig.3 Stacking force chain diagram under natural working conditions

如圖4所示，由于堆積體的作用，壩體左下角出現(xiàn)應(yīng)力集中，最大主應(yīng)力為3.2 MPa。

圖4 攔擋壩壩體的大主應(yīng)力等值線圖Fig.4 Contour map of high principal stress of retaining dam

3.2 降雨工況

在降雨條件下，降雨過程中及降雨后雨水滲透土體邊坡，會引起坡體非飽和帶的土體基質(zhì)吸力降低，由于雨水向下滲透以及地下水位的升高，會導(dǎo)致巖土體的容重發(fā)生變化，同時(shí)由于水對巖土體的軟化作用，使得巖土體的物理力學(xué)參數(shù)降低[16-18]，此外水會對土體產(chǎn)生滲流力與浮力的作用，因此邊坡體內(nèi)的滲流場的變化會引起應(yīng)力場的變化，同時(shí)滲流場的變化也會相應(yīng)地影響到邊坡體中的應(yīng)力場，水的流動(dòng)與土體的變形是相互制約相互作用的，在模擬中實(shí)現(xiàn)降雨中應(yīng)力場與滲流場的耦合過程較為復(fù)雜，因此為簡化降雨工況，將降雨過程簡化為對土石堆積體的細(xì)觀參數(shù)進(jìn)行折減，表2[19-26]統(tǒng)計(jì)了近年來8例滑坡中土體飽和狀態(tài)及天然狀態(tài)其物理力學(xué)參數(shù)的折減百分比，土體參數(shù)的折減范圍在4%～29%，將各案例中土體參數(shù)的折減百分比通過氣泡圖來展示(圖5)，可見折減百分比大部分在10%～20%，所統(tǒng)計(jì)滑坡土體參數(shù)折減的平均值在15.45%左右，對于宏觀土體應(yīng)對其內(nèi)聚力及內(nèi)摩擦角進(jìn)行折減，對于土體顆粒的細(xì)觀參數(shù)，根據(jù)文獻(xiàn)[11，27]的研究，并依據(jù)工程經(jīng)驗(yàn)，應(yīng)對離散顆粒間接觸的摩擦系數(shù)、接觸粘結(jié)張拉強(qiáng)度和接觸粘結(jié)強(qiáng)度折減15%。

表2 滑坡體材料折減參數(shù)

圖5 土體參數(shù)折減百分比氣泡圖Fig.5 Bubble chart of reduction percentage of soil parameters

為監(jiān)測壩體以及堆積體的位移情況，在壩體背面設(shè)置監(jiān)測線，沿堆積體坡面設(shè)置監(jiān)測點(diǎn)，由于監(jiān)測單獨(dú)顆粒具有一定的隨機(jī)性，因此設(shè)置數(shù)個(gè)顆粒團(tuán)作為一個(gè)監(jiān)測點(diǎn)，并取這幾個(gè)顆粒的平均位移作為測點(diǎn)的位移，測線及測點(diǎn)分布如圖6所示。

圖6 測線及測點(diǎn)設(shè)置Fig.6 Setup of measuring lines and points

降雨工況下壩體的位移等值線圖如圖7(a)所示，壩體位移較小，最大位移處位于壩體頂部，為3.8 mm，從壩體頂部至壩體底部位移逐漸減小，壩體背面的位移曲線如圖7(b)所示，壩體位移隨高程增加而增加，但降雨工況下位移相對較小。

圖7 降雨工況下壩體位移Fig.7 Dam displacement under rainfall condition

降雨工況下堆積體的位移如圖8(a)所示，最大位移處位于堆積體的頂部，為0.5 m。堆積體坡面沿高程方向的位移曲線如圖8(c)所示，堆積體的位移沿坡面從坡底向坡頂逐漸增加，位移與高程基本呈線性相關(guān)，出現(xiàn)這一現(xiàn)象的原因是由于降雨工況下堆積體的材料參數(shù)折減導(dǎo)致的不均勻沉降。由圖8(b) 可見，堆積體的位移主要朝下，在與風(fēng)化帶接觸界面角度的影響下，堆積體產(chǎn)生少量沿坡面方向的位移，但在壩體的阻擋作用下，未產(chǎn)生明顯的沿坡面方向位移。

3.3 強(qiáng)度折減

利用強(qiáng)度折減法對堆積體的力學(xué)參數(shù)進(jìn)行折減以尋求邊坡安全系數(shù)，對堆積體顆粒間接觸的摩擦系數(shù)、接觸粘結(jié)張拉強(qiáng)度和接觸粘結(jié)強(qiáng)度進(jìn)行折減，按照二分法依次折減，若堆積體出現(xiàn)明顯滑動(dòng)面或塌落，則判定為破壞，當(dāng)折減至臨界破壞狀態(tài)，則取安全系數(shù)為讓模型恰好出現(xiàn)破壞的折減系數(shù)。圖8為對模型采用二分法折減時(shí)堆積體的位移云圖，當(dāng)折減系數(shù)為2.0時(shí)，堆積體出現(xiàn)明顯塌落，并沿圖示曲線滑塌，則調(diào)整折減系數(shù)為1.50，堆積體并無明顯滑動(dòng)面，繼續(xù)折減直至堆積體恰好出現(xiàn)破壞，如圖9(d)和圖9(f)所示，當(dāng)折減系數(shù)為1.60時(shí)，堆積體恰好出現(xiàn)滑動(dòng)面，當(dāng)折減系數(shù)為1.59時(shí)，堆積體無滑動(dòng)面出現(xiàn)，因此判定堆積體的安全系數(shù)為1.60。

圖8 降雨工況下堆積體位移Fig.8 Displacement of accumulation body under rainfall condition

圖9 不同折減系數(shù)下堆積體位移云圖Fig.9 Displacement nephogram of accumulation under different reduction factors

當(dāng)折減系數(shù)為1.60時(shí)，堆積體顆粒的位移情況如圖10(b)所示，堆積體顆粒出現(xiàn)塌落，并沿圖示曲線滑塌，在壩體頂端滑出，但壩體無明顯位移，顆粒的位移方向如圖10(a)所示，滑塌區(qū)域的顆粒位移方向基本一致，位移最大的區(qū)域位于滑塌體中部，堆積體顆粒的力鏈如圖10(c)所示，滑塌區(qū)域頂部(如圖所示方框內(nèi))力鏈出現(xiàn)斷裂，可見滑塌區(qū)域與后方堆積體出現(xiàn)明顯的錯(cuò)動(dòng)；堆積體顆粒的接觸破壞模式如圖10(d)所示，不同的接觸狀態(tài)表明不同的接觸破壞方式。從圖中接觸粘結(jié)破壞模式可以看出，堆積體的滑動(dòng)面附近大部分顆粒接觸處于非粘結(jié)狀態(tài)和由剪切破壞導(dǎo)致的接觸非粘結(jié)狀態(tài)，產(chǎn)生的裂隙大部分位于滑動(dòng)面附近，與實(shí)際滑坡相符合。

注：圖例中0表示相鄰顆粒處于非粘結(jié)狀態(tài)，1表示接觸非粘結(jié)并且是由于張拉導(dǎo)致破壞，2表示接觸非粘結(jié)并且是由于剪切破導(dǎo)致產(chǎn)生，3表示接觸處于粘結(jié)狀態(tài)且粘結(jié)完好無損。圖10 折減系數(shù)為1.60時(shí)顆粒情況Fig.10 Particle size distribution with reduction factor of 1.60

堆積體坡面沿高程方向的顆粒位移情況圖11所示，坡面位移最大的區(qū)域位于堆積體頂部，為3.78 m，位移最小的區(qū)域位于堆積體坡面底部，為1.64 m，除堆積體底部與壩體接觸區(qū)域外，堆積體坡面的位移基本與高程成正比，結(jié)合圖10(b)，高填方邊坡在極端情況(如持續(xù)高強(qiáng)度降雨、地震等)下，潛在失穩(wěn)區(qū)域更大，發(fā)生失穩(wěn)時(shí)滑面更深，范圍更廣，更應(yīng)做好邊坡的支擋工作。

圖11 折減系數(shù)為1.60時(shí)坡面高程位移Fig.11 Elevation displacement of slope with reduction factor of 1.60

4 結(jié)論

本文基于連續(xù)-非連續(xù)耦合數(shù)值模擬方法，依據(jù)邊界墻耦合理論建立PFC-FLAC耦合模型，并通過強(qiáng)度折減法對邊坡進(jìn)行了穩(wěn)定性分析，得到結(jié)論如下：

1) 通過該方法驗(yàn)證了連續(xù)-非連續(xù)耦合分析在邊坡工程中的應(yīng)用可行性，運(yùn)用連續(xù)-非連續(xù)方法分析可以更直觀高效地評估邊坡穩(wěn)定性。

2) 通過強(qiáng)度折減法求解了某工程高填方邊坡的安全系數(shù)，并研究了其在降雨工況下和極端工況下的邊坡穩(wěn)定性情況，得出在降雨工況下邊坡較為穩(wěn)定，在極端工況下邊坡產(chǎn)生大范圍滑塌，滑塌區(qū)域沿壩體頂部滑出。

3) 高填方邊坡在極端工況下產(chǎn)生滑坡時(shí)，會產(chǎn)生更大的滑坡區(qū)域，更應(yīng)做好對高填方邊坡的支護(hù)工作。