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

張天軍，張 磊，，張 超，包若羽，尚宏波

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

0 引 言

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

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

1 平行鉆孔模型

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

1.1 生成塊狀單元

以巷道及鉆孔的半徑，長度為參數，生成相應8 個基本塊狀單元，如圖2 所示。模型上下分2層，其中下層深藍色、黃色(位于立方體背面)、紅色及綠色分別為巷道鉆孔交叉區域，巷道區域，鉆孔區域和無鉆孔煤層區域，記為組1 ～4. 上層青色、紫色、咖啡色與橙色分別為巷道鉆孔交叉頂部，巷道頂部，鉆孔頂部，無鉆孔煤層區域頂部，記為組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 鉆孔區域網格變形

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

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

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

其中 Z(NI，1，k)為節點NI，1，k的Z 坐標。

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

隨后對組3 的右側面內所有的網格做變形，即調整NI，j，k至NI，J，K共計(J -1)×(K -1)個節點的Y 坐標及Z 坐標，使網格均勻變化。節點NI，j，k的Y坐標以及Z 坐標如式(11)及式(12)。

之后對組3 的左側面網格做變形，即調整N1，1，1至N1，J，K共計J×K 個節點的坐標，使左面形成與巷道相切的弧面。節點N1，j，k的坐標如式(13)。

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

1.3 交叉區域網格變形

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

式中 R 為巷道半徑。

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

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

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

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

最后對組1 內部的所有網格變形，即調整內部節點N2，1，1至NI2，J-1，K-1，共計(I2 -1)×(J -1)×(K-1)個，節點Ni，J，k坐標如式21，調整完成如圖3(d)。

1.4 巷道區域網格變形

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

1.5 無鉆孔煤層區域網格變形

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

1.6 頂部區域網格變形

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

1.7 鏡像命令補全模型

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

2 復雜鉆孔模型

復雜模型主要指傾斜鉆孔模型及多孔模型。

2.1 傾斜鉆孔模型

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

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

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

2.2 多孔模型

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

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

圖3 網格變形過程Fig.3 Process of shaping grid

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

3 應用實例

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

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

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

圖6 鉆孔內應力狀態Fig.6 Stress states of borehole

(a)X 應力分布 (b)Y 應力分布(c)Z 應力分布 (d)塑性區域分布圖

4 結 論

針對FLAC3D建模復雜的問題，提出基于FISH語言的巷道鉆孔交叉網格變形建模法，編制了用于采礦工程中的巷道鉆孔交叉結構的模型建立程序。該方法及理論直白易懂，操作簡單易行;該程序實現了建模的自動化、參數化以及可重用化。減小了建模的時間與難度，提高建模的精度與效益。通過應用實例驗證了所述方法的可行性與有效性，有助于幫助研究人員深入了解復雜建模的本質，解決了巷道鉆孔交叉建模的難題，在采礦工程領域FLAC3D模型建立方面具有借鑒意義。

References

［1］ 陳育民，劉漢龍. 鄧肯-張本構模型在FLAC3D中的開發與實現［J］. 巖土力學，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］ 王松陽，李樹剛，王紅勝.水滲流對薄層復合頂板的弱化影響［J］.西安科技大學學報，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.，Fast Language Analysis of continua in three dimensions versin 3.0 user’s mannual［R］.USA:Itasca Consulting Group，Inc.，2005.

［4］ 廖秋林，曾錢幫，劉 彤，等. 基于ANSYS 平臺復雜地質體FLAC3D模型的自動生成［J］. 巖石力學與工程學報，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］ 李 翔，王金安，張少杰.復雜地質體三維數值建模方法研究［J］.西安科技大學學報，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］.西安:西安科技大學，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］ 馬長年，徐國元，江文武，等. 復雜開挖過程FLAC3D力學仿真代碼生成系統研究［J］.巖土力學，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前處理程序的開發及仿真效果檢驗［J］.巖石力學與工程學報，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］ 張 超，林柏泉，周 延，等.本煤層近水平瓦斯抽采鉆孔“強弱強”帶壓封孔技術研究［J］.采礦與安全工程學報，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］張 超，林柏泉，周 延，等.本煤層深孔定向靜態破碎卸壓增透技術研究與應用［J］.采礦與安全工程學報，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］王宏圖，江記記，王再清，等.本煤層單一順層瓦斯抽采鉆孔的滲流場數值模擬［J］.重慶大學學報，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 的采動巖體等效連續介質流固耦合模型及應用［J］.采礦與安全工程學報，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］王兆豐，周少華，李志強.瓦斯抽采鉆孔有效抽采半徑的數值計算方法［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］.采礦與安全工程學報，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］高富強.斷面形狀對巷道圍巖穩定性影響的數值模擬分析［J］.山東科技大學學報:自然科學版，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］張強勇，張緒濤，向 文，等.不同洞形與加載方式對深部巖體分區破裂影響的模型試驗研究［J］.巖石力學與工程學報，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］王 振，梁運培，金洪偉.防突鉆孔失穩的力學條件分析［J］.采礦與安全工程學報，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］姚向榮，程功林，石必明.深部圍巖遇弱結構瓦斯抽采鉆孔失穩分析與成孔方法［J］.煤炭學報，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.