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

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

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

0 引言

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

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

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

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

1 原接觸算法簡(jiǎn)介

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

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

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

2 基于高斯求積公式的勢(shì)接觸力計(jì)算方法

2.1 接觸檢測(cè)

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

(1)

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

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

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

2.2 接觸力計(jì)算

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

(2)

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

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

F=0.5KnLAihn

(3)

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

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

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

2.3 算法檢驗(yàn)

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

圖1 2塊體接觸過程

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

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

3.1 計(jì)算模型及方案

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

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

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

表1 各巖層參數(shù)

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

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

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

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

圖3 右巷支護(hù)情況

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3 θ的影響

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

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

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

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

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

4 結(jié)論

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

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

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