賈永杰

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

0 引 言

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

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

1.1 FDEM 本構(gòu)模型

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

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

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

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

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

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

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

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

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

2.1 模型與參數(shù)

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

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

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

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

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

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

2.2 各向同性圍巖

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

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

以圖4（a） 為例，由于巷道開(kāi)挖導(dǎo)致徑向應(yīng)力降低、切向應(yīng)力升高，進(jìn)而產(chǎn)生的X 型共軛剪切裂隙與牛雙建[11]和潘一山[3]等人通過(guò)室內(nèi)巷道開(kāi)挖模型試驗(yàn)得到的圍巖破裂形態(tài)和理論剪切滑移線場(chǎng)（圖6a） 非常相似；此外，巷道表面破碎塊體的翻轉(zhuǎn)大運(yùn)動(dòng)形態(tài)與潘一山[12]等人通過(guò)室內(nèi)模型試驗(yàn)結(jié)果非常相似，如圖5（b） 所示。在圍巖位移場(chǎng)方面，本文通過(guò)FDEM 數(shù)值模擬得到的圍巖位移場(chǎng)（圖6b） 與洛鋒等人基于室內(nèi)模型試驗(yàn)和FLAC3D 數(shù)值模擬提出的解析解非常相似（圖6a），亦可分為3 個(gè)區(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é)果可知，對(duì)于各向同性巖層且層理面參數(shù)與圍巖體一致時(shí)，褶皺的存在對(duì)圍巖破裂碎脹大變形影響不大，均呈現(xiàn)以X 型共軛剪切滑移破裂為主且伴隨少量拉伸裂隙的對(duì)稱性破裂碎脹大變形。當(dāng)然，對(duì)于不同的褶皺形態(tài)，圍巖碎裂形態(tài)和裂隙網(wǎng)絡(luò)形態(tài)存在微小差異，這是由于FDEM 數(shù)值模擬對(duì)網(wǎng)格的依賴性造成的。

3 結(jié) 語(yǔ)

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