考虑结构面黏聚力损伤变化的洞室块体稳定分析

骆 健， 章杨松，*， 李晓昭

(1. 南京理工大学理学院， 江苏 南京 210094; 2. 南京大学地球科学与工程学院， 江苏 南京 210093)

0 引言

岩体在形成和演化的过程中，由于受到各种地质作用和构造力的影响，会形成大小不一、形态复杂、分布错综复杂的岩体结构面，并与隧道开挖形成的临空面组成不同规模的岩石块体，这些块体的失稳或垮落会破坏岩体整体的稳定性。因此,在设计阶段，根据形成的结构面信息分析块体的稳定性，对优化设计和开挖支护方案具有重要的指导意义。

由石根华[1]提出的关键块体理论和P.M.Warburton[2]提出的块体可移动性的矢量法是岩石块体稳定性分析方法的重大突破，为用数值方法研究裂隙岩体稳定奠定了重要基础。研究裂隙岩体稳定的方法主要有连续介质法、非连续介质法和连续-非连续介质的混合法等[3]。黄达等[4]运用数值计算软件FLAC3D，通过连续介质法考虑块体的变形以及塑性区分布特征来评价块体的稳定性。张雨霆等[5-6]提出复杂岩石块体识别的单元重构-聚合法，并进一步采用连续介质法，运用数值分析软件FLAC3D中的Interface界面单元模拟岩石块体稳定性。Wu等[7]提供了一种DDA(非连续变形分析)的离散化方法来计算块体稳定性，此方法能够模拟出块体的掉落行为。

连续介质法的缺点是不能模拟块体的掉落行为；在模拟具有胶结结构面的块体失稳时，目前常用的非连续介质法无法分析其胶结面的渐进损伤过程及状况。针对以上问题，本文采用显式动力有限元模拟裂隙岩体的连续-非连续变形以及含胶结的结构面组成块体的局部失稳过程。选取甘肃北山某典型花岗岩研究区域测量的结构面数据以及搜索出的多种块体，运用有限元软件ABAQUS中的黏聚力模型(cohesive zone model，CZM)模拟具有填充和胶结的结构面的抗拉强度和黏聚强度。根据测量的结构面数据信息，分4种方案对胶结结构面的抗拉强度和黏聚强度进行赋值；运用Matlab程序对隧道洞室模型的inp文件进行修改，从而可以对洞室开挖边界施加相应的节点荷载，通过对反加于围岩边界荷载的逐渐释放来分析初始地应力和开挖诱导的二次地应力对块体稳定的影响。结合以上几点，综合考虑初始地应力、开挖诱导的二次地应力、结构面的黏聚强度和抗拉强度，以及它们之间的摩擦应力等因素来研究块体的稳定性。

1 黏聚力模型

(a) 法向

(b) 切向图1 双线性张力-位移关系Fig. 1 Bi-linear tension-displacement relationship

1.1 黏聚力单元损伤原理

在黏聚力单元损伤之前，认为其应力应变关系符合线弹性准则，如式(1)所示。2个相邻单元之间用类似胶材料的黏聚力单元连接，该黏聚力单元很薄，有时考虑到实际的工程问题，黏聚力单元厚度常常设置成0。

(1)

式中：t为单元应变矢量；tn、ts、tt分别为单元法向和2个切向的牵引力；K为单元刚度矩阵；ε为应变矢量；εn、εs、εt分别为单元法向和2个切向的应变。

矢量ε表达式如下：

εn=dn/To,εs=ds/To,εt=dt/To

。

(2)

式中：dn、ds、dt分别为单元法向和2个切向变形黏聚力单元分开的距离；To为黏聚力单元最初的厚度。

利用名义应力比的二次交互函数准则来确定胶结结构面是否会断裂，当临界应力比平方和等于1时，认为损伤完全。此准则用公式表达如下[9]：

(3)

1.2 损伤模式

研究表明，裂纹的扩展是一个复合型损伤演化过程，包括法向和2个切向的扩展。B-K准则[10]被广泛应用于预测裂缝扩展，其表达式见式(4)。该式用于当岩体完全沿第1、2剪切方向变形时，2个方向的临界断裂能相同的情况。

(4)

1.3 损伤演化

模拟材料的断裂损伤演化通常采用刚度退化准则，该准则表达式如下[11]：

(5)

(6)

零厚度的黏聚力单元模型如图2所示。图中的中间层表示具有零厚度的黏聚力单元层，当块体受到外力作用时，中间的黏聚力层会逐渐损伤；当其应力或应变满足定义的初始临界损伤准则时，则开始退化，黏聚力单元逐渐删除，单元间的黏聚力消失，而上下面之间的摩擦力从分析步开始只要面面接触就会一直起作用。根据摩尔-库仑理论，黏聚力单元层的强度τ=c+σtanφ，其中c为结构面的黏聚力强度，σtanφ为结构面之间的摩擦强度。

图2 零厚度的黏聚力单元模型Fig. 2 Cohesive element model of zero thickness

2 围岩块体的稳定性分析

在甘肃北山某花岗岩典型区域研究范围内进行块体的搜索，得到的块体如图3所示。其中，1、2、4号块体是四面体，3号块体是五面体，且由几何位置初步判断，1、2、3号块体可能出现失稳情况，4号块体是稳定块体。

1—4表示块体编号。

图3搜索出的块体示意图
Fig. 3 Sketch of searched blocks

2.1 模型建立

根据得到的块体结构面数据，以花岗岩典型区域地下开挖隧道为例，对开挖范围内搜索出的块体的稳定性进行研究。开挖隧道断面由边墙与拱顶组成，可以分别定义为直线段和圆弧段。考虑到模拟隧道开挖仅对一定范围内的岩体产生影响，根据工程经验，模型边界距开挖洞口边的距离超过5倍洞口直径时，开挖产生的影响就会变得很小。在本次模拟中，考虑到洞口直径为5 m左右，模型尺寸可取为40 m×60 m×60 m (长×宽×高)。根据北山岩体条件及核废料地下处置的基本要求，计算模型上表面位于地下600 m埋深处，同时上表面的竖向应力按实测应力取值； 下表面位于660 m埋深处。模型的边界条件为： 上边界为自由边界，下边界约束竖向(z向)位移，左右边界为侧向约束。单元网格的尺寸从模型边界到孔洞壁附近逐渐被细化，如图4所示。

图4 隧道模型的网格划分Fig. 4 Meshing of tunnel model

2.2 材料及结构面参数

模型顶部为应力边界条件，应力值为上覆岩石的质量，依据赵宏刚等[12]得出的岩石埋深与地应力的关系，取垂直主应力为16.2 MPa，x方向的侧压力系数为0.92，y方向的侧压力系数为0.65。根据北山某典型花岗岩采石场地区的统计资料，分别取表1所示的岩体物理力学参数，按照等效连续介质进行计算。

表1 岩体的物理力学参数Table 1 Physico-mechanical parameters of rock mass

结构面的强度包含抗剪强度和抗拉强度2部分，结合搜索出块体的结构面特征以及相应的资料[13]，参照GB/T 50218—2014《工程岩体分级标准》[14]，得出其抗拉强度在0～0.04 MPa，相应的黏聚强度在0～0.12 MPa，内摩擦角在19°～29°。

2.3 模型分析步

本次模拟分为2个分析步，每个分析步与时间步相对应。第1个分析步是初始地应力平衡，时间步是0～0.1； 第2个分析步是模拟隧道开挖，时间步是0.1～1.4。为了分析开挖诱导的二次地应力的影响，把初始平衡步得到的开挖边界应力加到隧道边界上，如图5所示。随着时间步的增加，开挖边界的节点力逐渐释放，当时间步增加到1时，节点力释放完成。释放过程中，由于中间岩体被挖出，原有的地应力平衡被打破，将由此产生应力二次重分布； 同时，原有的结构面与新产生的临空面形成关键块体，在地应力以及自身重力的作用下，此块体可能会出现失稳。

图5 隧道边界节点力Fig. 5 Joint force of tunnel boundary

2.4 模拟结果

考虑到不同位置的块体可能的失稳情况，分4种方案进行模拟,具体参数见表2。整个方案中，4号块体始终保持稳定。为了仅考虑摩擦力的影响，方案1把抗拉强度和黏聚强度赋予一个很小的数值(1×10-9MPa)，相应的摩擦因数见表2； 方案2、3、4分别考虑不同位置块体可能的失稳情况，通过调整不同块体结构面的抗拉强度和黏聚强度，确定影响块体稳定的主要因素。

表2 块体结构面参数的4种方案Table 2 Four schemes of block structural surface parameters

1)方案1模型不同时间步的位移云图如图6所示。在抗拉强度和黏聚强度接近于0，仅摩擦应力起作用的情况下，1、2、3号块体几乎同时掉落。图7示出开挖时间步为0.1～1.4时的摩擦应力变化情况。由图7可知，摩擦应力在起始阶段逐渐增大，达到峰值后逐渐降低并趋于0。对应图6来看，趋于0即块体已经完全脱离岩体结构面并向下滑落，顶部3号块体大约在时间步为0.5时最先脱离，另外2个块体大约在时间步为0.6时脱离。2号块体的摩擦应力始终相对较小，因此若想让2号块体保持稳定，提高块体结构面的抗拉强度是最有效的办法。对于3号块体，由于其位置的特殊性而最容易失稳，在接下来的探讨中，发现即使赋予该结构面参数区间内最大的抗拉强度、黏聚强度和摩擦因数，该块体还是无法保持稳定。

2)方案2模型的位移云图如图8所示。由图8可知： 开挖时间步为0.1时，3号块体已经失稳，而2号块体此时还处于稳定状态； 随着开挖进行，二次应力逐渐释放，在时间步为1时，应力释放完成，2号块体已经失稳，1号块体始终保持稳定状态，块体的位移曲线如图9所示。下面主要探讨3号块体结构面的损伤情况，黏聚力单元的损伤情况如图10所示。图10(a)示出时间步为0.1时，3号块体的结构面接近损伤完全(时间步为1时损伤破坏)，黏聚力单元的底边边角损伤最为严重； 时间步为0.3时，底边边角和两侧面部分已经完全损伤，并失效删除； 时间步为0.4时，已经完全失效删除。

(a) 时间步为0.3

(b) 时间步为1图6 方案1不同时间步的位移云图(单位： m)Fig. 6 Displacement nephograms of block with different time steps in scheme 1 (unit: m)

图7 方案1不同块体的摩擦应力的时间步曲线Fig. 7 Curves of relationship between time steps and friction stresses of different blocks in scheme 1

(a) 时间步为0.1

(b) 时间步为1图8 方案2不同时间步的位移云图(单位： m)Fig. 8 Displacement nephograms of block with different time steps in scheme 2 (unit: m)

图9 方案2块体位移随时间步的变化情况Fig. 9 Curves of relationship between time steps and friction stresses of different blocks in scheme 2

3)方案3模型的位移云图如图11所示。由图11可知: 降低1号块体结构面的摩擦因数而黏聚力不变，块体在时间步为1时出现失稳并逐渐滑落，此时由于2号块体结构面的黏聚力增加到0.1 MPa，摩擦因数为0.45，整个分析步过程中，2号块体都保持稳定； 3号块体此时已被赋予结构面参数区间的最大抗拉强度、黏聚强度和摩擦因数，但依然失稳。方案3块体位移曲线如图12所示，1号块体结构面的损伤情况如图13所示。

(a) 时间步为0.1

(b) 时间步为0.3图10 3号块体结构面随时间步变化的损伤情况Fig. 10 Damage of structural surfaces of block 3 variating with time steps

(a) 时间步为0.3

(b) 时间步为1图11 方案3不同时间步的位移云图(单位： m)Fig. 11 Displacement nephograms of block with different time steps in scheme 3 (unit: m)

4)方案4位移云图如14所示。把2号块体结构面的摩擦因数降到最小值0.34，此时黏聚力没有变，2号块体在时间步为0.8时出现失稳，到时间步为1时出现明显的滑落情况。值得注意的是，计算过程中1号块体结构面的摩擦因数由方案3的0.4提高到0.42，而此时的1号块体却比方案3滑落的早。由此可见，块体间的滑落是相互影响的，一个块体的滑落势必会影响其他块体。方案4块体的位移曲线如图15所示。

图12 方案3块体位移随时间步的变化情况Fig. 12 Curves of relationship between time steps and friction stresses of different blocks in scheme 3

(a) 时间步为0.3

(b) 时间步为0.9图13 1号块体结构面随时间步变化的损伤情况Fig. 13 Damage of structural surfaces of block 1 variating with time steps

(a) 时间步为0.5

(b) 时间步为1图14 方案4不同时间步的位移云图(单位： m)Fig. 14 Displacement nephograms of block with different time steps in scheme 4 (unit: m)

图15 方案4块体位移随时间步的变化情况Fig. 15 Curves of relationship between time steps and friction stresses of different blocks in scheme 4

3 结论与讨论

本文运用黏聚力模型，结合甘肃北山某花岗岩典型区域搜索出的块体，通过考虑具有填充和胶结的结构面的抗拉强度、黏聚强度以及摩擦应力等因素，分析不同位置块体的稳定情况。得到主要结论如下。

1)采用黏聚力单元不仅能模拟仅有摩擦应力作用的结构面情况，还能较好地模拟具有填充和胶结的结构面的抗拉强度和黏聚强度，以及结构面的滑移和脱开行为，可作为块体稳定分析计算的工具。

2)实际隧道开挖过程中，顶部的块体最容易失稳，因此在为隧道设计衬砌时应适当加强衬砌顶部的强度。此外，抗拉强度、黏聚强度和摩擦应力共同构成了结构面的强度，不同位置的块体滑落是相互影响的，一个块体的滑落势必会影响其他块体的稳定性。

3)采用黏聚力单元模拟胶结结构面的抗拉强度和黏聚强度，通过把初始平衡步得到的开挖边界应力反加到隧道边界上，并在开挖步逐渐释放边界应力，从而可以表征出黏聚力单元随二次地应力释放而逐渐损伤的过程，这为模拟结构面损伤随时间的渐进性发展提供了一个新思路。

本文的研究尚有许多不足之处，使用的黏聚力损伤模型为常用的双线性函数，与实际胶结结构面随外力具体损伤的契合度也有待改进。此外，模拟含有岩桥的块体稳定也是以后研究的重点与难点。