基于點-單元接觸模式的水平巖層運動連續(xù)-非連續(xù)方法模擬

王學濱，薛承宇，岑子豪

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

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

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

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

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

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

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

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

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

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

根據增量形式的線彈性本構方程，由子單元的應變增量Δeij(當i=j時，表示為Δekk)求解子單元的應力增量Δσij：

(1)

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

(2)

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

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

ft=σt-σ3<0，

(3)

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

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

(4)

(5)

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

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

勢函數(shù)

(6)

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

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

2 基于點-單元接觸模式的水平巖層運動計算方法

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

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

圖2 接觸單元上一條邊的接觸點

對于點-單元接觸模式，當1個接觸點嵌入1個單元時，該接觸點上的接觸力

F=0.5KnLAhn。

(7)

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

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

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

3 巖層運動過程模擬

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

表1 3個采場模型的部分參數(shù)

圖3 3個采場模型示意圖

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

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

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

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

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

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

3.2 110工法條件下側向巖層運動模擬

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

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

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

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

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

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

3.3 傳統(tǒng)工法條件下側向巖層運動模擬

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

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

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

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

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

4 結論

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

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