爆破作用下岩溶对边坡稳定性的数值模拟研究*

任高峰,王 鑫,周汉红,邓德志,廖 兵,向小东,胡世士,张聪瑞,SAMI Kakar,THEONESTE Ndayiragije

(1.武汉理工大学 资源与环境工程学院,武汉 430070;2.武汉市公安局 治安管理局,武汉 430077;3.中交四航局第二工程有限公司,广州 510230;4.葛洲坝兴山水泥有限公司,兴山 443000;5.中国葛洲坝集团 水泥有限公司,荆门 448000)

露天金属矿山大多采用爆破的方式进行开采,但是,在爆炸过程中,炸药产生的一部分能量会以地震波的形式向外传播[1],对爆区范围内的岩体安全、边坡稳定性造成影响,特别是边坡处于岩溶发育区时更容易造成边坡失稳[2]。

数值模拟计算因其成本低、求解速度快备受国内外研究人员青睐,众多学者也采取数值模拟的方法开展各类研究[3,4]。杜学领等通过建立包含岩溶结构的开采模型[5],研究岩溶结构对开采扰动、开采速度、开采方向等方面造成的影响,得出岩溶结构造成顶板应力突变的位置在岩溶与非岩溶边界区域更加明显。黄静波等运用MIDAS/GTS数值模拟软件对溶洞围岩中爆破地震波的传播规律进行分析[6],得出岩溶隧道受质点振速而产生溶洞围岩扰动,而避免产生裂隙的安全振速范围是小于6 cm/s;王少辉等运用MIDAS软件进行数值模拟得出在爆破荷载作用下[7],隧道溶洞顶壁、隧道拱顶以及拱脚更容易发生应力集中,同时这些部分更容易出现不稳定岩体掉落和围岩裂隙扩展现象,并且在爆破工艺中,由于工艺所造成的干扰损伤[8],往往会在溶洞顶壁和隧道的轮廓交会处,而爆点受力集中的可能性就越大,岩壁裂隙扩展趋势也随着质点振速的增加而增加。

综上所述,部分学者在研究岩溶区对工程结构造成的影响方面已取得了一系列成果[9,10],在前人研究的基础上,论文以兴山县塘垭石灰岩矿+1014 m平台出露岩溶为研究对象,采用ANSYS/LS-DYNA有限元分析程序进行数值模拟,建立现场炮孔布置方式计算模型,并分别在模型有无溶洞位置设置应力监测单元,重点对比分析了有无溶洞对边坡应力[11]、爆破振速[12]、边坡岩体损伤的差异[13],并提出了岩溶处理的措施,为岩溶区的支护设计提供了理论基础和技术支撑。并提出有效的支护处理措施,增强岩溶区边坡岩土体的稳定性,保障边坡的长期稳定[14]。

1 工程概况

塘垭水泥岩矿区位于兴山县城西南约28 km处,属兴山县昭君镇青华村管辖,矿区面积16 000 m2。矿山为山坡露天矿山,设计范围内最高开采标高+1062 m,最低开采标高+880 m,根据矿山地形条件、矿层赋存情况结合矿区开采现状,设计开采方法为自上而下分层开采,台段高12 m。如图1所示。

图1 塘垭水泥灰岩矿露天边坡Fig. 1 Open-pit slope of Tangya cement limestone mine

一般矿山地表岩溶裂隙率是0.05%～1.3 %,而塘垭矿山最大达到了9.23 %,塘垭矿山矿区方向的岩溶率0.61%,其岩溶裂隙在深部钻孔中最大长度达到了7.8 m,矿区平均的岩溶裂隙率2.84%,说明矿区整体岩溶不发育。三峡大学在经过调查评估后,证实前期关于岩溶的结论基本正确。即塘垭矿山岩溶总体不发育,该地区岩溶也为边坡的滑塌提供了一定的自由面。

矿山从坡脚+1014 m平台到坡顶+1132 m平台共发现12个溶洞,呈溶缝或落水洞型,宽度在1 m以上的4个,占30%,宽度在2 m以上的4个,占30%。

坡体内分布有岩溶裂隙和小型溶洞,主要集中于中北部,部分裂隙贯通至地表或坡面,这些裂隙提供了一条通道,使大气层的降水渗透到补给地下水中,使岩体的力学性能发生了劣化。塘垭石灰岩矿深部扩界开采过程中可能遭遇岩溶较发育地段,存在诱发局部崩陷可能,增加了矿山日常生产可能产生的一系列不安全因素。因此,研究岩溶对爆破地震波的传播影响,判定隧道爆破过程中空洞效应对围岩的影响程度以及影响范围,以指导边坡稳定性的设计尤为关键[15],采用爆破动荷载数值模拟手段,研究在爆破工况下岩溶对边坡稳定性的影响,为岩溶附近区域的边坡岩体设计支护提供有效的理论依据。

2 数值计算模型及参数

2.1 模型的构建

根据边坡的地形平面图及勘探线地质剖面图,进行三维坐标变化,建立地质体之间的关系,并按照边坡开挖顺序,采用1∶1的比例处理成三维实体模型如图2所示,建立该模型重点是基于主要断层对边坡稳定性的影响。

图2 三维实体模型Fig. 2 Three-dimensional solid model

在数值模拟模型的建立过程中,首先对该模型进行以下几个假设:(1)假定岩石介质具有连续性;(2)假定爆炸产物的膨胀是绝热工艺;(3)爆轰气体的渗流效应不在考虑范围。

为了简化模型,减少单元数量,提高计算速度,选取高边坡北侧+1026 ～+1040 m平台结构断面为原型,并依照ANSYS数值模拟软件的基本原理,在实际计算中将z方向长度设置为1个单位进行建模,然后在后处理当中设置z方向位移约束。

探究爆破对边坡的影响,假设z方向无限大。所建立的模型大小为13 800 cm×235 800 cm×1 cm,炮孔半径为4.5cm,排间距为380 cm,假设溶洞为空气填充的均质空材料,尺寸为300 cm×300 cm×1 cm数值模型采用单位制。台阶后部,底部和z方向设定无反射边界。

2.2 爆破荷载

1)爆破介质模型

为了研究爆破荷载作用下的岩石材料变形和破坏机理,本文在LS-DYNA中采用比较具有代表性的*MAT_Johnson_Holmquist_Concret岩石材料模型,简称H-J-C模型,其在爆破作用下对岩石与混凝土材料的高应变和大变形的物理性能描述较为准确,岩石物理力学参数如表1所示。

表1 岩石物理力学参数Table 1 Physical and mechanical parameters of rock

2)炸药模型

LS-DYNA中包含两种炸药爆轰模型:高能燃烧模型和点火生成模型。本次模拟用到高能燃烧模型(MAT_HIGH_EXPLOSIVE_BURN)。高能燃烧模型根据炸药上各点距起爆点的距离和炸药爆速来确定每点的起爆时间。炸药参数如表2所示。

表2 炸药参数Table 2 Explosives parameters

3)状态方程

在高压状态下,固体通常会表现出流体特性,故更新静水压力需要引入状态方程。用状态方程来表述压强、内能和密度等参数之间的关联。材料模型只给出塑性阶段材料的本构,状态方程给出了高压状态下弹性材料的本构。本次模拟选用JWL方程,该状态方程是目前模拟计算中使用最为广泛的一类,属于不显含产物组成的状态方程。

1956年,美国工程师Lee根据H Jones和M L Wilkins的研究理论上提出了一种新的状态方程,即JWL方程,该方程是一种半经验状态方程。JWL等熵方程形式如下

Ps=Ae-R1V+Be-R2V+CV-(ω+1)

(1)

式中:V=v/v0=ρ0/σ为相对比容,其中v代表现时比容,v0代表初始比容;A、B、C以及R1、R2、ω分别为材料参数,在计算中材料参数根据实际情况进行取值。炸药状态方程(JWL)的参数如表3所示。

表3 炸药状态方程(JWL)的参数Table 3 Parameters of the Explosives Equation of State(JWL)

模拟所用参数的数据,一部分来自于矿山的地质资料,还有一部分来自于具体的推演计算。在进行下一步工作,即将各项参数在后处理中赋值的过程前要进行单位换算,将原单位换算至cm-g-μs单位制中,便于之后读取数据。

2.3 监测点的布置

将监控点布置在LS-Prepost后处理软件中,是得出数值模拟结果的一种常用方法,在选定监测点后,输出各监测点的数据进行理论分析,由局部分析总体,从而达到理论分析的目的。监测点的布置要遵守一定的原则,使其能够更好地反映总体的模拟结果。

在考虑距离炮孔方面,监测点要选择距离炮孔位置适当的地方,由于爆破模拟属于大应力模拟,故在距离炮孔过近或者较远的地方布置监测点都不能反映较好的模拟结果。

在考虑自由面方面,监测点不应选在自由面附近,由于自由面会对应力波产生折射和反射,因此,靠近自由面的地方可能会发生较大的震荡。

本次模拟主要是为了验证边坡爆破过程中,溶洞对附近边坡造成的影响,因此考虑上述选点原则,应将监测点设置在溶洞附近,远离自由面的地方,具体选点位置如边坡模型构建如图3、图4所示。

图3 有溶洞和无溶洞模型建立Fig. 3 Establishment of model with and without karst cave

图4 有溶洞和无溶洞边坡监测点提取单元位置示意图Fig. 4 Extraction unit location diagram of slope monitoring points with and without karst cave

3 模拟结果及分析

3.1 岩溶对边坡应力影响分析

应力云图表示的是不同时刻下应力的分布,为了使云图显示更清楚,在后处理中设置显示的最大值为7 MPa,即在上图显示中应力最小值在红色区域处,为7 MPa。

由图5可知,在11 ms时应力波传到了溶洞附近,并在溶洞尖端位置产生了应力集中现象,而同一时刻的无溶洞边坡中,应力的传播是规律的,在同位置处并没有产生上述的应力集中现象。

图5 第11 ms时边坡应力云图Fig. 5 Slope stress nephogram at 11 ms

从无溶洞第11 ms时边坡应力云图中可以看出,在无溶洞边坡中,第一次爆破产生的应力波已经穿过了观察区域,而在有溶洞第11 ms时边坡应力云图中,由于溶洞的存在,溶洞尖端再次出现了应力集中现象,且应力值不小于7.6 MPa。

由图6可知,观察溶洞区域边坡受力情况,在15 ms时刻,在完整边坡第二次爆破产生的应力波到达观察区域之前,含溶洞边坡已经在溶洞附近产生了应力集中区,且应力值已经达到的所设置的显示阈值。通过观察后续云图,该区域应力仍然在不断增加。

图6 第15 ms时边坡应力云图Fig. 6 Slope stress nephogram at 15 ms

综上所述,溶洞的存在,使得台阶出现受力改变,出现应力集中,根据前期开展的岩石单轴抗压强度实验,所取岩样平均单轴抗压强度是94.71 MPa,所以台阶受力产生的应力值较小,仍不足以破坏局部边坡岩体。

3.2 岩溶对爆破振速影响分析

选取的监测点在空间位置上两两近似相同,根据其是否穿过溶洞对其进行编号,含溶洞边坡三个监测点分别为A1、B1、C1,完整边坡相对应位置的三个监测点分别为A0、B0、C0。

为方便描述以相同空间位置的两个监测点为一组,分为第一组(A1、A0)、第二组(B1、B0)、第三组(C1、C0),并通过计算,分别提取其在x和y方向上的振速时程曲线,如图7～图12所示。

图7 A点x方向振动速度Fig. 7 Vibration velocity of point A in x direction

图8 A点y方向振动速度Fig. 8 A point y direction vibration velocit

图9 B点x方向振动速度Fig. 9 Vibration velocity of point B in x direction

图10 B点y方向振动速度Fig. 10 B point y direction vibration velocity

图11 C点x方向振动速度Fig. 11 C point x direction vibration velocity

图12 C点y方向振动速度Fig. 12 C point y direction vibration velocity

三组监测点的x方向最大振动速度在2.2 cm/s上下,且有无溶洞情况下并没有明显差别。其中第一组和第三组在完整边坡下的x方向振动速度略大于含溶洞边坡,比较第一组和第三组数据发现,A0和C0之间的振动速度相差约为0.15 cm/s,而A1和C1两监测点之间的振动速度差值约0.1 cm/s,分析空间位置关系,可以认为,溶洞的存在阻碍了震动波的传播,导致上述振动速度差值的变化。

就y方向振动速度而言,组内比较第一二组,振动速度差值约为0.3 cm/s,第三组的差值约为0.7 cm/s,含溶洞边坡监测点的振动速度均大于同监测点完整边坡的振动速度。溶洞的存在增加了新的自由面,尤其是y向上,空间位置上本来就存在较近的台阶坡面,振动速度明显高于x方向因此,溶洞对振动速度的影响,在y方向的表达更为明显。第一三组的A1与C1的差值为1.01 cm/s,而A0和C0点的差值为0.63 cm/s,这可能是由于溶洞的存在,增加了自由面,使得前者的空间位置不同导致的振动速度差值明显增加。

综上所述,有无溶洞的存在,对相同空间位置监测点x方向爆破振动速度的影响较小,但不同空间位置的振动速度,仍然受到影响,具体表现为,溶洞的存在使得空间位置的不同的监测点之间的振动差值变大。而y方向振动速度受溶洞存在的影响则表现更加明显,具体表现为,同一空间位置的监测点,由于溶洞的存在,产生了新自由面,加剧了该点的震动;不同空间位置的监测点,由于溶洞的存在,振动速度的差值变大,尤其是空间位置靠近边坡坡面的监测点,对这种情况的反应更加明显。

溶洞的存在,使溶洞两侧的爆破振动速度出现差异,具体分析结果:表明溶洞的存在,较无溶洞的同一空间位置监测点振速会出现差异,其中y方向振动速度受溶洞存在的影响则表现更加明显;由于振速的差值较小,故对边坡的稳定性影响不大,在距离爆源相对较远的情况下,按照《爆破安全规程》(GB6722—2014)规程,该边坡爆破振动安全允许质点振动速度是8～12,将模拟结果与上述安全范围相对比,能够看出所有数据显示的振速峰值均在安全允许振速的范围以内。

3.3 岩溶对边坡岩体损伤影响分析

在ls-prepost软件中,损伤云图是描述岩石破坏最直接的表示方式。由于岩石成岩的是一个复杂的过程,初期破坏是不可避免的。含有初始损伤的脆性岩石在冲击载荷的作用下将会产生两方面的影响:一方面会导致材料刚度的劣化(等效模量的降低);另一方面也会使应力波能量耗散。

根据应变等价理论和能量理论,一般认为初始破坏为零,从而分析岩体破坏特征。但是对于相对完整的岩石而言,由于天然裂隙的存在,裂隙岩体实际上是具有初始损伤的材料,对其进行损伤演化规律研究时需考虑初始损伤。因此将损伤显示上限设置为0.1,得到不同时刻下模型的损伤云图如图13、图14所示,由于设置了损伤上限为0.1,故图中红色部分表达的意义是,该区域最小损伤为0.1。

图13 边坡岩石19 ms的损伤情况Fig. 13 Damage of slope rock in 19 ms

图14 边坡岩石22 ms的损伤情况Fig. 14 The damage of slope rock in 22 ms

从图13和图14可以看出,在爆破次数不断增加的过程中,岩体内部的损伤也在逐步加剧,且损伤在有溶洞的影响下更为明显,岩溶区半径越大,爆破次数越多,损伤就更为严重,且在实际模拟过程中发现,有溶洞区域爆破可能会导致临近岩石区域产生空洞,加剧岩体损伤。

3.4 岩溶对边坡有效应力影响分析

依照模拟方案,在分析有效应力时,利用后处理器Ls-prepost在模型中设置的监测点与振动速度的分析相同,以相同空间位置的两个监测点为一组,分为第一组(A1、A0)、第二组(B1、B0)、第三组(C1、C0),分别分析完整边坡和含溶洞边坡有效应力变化情况,如图15～17所示。

图15 溶洞边坡监测点有效应力(第一组)Fig. 15 Effective stress of karst cave slope monitoring points(first group)

图16 溶洞边坡监测点有效应力(第二组)Fig. 16 Effective stress of karst cave slope monitoring points(second group)

图17 溶洞边坡监测点有效应力(第三组)Fig. 17 Effective stress at monitoring points of karst cave slope(third group)

由于自由面的存在,有效应力的分布也显现出震荡,这种震荡在含溶洞边坡的B1和C1监测点表现得更为明显,分别读出各监测点的最大有效应力值,列表如表4。

表4 不同监测点有效应力峰值Table 4 Effective stress peaks at different monitoring points

根据表4中数据,有溶洞边坡模型的A1、B1、C1单元处的有效应力峰值均大于同组的无溶洞边坡模型的监测点的值,产生的应力差值大于7.6 MPa,原因是由于溶洞的存在,溶洞尖端会出现多次应力集中现象,根据前期开展的岩石单轴抗压强度实验,所取岩样单轴抗压强度平均值为94.71 MPa,因此模拟得出的有效应力峰值总体均低于岩体的单轴抗压强度。

4 岩溶支护处理措施分析

为保障边坡的长期稳定,消除安全隐患,必须根据实际情况,结合溶洞对边坡稳定性的影响,通过有针对性的支护处理措施,增强岩溶区边坡岩土体的稳定性,常用的措施有[16]:

(1)在岩溶边坡开挖施工过程中,为了降低坡体的承载能力,可以对坡体进行削减,将坡体的高度和坡度降低,并且在施工过程中,削减的原则是对相对稳定的坡体进行保留,对不稳定的坡体进行削减。

(2)对于岩溶边坡不稳定的岩体,可以通过修建挡土墙、护墙等人工加固边坡方式,从而提高边坡的稳定性。对于低抗滑性能的岩体,可用钢筋水泥防滑桩进行加固岩体。

(3)对于有裂隙和结构面软弱的岩质边坡,可采用预应力锚杆或者锚索和裂缝注浆等方法对其进行加固处理,从而增强裂缝的抗剪强度。

根据本论文数值模拟结果,对岩溶边坡进行有针对性的治理,在溶洞尖端产生应力集中的地方采用喷锚网支护,即“锚杆(锚索)+钢筋网+喷射混凝土”。利用锚杆(锚索)锚固于溶洞尖端部分,用于岩体整体稳定性的加固;利用钢筋网和喷射混凝土对溶洞边坡表层松动的岩块进行加固,该联合支护方法对增强整个岩土体的稳定性比单独使用锚杆(锚索)支护、喷射混凝土支护效果更好[17],从而确保工程安全。

5 结论

采用ANSYS/LS-DYNA数值模拟软件,建立现场炮孔布置方式计算模型,分别对有无溶洞的爆破过程进行数值模拟计算,研究爆破作用下岩溶对边坡应力、爆破振速,边坡岩体损伤,边坡有效应力的影响,得出以下结论:

(1)在爆破作用下,岩溶区的岩体受力发生改变,产生应力集中,这种现象在溶洞尖端部位表现得更为明显。

(2)溶洞的存在,会产生新的自由面,加剧岩体的震动,且溶洞对振动速度的影响,在垂直方向的表达远超水平方向。

(3)岩体损伤范围与岩溶区尺寸变化规律相同,岩溶区半径越大,损伤就更为严重。因此在实际施工中,可适当填充溶洞以减小岩溶区域的范围来削弱安全隐患。

(4)为保证溶洞作业安全,可采用“锚杆(锚索)+钢筋网+喷射混凝土”建立溶洞作业安全防护体系从而确保工程安全。