地下巷道开挖过程往往造成围岩的变形和应力重分布, 对围岩应力的准确分析是进行巷道稳定性控制的前提.岩体中往往存在大量的不连续面, 如裂隙、节理、断层及层理面等.在地下巷道开挖过程中, 不连续面对围岩变形及稳定性有显著的影响^[1].目前, 评价节理对巷道围岩变形及稳定性的影响主要采用岩体分级法、不连续分析法及等效连续分析法.其中, 不连续分析法可对不同节理面对变形影响的方向效应进行有效的模拟, 且能模拟不连续岩体的大变形和转动问题, 也能进行等效连续分析.

目前, 已有多种数值模拟方法可用于分析和评价节理岩体中开挖工程的变形及稳定性评价.有限元法、边界元法及拉格朗日法在某些情况下可以通过接触单元模拟不连续面.但是, 上述方法不能自动识别计算过程中产生的新的不连续面, 只能模拟小变形或转动.对于地下开挖的大变形问题, 离散元法更加适用.目前, 由Cundall提出的离散单元法和由Lemos等改进的离散单元法是用来分析不连续块体应力变化的常用工具^[2-3].考虑到块体和节理间的相互作用, 离散元方法可以有效分析不同应力及应变条件下块体的力学行为, 可以较容易地计算出大变形和完整岩体及节理间的复杂本构行为.

离散单元法被用来评价地下硐室群的稳定性^[2-9], 绝大多数研究者都是使用UDEC (universal distinct element code)在二维状态下通过对具有不同尺寸的岩块进行力学试验, 在莫尔库伦理论的基础上求得岩体的力学参数, 提出适用于节理岩体的介质模型^[10-12], 并对节理岩体力学性质进行分析, 得到对应的等效力学参数^[13-17].

本文采用3DEC(three dimension distinct element code)软件建立三维离散模型, 研究节理面对巷道开挖围岩稳定性的影响、变形规律, 分析岩块和节理力学参数对巷道围岩破坏的影响和等效连续性特征.

1 工程背景及基础参数

程潮铁矿是我国著名的大型地下金属矿山, 矿体呈南东65°走向, 倾角约为30°, 区内矿体宽度大, 最宽处超过100 m, 是典型的缓倾斜厚大矿体.程潮铁矿以大理岩为主体, 是由闪长岩、铁矿体接触破碎带、大理岩和花岗岩共同组成的封闭性的裂隙岩溶综合性含水构造带.

目前, 影响程潮铁矿围岩稳定性的因素主要包括：1)矿体内沿着水平最大应力分布大量的节理；2)矿体内高地应力的存在；3)每种岩性内分布2~3组节理组, 且方向各异；4)节理密度大, 分布在4~22条/m；5)岩块及节理力学参数的各向异性及不确定性.由上述分析可知, 程潮铁矿围岩变形及稳定性是由节理控制, 因此, 有必要开展节理性质对围岩变形的影响规律研究.

1.1 岩石及节理力学参数

选取磁铁矿、花岗岩、矽卡岩、大理岩、闪长岩和硬石膏6种岩性129件试件, 进行岩石密度和含水率的测定、单轴试验和三轴试验.如表 1所示为所选用的各个岩性的物理力学参数, 其中ρ为密度, K为体积模量, G为剪切模量, C为内聚力, Φ为内摩擦角, σ_t为抗拉强度.节理的物理力学参数如表 2所示.表中, J_i为节理倾向, J_Φ为节理倾角, J_s为节理间距, r为节理迹长.不同岩性的节理面的力学属性如表 3所示.

由于2种相邻的岩性之间存在层面, 在模拟时应考虑层面的影响.2个相邻材料之间层面的力学参数如表 4所示, k_N为法向刚度, k_S为剪切刚度.

1.2 地应力场分布规律

如表 5所示为现场监测的应力测量结果.表中, h为开采深度, σ₁为最大水平主应力, σ₂为垂直第二主应力，γH为垂直方向上的应力, σ₃为最小水平主应力.根据测试结果进行计算可得, -395 m开采水平的最大主应力σ₁=1.60 γH，最小主应力为0.39 γH, 施加在模型上方的应力为13.96 MPa (3DEC中的z方向应力σ_zz), 最大水平主应力为22.34 MPa(3DEC中的x方向应力σ_xx), 同样, 最小水平主应力5.44 MPa (3DEC中的y方向应力σ_yy), 最大水平主应力方向为NW75°.在3DEC中, 最大水平主应力和最小水平主应力方向分别为x和y方向.

目前, 程潮铁矿采用锚喷支护技术.顶板的砂浆锚杆长度为2.0 m, 两帮为1.8 m, 底板未支护, 锚杆间距为1.0 m.钢筋网度为0.2 m, 钢筋直径为6 mm, 喷层厚度为50 mm.巷道的支护断面图如图 1所示.模型中锚杆的杨氏模量为210 GPa, 杆体横截面积为2.55×10^-4 m², 抗拉能力为240 MN, 单位长度的灌浆黏结力为100 MN/m, 单位长度的灌浆刚度为2×10³ MN/m².

2 模型的建立与参数校核 2.1 模型构建

基于3DEC建立数值计算模型, 如图 2所示, 其中位于模型中心的是-395 m水平巷道, 模型的尺寸为40 m×10 m×40 m, 模型中包含的岩性有闪长岩、大理岩、花岗岩、磁铁矿.

2.2 模型试算

为了系统研究力学参数对巷道位移的影响, 需要进行多种不同方案的研究, 根据程潮铁矿现场地应力环境、各个岩性的力学参数、节理参数等实测数值, 提出试算方案：

1) 方案1：模型中无节理, 所选用的岩石力学参数为程潮铁矿岩样的实测值, 在此将岩样的实测值定义为岩石力学参数的标准值.

2) 方案2：模型中有节理, 所选用的岩石力学参数和节理参数为程潮铁矿岩样的实测值, 在此将岩样和节理的力学参数实测值定义为标准值.

3) 方案3：模型中有节理、有支护, 岩石和节理的力学参数均为标准值.

根据方案1的计算结果, 分析巷道开挖后模型的最大水平应力和垂直应力分布情况, 以及模型主应力矢量的分布情况.文中所有的结果图都取自垂直于y=0的面.如图 3所示为最大水平主应力和垂直应力的分布情况, 如图 4所示为巷道开挖后达到平衡后主应力矢量的分布情况.

图 3(a)为巷道附近的最大水平主应力分布情况, 在巷道的两帮出现0值, 最大水平主应力分布情况符合地应力的分布规律, 最高的最大水平主应力出现在巷道顶板和底板处.图 3(b)为巷道周围垂直应力分布情况, 0值出现在巷道顶板和底板处, 最大的垂直应力出现在巷道两帮, 模型从上到下垂直应力的分布符合所施加的地应力.

由计算结果可知, 应力分布云图符合预期, 效果良好, 表明了用来模拟巷道开挖的3DEC源文件的正确性.

2.2 与现场实测值的对比分析

对程潮铁矿巷道围岩表面和深部的位移变化进行监测, 通过对数值模拟和现场监测结果的对比分析, 验证模型建立的正确性.多点位移计采用自设式和JSS30A型10 m电子数显收敛计, 对-395 m水平进行长期监测, 监测点高度设置在距离巷道底板1.2 m左右, 监测点布置如图 5所示.

由于现场岩体存在大量的节理, 且存在支护, 因此, 在与现场监测结果进行对比分析时, 选取方案3进行模拟.巷道监测所得到的结果如图 6和7所示, 图中, S为位移, T为观测天数.由图可知, 巷道变形的最大位移出现在巷道两帮的边缘外侧, 越往内侧, 变形值越小.如图 8所示为现场监测结果和不同数值模拟结果的对比图, 图中H为监测点深度.

总体来说, 现场监测的结果和数值模拟的结果匹配较好, 右帮位移的数值模拟结果平均值约为现场监测的86%左右, 左帮位移数值模拟结果平均值约为现场监测的79%左右, 且监测深度越大, 数值模拟结果与现场监测结果越接近, 表明此次模拟选取的力学参数与建立的模型都较为合理.

3 节理力学参数敏感性分析 3.1 单因素影响程度分析

影响巷道围岩稳定性的因素较多, 为了找出各个因素的重要程度, 提出不同的计算方案, 对比分析巷道位移值, 确定因素的影响程度排序, 具体计算方案如下：

1) 方案4：模型中无节理、无支护, 在一定范围内改变内聚力的数值, 其他不变；

2) 方案5：模型中无节理、无支护, 在一定范围内改变抗拉强度的数值, 其他不变；

3) 方案6：模型中无节理、无支护, 在一定范围内改变内摩擦角的数值, 其他不变；

4) 方案7：模型中有节理、无支护, 在一定范围内改变JKN的值, 其他不变；

5) 方案8：模型中有节理、无支护, 在一定范围内改变JKS的值, 其他不变；

6) 方案9：模型中有节理、无支护, 在一定范围内改变节理的内摩擦角, 其他不变.

计算时所采用的力学参数标准值和因素的改变范围见表 6所示, 表中Φ_j为节理内摩擦角.

如图 9和10所示, 为改变不同力学参数对巷道开挖附近的最大位移的影响.图中, S_lmax为左侧最大位移, S_rmax为右侧最大位移, S_tmax为顶板最大位移, S_bmax为底板最大位移.

由图可知, 力学参数的增大导致最大位移的变小.由于不同类型的力学参数具有不同的量纲, 不同参数之间没有共性, 很难直接进行比较, 因此将这些指标按照某种标准归一化到无量纲区间.本文按照如下的公式将各参数无量纲化：

$L = \frac{{\left( {{S_{rp}}/\Delta {S_{rp}}} \right)}}{{\left( {{\Delta _p}/rp} \right)}} \times 100\% $

(1)

式中：L为因素的敏感性；S_rp为基准参数下巷道位移对值；ΔS_rp为由参数改变引起的巷道位移变化；Δ_p为力学参数的变化量；rp为基准值.

在计算过程中, 对于每种参数, 选取其左帮最大位移进行比较, 由计算结果可知, 改变C、σ_t、Φ、JKN、JKS、Φ_j, 位移该变量分别为70.5%、4.7%、141.5%、62.1%、66.8%、52.1%.由此可知, 各个岩体力学参数对巷道位移变化的影响程度由大到小依次是：Φ＞C＞JKS＞JKN＞Φ_j＞σ_t.

3.2 多因素影响程度分析

为了研究模型的极限值, 所有完整岩石和节理力学参数都在一定范围内同时增大或者减小, 根据相应最大位移变化的数值确定参数极限值, 参数增大或者减小后模型分别为极限塑性模型和极限弹性模型.完整岩石的力学参数分别增大和减小了20%, 节理和不连续面的力学参数分别增大和减小了30%, 力学参数增大的模型, 为极限塑性模型, 力学参数减小的模型, 为极限弹性模型.设计了如下计算方案：

方案10：模型无节理、无支护, 模型选用的岩石力学参数为标准值的80%；

方案11：模型无节理、无支护, 模型选用的岩石力学参数为标准值的120%；

方案12：模型有节理、无支护, 模型选用的岩石力学参数为标准值的80%, 节理力学参数为标准值的70%；

方案13：模型有节理、无支护, 模型选用的岩石力学参数为标准值的120%, 节理力学参数为标准值的130%；

方案14：模型有节理、有支护, 模型选用的岩石力学为标准值的80%, 节理力学参数为标准值的70%；

方案15：模型有节理、有支护, 模型选用的岩石力学参数为标准值的120%, 节理力学参数为标准值的130%.

通过计算可以得到各个方案中左帮、右帮、顶板和底板的最大位移, 将这些位移值汇总到图 11.由图 11可知, 当模型中没有节理时, 改变完整岩石的力学参数, 塑性模型的巷道最大位移是弹性模型的2.5~2.8倍.当模型中有节理存在时, 改变各种岩性和节理的力学参数, 塑性模型的巷道最大位移是弹性模型的1.8~4.4倍.

4 岩体等效连续性特征

节理的存在使得岩体的力学特征表现出尺寸效应和各向异性特征, 降低了整个岩体的强度.节理对整个岩体完整性和岩石力学参数的影响取决于裂隙网络中完整岩石和不连续面的属性值之间的差值^[18].为了研究岩体的等效连续性力学行为, 需要为模型选取合适的岩石力学参数, 以便等效模拟模型中含有节理时产生的力学行为.根据裂隙张量理论及裂隙张量和张量的第一不变量与岩石变形模量、剪切模量、内聚力和内摩擦角之间的关系, 首先计算出模型中各岩性的裂隙张量和第一不变量, 然后根据关系曲线计算出等效的岩石变形模量、剪切模量、内聚力和内摩擦角等参数值.假设节理是一个薄圆盘, 三维状态下第k组的节理裂隙张量的表达形式^[19-20]为：

$F_{ij}^{\left( k \right)} = 2{\rm{ \mathit{ π} }}\rho \int_0^\infty {\iint_{\mathit{\Omega }/2} {{r^3}{\mathit{\boldsymbol{n}}_i}{\mathit{\boldsymbol{n}}_j}f\left( {\mathit{\boldsymbol{n}},r} \right){\rm{d}}\mathit{\Omega }{\rm{d}}r}} .$

(2)

式中：ρ为节理的密度；n为直于节理面的单位矢量；f(n, r)为节理组n和r的概率密度函数；Ω/2为与节理面的单位半球相关的立体角；n_i和n_j(i, j=x, y, z)为直角坐标系中的矢量分量.

如果节理的方向和迹长是无关联的, 则三维状态下的裂隙张量公式可以表述为

$F_{ij}^{\left( k \right)} = 2{\rm{ \mathit{ π} }}\rho \int_0^\infty {{r^3}f\left( r \right){\rm{d}}r} \iint_{\mathit{\Omega }/2} {{\mathit{\boldsymbol{n}}_i}{\mathit{\boldsymbol{n}}_j}f\left( \mathit{\boldsymbol{n}} \right){\rm{d}}\mathit{\Omega }}.$

(3)

式中：f(n)为法向矢量n的概率密度函数；f(r)为迹长r的概率密度函数.

因此, 第k组节理的定向裂隙张量分量为

$F_{l}^{\left( k \right)} = 2{\rm{ \mathit{ π} }}\rho \int_0^\infty {{r^3}f\left( r \right){\rm{d}}r} \iint_{\mathit{\Omega }/2} {{\mathit{\boldsymbol{n}}_i}{\mathit{\boldsymbol{n}}_j}f\left( \mathit{\boldsymbol{n}} \right){\rm{d}}\mathit{\Omega }}.$

(4)

式中：l为方向矢量.

如果岩体中包含不止一组节理, 则整个岩体的裂隙张量可以表示为

${F_{ij}} = \sum\limits_{k = 1}^N {F_{ij}^{\left( k \right)}} .$

(5)

式中：N为岩体中的节理组数.

若F_ij=F_ji时, 裂隙张量F有3个主应力分量：F₁、F₂和F₃, 对应的裂隙张量的不变量为

$\left. {\begin{array}{*{20}{l}} {{I_{1\left( F \right)}} = {F_1} + {F_2} + {F_3},} \\ {{I_{2\left( F \right)}} = - \left( {{F_1}{F_2} + {F_2}{F_3} + {F_3}{F_1}} \right),} \\ {{I_{3\left( F \right)}} = {F_1}{F_2}{F_3}.} \end{array}} \right\}$

(6)

根据式(6)计算出各节理组的裂隙张量的分量和不变量, 再根据裂隙张量理论以及裂隙张量与张量的第一不变量与岩石变形模量、剪切模量、内聚力和内摩擦角之间的关系, 计算得出等效的岩石变形模量、剪切模量、内聚力和内摩擦角等参数值.

根据裂隙张量理论, 可计算出等效连续模型的力学参数, 为非等效连续模型中的力学参数的40%, 设非等效连续模型中的力学参数为标准值.采用3DEC软件研究不同岩体力学参数与标准值的比例, 具体方案如下：

方案16：无支护, 岩石力学参数为标准值的70%；

方案17, 无支护, 岩石力学参数为标准值的60%；

方案18：无支护, 岩石力学参数为标准值的50%；

方案19：无支护, 岩石力学参数为标准值的40%；

方案20, 无支护, 岩石力学参数为标准值的30%；

方案21：有支护, 岩石力学参数为标准值的80%；

方案22：有支护, 岩石力学参数为标准值的70%；

方案23：有支护, 岩石力学参数为标准值的60%；

方案24：有支护, 岩石力学参数为标准值的50%；

方案25：支护, 岩石力学参数为标准值的40%；

方案26：有支护, 岩石力学参数为标准值的30%.

如图 12所示为计算所得到的结果, 图中Δ_i为岩石的力学参数相对于标准值的比例.

如表 7所示为不同方案的巷道最大位移的对比结果.对于有支护和没有支护2种情况下的模型, 当等效连续模型使用的参数为节理岩体参数的40%~50%时, 等效连续模型和非等效连续模型会达到比较相似的结果.

由图 12和13的对比可知, 支护有效地抑制了顶板的下沉, 底板的位移则减小幅度较小.支护有效地控制了巷道的变形, 可增强巷道的稳定性.

5 结论

(1) 通过现场监测和数值模拟结果对比可知, 数值模拟结果中巷道左、右两帮与监测结果且监测深度越深, 数值模拟结果与现场监测结果接近, 表明此次模拟选取的力学参数与建立的模型都较为合理.

(2) 各参数对巷道位移的影响程度, 从大到小依次如下：岩石的内摩擦角, 岩石的内聚力, 节理的剪切刚度, 节理的法向刚度, 节理的内摩擦角, 岩石的抗拉强度.

(3) 当模型中未添加节理, 仅仅改变岩石力学参数时, 塑性模型的巷道最大位移是弹性模型的2.5~2.8倍.对于模型中有节理时, 同时改变各种岩性和节理的力学参数, 塑性模型的巷道最大位移是弹性模型的1.8~4.4倍.

(4) 当等效连续模型使用的参数为节理岩体参数的40%~50%时, 等效连续模型和非等效连续模型会达到比较相似的结果.

(5) 增加支护使变形减小了约43%~51%, 顶板位移减小了约29%~39%, 底板则只减小了9%~10%, 主要原因是只在两帮及顶板施加了支护, 而底板没有支护.由此可见, 支护提高了巷道的稳定性.