王學(xué)濱，薛承宇，岑子豪

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

采場(chǎng)上覆巖層的開裂和破斷在力學(xué)上可歸結(jié)為連續(xù)介質(zhì)向非連續(xù)介質(zhì)轉(zhuǎn)化或非連續(xù)介質(zhì)進(jìn)一步演化。從力學(xué)角度研究此規(guī)模較大問題(例如，幾公里范圍的采場(chǎng))的難度較大，其中涉及連續(xù)介質(zhì)力學(xué)、接觸力學(xué)、非線性斷裂力學(xué)和計(jì)算力學(xué)等分支。正確模擬采場(chǎng)上覆巖層運(yùn)動(dòng)過程對(duì)于巖層穩(wěn)定性控制和有關(guān)災(zāi)害預(yù)防具有極其重要的意義。目前，尚缺乏有效的求解手段。

顯然，欲模擬連續(xù)介質(zhì)向非連續(xù)介質(zhì)轉(zhuǎn)化或非連續(xù)介質(zhì)進(jìn)一步演化，有限元方法和有限差分方法等為代表的連續(xù)方法無能為力。在以離散元方法和非連續(xù)變形分析方法等為代表的非連續(xù)方法中，即使在彈性階段，通常都需要接觸面剛度系數(shù)等參數(shù)；巖層被預(yù)先根據(jù)經(jīng)驗(yàn)劃分為不同尺寸的塊體，這與實(shí)際情況有一定差距。兼具連續(xù)方法和非連續(xù)方法優(yōu)勢(shì)的連續(xù)-非連續(xù)方法被認(rèn)為在模擬連續(xù)介質(zhì)向非連續(xù)介質(zhì)轉(zhuǎn)化和非連續(xù)介質(zhì)進(jìn)一步演化方面具有突出的能力，正在快速發(fā)展。目前，已應(yīng)運(yùn)而生了許多著名的程序或軟件，例如有限元與離散元耦合方法[1-2]和拉格朗日元與離散元耦合方法[3]。目前，采用拉格朗日元與離散元耦合方法已經(jīng)開展了一系列研究[4-9]，效果良好。

本研究介紹了自主開發(fā)的拉格朗日元與離散元耦合方法的基本原理，以及以此為基礎(chǔ)最新發(fā)展的水平巖層運(yùn)動(dòng)計(jì)算方法，并通過3個(gè)采場(chǎng)算例，展現(xiàn)了該方法的優(yōu)勢(shì)。

1 拉格朗日元與離散元耦合方法的基本原理

拉格朗日元與離散元耦合方法本質(zhì)上是連續(xù)方法中拉格朗日元方法與非連續(xù)方法中離散元方法的耦合方法。拉格朗日元方法被用于求解彈性體的應(yīng)力和應(yīng)變問題，離散元方法被用于求解接觸和摩擦問題。同時(shí)，還需引入強(qiáng)度理論和虛擬裂紋模型以處理開裂和軟化問題(軟化程度由斷裂能控制)。

該方法的適用性包括兩方面：其一，連續(xù)介質(zhì)向非連續(xù)介質(zhì)轉(zhuǎn)化研究，例如，完整的連續(xù)介質(zhì)由于應(yīng)力超標(biāo)向破裂、分離的非連續(xù)介質(zhì)轉(zhuǎn)化模擬；其二，非連續(xù)介質(zhì)的進(jìn)一步演化研究，例如，若干疊合在一起的巖層(通過嵌入發(fā)生相互作用)運(yùn)動(dòng)模擬。本研究的算例即是如此。

該方法由4個(gè)模塊組成，分別是應(yīng)力-應(yīng)變模塊、開裂模塊、接觸-摩擦模塊和運(yùn)動(dòng)模塊。該方法流程圖如圖1所示。應(yīng)當(dāng)指出，在1個(gè)時(shí)步內(nèi)，上述4個(gè)模塊依次被執(zhí)行1遍，循環(huán)的次數(shù)被稱為時(shí)步數(shù)目N。當(dāng)N達(dá)到預(yù)設(shè)的結(jié)束時(shí)步數(shù)目時(shí)，循環(huán)結(jié)束，輸出單元和節(jié)點(diǎn)等信息。

圖1 拉格朗日元與離散元耦合方法的流程圖

應(yīng)力-應(yīng)變模塊被用于計(jì)算單元的應(yīng)力和應(yīng)變，其具體原理同F(xiàn)LAC(連續(xù)介質(zhì)快速拉格朗日分析)。在彈性階段，對(duì)于連續(xù)介質(zhì)，只引入兩個(gè)彈性參數(shù)，無需引入影響應(yīng)力、應(yīng)變的有關(guān)剛度參數(shù)，因而，介質(zhì)是真正的連續(xù)介質(zhì)。

根據(jù)增量形式的線彈性本構(gòu)方程，由子單元的應(yīng)變?cè)隽喀ij(當(dāng)i=j時(shí)，表示為Δekk)求解子單元的應(yīng)力增量Δσij：

(1)

式中：K為體積模量，Pa；G為剪切模量，Pa；δij為Kronecker符號(hào)，即

(2)

開裂模塊被用于處理介質(zhì)的拉裂和剪裂。利用最大主應(yīng)力準(zhǔn)則判斷介質(zhì)是否發(fā)生拉裂，利用莫爾-庫(kù)侖準(zhǔn)則判斷介質(zhì)是否發(fā)生剪裂。同時(shí)引入I型和II型斷裂能。在必要時(shí)，單元可以沿對(duì)角線開裂。

介質(zhì)發(fā)生拉裂的條件為：

ft=σt-σ3<0，

(3)

式中：ft是拉裂屈服函數(shù)；σt為單軸抗拉強(qiáng)度；σ3為最大主應(yīng)力。

介質(zhì)發(fā)生剪裂的條件為：

(4)

(5)

式中：fs是剪裂屈服函數(shù)；σ1為最小主應(yīng)力；c為黏聚力；φ為內(nèi)摩擦角。

接觸-摩擦模塊被用于處理單元嵌入后的相互作用。引入勢(shì)的概念計(jì)算接觸力[10-11]，以避免法向接觸力與接觸面積無關(guān)的弊端，并使針對(duì)不同的接觸類型的計(jì)算具有統(tǒng)一性。

勢(shì)函數(shù)

(6)

式中：Kn為法向剛度系數(shù)，Pa/m；hmin(Q)為單元內(nèi)任一點(diǎn)Q距單元邊界的最短距離，m；Ω為互嵌區(qū)域。

運(yùn)動(dòng)模塊被用于計(jì)算節(jié)點(diǎn)的速度和位移。彈性力由單元應(yīng)力計(jì)算得到；不平衡力是彈性力和外力(例如，自重和外部載荷)等力的合力；阻尼力與不平衡力呈線性關(guān)系，其比例系數(shù)為局部自適應(yīng)阻尼系數(shù)。計(jì)算模式包括動(dòng)力計(jì)算模式和準(zhǔn)靜力計(jì)算模式(模擬時(shí)間較長(zhǎng)的開采過程)。在連續(xù)介質(zhì)向非連續(xù)介質(zhì)轉(zhuǎn)化或非連續(xù)介質(zhì)進(jìn)一步演化過程中，一些不能再開裂的煤層單元會(huì)與周圍單元發(fā)生相互作用，能承受超過自身承載能力的高應(yīng)力。為此，提出了針對(duì)煤層離散單元聚集體的應(yīng)力跌落方法，或者將煤層視為連續(xù)介質(zhì)，僅對(duì)煤層單元在必要時(shí)進(jìn)行應(yīng)力跌落。

2 基于點(diǎn)-單元接觸模式的水平巖層運(yùn)動(dòng)計(jì)算方法

拉格朗日元與離散元耦合方法僅適于模擬幾萬個(gè)單元數(shù)的采場(chǎng)模型的巖層運(yùn)動(dòng)過程，并不適于模擬10萬個(gè)單元以上的情況，主要是受到微機(jī)計(jì)算速度的限制。例如，當(dāng)巖層較為破碎時(shí)，由于存在大量的接觸檢測(cè)等原因，原接觸-摩擦模塊的用時(shí)較長(zhǎng)。為此，需要研究水平巖層運(yùn)動(dòng)的并行計(jì)算方法，例如，借助圖形處理器(GPU)進(jìn)行加速。以接觸-摩擦模塊中的點(diǎn)-單元接觸模式為例進(jìn)行介紹。

如圖2所示，采用點(diǎn)-單元接觸模式取代單元-單元接觸模式，以起到精簡(jiǎn)算法和適于GPU并行的目的。點(diǎn)即接觸點(diǎn)，布置在單元的邊上。在建模完成且開始計(jì)算前，建立盒子網(wǎng)絡(luò)覆蓋模型。盒子為寬度一致的正方形。在1個(gè)時(shí)步內(nèi)，首先，將各接觸點(diǎn)及與其可能發(fā)生接觸的單元加入盒子；然后，對(duì)于每個(gè)盒子，將其中任一接觸點(diǎn)和任一單元作為1組，對(duì)所有組合進(jìn)行嵌入判斷(每當(dāng)有1組被判定為發(fā)生嵌入，則生成1個(gè)包含該組接觸點(diǎn)和單元的接觸對(duì))；最后，計(jì)算接觸對(duì)中接觸點(diǎn)的接觸力，將接觸力及其反作用力進(jìn)行分配，令接觸對(duì)消亡。應(yīng)當(dāng)指出，在巖層運(yùn)動(dòng)過程中，當(dāng)某一接觸點(diǎn)與任一單元發(fā)生嵌入時(shí)，則1組包含該接觸點(diǎn)和該單元的接觸對(duì)生成；否則，接觸對(duì)未生成。盒子寬度影響接觸對(duì)數(shù)量和計(jì)算效率，盒子寬度過小會(huì)使部分接觸無法被檢索到，盒子寬度過大會(huì)增加不必要的計(jì)算量?？紤]到單元會(huì)變形，根據(jù)經(jīng)驗(yàn)選取盒子寬度。在接觸-摩擦模塊中，由設(shè)備端計(jì)算得出的單元和節(jié)點(diǎn)信息，進(jìn)行接觸判斷，生成接觸對(duì)并計(jì)算接觸力，再將計(jì)算出的接觸力分配給相應(yīng)節(jié)點(diǎn)。給每個(gè)盒子分配1個(gè)線程，各線程同時(shí)計(jì)算，大大提高了計(jì)算效率。

圖2 接觸單元上一條邊的接觸點(diǎn)

對(duì)于點(diǎn)-單元接觸模式，當(dāng)1個(gè)接觸點(diǎn)嵌入1個(gè)單元時(shí)，該接觸點(diǎn)上的接觸力

F=0.5KnLAhn。

(7)

式中：L為該接觸點(diǎn)所在邊的長(zhǎng)度；A為該接觸點(diǎn)對(duì)應(yīng)的求積系數(shù)；h為該接觸點(diǎn)到被嵌入單元邊界的最短距離；n為垂直于該接觸點(diǎn)所在邊的單位矢量，方向指向該接觸點(diǎn)所在單元內(nèi)部。將F分配至該接觸點(diǎn)所在邊兩端節(jié)點(diǎn)，并將F的反力分配至被嵌入單元的4個(gè)節(jié)點(diǎn)。

顯然，新接觸力的計(jì)算比原接觸力的積分計(jì)算簡(jiǎn)單。而且，當(dāng)邊上的接觸點(diǎn)數(shù)量足夠多時(shí)，兩種方法將具有相同的精度。目前，在1條邊上布置了4個(gè)接觸點(diǎn)，接觸點(diǎn)的坐標(biāo)可由其所在邊兩端節(jié)點(diǎn)的坐標(biāo)求得。

在水平巖層運(yùn)動(dòng)計(jì)算方法中，要求各巖層的厚度不變。在模型建立過程中，由下至上依次建立各巖層。通常，在模型的左、右兩側(cè)施加水平約束，在模型的下端面施加垂直約束，在模型的上端面施加垂直向下的均布載荷p。計(jì)算條件為平面應(yīng)變、大變形。對(duì)于各巖層，主要的力學(xué)參數(shù)包括5個(gè)：彈性模量、泊松比、黏聚力、內(nèi)摩擦角和抗拉強(qiáng)度。上述主要參數(shù)均可通過實(shí)驗(yàn)測(cè)得。對(duì)于煤層，還需要應(yīng)力跌落系數(shù)。單元尺寸過大會(huì)導(dǎo)致計(jì)算結(jié)果精度較低，即模擬效果不好；單元尺寸過小會(huì)導(dǎo)致計(jì)算效率較低。根據(jù)經(jīng)驗(yàn)，將單元邊長(zhǎng)取為0.5 m較為適合。這樣，不需要引入I型和II型斷裂能。應(yīng)當(dāng)指出，各巖層的相互作用通過相互嵌入實(shí)現(xiàn)，無需給任意兩個(gè)相接觸的巖層提供接觸剛度。Kn為材料屬性，而非界面屬性，其大小與彈性模量成正比。比例系數(shù)的提高有利于減小嵌入量造成的誤差。

3 巖層運(yùn)動(dòng)過程模擬

本研究模擬了3個(gè)采場(chǎng)模型的巖層運(yùn)動(dòng)，示意圖見圖3。模型1是新巨龍礦1302N工作面走向模型[15]，采用逐漸刪除煤層單元的方式模擬煤層開采，推進(jìn)方向位于紙面內(nèi)且向右，每隔10 000個(gè)時(shí)步開采5 m。模型2是檸條塔煤礦S1201工作面?zhèn)认蚰Ｐ蚚16]，當(dāng)N=50 000～950 000時(shí)，對(duì)右巷右方待采煤層的頂、底板進(jìn)行卸壓，卸壓完畢相當(dāng)于煤層開采完畢，左巷和右巷的高度等于煤層厚度，對(duì)右巷頂板進(jìn)行了切縫，切縫角度為5°，同時(shí)，對(duì)右巷圍巖進(jìn)行等效支護(hù)(包括錨桿、錨索和液壓支柱)，推進(jìn)方向垂直于紙面。模型3是柳巷煤礦30105工作面?zhèn)认蚰Ｐ蚚17]，一次性開挖煤層，煤層中的左巷和右巷的高度小于煤層厚度，推進(jìn)方向垂直于紙面。限于篇幅，僅在表1中給出了部分所需參數(shù)。

表1 3個(gè)采場(chǎng)模型的部分參數(shù)

圖3 3個(gè)采場(chǎng)模型示意圖

3.1 傳統(tǒng)工法條件下走向巖層運(yùn)動(dòng)模擬

當(dāng)N=80 000時(shí)(工作面推進(jìn)距離D=30 m)(圖4(a))，采空區(qū)上方的粉砂巖層3存在一個(gè)正V字形σ3高值區(qū)，采空區(qū)下方存在一個(gè)倒V字形σ3高值區(qū)，兩個(gè)高值區(qū)均為拉應(yīng)力；開切眼后方與工作面煤壁前方的煤層分別存在σ3較高區(qū)域，均為壓應(yīng)力，擠壓程度較輕，前者寬度約為28 m，后者寬度約為15 m；關(guān)鍵層Ⅰ、Ⅱ之間的各巖層的σ3分布不均勻(例如，關(guān)鍵層Ⅱ下方的細(xì)粒砂巖層5的中性層上方的σ3較低，擠壓較為嚴(yán)重)；關(guān)鍵層Ⅱ、Ⅲ的σ3分布較為均勻。

當(dāng)N=160 000時(shí)(D=70m)(圖4(b))，開切眼后方和工作面煤壁前方的煤層的σ3高值區(qū)的寬度分別為30和12 m；頂板存在大V字形的σ3高值區(qū)，為拉應(yīng)力。開切眼附近出現(xiàn)1條向右上方發(fā)展的傾斜拉、剪復(fù)合裂紋(圖中灰色、黑色線段分別代表拉裂紋、剪裂紋)，前期以剪為主，后期以拉為主，已貫穿關(guān)鍵層Ⅰ和細(xì)粒砂巖層2；采空區(qū)上方出現(xiàn)1條由下至上發(fā)展的拉裂紋，已貫穿關(guān)鍵層Ⅰ；關(guān)鍵層Ⅰ與其上巖層出現(xiàn)離層；開切眼后方的細(xì)粒砂巖層2出現(xiàn)1條由上至下發(fā)展的剪裂紋；工作面煤壁上方的泥巖層2出現(xiàn)若干條剪裂紋。

當(dāng)N=240 000時(shí)(D=110 m)(圖4(c))，采空區(qū)上方的關(guān)鍵層Ⅰ發(fā)生局部垮落，垮落區(qū)域尺寸達(dá)60 m；關(guān)鍵層Ⅰ與其上巖層的離層十分明顯；開切眼前方懸露的關(guān)鍵層Ⅰ出現(xiàn)若干條剪裂紋，工作面煤壁后方懸露的關(guān)鍵層Ⅰ已出現(xiàn)若干條由上至下發(fā)展的近似平行的傾斜剪裂紋；采空區(qū)上方的離層向上發(fā)展至距離煤層上表面約101 m處；開切眼后方與工作面煤壁前方的若干巖層出現(xiàn)多條剪裂紋；采空區(qū)上方的關(guān)鍵層Ⅱ的裂紋較為稀疏；關(guān)鍵層Ⅲ出現(xiàn)1條由下至上發(fā)展和3條由上至下發(fā)展的拉裂紋，其中，2條由上至下發(fā)展的拉裂紋貫穿關(guān)鍵層Ⅲ，可能是由于模型水平方向仍不足夠長(zhǎng)且在模型左、右兩端施加了水平約束。

當(dāng)N=450 000時(shí)(D=110 m)(圖4(d))，采空區(qū)上方的大部分離層已閉合；模型上表面發(fā)生明顯沉降；開切眼前方的采空區(qū)空間縮??；工作面煤壁后方懸露的關(guān)鍵層Ⅰ已垮落；若干巖層出現(xiàn)更多拉裂紋和剪裂紋，越向上，裂紋越向模型左、右邊界發(fā)展，這導(dǎo)致裂紋的分布在整體上呈上底長(zhǎng)下底短的梯形。此后，裂紋的分布形態(tài)基本不變，各巖層處于穩(wěn)定狀態(tài)。

圖4 模型1的σ3的時(shí)空分布(單位：Pa)

3.2 110工法條件下側(cè)向巖層運(yùn)動(dòng)模擬

當(dāng)N=40 000時(shí)(巷道的開挖、支護(hù)及頂板切縫均已完成，尚未開采)(圖5(a))，兩巷的頂、底板均出現(xiàn)垂直拉裂紋，與左巷相比，右巷頂板的裂紋更少，這是由于對(duì)右巷進(jìn)行了支護(hù)；拉裂紋尖端、切縫尖端和切縫右下方存在σ3集中；兩巷兩幫的σ3較高，左巷左、右?guī)偷摩?較高區(qū)域的寬度分別為13和14 m；右巷左、右?guī)偷摩?較高區(qū)域的寬度分別為4和14 m；在粉砂巖層與石英砂巖層中，兩個(gè)相連的V字形σ3高值區(qū)有形成的趨勢(shì)，但切縫的存在使右V字形σ3高值區(qū)右半部分基本消失；兩巷底板均存在倒V字形σ3高值區(qū)；遠(yuǎn)離巷道的煤層的σ3最低，這是由于煤層的泊松比最大，導(dǎo)致其水平變形最大，且模型左、右邊界被施加水平約束。應(yīng)當(dāng)指出，上述σ3高值區(qū)內(nèi)均為拉應(yīng)力。

當(dāng)N=120 000時(shí)(圖5(b))，兩巷之間的煤柱上方和左巷左方的粉砂巖層和石英砂巖層存在σ3集中，擠壓較為嚴(yán)重；左巷左方煤層和煤柱上方的中砂巖層和松散巖層處于拉應(yīng)力區(qū)，出現(xiàn)由上至下發(fā)展的拉裂紋；遠(yuǎn)離右巷的右方粉砂巖層和石英砂巖層處于拉應(yīng)力區(qū)，出現(xiàn)由下至上發(fā)展的拉裂紋。

當(dāng)N=800 000時(shí)(圖5(c))，在石英砂巖層、中砂巖層和松散巖層中，煤柱和左巷左方煤層的上方出現(xiàn)更多拉裂紋，且各巖層下部分出現(xiàn)成片的剪裂紋；右巷右方的巖層出現(xiàn)眾多裂紋：由下至上形成了許多近似等間距的拉裂紋，由上至下形成了成片的剪裂紋。在粉砂巖層中，煤柱上方出現(xiàn)成片的剪裂紋；右巷右方出現(xiàn)若干近似平行的拉裂紋。而且煤柱下方的底板出現(xiàn)成片的剪裂紋；右巷右方的巖層發(fā)生少許沉降；由切縫尖端已經(jīng)發(fā)展出裂紋。

當(dāng)N=1 200 000時(shí)(待采煤層已被開采，采空區(qū)已形成)(圖5(d))，由切縫尖端發(fā)展出的裂紋促進(jìn)了采空區(qū)上方的部分巖層垮落，右巷右側(cè)形成矸石幫，使巷道得以保留；煤體向右巷的擠入現(xiàn)象明顯，右巷頂板發(fā)生一定程度的下沉；在石英砂巖層、中砂巖層和松散巖層中，左巷左方煤層上方的裂紋變化不大，采空區(qū)上方出現(xiàn)更多的拉裂紋和剪裂紋，上述兩個(gè)區(qū)域之間的裂紋有所發(fā)展；采空區(qū)上方的各巖層發(fā)生明顯沉降。

圖5 模型2的σ3的時(shí)空分布(單位：Pa)

文獻(xiàn)[12]采用非連續(xù)方法UDEC呈現(xiàn)了頂板的垮落與離層。模型2的結(jié)果不僅出現(xiàn)了上述現(xiàn)象，還呈現(xiàn)了拉、剪裂紋的時(shí)空分布。

3.3 傳統(tǒng)工法條件下側(cè)向巖層運(yùn)動(dòng)模擬

當(dāng)N=40 000時(shí)(圖6(a))，采空區(qū)左方和右方的細(xì)砂巖層1分別存在一個(gè)σ3低值區(qū)，擠壓較為嚴(yán)重，寬度分別為30和32 m。上述兩個(gè)區(qū)域上方的細(xì)砂巖層1出現(xiàn)3條較長(zhǎng)的由上至下發(fā)展的拉裂紋；采空區(qū)上方的細(xì)砂巖層1出現(xiàn)7條由下至上發(fā)展的拉裂紋；細(xì)砂巖層1、泥質(zhì)砂巖層3及其之間的巖層出現(xiàn)大量?jī)A斜拉裂紋和少量剪裂紋；采空區(qū)上方部分巖層發(fā)生少許冒落。

當(dāng)N=800 000時(shí)(左、右巷道開挖已完成)(圖6(b))，采空區(qū)左方和右方的細(xì)砂巖層1與泥質(zhì)砂巖層3之間的巖層出現(xiàn)大量拉裂紋和剪裂紋：越靠近模型的左、右邊界，裂紋越垂直；采空區(qū)上方的巖層出現(xiàn)密集的拉裂紋和剪裂紋，而且，成片的剪裂紋出現(xiàn)在某些巖層(例如粗砂巖層)的上部；頂板已冒落，大部分采空區(qū)已閉合；采空區(qū)左方和右方的巷道頂板均有所沉降；左巷右?guī)秃陀蚁镒髱偷拿后w向采空區(qū)鼓出；沙土層上表面明顯彎曲，伴隨著少量拉裂。

當(dāng)N=1 200 000時(shí)(圖6(c))，左巷頂板發(fā)生明顯沉降，煤體向兩巷內(nèi)的擠入現(xiàn)象明顯；兩巷的頂、底板以及兩幫均出現(xiàn)少量剪裂紋，左巷右?guī)秃陀蚁镒髱偷拿后w向采空區(qū)的鼓出更加明顯；采空區(qū)上方的冒落帶和裂縫帶呈上底短下底長(zhǎng)的梯形，而模型中裂紋的分布在整體上呈上底長(zhǎng)下底短的梯形。此后，模型基本處于平衡狀態(tài)。

圖6 模型3的σ3的時(shí)空分布(單位：Pa)

采用連續(xù)方法FLAC3D，文獻(xiàn)[13-14]均呈現(xiàn)了采空區(qū)上方的上底短下底長(zhǎng)的梯形塑性區(qū)。模型1和模型3的結(jié)果不僅呈現(xiàn)了采空區(qū)上方的上底短下底長(zhǎng)的梯形的冒落帶與裂隙帶，還呈現(xiàn)上底長(zhǎng)下底短的梯形裂紋分布形態(tài)，這與文獻(xiàn)[15-17]中的結(jié)果(導(dǎo)水裂隙帶中的拉裂紋集中在采空區(qū)上方和開采邊界，而剪裂紋多分布在開采邊界外側(cè))定性相符。

4 結(jié)論

在拉格朗日元與離散元耦合方法的基礎(chǔ)上，發(fā)展了水平巖層運(yùn)動(dòng)計(jì)算方法，采用點(diǎn)-單元接觸模式取代單元-單元接觸模式，以起到精簡(jiǎn)算法和適于GPU并行的目的。在傳統(tǒng)工法條件下走向和側(cè)向巖層運(yùn)動(dòng)模擬中，隨著時(shí)步數(shù)目的增加，若干巖層不斷開裂，巖層之間的離層在水平和垂直方向上不斷發(fā)育，上覆巖層破斷會(huì)導(dǎo)致下方的離層和采空區(qū)閉合。最終，裂紋的分布在整體上呈上底長(zhǎng)下底短的梯形，這與導(dǎo)水裂隙帶等的現(xiàn)場(chǎng)觀測(cè)結(jié)果一致。在110工法條件下側(cè)向巖層運(yùn)動(dòng)模擬中，采空區(qū)形成后，由右巷頂板切縫尖端發(fā)展出的裂紋促進(jìn)了采空區(qū)上方部分巖層垮落，右巷右側(cè)形成了矸石幫，使其得以保留。

本方法適于模擬采場(chǎng)模型的巖層運(yùn)動(dòng)過程，例如，巖層的變形、開裂、離層、破斷和冒落等，還適于研究多種載荷(靜水壓力和沖擊波等)條件下巷道圍巖模型的破壞，例如，圍巖開裂和坍塌以及巖塊彈射等，具有需要參數(shù)少且參數(shù)容易獲取的特點(diǎn)。