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

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

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

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

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

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

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

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

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

1.2 強度折減法

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

(1)

(2)

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

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

2 工程應用

2.1 工程概況

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

2.2 計算模型與條件

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

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

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

3 結(jié)果分析

3.1 天然工況

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

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

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

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

3.2 降雨工況

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

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

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

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

圖6 測線及測點設置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)所示，堆積體的位移沿坡面從坡底向坡頂逐漸增加，位移與高程基本呈線性相關，出現(xiàn)這一現(xiàn)象的原因是由于降雨工況下堆積體的材料參數(shù)折減導致的不均勻沉降。由圖8(b) 可見，堆積體的位移主要朝下，在與風化帶接觸界面角度的影響下，堆積體產(chǎn)生少量沿坡面方向的位移，但在壩體的阻擋作用下，未產(chǎn)生明顯的沿坡面方向位移。

3.3 強度折減

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

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

注：圖例中0表示相鄰顆粒處于非粘結(jié)狀態(tài)，1表示接觸非粘結(jié)并且是由于張拉導致破壞，2表示接觸非粘結(jié)并且是由于剪切破導致產(chǎn)生，3表示接觸處于粘結(jié)狀態(tài)且粘結(jié)完好無損。圖10 折減系數(shù)為1.60時顆粒情況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ù)高強度降雨、地震等)下，潛在失穩(wěn)區(qū)域更大，發(fā)生失穩(wěn)時滑面更深，范圍更廣，更應做好邊坡的支擋工作。

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

4 結(jié)論

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

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

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

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