基金项目:国家自然科学基金(51174158,51374168)
通讯作者:张天军(1971-),男,陕西临潼人,教授,E-mail:tianjun_zhang@126.com
1.西安科技大学 理学院,陕西 西安 710054; 2.西安科技大学 能源学院,陕西 西安 710054
(1.College of Sciences,Xi'an University of Science and Technology,Xi'an 710054,China; 2.College of Energy Science and Engineering,Xi'an University of Science and Technology,Xi'an 710054,China)
FLAC3D; roadway borehole intersecting; gas drainage borehole; Numerical Simulation; borehole modeling
DOI: 10.13800/j.cnki.xakjdxxb.2015.0603
针对FLAC3D建模难度高,效率低的问题,提出一种具有复用性、参数化特点的网格变形建模法,且依照该方法编写了基于FISH语言建模程序。该程序首先用块状单元生成模型网格,接着根据巷道钻孔交叉的几何关系计算各节点坐标,通过程序依次调整模型各部分网格形状,从而准确地建立巷道钻孔交叉计算模型。通过本煤层瓦斯抽采钻孔的计算实例,研究含有巷道影响下的钻孔应力。计算结果能够较为准确地反应钻孔的应力状态,验证了所述方法的可行性与准确性。该方法原理简单,能够快速地建立不同参数的大量模型,有助于理解FLAC3D巷道钻孔交叉几何建模原理,为深入研究复杂建模机制提供参考,降低研究巷道钻孔交叉问题的难度,提高抽采钻孔问题研究效率。
A mesh morphing method is proposed and a program based on FISH is written,which has parameterization and reusable features,as there are difficulty and inefficient in FLAC3D modeling.The grid is generated by brick unit in this program.Then the grid position is calculated according to the geometrical relationship of roadway borehole intersecting.And the grid is adjusted by the program.In this way,an accurate model of roadway borehole intersecting is obtained.The stress of borehole under the influence of roadway is studied by a calculation case of horizontal gas drainage borehole in mining seam.The stress state of borehole is accurately reflected according to the result.Feasibility and accuracy are verified.This method is easy to understand and can quickly build a large number of models with different parameters.The difficulty of study on the roadway borehole intersecting is reduced.This method contributes to understanding of the principle of modeling in FLAC3D and provides researchers with reference on depth of complex modeling mechanism.
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(a),1(b)表示,简化后利用网格变形建模法建立的FLAC3D计算模型如图1(c)所示。平行钻孔的网格变形建模法步骤如下。
以巷道及钻孔的半径,长度为参数,生成相应8个基本块状单元,如图2所示。模型上下分2层,其中下层深蓝色、黄色(位于立方体背面)、红色及绿色分别为巷道钻孔交叉区域,巷道区域,钻孔区域和无钻孔煤层区域,记为组1~4.上层青色、紫色、咖啡色与橙色分别为巷道钻孔交叉顶部,巷道顶部,钻孔顶部,无钻孔煤层区域顶部,记为组5~8.组1,2,5,6长度等于巷道的半径R,组3,4,7,8长度等于钻孔的长度l,组1,3,5,7的厚度等于钻孔的半径r.
首先对钻孔区域,即组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)。
dy=r/J,(1)
Y(NI,j,K)=dy·j.(2)
其中 r为钻孔半径; Y(NI,j,K)为节点NI,j,K的Y坐标。
其次沿组3的右侧面前部棱进行变形,将组3网格的Z方向尺寸更改为钻孔半径,即修改NI11至NI1K共计K个节点的Z坐标。棱上相邻节点的间距dz的计算公式为式(3),节点Ni的Z坐标计算如式(4)。
dz=r/K,(3)
Z(NI,j,K)=dz·k.(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)。
Y(NI,j,K)=cosθj·r,(5)
Z(NI,j,K)=sinθJ·r,(6)
θj=(π)/2-(π)/4·j/J,(7)
Y(NI,J,K)=sinθk·r,(8)
Z(NI,J,K)=cosθk·r,(9)
θk=(π)/4·k/K,(10)
随后对组3的右侧面内所有的网格做变形,即调整NI,j,k至NI,J,K共计(J-1)×(K-1)个节点的Y坐标及Z坐标,使网格均匀变化。节点NI,j,k的Y坐标以及Z坐标如式(11)及式(12)。
Y(NI,j,k)=Y(NI,1,k)+j/J·(Y(NI,J,K)-Y(NI,1,k)),(11)
Z(NI,j,k)=Z(NI,j,1)+k/K·(Z(NI,J,K)-Z(NI,j,1)),(12)
之后对组3的左侧面网格做变形,即调整N1,1,1至N1,J,K共计J×K个节点的坐标,使左面形成与巷道相切的弧面。节点N1,j,k的坐标如式(13)。
{X(N1,j,k)=(R2-Z(N1,j,k)2)1/2
Y(N1,j,k)=Y(NI,j,k)
Z(N1,j,k)=Z(NI,j,k)(13)
最后对组3内部的所有网格做变形,修改其余节点的坐标,即调整N2,1,1至NI-1,J,K共计(I-2)×J×K个节点的坐标,使组3的整体形状成为1/4圆柱,如图3(b),节点Ni,j,k的坐标如式(14)。
{X(Ni,j,k=X(N1,j,k)+i/I·(X(NI,j,k-X(N1,j,k))
Y(N1,j,k)=Y(NI,j,k)
Z(N1,j,k)=Z(NI,j,k).(14)
完成钻孔区域后,对巷道钻孔交叉区域,即组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)。
{X(Ni,j,K)=R·sinθi,j
Y(Ni,j,K)=j/JY(NI2,J,K)
Z(Ni,j,K)=R·cosθi,j,(15)
θi,j=(π)/2-arctan((Z(Ni,j,K))/(X(Ni,j,K))),(16)
式中 R为巷道半径。
其次对组1前面做变形,即调整N1,1,1至NI2,1,K,共计I2×K个节点的坐标,节点Ni,1,k的坐标如式(17)。
{X(Ni,1,k)=i/I·X(NI2,1,k)
Z(Ni,1,k)=k/K·Z(Ni,1,K),(17)
接着对组1底面做变形,即调整N1,1,1至NI2,J,1,共计I2×J个节点的坐标,节点Ni,j,1的坐标如式(18)。
{X(Ni,j,1)=i/I·X(NI2,j,1)
Y(Ni,j,1)=k/K·Y(NI2,j,1),(18)
随后对组1左面做变形,即调整N1,1,1至N1,J,K,共计J×K个节点的坐标,节点N1,j,k坐标如式(19)。
{Y(N1,j,k)=i/I·Y(N1,J,k)
Z(N1,j,k)=k/K·Z(N1,j,K),(19)
然后对组1后面做变形,即调整N1,J,1至NI2,J,K,共计I2×K节点的坐标,节点Ni,J,k坐标如式(20)。
{X(Ni,J,k)=i/I·X(NI2,j,1)
Y(Ni,J,k)=Y(N1,J,1)+i/I·(Y(NI2,J,k-Y(N1,J,1))
Y(Ni,J,k)=k/K+Ri/I·(Z(NI2,J,k-k/K·R),(20)
最后对组1内部的所有网格变形,即调整内部节点N2,1,1至NI2,J-1,K-1,共计(I2-1)×(J-1)×(K-1)个,节点Ni,J,k坐标如式21,调整完成如图3(d)。
{X(Ni,J,k)=i/I·X(NI2,j,1)
Y(Ni,J,k)=Y(N1,J,1)+i/I·(Y(NI2,J,k-Y(N1,J,1))
Y(Ni,J,k)=Z(N1,J,1)+i/I·(Z(NI2,J,k)-Z(N1,J,1)).(21)
巷道区域即为组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)。
{X(Ni,j,k)=X(Ni,1,k)
Y(Ni,j,k)=Y(Ni,1,k)+j/J((R+L)-Y(Ni,1,k))
Z(Ni,j,k)=Z(Ni,1,k).(22)
无钻孔煤层区域即为组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)。
{X(Ni,j,k)=X(Ni,1,K)
Y(Ni,j,k)=Y(N1,j,k)
Z(Ni,j,k)=Z(N1,1,k).(23)
完成巷道与钻孔的整体变形后,顶部区域仍然存在网格不规则的情况,如图3(i)。按照巷道钻孔交叉区域,巷道区域,钻孔区域及无钻孔煤层区域的顶部修正顶部4个组,如图3(j)。节点Ni,j,k坐标如式(24)。
{X(Ni,j,k)=k/(K2)·(X(Ni,j,K2)-X(Ni,j,1))+X(Ni,j,1)
Y(Ni,j,k)=k/(K2)·(Y(Ni,j,K2)-Y(Ni,j,1))+Y(Ni,j,1)
Z(Ni,j,k)=k/(K2)·(Z(Ni,j,K2)-Z(Ni,j,1))+Z(Ni,j,1).(24)1.7 镜像命令补全模型
经过以上6步,建立了1/4巷道钻孔交叉模型,使用FLAC3D中的镜像命令,即可完成整个模型的建立。
复杂模型主要指倾斜钻孔模型及多孔模型。
倾斜钻孔的网格变形建模首先根据需要建立相应的平行钻孔,如图4(a),接着对钻孔部分进行网格变形。
经过平行钻孔7步后,对钻孔区域,无钻孔煤层区域及两其顶部区域,即组3,组4,组7,组8沿Y轴及Z轴进行平移,得到所需钻孔角度,如图4(b)。节点Ni,j,k坐标如式(25)。
{Y(Ni,j,k)=Y(Ni,j,k)+tanρ1·(X(Ni,j,k)-R)
Z(Ni,j,k)=Z(Ni,j,k)+cotρ2·(X(Ni,j,k)-R).(25)
式中 ρ1为钻孔在XOY平面内投影与X轴的夹角; ρ2为钻孔与Z轴夹角。
建立单孔模型后,对平行钻孔可以直接使用镜像命令复制,形成多孔模型。而对于倾斜钻孔,由于其模型外形不规则,不能使用前述方法。采用导出移动导入法,能够较好的解决此问题。首先使用export命令导出模型,其次对模型进行全坐标移动,移动后的坐标如式(26)。再次使用import命令导入模型。最终得到多孔模型如图4(c)。
Y(Ni,j,k)=Y(Ni,j,k)+P.(26)
式中 P为模型沿Y轴移动量,大小等于模型的Y方向尺寸。
建模对象为近水平煤层,瓦斯抽采钻孔使用斜向孔,孔长为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]是相吻合的,模型能够较好的反应巷道钻孔交叉情况下的力学状态。
针对FLAC3D建模复杂的问题,提出基于FISH语言的巷道钻孔交叉网格变形建模法,编制了用于采矿工程中的巷道钻孔交叉结构的模型建立程序。该方法及理论直白易懂,操作简单易行; 该程序实现了建模的自动化、参数化以及可重用化。减小了建模的时间与难度,提高建模的精度与效益。通过应用实例验证了所述方法的可行性与有效性,有助于帮助研究人员深入了解复杂建模的本质,解决了巷道钻孔交叉建模的难题,在采矿工程领域FLAC3D模型建立方面具有借鉴意义。