基于顆粒離散元的礦柱群連鎖失穩(wěn)機理分析

周子龍，柯昌濤，王亦凡，周 靜

(中南大學資源與安全工程學院，湖南 長沙 410083)

0 引言

房柱采礦法長期以來作為地下采礦的主要方法，礦石采出后在地下產(chǎn)生了大量空區(qū)，這些空區(qū)依靠礦柱群支撐頂板以上的巖體載荷，一旦礦柱群失穩(wěn)，就會造成頂板大規(guī)模坍塌的地質災害，甚至還會形成礦難。實際上，由于受地下水、腐蝕性化學物質以及物理風化等持續(xù)的侵蝕，礦柱的功能在不斷退化，強度在持續(xù)減弱，最終會使礦柱群處于瀕臨崩塌的臨界狀態(tài)，一個微小的擾動就會誘發(fā)大規(guī)模失穩(wěn)。國內(nèi)外已經(jīng)發(fā)生了很多礦柱群大規(guī)模失穩(wěn)的災害，造成了重大的人員傷亡、嚴重的設備損毀以及生產(chǎn)被迫中斷的災難性后果[1-5]。最典型的是1960年發(fā)生在南非的礦柱群大規(guī)模失穩(wěn)的事故，近2 km2的空區(qū)坍塌，437個礦工被災難奪去生命。一系列類似的災難，引發(fā)了國內(nèi)外學術界和工程界對礦柱群大規(guī)模失穩(wěn)機理的研究和探討，例如：ZIPF R K[1]對發(fā)生過礦柱連鎖失穩(wěn)的大量案例進行分析后發(fā)現(xiàn)，這些礦山在礦柱寬高比、開挖率、邊界礦柱等方面具有相似的特征，即煤礦寬高比普遍小于3.0，金屬非金屬礦山普遍小于1.0，開挖率普遍高于60%，并且沒有邊界礦柱的存在。除了寬高比、開挖率、邊界礦柱，BE′RESTA P[2]提出頂板的剛度對礦柱群失穩(wěn)有一定的影響。POULSEN B A[6]對礦柱失穩(wěn)過程中的載荷轉移規(guī)律進行了定量的假設，并用數(shù)學方法對礦柱群連鎖失穩(wěn)的風險進行了定量的評估，引入概率方法預測礦柱群連鎖失穩(wěn)的概率；丘帆等[7]用Voronoi圖的方法估算礦柱承擔的頂板載荷面積，從而估算礦柱應力和評價系統(tǒng)穩(wěn)定性；MA Haitao等[8]研究了采空區(qū)大規(guī)模坍塌的多米諾效應；周子龍等[9]應用重整化群理論計算礦柱群的臨界失穩(wěn)概率，并且開展了雙礦柱的力學實驗和數(shù)值模擬研究，揭示了雙礦柱協(xié)同變形和強度特征[10-11]；REED GUY[12]提出了煤礦礦柱群和頂板的系統(tǒng)穩(wěn)定性標準；LI Chong等[13]則嘗試從礦柱群系統(tǒng)能量的角度出發(fā)解釋多米諾失穩(wěn)的機理。目前國內(nèi)外的研究集中于礦柱大規(guī)模失穩(wěn)的評價與案例分析，而對礦柱群失穩(wěn)的動態(tài)過程和發(fā)生機制則研究的較少。因此本文利用PFC2D顆粒離散元軟件建立了礦柱群數(shù)值模型，模擬礦柱群的漸進式失穩(wěn)過程，揭示礦柱群連鎖失穩(wěn)的誘導因素和發(fā)生機制，具有實際意義和工程應用價值。

1 數(shù)值模型

某礦采用房柱法開采，形成了面積達1.2×104m2的采空區(qū)。空區(qū)埋深為h=500 m，地表地勢平坦，圍巖密度為ρ=3 169 kg/m3，垂直地應力可估算為σy=ρgh=15.5 MPa(g=9.8 m/s2)，由于空區(qū)埋深大，水平地應力不可忽視，水平地應力可根據(jù)側壓系數(shù)估算為σx=kσy=17.05 MPa，式中k為側壓系數(shù)取1.1。礦柱寬W=6 m，高H=6 m，礦房寬L=9 m，開挖率為e=84%。根據(jù)以上數(shù)據(jù)建立PFC2D數(shù)值計算模型見圖1。模型長160 m，高56 m，頂板距離模型上邊界30 m，底板距離模型下邊界20 m，施加恒定的水平和垂直應力邊界。

圖1 PFC2D多礦柱數(shù)值模型及恒定應力邊界Fig.1 Numerical model of multiple pillars with PFC2D and its constant stress boundary condition

1.1 黏結模型及微觀參數(shù)

PFC中最基本的單元是球(三維)或單位厚度的圓盤(二維)，稱為顆粒，通過給顆粒賦予微觀物理力學參數(shù)，來模擬材料的宏觀力學行為。通過在顆粒的接觸處施加黏結模型將顆粒散體黏結為一個具有強度的宏觀物體，PFC中存在兩種黏結模型，接觸黏結和平行黏結。接觸黏結模型作用于顆粒接觸點處，等效于在顆粒的接觸點處用膠水黏結，在黏結點處只能承受力而無法承受力矩的作用。平行黏結模型是更接近材料微觀黏結本質的模型，圖2所示是平行黏結模型和顆粒接觸的示意圖，表1列出了顆粒和平行黏結微觀參數(shù)。

圖2 PFC2D顆粒接觸及平行黏結模型示意圖Fig.2 Schematic of particle and parallel bond in PFC2D

顆粒參數(shù)Rmin最小半徑0.15 mmRmax/Rmin最大最小半徑比1.66ρ顆粒密度3 169 kg/m3Ec彈性模量50 GPakn/ks法向切向剛度比2.5μ摩擦系數(shù)0.1平行黏結參數(shù)λp黏結半徑乘積系數(shù)1.0Epc黏結彈性模量50 GPakpn/kps黏結法向切向剛度比2.5σpn黏結法向強度30 MPacp黏結內(nèi)聚力30 MPaαp黏結內(nèi)摩擦角10°

平行黏結模型等效于一個充填在兩個顆粒之間的膠結體，能夠承受力和力矩的作用。在PFC2D中，平行黏結等效于一個單位厚度的長方體，其三個維度的尺寸為：

Rp=λpmin(RA,RB)

Lp=(RA+RB)/2t=1

(1)

顆粒在運動過程中，平行黏結內(nèi)部產(chǎn)生應力，當黏結體內(nèi)部應力超過強度極限時，平行黏結破裂，當微觀裂紋貫通時可以模擬宏觀破裂。給模型安裝平行黏結，在現(xiàn)場取得50 mm×100 mm長方體巖石取樣，測得單軸壓縮應力應變曲線，用于校正數(shù)值模型微觀參數(shù)。在PFC2D中建立與巖石試樣相同尺寸的數(shù)值模型，通過不斷調整微觀力學參數(shù)，使得模擬得到的單軸壓縮應力應變曲線與實驗獲得的曲線較好地吻合，對應的微觀力學參數(shù)值列于表1。

1.2 計算步驟

模型計算步驟如下：(1)建立模型，施加微觀力學參數(shù)，計算使得顆粒達到平衡狀態(tài)；(2)在模型頂部施加垂直應力σy=15.5 MPa，計算后達到平衡；(3)通過伺服機制移動垂直墻體使模型達到水平應力σx=17.05 MPa；(4)給模型施加平行黏結，開挖礦房顆粒，在1#到7#礦柱中建立應力監(jiān)測圓監(jiān)測礦柱應力，計算達到平衡狀態(tài)；(5)刪除4#礦柱，測得1#～3#和5#～7#礦柱垂直應力曲線。

2 結果分析與討論

2.1 礦柱群連鎖失穩(wěn)過程分析

圖3(a)是4#礦柱失穩(wěn)后，1#～7#礦柱垂直應力的變化情況，圖3(b)是各礦柱應力峰值時刻的裂紋發(fā)展情況，在A～F的各個時刻各礦柱的垂直應力列于表2。在A時刻，各礦柱都處于穩(wěn)定狀態(tài)，此時1#～7#礦柱垂直應力分別為37 MPa、38 MPa、39 MPa、41 MPa、39 MPa、38 MPa、37 MPa，礦柱中無裂紋出現(xiàn)。

圖3 礦柱應力變化及裂紋擴展Fig.3 Vertical stress variation and crack growth of the pillars

不同時刻礦柱垂直應力/MPa1#2#3#4#5#6#7#A37383941393837B374049-494037C384511-115038D384913-152240E433415-26849F49345-20713

4#礦柱失穩(wěn)(刪除顆粒單元)后，3#和5#礦柱最先受到影響，應力開始增加，3#礦柱和5#礦柱的應力增加速度大致相等，在B時刻均達到強度峰值，分別為49 MPa、49 MPa，3#礦柱呈東北向西南的45°剪切破壞，5#礦柱呈西北到東南的45°剪切破壞。在3#和5#礦柱應力增加過程中，1#、7#、2#、6#礦柱應力基本保持不變，3#和5#礦柱應力接近峰值的時候，2#和6#礦柱應力開始增加，6#礦柱在C時刻達到應力峰值50 MPa，礦柱呈西北到東南的45°剪切破壞。隨后在D時刻，2#礦柱達到強度峰值49 MPa，礦柱呈東北向西南的45°剪切破壞。6#礦柱達到峰值破壞時，1#和7#礦柱應力開始增加，7#礦柱在E時刻達到應力峰值49 MPa，礦柱呈西北到東南的45°剪切破壞。到F時刻，1#礦柱應力達到峰值49 MPa，礦柱呈東北向西南的45°剪切破壞，至此全部礦柱發(fā)生破壞。

根據(jù)應力曲線峰值，1#～3#，5#～7#礦柱的強度分別為49 MPa、49 MPa、49 MPa、49 MPa、50 MPa、49 MPa，根據(jù)計算公式：安全系數(shù)=礦柱強度/礦柱應力，礦柱群處于穩(wěn)定狀態(tài)即圖中A時刻的安全系數(shù)分別為：1.32、1.29、1.26、1.26、1.32、1.32。在A時刻時，礦柱的剩余強度只有10 MPa，4#礦柱失穩(wěn)后，轉移到各個礦柱上的應力使它們的載荷上升達到強度值，安全系數(shù)下降到1.0，造成鄰近礦柱的連鎖失穩(wěn)。另外，礦柱在破壞過程中的裂紋發(fā)生發(fā)展有一些共同的特點：5#、6#、7#礦柱的剪切裂紋都是從西北往東南的45°方向，而1#、2#、3#礦柱的剪切裂紋卻是從東北往西南的45°方向，這可能是由于頂板的沉降使礦柱受到彎矩的作用導致，1#、2#、3#礦柱的彎矩方向與5#、6#、7#礦柱彎矩方向相反。

2.2 不同的礦柱安全系數(shù)誘發(fā)連鎖失穩(wěn)結果分析

根據(jù)2.1節(jié)的分析，礦柱的初始安全系數(shù)分布在1.26～1.32，通過移除空區(qū)中間位置的礦柱，隨即誘發(fā)了連鎖失穩(wěn)。在本節(jié)通過提高平行黏結模型的微觀強度參數(shù)，增強礦柱的強度，保持其它一切變形參數(shù)不變，然后移除4#礦柱，觀察礦柱是否發(fā)生大規(guī)模失穩(wěn)。開展了兩組數(shù)值模擬研究，微觀參數(shù)的取值見表3。

表3 不同安全系數(shù)的平行黏結微觀參數(shù)

模型1的礦柱應力變化曲線見圖4(a)，礦柱破壞裂紋擴展見圖4(b)，在不同時刻各礦柱垂直應力值列于表4。

圖4 礦柱應力變化及裂紋擴展Fig.4 Vertical stress variation and crack growth of pillars

不同時刻礦柱垂直應力/MPa1#2#3#4#5#6#7#A37383941393837B374354—544337C38495—125441D40525—204543E50285—301052F55245—25732

圖4(b)礦柱的破壞順序與圖3(b)一樣，圖4(a)的應力變化曲線與圖3(a)也十分相似。在A時刻的穩(wěn)定狀態(tài)時，1#～7#礦柱的垂直應力分別為37 MPa、38 MPa、39 MPa、41 MPa、39 MPa、38 MPa、37 MPa。刪除4#礦柱后，在B時刻，3#和5#礦柱應力達到峰值54 MPa，3#和5#礦柱破壞后，2#和6#礦柱應力開始大幅增加，分別在D時刻和C時刻達到峰值52 MPa、54 MPa，2#和6#礦柱破壞后，1#和7#礦柱應力開始增加，分別在F時刻，E時刻達到峰值強度55 MPa、52 MPa。根據(jù)應力曲線峰值，1#～3#，5#～7#礦柱的強度分別55 MPa、52 MPa、54 MPa、54 MPa、54 MPa、52 MPa，計算得到礦柱的初始安全系數(shù)分別為：1.49、1.37、1.38、1.38、1.42、1.41，安全系數(shù)分布在1.37～1.49。模型1的結果與圖3相比，微觀參數(shù)增大后，礦柱強度提高，礦柱安全系數(shù)也相應提高，因此礦柱能夠承載更多的轉移載荷，模型1中礦柱承載的轉移載荷增大為18 MPa，相比圖3的10 MPa有所提高，但是礦柱仍然發(fā)生了大規(guī)模的連鎖失穩(wěn)，可知此模型的載荷轉移量大于18 MPa，要想避免連鎖失穩(wěn)，必須進一步提高礦柱安全系數(shù)值。

模型2相對模型1具有更高的礦柱強度，其它所有參數(shù)保持一致。得到的應力變化曲線見圖5。

圖5 礦柱應力變化曲線Fig.5 Vertical stress variation of pillars

由圖5可知，4#礦柱失穩(wěn)后，3#和5#礦柱的應力分別增加到54 MPa和55 MPa后達到穩(wěn)定的狀態(tài)，2#和6#礦柱最終應力增加到44 MPa后達到穩(wěn)定狀態(tài)，1#和7#礦柱應力基本保持不變，礦柱群的連鎖失穩(wěn)沒有發(fā)生。由此推斷，在圖1所示的礦柱和頂板結構形狀和地應力狀態(tài)以及表1所示的巖石變形微觀力學參數(shù)條件下，礦柱群發(fā)生連鎖失穩(wěn)的臨界安全系數(shù)就位于模型1和模型2之間。

2.3 不同的礦柱破壞位置誘發(fā)連鎖失穩(wěn)結果分析

2.1節(jié)和2.2節(jié)模擬了采空區(qū)中央位置的礦柱失穩(wěn)誘發(fā)連鎖失穩(wěn)的情況，以及誘發(fā)連鎖失穩(wěn)的礦柱安全系數(shù)臨界條件。通常，由于采空區(qū)頂板的沉降，中央位置的頂板沉降最大，使整個空區(qū)頂板呈現(xiàn)出一個類似漏斗的形狀，由此造成空區(qū)中央位置的礦柱應力最大，越靠近空區(qū)邊緣，礦柱應力越小，空區(qū)中央位置的礦柱失穩(wěn)帶來的擾動效應一般來說也最嚴重。實際工程中的采空區(qū)礦柱形狀往往不規(guī)則，礦柱布局也不均勻，礦柱應力大小不一，因此不一定是空區(qū)中央的礦柱最先失穩(wěn)，本節(jié)將探討空區(qū)邊緣的礦柱失穩(wěn)是否會誘發(fā)礦柱連鎖失穩(wěn)的問題。基于圖1所示的數(shù)值模型和表1的顆粒微觀力學參數(shù)，模型達到平衡狀態(tài)時，移除邊緣的1#礦柱，得到的應力變化曲線如圖6所示。

圖6 邊界礦柱失穩(wěn)后應力變化曲線Fig.6 Vertical stress variation and crack growth of pillar 2#～6# after pillar 1# fails

A時刻礦柱處于穩(wěn)定狀態(tài)，1#～7#礦柱的垂直應力分別為37 MPa、38 MPa、39 MPa、41 MPa、39 MPa、38 MPa、37 MPa，根據(jù)2.1的分析，礦柱初始安全系數(shù)分別為：1.29、1.26、1.26、1.26、1.32、1.32。1#礦柱移除后，2#礦柱的應力上升后又小幅下降最終穩(wěn)定在52 MPa，最終安全系數(shù)下降到0.94。2#礦柱安全系數(shù)小于1.0，但是圖6(b)B時刻的裂紋圖表明2#礦柱中并沒有形成貫通的裂紋，2#礦柱并沒有發(fā)生破壞，可能是因為兩個方面的原因：一是礦柱內(nèi)應力分布不均勻，局部應力很高，造成了礦柱的局部破壞，如圖6(b)，2#礦柱左下角和右上角有少量裂紋擴展，但并沒有貫通；二是邊界礦柱對頂板的端部支撐作用，此時頂板類似懸臂梁，沉降得到有效控制，沒有對2#礦柱形成持續(xù)的加載，圖6(a)也證明2#礦柱的應力最終維持平穩(wěn)狀態(tài)，因此沒有造成2#礦柱的破壞。與此不同的是，中央位置的礦柱失穩(wěn)，由于離邊界礦柱遠，邊界礦柱對頂板的沉降發(fā)揮不了有效的約束作用，頂板的沉降對礦柱形成了持續(xù)的加載，直到礦柱破壞。2#礦柱未失穩(wěn)，發(fā)揮了類似邊界礦柱的作用，因此3#～7#礦柱也保持穩(wěn)定，3#礦柱應力最終增加到44 MPa后穩(wěn)定，4#礦柱應力增加到43 MPa后穩(wěn)定，而5#、6#、7#礦柱的應力基本不變，沒有受到1#礦柱失穩(wěn)的影響。

2.4 不同剛度對礦柱群連鎖失穩(wěn)的影響

在PFC2D中，顆粒法向剛度與彈性模量成正比，因此要研究剛度對礦柱群連鎖失穩(wěn)的影響，只需要改變顆粒彈性模量Ec的值就能實現(xiàn)不同的剛度。在圖1數(shù)值模型和表1微觀力學參數(shù)的基礎上，分別將彈性模量縮小10倍和增大10倍，不改變法向剛度切向剛度的比值，其它參數(shù)也保持不變?nèi)绫?所示，分別建立兩組模型，按1.2節(jié)的步驟開展計算。

表5 不同頂板顆粒剛度參數(shù)

在模型1的軟弱頂板條件下，得到的礦柱應力變化曲線見圖7，各峰值時刻各礦柱的垂直應力值見表6。

圖7 礦柱應力變化曲線Fig.7 Vertical stress variation of pillar 1#～7# after 4# fails

所得的圖7礦柱破壞的順序與圖3基本一致，呈現(xiàn)出漸進的破壞順序，4#礦柱失穩(wěn)瞬間，5#和3#礦柱首先在B時刻破壞，6#和2#礦柱分別在C時刻和D時刻破壞，最后破壞的是位于空區(qū)兩邊的1#和7#礦柱，破壞順序呈現(xiàn)出典型的“多米諾骨牌效應”。由表6可知，1#～3#、5#～7#礦柱的強度分別45 MPa、50 MPa、49 MPa、49 MPa、49 MPa、50 MPa。

表6 軟弱頂板條件下礦柱垂直應力變化值

模型2的堅硬頂板條件下所得的計算結果如圖8所示。各時刻礦柱垂直應力值見表7。

圖8 礦柱應力變化曲線Fig.8 Vertical stress variation of pillar 1#～7# after 4# fails

不同時刻礦柱垂直應力/MPa1#2#3#4#5#6#7#A37383941393837B464750—504746C404350—104740D384345—125342E38509—101546F404710—131347

由圖8可知，在堅硬頂板條件下，所得和結果與軟弱頂板條件完全不同，圖8的應力曲線變化方式與圖7有本質的區(qū)別。在刪除4#礦柱的瞬間，即A時刻，1#～3#、5#～7#礦柱的應力都明顯增加，在B時刻同時達到峰值，3#和5#礦柱應力增加了11 MPa，2#和6#礦柱應力增加了9 MPa，1#和7#礦柱應力也增加了11 MPa，說明在頂板剛度非常大的情況下，4#礦柱的失穩(wěn)瞬間使得頂板的載荷均勻地轉移給剩余的礦柱。B時刻之后，1#～3#、5#～7#礦柱的應力幾乎不再增長，5#礦柱率先失穩(wěn)，應力迅速下降，1#和3#礦柱的應力經(jīng)過一段平穩(wěn)的時期也開始迅速下降，2#、6#、7#礦柱應力有微小的增加后也迅速破壞，因此總結圖8中各礦柱的破壞特征，可以發(fā)現(xiàn)礦柱群失穩(wěn)的順序不再具有圖7漸近式失穩(wěn)的特征，因此在堅硬頂板條件下礦柱群失穩(wěn)不再顯著地表現(xiàn)出“多米諾骨牌效應”。由此可以證明，頂板剛度在礦柱群大規(guī)模失穩(wěn)中確實發(fā)揮了重要的作用，這種作用主要體現(xiàn)在影響載荷轉移上，當頂板較軟弱時，礦柱的載荷轉移從失穩(wěn)誘發(fā)礦柱開始，離失穩(wěn)誘發(fā)礦柱越遠，受到載荷轉移的影響就越遲；但是當頂板是剛度很大的硬巖時，在失穩(wěn)瞬間，載荷幾乎同時均勻地轉移到剩余穩(wěn)定的礦柱中。值得一提的是，PFC模擬巖石材料時的微觀顆粒彈性模量不可能達到500 GPa之大，因此在礦柱群失穩(wěn)過程中，載荷均勻轉移到所有礦柱上的情況較少，圖7的礦柱應力轉移模式是最常見的。

3 結論

結合工程實例，基于PFC2D顆粒流離散元建立了礦柱群的數(shù)值模型，模擬了礦柱群大規(guī)模的連鎖失穩(wěn)的發(fā)生機理，以及載荷轉移過程，研究了不同的礦柱安全系數(shù)、不同的誘發(fā)位置、不同頂板剛度對礦柱群大規(guī)模失穩(wěn)的影響。得到了以下結論：

(1)礦柱群連鎖失穩(wěn)是一個載荷轉移的過程，連鎖失穩(wěn)從最先破壞的失穩(wěn)誘發(fā)礦柱位置開始，逐漸向其它礦柱轉移載荷，由此引發(fā)連鎖失穩(wěn)，表現(xiàn)出典型的“多米諾骨牌效應”。

(2)沒有足夠的安全系數(shù)儲備是礦柱群連鎖失穩(wěn)的主要原因，安全系數(shù)小使得礦柱沒有足夠的剩余強度承擔轉移過來的載荷，而充足的安全系數(shù)可以使礦柱發(fā)揮邊界礦柱的作用，有效阻擋失穩(wěn)擾動向外擴散。

(3)礦柱群連鎖失穩(wěn)是通過一個或幾個礦柱誘發(fā)產(chǎn)生的，誘發(fā)的位置對連鎖失穩(wěn)產(chǎn)生了顯著的影響，空區(qū)中央位置是礦柱應力最大的位置，中央位置失穩(wěn)誘發(fā)容易產(chǎn)生連鎖失穩(wěn)，而離邊界礦柱近的位置失穩(wěn)不容易誘發(fā)連鎖失穩(wěn)，因為邊界礦柱對頂板具有強大的支撐約束作用。

(4)頂板剛度影響載荷轉移的順序，軟弱的頂板使得載荷從失穩(wěn)礦柱位置開始逐漸向外擴散轉移，堅硬的剛性頂板使得載荷在誘發(fā)礦柱失穩(wěn)瞬間，迅速均勻地轉移到所有未失穩(wěn)的礦柱上。