考慮礦柱流變效應(yīng)的采空區(qū)力學(xué)模型分析

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

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

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

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

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

1.1 工程背景

某鎢礦采用房柱法開采。由于多年的開采和早年開采不規(guī)范,已經(jīng)形成了大大小小的采空區(qū),并且出現(xiàn)了礦柱片落、劈裂和頂板冒頂?shù)默F(xiàn)象。迫切需要對某鎢礦2線0沿東采空區(qū)進(jìn)行長期穩(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。礦柱設(shè)計為5 m×5 m方柱,礦房跨度10 m。

1.2 影響因素敏感性分析

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

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

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

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

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

2.1 礦柱蠕變模型

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

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

Burgers蠕變模型本構(gòu)方程:

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

化簡有:

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

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

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

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

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

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

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

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

構(gòu)造近似解的形式,設(shè)為

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

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

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

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

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

階段二的撓度方程形式設(shè)為

簡支邊條件下的分布函數(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)時間解析解推導(dǎo)

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

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

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

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

問題轉(zhuǎn)化成解此微分方程的解,式（13）解得:

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

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

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

式中,σ為礦柱中應(yīng)力,MPa。

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

階段一的破壞條件為

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

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

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

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

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

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

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

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

3.1 應(yīng)力分析

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

圖8 開采完豎直應(yīng)力云圖Fig.8 Vertical stress nephogram after mining

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

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

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

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

3.2 位移分析

由圖12知,模型豎直位移云圖上,礦房開采引起的應(yīng)力變化相互疊加,導(dǎo)致周圍巖體產(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)測點(diǎn)和5#監(jiān)測點(diǎn)、2#監(jiān)測點(diǎn)和4#監(jiān)測點(diǎn)符合對稱的結(jié)果。從整體上來看,采場中心的撓度到采場邊緣呈現(xiàn)遞減的現(xiàn)象,符合之前假設(shè)的勢函數(shù)形式。5個監(jiān)測點(diǎn)在開采完27.70個月內(nèi),礦巖在上覆巖層持續(xù)應(yīng)力作用下顯現(xiàn)出顯著的時間效應(yīng),先變形速率快,礦柱的頂部位移一直呈上升趨勢,后變形速率減小,曲線趨于平穩(wěn)。最后在118.06個月時,模型停止計算。

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

4 結(jié) 論

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

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

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