李昊城，宋選民，朱德福,2

(1.太原理工大學 原位改性采礦教育部重點實驗室，山西 太原 030024；2.中國礦業(yè)大學 煤炭資源與安全開采國家重點實驗室，江蘇 徐州 221116)

0 引 言

西北礦區(qū)煤炭資源儲量豐富，可采煤層數(shù)量多，且層間距較小。由于早期開采技術相對落后，上組煤層多采用房式開采，遺留大面積房式采空區(qū)，據(jù)統(tǒng)計，僅鄂爾多斯境內的房式采空區(qū)面積就達307.61 km2[1]。房式采空區(qū)的獨特性在于其上覆巖層自重均由遺留煤柱群來承擔，煤柱的穩(wěn)定性決定了整個采空區(qū)的穩(wěn)定狀態(tài)。下組煤層開采過程中，煤柱上方應力重新分布，造成煤柱的最大主應力增大，當上覆載荷大于煤柱極限強度時，煤柱塑性破壞區(qū)將由四周逐漸向中心擴展，直至發(fā)生整體塑性破壞從而失去支撐頂板的能力[2]，此過程易誘發(fā)煤層頂板聯(lián)動垮落形成強礦壓致災[3]。因此，評價房式煤柱穩(wěn)定性對下組煤層群安全開采具有重要的意義。房式采煤法是柱式采煤法的一種，在煤層每隔一定距離開采煤房，煤房之間留有一定寬度的煤柱支撐頂板，煤房及煤柱的尺寸設計取決于采高、埋深、煤層頂?shù)装辶W特性等[4]。房、柱分布規(guī)律分為平行分布(棋盤形)、“品”字形分布兩類，煤柱為采空區(qū)內的唯一承載結構[5-6]。影響房式采空區(qū)煤柱穩(wěn)定性的因素主要有煤柱分布方式、上覆載荷及力學、幾何特性參數(shù)等[7]。Voronoi圖法、模糊理論、重整化群理論廣泛應用于不規(guī)則房式煤柱群、條帶煤柱、非煤礦柱的穩(wěn)定性研究中[8-10]，認為影響煤柱群失穩(wěn)的主要因素為柱群的分布形式、各柱體間尺寸及強度的均勻性。為具體研究煤柱群中某一煤柱失穩(wěn)破壞特征，基于煤柱群穩(wěn)定性的研究基礎，結合煤柱上覆載荷分布及傳遞特性發(fā)展出了概率逼近法、邏輯回歸法及蒙特卡洛法等數(shù)理統(tǒng)計方法[11-13]，這類方法考慮了煤柱尺寸及力學參數(shù)的不確定性及其他因素，在房式煤柱的穩(wěn)定性評估中應用廣泛。隨后尖點突變失穩(wěn)模型將拓撲學、奇點理論等非線性理論引入到煤柱失穩(wěn)的機理中，得出突變失穩(wěn)對應的控制變量變化路徑及臨界條件[14]。由于房式采空區(qū)中煤柱的穩(wěn)定性是一個動態(tài)演化的過程且受眾多因素的影響，需充分考慮煤巖單元塑性弱化、能量傳遞及損傷發(fā)展等非線性變化過程。元胞自動推理機可利用簡單單元之間的相互作用規(guī)則，得到整個離散動力系統(tǒng)自組織演化的非線性數(shù)學模型[15]，從而用于巖石類材料受載后動態(tài)演化過程的研究中[16-18]。上述研究成果從煤柱系統(tǒng)穩(wěn)定性的角度分析并研究了煤柱的長期承載及失穩(wěn)特征，而針對煤柱分形單元的損傷演化涉及較少，基于上述研究，采用元胞自動機的方法建立二維8Moore元胞房式煤柱穩(wěn)定性評價模型，結合能量傳遞及元胞儲能極限準則，動態(tài)演化煤柱元胞破壞失穩(wěn)規(guī)律，通過有限差分模擬程序揭示房式煤柱單元失穩(wěn)時序，確定煤柱內元胞動態(tài)失穩(wěn)模式煤柱穩(wěn)定特征，并結合現(xiàn)場試驗評價模型的可靠性。

1 房式煤柱元胞自動機模型

1.1 元胞上覆均布載荷計算

選取品字形分布房式采空區(qū)煤柱群中1個煤柱作為研究對象，煤柱上方的受載情況如圖1所示。

h—煤柱高度；H—開采深度；ω—煤柱截面邊長；b—煤柱間距即煤房寬度

結合圖1，考慮房式采空區(qū)煤柱埋深較淺，采出寬度較小，采空區(qū)矸石不承載且每個煤柱承載特性相同等因素，采用有效面積理論對元胞上覆均布載荷進行估算，煤柱的水平截面形狀及從屬面積形狀均為矩形，則煤柱元胞承受的均布載荷Lv可由式(1)表示[19]。

Lv=γH(ω+b)2/ω2

(1)

式中：γ為上覆巖層的平均容重。

1.2 元胞模型建立

采用8Moore元胞自動機模型將被研究的采空區(qū)煤柱沿其橫截面離散成15×15的矩形網格，每個網格單元均為1×1的正方形，代表1個元胞，每個元胞的位置用(xi，yi)(i，j=1，2,…，15)表示。由于每個元胞單元的邊長均為1，可用元胞吸收的廣義能量代表系統(tǒng)施加的廣義應力，且二者在量綱上保持統(tǒng)一。以E(xi，yi，t)代表各位置元胞單元t時刻存儲的廣義能量，等效于該元胞此刻存儲的彈性應變能，則E(xi，yi，0)代表元胞初始廣義能量。本模型將煤柱模型抽象為受損和完好2個部分，且規(guī)定已損部分不再具有承擔載荷的能力，全部載荷將由未受損的部分承擔。

1.3 元胞儲能極限

假設煤柱為各向同性材料，除四角元胞外，位于模型內部的元胞均受到來自周圍元胞群的約束作用，這種約束作用來自相連元胞之間的黏聚力以及四周元胞群上下表面所受摩擦力，定義元胞儲能極限為使周圍元胞群約束作用喪失時的元胞能量。在模型中，將這種約束作用等效于圍壓，不同圍壓對應不同的儲能極限，元胞所受的圍壓狀態(tài)及大小與元胞所處位置有關，由此定義二維元胞模型儲能極限如圖2所示。

圖2 元胞儲能極限分布

由于模型采用15×15正方形元胞矩陣，所以整體元胞儲能極限沿對角線對稱分布，因此僅分析圖2中黑色線框區(qū)域。圖2中，一級元胞：位于模型四角的4個元胞，處于完全單軸壓縮狀態(tài)，所以儲能極限最小，表示為U1；二級元胞：位于邊緣一周除四角外的元胞，近似處于單軸壓縮狀態(tài)，但每個元胞兩邊存在較小約束，儲能極限大于一級元胞，表示為U2；三級元胞：位于模型對角線上的元胞單元，處于等圍壓狀態(tài)，研究表明，隨著圍壓的增高，巖石的強度及儲能極限大幅增長，總體呈冪指數(shù)形式[20]。本模型隨著圍壓的逐漸增大，對應三級元胞的儲能極限也隨之呈冪指數(shù)增加，可由式(2)表示。

U3(xi,yi)=U1+σik

(2)

式中：σi為各位置三級元胞破壞所需的極限圍壓，(i=2，3，…，8)；k為與煤巖物理力學性質相關的常數(shù)。

對于四級元胞：處于雙向不等圍壓下其他位置的元胞，王云飛等[21]通過室內試驗及顆粒流與Fish語言程序對煤巖進行不同圍壓組合試驗，得出在最小主應力一定的情況下，隨著中間主應力的增加，煤巖的彈性儲能極限也隨之線性增加，則在元胞模型中，在每一圈的四級元胞向中間移進的過程中，各四級元胞破壞所需的圍壓比n(中間主應力σ2/最小主應力σ3)隨著增加，對應的儲能極限也隨之線性增加，可由式(3)表示。

U4(xi,yj)=U3(xi,yi)+kn

(3)

式中：U4(xi，yj)為各位置四級元胞儲能極限(i≠j≠1)；U3(xi，yi)為各位置三級元胞儲能極限；n為四級元胞破壞的極限圍壓比。

三級元胞中位于模型中心的元胞處于最大等圍壓狀態(tài)，且所受的約束力最大，儲能極限值也最大。根據(jù)模型的對稱性及儲能極限分布規(guī)律得到模型中所有元胞的破壞極限閾值U(xi，yi)，(i，j=1，2，3，…，15)。

1.4 元胞模型能量輸入及傳遞規(guī)律

煤柱上方為均布載荷，因此系統(tǒng)廣義能量的輸入方式為恒定加載，每一時步對系統(tǒng)所有未破損元胞施加相同的廣義能量，模擬煤柱上方的恒定均布載荷，隨后元胞系統(tǒng)開始自動進行演化，在任意時步t內，計算該時步內每個元胞所存儲的廣義應變能E(xi，yi，t)若滿足式(4)

E(xi,yi,t)≥U(xi,yi)

(4)

則認為該元胞破壞且不再具有承載能力，破壞元胞的一部分能量以等比例、等概率的形式傳遞給周圍未破壞的鄰居元胞，剩余的能量根據(jù)熱力學第一定律以一定的能量耗散率損耗，系統(tǒng)滿足能量守恒定律。

所有元胞按照上述規(guī)則進行演化，形成雪崩式的能量傳遞及元胞破損，直到某一時步t，輸入廣義能量后，系統(tǒng)內剩余各元胞的能量均小于其對應的儲能極限閾值時，演化終止，結合破壞統(tǒng)計量D表示元胞模型的破壞程度

D=At/A

(5)

其中：At為演化終止時已破壞元胞的數(shù)量，即煤柱橫截面上宏觀破壞的表征；A為系統(tǒng)內元胞總數(shù)量(15×15)。當損傷量D達到一定值時，認為煤柱失穩(wěn)。

2 演化結果分析

由于模型具有對稱性，演化結果分析依舊以圖2中黑色線框內的元胞為主。為便于描述演化過程及結果，假設所有元胞的初始能量均為0，并規(guī)定每個元胞破壞后向周圍各元胞傳遞的能量為定值。由于模型破壞演化的形式主要取決于不同位置元胞所受約束力大小，若假定每個元胞在x、y軸方向上對其他元胞所提供的約束力均為σt，且隨著軸向元胞數(shù)量增多，約束力線性疊加。對于對角線上的元胞，隨著x的增加，對應元胞所受的最大約束(極限圍壓)也隨之線性增加，如式(6)所示。

σi=xσt

(6)

式中：σi為位于(xi，yi)處的三級元胞的最大約束(極限圍壓)。

而對于同圈的四級元胞，其圍壓比也會隨著x近似線性變化，如式(7)所示。

n=k1x

(7)

式中：k1為煤巖單元圍壓變化系數(shù)，隨著y增加逐漸減小。

將式(6)及式(7)分別代入式(2)和式(3)中，可得式(8)，即

(8)

將式(8)統(tǒng)一到同一坐標中，其中橫坐標為煤柱二維元胞模型的x軸，縱坐標為各元胞儲能極限值U，如圖3所示。

圖3 Ⅰ型演化元胞儲能極限變化

其中，一級元胞儲能極限最低，在施加廣義能量后最先破壞，隨后二級元胞逐漸全部破壞，而三級與四級元胞的破壞有2種情況，對應以下2種演化模式。

Ⅰ型演化：由圖3可知，每圈所有四級元胞儲能極限均小于其緊鄰內圈三級元胞的儲能極限，在能量加載過程中，儲能極限閾值較小的元胞容易先行破壞。如圖4所示，由于對角三級元胞的儲能極限均大于其緊鄰外圈四級元胞的極限閾值，且假設的每個元胞破壞后向周圍元胞傳遞能量的值相同，在演化過程中只有當外圈四級元胞先開始破壞或全部破壞時，緊鄰內圈的三級元胞才開始破壞。而在同一圈中只有當三級元胞破壞后，四級元胞才開始破壞，因此煤柱的破壞過程呈逐層演化形式，其彈性核區(qū)逐漸向矩形轉化，演化特性如圖4b、圖4c、圖4d所示，并在圖4b、圖4c、圖4d之間不斷轉化；隨著演化逐漸轉向內部，三級元胞的儲能極限不斷增加，直至某一步時演化停止，煤柱趨于穩(wěn)定，最終的彈性核區(qū)如圖4d所示呈矩形。

圖4 Ⅰ型演化過程

演化初期元胞的損傷較劇烈，損傷當量D隨著時步先急劇增加后逐漸平緩，代表煤柱逐漸趨于穩(wěn)定。由于初期破壞主要為一級元胞與二級元胞，且二級元胞儲能極限相同，邊緣煤柱破壞較快，因此初期演化破壞較劇烈，損傷當量D也隨之急劇增加。當邊緣元胞全部破壞后，演化向內部轉移，破壞主要以四級元胞為主，三級元胞破壞較少。

Ⅱ型演化如圖5所示，每圈的四級元胞中均存在1個與其緊鄰內圈三級元胞儲能極限相等的元胞，稱之為等效元胞，位于等效元胞前的四級元胞儲能極限小于該三級元胞儲能極限，之后的四級元胞極限值大于該三級元胞儲能極限值。

圖5 Ⅱ型演化元胞儲能極限變化

Ⅱ型演化元胞的破壞特征同樣服從極限值較小的元胞先破壞的規(guī)律，演化也是先由一級元胞開始破壞，隨后二級元胞逐漸全部破壞，并誘發(fā)緊鄰內圈的三級元胞及四級元胞逐層破壞。元胞破壞演化特征如圖6a、圖6b、圖6c所示，與Ⅰ型演化初期(一級、二級及可能的三級和四級元胞破壞)的演化過程相一致，但由于Ⅱ型演化每圈四級元胞中等效元胞的存在，位于等效元胞之后的四級元胞儲能極限值大于該元胞儲能極限值，且四級元胞儲能能力由外向里逐漸遞增，因此Ⅱ型演化四級元胞整體儲能極限要大于Ⅰ型演化。當2種演化中三級元胞儲能極限及外部施加廣義能量相同時，Ⅱ型演化四級元胞的損傷破壞范圍較Ⅰ型演化小，則演化終止時Ⅱ型演化的彈性核區(qū)形態(tài)特征如圖6d所示近似呈四周凸起的多邊形。

圖6 Ⅱ型演化過程

由此可見，Ⅱ型演化的演化過程中也是初期損傷較多，損傷當量先急劇增加，后逐漸趨于緩和。但由于Ⅱ型演化四級元胞破壞范圍較小，最終穩(wěn)定后彈性核區(qū)總面積要大于Ⅰ型演化，因此穩(wěn)定性相比Ⅰ型演化更好。

3 數(shù)值模擬

為研究房式采空區(qū)下單個煤柱破壞演化的特性，結合石圪臺煤礦3-1-1煤層房式開采現(xiàn)狀，建立FLAC3D數(shù)值模擬模型，模型外界尺寸為10 m×10 m×38 m，煤柱尺寸為6 m×6 m×2 m，如圖7所示。

圖7 煤柱數(shù)值計算模型

結合現(xiàn)場現(xiàn)狀及巖石力學試驗得出巖層各分層的巖性，見表1[10]。

表1 煤巖力學參數(shù)

模型下表面固支，頂?shù)装鍘r層四周表面施加法向位移約束，煤柱四周無約束，模型上表面施加1.0 MPa垂直應力，側壓系數(shù)為1.2。塑性區(qū)演化結果如圖8所示。

由圖8a可知，由于煤柱的邊緣處于單向受力狀態(tài)，所以煤柱4個角先出現(xiàn)剪切破壞，隨后應力向煤柱內部轉移，使邊緣煤柱單元依次塑性破壞，并致使緊鄰內圈煤柱的對角單元破壞，如圖8b所示一級元胞及二級元胞全部破壞。由圖8c可知，在圖8b向圖8c過渡期間，對角單元不再產生破壞，而其緊鄰外圈單元從兩端開始向里逐漸破壞，這與元胞模型Ⅱ型演化圖6c對應；由圖8d可知，模擬計算至平衡，煤柱的彈性核區(qū)形態(tài)呈四周凸起的多邊形，與Ⅱ型演化最終彈性核區(qū)形態(tài)特征相對應。

圖8 數(shù)值計算模擬結果

數(shù)值模擬結果表明：煤柱的破壞主要為剪切破壞，且同樣由外向里漸進式破壞，其彈性核區(qū)分布變化規(guī)律與元胞模型的Ⅱ型演化一致，由此揭示房式采空區(qū)煤柱的穩(wěn)定演化服從Ⅱ型演化，同時驗證了煤柱元胞自動機模型的可靠性。

4 現(xiàn)場實例

烏蘭集團石圪臺煤礦3-1-1號煤層埋深較淺，為確保生態(tài)環(huán)境不遭到破壞，采用房式采煤法開采3-1-1號煤層，煤柱及煤房寬度均為6 m，經過多年后形成房式采空區(qū)，遺留大量煤柱，通過對現(xiàn)場1435房采工作面采空區(qū)進行表土層剝離揭露，揭露后處于穩(wěn)定狀態(tài)下的煤柱如圖9所示。

圖9 房式采空區(qū)表土層剝離試驗

由圖9a可知，揭露后的房式煤柱呈“品”字形分布，煤柱均處于穩(wěn)定狀態(tài)，且煤柱尺寸基本相同。選取其中穩(wěn)定性較好的煤柱(在圖9a中紅色線框標出)并放大如圖9b所示，此時煤柱四周存在塑性區(qū)并有少量煤塊脫落，表明一級元胞及二級元胞已全部破壞，煤柱水平截面彈性核區(qū)的形態(tài)特征與Ⅱ型演化結論及數(shù)值模擬結果一致，大致呈四周凸起的多邊形(圖9b中紅色虛線所示)。表明提出的房式煤柱穩(wěn)定性元胞演化模型及數(shù)值模擬結果的可靠性強，可為后續(xù)房式采空區(qū)煤柱穩(wěn)定性分析提供新的研究思路。

5 結 論

1)基于房式采空區(qū)承載特性，采用元胞自動推理機理論，建立了二維煤柱元胞模型，結合圍壓對煤巖儲能極限影響規(guī)律定義了元胞模型的儲能極限分布規(guī)律，推演出房式煤柱的2種演化模式，其中Ⅰ型演化最終彈性核區(qū)近似呈矩形，Ⅱ型演化最終彈性核區(qū)近似呈四周凸起的多邊形。

2)Ⅱ型演化模式每圈的四級元胞中存在與緊鄰內圈三級元胞儲能極限相同的等效元胞，其整體元胞儲能極限高于Ⅰ型演化模式，在同樣演化規(guī)則及外部施加相同廣義能量的條件下，Ⅱ型演化模式的穩(wěn)定性相對較好。

3)采用FLAC3D模擬了煤柱、煤房寬度均為6 m時的淺埋房式采空區(qū)煤柱，煤柱整體穩(wěn)定性較好，煤柱塑性區(qū)及彈性核區(qū)的演化規(guī)律與Ⅱ型演化相一致。

4)通過對石圪臺煤礦3-1-1號煤層房式采空區(qū)煤柱進行揭露試驗得出，煤柱均處于穩(wěn)定狀態(tài)，且煤柱水平截面彈性核區(qū)形態(tài)特征與Ⅱ型演化相對應，由此驗證了房式煤柱元胞演化模型及數(shù)值模擬結果的正確性。