王 震，任高峰，張聰瑞，張健峰

（1.武漢理工大學(xué) 資源與環(huán)境工程學(xué)院，湖北 武漢 430070；2.礦物資源加工與環(huán)境湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430070）

地下開采形成的采空區(qū)在我國廣泛、大體量存在，已成為礦區(qū)的重大安全隱患。采空區(qū)頂板長時(shí)間大面積懸露，在動(dòng)力擾動(dòng)、次生應(yīng)力場、巖體儲能條件等組合因素影響下易發(fā)生失穩(wěn)斷裂，將采空區(qū)的空氣瞬間壓縮，形成巨大的沖擊氣浪，造成極大的破壞[1-4]。因此，采空區(qū)失穩(wěn)誘發(fā)的動(dòng)力學(xué)災(zāi)害是確保礦區(qū)安全以及深部開采穩(wěn)步實(shí)施過程中必須解決的問題。針對動(dòng)力擾動(dòng)下采空區(qū)圍巖的力學(xué)行為，國內(nèi)外的學(xué)者做了大量的研究工作。Zou C[5]等研究了單裂紋石膏礦試樣在動(dòng)載荷和準(zhǔn)靜載荷作用下的力學(xué)性質(zhì)和破裂行為，研究發(fā)現(xiàn)應(yīng)變率越大，石膏礦抗壓強(qiáng)度越大，且在動(dòng)載荷作用下的抗壓強(qiáng)度高于靜載荷。周小平[6]根據(jù)損傷斷裂力學(xué)知識建立了巖石處于卸荷條件全過程應(yīng)力應(yīng)變關(guān)系。王軍[7]通過運(yùn)用FLAC3D數(shù)值模擬、理論分析等方法對不同構(gòu)造應(yīng)力場中采空區(qū)上覆巖層的破壞規(guī)律進(jìn)行了研究。王開等[8]建立了頂板采空區(qū)堅(jiān)硬頂板斷裂懸臂巖梁力學(xué)模型，推導(dǎo)了采空區(qū)上覆巖層載荷和控頂距關(guān)系式。Gao F等[9]通過數(shù)值方式，研究了采空區(qū)上方巖層的應(yīng)力分布以及與采空區(qū)相連巷道的動(dòng)態(tài)載荷，揭露了巷道的失穩(wěn)機(jī)理。魏明堯等[10]利用FLAC3D分析了動(dòng)力擾動(dòng)下頂板的穩(wěn)定性以及頂板失穩(wěn)的動(dòng)力學(xué)特性。針對采空區(qū)失穩(wěn)產(chǎn)生的沖擊氣浪災(zāi)害特性，任高峰等[11]將采空區(qū)失穩(wěn)簡化為打氣筒模型，利用高速攝像和紋影技術(shù)，記錄和分析了沖擊氣浪在不同巷道類型的傳播規(guī)律以及變化特性。錢兆明[12]通過結(jié)合采空區(qū)坍塌的空氣壓縮模型，利用相似實(shí)驗(yàn)，分析了采空區(qū)坍塌造成的高壓氣浪在巷道中的傳遞規(guī)律。然而，這些研究缺乏對采空區(qū)失穩(wěn)過程中時(shí)空能量規(guī)律的探索，對于失穩(wěn)的演化過程及誘發(fā)多災(zāi)種連鎖反應(yīng)也缺乏深入分析。針對采空區(qū)失穩(wěn)下衍生多災(zāi)種連鎖反應(yīng)現(xiàn)象，采用數(shù)值模擬等研究方法，分析采空區(qū)失穩(wěn)誘發(fā)多災(zāi)種連鎖反應(yīng)過程中頂板及沖擊氣浪災(zāi)害的動(dòng)力學(xué)行為特性，研究結(jié)果可供礦山的防災(zāi)減災(zāi)提供一定的理論支持。

1 采空區(qū)流固耦合建模

建立頂板整體切落的失穩(wěn)流固耦合模型時(shí)，將其類比于“打氣筒”模型，將頂板視為整體?；谝陨霞僭O(shè)，頂板下落過程中的受力為自身重力和頂板下部被壓縮的氣體的反向阻力，該阻力隨著頂板下落高度的增大而逐漸增大，采空區(qū)內(nèi)的氣體瞬間受壓，流速急劇變大，形成能量巨大的沖擊氣浪對采場內(nèi)人員以及設(shè)備造成極大危害[13]。采空區(qū)頂板整體切落簡化圖形如圖1。

圖1 采空區(qū)頂板整體切落簡化圖形

完全氣體假設(shè)下，其狀態(tài)方程滿足：

式中：p 為壓力，Pa；ρ為流體密度，kg/m3；R 為氣體常數(shù)，取值8.314 4 J/（mol·K-1）；T為熱力學(xué)溫度，K。

根據(jù)等壓比熱公式，采空區(qū)內(nèi)和采空區(qū)坍塌時(shí)巷道風(fēng)流入口的溫度比值滿足：

式中：T2為巷道接口處的氣體溫度，K；T1為采空區(qū)內(nèi)部溫度，K；k為空氣的比熱比；p0為采空區(qū)內(nèi)初始壓力，Pa；p1為采空區(qū)坍塌時(shí)的巷道接口處壓力，Pa。

得到?jīng)_擊氣浪離開坍塌中的采空區(qū)進(jìn)入巷道的流速：

式中：voutflow為巷道與采空區(qū)接口處的空氣流速，m/s；vRoof為頂板下落速度，m/s；pgoof為采空區(qū)內(nèi)的氣體壓力，Pa；poutlet為采空區(qū)與巷道接口處的氣體壓力，Pa。

2 采空區(qū)頂板整體失穩(wěn)災(zāi)害分析

選擇大柳塔礦活雞兔井1-2煤層頂板[14]作為切落模型的預(yù)測對象，采空區(qū)的幾何尺寸和巖層密度等參數(shù)按照該對象的實(shí)際工程情況設(shè)置。實(shí)際工程情況采空區(qū)參數(shù)見表1。

表1 實(shí)際工程情況采空區(qū)參數(shù)表

2.1 采空區(qū)風(fēng)速變化規(guī)律

根據(jù)流固耦合模型，利用Matlab計(jì)算頂板整體切落采空區(qū)失穩(wěn)過程中的風(fēng)速，繪制的采空區(qū)內(nèi)風(fēng)速隨頂板下落高度變化規(guī)律圖如圖2。

從圖2可以看出，采空區(qū)內(nèi)風(fēng)速的變化大致可分為2個(gè)階段：風(fēng)速增大階段和風(fēng)速下降階段。

1）采空區(qū)內(nèi)風(fēng)速增大階段（0～A）。頂板開始下落時(shí)加速度較大，隨后遞減，設(shè)頂板下落高度為Z，下落時(shí)間為t。從頂板開始下落Z=0 m，t=0 s開始算起，當(dāng)t=1.096 s時(shí)，空區(qū)的風(fēng)速從0 m/s增大到峰值 6.85 m/s，此時(shí)對應(yīng)的下落高度為 Z=4.307 7 m。

圖2 采空區(qū)風(fēng)速隨頂板下落高度變化規(guī)律

2）采空區(qū)內(nèi)風(fēng)速減小階段（A～B）。當(dāng)頂板繼續(xù)下落時(shí)，因采空區(qū)內(nèi)空氣阻力對頂板的作用增大，頂板開始做減速運(yùn)動(dòng)，采空區(qū)風(fēng)速也隨之開始減小。從頂板下落時(shí)間t=1.096 s算起，頂板下落速度減緩，采空區(qū)風(fēng)速開始減小，此后頂板繼續(xù)下落，直至觸底，此時(shí)頂板下落高度Z=5 m，采空區(qū)風(fēng)速降低至 6.38 m/s，總體下落時(shí)間 t=1.200 s，頂板下落過程結(jié)束。

2.2 采空區(qū)風(fēng)壓變化規(guī)律

利用Matlab計(jì)算頂板整體切落采空區(qū)失穩(wěn)過程中的風(fēng)壓，繪制的風(fēng)壓隨頂板下落高度變化規(guī)律圖如圖3。

圖3 采空區(qū)內(nèi)風(fēng)壓隨頂板下落高度變化規(guī)律

從圖3可以看出，采空區(qū)內(nèi)風(fēng)壓的變化趨勢并非是簡單的線性，可分為緩慢上升階段和指數(shù)上升階段。

1）緩慢上升階段（0～A）。從 t=0 s 到 t=1.052 s的下落過程中，采空區(qū)內(nèi)的壓強(qiáng)由初始狀態(tài)的大氣壓 0.101 MPa，升高到了 0.356 8 MPa。此過程中頂板的下落高度Z=4 m，這個(gè)階段下氣體壓強(qiáng)的增幅并不是十分明顯。

2）指數(shù)上升階段（A～B）。頂板從Z=4 m繼續(xù)下落，采空區(qū)內(nèi)風(fēng)壓繼續(xù)增大，當(dāng)頂板下落高度到達(dá)Z=5 m時(shí)，整個(gè)頂板坍塌過程結(jié)束，結(jié)束的時(shí)間t=1.200 s，采空區(qū)風(fēng)壓達(dá)到整個(gè)過程的最大值0.961 7 MPa，由此可見，第2階段頂板下落時(shí)間不到第1階段的1/5，然而風(fēng)壓的增幅卻是第1階段的2.36倍，這也展示了2階段的差別。

2.3 沖擊氣浪災(zāi)害分析

2.3.1 沖擊氣浪縱向流速演化規(guī)律

采用采空區(qū)頂板整體切落模型，利用伯努利能量方程，考慮沖擊氣浪在巷道中傳播的沿程損失，計(jì)算距離巷道入口處距離L處的沖擊氣浪流速隨頂板下落高度分布（圖4）。

圖4 巷道距離 L=10、30、50、80、120、180 m 處沖擊氣浪流速隨頂板下落高度分布圖

由此可見，即使距離巷道入口的距離已經(jīng)達(dá)到180 m，沖擊氣浪仍然以亞音速狀態(tài)傳播，仍然具有極大的破壞力。根據(jù)圖4可得巷道內(nèi)最大風(fēng)速隨巷道距離的增大呈負(fù)指數(shù)趨勢減小，對沖擊氣浪最大流速衰減規(guī)律進(jìn)行擬合，得到擬合曲線（圖5）和擬合公式Vmax（L）=819.2e-0.6902×10-2L

圖5 最大流速衰減擬合曲線

2.3.2 沖擊氣浪橫截面流速分布規(guī)律

FLUENT模擬之前，選擇巷道各橫截面處的測點(diǎn)用于數(shù)據(jù)分析，按照從巷道出口到巷道與采空區(qū)相連入口，共70 m，每5 m選取橫截面。巷道出口為橫截面編號0，依次遞增至編號14。各橫截面上按照對角線的四等分線取點(diǎn)，設(shè)取點(diǎn)的標(biāo)號在截面0上分別為p1～p9，截面1上分為p10～p18，以此類推。取點(diǎn)分布如圖6。各點(diǎn)之間的相對坐標(biāo)見表2。

圖6 與采空區(qū)相連巷道FLUENT模擬各截面取點(diǎn)分布圖

表2 巷道各截面相對坐標(biāo)

利用FLUENT數(shù)值模擬在不同的時(shí)間間隔下取巷道內(nèi)的沖擊氣浪流場圖。采空區(qū)頂板垮落t=0.1、0.2、0.3、0.4、0.5 s 巷道內(nèi)沖擊氣浪流速分布規(guī)律如圖7。圖中下方長方形區(qū)域代表的是采空區(qū)，中間連接部分代表的是與采空區(qū)相連的巷道，在模擬時(shí)設(shè)置為70 m。上方面積略小的長方形代表的是從巷道中排出的沖擊氣浪的作用區(qū)域。

圖7 采空區(qū)頂板垮落 t=0.1、0.2、0.3、0.4、0.5 s巷道內(nèi)沖擊氣浪流速分布規(guī)律圖

由圖 7 可知，頂板下落時(shí)間在 t=0.1、0.2、0.3、0.4、0.5 s 時(shí)，沖擊氣浪最大值分別為 4.5、22、60、130、260 m/s，大流速集中在采空區(qū)與巷道入口相連處。

采空區(qū)頂板垮落 t=0.6、0.7、0.8、0.9、1.0 s 巷道內(nèi)沖擊氣浪流速分布規(guī)律如圖8。由圖8可知，從t=0.6 s開始，巷道內(nèi)沖擊氣浪流速的分布較之前發(fā)生了變化，在距離采空區(qū)出口處較近的截面，靠近巷道左側(cè)壁相比與右側(cè)壁具有更高的流速，左側(cè)的最高流速達(dá)到了音速水平，右側(cè)的流速衰減至100 m/s以下，并且巷道內(nèi)的流速沿巷道縱向出現(xiàn)了左側(cè)高流速區(qū)和右側(cè)低流速區(qū)。

根據(jù)圖6取點(diǎn)方式，結(jié)合流速分區(qū)的現(xiàn)象，發(fā)現(xiàn)在巷道左側(cè)底角的沖擊氣浪流速最大，比較不同截面左側(cè)底角測點(diǎn)的流速，（截面0中為p9點(diǎn)，其他截面依次類推），并繪制了不同截面同一測點(diǎn)的流速分布圖（圖 9）。

通過圖9比較發(fā)現(xiàn)，靠近采空區(qū)與巷道入口處的p126點(diǎn)的沖擊氣浪流速最大，最高風(fēng)速達(dá)到了950 m/s，隨距離采空區(qū)的距離的增大，沖擊氣浪在傳播過程中產(chǎn)生沿程能量損失，最大流速逐步減小。

3 結(jié)論

1）采空區(qū)頂板災(zāi)害發(fā)生時(shí)，隨著頂板下落高度的增加，采空區(qū)內(nèi)風(fēng)速變化經(jīng)歷了加速以及減速階段。壓強(qiáng)變化分為緩慢上升階段和指數(shù)上升階段，其中緩慢上升階段占總時(shí)間比例大。

圖8 采空區(qū)頂板垮落 t=0.6、0.7、0.8、0.9、1.0 s巷道內(nèi)沖擊氣浪流速分布規(guī)律圖

圖9 巷道橫截面各測點(diǎn)的沖擊氣浪流速

2）沖擊氣浪災(zāi)害演化過程中，最大流速沿巷道軸向衰減與巷道距離大致為指數(shù)關(guān)系，沖擊氣浪流速在傳播過程中不斷下降，然而不同距離處沖擊氣浪隨時(shí)間的變化趨勢一致。

3）沖擊氣浪在巷道截面的二維分布有流速分區(qū)現(xiàn)象，流速沿巷道軸線可分為左側(cè)高流速區(qū)和右側(cè)低流速區(qū)，在模擬巷道中，各截面左側(cè)底角沖擊氣浪的流速分布不均，巷道入口左側(cè)底角處流速最大。