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

劉天成，王學濱，岑子豪，杜 軒，薛承宇

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

0 引言

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

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

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

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

1 原接觸算法簡介

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

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

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

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

2.1 接觸檢測

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

(1)

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

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

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

2.2 接觸力計算

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

(2)

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

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

F=0.5KnLAihn

(3)

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

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

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

2.3 算法檢驗

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

圖1 2塊體接觸過程

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

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

3.1 計算模型及方案

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

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

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

表1 各巖層參數

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

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

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

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

圖3 右巷支護情況

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

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

3.2 方案3的結果及分析

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

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

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

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

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

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

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

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

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

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

對于煤柱，當N=800 000～1 200 000時，支承壓力總體下降，這是由于采空區形成前后，部分巖層垮落導致了煤柱的部分應力轉移至底板。

3.3 θ的影響

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

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

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

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

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

4 結論

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

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

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