矩形巷道礦震應(yīng)力波擴(kuò)散規(guī)律譜元法分析*

唐杰兵,鞠文君,焦建康,陳法兵

(1.煤炭科學(xué)研究總院 開采研究分院，北京 100013；2.中煤科工開采研究院有限公司，北京 100013；3.天地科技股份有限公司 開采設(shè)計(jì)事業(yè)部，北京 100013)

0 引言

礦震應(yīng)力波觸發(fā)巷道沖擊是一種常見的現(xiàn)象[1]，研究礦震應(yīng)力波傳播特性的方法較多，由于應(yīng)力波在介質(zhì)中傳播是一個(gè)入射、折射、反射、散射及它們之間相互干擾等的復(fù)雜過程，當(dāng)傳播過程中遇到異質(zhì)性巖體，應(yīng)力波會發(fā)生疊加、產(chǎn)生角波、P波轉(zhuǎn)化為S波以及振幅增大等復(fù)雜現(xiàn)象[2]，在理論分析上，采用數(shù)學(xué)、力學(xué)等學(xué)科較難建立相應(yīng)的沖擊地壓復(fù)雜現(xiàn)象的數(shù)學(xué)模型，較難捕捉到應(yīng)力波傳播過程特征[3-4]。實(shí)驗(yàn)室試驗(yàn)也是如此，霍普金森桿試驗(yàn)研究動載應(yīng)力波的比較多[5]，也有相似材料模擬的動載應(yīng)力波研究[6-7](震源激發(fā)采用爆炸物或擺錘)，其無法真實(shí)模擬原巖應(yīng)力條件和介質(zhì)的力學(xué)環(huán)境，或離心試驗(yàn)?zāi)Ｐ蚚8](人造重力場)框架周邊橡膠墊層只能減輕邊界效應(yīng)，所以這些現(xiàn)場或?qū)嶒?yàn)室試驗(yàn)無法解決應(yīng)力波的邊界影響問題，且不便進(jìn)行參數(shù)敏感性分析，現(xiàn)場或?qū)嶒?yàn)室試驗(yàn)中只能設(shè)置有限的監(jiān)測點(diǎn)，震源模擬的可靠性缺乏可依靠的研究。在定量化研究中，主要的研究手段還是數(shù)值模擬方法，可用于應(yīng)力波傳播的數(shù)值模擬軟件：FALC2D[9]，UDEC[10]，PFC[11]，ABAQUS[12]等，可以合理地設(shè)置應(yīng)力波的吸收邊界條件及數(shù)據(jù)監(jiān)測等，但是這些數(shù)值模擬及理論分析中常用的是簡諧二維波，SPECFEM2D程序包譜元法數(shù)值模擬在模擬應(yīng)力波傳播上的應(yīng)用，彌補(bǔ)了上述幾種數(shù)值模擬、現(xiàn)場或?qū)嶒?yàn)室試驗(yàn)以及理論分析的不足，譜元法的破裂震源更接近于現(xiàn)場實(shí)際情況，且SPECFEM2D可用于模擬煤礦巷道等小尺度范圍內(nèi)的礦震應(yīng)力波擴(kuò)散[13]。沖擊地壓巷道支護(hù)設(shè)計(jì)需要同時(shí)考慮動靜作用，以往的研究側(cè)重于巷道圍巖的靜載分布，對巷道周圍應(yīng)力波擴(kuò)散規(guī)律研究還不多見。本文采用SPECFEM2D程序包譜元法數(shù)值模擬方法，分析應(yīng)力波在巷道周圍擴(kuò)散過程中的波場快照、質(zhì)點(diǎn)最大振動速度PPV在巷道圍巖表面以及巷道圍巖分布規(guī)律以及應(yīng)力波記錄等。研究結(jié)果對于沖擊巷道危險(xiǎn)性分析及其合理支護(hù)設(shè)計(jì)具有借鑒意義。

1 譜元法及其計(jì)算流程

譜元法(Spectral-element method,SEM)最初是Patera 在1984年流體力學(xué)計(jì)算中提出[14]，Gharti等[15]在1998年把拉格朗日譜元法用于地震波模擬，推動了應(yīng)力波的數(shù)值模擬研究。譜元法的二維或三維彈性波波動方程強(qiáng)形式可以轉(zhuǎn)化為弱形式，使用基于高斯-勒讓德-洛巴托多項(xiàng)式(Gauss-Lobatto-Legendre,GLL)插值點(diǎn)[16]的拉格朗日高階正交多項(xiàng)式基函數(shù)獲得對角質(zhì)量矩陣[17]，大大簡化了算法實(shí)現(xiàn)和存儲量，并在時(shí)間域內(nèi)采用2階精度的顯式Newmark中心差分遞推格式求解，有利于計(jì)算機(jī)并行計(jì)算(Parallel computing,MPI)的實(shí)現(xiàn)[18]。譜元法能用較少的單元獲得與其他方法相同的精度，且具有最優(yōu)張量乘積技術(shù)的逐元算法，可處理復(fù)雜地質(zhì)結(jié)構(gòu)模型。與有限元法的區(qū)別在于基函數(shù)的選取和質(zhì)量矩陣中積分計(jì)算方法的不同，它兼具了有限元的處理邊界和結(jié)構(gòu)的靈活性和譜方法[19]的快速收斂特性。SPECFEM2D譜元法程序包可采用4節(jié)點(diǎn)或9節(jié)點(diǎn)網(wǎng)格單元對各向同性或各向異性介質(zhì)進(jìn)行彈性波場模擬、伴隨波場模擬以及全波形反演等，該程序包能模擬聲波和彈性波在局部、區(qū)域范圍內(nèi)各種介質(zhì)如流體、彈性體、粘彈性體、各向異性及多孔介質(zhì)中的擴(kuò)散過程。另外，SPECFEM2D/3D程序包具有與其他程序進(jìn)行耦合分析的擴(kuò)展性。

SPECFEM2D的非結(jié)構(gòu)化網(wǎng)格可由外部程序CUBIT，Gmsh，GID等軟件生成，由程序包自帶的xmeshfem2D命令也可以進(jìn)行簡單的網(wǎng)格劃分，并由SCOTH 進(jìn)行網(wǎng)格節(jié)點(diǎn)的負(fù)載分配，并通過轉(zhuǎn)換程序或程序包自帶的xmeshfem2D命令生成SPECFEM2D的網(wǎng)格文件、節(jié)點(diǎn)文件、材料屬性文件、應(yīng)力波吸收單元文件等程序運(yùn)行所需的參數(shù)輸入，然后程序初始化，得到離散方程組中的質(zhì)量矩陣、剛度矩陣、吸收邊界條件矩陣及震源的等效力項(xiàng)，然后用xspcfem2D命令進(jìn)行迭代求解，計(jì)算流程如圖1所示。

圖1 SPECFEM2D數(shù)值模擬流程Fig.1 Procedure of SPECFEM2Dnumerical simulation

2 譜元法數(shù)值模型建立

圖2 網(wǎng)格建立Fig.2 Grid establishment

圖3 SEM網(wǎng)格分塊Fig.3 SEM grid database partition

表1 材料屬性Table 1 Material properties

3 譜元法數(shù)值模擬結(jié)果分析

如圖4所示，淺色代表應(yīng)力波振幅為正值，深色代表應(yīng)力波振幅為負(fù)值，方框部分為應(yīng)力波監(jiān)測測點(diǎn)，為35 m×35 m的矩形范圍，每個(gè)測點(diǎn)相隔1 m，巷道內(nèi)部無測點(diǎn)，無巷道時(shí)布置1 332個(gè)測點(diǎn)，有巷道時(shí)布置1 310個(gè)測點(diǎn)，展示了應(yīng)力波在0.028,0.06,0.076和0.143 s擴(kuò)散時(shí)的波場快照。可以看出，CPML邊界應(yīng)力波吸收層達(dá)到了效果。在0.028 s時(shí)，震源開始擴(kuò)散；0.06 s時(shí)，P波和S波開始分離；0.076 s時(shí)，P波經(jīng)過巷道，由于巷道與周圍圍巖體的波阻抗不同，應(yīng)力波在巷道周圍擴(kuò)散，值得注意的是當(dāng)應(yīng)力波不在震源與巷道的直線方向上時(shí)，巷道的存在并不會引起應(yīng)力波擴(kuò)散的變化，因此，巷道周圍的應(yīng)力波擴(kuò)散樣式與震源類型和它與巷道的相對位置有關(guān)；0.143 s時(shí)，S波經(jīng)過巷道，可見在S波經(jīng)過巷道后，應(yīng)力波在巷道周圍擴(kuò)散的比在P波經(jīng)過巷道時(shí)劇烈，S波攜帶了大量的能量，是礦震應(yīng)力波導(dǎo)致巷道破壞的主要原因。

圖4 應(yīng)力波在有無巷道時(shí)擴(kuò)散過程波場快照Fig.4 Wave field snapshot of stress wave diffusion process with or without roadway

圖5 無巷道時(shí)巷道周圍PPV分布規(guī)律Fig.5 Distribution laws of PPV around roadway without roadway

圖6 有巷道時(shí)巷道周圍PPV分布規(guī)律Fig.6 Distribution laws of PPV around roadway

如圖7所示，在巷道圍巖表面布置23個(gè)測點(diǎn)，在周圍圍巖及內(nèi)部布置28個(gè)測點(diǎn)，箭頭從1開始依次表示了監(jiān)測點(diǎn)的序列號，每個(gè)測點(diǎn)相隔1 m。如圖8所示，巷道圍巖表面PPV最大值分布規(guī)律為：左幫>頂板>右?guī)?底板，分別為測點(diǎn)5、測點(diǎn)18、測點(diǎn)15和測點(diǎn)7。較無巷道時(shí)，PPV值增大倍數(shù)分別為2.4，1.7，1.7和1.4。如圖9所示，巷道圍巖內(nèi)部PPV在迎波側(cè)較無巷道時(shí)應(yīng)力波擴(kuò)散的PPV值大，背波側(cè)除了巷道圍巖表面的監(jiān)測點(diǎn)PPV值大于無巷道時(shí)的PPV值，其他均小于無巷道時(shí)的PPV值，可見，在背波側(cè)出現(xiàn)了1個(gè)應(yīng)力波擴(kuò)散的屏蔽區(qū)。巷道圍巖界面上入射壓縮波在自由表面反射成卸載波、再與入射卸載波相互作用，將形成拉應(yīng)力所產(chǎn)生的卸載效應(yīng)。在巷道圍巖高靜載接近煤巖體的極限強(qiáng)度時(shí)，極小的應(yīng)力波擾動引起的靜載增量都可能導(dǎo)致沖擊地壓的發(fā)生。

圖7 巷道圍巖測點(diǎn)布置示意Fig.7 Schematic diagram for layout of monitoring points in surrounding rock of roadway

圖8 巷道圍巖表面PPV分布規(guī)律Fig.8 Distribution law of PPV on surrounding rock surface of roadway

圖9 巷道圍巖內(nèi)部PPV分布規(guī)律Fig.9 Distribution laws of PPV in surrounding rock of roadway

由圖10所示，水平方向X和垂直方向Z的應(yīng)力波的S波幅值均遠(yuǎn)大于P波幅值。巷道左幫監(jiān)測點(diǎn)5表明頂?shù)装錢方向應(yīng)力波振幅放大明顯，尾波波動幅度比較大，應(yīng)力波持時(shí)顯著增大，Z方向應(yīng)力波與無巷道存在時(shí)相差較小。巷道底板監(jiān)測點(diǎn)11表明Z方向應(yīng)力波振幅放大明顯，尾波波動幅度比較大，應(yīng)力波持時(shí)顯著增大，X方向應(yīng)力波與無巷道存在時(shí)相差較小。受篇幅限制，在同時(shí)分析了巷道周圍其他監(jiān)測點(diǎn)應(yīng)力波幅值大小及衰減快慢后，發(fā)現(xiàn)應(yīng)力波在巷道周圍擴(kuò)散的速度矢量方向有一定的規(guī)律，也即巷道頂?shù)装逯饕芩椒较驊?yīng)力波作用，巷道兩幫主要受垂直方向應(yīng)力波作用，巷道兩幫徑向約束的最小主應(yīng)力較小，幅值較大的垂直方向應(yīng)力波與幫部垂直集中應(yīng)力疊加，其表現(xiàn)形式為巷道兩幫圍巖體向開挖空間拋射，出現(xiàn)層裂現(xiàn)象以及巷道肩窩處剪切滑移等，這與某礦現(xiàn)場出現(xiàn)的礦震應(yīng)力波破壞情況相吻合，如圖11所示。

圖10 巷道圍巖表面測點(diǎn)5和測點(diǎn)11 X及Z方向應(yīng)力波波形Fig.10 Stress wave waveforms in X and Z directions of No.5 and No.11 measuring points on surrounding rock surface of roadway

圖11 礦震應(yīng)力波致災(zāi)現(xiàn)場圖示Fig.11 Field diagram of disaster caused by mine earthquake stress wave

4 結(jié)論

1)質(zhì)點(diǎn)最大振動速度PPV>1 m/s在巷道迎波側(cè)的分布范圍較大，在背波側(cè)形成了一定范圍的應(yīng)力波屏蔽區(qū)。

2)礦震應(yīng)力波在遇到不同波阻抗的圍巖界面時(shí)，產(chǎn)生應(yīng)力波的疊加現(xiàn)象，巷道圍巖表面PPV最大值發(fā)生左幫上半部分，是無巷道影響時(shí)的2.4倍，深入巷道左幫、底板圍巖內(nèi)部的測點(diǎn)有巷道時(shí)PPV值反而沒有無巷道時(shí)大。

3)巷道兩幫垂直方向礦震應(yīng)力波幅值較水平方向大，頂?shù)装逅椒较驊?yīng)力波幅值較垂直方向大。且幅值較大的垂直或水平方向應(yīng)力波尾波幅值較大，持續(xù)時(shí)間長。

4)PPV值在巷道圍巖內(nèi)形成了升高區(qū)和屏蔽區(qū)，在巷道圍巖集中應(yīng)力區(qū)，較小的應(yīng)力波擾動極可能誘發(fā)沖擊顯現(xiàn)，數(shù)值模擬有助于加深理解巷道沖擊地壓的礦震應(yīng)力波誘發(fā)機(jī)理，并為下一步合理的沖擊巷道支護(hù)設(shè)計(jì)提供了參考。