基于数值模拟技术的井筒位置优化

吉学文，刘 辉，潘张伟，王 志，扶海鹰

(1.大兴安岭金欣矿业有限公司， 黑龙江大兴安岭市 165100;2.驰宏科技工程股份有限公司， 湖南长沙 410000)

1 工程背景

大兴安岭岔路口钼铅锌多金属矿为斑岩型矿床，主要围岩是花岗斑岩、石英斑岩、闪长玢岩等。钼矿体总体呈拉长的穹窿状，主体隐伏，现已控制矿脉长度2600 m，剖面上控制宽度360～1260 m，最大宽达1760 m;矿体延深大于1500 m，垂直厚度一般200～900 m，最大为1097.48 m。设计生产规模为1650万t/a，矿体按上下分区同时开采，采用大直径深孔空场嗣后充填法和分段空场嗣后充填法自下而上开采，上采区开采范围为+20～+440 m之间的矿体，上采区首采中段为140 m中段，下采区首采中段为－460 m中段，上下采区形成独立的回风系统，1#风井为上采区的专用回风井筒，设计服务年限为20 a。

2 设计方案

按照井巷布置有关规定，若将1#回风井位置设置在全矿区地表移动界线范围之外，距离上采区矿体约为1500 m，同时每个中段回风井石门也长达1500 m。这对于1#回风井仅服务于上采区矿体开采的回风，极大增加了基建投资及通风运营成本。为优化工程布置，根据矿体赋存状况、采矿工艺要求及地表地形，在地表初步确定了回风井的3个位置，3个方案的回风井距离矿体300～500 m，按照采矿设计规范及建构筑物设计要求，通过FLAC3D三维数值模拟方法，对初步确定的回风井3个优选方案在不同时期(15 a，20 a)的稳定性进行验证比较分析。

通过FLAC3D分析研究井筒位置的可行性，并与地表移动界限外井筒位置比较，确定合理可行的井筒位置。

3 数值模拟分析

有限差分法在计算机数值模拟中已被广泛运用。采用拉格朗日法的有限差分数值模拟计算分析软件已广泛应用于岩土工程中，弥补了类比法和解析法难以实现大范围矿床开采对井巷工程稳定性的分析的不足，同时也更好地考虑介质的各向异性、非均质特性及其随时间的变化、复杂边界条件和介质不连续性等复杂地质条件。

3.1 基本假设及模型建立

由于矿山地质结构的复杂性，影响采场稳定因素的多样性，在模型建立和计算过程中采用以下简化和假设:

(1)阶段运输巷道、主井、通风井等大量工程对整个矿山开采稳定性有一定的影响，但这些工程对宏观矿区开采的稳定性影响较小，因此，本次研究不考虑这些工程的影响;

(2)在建模中忽略模拟范围内的节理、裂隙和断层，在模拟过程中，选择具有代表性的3个岩层(风化带岩组，微风化带岩组和未风化带岩组)建立围岩;

(3)在模拟计算中，忽略爆破震动、地震波等动力因素及地下水对矿岩稳定性的影响。

本次模拟过程严格依照设计开采进度计划，先进行风井的开挖，然后以年为单位进行各个中段的开采和充填，记录下开采15 a和20 a对1#回风井的稳定性影响。经过折减后的矿岩物理力学参数见表1。

3.2 本构模型与边界条件

莫尔－库仑模型是一种弹塑性本构模型，它能反应岩土材料剪切破坏特性，并广泛应用于岩土实际破坏分析中，剪切屈服函数和拉应力屈服函数组成了莫尔－库仑模型破坏包络线。

表1 矿岩物理力学参数

采用3Dmine和FLAC3D结合的方式建立数值计算模型。先以位置1进行建模，根据矿区地表等高线模型图及整个矿体的空间形态，模型的计算域长×宽取为3200 m×2300 m，高度取－940 m标高到地表。采用3Dmine软件对矿体实体模型进行块处理，同时加上地表、矿体和1#回风井的约束条件，最终生产模型的FLAC3D命令流文件，在块体大小的划分中，由于模拟范围大，同时开挖步骤复杂、中段高度为60 m，最终选用的最大单元格大小为50 m×50 m ×60 m，次级单元格大小为12.5 m ×12.5 m ×20 m，FLAC3D软件读取命令流后生成计算模型，位置1的模型共235989个单元。

模型四周采用可滚动的轴支座约束侧向变形，底面采用固定支座边界类型约束竖直方向变形，顶面为自由面。

模拟计算主要记录3个不同风井位置条件下应力场、位移场和塑性区的变化过程，通过对比分析选取合理的风井位置。

3.3 模拟结果

以位置1进行分析，主要对模型内的位移变化、应力值变化规律和塑性区分布大小3个方面进行分析。

(1)变形规律分析。开采15 a对回风井拟选1#位置地表造成的沉降约为65 mm，开采20 a的最终沉降值约为85 mm。图1中，地表最大位移区域在15 a后出现向右下角移动的趋势，是由于在模拟开挖过程中，15 a后开挖的80 m和20 m中段的矿体主要集中在风井位置的右下角。距井筒腰部最近的矿体约为50 m，见图2，由于该矿体为独立矿体，开挖跨度小，开采对井筒的影响不大;下采区矿体开采20 a，由于中间还隔有近400 m厚度的未采动区域，由计算结果可知，此条件下下采区的采动对1#回风井的稳定性影响不大;开采15 a和20 a井筒底部的沉降值分别约为80 mm和120 mm。

图1 20 a地表沉降

图2 剖面竖向位移云图

(2)应力值变化规律。井筒开挖，最大主应力出现在井筒的底部，为9.82 MPa，整个开采过程到20 a时，经过开挖和回填，模型达到新的力学平衡，内部应力重新进行分配，最大主应力出现在开采矿体的顶板和矿柱上，如图3所示，最大主应力为15.32 MPa，但是井筒底部受到扰动影响，应力变化为8.77 MPa，出现了一定的应力释放。

(3)塑性区变化规律分析。在整个开采过程中，矿体及围岩的剪切破坏和拉伸破坏都主要集中在开采空区的顶板和矿柱区域，随着开采范围的不断扩大，胶结充填体边缘也出现了部分拉伸破坏，但是靠近1#回风井位置并未出现塑性区。

同理，对其它2个方案进行分析，得到以下数值模拟结论:

图3 最大主应力云图

(1)开采到20 a时，影响1#回风井稳定性的主要因素是上采区的开采扰动，由于上下采区之间存在厚大的未采动区域，极大地减弱了下部开采对回风井的影响;

(2)模拟结果对比见表2，同时根据水平位移，可计算得出3个位置在不同年份条件下的倾斜、曲率和水平变形(见表3)，均符合《有色金属采矿设计规范(GB50771－2012)》要求，结合井筒与矿体的三维空间位置，位置3变形位移量较小，符合国家建构筑物规范，位置3为最佳选择。

表2 3个拟选位置的数值模拟结果

表3 数值模拟计算移动变形值

4 井筒优化结果对比

按照有关矿山设计规定，通风井设计在地表移动线以外，风井距离矿体较远，与风井优化位置3相比，每个中段增加巷道及通风路线1000 m，共有6各中段与通风井连接，增加了工程投资5280万元;增加了通风机功率200 kW，年增加了电耗共1584000度。

5 结论

(1)根据模拟结果，推荐井筒位置3为最佳方案，减少了通风井巷工程量，节约了工程投资及生产成本。

(2)采用FLAC3D软件对井筒位置进行模拟计算，得到了整个矿区的采动应力场及位移场的分布规律，说明了该方法的可行性。

(3)由模拟结果推荐井筒位置3的方案，但数值模拟的计算结果都是在现有的岩体力学参数、原岩应力场及简化的开采工艺过程条件下进行的，由于岩石条件的多变性、岩石力学性质试验取样的局限性、岩体结构面的复杂性，计算所需的岩体力学参数难以准确确定。因此，建议在生产期，设置位移监测点，监测井筒位移的变化。

［1］王发芝.冬瓜山铜矿深部采场充填技术［J］.矿业研究与开发，2008，28(8):4-5.

［2］李克顺，吴 姗，宋卫东.大冶铁矿充填法采场结构参数校核与建议［J］.金属矿山，2010(9):16-22.

［3］王金安，赵志宏，侯志鹰.浅埋坚硬覆岩下开采地表塌陷机理研究［J］.煤炭学报，2007，32(10):1051-1056.

［4］《采矿设计手册》编辑部.采矿设计手册(矿床开采卷)［M］.北京:中国建筑工业出版社，1987，895-899.

［5］GB50771－2012.有色金属采矿设计规范［S］.

［6］吉学文，唐绍辉，李爱兵.某地下矿山采场与围岩稳定性三维有限元模拟分析［J］.矿业研究与开发，2003，23(2):7-9.

［7］乔 兰，王双红，蔡美峰.某地下矿山岩层及地表移动规律的有限元模拟研究［J］.金属矿山，2000(4):23-25.