含褶皺巷道圍巖破裂碎脹大變形機理FDEM 數(shù)值模擬研究

賈永杰

（山西工程職業(yè)學院 資源與安全工程系，山西 太原 030009）

0 引 言

當前對含褶皺巷道軟弱圍巖破裂碎脹大變形機理的研究還存在諸多不足之處，如大都采用連續(xù)性數(shù)值模擬方法，該方法僅能模擬巷道開挖后圍巖的彈塑性連續(xù)變形，無法模擬圍巖破裂過程、破碎塊體間的接觸擠壓和塊體的宏觀大運動過程；此外，未考慮地應力側(cè)壓系數(shù)、褶皺幾何形態(tài)和巖層特性的影響。因此，本文將采用有限元-離散元耦合數(shù)值模擬方法（FDEM） 模擬含褶皺巷道圍巖破裂碎脹大變形機制，并研究地應力側(cè)壓系數(shù)、褶皺幾何形態(tài)和不同巖層特性的影響。相比于傳統(tǒng)的單一連續(xù)性或非連續(xù)性數(shù)值方法，F(xiàn)DEM 能夠模擬完整巖石材料受載后的彈塑性連續(xù)變形及開裂失效過程，也能夠模擬破裂塊體間的接觸擠壓作用及塊體的宏觀大運動過程，且裂隙網(wǎng)絡數(shù)目不受限、計算效率較高，模擬結(jié)果正確性已經(jīng)過室內(nèi)試驗和工程現(xiàn)場監(jiān)測所驗證。

1 FDEM 數(shù)值模擬基本原理

1.1 FDEM 本構(gòu)模型

FDEM 基本原理：將巖石材料劃分為三角形單元和初始無厚度的四邊形節(jié)理單元，三角形單元在外荷載或其他荷載作用下僅發(fā)生線彈性變形，四邊形節(jié)理單元則經(jīng)歷峰前線彈性變形→峰后塑性損傷→斷裂失效全過程，失效后的節(jié)理單元不再進入計算程序并產(chǎn)生一條微裂隙，且兩側(cè)的三角形單元由粘結(jié)關(guān)系轉(zhuǎn)換為接觸關(guān)系。由于FDEM 基本原理能在諸多文獻中輕易找到，本文不再贅述，僅介紹改進后的四邊形節(jié)理單元本構(gòu)模型。

本文對四邊形節(jié)理單元粘結(jié)應力本構(gòu)模型進行了修正，如下：

式中：σn、τ 分別為法向應力和切向應力；o、分別為拉伸/壓縮位移和剪切位移；op、ot、sp和st分別為拉伸峰值位移和拉伸極限位移、剪切峰值位移和剪切極限位移；ot和st的值采用Deng 等[18]提出的相對值形式，即相對于網(wǎng)格尺寸，而非絕對位移形式，如此可消除I 型和II 型斷裂能GI、GII 對網(wǎng)格尺寸的依賴性，即不同的網(wǎng)格尺寸可采用相同的斷裂能；ft、c分別為抗拉強度和粘聚力；z為峰后軟化函數(shù)。

1.2 層狀巖體FDEM 數(shù)值模擬原理

可采用斷續(xù)裂隙網(wǎng)絡法（DNF） 和彌散方法模擬（smeared method），而DFN 難以模擬真實層狀巖體特性[9]，因而在后續(xù)的研究中更多地采用彌散方法[13]，其基本原理概述如下：將層狀巖體的層理面通過直接建模的方式表征（圖1a），平行于層理面的四邊形節(jié)理單元巖石強度最低（包括抗拉強度ft、粘聚力c、I 型斷裂能GI 和II 型斷裂能GII）、垂直于層理面的四邊形節(jié)理單元巖石強度最高（圖1b），與層理面呈其他角度的四邊形節(jié)理單元強度采用線性計算公式：

圖1 層狀巖體FDEM 數(shù)值模擬原理Fig.1 FDEM numerical simulation principle of layered rock mass

式中：ft,γ、cγ、GI,γ、GII,γ為與層理面呈γ 角度的四邊形節(jié)理單元抗拉強度、粘聚力、I 型斷裂能和II型斷裂能；ft,min、cmin、GI,min、GII,min為平行于層理面的四邊形節(jié)理單元抗拉強度、粘聚力、I 型斷裂能和II 型斷裂能；ft,max、cmax、GI,max、GII,max為垂直于層理面的四邊形節(jié)理單元相應參數(shù)。

實際上，彌散方法是一種等效方法，即真實巖體若干層理面內(nèi)的巖石被表征為FDEM 中的兩相鄰層理面內(nèi)的網(wǎng)格，這些若干層內(nèi)的真實巖石宏觀力學特性可通過FDEM 四邊形節(jié)理單元強度參數(shù)的線性變化來反映，避免在FDEM 中劃分非常微小的層理間距和網(wǎng)格尺寸，提高了建模效率和計算效率。

2 含褶皺巖體破裂碎脹大變形機理FDEM數(shù)值模擬

2.1 模型與參數(shù)

平直巖層經(jīng)水平地應力擠壓作用后將發(fā)生彎曲變形，F(xiàn)ereshtenejad 等[1]指出可采用正弦/余弦函數(shù)表示簡單的褶皺構(gòu)造。

為了消除異形斷面的影響，本文以TBM 掘進圓形巷道為例。如圖2（a） 所示，在本文中，將波長恒定為3.6 m，不同形狀的褶皺通過式（4）中的a值體現(xiàn)。將a值分別設定為0（平直層理面）、0.5 m 和1.5 m，建模過程如下：根據(jù)公式（4） 獲取圖2（a） 所示9 個綠色拐點的坐標，在Gmsh 軟件中采用樣條曲線Spline 命令將這9 個點相連獲得一條余弦曲線，同一層面的余弦曲線根據(jù)相同的方法繪制，首尾相連獲得一條褶皺曲線，將該層的褶皺曲線向上、下平移，確保翼部間距t=1.0 m，由此形成整個褶皺構(gòu)造數(shù)值模型。

巷道開挖位置如圖2（a） 所示，圓形巷道直徑3.0 m，模型為邊長80 m 的正方形，以a=1.5 m為例，其模型如圖3（b） 所示，為減小網(wǎng)格數(shù)目、提高計算效率，將模型劃分為3 個區(qū)域[13]：遠場區(qū)、網(wǎng)格細化區(qū)和巷道區(qū)（亦稱為核心材料區(qū)）。褶皺構(gòu)造僅在巷道區(qū)和網(wǎng)格細化區(qū)內(nèi)體現(xiàn)，遠場區(qū)內(nèi)不設置褶皺構(gòu)造，但需確保所有的裂隙網(wǎng)絡在網(wǎng)格細化區(qū)內(nèi)擴展。本文中網(wǎng)格細化區(qū)為邊長約30 m 的正方形，巷道周邊網(wǎng)格尺寸h=0.1 m。巷道開挖模擬分為兩個階段：地應力施加階段和開挖模擬階段。在地應力施加階段，根據(jù)所需地應力計算出所有節(jié)點的相應節(jié)點力，而后將節(jié)點力反向施加至對應節(jié)點上，在該階段不插入四邊形節(jié)理單元，僅采用三角形單元，且模型邊界自由，系統(tǒng)將產(chǎn)生巨大動能，在粘滯阻尼和臨界遲滯阻尼[10]作用下，模型達到平衡，獲得所需地應力；在開挖模擬階段，插入初始零厚度的四邊形節(jié)理單元，固定模型邊界以保持第一階段獲得的地應力；由于四邊形節(jié)理單元的插入，三角形單元間會發(fā)生非常微小的嵌入，模型再次產(chǎn)生較大的動能，但在阻尼作用下很快達到平衡狀態(tài)，隨后可激活巷道開挖模擬程序。

圖3 真實巖層褶皺構(gòu)造及FDEM 數(shù)值模型Fig.3 Real rock fold structure and FDEM numerical model

輸入表1 參數(shù)，對于層狀巖體可直接取表1 中的參數(shù)；對于各向同性巖體，取表1 中垂直層理面的參數(shù)值，該參數(shù)的可靠性已在文獻[13]中得到了驗證。

表1 FDEM 模擬參數(shù)[13]Table 1 Parameters of FDEM simulation

2.2 各向同性圍巖

在各向同性巖體中，各層巖體物理力學特性相同，且層理面力學參數(shù)與完整巖石力學參數(shù)相同。將地應力側(cè)壓系數(shù)λ 分別設定為0.5（水平地應力σh=14 MPa、垂直地應力σv=28 MPa） 和2.0（σh=28 MPa、σv=14 MPa），得到如圖4～圖5 所示的模擬結(jié)果。

根據(jù)模擬結(jié)果可知，在不同地應力側(cè)壓系數(shù)和褶皺形態(tài)下，巷道均因圍巖的破裂碎脹產(chǎn)生嚴重大變形災害，其機理可表述為：隨著巷道核心材料的逐步軟化（在真實巷道掘進中體現(xiàn)為隨著掌子面的推進），巷道表面圍巖產(chǎn)生徑向應力降低、切向應力升高現(xiàn)象，即為卸圍壓、升軸壓的三軸壓縮力學模型，當升高的切向應力超過該圍壓下的巖體強度時，圍巖發(fā)生X 型共軛剪切破裂并伴隨著拉伸斷裂，最大切向集中應力向深處的完整巖體轉(zhuǎn)移，直至開挖完畢后切向應力在圍巖深處與巖體強度達到極限平衡狀態(tài)，裂紋不再往更深處擴展；破裂后的塊體一方面發(fā)生彈性變形恢復，產(chǎn)生體積膨脹，另一方面主剪切帶發(fā)生滑移剪脹效應、塊體間不嚙合產(chǎn)生大量空隙，更重要的是破碎塊體間由于運動速度的不一致性產(chǎn)生相互擠壓導致塊體發(fā)生翻轉(zhuǎn)大運動并向巷道空間內(nèi)移動，使巷道斷面積急劇縮小，上述過程為圍巖體的破裂碎脹大變形災變過程。

以圖4（a） 為例，由于巷道開挖導致徑向應力降低、切向應力升高，進而產(chǎn)生的X 型共軛剪切裂隙與牛雙建[11]和潘一山[3]等人通過室內(nèi)巷道開挖模型試驗得到的圍巖破裂形態(tài)和理論剪切滑移線場（圖6a） 非常相似；此外，巷道表面破碎塊體的翻轉(zhuǎn)大運動形態(tài)與潘一山[12]等人通過室內(nèi)模型試驗結(jié)果非常相似，如圖5（b） 所示。在圍巖位移場方面，本文通過FDEM 數(shù)值模擬得到的圍巖位移場（圖6b） 與洛鋒等人基于室內(nèi)模型試驗和FLAC3D 數(shù)值模擬提出的解析解非常相似（圖6a），亦可分為3 個區(qū)域：第I 破裂帶、第II 破裂帶和第III 破裂帶，且第I 破裂帶呈楔形三角形形態(tài)。

圖6 巷道圍巖X 型共軛剪切滑移解析解及FDEM 數(shù)值模擬結(jié)果Fig.6 Analytical solution of X-type conjugate shear slip of roadway surrounding rock and FDEM numerical simulation results

根據(jù)本節(jié)模擬結(jié)果可知，對于各向同性巖層且層理面參數(shù)與圍巖體一致時，褶皺的存在對圍巖破裂碎脹大變形影響不大，均呈現(xiàn)以X 型共軛剪切滑移破裂為主且伴隨少量拉伸裂隙的對稱性破裂碎脹大變形。當然，對于不同的褶皺形態(tài)，圍巖碎裂形態(tài)和裂隙網(wǎng)絡形態(tài)存在微小差異，這是由于FDEM 數(shù)值模擬對網(wǎng)格的依賴性造成的。

3 結(jié) 語

本文采用有限元- 離散元耦合數(shù)值模擬方法FDEM 揭示了不同地應力側(cè)壓系數(shù)和不同褶皺形態(tài)下的圍巖破裂特征和破裂碎脹大變形機理。在不考慮褶皺形成后殘余構(gòu)造應力的影響下，對于各向同性圍巖體而言，褶皺的影響較微弱。巷道開挖卸荷引起的切向集中應力超過巖體強度時發(fā)生X 型共軛剪切破壞并伴隨少量拉伸裂隙，破碎巖塊沿著主剪切裂隙帶發(fā)生滑移剪脹，此外，大量巖塊在相互接觸擠壓下產(chǎn)生向巷道空間的宏觀大運動（包括平移和翻轉(zhuǎn)大運動） 并產(chǎn)生大量空隙，造成破碎巖體的體積膨脹并引起巷道斷面積的急劇縮小，是破裂碎脹擠壓大變形致災機理。