巷道与钻孔交叉三维数值建模方法＊

张天军，张 磊，，张 超，包若羽，尚宏波

(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.