巷道與鉆孔交叉三維數(shù)值建模方法＊

張?zhí)燔?，?磊，，張 超，包若羽，尚宏波

(1.西安科技大學(xué) 理學(xué)院，陜西 西安710054;2.西安科技大學(xué) 能源學(xué)院，陜西 西安710054)

0 引 言

FLAC3D是由美國ITASCA 國際咨詢與軟件開發(fā)公司開發(fā)的一種采用拉格朗日有限差分法求解的連續(xù)介質(zhì)力學(xué)分析軟件。由于其對地質(zhì)材料的模擬程度高，能夠分析材料在強(qiáng)度極限以及屈服極限時產(chǎn)生的破壞和失穩(wěn)。尤其在模擬大變形時，可以較為真實地還原所發(fā)生的破壞和塑形流動。因此在巖土工程界具有深遠(yuǎn)的影響，在我國被廣泛地應(yīng)用于建筑、交通、采礦、地質(zhì)等行業(yè)［1-3］。FLAC3D 軟件提供了13 種基本形狀網(wǎng)格，具有良好的FISH 語言接口，但其復(fù)雜的建模方式一直困擾著研究人員。許多學(xué)者通過引入ANSYS［4］，SURPAC［5］，AutoCAD［6］等三維建 模軟件來簡化建模工作。在此種情況下，研究人員通過自主開發(fā)仿真代碼生成系統(tǒng)［7］，前處理程序［8］來提高建模效率。相關(guān)研究的核心都集中在如何建立滿足需求的模型，而在瓦斯抽采鉆孔參數(shù)選型設(shè)計時，需要生成大量模型，以上方法因借助輔助軟件生成模型，需耗費大量時間建立三維模型。網(wǎng)格變形建模法是基于FISH 語言編寫的建模程序，具有易修改巷道、鉆孔等尺寸，易建立多組模型等特點。為此深入了解FLAC3D的網(wǎng)格變形建模思想極為重要，是進(jìn)行高精度計算、二次開發(fā)、自主研發(fā)以及批量建模的必備條件。

瓦斯抽采是一種廣泛應(yīng)用于煤礦生產(chǎn)當(dāng)中的重要技術(shù)［9-10］，其核心技術(shù)在于維持鉆孔穩(wěn)定，有效地封堵鉆孔。通過數(shù)值計算研究瓦斯抽采鉆孔穩(wěn)定性是一種高效、經(jīng)濟(jì)、安全的手段［11-14］。高富強(qiáng)［15］通過數(shù)值模擬研究了圓形、拱形以及矩形巷道的穩(wěn)定性，得出圓形巷道穩(wěn)定性最高的結(jié)論。張強(qiáng)勇［16］研究表明，深部巷道圍巖的破裂形狀與洞型無關(guān)，不同斷面形狀的應(yīng)力分布趨勢相近。由于圓形巷道存在幾何對稱特性，采用對稱建模法，可快速的建立計算模型。在深部開采的條件下，圓形巷道能夠在滿足計算精度要求的情況下，較為便捷的建立計算模型。文中通過FISH 語言編寫程序，建立圓形巷道與鉆孔交叉模型，探究基于參數(shù)的網(wǎng)格建模方法，提出了網(wǎng)格變形建模法，該方法能夠完成巷道與鉆孔交叉條件下的高效建模。

1 平行鉆孔模型

文中以本煤層瓦斯抽采鉆孔為建模對象，抽采鉆孔等間距地分布布置于巷道一側(cè)。鉆孔詳細(xì)布置方式如圖1(a)，1(b)表示，簡化后利用網(wǎng)格變形建模法建立的FLAC3D計算模型如圖1(c)所示。平行鉆孔的網(wǎng)格變形建模法步驟如下。

1.1 生成塊狀單元

以巷道及鉆孔的半徑，長度為參數(shù)，生成相應(yīng)8 個基本塊狀單元，如圖2 所示。模型上下分2層，其中下層深藍(lán)色、黃色(位于立方體背面)、紅色及綠色分別為巷道鉆孔交叉區(qū)域，巷道區(qū)域，鉆孔區(qū)域和無鉆孔煤層區(qū)域，記為組1 ～4. 上層青色、紫色、咖啡色與橙色分別為巷道鉆孔交叉頂部，巷道頂部，鉆孔頂部，無鉆孔煤層區(qū)域頂部，記為組5 ～8.組1，2，5，6 長度等于巷道的半徑R，組3，4，7，8 長度等于鉆孔的長度l，組1，3，5，7 的厚度等于鉆孔的半徑r.

圖1 巷道鉆孔交叉模型Fig.1 Model of roadway borehole intersecting

圖2 基本塊狀單元Fig.2 Brick units

1.2 鉆孔區(qū)域網(wǎng)格變形

首先對鉆孔區(qū)域，即組3 進(jìn)行網(wǎng)格變形，假設(shè)組3 含有網(wǎng)格為I×J×K 個，即沿X 方向有I 個節(jié)點，沿Y 方向有J 個節(jié)點，沿Z 方向有K 個節(jié)點，每個節(jié)點用Ni，j，k表示，其中1≤i≤I，1≤j≤J，1≤k≤K，i，j，k 均為整數(shù)。沿組3 的右側(cè)面底部棱進(jìn)行變形，將組3 網(wǎng)格Y 方向尺寸更改為鉆孔半徑，即修改NI，1，K至NI，J，K共計J 個節(jié)點的Y 坐標(biāo)，變形前如圖3(a)所示。相鄰節(jié)點的間距dy的計算公式如式(1)，節(jié)點NI，j，K的Y 坐標(biāo)計算如式(2)。

其中 r 為 鉆 孔 半 徑;Y(NI，j，K)為 節(jié) 點NI，j，K的Y坐標(biāo)。

其次沿組3 的右側(cè)面前部棱進(jìn)行變形，將組3網(wǎng)格的Z 方向尺寸更改為鉆孔半徑，即修改NI11至NI1K共計K 個節(jié)點的Z 坐標(biāo)。棱上相鄰節(jié)點的間距dz的計算公式為式(3)，節(jié)點Ni的Z 坐標(biāo)計算如式(4)。

其中 Z(NI，1，k)為節(jié)點NI，1，k的Z 坐標(biāo)。

接著分別沿組3 的右側(cè)面頂部棱及后部棱進(jìn)行變形，通過修改節(jié)點NI，1，K至NI，J，K以及NI，J，1至NI，J，K，共計J +K 個節(jié)點的Y 坐標(biāo)以及Z 坐標(biāo)，將頂部及后部近似為半徑r 的1/4 圓周。節(jié)點NI，j，K的Y 坐標(biāo)以及Z 坐標(biāo)如式(5)及式(6)，節(jié)點NI，J，k坐標(biāo)如式(8)及式(9)。

隨后對組3 的右側(cè)面內(nèi)所有的網(wǎng)格做變形，即調(diào)整NI，j，k至NI，J，K共計(J -1)×(K -1)個節(jié)點的Y 坐標(biāo)及Z 坐標(biāo)，使網(wǎng)格均勻變化。節(jié)點NI，j，k的Y坐標(biāo)以及Z 坐標(biāo)如式(11)及式(12)。

之后對組3 的左側(cè)面網(wǎng)格做變形，即調(diào)整N1，1，1至N1，J，K共計J×K 個節(jié)點的坐標(biāo)，使左面形成與巷道相切的弧面。節(jié)點N1，j，k的坐標(biāo)如式(13)。

最后對組3 內(nèi)部的所有網(wǎng)格做變形，修改其余節(jié)點的坐標(biāo)，即調(diào)整N2，1，1至NI-1，J，K共計(I -2)×J×K 個節(jié)點的坐標(biāo)，使組3 的整體形狀成為1/4圓柱，如圖3(b)，節(jié)點Ni，j，k的坐標(biāo)如式(14)。

1.3 交叉區(qū)域網(wǎng)格變形

完成鉆孔區(qū)域后，對巷道鉆孔交叉區(qū)域，即組1 進(jìn)行網(wǎng)格變形，假設(shè)組1 含有網(wǎng)格為I2 ×J ×K個，則組1 中的節(jié)點可用Ni，j，k表示，其中1≤i≤I2，1≤j≤J，1≤k≤K，i，j，k 均為整數(shù)。首先對組1 頂面做變形，即調(diào)整N1，1，K至NI2，J，K，共計I2 ×J 個節(jié)點的坐標(biāo)，變形前如圖3(c)，節(jié)點Ni，j，k的調(diào)整后坐標(biāo)如式(15)。

式中 R 為巷道半徑。

其次 對 組1 前 面 做 變 形，即 調(diào) 整N1，1，1至NI2，1，K，共計I2 ×K 個節(jié)點的坐標(biāo)，節(jié)點Ni，1，k的坐標(biāo)如式(17)。

接著 對 組1 底 面 做 變 形，即 調(diào) 整N1，1，1至NI2，J，1，共計I2 ×J 個節(jié)點的坐標(biāo)，節(jié)點Ni，j，1的坐標(biāo)如式(18)。

隨后 對 組1 左 面 做 變 形，即 調(diào) 整N1，1，1至N1，J，K，共計J ×K 個節(jié)點的坐標(biāo)，節(jié)點N1，j，k坐標(biāo)如式(19)。

然后 對 組1 后 面 做 變 形，即 調(diào) 整N1，J，1至NI2，J，K，共計I2 ×K 節(jié)點的坐標(biāo)，節(jié)點Ni，J，k坐標(biāo)如式(20)。

最后對組1 內(nèi)部的所有網(wǎng)格變形，即調(diào)整內(nèi)部節(jié)點N2，1，1至NI2，J-1，K-1，共計(I2 -1)×(J -1)×(K-1)個，節(jié)點Ni，J，k坐標(biāo)如式21，調(diào)整完成如圖3(d)。

1.4 巷道區(qū)域網(wǎng)格變形

巷道區(qū)域即為組2 區(qū)域，假設(shè)組2 內(nèi)的網(wǎng)格數(shù)量為I2 ×J2 ×K，如圖3(e)。用Ni，j，k表示其中的節(jié)點，其中1≤i≤I2，1≤j≤J2，1≤k≤K，i，j，k 均為整數(shù)。對巷道區(qū)域網(wǎng)格做變形，僅需按照巷道交叉區(qū)域的后面修正巷道區(qū)域內(nèi)的節(jié)點，將整個區(qū)域變形為1/4 圓柱狀，變形后如圖3(f)，節(jié)點Ni，j，k坐標(biāo)如式(22)。

1.5 無鉆孔煤層區(qū)域網(wǎng)格變形

無鉆孔煤層區(qū)域即為組4 內(nèi)區(qū)域，假設(shè)組4 內(nèi)的網(wǎng)格數(shù)量為I×J2 ×K，用Ni，j，k表示其中的節(jié)點，其中1≤i≤I，1≤j≤J2，1≤k≤K，i，j，k 均為整數(shù)。變形前如圖3(g)，對此區(qū)域變形需要參照巷道區(qū)域以及鉆孔區(qū)域修改組4 內(nèi)部節(jié)點，修正邊界存在的不規(guī)則形狀，變形后如圖3(h)，節(jié)點坐標(biāo)如式(23)。

1.6 頂部區(qū)域網(wǎng)格變形

完成巷道與鉆孔的整體變形后，頂部區(qū)域仍然存在網(wǎng)格不規(guī)則的情況，如圖3(i)。按照巷道鉆孔交叉區(qū)域，巷道區(qū)域，鉆孔區(qū)域及無鉆孔煤層區(qū)域的頂部修正頂部4 個組，如圖3(j)。節(jié)點Ni，j，k坐標(biāo)如式(24)。

1.7 鏡像命令補全模型

經(jīng)過以上6 步，建立了1/4 巷道鉆孔交叉模型，使用FLAC3D中的鏡像命令，即可完成整個模型的建立。

2 復(fù)雜鉆孔模型

復(fù)雜模型主要指傾斜鉆孔模型及多孔模型。

2.1 傾斜鉆孔模型

傾斜鉆孔的網(wǎng)格變形建模首先根據(jù)需要建立相應(yīng)的平行鉆孔，如圖4(a)，接著對鉆孔部分進(jìn)行網(wǎng)格變形。

經(jīng)過平行鉆孔7 步后，對鉆孔區(qū)域，無鉆孔煤層區(qū)域及兩其頂部區(qū)域，即組3，組4，組7，組8 沿Y 軸及Z 軸進(jìn)行平移，得到所需鉆孔角度，如圖4(b)。節(jié)點Ni，j，k坐標(biāo)如式(25)。

式中 ρ1為鉆孔在XOY 平面內(nèi)投影與X 軸的夾角;ρ2為鉆孔與Z 軸夾角。

2.2 多孔模型

建立單孔模型后，對平行鉆孔可以直接使用鏡像命令復(fù)制，形成多孔模型。而對于傾斜鉆孔，由于其模型外形不規(guī)則，不能使用前述方法。采用導(dǎo)出移動導(dǎo)入法，能夠較好的解決此問題。首先使用export 命令導(dǎo)出模型，其次對模型進(jìn)行全坐標(biāo)移動，移動后的坐標(biāo)如式(26)。再次使用import命令導(dǎo)入模型。最終得到多孔模型如圖4(c)。

式中 P 為模型沿Y 軸移動量，大小等于模型的Y方向尺寸。

圖3 網(wǎng)格變形過程Fig.3 Process of shaping grid

圖4 鉆孔變形與移動Fig.4 Shaping and moving borehole

3 應(yīng)用實例

建模對象為近水平煤層，瓦斯抽采鉆孔使用斜向孔，孔長為100 m，鉆孔與X 軸夾角為30°，鉆孔間距為1.0 m，布孔方式如圖1(b)所示，計算模型如圖1(c)所示。

根據(jù)圖5 中位移曲線能夠看出，鉆孔受到了巷道的影響，在0 ～4 m 的范圍內(nèi)產(chǎn)生較大的位移，4～20 m 的范圍內(nèi)，位移逐漸減小，而在20 m 外的鉆孔位移均保持一致。從圖6 中的應(yīng)力分布可以明顯地看出卸壓區(qū)、峰后應(yīng)力集中區(qū)、峰前應(yīng)力集中區(qū)、原巖應(yīng)力區(qū)3 個不同力學(xué)狀態(tài)的區(qū)域，在孔口處出現(xiàn)大范圍的塑性區(qū)域，煤巖體破裂。1 ～2倍巷道半徑處，出現(xiàn)較明顯的應(yīng)力集中現(xiàn)象，這與當(dāng)前的理論研究［17-18］是相吻合的，模型能夠較好的反應(yīng)巷道鉆孔交叉情況下的力學(xué)狀態(tài)。

圖5 鉆孔壁位移Fig.5 Displacement of borehole wall

圖6 鉆孔內(nèi)應(yīng)力狀態(tài)Fig.6 Stress states of borehole

(a)X 應(yīng)力分布 (b)Y 應(yīng)力分布(c)Z 應(yīng)力分布 (d)塑性區(qū)域分布圖

4 結(jié) 論

針對FLAC3D建模復(fù)雜的問題，提出基于FISH語言的巷道鉆孔交叉網(wǎng)格變形建模法，編制了用于采礦工程中的巷道鉆孔交叉結(jié)構(gòu)的模型建立程序。該方法及理論直白易懂，操作簡單易行;該程序?qū)崿F(xiàn)了建模的自動化、參數(shù)化以及可重用化。減小了建模的時間與難度，提高建模的精度與效益。通過應(yīng)用實例驗證了所述方法的可行性與有效性，有助于幫助研究人員深入了解復(fù)雜建模的本質(zhì)，解決了巷道鉆孔交叉建模的難題，在采礦工程領(lǐng)域FLAC3D模型建立方面具有借鑒意義。

References

［1］ 陳育民，劉漢龍. 鄧肯-張本構(gòu)模型在FLAC3D中的開發(fā)與實現(xiàn)［J］. 巖土力學(xué)，2007，28(10):2 123 -2 126.CHEN Yu-min，LIU Han-long. Development and implementation of Duncan-Chang constitutive model in FLAC3D［J］.Rock and Soil Mechanics，2007，28(10):2 123 -2 126.

［2］ 王松陽，李樹剛，王紅勝.水滲流對薄層復(fù)合頂板的弱化影響［J］.西安科技大學(xué)學(xué)報，2012，32(6):686-690.WANG Song-yang，LI Shu-gang，WANG Hong-sheng.Weaken influence of water seepage on thin-layer compound roof［J］. Journal of Xi’an University of Science and Technology，2012，32(6):686 -690，717.

［3］ Itasca Consulting Group，Inc.，F(xiàn)ast Language Analysis of continua in three dimensions versin 3.0 user’s mannual［R］.USA:Itasca Consulting Group，Inc.，2005.

［4］ 廖秋林，曾錢幫，劉 彤，等. 基于ANSYS 平臺復(fù)雜地質(zhì)體FLAC3D模型的自動生成［J］. 巖石力學(xué)與工程學(xué)報，2005，24(6):1 010 -1 013.LIAO Qiu-lin，ZENG Qian-bang，LIU Tong，et al. Automatic model generation of complex geologicbody with FLAC3Dbased on ANSYS platform［J］. Chinese Journal of Rock Mechanics and Engineering，2005，24(6):1 010 -1 013.

［5］ 李 翔，王金安，張少杰.復(fù)雜地質(zhì)體三維數(shù)值建模方法研究［J］.西安科技大學(xué)學(xué)報，2012，32(6):676-681.LI Xiang，WANG Jin-an，ZHANG Shao-jie.3D modeling method of complicated geological body［J］. Journal of Xi’an University of Science and Technology，2012，32(6):676 -681.

［6］ 高永剛.基于AutoCAD 的FLAC3D模型快速建模方法研究［D］.西安:西安科技大學(xué)，2012.GAO Yong-gang. Study on rapid modeling method of FLAC3Dmodel based on AutoCAD file［D］.Xi’an:Xi’an University of Science and Technology，2012.

［7］ 馬長年，徐國元，江文武，等. 復(fù)雜開挖過程FLAC3D力學(xué)仿真代碼生成系統(tǒng)研究［J］.巖土力學(xué)，2012，33(8):2 536 -2 542.MA Chang-nian，XU Guo-yuan，JIANG Wen-wu，et al.Study of generating code system of FLAC3Dfor simulating complicated excavating process［J］. Rock and Soil Mechanics，2012，33(8):2 536 -2 542.

［8］ 胡 斌，張倬元，黃潤秋，等.FLAC3D前處理程序的開發(fā)及仿真效果檢驗［J］.巖石力學(xué)與工程學(xué)報，2002，21(9):1 387 -1 391.HU Bin，ZHANG Zhuo-yuan，HUANG Run-qiu，et al.Development of pre-processing package for FLAC3Dand verification of its simulating effects［J］.Chinese Journal of Rock Mechanics and Engineering，2002，21(9):1 387 -1 391.

［9］ 張 超，林柏泉，周 延，等.本煤層近水平瓦斯抽采鉆孔“強(qiáng)弱強(qiáng)”帶壓封孔技術(shù)研究［J］.采礦與安全工程學(xué)報，2013，30(6):935 -939.ZHANG Chao，LIN Bai-quan，ZHOU Yan，et al. Strongweak-strong borehole pressurized sealing technology for horizontal gas drainage borehole in mining seam［J］.Journal of Mining ＆ Safety Engineering，2013，30(6):935 -939.

［10］張 超，林柏泉，周 延，等.本煤層深孔定向靜態(tài)破碎卸壓增透技術(shù)研究與應(yīng)用［J］.采礦與安全工程學(xué)報，2013，30(4):600 -604.ZHANG Chao，LIN Bai-quan，ZHOU Yan，et al. Deephole directional static cracking technique for pressure relief and permeability improvement in mining-coal bed［J］. Journal of Mining ＆ Safety Engineering，2013，30(4):600 -604.

［11］王宏圖，江記記，王再清，等.本煤層單一順層瓦斯抽采鉆孔的滲流場數(shù)值模擬［J］.重慶大學(xué)學(xué)報，2011，34(4):24 -29.WANG Hong-tu，JIANG Ji-ji，WANG Zai-qing，et al.Numerical simulation of seepage field of gas extraction drilling of single bedding of mining-coal bed［J］.Journal of Chongqing University，2011，34(4):24 -29.

［12］白國良.基于FLAC3D 的采動巖體等效連續(xù)介質(zhì)流固耦合模型及應(yīng)用［J］.采礦與安全工程學(xué)報，2010，27(1):106 -110.BAI Guo-liang.Fluid-solid coupling model of equivalent continuous medium in overburden rock based on FLAC3Dand its application［J］. Journal of Mining and Safety Engineering，2010，27(1):106 -110.

［13］王兆豐，周少華，李志強(qiáng).瓦斯抽采鉆孔有效抽采半徑的數(shù)值計算方法［J］.煤炭工程，2011(6):82 -84.WANG Zhao-feng，ZHOU Shao-hua，LI Zhi-qiang. Numerical calculation method of effective drainage radius for gas drainage borehole［J］. Coal Engineering，2011(6):82 -84.

［14］魯 義，申宏敏，秦波濤，等.順層鉆孔瓦斯抽采半徑及布孔間距研究［J］.采礦與安全工程學(xué)報，2015，32(1):156 -162.LU Yi，SHEN Hong-min，QIN Bo-tao，et al.Gas drainage radius and borehole distance along seam［J］. Journal of Mining and Safety Engineering，2015，32(1):156 -162.

［15］高富強(qiáng).斷面形狀對巷道圍巖穩(wěn)定性影響的數(shù)值模擬分析［J］.山東科技大學(xué)學(xué)報:自然科學(xué)版，2007，26(2):43 -46.GAO Fu-qiang.The effect of numerical simulation analysis of cross-section shapes of roadways on surrounding rock stability［J］.Journal of Shandong University of Science and Technology:Natural Science Edition，2007，26(2):43 -46.

［16］張強(qiáng)勇，張緒濤，向 文，等.不同洞形與加載方式對深部巖體分區(qū)破裂影響的模型試驗研究［J］.巖石力學(xué)與工程學(xué)報，2013，32(8):1 564 -1 571.ZHANG Qiang-yong，ZHANG Xu-tao，XIANG Wen，et al.Model test study of zonal disintegration in deep rock mass under different cavern shapes and loading conditions［J］.Chinese Journal of Rock Mechanics and Engineering，2013，32(8):1 564 -1 571.

［17］王 振，梁運培，金洪偉.防突鉆孔失穩(wěn)的力學(xué)條件分析［J］.采礦與安全工程學(xué)報，2008，25(4):444 -448.WANG Zhen，LIANG Yun-pei，JIN Hong-wei. Analysis of mechanics conditions for instability of outburst-preventing borehole［J］. Journal of Mining ＆ Safety Engineering，2008，25(4):444 -448.

［18］姚向榮，程功林，石必明.深部圍巖遇弱結(jié)構(gòu)瓦斯抽采鉆孔失穩(wěn)分析與成孔方法［J］.煤炭學(xué)報，2010，35(12):2 073 -2 081.YAO Xiang-rong，CHENG Gong-lin，SHI Bi-ming.Analysis on gas extraction drilling instability and control method of pore-forming in deep surrounding-rock with weak structure［J］.Journal of China Coal Society，2010，35(12):2 073 -2 081.