王文旭, 李示波
(1.中国冶金矿业总公司,北京 100011;2.河北联合大学,河北 唐山 063009)
河北省的冀东地区,地表地形多为山脉、丘陵,属于中低山区,中等切割,地形复杂,山峦起伏。山区内铁矿资源较为丰富,其矿床类型为鞍山式沉积变质铁矿。地下铁矿体埋藏深度不一,埋藏较浅的矿体上部先露天开采,随着开采深度增加目前多数转为地下开采。有用矿体被采出以后,开采区域周围的岩体原始应力平衡状态受到破坏,应力重新分布,采空区围岩崩落,岩层将移动变形。此影响如果发展到地表,将产生连续或非连续的地表变形,引起一系列环境岩土工程问题,给矿区人民生产生活带来巨大的影响和损失,同时给矿区未来工程建设留下巨大隐患。由于地下开采对岩体造成扰动,其对上覆岩体稳定性的影响可诱发上部边坡体的滑移,会引起山体受采动影响稳定性降低,对山体边坡安全构成危害;或引起山体变形,对山体上的构筑物产生破坏[1]。即使是地下较深的采空区,地表也会引起显著移动变形[2]。在山区,地下开采次生应力场对山体边坡岩体力学场的影响机理的研究还很少。因此,开展山体下采空区地表沉陷变形的研究,对于进一步认识开采沉陷基本规律,为矿层的合理开采及上部边坡的安全性评价提供依据,有着重要意义。本文以唐山北部山区某铁矿地下开采的实际情况为实例,分析了地下开采对山体边坡的影响。
矿区位于华北地台燕山台褶带马兰峪复式背斜迁西穹褶皱西部。矿区出露地层主要为下太古界迁西群马兰峪组片麻岩系,主要为中细粒紫苏黑云角闪斜长片麻岩、角闪斜长片麻岩、黑云角闪斜长片麻岩等。其次为南北两侧出露的长城系常州沟组地层,其岩性为底砾岩、砾岩、石英砂岩、石英岩。第四系地层分布于山间沟谷地带,主要为冲积洪积层,次为残坡积层及范围不大的坡积层。矿体埋藏于山坡下,为一隐伏矿体,倾角60°~65°,矿体较厚。矿石自然类型为石英岩型磁铁矿石,矿石呈粒状变晶结构或粒状花岗变晶结构,条带状、片麻状构造。
矿体上盘山坡上有一条供电线路,采用混凝土电线杆架设。地下开采引起的地表移动变形,会引起供电线路倾斜,倾斜变形严重时,可引起供电事故。因此,有必要对地下矿体开采引起的山体变形进行计算分析,确定其影响的程度,论证其变形是否在可接受的范围,以便能够及时采取必要的对策措施,合理开采。
地表变形的具体情况,可通过FLAC软件进行数值模拟。FLAC主要适用于模拟计算岩土体的力学变形情况和岩土体达到屈服极限以后所产生的塑性流动情况。其所采用的显式拉格朗日快速算法,特别适合模拟大变形和扭曲,能使计算结果更趋于准确。FLAC为解决工程地质问题提供了强有力的理想工具[3]。本文中的地下铁矿,矿体和山脉走向近乎平行,是较为典型的平面应变问题。因此,采用二维数值模拟,能够模拟出地下开采对上部山体的采动影响效果,对地下开采引起的山体边坡岩体和采场围岩力学变化规律以及稳定性进行分析研究,以便给工程施工提供理论依据,减少次生灾害发生。首先利用FLAC建立开采模型,垂直矿体走向取一个剖面,对模型施加边界条件并进行求解后,获取各个单元的弹塑性变化图和位移等值线图。
矿体上盘存在供电线路,地表架设有电线杆,地表变形对它的影响主要是地表倾斜变形的影响,这就需要计算地表下沉引起的倾斜大小,也就是需要求出这一点下沉变形的一阶导数。
由于有限差分是按照一定网度划分出网格,即FLAC得出的数值解是一些离散的点的数值。对离散点求导数,显然不能用连续函数的求导公式进行计算,此时,使用数值微分方法比较方便、准确[4]。
由数值计算方法可知,若函数f(x)由表格给出,求f(x)在结点上的微商,通常称为数值微分。最简单的数值微分公式是用差商代替微商。
(1)
显然,从理论上看,h越小,微分近似越精确。但从计算误差上看,步长h越小,f(x0)与f(x0+h)越接近,根据误差理论,两个相近的数相减,将会大量丢失有效数字。为了克服上述弊病,一个很自然的想法是求式(1)误差渐进展开式,再利用外推法提高精度。根据Taylor公式:
(2)
和
(3)
两式相减,除以2h,移项得:
(4)
令G(h)为上式等号左端的第二项,在Richardson外推法中,取q=1/2,则对于一阶导数可建立外推公式:
(5)
(m=1,2,…)
其中:f′(x)-Gm+1(h)=O(h2(m+1)) 。
根据此外推公式,在步长选择较为合适时,可获得较高的精度[6]。
计算时,步长可取为测点间距,即可按照所得的Gm+1(h)值计算出相应的一阶导数,其中m=1,即外推一次即可。忽略高阶无穷小,相应的Gm+1(h)即可取为一阶导数。
FLAC模拟的竖直位移等值线,供电线路附近取5个点,记录各点的竖直位移y值及水平位移x。由模拟的历史记录点得到各点竖直位移见表1。用数值微分计算由于数值位移引起的电线杆附近地表倾斜变形,见表2。
表1 地表竖直位移
表2 相应点一阶导数
同时,由于山体边坡倾角较大,各点水平位移x的差异也会引起边坡的倾斜变形,因此水平变形引起的倾斜值应一并计算在倾斜变形内。计算方法仍选用上述的数值微分方法。在FLAC历史记录点上,取各点的计算后水平位移,计算后得水平位移引起的地表倾斜变形值为-0.0016。则电线杆附近地表倾斜变形值为0.0051。
由上面的模拟和数值微分计算结果可以看出,地表供电线路位置倾斜变形最大值为0.51%。根据《架空送电线路运行规程》规定,线杆横线路方向的最大偏斜为1.0 %的最大允许限差,上述的计算最大倾斜变形相对于高压线路是横线路方向,其值小于1.0 %。供电线路处于安全状态。
该矿山于1999年12月开始对山坡下矿体进行开采,经过1年多的开采,地表未产生大的变形,农田和树木生长正常,未收到采矿影响;并且地表变形未能引起供电线路有较大倾斜变形,运行正常。
1)应力、应变分析求解岩体移动变形的方法是目前较为严格的分析方法之一。对于FLAC数值模拟结果,进一步通过数值微分方法进行分析处理。数值微分方法能够将该点附近的许多点的模拟数据一并用于计算,可获得更为精确的地表变形结果。通过求解一阶导数可计算出倾斜变形大小,还可以采用数值微分求解二阶导数计算出地表的弯曲变形,为安全评价提供定量的分析依据。
2)在山坡地形,地下开采对地表倾斜变形的影响,不仅取决于竖直位移引起的倾斜变形,还与水平位移引起的倾斜变形有关,两个方向的位移引起的倾斜变形应一并计算在内。
[1] 纪万斌.塌陷与灾害[M].北京:地震出版社,1998.
[2] 韩放,谢芳,王金安. 露天转地下开采岩体稳定性三维数值模拟[J]. 北京科技大学学报, 2006,28(6):509-514.
[3] Wen-Xiu Li,LeiWen,Xiao-Min Liu.Ground movements caused by deep underground mining in Guan-Zhuang iron mine[J]. International Journal of Applied Earth Observation and Geoinformation, 2010,12:175-182.
[4] 黄志安,张英华,李示波,等. FLAC在确定沙曲矿裂隙带上下界中的应用[J]. 矿业研究与开发,2006,26(1):21-23.
[5] 黄志安,李示波,赵永祥,等. FLAC和数值分析在矿山地表沉降预测中的应用[J]. 有色金属,2005,57(3):95-98.
[6] 刘钦圣, 张晓丹, 王兵团. 数值计算方法教程[M]. 北京: 冶金工业出版社, 1998.