寇向宇,贾明涛,王李管,吴 霞

(1.中南大学资源与安全工程学院,湖南 长沙 410083)

(2.中南大学数字矿山研究中心,湖南 长沙 410083)

0 引言

FLAC3D是一种基于连续介质快速拉格朗日差分法的三维空间分析软件,在采矿与岩土工程中广泛应用。但 FLAC3D的前处理能力不强,在建立计算模型时主要采用键入数据/命令行文件方式,虽然可以采用内置 Fish语言编程实现复杂模型的建立,但其独特的源代码表达方式、加之实际采矿工程介质和工程形体的复杂性,由 FLAC3D直接构建复杂模型十分困难[1～2]。

三维矿床建模软件可以利用地质钻孔、平剖面图或其他测量数据,生成复杂地质和工程对象界(表)面的 T IN模型,并在此模型控制下,将其离散为可以逼近形体表面的规则六面体单元,从而用于对地质体属性(如品位)等的估值。

侯恩科等[3～4]介绍了将矿床建模软件与FLAC3D软件进行耦合的应用实例,但由于六面体单元在 2类软件中的应用目的不同,关于如何在块段建模过程中对模型进行优化以使之满足力学分析要求,在这些文献中并没有详细讨论。

本文从 FLAC3D网格划分的原则出发,根据矿床建模软件 SURPAC构建块段模型的方法,对用于力学分析的块段模型构建的步骤、原则等进行研究,使建立的模型得到优化,使其能够直接作为 FLAC3D分析的前处理模型进行模拟计算。

1 FLAC3D网格划分原则及 SURPAC块段模型构建方法

1.1 FLAC3D网格划分原则

有限元网格划分是进行数值模拟分析至关重要的一步,它直接影响着后续数值计算分析结果的精确性。网格划分涉及单元的形状及其网格的密度、数量、质量,单元阶次,网格分界面和分界点,位移协调性,网格布局,节点和单元编号等[5]。

网格数量的多少将影响计算结果的精度和计算规模的大小。网格疏密是指在结构不同部位采用大小不同的网格,这是为了适应计算数据的分布特点。在计算数据变化梯度较大的部位(如应力集中处),为了较好地反映数据变化规律,需要采用比较密集的网格。而在计算数据变化梯度较小的部位,为减小模型规模,则应划分相对稀疏的网格。选用高阶单元可提高计算精度,因为高阶单元的曲线或曲面边界能够更好地逼近结构的曲线和曲面边界,且高次插值函数可更高精度地逼近复杂场函数,所以当结构形状不规则、应力分布或变形很复杂时可以选用高阶单元。但高阶单元的节点数较多,在网格数量相同的情况下由高阶单元组成的模型规模要大得多,因此在使用时应权衡考虑计算精度和时间。网格质量好坏将影响计算精度。结构中的一些特殊界面和特殊点应分为网格边界或节点以便定义材料特性、物理特性、载荷和位移约束条件。当结构形状对称时,其网格也应划分对称网格,以使模型表现出相应的对称特性(如集中质矩阵对称)。不对称布局会引起一定误差。节点和单元的编号影响结构总刚矩阵的带宽和波前数,因而影响计算时间和存储容量的大小,因此合理的编号有利于提高计算速度。

1.2 SURPAC块段模型构建方法

块段模型是一种分割矿体、模拟矿体内部特征的有效数据模型,就是对目标地质体进行三维栅格化处理[6],把矿体划分为一系列小长方体单元,其方法是,首先将研究的范围形成最小包络长方体,将其定义为原型；然后根据地质勘探网度、采矿方法、地质条件以及地质统计学等方面对块段要求来确定单元块尺寸,以此对原型进行三维栅格化；然后对地质体进行块段划分时,应选用最合适的、且最具代表性的矿体平剖面图建立地质体表面模型(包括矿体、地层和断面等模型)对实体模型进行边界约束,在边界处进一步细分,以逼近地质体的空间形态；每个小长方体单元的内部属性相同,且该属性代表了矿体该位置的内部特征,这样所有小长方体单元的属性变化规律就是整个矿体的特征变化规律。因此,在建立好块段模型之后,就可以进行相应的统计分析工作了。

2 用于力学分析的 SURPAC块段模型优化建模技术

为使 S URPAC建立的块段模型能够直接用于FLAC3D力学分析,保证力学分析的效率和精度,根据前文分析,块段模型建立时应遵从以下步骤和方法。

2.1 模型分区

根据网格划分原则,为了适应计算数据的分布特点,较好地反映数据变化规律,应力可能集中的地方应密,其他地方可稀疏一些。网格数量也要考虑精度和计算规模而确定。因此,考虑增加网格的经济性和分析数据的类型,按照不同的划分方法和力学参数建立原始块段模型,并根据三维模型进行约束,最后,利用编制的接口程序,将块段模型转换为 FLAC3D分析前处理模型,然后对比不同网格划分方法的模拟结果,以确定一个合适的值。这就是为满足疏密、密度等的要求,在建模时,可根据实际情况分成相关区段,外围区段单元大、往里逐渐减小,只对最里面有形状变化的区段细分。这样既保证计算精度,又保证单元数量不会太多,且虚拟节点不会很多,找不全的问题,进而提高了分析效率,实现了网格划分及其优化。

2.2 边界单元细分

边界单元细分一是为逼近几何边界、二是达到疏密网格的目的,从而反映应力集中区的应力变化。在通过 SURPAC软件输出网格时,要进行块段模型的建立,在建立的过程中,块段模型选取地质体线框模型进行约束,块段模型的单元块在边界处将自动按照一定的准则发生细分(变块技术)。同时在建模时细分次数也可超过 2次以上,导入 FLAC3D中使一个单元跟相邻单元的边长之比小于 1∶1或 1∶2,这样对复杂形态边界的逼近能力更高,效果更好。另外,在细分后导入 FLAC3D必须利用Attach FAce或 Gen merge命令找出并处理这些虚节点。因为当模型较大时,如细分次数过多,会导致接触面虚节点数量多,FLAC3D可能找不全这些点,导致计算结果出错。

2.3 结果验证

为了验证这种网格划分实现前处理模型优化方法的正确性,根据二维弹性力学的原理[7],设计了一个边长 48 m×48 m的板中有一个半径为 3 m圆孔,在左右两边受均布拉力,大小为 1 MPa,然后根据对模型不同的细分方法得出疏密不同的单元体网格,进行分析,再与理论值作比较,最终得出最优的网格划分方法,进而也验证了这种优化方法在力学分析中的重要性和正确性。所以以此方法做出理论应力分布图(11200Zones and 22402 Grid),远近疏密一致的应力分布图 (9988zones and 15448 grid points),疏密程度不一致的应力分布图 (9256zones and 18544 grid points),比较应力的主要分布,显示出不同网格划分对数值模拟效果的影响,如图 1～图 3所示。图中显示,理论上的应力应该集中在中心的两侧,而在不同位置的不同疏密程度模拟出来的效果不同,并不是网格越密效果越好,也不是网格越疏效果越好,只有在应力相对应该集中的地方密一些,其他地方疏一些,图 2在两侧时应力出现了发散,而图 3模拟的效果最接近于理论值,效果最好,因此得出了网格划分的最优选择,图 4为网格划分的显示,进而也得到了前处理模型的优化。

图 1 理论应力分布图

图 2 疏密程度一致应力分布图

图 3 疏密程度不一致应力分布图

图 4 最优网格划分效果图

3 工程应用

3.1 工程概况

云南某铜矿经过数十年的露天及地下大规模工业开采,井下采空区越来越多、越来越大,且许多早期、浅部空区已经发生垮塌,多个中段上下相邻或同一水平相邻的空区发生贯通,甚至部分已经贯通地表,对地表各种工业和社会活动,以及生态环境造成一定的威胁,存在较大的隐患。如果没有有效的采空区稳定性及地表移动规律的评价研究,则采空区垮落不仅加大矿石回采的贫化与损失,也必然对矿井及矿区相关工农业的安全生产造成严重的威胁。为保证矿山有序的正常安全生产以及后续开采的安全,对该矿采空区稳定性及由此引发的地表移动问题进行研究具有非常重要的现实意义。因此,运用FLAC3D数值方法模拟了影响采场稳定性的空区围岩应力及地表位移变化规律,并在模型建立时以不同网格划分方法进行对比,实现前处理模型的优化并进行数值模拟,进而有效地提高了采空区稳定性分析和评价结果的准确性和可靠性。

3.2 力学模拟及其分析

在 SURPAC中,在地表三维模型的约束下,建立原始块段模型。其次,在各中段空区三维表面模型的约束下,建立相应的空区块段模型,然后,把所有块段模型合并,形成一个统一的模型,该模型内单元块数目为 638762块。最后,利用编制的接口程序,将 SURPAC块段模型转换为 FLAC3D分析前处理模型,该模型节点数为 810780个。

在前处理模型中,将各采空区所包含的单元设置为不同的组 (Group),赋力学参数、初始地应力、边界条件,就可以在分析过程中,通过指定某个组的移出,以达到模拟开挖,考察在采空区形成后,周边围岩应力及位移分布情况的目的,最后得出结果。本文选取矿区某中段的空区作为实例进行分析比较,按照相关力学参数进行赋值,做出空区围岩应力分布图和地表位移分布图(图 5、图 6)。

图 5 中段开采后最小主应力空间分布图

图 6 中段开采后地表位移图

图 5为第一步开挖后模型内的最小主应力分布区及其大小。由图 5可见,在采空区周边围岩中出现了拉应力,多在0.3MPa以下,少数单元拉应力较大,接近但小于 0.8 MPa。由图 6可以看出,地表位移量整体上不是很大,在 1.5 cm以内；最大位移发生在露天坑坑底中心位置,随着离坑底中心渐远,整体位移逐渐减小；然后按采场模拟开挖顺序对各中段进行开挖,直至结束,并且分析各中段开挖后模型内的最小主应力分布区及其大小,以及采场中各应力的变化情况以及开采后的地表位移变化,实现对该矿区采空区稳定性的准确模拟,以制定措施对井下空区进行处理,排除矿区生产的安全隐患,保障矿区的正常生活和作业。

4 结论

为了实现数值模拟的效果更加准确,研究了有限元网格划分原则和 SURPAC构建块段模型的方法,然后建立了最优的 FLAC3D前处理模型,并构建了一个实例,依据弹性力学的原理,判定其有效性和准确性,同时也验证此种前处理模型优化方法的正确性,从而实现了 FLAC3D的前处理模型优化。以此在云南某铜矿进行设计实践,通过按照矿区实际开采顺序进行模拟,做出采空区的应力分布和地表的位移移动情况,判断空区的稳定性和地表生产生活的安全性,并制定一系列措施进行治理,排除了矿区生产的安全隐患,保障了矿区的正常生活和安全作业,在实际中验证了此模拟的有效性和精确性,同时也验证了该前处理模型优化方法的可行性和必要性。

[1] 邓红卫,朱和玲,周科平,等.基于 FLAC3D数值模拟的前后处理优化研究 [J].矿业研究与开发,2008,28 (2)：60-62.

[2] 胡斌,张倬元,黄润秋,等.FLAC3D前处理程序的开发及仿真效果检验 [J].岩石力学与工程学报,2002, 21(9)：1387-1391.

[3] 刘科伟,李夕兵,宫凤强,等.基于 CALS及 Surpac-FLAC3D耦合技术的复杂空区稳定性分析[J].岩石力学与工程学报,2008,27(9)：1924-1931.

[4] 侯恩科,吴立新,李建民,等.三维地学模拟与数值模拟的耦合方法研究[J].煤炭学报,2002,27(4)：388 -392.

[5] 杜平安.有限元网格划分的基本原则[J].机械设计与制造,2000：34-36.

[6] 冯兴隆,王李管,毕林,等.基于三维模拟技术的矿岩可崩性评价[J].煤炭学报,2008,33(9)：971-976.

[7] 徐芝纶.弹性力学简明教程 [M].人民教育出版社, 1983：200-203.