蔡 泳 鄧國梁 甄利兵 成 墻 黃金涌

（1.煤礦災害動力學與控制國家重點實驗室,重慶 400044;2.重慶大學資源與安全學院,重慶 400044;3.貴州紫金礦業(yè)股份有限公司,貴州 貞豐 562200）

采礦方法中,房柱法劃分礦塊為礦柱和礦房,回采礦房而預留下礦柱,形成大量的采空區(qū)[1]。如果不及時處理采空區(qū),可能造成人員上的傷亡和財產(chǎn)損失[2-5]。因此采空區(qū)自穩(wěn)時間的預估顯得尤其重要。在穩(wěn)定時間內,對采空區(qū)進行處理,如支護、充填等工作,可避免安全事故的發(fā)生。采空區(qū)力學分析能有效地預估采空區(qū)穩(wěn)定時間。多數(shù)力學分析將礦柱簡化為winkler彈性基礎和蠕變模型[6-8],將頂板簡化成彈性薄板、彈性巖梁或彈性厚板組建力學模型分析[9-11]。Zhao等[12]基于巖梁模型,確定了頂板巖塊轉角與撓度的關系。孫琦等[13]將采空區(qū)頂板簡化成彈性薄板,將礦柱簡化成Kelvin-Voigt模型,考慮水平應力和垂直應力的共同作用,建立起多向荷載作用下的采空區(qū)力學模型,經(jīng)過理論計算與現(xiàn)場實際監(jiān)測結果對比分析,驗證了模型的實用性。樓小明等[14]將礦柱簡化成Hooke-Kelvin-kelvin蠕變模型,頂板簡化為彈性薄板,建立力學模型應用到工程中,并指出礦柱的流變變形是采空區(qū)失穩(wěn)的主要原因。Wang等[15]基于Reissner厚板理論建立采空區(qū)力學模型,通過互等定理推導出在不同的邊界條件頂板撓度隨時間變化的表達式,借助數(shù)值模擬軟件驗證解析解的合理性。謝學斌等[16]將采空區(qū)頂板和礦柱簡化成彈性薄板和Burgers蠕變模型,構建出三維蠕變模型,借助尖點突變理論對采空區(qū)失穩(wěn)傾向做出定量分析,并分析了各個影響因子對系統(tǒng)穩(wěn)定性的影響。上述礦柱—頂板力學模型建立過程中,大多是將礦柱視為流變模型,頂板視為彈性模型。實際上,采空區(qū)自穩(wěn)時間是受礦柱和頂板共同流變作用的影響。若對采空區(qū)影響因素進行敏感性分析,對主控因素考慮其流變性作用,對非主控因素理想化考慮,會縮減力學模型計算量的同時,還能保證模型精確度。

鑒于此,本研究以某鎢礦2線0沿東采空區(qū)為工程背景,根據(jù)因素敏感性分析結果,構建采空區(qū)三維流變力學模型,得出采空區(qū)自穩(wěn)時間的表達式,并利用Flac3D對采空區(qū)蠕變過程中應力、位移進行了分析。

1 采空區(qū)影響因素敏感性分析

1.1 工程背景

某鎢礦采用房柱法開采。由于多年的開采和早年開采不規(guī)范,已經(jīng)形成了大大小小的采空區(qū),并且出現(xiàn)了礦柱片落、劈裂和頂板冒頂?shù)默F(xiàn)象。迫切需要對某鎢礦2線0沿東采空區(qū)進行長期穩(wěn)定性分析,采空區(qū)示意圖如圖1所示。

圖1 2線0沿東采空區(qū)示意Fig.1 Schematic diagram of goaf along E0 of line 2

取某鎢礦采場基本參數(shù),分析采場的長×寬為80 m×80 m,礦柱高度6 m,埋深200 m。礦柱設計為5 m×5 m方柱,礦房跨度10 m。

1.2 影響因素敏感性分析

參考李令鑫等[17]對采空區(qū)安全系數(shù)的定義,將采空區(qū)的影響確定為6因素5水平,進行正交極差分析,所得結果如表1所示。

表1 各因素極差分析結果Table 1 Range analysis results of various factors

通過極差分析可得各個影響因素的主次關系:礦柱的寬高比＞礦柱的單軸抗壓強度＞頂板上覆巖層的容重＞采空區(qū)埋深＞礦房跨度＞頂板厚度。

2 采空區(qū)長期穩(wěn)定性力學分析

采空區(qū)的穩(wěn)定與其結構參數(shù)和時間有關,考慮時間效應的作用下,礦柱能承受最大應力值應為巖石的長期強度[18]。根據(jù)因素敏感性分析結果,在采空區(qū)自穩(wěn)時間的力學模型構建中,礦柱若簡單視為彈性元件,可能會對結果造成較大的偏差。同時,頂板厚度影響性最小,力學模型的構建中可以將頂板簡化為彈性薄板。

2.1 礦柱蠕變模型

經(jīng)過文獻[19]研究,對比不同應力下的蠕變實驗曲線與蠕變模型發(fā)現(xiàn)Burgers模型能夠較好表征硬質礦柱的蠕變行為。Burgers流變模型如圖2所示。

圖2 巖石Burgers流變模型Fig.2 Burgers creep model of rock

Burgers蠕變模型本構方程:

式中,η1、η2為牛頓黏性系數(shù);k1、k2為彈性元件的彈性系數(shù);σ為模型的總應力,MPa;ε為模型的總應變。

化簡有:

2.2 采空區(qū)力學模型

根據(jù)因素敏感性分析,采空區(qū)的力學模型可以簡化礦柱為Burgers模型,而頂板可以簡化成彈性薄板。簡化后的采空區(qū)力學模型如圖3所示。

圖3 采空區(qū)簡化模型Fig.3 simplified model of goaf

隨著開采的進行,頂板邊界逐漸破壞,采空區(qū)的邊界條件分為3個階段:固支邊、簡支邊和自由邊。

（1）階段一:固支邊位移方程。當頂板處于固支邊時,則頂板的控制方程為:

式中,D為頂板抗彎剛度,kN/m2;w為頂板撓度,m;ξ為將礦柱中的集中力均布化系數(shù);σ為礦柱中應力,MPa;q為上覆巖層作用在頂板上的壓力,MPa;n為采場內礦柱個數(shù);A為礦柱的橫截面積,m2;a、b為矩形采場的長邊和短邊,m;E為頂板巖石彈性模量,GPa;h為頂板厚度,m。

聯(lián)立式（1）、式（3）消除式中σ,得到采空區(qū)頂板—礦柱系統(tǒng)的控制方程。

對式（4）,采用伽遼金法求其弱解。階段一時。頂板邊界處于固定,撓度和轉角為0:

構造近似解的形式,設為

式中,w0（t）為頂板中心點隨時間變化的撓度,m;φ（x,y）為撓度隨位置變化的分布函數(shù)。當t=0時,固支邊條件下采場的撓度如圖4所示。

圖4 t=0固支邊時采場撓度Fig.4 Stope deflection with fixed support when t=0

近似解滿足邊界條件,加權余數(shù)表示為

式中,w0為采場中心撓度,m;J1、J2、J3為化簡后系數(shù)。

（2）階段二:簡支邊位移方程。階段一采場位移達到破壞條件,此時,頂板邊界出現(xiàn)了轉角,處于簡支邊,邊界的撓度和彎矩為0。

階段二的撓度方程形式設為

簡支邊條件下的分布函數(shù)φ（x,y）如圖5所示。

圖5 t=0簡支邊時采場撓度Fig.5 Stope deflection with simply supported edge when t=0

同理,階段二采場中心撓度的微分方程為

頂板繼續(xù)持續(xù)破壞,表面裂隙形成并貫通,頂板已經(jīng)完全失去了承載能力（D=0）。頂板的邊界處于自由邊,階段三狀態(tài)。

2.3 采空區(qū)自穩(wěn)時間解析解推導

采空區(qū)的失穩(wěn)是一個逐漸發(fā)展的過程,自穩(wěn)時間可以看作階段一、階段二的時間總和。

采空區(qū)處在階段一時,采場的撓度表示為

表達式滿足邊界條件,代入某鎢礦基本數(shù)據(jù)和礦柱流變數(shù)據(jù)（表2）并根據(jù)式（7）有:

表2 礦柱流變參數(shù)Table 2 Rheological parameters of pillars

問題轉化成解此微分方程的解,式（13）解得:

接收現(xiàn)場數(shù)據(jù)，并對數(shù)據(jù)進行分析，根據(jù)特定的計算方式給出相應的指令，這一部分稱之為運算部分。邏輯控制部分：與現(xiàn)場設備相連接，接收運算部分的動作指令，判斷現(xiàn)場設備的實際運行狀況并將新的動作指令分別下發(fā)給各終端，使其按照要求單獨動作或者聯(lián)動動作。這2個部分互有區(qū)別，不應混為一談。

式中,A1、A2為待定系數(shù),根據(jù)邊界條件確定。

t=0時,根據(jù)前人的研究[20]和Burgers蠕變方程,采場礦柱剛受壓時,有初始位移和初始蠕變速度。

式中,σ為礦柱中應力,MPa。

求解出待定系數(shù),最終得出階段一時采場中心撓度隨時間的表達式:

階段一的破壞條件為

解得階段一采空區(qū)自穩(wěn)時間為32.14個月。

同理,也可以得到階段二頂板中心撓度的表達式和破壞表達式。

式中,[σt]為允許抗拉強度,MPa,ν為泊松比。

代入數(shù)據(jù),階段二的破壞時間為88.03個月。則采空區(qū)總的自穩(wěn)時間120.22個月。根據(jù)階段一和階段二的表達式,頂板中心點的撓度隨時間變化的趨勢如圖6所示。

圖6 頂板中心撓度隨時間變化的曲線Fig.6 Curve of central deflection of roof with time

3 采空區(qū)長期穩(wěn)定性數(shù)值分析

Flac3D是ltasca公司研發(fā)的連續(xù)介質力學分析軟件[21],被大量用于巖土工程和采礦工程中。建立計算模型如圖7所示,礦柱設為Burgers本構模型,頂板和底板為Mohr-Coulomb模型。設立5個礦柱頂部監(jiān)測點,觀察采空區(qū)失穩(wěn)的破壞規(guī)律。蠕變計算一步設為150 h,計算至失穩(wěn)后停止計算。

圖7 Flac3D計算模型示意Fig.7 Schematic diagram of Flac3D model calculation

3.1 應力分析

開挖礦體對地下初始應力場造成擾動,導致地下的應力場重新分布。由圖8和圖9知,模型豎直應力云圖上,豎直應力整體上呈現(xiàn)隨著深度的增加而增加。采空區(qū)形成后,礦體采空部分形成了應力釋放,礦柱中心部分有明顯的應力集中。礦柱頂端的壓應力為11 MPa接近理論計算的應力值。礦柱上部形成明顯的拱形等值線,礦柱頂板中間出現(xiàn)了較小的拉應力。而隨著蠕變迭代步數(shù)的進行,觀察到礦柱上方的豎直拉應力和影響區(qū)域都在不斷增加。

圖8 開采完豎直應力云圖Fig.8 Vertical stress nephogram after mining

圖9 蠕變過程中豎直應力云圖Fig.9 Vertical stress nephogram in creep process

由圖10知,模型最大主應力云圖上,采空區(qū)頂板上主要分布拉應力,最大拉應力出現(xiàn)在2個礦柱之間的中心的頂板上,未達到巖石的抗拉強度,此時頂板關鍵層并不會發(fā)生破壞。單獨通過拉應力值的大小判斷頂板是否會發(fā)生破斷并不合理,還需要結合頂板下沉量綜合分析。

由圖11知,模型最小主應力云圖上,最小主應力基本上呈層狀分布,礦柱支撐位置出現(xiàn)了應力集中,礦柱支撐力分布表現(xiàn)出中間礦柱支撐力大,而邊緣礦柱支撐力小的現(xiàn)象。

圖11 蠕變過程中最小主應力云圖Fig.11 Minimum principal stress nephogram in creep process

3.2 位移分析

由圖12知,模型豎直位移云圖上,礦房開采引起的應力變化相互疊加,導致周圍巖體產(chǎn)生的位移也隨之變化。此時,采空區(qū)頂板最大下沉量為18.93 mm,小于臨界破壞條件21.30 mm。而采空區(qū)底板出現(xiàn)了正向位移,最大位移量為2.59 mm。在采空區(qū)的上部形成了規(guī)律的位移等值線,并隨著埋深的增加而等值線的大小逐漸減小。礦柱的位移規(guī)律為:礦柱上半部表現(xiàn)為向下的位移,而下半部為向上的位移,表現(xiàn)出兩端受壓的狀態(tài)。

圖12 蠕變過程豎直位移云圖Fig.12 Vertical displacement nephogram of creep process

對5個礦柱頂部中心位置監(jiān)測位移,由圖13知,1#監(jiān)測點和5#監(jiān)測點、2#監(jiān)測點和4#監(jiān)測點符合對稱的結果。從整體上來看,采場中心的撓度到采場邊緣呈現(xiàn)遞減的現(xiàn)象,符合之前假設的勢函數(shù)形式。5個監(jiān)測點在開采完27.70個月內,礦巖在上覆巖層持續(xù)應力作用下顯現(xiàn)出顯著的時間效應,先變形速率快,礦柱的頂部位移一直呈上升趨勢,后變形速率減小,曲線趨于平穩(wěn)。最后在118.06個月時,模型停止計算。

圖13 5個監(jiān)測點豎直位移隨時間的曲線Fig.13 Curve of vertical displacement of five monitoring points with time

4 結 論

（1）采用敏感性分析方法得出了各個影響因素的敏感性,根據(jù)敏感性分析結果將礦柱簡化Burgers蠕變模型,頂板簡化為彈性薄板,建立了采空區(qū)三維蠕變力學模型,結合Burgers蠕變方程,推導出頂板位移隨時間變化的表達式。

（2）根據(jù)新構建的力學模型對某鎢礦采空區(qū)自穩(wěn)時間進行了預測,與數(shù)值模擬軟件對比分析,驗證了模型的合理性和可靠性。

（3）利用Flac3D軟件對某鎢礦采場模擬分析,隨著蠕變迭代步數(shù)的進行,觀察到礦柱上方的豎直拉應力和影響區(qū)域都在不斷增加,頂板的位移變化速率呈先增大后減小的趨勢。