復(fù)雜巖質(zhì)邊坡隨機(jī)節(jié)理模擬及離散元穩(wěn)定性分析

吳宇欣,巫錫勇,凌斯祥,劉禮釗

(西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院,成都 611756)

2008年汶川地震后，震區(qū)斜坡巖土體受到很大擾動，加之地震區(qū)余震頻發(fā)，汶川縣地質(zhì)災(zāi)害更加頻繁且嚴(yán)重。地震導(dǎo)致邊坡巖體節(jié)理裂隙發(fā)育，巖體強(qiáng)度隨之降低，進(jìn)而影響邊坡的穩(wěn)定性，巖質(zhì)邊坡的穩(wěn)定性問題已成為影響和制約公路和鐵路發(fā)展的關(guān)鍵技術(shù)問題之一[1]。汶川地震后震區(qū)邊坡的研究成果豐富，冉濤[2]對川西交通走廊的巖質(zhì)邊坡失穩(wěn)模式及破壞機(jī)理進(jìn)行了分析，總結(jié)出典型的6種斜坡破壞-失穩(wěn)模式；許強(qiáng)[3]對汶川誘發(fā)的大型滑坡崩塌災(zāi)害進(jìn)行了研究，發(fā)現(xiàn)其主要具有震裂潰屈、臨空拋射和碎屑流化等獨特的動力學(xué)特征；牛家永[4]對地震作用下邊坡的動力響應(yīng)規(guī)律進(jìn)行了研究，證明地震波對復(fù)雜結(jié)構(gòu)面邊坡的影響特征。某擬建鐵路穿越汶川地震核心影響區(qū)，在調(diào)查了鐵路沿線隧道進(jìn)出口邊坡后，發(fā)現(xiàn)漁子溪北岸某隧道進(jìn)口巖質(zhì)邊坡由花崗巖組成，其巖體節(jié)理發(fā)育且空間分布復(fù)雜，具有變形破壞風(fēng)險，需對其邊坡穩(wěn)定性進(jìn)行詳細(xì)研究，以保證隧道進(jìn)口安全。

邊坡巖體結(jié)構(gòu)模型的建立很大程度上將影響其數(shù)值模擬計算[5]。陳永明[6]建立貫穿的平行節(jié)理模型，但沒有考慮到巖橋以及節(jié)理接觸關(guān)系對巖體強(qiáng)度的影響；李源亮[7]用Beacher模型對巖質(zhì)節(jié)理邊坡進(jìn)行建模，解決的實際邊坡節(jié)理離散的問題，但不能準(zhǔn)確反映各組節(jié)理模擬程度與分布類型。蒙特卡洛法[8]是結(jié)構(gòu)面網(wǎng)絡(luò)模擬的重要方法，其根據(jù)節(jié)理的幾何參數(shù)分布規(guī)律，用跡長、傾角、間距的概率分布模擬出符合實際情況的復(fù)雜節(jié)理邊坡模型，可根據(jù)統(tǒng)計數(shù)據(jù)模擬出巖體內(nèi)部無法觀測的結(jié)構(gòu)面網(wǎng)絡(luò)。

目前用于研究邊坡穩(wěn)定性的數(shù)值模擬方法主要有流形元法(NMM)[9]、不連續(xù)變形分析法(DDA)[10]、有限元法(FEM)[11-12]、離散元法(DEM)[13-14]，有限元法雖然被廣泛應(yīng)用，但不能考慮塊體完全脫離、旋轉(zhuǎn)等情況。NMM能夠?qū)r體連續(xù)與非連續(xù)變形進(jìn)行統(tǒng)一，但目前缺少實用高效的流形元算法。DDA可以用來模擬巖體的裂隙形變和結(jié)構(gòu)的大位移、大轉(zhuǎn)動，但其缺少部分材料和節(jié)理的本構(gòu)關(guān)系，且DDA 通常假定整個塊體為常應(yīng)變塊體，其應(yīng)力場與坡體實際應(yīng)力場不同。離散單元法常用于不連續(xù)巖體變形與穩(wěn)定性分析，且適合模擬節(jié)理系統(tǒng)或者離散顆粒組合體在準(zhǔn)靜態(tài)或者動態(tài)下的變形過程?？紤]到巖質(zhì)節(jié)理邊坡形變的復(fù)雜性與塊體的離散性，離散元軟件中 UDEC(universal distinct element code)軟件[14]能準(zhǔn)確反映坡體材料性質(zhì)本構(gòu)關(guān)系和節(jié)理、裂隙等性質(zhì)的本構(gòu)關(guān)系，能直接計算和模擬結(jié)構(gòu)材料的具體破壞形式和過程，還能滿足全部的平衡方程和邊界條件，在實際工程中受到廣泛應(yīng)用。

本文的研究對象為地處地震帶的汶川縣某擬建鐵路花崗巖邊坡，為了針對花崗巖受風(fēng)化、卸荷等外力地質(zhì)作用產(chǎn)生產(chǎn)狀不穩(wěn)定、組合形式多變的復(fù)雜節(jié)理，采用蒙特卡洛隨機(jī)模擬理論，通過Matlab數(shù)學(xué)分析軟件，模擬巖體隨機(jī)節(jié)理網(wǎng)絡(luò)，并利用赤平投影和UDEC離散元數(shù)值模擬該邊坡進(jìn)行自然狀態(tài)和地震作用下的穩(wěn)定性及變形破壞模式分析。研究成果對該邊坡的鐵路建設(shè)及安全運(yùn)營提供了可靠建議，蒙特卡洛隨機(jī)節(jié)理網(wǎng)絡(luò)結(jié)合UDEC數(shù)值分析的方法對地震高發(fā)地區(qū)節(jié)理發(fā)育邊坡的穩(wěn)定性和破壞模式研究提供了新思路和參考。

1 基本原理及方法

1.1 赤平投影法

赤平極射投影(簡稱赤平投影)是以球體作為投影工具，用以表示物體的幾何要素或點、直線、平面的空間方向和它們之間的角距關(guān)系的一種平面投影。在巖體工程地質(zhì)力學(xué)研究和實踐中，用于表示巖體的結(jié)構(gòu)面、巖體的滑移方向、滑動力和抗滑力等，赤平投影圖常用于結(jié)構(gòu)面發(fā)育巖體邊坡穩(wěn)定性分析。本文采用蔣爵光[15]提出的赤平投影分析方法分析邊坡穩(wěn)定性，該方法將坡面與結(jié)構(gòu)面的關(guān)系在赤平投影上表現(xiàn)出來，結(jié)構(gòu)面運(yùn)動方向位于坡面投影外的塊體有向臨空面運(yùn)動的趨勢，反之則相對穩(wěn)定。

1.2 巴頓法

為避免庫倫-摩爾準(zhǔn)則在抗剪強(qiáng)度計算上的誤差，綜合現(xiàn)場回彈試驗，基于巴頓法(Barton)[16]對其進(jìn)行簡單的等效處理，得到等效的莫爾-庫侖準(zhǔn)則，最終得到等效參數(shù)C和φ：

(1)

(2)

(3)

則等效粘聚力用下面的公式表示：

Ci=τ-σn·tan(φi)

(4)

式中,JRC為結(jié)構(gòu)面粗糙度系數(shù)(Joint Roughness Coefficient)；JCS為結(jié)構(gòu)面壁面抗壓強(qiáng)度(Joint Compression Strength)；φb為基本摩擦角；σn為結(jié)構(gòu)面上的有效正應(yīng)力。

1.3 蒙特卡洛法

蒙特卡洛(Monte-Carlo)法是常用的統(tǒng)計模擬方法和隨機(jī)抽樣技術(shù)，具有廣泛的應(yīng)用范圍，是目前處理隨機(jī)問題相對精確的方法。其一般思想是：通過一些抽樣試驗的方法，來獲得隨機(jī)變量或者隨機(jī)事件的數(shù)字特征，從而達(dá)到求解問題的目的，其理論基礎(chǔ)為大數(shù)定律，即在隨機(jī)試驗中，隨著試驗次數(shù)的增加，試驗結(jié)果的平均值越接近于某個確定的值的規(guī)律。

2 地質(zhì)特征分析

2.1 地質(zhì)概況

以四川省汶川縣映秀鎮(zhèn)某擬建鐵路隧道進(jìn)口邊坡為研究對象(圖1)，該邊坡高約176 m，坡度為40°～60°。斜坡表面地帶第四系粉質(zhì)黏土、塊石土覆蓋層厚度較大，下部基巖為花崗巖，在擬建隧道洞口附近巖石破碎堆積，洞口所在坡面有散落巖塊，G350國道穿越坡體南部。邊坡為晉寧-澄江期第四期巖漿巖(γο2(4))斜長花崗巖構(gòu)成，處于地震高發(fā)區(qū)，由于長期風(fēng)化、卸荷等外界地質(zhì)應(yīng)力作用，邊坡節(jié)理發(fā)育，節(jié)理產(chǎn)狀不穩(wěn)定，且組合形式多樣，與常見的層狀結(jié)構(gòu)面邊坡在地質(zhì)構(gòu)造方面有很大差異。通過實地調(diào)查與勘探，邊坡由外及里分為強(qiáng)風(fēng)化帶、中風(fēng)化帶、弱風(fēng)化帶。根據(jù)《工程巖體分級標(biāo)準(zhǔn)》[17]，邊坡發(fā)育的結(jié)構(gòu)面以Ⅳ類、Ⅴ類硬性結(jié)構(gòu)面(節(jié)理)為主，邊坡存在3組優(yōu)勢節(jié)理，其幾何特征參數(shù)見表1。

表1 巖體結(jié)構(gòu)面幾何特征參數(shù)

圖1 邊坡巖體滑移

2.2 赤平投影分析

根據(jù)上文蔣爵光的理論分析，由實地調(diào)查邊坡巖體節(jié)理數(shù)據(jù)，制作赤平投影圖(圖2)可知，該邊坡巖體的3組節(jié)理切割巖體形成了沿單滑面滑動的滑塌體1、2、3，沿雙滑面滑動的滑塌體12、13、23，及墜落體G，各個滑塌體滑動方向如圖2所示。分析滑塌體穩(wěn)定性，從赤平投影圖中可以看出：單滑面滑塌體2、3以及雙滑面滑塌體12、13、23在坡面投影之外，單滑面滑塌體2、 3的滑動方向分別為316°和324°，雙滑面滑塌體12、13和23的滑動方向分別為301°、310°和243°，其中滑塌體2、3、12、13滑動方向為近坡向，所以易在邊坡面上形成滑塌體。說明該斜坡易沿坡面傾向向下失穩(wěn)破壞或形變。

圖2 斜坡節(jié)理赤平投影圖

從地震的震裂效應(yīng)可知，由于斜坡的坡頂和坡面存在自由邊界，沒有約束，所以在地震壓縮波和剪切波的作用下，被切割的巖體沿節(jié)理面錯動，陡傾節(jié)理裂隙卸荷拉張，緩傾節(jié)理向臨空面滑移，拉張裂隙和壓剪裂隙擴(kuò)展，兩組節(jié)理逐漸發(fā)展形成失穩(wěn)的鍥形體，組成經(jīng)典滑移-拉裂組合(圖3)?；泼鎺r石破碎程度較大，滑移過程中受到極強(qiáng)的拉張剪切作用，由此可以反演推測斜坡是滑移-拉裂破壞。

圖3 斜坡破壞模式

3 基于蒙特卡洛法的隨機(jī)節(jié)理網(wǎng)絡(luò)模型建立

結(jié)合地質(zhì)調(diào)查數(shù)據(jù)和統(tǒng)計資料分析結(jié)果，在邊坡結(jié)構(gòu)面網(wǎng)絡(luò)模擬中，表面強(qiáng)風(fēng)化層的Ⅳ、Ⅴ級結(jié)構(gòu)面跡長和傾角的幾何參數(shù)分布規(guī)律如圖4～6，經(jīng)過函數(shù)曲線擬合，得出J1、J2和J3的傾角和跡長概率分布類型及其均值和標(biāo)準(zhǔn)差，見表2。

表2 蒙特卡洛模擬巖體結(jié)構(gòu)面幾何特征參數(shù)

圖4 J1傾角與跡長概率分布圖

圖5 J2傾角與跡長概率分布圖

圖6 J3傾角與跡長概率分布圖

本文通過matlab實現(xiàn)基于蒙特卡洛法的復(fù)雜節(jié)理模擬，編程的方法如下：

(1) 確定隨機(jī)結(jié)構(gòu)面模型大小(長度和寬度)，依據(jù)結(jié)構(gòu)面間距均值以及結(jié)構(gòu)面面密度公式計算出該組結(jié)構(gòu)面的條數(shù)(n)。

(2) 假定模型區(qū)域結(jié)構(gòu)面跡線中點服從均勻分布，由蒙特卡洛法在模型區(qū)域內(nèi)均勻抽取該組結(jié)構(gòu)面n條跡線的中心點坐標(biāo)(xc，yc)。

(3) 確定該組結(jié)構(gòu)面傾角、跡長的概率密度分布函數(shù)，同樣由蒙特卡洛法隨機(jī)生成對應(yīng)的傾角(θ)、跡長(l)的具體隨機(jī)數(shù)值。

(4) 根據(jù)中心點坐標(biāo)、傾角和跡長這3個參數(shù)就可以確定一條結(jié)構(gòu)面跡線。重復(fù)以上步驟即可獲得某剖面的二維結(jié)構(gòu)面網(wǎng)絡(luò)圖。

根據(jù)表2數(shù)據(jù)，代入matlab計算得到圖7結(jié)構(gòu)面網(wǎng)絡(luò)圖。

圖7 Matlab隨機(jī)結(jié)構(gòu)面模擬

4 UDEC數(shù)值模擬分析

4.1 計算模型及參數(shù)

將上文通過Matlab創(chuàng)建的隨機(jī)網(wǎng)絡(luò)導(dǎo)入UDEC程序，并根據(jù)邊坡的最易滑動面的高程數(shù)據(jù)，建立隨機(jī)節(jié)理邊坡數(shù)值計算模型。該模型相較于傳統(tǒng)的貫穿節(jié)理模型，能更好地模擬出實際情況下節(jié)理的空間分布特征，使得數(shù)值模擬計算更為精確。UDEC數(shù)值模擬的模型高60 m，選取4個監(jiān)測點A、B、C、D，對計算過程中斜坡的位移和速度進(jìn)行監(jiān)測，邊坡結(jié)構(gòu)面網(wǎng)絡(luò)圖模擬結(jié)果見圖8。

圖8 數(shù)值計算模型及監(jiān)測點分布圖

利用二維離散元程序UDEC對邊坡模型進(jìn)行數(shù)值模擬，自然工況下巖體及結(jié)構(gòu)面均采用彈塑性材料，采用摩爾-庫侖屈服準(zhǔn)則，兩側(cè)邊界水平方向約束，底部邊界豎直方向約束；地震工況下巖體及結(jié)構(gòu)面均采用彈塑性材料，采用塑性模型中的摩爾-庫侖模型，兩側(cè)邊界采用自由邊界，底部邊界采用黏滯(不反射)邊界。

根據(jù)《中國地震動參數(shù)區(qū)劃圖》(GB 18306-2015)，該邊坡處于Ⅷ烈度區(qū)，本模型施加的地震加速度為正弦模擬地震波，最大幅值為0.20 g，在模型底部輸入，按照比例放縮，相當(dāng)于烈度為Ⅷ度的地震作用，作用時間0.1 s。模型計算中的阻尼選取為瑞利阻尼，自然頻率的0.1%。

根據(jù)上文巴頓(Barton)法，取得現(xiàn)場巖體力學(xué)參數(shù)，結(jié)合《巖石力學(xué)參數(shù)手冊》，計算得出修正后的邊坡巖體力學(xué)參數(shù)，見表3，計算得出結(jié)構(gòu)面力學(xué)參數(shù)，見表4。

表3 巖體物理力學(xué)參數(shù)

表4 結(jié)構(gòu)面物理力學(xué)參數(shù)

4.2 計算結(jié)果分析

4.2.1 自然狀態(tài)下邊坡穩(wěn)定性分析

自然狀態(tài)下(圖9)，位移最大處位于斜坡頂部，位移量4.5 mm，沿豎向位移量逐漸減小，在斜坡底部達(dá)到最小值。從圖中可以明顯看出，在部分節(jié)理面出現(xiàn)位移量的較大變化，說明在節(jié)理切割巖體，節(jié)理面較巖體強(qiáng)度降低，易發(fā)生相對位移。從圖(圖10)可以看出，在計算8 000步后，最大不平衡力趨近于0，說明坡體穩(wěn)定，不易發(fā)生破壞或失穩(wěn)。

圖9 自然狀態(tài)邊坡位移云圖

圖10 自然狀態(tài)邊坡最大不平衡力圖

從位移時程曲線(圖11)可以看出，隨著計算步長增大，各監(jiān)測點水平和豎直方向位移增大，最終收斂，說明斜坡是穩(wěn)定的。水平方向位移時程曲線顯示，坡體穩(wěn)定后，坡體淺表的4個監(jiān)測點，從坡角到坡頂水平位移值逐漸減小，而豎直方向的位移值從坡角到坡頂水平位移值逐漸增大，說明越往坡頂豎直方向位移越大水平方向位移越小，而越向坡腳則豎直方向位移越小水平方向位移越大。該坡體頂部以豎向位移為主，坡腳處則以水平方向位移為主，符合典型的邊坡滑移模式特征——坡頂豎向位移為主，坡腳水平位移為主，在坡體內(nèi)部形成潛在的滑動面。

圖11 自然狀態(tài)下邊坡監(jiān)測點的位移變化

從自然狀態(tài)下應(yīng)力云圖(圖12)可以看出，(UDEC中，拉應(yīng)力為正值，壓應(yīng)力為負(fù)值)。壓應(yīng)力最大值出現(xiàn)在斜坡底部，最小值出現(xiàn)在坡頂和坡面，壓應(yīng)力沿豎直方向逐漸增大，在自然狀態(tài)下拉應(yīng)力在斜坡中幾乎沒有出現(xiàn)，邊坡在自然狀態(tài)下整體處于受壓狀態(tài)，斜坡整體穩(wěn)定，破壞的風(fēng)險較小。

圖12 自然狀態(tài)下邊坡應(yīng)力云圖

4.2.2 地震作用下邊坡穩(wěn)定性分析

由地震作用下邊坡位移圖(圖13)可知，斜坡表面發(fā)生明顯形變，斜坡整體沿坡面向下以及臨空面運(yùn)動，與節(jié)理赤平投影分析結(jié)果吻合，即節(jié)理發(fā)育的淺表層巖體向臨空面卸荷滑移。斜坡表面節(jié)理發(fā)育處，巖體沿豎直方向節(jié)理張裂，并沿較水平節(jié)理滑移，最大位移達(dá)3 m。危巖體后緣拉裂，沿下緣滑移，節(jié)理對坡體破壞的影響明顯，與實地調(diào)查結(jié)果(圖7)相似。符合滑移-拉裂示意圖(圖14)，后緣結(jié)構(gòu)面逐漸發(fā)育，直至發(fā)生后緣拉裂破壞，后緣裂隙延伸至下部裂隙，發(fā)生滑移-拉裂破壞。從位移云圖也可以看出，斜坡中部存在潛在的貫通的滑移面，在地震作用下易發(fā)生斜坡整體的滑移破壞。

圖13 地震狀態(tài)下邊坡位移云圖

圖14 滑移-拉裂示意圖

在地震作用后(圖15)，斜坡應(yīng)力發(fā)生明顯變化，坡體中上部出現(xiàn)拉應(yīng)力，尤其坡頂及坡面，使得坡頂和坡面發(fā)生拉裂，在坡角出現(xiàn)應(yīng)力集中，坡腳有向上的位移趨勢，坡體中下部處于受壓狀態(tài)。對比自然狀態(tài)下應(yīng)力云圖(圖12)，斜坡整體應(yīng)力明顯增大，壓應(yīng)力減小，拉應(yīng)力增大，說明地震作用賦予坡體拉應(yīng)力，使斜坡處于不穩(wěn)定狀態(tài)，坡面受拉沿節(jié)理發(fā)生破壞，向臨空面滑移。

圖15 地震作用下斜坡監(jiān)測點的位移變化圖

5 結(jié)論

本文結(jié)合實際調(diào)查資料，利用赤平投影和離散元數(shù)值模擬方法，對汶川縣某擬建鐵路斜坡在自然狀態(tài)和地震作用下的穩(wěn)定性及破壞模式進(jìn)行分析，得出以下結(jié)論：

(1) 基于蒙特卡洛法，通過Matlab模擬隨機(jī)復(fù)雜節(jié)理，能準(zhǔn)確地模擬出節(jié)理產(chǎn)狀在一定概率型下的模型和結(jié)構(gòu)面的復(fù)雜接觸關(guān)系，以及符合自然情況下節(jié)理的分布情況。

(2) 自然狀況下，由赤平投影分析得出坡體有向臨空面卸荷滑移的趨勢，節(jié)理組成典型的滑移-拉裂破壞模式；UDEC離散元數(shù)值模擬分析得出斜坡有向臨空面和向下的微小位移，但斜坡整體相對穩(wěn)定。赤平投影和UDEC離散元數(shù)值模擬都分析得出自然狀況下坡體易向臨空面位移，但UDEC離散元數(shù)值模擬更能定量分析出位移大小，赤平投影則能從地質(zhì)角度定性分析斜坡的運(yùn)動趨勢。

(3) 地震作用下，坡面有較大向臨空面的位移趨勢，與自然狀態(tài)相比壓應(yīng)力明顯減小拉應(yīng)力增大，說明地震波對斜坡有拉應(yīng)力作用，斜坡中部表面節(jié)理發(fā)育處破壞最嚴(yán)重，巖體沿豎向節(jié)理拉裂，并沿近水平向節(jié)理滑移，最終形成斜坡滑移-拉裂破壞。