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

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

0 引 言

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

1 合成巖體模型

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

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

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

圖1 巖石的微觀結構Fig.1 Rock microstructure

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

τe=c+σetanφ

(1)

其中：σe和τe分別為黏結接觸法向和切向應力；c、φ分別為黏結的黏聚力、內摩擦角。黏結接觸受力大于剪切強度或者抗拉強度時，子接觸破壞形成微裂隙。

在未黏結接觸中，其剪切力服從庫倫摩擦準則：

τe=σetanφr

(2)

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

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

1.2 光滑節(jié)理模型

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

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

1.3 裂隙網絡模型

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

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

x0=lR1

(3)

y0=lR2

(4)

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

(5)

式中：μ為裂隙的平均角度；σ為傾角的標準差。

馮·米塞斯分布如下：

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

(6)

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

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

L=-uln(1-R)

(7)

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

2 合成巖體模型建立

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

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

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

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

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

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

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

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

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

2.3 建立合成巖體模型

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

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

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

3 模擬結果

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

3.1 裂隙密度對巖柱強度的影響

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

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

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

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

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

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

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

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

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

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

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

2)巖柱強度隨裂隙密度P21增加而降低。裂隙密度較小時，巖柱成“X”形壓剪破壞；裂隙密度較大時，巖柱“X”形壓剪破壞逐漸消失；巖柱產生的微裂隙與裂隙密度成負相關。

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