考慮礦柱流變效應的采空區力學模型分析

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

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

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

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

1 采空區影響因素敏感性分析

1.1 工程背景

某鎢礦采用房柱法開采。由于多年的開采和早年開采不規范,已經形成了大大小小的采空區,并且出現了礦柱片落、劈裂和頂板冒頂的現象。迫切需要對某鎢礦2線0沿東采空區進行長期穩定性分析,采空區示意圖如圖1所示。

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

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

1.2 影響因素敏感性分析

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

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

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

2 采空區長期穩定性力學分析

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

2.1 礦柱蠕變模型

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

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

Burgers蠕變模型本構方程:

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

化簡有:

2.2 采空區力學模型

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

圖3 采空區簡化模型Fig.3 simplified model of goaf

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

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

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

聯立式（1）、式（3）消除式中σ,得到采空區頂板—礦柱系統的控制方程。

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

構造近似解的形式,設為

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

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

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

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

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

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

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

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

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

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

2.3 采空區自穩時間解析解推導

采空區的失穩是一個逐漸發展的過程,自穩時間可以看作階段一、階段二的時間總和。

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

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

表2 礦柱流變參數Table 2 Rheological parameters of pillars

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

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

式中,A1、A2為待定系數,根據邊界條件確定。

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

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

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

階段一的破壞條件為

解得階段一采空區自穩時間為32.14個月。

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

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

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

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

3 采空區長期穩定性數值分析

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

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

3.1 應力分析

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

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

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

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

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

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

3.2 位移分析

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

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

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

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

4 結 論

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

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

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