許永祥，何振江，鄭興博，周昌臺

(1.中煤科工開采研究院有限公司,北京 100013;2.煤炭科學(xué)研究總院 開采研究院,北京 100013; 3.大唐陜西府谷煤電有限責(zé)任公司,陜西 榆林 719000;4.山東新巨龍能源責(zé)任有限公司，山東 菏澤 274918； 5.深圳大學(xué) 深地科學(xué)與綠色能源研究院,廣東 深圳 518060)

0 引 言

地下工程中，巖柱或煤柱作為主要支承結(jié)構(gòu)在地下煤礦開采和金屬礦床開采等工程中廣泛存在。無論是長期永久使用還是短期臨時(shí)使用，巖柱或煤柱都要有足夠的穩(wěn)定性，否則巖柱發(fā)生破壞可能會(huì)造成十分嚴(yán)重的后果[1-2]。地下空間開挖打破原有應(yīng)力平衡狀態(tài)，將原有三向受力狀態(tài)轉(zhuǎn)為二向甚至單向受力狀態(tài)，應(yīng)力轉(zhuǎn)移到相鄰支承結(jié)構(gòu)，導(dǎo)致地下空間結(jié)構(gòu)中局部地方出現(xiàn)應(yīng)力集中、大變形、甚至垮塌現(xiàn)象，圍巖結(jié)構(gòu)穩(wěn)定性顯著降低[3-4]。因此，巖柱強(qiáng)度問題直接關(guān)系到巖柱穩(wěn)定性，關(guān)系到地下空間工程安全。地下空間環(huán)境復(fù)雜，巖柱強(qiáng)度受到多種因素影響，包括巖柱尺寸、巖體巖性、巖體中節(jié)理裂隙分布特征等因素。由于所受影響因素復(fù)雜，目前尚未有確定巖柱強(qiáng)度的統(tǒng)一標(biāo)準(zhǔn)與方法。在工程實(shí)踐中，一般采用經(jīng)驗(yàn)類比法和簡化條件下的理論分析[5-6]。經(jīng)驗(yàn)類比法是以實(shí)際工程為背景，考慮巖柱形狀、尺寸和不連續(xù)等特性，基于巖樣實(shí)驗(yàn)室?guī)r石物理力學(xué)試驗(yàn)估算的強(qiáng)度。常用經(jīng)驗(yàn)類比法包括形狀影響經(jīng)驗(yàn)公式法、尺寸影響經(jīng)驗(yàn)公式法、經(jīng)驗(yàn)巖體強(qiáng)度破壞準(zhǔn)則和結(jié)構(gòu)不連續(xù)經(jīng)驗(yàn)分析。理論分析多以彈塑性理論為基礎(chǔ)，認(rèn)為巖柱從邊緣到核心是由塑性屈服區(qū)向彈性區(qū)過渡的，承載能力從外向內(nèi)逐漸增大。雖然這2種方法為實(shí)際工程解決了很多問題，但均不能詳細(xì)地考慮具體實(shí)際工況，不能很好地預(yù)測巖柱的破壞強(qiáng)度及破壞形式。隨著計(jì)算機(jī)技術(shù)發(fā)展，數(shù)值模擬技術(shù)在解決工程實(shí)際復(fù)雜問題中應(yīng)用越來越廣泛。對于巖柱破壞強(qiáng)度和破壞形式的研究，常用的數(shù)值模擬方法是離散元法，該方法最早由CUNDALL[7]提出：該方法將巖體視為由節(jié)理分割形成的離散塊體組合，并通過接觸相互連接，每個(gè)塊體服從牛頓第二定律。近年來，在傳統(tǒng)離散元基礎(chǔ)上，POTYONDY[8]提出黏結(jié)顆粒模型(Bonded particle model)，成功模擬巖石多種力學(xué)特征，包括變形、開裂、膨脹和聲發(fā)射等。在此基礎(chǔ)之上，合成巖體模型技術(shù)(Synthetic Rock Mass,SRM)[9]被提出并用以構(gòu)造巖體結(jié)構(gòu)模型，解決實(shí)際工程復(fù)雜問題。周喻等[10]基于等效巖體技術(shù)分析了斷續(xù)雙節(jié)理巖石試件破裂機(jī)制，并結(jié)合室內(nèi)試驗(yàn)研究成果，驗(yàn)證了等效巖體技術(shù)在斷續(xù)節(jié)理巖體力學(xué)特性研究中的適宜性和可靠性；朱萬成等[11]從節(jié)理巖體表征單元體的力學(xué)意義出發(fā)，研究了節(jié)理巖體參數(shù)確定方法，得出節(jié)理巖體表征單元尺寸(Representative Elementary Volume,REV)及其影響因素；王培濤等[12]研究了應(yīng)用于顆粒流的節(jié)理邊坡巖體強(qiáng)度參數(shù)獲取方法和掛幫礦回采中巖體變形、失穩(wěn)、破壞過程，得到的巖體破壞模式與現(xiàn)場勘測結(jié)果以及相似試驗(yàn)結(jié)果相似。但以往研究較少涉及巖柱穩(wěn)定性，基于此，針對裂隙分布特征對巖柱穩(wěn)定性影響問題，闡述合成巖體模型、建模方法和微觀參數(shù)的選擇與標(biāo)定，通過模擬結(jié)果分析巖體中裂隙密度和裂隙傾角對巖柱破壞強(qiáng)度及破壞形式的影響，為巖柱穩(wěn)定性分析提供思路和方法。

1 合成巖體模型

合成巖體模型是由IVARS等[9]提出的數(shù)值模擬技術(shù)，主要包含離散元(Discrete Element Model,DEM)和節(jié)理網(wǎng)絡(luò)模型(Discrete Fracture Network,DFN)2種模擬方法。復(fù)雜結(jié)構(gòu)巖體可通過黏結(jié)顆粒模型(使用平直節(jié)理黏結(jié)模型)、節(jié)理模型(采用光滑節(jié)理模型)和節(jié)理網(wǎng)絡(luò)模型(DFN)模擬巖塊、裂隙以及節(jié)理網(wǎng)絡(luò)。

1.1 平直節(jié)理黏結(jié)模型

顆粒流離散元模型由一系列直徑不同的剛性顆粒組成，通過顆粒與顆粒間接觸的相互作用來模擬顆粒流介質(zhì)力學(xué)特性。顆粒間接觸的相互作用服從初始設(shè)置的本構(gòu)模型。目前模擬巖石最常用的接觸類型有3種：接觸黏結(jié)模型(Contact-Bonded Model,CBM)、平行黏結(jié)模型(Parallel-Bonded Model,PBM)和平直節(jié)理黏結(jié)模型(Flat-Joint-Bonded Model,FJBN)。以往研究[8,11]表明前2種接觸模型較難得到合理的壓拉強(qiáng)度比，而平直節(jié)理黏結(jié)模型能夠較好解決這個(gè)問題，除此之外，還能更好地模擬巖石的微觀結(jié)構(gòu)(圖1)。因此，宜采用平直節(jié)理黏結(jié)模型模擬巖塊。

圖1 巖石的微觀結(jié)構(gòu)Fig.1 Rock microstructure

典型的二維平直節(jié)理黏結(jié)接觸模型包含剛性球和接觸(圖2)。不同于接觸黏結(jié)模型和平行黏結(jié)模型，平直節(jié)理黏結(jié)模型的接觸可以分割成均長黏結(jié)的子接觸和未黏結(jié)的子接觸[13-16]。這2種接觸服從不同的本構(gòu)關(guān)系，在黏結(jié)接觸中，其剪切力服從摩爾-庫倫準(zhǔn)則：

τe=c+σetanφ

(1)

其中：σe和τe分別為黏結(jié)接觸法向和切向應(yīng)力；c、φ分別為黏結(jié)的黏聚力、內(nèi)摩擦角。黏結(jié)接觸受力大于剪切強(qiáng)度或者抗拉強(qiáng)度時(shí)，子接觸破壞形成微裂隙。

在未黏結(jié)接觸中，其剪切力服從庫倫摩擦準(zhǔn)則：

τe=σetanφr

(2)

式中：φr為殘余摩擦角。

μ′—顆粒間切向摩擦因數(shù)；σt—黏結(jié)的抗拉強(qiáng)度； gc—顆粒表面最小間距，用于判斷是否接觸圖2 平直節(jié)理黏結(jié)模型Fig.2 Flat joint model

1.2 光滑節(jié)理模型

為模擬裂隙力學(xué)性質(zhì)及其破壞行為，文獻(xiàn)[7]提出了光滑節(jié)理模型。當(dāng)光滑節(jié)理模型嵌入黏結(jié)模型時(shí)，平直節(jié)理黏結(jié)模型將會(huì)被光滑節(jié)理模型代替，周圍被賦予光滑節(jié)理模型的球體的力學(xué)行為不再受顆粒黏結(jié)接觸方向影響。2個(gè)相鄰顆粒不發(fā)生沿顆粒表面的繞行行為，只發(fā)生相對滑動(dòng)，從而模擬裂隙的力學(xué)行為(圖3)。

節(jié)理單位法向向量；節(jié)理單位切向向量；光滑節(jié)理接觸法向量；θp—節(jié)理傾角；兩顆粒平均半徑圖3 光滑節(jié)理模型Fig.3 Smooth joint model

1.3 裂隙網(wǎng)絡(luò)模型

成巖作用和長期地質(zhì)作用下，巖體中通常存在大量隨機(jī)裂隙。裂隙網(wǎng)絡(luò)模擬技術(shù)基于蒙特卡洛方法，認(rèn)為裂隙參數(shù)(長度、傾角、密度、位置等)服從某種隨機(jī)分布的隨機(jī)變量，通過計(jì)算機(jī)模擬生成相關(guān)隨機(jī)數(shù)，進(jìn)而生成相應(yīng)隨機(jī)裂隙，典型的裂隙分布函數(shù)有均勻分布、負(fù)指數(shù)分布、正態(tài)分布、對數(shù)分布[17-18]。這些方法生成隨機(jī)網(wǎng)絡(luò)裂隙一般包括裂隙中心點(diǎn)坐標(biāo)、裂隙傾角以及裂隙跡長的確定。

1)裂隙中心點(diǎn)坐標(biāo)的確定。以往研究表明裂隙中心坐標(biāo)服從均勻分布，即裂隙中心點(diǎn)坐標(biāo)(x0,y0)以相同概率出現(xiàn)在研究區(qū)域內(nèi)。假設(shè)研究區(qū)域?yàn)檫呴L為l的正方形，R1和R2為服從(0,1)均勻分布的隨機(jī)數(shù)，則中心點(diǎn)坐標(biāo)可以通過式(3)、式(4)獲得

x0=lR1

(3)

y0=lR2

(4)

2)裂隙傾角的確定。在二維應(yīng)用過程中，裂隙傾角θ服從可以通過包裝正態(tài)分布(wrapped normal distribution)和馮·米塞斯分布模擬。包裝正態(tài)分布函數(shù)如下：

(5)

式中：μ為裂隙的平均角度；σ為傾角的標(biāo)準(zhǔn)差。

馮·米塞斯分布如下：

f(θ)=ekcos(θ-μ)/[2πI0(k)]

(6)

式中：I0(k)為貝葉斯公式；k為方位角集中程度。

3)裂隙跡長的確定。現(xiàn)場研究表明，裂隙的跡長L服從負(fù)指數(shù)分布和冪律分布。假設(shè)跡長服從負(fù)指數(shù)分布，跡長均長為u，則裂隙跡長見式(7)。

L=-uln(1-R)

(7)

式中：R為服從(0,1)均勻分布的隨機(jī)數(shù)。

2 合成巖體模型建立

離散元模型細(xì)觀參數(shù)與模型宏觀力學(xué)性質(zhì)沒有直接關(guān)系，需通過標(biāo)定過程來確定與宏觀力學(xué)性質(zhì)對應(yīng)的細(xì)觀力學(xué)參數(shù)。此過程通常采用數(shù)值模擬試驗(yàn)(單軸抗壓、單軸抗拉或者劈裂試驗(yàn)等)來建立彼此間關(guān)系。以陜西西王寨煤礦為研究背景，地質(zhì)勘探表明巖體完整性差，內(nèi)、外生裂隙較發(fā)育，有必要研究裂隙對巖柱穩(wěn)定性影響。以灰?guī)r為研究對象，其力學(xué)性質(zhì)見表1。

表1 宏觀力學(xué)參數(shù)Table 1 Macro-mechanical parameters

2.1 平直節(jié)理黏結(jié)模型

采用單軸抗壓和單軸抗拉試驗(yàn)標(biāo)定巖體宏觀力學(xué)參數(shù)(單軸抗壓強(qiáng)度、抗拉強(qiáng)度、楊氏模量和泊松比)，即用數(shù)值模擬得到的宏觀力學(xué)參數(shù)去匹配實(shí)驗(yàn)室所得的相應(yīng)力學(xué)參數(shù)。如圖4所示，試件寬度50 mm，高100 mm。試件上下兩端為剛性墻體，模擬加載壓盤，單軸抗壓加載過程中，底壓盤保持不動(dòng)，而上壓盤以一定速度向下移動(dòng)，直至試件破壞，加載過程記錄試件應(yīng)力-應(yīng)變關(guān)系，以獲取試件力學(xué)參數(shù)(單軸抗壓強(qiáng)度、楊氏模量和泊松比)。單軸抗拉試驗(yàn)中，定義上下兩端的部分球體為端部，在抗拉過程中，端部分別以一定速度向外側(cè)拉，從而獲取抗拉強(qiáng)度。

圖4 巖石單軸壓縮實(shí)驗(yàn)室結(jié)果和數(shù)值模擬結(jié)果Fig.4 Uniaxial compression tests of intact rock from experimental and numerical results

標(biāo)定PFC模型的力學(xué)參數(shù)與試驗(yàn)結(jié)果對比見表2，誤差基本小于1%，標(biāo)定的符合宏觀巖體特性的顆粒流細(xì)觀力學(xué)參數(shù)見表2。另外，實(shí)驗(yàn)室和數(shù)值模擬單軸壓縮應(yīng)力-應(yīng)變曲線對比(圖4)表明：數(shù)值模擬結(jié)果高度吻合實(shí)驗(yàn)室試驗(yàn)結(jié)果，可應(yīng)用于下一步數(shù)值建模分析。

表2 合成巖體模型細(xì)觀力學(xué)參數(shù)Table 2 Micro-parameters of SRM model

2.2 光滑節(jié)理模型參數(shù)

完整巖塊力學(xué)參數(shù)標(biāo)定后，加入光滑節(jié)理模型，可模擬裂隙巖體。采用試錯(cuò)方法，標(biāo)定光滑節(jié)理模型力學(xué)參數(shù)。光滑節(jié)理微觀力學(xué)參數(shù)如下：法向剛度200 GPa/m，切向剛度100 GPa/m，黏聚力0.5 MPa，內(nèi)摩擦角30°。

2.3 建立合成巖體模型

巖體合成模型技術(shù)如圖5所示，將裂隙網(wǎng)絡(luò)模型建立的裂隙網(wǎng)絡(luò)嵌入由平直節(jié)理黏結(jié)模型建立的巖石模型中，并將光滑節(jié)理模型的參數(shù)賦予裂隙，從而建立起巖體模型[19]。

圖5 合成巖體模型Fig.5 Synthetic rock mass model

巖柱模型相對于實(shí)驗(yàn)室?guī)r樣較大，需要顆粒數(shù)目較多，采用連續(xù)-非連續(xù)周期邊界技術(shù)建立平直節(jié)理黏結(jié)模型，以降低模型生成時(shí)間。

3 模擬結(jié)果

采用2 m(直徑)×4 m(高度)巖柱模型研究裂隙密度和裂隙傾角對巖柱強(qiáng)度的影響，模擬過程和結(jié)果分析如下。

3.1 裂隙密度對巖柱強(qiáng)度的影響

為研究裂隙密度對巖柱強(qiáng)度的影響，采取控制變量法，即裂隙的其他參數(shù)保持不變，裂隙位置和裂隙傾角采用均勻分布，裂隙長度采用冪次分布。圖6為裂隙密度P21(單位面積內(nèi)累計(jì)裂縫長度)由1.0 m/m2增加至5.0 m/m2的合成巖體巖柱模型。

圖6 不同裂隙密度的合成巖體巖柱模型Fig.6 Composite rock pillar models with different fracture density

為更好反映裂隙隨機(jī)性和不確定性，裂隙網(wǎng)絡(luò)隨機(jī)因子設(shè)置8次，獲取相同參數(shù)條件下不同裂隙結(jié)構(gòu)，然后將巖柱模型進(jìn)行單軸加載獲取巖柱單軸抗壓強(qiáng)度。不同裂隙密度條件下，單軸抗壓強(qiáng)度和變形模量隨裂隙密度增強(qiáng)而減弱。如圖7所示，由于裂隙結(jié)構(gòu)不確定性和隨機(jī)性，裂隙密度相同的條件下單軸抗壓強(qiáng)度和變形模量在一定區(qū)間內(nèi)波動(dòng)，其平均值與裂隙密度呈負(fù)相關(guān)，因而裂隙的存在可顯著降低巖柱強(qiáng)度。

圖7 裂隙密度P21對巖柱抗壓強(qiáng)度和變形模量的影響Fig.7 Effect of fracture density P21 on compressive strength and deformation modulus of rock pillars

裂隙密度P21由1 m/m2增加至5 m/m2過程中，巖柱平均強(qiáng)度由103.5 MPa 降低至69.5 MPa，巖柱平均強(qiáng)度降低34.0 MPa，巖柱平均變形模量由66.7 GPa降低至45.7 GPa。巖柱破壞形式也受裂隙密度影響，如圖8所示，裂隙密度較小時(shí)，巖柱易形成“X”形壓剪破壞，并且形成較小的破壞塊體；隨裂隙密度增大，“X”形壓剪破壞逐漸消失，并形成較大的破壞塊體。此外，如圖9所示，隨裂隙密度增大，外界作用力產(chǎn)生的微裂隙反而減小。這種現(xiàn)象是由于在裂隙較多的情況下，更容易使本身具有的裂隙相互貫通，形成破壞，而裂隙較少時(shí)，巖柱的巖塊部分更容易發(fā)生破壞，形成小塊體和微裂隙。

圖8 裂隙密度對巖柱破壞形式的影響Fig.8 Effect of fracture density on failure modes of rock pillars

圖9 裂隙密度對微裂隙的影響Fig.9 Effect of fracture density on micro-fractures

3.2 裂隙傾角對巖柱強(qiáng)度的影響

類似的，研究裂隙傾角對巖柱強(qiáng)度影響過程中，其他參數(shù)保持不變，裂隙密度P21設(shè)置為3 m/m2，裂隙位置采用均勻分布，裂隙長度采用冪次分布，裂隙傾角由0°增加至90°。不同角度、相同裂隙網(wǎng)絡(luò)參數(shù)條件下，生成8個(gè)隨機(jī)模型如圖10所示。

裂隙角度小于裂隙摩擦角時(shí)，裂隙的存在對巖柱強(qiáng)度沒有影響，如圖11所示。裂隙角度大于裂隙摩擦角時(shí)，巖柱強(qiáng)度隨傾角增加而降低，裂隙角度為60°時(shí)，達(dá)到最小值，然后隨裂隙角度增加而降低；裂隙角度與加載方向一致時(shí)，巖柱強(qiáng)度又重新達(dá)到最大值；巖柱變形模量隨裂隙傾角增加而逐漸增加。

圖10 裂隙不同傾角的合成巖體巖柱模型Fig.10 Composite rock pillar models with different dip angles

圖11 裂隙傾角對巖柱抗壓強(qiáng)度和變形模量的影響Fig.11 Effect of joint orientation on compressive strength and deformation modulus of rock pillars

巖柱破壞形式也受裂隙密度影響，如圖12所示，裂隙傾角小于臨界角φ時(shí)，巖柱破壞形式成“X”形壓剪破壞，形成較小的破壞塊體；隨裂隙傾角增大至45°和60°，巖柱“X”形破壞逐漸消失，并形成剪切滑移破壞，形成的破壞塊體較大；隨裂隙傾角繼續(xù)增大至大于45°+φ/2時(shí)，巖柱破壞形式又重新形成“X”形破壞，形成較小的破壞塊體。另一方面，如圖13所示，微裂隙數(shù)量也隨裂隙角度成“U”形變化。裂隙傾角大于臨界角度時(shí)，原裂隙端部更容易產(chǎn)生翼型張拉微裂隙，更容易發(fā)生貫通，形成破壞。

圖12 裂隙傾角對巖柱破壞形式的影響Fig.12 Effect of fracture dip on failure mode of rock pillars

圖13 裂隙傾角對微裂隙的影響 Fig.13 Effect of fracture dip on micro-fractures

4 結(jié) 論

1)通過節(jié)理網(wǎng)絡(luò)模型和離散元構(gòu)建合成巖柱模型，研究裂隙參數(shù)對巖柱強(qiáng)度的影響，可為確定巖柱強(qiáng)度提供思路。

2)巖柱強(qiáng)度隨裂隙密度P21增加而降低。裂隙密度較小時(shí)，巖柱成“X”形壓剪破壞；裂隙密度較大時(shí)，巖柱“X”形壓剪破壞逐漸消失；巖柱產(chǎn)生的微裂隙與裂隙密度成負(fù)相關(guān)。

3)巖柱強(qiáng)度和微裂隙隨裂隙傾角變化(0°到90°)而呈“U”形變化，裂隙傾角為0°和90°時(shí)達(dá)到最大值，在60°達(dá)到最小值；巖柱變形模量隨裂隙傾角增加而增加；裂隙傾角小于φ或大于45°+φ/2時(shí)，巖柱易形成“X”形壓剪破壞，而裂隙傾角介于φ與45°+φ/2之間時(shí)，巖柱易形成剪切滑移破壞。