切頂卸壓自成巷的連續(xù)-非連續(xù)方法模擬及切縫傾角的影響*

劉天成，王學(xué)濱，岑子豪，杜 軒，薛承宇

(1.遼寧工程技術(shù)大學(xué) 力學(xué)與工程學(xué)院，遼寧 阜新 123000；2.遼寧工程技術(shù)大學(xué) 計算力學(xué)研究所，遼寧 阜新 123000)

0 引言

沿空留巷技術(shù)[1-3]取消了傳統(tǒng)的區(qū)段煤柱，使每條巷道服務(wù)2個工作面。切頂卸壓自成巷技術(shù)通過在巷道頂板切縫，并利用礦山壓力和巖體的碎脹特性實現(xiàn)無煤柱開采[4-11]，近年來得到了快速發(fā)展，具有煤炭采出率高、巷道掘進率低及作業(yè)安全性高等優(yōu)勢。

目前，科技人員多采用連續(xù)方法和非連續(xù)方法開展切頂卸壓自成巷研究[12-13]。韓昌良等[12]采用FLAC3D分析了不切頂與切頂?shù)牟町悾芯堪l(fā)現(xiàn)，卸壓對最終采動應(yīng)力的改善很小，但能大幅降低巷道變形量；高玉兵等[13]采用UDEC考察切縫傾角的影響，研究發(fā)現(xiàn)，切縫傾角偏小將造成矸石對切頂短臂結(jié)構(gòu)產(chǎn)生向下的摩擦作用，且無法對其產(chǎn)生斜撐作用，而切縫傾角偏大將造成切頂短臂結(jié)構(gòu)的重量偏大。所以，切縫傾角以10°～20°為宜。在UDEC中，塊體尺寸被預(yù)先假定，且某些巖層的塊體尺寸較大，這將導(dǎo)致開裂路徑有限，從而不能很好地描述裂紋的時空分布。

連續(xù)方法雖在一定程度上適于模擬連續(xù)介質(zhì)的變形、破壞問題，但一般不適于模擬非連續(xù)問題。非連續(xù)方法雖然適于模擬非連續(xù)問題，但對于應(yīng)力和應(yīng)變的描述較為粗糙。近些年來，連續(xù)-非連續(xù)方法應(yīng)運而生，具有廣闊的應(yīng)用前景。

以作者自主研發(fā)的適于模擬連續(xù)介質(zhì)向非連續(xù)介質(zhì)轉(zhuǎn)化或非連續(xù)介質(zhì)進一步演化的拉格朗日元與離散元耦合連續(xù)-非連續(xù)方法為基礎(chǔ)[14]，提出基于高斯求積公式的勢接觸力計算方法以提高計算效率，并模擬不同切縫傾角時檸條塔煤礦S1201工作面的切頂卸壓自成巷過程。

1 原接觸算法簡介

為便于下文闡述，作如下定義和約定：若某條棱(單元的邊)被2個單元所擁有，則稱為內(nèi)棱；若僅被1個單元所擁有，則稱為外棱；若某個節(jié)點被某條棱所擁有，則稱該棱為該節(jié)點的關(guān)聯(lián)棱；若某節(jié)點的關(guān)聯(lián)棱均為內(nèi)棱，則稱該節(jié)點為內(nèi)節(jié)點，否則稱為外節(jié)點；若某單元所擁有的節(jié)點均為內(nèi)節(jié)點，則稱該單元為內(nèi)單元，否則稱為外單元；將計算模型所在空間劃分為尺寸一致的正方形網(wǎng)絡(luò)，稱每個網(wǎng)格為盒子。

拉格朗日元和離散元耦合的連續(xù)-非連續(xù)方法包括4個模塊：應(yīng)力-應(yīng)變模塊、開裂模塊、接觸-摩擦模塊和運動模塊。其中，接觸算法[15]是基于勢接觸理論。

在原接觸算法中，采用單元-單元接觸模式。該算法包括接觸檢測和接觸力計算2個部分。在接觸檢測部分，在每個盒子內(nèi)，通過判斷每2個外單元是否存在相交的外棱實現(xiàn)嵌入判斷，需進行大量的矢量計算。在接觸力計算部分，對于每2個相嵌入的外單元，通過計算勢函數(shù)在其互嵌區(qū)域邊界上的積分求得接觸力，并將該接觸力分配至相應(yīng)節(jié)點。其中，勢函數(shù)在積分域上分段線性，為此，需確定各區(qū)段端點位置并計算它們的勢，這導(dǎo)致計算效率低下。

2 基于高斯求積公式的勢接觸力計算方法

2.1 接觸檢測

在本文提出算法中，采用點-單元接觸模式。在每條外棱上布置4個接觸點，其編號i為1～4，其位置矢量ri滿足式(1)：

(1)

式中：r0和r5分別為每條外棱上兩端外節(jié)點的位置矢量；xi為4點高斯-勒讓德求積公式的求積節(jié)點，-x1=x4≈0.861 1，-x2=x3≈0.340 0。

將所有接觸點和外單元加入盒子。對每個盒子內(nèi)的接觸點和外單元兩兩進行嵌入判斷。當(dāng)某一接觸點P嵌入某一外單元γ時，稱P與γ的組合為1個接觸對。

提出算法的優(yōu)勢在于：P只會被加入1個盒子。所以，涉及P的接觸對只生成于該盒子的嵌入判斷，而不會有重復(fù)的接觸對生成。

2.2 接觸力計算

γ上的勢函數(shù)φ(P)與原算法中的一致，如式(2)所示：

(2)

式中：Kn為法向接觸剛度，Pa/m；h(P)為點P到γ邊界最短距離函數(shù)。

對于每個接觸對，P上的接觸力F如式(3)所示：

F=0.5KnLAihn

(3)

式中：L為P所在外棱長度，m；Ai為4點高斯-勒讓德求積公式的求積系數(shù)，A1=A4≈0.347 9，A2=A3=1-A1；h為P到γ邊界的最短距離，m；n為垂直于P所在外棱的單位矢量，方向指向P所在單元內(nèi)部。根據(jù)靜力等效原則，將F分配至P所在外棱兩端外節(jié)點，并將F的反力分配至γ的4個節(jié)點。

提出算法的優(yōu)勢在于：通過采用4點高斯-勒讓德求積公式計算勢函數(shù)在2個外單元互嵌區(qū)域邊界上的積分，提高了計算效率。

應(yīng)當(dāng)指出，提出算法是以4個求積節(jié)點為例進行闡述，易被推廣至求積節(jié)點更多的情形。隨著求積節(jié)點的增多，該算法的精度將逼近原算法的精度。

2.3 算法檢驗

模型由上部正方形塊體和下部矩形塊體2部分構(gòu)成，如圖1所示，二者間距為0。上部正方形塊體邊長為0.04 m，被剖分為20×20個正方形單元。下部矩形塊體長度為0.2 m，高度為0.09 m，被剖分為100×45個正方形單元。正方形塊體上邊界受到0.1 MPa的壓應(yīng)力，矩形塊體下邊界被施加法向約束。無重力作用，計算在平面應(yīng)力、大變形條件下進行。各種計算參數(shù)取值：Kn取1×1010Pa/m，面密度ρ取2 700 kg/m2，摩擦系數(shù)f取0.1，時間步長Δt取7.794 2×10-7s，局部自適應(yīng)阻尼系數(shù)α取0.2，彈性模量E取1 GPa，泊松比μ取0.2。

圖1 2塊體接觸過程

圖1(a)給出了正方形塊體下邊界受到的平均法向接觸力隨時間t的演變規(guī)律；圖1(b)給出了模型平衡后垂直應(yīng)力的分布，正、負分別代表拉應(yīng)力、壓應(yīng)力。由此可發(fā)現(xiàn)：在接觸過程中，正方形塊體下邊界受到的平均法向接觸力存在一定的波動，隨著t的增加，波動幅度逐漸減小，直至保持4 000.34 N不變，這與理論值4 000 N相一致。由此說明本文提出算法具有較高的計算精度。另外，正方形塊體的垂直應(yīng)力分布基本均勻；對于矩形塊體，接觸面附近擠壓嚴重，遠離接觸面的上邊界的垂直應(yīng)力為正，這與常識相符。

3 切頂卸壓自成巷過程模擬

3.1 計算模型及方案

模擬不同切縫傾角θ時檸條塔煤礦S1201工作面的切頂卸壓自成巷過程。采場的具體布置及巷道支護方式等內(nèi)容參見文獻[16]。采用錨桿進行原始支護。采用恒阻錨索進行補強支護，采用液壓支柱進行臨時支護。

模型長度為300 m，高度為80 m，被剖分為600×160個單元，工作面推進方向垂直于紙面。以模型左下角為原點O，水平向右為x軸正方向，豎直向上為y軸正方向，模型的左、右及下端面被施加了法向約束，如圖2所示。該模型共包含7個巖層，其中石英砂巖層和粉砂巖層被剖分為四邊形單元，而其他巖層被剖分為正方形單元，邊長為0.5 m。各巖層參數(shù)見表1，表中數(shù)據(jù)主要取自文獻[16]。

圖2 巷道開挖后、煤層開挖前的計算模型

表1 各巖層參數(shù)

為避免四邊形單元不能進一步開裂可能帶來的高應(yīng)力問題，對于煤層和脫離巖層的單元，采用文獻[17]的基于應(yīng)力球量不變假設(shè)的應(yīng)力跌落方法處理，即當(dāng)其應(yīng)力狀態(tài)滿足莫爾-庫侖準則時，將其應(yīng)力從初始屈服面跌至殘余屈服面。煤層和脫離巖層的單元的應(yīng)力跌落系數(shù)β分別取為0.25，0.99。應(yīng)當(dāng)指出，在計算過程中，不允許煤層單元開裂，即將其視為連續(xù)介質(zhì)，而對于其他單元，在脫離巖層前后，介質(zhì)由連續(xù)介質(zhì)向非連續(xù)介質(zhì)轉(zhuǎn)化，煤層與脫離巖層單元的β取值有所差異。對于脫離巖層的單元，β取值較高是為了避免對采場的過大擾動。

其他參數(shù)：Kn取為2個相接觸單元的彈性模量E均值的10 m-1倍，f取0.3，重力加速度g取9.8 m/s2，α取0.5，2個巖層間界面黏聚力峰值取為其黏聚力c均值的0.77倍。計算在平面應(yīng)變、大變形條件下進行。

計算過程：當(dāng)時步數(shù)目N=0～5 000時，模型逐漸達到靜力平衡。模型的上端面被施加了2 MPa的壓應(yīng)力，以模擬80 m厚的上覆巖層[16]。隨后，開挖2條巷道，左巷為輔運巷，右巷為膠運巷，寬度均為6 m，高度均為4 m；左巷左壁的x坐標為77 m，2條巷道之間留有寬25 m的煤柱。

當(dāng)N=20 000時，對右巷進行支護，如圖3所示。將部分煤巖的強度參數(shù)提高0.5倍，以近似模擬錨桿的支護作用。共有2個錨桿支護區(qū)域，分別為煤幫2 m寬的范圍和頂板2 m高、6.5 m寬的范圍。在巷道頂板上施加2個相對的集中力，以近似模擬單列錨索的支護作用。共有2列錨索，錨固力分別為0.087 5 MN(左)和0.35 MN(右)。在巷道頂、底板表面的局部施加相對的應(yīng)力，其與液壓支柱的變形量(即頂、底板高度差的變化量)呈線性關(guān)系，系數(shù)為2 GPa/m，以近似模擬單列液壓支柱的支護作用。共有2列液壓支柱，初始值(初撐力)分別為0.206 MPa(左)和2.52 MPa(右)，最大值(工作阻力)分別為0.7 MPa(左)和4.44 MPa(右)。應(yīng)當(dāng)指出，為了直觀研究左巷的破壞，不對左巷進行支護。

圖3 右巷支護情況

當(dāng)N=35 000時，在右巷的右上角切縫。切縫高度取為9 m。設(shè)計5個計算方案：θ分別取0°,5°,10°,15°，20°。

當(dāng)N=50 000～950 000時，以逐漸卸壓的方式開挖右巷右側(cè)的煤層，共計186 m。

3.2 方案3的結(jié)果及分析

θ=10°時(方案3)最大主應(yīng)力σ3及裂紋區(qū)段的時空分布如圖4所示，正、負分別代表拉應(yīng)力、壓應(yīng)力，灰色和黑色線段分別代表拉裂紋區(qū)段和剪裂紋區(qū)段[18]。

由圖4(a)可知，當(dāng)N=40 000時，巷道的開挖、支護及切縫均已完成，回采尚未開始。兩巷的頂、底板均出現(xiàn)了垂直拉裂紋，相比之下，右巷頂板的裂紋更少，這是由于對右巷進行了支護；切縫尖端出現(xiàn)了σ3集中，且由此發(fā)育出了少量的拉裂紋；煤體發(fā)生了局部破壞，煤體的應(yīng)力分布不連續(xù)。

由圖4(b)～圖4(d)可知，當(dāng)N=120 000～800 000時，右巷右側(cè)的煤層的頂、底板處于逐漸卸壓狀態(tài)。若干巖層出現(xiàn)了越來越多的拉裂紋和剪裂紋，并向下運動；切縫尖端發(fā)展出了較長的低角度拉裂紋，并逐漸向右上方發(fā)展。應(yīng)當(dāng)指出，當(dāng)N=240 000時，左巷直接頂拉裂紋附近已發(fā)展出了一定數(shù)量的剪裂紋。

由圖4(e)～圖4(f)可知，當(dāng)N=1 000 000～1 200 000時，采空區(qū)已形成。由切縫尖端發(fā)展出的低角度拉裂紋促進了采空區(qū)上方部分巖層冒落，模型上表面出現(xiàn)了明顯下沉；左巷頂板發(fā)生了少量冒落。

方案3的煤體支承壓力-x曲線如圖5所示。支承壓力是指初始中心縱坐標為16.75 m的一行煤層單元的垂直應(yīng)力的絕對值。由圖5可以看出，左側(cè)煤層和煤柱的支承壓力的峰值不位于二者的邊緣，這意味著它們均發(fā)生了破壞。應(yīng)當(dāng)指出，煤柱的支承壓力分布呈單峰特性，這意味著煤柱已完全破壞。

圖5 煤體的支承壓力-x曲線(方案3)

對于左側(cè)煤層，當(dāng)N=120 000～240 000時，支承壓力的峰值左移，這是由于向煤層轉(zhuǎn)移的巖層壓力使其破壞區(qū)尺寸增大；支承壓力的峰值增加，這是由于越接近模型左邊界σ3(圍壓)越大。

對于左側(cè)煤層，當(dāng)N=240 000～1 200 000時，支承壓力的峰值減小，這是由于左巷上方逐漸擴展的裂紋阻隔了巖層壓力向左傳遞(圖4(c))。另外，支承壓力的峰值左移。眾所周知，煤層右部巖層壓力的降低將引起其與頂、底板摩擦力的降低，這將導(dǎo)致煤層的圍壓有所釋放，進而導(dǎo)致支承壓力的降低。若出現(xiàn)煤層的支承壓力不足以抵抗巖層壓力的情況，則支承壓力峰值會左移。

圖4 σ3及裂紋區(qū)段的時空分布(方案3)

對于煤柱，當(dāng)N=120 000～800 000時，支承壓力總體上升，這可能是由于對煤層單元進行了應(yīng)力跌落，σ3(圍壓)有所提高。

對于煤柱，當(dāng)N=800 000～1 200 000時，支承壓力總體下降，這是由于采空區(qū)形成前后，部分巖層垮落導(dǎo)致了煤柱的部分應(yīng)力轉(zhuǎn)移至底板。

3.3 θ的影響

N=1 200 000時方案1～2和方案4～5的σ3的分布如圖6所示，具體含義同第3.2節(jié)。由圖6可以看出，θ越大，右巷頂板的完整性越好，這是由于當(dāng)θ較大時切縫尖端發(fā)展出的裂紋促進了采空區(qū)頂板冒落，從而阻隔了巖層壓力向右巷頂板傳遞；θ越大，左巷頂板的冒落越嚴重，甚至，在方案5中，左巷已被冒落的巖塊填滿；θ=0～5°時，左巷上方裂隙較發(fā)育，這說明巖層壓力向左傳遞未被較好地阻隔住。綜合考慮，θ以10°為宜。

圖6 方案1～2和方案4～5的σ3及裂紋區(qū)段的分布(N=1 200 000)

N=1 200 000時方案1～5的煤體支承壓力-x曲線如圖7所示。由圖7可以看出，當(dāng)θ=0～10°時(方案1～3)，支承壓力-x曲線基本一致；當(dāng)θ=15～20°時(方案4～5)，煤柱中部的支承壓力-x曲線呈上凹，且θ越大，煤柱中部的支承壓力越低，煤柱兩幫和左巷左幫的支承壓力越高。

圖7 方案1～5的煤體支承壓力-x曲線(N=1 200 000)

當(dāng)θ較大時，右巷頂板懸露的巖層尺寸較大，使煤柱正上方的直接頂呈上凸的拱形，這會使本應(yīng)向煤柱中部傳遞的較大巖層壓力向別處轉(zhuǎn)移。自然地，左拱角對左巷的直接頂會產(chǎn)生額外的水平方向巖層壓力，這將使該處發(fā)生大規(guī)模剪裂，進而嚴重冒落。與此同時，左巷中冒落的巖塊會對左巷兩幫產(chǎn)生一定的圍壓，這會提升兩幫的支承壓力。

4 結(jié)論

1)以自主研發(fā)的適于模擬連續(xù)介質(zhì)向非連續(xù)介質(zhì)轉(zhuǎn)化或非連續(xù)介質(zhì)進一步演化的拉格朗日元與離散元耦合連續(xù)-非連續(xù)方法為基礎(chǔ)，提出基于高斯求積公式的勢接觸力計算方法，以提高計算效率。

2)利用改進后的連續(xù)-非連續(xù)方法較好地呈現(xiàn)不同切縫傾角時檸條塔煤礦S1201工作面的切頂卸壓自成巷過程，其中，考慮了右巷(膠運巷)的錨桿原始支護、錨索補強支護和液壓元件的臨時支護。具體結(jié)果包括剪裂紋和拉裂紋的時空分布和煤層支承壓力的演化過程。

3)切縫傾角越大，右巷頂板的完整性越好，這是由于當(dāng)傾角較大時切縫尖端發(fā)展出的裂紋促進了采空區(qū)頂板冒落，從而阻隔了巖層壓力向右巷頂板傳遞；左巷頂板的冒落越嚴重。切縫傾角以10°為宜。