月牙山墜落式危巖穩(wěn)定性離散元數(shù)值分析

張 研，吳康麗，王 偉，高 宇

(1.桂林理工大學 土木與建筑工程學院， 廣西 桂林 541004；2.廣西巖土力學與工程重點實驗室， 廣西 桂林 541004)

隨著經(jīng)濟的發(fā)展，我國西南山區(qū)建設速度呈逐年增長趨勢，但由于自然地理條件、地質構造和地層巖性等復雜因素的影響，使得地質災害成為威脅西南山區(qū)居民生命財產(chǎn)、城鎮(zhèn)建設和交通運輸安全的首要因素。危巖崩塌屬于三個地質災害種類中重要的一種，其破壞性雖比地震、火山要小，但由于其發(fā)生的頻率和廣度遠超其它災害，故其危害性不容小覷[1]。巖體的穩(wěn)定性常常由巖塊之間的裂隙、節(jié)理、斷層等不連續(xù)結構面所決定?？v橫交錯的結構面將巖體進行切割，形成我們?nèi)粘Ｋ姷牟煌耆B續(xù)體，并可逐漸形成為危巖體[2]。

目前，針對危巖穩(wěn)定性分析及演化所采用的研究方法主要有實驗法、解析法和數(shù)值模擬法。其中，實驗法最為直接，但是具有周期長、成本較高和偶然性大等諸多缺點。又因危巖穩(wěn)定性受多種因素影響，解析法無法將所有影響因素進行系統(tǒng)分析，故其計算結果代表性不佳。隨著計算機技術和理論的發(fā)展，使得數(shù)值模擬技術已經(jīng)廣泛應用于危巖穩(wěn)定性問題研究[3]。目前，主要有有限元法、有限差分法、離散元法等。有限元法和有限差分法主要以連續(xù)介質為基礎，而危巖的崩塌過程大多是不連續(xù)的，這使得在實際工程中無法求解有初始間隙的相互接觸體；并且以單元內(nèi)應力作為屈服依據(jù)，可能會導致在計算位移和應力時出現(xiàn)計算失真[4]。離散法允許巖體之間出現(xiàn)滑動、平移、轉動和巖體斷裂等復雜過程,可用于模擬巖塊產(chǎn)生大變形、節(jié)理面發(fā)生滑動甚至滑落等危巖穩(wěn)定性問題時，不會對模型產(chǎn)生較大影響而產(chǎn)生誤差，此外還具有識別新的接觸面、重構計算模型進行分析等諸多優(yōu)點[5-6]。

本文通過對桂林月牙山進行現(xiàn)場地質勘察，共發(fā)現(xiàn)10處群體危巖，其中W1危巖為墜落式危巖。從巖體規(guī)模和所處位置的高度來看，其體積較山體北側的其他危巖體都要大，一旦發(fā)生危巖崩塌，W1墜落式危巖體所造成的影響將為最大。本文選取W1墜落式危巖為研究對象，運用UDEC離散元軟件對其穩(wěn)定性進行定性評價，揭示了危巖體的變形發(fā)展趨勢，為合理選取工程防治措施提供參考。

1 月牙山區(qū)域地質概況

1.1 工程概況

桂林月牙山位于七星景區(qū)內(nèi)，為典型的喀斯特巖溶地貌，山體主要由石灰?guī)r組成，巖體長期裸露，受降雨等外界環(huán)境侵蝕，巖體表面裂隙發(fā)育豐富(見圖1)。隨著時間推移，風化、水侵蝕等作用使得裂隙逐漸擴大、加深，在危巖體自身重力或外力影響下，隨時會發(fā)生崩塌落石災害。

圖1月牙山巖體裂隙

1.2 危巖分布概況

月牙山主體由硬質的石灰?guī)r巖體、巖塊組成，經(jīng)過初步的野外地質勘察，共發(fā)現(xiàn)危巖體10處(見圖2)，分別有傾倒式危巖、滑移式危巖、墜落式危巖，其中傾倒式危巖2處，滑移式危巖5處，墜落式危巖3處，月牙山危巖體具體分布情況見表1。

圖2月牙山初步勘察危巖分布概況

表1 月牙山野外初步勘察危巖體資料表

2 月牙山危巖成因分析

2.1 地形地貌特征

桂林地區(qū)為典型的喀斯特地貌，受巖溶地質作用十分明顯，山體巖性多為上泥盆統(tǒng)融縣組(D3r)石灰?guī)r，石灰?guī)r屬碳酸鹽巖，堅硬、質純、厚層，由于燕山運動的原因，境內(nèi)山體多高聳直立、陡峭挺拔，山體上植被發(fā)育茂密，亦有石灰?guī)r大量裸露，容易發(fā)生危巖崩塌災害[7-8]。

2.2 水文氣象條件

桂林所處緯度較低，接近北回歸線，屬亞熱帶季風性氣候。全年氣候溫和，降雨量大且集中于夏季，夏長冬短且雨熱相隨，年平均相對濕度為73%～79%。該地區(qū)內(nèi)地表河徑流豐富，集雨面積較大且達到100 km2及以上的河流有六十余條，地下河流的分布亦相當豐富，平均每年的水分蒸發(fā)量約為1 490 mm～1 905 mm。因此，該地區(qū)大部分石灰?guī)r山體長期受到降雨、霧氣和地下水活動的影響作用，產(chǎn)生危巖崩塌的可能性更大。

2.3 人為因素影響

月牙山位于著名的七星景區(qū)內(nèi)，游客數(shù)目巨大，頻繁的人類活動對自然生態(tài)和地理環(huán)境造成一定程度的影響。工業(yè)和城市化進程的不斷發(fā)展也給空氣帶來了越來越多的工業(yè)廢氣和汽車尾氣，多種物質經(jīng)過一系列復雜的化學反應后會產(chǎn)生SO2-和NO3-等酸根離子，這些離子與空氣中的水蒸氣結合后形成酸雨。酸雨的出現(xiàn)，使得巖石中的一些礦物在酸性條件下發(fā)生化學分解，導致原有的巖石結構破壞，從而降低巖石的力學性能，加大危巖崩塌、滑落的風險[9]。

3 離散元法基本原理

目前在工程領域廣為應用的數(shù)值模擬方法主要有有限單元法、有限差分法和離散單元法。有限元法(如ANSYS、ABAQUS)雖然應用最廣，但不適用于大變形、不連續(xù)體等問題的模擬；有限差分法(如FLAC3D)存在計算邊界和單元網(wǎng)格劃分的隨意性；不連續(xù)變形分析法(如DDA)亦有參數(shù)選取的局限性[10-11]。離散元法除具有將巖體視為非連續(xù)結構體、允許發(fā)生接觸位移、模擬過程更加真實等優(yōu)點之外，還包含豐富的材料本構關系和接觸本構關系[12]。是研究和分析非連續(xù)巖體穩(wěn)定性問題的優(yōu)選方法。

基于牛頓第二定律建立起來的離散元法與其他數(shù)值模擬方法相比，不需要復雜的數(shù)學推理過程，只需滿足物理方程和運動方程：

(1) 物理方程。假定塊體之間的法向作用力Fn和法向位移Wn成正比，則表達式為：

Fn=Kn·Wn

(1)

式中：Kn為界面的法向剛度系數(shù)，GPa/m。

塊體與塊體之間剪切力用剪切增量ΔFs表示，則表達式為：

ΔFs=Ks·Ws

(2)

式中：Ks為界面的切向剛度系數(shù)，GPa/m；Ws為塊體之間的相對位移，m。

(2) 運動方程。由于塊體之間是相互接觸的，因此每個塊體上受到不同角度的力。依據(jù)公式分別計算出作用在塊體X方向、Y方向上的合力FX、FY和合力矩M。

(3)

(4)

μ=F合/m

(5)

式中：X0為質點橫坐標；Y0為質心縱坐標；μ為加速度，m/s2；F合為合力，N；m為塊體質量，kg。

美國ITASCA軟件公司開發(fā)的二維離散元法UDEC軟件利用顯式差分法基本原理，可以為不穩(wěn)定物理過程提供穩(wěn)定解和為分析精度提供基本技術保障[13]。應用離散元法模擬分析時，巖體被當成僅由塊體和接觸面組成的模型，巖體整體模型設置有塊體模型和節(jié)理模型兩種，可以根據(jù)實際情形進行相應設置。UDEC離散元法程序軟件中有兩個較為常用的裂隙命令：Jset和Crack。Jset命令可以根據(jù)現(xiàn)場勘查所得的巖體節(jié)理參數(shù)資料賦值建立巖體裂隙，一般用于大量節(jié)理組的建立，并且可根據(jù)實際模型情況設置節(jié)理相關參數(shù)(如起始坐標、傾角、間距、巖橋等)，準確建立節(jié)理在模型中分布的位置和情況。Crack命令通過對任意兩點坐標進行設置，兩點坐標的直線連線即為巖體中的裂隙，這一命令常適用于建立單一節(jié)理面巖體的裂隙[14]。本文UDEC軟件進行危巖體穩(wěn)定性分析，具體實現(xiàn)步驟如圖3所示。

圖3 UDEC模擬計算流程圖

4 工程實例分析

4.1 模型的建立

本文塊體模型選用巖石力學界廣為使用的摩爾-庫侖本構模型[15]，掌握模型巖體實際的密度、體積模量、剪切模量、內(nèi)摩擦角、黏聚力、剪脹角以及抗拉強度，并且根據(jù)這幾個參數(shù)對模型巖體進行賦值建立初步塊體模型。同時采用Jset命令建立模型的節(jié)理。選取W1墜落式危巖體為對象進行研究，危巖算例高約20 m，見圖4，巖體本構模型采用摩爾-庫侖模型，巖體頂部由于當?shù)貧夂蚨嘤甑年P系出現(xiàn)溶痕，漸漸發(fā)展成裂隙、裂縫(圖4中黑粗線段表示裂隙)，巖體右側下部受到差異風化的作用或是下部受到地下水侵蝕的作用逐漸被掏空，形成了此墜落式危巖模型。

圖4山體模型南北向橫截面圖

在模型中，模型底部邊界X軸和Y軸需要固定X、Y軸速度位移都為0，即：Xvel=0，Yvel=0；模型的左側邊界X方向速度需要邊界約束，即：Xvel=0。

4.2 模擬參數(shù)確定

本文結合現(xiàn)場實地勘察和已有的桂林地區(qū)巖層研究成果獲得了相關的巖石力學性質指標，結合地質勘察經(jīng)驗和監(jiān)測分析，依據(jù)文獻合理的選取參數(shù)數(shù)值[7-8,16]，見表2、表3。

表2 UDEC巖體模型石灰?guī)r參數(shù)取值

表3 UDEC巖體模型節(jié)理參數(shù)取值

4.3 數(shù)值模擬結果分析

運算前將UDEC模型初始應力場調整為穩(wěn)定的應力狀態(tài)，山體的橫截面右側和山體底部設置為固定邊界，在自然條件下僅僅受到重力作用，UDEC模型中重力的作用通過設置重力加速度達到，重力加速度數(shù)值取10 m/s2。在初始應力場中，最大主應力與最小主應力隨著深度的變化而變化，在地勢陡峭的地方主應力數(shù)值最大，而與巖體平行的地方主應力數(shù)值最小，幾乎為零?？拷Ｐ投秆碌撞孔畲笾鲬妥钚≈鲬饾u發(fā)生偏轉，由豎直轉為水平或由平行轉為垂直。

模擬運算開始(見圖5、圖6、圖7)，隨著時間的推移，迭代步數(shù)的增加，在受重力的作用下，W1危巖體與山體之間的裂縫、裂隙逐漸變寬。由于此時危巖體與山體間出現(xiàn)摩擦和拉扯等力的作用，從初期應力場圖可以看出應力變大且集中于裂隙下部或是危巖體與山體連接部，由初期位移矢量圖可以看出裂隙、裂縫變寬后危巖體整體位移明顯，巖體上部位移矢量值最大，巖體中部次之，位移矢量值大小由上至下依次減小。再由初期速度矢量圖可以看出危巖體整體下落趨勢趨于明顯，尤其巖體上部加速度數(shù)值最大，巖體中部加速度數(shù)值次之，此時危巖體的潛在威脅變大。

圖5 山體初期應力場圖

圖6 山體初期位移矢量圖

圖7山體初期速度矢量圖

隨著計算的進行(見圖8、圖9、圖10)，迭代步數(shù)逐漸變大，后期山體與危巖體的裂隙、裂縫變得越來越寬，此時巖體之間的摩擦、拉扯現(xiàn)象更為劇烈，力的作用也更為明顯。通過山體后期應力場圖可以看出裂隙下部或者是山體與危巖體連接部的應力增大和集中現(xiàn)象也變得越來越明顯，山體后期的位移矢量圖中可以清晰的看到位移矢量方向較山體初期位移矢量圖趨于向下的態(tài)勢明顯，而且危巖體上、中、下部位移矢量方向凌亂，與之前的初期位移矢量圖相比位移矢量方向不均勻。山體后期的速度矢量圖可以看出危巖體整體下落趨勢已經(jīng)十分明顯，巖體上、中、下部的速度矢量方向凌亂但整體方向向下，危巖體此時的潛在威脅變得更大。

圖8 山體后期應力場圖

圖9山體后期位移矢量圖

圖10山體后期速度矢量圖

根據(jù)極限平衡法進行危巖體穩(wěn)定性計算的結果為：W1巖體的穩(wěn)定系數(shù)為1.46且小于1.5，說明該巖體處于欠穩(wěn)定狀態(tài)。當條件具備時，W1巖體極易發(fā)生墜落式崩塌破壞。由此可見，上述UDEC所模擬巖體變化發(fā)展趨勢結果與該結論較為一致，且與實際相吻合，故可將計算結果應用于該危巖體穩(wěn)定性的判定。同時，針對月牙山的其它危巖體穩(wěn)定性也應及時做出判斷和治理。

5 結 論

本文采用UDEC軟件對月牙山墜落式危巖體(W1)穩(wěn)定性進行數(shù)值分析，得出如下結論：

(1) 巖體的穩(wěn)定性主要取決于結構面，危巖體在重力的影響下會使裂隙、裂縫下部以及山體與危巖體的連接處產(chǎn)生應力增大和應力不斷集中的現(xiàn)象。

(2) W1危巖體已處于欠穩(wěn)定狀態(tài)，受危巖體所處的高度(重力勢能)、裂隙深度的發(fā)展和外界因素等方面影響，隨著時間的推移，該巖體終將會發(fā)生破壞。

(3) 將UDEC離散元法數(shù)值模擬分析軟件應用于危巖崩塌模擬，結合工程實例經(jīng)驗和先前學者的研究，科學合理地選取巖體力學參數(shù)和調整軟件相應設置。通過與極限平衡法計算結果進行對比，表明該方法對危巖穩(wěn)定性的分析結果準確，可為日后危巖的防治提供有價值的參考數(shù)據(jù)。