譚寶會，梁 博，賈凱躍，李明潤

(西南科技大學 環(huán)境與資源學院，四川 綿陽 621010)

0 引言

下向分層膠結充填法因具有安全、高效、低貧損等優(yōu)點而被廣泛應用于貴重有色金屬礦山和稀有金屬礦山地下開采中，但當某些特定因素發(fā)生改變后(如礦石價格下跌、品位降低等)，礦山可能不再適合采用生產(chǎn)成本較高的膠結充填法進行采礦，轉而考慮采用生產(chǎn)成本較低的無底柱分段崩落法。此時崩落法采場頂板以上通常至少包含了兩種屬性明顯不同的材料，其中一種是人工構筑的膠結充填體，另一種是充填體上部的天然圍巖，這兩種材料相互疊加構成了一種特殊的人造-天然復合巖層。崩落法采場在首采分段采礦過程中需要誘導頂板復合巖層冒落，形成正常采礦所需的散體覆蓋層。因此研究復合巖層的冒落機理及模式，對于制訂崩落法采礦方案及安全管控措施有著重要意義。但崩落法采空區(qū)頂板的冒落通常發(fā)生在無法直接觀測的崩落區(qū)內，這就增加了冒落機理、冒落模式及冒落發(fā)展過程的研究難度。

近年來，數(shù)值模擬方法的快速發(fā)展為頂板圍巖的冒落機理研究提供了一條便捷途徑。數(shù)值模擬方法可以很好地分析地下采礦過程中巖體的位移、應力及應變的分布模式[1-3]。因此有學者將連續(xù)介質數(shù)值模擬方法[如有限元方法(FEM)、有限差分方法(FDM)]應用于地下采礦引發(fā)的頂板冒落機理研究，但是地下采礦誘發(fā)頂板巖體冒落為松散巖塊，冒落前后巖體力學特性將發(fā)生顯著改變，因此利用FEM、FDM等方法模擬頂板冒落時，將巖體定義為一種連續(xù)的彈塑性材料并利用彈塑性理論將頂板的冒落歸因于材料的屈服是有一定缺陷的。而離散元方法(DEM)、非連續(xù)變形分析(DDA)以及數(shù)值流形方法(NMM)在模擬巖體不連續(xù)的大變形、大位移方面具有顯著優(yōu)勢，被廣泛應用于模擬地下開挖引發(fā)的頂板冒落行為[4-7]。DDA方法是一種對不連續(xù)塊體系統(tǒng)的靜態(tài)和動態(tài)力學行為進行計算的數(shù)值模擬方法，其與DEM相比，共同之處在于模擬對象都是離散的塊體單元，但由于計算收斂性等因素，DDA方法在模擬不連續(xù)塊體的大位移、大變形等方面(如巖體冒落、邊坡滑坡等)更具優(yōu)勢[8-9]。

綜合比較當前可利用的數(shù)值模擬工具，采用DDA或者DEM等方法對巖體的冒落機理進行研究更為合適。本文以某礦山下向分層膠結充填法轉無底柱分段崩落法為工程背景，采用DDA方法對崩落法采場頂板膠結充填體及上部圍巖構成的復合巖層冒落機理、模式及過程展開研究，以期為崩落法采場覆蓋層形成方案的設計和實施提供指導。

1 DDA方法的改進

1.1 DDA子塊體開裂算法的優(yōu)化

在以往的預離散子塊體DDA開裂模擬方法中，通過采用虛擬節(jié)理對子塊體進行黏結以模擬連續(xù)體的變形，進而根據(jù)虛擬節(jié)理面上的塊間接觸力判斷沿預設節(jié)理面的拉伸或剪切破壞情況；該算法能夠較好地模擬巖體的開裂路徑和破壞形態(tài)，但受預設虛擬節(jié)理面方向對塊間接觸力大小的影響，得到的巖體開裂強度與實際可能存在較大差異。為解決這一問題，倪克松等[10]對DDA方法中的子塊體開裂模擬算法進行了優(yōu)化，在該算法中不再根據(jù)塊間接觸力判斷開裂，而改進為根據(jù)鄰近子塊體的應力狀態(tài)判斷開裂，并仍假定裂紋沿虛擬節(jié)理面產(chǎn)生。研究表明，改進后的開裂算法具有較高的準確性，并減小了虛擬節(jié)理分布對開裂破壞路徑模擬結果的影響。

1.2 DDA塊體切角功能的實現(xiàn)

一般來說，在DDA模型中劃分的子塊體是最小的塊體單元，無法二次開裂，因此當這些塊體帶有銳角時在冒落或移動過程中就極易與周邊塊體相互擠壓而卡死，從而影響塊體運動行為的發(fā)展。為解決這一問題，在DDA前處理程序中添加了塊體的切角功能，假設模型中某個塊體為銳角三角形，其中一條邊長為L，通過設定切角比例k，便可在該條邊的首尾兩端各切去kL/2，從而消除了塊體中的尖角(見圖1)。

圖1 改進的DDA塊體切角功能

實驗結果表明，切角后的塊體系統(tǒng)發(fā)生不正?？ㄋ赖母怕蚀蟠蠼档汀４送?，切角后的模型在一定程度上也更符合真實的巖體結構，但要確定合理的切角比例則需要進行大規(guī)模的巖體結構信息調查，這在實際操作中難度較大，但從減少塊體系統(tǒng)卡死的概率來說，只要進行適當切角，其作用就會非常明顯。

2 數(shù)值模型構建

2.1 模型的建立及數(shù)值控制參數(shù)設定

選取垂直礦體走向的典型剖面建立DDA模型，模型從下至上按照實際情況劃分為4種不同屬性的材料(不包含模型邊框)，其中最底層為10 m厚的礦體，該層代表實際中崩落法采場的首采分段?？紤]到實際中礦體的強度相對較大，對頂板冒落所產(chǎn)生的影響較小，因此為了盡可能地縮減模型中塊體的數(shù)量，將礦體材料的網(wǎng)格簡化為9個矩形剛性塊體，模擬時對位于中部的7個塊體從左至右依次開挖，以代表實際中沿進路方向的退采，每個開挖塊體的寬度為10 m，大約等同于實際中的4～5個回采步距。

礦體之上為30 m厚的膠結充填體，共包含6個充填分層。膠結充填體中包含了兩種不同強度的不連續(xù)面：一種為上下分層交界處的水平不連續(xù)面，其代表的是實際中膠結充填體上下分層的接觸面，該不連續(xù)面的強度要小于膠結充填體的強度；另一種不連續(xù)面為膠結充填體內部的虛擬節(jié)理，其強度可取膠結充填體試樣的實測強度。

本次建模采用了蜂巢狀網(wǎng)格作為膠結充填體的虛擬節(jié)理，以便更真實地模擬膠結充填體的破壞形態(tài)，同時考慮到當前版本的DDA程序尚無法實現(xiàn)子塊體的二次開裂，這在一定程度上限制了膠結充填體剪切滑移破壞的發(fā)展，為了彌補這一缺點，在膠結充填體中增加了一組垂直方向上的貫通虛擬節(jié)理。

膠結充填體之上為50 m厚的片麻巖，片麻巖中包含了兩組優(yōu)勢結構面，其傾角及傾向分別為88°∠56°和3°∠52°；由于這兩組節(jié)理均具有較好的連續(xù)性，因此在模型中將這兩組節(jié)理均設定為貫通節(jié)理，同時為了縮減模型中的塊體數(shù)量，在保持這兩組節(jié)理走向、傾角及接觸面力學參數(shù)貼近真實的基礎上，將節(jié)理間距擴大10倍，這種做法雖然會在某些方面降低模擬結果的真實程度(如巖體開裂后的碎脹系數(shù)等)，但并不會從本質上影響巖層的失穩(wěn)模式。

片麻巖之上為10 m厚的第四系表土，由于第四系表土強度極低，十分松散，極易發(fā)生破壞，因此為進一步縮減模型中的塊體數(shù)量，僅在模型中設定了10 m厚的第四系表土為代表，同時采用2組互相正交的虛擬節(jié)理來對其強度進行賦值并表征破壞。最終模型中塊體數(shù)量維持在4.0×103個左右，模型幾何尺寸及測點分布情況如圖2所示。

模型中各材料的關鍵物理力學參數(shù)是根據(jù)礦山提供的地質報告以及一系列相關的物理實驗獲得的。為了使模型中各材料的破壞更貼近真實情況，設定當模型中某個不連續(xù)面發(fā)生破壞后，該不連續(xù)面上的黏聚力、抗拉強度、摩擦角等強度參數(shù)均會發(fā)生衰減[11]。模型中各材料的物理力學參數(shù)見表1。

圖2 DDA模型幾何尺寸及測點布置

表1 模型中各巖層的關鍵物理力學參數(shù)

DDA方法所用到的數(shù)值控制參數(shù)是由用戶自定義的，合理的數(shù)值控制參數(shù)通常要經(jīng)過大量的前期試算才能得到。采用合理的數(shù)值控制參數(shù)在進行求解計算時，每一時步的計算會在最小次數(shù)的開-合迭代之后達到收斂[12]。本次模擬使用的數(shù)值控制參數(shù)見表2。

表2 模型求解所用的數(shù)值控制參數(shù)

2.2 模擬方案設定及初始應力平衡求解

為了使模擬更貼近現(xiàn)實，模型在開挖之前需要先在靜態(tài)條件下求解初始應力平衡，即在模型中先實現(xiàn)地應力的還原再進行開挖。經(jīng)現(xiàn)場地應力測試，該礦區(qū)地應力以水平應力為主，垂向壓力基本為上覆巖體自重，各測點的水平應力與垂向應力之比(側壓比)λ大約維持在1～2，且最大水平應力方向基本垂直于礦體走向。為研究水平應力對復合巖層冒落發(fā)展的影響，本次數(shù)值模擬構建了兩組模擬方案，即方案1和方案2，方案1側壓比λ=1，方案2側壓比λ=2。

圖3以模型中1#測點為例展示了兩組模擬方案的初始應力平衡求解過程。1#測點位置對應的實際埋深在130 m左右，由理論計算可知，該測點處的垂向壓力σv=3 MPa。因方案1側壓比為1，故該測點的水平應力σh=3 MPa；而方案2側壓比為2，故該測點的水平應力σh=6 MPa。

圖3 模型求解初始應力平衡

從圖3可以看出，兩組模型在物理時間t=20 s時(共計算20 000時步)達到了初始應力平衡。在方案1中，1#測點處的垂向應力達到3 MPa，水平應力也達到3 MPa。在方案2中，1#測點的垂向應力達到3 MPa，水平應力基本達到6 MPa?？梢妰山M方案的應力條件基本達到預期狀態(tài)，此時便可以從左至右依次移除底部所設定的7個塊體，執(zhí)行開挖后程序自動轉為動態(tài)條件進行計算。模型中第1步開挖設定在20 000時步，第二步開挖設定在60 000時步(t=24 s)，之后每隔10萬時步進行一次開挖，兩次開挖之間保持較長時間間隔的目的是使每次開挖引發(fā)的頂板冒落都能夠得到充分發(fā)展。

3 模擬方案及結果分析

圖4給出了模擬方案1中復合巖層隨著開挖逐漸發(fā)生失穩(wěn)冒落的動態(tài)演化過程。由圖4(a)可知，在第1步開挖后，形成了跨度為10 m的采空區(qū)，此時空區(qū)頂板圍巖幾乎未發(fā)生任何開裂和破壞。由圖4(b)可知，在第2步開挖后，空區(qū)跨度達到了20 m，此時頂板膠結充填體開始發(fā)生破壞，膠結充填體沿空區(qū)兩側逐層向上破壞，很快形成貫通破壞裂隙，延伸至第五分層時才中止，此時上部片麻巖也發(fā)生了部分開裂，最終第一層膠結充填體巖梁發(fā)生了彎曲離層，但未冒落，而是在巖梁兩端摩擦阻力及水平應力作用下形成了穩(wěn)定的裂隙梁，在此過程中有個別散塊冒落至空區(qū)底板。由圖4(c)可知，在第3步開挖后，空區(qū)跨度達到了30 m，此時空區(qū)上方的膠結充填體、片麻巖和第四系表土瞬間完全垮塌冒落，以柱狀形式整體滑落至空區(qū)底板，此階段可視為空區(qū)頂板突發(fā)大規(guī)模冒落。由圖4(d)可知，第4步開挖后，新開挖空區(qū)上部的膠結充填體及覆巖呈現(xiàn)為臨空的懸臂梁形式，懸臂梁在自重作用下于支點處發(fā)生折斷并冒落至空區(qū)。為節(jié)省計算時間，在第4步開挖引發(fā)的冒落完成后終止了計算，可以推斷的是在后續(xù)退采中新開挖空區(qū)上方的復合巖層將周期性地發(fā)生懸臂梁折斷冒落。

(a)開挖跨度10 m

圖5為模擬方案2中復合巖層隨著開挖逐漸發(fā)生失穩(wěn)冒落的動態(tài)演化過程。

(a)開挖跨度10 m

由圖5(a)可知，在第1步開挖后，空區(qū)跨度為10 m，此時空區(qū)頂板各巖層幾乎未發(fā)生任何開裂和破壞。由圖5(b)可知，在第2步開挖后，空區(qū)跨度達到了20 m，此時僅有空區(qū)頂板膠結充填體的第一分層產(chǎn)生了少許離層裂紋。由圖5(c)可知，在第3步開挖后，空區(qū)跨度達到了30 m，此時第一分層膠結充填體冒落至空區(qū)底板，同時上部分層中也開始產(chǎn)生少量裂紋。由圖5(d)可知，在第4步開挖后，空區(qū)跨度達到了40 m，此時第二分層膠結充填體發(fā)生了彎曲離層，第3-第5分層兩端發(fā)生貫通式開裂，但并未滑落，而是形成了穩(wěn)定的裂隙梁，同時空區(qū)上方對應的片麻巖和第四系表土中產(chǎn)生大量裂紋。由圖5(e)可知，在第5步開挖后，空區(qū)跨度達到了50 m，此時膠結充填體和頂板圍巖突發(fā)大規(guī)模整體滑落。由圖5(f)可知，第6步開挖后，新開挖空區(qū)上部復合巖層以臨空懸臂梁折斷的形式發(fā)生破壞。在模型第6步開挖引發(fā)的冒落活動穩(wěn)定后，終止了計算，可以推斷的是在后續(xù)退采中新開挖空區(qū)頂板圍巖將周期性地發(fā)生懸臂梁折斷垮塌。

對比方案1和方案2的冒落過程可知，在該礦山膠結充填體和頂板圍巖的強度和結構特征條件下，崩落法采場空區(qū)頂板復合巖層頂板很容易發(fā)生柱狀整體陷落，從而引發(fā)直通地表的大規(guī)模冒落，而且水平應力顯著提升了各巖層的抗破壞能力，增加了頂板的臨界冒落跨度，但采空區(qū)的跨度越大，突發(fā)冒落的規(guī)模也就越大，更容易引起冒落沖擊災害。

為深入分析空區(qū)頂板的冒落機理，圖6以方案1為例給出了模型計算過程中的最大剪應力分布云圖，以及空區(qū)上方巖層整體冒落后模型中的破壞裂紋分布模式。

(a)第60 600時步最大剪應力云圖

由圖6可知，模型開挖后空區(qū)兩側發(fā)生了明顯的剪應力集中，在初期時剪應力主要集中在膠結充填體第一分層巖梁兩端[見圖6(a)]，率先致使第一分層巖梁兩端發(fā)生剪切破壞；隨后，集中剪應力轉移至第二分層，使第二分層巖梁兩端發(fā)生剪切破壞，以此類推，巖梁兩端的剪切破壞逐層向上發(fā)展[見圖6(b)、圖6(c)]。由于膠結充填體強度較低，剪切裂紋很快就貫通了整個膠結充填體，同時膠結充填體在自重作用下發(fā)生彎曲下沉，致使巖梁中部產(chǎn)生拉應力集中，從而發(fā)生拉伸破壞[見圖6(d)]，嚴重破壞了膠結充填體的完整性。當巖梁兩端所受的水平擠壓力和摩擦力不足以阻止膠結充填體滑落時，膠結充填體便會發(fā)生冒落。失去膠結充填體的支撐后，上部片麻巖在重力作用下沿陡立的貫通節(jié)理發(fā)生整體剪切下滑，同時強度較低的片麻巖也隨之陷落?？傊?，由于膠結充填體和頂板圍巖強度均較低，空區(qū)頂板圍巖往往會在較短的時間內發(fā)生整體垮塌，從而造成大規(guī)模冒落。

為了展現(xiàn)開挖引發(fā)的應力重分布規(guī)律，圖7以模擬方案1為例給出了模型中1#-4#測點的垂向應力演化曲線。由圖7可知，在模型連續(xù)開挖到40 m跨度的過程中，當測點位于空區(qū)一側時，測點處將產(chǎn)生垂向應力集中現(xiàn)象，且距離空區(qū)越近，垂向應力集中現(xiàn)象越明顯；而隨著開挖的進行，一旦測點的位置由空區(qū)一側轉為空區(qū)頂板之上后，其垂向應力被立即釋放，直至測點所在塊體冒落至空區(qū)底板，上部冒落巖體逐漸堆積在該塊體上，該測點的垂向壓力才逐漸出現(xiàn)回升，此時該測點上的垂向應力即為上部散體產(chǎn)生的垂壓。在整個開挖過程中，3#測點和4#測點始終位于空區(qū)一側。由圖7(c)可以看出，3#測點在前兩步開挖后垂向應力基本未發(fā)生改變，當?shù)谌介_挖后，垂向應力由3.0 MPa增加到3.5 MPa左右，增加幅度為17%左右，可認為此時該測點受到了采動壓力的影響；由圖7(d)可以看出，4#測點的垂向應力自始至終基本不變，表明在開挖過程中該測點未受到采動壓力的影響。

(a)1#測點

數(shù)值模擬結果表明，采空區(qū)頂板發(fā)生冒落的根本原因是開挖造成模型中應力重新分布；具體來說是開挖后空區(qū)上方垂向應力向空區(qū)兩側轉移，致使空區(qū)兩側產(chǎn)生剪應力集中，膠結充填體在空區(qū)兩端的集中剪應力作用下逐層向上發(fā)生剪切破壞，同時橫梁彎曲下沉時在中部產(chǎn)生拉應力集中，引發(fā)拉伸破壞，頂板圍巖在失去膠結充填體的支撐后，便在自重作用下沿軟弱的陡立結構面發(fā)生柱狀剪切滑移破壞。

由于該礦區(qū)礦體埋藏較淺，水平應力較大，同時膠結充填體強度較低而頂板圍巖富含弱陡立節(jié)理，在這些因素的共同作用下，崩落法采場回采過程中頂板復合巖層易突發(fā)直通地表的大規(guī)模冒落。此時，氣浪沖擊危害和沖擊振動危害將是最可能致災的兩大危險源。工程實踐證明[13-15]，在崩落法首采分段退采過程中均勻擴展采空區(qū)，避免在空區(qū)內部留支撐礦柱，可促使空區(qū)頂板及時發(fā)生冒落；同時在空區(qū)底板預留足夠厚度的礦石散體墊層，既可以有效緩沖頂板冒落的沖擊振動作用，還可封堵出礦口以預防氣浪沖擊。

4 結論

a.數(shù)值模擬結果表明，隨著采空區(qū)的不斷擴大，空區(qū)上方垂向應力逐漸向空區(qū)兩側轉移，致使空區(qū)兩側產(chǎn)生剪應力集中，膠結充填體在空區(qū)兩端的集中剪應力作用下逐層向上發(fā)生剪切破壞；由于膠結充填體的強度較低，剪切破壞很快貫通整個膠結充填體，同時在自重作用下巖梁發(fā)生彎曲下沉，致使巖梁中部產(chǎn)生拉應力集中，從而發(fā)生拉伸破壞，嚴重破壞了膠結充填體的完整性；當巖梁兩端所受的水平擠壓力和摩擦力不足以阻止膠結充填體滑落時，膠結充填體便發(fā)生冒落。

b.由于礦山片麻巖中富含貫通的光滑陡立節(jié)理，當膠結充填體發(fā)生失穩(wěn)下滑后，上部片麻巖便失去支撐，極易在重力作用下沿陡立節(jié)理發(fā)生整體剪切滑移破壞，同時由于第四系表土的抗剪強度極低，當其臨空后也極易發(fā)生整體陷落。

c.研究結果表明，由于該礦區(qū)礦體埋藏較淺，水平應力較大，同時膠結充填體強度較低而頂板圍巖又富含連續(xù)性較好的光滑陡立節(jié)理，因此在崩落法回采過程中采空區(qū)頂板復合巖層易突發(fā)直通地表的柱狀大規(guī)模冒落，此時氣浪沖擊和沖擊振動將是復合巖層冒落過程中最可能致災的兩大危險源，在實際生產(chǎn)過程中，建議采取均勻擴展采空區(qū)和預留礦石散體緩沖墊層等措施進行防控。