蔡 泳 鄧國(guó)梁 甄利兵 成 墻 黃金涌

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

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

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

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

1.1 工程背景

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

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

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

1.2 影響因素敏感性分析

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

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

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

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

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

2.1 礦柱蠕變模型

經(jīng)過文獻(xiàn)[19]研究,對(duì)比不同應(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)變。

化簡(jiǎn)有:

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

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

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

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

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

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

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

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

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

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

圖4 t=0固支邊時(shí)采場(chǎng)撓度Fig.4 Stope deflection with fixed support when t=0

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

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

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

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

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

圖5 t=0簡(jiǎn)支邊時(shí)采場(chǎng)撓度Fig.5 Stope deflection with simply supported edge when t=0

同理,階段二采場(chǎng)中心撓度的微分方程為

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

2.3 采空區(qū)自穩(wěn)時(shí)間解析解推導(dǎo)

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

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

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

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

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

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

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

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

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

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

階段一的破壞條件為

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

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

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

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

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

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

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

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

3.1 應(yīng)力分析

開挖礦體對(duì)地下初始應(yīng)力場(chǎng)造成擾動(dòng),導(dǎo)致地下的應(yīng)力場(chǎng)重新分布。由圖8和圖9知,模型豎直應(yīng)力云圖上,豎直應(yīng)力整體上呈現(xiàn)隨著深度的增加而增加。采空區(qū)形成后,礦體采空部分形成了應(yīng)力釋放,礦柱中心部分有明顯的應(yīng)力集中。礦柱頂端的壓應(yīng)力為11 MPa接近理論計(jì)算的應(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個(gè)礦柱之間的中心的頂板上,未達(dá)到巖石的抗拉強(qiáng)度,此時(shí)頂板關(guān)鍵層并不會(huì)發(fā)生破壞。單獨(dú)通過拉應(yīng)力值的大小判斷頂板是否會(huì)發(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)生的位移也隨之變化。此時(shí),采空區(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

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

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

4 結(jié) 論

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

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

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