基于节理网络有限元的岩质边坡稳定性分析

，，，，

(1.成都理工大学 地质灾害防治与地质环境保护国家重点实验室，成都 610059；2.广西交通设计集团有限公司，南宁 530011； 3.中国电建集团 西北勘测设计研究院有限公司，西安 710065；4.国网四川省电力公司 经济技术研究院，成都 610041)

1 研究背景

岩质边坡稳定性历来是岩土工程领域的热门研究内容。早期人们基于地质条件对边坡稳定性进行定性评价；随后极限平衡法的引入使岩质边坡稳定性研究步入定量评价阶段；20世纪70年代以来，计算机技术的迅速发展使得以有限元、离散元等为代表的数值分析方法在边坡稳定性研究中得到大量运用[1]。

以往研究中，无论是极限平衡法还是有限元等连续变形数值分析方法，多将岩质边坡概化成连续均质体，未能考虑岩体是各种结构面切割而成的非均质地质体这一事实；离散元、DDA[2]等非连续变形方法将岩体抽象为离散刚性单元的集合体[3]，适用于分析贯通性结构面网络切割的岩体；然而岩体中的结构面——尤其是Ⅳ类、Ⅴ类结构面通常是随机、断续发育，结构面末端多终止于完整岩石，结构面间以岩桥衔接，因此离散元等将岩体抽象成独立、分离块体组合的做法也不能很好地反映真实岩体中结构面非贯通、随机发育的特性。

耦合了节理网络的有限元分析方法[4](jointed finite element method)则很好地弥补了以上方法的不足。这种方法基于野外调查获得结构面特征参数，采用Monte-Carto法模拟生成岩体中的结构面网络，在此基础上建立有限元计算模型，可分析边坡应力、变形及塑性区特征，进行有限元强度折减计算还能给出边坡的稳定性系数、分析边坡的变形破坏机制。因此，节理网络有限元法在岩石工程领域具有广阔的应用前景，目前在国外已得到较广泛的运用，如Hammah等[5]、Ma等[6]将其应用于岩质边坡稳定性分析，SatcÖ[7]等将其运用于地下工程，国内在这方面的相关研究及应用[8-9]处于起步阶段。

本文以地处旱寒地区的新疆某抽水蓄能电站库区花岗岩边坡为研究对象。基于现场调查，利用节理网络有限元方法对该边坡进行稳定性计算、变形破坏模式分析。研究成果对该边坡的处置具有指导意义、对旱寒地区花岗岩边坡变形破坏模式研究具有一定的参考价值，研究方法为节理岩质边坡稳定性分析提供了新的思路与参考。

2 基本原理及方法

2.1 节理网络模拟

岩体中的结构面成因不同、规模各异，是成岩建造、构造及浅表生作用的综合产物。早期人们提出一系列简化的节理网络模型(图1)模拟岩体中的结构面，如平行确定(Parallel Deterministic)模型、平行统计(Parallel Statistical)模型、交错(Cross Jointed)模型、Voronoi泰森多边形模型等。平行确定模型及平行统计模型适用于模拟结构面产状稳定、在空间上呈平行状态的节理网络，前者的特点是节理平行、等间距，后者的特点是节理平行但间距不等；交错模型适用于模拟连贯岩层层面及层内节理构成的结构面网络；Voronoi模型则随机地将二维平面离散成一个个非重叠的凸多边形，多边形边界即为离散的节理网络，适用于模拟风化剧烈、节理密布相互交切且产状凌乱的岩体。

图1 4种常见节理网络模型Fig.1 Four common models of joint network

上述简化模型的运用均受一定的限制,Baecher模型[10-11]则较为灵活，可以考虑任意结构面位置、方位、规模的结合，因而其不仅能模拟简单贯通结构面网络，更能快速模拟一系列复杂的非贯通、随机结构面网络。该模型假定结构面中心在空间上的分布服从泊松点过程，结构面产状、长度、发育密度等服从一定的统计规律，在野外结构面实测的基础上，用统计学方法获取相关参数，通过Baecher算法生成结构面网络。

在Baecher模型中，需获取或指定如下参数。

(1)结构面方位(产状)：一般可通过Fisher概率函数或定义倾向/倾角以确定空间结构面的方位。采用倾向/倾角法，一旦平均倾角在0°或90°附近，或者平均倾向接近0°或360°时，则可能无法得到理想的结构面分布结果。采用Fisher函数定义结构面方位则不受限制，结构面产状用其法线矢量的方向表示，采用球面坐标时某组结构面的法向量将围绕其平均法线按Fisher概率函数分布[12],即

(1)

式中：φ′为概率密度函数；θ为各结构面矢量与结构面组平均矢量的角度差(°)；K称“Fisher常量”，反映了结构面产状的离散程度，K值越大，结构面产状离散性越小。

(2)结构面长度：实测结构面长度一般服从某种分布函数规律，可用正态分布、均匀分布、指数分布等概率密度函数描述其长度规律。

(3)结构面密度：结构面密度反映了岩体内节理的发育程度，常用单位面积内结构面条数(1/m2)、单位面积内结构面总迹长(1/m)、单位面积内结构面总迹长与面积平方根之比等。

此外，节理网络有限元分析中对节理末端条件的定义极为重要，末端条件规定了节理的运动行为。当末端被指定为“闭合”时，采用1个有限元节点代表节理末端，此时节理末端不能发生相对运动；当节理末端指定为“张开”时，采用2个节点代表节理的末端，节理两侧单元可以发生相对运动(滑动、拉开等)。图2为一个“收缩”了的节理网络有限元模型，结构面左下端定义为“闭合”，则结构面端点只有一个有限元节点1；右上端为“张开”，则末端共用2号和3号2个节点，允许2个节点发生相对运动可以模拟节理剪切错动、拉张等运动行为。

图2 节理末端条件Fig.2 End condition of joint

2.2 有限元强度折减法

近年来，有限元法强度折减法得到了大量的运用[13]。强度折减计算中，定义岩土体剪切强度折减系数SRF，按式(2)对岩土体剪切强度参数进行折减。

(2)

式中：c与c′、φ与φ′分别表示折减前后的内聚力、折减前后的内摩擦角。

本文采用计算位移作为收敛判据。具体而言：先利用原始强度参数(SRF=1)进行计算，一旦计算结果收敛，则增大折减系数，利用折减后的强度参数再次进行计算，如此不断地对强度参数进行折减，直到在指定的迭代步数及容差范围内计算不收敛则结束，最终的折减系数即是边坡的稳定性系数。对于稳定性系数<1的边坡，采用原始强度参数进行计算时必然是不收敛的；对此，从SRF=1开始逐步减小折减系数(相当于不断增大材料的强度参数)，直至计算结果收敛，最终这个<1的折减系数即为稳定性系数。因此从后一种情况的角度出发，人们常说的“强度折减法”并不严谨，只是由于实际评估中岩土体稳定性系数一般≥1，因而“强度折减”这一说法得以流传。

2.3 基于节理网络有限元进行稳定性分析的步骤

主要包括现场调查、结构面网络模拟及有限元分析几个步骤。

(1)现场调查及测量。现场工作阶段应查明边坡工程地质条件、风化卸荷特征、变形迹象等；查明边坡岩体中结构面类型、规模及力学性状，采用测窗法、测线法等对坡体结构面进行测量统计。

(2)建立结构面概率模型。将实测结构面按类型、规模进行分类，分别根据优势方位对各级结构面进行分组，统计各组结构面产状、长度、密度等特征参数，建立其概率模型。

(3)节理网络模拟。在正确认识边坡地质条件的基础上，对地质原型进行概化；根据边坡结构面的发育特点，选择合适的模型生成三维结构面网络，沿指定方向剖切得到二维结构面网络。

(4)对耦合了节理网络的有限元模型进行计算，结合现场调查对计算成果进行综合分析。

3 应用实例

3.1 工程概况及基本地质条件

某拟建抽水蓄能电站位于新疆东部，工程区属典型的旱寒地区，气候具有“降雨稀少、蒸发量大，年均气温低、昼夜温差大”的特点。

电站由上、下2个蓄水水库及其他设施组成，2个水库间通过永久道路连接；沿线花岗岩高边坡众多，受风化、卸荷作用的强烈影响，边坡岩体松动。这些边坡普遍存在稳定性问题，其中某段边坡更是发生了较大规模的崩塌(图3)，为此有必要对该边坡的稳定性、失稳模式及诱发因素进行研究，为工程后续处置提供依据。

图3 坡脚崩落的巨石Fig.3 Huge cavings at the foot of slope

该边坡高约170 m，坡向330°，平均坡度约40°，上部陡坡向外微凸，中下部较平直。坡体主要由华力西期(γ42b)中粒-粗粒块状花岗岩构成，坡脚堆积第四系(Q4)崩坡积、冲洪积块碎石混杂物。边坡区发育的结构面类型单一，主要以Ⅳ类、Ⅴ类硬性结构面(节理)为主，发育3组优势节理，最发育一组(J1)产状280°～355°∠60°～84°；次发育组(J2)产状140°～190°∠40°～70°；发育最少的缓倾角节理组(J3)产状300°～340°∠8°～25°。除J3组节理外，其余2组节理在边坡不同风化带内的性状(密度、张开度、填充及力学特性)差异较大；浅表岩体受强烈风化及卸荷的影响，陡倾角裂隙极为发育，一般张开度0.2～2 cm，充填岩屑，最大卸荷裂隙宽约15 cm，充填岩块。结合附近场区勘探平硐调查结果将边坡划分为强风化带(厚20～30 m)、弱风化带(厚15～20 m)及微新岩体，风化界线大致平行坡面。

由于该花岗岩边坡节理极为发育，节理对边坡稳定性的控制作用极为显著，因此采用节理网络有限元方法对其稳定性进行分析。

3.2 节理网络模拟及有限元计算

现场采用测窗法[14]测量120余条结构面的几何参数，并对其张开度、填充特征等进行精细描述。通过实测节理产状编制结构面等密图，将边坡区结构面划分为3组，统计得到各组结构面的几何特征参数(表1)。表1中结构面迹长服从正态分布。不同风化带内结构面倾向及倾角均值相同，但其产状离散程度、平均迹长及密度有所差异。

表1 结构面几何特征参数Table 1 Geometry parameters of discontinuities

注：强风带内节理末端条件在地表“张开”，其余全为“闭合”

采用Baecher模型生成三维节理网络，按边坡倾向剖切得到的二维节理网络模型见图4。

图4 生成的节理网络Fig.4 Generated joint network model

3.3 有限元建模及计算

在结构面网络模拟的基础上，对该边坡进行有限元分析。计算模型两侧采用滑动边界，底边采用约束边界，地表为自由面。计算中只考虑坡体自重应力场，采用三角形网格单元进行离散化。岩土体材料及结构面均采用摩尔-库伦屈服准则，计算材料参数见表2、表3。由于边坡区位于旱寒地区，降雨稀少，因此不考虑暴雨工况，仅考虑天然及地震工况。按相关规范取地震水平加速度0.15g并折减0.25倍，即实际取水平地震加速度0.037 5g，不考虑竖向加速度。

表2 岩体参数取值Table 2 Parameters of intact rock

表3 结构面参数取值Table 3 Parameters of discontinuities

3.4 计算结果分析

3.4.1 总体稳定性

计算得该边坡天然、地震工况下的稳定性系数分别为1.57和1.46，计算工况下处于整体稳定状态。

图5(a)反映了天然工况下对材料剪切强度折减至边坡处于临界状态时(折减系数SRF=1.57)，在中上部凸出的强风化岩体内形成了一个几近贯通的最大剪应变集中带，这个弧形带指示了边坡整体失稳时的潜在滑带。与最大剪应变集中带相对应，潜在滑体发生了显著的位移，位移云图(图5(b))总体呈追踪陡倾结构面和缓倾结构面发展的“锯齿”状，位移矢量(红色箭头)指示了边坡的变形方向，临界状态下强风化带内的陡倾节理几乎全部屈服(图5(b)中红色直线为屈服节理，黑色直线为未屈服节理)，屈服的缓倾节理相对较少。

图5 最大剪应变、最大总位移云图(SRF=1.57)Fig.5 Contours of maximum shear strain and maximum total displacement (SRF=1.57)

3.4.2 变形破坏演化模式

为研究该边坡变形破坏的演化模式，以图5(b)中的A区为例进行重点分析，放大的A区最大总位移云图见图6。

图6 A区最大总位移云图(SRF=1.57)Fig.6 Contours of maximum total displacement of zone A in Fig.5 (b) when SRF=1.57

原始状态下(即SRF=1时)，A区岩体未发生显著的变形，仅部分节理屈服，且屈服的陡倾角节理多于缓倾角节理；随着强度折减系数逐渐增大，绝大多数节理逐步发生屈服；被切割的岩块沿节理面剪切错动(节理两侧箭头表示相对错动方向)，图中缓倾节理①上方岩块向临空方向滑移，伴随这一过程的发生，与节理①交切的陡倾节理被拉开成开口下大上小的“楔形”，其演化趋势是以“滑移-拉裂”[15]的形式剪出失稳。

因此，“滑移-拉裂”是该边坡变形破坏的主要演化模式。具体而言，在自重等因素的长期作用下，边坡岩体将沿缓倾坡外的J3节理组发生剪切滑移，伴随这一过程的发生，陡倾的J1及J2节理组发生拉裂为主、剪切为次的拉-剪破坏，应力逐渐向完整岩桥集中，随着岩桥逐步破裂、贯通，滑面形成。

现场调查的结果证实了这类变形破坏模式的存在。图7(a)中的危岩体边界主要受由2组节理控制，后缘陡倾节理产状312°∠70°，底滑控制节理产状330°∠10°，底滑面与陡倾节理相交处普遍被拉开，充填岩屑。图7(b)中岩块滑移使陡倾结构面显著拉开2～5 cm，拉裂缝开口下大上小，节理交切处明显架空，充填岩屑及岩块。

图7 典型滑移-拉裂节理组合Fig.7 Typical joint combination of shearing and cracking

3.5 稳定性综合评价及防治措施建议

现场调查表明该花岗岩边坡无贯穿的长大控制性结构面，坡体主要发育Ⅳ类、Ⅴ类硬性结构面；有限元强度折减计算结果表明该边坡的整体稳定性较好。

该花岗岩边坡位于旱寒地区，昼夜温差大，气温骤变及由此产生的冻、融为主导的物理风化作用十分强烈[16]。强烈风化使岩体产生新裂隙的同时不断扩大已有裂隙，岩体完整性遭到不断破坏，在裂隙“冰劈”(图8)、地震等作用下易发生崩塌灾害；冰劈诱发的崩塌多发生于夜晚，显然这与裂隙水夜冻昼融有关。

综上，该边坡具有整体稳定性好、局部强风化松动块体稳定性差的特点。针对其现状稳定性、变形演化特点及影响诱发因素，建议采取如下防治措施：

(1)人工清坡，翘落或采取小型爆破清除坡面松动危岩。

(2)对难以清除的较大规模不稳定块体，采取后缘裂缝灌浆的措施，以防“冰劈”作用下块体失稳；具备施工条件的，宜对大的危险块体采取锚固措施。

(3)坡面重点区域挂网喷混凝土，以减弱坡体浅表物理风化作用，防止坡面掉块；下方公路内侧设置多道被动拦截网。

4 结 论

(1)岩质边坡稳定性分析中，传统的做法多是将边坡概化成连续均质体或离散的刚性块体集合。前者过于简化，忽略了岩体中结构面的影响；后者则在一定程度上夸大了结构面的控制效应，难以反映真实岩体中结构面非贯通、随机发育的特点。

(2)节理网络有限元法为岩质边坡稳定性分析提供了新的有效手段，该法基于现场结构面测量统计，采用节理网络模拟的方法可以高度还原边坡中结构面的发育特征，在此基础上进行有限元分析，结果更加真实可靠。

(3)采用节理网络有限元分析方法，不仅能了解边坡应力场、位移场、塑性区分布等特征，进行强度折减计算，还能给出边坡的稳定性系数、分析边坡变形破坏的范围及演化模式，这为人们评估边坡的现状、预测边坡的未来提供了有力帮助。

(4)本文以地处旱寒地区的某花岗岩边坡为例进行了节理网络有限元分析。结果表明该边坡处于整体稳定状态，浅表节理切割、强烈风化的松动岩体稳定性较差，易在外界因素诱发下发生滑移-拉裂式失稳。