基于FLAC-3D 砂岩块体超低摩擦鞭梢效应研究*

李利萍，唐 垒，潘一山，邰英楼，张海涛

（1.辽宁工程技术大学 力学与工程学院，辽宁 阜新 123000；2.辽宁大学 物理学院，沈阳 110036）

引 言

矿山资源是我国长期以来高速发展的有力保障，随着我国社会发展需要、经济体系不断壮大，原煤需求量也持续增长[1].现今，深部开采已成为主战场，通常把煤炭开采深度大于等于800 m 的矿井称作深部矿井，迄今为止，我国已经探明煤炭资源量超过5 万亿吨，其中57%的煤炭资源为深部资源.煤炭开采深度逐年递增，其速度最高可达25 m/年.已进入深部开采的矿山有47 座，开采深度最高达到1 501 m.随着开采深度增加，深部岩体应力状态越趋复杂、岩体结构愈加破碎，与浅部特征截然不同，冲击地压等动力灾害更加剧烈和频繁[2-3].促进我国煤炭资源高效开采、安全开采、环境友好型开采，已成为深部开采的首要问题.冲击地压现场观测表明，仍有一些特殊的动力现象无法得到合理解释，如沈阳焦煤顶板事故调查报告中指出：西三上采区702 工作面上出口204～214 m 范围内巷道上帮煤体向下帮煤壁侧整体滑移3.0 m，巷道顶部留有宽1.7～1.8 m、高0.5～0.6 m 的空间，极有可能是煤岩体的超低摩擦鞭梢效应诱发所致.

在国外，Kurlenya 和Oparin 等[4-6]发现并提出了岩体超低摩擦效应，同时还通过试验证明了超低摩擦效应真实存在.在国内，王明洋等[7-8]分析得知了深部岩体介质变形过程中诱发块系岩体超低摩擦现象的机理；吴昊等[9-10]基于块系岩体动力模型，计算得到了工作块体水平位移公式；何满潮等[11]采用二维数字图像技术对块系花岗岩进行了超低摩擦效应实验研究；王德荣等[12]研制出深部岩体动力响应试验系统，并分析了界面摩擦减弱现象；李利萍等[13-14]研究分析了冲击扰动对冲击地压型超低摩擦作用机理及围压与冲击扰动组合作用对超低摩擦效应影响的规律；唐红梅等[15]构建了双裂缝主控结构面危岩的力学模型和断裂力学模型，按照最大周向应力准则，基于危岩自重、裂隙水压力及地震荷载作用，得出了双裂缝主控结构面危岩的断裂稳定系数表达式；钟万勰等[16]从耦合振子振动分析出发，指出了高层建筑顶部构造物最易发生鞭梢效应的条件，及消除或减小鞭梢效应的一种可能性；李桂青等[17]提出用质量泵控制鞭梢效应发生的理论，将适宜的质量泵安装在主体结构顶部便可能从根本上降低或消除这种不利影响；杨佑发等[18]详细分析了鞭梢效应在高层建筑结构抗震中的作用；季金铭等[19]研究高陡坡在爆破作用下动力稳定问题时发现台阶处存在鞭梢效应，并且具有方向性；郜宁静等[20]通过ABAQUS 建立塔柱结构模型，利用动力时程法分析了地震作用下塔柱顶部机房的鞭梢效应.

综上所述，国内外学者在研究岩体超低摩擦效应上多集中于岩体接触界面，鲜有从深部岩体断续、破碎等结构特征出发的.随着开采深度的增加，岩体断续结构愈加明显、越趋破碎，深部岩体在结构上与鞭梢效应现象的外部结构特征极大吻合，故将鞭梢效应外部结构特征用于描述深部开采中某些特殊动力响应现象.本文首次用鞭梢效应补充、细化超低摩擦效应发生机理，通过试验与数值模拟结合的方式去研究深部岩体超低摩擦鞭梢效应发生机制，为冲击地压防治提供参考.

1 模型的建立及前期准备

1.1 理论模型的建立

随着开采深度的逐步增加，岩体结构随之变得更加破碎，在深部岩体中存在诸多具有力学特性的结构（节理、软弱夹层、断层等非连续结构面），因此参照Kurlenya 和Oparin 等[4-6]提出的一维块体模型，我们将深部岩体中由于节理、夹层、断层等形成的块状岩体视作随机质量刚体，并建立了砂岩块体超低摩擦鞭梢效应理论基准模型.块体模型共有5 个力学性质相同的砂岩块体，其中包含四个标准块体和一个工作块体，并按顺序从上至下依次编号为1～5 号块体，其中标准块体尺寸均为100 mm × 100 mm × 100 mm，工作块体尺寸为lmm×lmm ×lmm.若l=100，5 个标准块体依次叠放且以块体3 为工作块体，构成超低摩擦鞭梢效应理论基准模型，1 号块体上表面受垂直方向应力波扰动FV，工作块体左侧表面受水平冲击F，如图1 所示.

图1 超低摩擦鞭梢效应理论基准模型Fig.1 The theoretical model for the ultra low-friction whiplash effect of sandstone blocks

1.2 前期准备

1.2.1 模型网格单元划分

超低摩擦鞭梢效应理论基准模型由5 个尺寸均为100 mm×100 mm×100 mm 的标准块体构成，质量均为2.75 kg，每个块体由5×5×5 个立方体单元构成，其中超低摩擦鞭梢效应理论基准模型的工作块体为块体3，块体之间则用Interface 命令建立接触面，然后在接触面上赋予相应力学性质.具体模拟基准模型如图2 所示.

图2 超低摩擦鞭梢效应理论基准模型模拟Fig.2 The numerical simulation model for the ultra low-friction whiplash effect of sandstone blocks

1.2.2 确定本构模型

本文的研究对象为砂岩块体，在数值模拟软件FLAC-3D 中Mohr-Coulomb 塑性模型是岩土力学通用模型，适合松散或胶结的粒状结构（如土体、岩石、混凝土等）.因此，本文选用Mohr-Coulomb 塑性模型作为本构模型.

1.2.3 边界条件及加载条件的设置

针对FLAC-3D 动力计算部分，模型的边界会反射入射波，进而影响模拟结果真实性.因此，FLAC-3D 提供了静态边界和自由场边界来减少模型边界上波的反射.考虑到实际工况，本文模拟选用静态边界.其中，模型系统受重力、垂直方向应力波扰动（正弦式应力波）和水平冲击组合作用.

1.2.4 单轴压缩试验及力学参数的确定

本文中的砂岩块体试件由阜新孙家湾煤矿采购，将采得的大块岩样运至加工厂严格按照岩石力学试件切割要求进行加工，形成标准柱形试件煤样（直径为50 mm、高为100 mm），再利用YAW-2000 岩石压力试验机（如图3 所示）进行单轴压缩试验.得到相关力学参数后取平均值降低参数误差，同时计算数值模拟所需其他参数后一并列入表内，其中表1 为砂岩块体模型力学参数，表2 为砂岩块体间接触面参数.

表1 砂岩块体力学参数Table 1 Physical parameters of sandstone blocks

表2 砂岩块体间接触面力学参数Table 2 Mechanical parameters of contact surfaces between sandstone blocks

图3 砂岩块体单轴压缩试验Fig.3 Uniaxial compression tests of sandstone blocks

在FLAC-3D 中，K，G，E，υ 等参数之间有着如下关系[21]：

其中，Δzmin是指接触面法向连接区域上最小尺寸，kn为法向刚度，ks为剪切刚度.经计算可确定砂岩块体体积模量为3.442 GPa，剪切模量为2.175 GPa，kn，ks皆为6 342 GPa/m.

1.2.5 确定最大静摩擦力

利用FLAC-3D 软件模拟模型仅在水平冲击作用下工作块体的动力响应情况，进行对比分析，可确定出工作块体最大静摩擦力.由图4 可知，当工作块体所受水平静力为71 N 时（图中折线），监测点水平位移Shd产生最大滑移34 µm；而当工作块体所受水平静力为72 N 时（图中阶梯线），监测点水平位移在17.5 ms 处产生最大滑移为82 mm.将两种工况对比可知，72 N 为引起工作块体严重滑移失稳的临界力，即得最大静摩擦力为71 N.

图4 不同水平静力作用下砂岩工作块体水平位移变化曲线Fig.4 Horizontal displacement variation curves of the sandstone working block under different horizontal static forces

1.2.6 数值模拟可行性分析

针对数值模拟结果的可靠性，本文在正式开展超低摩擦鞭梢效应模拟试验之前，设计了一组对比试验说明其可行性.

试验部分通过自主研制的块系模型加载试验装置（如图5 所示）对砂岩块体进行垂向应力波扰动及水平冲击共同作用下工作块体动力响应研究.具体工况如下：1 号块体上表面施加垂直向下应力波扰动（频率3 Hz、振幅200 N），同时在工作块体左侧施加水平向右水平冲击（20 N）.采用数据采集仪（DH-5923N）和三向加速度计（CA-YD-193A）对监测点进行监测，其中监测点为工作块体右侧表面形心处.

图5 块系模型加载试验装置Fig.5 The block model loading test device

模拟部分采用命令流建立对应模型后，将试验部分垂直方向应力波扰动及水平冲击通过编写FISH 函数的方式在FLAC-3D 完成加载模拟，其中垂直方向应力波扰动通过半正弦函数进行模拟.

对比图6（a）和（b）可知，试验所得数据与模拟所得数据整体趋势基本一致，工作块体水平位移皆在剧烈动力响应后逐渐稳定，最终进入平衡状态.工作块体水平位移幅值、波动特性及响应时间仍有些许差异，但对于不同结构特征下工作块体动力响应特征整体趋势及程度的分析影响甚微，仍然具有实际意义.从试验数据与模拟数据对比的角度分析，验证了本文超低摩擦鞭梢效应理论基准模型模拟的可行性及垂直方向扰动、水平冲击等外载施加的准确性.经分析得知，造成试验数据与模拟数据存在一定差距的原因大致有以下3 点：① 针对块体间接触面力学性质的模拟，由于部分参数难以确定导致模拟数据与试验数据只能做到尽可能接近，并不能完全模拟.② 针对试件的模拟，试验所用试件无法做到模拟试件力学性质、体积形态的高度统一，只能通过一些手段降低块体间的差异.③ 针对应力环境的模拟，轴压、应力波扰动的施加是持续的，相对比较稳定；而水平冲击是瞬时的，试验中水平冲击控制存在一个过程，先是计算机下达冲击指令，冲头会在接触块体冲击面后进行瞬间冲击，然后冲头停止运动，此刻岩体极易发生回落碰触冲头后再次弹出，数值模拟中水平冲击则是通过FISH 函数下达指令，表现为瞬间出现、瞬间消失，故在试验中时间响应比数值模拟要慢.

图6 垂直扰动及水平静力作用下工作块体水平位移动态响应曲线：(a) 模拟结果；(b) 实验数据Fig.6 Dynamic response curves of the horizontal displacement of the working block under vertical disturbance and horizontal static force:(a) the simulated results; (b) the test data

2 工作块体尺寸对超低摩擦鞭梢效应影响分析

2.1 数值模拟方案设计

尺寸大小是引起鞭梢效应的重要因素.故本节数值模拟研究计算方案基于超低摩擦鞭梢效应理论基准模型，工作块体始终设为模型中第3 号块体.模型受垂直方向应力波扰动和水平冲击组合作用，其中应力波扰动频率3 Hz，振幅200 N、水平冲击20 N.取边长不同的5 个立方体试件作为工作块体，研究分析尺寸大小对超低摩擦鞭梢效应影响机制，其中工作块体尺寸分别取为标准块体尺寸的5/5，4/5，3/5，2/5，1/5，也即边长100 mm，80 mm，60 mm，40 mm，20 mm 的立方体试件.以工作块体水平位移及加速度作为参考指标，对比分析工作块体不同尺寸作用下工作块体动力响应规律，进而得出模型系统产生超低摩擦鞭梢效应强度随工作块体尺寸变化规律.

2.2 尺寸对工作块体水平位移Shd 影响分析

通过history 命令监测工作块体右侧表面形心处位移及加速度变化情况，再对输出数据加工处理后绘制工作块体水平位移时程曲线，如图7 所示.经分析可得如下结论：

1) 随着工作块体尺寸逐渐变小，工作块体动力响应越趋明显，岩体超低摩擦鞭梢效应强度越大.图7（a）～（e）5 种工况下砂岩工作块体动力响应整体趋势一致，其水平位移时程曲线主要经历两个阶段（动力响应阶段及稳定阶段）.其中，动力响应阶段表现为工作块体监测点水平位移发生突增；稳定阶段则为工作块体动力响应结束后，其水平位移达到幅值并趋于恒定而处于稳定状态.

事实上，图7 中工作块体水平位移响应幅值分别为3.49 mm，4.82 mm，18.12 mm，77.29 mm，321.7 mm，可知随着工作块体尺寸逐渐减小，模型系统产生块体超低摩擦鞭梢效应现象越明显，尤其是在工作块体边长为标准块体的2/5 和1/5 时，模型系统动力响应最为剧烈.

2) 对比已研究的5 种工况发现，边长为40 mm 时，工作块体超低摩擦鞭梢效应显著且数据相对可靠.图7（e）工况在模拟计算过程工作块体水平位移响应最为突出，但其幅值高达321.7 mm，远大于块体自身尺寸，且该工况下动力模拟计算不收敛，数据误差较大.

图7 不同尺寸工作块体水平位移时程曲线：(a) 边长100 mm 立方体试件；(b) 边长80 mm 立方体试件；(c) 边长60 mm立方体试件；(d) 边长40 mm 立方体试件；(e) 边长20 mm 立方体试件Fig.7 Time history curves of horizontal displacements of working blocks with different sizes: (a) the cubic specimen with a side length of 100 mm; (b) the cubic specimen with a side length of 80 mm; (c) the cubic specimen with a side length of 60 mm; (d) the cubic specimen with a side length of 40 mm;(e) the cubic specimen with a side length of 20 mm

2.3 尺寸对工作块体加速度影响分析

同理可得监测点加速度a的变化特征曲线，如图8 所示.其中虚线为工作块体水平方向加速度变化曲线，实线为工作块体垂直方向（z向）加速度变化曲线.通过对比分析工作块体加速度特征曲线可得如下规律：

1) 工作块体尺寸越小，其动力响应则越趋明显，也即超低摩擦鞭梢效应强度越大.纵观图8 中实线（工作块体垂直方向加速度变化曲线）可知，在水平冲击及垂直方向应力波扰动不变的情况下，工作块体尺寸越小会直接引起其质量越小，使得工作块体加速度幅值在逐步增大，其垂直方向加速度幅值依次为2.186 mm/ms2，3.511 mm/ms2，6.756 mm/ms2，20.140 mm/ms2，24.580 mm/ms2，水平方向加速度幅值依次为1.152 mm/ms2，1.618 mm/ms2，3.599 mm/ms2，13.470 mm/ms2，70.030 mm/ms2.

图8 工作块体加速度幅值变化曲线：(a) 边长100 mm 立方体试件；(b) 边长80 mm 立方体试件；(c) 边长60 mm 立方体试件；(d) 边长40 mm 立方体试件；(e) 边长20 mm 立方体试件Fig.8 Change curves of acceleration amplitudes of working blocks: (a) the cubic specimen with a side length of 100 mm; (b) the cubic specimen with a side length of 80 mm; (c) the cubic specimen with a side length of 60 mm; (d) the cubic specimen with a side length of 40 mm; (e) the cubic specimen with a side length of 20 mm

2) 工作块体尺寸为标准块体尺寸2/5 时，极易诱发超低摩擦鞭梢效应.对比图8 中虚线（工作块体水平方向加速度变化曲线）可知，图8（a）、（b）、（c）3 个工况加速度幅值变化甚微，图8（d）、（e）两种工况有着明显动力响应，但图8（e）工况在模拟计算过程中模型系统最大不平衡力与典型内力的比率处于10-1，远大于10-5，计算过程收敛性极差，存在的误也大大增加，反观图8（d）工况数据较好.

3 工作块体位置对超低摩擦鞭梢效应影响分析

3.1 数值模拟方案设计

工作块体位置是影响超低摩擦鞭梢效应的重要因素.由图7、8 分析可知，工作块体边长为标准块体边长2/5 时（也即边长为40 mm），工作块体超低摩擦鞭梢效应尤其显著，数据相对较好.故在设计工作块体位置对超低摩擦鞭梢效应影响模拟试验时，采用工作块体边长为40 mm 的立方体试件作为基准.在该部分研究中设有4 种工况，即工作块体分别为模型中2 号、3 号、4 号及5 号块体，模型受垂直方向应力波扰动和水平冲击组合作用，其中应力波扰动频率3 Hz、振幅200 N、水平冲击20 N，水平冲击作用于工作块体左侧表面，监测点设于工作块体右侧表面形心，水平冲击及监测点随工作块体位置变化而变化.工况具体指代模型如图9 所示.

3.2 位置对工作块体水平位移影响分析

将边长为40 mm 的立方体试件分别置于不同位置，构成工况1～4 具体模型图及相应工况模型下工作块体水平位移响应曲线，如图9（a）～（d）所示.分析可得：

1) 工作块体水平位移幅值随工作块体与扰动源之间距离的增大呈先增后减的趋势.纵观图9（a）～（d），4 种工况在图中都呈现出先增加后稳定的趋势，且4 种工况中工作块体对应水平位移变化曲线整体趋势一致，先后都经历了动力响应阶段及稳定阶段.4 种工况下工作块体水平位移幅值分别为9.79 mm，17.78 mm，21.03 mm，16.04 mm，说明工作块体与扰动源之间距离在一定范围内越大，其水平位移幅值越大，超低摩擦鞭梢效应强度越大；超出临界范围，工作块体水平位移则会随着距离的持续增大而呈递减趋势，超低摩擦鞭梢效应减弱.同时说明，距离扰动源最近的岩体一般不是动力响应最强烈的那部分.

图9 工作块体不同位置下水平位移时程曲线：(a) 工况1；(b) 工况2；(c) 工况3；(d) 工况4Fig.9 Time history curves of horizontal displacements at different heights of the working block: (a) working condition 1;(b) working condition 2; (c) working condition 3; (d) working condition 4

2) 位于4 号的工作块体水平位移响应幅值最大，其中工作块体形心与扰动源之间距离占整体模型高度的16/21.将工况1～4 在内的4 种工况下工作块体水平位移幅值Shda绘制点线图，如图10 所示，可直观地得知工作块体水平位移随着工作块体位置的远离呈递增趋势.但由能量耗散及扰动的反射可知：在工作块体远离扰动源的过程中存在一个临界距离使得工作块体超低摩擦鞭梢效应强度最大，超过这一临界值后，工作块体水平位移会随着工作块体与扰动源距离的持续增加而表现为递减趋势，本文模拟计算中临界距离趋于0.32 m.从整体角度看，工作块体水平位移幅值Shda随着距离的增加呈现出先增后减的趋势，类似于上凸型抛物线.

图10 不同工况下工作块体水平位移幅值变化曲线Fig.10 Variation curve of horizontal displacement amplitude of working block under different working conditions

3.3 位置对工作块体加速度影响分析

对工作块体受力分析如图11 所示，其中G0为单个标准块体重力位置，n为工作块体上方标准块体个数，Gs为工作块体自身重力，f1，f2分别为工作块体上下接触界面摩擦力，N为支持力.

由图12 可知，工况1～4 垂直方向加速度幅值分别为13.290 mm/ms2，13.470 mm/ms2，11.330 mm/ms2，2.363 mm/ms2，水平方向加速度幅值分别为3.182 mm/ms2，3.241 mm/ms2，5.626 mm/ms2，3.837 mm/ms2，绘图可得图13.基于图11 中工作块体受力情况分析及图13 可得如下规律：

图11 工作块体受力示意图Fig.11 The force diagram for the working block

图12 工作块体加速度幅值变化曲线：(a) 工况1；(b) 工况2；(c) 工况3；(d) 工况4Fig.12 Change curves of acceleration amplitudes of the working block: (a) working condition 1;(b) working condition 2; (c) working condition 3; (d) working condition 4

图13 工作块体加速度幅值变化曲线Fig.13 Change curves of acceleration amplitudes of the working block

1) 随着工作块体与扰动源距离的增加，其z向加速度幅值逐渐降低，x向加速度幅值先升后降.即工作块体在远离扰动源过程中，其z向加速度幅值并未随着上方标准块体对其产生的压力逐渐增大而增大，而是呈现出先增后减的趋势，说明这种鞭梢结构特征影响着工作块体z向受力.

2) 工作块体水平加速度幅值随其与扰动源的距离增大呈先增后减的趋势，工况3 为试验过程中的一个临界工况，动力响应较为显著.在工况3 的模拟模型下，工作块体所产生的超低摩擦鞭梢效应强度最大.

4 结 论

1）系统产生超低摩擦鞭梢效应难易程度与工作块体尺寸密切相关.工作块体尺寸越小，其动力响应则越趋明显，也即超低摩擦鞭梢效应强度越大.就本文所建模型而言，工作块体尺寸为标准块体尺寸2/5 时极易诱发超低摩擦鞭梢效应.

2）系统产生超低摩擦鞭梢效应难易程度与工作块体所处位置密切相关.超低摩擦鞭梢效应强度随工作块体与扰动源之间距离增大呈先增后减趋势，也即上凸型抛物线关系，如图10 所示.临界值内，距离越大则能引起超低摩擦鞭梢效应强度越大，当超过这一临界值则会随之减小，文中临界工况为工作块体形心与扰动源之间距离占整体模型高度的16/21.

3）块体尺寸及块体位置是岩体超低摩擦鞭梢效应的两个重要因素，深部开采时应多注意岩体整体结构特征.